python - 将指针定义为类的方法？

Question

我正在尝试加速我的 cython 代码。我遇到了this link作者描述了如何使用指针代替numpy数组来提高 cython 代码的速度。在我的宇宙学类(class)中，瓶颈是Da函数。我对C中的指针不太熟悉，如果有人给我一个想法，我将不胜感激:

是否可以将类的方法定义为指针，例如在我的例子中将 np.ndarray[double, ndim=1] Da 转换为 double* Da？

from __future__ import division
import numpy as np
cimport numpy as np
cimport cython
import copy
cdef extern from "gsl/gsl_math.h":
    ctypedef struct gsl_function:
        double (* function) (double x, void * params)
        void * params

cdef extern from "gsl/gsl_integration.h":
    ctypedef struct gsl_integration_workspace
    gsl_integration_workspace *  gsl_integration_workspace_alloc(size_t n)
    void  gsl_integration_workspace_free(gsl_integration_workspace * w)
    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 double func_callback(double x, void* params): 
     return (<cosmology>params).__angKernel(x) 

cdef class cosmology(object):
    cdef public double omega_m, omega_l, h, w, omega_r, G, v_c
    cdef object omega_c
    def __init__(self,double omega_m = 0.3, double omega_l = 0.7, double h = 0.7, double w = -1, double omega_r = 0., double G = std_G):

        self.omega_m = omega_m
        self.omega_l = omega_l
        self.omega_r = omega_r
        self.omega_c = (1. - omega_m - omega_l)
        self.h = h
        self.w = w
        self.G = G
        self.v_c = v_c

    def __copy__(self):

        return cosmology(omega_m = self.omega_m, omega_l = self.omega_l, h = self.h, w = self.w, omega_r = self.omega_r, G = self.G)

    property H0:
       def __get__(self):
           return 100*self.h  #km/s/MPC

    cpdef double a(self, double z):
        return 1./(1.+z)

    cpdef double E(self, double a):
        return (self.omega_r*a**(-4) + self.omega_m*a**(-3) + self.omega_c*a**(-2) + self.omega_l)**0.5

    @cython.boundscheck(False)
    @cython.wraparound(False)
    @cython.nonecheck(False)
    cdef double __angKernel(self, double x):
         """Integration kernel for angular diameter distance computation.
         """
         return self.E(x**-1)**-1

    @cython.boundscheck(False)
    @cython.wraparound(False)
    @cython.nonecheck(False)
    cpdef np.ndarray[double, ndim=1] Da(self, np.ndarray[double, ndim=1] z, double z_ref=0):
          cdef gsl_integration_workspace* w =gsl_integration_workspace_alloc(1000)

          cdef gsl_function F
          F.function = &func_callback
          F.params = <void*>self

          cdef double result = 3, error = 5
          cdef double err, rk, zs, omc
          omc=self.omega_c
          cdef np.ndarray[double,ndim=1] d = np.ones_like(z, dtype=np.float64, order='C')
          cdef int i, num
          num = len(z)
          for i in range(num):
              zs=z[i]
              if zs < 0:
                 raise ValueError("Redshift z must not be negative")
              if zs < z_ref:
                 raise ValueError("Redshift z must not be smaller than the reference redshift")
              gsl_integration_qags(&F, z_ref+1, zs+1, 0, 1e-7, 1000, w, &result, &error)
              d[i], err = result, error 
              # check for curvature
              rk = (fabs(omc))**0.5
              if (rk*d[i] > 0.01):
                 if omc > 0:
                    d[i] = sinh(rk*d[i])/rk
                 if omc < 0:
                    d[i] = sin(rk*d[i])/rk

          gsl_integration_workspace_free(w) 
          return d/(1.+z)

提前致谢。

Answer 1

自从我在 cython 中开发以来已经有一段时间了，但是如果我没记错的话，我相信您可以按如下方式声明该函数:

ctypedef double* ( * Da)(double* z, double z_ref, int length)

此函数将返回 double 类型的数组，并允许您将 double 数组作为 z 传递。这是一个函数指针，所以可能不完全是您想要的。

ctypedef double* Da(double* z, double z_ref, int length)

这将完成同样的事情，但作为一个常规函数，而不仅仅是一个函数指针。函数和函数指针之间的区别是您必须为函数指针分配一个指向的函数。