To create a CosmoSlik script which runs a chain, you start with the following boilerplate code, which you should put in some file, e.g. myscipt.py
,
from cosmoslik import *
class myscript(SlikPlugin):
def __init__(self):
super().__init__()
# your initialization code will go here
def __call__(self):
# you'll return negative log-likelihood here
This script is like the "ini" file other similar codes use to define a particular run, but is much more flexible because its just Python code which can do arbitrarily complex things, as we'll see below.
Lets code up a simple example that samples a 2D unit Gaussian likelihood. In the initialization section we define the parameters to be sampled by the MCMC chain, and in the likelihood section we use these parameters to return a likelihood (by convention, CosmoSlik wants the negative log-likelihood). Finally, we set which MCMC sampler to use. The new script looks like this:
from cosmoslik import *
class myscript(SlikPlugin):
def __init__(self):
super().__init__()
# define sampled parameters
self.a = param(start=0, scale=1)
self.b = param(start=0, scale=1)
# set the sampler
self.sampler = samplers.metropolis_hastings(
self,
num_samples=1e5,
output_file="myscript.chain",
)
def __call__(self):
# compute the likelihood
return (self.a**2 + self.b**2)/2
You can now immediately run this chain from the command line by doing,
$ cosmoslik -n 8 myscript.py
The -n 8
option runs 8 chains in parallel with MPI (you will need mpi4py
installed) and automatically updates the proposal covariance. The resulting 8 chains are all written to the single file specified in the script, in our case myscript.chain
. You can read the chain file (including while the chain is running to see its progress) via
>>> chain = load_chain("myscript.chain")
Alternatively, you could run this script from an interactive session and get a chain directly via
>>> chain = run_chain("myscript.py", nchains=8)
Or you could have skipped the script file entirely and just run
>>> chain = run_chain(myscript, nchains=8)
assuming you defined myscript
in your interactive session.
load_chain
and run_chain
return Chain
or Chains
objects which have a number of useful methods for manipulating, post-processing, and plotting results. For example, since we ran 8 chains, we might want to concatenate them into a single chain, first burning-off some samples at the beginning of each chain, and then create a "triangle" plot from the result. This would be done with,
>>> chain.burnin(500).join().likegrid()
The power in using Python to define scripts is that we can do lots of advanced things that can't be done in "ini" files or other mini-languages. This means a single script can in fact be a powerful multi-purpose tool that runs several different chains. For example, we can decide on the fly how many parameters to sample. Consider the following script which samples an ndim
-dimensional Gaussian where ndim
is a parameter we will pass in,
from cosmoslik import *
class myscript(SlikPlugin):
def __init__(self, ndim, num_samples=1000):
super().__init__()
# save the ndim parameter so we can use it below in __call__ too
ndim = self.ndim = int(ndim)
# dynamically create the sampled parameters we need
for i in range(ndim):
self["param%i"%i] = param(start=0, scale=1)
self.sampler = samplers.metropolis_hastings(
self,
num_samples=num_samples,
output_file="myscript_ndim_%i.chain"%ndim
)
def __call__(self):
return sum([self["param%i"%i]**2 for i in range(self.ndim)])/2
Let's break down some new things we did here:
- We gave the
__init__
function some arguments,ndim
andnum_samples
- We defined the sampled parameters for this chain dynamically with a
for
loop. - We used
self
as a dictionary, i.e.self["param"]
in place ofself.param
, which allows us to create the numbered parameters
Additionally, simply by having written this script, CosmoSlik creates a nice wrapper you can use to call this script from the command line and specify parameters. You can see the "help" for it via:
$ cosmoslik myscript.py -h
usage: cosmoslik myscript.py [-h] [--num_samples NUM_SAMPLES] ndim
positional arguments:
ndim
optional arguments:
-h, --help show this help message and exit
--num_samples NUM_SAMPLES
default: 1000
You can then run e.g. a a 10-dimensional chain with 10000 steps via:
$ cosmoslik myscript.py 10 --num_samples 10000
Having seen the basic machinery of how to write CosmoSlik scripts and run them, lets see how to run a real cosmology chain. As an example, we'll run a Planck chain, using CAMB to compute the CMB Cl's. The script file looks like this,
class planck(SlikPlugin):
def __init__(self):
super().__init__()
# define sampled cosmological parameters
param = param_shortcut('start','scale')
self.cosmo = SlikDict(
logA = param(3.108, 0.03),
ns = param(0.962, 0.006),
ombh2 = param(0.02221, 0.0002),
omch2 = param(0.1203, 0.002),
theta = param(0.0104, 0.00003),
tau = param(0.055, 0.01),
)
# sample the Planck calibration as well
self.calPlanck = param(1,0.0025,gaussian_prior=(1,0.0025))
# load CAMB to compute Cls
self.camb = models.camb()
# load the Planck likelihood files
self.highlTT = likelihoods.clik(clik_file='plik_lite_v18_TT.clik')
self.lowlTT = likelihoods.clik(clik_file='commander_rc2_v1.1_l2_29_B.clik')
self.sampler = samplers.metropolis_hastings(
self,
num_samples = 1e7,
output_file = 'planck.chain',
)
def __call__(self):
# we sample logA but CAMB needs As
self.cosmo.As = exp(self.cosmo.logA)*1e-10
# compute Cls
self.cls = self.camb(**self.cosmo)
# the two Planck likelihoods read the calibration from `A_Planck` and
# `A_planck`, so set those based on our sampled parameter
self.highlTT.A_Planck = self.lowlTT.A_planck = self.calPlanck
# compute likelihood
return lsum(
lambda: self.highlTT(self.cls),
lambda: self.lowlTT(self.cls)
)
Some new things here are:
- "Attaching" sampled parameters not directly to
self
but to a sub-attribute, in this case,self.cosmo
. CosmoSlik will find all sampled parameters if they are attached to anySlikDict
s attached toself
(recursively, any number ofSlikDict
s deep). You can use this to organize parameters into convenient subgroups. - Using the
lsum
function. This is a convenience function which simply sums up its arguments in order, but if one of them returnsinf
, it doesn't waste time evaluating the rest.