楔子
Cython 的兩個優秀的品質就是它的廣度和成熟度，可以編譯所有的 Python 代碼，並且將 C 的速度代入了 Python，並且還能輕松的和 C、C++ 集成。而本篇文章的任務就是完善 Cython 的功能，並介紹 Cython 的陣列特性，比如：對 Numpy 數組的深入支持。
我們已經知道，Cython 可以很好的支持 list、dict、set、tuple 等內置容器，這些容器非常容易使用，可以包含指向任意 Python 對象的變量，並且對 Python 對象的查找、分配、檢索都進行了高度的優化。盡管 list 類型和 dict 類型之間實現存儲和檢索的方式有很大差異，但是從實現角度來講，所有的容器都有一個共同點：它們存儲的都是對 Python 對象的一個引用。如果我們有一個包含100萬個整型的列表，那么在 C 中，每一個元素都是一個 Python *。雖然整型在底層對應 PyLongObject，但是會轉成 PyObject * 存在列表里面，因為 C 中的數組只能存放相同類型的元素（指針也是相同類型的指針），然后用的時候再轉回對應的類型，而 PyObject 里面的 ob_type 成員便是指向對應類型對象的指針。所以 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 數組的高性能是有目共睹的。
有效地使用緩沖區通常是從 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（數組） 是 Python 標准庫中提供的數據結構，它不是 Python 內置的
# 數組不同於列表，因為數組里面存儲的都是連續的整數塊，它的數據存儲要比列表緊湊得多
# 因此一些操作也可以更快的執行
np_arr = np.asarray(arr)
print(np_arr) # [0 1 2 3 4 5]
np_arr[0] = 123
print(arr) # array('i', [123, 1, 2, 3, 4, 5])
Python 提供的數組使用的不是特別多，而 Numpy 中的數組使用的則是非常廣泛，並且支持的操作非常豐富。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])
# 不進行拷貝，此時會共享緩沖區
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 的緩沖區協議本質上是一個結構體，它為多維數組定義了一個非常靈活的接口，我們看一下底層實現，源碼位於 Include\cpython\object.h 中（Python 3.9），小於 Python 3.9 的話則位於 Include\object.h 中。
/* buffer interface */
typedef struct bufferinfo {
void *buf;
PyObject *obj; /* owned reference */
Py_ssize_t len;
Py_ssize_t itemsize; /* This is Py_ssize_t so it can be
pointed to by strides in simple case.*/
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) 的數組，那么底層緩沖區的維度就是 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 來實現跨步訪問，但仍然共享同一個緩沖區，。
import numpy as np
arr1 = np.array(range(10), dtype="int64")
arr2 = arr1[:: 2]
print(arr2) # [0 2 4 6 8]
arr1[0] = 111
print(arr2) # [111 2 4 6 8]
arr1[:] = 0
print(arr2) # [0 0 0 0 0]
既然共享同一個緩沖區，那么改變 arr1 是會影響 arr2 的。
回歸正題，以 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，是因為后面會用到。
下面我們來實現緩沖區協議，由於會涉及到 C 和 Python / C API，所以可能比較枯燥，而且需要你對 C 語言 和 CPython 底層有一定的了解。當然也可以跳過，因為 Cython 允許我們在不了解協議的情況下使用緩沖區，從而在不復制的情況下高效的訪問底層數據、減少開銷，有這一點就足夠了。
實現緩沖區協議（准備工作）
我們說 Numpy 的數組已經實現了緩沖區協議，我們可以直接使用而不需要關注背后的細節，因為緩沖區協議本質上是一個 Python / C API 實體，如果我們想直接操作它那么必須要深入了解 Python / C API。而 Python 則是抽象了 C 級的混亂，給我們提供了一個簡潔、美麗、干凈的 Python 接口。
但有些時候我們不得不親自實現，假設我們有一個 C 對象，需要對其進行包裝然后讓 Python 可用，那么或許你首先會想到用 Cython 來實現（后面會介紹）。當然這是絕對推薦的，因為這是非常正確的做法，Cython 的目的之一就是讓我們能夠更方便地調用現有的 C 庫，此外除了 Cython 你還可以使用 f2py、swig 等等，這些工具的存在仍然可以讓我們不用直接和緩沖區協議打交道。
但話雖如此，可實際了解一下緩沖區協議也是非常有必要的，而想要了解緩沖區協議的最好方式就是手動實現它。一旦了解如何實現緩沖區協議，那么你可以很輕松地將其它語言（比如 Julia）的數據結構，以靈活且優雅的方式暴露給 Python。
下面會涉及到 C 和 Python / C API，需要你對 C 語言 和 CPython 有一定的了解。
// 文件名：my_array.c
#include <stdio.h>
#include <stdlib.h>
// 定義一個一維的數組
typedef struct {
int *arr;
int length;
} MyArray;
// 初始化
void initial_MyArray(MyArray *array, int length) {
// length >= 0
array->length = length;
if (length == 0) {
array->arr = NULL;
} else {
array->arr = (int *) malloc(sizeof(int) * length);
for (int i = 0; i < length; i++) {
array->arr[i] = i;
}
}
}
// 釋放內存
void dealloc_MyArray(MyArray *array) {
if (array->arr != NULL) free(array->arr);
array->arr = NULL;
}
// 將數組進行格式化方便打印
char *change_array_to_string(MyArray *array) {
// 我們模仿 Python 中的列表，將 C 數組格式化成字符串，例如 "[1, 2, 3, 4, 5]"
// 但是最多顯示 10 個
int max = array->length < 10 ? array->length : 10;
char *output = (char *) malloc(100);
// 長度為 0，直接返回 "[]"
if (max == 0) {
sprintf(&output[0], "[]");
return output;
}
// 數據長度不為 0，那么先寫入 "["
int pos = sprintf(&output[0], "[");
// 然后依次寫入數字，由於最后一個數的后面不需要 ", " 所以先設置前 n - 1 個數
for (int i = 0; i < max - 1; i++) {
pos += sprintf(&output[pos], "%d, ", array->arr[i]);
}
// 然后再將最后一個元素設置進去
pos += sprintf(&output[pos], "%d", array->arr[max - 1]);
// 如果數組的長度比 10 大，多余的元素用 ", ..." 顯示
if (array->length > 10) pos += sprintf(&output[pos], ", ...");
// 寫入最后一個 "]"
sprintf(&output[pos], "]");
return output;
}
// 打印數組
void print_MyArray(MyArray *array) {
char *s = change_array_to_string(array);
printf("%s\n", s);
free(s);
}
我們上面定義了一個結構體 MyArray，里面包含了一個數組以及該數組的長度信息，同時還定義了數組的初始化、銷毀、格式化、打印等函數。當然上述這段 C 代碼還是過於簡單，以及錯誤檢測我們也沒有做，像初始化的時候傳遞的 length 小於 0 的話會造成段錯誤，不過這不是我們的重點，目前能用就行。
#include "my_array.c"
int main() {
MyArray array;
initial_MyArray(&array, 0);
print_MyArray(&array);
initial_MyArray(&array, 10);
print_MyArray(&array);
initial_MyArray(&array, 20);
print_MyArray(&array);
dealloc_MyArray(&array);
}
來編譯測試一下：
D:\C\matsuri>gcc -std=c99 main.c -o a.exe
D:\C\matsuri>a.exe
[]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ...]
以上我們就實現了一個類似於 Python 的 range，而且還是 low 版的，但是使用的 C 而不是 Python。下面我們就來將 MyArray 包裝一下，交給 Python 使用，這里會涉及到 Python / C API，需要事先對其有相應的了解。
#include <Python.h>
#include "my_array.c"
// Python 中的對象在 C 中都嵌套了 PyObject
typedef struct {
PyObject_HEAD
MyArray array;
} PyMyArray;
// 初始化函數 __init__
static int
PyMyArray_init(PyMyArray *self, PyObject *args, PyObject *kwargs) {
if (self->array.arr != NULL) {
dealloc_MyArray(&self->array);
}
int length = 0;
static char *kwlist[] = {"length", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|i", kwlist, &length)) {
return -1;
}
// 因為給 Python 調用，所以這里額外對 length 進行一下檢測
if (length < 0) {
PyErr_SetString(PyExc_ValueError, "argument 'length' can not be negative");
return -1;
}
initial_MyArray(&self->array, length);
return 0;
}
// 析構函數
static void
PyMyArray_dealloc(PyMyArray *self) {
dealloc_MyArray(&self->array);
Py_TYPE(self)->tp_free((PyObject *) self);
}
// __str__ 函數
static PyObject *
PyMyArray_str(PyMyArray *self) {
char *s = change_array_to_string(&self->array);
PyObject *ret = PyUnicode_FromString(s);
free(s);
return ret;
}
static PyTypeObject PyMyArrayType = {
PyVarObject_HEAD_INIT(NULL, 0)
"py_my_array.PyMyArray", /* tp_name */
sizeof(PyMyArray), /* tp_basicsize */
0, /* tp_itemsize */
(destructor) PyMyArray_dealloc, /* tp_dealloc */
0, /* tp_print */
0, /* tp_getattr */
0, /* tp_setattr */
0, /* tp_reserved */
(reprfunc) PyMyArray_str, /* tp_repr */
0, /* tp_as_number */
0, /* tp_as_sequence */
0, /* tp_as_mapping */
0, /* tp_hash */
0, /* tp_call */
(reprfunc) PyMyArray_str, /* tp_str */
0, /* tp_getattro */
0, /* tp_setattro */
0, /* tp_as_buffer */
Py_TPFLAGS_DEFAULT, /* tp_flags */
"PyMyArray object", /* tp_doc */
0, /* tp_traverse */
0, /* tp_clear */
0, /* tp_richcompare */
0, /* tp_weaklistoffset */
0, /* tp_iter */
0, /* tp_iternext */
0, /* tp_methods */
0, /* tp_members */
0, /* tp_getset */
0, /* tp_base */
0, /* tp_dict */
0, /* tp_descr_get */
0, /* tp_descr_set */
0, /* tp_dictoffset */
(initproc) PyMyArray_init, /* tp_init */
};
static PyModuleDef py_my_array_module = {
PyModuleDef_HEAD_INIT,
"py_my_array",
"this is a module named py_my_array",
-1,
0,
NULL,
NULL,
NULL,
NULL
};
PyMODINIT_FUNC
PyInit_py_my_array(void) {
PyObject *m;
PyMyArrayType.tp_new = PyType_GenericNew;
if (PyType_Ready(&PyMyArrayType) < 0) return NULL;
m = PyModule_Create(&py_my_array_module);
if (m == NULL) return NULL;
Py_XINCREF(&PyMyArrayType);
PyModule_AddObject(m, "PyMyArray", (PyObject *) &PyMyArrayType);
return m;
}
這里我們算是用 Python / C API 編寫了一個擴展模塊，相信你此刻應該明白為什么要 Cython 了，因為 Python / C API 實在是太痛苦了。
然后我們來編譯成擴展模塊：
from distutils.core import setup, Extension
setup(ext_modules=[Extension("py_my_array", [r"D:\C\matsuri\py_my_array.c"])])
編譯命令和 Cython 是一致的，然后我們來導入測試一下：
import py_my_array
print(py_my_array) # <module 'py_my_array' from ...\\...\\py_my_array.cp38-win_amd64.pyd'>
arr = py_my_array.PyMyArray()
print(arr) # []
arr = py_my_array.PyMyArray(5)
print(arr) # [0, 1, 2, 3, 4]
arr = py_my_array.PyMyArray(20)
print(arr) # [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, ...]
arr = py_my_array.PyMyArray(-1)
"""
arr = py_my_array.PyMyArray(-1)
ValueError: argument 'length' can not be negative
"""
下面我們來給 Numpy 使用，看看會有什么結果：
import numpy as np
import py_my_array
array = py_my_array.PyMyArray(5)
print(array) # [0, 1, 2, 3, 4]
np_arr = np.asarray(array)
# 看到這里相信你已經有所察覺了，因為 Numpy 的數組在打印的時候，元素之間是沒有逗號的
print(np_arr) # [0, 1, 2, 3, 4]
# 而且 shape 為 空元組，很明顯這也是不正常的，應該是 (5,) 才對
print(np_arr.shape) # ()
# 然后我們來修改一下元素
try:
np_arr[0] = 123
except IndexError as e:
# 發現無法修改
print(e) # too many indices for array
# 我們再手動創建一個數組
np_arr = np.array([0, 1, 2, 3, 4])
# 此時一切正常
print(np_arr) # [0 1 2 3 4]
print(np_arr.shape) # (5,)
np_arr[0] = 123
print(np_arr) # [123 1 2 3 4]
我們看到 Numpy 在處理 py_my_array.PyMyArray(5) 的時候，表現非常怪異，原因就在於 PyMyArray 實例化得到的不具備列表的性質，它只是在打印的時候以列表的形式進行打印的。而 Numpy 不知道該如何處理這個對象，於是它創建了一個標量，其值等於這個對象。至於其它對象也是一樣的，舉個例子：
import numpy as np
class A:
def __str__(self):
return "神樂七奈"
arr = np.asarray(A())
print(arr) # 神樂七奈
print(type(arr)) # <class 'numpy.ndarray'>
顯然這不是我們想要的，而解決這一點就是實現緩沖區協議。
使用 C 實現緩沖區協議
再來看一下結構體定義：
/* buffer interface */
typedef struct bufferinfo {
void *buf;
PyObject *obj; /* owned reference */
Py_ssize_t len;
Py_ssize_t itemsize; /* This is Py_ssize_t so it can be
pointed to by strides in simple case.*/
int readonly;
int ndim;
char *format;
Py_ssize_t *shape;
Py_ssize_t *strides;
Py_ssize_t *suboffsets;
void *internal;
} Py_buffer;
buffer interface 所做的事情就是允許你定義一個函數，在函數中構造這樣的一個指向特定數據的結構體。從結構體定義中我們可以看到，數組可以是多維的，並且支持以指定步長訪問。當然對於我們的一維數組而言，我們只會使用一部分簡單的功能。下面我們來創建擴展模塊 py_my_array2，之前的代碼直接搬過來，然后再加上一些緩沖區協議相關的鈎子即可。
#include <Python.h>
#include "my_array.c"
// Python 中的對象在 C 中都嵌套了 PyObject
typedef struct {
PyObject_HEAD
MyArray array;
} PyMyArray;
// 初始化函數 __init__
static int
PyMyArray_init(PyMyArray *self, PyObject *args, PyObject *kwargs) {
if (self->array.arr != NULL) {
dealloc_MyArray(&self->array);
}
int length = 0;
static char *kwlist[] = {"length", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwargs, "|i", kwlist, &length)) {
return -1;
}
// 因為給 Python 調用，所以這里額外對 length 進行一下檢測
if (length < 0) {
PyErr_SetString(PyExc_ValueError, "argument 'length' can not be negative");
return -1;
}
initial_MyArray(&self->array, length);
return 0;
}
// 析構函數
static void
PyMyArray_dealloc(PyMyArray *self) {
dealloc_MyArray(&self->array);
Py_TYPE(self)->tp_free((PyObject *) self);
}
// __str__ 函數
static PyObject *
PyMyArray_str(PyMyArray *self) {
char *s = change_array_to_string(&self->array);
PyObject *ret = PyUnicode_FromString(s);
free(s);
return ret;
}
// 定義 buffer interface 支持的函數
static int
PyMyArray_getbuffer(PyObject *obj, Py_buffer *view, int flags) {
if (view == NULL) {
PyErr_SetString(PyExc_ValueError, "NULL view in getbuffer");
return -1;
}
PyMyArray* self = (PyMyArray*)obj;
view->obj = (PyObject*)self;
view->buf = (void*)self->array.arr;
view->len = self->array.length * sizeof(int);
view->readonly = 0;
view->itemsize = sizeof(int);
view->format = "i";
view->ndim = 1;
view->shape = &self->array.length;
view->strides = &view->itemsize;
view->suboffsets = NULL;
view->internal = NULL;
Py_INCREF(self);
return 0;
}
// 將上面的函數放入到 PyBufferProcs 結構體中
static PyBufferProcs PyMyArray_as_buffer = {
(getbufferproc)PyMyArray_getbuffer,
(releasebufferproc)0
};
static PyTypeObject PyMyArrayType = {
PyVarObject_HEAD_INIT(NULL, 0)
"py_my_array2.PyMyArray", /* tp_name */
sizeof(PyMyArray), /* tp_basicsize */
0, /* tp_itemsize */
(destructor) PyMyArray_dealloc, /* tp_dealloc */
0, /* tp_print */
0, /* tp_getattr */
0, /* tp_setattr */
0, /* tp_reserved */
(reprfunc) PyMyArray_str, /* tp_repr */
0, /* tp_as_number */
0, /* tp_as_sequence */
0, /* tp_as_mapping */
0, /* tp_hash */
0, /* tp_call */
(reprfunc) PyMyArray_str, /* tp_str */
0, /* tp_getattro */
0, /* tp_setattro */
// 指定 tp_as_buffer
&PyMyArray_as_buffer, /* tp_as_buffer */
Py_TPFLAGS_DEFAULT, /* tp_flags */
"PyMyArray object", /* tp_doc */
0, /* tp_traverse */
0, /* tp_clear */
0, /* tp_richcompare */
0, /* tp_weaklistoffset */
0, /* tp_iter */
0, /* tp_iternext */
0, /* tp_methods */
0, /* tp_members */
0, /* tp_getset */
0, /* tp_base */
0, /* tp_dict */
0, /* tp_descr_get */
0, /* tp_descr_set */
0, /* tp_dictoffset */
(initproc) PyMyArray_init, /* tp_init */
};
static PyModuleDef py_my_array_module = {
PyModuleDef_HEAD_INIT,
"py_my_array2",
"this is a module named py_my_array2",
-1,
0,
NULL,
NULL,
NULL,
NULL
};
PyMODINIT_FUNC
PyInit_py_my_array2(void) {
PyObject *m;
PyMyArrayType.tp_new = PyType_GenericNew;
if (PyType_Ready(&PyMyArrayType) < 0) return NULL;
m = PyModule_Create(&py_my_array_module);
if (m == NULL) return NULL;
Py_XINCREF(&PyMyArrayType);
PyModule_AddObject(m, "PyMyArray", (PyObject *) &PyMyArrayType);
return m;
}
然后編譯成擴展，方式不變，編譯之后導入測試一下：
import numpy as np
import py_my_array2
array = py_my_array2.PyMyArray(5)
print(array) # [0, 1, 2, 3, 4]
np_arr = np.asarray(array)
# 我們看到此刻是正常打印的
print(np_arr) # [0 1 2 3 4]
print(np_arr.shape) # (5,)
# 那么兩者是不是共享內存的呢
np_arr[0] = 123
print(array) # [123, 1, 2, 3, 4]
print(np_arr) # [123 1 2 3 4]
顯然此時就萬事大吉了，因為實現了緩沖區協議，Numpy 知道了緩沖區數據，因此會在此基礎之上建一個 view，因此 array 和 np_arr 是共享內存的。
因此核心就在於緩沖區協議的理解，它本質上就是一個結構體，內部的成員描述了緩沖區數據的所有信息。而我們只需要定義一個函數，然后根據數組初始化這些信息即可，最后構建 PyBufferProcs 實例作為 tp_as_buffer 成員的值。
我們目前實現的緩沖區協議非常簡單，但是至少我們理解它是如何工作的了，至於想要實現更復雜的功能，我們會使用 Cython。
使用 Cython 實現緩沖區協議
下面我們來看看如何使用 Cython 實現緩沖區協議，看完之后你會發現使用 Cython 真是一個幸運的選擇：
from cpython cimport Py_buffer
from cpython.mem cimport PyMem_Malloc, PyMem_Free
cdef class Matrix:
cdef Py_ssize_t shape[2] # 數組的形狀
cdef Py_ssize_t strides[2] # 數組的 stride
cdef float *array
def __cinit__(self, row, col):
self.shape[0] = <Py_ssize_t> row
self.shape[1] = <Py_ssize_t> col
self.strides[1] = sizeof(float)
self.strides[0] = self.strides[1] * self.shape[1]
self.array = <float *> PyMem_Malloc(self.shape[0] * self.shape[1] * sizeof(float))
def set_item_by_index(self, int index, float value):
"""留一個接口，用來設置元素"""
if index >= self.shape[0] * self.shape[1] or index < 0:
raise ValueError("索引無效")
self.array[index] = value
def __getbuffer__(self, Py_buffer *buffer, int flags):
"""自定義緩沖區需要實現 __getbuffer__ 方法"""
cdef int i;
for i from 0 <= i < self.shape[0] * self.shape[1]:
self.array[i] = float(i)
buffer.obj = self
buffer.buf = <void *> self.array
buffer.len = self.shape[0] * self.shape[1] * sizeof(float)
buffer.readonly = 0
buffer.itemsize = sizeof(float)
buffer.format = "f"
buffer.ndim = 2
buffer.shape = self.shape
buffer.strides = self.strides
buffer.suboffsets = NULL
def dealloc(self):
if self.array != NULL:
PyMem_Free(<void *> self.array)
在 Cython 中我們只需要實現一個相應的魔法方法即可，真的是非常方法，當然我們為了驗證是否共享內存，專門定義了一個方法。
import pyximport
pyximport.install(language_level=3)
import cython_test
import numpy as np
m = cython_test.Matrix(5, 4)
np_m = np.asarray(m)
print(m) # <cython_test.Matrix object at 0x7f96ba55a3f0>
print(np_m)
"""
[[ 0. 1. 2. 3.]
[ 4. 5. 6. 7.]
[ 8. 9. 10. 11.]
[12. 13. 14. 15.]
[16. 17. 18. 19.]]
"""
# 而且也是共享內存的
m.set_item_by_index(13, 666.666)
print(np_m)
"""
[[ 0. 1. 2. 3. ]
[ 4. 5. 6. 7. ]
[ 8. 9. 10. 11. ]
[ 12. 666.666 14. 15. ]
[ 16. 17. 18. 19. ]]
"""
結果上沒有任何問題。
以上就是緩沖區協議，其實在日常工作中我們不需要直接面對它，但了解一下總是好的。然后我們來介紹一下 Python 的 memoryview，從而引出 Cython 的 memoryview，這樣我們可以在不和緩沖區協議直接打交道的前提下使用緩沖區協議。
memoryview（內存視圖）
一個內置的 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。
b = bytearray(b"my name is satori")
m = memoryview(b)
# 此時 m 和 b 共享內存
print(b) # bytearray(b'my name is satori')
m[0] = ord("M")
print(b) # bytearray(b'My name is satori')
b[1] = ord("Y")
print(chr(m[1])) # Y
當然我們還可以傳入一個 Numpy 的 ndarray，只要是實現了新緩沖區協議的，都可以傳遞到 memoryview 中。
import numpy as np
np_mv = memoryview(np.ones((10, 20, 30)))
# 查看維度
print(np_mv.ndim) # 3
# 查看 shape
print(np_mv.shape) # (10, 20, 30)
# strides 屬性表示某個維度中, 一個元素和下一個元素之間差了多少個字節
print(np_mv.strides) # (4800, 240, 8)
# 查看每個元素的大小
print(np_mv.itemsize) # 8
# 查看整個數組占用的字節數, 等於 itemsize * 元素的個數
print(np_mv.nbytes) # 48000
# 查看數組類型
print(np_mv.format) # d
結構化數據也是支持的，首先我們來創建一個 Numpy 的 dtype，其中 a 和 b 的類型分別是 int8 和 字符串。
import numpy as np
dt = np.dtype([("a", "int8"), ("b", "U")])
print(dt) # [('a', 'i1'), ('b', '<U')]
print(np.empty(5, dtype=dt)) # [(0, '') (0, '') (0, '') (0, '') (0, '')]
structured_mv = memoryview(np.empty(5, dtype=dt))
print(structured_mv.format) # T{b:a:=0w:b:}
"""
這里的 format(格式化字符串) 來自標准庫模塊 struct 的規范, 對於結構化類型來說是相當神秘的
但是我們將 memoryview 的格式化字符串的細節留給官方文檔吧, 我們不需要與它們直接打交道
因為緩沖區和 memoryview 對象可以使用 簡單的標量類型 以及 用戶自定義的結構化類型
"""
那么問題來了，memoryview 對象和緩沖區如何轉換到 Cython 中呢？考慮到 Cython 是專門用來連接 Python 和 C 的，所以它一定非常適合在 C 級別使用 memoryview 對象和緩沖區協議。
類型化 memoryview
Cython 有一個 C 級類型，即：類型化 memoryview，它在概念上和 Python 中的 memoryview 重疊、並且在其基礎上展開。正如其名所示，類型化 memoryview 用於查看（共享）來自緩沖區對象的數據。並且類型化 memoryview 在 C 級別操作，所以它有着最小的 Python 開銷，因此非常高效，而且比直接使用 C 級緩沖區更方便。此外類型化 memoryview 是為了和緩沖區一起工作而被設計的，因此它可以有效支持任何來自緩沖區的對象，從而允許在不復制的情況下共享緩沖區數據。
假設我們想在 Cython 中有效地處理一維數據的緩沖區，而不關心如何在 Python 級別創建數據，我們只是想以一種有效的方式訪問它。
def summer(double[:] numbers):
cdef double res = 0
cdef double number
# memoryview 對象可以像 Python 迭代器一樣進行遍歷
for number in numbers:
res += number
return res
double[:] numbers 聲明了 numbers 是一個類型化 memoryview 對象，而 double 指定了該 memoryview 對象的基本類型，[:] 則表明這是一個一維的 memoryview 對象。
當我們從 Python 中調用 summer 函數時，會傳入一個 Python 對象，該對象會在調用函數的時候隱式的分配給參數 numbers。我們可以提供一個類型化 memoryview 對象，如果提供的不是 memoryview 對象，那么看該對象是否支持緩沖區協議，如果支持緩沖區協議，那么將會提供一個 C 級緩沖區來構建 memoryview；如果不支持該緩沖區協議（沒有提供相應的緩沖區），那么會引發一個 ValueError。
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
# 必須傳遞支持緩存區協議的對象
print(cython_test.summer(np.array([1.2, 2.3, 3.4, 4.5]))) # 11.4
# 可以直接傳入數組，也可以傳入memoryview對象
print(cython_test.summer(memoryview(np.array([1.2, 2.3, 3.4, 4.5])))) # 11.4
# 但是傳遞列表是不行的，因為它不支持緩沖區協議
print(cython_test.summer([1.2, 2.3, 3.4, 4.5]))
"""
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'
"""
不過當我們在編譯類型化 memoryview 對象時，Cython 本質上還是將它當成通用的 Python 迭代器來看待的，因此我們可以做的更好。
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
print(cython_test.summer(np.array([1.2, 2.3, 3.4, 4.5]))) # 11.4
這個版本會有更好的性能：對於百萬元素數組來說大約是 1 毫秒，因為我們用了一個有類型的整數去作為索引，那而這種索引訪問類型化 memoryview 時，Cython 會生成繞過 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
我們通過上下文管理的方式，在其內部對 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 函數的性能便在一個數量級了。我們編譯成擴展模塊測試一下吧：
可以看到，雖然沒有 Numpy 中的 sum 函數那么厲害，但至少在同一水平線上，反正都甩開內置的 sum 函數一條街。
那么到目前為止，我們都了解到了什么呢？首先我們知道如何在 Cython 中聲明一個簡單的類型化 memoryview，以及對它進行索引、訪問內部緩沖區的數據。並且還可以通過使用 boundscheck 和 wraparound 關閉邊界檢查，來生成更加高效的代碼，但前提是我們能確保不會出現索引越界，否則還是不要關閉檢查。因為為了安全，這些都是值得的。
但是還沒完，我們還有更多的細節要講。
聲明類型化 memoryview
當我們聲明一個類型化 memoryview 時，我們可以控制很多的屬性。
元素類型
類型化 memoryview 的元素類型可以任意，在 Cython 中凡是能拿來做變量類型聲明的都可以，因此也可以是 ctypedef 起的別名。
維度
類型化 memoryview 最多可以有7個維度，我們之前聲明一個一維的，使用的是
double[:]這種形式，如果是 3 位，那么寫成double[:, :, :]即可。當然類型不一定是 double，也可以是其它的。
C 和 Fortran 的連續性
通過指定數據打包約束的 C、Fortran 類型化內存視圖是一個非常重要的特例，C 連續和 Fortran 連續都意味着緩沖區在內存中是連續的。但如果是多維度，C 連續的 memoryview 的最后一個維度是連續的，而 Fortran 連續的 memoryview 的第一個維度是連續的。如果可能的話，從性能的角度上來說，將數組聲明為 C 或者 Fortran 連續是有利的，因為這使得 Cython 可以生成更快的代碼。如果不是 C 連續，也不是 Fortran 連續，那么我們稱之為 full strided。
我們通過 Numpy 來對比一下 C 連續和 Fortran 連續的區別。
import numpy as np
arr = np.arange(16)
print(arr.reshape((4 , 4), order="C"))
# 默認是 C 連續, 即 order="C", 最后一個維度是連續的
"""
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]
[12 13 14 15]]
"""
print(arr.reshape((4, 4), order="F"))
# 如果 Fortran, 即 order="F", 那么第一個維度是連續的
"""
[[ 0 4 8 12]
[ 1 5 9 13]
[ 2 6 10 14]
[ 3 7 11 15]]
"""
# 所以 C 連續的數組在轉置之后就會變成 Fortran 連續
直接或間接訪問
直接訪問是默認的，它涵蓋了幾乎所有的情況，它指定對應的維度可以通過索引的方式直接訪問底層數據。如果將某個維度指定為間接訪問，則底層緩沖區將存儲一個指向數組剩余部分的指針，該指針必須在訪問時解引用（因此是間接訪問）。而 Numpy 不支持間接訪問，所以我們不使用這種訪問規范，因此直接訪問是默認的。事實上，官方一般設置為默認的都是最好的。
下面我們來舉例說明：
import numpy as np
# 這是最靈活的聲明方式，可以從任何一個元素為 int 的二維類型化 memoryview 對象中獲取緩沖區
def func(int[:, :] ages):
# 直接打印是一個 memoryview 對象，這里轉成 ndarray, 這里說一句: Numpy 中數組的類型是 <class 'ndarray'>, 為了方便有時叫它 array
print(np.array(ages))
然后我們測試一下：
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
# 這里類型要匹配，C 中的 int 對應 numpy 中的 int32，而 numpy 的整型默認是 int64
# 因為 memoryview 對類型要求很嚴格
arr = np.random.randint(1, 10, (3, 3), dtype="int32")
cython_test.func(arr)
"""
[[7 9 2]
[9 6 9]
[9 9 1]]
"""
# 也可以對 arr 進行切片
cython_test.func(arr[1: 3, 0: 2])
"""
[[9 6]
[9 9]]
"""
當我們對數據進行索引的時候，Cython 會自動生成索引代碼來考慮數據的跨步訪問。但如果我們願意用一些靈活性來換取速度的話，那么 C 連續或者 Fortran 連續的類型化內存視圖可以更有效的建立索引。
聲明一個 C 連續的類型化內存視圖，需要對最后一個維度進行修改。前 n - 1 個維度不變，還是一個冒號，最后一個維度換成兩個冒號並跟一個數字 1。舉個栗子：之前的聲明是 double [:, :, :]，如果想要 C 連續，那么應該改成 double [:, :, :: 1] ，表示最后一個維度具有統一的步長。
import numpy as np
def func(int[:, :: 1] ages):
print(np.array(ages))
然后我們講一個 C 連續的數組賦值給這里的 ages，C 連續在 Numpy 中是所有數組創建函數的默認布局。
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
arr = np.arange(16, dtype="int32").reshape((4, 4))
cython_test.func(arr)
"""
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]
[12 13 14 15]]
"""
但如果賦值一個 Fortran 連續，或者指定步長篩選之后的數組就會引發一個運行的 ValueError。
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
arr = np.arange(16, dtype="int32")
try:
cython_test.func(arr.reshape((4, 4), order="F"))
except ValueError as e:
print(e) # ndarray is not C-contiguous
try:
cython_test.func(arr.reshape((4, 4))[:, :: 2])
except ValueError as e:
print(e) # ndarray is not C-contiguous
如果我們聲明一個 Fortran 連續的數組，那么很明顯，我們需要將第一個維度指定為 :: 1。
import numpy as np
def func(int[::1, :] ages):
print(np.array(ages))
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
arr = np.arange(16, dtype="int32")
try:
cython_test.func(arr.reshape((4, 4)))
except ValueError as e:
print(e) # ndarray is not Fortran contiguous
cython_test.func(arr.reshape((4, 4), order="F"))
"""
[[ 0 4 8 12]
[ 1 5 9 13]
[ 2 6 10 14]
[ 3 7 11 15]]
"""
try:
cython_test.func(arr.reshape((4, 4), order="F")[:: 2])
except ValueError as e:
print(e) # ndarray is not Fortran contiguous
一個多維數組要么是 C 連續，要么是 Fortran 連續，但不能同時是 C 連續和 Fortran 連續。但是一維數組除外，一維數組既可以是 C 連續、也可以是 Fortran 連續。
import numpy as np
def func(int[::1] ages):
print(np.array(ages))
所以我們目前已經介紹了三種類型化內存視圖，分別是：C 連續、Fortran 連續、fully strided。常見的情況下，所有數組都是 C 連續的，這是最常見的內存布局。特別是在需要和外部 C、C++ 庫進行交互的時候，這種布局就顯得尤為重要。並且在很多情況下，如果傳遞了非 C 連續的數組，比如：full strided 或者 Fortran 連續，將會引發一個 ValueError。
但如果你的應該是以 Fortran 為中心的，那么應該將數組聲明為 Fortran 連續，這樣會更好一些。
而 Numpy 提供了 ascontiguousarray 和 asfortranarray 這兩個轉換函數，可以接收一個數組並返回一個 C 連續或者 Fortran 連續的數組。當然，如果數組本身就是 C 連續或者 Fortran 連續，那么就不會轉化了，這樣更有效率。
import numpy as np
arr = np.arange(16).reshape((4, 4))
print(arr.flags["C_CONTIGUOUS"]) # True
print(arr.flags["F_CONTIGUOUS"]) # False
arr = np.asfortranarray(arr)
print(arr.flags["C_CONTIGUOUS"]) # False
print(arr.flags["F_CONTIGUOUS"]) # True
# 這兩個方法底層都調用了 np.array, 所以下面這種做法也是可以的
arr2 = np.arange(16).reshape((4, 4))
print(arr2.flags["F_CONTIGUOUS"]) # False
# 這里面有一個 copy 參數, 表示是否拷貝, 默認為 True, 這樣得到的新數組和原數組之前沒有任何關系
# 但這里我們顯然無需拷貝, 因此指定 copy=False 會更有效率
arr2 = np.array(arr2, order="F", copy=False)
print(arr2.flags["F_CONTIGUOUS"]) # True
# 如果要將一個數組改成 C 連續或者 Fortran 連續, 推薦上面兩種做法
# 另外, 我們說一維數組既是 C 連續又是 Fortran 連續
arr3 = np.arange(16)
print(arr3.flags["C_CONTIGUOUS"]) # True
print(arr3.flags["F_CONTIGUOUS"]) # True
但是這並不代表我們在 Cython 中聲明 memoryview 對象時，一定非要聲明成 C 連續或者 Fortran 連續，當你接收的數組的連續性不確定時，應該采用 full strided 類型，也就是聲明的時候不指定 :: 1。因為一旦指定連續，不管是 C 連續、還是 Fortran 連續，那么你的數組必須要滿足，否則就會報出我們之前出現的錯誤：ndarray is not C-contiguous 或者 ndarray is not Fortran contiguous。這個時候你就需要先進行轉化，通過手動創建一個 C 連續或者 Fortran 連續的數組（或者說 memoryview）的副本，說白了就是將那些指定步長訪問的數組對應的元素拷貝一份，建立一個新的連續數組，但是這會帶來額外的開銷，甚至會超過連續訪問帶來的性能收益，舉個栗子：
import numpy as np
arr = np.arange(16).reshape((4, 4))
# 默認情況下 copy=True，創建數組時會將元素拷貝一份，此時創建的數組默認是 C 連續
# 但如果是在不拷貝的情況下根據已有的數組創建數組，那么元素不會拷貝，而是共享緩沖區
# 此時創建的數組的連續性和已有的數組保持一致，因此，由於 arr[:: 2] 不是連續的，所以 arr1 也不是
arr1 = np.array(arr[:: 2], copy=False)
arr1[0, 0] = 111
# arr 和 arr1 共享緩沖區，修改 arr1 會改變 arr
print(arr[0, 0]) # 111
# 但 arr1 不是 C 連續，因為它指定了步長，實現了跨步訪問，所以它就不再具備連續性
print(arr1.flags["C_CONTIGUOUS"]) # False
# 還是相同的方式，但是我創建的時候強行指定讓 arr2 連續
arr2 = np.array(arr[:: 2], copy=False, order="C")
print(arr2.flags["C_CONTIGUOUS"]) # True
arr2[1, 1] = 222
# 但是我們看到 arr 並沒有被改變
# 原因就在於將一個不是連續的數組變成連續的數組，會將不是連續的數組中對應的元素拷貝一份
# 以此來構建一個連續的數組，所以此時指定 copy=False 是沒有用的
print(arr[1, 1]) # 5
# 所以在不指定連續性的情況下，指定 copy=False 則是和原來的數組共享緩沖區
# 但指定了連續性之后，只能拷貝一份，會有額外的性能開銷
對於類型化 memoryview，我們傳遞 None 也是合法的，因此需要有一步檢測，或者使用 not None 子句聲明。
混合類型
還記得我們之前提到的混合類型嗎？假設我希望某一個參數既可以是接收 list 對象，也可以接收 dict 對象，那么可以這么做。
cdef fused list_dict:
list
dict
cpdef func(list_dict var):
return var
而類型化內存的視圖的類型也可以是混合類型，這樣可以保證更強的泛化能力和靈活性。但是很明顯，所謂的混合類型無非就是創建了多個版本。
cimport cython
cpdef cython.floating generic_summer(cython.floating[:] m):
cdef cython.floating f, s = 0.0
for f in m:
s += f
return s
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
print(
cython_test.generic_summer(np.array([1, 2, 3], dtype="float64"))
) # 6.0
print(
cython_test.generic_summer(np.array([1, 2, 3], dtype="float32"))
) # 6.0
類型化內存視圖對元素類型的要求是很嚴格的，但是我們通過混合類型的方式可以同時接受 float32 和 float64，也就是 C 中的 float 和double。
但是類型化 memoryview 對混合類型目前支持的其實不是特別友好，個人不太建議使用。
使用類型化 memoryview
一旦聲明了類型化 memoryview，就必須給它分配一個支持緩沖區的對象，然后兩者共享底層緩沖區。
那么問題來了，類型化 memoryview 都支持哪些操作呢？
首先我們可以像 Numpy 一樣，對類型化 memoryview 進行訪問和修改。
import numpy as np
cpdef array(int[:, :] numbers):
print("----------")
print(np.array(numbers))
numbers[0, 0] = 3
print("----------")
print(np.array(numbers))
print("----------")
print(np.array(numbers[1: 3, : 2]))
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
arr = np.random.randint(1, 10, (3, 3))
cython_test.array(arr)
"""
----------
[[5 8 1]
[7 9 9]
[4 6 2]]
----------
[[3 8 1]
[7 9 9]
[4 6 2]]
----------
[[7 9]
[4 6]]
"""
正如我們之前說的，類型化內存視圖可以建立高效的索引，特別是當我們通過 boundscheck、wraparound 關閉檢查的時候。
from cython cimport boundscheck, wraparound
cpdef summer(int[:, :] numbers):
cdef int N, M, i, j
cdef long s=0
# 這里類型化 memoryview 的 shape 是一個含有 8 個元素的元組
# 但我們這里只有兩個維度, 所以截取前兩位, 至於后面的元素都是 0
N, M = numbers.shape[: 2]
with boundscheck(False), wraparound(False):
for i in range(N):
for j in range(M):
s += numbers[i, j]
return s
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
arr = np.random.randint(1, 10, (300, 300))
print(np.sum(arr)) # 450832
print(cython_test.summer(arr)) # 450832
另外類型化 memoryview 和 Numpy 中的 array 一樣，也支持 ... 這個語法糖，表示某個維度上、或者整體的全部篩選。
import numpy as np
cdef int[:, :] m = np.zeros((2, 2), dtype="int32")
print("*" * 10)
# 直接打印的話是一個 memoryview 對象, 需要再轉成 array 進行打印
print(np.array(m))
# 通過 ... 表示全局修改
# 因此 m[...] 等價於 m[:]
m[...] = 123
print("*" * 10)
print(np.array(m))
# 在某一個維度上使用 ..., 可以實現某個維度上的全局修改
# 等價於 m[0, :] = 456
m[0, ...] = 456
print("*" * 10)
print(np.array(m))
import pyximport
pyximport.install(language_level=3)
import cython_test
"""
**********
[[0 0]
[0 0]]
**********
[[123 123]
[123 123]]
**********
[[456 456]
[123 123]]
"""
因此在用法上，類型化 memoryview 和 Numpy 的 array 是一致的，當然我們也可以指定步長等等。注意：類型化 memoryview 的右邊只能跟標量。
但是類型化 memoryview 和 Numpy 的 array 還是有些差別的，類型化 memoryview 不支持通用函數，所以在賦值的時候除了簡單的標量賦值之外，無法進行廣播操作。
import numpy as np
arr = np.arange(9).reshape((3, 3))
print(arr)
"""
[[0 1 2]
[3 4 5]
[6 7 8]]
"""
# 會將所有元素都賦值為 123, 因為是標量賦值, 所以類型化 memoryview 是支持的
arr[:] = 123
print(arr)
"""
[[123 123 123]
[123 123 123]
[123 123 123]]
"""
# 這里就涉及到了廣播, 因為 (3, 3) 和 (3,) 兩個維度明顯不一致
# 所以會將 arr 的每一行都替換成 [11 22 33], 但這個是 Numpy 的 array 的功能, memoryview 是不支持的, 因為它在廣播的時候右邊只能跟標量
arr[:] = np.array([11, 22, 33])
print(arr)
"""
[[11 22 33]
[11 22 33]
[11 22 33]]
"""
可能有人覺得這也太不方便了吧，但辦法總比困難多，我們可以根據類型化 memoryview 拿到對應的 array，對這個 array 進行操作不就行了。那么這么做的效率會不會低呢？答案是不會的，因為類型化 memoryview 和 array 之間是共享內存的，這么做不會有任何的性能損失。正如 torch 里面的 tensor 一樣，它和 Numpy 的 array 之間也是共享內存的，由於 Numpy 的 API 用起來非常方便，也已經習慣了，加上個人也懶得使用 tensor 的一些操作，所以我都會先將 tensor 轉成 array，對 array 操作之后再轉回 tensor。雖然多了兩次轉化，但還是那句話，它們是共享內存的，所以完全沒問題。
from cython cimport boundscheck, wraparound
import numpy as np
cdef int[:, :] m = np.arange(9).reshape((3, 3))
# 這里一定要指定 copy=False, 否則在拿到緩沖區里面的內容時, 還會拷貝一份新的出來
# 這么做的話, 就會有性能損失了, 因為兩者本來就是共享內存的, 獲取出來之后直接操作就可以了, 為什么要再創建一個副本呢
# 另外, 如果拷貝了一個副本, 你再操作完畢之后還要賦值回去
np.array(m, copy=False)[:] = [1, 2, 3]
# 以上我們就實現了修改, 這里我們再來打印一下看看, 這里打印就創建一個新的吧, 看看數據有沒有修改
print(np.array(m))
import pyximport
pyximport.install(language_level=3)
import cython_test
"""
[[1 2 3]
[1 2 3]
[1 2 3]]
"""
因此我們可以把類型化 memoryview 看成是非常靈活 Cython 空間，可以有效地共享、索引定位、以及修改同質數據。它具有很多 Numpy array 的特性，特別是通過索引定位數據。至於那些沒有特性，也很容易被兩者之間轉換的高效性所掩蓋。
但是類型化 memoryview 實際上是超越了緩沖區協議，所以它還有額外的特性，那就是它還可以對 C 一級的數組進行 view。
view C 級數組
C 級數組可以是在堆上動態分配的，也可以是在棧上分配的。如果要 view C 一級的數組，那么該數組賦值給 memoryview 即可。如果數組的大小是固定的（或完整的），賦值的右邊直接寫數組名即可，因為 Cython 有足夠的信息來跟蹤這個 C 數組的大小。
cdef int a[3][5][7]
cdef int[:, :, :: 1] m = a
m[...] = 123
# 這里我們將 memoryview 聲明是 C 連續的, 因為 C 里面的數組當然是 C 連續的
# 然后將其賦值為 123
import numpy as np
# 轉成 Numpy 中的數組
arr = np.array(m, copy=False)
print("*" * 10)
print(np.sum(arr), 3 * 5 * 7 * 123)
print("*" * 10)
print(arr[:, 1: 3])
import pyximport
pyximport.install(language_level=3)
import cython_test
"""
**********
12915 12915
**********
[[[123 123 123 123 123 123 123]
[123 123 123 123 123 123 123]]
[[123 123 123 123 123 123 123]
[123 123 123 123 123 123 123]]
[[123 123 123 123 123 123 123]
[123 123 123 123 123 123 123]]]
"""
以上 C 數組是在棧區分配的，我們也可以改為堆區分配。
from libc.stdlib cimport malloc, free
cdef int *a = <int *>malloc(3 * 5 * 7 * sizeof(int))
# 但是很明顯, 這樣做的話, 形狀的信息就丟失了, Cython 只知道這個一個 int *
# 所以直接將 a 賦值過去的話, 編譯時會出現 "Cannot convert long * to memoryviewslice"
# 因此我們在賦值給 memoryview 的時候, 必須給 Cython 提供更多的信息
# 而 <int[:3, :5, :7]> a 會告訴 Cython 這是一個三維數組, 維度分別是 3 5 7; 當然我們這里變成 7 5 3 也是可以的, 因為形狀是由我們決定的
cdef int[:, :, :: 1] m = <int[:3, :5, :7]> a
m[...] = 123
import numpy as np
arr = np.array(m, copy=False)
print("*" * 10)
print(np.sum(arr), 3 * 5 * 7 * 123)
print("*" * 10)
print(arr[:, 1: 3])
# 最后將 a 給釋放掉
free(a)
導入這個模塊之后，打印的結果和之前是一樣的。
在 C 級別，光靠一個頭指針是沒有辦法確定動態分配的 C 數組的長度，而這一點是需要由我們來確定的。因此將一個 C 數組轉成類型化 memoryview 時，如果數據不正確，那么可能會導致緩沖區溢出、段錯誤，或者數據損壞等等。
到此我們就算介紹完了類型化memoryview 的特性，並展示了它如何在支持新緩沖區協議的 Python 對象和 C 級數組中使用它們。如果 Cython 函數中有一個類型化 memoryview 參數，那么可以通過傳遞一個支持新緩沖區協議的 Python 對象或者 C 數組作為參數進行調用。
另外，當在函數中想返回一個類型化 memoryview 時，Cython 會直接將緩沖區內容（沒有拷貝）並轉化成 Python 的 memoryview，比如在函數中返回一個由 C 數組構建的類型化 memoryview，假設叫 m。如果這個 C 數組是在堆區通過 malloc 動態分配的，那么返回 m 沒有任何問題；但如果它是在棧區分配的，在函數結束后就會被銷毀，而我們說 Cython 又不會對緩沖區內容進行拷貝，因此此時就會出現錯誤。
不過即便是 C 數組在堆區申請，也是存在問題的。那就是當我們不再使用 memoryview 的時候，誰來釋放這個在堆上申請的數組呢？如何正確地管理它的內存呢？不過在探討這個問題之前，我們需要先來說一下另一種 Numpy。
另一種 Numpy
什么叫做另一種 Numpy，因為我們在 Cython 中除了 import numpy 之外，還可以 cimport numpy。
在類型化 memoryview 之前，Cython 可以使用不同的語法來很好地處理 Numpy 數組，這便是 原始緩沖區語法。盡管它已經被類型化 memoryview 所取代，但我們依舊可以正常使用它。
# 這里一定要使用 cimport
cimport numpy as np # 或者 from numpy cimport ndarray
# 如果是 import numpy as np, 那么不好意思, np.ndarray 是無法作為參數類型和返回值的
# 編譯時會報錯: 'np' is not a cimported module
# 如果是 from numpy import ndarray, 編譯時同樣報錯: 'ndarray' is not a type identifier
cpdef np.ndarray func(np.ndarray array):
return array
但是注意：此時我們不能夠直接導入，因為它依賴 numpy 的一個頭文件，所以我們需要通過編譯的方式。
from pathlib import Path
import numpy as np
from distutils.core import Extension, setup
from Cython.Build import cythonize
ext = Extension("cython_test",
["cython_test.pyx"],
# 因為我們 cimport numpy 時會使用 Numpy 提供的一個頭文件, 叫做 arrayobject.h
# 但是很明顯, 我們並沒有指定它的位置, 因此需要通過 include_dirs 告訴 cython 編譯器去哪里找這個頭文件
# 如果沒有這個參數, 那么編譯時會報錯：fatal error: numpy/arrayobject.h: No such file or directory
include_dirs=[str(Path(np.__file__).parent / "core" / "include")])
# 當然 numpy 也為我們封裝了一個函數, 直接通過 np.get_include() 即可獲取改路徑
# 對於我當前的環境來說就是 C:\python38\lib\site-packages\numpy\core\include
setup(
ext_modules=cythonize(ext, language_level=3),
)
然后導入測試一下：
import numpy as np
import cython_test
print(cython_test.func(np.array([[1, 2], [3, 4]])))
"""
[[1 2]
[3 4]]
"""
print(cython_test.func(np.array([["xx", None], [(1, 2), {1, 2}]])))
"""
[['xx' None]
[(1, 2) {1, 2}]]
"""
測試是沒有問題的，並且接收的就是 array，返回的也是 array。並且我們看到，這對 array 的類型沒有任何限制，但如果我們希望限制的 array 的類型、甚至是維度，這個時候這怎么做呢？
Cython為numpy提供了專門的方法：比如希望參數是一個接收類型為 int64、維度為 2 的數組，就可以使用 ndarray[long, ndim=2] 這種方式，我們演示一下。
cimport numpy as np
# 類型需要統一
# long 對應 np.int64, int 對應 np.int32, short 對應 np.int16, char 對應 int8
# unsigned long 對應 np.uint64, 其它同理
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)
# 除了作為函數參數和返回值之外, 還可以用來聲明普通的靜態變量
# 比如: cdef np.ndarray[double, ndim=2] arr
編譯之后並且導入進行測試：
import numpy as np
import cython_test
# 這里我們傳遞的時候, 參數和維度一定要匹配, 默認是 int64
cython_test.func1(np.array([[1, 2], [3, 4]]))
"""
[[1 2]
[3 4]]
"""
try:
cython_test.func1(np.array([1, 2, 3, 4]))
except ValueError as e:
print(e) # Buffer has wrong number of dimensions (expected 2, got 1)
try:
cython_test.func2(np.array([1, 2, 3, 4]))
except ValueError as e:
print(e) # Buffer dtype mismatch, expected 'double' but got 'long'
cython_test.func2(np.array([1, 2, 3, 4], dtype="float64"))
"""
[1. 2. 3. 4.]
"""
cython_test.func3(np.array(["a", "b", object], dtype="O"))
"""
['a' 'b' <class 'object'>]
"""
但以上是原始緩沖區語法，現在我們更推薦類型化 memoryview，雖然它比 Numpy 中的 array 少了許多功能，但我們說這兩者之間是可以高效轉換的。並且如果是通過 np 來調用的話，那么兩者是等價的。舉個栗子：
import numpy as np
def func(int[:, :: 1] m):
# m 本身沒有 sum 方法, 但是我們可以將它傳遞給 np.sum
print(np.sum(m))
print(np.sum(m, axis=0))
print(np.sum(m, axis=1))
import pyximport
pyximport.install(language_level=3)
import numpy as np
import cython_test
cython_test.func(np.array([[1, 2], [3, 4]], dtype="int32"))
"""
10
[4 6]
[3 7]
"""
因此這種聲明方式是更加推薦的，而且也更加清晰和簡潔，以及我們可以通過 pyximport 進行導入了。之前的方式由於依賴一個頭文件，你必須要告訴 cython 編譯器頭文件去哪里找，而直接導入、自動編譯的話會找不到頭文件。但是現在不需要了，因為我們根本沒有使用 cimport numpy。
但是以上這些說實話都不能算是優點，所以肯定還有其它的優點，那么都有哪些呢？
1. 類型化 memoryview 支持的對象種類非常多，只要是支持緩沖區的對象均可，比如：Numpy array、Python memoryview 對象、array.array 對象、支持緩沖區協議任意類型的對象，並且它也適用於 C 數組。所以它比原始緩沖區語法更加通用，原始緩沖區語法只適用於 Numpy array。
2. 類型化 memoryview 有着更多的選擇來控制數據到底是 C 連續還是 Fortran 連續，以及直接或者間接訪問。並且一些選項可以按照維度逐個控制，而 Numpy 的原始緩沖區語法不提供這種級別的控制。
3. 在任何情況下，類型化 memoryview 都有着超越原始緩沖區語法的性能，這一點才是我們最關注的，至於語法的間接程度其實沒有太大的意義，個人這么認為。
包裝 C、C++ 數組
回到我們之前的問題，當我們 C 的數組在堆上分配，那么返回之后要如何釋放堆區申請的內存呢？我們舉個栗子：
// heap_malloc.h
float *make_matrix_c(int nrows, int ncols);
// heap_malloc.c
#include <stdlib.h>
float *make_matrix_c(int nrows, int ncols) {
float *matrix = (float *) malloc(nrows * ncols * sizeof(float));
return matrix;
}
以上是返回一個在堆上分配的 C 數組：
import numpy as np
cdef extern from "heap_malloc.h":
float *make_matrix_c(int nrows, int ncols)
def make_matrix(int nrows, int ncols):
cdef float[:, :: 1] m = <float[:nrows, :ncols]> make_matrix_c(nrows, ncols)
# 因為元素都未初始化, 所以里面的值都是不確定的, 雖然這無傷大雅, 但是更優雅的處理方式是將值都初始化為零值
m[...] = 0.0
# 如果直接返回 m, 那么 Python 在接收到之后會打印 <MemoryView of 'ndarray' object>
# 所以這里我們轉成 array 之后返回, 當然這里返回 m、然后在 Python 中接收到之后再將 m 轉成 array 也是可以的
# 只不過畢竟涉及到類型化 memoryview, 它是 Cython 的東西, 所以這一步還是在 Cython 中直接完成會比較好, 返回給 Python 的就是一個普通的 array
return np.array(m, copy=False)
顯然這個例子已經無需解釋了，然后我們直接編譯：
from distutils.core import Extension, setup
from Cython.Build import cythonize
ext = Extension("cython_test",
["cython_test.pyx", "heap_malloc.c"])
setup(
ext_modules=cythonize(ext, language_level=3),
)
編譯完成之后將文件移到當前目錄，導入測試：
import cython_test
arr1 = cython_test.make_matrix(3, 4)
print(arr1.__class__) # <class 'numpy.ndarray'>
print(arr1)
"""
[[0. 0. 0. 0.]
[0. 0. 0. 0.]
[0. 0. 0. 0.]]
"""
arr2 = cython_test.make_matrix(2, 5)
print(arr2)
"""
[[0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0.]]
"""
整體沒有任何問題，很 happy，但是 arr1 和 arr2 對應的 C 數組都是在堆區，那這兩個堆區的數組咋辦？所以目前這個 make_matrix 是存在致命缺陷的，因為它有內存泄露的風險，而這個風險對於任何一個程序而言都是具有致命性的。
使用 Cython 正確（自動）地管理 C 數組
首先我們是需要對內存負責的，但很明顯當和 C 共享數組的時候，合適的解決內存問題就變成了一件很棘手的問題，因為 C 語言沒有自動管理內存的特性。通常在這種情況下，最干凈利索的做法就是復制數據，來澄清各自對數據的所有權。比如：在返回 array 的時候，我們不加 copy=False，這樣就會創建一個新的副本，所以你的是你的，我的是我的，兩者之間沒有關系。
from libc.stdlib cimport free
import numpy as np
cdef extern from "heap_malloc.h":
float *make_matrix_c(int nrows, int ncols)
def make_matrix(int nrows, int ncols):
cdef float[:, :: 1] m = <float[:nrows, :ncols]> make_matrix_c(nrows, ncols)
m[...] = 0.0
matrix = np.array(m)
# 這里直接就可以釋放了, 因為此時的 matrix 里面的元素是獨立的
free(m)
return matrix
但很明顯，如果數據量非常大，我們這么做是不是很影響效率呢？所以它雖然是一個解決問題的辦法，但不是最好的辦法。
應該有一種機制，在 Cython 中返回數組之前應該將這個數組和某個函數關聯起來，一旦當數組被回收，那么要觸發相應的函數，而在這個函數里面執行釋放堆區內存的邏輯。那么我們看看這是如何實現的。
首先 Numpy 中的數組是在底層是一個 PyArrayObject，而 Numpy / C API 在此基礎之上定義了一個 base 屬性，而這個 base 屬性就是為了這個目的而存在的。如果你正在使用 C API 來構建一個數組，並在堆區手動申請內存，那么你應該使用 PyArray_SetBaseObject 函數將 base 屬性設置在內存申請在堆區的數組中。但是為了達到同樣的目的，我們會使用 Cython 提供的函數，而不是使用 PyArray_SetBaseObject。
下面我們修改一下 pyx 文件，C 的頭文件和源文件不需要做改動。
import numpy as np
from libc.stdlib cimport free
cimport numpy as c_np # 用於訪問 Numpy 的 C API
cdef extern from "heap_malloc.h":
float *make_matrix_c(int nrows, int ncols)
cdef class _finalizer:
# 定義一個類, 這個 void * 顯然是申請在堆區的數組
cdef void *_data
def __dealloc__(self):
# 析構函數, 這里再打印一句話
print("數組被回收了")
if self._data != NULL:
free(self._data)
cdef void *set_base(c_np.ndarray arr, void *c_arr):
# 函數接收兩個參數: ndarray, 堆區申請的數組轉成的 void * 指針
# 實例化一個類
cdef _finalizer f = _finalizer()
# 將數組轉成 void *, 交給 f._data
f._data = c_arr
# 重點來了, 這里便是底層對應的 PyArray_SetBaseObject, 但這里我們可以很輕松的進行調用
# 因為 Cython 幫我們做好了這一切, 調用這個函數將 ndarray 和 f 關聯起來, 而 f 里面存儲了指向堆區的指針
# 當 arr 被銷毀時, 會觸發 f 內部的 __dealloc__ 方法, 在這個方法將指針指向的內存釋放掉
c_np.set_array_base(arr, f)
def make_matrix(int nrows, int ncols):
# 獲取堆區申請的指針
cdef float *c_arr = make_matrix_c(nrows, ncols)
cdef float[:, :: 1] m = <float[:nrows, :ncols]> c_arr
m[...] = 0.0
# 這里通過 cdef c_np.ndarray 的方式靜態聲明
cdef c_np.ndarray matrix = np.array(m, copy=False)
# 調用 set_base, 傳入 ndarray 和 指向堆區的指針
set_base(matrix, <void *> c_arr)
return matrix
然后我們來重新編譯一下：
import numpy as np
from distutils.core import Extension, setup
from Cython.Build import cythonize
ext = Extension("cython_test",
["cython_test.pyx", "heap_malloc.c"],
# 在當前路徑和 np.get_include() 中尋找頭文件
include_dirs=[".", np.get_include()])
setup(
ext_modules=cythonize(ext, language_level=3),
)
說實話個人覺得還是比較麻煩的，如果 C 返回的數組不大的話，還是直接拷貝一份吧，這樣絕對是最方便的選擇。話說只是一個數組而已，如果你不是對性能要求到極致的話，那么直接通過拷貝的方式絕對是最方便、最穩妥的選擇。
總結
到目前為止，我們介紹了 Python 中的各種可以轉成類型化 memoryview 的對象，但是 Numpy array 絕對是最具普遍性、靈活性以及表現力。而且除了 Python 對象之外，類型化 memoryview 還可以使用 C 級數組，不管是棧上分配，還是堆上分配。
而核心就在於 Cython 的類型化 memoryview，它提供了一個一致的抽象，這個抽象適用於所有支持緩沖區協議的對象，並為我們提供了高效的 C 級訪問緩沖區，並且它們是共享內存的。
