How to extend NumPy

那是动态的
和随机混乱。在两者之间的艺术。
- John A. Locke
科学是一个微分方程。宗教是一种边界条件。
- Alan Turing

Writing an extension module

虽然ndarray对象被设计为允许在Python中快速计算,但它也被设计为通用的并满足各种各样的计算需求。因此,如果绝对速度是必不可少的,则不能替换专门针对您的应用程序和硬件的精心制作的编译循环。这是numpy包含f2py的原因之一,因此可以使用易于使用的将(简单)C / C ++和(任意)Fortran代码直接转换为Python的机制。你鼓励使用和改进这种机制。本节的目的不是记录此工具,而是记录编写此工具依赖的扩展模块的更基本步骤。

当扩展模块被编写,编译并安装到Python路径(sys.path)中的某个地方时,代码可以被导入到Python中,就像它是一个标准的python文件一样。它将包含已在C代码中定义和编译的对象和方法。在Python中执行此操作的基本步骤有详细的文档,您可以在www.python.org在Python在线文档中找到更多信息。

除了Python C-API之外,还有一个用于NumPy的完整而丰富的C-API,允许在C级进行复杂的操作。然而,对于大多数应用程序,通常仅使用几个API调用。如果所有你需要做的是提取一个指针到内存以及一些形状信息传递到另一个计算例程,那么你将使用非常不同的调用,那么如果你正在尝试创建一个新的数组类型或添加一个新的数据类型为ndarrays。本章介绍了最常用的API调用和宏。

Required subroutine

只有一个函数必须在C代码中定义,以便Python将其用作扩展模块。该函数必须调用init {name},其中{name}是Python中模块的名称。此函数必须声明,以便在例程外部的代码中可见。除了添加你想要的方法和常量之外,这个子程序还必须包含对import_array()和/或import_ufunc()的调用,这取决于需要哪个C-API。一旦任何C-API子例程被实际调用,忘记放置这些命令将显示为丑陋的分段故障(崩溃)。实际上可能在单个文件中有多个init {name}函数,在这种情况下,多个模块将由该文件定义。但是,有一些技巧来使它正常工作,这里不包括。

最小的init{name}方法如下:

PyMODINIT_FUNC
init{name}(void)
{
   (void)Py_InitModule({name}, mymethods);
   import_array();
}

mymethods必须是一个数组(通常静态声明)的PyMethodDef结构,它包含方法名称,实际的C函数,一个指示方法是否使用关键字参数的变量,以及docstrings。这些在下一节中解释。如果你想向模块添加常量,那么你存储从Py_InitModule返回的值,它是一个模块对象。向模块添加项最常用的方法是使用PyModule_GetDict(module)获取模块字典。使用模块字典,你可以手动添加任何你喜欢的模块。向模块添加对象的一种更简单的方法是使用三个附加的Python C-API调用中的一个,不需要单独提取模块字典。这些都在Python文档中进行了说明,为了方便起见,这里重复一下:

int PyModule_AddObject(PyObject* module, char* name, PyObject* value)
int PyModule_AddIntConstant(PyObject* module, char* name, long value)
int PyModule_AddStringConstant(PyObject* module, char* name, char* value)

所有这三个函数都需要模块对象(Py_InitModule的返回值)。名称是标记模块中的值的字符串。根据调用的函数,参数是一个通用对象(PyModule_AddObject窃取对它的引用),整数常量或字符串常量。

Defining functions

传递给Py_InitModule函数的第二个参数是一个结构,它使得容易定义模块中的函数。在上面给出的示例中,mymethods结构将在文件中(通常在init {name}子例程之前)定义为:

static PyMethodDef mymethods[] = {
    { nokeywordfunc,nokeyword_cfunc,
      METH_VARARGS,
      Doc string},
    { keywordfunc, keyword_cfunc,
      METH_VARARGS|METH_KEYWORDS,
      Doc string},
    {NULL, NULL, 0, NULL} /* Sentinel */
}

mymethods数组中的每个条目都是一个PyMethodDef结构,包含1)Python名称,2)实现该函数的C函数,3)指示关键字是否被该函数接受的标志,以及4)函数的docstring。可以通过向该表添加更多条目来为单个模块定义任何数量的函数。最后一个条目必须全部为NULL,如图所示作为一个哨兵。Python查找此条目以了解模块的所有函数已定义。

完成扩展模块必须完成的最后一件事是实际写入执行所需功能的代码。有两种类型的函数:那些不接受关键字参数的函数,以及那些。

Functions without keyword arguments

不接受关键字参数的函数应写为:

static PyObject*
nokeyword_cfunc (PyObject *dummy, PyObject *args)
{
    /* convert Python arguments */
    /* do function */
    /* return something */
}

虚拟参数不在此上下文中使用,可以安全地忽略。args参数包含作为元组传递到函数的所有参数。你可以做任何你想要的在这一点,但通常最简单的方式来管理输入参数是调用PyArg_ParseTuple(args,format_string,addresses_to_C_variables ...)或PyArg_UnpackTuple (tuple,“name”,min,max,...)。有关如何使用第一个函数的好描述包含在Python C-API参考手册的第5.5节(解析参数和构建值)中。你应该特别注意使用转换函数在Python对象和C对象之间转换的“O&”格式。所有其他格式函数可以(大多数)被认为是这个一般规则的特殊情况。在NumPy C-API中定义了几个可能有用的转换器函数。特别地,PyArray_DescrConverter函数对于支持任意数据类型规范非常有用。此函数将任何有效的数据类型Python对象转换为PyArray_Descr *对象。记住要传入应该填写的C变量的地址。

有很多关于如何在NumPy源代码中使用PyArg_ParseTuple的例子。标准用法如下:

PyObject *input;
PyArray_Descr *dtype;
if (!PyArg_ParseTuple(args, "OO&", &input,
                      PyArray_DescrConverter,
                      &dtype)) return NULL;

请务必记住,当使用“O”格式字符串时,您将获得对对象的借用引用。然而,转换器功能通常需要某种形式的存储器处理。在此示例中,如果转换成功,dtype将保存对PyArray_Descr *对象的新引用,而输入将保存借用的引用。因此,如果此转换与另一个转换(例如整数)混合,并且数据类型转换成功,但整数转换失败,那么您将需要在返回之前将引用计数释放到数据类型对象。A typical way to do this is to set dtype to NULL before calling PyArg_ParseTuple and then use Py_XDECREF on dtype before returning.

处理输入参数后,实际执行工作的代码将被写入(可能需要调用其他函数)。C函数的最后一步是返回一些东西。如果遇到错误,则应返回NULL(确保实际设置了错误)。如果不应返回任何值,则递增Py_None并返回。如果单个对象应该返回,那么它将返回(确保您自己的引用它首先)。如果多个对象应该返回,那么你需要返回一个元组。Py_BuildValue(format_string,c_variables ...)函数使得从C变量构建Python对象的元组变得容易。特别注意格式字符串中'N'和'O'之间的区别,或者你可以很容易地创建内存泄漏。'O'格式字符串递增其对应的PyObject * C变量的引用计数,而“N”格式字符串对相应的PyObject * C变量的引用。你应该使用'N',如果你已经创建了一个对象的引用,只是想给该元组的引用。你应该使用'O',如果你只有一个借用引用一个对象,需要创建一个提供元组。

Functions with keyword arguments

这些函数非常类似于没有关键字参数的函数。唯一的区别是函数签名是:

static PyObject*
keyword_cfunc (PyObject *dummy, PyObject *args, PyObject *kwds)
{
...
}

kwds参数包含一个Python字典,其键是关键字参数的名称,其值是相应的关键字参数值。这本字典可以处理,但你认为合适。然而,处理它的最简单的方法是用PyArg_ParseTupleAndKeywords(args,kwds,format_string())的调用替换PyArg_ParseTuple(args,format_string,addresses ...) ,char * kwlist [],addresses ...)。此函数的kwlist参数是提供预期关键字参数的字符串的NULL终结数组。format_string中的每个条目都应该有一个字符串。如果传入无效的关键字参数,则使用此函数将引发TypeError。

有关此函数的更多帮助,请参阅Python文档中扩展和嵌入教程的第1.8节(扩展函数的关键字参数)。

Reference counting

编写扩展模块时最大的困难是引用计数。这是f2py,编织,Cython,ctypes等的普及的一个重要原因。如果你错误处理引用计数你可以得到问题从内存泄漏到分割故障。我知道处理参考计数正确的唯一策略是血,汗和眼泪。首先,你强制它进入你的头,每个Python变量有一个引用计数。然后,你明确了解每个函数对对象的引用计数,所以你可以正确地使用DECREF和INCREF当你需要它们。参考计数可以真正测试你对你的编程工艺的耐心和勤奋的数量。尽管粗糙的描述,大多数情况下的引用计数是相当简单的,最常见的困难是没有使用DECREF对象,从一个例程由于一些错误提前退出。第二,是不拥有对传递给将窃取引用的函数或宏的对象的引用的常见错误(,例如 PyTuple_SET_ITEM,大多数函数需要PyArray_Descr对象)。

通常,当创建一个变量时,你会得到一个新的引用,或者是一些函数的返回值(但是有一些突出的异常,例如从一个元组或字典中获取一个项目)。当你拥有引用时,你负责确保当变量不再需要时调用Py_DECREF(var),并且没有其他函数“引用”它的引用。此外,如果你传递一个Python对象到一个函数,将“窃取”引用,那么你需要确保你自己的(或使用Py_INCREF来获得自己的引用)。你也会遇到借用参考的概念。借用引用的函数不会改变对象的引用计数,并且不希望“保持”引用。它只是暂时使用对象。当你使用PyArg_ParseTuplePyArg_UnpackTuple你收到一个借用引用的元组中的对象,不应该改变它们的函数内的引用计数。随着练习,你可以学习得到引用计数权,但它可以是令人沮丧的起初。

引用计数错误的一个常见来源是Py_BuildValue函数。请注意“N”格式字符和“O”格式字符之间的差异。如果你在你的子程序中创建一个新对象(例如输出数组),并且你将它返回一个返回值的元组,那么你应该在Py_BuildValue“O”字符将引用计数增加1。这将为调用者留下一个全新的数组的两个引用计数。当变量被删除并且引用计数减1时,仍然会有额外的引用计数,并且数组永远不会被释放。你会有一个引用计数诱导内存泄漏。使用“N”字符将避免这种情况,因为它将返回一个具有单个引用计数的对象(在元组内部)。

Dealing with array objects

大多数NumPy的扩展模块将需要访问一个ndarray对象(或它的一个子类)的内存。最简单的方法不需要你知道NumPy的内部。方法是

  1. 确保你正在处理一个良好的数组(对齐,以机器字节顺序和单段),具有正确的类型和数量的维度。

    1. 通过使用PyArray_FromAny或基于它构建的宏来转换某些Python对象。
    2. 通过使用PyArray_NewFromDescr或基于它的更简单的宏或函数构造所需形状和类型的新ndarray。
  2. 获取数组的形状和指向其实际数据的指针。

  3. 将数据和形状信息传递到实际执行计算的子例程或其他代码段。

  4. 如果你正在编写算法,那么我建议你使用数组中包含的stride信息来访问数组的元素(PyArray_GETPTR宏使这个无痛)。然后,您可以放宽您的要求,以免强制单段数据组和可能导致的数据复制。

以下各小节将介绍这些子主题。

Converting an arbitrary sequence object

从任何可以转换为数组的Python对象获取数组的主例程是PyArray_FromAny此函数对于许多输入参数非常灵活。几个宏使得更容易使用基本功能。PyArray_FROM_OTF可以说是这些宏中最常用的最有用的。它允许你将任意的Python对象转换为特定内置数据类型的数组(,例如 float),同时指定一组特定的需求(,例如 ,和可写)。语法是

PyObject *PyArray_FROM_OTF(PyObject* obj, int typenum, int requirements)

从任何可以转换为数组的Python对象obj返回一个ndarray。返回的数组中的维数由对象确定。返回的数组的所需数据类型在typenum中提供,它应该是枚举类型之一。返回的数组的要求可以是标准数组标志的任意组合。下面更详细地解释这些参数中的每一个。你在成功时收到数组的新引用。失败时,返回NULL,并设置异常。

obj

对象可以是可转换为ndarray的任何Python对象。如果对象已经是满足要求的ndarray(的子类),则返回新的引用。否则,构造一个新的数组。除非使用数组接口,否则obj的内容将复制到新数组中,以便不必复制数据。可以转换为数组的对象包括:1)任何嵌套序列对象,2)暴露数组接口的任何对象,3)任何具有__array__方法的对象(应返回一个ndarray) 4)任何标量对象(变成零维数组)。否则符合要求的ndarray的子​​类将被传递。如果你想确保基类ndarray,那么在需求标志中使用NPY_ENSUREARRAY只有在必要时才进行复印。如果您想要保证副本,请将NPY_ENSURECOPY传递到需求标志。

typenum

枚举类型之一或NPY_NOTYPE,如果数据类型应该从对象本身确定。可以使用基于C的名称:

或者,位宽名称可以在平台上支持使用。例如:

只有在不损失精度的情况下,对象才会被转换为所需类型。否则将返回NULL并引发错误。在需求标志中使用NPY_FORCECAST可覆盖此行为。

要求

用于ndarray的存储器模型允许在每个维度中的任意步幅前进到数组的下一个元素。然而,通常,您需要与期望C连续或Fortran连续内存布局的代码进行交互。此外,一个ndarray可能不对齐(元素的地址不是元素的大小的整数倍),这可能会导致程序崩溃(或至少工作更慢)如果你尝试和解引用指针到数组数据。这两个问题都可以通过将Python对象转换为一个数组,为您的具体使用更“行为良好”解决。

需求标志允许指定什么类型的数组是可接受的。如果传入的对象不满足此要求,那么将进行复制,以使返回的对象满足要求。这些ndarray可以使用非常通用的指针指向内存。此标志允许指定返回的数组的所需属性。所有的标志都在详细的API章节中解释。最常需要的标志是NPY_ARRAY_IN_ARRAYNPY_OUT_ARRAYNPY_ARRAY_INOUT_ARRAY

NPY_ARRAY_IN_ARRAY

等效于NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED这种标志的组合对于必须以C连续顺序并对齐的数组是有用的。这些类型的数组通常是一些算法的输入数组。

NPY_ARRAY_OUT_ARRAY

等效于NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE这种标志的组合对于指定以C连续顺序,对齐并且也可以写入的数组是有用的。这样的数组通常作为输出返回(尽管通常这种输出数组是从头开始创建的)。

NPY_ARRAY_INOUT_ARRAY

等效于NPY_ARRAY_C_CONTIGUOUS | NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE | NPY_ARRAY_UPDATEIFCOPY这种标志的组合对于指定将用于输入和输出的数组是有用的。如果需要复制,那么当删除临时(通过在接口例程结束时使用Py_DECREF)时,临时数组将被复制回传入的原数组中。使用NPY_ARRAY_UPDATEIFCOPY标志要求输入对象已是数组(因为其他对象不能以此方式自动更新)。如果发生错误,则使用设置了NPY_ARRAY_UPDATEIFCOPY标志的数组上的PyArray_DECREF_ERR(obj)这将删除数组,而不会将内容复制回原始数组。

可作为附加要求OR'的其他有用标志是:

NPY_ARRAY_FORCECAST

投射到所需类型,即使它不能在不丢失信息的情况下完成。

NPY_ARRAY_ENSURECOPY

确保生成的数组是原始数据的副本。

NPY_ARRAY_ENSUREARRAY

确保生成的对象是实际的ndarray,而不是子类。

注意

数组是否被字节交换由数组的数据类型确定。本地字节序数组总是由PyArray_FROM_OTF请求,因此在需求参数中不需要NPY_ARRAY_NOTSWAPPED标志。也没有办法从此例程获取字节交换数组。

Creating a brand-new ndarray

通常,必须从扩展模块代码中创建新的数组。也许一个输出数组是需要的,你不想让调用者提供它。也许只需要一个临时数组来进行中间计算。无论需要什么数据类型的需要有简单的方法来获得一个ndarray对象。执行此操作的最通用的函数是PyArray_NewFromDescr所有数组创建函数都会遍历这个重复使用的代码。由于它的灵活性,它可能有点混乱使用。因此,存在更容易使用的更简单的形式。

PyObject *PyArray_SimpleNew(int nd, npy_intp* dims, int typenum)

此函数分配新内存并将其放在具有nd维的ndarray中,其形状由dims t2指向的至少项的数组决定>。数组的内存未初始化(除非typenum为NPY_OBJECT,在这种情况下,数组中的每个元素都设置为NULL)。typenum参数允许指定任何内置数据类型,例如NPY_FLOATNPY_LONG如果需要,可以使用PyArray_FILLWBYTE(return_object,0)将数组的内存设置为零。

PyObject *PyArray_SimpleNewFromData(int nd, npy_intp* dims, int typenum, void* data)

有时,你想将其他地方分配的内存包装成一个ndarray对象供下游使用。这个例程使得直接这样做。前三个参数与PyArray_SimpleNew中的相同,最后一个参数是指向一个连续内存块的指针,它应该使用它作为它的数据缓冲区,将以C样式连续方式解释。返回对ndarray的新引用,但ndarray将不拥有其数据。当这个ndarray被释放时,指针不会被释放。

你应该确保所返回的数组存在时所提供的内存不会被释放。处理这个的最简单的方法是如果数据来自另一个引用计数的Python对象。在传递指针后,应增加此对象的引用计数,返回的ndarray的基本成员应指向拥有数据的Python对象。然后,当解除分配ndarray时,基本成员将适当地DECREF。如果你想要释放内存,只要ndarray释放,然后简单地设置OWNDATA标志返回的ndarray。

Getting at ndarray memory and accessing elements of the ndarray

如果obj是ndarray(PyArrayObject *),则ndarray的数据区域由void *指针指向PyArray_DATA(obj)或char *指针PyArray_BYTES(obj)。请记住(通常)此数据区可能不根据数据类型进行对齐,它可能表示字节交换数据,和/或它可能不可写。如果数据区域以原始字节顺序对齐,那么如何获取数组的特定元素只由npy_intp变量的数组决定,PyArray_STRIDES(obj)。特别地,整数的这个c数组示出了必须将多少字节添加到当前元素指针以获得每个维度中的下一个元素。对于小于4维的数组,有PyArray_GETPTR{k}(obj,...)宏,其中{k}是使用数组宽度的整数1,2,3或4 。参数....在数组中表示{k}非负整数索引。例如,假设E是三维参数。指向元素E[i,j,k]的(void *)指针被获得为PyArray_GETPTR3(E,i,j,k)。

如前所述,C风格的连续数组和Fortran风格的连续数组有特殊的跨步模式。两个数组标志(NPY_C_CONTIGUOUS和:cdata`NPY_F_CONTIGUOUS`)表示特定数组的步幅模式是否匹配C风格的连续或Fortran风格的连续或两者。可以分别使用PyArray_ISCONTIGUOUS(obj)和PyArray_ISFORTRAN(obj)来测试步幅模式是否匹配标准C或Fortran。大多数第三方库期望连续数组。但是,通常不难支持通用大步。我鼓励你尽可能使用自己的代码中的跨越信息,并保留单段需求来包装第三方代码。使用与ndarray一起提供的步幅信息,而不是要求连续的步幅减少否则必须进行的复制。

Example

以下示例显示了您可以如何编写接受两个输入参数(将转换为数组)和输出参数(必须是数组)的包装器。该函数返回None并更新输出数组。

static PyObject *
example_wrapper(PyObject *dummy, PyObject *args)
{
    PyObject *arg1=NULL, *arg2=NULL, *out=NULL;
    PyObject *arr1=NULL, *arr2=NULL, *oarr=NULL;

    if (!PyArg_ParseTuple(args, "OOO!", &arg1, &arg2,
        &PyArray_Type, &out)) return NULL;

    arr1 = PyArray_FROM_OTF(arg1, NPY_DOUBLE, NPY_IN_ARRAY);
    if (arr1 == NULL) return NULL;
    arr2 = PyArray_FROM_OTF(arg2, NPY_DOUBLE, NPY_IN_ARRAY);
    if (arr2 == NULL) goto fail;
    oarr = PyArray_FROM_OTF(out, NPY_DOUBLE, NPY_INOUT_ARRAY);
    if (oarr == NULL) goto fail;

    /* code that makes use of arguments */
    /* You will probably need at least
       nd = PyArray_NDIM(<..>)    -- number of dimensions
       dims = PyArray_DIMS(<..>)  -- npy_intp array of length nd
                                     showing length in each dim.
       dptr = (double *)PyArray_DATA(<..>) -- pointer to data.

       If an error occurs goto fail.
     */

    Py_DECREF(arr1);
    Py_DECREF(arr2);
    Py_DECREF(oarr);
    Py_INCREF(Py_None);
    return Py_None;

 fail:
    Py_XDECREF(arr1);
    Py_XDECREF(arr2);
    PyArray_XDECREF_ERR(oarr);
    return NULL;
}