Benchmark ODE solver: GSL V/s Boost Odeint library

For our neural simulator, MOOSE, we use GNU Scientific Library (GSL) for random number generation, for solving system of non-linear equations, and for solving ODE system.

Recently I checked the performance of GSL ode solver V/s Boost ode solver; both using Runge-Kutta 4. The boost-odeint outperformed GSL by approximately by a factor of 4. The numerical results were same. Both implementation were compiled with -O3 switch.

Below are the numerical results. In second subplot, a molecule is changing is concentration ( calcium ) and on the top, another molecule is produced (CaMKII). This network has more than 30 reactions and approximately 20 molecules.

GSL took approximately 39 seconds to simulate the system for 1 year. While Boost-odeint took only 8.6 seconds. (This turns out to be a corner case)

Update: If I let both solvers choose the step size by themselves for given accuracy, boost usually outperforms the GSL solver by a factor of 1.2x to 2.9x. These tests were done during development of a signaling network which has many many reactions. Under no circumstances I have tested, BOOSE ODE solver was slower than GSL solver.


Multi-dimentional root fiding using Newton-Raphson method

In short, It does the same thing as GNU-gsl `gnewton` solver for finding roots of a system of non-linear system.

The standard caveats apply when finding the multi-dimensional roots.

The templated class (must be compiled with  -std=c++11). It uses boost::numeric::ublas library.

The code is available here:

The header file contains the class, file `test_root_finding.cpp` shows how to use it.

Getting caught in a loop while traversing a (directed acyclic) graph

How to detect if you are getting caught inside a loop when you are traversing a directed graph. This is equivalent of checking if you graph is a Directed Acyclic Graph (DAG).

Assume that each vertex of graph has a label such as A, B, C, etc. While you are traversing the graph, you keep your labels of traversed vertices on a stack or in a vector. If you construct a string out of this traversal, you will get a pattern like the following one.


Whenever you are caught in a loop, a path will repeat itself. We have D,E,F repeating itself. Now there are two ways to detect this. One is a classical graph based method which is efficient O(E+V). The second one, about which I am going to talk about now, is not as efficient as first but but it is very easy to code. It is based on regular expressions.

All we have to do is  construct a string out of the data-structure which holds the vertex traversal order and match it against a regular-pattern which looks for the repeated substrings.

Here is one such pattern written using boost regex library.

 boost::regex rx ("^([0-9a-zA-Z_@]+,)+(?[0-9a-zA-Z@_,]+)(g{cycle})$", boost::regex::perl);

Care has to be taken while using this regex. First, I have special character @ in the label of vertex. And this pattern assumes that if there is a single repetition such as A,B,A,B or A,B,B etc, then it must lead to a loop. In my case it is justified. You should think if this is the behavior you want. And again, it is not as efficient as standard graph algorithms to detect loop (which uses Depth First Search). But for writing quick tests, it is a good approach.

I have used boost regex library. I am already using boost graph library therefore I need not look for other libraries. I have heard a some good things about GNU-regex library but have never used myself. Also, my gcc (4.6) does not support 2011C++ fully e.g. regex is not implemented fully. Have a look at its status page before writing any regex using header <regex>.