1

This question builds upon the one here

I have the following cython code that imports a file (cGslInteg.pxd) containing the cython declarations for the relevant parts of the gsl/gsl_integration.h and gsl/gsl_math.h header files and then defines a class to compute an integral:

from libc.math cimport log
from libc.math cimport sqrt
from cGslInteg cimport *

ctypedef double * double_ptr

cdef double integrand(double x, void * params):
    """integrand implemented outside class """

    cdef double alpha = (<double_ptr> params)[0]
    cdef double f = log(alpha*x) / sqrt(x)
    return f


cdef class CalcSomething(object):

    def integrate(self, double a, double b, double alpha):
        """This does the integral"""

        cdef double result, error

        cdef gsl_integration_workspace * W
        W = gsl_integration_workspace_alloc(1000)

        cdef gsl_function F

        F.function = &integrand          # works 
        # F.function = &self._integrand  # doesn't work

        cdef double params[1]
        params[0] = alpha
        F.params = &params

        cdef double epsabs = 0.
        cdef double epsrel = 1e-7
        gsl_integration_qags(&F, a, b, epsabs, epsrel, 1000, W, &result, &error)
        gsl_integration_workspace_free(W)

        return result


    cdef double _integrand(self, double x, void * params):
        """integrand implemented inside class """

        cdef double alpha = (<double_ptr> params)[0]
        cdef double f = log(alpha*x) / sqrt(x)
        return f

The above code compiles and runs correctly as the integrand function outside the class is used (F.function = &integrand). However, changing this line to the one commented out beneath (F.function = &self._integrand) in order to use instead the integrand function inside the class, brings the following compilation error:

cython_class_gsl.pyx:31:21: Cannot assign type 'double (*)(CalcSomething, double, void *)' to 'double (*)(double, void *)'

Which makes sense because my declaration for gsl_function in cGslInteg.pxd is:

cdef extern from "gsl/gsl_math.h":

    # Definition of an arbitrary function with parameters
    ctypedef struct gsl_function:
        double (* function) (double x, void * params)
        void * params

My question is: can I re-declare gsl_function so it expects something with type e.g. double (*)(PythonObject, double, void *) or can I wrap self._integrand so it appears to have type double (*)( double, void *)?

EDIT: Some extra notes to make this closer to the real-life problem/issue:

  1. Any wrapper of self._integrand would have to be defined within the CalcSomething class.
  2. The integration parameter alpha cannot be a class variable (i.e. it must be passed within the params argument to the integrand)
Community
  • 1
  • 1
alexabate
  • 143
  • 2
  • 10
  • I've edited my answer to address note 2. I don't understand why note 1 "Any wrapper of self._integrand would have to be defined within the CalcSomething class." is a requirement? What stops you from declaring the wrapper function outside the class? – DavidW Oct 07 '15 at 06:28
  • It is probably just an implementation choice, but I'm interested to understand if it is possible. As was described (maybe poorly/unclearly) in my original question, the code **does** compile and run correctly when the integrand function is outside the class, so there's no need to wrap in this case. – alexabate Oct 07 '15 at 18:17
  • I think "not easily". Somehow ctypes finds a way (see https://docs.python.org/2/library/ctypes.html#callback-functions) that does manage to do something very similar. (I think cffi can do this too). How it does this looks __very__ involved but you might be able to mix Cython and ctypes a little to accomplish what you want with a bit of effort. – DavidW Oct 07 '15 at 20:14
  • I think I will leave finding a solution to note 1. for now: sounds like it will be more of a mess than having everything self-contained within the one class – alexabate Oct 07 '15 at 23:20

1 Answers1

1

This is what params is for. You'd do something like this (hugely cut down example...)

from libc.math cimport log,sqrt

cdef class CalcSomething:   
    cdef double integrand(self, double x, double alpha):
        return log(alpha*x)/sqrt(x)

cdef double wrapped_integrand_func(double x, void* params):
    cdef object params_as_object = <object>params
    cdef CalcSomething instance = params_as_object[0] # we need to know the type to use cdef functions
    alpha = params_as_object[1]
    return instance.integrand(x,alpha)

def test(alpha):
    cdef object o = (CalcSomething(),alpha) # store Calc something and alpha as a tuple
    cdef void* o_ptr = <void*>o # be careful with reference counts here - o_ptr doesn't keep o alive!
    print(wrapped_integrand_func(4.3,o_ptr))

This obviously isn't a general solution to the problem is wrapping a bound method but in this case it's exactly why gsl lets you pass a void*.

The one caveat is that is that you have to ensure that o remains alive for as long as you want to use o_ptr.

Update (to match slightly edited question):

The integration parameter alpha cannot be a class variable (i.e. it must be passed within the params argument to the integrand)

This is fairly easily done. In my original code I used a class variable, but I've now passed them both in separately as a tuple.

DavidW
  • 29,336
  • 6
  • 55
  • 86
  • Maybe I've done something silly but, in the case where passing `alpha` I change my `integrate` function to have the following lines: `cdef object o = (self, alpha); cdef void* p = o; cdef gsl_function F; F.function = &wrapped_integrand_func; F.params = p;` but I get a runtime error of `cython_class_gsl.CalcSomething' object has no attribute 'integrand'" in 'cython_class_gsl.wrapped_integrand_func` – alexabate Oct 07 '15 at 23:13
  • Can't work out why `wrapped_integrand_func` can't see the `integrand` method in this case? (They are `cdef` in both) – alexabate Oct 07 '15 at 23:16
  • Oh defining `integrand` as a `cpdef` type works, but unsure why `cdef` didn't work after the minor update to `wrapped_integrand_func` to include `alpha` in the `void * params` – alexabate Oct 07 '15 at 23:17
  • That's my fault (I didn't test my edit properly...). In order to use `cdef` functions Cython needs to be told the type of the object (it can't find it an runtime like Python). I've fixed that now (`cdef CalcSomething instance = params_as_object[0]`) – DavidW Oct 08 '15 at 07:08