Notebook

Parallel MULTINEST with 3ML¶

J. Michael Burgess¶

MULTINEST¶

MULTINEST is a Bayesian posterior sampler that has two distinct advantages over traditional MCMC:

Recovering multimodal posteriors
- In the case that the posterior is does not have a single maximum, traditional MCMC may miss other modes of the posterior
Full marginal likelihood computation
- This allows for direct model comparison via Bayes factors

To run the MULTINEST sampler in 3ML, one must have the foloowing software installed:

MULTINEST (http://xxx.lanl.gov/abs/0809.3437) (git it here: https://github.com/JohannesBuchner/MultiNest)
pymultinest (https://github.com/JohannesBuchner/PyMultiNest)

Parallelization¶

MULTINEST can be run in a single instance, but it can be incredibly slow. Luckily, it can be built with MPI support enabling it to be run on a multicore workstation or cluster very effeciently.

There are multiple ways to invoke the parallel run of MULTINEST in 3ML: e.g., one can write a python script with all operations and invoke:

$> mpiexec -n <N> python my3MLscript.py

However, it is nice to be able to stay in the Jupyter environment with ipyparallel which allow the user to easily switch bewteen single instance, desktop cores, and cluster environment all with the same code.

Setup¶

The user is expected to have and MPI distribution installed (open-mpi, mpich) and have compiled MULTINEST against the MPI library. Additionally, the user should setup and ipyparallel profile. Instructions can be found here: http://ipython.readthedocs.io/en/2.x/parallel/parallel_mpi.html

Initialize the MPI engine¶

Details for luanching ipcluster on a distributed cluster are not covered here, but everything is the same otherwise.

In the directory that you want to run 3ML in the Jupyter notebook launch and ipcontroller:

$> ipcontroller start --profile=mpi --ip='*'

Next, launch MPI with the desired number of engines:

$> mpiexec -n <N> ipengine --mpi=mpi4py --profile=mpi

Now, the user can head to the notebook and begin!

Running 3ML¶

First we get a client and and connect it to the running profile

In [1]:

from ipyparallel import Client
rc = Client(profile='mpi')
# Grab a view
view = rc[:]

# Activate parallel cell magics
view.activate()

Import 3ML and astromodels to the workers

In [2]:

with view.sync_imports():
    import threeML
    import astromodels

WARNING CppInterfaceNotAvailable: The cthreeML package is not installed. You will not be able to use plugins which require the C/C++ interface (currently HAWC)

Configuration read from /Users/jburgess/.threeML/threeML_config.yml

WARNING NaimaNotAvailable: The naima package is not available. Models that depend on it will not be available

importing threeML on engine(s)
importing astromodels on engine(s)

WARNING CannotImportPlugin: Could not import plugin /usr/local/lib/python2.7/site-packages/threeML-0.2.0-py2.7.egg/threeML/plugins/FermiLATLike.py. Do you have the relative instrument software installed and configured?


WARNING CannotImportPlugin: Could not import plugin /usr/local/lib/python2.7/site-packages/threeML-0.2.0-py2.7.egg/threeML/plugins/HAWCLike.py. Do you have the relative instrument software installed and configured?


WARNING CannotImportPlugin: Could not import plugin /usr/local/lib/python2.7/site-packages/threeML-0.2.0-py2.7.egg/threeML/plugins/SherpaLike.py. Do you have the relative instrument software installed and configured?

Now we set up the analysis in the normal way except the following two caveats:

we must call the threeML module explicity because ipyparallel does not support from <> import *
we use the %%px cell magic (or %px line magic) to perfrom operations in the workers

In [3]:

%%px

# Make GBM detector objects
src_selection = "0.-10."

nai0 = threeML.FermiGBM_TTE_Like('NAI0',
                         "glg_tte_n0_bn080916009_v01.fit",
                         "-10-0,100-200", # background selection
                         src_selection,          # source interval
                         rspfile="glg_cspec_n0_bn080916009_v07.rsp")

nai3 = threeML.FermiGBM_TTE_Like('NAI3',"glg_tte_n3_bn080916009_v01.fit",
                         "-10-0,100-200",
                         src_selection,
                         rspfile="glg_cspec_n3_bn080916009_v07.rsp")

nai4 = threeML.FermiGBM_TTE_Like('NAI4',"glg_tte_n4_bn080916009_v01.fit",
                         "-10-0,100-200",
                         src_selection,
                         rspfile="glg_cspec_n4_bn080916009_v07.rsp")


bgo0 = threeML.FermiGBM_TTE_Like('BGO0',"glg_tte_b0_bn080916009_v01.fit",
                         "-10-0,100-200",
                         src_selection,
                         rspfile="glg_cspec_b0_bn080916009_v07.rsp")

# Select measurements
nai0.set_active_measurements("10.0-30.0", "40.0-950.0")
nai3.set_active_measurements("10.0-30.0", "40.0-950.0")
nai4.set_active_measurements("10.0-30.0", "40.0-950.0")
bgo0.set_active_measurements("250-43000")


# Set up 3ML likelihood object
triggerName = 'bn080916009'
ra = 121.8
dec = -61.3


data_list = threeML.DataList( nai0,nai3,nai4,bgo0 )

band = astromodels.Band()


GRB = threeML.PointSource( triggerName, ra, dec, spectral_shape=band )

model = threeML.Model( GRB )


# Set up Bayesian details

bayes = threeML.BayesianAnalysis(model, data_list)

band.K.prior     = astromodels.Log_uniform_prior(lower_bound=1E-2, upper_bound=5)
band.xp.prior    = astromodels.Log_uniform_prior(lower_bound=1E2, upper_bound=2E3)
band.alpha.prior = astromodels.Uniform_prior(lower_bound=-1.5,upper_bound=0.)
band.beta.prior  = astromodels.Uniform_prior(lower_bound=-3.,upper_bound=-1.5)

Out[3]:

<AsyncResult: execute>

Finally we call MULTINEST. If all is set up properly, MULTINEST will gather the distributed objects and quickly sample the posterior

In [4]:

%px samples = bayes.sample_multinest(n_live_points=400,resume=False)

Out[4]:

<AsyncResult: execute>

Viewing the results¶

Now we need to bring the BayesianAnalysis object back home. Unfortunately, not all objects can be brought back. So you must save figures to the workers. Future implementations of 3ML will allow for saving of the results to a dedicated format which can then be viewed on the local machine. More soon!

In [17]:

# Execute commands that allow for saving figures
# grabbing samples, etc

%%px --targets ::1

samples = bayes.raw_samples()
f=bayes.get_credible_intervals()
bayes.corner_plot(plot_contours=True, plot_density=False)

Out[17]:

<AsyncResult: execute>

In [22]:

# Bring the raw samples local
raw_samples=view['samples'][0]

In [23]:

raw_samples['bn080916009.spectrum.main.Band.K']

Out[23]:

array([ 0.02768163,  0.0326805 ,  0.03299265, ...,  0.03026221,
        0.03026221,  0.03110709])

In [ ]: