1

I would like to integrate locally defined python functions with the gsl libraries. To do that, i have implemented the following code with Cython (example with Gauss-Legendre quadrature) :

pxd file :

cdef extern from "gsl/gsl_math.h":
ctypedef struct gsl_function:
double (* function) (double x, void * params) nogil
void * params

cdef extern from "gsl/gsl_integration.h":
gsl_integration_glfixed_table * gsl_integration_glfixed_table_alloc(size_t n) nogil 
double gsl_integration_glfixed(gsl_function *f, double a, double b, gsl_integration_glfixed_table * t) nogil
void  gsl_integration_glfixed_table_free(gsl_integration_glfixed_table *t) nogil 

cdef double int_gsl_GaussLegendre(double func(double, void *) nogil, void * p, double xmin, double xmax) nogil

and the pyx file :

cdef size_t size_GL=1000

cdef double int_gsl_GaussLegendre(double func(double, void *) nogil, void * p, double xmin, double xmax) nogil:
cdef double result, error;
cdef gsl_integration_glfixed_table * W 
cdef gsl_function F
W = gsl_integration_glfixed_table_alloc(size_GL)
F.function = func
F.params = p
result = gsl_integration_glfixed(&F, xmin, xmax, W)    
gsl_integration_glfixed_table_free(W)
return result

This code work for any C function declared within my Cython code. Of course, this will fail when i pass as argument a python function. My script in python :

def gsl_integral(py_func, xmin, xmax, args=()):
cdef size_t sizep = <int>(args.size)
cdef double[:] params = np.empty(sizep, dtype=np.double)
for i in range(0,sizep):
    params[i]=args[i]
cdef gsl_function F
F.function = py_func
F.params = params

which return : "Cannot convert Python object to 'double (*)(double, void *) nogil'"

If i use instead :

def gsl_integral(py_func, xmin, xmax, args=()):
cdef size_t sizep = <int>(args.size)
cdef double[:] params = np.empty(sizep, dtype=np.double)
for i in range(0,sizep):
    params[i]=args[i]
cdef gsl_function F
F.function = <double *>py_func
F.params = params

which return: "Python objects cannot be cast to pointers of primitive types"

I have seen that i could wrap my cython function into a class (Pass c function as argument to python function in cython), but i'm not quite sure to understand how to do it in this situation (plus the example doesn't work me.) As a workaround, i have been passing to Cython two array, x and f(x), the latter estimated with my local python f, in order to define a gsl spline that i can latter integrate, but this not elegant at all. Is there any other ways? I would like to use GSL integration without the GIL

Many thanks, Romain

3
  • Remember that for Python to see your modules you need an entry point (function) with a cpdef declaration which can talk to C/C++ as well as Python. Easy fix... Just write a function with a cpdef that does the C work and returns the Python items.
    – Matt
    Mar 6, 2022 at 19:45
  • Hello Matt, thanks for your answer! I have badly explain myself : I want to integrate purely python functions in .py file by importing the cython code. I'm going try David answers. Cheers, Romain
    – rpaviot
    Mar 8, 2022 at 15:12
  • There is already a GSL package you can install with pip install pygsl -> may want to check that out instead of wrapping functions yourself (although that can be rewarding too) :)
    – Matt
    Mar 8, 2022 at 16:17

1 Answer 1

0

There's a few fundamental issues with your basic premise:

  • C function pointers do not have any state. Python callable objects do have state. Therefore, C function pointers do not have space to store the information required to call a Python function.
  • If you're calling a Python callable object, you cannot release the GIL, since you must have the GIL to execute that object.

This answer outlines your options (the C++ std::function option doesn't apply here). Essentially you can either

  • use ctypes to generate a function pointer from the Python function (it does this with hacky code-generation at runtime - it isn't possible in standard C),
  • or write a cdef function with an appropriate signature to pass as the function pointer, and pass your Python callable to that as part of the void* params.

I recommend the latter, but neither is a great options.

1
  • Hello David, i'm going to try both method to see which one is faster (and safer). Many thanks ! Cheers, Romain
    – rpaviot
    Mar 8, 2022 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.