Friday, 18 September 2015

Basic statistics using higher-order functions in C++

I do a lot of individual-based simulations, and often I want to calculate some population statistics 'on the fly'. I found that it can be helpful to use C/C++'s basic functional capabilities.

A higher-order function is a function that takes other functions as arguments (or returns a function as a result). In C/C++, functions can be passed to other functions, but the notation can be a bit awkward, as one needs to pass a reference to a function. If the function is a method of some class, then the notation gets even more involved. You can make your life easier by using a typedef.

The following code snippet shows the header file of a simple example. The goal is to calculate some statistics on a population of "Critters". These Critters have multiple "traits", and the traits are accessed by methods of the class Critter of signature "double Critter::trait() const". Suppose that we want to calculate the mean of the trait "happiness". It's trivial to write a function that does this, but then we might also get interested in the average "wealth". The function that calculates the average wealth is identical to the function that calculates average happiness, except that happiness is replaced by wealth. We can get rid of this redundancy by defining the typedef Trait as a method of Critter, and writing a function that takes the average of an arbitrary Trait.

Let us now look at the source file. The most important things to notice are...
(1) whenever you pass a member "trait" (e.g. wealth) of Critter to a function, you should pass it as "&Critter::trait" (i.e. pass the reference).
(2) when you want to evaluate Trait x of Critter alice, you'll need to de-reference x, and call the resulting function: "(alice.*x)()"

If you want to play with this example, put the header in a file called "main.hpp", and the source in a file called "main.cpp", and compile main.cpp by typing "g++ main.cpp -o critters" in your terminal (I assume that you are using Linux and have the gcc compiler installed).

Friday, 10 July 2015

Blob plots, like violin plots, but more precise

Violin plots are beautiful, useful and much clearer than drawing several histograms on the same plot.
In my opinion, the smoothing function that they implement with python does not really compensate the loss of precision with the visual appealing.
Why not plotting the whole distribution, then?

I modified the violin plots source code that I found here HERE

And produced... let's called them blob plots:


Here is the code.

Feel free to take it and do whatever you want with it.
If you have ideas to improve it let me know!
If you like the result, go say thanks to the violin plot maker (I did very little)

Friday, 26 June 2015

A simple class wrapper for parallelism in C++

Concurrency can be extremely complicated, and causes problems that will haunt you in your dreams. The classical libraries in C/C++ don't protect you from this horror in any way, and you will have to figure it out yourself. Parallelism is supposed to be a lot easier, but C/C++ does not have standard libraries like---for instance---Pythons parallelism modules. The boost libraries, however, provide an extensive interface for concurrency and parallelism.
If you don't want to use boost, don't panic. There are other options. The POSIX pthread library provides a down-to-the-metal concurrency model. It took me a while to find out what it is all about, and I haven't successfully applied all it's possibilities. What I have managed to apply, is a so-called "worker-pool" model. This is one of the easier concurrency applications (it falls under the parallelism category) of the pthread library, but can be quite useful. Here I will demonstrate a "wrapper" that C++ classes can inherit.

Suppose that you have a class Fun, that needs to do some computations. We declare Fun in---say---Fun.hpp:
Fun inherits the worker pool functionality from the class "WorkerPoolWrapper" (imported from WorkerPoolWrapper.hpp). The class WorkerPoolWrapper has a "pure virtual" member executeJob(void* ). You, the user, must specify what you want to compute in your own definition of the executeJob method. Besides implementing executeJob, you must also initiate the worker pool, and somewhere before Fun's destuctor returns, the threads must be "joined", i.e., all worker threads must finish their execution. In this case, I use the constructor and destructor of Fun to accomplish these things:
The methods initWorkerPool and waitForWorkerPoolToExit are inherited from WorkerPoolWrapper. Lets use Fun to compute the number of primes pi(n) below a number n. We overload operator() as follows:
Notice that this implementation of pi(n) is not the most efficient one. It checks for every integer i between 0 and n whether i is prime or not. This prime test is performed by executeJob. In the background, WorkerPoolWrapper has a list of jobs that have to be executed. Jobs can be added to this job list using addNewJob(void* ). Once executed, the result of a job must somehow be stored in the job again. Above, the number pi is the sum of the array ns, which makes sense when we look at the implementation of executeJob:
Hence, executeJob transforms the number i pointed to by job into 0 if i is composit, or 1 if i is prime, such that the sum of the i's equals pi(n). Before we gathered the results in Fun::operator(), we called syncWorkerThreads(). This method lets the program halt until every job in the job list has been executed.
Using the functionoid Fun now works as follows:
The class WorkerPoolWrapper is declared here:
And the members are defined here: The credits for combining pthread with C++ classes go to Jeremy Friesner.