如何更好地使用 Cython 更快地求解微分方程？
How can I use Cython well to solve a differential equation faster?
我想缩短 Scipy 的 odeint 求解微分的时间
方程。
为了练习，我使用了 Python in scientific computations 中的示例作为模板。因为 odeint 将函数 f 作为参数，所以我将此函数编写为静态类型的 Cython 版本并希望
odeint 的 运行ning 时间将显着减少。
函数 f 包含在名为 ode.pyx 的文件中，如下所示：
import numpy as np
cimport numpy as np
from libc.math cimport sin, cos
def f(y, t, params):
cdef double theta = y[0], omega = y[1]
cdef double Q = params[0], d = params[1], Omega = params[2]
cdef double derivs[2]
derivs[0] = omega
derivs[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)
return derivs
def fCMath(y, double t, params):
cdef double theta = y[0], omega = y[1]
cdef double Q = params[0], d = params[1], Omega = params[2]
cdef double derivs[2]
derivs[0] = omega
derivs[1] = -omega/Q + sin(theta) + d*cos(Omega*t)
return derivs
然后我创建一个文件setup.py来编译函数：
from distutils.core import setup
from Cython.Build import cythonize
setup(ext_modules=cythonize('ode.pyx'))
求解微分方程的脚本（也包含Python
f) 的版本称为 solveODE.py，看起来像：
import ode
import numpy as np
from scipy.integrate import odeint
import time
def f(y, t, params):
theta, omega = y
Q, d, Omega = params
derivs = [omega,
-omega/Q + np.sin(theta) + d*np.cos(Omega*t)]
return derivs
params = np.array([2.0, 1.5, 0.65])
y0 = np.array([0.0, 0.0])
t = np.arange(0., 200., 0.05)
start_time = time.time()
odeint(f, y0, t, args=(params,))
print("The Python Code took: %.6s seconds" % (time.time() - start_time))
start_time = time.time()
odeint(ode.f, y0, t, args=(params,))
print("The Cython Code took: %.6s seconds ---" % (time.time() - start_time))
start_time = time.time()
odeint(ode.fCMath, y0, t, args=(params,))
print("The Cython Code incorpoarting two of DavidW_s suggestions took: %.6s seconds ---" % (time.time() - start_time))
然后我运行:
python setup.py build_ext --inplace
python solveODE.py
在终端中。
python版本的时间约为0.055秒，
而 Cython 版本大约需要 0.04 秒。
有人对我解决问题的尝试有什么建议吗？
微分方程，最好不要用 Cython 修改 odeint 例程本身？
编辑
我将 DavidW 的建议合并到 ode.pyx 和 solveODE.py 这两个文件中 运行 使用这些建议的代码只用了大约 0.015 秒。
最简单的更改（这可能会给您带来很多好处）是使用 C 数学库 sin 和 cos 对单个数字而不是数字进行运算。调用 numpy 并确定它不是数组所花费的时间相当昂贵。
from libc.math cimport sin, cos
# later
-omega/Q + sin(theta) + d*cos(Omega*t)
我很想为输入分配一个类型 d（其他输入中的 none 无需更改界面即可轻松键入）：
def f(y, double t, params):
我想我也只是 return 一个列表，就像您在 Python 版本中所做的那样。我不认为你通过使用 C 数组获得很多。
tldr；使用 numba.jit 实现 3 倍加速...
我没有太多使用 cython 的经验，但我的机器似乎与您严格 python 版本的计算时间相似，因此我们应该能够大致比较苹果与苹果。我使用 numba 来编译函数 f（我稍微重写了它以使其在编译器中更好地发挥作用）。
def f(y, t, params):
return np.array([y[1], -y[1]/params[0] + np.sin(y[0]) + params[1]*np.cos(params[2]*t)])
numba_f = numba.jit(f)
用 numba_f 代替你的 ode.f 给我这个输出...
The Python Code took: 0.0468 seconds
The Numba Code took: 0.0155 seconds
然后我想知道我是否可以复制 odeint 并使用 numba 进行编译以进一步加快速度...（我做不到）
这是我的 Runge-Kutta 数值微分方程积分器：
#function f is provided inline (not as an arg)
def runge_kutta(y0, steps, dt, args=()): #improvement on euler's method. *note: time steps given in number of steps and dt
Y = np.empty([steps,y0.shape[0]])
Y[0] = y0
t = 0
n = 0
for n in range(steps-1):
#calculate coeficients
k1 = f(Y[n], t, args) #(euler's method coeficient) beginning of interval
k2 = f(Y[n] + (dt * k1 / 2), t + (dt/2), args) #interval midpoint A
k3 = f(Y[n] + (dt * k2 / 2), t + (dt/2), args) #interval midpoint B
k4 = f(Y[n] + dt * k3, t + dt, args) #interval end point
Y[n + 1] = Y[n] + (dt/6) * (k1 + 2*k2 + 2*k3 + k4) #calculate Y(n+1)
t += dt #calculate t(n+1)
return Y
天真的循环函数通常是编译后最快的，尽管这可能会被重新构造以获得更好的速度。我应该注意，这给出了与 odeint 不同的答案，在大约 2000 步之后偏离了 .001，并且在 3000 之后完全不同。对于函数的 numba 版本，我只是替换了 f 和 numba_f，并添加了带有 @numba.jit 的编译作为装饰器。在这种情况下，正如预期的那样，纯 python 版本非常慢，但 numba 版本并不比带有 odeint 的 numba 快（再次，ymmv）。
using custom integrator
The Python Code took: 0.2340 seconds
The Numba Code took: 0.0156 seconds
这里有一个提前编译的例子。我在这台计算机上没有必要的工具链来编译，我没有管理员安装它，所以这给了我一个错误，我没有所需的编译器，但它应该可以正常工作。
import numpy as np
from numba.pycc import CC
cc = CC('diffeq')
@cc.export('func', 'f8[:](f8[:], f8, f8[:])')
def func(y, t, params):
return np.array([y[1], -y[1]/params[0] + np.sin(y[0]) + params[1]*np.cos(params[2]*t)])
cc.compile()
如果其他人用其他模块回答这个问题，我不妨附和一下：
我是 JiTCODE, which accepts an ODE written in SymPy symbols and then converts this ODE to C code for a Python module, compiles this C code, loads the result and uses this as a derivative for SciPy’s ODE 的作者。翻译成 JiTCODE 的示例如下所示：
from jitcode import jitcode, provide_basic_symbols
import numpy as np
from sympy import sin, cos
import time
Q = 2.0
d = 1.5
Ω = 0.65
t, y = provide_basic_symbols()
f = [
y(1),
-y(1)/Q + sin(y(0)) + d*cos(Ω*t)
]
initial_state = np.array([0.0,0.0])
ODE = jitcode(f)
ODE.set_integrator("lsoda")
ODE.set_initial_value(initial_state,0.0)
start_time = time.time()
data = np.vstack(ODE.integrate(T) for T in np.arange(0.05, 200., 0.05))
end_time = time.time()
print("JiTCODE took: %.6s seconds" % (end_time - start_time))
这需要 0.11 秒，与基于 odeint 的解决方案相比非常慢，但这不是由于实际集成而是结果处理方式造成的：虽然 odeint 直接在内部有效地创建一个数组，这是通过 Python 此处完成的。根据您的操作，这可能是一个关键的缺点，但是对于较粗略的采样或较大的微分方程，这很快就会变得无关紧要。
所以，让我们删除数据收集，只看集成，将最后几行替换为以下内容：
ODE = jitcode(f)
ODE.set_integrator("lsoda", max_step=0.05, nsteps=1e10)
ODE.set_initial_value(initial_state,0.0)
start_time = time.time()
ODE.integrate(200.0)
end_time = time.time()
print("JiTCODE took: %.6s seconds" % (end_time - start_time))
请注意，我设置 max_step=0.05 以强制积分器至少执行与您的示例中一样多的步骤，并确保唯一的区别是积分结果不会存储到某个数组中。运行时间为 0.010 秒。
NumbaLSODA 需要 0.00088 秒（比 Cython 快 17 倍）。
from NumbaLSODA import lsoda_sig, lsoda
import numba as nb
import numpy as np
import time
@nb.cfunc(lsoda_sig)
def f(t, y_, dy, p_):
p = nb.carray(p_, (3,))
y = nb.carray(y_, (2,))
theta, omega = y
Q, d, Omega = p
dy[0] = omega
dy[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)
funcptr = f.address # address to ODE function
y0 = np.array([0.0, 0.0])
data = np.array([2.0, 1.5, 0.65])
t = np.arange(0., 200., 0.05)
start_time = time.time()
usol, success = lsoda(funcptr, y0, t, data = data)
print("NumbaLSODA took: %.8s seconds ---" % (time.time() - start_time))
结果
NumbaLSODA took: 0.000880 seconds ---
我想缩短 Scipy 的 odeint 求解微分的时间 方程。
为了练习，我使用了 Python in scientific computations 中的示例作为模板。因为 odeint 将函数 f 作为参数，所以我将此函数编写为静态类型的 Cython 版本并希望
odeint 的 运行ning 时间将显着减少。
函数 f 包含在名为 ode.pyx 的文件中，如下所示：
import numpy as np
cimport numpy as np
from libc.math cimport sin, cos
def f(y, t, params):
cdef double theta = y[0], omega = y[1]
cdef double Q = params[0], d = params[1], Omega = params[2]
cdef double derivs[2]
derivs[0] = omega
derivs[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)
return derivs
def fCMath(y, double t, params):
cdef double theta = y[0], omega = y[1]
cdef double Q = params[0], d = params[1], Omega = params[2]
cdef double derivs[2]
derivs[0] = omega
derivs[1] = -omega/Q + sin(theta) + d*cos(Omega*t)
return derivs
然后我创建一个文件setup.py来编译函数：
from distutils.core import setup
from Cython.Build import cythonize
setup(ext_modules=cythonize('ode.pyx'))
求解微分方程的脚本（也包含Python
f) 的版本称为 solveODE.py，看起来像：
import ode
import numpy as np
from scipy.integrate import odeint
import time
def f(y, t, params):
theta, omega = y
Q, d, Omega = params
derivs = [omega,
-omega/Q + np.sin(theta) + d*np.cos(Omega*t)]
return derivs
params = np.array([2.0, 1.5, 0.65])
y0 = np.array([0.0, 0.0])
t = np.arange(0., 200., 0.05)
start_time = time.time()
odeint(f, y0, t, args=(params,))
print("The Python Code took: %.6s seconds" % (time.time() - start_time))
start_time = time.time()
odeint(ode.f, y0, t, args=(params,))
print("The Cython Code took: %.6s seconds ---" % (time.time() - start_time))
start_time = time.time()
odeint(ode.fCMath, y0, t, args=(params,))
print("The Cython Code incorpoarting two of DavidW_s suggestions took: %.6s seconds ---" % (time.time() - start_time))
然后我运行:
python setup.py build_ext --inplace
python solveODE.py
在终端中。
python版本的时间约为0.055秒， 而 Cython 版本大约需要 0.04 秒。
有人对我解决问题的尝试有什么建议吗？ 微分方程，最好不要用 Cython 修改 odeint 例程本身？
编辑
我将 DavidW 的建议合并到 ode.pyx 和 solveODE.py 这两个文件中 运行 使用这些建议的代码只用了大约 0.015 秒。
最简单的更改（这可能会给您带来很多好处）是使用 C 数学库 sin 和 cos 对单个数字而不是数字进行运算。调用 numpy 并确定它不是数组所花费的时间相当昂贵。
from libc.math cimport sin, cos
# later
-omega/Q + sin(theta) + d*cos(Omega*t)
我很想为输入分配一个类型 d（其他输入中的 none 无需更改界面即可轻松键入）：
def f(y, double t, params):
我想我也只是 return 一个列表，就像您在 Python 版本中所做的那样。我不认为你通过使用 C 数组获得很多。
tldr；使用 numba.jit 实现 3 倍加速...
我没有太多使用 cython 的经验，但我的机器似乎与您严格 python 版本的计算时间相似，因此我们应该能够大致比较苹果与苹果。我使用 numba 来编译函数 f（我稍微重写了它以使其在编译器中更好地发挥作用）。
def f(y, t, params):
return np.array([y[1], -y[1]/params[0] + np.sin(y[0]) + params[1]*np.cos(params[2]*t)])
numba_f = numba.jit(f)
用 numba_f 代替你的 ode.f 给我这个输出...
The Python Code took: 0.0468 seconds
The Numba Code took: 0.0155 seconds
然后我想知道我是否可以复制 odeint 并使用 numba 进行编译以进一步加快速度...（我做不到）
这是我的 Runge-Kutta 数值微分方程积分器：
#function f is provided inline (not as an arg)
def runge_kutta(y0, steps, dt, args=()): #improvement on euler's method. *note: time steps given in number of steps and dt
Y = np.empty([steps,y0.shape[0]])
Y[0] = y0
t = 0
n = 0
for n in range(steps-1):
#calculate coeficients
k1 = f(Y[n], t, args) #(euler's method coeficient) beginning of interval
k2 = f(Y[n] + (dt * k1 / 2), t + (dt/2), args) #interval midpoint A
k3 = f(Y[n] + (dt * k2 / 2), t + (dt/2), args) #interval midpoint B
k4 = f(Y[n] + dt * k3, t + dt, args) #interval end point
Y[n + 1] = Y[n] + (dt/6) * (k1 + 2*k2 + 2*k3 + k4) #calculate Y(n+1)
t += dt #calculate t(n+1)
return Y
天真的循环函数通常是编译后最快的，尽管这可能会被重新构造以获得更好的速度。我应该注意，这给出了与 odeint 不同的答案，在大约 2000 步之后偏离了 .001，并且在 3000 之后完全不同。对于函数的 numba 版本，我只是替换了 f 和 numba_f，并添加了带有 @numba.jit 的编译作为装饰器。在这种情况下，正如预期的那样，纯 python 版本非常慢，但 numba 版本并不比带有 odeint 的 numba 快（再次，ymmv）。
using custom integrator
The Python Code took: 0.2340 seconds
The Numba Code took: 0.0156 seconds
这里有一个提前编译的例子。我在这台计算机上没有必要的工具链来编译，我没有管理员安装它，所以这给了我一个错误，我没有所需的编译器，但它应该可以正常工作。
import numpy as np
from numba.pycc import CC
cc = CC('diffeq')
@cc.export('func', 'f8[:](f8[:], f8, f8[:])')
def func(y, t, params):
return np.array([y[1], -y[1]/params[0] + np.sin(y[0]) + params[1]*np.cos(params[2]*t)])
cc.compile()
如果其他人用其他模块回答这个问题，我不妨附和一下：
我是 JiTCODE, which accepts an ODE written in SymPy symbols and then converts this ODE to C code for a Python module, compiles this C code, loads the result and uses this as a derivative for SciPy’s ODE 的作者。翻译成 JiTCODE 的示例如下所示：
from jitcode import jitcode, provide_basic_symbols
import numpy as np
from sympy import sin, cos
import time
Q = 2.0
d = 1.5
Ω = 0.65
t, y = provide_basic_symbols()
f = [
y(1),
-y(1)/Q + sin(y(0)) + d*cos(Ω*t)
]
initial_state = np.array([0.0,0.0])
ODE = jitcode(f)
ODE.set_integrator("lsoda")
ODE.set_initial_value(initial_state,0.0)
start_time = time.time()
data = np.vstack(ODE.integrate(T) for T in np.arange(0.05, 200., 0.05))
end_time = time.time()
print("JiTCODE took: %.6s seconds" % (end_time - start_time))
这需要 0.11 秒，与基于 odeint 的解决方案相比非常慢，但这不是由于实际集成而是结果处理方式造成的：虽然 odeint 直接在内部有效地创建一个数组，这是通过 Python 此处完成的。根据您的操作，这可能是一个关键的缺点，但是对于较粗略的采样或较大的微分方程，这很快就会变得无关紧要。
所以，让我们删除数据收集，只看集成，将最后几行替换为以下内容：
ODE = jitcode(f)
ODE.set_integrator("lsoda", max_step=0.05, nsteps=1e10)
ODE.set_initial_value(initial_state,0.0)
start_time = time.time()
ODE.integrate(200.0)
end_time = time.time()
print("JiTCODE took: %.6s seconds" % (end_time - start_time))
请注意，我设置 max_step=0.05 以强制积分器至少执行与您的示例中一样多的步骤，并确保唯一的区别是积分结果不会存储到某个数组中。运行时间为 0.010 秒。
NumbaLSODA 需要 0.00088 秒（比 Cython 快 17 倍）。
from NumbaLSODA import lsoda_sig, lsoda
import numba as nb
import numpy as np
import time
@nb.cfunc(lsoda_sig)
def f(t, y_, dy, p_):
p = nb.carray(p_, (3,))
y = nb.carray(y_, (2,))
theta, omega = y
Q, d, Omega = p
dy[0] = omega
dy[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)
funcptr = f.address # address to ODE function
y0 = np.array([0.0, 0.0])
data = np.array([2.0, 1.5, 0.65])
t = np.arange(0., 200., 0.05)
start_time = time.time()
usol, success = lsoda(funcptr, y0, t, data = data)
print("NumbaLSODA took: %.8s seconds ---" % (time.time() - start_time))
结果
NumbaLSODA took: 0.000880 seconds ---