eb.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run notebook_setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Eclipsing binary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import exoplanet as xo\n",
    "import lightkurve as lk\n",
    "import pymc3 as pm\n",
    "import theano.tensor as tt\n",
    "import corner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lit_period = 4.541710\n",
    "lit_t0 = 896.19699 + 2456000 - 2454833"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lcf = lk.search_lightcurvefile(\"EPIC 203868608\").download(quality_bitmask=0)\n",
    "\n",
    "hdr = lcf.header(1)\n",
    "texp = hdr[\"FRAMETIM\"] * hdr[\"NUM_FRM\"]\n",
    "texp /= 60.0 * 60.0 * 24.0\n",
    "\n",
    "lc = lcf.PDCSAP_FLUX\n",
    "lc = lc.remove_nans()\n",
    "\n",
    "# lc = lc.flatten(window_length=401)\n",
    "\n",
    "# lc = lc.to_corrector(\"sff\").correct(windows=20)\n",
    "\n",
    "lc = lc.normalize()\n",
    "\n",
    "x = np.ascontiguousarray(lc.time, dtype=np.float64)\n",
    "y = np.ascontiguousarray(1e3 * (lc.flux - 1), dtype=np.float64)\n",
    "yerr = np.ascontiguousarray(1e3 * lc.flux_err, dtype=np.float64)\n",
    "\n",
    "# plt.plot(x, y)\n",
    "plt.scatter((x - lit_t0 + 0.5 * lit_period) % lit_period - 0.5 * lit_period, y, c=x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_fold = (x - lit_t0 + 0.5 * lit_period) % lit_period - 0.5 * lit_period\n",
    "m = np.abs(x_fold) > 0.1\n",
    "bins = np.linspace(-0.5 * lit_period, 0.5 * lit_period, 200)\n",
    "\n",
    "num, _ = np.histogram(x_fold[m], bins, weights=y[m])\n",
    "denom, _ = np.histogram(x_fold[m], bins)\n",
    "denom[denom == 0] = 1.0\n",
    "\n",
    "ind = np.argmin(num / denom)\n",
    "dt = 0.5 * (bins[ind] + bins[ind + 1])\n",
    "ecosw = 0.25 * np.pi * (2 * dt / lit_period - 1)\n",
    "sqrtecosw_guess = np.sign(ecosw) * np.sqrt(np.abs(ecosw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_model(mask):\n",
    "\n",
    "    with pm.Model() as model:\n",
    "\n",
    "        mean = pm.Normal(\"mean\", mu=0.0, sd=5.0)\n",
    "        flux_ratio = pm.Lognormal(\"flux_ratio\", mu=np.log(50.0), sigma=10.0)\n",
    "        \n",
    "        u1 = xo.QuadLimbDark(\"u1\")\n",
    "        u2 = xo.QuadLimbDark(\"u2\")\n",
    "        \n",
    "        # Following Maxted (2016):\n",
    "        #  r1 = (R1 + R2) / a; r2 = R2 / R1  ->  a / R1 = (1 + r2) / r1\n",
    "        r1 = pm.Lognormal(\"r1\", mu=np.log(0.2), sigma=50)\n",
    "        r2 = pm.Lognormal(\"r2\", mu=np.log(1.0), sigma=50)\n",
    "        a = pm.Deterministic(\"a\", (1 + r2) / r1)\n",
    "\n",
    "        period = pm.Lognormal(\"period\", mu=np.log(lit_period), sigma=1.0)\n",
    "        t0 = pm.Normal(\"t0\", mu=lit_t0, sigma=1.0)\n",
    "        b = xo.ImpactParameter(\"b\", ror=r2, testval=0.5)\n",
    "                   \n",
    "        hk = xo.UnitDisk(\"hk\", testval=np.array([sqrtecosw_guess, 0.0]))\n",
    "        ecc = pm.Deterministic(\"ecc\", tt.sum(hk ** 2))\n",
    "        omega = pm.Deterministic(\"omega\", tt.arctan2(hk[1], hk[0]))\n",
    "\n",
    "        orbit = xo.orbits.KeplerianOrbit(period=period, t0=t0, ecc=ecc, omega=omega, b=b, a=a)\n",
    "\n",
    "        lc = xo.SecondaryEclipseLightCurve(u1, u2, 1e-3 * flux_ratio)\n",
    "\n",
    "        def lc_model(t):\n",
    "            return 1e3 * lc.get_light_curve(orbit=orbit, r=r2, t=t, texp=texp)[:, 0]\n",
    "\n",
    "        def full_model(t):\n",
    "            return mean + lc_model(t)\n",
    "\n",
    "        sigma = pm.InverseGamma(\"sigma\", testval=1.0, **xo.estimate_inverse_gamma_parameters(0.1, 2.0))\n",
    "        S_tot = pm.InverseGamma(\"S_tot\", testval=2.5, **xo.estimate_inverse_gamma_parameters(1.0, 5.0))\n",
    "        ell = pm.InverseGamma(\"ell\", testval=2.0, **xo.estimate_inverse_gamma_parameters(1.0, 5.0))\n",
    "        kernel = xo.gp.terms.SHOTerm(S_tot=S_tot, w0=2 * np.pi / ell, Q=1. / 3)\n",
    "        gp = xo.gp.GP(kernel, x[mask], yerr[mask] ** 2 + sigma ** 2, mean=full_model)\n",
    "\n",
    "        gp.marginal(\"obs\", observed=y[mask])\n",
    "        \n",
    "        map_soln = model.test_point\n",
    "        map_soln = xo.optimize(map_soln, [mean, u1, u2, r1, r2, flux_ratio, b])\n",
    "        map_soln = xo.optimize(map_soln, [mean, sigma, S_tot, ell])\n",
    "        map_soln = xo.optimize(map_soln, [t0, period])\n",
    "        map_soln = xo.optimize(map_soln, [hk])\n",
    "        map_soln = xo.optimize(map_soln)\n",
    "\n",
    "        model.gp = gp\n",
    "        model.lc_model = lc_model\n",
    "        model.full_model = full_model\n",
    "        model.x = x[mask]\n",
    "        model.y = y[mask]\n",
    "        model.yerr = yerr[mask]\n",
    "        \n",
    "    return model, map_soln\n",
    "\n",
    "def sigma_clip():\n",
    "    mask = np.ones(len(x), dtype=bool)\n",
    "    num = len(mask)\n",
    "    \n",
    "    for i in range(10):\n",
    "        model, map_soln = build_model(mask)\n",
    "\n",
    "        with model:\n",
    "            mdl = xo.eval_in_model(model.full_model(x[mask]) + model.gp.predict(), map_soln)\n",
    "\n",
    "        resid = y[mask] - mdl\n",
    "        sigma = np.sqrt(np.median((resid - np.median(resid)) ** 2))\n",
    "        mask[mask] = np.abs(resid - np.median(resid)) < 7 * sigma \n",
    "        print(num, mask.sum())\n",
    "        if num == mask.sum():\n",
    "            break\n",
    "        num = mask.sum()\n",
    "        \n",
    "    return model, map_soln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model, map_soln = sigma_clip()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_soln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with model:    \n",
    "    gp_pred = xo.eval_in_model(model.gp.predict(), map_soln) + map_soln[\"mean\"]\n",
    "    lc = xo.eval_in_model(model.lc_model(model.x), map_soln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2, sharex=True, figsize=(12, 7))\n",
    "\n",
    "ax1.plot(model.x, model.y, \"k.\", alpha=0.2)\n",
    "ax1.plot(model.x, gp_pred, color=\"C1\", lw=1)\n",
    "\n",
    "ax2.plot(model.x, model.y - gp_pred, \"k.\", alpha=0.2)\n",
    "ax2.plot(model.x, lc, color=\"C2\", lw=1)\n",
    "ax2.set_xlim(model.x.min(), model.x.max())\n",
    "\n",
    "ax1.set_title(\"EPIC 203868608\", fontsize=14)\n",
    "ax1.set_ylabel(\"raw flux [ppt]\")\n",
    "ax2.set_ylabel(\"de-trended flux [ppt]\")\n",
    "ax2.set_xlabel(\"time [KBJD]\")\n",
    "\n",
    "fig.subplots_adjust(hspace=0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots(1, figsize=(12, 3.5))\n",
    "\n",
    "x_fold = (model.x - map_soln[\"t0\"]) % map_soln[\"period\"] / map_soln[\"period\"]\n",
    "inds = np.argsort(x_fold)\n",
    "\n",
    "ax1.plot(x_fold[inds], model.y[inds] - gp_pred[inds], \"k.\", alpha=0.2)\n",
    "ax1.plot(x_fold[inds] - 1, model.y[inds] - gp_pred[inds], \"k.\", alpha=0.2)\n",
    "ax2.plot(x_fold[inds], model.y[inds] - gp_pred[inds], \"k.\", alpha=0.2, label=\"data!\")\n",
    "ax2.plot(x_fold[inds] - 1, model.y[inds] - gp_pred, \"k.\", alpha=0.2)\n",
    "\n",
    "yval = model.y[inds] - gp_pred\n",
    "bins = np.linspace(0, 1, 75)\n",
    "num, _ = np.histogram(x_fold[inds], bins, weights=yval)\n",
    "denom, _ = np.histogram(x_fold[inds], bins)\n",
    "ax2.plot(0.5 * (bins[:-1] + bins[1:]) - 1, num / denom, \".w\")\n",
    "\n",
    "args = dict(lw=1)\n",
    "\n",
    "ax1.plot(x_fold[inds], lc[inds], \"C1\", **args)\n",
    "ax1.plot(x_fold[inds] - 1, lc[inds], \"C1\", **args)\n",
    "    \n",
    "ax1.set_xlim(-1, 1)\n",
    "ax1.set_ylabel(\"de-trended flux [ppt]\")\n",
    "ax1.set_xlabel(\"phase\")\n",
    "ax1.set_title(\"EPIC 203868608\", fontsize=14);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(203868608)\n",
    "with model:\n",
    "    trace = pm.sample(\n",
    "        tune=2000,\n",
    "        draws=2000,\n",
    "        start=map_soln,\n",
    "        cores=2,\n",
    "        chains=2,\n",
    "        step=xo.get_dense_nuts_step(start=map_soln, target_accept=0.9),\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pm.summary(trace)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = pm.trace_to_dataframe(trace, varnames=[\"flux_ratio\", \"r1\", \"r2\", \"b\", \"ecc\"])\n",
    "samples[\"flux_ratio\"] *= 1e-3\n",
    "m = samples[\"r2\"] < 2.0\n",
    "corner.corner(samples[m]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = pm.trace_to_dataframe(trace, varnames=[\"hk\", \"ecc\", \"omega\"])\n",
    "corner.corner(samples, labels=[\"$\\sqrt{e}\\sin\\omega$\", \"$\\sqrt{e}\\cos\\omega$\", \"$e$\", \"$\\omega$\"]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}