使用Python作为胶水

一个基本的扩展模块创建源¶

引入f2py的最简单方法可能是提供一个简单示例。这是名为的文件中包含的子例程之一 add.f：

C
      SUBROUTINE ZADD(A,B,C,N)
C
      DOUBLE COMPLEX A(*)
      DOUBLE COMPLEX B(*)
      DOUBLE COMPLEX C(*)
      INTEGER N
      DO 20 J = 1, N
         C(J) = A(J)+B(J)
 20   CONTINUE
      END

该例程仅将元素添加到两个连续的数组中，然后将结果放在第三个数组中。所有三个阵列的内存必须由调用例程提供。f2py可以自动生成此例程的一个非常基本的接口：

f2py -m add add.f

假设您的搜索路径设置正确，您应该能够运行此命令。此命令将在当前目录中生成一个名为addmodule.c的扩展模块。现在可以像其他任何扩展模块一样从Python编译和使用此扩展模块。

创建编译的扩展模块¶

您还可以让f2py编译add.f并编译其产生的扩展模块，仅留下可以从Python导入的共享库扩展文件：

f2py -c -m add add.f

此命令在当前目录中保留一个名为add。{ext}的文件（其中{ext}是平台上python扩展模块的适当扩展名-pyd 等）。然后可以从Python导入此模块。它将为add中的每个子例程包含一个方法（zadd，cadd，dadd，sadd）。每个方法的文档字符串包含有关如何调用模块方法的信息：

>>> import add
>>> print add.zadd.__doc__
zadd - Function signature:
  zadd(a,b,c,n)
Required arguments:
  a : input rank-1 array('D') with bounds (*)
  b : input rank-1 array('D') with bounds (*)
  c : input rank-1 array('D') with bounds (*)
  n : input int

完善基本接口¶

默认接口是将fortran代码非常直观地翻译成Python。现在，Fortran数组参数必须是NumPy数组，整数参数应该是整数。接口将尝试将所有参数转换为所需的类型（和形状），如果不成功，则发出错误。但是，由于它对参数的语义一无所知（例如C是输出，n应该真正匹配数组的大小），因此有可能以可能导致Python崩溃的方式滥用此函数。例如：

>>> add.zadd([1,2,3], [1,2], [3,4], 1000)

在大多数系统上都会导致程序崩溃。在幕后，这些列表将被转换为适当的数组，但随后告知底层的add循环以循环的方式超出已分配内存的边界。

为了改善界面，应提供指令。这是通过构造接口定义文件来完成的。通常最好从f2py可以生成的接口文件开始（从中获取其默认行为）。要使f2py生成接口文件，请使用-h选项：

f2py -h add.pyf -m add add.f

此命令将文件add.pyf保留在当前目录中。该文件对应于zadd的部分是：

subroutine zadd(a,b,c,n) ! in :add:add.f
   double complex dimension(*) :: a
   double complex dimension(*) :: b
   double complex dimension(*) :: c
   integer :: n
end subroutine zadd

通过放置intent指令并检查代码，可以相当多地清理接口，直到Python模块方法更易于使用和更强大为止。

subroutine zadd(a,b,c,n) ! in :add:add.f
   double complex dimension(n) :: a
   double complex dimension(n) :: b
   double complex intent(out),dimension(n) :: c
   integer intent(hide),depend(a) :: n=len(a)
end subroutine zadd

Intent指令intent（out）用于告诉f2py这c是一个输出变量，应该在传递给基础代码之前由接口创建。intent（hide）指令告诉f2py不允许用户指定变量n，而是从的大小获取它a。将取决于（a）指令是必要告诉f2py n的值取决于输入（所以它不会尝试创建变量n，直至创建一个变量）。

修改后add.pyf，可以通过add.f95和编译生成新的python模块文件add.pyf：

f2py -c add.pyf add.f95

新接口具有文档字符串：

>>> import add
>>> print add.zadd.__doc__
zadd - Function signature:
  c = zadd(a,b)
Required arguments:
  a : input rank-1 array('D') with bounds (n)
  b : input rank-1 array('D') with bounds (n)
Return objects:
  c : rank-1 array('D') with bounds (n)

现在，可以以更强大的方式调用该函数：

>>> add.zadd([1,2,3],[4,5,6])
array([ 5.+0.j,  7.+0.j,  9.+0.j])

注意自动转换为正确的格式。

在Fortran源代码中插入指令¶

还可以通过将变量指令作为特殊注释放入原始fortran代码中来自动生成nice接口。因此，如果我将源代码修改为包含：

C
      SUBROUTINE ZADD(A,B,C,N)
C
CF2PY INTENT(OUT) :: C
CF2PY INTENT(HIDE) :: N
CF2PY DOUBLE COMPLEX :: A(N)
CF2PY DOUBLE COMPLEX :: B(N)
CF2PY DOUBLE COMPLEX :: C(N)
      DOUBLE COMPLEX A(*)
      DOUBLE COMPLEX B(*)
      DOUBLE COMPLEX C(*)
      INTEGER N
      DO 20 J = 1, N
         C(J) = A(J) + B(J)
 20   CONTINUE
      END

然后，我可以使用以下命令编译扩展模块：

f2py -c -m add add.f

函数add.zadd的结果签名与之前创建的签名完全相同。如果原始源代码包含和A(N)而不是and ，则只需将注释行放在源代码中，就可以获得（几乎）相同的接口 。唯一的区别是这是一个可选输入，默认为。A(*)BCINTENT(OUT) :: CNA

过滤例子¶

为了与要讨论的其他方法进行比较。这是使用固定平均滤波器对双精度浮点数的二维数组进行滤波的函数的另一个示例。从此示例中应该清楚知道使用Fortran索引到多维数组的优点。

      SUBROUTINE DFILTER2D(A,B,M,N)
C
      DOUBLE PRECISION A(M,N)
      DOUBLE PRECISION B(M,N)
      INTEGER N, M
CF2PY INTENT(OUT) :: B
CF2PY INTENT(HIDE) :: N
CF2PY INTENT(HIDE) :: M
      DO 20 I = 2,M-1
         DO 40 J=2,N-1
            B(I,J) = A(I,J) +
     $           (A(I-1,J)+A(I+1,J) +
     $            A(I,J-1)+A(I,J+1) )*0.5D0 +
     $           (A(I-1,J-1) + A(I-1,J+1) +
     $            A(I+1,J-1) + A(I+1,J+1))*0.25D0
 40      CONTINUE
 20   CONTINUE
      END

可以使用以下代码编译此代码并将其链接到名为filter的扩展模块中：

f2py -c -m filter filter.f

这将在当前目录中生成一个名为filter.so的扩展模块，该模块具有一个名为dfilter2d的方法，该方法返回输入的过滤版本。

从Python中调用f2py ¶

f2py程序是用Python编写的，可以在您的代码内部运行以在运行时编译Fortran代码，如下所示：

from numpy import f2py
with open("add.f") as sourcefile:
    sourcecode = sourcefile.read()
f2py.compile(sourcecode, modulename='add')
import add

源字符串可以是任何有效的Fortran代码。如果要保存扩展模块源代码，则可以通过source_fn关键字为编译功能提供合适的文件名。

自动扩展模块生成¶

如果要分发f2py扩展模块，则只需包括.pyf文件和Fortran代码。NumPy中的distutils扩展允许您完全根据此接口文件定义扩展模块。setup.py允许分发add.f模块的有效文件（作为软件包的一部分， f2py_examples以便作为加载f2py_examples.add）是：

def configuration(parent_package='', top_path=None)
    from numpy.distutils.misc_util import Configuration
    config = Configuration('f2py_examples',parent_package, top_path)
    config.add_extension('add', sources=['add.pyf','add.f'])
    return config

if __name__ == '__main__':
    from numpy.distutils.core import setup
    setup(**configuration(top_path='').todict())

使用以下命令即可轻松安装新软件包：

pip install .

假设您具有使用您所使用的Python版本写入主sitepackage目录的适当权限。为了使生成的包起作用，您需要创建一个名为的文件__init__.py （与相同的目录中add.pyf）。请注意，扩展模块完全根据add.pyf和add.f文件定义。.pyf文件到.c文件的转换由numpy.disutils处理。

结论¶

接口定义文件（.pyf）是您如何微调Python和Fortran之间的接口。在numpy / f2py / docs目录中找到关于f2py的体面文档，无论您在系统上安装了NumPy的位置如何（通常在站点软件包下）。在https://scipy-cookbook.readthedocs.io的“与其他语言交互”标题下，还提供了有关使用f2py的更多信息（包括如何使用它包装C代码）。

链接已编译代码的f2py方法是目前最复杂和集成的方法。它允许将Python与已编译的代码完全分开，同时仍允许扩展模块的单独分发。唯一的缺点是它需要存在Fortran编译器才能使用户安装代码。但是，由于存在免费编译器g77，gfortran和g95以及高质量的商用编译器，因此这种限制并不是特别繁重。我认为，Fortran仍然是编写快速清晰的代码以进行科学计算的最简单方法。它以最直接的方式处理复数和多维索引。但是请注意，某些Fortran编译器将无法优化代码以及良好的手写C代码。

复杂另外在用Cython ¶

这是一个名为Cython模块的一部分，该模块add.pyx实现了我们之前使用f2py实现的复杂加法功能：

cimport cython
cimport numpy as np
import numpy as np

# We need to initialize NumPy.
np.import_array()

#@cython.boundscheck(False)
def zadd(in1, in2):
    cdef double complex[:] a = in1.ravel()
    cdef double complex[:] b = in2.ravel()

    out = np.empty(a.shape[0], np.complex64)
    cdef double complex[:] c = out.ravel()

    for i in range(c.shape[0]):
        c[i].real = a[i].real + b[i].real
        c[i].imag = a[i].imag + b[i].imag

    return out

此模块显示了如何使用该cimport语句从numpy.pxdCython附带的标头中加载定义。看起来NumPy被导入了两次；cimport仅使NumPy C-API可用，而常规规则则import在运行时导致Python样式的导入，并可以调用熟悉的NumPy Python API。

该示例还演示了Cython的“类型化内存视图”，类似于C级别的NumPy数组，从某种意义上来说，它们是形状和跨步的数组，它们知道自己的范围（不同于通过裸指针寻址的C数组）。语法表示具有任意跨度的双精度的一维数组（向量）。一个连续的int数组将是，而一个浮点数矩阵将是。double complex[:]int[::1]float[:, :]

显示为cython.boundscheck带有注释的装饰器，该装饰器根据功能打开或关闭内存视图访问的边界检查。我们可以使用它来进一步加快代码速度，但要牺牲安全性（或在进入循环之前进行手动检查）。

除了视图语法外，该函数可立即由Python程序员读取。变量的静态类型i是隐式的。除了视图语法，我们还可以使用Cython的特殊NumPy数组语法，但最好使用视图语法。

在用Cython图像滤波器¶

我们使用Fortran创建的二维示例在Cython中也很容易编写：

cimport numpy as np
import numpy as np

np.import_array()

def filter(img):
    cdef double[:, :] a = np.asarray(img, dtype=np.double)
    out = np.zeros(img.shape, dtype=np.double)
    cdef double[:, ::1] b = out

    cdef np.npy_intp i, j

    for i in range(1, a.shape[0] - 1):
        for j in range(1, a.shape[1] - 1):
            b[i, j] = (a[i, j]
                       + .5 * (  a[i-1, j] + a[i+1, j]
                               + a[i, j-1] + a[i, j+1])
                       + .25 * (  a[i-1, j-1] + a[i-1, j+1]
                                + a[i+1, j-1] + a[i+1, j+1]))

    return out

该二维平均滤波器运行很快，因为循环位于C中，并且指针计算仅在需要时执行。如果上面的代码被编译为module image，则img可以使用以下代码非常快速地过滤出2-d图像：

import image
out = image.filter(img)

关于代码，有两点需要注意：首先，不可能将内存视图返回给Python。而是out先创建一个NumPy数组，然后使用b该数组的视图进行计算。其次，b键入视图。这意味着二维数组具有连续的行，即C矩阵顺序。明确指定顺序可以加快某些算法的速度，因为它们可以跳过跨步计算。double[:, ::1]

结论¶

Cython是几种科学的Python库（包括Scipy，Pandas，SAGE，scikit-image和scikit-learn）以及XML处理库LXML的首选扩展机制。语言和编译器维护良好。

使用Cython有几个缺点：

1. 在编写自定义算法时，有时在包装现有的C库时，需要对C有一定的了解。特别是，在使用C内存管理（malloc和朋友）时，很容易引入内存泄漏。但是，只需编译一个重命名为的Python模块就.pyx 可以加快速度，而添加一些类型声明可以在某些代码中显着提高速度。

2. 很容易在Python和C之间失去清晰的分隔，这使得将C代码重用于其他与Python不相关的项目变得更加困难。

3. Cython生成的C代码很难读取和修改（通常会在烦人但无害的警告下进行编译）。

Cython生成的扩展模块的一大优势是它们易于分发。总而言之，Cython是一种非常强大的工具，可用于粘合C代码或快速生成扩展模块，因此不应忽视。这对于不能或不会编写C或Fortran代码的人特别有用。

有一个共享库¶

对于可以与特定于平台的ctypes一起使用的共享库有一些要求。本指南假定您已熟悉在系统上建立共享库的过程（或只是拥有一个可用的共享库）。要记住的项目是：

• 共享库必须以特殊方式编译（例如，将-shared标志与gcc一起使用）。

• 在某些平台（例如 Windows）上，共享库需要一个.def文件，该文件指定要导出的功能。例如，mylib.def文件可能包含：

  LIBRARY mylib.dll
  EXPORTS
  cool_function1
  cool_function2
  

  或者，您可以__declspec(dllexport)在函数的C定义中使用存储类说明符 ，从而无需该.def文件。

Python distutils中没有以跨平台方式创建标准共享库（扩展模块是Python理解的“特殊”共享库）的标准方法。因此，在撰写本书时，ctypes的一大缺点是，很难以跨平台的方式分发使用ctypes并包含您自己的代码的Python扩展，该扩展应在用户系统上编译为共享库。 。

加载共享库¶

一种简单但可靠的加载共享库的方法是获取绝对路径名，并使用ctypes的cdll对象加载它：

lib = ctypes.cdll[<full_path_name>]

但是，在Windows上，访问cdll方法的属性将通过在当前目录或PATH中找到的名称加载第一个DLL。加载绝对路径名对于跨平台工作需要一些技巧，因为共享库的扩展会有所不同。有ctypes.util.find_library可用的实用程序可以简化查找要加载的库的过程，但并非万无一失。复杂的是，不同的平台具有共享库使用的默认扩展名（例如，.dll – Windows，.so – Linux，.dylib – Mac OS X）。如果您使用ctypes包装需要在多个平台上工作的代码，则也必须考虑到这一点。

NumPy提供了一个称为ctypeslib.load_library（名称，路径）的便捷功能 。此函数使用共享库的名称（包括“ lib”之类的任何前缀，但不包括扩展名）和共享库可位于的路径。它返回一个ctypes库对象，OSError如果找不到该库，则引发；ImportError如果ctypes模块不可用，则引发。（Windows用户：load_library始终使用cdecl调用约定来加载使用加载的ctypes库对象 。有关 在其他调用约定下加载库的方法，请参见ctypes.windll和/或下的ctypes文档ctypes.oledll。）

共享库中的函数可用作ctypes库对象的属性（从返回ctypeslib.load_library），也可以用作使用lib['func_name']语法的项。如果函数名称包含Python变量名称中不允许的字符，则后一种检索函数名称的方法特别有用。

转换参数¶

Python int / long，字符串和unicode对象会根据需要自动转换为等效的ctypes参数。None对象也会自动转换为NULL指针。所有其他Python对象必须转换为特定于ctypes的类型。围绕此限制有两种方法可以使ctype与其他对象集成。

1. 不要设置函数对象的argtypes属性，也不要_as_parameter_为要传递的对象定义 方法。该 _as_parameter_方法必须返回一个Python int，它将直接传递给函数。

2. 将argtypes属性设置为一个列表，该列表的条目包含带有名为from_param的类方法的对象，该方法知道如何将您的对象转换为ctypes可以理解的对象（int / long，字符串，unicode或具有该_as_parameter_属性的对象）。

NumPy会优先使用第二种方法和第二种方法，因为它可以更安全。ndarray的ctypes属性返回一个对象，该对象的_as_parameter_属性返回一个整数，该整数表示与之关联的ndarray的地址。结果，可以将此ctypes属性对象直接传递给期望指向ndarray中数据指针的函数。调用者必须确保ndarray对象的类型，形状正确，并且设置了正确的标志，否则，如果传入不适当数组的数据指针，则可能导致崩溃。

为了实现第二种方法，NumPy ndpointer在numpy.ctypeslib模块中提供了class-factory函数。这个classfactory函数产生一个合适的类，该类可以放在ctypes函数的argtypes属性条目中。该类将包含from_param方法，ctypes将使用该方法将传递给函数的任何ndarray转换为ctypes识别的对象。在此过程中，转换将检查用户在调用中指定的ndarray的任何属性ndpointer。可以检查的ndarray的方面包括所传递的任何数组上的标志的数据类型，维数，形状和/或状态。from_param方法的返回值是数组的ctypes属性（因为它包含_as_parameter_ 指向数组数据区域的属性）可以直接由ctypes使用。

ndarray的ctypes属性还具有其他属性，当将有关数组的其他信息传递给ctypes函数时，这些属性可能很方便。属性data， shape和stride可以提供与数组的数据区域，形状和步幅相对应的ctype兼容类型。data属性返回一个c_void_p代表指向数据区域的指针的。shape和stride属性每个都返回一个ctypes整数数组（如果为0-d数组，则返回None表示NULL指针）。数组的基本ctype是与平台上的指针大小相同的ctype整数。还有一些方法 data_as({ctype})，和shape_as(<base ctype>)strides_as(<base ctype>)。这些返回数据作为您选择的ctype对象，并使用您选择的基础基本类型返回shape / strides数组。为了方便起见，该ctypeslib模块还包含c_intp一个ctypes整数数据类型，其大小c_void_p与平台上的大小相同 （如果未安装ctypes，则其值为None）。

调用函数¶

该函数可以作为已加载共享库的属性或项目来访问。因此，如果./mylib.so有一个名为的函数 cool_function1，我可以通过以下方式访问此函数：

lib = numpy.ctypeslib.load_library('mylib','.')
func1 = lib.cool_function1  # or equivalently
func1 = lib['cool_function1']

在ctypes中，默认情况下，函数的返回值设置为'int'。可以通过设置函数的restype属性来更改此行为。如果该函数没有返回值（“ void”），请对None类型使用None：

func1.restype = None

如前所述，您还可以设置函数的argtypes属性，以使ctypes在调用函数时检查输入参数的类型。使用ndpointer工厂函数来生成现成的类，以检查新函数的数据类型，形状和标志。该ndpointer功能具有签名

ndpointer（dtype = None，ndim = None，shape = None，flags = None ）¶

None不检查具有该值的关键字参数。指定关键字可在转换为ctypes兼容对象时强制检查ndarray的这一方面。dtype关键字可以是任何理解为数据类型对象的对象。ndim关键字应为整数，shape关键字应为整数或整数序列。flags关键字指定传入的任何数组上所需的最小标志。可以将其指定为以逗号分隔的要求字符串，指示要求位或的“整数”或从变量的flags属性返回的flags对象。具有必要要求的阵列。

在argtypes方法中使用ndpointer类可以使使用ctypes和ndarray的数据区域调用C函数的安全性大大提高。您可能仍然希望将函数包装在其他Python包装器中，以使其易于使用（隐藏一些明显的参数并使某些参数输出参数）。在此过程中，requiresNumPy中的函数对于从给定输入返回正确类型的数组可能很有用。

完整的示例¶

在此示例中，我将展示如何使用ctypes来实现以前使用其他方法实现的加法函数和过滤器函数。第一，它实现了算法的C代码所包含的功能zadd，dadd，sadd，cadd，和dfilter2d。该zadd函数是：

/* Add arrays of contiguous data */
typedef struct {double real; double imag;} cdouble;
typedef struct {float real; float imag;} cfloat;
void zadd(cdouble *a, cdouble *b, cdouble *c, long n)
{
    while (n--) {
        c->real = a->real + b->real;
        c->imag = a->imag + b->imag;
        a++; b++; c++;
    }
}

具有类似的代码cadd，dadd以及sadd用于处理复杂的浮点，双精度和浮点数据类型，分别为：

void cadd(cfloat *a, cfloat *b, cfloat *c, long n)
{
        while (n--) {
                c->real = a->real + b->real;
                c->imag = a->imag + b->imag;
                a++; b++; c++;
        }
}
void dadd(double *a, double *b, double *c, long n)
{
        while (n--) {
                *c++ = *a++ + *b++;
        }
}
void sadd(float *a, float *b, float *c, long n)
{
        while (n--) {
                *c++ = *a++ + *b++;
        }
}

该code.c文件还包含以下功能dfilter2d：

/*
 * Assumes b is contiguous and has strides that are multiples of
 * sizeof(double)
 */
void
dfilter2d(double *a, double *b, ssize_t *astrides, ssize_t *dims)
{
    ssize_t i, j, M, N, S0, S1;
    ssize_t r, c, rm1, rp1, cp1, cm1;

    M = dims[0]; N = dims[1];
    S0 = astrides[0]/sizeof(double);
    S1 = astrides[1]/sizeof(double);
    for (i = 1; i < M - 1; i++) {
        r = i*S0;
        rp1 = r + S0;
        rm1 = r - S0;
        for (j = 1; j < N - 1; j++) {
            c = j*S1;
            cp1 = j + S1;
            cm1 = j - S1;
            b[i*N + j] = a[r + c] +
                (a[rp1 + c] + a[rm1 + c] +
                 a[r + cp1] + a[r + cm1])*0.5 +
                (a[rp1 + cp1] + a[rp1 + cm1] +
                 a[rm1 + cp1] + a[rm1 + cp1])*0.25;
        }
    }
}

该代码相对于Fortran等效代码可能具有的优势是，它可以任意跨步（即非连续数组），并且根据编译器的优化能力，其运行速度也可能更快。但是，它显然比中的简单代码复杂filter.f。此代码必须编译到共享库中。在我的Linux系统上，使用以下命令完成此操作：

gcc -o code.so -shared code.c

它将在当前目录中创建一个名为code.so的shared_library。在Windows上，不要忘记__declspec(dllexport)在每个函数定义之前的行的void前面添加，或编写一个