绘制微分方程时的 Bokeh JS 回调

Question

我对 Python 和 Whosebug 比较陌生，如果有人问过这个问题，我深表歉意，但我已经搜索了很长时间，但无法真正找到解决我问题的方法试图做。

问题：

我一直在尝试创建一个非常基本的 COVID-19 流行病模型来熟悉 Python。我的目的是生成一个基本的 SIR 模型，用于计算人群中的易感、感染和移除个体。到目前为止，一切都很好。

我的问题是我希望绘图有一个交互式滑块来改变微分方程中的一个常数。

我正在使用 Bokeh 并尝试为此使用 Javascript 回调，但是我在使用 Javascript 时遇到了困难。到目前为止，我看到的所有示例都使用线方程，其中 y 是 x 的函数，并且代码相对简单。就我而言，由于它是一个微分方程组，我不确定我应该如何处理这个问题。

我也尝试过使用scipy，但我仍然遇到同样的问题。

代码如下。任何 help/feedback/suggestions 将不胜感激。

谢谢！

from bokeh.layouts import column, row
from bokeh.models import CustomJS, Slider
from bokeh.plotting import ColumnDataSource, figure, output_file, show

t = []

for i in range(200):
    t.append(i)

slst = []
ilst = []
rlst = []

s = 489599/489609
i = 10/489609
r = 0/489609
bet = 0.28
gam = 0.189

for f in range(200):
    ds = (-bet * (s * i))
    di = ((bet * (s * i)) - (gam * i))
    dr = gam * i
    s = s + (ds)
    slst.append(s)
    i = i + (di)
    ilst.append(i)
    r = r + (dr)
    rlst.append(r)

source = ColumnDataSource(data=dict(x=t, y=[], s=slst, i=ilst))

plot = figure(plot_width=400, plot_height=400)
plot.line('x', 'y', source=source, line_width=3, line_alpha=0.6)

callback = CustomJS(args=dict(source=source), code="""

                    ????????????

    """)

slider = Slider(start=0.1, end=4, value=1, step=.1, title="Beta ")
slider.js_on_change('value', callback)

layout = column(slider, plot)

show(layout)

Answer 1

这是我的想法。有趣的是，随着 beta 值较高，易感线低于 0。也许我在将您的代码移植到 JavaScript 时犯了一个错误 - 如果是这样请纠正我。

from bokeh.core.property.instance import Instance
from bokeh.io import save
from bokeh.layouts import column
from bokeh.model import Model
from bokeh.models import CustomJS, Slider, Callback
from bokeh.plotting import ColumnDataSource, figure

source = ColumnDataSource(data=dict(t=[], s=[], i=[], r=[]))

plot = figure(plot_width=400, plot_height=400)
plot.line('t', 's', source=source, line_width=3, line_alpha=0.6)
plot.line('t', 'i', source=source, line_width=3, line_alpha=0.6, color='orange')
plot.line('t', 'r', source=source, line_width=3, line_alpha=0.6, color='green')

callback = CustomJS(args=dict(source=source), code="""\
    const N = 200;
    let s = 489599 / 489609;
    let i = 10 / 489609;
    let r = 0 / 489609;
    const bet = cb_obj.value;
    const gam = 0.189;
    const tlst = source.data.t = [];
    const slst = source.data.s = [];
    const ilst = source.data.i = [];
    const rlst = source.data.r = [];
    for (let t = 0; t < N; ++t) {
        s -= bet * s * i;
        i += bet * s * i - gam * i;
        r += gam * i;
        tlst.push(t);
        slst.push(s);
        ilst.push(i);
        rlst.push(r);
    }
    source.change.emit();
""")

slider = Slider(start=0.1, end=4, value=1, step=.1, title="Beta ")
slider.js_on_change('value', callback)


class IdleDocObserver(Model):
    """Work around https://github.com/bokeh/bokeh/issues/4272."""
    on_idle = Instance(Callback)

    # language=TypeScript
    __implementation__ = """\
        import {View} from "core/view"
        import {Model} from "model"
        import * as p from "core/properties"

        export class IdleDocObserverView extends View {}

        export namespace IdleDocObserver {
            export type Attrs = p.AttrsOf<Props>
            export type Props = Model.Props & {on_idle: p.Property<any>}
        }

        export interface IdleDocObserver extends IdleDocObserver.Attrs {}

        export class IdleDocObserver extends Model {
            static init_IdleDocObserver(): void {
                this.prototype.default_view = IdleDocObserverView
                this.define<IdleDocObserver.Props>({on_idle: [p.Any]})
            }

            _doc_attached(): void {
                super._doc_attached()
                const execute = () => this.on_idle!.execute(this)
                if (this.document!.is_idle)
                    execute();
                else
                    this.document!.idle.connect(execute);
            }
        }
    """


idle_doc_observer = IdleDocObserver(on_idle=CustomJS(args=dict(callback=callback, slider=slider),
                                                     code="callback.execute(slider);"))

layout = column(slider, plot)
save([idle_doc_observer, layout])

Answer 2

根据我的评论，可以在 javascript 实现中使用 RK4 将 ODE 集成到 html 文档中。在 bokeh 中似乎没有办法在回调、实用函数和常见计算之外实现 javascript 函数。因此，为避免代码重复，必须使一个回调足够通用，以便它可以为所有滑块更改事件提供服务。（或者，可以实现一个 "recompute" 按钮。）

为了看起来更专业，制作 2 个图，一个用于所有组件，一个用于 I 单独。

# Set up the plots and their data source
source = ColumnDataSource(data=dict(T=[], S=[], I=[], R=[]))

SIR_plot = figure(plot_width=400, plot_height=400)
SIR_plot.line('T', 'S', source=source, legend_label="S", line_width=3, line_alpha=0.6, color='blue')
SIR_plot.line('T', 'I', source=source, legend_label="I", line_width=3, line_alpha=0.6, color='orange')
SIR_plot.line('T', 'R', source=source, legend_label="R", line_width=3, line_alpha=0.6, color='green')

I_plot = figure(plot_width=400, plot_height=400)
I_plot.line('T', 'I', source=source, line_width=3, line_alpha=0.6, color='orange')

接下来为可能想要影响的参数设置 4 个滑块

# declare the interactive interface elements
trans_rate = Slider(start=0.01, end=0.4, value=0.3, step=.01, title="transmission rate ")
recov_rate = Slider(start=0.01, end=0.4, value=0.1, step=.01, title="recovery rate")

I_init = Slider(start=0.01, end=0.1, value=0.05, step=.002, title="initial infected [proportion] ")
max_time = Slider(start=10, end=200, value=50, step=1, title="time range [days] ")

现在作为对 Eugene Pakhomov 答案的主要更改，对所有滑块进行一次回调（参见散景画廊滑块演示）并使用矢量化 RK4 方法进行 ODE 积分

callback = CustomJS(args=dict(source=source, I_init=I_init, max_time=max_time, 
                              trans_rate=trans_rate, recov_rate=recov_rate), 
                    code="""\
    let i = I_init.value;
    let s = 1-i;
    let r = 0;
    const bet = trans_rate.value;
    const gam = recov_rate.value;
    let tf = max_time.value;
    const dt = 0.1;
    const tlst = source.data.T = [0];
    const slst = source.data.S = [s];
    const ilst = source.data.I = [i];
    const rlst = source.data.R = [r];

    function odefunc(t,sir) {
        let tr = bet*sir[0]*sir[1];
        let rc = gam*sir[1];
        return [-tr, tr-rc, rc];
    }
    let sir = [s,i,r];
    for (let t = 0; t < tf; t+=dt) {
        sir = RK4Step(t,sir,dt);
        tlst.push(t+dt);
        slst.push(sir[0]);
        ilst.push(sir[1]);
        rlst.push(sir[2]);
    }
    source.change.emit();

    function axpy(a,x,y) { 
        // returns a*x+y for arrays x,y of the same length
        var k = y.length >>> 0;
        var res = new Array(k);
        while(k-->0) { res[k] = y[k] + a*x[k]; }
        return res;
    }

    function RK4Step(t,y,h) {
        var k0 = odefunc(t      ,               y );
        var k1 = odefunc(t+0.5*h, axpy(0.5*h,k0,y));
        var k2 = odefunc(t+0.5*h, axpy(0.5*h,k1,y));
        var k3 = odefunc(t+    h, axpy(    h,k2,y));
        // ynext = y+h/6*(k0+2*k1+2*k2+k3);
        return axpy(h/6,axpy(1,k0,axpy(2,k1,axpy(2,k2,k3))),y);
    }

""")
trans_rate.js_on_change('value', callback)
recov_rate.js_on_change('value', callback)

I_init.js_on_change('value', callback)
max_time.js_on_change('value', callback)

最后，以某种布局将所有内容串在一起

# generate the layout

parameters_panel = column(trans_rate, recov_rate)
initials_panel = column(I_init,max_time)

plots = row(SIR_plot, I_plot)
inputs = row(parameters_panel, initials_panel)

simulation = column(plots, inputs)

show(simulation)

为了避免最初的空白图，我参考了 Eugene Pakhomov 的回答，因为它是在第一个滑块移动后出现的图。

绘制微分方程时的 Bokeh JS 回调

Bokeh JS callback when plotting differential equations

javascript

python

differential-equations

bokeh

bokehjs