Calling C from within Python, Part 2: Calling and Callbacks

Introduction

In the previous blog post, we have written a Runge-Kutta fourth order solver in plain C, which can be used to numerically solve systems of differential equations. We have also compiled it as a shared object in order to make it available from within other languages such as Python. In this post, we will do exactly this, using the ctypes library, which is a "foreign function library for Python", which "provides C compatible data types, and allows calling functions in DLLs or shared libraries."

Loading the Shared Object with ctypes

First, we have to load the ctypes library which allows us to interact with the shared object:

# Import ctypes
from ctypes import *

# Load the shared object
libname = "./rk4.so"
c_lib = CDLL(libname)

Setting up the Function

Next, we get the function from our shared object, and provide annotations for the return type, as well as the types of the passed function arguments. From my (limited) experience, annotating the argument types does not seem to be strictly required, as ctypes seems to be able to infer them by itself just fine. In contrast, however, the return type is always implicitly assumed to be a simple int. In our case, the solver will return a double**, so we definitely have to change the annotation manually. As it is (in my opinion) good practice, we will also provide explicit annotations for the argument types.

# The function from the SO we want to call
solve = c_lib.solve                                                     # Get the function
solve.restype = POINTER(POINTER(c_double))                              # Annotate return type (double**)
solve.argtypes = [                                                      # Annotate argument types:
    CFUNCTYPE(None, c_double, POINTER(c_double), POINTER(c_double)),    # - void (*f) (double, double*, double*)
    POINTER(c_double),                                                  # - double*
    c_double,                                                           # - double
    c_double,                                                           # - double
    c_int]                                                              # - int

Preparing the Parameters

Now, we prepare the parameters which we want to pass to the function. For now, we ignore the function pointer and just focus on all the "simple" parameters.

# The parameters we want to pass
y0 = (c_double*3)(1,2,3)            # An array y0 of three doubles (1.0, 2.0, 3.0)...
y0 = cast(y0, POINTER(c_double))    # ... which we cast to double* y0
dt = c_double(0.0001)               # double dt = 0.0001
tmax = c_double(10.0)               # double tmax = 10.0
d = c_int(3)                        # int d = 3

What's now missing is only the first parameter, which is a pointer to the function we want to numerically solve. It has the signature void (*f)(double, double*, double*). If we want to pass a function defined in Python (which sometimes is called a "callback"), we have to cast this function to a C function type with the required argument and return types. We do so by using the followig convenient decorator, in which the first element of the tuple encodes the function's return type (None is used for void), and all the subsequent elements corespond to the types of the function arguments.

#Our Python function to be called from within C
@CFUNCTYPE(None, c_double, POINTER(c_double), POINTER(c_double))
def lorenz(t, y, dydt):
    dydt[0] = 10.0*(y[1]-y[0])
    dydt[1] = y[0]*(28.0-y[2]) - y[1]
    dydt[2] = y[0]*y[1] - (8.0/3.0)*y[2]

Now we have prepared all parameters and are ready to call the solver.

Calling the Function and Working with its Output

The call to the function itself is now extremely straight-forward. All we have to do is:

ret = solve(lorenz, y0, dt, tmax, d)

Now, we can do whatever we want with the returned results. Accessing the contents the double** points to is straightforward using the square-bracket syntax:

# Do stuff with returned data
print(ret[0][0:100])        # Print the first 100 timepoints
print(ret[1][0:100])        # Print the first 100 values of x(t)

Cleanup!

Finally, when we're done we need to make sure to free all memory we have allocated. For this, in the last blog post, we have made sure to provide a method in our shared object that takes care of this. We just need to call it via the ctypes library:

# Proper deallocaton of memory in C from within Python
dealloc = c_lib.dealloc                             # The method to call
dealloc.argtypes = [POINTER(POINTER(c_double)),     # Its argument types:
                    c_int]                          # - (double**, int)
dealloc.restype = None                              # Its return type (void)
dealloc(ret, 3)                                     # Calling it

It is important to only call this method once, otherwise we'd cause a double free, causing the abortion of the software.

Outlook

One thing has become immediately obvious: While it is generally easy and straghtforward to interact with a SO from within Python, there are still some little details to be taken care of, in particular the type annotations, the preparation of the function parameters and the proper deallocation of memory. These things could easily be packaged away into a little Python library for the user to interact with, which takes care of all of these things "under the hood." Thus, this would be an immediately obvious next step for this project, which -- for completeness sake -- I will likely implement in a future blog post.

Second, this little project taught me a more general important thing: It generally is a really good idea to use every (personal) project to learn a completely new skill. This not only naturally expands one's skill set, but also makes every project significantly more exciting and enjoyable. Thus, I'll definitely continue doing so in the future :-).

Update (16th October 2022)

The little "convenience library" I mentioned in the previous section is now availabe on my github. Check it out!