『Python CoolBook』Cython_高效數組操做

數組運算加速是相當科學計算重要的領域，本節咱們以一個簡單函數爲例，使用C語言爲python數組加速。

1、Cython

本函數爲一維數組修剪最大最小值

version1

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef clip(double[:] a, double min, double max, double[:] out):
    '''
    Clip the values in a to be between min and max. Result in out
    '''
    if min > max:
        raise ValueError("min must be <= max")
    if a.shape[0] != out.shape[0]:
        raise ValueError("input and output arrays must be the same size")
    for i in range(a.shape[0]):
        if a[i] < min:
            out[i] = min
        elif a[i] > max:
            out[i] = max
        else:
            out[i] = a[i]

利用Cython類型的內存視圖，極大的簡化了數組的操做。

• cpdef clip() 聲明瞭 clip() 同時爲C級別函數以及Python級別函數。 在Cython中，這個是很重要的，由於它表示此函數調用要比其餘Cython函數更加高效 （好比你想在另一個不一樣的Cython函數中調用clip()）。 
•  類型參數 double[:] a 和 double[:] out 聲明這些參數爲一維的雙精度數組。 做爲輸入，它們會訪問任何實現了內存視圖接口的數組對象，這個在PEP 3118有詳細定義。 包括了NumPy中的數組和內置的array庫。
• clip() 定義以前的兩個裝飾器能夠優化下性能：
  • @cython.boundscheck(False) 省去了全部的數組越界檢查， 當你知道下標訪問不會越界的時候可使用它
  • @cython.wraparound(False) 消除了相對數組尾部的負數下標的處理（相似Python列表）

version2_條件表達式

任什麼時候候處理數組時，研究並改善底層算法一樣能夠極大的提示性能

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef clip(double[:] a, double min, double max, double[:] out):
    if min > max:
        raise ValueError("min must be <= max")
    if a.shape[0] != out.shape[0]:
        raise ValueError("input and output arrays must be the same size")
    for i in range(a.shape[0]):
        out[i] = (a[i] if a[i] < max else max) if a[i] > min else min

 version3_釋放GIL

釋放GIL，這樣多個線程能並行運行，要這樣作的話，須要修改代碼，使用 with nogil:

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef clip(double[:] a, double min, double max, double[:] out):
    if min > max:
        raise ValueError("min must be <= max")
    if a.shape[0] != out.shape[0]:
        raise ValueError("input and output arrays must be the same size")
    with nogil:
        for i in range(a.shape[0]):
            out[i] = (a[i] if a[i] < max else max) if a[i] > min else min

編寫setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_modules = [
    Extension('sample',
              ['sample.pyx'])
]

setup(
  name = 'Sample app',
  cmdclass = {'build_ext': build_ext},
  ext_modules = ext_modules
)

使用 python3 setup.py build_ext --inplace 來構建它

效率測試示意以下，

>>> import sample >>> import numpy >>> b = numpy.random.uniform(-10,10,size=1000000) >>> c = numpy.zeros_like(b)
>>> import timeit
>>> timeit.timeit('numpy.clip(b,-5,5,c)','from __main__ import b,c,numpy',number=1000)
>>> timeit.timeit('sample.clip(b,-5,5,c)','from __main__ import b,c,sample', ... number=1000)

其中使用numpy本身的clip對比試驗，

2.6287411409430206  # numpy
2.8034782900940627 # v1
2.7247575907967985 # v2
2.6071253868285567 # v3

版本三近似於numpy的實現效果，其餘版本差一些(每次試驗結果都會略有差別，這裏只是粗略的比較一下)。

二維數組處理版本參考：

@cython.boundscheck(False)
@cython.wraparound(False)
cpdef clip2d(double[:,:] a, double min, double max, double[:,:] out):
    if min > max:
        raise ValueError("min must be <= max")
    for n in range(a.ndim):
        if a.shape[n] != out.shape[n]:
            raise TypeError("a and out have different shapes")
    for i in range(a.shape[0]):
        for j in range(a.shape[1]):
            if a[i,j] < min:
                out[i,j] = min
            elif a[i,j] > max:
                out[i,j] = max
            else:
                out[i,j] = a[i,j]

2、本身寫接口

點擊進入項目

sample.c中添加

/* n:longth of array */
void clip(double *a, int n, double min, double max, double *out) {
  double x;
  for (; n >= 0; n--, a++, out++) {
    x = *a;

    *out = x > max ? max : (x < min ? min : x);
  }
}

pysample.c中添加

// void clip(double *a, int n, double min, double max, double *out);
static PyObject *py_clip(PyObject *self, PyObject *args){
  PyObject *a, *out;
  int min, max;
  if(!PyArg_ParseTuple(args, "OiiO", &a, &min, &max, &out)){  //py數組對象暫記
    return NULL;
  }

  // printf("%i, %i\n", min, max);
  Py_buffer view_a, view_out;  //py數組對象接收對象
  if (PyObject_GetBuffer(a, &view_a,
      PyBUF_ANY_CONTIGUOUS | PyBUF_FORMAT) == -1) {
    return NULL;
  }
  if (PyObject_GetBuffer(out, &view_out,
      PyBUF_ANY_CONTIGUOUS | PyBUF_FORMAT) == -1) {
    return NULL;
  }
  clip(view_a.buf, view_a.shape[0], min, max, view_out.buf);
  PyBuffer_Release(&view_a);
  PyBuffer_Release(&view_out);
  return Py_BuildValue("");
}

函數登記處添加

{"clip", py_clip, METH_VARARGS, "clip array"},

則可，實際測試發現，手動接口效率也很高，和numpy同一水平。