《Cython系列》7. Cython、numpy、以及类型化memoryview

楔子

Cython的两个优秀的品质就是它的广度和成熟度，可以编译所有的Python代码，并且将C的速度代入了Python，并且还能轻松的和C、C++集成。而本篇文章的任务就是完善Cython的功能，并介绍Cython的阵列特性，比如：对numpy数组的深入支持。

我们已经知道，Cython可以很好的支持list、dict、set、tuple等内置容器，这些容器非常容易使用，可以包含任意的Python对象，并且对对象的查找、分配、检索都进行了高度的优化。尽管list类型和dict类型之间实现存储和检索的方式有很大差异，但是从实现角度来讲，所有的容器都有一个共同点：它们存储的都是对Python对象的一个引用。如果我们有一个包含100万个整型的列表，那么在C中，每一个元素都是一个Python *，虽然整型在底层对应PyLongObject，但是会转成PyObject *存在列表里面，因为C中的数组只能存放相同类型的元素(指针也是相同类型的指针)，然后用的时候再转回对应的类型，因为PyObject里面存储了指向对应类型的指针。所以Python中的列表可以存储任意类型的变量，因为它们都是一个PyObject *。但我们说，由于每次运算的时候都要进行类型转化，所以效率低，尽管我们知道都是整型，但是Python不会做这种假设。

因此对于同质容器(只包含一种固定的类型)，可以在存储开销和性能方面做得更好，比如numpy中的数组。numpy中的数组就有一个准确的类型，这样的话在存储和计算的时候可以表现的更加优秀，我们举个栗子：

import numpy as np

# 创建一个整型数组
# 可以指定dtype为: int, int8, int16, int32, int64
# 或者是: uint, uint8, uint16, uint32, uint64
arr1 = np.zeros((3, 3), dtype="uint32")
print(arr1)
"""
[[0 0 0]
 [0 0 0]
 [0 0 0]]
"""
try:
    arr1[0, 0] = "xx"
except Exception as e:
    print(e)  # invalid literal for int() with base 10: 'xx'
# 我们看到报错了，因为已经规定了arr是一个整型数组，那么在存储和计算都是按照整型来处理的
# 既然是整型数组，那么赋值一个字符串是不允许的

arr1[0, 0] = -123
print(arr1)
"""
[[4294967173          0          0]
 [         0          0          0]
 [         0          0          0]]
"""
# 因为是uint32, 只能存储正整数, 所以结果是 uint32的最大值 - 123
print((2 << 31) - 123)  # 4294967173


# 创建一个浮点数组, 可以指定dtype为: float, float16, float32, float64
arr2 = np.zeros((3, 3), dtype="float")
print(arr2)
"""
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
"""


# 创建一个字符串类型, dtype可以为: U, str
# 如果是U, 那么还可以加上一个数值, 比如: U3, 表示最多存储3个字符。
# 并且还可以通过<U或者>U的方式来指定是小端存储或者大端存储
arr3 = np.zeros((3, 3), dtype="U3")
print(arr3)
"""
[['' '' '']
 ['' '' '']
 ['' '' '']]
"""
arr3[0, 0] = "古明地觉"
print(arr3)
"""
[['古明地' '' '']
 ['' '' '']
 ['' '' '']]
"""
# 我们看到被截断了，并且截断是按照字符来的，不是按照字节


# 创建一个Python对象, 注意：没有tuple、list、dict等类型
# 特定的类型只有整型、浮点型、字符串，至于其它的类型统统用object表示
# 可以指定dtype="O"
arr4 = np.zeros((3, 3), dtype="O")
print(arr4)
"""
[[0 0 0]
 [0 0 0]
 [0 0 0]]
"""
# 虽然打印的也是0，但它是一个object类型
print(arr4.dtype)  # object
# 我们可以使用empty创建
print(np.empty((3, 3), dtype="O"))
"""
[[None None None]
 [None None None]
 [None None None]]
"""

而实现同质容器是通过缓冲区的方式，它允许我们将连续的简单数据使用单个数据类型表示，而支持新缓冲区协议的numpy数组在Python中则是使用最广泛的数组类型。而numpy数组的速度有多快不用我多说，所以我们可以将缓冲区视为简化的numpy数组。

有效地使用缓冲区通常是从Cython代码中获得C性能的关键，而幸运的是，Cython使处理缓冲区变得非常的容易，它对新缓冲区协议和numpy数组有着一流的支持。

新缓冲区协议的威力

新缓冲区协议是一个C级协议，Python对象可以实现该协议，并且不会影响它们在Python级别的接口。所有的Python3版本都支持该协议，它定义了一个具有数据缓冲区和元数据的C级结构体，并用它来描述缓冲区的布局、数据类型和读写权限，并且还定义了支持协议的对象所必须实现的api。

新缓冲区协议最重要的特性就是它能以不同的方式来表示相同的底层数组，它允许numpy数组、几个Python的内置类型、和Cython级别的数组共享相同的数据，而无需赋值。并且使用Cython，我们还可以轻松地扩展缓冲区协议去处理来自外部库的数据。

关于协议的细节，在Python的C api参数手册中有详细的说明，但是我们不需要在这里讨论，因为它很复杂。但幸运的是，Cython允许在不了解协议的情况下使用缓冲区，可以在不复制的情况下高效的访问底层数据，从而减少开销，有这一点就足够了。

那都有哪些类型实现了新缓冲区协议呢？

numpy的ndarray

Python中最知名、使用最广泛的numpy包中有一个ndarray对象，它是一个支持新缓冲区协议的有效Python缓冲区。

Python2中的str

Python2中的str也实现了该协议，但是Python3的str则没有。

Python中的bytes和bytearray

既然Python2中的str实现了该协议，那么代表Python3的bytes也实现了。当然还有bytearray。

标准库array中的array

Python标准库中有一个array，里面的array也实现了该协议，但是我们用的不多。

标准库ctypes中的array

这个我们用的也不多。

其它的第三方数据类型

比如第三方库PIL，用于处理图片的，将图片读取进来得到的对象的类型也实现了新缓冲区协议。当然这个很好理解，因为它们读取进来可以直接转成numpy的ndarray。

内存视图的类型

还有一个内置的Python类型是memoryview(内存视图)，它存在的唯一目的就是在Python级别表示C级的缓冲区。我们可以像memoryview传递一个实现了新缓冲区协议的对象(比如：bytes对象)，来创建一个memoryview对象。

b = b"komeiji satori"
m = memoryview(b)

# 此时和m和b之间是共享内存的
print(m)  # <memory at 0x000001C3A54EFA00>

# 可以通过索引、切片的方式访问
print(m[0], m[-1])  # 107 105
# 通过索引获取得到的是数字，切片获取得到的还是一个memoryview对象
print(f"{m[0]:c}", f"{m[-1]:c}")  # k i
print(m[0: 2])  # <memory at 0x00000229D637F880>
print(m[0: 2][-1], f"{m[0: 2][-1]:c}")  # 111 o

我们可以通过切片的方式对memoryview进行任意截取，通过这种方式，灵活性就变得非常高。

那么能不能修改呢？答案是：因为bytes对象不可以修改，所以memoryview也不可以修改，也就是只读的。

b = b"komeiji satori"
m = memoryview(b)

print(m.readonly)  # True
try:
    m[0] = "K"  
except Exception as e:
    print(e)  # cannot modify read-only memory

bytes对应的缓冲区是不可以修改的，所以如果想修改的话，我们应该使用bytearray。

当然我们还可以传入一个numpy的ndarray等等，这里就不说了，memoryview我个人用的很少。

类型化的内存视图

Cython有一个C级类型，即：类型化的memoryview，它在概念上和Python中的memoryview重叠、并且在其基础上展开。正如其名所示，类型化的memoryview用于查看(共享)来自缓冲区对象的数据。并且类型化memoryview在C级别操作，所以它有着最小的Python开销，因此非常高效，而且比直接使用C级缓冲区跟方便。此外类型化的memoryview是为了和缓冲区一起工作而被设计的，因此它可以有效支持任何来自缓冲区的对象，从而允许在不复制的情况下共享缓冲区数据。

举个栗子

假设我们想在Cython中有效地处理一维数据的缓冲区，而不关心如何在Python级别创建数据，我们只是想以一种有效的方式访问它。

def summer(double[:] numbers):
    """
    code
    """

double[:] numbers声明了numbers是一个类型化的memoryview，而double指定了该memoryview的基本类型，[:]则表明这是一个一维的memoryview。

当我们从Python中调用summer函数时，会传入一个Python对象，该对象会在调用函数的时候隐式的分配给参数numbers。当分配个类型化的memoryview时，该memoryview会尝试访问该对象的数据缓冲区。如果该对象没有提供相应的缓冲区(不支持该缓冲区协议)，那么会引发一个ValueError；如果支持的话新缓冲区协议，那么将会提供一个C级缓冲区供memoryview使用。

def summer(double[:] numbers):
    cdef double s = 0
    cdef double number
    # memoryview可以像Python迭代器一样进行遍历
    for number in numbers:
        s += number
    return s

>>> import cython_test
>>> import numpy as np
>>> # 必须传递支持缓存区协议的对象
>>> cython_test.summer(np.array([1.2, 2.3, 3.4, 4.5]))
11.4
>>> # 可以直接传入数组，也可以传入memoryview对象
>>> cython_test.summer(memoryview(np.array([1.2, 2.3, 3.4, 4.5])))
11.4
>>> # 但是传递列表是不行的，因为它不支持缓冲区协议
>>> cython_test.summer([1.2, 2.3, 3.4, 4.5])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "cython_test.pyx", line 1, in cython_test.summer
    def summer(double[:] numbers):
  File "stringsource", line 658, in View.MemoryView.memoryview_cwrapper
  File "stringsource", line 349, in View.MemoryView.memoryview.__cinit__
TypeError: a bytes-like object is required, not 'list'
>>> 

C级访问类型化memoryview数据

类型化memoryview是为C风格的访问而设计的，没有开销，因此也可以用另一种方式去遍历number。

def summer(double[:] numbers):
    cdef double s = 0
    cdef int i, N
    # 调用shape拿到其长度
    N = numbers.shape[0]
    for i in range(N):
        s += numbers[i]
    return s

>>> import cython_test
>>> import numpy as np
>>> 
>>> cython_test.summer(np.array([1.2, 2.3, 3.4, 4.5]))
11.4

这个版本会有更好的性能：对于百万元素数组来说大约是1毫秒，因为它直接在Python/C api类型化的情况下生成一个整型索引，这是我们一个加速的来源，当然我们还可以做得更好。

用安全换取性能

每次我们访问memoryview对象时，Cython都会检测索引是否越界内，如果越界，那么Cython将引发一个IndexError。而且Cython也允许我们像Python一样通过负数索引对memoryview对象进行访问。

对于我们上面的summer函数，我们在访问内部memoryview对象之前就已经获取了它的元素个数，因为我们在遍历的时候用于不会越界，所以我们可以指示Cython关闭这些检查以获取更高的性能。而关闭检查可以使用上下文的方式：

from cython cimport boundscheck, wraparound

def summer(double[:] numbers):
    cdef double s = 0
    cdef int i, N
    N = numbers.shape[0]
    # 关闭检查
    with boundscheck(False), wraparound(False):
        for i in range(N):
            s += numbers[i]
        return s

>>> import cython_test
>>> import numpy as np
>>> 
>>> cython_test.summer(np.array([1.2, 2.3, 3.4, 4.5]))
11.4
>>> 

我们通过上下文管理的方式，在其内部对memoryview对象进行索引访问的时候，会关闭检测，当然这些修改仅在上下文管理器内部有效。得到的结果是性能会有小幅度提高(当然我们这里数据很少，看不出来)，代码生成效率也更高。但效率提升的后果就是必须由我们来确保索引不会越界、并且不可以使用负数索引，否则的话可能会导致段错误(非常危险，不管程序崩溃，解释器就直接退出了)。因此如果没有百分之百的把握，不要关闭检查。

如果我们想关闭整个函数内部的检查，那么可以使用装饰器的形式。

from cython cimport boundscheck, wraparound

@boundscheck(False)
@wraparound(False)
def summer(double[:] numbers):
    cdef double s = 0
    cdef int i, N
    N = numbers.shape[0]
    for i in range(N):
        s += numbers[i]
    return s

如果想关闭全局的边界检测，那么可以使用注释的形式。

# cython: boundscheck=False
# cython: wraparound=False

所以关闭边界检测有多种方式，但是不同的方式对应不同的作用域。

所以有了以上这些方法，我们的summer函数的性能，和numpy中sum函数的性能是几乎等价的。

但是还没完，我们还有更多的细节要将。

声明类型化memoryview

当我们声明一个类型化memoryview时，我们可以控制很多的属性。

元素类型

类型化memoryview的元素类型可以是int、long、float、double、complex等数值标量类型；也可以是ctypedef起的别名；也可以是使用cdef struct声明的结构体类型。

维度

类型化memoryview最多可以有7个维度，我们之前声明一个一维的，使用的是double[:]这种形式，如果是3位，那么写成double[:, :, :]即可。当然类型不一定是double，也可以是其它的。

直接或间接访问

直接访问是默认的，它涵盖了几乎所有的情况，它指定对应的维度可以通过索引的方式直接访问底层数据。如果将某个维度指定为间接访问，则底层缓冲区将存储一个指向数组剩余部分的指针，该指针必须在访问时解引用(因此是间接访问)。而numpy不支持间接访问，所以我们不使用这种访问规范，因此直接访问是默认的。事实上，官方一般设置为默认的都是最好的。

下面举例说明：

import numpy as np

# 这是最灵活的声明方式，可以从任何int类型的二维ndarray中获取缓冲区
def func(int[:, :] ages):
    # 直接打印是一个memoryview对象，这里转成ndarray
    print(np.array(ages))

>>> import cython_test
>>> import numpy as np
>>> # 这里类型要匹配，C中的int对应numpy中的int32，而numpy的整型默认是int64
>>> # 因为memoryview对类型要求很严格
>>> arr = np.random.randint(1, 10, (3, 3), dtype="int32")
>>> cython_test.func(arr)
[[9 3 8]
 [8 3 2]
 [1 1 4]]
>>> # 这里的order="F"表示按照Fortran的规则存储
>>> arr1 = np.array([[1, 2, 3], [2, 3, 3], [3, 2, 1]], dtype="int32", order="F")
>>> cython_test.func(arr1)
[[1 2 3]
 [2 3 3]
 [3 2 1]]
>>> 

numpy的ndarray作为变量类型

我们将numpy的ndarray作为类型的时候，Cython赋予了它额外的属性。

# 这里一定要使用cimport
cimport numpy as np  # 或者from numpy cimport ndarray


def func1(np.ndarray array):
    print(array)

能够作为Cython中的类型，那么它要么是扩展类型、要么是内置类型，所以这里一定要使用cimport，如果是import，那么对应的ndarray不可以作为类型。然后我们进行编译，但是注意：编译的时候需要多指定一个参数

from distutils.core import Extension, setup
from Cython.Build import cythonize
import numpy as np

ext = Extension("cython_test", ["cython_test.pyx"])

setup(
    ext_modules=cythonize(ext, language_level=3),
    # 如果没有这个参数，那么编译是会报错：fatal error: numpy/arrayobject.h: No such file or directory
    # 所以我们必须把响应的头文件包含进来
    include_dirs=[np.get_include()]  
)

上面的函数表示接收一个nd.array类型的数组，然后打印，至于这个数组里面的元素是什么类型则没有要求。

>>> import cython_test
>>> import numpy as np
>>> 
>>> cython_test.func1(np.array([1, 2, 3]))
[1 2 3]
>>> cython_test.func1(np.array([[1, 2, 3]]))
[[1 2 3]]
>>> cython_test.func1(np.array(["a", "b", "c"]))
['a' 'b' 'c']
>>> 

我们可以传递任意的array，但如果我们希望限制array内部元素的类型以及维度，这个时候该怎么做呢？

Cython为numpy提供了专门的方法：比如希望参数是一个接收类型为int64、维度为2的数组，就可以使用ndarray[long, ndim=2]这种方式，我们演示一下。

cimport numpy as np


def func1(np.ndarray[long, ndim=2] array):
    print(array)


def func2(np.ndarray[double, ndim=1] array):
    print(array)


def func3(np.ndarray[object, ndim=1] array):
    print(array)

>>> import cython_test
>>> import numpy as np
>>> # 类型需要匹配，参数是long，所以这里是int64
>>> cython_test.func1(np.array([[1, 2], [3, 4]]))
[[1 2]
 [3 4]]
>>> # 维度必须匹配
>>> cython_test.func1(np.array([1, 2, 3, 4]))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "cython_test.pyx", line 4, in cython_test.func1
    def func1(np.ndarray[long, ndim=2] array):
ValueError: Buffer has wrong number of dimensions (expected 2, got 1)
>>> cython_test.func3(np.array(["a", "b", "c"], dtype="O"))
['a' 'b' 'c']
>>> 

总结

这一张的内容我们只介绍了一部分，但是个人觉得对于日常工作来说差不多了，有兴趣深入的话可以参考相应的书籍。