{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from pycheops import Dataset, StarProperties\n",
    "from uncertainties import ufloat\n",
    "# Temporary(?) fix for problem with coordinate look-up for HD stars in astropy \n",
    "from astropy.coordinates.name_resolve import sesame_database\n",
    "sesame_database.set('simbad');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create and instance of the dataset\n",
    "\n",
    "*Run pycheops.core.setup_config once first to create the configuration file and data cache directory*\n",
    "\n",
    "```python\n",
    ">>> from pycheops.core import setup_config\n",
    ">>> setup_config()\n",
    "```\n",
    "\n",
    "This step can be used to download the data from CHEOPsim ftp site using the run number, e.g., \n",
    "```python\n",
    ">>> dataset = Dataset().from_simulation(7649, target='HD97658')\n",
    "```\n",
    "\n",
    "The cached version is then loaded using the program number, e.g., \n",
    "```python\n",
    ">>> dataset = Dataset('PR900048_TG001501')\n",
    "```\n",
    "\n",
    "The lightcurve for analysis is then selected using the get_lightcurve method. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/pflm/pycheops_data/CH_PR900048_TG001501.zip already downloaded\n",
      "Writing Imagette data to .tgz file...\n",
      "Writing Lightcurve data to .tgz file...\n",
      ".. CH_PR900048_TG001501_TU2020-04-03T04-11-52_SCI_COR_Lightcurve-DEFAULT_V0000.fits - done\n",
      ".. CH_PR900048_TG001501_TU2020-04-03T04-11-52_SCI_COR_Lightcurve-OPTIMAL_V0000.fits - done\n",
      ".. CH_PR900048_TG001501_TU2020-04-03T04-11-52_SCI_COR_Lightcurve-RINF_V0000.fits - done\n",
      ".. CH_PR900048_TG001501_TU2020-04-03T04-11-52_SCI_COR_Lightcurve-RSUP_V0000.fits - done\n",
      "Found archive tgzfile /Users/pflm/pycheops_data/PR900048_TG001501.tgz\n",
      " PI name     : CHEOPSim\n",
      " OBS ID      : 7649\n",
      " Target      : HD97658\n",
      " Coordinates : 11:14:33.16 +25:42:37.4\n",
      "Light curve data loaded from  /Users/pflm/pycheops_data/PR900048_TG001501-OPTIMAL.fits\n",
      "C_cut = 97317724\n",
      "N(C > C_cut) = 0\n",
      "Mean counts = 97182778.5\n",
      "Median counts = 97196778.9\n",
      "RMS counts = 39092.5 [402 ppm]\n",
      "Median standard error = 97196778.9 [104 ppm]\n"
     ]
    }
   ],
   "source": [
    "dataset = Dataset().from_simulation(7649, target='HD97658')\n",
    "#dataset = Dataset('PR900048_TG001501')\n",
    "_ = dataset.get_lightcurve(aperture='OPTIMAL')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate stellar density and limb darkening from parameers in SWEET-Cat \n",
    "\n",
    "Prior on mean stellar density is helpful to constrain the width of the eclipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Identifier : HD97658\n",
      "Coordinates: 11:14:33.16 +25:42:37.4\n",
      "T_eff :  5137 +/- 36 K\n",
      "log g :  4.47 +/- 0.09 [cgs]\n",
      "[M/H] : -0.35 +/- 0.02\n",
      "log rho :  0.06 +/- 0.16 [solar]\n",
      "h_1 : 0.708 +/- 0.010\n",
      "h_2 : 0.439 +/- 0.051\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pHD97658 = StarProperties(dataset.target)\n",
    "print(pHD97658)\n",
    "logrhoprior = pHD97658.logrho\n",
    "# logrhoprior = ufloat(0.0593, 0.1556) # for offline user"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Least-squares fit\n",
    "\n",
    "Light curve for a planet of radius $R_p$ transiting a star of radius $R_s$ with impact parameter $b$ is calculated using the power-2 limb-darkening law, \n",
    "$$I_{\\lambda}(\\mu) = 1 - c(1 - \\mu^{\\alpha}),$$\n",
    "and is described by the following parameters.\n",
    "\n",
    "- T_0:   time of mid-transit\n",
    "- P:     orbital period\n",
    "- D:     depth, $(R_p/R_s)^2 = k^2$\n",
    "- W:     width, $(R_s/a)\\sqrt{(1+k)^2 - b^2}/\\pi$\n",
    "- S:     shape, $[(1-k)^2-b^2]/[(1+k)^2 - b^2]$\n",
    "- c:     flux scaling factor\n",
    "- f_c:   $\\sqrt{e}\\cos(\\omega)$\n",
    "- f_s:   $\\sqrt{e}\\sin(\\omega)$\n",
    "- h_1:   $I_{\\lambda}(0.5) = 1 - c(1-2^{-\\alpha})$\n",
    "- h_2:   $I_{\\lambda}(0.5) - I_{\\lambda}(0) = c\\cdot2^{-\\alpha}$\n",
    "\n",
    "Model options available here are \n",
    " - linear or quadratic decorrelation against position (dfdx, dfdy, df2dx2, d2fdxdy, d2fdx2, d2fdy2)\n",
    " - linear decorrelation against roll angle (dfdsinphi, dfdcosphi)\n",
    " - linear or quadratic time trend (dfdt, d2fdt2)\n",
    " \n",
    "The stellar density can be included as a prior in the least-squares fit using the keyword logrhoprior.\n",
    "\n",
    "Parameters and their bounds/priors are  specified as follows.\n",
    "\n",
    "- as a floating point number for fixed parameters, e.g., \"P=1.234\".\n",
    "- as a tuple for a free parameter with bounds, e.g., \"h_2 = (0.3,0.7)\" or \"c = (0.9, 1.0, 1.1)\". The median value of the tuple is used as the initial value in the fit.\n",
    "- as a ufloat floating-point value with a standard error, e.g., \"P=ufloat(1.2345,0.0001)\", to specify a free parameter with a Gaussian prior.\n",
    "- as an lmfit Parameter, e.g., \"lmfit.Parameter(value=0.01,min=0.001,max=0.1,vary=True)\"\n",
    "\n",
    "The output from this method is an lmfit [MinimizerResult](https://lmfit.github.io/lmfit-py/fitting.html#minimizerresult-the-optimization-result) object.\n",
    "\n",
    "The limb-darkening is not well constrained by this light curve so we fix the parameter h_2 at its nominal value (\".n\") and fit the parameter h_1 with a prior.\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Fit Statistics]]\n",
      "    # fitting method   = leastsq\n",
      "    # function evals   = 126\n",
      "    # data points      = 580\n",
      "    # variables        = 7\n",
      "    chi-square         = 1534.45546\n",
      "    reduced chi-square = 2.67793274\n",
      "    Akaike info crit   = 578.283592\n",
      "    Bayesian info crit = 608.824789\n",
      "[[Variables]]\n",
      "    T_0:     0.89055878 +/- 4.4763e-04 (0.05%) (init = 0.89)\n",
      "    P:       9.48926400 +/- 1.6364e-05 (0.00%) (init = 9.489264)\n",
      "    D:       7.4000e-04 +/- 6.6703e-05 (9.01%) (init = 0.001244335)\n",
      "    W:       0.01426471 +/- 1.4973e-04 (1.05%) (init = 0.02208278)\n",
      "    S:       0.87971847 +/- 0.04220332 (4.80%) (init = 0.8283019)\n",
      "    f_c:     0 (fixed)\n",
      "    f_s:     0 (fixed)\n",
      "    h_1:     0.70028168 +/- 0.01532240 (2.19%) (init = 0.708)\n",
      "    h_2:     0.439 (fixed)\n",
      "    c:       1.00009341 +/- 8.5132e-06 (0.00%) (init = 1)\n",
      "    k:       0.02720288 +/- 0.00122603 (4.51%) == 'sqrt(D)'\n",
      "    bsq:     0.15050533 +/- 0.31741224 (210.90%) == '((1-k)**2-S*(1+k)**2)/(1-S)'\n",
      "    b:       0.38795017 +/- 0.40908893 (105.45%) == 'sqrt(((1-k)**2-S*(1+k)**2)/(1-S))'\n",
      "    aR:      21.2238985 +/- 0.22277322 (1.05%) == 'sqrt((1+k)**2-b**2)/W/pi'\n",
      "    logrho:  0.15369712 +/- 0.01367551 (8.90%) == 'log10(4.3275e-4*((1+k)**2-b**2)**1.5/W**3/P**2)'\n",
      "    e:       0.00000000 +/- 0.00000000  == 'f_c**2 + f_s**2'\n",
      "    q_1:     0.31472100 +/- 0.00000000 (0.00%) == '(1-h_2)**2'\n",
      "    q_2:     0.46574274 +/- 0.02731265 (5.86%) == '(h_1-h_2)/(1-h_2)'\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(D, S)   = -0.980\n",
      "    C(W, S)   = -0.737\n",
      "    C(D, W)   =  0.711\n",
      "    C(T_0, S) = -0.232\n",
      "    C(T_0, D) =  0.217\n",
      "    C(D, c)   =  0.136\n",
      "    C(W, h_1) = -0.128\n",
      "    C(W, c)   =  0.114\n",
      "[[Priors]]\n",
      "    T_0:    0.89000000 +/- 1.0000e-03\n",
      "    P:      9.48926400 +/- 1.0000e-05\n",
      "    h_1:    0.70800000 +/- 0.01000000\n",
      "    logrho: 0.05926000 +/- 0.15557188\n",
      "\n"
     ]
    }
   ],
   "source": [
    "result = dataset.lmfit_transit(P=ufloat(9.489264,0.00001), T_0=ufloat(0.890,0.001),\n",
    "                               h_1=pHD97658.h_1, h_2=pHD97658.h_2.n,\n",
    "                              logrhoprior=logrhoprior)\n",
    "print(dataset.lmfit_report())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 59456 characters]\n",
      "text/plain": [
       "<Figure size 648x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = dataset.plot_lmfit(title=dataset.target,figsize=(9,6),fontsize=12)\n",
    "fig.savefig(\"HD97658_lmfit.eps\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling the posterior\n",
    "Sampling the posterior probability distribution using the affine-invariant sampler [emcee](https://emcee.readthedocs.io/en/stable/). \n",
    "\n",
    "*Using a much lower number of burn-in and sampling steps here than would normally be used - just to see how it works*\n",
    "\n",
    "The $\\pm$1-sigma limits on parameters from the prior are indicated in the corner plot using dashed lines.\n",
    "\n",
    "This return value from emcee_transit is an instance of an [EnsembleSampler](https://emcee.readthedocs.io/en/latest/user/sampler/) object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running burn-in ..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " 50%|█████     | 16/32 [00:02<00:02,  7.06it/s]/Users/pflm/anaconda3/lib/python3.7/site-packages/asteval/asteval.py:776: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return func(*args, **keywords)\n",
      "100%|██████████| 32/32 [00:04<00:00,  7.03it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running sampler ..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "100%|██████████| 128/128 [00:18<00:00,  7.07it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Fit Statistics]]\n",
      "    # fitting method   = emcee\n",
      "    # function evals   = 4096\n",
      "    # data points      = 576\n",
      "    # variables        = 8\n",
      "    chi-square         = 1537.41443\n",
      "    reduced chi-square = 2.70671554\n",
      "    Akaike info crit   = -8352.22281\n",
      "    Bayesian info crit = -8317.37395\n",
      "[[Variables]]\n",
      "    T_0:        0.89102972 +/- 5.2278e-04 (0.06%) (init = 0.89)\n",
      "    P:          9.48924895 +/- 8.1155e-06 (0.00%) (init = 9.489264)\n",
      "    D:          8.0593e-04 +/- 5.3491e-05 (6.64%) (init = 0.001244335)\n",
      "    W:          0.01426424 +/- 1.5403e-04 (1.08%) (init = 0.02208278)\n",
      "    S:          0.83270305 +/- 0.04149386 (4.98%) (init = 0.8283019)\n",
      "    f_c:        0 (fixed)\n",
      "    f_s:        0 (fixed)\n",
      "    h_1:        0.69785553 +/- 0.00932562 (1.34%) (init = 0.708)\n",
      "    h_2:        0.439 (fixed)\n",
      "    c:          1.00009784 +/- 6.1286e-06 (0.00%) (init = 1)\n",
      "    k:          0.02838899 +/- 0.00122603 (4.32%) == 'sqrt(D)'\n",
      "    bsq:        0.37881499 +/- 0.31741224 (83.79%) == '((1-k)**2-S*(1+k)**2)/(1-S)'\n",
      "    b:          0.61547948 +/- 0.40908893 (66.47%) == 'sqrt(((1-k)**2-S*(1+k)**2)/(1-S))'\n",
      "    aR:         18.3849554 +/- 0.22277322 (1.21%) == 'sqrt((1+k)**2-b**2)/W/pi'\n",
      "    logrho:    -0.03338925 +/- 0.01367551 (40.96%) == 'log10(4.3275e-4*((1+k)**2-b**2)**1.5/W**3/P**2)'\n",
      "    e:          0.00000000 +/- 0.00000000  == 'f_c**2 + f_s**2'\n",
      "    q_1:        0.31472100 +/- 0.00000000 (0.00%) == '(1-h_2)**2'\n",
      "    q_2:        0.46141806 +/- 0.02731265 (5.92%) == '(h_1-h_2)/(1-h_2)'\n",
      "    log_sigma: -8.96257964 +/- 0.13115007 (1.46%) (init = -10)\n",
      "    sigma_w:    128.115335 == 'exp(log_sigma)*1e6'\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(D, S)         = -0.921\n",
      "    C(T_0, c)       =  0.725\n",
      "    C(W, S)         = -0.711\n",
      "    C(D, W)         =  0.597\n",
      "    C(P, W)         =  0.440\n",
      "    C(W, c)         = -0.427\n",
      "    C(T_0, P)       = -0.405\n",
      "    C(S, log_sigma) = -0.356\n",
      "    C(D, log_sigma) =  0.355\n",
      "    C(S, c)         =  0.324\n",
      "    C(h_1, c)       = -0.316\n",
      "    C(P, log_sigma) = -0.294\n",
      "    C(P, c)         = -0.290\n",
      "    C(P, D)         = -0.283\n",
      "    C(P, h_1)       = -0.240\n",
      "    C(P, S)         =  0.216\n",
      "    C(W, h_1)       = -0.204\n",
      "    C(D, c)         = -0.181\n",
      "    C(S, h_1)       =  0.163\n",
      "    C(D, h_1)       = -0.133\n",
      "    C(W, log_sigma) =  0.130\n",
      "    C(T_0, D)       =  0.126\n",
      "[[Priors]]\n",
      "    T_0:       0.89000000 +/- 1.0000e-03\n",
      "    P:         9.48926400 +/- 1.0000e-05\n",
      "    h_1:       0.70800000 +/- 0.01000000\n",
      "    logrho:    0.05926000 +/- 0.15557188\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "sampler = dataset.emcee_transit(burn=32, steps=32, nwalkers=32, thin=4)\n",
    "print(dataset.emcee_report())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 202240 characters]\n",
      "text/plain": [
       "<Figure size 849.6x849.6 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cornerplot = dataset.corner_plot(['k','logrho','T_0','S','h_1'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 62324 characters]\n",
      "text/plain": [
       "<Figure size 648x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = dataset.plot_emcee(title=dataset.target,figsize=(9,6),fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add correlated noise\n",
    "\n",
    "Now use \"add_shoterm=True\" so that correlated noise is modelled using the SHOTerm kernel from [celerite](https://celerite.readthedocs.io/en/stable/). By default the parameter quality factor is fixed at $Q=1/\\sqrt{2}$ so that the kernel is the one commonly used to model granulation noise\n",
    "$$ k(\\tau) = \\omega_0 e^{-\\omega_0\\tau/\\sqrt{2}}\\cos\\left(-\\omega_0\\tau/\\sqrt{2} - \\pi/4\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running burn-in ..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 128/128 [00:39<00:00,  3.20it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running sampler ..\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "100%|██████████| 256/256 [01:18<00:00,  3.35it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[Fit Statistics]]\n",
      "    # fitting method   = emcee\n",
      "    # function evals   = 16384\n",
      "    # data points      = 576\n",
      "    # variables        = 10\n",
      "    chi-square         = 1579.16358\n",
      "    reduced chi-square = 2.79004167\n",
      "    Akaike info crit   = -8343.27088\n",
      "    Bayesian info crit = -8299.70981\n",
      "[[Variables]]\n",
      "    T_0:         0.89092681 +/- 5.8843e-04 (0.07%) (init = 0.89)\n",
      "    P:           9.48926165 +/- 1.1423e-05 (0.00%) (init = 9.489264)\n",
      "    D:           8.2969e-04 +/- 9.8298e-05 (11.85%) (init = 0.001244335)\n",
      "    W:           0.01439629 +/- 1.7263e-04 (1.20%) (init = 0.02208278)\n",
      "    S:           0.83398169 +/- 0.03401786 (4.08%) (init = 0.8283019)\n",
      "    f_c:         0 (fixed)\n",
      "    f_s:         0 (fixed)\n",
      "    h_1:         0.70162736 +/- 0.01243761 (1.77%) (init = 0.708)\n",
      "    h_2:         0.439 (fixed)\n",
      "    c:           1.00008598 +/- 2.4863e-05 (0.00%) (init = 1)\n",
      "    k:           0.02880439 +/- 0.00122603 (4.26%) == 'sqrt(D)'\n",
      "    bsq:         0.36443327 +/- 0.31741224 (87.10%) == '((1-k)**2-S*(1+k)**2)/(1-S)'\n",
      "    b:           0.60368309 +/- 0.40908893 (67.77%) == 'sqrt(((1-k)**2-S*(1+k)**2)/(1-S))'\n",
      "    aR:          18.4196288 +/- 0.22277322 (1.21%) == 'sqrt((1+k)**2-b**2)/W/pi'\n",
      "    logrho:     -0.03093553 +/- 0.01367551 (44.21%) == 'log10(4.3275e-4*((1+k)**2-b**2)**1.5/W**3/P**2)'\n",
      "    e:           0.00000000 +/- 0.00000000  == 'f_c**2 + f_s**2'\n",
      "    q_1:         0.31472100 +/- 0.00000000 (0.00%) == '(1-h_2)**2'\n",
      "    q_2:         0.46814147 +/- 0.02731265 (5.83%) == '(h_1-h_2)/(1-h_2)'\n",
      "    log_S0:     -13.6267180 +/- 2.82540367 (20.73%) (init = -11)\n",
      "    log_omega0:  0.53973210 +/- 1.13765591 (210.78%) (init = 6)\n",
      "    log_Q:      -0.3465736 (fixed)\n",
      "    log_sigma:  -8.93508858 +/- 0.05635832 (0.63%) (init = -10)\n",
      "    sigma_w:     131.686221 == 'exp(log_sigma)*1e6'\n",
      "[[Correlations]] (unreported correlations are < 0.100)\n",
      "    C(W, S)                  = -0.522\n",
      "    C(D, S)                  = -0.479\n",
      "    C(log_S0, log_omega0)    = -0.442\n",
      "    C(D, W)                  =  0.341\n",
      "    C(T_0, W)                = -0.331\n",
      "    C(W, log_S0)             =  0.312\n",
      "    C(D, log_S0)             =  0.304\n",
      "    C(D, h_1)                = -0.271\n",
      "    C(S, log_S0)             = -0.265\n",
      "    C(c, log_omega0)         =  0.251\n",
      "    C(D, log_omega0)         =  0.243\n",
      "    C(P, D)                  = -0.207\n",
      "    C(S, c)                  =  0.205\n",
      "    C(c, log_S0)             = -0.196\n",
      "    C(h_1, log_omega0)       = -0.194\n",
      "    C(W, h_1)                = -0.153\n",
      "    C(T_0, D)                =  0.138\n",
      "    C(log_omega0, log_sigma) = -0.134\n",
      "    C(P, S)                  =  0.131\n",
      "    C(W, c)                  = -0.120\n",
      "    C(P, c)                  =  0.117\n",
      "    C(T_0, h_1)              = -0.111\n",
      "[[Priors]]\n",
      "    T_0:        0.89000000 +/- 1.0000e-03\n",
      "    P:          9.48926400 +/- 1.0000e-05\n",
      "    h_1:        0.70800000 +/- 0.01000000\n",
      "    logrho:     0.05926000 +/- 0.15557188\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "sampler = dataset.emcee_transit(add_shoterm=True, burn=128, steps=64, nwalkers=64, thin=4)\n",
    "print(dataset.emcee_report())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:Too few points to create valid contours\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 188828 characters]\n",
      "text/plain": [
       "<Figure size 849.6x849.6 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cornerplot = dataset.corner_plot(['log_sigma','log_omega0','log_S0','D','W'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 75712 characters]\n",
      "text/plain": [
       "<Figure size 648x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = dataset.plot_emcee(title=dataset.target,figsize=(9,6),fontsize=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "(c) 2019 Pierre Maxted, Keele Univerity, UK"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
