Wrap C functions in Python

C_and_Python.md

Pi Calculation in C vs Python

This is another example which we show how we can wrap a function in C to gain more performance. We use the integral of $\int_{0}^{1} \frac{4}{1+x^{2}} dx$ to calculate $\pi$.

This is the C code:

/*
    Code to compute the number Pi = 3.14

    Integral 4/(1+x^2) from 0 to 1 is 3.14

    gcc -shared -o pi.so -fPIC pi.c
*/
#include <stdio.h>
#include <time.h>


float pi(int num_steps);

float pi(int num_steps)
{
     double step, x, sum;
     
     step = 1.0/(double)num_steps;
     sum = 0;
     for (int i=0;i<num_steps;i++)
     {
          x = (i+0.5)*step;
          sum = sum + 4.0/(1.0+x*x);
     }
     return step * sum;
}

int main()
{
    float result;
    double start, end;
    start = clock();
    result = pi(10000000);
    end = clock();

    printf("Pi = %f in %f secs\n", result, (float)(end-start)/CLOCKS_PER_SEC);
    return 0;
}

You need to compile this C function with $> gcc -shared -o pi.so -fPIC pi.c. Then you can use the following python wrapper.

"""
    Compute Pi with a C function in Python

    1. First, compile a shared library
        $> gcc -shared -o pi.so -fPIC pi.c

    2. Then, execute this script with
        $> python Pi.py
"""
import ctypes
import numpy as np


def pypi(npoints):
    ''' Python implementation of Pi calculator
    '''
    step = 1.0/npoints
    sum = 0

    for i in range(npoints):
        x = (i+0.5)*step
        sum = sum + 4.0/(1.0+x*x)

    return sum * step


class Pi:

    def __init__(self):
        self.lib = ctypes.CDLL(r'/Users/mehdi/github/ThirtyDaysOfC/pi.so')
        self.lib.pi.restype = ctypes.c_float
        self.lib.pi.argtypes = (ctypes.c_int, ) #_argtypes_ must be a sequence of types

    def __call__(self, npoints):
        return self.lib.pi(ctypes.c_int(npoints))


if __name__ == '__main__':
    from time import time

    # initialize the C Object
    pi = Pi()

    times = {'C':[],
            'Python':[]}

    np_grid = np.logspace(1, 7)
    for np_float in np_grid:

        npoints = int(np_float)

        t0 = time()
        pi_c = pi(npoints)
        t1 = time()
        pi_py = pypi(npoints)
        t2 = time()

        print((abs(pi_c - np.pi) < 1.0e6), (abs(pi_c - pi_py) < 1.0e6))

        #print('Pi (C) = %f in %f secs'%(pi_c, t1-t0))
        #print('Pi (Python) = %f in %f secs'%(pi_py, t2-t1))
        #print('C implmentation is %.1f times faster than Python'%((t2-t1)/(t1-t0)))

        times['C'].append(t1-t0)
        times['Python'].append(t2-t1)

    import matplotlib
    matplotlib.use('TKAgg')
    import matplotlib.pyplot as plt

    plt.figure()
    for language in ['C', 'Python']:
        plt.plot(np_grid, times[language], marker='o', label=language)

    plt.xscale('log')
    plt.yscale('log')
    plt.grid(ls=':', color='grey', alpha=0.5)

    plt.tick_params(direction='in', right=True,
                    top=True, which='both', axis='both')

    plt.xlabel('Nsteps', fontsize=14)
    plt.ylabel('time [sec]', fontsize=14)

    plt.legend(fontsize=14)
    #plt.show()
    plt.savefig('pi_c_vs_python.png',
                dpi=300, bbox_inches='tight')

mandelbrot.md

Mandelbrot in C

This note demonstrates how we can compile a code written in C for Python. The example is based on the IPython Interactive Computing and Visualization Cookbook, Second Edition (2018), by Cyrille Rossant

• Step 1: have a C function. In this example, we cover Mandelbrot, the C code, mandelbrot.c, is:

#include "stdio.h"
#include "stdlib.h"

// function declaration
void mandelbrot(int size, int iterations, int *col);

void mandelbrot(int size, int iterations, int *col)
{
    // Variable declarations.
    int i, j, n, index;
    double cx, cy;
    double z0, z1, z0_tmp, z0_2, z1_2;

    // Loop within the grid.
    for (i = 0; i < size; i++)
    {
        cy = -1.5 + (double)i / size * 3;
        for (j = 0; j < size; j++)
        {
            // We initialize the loop of the system.
            cx = -2.0 + (double)j / size * 3;
            index = i * size + j;
            // Let's run the system.
            z0 = 0.0;
            z1 = 0.0;
            for (n = 0; n < iterations; n++)
            {
                z0_2 = z0 * z0;
                z1_2 = z1 * z1;
                if (z0_2 + z1_2 <= 100)
                {
                    // Update the system.
                    z0_tmp = z0_2 - z1_2 + cx;
                    z1 = 2 * z0 * z1 + cy;
                    z0 = z0_tmp;
                    col[index] = n;
                }
                else
                {
                    break;
                }
            }
        }
    }
}

• Step 2: compile the C code into a shared library with $> gcc -shared -Wl,-soname,mandelbrot -o mandelbrot.so -fPIC mandelbrot.c

• Step 3: write a wrapper object, MandelBrot defined in Mandelbrot.py to access the library inside Python with Ctypes:

import ctypes
import numpy as np
from numpy.ctypeslib import ndpointer

class MandelBrot:
    
    def __init__(self):
        
        # path to the shared library
        self.lib = ctypes.CDLL(r'/home/medirz90/Downloads/github/ThirtyDaysOfC/mandelbrot.so')
        
        self.lib.mandelbrot.restype = None
        self.lib.mandelbrot.argtypes = (ctypes.c_int, 
                                   ctypes.c_int,
                                   ndpointer(ctypes.c_int))
    
    def __call__(self, size, iterations):
        
        # placeholder for the output
        col = np.empty((size, size), dtype=np.int32)
        
        self.lib.mandelbrot(ctypes.c_int(size),
                            ctypes.c_int(iterations),
                            col)
        return col

• Step 4: Finally, you can import the MandelBrot object in Python with >>> from Mandelbrot import MandelBrot. You can run the following example, test_mandelbrot.py with $> python test_mandelbrot.py:

import numpy as np
import matplotlib.pyplot as plt

from time import time
from Mandelbrot import MandelBrot

def test():
    t0 = time()

    size = 400
    iterations = 100

    # call C
    MB = MandelBrot()
    col = MB(size, iterations)

    print('done in {:.2f} sec'.format(time()-t0))

    # plot
    fig, ax = plt.subplots(1, 1, figsize=(10, 10))
    ax.imshow(np.log(col), cmap=plt.cm.hot)
    ax.set_axis_off()
    plt.savefig('mandelbrot.png')
    
    return 0

if __name__ == '__main__':
    test()

Then, you will get something like: