Python is particular suitable for rapid prototyping and the testing the suitability of algorithms. In scientific computing, one often encounters two problems: Python is too slow. One needs to use available code and external libraries. This article intends to give an overview of different techniques to solve these problems.

Intro

In the following, we use a simple example to demonstrate the usage of Cython, the Scipy module Weave, Ctypes and the Python C-API. The article should give an impression of the technique but not intends to replace howtos and the corresponding documentation.

Application Example

In this example, the computational expensive part of a photon generating Monte Carlo algorithm will be optimized. In this Monte Carlo experiment, we roll a dice in each timestep to discriminate between fluorescence (=0) and FRET (=1) as well as none of both. The probability for fluorescence is time-independent while the FRET probability has a time dependency. The pure python implementation looks as follows (python.py):

import random

def getPhoton(p0,pvar,start,seed):
  """
  p0:    probability for donor fluorescence
  pvar:  probability trajectory  for FRET + probability for donor fluorescence
  start: start position on trajectory
  seed:  random number generator seed
  """
  random.seed(seed)
  for curndx in range(start,len(pvar)):
    rnd=random.random()
    if rnd < p0:
      return 0
    if rnd < pvar[curndx]:
      return 1

  return -1

The function parameters are the fluorescence probability p0, the FRET probility pvar (list, tupel or numpy array), the starting time and the random number generator seed for the initialization of the RNG.

Cython

We now wanna improve the performance of the code. Cython converts python code in C, which is then compiled an used as Python module. First of all the unmodified python code (cython.pyx):

import random

def getPhoton(p0,pvar,start,seed):
  random.seed(seed)
  for curndx in range(start,len(pvar)):
    rnd=random.random()
    if rnd < p0:
      return 0
    if rnd < pvar[curndx]:
      return 1

    return -1

Translated to C-code by

cython cython.pyx

which is then compiled into a Python module via

gcc -shared -O3 -fPIC -I/usr/include/python2.7 -o cython.so cython.c

Static Typedeclaration

Unmodified Cython code is often not much faster as the original. One reason is the dynamic typing of Python. Cython allows static typing by extending the Python syntax. Therefore, by simply adding the static typing information, a significant speed improvement can be achieved. (ctypes-statictypedec.pyx):

cimport numpy as np
from numpy import random

def getPhoton(double p0,np.ndarray pvar,long start,long seed):
  cdef double rnd
  cdef long end,curndx

  end = len(pvar)  
  random.seed(seed)

  for curndx in range(start,end):
    rnd=random.random()
    if rnd < p0:
      return 0
    if rnd < pvar[curndx]:
      return 1

  return -1

As in C, standard datatypes such as int, long, can be declared. For special data types such as Numpy arrays, a special cython module has to be imported via the cimport line.

Including external libraries and/or code

Cython is pretty flexible when it comes to the inclusion of external code. In the following, we wanna include a SIMD unit optimized Mersenne Twister RNG. To achieve this, we copy the required signatures from the SFMT.h and change them as follows (cython-sfmt.pyx):

cdef extern from "SFMT.h":
  void init_gen_rand(int seed)
  double genrand_res53() 

cimport numpy as np

def getPhoton(double p0,np.ndarray pvar,long start,long seed):
  cdef double rnd
  cdef long end,curndx

  end = len(pvar)  
  init_gen_rand(seed)

  for curndx in range(start,end):
    rnd=genrand_res53()
    if rnd < p0:
      return 0
    if rnd < pvar[curndx]:
      return 1

  return -1

The cdef in the first 3 lines allows usage of the functions in the python code. We can now link the objectfiles of our compiled cython and SFMT code to a Python-C module:

gcc -shared -fPIC -o cython-sfmt.so SFMT.o cython-sfmt.o

In summary, Cython is particular suitable to optimize existing python code and to interface or embed external libraries and C-code.

Weave

OK, we now assume that we have the code already as a C-code snippet, which we wanna use from Python. A simple way to do this is using the weave module and its blitz type converter (weave.py):

from scipy.weave import inline
from scipy.weave import converters

def getPhoton(pconst,u,startpos,rngseed):    
  code = """   
         double rnd;
         int retval,i;
         retval=-1;
         srand(rngseed);

         for (i=startpos;i<Nu[0];i++)
         {
           rnd = rand()/((float)RAND_MAX);
           if (rnd < pconst) {
             retval=0;
             break;
           }
           if (rnd < u(i)){
             retval=1;
             break;
           }
         }
         return_val = retval;
         """

  return inline(code, ['pconst','u','startpos','rngseed'],
                type_converters=converters.blitz,compiler='gcc',
                headers=['<stdlib.h>',],)

As we can see, the getPhoton function exclusively consists of a multiline string assignment containing the C-code to the code Variable which is then compiled “on-the-fly” by the inline function. Parameters of the inline function are the codestring and the list of function parameters. Further, typeconverter, compiler and headers are specified here. Compilation is done during run-time and is cached, such that the code has to be compiled only once.

In my opinion, weave is particular useful for short code snippets. I found it a bit tricky to include external libraries and further, syntax highlighting of the c-part is typically not possible, as it is considered as Python string by the code parsers.

Ctypes

Lets now assume, that the photon generator code is in an external library. The ctypes module from the Python standard library is suitable for this task. Ctypes allows direct calls to shared library functions and a proper conversion of Python objects into their C equivalents. We now look into a very simple example. Ctypes also allows mapping of structures, specification of function in and output variables and much more.

from ctypes import cdll,c_double,byref,c_long,POINTER
import ctypes

myclib = cdll.LoadLibrary("libfretext.so")

def getPhoton(p0,pvar,start,seed):
  nvar = len(pvar)
  doublepointer = ctypes.POINTER(ctypes.c_double)

  c_pconst   = c_double(p0)
  c_pvarlen  = c_long(nvar)
  c_pvarp    = pvar.ctypes.data_as(doublepointer)
  c_startpos = c_long(start)
  c_rngseed  = c_long(seed)

  retval = myclib.tryGetPhoton(c_pconst,c_pvarlen, c_pvarp,c_startpos, c_rngseed)
  return retval

First of all, we “import” the library “libfretext.so” via the LoadLibrary function and render them accessible as myclib. The sole purpose of the following function getPhoton is to convert Python objects into their corresponding C types. The converted objects are then used for the function call.

Ctypes is very useful to call functions in available shared libraries. It allows encapsulation of available functions without modification of the existing library.

Python / Numpy C-API

The most detailed interfacing capabilities are available by using the Python C-API. The API can be used to write modules in C or even to embed the Python interpreter in a C-application. Our example looks like:

#include <Python.h>
#include <numpy/arrayobject.h>
#include "cext.h"
#include "SFMT.h"

static PyObject *getPhoton(PyObject *self, PyObject *args)
{
  PyArrayObject *pvar;
  long startpos,rngseed;
  double pconst,*pvarp,rnd;
  npy_intp *pvarlen;
  int retval,i,j;

  if (!PyArg_ParseTuple(args,"dOll",&pconst,&pvar,&startpos,&rngseed)) return NULL;

  /*Checking if array is C-Contiguous, double and then initializing pointers*/

  if (!PyArray_ISCONTIGUOUS(pvar) || (sizeof(double) != PyArray_ITEMSIZE(pvar))) return NULL;
  pvarp=(double*)PyArray_DATA(pvar);
  pvarlen = PyArray_DIMS(pvar);

  retval=-1;

  init_gen_rand(rngseed);
  for (i=startpos;i<pvarlen[0];i++){
    rnd = genrand_res53();

    if (rnd < pconst) {
      retval=0;
      break;
    }
    if (rnd < pvarp[i]){
      retval=1;
      break;
    }
  }

  return Py_BuildValue("i",retval);
}

static PyMethodDef fretmethods[] = {
  {"getPhoton", getPhoton,METH_VARARGS,"Generation attempt of a photon in plain C-Code"},
  {NULL,NULL,0,NULL}
};

#ifndef PyMODINIT_FUNC    /* declarations for DLL import/export */
#define PyMODINIT_FUNC void
#endif

PyMODINIT_FUNC initcext(void)
{
  (void) Py_InitModule("cext",fretmethods);
  import_array();
  printf("-> FRET numpy extension initialized using SF Mersenne Twister\n");
}

The shown code consists of 3 blocks: The function getPhoton, which processes and returns python objects, the PyMethodDef block to list the functions of the module and the initcext function, which is called upon module initialization. In the above example, parts of the Numpy C-API is used additionally.

Compilation and linkeage:

gcc -I/usr/lib64/python2.7/site-packages/numpy/core/include/ -I/usr/include/python2.7 -c -O3 -fPIC -D_FORTIFY_SOURCE=0 -DMEXP=607 -msse2 -DHAVE_SSE2 SFMT.c
gcc -I/usr/lib64/python2.7/site-packages/numpy/core/include/ -I/usr/include/python2.7 -c -O3 -fPIC -D_FORTIFY_SOURCE=0 -DMEXP=607 -msse2 -DHAVE_SSE2 cext.c
gcc -shared -fPIC -Wl,-soname,cext.so -o cext.so SFMT.o cext.o

If you want the maximum of control, flexibility and speed, you should familiarize yourself with the C-API. The API however requires a solid C background and knowledge of Python interna and their implementation in C. As an example, one has to take care for incrementing, and decrementing python object reference counters in the C-code.

  1. http://www.scipy.org/PerformancePython
  2. http://technicaldiscovery.blogspot.de/2011/06/speeding-up-python-numpy-cython-and.html
  3. http://www.scipy.org/Weave
  4. http://www.cython.org
  5. http://docs.python.org/library/ctypes.html
  6. http://docs.python.org/c-api/
  7. http://docs.scipy.org/doc/numpy/reference/c-api.html
  8. http://docs.python.org/extending/
  9. http://www.swig.org/
  10. http://riverbankcomputing.co.uk/static/Docs/sip4/index.html
  11. http://www.boost.org/doc/libs/1_49_0/libs/python/doc
  12. http://thejosephturner.com/blog/2011/06/02/embedding-python-in-c-applications-with-boostpython-introduction/