{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# 🚀 Ariel Exoplanet Spectra Challenge – Quick Guide  \n\n## 1 The big picture 🌌  \n\n* **Mission**   The European Space Agency’s **Ariel** telescope (launch ≈ 2029) will observe ≈ 1 000 exoplanets as they **transit** their host stars.  \n\n### The science behind transmission spectroscopy\n\nWhen an **exoplanet transits** its star we normally see a tiny dip in total brightness, but a **sliver of starlight also skims through the planet’s atmosphere** on the way to us.  \nMolecules in that gas absorb light at their own fingerprint wavelengths, so the planet looks ever-so-slightly larger at those colours.\n\nComparing the star’s spectrum *in transit* vs *out of transit* yields a **transmission spectrum** that reveals which molecules are present.\n\n---\n\n### 1.2 What is “Transit Depth” and a Transmission Spectrum?\n\nWhen an exoplanet passes in front of its host star (a **transit**), it blocks a small fraction of the star’s light. The **transit depth** is simply the relative drop in stellar flux at mid‐transit:\n\n$$\n\\delta = \\frac{\\Delta F}{F_\\star} \\approx \\left(\\frac{R_p}{R_\\star}\\right)^2\n$$\n\nwhere  \n- $F_\\star$ is the out‐of‐transit stellar flux,  \n- $\\Delta F$ is the decrease in flux at mid‐transit,  \n- $R_p$ is the planet radius, and  \n- $R_\\star$ is the stellar radius.  \n\n*(See Perryman 2018, Chap. 6 §6.12 for a detailed derivation.)*\n\n---\n\n### 1.3 Transmission Spectrum\n\nA **transmission spectrum** measures how the transit depth varies with wavelength, \\(\\delta(\\lambda)\\).  Different molecules in the planet’s atmosphere absorb starlight more strongly at certain wavelengths, so:\n\n- **At wavelengths where the atmosphere is opaque**, the effective \\(R_p(\\lambda)\\) is larger ⇒ **deeper** transit.\n- **At transparent wavelengths**, the atmosphere contributes little ⇒ **shallower** transit.\n\nPlotting \\(\\delta(\\lambda)\\) thus reveals absorption features of molecules like H₂O, CH₄, CO₂, etc.  By fitting these features, we can infer atmospheric composition, temperature, and clouds.  \n\n*For an accessible review of transit spectroscopy, see Winn (2010; _Exoplanet transits and atmospheres_, arXiv:1001.2010) or Perryman (2018, Chap. 11 §11.6).*\n\n\n### 1.4;· Geometry & signal size\n* Only a few **scale heights** of atmosphere back-light the star:  \n  $ H = \\dfrac{kT}{\\mu g}$.\n* Extra area blocked ≈ $ 2\\pi R_p H $.  \n  For a hot Jupiter ( $ R_p ≈ 1\\,R_\\text{Jup}$, $ H ≈ 1000 $ km) the added dimming is **≈ 100–300 ppm**, demanding highly stable detectors.\n\n\n- But what is scale height (H)??\n\n| Symbol | Meaning | Typical hot-Jupiter value | Typical Earth value |\n|--------|---------|---------------------------|---------------------|\n| $k$  | Boltzmann constant $ 1.38 \\times 10^{-23}\\,\\text{J K}^{-1}$ | — | — |\n| $T$  | Atmospheric temperature | $ 1500\\text{–}2500\\ \\text{K}$ (irradiated) | $ 288\\ \\text{K}$ |\n| $\\mu$| Mean molecular mass | $ \\approx 2.3\\,m_{\\text{p}}$ for H/He gas | $\\approx 29\\,m_{\\text{p}}$ ($N_2$+ $O_2$) |\n| $g$  | Surface gravity | $ 10\\text{–}30\\ \\text{m s}^{-2}$ | $ 9.8\\ \\text{m s}^{-2}$ |\n\n* **Definition:** the vertical distance over which pressure (and density) fall by $ e \\approx 2.718 $. Hot, low-gravity planets therefore have thick, “puffed-up” atmospheres\n* **Why “puffed-up” atmospheres?** High \\(T\\) and low $\\mu$ increase \\(H\\); low gravity $g$ further inflates it.  \n  *Hot Jupiters easily reach $ H \\sim 500\\text{–}1500 $ km — orders of magnitude larger than Earth’s \\(H $\\approx$ 8\\) km.*\n\n\n- **Extra geometric area blocked**\n\nWe start with the opaque planetary disk of radius $R_p$.  \nAdd an atmospheric shell of thickness \\(H\\):\n\n$$\nA_{\\text{full}} = \\pi \\bigl(R_p + H\\bigr)^{2}\n                 = \\pi R_p^{2} + 2\\pi R_p H + \\pi H^{2}.\n$$\n\nBecause $ H \\ll R_p$ (e.g. $ H/R_p \\sim 10^{-2}$), the $ H^{2}$ term is negligible, so\n\n$$\n\\boxed{\\Delta A \\;\\approx\\; 2\\pi R_p H}.\n$$\n\n> **Intuition:** the atmosphere acts like a thin ring whose width is the scale height and whose circumference is $ 2\\pi R_p$.\n\n\n- Converting area to observable “transit depth”\n\nThe observable quantity is the *fractional* stellar flux drop,\n\n$$\n\\delta = \\frac{A}{A_\\star}\n        = \\frac{\\pi R_p^{2}}{\\pi R_\\star^{2}}\n        = \\left(\\frac{R_p}{R_\\star}\\right)^{2}.\n$$\n\nThe **atmospheric contribution** is then\n\n$$\n\\Delta\\delta\n  = \\frac{\\Delta A}{\\pi R_\\star^{2}}\n  \\approx \\frac{2 R_p H}{R_\\star^{2}}.\n$$\n\n\n- **Example calculation — *canonical hot Jupiter***\n\n| Parameter | Value | Comment |\n|-----------|-------|---------|\n| $R_p$     | $1.0\\,R_{\\text{Jup}} = 7.15\\times10^{7}\\ \\text{m}$ | Well-studied targets (HD 209458 b, WASP-39 b, …) |\n| $H$       | $1000\\ \\text{km} = 1.0\\times10^{6}\\ \\text{m}$ | For $T\\simeq1500\\ \\text{K},\\; g\\simeq15\\ \\text{m\\,s}^{-2}$ |\n| $R_\\star$ | $1.0\\,R_{\\odot} = 6.96\\times10^{8}\\ \\text{m}$ | Sun-like host |\n\n$$\n\\Delta\\delta\n  \\approx\n  \\frac{2 \\times 7.15\\times10^{7}\\ \\text{m} \\times 1.0\\times10^{6}\\ \\text{m}}\n       {(6.96\\times10^{8}\\ \\text{m})^{2}}\n  \\approx\n  1.5\\times10^{-4}\n  = 150\\ \\text{ppm}.\n$$\n\n* **ppm = parts per million** → $150\\,\\text{ppm} = 0.015\\%$ of the star’s light.  \n* Such minute signals require ultra-stable space telescopes (HST, JWST, CHEOPS, PLATO, Ariel).\n\n> **Rule of thumb:** each scale height adds roughly $(2 R_p H / R_\\star^{2})$ to the transit depth.  \n> Measuring how $\\Delta\\delta$ varies with wavelength—because molecular absorption changes the effective $H$—is the essence of transmission spectroscopy.\n\n---\n\n### 1.5 Nuances & caveats\n\n* **Multiple scale heights:** we may detect absorption from several scale heights, boosting $\\Delta\\delta$.  \n* **Star size:** smaller stars (M-dwarfs) make $\\Delta\\delta$ larger because $R_\\star$ is smaller.  \n* **Clouds/hazes:** high-altitude aerosols can truncate the effective height, flattening the spectrum.  \n* **Instrument stability:** we need relative precision of $\\sim10$ ppm between wavelengths—hence sophisticated calibration and systematics removal.\n\n---\n\n#### 1.6 Molecular fingerprints\n| Species | Strong bands (µm) | Notes |\n|---------|-------------------|-------|\n| H₂O     | 1.4, 1.9          | Dominant in many hot atmospheres |\n| CO₂     | 4.3               | First seen by JWST in WASP-39 b |\n| Na, K   | 0.59, 0.77        | Narrow optical lines |\n| CO      | 2.3, 4.6          | Traced with high-res ground spectra |\n\nEach species’ unique pattern is matched to databases such as **HITRAN** to retrieve composition, temperature and pressure.\n\n---\n\n#### 1.7 How we observe it\n\n| Facility | Wavelength range | Milestone results |\n|----------|------------------|-------------------|\n| **HST / WFC3** | 1.1–1.7 µm | First robust H₂O detections (HD 209458 b) |\n| **Ground** (e.g. VLT/ESPRESSO) | 0.4–2.5 µm (high-res) | Resolves individual CO, H₂O lines; measures winds |\n| **JWST** (NIRISS/NIRSpec/MIRI) | 0.6–12 µm | First CO₂ detection (WASP-39 b, 2022) |\n\nJWST’s broad simultaneous coverage lets us see multiple molecules and clouds in one pass.\n\n---\n\n#### 1.8 What we learn\n* **Elemental ratios** (e.g. C/O) constrain planet-formation zones.\n* **Clouds & hazes** flatten spectra; altitude and particle size affect albedo.\n* **Temperature–pressure profiles** inform circulation and photochemistry.\n\n---\n\n#### 1.9 Key challenges\n* **Stellar activity** (spots, faculae) can mimic or mask signals.  \n* **Clouds/hazes** may obscure molecular lines → need multi-band data.  \n* **Telluric absorption** for ground telescopes requires meticulous correction.\n\nWith larger telescopes and ultra-stable instruments, transmission spectroscopy is moving from hot Jupiters to **smaller, potentially habitable worlds**, giving us a first remote sniff of alien atmospheres.\n\n\n* **The task**   Given extremely noisy light‑curves that mimic Ariel observations, **reconstruct the planet’s transmission spectrum** (how much starlight is blocked at each wavelength) **and say how uncertain we are**.\n\n---\n","metadata":{}},{"cell_type":"markdown","source":"\n## 2 What is inside the dataset 📦  \n\n| Folder / file | What it contains | Shape / dtype (typical) | Purpose in analysis |\n|---------------|-----------------|-------------------------|---------------------|\n| `AIRS-CH0_signal_*.parquet` | Near‑infra‑red detector images for each *visit* | 11 250 frames × 32 px × 356 px (`uint16`) | Holds the wavelength‑dependent signal (1.95–3.90 µm). |\n| `FGS1_signal_*.parquet` | Optical photometry frames for the same visit | 135 000 frames × 32 px × 32 px (`uint16`) | Broadband “white‑light” reference (0.6–0.8 µm). |\n| `dark / flat / dead / linear_corr / read` | Calibration maps | various | Used to correct detector artefacts before extraction. |\n| `train.csv` | **Ground‑truth spectrum** (only for training planets) | #planets × 283 (`float32`) | Supervision target during model training. |\n| `wavelengths.csv` | Central λ for each spectral bin (µm) | 283 (`float64`) | Axis for the spectrum. |\n| `*_star_info.csv` | Stellar & orbital parameters (R$\\odot$, M$\\odot$, Ts, P, inclination…) | #planets × 9 | Helpful priors for transit modelling. |\n| `adc_info.csv` | 2 numbers per instrument: **gain** & **offset** | 2 rows | Restores full dynamic range: `raw*gain + offset`. |\n| `sample_submission.csv` | Template for predictions | #test_planets × (1 + 2×283) | Must match when you upload results. |\n\n🔑 **Key fact** – Each planet can have ***multiple visits***, so you may see several `*_signal_1.parquet`, `*_signal_2.parquet`, … files in the same folder.\n\n---\n\n## 3 What should be the output 🧮  \n\nFor every **test‑set planet** and for **each of the 283 wavelength bins** you have to predict two numbers:\n\n* **μ<sub>i</sub>** – your best estimate of the true transit depth (in parts‑per‑million, ppm).  \n* **σ<sub>i</sub>** – a \\(1σ\\) uncertainty that says *“I believe the real answer lies within μ±σ with ~68 % probability.”*\n\nThe submission file therefore has:\n\n| Column | Meaning |\n|--------|---------|\n| `planet_id` | unique identifier |\n| `mu_1 … mu_283` | mean spectrum values |\n| `sigma_1 … sigma_283` | corresponding uncertainties |\n\nTotal columns = **1 + 283 + 283 = 567**.\n\n---\n\n## 4 How your predictions are scored 🔍  \n\n1. **Gaussian log‑likelihood** (GLL) is computed for every wavelength point **i** and planet **p**:\n\n$$\n\\text{GLL}_{p,i}\n= -\\tfrac12 \\Bigl[\n\\ln\\!\\bigl(2\\pi\\sigma_{p,i}^2\\bigr)\n+ \\frac{\\bigl(y_{p,i}-\\mu_{p,i}\\bigr)^2}{\\sigma_{p,i}^2}\n\\Bigr]\n$$\n\n* $ y_{p,i}$ = hidden ground‑truth depth  \n* $\\mu_{p,i}$ = our prediction  \n* $\\sigma_{p,i}$ = our uncertainty  \n\n2.   Sum over wavelengths **and** over all test planets to get a total **L**.\n\n3.   Convert **L** to a [0, 1] score so that *1 = perfection* and *0 = very poor*:\n\n$$\n\\text{score}\n= \\frac{L_{\\text{ref}} - L}{L_{\\text{ref}} - L_{\\text{ideal}}}\n$$\n\n* $\\mathbf{L_{\\text{ideal}}}$ – GLL achieved by a *perfect* model that recovers the exact truth with an unrealistically small uncertainty (10 ppm for AIRS, 1 ppm for FGS1).  \n* $\\mathbf{L_{\\text{ref}}}$ – GLL of a lazy baseline that always predicts the **mean** training spectrum with its **variance**.\n\n4.   **Instrument weighting** – Every AIRS wavelength bin contributes\n\n$\nw_{\\text{AIRS}} \\approx \\frac{1.95}{282} \\; \\approx \\; 0.0069,\n$\n\nwhile the single FGS1 broadband point is given **double extra weight**\n\n$\nw_{\\text{FGS1}} = 2 \\times 0.2 = 0.4.\n$\n\nThese factors are applied before summing into **L**.","metadata":{}},{"cell_type":"markdown","source":"## 5. Lets observe the data","metadata":{}},{"cell_type":"code","source":"!pip install -q batman-package","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T07:36:31.204191Z","iopub.execute_input":"2025-08-03T07:36:31.206252Z","iopub.status.idle":"2025-08-03T07:39:55.682602Z","shell.execute_reply.started":"2025-08-03T07:36:31.206171Z","shell.execute_reply":"2025-08-03T07:39:55.68126Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import random\nimport time\nimport logging\nfrom tqdm import tqdm\nfrom pathlib import Path\nimport pathlib, json, gc, warnings\nimport numpy as np, pandas as pd\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nfrom scipy.ndimage import median_filter\nwarnings.filterwarnings(\"ignore\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:06.533326Z","iopub.execute_input":"2025-07-22T12:26:06.533684Z","iopub.status.idle":"2025-07-22T12:26:09.957088Z","shell.execute_reply.started":"2025-07-22T12:26:06.533649Z","shell.execute_reply":"2025-07-22T12:26:09.956186Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 🔍 point this path to the folder that contains 'train.csv', etc.\nDATA_DIR = pathlib.Path(\"/kaggle/input/ariel-data-challenge-2025/\")\nWORK_DIR    = pathlib.Path(\"/kaggle/working\")\nassert DATA_DIR.exists()\n\n\nTRAIN_ROOT  = DATA_DIR/\"train\"\nTEST_ROOT   = DATA_DIR/\"test\"\n\nSAMPLE_COLS = pd.read_csv(DATA_DIR/\"sample_submission.csv\", nrows=0).columns\n\n\nwav_df = pd.read_csv(DATA_DIR / \"wavelengths.csv\", header=None)\nwav = wav_df.iloc[1].astype(float).values   # (283,)\nadc = pd.read_csv(DATA_DIR / \"adc_info.csv\", index_col=0)\n\ndisplay(adc)\n\ntrain_truth = (\n    pd.read_csv(DATA_DIR / \"train.csv\")\n      .set_index(\"planet_id\")\n      .astype(float)              # every cell is now a plain Python float\n)\n\nn_planets    = len(train_truth)\nprint(f\"Number of training planets: {n_planets}\")\n\n\nstar_info = pd.read_csv(DATA_DIR/\"train_star_info.csv\")\ntest_star_info   = pd.read_csv(DATA_DIR/\"test_star_info.csv\")\nstar_info.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:09.958783Z","iopub.execute_input":"2025-07-22T12:26:09.959332Z","iopub.status.idle":"2025-07-22T12:26:10.273616Z","shell.execute_reply.started":"2025-07-22T12:26:09.959305Z","shell.execute_reply":"2025-07-22T12:26:10.272795Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_planets = star_info[\"planet_id\"].astype(str).tolist()\n\n# ─── 1) Seed the RNG for reproducibility ────────────────────────────────────\nrandom.seed(42)\n\n# ─── 2) Pick one planet (always the same) ──────────────────────────────────\npid = random.choice(train_planets)\n# Or, to avoid randomness altogether, just do:\n# pid = train_planets[0]\n\nprint(\"chosen planet :\", pid)\n\n# ─── 3) Point at its data folder and list the first few files ─────────────\nplanet_dir = DATA_DIR / \"train\" / pid\nprint(list(planet_dir.glob(\"*\"))[:10])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:10.274436Z","iopub.execute_input":"2025-07-22T12:26:10.274775Z","iopub.status.idle":"2025-07-22T12:26:10.291767Z","shell.execute_reply.started":"2025-07-22T12:26:10.274745Z","shell.execute_reply":"2025-07-22T12:26:10.290894Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 6. Zooming into one **signal** Parquet file 🔬  \n\n#### What the numbers are  \n\n* Each **row** = a *detector frame* (one read of the sensor).  \n* Each **column** = one *pixel* of that frame, flattened.  \n* **Data type** = `uint16` → raw analogue-to-digital counts.  \n* **Conversion** to physical units (photo-electrons) requires  \n\n$$\n\\text{electrons} = \\text{raw\\_count} \\times \\text{gain} + \\text{offset},\n$$\n\nwhere the gain & offset come from **`adc_info.csv`**.\n","metadata":{}},{"cell_type":"code","source":"# The *row label* is the FGS1 offset; the real numbers sit in row 0\ngain_fgs   = float(adc.iloc[0][\"FGS1_adc_gain\"])\noffset_fgs = float(adc.index[0])                      # index label itself\ngain_air   = float(adc.iloc[0][\"AIRS-CH0_adc_gain\"])\noffset_air = float(adc.iloc[0][\"AIRS-CH0_adc_offset\"])\n\nprint(\"FGS1  gain / offset :\", gain_fgs, offset_fgs)\nprint(\"AIRS  gain / offset :\", gain_air, offset_air)\n\n# ───────────────────────────────────────────────────────────────\n# AIRS-CH0  (near-infrared spectrograph)\n# ───────────────────────────────────────────────────────────────\nairs_path = next(planet_dir.glob(\"AIRS-CH0_signal_*.parquet\"))\nairs_raw  = pd.read_parquet(airs_path).values.astype(\"float32\")   # (11250, 11392)\n\nairs_phys = airs_raw * gain_air + offset_air                      # counts → electrons\ncube_air  = airs_phys.reshape(-1, 32, 356)                        # (time, row, λ)\n\nprint(\"AIRS cube:\", cube_air.shape)                               # sanity: (11250, 32, 356)\n\n# ───────────────────────────────────────────────────────────────\n# FGS1  (optical broadband photometer)\n# ───────────────────────────────────────────────────────────────\nfgs_path = next(planet_dir.glob(\"FGS1_signal_*.parquet\"))\nfgs_raw  = pd.read_parquet(fgs_path).values.astype(\"float32\")     # (135000, 1024)\n\nfgs_phys = fgs_raw * gain_fgs + offset_fgs\ncube_fgs = fgs_phys.reshape(-1, 32, 32)                           # (time, row, col)\n\nprint(\"FGS cube :\", cube_fgs.shape)                               # sanity: (135000, 32, 32)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:10.293719Z","iopub.execute_input":"2025-07-22T12:26:10.29398Z","iopub.status.idle":"2025-07-22T12:26:15.3282Z","shell.execute_reply.started":"2025-07-22T12:26:10.293959Z","shell.execute_reply":"2025-07-22T12:26:15.327354Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pyarrow.parquet as pq\n\ndef first_file(glob_pattern):\n    \"\"\"Return the first file that matches the given pattern.\"\"\"\n    files = list(DATA_DIR.glob(glob_pattern))\n    if not files:\n        raise FileNotFoundError(f\"No files match {glob_pattern}\")\n    return files[0]\n\n# Grab one AIRS and one FGS Parquet file (any planet, any visit)\nairs_file = first_file(\"train/*/AIRS-CH0_signal_*.parquet\")\nfgs_file  = first_file(\"train/*/FGS1_signal_*.parquet\")\n\ndef parquet_shape(path):\n    \"\"\"Return (#rows, #columns) without loading the whole file.\"\"\"\n    pf = pq.ParquetFile(path)\n    return pf.metadata.num_rows, pf.metadata.num_columns\n\nairs_rows, airs_cols = parquet_shape(airs_file)\nfgs_rows,  fgs_cols  = parquet_shape(fgs_file)\n\n# Derive the 2-D geometry that the columns flatten into\nairs_geometry = f\"32 × {airs_cols // 32}\"           # 32 spatial rows\nfgs_side      = int(fgs_cols ** 0.5)                # expect 32 for 1 024\nfgs_geometry  = f\"{fgs_side} × {fgs_side}\"\n\nsummary = pd.DataFrame({\n    \"Instrument\":          [\"AIRS-CH0\", \"FGS1\"],\n    \"Frames (rows)\":       [airs_rows,   fgs_rows],\n    \"Flattened columns\":   [airs_cols,   fgs_cols],\n    \"Un-flattened geometry\": [airs_geometry, fgs_geometry],\n    \"Cadence\":             [\"0.8 s\", \"0.1 s\"]  # fixed by instrument design\n})\n\nsummary","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:15.32924Z","iopub.execute_input":"2025-07-22T12:26:15.329481Z","iopub.status.idle":"2025-07-22T12:26:23.467Z","shell.execute_reply.started":"2025-07-22T12:26:15.329462Z","shell.execute_reply":"2025-07-22T12:26:23.466Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 7 What do these pixels *mean* physically? 🌈  \n\n### AIRS-CH0  \n\n* The **columns** run along the **dispersion axis** → a pixel column ~ one wavelength bin.  \n* The **rows** are the **spatial axis** (height of the spectral trace).  \n* A typical reduction step: **optimal extraction** collapses rows with weights → a 1-D spectrum per frame.\n\n### FGS1  \n\n* No dispersion: it’s an “imaging photometer”.  \n* Summing all 32 × 32 pixels gives a **broad-band light curve** used for systematics correction.","metadata":{}},{"cell_type":"markdown","source":"## 8 Making sense mathematically 📐  \n\n* A single AIRS pixel time-series $F_{t,r,c}$ is the sum of  \n\n$$\nF = S_{\\star} \\;\\bigl[1 - \\delta(\\lambda_c)\\,\\mathsf{T}(t)\\bigr] + N_\\text{det} + N_\\text{astro},\n$$\n\n\nBelow we expand the compact formula that appeared earlier and show how each\nterm plugs into a practical reduction pipeline.\n\n---\n\n### 8.1 The forward model for one detector pixel\n\nFor a given **time index `t`**, spatial **row `r`** (cross-dispersion), and\nspectral **column `c`** (≈ wavelength bin) the **measured** count in the\nAIRS-CH0 cube is\n\n$$\n\\boxed{\nF_{t,r,c}\n\\;=\\;\nS_\\star \\,\\Bigl[1 - \\delta(\\lambda_c)\\,\\,\\mathsf{T}(t)\\Bigr]\n\\;+\\;\nN_\\text{det}(t,r,c)\n\\;+\\;\nN_\\text{astro}(t)\n}\n\\tag{1}\n$$\n\n| Term  | Description (and how we deal with it) |\n|-------|---------------------------------------|\n| $S_\\star$ | **Stellar baseline** — the mean (out-of-transit) flux; essentially a multiplicative constant. |\n| $\\delta(\\lambda_c)$ | **Transit depth** (what we want) at wavelength \\(λ_c\\). In ideal noiseless data: $\\delta = (R_p/R_\\star)^2$ but in reality it varies with \\(λ\\) because different molecules absorb different bands. |\n| $\\mathsf{T}(t)$ | **Transit shape**: a time-varying function that is 0 out of transit, rises to 1 at mid-transit, and is analytically given by the *Mandel & Agol* (2002) model once you know the orbital parameters. |\n| $N_\\text{det}(t,r,c)$ | **Detector/systematics noise** — thermal drifts, pointing jitter, inter-pixel sensitivity, ramp effect, etc.  Typically *smooth* in time and/or correlated across pixels, so we model it with low-order polynomials, Gaussian Processes, or Pixel-Level-Decorrelation (PLD) basis vectors and **subtract it out**. |\n| $N_\\text{astro}(t)$ | **Astrophysical noise** — e.g. star-spot–driven brightness changes.  Often traced by the broadband FGS1 light-curve or by first principal components of the AIRS cube. |\n\n---\n\n### 8.2 Normalise & linearise the equation\n\nDivide equation (1) by an estimate of the stellar baseline\n$\\hat S_\\star$ (obtained by averaging many out-of-transit frames in the\nbroadband **FGS1** channel or with a polynomial fit):\n\n$$\n\\frac{F_{t,r,c}}{\\hat S_\\star}\n\\;\\approx\\;\n1 - \\delta(\\lambda_c)\\,\\mathsf{T}(t)\n\\;+\\;\n\\underbrace{\\frac{N_\\text{det}}{\\hat S_\\star}}_{\\text{systematics}}\n\\;+\\;\n\\underbrace{\\frac{N_\\text{astro}}{\\hat S_\\star}}_{\\text{stellar var.}}\n$$\n\nThis puts everything on a **dimensionless, near-unity scale** (ppm-level\ndepartures are easier to model numerically).\n\n---\n\n### 8.3 Fitting strategy in practice\n\n1. **Build a design-matrix `X`**  \n   * Columns:  \n     * polynomial/time basis (1, t, t², …),  \n     * roll angle, temperature sensors,  \n     * first _k_ PCA components of the detector cube (PLD),  \n     * normalised FGS1 white-light flux, …  \n   * Goal: capture $N_\\text{det}$ and large-scale $N_\\text{astro}$.\n\n2. **Joint regression**  \n   At each wavelength bin \\(c\\) (or all bins simultaneously in a multi-output\n   GP) solve\n\n   $$\n   y_t \\;=\\; 1 - \\delta_c\\,\\mathsf{T}(t) + X_t\\,\\beta_c + \\epsilon_t,\n   $$\n   where $y_t = F_{t,\\,\\text{collapsed},\\,c}/\\hat S_\\star$.\n\n   * $\\beta_c$ are nuisance coefficients (systematics).  \n   * $\\epsilon_t$ is (often) assumed Gaussian white noise.\n\n3. **Extract $\\delta_c$**  \n   * Point estimate → least-squares or MCMC posterior mean.  \n   * Uncertainty $\\sigma_c$ → covariance of the estimator or full posterior\n     width.\n\n4. **Stack over visits**  \n   If a planet has multiple visits $v=1…V$, repeat\n   Steps 1-3 per visit and combine with a hierarchical model:\n\n   $$\n   \\delta_c\n   \\sim \\mathcal N\\!\\bigl(\\bar\\delta_c,\\; \\sigma_{v,c}^2\\bigr)\n   ,\\qquad\n   \\bar\\delta_c\n   \\sim \\text{prior}.\n   $$\n\n   The posterior mean of $\\bar\\delta_c$ is the **μ** you submit;\n   its posterior σ becomes your **σ** column.\n\n---\n\n### 8.4 Where does the *spectrum* live?\n\nBecause $\\delta(\\lambda_c)$ is *inside* the product with $\\mathsf{T}(t)$,\nall **wavelength channels share the *same* transit shape**\nand differ only in amplitude.  This is what lets us solve for a\n283-element vector of depths from hundreds of thousands of time samples.\n\n* Out-of-transit frames ($\\mathsf{T}=0$) ⇒ just set the baseline.  \n* In-transit frames $(\\mathsf{T}>0$) ⇒ encode the differential depth.\n\nThe job of the reduction pipeline is to **remove everything else** so that the\nonly remaining difference between out- and in-transit parts of the light\ncurve is that tiny amplitude $\\delta(\\lambda_c)$.\n\n---\n\n### 8.5 Connecting to the submission file  \n\nFor every wavelength bin \\(c\\):\n\n| Column in submission | Quantity | Derived from |\n|----------------------|----------|--------------|\n| `mu_c`     | $\\bar\\delta(\\lambda_c)$ | posterior mean of the hierarchical fit |\n| `sigma_c`  | $\\text{StdDev}\\bigl[\\bar\\delta(\\lambda_c)\\bigr]$ | posterior width or propagated uncertainty |\n\nA perfectly calibrated pipeline will yield\n$\\text{(truth)} \\in [\\mu_c - \\sigma_c,\\, \\mu_c + \\sigma_c]$\n~68 % of the time, maximising the Gaussian log-likelihood score.\n\n---\n\n> **Key intuition**:  \n> *Everything* you measure is the star, minus a *tiny* wavelength-dependent dip\n> that slides on and off with the transit shape. Strip away detector quirks and\n> stellar wiggles, and that dip is the planet’s fingerprint.\n","metadata":{}},{"cell_type":"markdown","source":"## 9. Visualization","metadata":{}},{"cell_type":"code","source":"# ---------- ORIGINAL DATA ----------\nwhite_light = cube_fgs.sum(axis=(1, 2))          # integrate all pixels\ntime_s      = np.arange(len(white_light)) * 0.1  # cadence 0.1 s\n\n# ---------- 1) BIN TO 2-SECOND CADENCE ----------\nbin_size   = 20                # 20 × 0.1 s  = 2 s\nn_bins     = len(white_light) // bin_size\nwl_binned  = white_light[: n_bins*bin_size].reshape(n_bins, bin_size).mean(axis=1)\ntime_b     = time_s[: n_bins*bin_size].reshape(n_bins, bin_size).mean(axis=1)\n\n# ---------- 2) RUNNING MEDIAN (±30 s window) ----------\nwl_med = median_filter(wl_binned, size=15)       # 15 × 2 s ≈ 30 s\n\n# ---------- 3) NORMALISE ----------\nwl_norm = wl_binned / wl_binned.mean()\nwl_med  = wl_med    / wl_binned.mean()\n\n# ---------- 4) PLOT ----------\nplt.figure(figsize=(11, 3.5))\n\n# lightly-coloured line for every 2-s bin\nplt.plot(time_b / 60, wl_norm,\n         lw=0.6, alpha=0.3, label=\"2 s bins\")\n\n# bold line for the running median\nplt.plot(time_b / 60, wl_med,\n         lw=1.2, label=\"30 s running median\")\n\n# cosmetic tweaks\nplt.axhline(1.0, color=\"k\", lw=0.8, ls=\"--\", alpha=0.5)\nplt.xlabel(\"Time [minutes]\")\nplt.ylabel(\"Relative flux\")\nplt.title(f\"FGS1 white-light curve  (planet {pid})\")\nplt.ylim(0.985, 1.015)          # tight y-range shows ~1–2 % variations\nplt.legend(frameon=False)\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:23.468003Z","iopub.execute_input":"2025-07-22T12:26:23.468372Z","iopub.status.idle":"2025-07-22T12:26:23.88239Z","shell.execute_reply.started":"2025-07-22T12:26:23.46834Z","shell.execute_reply":"2025-07-22T12:26:23.881482Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 📊 9.1 Interpreting the FGS1 White-Light Curve\n\n| Element on the plot | Physical meaning | Typical next step |\n|---------------------|------------------|-------------------|\n| **Blue mist (“2 s bins”)** | Broadband flux averaged over 2-second windows. Values are *relative* (1 = mean stellar flux). Scatter ≈ photon + detector noise (~300-400 ppm). | Estimate white-noise level; propagate into AIRS channel uncertainties. |\n| **Orange line (“30 s running median”)** | Low-frequency trend (systematics + stellar variability) with high-frequency noise suppressed. | Use as a basis vector or GP component when de-trending AIRS data. |\n| **Dashed line at 1.0** | Reference “no variation” level. | — |\n| **Thin negative spikes** | Cosmic-ray hits or telemetry glitches — sudden frame losses. | Flag & mask those times so they don’t bias fits. |\n| **Broad hump (~100–160 min)** | Classic **detector ramp**: charge-trapping causes counts to creep upward then relax. Amplitude ≈ 0.1 % (1000 ppm). | Fit & divide out (polynomial, exponential, or GP). |\n| **Slow baseline tilt** | Long-term drift / thermal settling. | Include in same de-trending model. |\n\n#### Where is the transit dip?\n* Typical exoplanet depth is **tens–hundreds ppm** → hidden beneath the orange curve at this plot scale.\n* After removing ramp & drift and zooming to ±500 ppm, a symmetric “U-shape” should appear.\n\n#### Why we care about this curve\n1. **Systematics tracer** – Detector & pointing trends are common-mode; the FGS1 channel is a proxy to remove them from AIRS.\n2. **Stellar variability monitor** – Spots/faculae modulate both instruments alike; subtracting this curve helps isolate the planet signal.\n3. **Transit timing anchor** – The white-light dip pins down ingress/egress times that apply to every wavelength bin.\n\n> **Key takeaway:**  \n> The data are dominated by the star’s steady light, plus slow instrumental trends and occasional glitches.  \n> Once those are modelled and removed, the tiny wavelength-dependent transit depth (the planet’s fingerprint) can be measured.","metadata":{}},{"cell_type":"code","source":"# ─── 1) Cache all calibration maps once ──────────────────────────────────────\n_CAL_DIR     = None\n_DARK        = {}\n_FLAT        = {}\n_DEAD_MASK   = {}\n_LIN_COEFFS  = {}\n_READ_NOISE  = {}\n\ndef load_all_calibrations(sample_root: Path, instr=\"AIRS-CH0\"):\n    global _CAL_DIR\n    # pick first planet folder\n    sample_planet = next(p for p in sample_root.iterdir() if p.is_dir())\n    cal_dir = sample_planet / f\"{instr}_calibration_0\"\n    _CAL_DIR = cal_dir\n\n    _DARK[instr]      = pd.read_parquet(cal_dir/\"dark.parquet\").to_numpy(np.float32).reshape(32,356)\n    _FLAT[instr]      = pd.read_parquet(cal_dir/\"flat.parquet\").to_numpy(np.float32).reshape(32,356)\n    _DEAD_MASK[instr] = pd.read_parquet(cal_dir/\"dead.parquet\").to_numpy().astype(bool).reshape(32,356)\n    lin = pd.read_parquet(cal_dir/\"linear_corr.parquet\").to_numpy(np.float32)\n    _LIN_COEFFS[instr] = lin.reshape(-1,32,356)[::-1]    # reverse for Horner’s\n    _READ_NOISE[instr] = pd.read_parquet(cal_dir/\"read.parquet\").to_numpy(np.float32).reshape(32,356)\n\nload_all_calibrations(TRAIN_ROOT, instr=\"AIRS-CH0\")\n\n# ─── 2) Calibrate raw cube ────────────────────────────────────────────────────\ndef calibrate_cube(planet_dir: Path, instr=\"AIRS-CH0\"):\n    \"\"\"\n    1) read raw → apply ADC\n    2) subtract dark, divide flat\n    3) mask dead pixels\n    4) Horner’s nonlinearity correction\n    5) compute variance = shot + read_noise^2\n    returns cube (n_frames,32,356) and var (same)\n    \"\"\"\n    sig = next(planet_dir.glob(f\"{instr}_signal_*.parquet\"))\n    raw = pd.read_parquet(sig).to_numpy(np.float32)\n    arr = raw * gain_air + offset_air\n    cube = arr.reshape(-1,32,356)\n\n    # dark & flat\n    cube -= _DARK[instr][None]\n    safe_flat = np.where(_FLAT[instr]==0, np.nan, _FLAT[instr])\n    cube /= safe_flat[None]\n\n    # mask\n    cube[:, _DEAD_MASK[instr]] = np.nan\n\n    # nonlinearity via Horner’s\n    corr = np.zeros_like(cube)\n    x = cube\n    for coeff in _LIN_COEFFS[instr]:\n        corr = corr * x + coeff[None]\n    cube = corr\n\n    # variance\n    var = np.abs(cube) + _READ_NOISE[instr][None]**2\n\n    return cube, var\n\n# ─── 3) Optimal white-light extraction ───────────────────────────────────────\ndef extract_white_light(cube: np.ndarray, var: np.ndarray):\n    \"\"\"\n    PSF = median over time → normalize → optimal extraction:\n    flux = sum( PSF * cube/var ) / sum( PSF^2/var )\n    sigma = 1/sqrt(denominator)\n    returns time_s (0.8s cadence), flux, sigma\n    \"\"\"\n    psf = np.nanmedian(cube, axis=0)\n    psf /= np.nansum(psf)\n\n    num = np.nansum(psf[None,:,:] * cube / var, axis=(1,2))\n    den = np.nansum(psf[None,:,:]**2 / var, axis=(1,2))\n\n    flux  = num/den\n    sigma = 1.0/np.sqrt(den)\n    time_s = np.arange(cube.shape[0]) * 0.8\n    return time_s, flux, sigma\n\n# ─── 4) Correlated-double sampling & binning ─────────────────────────────────\ndef get_cds(signal: np.ndarray) -> np.ndarray:\n    \"\"\"(end - start) reading pairs\"\"\"\n    return signal[1::2] - signal[::2]\n\ndef bin_obs(cds_signal: np.ndarray, binning: int) -> np.ndarray:\n    \"\"\"sum every `binning` ramps into one frame\"\"\"\n    n = cds_signal.shape[0] // binning\n    b = np.zeros((n, *cds_signal.shape[1:]), dtype=cds_signal.dtype)\n    for i in range(n):\n        b[i] = np.nansum(cds_signal[i*binning:(i+1)*binning], axis=0)\n    return b\n\n# ─── 5) Plot 10 random train white-light curves ──────────────────────────────\ndef plot_sample_train_curves(n_samples=20, binning=30):\n    random.seed(0)\n    planets = random.sample([d for d in TRAIN_ROOT.iterdir() if d.is_dir()], n_samples)\n\n    plt.figure(figsize=(8,5))\n    for pd_dir in tqdm(planets, desc=\"plotting\"):\n        cube, var        = calibrate_cube(pd_dir)\n        time_s, flux_wl, _ = extract_white_light(cube, var)\n\n        # build CDS+binned cube to sum pixels\n        cds   = get_cds(cube)           # (n_ramps,32,356)\n        binned= bin_obs(cds, binning)   # (n_bins,32,356)\n\n        # white-light = sum over all pixels in each binned frame\n        wl_curve = binned.sum(axis=(1,2))\n\n        # normalize to first+last 10% OOT baseline\n        n = len(wl_curve)\n        oot = np.concatenate([wl_curve[:n//10], wl_curve[-n//10:]])\n        baseline = np.nanmedian(oot)\n        wl_norm  = wl_curve / baseline\n\n        frames = np.arange(len(wl_norm))\n        plt.plot(frames, wl_norm, alpha=0.5, lw=1)\n\n    plt.xlabel(\"Time (binned frame index)\")\n    plt.ylabel(\"Normalized flux\")\n    plt.title(f\"White-light curves from {n_samples} random training planets\")\n    plt.ylim(0.9, 1.02)\n    plt.grid(alpha=0.3)\n    plt.tight_layout()\n    plt.show()\n\n# ─── Run ─────────────────────────────────────────────────────────────────────\nplot_sample_train_curves(n_samples=100, binning=30)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:26:23.88349Z","iopub.execute_input":"2025-07-22T12:26:23.883755Z","iopub.status.idle":"2025-07-22T12:42:02.301596Z","shell.execute_reply.started":"2025-07-22T12:26:23.883735Z","shell.execute_reply":"2025-07-22T12:42:02.300693Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### White-light curves from 100 random training planets\n\nThe overplotted curves above show the **normalized white-light flux** extracted from 100 randomly chosen planets in the training set, after full calibration and time-binning. Each colored line is one planet’s transit observation:\n\n- **Horizontal axis**: binned frame index (each bin combines 30 CDS ramps → ∼24 s per point).  \n- **Vertical axis**: normalized total flux in the frame, scaled so that the out-of-transit baseline sits at ≃ 1.0.  \n\n### What we can read from this figure\n\n1. **Transit depths** vary from a few per‐thousand (1.0 – flux dips of ∼0.001) up to several percent (flux dips > 0.02).  \n2. **Transit durations and shapes** differ from object to object: some curves are wide and shallow; others are narrow and deep.  \n3. Because we normalize each curve to its own out-of-transit baseline (using the first / last 10 % of points), all baselines line up at 1.0. The spread you see around 1.0 before and after ingress/egress is residual instrumental noise and systematics.\n\nThis single-panel view gives an at-a-glance sense of the diversity in both **depth** and **duration** across the sample, and confirms that our calibration+binning pipeline is recovering clean, transit-shaped signatures on top of a flat baseline.","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(8,3))\nplt.imshow(cube_air[0], aspect=\"auto\", origin=\"lower\", cmap=\"viridis\")\nplt.colorbar(label=\"e‑ counts\")\nplt.xlabel(\"Dispersion (λ bins)\"); plt.ylabel(\"Spatial rows\")\nplt.title(\"AIRS‑CH0 raw detector frame\")\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:42:02.302875Z","iopub.execute_input":"2025-07-22T12:42:02.303202Z","iopub.status.idle":"2025-07-22T12:42:02.627331Z","shell.execute_reply.started":"2025-07-22T12:42:02.303174Z","shell.execute_reply":"2025-07-22T12:42:02.625979Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 🖼️ What you are looking at – AIRS-CH0 raw detector frame\n\n| Plot feature | Physical meaning |\n|--------------|------------------|\n| **Horizontal bright stripe (rows ≈ 13 – 18)** | The *spectral trace* of the target star.  Each column is a different wavelength bin; brightness encodes the number of photo-electrons captured in that sub-exposure. |\n| **Dispersion axis (x-axis, 0 – 356 px)** | Left → right maps onto **increasing wavelength** (1.95 → 3.90 µm).  The stripe gets noticeably brighter toward the right because the star is redder (more flux at longer λ) and detector QE may rise. |\n| **Spatial axis (y-axis, 0 – 31 px)** | Rows perpendicular to dispersion.  The star’s PSF spans ~4–5 pixels (FWHM) and is centered near row 16. |\n| **Color bar labelled “e-counts”** | Pixel values after ADC conversion (`raw*gain + offset`).  They are **negative** here because the constant electronic *offset* (≈ 1000 e⁻) has been subtracted but the dark current hasn’t been added back yet.  The relative scale still shows where the signal is strongest. |\n| **Speckled background** | Read-noise + dark current + low-level flat-field variations.  It dominates outside the spectral trace. |\n| **Isolated reddish pixels** | Hot pixels or cosmic-ray hits—individual detector defects.  These are masked later. |\n\n#### How to use this frame in the pipeline 🔧  \n1. **Calibration** – subtract dark, divide by flat, interpolate over hot pixels, linearise.  \n2. **Optimal extraction** – for each column, weight rows by the spatial profile (PSF) to sum the stripe into one **flux vs. wavelength** data point.  \n3. **Repeat for every time-step** → build a time–wavelength cube whose rows you saw in the white-light plot earlier.  \n4. **Systematics + transit fitting** – model out the “noise floor” and solve for the tiny wavelength-dependent transit depth.\n\n> **In one sentence:**  \n> The image is a single NIR snapshot where the bright horizontal band is the star’s dispersed light; extracting that band over time is the first step toward recovering the planet’s transmission spectrum.","metadata":{}},{"cell_type":"markdown","source":"### Train‑set labels (ground‑truth spectra) ","metadata":{}},{"cell_type":"code","source":"# ─── Select 3 reproducible planets ───────────────────────────────────────────\nvalid_planets = train_truth.index.intersection(star_info[\"planet_id\"].astype(int))\nrandom.seed(80)\npids = random.sample(list(valid_planets), 3)\nprint(\"Chosen planets:\", pids)\n\n# ─── Prepare molecular features ─────────────────────────────────────────────\nfeatures = {\n    2.7:  \"H₂O\",\n    2.97: \"NH₃\",\n    3.3:  \"CH₄\",\n    3.4:  \"C₂H₂\",\n    3.5:  \"HCN\",\n}\n\n# ─── Plot each spectrum in its own subplot ─────────────────────────────────\nfig, axes = plt.subplots(3, 1, figsize=(10, 12), sharex=True)\n\nfor ax, pi in zip(axes, pids):\n    # Extract AIRS spectrum (skip bin 0)\n    spec      = train_truth.loc[pi].values\n    wav_airs  = wav[1:]\n    depth_ppm = spec[1:] * 1e6\n\n    # Plot spectrum and noise floor\n    ax.plot(wav_airs, depth_ppm, \"-o\", ms=4, lw=1,\n            color=\"navy\", label=\"True spectrum\")\n    ax.fill_between(wav_airs,\n                    depth_ppm - 50,\n                    depth_ppm + 50,\n                    color=\"gray\", alpha=0.1,\n                    label=\"±50 ppm noise\")\n\n    # Formatting\n    ax.set_ylabel(\"Transit depth [ppm]\", fontsize=12)\n    ax.set_title(f\"Planet {pid}\", fontsize=14)\n    ax.grid(alpha=0.3)\n\n    # Mixed transform: x in data, y in axes fraction\n    trans = ax.get_xaxis_transform()\n\n    # Place labels a bit to the right of each dotted line\n    for wl0, label in features.items():\n        if wav_airs.min() < wl0 < wav_airs.max():\n            # draw the dotted line\n            ax.axvline(wl0,\n                       color=\"firebrick\", linestyle=\"--\", alpha=0.7,\n                       label=\"_nolegend_\")\n            # determine text position offset\n            x_text = wl0 + 0.03          # shift right by 0.03 μm\n            # if near right edge, shift left instead\n            if x_text > wav_airs.max():\n                x_text = wl0 - 0.03\n                ha = \"right\"\n            else:\n                ha = \"left\"\n            # place text just inside the top (y=0.98)\n            ax.text(x_text, 0.4,\n                    f\"{label} ({wl0:.2f} µm)\",\n                    color=\"firebrick\",\n                    fontsize=10,\n                    rotation=90,\n                    ha=ha, va=\"top\",\n                    transform=trans,\n                    clip_on=True,\n                    label=\"_nolegend_\")\n\n# ─── Final formatting ───────────────────────────────────────────────────────\naxes[-1].set_xlabel(\"Wavelength [µm]\", fontsize=12)\naxes[0].legend(loc=\"upper right\", fontsize=10)\nplt.tight_layout()\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:42:02.628638Z","iopub.execute_input":"2025-07-22T12:42:02.629482Z","iopub.status.idle":"2025-07-22T12:42:03.471069Z","shell.execute_reply.started":"2025-07-22T12:42:02.629444Z","shell.execute_reply":"2025-07-22T12:42:03.470065Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"### 🖼️  Reading the AIRS-CH0 ground-truth spectrum for planet **1124834224**\n\n| Plot element | What it means | Typical interpretation |\n|--------------|---------------|------------------------|\n| **X-axis (1.95 – 3.90 µm)** | Infra-red wavelengths sampled by Ariel’s AIRS Channel 0. | Each point is a ~7 nm-wide bin (R≈100). |\n| **Y-axis (≈ 4430 – 4500 ppm)** | Extra stellar light blocked by the planet, in **parts-per-million**. (4 500 ppm ≃ 0.45 % of the star’s flux.) | Higher value ⇒ the planet’s apparent radius is larger at that λ. |\n| **Deep trough at ~2.20 µm** | Atmosphere is *transparent* here → smaller apparent radius. | Likely a “window” between strong H<sub>2</sub>O / CO bands. |\n| **Steep rise 2.25 → 2.50 µm** | Rapid increase in opacity. | Water-vapour continuum plus onset of H<sub>2</sub>O band head. |\n| **Broad plateau & local maxima 2.7–3.0 µm** | Strong absorption hump. | Hot-Jupiter spectra often show an H<sub>2</sub>O peak here; cooler planets might show CH<sub>4</sub>. |\n| **Small wiggles atop the plateau** | Real molecular sub-structure (richer line list) **plus** detector/read noise at the ~20 ppm level. | You’ll need to recover these tiny modulations to score well. |\n| **Slow decline beyond 3.2 µm** | Opacity drops; probing deeper, cooler layers. | Can indicate cloud top or decreasing H<sub>2</sub>O cross-section. |\n\n#### Why the exact shape matters\n\n* **Amplitude** (4 430 → 4 500 ppm) encodes **scale height** → temperature, mean molecular weight, gravity.  \n* **Precise wavelength of peaks/troughs** fingerprint specific molecules.  \n* **Slope of the continuum** hints at hazes/clouds and pressure levels.\n\n---\n\n> **In short:** this plot is the planet’s “barcode” in the near-IR.  \n> Your model must predict the same bumps and dips (within well-calibrated error bars) for the unseen test planets to score highly in the competition.","metadata":{}},{"cell_type":"markdown","source":"## 8. A physics based baseline","metadata":{}},{"cell_type":"markdown","source":"## Physics & Math Behind Each Pipeline Step\n\n### 1) ADC Gain & Offset  \nEach pixel records a raw digital count proportional to accumulated charge. We convert to physical electrons by  \n$$\n\\text{electrons} \\;=\\; \\text{raw\\_counts}\\times\\mathrm{gain} \\;+\\;\\mathrm{offset}\n$$  \nrestoring the true signal level before calibration.\n\n### 2) Dark-Current Subtraction  \nEven with the shutter closed, detectors register thermal and bias currents. A “dark frame” $D_{r,c}$ measures this, so we correct:  \n$$\nF'_{t,r,c} \\;=\\; F_{t,r,c} \\;-\\; D_{r,c}\n$$  \nremoving stable bias and thermal noise.\n\n### 3) Flat-Field Correction  \nPixels vary in sensitivity. A “flat field” $F_{f,r,c}$ maps relative pixel gains under uniform illumination. We normalize:  \n$$\nF''_{t,r,c} \\;=\\; \\frac{F'_{t,r,c}}{F_{f,r,c}}\n$$  \nso every pixel responds uniformly to the same incident flux.\n\n### 4) Dead & Hot Pixel Masking  \nSome pixels never respond (“dead”) or saturate (“hot”). We apply a boolean mask $M_{r,c}$, setting those pixels to NaN so they don’t bias later steps.\n\n### 5) Non-linearity Correction (Horner’s Method)  \nAs pixels approach saturation, their response deviates from linear. We have per-pixel polynomial coefficients $a_0,a_1,\\dots,a_5$ and apply Horner’s scheme:  \n$$\nP(x) = a_5x^5 + a_4x^4 + \\dots + a_1x + a_0\n$$  \nimplemented as  \n$$\ny = (((a_5\\cdot x + a_4)\\cdot x + a_3)\\cdot x + \\dots + a_0)\\,,\n$$  \nto linearize each pixel’s output.\n\n### 6) Variance Map (Shot + Read Noise)  \nPhoton arrival is Poisson, so variance $\\approx$ signal. Electronic read noise adds $\\sigma_{\\rm read}^2$. We compute per-pixel variance  \n$$\n\\sigma^2_{t,r,c} = \\bigl|F''_{t,r,c}\\bigr| \\;+\\;\\sigma_{\\rm read,\\;r,c}^2.\n$$\n\n### 7) White-Light Extraction (Optimal PSF-Weighted)  \nTo get a high-SNR broadband flux curve, we weight each pixel by its average brightness pattern (PSF) and inverse variance:  \n1. Compute the median image  \n   $$\n   \\mathrm{PSF}_{r,c} = \\frac{\\mathrm{median}_t\\,F''_{t,r,c}}{\\sum_{r,c}\\mathrm{median}_t\\,F''_{t,r,c}}\n   $$  \n2. Form the weighted sum  \n   $$\n   F_{\\rm wl}(t) \n   = \\frac{\\sum_{r,c}\\mathrm{PSF}_{r,c}\\,F''_{t,r,c}/\\sigma^2_{t,r,c}}\n          {\\sum_{r,c}\\mathrm{PSF}_{r,c}^2/\\sigma^2_{t,r,c}},\n   $$  \n   with uncertainty  \n   $$\n   \\sigma_{\\rm wl}(t)\n   = \\biggl(\\sum_{r,c}\\frac{\\mathrm{PSF}_{r,c}^2}{\\sigma^2_{t,r,c}}\\biggr)^{-1/2}.\n   $$\n\n### 8) Spectral “Box” Extraction  \nTo build a light curve in each wavelength channel, we sum a small band of rows ($\\pm3$ pixels around the PSF center):  \n$$\nS_{t,\\lambda} = \\sum_{r=r_0-3}^{r_0+3} F''_{t,r,\\lambda}.\n$$  \nThis collapses the 2D cube $(r,c)$ into a 2D time–wavelength array.\n\n### 9) Common-Mode Removal  \nInstrumental/stellar systematics often affect all wavelengths simultaneously. We remove the median at each time step:  \n$$\nS'_{t,\\lambda} = \\frac{S_{t,\\lambda}}{\\mathrm{median}_\\lambda\\,S_{t,\\lambda}},\n$$  \nflattening out shared trends and isolating wavelength-dependent transit signals.\n\n### 10) Phase-Folding & Masking  \nKnowing the orbital period $P$ and mid-transit time $t_0$, we compute phase  \n$$\n\\phi = \\frac{(t - t_0)\\bmod P}{P} \\;-\\; 0.5,\\quad \\phi\\in[-0.5,0.5],\n$$  \nand define “in-transit” frames $|\\phi|<\\Delta\\phi$ (e.g.\\ $\\Delta\\phi=0.02$) vs.\\ “out-of-transit.”\n\n### 11) Channel-by-Channel Depth Measurement  \nFor each $\\lambda$-bin we:\n1. Detrend the out-of-transit baseline (e.g.\\ linear fit),  \n2. Normalize and sigma-clip outliers,  \n3. Bin the folded light curve into uniform phase bins,  \n4. Compute the transit depth  \n   $$\n   \\delta(\\lambda) \n   = 1 \\;-\\;\\mathrm{median}\\bigl\\{F_{\\rm binned}(\\phi)\\bigr\\}_{|\\phi|<\\Delta\\phi},\n   $$  \n   and its uncertainty from the scatter within those in-transit bins.\n\n---\n\n**In summary:**  \nWe convert raw detector counts into a fully calibrated, time- and wavelength-dependent flux cube; extract both broadband and spectrally resolved light curves; remove common systematics; fold on the known orbit; and measure the tiny, wavelength-dependent dip $\\delta(\\lambda)$ imprinted by the planet’s atmosphere.  ","metadata":{}},{"cell_type":"code","source":"# ─── 5) Spectral “box” extraction ──────────────────────────────────────────────\ndef box_extract(cube: np.ndarray, center: int = 16, half: int = 3) -> np.ndarray:\n    \"\"\"\n    Sum ±half rows around `center` in the 32×356 cube → shape (n_frames, 356).\n    \"\"\"\n    r0, r1 = center-half, center+half+1\n    return cube[:, r0:r1, :].sum(axis=1)\n\n# ─── 6) Common‐mode removal ────────────────────────────────────────────────────\ndef remove_common_mode(spec_cube: np.ndarray) -> np.ndarray:\n    \"\"\"\n    spec_cube: (n_frames, λ_bins)\n    divide out the common‐mode (per‐frame median) → same shape\n    \"\"\"\n    cm = np.nanmedian(spec_cube, axis=1)\n    return spec_cube / cm[:, None]\n\n# ─── 7) Phase‐fold & masks ─────────────────────────────────────────────────────\ndef phase_fold_masks(time_s: np.ndarray, pid: int, dur_phase: float = 0.02):\n    \"\"\"\n    time_s [s] → days; look up P,t0 for this planet → phase in (–0.5…+0.5)\n    returns (phase, in_mask, out_mask)\n    \"\"\"\n    time_d = time_s / 86400.0\n    row    = star_info.query(\"planet_id==@pid\").iloc[0]\n    P      = float(row[\"P\"])\n    t0     = float(row.get(\"t0\", 0.5*P))\n    phase  = ((time_d - t0 + 0.5*P) % P)/P - 0.5\n\n    in_m = np.abs(phase) < dur_phase\n    if not in_m.any():  # widen if no points\n        in_m = np.abs(phase) < (dur_phase*2)\n    return phase, in_m, ~in_m\n\n# ─── 8) Depth & σ measurement ─────────────────────────────────────────────────\ndef measure_depths(spec_cube, flux_wl, in_mask, out_mask, eps=1e-8):\n    \"\"\"\n    Compute transit depth channel-by-channel as 1 - (mean in / mean out).\n    Adds a tiny eps when dividing to avoid NaNs.\n    \"\"\"\n    # make sure flux_wl never zero\n    flux_safe = flux_wl.copy()\n    flux_safe[np.abs(flux_safe) < eps] = np.nanmedian(flux_safe)\n\n    # relative spectrum\n    rel_cube = spec_cube / flux_safe[:, None]\n\n    # compute in/out means\n    mu_in  = np.nanmean(rel_cube[in_mask],  axis=0)\n    mu_out = np.nanmean(rel_cube[out_mask], axis=0)\n\n    # avoid zero out:\n    mu_out_safe = mu_out.copy()\n    mu_out_safe[np.abs(mu_out_safe) < eps] = eps\n\n    depth = 1.0 - mu_in / mu_out_safe\n    depth = np.clip(depth, 0.0, 1.0)\n\n    sigma = np.nanstd(rel_cube[in_mask], axis=0) / np.sqrt(np.sum(in_mask))\n\n    return depth, sigma\n\ndef build_submission(default_sigma=2e-3):\n    sample = pd.read_csv(DATA_DIR/\"sample_submission.csv\", nrows=0).columns\n    rows   = []\n    star_info = pd.read_csv(DATA_DIR/\"test_star_info.csv\")\n\n    for d in tqdm(sorted(TEST_ROOT.iterdir()), desc=\"planets\"):\n        pid = int(d.name)\n\n        # 1) CALIBRATE & EXTRACT WHITE LIGHT\n        cube, var        = calibrate_cube(d)\n        t_s, flux_wl, _  = extract_white_light(cube, var)\n\n        # 2) SPECTRAL BOX & NORMALIZE\n        spec_cube = box_extract(cube)\n        spec_rel  = spec_cube / flux_wl[:, None]\n\n        # 3) COMMON‐MODE CORRECTION\n        corr_cube = remove_common_mode(spec_rel)\n\n        # 4) ORBITAL PARAMETERS\n        rowp = star_info[star_info[\"planet_id\"] == pid].iloc[0]\n        P   = float(rowp[\"P\"])\n        dur = float(rowp.get(\"dur\", 0.05 * P))\n        t0  = float(rowp.get(\"t0\", 0.5 * P))\n\n        # 5) PHASE\n        time_d = t_s / 86400.0\n        phase  = ((time_d - t0 + 0.5*P) % P) / P - 0.5\n\n        # 6) IN/OUT MASKS WITH FALLBACK\n        half_w = 0.5 * dur / P\n        in_m   = np.abs(phase) < half_w\n        if in_m.sum() == 0:\n            n = len(phase)\n            lo, hi = int(0.35*n), int(0.65*n)\n            in_m = np.zeros(n, bool)\n            in_m[lo:hi] = True\n        out_m = ~in_m\n\n        # 7) MEASURE DEPTHS & UNCERTAINTIES\n        depths, sigmas = measure_depths(corr_cube, phase, in_m, out_m)\n\n        # 8) ASSEMBLE ROW WITH SANITY‐CHECK\n        row = {\"planet_id\": pid, \"wl_1\": 0.0, \"sigma_1\": default_sigma}\n        for idx, (δ, σ) in enumerate(zip(depths, sigmas), start=2):\n            # if δ is NaN or negative, set to zero\n            if not np.isfinite(δ) or δ < 0:\n                δ = 0.0\n            # if σ is NaN or non-positive, set to default_sigma\n            if not np.isfinite(σ) or σ <= 0:\n                σ = default_sigma\n\n            row[f\"wl_{idx}\"]    = float(δ)\n            row[f\"sigma_{idx}\"] = float(σ)\n\n        rows.append(row)\n\n    submission = pd.DataFrame(rows)[sample]\n    submission.to_csv(WORK_DIR/\"submission.csv\", index=False)\n    print(\"✅ Written submission.csv, shape:\", submission.shape)\n\n# ─── run it! ────────────────────────────────────────────────────────────────\nbuild_submission()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-22T12:42:03.473577Z","iopub.execute_input":"2025-07-22T12:42:03.473879Z","iopub.status.idle":"2025-07-22T12:42:13.200285Z","shell.execute_reply.started":"2025-07-22T12:42:03.473859Z","shell.execute_reply":"2025-07-22T12:42:13.19924Z"}},"outputs":[],"execution_count":null}]}