c - 使用 cython_gsl 进行集成(加上将 numpy.array 寻址到 cython)

Question

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

我的类(class)如下:

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 有 an example其中集成是在全局函数上完成的。

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 参数中发送类指针。