MCPlusPlus
A Monte Carlo C++ code for radiative transport
MCPlusPlus Documentation

Table of Contents

Introduction

Numerous areas of science employ Monte Carlo (MC) methods to simulate complex processes. Light propagation in random media, often referred to as photon migration, is an area of physics in which such methods are of great importance. Prime examples include radiative transfer in highly scattering materials such as biological tissues, clouds, and pharmaceuticals. There, MC simulation is generally considered the gold standard of modeling and is used to investigate complex systems and processes, to validate simpler models, as well as to evaluate data.

Given the complexity of the radiative transfer equation (RTE) – a widely used analytical model for light diffusion through turbid media – approximations such as the diffusion approximation (DA) are most often used. However these simplifications do not perform equally well in every situation and they have their shortcomings. On the other hand the RTE directly follows from energy conservation. Its derivation is phenomenological, where the complex wave propagation is reduced to a Poissonian random walk with exponentially distributed steps between scattering events, with correlations in directionality given by the phase function. With the computing power available nowadays, it is possible to implement a Monte Carlo method to find exact solutions to the RTE, the only approximation being that the number of simulated photons is finite. This is especially true because this kind of simulations have an inherently parallel nature and so they can leverage nowadays multi-threading CPUs and GPUs.

Some excerpts rearranged from [3] [1].

Motivation

MCPlusPlus is a Monte Carlo implementation of photon transport in turbid layered media. As its name suggests, MCPlusPlus is entirely written in C++. We think that the Object-oriented programming (OOP) paradigm and the abstraction layer that it provides apply very well to the problem described above. OOP allows to encapsulate concepts as reusable objects, thus providing ease of maintenance, increased code understandability and modular architecture. With this powerful tool, we can implement a software package that is much more flexible than the most widely used MC codes, yet having equal to superior performance. The OOP makes it easy to expand and customize the code, that's why we can provide some unique features such as:

Download

mcplusplus_2016.03.04_amd64.deb (211 KiB)

mcplusplus-2016.03.04.tar.gz (183 KiB)

License

This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Authors: Giacomo Mazzamuto and Lorenzo Pattelli, European Laboratory for Non-Linear Spectroscopy (LENS)

How to install

Extract the source tarball and create a new folder that will be used for compilation. This is called the build folder. Open a terminal in the build folder and issue the following command:

$ cmake [ADDITIONAL_DEFINES] pathToSourceFolder

where pathToSourceFolder is the absolute or relative path to the extracted source folder. This will prepare the files needed for compilation.

If you wish to specify an installation path different from your system's default (for example you want to install in /home/yourname/local), specify the following define in the command above:

-DCMAKE_INSTALL_PREFIX=/home/yourname/local

If you have libraries and headers installed in non-standard paths, for example in /home/yourname/local/include or /home/yourname/local/lib, you can add those paths to the cmake command line like this:

-DCMAKE_LIBRARY_PATH=/home/yourname/local/lib -DCMAKE_CXX_FLAGS=-I/home/yourname/local/include/

Once cmake has been run, you can build the project and install it in the specified folder with the following command:

$ make install

If you installed in a local folder, in order to be able to use the software you'll have to export the following environment variables in every shell that you open:

export PATH=/home/yourname/local/bin:$PATH
export LD_LIBRARY_PATH=/home/yourname/local/lib:$LD_LIBRARY_PATH
export PYTHONPATH=/home/yourname/local/lib:$PYTHONPATH

assuming, again, your installation prefix to be /home/yourname/local/. You can also place the lines above in your shell's configuration file.

Build Options

The following options can be specified in the cmake command line preceded by -D, for example

-DBUILD_PYTHON_BINDINGS=TRUE

Boolean options:

Tutorials

These simple tutorials are meant to showcase the easy and straightforward interface provided by MCPlusPlus.

Python Tutorial

To take advantage of Python scripting MCPlusPlus is compiled by default with the -DBUILD_PYTHON_BINDINGS flag set to TRUE. Remember to set the PYTHONPATH environment variable as described above.

Create a file example.py with the contents shown below. You can also find this file in the examples folder of the source tarball.

#!/usr/bin/env python3
from pymcplusplus import *
# define a sample consisting of two layers of different materials
sample = Sample()
mat = Material()
mat.n = 1.5
mat.g = 0
mat.ls = 1000
# add a layer of material mat and thickness 1000um
sample.addLayer(mat, 1000)
mat2 = Material()
mat2.n = 1.3
mat.g = 0.5
mat.ls = 700
sample.addLayer(mat2, 500)
externalMaterial = Material()
externalMaterial.n = 1
sample.setSurroundingEnvironment(externalMaterial)
# define a photon source
source = PencilBeamSource()
# define the main simulation object, simulate 1e6 photons using 8 parallel
# threads. Note that each thread will use increasing seeds starting from 0.
sim = Simulation()
sim.setSample(sample)
sim.setSource(source)
sim.setNPhotons(1000000)
sim.setNThreads(8)
sim.setSeed(0)
sim.setOutputFileName("example.h5")
# define and add several histograms to the simulation
# 1) a histogram of the exit times
hist = Histogram()
hist.setDataDomain(DATA_TIMES)
hist.setPhotonTypeFlags(FLAG_TRANSMITTED)
hist.setMax(1000)
hist.setBinSize(2)
hist.setName("times")
hist.addMomentExponent(2) # compute spatial variance
sim.addHistogram(hist)
# 2) a histogram of the exit distances from the center
hist = Histogram()
hist.setDataDomain(DATA_POINTS)
hist.setPhotonTypeFlags(FLAG_TRANSMITTED)
hist.setMax(1000000)
hist.setBinSize(50)
hist.setName("points")
sim.addHistogram(hist)
# 3) a bivariate histogram of the exit distances as a function of time
hist = Histogram()
hist.setDataDomain(DATA_POINTS, DATA_TIMES)
hist.setPhotonTypeFlags(FLAG_TRANSMITTED)
hist.setMax(100000, 1000)
hist.setBinSize(500, 2)
hist.setName("points_vs_times")
sim.addHistogram(hist)
sim.run()

Once you have created it, you probably need to grant it execution permissions using

$ chmod +x example.py

Then run the simulation by simply issuing

$ ./example.py

Here is an example of what the output may look like:

[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 0
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 1
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 2
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 3
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 7
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 5
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 6
[2015-02-19 17:17:03][MCPP::Simulation] starting... Number of walkers = 125000000, original seed = 4
[2015-02-19 17:17:04][MCPP::Simulation]
================
transmitted: 454581567
ballistic: 225229317
reflected: 278973129
back-reflected: 41231183
Completed in 6871 seconds
================

While the simulation runs, you can query its progress by sending a USR2 signal to the process. To do this, open another shell and issue:

$ killall -USR2 example.py

or

$ kill -USR2 [pid of your process]

Once the simulation has ended, the histogrammed data is saved in an H5 file named example.h5 as was specified in the script. This file contains several datasets, you can list them with the command

$ h5ls example.h5
photon-counters Dataset {4}
points Dataset {20001, 2}
points_vs_times Dataset {201, 502}
sample Dataset {4, 4/Inf}
times Dataset {501, 3}

Notice the three datasets containg the histograms that were defined in the Python script: points, times and points_vs_times. You can inspect the file contents using a graphical viewer:

$ hdfview example.h5

or you can extract single datasets as text files:

$ h5totxt -d [dataset name] example.h5 > dataset.txt

In MATLAB you can directly load a dataset with:

hdf5read('example.h5', 'dataset_name')

C++ Tutorial

#include <stdio.h>
#include <MCPlusPlus/simulation.h>
using namespace MCPP;
int main(int argc, char *argv[])
{
// define a sample consisting of two layers of different materials
Sample *sample = new Sample();
Material *mat = new Material();
mat->n = 1.5;
mat->g = 0;
mat->ls = 1000;
sample->addLayer(mat,1000); //add a layer of material mat and thickness 1000um
Material *mat2 = new Material();
mat2->n = 1.3;
mat2->g = 0.5;
mat2->ls = 700;
sample->addLayer(mat2,700);
Vacuum *vacuum = new Vacuum();
sample->setSurroundingEnvironment(vacuum);
// define a photon source
Source *source = new PencilBeamSource(0);
// define the main simulation object, simulate 1e6 photons using 8 parallel threads.
// note that each thread will use increasing seeds starting from 0
Simulation *sim = new Simulation(0);
sim->setSample(sample);
sim->setSource(source);
sim->setNPhotons(1000000);
sim->setNThreads(8);
sim->setSeed(0);
sim->setOutputFileName("example.h5");
// define and add several histograms to the simulation
Histogram *hist;
// 1) a histogram of the exit times
hist = new Histogram(0);
hist->setDataDomain(DATA_TIMES);
hist->setPhotonTypeFlags(FLAG_TRANSMITTED);
hist->setMax(1000);
hist->setBinSize(2);
hist->setName("times");
hist->addMomentExponent(2); //compute spatial variance
sim->addHistogram(hist);
// 2) a histogram of the exit distances from the center
hist = new Histogram(0);
hist->setDataDomain(DATA_POINTS);
hist->setPhotonTypeFlags(FLAG_TRANSMITTED);
hist->setMax(1000000);
hist->setBinSize(50);
hist->setName("points");
sim->addHistogram(hist);
// 3) a bivariate histogram of the exit distances as a function of time
hist = new Histogram(0);
hist->setDataDomain(DATA_POINTS, DATA_TIMES);
hist->setPhotonTypeFlags(FLAG_TRANSMITTED);
hist->setMax(100000,1000);
hist->setBinSize(500,2);
hist->setName("points_vs_times");
sim->addHistogram(hist);
sim->run();
return 0;
}

Create a file example.py with the contents shown above. You can find this file also in the examples folder of the source tarball. Then compile it:

g++ -I/home/yourname/local/include -L/home/yourname/local/lib example.cpp -o example /home/yourname/local/lib/libMCPlusPlus.a -lhdf5 -lhdf5_cpp -lboost_thread -lboost_system

assuming again MCPlusPlus installed in /home/yourname/local/, and run it:

$ ./example

See the last part of Python Tutorial for an explanation of how to examine the output data.