авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ

Министерство образования и науки Российской Федерации

Федеральное государственное бюджетное образовательное учреждение высшего

профессионального образования «Нижегородский государственный

университет

им. Н.И. Лобачевского»

Факультет вычислительной математики и кибернетики

Кафедра: Математического обеспечения ЭВМ

Направление: «Прикладная математика и информатика»

ВЫПУСКНАЯ КВАЛИФИКАЦИОННАЯ РАБОТА БАКАЛАВРА Математическое моделирование динамики клеточных культур на GPU Заведующий кафедрой:

доктор физ.-мат. наук Стронгин Р.Г.

«_»2012 г.

Выполнил: студент группы 8402 Вильдеманов А.В.

Подпись Научный руководитель:

д.ф.-м.н., доцент каф. МО ЭВМ Иванченко М.В.

Подпись Нижний Новгород Содержание Содержание............................................................................................................................ Введение................................................................................................................................. Описание модели......................................................................................................... 1.

Детерминистическая модель...................................................................................... 2.

Метод Рунге-Кутты................................................................................................. 2. Стохастическая модель............................................................................................... 3.

Стохастическая химическая кинетика................................................................... 3. 3.1.1 Модификация стохастического алгоритма..................................................... Моделирование диффузии........................................................................................ 4.

Технология CUDA..................................................................................................... 5.

Параллельные алгоритмы......................................................................................... 6.

Описание реализации программной системы......................................................... 7.

Структура программной схемы............................................................................ 7. Описание основных методов программ.............................................................. 7. Результаты тестирования программной системы............................................... 7. Заключение........................................................................................................................... Список литературы.............................................................................................................. Приложение 1. Исходный код основных структурных модулей программной системы........................................................................................................................................................... Введение Одна из центральный задач современной регенеративной медицины является искусственное формирование тканей на основе клеточных культур. Например, актуальным направлением исследований является разработка комплекса «Искусственная печень» на базе культур – гепатоцитов.

Печень – жизненно важный орган, выполняющий большое количество физиологических функций, в том числе, метаболизма, синтеза, и дезинтоксикации.

Заболевания печени ослабляют её функциональные свойства, что может привести к острым состояниям, требующим интенсивной терапии, трансплантации, либо к летальному исходу.

Одним из ключевых этапов лечения пациентов с тяжелой патологией печени является замещение утраченной или недостаточной печеночной функции до восстановления ее необходимого объема в результате регенерации, либо до пересадки донорского органа. С этим связаны активные работы в области биомедицинских технологий, ведущиеся в мире.

Сложность создания эффективной системы «искусственная печень» состоит в многофункциональности человеческой печени.

Достижение целевых характеристик невозможно без создания эффективных и гибких инструментов проектирования. В первую очередь требуется разработка методов и программного обеспечения для построения математической модели эволюции клеточного материала.

С помощью математической модели возможно улучшать клинические и биохимические параметры систем основанных на клеточном взаимодействии, например, в биореакторах.

В упрощенных популяционных моделях могут возникать сложные хаотические пространственно-неоднородные решения. Кроме того, усреднение характеристик, таких как состава клеточных популяций, не позволяет учесть редкие, но принципиально значимые в случае небольших концентраций, события, что приводит к неверным результатам. Поэтому моделирование комплексной задачи роста клеточных систем изначально требует реалистичной модели.

1. Описание модели Любой живой организм состоит из множества различных видов клеток. В рассматриваемой модели нас будут интересовать три вида клеток, разделяемых по степени зрелости: стволовые, неразвитые и развитые клетки. Изменение числа клеток происходит посредством деления или гибели. В результате деления клетка может разделиться пополам или перейти из одного типа в другой. Так же для клеточных культур характерна зависимость от концентрации и от положения в пространстве. Динамику развития такой модели можно описать детерминистическими или стохастическими распределенными уравнениями.



В модель рассматриваются следующие типы клеток:

Стволовые (S) Не зрелые (T) Развитые, 8 типов (Di, i = 1…8) Поскольку вследствие укорачивания теломеразы у развитых клеток, их число делений ограничено (порядка 7), удобно разбивать их на 8 групп, идентичных по физиологическим свойствам, каждая из которых (Di, i = 1…8) отвечает клеткам прошедшим i-1 деление.

В каждый момент времени состояние каждой из популяции можно охарактеризовать ( ), где () количеством клеток. Заметим, что хотя и должен принимать целые значения, но при построении детерминистической модели рассматриваются средние значения по популяциям и поэтому величина может быть любым действительным числом, непрерывно меняющемся во времени.

Различные типы клеток различаются по времени жизни, скорости и сценариям деления. В модели возможны следующие реакции:

(симметричное деление – одна клетка делется на две себе подобных;

воспроизводство стволовых клеток) (асимметричное деление – клетка делится на одну своего вида и на одну другого;

рождение неразвитых клеток) (симметричное деление;

переход стволовых клеток в неразвитые) (симметричное деление;

переход неразвитых клеток в развитые первого типа) (симметричное деление) (гибель клетки) (гибель клетки) (гибель клетки) 2. Детерминистическая модель В качестве детерминистической модели была составлена система дифференциальных уравнений, которая учитывает все описанные выше процессы.

( ) ( ) ( ) (1) ( ) ( ) ( ) { где – отвечает за деление клетки типа S на две клетки типа S;

( ) – отвечает за деление клетки типа S на две клетки типа T в зависимости от количества клеток типа T(D) в популяции;

– отвечает за деление клетки типа S на две клетки типа S и T;

– отвечает за деление клетки типа T на две клетки типа ;

– отвечает за деление клетки типа на две клетки типа ;

– коэффициент смертности;

Данная система дифференциальных уравнений относится к стандартному классу линейных уравнений первого порядка.

Решение данной системы аналитически не представляется возможным, поэтому необходимо использовать численные методы. В качестве численного метода был взят метод Рунге-Кутты [11] четвертого порядка, явная схема с постоянным шагом.

2.1 Метод Рунге-Кутты Методы Рунге-Кутты являются семейством численных алгоритмов для решения обыкновенных дифференциальных уравнений и их систем.

Методы Рунге-Кутты является модификацией метода Эйлера, они представляют собой схемы второго и более порядка точности. На практике наиболее часто пользуются схемой четвертого порядка.

Свойства:

Метод является одноступенчатым: чтобы найти ( ) нужно знать только ( ).

( ), нужны только вычисление Не требует вычисление производных от самой функции.

Более точный, чем метод Эйлера, потому что учитывает поле ( ) не только в точке ( ), но и её окрестности.

Рассмотрим задачу Коши:

( ) { (2.1) () где и могут быть как скалярами (для уравнения) так и векторами (для систем уравнений).

Метод Рунге-Кутты четвертого порядка (явный) имеет следующий вид ( ) ( ) (2.2) ( ) ( ) ( ) { Четвертый порядок метода характеризуется тем, что при уменьшении шага в 2 раза ошибка между точным и вычисленным значением уменьшится в 16 раз.

3. Стохастическая модель Стохастическое описание является наиболее точным способом отобразить эволюцию системы, оно детально отслеживает популяционный состав после каждой реакции.

Стохастическая модель предназначена для моделирования пространственно-распределенной динамики с использованием многомерных марковских процессов рождения-гибели (алгоритм Gillespie [2]). Отличие стохастической модели от детерминистического популяционно-видового описания заключается в высоком разрешении клеточных процессов: деления рассматриваются на уровне клеток, в отличие от усредненных характеристик клеточных популяций.

3.1 Стохастическая химическая кинетика Изложим численный алгоритм моделирования стохастической динамики следуя алгоритму Gillespie [2], который изначально применялся для моделирования химических реакций.

У нас имеется:

– популяций (химических веществ).

реакций { }.

Например:.

размер популяции.

вектор изменения в – реакции.

( ).

Например:

В каждой популяции может произойти реакций.

() кинетический коэффициент реакции.

() вероятность, того что произойдет реакция в следующий бесконечно малый промежуток времени [ ), для определенного ( ).

Например:

( ). Константа Реакция заключается в том, что отвечает за интенсивность реакции,а дает вероятность того, что реакция произойдет в бесконечно малый отрезок времени. Для клеток типа вероятность реакции в () момент времени будет равна. Кинетический коэффициент реакции.

() кинетический коэффициент того, что произойдет любая реакция.

() () (3.1) Вероятность того, что любая реакция произошла в первый раз в момент времени [ ).

( ) вероятность того что до ничего не происходило.

( ) () (3.2) (3.3) ( ) () | | () (3.4) () (3.5) ] не произошла не одна Вероятность (3.5) показывает то, что в течение времени[ реакция. Вероятность первой реакции в промежуток времени [ ):

() () (3.6) Плотность вероятности возникновения первой реакции:

() () () (3.7) () (3.8) Рис1. Схематический график плотности вероятности Теперь найдем время первой реакции:

() [ ] (3.9) ( ) ( ) (3.10) () () (3.11) (3.12) () (3.13) () () () (3.14) Время первой реакции равно:

(4) Чтобы определить какая реакция произошла в момент времени необходимо [ ] и из условия:

генерировать вторую случайную величину (5) Мы находим, где номер реакции.

Схема стохастического алгоритма:

0. Инициализируем время и начальное состояние системы.

1. Для состояния и времени вычисляем все кинетические коэффициенты ( ) и их сумму ( ).

2. Генерируем и используя (4) и (5).

3. Изменяем время и состояния для следующей реакции и.

4. Записываем ( ) и переходим на Шаг 1 или заканчиваем при условии, что дошло до границы интегрирования.

Рис 2. Пример моделируемой динамики клеточных культур (стволовые клетки (S), незрелые (T) и зрелые (D) гепатоциты), Алгоритм Рунге Кутта и стохастический алгоритм.

Модификация стохастического алгоритма 3.1. Для больших популяций стохастический алгоритм сходится очень медленно, так как ( ). Чтобы увеличить скорость интегрирования шаг интегрирования при увеличении в работе Gillespie [2] был предложен, так называемый метод “первой семьи”.

“семей” { }, в каждой из которой находится Разобьем реакций на реакций { }. Для каждой семьи пересчитаем кинетические коэффициенты и () ( ). Генерируем время найдем их сумму, через которое произойдет следующая реакция и пару ( ), которая идентифицирует произошедшую реакцию. Для этого потребуется случайных чисел, распределенных по нормальному закону и принадлежащих отрезку [ ]:

)( ) ( (6) () где мы берем { } { (7) { } В итоге находим () () (8) Зная реакцию и ее время, выполняем Шаг 3 и Шаг 4 стохастического алгоритма.





Данная модификация существенно увеличивает скорость интегрирования (рис 3), так как на ( ), а только те, которые затронула каждом шаге мы не пересчитываем все коэффициенты реакция. При больших размерностях необходимо брать большее количество семей, чтобы получить наилучшее ускорение (рис 3).

L= 6 L= L= 8 16 32 64 128 256 512 Размерность Рис 3. График отношения «семейного» к стохастическому алгоритмов.

(Windows 7 Enterprise 64-bit, CPU Intel Core i5 M540 @ 2.53GHz RAM 4.00 GB) Из графика (рис. 3) видно, что оптимальная размерность семьи равняется 8, так как при размерности пространства равным 16, разбиение на предложенное количество семей дало одинаковое ускорение. Алгоритм работает одинаковое количество времени при разбиение на 4 и 8 семей при размерности 32, а при разбиение на 2 семьи в два раза меньше. По индукции получаем, что оптимальное разбиение будет равняться отношению размерности пространства к оптимальной размерности семьи, то есть к 8.

4. Моделирование диффузии Чтобы учесть влияние соседних популяций друг на друга, введем диффузионное взаимодействие между ними. Это взаимодействие возникает в следствие диффузии сложных молекулярных и питательных веществ (цитоминов, хемокинов) в межклеточном пространстве и самих клеток. Она является выравнивающим фактором, модифицируя кинетические коэффициенты деления/смерти. Диффузионное взаимодействие в нашей модели полагается локальным.

Для одномерной модели схема взаимодействия и математическая модель выглядят следующим образом:

( ) ( ) (9) Для двумерной модели:

( ) ( ) ( ) (10) Здесь коэффициенты диффузии по соответственно.

Заметим, что в эксперименте может реализоваться и противоположная ситуация, когда присутствует сильное перемешивание среды. Тогда связь приобретает глобальный характер и в качестве математической модели следует использовать:

( ) (11) 5. Технология CUDA CUDA – это архитектура параллельных вычислений от компании NVIDIA, позволяющая увеличить вычислительную производительность благодаря использованию графических процессоров (GPU). Разработчики программного обеспечения, ученые и исследователи используют CUDA в таких областях, как обработка видео, астрофизика, вычислительная биология и химия, моделирование, томография, трассировка лучей и многое другое. Все эти задачи хорошо ложатся на GPU, потому что в этих задачах большой поток данных многократно обрабатывается одними и теме же командами, то есть представляют собой SIMD-модель (Single Instruction Multiple Data). Компания NVidia использует термин SIMT (Single Instruction, Multiple Thread).

CUDA состоит из концепции, что GPU выступает в роли массивно-параллельного сопроцессора к CPU. Программа на CUDA выполняется как на CPU (последовательный код) так и на GPU (для массивно-параллельные вычисления), как набор одновременно выполняющихся нитей(потоков). Нити разбиваются на группы по 32 нити, называемые warp`ами. Только нити в пределах одного warp`а выполняются физически одновременно.

Нити разных warp`ов могут находиться на разных стадиях выполнения программы, при этом управление warp`ами осуществляет сам GPU.

Все запущенные на выполнение нити организованы в следующую иерархию. Верхний уровень иерархии – сетка (grid) – соответствует всем нитям, выполняющим данное ядро.

Верхний уровень представляет из себя одномерный или двумерный массив блоков (block).

Каждый блок это одномерный, двумерный или трехмерный массив нитей. При этом все блоки одной сетки имеют одинаковый размер и размерность.

Графические процессоры являются параллельными архитектурами, которые обладают своей памятью (DRAM), содержат ряд потоковых мультипроцессоров (SM, streaming multiprocessor), каждый из которых способен выполнять до 1024 нитей одновременно.

Каждый мультипроцессор работает независимо от остальных. Когда CUDA запускает ядро на исполнение, то блоки сетки выполняются на имеющихся мультипроцессорах. При этом каждый блок целиком выполняется на одном из потоковых мультипроцессоров.

Мультипроцессор способен одновременно выполнить до восьми блоков.

Чисто физически память GPU можно разделить на DRAM и на память, размещенную в потоковых мультипроцессорах. Однако классификация памяти не ограничивается её чисто физическим расположением. Каждый потоковый мультипроцессор содержит 8192 или 32-битовых регистра. Имеющиеся регистры разделяются между нитями блока на этапе компиляции. Регистры обладают максимальной скоростью доступа, но нить не имеет доступа к регистрам других нитей. Если регистров не хватает то для размещения локальных данных нити используют локальную память, размещенная в DRAM. Доступ к локальной памяти характеризуется высокой латентностью – от 400 до 600 тактов. Разделяемая (shared) память – память расположена непосредственно в потоковом мультипроцессоре. Она выделяется на уровне блоков – каждый блок получает в свое распоряжение одно и тоже количество разделяемой памяти. На каждом мультипроцессоре содержится 16Кбайт разделяемой памяти. Глобальная память – это обычная DRAM-память, которая доступна всем нитям.

6. Параллельные алгоритмы Как численный метод интегрирования Рунге-Кутты, так и стохастический алгоритм выполняют одни и те же операции над разными данными, что удовлетворяет критерию SIMD (single instruction, multiple data). Поэтому параллельная реализация была разработана под SIMD-процессоры, а именно под графические ускорители NVIDIA.

При примитивном распараллеливании алгоритма Рунге-Кутты под с GPU использованием технологии CUDA было получено ускорение в 11 раз по сравнению с последовательной версией на центральном процессоре (Рис 4).

Рис 4. Тестирование на OS Windows 7 Enterprise 64-bit, Intel Xeon E5520 2.27GHz ( CPUs), nVidia Tesla C Грамотное использование особенностей архитектуры GPU, а именно разделяемой памяти, доступ к которой в 100 раз быстрее, чем к общей, позволило получить ускорение в 45 раз. Отмеченным недостатком использования этого типа памяти является то, что к ней могут обращаться только нити одного блока (алгоритм делится на блоки, которые выполняются независимо друг от друга на потоковых мультипроцессорах на GPU, блоки же делятся на нити, а нити уже параллельно выполняют операции). Также отметим, что технология CUDA предоставляет математический набор функций пониженной точности, но обеспечивают большее быстродействие. Ускорение, полученное с использованием этих функций, составило 50 раз по отношению к одному CPU (рис 4).

При реализации диффузии мы имеем зависимость по данным, и чтобы потоки брали нужную информацию на границах блоков, а не измененную (как известно блоки могут исполнятся с разной скоростью и разделяемая память не доступна другим блокам) реализовано чтение из глобальной памяти, а запись в разделяемую, тем самым мы обходим проблему чтения не верных данных и получаем оптимальный по памяти алгоритм (рис 5).

Рис 5. Схема работы с памятью в задаче моделирования диффузии.

Так как система ОДУ со временем приходит к своему состоянию равновесию и слабо реагирует на погрешности счета, то можно воспользоваться математическими CUDA функциями одинарной точности, но меньшим временем выполнения.

7. Описание реализации программной системы 7.1 Структура программной схемы Для реализации поставленной задачи была разработана программная система на алгоритмическом языке C++ с использованием алгоритмического языка параллельного программирования на GPU – CUDA C. Разработка приложений производилась в визуальной среде разработки Microsoft Visual Studio 2008. Было написано около 2000 строк кода на С++ и 1200 строк кода на CUDA C. Для обработки данных и получения изображений (видео) использовался скриптовый язык Python.

7.2 Описание основных методов программ Программа представляет собой набор методов, позволяющих получать требуемый результат. Кроме того, все основные алгоритмы могут выполняться как на центральном процессоре, так и на графическом (рис 6). Методы, использующие вызовы вычислений на GPU, главным образом, запускают эти вычисления, обрабатывают и сохраняют полученные результаты в файл.

Чтение данных Алгоритм Рунге-Кутты CPU + Стохастический алгоритм GPU Обработка CPU данных Рис 6. Схема работы и взаимодействия программной системы.

Приложение считывает из файла информацию о моделирование: размер системы, размер шага интегрирования, время моделирования, коэффициенты диффузии системы, количество «семей» (для семейного алгоритма), имя файла с коэффициентами системы.

Методы реализованные в программе:

1) void parserCommandLine( int argc, char * argv[] ) – считывает параметры с командной строки, задает параметры системы из файла, который был подан на вход.

2) double Runge(float* x0) – метод Рунге-Кутта на CPU, принимает на вход начальное состояние системы и результаты записывает в файл или на экран (опционо) 3) double stochasticAlgorithm(float* x0 ) – алгоритм Gillespie на CPU.

4) void familyStochasticAlgorithm(float *x0) – «семейный» алгоритм Gillespie на CPU.

5) double GPU_Runge(float* x0) – метод Рунге-Кутта на GPU, в данном методе используется CUDA ядро: global void calculate(float* x_,float* k1,float* k2,float* для одномерного случая и k3,float* k4, int size, int n, int blockSize, float d, float step) – global void calculate_2D(float* x_,float* k1,float* k2,float* k3,float* k4, int size, int n, int blockSize, float d, float step) – для двумерной области.

6) double GPU_stochasticAlgorithm( float *x0) – стохастический алгоритм на GPU, в данном методе используется CUDA ядра: global void diffusion(float *x, float *y, int – алгоритм рассчета диффузии, global void size, float d, float sum) – метод пересчета кинетических propensityFunction(float *a, float *y_, int n, int m) коэффициентов. Ядра имеются для одномерного и двумерного случая.

7) void GPU_familyAlgorithm(float *x0) – семейный стохастический алгоритм на GPU, в данном методе используется CUDA ядра аналогичные CUDA ядрам в стахостическом алгоритме.

8) float cuReduce ( int n, float *in_data, float *out_data ) – метод редуцирования массива кинетических коэффициентов для стохастических алгоритмов на GPU.

Использует CUDA ядро: template unsigned int blockSizeglobal void reduce4(float *g_idata, float *g_odata, unsigned int n) 7.3 Результаты тестирования программной системы Проверка корректности математической модели проходила путем сравнения с физиологическими данными (Hussain et al. 2005, Hoeme et al. 2010). Характерным сценарием развития клеточных культур на плоской подложке, наблюдаемый в биологическом эксперименте, является интенсивное деление стволовых клеток, дифференцировка в незрелые гепатоциты, а затем – в зрелые. При этом достигается тысячекратное увеличение числа клеток. Характерный временной масштаб этой стадии составляет несколько часов.

Затем происходит насыщение роста зрелых гепатоцитов, как за счет конкуренции взрослых клеток за питательные ресурсы, так и за счет исчерпания пула стволовых клеток. На характерном времени порядка суток численность взрослых гепатоцитов начинает сокращаться. Начальные условия, использованные в численном эксперименте, соответствовали физиологически реалистичным: в начале эксперимента в «пробирке»

находится около 106 клеток, из них 60% стволовых клеток и 40% неразвитых. Полученные данные: происходит деление и дифференцировка с образованием порядка 4*106 развитых клеток, после определенного времени (порядка суток, единиц времени математической модели) клетки начинают погибать (рис 7).

Рис 7. Пример моделируемой динамики клеточных культур (стволовые клетки (S), незрелые (T) и зрелые (D) гепатоциты).

Полученные численные результаты качественно согласуются с экспериментальными наблюдениями, что подтверждает корректность используемой модели.

Заключение Разработаны методы математического моделирования клеточных культур, состоящих из клеток различной степени зрелости (стволовые, незрелые, зрелые), реализующие эти методы с помощью средств последовательного и параллельного программирования.

Математическая модель детализирует процессы деления и дифференцировки клеток различных типов зрелости, содержит как реализацию в приближении непрерывного времени, так и реализацию на уровне единичных клеточных событий.

Результатом проделанной работы является программный комплекс позволяющий получить картину динамики развития клеточных культур как на уровне отдельных молекул, так и на общем уровне. Комплекс позволяет быстро моделировать динамику эксперимента, который может в реальности проходить несколько суток. Это позволит прогнозировать результаты и существенно улучшить планирование эксперимента. ПО было реализовано на алгоритмических языках C++ и CUDA. Суммарное число строк кода около 2000, из них строк на CUDA и около 500 строк вспомогательных функций.

Работа выполнена при поддержке компании «Т-платформы» и ФЦП Кадры, контракт 02.740.11. Список литературы 1. Franziska Michor, Martin A. Nowak, Steven A. Frank and Yoh Iwasa Stochastic elimination of cancer cells [Журнал] // The Royal Society 2003 г.. – 270. – стр. 2017- 2. Gillespie Daniel T. Stochastic Simulation of Chemical Kinetics [Журнал] // Annm. Rev. Phys.

Chem. 2007 г.. – 58. – стр. 35- 3. Green J.E.F., Waters S.L., Shakesheff K.M., Byrne H.M. A Mathematical Model of Liver Cell Aggregation In Vitro [Журнал] // Bulletin of Mathematical Biology 2009 г.. – 71. – стр.906– 4. Hoehme Stefan, Drasdo Dirk A cell-based simulation software for multi-cellular systems [Журнал] // bioinformatics 2010 г.. - 26. – стр. 2641– 5. Kirouac Daniel C, Madlambayan Gerard J, Mei Yu1, Sykes Edward A, Caryn Ito3 and Peter W Zandstra Cell–cell interaction networks regulate blood stemand progenitor cell fate [Журнал] // Molecular Systems Biology 2009г.. – 293. – стр. 1- 6. O`Dea R. D., Waters S. L., Byrne H. M. A multiphase model for tissue construct growth in a perfusion bioreactor [Журнал] // Mathematical Medicine and Biology 2010 г.. – 27 – стр.

7. Боресков А.В., Харламов А.А. Основы работы с технологией CUDA [Книга]. – М.: ДМК Пресс, 2010. – 232 с.: ил. ISBN 978-5-94074-578- 8. Вильдеманов А.В., Иванченко М.В. Моделирование пролиферации в ансамблях живых клеток [Публикация] // Молодежная Радиофизическая конференция 2011г.. - стр 85-87.

9. Вильдеманов А.В., Иванченко М.В. Моделирование размножения клеточных культур [Публикация] // Форум HPC-2011 2011 г.. - стр. 51- 10. Кузенков О.А., Рябова Е.А., Круподерова К.Р. Математические модели процессов отбора. Учебное пособие [Книга]. Нижний Новгород: Нижегородский госуниверситет, 2010. 133 с.

11. Метод Рунге – Кутты [В Интернете].- ru.wikipedia.org/wiki/Метод_Рунге_-_Кутты.

Приложение 1. Исходный код основных структурных модулей программной системы Метод Рунге-Кутта:

// размер поверхности 1 extern int N ;

extern double step ;

extern float finish ;

// печать extern bool print ;

extern bool save ;

// коэффициент диффузии extern float d ;

// коэффициент диффузии в глубь extern float dd ;

// размер системы extern int size ;

extern float bst, btd, bdd, btt1, btt2, bss, dead1, dead2, dead3 ;

extern int rang ;

// количество реакций extern int m ;

void F(float* x,float* k) { float *y = new float[size];

#ifdef ENABLE_2D diffusion_2d(x,y,d,dd,size,rang,N);

#else diffusion(x,y,d,size,rang);

#endif for ( int i = 0;

i size;

i += rang ) { k[i+0] = x[i+0] * (-btt1/(y[i+1]+1.)-btt2/(y[i+2]+1.)-dead1*y[i])+bss*y[i]/(y[i]+1.);

k[i+1] = x[i+1] * (-btd-dead2*y[i+1])+y[i+0]*bst+2.*y[i]*(btt1/(y[i+1]+1.)+btt2/(y[i+2]+1.));

k[i+2] = x[i+2] * ( -dead3*y[i+2]-bdd)+2.*btd*y[i+1];

for ( int j = 1;

j 7;

j++ ) k[i+2+j] = x[i+2+j] * ( -dead3*y[i+2+j] - bdd) + 2.*bdd * y[i+1+j];

k[i+2+7] = x[i+2+7] * ( -dead3*y[i+2+7])+2.*bdd*y[i+1+7];

} delete[] y;

} double Runge(float* x0) { printf("Runge-Kutta Algorithm\n\n");

printf("Size is %d\n",N);

printf("final time is %f\n\n",finish);

float t=0;

// current time int stepSave = 10;

// для шага слива в фаил int temp = 0;

float eps = 0.00000001;

int j = 0;

// position of population int iter = 0;

float* x = new float [size];

// array of variables for ( int i = 0;

i size;

i++) { x[i] = x0[i];

} float* k1 = new float [size];

float* k2 = new float [size];

float* k3 = new float [size];

float* k4 = new float [size];

float* xx = new float [size];

clock_t stime;

clock_t ftime;

stime = clock();

while ( t = finish ) { t = t + step;

iter++;

if ( save ) { { Save(t, x, size);

} } if( print ) MyPrint( t, x, size );

//--------------------------------------- // расчитать k F ( x, k1 );

for ( int i = 0;

i size;

i++ ) { xx[i] = x[i] + k1[i] * (step / 2.0);

} //--------------------------------------- // расчитать k F( xx, k2 );

for ( int i = 0;

i size;

i++ ) { xx[i] = x[i] + k2[i] * (step / 2.0);

} //--------------------------------------- // расчитать k F( xx, k3 );

for ( int i = 0;

i size;

i++ ) { xx[i] = x[i] + k3[i] * step;

} //--------------------------------------- // расчитать k F( xx, k4 );

for ( int i = 0;

i size;

i++ ) { x[i] = x[i] + step / 6.0 * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]);

// расчитать решение на новом шаге //--------------------------------------- } } ftime = clock();

double t_time = (double)(ftime - stime)/CLOCKS_PER_SEC;

printf("\n\ntime = %f\n",t_time );

double time_iter = (t_time*1000)/(double)iter;

printf("iter = %d\ntime iter = %f msec\n\n\n", iter, time_iter);

delete[] k1;

delete[] k2;

delete[] k3;

delete[] k4;

delete[] xx;

delete[] x;

return t_time;

} Стохастический Метод:

double stochasticAlgorithm(float* x0 ) { printf("Stochastic Algorithm\n\n");

printf("Size is %d\n",N);

printf("final time is %f\n\n",finish);

// счетчик итераций int iter = 0;

clock_t stime;

clock_t ftime;

double t = 0;

// шаг по времени float h;

// обший кинетический коэффициент float a0 = 0;

// локальный помер реакции int num = -1;

// координата в пространстве int coord = -1;

// номер в массиве реакций, где произошла реакция int number;

#ifdef ENABLE_2D int countReaction = m*N*N;

#else int countReaction = m*N;

#endif // вектор изменений float * v = new float[rang];

float * a = new float[countReaction];

// propensity function // вектор переменных float * x = new float[size];

float * y = new float[size];

for ( int i = 0;

i size;

i++ ) { x[i] = x0[i];

} stime = clock();

//--------------------------------------- int ff = 5;

double stepSave = 0.0;

while ( t = finish ) { if ( print ) { MyPrint(t,x,size);

} if ( save ) { if ( (stepSave - t)*(stepSave - t) EPS ) { stepSave += step;

Save(t, x, size);

} } #ifdef ENABLE_2D diffusion_2d(x,y,d,dd,size,rang,N);

#else diffusion(x,y,d,size,rang);

#endif changePropensityFunction(a,y,size);

//p a0 = 0;

for(int i = 0;

i countReaction;

i++) a0 += a[i];

number = numberReaction( a, countReaction, a0 );

coord = number/m;

num = number - m*coord;

changeVector(v, num);

if ( num -1 ) { int k = rang*coord;

for(int i = k;

i k+rang;

i++) { x[i] += v[i-k];

if(x[i] 0.0001) x[i] = 0;

} } //-------------------------------------- h = timeFirstReaction( a0 );

t = t + h;

iter++;

if ( iter % 1000000 == 0 ) printf("iter = %f\n",t );

}// end while ftime = clock();

double t_time = (double)(ftime - stime)/CLOCKS_PER_SEC;

printf("\n\ntime = %f\n",t_time );

double time_iter = (t_time*1000)/(double)iter;

printf("iter = %d\ntime iter = %f msec\n\n\n", iter, time_iter);

delete[] a;

delete[] v;

delete[] x;

delete[] y;

return t_time;

} «Семейный» алгоритм:

void familyStochasticAlgorithm(float *x0) { printf("Stochastic Family Algorithm\n\n");

printf("Size is %d\n",N);

printf("Family is %d\n\n",L);

double stime;

double ftime;

int iter = 0;

double t = 0;

double minTime;

// номер семьи int l = -1;

// локальный помер реакции int num = -1;

// координата в пространстве int coord = -1;

// номер в массиве реакций, где произошла реакция int number = -1;

double r;

double sum;

#ifdef ENABLE_2D int numberReaction = m*N*N;

float * a0_l = new float[L*L];

// family propensity function double * time = new double[L*L];

#else int numberReaction = m*N;

float * a0_l = new float[L];

// family propensity function double * time = new double[L];

#endif int sizeFamily = N*m/L;

// вектор изменений float *v = new float [rang];

float *a = new float[numberReaction];

// propensity function // вектор переменных float *x = new float[size];

float *y = new float[size];

for ( int i = 0;

i size;

i++ ) { x[i] = x0[i];

} double stepSave = 0.0;

stime = clock();

//--------------------------------------- while( t = finish) { if( print) { MyPrint(t,x,size);

} if ( save ) { if ( (stepSave - t)*(stepSave - t) EPS ) { stepSave += step;

Save(t, x, size);

} } #ifdef ENABLE_2D if ( iter == 24939) printf("");

if ( l != -1 ) { diffusion_2d(x,y,d,dd,size,rang,N);

// diffusion_2d(x,y,coord);

changePropensityFunction_2d(a,y,coord);

// changePropensityFunction(a,y,size);

} else { diffusion_2d(x,y,d,dd,size,rang,N);

changePropensityFunction(a,y,size);

} // ---------------------------------------- // // (xT, yT) // | | || // | | || // | | || // ----------------------- // | | || // | | || // | | || // ----------------------- // | | || // | | || // | | || // ------------------------ (xB, yB) int xTop = 0, yTop = 0, xBottom = L, yBottom = L;

if( l != -1 ) { xTop = l%L;

yTop = l/L;

xBottom = xTop + 1;

yBottom = yTop + 1;

if ( coord != 0 ) { if ( (coord+1) % (N/L) == 0 && (1-(coord+1)%N == 0) ){ xBottom++;

} if ( coord % (N/L) == 0 && (1-(coord%N == 0 ))){ xTop--;

} if ( (coord/N) % (N/L) == 0 && (1- coord/N == 0) ){ yTop--;

} if ( (coord/N + 1) % (N/L) == 0 && (1- (coord/N +1 == N ) )){ yBottom++;

} } } for ( int i = yTop;

i yBottom;

i++) for(int j = xTop;

j xBottom;

j++) { int index = i*L+j;

if ( index = L*L ) printf("error access (a0_l[%d])\n",index);

a0_l[index] = 0;

int startIndexFamily = i*sizeFamily*N + j*sizeFamily;

for ( int iL = 0;

iL N/L;

iL++) for ( int jL = startIndexFamily;

jL startIndexFamily + sizeFamily;

jL++) a0_l[index] += a[iL*N*m+jL];

} for(int i = 0;

iL*L;

i++) { r = MyRand();

if(a0_l[i]) time[i] = -log(r)/(a0_l[i]);

} l = 0;

minTime = time[0];

for(int i = 1;

i L*L;

i++) { if(minTime time[i]) { minTime = time[i];

l = i;

} } int index = (l/L)*sizeFamily*N + (l%L)*sizeFamily;

//l*sizeFamily;

sum = 0;

num = -1;

r = MyRand();

bool flag = true;

for(int i = 0;

i N/L;

i++) { if (flag) { for(int j = index;

j index+sizeFamily;

j++) { sum += a[i*N*m+j];

number = i*N*m+j;

if(a0_l[l] * r sum ) { flag = false;

break;

} } } else { break;

} } coord = number/ m;

num = number - m*coord;

changeVector(v, num);

//-------------------------------------- t = t + minTime;

iter++;

if(num -1 ) { int k = rang*coord;

for(int i = k;

i k+ rang;

i++) { x[i] += v[i-k];

if(x[i] 0.0001) x[i] = 0;

} } #else if ( l != -1 ) { diffusion(x,y,d,size,N,rang,coord);

changePropensityFunction(a,y,N,coord);

} else { diffusion(x,y,d,size,rang);

changePropensityFunction(a,y,size);

} int startFamily = 0, endFamily = L;

if( l != -1 ) { startFamily = ( m*coord)/(sizeFamily);

endFamily = startFamily +1;

if ( m*coord % sizeFamily == 0 ) { startFamily = ( m*(coord-1))/(sizeFamily);

endFamily = ( m*coord)/(sizeFamily) + 1;

} if ( m*(coord+1) % sizeFamily == 0 ) { startFamily = ( m*coord)/(sizeFamily);

endFamily = ( m*(coord+1))/(sizeFamily)+1;

} if ( coord == 0 ) { startFamily = 0;

endFamily = 1;

} if ( coord == N-1 ) { startFamily = L-1;

endFamily = L;

} } for(int i = startFamily;

i endFamily;

i++) { a0_l[i] = 0;

} for ( int k = startFamily;

k endFamily;

k++) { for ( int i = k*sizeFamily;

i sizeFamily*(k+1);

i++) a0_l[k] += a[i];

} for(int i = 0;

iL;

i++) { r = MyRand();

if(a0_l[i]) time[i] = -log(r)/(a0_l[i]);

} l = 0;

minTime = time[0];

for(int i = 1;

i L;

i++) { if(minTime time[i]) { minTime = time[i];

l = i;

} } int index = l*sizeFamily;

sum = 0;

num = -1;

r = MyRand();

for(int i = index;

i sizeFamily*(l+1);

i++) { sum += a[i];

if(a0_l[l] * r sum ) { number = i;

break;

} } coord = number/ m;

num = number - m*coord;

changeVector(v, num);

//-------------------------------------- t = t + minTime;

iter++;

if(num -1 ) { int k = rang*coord;

for(int i = k;

i k+ rang;

i++) { x[i] += v[i-k];

if(x[i] 0.0001) x[i] = 0;

} } #endif } ftime = clock();

double t_time = (double)(ftime - stime)/CLOCKS_PER_SEC;

printf("\n\ntime = %f\n",t_time );

double time_iter = (t_time*1000)/(double)iter;

printf("iter = %d\ntime iter = %f msec\n\n\n", iter, time_iter);

delete[] a;

delete[] a0_l;

delete[] time;

delete[] v;

delete[] x;

delete[] y;

} Метод Рунге-Кутта в виде функций CUDA:

global void calculate(float* x_,float* k1,float* k2,float* k3,float* k4, int size, int rang, int blockSize, float d, float step) { extern shared float x[];

extern shared float y[];

extern shared float k[];

unsigned int tid = threadIdx.x;

unsigned int i = blockDim.x*blockIdx.x +threadIdx.x;

if ( i size ) { x[tid] = x_[i];

syncthreads();

diff(x,x_,(y + blockDim.x),i,tid,d,rang,size);

syncthreads();

int l = (blockDim.x*blockIdx.x)/rang+tid;

if( tid blockDim.x/rang) { if ( blockIdx.x*blockDim.x+tid*rang size ) coefficient(k1,(x),(y + blockDim.x),tid,l,rang);

} syncthreads();

float b = fdividef( step, 2.0f);

k[tid + 2*blockDim.x] = x[tid] + b*k1[i];

//-------------------------------------------- syncthreads();

diff1(x_,(y + blockDim.x),(k+ 2*blockDim.x),k1,i,tid,d,rang,size, b);

syncthreads();

if( tid blockDim.x/rang) { if ( blockIdx.x*blockDim.x+tid*rang size ) coefficient(k2,(x),(y + blockDim.x),tid,l,rang);

} syncthreads();

k[tid+ 2*blockDim.x] = x[tid] + b*k2[i];

//-------------------------------------------- syncthreads();

diff1(x_,(y + blockDim.x),(k+ 2*blockDim.x),k2,i,tid,d,rang,size, b);

syncthreads();

if( tid blockDim.x/rang) { if ( blockIdx.x*blockDim.x+tid*rang size ) coefficient(k3,(x),(y + blockDim.x),tid,l,rang);

} syncthreads();

b = step;

k[tid+ 2*blockDim.x] = x[tid] + b*k3[i];

//-------------------------------------------- syncthreads();

diff1(x_,(y + blockDim.x),(k+ 2*blockDim.x),k3,i,tid,d,rang,size, b);

syncthreads();

if( tid blockDim.x/rang) { if ( blockIdx.x*blockDim.x+tid*rang size ) coefficient(k4,(x),(y + blockDim.x),tid,l,rang);

} syncthreads();

x_[i] += step / 6.0f * (k1[i] + 2.0f * k2[i] + 2.0f * k3[i] + k4[i]);

} } Стохастический алгоритм на CUDA:

int threads = (N 512) ? N : 512;

cudaEvent_t s,f;

cudaEventCreate ( &s );

cudaEventCreate ( &f );

cudaEventRecord ( s, 0 );

while ( t = finish ) { int blocks = size/threads;

int sizeShared = threads*sizeof(float);

GPU_diffusion blocks,threads,sizeShared(x_gpu,y_gpu,d,rang,size);

blocks = N/threads;

sizeShared = rang*threads*sizeof(float);

propensityFunction blocks,threads,sizeShared(a_gpu,x_gpu,rang,m);

cudaMemcpy ( a,a_gpu,countReaction*sizeof( float ), cudaMemcpyDeviceToHost );

a0 = 0;

a0 = cuReduce(countReaction, a_gpu, sums);

r = MyRand();

number = numberReaction( a, countReaction, a0 );

coord = number/m;

num = number - m*coord;

changeVector(v, num);

cudaMemcpy ( v_gpu, v, rang * sizeof ( float ), cudaMemcpyHostToDevice );

////-------------------------------------- if ( print ) { cudaMemcpy ( x, x_gpu, size * sizeof ( float ), cudaMemcpyDeviceToHost );

MyPrint(t,x,size);

} if ( save ) { if ( (stepSave - t)*(stepSave - t) EPS ) { cudaMemcpy ( x, x_gpu, size * sizeof ( float ), cudaMemcpyDeviceToHost );

stepSave += step;

Save(t, x, size);

} } h = timeFirstReaction( a0 );

//--------------------------------------------------------------------- iter++;

blocks = size/threads;

final1,1(x_gpu,rang,v_gpu,coord);

t = t + h;

} «Семейный» алгоритм на CUDA:

//--------------------------------------- while( t = finish) { int blocks = size/threads;

int sizeShared = threads*sizeof(float);

cudaMemcpy ( x_gpu, x, size * sizeof ( float ), cudaMemcpyHostToDevice );

GPU_diffusion blocks,threads,sizeShared(x_gpu,y_gpu,d,rang,size);

blocks = N/threads;

sizeShared = rang*threads*sizeof(float);

propensityFunction blocks,threads,sizeShared(a_gpu,y_gpu,rang,m);

sizeShared = (2*L+rang)*sizeof(float);

cudaMemcpy ( a,a_gpu,countReaction*sizeof( float ), cudaMemcpyDeviceToHost );

for ( int i = 0;

i L;

i++ ) { a0_l[i] = 0;

for ( int j = i*sizeFamily;

j sizeFamily*(i+1);

j++) a0_l[i] += a[j];

r = MyRand();

time[i] = -log(r)/(a0_l[i]);

} float maxTime;

int l;

maxTime = time[0];

for( int j = 1;

j L;

j++ ) { if ( maxTime time[j] ) { maxTime = time[j];

l = j;

} } for ( int i = 0;

i L;

i++ ) while( maxTime = time[i] ) { float sum = 0;

int number = -1;

r = MyRand();

for ( int j = i*sizeFamily;

j sizeFamily*(i+1);

j++) { sum += a[j];

number = j;

if(a0_l[i] * r sum ) { break;

} } if(number -1 ) { int coord = number/m;

int num = number - m*coord;

changeVector(v, num);

int k = rang*coord;

for(int j = k;

j k+rang;

j++) { x[j] += v[j-k];

if(x[j] 0.0001) x[j] = 0;

} } changePropensityFunction(a,x,N,coord);

for ( int j = i*sizeFamily;

j sizeFamily*(i+1);

j++) a0_l[i] += a[j];

r = MyRand();

time[i] += -log(r)/(a0_l[i]);

} if( print) { cudaMemcpy ( x, x_gpu, size * sizeof ( float ), cudaMemcpyDeviceToHost );

MyPrint(t,x,size);

} if ( save ) { if ( (stepSave - t)*(stepSave - t) 0.0001) { stepSave += step;

Save(t, x, size);

} } t += maxTime;

} Файл ресурсов:

void diffusion(float * x,float * y,float d, int n, int rang ) { for ( int i = 0;

i rang;

i++) { y[i] =(1-d/2.)*x[i]+0.5*d*x[i+rang];

} //--------------------------------------- for ( int i = rang;

i n-rang;

i++) { y[i] = (1-d)*x[i]+0.5*d*(x[i-rang] + x[i+rang]);

} //--------------------------------------- for ( int i = n-rang;

i n;

i++) { y[i] = (1-d/2.)*x[i]+0.5*d*(x[i-rang]);

} } void diffusion_2d(const float * x,float * y, float d, float dd, int n, int rang, int N ) { for ( int i = 0;

i rang;

i++) // левый верхний { y[i] = ( 1 - d/2.f)*(1 - dd/2.f ) * x[i] + 0.5*d*x[i+rang] + 0.5*dd*x[N*rang+i];

} //--------------------------------------- for ( int i = rang*N-rang;

i rang*N;

i++) // правый верхний { y[i] = ( 1 - d/2.f)*(1 - dd/2.f ) * x[i] + 0.5*d*x[i-rang] + 0.5*dd*x[i+N*rang];

// } //--------------------------------------- for ( int i = rang;

i N*rang-rang;

i++) // верх { y[i] = ( 1 - d)*(1 - dd/2.f ) * x[i] + 0.5*d*x[i-rang] + 0.5*d*x[i+rang] + 0.5*dd*x[N*rang+i];

} //--------------------------------------- for ( int i = n-rang;

i n;

i++) // правый нижний { y[i] = ( 1 - d/2.f)*(1 - dd/2.f ) * x[i] + 0.5*d*x[i-rang] + 0.5*dd*x[i-N*rang];

} //--------------------------------------- for ( int i = n - N*rang;

i n-N*rang + rang;

i++) // левый нижний { y[i] = ( 1 - d/2.f)*(1 - dd/2.f ) * x[i] + 0.5*d*x[i+rang] + 0.5*dd*x[i-N*rang];

} //--------------------------------------- for ( int i = n - N*rang+rang;

i n - rang;

i++) // низ { y[i] = ( 1 - d)*(1 - dd/2.f ) * x[i] + 0.5*d*x[i+rang] + 0.5*d*x[i-rang] + 0.5*dd*x[i-N*rang];

} //--------------------------------------- for ( int j = N*rang;

j n-N*rang ;

j+= N*rang) // левый for ( int i = j;

i j+rang;

i++ ) { y[i] = ( 1 - d/2.)*(1 - dd ) * x[i] + 0.5*d*x[i+rang] + 0.5*dd*x[N*rang+i]+ 0.5*dd*x[i-N*rang];

} //--------------------------------------- for ( int j = 2*N*rang - rang;

j n - N*rang ;

j+= N*rang) // правый for ( int i = j;

i j+rang;

i++ ) { y[i] = ( 1 - d/2.)*(1 - dd ) * x[i] + 0.5*d*x[i-rang] + 0.5*dd*x[N*rang+i]+ 0.5*dd*x[i-N*rang];

} //--------------------------------------- for ( int j = 1;

j N -1;

j++) // мякоть for ( int i = j*N*rang+rang;

i (j+1)*N*rang-rang;

i++) { y[i] = ( 1 - d)*(1 - dd ) * x[i] + 0.5*d*x[i-rang] + 0.5*d*x[i+rang] + 0.5*dd*x[N*rang+i]+ 0.5*dd*x[i-N*rang];

} } void PropensityFunction(float * a, const float * y, int ind, int j ) { a[0+ind] = bst*y[j];

// S - S + T a[1+ind] = (btt1*y[j]/(y[1+j]+1)+ btt2*y[j]/(y[2+j]+1));

// S - T + T a[2+ind] = bss*y[j]/(y[j]+1);

// S - S + S a[3+ind] = y[1+j]*btd;

// T - D1 + D a[4+ind] = bdd*y[2+j];

// D1 - D2 + D a[5+ind] = bdd*y[3+j];

// D2 - D3 + D a[6+ind] = bdd*y[4+j];

// D3 - D4 + D a[7+ind] = bdd*y[5+j];

// D4 - D5 + D a[8+ind] = bdd*y[6+j];

// D5 - D6 + D a[9+ind] = bdd*y[7+j];

// D6 - D7 + D a[10+ind] = bdd*y[8+j];

// D7 - D8 + D a[11+ind] = dead1*(y[j]*y[j]);

// S - a[12+ind] = dead2*(y[1+j]*y[1+j]);

// T - a[13+ind] = dead3*(y[2+j]*y[2+j]);

// D1 - a[14+ind] = dead3*(y[3+j]*y[3+j]);

// D2 - a[15+ind] = dead3*(y[4+j]*y[4+j]);

// D3 - a[16+ind] = dead3*(y[5+j]*y[5+j]);

// D4 - a[17+ind] = dead3*(y[6+j]*y[6+j]);

// D5 - a[18+ind] = dead3*(y[7+j]*y[7+j]);

// D6 - a[19+ind] = dead3*(y[8+j]*y[8+j]);

// D7 - a[20+ind] = dead3*(y[9+j]*y[9+j]);

// D8 - } void changePropensityFunction(float * a, const float * y, int size ) { int j = 0;

int ind = 0;

while ( j size ) { PropensityFunction(a,y,ind,j);

j += rang;

ind += m;

} } void changeVector(float *v, int num) { v[0] = 0;

// S v[1] = 0;

// T v[2] = 0;

// D v[3] = 0;

// D v[4] = 0;

// D v[5] = 0;

// D v[6] = 0;

// D v[7] = 0;

// D v[8] = 0;

// D v[9] = 0;

// D switch( num ) { case 0:

v[1] = 1;

// S - S + T break;

case 1:

v[0] = -1;

// v[1] = 2;

// S - T + T break;

case 2:

v[0] = 1;

// S - S + S break;

case 3:

v[1] = -1;

// T - D1 + D v[2] = 2;

// break;

case 4:

v[2] = -1;

// D1 - D2 + D v[3] = 2;

// break;

case 5:

v[3] = -1;

// D2 - D3 + D v[4] = 2;

break;

case 6:

v[4] = -1;

// D3 - D4 + D v[5] = 2;

break;

case 7:

v[5] = -1;

// D4 - D5 + D v[6] = 2;

break;

case 8:

v[6] = -1;

// D5 - D6 + D v[7] = 2;

break;

case 9:

v[7] = -1;

// D6 - D7 + D v[8] = 2;

break;

case 10:

v[8] = -1;

// D7 - D8 + D v[9] = 2;

break;

case 11:

v[0] = -1;

// S - break;

case 12:

v[1] = -1;

// T - break;

case 13:

v[2] = -1;

// D1 - break;

case 14:

v[3] = -1;

// D2 - break;

case 15:

v[4] = -1;

// D3 - break;

case 16:

v[5] = -1;

// D4 - break;

case 17:

v[6] = -1;

// D5 - break;

case 18:

v[7] = -1;

// D6 - break;

case 19:

v[8] = -1;

// D7 - break;

case 20:

v[9] = -1;

// D8 - break;

} } int numberReaction(float *a, int n, float a0) { int j = -1;

float sum = 0;

float r = MyRand();

for(int i = 0;

i n;

i++) { sum += a[i];

j = i;

if(a0*r sum) { break;

} } return j;

}

 

Похожие работы:





 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.