Написание собственной ufunc
Создание новой универсальной функции
Перед чтением этого раздела рекомендуется ознакомиться с основами расширений C для Python, прочитав/пробежав глазами учебные материалы в разделе 1 Расширения и интеграция интерпретатора Python и в Как расширить NumPy.
Модуль umath — это сгенерированный компьютером модуль C, который создаёт множество ufunc. Он предоставляет множество примеров того, как создать универсальную функцию. Также несложно создать собственную ufunc, которая будет использовать механизм ufunc. Предположим, у вас есть функция, которую вы хотите применять к своим входным данным поэлементно. Создав новую ufunc, вы получите функцию, которая обрабатывает
- векторизацию
- многомерные циклы
- автоматическое преобразование типов с минимальным использованием памяти
- массивы выходных данных (по желанию)
Создание собственной ufunc не сложно. Всё, что требуется, — это цикл по одномерному массиву для каждого типа данных, который вы хотите поддерживать. Каждый одномерный цикл должен иметь определённую сигнатуру, и можно использовать только ufunc для типов данных фиксированного размера. Приведен ниже вызов функции для создания новой ufunc, работающей с встроенными типами данных. Для регистрации ufunc для пользовательских типов данных используется другой механизм.
В следующих разделах мы приводим примеры кода, которые можно легко модифицировать для создания собственных ufunc. Примеры постепенно усложняются и представляют собой различные версии функции logit, часто используемой в статистическом моделировании. Logit также интересна, поскольку благодаря особенностям стандартов IEEE (в частности, IEEE 754) все создаваемые ниже функции logit автоматически ведут себя следующим образом.
>>> logit(0) -inf >>> logit(1) inf >>> logit(2) nan >>> logit(-2) nan
Это замечательно, потому что автору функции не нужно вручную распространять inf или nan.
Пример расширения без ufunc
Для сравнения и общего ознакомления читателя мы предоставляем простое реализацию расширения C для функции logit, не использующей NumPy.
Для этого нам нужны два файла. Первый — файл 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. */
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;
}
Чтобы использовать файл setup.py, поместите файлы setup.py и spammodule.c в одну папку. Затем выполните python setup.py build для построения модуля для импорта или python 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])
После импорта модуля spam в Python вы можете вызвать logit через spam.logit. Обратите внимание, что функцию, используемую выше, нельзя применять к массивам NumPy напрямую. Для этого необходимо вызвать numpy.vectorize. Например, если интерпретатор Python открыт в файле, содержащем библиотеку spam, или spam был установлен, можно выполнить следующие команды:
>>> 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])
РЕЗУЛЬТИРУЮЩАЯ ФУНКЦИЯ LOGIT НЕ БЫСТРА! numpy.vectorize просто перебирает spam.logit. Цикл выполняется на уровне C, но массив NumPy постоянно анализируется и восстанавливается. Это дорогостоящая операция. При сравнении numpy.vectorize(spam.logit) и функций logit ufuncs, созданных ниже, функции logit ufuncs оказались примерно в 4 раза быстрее. Конечно, возможны и другие значения ускорения, в зависимости от характера функции.
Пример NumPy ufunc для одного типа данных
Для простоты мы приведём ufunc для одного типа данных — 'f8' double. Как и в предыдущем разделе, сначала мы приведём файл .c, а затем файл setup.py, используемый для создания модуля, содержащего ufunc.
Место в коде, соответствующее фактическим вычислениям для ufunc, отмечено /*BEGIN main ufunc computation*/ и /*END main ufunc computation*/. Код между этими строками — основная часть, которую нужно изменить, чтобы создать собственную 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};
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;
}
Это файл setup.py для приведенного выше кода. Как и раньше, модуль можно построить, вызвав python setup.py build в командной строке, или установить в site-packages, выполнив 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])
Пример NumPy ufunc с несколькими типами данных
Наконец, мы приведём пример полной ufunc с внутренними циклами для half-floats, floats, doubles и long doubles. Как и в предыдущих разделах, сначала мы приведём файл .c, а затем соответствующий файл setup.py.
Места в коде, соответствующие фактическим вычислениям для ufunc, отмечены /*BEGIN main ufunc computation*/ и /*END main ufunc computation*/. Код между этими строками — основная часть, которую нужно изменить, чтобы создать собственную 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};
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;
}
Это файл setup.py для приведенного выше кода. Как и раньше, модуль можно построить, вызвав python setup.py build в командной строке, или установить в site-packages, выполнив 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])
Пример NumPy ufunc с несколькими аргументами/значениями возврата
Наш последний пример — ufunc с несколькими аргументами. Это модификация кода для logit ufunc для данных с одним типом. Мы вычисляем (A*B, logit(A*B)).
Мы приводим только код C, так как файл setup.py точно такой же, что и в примере Пример NumPy ufunc для одного типа данных, за исключением строки
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. */
static void double_logitprod(char **args, npy_intp *dimensions,
npy_intp* steps, void* data)
{
npy_intp i;
npy_intp n = dimensions[0];
char *in1 = args[0], *in2 = args[1];
char *out1 = args[2], *out2 = args[3];
npy_intp in1_step = steps[0], in2_step = steps[1];
npy_intp out1_step = steps[2], out2_step = steps[3];
double tmp;
for (i = 0; i < n; i++) {
/*BEGIN main ufunc computation*/
tmp = *(double *)in1;
tmp *= *(double *)in2;
*((double *)out1) = tmp;
*((double *)out2) = log(tmp/(1-tmp));
/*END main ufunc computation*/
in1 += in1_step;
in2 += in2_step;
out1 += out1_step;
out2 += out2_step;
}
}
/*This a pointer to the above function*/
PyUFuncGenericFunction funcs[1] = {&double_logitprod};
/* These are the input and return dtypes of logit.*/
static char types[4] = {NPY_DOUBLE, NPY_DOUBLE,
NPY_DOUBLE, NPY_DOUBLE};
static void *data[1] = {NULL};
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, 2, 2,
PyUFunc_None, "logit",
"logit_docstring", 0);
d = PyModule_GetDict(m);
PyDict_SetItemString(d, "logit", logit);
Py_DECREF(logit);
return m;
}
Пример NumPy ufunc со структурированным типом данных аргументов
Этот пример демонстрирует, как создать ufunc для структурированного массива типа данных. В качестве примера мы показываем простую ufunc для сложения двух массивов с типом 'u8,u8,u8'. Процесс немного отличается от других примеров, так как вызов PyUFunc_FromFuncAndData не полностью регистрирует ufunc для пользовательских типов данных и типов данных структурированных массивов. Нам также нужно вызвать PyUFunc_RegisterLoopForDescr, чтобы завершить настройку ufunc.
Мы приводим только код C, так как файл setup.py точно такой же, что и в примере Пример NumPy ufunc для одного типа данных, за исключением строки
config.add_extension('npufunc', ['single_type_logit.c'])
которая заменяется на
config.add_extension('npufunc', ['add_triplet.c'])
Файл C приведён ниже.
#include "Python.h"
#include "math.h"
#include "numpy/ndarraytypes.h"
#include "numpy/ufuncobject.h"
#include "numpy/npy_3kcompat.h"
/*
* add_triplet.c
* This is the C code for creating your own
* NumPy ufunc for a structured array dtype.
*
* Details explaining the Python-C API can be found under
* 'Extending and Embedding' and 'Python/C API' at
* docs.python.org .
*/
static PyMethodDef StructUfuncTestMethods[] = {
{NULL, NULL, 0, NULL}
};
/* The loop definition must precede the PyMODINIT_FUNC. */
static void add_uint64_triplet(char **args, npy_intp *dimensions,
npy_intp* steps, void* data)
{
npy_intp i;
npy_intp is1=steps[0];
npy_intp is2=steps[1];
npy_intp os=steps[2];
npy_intp n=dimensions[0];
uint64_t *x, *y, *z;
char *i1=args[0];
char *i2=args[1];
char *op=args[2];
for (i = 0; i < n; i++) {
x = (uint64_t*)i1;
y = (uint64_t*)i2;
z = (uint64_t*)op;
z[0] = x[0] + y[0];
z[1] = x[1] + y[1];
z[2] = x[2] + y[2];
i1 += is1;
i2 += is2;
op += os;
}
}
/* This a pointer to the above function */
PyUFuncGenericFunction funcs[1] = {&add_uint64_triplet};
/* These are the input and return dtypes of add_uint64_triplet. */
static char types[3] = {NPY_UINT64, NPY_UINT64, NPY_UINT64};
static void *data[1] = {NULL};
static struct PyModuleDef moduledef = {
PyModuleDef_HEAD_INIT,
"struct_ufunc_test",
NULL,
-1,
StructUfuncTestMethods,
NULL,
NULL,
NULL,
NULL
};
PyMODINIT_FUNC PyInit_struct_ufunc_test(void)
{
PyObject *m, *add_triplet, *d;
PyObject *dtype_dict;
PyArray_Descr *dtype;
PyArray_Descr *dtypes[3];
m = PyModule_Create(&moduledef);
if (m == NULL) {
return NULL;
}
import_array();
import_umath();
/* Create a new ufunc object */
add_triplet = PyUFunc_FromFuncAndData(NULL, NULL, NULL, 0, 2, 1,
PyUFunc_None, "add_triplet",
"add_triplet_docstring", 0);
dtype_dict = Py_BuildValue("[(s, s), (s, s), (s, s)]",
"f0", "u8", "f1", "u8", "f2", "u8");
PyArray_DescrConverter(dtype_dict, &dtype);
Py_DECREF(dtype_dict);
dtypes[0] = dtype;
dtypes[1] = dtype;
dtypes[2] = dtype;
/* Register ufunc for structured dtype */
PyUFunc_RegisterLoopForDescr(add_triplet,
dtype,
&add_uint64_triplet,
dtypes,
NULL);
d = PyModule_GetDict(m);
PyDict_SetItemString(d, "add_triplet", add_triplet);
Py_DECREF(add_triplet);
return m;
}
Возвращаемый объект ufunc — вызываемый объект Python. Он должен быть помещён в (модульный) словарь под тем же именем, что и в аргументе name при создании ufunc. Следующий пример взят из модуля umath
static PyUFuncGenericFunction atan2_functions[] = {
PyUFunc_ff_f, PyUFunc_dd_d,
PyUFunc_gg_g, PyUFunc_OO_O_method};
static void* atan2_data[] = {
(void *)atan2f,(void *) atan2,
(void *)atan2l,(void *)"arctan2"};
static char atan2_signatures[] = {
NPY_FLOAT, NPY_FLOAT, NPY_FLOAT,
NPY_DOUBLE, NPY_DOUBLE, NPY_DOUBLE,
NPY_LONGDOUBLE, NPY_LONGDOUBLE, NPY_LONGDOUBLE
NPY_OBJECT, NPY_OBJECT, NPY_OBJECT};
...
/* in the module initialization code */
PyObject *f, *dict, *module;
...
dict = PyModule_GetDict(module);
...
f = PyUFunc_FromFuncAndData(atan2_functions,
atan2_data, atan2_signatures, 4, 2, 1,
PyUFunc_None, "arctan2",
"a safe and correct arctan(x1/x2)", 0);
PyDict_SetItemString(dict, "arctan2", f);
Py_DECREF(f);
...
© 2005–2022 NumPy Developers
Licensed under the 3-clause BSD License.
https://numpy.org/doc/1.21/user/c-info.ufunc-tutorial.html
