Cython中的Memoryview切片

安装和环境设置

从下面的地址下载并安装Cython 0.16：

Cython的Windows安装版下载地址

在c:\Python27\Lib\distutils路径下创建distutils.cfg并写入如下内容，这样Cython编写的扩展库将采用Python(x,y)自带的MinGW32进行编译：

[build]
compiler = mingw32

[build_ext]
compiler = mingw32

在Spyder中编写如下三个文件，并保存到同一个目录之下：

• memview_test.py用来测试编译之后的扩展库的测试程序
• memview.pyxCython源程序
• setup.py用于编译Cython源程序

下面是setup.py中的内容，它用编译memview.pyx，只需要输入setup.py build_ext --inplace即可将memview.pyx编译为memview.pyd。

import numpy as np
from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
    cmdclass = {'build_ext': build_ext},
    ext_modules = [
        Extension("memview", ["memview.pyx"],
                 extra_compile_args=['-march=i486']) ❶
    ],
    include_dirs=[np.get_include()]
)

❶通过extra_compile_args参数指定编译时传递给gcc的额外参数，如果没有-march=i486参数，在连接时会报错。

每次修改完pyx文件之后，都执行setup.py进行编译连接是一件很繁琐的事情，因此我们在memview_test.py中使用pyximport实现自动编译：

import numpy as np
import pyximport
from distutils.extension import Extension
exts = [Extension("memview", ["memview.pyx"],
                 extra_compile_args=['-march=i486'])]

pyximport.install(setup_args=dict(
    include_dirs=[np.get_include()],
    ext_modules=exts))

import memview

为了测试Cython编译环境是否正常，在memview.pyx中添加一行：print "hello world"。然后切换到memview_test.py按F5，弹出如下的执行配置窗口，选择下图中红框标示的项目。这样每次都在新的进程中运行程序，程序运行结束之后，进入交互模式，用户可以输入命令查看各个对象的内容。

如果一切正常，就可以看到新创建的控制台窗口中显示“hello world”。

Memoryview切片

Memoryview切片(Memoryview slices)是Cython中的一种特殊类型，通过它可以高效地访问支持buffer接口的Python对象内部的数据区，例如NumPy中的ndarray对象。下面我们通过一个例子说明它的用法。在memview.pyx中添加如下代码：

def memview_sum(int[:] a): ❶
    cdef int i
    cdef int s = 0
    for i in range(a.shape[0]): ❷
        s += a[i]
    return s

❶参数a是一个一维整数切片类型，可以将与此切片类型一致的数组传递给它。❷和NumPy数组一样，它的shape属性为其各个轴的长度。在命令行中输入cython.py -a memview.pyx将生成一个memview.html文件，它显示Cython源代码和编译之后的C语言代码之间的关系，并用颜色的深浅表示所生成的代码的长短，如下图所示。

分析所生成的C语言代码可知，Cython将第7行的s += a[i]翻译成了一大段C语言代码。这段代码可以处理下标越界以及下标为负数的情况。由于需要在循环中对每个元素进行判断，因此这些代码会降低运算速度。可以使用Cython提供的wraparound和boundscheck修饰器关闭这两项功能。

cimport cython

@cython.boundscheck(False)
@cython.wraparound(False)
def memview_sum(int[:] a):
    # ...

使用修饰器对整个函数体有效，如果只希望对某一段程序有效的话，可以使用with关键字：

for i in range(a.shape[0]):
with cython.boundscheck(False):
    with cython.wraparound(False):
        s += a[i]

关闭这两个选项之后，输出的代码如下：

__pyx_t_3 = __pyx_v_i;
__pyx_v_s = (__pyx_v_s + (*((int *) ( /* dim=0 */ (__pyx_v_a.data + __pyx_t_3 * __pyx_v_a.strides[0]) ))));

按F5运行memview_test.py，运行结束之后，在控制台中输入：

>>> a = np.arange(11)
>>> memview.memview_sum(a)
55
>>> memview.memview_sum(a[::2]) ❶
30

❶由于在C语言代码中使用data和strides属性访问数组中的数据，因此即使对于元素不是连续存储的数组也能正常运算。

如果希望数据在内存中是连续存储的，那么可以用int[::1] a定义：

def memview_sum2(int[::1] a):
    cdef int i
    cdef int s = 0
    for i in range(a.shape[0]):
        s += a[i]
    return s

在控制台中的执行结果为：

>>> memview.memview_sum2(a)
55
>>> memview.memview_sum2(a[::2])
...
ValueError: ndarray is not C-contiguous

多维内存视图切片类型可以用如下方式定义：

cdef int[:, :] # 二维切片
cdef int[:, ::1] # C语言连续(C-contiguous)的二维切片
cdef int[::1, :] # Fortran语言连续(Fortran-contiguous)的二维切片

内存视图对象

当将Cython的内存视图切片类型返回到Python中时，它就变成了一个内存视图对象。下面看一个例子。在memview.pyx中添加如下代码：

def memview_object(int[:, :] a): ❶
    return a[::2, ::2] ❷

❶参数a是一个二维切片类型，❷它除了支持整数下标之外，还可以通过切片下标生成新的切片对象。我们直接将新生成的切片对象返回。在控制台中执行如下代码：

>>> b = np.arange(24).reshape(6, 4)
>>> mv = memview.memview_object(b)
>>> mv

>>> dir(mv)
[..., 'T', 'base', 'copy', 'copy_fortran', 'is_c_contig', 'is_f_contig', 'itemsize',
'nbytes', 'ndim', 'shape', 'size', 'strides', 'suboffsets']
>>> mv.shape
(3, 2)

Cython中的切片类型和Python中的MemoryView对象拥有许多相同的属性，例如：itemsize, nbytes, ndim, shape, size, strides等。

MemoryView对象的这些属性和NumPy数组十分类似。我们可以通过numpy.asarray()将MemoryView转换为NumPy数组：

>>> c = np.asanyarray(mv)
>>> c
array([[ 0, 2],
     [ 8, 10],
     [16, 18]])
>>> c[0,0] = 100
>>> b[0,0]
100

新创建的数组c和原来的数组b共享内存，因此修改c[0,0]会同时修改b[0,0]。

Cython数组

当在Cython中调用切片对象的copy()和copy_fortran()时，将创建一个Cython数组，并在此数组上创建一个切片对象。在Python中调用MemoryView对象的copy()方法也与此类似。

>>> mv.base
array([[ 0, 1, 2, 3],
     [ 4, 5, 6, 7],
     [ 8, 9, 10, 11],
     [12, 13, 14, 15],
     [16, 17, 18, 19],
     [20, 21, 22, 23]])
>>> mv2 = mv.copy()
>>> mv2

>>> mv2.base

在上面的程序中，mv是一个NumPy数组的内存视图，调用mv.copy()之后得到一个新的内存视图mv2，它是一个Cython数组的内存视图。

我们也可以在Cython中直接创建Cython数组：

cimport cython.view
def cython_array(int w, int h):
    cdef int [:,:] m = cython.view.array(shape=(h, w), ❶
                             itemsize=sizeof(int), format="i")
    cdef int i, j
    for i in range(h):
        for j in range(w):
            m[i, j] = i+j ❷
    return m

❶创建一个形状为(h, w)的整型Cython数组，并用一个内存视图切片m对它进行存取。❷通过切片m我们可以设置Cython数组的内容。在memview_test.py中添加如下的程序将创建的数组用imshow()显示：

import pylab as pl
m = memview.cython_array(400, 250)
pl.imshow(m)
pl.show()

程序的输出如下图所示：

绘制Mandelbrot集

下面是使用Cython数组计算Mandelbrot集的程序：

cdef int mand_iter(double cx, double cy):
    cdef int i
    cdef double zx, zy, zx2, zy2
    zx = 0
    zy = 0
    for i in range(256):
        zx2 = cx + zx**2 - zy**2
        zy2 = cy + 2*zx*zy
        zx = zx2
        zy = zy2
        if zx**2+zy**2 > 1000:
            return i
    return i

@cython.boundscheck(False)
@cython.wraparound(False)
def mand(int w, int h, double xmin, double xmax, double ymin, double ymax):
    cdef int i, j
    cdef double x, y
    cdef int[:,:] m = cython.view.array(shape=(h, w),
        itemsize=sizeof(int), format="i")

    for i in range(h):
        for j in range(w):
            x = j*(xmax-xmin)/w + xmin
            y = i*(ymax-ymin)/h + ymin
            m[i, j] = mand_iter(x, y)
    return m

在memview_test.py中用下面的程序调用mand()：

import pylab as pl
import matplotlib.cm as cm

xmin = -2.0
xmax = 1.0
ymin = -1.2
ymax = 1.2
m = memview.mand(400, 400, xmin, xmax, ymin, ymax)
pl.imshow(m, extent=[xmin, xmax, ymin, ymax], cmap=cm.Blues_r)

程序的输出如下图所示：