{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","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":13142515,"sourceType":"datasetVersion","datasetId":8321298},{"sourceId":247163208,"sourceType":"kernelVersion"}],"dockerImageVersionId":31040,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ariel '25 Baseline: FGS1->Features->LGBM\n\n* Wow - 263GB of data!\n* Just a single transit has 135,000 32x32 FGS1 frames!\n* Hard to figure out where to start....\n\nThe trick to machine learning seems to be making a new problem look like one you already know. So - let's try some feature engineering to **turn this all into tabular data and feed it into LGBM.**\n\n*Friendly Reminder:* If re-using large parts of this work in a public notebook - **please credit where you found the code**\n\n## Approach\n\nWe focus on **just the FGS1 data** - we do ADC correction but **ignore other calibration (dead-pixels, etc.)** for now.\n\nSince we have 283 spectral targets, we train *283 separate LGBM models per CV fold.* That's **1415 LGBM models total!**\n\nMultiple signals create multiple rows of tabular data for both train and test.  Planets with multiple signals result in multiple predictions (which are averaged together).\n\n# Engineered Features\n\nI spent some quality time with Claude to come up with a pipeline extracts **175 engineered features** from the FGS1 data:\n\n| Category | Count | Purpose |\n|----------|-------|---------|\n| **Global Flux Statistics** | 19 | Overall brightness characteristics & fundamental transit depth measurement |\n| **Rolling Statistics** | 54 | Multi-timescale transit detection across 6 window sizes (5 sec to 8.3 min) |\n| **Transit Detection** | 20 | Specialized transit identification, timing, duration & shape analysis |\n| **Frequency Domain** | 16 | Signal/noise separation & periodic pattern detection via FFT & autocorrelation |\n| **Spatial Features** | 30 | Star position/shape changes during transit across 5 key timepoints |\n| **Gradient & Change** | 36 | Rate of brightness change & 4-segment temporal analysis |\n\nThese features are combined with the star info features **(Rs, Ms, Ts, Mp, e, P, sma, i)** for training the LGBM models.\n\n## LB Score Estimate + Sigma Optimization\n\nFirst, we use **residual errors** from our train process to try predict those pesky **uncertainties (sigma) columns**.\n\nThen we use the competition metric to **estimate our LB Score**.\n\nFinally, we run a **per-wavelength optimization loop to tune sigma values** - and then generate an updated LB estimate.\n* Note: Sigma optimization hasn't produced the gains it has for other notebooks.  Unclear if this is a difference in implementation - or the residual error approach is already optimal.\n\n## Timing\n\n- **Feature extraction**: ~2 hours for training data\n- **Model training**: Few minutes for 1415 models  \n- **Scoring**: ~4 hours total (includes feature extraction on train+test data)","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport lightgbm as lgb\nfrom sklearn.model_selection import KFold\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error, r2_score\nfrom tqdm import tqdm\nimport time\nimport os\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport glob\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom sklearn.preprocessing import StandardScaler, RobustScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.manifold import TSNE\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.neighbors import NearestNeighbors\nfrom scipy.spatial.distance import cdist\nfrom scipy.stats import energy_distance\nfrom scipy import signal, stats\nfrom scipy.stats import skew, kurtosis\nfrom scipy.optimize import minimize_scalar\nfrom scipy.stats import norm\nfrom sklearn.model_selection import GroupKFold\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport sys","metadata":{"trusted":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2025-09-23T11:07:54.441167Z","iopub.execute_input":"2025-09-23T11:07:54.441567Z","iopub.status.idle":"2025-09-23T11:08:01.869092Z","shell.execute_reply.started":"2025-09-23T11:07:54.441523Z","shell.execute_reply":"2025-09-23T11:08:01.868193Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# =========================================================\n# Ariel 2025 — FGS1 특징 (베이직셋) + LightGBM (CPU only)\n#  + Post-Process Patch (튜닝 가능):\n#    - 파장방향 Savitzky–Golay 스무딩 + α 블렌딩\n#    - OOF 잔차 기반 σ(타깃 불확실도) 재보정(+ clip, scale)\n#  + Template 정렬(Nearest Neighbors on PCA of AIRS) + PCA 서브스페이스 회귀 블렌드\n#  + (옵션) 변동성 기반 클러스터링 + GMM 증강\n#  + Wide 튜너로 후보정 하이퍼파라미터 탐색 + OOF GLL 스코어링\n#  + 시각화/캐시/서브미션 유틸\n# =========================================================\n\nimport os, re, glob, time, warnings, random, itertools\nimport numpy as np\nimport pandas as pd\nfrom pathlib import Path\nfrom tqdm import tqdm\nimport matplotlib\nmatplotlib.use(\"Agg\")\nimport matplotlib.pyplot as plt\nfrom scipy.signal import find_peaks, savgol_filter\nfrom scipy.stats import skew, kurtosis, norm\nfrom sklearn.model_selection import KFold\nimport lightgbm as lgb\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.mixture import GaussianMixture\nfrom sklearn.decomposition import PCA\nfrom sklearn.linear_model import Ridge\nfrom sklearn.neighbors import NearestNeighbors  # ★ 누락 보완\nwarnings.filterwarnings(\"ignore\")\ntry:\n    from scipy.fft import dct as _dct\nexcept Exception:\n    from scipy.fftpack import dct as _dct\n\n\n# --- LightGBM 무음 로거 ---\nclass _SilentLogger:\n    def info(self, msg: str):   pass\n    def warning(self, msg: str): pass\ntry:\n    lgb.register_logger(_SilentLogger())\nexcept Exception as e:\n    print(\"[LGBM] register_logger failed:\", e)\n\n# ----------------------- Config -----------------------\nBASE_PATH = \"/kaggle/input/ariel-data-challenge-2025\"\nWORK_DIR  = \"/kaggle/working\"\nPLOT_DIR  = f\"{WORK_DIR}/oof_plots\"\nos.makedirs(WORK_DIR, exist_ok=True); os.makedirs(PLOT_DIR, exist_ok=True)\n\nUSE_CACHED_XY = False\nCACHED_TRAIN_XY = \"/kaggle/input/fg1-train-xy/fgs1_train_xy.csv\"\n\nSEED = 42\nrandom.seed(SEED); np.random.seed(SEED)\n\n# ★ Train만 fraction 조절 (test는 항상 100%)\nTRAIN_FRAC = 1\nMAX_SIGNALS_PER_PLANET = 1\n\n# ★ 밝은 픽셀 파생 (원하면 False로 끄기)\nUSE_BRIGHTPIXEL = True\nBRIGHTPIX_RELAX = 0.75\n\n# ★ Post-process 초기 하이퍼파라미터 (튜닝 전 기본값)\nPOST_SMOOTH_WINDOW = 13     # Savitzky–Golay window (odd)\nPOST_SMOOTH_POLY   = 2\nPOST_BLEND_ALPHA   = 0.65   # 최종 = alpha*raw + (1-alpha)*smooth\nSIGMA_CLIP_FGS     = (0.85, 1.30)\nSIGMA_CLIP_AIRS    = (0.92, 1.22)\nSIGMA_FINAL_SCALE  = 1.04\n\n# LightGBM (CPU only)\nN_SPLITS = 5\nEARLY_STOP = 200\nBASE_LGB_PARAMS = dict(\n    num_leaves=31,\n    max_depth=-1,\n    learning_rate=0.05,\n    n_estimators=3000,\n    subsample=0.9,\n    colsample_bytree=0.9,\n    reg_alpha=1e-2,\n    reg_lambda=1e-1,\n    random_state=SEED,\n    n_jobs=-1,\n    num_threads=-1\n)\n\n# ----------------------- Load -----------------------\ntrain_df     = pd.read_csv(f\"{BASE_PATH}/train.csv\")\ntrain_star   = pd.read_csv(f\"{BASE_PATH}/train_star_info.csv\")\ntest_star    = pd.read_csv(f\"{BASE_PATH}/test_star_info.csv\")\nsample_sub   = pd.read_csv(f\"{BASE_PATH}/sample_submission.csv\")\n\nWL_COLS = sorted([c for c in train_df.columns if re.fullmatch(r\"(wl|wavelength)_\\d+\", c)],\n                 key=lambda x: int(re.findall(r\"\\d+\", x)[0]))\nassert len(WL_COLS) == 283\n\nwldf = pd.read_csv(f\"{BASE_PATH}/wavelengths.csv\")\nwavelengths = wldf[WL_COLS].iloc[0].to_numpy(float)\n\narr = train_df[WL_COLS].to_numpy(dtype=float)\nNAIVE_MEAN  = float(np.nanmean(arr))\nNAIVE_SIGMA = float(np.nanstd(arr) + 1e-12)\n\n# ----------------------- Utils -----------------------\ndef add_interactions_to_X(X, small_feats, large_feats, eps=1e-9):\n    new_cols = []\n    for a, b in itertools.permutations(large_feats, 2):\n        if a in X.columns and b in X.columns:\n            c1 = f\"{a}__minus__{b}\"; X[c1] = X[a] - X[b]; new_cols.append(c1)\n            c2 = f\"{a}__div__{b}\"  ; X[c2] = X[a] / (X[b].abs() + eps); new_cols.append(c2)\n    for a, b in itertools.permutations(small_feats, 2):\n        if a in X.columns and b in X.columns:\n            c3 = f\"{a}__mul__{b}\"  ; X[c3] = X[a] * X[b]; new_cols.append(c3)\n            c4 = f\"{a}__div__{b}\"  ; X[c4] = X[a] / (X[b].abs() + eps); new_cols.append(c4)\n    for a in small_feats:\n        for b in large_feats:\n            if a in X.columns and b in X.columns:\n                c5 = f\"{a}__mul__{b}\"; X[c5] = X[a] * X[b]; new_cols.append(c5)\n    return new_cols\n\ndef get_mean_feature_importance(trained, kind=\"gain\"):\n    feats = trained[\"feature_columns\"]\n    imp = np.zeros(len(feats), dtype=np.float64); n = 0\n    for fold_models in trained[\"models\"]:\n        for m in fold_models:\n            try:\n                arr = m.booster_.feature_importance(importance_type=kind)\n            except Exception:\n                arr = m.feature_importances_\n            imp += np.asarray(arr, dtype=np.float64); n += 1\n    imp = imp / max(n, 1)\n    fi = pd.DataFrame({\"feature\": feats, \"importance\": imp}).sort_values(\n        \"importance\", ascending=False).reset_index(drop=True)\n    return fi\n\ndef build_derived_from_topk(X, fi_df, topk=20, mean_threshold=1.0, eps=1e-9):\n    top = [f for f in fi_df[\"feature\"].tolist() if f in X.columns][:topk]\n    means = X[top].mean(numeric_only=True)\n    small = [f for f in top if means[f] < mean_threshold]\n    large = [f for f in top if means[f] >= mean_threshold]\n    new_cols = add_interactions_to_X(X, small, large, eps=eps)\n    info = {\"small_feats\": small, \"large_feats\": large, \"eps\": eps, \"topk\": top, \"new_cols\": new_cols}\n    print(f\"[derive] top{len(top)} -> small={len(small)}, large={len(large)}, new_cols={len(new_cols)}\")\n    return X, info, new_cols\n\ndef apply_derived_info_to_X(X, info):\n    return add_interactions_to_X(X, info[\"small_feats\"], info[\"large_feats\"], eps=info.get(\"eps\", 1e-9))\n\ndef extract_transit_shape_features(total_flux, baseline_win=5000, local_win=600):\n    f = {}\n    x = np.asarray(total_flux, dtype=np.float64); n = len(x)\n    if n < 10:\n        return {k:0.0 for k in [\n            \"shape_depth\",\"shape_snr\",\"shape_halfdur\",\"shape_asym\",\n            \"shape_curvature\",\"shape_vu_ratio\",\"shape_local_slope_l\",\"shape_local_slope_r\"\n        ]}\n    s = pd.Series(x)\n    base = s.rolling(window=min(baseline_win, max(11, n//10)), center=True).median()\n    base = base.fillna(method=\"bfill\").fillna(method=\"ffill\").to_numpy()\n    d = x - base\n    idx0 = int(np.argmin(d)); depth = float(-d[idx0]); sigma = float(np.nanstd(d))\n    snr = depth / (sigma + 1e-12)\n    w = int(min(local_win, n-1)); lo, hi = max(0, idx0 - w//2), min(n, idx0 + w//2)\n    local = d[lo:hi]\n    if local.size < 5:\n        halfdur = asym = curv = vu = sl_l = sl_r = 0.0\n    else:\n        half = -depth/2.0; rel_idx0 = idx0 - lo\n        left = local[:rel_idx0+1]; right = local[rel_idx0:]\n        Li = np.where(left >= half)[0]; li = int(Li[-1]) if Li.size > 0 else 0\n        Ri = np.where(right >= half)[0]; ri = int(Ri[0]) if Ri.size > 0 else (len(right)-1)\n        halfdur = float((ri + (rel_idx0 - li)))\n        left_w = float(rel_idx0 - li); right_w = float(ri)\n        asym = float((right_w - left_w) / (left_w + right_w + 1e-12))\n        rad = max(5, min(60, local.size//10))\n        a_lo = max(0, rel_idx0 - rad); a_hi = min(local.size, rel_idx0 + rad + 1)\n        yy = local[a_lo:a_hi]; xx = np.arange(len(yy), dtype=np.float64)\n        if len(yy) >= 6:\n            coef = np.polyfit(xx, yy, 2); curv = float(coef[0])\n        else:\n            curv = 0.0\n        d1 = np.diff(local)\n        vu = float(np.median(np.abs(d1)) / (depth + 1e-12))\n        sl_win = max(3, min(40, local.size//20))\n        sl_l = float(np.median(np.diff(local[max(0, rel_idx0-2*sl_win):rel_idx0])))\n        sl_r = float(np.median(np.diff(local[rel_idx0:rel_idx0+2*sl_win])))\n    f.update({\n        \"shape_depth\": depth, \"shape_snr\": snr, \"shape_halfdur\": halfdur, \"shape_asym\": asym,\n        \"shape_curvature\": curv, \"shape_vu_ratio\": vu, \"shape_local_slope_l\": sl_l, \"shape_local_slope_r\": sl_r,\n    })\n    return f\n\ndef extract_diff_volatility_features(total_flux):\n    f = {}\n    x = np.asarray(total_flux, dtype=np.float64)\n    if x.size < 5:\n        return {k:0.0 for k in [\n            \"d1_std\",\"d1_medabs\",\"d1_p95\",\"d1_p99\",\"d2_std\",\n            \"d1_roll100_std_mean\",\"d1_roll1000_std_max\",\"d1_iqr\",\n            \"vol_range_ratio\"\n        ]}\n    d1 = np.diff(x); d2 = np.diff(d1)\n    f[\"d1_std\"] = float(np.std(d1)); f[\"d1_medabs\"] = float(np.median(np.abs(d1)))\n    f[\"d1_p95\"] = float(np.percentile(np.abs(d1), 95)); f[\"d1_p99\"] = float(np.percentile(np.abs(d1), 99))\n    f[\"d2_std\"] = float(np.std(d2)); f[\"d1_iqr\"] = float(np.percentile(d1, 75) - np.percentile(d1, 25))\n    s = pd.Series(d1); r100 = s.rolling(100, center=True).std().dropna(); r1000 = s.rolling(1000, center=True).std().dropna()\n    f[\"d1_roll100_std_mean\"] = float(r100.mean()) if len(r100) else 0.0\n    f[\"d1_roll1000_std_max\"] = float(r1000.max()) if len(r1000) else 0.0\n    rng = float(np.max(x) - np.min(x)); mu = float(np.mean(x))\n    f[\"vol_range_ratio\"] = float(rng / (mu + 1e-12))\n    return f\n\ndef _fft_autocorr(x):\n    n = len(x); x = np.asarray(x, dtype=np.float64); x = x - np.mean(x)\n    nfft = 1 << (n - 1).bit_length()\n    f = np.fft.rfft(x, n=2*nfft)\n    ac = np.fft.irfft(f * np.conj(f))[:n]\n    ac /= (ac[0] + 1e-12)\n    return ac\n\ndef extract_autocorr_rich_features(total_flux, maxlag=4000):\n    f = {}\n    x = np.asarray(total_flux, dtype=np.float64)\n    if x.size < 20:\n        return {k: 0.0 for k in [\n            \"acf_first_zero_lag\",\"acf_peak1_lag\",\"acf_peak1_val\",\n            \"acf_peak2_lag\",\"acf_peak2_val\",\"acf_int_time\"\n        ]}\n    ac = _fft_autocorr(x); m = int(min(maxlag, len(ac)-1)); a = ac[1:m+1]\n    zc = np.where(a <= 0)[0]; f[\"acf_first_zero_lag\"] = float(zc[0]+1) if zc.size else float(m)\n    peaks, _ = find_peaks(a, height=0.05, distance=5)\n    if peaks.size > 0:\n        p_sorted = peaks[np.argsort(a[peaks])[::-1]]\n        p1 = int(p_sorted[0]); f[\"acf_peak1_lag\"] = float(p1+1); f[\"acf_peak1_val\"] = float(a[p1])\n        if len(p_sorted) > 1:\n            p2 = int(p_sorted[1]); f[\"acf_peak2_lag\"] = float(p2+1); f[\"acf_peak2_val\"] = float(a[p2])\n        else:\n            f[\"acf_peak2_lag\"] = 0.0; f[\"acf_peak2_val\"] = 0.0\n    else:\n        f[\"acf_peak1_lag\"] = 0.0; f[\"acf_peak1_val\"] = 0.0; f[\"acf_peak2_lag\"] = 0.0; f[\"acf_peak2_val\"] = 0.0\n    pos = a[a > 0]; f[\"acf_int_time\"] = float(np.sum(pos))\n    return f\n\ndef extract_centroid_track_features(signal_3d, step=50, max_frames=20000):\n    f = {}\n    arr = np.asarray(signal_3d, dtype=np.float32); T, H, W = arr.shape if arr.ndim==3 else (0,0,0)\n    if T == 0:\n        return {k: 0.0 for k in [\n            \"centroid_x_std\",\"centroid_y_std\",\"centroid_drift\",\n            \"centroid_speed_mean\",\"centroid_speed_std\",\"centroid_flux_corr\"\n        ]}\n    idx = np.arange(0, T, step, dtype=int); \n    if idx.size > max_frames: idx = idx[:max_frames]\n    xs = np.arange(W, dtype=np.float32)[None, :, None]\n    ys = np.arange(H, dtype=np.float32)[:, None, None]\n    cx, cy, tf = [], [], []\n    for t in idx:\n        frame = arr[t]; s = float(frame.sum()); tf.append(s)\n        if s <= 0:\n            if cx: cx.append(cx[-1]); cy.append(cy[-1])\n            else:  cx.append(W/2.0);  cy.append(H/2.0)\n            continue\n        cx.append(float((frame * xs.squeeze(2)).sum() / s))\n        cy.append(float((frame * ys.squeeze(2)).sum() / s))\n    cx = np.asarray(cx); cy = np.asarray(cy); tf = np.asarray(tf)\n    if cx.size < 3:\n        return {k: 0.0 for k in [\n            \"centroid_x_std\",\"centroid_y_std\",\"centroid_drift\",\n            \"centroid_speed_mean\",\"centroid_speed_std\",\"centroid_flux_corr\"\n        ]}\n    f[\"centroid_x_std\"] = float(np.std(cx)); f[\"centroid_y_std\"] = float(np.std(cy))\n    dx = cx - cx[0]; dy = cy - cy[0]; f[\"centroid_drift\"] = float(np.hypot(dx[-1], dy[-1]))\n    vx = np.diff(cx); vy = np.diff(cy); speed = np.hypot(vx, vy)\n    f[\"centroid_speed_mean\"] = float(np.mean(speed)); f[\"centroid_speed_std\"]  = float(np.std(speed))\n    if np.std(tf) > 0 and np.std(cx) > 0:\n        f[\"centroid_flux_corr\"] = float(np.corrcoef(tf, cx)[0,1])\n    else:\n        f[\"centroid_flux_corr\"] = 0.0\n    return f\n\ndef extract_segment_features(total_flux, n_segments=4):\n    f = {}\n    x = np.asarray(total_flux, dtype=np.float64); n = len(x)\n    if n < n_segments:\n        return {k: 0.0 for k in [\n            \"seg_mean_max\",\"seg_mean_min\",\"seg_std_mean\",\"seg_min_min\",\n            \"seg_range_max\",\"seg_mean_slope\",\"seg_mean_cv\"\n        ]}\n    segs = np.array_split(x, n_segments)\n    means = np.array([s.mean() for s in segs]); stds  = np.array([s.std()  for s in segs])\n    mins  = np.array([s.min()  for s in segs]); maxs  = np.array([s.max()  for s in segs])\n    rngs  = maxs - mins; xs = np.arange(n_segments, dtype=np.float64)\n    A = np.vstack([xs, np.ones_like(xs)]).T; slope = float(np.linalg.lstsq(A, means, rcond=None)[0][0])\n    f.update({\n        \"seg_mean_max\": float(means.max()), \"seg_mean_min\": float(means.min()),\n        \"seg_std_mean\": float(stds.mean()), \"seg_min_min\": float(mins.min()),\n        \"seg_range_max\": float(rngs.max()), \"seg_mean_slope\": slope,\n        \"seg_mean_cv\": float(stds.mean() / (means.mean() + 1e-12)),\n    })\n    return f\n\ndef _robust_std(x, axis=0):\n    med = np.median(x, axis=axis, keepdims=True)\n    mad = np.median(np.abs(x - med), axis=axis)\n    return 1.4826 * mad + 1e-12\n\ndef load_cached_xy_csv(path, wl_cols=None, index_candidates=(\"planet_id\",\"id\",\"pid\")):\n    df = pd.read_csv(path)\n    if wl_cols is None:\n        y_cols = [c for c in df.columns if re.fullmatch(r\"(wl|wavelength)_\\d+\", c)]\n        if not y_cols: raise ValueError(\"Cached CSV에서 타깃 컬럼(wl_*)을 찾지 못했습니다.\")\n        wl_cols = sorted(y_cols, key=lambda x: int(re.findall(r\"\\d+\", x)[0]))\n    else:\n        y_cols = list(wl_cols)\n    idx_col = next((c for c in index_candidates if c in df.columns), None)\n    if idx_col: df[idx_col] = df[idx_col].astype(str); df = df.set_index(idx_col)\n    else:       df.index = df.index.astype(str)\n    num_cols = df.select_dtypes(include=[np.number]).columns\n    x_cols = [c for c in num_cols if c not in y_cols]\n    X = df[x_cols].astype(np.float32).copy()\n    y = df[y_cols].astype(np.float32).copy()\n    return X, y, wl_cols\n\ndef adaptive_sg_blend(\n    mu_raw,\n    *,\n    window=13,              # mid 구간 window\n    poly=2,                 # SG poly\n    alpha=0.65,             # mid 구간 α\n    w_lo=9,  w_hi=25,       # oscillatory / flat window\n    a_lo=0.85, a_hi=0.45    # oscillatory / flat α\n):\n    \"\"\"\n    스펙트럼 별 roughness = std(diff(mu))로 flat/mid/oscillatory 분기하여\n    서로 다른 (W, α) 적용. 결과는 >=0 로 clip.\n    \"\"\"\n    rough = np.std(np.diff(mu_raw, axis=1), axis=1)\n    q1, q2 = np.percentile(rough, [33, 66])\n    out = np.empty_like(mu_raw)\n    for i, r in enumerate(rough):\n        if r <= q1:            # flat → 크게/많이 스무딩\n            W, A = w_hi, a_hi\n        elif r <= q2:          # mid\n            W, A = window, alpha\n        else:                  # oscillatory → 작게/덜 스무딩\n            W, A = w_lo, a_lo\n        W = max(3, int(W) | 1); P = max(1, min(int(poly), W - 1))\n        sm = savgol_filter(mu_raw[i], window_length=W, polyorder=P, mode=\"interp\")\n        out[i] = A * mu_raw[i] + (1.0 - A) * sm\n    return np.clip(out, 0.0, None)\n\ndef _odd_window(n, prefer):\n    w = int(prefer); w = max(3, min(w, (n if n%2==1 else n-1)))\n    if w % 2 == 0: w -= 1\n    w = max(3, w)\n    return w\n\ndef apply_adc_correction(signal_array, instrument, adc_info_path=f\"{BASE_PATH}/adc_info.csv\"):\n    adc_df = pd.read_csv(adc_info_path)\n    gain = adc_df.at[0, f\"{instrument}_adc_gain\"]\n    offset = adc_df.at[0, f\"{instrument}_adc_offset\"]\n    return signal_array.astype(np.float32) / gain + offset\n\ndef extract_brightpixel_features(signal_3d, total_flux=None, relax=0.30):\n    if total_flux is None:\n        total_flux = np.sum(signal_3d, axis=(1, 2))\n    bidx = int(np.argmax(total_flux))\n    frame = signal_3d[bidx]\n    max_val = float(np.max(frame))\n    if max_val <= 0: mask = (frame == max_val)\n    else:\n        thr = max_val * (1.0 - relax)\n        mask = frame >= thr\n    vals = frame[mask]\n    if vals.size == 0: vals = np.array([0.0], dtype=np.float32)\n    return {\n        \"brightpix_npix\": float(mask.sum()),\n        \"brightpix_brightframe_mean\": float(np.mean(vals)),\n        \"brightpix_brightframe_max\":  float(np.max(vals)),\n        \"brightpix_brightframe_min\":  float(np.min(vals)),\n    }\n\ndef _rolling_median(x, win):\n    import pandas as pd, numpy as np\n    w = int(max(11, min(win, len(x)//2)*1) )\n    if w % 2 == 0: w += 1\n    s = pd.Series(np.asarray(x, float))\n    b = s.rolling(w, center=True).median()\n    return b.fillna(method=\"bfill\").fillna(method=\"ffill\").to_numpy()\n\ndef _detrend_median(x, base_win=5000):\n    x = np.asarray(x, float)\n    base = _rolling_median(x, min(base_win, max(11, len(x)//10)))\n    return x - base\n\ndef _transit_shape_physics(detrended):\n    from scipy.signal import savgol_filter\n    z = savgol_filter(detrended, max(31, (401 if len(detrended)>401 else (len(detrended)//5*2+1))), 3, mode=\"interp\")\n    i0 = int(np.argmin(z)); depth = float(-z[i0])\n    half = -depth/2.0\n    L = np.where(z[:i0] >= half)[0]; R = np.where(z[i0:] >= half)[0]\n    l = int(L[-1]) if L.size else 0\n    r = int(R[0])  if R.size else 0\n    T14 = float(np.sum(z < 0))\n    T23 = float(r + (i0 - l))\n    d1 = np.gradient(z); d2 = np.gradient(d1)\n    ingress = float(-np.min(d1[l:i0])) if i0>l else 0.0\n    egress  = float(np.max(d1[i0:i0+r+1])) if r>0 else 0.0\n    curv    = float(np.median(d2[max(0,i0-50):i0+51]))\n    # 바깥 소음\n    oot_std = float(np.std(np.concatenate([z[:max(1,l)], z[i0+r+1:]]))) + 1e-12\n    return {\n        \"phys_depth\": depth,\n        \"phys_T14\": T14,\n        \"phys_T23\": T23,\n        \"phys_T23_over_T14\": T23 / (T14 + 1e-12),\n        \"phys_ingress_slope\": ingress,\n        \"phys_egress_slope\": egress,\n        \"phys_bottom_curvature\": curv,\n        \"phys_snr_depth\": depth / oot_std,\n        \"phys_energy\": depth * T14\n    }\n\ndef _rednoise_features(x, bin_sizes=(30,60,120,300)):\n    x = np.asarray(x, float)\n    z = x - np.median(x)\n    out = {}\n    base_std = np.std(z) + 1e-12\n    eff_bins = [b for b in bin_sizes if len(z)//b >= 2]\n    for B in eff_bins:\n        n = len(z)//B\n        xb = z[:n*B].reshape(n, B).mean(1)\n        out[f\"rn_bin{B}_std\"]  = float(np.std(xb))\n        out[f\"rn_beta{B}\"]     = float(np.std(xb) / (base_std/np.sqrt(B)))\n    # PSD 기울기(저주파 지배 판단)\n    fx  = np.fft.rfft(z); ps = (np.abs(fx)**2)[1:]\n    fr  = np.fft.rfftfreq(len(z))[1:]\n    m = (fr>0) & np.isfinite(ps)\n    if m.sum()>10:\n        X = np.vstack([np.log10(fr[m]), np.ones(m.sum())]).T\n        a,_ = np.linalg.lstsq(X, np.log10(ps[m] + 1e-20), rcond=None)[0]\n        out[\"rn_psd_slope\"] = float(a)      # -1~-2면 적색 잡음 우세\n    else:\n        out[\"rn_psd_slope\"] = 0.0\n    return out\n\ndef _dct_lowfreq_coeffs(x, m=256, k=16):\n    \"\"\"균일 리샘플 후 DCT 저주파 k개(DC 포함)\"\"\"\n    x = np.asarray(x, float)\n    t = np.linspace(0, 1, len(x))\n    grid = np.linspace(0, 1, m)\n    y = np.interp(grid, t, x)\n    y = (y - y.mean())/(y.std()+1e-12)\n    c = _dct(y, type=2, norm=\"ortho\")[:k]\n    return {f\"dct_lf_{i:02d}\": float(c[i]) for i in range(k)}\n\ndef _wavelet_energy(x, wave=\"db4\", levels=(2,3,4,5)):\n    out = {}\n    try:\n        import pywt\n        x = (np.asarray(x, float) - np.mean(x))/(np.std(x)+1e-12)\n        L = max(levels) if len(x) > 2**max(levels) else int(np.floor(np.log2(len(x))) - 1)\n        L = max(1, min(L, max(levels)))\n        coeffs = pywt.wavedec(x, wave, level=L)\n        for Lv in levels:\n            if Lv <= L:\n                cD = coeffs[Lv]\n                out[f\"wavE_L{Lv}\"] = float(np.mean(cD**2))\n        if L >= 2:\n            lowE  = float(np.mean(coeffs[-1]**2))\n            highE = float(np.mean(np.hstack(coeffs[1:-1])**2) + 1e-12)\n            out[\"wavE_low_over_high\"] = lowE / highE\n        else:\n            out[\"wavE_low_over_high\"] = 0.0\n    except Exception:\n        # PyWavelets 없으면 생략\n        out[\"wavE_low_over_high\"] = 0.0\n    return out\n\ndef _centroid_regression_residuals(signal_3d, step=100):\n    \"\"\"cx, cy, t, t^2 회귀 후 잔차 표준편차/딥\"\"\"\n    arr = np.asarray(signal_3d, np.float32)\n    T,H,W = arr.shape\n    idx = np.arange(0, T, step, dtype=int)\n    xs = np.arange(W, dtype=np.float32)[None, :, None]\n    ys = np.arange(H, dtype=np.float32)[:, None, None]\n    cx, cy, tf = [], [], []\n    for t in idx:\n        f = arr[t]; s = float(f.sum()); tf.append(s)\n        if s <= 0:\n            cx.append(cx[-1] if cx else W/2.0); cy.append(cy[-1] if cy else H/2.0)\n        else:\n            cx.append(float((f*xs.squeeze(2)).sum()/s))\n            cy.append(float((f*ys.squeeze(2)).sum()/s))\n    cx, cy, tf = np.asarray(cx), np.asarray(cy), np.asarray(tf)\n    if cx.size < 5:\n        return {\"cr_res_std\": 0.0, \"cr_res_depth\": 0.0}\n    n = len(tf); t = np.linspace(-1,1,n)\n    A = np.vstack([np.ones(n), cx, cy, t, t**2]).T\n    coef = np.linalg.lstsq(A, tf, rcond=None)[0]\n    res  = tf - A@coef\n    return {\"cr_res_std\": float(np.std(res)), \"cr_res_depth\": float(np.max(np.median(tf)-res))}\n\ndef extract_global_flux_features(total_flux):\n    f = {}\n    f['global_flux_mean'] = np.mean(total_flux); f['global_flux_std']  = np.std(total_flux)\n    f['global_flux_min']  = np.min(total_flux);  f['global_flux_max']  = np.max(total_flux)\n    f['global_flux_range'] = f['global_flux_max'] - f['global_flux_min']\n    f['global_flux_skew']   = skew(total_flux); f['global_flux_kurtosis'] = kurtosis(total_flux)\n    f['global_flux_cv'] = f['global_flux_std'] / (f['global_flux_mean'] + 1e-12)\n    for p in [1,5,10,25,50,75,90,95,99]:\n        f[f'global_flux_p{p}'] = np.percentile(total_flux, p)\n    f['global_flux_depth'] = f['global_flux_mean'] - f['global_flux_min']\n    f['global_flux_depth_ratio'] = f['global_flux_depth'] / (f['global_flux_mean'] + 1e-12)\n    return f\n\ndef extract_rolling_statistics_features(total_flux):\n    f = {}\n    s = pd.Series(total_flux)\n    for window in [50,100,500,1000,2000,5000]:\n        if window < len(total_flux):\n            rm = s.rolling(window, center=True).mean().dropna()\n            rs = s.rolling(window, center=True).std().dropna()\n            rmin = s.rolling(window, center=True).min().dropna()\n            rmax = s.rolling(window, center=True).max().dropna()\n            f[f'rolling{window}_mean_min']  = float(rm.min())\n            f[f'rolling{window}_mean_max']  = float(rm.max())\n            f[f'rolling{window}_mean_std']  = float(rm.std())\n            f[f'rolling{window}_std_mean']  = float(rs.mean())\n            f[f'rolling{window}_std_max']   = float(rs.max())\n            f[f'rolling{window}_deepest_dip'] = float(rmin.min())\n            f[f'rolling{window}_highest_peak'] = float(rmax.max())\n    return f\n\ndef extract_transit_detection_features(total_flux):\n    f = {}\n    baseline = pd.Series(total_flux).rolling(window=5000, center=True).median().fillna(method='bfill').fillna(method='ffill')\n    detrended = total_flux - baseline\n    f['detrended_min'] = float(np.min(detrended))\n    f['detrended_std'] = float(np.std(detrended))\n    f['detrended_skew'] = float(skew(detrended))\n    f['detrended_neg_excursions']  = int(np.sum(detrended < -2*np.std(detrended)))\n    f['detrended_deep_excursions'] = int(np.sum(detrended < -3*np.std(detrended)))\n    thr = np.mean(total_flux) - 1.0*np.std(total_flux)\n    below = total_flux < thr\n    if np.any(below):\n        diff = np.diff(np.concatenate(([False], below, [False])).astype(int))\n        starts = np.where(diff == 1)[0]; ends = np.where(diff == -1)[0]\n        durations = ends - starts\n        f['longest_dip_duration'] = int(np.max(durations)) if len(durations)>0 else 0\n        f['num_dip_periods']      = int(len(durations))\n        f['total_dip_time']       = int(np.sum(durations))\n        f['avg_dip_duration']     = float(np.mean(durations)) if len(durations)>0 else 0.0\n        deepest_idx = int(np.argmin(total_flux))\n        f['deepest_time_fraction'] = float(deepest_idx / len(total_flux))\n    else:\n        f.update(dict(longest_dip_duration=0, num_dip_periods=0, total_dip_time=0, avg_dip_duration=0.0, deepest_time_fraction=0.5))\n    return f\n\ndef extract_frequency_features(total_flux):\n    f = {}\n    x = total_flux - np.mean(total_flux)\n    fft_flux = np.fft.fft(x); fft_power = np.abs(fft_flux)\n    freqs = np.fft.fftfreq(len(x))\n    half = slice(1, len(x)//2)\n    ps = np.abs(fft_power[half]); fr = np.abs(freqs[half])\n    if len(ps)>0:\n        f['fft_peak_power']  = float(np.max(ps))\n        f['fft_total_power'] = float(np.sum(ps))\n        f['fft_mean_power']  = float(np.mean(ps))\n        f['fft_std_power']   = float(np.std(ps))\n        low = np.abs(fr) < 0.01\n        f['fft_low_freq_power']  = float(np.sum(ps[low]))\n        f['fft_low_freq_ratio']  = float(f['fft_low_freq_power']/(f['fft_total_power']+1e-12))\n    else:\n        f.update(dict(fft_peak_power=0, fft_total_power=0, fft_mean_power=0, fft_std_power=0, fft_low_freq_power=0, fft_low_freq_ratio=0))\n    return f\n\ndef extract_gradient_features(total_flux):\n    f = {}\n    d1 = np.diff(total_flux); d2 = np.diff(d1)\n    f['flux_diff1_mean'] = float(np.mean(d1)); f['flux_diff1_std'] = float(np.std(d1))\n    f['flux_diff1_min']  = float(np.min(d1));  f['flux_diff1_max'] = float(np.max(d1))\n    f['flux_diff2_std']  = float(np.std(d2))\n    return f\n\ndef extract_spatial_features(signal_data):\n    f = {}\n    idxs = [0, len(signal_data)//4, len(signal_data)//2, 3*len(signal_data)//4, -1]\n    for i, idx in enumerate(idxs):\n        frame = signal_data[idx]; tot = np.sum(frame)\n        center = frame[12:20,12:20].sum() / (tot+1e-12)\n        f[f'frame{i}_spatial_mean'] = float(np.mean(frame))\n        f[f'frame{i}_spatial_std']  = float(np.std(frame))\n        f[f'frame{i}_concentration'] = float(center)\n    return f\n\ndef extract_sustained_slope_features(\n    total_flux,\n    base_win=5000,\n    sg_win=401,\n    k_mad=5.0,\n    hysteresis=0.6,\n    min_len=80\n):\n    \"\"\"\n    detrend → SG smoothing → |d1| 히스테리시스 run 추출\n    \"\"\"\n    x = np.asarray(total_flux, dtype=np.float64)\n    n = len(x)\n    outs = {k:0.0 for k in [\n        \"ss_n_runs\",\"ss_len_min\",\"ss_len_median\",\"ss_len_mean\",\"ss_len_max\",\n        \"ss_len_sum\",\"ss_len_max_frac\",\"ss_run_mean_absd1_mean\",\"ss_run_peak_absd1_mean\",\n        \"ss_len_at_dip\",\"ss_dist_to_prev_run\",\"ss_dist_to_next_run\",\"ss_n_up\",\"ss_n_down\"\n    ]}\n    if n < 50:\n        return outs\n    # detrend + smoothing\n    s = pd.Series(x)\n    bw = min(base_win, max(11, n//10))\n    base = s.rolling(bw, center=True).median().fillna(method=\"bfill\").fillna(method=\"ffill\").to_numpy()\n    d = x - base\n    W = int(sg_win) | 1\n    if W >= n: W = n-1-(n%2==0)\n    W = max(31, W)\n    y = savgol_filter(d, window_length=W, polyorder=3, mode=\"interp\")\n    # thresholds\n    d1 = np.diff(y)\n    absd1 = np.abs(d1)\n    mad = np.median(np.abs(d1 - np.median(d1))) + 1e-12\n    thr_hi = 1.4826 * mad * k_mad\n    thr_lo = thr_hi * float(hysteresis)\n    mask_hi = absd1 >= thr_hi\n    mask_lo = absd1 >= thr_lo\n    # runs by hysteresis\n    runs = []\n    i, L = 0, len(d1)\n    while i < L:\n        if mask_hi[i]:\n            s_idx = i; i += 1\n            while i < L and mask_lo[i]: i += 1\n            e_idx = i\n            if e_idx - s_idx >= min_len:\n                runs.append((s_idx, e_idx))\n        else:\n            i += 1\n    if not runs:\n        return outs\n    lens = np.array([e - s for s, e in runs], dtype=np.float64)\n    mean_abs = np.array([absd1[s:e].mean() for s, e in runs], dtype=np.float64)\n    peak_abs = np.array([absd1[s:e].max()  for s, e in runs], dtype=np.float64)\n    outs.update({\n        \"ss_n_runs\": float(len(runs)),\n        \"ss_len_min\": float(lens.min()),\n        \"ss_len_median\": float(np.median(lens)),\n        \"ss_len_mean\": float(lens.mean()),\n        \"ss_len_max\": float(lens.max()),\n        \"ss_len_sum\": float(lens.sum()),\n        \"ss_len_max_frac\": float(lens.max() / (n + 1e-12)),\n        \"ss_run_mean_absd1_mean\": float(mean_abs.mean()),\n        \"ss_run_peak_absd1_mean\": float(peak_abs.mean()),\n    })\n    dip = int(np.argmin(y))\n    len_at_dip = 0; prev_end = None; next_start = None\n    for s_idx, e_idx in runs:\n        if s_idx <= dip < e_idx: len_at_dip = e_idx - s_idx\n        if e_idx <= dip: prev_end = e_idx\n        if next_start is None and s_idx > dip: next_start = s_idx\n    outs[\"ss_len_at_dip\"] = float(len_at_dip)\n    outs[\"ss_dist_to_prev_run\"] = float(dip - prev_end) if prev_end is not None else 0.0\n    outs[\"ss_dist_to_next_run\"] = float(next_start - dip) if next_start is not None else 0.0\n    up = sum(np.sign(d1[s:e]).mean() >= 0 for s, e in runs)\n    downs = len(runs) - up\n    outs[\"ss_n_up\"] = float(up); outs[\"ss_n_down\"] = float(downs)\n    return outs\n\ndef extract_enhanced_fgs1_features(signal_3d, verbose=False):\n    \"\"\"\n    기존 베이직셋 + (패치) 트랜짓물리/적색노이즈/DCT/웨이블릿/센트로이드회귀 잔차\n    - 기존 코드의 apply_adc_correction, 기타 유틸이 같은 파일에 있다고 가정\n    \"\"\"\n    # --- ADC 보정 & 총광도 ---\n    signal_3d = apply_adc_correction(signal_3d, instrument='FGS1')\n    total_flux = np.sum(signal_3d, axis=(1,2))\n\n    feats = {}\n\n    # ===== 기존 베이직 특징들 =====\n    feats.update(extract_global_flux_features(total_flux))\n    feats.update(extract_rolling_statistics_features(total_flux))\n    feats.update(extract_transit_detection_features(total_flux))\n    feats.update(extract_frequency_features(total_flux))\n    feats.update(extract_gradient_features(total_flux))\n    feats.update(extract_transit_shape_features(total_flux))\n    feats.update(extract_sustained_slope_features(total_flux, base_win=5000, sg_win=401,\n                                                  k_mad=5.0, hysteresis=0.6, min_len=80))\n    feats.update(extract_diff_volatility_features(total_flux))\n    feats.update(extract_autocorr_rich_features(total_flux))\n    feats.update(extract_centroid_track_features(signal_3d, step=50, max_frames=20000))\n    feats.update(extract_segment_features(total_flux, n_segments=4))\n    feats.update(extract_spatial_features(signal_3d))\n    if USE_BRIGHTPIXEL:\n        feats.update(extract_brightpixel_features(signal_3d, total_flux=total_flux, relax=BRIGHTPIX_RELAX))\n\n    # ====== (패치) 고효과 파생 ======\n    detr = _detrend_median(total_flux, base_win=5000)\n    feats.update(_transit_shape_physics(detr))\n    feats.update(_rednoise_features(total_flux, bin_sizes=(30,60,120,300)))\n    feats.update(_dct_lowfreq_coeffs(detr, m=256, k=16))\n    feats.update(_wavelet_energy(total_flux, wave=\"db4\", levels=(2,3,4,5)))\n    feats.update(_centroid_regression_residuals(signal_3d, step=100))\n\n    return feats\n\ndef filter_similar_derived(X, derived_cols, fi_df=None, thresh=0.90, max_check=2000):\n    cols = [c for c in derived_cols if c in X.columns]\n    const = [c for c in cols if float(np.std(X[c].values)) <= 1e-12]\n    if const: X.drop(columns=const, inplace=True, errors=\"ignore\")\n    cols = [c for c in cols if c not in const]\n    if not cols:\n        print(f\"[derive-dup] const={len(const)} -> no derived left\")\n        return []\n    if fi_df is not None and \"feature\" in fi_df and \"importance\" in fi_df:\n        imp = fi_df.set_index(\"feature\")[\"importance\"]\n        cols = sorted(cols, key=lambda c: -float(imp.get(c, 0.0)))\n    else:\n        cols = sorted(cols, key=lambda c: -float(np.nanvar(X[c].values)))\n    sub = cols[:min(len(cols), max_check)]\n    C = np.abs(pd.DataFrame(X[sub], copy=False).corr().to_numpy())\n    keep, dropped = [], set()\n    for i, ci in enumerate(sub):\n        if ci in dropped: continue\n        keep.append(ci)\n        high = np.where(C[i, i+1:] >= thresh)[0]\n        for h in high:\n            cj = sub[i+1+h]; dropped.add(cj)\n    if dropped: X.drop(columns=list(dropped), inplace=True, errors=\"ignore\")\n    kept_full = [c for c in derived_cols if c in X.columns]\n    print(f\"[derive-dup] checked={len(sub)}, kept={len(keep)}, dropped≈{len(dropped)}, const={len(const)}\")\n    return kept_full\n\ndef load_fgs1_signals_for_planet(planet_id_str, split=\"train\"):\n    base = f\"{BASE_PATH}/{split}/{int(planet_id_str)}\"\n    paths = sorted(glob.glob(f\"{base}/FGS1_signal_*.parquet\"))\n    if MAX_SIGNALS_PER_PLANET is not None:\n        paths = paths[:MAX_SIGNALS_PER_PLANET]\n    return paths\n\ndef build_features(split=\"train\", star_info_df=None, targets_df=None, frac=1.0):\n    sid = star_info_df.copy()\n    sid['planet_id'] = sid['planet_id'].astype(int).astype(str)\n    if split==\"train\" and (0 < frac < 1.0):\n        n = max(1, int(len(sid)*frac))\n        sid = sid.sample(n=n, random_state=SEED).sort_values(\"planet_id\")\n    feats = []\n    for _, row in tqdm(sid.iterrows(), total=len(sid), desc=f\"Extract {split} FGS1 features\"):\n        pid = row['planet_id']\n        for j, p in enumerate(load_fgs1_signals_for_planet(pid, split=split)):\n            try:\n                df = pd.read_parquet(p)\n                signal_3d = df.values.reshape(135000, 32, 32)\n                f = extract_enhanced_fgs1_features(signal_3d, verbose=False)\n                f['planet_id'] = f\"{pid}_{j}\"\n                feats.append(f)\n            except Exception as e:\n                print(f\"[WARN] {split} {pid} signal{j} fail:\", e)\n    X = pd.DataFrame(feats).set_index(\"planet_id\").sort_index()\n    meta = star_info_df.copy(); meta['planet_id'] = meta['planet_id'].astype(int).astype(str)\n    exp = []\n    for pid in X.index:\n        base_id = pid.split(\"_\")[0]\n        m = meta[meta['planet_id']==base_id].copy()\n        if not m.empty:\n            m['planet_id'] = pid; exp.append(m)\n    if len(exp)>0:\n        meta_expand = pd.concat(exp, ignore_index=True).set_index(\"planet_id\").sort_index()\n        X = X.join(meta_expand.select_dtypes(include=[np.number]).astype(np.float32), how=\"left\")\n    X = X.select_dtypes(include=[np.number]).fillna(0).astype(np.float32)\n    if split==\"train\":\n        tgt = targets_df.copy()\n        tgt['planet_id'] = tgt['planet_id'].astype(int).astype(str)\n        tgt = tgt.set_index(\"planet_id\")[WL_COLS].astype(np.float32)\n        y_rows = []\n        for pid in X.index:\n            base_id = pid.split(\"_\")[0]\n            if base_id in tgt.index:\n                r = tgt.loc[base_id].copy(); r.name = pid; y_rows.append(r)\n        y = pd.DataFrame(y_rows).sort_index().astype(np.float32)\n        common = X.index.intersection(y.index)\n        X = X.loc[common]; y = y.loc[common]\n        print(f\"[SHAPE] X:{X.shape}, y:{y.shape}\")\n        return X, y\n    else:\n        print(f\"[SHAPE] X_test:{X.shape}\")\n        return X\n\n# ----------------------- Train / Post-calibration / Predict -----------------------\ndef train_lgbm_cv(X, y):\n    print(\"[LGBM] Device: CPU\")\n    params = BASE_LGB_PARAMS.copy()\n    models = [[] for _ in range(y.shape[1])]\n    oof_pred = np.zeros_like(y.values, dtype=np.float32)\n    kf = KFold(n_splits=N_SPLITS, shuffle=True, random_state=SEED)\n    for t_idx, target in enumerate(WL_COLS):\n        y_t = y[target].values\n        fold_models = []\n        for fold, (tr, va) in enumerate(kf.split(X), 1):\n            m = lgb.LGBMRegressor(**params)\n            m.fit(\n                X.iloc[tr], y_t[tr],\n                eval_set=[(X.iloc[va], y_t[va])],\n                eval_metric=\"l2\",\n                callbacks=[lgb.early_stopping(EARLY_STOP, verbose=False)]\n            )\n            pred_va = m.predict(X.iloc[va]).astype(np.float32)\n            oof_pred[va, t_idx] = pred_va\n            fold_models.append(m)\n        models[t_idx] = fold_models\n        if (t_idx+1) % 25 == 0:\n            print(f\"[{t_idx+1}/{y.shape[1]}] targets trained\")\n    resid_base = y.values - oof_pred\n    mad   = np.median(np.abs(resid_base - np.median(resid_base, axis=0, keepdims=True)), axis=0)\n    target_unc = 1.4826 * mad\n    return {\n        \"models\": models,\n        \"feature_columns\": list(X.columns),\n        \"target_columns\": WL_COLS,\n        \"target_uncertainties\": target_unc.astype(np.float32),\n        \"oof_predictions\": oof_pred,\n        \"y_true\": y.values,\n        \"index\": list(y.index),\n    }\n\ndef _apply_smooth_blend(mu, window=13, poly=2, alpha=0.65):\n    W = _odd_window(mu.shape[1], window)\n    smooth = savgol_filter(mu, window_length=W, polyorder=poly, axis=1, mode=\"interp\")\n    mu2 = alpha * mu + (1.0 - alpha) * smooth\n    mu2 = np.clip(mu2, 0.0, None)\n    return mu2\n\ndef fit_postcalibration(trained, y_true,\n                        smooth_window=POST_SMOOTH_WINDOW,\n                        smooth_poly=POST_SMOOTH_POLY,\n                        blend_alpha=POST_BLEND_ALPHA,\n                        clip_fgs=SIGMA_CLIP_FGS, clip_airs=SIGMA_CLIP_AIRS,\n                        final_sigma_scale=SIGMA_FINAL_SCALE):\n    oof = trained[\"oof_predictions\"].astype(np.float64)\n    base_unc = trained[\"target_uncertainties\"].astype(np.float64)\n    bias = np.mean(y_true - oof, axis=0)  # y - pred\n    oof_corr = oof + bias\n    oof_final = _apply_smooth_blend(oof_corr, window=smooth_window, poly=smooth_poly, alpha=blend_alpha)\n    resid = y_true - oof_final\n    mad = np.median(np.abs(resid - np.median(resid, axis=0, keepdims=True)), axis=0)\n    robust_std = 1.4826 * mad + 1e-12\n    scale = robust_std / (base_unc + 1e-12)\n    if scale.size >= 1: scale[0]  = np.clip(scale[0],  clip_fgs[0],  clip_fgs[1])\n    if scale.size >= 2: scale[1:] = np.clip(scale[1:], clip_airs[0], clip_airs[1])\n    sigma_calib = base_unc * scale * final_sigma_scale\n    trained[\"post_bias\"] = bias.astype(np.float32)\n    trained[\"post_alpha\"] = float(blend_alpha)\n    trained[\"post_smooth_window\"] = int(smooth_window)\n    trained[\"post_smooth_poly\"] = int(smooth_poly)\n    trained[\"calibrated_uncertainties\"] = sigma_calib.astype(np.float32)\n    trained[\"oof_predictions_post\"] = oof_final.astype(np.float32)\n    return trained\n\ndef apply_postprocess_to_preds(mu_raw, trained):\n    bias = trained.get(\"post_bias\", np.zeros(mu_raw.shape[1], dtype=np.float32))\n    alpha = trained.get(\"post_alpha\", POST_BLEND_ALPHA)\n    window = trained.get(\"post_smooth_window\", POST_SMOOTH_WINDOW)\n    poly   = trained.get(\"post_smooth_poly\", POST_SMOOTH_POLY)\n    mu_corr = mu_raw + bias\n    mu_post = _apply_smooth_blend(mu_corr, window=window, poly=poly, alpha=alpha)\n    return mu_post\n\ndef gll_score_numpy(y_true, y_pred, sigma_pred,\n                    naive_mean, naive_sigma,\n                    fsg_sigma_true=1e-6, airs_sigma_true=1e-5, fgs_weight=0.4):\n    sigma_pred = np.clip(sigma_pred, 1e-15, None)\n    n_samples, n_waves = sigma_pred.shape\n    sigma_true = np.append([fsg_sigma_true], np.full(n_waves-1, airs_sigma_true))\n    sigma_true = np.tile(sigma_true, (n_samples, 1))\n    weights = np.append([fgs_weight], np.ones(n_waves-1))\n    weights = np.tile(weights, (n_samples, 1))\n    gll_pred  = norm.logpdf(y_true, loc=y_pred, scale=sigma_pred)\n    gll_true  = norm.logpdf(y_true, loc=y_true, scale=sigma_true)\n    gll_naive = norm.logpdf(y_true, loc=naive_mean, scale=naive_sigma)\n    ind = (gll_pred - gll_naive) / (gll_true - gll_naive + 1e-12)\n    return float(np.clip(np.average(ind, weights=weights), 0.0, 1.0))\n\ndef train_pca_coefficient_regressor(X_train: pd.DataFrame, Y_train: np.ndarray,\n                                    k=24, alpha=1.0, seed=42):\n    y_mean = Y_train.mean(0)\n    y_std  = Y_train.std(0) + 1e-12\n    Y_std = (Y_train - y_mean) / y_std\n    k_eff = int(min(k, max(4, Y_std.shape[1]//4)))\n    pca = PCA(n_components=k_eff, random_state=seed)\n    C = pca.fit_transform(Y_std)\n    reg = Ridge(alpha=alpha)\n    reg.fit(X_train.values, C)\n    return {\"pca\": pca, \"reg\": reg, \"y_mean\": y_mean.astype(np.float32), \"y_std\": y_std.astype(np.float32)}\n\ndef predict_with_pca_regressor(pack, X: pd.DataFrame):\n    C_hat = pack[\"reg\"].predict(X.values)\n    Y_std_hat = pack[\"pca\"].inverse_transform(C_hat)\n    Y_hat = Y_std_hat * pack[\"y_std\"] + pack[\"y_mean\"]\n    return np.clip(Y_hat, 0.0, None).astype(np.float32)\n\ndef _xcorr_shift(a, b, max_shift=3):\n    best_d, best_corr = 0, -1e18\n    for d in range(-max_shift, max_shift+1):\n        if d < 0:   x, y = a[:len(a)+d], b[-d:]\n        elif d > 0: x, y = a[d:], b[:len(b)-d]\n        else:       x, y = a, b\n        if len(x) < 10: continue\n        c = np.corrcoef(x, y)[0, 1]\n        if np.isfinite(c) and c > best_corr:\n            best_corr, best_d = c, d\n    return best_d\n\ndef _shift_apply_air(mu_row, d):\n    out = mu_row.copy()\n    if d == 0: return out\n    a = out[1:]\n    if d > 0:\n        out[1+d:] = a[:-d]; out[1:1+d] = a[:1]\n    else:\n        d = -d\n        out[1:-d] = a[d:]; out[-d:] = a[-1]\n    return out\n\ndef build_template_library(Y_train, k_pca=12, n_neighbors=8, seed=42):\n    pca = PCA(n_components=k_pca, random_state=seed)\n    Z = pca.fit_transform(Y_train[:, 1:])  # AIRS만 사용\n    nn = NearestNeighbors(n_neighbors=n_neighbors, metric='euclidean')\n    nn.fit(Z)\n    return {\"pca\": pca, \"nn\": nn, \"Y\": Y_train.astype(np.float32)}\n\ndef align_with_templates(mu_pred_all, lib, max_shift=3):\n    pca, nn, Y = lib[\"pca\"], lib[\"nn\"], lib[\"Y\"]\n    Zp = pca.transform(mu_pred_all[:, 1:])\n    _, idx = nn.kneighbors(Zp, return_distance=True)\n    out = np.empty_like(mu_pred_all)\n    for i in range(mu_pred_all.shape[0]):\n        tpl = Y[idx[i]].mean(axis=0)\n        d = _xcorr_shift(mu_pred_all[i, 1:], tpl[1:], max_shift=max_shift)\n        out[i] = _shift_apply_air(mu_pred_all[i], d)\n    return out\n\ndef predict_with_uncertainty(trained, X):\n    models = trained[\"models\"]; feats = trained[\"feature_columns\"]\n    X = X.reindex(columns=feats, fill_value=0)\n    n_targets = len(models)\n    preds = np.zeros((len(X), n_targets), dtype=np.float32)\n    for t in range(n_targets):\n        fold_models = models[t]\n        fold_preds = [m.predict(X).astype(np.float32) for m in fold_models]\n        preds[:, t] = np.mean(fold_preds, axis=0)\n    # 후보정\n    mu = apply_postprocess_to_preds(preds, trained)\n    # 템플릿 정렬\n    if \"template_lib\" in trained:\n        mu = align_with_templates(mu, trained[\"template_lib\"], max_shift=trained.get(\"max_shift\", 3))\n    # PCA 서브스페이스 블렌드\n    if \"pca_pack\" in trained:\n        mu_pca = predict_with_pca_regressor(trained[\"pca_pack\"], X)\n        rough = np.std(np.diff(mu[:, 1:], axis=1), axis=1, keepdims=True)  # AIRS roughness\n        q1, q2 = np.quantile(rough, [0.2, 0.8])\n        w = np.clip((rough - q1) / (q2 - q1 + 1e-12), 0, 1) * (0.6 - 0.3) + 0.3  # 0.3~0.6\n        mu = (1 - w) * mu + w * mu_pca\n    # SG 한번 더(미세 계단 제거)\n    mu = adaptive_sg_blend(mu, window=POST_SMOOTH_WINDOW, poly=POST_SMOOTH_POLY, alpha=POST_BLEND_ALPHA)\n    # σ\n    sigmas_base = trained.get(\"calibrated_uncertainties\", trained[\"target_uncertainties\"])\n    sigmas = np.tile(sigmas_base, (len(X), 1))\n    return mu.astype(np.float32), sigmas.astype(np.float32)\n\ndef make_submission(trained, test_star_info, out_path=\"submission.csv\"):\n    X_test = build_features(split=\"test\", star_info_df=test_star_info, targets_df=None, frac=1.0)\n    if \"derived_info\" in trained:\n        apply_derived_info_to_X(X_test, trained[\"derived_info\"])\n    y_pred, sigma_pred = predict_with_uncertainty(trained, X_test)\n    base_ids = [idx.split(\"_\")[0] for idx in X_test.index]\n    df_mu    = pd.DataFrame(y_pred, index=base_ids, columns=WL_COLS)\n    df_sigma = pd.DataFrame(sigma_pred, index=base_ids, columns=[f\"sigma_{i+1}\" for i in range(y_pred.shape[1])])\n    mu_mean    = df_mu.groupby(df_mu.index).mean()\n    sigma_mean = df_sigma.groupby(df_sigma.index).mean()\n    sub = pd.concat([mu_mean, sigma_mean], axis=1).reset_index().rename(columns={\"index\":\"planet_id\"})\n    sub[\"planet_id\"] = sub[\"planet_id\"].astype(int)\n    cols = list(sample_sub.columns)\n    sub = sub.set_index(\"planet_id\").reindex(sample_sub[\"planet_id\"]).reset_index()\n    sub = sub[[\"planet_id\"] + cols[1:]].fillna(0)\n    sub.to_csv(out_path, index=False)\n    print(f\"[SAVE] submission -> {out_path}  shape={sub.shape}\")\n    return sub\n\n#--------------------cluster, GMM--------------------\ndef _select_volatility_columns(X: pd.DataFrame):\n    cols = []\n    prefixes = (\n        \"d1_\", \"d2_\", \"d3_\", \"diff_\", \"delta_\", \"flux_diff\",\n        \"vol_\", \"std_\", \"var_\", \"mad_\", \"iqr_\", \"rms_\", \"zscore_\",\n        \"acf_\", \"pacf_\", \"autocorr_\", \"lag_\", \"xcorr_\",\n        \"fft_\", \"spec_\", \"psd_\", \"welch_\", \"stft_\", \"hilbert_\", \"envelope_\", \"band_\",\n        \"seg_\", \"segment_\", \"cpd_\", \"changepoint_\", \"ruptures_\", \"bkpt_\",\n        \"rolling\", \"roll_\", \"window_\", \"mov_\", \"ma_\", \"sma_\", \"ema_\", \"ewm_\",\n        \"detrended_\", \"resid_\", \"noise_\", \"snr_\",\n        \"shape_\", \"curv_\", \"curvature_\", \"asym_\", \"half\", \"fwhm_\",\n        \"ss_\", \"centroid_\",\n    )\n    extra_exact = {\n        \"global_flux_cv\", \"shape_vu_ratio\", \"shape_halfdur\", \"shape_asym\", \"shape_curvature\",\n        \"d1_roll100_std_mean\", \"d1_roll1000_std_max\", \"d1_iqr\",\n        \"global_flux_std\", \"global_flux_range\", \"global_flux_depth\", \"global_flux_depth_ratio\",\n        \"global_flux_skew\", \"global_flux_kurtosis\", \"vol_range_ratio\",\n        \"acf_int_time\", \"acf_peak1_val\", \"acf_peak2_val\",\n        \"fft_peak_power\", \"fft_total_power\", \"fft_low_freq_ratio\",\n        \"ss_n_runs\", \"ss_len_max\", \"ss_len_sum\", \"ss_len_max_frac\",\n        \"ss_run_mean_absd1_mean\", \"ss_run_peak_absd1_mean\",\n    }\n    for c in X.columns:\n        if c in extra_exact: cols.append(c); continue\n        if any(c.startswith(p) for p in prefixes): cols.append(c)\n    cols = [c for c in cols if c in X.columns]\n    if len(cols) < 5: cols = X.columns.tolist()\n    return sorted(list(dict.fromkeys(cols)))\n\ndef cluster_and_gmm_augment(\n    X: pd.DataFrame,\n    y: pd.DataFrame,\n    *,\n    k_clusters: int = 3,\n    aug_mult: float = 0.5,\n    gm_components: int = 3,\n    y_pca_dim: int = 24,\n    random_state: int = 42,\n    min_cluster_size: int = 30\n):\n    rng = np.random.RandomState(random_state)\n    X = X.copy(); y = y.copy()\n    vol_cols = _select_volatility_columns(X)\n    X_vol = X[vol_cols].values.astype(np.float32)\n    vol_scaler = StandardScaler(); X_vol_scaled = vol_scaler.fit_transform(X_vol)\n    y_scaler = StandardScaler(with_mean=True, with_std=True)\n    y_std = y_scaler.fit_transform(y.values.astype(np.float32))\n    y_pca_dim_eff = int(min(y_pca_dim, y_std.shape[1]-1, max(4, y_std.shape[1]//4)))\n    pca = PCA(n_components=y_pca_dim_eff, random_state=random_state)\n    y_pca = pca.fit_transform(y_std)\n    y_pca_scaler = StandardScaler(); y_pca_scaled = y_pca_scaler.fit_transform(y_pca)\n    gmm_cluster = GaussianMixture(\n        n_components=k_clusters, covariance_type=\"full\",\n        random_state=random_state, reg_covar=1e-6\n    )\n    gmm_cluster.fit(X_vol_scaled)\n    probs = gmm_cluster.predict_proba(X_vol_scaled)\n    labels = probs.argmax(axis=1)\n    X_synth_rows, y_synth_rows, synth_index = [], [], []\n    all_cols = X.columns.tolist()\n    other_cols = [c for c in all_cols if c not in vol_cols]\n    X_df = X.copy()\n    for c in range(k_clusters):\n        mask = (labels == c); n_c = int(mask.sum())\n        if n_c < min_cluster_size or aug_mult <= 0: continue\n        n_aug = int(max(1, round(n_c * aug_mult)))\n        Z = np.hstack([X_vol_scaled[mask], y_pca_scaled[mask]])\n        cov_type = \"diag\" if n_c < 200 else \"full\"\n        gm = GaussianMixture(\n            n_components=min(gm_components, max(1, n_c // 20)),\n            covariance_type=cov_type, random_state=random_state, reg_covar=1e-6\n        )\n        gm.fit(Z)\n        Z_synth, _ = gm.sample(n_aug)\n        x_vol_synth_scaled = Z_synth[:, :X_vol_scaled.shape[1]]\n        y_pca_synth_scaled = Z_synth[:, X_vol_scaled.shape[1]:]\n        x_vol_synth = vol_scaler.inverse_transform(x_vol_synth_scaled)\n        y_pca_synth = y_pca_scaler.inverse_transform(y_pca_synth_scaled)\n        y_std_synth = pca.inverse_transform(y_pca_synth)\n        y_synth = y_scaler.inverse_transform(y_std_synth)\n        y_synth = np.clip(y_synth, 0.0, None)\n        cluster_stats = X_df[mask].agg([\"mean\", \"std\"])\n        other_means = cluster_stats.loc[\"mean\", other_cols].values.astype(np.float32)\n        other_stds  = cluster_stats.loc[\"std\",  other_cols].values.astype(np.float32)\n        jitter = rng.normal(loc=0.0, scale=np.maximum(other_stds*0.05, 1e-6), size=(n_aug, len(other_cols)))\n        for i in range(n_aug):\n            row = np.empty(len(all_cols), dtype=np.float32)\n            row[[all_cols.index(c) for c in vol_cols]] = x_vol_synth[i]\n            row[[all_cols.index(c) for c in other_cols]] = (other_means + jitter[i])\n            X_synth_rows.append(row)\n        y_synth_rows.append(y_synth)\n        synth_index += [f\"synth_c{c}_{i}\" for i in range(n_aug)]\n    if not X_synth_rows:\n        info = {\n            \"labels\": labels, \"cluster_sizes\": [int((labels==i).sum()) for i in range(k_clusters)],\n            \"X_aug\": X, \"y_aug\": y, \"vol_cols\": vol_cols\n        }\n        print(\"[AUG] No augmentation performed (clusters too small or aug_mult<=0).\")\n        return labels, gmm_cluster, info\n    X_synth = pd.DataFrame(X_synth_rows, columns=all_cols, index=synth_index).astype(np.float32)\n    y_synth = pd.DataFrame(np.vstack(y_synth_rows), columns=y.columns, index=synth_index).astype(np.float32)\n    X_aug = pd.concat([X, X_synth], axis=0); y_aug = pd.concat([y, y_synth], axis=0)\n    info = {\n        \"labels\": labels, \"cluster_sizes\": [int((labels==i).sum()) for i in range(k_clusters)],\n        \"n_synth\": len(X_synth), \"X_aug\": X_aug, \"y_aug\": y_aug, \"vol_cols\": vol_cols\n    }\n    print(f\"[AUG] clusters={info['cluster_sizes']}, synth={info['n_synth']} rows\")\n    return labels, gmm_cluster, info\n\n# ----------------------- Viz -----------------------\ndef plot_oof_lines(y_true, y_pred, wavelengths, ids, out_dir, batch_size=10, show=False):\n    Path(out_dir).mkdir(parents=True, exist_ok=True)\n    N = len(ids); batches = (N + batch_size - 1)//batch_size\n    for b in range(batches):\n        idxs = slice(b*batch_size, min((b+1)*batch_size, N))\n        n = idxs.stop - idxs.start; ncols=5; nrows=int(np.ceil(n/ncols))\n        fig, axes = plt.subplots(nrows, ncols, figsize=(18, 2.8*nrows), sharex=True)\n        axes = np.atleast_1d(axes).ravel()\n        for i, k in enumerate(range(idxs.start, idxs.stop)):\n            ax = axes[i]\n            yt = y_true[k]; yp = y_pred[k]\n            ax.plot(wavelengths, yt, lw=1.5, label=\"truth\")\n            ax.plot(wavelengths, yp, lw=1.2, alpha=0.9, label=\"pred\")\n            ax.set_title(f\"{ids[k]}\", fontsize=9); ax.grid(True, alpha=0.3)\n            if i//ncols == nrows-1: ax.set_xlabel(\"Wavelength (µm)\")\n            if i % ncols == 0:      ax.set_ylabel(\"Transit depth\")\n        for j in range(n, nrows*ncols): axes[j].set_visible(False)\n        if b==0: axes[0].legend(loc=\"best\", fontsize=8)\n        fig.suptitle(f\"OOF spectra (batch {b+1}/{batches})\", y=0.98, fontsize=12)\n        fig.tight_layout()\n        out = Path(out_dir)/f\"oof_batch_{b+1:03d}.png\"\n        fig.savefig(out, dpi=150)\n        if show: plt.show()\n        plt.close(fig)\n\ndef plot_bias_sigma(trained, wavelengths, out_dir):\n    Path(out_dir).mkdir(parents=True, exist_ok=True)\n    bias = trained.get(\"post_bias\", np.zeros_like(wavelengths, dtype=np.float32))\n    sig_base = trained[\"target_uncertainties\"]\n    sig_cal  = trained.get(\"calibrated_uncertainties\", sig_base)\n    scale = sig_cal / (sig_base + 1e-12)\n    fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)\n    ax = axes[0]\n    ax.plot(wavelengths, bias, lw=1.2)\n    ax.axhline(0, ls=\"--\", c=\"k\", lw=0.8)\n    ax.set_ylabel(\"bias (y - oof)\")\n    ax.set_title(\"λ-wise bias / uncertainty / scale\")\n    ax.grid(True, alpha=0.3)\n    ax = axes[1]\n    ax.plot(wavelengths, sig_base, lw=1.0, label=\"base_unc\")\n    ax.plot(wavelengths, sig_cal,  lw=1.0, label=\"calibrated_unc\")\n    ax.set_ylabel(\"σ\"); ax.legend(); ax.grid(True, alpha=0.3)\n    ax = axes[2]\n    ax.plot(wavelengths, scale, lw=1.0)\n    ax.axhline(1.0, ls=\"--\", c=\"k\", lw=0.8)\n    ax.set_ylabel(\"scale (=σ_cal/σ_base)\"); ax.set_xlabel(\"Wavelength (µm)\")\n    ax.grid(True, alpha=0.3)\n    fig.tight_layout()\n    out = Path(out_dir)/\"postcalib_bias_sigma.png\"\n    fig.savefig(out, dpi=150); plt.close(fig)\n\ndef plot_residual_std_vs_sigma(y_true, trained, wavelengths, out_dir):\n    Path(out_dir).mkdir(parents=True, exist_ok=True)\n    mu = trained.get(\"oof_predictions_post\", trained[\"oof_predictions\"])\n    resid = y_true - mu\n    med = np.median(resid, axis=0, keepdims=True)\n    mad = np.median(np.abs(resid - med), axis=0)\n    robstd = 1.4826 * mad\n    sig_cal = trained.get(\"calibrated_uncertainties\", trained[\"target_uncertainties\"])\n    fig = plt.figure(figsize=(12,4.2)); ax = plt.gca()\n    ax.plot(wavelengths, robstd, lw=1.1, label=\"residual robust std\")\n    ax.plot(wavelengths, sig_cal, lw=1.1, label=\"pred σ (calibrated)\")\n    ax.set_title(\"Residual robust std vs predicted σ\"); ax.set_xlabel(\"Wavelength (µm)\"); ax.set_ylabel(\"amplitude\")\n    ax.grid(True, alpha=0.3); ax.legend()\n    fig.tight_layout(); out = Path(out_dir)/\"resid_vs_sigma.png\"\n    fig.savefig(out, dpi=150); plt.close(fig)\n\ndef plot_scatter_and_residuals(y_true, y_pred, wavelengths, out_dir, picks=(0,50,150,250), show=False):\n    Path(out_dir).mkdir(parents=True, exist_ok=True)\n    for j in picks:\n        fig = plt.figure(figsize=(12,4))\n        ax1 = plt.subplot(1,2,1)\n        ax1.scatter(y_true[:,j], y_pred[:,j], s=6, alpha=0.6)\n        lims = [min(ax1.get_xlim()[0], ax1.get_ylim()[0]), max(ax1.get_xlim()[1], ax1.get_ylim()[1])]\n        ax1.plot(lims, lims, 'k--', lw=1); ax1.set_xlim(lims); ax1.set_ylim(lims)\n        ax1.set_title(f\"y_true vs pred @λ{j+1} ({wavelengths[j]:.3f} µm)\")\n        ax1.set_xlabel(\"y_true\"); ax1.set_ylabel(\"y_pred\"); ax1.grid(True, alpha=0.3)\n        ax2 = plt.subplot(1,2,2)\n        res = y_true[:,j] - y_pred[:,j]\n        ax2.hist(res, bins=40, alpha=0.8)\n        ax2.set_title(f\"Residuals @λ{j+1}  mean={res.mean():.3e}, std={res.std():.3e}\")\n        ax2.grid(True, alpha=0.3)\n        fig.tight_layout()\n        out = Path(out_dir)/f\"oof_scatter_residuals_lambda_{j+1:03d}.png\"\n        fig.savefig(out, dpi=150)\n        if show: plt.show()\n        plt.close(fig)\n\ndef save_cached_xy_csv(X: pd.DataFrame, y: pd.DataFrame, path: str, wl_cols=None, index_name=\"planet_id\"):\n    dfX = X.copy(); dfY = y.copy() if y is not None else None\n    idx = dfX.index.astype(str); dfX.index = idx\n    if dfY is not None:\n        dfY = dfY.copy(); dfY.index = dfY.index.astype(str)\n        if wl_cols is not None: dfY = dfY[wl_cols]\n        common = dfX.index.intersection(dfY.index); dfX = dfX.loc[common]; dfY = dfY.loc[common]\n        out = pd.concat([dfX, dfY], axis=1)\n    else:\n        out = dfX\n    out.to_csv(path, index=True)\n    print(f\"[CACHE] saved X,y -> {path}  shape={out.shape}\")\n\n# ----------------------- Postprocess 튜닝 -----------------------\ndef _cv_eval_postprocess(\n    oof_raw, y_true,\n    smooth_window, smooth_poly, alpha,\n    clip_fgs, clip_airs, final_sigma_scale,\n    naive_mean, naive_sigma, fgs_weight=0.4,\n    n_splits=5, seed=42\n):\n    kf = KFold(n_splits=n_splits, shuffle=True, random_state=seed)\n    scores = []\n    for tr, va in kf.split(oof_raw):\n        oof_tr, y_tr = oof_raw[tr], y_true[tr]\n        oof_va, y_va = oof_raw[va], y_true[va]\n        bias = np.mean(y_tr - oof_tr, axis=0)\n        resid_tr_base = y_tr - oof_tr\n        base_unc = _robust_std(resid_tr_base, axis=0)\n        oof_tr_post = _apply_bias_smooth_blend(oof_tr, bias, smooth_window, smooth_poly, alpha)\n        resid_tr_post = y_tr - oof_tr_post\n        post_unc = _robust_std(resid_tr_post, axis=0)\n        scale = post_unc / (base_unc + 1e-12)\n        if scale.size >= 1: scale[0]  = np.clip(scale[0],  clip_fgs[0],  clip_fgs[1])\n        if scale.size >= 2: scale[1:] = np.clip(scale[1:], clip_airs[0], clip_airs[1])\n        sigma_vec = base_unc * scale * final_sigma_scale\n        oof_va_post = _apply_bias_smooth_blend(oof_va, bias, smooth_window, smooth_poly, alpha)\n        sigma_val = np.tile(sigma_vec, (len(va), 1))\n        sc = gll_score_numpy(y_va, oof_va_post, sigma_val, naive_mean, naive_sigma, fgs_weight=fgs_weight)\n        scores.append(sc)\n    return float(np.mean(scores)), float(np.std(scores))\n\ndef _apply_bias_smooth_blend(mu, bias, window, poly, alpha):\n    mu_corr = mu + bias\n    W = _odd_window(mu_corr.shape[1], window)\n    poly = int(min(int(poly), W - 1))\n    sm  = savgol_filter(mu_corr, window_length=W, polyorder=poly, axis=1, mode=\"interp\")\n    out = alpha * mu_corr + (1.0 - alpha) * sm\n    return np.clip(out, 0.0, None)\n\ndef tune_postprocess_hparams(\n    trained, y_true, naive_mean, naive_sigma,\n    init=dict(window=13, poly=2, alpha=0.65,\n              clip_fgs=(0.85,1.30), clip_airs=(0.92,1.22), sigma_scale=1.04),\n    n_splits=5, seed=42, verbose=True\n):\n    oof_raw = trained[\"oof_predictions\"].astype(np.float32)\n    win_grid  = [7, 9, 11, 13, 15, 17]\n    poly_grid = [2, 3]\n    alp_grid  = [0.50, 0.60, 0.65, 0.70, 0.75]\n    best = dict(**init); best_score = -1.0\n    for w in win_grid:\n        for p in poly_grid:\n            if p >= w:  continue\n            for a in alp_grid:\n                sc, sd = _cv_eval_postprocess(\n                    oof_raw, y_true, w, p, a,\n                    best[\"clip_fgs\"], best[\"clip_airs\"], best[\"sigma_scale\"],\n                    naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n                )\n                if sc > best_score:\n                    best_score = sc; best.update(window=w, poly=p, alpha=a)\n                    if verbose:\n                        print(f\"[Stage1] win={w}, poly={p}, alpha={a} -> GLL={sc:.6f} (±{sd:.6f})\")\n    fgs_lo_grid = [0.80, 0.85, 0.90]; fgs_hi_grid = [1.20, 1.25, 1.30, 1.35]\n    for lo in fgs_lo_grid:\n        for hi in fgs_hi_grid:\n            if lo >= hi: continue\n            sc, sd = _cv_eval_postprocess(\n                oof_raw, y_true, best[\"window\"], best[\"poly\"], best[\"alpha\"],\n                (lo, hi), best[\"clip_airs\"], best[\"sigma_scale\"],\n                naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n            )\n            if sc > best_score:\n                best_score = sc; best[\"clip_fgs\"] = (lo, hi)\n                if verbose: print(f\"[Stage2a] FGS clip={best['clip_fgs']} -> GLL={sc:.6f} (±{sd:.6f})\")\n    air_lo_grid = [0.90, 0.92, 0.94]; air_hi_grid = [1.18, 1.20, 1.22, 1.24]\n    for lo in air_lo_grid:\n        for hi in air_hi_grid:\n            if lo >= hi: continue\n            sc, sd = _cv_eval_postprocess(\n                oof_raw, y_true, best[\"window\"], best[\"poly\"], best[\"alpha\"],\n                best[\"clip_fgs\"], (lo, hi), best[\"sigma_scale\"],\n                naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n            )\n            if sc > best_score:\n                best_score = sc; best[\"clip_airs\"] = (lo, hi)\n                if verbose: print(f\"[Stage2b] AIRS clip={best['clip_airs']} -> GLL={sc:.6f} (±{sd:.6f})\")\n    for s in [1.00, 1.02, 1.04, 1.06]:\n        sc, sd = _cv_eval_postprocess(\n            oof_raw, y_true, best[\"window\"], best[\"poly\"], best[\"alpha\"],\n            best[\"clip_fgs\"], best[\"clip_airs\"], s,\n            naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n        )\n        if sc > best_score:\n            best_score = sc; best[\"sigma_scale\"] = s\n            if verbose: print(f\"[Stage2c] sigma_scale={s} -> GLL={sc:.6f} (±{sd:.6f})\")\n    if verbose:\n        print(\"\\n[BEST] \",\n              f\"win={best['window']}, poly={best['poly']}, alpha={best['alpha']}, \",\n              f\"FGS={best['clip_fgs']}, AIRS={best['clip_airs']}, scale={best['sigma_scale']}, \",\n              f\"GLL={best_score:.6f}\")\n    return best, best_score\n\ndef tune_postprocess_hparams_wide(\n    trained, y_true, naive_mean, naive_sigma,\n    n_splits=5, seed=42, verbose=True,\n    extra_trials=300,\n    clip_random_trials=800,\n    grid_max_pairs=400,\n    fgs_lo_range=(0.45, 1.20, 0.05),\n    fgs_hi_range=(1.00, 2.00, 0.05),\n    air_lo_range=(0.70, 1.05, 0.03),\n    air_hi_range=(1.00, 1.70, 0.05),\n    sigma_scale_range=(1.50, 3.0, 0.02),\n):\n    best, best_score = tune_postprocess_hparams(\n        trained, y_true, naive_mean, naive_sigma,\n        n_splits=n_splits, seed=seed, verbose=verbose\n    )\n    oof_raw = np.asarray(trained[\"oof_predictions\"], dtype=np.float32)\n    def _step_sample(rng):\n        lo, hi, step = rng; n = int(np.floor((hi - lo) / step))\n        k = np.random.randint(0, n + 1); return float(lo + k * step)\n    # 1) 풀 랜덤\n    for _ in range(int(extra_trials)):\n        w = int(np.random.choice([7, 9, 11, 13, 15, 17, 19]))\n        p = int(np.random.choice([2, 3])); \n        if p >= w:  continue\n        a = float(np.random.uniform(0.50, 0.80))\n        fgs_lo = _step_sample(fgs_lo_range); fgs_hi = _step_sample(fgs_hi_range)\n        air_lo = _step_sample(air_lo_range); air_hi = _step_sample(air_hi_range)\n        if not (fgs_lo < fgs_hi and air_lo < air_hi): continue\n        s = _step_sample(sigma_scale_range)\n        sc, sd = _cv_eval_postprocess(\n            oof_raw, y_true, w, p, a,\n            (fgs_lo, fgs_hi), (air_lo, air_hi), s,\n            naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n        )\n        if sc > best_score:\n            best_score = sc\n            best.update(window=w, poly=p, alpha=a,\n                        clip_fgs=(fgs_lo, fgs_hi),\n                        clip_airs=(air_lo, air_hi),\n                        sigma_scale=s)\n            if verbose:\n                print(f\"[WIDE-RAND] win={w}, poly={p}, alpha={a}, \"\n                      f\"FGS={best['clip_fgs']}, AIRS={best['clip_airs']}, scale={s} \"\n                      f\"-> GLL={sc:.6f} (±{sd:.6f})\")\n    # 1.5) 클립만 랜덤\n    for _ in range(int(clip_random_trials)):\n        f_lo = _step_sample(fgs_lo_range); f_hi = _step_sample(fgs_hi_range)\n        a_lo = _step_sample(air_lo_range); a_hi = _step_sample(air_hi_range)\n        if not (f_lo < f_hi and a_lo < a_hi): continue\n        sc, sd = _cv_eval_postprocess(\n            oof_raw, y_true, best[\"window\"], best[\"poly\"], best[\"alpha\"],\n            (f_lo, f_hi), (a_lo, a_hi), best[\"sigma_scale\"],\n            naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n        )\n        if sc > best_score:\n            best_score = sc; best[\"clip_fgs\"] = (f_lo, f_hi); best[\"clip_airs\"] = (a_lo, a_hi)\n            if verbose: print(f\"[WIDE-CLIP-RAND] FGS={best['clip_fgs']}, AIRS={best['clip_airs']} -> GLL={sc:.6f} (±{sd:.6f})\")\n    # 2) 코스 그리드\n    def _linspace_with_step(lo, hi, step, max_points=12):\n        num = int(np.floor((hi - lo) / step)) + 1\n        if num <= max_points: return [lo + i * step for i in range(num)]\n        stride = int(np.ceil(num / max_points))\n        return [lo + i * step for i in range(0, num, stride)]\n    fgs_lo_grid = _linspace_with_step(*fgs_lo_range, max_points=10)\n    fgs_hi_grid = _linspace_with_step(*fgs_hi_range, max_points=10)\n    air_lo_grid = _linspace_with_step(*air_lo_range, max_points=10)\n    air_hi_grid = _linspace_with_step(*air_hi_range, max_points=10)\n    tried = 0\n    for _ in range(int(grid_max_pairs)):\n        f_lo = float(np.random.choice(fgs_lo_grid)); f_hi = float(np.random.choice(fgs_hi_grid))\n        a_lo = float(np.random.choice(air_lo_grid)); a_hi = float(np.random.choice(air_hi_grid))\n        if not (f_lo < f_hi and a_lo < a_hi): continue\n        tried += 1\n        sc, sd = _cv_eval_postprocess(\n            oof_raw, y_true, best[\"window\"], best[\"poly\"], best[\"alpha\"],\n            (f_lo, f_hi), (a_lo, a_hi), best[\"sigma_scale\"],\n            naive_mean, naive_sigma, n_splits=n_splits, seed=seed\n        )\n        if sc > best_score:\n            best_score = sc; best[\"clip_fgs\"] = (f_lo, f_hi); best[\"clip_airs\"] = (a_lo, a_hi)\n            if verbose: print(f\"[WIDE-GRID] FGS={best['clip_fgs']}, AIRS={best['clip_airs']} -> GLL={sc:.6f} (±{sd:.6f})\")\n    if verbose:\n        print(\"\\n[WIDE-BEST] \"\n              f\"win={best['window']}, poly={best['poly']}, alpha={best['alpha']}, \"\n              f\"FGS={best['clip_fgs']}, AIRS={best['clip_airs']}, scale={best['sigma_scale']}, \"\n              f\"GLL={best_score:.6f}\")\n    return best, best_score\n\n# ----------------------- Main -----------------------\ndef _nrows(a):\n    try: return a.shape[0]\n    except Exception: return len(a)\n\ndef _take_rows(a, mask):\n    import numpy as _np\n    if hasattr(a, \"iloc\"): return a.iloc[mask]\n    else: return a[_np.asarray(mask)]\n\nif __name__ == \"__main__\":\n    t0 = time.time()\n    if USE_CACHED_XY:\n        print(\"== Load TRAIN (cached X,y) ==\")\n        X, y, WL_COLS = load_cached_xy_csv(CACHED_TRAIN_XY)  # 자동 탐지\n        print(f\"[SHAPE] X:{X.shape}, y:{y.shape}\")\n        arr_y = y.values.astype(np.float32)\n        NAIVE_MEAN  = float(np.nanmean(arr_y))\n        NAIVE_SIGMA = float(np.nanstd(arr_y) + 1e-12)\n    else:\n        print(\"== Build TRAIN features (from images) ==\")\n        X, y = build_features(split=\"train\", star_info_df=train_star, targets_df=train_df, frac=float(TRAIN_FRAC))\n        arr_y = y.values.astype(np.float32)\n        NAIVE_MEAN  = float(np.nanmean(arr_y))\n        NAIVE_SIGMA = float(np.nanstd(arr_y) + 1e-12)\n\n    # (A) Baseline 학습 → 중요도\n    print(\"== Baseline train (for FI) ==\")\n    baseline_trained = train_lgbm_cv(X, y)\n    fi = get_mean_feature_importance(baseline_trained, kind=\"gain\")\n\n    # (B) Top-K 기반 파생 생성\n    TOPK_DERIVE = 50; MEAN_THRESH = 1.0\n    X_derived = X.copy()\n    X_derived, DERIVED_INFO, DERIVED_COLS = build_derived_from_topk(\n        X_derived, fi, topk=TOPK_DERIVE, mean_threshold=MEAN_THRESH, eps=1e-9\n    )\n    X_derived.replace([np.inf, -np.inf], np.nan, inplace=True); X_derived.fillna(0, inplace=True)\n    DERIVED_COLS = filter_similar_derived(\n        X_derived, DERIVED_COLS, fi_df=fi, thresh=0.90, max_check=5000\n    )\n    DERIVED_INFO[\"new_cols\"] = DERIVED_COLS\n\n    # (B-1) (옵션) GMM 증강\n    USE_GMM_AUG  = True\n    K_CLUSTERS   = 4\n    AUG_MULT     = 0.30\n    GM_COMPONENTS= 2\n    Y_PCA_DIM    = 8\n    if USE_GMM_AUG:\n        print(\"== Cluster by volatility & GMM-augment (on derived) ==\")\n        labels, clust_model, aug = cluster_and_gmm_augment(\n            X_derived, y, k_clusters=K_CLUSTERS, aug_mult=AUG_MULT,\n            gm_components=GM_COMPONENTS, y_pca_dim=Y_PCA_DIM, random_state=SEED\n        )\n        X_train, y_train = aug[\"X_aug\"], aug[\"y_aug\"]\n    else:\n        X_train, y_train = X_derived, y\n\n    # (B-2) 파생(증강) 포함 재학습\n    print(\"== Retrain with derived features ==\")\n    trained = train_lgbm_cv(X_train, y_train)\n    trained[\"derived_info\"] = DERIVED_INFO; trained[\"derived_cols\"] = DERIVED_COLS\n\n    # (C) 후보정(초기값) — 이후 WIDE 튜너로 갱신 예정\n    trained = fit_postcalibration(\n        trained, y_train.values,\n        smooth_window=POST_SMOOTH_WINDOW, smooth_poly=POST_SMOOTH_POLY, blend_alpha=POST_BLEND_ALPHA,\n        clip_fgs=SIGMA_CLIP_FGS, clip_airs=SIGMA_CLIP_AIRS, final_sigma_scale=SIGMA_FINAL_SCALE\n    )\n\n    # (D) Template 라이브러리 & PCA 서브스페이스 회귀\n    print(\"== Build template library & PCA subspace regressor ==\")\n    trained[\"template_lib\"] = build_template_library(y_train.values, k_pca=12, n_neighbors=8, seed=SEED)\n    trained[\"pca_pack\"]     = train_pca_coefficient_regressor(X_train, y_train.values, k=24, alpha=1.0, seed=SEED)\n    trained[\"max_shift\"]    = 3  # AIRS 정렬 허용 시프트\n\n    # (E) 폭넓은 후보정 튜닝 → 최종 보정 반영\n    print(\"== Tune postprocess (WIDE) ==\")\n    oof_base = trained[\"oof_predictions\"]\n    n_oof = _nrows(oof_base)\n    orig_n = len(y); aug_n  = len(y_train)\n    if n_oof == aug_n:\n        y_for_trained = np.asarray(y_train.values if hasattr(y_train, \"values\") else y_train)\n        mask_orig = np.arange(n_oof) < orig_n\n    elif n_oof == orig_n:\n        y_for_trained = np.asarray(y.values if hasattr(y, \"values\") else y)\n        mask_orig = np.ones(n_oof, dtype=bool)\n    else:\n        raise ValueError(f\"Length mismatch: OOF={n_oof}, y={orig_n}, y_train={aug_n}\")\n\n    best_params, best_cv = tune_postprocess_hparams_wide(\n        trained, y_for_trained, NAIVE_MEAN, NAIVE_SIGMA,\n        n_splits=N_SPLITS, seed=SEED, verbose=True,\n        extra_trials=300, clip_random_trials=800, grid_max_pairs=400,\n        fgs_lo_range=(0.45, 1.20, 0.05),\n        fgs_hi_range=(1.00, 2.00, 0.05),\n        air_lo_range=(0.70, 1.05, 0.03),\n        air_hi_range=(1.00, 1.70, 0.05),\n        sigma_scale_range=(1.50, 3.0, 0.02),\n    )\n\n    trained = fit_postcalibration(\n        trained, y_for_trained,\n        smooth_window=best_params[\"window\"], smooth_poly=best_params[\"poly\"], blend_alpha=best_params[\"alpha\"],\n        clip_fgs=best_params[\"clip_fgs\"], clip_airs=best_params[\"clip_airs\"],\n        final_sigma_scale=best_params[\"sigma_scale\"],\n    )\n\n    # (F) 원본 구간만 OOF 점수\n    oof_post = trained[\"oof_predictions_post\"]\n    oof_mu_post = _take_rows(oof_post, mask_orig)\n    calib_sig = np.asarray(trained[\"calibrated_uncertainties\"]).reshape(-1)\n    oof_sigma = np.tile(calib_sig, (mask_orig.sum(), 1))\n    oof_score = gll_score_numpy(\n        np.asarray(y.values), np.asarray(oof_mu_post), oof_sigma, NAIVE_MEAN, NAIVE_SIGMA, fgs_weight=0.4\n    )\n    print(f\"[OOF/Post tuned] GLL: {oof_score:.6f}\")\n\n    # (G) 개선 정도 점검 + 플롯 저장\n    mu_oof_improved, _ = predict_with_uncertainty(trained, X_train)\n    delta = float(np.mean(np.abs(mu_oof_improved - trained[\"oof_predictions_post\"])))\n    print(f\"[OOF] improved-vs-post mean|Δ| = {delta:.3e}\")\n\n    try:\n        plot_oof_lines(y_for_trained, mu_oof_improved, wavelengths, trained[\"index\"], PLOT_DIR, batch_size=12, show=False)\n        plot_bias_sigma(trained, wavelengths, PLOT_DIR)\n        plot_residual_std_vs_sigma(y_for_trained, trained, wavelengths, PLOT_DIR)\n        plot_scatter_and_residuals(y_for_trained, mu_oof_improved, wavelengths, PLOT_DIR, picks=(0,50,150,250), show=False)\n        print(f\"[PLOTS] saved under {PLOT_DIR}\")\n    except Exception as e:\n        print(\"[PLOTS] failed:\", e)\n\n    # (H) 캐시 저장\n    cache_path = f\"{WORK_DIR}/fgs1_train_xy.csv\"\n    save_cached_xy_csv(X_train, y_train, cache_path, wl_cols=WL_COLS)\n\n    print(f\"[DONE] total {time.time()-t0:.1f}s\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-23T11:08:01.871469Z","iopub.execute_input":"2025-09-23T11:08:01.872085Z","iopub.status.idle":"2025-09-23T11:25:04.838316Z","shell.execute_reply.started":"2025-09-23T11:08:01.872058Z","shell.execute_reply":"2025-09-23T11:25:04.837168Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ===== Imports (필수) =====\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom scipy.spatial.distance import cdist\nfrom sklearn.preprocessing import StandardScaler, RobustScaler\nfrom sklearn.decomposition import PCA\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.metrics import roc_auc_score\nfrom sklearn.neighbors import NearestNeighbors\n\n# --- 커널 MMD (Gaussian RBF) ---\ndef _rbf_kernel(X, Y=None, gamma=None):\n    X = np.asarray(X, dtype=np.float64)\n    Y = X if Y is None else np.asarray(Y, dtype=np.float64)\n    if gamma is None:\n        # median heuristic (0 방지)\n        D = cdist(X, Y, metric=\"euclidean\")\n        D = D[np.isfinite(D) & (D > 0)]\n        med = np.median(D) if D.size else 1.0\n        gamma = 1.0 / (2.0 * (med**2 + 1e-12))\n    return np.exp(-gamma * cdist(X, Y, metric=\"sqeuclidean\"))\n\ndef mmd_rbf(X, Y, gamma=None):\n    X = np.asarray(X, dtype=np.float64)\n    Y = np.asarray(Y, dtype=np.float64)\n    n, m = X.shape[0], Y.shape[0]\n    if n < 2 or m < 2:\n        return 0.0\n    Kxx = _rbf_kernel(X, None, gamma)\n    Kyy = _rbf_kernel(Y, None, gamma)\n    Kxy = _rbf_kernel(X, Y,  gamma)\n    mmd2 = (Kxx.sum() - np.trace(Kxx)) / (n*(n-1) + 1e-12) \\\n         + (Kyy.sum() - np.trace(Kyy)) / (m*(m-1) + 1e-12) \\\n         - 2.0 * Kxy.mean()\n    return float(max(mmd2, 0.0))\n\n# --- 증강 마스크 감지: 인덱스가 'synth_'로 시작하면 증강으로 간주 ---\ndef _is_aug_mask(df_like, n_all: int = None):\n    \"\"\"\n    인덱스가 'synth_'로 시작하는지 보고 증강샘플 마스크를 만든다.\n    - df_like: DataFrame/Series/ndarray 모두 허용\n    - n_all  : 전체 길이(옵션). None이면 len(df_like)로 추정\n    \"\"\"\n    try:\n        if hasattr(df_like, \"index\"):\n            idx = pd.Index(getattr(df_like, \"index\"))\n            mask = pd.Series(idx.astype(str), dtype=\"string\") \\\n                     .str.startswith(\"synth_\") \\\n                     .fillna(False) \\\n                     .to_numpy(dtype=bool)\n            return mask\n    except Exception:\n        pass\n    if n_all is None:\n        try:\n            n_all = len(df_like)\n        except Exception:\n            n_all = 0\n    return np.zeros(int(n_all), dtype=bool)\n\ndef diagnose_aug(\n    X, y, X_train, y_train,\n    *, use_y=True, robust=True, whiten=True, max_n=4000, random_state=42\n):\n    \"\"\"\n    - PCA 2D로 원본/증강 분포를 바로 띄움(plt.show)\n    - 원본 vs 증강 분리 AUC, RBF-MMD^2, 최근접거리 히스토그램 표시\n    \"\"\"\n    rng = np.random.default_rng(random_state)\n\n    def _to_matrix(A):\n        if hasattr(A, \"values\"): A = A.values\n        A = np.asarray(A)\n        if A.ndim == 1: A = A.reshape(-1, 1)\n        return A.astype(np.float32)\n\n    Y_all_df  = (y_train if use_y else X_train)\n    Y_orig_df = (y if use_y else X)\n\n    Y_all  = _to_matrix(Y_all_df)\n    Y_orig = _to_matrix(Y_orig_df)\n\n    is_aug = _is_aug_mask(Y_all_df, len(Y_all))\n    has_aug = bool(is_aug.any())\n\n    # ----- 스케일 & PCA -----\n    scaler = RobustScaler() if robust else StandardScaler()\n    Ys = scaler.fit_transform(Y_all)\n\n    n_comp = int(max(2, min(10, Ys.shape[1], Ys.shape[0]-1)))\n    pca = PCA(n_components=n_comp, whiten=whiten, random_state=random_state)\n    Z = pca.fit_transform(Ys)\n    evr = pca.explained_variance_ratio_\n    print(f\"[PCA] n_comp={pca.n_components_}, EVR top2={evr[:2].round(4)}, cum{n_comp}={evr[:n_comp].sum():.4f}\")\n    print(f\"[diagnose_aug] n_all={len(Y_all)}, n_orig={len(Y_orig)}, n_synth={int(is_aug.sum())}\")\n\n    # ----- 2D 시각화 -----\n    plt.figure(figsize=(6.8,5.2))\n    idx = np.arange(len(Z))\n    if len(idx) > max_n:\n        idx = rng.choice(idx, size=max_n, replace=False)\n    if has_aug:\n        base = idx[~is_aug[idx]]\n        synth = idx[is_aug[idx]]\n        plt.scatter(Z[base,0],  Z[base,1],  s=8, alpha=0.65, label=\"original\")\n        plt.scatter(Z[synth,0], Z[synth,1], s=10, alpha=0.75, marker='x', label=\"synthetic\")\n    else:\n        plt.scatter(Z[idx,0], Z[idx,1], s=8, alpha=0.7, label=\"original (no synthetic)\")\n    plt.title(f\"PCA 2D of {'y' if use_y else 'X'}  — robust={robust}, whiten={whiten}\")\n    plt.legend()\n    plt.tight_layout()\n    plt.show()\n\n    # ----- 원본/증강 분리 가능도 -----\n    if has_aug and np.unique(is_aug.astype(int)).size >= 2:\n        k = int(min(20, Z.shape[1]))\n        clf = LogisticRegression(max_iter=300, solver=\"lbfgs\")\n        clf.fit(Z[:,:k], is_aug.astype(int))\n        auc = roc_auc_score(is_aug.astype(int), clf.predict_proba(Z[:,:k])[:,1])\n        print(f\"[AUC(original vs synthetic) @ PCA-{k}] = {auc:.3f}\")\n    else:\n        print(\"[AUC] skipped (no synthetic or one class)\")\n\n    # ----- 분포 거리 & 최근접 거리 -----\n    if has_aug:\n        def _sub(A, n=max_n):\n            if len(A) <= n: return A\n            return A[rng.choice(len(A), n, replace=False)]\n        A = _sub(Ys[~is_aug])\n        B = _sub(Ys[is_aug])\n        if len(A) > 10 and len(B) > 10:\n            mmd = mmd_rbf(A, B)\n            mean_dist = float(np.linalg.norm(A.mean(axis=0) - B.mean(axis=0)))\n            print(f\"[Distance] MMD^2={mmd:.6f} (↓유사), MeanVecDist={mean_dist:.6f}\")\n        else:\n            print(\"[Distance] skipped (too few samples)\")\n\n        if (~is_aug).any():\n            nn = NearestNeighbors(n_neighbors=1, metric=\"euclidean\").fit(Ys[~is_aug])\n            dist, _ = nn.kneighbors(Ys[is_aug])\n            plt.figure(figsize=(6.8,4.2))\n            plt.hist(dist.ravel(), bins=40, alpha=0.85, density=True)\n            plt.title(\"Nearest-original distance of synthetic samples\")\n            plt.xlabel(\"euclidean distance\"); plt.ylabel(\"density\")\n            plt.tight_layout()\n            plt.show()\n        else:\n            print(\"[NN] skipped (no original rows)\")\n    else:\n        print(\"[NN/Distance] skipped (no synthetic)\")\n    \n# ===== 사용 예시 =====\nif len(X_train) > len(X):\n    diagnose_aug(X, y, X_train, y_train, use_y=True,  robust=True,  whiten=True)\n    diagnose_aug(X, y, X_train, y_train, use_y=False, robust=True,  whiten=True)\nelse:\n    print(\"[diagnose_aug] skipped: no synthetic rows\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-23T11:25:04.839863Z","iopub.execute_input":"2025-09-23T11:25:04.840277Z","iopub.status.idle":"2025-09-23T11:25:05.747396Z","shell.execute_reply.started":"2025-09-23T11:25:04.840242Z","shell.execute_reply":"2025-09-23T11:25:05.746373Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# === Visualization ===\ntry:\n    ids_for_plot = trained.get(\"index\", [f\"pid_{i}\" for i in range(oof_mu_post.shape[0])])\n\n    # 개별 스펙트럼 라인 (배치 저장)\n    plot_oof_lines(\n        y_true=y.values,\n        y_pred=oof_mu_post,\n        wavelengths=wavelengths,\n        ids=ids_for_plot,\n        out_dir=PLOT_DIR,\n        batch_size=12,\n        show=True\n    )\n\n    # 몇 개 λ에서 산포/잔차 분포\n    plot_scatter_and_residuals(\n        y_true=y.values,\n        y_pred=oof_mu_post,\n        wavelengths=wavelengths,\n        out_dir=PLOT_DIR,\n        picks=(0, 50, 150, 250),\n        show=True\n    )\n\n    # λ별 바이어스/σ/스케일 추적\n    plot_bias_sigma(trained, wavelengths, PLOT_DIR)\n\n    # λ별 잔차 robust std와 예측 σ 비교\n    plot_residual_std_vs_sigma(y.values, trained, wavelengths, PLOT_DIR)\n\n    print(f\"[PLOTS] saved to: {PLOT_DIR}\")\nexcept Exception as e:\n    print(\"[PLOTS] failed:\", e)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-23T11:25:05.749723Z","iopub.execute_input":"2025-09-23T11:25:05.75009Z","iopub.status.idle":"2025-09-23T11:25:23.31989Z","shell.execute_reply.started":"2025-09-23T11:25:05.750063Z","shell.execute_reply":"2025-09-23T11:25:23.318823Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# == Make submission ==\nsubmission_path = f\"{WORK_DIR}/submission.csv\"  # Kaggle는 /kaggle/working/submission.csv 를 찾음\nsubmission = make_submission(trained, test_star, out_path=submission_path)\nprint(f\"[DONE] saved: {submission_path}  shape={submission.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-23T11:25:23.320626Z","iopub.execute_input":"2025-09-23T11:25:23.320988Z","iopub.status.idle":"2025-09-23T11:25:41.650883Z","shell.execute_reply.started":"2025-09-23T11:25:23.320955Z","shell.execute_reply":"2025-09-23T11:25:41.649855Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"submission","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-23T11:25:41.652283Z","iopub.execute_input":"2025-09-23T11:25:41.652653Z","iopub.status.idle":"2025-09-23T11:25:41.686059Z","shell.execute_reply.started":"2025-09-23T11:25:41.652619Z","shell.execute_reply":"2025-09-23T11:25:41.685039Z"}},"outputs":[],"execution_count":null}]}