1

首先，请原谅我对编程术语的不当使用和糟糕（甚至完全错误）的编程实践；我仍在努力寻找自己的脚 :-)

总之，我正在尝试在 Numerical Recipes，第 3 版，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 函数？

4

0 回答 0