[Python] How to efficiently wrap GSL functions with structs and pointers using Numba?

I'm trying to wrap some GSL (GNU Scientific Library) functions using Numba. I'd like to avoid C glue code wrappers and to be able to cache and extend the wrapped functions.

Here’s a simplified structure of the site-package:

numba_gsl/
├── numba_gsl/
│ ├── gsl_integration.c # Minimal C wrapper for GSL integration
│ ├── integration.py # Python-side ctypes + numba wrapper
├── setup.py # build script


gsl_integration.c exposes a GSL function:

// gsl_integration.c (must be compiled)
#include <gsl/gsl_integration.h>
#include <stdint.h>

typedef double (*func_ptr)(double x, void* params);

double qag(
func_ptr f,
void* user_data,
int which,
double a,
double b,
double epsabs,
double epsrel,
int limit,
int key) {
gsl_integration_workspace* w = gsl_integration_workspace_alloc(limit);
double result, error;

gsl_function gsl_func;
gsl_func.function = f;
gsl_func.params = user_data;

gsl_integration_qag(&gsl_func, a, b, epsabs, epsrel, limit, key, w, &result, &error);
gsl_integration_workspace_free(w);

return result;
}


In integration.py I load the compiled shared library with ctypes, define argument types, and expose a jitted function that calls into it by passing function pointers obtained via Numba's cfunc.

# integration.py
import ctypes as ct
from numba import njit

def get_extension_path(lib_name: str) -> str:
search_path = Path(__file__).parent
pattern = f"*{lib_name}.*"
matches = search_path.glob(pattern)
try:
return str(next(matches))
except StopIteration:
return None

# Load the shared GSL integration library (.so)
_lib_path = get_extension_path('gsl_integration')
_lib = ct.CDLL(_lib_path)

# Define ctypes function prototype for GSL wrapper:
# double qag(void *f, void *params, double a, double b, double epsabs, double epsrel, int limit, int key)
qag_c = _lib.qag
qag_c.argtypes = [ct.c_void_p, ct.c_void_p,
ct.c_double, ct.c_double,
ct.c_double, ct.c_double,
ct.c_int, ct.c_int]
qag_c.restype = ct.c_double

@njit
def qag(func_ptr: int,
a: float,
b: float,
epsabs: float = 1.49e-08,
epsrel: float = 1.49e-08,
limit: int = 50,
key: int = 1,
params_ptr: int = 0) -> float:
"""GSL integration wrapper with finite intervalls."""
return qag_c(func_ptr, params_ptr, a, b, epsabs, epsrel, limit, key)


Here is a build-script:

# setup.py
from setuptools import setup, Extension, find_packages

ext = Extension(
"numba_gsl.gsl_integration",
sources=["numba_gsl/gsl_integration.c"],
libraries=["gsl", "gslcblas"],
extra_compile_args=["-O3"]
)

setup(
name="numba_gsl",
version="0.1.0",
description="Numba-friendly wrappers for GSL routines",
author="No Name",
packages=find_packages(),
ext_modules=[ext],
install_requires=["numba", "llvmlite"],
python_requires=">=3.12",
)


And an example:

# example.py
import math
from numba import cfunc, types
from numba_gsl.integration import qag
from scipy.integrate import quad

# Define the integrand as a Numba cfunc with the proper signature:
@cfunc(types.float64(types.float64, types.voidptr))
def sin_over_x(x, _):
return math.sin(x) / x if x != 0.0 else 1.0

def py_func(x):
return math.sin(x) / x if x != 0 else 1.0

func_ptr = sin_over_x.address
qag_res = qag(func_ptr, a=1e-8, b=3.14)
scipy_res = quad(py_func, a=1e-8, b=3.14)[0]

print("numba_gsl quad result:", qag_res)
print("scipy quad result:", scipy_res)
# numba_gsl quad result: 1.8519366381423115
# scipy quad result: 1.8519366381423115


Is there a (better) way to wrap complex GSL structs and pass them along with function pointers to GSL routines with Numba, perhaps using:

• llvmlite.binding.load_library_permanently
• numba.types.ExternalFunction
• numba.types.Record.make_c_struct

Continue reading...