优化我的Cython/Numpy代码？至今仅提高了30%的性能

Question

我是不是还有什么事情没做，能让这个过程快一点？我正在尝试实现一本叫《调音音色谱尺度》的书里描述的算法。另外，如果其他方法都不行，我能不能把这部分代码用C语言写，然后在Python里调用它？

import numpy as np
cimport numpy as np

# DTYPE = np.float
ctypedef np.float_t DTYPE_t

np.seterr(divide='raise', over='raise', under='ignore', invalid='raise')

"""
I define a timbre as the following 2d numpy array:
[[f0, a0], [f1, a1], [f2, a2]...] where f describes the frequency
of the given partial and a is its amplitude from 0 to 1. Phase is ignored.
"""

#Test Timbre
# cdef np.ndarray[DTYPE_t,ndim=2] t1 = np.array( [[440,1],[880,.5],[(440*3),.333]])

# Calculates the inherent dissonance of one timbres of the above form
# using the diss2Partials function
cdef DTYPE_t diss1Timbre(np.ndarray[DTYPE_t,ndim=2] t):
    cdef DTYPE_t runningDiss1
    runningDiss1 = 0.0
    cdef unsigned int len = np.shape(t)[0]
    cdef unsigned int i
    cdef unsigned int j
    for i from 0 <= i < len:
        for j from i+1 <= j < len:
            runningDiss1 += diss2Partials(t[i], t[j])
    return runningDiss1

# Calculates the dissonance between two timbres of the above form 
cdef DTYPE_t diss2Timbres(np.ndarray[DTYPE_t,ndim=2] t1, np.ndarray[DTYPE_t,ndim=2] t2):
    cdef DTYPE_t runningDiss2
    runningDiss2 = 0.0
    cdef unsigned int len1 = np.shape(t1)[0]
    cdef unsigned int len2 = np.shape(t2)[0]
    runningDiss2 += diss1Timbre(t1)
    runningDiss2 += diss1Timbre(t2)
    cdef unsigned int i1
    cdef unsigned int i2
    for i1 from 0 <= i1 < len1:
        for i2 from 0 <= i2 < len2:
            runningDiss2 += diss2Partials(t1[i1], t2[i2])
    return runningDiss2

cdef inline DTYPE_t float_min(DTYPE_t a, DTYPE_t b): return a if a <= b else b

# Calculates the dissonance of two partials of the form [f,a]
cdef DTYPE_t diss2Partials(np.ndarray[DTYPE_t,ndim=1] p1, np.ndarray[DTYPE_t,ndim=1] p2):
    cdef DTYPE_t f1 = p1[0]
    cdef DTYPE_t f2 = p2[0]
    cdef DTYPE_t a1 = abs(p1[1])
    cdef DTYPE_t a2 = abs(p2[1])

    # In order to insure that f2 > f1:
    if (f2 < f1):
        (f1,f2,a1,a2) = (f2,f1,a2,a1)

    # Constants of the dissonance curves
    cdef DTYPE_t _xStar
    _xStar = 0.24
    cdef DTYPE_t _s1
    _s1 = 0.021
    cdef DTYPE_t _s2
    _s2 = 19
    cdef DTYPE_t _b1
    _b1 = 3.5
    cdef DTYPE_t _b2
    _b2 = 5.75

    cdef DTYPE_t a = float_min(a1,a2)
    cdef DTYPE_t s = _xStar/(_s1*f1 + _s2)
    return (a * (np.exp(-_b1*s*(f2-f1)) - np.exp(-_b2*s*(f2-f1)) ) )

cpdef dissTimbreScale(np.ndarray[DTYPE_t,ndim=2] t,np.ndarray[DTYPE_t,ndim=1] s):
    cdef DTYPE_t currDiss
    currDiss = 0.0;
    cdef unsigned int i
    for i from 0 <= i < s.size:
        currDiss += diss2Timbres(t, transpose(t,s[i]))
    return currDiss

cdef np.ndarray[DTYPE_t,ndim=2] transpose(np.ndarray[DTYPE_t,ndim=2] t, DTYPE_t ratio):
    return np.dot(t, np.array([[ratio,0],[0,1]]))

代码链接：Cython代码

Answer 1

我建议你先分析一下你的代码，看看哪个函数耗时最长。如果是 diss2Timbres，那么你可以试试“numexpr”这个包，它可能会对你有帮助。

我曾经比较过 Python/Cython 和 Numexpr 在我一个函数上的表现（链接到 SO）。根据数组的大小，numexpr 的表现比 Cython 和 Fortran 都要好。

注意：我刚发现这篇帖子其实很旧了……

Answer 2

在你的代码中：

for i from 0 <= i < len:
    for j from i+1 <= j < len:
        runningDiss1 += diss2Partials(t[i], t[j])
return runningDiss1

每次查找数组时都会进行边界检查。你可以在函数前面加上这个装饰器 @cython.boundscheck(False)，这样就可以关闭这个检查。然后在使用 i 和 j 作为索引之前，把它们转换成无符号整数类型。想了解更多信息，可以查看 Cython与Numpy的教程。

Answer 3

我注意到了一些事情：

1. 用 t1.shape[0] 替代 np.shape(t1)[0]，在其他地方也一样。
2. 不要把 len 当作变量名，因为它是Python的内置函数（这不是为了速度，而是为了良好的编程习惯）。可以用 L 或其他类似的名字。
3. 除非真的需要，不要把两个元素的数组传给函数。因为每次传数组时，Cython都会检查这个数组的缓冲区。所以，当你用 diss2Partials(t[i], t[j]) 时，应该改成 diss2Partials(t[i,0], t[i,1], t[j,0], t[j,1])，并相应地重新定义 diss2Partials。
4. 不要使用 abs，或者至少不要用Python的那个。因为它需要把你的C语言的双精度浮点数转换成Python的浮点数，调用abs函数，然后再转换回C的双精度浮点数。可能更好的是像你做的 float_min 那样，创建一个内联函数。
5. 调用 np.exp 的过程和使用 abs 类似。把 np.exp 改成 exp，并在文件顶部添加 from libc.math cimport exp。
6. 完全去掉 transpose 函数。 np.dot 会让程序变慢，但这里其实并不需要矩阵乘法。重写你的 dissTimbreScale 函数，创建一个空矩阵，比如叫 t2。在当前循环之前，把 t2 的第二列设置为 t 的第二列（最好用循环，但你也可以用Numpy操作）。然后在当前循环内部，添加一个循环，把 t2 的第一列设置为 t 的第一列乘以 s[i]。这就是你之前矩阵乘法的实际作用。最后，把 t2 作为第二个参数传给 diss2Timbres，而不是用 transpose 函数返回的那个。

先做1到5，因为这些比较简单。第6个可能需要更多时间、精力和一些实验，但我怀疑这会显著提高你的速度。