Question

有时候自动化解决方案不起作用，你需要自己写一些定制的C或Fortran代码。这是有可能的，因为编译手段没有找到一些潜在的优化，或者因为你想要利用在Python中没有的库或语言特色。在所有这些情况中，你会需要使用外部函数接口，让你去访问用其他语言编写和编译的代码。

在本章其余部分，我们试图用一个外部库来解决2阶扩散方程，就以我们在第6章^[2]中所做的相同的方式。显示在例7-20中的这个库代码，能够代表你已经安装的库，或者一些你已经编写的代码。我们要看的方法发挥了巨大的作用，来提取出你的部分代码并把它们挪到另一种语言中去，以便做很目标化的基于语言的优化。

例7-20　解决2阶扩散问题的C代码示例

void evolve(double in[][512], double out[][512], double D, double dt) {
  int i, j;
   double laplacian;
   for (i=1; i<511; i++) {
    for (j=1; j<511; j++) {
     laplacian = in[i+1][j] + in[i-1][j] + in[i][j+1] + in[i][j-1]\
      - 4 * in[i][j];
     out[i][j] = in[i][j] + D * dt * laplacian;
    }
   }
}

为了使用这个代码，我们必须把它编译成一个创建了.so文件的共享模块。我们能使用gcc（或任何其他C编译器）通过以下步骤来做：

$ gcc -O3 -std=gnu99 -c diffusion.c $ gcc -shared -o diffusion.so diffusion.o

我们能把这个最终的共享库文件放置于可以被我们的Python代码所访问的任何地方，但是标准*nix组织把共享库存放于/usr/lib和/usr/local/lib。

7.13.1　ctypes

cPython^[3]中最基本的外部函数接口是通过ctypes模块。这个模块的特性就是具有很多限制性——你要负责做一切，并且需要用一段时间来确认你是按顺序做了这一切。这个额外级别的复杂性在我们的ctypes扩散版本中得到了证明，显示于例7-21中。

例7-21　ctypes 2阶扩散代码

import ctypes

grid_shape = (512, 512)
_diffusion = ctypes.CDLL("../diffusion.so") # ❶

# Create references to the C types that we will need to simplify future code
TYPE_INT = ctypes.c_int
TYPE_DOUBLE = ctypes.c_double
TYPE_DOUBLE_SS = ctypes.POINTER(ctypes.POINTER(ctypes.c_double))

# Initialize the signature of the evolve function to:
# void evolve(int, int, double**, double**, double, double)
_diffusion.evolve.argtypes = [
  TYPE_INT,
  TYPE_INT,
  TYPE_DOUBLE_SS,
  TYPE_DOUBLE_SS,
  TYPE_DOUBLE,
  TYPE_DOUBLE,
]
_diffusion.evolve.restype = None

def evolve(grid, out, dt, D=1.0):
  # First we convert the Python types into the relevant C types
  cX = TYPE_INT(grid_shape[0])
  cY = TYPE_INT(grid_shape[1])
  cdt = TYPE_DOUBLE(dt)
  cD = TYPE_DOUBLE(D)
  pointer_grid = grid.ctypes.data_as(TYPE_DOUBLE_SS) # ❷
  pointer_out = out.ctypes.data_as(TYPE_DOUBLE_SS)

  # Now we can call the function
  _diffusion.evolve(cX, cY, pointer_grid, pointer_out, cD, cdt) # ❸

❶ 这类似于导入diffusion.so库。

❷ grid和out都是numpy数组。

❸ 我们最终做完所有必要的步骤并能直接调用C函数。

我们做的第一件事就是“导入”我们的共享库。通过ctypes.CDLL调用来做。在这行中，我们可以声明Python能够访问的任何共享库（例如，ctypes-opencv模块装载libcv.so库）。由此，我们得到了一个_diffusion对象，包含了共享库所含有的所有成员。在这个例子中，diffusion.so只包含了一个函数——evolve，这不是一个对象的属性。如果diffusion.so中含有许多函数和属性，我们能通过_diffusion对象来访问它们全部。

无论怎样，即使_diffusion对象里有可调用的evolve函数，却不知道怎样来使用它。C是静态类型的，并且函数有很具体的签名。为了恰当地用evolve函数来工作，我们必须显式地设置输入参数类型和返回类型。当用Python接口串行地开发库，或者当处理一个快速变化的库时，这就会变得很枯燥。而且，既然ctypes不能检查你所给的是否是正确的类型，一旦你犯错，你的代码就可能默默地失败或发生段错误。

而且除了要设置参数和函数对象的返回类型，我们也需要小心地转换使用的数据（这叫“类型转换”）。我们传送给函数的每一个参数必须很小心地转换为本地的C类型。有时这事会变得相当诡异，因为Python的变量类型很灵活。例如，如果我们有num1 = 1e5，我们就知道这是一个Python的浮点数，因此我们应该使用ctypes.c_float。另一方面，对于num2 = 1e300来说，我们就不得不使用ctype.c_double，因为它对标准C浮点数会有溢出。

那就是说，numpy给它的数组提供了一个.ctypes属性来使它易与ctypes兼容。如果numpy没有提供这个功能，我们将不得不把一个ctypes数组初始化成正确类型，接着找到我们原始数据的位置并让我们的新ctypes对象指向那儿。

　警告　

当你正在把一个对象转换为一个ctype对象时，除非那个对象实现了一个缓存（就如array模块、numpy arrays、cStringIO等那样），不然你的数据将会被拷贝进新对象中。在把一个整数转换为浮点数的情况下，这并不会对你的代码性能有什么影响。无论如何，如果你正在转换一个很长的Python list，这会导致很大的性能惩罚！在这些情况下，使用array模块，或numpy array，或者干脆用struct模块构建你自己的缓存对象，这些都会有帮助。但是这样做会伤害你的代码可读性，因为这些对象一般没有它们的原生Python对应物那样灵活。

如果你必须要给库传送一个复杂的数据结构，情况就会变得愈加复杂。例如，如果你的库期待一个代表空间中的点的C结构，有x和y属性，你不得不定义：

from ctypes import Structure
class cPoint(Structure):
  _fields_ = ("x", c_int), ("y", c_int)

在这个点上，你能通过初始化一个cPoint对象（例如，point = cPoint(10, 5)）来开始创建一个兼容C的对象。这不是一项有可怕工作量的工作，但是却会变得乏味，并且产生一些代码碎片。如果一个新版本库发布了，稍微改了下结构体，会发生什么事呢？这会让你的代码很难维护，一般会导致僵硬的代码，开发者会决定从不去更新正在被使用的低层库。

基于这些理由，如果你已经很好地理解C并且想要能够微调接口的每一方面，使用ctypes模块是很正确的。它有很好的可移植性，因为它是标准库的一部分，如果你的任务简单，它就提供简单的解决方案。要仔细一点，因为ctypes解决方案的复杂性（类似低层解决方案）会很快变得难以管理。

7.13.2　cffi

意识到ctypes有时会相当难用，cffi意图简化程序员所使用的许多标准运算。它用一个内部的C解析器来做，能够理解函数和结构体定义。

作为结果，我们能仅仅写出C代码来定义我们希望使用的库的结构体，接着cffi将会为我们做所有重量级的工作：导入模块并确认我们给结果函数声明了正确的类型。事实上，如果有库的源码，这个工作几乎微不足道，因为头文件（用.h结尾的文件）会包含我们所需的所有相关定义。例7-22演示了cffi版本的2阶扩散代码。

例7-22　cffi的2阶扩散代码

from cffi import FFI

ffi = FFI()
ffi.cdef(r'''
   void evolve(
  int Nx, int Ny,
  double **in, double **out,
  double D, double dt
   ); # ❶
''')
lib = ffi.dlopen("../diffusion.so")

def evolve(grid, dt, out, D=1.0):
  X, Y = grid_shape
  pointer_grid = ffi.cast('double**', grid.ctypes.data) # ❷
  pointer_out = ffi.cast('double**', out.ctypes.data)
  lib.evolve(X, Y, pointer_grid, pointer_out, D, dt)

❶ 这个定义的内容能被正常地从你正在使用的库手册中，或者通过查看库的头文件获取。

❷ 尽管我们还是需要把非本地Python对象做类型转换来使用我们的C模块，语法对那些有C经验的人来说是很熟悉的。

在前面的代码中，我们可以把cffi的初始化过程看作两个步骤。首先，我们创建了一个FFI对象并且给出我们所需的全局C声明。除了函数签名以外，还可以包括数据类型。接着，我们使用dlopen把一个共享库导入它自己的名字空间，这是一个FFI的子空间。这意味着我们能够装载两个具有相同evolve函数的库，分别赋给变量lib1和lib2，然后独立地使用它们（这对于调试和剖析很给力）。

除了简单地导入C共享库以外，cffi允许你只写C代码，然后使用verify函数来即时编译。这有很多即时收益——你能够简单地把你的一小部分代码用C来重写，而不去调用独立的C库的庞大的机制。可做替换的是，如果有一个你希望使用的库，但是要求用一些C的胶水代码来让接口完美工作，你可以按例7-23显示的那样只是把它和你的cffi代码内联起来从而让一切都处于一个中心化的位置上。除此之外，既然代码是即时编译的，你可以给你需要编译的每块代码声明编译指令。需要注意的是，无论如何，当每次verify函数运行去做实际的编译时， 这种编译方式会有首次惩罚。

例7-23　内联2阶扩散代码的cffi

ffi = FFI()
ffi.cdef(r'''
  void evolve(
     int Nx, int Ny,
     double **in, double **out,
     double D, double dt
  );
''')
lib = ffi.verify(r'''
void evolve(int Nx, int Ny,
    double in[][Ny], double out[][Ny],
    double D, double dt) {
  int i, j;
  double laplacian;
  for (i=1; i<Nx-1; i++) {
    for (j=1; j<Ny-1; j++) {
      laplacian = in[i+1][j] + in[i-1][j] + in[i][j+1] + in[i][j-1]\
        - 4 * in[i][j];
      out[i][j] = in[i][j] + D * dt * laplacian;
    }
  }
}
''', extra_compile_args=["-O3",]) # ❶

❶ 既使我们正在即时编译这段代码，我们也可以提供相关的编译标志。

Verify功能的另一个好处是它与复杂的cdef声明交互得很好。例如，如果我们正要使用一个具有很复杂结构体的库，但是只是想使用它的一部分，我们能够使用部分结构定义。为此，我们在ffi.cdef中的结构体定义中添加一个…，并且在后面的verify中#include相关的头文件。

例如，假设我们正使用一个有complicated.h头文件的库，其中包含了一个类似于下面的结构体：

struct Point {
  double x;
  double y;
  bool isActive;
  char *id;
  int num_times_visited;
}

如果我们只需要关心x和y的属性，我们能写一些简单的cffi代码，只关心下面这些值：

from cffi import FFI

ffi = FFI()
ffi.cdef(r"""
  struct Point {
    double x;
    double y;
    ...;
  };
  struct Point do_calculation();
""")
lib = ffi.verify(r"""
  #include <complicated.h>
""")

我们能运行comlicated.h库的do_calculation函数，并且返回给我们一个Point对象，它的x和y属性是可以访问的。这种移植性很惊人，因为代码在具有不同的Point实现的系统上都能很好地工作，或者当complicated.h的新版本出来后也能很好地工作，只要它们都有x和y属性就可以。

所有这些优点让cffi成为在Python中用C代码来工作时的一个很优秀的工具。它比ctypes简单很多，然而当直接和外来函数接口一起用时，还是给予了你想要的等量的细粒度控制。

7.13.3　f2py

对于许多科学应用来说，Fortran还是一个黄金标准。尽管它作为一种通用语言的日子已经结束了，但它还是有很多优点来让编写矢量运算变得容易，并且运行相当快。另外，有很多性能数学库用Fortran编写（LAPACK、BLAS等），并且能够使用在你的Python代码中的性能关键部分。

对于这种情况，f2py提供了一个非常简单的把Fortran代码导入Python的方法。归因于Fortran的显式类型，这个模块能做得很简单。既然解析和理解类型是简单的，f2py就能容易地创建一个CPython模块，该模块依靠C中的本地外部函数支持来使用Fortran代码。这意味着当你要使用f2py时，你其实在自动生成一个知道怎样使用Fortran代码的C模块！这样的结果就是，许多在ctypes和cffi的方案中与生俱来的混淆就不存在了。

在例7-24中，我们能看到一些用于解决扩散方程的与f2py兼容的代码。事实上，所有本地Fortran代码都与f2py兼容。无论如何，函数参数注解（由!f2py为前缀的语句）简化了最终的Python模块并且会生成易于使用的接口。注解隐式地告诉f2py我们是否意图让一个参数只用于输出还是只用于输入，是让我们就地修改还是完全隐式地修改。隐式类型对于vectors的大小尤其有用：当Fortran可能需要让那些数字变显式时，我们的Python代码已经有这个信息在手上了。当我们把类型设为“隐式”时，f2py能够自动为我们填充那些值，基本上在最终的Python接口中是让它们保持隐式的。

例7-24　使用f2py注解的Fortran 2阶扩散代码

SUBROUTINE evolve(grid, next_grid, D, dt, N, M)
  !f2py threadsafe
  !f2py intent(in) grid
  !f2py intent(inplace) next_grid
  !f2py intent(in) D
  !f2py intent(in) dt
  !f2py intent(hide) N
  !f2py intent(hide) M
  INTEGER :: N, M
  DOUBLE PRECISION, DIMENSION(N,M) :: grid, next_grid
  DOUBLE PRECISION, DIMENSION(N-2, M-2) :: laplacian
  DOUBLE PRECISION :: D, dt

  laplacian = grid(3:N, 2:M-1) + grid(1:N-2, 2:M-1) + &
        grid(2:N-1, 3:M) + grid(2:N-1, 1:M-2) - 4 * grid(2:N-1, 2:M-1)
   next_grid(2:N-1, 2:M-1) = grid(2:N-1, 2:M-1) + D * dt * laplacian
END SUBROUTINE evolve

我们运行下面的命令来把代码构建成一个Python模块：

$ f2py -c -m diffusion --fcompiler=gfortran --opt='-O3' diffusion.f90

这会创建一个能被直接导入Python的diffusion.so文件。

如果我们交互性地运行结果模块，我们能看到f2py带给我们的好处，多亏了我们的注解和f2py解析Fortran代码的能力：

In [1]: import diffusion

In [2]: diffusion?
Type:    module
String form: <module 'diffusion' from 'diffusion.so'>
File:    .../examples/compilation/f2py/diffusion.so
Docstring:
This module 'diffusion' is auto-generated with f2py (version:2).
Functions:
  evolve(grid,next_grid,d,dt)
.

In [3]: diffusion.evolve?
Type:    fortran
String form: <fortran object>
Docstring:
evolve(grid,next_grid,d,dt)

Wrapper for evolve.

Parameters
grid : input rank-2 array('d') with bounds (n,m)
next_grid : rank-2 array('d') with bounds (n,m)
d : input float
dt : input float

这个代码显示出f2py生成的结果是自动文档化的，而且接口相当简化。例如，与其我们来抽取出vectors的大小，f2py已经能发现自动寻找这些信息的方法并且在最终的接口中只是把它隐藏了起来。事实上，最终的evolve函数签名看起来和我们在例6-14中所写的纯Python版本一模一样。

我们必须要仔细的唯一事情就是在内存中numpy array的顺序。既然绝大多数我们所用的numpy和Python集中于从C演变过来的代码，我们对于内存中的数据顺序一直使用C的惯例（称为行优先顺序）。Fortran使用了一种不同的惯例（列优先顺序），这样我们必须要确保让我们的vectors遵守惯例。这些顺序仅仅表明了对于一个2维数组，它的列或行在内存中是否是紧挨的^[4]。幸运的是，这仅仅意味着当我们声明vectors时，给numpy指明order=’F’参数就行。

　备忘　

行优先顺序和列优先顺序的差异意味着矩阵[ [1, 2], [3, 4] ]在内存中按行优先顺序排列成[1, 2, 3, 4]，按列优先顺序排列成[1, 3, 2, 4]。这个差异只是习惯上的，当使用恰当时对性能没有任何真正影响。

这就产生了下面我们使用Fortran子例程的代码。除了导入f2py衍生库和让我们的数据遵照显式的Fortran顺序之外，这个代码看上去和我们在例6-14中所使用的简直一模一样。

from diffusion import evolve

def run_experiment(num_iterations):
  next_grid = np.zeros(grid_shape, dtype=np.double, order='F') # ❶
  grid = np.zeros(grid_shape, dtype=np.double, order='F')

  # ... standard initialization ...

  for i in range(num_iterations):
    evolve(grid, next_grid, 1.0, 0.1)
    grid, next_grid = next_grid, grid

❶ Fortran在内存中以不同的顺序排列数字，所以我们必须记得把我们的numpy array设置成使用这个标准。

7.13.4　CPython模块

最后，我们总是能一路走到CPython API的层面，并且写一个CPython模块。这要求我们以和开发CPython一样的方式去写代码，并且需要小心对待我们的代码和CPython实现之间的所有交互。

这样会具有很优秀的移植性优势，取决于Python版本。我们不需要任何外部模块或库，只需一个C编译器和Python！无论怎样，这样并不一定能很好地扩展到Python的新版本中去。例如，让用Python 2.7写的CPython模块和Python 3一起工作。

然而为移植性付出了巨大的代价，你要负责你的Python代码和模块之间接口的所有方方面面。甚至只是为最简单的任务也要去写上几十行的代码。例如，为了和来自例7-20中的扩散库对接，我们必须写28行代码只是去读一个函数的参数并解析它（例7-25）。当然，这确实意味着你有惊人的细粒度去控制正在发生的事情。如此可以一直走下去直到能够手动为Python的垃圾收集去更改引用计数（当创建处理本地Python类型的CPython模块时，这是许多痛楚的根源）。就因为这样，最终的代码趋向于比其他的接口方式快几分钟。

　警告　

总而言之，这个方法应该留作最后的解决手段。虽然它为编写CPython模块提供了很多信息，但最终的代码不如其他潜在的方法那么可重用和可维护。在模块中做很细微的改动常常需要完全重写。事实上，我们包含了模块代码，并且需要setup.py来编译它（例7-26）作为警示。

例7-25　与2阶扩散库对接的CPython模块

// python_interface.c
// - cpython module interface for diffusion.c
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION

#include <Python.h>
#include <numpy/arrayobject.h>
#include "diffusion.h"

/* Docstrings */
static char module_docstring[] =
  "Provides optimized method to solve the diffusion equation";
static char cdiffusion_evolve_docstring[] =
  "Evolve a 2D grid using the diffusion equation";

PyArrayObject* py_evolve(PyObject* self, PyObject* args) {
  PyArrayObject* data;
  PyArrayObject* next_grid;
  double dt, D=1.0;

  /* The "evolve" function will have the signature:
  *     evolve(data, next_grid, dt, D=1)
  */
  if (!PyArg_ParseTuple(args, "OOd|d", &data, &next_grid, &dt, &D)) {
    PyErr_SetString(PyExc_RuntimeError, "Invalid arguments");
    return NULL;
  }

  /* Make sure that the numpy arrays are contiguous in memory */
  if (!PyArray_Check(data) || !PyArray_ISCONTIGUOUS(data)) {
    PyErr_SetString(PyExc_RuntimeError,"data is not a contiguous array.");
    return NULL;
  }
  if (!PyArray_Check(next_grid) || !PyArray_ISCONTIGUOUS(next_grid)) {
    PyErr_SetString(PyExc_RuntimeError,"next_grid is not a contiguous array.");
    return NULL;
  }

  /* Make sure that grid and next_grid are of the same type and have the same
  * dimensions
  */
  if (PyArray_TYPE(data) != PyArray_TYPE(next_grid)) {
    PyErr_SetString(PyExc_RuntimeError,
            "next_grid and data should have same type.");
    return NULL;
  }
  if (PyArray_NDIM(data) != 2) {
    PyErr_SetString(PyExc_RuntimeError,"data should be two dimensional");
    return NULL;
  }
  if (PyArray_NDIM(next_grid) != 2) {
    PyErr_SetString(PyExc_RuntimeError,"next_grid should be two dimensional");
    return NULL;
  }
  if ((PyArray_DIM(data,0) != PyArrayDim(next_grid,0)) ||
    (PyArray_DIM(data,1) != PyArrayDim(next_grid,1))) {
    PyErr_SetString(PyExc_RuntimeError,
            "data and next_grid must have the same dimensions");
    return NULL;
  }

  /* Fetch the size of the grid we are working with */
  const int N = (int) PyArray_DIM(data, 0);
  const int M = (int) PyArray_DIM(data, 1);

  evolve(
    N,
    M,
    PyArray_DATA(data),
    PyArray_DATA(next_grid),
    D,
    dt
  );

  Py_XINCREF(next_grid);
  return next_grid;
}

/* Module specification */
static PyMethodDef module_methods[] = {
/* { method name , C function , argument types , docstring } */
   { "evolve"   , py_evolve  , METH_VARARGS  , cdiffusion_evolve_docstring } ,
   { NULL     , NULL     , 0         , NULL }
};

/* Initialize the module */
PyMODINIT_FUNC initcdiffusion(void)
{
  PyObject *m = Py_InitModule3("cdiffusion", module_methods, module_docstring);
  if (m == NULL)
    return;
  /* Load `numpy` functionality. */
  import_array();
}

为了构建这个代码，我们需要创建一个setup.py脚本，使用distutils模块来确定构建代码的方式，这样才是兼容Python（例7-26）。除了标准distutils模块，numpy还提供了它自己的模块用以帮助在你的CPython模块中加入与numpy的整合。

例7-26　用于CPython模块扩散接口的setup文件

"""
setup.py for cpython diffusion module. The extension can be built by running

  $ python setup.py build_ext --inplace

which will create the __cdiffusion.so__ file, which can be directly imported into
Python.
"""

from distutils.core import setup, Extension
import numpy.distutils.misc_util

__version__ = "0.1"

cdiffusion = Extension(
  'cdiffusion',
  sources = ['cdiffusion/cdiffusion.c', 'cdiffusion/python_interface.c'],
  extra_compile_args = ["-O3", "-std=c99", "-Wall", "-p", "-pg", ],
  extra_link_args = ["-lc"],
)

setup (
  name = 'diffusion',
  version = __version__,
  ext_modules = [cdiffusion,],
  packages = ["diffusion", ],
  include_dirs = numpy.distutils.misc_util.get_numpy_include_dirs(),
)

生成的结果是一个cdiffusion.so文件，能被Python直接导入，而且使用起来相当简单。既然我们已经完全控制了最终函数的签名以及我们的C代码与库的精确交互方式，我们能够（以一些艰难的工作）创建出一个容易使用的模块。

from cdiffusion import evolve

def run_experiment(num_iterations):
  next_grid = np.zeros(grid_shape, dtype=np.double)
  grid = np.zeros(grid_shape, dtype=np.double)

  # ... standard initialization ...

  for i in range(num_iterations):
    evolve(grid, next_grid, 1.0, 0.1)
    grid, next_grid = next_grid, grid