{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallel MULTINEST with 3ML\n",
    "###### J. Michael Burgess\n",
    "\n",
    "## MULTINEST \n",
    "MULTINEST is a Bayesian posterior sampler that has two distinct advantages over traditional MCMC:\n",
    "* Recovering multimodal posteriors\n",
    "    * In the case that the posterior is does not have a single maximum, traditional MCMC\n",
    "        may miss other modes of the posterior\n",
    "* Full marginal likelihood computation\n",
    "    * This allows for direct model comparison via Bayes factors\n",
    "\n",
    "\n",
    "To run the MULTINEST sampler in **3ML**, one must have the foloowing software installed:\n",
    "* MULTINEST (http://xxx.lanl.gov/abs/0809.3437) (*git* it here: https://github.com/JohannesBuchner/MultiNest)\n",
    "* pymultinest (https://github.com/JohannesBuchner/PyMultiNest)\n",
    "\n",
    "## Parallelization\n",
    "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. \n",
    "\n",
    "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:\n",
    "\n",
    "```bash\n",
    "$> mpiexec -n <N> python my3MLscript.py\n",
    "\n",
    "```\n",
    "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.\n",
    "\n",
    "### Setup\n",
    "\n",
    "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\n",
    "\n",
    "### Initialize the MPI engine\n",
    "Details for luanching ipcluster on a distributed cluster are not covered here, but everything is the same otherwise.\n",
    "\n",
    "In the directory that you want to run 3ML in the Jupyter notebook launch and ipcontroller:\n",
    "\n",
    "```bash\n",
    "$> ipcontroller start --profile=mpi --ip='*'\n",
    "\n",
    "```\n",
    "Next, launch MPI with the desired number of engines:\n",
    "\n",
    "```bash\n",
    "$> mpiexec -n <N> ipengine --mpi=mpi4py --profile=mpi\n",
    "\n",
    "```\n",
    "\n",
    "Now, the user can head to the notebook and begin!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running 3ML\n",
    "First we get a client and and connect it to the running profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from ipyparallel import Client\n",
    "rc = Client(profile='mpi')\n",
    "# Grab a view\n",
    "view = rc[:]\n",
    "\n",
    "# Activate parallel cell magics\n",
    "view.activate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import 3ML and astromodels to the workers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING CppInterfaceNotAvailable: The cthreeML package is not installed. You will not be able to use plugins which require the C/C++ interface (currently HAWC)\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Configuration read from /Users/jburgess/.threeML/threeML_config.yml\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING NaimaNotAvailable: The naima package is not available. Models that depend on it will not be available\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "importing threeML on engine(s)\n",
      "importing astromodels on engine(s)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "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?\n",
      "\n",
      "\n",
      "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?\n",
      "\n",
      "\n",
      "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?\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with view.sync_imports():\n",
    "    import threeML\n",
    "    import astromodels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we set up the analysis in the normal way except the following two caveats:\n",
    "* we must call the threeML module explicity because ipyparallel does not support from <> import \\*\n",
    "* we use the %%px cell magic (or %px line magic) to perfrom operations in the workers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AsyncResult: execute>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%px\n",
    "\n",
    "# Make GBM detector objects\n",
    "src_selection = \"0.-10.\"\n",
    "\n",
    "nai0 = threeML.FermiGBM_TTE_Like('NAI0',\n",
    "                         \"glg_tte_n0_bn080916009_v01.fit\",\n",
    "                         \"-10-0,100-200\", # background selection\n",
    "                         src_selection,          # source interval\n",
    "                         rspfile=\"glg_cspec_n0_bn080916009_v07.rsp\")\n",
    "\n",
    "nai3 = threeML.FermiGBM_TTE_Like('NAI3',\"glg_tte_n3_bn080916009_v01.fit\",\n",
    "                         \"-10-0,100-200\",\n",
    "                         src_selection,\n",
    "                         rspfile=\"glg_cspec_n3_bn080916009_v07.rsp\")\n",
    "\n",
    "nai4 = threeML.FermiGBM_TTE_Like('NAI4',\"glg_tte_n4_bn080916009_v01.fit\",\n",
    "                         \"-10-0,100-200\",\n",
    "                         src_selection,\n",
    "                         rspfile=\"glg_cspec_n4_bn080916009_v07.rsp\")\n",
    "\n",
    "\n",
    "bgo0 = threeML.FermiGBM_TTE_Like('BGO0',\"glg_tte_b0_bn080916009_v01.fit\",\n",
    "                         \"-10-0,100-200\",\n",
    "                         src_selection,\n",
    "                         rspfile=\"glg_cspec_b0_bn080916009_v07.rsp\")\n",
    "\n",
    "# Select measurements\n",
    "nai0.set_active_measurements(\"10.0-30.0\", \"40.0-950.0\")\n",
    "nai3.set_active_measurements(\"10.0-30.0\", \"40.0-950.0\")\n",
    "nai4.set_active_measurements(\"10.0-30.0\", \"40.0-950.0\")\n",
    "bgo0.set_active_measurements(\"250-43000\")\n",
    "\n",
    "\n",
    "# Set up 3ML likelihood object\n",
    "triggerName = 'bn080916009'\n",
    "ra = 121.8\n",
    "dec = -61.3\n",
    "\n",
    "\n",
    "data_list = threeML.DataList( nai0,nai3,nai4,bgo0 )\n",
    "\n",
    "band = astromodels.Band()\n",
    "\n",
    "\n",
    "GRB = threeML.PointSource( triggerName, ra, dec, spectral_shape=band )\n",
    "\n",
    "model = threeML.Model( GRB )\n",
    "\n",
    "\n",
    "# Set up Bayesian details\n",
    "\n",
    "bayes = threeML.BayesianAnalysis(model, data_list)\n",
    "\n",
    "band.K.prior     = astromodels.Log_uniform_prior(lower_bound=1E-2, upper_bound=5)\n",
    "band.xp.prior    = astromodels.Log_uniform_prior(lower_bound=1E2, upper_bound=2E3)\n",
    "band.alpha.prior = astromodels.Uniform_prior(lower_bound=-1.5,upper_bound=0.)\n",
    "band.beta.prior  = astromodels.Uniform_prior(lower_bound=-3.,upper_bound=-1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we call MULTINEST. If all is set up properly, MULTINEST will gather the distributed objects and quickly sample the posterior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AsyncResult: execute>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%px samples = bayes.sample_multinest(n_live_points=400,resume=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewing the results\n",
    "\n",
    "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!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AsyncResult: execute>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute commands that allow for saving figures\n",
    "# grabbing samples, etc\n",
    "\n",
    "%%px --targets ::1\n",
    "\n",
    "samples = bayes.raw_samples()\n",
    "f=bayes.get_credible_intervals()\n",
    "bayes.corner_plot(plot_contours=True, plot_density=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Bring the raw samples local\n",
    "raw_samples=view['samples'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.02768163,  0.0326805 ,  0.03299265, ...,  0.03026221,\n",
       "        0.03026221,  0.03110709])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "raw_samples['bn080916009.spectrum.main.Band.K']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}