元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

Julia tips #9: CVODE (Sundials.jl) の数値積分が終わらない時の回避方法

Julia

追記
どうやら、せっかく設定したhminは内部でフラグチェックが行われないようです。
solve関数のmaxitersオプションで対応するしかないかもしれません。 また後程検証してみます。

以下は検証が十分ではありません。
ソルバ自体の動作に支障はありませんが、期待通りに動いているかは検証する必要があります。
そのことをご理解いただいたうえでご覧ください。

2018/03/08当時の環境
- Julia 0.6.2 - Sundials 1.2.0

事の発端

ODEモデルのパラメータ推定をしていたら、非線形なモデルだとどうしても途中でソルバがフリーズしてしまう...。

原因はソルバが無限に数値積分のステップを続けてしまうことでした。

数値計算に疎いのでどうしてそんなことになってしまうのかはちゃんとわかっていないのですが、stiffなODEを解こうとすると、数値積分時に細かいステップを刻まないといけなくなるのが原因です。

「DifferentialEquationsから簡単にSundialsが使える！」と喜んだのもつかの間。

DifferentialEquations.jlでもこの件について触れられていました。

Frequently Asked Questions · DifferentialEquations.jl

公式によると、現状アルゴリズムを変えたり、ソルバを変えたりするしかないようです。

しかし、stiffなODEが解けないような状況においてCVODE以外ではお話にならないのです（経験上）。

そもそもODEが解けないようなパラメータは捨ててしまえばよいので、途中で計算を止めてなんとか捨てたいわけです。

Sundials? CVODE?

Sundialsは常微分方程式のソルバ、CVODEを作っているところです（他に幾つも作ってる）。多用しないのでわかりませんが、良く名前を聞くのでシェアは大きいと思います、MATLABにもあります。

stiffな問題をCVODEに解かせたら他のソルバより圧倒的に高速で数値積分してくれます。
stiffでなくてもめっちゃ速いです。過去記事中のリンク先でJuliaのソルバで速度比較をした方がいます。

hotoke-x.hatenablog.com

ちなみにCVODEとよく比較されるソルバにLSODAっていうのがあります。自分の場合はCVODEの方が5倍くらい速かったです。問題にもよると思うので、CVODEがダメだったら試す価値はあるかもしれません。JuliaならPkg.add("LSODA")ってやればDiffEq.jlのインターフェースから使えるようになります。 以下のブログ記事ではLSODAも良いぞ！って言ってますね。

Solving Stiff ODEs

CVODEが途中で止まる

最強のCVODEといえど限界はあります。積分のステップが小さすぎると計算が終わらなくなります。途中でエラー吐けよって思うんですが、デフォルト設定だと無限にステップを刻めることになっています。こいつをちょちょいといじれば良いのです。

以下のhminというオプションをいじる（CVODE公式ドキュメント, p37より）↓ f:id:hotoke-X:20180308031520p:plain

The default value is 0.0となっていますね。

これがいかんわけです。

Sundials.jlを書き換える（自己責任で）

DiffEq.jlのインターフェースには残念なことに、hminを設定するオプションがありません。しかし、Sundials.jl自体ではちゃんとラップしてくれているのです。じゃあ直接Sundials.jlを使えばいいんですが、exampleを見てるとわかるんですがjuliaの可読性が損なわれます。何より面倒です。

なので、DiffEq.jlのインターフェースから使えるようにSundials.jlを少しだけいじりました。

自分でやる場合は自己責任でお願いします。 ...といっても、ミスっても書き換えたところを戻せば元通りです（2敗）

DiffEq.jlのインターフェースとなるsolve.jlをSundialsのパッケージディレクトリから探します。 f:id:hotoke-X:20180308033027p:plain

hminを引数としてsolveに追加（後ろにkwargs...があるので多分これやらなくても動くけど）。

f:id:hotoke-X:20180308033609p:plain

以下のようにhminをCVODEの設定として反映させる（ハイライト部分）。

f:id:hotoke-X:20180308033746p:plain

これで完了。

SundialsはCのライブラリなんですが、Sundilas.jlの作者さんがきれいにラップしてくれているおかげでインターフェース部分を少し書き換えるだけで良いです。作者さんありがとうございます。

後は動作確認したらよいです。

パッケージを書き換えたので、初めにプリコンパイルが走ります。

（以下ではCVODEが固まるほどのstiffなODEがどんなものかわからないので、ソルバが動くかどうかだけ確かめています。）

using DifferentialEquations
using Sundials
using Plots

function parameterized_lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
alg = CVODE_BDF()
prob = ODEProblem(parameterized_lorenz, u0, tspan, p)
sol = solve(prob, hmin=1.0e-8)

plot(sol,vars=(1,2,3))

以下のようにグラフがプロットされればとりあえず正常に動作しています。

f:id:hotoke-X:20180308043930p:plain

しかし、本当に設定が反映されたのだろうか...

（どなたか、CVODEをフリーズさせる程のめちゃめちゃ硬いODEを教えてください）