TESS Atlas fit for TOI 349

Version: 0.1.1

Note: This notebook was automatically generated as part of the TESS Atlas project. More information can be found on GitHub: github.com/dfm/tess-atlas

In this notebook, we do a quicklook fit for the parameters of the TESS Objects of Interest (TOI) in the system number 349. To do this fit, we use the exoplanet library and you can find more information about that project at exoplanet.dfm.io.

From here, you can scroll down and take a look at the fit results, or you can:


There are many caveats associated with this relatively simple "quicklook" type of analysis that should be kept in mind. Here are some of the main things that come to mind:

  1. The orbits that we fit are constrained to be circular. One major effect of this approximation is that the fit will significantly overestimate the confidence of the impact parameter constraint, so the results for impact parameter shouldn't be taken too seriously.

  2. Transit timing variations, correlated noise, and (probably) your favorite systematics are ignored. Sorry!

  3. This notebook was generated automatically without human intervention. Use at your own risk!

Table of Contents

  1. Getting started
  2. Data & de-trending
  3. Removing stellar variability
  4. Transit model in PyMC3 & exoplanet
  5. Sampling
  6. Posterior constraints
  7. Attribution

Getting started

To get going, we'll need to make out plots show up inline and install a few packages:

In [1]:
%matplotlib inline
!pip install -q -U lightkurve fbpca exoplanet corner pymc3

Then we'll set up the plotting styles and do all of the imports:

In [2]:
import matplotlib.pyplot as plt

from matplotlib import rcParams
rcParams["savefig.dpi"] = 100
rcParams["figure.dpi"] = 100
rcParams["font.size"] = 16
rcParams["text.usetex"] = False
rcParams["font.family"] = ["sans-serif"]
rcParams["font.sans-serif"] = ["cmss10"]
rcParams["axes.unicode_minus"] = False

import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)
warnings.filterwarnings("ignore", category=FutureWarning)
import logging
logger = logging.getLogger("theano.gof.compilelock")

import corner
import numpy as np
import pandas as pd
import lightkurve as lk
import matplotlib.pyplot as plt

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

Next, we grab the TOI list from ExoFOP to get the information about the system:

In [3]:
toi_num = 349

# Get the table of TOI info from ExoFOP
tois = pd.read_csv("https://exofop.ipac.caltech.edu/tess/download_toi.php?sort=toi&output=csv")

# Select all of the rows in the TOI table that are associated with this target
toi = tois[tois["TOI"] == toi_num + 0.01].iloc[0]
tic = toi['TIC ID']
tois = tois[tois["TIC ID"] == tic].sort_values("TOI")

# Extract the planet periods
periods = np.array(tois["Period (days)"], dtype=float)

# Convert the phase to TBJD from BJD
t0s = np.array(tois["Epoch (BJD)"], dtype=float) - 2457000

# Convert the depth to parts per thousand from parts per million
depths = 1e-3 * np.array(tois["Depth (ppm)"], dtype=float)

# Convert the duration to days from hours
durations = np.array(tois["Duration (hours)"], dtype=float) / 24.0

# Extract the stellar radius from the table
toi_r_star = toi['Stellar Radius (R_Sun)']
toi_r_star_err = toi['Stellar Radius (R_Sun) err']
toi_logg_star = toi['Stellar log(g) (cm/s^2)']
toi_logg_star_err = toi['Stellar log(g) (cm/s^2) err']

# If there is no entry in the table (does this ever happen?)
if not (np.isfinite(toi_r_star) and np.isfinite(toi_r_star_err)):
    raise ValueError("no radius measurement in table")
if not (np.isfinite(toi_logg_star) and np.isfinite(toi_logg_star_err)):
    raise ValueError("no logg measurement in table")

# These are the letters that will be used to identify each candidate
# (are we being a bit optimistic?)
letters = "bcdefghijklmnopqrstuvwxyz"[:len(periods)]

Then we use the lightkurve library to download and de-trend the time series using pixel-level decorrelation (PLD). We read in target pixel files (TPFs) for each of the campaigns in which TOI 349 was observed. To remove systematic noise, we mask out known transits and perform second order PLD. The noise-corrected light curves are stitched together to create a single contiguous light curve.

In [4]:
# Download the target pixel files
sr = lk.search_targetpixelfile('TIC %i' % tic)
tpf_collection = sr.download_all(quality_bitmask="hardest")
if tpf_collection is None or not len(tpf_collection):
    raise ValueError("the TESS atlas only works for TPF files")
# Extract the exposure time associated with the TPF
hdr = tpf_collection[0].hdu[1].header
texp = hdr["FRAMETIM"] * hdr["NUM_FRM"]
texp /= 60.0 * 60.0 * 24.0

# This function can be used to estimate which data points are in transit
# for known phase, period, and duration
def get_transit_mask(t, t0, period, duration):
    hp = 0.5*period
    return np.abs((t-t0+hp) % period - hp) < 0.5*duration

# Run PLD on each TPF to extract the light curves
lc_collection = []
for tpf in tpf_collection:
    mask = np.zeros_like(tpf.time, dtype=bool)
    for i in range(len(periods)):
        mask |= get_transit_mask(tpf.time, t0s[i], periods[i], 5*durations[i])
    pld = tpf.to_corrector("pld")
        lc = pld.correct(aperture_mask="pipeline", cadence_mask=~mask, use_gp=False, pld_order=2)
    except ValueError:
        lc = pld.correct(aperture_mask="pipeline", cadence_mask=~mask, use_gp=False, pld_order=1)

# Normalize and stitch the sectors
lc = lc_collection[0]
if len(lc_collection) > 1:
    lc = lc.append([next_lc for next_lc in lc_collection[1:]])
# Remove outliers
_, outliers = lc.remove_outliers(return_mask=True)
mask = np.zeros_like(lc.time, dtype=bool)
for i in range(len(periods)):
    mask |= get_transit_mask(lc.time, t0s[i], periods[i], 2*durations[i])
outliers[mask] = False
lc = lc[~outliers]

Removing stellar variability

Next up, we remove stellar variability using a Gaussian Processes model fit to the out of transit data.

In [5]:
# Extract the data and convert to parts per thousand
x = np.ascontiguousarray(lc.time, dtype=np.float64)
y = np.ascontiguousarray((lc.flux - 1.0) * 1e3, dtype=np.float64)
yerr = np.ascontiguousarray(lc.flux_err * 1e3, dtype=np.float64)

# Compute the transit mask
mask = np.zeros_like(x, dtype=bool)
for i in range(len(periods)):
    mask |= get_transit_mask(x, t0s[i], periods[i], 5*durations[i])

# Temporarily increase the in transit error bars substantially
diag = np.array(yerr**2)
diag[mask] += 10000.0

# Build a GP model
with pm.Model() as model:
    logs2 = pm.Normal("logs2", mu=np.log(1e-4*np.var(y)), sd=10)
    logsigma = pm.Normal("logsigma", mu=np.log(np.std(y)), sd=10)
    logrho = pm.Normal("logrho", mu=np.log(10.0), sd=10.0)
    kernel = xo.gp.terms.Matern32Term(log_sigma=logsigma, log_rho=logrho)
    gp = xo.gp.GP(kernel, x, diag + tt.exp(logs2), J=2)
    pm.Potential("loglike", gp.log_likelihood(y))
    map_soln = xo.optimize()
    pred = xo.utils.eval_in_model(gp.predict(), map_soln)

# Flatten the light curve
y -= pred
optimizing logp for variables: ['logrho', 'logsigma', 'logs2']
message: Optimization terminated successfully.
logp: -51303.11807668808 -> -50001.299705201694

Transit model in PyMC3 & exoplanet

Here's how we set up the transit model using exoplanet and PyMC3. For more information about how to use these libraries take a look at the docs that are linked above. In this model, the parameters that we're fitting are:

  • mean: the mean (out-of-transit) flux of the star,
  • r_star: the radius of the star (with the prior from the TOI list),
  • logg_star: the surface gravity of the star (with the prior from the TOI list),
  • u: the quadratic limb darkening parameters, parameterized following Kipping (2013)
  • t0: the time of a reference transit for each planet,
  • logP: the log of the obribatl periods,
  • r: the planet radius ratios (relative to the star),
  • b: the impact parameter in units of the stellar radius, b and r are both parameterized following Espinoza (2018), and
  • logs2: a jitter parameter that captures excess noise or underrestimated error bars.

A few key assumptions include:

  • The orbits are assumed to be circular so the constraints on impact parameter (which would be severely degenerate with eccentricity) will be tighter than they should be.
  • The noise is assumed to be Gaussian and independent. This means that all correlated noise should be removed in advance. Since we flattened the light curve using a Gaussian process above, this should be not totally unreasonable.
  • We are neglecting transit times (the ephemeris is assumed to be linear) which should be sufficient for most cases with the short TESS baseline, but transit timing variations could be important for some targets.
In [6]:
# factor * 10**logg / r_star = rho
factor = 5.141596357654149e-05

def build_model(x, y, yerr, periods, t0s, depths, mask=None, start=None):
    """Build an exoplanet model for a dataset and set of planets
        x: The time series (in days); this should probably be centered
        y: The relative fluxes (in parts per thousand)
        yerr: The uncertainties on ``y``
        periods: The periods of the planets (in days)
        t0s: The phases of the planets in the same coordinates as ``x``
        depths: The depths of the transits in parts per thousand
        mask: A boolean mask with the same shape as ``x`` indicating which
            data points should be included in the fit
        start: A dictionary of model parameters where the optimization
            should be initialized
        A PyMC3 model specifying the probabilistic model for the light curve

    if mask is None:
        mask = np.ones(len(x), dtype=bool)
    periods = np.atleast_1d(periods)
    t0s = np.atleast_1d(t0s)
    depths = np.atleast_1d(depths)
    n_planets = len(periods)
    with pm.Model() as model:
        # Extract the un-masked data points
        model.x = x[mask]
        model.y = y[mask]
        model.yerr = (yerr + np.zeros_like(x))[mask]
        model.mask = mask

        # The baseline (out-of-transit) flux for the star in ppt. This
        # should be close to one because of how we normalized the data
        mean = pm.Normal("mean", mu=0.0, sd=10.0)
        logg_star = pm.Normal("logg_star", mu=toi_logg_star, sd=toi_logg_star_err)
        r_star = pm.Bound(pm.Normal, lower=0.0)("r_star", mu=toi_r_star, sd=toi_r_star_err)
        rho_star = pm.Deterministic("rho_star", factor * 10**logg_star / r_star)

        # The time of a reference transit for each planet
        t0 = pm.Normal("t0", mu=t0s, sd=1.0, shape=n_planets)

        # The log period; also tracking the period itself
        logP = pm.Normal("logP", mu=np.log(periods), sd=0.1, shape=n_planets)
        period = pm.Deterministic("period", tt.exp(logP))

        # The Kipping (2013) parameterization for quadratic limb darkening paramters
        u = xo.distributions.QuadLimbDark("u")

        # The Espinoza (2018) parameterization for the joint radius ratio and
        # impact parameter distribution
        r, b = xo.distributions.get_joint_radius_impact(
            min_radius=0.001, max_radius=1.0,
        r_pl = pm.Deterministic("r_pl", r * r_star)

        # This shouldn't make a huge difference, but I like to put a uniform
        # prior on the *log* of the radius ratio instead of the value. This
        # can be implemented by adding a custom "potential" (log probability).
        pm.Potential("r_prior", -pm.math.log(r))

        # Set up a Keplerian orbit for the planets
        model.orbit = xo.orbits.KeplerianOrbit(
            period=period, t0=t0, b=b, r_star=r_star, rho_star=rho_star)
        # Compute the model light curve using starry
        model.light_curves = xo.StarryLightCurve(u).get_light_curve(
            orbit=model.orbit, r=r_pl, t=model.x)
        model.light_curve = pm.math.sum(model.light_curves, axis=-1) * 1e3 + mean

        # Jitter and likelihood function
        logs2 = pm.Normal("logs2", mu=np.log(np.mean(model.yerr)), sd=10)
        pm.Normal("obs", mu=model.light_curve, sd=tt.sqrt(model.yerr**2+tt.exp(logs2)),

        # 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 = start        
        map_soln = xo.optimize(start=map_soln, vars=[logs2, mean])
        map_soln = xo.optimize(start=map_soln, vars=[model.rb, mean])
        map_soln = xo.optimize(start=map_soln, vars=[logg_star])
        map_soln = xo.optimize(start=map_soln, vars=[logP, t0, mean])
        map_soln = xo.optimize(start=map_soln, vars=[model.rb, mean])
        map_soln = xo.optimize(start=map_soln)
        model.map_soln = map_soln
    return model

def build_model_sigma_clip(x, y, yerr, periods, t0s, depths, sigma=5.0, maxiter=10, start=None):
    ntot = len(x)
    for i in range(maxiter):
        print("*** Sigma clipping round {0} ***".format(i+1))
        # Build the model
        model = build_model(x, y, yerr, periods, t0s, depths, start=start)
        start = model.map_soln

        # Compute the map prediction
        with model:
            mod = xo.utils.eval_in_model(model.light_curve, model.map_soln)
        # Do sigma clipping
        resid = y - mod
        rms = np.sqrt(np.median(resid**2))
        mask = np.abs(resid) < sigma * rms
        if ntot == mask.sum():
        ntot = mask.sum()

    return model

Using the above function, we'll generate a probabilistic model for the light curve and plot the maximum a posteriori fit.

In [7]:
model = build_model_sigma_clip(x, y, yerr, periods, t0s, depths)

with model:
    mean = model.map_soln["mean"]
    light_curves = xo.utils.eval_in_model(model.light_curves, model.map_soln)

plt.plot(model.x, model.y - mean, "k", label="data")
for n, l in enumerate(letters):
    plt.plot(model.x, 1e3 * light_curves[:, n], label="planet {0}".format(l), zorder=100-n)

plt.xlabel("time [days]")
plt.ylabel("flux [ppt]")
plt.title("initial fit")
plt.xlim(model.x.min(), model.x.max())
*** Sigma clipping round 1 ***
optimizing logp for variables: ['mean', 'logs2']
message: Optimization terminated successfully.
logp: -33275.3297112618 -> -33112.91919144288
optimizing logp for variables: ['mean', 'rb_radiusimpact__']
message: Optimization terminated successfully.
logp: -33112.91919144288 -> -31102.284421709057
optimizing logp for variables: ['logg_star']
message: Optimization terminated successfully.
logp: -31102.284421709057 -> -29221.035556821444
optimizing logp for variables: ['mean', 't0', 'logP']
message: Desired error not necessarily achieved due to precision loss.
logp: -29221.035556821444 -> -29173.812203185375
optimizing logp for variables: ['mean', 'rb_radiusimpact__']
message: Desired error not necessarily achieved due to precision loss.
logp: -29173.812203185375 -> -29121.445773680967
optimizing logp for variables: ['logs2', 'rb_radiusimpact__', 'u_quadlimbdark__', 'logP', 't0', 'r_star_lowerbound__', 'logg_star', 'mean']
message: Desired error not necessarily achieved due to precision loss.
logp: -29121.445773680967 -> -28667.6318455909
*** Sigma clipping round 2 ***
optimizing logp for variables: ['mean', 'logs2']
message: Desired error not necessarily achieved due to precision loss.
logp: -28667.6318455909 -> -27847.40936259271
optimizing logp for variables: ['mean', 'rb_radiusimpact__']
message: Desired error not necessarily achieved due to precision loss.
logp: -27847.40936259271 -> -27836.15165877355
optimizing logp for variables: ['logg_star']
message: Optimization terminated successfully.
logp: -27836.15165877355 -> -27820.36699970267
optimizing logp for variables: ['mean', 't0', 'logP']
message: Desired error not necessarily achieved due to precision loss.
logp: -27820.36699970267 -> -27819.91315804824
optimizing logp for variables: ['mean', 'rb_radiusimpact__']
message: Optimization terminated successfully.
logp: -27819.91315804824 -> -27818.751121191275
optimizing logp for variables: ['logs2', 'rb_radiusimpact__', 'u_quadlimbdark__', 'logP', 't0', 'r_star_lowerbound__', 'logg_star', 'mean']
message: Desired error not necessarily achieved due to precision loss.
logp: -27818.751121191275 -> -27746.61610774118
*** Sigma clipping round 3 ***
optimizing logp for variables: ['mean', 'logs2']
message: Desired error not necessarily achieved due to precision loss.
logp: -27746.616107741185 -> -27746.616107741185
optimizing logp for variables: ['mean', 'rb_radiusimpact__']
message: Desired error not necessarily achieved due to precision loss.
logp: -27746.61610774118 -> -27746.61610774118
optimizing logp for variables: ['logg_star']
message: Desired error not necessarily achieved due to precision loss.
logp: -27746.61610774118 -> -27746.61610774118
optimizing logp for variables: ['mean', 't0', 'logP']
message: Desired error not necessarily achieved due to precision loss.
logp: -27746.61610774118 -> -27746.61610774118
optimizing logp for variables: ['mean', 'rb_radiusimpact__']
message: Desired error not necessarily achieved due to precision loss.
logp: -27746.61610774118 -> -27746.61610774118
optimizing logp for variables: ['logs2', 'rb_radiusimpact__', 'u_quadlimbdark__', 'logP', 't0', 'r_star_lowerbound__', 'logg_star', 'mean']
message: Desired error not necessarily achieved due to precision loss.
logp: -27746.61610774118 -> -27746.616107740945


Now we use PyMC3 to sample the posterior density for the parameters of this model.

In [8]:
sampler = xo.PyMC3Sampler(window=50, start=50, finish=500)
with model:
    burnin = sampler.tune(tune=3000, start=model.map_soln,
                          chains=2, cores=2)
    trace = sampler.sample(draws=1000, chains=2, cores=2)
Sampling 2 chains: 100%|██████████| 104/104 [00:43<00:00,  2.39draws/s]
Sampling 2 chains: 100%|██████████| 104/104 [00:38<00:00,  1.05s/draws]
Sampling 2 chains: 100%|██████████| 204/204 [01:08<00:00,  1.41draws/s]
Sampling 2 chains: 100%|██████████| 404/404 [00:09<00:00, 42.47draws/s]
Sampling 2 chains: 100%|██████████| 804/804 [00:18<00:00, 18.57draws/s]
Sampling 2 chains: 100%|██████████| 4404/4404 [01:35<00:00, 29.65draws/s]
Sampling 2 chains: 100%|██████████| 1004/1004 [00:27<00:00, 19.45draws/s]
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [logs2, rb, u, logP, t0, r_star, logg_star, mean]
Sampling 2 chains: 100%|██████████| 2000/2000 [00:47<00:00, 28.79draws/s]

We can use the summary function to check convergence and report posterior means and variances for the parameters.

In [9]:
pm.summary(trace, varnames=["mean", "u", "period", "t0", "r", "b", "logs2", "r_star", "logg_star"])
mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
mean -0.061932 0.008287 1.701651e-04 -0.078196 -0.045428 2241.154899 1.000339
u__0 0.273758 0.155750 4.684446e-03 0.014131 0.547467 959.602046 0.999598
u__1 0.156244 0.221811 6.984692e-03 -0.242528 0.533786 899.910171 0.999684
period__0 2.652187 0.000045 9.041861e-07 2.652098 2.652267 2216.968574 1.000088
t0__0 1354.561111 0.000261 5.411704e-06 1354.560605 1354.561608 2139.337422 0.999616
r__0 0.102124 0.000798 2.364763e-05 0.100653 0.103640 1094.107212 0.999914
b__0 0.642021 0.019042 4.968625e-04 0.605791 0.679404 1467.950270 0.999680
logs2 -2.400465 0.146533 3.196056e-03 -2.685791 -2.115109 1923.103174 0.999912
r_star 1.674044 0.971972 2.169671e-02 0.005019 3.382671 1572.841060 0.999523
logg_star 3.850623 0.372392 1.325287e-02 3.041853 4.356822 575.830353 0.999545

And then plot the posterior constraints on the folded transit models.

In [10]:
with model:
    light_curves = np.empty((500, len(model.x), len(periods)))
    func = xo.utils.get_theano_function_for_var(model.light_curves)
    for i, sample in enumerate(xo.utils.get_samples_from_trace(
            trace, size=len(light_curves))):
        light_curves[i] = func(*xo.utils.get_args_for_theano_function(sample))

for n, letter in enumerate(letters):

    # Compute the GP prediction
    mean_mod = np.median(trace["mean"][:, None])

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

    # Compute the median of posterior estimate of the contribution from
    # the other planet. Then we can remove this from the data to plot
    # just the planet we care about.
    inds = np.arange(len(periods)) != n
    others = np.median(1e3*np.sum(light_curves[:, :, inds], axis=-1), axis=0)

    # Plot the folded data
    x_fold = (model.x - t0 + 0.5*p) % p - 0.5*p
    plt.plot(x_fold, model.y - mean_mod - others, ".k", label="data", zorder=-1000)

    # Plot the folded model
    inds = np.argsort(x_fold)
    inds = inds[np.abs(x_fold)[inds] < 0.3]
    pred = 1e3 * light_curves[:, inds, n]
    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,

    # Annotate the plot with the planet's period
    txt = "period = {0:.4f} +/- {1:.4f} d".format(
        np.mean(trace["period"][:, n]), np.std(trace["period"][:, n]))
    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.title("TOI {0}{1}".format(toi_num, letter));
    plt.xlim(-0.3, 0.3)

Posterior constraints

Given the MCMC sampling, we can make some plots summarizing the constraints on the key parameters. First up, we plot the physical radius of the planet using the stellar radius constraint from the TOI list, and the impact parameter (remember that we're using circular orbits).

In [11]:
# Convert to Earth radii
r_pl = trace["r_pl"] * 109.07637070600963
samples = np.concatenate((r_pl, trace["b"]), axis=-1)

labels = ["$R_{{\mathrm{{Pl}},{0}}}$ [$R_\oplus$]".format(i) for i in letters]
labels += ["impact param {0}".format(i) for i in letters]

corner.corner(samples, labels=labels,
              show_titles=True, title_kwargs=dict(fontsize=10));

The other most interesting parameters are the period and transit times.

In [12]:
labels = ["$P_{{{0}}}$ [days]".format(i) for i in letters]
labels += ["$t0_{{{0}}}$ [TBJD]".format(i) for i in letters]
samples = np.concatenate((trace["period"], trace["t0"]), axis=-1)
corner.corner(samples, labels=labels,
              show_titles=True, title_fmt=".5f",

Finally, here are the posterior constraints on the stellar properties. These won't be exactly the same as the prior because the duration of the transits (and the assumption of circular orbits )

In [13]:
labels = ["$R_\mathrm{star}$ [$R_\odot$]", "$\log g$ [cm/s$^2$]",
          r"$\rho_\mathrm{star}$ [g/cm$^3$]"]
samples = np.vstack((trace["r_star"], trace["logg_star"], trace["rho_star"])).T
corner.corner(samples, labels=labels,


If you use these results or this code, please consider citing the relevant sources. First, you can cite the lightkurve package:

  author       = {Geert Barentsen and
                  Christina Hedges and
                  Zé Vinícius and
                  Nicholas Saunders and
                  gully and
                  Oliver Hall and
                  Sheila Sagear and
                  Tom Barclay and
                  KenMighell and
                  Keaton Bell and
                  Johnny Zhang and
                  Emma Turtelboom and
                  Zach Berta-Thompson and
                  Peter Williams and
                  Jose A Lerma III and
                  Guy Davies and
                  Brennan Vincello and
                  Anand Sundaram},
  title        = {KeplerGO/lightkurve: Lightkurve v1.0b30},
  month        = mar,
  year         = 2019,
  doi          = {10.5281/zenodo.2611871},
  url          = {https://doi.org/10.5281/zenodo.2611871}

You can also cite the exoplanet project and its dependencies using the following acknowledgement:

In [14]:
with model:
    txt, bib = xo.citations.get_citations_for_model()
This research made use of \textsf{exoplanet} \citep{exoplanet} and its
dependencies \citep{exoplanet:astropy13, exoplanet:astropy18,
exoplanet:espinoza18, exoplanet:exoplanet, exoplanet:kipping13,
exoplanet:luger18, exoplanet:pymc3, exoplanet:theano}.

and BibTeX entries:

In [15]:
  author = {Dan Foreman-Mackey and
            Geert Barentsen and
            Tom Barclay},
   title = {dfm/exoplanet: exoplanet v0.1.5},
   month = mar,
    year = 2019,
     doi = {10.5281/zenodo.2587222},
     url = {https://doi.org/10.5281/zenodo.2587222}

    title={Probabilistic programming in Python using PyMC3},
   author={Salvatier, John and Wiecki, Thomas V and Fonnesbeck, Christopher},
  journal={PeerJ Computer Science},
publisher={PeerJ Inc.}

    title="{Theano: A {Python} framework for fast computation of mathematical
   author={{Theano Development Team}},
  journal={arXiv e-prints},

   author = {{Kipping}, D.~M.},
    title = "{Efficient, uninformative sampling of limb darkening coefficients
              for two-parameter laws}",
  journal = {\mnras},
     year = 2013,
    month = nov,
   volume = 435,
    pages = {2152-2160},
      doi = {10.1093/mnras/stt1435},
   adsurl = {http://adsabs.harvard.edu/abs/2013MNRAS.435.2152K},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

   author = {{Espinoza}, N.},
    title = "{Efficient Joint Sampling of Impact Parameters and Transit Depths
              in Transiting Exoplanet Light Curves}",
  journal = {Research Notes of the American Astronomical Society},
     year = 2018,
    month = nov,
   volume = 2,
   number = 4,
    pages = {209},
      doi = {10.3847/2515-5172/aaef38},
   adsurl = {http://adsabs.harvard.edu/abs/2018RNAAS...2d.209E},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

   author = {{Astropy Collaboration} and {Robitaille}, T.~P. and {Tollerud},
             E.~J. and {Greenfield}, P. and {Droettboom}, M. and {Bray}, E. and
             {Aldcroft}, T. and {Davis}, M. and {Ginsburg}, A. and
             {Price-Whelan}, A.~M. and {Kerzendorf}, W.~E. and {Conley}, A. and
             {Crighton}, N. and {Barbary}, K. and {Muna}, D. and {Ferguson}, H.
             and {Grollier}, F. and {Parikh}, M.~M. and {Nair}, P.~H. and
             {Unther}, H.~M. and {Deil}, C. and {Woillez}, J. and {Conseil}, S.
             and {Kramer}, R. and {Turner}, J.~E.~H. and {Singer}, L. and
             {Fox}, R. and {Weaver}, B.~A. and {Zabalza}, V. and {Edwards},
             Z.~I. and {Azalee Bostroem}, K. and {Burke}, D.~J. and {Casey},
             A.~R. and {Crawford}, S.~M. and {Dencheva}, N. and {Ely}, J. and
             {Jenness}, T. and {Labrie}, K. and {Lim}, P.~L. and
             {Pierfederici}, F. and {Pontzen}, A. and {Ptak}, A. and {Refsdal},
             B. and {Servillat}, M. and {Streicher}, O.},
    title = "{Astropy: A community Python package for astronomy}",
  journal = {\aap},
     year = 2013,
    month = oct,
   volume = 558,
    pages = {A33},
      doi = {10.1051/0004-6361/201322068},
   adsurl = {http://adsabs.harvard.edu/abs/2013A%26A...558A..33A},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

   author = {{Astropy Collaboration} and {Price-Whelan}, A.~M. and
             {Sip{\H o}cz}, B.~M. and {G{\"u}nther}, H.~M. and {Lim}, P.~L. and
             {Crawford}, S.~M. and {Conseil}, S. and {Shupe}, D.~L. and
             {Craig}, M.~W. and {Dencheva}, N. and {Ginsburg}, A. and
             {VanderPlas}, J.~T. and {Bradley}, L.~D. and
             {P{\'e}rez-Su{\'a}rez}, D. and {de Val-Borro}, M.
             and {Aldcroft}, T.~L. and {Cruz}, K.~L. and {Robitaille}, T.~P.
             and {Tollerud}, E.~J. and {Ardelean}, C. and {Babej}, T. and
             {Bach}, Y.~P. and {Bachetti}, M. and {Bakanov}, A.~V. and
             {Bamford}, S.~P. and {Barentsen}, G. and {Barmby}, P. and
             {Baumbach}, A. and {Berry}, K.~L.  and {Biscani}, F. and
             {Boquien}, M. and {Bostroem}, K.~A. and {Bouma}, L.~G. and
             {Brammer}, G.~B. and {Bray}, E.~M. and {Breytenbach}, H. and
             {Buddelmeijer}, H. and {Burke}, D.~J. and {Calderone}, G. and
             {Cano Rodr{\'{\i}}guez}, J.~L. and {Cara}, M. and {Cardoso},
             J.~V.~M. and {Cheedella}, S. and {Copin}, Y. and {Corrales}, L.
             and {Crichton}, D. and {D'Avella}, D. and {Deil}, C. and
             {Depagne}, {\'E}. and {Dietrich}, J.~P. and {Donath}, A. and
             {Droettboom}, M. and {Earl}, N. and {Erben}, T. and {Fabbro}, S.
             and {Ferreira}, L.~A. and {Finethy}, T. and {Fox}, R.~T. and
             {Garrison}, L.~H. and {Gibbons}, S.~L.~J. and {Goldstein}, D.~A.
             and {Gommers}, R. and {Greco}, J.~P. and {Greenfield}, P. and
             {Groener}, A.~M. and {Grollier}, F. and {Hagen}, A. and {Hirst},
             P. and {Homeier}, D. and {Horton}, A.~J. and {Hosseinzadeh}, G.
             and {Hu}, L. and {Hunkeler}, J.~S. and {Ivezi{\'c}}, {\v Z}. and
             {Jain}, A. and {Jenness}, T. and {Kanarek}, G. and {Kendrew}, S.
             and {Kern}, N.~S. and {Kerzendorf}, W.~E. and {Khvalko}, A. and
             {King}, J. and {Kirkby}, D. and {Kulkarni}, A.~M. and {Kumar}, A.
             and {Lee}, A.  and {Lenz}, D.  and {Littlefair}, S.~P. and {Ma},
             Z. and {Macleod}, D.~M. and {Mastropietro}, M. and {McCully}, C.
             and {Montagnac}, S. and {Morris}, B.~M. and {Mueller}, M. and
             {Mumford}, S.~J. and {Muna}, D. and {Murphy}, N.~A. and {Nelson},
             S. and {Nguyen}, G.~H. and {Ninan}, J.~P. and {N{\"o}the}, M. and
             {Ogaz}, S. and {Oh}, S. and {Parejko}, J.~K.  and {Parley}, N. and
             {Pascual}, S. and {Patil}, R. and {Patil}, A.~A.  and {Plunkett},
             A.~L. and {Prochaska}, J.~X. and {Rastogi}, T. and {Reddy Janga},
             V. and {Sabater}, J.  and {Sakurikar}, P. and {Seifert}, M. and
             {Sherbert}, L.~E. and {Sherwood-Taylor}, H. and {Shih}, A.~Y. and
             {Sick}, J. and {Silbiger}, M.~T. and {Singanamalla}, S. and
             {Singer}, L.~P. and {Sladen}, P.~H. and {Sooley}, K.~A. and
             {Sornarajah}, S. and {Streicher}, O. and {Teuben}, P. and
             {Thomas}, S.~W. and {Tremblay}, G.~R. and {Turner}, J.~E.~H. and
             {Terr{\'o}n}, V.  and {van Kerkwijk}, M.~H. and {de la Vega}, A.
             and {Watkins}, L.~L. and {Weaver}, B.~A. and {Whitmore}, J.~B. and
             {Woillez}, J.  and {Zabalza}, V. and {Astropy Contributors}},
    title = "{The Astropy Project: Building an Open-science Project and Status
              of the v2.0 Core Package}",
  journal = {\aj},
     year = 2018,
    month = sep,
   volume = 156,
    pages = {123},
      doi = {10.3847/1538-3881/aabc4f},
   adsurl = {http://adsabs.harvard.edu/abs/2018AJ....156..123A},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}

   author = {{Luger}, R. and {Agol}, E. and {Foreman-Mackey}, D. and {Fleming},
             D.~P. and {Lustig-Yaeger}, J. and {Deitrick}, R.},
    title = "{starry: Analytic Occultation Light Curves}",
  journal = {\aj},
     year = 2019,
    month = feb,
   volume = 157,
    pages = {64},
      doi = {10.3847/1538-3881/aae8e5},
   adsurl = {http://adsabs.harvard.edu/abs/2019AJ....157...64L},
  adsnote = {Provided by the SAO/NASA Astrophysics Data System}


This notebook was run with the following conda environment:

In [16]:
!conda env export
name: tessatlas
  - defaults
  - astropy=3.1.2=py37h7b6447c_0
  - atomicwrites=1.3.0=py37_1
  - attrs=19.1.0=py37_1
  - blas=1.0=mkl
  - ca-certificates=2019.1.23=0
  - certifi=2019.3.9=py37_0
  - cycler=0.10.0=py37_0
  - dbus=1.13.6=h746ee38_0
  - expat=2.2.6=he6710b0_0
  - fontconfig=2.13.0=h9420a91_0
  - freetype=2.9.1=h8a8886c_1
  - glib=2.56.2=hd408876_0
  - gst-plugins-base=1.14.0=hbbd80ab_1
  - gstreamer=1.14.0=hb453b48_1
  - icu=58.2=h9c2bf20_1
  - intel-openmp=2019.3=199
  - jpeg=9b=h024ee3a_2
  - kiwisolver=1.0.1=py37hf484d3e_0
  - libedit=3.1.20181209=hc058e9b_0
  - libffi=3.2.1=hd88cf55_4
  - libgcc-ng=8.2.0=hdf63c60_1
  - libgfortran-ng=7.3.0=hdf63c60_0
  - libpng=1.6.36=hbc83047_0
  - libstdcxx-ng=8.2.0=hdf63c60_1
  - libuuid=1.0.3=h1bed415_2
  - libxcb=1.13=h1bed415_1
  - libxml2=2.9.9=he19cac6_0
  - matplotlib=3.0.3=py37h5429711_0
  - mkl=2019.3=199
  - mkl_fft=1.0.10=py37ha843d7b_0
  - mkl_random=1.0.2=py37hd81dba3_0
  - more-itertools=6.0.0=py37_0
  - ncurses=6.1=he6710b0_1
  - numpy=1.16.2=py37h7e9f1db_0
  - numpy-base=1.16.2=py37hde5b4d6_0
  - openssl=1.1.1b=h7b6447c_1
  - pandas=0.24.2=py37he6710b0_0
  - pcre=8.43=he6710b0_0
  - pip=19.0.3=py37_0
  - pluggy=0.9.0=py37_0
  - psutil=5.6.1=py37h7b6447c_0
  - py=1.8.0=py37_0
  - pyparsing=2.3.1=py37_0
  - pyqt=5.9.2=py37h05f1152_2
  - pytest=4.3.1=py37_0
  - pytest-arraydiff=0.3=py37h39e3cac_0
  - pytest-astropy=0.5.0=py37_0
  - pytest-doctestplus=0.3.0=py37_0
  - pytest-openfiles=0.3.2=py37_0
  - pytest-remotedata=0.3.1=py37_0
  - python=3.7.3=h0371630_0
  - python-dateutil=2.8.0=py37_0
  - pytz=2018.9=py37_0
  - qt=5.9.7=h5867ecd_1
  - readline=7.0=h7b6447c_5
  - scipy=1.2.1=py37h7c811a0_0
  - setuptools=40.8.0=py37_0
  - sip=4.19.8=py37hf484d3e_0
  - six=1.12.0=py37_0
  - sqlite=3.27.2=h7b6447c_0
  - tk=8.6.8=hbc83047_0
  - tornado=6.0.2=py37h7b6447c_0
  - wheel=0.33.1=py37_0
  - xz=5.2.4=h14c3975_4
  - zlib=1.2.11=h7b6447c_3
  - pip:
    - alembic==1.0.8
    - asn1crypto==0.24.0
    - astroquery==0.3.9
    - async-generator==1.10
    - autograd==1.2
    - backcall==0.1.0
    - beautifulsoup4==4.7.1
    - bleach==3.1.0
    - bs4==0.0.1
    - cffi==1.12.2
    - chardet==3.0.4
    - corner==2.0.1
    - cryptography==2.6.1
    - decorator==4.4.0
    - defusedxml==0.5.0
    - entrypoints==0.3
    - exoplanet==0.1.5
    - fbpca==1.0
    - future==0.17.1
    - h5py==2.9.0
    - html5lib==1.0.1
    - idna==2.8
    - ipykernel==5.1.0
    - ipython==7.4.0
    - ipython-genutils==0.2.0
    - ipywidgets==7.4.2
    - jedi==0.13.3
    - jeepney==0.4
    - jinja2==2.10
    - joblib==0.12.5
    - jsonschema==3.0.1
    - jupyter==1.0.0
    - jupyter-client==5.2.4
    - jupyter-console==6.0.0
    - jupyter-core==4.4.0
    - jupyterhub==0.9.5
    - jupyterlab==0.35.4
    - jupyterlab-server==0.2.0
    - keyring==19.0.1
    - lightkurve==1.0b30
    - mako==1.0.8
    - markupsafe==1.1.1
    - mistune==0.8.4
    - nbconvert==5.4.1
    - nbformat==4.4.0
    - notebook==5.7.7
    - oktopus==0.1.2
    - pamela==1.0.0
    - pandocfilters==1.4.2
    - parso==0.3.4
    - patsy==0.5.1
    - pexpect==4.6.0
    - pickleshare==0.7.5
    - prometheus-client==0.6.0
    - prompt-toolkit==2.0.9
    - ptyprocess==0.6.0
    - pycparser==2.19
    - pygments==2.3.1
    - pymc3==3.6
    - pyrsistent==0.14.11
    - python-editor==1.0.4
    - python-oauth2==1.1.0
    - pyzmq==18.0.1
    - qtconsole==4.4.3
    - requests==2.21.0
    - secretstorage==3.1.1
    - send2trash==1.5.0
    - soupsieve==1.9
    - sqlalchemy==1.3.1
    - terminado==0.8.1
    - testpath==0.4.2
    - theano==1.0.4
    - tqdm==4.31.1
    - traitlets==4.3.2
    - urllib3==1.24.1
    - wcwidth==0.1.7
    - webencodings==0.5.1
    - widgetsnbextension==3.4.2
prefix: /mnt/home/dforeman/miniconda3/envs/tessatlas