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

rpaviot
  • 11
  • 2
  • 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 06 '22 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 08 '22 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 08 '22 at 16:17

1 Answers1

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.

DavidW
  • 29,336
  • 6
  • 55
  • 86
  • Hello David, i'm going to try both method to see which one is faster (and safer). Many thanks ! Cheers, Romain – rpaviot Mar 08 '22 at 15:13