Cython 代码比 Python / Numpy 代码慢 3-4 倍？答案

【问题标题】：Cython code 3-4 times slower than Python / Numpy code?Cython 代码比 Python / Numpy 代码慢 3-4 倍？
【发布时间】：2012-10-14 05:10:34
【问题描述】：

我正在尝试将我的 Python / Numpy 代码转换为 Cython 代码以加快速度。但是，Cython 比 Python / Numpy 代码慢得多（3-4 倍）。我是否正确使用 Cython？我是否在 Cython 代码中将参数正确传递给 myc_rb_etc() ？当我调用集成函数时呢？预先感谢您的帮助。这是我的 Python / Numpy 代码：

from pylab import * 
import pylab as pl
from numpy import *
import numpy as np
from scipy import integrate

def myc_rb_e2f(y,t,k,d):

    M = y[0]
    E = y[1]
    CD = y[2]
    CE = y[3]
    R = y[4]
    RP = y[5] 
    RE = y[6]

    S = 0.01
    if t > 300:
        S = 5.0
    #if t > 400
        #S = 0.01

    t1 = k[0]*S/(k[7]+S);
    t2 = k[1]*(M/(k[14]+M))*(E/(k[15]+E));
    t3 = k[5]*M/(k[14]+M);
    t4 = k[11]*CD*RE/(k[16]+RE);
    t5 = k[12]*CE*RE/(k[17]+RE);
    t6 = k[2]*M/(k[14]+M);
    t7 = k[3]*S/(k[7]+S);
    t8 = k[6]*E/(k[15]+E);
    t9 = k[13]*RP/(k[18]+RP);
    t10 = k[9]*CD*R/(k[16]+R);
    t11 = k[10]*CE*R/(k[17]+R);

    dM = t1-d[0]*M
    dE = t2+t3+t4+t5-k[8]*R*E-d[1]*E
    dCD = t6+t7-d[2]*CD
    dCE = t8-d[3]*CE
    dR = k[4]+t9-k[8]*R*E-t10-t11-d[4]*R
    dRP = t10+t11+t4+t5-t9-d[5]*RP
    dRE = k[8]*R*E-t4-t5-d[6]*RE

    dy = [dM,dE,dCD,dCE,dR,dRP,dRE]

    return dy

t = np.zeros(10000)
t = np.linspace(0.,3000.,10000.)

# Initial concentrations of [M,E,CD,CE,R,RP,RE]
y0 = np.array([0.,0.,0.,0.,0.4,0.,0.25])
E_simulated = np.zeros([10000,5000])
E_avg = np.zeros([10000])
k = np.zeros([19])
d = np.zeros([7])

for i in range (0,5000):
    k[0] = 1.+0.1*randn(1)
    k[1] = 0.15+0.05*randn(1)
    k[2] = 0.2+0.05*randn(1)
    k[3] = 0.2+0.05*randn(1)
    k[4] = 0.35+0.05*randn(1)
    k[5] = 0.001+0.0001*randn(1)
    k[6] = 0.5+0.05*randn(1)
    k[7] = 0.3+0.05*randn(1)
    k[8] = 30.+5.*randn(1)
    k[9] = 18.+3.*randn(1)
    k[10] = 18.+3.*randn(1)
    k[11] = 18.+3.*randn(1)
    k[12] = 18.+3.*randn(1)
    k[13] = 3.6+0.5*randn(1)
    k[14] = 0.15+0.05*randn(1)
    k[15] = 0.15+0.05*randn(1)
    k[16] = 0.92+0.1*randn(1)
    k[17] = 0.92+0.1*randn(1)
    k[18] = 0.01+0.001*randn(1)
    d[0] = 0.7+0.05*randn(1)
    d[1] = 0.25+0.025*randn(1)
    d[2] = 1.5+0.05*randn(1)
    d[3] = 1.5+0.05*randn(1)
    d[4] = 0.06+0.01*randn(1)
    d[5] = 0.06+0.01*randn(1)
    d[6] = 0.03+0.005*randn(1)
    r = integrate.odeint(myc_rb_e2f,y0,t,args=(k,d))
    E_simulated[:,i] = r[:,1]

for i in range(0,10000):
    E_avg[i] = sum(E_simulated[i,:])/5000.

pl.plot(t,E_avg,'-ro')
pl.show()

下面是转换成 Cython 的代码：

cimport numpy as np
import numpy as np
from numpy import *
import pylab as pl
from pylab import * 
from scipy import integrate

def myc_rb_e2f(y,t,k,d):

    cdef double M = y[0]
    cdef double E = y[1]
    cdef double CD = y[2]
    cdef double CE = y[3]
    cdef double R = y[4]
    cdef double RP = y[5] 
    cdef double RE = y[6]

    cdef double S = 0.01
    if t > 300.0:
        S = 5.0
    #if t > 400
        #S = 0.01

    cdef double t1 = k[0]*S/(k[7]+S)
    cdef double t2 = k[1]*(M/(k[14]+M))*(E/(k[15]+E))
    cdef double t3 = k[5]*M/(k[14]+M)
    cdef double t4 = k[11]*CD*RE/(k[16]+RE)
    cdef double t5 = k[12]*CE*RE/(k[17]+RE)
    cdef double t6 = k[2]*M/(k[14]+M)
    cdef double t7 = k[3]*S/(k[7]+S)
    cdef double t8 = k[6]*E/(k[15]+E)
    cdef double t9 = k[13]*RP/(k[18]+RP)
    cdef double t10 = k[9]*CD*R/(k[16]+R)
    cdef double t11 = k[10]*CE*R/(k[17]+R)

    cdef double dM = t1-d[0]*M
    cdef double dE = t2+t3+t4+t5-k[8]*R*E-d[1]*E
    cdef double dCD = t6+t7-d[2]*CD
    cdef double dCE = t8-d[3]*CE
    cdef double dR = k[4]+t9-k[8]*R*E-t10-t11-d[4]*R
    cdef double dRP = t10+t11+t4+t5-t9-d[5]*RP
    cdef double dRE = k[8]*R*E-t4-t5-d[6]*RE

    dy = [dM,dE,dCD,dCE,dR,dRP,dRE]

    return dy


def main():
    cdef np.ndarray[double,ndim=1] t = np.zeros(10000)
    t = np.linspace(0.,3000.,10000.)
    # Initial concentrations of [M,E,CD,CE,R,RP,RE]
    cdef np.ndarray[double,ndim=1] y0 = np.array([0.,0.,0.,0.,0.4,0.,0.25])
    cdef np.ndarray[double,ndim=2] E_simulated = np.zeros([10000,5000])
    cdef np.ndarray[double,ndim=2] r = np.zeros([10000,7])
    cdef np.ndarray[double,ndim=1] E_avg = np.zeros([10000])
    cdef np.ndarray[double,ndim=1] k = np.zeros([19])
    cdef np.ndarray[double,ndim=1] d = np.zeros([7])
    cdef int i
    for i in range (0,5000):
        k[0] = 1.+0.1*randn(1)
        k[1] = 0.15+0.05*randn(1)
        k[2] = 0.2+0.05*randn(1)
        k[3] = 0.2+0.05*randn(1)
        k[4] = 0.35+0.05*randn(1)
        k[5] = 0.001+0.0001*randn(1)
        k[6] = 0.5+0.05*randn(1)
        k[7] = 0.3+0.05*randn(1)
        k[8] = 30.+5.*randn(1)
        k[9] = 18.+3.*randn(1)
        k[10] = 18.+3.*randn(1)
        k[11] = 18.+3.*randn(1)
        k[12] = 18.+3.*randn(1)
        k[13] = 3.6+0.5*randn(1)
        k[14] = 0.15+0.05*randn(1)
        k[15] = 0.15+0.05*randn(1)
        k[16] = 0.92+0.1*randn(1)
        k[17] = 0.92+0.1*randn(1)
        k[18] = 0.01+0.001*randn(1)
        d[0] = 0.7+0.05*randn(1)
        d[1] = 0.25+0.025*randn(1)
        d[2] = 1.5+0.05*randn(1)
        d[3] = 1.5+0.05*randn(1)
        d[4] = 0.06+0.01*randn(1)
        d[5] = 0.06+0.01*randn(1)
        d[6] = 0.03+0.005*randn(1)
        r = integrate.odeint(myc_rb_e2f,y0,t,args=(k,d))
        E_simulated[:,i] = r[:,1]
    for i in range(0,10000):
        E_avg[i] = sum(E_simulated[i,:])/5000.
    pl.plot(t,E_avg,'-ro')
    pl.show()

以下是我的 Python / Numpy 代码中 cProfile 中的一些 pstats：

ncalls tottime percall cumtime percall

5000 82.505 0.017 236.760 0.047 {scipy.integrate._odepack.odeint}

1 1.504 1.504 238.949 238.949 myc_rb_e2f.py:1(<module>)

5000 0.025 0.000 236.855 0.047 C:\Python27\lib\site-packages\scipy\integrate\odepack.py:18(odeint)

12291237 154.255 0.000 154.255 0.000 myc_rb_e2f.py:7(myc_rb_e2f)

【问题讨论】：

分析器对时间花在哪里说了些什么？
@Zhenya - 我以前从未分析过。你能告诉我该怎么做吗？
最好将系数系列1., 0.15, 0.2, ...放入一个numpy数组（a），将0.1, 0.05, 0.05, ...放入b，然后使用k = a + b * numpy.random.randn(18)。 d 也是如此。
@Zack 要分析您的代码，您可以使用python -m cProfile -o profile.dat <your-script> 运行它，然后使用python -m pstats profile.dat 后跟命令stats 开始了解配置文件数据。要查看 pstats 中的其他命令帮助，请使用命令 help。
我在脚本上运行了 cProfile -- 我复制了一些我认为对你们有用的输出行：1 1.530 1.530 240.323 240.323 myc_rb_e2f.py:1(<module>)12267587 154.889 0.00 154.88 0.00 myc_rb_e2f.py:7(myc_rb_e2f)5000 0.024 0.000 238.539 0.048 odepack.py:18(odeint)5000 83.554 0.017 238.443 0.048{scipy.integrate._odepack.odeint} 对我来说，它看起来就像我唯一需要时间的是集成电话。请记住，这是 Python / Numpy 代码的配置文件，而不是 Cython 代码。

标签： python performance numpy scipy cython

【解决方案1】：

更改函数定义以包含参数的类型：

def myc_rb_e2f(np.ndarray[double,ndim=1]y, double t, np.ndarray[double, ndim=1] k, np.ndarray[double, ndim=1] d):

这将使运行时间比 numpy 实现提高约 3 倍，比初始 cython 实现提高 6-7 倍。仅供参考，我降低了迭代次数，这样我就不必在测试时永远等待它完成。它应该随着您期望的迭代次数而扩大。

[pkerp@plastilin so]$ time python run_numpy.py

real    0m47.572s
user    0m45.702s
sys     0m0.049s

[pkerp@plastilin so]$ time python run_cython1.py

real    1m14.851s
user    1m12.308s
sys     0m0.135s

[pkerp@plastilin so]$ time python run_cython2.py

real    0m15.774s
user    0m14.115s
sys     0m0.105s

编辑：

此外，您不必每次想要从myc_rb_e2f 返回结果时都创建一个新数组。您可以在main 中声明一个结果数组，在每次调用时将其传入，然后填写。这样可以为您节省大量不必要的分配。这比之前的最佳运行时间减半：

[pkerp@plastilin so]$ time python run_cython3.py

real    0m6.165s
user    0m4.818s
sys     0m0.152s

还有代码：

cimport numpy as np
import numpy as np
from numpy import *
import pylab as pl
from pylab import * 
from scipy import integrate

def myc_rb_e2f(np.ndarray[double,ndim=1]y, double t, np.ndarray[double, ndim=1] k, np.ndarray[double, ndim=1] d, np.ndarray[double, ndim=1] res):

    cdef double S = 0.01
    if t > 300.0:
        S = 5.0
    #if t > 400
        #S = 0.01

    cdef double t1 = k[0]*S/(k[7]+S)
    cdef double t2 = k[1]*(y[0]/(k[14]+y[0]))*(y[1]/(k[15]+y[1]))
    cdef double t3 = k[5]*y[0]/(k[14]+y[0])
    cdef double t4 = k[11]*y[2]*y[6]/(k[16]+y[6])
    cdef double t5 = k[12]*y[3]*y[6]/(k[17]+y[6])
    cdef double t6 = k[2]*y[0]/(k[14]+y[0])
    cdef double t7 = k[3]*S/(k[7]+S)
    cdef double t8 = k[6]*y[1]/(k[15]+y[1])
    cdef double t9 = k[13]*y[5]/(k[18]+y[5])
    cdef double t10 = k[9]*y[2]*y[4]/(k[16]+y[4])
    cdef double t11 = k[10]*y[3]*y[4]/(k[17]+y[4])

    cdef double dM = t1-d[0]*y[0]
    cdef double dE = t2+t3+t4+t5-k[8]*y[4]*y[1]-d[1]*y[1]
    cdef double dCD = t6+t7-d[2]*y[2]
    cdef double dCE = t8-d[3]*y[3]
    cdef double dR = k[4]+t9-k[8]*y[4]*y[1]-t10-t11-d[4]*y[4]
    cdef double dRP = t10+t11+t4+t5-t9-d[5]*y[5]
    cdef double dRE = k[8]*y[4]*y[1]-t4-t5-d[6]*y[6]

    res[0] = dM
    res[1] = dE
    res[2] = dCD
    res[3] = dCE
    res[4] = dR
    res[5] = dRP
    res[6] = dRE

    return res


def main():
    cdef np.ndarray[double,ndim=1] t = np.zeros(467)
    cdef np.ndarray[double,ndim=1] results = np.zeros(7)
    t = np.linspace(0.,3000.,467.)
    # Initial concentrations of [M,E,CD,CE,R,RP,RE]
    cdef np.ndarray[double,ndim=1] y0 = np.array([0.,0.,0.,0.,0.4,0.,0.25])
    cdef np.ndarray[double,ndim=2] E_simulated = np.zeros([467,554])
    cdef np.ndarray[double,ndim=2] r = np.zeros([467,7])
    cdef np.ndarray[double,ndim=1] E_avg = np.zeros([467])
    cdef np.ndarray[double,ndim=1] k = np.zeros([19])
    cdef np.ndarray[double,ndim=1] d = np.zeros([7])
    cdef int i
    for i in range (0,554):
        k[0] = 1.+0.1*randn(1)
        k[1] = 0.15+0.05*randn(1)
        k[2] = 0.2+0.05*randn(1)
        k[3] = 0.2+0.05*randn(1)
        k[4] = 0.35+0.05*randn(1)
        k[5] = 0.001+0.0001*randn(1)
        k[6] = 0.5+0.05*randn(1)
        k[7] = 0.3+0.05*randn(1)
        k[8] = 30.+5.*randn(1)
        k[9] = 18.+3.*randn(1)
        k[10] = 18.+3.*randn(1)
        k[11] = 18.+3.*randn(1)
        k[12] = 18.+3.*randn(1)
        k[13] = 3.6+0.5*randn(1)
        k[14] = 0.15+0.05*randn(1)
        k[15] = 0.15+0.05*randn(1)
        k[16] = 0.92+0.1*randn(1)
        k[17] = 0.92+0.1*randn(1)
        k[18] = 0.01+0.001*randn(1)
        d[0] = 0.7+0.05*randn(1)
        d[1] = 0.25+0.025*randn(1)
        d[2] = 1.5+0.05*randn(1)
        d[3] = 1.5+0.05*randn(1)
        d[4] = 0.06+0.01*randn(1)
        d[5] = 0.06+0.01*randn(1)
        d[6] = 0.03+0.005*randn(1)
        r = integrate.odeint(myc_rb_e2f,y0,t,args=(k,d,results))
        E_simulated[:,i] = r[:,1]
    for i in range(0,467):
        E_avg[i] = sum(E_simulated[i,:])/554.
    #pl.plot(t,E_avg,'-ro')
    #pl.show()

if __name__ == "__main__":
    main()

【讨论】：

首先，非常感谢。在我发布我的问题之前，我尝试将参数的类型包括到myc_rb_e2f()，但我得到了不好的结果，因为我传递的是np.ndarray[double,ndim=1] t而不是double t。更改为 double t 给了我正确的结果 - 为什么将 t 作为双精度而不是数组传递？此外，您的加速也有很大帮助。随着您的编辑，我的代码从 100 秒变为 40 秒。
好问题。我实际上并没有真正看你的程序是做什么的。我刚刚看到您进行了if t > 300.0 检查并假设它是双重的。既然您提到了它，它似乎确实被声明为一个数组。恐怕我无法帮助您解决这种差异，因为这取决于您实际尝试做什么。
另外，你是如何编译你的 cython 代码的？您是否包含优化标志？这些可以进一步改善它。
这段代码所做的就是使用随机生成的参数k 和d 生成合成数据。我多次求解 ODE 的系统，并取 dE 的时间过程的平均值 E_avg 并绘制它。我知道情节应该是什么样子，当我使用double t 作为参数时，它看起来很完美。然而，当我试图将t 作为数组传递给myc_rb_e2f() 时，我得到了错误的情节。我有一个 setup.py 文件。所以从 IPython，我做：%run setup.py build_ext --inplace;import myscript;myscript.main(); 编译并运行我的代码。
是的，这是因为 odeint 函数会针对 t 中的每个点评估您的 myc_rb_e2f 函数。因此它会遍历这些值并将双精度值传递给您的 myc_rb_e2f 函数。

【解决方案2】：

免责声明：我是子提及工具的核心开发人员之一

作为替代方案，将myc_rb_e2f 函数提取到单个文件并使用pythran 编译它会产生x5 的加速。

添加了提取的python代码（只有一个类型签名）：

#pythran export myc_rb_e2f(float[], float, float[], float[])
def myc_rb_e2f(y,t,k,d):

    M = y[0]
    E = y[1]
    CD = y[2]
    CE = y[3]
    R = y[4]
    RP = y[5] 
    RE = y[6]

    S = 0.01
    if t > 300:
        S = 5.0
    #if t > 400
        #S = 0.01

    t1 = k[0]*S/(k[7]+S);
    t2 = k[1]*(M/(k[14]+M))*(E/(k[15]+E));
    t3 = k[5]*M/(k[14]+M);
    t4 = k[11]*CD*RE/(k[16]+RE);
    t5 = k[12]*CE*RE/(k[17]+RE);
    t6 = k[2]*M/(k[14]+M);
    t7 = k[3]*S/(k[7]+S);
    t8 = k[6]*E/(k[15]+E);
    t9 = k[13]*RP/(k[18]+RP);
    t10 = k[9]*CD*R/(k[16]+R);
    t11 = k[10]*CE*R/(k[17]+R);

    dM = t1-d[0]*M
    dE = t2+t3+t4+t5-k[8]*R*E-d[1]*E
    dCD = t6+t7-d[2]*CD
    dCE = t8-d[3]*CE
    dR = k[4]+t9-k[8]*R*E-t10-t11-d[4]*R
    dRP = t10+t11+t4+t5-t9-d[5]*RP
    dRE = k[8]*R*E-t4-t5-d[6]*RE

    dy = [dM,dE,dCD,dCE,dR,dRP,dRE]

    return dy

编译：

> pythran myc.py

并用from myc import myc_rb_e2f 代替函数定义导入。

原始代码在我笔记本电脑上的103.46s 中运行，调用pythran 编译例程的代码在22.23s 中运行。

【讨论】：