C语言 使用cython_gsl进行集成(加上寻址numpy.数组为cython)

Question

我想提高我的代码的速度，所以我将我的原始Python代码转移到Cython。代码的瓶颈是我用scipy.integral.quad做的集成部分，然后在二维数组的每一行上连接一堆numpy.array和sum。所以我试图重写我在Cython中使用的类。

我的班级如下:

import numpy as np
import cosmology
from cython_gsl cimport *
ctypedef double * double_ptr
ctypedef void * void_ptr

cdef class Cosmology(object):
    cdef double omega_c
    def __init__(self, double omega_m=0.3, double omega_lam=0.7):
        # no quintessence, no radiation in this universe!
        self.omega_m = omega_m
        self.omega_lam = omega_lam
        self.omega_c = (1. - omega_m - omega_lam)
    cdef double a(self, double z):
        return 1./(1+z)
    cdef double E(self, double a):
        return (self.omega_m*a**(-3) + self.omega_c*a**(-2) + self.omega_lam)**0.5
    cdef double __angKernel(self, double x, void * params) nogil:
        cdef double alpha, f
        alpha = (<double_ptr> params)[0]
        f=self.E(alpha*x**-1)**-1
        return f
    cdef double Da(self, int size, double *z, double z_ref=0):
        cdef gsl_integration_workspace *w
        cdef double result, error, expected, alpha
        w = gsl_integration_workspace_alloc (1000)
        cdef gsl_function F
        if size>1:
           da = np.zeros(size)
           for i in range(size):
               da[i] = self.Da(1, &z[i], z_ref)
           return da
        else:
           if z[0] < 0:
                    raise ValueError("Redshift z must not be negative")
           if z[0] < z_ref:
                    raise ValueError("Redshift z must not be smaller than the reference redshift")
           expected = -4.0
           alpha = 1
           F.function = &self.__angKernel
           F.params = &alpha
           gsl_integration_qags (&F, z_ref+1, z[0]+1, 0, 1e-7, 1000, w, &result, &error)
           d=result
           rk = (abs(self.omega_c))**0.5
           if (rk*d > 0.01):
              if self.omega_c > 0:
                 d = sinh(rk*d)/rk
              if self.omega_c < 0:
                 d = sin(rk*d)/rk
           return d/(1+z[0])

cdef class halo_positions(object):
    cdef double *x 
    cdef double *y 
    def __init__(self, double[:] positions):
        self.x = &positions[0]
        self.y = &positions[1]

我试图使用cython_gsl集成，之后我试图实现我的集成和构建cnfw_model.pyx代码。我得到了这个错误信息:

Error compiling Cython file:
------------------------------------------------------------
...
                    raise ValueError("Redshift z must not be smaller than the reference redshift")
           expected = -4.0
           alpha = 1
           F.function = &self.__angKernel
                       ^
------------------------------------------------------------
cnfw_model.pyx:102:24: Cannot assign type 'double (*)(Cosmology, double, void *) nogil' to 'double (*)(double, void *) nogil'
Error compiling Cython file:
------------------------------------------------------------
...
    cdef double __angKernel(self, double x, void * params) nogil:
        cdef double alpha, f
        alpha = (<double_ptr> params)[0]
        """Integration kernel for angular diameter distance computation.
        """
        f=self.E(alpha*x**-1)**-1
               ^
------------------------------------------------------------
cnfw_model.pyx:73:16: Calling gil-requiring function not allowed without gil

我不知道应该如何定义function F来避免这个问题。

second:在class halo_positions中，我想更改它以获得2xN维度数组并将其划分为两个1xN维度数组并将其分配给x和y。

Answer 1

似乎与c++中也发生的问题相似(我不使用python编写代码，但问题似乎是相同的)。不能将成员函数转换为全局函数。记住，gsl_function是下面的c_struct

struct gsl_function_struct 
{
   double (* function) (double x, void * params);
   void * params;
};
typedef struct gsl_function_struct gsl_function ;

你可以看到gsl。Function是指向全局函数的指针(成员函数有一个保存this指针的额外形参)。

cython_gsl的github有一个在全局函数上完成集成的示例。

from cython_gsl cimport *
ctypedef double * double_ptr
ctypedef void * void_ptr
cdef double normal(double x, void * params) nogil:
  cdef double mu = (<double_ptr> params)[0]
  cdef double sigma = (<double_ptr> params)[1]
  return gsl_ran_gaussian_pdf(x - mu, sigma)
def cdf_numerical(double x, double mu, double sigma):
  cdef double alpha, result, error, expected
  cdef gsl_integration_workspace * W
  W = gsl_integration_workspace_alloc(1000)
  cdef gsl_function F
  cdef double params[1]
  cdef size_t neval
  params[0] = mu
  params[1] = sigma
  F.function = &normal
  F.params = params
  gsl_integration_qag(&F, -10, x, 1e-2, 1e-2, 1000, GSL_INTEG_GAUSS15, W, &result, &error)
  gsl_integration_workspace_free(W)
  return result

解决这个问题的一种方法是创建一个全局函数来包装成员函数。在这种情况下，你需要在void形参中发送类指针。