# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.6.0
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + nbsphinx="hidden"
import lightkurve as lk

# %matplotlib inline

# + nbsphinx="hidden"
# %run notebook_setup
# -

# # Fitting TESS data

# In this tutorial, we will reproduce the fits to the transiting planet in the Pi Mensae system discovered by [Huang et al. (2018)](https://arxiv.org/abs/1809.05967).
# The data processing and model are similar to the [Putting it all together](./together.ipynb) case study, but with a few extra bits like aperture selection and de-trending.
#
# To start, we need to download the target pixel file:

# +
import numpy as np
import lightkurve as lk
import matplotlib.pyplot as plt

lc_file = lk.search_lightcurvefile("TIC 261136679", sector=1).download(
    quality_bitmask="hardest"
)
lc = lc_file.PDCSAP_FLUX.remove_nans().normalize().remove_outliers()
time = lc.time
flux = lc.flux
m = lc.quality == 0
with lc_file.hdu as hdu:
    hdr = hdu[1].header

texp = hdr["FRAMETIM"] * hdr["NUM_FRM"]
texp /= 60.0 * 60.0 * 24.0

ref_time = 0.5 * (np.min(time) + np.max(time))
x = np.ascontiguousarray(time[m] - ref_time, dtype=np.float64)
y = np.ascontiguousarray(1e3 * (flux[m] - 1.0), dtype=np.float64)

plt.plot(x, y, ".k")
plt.xlabel("time [days]")
plt.ylabel("relative flux [ppt]")
_ = plt.xlim(x.min(), x.max())
# -

# ## Transit search
#
# Now, let's use [the box least squares periodogram from AstroPy](http://docs.astropy.org/en/latest/timeseries/bls.html)
# (Note: you'll need AstroPy v3.1 or more recent to use this feature) to estimate the period, phase, and depth of the transit.

# +
from astropy.timeseries import BoxLeastSquares

period_grid = np.exp(np.linspace(np.log(1), np.log(15), 50000))

bls = BoxLeastSquares(x, y)
bls_power = bls.power(period_grid, 0.1, oversample=20)

# Save the highest peak as the planet candidate
index = np.argmax(bls_power.power)
bls_period = bls_power.period[index]
bls_t0 = bls_power.transit_time[index]
bls_depth = bls_power.depth[index]
transit_mask = bls.transit_mask(x, bls_period, 0.2, bls_t0)

fig, axes = plt.subplots(2, 1, figsize=(10, 10))

# Plot the periodogram
ax = axes[0]
ax.axvline(np.log10(bls_period), color="C1", lw=5, alpha=0.8)
ax.plot(np.log10(bls_power.period), bls_power.power, "k")
ax.annotate(
    "period = {0:.4f} d".format(bls_period),
    (0, 1),
    xycoords="axes fraction",
    xytext=(5, -5),
    textcoords="offset points",
    va="top",
    ha="left",
    fontsize=12,
)
ax.set_ylabel("bls power")
ax.set_yticks([])
ax.set_xlim(np.log10(period_grid.min()), np.log10(period_grid.max()))
ax.set_xlabel("log10(period)")

# Plot the folded transit
ax = axes[1]
x_fold = (x - bls_t0 + 0.5 * bls_period) % bls_period - 0.5 * bls_period
m = np.abs(x_fold) < 0.4
ax.plot(x_fold[m], y[m], ".k")

# Overplot the phase binned light curve
bins = np.linspace(-0.41, 0.41, 32)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y)
denom[num == 0] = 1.0
ax.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, color="C1")

ax.set_xlim(-0.3, 0.3)
ax.set_ylabel("de-trended flux [ppt]")
_ = ax.set_xlabel("time since transit")
# -

# ## The transit model in PyMC3
#
# The transit model, initialization, and sampling are all nearly the same as the one in [Putting it all together](./together.ipynb).

# +
import exoplanet as xo
import pymc3 as pm
import theano.tensor as tt

import pymc3_ext as pmx
from celerite2.theano import terms, GaussianProcess


def build_model(mask=None, start=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    with pm.Model() as model:

        # Parameters for the stellar properties
        mean = pm.Normal("mean", mu=0.0, sd=10.0)
        u_star = xo.QuadLimbDark("u_star")

        # Stellar parameters from Huang et al (2018)
        M_star_huang = 1.094, 0.039
        R_star_huang = 1.10, 0.023
        BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=3)
        m_star = BoundedNormal(
            "m_star", mu=M_star_huang[0], sd=M_star_huang[1]
        )
        r_star = BoundedNormal(
            "r_star", mu=R_star_huang[0], sd=R_star_huang[1]
        )

        # Orbital parameters for the planets
        period = pm.Lognormal("period", mu=np.log(bls_period), sd=1)
        t0 = pm.Normal("t0", mu=bls_t0, sd=1)
        r_pl = pm.Lognormal(
            "r_pl",
            sd=1.0,
            mu=0.5 * np.log(1e-3 * np.array(bls_depth))
            + np.log(R_star_huang[0]),
        )
        ror = pm.Deterministic("ror", r_pl / r_star)
        b = xo.distributions.ImpactParameter("b", ror=ror)

        ecs = pmx.UnitDisk("ecs", testval=np.array([0.01, 0.0]))
        ecc = pm.Deterministic("ecc", tt.sum(ecs ** 2))
        omega = pm.Deterministic("omega", tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.kipping13("ecc_prior", fixed=True, observed=ecc)

        # Transit jitter & GP parameters
        sigma_lc = pm.Lognormal("sigma_lc", mu=np.log(np.std(y[mask])), sd=10)
        rho_gp = pm.Lognormal("rho_gp", mu=0, sd=10)
        sigma_gp = pm.Lognormal("sigma_gp", mu=np.log(np.std(y[mask])), sd=10)

        # Orbit model
        orbit = xo.orbits.KeplerianOrbit(
            r_star=r_star,
            m_star=m_star,
            period=period,
            t0=t0,
            b=b,
            ecc=ecc,
            omega=omega,
        )

        # Compute the model light curve
        light_curves = pm.Deterministic(
            "light_curves",
            xo.LimbDarkLightCurve(u_star).get_light_curve(
                orbit=orbit, r=r_pl, t=x[mask], texp=texp
            )
            * 1e3,
        )
        light_curve = tt.sum(light_curves, axis=-1) + mean
        resid = y[mask] - light_curve

        # GP model for the light curve
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1 / np.sqrt(2))
        gp = GaussianProcess(kernel, t=x[mask], yerr=sigma_lc)
        gp.marginal("gp", observed=resid)
        pm.Deterministic("gp_pred", gp.predict(resid))

        # Fit for the maximum a posteriori parameters, I've found that I can get
        # a better solution by trying different combinations of parameters in turn
        if start is None:
            start = model.test_point
        map_soln = pmx.optimize(start=start, vars=[sigma_lc, sigma_gp, rho_gp])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[period, t0])
        map_soln = pmx.optimize(start=map_soln, vars=[u_star])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[ecs])
        map_soln = pmx.optimize(start=map_soln, vars=[mean])
        map_soln = pmx.optimize(
            start=map_soln, vars=[sigma_lc, sigma_gp, rho_gp]
        )
        map_soln = pmx.optimize(start=map_soln)

    return model, map_soln


model0, map_soln0 = build_model()


# -

# Here's how we plot the initial light curve model:

# +
def plot_light_curve(soln, mask=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)

    fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)

    ax = axes[0]
    ax.plot(x[mask], y[mask], "k", label="data")
    gp_mod = soln["gp_pred"] + soln["mean"]
    ax.plot(x[mask], gp_mod, color="C2", label="gp model")
    ax.legend(fontsize=10)
    ax.set_ylabel("relative flux [ppt]")

    ax = axes[1]
    ax.plot(x[mask], y[mask] - gp_mod, "k", label="de-trended data")
    for i, l in enumerate("b"):
        mod = soln["light_curves"][:, i]
        ax.plot(x[mask], mod, label="planet {0}".format(l))
    ax.legend(fontsize=10, loc=3)
    ax.set_ylabel("de-trended flux [ppt]")

    ax = axes[2]
    mod = gp_mod + np.sum(soln["light_curves"], axis=-1)
    ax.plot(x[mask], y[mask] - mod, "k")
    ax.axhline(0, color="#aaaaaa", lw=1)
    ax.set_ylabel("residuals [ppt]")
    ax.set_xlim(x[mask].min(), x[mask].max())
    ax.set_xlabel("time [days]")

    return fig


_ = plot_light_curve(map_soln0)
# -

# As in [Putting it all together](./together.ipynb), we can do some sigma clipping to remove significant outliers.

# +
mod = (
    map_soln0["gp_pred"]
    + map_soln0["mean"]
    + np.sum(map_soln0["light_curves"], axis=-1)
)
resid = y - mod
rms = np.sqrt(np.median(resid ** 2))
mask = np.abs(resid) < 5 * rms

plt.figure(figsize=(10, 5))
plt.plot(x, resid, "k", label="data")
plt.plot(x[~mask], resid[~mask], "xr", label="outliers")
plt.axhline(0, color="#aaaaaa", lw=1)
plt.ylabel("residuals [ppt]")
plt.xlabel("time [days]")
plt.legend(fontsize=12, loc=3)
_ = plt.xlim(x.min(), x.max())
# -

# And then we re-build the model using the data without outliers.

model, map_soln = build_model(mask, map_soln0)
_ = plot_light_curve(map_soln, mask)

# Now that we have the model, we can sample:

np.random.seed(261136679)
with model:
    trace = pmx.sample(
        tune=2500,
        draws=2000,
        start=map_soln,
        cores=2,
        chains=2,
        initial_accept=0.8,
        target_accept=0.95,
    )

with model:
    summary = pm.summary(
        trace,
        var_names=[
            "omega",
            "ecc",
            "r_pl",
            "b",
            "t0",
            "period",
            "r_star",
            "m_star",
            "u_star",
            "mean",
        ],
    )
summary

# ## Results
#
# After sampling, we can make the usual plots.
# First, let's look at the folded light curve plot:

# +
# Compute the GP prediction
gp_mod = np.median(trace["gp_pred"] + trace["mean"][:, None], axis=0)

# Get the posterior median orbital parameters
p = np.median(trace["period"])
t0 = np.median(trace["t0"])

# Plot the folded data
x_fold = (x[mask] - t0 + 0.5 * p) % p - 0.5 * p
plt.plot(x_fold, y[mask] - gp_mod, ".k", label="data", zorder=-1000)

# Overplot the phase binned light curve
bins = np.linspace(-0.41, 0.41, 50)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y[mask])
denom[num == 0] = 1.0
plt.plot(
    0.5 * (bins[1:] + bins[:-1]), num / denom, "o", color="C2", label="binned"
)

# Plot the folded model
inds = np.argsort(x_fold)
inds = inds[np.abs(x_fold)[inds] < 0.3]
pred = trace["light_curves"][:, inds, 0]
pred = np.percentile(pred, [16, 50, 84], axis=0)
plt.plot(x_fold[inds], pred[1], color="C1", label="model")
art = plt.fill_between(
    x_fold[inds], pred[0], pred[2], color="C1", alpha=0.5, zorder=1000
)
art.set_edgecolor("none")

# Annotate the plot with the planet's period
txt = "period = {0:.5f} +/- {1:.5f} d".format(
    np.mean(trace["period"]), np.std(trace["period"])
)
plt.annotate(
    txt,
    (0, 0),
    xycoords="axes fraction",
    xytext=(5, 5),
    textcoords="offset points",
    ha="left",
    va="bottom",
    fontsize=12,
)

plt.legend(fontsize=10, loc=4)
plt.xlim(-0.5 * p, 0.5 * p)
plt.xlabel("time since transit [days]")
plt.ylabel("de-trended flux")
_ = plt.xlim(-0.15, 0.15)
# -

# And a corner plot of some of the key parameters:

# +
import corner
import astropy.units as u

varnames = ["period", "b", "ecc", "r_pl"]
samples = pm.trace_to_dataframe(trace, varnames=varnames)

# Convert the radius to Earth radii
samples["r_pl"] = (np.array(samples["r_pl"]) * u.R_sun).to(u.R_earth).value

_ = corner.corner(
    samples[["period", "r_pl", "b", "ecc"]],
    labels=[
        "period [days]",
        "radius [Earth radii]",
        "impact param",
        "eccentricity",
    ],
)
# -

# These all seem consistent with the previously published values.

# ## Citations
#
# As described in the [citation tutorial](https://docs.exoplanet.codes/en/stable/tutorials/citation/), we can use [citations.get_citations_for_model](https://docs.exoplanet.codes/en/stable/user/api/#exoplanet.citations.get_citations_for_model) to construct an acknowledgement and BibTeX listing that includes the relevant citations for this model.

with model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)

print("\n".join(bib.splitlines()[:10]) + "\n...")

348fbbfab1bb2b8ea575e848a692748e45017faf
<mask token>


def build_model(mask=None, start=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    with pm.Model() as model:
        mean = pm.Normal('mean', mu=0.0, sd=10.0)
        u_star = xo.QuadLimbDark('u_star')
        M_star_huang = 1.094, 0.039
        R_star_huang = 1.1, 0.023
        BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=3)
        m_star = BoundedNormal('m_star', mu=M_star_huang[0], sd=M_star_huang[1]
            )
        r_star = BoundedNormal('r_star', mu=R_star_huang[0], sd=R_star_huang[1]
            )
        period = pm.Lognormal('period', mu=np.log(bls_period), sd=1)
        t0 = pm.Normal('t0', mu=bls_t0, sd=1)
        r_pl = pm.Lognormal('r_pl', sd=1.0, mu=0.5 * np.log(0.001 * np.
            array(bls_depth)) + np.log(R_star_huang[0]))
        ror = pm.Deterministic('ror', r_pl / r_star)
        b = xo.distributions.ImpactParameter('b', ror=ror)
        ecs = pmx.UnitDisk('ecs', testval=np.array([0.01, 0.0]))
        ecc = pm.Deterministic('ecc', tt.sum(ecs ** 2))
        omega = pm.Deterministic('omega', tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.kipping13('ecc_prior', fixed=True, observed=ecc)
        sigma_lc = pm.Lognormal('sigma_lc', mu=np.log(np.std(y[mask])), sd=10)
        rho_gp = pm.Lognormal('rho_gp', mu=0, sd=10)
        sigma_gp = pm.Lognormal('sigma_gp', mu=np.log(np.std(y[mask])), sd=10)
        orbit = xo.orbits.KeplerianOrbit(r_star=r_star, m_star=m_star,
            period=period, t0=t0, b=b, ecc=ecc, omega=omega)
        light_curves = pm.Deterministic('light_curves', xo.
            LimbDarkLightCurve(u_star).get_light_curve(orbit=orbit, r=r_pl,
            t=x[mask], texp=texp) * 1000.0)
        light_curve = tt.sum(light_curves, axis=-1) + mean
        resid = y[mask] - light_curve
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1 / np.sqrt(2))
        gp = GaussianProcess(kernel, t=x[mask], yerr=sigma_lc)
        gp.marginal('gp', observed=resid)
        pm.Deterministic('gp_pred', gp.predict(resid))
        if start is None:
            start = model.test_point
        map_soln = pmx.optimize(start=start, vars=[sigma_lc, sigma_gp, rho_gp])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[period, t0])
        map_soln = pmx.optimize(start=map_soln, vars=[u_star])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[ecs])
        map_soln = pmx.optimize(start=map_soln, vars=[mean])
        map_soln = pmx.optimize(start=map_soln, vars=[sigma_lc, sigma_gp,
            rho_gp])
        map_soln = pmx.optimize(start=map_soln)
    return model, map_soln


<mask token>


def plot_light_curve(soln, mask=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
    ax = axes[0]
    ax.plot(x[mask], y[mask], 'k', label='data')
    gp_mod = soln['gp_pred'] + soln['mean']
    ax.plot(x[mask], gp_mod, color='C2', label='gp model')
    ax.legend(fontsize=10)
    ax.set_ylabel('relative flux [ppt]')
    ax = axes[1]
    ax.plot(x[mask], y[mask] - gp_mod, 'k', label='de-trended data')
    for i, l in enumerate('b'):
        mod = soln['light_curves'][:, i]
        ax.plot(x[mask], mod, label='planet {0}'.format(l))
    ax.legend(fontsize=10, loc=3)
    ax.set_ylabel('de-trended flux [ppt]')
    ax = axes[2]
    mod = gp_mod + np.sum(soln['light_curves'], axis=-1)
    ax.plot(x[mask], y[mask] - mod, 'k')
    ax.axhline(0, color='#aaaaaa', lw=1)
    ax.set_ylabel('residuals [ppt]')
    ax.set_xlim(x[mask].min(), x[mask].max())
    ax.set_xlabel('time [days]')
    return fig


<mask token>

<mask token>
with lc_file.hdu as hdu:
    hdr = hdu[1].header
<mask token>
texp /= 60.0 * 60.0 * 24.0
<mask token>
plt.plot(x, y, '.k')
plt.xlabel('time [days]')
plt.ylabel('relative flux [ppt]')
<mask token>
ax.axvline(np.log10(bls_period), color='C1', lw=5, alpha=0.8)
ax.plot(np.log10(bls_power.period), bls_power.power, 'k')
ax.annotate('period = {0:.4f} d'.format(bls_period), (0, 1), xycoords=
    'axes fraction', xytext=(5, -5), textcoords='offset points', va='top',
    ha='left', fontsize=12)
ax.set_ylabel('bls power')
ax.set_yticks([])
ax.set_xlim(np.log10(period_grid.min()), np.log10(period_grid.max()))
ax.set_xlabel('log10(period)')
<mask token>
ax.plot(x_fold[m], y[m], '.k')
<mask token>
ax.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, color='C1')
ax.set_xlim(-0.3, 0.3)
ax.set_ylabel('de-trended flux [ppt]')
<mask token>


def build_model(mask=None, start=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    with pm.Model() as model:
        mean = pm.Normal('mean', mu=0.0, sd=10.0)
        u_star = xo.QuadLimbDark('u_star')
        M_star_huang = 1.094, 0.039
        R_star_huang = 1.1, 0.023
        BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=3)
        m_star = BoundedNormal('m_star', mu=M_star_huang[0], sd=M_star_huang[1]
            )
        r_star = BoundedNormal('r_star', mu=R_star_huang[0], sd=R_star_huang[1]
            )
        period = pm.Lognormal('period', mu=np.log(bls_period), sd=1)
        t0 = pm.Normal('t0', mu=bls_t0, sd=1)
        r_pl = pm.Lognormal('r_pl', sd=1.0, mu=0.5 * np.log(0.001 * np.
            array(bls_depth)) + np.log(R_star_huang[0]))
        ror = pm.Deterministic('ror', r_pl / r_star)
        b = xo.distributions.ImpactParameter('b', ror=ror)
        ecs = pmx.UnitDisk('ecs', testval=np.array([0.01, 0.0]))
        ecc = pm.Deterministic('ecc', tt.sum(ecs ** 2))
        omega = pm.Deterministic('omega', tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.kipping13('ecc_prior', fixed=True, observed=ecc)
        sigma_lc = pm.Lognormal('sigma_lc', mu=np.log(np.std(y[mask])), sd=10)
        rho_gp = pm.Lognormal('rho_gp', mu=0, sd=10)
        sigma_gp = pm.Lognormal('sigma_gp', mu=np.log(np.std(y[mask])), sd=10)
        orbit = xo.orbits.KeplerianOrbit(r_star=r_star, m_star=m_star,
            period=period, t0=t0, b=b, ecc=ecc, omega=omega)
        light_curves = pm.Deterministic('light_curves', xo.
            LimbDarkLightCurve(u_star).get_light_curve(orbit=orbit, r=r_pl,
            t=x[mask], texp=texp) * 1000.0)
        light_curve = tt.sum(light_curves, axis=-1) + mean
        resid = y[mask] - light_curve
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1 / np.sqrt(2))
        gp = GaussianProcess(kernel, t=x[mask], yerr=sigma_lc)
        gp.marginal('gp', observed=resid)
        pm.Deterministic('gp_pred', gp.predict(resid))
        if start is None:
            start = model.test_point
        map_soln = pmx.optimize(start=start, vars=[sigma_lc, sigma_gp, rho_gp])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[period, t0])
        map_soln = pmx.optimize(start=map_soln, vars=[u_star])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[ecs])
        map_soln = pmx.optimize(start=map_soln, vars=[mean])
        map_soln = pmx.optimize(start=map_soln, vars=[sigma_lc, sigma_gp,
            rho_gp])
        map_soln = pmx.optimize(start=map_soln)
    return model, map_soln


<mask token>


def plot_light_curve(soln, mask=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
    ax = axes[0]
    ax.plot(x[mask], y[mask], 'k', label='data')
    gp_mod = soln['gp_pred'] + soln['mean']
    ax.plot(x[mask], gp_mod, color='C2', label='gp model')
    ax.legend(fontsize=10)
    ax.set_ylabel('relative flux [ppt]')
    ax = axes[1]
    ax.plot(x[mask], y[mask] - gp_mod, 'k', label='de-trended data')
    for i, l in enumerate('b'):
        mod = soln['light_curves'][:, i]
        ax.plot(x[mask], mod, label='planet {0}'.format(l))
    ax.legend(fontsize=10, loc=3)
    ax.set_ylabel('de-trended flux [ppt]')
    ax = axes[2]
    mod = gp_mod + np.sum(soln['light_curves'], axis=-1)
    ax.plot(x[mask], y[mask] - mod, 'k')
    ax.axhline(0, color='#aaaaaa', lw=1)
    ax.set_ylabel('residuals [ppt]')
    ax.set_xlim(x[mask].min(), x[mask].max())
    ax.set_xlabel('time [days]')
    return fig


<mask token>
plt.figure(figsize=(10, 5))
plt.plot(x, resid, 'k', label='data')
plt.plot(x[~mask], resid[~mask], 'xr', label='outliers')
plt.axhline(0, color='#aaaaaa', lw=1)
plt.ylabel('residuals [ppt]')
plt.xlabel('time [days]')
plt.legend(fontsize=12, loc=3)
<mask token>
np.random.seed(261136679)
with model:
    trace = pmx.sample(tune=2500, draws=2000, start=map_soln, cores=2,
        chains=2, initial_accept=0.8, target_accept=0.95)
with model:
    summary = pm.summary(trace, var_names=['omega', 'ecc', 'r_pl', 'b',
        't0', 'period', 'r_star', 'm_star', 'u_star', 'mean'])
summary
<mask token>
plt.plot(x_fold, y[mask] - gp_mod, '.k', label='data', zorder=-1000)
<mask token>
plt.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, 'o', color='C2', label=
    'binned')
<mask token>
plt.plot(x_fold[inds], pred[1], color='C1', label='model')
<mask token>
art.set_edgecolor('none')
<mask token>
plt.annotate(txt, (0, 0), xycoords='axes fraction', xytext=(5, 5),
    textcoords='offset points', ha='left', va='bottom', fontsize=12)
plt.legend(fontsize=10, loc=4)
plt.xlim(-0.5 * p, 0.5 * p)
plt.xlabel('time since transit [days]')
plt.ylabel('de-trended flux')
<mask token>
with model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)
print('\n'.join(bib.splitlines()[:10]) + '\n...')

<mask token>
lc_file = lk.search_lightcurvefile('TIC 261136679', sector=1).download(
    quality_bitmask='hardest')
lc = lc_file.PDCSAP_FLUX.remove_nans().normalize().remove_outliers()
time = lc.time
flux = lc.flux
m = lc.quality == 0
with lc_file.hdu as hdu:
    hdr = hdu[1].header
texp = hdr['FRAMETIM'] * hdr['NUM_FRM']
texp /= 60.0 * 60.0 * 24.0
ref_time = 0.5 * (np.min(time) + np.max(time))
x = np.ascontiguousarray(time[m] - ref_time, dtype=np.float64)
y = np.ascontiguousarray(1000.0 * (flux[m] - 1.0), dtype=np.float64)
plt.plot(x, y, '.k')
plt.xlabel('time [days]')
plt.ylabel('relative flux [ppt]')
_ = plt.xlim(x.min(), x.max())
<mask token>
period_grid = np.exp(np.linspace(np.log(1), np.log(15), 50000))
bls = BoxLeastSquares(x, y)
bls_power = bls.power(period_grid, 0.1, oversample=20)
index = np.argmax(bls_power.power)
bls_period = bls_power.period[index]
bls_t0 = bls_power.transit_time[index]
bls_depth = bls_power.depth[index]
transit_mask = bls.transit_mask(x, bls_period, 0.2, bls_t0)
fig, axes = plt.subplots(2, 1, figsize=(10, 10))
ax = axes[0]
ax.axvline(np.log10(bls_period), color='C1', lw=5, alpha=0.8)
ax.plot(np.log10(bls_power.period), bls_power.power, 'k')
ax.annotate('period = {0:.4f} d'.format(bls_period), (0, 1), xycoords=
    'axes fraction', xytext=(5, -5), textcoords='offset points', va='top',
    ha='left', fontsize=12)
ax.set_ylabel('bls power')
ax.set_yticks([])
ax.set_xlim(np.log10(period_grid.min()), np.log10(period_grid.max()))
ax.set_xlabel('log10(period)')
ax = axes[1]
x_fold = (x - bls_t0 + 0.5 * bls_period) % bls_period - 0.5 * bls_period
m = np.abs(x_fold) < 0.4
ax.plot(x_fold[m], y[m], '.k')
bins = np.linspace(-0.41, 0.41, 32)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y)
denom[num == 0] = 1.0
ax.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, color='C1')
ax.set_xlim(-0.3, 0.3)
ax.set_ylabel('de-trended flux [ppt]')
_ = ax.set_xlabel('time since transit')
<mask token>


def build_model(mask=None, start=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    with pm.Model() as model:
        mean = pm.Normal('mean', mu=0.0, sd=10.0)
        u_star = xo.QuadLimbDark('u_star')
        M_star_huang = 1.094, 0.039
        R_star_huang = 1.1, 0.023
        BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=3)
        m_star = BoundedNormal('m_star', mu=M_star_huang[0], sd=M_star_huang[1]
            )
        r_star = BoundedNormal('r_star', mu=R_star_huang[0], sd=R_star_huang[1]
            )
        period = pm.Lognormal('period', mu=np.log(bls_period), sd=1)
        t0 = pm.Normal('t0', mu=bls_t0, sd=1)
        r_pl = pm.Lognormal('r_pl', sd=1.0, mu=0.5 * np.log(0.001 * np.
            array(bls_depth)) + np.log(R_star_huang[0]))
        ror = pm.Deterministic('ror', r_pl / r_star)
        b = xo.distributions.ImpactParameter('b', ror=ror)
        ecs = pmx.UnitDisk('ecs', testval=np.array([0.01, 0.0]))
        ecc = pm.Deterministic('ecc', tt.sum(ecs ** 2))
        omega = pm.Deterministic('omega', tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.kipping13('ecc_prior', fixed=True, observed=ecc)
        sigma_lc = pm.Lognormal('sigma_lc', mu=np.log(np.std(y[mask])), sd=10)
        rho_gp = pm.Lognormal('rho_gp', mu=0, sd=10)
        sigma_gp = pm.Lognormal('sigma_gp', mu=np.log(np.std(y[mask])), sd=10)
        orbit = xo.orbits.KeplerianOrbit(r_star=r_star, m_star=m_star,
            period=period, t0=t0, b=b, ecc=ecc, omega=omega)
        light_curves = pm.Deterministic('light_curves', xo.
            LimbDarkLightCurve(u_star).get_light_curve(orbit=orbit, r=r_pl,
            t=x[mask], texp=texp) * 1000.0)
        light_curve = tt.sum(light_curves, axis=-1) + mean
        resid = y[mask] - light_curve
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1 / np.sqrt(2))
        gp = GaussianProcess(kernel, t=x[mask], yerr=sigma_lc)
        gp.marginal('gp', observed=resid)
        pm.Deterministic('gp_pred', gp.predict(resid))
        if start is None:
            start = model.test_point
        map_soln = pmx.optimize(start=start, vars=[sigma_lc, sigma_gp, rho_gp])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[period, t0])
        map_soln = pmx.optimize(start=map_soln, vars=[u_star])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[ecs])
        map_soln = pmx.optimize(start=map_soln, vars=[mean])
        map_soln = pmx.optimize(start=map_soln, vars=[sigma_lc, sigma_gp,
            rho_gp])
        map_soln = pmx.optimize(start=map_soln)
    return model, map_soln


model0, map_soln0 = build_model()


def plot_light_curve(soln, mask=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
    ax = axes[0]
    ax.plot(x[mask], y[mask], 'k', label='data')
    gp_mod = soln['gp_pred'] + soln['mean']
    ax.plot(x[mask], gp_mod, color='C2', label='gp model')
    ax.legend(fontsize=10)
    ax.set_ylabel('relative flux [ppt]')
    ax = axes[1]
    ax.plot(x[mask], y[mask] - gp_mod, 'k', label='de-trended data')
    for i, l in enumerate('b'):
        mod = soln['light_curves'][:, i]
        ax.plot(x[mask], mod, label='planet {0}'.format(l))
    ax.legend(fontsize=10, loc=3)
    ax.set_ylabel('de-trended flux [ppt]')
    ax = axes[2]
    mod = gp_mod + np.sum(soln['light_curves'], axis=-1)
    ax.plot(x[mask], y[mask] - mod, 'k')
    ax.axhline(0, color='#aaaaaa', lw=1)
    ax.set_ylabel('residuals [ppt]')
    ax.set_xlim(x[mask].min(), x[mask].max())
    ax.set_xlabel('time [days]')
    return fig


_ = plot_light_curve(map_soln0)
mod = map_soln0['gp_pred'] + map_soln0['mean'] + np.sum(map_soln0[
    'light_curves'], axis=-1)
resid = y - mod
rms = np.sqrt(np.median(resid ** 2))
mask = np.abs(resid) < 5 * rms
plt.figure(figsize=(10, 5))
plt.plot(x, resid, 'k', label='data')
plt.plot(x[~mask], resid[~mask], 'xr', label='outliers')
plt.axhline(0, color='#aaaaaa', lw=1)
plt.ylabel('residuals [ppt]')
plt.xlabel('time [days]')
plt.legend(fontsize=12, loc=3)
_ = plt.xlim(x.min(), x.max())
model, map_soln = build_model(mask, map_soln0)
_ = plot_light_curve(map_soln, mask)
np.random.seed(261136679)
with model:
    trace = pmx.sample(tune=2500, draws=2000, start=map_soln, cores=2,
        chains=2, initial_accept=0.8, target_accept=0.95)
with model:
    summary = pm.summary(trace, var_names=['omega', 'ecc', 'r_pl', 'b',
        't0', 'period', 'r_star', 'm_star', 'u_star', 'mean'])
summary
gp_mod = np.median(trace['gp_pred'] + trace['mean'][:, None], axis=0)
p = np.median(trace['period'])
t0 = np.median(trace['t0'])
x_fold = (x[mask] - t0 + 0.5 * p) % p - 0.5 * p
plt.plot(x_fold, y[mask] - gp_mod, '.k', label='data', zorder=-1000)
bins = np.linspace(-0.41, 0.41, 50)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y[mask])
denom[num == 0] = 1.0
plt.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, 'o', color='C2', label=
    'binned')
inds = np.argsort(x_fold)
inds = inds[np.abs(x_fold)[inds] < 0.3]
pred = trace['light_curves'][:, inds, 0]
pred = np.percentile(pred, [16, 50, 84], axis=0)
plt.plot(x_fold[inds], pred[1], color='C1', label='model')
art = plt.fill_between(x_fold[inds], pred[0], pred[2], color='C1', alpha=
    0.5, zorder=1000)
art.set_edgecolor('none')
txt = 'period = {0:.5f} +/- {1:.5f} d'.format(np.mean(trace['period']), np.
    std(trace['period']))
plt.annotate(txt, (0, 0), xycoords='axes fraction', xytext=(5, 5),
    textcoords='offset points', ha='left', va='bottom', fontsize=12)
plt.legend(fontsize=10, loc=4)
plt.xlim(-0.5 * p, 0.5 * p)
plt.xlabel('time since transit [days]')
plt.ylabel('de-trended flux')
_ = plt.xlim(-0.15, 0.15)
<mask token>
varnames = ['period', 'b', 'ecc', 'r_pl']
samples = pm.trace_to_dataframe(trace, varnames=varnames)
samples['r_pl'] = (np.array(samples['r_pl']) * u.R_sun).to(u.R_earth).value
_ = corner.corner(samples[['period', 'r_pl', 'b', 'ecc']], labels=[
    'period [days]', 'radius [Earth radii]', 'impact param', 'eccentricity'])
with model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)
print('\n'.join(bib.splitlines()[:10]) + '\n...')

import lightkurve as lk
import numpy as np
import lightkurve as lk
import matplotlib.pyplot as plt
lc_file = lk.search_lightcurvefile('TIC 261136679', sector=1).download(
    quality_bitmask='hardest')
lc = lc_file.PDCSAP_FLUX.remove_nans().normalize().remove_outliers()
time = lc.time
flux = lc.flux
m = lc.quality == 0
with lc_file.hdu as hdu:
    hdr = hdu[1].header
texp = hdr['FRAMETIM'] * hdr['NUM_FRM']
texp /= 60.0 * 60.0 * 24.0
ref_time = 0.5 * (np.min(time) + np.max(time))
x = np.ascontiguousarray(time[m] - ref_time, dtype=np.float64)
y = np.ascontiguousarray(1000.0 * (flux[m] - 1.0), dtype=np.float64)
plt.plot(x, y, '.k')
plt.xlabel('time [days]')
plt.ylabel('relative flux [ppt]')
_ = plt.xlim(x.min(), x.max())
from astropy.timeseries import BoxLeastSquares
period_grid = np.exp(np.linspace(np.log(1), np.log(15), 50000))
bls = BoxLeastSquares(x, y)
bls_power = bls.power(period_grid, 0.1, oversample=20)
index = np.argmax(bls_power.power)
bls_period = bls_power.period[index]
bls_t0 = bls_power.transit_time[index]
bls_depth = bls_power.depth[index]
transit_mask = bls.transit_mask(x, bls_period, 0.2, bls_t0)
fig, axes = plt.subplots(2, 1, figsize=(10, 10))
ax = axes[0]
ax.axvline(np.log10(bls_period), color='C1', lw=5, alpha=0.8)
ax.plot(np.log10(bls_power.period), bls_power.power, 'k')
ax.annotate('period = {0:.4f} d'.format(bls_period), (0, 1), xycoords=
    'axes fraction', xytext=(5, -5), textcoords='offset points', va='top',
    ha='left', fontsize=12)
ax.set_ylabel('bls power')
ax.set_yticks([])
ax.set_xlim(np.log10(period_grid.min()), np.log10(period_grid.max()))
ax.set_xlabel('log10(period)')
ax = axes[1]
x_fold = (x - bls_t0 + 0.5 * bls_period) % bls_period - 0.5 * bls_period
m = np.abs(x_fold) < 0.4
ax.plot(x_fold[m], y[m], '.k')
bins = np.linspace(-0.41, 0.41, 32)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y)
denom[num == 0] = 1.0
ax.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, color='C1')
ax.set_xlim(-0.3, 0.3)
ax.set_ylabel('de-trended flux [ppt]')
_ = ax.set_xlabel('time since transit')
import exoplanet as xo
import pymc3 as pm
import theano.tensor as tt
import pymc3_ext as pmx
from celerite2.theano import terms, GaussianProcess


def build_model(mask=None, start=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    with pm.Model() as model:
        mean = pm.Normal('mean', mu=0.0, sd=10.0)
        u_star = xo.QuadLimbDark('u_star')
        M_star_huang = 1.094, 0.039
        R_star_huang = 1.1, 0.023
        BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=3)
        m_star = BoundedNormal('m_star', mu=M_star_huang[0], sd=M_star_huang[1]
            )
        r_star = BoundedNormal('r_star', mu=R_star_huang[0], sd=R_star_huang[1]
            )
        period = pm.Lognormal('period', mu=np.log(bls_period), sd=1)
        t0 = pm.Normal('t0', mu=bls_t0, sd=1)
        r_pl = pm.Lognormal('r_pl', sd=1.0, mu=0.5 * np.log(0.001 * np.
            array(bls_depth)) + np.log(R_star_huang[0]))
        ror = pm.Deterministic('ror', r_pl / r_star)
        b = xo.distributions.ImpactParameter('b', ror=ror)
        ecs = pmx.UnitDisk('ecs', testval=np.array([0.01, 0.0]))
        ecc = pm.Deterministic('ecc', tt.sum(ecs ** 2))
        omega = pm.Deterministic('omega', tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.kipping13('ecc_prior', fixed=True, observed=ecc)
        sigma_lc = pm.Lognormal('sigma_lc', mu=np.log(np.std(y[mask])), sd=10)
        rho_gp = pm.Lognormal('rho_gp', mu=0, sd=10)
        sigma_gp = pm.Lognormal('sigma_gp', mu=np.log(np.std(y[mask])), sd=10)
        orbit = xo.orbits.KeplerianOrbit(r_star=r_star, m_star=m_star,
            period=period, t0=t0, b=b, ecc=ecc, omega=omega)
        light_curves = pm.Deterministic('light_curves', xo.
            LimbDarkLightCurve(u_star).get_light_curve(orbit=orbit, r=r_pl,
            t=x[mask], texp=texp) * 1000.0)
        light_curve = tt.sum(light_curves, axis=-1) + mean
        resid = y[mask] - light_curve
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1 / np.sqrt(2))
        gp = GaussianProcess(kernel, t=x[mask], yerr=sigma_lc)
        gp.marginal('gp', observed=resid)
        pm.Deterministic('gp_pred', gp.predict(resid))
        if start is None:
            start = model.test_point
        map_soln = pmx.optimize(start=start, vars=[sigma_lc, sigma_gp, rho_gp])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[period, t0])
        map_soln = pmx.optimize(start=map_soln, vars=[u_star])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[ecs])
        map_soln = pmx.optimize(start=map_soln, vars=[mean])
        map_soln = pmx.optimize(start=map_soln, vars=[sigma_lc, sigma_gp,
            rho_gp])
        map_soln = pmx.optimize(start=map_soln)
    return model, map_soln


model0, map_soln0 = build_model()


def plot_light_curve(soln, mask=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
    ax = axes[0]
    ax.plot(x[mask], y[mask], 'k', label='data')
    gp_mod = soln['gp_pred'] + soln['mean']
    ax.plot(x[mask], gp_mod, color='C2', label='gp model')
    ax.legend(fontsize=10)
    ax.set_ylabel('relative flux [ppt]')
    ax = axes[1]
    ax.plot(x[mask], y[mask] - gp_mod, 'k', label='de-trended data')
    for i, l in enumerate('b'):
        mod = soln['light_curves'][:, i]
        ax.plot(x[mask], mod, label='planet {0}'.format(l))
    ax.legend(fontsize=10, loc=3)
    ax.set_ylabel('de-trended flux [ppt]')
    ax = axes[2]
    mod = gp_mod + np.sum(soln['light_curves'], axis=-1)
    ax.plot(x[mask], y[mask] - mod, 'k')
    ax.axhline(0, color='#aaaaaa', lw=1)
    ax.set_ylabel('residuals [ppt]')
    ax.set_xlim(x[mask].min(), x[mask].max())
    ax.set_xlabel('time [days]')
    return fig


_ = plot_light_curve(map_soln0)
mod = map_soln0['gp_pred'] + map_soln0['mean'] + np.sum(map_soln0[
    'light_curves'], axis=-1)
resid = y - mod
rms = np.sqrt(np.median(resid ** 2))
mask = np.abs(resid) < 5 * rms
plt.figure(figsize=(10, 5))
plt.plot(x, resid, 'k', label='data')
plt.plot(x[~mask], resid[~mask], 'xr', label='outliers')
plt.axhline(0, color='#aaaaaa', lw=1)
plt.ylabel('residuals [ppt]')
plt.xlabel('time [days]')
plt.legend(fontsize=12, loc=3)
_ = plt.xlim(x.min(), x.max())
model, map_soln = build_model(mask, map_soln0)
_ = plot_light_curve(map_soln, mask)
np.random.seed(261136679)
with model:
    trace = pmx.sample(tune=2500, draws=2000, start=map_soln, cores=2,
        chains=2, initial_accept=0.8, target_accept=0.95)
with model:
    summary = pm.summary(trace, var_names=['omega', 'ecc', 'r_pl', 'b',
        't0', 'period', 'r_star', 'm_star', 'u_star', 'mean'])
summary
gp_mod = np.median(trace['gp_pred'] + trace['mean'][:, None], axis=0)
p = np.median(trace['period'])
t0 = np.median(trace['t0'])
x_fold = (x[mask] - t0 + 0.5 * p) % p - 0.5 * p
plt.plot(x_fold, y[mask] - gp_mod, '.k', label='data', zorder=-1000)
bins = np.linspace(-0.41, 0.41, 50)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y[mask])
denom[num == 0] = 1.0
plt.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, 'o', color='C2', label=
    'binned')
inds = np.argsort(x_fold)
inds = inds[np.abs(x_fold)[inds] < 0.3]
pred = trace['light_curves'][:, inds, 0]
pred = np.percentile(pred, [16, 50, 84], axis=0)
plt.plot(x_fold[inds], pred[1], color='C1', label='model')
art = plt.fill_between(x_fold[inds], pred[0], pred[2], color='C1', alpha=
    0.5, zorder=1000)
art.set_edgecolor('none')
txt = 'period = {0:.5f} +/- {1:.5f} d'.format(np.mean(trace['period']), np.
    std(trace['period']))
plt.annotate(txt, (0, 0), xycoords='axes fraction', xytext=(5, 5),
    textcoords='offset points', ha='left', va='bottom', fontsize=12)
plt.legend(fontsize=10, loc=4)
plt.xlim(-0.5 * p, 0.5 * p)
plt.xlabel('time since transit [days]')
plt.ylabel('de-trended flux')
_ = plt.xlim(-0.15, 0.15)
import corner
import astropy.units as u
varnames = ['period', 'b', 'ecc', 'r_pl']
samples = pm.trace_to_dataframe(trace, varnames=varnames)
samples['r_pl'] = (np.array(samples['r_pl']) * u.R_sun).to(u.R_earth).value
_ = corner.corner(samples[['period', 'r_pl', 'b', 'ecc']], labels=[
    'period [days]', 'radius [Earth radii]', 'impact param', 'eccentricity'])
with model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)
print('\n'.join(bib.splitlines()[:10]) + '\n...')

# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.6.0
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + nbsphinx="hidden"
import lightkurve as lk

# %matplotlib inline

# + nbsphinx="hidden"
# %run notebook_setup
# -

# # Fitting TESS data

# In this tutorial, we will reproduce the fits to the transiting planet in the Pi Mensae system discovered by [Huang et al. (2018)](https://arxiv.org/abs/1809.05967).
# The data processing and model are similar to the [Putting it all together](./together.ipynb) case study, but with a few extra bits like aperture selection and de-trending.
#
# To start, we need to download the target pixel file:

# +
import numpy as np
import lightkurve as lk
import matplotlib.pyplot as plt

lc_file = lk.search_lightcurvefile("TIC 261136679", sector=1).download(
    quality_bitmask="hardest"
)
lc = lc_file.PDCSAP_FLUX.remove_nans().normalize().remove_outliers()
time = lc.time
flux = lc.flux
m = lc.quality == 0
with lc_file.hdu as hdu:
    hdr = hdu[1].header

texp = hdr["FRAMETIM"] * hdr["NUM_FRM"]
texp /= 60.0 * 60.0 * 24.0

ref_time = 0.5 * (np.min(time) + np.max(time))
x = np.ascontiguousarray(time[m] - ref_time, dtype=np.float64)
y = np.ascontiguousarray(1e3 * (flux[m] - 1.0), dtype=np.float64)

plt.plot(x, y, ".k")
plt.xlabel("time [days]")
plt.ylabel("relative flux [ppt]")
_ = plt.xlim(x.min(), x.max())
# -

# ## Transit search
#
# Now, let's use [the box least squares periodogram from AstroPy](http://docs.astropy.org/en/latest/timeseries/bls.html)
# (Note: you'll need AstroPy v3.1 or more recent to use this feature) to estimate the period, phase, and depth of the transit.

# +
from astropy.timeseries import BoxLeastSquares

period_grid = np.exp(np.linspace(np.log(1), np.log(15), 50000))

bls = BoxLeastSquares(x, y)
bls_power = bls.power(period_grid, 0.1, oversample=20)

# Save the highest peak as the planet candidate
index = np.argmax(bls_power.power)
bls_period = bls_power.period[index]
bls_t0 = bls_power.transit_time[index]
bls_depth = bls_power.depth[index]
transit_mask = bls.transit_mask(x, bls_period, 0.2, bls_t0)

fig, axes = plt.subplots(2, 1, figsize=(10, 10))

# Plot the periodogram
ax = axes[0]
ax.axvline(np.log10(bls_period), color="C1", lw=5, alpha=0.8)
ax.plot(np.log10(bls_power.period), bls_power.power, "k")
ax.annotate(
    "period = {0:.4f} d".format(bls_period),
    (0, 1),
    xycoords="axes fraction",
    xytext=(5, -5),
    textcoords="offset points",
    va="top",
    ha="left",
    fontsize=12,
)
ax.set_ylabel("bls power")
ax.set_yticks([])
ax.set_xlim(np.log10(period_grid.min()), np.log10(period_grid.max()))
ax.set_xlabel("log10(period)")

# Plot the folded transit
ax = axes[1]
x_fold = (x - bls_t0 + 0.5 * bls_period) % bls_period - 0.5 * bls_period
m = np.abs(x_fold) < 0.4
ax.plot(x_fold[m], y[m], ".k")

# Overplot the phase binned light curve
bins = np.linspace(-0.41, 0.41, 32)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y)
denom[num == 0] = 1.0
ax.plot(0.5 * (bins[1:] + bins[:-1]), num / denom, color="C1")

ax.set_xlim(-0.3, 0.3)
ax.set_ylabel("de-trended flux [ppt]")
_ = ax.set_xlabel("time since transit")
# -

# ## The transit model in PyMC3
#
# The transit model, initialization, and sampling are all nearly the same as the one in [Putting it all together](./together.ipynb).

# +
import exoplanet as xo
import pymc3 as pm
import theano.tensor as tt

import pymc3_ext as pmx
from celerite2.theano import terms, GaussianProcess


def build_model(mask=None, start=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    with pm.Model() as model:

        # Parameters for the stellar properties
        mean = pm.Normal("mean", mu=0.0, sd=10.0)
        u_star = xo.QuadLimbDark("u_star")

        # Stellar parameters from Huang et al (2018)
        M_star_huang = 1.094, 0.039
        R_star_huang = 1.10, 0.023
        BoundedNormal = pm.Bound(pm.Normal, lower=0, upper=3)
        m_star = BoundedNormal(
            "m_star", mu=M_star_huang[0], sd=M_star_huang[1]
        )
        r_star = BoundedNormal(
            "r_star", mu=R_star_huang[0], sd=R_star_huang[1]
        )

        # Orbital parameters for the planets
        period = pm.Lognormal("period", mu=np.log(bls_period), sd=1)
        t0 = pm.Normal("t0", mu=bls_t0, sd=1)
        r_pl = pm.Lognormal(
            "r_pl",
            sd=1.0,
            mu=0.5 * np.log(1e-3 * np.array(bls_depth))
            + np.log(R_star_huang[0]),
        )
        ror = pm.Deterministic("ror", r_pl / r_star)
        b = xo.distributions.ImpactParameter("b", ror=ror)

        ecs = pmx.UnitDisk("ecs", testval=np.array([0.01, 0.0]))
        ecc = pm.Deterministic("ecc", tt.sum(ecs ** 2))
        omega = pm.Deterministic("omega", tt.arctan2(ecs[1], ecs[0]))
        xo.eccentricity.kipping13("ecc_prior", fixed=True, observed=ecc)

        # Transit jitter & GP parameters
        sigma_lc = pm.Lognormal("sigma_lc", mu=np.log(np.std(y[mask])), sd=10)
        rho_gp = pm.Lognormal("rho_gp", mu=0, sd=10)
        sigma_gp = pm.Lognormal("sigma_gp", mu=np.log(np.std(y[mask])), sd=10)

        # Orbit model
        orbit = xo.orbits.KeplerianOrbit(
            r_star=r_star,
            m_star=m_star,
            period=period,
            t0=t0,
            b=b,
            ecc=ecc,
            omega=omega,
        )

        # Compute the model light curve
        light_curves = pm.Deterministic(
            "light_curves",
            xo.LimbDarkLightCurve(u_star).get_light_curve(
                orbit=orbit, r=r_pl, t=x[mask], texp=texp
            )
            * 1e3,
        )
        light_curve = tt.sum(light_curves, axis=-1) + mean
        resid = y[mask] - light_curve

        # GP model for the light curve
        kernel = terms.SHOTerm(sigma=sigma_gp, rho=rho_gp, Q=1 / np.sqrt(2))
        gp = GaussianProcess(kernel, t=x[mask], yerr=sigma_lc)
        gp.marginal("gp", observed=resid)
        pm.Deterministic("gp_pred", gp.predict(resid))

        # Fit for the maximum a posteriori parameters, I've found that I can get
        # a better solution by trying different combinations of parameters in turn
        if start is None:
            start = model.test_point
        map_soln = pmx.optimize(start=start, vars=[sigma_lc, sigma_gp, rho_gp])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[period, t0])
        map_soln = pmx.optimize(start=map_soln, vars=[u_star])
        map_soln = pmx.optimize(start=map_soln, vars=[r_pl])
        map_soln = pmx.optimize(start=map_soln, vars=[b])
        map_soln = pmx.optimize(start=map_soln, vars=[ecs])
        map_soln = pmx.optimize(start=map_soln, vars=[mean])
        map_soln = pmx.optimize(
            start=map_soln, vars=[sigma_lc, sigma_gp, rho_gp]
        )
        map_soln = pmx.optimize(start=map_soln)

    return model, map_soln


model0, map_soln0 = build_model()


# -

# Here's how we plot the initial light curve model:

# +
def plot_light_curve(soln, mask=None):
    if mask is None:
        mask = np.ones(len(x), dtype=bool)

    fig, axes = plt.subplots(3, 1, figsize=(10, 7), sharex=True)

    ax = axes[0]
    ax.plot(x[mask], y[mask], "k", label="data")
    gp_mod = soln["gp_pred"] + soln["mean"]
    ax.plot(x[mask], gp_mod, color="C2", label="gp model")
    ax.legend(fontsize=10)
    ax.set_ylabel("relative flux [ppt]")

    ax = axes[1]
    ax.plot(x[mask], y[mask] - gp_mod, "k", label="de-trended data")
    for i, l in enumerate("b"):
        mod = soln["light_curves"][:, i]
        ax.plot(x[mask], mod, label="planet {0}".format(l))
    ax.legend(fontsize=10, loc=3)
    ax.set_ylabel("de-trended flux [ppt]")

    ax = axes[2]
    mod = gp_mod + np.sum(soln["light_curves"], axis=-1)
    ax.plot(x[mask], y[mask] - mod, "k")
    ax.axhline(0, color="#aaaaaa", lw=1)
    ax.set_ylabel("residuals [ppt]")
    ax.set_xlim(x[mask].min(), x[mask].max())
    ax.set_xlabel("time [days]")

    return fig


_ = plot_light_curve(map_soln0)
# -

# As in [Putting it all together](./together.ipynb), we can do some sigma clipping to remove significant outliers.

# +
mod = (
    map_soln0["gp_pred"]
    + map_soln0["mean"]
    + np.sum(map_soln0["light_curves"], axis=-1)
)
resid = y - mod
rms = np.sqrt(np.median(resid ** 2))
mask = np.abs(resid) < 5 * rms

plt.figure(figsize=(10, 5))
plt.plot(x, resid, "k", label="data")
plt.plot(x[~mask], resid[~mask], "xr", label="outliers")
plt.axhline(0, color="#aaaaaa", lw=1)
plt.ylabel("residuals [ppt]")
plt.xlabel("time [days]")
plt.legend(fontsize=12, loc=3)
_ = plt.xlim(x.min(), x.max())
# -

# And then we re-build the model using the data without outliers.

model, map_soln = build_model(mask, map_soln0)
_ = plot_light_curve(map_soln, mask)

# Now that we have the model, we can sample:

np.random.seed(261136679)
with model:
    trace = pmx.sample(
        tune=2500,
        draws=2000,
        start=map_soln,
        cores=2,
        chains=2,
        initial_accept=0.8,
        target_accept=0.95,
    )

with model:
    summary = pm.summary(
        trace,
        var_names=[
            "omega",
            "ecc",
            "r_pl",
            "b",
            "t0",
            "period",
            "r_star",
            "m_star",
            "u_star",
            "mean",
        ],
    )
summary

# ## Results
#
# After sampling, we can make the usual plots.
# First, let's look at the folded light curve plot:

# +
# Compute the GP prediction
gp_mod = np.median(trace["gp_pred"] + trace["mean"][:, None], axis=0)

# Get the posterior median orbital parameters
p = np.median(trace["period"])
t0 = np.median(trace["t0"])

# Plot the folded data
x_fold = (x[mask] - t0 + 0.5 * p) % p - 0.5 * p
plt.plot(x_fold, y[mask] - gp_mod, ".k", label="data", zorder=-1000)

# Overplot the phase binned light curve
bins = np.linspace(-0.41, 0.41, 50)
denom, _ = np.histogram(x_fold, bins)
num, _ = np.histogram(x_fold, bins, weights=y[mask])
denom[num == 0] = 1.0
plt.plot(
    0.5 * (bins[1:] + bins[:-1]), num / denom, "o", color="C2", label="binned"
)

# Plot the folded model
inds = np.argsort(x_fold)
inds = inds[np.abs(x_fold)[inds] < 0.3]
pred = trace["light_curves"][:, inds, 0]
pred = np.percentile(pred, [16, 50, 84], axis=0)
plt.plot(x_fold[inds], pred[1], color="C1", label="model")
art = plt.fill_between(
    x_fold[inds], pred[0], pred[2], color="C1", alpha=0.5, zorder=1000
)
art.set_edgecolor("none")

# Annotate the plot with the planet's period
txt = "period = {0:.5f} +/- {1:.5f} d".format(
    np.mean(trace["period"]), np.std(trace["period"])
)
plt.annotate(
    txt,
    (0, 0),
    xycoords="axes fraction",
    xytext=(5, 5),
    textcoords="offset points",
    ha="left",
    va="bottom",
    fontsize=12,
)

plt.legend(fontsize=10, loc=4)
plt.xlim(-0.5 * p, 0.5 * p)
plt.xlabel("time since transit [days]")
plt.ylabel("de-trended flux")
_ = plt.xlim(-0.15, 0.15)
# -

# And a corner plot of some of the key parameters:

# +
import corner
import astropy.units as u

varnames = ["period", "b", "ecc", "r_pl"]
samples = pm.trace_to_dataframe(trace, varnames=varnames)

# Convert the radius to Earth radii
samples["r_pl"] = (np.array(samples["r_pl"]) * u.R_sun).to(u.R_earth).value

_ = corner.corner(
    samples[["period", "r_pl", "b", "ecc"]],
    labels=[
        "period [days]",
        "radius [Earth radii]",
        "impact param",
        "eccentricity",
    ],
)
# -

# These all seem consistent with the previously published values.

# ## Citations
#
# As described in the [citation tutorial](https://docs.exoplanet.codes/en/stable/tutorials/citation/), we can use [citations.get_citations_for_model](https://docs.exoplanet.codes/en/stable/user/api/#exoplanet.citations.get_citations_for_model) to construct an acknowledgement and BibTeX listing that includes the relevant citations for this model.

with model:
    txt, bib = xo.citations.get_citations_for_model()
print(txt)

print("\n".join(bib.splitlines()[:10]) + "\n...")
