2016-02-11 74 views
0

我正在嘗試爲cython和numpy獲取FLOPS基準。爲此我寫了一個cython程序。這裏是:Cython性能基準

cimport numpy as np 
import numpy as np 
import time 

cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def numpybenchmark(): 

    cdef np.ndarray[np.float64_t, ndim=2] m1 = np.random.rand(3,3) 
    cdef np.ndarray[np.float64_t, ndim=1] m2 = np.random.rand(3) 
    cdef np.ndarray[np.float64_t, ndim=1] res 

    cdef int niters = 10000000 
    cdef int x 

    t1 = time.time() 
    for x in range(niters): 
     res = np.dot(m1, m2) 
    t2 = time.time() 
    cdef double numopsperloop = 9. + 6. 
    cdef double totalops = numopsperloop * float(niters) 
    cdef double mflops = totalops/(t2-t1)/1024./1024. 
    print 'Computed MFLops is: ' + str(mflops) 

在我的機器上我測量「計算MFLops是:7.42390102416」。我的機器具有Intel Core i7-6700HQ CPU @ 2.6 GHz並運行Windows 10.

如果要在機器上運行它,請將代碼保存在名爲「benchmark.pyx」的文件中。然後創建一個名爲「setup.py」與文件,內容如下:

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

setup(                       
    ext_modules = cythonize("benchmark.pyx"), 
    include_dirs=[numpy.get_include()]               
) 

那麼你應該能夠「蟒蛇setup.py build_ext --inplace」進行編譯。在Windows中,遇到可怕的「無法找到vcvarsall.bat」錯誤並且不得不花費大量精力來解決這個問題時,它可能會更加困難。

這個表現對我來說似乎很差。我想知道是否有人可以在他們的平臺上運行它並告訴我你得到了什麼?或者指出我的代碼中對性能有不利影響的任何明顯錯誤?

謝謝!

+1

這似乎是它會更適合codereview.stackexchange.com –

回答

1

Cython實際上並沒有消除np.dot上的任何Python調用開銷。這包括(請注意,名單並不詳盡,並可能地方稍有不妥,但它給人的要點):

  1. 尋找np.dot打電話:

    • 在字典查找全局名字空間爲np
    • np的字典查詢爲dot的名字空間。 (請注意,以上所有內容都可以通過在功能中執行dot = np.dot來消除,然後調用dot
    • 關於dot__call__的字典查詢。 (這可通過更快的機制來完成,如果點是一個C/Fortran的編譯函數)
  2. 包裝製備的論據np.dot

    • 創建包含傳遞給np.dot
    • 兩個參數的元組
    • 增加每個參數的引用計數。
  3. np.dot接着處理所述參數...

    • 解包元組的每個參數
    • 檢查在numpy的陣列。
    • 檢查每個numpy數組的dtype是否相同,並基於dtype選擇哪個BLAS函數來調用。
    • 檢查數組尺寸並確保它們匹配。
  4. ...的輸出參數分配空間...

    • 分配新對象np.ndarray
    • 遞增該
    • 分配空間的引用計數爲ndarray內的物理陣列
  5. ...調用BLAS操作,給出您在浮點運算...

  6. ...而遞減的它傳遞的輸入參數引用計數(檢查,看是否有應被釋放,但沒有人會是)

  7. 您的通話功能然後必須:

    • 檢查是否異常,通過np.dot
    • 凸起接收輸出陣列(可能還有一些的refcount這裏雜耍)
    • 遞減引用計數以前的內容res
    • 釋放以前的內容res,記住它至少需要2步處理,因爲陣列與ndarray持有者分開持有。

如果你想使這個最(可能除了分配)微不足道相比,矩陣向量乘法,那麼你需要做顯著更大的陣列的測量。您可以使用np.dot中的out可選參數擺脫分配。如果您想使其全部消失,則可以使用scipy Cython BLAS interface直接調用BLAS功能。

+0

謝謝大衛!很高興知道。 – jungleman007

1

在仔細閱讀DavidW的文章並做了一些實驗之後,我找到了避免所有numpy開銷的方法。它涉及使用指針,特別是不會將numpy數組傳遞給循環內的函數。

這裏是全碼:

cimport numpy as np 
import numpy as np 
import time 


cdef matrixdotvector(double* mat, int numrows, int numcols, double* vec, double* outputvec): 
    outputvec[0] = mat[0+0*numcols] * vec[0] + mat[1+0*numcols] * vec[1] + mat[2+0*numcols] * vec[2] 
    outputvec[1] = mat[0+1*numcols] * vec[0] + mat[1+1*numcols] * vec[1] + mat[2+1*numcols] * vec[2] 
    outputvec[2] = mat[0+2*numcols] * vec[0] + mat[1+2*numcols] * vec[1] + mat[2+2*numcols] * vec[2] 

cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def numpybenchmark(): 

    cdef np.ndarray[np.float64_t, ndim=2] m1 = np.random.rand(3,3) 
    cdef np.ndarray[np.float64_t, ndim=1] m2 = np.transpose(np.random.rand(3)) 
    cdef np.ndarray[np.float64_t, ndim=1] res 

    cdef int niters = 10000000 
    cdef int x 

    t1 = time.time() 
    for x in range(niters): 
     res = np.dot(m1, m2) 
    t2 = time.time() 
    cdef double numopsperloop = 9. + 6. 
    cdef double totalops = numopsperloop * float(niters) 
    cdef double mflops = totalops/(t2-t1)/1024./1024. 
    print 'Computed MFLops is: ' + str(mflops) 

cimport cython 
@cython.boundscheck(False) 
@cython.wraparound(False) 
def numpybenchmark2(): 

    cdef int numrows = 3 
    cdef int numcols = 3 
    cdef np.ndarray[np.float64_t, ndim=2] m1 = np.random.rand(3,3) 
    cdef np.ndarray[np.float64_t, ndim=1] m2 = np.transpose(np.random.rand(3)) 
    cdef np.ndarray[np.float64_t, ndim=1] res = np.zeros(3) 

    cdef int niters = 10000000 
    cdef int x 

    t1 = time.time() 
    for x in range(niters): 
     matrixdotvector(&m1[0,0], numrows, numcols, &m2[0], &res[0]) 
    t2 = time.time() 

    assert (np.linalg.norm(np.dot(m1,m2) - res) < 1.0e-6), "Arrays do not match" 

    cdef double numopsperloop = 9. + 6. 
    cdef double totalops = numopsperloop * float(niters) 
    cdef double mflops = totalops/(t2-t1)/1024./1024. 
    print 'Computed MFLops is: ' + str(mflops) 

numpybenchmark之間的大的差異()和numpybenchmark2()是,我通過使指針指向numpy的數據陣列中numpybenchmark2避免所有numpy的開銷()(而不是傳遞鍵入的numpy對象,這也很慢)。我通過展開並在代碼中重新實現它來避免np.dot計算的開銷。

所以,現在我得到的基準測試結果是:

在[13]:benchmark.numpybenchmark() 已計算MFLOPS是:7.3412268815

在[14]:benchmark.numpybenchmark2() 計算MFLops是:1521.81908107

所以這是一個相當大的增加。老實說,這不是一種「pythonic」的方式,但它是快速尖叫,所以在某些情況下可能有用。由於matrixdotvector()中的代碼看起來像C代碼,所以我們可以說這應該全部用C編寫。就我個人而言,我發現使用類C的Cython代碼而不是直接跳入C的原型會更快。

無論如何,對於正在學習cython的人來說,這篇文章可能有用。