/* I assume that you have a class called RngClass that
* can produce standard normal deviates (accessed with
* RngClass::Normal(0,1)). Add this member to your class
* RngClass and off you go. Note that when alpha != 0.0,
* m is NOT the mean and s is also NOT the standard deviation.
*/#include <cmath> // for sqrt, probably already done in the headerdouble RngClass::SkewNormal(double m, double s, double alpha){/* If s < 0.0, give an error. I throw an exception, but
* you may want to do something else. Note that I have defined
* the object rngException myself. Take a look online at the struct
* std::exception for more info.
*/if( s <0.0)throw rngException();/* if s == 0, then just return m
*/if( s ==0.0)return m;/* now handle non-trivial cases
*/double Z1 = Normal(0.0,1.0);/* if alpha == 0.0, then just return Normal(m,s)
*/if( alpha ==0.0)return s * Z1 + m;else{/* X = 1/sqrt{1+alpha^2} alpha |Z1| + Z2
*/double Z2 = Normal(0.0,1.0);double absZ1 =(Z1 >0.0? Z1 :-Z1);return m + s *((alpha*absZ1 + Z2)/sqrt(1.0+ alpha*alpha));}}
Make and makefiles
... and for those of you that get sweaty palms when you have to run gcc to compile your code or if you completely panic when you open the makefile that someone else kindly stitched together for your project and you see :
%.o :%.cpp
g++-g -o $@-c $<
... here is a nice introduction that guides you through the fundamentals of preprosessing, compiling, assembling and linking, setting your library and include paths right and how to write those darn makefiles in an understandable way: link
An example of a makefile that 1) puts the .o files in a seperate folder, 2) handles dependencies automatically and 3) compiles all .cpp files in a source folder:
# An "omnipotent" makefile for g++ that builds the# dependecy tree "automagically".# I assume that you have a subfolder ./src that contains# all your .cpp and .hpp files, and that you made a subfolder# ./obj that will contain the .o files.# Hint: the compiler option -MMD does a lot of magic.# -----------------------------------------------------------# the compiler, in my case the C++ compiler
CC = g++# flags for compiling .cpp files
CC_FLAGS =-O3 -Wall -c
# flags for linking the .o files
LD_FLAGS =-O3 -Wall
# all my .cpp files are saved in a folder ./src
CPP_FILES :=$(wildcard src/*.cpp)# and I want my .o files to be stored in a ./obj folder
OBJ_FILES :=$(addprefix obj/,$(notdir$(CPP_FILES:.cpp=.o)))# the target "all" is my main executable called# "hawaiitoaster" (make, make all)
all: hawaiitoaster
# a rule for cleaning up (make clean)
clean:
rm -f obj/*.o obj/*.d hawaiitoaster
# the rule for the main executable
hawaiitoaster:$(OBJ_FILES)$(CC)$(LD_FLAGS)$^-o $@# a pattern for compiling .cpp files# (and hence, making the .o files)
obj/%.o: src/%.cpp $^$(CC)$(CC_FLAGS)$<-o $@# now include dependencies. "gcc -c -MMD src/foo.cpp -o obj/foo.o"# will make a file foo.d (and foo.o). The foo.d file contains the# line "obj/foo.o: src/foo.cpp src/foo.hpp src/bar.hpp", depending# on whatever has been "#include"d in the header "foo.hpp".
CC_FLAGS +=-MMD
-include$(OBJ_FILES:.o=.d)
NB: When you copy this makefile, the indentation might consist of spaces. make wants tabs! Or download the makefile here
Multi-threading
Our machines have 8, 12, 16 or even 20 cores, nowadays! Together we have way more cores than (PhD) students. This means that you can run a gazillion simulations simultaneously, or add some multi-threading to your c++ program.
Concurrency can be extremely complicated, and causes problems that will haunt you in your dreams. [TODO: better link about the basic problems] 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 [TODO: link to python wiki page]. 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 explain 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();// destructorint 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...};
Fun inherits the worker pool functionality from "WorkerPoolWrapper" by writing " : public WorkerPoolWrapper" just after "class Fun". The class "Fun has a "pure virtual" member executeJob(void* ). This means that no object of class WorkerPoolWrapper can be created. 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. 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// ...
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 continuedint Fun::operator()(int n){int* ns =newint[n];// make array with integers 1,...,nfor(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 completedint pi =0;// gather (transformed) data and compute final resultfor(int i =0; i < n;++i ) pi += ns[i];delete[] ns;// free allocated spacereturn pi;// return answer}// ...
Notice that this implementation of $\pi(n)$ is not the most efficient one. It checks for every integer $i$ between $0$ and $n$ if $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, 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 continuedvoid 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}
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
Programming in C++
Custom Random Number Generators
Have you always wanted to sample from a Skew Normal Distribution? I have never been able to find one online, so I made one my self: See also
http://en.wikipedia.org/wiki/Skew_normal_distribution
Make and makefiles
... and for those of you that get sweaty palms when you have to run gcc to compile your code or if you completely panic when you open the makefile that someone else kindly stitched together for your project and you see :
... here is a nice introduction that guides you through the fundamentals of preprosessing, compiling, assembling and linking, setting your library and include paths right and how to write those darn makefiles in an understandable way: link
An example of a makefile that 1) puts the .o files in a seperate folder, 2) handles dependencies automatically and 3) compiles all .cpp files in a source folder:
NB: When you copy this makefile, the indentation might consist of spaces. make wants tabs! Or download the makefile here
Multi-threading
Our machines have 8, 12, 16 or even 20 cores, nowadays! Together we have way more cores than (PhD) students. This means that you can run a gazillion simulations simultaneously, or add some multi-threading to your c++ program.
Concurrency can be extremely complicated, and causes problems that will haunt you in your dreams. [TODO: better link about the basic problems] 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 [TODO: link to python wiki page]. 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 explain 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 "WorkerPoolWrapper" by writing " : public WorkerPoolWrapper" just after "class Fun". The class "Fun has a "pure virtual" member executeJob(void* ). This means that no object of class WorkerPoolWrapper can be created. 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. 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$ if $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, 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 source for WorkerPoolWrapper can be downloaded here: WorkerPoolWrapper.hpp and here: WorkerPoolWrapper.cpp