【发布时间】：2015-06-09 17:24:50
【问题描述】：
我是第一次尝试 cython。并尝试将函数从使用纯 numpy 转换为 cython
这是两个函数：
from __future__ import division
import numpy as np
cimport numpy as np
DTYPEf = np.float64
ctypedef np.float64_t DTYPEf_t
DTYPEi = np.int64
ctypedef np.int64_t DTYPEi_t
DTYPEu = np.uint8
ctypedef np.uint8_t DTYPEu_t
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
def twodcitera(np.ndarray[DTYPEf_t, ndim=3] data, int res, int indexl, int indexu, float radius1, float radius2, output, float height1, float height2 ):
'''
Function to return correlation for fixed radius using Cython
'''
cdef float sum_mask = 0
cdef int i,j,k
cdef int a, b, c
cdef np.ndarray[DTYPEi_t, ndim=3] x
cdef np.ndarray[DTYPEi_t, ndim=3] y
cdef np.ndarray[DTYPEi_t, ndim=3] z
cdef np.ndarray[DTYPEu_t, ndim=3, cast=True] R
a,b,c = res//2,res//2,res//2
x,y,z = np.ogrid[-a:a,-b:b,-c:c]
for i in xrange(indexl,indexu):
for j in xrange(1):
for k in xrange(1):
R = np.roll(np.roll(np.roll(np.logical_and(np.logical_or(np.logical_and(z>height1,z<=height2), np.logical_and(z<-height1,z>=-height2)), np.logical_and(x**2 + y**2<= radius2**2, x**2 + y**2 > radius1**2)), (i-a), axis =0), (j-a), axis =1), (k-a), axis =2)
sum_mask += (data[i][j][k] * np.average(data[R]))
output.put(sum_mask)
对于 numpy 实现：
def no_twodcitera(data, res, indexl, indexu, radius1, radius2, output, height1, height2 ):
'''
Function to return correlation for fixed radius
'''
a,b,c = res/2,res/2,res/2
x,y,z = np.ogrid[-a:a,-b:b,-c:c]
sum_mask = 0
for i in xrange(indexl,indexu):
for j in xrange(1):
for k in xrange(1):
R = np.roll(np.roll(np.roll(np.logical_and(np.logical_or(np.logical_and(z>height1,z<=height2), np.logical_and(z<-height1,z>=-height2)), np.logical_and(x**2 + y**2<= radius2**2, x**2 + y**2 > radius1**2)), (i-a), axis =0), (j-a), axis =1), (k-a), axis =2)
sum_mask += (data[i][j][k] * np.average(data[R]))
output.put(sum_mask)
这两个函数实际上给了我相同的完成时间。
%timeit -n200 -r10 twodcitera(dd, tes_res,in1,in2,r[k],r[k+1], output, r[l], r[l+1])
200 loops, best of 10: 1.57 ms per loop
%timeit -n200 -r10 no_twodcitera(dd, tes_res,in1,in2,r[k],r[k+1], output, r[l], r[l+1])
200 loops, best of 10: 1.57 ms per loop
我想知道我做错了什么，或者在尝试实现 cython 时我没有正确理解。输入是：
dd = np.random.randn(64,64,64)
res = 64
r = np.arange(0,21,2)
in1 = 0
in2 = 1
l = 5
k = 7
output = mp.Queue()
谢谢你在这里指出我的误解。
【问题讨论】：
-
j 和 k 的 xrange 保留为 1 仅用于测试目的，最终它将是 xrange(res) 中的 j 和 xrange(res) 中的 k
-
您是否尝试使用 cython -a 运行代码？ docs.cython.org/src/quickstart/…
-
还有什么是 in1,in2 等..
-
所以 r = np.arange(0,21,2)。 in1 和 in2 以及遍历数组的两个索引，例如 in1 = 0，in2 =1。 Radius1 和 radius2 是圆形壳的半径，而 height1 和 height2 是圆柱壳的长度。我已将代码并行化以在多个内核上运行，因此使用了输出，即 output=mp.Queue()。希望这有助于更好地理解代码。
-
np.roll,np.logical_and等：所有这些都是慢速 python 函数，如果将它们转换为 cython，则不会有太大的加速。如果您真的想要加速，请将其表达为对 R 数组的每个索引进行循环以执行相同的操作。此外，data[i][j][k]对 Cython 速度来说很糟糕，请改用data[i,j,k]。
标签： python numpy optimization cython