为什么代码加速不适用于Cython？

Question

我需要将此代码加速到4毫秒。

import numpy as np

def return_call(data):
num = int(data.shape[0] / 4096)
buff_spectrum  = np.empty(2048,dtype= np.uint64)
buff_detect =  np.empty(2048,dtype= np.uint64)
end_spetrum = np.empty(num*1024,dtype=np.uint64)
end_detect = np.empty(num*1024,dtype= np.uint64)
_data = np.reshape(data,(num,4096))
for _raw_data_spec in _data:
raw_data_spec = np.reshape(_raw_data_spec,(2048,2))
for i in range(2048):
buff_spectrum[i] = (np.int16(raw_data_spec[i][0])<<17)|(np.int16(raw_data_spec[i][1] <<1))>>1
buff_detect[i] = (np.int16(raw_data_spec[i][0])>>15)
for i in range (511,-1,-1):
if buff_spectrum[i+1024] != 0:
end_spetrum[i]=(np.log10(buff_spectrum[i+1024]))
end_detect[i]=buff_detect[i+1024]
else:
end_spetrum[i] =0
end_detect[i] = 0
for i in range(1023, 511, -1):
if buff_spectrum[i+1024] != 0:
end_spetrum[i] = (np.log10(buff_spectrum[i + 1024]))
end_detect[i] = buff_detect[i + 1024]
else:
end_spetrum[i] = 0
end_detect[i] = 0
return end_spetrum, end_detect

我决定使用Cython来完成这项任务。但我没有得到任何加速。

import numpy as np
cimport numpy

ctypedef signed short DTYPE_t
cpdef return_call(numpy.ndarray[DTYPE_t, ndim=1] data):
cdef int i
cdef int num = data.shape[0]/4096
cdef numpy.ndarray _data
cdef numpy.ndarray[unsigned long long, ndim=1] buff_spectrum  = np.empty(2048,dtype= np.uint64)
cdef numpy.ndarray[ unsigned long long, ndim=1] buff_detect =  np.empty(2048,dtype= np.uint64)
cdef numpy.ndarray[double , ndim=1] end_spetrum = np.empty(num*1024,dtype= np.double)
cdef numpy.ndarray[double , ndim=1] end_detect = np.empty(num*1024,dtype= np.double)
_data = np.reshape(data,(num,4096))
for _raw_data_spec in _data:
raw_data_spec = np.reshape(_raw_data_spec,(2048,2))
for i in range(2048):
buff_spectrum[i] = (np.uint16(raw_data_spec[i][0])<<17)|(np.uint16(raw_data_spec[i][1] <<1))>>1
buff_detect[i] = (np.uint16(raw_data_spec[i][0])>>15)
for i in range (511,-1,-1):
if buff_spectrum[i+1024] != 0:
end_spetrum[i]=(np.log10(buff_spectrum[i+1024]))
end_detect[i]=buff_detect[i+1024]
else:
end_spetrum[i] =0
end_detect[i] = 0
for i in range(1023, 511, -1):
if buff_spectrum[i+1024] != 0:
end_spetrum[i] = (np.log10(buff_spectrum[i + 1024]))
end_detect[i] = buff_detect[i + 1024]
else:
end_spetrum[i] = 0
end_detect[i] = 0
return end_spetrum, end_detect

我达到的最大速度是80毫秒，但我需要更快。由于您需要几乎实时地处理铁的数据告诉我原因。实现预期结果是否现实。我还附上了测试文件的代码。

import numpy as np
import example_original
import example_cython
data = np.empty(8192*2, dtype=np.int16)
import time
startpy = time.time()

example_original.return_call(data)
finpy = time.time() -startpy
startcy = time.time()
k,r = example_cython.return_call(data)
fincy = time.time() -startcy
print( fincy, finpy)
print('Cython is {}x faster'.format(finpy/fincy))

Answer 1

我认为这可能是因为您的python代码几乎没有python操作，而且都是numpy操作。numpy代码的很大一部分是用C编写的。其中一些是用Fortran编写的。很多都是用Python编写的。写得好的numpy代码在速度上与C代码相当。

Answer 2

raw_data_spec = np.reshape(_raw_data_spec,(2048,2))

未键入raw_data_spec。在函数的开头添加一个定义。我推荐更新的memoryview语法(但如果你愿意，可以使用旧的numpy语法(：

cdef DTYPE_t[:,:] raw_data_spec

---

这条线(你已经确定为瓶颈(一团糟：

buff_spectrum[i] = (np.int16(raw_data_spec[i][0])<<17)|(np.int16(raw_data_spec[i][1] <<1))>>1

1. 在一个步骤中进行索引，而不是两个步骤：raw_data_spec[i, 0](注意一个大括号和一个逗号(。

2. 重新考虑转换为16位整数。将16位整数移位17位真的有意义吗？

3. 您可能根本不需要强制转换，因为数据已知为DTYPE_t，但如果您确实想要强制转换，则使用尖括号：<numpy.uint16_t>(raw_data_spec[i, 0])

---

考虑关闭boundscheck和wraparound请自己验证这样做是安全的，并且当索引超出数组末尾或使用负索引时，您不会依赖异常来告诉您。只有经过思考才能做到这一点——而不是自动地以"货物崇拜"的方式。

cimport cython    
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef return_call(numpy.ndarray[DTYPE_t, ndim=1] data):

---

放弃对np.log10的调用。这是对单个元素的一个完整的Python调用，最终效率低下。您可以使用C标准库的数学函数：

from libc.math cimport log10

则用CCD_ 9代替CCD_。

Answer 3

我对Cython没有太多经验，所以这只是Cython计时的一个例子。

示例

import numpy as np
import numba as nb
@nb.njit(cache=True)
def return_call(data_in):
#If the input is not contigous the reshape will fail
#-> make a c-contigous copy if the array isn't c-contigous
data=np.ascontiguousarray(data_in)
num = int(data.shape[0] / 4096)
buff_spectrum  = np.zeros(2048,dtype= np.uint64)
buff_detect =  np.zeros(2048,dtype= np.uint64)
end_spetrum = np.zeros(num*1024,dtype=np.float64)
end_detect = np.zeros(num*1024,dtype= np.float64)
_data = np.reshape(data,(num,4096))
#for _raw_data_spec in _data: is not supported
#but the followin works
for x in range(_data.shape[0]):
raw_data_spec = np.reshape(_data[x],(2048,2))
for i in range(2048):
buff_spectrum[i] = (np.int16(raw_data_spec[i][0])<<17)|(np.int16(raw_data_spec[i][1] <<1))>>1
buff_detect[i] = (np.int16(raw_data_spec[i][0])>>15)
for i in range (511,-1,-1):
if buff_spectrum[i+1024] != 0:
end_spetrum[i]=(np.log10(buff_spectrum[i+1024]))
end_detect[i]=buff_detect[i+1024]
else:
end_spetrum[i] =0
end_detect[i] = 0
for i in range(1023, 511, -1):
if buff_spectrum[i+1024] != 0:
end_spetrum[i] = (np.log10(buff_spectrum[i + 1024]))
end_detect[i] = buff_detect[i + 1024]
else:
end_spetrum[i] = 0
end_detect[i] = 0
return end_spetrum, end_detect

计时

data = np.random.rand(8192*2)*20
data=data.astype(np.int16)
#with compilation
%timeit end_spetrum, end_detect=return_call(data)
#32.7 µs ± 5.61 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)
#without compilation
%timeit end_spetrum, end_detect=return_call_orig(data)
#106 ms ± 448 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)