Charge Carrier Trapping Experiment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Charge Carrier Trapping Experiment\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, DataFrame, concat\n",
    "from time   import time\n",
    "\n",
    "from exoparams import PlanetParams\n",
    "from astropy import units as u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test_data = read_csv('test_data.dat')\n",
    "test_data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Check if the data makes sense as is**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = figure(figsize=(10,10))\n",
    "errorbar(test_data['DeltaPhase'], test_data['Flux'] , test_data['Sigma'], fmt='o')\n",
    "\n",
    "fig = figure(figsize=(10,10))\n",
    "errorbar(np.arange(test_data['DeltaPhase'].size), test_data['Flux'] , test_data['Sigma'], fmt='o')\n",
    "# errorbar(h38_v2_orbitPhased['DeltaPhase'], h38_v2_orbitPhased['Flux'] , h38_v2_orbitPhased['Sigma'], fmt='o')\n",
    "# errorbar(w67_v1_orbitPhased['DeltaPhase'], w67_v1_orbitPhased['Flux'] , w67_v1_orbitPhased['Sigma'], fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for k in test_data['OrbitNumber'].unique():\n",
    "    orbitNNow = test_data['OrbitNumber'] == k\n",
    "    errorbar(test_data['DeltaPhase'][orbitNNow]    , \\\n",
    "             test_data['Flux'][orbitNNow]          , \\\n",
    "             test_data['Sigma'][orbitNNow], fmt='o')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def zhou_model(params):\n",
    "    # Zhou et al. 2017\n",
    "    # The exponential ramp timescale is detector dependennt, and therfore uniform across all observations\n",
    "    # But the difference from orbit to orbit is predicted to be related \n",
    "    #   ONLY to the inital number of charge traps populated at the start of the each ramp\n",
    "    \n",
    "    # BIG ASSUMPTION\n",
    "    \n",
    "    flux    = ydata.copy() # I assume that what Zhou means by `flux` is either the WLC or avg WLC value\n",
    "    # flux    = ydata.copy() / 128 # I assume that what Zhou means by `flux` is either the WLC or avg WLC value\n",
    "    \n",
    "    E0fast  = params[0] # Orbit 0; Start with per frame; may want per pixel\n",
    "    E0slow  = params[1] # Orbit 0; Start with per frame; may want per pixel\n",
    "    \n",
    "    # Separate out the delta-E0 components per orbit\n",
    "    # Keep dE0fast[0] and dE0slow[0] == 0.0 because they correspond to E0fast and E0slow (initial)\n",
    "    dE0fast     = np.zeros(nOrbits)\n",
    "    dE0slow     = np.zeros(nOrbits)\n",
    "    for k in range(1, nOrbits):\n",
    "        print(k,2*k, 2*k+1,len(params))\n",
    "        dE0fast[k] = params[2*k]\n",
    "        dE0slow[k] = params[2*k+1]\n",
    "    \n",
    "    # From Table 3 of Zhou et al. 2017\n",
    "    ETotFast  = 270.6\n",
    "    etaFast   = 0.006863\n",
    "    tauFast   = 224.8\n",
    "    \n",
    "    ETotSlow  = 1320.0\n",
    "    etaSlow   = 0.01311\n",
    "    tauSlow   = 2.45e4\n",
    "    \n",
    "    coeffFast0= (etaFast * flux / ETotFast + tauFast**-1)\n",
    "    coeffSlow0= (etaSlow * flux / ETotSlow + tauSlow**-1)\n",
    "    \n",
    "    coeffFast1= etaFast*flux / coeffFast0\n",
    "    coeffSlow1= etaSlow*flux / coeffSlow0\n",
    "    \n",
    "    Efast     = zeros(orbit_phase.shape)\n",
    "    Eslow     = zeros(orbit_phase.shape)\n",
    "    for k in range(nOrbits):\n",
    "        orbitNow         = where(orbitNumber == k)[0]\n",
    "        Efast[orbitNow]  = coeffFast1 + (E0fast + dE0fast[k] - coeffFast1)*exp(-coeffFast0 * tphase[orbitNow])\n",
    "        Eslow[orbitNow]  = coeffSlow1 + (E0slow + dE0slow[k] - coeffSlow1)*exp(-coeffSlow0 * tphase[orbitNow])\n",
    "    \n",
    "    dEFastDtP = etaFast * flux * (ETotFast - Efast) / ETotFast\n",
    "    dEFastDtN = -Efast / tauFast\n",
    "    \n",
    "    dESlowDtP = etaSlow * flux * (ETotSlow - Eslow) / ETotSlow\n",
    "    dESlowDtN = -Eslow / tauSlow\n",
    "    \n",
    "    lambda phase: 1 - dEFastDtP - dESlowDtP - dEFastDtP - dESlowDtP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PyMultiNest Demo"
   ]
  },
  {
   "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",
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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 single_exponential_model(cube):\n",
    "    alpha = cube[0]\n",
    "    beta  = cube[1]\n",
    "    gamma = cube[2]\n",
    "    \n",
    "    return lambda xdata: alpha - beta*exp(-gamma*xdata)\n",
    "\n",
    "def double_exponential_model(cube):\n",
    "    alpha   = cube[0]\n",
    "    beta    = cube[1]\n",
    "    gamma   = cube[2]\n",
    "    delta   = cube[3]\n",
    "    epsilon = cube[4]\n",
    "    \n",
    "    return lambda xdata: alpha - beta*exp(-gamma*xdata) - delta*exp(-epsilon*xdata)\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Generate Fake Data for Algorithm Testing **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "\n",
    "param0_test= 1#0.05\n",
    "param1_test= .1#5*pi\n",
    "param2_test= 10.0\n",
    "yunc_test  = 0.01\n",
    "nPts_test  = int(50)\n",
    "nThPts_test= int(1e3)\n",
    "\n",
    "xmin_test  = -0.0#*pi\n",
    "xmax_test  =  1.0#*pi\n",
    "dx_test    = 0.01*(xmax_test - xmin_test)\n",
    "\n",
    "model_test = single_exponential_model\n",
    "# model_test = sine_wave\n",
    "# model_test = straight_line\n",
    "\n",
    "yuncs_test = np.random.normal(yunc_test, 1e-2 * yunc_test, nPts_test)\n",
    "thdata_test= np.linspace(xmin_test-dx_test, xmax_test+dx_test, nThPts_test)\n",
    "\n",
    "xdata_test = np.random.uniform(xmin_test, xmax_test, nPts_test)\n",
    "xdata_test = sort(xdata_test)\n",
    "\n",
    "ydata_test = model_test([param0_test,param1_test,param2_test])(xdata_test)\n",
    "\n",
    "yerr_test  = np.random.normal(0, yuncs_test, nPts_test)\n",
    "zdata_test = ydata_test + yerr_test\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata_test, model_test([param0_test,param1_test,param2_test])(thdata_test))\n",
    "errorbar(xdata_test, zdata_test, yunc_test*ones(zdata_test.size), fmt='o')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Single Exponential Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nThPts    = int(1e3)\n",
    "model_SEM = single_exponential_model\n",
    "\n",
    "xdata = test_data['DeltaPhase']\n",
    "ydata = test_data['Flux']\n",
    "yuncs = test_data['Sigma']\n",
    "\n",
    "xmin, xmax = xdata.min(), xdata.max()\n",
    "dx         = (xmax - xmin)/100\n",
    "thdata_SEM = np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "param0_SEM_init= 1.0  # by defintion\n",
    "param1_SEM_init= (ydata.max() - ydata.min())#/100\n",
    "param2_SEM_init= round(5/(xdata.max() - xdata.min()))\n",
    "\n",
    "print(param0_SEM_init, param1_SEM_init, param2_SEM_init)\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata_SEM, model_SEM([param0_SEM_init,param1_SEM_init,param2_SEM_init])(thdata_SEM))\n",
    "errorbar(xdata, ydata, yuncs, 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\n",
    "# parameters = [\"amp\", \"period\"]\n",
    "\n",
    "# model = straight_line\n",
    "# parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "model_SEM       = single_exponential_model\n",
    "parameters_SEM  = ['max', 'amp1', 'scale1']\n",
    "\n",
    "def myprior_SEM(cube, ndim, nparams):\n",
    "    cube0_width = 1e-3\n",
    "    cube[0] = cube[0] * cube0_width  + (1 - 0.5*cube0_width)# - 10# U(0,2)\n",
    "    cube[1] = cube[1] # - 10# U(0,1) -- default\n",
    "    cube[2] = cube[2] * 1e4   - 5e3# - 1000 # U(0,2000)\n",
    "\n",
    "def myloglike_SEM(cube, ndim, nparams):\n",
    "    chi = 1.\n",
    "    # print \"cube\", [cube[i] for i in range(ndim)], cube\n",
    "    # for i in range(ndim):\n",
    "    #     chi *= -0.5 * ((cube[i] - 0.2) / 0.1)**2#math.cos(cube[i] / 2.) * math.sin(cube[i] / 2.)\n",
    "    # print \"returning\", math.pow(2. + chi, 5)\n",
    "    modelNow = model_SEM(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": 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_SEM = len(parameters_SEM)\n",
    "\n",
    "planetName           = 'HAT38'\n",
    "visitName            = 'visit1'\n",
    "modelName            = 'single_exponential_model'\n",
    "outputfiles_basename = 'chains/' + planetName + '-' + visitName + '-' + modelName + '-'\n",
    "\n",
    "start = time()\n",
    "\n",
    "plt.figure(figsize=(5*n_params_SEM, 5*n_params_SEM))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params_SEM, outputfiles_basename=outputfiles_basename); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike_SEM, myprior_SEM, n_params_SEM, \\\n",
    "                importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "                sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename=outputfiles_basename);\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop();\n",
    "print('SEM took', time() - start, 'seconds')\n",
    "# lets analyse the results\n",
    "a_SEM = pymultinest.Analyzer(n_params = n_params_SEM, outputfiles_basename=outputfiles_basename);\n",
    "s_SEM = a_SEM.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_SEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_SEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_SEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_SEM, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\t%.15e +- %.15e\" % ( s_SEM['nested sampling global log-evidence'], \\\n",
    "                                             s_SEM['nested sampling global log-evidence error'] ))\n",
    "print(\"Global Evidence:\\t%.3f +- %.3f\" % ( s_SEM['nested sampling global log-evidence'], \\\n",
    "                                             s_SEM['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_SEM = pymultinest.PlotMarginalModes(a_SEM)\n",
    "plt.figure(figsize=(5*n_params_SEM, 5*n_params_SEM))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params_SEM):\n",
    "    plt.subplot(n_params_SEM, n_params_SEM, n_params_SEM * i + i + 1)\n",
    "    p_SEM.plot_marginal(i, with_ellipses = False, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_SEM[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params_SEM, n_params_SEM, n_params_SEM * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p_SEM.plot_conditional(i, j, with_ellipses = False, with_points = False, grid_points=30)\n",
    "        plt.xlabel(parameters_SEM[i])\n",
    "        plt.ylabel(parameters_SEM[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 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_SEM):\n",
    "    # print(5*n_params, 1, i+1)\n",
    "    plt.subplot(5*n_params_SEM, 1, i+1)\n",
    "    p_SEM.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_SEM[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_SEM.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters_SEM[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_SEM.analyser.get_best_fit()['parameters'], [param0_SEM_init, param1_SEM_init, param2_SEM_init]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "plot(thdata_SEM, model_SEM([param0_SEM_init,param1_SEM_init, param2_SEM_init])(thdata_SEM), label='Initial Model')\n",
    "errorbar(xdata, ydata, yuncs, fmt='o', label='Data')\n",
    "plot(thdata_SEM, model_SEM(p_SEM.analyser.get_best_fit()['parameters'])(thdata_SEM), label='PMN Model')\n",
    "legend(loc=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_SEM.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unrestricted Double Exponential Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nThPts= int(1e3)\n",
    "model_UDEM = double_exponential_model\n",
    "\n",
    "xdata = test_data['DeltaPhase']\n",
    "ydata = test_data['Flux']\n",
    "yuncs = test_data['Sigma']\n",
    "\n",
    "xmin, xmax = xdata.min(), xdata.max()\n",
    "dx         = (xmax - xmin)/100\n",
    "thdata_UDEM = np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "param0_UDEM_init  = 1.0  # by defintion\n",
    "param1_UDEM_init  = 0.5*(ydata.max() - ydata.min())#/100\n",
    "param2_UDEM_init  = round(5/(xdata.max() - xdata.min()))\n",
    "param3_UDEM_init  = 0.5*(ydata.max() - ydata.min())#/100\n",
    "param4_UDEM_init  = round(5/(xdata.max() - xdata.min()))\n",
    "\n",
    "print(param0_UDEM_init, param1_UDEM_init, param2_UDEM_init, param3_UDEM_init, param4_UDEM_init)\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata_UDEM, model_UDEM([param0_UDEM_init,param1_UDEM_init,param2_UDEM_init, param3_UDEM_init, param4_UDEM_init])(thdata_UDEM))\n",
    "errorbar(xdata, ydata, yuncs, 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\n",
    "# parameters = [\"amp\", \"period\"]\n",
    "\n",
    "# model = straight_line\n",
    "# parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "model_UDEM      = double_exponential_model\n",
    "parameters_UDEM = ['max', 'amp1', 'scale1', 'amp2', 'scale2']\n",
    "\n",
    "# def myprior_RDEM(cube, ndim, nparams):\n",
    "#     cube[0] = cube[0] * 1e-3  + (1 - 1e-3/2)# - 10# U(0,2)\n",
    "#     cube[1] = -cube[1] * 5e-3 + 5e-4 # - 10# U(0,1) -- default\n",
    "#     cube[2] = cube[2] * 1e4   - 5e3# - 1000 # U(0,2000)\n",
    "#     cube[3] = cube[3] * 5e-3  + 5e-4# - 10# U(0,1) -- default\n",
    "#     cube[4] = cube[4] * 1e4   - 5e3# - 1000 # U(0,2000)\n",
    "\n",
    "def myprior_UDEM(cube, ndim, nparams):\n",
    "    cube[0] = cube[0] * 1e-2  + (1 - 1e-2/2)# - 10# U(0,2)\n",
    "    cube[1] = cube[1] * 2     - 2/2 # - 10# U(0,1) -- default\n",
    "    cube[2] = cube[2] * 1e4   - 1e4/2# - 1000 # U(0,2000)\n",
    "    cube[3] = cube[3] * 2     - 2/2# - 10# U(0,1) -- default\n",
    "    cube[4] = cube[4] * 1e4   - 1e4/2# - 1000 # U(0,2000)\n",
    "\n",
    "def myloglike_UDEM(cube, ndim, nparams):\n",
    "    modelNow = model_UDEM(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "start = time()\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params_UDEM = len(parameters_UDEM)\n",
    "\n",
    "savedir              = 'chains'\n",
    "planetName           = 'HAT38'\n",
    "visitName            = 'visit1'\n",
    "modelName            = 'unrestricted_double_exponential_model'\n",
    "outputfiles_basename = savedir + '/' + planetName + '-' + visitName + '-' + modelName + '-'\n",
    "\n",
    "plt.figure(figsize=(5*n_params_UDEM, 5*n_params_UDEM))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params_UDEM, outputfiles_basename=outputfiles_basename)\n",
    "progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike_UDEM, myprior_UDEM, n_params_UDEM, importance_nested_sampling = False, resume = False, \\\n",
    "                verbose = True, sampling_efficiency = 'model', n_live_points = 1000, \\\n",
    "                outputfiles_basename=outputfiles_basename)\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a_UDEM = pymultinest.Analyzer(n_params = n_params_UDEM, outputfiles_basename=outputfiles_basename)\n",
    "s_UDEM = a_UDEM.get_stats()\n",
    "print('UDEM took', time() - start, 'seconds')\n",
    "\n",
    "# fig = plt.gcf()\n",
    "# axs = fig.get_axes()\n",
    "\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim()\n",
    "#     # ax.set_xscale(\"log\", nonposx='clip')\n",
    "#     # ax.set_yscale(\"log\", nonposy='clip')"
   ]
  },
  {
   "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_UDEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_UDEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_UDEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_UDEM, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\t%.15e +- %.15e\" % ( s_UDEM['nested sampling global log-evidence'], \\\n",
    "                                               s_UDEM['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_UDEM = pymultinest.PlotMarginalModes(a_UDEM)\n",
    "plt.figure(figsize=(5*n_params_UDEM, 5*n_params_UDEM))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params_UDEM):\n",
    "    plt.subplot(n_params_UDEM, n_params_UDEM, n_params_UDEM * i + i + 1)\n",
    "    p_UDEM.plot_marginal(i, with_ellipses = False, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_UDEM[i])\n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params_UDEM, n_params_UDEM, n_params_UDEM * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        # p_UDEM.plot_conditional(i, j, with_ellipses=False, with_points=False, grid_points=30)\n",
    "        p_UDEM.plot_conditional(i, j, with_ellipses=False , with_points=False , grid_points=30, \\\n",
    "                                only_interpolate=False, use_log_values=False, marginalization_type='sum')\n",
    "        plt.xlabel(parameters_UDEM[i])\n",
    "        plt.ylabel(parameters_UDEM[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "axes_colors = rcParams['axes.prop_cycle'].by_key()['color']\n",
    "nColors     = len(axes_colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "minLogE, maxLogE = min(a_UDEM.get_equal_weighted_posterior().T[-1]), max(a_UDEM.get_equal_weighted_posterior().T[-1])\n",
    "rangeLogE        = maxLogE - minLogE\n",
    "minLogE, maxLogE, rangeLogE, nColors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from astroML.plotting import hist\n",
    "from statsmodels.robust import scale\n",
    "hist(a_UDEM.get_equal_weighted_posterior().T[-1], bins='blocks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nSig      = 10\n",
    "mad_logE  = scale.mad(a_UDEM.get_equal_weighted_posterior().T[-1])\n",
    "med_logE  = median(a_UDEM.get_equal_weighted_posterior().T[-1])\n",
    "\n",
    "madBins   = [med_logE - nSig*mad_logE for nSig in range(nColors)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "fig = figure(figsize=(15,15));\n",
    "\n",
    "logEchain = a_UDEM.get_equal_weighted_posterior().T[-1]\n",
    "mad_logE  = scale.mad(a_UDEM.get_equal_weighted_posterior().T[-1])\n",
    "med_logE  = median(a_UDEM.get_equal_weighted_posterior().T[-1])\n",
    "\n",
    "madBins   = [med_logE - nSig*mad_logE for nSig in range(nColors+1)]\n",
    "\n",
    "for k in range(5):\n",
    "    ax = fig.add_subplot(5,1,k+1);\n",
    "    for nSig in range(nColors):\n",
    "        for klogE in range(logEchain.size):\n",
    "            if logEchain[klogE] > madBins[nSig] or logEchain[klogE] < madBins[nSig+1]:\n",
    "                ax.plot(a_UDEM.get_equal_weighted_posterior().T[k],\\\n",
    "                        logEchain,'o', \\\n",
    "                        color = axes_colors[nSig], alpha=0.1);\n",
    "\n",
    "fig.canvas.draw();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 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_UDEM):\n",
    "    # print(5*n_params, 1, i+1)\n",
    "    plt.subplot(5*n_params_UDEM, 1, i+1)\n",
    "    p_UDEM.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_UDEM[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_UDEM.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters_UDEM[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_UDEM.analyser.get_best_fit()['parameters'], [param0_UDEM_init, param1_UDEM_init, param2_UDEM_init, param3_UDEM_init, param4_UDEM_init]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "plot(thdata_UDEM, model_UDEM([param0_UDEM_init,param1_UDEM_init, param2_UDEM_init, param3_UDEM_init, param4_UDEM_init])(thdata_UDEM), label='Initial Model')\n",
    "errorbar(xdata, ydata, yuncs, fmt='o', label='Data')\n",
    "plot(thdata_UDEM, model_UDEM(p_UDEM.analyser.get_best_fit()['parameters'])(thdata_UDEM), label='PMN UDEM Model')\n",
    "legend(loc=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_UDEM.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Restricted Double Exponential Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nThPts= int(1e3)\n",
    "model_RDEM = double_exponential_model\n",
    "\n",
    "xdata = test_data['DeltaPhase']\n",
    "ydata = test_data['Flux']\n",
    "yuncs = test_data['Sigma']\n",
    "\n",
    "xmin, xmax = xdata.min(), xdata.max()\n",
    "dx         = (xmax - xmin)/100\n",
    "thdata_RDEM = np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "param0_RDEM_init  = 1.0  # by defintion\n",
    "param1_RDEM_init  = 0.5*(ydata.max() - ydata.min())#/100\n",
    "param2_RDEM_init  = round(5/(xdata.max() - xdata.min()))\n",
    "param3_RDEM_init  = 0.5*(ydata.max() - ydata.min())#/100\n",
    "param4_RDEM_init  = round(5/(xdata.max() - xdata.min()))\n",
    "\n",
    "print(param0_RDEM_init, param1_RDEM_init, param2_RDEM_init, param3_RDEM_init, param4_RDEM_init)\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata_RDEM, model_RDEM([param0_RDEM_init,param1_RDEM_init,param2_RDEM_init, param3_RDEM_init, param4_RDEM_init])(thdata_RDEM))\n",
    "errorbar(xdata, ydata, yuncs, fmt='o')"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "def show(filepath): \n",
    "    \"\"\" open the output (pdf) file for the user \"\"\"\n",
    "    if os.name == 'mac' or platform == 'darwin': subprocess.call(('open', filepath))\n",
    "    elif os.name == 'nt' or platform == 'win32': os.startfile(filepath)\n",
    "    elif platform.startswith('linux') : subprocess.call(('xdg-open', filepath))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# our probability functions\n",
    "# Taken from the eggbox problem.\n",
    "# model = sine_wave\n",
    "# parameters = [\"amp\", \"period\"]\n",
    "\n",
    "# model = straight_line\n",
    "# parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "model_RDEM      = double_exponential_model\n",
    "parameters_RDEM = ['max', 'amp1', 'scale1', 'amp2', 'scale2']\n",
    "\n",
    "def myprior_RDEM(cube, ndim, nparams):\n",
    "    cube[0] = cube[0] * 1e-3  + (1 - 1e-3/2)# - 10# U(0,2)\n",
    "    cube[1] = -cube[1] * 5e-3 + 5e-4 # - 10# U(0,1) -- default\n",
    "    cube[2] = cube[2] * 1e4   - 5e3# - 1000 # U(0,2000)\n",
    "    cube[3] = cube[3] * 5e-3  + 5e-4# - 10# U(0,1) -- default\n",
    "    cube[4] = cube[4] * 1e4   - 5e3# - 1000 # U(0,2000)\n",
    "\n",
    "def myloglike_RDEM(cube, ndim, nparams):\n",
    "    chi = 1.\n",
    "    # print \"cube\", [cube[i] for i in range(ndim)], cube\n",
    "    # for i in range(ndim):\n",
    "    #     chi *= -0.5 * ((cube[i] - 0.2) / 0.1)**2#math.cos(cube[i] / 2.) * math.sin(cube[i] / 2.)\n",
    "    # print \"returning\", math.pow(2. + chi, 5)\n",
    "    modelNow = model_RDEM(cube)(xdata)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "start = time()\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params_RDEM = len(parameters_RDEM)\n",
    "\n",
    "savedir              = 'chains'\n",
    "planetName           = 'HAT38'\n",
    "visitName            = 'visit1'\n",
    "modelName            = 'restricted_double_exponential_model'\n",
    "outputfiles_basename = savedir + '/' + planetName + '-' + visitName + '-' + modelName + '-'\n",
    "\n",
    "plt.figure(figsize=(5*n_params_RDEM, 5*n_params_RDEM))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params_RDEM, outputfiles_basename=outputfiles_basename); progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike_RDEM, myprior_RDEM, n_params_RDEM, importance_nested_sampling = False, resume = False, verbose = True, \\\n",
    "            sampling_efficiency = 'model', n_live_points = 1000, outputfiles_basename=outputfiles_basename)\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a_RDEM = pymultinest.Analyzer(n_params = n_params_RDEM, outputfiles_basename=outputfiles_basename)\n",
    "s_RDEM = a_RDEM.get_stats()\n",
    "print('RDEM took', time() - start, 'seconds')\n",
    "\n",
    "# fig = plt.gcf()\n",
    "# axs = fig.get_axes()\n",
    "\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim()\n",
    "#     # ax.set_xscale(\"log\", nonposx='clip')\n",
    "#     # ax.set_yscale(\"log\", nonposy='clip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a_RDEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_RDEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_RDEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_RDEM, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\t%.15e +- %.15e\" % ( s_RDEM['nested sampling global log-evidence'], \\\n",
    "                                               s_RDEM['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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_RDEM = pymultinest.PlotMarginalModes(a_RDEM)\n",
    "plt.figure(figsize=(5*n_params_RDEM, 5*n_params_RDEM))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params_RDEM):\n",
    "    plt.subplot(n_params_RDEM, n_params_RDEM, n_params_RDEM * i + i + 1)\n",
    "    p_RDEM.plot_marginal(i, with_ellipses = False, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_RDEM[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params_RDEM, n_params_RDEM, n_params_RDEM * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p_RDEM.plot_conditional(i, j, with_ellipses = False, with_points = False, grid_points=30)\n",
    "        plt.xlabel(parameters_RDEM[i])\n",
    "        plt.ylabel(parameters_RDEM[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 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_RDEM):\n",
    "    # print(5*n_params, 1, i+1)\n",
    "    plt.subplot(5*n_params_RDEM, 1, i+1)\n",
    "    p_RDEM.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_RDEM[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_RDEM.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters_RDEM[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_RDEM.analyser.get_best_fit()['parameters'], [param0_RDEM_init, param1_RDEM_init, param2_RDEM_init, param3_RDEM_init, param4_RDEM_init]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "plot(thdata_RDEM, model_RDEM([param0_RDEM_init,param1_RDEM_init, param2_RDEM_init, param3_RDEM_init, param4_RDEM_init])(thdata_RDEM), label='Initial Model')\n",
    "errorbar(xdata, ydata, yuncs, fmt='o', label='Data')\n",
    "plot(thdata_RDEM, model_RDEM(p_RDEM.analyser.get_best_fit()['parameters'])(thdata_RDEM), label='PMN Model')\n",
    "legend(loc=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_RDEM.analyser.get_stats()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare Unrestricted Double, Restricted Double, and Single Exponential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a_SEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_SEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_SEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_SEM, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"SEM Global Evidence:\\t\\t%.3f +- %.3f\" % ( s_SEM['nested sampling global log-evidence'], \\\n",
    "                                                   s_SEM['nested sampling global log-evidence error'] ))\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a_UDEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_UDEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_UDEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_UDEM, f, indent=2)\n",
    "\n",
    "# print()\n",
    "# print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"UDEM Global Evidence:\\t\\t%.3f +- %.3f\" % ( s_UDEM['nested sampling global log-evidence'], \\\n",
    "                                                   s_UDEM['nested sampling global log-evidence error'] ))\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a_RDEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_RDEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_RDEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_RDEM, f, indent=2)\n",
    "\n",
    "# print()\n",
    "# print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"RDEM Global Evidence:\\t\\t%.3f +- %.3f\" % ( s_RDEM['nested sampling global log-evidence'], \\\n",
    "                                                   s_RDEM['nested sampling global log-evidence error'] ))\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata_UDEM, model_SEM([param0_SEM_init,param1_SEM_init, param2_SEM_init])(thdata_SEM), \\\n",
    "     '.', label='Initial SEM Model')\n",
    "plot(thdata_UDEM, \\\n",
    " model_RDEM([param0_RDEM_init,param1_RDEM_init, param2_RDEM_init, param3_RDEM_init, param4_RDEM_init])(thdata_RDEM), \\\n",
    " '--', label='Initial DEM Model')\n",
    "errorbar(xdata, ydata, yuncs, fmt='o', label='Data')\n",
    "plot(thdata_SEM, model_SEM(p_SEM.analyser.get_best_fit()['parameters'])(thdata_SEM), label='PMN SEM Model')\n",
    "plot(thdata_UDEM, model_UDEM(p_UDEM.analyser.get_best_fit()['parameters'])(thdata_UDEM), label='PMN UDEM Model')\n",
    "plot(thdata_RDEM, model_RDEM(p_RDEM.analyser.get_best_fit()['parameters'])(thdata_RDEM), label='PMN RDEM Model')\n",
    "legend(loc=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "figure(figsize=(20,20))\n",
    "from numpy.polynomial import polynomial\n",
    "def time_polynomial(params):\n",
    "    # modelOut = np.zeros(tdata.size)\n",
    "    # for kc, coeff in enumerate(params):\n",
    "    #     modelOut += coeff * tdata**kc\n",
    "    if len(params):\n",
    "        return lambda tdata: polynomial.polyval(tdata, params)\n",
    "    else:\n",
    "        return lambda tdata: zeros(tdata.size)\n",
    "\n",
    "def orbital_polynomial(params):\n",
    "    # modelOut = np.zeros(xdata.size)\n",
    "    # for kc, coeff in enumerate(params):\n",
    "    #     modelOut += coeff * xdata**kc\n",
    "    # return modelOut\n",
    "    if len(params):\n",
    "        return lambda odata: polynomial.polyval(odata, params)\n",
    "    else:\n",
    "        return lambda odata: zeros(odata.size)\n",
    "\n",
    "def wavelength_polynomial(params):\n",
    "    # modelOut = np.zeros(ldata.size)\n",
    "    # for kc, coeff in enumerate(params):\n",
    "    #     modelOut += coeff * ldata**kc\n",
    "    # return modelOut\n",
    "    if len(params):\n",
    "        return lambda ldata: polynomial.polyval(ldata, params)\n",
    "    else:\n",
    "        return lambda ldata: zeros(ldata.size)\n",
    "\n",
    "def polynomial_model(params):\n",
    "    params_list= list(params.copy())[::-1]\n",
    "    timeParams  = array([params_list.pop() for _ in range(nTimeCoeffs)])\n",
    "    orbitParams = array([params_list.pop() for _ in range(nOrbitCoeffs)])\n",
    "    waveParams  = array([params_list.pop() for _ in range(nWaveCoeffs)])\n",
    "    return lambda tdata, odata, ldata: time_polynomial(timeParams)(tdata)      + \\\n",
    "                                       orbital_polynomial(orbitParams)(odata)  + \\\n",
    "                                       wavelength_polynomial(waveParams)(ldata)\n",
    "\n",
    "tdata, xdata, ldata = np.random.uniform(-10,10,(3,100))\n",
    "tdata.sort()\n",
    "xdata.sort()\n",
    "ldata.sort()\n",
    "\n",
    "# tdata, xdata, ldata = [np.linspace(-10,10,100) for _ in range(3)]\n",
    "for nTimeCoeffs in range(4):\n",
    "    for nOrbitCoeffs in range(4):\n",
    "        for nWaveCoeffs in range(4):\n",
    "            params = np.random.uniform(-20,20,nTimeCoeffs+nOrbitCoeffs+nWaveCoeffs)\n",
    "            plot(tdata, polynomial_model(params)(tdata,xdata,ldata),'.', alpha=0.5, mew=0)\n",
    "            plot(xdata, polynomial_model(params)(tdata,xdata,ldata),'.', alpha=0.5, mew=0)\n",
    "            plot(ldata, polynomial_model(params)(tdata,xdata,ldata),'.', alpha=0.5, mew=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nThPts= int(1e3)\n",
    "model_Poly  = polynomial_model\n",
    "\n",
    "xdata = test_data['DeltaPhase']\n",
    "ydata = test_data['Flux']\n",
    "yuncs = test_data['Sigma']\n",
    "\n",
    "nTimeCoeffs     = 2\n",
    "nOrbitCoeffs    = 3\n",
    "nWaveCoeffs     = 0\n",
    "\n",
    "h38PlanetPhase  = test_data_input_input['Phase']\n",
    "h38HSTPhase     = test_data['DeltaPhase']\n",
    "\n",
    "xmin, xmax = xdata.min(), xdata.max()\n",
    "dx         = (xmax - xmin)/100\n",
    "thdata_Poly = np.linspace(xmin-dx, xmax+dx, nThPts)\n",
    "\n",
    "param0_Poly_init  = 1.0  # by defintion\n",
    "param1_Poly_init  = 1.0\n",
    "param2_Poly_init  = 1.0\n",
    "param3_Poly_init  = 1.0\n",
    "param4_Poly_init  = 1.0\n",
    "\n",
    "print(param0_Poly_init, param1_Poly_init, param2_Poly_init, param3_Poly_init, param4_Poly_init)\n",
    "\n",
    "figure(figsize=(10,10))\n",
    "plot(thdata_Poly, model_Poly([param0_Poly_init,param1_Poly_init,param2_Poly_init, param3_Poly_init, param4_Poly_init])(thdata_Poly))\n",
    "errorbar(xdata, ydata, yuncs, 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\n",
    "# parameters = [\"amp\", \"period\"]\n",
    "\n",
    "# model = straight_line\n",
    "# parameters = [\"offset\", \"slope\"]\n",
    "\n",
    "nTimeCoeffs     = 2\n",
    "nOrbitCoeffs    = 3\n",
    "nWaveCoeffs     = 0\n",
    "\n",
    "h38PlanetPhase  = test_data_input_input['Phase']\n",
    "h38HSTPhase     = test_data['DeltaPhase']\n",
    "\n",
    "model_Poly      = polynomial_model\n",
    "parameters_Poly = ['timeIntercept', 'timeSlope', 'orbitIntercept', 'orbitSlope', 'orbitQuadratic']\n",
    "\n",
    "cubeKWith = 1e3\n",
    "def myprior_Poly(cube, ndim, nparams):\n",
    "    for k in len(cube):\n",
    "        cube[k] = cube[k] * cubeKWith - 0.5*cubeKWith\n",
    "\n",
    "def myloglike_Poly(cube, ndim, nparams):\n",
    "    chi = 1.\n",
    "    # print \"cube\", [cube[i] for i in range(ndim)], cube\n",
    "    # for i in range(ndim):\n",
    "    #     chi *= -0.5 * ((cube[i] - 0.2) / 0.1)**2#math.cos(cube[i] / 2.) * math.sin(cube[i] / 2.)\n",
    "    # print \"returning\", math.pow(2. + chi, 5)\n",
    "    modelNow = model_Poly(cube)(times, HSTPhase, 0)\n",
    "    return -0.5*((modelNow - ydata)**2. / yuncs**2.).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if not os.path.exists(\"chains\"): os.mkdir(\"chains\")\n",
    "\n",
    "start = time()\n",
    "# number of dimensions our problem has\n",
    "# parameters = [\"x\", \"y\"]\n",
    "n_params_Poly = len(parameters_Poly)\n",
    "\n",
    "savedir              = 'chains'\n",
    "planetName           = 'HAT38'\n",
    "visitName            = 'visit1'\n",
    "modelName            = 'polynomial_model'\n",
    "outputfiles_basename = savedir + '/' + planetName + '-' + visitName + '-' + modelName + '-'\n",
    "\n",
    "plt.figure(figsize=(5*n_params_Poly, 5*n_params_Poly))\n",
    "# we want to see some output while it is running\n",
    "progress = pymultinest.ProgressPlotter(n_params = n_params_Poly, outputfiles_basename=outputfiles_basename)\n",
    "\n",
    "progress.start()\n",
    "# threading.Timer(2, show, [\"chains/2-phys_live.points.pdf\"]).start() # delayed opening\n",
    "# run MultiNest\n",
    "pymultinest.run(myloglike_Poly, myprior_Poly, n_params_Poly, importance_nested_sampling = False, \\\n",
    "                resume = False, verbose = True, sampling_efficiency = 'model', n_live_points = 1000, \\\n",
    "                outputfiles_basename=outputfiles_basename)\n",
    "\n",
    "# ok, done. Stop our progress watcher\n",
    "progress.stop()\n",
    "\n",
    "# lets analyse the results\n",
    "a_Poly = pymultinest.Analyzer(n_params = n_params_Poly, outputfiles_basename=outputfiles_basename)\n",
    "s_Poly = a_Poly.get_stats()\n",
    "print('Polynomial took', time() - start, 'seconds')\n",
    "\n",
    "# fig = plt.gcf()\n",
    "# axs = fig.get_axes()\n",
    "\n",
    "# for ax in axs:\n",
    "#     ax.set_ylim()\n",
    "#     # ax.set_xscale(\"log\", nonposx='clip')\n",
    "#     # ax.set_yscale(\"log\", nonposy='clip')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "# store name of parameters, always useful\n",
    "with open('%sparams.json' % a_RDEM.outputfiles_basename, 'w') as f:\n",
    "    json.dump(parameters_RDEM, f, indent=2)\n",
    "# store derived stats\n",
    "with open('%sstats.json' % a_RDEM.outputfiles_basename, mode='w') as f:\n",
    "    json.dump(s_RDEM, f, indent=2)\n",
    "\n",
    "print()\n",
    "print(\"-\" * 30, 'ANALYSIS', \"-\" * 30)\n",
    "print(\"Global Evidence:\\t%.15e +- %.15e\" % ( s_RDEM['nested sampling global log-evidence'], \\\n",
    "                                               s_RDEM['nested sampling global log-evidence error'] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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_RDEM = pymultinest.PlotMarginalModes(a_RDEM)\n",
    "plt.figure(figsize=(5*n_params_RDEM, 5*n_params_RDEM))\n",
    "#plt.subplots_adjust(wspace=0, hspace=0)\n",
    "for i in range(n_params_RDEM):\n",
    "    plt.subplot(n_params_RDEM, n_params_RDEM, n_params_RDEM * i + i + 1)\n",
    "    p_RDEM.plot_marginal(i, with_ellipses = False, with_points = False, grid_points=50)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_RDEM[i])\n",
    "    \n",
    "    for j in range(i):\n",
    "        plt.subplot(n_params_RDEM, n_params_RDEM, n_params_RDEM * j + i + 1)\n",
    "        #plt.subplots_adjust(left=0, bottom=0, right=0, top=0, wspace=0, hspace=0)\n",
    "        p_RDEM.plot_conditional(i, j, with_ellipses = False, with_points = False, grid_points=30)\n",
    "        plt.xlabel(parameters_RDEM[i])\n",
    "        plt.ylabel(parameters_RDEM[j])\n",
    "\n",
    "# plt.savefig(\"chains/marginals_multinest.pdf\") #, bbox_inches='tight')\n",
    "# show(\"chains/marginals_multinest.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# 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_RDEM):\n",
    "    # print(5*n_params, 1, i+1)\n",
    "    plt.subplot(5*n_params_RDEM, 1, i+1)\n",
    "    p_RDEM.plot_modes_marginal(i, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Probability\")\n",
    "    plt.xlabel(parameters_RDEM[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_RDEM.plot_modes_marginal(i, cumulative = True, with_ellipses = True, with_points = False)\n",
    "    plt.ylabel(\"Cumulative probability\")\n",
    "    plt.xlabel(parameters_RDEM[i])\n",
    "    # plt.savefig(outfile, format='pdf', bbox_inches='tight')\n",
    "    # plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_RDEM.analyser.get_best_fit()['parameters'], [param0_RDEM_init, param1_RDEM_init, param2_RDEM_init, param3_RDEM_init, param4_RDEM_init]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "figure(figsize=(10,10))\n",
    "plot(thdata_RDEM, model_RDEM([param0_RDEM_init,param1_RDEM_init, param2_RDEM_init, param3_RDEM_init, param4_RDEM_init])(thdata_RDEM), label='Initial Model')\n",
    "errorbar(xdata, ydata, yuncs, fmt='o', label='Data')\n",
    "plot(thdata_RDEM, model_RDEM(p_RDEM.analyser.get_best_fit()['parameters'])(thdata_RDEM), label='PMN Model')\n",
    "legend(loc=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p_RDEM.analyser.get_stats()"
   ]
  }
 ],
 "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
}