{"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":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"#### Everything about the data can be found here: https://www.kaggle.com/code/kuntalpal/ariel-data-neurips-2025","metadata":{}},{"cell_type":"markdown","source":"\n# Time for some ML\n\n## 11. What features to use?? \n\n- ### Calibrated pixel cube \n $$\n I_{\\mathrm{cal}}(x,y,t) =\n \\frac{\\bigl(I_{\\mathrm{raw}}(x,y,t)\\times \\mathrm{gain} + \\mathrm{offset}\\bigr)\\;-\\;D(x,y)}\n {F(x,y)}\n \\;\\times\\;\\mathrm{LinCorr}\\bigl(x,y,I_{\\mathrm{raw}}\\bigr)\n $$\n\n where\n• $D$ = dark map (thermal + bias), \\\n• $F$ = flat‑field response, \\\n• $\\mathrm{LinCorr}$ = polynomial linearisation term.\n\n- ### Transit‑phase mask \n\n Compute a boolean mask using the Mandel–Agol analytic model with orbital elements supplied in `*_star_info.csv`.\n\n\n- ### Frame-level (pixel-domain) statistics\n\n| Feature | Formula | Physics intuition |\n|---------------------------|-------------------------------------------------------------------------------------------------------------------|----------------------------------------------------|\n| Total flux $F_{\\rm tot}(t)$ | $\\displaystyle \\sum_{x,y} I_{\\rm cal}(x,y,t)$ | Broadband stellar photon count |\n| Centroid $(x_c, y_c)$ | $\\displaystyle x_c(t)=\\frac{\\sum_{x,y}x\\,I_{\\rm tot}(t)}{\\!F_{\\rm tot}(t)},\\quad y_c(t)=\\frac{\\sum_{x,y}y\\,I_{\\rm tot}(t)}{\\!F_{\\rm tot}(t)}$ | Pointing jitter / PSF wander |\n| Centroid scatter | $\\sigma_{x_c},\\;\\sigma_{y_c}$ over the whole visit | Jitter couples to wavelength-dependent throughput |\n| Background mean $B(t)$ | mean of outer two detector rows/cols | Thermal & zodiacal background |\n| Spatial RMS $\\mathrm{RMS}_{xy}(t)$ | $\\sqrt{\\langle\\,I-\\widetilde{I}\\rangle^2_{x,y}}$ | Hot pixels, Poisson noise proxy |\n| Flux derivative $\\dot F(t)$ | $\\displaystyle \\frac{F_{\\rm tot}(t+1)-F_{\\rm tot}(t)}{\\Delta t}$ | Rapid thermal drifts / reaction wheel hits |\n\nFor each scalar time-series above, store **mean**, **standard deviation**, **skewness** computed *separately* for the out-of-transit (OOT) and in-transit (IN) windows.\n\n\n- ### Spectrophotometric light-curve features\n\n\n - **AIRS‑CH0** (infra‑red, 282 bins) \n $$L_{\\lambda}(t) \\;=\\; \\sum_{y=0}^{31} I_{\\rm cal}\\bigl(x=\\lambda,\\;y,\\;t\\bigr)$$\n\n - **FGS1** (broadband) \n $$L_{\\rm FGS}(t) \\;=\\; \\sum_{x,y} I_{\\rm cal}(x,y,t)$$\n\n\n - Transit‑depth measures\n\n| Quantity | Definition |\n|--------------------|----------------------------------------------------------------------------------------------------------------------------------------------------------|\n| Analytic depth | $$d_{\\rm phys}(\\lambda) = \\frac{R_p^2(\\lambda)}{R_*^2}$$ |\n| Empirical depth | $$d_{\\rm emp}(\\lambda) = \\dfrac{\\mathrm{median}_{\\rm OOT}L(\\lambda) - \\mathrm{median}_{\\rm IN}L(\\lambda)}{\\mathrm{median}_{\\rm OOT}L(\\lambda)}$$ |\n| Residual target | $$r(\\lambda) = d_{\\rm emp}(\\lambda) - d_{\\rm phys}(\\lambda)$$ |\n\n - Correlated‑noise & trend metrics\n\n| Metric | Formula |\n|-------------------------------------|----------------------------------------------------------------------------------------------------------------------------------|\n| White‑trend slope $a_{\\lambda}$ | Fit $L_{\\lambda}(t) = a_{\\lambda}\\,t + b_{\\lambda}$ on OOT data |\n| Correlated‑noise RMS $\\sigma_{\\rm corr,\\lambda}$ | RMS of $L_{\\lambda}(t)$ minus fitted trend (OOT) |\n| Pont $\\beta$‑factor $\\beta_{\\lambda}$ | $$\\beta_{\\lambda} = \\sqrt{\\frac{\\sigma_N^2}{\\sigma_1^2}}$$ comparing std‑dev before/after $N$‑point binning |\n\nStore **$a_{\\lambda}$**, **$\\sigma_{\\rm corr,\\lambda}$**, **$\\beta_{\\lambda}$** for each $\\lambda$ plus their visit‑level **mean** & **std**.\n\n - Broadband aggregates\n\nSystematics often act coherently across wavelength. For any per‑channel statistic $X_{\\lambda}$, compute:\n\n$$\n\\mu_{X} = \\frac{1}{N_{\\lambda}} \\sum_{\\lambda} X_{\\lambda}, \n\\qquad\n\\sigma_{X} = \\sqrt{\\frac{1}{N_{\\lambda} - 1} \\sum_{\\lambda} \\bigl(X_{\\lambda} - \\mu_{X}\\bigr)^{2}}.\n$$\n\nApply to $X \\in \\{d_{\\rm emp},\\,a,\\,\\sigma_{\\rm corr},\\,\\beta\\}$.\n\n- ### Astrophysical & Geometric Metadata \n\nThese **global parameters** describe the star–planet system. We transform or scale each so that the ML model can learn more effectively.\n\n| CSV Field | Stored Value | Interpretation & Physics Rationale |\n|-----------|---------------------------------------|-------------------------------------------------------------------------------|\n| **`Rs`** | $\\ln\\bigl(R_s/R_\\odot\\bigr)$ | Stellar radius in Solar radii (log‐scaled). Sets the transit **depth** scale: $d_{\\rm phys}\\propto\\frac{R_p^2}{R_s^2}\\,. $ |\n| **`Ms`** | $\\ln\\bigl(M_s/M_\\odot\\bigr)$ | Stellar mass in Solar masses (log‐scaled). Determines surface gravity $g_s$ and thus the planet’s **scale height** $H = \\frac{k_B T_{\\rm eq}}{\\mu m_p g_s}\\,.$ |\n| **`Ts`** | $T_s/5000\\,$K | Stellar effective temperature (scaled to $\\sim5000\\,$K). Influences equilibrium temperature of the planet and limb‐darkening. |\n| **`Mp`** | $\\ln\\bigl(M_p/M_\\oplus\\bigr)$ | Planet mass in Earth masses (log‐scaled). Enters the scale height via $g_p\\propto M_p/R_p^2$. |\n| **`P`** | $\\log_{10}\\bigl(P/\\mathrm{days}\\bigr)$ | Orbital period in days (log‐scaled). Sets the **transit duration** $T_{14}\\approx\\frac{P}{\\pi}\\arcsin\\!\\Bigl(\\frac{R_s}{a}\\Bigr)\\,.$ |\n| **`sma`** | Raw | Semi-major axis in units of $R_s$. Controls ingress/egress timing and equilibrium temperature $T_{\\rm eq}\\propto\\Bigl(\\frac{R_s}{2\\,a}\\Bigr)^{1/2}T_s\\,.$ |\n| **`i`** | $\\cos(i\\,\\pi/180)$ | Orbital inclination in degrees, converted via cosine so $i\\approx90^\\circ$ (central transit) maps to 0 smoothly. |\n| **`e`** | Raw | Orbital eccentricity. Affects transit asymmetry and time‐varying star–planet separation: $r(\\theta)=\\frac{a(1-e^2)}{1+e\\cos\\theta}\\,.$ |\n\n---\n\n### Why these matter\n\n- **Transit geometry** \n - Impact parameter: \n $$b = \\frac{a\\cos i}{R_s}\\,,\\quad a = \\texttt{sma}\\,\\times R_s.$$ \n - Duration & shape set by $P$, $sma$, $R_s$, and $i$.\n\n- **Atmospheric scale height** $H$ \n $$H = \\frac{k_B\\,T_{\\rm eq}}{\\mu\\,m_p\\,g_p},\\quad g_p = \\frac{G\\,M_p}{R_p^2},$$ \n so $M_s$, $T_s$, and $M_p$ control the amplitude of spectral features.\n\n- **Stellar limb-darkening & temperature** \n - $T_s$ sets the baseline spectral energy distribution; helps the model separate instrumental trends from real color changes.\n\nBy **log‐scaling** radii/masses and **cosine‐transforming** inclination, we avoid large numerical ranges and boundary effects, speeding up learning. \n","metadata":{}},{"cell_type":"markdown","source":"### Target & loss for ML training\n\n$$\n\\mathcal{L}\n= \\sum_{\\lambda}\n\\frac{\\bigl(r_{\\text{pred}}(\\lambda) - r_{\\text{true}}(\\lambda)\\bigr)^{2}}\n {2\\,\\sigma_{\\text{ML}}^{2}(\\lambda)}\n\\;+\\;\n\\frac{1}{2}\\,\\ln\\sigma_{\\text{ML}}^{2}(\\lambda)\n$$\n\n$$\n\\mu_{\\text{user}}(\\lambda)\n= d_{\\text{phys}}(\\lambda) + r_{\\text{pred}}(\\lambda),\n\\qquad\n\\sigma_{\\text{user}}(\\lambda)\n= \\sigma_{\\text{ML}}(\\lambda).\n$$\n\n---\n\n### Summary checklist\n\n- **Pixel-domain stats:** $F_{\\text{tot}},\\, x_c,\\, y_c,\\, B,\\, \\mathrm{RMS}_{xy},\\, \\dot{F}$ (mean/std **IN** & **OOT**)\n- **Per-channel metrics:** $d_{\\text{phys}},\\, d_{\\text{emp}},\\, a_\\lambda,\\, \\sigma_{\\text{corr},\\lambda},\\, \\beta_\\lambda$\n- **Broadband aggregates:** $\\mu_X,\\, \\sigma_X$\n- **Astro metadata:** $(R_\\ast,\\, M_\\ast,\\, T_\\ast,\\, M_p,\\, P,\\, a/R_\\ast,\\, \\cos i,\\, e)$\n\n*Total: ≈ 500 scalar features — ready for LightGBM or extension to deep models.*\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nfrom sklearn.linear_model import LinearRegression\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Calibrated Pixel Cube\ndef calibrate_cube(raw_cube, gain, offset, dark_map, flat_map, lin_corr_map):\n \"\"\"\n Restore calibrated pixel values.\n \n Parameters\n ----------\n raw_cube : np.ndarray\n Raw pixel cube, shape (T, H, W), dtype uint16.\n gain : float\n Detector gain.\n offset : float\n Detector offset.\n dark_map : np.ndarray\n Dark current + bias map, shape (H, W).\n flat_map : np.ndarray\n Flat-field response map, shape (H, W).\n lin_corr_map : np.ndarray\n Linearity correction factors, shape (H, W).\n \n Returns\n -------\n I_cal : np.ndarray\n Calibrated pixel cube, shape (T, H, W), dtype float.\n \"\"\"\n # Recover dynamic range\n I = raw_cube.astype(float) * gain + offset\n # Subtract dark current, divide flat field\n I = (I - dark_map[None, :, :]) / flat_map[None, :, :]\n # Apply linearity correction\n I *= lin_corr_map[None, :, :]\n return I\n\n# Transit-Phase Mask (approximate)\ndef transit_phase_mask(times, P, sma, inc_deg):\n \"\"\"\n Compute in-transit and out-of-transit masks without external packages.\n \n Parameters\n ----------\n times : np.ndarray\n Time stamps of observations, shape (T,), same units as P.\n P : float\n Orbital period (days).\n sma : float\n Semi-major axis in units of stellar radii (R*).\n inc_deg : float\n Orbital inclination in degrees.\n \n Returns\n -------\n in_transit : np.ndarray\n Boolean mask for in-transit frames, shape (T,).\n oot : np.ndarray\n Boolean mask for out-of-transit frames, shape (T,).\n \"\"\"\n # Compute orbital phase in [0,1)\n phase = ((times - times[0]) / P) % 1.0\n \n # Approximate full transit duration (T14) using small-planet formula\n i = np.deg2rad(inc_deg)\n b = sma * np.cos(i) # impact parameter\n arg = np.clip(np.sqrt(max(0, (1.0) - b**2)) / (sma * np.sin(i)), -1, 1)\n T14 = (P / np.pi) * np.arcsin(arg)\n \n # Fractional duration of transit\n frac = T14 / P\n \n # Center transit at phase = 0.5\n in_transit = (phase >= 0.5 - frac/2) & (phase <= 0.5 + frac/2)\n oot = ~in_transit\n return in_transit, oot\n\n# Example usage (replace raw_cube, times, star parameters as needed):\n# I_cal = calibrate_cube(raw_cube, gain, offset, dark_map, flat_map, lin_corr_map)\n# in_tr, oot = transit_phase_mask(times, P, sma, inclination)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T08:36:49.67692Z","iopub.execute_input":"2025-08-03T08:36:49.677628Z","iopub.status.idle":"2025-08-03T08:36:49.69068Z","shell.execute_reply.started":"2025-08-03T08:36:49.677586Z","shell.execute_reply":"2025-08-03T08:36:49.689608Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Frame-level (pixel-domain) statistics\ndef frame_stats(I_cal, in_transit, oot):\n \"\"\"\n Compute per-frame and aggregated statistics:\n - Total flux\n - Centroids x_c, y_c\n - Centroid scatter\n - Background mean (border pixels)\n - Spatial RMS\n - Flux derivative\n \n Aggregates mean, std, skew for IN and OOT segments.\n \n Parameters\n ----------\n I_cal : np.ndarray, shape (T, H, W)\n Calibrated pixel cube.\n in_transit : np.ndarray, bool, shape (T,)\n In-transit mask.\n oot : np.ndarray, bool, shape (T,)\n Out-of-transit mask.\n \n Returns\n -------\n stats : dict\n Feature dictionary with keys like 'F_IN_mean', 'F_OOT_std', etc.\n \"\"\"\n T, H, W = I_cal.shape\n flat = I_cal.reshape(T, -1)\n \n # 1. Total flux\n F = flat.sum(axis=1)\n \n # 2. Centroids\n ys, xs = np.indices((H, W))\n x_c = (I_cal * xs[None]).sum(axis=(1,2)) / F\n y_c = (I_cal * ys[None]).sum(axis=(1,2)) / F\n \n # 3. Background mean (outer border)\n border = np.concatenate([\n I_cal[:, 0, :], I_cal[:, -1, :], \n I_cal[:, :, 0], I_cal[:, :, -1]\n ], axis=0).reshape(4, T, -1)\n B = border.mean(axis=(0,2))\n \n # 4. Spatial RMS\n med = np.median(flat, axis=1)\n RMS = np.sqrt(((flat - med[:, None])**2).mean(axis=1))\n \n # 5. Flux derivative\n dotF = np.concatenate(([0], np.diff(F))) # assume uniform Δt\n \n # Helper for aggregation\n def agg(arr):\n return arr.mean(), arr.std(), skew(arr)\n \n stats = {}\n for name, arr in [('F', F),\n ('x_c', x_c),\n ('y_c', y_c),\n ('B', B),\n ('RMS', RMS),\n ('dotF', dotF)]:\n m_in, s_in, k_in = agg(arr[in_transit])\n m_oot, s_oot, k_oot = agg(arr[oot])\n stats[f'{name}_IN_mean'] = m_in\n stats[f'{name}_IN_std'] = s_in\n stats[f'{name}_IN_skew'] = k_in\n stats[f'{name}_OOT_mean'] = m_oot\n stats[f'{name}_OOT_std'] = s_oot\n stats[f'{name}_OOT_skew'] = k_oot\n \n return stats\n\n\n# Example usage:\n# pixel_feats = frame_stats(I_cal, in_tr, oot)\n","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\ndef spectro_lightcurve_features(I_cal_airs, I_cal_fgs, in_transit, oot, d_phys_airs=None, N_beta=10):\n \"\"\"\n Extract spectrophotometric light-curve features for AIRS-CH0 and FGS1.\n \n Parameters\n ----------\n I_cal_airs : np.ndarray, shape (T, H, W)\n Calibrated cube for AIRS-CH0.\n I_cal_fgs : np.ndarray, shape (T, H, W')\n Calibrated cube for FGS1 (H'=W'=32).\n in_transit : np.ndarray, bool, shape (T,)\n In-transit mask.\n oot : np.ndarray, bool, shape (T,)\n Out-of-transit mask.\n d_phys_airs : np.ndarray or None, shape (W,)\n Analytic physics depth per AIRS channel. If None, zeros used.\n N_beta : int\n Bin size for beta factor calculation.\n \n Returns\n -------\n feats : dict\n Dictionary containing per-channel arrays and aggregated stats:\n - 'd_emp': empirical depth array for AIRS shape (W,)\n - 'r': residual array for AIRS shape (W,)\n - 'a': white-trend slopes array for AIRS shape (W,)\n - 'sigma_corr': correlated-noise RMS array shape (W,)\n - 'beta': Pont beta-factor array shape (W,)\n - '_mean', '_std' for each feat in [d_emp, r, a, sigma_corr, beta]\n - 'd_emp_fgs', 'r_fgs' for FGS1 empirical depth & residual\n \"\"\"\n T = I_cal_airs.shape[0]\n # 1) Light curves\n L_airs = I_cal_airs.sum(axis=1) # shape (T, W)\n L_fgs = I_cal_fgs.sum(axis=(1,2)) # shape (T,)\n \n # 2) Analytic physics depth\n W = L_airs.shape[1]\n if d_phys_airs is None:\n d_phys_airs = np.zeros(W)\n \n # 3) Empirical depth & residual for AIRS\n med_oot = np.median(L_airs[oot], axis=0)\n med_in = np.median(L_airs[in_transit], axis=0)\n d_emp = (med_oot - med_in) / med_oot\n r = d_emp - d_phys_airs\n \n # 4) Trend & correlated noise per channel\n t = np.arange(T)[:, None]\n a = np.zeros(W)\n sigma_corr = np.zeros(W)\n beta = np.zeros(W)\n for i in range(W):\n # slope fit on OOT\n lr = LinearRegression().fit(t[oot], L_airs[oot, i])\n a[i] = lr.coef_[0]\n resid = L_airs[:, i] - lr.predict(t)\n sigma_corr[i] = np.sqrt(np.mean(resid[oot]**2))\n # beta factor\n M = (T // N_beta) * N_beta\n binned = resid[:M].reshape(-1, N_beta).mean(axis=1)\n beta[i] = np.std(binned) / np.std(resid[:M])\n \n # 5) Empirical depth & residual for FGS1\n med_oot_fgs = np.median(L_fgs[oot])\n med_in_fgs = np.median(L_fgs[in_transit])\n d_emp_fgs = (med_oot_fgs - med_in_fgs) / med_oot_fgs\n r_fgs = d_emp_fgs # if no physics model for FGS, residual = empirical\n \n # 6) Aggregates\n def agg(arr):\n return arr.mean(), arr.std()\n \n feats = {\n 'd_emp': d_emp,\n 'r': r,\n 'a': a,\n 'sigma_corr': sigma_corr,\n 'beta': beta,\n 'd_emp_fgs': d_emp_fgs,\n 'r_fgs': r_fgs,\n }\n \n for name in ['d_emp', 'r', 'a', 'sigma_corr', 'beta']:\n m, s = agg(feats[name])\n feats[f'{name}_mean'] = m\n feats[f'{name}_std'] = s\n \n return feats\n\n# Example usage:\n# feats_spec = spectro_lightcurve_features(I_cal_airs, I_cal_fgs, in_tr, oot, d_phys_airs=None)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-08-03T09:16:11.493297Z","iopub.execute_input":"2025-08-03T09:16:11.498512Z","iopub.status.idle":"2025-08-03T09:16:13.725405Z","shell.execute_reply.started":"2025-08-03T09:16:11.498381Z","shell.execute_reply":"2025-08-03T09:16:13.72442Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}