{"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,"isSourceIdPinned":false,"sourceType":"competition"},{"sourceId":9629432,"sourceType":"datasetVersion","datasetId":5846888}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# はじめに（ノートブック冒頭の説明文）\n\n## 🌌 NeurIPS Ariel Data Challenge 2025 への挑戦：物理に素直な前処理 × シンプルなトランジット推定\n\nこのノートブックでは、**Ariel 衛星の模擬観測データ**（FGS1 + AIRS-CH0）から、\n**惑星トランジットによる極小の減光（トランジット深さ）を安定に推定**するための、\n\\*\\*“物理に基づく前処理 → 低自由度フィッティング”\\*\\*というミニマルな方針をまとめます。\n\n> 目標：生データ（ADU）から**科学的に妥当なキャリブレーション**を行い、\n> 時系列の滑らかさを指標に**トランジット深さ $s$** を1スカラーで推定するベースラインを作ること。\n\n---\n\n## 🎯 問題の要点（何を解くか）\n\n* 入力：FGS1（ガイドカメラ）と AIRS-CH0（分光器）の**生信号**（ADU, 未フラックス較正）。\n* キャリブレーション：ダーク電流・フラット・非線形性・死/ホットピクセルなど**検出器依存の系統**を補正。\n* 読み出し：\\*\\*Correlated Double Sampling（CDS）\\*\\*で各露光の開始/終了の差を取り、**オフセットを打ち消す**。\n* 目的：**トランジット中の微小な減光**を時系列から取り出す（オーダーは \\~100 ppm レベル）。\n* 出力：各惑星について**1つのトランジット深さ $s$** と、その共通シグマ（提出形式に合わせて複製）。\n\n---\n\n## 🧭 このノートのアプローチ（全体像）\n\n本ワークは、2024 年版の校正レシピを踏襲しつつ、2025 データ仕様に合わせて**一体化した前処理クラス**を設計。\nデータ縮約を徹底し、**軽量な最適化**でスカラー $s$ を推定します。\n\n1. **ADC 逆変換**：ADU → アナログ相当へ（gain/offset で戻す）。\n2. **死/ホットピクセルのマスク**：ダークからホット推定（sigma-clip）、dead/ hot を mask。\n3. **（任意）非線形補正**：各画素ごとの逆多項式で応答の非線形を補正。\n4. **ダーク減算**：**積分時間に比例**させてスケールし、フレーム毎に減算。\n5. **CDS**：各露光の（終−初）で**加法的バイアスを打ち消す**。\n6. **時間ビニング**：AIRS=30枚、FGS1=360枚（30×12）で平均し、**時系列長を \\~187** に圧縮。\n7. **フラット補正**：**ビニング後**にフラット/デッドを適用し、相対感度のむらを補正。\n8. **空間縮約**：\n\n   * AIRS-CH0：**中心行（12行）の平均**で空間ノイズを下げ、**波長次元（282）を保持**。\n   * FGS1：**中心 12×12** を平均して **1 チャンネル**化。\n9. **結合**：最終形状を **(N, 187, 283)**（FGS1=1 + AIRS=282）に整形。\n\nその上で、**非 ML** のミニマルモデルで $s$ を推定：\n\n* **トランジット位相検出**：平滑化した光度曲線からディップの開始/終了（phase1/phase2）を自動検出。\n* **目的関数**：トランジット区間の信号を $(1+s)$ で持ち上げたとき、**全体に低次数多項式が最も滑らかに当たる $s$** を選ぶ。\n  直感的には、**“ディップを埋めて滑らかさ（非連続の小ささ）を最大化”** する $s$ が真の深さに対応。\n\n---\n\n## 🔬 なぜこの順番・設計か（意図）\n\n* **CDS前にダークを適切にスケールして引く**：加法バイアスを残すと**相対比**の分母が歪み、深さがバイアスします。\n* **フラットはビニング後**：S/N が上がった段階で割ると**数値安定**かつ**欠損伝播が穏やか**。\n* **AIRS は波長軸（282）を保持**し、空間だけ平均：分光情報を落とさず、**位相検出のロバスト性**を確保。\n* **FGS1 はガイド的に 1ch 化**：不要な空間ゆらぎを抑えて**参照光度**に。\n* **スカラ $s$**：本コンペの提出形式に沿い、**軽量・頑健**・**解釈容易**。後で **$s(\\lambda)$** へ拡張も可能。\n\n---\n\n## 🔗 参考にした議論・ノートブック（謝辞）\n\n* **Calibrating and Binning Ariel Data**（2024→2025 反映版）\n  [https://www.kaggle.com/code/gordonyip/calibrating-and-binning-ariel-data](https://www.kaggle.com/code/gordonyip/calibrating-and-binning-ariel-data)\n  → 校正の順序（ADC/ダーク/フラット/非線形/マスク）と CDS・ビニングの実務的レシピ。\n\n* **NeurIPS Non-ML Transit Curve Fitting**\n  [https://www.kaggle.com/code/vitalykudelya/neurips-non-ml-transit-curve-fitting](https://www.kaggle.com/code/vitalykudelya/neurips-non-ml-transit-curve-fitting)\n  → **低自由度フィット**でのトランジット深度推定という思想、**位相検出 + 最小化**の枠組み。\n\nこれらを踏まえ、本ノートでは**クラス化して一気通貫の再現性**を重視し、\n**I/O と形状（(N,187,283)）を厳密に保証**する実装へとまとめています。\n\n---\n\n## 🧪 このノートの流れ（各セルの役割）\n\n* **Cell 1**：`pqdm` をオフラインインストール（並列化オプション）。\n* **Cell 2**：インポート・`Config` 定義（パス、切り出し範囲、ビニング、トグル）。\n* **Cell 3**：`UnifiedArielPreprocessor`（前処理の**全工程を一括**で実施）。\n* **Cell 4**：`TransitModel`（位相検出 → $s$ 最適化）。\n* **Cell 5**：`SubmissionGenerator`（提出ファイル `submission.csv` の生成）。\n* **Cell 6**：実行ブロック（前処理→推定→提出）。\n\n---\n\n## 🚀 今後の発展（やってみる価値が高い拡張）\n\n* $s$ を **波長依存 $s(\\lambda)$** に拡張（物理的拘束＋滑らかさ正則化）。\n* **システマティックモデル**（温度ドリフト・姿勢等）を外生変数で同時回帰。\n* **ロバスト回帰**／**分位点損失**で外れ値耐性を強化。\n* **FGS1 の使い方**：参照光のトレンド除去により AIRS の安定化。\n* **ベイズ推定**（位相と深さの同時事後）で不確かさ伝播を厳密化。\n\n---\n","metadata":{}},{"cell_type":"code","source":"!pip install --no-index --find-links=/kaggle/input/ariel-2024-pqdm pqdm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-07T09:26:34.442841Z","iopub.execute_input":"2025-09-07T09:26:34.443499Z","iopub.status.idle":"2025-09-07T09:26:38.357062Z","shell.execute_reply.started":"2025-09-07T09:26:34.443472Z","shell.execute_reply":"2025-09-07T09:26:38.355832Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# セットアップ解説\n\n## 🔧 Cell 1 で何をしているか\n\n* `pqdm` を**オフライン**でインストールしています。\n  Kaggle 環境のネットワーク制限でも、提供ホイールから**並列実行**を有効化できるようにするためです。\n  並列化は **前処理（惑星ごと）** の速度を大きく改善します（CPU ワーカー × I/O 並列）。\n\n> もし `pqdm` が使えない（例：環境差異）場合は、コード側で自動的に**逐次実行**へフォールバックします。\n\n---\n\n## 📦 Cell 2 の役割：ライブラリ読込と Config 定義\n\nここでは**必要ライブラリの読み込み**と、前処理・モデリングの振る舞いを決める**設定クラス `Config`** を定義します。\nこの `Config` を変えるだけで、**パス・切り出し・校正 ON/OFF・ビニング幅・並列度**などを一括調整できます。\n\n### `Config` の主な項目（ざっくり）\n\n* **DATA\\_PATH / DATASET / N\\_JOBS**\n  データ格納場所・対象（`train`/`test`）・並列ジョブ数。\n* **ADC パラメータ（フォールバック）**\n  `adc_info.csv` が無い場合に使うゲイン/オフセット（ADU→アナログ相当の逆変換）。\n* **スペクトル切り出し**\n  `CUT_INF=39, CUT_SUP=321` → AIRS-CH0 の**波長側 282 チャンネル**を採用（321 は**排他的**上限）。\n* **パイプライン・トグル**\n  `DO_MASK/DO_NL_CORR/DO_DARK/DO_FLAT` で**個別の校正工程**を ON/OFF。\n  まずは `DO_MASK=True, DO_DARK=True, DO_FLAT=True, DO_NL_CORR=False` が堅い出発点。\n* **露光パターン（CDS 用）**\n  `AIRS_DT_BASE=0.1, AIRS_DT_INC=4.5` / `FGS_DT_BASE=0.1, FGS_DT_INC=0.1`\n  → 奇数フレームで**追加積分時間**を持つ前提（終−初の CDS に合わせたスケーリング）。\n* **時間ビニング幅**\n  `AIRS_BIN=30`, `FGS_BIN=30*12`（=360）。CDS 後の**平均化**で時系列長を \\~187 に圧縮。\n* **形状と ROI**\n  センサーごとの生形状、非線形係数形状、**中心 12 行/12×12** の空間 ROI。\n* **モデル関連**\n  位相検出用スライス、ギャップ `DELTA`、多項式次数、スケール・既定シグマ。\n\n---\n\n## 🧭 形状の約束（とくに重要）\n\n最終的にモデルへ渡すテンソルは **(N, 187, 283)**。\n\n* `187`：CDS + ビニング後の**時間**長\n* `283`：**FGS1=1**（中心 12×12 平均）＋ **AIRS-CH0=282**（波長 282 を保持）\n\n中間形状は下記の流れで変わります：\n\n| 段階               | FGS1                 | AIRS-CH0                  |\n| ---------------- | -------------------- | ------------------------- |\n| **生信号**          | (135000, 32, 32)     | (11250, 32, 356)          |\n| **スペクトル切り出し**    | —                    | (11250, 32, **282**)      |\n| **CDS（終−初）**     | (\\~67500, 32, 32)    | (\\~5625, 32, 282)         |\n| **ビニング**         | (1, **187**, 32, 32) | (1, **187**, **282**, 32) |\n| **フラット補正**       | (同上)                 | (同上)                      |\n| **空間縮約（ROI 平均）** | (1, 187, **1**)      | (1, 187, **282**)         |\n| **結合**           | -                    | **(N, 187, 283)**         |\n\n> **注意**：AIRS は**波長（W=282）を温存**し、**空間 X=32 行**側を中心 12 行で平均します。\n> FGS1 は中心 12×12 を平均して **1 チャンネル化**します。\n\n---\n\n## 🎚️ 前処理トグルの使い方（おすすめの流儀）\n\n* **DO\\_MASK**：`True` 推奨。dead/hot を**マスク**するだけでもノイズが一段下がる。\n* **DO\\_NL\\_CORR**：まずは `False` で開始。学習や最適化が安定したら `True` で再検証。\n* **DO\\_DARK**：`True` 推奨。**積分時間でスケール**してから減算するのがコツ。\n* **DO\\_FLAT**：`True` 推奨。**ビニング後**に適用することで数値安定・欠損伝播の抑制。\n\n---\n\n## 🧱 よくある落とし穴と対策\n\n* **軸の取り違え**：AIRS のレイアウトは「(T, 32, 282) → (1, T, **282**, **32**)」に**transpose**。\n* **CDS 入力次元**：`_cds` は **3D/4D 両対応**（(T,W,X) or (B,T,W,X)）。\n* **masked array の挙動**：CDS 前に `np.asarray` で**ndarray に戻す**と安全。\n* **並列実行の例外**：`pqdm` は**例外オブジェクト**を返すことあり → `process_all()` 内で検出して即時に再スロー。\n* **時系列長の端数**：ビニングは**切り捨て**。187 に満たないケースは設定見直しを。\n\n---\n\n## ⚙️ 実行のヒント\n\n* **まずは逐次実行**で形状と動作を確認 → 問題なければ `N_JOBS` を上げて並列化。\n* **ログ**：惑星 ID ごとに失敗箇所が分かるよう `RuntimeError` に ID を添付済み。\n* **再現性**：`Config` だけで全体の再実行条件を固定できる構成にしています。\n\n---\n\n","metadata":{}},{"cell_type":"code","source":"import os, glob, itertools\nimport numpy as np\nimport pandas as pd\n\nfrom tqdm import tqdm\nfrom astropy.stats import sigma_clip\nfrom scipy.signal import savgol_filter\nfrom scipy.optimize import minimize\n\n# ---- Optional parallel\ntry:\n    from pqdm.threads import pqdm\n    PQDM_AVAILABLE = True\nexcept Exception:\n    PQDM_AVAILABLE = False\n\n\nclass Config:\n    # paths\n    DATA_PATH  = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET    = 'test'      # 'train' or 'test'\n    N_JOBS     = 4           # for parallel preprocessing if pqdm is available\n\n    # ADC (fallback defaults if adc_info.csv missing)\n    ADC_GAIN_DEFAULT   = 0.4369\n    ADC_OFFSET_DEFAULT = -1000.0\n\n    # spectral crop (keep 39..320 -> 282 bands)\n    CUT_INF = 39\n    CUT_SUP = 321  # exclusive\n\n    # pipeline toggles\n    DO_MASK     = True\n    DO_NL_CORR  = False\n    DO_DARK     = True\n    DO_FLAT     = True\n\n    # timing pattern for CDS pairs\n    AIRS_DT_BASE, AIRS_DT_INC = 0.1, 4.5\n    FGS_DT_BASE,  FGS_DT_INC  = 0.1, 0.1\n\n    # binning\n    AIRS_BIN = 30\n    FGS_BIN  = 30 * 12\n\n    # shapes\n    AIRS_RAW_SHAPE = (11250, 32, 356)\n    FGS_RAW_SHAPE  = (135000, 32, 32)\n    LIN_AIRS_SHAPE = (6, 32, 356)\n    LIN_FGS_SHAPE  = (6, 32, 32)\n\n    # ROI (center crop used before spatial aggregation)\n    ROI_ROWS = slice(10, 22)     # 12 rows\n    ROI_COLS = slice(10, 22)     # 12 cols for FGS\n\n    # model knobs (kept same for your TransitModel)\n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA    = 7\n    MODEL_POLYNOMIAL_DEGREE     = 3\n    SCALE = 0.95\n    SIGMA = 0.0009\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-07T09:26:38.359047Z","iopub.execute_input":"2025-09-07T09:26:38.359437Z","iopub.status.idle":"2025-09-07T09:26:38.368863Z","shell.execute_reply.started":"2025-09-07T09:26:38.359392Z","shell.execute_reply":"2025-09-07T09:26:38.36776Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n## 📚 このセルでやること\n\n* 解析に必要な **標準ライブラリ / 科学計算ライブラリ** を読み込みます。\n* 前処理とモデルの振る舞いを一括制御する **`Config` クラス** を定義します。\n  以降のセルは **`Config` を差し替えるだけ**で同じコードを繰り返し再現できます。\n\n---\n\n## 🧩 依存ライブラリの役割（ざっくり）\n\n* **NumPy / Pandas**：数値計算・Parquet/CSV 読み込み・配列操作\n* **tqdm**：ループ進捗の可視化\n* **Astropy (`sigma_clip`)**：ダークからのホットピクセル検出（統計的外れ値マスク）\n* **SciPy**：\n\n  * `savgol_filter`：光度曲線の**スムージング**（位相検出を安定化）\n  * `minimize`：トランジット深さ $s$ の**最適化**\n* **pqdm（任意）**：惑星ごとの前処理を**並列**で高速化（使えない場合は自動で逐次にフォールバック）\n\n---\n\n## ⚙️ `Config` の主な項目と意味\n\n* **パス / 実行**\n\n  * `DATA_PATH`：データのルートディレクトリ\n  * `DATASET`：`'train'` or `'test'`\n  * `N_JOBS`：並列ジョブ数（`pqdm` 使用時）\n* **ADC 逆変換（フォールバック）**\n\n  * `ADC_GAIN_DEFAULT`, `ADC_OFFSET_DEFAULT`：`adc_info.csv` が無いときの既定値\n* **スペクトル切り出し**\n\n  * `CUT_INF=39`, `CUT_SUP=321`（**321 は上限排他**）→ AIRS の波長 282 チャンネルを採用\n* **パイプライン・トグル**\n\n  * `DO_MASK` / `DO_NL_CORR` / `DO_DARK` / `DO_FLAT`\n  * 推奨スタート：`True, False, True, True`（非線形は後から検証）\n* **CDS 用の積分時間パターン**\n\n  * `AIRS_DT_BASE, AIRS_DT_INC = 0.1, 4.5`\n  * `FGS_DT_BASE,  FGS_DT_INC  = 0.1, 0.1`\n    → 奇数フレームが「終了リード」で**追加積分**される前提。\n* **時間ビニング**\n\n  * `AIRS_BIN=30`, `FGS_BIN=30*12`（=360）→ どちらも**約 187 ステップ**へ縮約\n* **形状と非線形係数形状**\n\n  * `AIRS_RAW_SHAPE=(11250, 32, 356)` / `FGS_RAW_SHAPE=(135000, 32, 32)`\n  * `LIN_AIRS_SHAPE=(6, 32, 356)` / `LIN_FGS_SHAPE=(6, 32, 32)`\n* **空間 ROI**\n\n  * `ROI_ROWS = slice(10, 22)`（中心 12 行）\n  * `ROI_COLS = slice(10, 22)`（FGS1 の中心 12 列）\n* **モデル（後段の位相検出・最適化）**\n\n  * `MODEL_PHASE_DETECTION_SLICE = slice(30, 140)`\n  * `MODEL_OPTIMIZATION_DELTA = 7`（端の保護）\n  * `MODEL_POLYNOMIAL_DEGREE = 3`（滑らかさの尺度）\n  * `SCALE = 0.95`, `SIGMA = 0.0009`（提出用のスケール/共通シグマ）\n\n---\n\n## 🎛️ チューニングのコツ\n\n* **AIRS の波長窓**（`CUT_*`）は、S/N や系統の見え方に応じて再検討可。\n* **ビニング幅**は S/N と時間分解能のトレードオフ。`AIRS_BIN` と `FGS_BIN` を同時に変えると整合が取りやすい。\n* **DO\\_NL\\_CORR** を有効化する場合は、**線形性補正の係数**が信号域（飽和近傍を避ける等）に適しているかを確認。\n\n---\n\n## 🧪 確認ポイント\n\n* `Config` を `cfg = Config()` で生成できる\n* `DATA_PATH/DATASET` の中に `*_star_info.csv` とセンサー別 Parquet 群が存在\n* 以降、**Cell 3** の前処理クラスが `cfg` を受け取り、その設定に従って動作します。\n","metadata":{}},{"cell_type":"code","source":"class UnifiedArielPreprocessor:\n    \"\"\"\n    One-stop class that:\n      - Reads per-planet AIRS-CH0 & FGS1\n      - ADC invert\n      - mask hot/dead\n      - (optional) non-linearity inverse poly\n      - dark subtraction scaled by integration time\n      - CDS (end - start)\n      - time binning\n      - flat-field division (per-planet, safe masking)\n      - spatial ROI + aggregation\n      - Concatenate into (N, 187, 283): [FGS1(187,1) | AIRS(187,282)]\n    \"\"\"\n\n    def __init__(self, cfg: Config):\n        self.cfg = cfg\n        # optional adc table\n        adc_csv = os.path.join(cfg.DATA_PATH, 'adc_info.csv')\n        self.adc_info = pd.read_csv(adc_csv) if os.path.exists(adc_csv) else None\n\n        # planet IDs from star_info\n        ids_csv = os.path.join(cfg.DATA_PATH, f'{cfg.DATASET}_star_info.csv')\n        self.planet_ids = pd.read_csv(ids_csv, index_col='planet_id').index.astype(int).tolist()\n\n    # ---------- utilities ----------\n    @staticmethod\n    def _adc_convert(signal, gain, offset):\n        signal = np.asarray(signal, dtype=np.float64)\n        signal /= gain\n        signal += offset\n        return signal\n\n    def _adc_for(self, sensor: str):\n        if self.adc_info is not None:\n            g = float(self.adc_info[f\"{sensor}_adc_gain\"].iloc[0])\n            o = float(self.adc_info[f\"{sensor}_adc_offset\"].iloc[0])\n            return g, o\n        return self.cfg.ADC_GAIN_DEFAULT, self.cfg.ADC_OFFSET_DEFAULT\n\n    def _dt_pattern(self, sensor: str, length: int):\n        if sensor == 'AIRS-CH0':\n            base, inc = self.cfg.AIRS_DT_BASE, self.cfg.AIRS_DT_INC\n        else:\n            base, inc = self.cfg.FGS_DT_BASE, self.cfg.FGS_DT_INC\n        dt = np.ones(length, dtype=np.float64) * base\n        dt[1::2] += inc\n        return dt\n\n    @staticmethod\n    def _mask_hot_dead(signal, dead, dark):\n        \"\"\"Dead & hot masking (hot from sigma-clip on dark).\"\"\"\n        signal = np.asarray(signal)  # avoid masked-array surprises later\n        hot_mask = sigma_clip(dark, sigma=5, maxiters=5).mask\n        hot  = np.tile(hot_mask,  (signal.shape[0], 1, 1))\n        dead = np.tile(dead,      (signal.shape[0], 1, 1))\n        sig = np.ma.masked_where(dead, signal)\n        sig = np.ma.masked_where(hot,  sig)\n        return sig\n\n    @staticmethod\n    def _apply_linear_corr(linear_corr, clean_signal):\n        \"\"\"Inverse polynomial per pixel across time.\"\"\"\n        lin = np.flip(linear_corr, axis=0)  # highest degree first for np.poly1d\n        out = np.asarray(clean_signal, dtype=np.float64).copy()\n        T, X, W = out.shape\n        for x in range(X):\n            for w in range(W):\n                poly = np.poly1d(lin[:, x, w])\n                out[:, x, w] = poly(out[:, x, w])\n        return out\n\n    @staticmethod\n    def _clean_dark(signal, dead, dark, dt):\n        signal = np.asarray(signal)\n        dark_m = np.ma.masked_where(dead, dark)\n        dark_m = np.tile(dark_m, (signal.shape[0], 1, 1))\n        return signal - dark_m * dt[:, np.newaxis, np.newaxis]\n\n    @staticmethod\n    def _cds(signal):\n        \"\"\"\n        Correlated Double Sampling:\n        - If signal is 4D: (B, T, W, X) → return (B, T/2, W, X)\n        - If signal is 3D: (T, W, X)     → return (T/2, W, X)\n        \"\"\"\n        if signal.ndim == 4:\n            return signal[:, 1::2, :, :] - signal[:, ::2, :, :]\n        elif signal.ndim == 3:\n            return signal[1::2, :, :] - signal[0::2, :, :]\n        else:\n            raise ValueError(f\"_cds expects 3D or 4D, got shape {signal.shape}\")\n\n    @staticmethod\n    def _time_bin_cds(cds_btxw, binning):\n        \"\"\"Average over time windows; drop remainder. cds_btxw: (B, T, W, X)\"\"\"\n        cds_btxw = np.asarray(cds_btxw)\n        B, T, W, X = cds_btxw.shape\n        n_bins = T // binning\n        out = np.empty((B, n_bins, W, X), dtype=cds_btxw.dtype)\n        for j in range(n_bins):\n            out[:, j, :, :] = cds_btxw[:, j*binning:(j+1)*binning, :, :].mean(axis=1)\n        return out\n\n    @staticmethod\n    def _flat_field(signal_btxw, flat_wx, dead_wx):\n        \"\"\"\n        signal_btxw: (B, T, W, X)\n        flat_wx, dead_wx: (W, X)\n        \"\"\"\n        signal_btxw = np.asarray(signal_btxw)\n        B, T, W, X = signal_btxw.shape\n        flat2 = np.ma.masked_where(dead_wx, flat_wx)\n        flat2 = np.tile(flat2, (B, T, 1, 1))  # (B,T,W,X)\n        divided = signal_btxw / flat2\n        return np.asarray(divided)  # back to plain ndarray\n\n    # ---------- per-sensor processing ----------\n    def _load_calibration(self, base, sensor):\n        dark = pd.read_parquet(f\"{base}/{sensor}_calibration_0/dark.parquet\").to_numpy()\n        dead = pd.read_parquet(f\"{base}/{sensor}_calibration_0/dead.parquet\").to_numpy()\n        flat = pd.read_parquet(f\"{base}/{sensor}_calibration_0/flat.parquet\").to_numpy()\n        lin  = pd.read_parquet(f\"{base}/{sensor}_calibration_0/linear_corr.parquet\").to_numpy()\n        return dark, dead, flat, lin\n\n    def _load_signal(self, base, sensor):\n        sig = pd.read_parquet(f\"{base}/{sensor}_signal_0.parquet\").to_numpy()\n        if sensor == 'AIRS-CH0':\n            sig = sig.reshape(self.cfg.AIRS_RAW_SHAPE)  # (11250, 32, 356)\n        else:\n            sig = sig.reshape(self.cfg.FGS_RAW_SHAPE)   # (135000, 32, 32)\n        return sig\n\n    def _calibrate_sensor_full(self, planet_id, sensor):\n        base = f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}\"\n\n        # raw signal & calibration\n        signal = self._load_signal(base, sensor)              # (T, 32, Worig)\n        dark, dead, flat, lin = self._load_calibration(base, sensor)\n\n        # reshape cal maps\n        if sensor == 'AIRS-CH0':\n            # maps originally (32,356)\n            dark = dark.reshape(self.cfg.AIRS_RAW_SHAPE[1:])\n            dead = dead.reshape(self.cfg.AIRS_RAW_SHAPE[1:])\n            flat = flat.reshape(self.cfg.AIRS_RAW_SHAPE[1:])\n            lin  = lin.astype(np.float64).reshape(self.cfg.LIN_AIRS_SHAPE)\n\n            # crop wavelengths to 39..320 (→ 282)\n            signal = signal[:, :, self.cfg.CUT_INF:self.cfg.CUT_SUP]  # (T, 32, 282)\n            dark   = dark[:,  self.cfg.CUT_INF:self.cfg.CUT_SUP]      # (32, 282)\n            dead   = dead[:,  self.cfg.CUT_INF:self.cfg.CUT_SUP]\n            flat   = flat[:,  self.cfg.CUT_INF:self.cfg.CUT_SUP]\n            lin    = lin[:, :, self.cfg.CUT_INF:self.cfg.CUT_SUP]\n            base_dt, inc = self.cfg.AIRS_DT_BASE, self.cfg.AIRS_DT_INC\n        else:\n            # FGS maps (32,32)\n            dark = dark.reshape(self.cfg.FGS_RAW_SHAPE[1:])\n            dead = dead.reshape(self.cfg.FGS_RAW_SHAPE[1:])\n            flat = flat.reshape(self.cfg.FGS_RAW_SHAPE[1:])\n            lin  = lin.astype(np.float64).reshape(self.cfg.LIN_FGS_SHAPE)\n            base_dt, inc = self.cfg.FGS_DT_BASE, self.cfg.FGS_DT_INC\n\n        # ADC invert\n        g, o = self._adc_for(sensor)\n        signal = self._adc_convert(signal, g, o)\n\n        # dt pattern\n        dt = self._dt_pattern(sensor, len(signal))\n\n        # clip\n        signal = np.clip(signal, 0, None)\n\n        # mask (returns masked array; convert to ndarray right before CDS)\n        if self.cfg.DO_MASK:\n            signal = self._mask_hot_dead(signal, dead, dark)\n\n        # NL corr\n        if self.cfg.DO_NL_CORR:\n            signal = self._apply_linear_corr(lin, signal)\n\n        # dark subtraction\n        if self.cfg.DO_DARK:\n            signal = self._clean_dark(signal, dead, dark, dt)\n\n        # ensure ndarray (not masked) for CDS\n        signal = np.asarray(signal)\n\n        # CDS: (T, 32, W) → (T/2, 32, W)\n        cds = self._cds(signal)\n\n        # To (B, T, W, X) layout for both sensors\n        # We define W = spectral axis; X = spatial rows\n        if sensor == 'AIRS-CH0':          # cds: (T', 32, 282)\n            cds_btxw = cds[np.newaxis, ...].transpose(0, 1, 3, 2)  # (1, T', W=282, X=32)\n        else:                              # FGS1: (T', 32, 32)\n            # treat W=32, X=32\n            cds_btxw = cds[np.newaxis, ...]                           # (1, T', 32, 32)\n\n        # bin time (keep (B,T,W,X))\n        binning = self.cfg.AIRS_BIN if sensor == 'AIRS-CH0' else self.cfg.FGS_BIN\n        binned = self._time_bin_cds(cds_btxw, binning)  # (1, 187, W, X)\n\n        # flat-field at end (safer numerically after bin)\n        if self.cfg.DO_FLAT:\n            if sensor == 'AIRS-CH0':\n                # flat is (32, 282); need (W=282, X=32)\n                flat_wx = flat.T\n                dead_wx = dead.T\n            else:\n                flat_wx = flat\n                dead_wx = dead\n            binned = self._flat_field(binned, flat_wx, dead_wx)  # (1,187,W,X)\n\n        return np.asarray(binned)  # (1, 187, W, X)\n\n    def _roi_and_aggregate(self, binned, sensor):\n        \"\"\"\n        Reduce spatial dims to (1, 187, K):\n          K = 282 for AIRS-CH0 (mean over spatial X ROI; keep spectral W)\n          K = 1   for FGS1    (mean over 12x12 ROI)\n        \"\"\"\n        binned = np.asarray(binned)\n        B, T, W, X = binned.shape  # B=1\n        if sensor == 'AIRS-CH0':\n            # center rows in X (spatial), mean over them → keep W=282\n            x_roi = binned[:, :, :, self.cfg.ROI_ROWS]              # (1,187,282,12)\n            mean_over_x = np.nanmean(x_roi, axis=3)                 # (1,187,282)\n            return mean_over_x.reshape(1, T, W)\n        else:\n            # FGS1: 12x12 center then mean → scalar per time\n            roi = binned[:, :, self.cfg.ROI_ROWS, self.cfg.ROI_COLS]            # (1,187,12,12)\n            mean_scalar = np.nanmean(roi.reshape(B, T, -1), axis=2)             # (1,187)\n            return mean_scalar.reshape(1, T, 1)\n\n    # ---------- public API ----------\n    def process_all(self):\n        \"\"\"Return (N, 187, 283) by concatenating FGS1(187,1) and AIRS(187,282).\"\"\"\n\n        def run_one(pid, sensor):\n            binned  = self._calibrate_sensor_full(pid, sensor)   # (1, 187, W, X)\n            reduced = self._roi_and_aggregate(binned, sensor)    # (1,187,K)\n            arr = np.asarray(reduced[0])                         # (187,K) or (187,)\n            if arr.ndim == 1:                                    # ensure (187,1)\n                arr = arr.reshape(arr.shape[0], 1)\n            return arr\n\n        # parallel or sequential\n        if PQDM_AVAILABLE:\n            fgs_args  = [dict(pid=p, s='FGS1')     for p in self.planet_ids]\n            airs_args = [dict(pid=p, s='AIRS-CH0') for p in self.planet_ids]\n            fgs_list  = pqdm(fgs_args,  lambda d: run_one(d['pid'], d['s']), n_jobs=self.cfg.N_JOBS)\n            airs_list = pqdm(airs_args, lambda d: run_one(d['pid'], d['s']), n_jobs=self.cfg.N_JOBS)\n        else:\n            fgs_list  = [run_one(p, 'FGS1')     for p in tqdm(self.planet_ids, desc='FGS1')]\n            airs_list = [run_one(p, 'AIRS-CH0') for p in tqdm(self.planet_ids, desc='AIRS-CH0')]\n\n        # detect and raise worker exceptions early (pqdm can return Exceptions)\n        for name, lst in (('FGS1', fgs_list), ('AIRS-CH0', airs_list)):\n            for idx, item in enumerate(lst):\n                if isinstance(item, Exception):\n                    pid = self.planet_ids[idx] if idx < len(self.planet_ids) else None\n                    raise RuntimeError(f\"{name} worker failed for planet_id={pid}\") from item\n\n        # force proper rank before stacking\n        fgs_list  = [a if a.ndim == 2 else a.reshape(a.shape[0], 1) for a in fgs_list]   # (187,1)\n        airs_list = [a if a.ndim == 2 else a.reshape(a.shape[0], 1) for a in airs_list]  # (187,282)\n\n        fgs_arr  = np.stack(fgs_list)   # (N, 187, 1)\n        airs_arr = np.stack(airs_list)  # (N, 187, 282)\n\n        # sanity checks\n        assert fgs_arr.ndim == 3 and fgs_arr.shape[1] == 187 and fgs_arr.shape[2] == 1, \\\n            f\"FGS bad shape: {fgs_arr.shape}\"\n        assert airs_arr.ndim == 3 and airs_arr.shape[1] == 187 and airs_arr.shape[2] >= 100, \\\n            f\"AIRS bad shape: {airs_arr.shape}\"\n\n        return np.concatenate([fgs_arr, airs_arr], axis=2)  # (N, 187, 283)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-07T09:26:38.370164Z","iopub.execute_input":"2025-09-07T09:26:38.370799Z","iopub.status.idle":"2025-09-07T09:26:38.413358Z","shell.execute_reply.started":"2025-09-07T09:26:38.370769Z","shell.execute_reply":"2025-09-07T09:26:38.412338Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# `UnifiedArielPreprocessor`：前処理フルパイプライン）\n\n## 🎛 クラスの目的\n\n`UnifiedArielPreprocessor` は、**AIRS-CH0 / FGS1 の生データを、提出に使える軽量特徴 (N, 187, 283)** に一気通貫で変換します。\n工程は **ADC 逆変換 → マスク →（任意）非線形補正 → ダーク減算 → CDS → 時間ビニング → フラット補正 → 空間縮約 → 結合**。\n\n---\n\n## 🧱 入出力の形状（最重要）\n\n* **入力**\n\n  * AIRS-CH0 生: `(11250, 32, 356)` → **スペクトル切出し** `(T, 32, 282)`\n  * FGS1 生: `(135000, 32, 32)`\n* **中間（CDS → ビニング）**\n\n  * AIRS: `(1, 187, W=282, X=32)`\n  * FGS : `(1, 187, W=32,  X=32)`\n* **空間縮約（ROI 平均）**\n\n  * AIRS: `(1, 187, 282)`（**波長**を保持、**空間**を平均）\n  * FGS : `(1, 187, 1)`（中心 **12×12** を平均）\n* **出力**\n\n  * 結合: **`(N, 187, 283)` = FGS(1) + AIRS(282)**\n\n> ここで **W=波長**, **X=空間行**。AIRS は **W を残し X を平均**、FGS は **12×12 を平均して 1ch** にします。\n\n---\n\n## 🔩 主なメソッドと役割\n\n### `_load_signal / _load_calibration`\n\n* 各惑星・各センサーの **signal / dark / dead / flat / linear\\_corr** を Parquet から読み込み。\n* AIRS は **波長 39..320（排他上限321）** に切り出して **282 チャンネル**にします。\n\n### `_adc_convert`（ADC 逆変換）\n\n* ADU → アナログ相当へ。`adc_info.csv` があれば**実測ゲイン/オフセット**を、なければ `Config` の**既定値**を使用。\n\n### `_mask_hot_dead`\n\n* **dead** はそのままマスク。**hot** はダークを `sigma_clip` で抽出してマスク。\n* NumPy の **masked array** を返します（後続で `np.asarray` へ戻して安定化）。\n\n### `_apply_linear_corr`（任意）\n\n* 画素ごとの **逆多項式係数**で **非線形応答**を補正。\n* まずは OFF（`DO_NL_CORR=False`）で安定動作を確認 → 効果検証の段階で ON がおすすめ。\n\n### `_clean_dark`\n\n* **積分時間 `dt`**（偶奇で `base + inc`）に合わせて **dark × dt** を減算。\n* 加法バイアスを正しく除くことで、**相対比（トランジット深さ）の分母の歪み**を避けます。\n\n### `_cds`\n\n* 読み出し **開始/終了の差（終−初）** を取ることで、**オフセットやゆっくりしたドリフト**を打消し。\n* **3D/4D 両対応**（(T,W,X) or (B,T,W,X)）で安全。\n\n### `_time_bin_cds`\n\n* 時間方向を **非重複平均**。AIRS=30, FGS=360（=30×12）で **\\~187 ステップ**に縮約。\n* 端数は切り捨て（提出形状を安定確保）。\n\n### `_flat_field`\n\n* **(W, X)** 形状のフラットを **デッドでマスク**し、(B,T,W,X) 全体にタイリングして除算。\n* AIRS は `flat.T`（(32,282) → (282,32)）で **W,X の向きを一致**させます。\n* **ビニング後**に適用することで、**数値安定**＆**欠損伝播の抑制**。\n\n### `_roi_and_aggregate`\n\n* AIRS：**中心 12 行**（Config の `ROI_ROWS`）平均 → **波長 282 を保持**。\n* FGS：**中心 12×12**（`ROI_ROWS`×`ROI_COLS`）平均 → **1 チャンネル化**。\n* どちらも **(1, 187, K)** に正規化して返します（K=282 or 1）。\n\n### `process_all`（公開 API）\n\n* 全惑星に対して、AIRS/FGS を \\*\\*順次 or 並列（pqdm）\\*\\*で処理。\n* 並列実行では **例外オブジェクト**が返る場合があるため、**検出して惑星ID付きで即時再スロー**。\n* 最後に **`np.stack → np.concatenate(axis=2)`** で **(N, 187, 283)** を構築。\n\n---\n\n## 🛡️ 安全装置（落とし穴対策）\n\n* **AIRS の transpose**：`(T,32,282) → (1,T,282,32)` の順番を**固定**（重複軸指定のミスを排除）。\n* **masked array → ndarray**：CDS 前に `np.asarray` で**強制変換**。\n* **形状アサート**：結合前に `(N,187,1)` / `(N,187,≥100)` をチェックして**早期に形状バグを発見**。\n* **例外の見える化**：`RuntimeError(f\"{sensor} worker failed for planet_id=...\")` で**データ側の異常箇所を特定**。\n\n---\n\n## 🧪 検証の観点\n\n* AIRS/FGS それぞれで **CDS の平均が \\~0 付近**か（大きなオフセット残りがないか）。\n* フラット適用後、**空間方向のばらつき**が減っているか。\n* ROI 平均後、**残差の時間相関**が減少しているか。\n* 出力 `(N,187,283)` の **NaN 率**や**外れ値**が過大でないか。\n\n---\n\n## 🔄 カスタマイズ例\n\n* **ROI を広げる/狭める**：`ROI_ROWS`, `ROI_COLS` を変更（S/N vs. 系統のバランス）。\n* **ビニング幅の調整**：`AIRS_BIN`, `FGS_BIN` を同時に変更して**187 付近を維持**。\n* **非線形補正の有効化**：`DO_NL_CORR=True`。ただし **係数の信頼領域**（飽和近傍など）を要確認。\n* **フラットの順序**：理論的には前後で等価な場合もありますが、**本実装ではビニング後**が安定。\n\n---\n\n","metadata":{}},{"cell_type":"code","source":"class TransitModel:\n    def __init__(self, config: Config):\n        self.cfg = config\n\n    def _phase_detector(self, signal_1d):\n        sl = self.cfg.MODEL_PHASE_DETECTION_SLICE\n        min_idx = int(np.argmin(signal_1d[sl]) + sl.start)\n\n        pre = signal_1d[:min_idx]\n        post = signal_1d[min_idx:]\n\n        g1 = np.gradient(pre);  g1 /= max(g1.max(), 1e-12)\n        g2 = np.gradient(post); g2 /= max(g2.max(), 1e-12)\n\n        phase1 = int(np.argmin(g1))\n        phase2 = int(np.argmax(g2) + min_idx)\n        return phase1, phase2\n\n    def _objective(self, s, signal, phase1, phase2):\n        delta = self.cfg.MODEL_OPTIMIZATION_DELTA\n        power = self.cfg.MODEL_POLYNOMIAL_DEGREE\n\n        if phase1 - delta <= 0 or phase2 + delta >= len(signal) or (phase2 - delta) - (phase1 + delta) < 5:\n            delta = 2\n\n        y = np.concatenate([\n            signal[:phase1 - delta],\n            signal[phase1 + delta:phase2 - delta] * (1 + s),\n            signal[phase2 + delta:]\n        ])\n        x = np.arange(len(y))\n        coeffs = np.polyfit(x, y, deg=power)\n        poly = np.poly1d(coeffs)\n        return np.mean(np.abs(poly(x) - y))\n\n    def predict(self, single_preprocessed_signal):\n        # drop FGS1 column, average over wavelengths → (187,)\n        signal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)\n        signal_1d = savgol_filter(signal_1d, 20, 2)\n\n        p1, p2 = self._phase_detector(signal_1d)\n        p1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, p1)\n        p2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, p2)\n\n        res = minimize(\n            fun=self._objective,\n            x0=[0.0001],\n            args=(signal_1d, p1, p2),\n            method='Nelder-Mead'\n        )\n        return float(res.x[0])\n\n    def predict_all(self, preprocessed_signals):\n        preds = [self.predict(sp) for sp in tqdm(preprocessed_signals, desc='Model')]\n        return np.array(preds) * self.cfg.SCALE\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-07T09:26:38.414496Z","iopub.execute_input":"2025-09-07T09:26:38.414783Z","iopub.status.idle":"2025-09-07T09:26:38.432902Z","shell.execute_reply.started":"2025-09-07T09:26:38.414762Z","shell.execute_reply":"2025-09-07T09:26:38.431903Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# `TransitModel`：トランジット深さ $s$ の最小化推定\n\n## 🎯 目的\n\n前処理で得た **(N, 187, 283)** の特徴から、各惑星について **1 つのトランジット深さ $s$** を推定します。\nここでの方針は **“非 ML の最小自由度”**：\n\n* \\*\\*トランジット区間（in-transit）\\*\\*を見つける\n* その区間の信号を $(1+s)$ 倍して**ディップ（落ち込み）を埋める**\n* すると、全区間は**スムーズ**につながるはず\n* **“全体が最も滑らかになる $s$”** を最適化で選ぶ\n\n---\n\n## 🧠 理解のコア（直感）\n\n* 光度曲線 $L(t)$ はトランジット中だけ**わずかに下がる**。\n* もし真の深さ $s^\\*$ を知っていて、in-transit 部分を $L(t)\\times (1+s^\\*)$ に持ち上げたら、\n  **カーブの前後（out-of-transit）と段差が消え、全体は滑らか**になります。\n* 逆に、間違った $s$ だと**段差（不連続）や歪み**が残ります。\n* その“滑らかさ”を、\\*\\*低次数多項式（ここでは 3 次）\\*\\*の当たり具合（残差）で測ります。\n\n---\n\n## 🔎 ステップ内訳\n\n### 1) 1次元化と平滑化\n\n```python\nsignal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)  # (187,), FGS(1列)は除き、AIRS 282列の平均\nsignal_1d = savgol_filter(signal_1d, 20, 2)                  # Savitzky–Golayで軽くスムージング\n```\n\n* **FGS 列（1 列）を除外**して、**AIRS の 282 波長を平均**した時系列を使います（ベースライン想定）。\n* Savitzky–Golay でノイズを抑え、**位相検出を安定**させています。\n\n  * ウィンドウ=20, 次数=2 は**過平滑しすぎず**にトレンド把握を狙う設定。\n\n> 拡張案：AIRS を重み付き平均したり、FGS をトレンド補正に用いたり、\n> 波長ごとに $s(\\lambda)$ を出すなども可能。\n\n---\n\n### 2) 位相検出（`_phase_detector`）\n\n```python\nsl = self.cfg.MODEL_PHASE_DETECTION_SLICE        # 例: slice(30,140)\nmin_idx = argmin(signal_1d[sl]) + sl.start       # まずディップ底の近傍を見つける\npre  = signal_1d[:min_idx]\npost = signal_1d[min_idx:]\n\ng1 = gradient(pre);  g1 /= max(g1.max(), 1e-12)  # 前半の下降エッジ\ng2 = gradient(post); g2 /= max(g2.max(), 1e-12)  # 後半の上昇エッジ\n\nphase1 = argmin(g1)                              # 下がり始め（最も強い負勾配付近）\nphase2 = argmax(g2) + min_idx                    # 上がり切り（最も強い正勾配付近）\n```\n\n* ディップの**底**（最小値）を中心に、前半/後半の**勾配**を見て **下降開始/上昇終了**を推定。\n* これで **in-transit 区間 $[phase1, phase2]$** が粗く決まります。\n\n> うまくいかない時のシグナル：\n>\n> * 勾配がノイジー → `savgol_filter` の窓やスライス `slice(30,140)` を調整\n> * 2峰性などで誤検出 → **探索スライス**を物理的に妥当な範囲へ絞る\n\n---\n\n### 3) 目的関数（`_objective`）\n\n```python\n# 端の保護（in-transitの両端からdelta点は使わない）\ndelta = self.cfg.MODEL_OPTIMIZATION_DELTA  # 例: 7\nif (phase1 - delta <= 0) or (phase2 + delta >= len(signal)) or ((phase2 - delta) - (phase1 + delta) < 5):\n    delta = 2\n\n# in-transit区間だけ (1+s) 倍に“持ち上げる”\ny = concat([\n    signal[:phase1 - delta],\n    signal[phase1 + delta : phase2 - delta] * (1 + s),\n    signal[phase2 + delta :]\n])\nx = arange(len(y))\n\n# 低次数多項式（3次）を当てたときの絶対誤差の平均\ncoeffs = polyfit(x, y, deg=power)  # power=3\npoly   = poly1d(coeffs)\nerror  = mean(abs(poly(x) - y))\n```\n\n* **端の $\\delta$** を捨てるのは、**エッジ効果や位相ずれ**の影響を減らすため。\n* in-transit を $(1+s)$ 倍し、**全体に 3 次多項式**をフィット。\n* **残差の平均絶対値**をエラーに採用 → ステップ（不連続）を嫌う性質にマッチ。\n\n  * L2（MSE）でなく L1（MAE）気味にするのは、**外れ値に少し強く**する意図。\n\n---\n\n### 4) 最適化（`predict`）\n\n```python\nres = minimize(\n    fun=self._objective,\n    x0=[1e-4],                      # 初期値 s≈0\n    args=(signal_1d, p1, p2),\n    method='Nelder-Mead'            # 派生なしのシンプルなシンプレックス\n)\ns_hat = float(res.x[0]) * self.cfg.SCALE  # スケール係数で微調整\n```\n\n* **Nelder–Mead** は導関数不要・ロバストで、**1 変数**最適化に最適。\n* 初期値は小さめ（ $10^{-4}$ 程度）で十分。\n* `SCALE` は全体的な**軽い補正ノブ**（例：0.95）。バリデーションで調整。\n\n---\n\n## 🧪 品質チェック（model 側）\n\n* **位相の妥当性**：`phase1 < phase2`、かつ in-transit 長が短すぎないか（`delta` で守る）。\n* **エラー曲線の形**：`s` 近傍で**凸**に見えるか（極小がはっきり）。\n* **平滑化の過不足**：平滑化しすぎると深さが**過小**、しなさすぎると位相が**不安定**。\n* **FGS の扱い**：今回は除外して AIRS 平均のみ。FGS を**トレンド除去**に使うのも有効。\n\n---\n\n## 🔄 発展（この枠組みを保ったまま強化）\n\n* **波長依存 $s(\\lambda)$**：AIRS 282 次元に対して、**スプライン等で滑らかさ正則化**を入れて同じ目的を最小化。\n* **物理モデルの併用**：非線形ベースライン（姿勢/温度）や**システマティック回帰**を同時に入れる。\n* **ロバスト損失/Huber**：外れ値耐性をさらに上げる。\n* **位相と深さの同時推定**：位相と $s$ を**同時に最適化**、あるいは**ベイズ推定**で事後分布へ。\n\n---\n\n## 📝 まとめ\n\n* **情報最小**の 1D 光度曲線から、**ディップを“埋める”という幾何学的直感**で $s$ を推定。\n* **位相検出 → in-transit の $(1+s)$ スケーリング → 低次数多項式の滑らかさ** という\n  **シンプルだが強いバイアス**を使い、**計算は軽く・解釈は明快**に保っています。\n* このベースラインは、**後で AIRS の波長情報を細かく使う拡張**や、**物理的トレンドの同時回帰**に自然に接続できます。\n\n次の **Cell 5** では、推定した $s$ を提出フォーマットに合わせて展開し、`submission.csv` を作成します。\n","metadata":{}},{"cell_type":"code","source":"class SubmissionGenerator:\n    def __init__(self, config: Config):\n        self.cfg = config\n        self.sample_submission = pd.read_csv(\n            f\"{self.cfg.DATA_PATH}/sample_submission.csv\",\n            index_col=\"planet_id\"\n        )\n\n    def create(self, predictions):\n        planet_ids = self.sample_submission.index\n        K = self.sample_submission.shape[1] // 2  # half are sigma columns\n        vals = np.repeat(predictions, K).reshape(len(predictions), -1)\n        vals = np.clip(vals, 0, None)\n        sigmas = np.ones_like(vals) * self.cfg.SIGMA\n\n        df = pd.DataFrame(\n            np.concatenate([vals, sigmas], axis=1),\n            columns=self.sample_submission.columns,\n            index=planet_ids\n        )\n        df.to_csv(\"submission.csv\")\n        return df\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-07T09:26:38.43508Z","iopub.execute_input":"2025-09-07T09:26:38.435411Z","iopub.status.idle":"2025-09-07T09:26:38.457242Z","shell.execute_reply.started":"2025-09-07T09:26:38.435377Z","shell.execute_reply":"2025-09-07T09:26:38.456377Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"\n\n## 📤 `SubmissionGenerator`：提出ファイルの生成\n\n### 何をしているか\n\n* Kaggle が配布する **`sample_submission.csv`** を読み込み、**列構成（順序・本数）** をそのまま踏襲します。\n* 予測した **トランジット深さ $s$** を、提出フォーマットの\\*\\*「深さ列」全てに複製\\*\\*して埋め込みます。\n  （今回のベースラインは **単一スカラー**なので、**すべての波長に同じ値**を入れます）\n* シグマ列は**定数 `SIGMA`** を一律で設定（ベースライン）。\n* 最終的に **`submission.csv`** を出力します。\n\n### 実装のポイント\n\n* `K = sample_submission.shape[1] // 2`\n  → 右半分が $\\sigma$ 列なので、左半分（K 列ぶん）に $s$ を複製します。\n* **安全策**：`np.clip(vals, 0, None)` で負値の深さを 0 に切り上げ（提出の健全性を担保）。\n* 将来拡張（例）\n\n  * **波長依存の $s(\\lambda)$** を学習したら、**AIRS の 282 列**に対応させて列ごとに値を入れる。\n  * $\\sigma$ を**推定由来**に変えて波長ごとに異なる不確かさを反映。\n\n---","metadata":{}},{"cell_type":"code","source":"cfg = Config()\nprep = UnifiedArielPreprocessor(cfg)\n\n# unified preprocessing → (N, 187, 283)\npreprocessed = prep.process_all()\nprint('Preprocessed shape:', preprocessed.shape)\n\n# model\nmodel = TransitModel(cfg)\npreds = model.predict_all(preprocessed)\n\n# submission\nsub = SubmissionGenerator(cfg).create(preds)\nsub.head()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-07T09:26:38.458229Z","iopub.execute_input":"2025-09-07T09:26:38.45853Z","iopub.status.idle":"2025-09-07T09:26:52.030723Z","shell.execute_reply.started":"2025-09-07T09:26:38.458509Z","shell.execute_reply":"2025-09-07T09:26:52.029828Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 🚀 実行ブロック：一気通貫で実行\n\n### 処理フロー\n\n1. `cfg = Config()`：設定を一括管理。\n2. `prep = UnifiedArielPreprocessor(cfg)`：前処理クラスのインスタンス化。\n3. `preprocessed = prep.process_all()`：\n\n   * **(N, 187, 283)** のテンソルを得る（FGS1=1 + AIRS=282）。\n   * 内部で **ADC → マスク →（任意）非線形 → ダーク → CDS → ビニング → フラット → ROI 平均** を実施。\n4. `model = TransitModel(cfg); preds = model.predict_all(preprocessed)`：\n\n   * 各惑星について **スカラー $s$** を最適化で推定（**段差が最小になる $s$**）。\n5. `sub = SubmissionGenerator(cfg).create(preds)`：\n\n   * 提出ファイル **`submission.csv`** を生成。`sub.head()` で先頭を確認。\n\n### 実行時の確認ポイント\n\n* `print('Preprocessed shape:', preprocessed.shape)` が **(N, 187, 283)** になっている。\n* 処理が重い場合は `cfg.N_JOBS` を下げるか、`pqdm` なし（逐次）で動作確認 → その後並列化。\n* 例外が出たら、メッセージに出る **planet\\_id** を手がかりに該当データを個別デバッグ。\n\n---\n\n## ✅ 結論（まとめ）\n\n* **本ノートの核**は、Ariel データに対する **科学的に素直な前処理**と、**シンプルで解釈容易な推定器**です。\n\n  * **前処理**：検出器の系統（ダーク・フラット・非線形・死/ホット）をケアし、**CDS とビニングで S/N を改善**。\n  * **特徴設計**：AIRS は**波長情報（282）を保持**し、空間方向のみを平均。FGS は**参照光**として 1ch 化。\n  * **推定**：トランジット区間を $(1+s)$ 倍して**全体を滑らかにする $s$** を選ぶ、**幾何学的でロバストな**非 ML ベースライン。\n\n* **強み**\n\n  * 実装が **短く・速く・再現性が高い**。\n  * 物理の観点で**何をどの順で補正しているかが明確**。\n  * そのまま **$s(\\lambda)$** や**トレンド同時回帰**などへ**拡張しやすい設計**。\n\n* **今後の伸ばし方**\n\n  1. **波長依存 $s(\\lambda)$**：AIRS 282 チャンネルへ**滑らかさ正則化**を入れたマルチ出力化。\n  2. **FGS の活用強化**：AIRS のトレンド補償（共通モード除去）や説明変数として併用。\n  3. **ロバスト損失 / Huber / Quantile**：外れ値耐性・非対称ノイズへの強化。\n  4. **ベイズ的枠組み**：位相・深さ・ベースラインを**同時に推定**し、**不確かさの伝播**を厳密化。\n  5. **高速化**：I/O 並列、メモリ節約、JIT（numba）や GPU 化で**前処理コスト**を削減。\n\n* **参考ノート**（再掲）\n\n  * Calibrating and Binning Ariel Data（前処理の実務レシピ）\n    [https://www.kaggle.com/code/gordonyip/calibrating-and-binning-ariel-data](https://www.kaggle.com/code/gordonyip/calibrating-and-binning-ariel-data)\n  * NeurIPS Non-ML Transit Curve Fitting（非 ML 最適化の思想）\n    [https://www.kaggle.com/code/vitalykudelya/neurips-non-ml-transit-curve-fitting](https://www.kaggle.com/code/vitalykudelya/neurips-non-ml-transit-curve-fitting)\n\n> このベースラインは「**堅い土台**」として、あなたの工夫（特徴量設計・目的関数・正則化・ベイズ化）を\n> 上に積みやすい構造になっています。ここから一緒に**精度をガンガン伸ばして**いきましょう！ 🌟\n","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}