Define a pointer as a method of a class?

Question

I am trying to speed up my cython code. I came across this link where the author has described how using pointers instead of numpy arrays can improve the speed of cython codes. In my cosmology class the bottleneck is Da function. I am not very familiar with pointers in C, I would appreciate if somebody give me an idea:

Is it possible to define a method of a class as a pointer for instance in my case convert np.ndarray[double, ndim=1] Da to something like 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)

Thanks in advance.

Answer 1

It has been a while since I developed in cython, but if memory serves me I believe you could declare the function as follows:

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

This function will return an array of type double and allow you to pass an array of doubles in as z. This is a function pointer, so maybe not quite what you want.

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

this will accomplish same thing but as a regular function, not just a function pointer. Difference between function and function pointer is you have to assign a function pointer a function to point to.