4

I want to cythonize a code in python to speed up the code. In the following you can see my attempt to make my python class understandable for cython:

import numpy as np
cimport numpy as np
ctypedef np.double_t DTYPE_T
cpdef double std_G,v_c
std_G=4.3e-9 # Newton's const   in Mpc (km/s)^2 M_sol^{-1}
v_c = 299792.458 #km/s

cdef extern from "math.h":
    double log(double) nogil
    double sqrt(double) nogil


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 do_callback(double x, void* params): 
     return (<ComovingDistMemoization>params).eval(x) 


cdef class ComovingDistMemoization(object):
     cdef list _memotable
     cdef cosmolgy
     def __cinit__(self, cosmology, memotable = None):

         if memotable is None:
            self._memotable = []

         self._memotable = memotable
         self.cosmology = cosmology
     def __call__(self, double z):


        if z in self._memotable:
            return self._memotable[z]

        def eval(z):

            return 1./sqrt(self.cosmology.hubble2(z))
        cdef gsl_integration_workspace* w =gsl_integration_workspace_alloc(1000)
        cdef gsl_function F            

        F.function = &do_callback 
        F.params = <void*>self 
        cdef double result = 3, error = 5
        cdef double y, err, dist
        gsl_integration_qags (&F, 0, z, 0, 1e-7, 1000, w, &result, &error) 
        y, err = result, error 

        gsl_integration_workspace_free(w) 

        dist = self.cosmology.v_c * y  #to get proper units, ie to put in the hubble length

        self._memotable[z] = dist

        return dist


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

        self.omega_m = omega_m
        self.omega_l = omega_l
        self.omega_r = omega_r
        self.h = h
        self.w = w
        self.G = G
        self.v_c = v_c

        self.comovingdist = ComovingDistMemoization(self)
    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

    def hubble2(self, double z):
        cdef double inv_a
        inv_a = 1.+z
        return (self.omega_r*inv_a**4 + self.omega_m*inv_a**3 + \
                  self.omega_l*(inv_a**(3*(1+self.w))) + (1 - self.omega_m - self.omega_l - self.omega_r)*inv_a**2)*self.H0**2

    property hubble_length:
        def __get__(self):
            return self.v_c / self.H0

    def rho_crit(self, double z):
        return 3.*self.hubble2(z)/(8*np.pi*self.G)


    def angulardist(self, double z, double z2 = None):

        if z2 is None:
            return self.comovingdist(z) / (1+z)

        return (self.comovingdist(z2) - self.comovingdist(z)) / (1+z2)

However one of the places that raises error and I couldn't find so far any substitution and my investigation just reached to this point that the __copy__(self) function is not supported by cython. First question: How could I make a copy of one method of a class in cython? I read that with pickle.Pickler, it is possible to provide a substitution for an instance of a class but I don't know how it should produce a copy of a class inside the class?? Update: Is it a right way to replace __copy__ with the following piece of code:

def __reduce__(self):
    return (self.__class__, (), self.__getstate__())
def __getstate__(self):
    return (self.omega_m, self.omega_l, self.omega_r, self.h, self.w, self.G, self.v_c)
def __setstate__(self, data):
    (self.omega_m, self.omega_l, self.omega_r, self.h, self.w, self.G, self.v_c) = data

My second question is about property, have I defined the property in cython or I also need to have a __set__ function as well?

My last question: when I call all the instances of cosmology class using ComovingDistMemoization they return zero. I am wondering whether the way I have called for instance cosmology class in ComovingDistMemoization class and vice versa caused the problem and can not pass info between them ?

Dalek
  • 4,168
  • 11
  • 48
  • 100
  • I know very little about cython, but your current `__copy__` code wouldn't work in regular Python because you're using a capital C when you try to call `Cosmology` while the class name starts with a lowercase c. I'm also not sure why you need the `__copy__` method at all, since you don't seem to be doing anything fancy in it (a regular Python class would function exactly the same without that method defined). – Blckknght Jul 08 '14 at 15:59
  • @Blckknght lol! you are right. I updated it. However I am still wondering whether it can be replaced with `pickle.Pickler` and how? Another point whether `__copy__` is eligible in cython and which one is faster `__copy__` or using pickle? – Dalek Jul 08 '14 at 16:08
  • @Dalek did the answer below solved this question of yours about property? – Saullo G. P. Castro Aug 14 '14 at 12:06
  • @SaulloCastro Unfortunately defining the property like this in cython didn't work! – Dalek Aug 14 '14 at 12:14

2 Answers2

1

Answering your second question, in your case you don't need the __set__ function, and it seems there is an error in the way you are using the __get__ concept. The attribute name to be obtained has to be set like:

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

def get_hubble_length(self):
    return self.v_c / self.H0
hubble_length = property(fget=get_hubble_length)

EDIT: I was not aware that @property can't be used with Cython

Saullo G. P. Castro
  • 56,802
  • 26
  • 179
  • 234
  • so they way I have defined it now is not correct, because I haven't used `__set__` function or there is another mistake? – Dalek Jul 09 '14 at 07:06
  • @Dalek you don't have to define `__set__` or `__get__` when using `property`... in Python you can only put it as a decorator that it will work – Saullo G. P. Castro Jul 09 '14 at 07:08
  • I came across [this example and the answer](http://stackoverflow.com/questions/7666873/cython-and-deepcopy-woes-with-referenced-methods-functions-any-alternative-id) that cython doesn't work properly with decorators. It also decrease the speed of the code, isn't it right? – Dalek Jul 09 '14 at 07:16
  • @Dalek you are right, I was not aware that Cython did not accept `@property`, I've updated the answer... hope it helps... – Saullo G. P. Castro Jul 09 '14 at 07:42
  • Thanks a lot. I think my last question is also not very difficult but I couldn't find any example of it in the web. How could I **calling** a class as an object in another class in cython (especially in my case)? – Dalek Jul 09 '14 at 07:54
  • sorry, mate that I tag you again! I am a bit stuck again with declaration of class as my attributes for another class. I appreciate if you would suggest something how I could solve the problem. Thanks in advance! – Dalek Jul 11 '14 at 16:33
  • @Dalek did this answer solved your problem with `property`? you could post this new problem as another question... – Saullo G. P. Castro Jul 11 '14 at 18:47
  • I tried your reply and got an error message but I tested what I have already written in this post and it works perfectly. Actually I did post another question ans so far no one replied. – Dalek Jul 11 '14 at 19:48
  • @Dalek did you try the approach proposed here (after the edit)? – Saullo G. P. Castro Sep 20 '14 at 12:07
  • 1
    when I use your approach it does not work in `cython`, it works in python but when I define the property by using `__get__`, it works in `cython`. – Dalek Sep 20 '14 at 12:32
1

The best way to go is to combine cosmology class with ComovingDistMemoization class and write the second class as a method for the first one. I renamed this class as comovingdist function for the cosmology class.

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 do_callback(double x, void* params): 
     return (<cosmology>params).eval(x) 

cdef class cosmology(object):
    cdef public double omega_m, omega_l, h, w, omega_r, G, v_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.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

    def hubble2(self, double z):
        cdef double inv_a
        inv_a = 1.+z
        return (self.omega_r*inv_a**4 + self.omega_m*inv_a**3 + \
                  self.omega_l*(inv_a**(3*(1+self.w))) + (1 - self.omega_m - self.omega_l - self.omega_r)*inv_a**2)*self.H0**2

    property hubble_length:
        def __get__(self):
            return self.v_c / self.H0

    cpdef double eval(self, double z):
         return 1./sqrt(self.hubble2(z))

    def comovingdist(self, double z):
        cdef gsl_integration_workspace* w =gsl_integration_workspace_alloc(1000)

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

        cdef double result = 3, error = 5
        cdef double y, err, dist
        gsl_integration_qags(&F, 0, z, 0, 1e-7, 1000, w, &result, &error)
        y, err = result, error 

        gsl_integration_workspace_free(w) 

        dist = self.v_c * y  #to get proper units, ie to put in the hubble length

        return dist


    def rho_crit(self, double z):
        return 3.*self.hubble2(z)/(8*np.pi*self.G)


    def angulardist(self, z, z2=None):
        if z2 is None:
            return self.comovingdist(z) / (1+z)
        return (self.comovingdist(z2) - self.comovingdist(z)) / (1+z2)

    def beta(self, z, zcluster):

        Ds = np.array([self.angulardist(zi) for zi in z])
        Dls = np.array([self.angulardist(zcluster, zi) for zi in z])

        Dls_over_Ds = np.zeros_like(Dls)
        Dls_over_Ds[Ds > 0] = Dls[Ds > 0] / Ds[Ds > 0]
        Dls_over_Ds[Dls <= 0] = 0

        return Dls_over_Ds


    def beta_s(self, z, zcluster):

        betainf = self.beta([1e6], zcluster)

        beta_s = self.beta(z, zcluster) / betainf

        return beta_s


    def RhoCrit_over_SigmaC(self,z, zcluster):

        Dl = self.angulardist(zcluster)
        r= (4*np.pi*self.G)*self.rho_crit(zcluster)*self.beta( z, zcluster )*Dl/self.v_c**2
        return r
Dalek
  • 4,168
  • 11
  • 48
  • 100