.ipynb_checkpoints/PyMultiNest Learning Tutorial-checkpoint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyMultiNest My Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import absolute_import, unicode_literals, print_function\n",
    "import pymultinest\n",
    "import math\n",
    "import os\n",
    "import threading, subprocess\n",
    "from sys import platform\n",
    "\n",
    "from pylab import *;ion()\n",
    "\n",
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Straight Line **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian1D(cube):\n",
    "    # cube = np.float32(cube)\n",
    "    return lambda y: np.exp(-0.5*(( (cube[0] - y) / cube[1])**2)) / sqrt(2*pi*cube[1]**2)\n",
    "\n",
    "def straight_line(cube):\n",
    "    offset = cube[0]\n",
    "    slope  = cube[1]\n",
    "    return lambda abscissa: offset + slope * abscissa\n",
    "\n",
    "def sine_wave(cube):\n",
    "    amp    = cube[0]\n",
    "    period = cube[1]\n",
    "    return lambda abscissa: amp*sin(2*pi / period * abscissa)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0= 0.5#0.05\n",
    "param1= 0.9#5*pi\n",
    "yunc  = 0.1\n",
    "nPts  = int(50)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -1#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, nPts)\n",
    "thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.random.uniform(xmin, xmax, nPts)\n",
    "xdata = sort(xdata)\n",
    "\n",
    "ydata = model([param0,param1])(xdata)\n",
    "\n",
    "yerr  = np.random.normal(0, yuncs, nPts)\n",
    "zdata = ydata + yerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata, model([param0,param1])(thdata))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# If you want to \"mess up\" the model (for model comparison), \n",
    "#   then uncomment here a *different* model than above\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    # assumes a uniform(0,1) prior in nparams-space\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    # Normal Log Likelihood\n",
    "    modelNow = model(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "            sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3))\n",
    "\n",
    "print('true\\t', [param0, param1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# figure(figsize=(10,10))\n",
    "# plot(thdata, model(p.analyser.get_best_fit()['parameters'])(thdata))\n",
    "# errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')\n",
    "# plot(thdata, model([param0,param1])(thdata))\n",
    "figure(figsize=(10,10))\n",
    "# plot(thdata, model([param0a,param1a])(thdata) + model([param0b,param1b])(thdata))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')\n",
    "modelAll = np.zeros(thdata.size)\n",
    "for m in p.analyser.get_stats()['modes']:\n",
    "    modelAll = modelAll + model(m['mean'])(thdata)\n",
    "    plot(thdata, model(m['mean'])(thdata))\n",
    "\n",
    "plot(thdata, modelAll)\n",
    "plot(thdata, model(p.analyser.get_best_fit()['parameters'])(thdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Sine Wave **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian1D(cube):\n",
    "    # cube = np.float32(cube)\n",
    "    return lambda y: np.exp(-0.5*(( (cube[0] - y) / cube[1])**2)) / sqrt(2*pi*cube[1]**2)\n",
    "\n",
    "def straight_line(cube):\n",
    "    offset = cube[0]\n",
    "    slope  = cube[1]\n",
    "    return lambda abscissa: offset + slope * abscissa\n",
    "\n",
    "def sine_wave(cube):\n",
    "    amp    = cube[0]\n",
    "    period = cube[1]\n",
    "    return lambda abscissa: amp*sin(2*pi / period * abscissa)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0= 0.5#0.05\n",
    "param1= 0.9#5*pi\n",
    "yunc  = 0.1\n",
    "nPts  = int(50)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -1#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, nPts)\n",
    "thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.random.uniform(xmin, xmax, nPts)\n",
    "xdata = sort(xdata)\n",
    "\n",
    "ydata = model([param0,param1])(xdata)\n",
    "\n",
    "yerr  = np.random.normal(0, yuncs, nPts)\n",
    "zdata = ydata + yerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata, model([param0,param1])(thdata))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# If you want to \"mess up\" the model (for model comparison), \n",
    "#   then uncomment here a *different* model than above\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    # assumes a uniform(0,1) prior in nparams-space\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    # Normal Log Likelihood\n",
    "    modelNow = model(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "            sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3))\n",
    "\n",
    "print('true\\t', [param0, param1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')\n",
    "modelAll = np.zeros(thdata.size)\n",
    "\n",
    "sumLogE = -\n",
    "for m in p.analyser.get_stats()['modes']:\n",
    "    sumLogE += modelAll + model(m['mean'])(thdata) \n",
    "for m in p.analyser.get_stats()['modes']:\n",
    "    modelAll = modelAll + model(m['mean'])(thdata)\n",
    "    plot(thdata, model(m['mean'])(thdata))\n",
    "\n",
    "plot(thdata, modelAll)\n",
    "plot(thdata, model(p.analyser.get_best_fit()['parameters'])(thdata))\n",
    "\n",
    "figure()\n",
    "modelAll = np.zeros(thdata.size)\n",
    "errorbar(xdata[1:], np.diff(zdata),yunc*ones(zdata.size-1), fmt='o')\n",
    "for km, mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    if km:\n",
    "        plot(thdata, model(mode['mean'])(thdata))\n",
    "        modelAll = modelAll + model(m['mean'])(thdata)\n",
    "plot(thdata, modelAll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Gaussian Normal**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian1Dp(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    height = cube[2]\n",
    "    return lambda y: height*np.exp(-0.5*(( (center - y) / width)**2))# / sqrt(2*pi*width**2)\n",
    "\n",
    "def gaussian1D(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    return lambda y: np.exp(-0.5*(( (center - y) / width)**2)) / sqrt(2*pi*width**2)\n",
    "\n",
    "def straight_line(cube):\n",
    "    offset = cube[0]\n",
    "    slope  = cube[1]\n",
    "    return lambda abscissa: offset + slope * abscissa\n",
    "\n",
    "def sine_wave(cube):\n",
    "    amp    = cube[0]\n",
    "    period = cube[1]\n",
    "    return lambda abscissa: amp*sin(2*pi / period * abscissa)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0= 0.5#0.05\n",
    "param1= 0.1#0.05\n",
    "# param2= 0.8\n",
    "\n",
    "yunc  = 0.1\n",
    "nPts  = int(100)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -0#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = gaussian1Dp; parameters = [\"center\", \"width\", \"height\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, nPts)\n",
    "thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.random.uniform(xmin, xmax, nPts)\n",
    "xdata = sort(xdata)\n",
    "\n",
    "ydata = model([param0,param1])(xdata)\n",
    "\n",
    "yerr  = np.random.normal(0, yuncs, nPts)\n",
    "zdata = ydata + yerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata, model([param0,param1])(thdata))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# Taken from the eggbox problem.\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    modelNow = model(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "            sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()\n",
    "\n",
    "# fig = gcf()\n",
    "# axs  = fig.get_axes()\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim(-16,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3))\n",
    "\n",
    "print('True\\t', [param0, param1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')\n",
    "modelAll = np.zeros(thdata.size)\n",
    "for m in p.analyser.get_stats()['modes']:\n",
    "    modelAll = modelAll + model(m['mean'])(thdata) #* model(m[''])\n",
    "    plot(thdata, model(m['mean'])(thdata))\n",
    "\n",
    "plot(thdata, modelAll)\n",
    "plot(thdata, model(p.analyser.get_best_fit()['parameters'])(thdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** MultiModel Gaussian Normal**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian1Dp(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    height = cube[2]\n",
    "    return lambda y: height*np.exp(-0.5*(( (center - y) / width)**2))# / sqrt(2*pi*width**2)\n",
    "\n",
    "def gaussian1D(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    return lambda y: np.exp(-0.5*(( (center - y) / width)**2)) / sqrt(2*pi*width**2)\n",
    "\n",
    "def straight_line(cube):\n",
    "    offset = cube[0]\n",
    "    slope  = cube[1]\n",
    "    return lambda abscissa: offset + slope * abscissa\n",
    "\n",
    "def sine_wave(cube):\n",
    "    amp    = cube[0]\n",
    "    period = cube[1]\n",
    "    return lambda abscissa: amp*sin(2*pi / period * abscissa)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0a= -0.5#0.05\n",
    "param0b= 0.5#0.05\n",
    "param1a= 0.1#5*pi\n",
    "param1b= 0.1#5*pi\n",
    "# param2= 0.8\n",
    "\n",
    "yunc  = 0.1\n",
    "nPts  = int(100)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -1#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = gaussian1Dp; parameters = [\"center\", \"width\", \"height\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, nPts)\n",
    "thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.linspace(xmin,xmax,nPts)\n",
    "# xdata = np.random.uniform(xmin, xmax, nPts)\n",
    "# xdata = sort(xdata)\n",
    "\n",
    "ydata = model([param0a,param1a])(xdata) + model([param0b,param1b])(xdata)\n",
    "\n",
    "yerr  = np.random.normal(0, yuncs, nPts)\n",
    "zdata = ydata + yerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata, model([param0a,param1a])(thdata) + model([param0b,param1b])(thdata))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# Taken from the eggbox problem.\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    cube[0] = cube[0]*2 - 1\n",
    "    cube[1] = cube[1]*2\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    modelNow = model(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "            sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()\n",
    "\n",
    "# fig = gcf()\n",
    "# axs  = fig.get_axes()\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim(-16,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3))\n",
    "\n",
    "print('True a\\t', [param0a, param1a])\n",
    "print('True b\\t', [param0b, param1b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')\n",
    "modelAll = np.zeros(thdata.size)\n",
    "for m in p.analyser.get_stats()['modes']:\n",
    "    modelAll = modelAll + model(m['mean'])(thdata)# * model(m[''])\n",
    "    plot(thdata, model(m['mean'])(thdata))\n",
    "\n",
    "plot(thdata, modelAll)\n",
    "plot(thdata, model(p.analyser.get_best_fit()['parameters'])(thdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** MultiModel Gaussian Normal with Height Free**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian1Dp(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    height = cube[2]\n",
    "    return lambda y: height*np.exp(-0.5*(( (center - y) / width)**2))# / sqrt(2*pi*width**2)\n",
    "\n",
    "def gaussian1D(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    return lambda y: np.exp(-0.5*(( (center - y) / width)**2)) / sqrt(2*pi*width**2)\n",
    "\n",
    "def straight_line(cube):\n",
    "    offset = cube[0]\n",
    "    slope  = cube[1]\n",
    "    return lambda abscissa: offset + slope * abscissa\n",
    "\n",
    "def sine_wave(cube):\n",
    "    amp    = cube[0]\n",
    "    period = cube[1]\n",
    "    return lambda abscissa: amp*sin(2*pi / period * abscissa)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0a= -0.5#0.05\n",
    "param0b= 0.5#0.05\n",
    "param1a= 0.1#5*pi\n",
    "param1b= 0.1#5*pi\n",
    "param2a= 0.8\n",
    "param2b= 0.8\n",
    "\n",
    "yunc  = 0.1\n",
    "nPts  = int(100)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -1#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "model = gaussian1Dp; parameters = [\"center\", \"width\", \"height\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, nPts)\n",
    "thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.linspace(xmin,xmax,nPts)\n",
    "# xdata = np.random.uniform(xmin, xmax, nPts)\n",
    "# xdata = sort(xdata)\n",
    "\n",
    "ydata = model([param0a,param1a,param2a])(xdata) + model([param0b,param1b,param2b])(xdata)\n",
    "\n",
    "yerr  = np.random.normal(0, yuncs, nPts)\n",
    "zdata = ydata + yerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata, model([param0a,param1a,param2a])(thdata) + model([param0b,param1b,param2b])(thdata))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# Taken from the eggbox problem.\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    cube[0] = cube[0]*2 - 1\n",
    "    cube[1] = cube[1]*2\n",
    "    cube[2] = cube[2]*2\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    modelNow = model(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "            sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()\n",
    "\n",
    "# fig = gcf()\n",
    "# axs  = fig.get_axes()\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim(-16,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3))\n",
    "\n",
    "print('True a\\t', [param0a, param1a])\n",
    "print('True b\\t', [param0b, param1b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "errorbar(xdata, zdata, yunc*ones(zdata.size), fmt='o')\n",
    "modelAll = np.zeros(thdata.size)\n",
    "for m in p.analyser.get_stats()['modes']:\n",
    "    modelAll = modelAll + model(m['mean'])(thdata)# * model(m[''])\n",
    "    plot(thdata, model(m['mean'])(thdata))\n",
    "\n",
    "plot(thdata, modelAll)\n",
    "plot(thdata, model(p.analyser.get_best_fit()['parameters'])(thdata))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 2D Gaussian Modeling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2D Gaussian Normal Mono-Modal**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nPts  = int(100)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -0#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "ymin  = -0#*pi\n",
    "ymax  =  1#*pi\n",
    "dy    = 0.1*(ymax - ymin)\n",
    "\n",
    "xdata = np.ones((nPts,nPts))*np.linspace(xmin,xmax,nPts)\n",
    "ydata = (np.ones((nPts,nPts))*np.linspace(ymin,ymax,nPts)).T\n",
    "\n",
    "fig = figure(figsize=(16,8))\n",
    "ax = fig.add_subplot(221)\n",
    "ax.imshow(xdata)\n",
    "ax.set_ylim(0,nPts)\n",
    "ax = fig.add_subplot(222)\n",
    "ax.imshow(ydata)\n",
    "ax.set_ylim(0,nPts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian2D(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    return lambda y,x: np.exp(-0.5*((( (center - y) / width)**2) + (( (center - x) / width)**2))) / sqrt(2*pi*width**2)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0= 0.5#0.05\n",
    "param1= 0.1#0.05\n",
    "# param2= 0.8\n",
    "\n",
    "yunc  = 0.1\n",
    "nPts  = int(100)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -0#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "ymin  = -0#*pi\n",
    "ymax  =  1#*pi\n",
    "dy    = 0.1*(ymax - ymin)\n",
    "\n",
    "model = gaussian2D; parameters = [\"center\", \"width\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, (nPts,nPts))\n",
    "# thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.ones((nPts,nPts))*np.linspace(xmin,xmax,nPts)\n",
    "ydata = (np.ones((nPts,nPts))*np.linspace(ymin,ymax,nPts)).T\n",
    "\n",
    "zmodel  = model([param0,param1])(ydata,xdata)\n",
    "zerr    = np.random.normal(0, yuncs, (nPts,nPts))\n",
    "zdata   = zmodel + zerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "# imshow(thdata, model([param0,param1])(thdata))\n",
    "imshow(zdata, extent=[xdata.min(), xdata.max(), ydata.min(), ydata.max()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# Taken from the eggbox problem.\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    modelNow = gaussian2D(cube)(ydata,xdata)\n",
    "    return -0.5*((modelNow - zdata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, \\\n",
    "                resume = False, verbose = True, sampling_efficiency = 'model', n_live_points = 1000, \\\n",
    "                outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()\n",
    "\n",
    "# fig = gcf()\n",
    "# axs  = fig.get_axes()\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim(-16,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3),'\\t', np.round(mode['local log-evidence'],3))\n",
    "\n",
    "print('True\\t', [param0, param1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "modelAll = np.zeros((nPts, nPts))\n",
    "fig=figure(figsize=(20,10))\n",
    "for km,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    modelAll = modelAll + model(mode['mean'])(ydata,xdata) #* model(m[''])\n",
    "    ax = fig.add_subplot(1,len(p.analyser.get_stats()['modes']), km+1)\n",
    "    ims = ax.imshow(model(mode['mean'])(ydata,xdata))\n",
    "    plt.colorbar(ims)\n",
    "\n",
    "fig = figure(figsize=(20,10))\n",
    "ax = fig.add_subplot(131)\n",
    "ims = ax.imshow(modelAll)\n",
    "plt.colorbar(ims)\n",
    "ax = fig.add_subplot(132)\n",
    "ims = ax.imshow(model(p.analyser.get_best_fit()['parameters'])(ydata,xdata))\n",
    "plt.colorbar(ims)\n",
    "ax = fig.add_subplot(133)\n",
    "# ims = ax.imshow(model([param0, param1])(ydata,xdata))\n",
    "ims = ax.imshow(zdata)\n",
    "plt.colorbar(ims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2D Gaussian Normal Multi-Modal**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def gaussian2D(cube):\n",
    "    center = cube[0]\n",
    "    width  = cube[1]\n",
    "    return lambda y,x: np.exp(-0.5*((( (center - y) / width)**2) + (( (center - x) / width)**2))) / sqrt(2*pi*width**2)\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "param0a= 0.75#0.05\n",
    "param1a= 0.05#0.05\n",
    "param0b= 0.25#0.05\n",
    "param1b= 0.05#0.05\n",
    "\n",
    "# param2= 0.8\n",
    "\n",
    "yunc  = 0.1\n",
    "nPts  = int(100)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -0#*pi\n",
    "xmax  =  1#*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "ymin  = -0#*pi\n",
    "ymax  =  1#*pi\n",
    "dy    = 0.1*(ymax - ymin)\n",
    "\n",
    "model = gaussian2D; parameters = [\"center\", \"width\"]\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, (nPts,nPts))\n",
    "# thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.ones((nPts,nPts))*np.linspace(xmin,xmax,nPts)\n",
    "ydata = (np.ones((nPts,nPts))*np.linspace(ymin,ymax,nPts)).T\n",
    "\n",
    "zmodel  = model([param0a,param1a])(ydata,xdata) + model([param0b,param1b])(ydata,xdata)\n",
    "zerr    = np.random.normal(0, yuncs, (nPts,nPts))\n",
    "zdata   = zmodel + zerr\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "imshow(zdata, extent=[xdata.min(), xdata.max(), ydata.min(), ydata.max()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# Taken from the eggbox problem.\n",
    "# model = sine_wave; parameters = [\"amp\", \"period\"]\n",
    "# model = gaussian1D; parameters = [\"center\", \"width\"]\n",
    "# model = straight_line; parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "def myprior(cube, ndim, nparams):\n",
    "    pass\n",
    "\n",
    "def myloglike(cube, ndim, nparams):\n",
    "    modelNow = gaussian2D(cube)(ydata,xdata)\n",
    "    return -0.5*((modelNow - zdata)**2. / yuncs**2.).sum()\n",
    "\n",
    "def null_loglike(cube, ndim, nparams):\n",
    "    modelNow = gaussian2D(cube)(ydata,xdata)\n",
    "    return -0.5*((modelNow)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params, outputfiles_basename='chains/2-'); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike, myprior, n_params, importance_nested_sampling = False, \\\n",
    "                resume = False, verbose = True, sampling_efficiency = 'model', n_live_points = 1000, \\\n",
    "                outputfiles_basename='chains/2-')\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a = pymultinest.Analyzer(n_params = n_params, outputfiles_basename='chains/2-')\n",
    "s = a.get_stats()\n",
    "\n",
    "# fig = gcf()\n",
    "# axs  = fig.get_axes()\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim(-16,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\n\\t%.15e +- %.15e\" % ( s['nested sampling global log-evidence'], s['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.clf()\n",
    "\n",
    "# Here we will plot all the marginals and whatnot, just to show off\n",
    "# You may configure the format of the output here, or in matplotlibrc\n",
    "# All pymultinest does is filling in the data of the plot.\n",
    "\n",
    "# Copy and edit this file, and play with it.\n",
    "\n",
    "p = pymultinest.PlotMarginalModes(a)\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params):\n",
    "    plt.subplot(n_params, n_params, n_params * i + i + 1)\n",
    "    p.plot_marginal(i, with_ellipses = True, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params, n_params, n_params * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p.plot_conditional(i, j, with_ellipses = False, with_points = True, grid_points=30)\n",
    "        plt.xlabel(parameters[i])\n",
    "        plt.ylabel(parameters[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")\n",
    "\n",
    "plt.figure(figsize=(5*n_params, 5*n_params))\n",
    "plt.subplot2grid((5*n_params, 5*n_params), loc=(0,0))\n",
    "for i in range(n_params):\n",
    "    #plt.subplot(n_params, n_params, i + 1)\n",
    "    # outfile = '%s-mode-marginal-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "    \n",
    "    # outfile = '%s-mode-marginal-cumulative-%d.pdf' % (a.outputfiles_basename,i)\n",
    "    p.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()\n",
    "\n",
    "print(\"Take a look at the pdf files in chains/\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "null_logE = -4.43430028425182973E+04\n",
    "\n",
    "print('best\\t', np.round(p.analyser.get_best_fit()['parameters'],3))\n",
    "for k,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    print('mode' + str(k) + '\\t', np.round(mode['mean'],3),'\\t', np.round(mode['local log-evidence'],3))\n",
    "\n",
    "print('True a\\t', [param0a, param1a])\n",
    "print('True b\\t', [param0b, param1b])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "modelAll = np.zeros((nPts, nPts))\n",
    "fig=figure(figsize=(20,10))\n",
    "for km,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    modelAll = modelAll + model(mode['mean'])(ydata,xdata)\n",
    "    ax = fig.add_subplot(1,len(p.analyser.get_stats()['modes']), km+1)\n",
    "    ims = ax.imshow(model(mode['mean'])(ydata,xdata))\n",
    "    plt.colorbar(ims)\n",
    "\n",
    "fig = figure(figsize=(20,10))\n",
    "ax = fig.add_subplot(131)\n",
    "ims = ax.imshow(modelAll)\n",
    "plt.colorbar(ims)\n",
    "ax = fig.add_subplot(132)\n",
    "ims = ax.imshow(model(p.analyser.get_best_fit()['parameters'])(ydata,xdata))\n",
    "plt.colorbar(ims)\n",
    "ax = fig.add_subplot(133)\n",
    "ims = ax.imshow(zdata)\n",
    "plt.colorbar(ims)\n",
    "\n",
    "# Residuals\n",
    "modelAll = np.zeros((nPts, nPts))\n",
    "for km,mode in enumerate(p.analyser.get_stats()['modes']):\n",
    "    if np.round(mode['local log-evidence'],3) > -400000.:\n",
    "        modelAll = modelAll + model(mode['mean'])(ydata,xdata)\n",
    "\n",
    "fig = figure(figsize=(20,10))\n",
    "ax  = fig.add_subplot(131)\n",
    "ims = ax.imshow(zdata - modelAll)\n",
    "plt.colorbar(ims)\n",
    "ax  = fig.add_subplot(132)\n",
    "ax.hist((zdata - modelAll).ravel(), bins=1000, normed=True);\n",
    "# plt.colorbar(ims)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p.analyser.get_stats()\n",
    "p.analyser.post_file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyMultiNest Learning Tutorial\n",
    "\n",
    "CCT = Charge Carrier Trapping - This is a test of comparing the Zhou et al 2017 results with a data driven analysis using multinest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "from pandas import read_csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "if not os.path.exists('chains/'):\n",
    "    os.mkdir('chains')\n",
    "\n",
    "print(os.path.exists('chains'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyMultiNest Solver Class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Initial Example**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "from pymultinest.solve import Solver\n",
    "from numpy import pi, sin\n",
    "\n",
    "class EggboxProblem(Solver):\n",
    "    def Prior(self, cube):\n",
    "        return cube * 10 * pi\n",
    "    def LogLikelihood(self, cube):\n",
    "        chi = (sin(cube)).prod()\n",
    "        return (2. + chi)**5\n",
    "\n",
    "solution = EggboxProblem(n_dims = 1)\n",
    "print(solution)\n",
    "solution = EggboxProblem(n_dims = 2)\n",
    "print(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**My PyMultiNest Test**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from pylab import *;ion()\n",
    "\n",
    "from pymultinest.solve import Solver,solve\n",
    "from numpy import pi, sin, cos, linspace\n",
    "\n",
    "def straight_line(cube):\n",
    "    offset = cube[0]\n",
    "    slope  = cube[1]\n",
    "    return lambda abscissa: offset + slope * abscissa\n",
    "\n",
    "def sine_wave(cube):\n",
    "    amp    = cube[0]\n",
    "    period = cube[1]\n",
    "    return lambda abscissa: amp*sin(2*pi / period * abscissa)\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "param0= 0.1#0.05\n",
    "param1= 0.1#5*pi\n",
    "yunc  = 0.025\n",
    "nPts  = int(10)\n",
    "nThPts= int(1e3)\n",
    "\n",
    "xmin  = -0.5*pi\n",
    "xmax  =  0.5*pi\n",
    "dx    = 0.1*(xmax - xmin)\n",
    "\n",
    "# model = sine_wave\n",
    "model = straight_line\n",
    "\n",
    "yuncs = np.random.normal(yunc, 1e-2 * yunc, nPts)\n",
    "thdata= np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "xdata = np.random.uniform(xmin, xmax, nPts)\n",
    "xdata = sort(xdata)\n",
    "\n",
    "ydata = model([param0,param1])(xdata)\n",
    "\n",
    "yerr  = np.random.normal(0, yuncs, nPts)\n",
    "zdata = ydata + yerr\n",
    "\n",
    "figure(figsize=(20,5))\n",
    "plot(thdata, model([param0, param1])(thdata))\n",
    "errorbar(xdata, zdata, yuncs, fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class ChisqFit(Solver):\n",
    "    def Prior(self, cube):\n",
    "        return cube\n",
    "    \n",
    "    def LogLikelihood(self, cube):\n",
    "        chisq = (-0.5*((model(cube)(xdata) - ydata)**2.) )#/ yuncs**2.\n",
    "        return chisq.prod()\n",
    "\n",
    "# solution = ChisqFit(n_dims = 2, resume=False, outputfiles_basename='./42-')\n",
    "\n",
    "solution = ChisqFit(n_dims = 2, n_params=None, n_clustering_params=None, wrapped_params=None        , \\\n",
    "                    importance_nested_sampling=True, multimodal=True, const_efficiency_mode=False   , \\\n",
    "                    n_live_points=400, evidence_tolerance=0.5, sampling_efficiency=0.8              , \\\n",
    "                    n_iter_before_update=100, null_log_evidence=-1e+90, max_modes=100               , \\\n",
    "                    mode_tolerance=-1e+90, outputfiles_basename='chains/1-', seed=-1, verbose=False , \\\n",
    "                    resume=False, context=0, write_output=True, log_zero=-1e+100, max_iter=0         , \\\n",
    "                    init_MPI=False, dump_callback=None)\n",
    "print(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Simplest Example**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "import pymultinest\n",
    "\n",
    "def prior(cube, ndim, nparams):\n",
    "    cube[0] = cube[0] * 2\n",
    "\n",
    "def loglikelihood(cube, ndim, nparams):\n",
    "    return ((cube[0] - 0.2) / 0.1)**2\n",
    "\n",
    "pymultinest.run(loglikelihood, prior, n_dims=1, max_iter=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyMultiNest Solve Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!/usr/bin/env python\n",
    "from __future__ import absolute_import, unicode_literals, print_function\n",
    "import numpy\n",
    "from numpy import pi, cos\n",
    "from pymultinest.solve import solve\n",
    "import os\n",
    "\n",
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "# probability function, taken from the eggbox problem.\n",
    "\n",
    "def myprior(cube):\n",
    "    return cube * 10 * pi\n",
    "\n",
    "def myloglike(cube):\n",
    "    chi = (cos(cube / 2.)).prod()\n",
    "    return (2. + chi)**5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# number of dimensions our problem has\n",
    "parameters = [\"x\", \"y\"]\n",
    "n_params = len(parameters)\n",
    "\n",
    "# run MultiNest\n",
    "result = solve(LogLikelihood=myloglike, Prior=myprior, \n",
    "    n_dims=n_params, outputfiles_basename=\"chains/3-\")\n",
    "\n",
    "print()\n",
    "print('evidence: %(logZ).1f +- %(logZerr).1f' % result)\n",
    "print()\n",
    "print('parameter values:')\n",
    "for name, col in zip(parameters, result['samples'].transpose()):\n",
    "    print('%15s : %.3f +- %.3f' % (name, col.mean(), col.std()))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}