Python scripts for plotting, analysing and grids of runs

See the ./python directory for scripts, which require Python 2.7+ or 3.4+. To configure your system to work with the script libraries add CosmoMC's python path to your environment variables, e.g. in ~/.bashrc (or ~/bash_profile, etc):
where COSOMC..PATH is the full path of wherever you installed CosmoMC. If you have problems on a Mac or need to install python, see Installing Python.

Scripts for plotting and analysing are described below. See also grid scripts and GetDist GUI documentation. For plotting from Planck chains see the Planck readme for how to download and install.

Plotting scripts

This page gives a few simple examples to get you started with plotting from Planck of CosmoMC chains. For more detailed documentation of the full plotting and sample analysis API see the
GetDist documentation.

The getdist.plots module (directly, or via planckStyle) is used to make plots from chain results using your own custom scripts.

The sample scripts (under batch2/outputs) make use of planckStyle. This can be interchanged with getdist.plots if you do not want to use the Planck style (or some extra functions used in the batch2/output samples). Both implement the functions

The main plotting class is defined in python/getdist/ For example you can do
import planckStyle

g = planckStyle.getSinglePlotter(chain_dir = './PLA')

roots = ['base_nnu_mnu_plikHM_TT_lowTEB', 'base_nnu_mnu_plikHM_TT_lowTEB_lensing', 'base_nnu_mnu_plikHM_TT_lowTEB_lensing_post_BAO']
g.plot_2d(roots, 'mnu', 'nnu', filled=True)
g.add_legend(['Planck', 'Planck+lensing', 'Planck+lensing+BAO'], legend_loc='upper right');
The chain_dir argument can be neglected if you have set up a default in config.ini. You can do "from pylab import *" and use standard matplotlib commands to customize the plots, and there are also many plot script options you can use to customize settings, and also change global getdist plot settings using For example:
import planckStyle

g = planckStyle.getSinglePlotter(chain_dir = './PLA', ratio=1)

roots = ['base_nnu_mnu_plikHM_TT_lowTEB', 'base_nnu_mnu_plikHM_TT_lowTEB_lensing', 'base_nnu_mnu_plikHM_TT_lowTEB_lensing_post_BAO']
g.settings.solid_contour_palefactor = 0.8
g.plot_2d(roots, 'mnu', 'nnu', filled=True, colors=['orange', 'darkred', 'green'], lims=[0, 1.3, 2, 4.2])
g.add_legend(['Planck', 'Planck+lensing', 'Planck+lensing+BAO'], legend_loc='upper right', colored_text=True);
Outputs of the two versions should look like this:

See sample scripts in batch2/outputs, the Plot Gallery and tutorial, and the full GetDist documentation. Note that if you have installed CosmoMC, you don't need to separately install the GetDist python package.

Analysis scripts

The GetDist program can be used to get means, variances, limits etc from all parameters in a chain. The python scripts allow you to do this dynamically, and also offer additional features such as being able to define new derived parameters.

Calculating derived parameters

For simple calculations like finding the mean and variance of new derived parameters you can use the functions in python/ For example, if you want to calculate the posterior mean and limits for σ8 Ωm0.6 from Planck you could write a python script
import getdist.plots as gplot
g = gplot.getSinglePlotter(chain_dir=r'./PLA')

samples = g.sampleAnalyser.samplesForRoot('base_plikHM_TT_lowTEB_lensing')

p = samples.getParams()

derived = p.sigma8 * p.omegam ** 0.6

print 'mean, err = ', samples.mean(derived), samples.std(derived)
print '95% limits: ', samples.twoTailLimits(derived, 0.95)
Here p.omegam is a vector of parameter values, similar p.sigma8; samples.mean and samples.std sum the samples with the corresponding weights to calculate the result.

See python/getdist/ and python/getdist/ for other functions you can use. Use GetDist or if you want to reproduce Planck results. The mcsamples module gives you code access to most getdist results, and an MCSamples instance can be obtained from a grid as above, or you can load a chain file directly.

For example if you want to do a power law fit in the variable Ωm and H0, you could do

import GetDistPlots as s

g = s.getSinglePlotter(chain_dir=r'./PLA')
samples = g.sampleAnalyser.samplesForRoot('base_plikHM_TT_lowTEB')

print samples.PCA(['omegam', 'H0'], 'LL', 'omegam')
which fits Ωm and H0, using log transforms (L), normalized so the exponent of Ωm is unity (as in GetDist PCA outputs). The output includes
Principle components
PC1 (e-value: 0.006151)
[0.042284]  (\Omega_m/0.314833)^{1.000000}
[0.042284]  (H_0/67.300590)^{2.971216}
          = 1.000011 +- 0.004690
ND limits:     0.984    1.017    0.983    1.019
which tells you that ΩmH02.97is constrained at the 0.4% level.

Samples can also be loaded directly from single chains, optionally with custom settings, e.g.

import getdist

samples = getdist.loadMCSamples(r'./PLA/base/plik_HM_TT_lowTEB/base_plikHM_TT_lowTEB', settings={'ignore_rows':0.3})

Adding and plotting new derived parameters

You can use the addDerived function of the MCSamples class to add a new parameter to existing chains, which can then be used like any of the original parameters. An example gives the idea (from batch2/outputs/, which adds a new parameter with name tag 'rsH' which is a function of the original parameters 'Hubble057' and 'rdrag':
p = samples.getParams()
rd_fid = 149.28
rsH = p.Hubble057 * p.rdrag / rd_fid
samples.addDerived(rsH, name='rsH', label=r'H(0.57) (r_{\mathrm{drag}}/r_{\mathrm{drag}}^{\rm fid})\, [{\rm km} \,{\rm s}^{-1}{\rm Mpc}^{-1}]')
The last line is needed to ensure caches etc are cleared, but only needs to be called once if you are adding several new parameters. After doing this, you can use 'rsH' as you would any of the original parameter names in the chain.

Analysing samples directly

Rather than reading chain files, you can also analyse samples directly, for example using the constructor
 from getdist import mcsamples
 samples = mcsamples.MCSamples(samples=sample_points, loglikes=loglikes, names=names)
where sample_points is a matrix of sample values, names are the parameter names (list of strings), and loglikes is (optionally) an array of corresponding -log(likelihood) values.

Plotting densities

The g.plot_2d (etc) functions will calculate and plot densities from samples, but you can also add densities that are calculated yourself for comparison, e.g.:
from getdist.densities import Density2D
import getdist.plots as gplot
import numpy as np

g = gplot.getSinglePlotter(chain_dir=r'./PLA')

xvalues = np.arange(85, 110, 0.3)
yvalues = np.arange(1000, 1500, 4)
x,y = np.meshgrid(xvalues, yvalues)
loglike = my_loglike_func(x,y)
density = Density2D(xvalues,yvalues, np.exp(-loglike / 2))
density.contours = np.exp(-np.array([1.509, 2.4477]) ** 2 / 2)
g.add_2d_contours(root, 'x', 'y', filled=True, density=density)

Other scripts

Multi-purpose use and utility cosmology scripts include

Installing Python

Instructions below are for installing python to run things natively on your machine. You could also use CosmoBox.

Python on Windows

See the Python 2.7 download, or install a package like Anaconda (Python 2.7 or 3.4+). Also install matplotlib and PySide standard packages using standard methods, or download builds from here.

Python on a Mac

There are several options, here's one: Check that python --version says 2.7.9. If not, edit your system path to remove other python installations.

Writing Python

While not needed for GetDist/Cosmomc scripts to run, you may also find an integrated development environment like PyDev or PyCharm useful (they are similar). This will let you run in one click, check syntax as you type, format nicely, etc.