Ошибка переноса GSL в Cython: невозможно преобразовать объект Python в 'double (*)(double, void *) nogil "
Прежде всего, прошу прощения за плохое использование терминологии программирования и за плохие (или даже совершенно неправильные) методы программирования; я все еще пытаюсь найти свои ноги:-)
Таким образом, я пытаюсь написать версию многомерной процедуры интеграции в Numeric Recipes, 3rd Edition, pgs. 198 - 199. на Cython, используя библиотеку численного интегрирования GSL. Я воспроизвел соответствующий фрагмент из книги NR ниже:
struct NRf3 {
double xsav, ysav;
double (*func3d)(const double, const double, const double);
double operator()(const double z)
{
return func3d(xsav, ysav, z);
}
}
struct NRf2 {
NRf3 f3;
double (*z1)(double, double);
double (*z2)(double, double);
NRf2(double zz1(double, double), double zz2(double, double)): z1(zz1), z2(zz2)} {}
double operator()(const double y)
{
f3.ysav = y;
return qgaus(f3, z1(f3.xsav, y), z2(f3.xsav, y));
}
}
struct NRf1 {
double (*y1)(double);
double (*y2)(double);
NRf2 f2;
NRf1(double yy1(double), yy2(double), z1(double, double), z2(double, double)) : y1(yy1), y2(yy2), f2(z1, z2) {}
double operator()(const double x){
f2.f3.xsav = x;
return qgaus(f2, y1(x), y2(x));
}
}
template <class T>
double quad3d(T &func, const double x1, const double x2, double y1(double), double y2(double), double z1(double, double), double z2(double, double))
{
NRf1 f1(y1, y2, z1, z2);
f1.f2.f3.func3d = func;
return qgaus(f1, x1, x2)
}
Этот скрипт применяет ранее написанную процедуру интеграции qgaus (который я позже заменю подпрограммой из GSL), которая принимает двойную функцию и два двойных, представляющих пределы интегрирования в качестве входных параметров.
Математически код выполняет интеграцию, как показано в формуле. (4.6.2) этой ссылки. Пределы интегрирования y[i] и z[i] имеют зависимости (x) и (x,y) соответственно, потому что интеграция выполняется в предварительно определенной области в плоскости (x,y) (см. Рис. 4.6.1). по той же ссылке выше).
Моя собственная попытка в Cython ниже (обратите внимание, что моя проблема, я думаю, не имеет ничего общего с интерфейсом GSL, но все еще находится на уровне правильного Cythoning):
from cython_gsl cimport *
from numpy import *
ctypedef double * double_ptr
ctypedef void * void_ptr
cdef double sigma = 0.045
cdef double mu = 0.37
## Test : a test function in x, y, z
cdef double gauss(double x, double y, double z, void * params) nogil:
cdef double mu = (<double_ptr> params)[0]
cdef double sigma = (<double_ptr> params)[1]
return mu * x + sigma * y + z
cdef class NRf3:
cdef double xsav, ysav
cdef double (*func3d)(double x, double y, double z, void * params)
def __init__(self, xsav, ysav):
self.xsav = xsav
self.ysav = ysav
def func(self, double z):
xsav = self.xsav
ysav = self.ysav
return func3d(xsav, ysav, z, void * params)
## Calling the gsl integrator qags
## Function: int gsl_integration_qags (const gsl_function * f, double a, double b, double epsabs, double epsrel, size_t limit, gsl_integration_workspace * workspace, double * result, double * abserr)
cdef class NRf2:
cdef NRf3 f3
cdef double z1
cdef double z2
cdef double paramsf2[1]
cdef size_t nevalf2
def __init__(self, double z1, double z2):
self.z1 = z1
self.z2 = z2
def __call__(self, double y):
cdef double resultf2, errorf2
cdef gsl_integration_workspace * W
W = gsl_integration_workspace_alloc(1000)
cdef gsl_function F2
z1 = self.z1
z2 = self.z2
f3.ysav = y
f3_temp = f3(f3.xsav, y)
F2.function = f3_temp
F2.params = self.paramsf2
gsl_integration_qags(&F2, z1, z2, 1e-2, 1e-2, 1000, W, &resultf2, &errorf2) # this line calls the GSL routine
gsl_integration_workspace_free(W)
return resultf2
cdef class NRf1:
cdef NRf2 f2
cdef double y1
cdef double y2
cdef double paramsf1[1]
cdef size_t nevalf1
def __init__(self, double yy1, double yy2, double z1, double z2):
self.y1 = yy1
self.y2 = yy2
self.z1 = z1
self.z2 = z2
self.f2 = f2(self.z1, self.z2)
def __call__(self, func, double x):
cdef double result, error
cdef gsl_integration_workspace * W
W = gsl_integration_workspace_alloc(1000)
cdef gsl_function F1
f2 = self.f2
f2.f3.func3d = func
f2.f3.xsav = x
y1 = self.y1
y2 = self.y2
F1.function = f2
F1.params = self.paramsf1
gsl_integration_qags(&F1, y1, y2, 1e-2, 1e-2, 1000, W, &result, &error)
gsl_integration_workspace_free(W)
return result
Компилируем этот скрипт с обычным
python min_example.py build_ext --inplace
как обнаружено на сайте Cython, выдает следующие ошибки:
Error compiling Cython file:
------------------------------------------------------------ ...
z1 = self.z1
z2 = self.z2
f3.ysav = y
f3_temp = f3(f3.xsav, y)
F2.function = f3_temp
^
------------------------------------------------------------
min_example.pyx:61:29: Cannot convert Python object to 'double
(*)(double, void *) nogil'
Error compiling Cython file:
------------------------------------------------------------ ...
f2.f3.func3d = func
f2.f3.xsav = x
y1 = self.y1
y2 = self.y2
F1.function = f2
^
------------------------------------------------------------
min_example.pyx:92:24: Cannot convert Python object to 'double
(*)(double, void *) nogil'
Ошибки относятся к строкам F2.function = f3_temp а также F1.function = f2 соответственно. В этом примере gsl_integration_qags это библиотека от gsl, но для целей этой проблемы правильность взаимодействия не имеет значения, так как (я думаю) проблема заключается в неправильном синтаксисе Cython, который я использовал.
Моя интуиция говорит мне, что это очень тривиальная ошибка новичка, но я не смог найти причину. Любой вклад (по проблеме, в мой код в целом) будет очень приветствоваться.
Изменить 1: чтобы отразить фактические сообщения об ошибках
Редактировать 2: отредактированный заголовок для отражения характера сообщений об ошибках
Изменить 3: я должен сказать, что я понимаю, что ошибка исходит от меня, пытаясь назначить экземпляр объекта Python, то есть в строках F2.function = f3_temp а также F1.function = f2к struct член, который является функцией, принимающей два аргумента. Поскольку это кажется довольно тривиальным переносом кода с C на Cython и при условии, что он тесно связан с исходной реализацией C, я просто задаюсь вопросом, как правильно написать эти строки.
Редактировать 4: Немного поиска поднимает эту тему. Я думаю, что моя путаница может быть устранена таким же образом: как передать член класса Python (который, если я его правильно создал) в функцию C?