Getting Started¶
Foreword¶
Python has a very nice way of handling errors in the execution. Instead of a segmentation fault as in C, when the code breaks, you have access to the whole stack of actions that lead to the error. This helps you pin-point which function was called, which line was responsible for the error.
It can however be lengthy, and to help everyone reading it, a messaging system was implemented in Monte Python. After a blank line, a summary of the actual error will be displayed. When reporting for an error, please attach the entire output, as this is priceless for debugging.
Input parameter file¶
An example of input parameter file is provided with the download
package, under the name example.param
. Input files are
organised as follows:
data.experiments = ['experiment1', 'experiment2', ...]
data.parameters['cosmo_name'] = [mean, min, max, sigma, scale, 'cosmo']
...
data.parameters['nuisance_name'] = [mean, min, max, sigma, scale, 'nuisance']
...
data.parameters['cosmo_name'] = [mean, min, max, sigma, scale, 'derived']
...
data.cosmo_arguments['cosmo_name'] = value
data.N = 10
data.write_step = 5
The first command is rather explicit. You will list there all the
experiments you want to take into account. Their name should coincide
with the name of one of the several sub-directories in the
montepython/likelihoods/
directory. Likelihoods will be explained in the
Likelihood class module
In data.parameters
, you can list all the cosmo and nuisance
parameter that you want to vary in the Markov chains. For each of them
you must give an array with six elements, in this order:
- mean value (your guess for the best fitting value, from which the first jump will start)
- minimum value (set to -1 or None for unbounded prior edge),
- maximum value (set to -1 or None for unbounded prior edge),
- sigma (your guess for the standard deviation of the posterior of this parameter, its square will be used as the variance of the proposal density when there is no covariance matrix including this parameter passed as an input),
- scale (most of the time, it will be 1, but occasionnaly you can use a rescaling factor for convenience, for instance {tt 1.e-9} if you are dealing with \(A_s\) or
0.01
if you are dealing with \(\omega_b\))- role (
cosmo
for MCMC parameters used by the Boltzmann code,nuisance
for MCMC parameters used only by the likelihoods, andderived
for parameters not directly varied by the MCMC algorithm, but to be kept in the chains for memory).
In data.cosmo_arguments
, you can pass to the Boltzmann code
any parameter that you want to fix to a non-default value
(cosmological parameter, precision parameter, flag, name of input file
needed by the Bolztmann code, etc.). The names and values should be
the same as in a Class input file, so the values can be numbers or a
strings, e.g:
data.cosmo_arguments['Y_He'] = 0.25
or
data.cosmo_arguments['Y_He'] = 'BBN'
data.cosmo_arguments['sBBN file'] = data.path['cosmo']+'/bbn/sBBN.dat'
All elements you input with a cosmo
, derived
or
cosmo_arguments
role will be interpreted by the cosmological
code (only Class so far). They are not coded anywhere inside Monte Python.
Monte Python takes parameter names, assigns values, and passes all of these to
Class as if they were written in a Class input file. The
advantages of this scheme are obvious. If you need to fix or vary
whatever parameter known by Class, you don’t need to edit Monte Python, you
only need to write these parameters in the input parameter file. Also,
Class is able to interpret input parameters from a Class input
file with a layer of simple logic, allowing to specify different
parameter combinations. Parameters passed from the parameter file of
Monte Python go through the same layer of logic.
If a cosmo
, derived
or cosmo_arguments
parameter is not understood by the Boltzmann code, Monte Python will stop
and return an explicit error message. A similar error will occur if
one of the likelihoods requires a nuisance
parameter that is
not passed in the list.
You may wish occasionally to use in the MCMC runs a new parameter
that is not a Class parameter, but can be mapped to one or
several Class parameters (e.g. you may wish to use in your chains
\(\log(10^{10}A_s)\) instead of \(A_s\)). There is a function,
in the module data
, that you can edit to define such
mappings: it is called update_cosmo_arguments
. Before calling CLASS, this
function will simply substitute in the list of arguments your
customized parameters by some Class parameters. Several exemple of
such mappings are already implemented, allowing you for instance to
use 'Omega_Lambda'
, 'ln10^{10}A_s'
or
'exp_m_2_tau_As'
in your chains. Looking at these examples,
the user can easily write new ones even without knowing python.
The last two lines of the input parameter file are the number of steps
you want your chain to contain (data.N
) and the number of
accepted steps the system should wait before writing it down to a file
(data.write_step
). Typically, you will need a rather low
number here, e.g. data.write_step = 5
or 10
. The
reason for not setting this parameter to one is just to save a bit of
time in writing on the disk.
In general, you will want to specify the number of steps in the
command line, with the option -N
(see section~ref{commands}).
This will overwrite the value passed in the input parameter file. The
value by default in the parameter file, data.N = 10
, is
intentionnaly low, simply to prevent doing any mistake while testing
the program on a cluster.
Output directory¶
You are assumed to use the code in the following way: for every set of
experiments and parameters you want to test, including different
priors, some parameters fixed, etcldots you should use one output
folder. This way, the folder will keep track of the exact calling of
the code, allowing you to reproduce the data at later times, or to
complete the existing chains. All important data are stored in your
folder/log.param
file.
Incidentaly, if you are starting the program in an existing folder,
already containing a log.param
file, then you do not even have
to specify a parameter file: the code will use it automatically. This
will avoid mixing things up. If you are using one anyway, the code
will warn you that it did not read it: it will always only use the
log.param
file.
In the folder montepython
, you can create a folder
chains
where you will organize your runs e.g. in the
following way:
montepython/chains/set_of_experiments1/model1
montepython/chains/set_of_experiments1/model2
...
montepython/chains/set_of_experiments2/model1
montepython/chains/set_of_experiments2/model2
...
The minimum amount of command lines for running Monte Python is an input file,
an output directory and a configuration file: if you have already
edited defaut.conf
or copied it to your own
my-machine.conf
, you may already try a mini-run with the
command
montepython]$ montepython/MontePython.py -conf my-machine.conf -p example.param -o test
Analyzing chains and plotting¶
Once you have accumulated a few chains, you can analyse the run to get
convergence estimates, best-fit values, minimum credible intervals, a
covariance matrix and some plots of the marginalised posterior
probability. You can run again Monte Python with the info
prefix
followed by the name of a directory or of several chains, e.g.
info chains/myrun/
or info chains/myrun/2012-10-26*
chains/myrun/2012-10-27*
. There is no need to pass an input file
with parameter names since they have all been stored in the
log.param
.
Information on the acceptance rate and minimum \(-\log{\cal
L}=\chi^2_{\rm eff}/2\) is written in chains/myrun/myrun.log
.
Information on the convergence (Gelman-Rubin test for each chain
parameter), on the best fit, mean and minimum credible interval for
each parameter at the 68.26%, 95.4%, 99.7% level are written in
horizontal presentation in chains/myrun/myrun.h_info
, and in
vertical presentation in chains/myrun/myrun.v_info
(without
99.7% in the vertical one). A latex file to produce a table with
parameter names, means and 68% errors in written in
chains/myrun/myrun.tex
.
The covariance matrix of the run is written in
chains/myrun/myrun.covmat
. It can be used as an input for the
proposal density in a future run. The first line, containing the
parameter name, will be read when the covariance matrix will be passed
in input. This means that the list of parameters in the input
covariance matrix and in the run don’t need to coincide: the code will
automatically eliminate, add and reorder parameters (see
mcmc.get_covariance_matrix()
). Note that the rescaling factors
passed in the input file are used internally during the run and also
in the presentation of results in the .h_info
,
.v_info
, .tex
files, but not in the covariance matrix
file, which refers to the true parameters.
The 1D posteriors and 2D posterior contours are plotted in
chains/myrun/plots/myrun_1D.pdf
and
chains/myrun/plots/myrun_triangle.pdf
. You will find in the
Parser module documentation a list of commands to customize the
plots.
When the chains are not very converged and the posterior probability
has local maxima, the code will fail to compute minimum credible
intervals and say it in a warning. The two solutions are either to
re-run and increase the number of samples, or maybe just to decrease
the number of bins with the --bins
option.
Global running strategy¶
In the current version of Monte Python, we deliberately choose not to use MPI
communication between instances of the code. Indeed the use of MPI
usually makes the installation step more complicated, and the gain is,
in our opinion, not worth it. Several chains are launched as
individual serial runs (if each instance of Monte Python is launched on
several cores, Class and the WMAP likelihood will parallelize since
they use OpenMP). They can be run with the same command since chain
names are created automatically with different numbers for each
chain: the chain names are in the form yyyy-mm-dd_N__i.txt
where yyyy
is the year, mm
the month, dd
the
day, N
the requested number of steps and i
the
smallest available integer at the time of starting a new run.
However the absence of communication between chains implies that the
proposal density cannot be updated automatically during the initial
stage of a run. Hence the usual strategy consists in launching a first
run with a poor (or no) covariance matrix, and a low acceptance rate;
then to analyze this run and produce a better covariance matrix; and
then to launch a new run with high acceptance rate, leading to nice
plots. Remember that in order to respect strictly markovianity and the
Metropolis Hastings algorithm, one should not mix up chains produced
with different covariance matrices: this is easy if one takes
advantage of the info
syntax, for example info
chains/myrun/2012-10-26_10000*
. However mixing runs that started from
very similar covariance matrices is harmless.
It is also possible to run on several desktops instead of a single
cluster. Each desktop should have a copy of the output folder and with
the same log.param
file, and after running the chains can be
grouped on a single machine and analyse. In this case, take care of
avoiding that chains are produced with the same name (easy to ensure
with either the -N
or --chain-number
options). This is
a good occasion to keep the desktops of your department finally busy.