Монада моделирования. Часть вторая. Первые шаги в 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#.