{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "b904abe7",
"metadata": {
"ploomber": {
"memory_usage": 198.515625,
"timestamp_end": 1691727395.9039,
"timestamp_start": 1691727391.274553
},
"tags": [
"exe",
"hide-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": []
}
],
"source": [
"! pip install tess-atlas -q\n",
"\n",
"from tess_atlas.utils import notebook_initalisations\n",
"\n",
"notebook_initalisations()"
]
},
{
"cell_type": "markdown",
"id": "d93383fa",
"metadata": {
"tags": [
"def"
]
},
"source": [
"# TESS Atlas fit for TOI 3114\n",
"\n",
"**Version: '0.2.1.dev314+ge1761e8'**\n",
"\n",
"**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](https://github.com/dfm/tess-atlas)\n",
"\n",
"In this notebook, we do a quicklook fit for the parameters of the TESS Objects of Interest (TOI) in the system number 3114.\n",
"To do this fit, we use the [exoplanet](https://exoplanet.dfm.io) library and you can find more information about that project at [exoplanet.dfm.io](https://exoplanet.dfm.io).\n",
"\n",
"From here, you can\n",
"- scroll down and take a look at the fit results\n",
"- open the notebook in Google Colab to run the fit yourself\n",
"- download the notebook\n",
"\n",
"\n",
"\n",
"## Caveats\n",
"\n",
"There are many caveats associated with this relatively simple \"quicklook\" type of analysis that should be kept in mind.\n",
"Here are some of the main things that come to mind:\n",
"\n",
"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.\n",
"\n",
"2. Transit timing variations, correlated noise, and (probably) your favorite systematics are ignored. Sorry!\n",
"\n",
"3. This notebook was generated automatically without human intervention. Use at your own risk!\n",
"\n",
"\n",
"## Getting started\n",
"\n",
"To get going, we'll add some _magic_, import some packages, and run some setup steps."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "835d0bbe",
"metadata": {
"ploomber": {
"memory_usage": 198.69140625,
"timestamp_end": 1691727395.987593,
"timestamp_start": 1691727395.90526
},
"pycharm": {
"name": "#%%\n"
},
"tags": [
"def",
"hide-cell",
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The autoreload extension is already loaded. To reload it, use:\n",
" %reload_ext autoreload\n",
"time: 3.37 ms (started: 2023-08-11 14:16:35 +10:00)\n"
]
}
],
"source": [
"%load_ext autoreload\n",
"%load_ext memory_profiler\n",
"%load_ext autotime\n",
"%autoreload 2\n",
"# %matplotlib inline\n",
"\n",
"import os\n",
"\n",
"import numpy as np\n",
"import pymc3_ext as pmx\n",
"from arviz import InferenceData\n",
"\n",
"from tess_atlas.analysis.eccenticity_reweighting import (\n",
" calculate_eccentricity_weights,\n",
")\n",
"from tess_atlas.analysis.model_tools import (\n",
" compute_variable,\n",
" get_untransformed_varnames,\n",
" sample_prior,\n",
")\n",
"from tess_atlas.data.inference_data_tools import (\n",
" get_optimized_init_params,\n",
" summary,\n",
" test_model,\n",
")\n",
"from tess_atlas.data.tic_entry import TICEntry\n",
"from tess_atlas.logger import get_notebook_logger\n",
"from tess_atlas.plotting import (\n",
" plot_diagnostics,\n",
" plot_eccentricity_posteriors,\n",
" plot_inference_trace,\n",
" plot_lightcurve,\n",
" plot_phase,\n",
" plot_posteriors,\n",
" plot_priors,\n",
" plot_raw_lightcurve,\n",
")\n",
"\n",
"TOI_NUMBER = 3114\n",
"logger = get_notebook_logger(outdir=f\"toi_{TOI_NUMBER}_files\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "37cf57a3",
"metadata": {
"ploomber": {
"memory_usage": 198.796875,
"timestamp_end": 1691727396.064007,
"timestamp_start": 1691727395.988083
},
"tags": [
"exe",
"remove-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 9.92 ms (started: 2023-08-11 14:16:36 +10:00)\n"
]
}
],
"source": [
"import theano\n",
"\n",
"from tess_atlas.utils import tabulate_global_environ_vars\n",
"\n",
"logger.info(\"Logging some settings for future reference\")\n",
"logger.info(\"GLOBAL ENVS:\\n\" + tabulate_global_environ_vars())\n",
"logger.info(f\"THEANO Config:\\n{theano.config}\")"
]
},
{
"cell_type": "markdown",
"id": "b0779ec5",
"metadata": {
"lines_to_next_cell": 0,
"tags": [
"def"
]
},
"source": [
"## Downloading Data\n",
"\n",
"Next, we grab some inital guesses for the TOI's parameters from [ExoFOP](https://exofop.ipac.caltech.edu/tess/) and download the TOI's lightcurve with [Lightkurve].\n",
"\n",
"We wrap the information in three objects, a `TIC Entry`, a `Planet Candidate` and finally a `Lightcurve Data` object.\n",
"\n",
"- The `TIC Entry` object holds one or more `Planet Candidate`s (each candidate associated with one TOI id number) and a `Lightcurve Data` for associated with the candidates. Note that the `Lightcurve Data` object is initially the same fopr each candidate but may be masked according to the candidate transit's period.\n",
"\n",
"- The `Planet Candidate` holds information on the TOI data collected by [SPOC] (eg transit period, etc)\n",
"\n",
"- The `Lightcurve Data` holds the lightcurve time and flux data for the planet candidates.\n",
"\n",
"[ExoFOP]: https://exofop.ipac.caltech.edu/tess/\n",
"[Lightkurve]: https://docs.lightkurve.org/index.html\n",
"[SPOC]: https://heasarc.gsfc.nasa.gov/docs/tess/pipeline.html\n",
"\n",
"Downloading the data (this may take a few minutes):"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "d18e41d5",
"metadata": {
"ploomber": {
"memory_usage": 214.73046875,
"timestamp_end": 1691727396.74997,
"timestamp_start": 1691727396.064554
},
"pycharm": {
"name": "#%%\n"
},
"tags": [
"exe"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 620 ms (started: 2023-08-11 14:16:36 +10:00)\n"
]
},
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
TOI
\n",
"
3114.01
\n",
"
\n",
" \n",
" \n",
"
\n",
"
TOI
\n",
"
3114.01
\n",
"
\n",
"
\n",
"
Classification
\n",
"
Planet Candidate
\n",
"
\n",
"
\n",
"
Period (days)
\n",
"
4.407111
\n",
"
\n",
"
\n",
"
Epoch (TBJD)
\n",
"
2331.60796
\n",
"
\n",
"
\n",
"
Depth (ppt)
\n",
"
9.14
\n",
"
\n",
"
\n",
"
Duration (days)
\n",
"
0.088917
\n",
"
\n",
"
\n",
"
Planet SNR
\n",
"
61.0
\n",
"
\n",
"
\n",
"
Single Transit
\n",
"
False
\n",
"
\n",
"
\n",
"
PE Pipeline
\n",
"
QLP
\n",
"
\n",
"
\n",
"
Comments
\n",
"
found in faint-star QLP search
\n",
"
\n",
" \n",
"
\n",
"
\n",
"More data on ExoFOP page"
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tic_entry = TICEntry.load(toi=TOI_NUMBER)\n",
"tic_entry"
]
},
{
"cell_type": "markdown",
"id": "1fc52577",
"metadata": {},
"source": [
"If the amount of lightcurve data availible is large we filter the data to keep only data around transits."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "6698b984",
"metadata": {
"ploomber": {
"memory_usage": 214.734375,
"timestamp_end": 1691727396.821079,
"timestamp_start": 1691727396.750513
},
"tags": [
"exe"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 1.08 ms (started: 2023-08-11 14:16:36 +10:00)\n"
]
}
],
"source": [
"if tic_entry.lightcurve.len > 1e5:\n",
" tic_entry.lightcurve.filter_non_transit_data(tic_entry.candidates)\n",
"else:\n",
" logger.info(\"Using the full lightcurve for analysis.\")"
]
},
{
"cell_type": "markdown",
"id": "ba99a5cc",
"metadata": {
"lines_to_next_cell": 2
},
"source": [
"Plot of the lightcurve:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "df6538fc",
"metadata": {
"lines_to_next_cell": 0,
"ploomber": {
"memory_usage": 234.9140625,
"timestamp_end": 1691727398.895523,
"timestamp_start": 1691727396.821544
},
"pycharm": {
"name": "#%%\n"
},
"tags": [
"exe",
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 2 s (started: 2023-08-11 14:16:36 +10:00)\n"
]
}
],
"source": [
"plot_lightcurve(tic_entry, save=True)\n",
"\n",
"# Some diagnostics\n",
"plot_raw_lightcurve(tic_entry, save=True)\n",
"plot_raw_lightcurve(tic_entry, zoom_in=True, save=True)"
]
},
{
"cell_type": "markdown",
"id": "23bf0381",
"metadata": {},
"source": [
"![](toi_3114_files/flux_vs_time.png)"
]
},
{
"cell_type": "markdown",
"id": "92e2eea0",
"metadata": {
"tags": [
"hide-cell"
]
},
"source": [
"Diagnostic plots of the raw lightcurve (not applying sigma clipping/other cleaning methods to remove outliers). Some things to consider:\n",
"- Do the initial fits from ExoFOP match the transits if visible?\n",
"- If this is marked as a single-transit event, is there only 1 transit visible?\n",
"\n",
"![](toi_3114_files/diagnostic_raw_lc_plot.png)\n",
"![](toi_3114_files/diagnostic_raw_lc_plot_zoom.png)"
]
},
{
"cell_type": "markdown",
"id": "145cb607",
"metadata": {},
"source": [
"## Fitting transit parameters\n",
"Now that we have the data, we can define a Bayesian model to fit it.\n",
"\n",
"### The probabilistic model\n",
"\n",
"We use the probabilistic model as described in [Foreman-Mackey et al 2017] to determine the best parameters to fit the transits present in the lightcurve data.\n",
"\n",
"More explicitly, the stellar light curve $l(t; \\vec{\\theta})$ is modelled with a Gaussian Process (GP).\n",
"A GP consists of a mean function $\\mu(t;\\vec{\\theta})$ and a kernel function $k_\\alpha(t,t';\\vec{\\theta})$, where $\\vec{\\theta}$ is the vector of parameters descibing the lightcurve and $t$ is the time during which the lightcurve is under observation\n",
"\n",
"The 8 parameters describing the lightcurve are\n",
"$$\\vec{\\theta} = \\{d_i, t0_i, tmax_i, b_i, r_i, f0, u1, u2\\},$$\n",
"where\n",
"* $d_i$ transit durations for each planet,\n",
"* $tmin_i$ time of first transit for each planet (reference time),\n",
"* $tmax_i$ time of the last transit for each planet (a second reference time),\n",
"* $b_i$ impact parameter for each planet,\n",
"* $r_i$ planet radius in stellar radius for each planet,\n",
"* $f0$ baseline relative flux of the light curve from star,\n",
"* $u1$ $u2$ two parameters describing the limb-darkening profile of star.\n",
"\n",
"Note: if the observed data only records a single transit,\n",
"we swap $tmax_i$ with $p_i$ (orbital periods for each planet).\n",
"\n",
"With this we can write\n",
"$$l(t;\\vec{\\theta}) \\sim \\mathcal{GP} (\\mu(t;\\vec{\\theta}), k_\\alpha(t,t';\\vec{\\theta}))\\ .$$\n",
"\n",
"Here the mean and kernel functions are:\n",
"* $\\mu(t;\\vec{\\theta})$: a limb-darkened transit light curve ([Kipping 2013])\n",
"* $k_\\alpha(t,t';\\vec{\\theta}))$: a stochastically-driven, damped harmonic oscillator ([SHOTterm])\n",
"\n",
"\n",
"Now that we have defined our transit model, we can implement it in python (toggle to show).\n",
"\n",
"[Foreman-Mackey et al 2017]: https://arxiv.org/pdf/1703.09710.pdf\n",
"[Kipping 2013]: https://arxiv.org/abs/1308.0009\n",
"[SHOTterm]: https://celerite2.readthedocs.io/en/latest/api/python/?highlight=SHOTerm#celerite2.terms.SHOTerm"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "e1632350",
"metadata": {
"ploomber": {
"memory_usage": 236.9375,
"timestamp_end": 1691727399.012105,
"timestamp_start": 1691727398.896129
},
"pycharm": {
"name": "#%%\n"
},
"tags": [
"def",
"hide-cell"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 39.9 ms (started: 2023-08-11 14:16:38 +10:00)\n"
]
}
],
"source": [
"import aesara_theano_fallback.tensor as tt\n",
"import exoplanet as xo\n",
"import numpy as np\n",
"import pymc3 as pm\n",
"import pymc3_ext as pmx\n",
"from celerite2.theano import GaussianProcess, terms\n",
"\n",
"DEPTH = \"depth\"\n",
"DURATION = \"dur\"\n",
"RADIUS_RATIO = \"r\"\n",
"TIME_START = \"tmin\"\n",
"TIME_END = \"tmax\"\n",
"ORBITAL_PERIOD = \"p\"\n",
"MEAN_FLUX = \"f0\"\n",
"LC_JITTER = \"jitter\"\n",
"GP_RHO = \"rho\"\n",
"GP_SIGMA = \"sigma\"\n",
"RHO_CIRC = \"rho_circ\" # stellar density at e=0\n",
"LIMB_DARKENING_PARAM = \"u\"\n",
"IMPACT_PARAM = \"b\"\n",
"\n",
"\n",
"def get_test_duration(min_durations, max_durations, durations):\n",
" largest_min_duration = np.amax(\n",
" np.array([durations, 2 * min_durations]), axis=0\n",
" )\n",
" smallest_max_duration = np.amin(\n",
" np.array([largest_min_duration, 0.99 * max_durations]), axis=0\n",
" )\n",
" return smallest_max_duration\n",
"\n",
"\n",
"def build_planet_transit_model(tic_entry):\n",
" t = tic_entry.lightcurve.time\n",
" y = tic_entry.lightcurve.flux\n",
" yerr = tic_entry.lightcurve.flux_err\n",
"\n",
" n = tic_entry.planet_count\n",
" tmins = np.array([planet.tmin for planet in tic_entry.candidates])\n",
" depths = np.array([planet.depth for planet in tic_entry.candidates])\n",
" durations = np.array([planet.duration for planet in tic_entry.candidates])\n",
" max_durations = np.array(\n",
" [planet.duration_max for planet in tic_entry.candidates]\n",
" )\n",
" min_durations = np.array(\n",
" [planet.duration_min for planet in tic_entry.candidates]\n",
" )\n",
" test_duration = get_test_duration(min_durations, max_durations, durations)\n",
"\n",
" with pm.Model() as my_planet_transit_model:\n",
" ## define planet parameters\n",
"\n",
" # 1) d: transit duration (duration of eclipse)\n",
" d_priors = pm.Bound(\n",
" pm.Lognormal, lower=min_durations, upper=max_durations\n",
" )(\n",
" name=DURATION,\n",
" mu=np.log(durations),\n",
" sigma=np.log(1.2),\n",
" shape=n,\n",
" testval=test_duration,\n",
" )\n",
"\n",
" # 2) r: radius ratio (planet radius / star radius)\n",
" r_priors = pm.Lognormal(\n",
" name=RADIUS_RATIO, mu=0.5 * np.log(depths * 1e-3), sd=1.0, shape=n\n",
" )\n",
" # 3) b: impact parameter\n",
" b_priors = xo.distributions.ImpactParameter(\n",
" name=IMPACT_PARAM, ror=r_priors, shape=n\n",
" )\n",
" planet_priors = [r_priors, d_priors, b_priors]\n",
"\n",
" ## define orbit-timing parameters\n",
"\n",
" # 1) tmin: the time of the first transit in data (a reference time)\n",
" tmin_norm = pm.Bound(\n",
" pm.Normal, lower=tmins - max_durations, upper=tmins + max_durations\n",
" )\n",
" tmin_priors = tmin_norm(\n",
" TIME_START, mu=tmins, sigma=0.5 * durations, shape=n, testval=tmins\n",
" )\n",
"\n",
" # 2) period: the planets' orbital period\n",
" p_params, p_priors_list, tmax_priors_list = [], [], []\n",
" for n, planet in enumerate(tic_entry.candidates):\n",
" # if only one transit in data we use the period\n",
" if planet.has_data_only_for_single_transit:\n",
" p_prior = pm.Pareto(\n",
" name=f\"{ORBITAL_PERIOD}_{planet.index}\",\n",
" m=planet.period_min,\n",
" alpha=2.0 / 3.0,\n",
" testval=planet.period,\n",
" )\n",
" p_param = p_prior\n",
" tmax_prior = planet.tmin\n",
" # if more than one transit in data we use a second time reference (tmax)\n",
" else:\n",
" tmax_norm = pm.Bound(\n",
" pm.Normal,\n",
" lower=planet.tmax - planet.duration_max,\n",
" upper=planet.tmax + planet.duration_max,\n",
" )\n",
" tmax_prior = tmax_norm(\n",
" name=f\"{TIME_END}_{planet.index}\",\n",
" mu=planet.tmax,\n",
" sigma=0.5 * planet.duration,\n",
" testval=planet.tmax,\n",
" )\n",
" p_prior = (tmax_prior - tmin_priors[n]) / planet.num_periods\n",
" p_param = tmax_prior\n",
"\n",
" p_params.append(p_param) # the param needed to calculate p\n",
" p_priors_list.append(p_prior)\n",
" tmax_priors_list.append(tmax_prior)\n",
"\n",
" p_priors = pm.Deterministic(ORBITAL_PERIOD, tt.stack(p_priors_list))\n",
" tmax_priors = pm.Deterministic(TIME_END, tt.stack(tmax_priors_list))\n",
"\n",
" ## define stellar parameters\n",
"\n",
" # 1) f0: the mean flux from the star\n",
" f0_prior = pm.Normal(name=MEAN_FLUX, mu=0.0, sd=10.0)\n",
"\n",
" # 2) u1, u2: limb darkening parameters\n",
" u_prior = xo.distributions.QuadLimbDark(\"u\")\n",
" stellar_priors = [f0_prior, u_prior]\n",
"\n",
" ## define k(t, t1; parameters)\n",
" jitter_prior = pm.InverseGamma(\n",
" name=LC_JITTER, **pmx.estimate_inverse_gamma_parameters(1.0, 5.0)\n",
" )\n",
" sigma_prior = pm.InverseGamma(\n",
" name=GP_SIGMA, **pmx.estimate_inverse_gamma_parameters(1.0, 5.0)\n",
" )\n",
" rho_prior = pm.InverseGamma(\n",
" name=GP_RHO, **pmx.estimate_inverse_gamma_parameters(0.5, 10.0)\n",
" )\n",
" kernel = terms.SHOTerm(sigma=sigma_prior, rho=rho_prior, Q=0.3)\n",
" noise_priors = [jitter_prior, sigma_prior, rho_prior]\n",
"\n",
" ## define the lightcurve model mu(t;paramters)\n",
" orbit = xo.orbits.KeplerianOrbit(\n",
" period=p_priors,\n",
" t0=tmin_priors,\n",
" b=b_priors,\n",
" duration=d_priors,\n",
" ror=r_priors,\n",
" )\n",
" star = xo.LimbDarkLightCurve(u_prior)\n",
" lightcurve_models = star.get_light_curve(orbit=orbit, r=r_priors, t=t)\n",
" lightcurve = 1e3 * pm.math.sum(lightcurve_models, axis=-1) + f0_prior\n",
" my_planet_transit_model.lightcurve_models = lightcurve_models\n",
" rho_circ = pm.Deterministic(name=RHO_CIRC, var=orbit.rho_star)\n",
"\n",
" # Finally the GP likelihood\n",
" residual = y - lightcurve\n",
" gp_kwargs = dict(diag=yerr**2 + jitter_prior**2, quiet=True)\n",
" gp = GaussianProcess(kernel, t, **gp_kwargs)\n",
" gp.marginal(name=\"obs\", observed=residual)\n",
" my_planet_transit_model.gp_mu = gp.predict(residual, return_var=False)\n",
"\n",
" # cache params\n",
" my_params = dict(\n",
" planet_params=planet_priors,\n",
" noise_params=noise_priors,\n",
" stellar_params=stellar_priors,\n",
" period_params=p_params,\n",
" )\n",
" return my_planet_transit_model, my_params"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "14030d1b",
"metadata": {
"lines_to_end_of_cell_marker": 2,
"ploomber": {
"memory_usage": 263.078125,
"timestamp_end": 1691727606.960698,
"timestamp_start": 1691727399.012704
},
"pycharm": {
"name": "#%%\n"
},
"tags": [
"exe"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 3min 27s (started: 2023-08-11 14:16:39 +10:00)\n"
]
}
],
"source": [
"planet_transit_model, params = build_planet_transit_model(tic_entry)\n",
"model_varnames = get_untransformed_varnames(planet_transit_model)\n",
"test_model(planet_transit_model)"
]
},
{
"cell_type": "markdown",
"id": "2dd3941a",
"metadata": {},
"source": [
"### Optimizing the initial point for sampling\n",
"We help out the sampler we try to find an optimized set of initial parameters to begin sampling from."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "b3bd6f74",
"metadata": {
"ploomber": {
"memory_usage": 263.078125,
"timestamp_end": 1691727607.044851,
"timestamp_start": 1691727606.961898
},
"tags": [
"exe",
"remove-output"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 1.01 ms (started: 2023-08-11 14:20:07 +10:00)\n"
]
}
],
"source": [
"if tic_entry.optimized_params is None:\n",
" init_params = get_optimized_init_params(planet_transit_model, **params)\n",
" tic_entry.save_data(optimized_params=init_params)\n",
"else:\n",
" init_params = tic_entry.optimized_params.to_dict()"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "8e6531eb",
"metadata": {
"ploomber": {
"memory_usage": 274.0,
"timestamp_end": 1691727609.164817,
"timestamp_start": 1691727607.045324
},
"tags": [
"exe"
]
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"time: 2.04 s (started: 2023-08-11 14:20:07 +10:00)\n"
]
},
{
"data": {
"text/html": [
"