Passing arguments to scipy.LowLevelCallable while using functions in C

Question

I am trying to use C defined function for numerical integration in SciPy. The example case given here (SciPy documentation) works fine.

In my case, the testlib.c file is

/* testlib.c */

#include <math.h>
#define PI 3.14159265358979323846

double factor(double phi, double r) {
    double val = (2*PI)/(pow(r, 2) + 3*cos(phi));
    return val;
}

//------------------------------------

double f2(int n, double *x, void *user_data) {
    double c = *(double *)user_data;
    double v1 = factor(c, 0.25); // value of phi defined inline but it is an argument
    return v1 + x[0] - x[1] ; /* corresponds to v1 + x - y    */
}

And, the test.py function calling the testlib.so file obtained after compilation is below,

import os, ctypes
from scipy import integrate, LowLevelCallable

lib = ctypes.CDLL(os.path.abspath('testlib.so'))

# define return type in .restype
lib.f2.restype = ctypes.c_double

# define argument type in .argtypes
lib.f2.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)

# additional argument, here a constant, casting needed
c = ctypes.c_double(1.0)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)

# pass extra argument
func = LowLevelCallable(lib.f2, user_data)

# quadrature in 3-dimensions
out=integrate.nquad(func, [[0, 10], [-10, 0]])

print(out)

# -----------------------------------------------
phi = 0.25   #  How to pass these to the C function 
esv = 1.25   
cv1 = 0.03   
cv2 = -0.15  
cv3 = 3.0

My question: How to pass additional arguments such as c to the function f2. In my case I have 5 such arguments, available as np.float64 in the calling py file.

I wonder if I can pass the arguments as array as user_data to the function f2 .

From documentation for nquad, found that arguments are to be passed as array and int n in the C function is the number of arguments passed.
Also, I'm Open to try other options such as cython, pyCapsule but no experience there in. Found very similar question using numba and jit, where no additional arguments are passed. Using numba and jit for integration: SE

For compiling testlib.c: $ gcc -shared -fPIC -o testlib.so testlib.c

Btw, you probably should compile testlib.so with optimization on (i.e `-O2` or `-O3` or whatever best) to get best performance. — ead, May 05 '20 at 08:02
As long as you only pass scalars dtype np.float64 There is an example here https://stackoverflow.com/a/60619681/4045774, which shows how to do this using Numba. You can pass the additional params as args, as you would do in pure Python. — max9111, May 05 '20 at 12:07

ead · Accepted Answer · 2020-05-05T10:00:43.983

They are more than one way to skin the cat, but if your are using ctypes there is possibility to stick with ctypes, for example:

You can create an array an initialize it with values, e.g.:

ptr_to_buffer=(ctypes.c_double*5)(phi,esv,cv1,cv2,cv3)
user_data = ctypes.cast(ptr_to_buffer, ctypes.c_void_p)

or if the data is already in a numpy-array (as I originally understood your question):

import numpy as np
a=np.array([1.0, 2.0, 3.0, 4.0, 5.0], np.float64)

import ctypes
ptr_to_buffer=(ctypes.c_double*5).from_buffer(a)
user_data = ctypes.cast(ptr_to_buffer, ctypes.c_void_p)

in this case user_data is a copy of a and doesn't share the memory, which is sometimes a good thing and sometimes not.

For bigger arrays one can also let user_data also to share the memory with the numpy-array:

user_data = ctypes.c_void_p(a.__array_interface__['data'][0])

as can be verified via:

ctypes.cast(user_data, ctypes.POINTER(ctypes.c_double))[0] = 42.0
print(a[0])
# 42.0 and not 1.0

For this variant you actually have to check, that the memory of the numpy-array is contiguous, for an example how this information can be obtained cab be looked up for example in numpy.ctypeslib.as_ctypes.

Maybe a less low-level way of obtaining the pointer is

user_data =  ctypes.cast(np.ctypeslib.as_ctypes(a), ctypes.c_void_p)

but still checks of the shape/strides are needed.

Passing arguments to scipy.LowLevelCallable while using functions in C

1 Answers1