NumPy用户指南 >使用NumPy C-API > 编写自己的ufunc
在阅读本文之前,通过阅读/略读扩展和嵌入Python解释器的第1节以及如何扩展NumPy中的教程,可能有助于熟悉Python C扩展的基础知识。
umath模块是计算机生成的C模块,可创建许多ufunc。它提供了许多有关如何创建通用函数的示例。创建使用ufunc机制的自己的ufunc也不难。假设您有一个函数要在其输入中逐个元素地进行操作。通过创建一个新的ufunc,您将获得一个处理
广播
N维循环
自动类型转换,使用最少的内存
可选输出数组
创建自己的ufunc并不困难。对于您要支持的每种数据类型,只需要一个1-d循环即可。每个1-d循环必须具有特定的签名,并且只能使用固定大小的数据类型的ufunc。下面给出了用于创建新的ufunc以处理内置数据类型的函数调用。使用不同的机制为用户定义的数据类型注册ufunc。
在接下来的几节中,我们将提供示例代码,可以轻松对其进行修改以创建自己的ufunc。示例依次是logit函数的更完整或更复杂的版本,这是统计建模中的常用函数。Logit也很有趣,因为由于IEEE标准(特别是IEEE 754)的魔力,下面创建的所有logit函数都会自动具有以下行为。
>>> logit(0)
-inf
>>> logit(1)
inf
>>> logit(2)
nan
>>> logit(-2)
nan
这很棒,因为函数编写者不必手动传播infs或nans。
为了比较读者并使读者更容易理解,我们提供了一个不使用numpy的logit C扩展的简单实现。
为此,我们需要两个文件。第一个是包含实际代码的C文件,第二个是用于创建模块的setup.py文件。
#include <Python.h> #include <math.h> /* * spammodule.c * This is the C code for a non-numpy Python extension to * define the logit function, where logit(p) = log(p/(1-p)). * This function will not work on numpy arrays automatically. * numpy.vectorize must be called in python to generate * a numpy-friendly function. * * Details explaining the Python-C API can be found under * 'Extending and Embedding' and 'Python/C API' at * docs.python.org . */ /* This declares the logit function */ static PyObject* spam_logit(PyObject *self, PyObject *args); /* * This tells Python what methods this module has. * See the Python-C API for more information. */ static PyMethodDef SpamMethods[] = { {"logit", spam_logit, METH_VARARGS, "compute logit"}, {NULL, NULL, 0, NULL} }; /* * This actually defines the logit function for * input args from Python. */ static PyObject* spam_logit(PyObject *self, PyObject *args) { double p; /* This parses the Python argument into a double */ if(!PyArg_ParseTuple(args, "d", &p)) { return NULL; } /* THE ACTUAL LOGIT FUNCTION */ p = p/(1-p); p = log(p); /*This builds the answer back into a python object */ return Py_BuildValue("d", p); } /* This initiates the module using the above definitions. */ #if PY_VERSION_HEX >= 0x03000000 static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, "spam", NULL, -1, SpamMethods, NULL, NULL, NULL, NULL }; PyMODINIT_FUNC PyInit_spam(void) { PyObject *m; m = PyModule_Create(&moduledef); if (!m) { return NULL; } return m; } #else PyMODINIT_FUNC initspam(void) { PyObject *m; m = Py_InitModule("spam", SpamMethods); if (m == NULL) { return; } } #endif
要使用setup.py文件,请将setup.py和spammodule.c放在同一文件夹中。然后python setup.py build将构建要导入的模块,或者setup.py install将模块安装到您的site-packages目录。
''' setup.py file for spammodule.c Calling $python setup.py build_ext --inplace will build the extension library in the current file. Calling $python setup.py build will build a file that looks like ./build/lib*, where lib* is a file that begins with lib. The library will be in this file and end with a C library extension, such as .so Calling $python setup.py install will install the module in your site-packages file. See the distutils section of 'Extending and Embedding the Python Interpreter' at docs.python.org for more information. ''' from distutils.core import setup, Extension module1 = Extension('spam', sources=['spammodule.c'], include_dirs=['/usr/local/lib']) setup(name = 'spam', version='1.0', description='This is my spam package', ext_modules = [module1])
将垃圾邮件模块导入python后,您可以通过spam.logit调用logit。请注意,上面使用的函数不能按原样应用于numpy数组。为此,我们必须对其调用numpy.vectorize。例如,如果在包含垃圾邮件库的文件中打开了python解释器或已经安装了垃圾邮件,则可以执行以下命令:
>>> import numpy as np
>>> import spam
>>> spam.logit(0)
-inf
>>> spam.logit(1)
inf
>>> spam.logit(0.5)
0.0
>>> x = np.linspace(0,1,10)
>>> spam.logit(x)
TypeError: only length-1 arrays can be converted to Python scalars
>>> f = np.vectorize(spam.logit)
>>> f(x)
array([ -inf, -2.07944154, -1.25276297, -0.69314718, -0.22314355,
0.22314355, 0.69314718, 1.25276297, 2.07944154, inf])
结果登录功能不快!numpy.vectorize只是在spam.logit上循环。循环是在C级别完成的,但是numpy数组不断地被解析和构建。这很贵。当作者将numpy.vectorize(spam.logit)与下面构造的logit函数进行比较时,logit函数几乎快了4倍。当然,取决于功能的性质,更大或更小的加速是可能的。
为简单起见,我们为单个dtype(“ f8”双精度)提供ufunc。与上一节一样,我们首先提供.c文件,然后提供用于创建包含ufunc的模块的setup.py文件。
代码中与ufunc的实际计算相对应的位置用/ * BEGIN主ufunc计算* /和/ * END主ufunc计算* /标记。这些行之间的代码是必须首先更改以创建自己的ufunc的主要内容。
#include "Python.h" #include "math.h" #include "numpy/ndarraytypes.h" #include "numpy/ufuncobject.h" #include "numpy/npy_3kcompat.h" /* * single_type_logit.c * This is the C code for creating your own * NumPy ufunc for a logit function. * * In this code we only define the ufunc for * a single dtype. The computations that must * be replaced to create a ufunc for * a different function are marked with BEGIN * and END. * * Details explaining the Python-C API can be found under * 'Extending and Embedding' and 'Python/C API' at * docs.python.org . */ static PyMethodDef LogitMethods[] = { {NULL, NULL, 0, NULL} }; /* The loop definition must precede the PyMODINIT_FUNC. */ static void double_logit(char **args, npy_intp *dimensions, npy_intp* steps, void* data) { npy_intp i; npy_intp n = dimensions[0]; char *in = args[0], *out = args[1]; npy_intp in_step = steps[0], out_step = steps[1]; double tmp; for (i = 0; i < n; i++) { /*BEGIN main ufunc computation*/ tmp = *(double *)in; tmp /= 1-tmp; *((double *)out) = log(tmp); /*END main ufunc computation*/ in += in_step; out += out_step; } } /*This a pointer to the above function*/ PyUFuncGenericFunction funcs[1] = {&double_logit}; /* These are the input and return dtypes of logit.*/ static char types[2] = {NPY_DOUBLE, NPY_DOUBLE}; static void *data[1] = {NULL}; #if PY_VERSION_HEX >= 0x03000000 static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, "npufunc", NULL, -1, LogitMethods, NULL, NULL, NULL, NULL }; PyMODINIT_FUNC PyInit_npufunc(void) { PyObject *m, *logit, *d; m = PyModule_Create(&moduledef); if (!m) { return NULL; } import_array(); import_umath(); logit = PyUFunc_FromFuncAndData(funcs, data, types, 1, 1, 1, PyUFunc_None, "logit", "logit_docstring", 0); d = PyModule_GetDict(m); PyDict_SetItemString(d, "logit", logit); Py_DECREF(logit); return m; } #else PyMODINIT_FUNC initnpufunc(void) { PyObject *m, *logit, *d; m = Py_InitModule("npufunc", LogitMethods); if (m == NULL) { return; } import_array(); import_umath(); logit = PyUFunc_FromFuncAndData(funcs, data, types, 1, 1, 1, PyUFunc_None, "logit", "logit_docstring", 0); d = PyModule_GetDict(m); PyDict_SetItemString(d, "logit", logit); Py_DECREF(logit); } #endif
这是上述代码的setup.py文件。和以前一样,可以通过在命令提示符下调用python setup.py build来构建模块,也可以通过python setup.py install将其安装到站点软件包中。
''' setup.py file for logit.c Note that since this is a numpy extension we use numpy.distutils instead of distutils from the python standard library. Calling $python setup.py build_ext --inplace will build the extension library in the current file. Calling $python setup.py build will build a file that looks like ./build/lib*, where lib* is a file that begins with lib. The library will be in this file and end with a C library extension, such as .so Calling $python setup.py install will install the module in your site-packages file. See the distutils section of 'Extending and Embedding the Python Interpreter' at docs.python.org and the documentation on numpy.distutils for more information. ''' def configuration(parent_package='', top_path=None): import numpy from numpy.distutils.misc_util import Configuration config = Configuration('npufunc_directory', parent_package, top_path) config.add_extension('npufunc', ['single_type_logit.c']) return config if __name__ == "__main__": from numpy.distutils.core import setup setup(configuration=configuration)
安装完上述文件后,可以将其导入并按以下方式使用。
>>> import numpy as np
>>> import npufunc
>>> npufunc.logit(0.5)
0.0
>>> a = np.linspace(0,1,5)
>>> npufunc.logit(a)
array([ -inf, -1.09861229, 0. , 1.09861229, inf])
最后,我们给出一个完整ufunc的示例,其中包含用于半浮点数,浮点数,双打和长双打的内部循环。如前几节所述,我们首先提供.c文件,然后提供相应的setup.py文件。
代码中与ufunc的实际计算相对应的位置用/ * BEGIN主ufunc计算* /和/ * END主ufunc计算* /标记。这些行之间的代码是必须首先更改以创建自己的ufunc的主要内容。
#include "Python.h" #include "math.h" #include "numpy/ndarraytypes.h" #include "numpy/ufuncobject.h" #include "numpy/halffloat.h" /* * multi_type_logit.c * This is the C code for creating your own * NumPy ufunc for a logit function. * * Each function of the form type_logit defines the * logit function for a different numpy dtype. Each * of these functions must be modified when you * create your own ufunc. The computations that must * be replaced to create a ufunc for * a different function are marked with BEGIN * and END. * * Details explaining the Python-C API can be found under * 'Extending and Embedding' and 'Python/C API' at * docs.python.org . * */ static PyMethodDef LogitMethods[] = { {NULL, NULL, 0, NULL} }; /* The loop definitions must precede the PyMODINIT_FUNC. */ static void long_double_logit(char **args, npy_intp *dimensions, npy_intp* steps, void* data) { npy_intp i; npy_intp n = dimensions[0]; char *in = args[0], *out=args[1]; npy_intp in_step = steps[0], out_step = steps[1]; long double tmp; for (i = 0; i < n; i++) { /*BEGIN main ufunc computation*/ tmp = *(long double *)in; tmp /= 1-tmp; *((long double *)out) = logl(tmp); /*END main ufunc computation*/ in += in_step; out += out_step; } } static void double_logit(char **args, npy_intp *dimensions, npy_intp* steps, void* data) { npy_intp i; npy_intp n = dimensions[0]; char *in = args[0], *out = args[1]; npy_intp in_step = steps[0], out_step = steps[1]; double tmp; for (i = 0; i < n; i++) { /*BEGIN main ufunc computation*/ tmp = *(double *)in; tmp /= 1-tmp; *((double *)out) = log(tmp); /*END main ufunc computation*/ in += in_step; out += out_step; } } static void float_logit(char **args, npy_intp *dimensions, npy_intp* steps, void* data) { npy_intp i; npy_intp n = dimensions[0]; char *in=args[0], *out = args[1]; npy_intp in_step = steps[0], out_step = steps[1]; float tmp; for (i = 0; i < n; i++) { /*BEGIN main ufunc computation*/ tmp = *(float *)in; tmp /= 1-tmp; *((float *)out) = logf(tmp); /*END main ufunc computation*/ in += in_step; out += out_step; } } static void half_float_logit(char **args, npy_intp *dimensions, npy_intp* steps, void* data) { npy_intp i; npy_intp n = dimensions[0]; char *in = args[0], *out = args[1]; npy_intp in_step = steps[0], out_step = steps[1]; float tmp; for (i = 0; i < n; i++) { /*BEGIN main ufunc computation*/ tmp = *(npy_half *)in; tmp = npy_half_to_float(tmp); tmp /= 1-tmp; tmp = logf(tmp); *((npy_half *)out) = npy_float_to_half(tmp); /*END main ufunc computation*/ in += in_step; out += out_step; } } /*This gives pointers to the above functions*/ PyUFuncGenericFunction funcs[4] = {&half_float_logit, &float_logit, &double_logit, &long_double_logit}; static char types[8] = {NPY_HALF, NPY_HALF, NPY_FLOAT, NPY_FLOAT, NPY_DOUBLE,NPY_DOUBLE, NPY_LONGDOUBLE, NPY_LONGDOUBLE}; static void *data[4] = {NULL, NULL, NULL, NULL}; #if PY_VERSION_HEX >= 0x03000000 static struct PyModuleDef moduledef = { PyModuleDef_HEAD_INIT, "npufunc", NULL, -1, LogitMethods, NULL, NULL, NULL, NULL }; PyMODINIT_FUNC PyInit_npufunc(void) { PyObject *m, *logit, *d; m = PyModule_Create(&moduledef); if (!m) { return NULL; } import_array(); import_umath(); logit = PyUFunc_FromFuncAndData(funcs, data, types, 4, 1, 1, PyUFunc_None, "logit", "logit_docstring", 0); d = PyModule_GetDict(m); PyDict_SetItemString(d, "logit", logit); Py_DECREF(logit); return m; } #else PyMODINIT_FUNC initnpufunc(void) { PyObject *m, *logit, *d; m = Py_InitModule("npufunc", LogitMethods); if (m == NULL) { return; } import_array(); import_umath(); logit = PyUFunc_FromFuncAndData(funcs, data, types, 4, 1, 1, PyUFunc_None, "logit", "logit_docstring", 0); d = PyModule_GetDict(m); PyDict_SetItemString(d, "logit", logit); Py_DECREF(logit); } #endif
这是上述代码的setup.py文件。和以前一样,可以通过在命令提示符下调用python setup.py build来构建模块,也可以通过python setup.py install将其安装到站点软件包中。
''' setup.py file for logit.c Note that since this is a numpy extension we use numpy.distutils instead of distutils from the python standard library. Calling $python setup.py build_ext --inplace will build the extension library in the current file. Calling $python setup.py build will build a file that looks like ./build/lib*, where lib* is a file that begins with lib. The library will be in this file and end with a C library extension, such as .so Calling $python setup.py install will install the module in your site-packages file. See the distutils section of 'Extending and Embedding the Python Interpreter' at docs.python.org and the documentation on numpy.distutils for more information. ''' def configuration(parent_package='', top_path=None): import numpy from numpy.distutils.misc_util import Configuration from numpy.distutils.misc_util import get_info #Necessary for the half-float d-type. info = get_info('npymath') config = Configuration('npufunc_directory', parent_package, top_path) config.add_extension('npufunc', ['multi_type_logit.c'], extra_info=info) return config if __name__ == "__main__": from numpy.distutils.core import setup setup(configuration=configuration)
安装完上述文件后,可以将其导入并按以下方式使用。
>>> import numpy as np
>>> import npufunc
>>> npufunc.logit(0.5)
0.0
>>> a = np.linspace(0,1,5)
>>> npufunc.logit(a)
array([ -inf, -1.09861229, 0. , 1.09861229, inf])
我们的最后一个示例是带有多个参数的ufunc。它是对具有单个dtype的数据的logit ufunc的代码的修改。我们计算(A * B,logit(A * B))。
我们仅给出C代码,因为对于一个dtype,setup.py文件与示例NumPy ufunc中的setup.py文件完全相同,除了该行
config.add_extension('npufunc', ['single_type_logit.c'])
被替换为
config.add_extension('npufunc', ['multi_arg_logit.c'])
C文件如下。生成的ufunc具有两个参数A和B。它返回一个元组,其第一个元素为A * B,第二个元素为logit(A * B)。请注意,它自动支持广播以及ufunc的所有其他属性。
#include "Python.h" #include "math.h" #include "numpy/ndarraytypes.h" #include "numpy/ufuncobject.h" #include "numpy/halffloat.h" /* * multi_arg_logit.c * This is the C code for creating your own * NumPy ufunc for a multiple argument, multiple * return value ufunc. The places where the * ufunc computation is carried out are marked * with comments. * * Details explaining the Python-C API can be found under * 'Extending and Embedding' and 'Python/C API' at * docs.python.org . * */ static PyMethodDef LogitMethods[] = { {NULL, NULL, 0, NULL} }; /* The loop definition must precede the PyMODINIT_FUNC.