воскресенье, 28 февраля 2010 г.

Потоки (Streams)

Вчера прочитал главу SICP, посвященную потокам (streams). Возникло острое желание переписать примеры на Common Lisp (CL). Тема очень интересна сама по себе. Еще подстегивало то, что некоторые активисты упорно пропагандируют идею оторванности CL от функциональной парадигмы (ФП), как бы это абсурдно ни звучало. Но когда я приступил к переписыванию кода на CL, мое первое обманчивое впечатление было таким, что потоки будет реализовать труднее, чем в Схеме. Например, в CL нет аналога схемовского DEFINE, который бы позволил определять переменные рекурсивно. Но, к счастью, все разрешилось удачным образом. Как и должно было быть, выручили макросы, могучая сила CL. Они же мне позволили в свое время добавить очень простой и удобный синтаксический сахар для монад по типу нотации do, о чем я писал ранее.

Итак, все начинается с ленивости. Для этого определяются конструкции DELAY и FORCE, но прежде нужна вспомогательная функция MEMO, которая возвращает функцию без аргументов, результат которой кешируется:

(defun memo (fun)
  (let ((already-run? nil)
        (result nil))
    #'(lambda ()
        (if (not already-run?)
            (progn
              (setf result (funcall fun))
              (setf already-run? t)
              result)
            result))))

Тут все по книге SICP, только там функция называлась MEMO-PROC.

Дальше определяем конструкции DELAY и FORCE, причем DELAY должен быть непременно макросом, чтобы он мог упрятать свой аргумент в лямбду, сделав его тем самым ленивым:

(defmacro delay (exp)
  `(memo #'(lambda() ,exp)))

(defun force (delayed-exp)
  (funcall delayed-exp))

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

(defmacro cons-stream (a b)
  `(cons ,a (delay ,b)))

(defun stream-car (stream)
  (car stream))

(defun stream-cdr (stream)
  (force (cdr stream)))

(defun stream-null (stream)
  (null stream))

(defparameter *empty-stream* nil)

Функция STREAM-CDR раскрывает ленивую CDR-часть потока, если она не была до того уже раскрыта. Как помним, на нижнем уровне за это отвечает функция MEMO, которая встраивается при создании каждой пары. В общем, это все уже описано в SICP. Поэтому останавливаться не буду.

Следующая функция STREAM-REF аналогична AREF, но работает уже с потоками.

(defun stream-ref (stream n)
  (loop for i from 0 to n
     for s = stream then (stream-cdr s)
     finally (return (stream-car s))))

Несмотря на присутствие в реализации совсем нефункционального LOOP, функция является чистой.

Далее идут отображения для потоков. Вполне в духе ФП.

(defun stream-map (fun stream)
  (if (stream-null stream) stream
      (cons-stream (funcall fun (stream-car stream))
                   (stream-map fun (stream-cdr stream)))))

(defun stream-map2 (fun s1 s2)
  (cond
    ((stream-null s1) s1)
    ((stream-null s2) s2)
    (t (cons-stream (funcall fun (stream-car s1) (stream-car s2))
                    (stream-map2 fun (stream-cdr s1) (stream-cdr s2))))))

Далее понадобится функция фильтрации потока. Я поменял имя STREAM-FILTER на более идиоматическое. К тому же задействовал LOOP. При этом функция остается чистой.

(defun stream-remove-if-not (test stream)
  (loop for s = stream then (stream-cdr s)
     when (stream-null s) do (return s)
     when (funcall test (stream-car s)) do
       (return
         (cons-stream (stream-car s)
                      (stream-remove-if-not test (stream-cdr s))))))

Чтобы опробовать новые возможности в деле, понадобятся еще итератор и функция вывода:

(defun stream-iter (fun stream)
  (loop for s = stream then (stream-cdr s)
       while (not (stream-null s))
       do (funcall fun (stream-car s))))

(defun print-stream (stream)
  (stream-iter #'(lambda (x) (format t "~%~a" x)) stream))

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

CL-USER> (print-stream (cons-stream 1 (cons-stream 2 (cons-stream 3 nil))))

1
2
3
NIL

Подходим к главному препятствию. В схеме мы могли легко определить поток рекурсивно. Например, поток единиц задавался бы так:

(define ones (cons-stream 1 ones))  ;; Scheme!

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

(defparameter *ones*
  (recurrent-let ((ones (cons-stream 1 ones)))
    ones))

Здесь RECURRENT-LET похож на LET с одной переменной. Отличие заключается в том, что к самой переменной можно обращаться внутри определения. То есть, определение может быть рекурсивным.

Сам макрос RECURRENT-LET достаточно прост:

(defmacro recurrent-let (((name value)) &body body)
  (let ((x (gensym)))
    `(let ((,x (cons nil nil)))
       (symbol-macrolet ((,name (force (car ,x))))
         (setf (car ,x) (delay ,value))
         ,@body))))

Вот, во что будет раскрыта внутренняя часть определения параметра *ONES*:

(LET ((#:G764 (CONS NIL NIL)))
  (SYMBOL-MACROLET ((ONES (FORCE (CAR #:G764))))
    (SETF (CAR #:G764) (DELAY (CONS-STREAM 1 ONES)))
    ONES))

Мы создаем некий объект с одним полем. Это поле содержит ленивое значение переменной. Любое (рекурсивное) обращение к переменной мы трактуем как попытку немедленно вычислить ее значение, если оно еще не вычислено. Трюк работает, потому что мы успеваем присвоить полю объекта ленивое значение прежде первого обращения к этой переменной. Это гарантирует использованный макрос DELAY.

Стоит заметить, что данное представление бесконечного потока единиц на самом деле занимает очень мало памяти, поскольку поток ссылается на самого себя.

Если присмотреться к определению макроса, то можно увидеть, что внутри определения рекурсивной переменной мы можем использовать еще один RECURRENT-LET, а также любую другую конструкцию, включая LET, FLET и LABELS. Это открывает путь к вложенным и более сложным взаимно-рекурсивным определениям, что и будет продемонстрировано далее.

В соответствии с SICP определим вспомогательные функции:

(defun add-streams (stream1 stream2)
  (stream-map2 #'+ stream1 stream2))

(defun scale-stream (stream factor)
  (stream-map #'(lambda (x) (* x factor)) stream))

Далее определим поток целых чисел:

(defun integers-starting-from (n)
  (cons-stream n (integers-starting-from (+ n 1))))

(defparameter *integers-alpha*
  (integers-starting-from 1))

Этот же самый поток мы можем определить иначе:

(defparameter *integers*
  (recurrent-let
      ((integers (cons-stream 1 (add-streams *ones* integers))))
     integers))

Естественно, без определения потока чисел Фибоначчи мое сообщение можно было бы считать неполным:

(defparameter *fibs*
  (recurrent-let
      ((fibs (cons-stream 0
                          (cons-stream 1
                                       (add-streams (stream-cdr fibs)
                                                    fibs)))))
    fibs))

Теперь можем узнать, каким будет 1001-ое число Фибоначчи:

CL-USER> (stream-ref *fibs* 1001)
70330367711422815821835254877183549770181269836358732742604905087154537118196933579742249494562611733487750449241765991088186363265450223647106012053374121273867339111198139373125598767690091902245245323403501

Как и следовало ожидать, ответ был немедленным.

Следуя SICP, далее приведу пример определения потока простых чисел. Но прежде мне понадобятся две простые утилиты.

(defun square (n) (* n n))
(defun divisible-p (x y) (= (mod x y) 0))

Сам поток определен ниже. Это также пример взаимно-рекурсивного определения.

(defparameter *primes*
  (recurrent-let
      ((primes
        (flet ((prime-p (n)
                 (labels ((iter (ps)
                            (cond ((> (square (stream-car ps)) n) t)
                                  ((divisible-p n (stream-car ps)) nil)
                                  (t (iter (stream-cdr ps))))))
                   (iter primes))))
          (cons-stream
           2
           (stream-remove-if-not #'prime-p (integers-starting-from 3))))))
    primes))

Еще один пример — решатель дифференциальных уравнений. Сначала определяем интеграл, где интегрируемая функция передается лениво:

(defun integral (delayed-integrand initial-value dt)
  (recurrent-let
      ((int (cons-stream initial-value
                         (let ((integrand (force delayed-integrand)))
                           (add-streams (scale-stream integrand dt)
                                        int)))))
    int))

Теперь сам решатель:

(defun solve (f y0 dt)
  (recurrent-let
      ((y (recurrent-let
              ((dy (stream-map f y)))
            (integral (delay dy) y0 dt))))
    y))

Можем проверить как и в SICP:

CL-USER> (stream-ref (solve #'(lambda (y) y) 1 0.001) 1000)
2.7169204

Как видим, примеры достаточно легко ложатся на CL, хотя и выглядят несколько иначе. При этом они вполне соответствуют духу ФП. Разве что, использование LOOP несколько необычно для этой области, но это прекрасный пример сочетания разных подходов к программированию. Основа остается несомненно функциональной. Неужели после этого может кто-то по-прежнему утверждать, что CL не является языком функционального программирования?

суббота, 6 февраля 2010 г.

Монада моделирования. Код

Открыл новый проект Aivika на SourceForge. Там лежат исходники моделирующей библиотеки для F#. Пожалуй, функцию memo нужно будет оптимизировать. Вероятно, из-за нее тормозит. Но есть задача важнее. Теперь буду пытаться продвигать идею на инглише.

пятница, 5 февраля 2010 г.

Монада моделирования. Часть вторая. Первые шаги в Common Lisp.

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

В общем, это все уже работает на F#, хотя и медленно. Рабочее название соответствующей библиотеки для F# – Aivika. Теперь я перенес базовую часть на Common Lisp, назвав новую библиотеку Salika. Соответствующий монадический макрос назвал как WITH-DYNAMICS-MONAD. Этот макрос полностью соответствует интерфейсу монадических макросов из пакета cl-monad-macros, о котором я писал в своем блоге ранее. Вкратце, такие макросы добавляют удобный синтаксический сахар по типу нотации do из хаскеля.

Возьмем простую динамическую систему, описанную на языке MapSim:

A = integ (-F, 100);
B = integ (F - G, 0);
C = integ (G, 0);

F = ka * A;
G = kb * B;

ka = 1;
kb = 1;

Функция integ задает интеграл. В первом параметре функции определяется производная по времени. Во втором – начальное значение интеграла. Интегрирование происходит на промежутке starttime <= t <= stoptime с шагом dt. Здесь A(t), B(t) и C(t) – интегралы. Их еще иногда называют в системной динамике резервуарами. Функции F(t) и G(t) называют дополнительными. В системе также заданы константы ka и kb.

В первой ссылке есть пример того, как эту систему можно задать на языке F#. Там получается достаточно близко к математической нотации. Во многом благодаря тому, что для типа монады моделирования легко определить арифметические операции и стандартные функции типа синуса и косинуса. Такой подход можно встретить, например, в книге “The Haskell School of Expression” автора Paul Hudak. В случае же Common Lisp я предпочитаю использовать монадические макросы напрямую, поскольку это порождает в целом более эффективный код. Думаю, что наглядность при этом страдает несильно.

Итак, такую систему мы можем описать на Common Lisp следующим образом.

(defun create-process ()
  (let* ((ka 1d0)
         (kb 1d0))
    (let* ((integ-a (make-instance 'integ-stock :initial-value 100d0))
           (integ-b (make-instance 'integ-stock :initial-value 0d0))
           (integ-c (make-instance 'integ-stock :initial-value 0d0)))
      (let* ((f (with-dynamics-monad
                  (let! ((a (current-value integ-a)))
                    (unit (* ka a)))))
             (g (with-dynamics-monad
                  (let! ((b (current-value integ-b)))
                    (unit (* kb b))))))
        (setf (outflow integ-a) f)
        (setf (inflow integ-b) f)
        (setf (outflow integ-b) g)
        (setf (inflow integ-c) g)
        (current-value integ-a)))))

Единственный момент – в примере возвращается лишь интеграл A, но как часть всей системы. Здесь монаду можно увидеть в том, как задаются динамические процессы для переменных F и G. Внутри макроса WITH-DYNAMICS-MONAD макроc LET! эквивалентен стрелке из нотации do, а макрос UNIT – функции return для монады. Получается, что F и G возвращают монады. Более того, выражение (current-value integ-a) – тоже монада. Это представления динамических процессов.

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

(defun run-process ()
  (with-dynamics-monad
    (run! (create-process) (make-specs :start-time 0.0 :stop-time 10.0 :method 'runge-kutta-4 :dt 0.1)))))

Среда лиспа проинтегрирует модель методом Рунге-Кутта четвертого порядка, начиная от точки 0 и заканчивая точкой 10 с основным шагом интегрирования 0,1.

SALIKA> (run-process)
(100.0d0 90.48374986516933d0 81.8730898966253d0 74.08184186894766d0
         67.03202849220888d0 60.65309299043932d0 54.88119294695767d0
         49.658561349146126d0 44.932928437803035d0 40.65699857475723d0
         36.787976893068794d0 33.28714099238066d0 30.119453392811955d0
         27.253210868708226d0 24.65972715266909d0 22.313045834254343d0
         ...)

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

Тем не менее, во время интегрирования таких уравнений внутри используется сугубо императивная функция мемоизации MEMO-PROCESS. Она аналогична функции memo из F#-вской библиотеки. Обе эти функции принимают динамический процесс, т.е. некоторое значение в монаде моделирования, и возвращают процесс-двойник, который кеширует все значения первого процесса в узлах интегрирования, причем вычисления происходят строго последовательно по узлам. Это открывает путь к построению гибридных моделей, поскольку мы здесь знаем, что первый, т.е. кешируемый процесс будет вычисляться в строго определенном порядке (если нет других ссылок на него). Такой процесс мы можем реализовать как итерационный. Это еще одна вещь наряду с недетерминизмом, которую будет трудно или почти невозможно воспроизвести в хаскеле.

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

Это все прекрасно, но у моего метода есть большой недостаток. Монада моделирования работает медленно, очень медленно. Например, MapSim компилирует симуляцию в эффективный байт-код, который выполняется гораздо быстрее. Или, быть может, просто MapSim быстр? Но монада моделирования позволяет достаточно просто и легко строить очень сложные гибридные модели. И пока я не определился в своем отношении к изобретению. Просто игрушка для ума или серьезная и полезная вещь?

Версия для F# достаточно зрелая и многое умеет, включая стохастику. Версия для Common Lisp реализует пока лишь интегралы и мемоизацию, причем нет оптимизации по типам. Может быть, потом выложу результаты в свободный доступ. Хотя мне кажется, что описанное в обеих частях достаточно легко воспроизвести, поскольку первая часть содержит полное описание монады моделирования на языке F#.

среда, 3 февраля 2010 г.

Монада моделирования

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

Но с этой реализацией идеи возникла проблема, которую я пока не знаю как красиво разрешить в хаскеле. Моделируемая динамическая система на самом деле может быть стохастической, т.е. недетерминированной. Каждый новый прогон модели может порождать новые результаты. То есть, запускающая имитацию функция не является чистой. Тут на помощь пришел язык F#, который я начал изучать чуть позже. В нем уже достаточно легко сочетать функциональный и императивный подходы. В итоге это привело к любопытному результату, оформленному в виде библиотеки на F# с рабочим названием Aivika, которую я здесь и попытаюсь описать.

Сначала приведу пример простой динамической системы, записанной на языке MapSim:

A = integ (-F, 100);
B = integ (F - G, 0);
C = integ (G, 0);

F = ka * A;
G = kb * B;

ka = 1;
kb = 1;

Здесь функция integ задает интеграл. Ее первым параметром является производная по времени, а вторым – начальное значение. То есть, A(t) является интегралом (еще говорят резервуаром), производная которого равна dA/dt = –ka*A(t). В начальный момент времени значение A(t0) равно 100. Это задача Коши.

Такая система может быть интегрирована методами Эйлера и Рунге-Кутта второго и четвертого порядков. Задается начальное время starttime, конечное stoptime и основной шаг интегрирования dt. Затем временная шкала разбивается на узлы: starttime, starttime + dt, starttime + 2*dt, …, stoptime. В случае методов Рунге-Кутта создаются еще дополнительные внутренние под-узлы. В итоге каждый узел сетки интегрирования однозначно задается парой целых чисел: номером итерации и номером фазы. Метод Эйлера имеет всего одну фазу. Метод Рунге-Кутта второго порядка – две фазы. Метод четвертого порядка – четыре фазы.

type Iteration = int
type Phase = int
type Method = Euler | RungeKutta2 | RungeKutta4

Здесь и далее я буду использовать F# для записи кода.

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

type Randomness = SimpleRnd | StrongRnd

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

type Specs = {
    spcStartTime: float; spcStopTime: float; spcDT: float;
    spcMethod: Method; spcRandomness: Randomness
}

Для нашего примера мы могли бы взять следующие параметры: начальное время – 0, конечное время – 10, основной шаг интегрирования – 0,1, метод – Рунге-Кутта четвертого порядка, а режим генератора – «сильная» псевдо-случайность, хотя последний никак не используется в этой модели.

let specs = {
    spcStartTime=0.0; spcStopTime=10.0; spcDT=0.1;
    spcMethod=RungeKutta4; spcRandomness=StrongRnd
}

Саму систему мы можем начать описывать с констант ka и kb.

let ka = 1.0
let kb = 1.0

Затем объявляем интегралы и задаем для них начальные значения.

let a = Integ 100.0
let b = Integ 0.0
let c = Integ 0.0

Каждый из трех интегралов имеет свойства Value, Inflow и Outflow. Эти свойства связаны дифференциальными уравнениями, одно из которых я покажу на примере интеграла a:

d/dt (a.Value) = (a.Inflow – a.Outflow)

Такие же уравнения имеют место для других интегралов. Свойство Value возвращает текущее значение интеграла, и оно доступно только по чтению. Свойства Inflow и Outflow можно менять. Они оба вместе задают производную. По-умолчанию, они представляют нулевые значения. То есть, если мы не зададим ни одно из свойств Inflow и Outflow, то производная будет равна нулю, а интеграл окажется константой.

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

let f = ka * a.Value
let g = kb * b.Value

И завершаем описание задачи собственно самими дифференциальными уравнениями:

a.Inflow <- -f
b.Inflow <- f – g
c.Inflow <- g

Удивительно, но это все – что нам нужно, чтобы немедленно запросить библиотеку промоделировать, скажем, переменную a и вернуть результат:

> runDynamics a.Value specs;;
val it : float [] =
  [|100.0; 90.48375; 81.87309014; 74.0818422; 67.03202889; 60.65309344;
    54.88119344; 49.65856187; 44.93292897; 40.65699912; 36.78797744;
    33.28714154; 30.11945393; 27.2532114; 24.65972767; 22.31304633;
    20.18968106; 18.26838054; 16.52991577; 14.95688766; 13.53355284;
    12.24566612; 11.08033792; 10.02590526; 9.071815051; 8.208518451;
    7.427375314; 6.720567711; 6.081021686; 5.50233646; 4.978720367;...|]

Более того, мы можем запросить результаты не только для одной переменной, но и для целой группы, попросив включить в список и модельное время:

> runDynamics Dynamics.ofArray [| time; a.Value; b.Value; c.Value |] specs;;
[|[|0.0; 100.0; 0.0; 0.0|];
  [|0.1; 90.48375; 9.048333333; 0.4679166667|];
  [|0.2; 81.87309014; 16.37454263; 1.752367234|];
  [|0.3; 74.0818422; 22.22445032; 3.693707481|];
  [|0.4; 67.03202889; 26.81268809; 6.155283021|];
  [|0.5; 60.65309344; 30.32640707; 9.020499487|]; ...|]

Поразительная краткость, не правда ли?

Все это стало возможным благодаря использованию монады моделирования, которая является вариацией стандартной монады Reader. Еще меня натолкнула на ее написание замечательная книга «The Haskell School of Expression» автора Paul Hudak.

Итак, мы моделируем некоторый динамический процесс, который в каждый момент времени возвращает некоторое значение типа 'a. Обозначим этот полиморфный тип как Dynamics<'a>.

[]
type Dynamics<'a> (iter: Iterator<'a>) =
    member d.Iterator = iter

let iterator (p: Dynamics<'a>) = p.Iterator

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

type Iterator<'a> = Run -> Specs -> Iteration -> Phase -> 'a

То есть, итератор – это просто функция. Здесь появился новый тип Run, c которым связан недетерминизм и не только. В упрощенном виде (после удаления некоторых оптимизаций) этот тип выглядит так.

type Run () =
    let mutable disposed = false
    let disposedEvent = new Event<_>()
    member x.Disposed = disposedEvent.Publish
    member x.IsDisposed = disposed
    interface IDisposable with
        member x.Dispose() =
            if not disposed then
                disposedEvent.Trigger()
                disposed <- true

Когда мы запускаем имитацию, мы создаем уникальный экземпляр типа Run. Затем во время интегрирования мы передаем это значение в итератор динамического процесса. Ориентируясь на такое значение, любая функция может создать некоторые внутренние данные, которые будут актуальны, пока существует значение типа Run. Например, это может быть таблица кеширования предыдущих значений для интеграла внутри объекта типа Integ. Когда имитация заканчивается, значение типа Run разрушается, но предварительно вызывается метод Dispose, который посылает сигнал Disposed. Одна из внутренних функций интеграла или уровнем ниже ловит этот сигнал, и разрушает таблицу кеширования. Так достигается следующее.

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

Сам запуск динамического процесса очень прост. Его вызов уже использовался в примерах выше.

let runDynamics d specs = [|

    use r = new Run ()
   
    let f1 = iterator d r specs
    let f2 = fun n -> f1 n 0
   
    for i = 0 to (Specs.iterations specs - 1) do yield f2 i
|]

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

Собственно, все. Дальше монада моделирования определяется просто. Ее функции return и bind почти такие же, как и в случае монады Reader.

[]
type DynamicsBuilder() =
    member d.Return (a) =
        new Dynamics<'a> (fun r s n ph -> a)
    member d.Bind (m, k) =
        new Dynamics<'b> (fun r s n ph ->
            let a = iterator m r s n ph
            let m' = k a
            iterator m' r s n ph)

Тогда мы можем определить свой workflow. Назовем его dynamics.  Он позволяет использовать синтаксический сахар в F#, похожий на нотацию do из хаскеля, где let! будет эквивалентен стрелке.

let dynamics = new DynamicsBuilder()

Более того, мы можем определить арифметические операции и основные функции типа синуса в соответствии с образцом:

type Dynamics<'a> with
    static member (+) (a: Dynamics, b: Dynamics) =
        Dynamics.lift2 (+) a b

Здесь функция lift2 достаточно идиоматична. Она берет некоторую функцию и два значения, обернутые в монаду. Значения извлекаются, к ним применяется функция, и результат снова заворачивается в монаду.

module Dynamics =
    let lift2 f m1 m2 = dynamics {
        let! a1 = m1
        let! a2 = m2
        return f a1 a2}

Но я использовал версию оптимальнее, которая уже учитывает внутреннее устройство монады:

module Dynamics =
    let lift2 f m1 m2 = Dynamics<_> (fun r s n ph ->
        let a1 = iterator m1 r s n ph
        let a2 = iterator m2 r s n ph
        in f a1 a2)

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

module Dynamics =
    val memo: Dynamics<'a> -> Dynamics<'a>

Здесь для краткости приведу лишь функцию integ – аналог одноименной функции MapSim, причем только для интегрирования по методу Рунге-Кутта второго порядка. Для других методов определены аналогичные внутренние функции integEuler и integRK4.

    let integ (f: Dynamics) (i: Dynamics) =
       
        // integEuler and integRK4
       
        let integRK2 (y: Dynamics<_>) (f: Dynamics<_>) (i: Dynamics<_>) r s =
            let vy n = iterator y r s n 0
            let vi n = iterator i r s n 0
            let k1 n = iterator f r s n 0
            let k2 n = iterator f r s n 1
            fun n ph ->
                match n, ph with
                    | 0, 0 -> vi 0
                    | n, 0 -> vy (n-1) + s.spcDT/2.0 * (k1 (n-1) + k2 (n-1))
                    | n, 1 -> vy n + s.spcDT * k1 n
                    | _ -> failwithf "integRK2: incorrect phase = %i" ph
       
        let rec y = Dynamics.memo z
        and z = Dynamics<_> (fun r s ->
            match s.spcMethod with
                | Euler -> integEuler y f i r s
                | RungeKutta2 -> integRK2 y f i r s
                | RungeKutta4 -> integRK4 y f i r s)
        in y

Но одной функции integ недостаточно в F#, потому что нужно разруливать как-то рекурсивные связи. Поэтому придумал трюк с отдельным типом Integ и тремя его свойствами Value, Inflow и Outflow. Сам тип Integ использует внутри функцию integ.

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

В заключение несколько слов о скорости имитации. Она ужасна. Мои тесты показали, что библиотека Aivika медленнее MapSim в десятки или сотни раз, хотя и там, и там – .NET. И автор один и тот же J.

Теперь вот думаю, а не написать ли мне подобное на Common Lisp – там некоторые вещи можно сделать быстрее за счет более умной кодогенерации.