CC BY
37
Таким образом, кава-фильтр типа «Зонтик» вызывает замедление кровотока в области его установки. Сильное уменьшение скорости течения становится причиной нового тромбообразования, поэтому необходимо задавать имплантату такие параметры, чтобы он был «прозрачен» для кровотока. С помощью объединения моделей глобальной циркуляции крови и эластичной стенки сосуда становится возможным наблюдение за влиянием кава-фильтров с различными параметрами на гемодинамику, а также определение их оптимальных размеров.
Работа поддержана РФФИ, грант 08-01-00159, 09-01-00115, а также программой Президиума РАН 21-П «Фундаментальные науки в медицине».
Литература
1. Есикова Н.Б., Мухин С.И., Соснин Н.В., Фаворский А.П., Хруленко А.Б. Математическое моделирование течения крови с кава-фильтрами. Препринт. - М.: МАКС Пресс, 2004.
2. Холодов А.С., Симаков С.С. Численное исследование содержания кислорода в крови человека при низкочастотных воздействиях // Математическое моделирование. - 2008. - Т. 20. - № 4. - С. 87-102.
3. Каро К., Педли Т., Шротер Р., Сид У. Механика кровообращения / Под ред. С.А. Регирер, В.М. Хаю-тина. - М.: Мир, 1981.
4. Rosar M.E., Peskin C.S. Fluid Flow in Collapsible Elastic Tubes: A Three-Dimensional Numerical Model // New York J. Math. - 2001. - № 7. - Р. 281-302.
5. Савельев В.С., Гологорский В.А., Кириенко А.И. и др. Флебология / Под ред. В.С. Савельева - М.: Медицина, 2001.
6. Cheng Tu and Charles Peskin. Stability and Instability in the Computation of Flows with Moving Immersed Boundaries: A Comparsion of Three Methods // SIAM J. SCI. STAT. COMPUT. - November 1992. - Vol. 13. - № 6. - Р. 1361-1376.
Иванов Юрий Александрович - Московский государственный университет им. М.В. Ломоносо-
ва, студент, yura-vtm@yandex.ru Добросердова Татьяна Константиновна - Московский государственный университет им. М.В. Ломоносова, студент, DobroserdovaTK@gmail.com
УДК 519.711
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ТУБЕРКУЛЕЗА В.Н. Леоненко, К.К. Логинов
В работе рассматриваются подходы, позволяющие повысить быстродействие моделирующей программы для индивидуум-ориентированной модели распространения туберкулеза. Описывается процесс оптимизации вычислительных алгоритмов, а также использование распределенных вычислений на основе системы МОЫС и параллельных вычислений с применением технологии ОрепМР. Приведены результаты вычислительных экспериментов. Ключевые слова: имитационное моделирование, эпидемиология, распределенные вычисления, параллельное программирование.
Введение
Одним из инструментов исследований закономерностей распространения туберкулеза является метод математического моделирования. В настоящее время разработан широкий класс детерминированных моделей, описывающих распространение туберкулеза и контроль за этим заболеванием [1, 2]. Вместе с тем, практически отсутствуют модели, опирающиеся на стохастический подход и индивидуально-ориентированное описание неинфицированных, инфицированных и больных индивидуумов.
Индивидуум-ориентированные модели позволяют детально провести формализацию принципов передачи и развития туберкулеза и исследовать динамику заболеваемости туберкулезом от различных факторов. Одними из первых в этом направлении являются модели, представленные в работах [3, 4]. В работе [3] предложена стохастическая модель в форме высокоразмерной системы рекуррентных уравнений с дискретным временем, учитывающая зависимость развития заболевания от возраста индивидуумов. Модель, построенная в работе [4], имеет следующие отличительные свойства: используется непрерывное время, продолжительность пребывания индивидуумов в различных стадиях болезни задается с помощью случайных величин с произвольными законами распределения, вероятность контактов индивидуумов считается пропорциональной произведению численностей когорт.
В данной работе рассматривается стохастическая модель [5], созданная на основе дифференциальной модели динамики туберкулеза [2]. В индивидуум-ориентированной модификации модели, в отличие
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ...
от предшествующих работ [3, 4], учитывается влияние продолжительности болезни на вероятность гибели индивидуумов. Введение индивидуальных параметров привело к тому, что объем обрабатываемых данных повысился и время, затрачиваемое на моделирование динамики популяции, существенно возросло. Основной целью данной работы является рассмотрение подходов, позволяющих повысить быстродействие моделирующих программ путем оптимизации исходного алгоритма и более эффективного применения имеющихся компьютерных мощностей и вычислительных технологий.
Описание модели
Кратко опишем особенности модификаций модели, рассматриваемых в данной работе (полное формальное описание приведено в [5]).
Базовая модель. В популяции выделяется шесть когорт населения: восприимчивые индивиды (5), инфицированные индивидуумы (Ь), невыявленные больные индивидуумы без бактериовыделения (В), невыявленные больные индивидуумы с бактериовыделением (В), выявленные больные индивидуумы без бактериовыделения (В0) и выявленные больные индивидуумы с бактериовыделением (В0). Система уравнений модели имеет следующий вид:
(*) = % (1) - и8, Ь^В (1) + / (1),
Хь V) = Хь (Г) + иаЬ (Г) + ивЬЬ (Г) + иВаь (Г) - ив,В (Г) + /2 (Г),
хВ (*) = *в (*) + и8 В (*) + иЬ В V) + ив В ^) - ив Ь ^) - ив в V) - ив В0 (1) + /з (1),
' ' ' ' ' ' ° (1)
Хв V) = Хв (() + ив в V) - ив, В ({) - ив, в0 (*) + /а (1),
ХВ„ V) = ХВ0 V) + ив, в0 V) + ив0, в0 (0 - ив0 ,в„ V) - иВ0 , Ь V) + /5 V X
Хв0 V) = Хв0 V) + ив,в0 V) + ив0,в0 V) -ив0,В0 V) + /6VXt = ^Х...^
Система (1) дополняется начальными данными при t = 0. Здесь х¥ Ц) - численность индивидуумов когорты У в момент времени Хг У) - численность индивидуумов когорты У , доживших от момента времени t -1 до момента времени t, случайные величины вида иУ г ^) описывают переход из когорты У
в когорту 7 в течение промежутка времени ^-1,t], У, 1 е {а,Ь,В,в,В0,в0}. Слагаемые /(^, I = 1,6 задают приток населения в соответствующие когорты за промежуток времени ^ -1, t ].
Для описания процессов инфицирования индивидуумов и спонтанной активации болезни используется схема, предложенная в работе [3]. Инфицирование происходит при одновременном посещении больными и здоровыми индивидуумами мест возможного контакта. Считается, что за единицу времени индивидуумы из когорт в , во посещают в совокупности ^{в) и в'о) мест соответственно, где
^(в) = X(в), ^(во) = XХвв^У*(во), Уй(в) ^ 0 , (во) > 0 - взаимно независимые, одинаково распределенные случайные величины, не зависящие от хв^), хв ^), причем Е(уд(в>) = г >0, Е(уд(во)) = к.г , к е (0;1). Вероятности того, что индивидуумы из когорт а, Ь посетят одно из этих мест, обозначим через Xх е (0; 1) и XЬ е (0; 1) соответственно. Вероятности заражения индивидуумов после контакта с больным индивидуумом обозначим через 5Х е (0; 1) и 5Ь е (0; 1) соответственно. Тогда вероятности того, что индивидуум из конкретной когорты а, Ь будет инфицирован в течение суток, при фиксированном значении | t задаются следующими формулами: ц; ) =1-(1-Хх5Х)t t , ц; ) =1-(1-0ь)(1-ХЬ5Ь)t t , где 6Ь - вероятность спонтанной активации болезни для индивидуумов когорты Ь. Данные величины определяют количество инфицированных на интервале ^ -1, t ].
Индивидуум-ориентированная модель. В популяцию вводятся когорты Я и Я0, которые используются для описания индивидуумов, находящихся в состоянии ремиссии. Считается, что в когорту Я поступают самоизлечившиеся индивидуумы из когорты В , а в когорту Я0 - вылеченные индивидуумы из когорты В0. Переходы индивидуумов из В и В0 в Ь заменяются на переходы из В в Я и из В0 в Я0.
Предполагается, что больные индивидуумы после перехода в состояние ремиссии остаются инфицированными и существует риск повторного развития заболевания. Это означает, что возможны переходы индивидуумов из Я в В и из Я0 в В0.
Для индивидуумов когорт В, в, В0, в0, Я и Я0 вводится параметр ДХ^) = тах{0,пх^) -стх^)}, t = 0,1,2,..., отражающий влияние продолжительности болезни индивидуума на вероятность его дожития от момента времени t -1 до г. Здесь пх ^) = уВ аВх ^) + ув авх ^) + уВ ®ВоХ^) + ув ^) - величина,
отражающая негативное влияние заболевания на продолжительность жизни индивидуума x, ст x (t) = 5R arx (t) + 5 ^ x (t) - величина, отражающая снижение эффекта воздействия заболевания на
продолжительность жизни индивидуума x за счет его пребывания в состоянии ремиссии,
УD, Ув, Уd , YBo, 5r , 5r - положительные константы, rndx (t), rnbx (t), (t), roB (t), roRx(t), ra^x (t) -
продолжительности пребывания индивидуума x в соответствующих когортах до момента времени t включительно. Вероятность дожития от момента времени t -1 до t для индивидуума x из когорты H
при заданной Дx (t) равна рH x (t) = рHe Ax (t 1)1, H e {D, B, D0, B0, R, R0} (в отличие от базовой модели, где
Рн,x (t) = Рн = const).
Повышение быстродействия работы моделирующих программ
Рассмотрим далее способы ускорения работы моделирующих программ, позволяющих исследовать динамику численностей когорт в зависимости от параметров модели.
Оптимизация вычислительных алгоритмов. Рассмотрим индивидуум-ориентированную модификацию модели. Для всех вычислительных алгоритмов, использованных в моделирующей программе, справедливо следующее.
1. Когорты S и L полностью характеризуются своими численностями.
2. Динамика когорт S и L моделируется розыгрышем общей численности этих когорт на каждом шаге с помощью заданных вероятностей дожития индивидуумов.
Основные сложности возникают при моделировании динамики когорт, члены которых обладают индивидуальными параметрами.
Пусть T = {D, B, D0, B0, R, R0}. Наиболее простым является итеративный вариант моделирования всех взаимодействий, выглядящий следующим образом. Алгоритм № 1.
1. Для каждой когорты множества T создается список, содержащий индивидуальные параметры членов когорты.
2. Дожитие и переходы индивидуумов когорт множества T в течение суток (т.е. в течение одной итерации цикла) разыгрываются для каждого из этих индивидуумов последовательно в процессе прохода по спискам параметров.
Заметим, что вследствие линейности процессов дожития и переходов между когортами для индивидуумов когорт D, B, D0, B0 судьбу каждого такого индивидуума можно разыгрывать независимо. Тогда алгоритм моделирования принимает следующий вид. Алгоритм № 2.
1. Для когорт R и R создаются списки, содержащие индивидуальные параметры членов когорт.
2. Процессы дожития и переходов индивидуумов когорт R и R0 в течение суток разыгрываются для каждого из этих индивидуумов последовательно путем прохода по спискам параметров.
3. Процессы дожития и переходов индивидуумов когорт D, B, D0, B0 вплоть до их гибели или выздоровления разыгрываются для каждого индивидуума с помощью функции розыгрыша судьбы отдельно взятого индивидуума. В процессе розыгрыша принадлежность индивидуума к заданной когорте в заданный момент времени регистрируется в массиве численностей индивидуумов множества T. В случае перехода рассматриваемого индивидуума в состояние ремиссии значение его параметра передается в качестве выходного параметра функции и помещается в массив параметров индивидуумов когорт R и R0 .
Формальное описание алгоритмов №№ 1, 2 для базовой модели и алгоритма № 1 для индивидуум-ориентированной модели приведено в работе [5].
Функция розыгрыша судьбы отдельно взятого индивидуума реализована в двух вариантах.
(а) Итеративный вариант. Судьба индивидуума разыгрывается итеративно, с шагом в одни сутки (на каждом шаге с помощью метода Монте-Карло моделируется, погибнет ли индивидуум на данном шаге, и если нет, то перейдет ли он в другую когорту).
(б) Вариант со скачками. На основании значений вероятностей гибели индивидуума и перехода его в другую когорту разыгрывается момент ближайшей смены статуса рассматриваемого индивидуума, после чего внутри функции происходит «скачок» во времени до этого момента с одновременным пересчетом значения индивидуального параметра.
Для сравнения быстродействия работы моделирующих программ проводился вычислительный эксперимент (табл. 1). Расчеты производились для двух наборов параметров на персональном компьютере Intel Core2 Quad 2,66 ГГц, объем оперативной памяти - 3 Гб. Отрезок моделирования - 5500 суток. Получаем, что моделирующая программа, выполненная на основе алгоритма №2 и функцией розыгрыша
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ.
судьбы в варианте (б), производит расчеты с первым набором параметров примерно на 70% быстрее, а со вторым набором - на 55% быстрее по сравнению с исходным вариантом алгоритма.
Заметим, что реализация розыгрыша судьбы индивидуума в виде отдельной функции также облегчает модификацию программного кода при необходимости изменения законов протекания инфекции (в перспективе возможна реализация выбора из нескольких модулей протекания инфекции и схем заражения для разных типов заболеваний).
Время работы на наборе № 1, с Время работы на наборе № 2, с
Алгоритм №1 1611,39 2468,01
Алгоритм №2 (а) 970,14 1636,69
Алгоритм №2 (б) 478,43 1111,56
Таблица 1. Быстродействие моделирующей программы с разными алгоритмами
Распределенные вычисления с использованием системы MONC
Система MONC [6] (сокращение от Monte Carlo), разработанная М.А. Марченко (Институт вычислительной математики и математической геофизики СО РАН, г. Новосибирск), является универсальной высокопроизводительной системой параллельных вычислений для методов Монте-Карло на базе сети персональных компьютеров. Она распределяет независимые копии задания по персональным компьютерам в сети, отдает команды на их исполнение, следит за ходом выполнения заданий и по завершению выполнения заданий выполняет копирование и осреднение файлов с результатами расчетов.
Был проведен вычислительный эксперимент с использованием системы MONC на комплексе из четырех компьютеров с ОС Windows XP. Исполняемым приложением была модификация моделирующей программы для индивидуум-ориентированной модели, доработанная с учетом требований, предъявляемых системой (см. [6]). Производилось нахождение оценок математических ожиданий численностей когорт по ста реализациям для двух различных наборов параметров (отрезок моделирования - 5500 суток). Сначала измерялось время выполнения расчетов на каждом из компьютеров по отдельности, затем общее количество моделируемых реализаций разделялось по компьютерам пропорционально их быстродействию. После окончания расчетов на всех машинах система MONC усредняла полученные данные и сохраняла результат.
Таким образом, в данном случае запуск модели с помощью системы MONC на комплексе из четырех компьютеров позволяет сократить время выполнения задания примерно на 66% для обоих наборов параметров по сравнению с запуском на самом быстром из имеющихся компьютеров (O94-11).
Сетевое имя Процессор Оперативная память Время работы на наборе № 1, с Время работы на наборе № 2, с
O94-31 Intel Pentium 4, 3,2 ГГц 512 Мб 5156,7 17307,2
O91-01 AMD Athlon 64 x2 Dual Core Processor 5200+, 2,61 ГГц 1 Гб 4958,23 16755,6
O94-09 Intel Pentium 4, 3 ГГц 512 Мб 5826,1 19228,4
O94-11 Intel Core 2 Duo E7200, 2,53ГГц 1 Гб 3679,15 12031,9
Таблица 2. Характеристики компьютеров, задействованных для расчетов с использованием системы МОЫО
Набор № 1 Набор № 2
Сетевое Кол-во реализаций Время выполнения расче- Кол-во реализа- Время выполнения расче-
имя тов, с ции тов, с
O94-31 23 1216 23 3885
O91-01 24 1233 24 4073
O94-09 21 1183 20 3850
O94-11 32 1171 33 4019
Итого 100 1233 100 4073
Таблица 3. Результаты вычислительного эксперимента по применению системы MONC
Параллельные вычисления с использованием технологии OpenMP
Технология OpenMP (Open Multi-Processing) представляет собой набор директив компилятора, функций и переменных окружения, позволяющих создавать на базе последовательных программ многопоточные приложения, предназначенные для выполнения на многопроцессорных системах с единой па-
мятью. Важным достоинством OpenMP является то, что исходная программа не переписывается заново, а модифицируется добавлением в текст соответствующих директив [5].
К моменту написания статьи реализована поддержка технологии OpenMP в моделирующей программе для базовой модели с алгоритмом №1. Проводился вычислительный эксперимент с двумя наборами параметров, расчеты выполнялись на компьютере Intel Core 2 Duo E7400 2.8ГГц с объемом оперативной памяти 2 Гб. Отрезок моделирования - 5500 суток.
Количество потоков Время работы на наборе № 1, с Время работы на наборе № 2, с
1 129,968 191,093
2 71,578 99,515
Таблица 4. Результаты испытаний программы с поддержкой технологии OpenMP
Проведение вычислений в одном потоке соответствует выполнению последовательной версии моделирующей программы. Таким образом, применение технологии OpenMP при выполнении программы на двухъядерном процессоре дает экономию времени в 45% для первого набора параметров и 48% для второго набора параметров по сравнению с последовательным вычислением.
Заключение
В работе были описаны подходы, которые позволили существенно увеличить скорость работы моделирующей программы. Это особенно актуально в связи с прогнозируемым усложнением модели (введением возрастной структуры в популяции и реализацией более детального механизма заражения). Также важной задачей является освоение моделирования на высокопроизводительных ЭВМ с распределенной памятью, в связи с чем планируется адаптировать программу под распределенные вычисления на основе технологии MPI.
Авторы благодарят своего научного руководителя профессора Н. В. Перцева (ОФ ИМ им. С. Л. Соболева СО РАН) за постановку задачи и обсуждение результатов работы.
Литература
1. Авилов К. К., Романюха А. А. Математические модели распространения и контроля туберкулеза (обзор) // Математическая биология и биоинформатика. - 2007. - Т. 2. - № 2. - С. 188-318.
2. Perelman M.I., Marchuk G.I., Borisov S.E., et. al. Tuberculosis epidemiology in Russia: the mathematical model and data analysis // Russ. J. Numer. Anal. Math. Modelling. - 2004. - V. 19. - № 4. - Р. 305-314.
3. Перцев Н.В., Романюха А.А., Касаткина В.С. Нелинейная стохастическая модель распространения туберкулеза // Системы управления и информационные технологии. - 2008. - № 1.2 (31). - С. 246250.
4. Перцев Н.В., Пичугин Б.Ю. Индивидуум-ориентированная стохастическая модель распространения туберкулеза // Сибирский журнал индустриальной математики. - 2009. - Т. 12. - № 2(38). - С. 97110.
5. Pertsev N.V., Leonenko V.N. Stochastic individual-based model of spread of tuberculosis // Russ. J. Numer. Anal. Math. Modelling. - 2009. - V. 24. - № 4. - Р. 341-360.
6. Марченко М.А. Комплекс программ MONC для распределенных вычислений методом Монте-Карло. [Электронный ресурс]. - Режим доступа: http://osmf.sscc.ru/~mam/monc_rus.htm, своб.
7. Антонов А.С. Параллельное программирование с использованием технологии OpenMP. - М.: Изд-во МГУ, 2009.
Леоненко Василий Николаевич - Омский филиал Института математики им. С. Л. Соболева СО РАН,
аспирант, VNLeonenko@yandex.ru Логинов Константин Константи- - Омский филиал Института математики им. С. Л. Соболева СО РАН, нович аспирант, kloginov85@mail.ru
