Friday, April 17, 2009

Parallel version of the Gillespie solver

A few months back I blogged about my implementation of the Gillespie SSA algorithm for solving stochastic differential equations. I have now improved the code allowing for parallel execution on SMP systems. If you don't know what this is all about, please refer the the first post for a quick intro to what a stochastic differential equation is and what it means to numerically solve a system of them.

The code of the solver is in this module. There is an example in the package which is also worth running. On a two core machine the simulation runs about twice as fast, a expected. The parallelization was implemented using the multiprocessing module, and ammounts to roughly 4 lines of code. I love multiprocessing. and it's also compatible with psyco which makes the simulations run twice as fast (on top of the parallelization gains). Probably my next shot at improving its performance will be using scipy.weave and migrate the core of the algorithm to C. Simpler optimizations schemes such as cython (my quick port), runs twice as slow as my python code. Cython can be quite tricky to optimize beyond Python performance, and generates horribly verbose C code (if you don't believe me, try running cython on a print "hello world", you will get a 1000+ lines hello.c !!!).

It is worth noting that the code for the solver and examples are being maintained as a subpackage of BIP, my Bayesian inference package. Even though it is a bit unrelated, Since I frequently have to do Bayesian inference on mathematical models (both deterministic and stochastic), I thought it would be a good idea to keep the solver to those models in the same package as the analysis tools.