Ctypes和NumPy

用ctypes加速計算

Ctypes是Python處理動態鏈接庫的標準擴展模塊，在Windows下使用它可以直接調用C語言編寫的DLL動態鏈接庫。由于對傳遞的參數沒有類型和越界檢查，因此如果編寫的代碼有問題的話，很可能會造成程序崩潰。當將數組數據使用指針傳遞時，出錯誤的風險將更加大。

為了讓程序更加安全，通常會用Python代碼對Ctypes調用進行包裝，在調用Ctypes之前，在Python級別對數據類型和越界進行檢查。這樣做會使得調用接口部分比其它的一些手工編寫的擴展模塊速度要慢，但是如果C語言的代碼段處理相當多的數據的話，接口調用部分的速度損失是可以忽略不計的。

用ctypes調用DLL

為了使用CTypes，你必須依次完成以下步驟：

• 編寫動態連接庫程序
• 載入動態連接庫
• 將Python的對象轉換為ctypes所能識別的參數
• 使用ctypes的參數調用動態連接庫中的函數

下面我們來看看如何用ctypes調用動態鏈接庫。

numpy對ctypes的支持

為了方便動態連接庫的載入，numpy提供了一個便捷函數ctypeslib.load_library。它有兩個參數，第一個參數是庫的文件名，第二個參數是庫所在的路徑。函數返回的是一個ctypes的對象。通過此對象的屬性可以直接到動態連接庫所提供的函數。

例如如果我們有一個庫名為test_sum.dll，其中提供了一個函數mysum :

double mysum(double a[], long n)
{
    double sum = 0;
    int i;
    for(i=0;i<n;i++) sum += a[i];
    return sum;
}

的話，我們可以使用如下語句載入此庫：

>>> from ctypes import *
>>> sum_test = np.ctypeslib.load_library("sum_test", ".")
>>> print sum_test.mysum
<_FuncPtr object at 0x037D7210>

要正確調用sum函數，還必須對其參數類型進行說明，下面的語句描述了sum函數的兩個參數的類型和返回值的類型進行描述：

>>> sum_test.mysum.argtypes = [POINTER(c_double), c_long]
>>> sum_test.mysum.restype = c_double

接下來就可以正常調用sum函數了：

>>> x = np.arange(1, 101, 1.0)
>>> sum_test.mysum(x.ctypes.data_as(POINTER(c_double)), len(x))
5050.0

每次調用sum都需要進行類型轉換時比較麻煩的事情，因此可以編寫一個Python的mysum函數，將C語言的mysum函數包裝起來:

def mysum(x):
    return sum_test.mysum(x.ctypes.data_as(POINTER(c_double)), len(x))

在上面的例子中，test_sum.mysum的參數值使用標準的ctypes類型聲明：用POINTER(c_double)聲明mysum函數的第一個參數是一個指向double的指針；然后調用數組x的x.ctypes.data_as函數將x轉換為一個指向double的指針類型。

由于數組的元素在內存中的存儲可以是不連續的，而且可以是多維數組，因此我們不能指望前面的mysum函數能夠處理所有的情況：

>>> x = np.arange(1,11,1.0)
>>>  mysum(x[::2])
15.0
>>> sum(x[::2])
25.0

由于x[::2]和x共同一塊內存空間，而x[::2]中的元素是不連續的，每個元素之間的間隔為16byptes(2個double的大小)。因此將它傳遞給mysum的話，實際上計算的是x數組中前5項的和：1+2+3+4+5=15，而實際上我們希望的結果是：1+3+5+7+9=25。

為了對傳遞的數組參數進行更加詳細的描述，numpy庫提供了ndpointer函數。ndpointer函數對restype和argtypes中的數組參數進行描述，他有如下4個參數：

• dtype : 數組的元素類型
• ndim : 數組的維數
• shape : 數組的形狀，各個軸的長度
• flags : 數組的標志

例如：

test_sum.mysum.argtypes = [
    np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags="C_CONTIGUOUS"),
    c_long
]

描述了sumfunc函數的參數為一個元素類型為double的、一維的、連續的元素按C語言規定排列的數組。

這時傳遞給mysum函數的第一個參數可以直接是數組，因此無需再編寫一個Python函數對其進行包裝：

>>> sum_test.mysum(x,len(x))
55.0
>>> sum_test.mysum(x[::2],len(x)/2)
ArgumentError: argument 1: <type 'exceptions.TypeError'>:
array must have flags ['C_CONTIGUOUS']

我們注意到如果參數數組不是連續空間的話，mysum函數的調用會拋出異常錯誤，提醒我們其參數需要C語言排列的連續數組。

如果我們希望它能夠處理多維、不連續的數組的話，就需要把數組的shape和strides屬性都傳遞給過去。假設我們想寫一個通用的mysum2函數，它可以對二維數組的所有元素進行求和。下面是C語言的程序：

double mysum2(double a[], int strides[], int shapes[])
{
    double sum = 0;
    int i, j, M, N, S0, S1;
    M = shape[0]; N=shape[1];
    S0 = strides[0] / sizeof(double);
    S1 = strides[1] / sizeof(double);

    for(i=0;i<M;i++){
        for(j=0;j<N;j++){
            sum += a[i*S0 + j*S1];
        }
    }
    return sum;
}

mysum2函數有3個參數，第一個參數a[]指向保存數組數據的內存塊；第二個參數astrides指向保存數組各個軸元素之間的間隔(以byte為單位)；第三個參數dims指向保存數組各個軸長度的數組。

由于strides保存的是以byte為單位的間隔長度，因此需要除以sizeof(double)計算出以double為單位的間隔長度S0和S1。這樣二維數組a中的第i行、第j列的元素可以通過a[i*S0 + j*S1]來存取。下面用ctypes對mysum2函數進行包裝：

sum_test.mysum2.restype = c_double
sum_test.mysum2.argtypes = [
    np.ctypeslib.ndpointer(dtype=np.float64, ndim=2),
    POINTER(c_int),
    POINTER(c_int)
]

def mysum2(x):
    return sum_test.mysum2(x, x.ctypes.strides, x.ctypes.shape)

在mysum2函數中，為了將數組x的strides和shape屬性傳遞給C語言的函數，可以使用x.ctypes中提供的strides和shape屬性。注意不能直接傳遞x.strides和x.shape，因為這些是python的tuple對象，而x.ctypes.shape得到的是ctypes包裝的整數數組：

>>> x = np.zeros((3,4), np.float)
>>> x.ctypes.shape
<numpy.core._internal.c_long_Array_2 object at 0x020B4DF0>
>>> s = x.ctypes.shape
>>> s[0]
3
>>> s[1]
4

可以看出x.ctypes.shape是一個有兩個元素的C語言長整型數組。雖然我們也可以在Python中通過下標讀取其各個元素的值，但是通常它們是作為參數傳遞給C語言函數用的。