地用Cython創建小數組需要

Question

我正在寫一個新的隨機數發生器，用於numpy的根據任意分佈，當我遇到這確實怪異的行爲來產生隨機數的時候一個巨大的量：

這是考驗。 PYX

setup.py

from distutils.core import setup 
from Cython.Build import cythonize 
setup(name = "simple cython func",ext_modules = cythonize('test.pyx'),) 

配置代碼

#!/usr/bin/python 
from __future__ import division 

import subprocess 
import timeit 

#Compile the cython modules before importing them 
subprocess.call(['python', 'setup.py', 'build_ext', '--inplace']) 

sstr=""" 
import test 
import numpy 
u=numpy.random.random(10) 
a=numpy.random.random(10) 
a=numpy.cumsum(a) 
a/=a[-1] 
r=numpy.empty(10,int) 
""" 

print "binary search: creates an array[N] and performs N binary searches to fill it:\n",timeit.timeit('numpy.searchsorted(a,u)',sstr) 
print "Simple replacement for binary search:takes the same args as np.searchsorted and similarly returns a new array. this performs only one trivial operation per element:\n",timeit.timeit('test.BSReplacement(a,u)',sstr) 

print "barebones function doing nothing:",timeit.timeit('test.BareBones(a,u,r)',sstr) 
print "Untyped inputs and doing N iterations:",timeit.timeit('test.UntypedWithLoop(a,u,r)',sstr) 
print "time for just np.empty()",timeit.timeit('numpy.empty(10,int)',sstr) 

二分查找執行需要執行的時間順序爲len(u)*Log(len(a))。簡單的cython函數按照len(u)的順序運行。兩者都返回一個len（u）的1D int數組。然而，即使這樣，在numpy庫中，沒有計算繁瑣的實現需要比全二分搜索更長的時間。 （這是寫在C：https://github.com/numpy/numpy/blob/202e78d607515e0390cffb1898e11807f117b36a/numpy/core/src/multiarray/item_selection.c看到PyArray_SearchSorted）

的結果是：

binary search: creates an array[N] and performs N binary searches to fill it: 
1.15157485008 
Simple replacement for binary search:takes the same args as np.searchsorted and similarly returns a new array. this performs only one trivial operation per element: 
3.69442796707 
barebones function doing nothing: 0.87496304512 
Untyped inputs and doing N iterations: 0.244267940521 
time for just np.empty() 1.0983929634 

爲什麼np.empty（）步走這麼多的時間？我能做些什麼來獲得我可以返回的空數組？

C函數執行此操作並運行一大堆理智檢查，並在內部循環中使用較長的算法。 （ⅰ去除除循環本身來回我的示例中，所有的邏輯）

---

更新

原來有兩個不同的問題：

1. 的np.empty（10）單獨的呼叫具有巨大的開銷並且花費盡可能多的時間來進行搜索排序以創建新的數組並且執行10個二進制搜索
2. 只聲明緩衝區語法np.ndarray[...]也有比接收無類型變量和迭代50次花費更多時間的大量開銷。

結果50次迭代：

binary search: 2.45336699486 
Simple replacement:3.71126317978 
barebones function doing nothing: 0.924916028976 
Untyped inputs and doing N iterations: 0.316384077072 
time for just np.empty() 1.04949498177 

Answer 1

創建np.empty在用Cython函數內部有一定的開銷，因爲你已經看到了。在這裏，您將看到有關如何創建空數組並把它傳遞給用Cython模塊，以填補與正確值的示例：

n=10：

numpy.searchsorted: 1.30574745517 
cython O(1): 3.28732016088 
cython no array declaration 1.54710909596 

n=100：

numpy.searchsorted: 4.15200545373 
cython O(1): 13.7273431067 
cython no array declaration 11.4186086744 

正如您已經指出的那樣，numpy版本的規模更好，因爲它是O(len(u)*long(len(a)))，此處的算法是O(len(u)*len(a)) ...

我也嘗試過使用Memoryview，基本上改變np.ndarray[double, ndim=1]double[:]，但在這種情況下第一個選項更快。

新.pyx文件是：

from __future__ import division 
import numpy as np 
cimport numpy as np 
cimport cython 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def JustLoop(np.ndarray[double, ndim=1] a, np.ndarray[double, ndim=1] u, 
      np.ndarray[int, ndim=1] r): 
    cdef int i,j 
    for j in range(u.shape[0]): 
     if u[j] < a[0]: 
      r[j] = 0 
      continue 

     if u[j] > a[a.shape[0]-1]: 
      r[j] = a.shape[0]-1 
      continue 

     for i in range(1, a.shape[0]): 
      if u[j] >= a[i-1] and u[j] < a[i]: 
       r[j] = i 
       break 

@cython.boundscheck(False) 
@cython.wraparound(False) 
def WithArray(np.ndarray[double, ndim=1] a, np.ndarray[double, ndim=1] u): 
    cdef np.ndarray[np.int_t, ndim=1] r=np.empty(u.shape[0],dtype=int) 
    cdef int i,j 
    for j in range(u.shape[0]): 
     if u[j] < a[0]: 
      r[j] = 0 
      continue 

     if u[j] > a[a.shape[0]-1]: 
      r[j] = a.shape[0]-1 
      continue 

     for i in range(1, a.shape[0]): 
      if u[j] >= a[i-1] and u[j] < a[i]: 
       r[j] = i 
       break 
    return r 

新.py文件：

import numpy 
import subprocess 
import timeit 

#Compile the cython modules before importing them 
subprocess.call(['python', 'setup.py', 'build_ext', '--inplace']) 
from test import * 

sstr=""" 
import test 
import numpy 
u=numpy.random.random(10) 
a=numpy.random.random(10) 
a=numpy.cumsum(a) 
a/=a[-1] 
a.sort() 
r = numpy.empty(u.shape[0], dtype=int) 
""" 

print "numpy.searchsorted:",timeit.timeit('numpy.searchsorted(a,u)',sstr) 
print "cython O(1):",timeit.timeit('test.WithArray(a,u)',sstr) 
print "cython no array declaration",timeit.timeit('test.JustLoop(a,u,r)',sstr) 

Answer 2

還有就是這個用Cython名單上可能有一些有益的建議進行討論： https://groups.google.com/forum/#!topic/cython-users/CwtU_jYADgM

雖然我嘗試在Cython外部分配小數組，但將它們傳遞給d在隨後的方法調用中重新使用它們。我明白，這並不總是一種選擇。