Combining C++ and Python: dynamic library

The simplest way to combine the unprecedented graphical potential of Python (matplotlib library) and a speed of numerical calculations of pure C/C++ languages is to compile the C-code as a dynamic library.

Here I give an example.

#include <cmath> // for modf and cos 
using namespace std;
extern "C" double value (double);  
double value (double step) {
double res, pos, fractpart, intpart;
int n;
fractpart = modf(M_PI / 2.0 / step, &intpart); // Compute a number of steps for numerical integration
n = intpart;
res = 0.0;
for ( int i=0; i < n; i++) {
        pos = i * step;                 // Eulear integration method
        res += step * cos(pos);
}
res += fractpart*step * cos(M_PI / 2.0);
return res;  
}

extern “C” is necessary to clarify which functions are visible outside the dynamical library. The python package ‘ctype’ was written for pure C-language. Therefore, the class structure is not inherited. The class structure should be duplicated in the Python code.

To compile the dynamical library we use following command (iOS):

g++ -dynamiclib -flat_namespace test_func.cpp -o test_fucn.so

For Linux the arguments slightly differ:

g++ -fPIC -shared test_func.cpp -o test_func.so

The Python script has following structure:

import numpy as np                           ## library to work with arrays
import matplotlib.pyplot as plt              ## plotting library
import ctypes
from ctypes import cdll, c_double
lib = cdll.LoadLibrary('./test_func.so')    ## use our library written in C
grid = np.linspace (0.01, 1.0, 50)           ## prepare a grid with stepsizes
res_list=[]
lib.value.restype = ctypes.c_double          ## Declare a type of return argument (important)
for i in range (0, len(grid)):               ## for each step size we call our function
        step = c_double(grid[i])
        res = lib.value(step)
        res_list.append(res)
plt.plot(grid, res_list)                        ## plot the resylt of numerical integration
plt.xlabel('Numerical integration step size')
plt.ylabel('Result of integration')
plt.savefig('plot.pdf')
plt.show()

The only complicated part of the code is how to return values.

If the code works correctly after we call Python scripts, it prepares following figure:

When the integration step is tiny (~0.1) we get the right result. If the integration step is large, the result fluctuates a lot.