楔子
Cython 的两个优秀的品质就是它的广度和成熟度，可以编译所有的 Python 代码，并且将 C 的速度引入到了 Python 中，此外还能轻松的和 C、C++ 集成。而本篇文章的任务就是完善 Cython 的功能，并介绍 Cython 的阵列特性，比如：对 Numpy 数组的深入支持。
我们已经知道，Cython 可以很好地支持列表、元组、字典等内置容器，这些容器非常容易使用，可以包含指向任意 Python 对象的变量，并且对 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 数组的高性能是有目共睹的。
有效地使用缓冲区通常是从 Cython 代码中获得 C 性能的关键，而幸运的是，Cython 使处理缓冲区变得非常的容易，它对缓冲区协议和 Numpy 数组有着一流的支持。
什么是缓冲区协议？
下面来说一下缓冲区协议，缓冲区协议是一个 C 级协议，它定义了一个具有数据缓冲区和元数据的 C 级结构体，并用它来描述缓冲区的布局、数据类型和读写权限，并且还定义了支持协议的对象所必须实现的 API。
而实现了该协议的 Python 对象之间可以向彼此暴露自身的原始字节数组，这在科学计算中非常有用，因为在科学计算中我们经常会使用诸如 Numpy 之类的包来存储和操作大型数组，因为要对数据做各种各样的变换，所以难免要涉及到数组的拷贝。而使用缓冲区协议，那么数组之间就可以不用拷贝了，而是共享同一份数据缓冲区，这些数组都可以看成是该缓冲区的一个视图，那么也意味着操作任何一个数组都会影响其它数组。
那么都有哪些类型实现了缓冲区协议呢？
Numpy 中的 ndarray
Python 中最知名、使用最广泛的 Numpy 包里面有一个 ndarray 对象，它是一个支持缓冲区协议的有效 Python 缓冲区。
Python2 中的 str
Python2 中的 str 也实现了该协议，但是 Python3 的 str 则没有。
Python3 中的 bytes 和 bytearray
既然 Python2 中的 str 实现了该协议，那么代表 Python3 的 bytes 也实现了，当然还有 bytearray。
标准库 array 中的 array
Python 标准库中有一个 array 模块，里面的 array 也实现了该协议，但是我们用的不多。
标准库 ctypes 中的 array
这个我们用的也不多。
其它的第三方数据类型
比如第三方库 PIL，用于处理图片的，将图片读取进来得到的对象也实现了缓冲区协议。当然这个很好理解，因为它们读取进来可以直接转成 Numpy 的 ndarray。
缓冲区协议最重要的特性就是它能以不同的方式来表示相同的底层数组，它允许 Numpy 数组、几个 Python 的内置类型、标准库的数组之间共享相同的数据，而无需再拷贝。当然 Cython 级别的数组也是可以的，并且使用 Cython，我们还可以轻松地扩展缓冲区协议去处理来自外部库的数据（后面说）。
我们举个例子，看看不同类型的数据之间如何共享内存：
import array import numpy as np """ 'b' signed integer 'B' unsigned integer 'u' Unicode character 'h' signed short 'H' unsigned short 'i' signed int 'I' unsigned int 'l' signed long 'L' unsigned long 'q' signed long long 'Q' unsigned long long 'f' float 'd' double """ arr = array.array("i", range(6)) print(arr) """ array('i', [0, 1, 2, 3, 4, 5]) """ # array（数组）是标准库array中提供的数据结构，它不是 Python 内置的 # 数组不同于列表，因为数组里面存储的都是连续的整数块 # 它的数据存储要比列表紧凑得多，因此一些操作也可以更快的执行 # 基于 arr 创建 Numpy 的数组 np_arr = np.asarray(arr) print(np_arr) """ [0 1 2 3 4 5] """ # 修改 Numpy 数组 np_arr[0] = 123 # arr 也被改变了，因为它们共享内存 print(arr) """ array('i', [123, 1, 2, 3, 4, 5]) """
Python 提供的数组使用的不是特别多，而 Numpy 的数组使用的则是非常广泛，并且支持的操作非常丰富。但是这两种数组都实现了缓冲区协议，因此可以共享同一份数据缓冲区，它们在转化的时候是不用复制原始数据的。所以 np_arr 在将第一个元素修改之后，打印 arr 也发生了变化。
当然 np.asarray 等价于不拷贝的 np.array：
import array import numpy as np arr = array.array("i", range(6)) # np.array 内部有一个 copy 参数，默认是 True # 也就是会将原始数组拷贝一份 np_arr1 = np.array(arr) np_arr1[0] = 123 # 此时 arr 是没有变化的，因为操作的不是同一个数组 print(arr) # array('i', [0, 1, 2, 3, 4, 5]) # 不进行拷贝，则共享缓冲区，等价于 asarray np_arr2 = np.array(arr, copy=False) np_arr2[0] = 123 # 因此结果变了 print(arr) # array('i', [123, 1, 2, 3, 4, 5])
问题来了，如果我们将 array 换成 list 的话会怎么样呢？
import numpy as np s = [1, 2, 3] np_arr = np.asarray(s) np_arr[0] = 123 print(s) # [1, 2, 3]
因为列表不支持、或者说没有实现缓冲区协议，所以 Numpy 没办法与之共享数据缓冲区，因而只能将数据拷贝一份。
可能有人觉得以现如今的硬件来说，根本不需要考虑内存占用方面的问题，但即便如此，共享内存也是非常有必要的。因为在科学计算中，大部分的经典算法都是采用编译型语言实现的，像我们熟知的 scipy 本质上就是基于 NetLib 实现的一些包装器，NetLib 才是负责提供大量高效算法的工具箱，而这些高效算法几乎都是采用 Fortran 和 C 编写的。Python 能够和这些编译库（NetLib）共享本地数据对于科学计算而言非常重要，这也是 Numpy 被开发的一个重要原因，更是缓冲区协议被提出、并添加到 Python 中的原因。
在这里我们提一下 PyPy，我们知道它是用 CPython 编写的 Python 解释器，它的速度要比 Python 快很多，但是对于使用 Python 进行科学计算的人来说却反而没什么吸引力。原因是在科学计算时所使用的算法实际上都是采用 Fortran 和 C 等语言编写的、并被编译成库的形式，Python 只是负责对这些库进行包装、提供一个友好的接口，因此这意味着 Python 能够与之进行很好的交互。而 PyPy 还无法做到这一点，因此现在用的解释器基本都是 CPython，至于 PyPy 引入 JIT（即时编译）所带来的性能收益实际上用处不大。
Python 能成为有效的科学计算平台，主要得益于缓冲区协议的实现和 Numpy。
缓冲区协议长什么样子？
Python 的缓冲区协议本质上是一个结构体，它为多维数组定义了一个非常灵活的接口，我们看一下底层实现，源码位于 object.h 中。
typedef struct bufferinfo { void *buf; PyObject *obj; Py_ssize_t len; Py_ssize_t itemsize; int readonly; int ndim; char *format; Py_ssize_t *shape; Py_ssize_t *strides; Py_ssize_t *suboffsets; void *internal; } Py_buffer;
以上就是缓冲区协议的底层定义，我们来解释一下里面的成员都代表什么含义，至于如何实现一会再说。
void *buf
实现了缓冲区协议的对象的内部缓冲区（指针），可以被多个不同的对象共享，只要这些对象都实现了缓冲区协议。
PyObject *obj
实现了缓冲区协议的对象（指针），比如 ndarray 对象、bytes 对象等等。
Py_ssize_t len
不要被名字误导了，这里表示缓冲区的总大小。比如一个 shape 为 (3, 4, 5) 的数组，存储的元素是 8 字节的 int64，那么这个 len 就等于 3 * 4 * 5 * 8。
Py_ssize_t itemsize
缓冲区存储的元素都是同质的，每一个元素都占用相同的字节，而 itemsize 就表示每个元素所占的大小。比如缓冲区存储的元素是 int64，那么 itemsize 就是 8。
int readonly
缓冲区是否是只读的，为 0 表示可读写，为 1 表示只读。
int ndim
维度，比如 shape 为 (3, 4, 5) 的数组，那么 ndim 就是 3。注意：如果 ndim 为 0，表示 buf 指向的缓冲区代表的只是一个标量，这种情况下，字段 shape, strides, suboffsets 都必须为 NULL。
而且维度最大不超过 64，但在 Numpy 里面支持的最大维度是 32。
char *format
格式化字符，用于描述缓冲区元素的类型，和 Python 标准库 struct 使用的 format 是一致的。比如：i 表示 C 的 int，L 表示 C 的 unsigned long 等等。
Py_ssize_t *shape
这个很好理解，等同于 Numpy array 的 shape，只不过在 C 里面是一个数组。
Py_ssize_t *strides
维度为 ndim 的数组，里面的值表示在某个维度下，从一个元素到下一个元素所需要跳跃的字节数。举个栗子，假设有一个 shape 为 (10, 20, 30) 的数组，里面的元素是 int64，那么 strides 就是 (4800, 240, 8)。
因为有三个维度：对于第一个维度而言，每一个元素都是 (20, 30) 的二维数组，所以当前元素和下一个元素的地址差了 20 * 30 * 8 = 4800 个字节；对于第二个维度而言，每一个元素都是 (30,) 的一维数组，所以当前元素和下一个元素的地址差了 30 * 8 = 240 个字节；对于第三个维度而言，每一个元素都是一个标量，所以当前元素和下一个元素的地址差了 8 个字节。
根据 strides 我们可以引出一个概念：full strided array，直接解释的话比较费劲，我们用代码说明。
import numpy as np arr1 = np.array(range(10), dtype="int64") print(arr1.strides) # (8,) arr2 = arr1[:: 2] print(arr2.strides) # (16,)
显然 arr1 和 arr2 是共享缓冲区的，也就是说它们底层的 buf 指向的是同一块内存，但它们的 strides 不同。因此 arr1 从一个元素到下一个元素需要跳跃 8 字节，但是 arr2 则是跳跃 16 个字节，否则就无法实现步长为 2 了。
假设把步长从 2 改成 3，那么 arr2 的 strides 显然就变成了 (24,)，所以此刻你应该对 Numpy 数组有更深的认识了。使用切片，本质上是通过改变 strides 来实现跨步访问，但仍然共享同一个缓冲区。
但 arr2 只有一个维度，所以 strides 的元素个数为 1，里面的 16 表示数组 arr2 从一个元素到下一个元素所跳跃的字节数。但是问题来了，arr2 里面的元素大小只有 8 字节，所以像这种元素大小和对应的 strides 不相等的数组，我们称之为 full strided 数组。
对于多维数组也是一样，我们举个例子：
import numpy as np arr = np.ones((10, 20, 30), dtype="int64") print(arr.strides) # (4800, 240, 8) arr2 = arr[:: 2] # arr2 是 full strided，因为在第一个维度中 # 一个元素到下一个元素应该需要 4800 个字节 # 但是 arr2 的 strides 的第一个元素是 9600 # 因为不相等，所以是 full strided print(arr2.strides) # (9600, 240, 8) arr3 = arr[:, :: 2] # arr3 是 full strided，因为在第二个维度中 # 一个元素到下一个元素应该需要 240 个字节 # 但是 arr3 的 strides 的第二个元素是 480， # 因为不相等，所以是 full strided print(arr3.strides) # (4800, 480, 8) arr4 = arr[:, :, :: 2] # arr4 是 full strided，因为在第三个维度中 # 一个元素到下一个元素应该需要 8 个字节 # 但是 arr4 的 strides 的第三个元素是 16 # 因为不相等，所以是 full strided print(arr4.strides) # (4800, 240, 16)
说白了，只要任意维度出现了数组的跨步访问、且步长不为 1，那么这个数组就是 full strided 数组。之所以要说这个 full strided，是因为后面会用到。
代码演示缓冲区协议
再来看一下缓冲区协议长什么样子？
typedef struct bufferinfo { void *buf; PyObject *obj; Py_ssize_t len; Py_ssize_t itemsize; int readonly; int ndim; char *format; Py_ssize_t *shape; Py_ssize_t *strides; Py_ssize_t *suboffsets; void *internal; } Py_buffer;
Py_buffer 内部的 obj 指向了实现缓冲区协议的对象，内部的 buf 则指向了缓冲区本身。而缓冲区本质上就是一个一维数组，负责存储具体的数据，可以被任意多个对象共享。
像 Numpy 的数组，在拷贝的时候只会将 Py_buffer 拷贝一份，但是内部的 buf 成员指向的缓冲区则不会拷贝。
import numpy as np #Py_buffer -> buf 指向了缓冲区 #Py_buffer -> shape 为 (6,) arr1 = np.array([3, 9, 5, 7, 6, 8]) #将 Py_buffer 拷贝一份 #同时 Py_buffer -> shape 变成了 (2, 3) #但是 Py_buffer -> buf 指向的缓冲区没有拷贝 arr2 = arr1.reshape((2, 3)) #然后在通过索引访问的时候 #可以认为Numpy为其创建了虚拟的索引轴 #由于 arr1 只有一个维度 #那么numpy会为其创建一个虚拟的索引轴 """ arr1 = [3 9 5 7 6 8]: index1: 0 1 2 3 4 5 buf: 3 9 5 7 6 8 """ #arr2 有两个维度，shape 是 (2, 3) #那么Numpy会为其创建两个虚拟的索引轴 """ arr2 = [[3 9 5] [7 6 8]]: index1: 0 0 0 1 1 1 index2: 0 1 2 0 1 2 buf: 3 9 5 7 6 8 """ #缓冲区中索引为 4 的元素被修改 arr2[1, 1] = 666 #但由于 arr1 和 arr2 共享一个缓冲区 #所以 print(arr1[4]) 也会打印 666 print(arr1[4]) # 666
所以缓冲区非常简单，它就是一个一维数组，由 buf 成员指向，而其它的成员则负责描述该如何使用这个缓冲区，可以理解为元信息。正如Numpy的数组，虽然多个数组底层共用一个缓冲区，数据也只有那一份，但是在Numpy的层面却可以表现出不同的维度，究其原因就是元信息不同。
Py_buffer的实现，也是Numpy诞生的一个重要原因。另外，类型对象内部有一个tp_as_buffer成员，它是一个函数指针，在函数内部负责对Py_buffer进行初始化。如果实现了该成员，那么其实例对象便支持缓冲区协议。并且实现了缓冲区协议的对象，不会直接操作缓冲区，而是会借助于 Py_buffer。
相信你现在肯定明白Py_buffer存在的意义了，就是共享内存，实现了缓冲区协议的对象可以直接向彼此暴露对应的缓冲区，比如bytes对象和ndarray对象。
import numpy as np #缓冲区是char类型的一维数组: #{'a', 'b', 'c', 'd', '\0'} b = b"abcd" #直接共享底层的缓冲区 #但是Numpy不知道如何使用这个缓冲区 #所以我们必须显式地指定 dtype #"S1" 表示按照单个字节来进行解析 arr1 = np.frombuffer(b, dtype="S1") print(arr1) # [b'a' b'b' b'c' b'd'] #"S2" 表示按照两个字节来进行解析 arr2 = np.frombuffer(b, dtype="S2") print(arr2) # [b'ab' b'cd'] #那么问题来了，按照三个字节解析是否可行呢？ #答案是不可行，缓冲区的大小不是3的整数倍 #而 "S4" 显然是可以的 arr3 = np.frombuffer(b, dtype="S4") print(arr3) # [b'abcd'] #按照 int8 进行解析 arr4 = np.frombuffer(b, dtype="int8") print(arr4) # [ 97 98 99 100] #按照 int16 进行解析 #显然 97 98 会被解析成一个整数 #99 100 会被解析成一个整数 """ 97 -> 01100001 98 -> 01100010 那么 97 98 组合起来就是 01100010_01100001 99 -> 01100011 100 -> 01100100 那么 99 100 组合起来就是 01100100_01100011 """ print(0b01100010_01100001) # 25185 print(0b01100100_01100011) # 25699 print( np.frombuffer(b, dtype="int16") ) #[25185 25699] #按照int32来解析，显然这4个int8表示一个int32 print( 0b01100100_01100011_01100010_01100001 ) #1684234849 print( np.frombuffer(b, dtype="int32") ) # [1684234849]
怎么样，是不是有点神奇呢？相信你在使用Numpy的时候应该会有更加深刻的认识了，这就是缓冲区协议的威力。哪怕是不同的对象，只要都实现了缓冲区协议，那么彼此之间就可以暴露底层的缓冲区，从而实现共享内存。
所以 np.frombuffer 就是直接根据对象的缓冲区来创建数组，然后它底层的buf成员也指向这个缓冲区。但它不知道该如何解析这个缓冲区，所以我们需要显式地指定dtype来告诉它，相当于告诉它一些元信息。
那么问题来了，我们能不能修改缓冲区呢？
import numpy as np b = b"abcd" arr = np.frombuffer(b, dtype="S1") try: arr[0] = b'A' except ValueError as e: print(e) """ assignment destination is read-only """
答案是不可以的，因为原始的 bytes 对象不可修改，所以缓冲区是只读的。如果想修改的话，可以使用 bytearray。
import numpy as np #可以理解为可变的 bytes 对象 b = bytearray(b"abcd") print(b) # bytearray(b'abcd') #修改 arr arr = np.frombuffer(b, dtype="S1") arr[0] = b'A' #再次打印 print(b) # bytearray(b'Abcd')
小结
到目前为止，我们就解释了什么是缓冲区协议，下面再来总结一下：
- 如果一个类型对象实现了tp_as_buffer，那么它的实例对象便支持缓冲区协议；
- tp_as_buffer 是一个函数指针，指向的函数内部负责初始化Py_buffer；
- Py_buffer 的 buf 成员指向的就是缓冲区，支持缓冲区协议的对象内部的数据都存在缓冲区里面，操作缓冲区数据都是通过 Py_buffer 操作的；
- 实现了缓冲区协议的多个对象可以共享同一个缓冲区，具体做法就是让内部的 buf 成员都指向同一个缓冲区。比如 Numpy 的数组进行切片的时候会得到新数组，而新数组和原数组是共享内存的，原因就是创建新数组的时候只是将 Py_buffer 拷贝了一份，但是 buf 成员指向的缓冲区却没有拷贝；
- 在共享缓冲区的时候，比如np.frombuffer(obj)，会直接调用obj的类型对象的tp_as_buffer成员指向的函数，拿到Py_buffer实例的buf成员指向的缓冲区。但我们说numpy不知道该怎么解析这个缓冲区，所以还需要指定dtype参数；
- 缓冲区存在的最大意义就是共享内存，Numpy的数组在切片的时候，只拷贝 Py_buffer 实例，至于Py_buffer里面的buf成员指向的缓冲区是不会拷贝的。比如数组有100万个元素，这些元素都存在缓冲区中，被Py_buffer里面的buf成员指向，拷贝的时候这100万个元素是不会拷贝的；
- Numpy数组的维度、shape，是借助于Py_buffer中的元信息体现的，至于存储元素的缓冲区，永远是一个一维数组，由 buf 成员指向。只是维度、shape，以及 strides 不同，访问缓冲区元素的方式也不同。但还是那句话，缓冲区本身很单纯，就是一个一维数组。
