使用 NumPy

原文： http://docs.cython.org/en/latest/src/tutorial/numpy.html

注意

Cython 0.16 引入了类型化的内存视图作为此处描述的 NumPy 集成的后续版本。它们比下面的缓冲区语法更容易使用，开销更少，并且可以在不需要 GIL 的情况下传递。它们应该优先于本页面中提供的语法。有关 NumPy 用户 ，请参阅 Cython。

您可以使用 Cython 中的 NumPy 与常规 Python 中的 NumPy 完全相同，但这样做会导致潜在的高速加速，因为 Cython 支持快速访问 NumPy 数组。让我们看一下如何使用一个简单的例子。

下面的代码使用过滤器对图像进行 2D 离散卷积（我相信你可以做得更好！，让它用于演示目的）。它既是有效的 Python 又是有效的 Cython 代码。我将它称为 Python 版本的convolve_py.py和 Cython 版本的convolve1.pyx - Cython 使用“.pyx”作为其文件后缀。

  1. import numpy as np
  3. def naive_convolve(f, g):
  4.     # f is an image and is indexed by (v, w)
  5.     # g is a filter kernel and is indexed by (s, t),
  6.     #   it needs odd dimensions
  7.     # h is the output image and is indexed by (x, y),
  8.     #   it is not cropped
  9.     if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
  10.         raise ValueError("Only odd dimensions on filter supported")
  11.     # smid and tmid are number of pixels between the center pixel
  12.     # and the edge, ie for a 5x5 filter they will be 2.
  13.     #
  14.     # The output size is calculated by adding smid, tmid to each
  15.     # side of the dimensions of the input image.
  16.     vmax = f.shape[0]
  17.     wmax = f.shape[1]
  18.     smax = g.shape[0]
  19.     tmax = g.shape[1]
  20.     smid = smax // 2
  21.     tmid = tmax // 2
  22.     xmax = vmax + 2 * smid
  23.     ymax = wmax + 2 * tmid
  24.     # Allocate result image.
  25.     h = np.zeros([xmax, ymax], dtype=f.dtype)
  26.     # Do convolution
  27.     for x in range(xmax):
  28.         for y in range(ymax):
  29.             # Calculate pixel value for h at (x,y). Sum one component
  30.             # for each pixel (s, t) of the filter g.
  31.             s_from = max(smid - x, -smid)
  32.             s_to = min((xmax - x) - smid, smid + 1)
  33.             t_from = max(tmid - y, -tmid)
  34.             t_to = min((ymax - y) - tmid, tmid + 1)
  35.             value = 0
  36.             for s in range(s_from, s_to):
  37.                 for t in range(t_from, t_to):
  38.                     v = x - smid + s
  39.                     w = y - tmid + t
  40.                     value += g[smid - s, tmid - t] * f[v, w]
  41.             h[x, y] = value
  42.     return h

这应编译为生成yourmod.so（对于 Linux 系统，在 Windows 系统上，它将是yourmod.pyd）。我们运行 Python 会话来测试 Python 版本（从.py -file 导入）和编译的 Cython 模块。

  1. In [1]: import numpy as np
  2. In [2]: import convolve_py
  3. In [3]: convolve_py.naive_convolve(np.array([[1, 1, 1]], dtype=np.int),
  4. ...     np.array([[1],[2],[1]], dtype=np.int))
  5. Out [3]:
  6. array([[1, 1, 1],
  7.  [2, 2, 2],
  8.  [1, 1, 1]])
  9. In [4]: import convolve1
  10. In [4]: convolve1.naive_convolve(np.array([[1, 1, 1]], dtype=np.int),
  11. ...     np.array([[1],[2],[1]], dtype=np.int))
  12. Out [4]:
  13. array([[1, 1, 1],
  14.  [2, 2, 2],
  15.  [1, 1, 1]])
  16. In [11]: N = 100
  17. In [12]: f = np.arange(N*N, dtype=np.int).reshape((N,N))
  18. In [13]: g = np.arange(81, dtype=np.int).reshape((9, 9))
  19. In [19]: %timeit -n2 -r3 convolve_py.naive_convolve(f, g)
  20. 2 loops, best of 3: 1.86 s per loop
  21. In [20]: %timeit -n2 -r3 convolve1.naive_convolve(f, g)
  22. 2 loops, best of 3: 1.41 s per loop

还没有那么大的差别;因为 C 代码仍然完全符合 Python 解释器的作用（例如，意味着为每个使用的数字分配了一个新对象）。查看生成的 html 文件，看看即使是最简单的语句，您需要快速得到什么。我们需要给 Cython 更多信息;我们需要添加类型。

添加类型

要添加类型，我们使用自定义 Cython 语法，因此我们现在正在破坏 Python 源兼容性。考虑一下这段代码（ 阅读评论！）：

  1. # tag: numpy_old
  2. # You can ignore the previous line.
  3. # It's for internal testing of the cython documentation.
  5. import numpy as np
  7. # "cimport" is used to import special compile-time information
  8. # about the numpy module (this is stored in a file numpy.pxd which is
  9. # currently part of the Cython distribution).
  10. cimport numpy as np
  12. # We now need to fix a datatype for our arrays. I've used the variable
  13. # DTYPE for this, which is assigned to the usual NumPy runtime
  14. # type info object.
  15. DTYPE = np.int
  17. # "ctypedef" assigns a corresponding compile-time type to DTYPE_t. For
  18. # every type in the numpy module there's a corresponding compile-time
  19. # type with a _t-suffix.
  20. ctypedef np.int_t DTYPE_t
  22. # "def" can type its arguments but not have a return type. The type of the
  23. # arguments for a "def" function is checked at run-time when entering the
  24. # function.
  25. #
  26. # The arrays f, g and h is typed as "np.ndarray" instances. The only effect
  27. # this has is to a) insert checks that the function arguments really are
  28. # NumPy arrays, and b) make some attribute access like f.shape[0] much
  29. # more efficient. (In this example this doesn't matter though.)
  30. def naive_convolve(np.ndarray f, np.ndarray g):
  31.     if g.shape[0] % 2 != 1 or g.shape[1] % 2 != 1:
  32.         raise ValueError("Only odd dimensions on filter supported")
  33.     assert f.dtype == DTYPE and g.dtype == DTYPE
  35.     # The "cdef" keyword is also used within functions to type variables. It
  36.     # can only be used at the top indentation level (there are non-trivial
  37.     # problems with allowing them in other places, though we'd love to see
  38.     # good and thought out proposals for it).
  39.     #
  40.     # For the indices, the "int" type is used. This corresponds to a C int,
  41.     # other C types (like "unsigned int") could have been used instead.
  42.     # Purists could use "Py_ssize_t" which is the proper Python type for
  43.     # array indices.
  44.     cdef int vmax = f.shape[0]
  45.     cdef int wmax = f.shape[1]
  46.     cdef int smax = g.shape[0]
  47.     cdef int tmax = g.shape[1]
  48.     cdef int smid = smax // 2
  49.     cdef int tmid = tmax // 2
  50.     cdef int xmax = vmax + 2 * smid
  51.     cdef int ymax = wmax + 2 * tmid
  52.     cdef np.ndarray h = np.zeros([xmax, ymax], dtype=DTYPE)
  53.     cdef int x, y, s, t, v, w
  55.     # It is very important to type ALL your variables. You do not get any
  56.     # warnings if not, only much slower code (they are implicitly typed as
  57.     # Python objects).
  58.     cdef int s_from, s_to, t_from, t_to
  60.     # For the value variable, we want to use the same data type as is
  61.     # stored in the array, so we use "DTYPE_t" as defined above.
  62.     # NB! An important side-effect of this is that if "value" overflows its
  63.     # datatype size, it will simply wrap around like in C, rather than raise
  64.     # an error like in Python.
  65.     cdef DTYPE_t value
  66.     for x in range(xmax):
  67.         for y in range(ymax):
  68.             s_from = max(smid - x, -smid)
  69.             s_to = min((xmax - x) - smid, smid + 1)
  70.             t_from = max(tmid - y, -tmid)
  71.             t_to = min((ymax - y) - tmid, tmid + 1)
  72.             value = 0
  73.             for s in range(s_from, s_to):
  74.                 for t in range(t_from, t_to):
  75.                     v = x - smid + s
  76.                     w = y - tmid + t
  77.                     value += g[smid - s, tmid - t] * f[v, w]
  78.             h[x, y] = value
  79.     return h

在构建完这个并继续我的（非正式）基准测试后，我得到：

  1. In [21]: import convolve2
  2. In [22]: %timeit -n2 -r3 convolve2.naive_convolve(f, g)
  3. 2 loops, best of 3: 828 ms per loop

高效索引

仍有瓶颈杀戮性能，那就是数组查找和分配。 [] -operator 仍然使用完整的 Python 操作 - 我们想要做的是直接以 C 速度访问数据缓冲区。

我们需要做的是输入ndarray对象的内容。我们使用特殊的“缓冲区”语法来执行此操作，必须告诉数据类型（第一个参数）和维度数量（“ndim”仅限关键字参数，如果未提供，则假定为一维）。

这些是必要的变化：

  1. ...
  2. def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g):
  3. ...
  4. cdef np.ndarray[DTYPE_t, ndim=2] h = ...

用法：

  1. In [18]: import convolve3
  2. In [19]: %timeit -n3 -r100 convolve3.naive_convolve(f, g)
  3. 3 loops, best of 100: 11.6 ms per loop

请注意这种变化的重要性。

Gotcha ：这种有效的索引仅影响某些索引操作，即那些具有完全ndim数量的类型化整数索引的索引操作。因此，如果没有输入v，则查找f[v, w]不会被优化。另一方面，这意味着您可以继续使用 Python 对象进行复杂的动态切片等，就像未键入数组时一样。

进一步调整索引

数组查找仍然受到两个因素的影响：

  1. 进行边界检查。

  2. 检查负指数并正确处理。上面的代码是明确编码的，因此它不使用负索引，并且（希望）总是在边界内访问。我们可以添加一个装饰器来禁用边界检查：

    1. ...
    2. cimport cython
    3. @cython.boundscheck(False) # turn off bounds-checking for entire function
    4. @cython.wraparound(False)  # turn off negative index wrapping for entire function
    5. def naive_convolve(np.ndarray[DTYPE_t, ndim=2] f, np.ndarray[DTYPE_t, ndim=2] g):
    6. ...

现在没有执行边界检查（并且，作为一个副作用，如果你’碰巧’访问越界，你将在最好的情况下崩溃你的程序，在最坏的情况下会损坏数据）。可以通过多种方式切换边界检查模式，有关详细信息，请参阅 编译器指令 。

此外，我们已禁用检查以包装负指数（例如 g [-1]给出最后一个值）。与禁用边界检查一样，如果我们尝试实际使用带有此禁用的负索引，则会发生错误。

函数调用开销现在开始发挥作用，因此我们将后两个示例与较大的 N 进行比较：

  1. In [11]: %timeit -n3 -r100 convolve4.naive_convolve(f, g)
  2. 3 loops, best of 100: 5.97 ms per loop
  3. In [12]: N = 1000
  4. In [13]: f = np.arange(N*N, dtype=np.int).reshape((N,N))
  5. In [14]: g = np.arange(81, dtype=np.int).reshape((9, 9))
  6. In [17]: %timeit -n1 -r10 convolve3.naive_convolve(f, g)
  7. 1 loops, best of 10: 1.16 s per loop
  8. In [18]: %timeit -n1 -r10 convolve4.naive_convolve(f, g)
  9. 1 loops, best of 10: 597 ms per loop

（这也是一个混合基准，因为结果数组是在函数调用中分配的。）

警告

速度需要一些成本。特别是将类型对象（如我们的示例代码中的f，g和h）设置为None会很危险。将这些对象设置为None是完全合法的，但您可以使用它们检查它们是否为 None。所有其他用途（属性查找或索引）都可能会破坏或损坏数据（而不是像在 Python 中那样引发异常）。

实际规则有点复杂但主要信息很明确：不要使用类型化对象而不知道它们没有设置为 None。

更通用的代码

有可能做到：

  1. def naive_convolve(object[DTYPE_t, ndim=2] f, ...):

即使用 object 而不是np.ndarray。在 Python 3.0 下，这可以允许您的算法使用任何支持缓冲区接口的库;并支持例如如果有人对 Python 2.x 感兴趣，可以轻松添加 Python Imaging Library。

但是这有一些速度损失（如果类型设置为np.ndarray，则编译时假设编译时间更多，特别是假设数据以纯步进模式而不是间接模式存储）。