From c075016bef6dec6ec46213e1dfbb56253da2dfae Mon Sep 17 00:00:00 2001 From: Matt Craig Date: Mon, 25 Nov 2024 10:23:24 -0600 Subject: [PATCH] Add draft of next notebook --- .../06b-transit-fit-template-fancy-plot.ipynb | 420 ++++++++++++++++++ 1 file changed, 420 insertions(+) create mode 100644 stellarphot/notebooks/photometry/06b-transit-fit-template-fancy-plot.ipynb diff --git a/stellarphot/notebooks/photometry/06b-transit-fit-template-fancy-plot.ipynb b/stellarphot/notebooks/photometry/06b-transit-fit-template-fancy-plot.ipynb new file mode 100644 index 00000000..527778fb --- /dev/null +++ b/stellarphot/notebooks/photometry/06b-transit-fit-template-fancy-plot.ipynb @@ -0,0 +1,420 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from itertools import product\n", + "\n", + "import numpy as np\n", + "\n", + "from astropy.timeseries import TimeSeries, aggregate_downsample\n", + "from astropy.time import Time\n", + "from astropy.table import Table, Column\n", + "from astropy import units as u\n", + "from matplotlib import pyplot as plt\n", + "\n", + "from stellarphot.transit_fitting import TransitModelFit, TransitModelOptions\n", + "from stellarphot.io import TOI\n", + "from stellarphot.plotting import plot_transit_lightcurve\n", + "from stellarphot.gui_tools.photometry_widget_functions import TessAnalysisInputControls, filter_by_dates\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 0. Get some data\n", + "\n", + "+ Select photometry file with relative flux\n", + "+ Select passband\n", + "+ Select TESS info file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "taic = TessAnalysisInputControls()\n", + "taic" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# 👉 File with photometry, including flux\n", + "photometry_file = taic.photometry_data_file\n", + "inp_photometry = taic.phot_data\n", + "\n", + "# 👉 File with exoplanet info in\n", + "tess_info_output_file = taic.tic_info_file\n", + "tess_info = TOI.model_validate_json(tess_info_output_file.read_text())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Get just the target star and some information about it" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "photometry = inp_photometry.lightcurve_for(1, flux_column=\"relative_flux\", passband=taic.passband).remove_nans()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### You may need to alter some of the settings here" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Fit settings\n", + "\n", + "+ Do any detrending by a covariate?\n", + "+ Which parameters are fixed?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# These affect the fitting that is done\n", + "\n", + "model_options = TransitModelOptions()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Find the OOT region and use it to get normalization factor" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "that_transit = tess_info.transit_time_for_observation(photometry.time)\n", + "start = that_transit - tess_info.duration / 2\n", + "mid = that_transit\n", + "end = that_transit + tess_info.duration / 2\n", + "\n", + "after_transit = (photometry[\"bjd\"] - 2400000 * u.day) > end\n", + "\n", + "outside_transit = (photometry[\"bjd\"] < start) | (photometry[\"bjd\"] > end)\n", + "\n", + "normalization_factor = np.nanmean(1 / photometry[\"relative_flux\"][outside_transit])\n", + "normalized_flux = Column(photometry[\"relative_flux\"] * normalization_factor, name=\"normalized_flux\")\n", + "norm_flux_error = Column(normalization_factor * photometry[\"relative_flux_error\"].value, name=\"normalized_flux_error\")\n", + "photometry.add_columns([normalized_flux, norm_flux_error])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Bin Data\n", + "\n", + "Need\n", + "* data table\n", + "* start\n", + "* end\n", + "* bin_size\n", + "\n", + "Data is binned twice because one finds means and the other errors\n", + "\n", + "**Also make times smaller**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "t_ob = Time(photometry[\"bjd\"], scale=\"tdb\", format=\"jd\")\n", + "ts = TimeSeries(\n", + " [\n", + " photometry[\"normalized_flux\"],\n", + " photometry[\"airmass\"],\n", + " photometry[\"xcenter\"],\n", + " photometry[\"sky_per_pix_avg\"],\n", + " photometry[\"width\"],\n", + " ],\n", + " time=t_ob,\n", + ")\n", + "ts2 = TimeSeries(\n", + " [Column(\n", + " data=photometry[\"normalized_flux_error\"],\n", + " name=\"normalized_flux_error\"\n", + " )],\n", + " time=t_ob\n", + ")\n", + "\n", + "first_time = photometry[\"bjd\"][0] - 2400000\n", + "last_time = photometry[\"bjd\"][-1] - 2400000\n", + "\n", + "\n", + "def add_quad(x):\n", + " try:\n", + " n = len(x)\n", + " except TypeError:\n", + " n = 1\n", + " return np.sqrt(np.nansum(x**2)) / n\n", + "\n", + "\n", + "binned = aggregate_downsample(ts, time_bin_size=model_options.bin_size * u.min)\n", + "binned2 = aggregate_downsample(ts2, time_bin_size=model_options.bin_size * u.min, aggregate_func=add_quad)\n", + "\n", + "binned[\"normalized_flux_error\"] = binned2[\"normalized_flux_error\"]\n", + "# binned_time = BinnedTimeSeries(photometry['bjd'], time_bin_start=first_time, time_bin_end=last_time, time_bin_size=bin_size)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model, fit, plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create the model " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Make the model\n", + "mod = TransitModelFit()\n", + "\n", + "# Setup the model\n", + "mod.setup_model(\n", + " binned_data=binned,\n", + " t0=mid.jd - 2400000, # midpoint, BJD\n", + " depth=tess_info.depth_ppt, # parts per thousand\n", + " duration=tess_info.duration.to(\"day\").value, # days\n", + " period=tess_info.period.to(\"day\").value, # days\n", + " model_options=model_options,\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Run the fit" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mod.fit()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Look at the results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plt.plot(mod.times, mod.data, \".\")\n", + "plt.plot(mod.times, mod.model_light_curve())\n", + "plt.vlines(start.jd - 2400000, 0.98, 1.02, colors=\"r\", linestyle=\"--\", alpha=0.5)\n", + "plt.vlines(end.jd - 2400000, 0.98, 1.02, colors=\"r\", linestyle=\"--\", alpha=0.5)\n", + "plt.title(\"Data and fit\")\n", + "plt.grid()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exclude data by date *if needed*\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bad_time = filter_by_dates(\n", + " phot_times=photometry[\"bjd\"],\n", + " use_no_data_before=Time(2400000, format=\"jd\", scale=\"tdb\"),\n", + " use_no_data_between=[\n", + " [\n", + " Time(2400000, format=\"jd\", scale=\"tdb\"),\n", + " Time(2400000, format=\"jd\", scale=\"tdb\"),\n", + " ]\n", + " ],\n", + " use_no_data_after=Time(2499999, format=\"jd\", scale=\"tdb\"),\n", + ")\n", + "\n", + "photometry = photometry[~bad_time]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mod.model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Make the big plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "plot_transit_lightcurve(\n", + " photometry,\n", + " mod,\n", + " tess_info,\n", + " model_options.bin_size * u.min\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Attempt to calculate BIC, but...this seems to have side effects on the rest of notebook " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def evaluate_fits(mod):\n", + " BICs = []\n", + " settings = []\n", + " all_trendable = mod._all_detrend_params\n", + " tf_sequence = product([True, False], repeat=len(all_trendable))\n", + " for fixed in tf_sequence:\n", + " this_summary = []\n", + " for param, fix in zip(all_trendable, fixed):\n", + " trend_mod = getattr(mod.model, f\"{param}_trend\")\n", + " if fix:\n", + " setattr(mod.model, f\"{param}_trend\", 0.0)\n", + " trend_mod.fixed = fix\n", + " this_summary.append(f\"{param}: {not fix}\")\n", + "\n", + " settings.append(\", \".join(this_summary))\n", + " mod.fit()\n", + " BICs.append(mod.BIC)\n", + " return Table(data=[settings, BICs], names=[\"Fit this param?\", \"BIC\"])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "bic_table = evaluate_fits(mod)\n", + "bic_table.sort(\"BIC\")\n", + "bic_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "stelldev-pyd2", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}