{"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"},{"sourceId":9629432,"sourceType":"datasetVersion","datasetId":5846888},{"sourceId":13123369,"sourceType":"datasetVersion","datasetId":8313229},{"sourceId":13121320,"sourceType":"datasetVersion","datasetId":8304548},{"sourceId":564250,"sourceType":"modelInstanceVersion","modelInstanceId":426187,"modelId":443655},{"sourceId":564252,"sourceType":"modelInstanceVersion","modelInstanceId":426188,"modelId":443656},{"sourceId":564240,"sourceType":"modelInstanceVersion","isSourceIdPinned":true,"modelInstanceId":426180,"modelId":443648}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Solution1","metadata":{}},{"cell_type":"code","source":"# install pqdm for parallel processing\n!pip install --no-index --find-links=/kaggle/input/ariel-2024-pqdm pqdm","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:22.827122Z","iopub.execute_input":"2025-09-24T08:28:22.827438Z","iopub.status.idle":"2025-09-24T08:28:29.388255Z","shell.execute_reply.started":"2025-09-24T08:28:22.827415Z","shell.execute_reply":"2025-09-24T08:28:29.386809Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport torch\nimport torch.nn.functional as F\nimport multiprocessing as mp\nimport torch.nn as nn\nimport os\nimport matplotlib.pyplot as plt\nimport itertools\n\nfrom tqdm import tqdm\nfrom astropy.stats import sigma_clip\nfrom scipy.optimize import minimize\nfrom torch.utils.data import DataLoader, TensorDataset, random_split\nfrom sklearn.preprocessing import StandardScaler\nfrom scipy.signal import savgol_filter\nfrom astropy.stats import sigma_clip","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:29.390657Z","iopub.execute_input":"2025-09-24T08:28:29.391365Z","iopub.status.idle":"2025-09-24T08:28:37.503432Z","shell.execute_reply.started":"2025-09-24T08:28:29.391319Z","shell.execute_reply":"2025-09-24T08:28:37.501683Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\n\nfrom tqdm import tqdm\nfrom pqdm.threads import pqdm\nimport itertools\n\nfrom scipy.optimize import minimize\nfrom sklearn.metrics import mean_squared_error\n\nfrom astropy.stats import sigma_clip\nfrom scipy.signal import savgol_filter\n\nimport time\n__t0 = time.perf_counter()\nROOT_PATH = \"/kaggle/input/ariel-data-challenge-2025\"\nMODE = \"test\"\n\nclass Config:\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n\n    SCALE = 0.95\n    SIGMA = 0.0009\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(28, 145)\n    MODEL_OPTIMIZATION_DELTA = 11\n    MODEL_POLYNOMIAL_DEGREE = 3\n    \n    N_JOBS = 3\n\ndef _phase_detector_signal(signal, cfg):\n    sl = cfg.MODEL_PHASE_DETECTION_SLICE\n    min_idx = int(np.argmin(signal[sl])) + sl.start\n    s1 = signal[:min_idx]; s2 = signal[min_idx:]\n    if s1.size < 3 or s2.size < 3:\n        return 0, len(signal) - 1\n    g1 = np.gradient(s1); g1_max = np.max(g1) if np.size(g1) else 0.0\n    g2 = np.gradient(s2); g2_max = np.max(g2) if np.size(g2) else 0.0\n    if g1_max != 0: g1 /= g1_max\n    if g2_max != 0: g2 /= g2_max\n    phase1 = int(np.argmin(g1)); phase2 = int(np.argmax(g2)) + min_idx\n    return phase1, phase2\n\ndef estimate_sigma_fgs(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_1 (для FGS1) длиной N_planets — мягкий множитель к cfg.SIGMA.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n    for single in preprocessed_data:\n        # фазы по AIRS белой кривой — так же, как в модели\n        air_white = savgol_filter(single[:, 1:].mean(axis=1), 20, 2)\n        p1, p2 = _phase_detector_signal(air_white, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(air_white) - delta - 1, p2)\n\n        fgs = single[:, 0]\n        oot = (fgs[: p1 - delta] if p1 - delta > 0 else np.empty(0, fgs.dtype))\n        if p2 + delta < fgs.size:\n            oot = np.concatenate([oot, fgs[p2 + delta :]])\n        inn = fgs[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(fgs))\n        # относительная неопределённость глубины (в тех же ед., что s)\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель: корень, и узкий клип, чтобы не рисковать\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.8, 1.25)  # ±20–25% от базовой σ\n\n    return k * cfg.SIGMA\n\ndef estimate_sigma_air(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_air длиной N_planets — мягкий множитель к cfg.SIGMA для всех AIRS-каналов.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n\n    for single in preprocessed_data:\n        # белая кривая AIRS на бинированных данных (после всех твоих весов по λ)\n        white = np.nanmean(single[:, 1:], axis=1)         # (n_bins,)\n        white_s = savgol_filter(white, 20, 2)             # для фаз\n\n        p1, p2 = _phase_detector_signal(white_s, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(white) - delta - 1, p2)\n\n        oot_left = white[: p1 - delta] if p1 - delta > 0 else np.empty(0, white.dtype)\n        oot_right = white[p2 + delta :] if (p2 + delta) < white.size else np.empty(0, white.dtype)\n        oot = np.concatenate([oot_left, oot_right]) if (oot_left.size + oot_right.size) else oot_left\n        inn = white[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(white))\n\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель вокруг медианы\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.90, 1.20)  # ±10%–20%\n\n    return k * cfg.SIGMA\n\n\nclass SignalProcessor:\n    def __init__(self, config):\n        self.cfg = config\n        self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\n        self.planet_ids = pd.read_csv(f'{self.cfg.DATA_PATH}/{self.cfg.DATASET}_star_info.csv', index_col='planet_id').index.astype(int)\n\n    def _apply_linear_corr(self, linear_corr, signal):\n\n        coeffs = np.flip(linear_corr, axis=0)      # shape: (D, X, Y), D — старшая степень сначала\n        x = signal.astype(np.float64, copy=False)  # считаем в float64 для стабильности\n        out = np.empty_like(x, dtype=np.float64)\n        out[...] = coeffs[0]  # broadcast (X,Y) -> (T,X,Y)\n        for k in range(1, coeffs.shape[0]):\n            np.multiply(out, x, out=out)  # in-place умножение\n            out += coeffs[k]              # broadcast (X,Y)\n\n        return out.astype(signal.dtype, copy=False)\n\n    def _calibrate_single_signal(self, planet_id, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n\n        signal = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_signal_0.parquet\").to_numpy()\n        dark = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dark.parquet\").to_numpy()\n        dead = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dead.parquet\").to_numpy()\n        flat = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/flat.parquet\").to_numpy()\n        linear_corr = pd.read_parquet(f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\").values.astype(np.float64).reshape(sensor_cfg[\"linear_corr_shape\"])\n\n        signal = signal.reshape(sensor_cfg[\"raw_shape\"])\n        gain = self.adc_info[f\"{sensor}_adc_gain\"].iloc[0]\n        offset = self.adc_info[f\"{sensor}_adc_offset\"].iloc[0]\n        signal = signal / gain + offset\n\n        hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n\n        if sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            linear_corr = linear_corr[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dark = dark[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dead = dead[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            flat = flat[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            hot = hot[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n\n\n        if sensor == \"FGS1\":\n            y0, y1, x0, x1 = 10, 22, 10, 22\n            signal = signal[:, y0:y1, x0:x1]\n            dark   = dark[y0:y1, x0:x1]\n            dead   = dead[y0:y1, x0:x1]\n            flat   = flat[y0:y1, x0:x1]\n            linear_corr = linear_corr[:, y0:y1, x0:x1]\n            hot    = hot[y0:y1, x0:x1]\n\n        #np.maximum(signal, 0, out=signal)\n\n        if sensor == \"FGS1\":\n            signal = self._apply_linear_corr(linear_corr, signal)\n        elif sensor == \"AIRS-CH0\":\n            sl = (slice(None), slice(10, 22), slice(None))     # T, Y, λ (уже порезан по λ выше)\n            signal[sl] = self._apply_linear_corr(linear_corr[:, 10:22, :], signal[sl])\n        else:\n            signal = self._apply_linear_corr(linear_corr, signal)\n\n\n\n\n        # === стало (in‑place, без 3D-временного буфера) ===\n        base_dt, increment = sensor_cfg[\"dt_pattern\"]\n        even_scale = base_dt\n        odd_scale  = base_dt + increment\n\n        signal[::2] -= dark * even_scale   # кадры 0,2,4,...\n        signal[1::2] -= dark * odd_scale   # кадры 1,3,5,...\n\n\n        # --- APPLY FLAT (ROI-aware, без reshape) ---\n        if sensor == \"FGS1\":\n            # после раннего ROI: signal, dark, dead, hot, flat уже (12,12)\n            flat_roi = flat.astype(signal.dtype, copy=False).copy()   # (12,12)\n            mask_roi = (dead | hot)                                   # (12,12)\n            # защита от нулей/NaN в flat\n            flat_roi[mask_roi | ~np.isfinite(flat_roi) | (flat_roi == 0)] = np.nan\n            signal /= flat_roi                                        # (T,12,12) / (12,12)\n\n        elif sensor == \"AIRS-CH0\":\n            # делим только в рабочем ROI по Y (10:22), весь спектр по λ\n            y0, y1 = 10, 22\n            flat_roi = flat[y0:y1, :].astype(signal.dtype, copy=False).copy()   # (12, λ)\n            mask_roi = (dead | hot)[y0:y1, :]                                   # (12, λ)\n            flat_roi[mask_roi | ~np.isfinite(flat_roi) | (flat_roi == 0)] = np.nan\n            signal[:, y0:y1, :] /= flat_roi                                     # (T,12,λ) / (12,λ)\n\n        else:\n            # запасной вариант для прочих сенсоров с «полным кадром»\n            flat2 = flat.astype(signal.dtype, copy=False).copy()\n            mask2 = (dead | hot)\n            flat2[mask2 | ~np.isfinite(flat2) | (flat2 == 0)] = np.nan\n            signal /= flat2\n        # --- END FLAT ---\n        \n        return signal\n\n    def _preprocess_calibrated_signal(self, calibrated_signal, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n        binning = sensor_cfg[\"binning\"]\n\n        if sensor == \"AIRS-CH0\":\n            signal_roi = calibrated_signal[:, 10:22, :]\n        elif sensor == \"FGS1\":\n            signal_roi = calibrated_signal[:, 10:22, 10:22]\n            signal_roi = signal_roi.reshape(signal_roi.shape[0], -1)\n        \n        mean_signal = np.nanmean(signal_roi, axis=1)\n\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        n_bins = cds_signal.shape[0] // binning\n        binned = np.array([\n            cds_signal[j*binning : (j+1)*binning].mean(axis=0) \n            for j in range(n_bins)\n        ])\n\n        # >>> НОВОЕ: винсоризация ПОСЛЕ биннинга (дёшево), только для AIRS\n        if sensor == \"AIRS-CH0\":\n            q_lo = np.nanpercentile(binned, 5.0, axis=1, keepdims=True)    # (n_bins, 1)\n            q_hi = np.nanpercentile(binned, 95.0, axis=1, keepdims=True)   # (n_bins, 1)\n            np.clip(binned, q_lo, q_hi, out=binned)\n\n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        if sensor == \"AIRS-CH0\":\n            # Инверсно-дисперсные веса по λ на бinned-рядe (n_bins, λ)\n            var = np.nanvar(binned, axis=0, ddof=1)                 # (λ, )\n            med = np.nanmedian(var)\n            # заменим невалидные/слишком маленькие дисперсии на медиану\n            safe_var = np.where(~np.isfinite(var) | (var <= 0), med if (np.isfinite(med) and med > 0) else 1.0, var)\n            w = 1.0 / safe_var\n\n            # защитный клип весов, чтобы один канал не доминировал\n            lo, hi = np.nanpercentile(w, 5.0), np.nanpercentile(w, 95.0)\n            if np.isfinite(lo) and np.isfinite(hi) and lo < hi:\n                w = np.clip(w, lo, hi)\n\n            # нормировка: сумма весов = числу каналов → mean == взвешенному mean\n            M = binned.shape[1]\n            s = np.nansum(w)\n            if np.isfinite(s) and s > 0:\n                w = w * (M / s)\n            else:\n                w = np.ones_like(w)\n\n            # применяем веса к каждому времени (broadcast по оси 0)\n            binned *= w[None, :]\n\n\n        return binned\n\n    def _process_planet_sensor(self, args):\n        planet_id, sensor = args['planet_id'], args['sensor']\n        calibrated = self._calibrate_single_signal(planet_id, sensor)\n        preprocessed = self._preprocess_calibrated_signal(calibrated, sensor)\n        return preprocessed\n\n    def process_all_data(self):\n        args_fgs1 = [dict(planet_id=planet_id, sensor=\"FGS1\") for planet_id in self.planet_ids]\n        preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        args_airs_ch0 = [dict(planet_id=planet_id, sensor=\"AIRS-CH0\") for planet_id in self.planet_ids]\n        preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        preprocessed_signal = np.concatenate(\n            [np.stack(preprocessed_fgs1), np.stack(preprocessed_airs_ch0)], axis=2\n        )\n        return preprocessed_signal\n    \n\nclass TransitModel:\n    def __init__(self, config):\n        self.cfg = config\n\n    def _phase_detector(self, signal):\n        search_slice = self.cfg.MODEL_PHASE_DETECTION_SLICE\n        min_index = np.argmin(signal[search_slice]) + search_slice.start\n        \n        signal1 = signal[:min_index]\n        signal2 = signal[min_index:]\n\n        grad1 = np.gradient(signal1)\n        grad1 /= grad1.max()\n        \n        grad2 = np.gradient(signal2)\n        grad2 /= grad2.max()\n\n        phase1 = np.argmin(grad1)\n        phase2 = np.argmax(grad2) + min_index\n\n        return phase1, phase2\n    \n    def _objective_function(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\n        coeffs = np.polyfit(x, y, deg=power)\n        poly = np.poly1d(coeffs)\n        error = np.abs(poly(x) - y).mean()\n        \n        return error\n\n    def predict(self, single_preprocessed_signal):\n        signal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)\n        signal_1d = savgol_filter(signal_1d, 23, 2)\n        \n        phase1, phase2 = self._phase_detector(signal_1d)\n\n        phase1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, phase1)\n        phase2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, phase2)    \n\n        result = minimize(\n            fun=self._objective_function,\n            x0=[0.0001],\n            args=(signal_1d, phase1, phase2),\n            method=\"Nelder-Mead\"\n        )\n        \n        return result.x[0]\n\n    def predict_all(self, preprocessed_signals):\n        predictions = [\n            self.predict(preprocessed_signal)\n            for preprocessed_signal in tqdm(preprocessed_signals)\n        ]\n        return np.array(predictions) * self.cfg.SCALE\n    \nclass SubmissionGenerator:\n    def __init__(self, config):\n        self.cfg = config\n        self.sample_submission = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", index_col=\"planet_id\")\n\n    def create(self, predictions, sigma_fgs=None, sigma_air=None):\n        planet_ids = self.sample_submission.index\n        n_mu = self.sample_submission.shape[1] // 2  # 283\n\n        preds = np.asarray(predictions, dtype=float).reshape(-1)\n        mu = np.tile(preds.reshape(-1, 1), (1, n_mu))\n        mu = np.clip(mu, 0, None)\n\n        sigmas = np.full_like(mu, self.cfg.SIGMA, dtype=float)\n        if sigma_fgs is not None:\n            sigma_fgs = np.asarray(sigma_fgs, dtype=float).reshape(-1)\n            sigmas[:, 0] = np.clip(sigma_fgs, 1e-6, 0.1)\n        if sigma_air is not None:\n            sigma_air = np.asarray(sigma_air, dtype=float).reshape(-1, 1)\n            sigmas[:, 1:] = np.clip(sigma_air, 1e-6, 0.1)\n\n        submission_df = pd.DataFrame(\n            np.concatenate([mu, sigmas], axis=1),\n            columns=self.sample_submission.columns,\n            index=planet_ids\n        )\n        #submission_df.to_csv(\"submission.csv\")\n        return submission_df\n\n\n\nconfig = Config()\n    \nsignal_processor = SignalProcessor(config)\npreprocessed_data = signal_processor.process_all_data()\n\nmodel = TransitModel(config)\npredictions = model.predict_all(preprocessed_data)\nsigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)  # новый шаг\nsigma_air_vec = estimate_sigma_air(preprocessed_data, config)\n\nsubmission_generator = SubmissionGenerator(config)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:37.504393Z","iopub.execute_input":"2025-09-24T08:28:37.504796Z","iopub.status.idle":"2025-09-24T08:28:44.792692Z","shell.execute_reply.started":"2025-09-24T08:28:37.504773Z","shell.execute_reply":"2025-09-24T08:28:44.7915Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_solution1 = submission_generator.create(predictions, sigma_fgs=sigma_fgs_vec, sigma_air=sigma_air_vec)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:44.795447Z","iopub.execute_input":"2025-09-24T08:28:44.795733Z","iopub.status.idle":"2025-09-24T08:28:44.80167Z","shell.execute_reply.started":"2025-09-24T08:28:44.79571Z","shell.execute_reply":"2025-09-24T08:28:44.800469Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_solution1.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:44.802773Z","iopub.execute_input":"2025-09-24T08:28:44.803477Z","iopub.status.idle":"2025-09-24T08:28:44.876648Z","shell.execute_reply.started":"2025-09-24T08:28:44.803441Z","shell.execute_reply":"2025-09-24T08:28:44.875642Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution2","metadata":{}},{"cell_type":"code","source":"from astropy.stats import sigma_clip\nfrom scipy.signal import savgol_filter\nROOT_PATH = \"/kaggle/input/ariel-data-challenge-2025\"\nMODE = \"test\"\nimport time\n__t0 = time.perf_counter()\n\nclass Config:\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n\n    SCALE = 0.9475\n    SIGMA = 0.00055\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA = 11 # 9\n    MODEL_POLYNOMIAL_DEGREE = 4\n    \n    N_JOBS = 3\n\ndef _phase_detector_signal(signal, cfg):\n    sl = cfg.MODEL_PHASE_DETECTION_SLICE\n    min_idx = int(np.argmin(signal[sl])) + sl.start\n    s1 = signal[:min_idx]; s2 = signal[min_idx:]\n    if s1.size < 3 or s2.size < 3:\n        return 0, len(signal) - 1\n    g1 = np.gradient(s1); g1_max = np.max(g1) if np.size(g1) else 0.0\n    g2 = np.gradient(s2); g2_max = np.max(g2) if np.size(g2) else 0.0\n    if g1_max != 0: g1 /= g1_max\n    if g2_max != 0: g2 /= g2_max\n    phase1 = int(np.argmin(g1)); phase2 = int(np.argmax(g2)) + min_idx\n    return phase1, phase2\n\ndef estimate_sigma_fgs(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_1 (для FGS1) длиной N_planets — мягкий множитель к cfg.SIGMA.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n    for single in preprocessed_data:\n        # фазы по AIRS белой кривой — так же, как в модели\n        air_white = savgol_filter(single[:, 1:].mean(axis=1), 20, 2)\n        p1, p2 = _phase_detector_signal(air_white, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(air_white) - delta - 1, p2)\n\n        fgs = single[:, 0]\n        oot = (fgs[: p1 - delta] if p1 - delta > 0 else np.empty(0, fgs.dtype))\n        if p2 + delta < fgs.size:\n            oot = np.concatenate([oot, fgs[p2 + delta :]])\n        inn = fgs[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(fgs))\n        # относительная неопределённость глубины (в тех же ед., что s)\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель: корень, и узкий клип, чтобы не рисковать\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.8, 1.25)  # ±20–25% от базовой σ\n\n    return k * cfg.SIGMA * 1.042\n\ndef estimate_sigma_air(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_air длиной N_planets — мягкий множитель к cfg.SIGMA для всех AIRS-каналов.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n\n    for single in preprocessed_data:\n        # белая кривая AIRS на бинированных данных (после всех твоих весов по λ)\n        white = np.nanmean(single[:, 1:], axis=1)         # (n_bins,)\n        white_s = savgol_filter(white, 20, 2)             # для фаз\n\n        p1, p2 = _phase_detector_signal(white_s, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(white) - delta - 1, p2)\n\n        oot_left = white[: p1 - delta] if p1 - delta > 0 else np.empty(0, white.dtype)\n        oot_right = white[p2 + delta :] if (p2 + delta) < white.size else np.empty(0, white.dtype)\n        oot = np.concatenate([oot_left, oot_right]) if (oot_left.size + oot_right.size) else oot_left\n        inn = white[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(white))\n\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель вокруг медианы\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.90, 1.20)  # ±10%–20%\n\n    return k * cfg.SIGMA * 0.78\n\n\nclass SignalProcessor:\n    def __init__(self, config):\n        self.cfg = config\n        self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\n        self.planet_ids = pd.read_csv(f'{self.cfg.DATA_PATH}/{self.cfg.DATASET}_star_info.csv', index_col='planet_id').index.astype(int)\n\n    def _apply_linear_corr(self, linear_corr, signal):\n\n        coeffs = np.flip(linear_corr, axis=0)      # shape: (D, X, Y), D — старшая степень сначала\n        x = signal.astype(np.float64, copy=False)  # считаем в float64 для стабильности\n        out = np.empty_like(x, dtype=np.float64)\n        out[...] = coeffs[0]  # broadcast (X,Y) -> (T,X,Y)\n        for k in range(1, coeffs.shape[0]):\n            np.multiply(out, x, out=out)  # in-place умножение\n            out += coeffs[k]              # broadcast (X,Y)\n\n        return out.astype(signal.dtype, copy=False)\n\n    def _calibrate_single_signal(self, planet_id, sensor):\n        \"\"\"\n        Калибровка single-node сигнала.\n        Политика масок: DEAD — маскируем, HOT — НЕ маскируем (оставляем в данных).\n        \"\"\"\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n    \n        # --- load ---\n        signal = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_signal_0.parquet\"\n        ).to_numpy()\n        dark = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dark.parquet\"\n        ).to_numpy()\n        dead = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dead.parquet\"\n        ).to_numpy()\n        flat = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/flat.parquet\"\n        ).to_numpy()\n        linear_corr = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\"\n        ).values.astype(np.float64).reshape(sensor_cfg[\"linear_corr_shape\"])\n    \n        # --- reshape & ADC ---\n        signal = signal.reshape(sensor_cfg[\"raw_shape\"])\n        gain = self.adc_info[f\"{sensor}_adc_gain\"].iloc[0]\n        offset = self.adc_info[f\"{sensor}_adc_offset\"].iloc[0]\n        signal = signal / gain + offset  # сохраняем твою формулу\n    \n        # HOT только для мониторинга, не для маскирования\n        hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    \n        # --- crop per sensor ---\n        if sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            linear_corr = linear_corr[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dark = dark[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dead = dead[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            flat = flat[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            hot = hot[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]  # только для логов\n    \n        if sensor == \"FGS1\":\n            y0, y1, x0, x1 = 10, 22, 10, 22\n            signal = signal[:, y0:y1, x0:x1]\n            dark   = dark[y0:y1, x0:x1]\n            dead   = dead[y0:y1, x0:x1]\n            flat   = flat[y0:y1, x0:x1]\n            linear_corr = linear_corr[:, y0:y1, x0:x1]\n            hot    = hot[y0:y1, x0:x1]  # только для логов\n    \n        # --- non-neg clamp before linearity corr (как у тебя) ---\n        np.maximum(signal, 0, out=signal)\n    \n        # --- linearity correction ---\n        if sensor == \"FGS1\":\n            signal = self._apply_linear_corr(linear_corr, signal)\n        elif sensor == \"AIRS-CH0\":\n            sl = (slice(None), slice(10, 22), slice(None))  # T, Y, λ\n            signal[sl] = self._apply_linear_corr(linear_corr[:, 10:22, :], signal[sl])\n        else:\n            signal = self._apply_linear_corr(linear_corr, signal)\n    \n        # --- dark subtraction с учётом паттерна интеграций ---\n        base_dt, increment = sensor_cfg[\"dt_pattern\"]\n        even_scale = base_dt\n        odd_scale  = base_dt + increment\n        signal[::2]  -= dark * even_scale\n        signal[1::2] -= dark * odd_scale\n    \n        # --- APPLY FLAT (HOT-KEEP: не включаем hot в маску!) ---\n        if sensor == \"FGS1\":\n            flat_roi = flat.astype(signal.dtype, copy=False).copy()      # (12,12)\n            bad = (dead) | ~np.isfinite(flat_roi) | (flat_roi == 0)      # ← ТОЛЬКО dead/invalid\n            flat_roi[bad] = np.nan\n            signal /= flat_roi\n    \n        elif sensor == \"AIRS-CH0\":\n            y0, y1 = 10, 22\n            flat_roi = flat[y0:y1, :].astype(signal.dtype, copy=False).copy()  # (12, λ)\n            bad = (dead[y0:y1, :]) | ~np.isfinite(flat_roi) | (flat_roi == 0)  # ← ТОЛЬКО dead/invalid\n            flat_roi[bad] = np.nan\n            signal[:, y0:y1, :] /= flat_roi\n    \n        else:\n            flat2 = flat.astype(signal.dtype, copy=False).copy()\n            bad2 = (dead) | ~np.isfinite(flat2) | (flat2 == 0)                  # ← ТОЛЬКО dead/invalid\n            flat2[bad2] = np.nan\n            signal /= flat2\n        # --- END FLAT ---\n    \n        # (опционально) логируем метрики hot/dead\n        if getattr(self.cfg, \"LOG_HOT_STATS\", False):\n            if not hasattr(self, \"stats\"):\n                self.stats = []\n            self.stats.append({\n                \"planet_id\": int(planet_id),\n                \"sensor\": sensor,\n                \"hot_frac\": float(np.mean(hot)),\n                \"dead_frac\": float(np.mean(dead)),\n            })\n    \n        return signal\n\n\n    def _preprocess_calibrated_signal(self, calibrated_signal, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n        binning = sensor_cfg[\"binning\"]\n\n        if sensor == \"AIRS-CH0\":\n            signal_roi = calibrated_signal[:, 10:22, :]\n        elif sensor == \"FGS1\":\n            signal_roi = calibrated_signal[:, 10:22, 10:22]\n            signal_roi = signal_roi.reshape(signal_roi.shape[0], -1)\n        \n        mean_signal = np.nanmean(signal_roi, axis=1)\n\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        n_bins = cds_signal.shape[0] // binning\n        binned = np.array([\n            cds_signal[j*binning : (j+1)*binning].mean(axis=0) \n            for j in range(n_bins)\n        ])\n\n        if sensor == \"AIRS-CH0\":\n            q_lo = np.nanpercentile(binned, 5.0, axis=1, keepdims=True)    # (n_bins, 1)\n            q_hi = np.nanpercentile(binned, 95.0, axis=1, keepdims=True)   # (n_bins, 1)\n            np.clip(binned, q_lo, q_hi, out=binned)\n\n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        if sensor == \"AIRS-CH0\":\n            var = np.nanvar(binned, axis=0, ddof=1)                 # (λ, )\n            med = np.nanmedian(var)\n            safe_var = np.where(~np.isfinite(var) | (var <= 0), med if (np.isfinite(med) and med > 0) else 1.0, var)\n            w = 1.0 / safe_var\n\n            lo, hi = np.nanpercentile(w, 5.0), np.nanpercentile(w, 95.0)\n            if np.isfinite(lo) and np.isfinite(hi) and lo < hi:\n                w = np.clip(w, lo, hi)\n\n            M = binned.shape[1]\n            s = np.nansum(w)\n            if np.isfinite(s) and s > 0:\n                w = w * (M / s)\n            else:\n                w = np.ones_like(w)\n\n            binned *= w[None, :]\n\n\n        return binned\n\n    def _process_planet_sensor(self, args):\n        planet_id, sensor = args['planet_id'], args['sensor']\n        calibrated = self._calibrate_single_signal(planet_id, sensor)\n        preprocessed = self._preprocess_calibrated_signal(calibrated, sensor)\n        return preprocessed\n\n    def process_all_data(self):\n        args_fgs1 = [dict(planet_id=planet_id, sensor=\"FGS1\") for planet_id in self.planet_ids]\n        preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        args_airs_ch0 = [dict(planet_id=planet_id, sensor=\"AIRS-CH0\") for planet_id in self.planet_ids]\n        preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        preprocessed_signal = np.concatenate(\n            [np.stack(preprocessed_fgs1), np.stack(preprocessed_airs_ch0)], axis=2\n        )\n        return preprocessed_signal\n    \n\nclass TransitModel:\n    def __init__(self, config):\n        self.cfg = config\n\n    def _phase_detector(self, signal):\n        search_slice = self.cfg.MODEL_PHASE_DETECTION_SLICE\n        min_index = np.argmin(signal[search_slice]) + search_slice.start\n        \n        signal1 = signal[:min_index]\n        signal2 = signal[min_index:]\n\n        grad1 = np.gradient(signal1)\n        grad1 /= grad1.max()\n        \n        grad2 = np.gradient(signal2)\n        grad2 /= grad2.max()\n\n        phase1 = np.argmin(grad1)\n        phase2 = np.argmax(grad2) + min_index\n\n        return phase1, phase2\n    \n    def _objective_function(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\n        coeffs = np.polyfit(x, y, deg=power)\n        poly = np.poly1d(coeffs)\n        error = np.abs(poly(x) - y).mean()\n        \n        return error\n\n    def predict(self, single_preprocessed_signal):\n        signal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)\n        signal_1d = savgol_filter(signal_1d, 23, 2)\n        \n        phase1, phase2 = self._phase_detector(signal_1d)\n\n        phase1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, phase1)\n        phase2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, phase2)    \n\n        result = minimize(\n            fun=self._objective_function,\n            x0=[0.0001],\n            args=(signal_1d, phase1, phase2),\n            method=\"Nelder-Mead\"\n        )\n        \n        return result.x[0]\n\n    def predict_all(self, preprocessed_signals):\n        predictions = [\n            self.predict(preprocessed_signal)\n            for preprocessed_signal in tqdm(preprocessed_signals)\n        ]\n        return np.array(predictions) * self.cfg.SCALE\nclass ResidualBlock(nn.Module):\n    def __init__(self, dim, p=0.2):\n        super().__init__()\n        self.fc1 = nn.Linear(dim, dim)\n        self.bn1 = nn.BatchNorm1d(dim)\n        self.fc2 = nn.Linear(dim, dim)\n        self.bn2 = nn.BatchNorm1d(dim)\n        self.relu = nn.ReLU()\n        self.dropout = nn.Dropout(p)\n\n    def forward(self, x):\n        identity = x\n        out = self.relu(self.bn1(self.fc1(x)))\n        out = self.dropout(out)\n        out = self.bn2(self.fc2(out))\n        return self.relu(out + identity)\n\n\nclass ResNetMLP(nn.Module):\n    def __init__(self, input_dim=3, hidden_dim=32, output_dim=1, num_blocks=3, dropout_rate=0.2):\n        super().__init__()\n        self.input_layer = nn.Linear(input_dim, hidden_dim)\n        self.blocks = nn.Sequential(*[ResidualBlock(hidden_dim, p=dropout_rate) for _ in range(num_blocks)])\n        self.output_layer = nn.Linear(hidden_dim, output_dim)\n\n    def forward(self, x):\n        x = self.input_layer(x)\n        x = self.blocks(x)\n        x = self.output_layer(x)\n        return x\n        \nresnet = ResNetMLP(num_blocks=80, dropout_rate=0.3)\nresnet.load_state_dict(torch.load(\"/kaggle/input/fgs1/pytorch/default/1/best_model.pth\"))\nresnet.eval()\nStarInfo = pd.read_csv(ROOT_PATH + f\"/{MODE}_star_info.csv\")\nStarInfo[\"planet_id\"] = StarInfo[\"planet_id\"].astype(int)\nPlanetIds = StarInfo[\"planet_id\"].tolist()\nStarInfo = StarInfo.set_index(\"planet_id\")\n\nclass SubmissionGenerator:\n    def __init__(self, config):\n        self.cfg = config\n        self.sample_submission = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", index_col=\"planet_id\")\n\n    def create(self, predictions1, predictions2, predictions, sigma_fgs=None, sigma_air=None):\n        planet_ids = self.sample_submission.index\n        n_mu = self.sample_submission.shape[1] // 2  # 283\n\n        preds = np.asarray(predictions, dtype=float).reshape(-1)\n        mu = np.tile(preds.reshape(-1, 1), (1, n_mu))\n        mu = np.clip(mu, 0, None)\n\n        sigmas = np.full_like(mu, self.cfg.SIGMA, dtype=float)\n        if sigma_fgs is not None:\n            sigma_fgs = np.asarray(sigma_fgs, dtype=float).reshape(-1)\n            sigmas[:, 0] = np.clip(sigma_fgs, 1e-6, 0.1)\n        if sigma_air is not None:\n            sigma_air = np.asarray(sigma_air, dtype=float).reshape(-1, 1)\n            sigmas[:, 1:] = np.clip(sigma_air, 1e-6, 0.1)\n\n        submission_df = pd.DataFrame(\n            np.concatenate([mu, sigmas], axis=1),\n            columns=self.sample_submission.columns,\n            index=planet_ids\n        )\n        submission_df.iloc[:, 0] = predictions1\n        submission_df.iloc[:, 1:283] = predictions2\n        #submission_df.to_csv(\"submission.csv\")\n        \n        return submission_df\n\n\n\nconfig = Config()\n    \nsignal_processor = SignalProcessor(config)\npreprocessed_data = signal_processor.process_all_data()\n\nmodel = TransitModel(config)\npredictions = model.predict_all(preprocessed_data)\nsigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)\nsigma_air_vec = estimate_sigma_air(preprocessed_data, config)\npredictions\npredictions_df = pd.DataFrame({\n    \"planet_id\": PlanetIds,\n    \"transit_depth\": predictions\n})\n\ninput_df = pd.merge(predictions_df, StarInfo, on=\"planet_id\", how=\"left\")\ninput_df[\"transit_depth\"] *= 10000\nfeatures = ['transit_depth','Rs','i']\nX = input_df[features].values.astype(np.float32)\nX_tensor = torch.tensor(X, dtype=torch.float32)\nwith torch.no_grad():\n    predictions1 = resnet(X_tensor).numpy()\npredictions1 /= 10000\n\nclass ResidualBlock2(nn.Module):\n    def __init__(self, dim, p=0.2):\n        super().__init__()\n        self.fc1 = nn.Linear(dim, dim)\n        self.bn1 = nn.BatchNorm1d(dim)\n        self.fc2 = nn.Linear(dim, dim)\n        self.bn2 = nn.BatchNorm1d(dim)\n        self.relu = nn.ReLU()\n        self.dropout = nn.Dropout(p)\n\n    def forward(self, x):\n        identity = x\n        out = self.relu(self.bn1(self.fc1(x)))\n        out = self.dropout(out)\n        out = self.bn2(self.fc2(out))\n        return self.relu(out + identity)\n\n\nclass ResNetMLP2(nn.Module):\n    def __init__(self, input_dim=3, hidden_dim=128, output_dim=282, num_blocks=3, dropout_rate=0.2):\n        super().__init__()\n        self.input_layer = nn.Linear(input_dim, hidden_dim)\n        self.blocks = nn.Sequential(*[ResidualBlock2(hidden_dim, p=dropout_rate) for _ in range(num_blocks)])\n        self.output_layer = nn.Linear(hidden_dim, output_dim)\n\n    def forward(self, x):\n        x = self.input_layer(x)\n        x = self.blocks(x)\n        x = self.output_layer(x)\n        return x\nresnet2 = ResNetMLP2(num_blocks=120, dropout_rate=0.3)\nresnet2.load_state_dict(torch.load(\"/kaggle/input/airs-deeper/pytorch/default/1/best_model deep1.pth\"))\nresnet2.eval()\n\nwith torch.no_grad():\n    predictions2 = resnet2(X_tensor).numpy()\npredictions2 /= 10000\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:44.878061Z","iopub.execute_input":"2025-09-24T08:28:44.878559Z","iopub.status.idle":"2025-09-24T08:28:50.9789Z","shell.execute_reply.started":"2025-09-24T08:28:44.878522Z","shell.execute_reply":"2025-09-24T08:28:50.978007Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_generator = SubmissionGenerator(config)\nsubmission_solution2 = submission_generator.create(predictions1, predictions2, predictions, sigma_fgs=sigma_fgs_vec, sigma_air=sigma_air_vec)\n\nsubmission_solution2.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:50.980045Z","iopub.execute_input":"2025-09-24T08:28:50.980533Z","iopub.status.idle":"2025-09-24T08:28:51.01953Z","shell.execute_reply.started":"2025-09-24T08:28:50.980497Z","shell.execute_reply":"2025-09-24T08:28:51.01828Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution 3","metadata":{}},{"cell_type":"code","source":"ROOT_PATH = \"/kaggle/input/ariel-data-challenge-2025\"\nMODE = \"test\"\n\n__t0 = time.perf_counter()\n\nclass Config:\n    FEATURES = ['transit_depth', 'Rs', 'i']\n    # FEATURES = ['transit_depth', 'Rs', 'Ms', 'Ts', 'Mp', 'e', 'P', 'sma', 'i']\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n    load_data = True  # Set to True to load from disk, False to generate\n    DEBUG = False\n\n    SCALE = 0.96\n    SIGMA = 0.00056\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA = 11 # 9\n    MODEL_POLYNOMIAL_DEGREE = 3\n    \n    N_JOBS = 3\n\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    row_id_column_name: str,\n    naive_mean: float,\n    naive_sigma: float,\n    fsg_sigma_true: float = 1e-6,\n    airs_sigma_true: float = 1e-5,\n    fgs_weight: float = 1,\n) -> float:\n    \"\"\"\n    This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n    we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n    x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n    the GLL value for one spectrum is the sum of all of them.\n\n    Inputs:\n        - solution: Ground Truth spectra (from test set)\n            - shape: (nsamples, n_wavelengths)\n        - submission: Predicted spectra and errors (from participants)\n            - shape: (nsamples, n_wavelengths*2)\n        naive_mean: (float) mean from the train set.\n        naive_sigma: (float) standard deviation from the train set.\n        fsg_sigma_true: (float) standard deviation from the FSG1 instrument for the test set.\n        airs_sigma_true: (float) standard deviation from the AIRS instrument for the test set.\n        fgs_weight: (float) relative weight of the fgs channel\n    \"\"\"\n\n    del solution[row_id_column_name]\n    del submission[row_id_column_name]\n\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\n    for col in submission.columns:\n        if not pandas.api.types.is_numeric_dtype(submission[col]):\n            raise ParticipantVisibleError(f'Submission column {col} must be a number')\n\n    n_wavelengths = len(solution.columns)\n    if len(submission.columns) != n_wavelengths * 2:\n        raise ParticipantVisibleError('Wrong number of columns in the submission')\n\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=10**-15, a_max=None)\n    sigma_true = np.append(\n        np.array(\n            [\n                fsg_sigma_true,\n            ]\n        ),\n        np.ones(n_wavelengths - 1) * airs_sigma_true,\n    )\n    y_true = solution.values\n\n    GLL_pred = scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred)\n    GLL_true = scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true))\n    GLL_mean = scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true))\n\n    # normalise the score, right now it becomes a matrix instead of a scalar.\n    ind_scores = (GLL_pred - GLL_mean) / (GLL_true - GLL_mean)\n\n    weights = np.append(np.array([fgs_weight]), np.ones(len(solution.columns) - 1))\n    weights = weights * np.ones_like(ind_scores)\n    submit_score = np.average(ind_scores, weights=weights)\n    return float(np.clip(submit_score, 0.0, 1.0))\n\ndef _phase_detector_signal(signal, cfg):\n    sl = cfg.MODEL_PHASE_DETECTION_SLICE\n    min_idx = int(np.argmin(signal[sl])) + sl.start\n    s1 = signal[:min_idx]; s2 = signal[min_idx:]\n    if s1.size < 3 or s2.size < 3:\n        return 0, len(signal) - 1\n    g1 = np.gradient(s1); g1_max = np.max(g1) if np.size(g1) else 0.0\n    g2 = np.gradient(s2); g2_max = np.max(g2) if np.size(g2) else 0.0\n    if g1_max != 0: g1 /= g1_max\n    if g2_max != 0: g2 /= g2_max\n    phase1 = int(np.argmin(g1)); phase2 = int(np.argmax(g2)) + min_idx\n    return phase1, phase2\n\ndef estimate_sigma_fgs(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_1 (для FGS1) длиной N_planets — мягкий множитель к cfg.SIGMA.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n    for single in preprocessed_data:\n        # фазы по AIRS белой кривой — так же, как в модели\n        air_white = savgol_filter(single[:, 1:].mean(axis=1), 20, 2)\n        p1, p2 = _phase_detector_signal(air_white, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(air_white) - delta - 1, p2)\n\n        fgs = single[:, 0]\n        oot = (fgs[: p1 - delta] if p1 - delta > 0 else np.empty(0, fgs.dtype))\n        if p2 + delta < fgs.size:\n            oot = np.concatenate([oot, fgs[p2 + delta :]])\n        inn = fgs[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(fgs))\n        # относительная неопределённость глубины (в тех же ед., что s)\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель: корень, и узкий клип, чтобы не рисковать\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.8, 1.25)  # ±20–25% от базовой σ\n\n    return k * cfg.SIGMA\n\ndef estimate_sigma_air(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_air длиной N_planets — мягкий множитель к cfg.SIGMA для всех AIRS-каналов.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n\n    for single in preprocessed_data:\n        # белая кривая AIRS на бинированных данных (после всех твоих весов по λ)\n        white = np.nanmean(single[:, 1:], axis=1)         # (n_bins,)\n        white_s = savgol_filter(white, 20, 2)             # для фаз\n\n        p1, p2 = _phase_detector_signal(white_s, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(white) - delta - 1, p2)\n\n        oot_left = white[: p1 - delta] if p1 - delta > 0 else np.empty(0, white.dtype)\n        oot_right = white[p2 + delta :] if (p2 + delta) < white.size else np.empty(0, white.dtype)\n        oot = np.concatenate([oot_left, oot_right]) if (oot_left.size + oot_right.size) else oot_left\n        inn = white[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(white))\n\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель вокруг медианы\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.90, 1.20)  # ±10%–20%\n\n    return k * cfg.SIGMA\n\nclass SignalProcessor:\n    def __init__(self, config):\n        self.cfg = config\n        self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\n        self.planet_ids = pd.read_csv(f'{self.cfg.DATA_PATH}/{self.cfg.DATASET}_star_info.csv', index_col='planet_id').index.astype(int)\n\n    def _apply_linear_corr(self, linear_corr, signal):\n\n        coeffs = np.flip(linear_corr, axis=0)      # shape: (D, X, Y), D — старшая степень сначала\n        x = signal.astype(np.float64, copy=False)  # считаем в float64 для стабильности\n        out = np.empty_like(x, dtype=np.float64)\n        out[...] = coeffs[0]  # broadcast (X,Y) -> (T,X,Y)\n        for k in range(1, coeffs.shape[0]):\n            np.multiply(out, x, out=out)  # in-place умножение\n            out += coeffs[k]              # broadcast (X,Y)\n\n        return out.astype(signal.dtype, copy=False)\n\n    def _calibrate_single_signal(self, planet_id, sensor):\n        \"\"\"\n        Калибровка single-node сигнала.\n        Политика масок: DEAD — маскируем, HOT — НЕ маскируем (оставляем в данных).\n        \"\"\"\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n    \n        # --- load ---\n        signal = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_signal_0.parquet\"\n        ).to_numpy()\n        dark = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dark.parquet\"\n        ).to_numpy()\n        dead = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dead.parquet\"\n        ).to_numpy()\n        flat = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/flat.parquet\"\n        ).to_numpy()\n        linear_corr = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\"\n        ).values.astype(np.float64).reshape(sensor_cfg[\"linear_corr_shape\"])\n    \n        # --- reshape & ADC ---\n        signal = signal.reshape(sensor_cfg[\"raw_shape\"])\n        gain = self.adc_info[f\"{sensor}_adc_gain\"].iloc[0]\n        offset = self.adc_info[f\"{sensor}_adc_offset\"].iloc[0]\n        signal = signal / gain + offset  # сохраняем твою формулу\n    \n        # HOT только для мониторинга, не для маскирования\n        hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    \n        # --- crop per sensor ---\n        if sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            linear_corr = linear_corr[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dark = dark[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dead = dead[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            flat = flat[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            hot = hot[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]  # только для логов\n    \n        if sensor == \"FGS1\":\n            y0, y1, x0, x1 = 10, 22, 10, 22\n            signal = signal[:, y0:y1, x0:x1]\n            dark   = dark[y0:y1, x0:x1]\n            dead   = dead[y0:y1, x0:x1]\n            flat   = flat[y0:y1, x0:x1]\n            linear_corr = linear_corr[:, y0:y1, x0:x1]\n            hot    = hot[y0:y1, x0:x1]  # только для логов\n    \n        # --- non-neg clamp before linearity corr (как у тебя) ---\n        np.maximum(signal, 0, out=signal)\n    \n        # --- linearity correction ---\n        if sensor == \"FGS1\":\n            signal = self._apply_linear_corr(linear_corr, signal)\n        elif sensor == \"AIRS-CH0\":\n            sl = (slice(None), slice(10, 22), slice(None))  # T, Y, λ\n            signal[sl] = self._apply_linear_corr(linear_corr[:, 10:22, :], signal[sl])\n        else:\n            signal = self._apply_linear_corr(linear_corr, signal)\n    \n        # --- dark subtraction с учётом паттерна интеграций ---\n        base_dt, increment = sensor_cfg[\"dt_pattern\"]\n        even_scale = base_dt\n        odd_scale  = base_dt + increment\n        signal[::2]  -= dark * even_scale\n        signal[1::2] -= dark * odd_scale\n    \n        # --- APPLY FLAT (HOT-KEEP: не включаем hot в маску!) ---\n        if sensor == \"FGS1\":\n            flat_roi = flat.astype(signal.dtype, copy=False).copy()      # (12,12)\n            bad = (dead) | ~np.isfinite(flat_roi) | (flat_roi == 0)      # ← ТОЛЬКО dead/invalid\n            flat_roi[bad] = np.nan\n            signal /= flat_roi\n    \n        elif sensor == \"AIRS-CH0\":\n            y0, y1 = 10, 22\n            flat_roi = flat[y0:y1, :].astype(signal.dtype, copy=False).copy()  # (12, λ)\n            bad = (dead[y0:y1, :]) | ~np.isfinite(flat_roi) | (flat_roi == 0)  # ← ТОЛЬКО dead/invalid\n            flat_roi[bad] = np.nan\n            signal[:, y0:y1, :] /= flat_roi\n    \n        else:\n            flat2 = flat.astype(signal.dtype, copy=False).copy()\n            bad2 = (dead) | ~np.isfinite(flat2) | (flat2 == 0)                  # ← ТОЛЬКО dead/invalid\n            flat2[bad2] = np.nan\n            signal /= flat2\n        # --- END FLAT ---\n    \n        # (опционально) логируем метрики hot/dead\n        if getattr(self.cfg, \"LOG_HOT_STATS\", False):\n            if not hasattr(self, \"stats\"):\n                self.stats = []\n            self.stats.append({\n                \"planet_id\": int(planet_id),\n                \"sensor\": sensor,\n                \"hot_frac\": float(np.mean(hot)),\n                \"dead_frac\": float(np.mean(dead)),\n            })\n    \n        return signal\n\n    def _preprocess_calibrated_signal(self, calibrated_signal, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n        binning = sensor_cfg[\"binning\"]\n\n        if sensor == \"AIRS-CH0\":\n            signal_roi = calibrated_signal[:, 10:22, :]\n        elif sensor == \"FGS1\":\n            signal_roi = calibrated_signal[:, 10:22, 10:22]\n            signal_roi = signal_roi.reshape(signal_roi.shape[0], -1)\n        \n        mean_signal = np.nanmean(signal_roi, axis=1)\n\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        n_bins = cds_signal.shape[0] // binning\n        binned = np.array([\n            cds_signal[j*binning : (j+1)*binning].mean(axis=0) \n            for j in range(n_bins)\n        ])\n\n        if sensor == \"AIRS-CH0\":\n            q_lo = np.nanpercentile(binned, 5.0, axis=1, keepdims=True)    # (n_bins, 1)\n            q_hi = np.nanpercentile(binned, 95.0, axis=1, keepdims=True)   # (n_bins, 1)\n            np.clip(binned, q_lo, q_hi, out=binned)\n\n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        if sensor == \"AIRS-CH0\":\n            var = np.nanvar(binned, axis=0, ddof=1)                 # (λ, )\n            med = np.nanmedian(var)\n            safe_var = np.where(~np.isfinite(var) | (var <= 0), med if (np.isfinite(med) and med > 0) else 1.0, var)\n            w = 1.0 / safe_var\n\n            lo, hi = np.nanpercentile(w, 5.0), np.nanpercentile(w, 95.0)\n            if np.isfinite(lo) and np.isfinite(hi) and lo < hi:\n                w = np.clip(w, lo, hi)\n\n            M = binned.shape[1]\n            s = np.nansum(w)\n            if np.isfinite(s) and s > 0:\n                w = w * (M / s)\n            else:\n                w = np.ones_like(w)\n\n            binned *= w[None, :]\n\n\n        return binned\n\n    def _process_planet_sensor(self, args):\n        planet_id, sensor = args['planet_id'], args['sensor']\n        calibrated = self._calibrate_single_signal(planet_id, sensor)\n        preprocessed = self._preprocess_calibrated_signal(calibrated, sensor)\n        return preprocessed\n\n    def process_all_data(self):\n        args_fgs1 = [dict(planet_id=planet_id, sensor=\"FGS1\") for planet_id in self.planet_ids]\n        preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        args_airs_ch0 = [dict(planet_id=planet_id, sensor=\"AIRS-CH0\") for planet_id in self.planet_ids]\n        preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        preprocessed_signal = np.concatenate(\n            [np.stack(preprocessed_fgs1), np.stack(preprocessed_airs_ch0)], axis=2\n        )\n        return preprocessed_signal\n    \n\nclass TransitModel:\n    def __init__(self, config):\n        self.cfg = config\n\n    def _phase_detector(self, signal):\n        search_slice = self.cfg.MODEL_PHASE_DETECTION_SLICE\n        min_index = np.argmin(signal[search_slice]) + search_slice.start\n        \n        signal1 = signal[:min_index]\n        signal2 = signal[min_index:]\n\n        grad1 = np.gradient(signal1)\n        grad1 /= grad1.max()\n        \n        grad2 = np.gradient(signal2)\n        grad2 /= grad2.max()\n\n        phase1 = np.argmin(grad1)\n        phase2 = np.argmax(grad2) + min_index\n\n        return phase1, phase2\n    \n    def _objective_function(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\n        coeffs = np.polyfit(x, y, deg=power)\n        poly = np.poly1d(coeffs)\n        error = np.abs(poly(x) - y).mean()\n        \n        return error\n\n    def predict(self, single_preprocessed_signal):\n        signal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)\n        signal_1d = savgol_filter(signal_1d, 23, 2)\n        \n        phase1, phase2 = self._phase_detector(signal_1d)\n\n        phase1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, phase1)\n        phase2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, phase2)    \n\n        result = minimize(\n            fun=self._objective_function,\n            x0=[0.0001],\n            args=(signal_1d, phase1, phase2),\n            method=\"Nelder-Mead\"\n        )\n        \n        return result.x[0]\n\n    def predict_all(self, preprocessed_signals):\n        predictions = [\n            self.predict(preprocessed_signal)\n            for preprocessed_signal in tqdm(preprocessed_signals)\n        ]\n        return np.array(predictions) * self.cfg.SCALE\n\nStarInfo = pd.read_csv(ROOT_PATH + f\"/{MODE}_star_info.csv\")\nStarInfo[\"planet_id\"] = StarInfo[\"planet_id\"].astype(int)\nPlanetIds = StarInfo[\"planet_id\"].tolist()\nStarInfo = StarInfo.set_index(\"planet_id\")\nclass SubmissionGenerator:\n    def __init__(self, config):\n        self.cfg = config\n        self.sample_submission = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", index_col=\"planet_id\")\n\n    def create(self, predictions1, predictions, sigma_fgs=None, sigma_air=None):\n        planet_ids = self.sample_submission.index\n        n_mu = self.sample_submission.shape[1] // 2  # 283\n\n        preds = np.asarray(predictions, dtype=float).reshape(-1)\n        mu = np.tile(preds.reshape(-1, 1), (1, n_mu))\n        mu = np.clip(mu, 0, None)\n\n        sigmas = np.full_like(mu, self.cfg.SIGMA, dtype=float)\n        if sigma_fgs is not None:\n            sigma_fgs = np.asarray(sigma_fgs, dtype=float).reshape(-1)\n            sigmas[:, 0] = np.clip(sigma_fgs, 1e-6, 0.1)\n        if sigma_air is not None:\n            sigma_air = np.asarray(sigma_air, dtype=float).reshape(-1, 1)\n            sigmas[:, 1:] = np.clip(sigma_air, 1e-6, 0.1)\n\n        submission_df = pd.DataFrame(\n            np.concatenate([mu, sigmas], axis=1),\n            columns=self.sample_submission.columns,\n            index=planet_ids\n        )\n        submission_df.iloc[:, 0] = predictions\n        submission_df.iloc[:, 1:283] = predictions1\n        submission_df.to_csv(\"submission.csv\")\n        \n        return submission_df\n\nclass ResidualBlock2(nn.Module):\n    def __init__(self, dim, p=0.2):\n        super().__init__()\n        self.fc1 = nn.Linear(dim, dim)\n        self.fc2 = nn.Linear(dim, dim)\n        self.relu = nn.ReLU()\n        self.dropout = nn.Dropout(p)\n\n    def forward(self, x):\n        identity = x\n        out = self.relu(self.fc1(x))\n        out = self.dropout(out)\n        out = self.fc2(out)\n        return self.relu(out + identity)\n\nclass ResNetMLP2(nn.Module):\n    def __init__(self, input_dim=3, hidden_dim=128, output_dim=282, num_blocks=3, dropout_rate=0.2):\n        super().__init__()\n        self.input_layer = nn.Linear(input_dim, hidden_dim)\n        self.blocks = nn.Sequential(*[ResidualBlock2(hidden_dim, p=dropout_rate) for _ in range(num_blocks)])\n        self.output_layer = nn.Linear(hidden_dim, output_dim)\n\n    def forward(self, x):\n        x = self.input_layer(x)\n        x = self.blocks(x)\n        x = self.output_layer(x)\n        return x\n\n# Execute the training\n# Train and get scalers\n# all_models, scaler_X, scaler_y = train_new_model()\n\ndef load_cv_models_and_scalers(directory):\n    \"\"\"\n    Loads all cross-validation models and scalers from the specified directory.\n    Args:\n        directory (str): Path to the directory containing model and scaler files.\n    Returns:\n        all_models (list): List of loaded models.\n        scaler_X: Loaded X scaler.\n        scaler_y: Loaded y scaler.\n    \"\"\"\n    import os\n    import joblib\n    # Load scalers\n    scaler_X = joblib.load(os.path.join(directory, 'scaler_X.joblib'))\n    scaler_y = joblib.load(os.path.join(directory, 'scaler_y.joblib'))\n\n    # Load all CV models\n    all_models = []\n    model_params = {\n        'input_dim': len(Config.FEATURES),  # Always use the current feature count\n        'hidden_dim': 256,\n        'output_dim': 282,\n        'num_blocks': 80,\n        'dropout_rate': 0.3\n    }\n    for fold in range(1, 11):\n        model = ResNetMLP2(**model_params).double()\n        model_path = os.path.join(directory, f'best_model_airs_cv_fold{fold}.pth')\n        model.load_state_dict(torch.load(model_path, map_location=torch.device('cpu')))\n        model.eval()\n        all_models.append(model)\n    return all_models, scaler_X, scaler_y\n    \nall_models, scaler_X, scaler_y = load_cv_models_and_scalers('/kaggle/input/ariel-2025-result/results_3_v25')\n\nconfig = Config()\nsignal_processor = SignalProcessor(config)\npreprocessed_data = signal_processor.process_all_data()\n\nmodel = TransitModel(config)\npredictions = model.predict_all(preprocessed_data)\nsigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)\nsigma_air_vec = estimate_sigma_air(preprocessed_data, config)\n\npredictions_df = pd.DataFrame({\n    \"planet_id\": PlanetIds,\n    \"transit_depth\": predictions\n})\n\ninput_df = pd.merge(predictions_df, StarInfo, on=\"planet_id\", how=\"left\")\nX = input_df[Config.FEATURES].values.astype(np.float64)\nX_scaled = scaler_X.transform(X)\nX_tensor = torch.tensor(X_scaled, dtype=torch.float64)\n\n# Generate average prediction from all CV models (on scaled X, then inverse transform)\nwith torch.no_grad():\n    preds_scaled = [model(X_tensor).numpy() for model in all_models]\npredictions1_scaled = np.mean(preds_scaled, axis=0)\npredictions1 = scaler_y.inverse_transform(predictions1_scaled)\n\nclass Config:\n    # FEATURES = ['transit_depth', 'Rs', 'i']\n    FEATURES = ['transit_depth', 'Rs', 'Ms', 'Ts', 'Mp', 'e', 'P', 'sma', 'i']\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n    load_data = True  # Set to True to load from disk, False to generate\n    DEBUG = False\n\n    SCALE = 0.96\n    SIGMA = 0.00055\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA = 11 # 9\n    MODEL_POLYNOMIAL_DEGREE = 3\n    \n    N_JOBS = 3\n\ndef load_cv_models_and_scalers(directory):\n    \"\"\"\n    Loads all cross-validation models and scalers from the specified directory.\n    Args:\n        directory (str): Path to the directory containing model and scaler files.\n    Returns:\n        all_models (list): List of loaded models.\n        scaler_X: Loaded X scaler.\n        scaler_y: Loaded y scaler.\n    \"\"\"\n    import os\n    import joblib\n    # Load scalers\n    scaler_X = joblib.load(os.path.join(directory, 'scaler_X.joblib'))\n    scaler_y = joblib.load(os.path.join(directory, 'scaler_y.joblib'))\n\n    # Load all CV models\n    all_models = []\n    model_params = {\n        'input_dim': len(Config.FEATURES),  # Always use the current feature count\n        'hidden_dim': 256,\n        'output_dim': 282,\n        'num_blocks': 80,\n        'dropout_rate': 0.3\n    }\n    for fold in range(1, 11):\n        model = ResNetMLP2(**model_params).double()\n        model_path = os.path.join(directory, f'best_model_airs_cv_fold{fold}.pth')\n        model.load_state_dict(torch.load(model_path, map_location=torch.device('cpu')))\n        model.eval()\n        all_models.append(model)\n    return all_models, scaler_X, scaler_y\n    \nall_models, scaler_X, scaler_y = load_cv_models_and_scalers('/kaggle/input/ariel-2025-result/results_9_v26')\nconfig = Config()\nsignal_processor = SignalProcessor(config)\npreprocessed_data = signal_processor.process_all_data()\n\nmodel = TransitModel(config)\npredictions = model.predict_all(preprocessed_data)\nsigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)\nsigma_air_vec = estimate_sigma_air(preprocessed_data, config)\n\npredictions_df = pd.DataFrame({\n    \"planet_id\": PlanetIds,\n    \"transit_depth\": predictions\n})\n\ninput_df = pd.merge(predictions_df, StarInfo, on=\"planet_id\", how=\"left\")\nX = input_df[Config.FEATURES].values.astype(np.float64)\nX_scaled = scaler_X.transform(X)\nX_tensor = torch.tensor(X_scaled, dtype=torch.float64)\n\n# Generate average prediction from all CV models (on scaled X, then inverse transform)\nwith torch.no_grad():\n    preds_scaled = [model(X_tensor).numpy() for model in all_models]\n\n# Correctly average the predictions\npredictions2_scaled = np.mean(preds_scaled, axis=0)\n\n# Correctly inverse transform the *averaged* predictions\npredictions2 = scaler_y.inverse_transform(predictions2_scaled)\n\n# Combine the predictions into a single array\nall_predictions = np.array([predictions1, predictions2])\n\n# Calculate the mean across the models (axis=0)\nfinal_predictions = np.mean(all_predictions, axis=0)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:28:51.021054Z","iopub.execute_input":"2025-09-24T08:28:51.021542Z","iopub.status.idle":"2025-09-24T08:29:39.526143Z","shell.execute_reply.started":"2025-09-24T08:28:51.021393Z","shell.execute_reply":"2025-09-24T08:29:39.524143Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_generator = SubmissionGenerator(config)\nsubmission_solution3 = submission_generator.create(final_predictions, predictions, sigma_fgs=sigma_fgs_vec, sigma_air=sigma_air_vec)\nsubmission_solution3.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:29:39.528335Z","iopub.execute_input":"2025-09-24T08:29:39.528734Z","iopub.status.idle":"2025-09-24T08:29:39.579596Z","shell.execute_reply.started":"2025-09-24T08:29:39.528694Z","shell.execute_reply":"2025-09-24T08:29:39.57825Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution4","metadata":{}},{"cell_type":"code","source":"# ======== パスとモード ========\nROOT_PATH = \"/kaggle/input/ariel-data-challenge-2025\"\nMODE = \"test\"\n__t0 = time.perf_counter()\n\n\n# ======== Config ========\nclass Config:\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n\n    SCALE = 0.928          # ※ #9 は今回は変更しない\n    SIGMA = 0.00056\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(27, 143)\n    MODEL_OPTIMIZATION_DELTA = 11\n    MODEL_POLYNOMIAL_DEGREE = 4\n    \n    N_JOBS = 3\n\n\n# ======== 位相補助 ========\ndef _phase_detector_signal(signal, cfg):\n    sl = cfg.MODEL_PHASE_DETECTION_SLICE\n    min_idx = int(np.argmin(signal[sl])) + sl.start\n    s1 = signal[:min_idx]; s2 = signal[min_idx:]\n    if s1.size < 3 or s2.size < 3:\n        return 0, len(signal) - 1\n    g1 = np.gradient(s1); g1_max = np.max(g1) if np.size(g1) else 0.0\n    g2 = np.gradient(s2); g2_max = np.max(g2) if np.size(g2) else 0.0\n    if g1_max != 0: g1 /= g1_max\n    if g2_max != 0: g2 /= g2_max\n    phase1 = int(np.argmin(g1)); phase2 = int(np.argmax(g2)) + min_idx\n    return phase1, phase2\n\n\n# ======== σ推定（#4: クリップ幅＆仕上げ係数を拡張） ========\ndef estimate_sigma_fgs(preprocessed_data, cfg):\n    \"\"\"FGS1用 σ のソフト推定（保守拡張版）\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n    for single in preprocessed_data:\n        air_white = savgol_filter(single[:, 1:].mean(axis=1), 20, 2)\n        p1, p2 = _phase_detector_signal(air_white, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(air_white) - delta - 1, p2)\n\n        fgs = single[:, 0]\n        oot = (fgs[: p1 - delta] if p1 - delta > 0 else np.empty(0, fgs.dtype))\n        if p2 + delta < fgs.size:\n            oot = np.concatenate([oot, fgs[p2 + delta :]])\n        inn = fgs[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(fgs))\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n\n    # --- #4: clipをやや緩める（0.8–1.25 → 0.85–1.30）---\n    k = np.clip(k, 0.85, 1.30)\n\n    sigma_fgs = k * cfg.SIGMA\n\n    # --- #4: 仕上げの全体係数（わずかに拡張）---\n    sigma_fgs *= 1.04\n    return sigma_fgs\n\n\ndef estimate_sigma_air(preprocessed_data, cfg):\n    \"\"\"AIRS用 σ のソフト推定（保守拡張版）\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n\n    for single in preprocessed_data:\n        white = np.nanmean(single[:, 1:], axis=1)\n        white_s = savgol_filter(white, 20, 2)\n\n        p1, p2 = _phase_detector_signal(white_s, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(white) - delta - 1, p2)\n\n        oot_left = white[: p1 - delta] if p1 - delta > 0 else np.empty(0, white.dtype)\n        oot_right = white[p2 + delta :] if (p2 + delta) < white.size else np.empty(0, white.dtype)\n        oot = np.concatenate([oot_left, oot_right]) if (oot_left.size + oot_right.size) else oot_left\n        inn = white[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(white))\n\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n\n    # --- #4: clipをやや緩める（0.90–1.20 → 0.92–1.22）---\n    k = np.clip(k, 0.92, 1.22)\n\n    sigma_air = k * cfg.SIGMA\n\n    # --- #4: 仕上げの全体係数（わずかに拡張）---\n    #sigma_air *= 1.04\n    return sigma_air\n\n\n# ======== 前処理パイプライン ========\nclass SignalProcessor:\n    def __init__(self, config):\n        self.cfg = config\n        self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\n        self.planet_ids = pd.read_csv(f'{self.cfg.DATA_PATH}/{self.cfg.DATASET}_star_info.csv', index_col='planet_id').index.astype(int)\n\n    def _apply_linear_corr(self, linear_corr, signal):\n        coeffs = np.flip(linear_corr, axis=0)\n        x = signal.astype(np.float64, copy=False)\n        out = np.empty_like(x, dtype=np.float64)\n        out[...] = coeffs[0]\n        for k in range(1, coeffs.shape[0]):\n            np.multiply(out, x, out=out)\n            out += coeffs[k]\n        return out.astype(signal.dtype, copy=False)\n\n    def _calibrate_single_signal(self, planet_id, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n    \n        signal = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_signal_0.parquet\"\n        ).to_numpy()\n        dark = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dark.parquet\"\n        ).to_numpy()\n        dead = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dead.parquet\"\n        ).to_numpy()\n        flat = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/flat.parquet\"\n        ).to_numpy()\n        linear_corr = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\"\n        ).values.astype(np.float64).reshape(sensor_cfg[\"linear_corr_shape\"])\n    \n        signal = signal.reshape(sensor_cfg[\"raw_shape\"])\n        gain = self.adc_info[f\"{sensor}_adc_gain\"].iloc[0]\n        offset = self.adc_info[f\"{sensor}_adc_offset\"].iloc[0]\n        signal = signal / gain + offset\n    \n        hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    \n        if sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            linear_corr = linear_corr[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dark = dark[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dead = dead[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            flat = flat[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            hot = hot[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n    \n        if sensor == \"FGS1\":\n            y0, y1, x0, x1 = 10, 22, 10, 22\n            signal = signal[:, y0:y1, x0:x1]\n            dark   = dark[y0:y1, x0:x1]\n            dead   = dead[y0:y1, x0:x1]\n            flat   = flat[y0:y1, x0:x1]\n            linear_corr = linear_corr[:, y0:y1, x0:x1]\n            hot    = hot[y0:y1, x0:x1]\n    \n        np.maximum(signal, 0, out=signal)\n    \n        if sensor == \"FGS1\":\n            signal = self._apply_linear_corr(linear_corr, signal)\n        elif sensor == \"AIRS-CH0\":\n            sl = (slice(None), slice(10, 22), slice(None))\n            signal[sl] = self._apply_linear_corr(linear_corr[:, 10:22, :], signal[sl])\n        else:\n            signal = self._apply_linear_corr(linear_corr, signal)\n    \n        base_dt, increment = sensor_cfg[\"dt_pattern\"]\n        even_scale = base_dt\n        odd_scale  = base_dt + increment\n        signal[::2]  -= dark * even_scale\n        signal[1::2] -= dark * odd_scale\n    \n        if sensor == \"FGS1\":\n            flat_roi = flat.astype(signal.dtype, copy=False).copy()\n            bad = (dead) | ~np.isfinite(flat_roi) | (flat_roi == 0)\n            flat_roi[bad] = np.nan\n            signal /= flat_roi\n    \n        elif sensor == \"AIRS-CH0\":\n            y0, y1 = 10, 22\n            flat_roi = flat[y0:y1, :].astype(signal.dtype, copy=False).copy()\n            bad = (dead[y0:y1, :]) | ~np.isfinite(flat_roi) | (flat_roi == 0)\n            flat_roi[bad] = np.nan\n            signal[:, y0:y1, :] /= flat_roi\n    \n        else:\n            flat2 = flat.astype(signal.dtype, copy=False).copy()\n            bad2 = (dead) | ~np.isfinite(flat2) | (flat2 == 0)\n            flat2[bad2] = np.nan\n            signal /= flat2\n    \n        return signal\n\n    def _preprocess_calibrated_signal(self, calibrated_signal, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n        binning = sensor_cfg[\"binning\"]\n\n        if sensor == \"AIRS-CH0\":\n            signal_roi = calibrated_signal[:, 10:22, :]\n        elif sensor == \"FGS1\":\n            signal_roi = calibrated_signal[:, 10:22, 10:22]\n            signal_roi = signal_roi.reshape(signal_roi.shape[0], -1)\n        \n        mean_signal = np.nanmean(signal_roi, axis=1)\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        n_bins = cds_signal.shape[0] // binning\n        binned = np.array([\n            cds_signal[j*binning : (j+1)*binning].mean(axis=0) \n            for j in range(n_bins)\n        ])\n\n        if sensor == \"AIRS-CH0\":\n            q_lo = np.nanpercentile(binned, 5.0, axis=1, keepdims=True)\n            q_hi = np.nanpercentile(binned, 95.0, axis=1, keepdims=True)\n            np.clip(binned, q_lo, q_hi, out=binned)\n\n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        if sensor == \"AIRS-CH0\":\n            var = np.nanvar(binned, axis=0, ddof=1)\n            med = np.nanmedian(var)\n            safe_var = np.where(~np.isfinite(var) | (var <= 0), med if (np.isfinite(med) and med > 0) else 1.0, var)\n            w = 1.0 / safe_var\n\n            lo, hi = np.nanpercentile(w, 5.0), np.nanpercentile(w, 95.0)\n            if np.isfinite(lo) and np.isfinite(hi) and lo < hi:\n                w = np.clip(w, lo, hi)\n\n            M = binned.shape[1]\n            s = np.nansum(w)\n            if np.isfinite(s) and s > 0:\n                w = w * (M / s)\n            else:\n                w = np.ones_like(w)\n\n            binned *= w[None, :]\n\n        return binned\n\n    def _process_planet_sensor(self, args):\n        planet_id, sensor = args['planet_id'], args['sensor']\n        calibrated = self._calibrate_single_signal(planet_id, sensor)\n        preprocessed = self._preprocess_calibrated_signal(calibrated, sensor)\n        return preprocessed\n\n    def process_all_data(self):\n        args_fgs1 = [dict(planet_id=planet_id, sensor=\"FGS1\") for planet_id in self.planet_ids]\n        preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        args_airs_ch0 = [dict(planet_id=planet_id, sensor=\"AIRS-CH0\") for planet_id in self.planet_ids]\n        preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        preprocessed_signal = np.concatenate(\n            [np.stack(preprocessed_fgs1), np.stack(preprocessed_airs_ch0)], axis=2\n        )\n        return preprocessed_signal\n    \n\n# ======== 1Dトランジット深さ推定 ========\nclass TransitModel:\n    def __init__(self, config):\n        self.cfg = config\n\n    def _phase_detector(self, signal):\n        search_slice = self.cfg.MODEL_PHASE_DETECTION_SLICE\n        min_index = np.argmin(signal[search_slice]) + search_slice.start\n        \n        signal1 = signal[:min_index]\n        signal2 = signal[min_index:]\n\n        grad1 = np.gradient(signal1)\n        grad1 /= grad1.max()\n        \n        grad2 = np.gradient(signal2)\n        grad2 /= grad2.max()\n\n        phase1 = np.argmin(grad1)\n        phase2 = np.argmax(grad2) + min_index\n\n        return phase1, phase2\n    \n    def _objective_function(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\n        coeffs = np.polyfit(x, y, deg=power)\n        poly = np.poly1d(coeffs)\n        error = np.abs(poly(x) - y).mean()\n        \n        return error\n\n    def predict(self, single_preprocessed_signal):\n        signal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)\n        signal_1d = savgol_filter(signal_1d, 23, 2)   # #2 は今回は触らない\n        \n        phase1, phase2 = self._phase_detector(signal_1d)\n\n        phase1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, phase1)\n        phase2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, phase2)    \n\n        result = minimize(\n            fun=self._objective_function,\n            x0=[0.0001],\n            args=(signal_1d, phase1, phase2),\n            method=\"Nelder-Mead\"\n        )\n        \n        return result.x[0]\n\n    def predict_all(self, preprocessed_signals):\n        predictions = [\n            self.predict(preprocessed_signal)\n            for preprocessed_signal in tqdm(preprocessed_signals)\n        ]\n        return np.array(predictions) * self.cfg.SCALE\n\n\n# ======== メタデータ ========\nStarInfo = pd.read_csv(ROOT_PATH + f\"/{MODE}_star_info.csv\")\nStarInfo[\"planet_id\"] = StarInfo[\"planet_id\"].astype(int)\nPlanetIds = StarInfo[\"planet_id\"].tolist()\nStarInfo = StarInfo.set_index(\"planet_id\")\n\n\n# ======== AIRS 出力（282次元）モデル ========\nclass ResidualBlock2(nn.Module):\n    def __init__(self, dim, p=0.2):\n        super().__init__()\n        self.fc1 = nn.Linear(dim, dim)\n        self.bn1 = nn.BatchNorm1d(dim)\n        self.fc2 = nn.Linear(dim, dim)\n        self.bn2 = nn.BatchNorm1d(dim)\n        self.relu = nn.ReLU()\n        self.dropout = nn.Dropout(p)\n\n    def forward(self, x):\n        identity = x\n        out = self.relu(self.bn1(self.fc1(x)))\n        out = self.dropout(out)\n        out = self.bn2(self.fc2(out))\n        return self.relu(out + identity)\n\n\nclass ResNetMLP2(nn.Module):\n    def __init__(self, input_dim=3, hidden_dim=128, output_dim=282, num_blocks=3, dropout_rate=0.2):\n        super().__init__()\n        self.input_layer = nn.Linear(input_dim, hidden_dim)\n        self.blocks = nn.Sequential(*[ResidualBlock2(hidden_dim, p=dropout_rate) for _ in range(num_blocks)])\n        self.output_layer = nn.Linear(hidden_dim, output_dim)\n\n    def forward(self, x):\n        x = self.input_layer(x)\n        x = self.blocks(x)\n        x = self.output_layer(x)\n        return x\n\n\n# ======== 提出生成 ========\nclass SubmissionGenerator:\n    def __init__(self, config):\n        self.cfg = config\n        self.sample_submission = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", index_col=\"planet_id\")\n\n    def create(self, predictions1, predictions, sigma_fgs=None, sigma_air=None):\n        planet_ids = self.sample_submission.index\n        n_mu = self.sample_submission.shape[1] // 2  # 283\n\n        preds = np.asarray(predictions, dtype=float).reshape(-1)\n        mu = np.tile(preds.reshape(-1, 1), (1, n_mu))\n        mu = np.clip(mu, 0, None)\n\n        sigmas = np.full_like(mu, self.cfg.SIGMA, dtype=float)\n        if sigma_fgs is not None:\n            sigma_fgs = np.asarray(sigma_fgs, dtype=float).reshape(-1)\n            sigmas[:, 0] = np.clip(sigma_fgs, 1e-6, 0.1)\n        if sigma_air is not None:\n            sigma_air = np.asarray(sigma_air, dtype=float).reshape(-1, 1)\n            sigmas[:, 1:] = np.clip(sigma_air, 1e-6, 0.1)\n\n        submission_df = pd.DataFrame(\n            np.concatenate([mu, sigmas], axis=1),\n            columns=self.sample_submission.columns,\n            index=planet_ids\n        )\n        # μの1列目(白色FGS)と AIRS 1:282 を元の仕様どおり再配置\n        submission_df.iloc[:, 0] = predictions\n        submission_df.iloc[:, 1:283] = predictions1\n\n        #submission_df.to_csv(\"submission.csv\")\n        return submission_df\n\n\n# =========================\n# 実行フロー\n# =========================\nif __name__ == \"__main__\":\n    config = Config()\n    \n    # 前処理\n    signal_processor = SignalProcessor(config)\n    preprocessed_data = signal_processor.process_all_data()\n\n    # 1D深さ\n    model = TransitModel(config)\n    predictions = model.predict_all(preprocessed_data)\n\n    # --- #4: 保守拡張版の σ 推定 ---\n    sigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)\n    sigma_air_vec = estimate_sigma_air(preprocessed_data, config)\n\n    # メタ特徴作成\n    predictions_df = pd.DataFrame({\n        \"planet_id\": PlanetIds,\n        \"transit_depth\": predictions\n    })\n    input_df = pd.merge(predictions_df, StarInfo, on=\"planet_id\", how=\"left\")\n    input_df[\"transit_depth\"] *= 10000\n    features = ['transit_depth','Rs','i']\n    X = input_df[features].values.astype(np.float32)\n    X_tensor = torch.tensor(X, dtype=torch.float32)\n\n    # AIRS 予測モデル読み込み\n    resnet2 = ResNetMLP2(num_blocks=80, dropout_rate=0.3)\n    resnet2.load_state_dict(torch.load(\"/kaggle/input/airs/pytorch/default/1/best_model_airs.pth\", map_location=\"cpu\"))\n    resnet2.eval()\n\n    with torch.no_grad():\n        predictions1 = resnet2(X_tensor).numpy()\n    predictions1 /= 10000\n\n    # --- #5: AIRS 出力の軽スムージング・アンサンブル ---\n    # 波長方向(軸=1)に Savitzky–Golay で極軽く平滑化し、0.7:0.3 でブレンド\n    predictions1_smooth = savgol_filter(predictions1, window_length=13, polyorder=2, axis=1)\n    predictions1 = 0.65 * predictions1 + 0.35 * predictions1_smooth\n\n    submission_generator = SubmissionGenerator(config)\n    submission_solution4 = submission_generator.create(\n        predictions1,\n        predictions,\n        sigma_fgs=sigma_fgs_vec,\n        sigma_air=sigma_air_vec\n    )","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:29:39.583801Z","iopub.execute_input":"2025-09-24T08:29:39.584304Z","iopub.status.idle":"2025-09-24T08:29:45.553794Z","shell.execute_reply.started":"2025-09-24T08:29:39.584151Z","shell.execute_reply":"2025-09-24T08:29:45.5527Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_solution4.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:29:45.554786Z","iopub.execute_input":"2025-09-24T08:29:45.555168Z","iopub.status.idle":"2025-09-24T08:29:45.576892Z","shell.execute_reply.started":"2025-09-24T08:29:45.555131Z","shell.execute_reply":"2025-09-24T08:29:45.57565Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution 5","metadata":{}},{"cell_type":"code","source":"ROOT_PATH = \"/kaggle/input/ariel-data-challenge-2025\"\nMODE = \"test\"\n\n__t0 = time.perf_counter()\n\nclass Config:\n    FEATURES = ['transit_depth', 'Rs', 'i']\n    # FEATURES = ['transit_depth', 'Rs', 'i', 'P', 'scaled_depth']\n    # FEATURES = ['transit_depth', 'Rs', 'Ms', 'Ts', 'Mp', 'e', 'P', 'sma', 'i']\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n    load_data = False  # Set to True to load from disk, False to generate\n    DEBUG = False\n\n    SCALE = 0.96\n    SIGMA = 0.00055\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA = 11 # 9\n    MODEL_POLYNOMIAL_DEGREE = 3\n    \n    N_JOBS = 3\n\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    row_id_column_name: str,\n    naive_mean: float,\n    naive_sigma: float,\n    fsg_sigma_true: float = 1e-6,\n    airs_sigma_true: float = 1e-5,\n    fgs_weight: float = 1,\n) -> float:\n    \"\"\"\n    This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n    we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n    x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n    the GLL value for one spectrum is the sum of all of them.\n\n    Inputs:\n        - solution: Ground Truth spectra (from test set)\n            - shape: (nsamples, n_wavelengths)\n        - submission: Predicted spectra and errors (from participants)\n            - shape: (nsamples, n_wavelengths*2)\n        naive_mean: (float) mean from the train set.\n        naive_sigma: (float) standard deviation from the train set.\n        fsg_sigma_true: (float) standard deviation from the FSG1 instrument for the test set.\n        airs_sigma_true: (float) standard deviation from the AIRS instrument for the test set.\n        fgs_weight: (float) relative weight of the fgs channel\n    \"\"\"\n\n    del solution[row_id_column_name]\n    del submission[row_id_column_name]\n\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\n    for col in submission.columns:\n        if not pandas.api.types.is_numeric_dtype(submission[col]):\n            raise ParticipantVisibleError(f'Submission column {col} must be a number')\n\n    n_wavelengths = len(solution.columns)\n    if len(submission.columns) != n_wavelengths * 2:\n        raise ParticipantVisibleError('Wrong number of columns in the submission')\n\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=10**-15, a_max=None)\n    sigma_true = np.append(\n        np.array(\n            [\n                fsg_sigma_true,\n            ]\n        ),\n        np.ones(n_wavelengths - 1) * airs_sigma_true,\n    )\n    y_true = solution.values\n\n    GLL_pred = scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred)\n    GLL_true = scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true * np.ones_like(y_true))\n    GLL_mean = scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true))\n\n    # normalise the score, right now it becomes a matrix instead of a scalar.\n    ind_scores = (GLL_pred - GLL_mean) / (GLL_true - GLL_mean)\n\n    weights = np.append(np.array([fgs_weight]), np.ones(len(solution.columns) - 1))\n    weights = weights * np.ones_like(ind_scores)\n    submit_score = np.average(ind_scores, weights=weights)\n    return float(np.clip(submit_score, 0.0, 1.0))\n\ndef _phase_detector_signal(signal, cfg):\n    sl = cfg.MODEL_PHASE_DETECTION_SLICE\n    min_idx = int(np.argmin(signal[sl])) + sl.start\n    s1 = signal[:min_idx]; s2 = signal[min_idx:]\n    if s1.size < 3 or s2.size < 3:\n        return 0, len(signal) - 1\n    g1 = np.gradient(s1); g1_max = np.max(g1) if np.size(g1) else 0.0\n    g2 = np.gradient(s2); g2_max = np.max(g2) if np.size(g2) else 0.0\n    if g1_max != 0: g1 /= g1_max\n    if g2_max != 0: g2 /= g2_max\n    phase1 = int(np.argmin(g1)); phase2 = int(np.argmax(g2)) + min_idx\n    return phase1, phase2\n\ndef estimate_sigma_fgs(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_1 (для FGS1) длиной N_planets — мягкий множитель к cfg.SIGMA.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n    for single in preprocessed_data:\n        # фазы по AIRS белой кривой — так же, как в модели\n        air_white = savgol_filter(single[:, 1:].mean(axis=1), 20, 2)\n        p1, p2 = _phase_detector_signal(air_white, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(air_white) - delta - 1, p2)\n\n        fgs = single[:, 0]\n        oot = (fgs[: p1 - delta] if p1 - delta > 0 else np.empty(0, fgs.dtype))\n        if p2 + delta < fgs.size:\n            oot = np.concatenate([oot, fgs[p2 + delta :]])\n        inn = fgs[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(fgs))\n        # относительная неопределённость глубины (в тех же ед., что s)\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель: корень, и узкий клип, чтобы не рисковать\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.8, 1.25)  # ±20–25% от базовой σ\n\n    return k * cfg.SIGMA\n\ndef estimate_sigma_air(preprocessed_data, cfg):\n    \"\"\"Возвращает вектор sigma_air длиной N_planets — мягкий множитель к cfg.SIGMA для всех AIRS-каналов.\"\"\"\n    sig_rel = []\n    delta = cfg.MODEL_OPTIMIZATION_DELTA\n    eps = 1e-12\n\n    for single in preprocessed_data:\n        # белая кривая AIRS на бинированных данных (после всех твоих весов по λ)\n        white = np.nanmean(single[:, 1:], axis=1)         # (n_bins,)\n        white_s = savgol_filter(white, 20, 2)             # для фаз\n\n        p1, p2 = _phase_detector_signal(white_s, cfg)\n        p1 = max(delta, p1)\n        p2 = min(len(white) - delta - 1, p2)\n\n        oot_left = white[: p1 - delta] if p1 - delta > 0 else np.empty(0, white.dtype)\n        oot_right = white[p2 + delta :] if (p2 + delta) < white.size else np.empty(0, white.dtype)\n        oot = np.concatenate([oot_left, oot_right]) if (oot_left.size + oot_right.size) else oot_left\n        inn = white[p1 + delta : max(p1 + delta, p2 - delta)]\n\n        if oot.size == 0 or inn.size == 0:\n            sig_rel.append(np.nan); continue\n\n        n_oot, n_in = len(oot), len(inn)\n        var_oot = np.nanvar(oot, ddof=1)\n        var_in  = np.nanvar(inn, ddof=1)\n        oot_mean = float(np.nanmean(oot)) if np.isfinite(np.nanmean(oot)) else float(np.nanmean(white))\n\n        sigma_rel = np.sqrt(var_oot / max(n_oot,1) + var_in / max(n_in,1)) / max(oot_mean, eps)\n        sig_rel.append(sigma_rel)\n\n    s = np.asarray(sig_rel, dtype=float)\n    mask = np.isfinite(s) & (s > 0)\n    med = float(np.nanmedian(s[mask])) if mask.any() else 1.0\n\n    # мягкий множитель вокруг медианы\n    k = np.ones_like(s)\n    if med > 0 and np.isfinite(med):\n        k[mask] = np.sqrt(s[mask] / med)\n    k = np.clip(k, 0.90, 1.20)  # ±10%–20%\n\n    return k * cfg.SIGMA\n\nclass SignalProcessor:\n    def __init__(self, config):\n        self.cfg = config\n        self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\n        self.planet_ids = pd.read_csv(f'{self.cfg.DATA_PATH}/{self.cfg.DATASET}_star_info.csv', index_col='planet_id').index.astype(int)\n\n    def _apply_linear_corr(self, linear_corr, signal):\n\n        coeffs = np.flip(linear_corr, axis=0)      # shape: (D, X, Y), D — старшая степень сначала\n        x = signal.astype(np.float64, copy=False)  # считаем в float64 для стабильности\n        out = np.empty_like(x, dtype=np.float64)\n        out[...] = coeffs[0]  # broadcast (X,Y) -> (T,X,Y)\n        for k in range(1, coeffs.shape[0]):\n            np.multiply(out, x, out=out)  # in-place умножение\n            out += coeffs[k]              # broadcast (X,Y)\n\n        return out.astype(signal.dtype, copy=False)\n\n    def _calibrate_single_signal(self, planet_id, sensor):\n        \"\"\"\n        Калибровка single-node сигнала.\n        Политика масок: DEAD — маскируем, HOT — НЕ маскируем (оставляем в данных).\n        \"\"\"\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n    \n        # --- load ---\n        signal = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_signal_0.parquet\"\n        ).to_numpy()\n        dark = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dark.parquet\"\n        ).to_numpy()\n        dead = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/dead.parquet\"\n        ).to_numpy()\n        flat = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/flat.parquet\"\n        ).to_numpy()\n        linear_corr = pd.read_parquet(\n            f\"{self.cfg.DATA_PATH}/{self.cfg.DATASET}/{planet_id}/{sensor}_calibration_0/linear_corr.parquet\"\n        ).values.astype(np.float64).reshape(sensor_cfg[\"linear_corr_shape\"])\n    \n        # --- reshape & ADC ---\n        signal = signal.reshape(sensor_cfg[\"raw_shape\"])\n        gain = self.adc_info[f\"{sensor}_adc_gain\"].iloc[0]\n        offset = self.adc_info[f\"{sensor}_adc_offset\"].iloc[0]\n        signal = signal / gain + offset  # сохраняем твою формулу\n    \n        # HOT только для мониторинга, не для маскирования\n        hot = sigma_clip(dark, sigma=5, maxiters=5).mask\n    \n        # --- crop per sensor ---\n        if sensor == \"AIRS-CH0\":\n            signal = signal[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            linear_corr = linear_corr[:, :, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dark = dark[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            dead = dead[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            flat = flat[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]\n            hot = hot[:, self.cfg.CUT_INF : self.cfg.CUT_SUP]  # только для логов\n    \n        if sensor == \"FGS1\":\n            y0, y1, x0, x1 = 10, 22, 10, 22\n            signal = signal[:, y0:y1, x0:x1]\n            dark   = dark[y0:y1, x0:x1]\n            dead   = dead[y0:y1, x0:x1]\n            flat   = flat[y0:y1, x0:x1]\n            linear_corr = linear_corr[:, y0:y1, x0:x1]\n            hot    = hot[y0:y1, x0:x1]  # только для логов\n    \n        # --- non-neg clamp before linearity corr (как у тебя) ---\n        np.maximum(signal, 0, out=signal)\n    \n        # --- linearity correction ---\n        if sensor == \"FGS1\":\n            signal = self._apply_linear_corr(linear_corr, signal)\n        elif sensor == \"AIRS-CH0\":\n            sl = (slice(None), slice(10, 22), slice(None))  # T, Y, λ\n            signal[sl] = self._apply_linear_corr(linear_corr[:, 10:22, :], signal[sl])\n        else:\n            signal = self._apply_linear_corr(linear_corr, signal)\n    \n        # --- dark subtraction с учётом паттерна интеграций ---\n        base_dt, increment = sensor_cfg[\"dt_pattern\"]\n        even_scale = base_dt\n        odd_scale  = base_dt + increment\n        signal[::2]  -= dark * even_scale\n        signal[1::2] -= dark * odd_scale\n    \n        # --- APPLY FLAT (HOT-KEEP: не включаем hot в маску!) ---\n        if sensor == \"FGS1\":\n            flat_roi = flat.astype(signal.dtype, copy=False).copy()      # (12,12)\n            bad = (dead) | ~np.isfinite(flat_roi) | (flat_roi == 0)      # ← ТОЛЬКО dead/invalid\n            flat_roi[bad] = np.nan\n            signal /= flat_roi\n    \n        elif sensor == \"AIRS-CH0\":\n            y0, y1 = 10, 22\n            flat_roi = flat[y0:y1, :].astype(signal.dtype, copy=False).copy()  # (12, λ)\n            bad = (dead[y0:y1, :]) | ~np.isfinite(flat_roi) | (flat_roi == 0)  # ← ТОЛЬКО dead/invalid\n            flat_roi[bad] = np.nan\n            signal[:, y0:y1, :] /= flat_roi\n    \n        else:\n            flat2 = flat.astype(signal.dtype, copy=False).copy()\n            bad2 = (dead) | ~np.isfinite(flat2) | (flat2 == 0)                  # ← ТОЛЬКО dead/invalid\n            flat2[bad2] = np.nan\n            signal /= flat2\n        # --- END FLAT ---\n    \n        # (опционально) логируем метрики hot/dead\n        if getattr(self.cfg, \"LOG_HOT_STATS\", False):\n            if not hasattr(self, \"stats\"):\n                self.stats = []\n            self.stats.append({\n                \"planet_id\": int(planet_id),\n                \"sensor\": sensor,\n                \"hot_frac\": float(np.mean(hot)),\n                \"dead_frac\": float(np.mean(dead)),\n            })\n    \n        return signal\n\n    def _preprocess_calibrated_signal(self, calibrated_signal, sensor):\n        sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n        binning = sensor_cfg[\"binning\"]\n\n        if sensor == \"AIRS-CH0\":\n            signal_roi = calibrated_signal[:, 10:22, :]\n        elif sensor == \"FGS1\":\n            signal_roi = calibrated_signal[:, 10:22, 10:22]\n            signal_roi = signal_roi.reshape(signal_roi.shape[0], -1)\n        \n        mean_signal = np.nanmean(signal_roi, axis=1)\n\n        cds_signal = mean_signal[1::2] - mean_signal[0::2]\n\n        n_bins = cds_signal.shape[0] // binning\n        binned = np.array([\n            cds_signal[j*binning : (j+1)*binning].mean(axis=0) \n            for j in range(n_bins)\n        ])\n\n        if sensor == \"AIRS-CH0\":\n            q_lo = np.nanpercentile(binned, 5.0, axis=1, keepdims=True)    # (n_bins, 1)\n            q_hi = np.nanpercentile(binned, 95.0, axis=1, keepdims=True)   # (n_bins, 1)\n            np.clip(binned, q_lo, q_hi, out=binned)\n\n        if sensor == \"FGS1\":\n            binned = binned.reshape((binned.shape[0], 1))\n\n        if sensor == \"AIRS-CH0\":\n            var = np.nanvar(binned, axis=0, ddof=1)                 # (λ, )\n            med = np.nanmedian(var)\n            safe_var = np.where(~np.isfinite(var) | (var <= 0), med if (np.isfinite(med) and med > 0) else 1.0, var)\n            w = 1.0 / safe_var\n\n            lo, hi = np.nanpercentile(w, 5.0), np.nanpercentile(w, 95.0)\n            if np.isfinite(lo) and np.isfinite(hi) and lo < hi:\n                w = np.clip(w, lo, hi)\n\n            M = binned.shape[1]\n            s = np.nansum(w)\n            if np.isfinite(s) and s > 0:\n                w = w * (M / s)\n            else:\n                w = np.ones_like(w)\n\n            binned *= w[None, :]\n\n\n        return binned\n\n    def _process_planet_sensor(self, args):\n        planet_id, sensor = args['planet_id'], args['sensor']\n        calibrated = self._calibrate_single_signal(planet_id, sensor)\n        preprocessed = self._preprocess_calibrated_signal(calibrated, sensor)\n        return preprocessed\n\n    def process_all_data(self):\n        args_fgs1 = [dict(planet_id=planet_id, sensor=\"FGS1\") for planet_id in self.planet_ids]\n        preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        args_airs_ch0 = [dict(planet_id=planet_id, sensor=\"AIRS-CH0\") for planet_id in self.planet_ids]\n        preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n\n        preprocessed_signal = np.concatenate(\n            [np.stack(preprocessed_fgs1), np.stack(preprocessed_airs_ch0)], axis=2\n        )\n        return preprocessed_signal\n    \n\nclass TransitModel:\n    def __init__(self, config):\n        self.cfg = config\n\n    def _phase_detector(self, signal):\n        search_slice = self.cfg.MODEL_PHASE_DETECTION_SLICE\n        min_index = np.argmin(signal[search_slice]) + search_slice.start\n        \n        signal1 = signal[:min_index]\n        signal2 = signal[min_index:]\n\n        grad1 = np.gradient(signal1)\n        grad1 /= grad1.max()\n        \n        grad2 = np.gradient(signal2)\n        grad2 /= grad2.max()\n\n        phase1 = np.argmin(grad1)\n        phase2 = np.argmax(grad2) + min_index\n\n        return phase1, phase2\n    \n    def _objective_function(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\n        coeffs = np.polyfit(x, y, deg=power)\n        poly = np.poly1d(coeffs)\n        error = np.abs(poly(x) - y).mean()\n        \n        return error\n\n    def predict(self, single_preprocessed_signal):\n        signal_1d = single_preprocessed_signal[:, 1:].mean(axis=1)\n        signal_1d = savgol_filter(signal_1d, 23, 2)\n        \n        phase1, phase2 = self._phase_detector(signal_1d)\n\n        phase1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, phase1)\n        phase2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, phase2)    \n\n        result = minimize(\n            fun=self._objective_function,\n            x0=[0.0001],\n            args=(signal_1d, phase1, phase2),\n            method=\"Nelder-Mead\"\n        )\n        \n        return result.x[0]\n\n    def predict_all(self, preprocessed_signals):\n        predictions = [\n            self.predict(preprocessed_signal)\n            for preprocessed_signal in tqdm(preprocessed_signals)\n        ]\n        return np.array(predictions) * self.cfg.SCALE\n\nStarInfo = pd.read_csv(ROOT_PATH + f\"/{MODE}_star_info.csv\")\nStarInfo[\"planet_id\"] = StarInfo[\"planet_id\"].astype(int)\nPlanetIds = StarInfo[\"planet_id\"].tolist()\nStarInfo = StarInfo.set_index(\"planet_id\")\nclass SubmissionGenerator:\n    def __init__(self, config):\n        self.cfg = config\n        self.sample_submission = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", index_col=\"planet_id\")\n\n    def create(self, predictions1, predictions, sigma_fgs=None, sigma_air=None):\n        planet_ids = self.sample_submission.index\n        n_mu = self.sample_submission.shape[1] // 2  # 283\n\n        preds = np.asarray(predictions, dtype=float).reshape(-1)\n        mu = np.tile(preds.reshape(-1, 1), (1, n_mu))\n        mu = np.clip(mu, 0, None)\n\n        sigmas = np.full_like(mu, self.cfg.SIGMA, dtype=float)\n        if sigma_fgs is not None:\n            sigma_fgs = np.asarray(sigma_fgs, dtype=float).reshape(-1)\n            sigmas[:, 0] = np.clip(sigma_fgs, 1e-6, 0.1)\n        if sigma_air is not None:\n            sigma_air = np.asarray(sigma_air, dtype=float).reshape(-1, 1)\n            sigmas[:, 1:] = np.clip(sigma_air, 1e-6, 0.1)\n\n        submission_df = pd.DataFrame(\n            np.concatenate([mu, sigmas], axis=1),\n            columns=self.sample_submission.columns,\n            index=planet_ids\n        )\n        submission_df.iloc[:, 0] = predictions\n        submission_df.iloc[:, 1:283] = predictions1\n        submission_df.to_csv(\"submission.csv\")\n        \n        return submission_df\n\nclass SEBlock(nn.Module):\n    def __init__(self, dim, reduction=16):\n        super().__init__()\n        self.fc1 = nn.Linear(dim, dim // reduction, bias=False)\n        self.fc2 = nn.Linear(dim // reduction, dim, bias=False)\n        self.act = nn.SiLU()   # Or nn.ReLU()\n\n    def forward(self, x):\n        # Compute channel-wise attention\n        w = x.mean(dim=0, keepdim=True)      # Global context (mean across batch)\n        w = self.act(self.fc1(w))\n        w = torch.sigmoid(self.fc2(w))\n        return x * w   # Rescale input\n\nclass AttentionBlock(nn.Module):\n    def __init__(self, dim, num_heads=4, dropout=0.1):\n        super().__init__()\n        self.attn = nn.MultiheadAttention(embed_dim=dim, num_heads=num_heads, batch_first=True)\n        self.norm = nn.LayerNorm(dim)\n        self.dropout = nn.Dropout(dropout)\n\n    def forward(self, x):\n        # Expect shape [batch, seq_len, dim], so expand if only [batch, dim]\n        if x.dim() == 2:\n            x = x.unsqueeze(1)   # -> [batch, 1, dim]\n        attn_out, _ = self.attn(x, x, x)\n        out = self.norm(x + self.dropout(attn_out))\n        return out.squeeze(1)    # Back to [batch, dim]\n\nclass ResidualBlock2(nn.Module):\n    def __init__(self, dim, p=0.2):\n        super().__init__()\n        self.fc1 = nn.Linear(dim, dim)\n        self.fc2 = nn.Linear(dim, dim)\n        self.relu = nn.ReLU()\n        self.dropout = nn.Dropout(p)\n\n    def forward(self, x):\n        identity = x\n        out = self.relu(self.fc1(x))\n        out = self.dropout(out)\n        out = self.fc2(out)\n        return self.relu(out + identity)\n\nclass ResNetMLP2(nn.Module):\n    def __init__(self, input_dim=3, hidden_dim=128, output_dim=282, num_blocks=3, dropout_rate=0.2):\n        super().__init__()\n        self.input_layer = nn.Linear(input_dim, hidden_dim)\n        self.blocks = nn.Sequential(*[ResidualBlock2(hidden_dim, p=dropout_rate) for _ in range(num_blocks)])\n        self.output_layer = nn.Linear(hidden_dim, output_dim)\n\n    def forward(self, x):\n        x = self.input_layer(x)\n        x = self.blocks(x)\n        x = self.output_layer(x)\n        return x\n\n\ndef load_cv_models_and_scalers(directory):\n    \"\"\"\n    Loads all cross-validation models and scalers from the specified directory.\n    Args:\n        directory (str): Path to the directory containing model and scaler files.\n    Returns:\n        all_models (list): List of loaded models.\n        scaler_X: Loaded X scaler.\n        scaler_y: Loaded y scaler.\n    \"\"\"\n    import os\n    import joblib\n    # Load scalers\n    scaler_X = joblib.load(os.path.join(directory, 'scaler_X.joblib'))\n    scaler_y = joblib.load(os.path.join(directory, 'scaler_y.joblib'))\n\n    # Load all CV models\n    all_models = []\n    model_params = {\n        'input_dim': len(Config.FEATURES),  # Always use the current feature count\n        'hidden_dim': 256,\n        'output_dim': 282,\n        'num_blocks': 35,\n        'dropout_rate': 0.1\n    }\n    for fold in range(1, 11):\n        model = ResNetMLP2(**model_params).double()\n        model_path = os.path.join(directory, f'best_model_airs_cv_fold{fold}.pth')\n        model.load_state_dict(torch.load(model_path, map_location=torch.device('cpu')))\n        model.eval()\n        all_models.append(model)\n    return all_models, scaler_X, scaler_y\n    \nall_models, scaler_X, scaler_y = load_cv_models_and_scalers('/kaggle/input/d/pankajiitr/ariel-2025-result/results_v37')\n\nconfig = Config()\nsignal_processor = SignalProcessor(config)\npreprocessed_data = signal_processor.process_all_data()\n\nmodel = TransitModel(config)\npredictions = model.predict_all(preprocessed_data)\nsigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)\nsigma_air_vec = estimate_sigma_air(preprocessed_data, config)\n\npredictions_df = pd.DataFrame({\n    \"planet_id\": PlanetIds,\n    \"transit_depth\": predictions\n})\n\ninput_df = pd.merge(predictions_df, StarInfo, on=\"planet_id\", how=\"left\")\ninput_df['scaled_depth'] = input_df['transit_depth'] / input_df['Rs']\n\nX = input_df[Config.FEATURES].values.astype(np.float64)\nX_scaled = scaler_X.transform(X)\nX_tensor = torch.tensor(X_scaled, dtype=torch.float64)\n\n# Generate average prediction from all CV models (on scaled X, then inverse transform)\nwith torch.no_grad():\n    preds_scaled = [model(X_tensor).numpy() for model in all_models]\npredictions1_scaled = np.mean(preds_scaled, axis=0)\npredictions1 = scaler_y.inverse_transform(predictions1_scaled)\n\nclass Config:\n    FEATURES = ['transit_depth', 'Rs', 'i', 'P']\n\n    DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n    DATASET = \"test\"\n    load_data = False  # Set to True to load from disk, False to generate\n    DEBUG = False\n\n    SCALE = 0.96\n    SIGMA = 0.00055\n    \n    CUT_INF = 39\n    CUT_SUP = 321\n    \n    SENSOR_CONFIG = {\n        \"AIRS-CH0\": {\n            \"raw_shape\": [11250, 32, 356],\n            \"calibrated_shape\": [1, 32, CUT_SUP - CUT_INF],\n            \"linear_corr_shape\": (6, 32, 356),\n            \"dt_pattern\": (0.1, 4.5), \n            \"binning\": 30\n        },\n        \"FGS1\": {\n            \"raw_shape\": [135000, 32, 32],\n            \"calibrated_shape\": [1, 32, 32],\n            \"linear_corr_shape\": (6, 32, 32),\n            \"dt_pattern\": (0.1, 0.1),\n            \"binning\": 30 * 12\n        }\n    }\n    \n    MODEL_PHASE_DETECTION_SLICE = slice(30, 140)\n    MODEL_OPTIMIZATION_DELTA = 11 # 9\n    MODEL_POLYNOMIAL_DEGREE = 3\n    \n    N_JOBS = 3\n\n\ndef load_cv_models_and_scalers(directory):\n    \"\"\"\n    Loads all cross-validation models and scalers from the specified directory.\n    Args:\n        directory (str): Path to the directory containing model and scaler files.\n    Returns:\n        all_models (list): List of loaded models.\n        scaler_X: Loaded X scaler.\n        scaler_y: Loaded y scaler.\n    \"\"\"\n    import os\n    import joblib\n    # Load scalers\n    scaler_X = joblib.load(os.path.join(directory, 'scaler_X.joblib'))\n    scaler_y = joblib.load(os.path.join(directory, 'scaler_y.joblib'))\n\n    # Load all CV models\n    all_models = []\n    model_params = {\n        'input_dim': len(Config.FEATURES),  # Always use the current feature count\n        'hidden_dim': 256,\n        'output_dim': 282,\n        'num_blocks': 35,\n        'dropout_rate': 0.1\n    }\n    for fold in range(1, 11):\n        model = ResNetMLP2(**model_params).double()\n        model_path = os.path.join(directory, f'best_model_airs_cv_fold{fold}.pth')\n        model.load_state_dict(torch.load(model_path, map_location=torch.device('cpu')))\n        model.eval()\n        all_models.append(model)\n    return all_models, scaler_X, scaler_y\n    \nall_models, scaler_X, scaler_y = load_cv_models_and_scalers('/kaggle/input/d/pankajiitr/ariel-2025-result/results_sv44')\nconfig = Config()\nsignal_processor = SignalProcessor(config)\npreprocessed_data = signal_processor.process_all_data()\n\nmodel = TransitModel(config)\npredictions = model.predict_all(preprocessed_data)\nsigma_fgs_vec = estimate_sigma_fgs(preprocessed_data, config)\nsigma_air_vec = estimate_sigma_air(preprocessed_data, config)\n\npredictions_df = pd.DataFrame({\n    \"planet_id\": PlanetIds,\n    \"transit_depth\": predictions\n})\n\ninput_df = pd.merge(predictions_df, StarInfo, on=\"planet_id\", how=\"left\")\nX = input_df[Config.FEATURES].values.astype(np.float64)\nX_scaled = scaler_X.transform(X)\nX_tensor = torch.tensor(X_scaled, dtype=torch.float64)\n\n# Generate average prediction from all CV models (on scaled X, then inverse transform)\nwith torch.no_grad():\n    preds_scaled = [model(X_tensor).numpy() for model in all_models]\n\n# Correctly average the predictions\npredictions2_scaled = np.mean(preds_scaled, axis=0)\n\n# Correctly inverse transform the *averaged* predictions\npredictions2 = scaler_y.inverse_transform(predictions2_scaled)\n\n# Combine the predictions into a single array\nall_predictions = np.array([predictions1, predictions2])\n\n# Calculate the mean across the models (axis=0)\nfinal_predictions = np.mean(all_predictions, axis=0)\n\nsubmission_generator = SubmissionGenerator(config)\nsubmission_solution5 = submission_generator.create(final_predictions, predictions, sigma_fgs=sigma_fgs_vec, sigma_air=sigma_air_vec)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:29:45.578549Z","iopub.execute_input":"2025-09-24T08:29:45.578858Z","iopub.status.idle":"2025-09-24T08:30:11.872885Z","shell.execute_reply.started":"2025-09-24T08:29:45.578838Z","shell.execute_reply":"2025-09-24T08:30:11.871847Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission_solution5.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:30:11.874222Z","iopub.execute_input":"2025-09-24T08:30:11.874565Z","iopub.status.idle":"2025-09-24T08:30:11.896403Z","shell.execute_reply.started":"2025-09-24T08:30:11.874542Z","shell.execute_reply":"2025-09-24T08:30:11.895039Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Solution_ML","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Soft Voting ","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\nweights = np.array([0.03, 0.40, 0.20, 0.20, 0.17], dtype=float)\nsolutions = [submission_solution1, submission_solution2, submission_solution3, submission_solution4, submission_solution5]\n\n# 可选：权重归一化\n#w_sum = weights.sum()\n#if not np.isclose(w_sum, 1.0):\n#    weights = weights / w_sum\n\n# 先加权平均\nweighted_avg = sum(w * df for w, df in zip(weights, solutions))\n\n# 再裁剪\nweighted_avg = weighted_avg.clip(lower=1e-7)\n\nweighted_avg.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:31:53.977236Z","iopub.execute_input":"2025-09-24T08:31:53.9776Z","iopub.status.idle":"2025-09-24T08:31:54.007451Z","shell.execute_reply.started":"2025-09-24T08:31:53.977577Z","shell.execute_reply":"2025-09-24T08:31:54.006012Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"weighted_avg.to_csv('submission.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T08:32:05.526714Z","iopub.execute_input":"2025-09-24T08:32:05.527063Z","iopub.status.idle":"2025-09-24T08:32:05.538Z","shell.execute_reply.started":"2025-09-24T08:32:05.527039Z","shell.execute_reply":"2025-09-24T08:32:05.536752Z"}},"outputs":[],"execution_count":null}]}