python - 使用类方法作为 GSL QAGS 的被积函数

Question

我想将类的成员函数用于 gsl_function 的全局函数，但我不知道该怎么做。由于我现在只是关于 C 的非常基础的知识，我知道我必须将作为被积函数的类实例发送到 void 参数，但从技术上讲我不能在 cython 中对其进行编码，这是一个 cython 示例。

from cython_gsl cimport *

ctypedef double * double_ptr
ctypedef void * void_ptr

cdef double foo(double x, void * params) nogil:
    cdef double alpha, f
    alpha = (<double_ptr> params)[0]
    f = log(alpha*x) / sqrt(x)
    return f


def main():
    cdef gsl_integration_workspace * w
    cdef double result, error, expected, alpha
    w = gsl_integration_workspace_alloc (1000)

    expected = -4.0
    alpha = 1

    cdef gsl_function F
    F.function = &foo
    F.params = &alpha

    gsl_integration_qags (&F, 0, 1, 0, 1e-7, 1000, w, &result, &error)
    print "result          = % .18f\n" % result
    print "estimated error          = % .18f\n" % error

假设我有以下类(class):

from math import *
cdef class Foo(object):
    def __init__(self, double a=1.2, double b=0.6):
        self.a = a
        self.b = b
    def _integrand(self,double x):
         cdef double self.a = (<double_ptr> params)[0]
         cdef double self.b = (<double_ptr> params)[1]
         return self.a*log(x)+self.b/x**3
    def whole(self, double upper_limit=10,double lower_limit=0):
        cdef gsl_integration_workspace * w
        cdef double result, error, expected, alpha
        w = gsl_integration_workspace_alloc (1000)

        expected = -4.0
        params[0] = self.a
        params[1] = self.b

        cdef gsl_function F
        F.function = &self._integrand
        F.params = params

        gsl_integration_qags (&F, lower_limit,upper_limit , 0, 1e-7, 1000, w, &result, &error)
        return result,error

如何将 Foo._integrand 转换为可由 gsl 使用的全局函数？

Answer 1

在 Python 中，您可以使用 functools.partial 像函数一样传递方法，预定义 self 参数:

from functools import partial

foo = Foo()
integrand_func = partial(Foo._integrand, foo)

但您不能在 C 中定义这样的偏函数。因此，我会将被积函数外部化并将其定义为类的一个参数，而不是额外的方法。请注意，这就像使用静态方法一样，但是 Cython does not support it yet .另一个重要提示是使用 log 形式的 math.h。请参阅下面的原型(prototype):

#cython: wraparound=False
#cython: boundscheck=False
#cython: cdivision=True
#cython: nonecheck=False
cdef extern from "math.h":
    double log(double x) nogil

ctypedef double (*function)(double x, void *params)

cdef struct integrandFoo_p:
    double *a
    double *b

cdef struct gsl_function:
    function function
    void *params

cdef double integrandFoo(double x, void *params):
    cdef integrandFoo_p *p=<integrandFoo_p *>params
    cdef double a, b
    a = p.a[0]
    b = p.b[0]
    return a*log(x)+b/(x*x*x)

cdef void trapzd(gsl_function *F, double lower_limit, double upper_limit,
                 int num, double *result, double *error):
    cdef int i
    cdef double x
    f = F[0].function
    p = F[0].params
    result[0] = 0.
    for i in range(num+1):
        x = lower_limit + (upper_limit - lower_limit)*i/num
        if i==0 or i==num:
            result[0] += f(x, p)*0.5
        else:
            result[0] += f(x, p)

    error[0] = 0.

cdef class Foo(object):
    cdef double a
    cdef double b
    cdef double result
    cdef double error
    cdef function f
    def __init__(self):
        self.a = 1.2
        self.b = 0.6
        self.f = <function>integrandFoo
    cdef void integrate(self, double lower_limit=1, double upper_limit=10):
        #cdef gsl_integration_workspace * w
        cdef double result, error, expected, alpha
        cdef integrandFoo_p p
        cdef gsl_function F

        #w = gsl_integration_workspace_alloc(1000)
        p.a = &self.a
        p.b = &self.b
        expected = -4.0

        F.function = self.f
        F.params = &p

        trapzd(&F, lower_limit, upper_limit, 1000, &self.result, &self.error)
        #gsl_integration_qags(&F, lower_limit, upper_limit, 0, 1e-7,
                             #1000, w, &result, &error)
def main():
    foo = Foo()
    print 'HERE', foo.result, foo.error
    foo.integrate()
    print 'HERE', foo.result, foo.error