ETA for C++

Simulations can take some time, and I'd like to know how long. This is easy, right? Yes, it is. I've done it lots of times, but every time I do, I curse myself for not using an old piece of code.
most likely, there is some standard, best way of doing this, but I haven't found it. Most recently, I did this: I made a simple object "EtaEstimator", that can be updated every (costly) time step and asked for an estimated time of "arrival" at any time. Here's the header:

// eta.hpp
#include <ctime>
#include <cmath> // floor
#include <iostream>

class EtaEstimator {
public:
    EtaEstimator(int ); 
    // constuction starts the clock. Pass the number of steps
    void update();
    void print(std::ostream & ) const;
private:
    double ct, etl; // cumulative time, estimated time left
    int n, N; // steps taken, total amount of steps
    clock_t tick; // time after update ((c) matlab)
    // statics...
    static const int secperday = 86400;
    static const int secperhour = 3600;
    static const int secperminute = 60;
};

std::ostream & operator<<(std::ostream & , const EtaEstimator & );

The members are straight forward too:

// eta.cpp
#include "eta.hpp"

EtaEstimator::EtaEstimator(int N) :
        ct(0.0), etl(0.0), n(0), N(N) {
    tick = clock();
}

void EtaEstimator::update() {
    clock_t dt = clock() - tick;
    tick += dt;
    ct += (double(dt)/CLOCKS_PER_SEC); // prevent integer division
    // CLOCKS_PER_SEC is defined in ctime
    ++n;
    etl = (ct/n) * (N-n);
}

void EtaEstimator::print(std::ostream & os) const {
    double etlprime = etl;
    int days = floor(etlprime / secperday);
    etlprime -= days * secperday;
    int hours = floor(etlprime / secperhour); 
    etlprime -= hours * secperhour;
    int minutes = floor(etlprime / secperminute);
    etlprime -= minutes * secperminute;
    int seconds = floor(etlprime);
    os << (days > 0 ? std::to_string(days) + " " : "")
       << hours << ":" 
       << (minutes < 10 ? "0" : "") << minutes << ":" 
       << (seconds < 10 ? "0" : "") << seconds;
}

std::ostream & operator<<(std::ostream & os,
        const EtaEstimator & eta) {
    eta.print(os);
    return os;
}

Typical usage of EtaEstimator would be the following:

#include <iostream>
#include "eta.hpp"

// about to do lots of work...
int N = 1000;
EtaEstimator eta(N);
for ( int n = 0; n < N; ++n ) {
    // do something very expensive
    eta.update()
    std::cout << "\rETA: " << eta << std::flush;
}
// ...

PS: std::to_string is a C++11 feature, and can be ignored by using something like

if ( days > 0 ) os << days << " "; // else nothing at all

Carnes's life span distribution

In a paper by Carnes et. al, a simple parametric, but realistic life span distribution is given, and here I show how you can sample from it. In addition, assuming a demographic equilibrium, the age of individuals will have a particular distribution. I will show what this distribution is, and again how to sample from it. Sampling ages instead of lifespans might be useful for initiating simulations. I model epidemics, and I want my disease free (a.k.a. virgin) population to have the 'right' age distribution.
The life span distribution has hazard $$\lambda(a) = e^{u_1 a + v_1} + e^{u_2 a + v_2}\,.$$ Typical parameters are given by $u_1 = 0.1$, $v_1 = -10.5$, $u_2 = -0.4$, and $v_2 = -8$, so that infants have a slightly increased hazard of dying, and after the age of 60, the hazard rapidly starts to grow, until it becomes exceedingly large around $a = 100$.
The survival function $S(a) = {\rm Pr}(A > a)$, where $A$ is the random variable denoting 'age at death' is given by $S(a) = e^{-\Lambda(a)}$, with $\Lambda(a) := \int_0^a \lambda(a')da'$ denoting the cumulative hazard. The cumulative hazard $\Lambda$ is easily calculated: $$\Lambda(a) = \frac{e^{v_1}}{u_1}(e^{u_1 a}-1) + \frac{e^{v_2}}{u_2}(e^{u_2 a}-1)\,,$$ but its inverse, or the inverse of the survival function is more difficult to calculate.
We need the inverse of $S$, because sampling random deviates typically involves uniformly sampling a number $u\in [0,1)$. The number $S^{-1}(u)$ is then the desired deviate.

In a future post, I will show how to use the GNU Scientific Library (GSL) to sample deviates from $A$.

Suppose that the birth rate $b$ in our population is constant. A PDE describing the population is given by $$\frac{\partial N(a,t)}{\partial t} + \frac{\partial N(a,t)}{\partial a} = -\lambda(a)N(a,t)\,,$$ where $N(a,t)$ is the number (density) of individuals of age $a$, alive at time $t$. The boundary condition (describing birth) is given by $$N(0,t) = b$$ When we assume that the population is in a demographic equilibrium, the time derivative with respect to $t$ vanishes, and we get an ODE for the age distribution: $$\frac{\partial N(a)}{\partial a} = -\lambda(a) N(a)\,,\quad N(0) = b\,,$$ where we omitted the variable $t$. This equation can be solved: $$\frac{1}{N}\frac{\partial N}{\partial a} = \frac{\partial \log(N)}{\partial a} = -\lambda(a) \implies N(a) = c \cdot e^{-\Lambda(a)}$$ for some constant $c$. Since $b = N(0) = c \cdot e^{-\Lambda(0)} = c$, we have $$N(a) = b\cdot e^{-\Lambda(a)}\,.$$ Hence, we now know the PDF of the age distribution (up to a constant). Unfortunately, we can't get a closed form formula for the CDF, let alone invert it. Therefore, when we want to sample, we need another trick. I've used a method from Numerical Recipies in C. It involves finding a dominating function of the PDF, with an easy, and easily invertible, primitive.
Let's just assume that $b = 1$, so that the objective PDF is $N(a) = e^{-\Lambda(a)}$. Please notice that $N$ is not a proper PDF, since, in general, the expectation won't be $1$. We need to find a simple, dominating function for $N$. A stepwise defined function might be a good choice, since the hazard is practically zero when the age is below $50$, and then increases rapidly. We first find a comparison cumulative hazard $\tilde{\Lambda}$ that is dominated by the actual cumulative hazard $\Lambda$. Many choices are possible, but one can take for instance $$\tilde{\Lambda}(a) = \left\{ \begin{array}{ll} 0 & \mbox{if } a < a_0 \\ \lambda(a^{\ast})\cdot (a-a_0) & \mbox{otherwise} \end{array}\right.$$ where $$a_0 = a^{\ast} - \frac{\Lambda(a^{\ast})}{\lambda(a^{\ast})}\,.$$ The constant $a_0$ is chosen such that the cumulative hazards $\Lambda$ and $\tilde{\Lambda}$ are tangent at $a^{\ast}$.
Since $\Lambda$ dominates $\tilde{\Lambda}$, the survival function $\tilde{S}$ defined by $\tilde{S}(a) = e^{-\tilde{\Lambda}(a)}$ dominates $S$. It is easy to find the inverse of $a\mapsto\int_0^a \tilde{S}(a')da'$, and hence we can easily sample random deviates from the age distribution corresponding to $\tilde{\Lambda}$. In order to sample from the desired age distribution, one can use a rejection method: (i) sample an age $a$ from the easy age distribution. (ii) compute the ratio $r = S(a)/\tilde{S}(a) \leq 1$. (iii) sample a deviate $u \sim \mbox{Uniform}(0,1)$. (iv) accept the age $a$ when $u \leq r$, and reject $a$ otherwise. Repeat these steps until an age $a$ was accepted.

The only thing we still need to do, is to find a good value for $a^{\ast}$. To make the sampler as efficient as possible, we want to minimize the probability that we have to reject the initially sampled age $a$ from $S$. This boils down to minimizing $$\int_0^{\infty} \tilde{S}(a)da = a_0 + \frac{1}{\lambda(a^{\ast})} = a^{\ast} + \frac{1 - \Lambda(a^{\ast})}{\lambda(a^{\ast})}\,.$$ The derivative of $a^{\ast} \mapsto \int_0^{\infty} \tilde{S}(a)da$ equals $$\frac{(\Lambda(a^{\ast}) - 1)\lambda'(a^{\ast})}{\lambda(a^{\ast})^2}$$ and thus, we find an extreme value for $\int_0^{\infty} \tilde{S}(a)da$ when $\Lambda(a^{\ast}) = 1$ or when $\lambda'(a^{\ast}) = \frac{d\lambda}{da^{\ast}} = 0$. The second condition can only correspond to a very small $a^{\ast}$, and therefore will not minimize $\int_0^{\infty} \tilde{S}(a)da$. Hence, we have to solve $\Lambda(a^{\ast}) = 1$. When we ignore the second term of $\Lambda$, we find that: $$a^{\ast} = \log(1 + u_1 \exp(-v_1))/u_1$$

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.

	// main.hpp
	#include <list>
	#include <iostream>
	#include <cmath> // for sqrt and pow
	// the population consists of "Critters"
	class Critter {
	public:
	Critter(int ); // constructor takes "age"
	// typical critter traits...
	double happiness() const;
	double wealth() const;
	private:
	int age;
	};
	/* Trait will be a const method of Critter
	* that takes no argument, and returns a double.
	*/
	typedef double (Critter::*Trait)() const;
	// higher-order functions take a Trait as an argument
	double calcMean(const std::list<Critter> & , Trait );
	double calcCorrelation(const std::list<Critter> & , Trait , Trait );
	int main();

view raw main_critters.hpp hosted with ❤ by GitHub

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)()"

	// main.cpp
	#include "main.hpp"
	// Critters constructor
	Critter::Critter(int age) : age(age) { /* empty */ }
	// methods for Critter (Traits)
	double Critter::happiness() const {
	return 10.0/0.9*pow(0.9, age/10.0); // some arbitrary computation
	}
	double Critter::wealth() const {
	return (age-50)*(age-50)*0.01 + age*0.1; // some arbitrary computation
	}
	// higher-order functions
	double calcMean(const std::list<Critter> & population, Trait x) {
	double Ex = 0.0;
	int n = population.size();
	std::list<Critter>::const_iterator it = population.begin();
	for ( ; it != population.end(); ++it ) {
	Ex += ((*it).*x)(); // x needs to be derefecenced before it can be called
	}
	Ex /= n; // assumes a non-empty population
	return Ex;
	}
	double calcCorrelation(const std::list<Critter> & population, Trait x, Trait y) {
	double Ex = 0.0; double Ey = 0.0;
	double Exx = 0.0; double Eyy = 0.0;
	double Exy = 0.0;
	int n = population.size();
	std::list<Critter>::const_iterator it = population.begin();
	for ( ; it != population.end(); ++it ) {
	Ex += ((*it).*x)(); Exx += ((*it).*x)() * ((*it).*x)();
	Ey += ((*it).*y)(); Eyy += ((*it).*y)() * ((*it).*y)();
	Exy += ((*it).*x)() * ((*it).*y)();
	}
	Ex /= n; Ey /= n; Exx /= n; Eyy /= n; Exy /= n;
	double sdx = sqrt(Exx - Ex*Ex); double sdy = sqrt(Eyy - Ey*Ey);
	double covxy = Exy - Ex*Ey;
	double corrxy = covxy/(sdx*sdy);
	return corrxy;
	}
	// main
	int main() {
	// make a population of Critters
	int N = 100;
	std::list<Critter> critters;
	for ( int n = 0; n < N; ++n ) {
	critters.push_back(Critter(n));
	}
	// calculate statistics
	double Eh = calcMean(critters, &Critter::happiness);
	double Ew = calcMean(critters, &Critter::wealth);
	double corrhw = calcCorrelation(critters, &Critter::happiness, &Critter::wealth);
	// and show them to the world
	std::cout << "mean wealth = " << Ew << std::endl;
	std::cout << "mean happiness = " << Eh << std::endl;
	std::cout << "correlation happiness and wealth = " << corrhw << std::endl;
	return 0;
	}

view raw main_critters.cpp hosted with ❤ by GitHub

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).

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.

	import math
	import matplotlib as mpl
	import matplotlib.pyplot as plt
	import numpy as np
	def blob_plot(ax,data,list_pos, extra=False):
	lextra=[]
	for d,p in zip(data,list_pos):
	#this makes the histogram that constitutes the blob plot
	m=min(d)
	M=max(d)
	nbins=20
	x = np.linspace(m,M,nbins) # support for histogram
	his,bins = np.histogram(d, bins=nbins)
	#scale the histogram to a proper size (you may have to fiddle with this), then place it at position 'pos'
	scale_blob=2.
	shift_his_plus_pos = [ p + h*scale_blob/float(len(d)) for h in his]
	shift_his_minus_pos = [ p - h*scale_blob/float(len(d)) for h in his]
	facecolor,alpha=color_alpha_blobs # assign color and transparency
	#this is the matplotlib function that does the trick
	ax.fill_betweenx(x,shift_his_minus_pos, shift_his_plus_pos, linewidth=0.0, facecolor= facecolor, alpha=alpha)
	#calculates median or mean, if you want
	if extra=='median': lextra.append( np.median(d) )
	elif extra=='mean': lextra.append( np.mean(d) )
	#and plots it
	if extra != False:
	color = 'orangered'
	ax.plot(list_pos,lextra, color=color_mean_or_median,linestyle='-',marker='D',markersize=5, lw=1.5)
	# MAKE UP SOME DATA
	list_positions=[0,1,2,3,4] #this is their position on the x-axis
	data_for_blobplot=[] # this will contain a list for each blob we want in the plot
	extra='median'
	for pos in list_positions:
	if pos <= 2: some_data = np.random.gumbel(pos,1, 1000) # generate some data
	else: some_data = np.random.normal(pos,1, 1000) # generate some other data
	data_for_blobplot.append(some_data) # and append it to the list
	# PLOTTING
	color_alpha_blobs=('midnightblue',0.7) #this is the color of the blobplot
	color_mean_or_median='orangered' #this is the colot of the median or mean, whatever you specify
	fig = plt.figure() # makes a figure
	ax = fig.add_subplot(111) # adds a single subplot
	#this is the blobplotting function
	blob_plot(ax, data_for_blobplot, list_pos=list_positions, extra=extra) #makes blob plots, if extra is not spcified, makes no median/mean
	ax.set_xlim(-1., len(list_positions)) #to display all of them
	ax.set_ylim(-2., 12.)
	plt.title('Beautiful blob plots',fontsize=18)
	plt.xlabel('$\mu$ (yup, that is Latex)', fontsize=14)
	plt.ylabel('distribution of my variable', fontsize=14)
	####################################################
	# So far this is all you need to make the blob plot
	# Below is if you want a legend for the medain/mean
	####################################################
	# draw temporary dot to use it for a legend, later set to invisible
	h = ax.scatter( [1],[1], marker='D', color=color_mean_or_median, label='extra')
	#draws the legend, with larger symbols
	leg=plt.legend(loc= 'upper left',scatterpoints=1)
	for legobj in leg.legendHandles:
	legobj.set_linewidth(10.0)
	h.set_visible(False)
	plt.savefig('blobplot.png') # saving to png, alternatively you can write 'blobplot.pdf', or other formats
	plt.show() #prints to screen your splendid blobplot

view raw blobplot.py hosted with ❤ by GitHub

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)

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.hpp
	#include "WorkerPoolWrapper.hpp"
	class Fun : public WorkerPoolWrapper {
	public:
	Fun(int ); // constructor, pass the number of threads that you want to use
	~Fun(); // destructor
	int operator()(int ); // suppose that this class is used for a functionoid
	// other public members and methods...
	protected:
	void executeJob(void* ); // pure virtual function in the base class
	// other protected members and methods...
	};

view raw Fun.hpp hosted with ❤ by GitHub

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:

	// Fun.cpp
	#include "Fun.hpp"
	Fun::Fun(int w) { initWorkerPool(w); } // upon construction, Fun initiates the worker pool with w workers
	Fun::~Fun() { waitForWorkerPoolToExit(); } // upon destruction, Fun joins the threads
	// ...

view raw Fun_part1.cpp hosted with ❤ by GitHub

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:

	// ... Fun.cpp continued
	int Fun::operator()(int n) {
	int* ns = new int[n]; // make array with integers 1,...,n
	for ( int i = 0; i < n; ++i ) ns[i] = i+1;
	for ( int* job = ns; job < ns+n; ++job ) addNewJob((void*) job); // add the integers as jobs to the job que
	syncWorkerThreads(); // wait for all jobs to be completed
	int pi = 0; // gather (transformed) data and compute final result
	for ( int i = 0; i < n; ++i ) pi += ns[i];
	delete[] ns; // free allocated space
	return pi; // return answer
	}
	// ...

view raw Fun_part2.cpp hosted with ❤ by GitHub

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:

	// ... Fun.cpp continued
	void Fun::executeJob(void* job) {
	int* i = (int*) job; // cast the job back to int*
	bool prime = true;
	int d = 2;
	while ( prime && d != (*i) ) prime = (*i) % (d++); // naive prime test
	*i = prime; // store the result in the job
	}

view raw Fun_part3.cpp hosted with ❤ by GitHub

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:

	// main.cpp
	// ...
	Fun primePi(10); // uses 10 cores!
	int pi = primePi(144169); // equals 13355
	// ... when primePi gets out of scope, the threads are joined

view raw main_workerpoolwrapper.cpp hosted with ❤ by GitHub

The class WorkerPoolWrapper is declared here:

	#ifndef WORKERPOOLWRAPPER_HPP
	#define WORKERPOOLWRAPPER_HPP
	#include <pthread.h>
	#include <list>
	class WorkerPoolWrapper {
	public:
	/* constuctor */
	WorkerPoolWrapper();
	/* destructor */
	virtual ~WorkerPoolWrapper();
	/* initiate private members and set the number of workers.
	* Returns true if the threads were successfully started,
	* false if there was an error starting a thread
	*/
	bool initWorkerPool(int );
	/* wait for all workers to finish. First call sendNoMoreJobsSignal
	* to make sure that they do.
	*/
	void waitForWorkerPoolToExit();
	/* Wait for all Workers to finish their jobs. The Workers
	* will then be waiting for new jobs. (So no joining...)
	*/
	void syncWorkerThreads();
	/* Send the worker pools a signal that there will be no more jobs.
	* This will make them exit their WorkerThreadEntry function.
	*/
	void sendNoMoreJobsSignal();
	/* add a new job to the list of jobs.
	*/
	void addNewJob(void* );
	protected:
	/* the member executeJob needs to be implemented by the class that
	* inherits this class. It should do the actual work that needs
	* to be done in parallel.
	*/
	virtual void executeJob(void* )=0;
	private:
	/* pthread is a C library, and does not know about classes. Therefore
	* we need the following construction
	*/
	void workerThreadEntry();
	static void* workerThreadEntryFunc(void* );
	/* mutex and conditional variable for the job list */
	pthread_mutex_t jobsMutex;
	pthread_cond_t jobsConditionalVar;
	/* mutex and conditional variable for synchronizing */
	pthread_mutex_t jobsyncMutex;
	pthread_cond_t jobsyncConditionalVar;
	/* the number of jobs and the job list */
	volatile int jobsyncCounter;
	std::list<void*> jobs;
	/* a flag for terminating the threads */
	volatile bool noMoreJobsFlag;
	/* an array of worker threads */
	pthread_t* workerThreads;
	pthread_attr_t attr;
	int numOfWorkerThreads;
	};
	#endif

view raw WorkerPoolWrapper.hpp hosted with ❤ by GitHub

And the members are defined here:

	#include "WorkerPoolWrapper.hpp"
	WorkerPoolWrapper::WorkerPoolWrapper() {
	numOfWorkerThreads = 0;
	noMoreJobsFlag = true;
	jobsyncCounter = 0;
	workerThreads = NULL;
	}
	bool WorkerPoolWrapper::initWorkerPool(int numOfWorkerThreads) {
	this->numOfWorkerThreads = numOfWorkerThreads;
	workerThreads = new pthread_t[numOfWorkerThreads];
	noMoreJobsFlag = false;
	jobsyncCounter = 0;
	pthread_attr_init(&attr);
	pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
	pthread_mutex_init(&jobsMutex,NULL);
	pthread_mutex_init(&jobsyncMutex,NULL);
	pthread_cond_init(&jobsConditionalVar,NULL);
	pthread_cond_init(&jobsyncConditionalVar,NULL);
	if ( numOfWorkerThreads < 1 ) return false;
	else {
	bool ok = true;
	for ( int i = 0; i < numOfWorkerThreads; i++ ) {
	ok = ok && ( pthread_create(workerThreads+i, &attr, workerThreadEntryFunc, this) == 0 );
	}
	return ok;
	}
	}
	WorkerPoolWrapper::~WorkerPoolWrapper() {
	delete[] workerThreads; // assumes that the threads have been joined!
	pthread_attr_destroy(&attr);
	pthread_mutex_destroy(&jobsMutex);
	pthread_mutex_destroy(&jobsyncMutex);
	pthread_cond_destroy(&jobsConditionalVar);
	pthread_cond_destroy(&jobsyncConditionalVar);
	}
	void WorkerPoolWrapper::waitForWorkerPoolToExit() {
	sendNoMoreJobsSignal(); // the user could have given this signal manually
	for ( int i = 0; i < numOfWorkerThreads; i++ ) {
	pthread_join(workerThreads[i], NULL);
	}
	}
	void WorkerPoolWrapper::syncWorkerThreads() {
	pthread_mutex_lock(&jobsyncMutex);
	while ( jobsyncCounter > 0 ) {
	pthread_cond_wait(&jobsyncConditionalVar, &jobsyncMutex);
	}
	pthread_mutex_unlock(&jobsyncMutex);
	}
	void WorkerPoolWrapper::sendNoMoreJobsSignal() {
	pthread_mutex_lock(&jobsMutex);
	noMoreJobsFlag = true;
	pthread_cond_broadcast(&jobsConditionalVar); // wake all waiting threads
	pthread_mutex_unlock(&jobsMutex);
	}
	void WorkerPoolWrapper::addNewJob(void* job) {
	pthread_mutex_lock(&jobsyncMutex);
	jobsyncCounter++;
	pthread_mutex_unlock(&jobsyncMutex);
	pthread_mutex_lock(&jobsMutex);
	jobs.push_back(job);
	pthread_cond_signal(&jobsConditionalVar); // wake a waiting thread
	pthread_mutex_unlock(&jobsMutex);
	}
	void* WorkerPoolWrapper::workerThreadEntryFunc(void* thisObject) {
	((WorkerPoolWrapper*) thisObject)->workerThreadEntry();
	return NULL;
	}
	void WorkerPoolWrapper::workerThreadEntry() {
	void* myJob = NULL;
	bool lookForAJob = true;
	while ( lookForAJob ) {
	if ( myJob != NULL ) {
	executeJob(myJob); // execute the job
	myJob = NULL; // reset my job
	pthread_mutex_lock(&jobsyncMutex);
	jobsyncCounter--;
	pthread_cond_signal(&jobsyncConditionalVar);
	pthread_mutex_unlock(&jobsyncMutex);
	}
	else {
	pthread_mutex_lock(&jobsMutex);
	if ( jobs.empty() ) {
	if ( noMoreJobsFlag ) lookForAJob = false;
	else {
	pthread_cond_wait(&jobsConditionalVar, &jobsMutex);
	if ( !jobs.empty() ) {
	myJob = jobs.front();
	jobs.pop_front();
	}
	}
	}
	else {
	myJob = jobs.front();
	jobs.pop_front();
	}
	pthread_mutex_unlock(&jobsMutex);
	}
	} // while ( lookForAJob )
	}

view raw WorkerPoolWrapper.cpp hosted with ❤ by GitHub

The credits for combining pthread with C++ classes go to Jeremy Friesner.