Back to the main page 1. Overview 2. What you need to know about C++ 3. Structure of the C++ template code 4. Arrays 5. Splines 6. Solving ODEs 7. Utils 8. Good coding practices 9. Installing GSL

This project is given as the course AST5220 Cosmology II at ITA at the University of Oslo. It was first given, and designed, by Hans Kristian Eriksen (who, due to his old age, of course wrote the templates in Fortran). In this updated version we also provide a template in C++. You can do the coding in whatever language you want, but we will only offer templates in C++ and Fortran. For each module there are plots you can compare your results to check that you are doing it correctly.

Can I do it in Python? You can, but I strongly suggest you don't. Python is great, but it is very useful to not just use that and get some exerience with a more low-level (and strongly typed) language. This will make you a better programmer. Many important codes you might encounter in cosmology are written in such languages like C or Fortran so its very useful to know such a language. Python is not known for speed and doing this project in Python will very likely lead to a very slow code (which means it will be hard to debug) unless you really know what you are doing. That said, if you are good, then it is very possible to write this code in Python and have the whole thing running in a few seconds (though for this you will need to use Cython or similar to speed up the slow parts).

You can download the code template from Github. Once you feel like you are ready to start coding then go to Milestone I.

C++ is a very rich and powerful language. However you do not need to be an expert to take this course. We will not use many fancy features. Just knowing the basics like for-loops, if-statements and how a simple class works is enough to get started. If you know basic C then there are just a few new features you will need to learn (plus a few things that are normal in C, but not in C++ - like the way memory should be handled). If you come from Python then many things look similar, just with a slightly different syntax, but many things are fundamentally different (you will need to get used to being a bit more precise about what types variables have, how you pass them around, differences beetween references, pointers and values) so I suggest reading a simple introduction to C++ first. You will need a C++11 compatible compiler to compile the template and below we list some of the basic features that you need to know about.

Types: C++ is a strongly typed language so you always need to specify which type your variables have. This is very different from languages like Python or Javascript where the types are implicit.

double x = 1.0; // Make a double precision floating point number
int n = 1; // Declare an integer
std::vector x(10); // Make an array with 10 elements of doubles

It's very important to note that 1 is not the same as 1.0. The former is an integer and the latter is a floating point number. C++ allows for implicit conversions whenever you use integers in floating points expressions, but there are some issues related to this. Consider the following code, what is x here?

double x = 1/3 * 3.0;

This evaluates to $x=0.0$ not $x=1.0$. The reason is that when parsing this it first enters an integer 1 then divide and then an integer 3. It therefore does integer division which gives 0 which multiplied by $9.0$ gives $0.0$. Always add the dot-zero, i.e. write $1.0/3.0$ and not $1/3$, if you are dealing with integers in floating point expressions to avoid getting integer division. If you have integer variables that you use in a floating point expression then cast them to double, i.e. double(variable).

auto: We will use the 'type' [auto] many places. This is just a simplified way of not having to write the full type of the things we are working with and let the compiler deduce the type for us (kind of like Python). For example [auto x_array = Utils::linspace(...)] means the compiler deduces the type of x_array from the signature of linspace which will be Vector of doubles. Likewise auto x = 1.0 the compiler deduces that since 1.0 is a double precision number the type of x will be double. This is just to avoid writing long names of types so if this confuses stick with using explicit types. For more info see this.

std::vector<double> give_me_a_vector(int n) {
  return std::vector<double>(n);
}

// The type of x will be deduced by the compiler to be std::vector<double>
auto x = give_me_a_vector(10); 

// ...but we could just as well have written it as
std::vector<double> x = give_me_a_vector(10);

Aliases: We will use a few aliases throughout the code. Instead of having to type std::vector<double> x; whenever you need to make a vector we simply define an alias Vector = std::vector<double> and write Vector x; Similary we will use the alias Vector2D = std::vector<Vector>. For more info see this or this.

using Vector = std::vector<double>;

// This is now the same as writing std::vector<double> x(3);
Vector x(3);

Lambda functions: We will use lambda functions a few places - anonymous functions you can define wherever you want and easily pass around and execute. If you come from Python then you might be familiar with this and it works in a similar way. For example for the ODE solver we supply the right hand side of the ODE system as a lambda function (or strictly speaking an ODEFunction = std::function<int(double,const double*, double*)>, but for our purposes its the same). For more see this, this or this.

int n = 10;

// Make a lambda function which needs n. Capture this value by reference [&]
// (We can also capture it by value [=]  so if n changes the function will not change)
auto scale_by_n = [&](double x){ return n * x; };

// This now gives us 5*n
int m = scale_by_n(5);

Focus on writing simple readable code, study the examples below for how to use the tools (see also Examples.cpp), look through the templates and come and ask me if there is something you don't understand how works and you will be fine.

The code runs from the file Main.cpp. It creates the objects we aree to implement and calls the solve routine for each object. There are four milstones corresponding to four classes we need to make:

BackgroundCosmology: This class deals with anything related to the background. It keeps the (background) cosmology parameters, has functions to get the Hubble function, density parameters as function of time and the conformal time. The only real work here is to solve the ODE for the conformal time and spline this up.

RecombinationHistory: This class contains everything related to the recombination (and reionization) history of the Universe. It solves for the free electron fraction $X_e$ by solving the Saha and Peebles equation and then integrates up the optical depth $\tau$ and provides these quantities (plus the visibility function) as function of time.

Perturbations: This class integrates the perturbations in the metric and in each of the energy components from eearly times till today. This means solving a system of $10-20$ coupled differential equations and use the result to compute source functions that are the input for the line-of-sight integration.

PowerSpectrum: This class does the line-of-sight integration of the source-function(s) to compute $\theta_\ell(k)$ (and similar for polarization) which then again is integrated to give the CMB power-spectrum. We also compute the matter power-spectrum here.

In addition to this we have some useful libraries that will come in handy namely Spline, ODESolver and Utils. More information about these can be found below. If you want to read more about what you are to implement read this paper (and then read it again). It goes through all the different parts we are going to do and will be our main reference for the numerical part of the course.

You will need good a lot of arrays (and arrays of arrays) to store the data. Make sure you know how to work with these: how to make them, how to do operations on them. One thing you should never do is to allocate and handle the memory yourself, that is stupid and will lead to memory leaks. Use a standard container like std::vector (or std::array if you want things allocated on the stack).

Often you will have to make an array with a sequence of points between a start and an end point. For this you have to choose a spacing and what type of spacing you want depends on the problem at hand. Two very common choices that cover most use-cases are linear and logarithmically spaced arrays. In the Utils namespace we provide a function that works just as in Python called linspace.

To make a two-dimensional array (a matrix) you can write Vector2D two_dim_array(nx, Vector(ny,0.0)); This can then be accessed as two_dim_array[ix][iy] for $0\leq ix \leq nx$ and $0\leq iy \leq ny$. Another common way of having a 2D array is to make a 1D vector of size nx*ny and associate the $(ix\cdot ny + iy)$ element of this array with the matrix element $(ix,iy)$. This latter thing require book-keeping (that is easy to get wrong) so stick with Vector2D.

#include <vector>
#include "Utils.h"

// In the code we have defined these two aliases 
// to avoid having to write the long names
using Vector   = std::vector<double>; 
using Vector2D = std::vector<Vector>; 

int main(){

  // Make a vector that is 100 elements long and init it to 0.0
  Vector test_1D_array(100, 0.0);

  // Make a 2D vector (matrix) that has 10 x 20 elements and init it to 0.0
  Vector2D test_2D_array(10, Vector(20, 0.0));

  // Make a lineary spaced array (here 1,2,3,...,100)
  double x_start = 1.0, x_end = 100.0; int nx = 100;
  auto x_array_lin = Utils::linspace(x_start, x_end, nx);

  // Make a log spaced array (here 1.0,10,100.0)
  double y_start = 1.0, y_end = 100.0; int ny = 3;
  auto y_array_log = exp(Utils::linspace(log(y_start), log(y_end), ny));

}

Bounds checks: Always have bounds-checks on! The most common mistake with arrays is to try to access them out of bounds (for example if you have an array x that is 5 element long and you try to write or read from x[5] or x[752] then anything can happen). In the Makefile I have added the define _GLIBCXX_DEBUG This turns on bounds-checking on std::vector. It does not slow the code down much so always keep it on. Only when you are really really sure that things are working and you want to optimize the code to try to get it to run faster can you consider turning this off.

Taking a set of discrete points $(x,y)$ and making that into a function that allows us to evaluate it at any point inbetween is something we will need frequently. Such an interpolation is called a spline. There are many different types of splines and the best choice depends on how the function you want to interpolate behaves and the accuracy you need. The simplest choice is just using the nearest $x$-point to the ones you want and the next simplest is a linear interpolation between neighboring points. We will here mainly use a so-called cubic spline (see e.g. Section 3.3. of Numerical Recipies for the algorithm). This is exact if the function we interpolate is a polynomial of degree $3$ or less and is accurate enough for most purposes. In C++ we can use the GSL library (or if you want you can write it from scratch using algorithms outlined in Numerical Recipies). We provide a simple wrapper so that you don't have to deal with memory management of the spline data for 1D and 2D splines.

Below we show an example for how to make a 1D and 2D spline.

#include <cmath>
#include "Spline.h"
#include "Utils.h"

void make_1D_spline(){

  // Make a x-array
  double x_start = 0.0, x_end = 2.0; int nx = 10;
  auto x_array = linspace(x_start, x_end, nx);

  // Fill in the y = f(x) array
  auto y_array = exp(x_array);

  // Spline it up and test it
  Spline f(x_array, y_array);
  std::cout << "For f(x) = e^x we have f(x = 1.0) = " << f(1.0) << "\n";
}

void make_2D_spline(){

  // Make a x-array
  double x_start = 0.0, x_end = 2.0; int nx = 10;
  auto x_array = Utils::linspace(x_start, x_end, nx);

  // Make a y-array
  double y_start = 0.0, y_end = 2.0; int ny = 20;
  auto y_array = Utils::linspace(y_start, y_end, ny);

  // Fill inn the z = f(x,y) = e^x + sin(y) array
  Vector2D z_array(nx, Vector(ny));
  for(size_t ix = 0; ix < x_array.size(); ix++){ 
    for(size_t iy = 0; iy < y_array.size(); iy++) 
      double x = x_array[ix];
      double y = y_array[iy];
      z_array[ix][iy] = exp(x) + sin(y);
    }
  }

  // Spline it up and test it
  Spline2D f(x_array, y_array, z_array);
  std::cout << "For f(x,y) = e^x + sin(y) we have f(x = 1.0, y = 0.5) = " << f(1.0,0.5) << "\n";
}

int main(){
  make_1D_spline();
  make_2D_spline();
}

It is also possible to create a 2D spline from a 1D (flattened) array where the (ix,iy) component located at index (ix + nx * iy). See Spline.h for more info.

Out of bounds: If you try to get the value of the function outside of the domain you created it from then it will by default return the closest value to the domain, i.e. if you created it with data on (xmin,xmax) and try to access f(xmax+1) then it will return f(xmax). You can get a warning if this happens by calling myspline.set_out_of_bounds_warning(true). To always have this on compile the code with the define _SPLINE_WARNINGS_ON. Reccomend you always have this on!

Boundary conditions: Note that the spline is created using so-called natural boundary conditions. This means that it assumes that the second derivative of the function you want to spline vanishes at the end-points. If this turns out to be a bad approximation for your function then the resulting spline will be a little inaccurate very close to the end-points. This is usually not a big issue (and if it ever turns out to be a problem we would just sample the function we are to spline on a slightly larger interval than the range we need), but its good to be aware of this.

How spline lookup is done: Whenever you try to lookup a value in a spline what it does it to do a binary search for the index $i$ in the supplied x array such that $x[i] \lt x \lt x[i+1]$ and from this evaluates f(x). The binary search is the slow part. The most common case of using splines is to evaluate the function in order (for example $f(1.0)$ then $f(1.01)$ then $f(1.02)$ instead of $f(1.0)$ then $f(20.0)$ then $f(-30.0)$ and so on). What GSL does to make the spline lookup faster is that it caches (stores) the last index it looked up. Then when it gets a new $x$ value it checks if the stored index satisfies $x[i] \lt x \lt x[i+1]$. If it does then it does not have to do the expensive binary search. Be aware that because of this it will be much faster to evaluate splines in a regular fashion than doing random lookups.

Spline with threads: The Spline class is made to be safe to use with OpenMP threads (each thread gets their own GSL accelerator so spline lookups will still be fast and there is no problem with collisions.).

Derivatives: From a cubic spline we can also get the first two dervatives (or x, y, xx, yy and xy for the case of a 2D spline) directly from the spline.

Spline f;
Spline2D g;

// ... create the splines

// Evaluate derivatives: f'(x) and f''(x)
f.deriv_x(x);
f.deriv_xx(x);

// Evaluate derivatives: g_x(x,y), g_y(x,y), g_xx(x,y), g_xy(x,y), g_yy(x,y)
g.deriv_x(x,y);
g.deriv_y(x,y);
g.deriv_xx(x,y);
g.deriv_xy(x,y);
g.deriv_yy(x,y);

For more info on how to use the Spline class see Examples.cpp or read the source Spline.h

Systems of ordinary differential equations is something that comes up a lot in physics, and it will be crucial to be able to solve such equations in this project. The equations of most interest to us is initial value problems on the form $\dot{\vec{y}} = \vec{f}(x,y)$ with some initial value $\vec{y}(x_{\rm ini}) = \vec{y}_{\rm ini}$. All ordinary differential equations can be reduced to this form. For example if we want to solve the second order ODE $y''(x) + \omega^2 y(x) = 0$ then by defining $y_0 = y$ and $y_1 = y_0'$ the system reduces to the two coupled equations $y_0' = y_1$ and $y_1' = -w^2 y_0$. On vector form with $\vec{y} = \pmatrix{y_0\\y_1}$ this can be written $\vec{y}'(x) = \vec{f}(x,\vec{y})$ with $\vec{f}(x,\vec{y}) = \pmatrix{y_1\\-w^2y_0}$. This type of trick works for higher orders as well (e.g. $y''' = x$ can be written as $y_0' = y_1$, $y_1' = y_2$, $y_2' = x$ where $y_0 = y$ is the function we are interested in). There are many methods for solving such equations (Eulers methods, Runge-Kutta, etc.) and we won't go into detail on exactly how this is done, but if you want to try to implement a solver yourself then see for example "Numerical Recipies" for an algorithm. In C++ the GSL library contains a decent solver which we provide a simple wrapper for. Below are some exampes:

#include <cmath>
#include <iostream>
#include <vector>
#include "ODESolver.h"
#include "Utils.h"

int main(){
  ODESolver ode;

  // Set up x-range we want to solve the ODE over
  double x_start = 0.0;
  double x_end   = 1.0;
  int n          = 100;
  Vector x_array = Utils::linspace(x_start, x_end, n);

  // The right hand side of the ODE system y' = y + 2*x with y(0) = 1 i.e. y(x) = e^x + x^2
  ODEFunction dydx = [&](double x, const double *y, double *dydx){
    // We only have 1 equation so y[0] is our y and dydx[0] is our y' that we need to define
    dydx[0] =  y[0] + 2*x;
    return GSL_SUCCESS;
  };

  // Initial condition: y0 = 1.0
  Vector y_ini{1.0};

  // Solve the ODE from x_start till x_end and store the solution at all the points in x_array
  ode.solve(dydx, x_array, y_ini);

  // Extract the solution y(x)
  // Many other methods to fetch the data, see ODESolver.h
  auto y_array = sol.get_data_by_component(0);

  // Output result together with analytical solution (x, y, yanal)
  for(int i = 0; i < n; i++)
    std::cout << x_array[i] << " " << y_array[i] << " " << exp(x_array[i]) + x_array[i]*x_array[i] << "\n";
}

#include <cmath>
#include <iostream>
#include <vector>
#include "ODESolver.h"
#include "Utils.h"

int main(){
  ODESolver ode;

  // Set up x-range we want to solve the ODE over
  double x_start = 0.0;
  double x_end   = 1.0;
  int n          = 100;
  Vector x_array = Utils::linspace(x_start, x_end, n);

  // The right hand side of the ODE system y0' = y1 ;  y1' = -y0
  ODEFunction dydx = [&](double x, const double *y, double *dydx){
    dydx[0] =  y[1];
    dydx[1] = -y[0];
    return GSL_SUCCESS;
  };

  // Initial condition: y0 = 1.0 ; y1 = 0.0
  Vector y_ini{1.0, 0.0};

  // Solve the ODE from x_start till x_end and store the solution at all the points in x_array
  ode.solve(dydx, x_array, y_ini);

  // Extract the components y0(x) and y1(x)
  // Many other methods to fetch the data, see ODESolver.h
  auto y0_array = sol.get_data_by_component(0);
  auto y1_array = sol.get_data_by_component(1);

  // Output result together with analytical solution (x, y0, y0anal)
  for(int i = 0; i < n; i++)
    std::cout << x_array[i] << " " << y0_array[i] << " " << y_ini[0] * cos(x_array[i]) + y_ini[1] * sin(x_array[i]) << "\n";
}

For more info on how to use the ODESolver see Examples.cpp or read the source ODESolver.h

In addition to Splines and the ODE solver I have added some useful tools in the Utils namespace. See Utils.h for what is availiable and Utils.cpp for the implementation.

Mathematical function of vectors: You can apply standard mathematical functions to a vector directly

auto exp_of_x_array = exp(x_array);
auto log_of_x_array = log(x_array);
auto sin_of_x_array = sin(x_array);

Linspace: To make a linear spaced array with n points from x_start till x_end use

auto x_array = Utils::linspace(x_start, x_end, n);

If you want logarithimically spaced arrays then you can write (for this last case we must have x_start $> 0$)

auto x_array = exp(Utils::linspace(log(x_start), log(x_end), n));

Timing: To time the code you can add Utils::StartTiming("mylabel"); before whatever you are to time and then add Utils::EndTiming("mylabel"); right after and it will print out the time used.

Root finding: If you need to find the solution to f(x) = y of a function f that you have as a spline then call Utils::binary_search_for_value(f, y, {xmin, xmax}); where (xmin,xmax) is the search interval ($f$ must lie on opposite sides of $y$ on these two points to guarantee that the method will succeed). If there are no solutions in this interval then it will throw an exception. If you omit the search interval then it will use the whole domain that f was splined over.

Bessel functions: To evaluate spherical Bessel functions $j_n(x)$ call Utils::j_ell(n, x). The standard is to use GSL's implementation here (which does not work that well for large x and n). If you have problems with GSL errors download the Complex Bessel library and compile with the define _COMPLEX_BESSEL included in the Makefile).

Units: In some of the milestones we (unfortunately) have to deal with physical units. I have added a (global) struct Constants containing all the physical constants needed in this project, see Utils.h

Here are some general suggestions on things to do and not to do to improve your code:

Compile and fix bugs often: When developing the code try to compile often and correct the things the compiler shows as errors right away. Writing huge chunks of code before compiling will often lead to a large number of errors. If you compile often you will (hopefully) have a much smaller list of errors to correct. And do fix the errors right away - they won't magically go away.

Naming: Use proper names for variables and functions, don't just call variables a,b,c,d,e,f,g,h,i,j,k or similar. That makes the code unreadable. Anybody should be able to read the code and understand for the code itself what it is doing without having to decode it. There are exceptions: its OK to use single variable names for the few things that we call by a single variable like $z$ for redshift, $c$ for the speed of light, $a$ for the scale-factor, (and $x$ for $\log(a)$ in this project), but that is about it.

// Bad naming
double b,c,d,e,x1,x2,...,jalla,test1,xxxx;
double f(){ ... }
double g(){ ... }

// Good naming
double optical_depth = 0.0, density_matter = 1.0;
double calc_conformal_time(){ ... }
double solve_for_electron_density(){ ... }

  

Comments: Document your code with comments where nessesary. For "standard" things you don't need comments as long as long as you use proper names, however if you do something special like applying an algorithm etc. then add comments telling the reader (and yourself in a few months) exactly what you are doing.

Structure: Avoid methods with hundreds of lines of code. Try to factor out larger things into smaller functions (just don't take this too far, not everything has to be its own function) that does concrete work (and give the function a good descriptive name). This will make the code easier to read and work with. Group things that belong together into classes or modules (but don't make everything a class, that defeats the purpose).

Don't explicitly allocate memory: If you need to make pointers (using 'new'): use unique or shared pointers. These work just like normal pointers however they are automatically freed when they go out of scope (are no longer used) so you avoid memory leaks. Don't explicitly allocate memory, it's a bad idea and it will lead to memory leaks and tons of time of debugging. And speaking of debugging: add the compile flag [ -fsanitize=address ] while developing. It might slow things down, but it will check and tell you at runtime if you try to access values out of scope or use things after freeing them etc. which are really useful to find mistakes.

double *x = new double[42]; // Bad way of dynamically allocating memory
delete[] x; //...because its super easy to forget this
x[43] = 752.0; //...and its easy to accidentally access it out of bounds

unique_ptr<double> x = make_unique<double>(42); // Better. Its freed automatically 
//...but still not a good idea to make an array like we do here. Contains no info about the size.

std::vector<double> x(100, 0.0); // Preferred way, freed automatically.
x[43] = 752.0; //...and we can get an error when we access it out of bounds (without slowing it down much)
  

Make reusable code if possible: If you need to do implement something that is generic (for example finding the root of a function) then write this as an independent method (i.e. without any dependencies of the code you are writing). This way you can directly (or with minor adjustments) reuse it for other things later on. That being said If you need to do something generic then there is a good chance there exists a good library that does this for you. For example you should not have to implement things like a ODE solver (unless you have a very good reason for doing this), people have done this before you and they have done it better so take advantage of this. Its a good idea to do a search to see if things exists before trying to implement it yourself.

Avoid too many globals: Global variables might seem handy, but try to avoid it. They can be accessed from any file and therefore be modified from anywhere. It makes it very hard to deal with and often lead to bugs. Instead aim to store variables into classes and modules where they belong. If you use global variables then try to make most of them const so that they don't get changed by accident in some other place.

Const: For variables you are just going to read from and not modify make them const. This avoids you accidentally modifying things that should not be modified. It also makes the code easier to read - the reader can more easily see what is being changed and what is not.

const int n = 100; // Tells us that in this function n is always 100
n = 10;            // Error: you are not allowed to change something that is const
  

The code template requires the GNU scientific library (GSL) to be installed. See this page for how to install it on a Windows machine. On Linux or a Mac you can either use a package manager or install it directly as follows:

• Go the the home directory:
  cd `$`HOME
• Make a local folder to hold libraries:
  mkdir local
• Enter this directory:
  cd local
• Download the code (if you don't have wget you need to get the file to this dir by other means):
  wget ftp://ftp.gnu.org/gnu/gsl/gsl-2.6.tar.gz
• Untar the code:
  tar -xvf gsl-2.6.tar.gz
• You should now have the gsl-2.6 folder. Enter it:
  cd gsl-2.6
• Run the configure script specifying where to install it:
  ./configure --prefix=`$`HOME/local
• Compile and install it:
  make ; make install
• In the CMB code Makefile change the include and lib paths to point to where the library was installed:
  INC=-I$\text{\$}$(HOME)/local/include
  LIBS=-L$\text{\$}$(HOME)/local/lib -lgsl -lgslcblas
• If this fails with "libgsl.so not found" then run the command:
  export LD_LIBRARY_PATH="$\text{\$}$LD_LIBRARY_PATH:$\text{\$}$HOME/local/lib"
  and try to run ./cmb again and it should work. To avoid having to run this command every time you open a new terminal open the `$`HOME/.bashrc file and add this line to the end of the file and it will load everytime you open a new window.