ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ТУБЕРКУЛЕЗА
В.Н. Леоненко, К.К. Логинов
УДК 519.711
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ТУБЕРКУЛЕЗА
В.Н. Леоненко, К.К. Логинов
В работе рассматриваются подходы, позволяющие повысить быстродействие моделирующей программы для индивидуум-ориентированной модели распространения туберкулеза. Описывается процесс оптимизации вычислительных алгоритмов, а также использование распределенных вычислений на основе системы MONC и параллельных вычислений с применением технологии OpenMP. Приведены результаты вычислительных экспериментов. Ключевые слова: имитационное моделирование, эпидемиология, распределенные вычисления, параллельное программирование.
Введение
Одним из инструментов исследований закономерностей распространения туберкулеза является метод математического моделирования. В настоящее время разработан широкий класс детерминированных моделей, описывающих распространение туберкулеза и контроль за этим заболеванием [1, 2]. Вместе с тем, практически отсутствуют модели, опирающиеся на стохастический подход и индивидуальноориентированное описание неинфицированных, инфицированных и больных индивидуумов.
Индивидуум-ориентированные модели позволяют детально провести формализацию принципов передачи и развития туберкулеза и исследовать динамику заболеваемости туберкулезом от различных факторов. Одними из первых в этом направлении являются модели, представленные в работах [3, 4]. В работе [3] предложена стохастическая модель в форме высокоразмерной системы рекуррентных уравнений с дискретным временем, учитывающая зависимость развития заболевания от возраста индивидуумов. Модель, построенная в работе [4], имеет следующие отличительные свойства: используется непрерывное время, продолжительность пребывания индивидуумов в различных стадиях болезни задается с помощью случайных величин с произвольными законами распределения, вероятность контактов индивидуумов считается пропорциональной произведению численностей когорт.
В данной работе рассматривается стохастическая модель [5], созданная на основе дифференциальной модели динамики туберкулеза [2]. В индивидуум-ориентированной модификации модели, в отличие
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
99
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ...
от предшествующих работ [3, 4], учитывается влияние продолжительности болезни на вероятность гибели индивидуумов. Введение индивидуальных параметров привело к тому, что объем обрабатываемых данных повысился и время, затрачиваемое на моделирование динамики популяции, существенно возросло. Основной целью данной работы является рассмотрение подходов, позволяющих повысить быстродействие моделирующих программ путем оптимизации исходного алгоритма и более эффективного применения имеющихся компьютерных мощностей и вычислительных технологий.
Описание модели
Кратко опишем особенности модификаций модели, рассматриваемых в данной работе (полное
формальное описание приведено в [5]).
Базовая модель. В популяции выделяется шесть когорт населения: восприимчивые индивиды (S),
инфицированные индивидуумы (L), невыявленные больные индивидуумы без бактериовыделения (D),
невыявленные больные индивидуумы с бактериовыделением (B), выявленные больные индивидуумы без
бактериовыделения (D0) и выявленные больные индивидуумы с бактериовыделением (B0). Система уравнений модели имеет следующий вид:
xS (t) xˆS (t) uS,LD (t) f1 (t),
xL (t) xˆL (t) uS,L (t) uD,L (t) uD0 ,L (t) uL,D (t) f2 (t),
xD (t) xˆD (t) uS,D (t) uL,D (t) uB,D (t) uD,L (t) uD,B (t) uD,D0 (t) f3 (t), xB (t) xˆB (t) uD,B (t) uB,D (t) uB,B0 (t) f4 (t),
(1)
xD0 (t) xˆD0 (t) uD,D0 (t) uB0 ,D0 (t) uD0 ,B0 (t) uD0 ,L (t) f5 (t),
xB0 (t) xˆB0 (t) uB,B0 (t) uD0 ,B0 (t) uB0 ,D0 (t) f6 (t), t 0,1,...,.
Система (1) дополняется начальными данными при t 0. Здесь xY (t) − численность индивидуумов
когорты Y в момент времени t, xˆY (t) − численность индивидуумов когорты Y , доживших от момента
времени t 1 до момента времени t , случайные величины вида uY ,Z (t) описывают переход из когорты Y
в когорту Z в течение промежутка времени (t 1,t] , Y , Z {S, L, D, B, D0 , B0} . Слагаемые fi (t) , i 1, 6 задают приток населения в соответствующие когорты за промежуток времени (t 1,t].
Для описания процессов инфицирования индивидуумов и спонтанной активации болезни используется схема, предложенная в работе [3]. Инфицирование происходит при одновременном посещении больными и здоровыми индивидуумами мест возможного контакта. Считается, что за единицу
времени индивидуумы из когорт
B,
B0
посещают в совокупности
(B) t
и
(B0 ) t
мест соответственно, где
(B) t
=
, =xB (t ) (B)
i =1 it
( B0 ) t
,xB0 (t ) (B0 )
i =1 it
(B) it
0,
(B0 ) it
0
− взаимно независимые, одинаково распреде-
ленные случайные величины, не зависящие от xB (t) , xB0 (t) , причем E(it(B) ) = r > 0 , E(it(B0 ) ) = k r , k (0;1) . Вероятности того, что индивидуумы из когорт S, L посетят одно из этих мест, обозначим через
S (0;1) и L (0;1) соответственно. Вероятности заражения индивидуумов после контакта с больным
индивидуумом обозначим через S (0;1) и L (0;1) соответственно. Тогда вероятности того, что индивидуум из конкретной когорты S, L будет инфицирован в течение суток, при фиксированном значении
t
задаются следующими формулами:
(S ) t
= 1 (1 S S )t( B) t(B0 )
,
(L) t
= 1 (1 L )(1 L L )t(B) t(B0 ) ,
где
L − вероятность спонтанной активации болезни для индивидуумов когорты L. Данные величины опреде-
ляют количество инфицированных на интервале (t 1,t].
Индивидуум-ориентированная модель. В популяцию вводятся когорты R и R0 , которые используются для описания индивидуумов, находящихся в состоянии ремиссии. Считается, что в когорту R поступают самоизлечившиеся индивидуумы из когорты D , а в когорту R0 − вылеченные индивидуумы
из когорты D0 . Переходы индивидуумов из D и D0 в L заменяются на переходы из D в R и из D0 в R0 . Предполагается, что больные индивидуумы после перехода в состояние ремиссии остаются инфицированными и существует риск повторного развития заболевания. Это означает, что возможны переходы индивидуумов из R в D и из R0 в D0 .
Для индивидуумов когорт D, B, D0, B0, R и R0 вводится параметр x (t) = max{0, x (t) x (t)} , t = 0,1, 2,, отражающий влияние продолжительности болезни индивидуума на вероятность его дожития
от момента времени t 1 до t . Здесь x (t) = D D,x (t) B B,x (t) D0 D0 ,x (t) B0 B0 ,x (t) − величина,
100
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
В.Н. Леоненко, К.К. Логинов
отражающая негативное влияние заболевания на продолжительность жизни индивидуума x , x (t) = R R,x (t) R0 R0 ,x (t) − величина, отражающая снижение эффекта воздействия заболевания на продолжительность жизни индивидуума x за счет его пребывания в состоянии ремиссии, D , B , D0 , B0 , R , R0 − положительные константы, D,x (t) , B,x (t) , D0 ,x (t) , B0 ,x (t) , R,x (t) , R0 ,x (t) − продолжительности пребывания индивидуума x в соответствующих когортах до момента времени t включительно. Вероятность дожития от момента времени t 1 до t для индивидуума x из когорты H
при заданной x (t) равна H ,x (t) = H ex (t1) , H {D, B, D0 , B0 , R, R0} (в отличие от базовой модели, где
H ,x (t) H 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 и R0 создаются списки, содержащие индивидуальные параметры членов когорт.
2. Процессы дожития и переходов индивидуумов когорт R и R0 в течение суток разыгрываются для каждого из этих индивидуумов последовательно путем прохода по спискам параметров. 3. Процессы дожития и переходов индивидуумов когорт D, B, D0, B0 вплоть до их гибели или выздоровления разыгрываются для каждого индивидуума с помощью функции розыгрыша судьбы отдельно взятого индивидуума. В процессе розыгрыша принадлежность индивидуума к заданной когорте в заданный момент времени регистрируется в массиве численностей индивидуумов множества T. В случае перехода рассматриваемого индивидуума в состояние ремиссии значение его параметра передается в качестве выходного параметра функции и помещается в массив параметров индивидуумов когорт R и R0 .
Формальное описание алгоритмов №№ 1, 2 для базовой модели и алгоритма № 1 для индивидуумориентированной модели приведено в работе [5].
Функция розыгрыша судьбы отдельно взятого индивидуума реализована в двух вариантах. (а) Итеративный вариант. Судьба индивидуума разыгрывается итеративно, с шагом в одни сутки (на каждом шаге с помощью метода Монте-Карло моделируется, погибнет ли индивидуум на данном шаге, и если нет, то перейдет ли он в другую когорту). (б) Вариант со скачками. На основании значений вероятностей гибели индивидуума и перехода его в другую когорту разыгрывается момент ближайшей смены статуса рассматриваемого индивидуума, после чего внутри функции происходит «скачок» во времени до этого момента с одновременным пересчетом значения индивидуального параметра.
Для сравнения быстродействия работы моделирующих программ проводился вычислительный эксперимент (табл. 1). Расчеты производились для двух наборов параметров на персональном компьютере Intel Core2 Quad 2,66 ГГц, объем оперативной памяти – 3 Гб. Отрезок моделирования – 5500 суток. Получаем, что моделирующая программа, выполненная на основе алгоритма №2 и функцией розыгрыша
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
101
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ...
судьбы в варианте (б), производит расчеты с первым набором параметров примерно на 70% быстрее, а со вторым набором – на 55% быстрее по сравнению с исходным вариантом алгоритма.
Заметим, что реализация розыгрыша судьбы индивидуума в виде отдельной функции также облегчает модификацию программного кода при необходимости изменения законов протекания инфекции (в перспективе возможна реализация выбора из нескольких модулей протекания инфекции и схем заражения для разных типов заболеваний).
Алгоритм №1 Алгоритм №2 (а) Алгоритм №2 (б)
Время работы на наборе № 1, с 1611,39 970,14 478,43
Время работы на наборе № 2, с 2468,01 1636,69 1111,56
Таблица 1. Быстродействие моделирующей программы с разными алгоритмами
Распределенные вычисления с использованием системы MONC
Система MONC [6] (сокращение от Monte Carlo), разработанная М.А. Марченко (Институт вычислительной математики и математической геофизики СО РАН, г. Новосибирск), является универсальной высокопроизводительной системой параллельных вычислений для методов Монте-Карло на базе сети персональных компьютеров. Она распределяет независимые копии задания по персональным компьютерам в сети, отдает команды на их исполнение, следит за ходом выполнения заданий и по завершению выполнения заданий выполняет копирование и осреднение файлов с результатами расчетов.
Был проведен вычислительный эксперимент с использованием системы MONC на комплексе из четырех компьютеров с ОС Windows XP. Исполняемым приложением была модификация моделирующей программы для индивидуум-ориентированной модели, доработанная с учетом требований, предъявляемых системой (см. [6]). Производилось нахождение оценок математических ожиданий численностей когорт по ста реализациям для двух различных наборов параметров (отрезок моделирования – 5500 суток). Сначала измерялось время выполнения расчетов на каждом из компьютеров по отдельности, затем общее количество моделируемых реализаций разделялось по компьютерам пропорционально их быстродействию. После окончания расчетов на всех машинах система MONC усредняла полученные данные и сохраняла результат.
Таким образом, в данном случае запуск модели с помощью системы MONC на комплексе из четырех компьютеров позволяет сократить время выполнения задания примерно на 66% для обоих наборов параметров по сравнению с запуском на самом быстром из имеющихся компьютеров (O94-11).
Сетевое имя
O94-31 O91-01
O94-09
Процессор
Intel Pentium 4, 3,2 ГГц AMD Athlon 64 x2 Dual Core
Processor 5200+, 2,61 ГГц Intel Pentium 4, 3 ГГц
Оперативная память
512 Мб 1 Гб
512 Мб
Время работы на наборе № 1, с 5156,7 4958,23
5826,1
Время работы на наборе № 2, с 17307,2 16755,6
19228,4
O94-11
Intel Core 2 Duo E7200, 2,53ГГц
1 Гб
3679,15
12031,9
Таблица 2. Характеристики компьютеров, задействованных для расчетов с использованием системы MONC
Сетевое имя O94-31 O91-01 O94-09 O94-11
Итого
Набор № 1 Кол-во реализаций Время выполнения расче-
тов, с 23 1216 24 1233 21 1183 32 1171
100 1233
Набор № 2 Кол-во реализа- Время выполнения расче-
ций тов, с 23 3885 24 4073 20 3850 33 4019
100 4073
Таблица 3. Результаты вычислительного эксперимента по применению системы MONC
Параллельные вычисления с использованием технологии OpenMP
Технология OpenMP (Open Multi-Processing) представляет собой набор директив компилятора, функций и переменных окружения, позволяющих создавать на базе последовательных программ многопоточные приложения, предназначенные для выполнения на многопроцессорных системах с единой па-
102
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
В.Н. Леоненко, К.К. Логинов
мятью. Важным достоинством OpenMP является то, что исходная программа не переписывается заново, а модифицируется добавлением в текст соответствующих директив [5].
К моменту написания статьи реализована поддержка технологии OpenMP в моделирующей программе для базовой модели с алгоритмом №1. Проводился вычислительный эксперимент с двумя наборами параметров, расчеты выполнялись на компьютере Intel Core 2 Duo E7400 2.8ГГц с объемом оперативной памяти 2 Гб. Отрезок моделирования – 5500 суток.
Количество потоков 1 2
Время работы на наборе № 1, с 129,968 71,578
Время работы на наборе № 2, с 191,093 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). – С. 246– 250.
4. Перцев Н.В., Пичугин Б.Ю. Индивидуум-ориентированная стохастическая модель распространения туберкулеза // Сибирский журнал индустриальной математики. – 2009. – Т. 12. – № 2(38). – С. 97– 110.
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
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
103
УДК 519.711
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАСПРОСТРАНЕНИЯ ТУБЕРКУЛЕЗА
В.Н. Леоненко, К.К. Логинов
В работе рассматриваются подходы, позволяющие повысить быстродействие моделирующей программы для индивидуум-ориентированной модели распространения туберкулеза. Описывается процесс оптимизации вычислительных алгоритмов, а также использование распределенных вычислений на основе системы MONC и параллельных вычислений с применением технологии OpenMP. Приведены результаты вычислительных экспериментов. Ключевые слова: имитационное моделирование, эпидемиология, распределенные вычисления, параллельное программирование.
Введение
Одним из инструментов исследований закономерностей распространения туберкулеза является метод математического моделирования. В настоящее время разработан широкий класс детерминированных моделей, описывающих распространение туберкулеза и контроль за этим заболеванием [1, 2]. Вместе с тем, практически отсутствуют модели, опирающиеся на стохастический подход и индивидуальноориентированное описание неинфицированных, инфицированных и больных индивидуумов.
Индивидуум-ориентированные модели позволяют детально провести формализацию принципов передачи и развития туберкулеза и исследовать динамику заболеваемости туберкулезом от различных факторов. Одними из первых в этом направлении являются модели, представленные в работах [3, 4]. В работе [3] предложена стохастическая модель в форме высокоразмерной системы рекуррентных уравнений с дискретным временем, учитывающая зависимость развития заболевания от возраста индивидуумов. Модель, построенная в работе [4], имеет следующие отличительные свойства: используется непрерывное время, продолжительность пребывания индивидуумов в различных стадиях болезни задается с помощью случайных величин с произвольными законами распределения, вероятность контактов индивидуумов считается пропорциональной произведению численностей когорт.
В данной работе рассматривается стохастическая модель [5], созданная на основе дифференциальной модели динамики туберкулеза [2]. В индивидуум-ориентированной модификации модели, в отличие
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
99
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ...
от предшествующих работ [3, 4], учитывается влияние продолжительности болезни на вероятность гибели индивидуумов. Введение индивидуальных параметров привело к тому, что объем обрабатываемых данных повысился и время, затрачиваемое на моделирование динамики популяции, существенно возросло. Основной целью данной работы является рассмотрение подходов, позволяющих повысить быстродействие моделирующих программ путем оптимизации исходного алгоритма и более эффективного применения имеющихся компьютерных мощностей и вычислительных технологий.
Описание модели
Кратко опишем особенности модификаций модели, рассматриваемых в данной работе (полное
формальное описание приведено в [5]).
Базовая модель. В популяции выделяется шесть когорт населения: восприимчивые индивиды (S),
инфицированные индивидуумы (L), невыявленные больные индивидуумы без бактериовыделения (D),
невыявленные больные индивидуумы с бактериовыделением (B), выявленные больные индивидуумы без
бактериовыделения (D0) и выявленные больные индивидуумы с бактериовыделением (B0). Система уравнений модели имеет следующий вид:
xS (t) xˆS (t) uS,LD (t) f1 (t),
xL (t) xˆL (t) uS,L (t) uD,L (t) uD0 ,L (t) uL,D (t) f2 (t),
xD (t) xˆD (t) uS,D (t) uL,D (t) uB,D (t) uD,L (t) uD,B (t) uD,D0 (t) f3 (t), xB (t) xˆB (t) uD,B (t) uB,D (t) uB,B0 (t) f4 (t),
(1)
xD0 (t) xˆD0 (t) uD,D0 (t) uB0 ,D0 (t) uD0 ,B0 (t) uD0 ,L (t) f5 (t),
xB0 (t) xˆB0 (t) uB,B0 (t) uD0 ,B0 (t) uB0 ,D0 (t) f6 (t), t 0,1,...,.
Система (1) дополняется начальными данными при t 0. Здесь xY (t) − численность индивидуумов
когорты Y в момент времени t, xˆY (t) − численность индивидуумов когорты Y , доживших от момента
времени t 1 до момента времени t , случайные величины вида uY ,Z (t) описывают переход из когорты Y
в когорту Z в течение промежутка времени (t 1,t] , Y , Z {S, L, D, B, D0 , B0} . Слагаемые fi (t) , i 1, 6 задают приток населения в соответствующие когорты за промежуток времени (t 1,t].
Для описания процессов инфицирования индивидуумов и спонтанной активации болезни используется схема, предложенная в работе [3]. Инфицирование происходит при одновременном посещении больными и здоровыми индивидуумами мест возможного контакта. Считается, что за единицу
времени индивидуумы из когорт
B,
B0
посещают в совокупности
(B) t
и
(B0 ) t
мест соответственно, где
(B) t
=
, =xB (t ) (B)
i =1 it
( B0 ) t
,xB0 (t ) (B0 )
i =1 it
(B) it
0,
(B0 ) it
0
− взаимно независимые, одинаково распреде-
ленные случайные величины, не зависящие от xB (t) , xB0 (t) , причем E(it(B) ) = r > 0 , E(it(B0 ) ) = k r , k (0;1) . Вероятности того, что индивидуумы из когорт S, L посетят одно из этих мест, обозначим через
S (0;1) и L (0;1) соответственно. Вероятности заражения индивидуумов после контакта с больным
индивидуумом обозначим через S (0;1) и L (0;1) соответственно. Тогда вероятности того, что индивидуум из конкретной когорты S, L будет инфицирован в течение суток, при фиксированном значении
t
задаются следующими формулами:
(S ) t
= 1 (1 S S )t( B) t(B0 )
,
(L) t
= 1 (1 L )(1 L L )t(B) t(B0 ) ,
где
L − вероятность спонтанной активации болезни для индивидуумов когорты L. Данные величины опреде-
ляют количество инфицированных на интервале (t 1,t].
Индивидуум-ориентированная модель. В популяцию вводятся когорты R и R0 , которые используются для описания индивидуумов, находящихся в состоянии ремиссии. Считается, что в когорту R поступают самоизлечившиеся индивидуумы из когорты D , а в когорту R0 − вылеченные индивидуумы
из когорты D0 . Переходы индивидуумов из D и D0 в L заменяются на переходы из D в R и из D0 в R0 . Предполагается, что больные индивидуумы после перехода в состояние ремиссии остаются инфицированными и существует риск повторного развития заболевания. Это означает, что возможны переходы индивидуумов из R в D и из R0 в D0 .
Для индивидуумов когорт D, B, D0, B0, R и R0 вводится параметр x (t) = max{0, x (t) x (t)} , t = 0,1, 2,, отражающий влияние продолжительности болезни индивидуума на вероятность его дожития
от момента времени t 1 до t . Здесь x (t) = D D,x (t) B B,x (t) D0 D0 ,x (t) B0 B0 ,x (t) − величина,
100
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
В.Н. Леоненко, К.К. Логинов
отражающая негативное влияние заболевания на продолжительность жизни индивидуума x , x (t) = R R,x (t) R0 R0 ,x (t) − величина, отражающая снижение эффекта воздействия заболевания на продолжительность жизни индивидуума x за счет его пребывания в состоянии ремиссии, D , B , D0 , B0 , R , R0 − положительные константы, D,x (t) , B,x (t) , D0 ,x (t) , B0 ,x (t) , R,x (t) , R0 ,x (t) − продолжительности пребывания индивидуума x в соответствующих когортах до момента времени t включительно. Вероятность дожития от момента времени t 1 до t для индивидуума x из когорты H
при заданной x (t) равна H ,x (t) = H ex (t1) , H {D, B, D0 , B0 , R, R0} (в отличие от базовой модели, где
H ,x (t) H 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 и R0 создаются списки, содержащие индивидуальные параметры членов когорт.
2. Процессы дожития и переходов индивидуумов когорт R и R0 в течение суток разыгрываются для каждого из этих индивидуумов последовательно путем прохода по спискам параметров. 3. Процессы дожития и переходов индивидуумов когорт D, B, D0, B0 вплоть до их гибели или выздоровления разыгрываются для каждого индивидуума с помощью функции розыгрыша судьбы отдельно взятого индивидуума. В процессе розыгрыша принадлежность индивидуума к заданной когорте в заданный момент времени регистрируется в массиве численностей индивидуумов множества T. В случае перехода рассматриваемого индивидуума в состояние ремиссии значение его параметра передается в качестве выходного параметра функции и помещается в массив параметров индивидуумов когорт R и R0 .
Формальное описание алгоритмов №№ 1, 2 для базовой модели и алгоритма № 1 для индивидуумориентированной модели приведено в работе [5].
Функция розыгрыша судьбы отдельно взятого индивидуума реализована в двух вариантах. (а) Итеративный вариант. Судьба индивидуума разыгрывается итеративно, с шагом в одни сутки (на каждом шаге с помощью метода Монте-Карло моделируется, погибнет ли индивидуум на данном шаге, и если нет, то перейдет ли он в другую когорту). (б) Вариант со скачками. На основании значений вероятностей гибели индивидуума и перехода его в другую когорту разыгрывается момент ближайшей смены статуса рассматриваемого индивидуума, после чего внутри функции происходит «скачок» во времени до этого момента с одновременным пересчетом значения индивидуального параметра.
Для сравнения быстродействия работы моделирующих программ проводился вычислительный эксперимент (табл. 1). Расчеты производились для двух наборов параметров на персональном компьютере Intel Core2 Quad 2,66 ГГц, объем оперативной памяти – 3 Гб. Отрезок моделирования – 5500 суток. Получаем, что моделирующая программа, выполненная на основе алгоритма №2 и функцией розыгрыша
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
101
ВЫЧИСЛИТЕЛЬНЫЕ АСПЕКТЫ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ...
судьбы в варианте (б), производит расчеты с первым набором параметров примерно на 70% быстрее, а со вторым набором – на 55% быстрее по сравнению с исходным вариантом алгоритма.
Заметим, что реализация розыгрыша судьбы индивидуума в виде отдельной функции также облегчает модификацию программного кода при необходимости изменения законов протекания инфекции (в перспективе возможна реализация выбора из нескольких модулей протекания инфекции и схем заражения для разных типов заболеваний).
Алгоритм №1 Алгоритм №2 (а) Алгоритм №2 (б)
Время работы на наборе № 1, с 1611,39 970,14 478,43
Время работы на наборе № 2, с 2468,01 1636,69 1111,56
Таблица 1. Быстродействие моделирующей программы с разными алгоритмами
Распределенные вычисления с использованием системы MONC
Система MONC [6] (сокращение от Monte Carlo), разработанная М.А. Марченко (Институт вычислительной математики и математической геофизики СО РАН, г. Новосибирск), является универсальной высокопроизводительной системой параллельных вычислений для методов Монте-Карло на базе сети персональных компьютеров. Она распределяет независимые копии задания по персональным компьютерам в сети, отдает команды на их исполнение, следит за ходом выполнения заданий и по завершению выполнения заданий выполняет копирование и осреднение файлов с результатами расчетов.
Был проведен вычислительный эксперимент с использованием системы MONC на комплексе из четырех компьютеров с ОС Windows XP. Исполняемым приложением была модификация моделирующей программы для индивидуум-ориентированной модели, доработанная с учетом требований, предъявляемых системой (см. [6]). Производилось нахождение оценок математических ожиданий численностей когорт по ста реализациям для двух различных наборов параметров (отрезок моделирования – 5500 суток). Сначала измерялось время выполнения расчетов на каждом из компьютеров по отдельности, затем общее количество моделируемых реализаций разделялось по компьютерам пропорционально их быстродействию. После окончания расчетов на всех машинах система MONC усредняла полученные данные и сохраняла результат.
Таким образом, в данном случае запуск модели с помощью системы MONC на комплексе из четырех компьютеров позволяет сократить время выполнения задания примерно на 66% для обоих наборов параметров по сравнению с запуском на самом быстром из имеющихся компьютеров (O94-11).
Сетевое имя
O94-31 O91-01
O94-09
Процессор
Intel Pentium 4, 3,2 ГГц AMD Athlon 64 x2 Dual Core
Processor 5200+, 2,61 ГГц Intel Pentium 4, 3 ГГц
Оперативная память
512 Мб 1 Гб
512 Мб
Время работы на наборе № 1, с 5156,7 4958,23
5826,1
Время работы на наборе № 2, с 17307,2 16755,6
19228,4
O94-11
Intel Core 2 Duo E7200, 2,53ГГц
1 Гб
3679,15
12031,9
Таблица 2. Характеристики компьютеров, задействованных для расчетов с использованием системы MONC
Сетевое имя O94-31 O91-01 O94-09 O94-11
Итого
Набор № 1 Кол-во реализаций Время выполнения расче-
тов, с 23 1216 24 1233 21 1183 32 1171
100 1233
Набор № 2 Кол-во реализа- Время выполнения расче-
ций тов, с 23 3885 24 4073 20 3850 33 4019
100 4073
Таблица 3. Результаты вычислительного эксперимента по применению системы MONC
Параллельные вычисления с использованием технологии OpenMP
Технология OpenMP (Open Multi-Processing) представляет собой набор директив компилятора, функций и переменных окружения, позволяющих создавать на базе последовательных программ многопоточные приложения, предназначенные для выполнения на многопроцессорных системах с единой па-
102
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
В.Н. Леоненко, К.К. Логинов
мятью. Важным достоинством OpenMP является то, что исходная программа не переписывается заново, а модифицируется добавлением в текст соответствующих директив [5].
К моменту написания статьи реализована поддержка технологии OpenMP в моделирующей программе для базовой модели с алгоритмом №1. Проводился вычислительный эксперимент с двумя наборами параметров, расчеты выполнялись на компьютере Intel Core 2 Duo E7400 2.8ГГц с объемом оперативной памяти 2 Гб. Отрезок моделирования – 5500 суток.
Количество потоков 1 2
Время работы на наборе № 1, с 129,968 71,578
Время работы на наборе № 2, с 191,093 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). – С. 246– 250.
4. Перцев Н.В., Пичугин Б.Ю. Индивидуум-ориентированная стохастическая модель распространения туберкулеза // Сибирский журнал индустриальной математики. – 2009. – Т. 12. – № 2(38). – С. 97– 110.
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
Научно-технический вестник Санкт-Петербургского государственного университета информационных технологий, механики и оптики, 2010, № 4(68)
103