如何监控SciPy.odeint的进程？

Question

SciPy 可以通过 scipy.integrate.odeint 或其他包求解 ode 方程，但它在函数完全求解后给出结果。但是，如果 ode 函数非常复杂，程序将花费大量时间（一两天）才能给出整个结果。那么我如何监控求解方程的步骤（当方程尚未完全求解时打印出结果）？

Answer 1

当我在谷歌上寻找答案时，我找不到满意的答案。因此，我使用tqdm项目通过概念验证解决方案提出了一个简单的要点。希望对您有帮助。

\n

编辑：版主要求我对上面链接中发生的情况进行解释。

\n

首先，我使用 scipy 的 OOP 版本odeint( solve_ivp)，但您可以将其改编回odeint. T0假设您希望时不时地进行集成，T1并且希望显示每 0.1% 的进度。您可以修改 ode 函数以采用两个额外参数，a pbar（进度条）和 a state（当前积分状态）。就像这样：

\n

def fun(t, y, omega, pbar, state):\n    # state is a list containing last updated time t:\n    # state = [last_t, dt]\n    # I used a list because its values can be carried between function\n    # calls throughout the ODE integration\n    last_t, dt = state\n    \n    # let\'s subdivide t_span into 1000 parts\n    # call update(n) here where n = (t - last_t) / dt\n    time.sleep(0.1)\n    n = int((t - last_t)/dt)\n    pbar.update(n)\n    \n    # we need this to take into account that n is a rounded number.\n    state[0] = last_t + dt * n\n    \n    # YOUR CODE HERE\n    dydt = 1j * y * omega\n    return dydt\n

\n

这是必要的，因为函数本身必须知道它所在的位置，但 scipy 并odeint没有真正将此上下文提供给函数。然后，您可以fun与以下代码集成：

\n

T0 = 0\nT1 = 1\nt_span = (T0, T1)\nomega = 20\ny0 = np.array([1], dtype=np.complex)\nt_eval = np.arange(*t_span, 0.25/omega)\n\nwith tqdm(total=1000, unit="\xe2\x80\xb0") as pbar:\n    sol = solve_ivp(\n        fun,\n        t_span,\n        y0,\n        t_eval=t_eval,\n        args=[omega, pbar, [T0, (T1-T0)/1000]],\n    )\n

\n

请注意， 中任何可变的（如列表）args都会被实例化一次，并且可以在函数内进行更改。我建议这样做而不是使用全局变量。

\n

这将显示一个进度条，如下所示：

\n

100%|\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x88\xe2\x96\x89| 999/1000 [00:13<00:00, 71.69\xe2\x80\xb0/s]\n

\n

Answer 2

您可以分割积分域并积分各个段，将前一个段的最后一个值作为下一个段的初始条件。在此期间，打印出您想要的任何内容。如有必要，请使用 numpy.concatenate 组装各个部分。

在三体太阳系模拟的标准示例中，替换代码

u0 = solsys.getState0();
t = np.arange(0, 100*365.242*day, 0.5*day);
%timeit u_res = odeint(lambda u,t: solsys.getDerivs(u), u0, t, atol = 1e11*1e-8, rtol = 1e-9)

输出：1 loop, best of 3: 5.53 s per loop

带有进度报告代码

def progressive(t,N):
    nk = [ int(n+0.5) for n in np.linspace(0,len(t),N+1) ]
    u0 = solsys.getState0();
    u_seg = [ np.array([u0]) ];
    for k in range(N):
        u_seg.append( odeint(lambda u,t: solsys.getDerivs(u), u0, t[nk[k]:nk[k+1]], atol = 1e11*1e-8, rtol = 1e-9)[1:] )
        print t[nk[k]]/day
        for b in solsys.bodies: print("%10s %s"%(b.name,b.x))
    return np.concatenate(u_seg)
%timeit u_res = progressive(t,20)

输出：1 loop, best of 3: 5.96 s per loop

显示控制台打印的开销仅为 8%。有了更实质性的 ODE 功能，报告开销的比例将显着减少。

也就是说，Python，至少其标准包，并不是工业规模的数字运算工具。始终使用具有强类型变量的编译版本，以尽可能减少解释开销。

• 使用一些经过大量开发和测试的软件包，例如 Sundials 或 julia-lang 框架 Differentialequations.jl，直接以适当的编译语言编码 ODE 函数。使用高阶方法获得更大的步长，从而获得更小的步长。测试使用隐式或指数/Rosenbrock 方法是否会进一步减少每个固定间隔的步骤数或 ODE 函数评估。加速比的差异可能是 10 到 100 倍。

• 使用上述 Python 包装器以及 ODE 函数的一些加速友好实现。

• 使用准源代码翻译工具 JITcode 将您的 Python ODE 函数翻译为 C 指令的意大利面条列表，然后给出一个编译函数，可以（几乎）直接从编译后的 odeint 的 FORTRAN 内核调用该函数。