Cython speed vs numpy

2019-07-21 20:35发布

站内问答 / Python
1231 1 6

I am trying out cython for first time. And tried to convert a function from using pure numpy to cython

Here are the two functions:

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)

And for numpy implementation:

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)

Both functions actually give me same time for completion.

%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

I was wondering what I was doing wrong or I haven't understood correctly when trying to implement cython. Inputs are:

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()

Thank you if you could point out my misunderstanding here.

1条回答
姐就是有狂的资本
2楼-- · 2019-07-21 21:11

Without knowing your input and output the following compiled for me following the cython guide If you explain how to create test input i could maybe provide more help.

EDIT: My first thought was that there is maybe something up with the cython compilation. But i could not find anything realy helpful. Therfore this answer is not really helpful for improving the speed problem. Anyway i leave it here for those interested in the testing and understanding.

Put the code into test.pyx

cimport cython
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


@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)

Create a make file setup.py and put

from distutils.core import setup
from Cython.Build import cythonize

setup(
    name = "testapp",
    ext_modules = cythonize('test.pyx'),  # accepts a glob pattern
    )

Go to the shell and compile it:

$python setup.py build_ext --inplace

Go to ipython and try to import:

from test import *

worked for me to run.

A speed test showed:

In [28]: %timeit -n200 -r10 no_twodcitera(dd, res,in1,in2,r[k],r[k+1], output, r[l], r[l+1])
200 loops, best of 10: 1.29 ms per loop

In [29]: %timeit -n200 -r10 test.twodcitera(dd, res,in1,in2,r[k],r[k+1], output, r[l], r[l+1])
200 loops, best of 10: 1.31 ms per loop

So the results are the same and there is not much difference. I furthermore conducted a cProfile study to see if there is something showing up in the runtime of the call stack. One has to confess, that cProfile gets hard to interpret when in comes to ms seconds speed! But lets give it a try.

In [34]: cProfile.run("""no_twodcitera(dd, res,in1,in2,r[k],r[k+1], output, r[l], r[l+1])""")
         82 function calls in 0.004 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.001    0.001    0.004    0.004 <ipython-input-27-663e142d15fb>:1(no_twodcitera)
        1    0.000    0.000    0.004    0.004 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 _methods.py:43(_count_reduce_items)
        1    0.000    0.000    0.000    0.000 _methods.py:53(_mean)
        1    0.000    0.000    0.000    0.000 function_base.py:436(average)
        1    0.000    0.000    0.000    0.000 index_tricks.py:151(__getitem__)
        3    0.000    0.000    0.002    0.001 numeric.py:1279(roll)
        1    0.000    0.000    0.000    0.000 numeric.py:394(asarray)
        4    0.000    0.000    0.000    0.000 numeric.py:464(asanyarray)
        1    0.000    0.000    0.000    0.000 queues.py:99(put)
        1    0.000    0.000    0.000    0.000 threading.py:299(_is_owned)
        1    0.000    0.000    0.000    0.000 threading.py:372(notify)
        1    0.000    0.000    0.000    0.000 threading.py:63(_note)
        1    0.000    0.000    0.000    0.000 {hasattr}
       18    0.000    0.000    0.000    0.000 {isinstance}
        1    0.000    0.000    0.000    0.000 {issubclass}
        5    0.000    0.000    0.000    0.000 {len}
        3    0.000    0.000    0.000    0.000 {math.ceil}
        1    0.000    0.000    0.000    0.000 {method 'acquire' of '_multiprocessing.SemLock' objects}
        2    0.000    0.000    0.000    0.000 {method 'acquire' of 'thread.lock' objects}
        1    0.000    0.000    0.000    0.000 {method 'append' of 'collections.deque' objects}
        3    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        1    0.000    0.000    0.000    0.000 {method 'mean' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 {method 'reduce' of 'numpy.ufunc' objects}
        1    0.000    0.000    0.000    0.000 {method 'release' of 'thread.lock' objects}
        3    0.002    0.001    0.002    0.001 {method 'take' of 'numpy.ndarray' objects}
        9    0.000    0.000    0.000    0.000 {numpy.core.multiarray.arange}
        5    0.000    0.000    0.000    0.000 {numpy.core.multiarray.array}
        3    0.000    0.000    0.000    0.000 {numpy.core.multiarray.concatenate}
        4    0.000    0.000    0.000    0.000 {range}
        1    0.000    0.000    0.000    0.000 {zip}



In [35]: cProfile.run("""test.twodcitera(dd, res,in1,in2,r[k],r[k+1], output, r[l], r[l+1])""")
         82 function calls in 0.003 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    0.003    0.003 <string>:1(<module>)
        1    0.000    0.000    0.000    0.000 _methods.py:43(_count_reduce_items)
        1    0.000    0.000    0.000    0.000 _methods.py:53(_mean)
        1    0.000    0.000    0.000    0.000 function_base.py:436(average)
        1    0.000    0.000    0.000    0.000 index_tricks.py:151(__getitem__)
        3    0.000    0.000    0.001    0.000 numeric.py:1279(roll)
        1    0.000    0.000    0.000    0.000 numeric.py:394(asarray)
        4    0.000    0.000    0.000    0.000 numeric.py:464(asanyarray)
        1    0.000    0.000    0.000    0.000 queues.py:99(put)
        1    0.000    0.000    0.000    0.000 threading.py:299(_is_owned)
        1    0.000    0.000    0.000    0.000 threading.py:372(notify)
        1    0.000    0.000    0.000    0.000 threading.py:63(_note)
        1    0.000    0.000    0.000    0.000 {hasattr}
       18    0.000    0.000    0.000    0.000 {isinstance}
        1    0.000    0.000    0.000    0.000 {issubclass}
        5    0.000    0.000    0.000    0.000 {len}
        3    0.000    0.000    0.000    0.000 {math.ceil}
        1    0.000    0.000    0.000    0.000 {method 'acquire' of '_multiprocessing.SemLock' objects}
        2    0.000    0.000    0.000    0.000 {method 'acquire' of 'thread.lock' objects}
        1    0.000    0.000    0.000    0.000 {method 'append' of 'collections.deque' objects}
        3    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
        1    0.000    0.000    0.000    0.000 {method 'mean' of 'numpy.ndarray' objects}
        1    0.000    0.000    0.000    0.000 {method 'reduce' of 'numpy.ufunc' objects}
        1    0.000    0.000    0.000    0.000 {method 'release' of 'thread.lock' objects}
        3    0.001    0.000    0.001    0.000 {method 'take' of 'numpy.ndarray' objects}
        9    0.000    0.000    0.000    0.000 {numpy.core.multiarray.arange}
        5    0.000    0.000    0.000    0.000 {numpy.core.multiarray.array}
        3    0.000    0.000    0.000    0.000 {numpy.core.multiarray.concatenate}
        4    0.000    0.000    0.000    0.000 {range}
        1    0.001    0.001    0.003    0.003 {test.twodcitera}
        1    0.000    0.000    0.000    0.000 {zip}

Sadly there is nothing popping up. I would conclude that the reason might be that numpy is already well implemented and most of the time is not lost in the nested loops. Furthermore cPython mostly benefits from static typing. Due to we are using numpy here this might not be a big benefit.

查看更多
0人赞 添加讨论(0) 举报
登录 后发表回答
相关问题
查看全部
相关文章
查看全部
收藏的人(6)