Skip to content

用于 NumPy 用户的 Cython

原文: http://docs.cython.org/en/latest/src/userguide/numpy_tutorial.html

本教程针对的是根本没有 Cython 经验的 NumPy 用户。如果您对 Cython 有一些了解,可以直接跳到“高效索引”部分。

这里考虑的主要场景是 NumPy 最终用途而不是 NumPy / SciPy 开发。原因是 Cython 不能(还)能够以高级方式支持关于维数的通用功能。对于 SciPy 开发而言,这种限制比更具体的“最终用户”功能要严格得多。有关详细信息,请参阅最后一节。

本教程的风格不适合所有人,所以你也可以考虑:

Cython 初识

Cython 是一个编译器,它将 Python风格 的代码文件编译为 C 代码。但是,Cython 不是 Python 到 C 的翻译机。也就是说,它不会接收您的完整程序并“将其转换为 C语言代码” - 相反,编译结果会充分利用 Python 运行时环境。一种看待它的方式可能是你的代码仍然是 Python,因为它在 Python 运行时环境中运行,而不是像你编译为机器码一样的编译了的解释过的Python字节码(但是增加了额外的语法以便于嵌入更快的类 C 代码)。

这有两个重要的结果:

  • 速度。提升多少在很大程度上取决于程序内容。典型的 Python 数值程序往往会获得很少的收益,因为大部分时间花在以高级方式使用的较低级别 C 上。然而,当添加类型信息时,for 循环风格的程序可以获得巨大的提升(并且因此可以作为现实的替代方案)。
  • 轻松调用 C 代码。 Cython 的目的之一是允许轻松包装 C 库。在 Cython 中编写代码时,您可以像使用 Python 代码一样轻松调用 C 代码。

只有很少有 Python 构造还没在Cython被支持,但是Cython的一个明确的目标是能够编译所有 Python 代码,你可以在 限制 这一节文档中看到与 Python 的差异。

您的 Cython 环境

使用 Cython 包含以下步骤:

  1. 写一个 .pyx 源文件
  2. 运行 Cython 编译器以生成 C 文件
  3. 运行 C 编译器以生成编译库
  4. 运行 Python 解释器并要求它导入模块

但是,有几个选项可以自动执行这些步骤:

  1. SAGE 数学软件系统为从交互式命令行或笔记本界面(如 Maple / Mathematica)使用 Cython 和 NumPy 提供了出色的支持。参见本文档
  2. Cython 可以作为 Jupyter 笔记本中的扩展,使得编译和使用 Cython 代码只需在单元格顶部只有一个%%cython。有关更多信息,请参阅 使用 Jupyter 笔记本
  3. pyximport 的一个版本随 Cython 一起提供,因此您可以动态地将 pyx 文件导入 Python 并自动编译(参见 使用 pyximport进行编译)。
  4. Cython 支持 setuptools,因此您可以非常轻松地创建自动化过程的构建脚本,这是 Cython 实现的库和包的首选方法。参见 Basic setup.py
  5. 手动编译(见下文)

注意: 如果使用除 SAGE 之外的其他交互式命令行环境(如 IPython 或 Python 本身),则在重新编译模块时重新启动该进程非常重要。再次发出“导入”声明是不够的。

安装

如果您已经有一个 C 编译器,那么就这样做:

pip install Cython

否则,请参见 安装页面

在撰写本文时,SAGE 附带了比本教程所要求的更早版本的 Cython。因此,如果使用 SAGE,您应该下载最新的 Cython 然后执行

$ cd path/to/cython-distro
$ path-to-sage/sage -python setup.py install

这将把最新的 Cython 安装到 SAGE 中。

手动编译

由于知道发生了什么总是很重要,我将在这里描述手动方法。第一次运行 Cython:

$ cython yourmod.pyx

这会创建yourmod.c,这是 Python 扩展模块的 C 源代码。一个有用的附加开关是-a,它将生成一个文档yourmod.html,它显示Cython 代码和逐行转换 C 代码的。

然后我们编译 C 文件。这可能会因您的系统而异,但 C 文件应该像 Python 一样构建。用于编写扩展的 Python 文档应该有一些细节。在 Linux 上,这通常意味着:

$ gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.7 -o yourmod.so yourmod.c

gcc应该可以访问 NumPy C 头文件,因此如果它们未安装在/usr/include/numpy或类似文件中,您可能需要为这些文件传递另一个选项。如果你写的话,你只需要提供 NumPy 标题:

cimport numpy

在您的 Cython 代码中。

这会在同一目录中创建yourmod.so,Python 可以使用普通的import yourmod语句导入该目录。

第一个 Cython 程序

您可以通过下载 Jupyter 笔记本轻松执行本教程的代码。

下面的代码在 numpy 中相当于这个函数:

def compute_np(array_1, array_2, a, b, c):
    return np.clip(array_1, 2, 10) * a + array_2 * b + c

假设 array_1array_2是整数类型的 2D NumPy 数组,abc 是三个 Python 整数。

这个函数使用 NumPy 并且已经非常快了,所以使用 Cython 再次执行它可能有点过分。这是出于演示目的。尽管如此,我们将证明我们以更加冗长的代价实现了比 NumPy 更好的速度和内存效率。

此代码在展开的两个维度上的循环上计算函数结果。它既是有效的 Python 又是有效的 Cython 代码。我将其称为 Python 版本的 compute_py.py 和 Cython 版本的 compute_cy.pyx - Cython 使用 .pyx 作为其文件后缀(但它也可以编译.py文件)。

import numpy as np

def clip(a, min_value, max_value):
    return min(max(a, min_value), max_value)

def compute(array_1, array_2, a, b, c):
    """
 This function must implement the formula
 np.clip(array_1, 2, 10) * a + array_2 * b + c

 array_1 and array_2 are 2D.
 """
    x_max = array_1.shape[0]
    y_max = array_1.shape[1]

    assert array_1.shape == array_2.shape

    result = np.zeros((x_max, y_max), dtype=array_1.dtype)

    for x in range(x_max):
        for y in range(y_max):
            tmp = clip(array_1[x, y], 2, 10)
            tmp = tmp * a + array_2[x, y] * b
            result[x, y] = tmp + c

    return result

对于Linux系统这应该被编译为 compute_cy.so(在 Windows 系统上,这将是 .pyd 文件)。我们运行 Python 进程来测试这两个 Python 版本(从.py -file 导入的版本)和编译的 Cython 模块。

In [1]: import numpy as np
In [2]: array_1 = np.random.uniform(0, 1000, size=(3000, 2000)).astype(np.intc)
In [3]: array_2 = np.random.uniform(0, 1000, size=(3000, 2000)).astype(np.intc)
In [4]: a = 4
In [5]: b = 3
In [6]: c = 9
In [7]: def compute_np(array_1, array_2, a, b, c):
 ...:     return np.clip(array_1, 2, 10) * a + array_2 * b + c
In [8]: %timeit compute_np(array_1, array_2, a, b, c)
103 ms ± 4.16 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [9]: import compute_py
In [10]: compute_py.compute(array_1, array_2, a, b, c)
1min 10s ± 844 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [11]: import compute_cy
In [12]: compute_cy.compute(array_1, array_2, a, b, c)
56.5 s ± 587 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

还没有那么大的差别;因为 C 代码仍然完全符合 Python 解释器的作用(例如,意味着为每个使用的数字分配了一个新对象)。

从命令行调用 Cython 时可以使用-a查看 Python 交互和生成的 C 代码,使用 Jupyter Notebook 时使用 %%cython -a ,或者在使用 setup.py 时使用 cythonize('compute_cy.pyx', annotate=True) 。查看生成的 html 文件,看看即使是最简单的语句也需要什么。你很快就明白了。我们需要给 Cython 更多信息;我们需要添加类型。

添加类型

要添加类型,我们使用自定义 Cython 语法,因此我们现在正在破坏 Python 源兼容性。这是 compute_typed.pyx阅读注释

import numpy as np

# 我需要为数组固定一个数据类型
# 为此,我使用DTYPE变量,赋值为正常的NumPy运行时类型信息对象。
DTYPE = np.intc

# cdef代表这个函数是普通的C函数(因此更快)
# 为了追求最优,我们给参数和返回值都设置了类型
cdef int clip(int a, int min_value, int max_value):
    return min(max(a, min_value), max_value)

def compute(array_1, array_2, int a, int b, int c):

    # 这个 "cdef" 关键词还可以用来在函数里设置变量类型。它只能在最顶层的锁进级别中使用,
    # 在其他地方允许使用它们会造成严重的问题。我们也很期待有方案可以解决这个问题。
    cdef Py_ssize_t x_max = array_1.shape[0]
    cdef Py_ssize_t y_max = array_1.shape[1]

    assert array_1.shape == array_2.shape
    assert array_1.dtype == DTYPE
    assert array_2.dtype == DTYPE

    result = np.zeros((x_max, y_max), dtype=DTYPE)

    # 给你所有的变量设置类型是很重要的。如果你不这么做,你会被警告,但是代码速度会更慢
    # 因为它们会被隐式的定义为Python对象类型
    # 对于 "tmp" 变量,我们希望使用和数组中存储的对象相同的类型,所以我们使用int
    # 因为它和np.intc相对。
    # 注意:这么做的一个严重的副作用就是如果"tmp"溢出了它的类型大小,它会像在C中一样
    # 环绕,而不是像Python一样引发错误。

    cdef int tmp

    # Py_ssize_t is the proper C type for Python array indices.
    cdef Py_ssize_t x, y

    for x in range(x_max):
        for y in range(y_max):

            tmp = clip(array_1[x, y], 2, 10)
            tmp = tmp * a + array_2[x, y] * b
            result[x, y] = tmp + c

    return result

../../_img/compute_typed_html.jpg

此时,查看为compute_cy.pyxcompute_typed.pyx生成的 C 代码。单击这些行以展开它们并查看相应的 C.

特别是看一下for-loops:在compute_cy.c中,这些是要设置的~20 行 C 代码,而在compute_typed.c中使用正常的 C for 循环。

在构建完这个并继续我的(非正式)基准测试后,我得到:

In [13]: %timeit compute_typed.compute(array_1, array_2, a, b, c)
26.5 s ± 422 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

因此添加类型确实使代码更快,但是没有接近 NumPy 的速度?

发生的事情是,在这段代码中花费的大部分时间用在以下行中,并且这些行的执行速度比纯 Python 要慢:

tmp = clip(array_1[x, y], 2, 10)
tmp = tmp * a + array_2[x, y] * b
result[x, y] = tmp + c

那么是什么让这几行代码比纯 Python 版本慢得多?

array_1array_2仍然是 NumPy 数组,所以 Python 对象,并期望 Python 整数作为索引。这里我们传递 C int 值。所以每次 Cython 到达这一行时,它必须将所有 C 整数转换为 Python int 对象。由于这几行代码经常被调用,它完全抵消了之前从range()创建的纯 C 循环的速度优势。

此外,tmp * a + array_2[x, y] * b返回一个 Python 整数,tmp是一个 C 整数,因此 Cython 必须再次进行类型转换。最后,这些类型的转换会加起来。并使我们的计算非常缓慢。但是这个问题可以通过使用内存视图轻松解决。

使用内存视图进行高效索引

仍有两个瓶颈会降低性能,即数组查找和分配,以及 C / Python 类型转换。 [] -运算符仍然使用完整的 Python 操作 - 我们想要做的是直接以 C 速度访问数据缓冲区。

我们需要做的是给ndarray对象的内容设置类型。我们使用内存视图执行此操作。 Cython 文档有一页专门讲内存视图

简而言之,内存视图是 C 结构,可以保存指向 NumPy 数组的数据的指针以及所有必要的缓冲区元数据,以提供有效和安全的访问:维度,步幅,项大小,项类型信息等...它们也支持切片,所以即使 NumPy 数组在内存中不连续,它们也能正常工作。它们可以由 C 整数索引,从而允许快速访问 NumPy 数组数据。

以下是如何声明整数的内存视图:

cdef int [:] foo         # 1维 内存视图
cdef int [:, :] foo      # 2维 内存视图
cdef int [:, :, :] foo   # 3维 内存视图
...                      # 你懂了吧

在我们的示例中,没有数据从 NumPy 数组复制到 memoryview。顾名思义,它只是内存的“视图”。因此,我们可以使用视图result_view进行有效索引,最后返回保存我们操作数据的真实 NumPy 数组result

以下是如何在我们的代码中使用它们:

compute_memview.pyx

import numpy as np

DTYPE = np.intc

cdef int clip(int a, int min_value, int max_value):
    return min(max(a, min_value), max_value)

def compute(int[:, :] array_1, int[:, :] array_2, int a, int b, int c):

    cdef Py_ssize_t x_max = array_1.shape[0]
    cdef Py_ssize_t y_max = array_1.shape[1]

    # array_1.shape is now a C array, no it's not possible
    # to compare it simply by using == without a for-loop.
    # To be able to compare it to array_2.shape easily,
    # we convert them both to Python tuples.
    assert tuple(array_1.shape) == tuple(array_2.shape)

    result = np.zeros((x_max, y_max), dtype=DTYPE)
    cdef int[:, :] result_view = result

    cdef int tmp
    cdef Py_ssize_t x, y

    for x in range(x_max):
        for y in range(y_max):

            tmp = clip(array_1[x, y], 2, 10)
            tmp = tmp * a + array_2[x, y] * b
            result_view[x, y] = tmp + c

    return result

让我们看看现在访问速度有多快。

In [22]: %timeit compute_memview.compute(array_1, array_2, a, b, c)
22.9 ms ± 197 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

请注意这种变化的重要性。我们现在比 Python 的解释版快 3081 倍,比 NumPy 快 4.5 倍。

内存视图也可以与切片一起使用,甚至可以与 Python 数组一起使用。查看 内存视图页面 ,看看它们能为您做些什么。

进一步调整索引

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

  1. 进行边界检查。
  2. 检查负指数并正确处理。上面的代码是明确编码的,因此它不使用负索引,并且(希望)总是在边界内访问。

使用装饰器,我们可以停用这些检查:

...
cimport cython
@cython.boundscheck(False)  # Deactivate bounds checking
@cython.wraparound(False)   # Deactivate negative indexing.
def compute(int[:, :] array_1, int[:, :] array_2, int a, int b, int c):
...

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

In [23]: %timeit compute_index.compute(array_1, array_2, a, b, c)
16.8 ms ± 25.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

我们比 NumPy 版本(6.2x)更快。 NumPy 编写得很好,但不会执行懒操作,导致在内存中进行大量的中间复制操作。我们的版本具有非常高的内存效率和缓存友好性,因为我们可以在一次运行中执行数据操作。

警告: 提高速度需要付出成本。特别是将类型对象(如我们的示例代码中的 array_1array_2result_view )设置为 None 会很危险。将这些对象设置为 None 是完全合法的,你唯一可以做的就是检查它们是否为 None。所有其他用途(属性查找或索引)都可能会破坏或损坏数据(而不是像在 Python 中那样引发异常)。

实际的规则有点复杂,但主旨很明确:不要在不知道是否为None的情况下,使用类型化的对象。

将 NumPy 数组声明为连续

为了获得额外的速度增益,如果您知道您提供的 NumPy 阵列在内存中是连续的,则可以将 内存视图 声明为连续的。

我们举例说明一个有 3 个维度的数组。如果你想给 Cython 提供数据是 C-连续的 这个信息,你必须像这样声明内存视图:

cdef int [:,:,::1] a

如果你想给 Cython 提供数据是 Fortran 连续的信息,你必须像这样声明内存视图:

cdef int [::1, :, :] a

(译者注:内存视图这一章这里提到了, ::1 表示全部维度连续,对于C意味着之前的所有,对于Fortran意味着之后的所有)

如果所有这些对你没有意义,你可以跳过这一部分,声明数组是连续的约束你的函数的用法,因为它拒绝数组切片作为输入。如果您仍然想了解连续数组的全部内容,可以在 StackOverflow 上看到这个答案。

为了给出数字,这里是通过将 memoryviews 声明为连续来获得的速度增益:

In [23]: %timeit compute_contiguous.compute(array_1, array_2, a, b, c)
11.1 ms ± 30.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

我们现在的速度比 NumPy 版快 9 倍,比纯 Python 版快 6300 倍!

使函数更清洁

声明类型会令您的代码显得非常啰嗦。如果您不介意由 Cython 来推断变量的 C 类型,可以使用文件顶部的infer_types=True编译器指令。它会为你节省很多敲代码的时间。

请注意,由于类型声明必须发生在顶部缩进级别,因此 Cython 不会推断在其他缩进级别中第一次声明的变量类型。它会改变我们代码的含义。这就是为什么我们仍然必须手动声明 tmpxy 变量的类型。

实际上,在使用融合类型时,手动给出 tmp 变量的类型会很有用。

# cython: infer_types=True
import numpy as np
cimport cython

DTYPE = np.intc

cdef int clip(int a, int min_value, int max_value):
    return min(max(a, min_value), max_value)

@cython.boundscheck(False)
@cython.wraparound(False)
def compute(int[:, ::1] array_1, int[:, ::1] array_2, int a, int b, int c):

    x_max = array_1.shape[0]
    y_max = array_1.shape[1]

    assert tuple(array_1.shape) == tuple(array_2.shape)

    result = np.zeros((x_max, y_max), dtype=DTYPE)
    cdef int[:, ::1] result_view = result

    cdef int tmp
    cdef Py_ssize_t x, y

    for x in range(x_max):
        for y in range(y_max):

            tmp = clip(array_1[x, y], 2, 10)
            tmp = tmp * a + array_2[x, y] * b
            result_view[x, y] = tmp + c

    return result

我们现在进行速度测试:

In [24]: %timeit compute_infer_types.compute(array_1, array_2, a, b, c)
11.5 ms ± 261 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

瞧,速度没变。

更泛型的代码

所有这些速度提升都很好,但添加类型会限制我们的代码。目前,这意味着我们的函数只能与具有 np.intc 类型的 NumPy 数组一起使用。是否可以使我们的代码适用于多种 NumPy 数据类型?

是的,借助称为融合类型的新功能。您可以在 文档 的这一部分了解更多相关信息。它类似于 C ++的模板。它在编译时生成多个函数声明,然后根据提供的参数类型在运行时选择正确的声明。通过比较 if 条件中的类型,还可以根据特定数据类型执行完全不同的代码路径。

在我们的例子中,由于我们不再访问 NumPy 的输入数组的 dtype,我们使用那些if-else语句来知道我们应该为输出数组使用什么 NumPy 数据类型。

在这种情况下,我们的函数现在适用于整数,双精度和浮点数。

# cython: infer_types=True
import numpy as np
cimport cython

ctypedef fused my_type:
    int
    double
    long long

cdef my_type clip(my_type a, my_type min_value, my_type max_value):
    return min(max(a, min_value), max_value)

@cython.boundscheck(False)
@cython.wraparound(False)
def compute(my_type[:, ::1] array_1, my_type[:, ::1] array_2, my_type a, my_type b, my_type c):

    x_max = array_1.shape[0]
    y_max = array_1.shape[1]

    assert tuple(array_1.shape) == tuple(array_2.shape)

    if my_type is int:
        dtype = np.intc
    elif my_type is double:
        dtype = np.double
    elif my_type is cython.longlong:
        dtype = np.longlong

    result = np.zeros((x_max, y_max), dtype=dtype)
    cdef my_type[:, ::1] result_view = result

    cdef my_type tmp
    cdef Py_ssize_t x, y

    for x in range(x_max):
        for y in range(y_max):

            tmp = clip(array_1[x, y], 2, 10)
            tmp = tmp * a + array_2[x, y] * b
            result_view[x, y] = tmp + c

    return result

我们可以检查输出类型是否正确:

>>>compute(array_1, array_2, a, b, c).dtype
dtype('int32')
>>>compute(array_1.astype(np.double), array_2.astype(np.double), a, b, c).dtype
dtype('float64')

We now do a speed test:

In [25]: %timeit compute_fused_types.compute(array_1, array_2, a, b, c)
11.5 ms ± 258 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

在编译时创建更多版本的函数。因此有意义的是,与以前一样使用整数执行此函数的速度不会改变。

使用多个线程

Cython 支持 OpenMP。它周围也有一些很好的包装器,比如函数 prange()。您可以在 使用并行 中查看有关 Cython 和并行性的更多信息。由于我们进行元素操作,因此我们可以轻松地在多个线程之间分配工作。重要的是不要忘记将正确的参数传递给编译器以启用 OpenMP。使用 Jupyter 笔记本时,你应该像这样使用单元格魔法指令:

%%cython --force
# distutils: extra_compile_args=-fopenmp
# distutils: extra_link_args=-fopenmp

必须释放 GIL(参见 释放 GIL),这就是我们声明 clip() 函数 nogil 的原因。

# tag: openmp
# You can ignore the previous line.
# It's for internal testing of the cython documentation.

# distutils: extra_compile_args=-fopenmp
# distutils: extra_link_args=-fopenmp

import numpy as np
cimport cython
from cython.parallel import prange

ctypedef fused my_type:
    int
    double
    long long

# We declare our plain c function nogil
cdef my_type clip(my_type a, my_type min_value, my_type max_value) nogil:
    return min(max(a, min_value), max_value)

@cython.boundscheck(False)
@cython.wraparound(False)
def compute(my_type[:, ::1] array_1, my_type[:, ::1] array_2, my_type a, my_type b, my_type c):

    cdef Py_ssize_t x_max = array_1.shape[0]
    cdef Py_ssize_t y_max = array_1.shape[1]

    assert tuple(array_1.shape) == tuple(array_2.shape)

    if my_type is int:
        dtype = np.intc
    elif my_type is double:
        dtype = np.double
    elif my_type is cython.longlong:
        dtype = np.longlong

    result = np.zeros((x_max, y_max), dtype=dtype)
    cdef my_type[:, ::1] result_view = result

    cdef my_type tmp
    cdef Py_ssize_t x, y

    # We use prange here.
    for x in prange(x_max, nogil=True):
        for y in range(y_max):

            tmp = clip(array_1[x, y], 2, 10)
            tmp = tmp * a + array_2[x, y] * b
            result_view[x, y] = tmp + c

    return result

我们可以以最小的努力获得可观的速度提升:

In [25]: %timeit compute_prange.compute(array_1, array_2, a, b, c)
9.33 ms ± 412 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

我们现在比纯 Python 版本快 7558 倍,比 NumPy 快 11.1 倍!

然后去哪儿?



回到顶部