使用 NumPy · Cython 3.0 中文文檔

# 使用 NumPy > 原文： [http://docs.cython.org/en/latest/src/tutorial/numpy.html](http://docs.cython.org/en/latest/src/tutorial/numpy.html) 注意 Cython 0.16 引入了類型化的內存視圖作為此處描述的 NumPy 集成的后續版本。它們比下面的緩沖區語法更容易使用，開銷更少，并且可以在不需要 GIL 的情況下傳遞。它們應該優先于本頁面中提供的語法。有關 NumPy 用戶，請參閱 [Cython。](../userguide/numpy_tutorial.html#numpy-tutorial) 您可以使用 Cython 中的 NumPy 與常規 Python 中的 NumPy 完全相同，但這樣做會導致潛在的高速加速，因為 Cython 支持快速訪問 NumPy 數組。讓我們看一下如何使用一個簡單的例子。下面的代碼使用過濾器對圖像進行 2D 離散卷積（我相信你可以做得更好！，讓它用于演示目的）。它既是有效的 Python 又是有效的 Cython 代碼。我將它稱為 Python 版本的`convolve_py.py`和 Cython 版本的`convolve1.pyx` - Cython 使用“.pyx”作為其文件后綴。 ```py import numpy as np def naive_convolve(f, g): # f is an image and is indexed by (v, w) # g is a filter kernel and is indexed by (s, t), # it needs odd dimensions # h is the output image and is indexed by (x, y), # it is not cropped if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1: raise ValueError("Only odd dimensions on filter supported") # smid and tmid are number of pixels between the center pixel # and the edge, ie for a 5x5 filter they will be 2. # # The output size is calculated by adding smid, tmid to each # side of the dimensions of the input image. vmax = f.shape[0] wmax = f.shape[1] smax = g.shape[0] tmax = g.shape[1] smid = smax // 2 tmid = tmax // 2 xmax = vmax + 2 * smid ymax = wmax + 2 * tmid # Allocate result image. h = np.zeros([xmax, ymax], dtype=f.dtype) # Do convolution for x in range(xmax): for y in range(ymax): # Calculate pixel value for h at (x,y). Sum one component # for each pixel (s, t) of the filter g. s_from = max(smid - x, -smid) s_to = min((xmax - x) - smid, smid + 1) t_from = max(tmid - y, -tmid) t_to = min((ymax - y) - tmid, tmid + 1) value = 0 for s in range(s_from, s_to): for t in range(t_from, t_to): v = x - smid + s w = y - tmid + t value += g[smid - s, tmid - t] * f[v, w] h[x, y] = value return h ``` 這應編譯為生成`yourmod.so`（對于 Linux 系統，在 Windows 系統上，它將是`yourmod.pyd`）。我們運行 Python 會話來測試 Python 版本（從`.py` -file 導入）和編譯的 Cython 模塊。 ```py In [1]: import numpy as np In [2]: import convolve_py In [3]: convolve_py.naive_convolve(np.array([[1, 1, 1]], dtype=np.int), ... np.array([[1],[2],[1]], dtype=np.int)) Out [3]: array([[1, 1, 1], [2, 2, 2], [1, 1, 1]]) In [4]: import convolve1 In [4]: convolve1.naive_convolve(np.array([[1, 1, 1]], dtype=np.int), ... np.array([[1],[2],[1]], dtype=np.int)) Out [4]: array([[1, 1, 1], [2, 2, 2], [1, 1, 1]]) In [11]: N = 100 In [12]: f = np.arange(N*N, dtype=np.int).reshape((N,N)) In [13]: g = np.arange(81, dtype=np.int).reshape((9, 9)) In [19]: %timeit -n2 -r3 convolve_py.naive_convolve(f, g) 2 loops, best of 3: 1.86 s per loop In [20]: %timeit -n2 -r3 convolve1.naive_convolve(f, g) 2 loops, best of 3: 1.41 s per loop ``` 還沒有那么大的差別;因為 C 代碼仍然完全符合 Python 解釋器的作用（例如，意味著為每個使用的數字分配了一個新對象）。查看生成的 html 文件，看看即使是最簡單的語句，您需要快速得到什么。我們需要給 Cython 更多信息;我們需要添加類型。 ## 添加類型要添加類型，我們使用自定義 Cython 語法，因此我們現在正在破壞 Python 源兼容性。考慮一下這段代碼（_ 閱讀評論！_）： ```py # tag: numpy_old # You can ignore the previous line. # It's for internal testing of the cython documentation. import numpy as np # "cimport" is used to import special compile-time information # about the numpy module (this is stored in a file numpy.pxd which is # currently part of the Cython distribution). cimport numpy as np # We now need to fix a datatype for our arrays. I've used the variable # DTYPE for this, which is assigned to the usual NumPy runtime # type info object. DTYPE = np.int # "ctypedef" assigns a corresponding compile-time type to DTYPE_t. For # every type in the numpy module there's a corresponding compile-time # type with a _t-suffix. ctypedef np.int_t DTYPE_t # "def" can type its arguments but not have a return type. The type of the # arguments for a "def" function is checked at run-time when entering the # function. # # The arrays f, g and h is typed as "np.ndarray" instances. The only effect # this has is to a) insert checks that the function arguments really are # NumPy arrays, and b) make some attribute access like f.shape[0] much # more efficient. (In this example this doesn't matter though.) def naive_convolve(np.ndarray f, np.ndarray g): if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1: raise ValueError("Only odd dimensions on filter supported") assert f.dtype == DTYPE and g.dtype == DTYPE # The "cdef" keyword is also used within functions to type variables. It # can only be used at the top indentation level (there are non-trivial # problems with allowing them in other places, though we'd love to see # good and thought out proposals for it). # # For the indices, the "int" type is used. This corresponds to a C int, # other C types (like "unsigned int") could have been used instead. # Purists could use "Py_ssize_t" which is the proper Python type for # array indices. cdef int vmax = f.shape[0] cdef int wmax = f.shape[1] cdef int smax = g.shape[0] cdef int tmax = g.shape[1] cdef int smid = smax // 2 cdef int tmid = tmax // 2 cdef int xmax = vmax + 2 * smid cdef int ymax = wmax + 2 * tmid cdef np.ndarray h = np.zeros([xmax, ymax], dtype=DTYPE) cdef int x, y, s, t, v, w # It is very important to type ALL your variables. You do not get any # warnings if not, only much slower code (they are implicitly typed as # Python objects). cdef int s_from, s_to, t_from, t_to # For the value variable, we want to use the same data type as is # stored in the array, so we use "DTYPE_t" as defined above. # NB! An important side-effect of this is that if "value" overflows its # datatype size, it will simply wrap around like in C, rather than raise # an error like in Python. cdef DTYPE_t value for x in range(xmax): for y in range(ymax): s_from = max(smid - x, -smid) s_to = min((xmax - x) - smid, smid + 1) t_from = max(tmid - y, -tmid) t_to = min((ymax - y) - tmid, tmid + 1) value = 0 for s in range(s_from, s_to): for t in range(t_from, t_to): v = x - smid + s w = y - tmid + t value += g[smid - s, tmid - t] * f[v, w] h[x, y] = value return h ``` 在構建完這個并繼續我的（非正式）基準測試后，我得到： ```py In [21]: import convolve2 In [22]: %timeit -n2 -r3 convolve2.naive_convolve(f, g) 2 loops, best of 3: 828 ms per loop ``` ## 高效索引仍有瓶頸殺戮性能，那就是數組查找和分配。 `[]` -operator 仍然使用完整的 Python 操作 - 我們想要做的是直接以 C 速度訪問數據緩沖區。我們需要做的是輸入`ndarray`對象的內容。我們使用特殊的“緩沖區”語法來執行此操作，必須告訴數據類型（第一個參數）和維度數量（“ndim”僅限關鍵字參數，如果未提供，則假定為一維）。這些是必要的變化： ```py ... def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g): ... cdef np.ndarray[DTYPE_t, ndim=2] h = ... ``` 用法： ```py In [18]: import convolve3 In [19]: %timeit -n3 -r100 convolve3.naive_convolve(f, g) 3 loops, best of 100: 11.6 ms per loop ``` 請注意這種變化的重要性。 _Gotcha_ ：這種有效的索引僅影響某些索引操作，即那些具有完全`ndim`數量的類型化整數索引的索引操作。因此，如果沒有輸入`v`，則查找`f[v, w]`不會被優化。另一方面，這意味著您可以繼續使用 Python 對象進行復雜的動態切片等，就像未鍵入數組時一樣。 ## 進一步調整索引數組查找仍然受到兩個因素的影響： 1. 進行邊界檢查。 2. 檢查負指數并正確處理。上面的代碼是明確編碼的，因此它不使用負索引，并且（希望）總是在邊界內訪問。我們可以添加一個裝飾器來禁用邊界檢查： ```py ... cimport cython @cython.boundscheck(False) # turn off bounds-checking for entire function @cython.wraparound(False) # turn off negative index wrapping for entire function def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g): ... ``` 現在沒有執行邊界檢查（并且，作為一個副作用，如果你'碰巧'訪問越界，你將在最好的情況下崩潰你的程序，在最壞的情況下會損壞數據）。可以通過多種方式切換邊界檢查模式，有關詳細信息，請參閱 [編譯器指令](../userguide/source_files_and_compilation.html#compiler-directives) 。此外，我們已禁用檢查以包裝負指數（例如 g [-1]給出最后一個值）。與禁用邊界檢查一樣，如果我們嘗試實際使用帶有此禁用的負索引，則會發生錯誤。函數調用開銷現在開始發揮作用，因此我們將后兩個示例與較大的 N 進行比較： ```py In [11]: %timeit -n3 -r100 convolve4.naive_convolve(f, g) 3 loops, best of 100: 5.97 ms per loop In [12]: N = 1000 In [13]: f = np.arange(N*N, dtype=np.int).reshape((N,N)) In [14]: g = np.arange(81, dtype=np.int).reshape((9, 9)) In [17]: %timeit -n1 -r10 convolve3.naive_convolve(f, g) 1 loops, best of 10: 1.16 s per loop In [18]: %timeit -n1 -r10 convolve4.naive_convolve(f, g) 1 loops, best of 10: 597 ms per loop ``` （這也是一個混合基準，因為結果數組是在函數調用中分配的。）警告速度需要一些成本。特別是將類型對象（如我們的示例代碼中的`f`，`g`和`h`）設置為`None`會很危險。將這些對象設置為`None`是完全合法的，但您可以使用它們檢查它們是否為 None。所有其他用途（屬性查找或索引）都可能會破壞或損壞數據（而不是像在 Python 中那樣引發異常）。實際規則有點復雜但主要信息很明確：不要使用類型化對象而不知道它們沒有設置為 None。 ## 更通用的代碼有可能做到： ```py def naive_convolve(object[DTYPE_t, ndim=2] f, ...): ``` 即使用 [`object`](https://docs.python.org/3/library/functions.html#object "(in Python v3.7)") 而不是`np.ndarray`。在 Python 3.0 下，這可以允許您的算法使用任何支持緩沖區接口的庫;并支持例如如果有人對 Python 2.x 感興趣，可以輕松添加 Python Imaging Library。但是這有一些速度損失（如果類型設置為`np.ndarray`，則編譯時假設編譯時間更多，特別是假設數據以純步進模式而不是間接模式存儲）。