CosmoSlik (Cosmology Sampler of Likelihoods) is:
SlikPlugins.SlikPlugins for analyzing cosmological datasets, such as the Planck likelihood, CAMB boltzmann solver, PICO interpolator, etc...CosmoSlik is functionally very similar to many Python MCMC samplers such as PyMC or emcee (in fact it contains a SlikPlugin to sample using emcee), as well as being similar to other cosmological samplers such as CosmoMC, MontePython, or Cosmosis. What distinguishes it is its modular structure, and the ease and power with which one can construct likelihoods, debug them, and analyze their state at any point in the chain.
In [1]:
from cosmoslik import *
As a first example, let's run a chain that samples a simple 2-d Gaussian likelihood. First I'll show you the code, then explain it:
In [2]:
class simple_gaussian(SlikPlugin):
def __init__(self):
super(SlikPlugin,self).__init__()
self.a = param(start=0,scale=1)
self.b = param(start=0,scale=1)
self.sampler = get_plugin('samplers.metropolis_hastings')(self,num_samples=10000)
def __call__(self):
return self.a**2/2 + self.b**2/2
First we create the likelihood by creating a class (simple_gaussian in this case) which is a subclass of SlikPlugin. All SlikPlugins have an __init__() function, which does any initialization, and a __call__() function, which executes the plugins main functionality. Here that functionality is returning the likelihood, but we will see in further examples __call__() functions which do other things, for example returning some model rather evaluating a likelihood.
In the initialization function, we define what parameters are going to be sampled by the MCMC sampler:
self.a = param(start=0,scale=1)
self.b = param(start=0,scale=1)
This means parameter "a" will start from 0, and be sampled with "scale"=1. The "scale" is a loosely defined quantity which tells the MCMC sampler the scale on which to vary this parameter, and may be used in slightly different ways by various samplers (for the default Metropolis-Hasting sampler, this is used as the initial proposal standard-deviation).
The other initialization step is to pick a sampler, and attach it to the likelihoood's sampler key:
self.sampler = get_plugin('samplers.metropolis_hastings')(self,num_samples=1000)
get_plugin is a convenience function to load one of the installed plugin, which reside in the cosmoslik_plugins Python package. To get a list of all available plugins you can run:
In [3]:
sorted(get_all_plugins().values())
Out[3]:
At this point its worth also noting that all SlikPlugin (such as the simple_gaussian that we created, as well as all plugins returned by get_plugin) are just slightly enhanced Python dictionaries which allow convenient access to their data via the "." notation, as well as accessing the data sub-plugins recursively. For example,
In [4]:
x = SlikPlugin()
In [5]:
x['a']=3
x.a
Out[5]:
In [6]:
x['b']=SlikPlugin(c=1)
x.b.c
Out[6]:
In [7]:
x['b.c']
Out[7]:
Now back to the simple Gaussian likelihood. The __call__() function by convention returns the negative log-likelihood. Keys which were set in __init__() as param(...) will be set to numerical values by the CosmoSlik sampler before calling the __call__() function. Thus the our simple likelihood is, in this case,
return self.a**2/2 + self.b**2/2
Running the chain is as simple as:
In [8]:
simple_gaussian_chain = run_chain(simple_gaussian)
This returns a Chain object, which is the default way CosmoSlik stores chains.
In [9]:
simple_gaussian_chain
Out[9]:
Like SlikPlugin, Chain is also a Python dictionary with some enhancements. It contains functions for plotting, analyzing, and manipulating data. Below are some examples.
In [10]:
simple_gaussian_chain.params()
Out[10]:
In [11]:
zip(simple_gaussian_chain.params(),simple_gaussian_chain.std())
Out[11]:
In [12]:
simple_gaussian_chain.acceptance()
Out[12]:
In [13]:
simple_gaussian_chain.plot()
In [14]:
simple_gaussian_chain.likegrid()
Coding up a non-trivial likelihood can often be a difficult process. One of the biggest strengths of CosmoSlik is the ability to manually step through a chain, examining and debugging the state of the SlikPlugin object at each step. To see this, let's first create a slightly more complicated likelihood, one that fits a line to some random data using least-squares.
In [15]:
class fit_line(SlikPlugin):
ndata = 30
def __init__(self):
super(SlikPlugin,self).__init__()
self.data = (3*arange(self.ndata) + 11) + randn(self.ndata)*10
self.a = param(start=0,scale=0.3)
self.b = param(start=0,scale=5)
self.sampler = get_plugin('samplers.metropolis_hastings')(self,num_samples=10000)
def get_model(self):
return self.a*arange(self.ndata) + self.b
def plot(self):
errorbar(arange(self.ndata),self.data,yerr=10*ones(self.ndata),ls='',marker='.',color='k')
plot(self.get_model())
def __call__(self):
return ((self.get_model() - self.data)**2/2/100).sum()
fit_line_chain = run_chain(fit_line)
Note that we created a separate function to get the "model" (just our straight line for the current values of a and b), as well as one for plotting. We'll see how they come in handy in a second. To get started stepping through the state of this chain, we first wrap a fit_line instance in a Slik object:
In [16]:
fit_line_instance = fit_line()
fit_line_slik=Slik(fit_line_instance)
It has a few methods, including one which automatically scans through the SlikPlugin and recursively through any sub-plugins to get an ordered list of the sampled parameters:
In [17]:
fit_line_slik.get_sampled().keys()
Out[17]:
The evalute function is what lets us call the likelihood for a given set of parameters. It returns a tuple of (lnl, instance) where lnl is the negative log likelihood as returned by the __call__ function, and instance is an instance of fit_line with the parameters evaluated at the values passed to evaluate.
In [18]:
lnl, fit_line_eval = fit_line_slik.evaluate(a=1,b=3)
In [19]:
lnl
Out[19]:
In [20]:
fit_line_eval
Out[20]:
We now call any of the functions on fit_line_eval, such as plot,
In [21]:
fit_line_eval.plot()
or something which is often useful, plot the best-fitting model from the chain:
In [22]:
fit_line_chain.best_fit()
Out[22]:
In [23]:
fit_line_slik.evaluate(**fit_line_chain.best_fit())[1].plot()
We can also manually run the first 100 steps of the chain by hand, calling plot at each step:
In [24]:
for sample,_ in zip(fit_line_slik.sample(),range(100)):
sample.extra.plot()
We can see from the plot above the chain requires some burnin this time, as we did not start at the best-fit parameters. This is also visible if we plot the first 50 accepted samples. Note that some of these samples have weight>1, hence why we see the x-axis below go above 50,
In [25]:
fit_line_chain.sample(slice(0,50)).plot()
Different sampler plugins each provide their own method of parallelization. However, there is one CosmoSlik-wide way to trivially run multiple chains in parallel using Python's multiprocessing module by passing the nchains keyword to run_chain. In this case the return value is a Chains object which is just a list of Chain objects,
In [26]:
simple_gaussian_parallel_chain = run_chain(simple_gaussian,nchains=4)
In [27]:
simple_gaussian_parallel_chain.plot()
Often you'll want to remove the first several samples from the chains ("burn in"), then join them into one,
In [28]:
simple_gaussian_parallel_chain.burnin(2000).join().likegrid()
Now for a more interesting example of parallelization with MPI. Again, each sampler plugin provides it's own method for parallelization, this example is using the default metropolis_hastings plugin. It will require both an MPI implementation and mpi4py installed on your system.
To use MPI we will have to run from outside of IPython, which means we write the simple_gaussian code to a CosmoSlik script file. This is just a Python file which conatins the likelihood plugin (you can think of this as the "ini" file some other codes use). Note we've also told the sampler to save the chain to the file output_file,
In [15]:
%%file simple_gaussian.py
from cosmoslik import *
class simple_gaussian(SlikPlugin):
def __init__(self):
super(SlikPlugin,self).__init__()
self.a = param(start=0,scale=0.1)
self.b = param(start=0,scale=10)
self.sampler = get_plugin('samplers.metropolis_hastings')(self,num_samples=10000,output_file='simple_gaussian.chain')
def __call__(self):
return self.a**2/2 + self.b**2/2
Now we simply run multiple MPI chains with mpiexec. The default metropolis_hastings sampler uses one process to coordinate the others, so create one more process than the number of chains you'd like. For example, to get 4 chains,
In [16]:
%%sh
mpiexec -n 5 python -m cosmoslik simple_gaussian.py
The chain is now saved to a file, which can be loaded with a utility function,
In [17]:
simple_gaussian_mpi_chains = utils.load_chain("simple_gaussian.chain")
In [18]:
simple_gaussian_mpi_chains.plot()
Note that we purposely gave a poor choice of scale in the script file. The metropolis_hastings sampler automatically does proposal updating when running with MPI, as seen above.
Now that we understand the general idea of running chains, let's show how to run a Planck chain. We'll give the full script then explain the pieces one-by-one.
In [33]:
from cosmoslik import *
In [34]:
import os.path as osp
param = param_shortcut('start','scale')
class planck(SlikPlugin):
def __init__(self, camspec_clik_file, model='lcdm'):
super(SlikPlugin,self).__init__(**all_kw(locals()))
self.cosmo = get_plugin('models.cosmology')(
logA = param(3.2),
ns = param(0.96),
ombh2 = param(0.0221),
omch2 = param(0.12),
tau = param(0.09,min=0,gaussian_prior=(0.085,0.015)),
theta = param(0.010413),
omnuh2 = 0.000645,
massive_neutrinos=1,
massless_neutrinos=2.046,
)
if 'neff' in model: self.cosmo.massless_neutrinos = param(3,.2)
if 'yp' in model: self.cosmo.Yp = param(.24,0.1)
if 'mnu' in model: self.cosmo.omnuh2 = param(0,0.001,range=(0,1))
self.camspec = get_plugin('likelihoods.clik_like')(
clik_file=camspec_clik_file,
A_ps_100=param(150,min=0),
A_ps_143=param(60,min=0),
A_ps_217=param(60,min=0),
A_cib_143=param(10,min=0),
A_cib_217=param(40,min=0),
A_sz=param(5,scale=1,range=(0,20)),
r_ps=param(0.7,range=(0,1)),
r_cib=param(0.7,range=(0,1)),
n_Dl_cib=param(0.8,scale=0.2,gaussian_prior=(0.8,0.2)),
cal_100=param(1,scale=0.001),
cal_217=param(1,scale=0.001),
xi_sz_cib=param(0.5,range=(-1,1),scale=0.2),
A_ksz=param(1,range=(0,5)),
Bm_1_1=param(0,gaussian_prior=(0,1),scale=1)
)
self.get_cmb = get_plugin('models.camb')()
self.bbn = get_plugin('models.bbn_consistency')()
self.hubble_theta = get_plugin('models.hubble_theta')()
self.priors = get_plugin('likelihoods.priors')(self)
self.sampler = get_plugin('samplers.metropolis_hastings')(
self,
num_samples=1000000,
output_file='chains/chain_%s_%s.chain'%(model,osp.basename(camspec_clik_file)),
proposal_cov='planck_lcdm.covmat',
proposal_scale=1,
print_level=1,
output_extra_params=['cosmo.Yp','cosmo.H0']
)
def __call__(self):
self.cosmo.As = exp(self.cosmo.logA)*1e-10
if 'yp' not in self.model: self.cosmo.Yp = self.bbn(**self.cosmo)
self.cosmo.H0 = self.hubble_theta.theta_to_hubble(**self.cosmo)
self.cmb_result = self.get_cmb(outputs=['cl_TT'],**self.cosmo)
return lsum(lambda: self.priors(self),
lambda: self.camspec(self.cmb_result))
In [35]:
planck_slik = Slik(planck('/home/marius/workspace/planck/clik/CAMspec_v6.2TN_2013_02_26_dist.clik','lcdm'))
In [36]:
lnl, planck_eval = planck_slik.evaluate(**{k:v.start for k,v in planck_slik.get_sampled().items()})
In [37]:
planck_eval.camspec(planck_eval.cmb_result)
Out[37]:
In [38]:
semilogy(planck_eval.cmb_result['cl_TT'])
Out[38]:
In [12]:
clik.clik
Out[12]:
In [38]: