How does one directly call a ctypes SciPy integrand?

Question

I have some numerical code to run in SciPy. It involves complicated rational functions of exponentials of polynomials, so it's pretty computationally expensive; ergo, I have written it in C and am calling it using ctypes. The most important use case is as an integrand for scipy.integrate.quad, but I also occasionally need to call it directly.

The "natural" signature of my function is

double f(double x, double y, double z){}

which the ctypes documentation suggests using, with the corresponding Python

import ctypes
so = ctypes.CDLL('/path/to/f.so')
f = so.f
f.restype = ctypes.c_double
f.argtypes = (ctypes.c_double, ctypes.c_double, ctypes.c_double)

Regardless, however, to call it as an integrand, SciPy requires a certain function signature

double f(int n, double args[n]){}

which is specified in the Python code as

import ctypes
so = ctypes.CDLL('/path/to/f.so')
f = so.f
f.restype = ctypes.c_double
f.argtypes = (ctypes.c_int, ctypes.c_double)

To pass the parameters y and z when performing the integral, they're passed to quad as a tuple called args.

scipy.integrate.quad(f, lowerlimit, upperlimit, args=(y,z))

This makes it unclear how to call f directly. My naive attempt was

f(3, (x,y,z))

but this yields a type error for argument 2. Several variations have similarly failed to work. It's not really surprising; ctypes is expecting a function call with exactly one integer argument followed by exactly one double argument.

I have absolutely no idea how quad gets y and z into f. I tried looking into the SciPy source, but I must admit I got lost trying to trace the calls from the Python into the C into the Fortran.

I could just write another function as a wrapper for either the direct calls or the integration, but it seems more elegant to use just one form, and at the very least, I'd like to understand how the SciPy call works.

How can I call the SciPy integrand form of f directly, passing all three parameters x, y, and z?

I'm using Python 3.4.3, NumPy 1.11.2, and SciPy 0.18.1.

Edit: Note that f can be called by changing its argtypes:

f.argtypes = (ctypes.c_int, 3*ctypes.c_double)
f(3, (3*ctypes.c_double)(x, y, z))

I'm still curious what it is that SciPy is doing, though. Setting the argtypes back and forth all the time is still inelegant and inconvenient, at best.

Edit 2: Note that this question, after the previous edit, is now basically a duplicate of this one, which has helpfully popped up in the right column.

score 2 · Answer 1 · answered Mar 07 '17 at 21:13

2

You cannot fix this at the Python level. The documentation https://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#faster-integration-using-ctypes states

Write an integrand function in C with the function signature double f(int n, double args[n])

In your situation, you must thus add a function at the C level that is

double f_for_scipy(int n, double args[n]) {
    return f(args[0], args[1], args[2]);
}

and feed that to quad instead.

answered Mar 07 '17 at 21:13

Pierre de Buyl

7,074
2
16
22

"You cannot fix this at the Python level." Why not? SciPy has to call the function, with all its parameters, somehow. – calavicci Mar 07 '17 at 21:49
1

You cannot fix it while keeping the compiled performance of quad. When the integrand is a compiled C function that is loaded with ctypes and passed to quad, the whole procedure is done with compiled code. If you pass a Python callable, this would also work but the performance would be less interesting. You should then use the "real" interface `f.argtypes = (ctypes.c_double, ctypes.c_double, ctypes.c_double)` and call `scipy.integrate.quad(lambda x, y, z: f(x, y, z), lowerlimit, upperlimit, args=(y,z))`. This means Python will convert the `x, y, z` argument properly at a cost. – Pierre de Buyl Mar 07 '17 at 21:59
I did find a way to do it in Python with no performance cost; see above edit. It's not pretty, though, and I'm still curious about just how SciPy is going about invoking the function. – calavicci Mar 07 '17 at 22:06
SciPy is "just" passing the function pointer up to the compiled quad routine. I am curious to know why your solution works, do you have a pointer to the ctypes doc for this? Do you have any method to check that the numerical result is correct? – Pierre de Buyl Mar 07 '17 at 22:18

score 0 · Answer 2 · edited May 23 '17 at 12:10

I don't know if this helps with the ctypes case, but when calling a Python integrad these scipy functions concatenate the free variable with the args. In other words

def bar(x,y,z):
   return np.sin(x*y*z)
In [43]: quad(bar,0,1, args=(1,2))
Out[43]: (0.7080734182735712, 7.861194120923578e-15)

When evaluating at 0.5, it does, (x,)+args:

In [49]: bar(*(.5,)+(1,2))
Out[49]: 0.8414709848078965

In [50]: bar(.5,1,2)
Out[50]: 0.8414709848078965

So with the c signature:

f_for_scipy(int n, double args[n])

n is number of arguments, and args[n] is the pointer that array of values.

I'm experimenting these sorts of calls using cython and its extension types

https://cython.readthedocs.io/en/latest/src/tutorial/cdef_classes.html

Passing a cython function vs a cython method to scipy.integrate

score 0 · Answer 3 · answered Mar 09 '17 at 17:09

You need to have the correct function signature, but type Python argtypes should be POINTER(c_double). Arrays in C decay to pointers in function arguments:

C Example (Windows)

#include <stdio.h>
__declspec(dllexport) double f(int n, double args[])
{
    double sum = 0;
    int nn;
    for(nn = 0; nn < n; ++nn)
    {
        printf("args[%d] = %ld\n",nn,args[nn]);
        sum += args[nn];
    }
    return sum;
}

ctypes Example

>>> from ctypes import *
>>> dll = CDLL('x')
>>> dll.f.restype = c_double
>>> dll.f.argtypes = c_int,POINTER(c_double)
>>> L=[1.1,2.2,3.3,4.4,5.5]
>>> args = (c_double*len(L))(*L)
>>> dll.f(len(args),args)
args[0] = 1.100000
args[1] = 2.200000
args[2] = 3.300000
args[3] = 4.400000
args[4] = 5.500000
16.5

How does one directly call a ctypes SciPy integrand?

3 Answers3