Создание собственной 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
Это замечательно, потому что автору функции не нужно вручную распространять infs или nans.
Пример расширения без 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 построит импортируемый модуль, или 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 основное вычисление 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};
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 с внутренними циклами для полуточных чисел, чисел с плавающей точкой, двойных и длинных двойных чисел. Как и в предыдущих разделах, мы сначала приводим файл .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};
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 с несколькими аргументами. Это модификация кода для ufunc logit для данных с одним типом данных. Мы вычисляем (A*B, logit(A*B)).
Мы приводим только код C, поскольку файл setup.py точно такой же, как и файл 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 со структурированным массивом dtype аргументов
Этот пример демонстрирует, как создать ufunc для структурированного массива dtype. В примере показана тривиальная ufunc для сложения двух массивов с dtype ‘u8,u8,u8’. Процесс немного отличается от других примеров, так как вызов PyUFunc_FromFuncAndData не полностью регистрирует ufunc для пользовательских типов данных и структурированных типов массивов. Нам также необходимо вызвать PyUFunc_RegisterLoopForDescr для завершения настройки ufunc.
Мы приводим только код C, поскольку файл setup.py точно такой же, как и файл 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–2020 NumPy Developers
Licensed under the 3-clause BSD License.
https://numpy.org/doc/1.19/user/c-info.ufunc-tutorial.html
