{"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":9806277,"sourceType":"datasetVersion","datasetId":5998787}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# install pqdm for parallel processing\n!pip install --no-index --find-links=/kaggle/input/ariel-2024-pqdm pqdm","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T00:14:02.781431Z","iopub.execute_input":"2025-09-18T00:14:02.78192Z","iopub.status.idle":"2025-09-18T00:14:08.614457Z","shell.execute_reply.started":"2025-09-18T00:14:02.781885Z","shell.execute_reply":"2025-09-18T00:14:08.613337Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Ariel Data Challenge 2025 — baseline analytics + safe ML calibration\n# ---------------------------------------------------------------\n# This script augments the user's analytical pipeline with a light-weight,\n# conservative ML layer (OOF-calibrated scale, low-rank spectral shape via PCA + Ridge,\n# and sigma calibration). If ML fails at any point, it gracefully falls back to the\n# original baseline behavior.\n#\n# Notes:\n# - Designed for Kaggle environment paths/layout.\n# - Keeps your preprocessing and TransitModel logic intact.\n# - Uses strong shrinkage (alpha/beta/gamma <= 0.3) to avoid degradation.\n# ---------------------------------------------------------------\n\nimport os\nimport time\nimport itertools\nimport numpy as np\nimport pandas as pd\n\nfrom tqdm import tqdm\nfrom pqdm.threads import pqdm\n\nfrom scipy.optimize import minimize\nfrom scipy.signal import savgol_filter, find_peaks\nfrom sklearn.linear_model import Ridge\nfrom sklearn.decomposition import PCA\nfrom sklearn.model_selection import KFold\nfrom sklearn.metrics import mean_squared_error\n\nfrom astropy.stats import sigma_clip\n\nfrom scipy.stats import skew, kurtosis\nfrom sklearn.preprocessing import StandardScaler\n\nimport matplotlib.pyplot as plt\n\n# ------------------------------\n# Config\n# ------------------------------\n\nclass Config:\n DATA_PATH = '/kaggle/input/ariel-data-challenge-2025'\n DATASET_TEST = 'test'\n DATASET_TRAIN = 'train'\n\n SCALE = 0.952\n SIGMA = 0.00055\n \n\n CUT_INF = 39\n CUT_SUP = 321\n \n def __init__(self):\n # Load wavelength information\n try:\n self.wavelengths_df = pd.read_csv(f\"{self.DATA_PATH}/wavelengths.csv\")\n # Extract actual wavelength values from the single row\n self.wavelengths = self.wavelengths_df.iloc[0].values\n print(f\"[INFO] Loaded {len(self.wavelengths)} wavelength points: {self.wavelengths[0]:.3f} - {self.wavelengths[-1]:.3f} µm\")\n except Exception as e:\n print(f\"[WARN] Could not load wavelengths.csv: {e}\")\n self.wavelengths = None\n self.wavelengths_df = None\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\n MODEL_POLYNOMIAL_DEGREE = 3\n\n N_JOBS = 3\n\n# ------------------------------\n# Utilities\n# ------------------------------\n\n\ndef smooth_spectrum(spec, win=7, poly=2):\n # spec: (n_ch,) or (n_bins, n_ch) — smooth along last axis\n if spec.ndim == 1:\n L = len(spec)\n win = min(win, L if L%2==1 else L-1)\n if win < 3:\n return spec\n return savgol_filter(spec, window_length=win, polyorder=poly, mode='mirror')\n else:\n out = np.empty_like(spec)\n for i in range(spec.shape[0]):\n out[i] = smooth_spectrum(spec[i], win=win, poly=poly)\n return out\ndef apply_spectral_smoothness_constraint(mu_predictions, sigma_predictions, lambda_smooth=0.1):\n \"\"\"Apply smoothness constraint across wavelengths for physically realistic transmission spectra\"\"\"\n \n n_samples, n_mu = mu_predictions.shape\n mu_smooth = mu_predictions.copy()\n \n for i in range(n_samples):\n spectrum = mu_predictions[i, :]\n weights = 1.0 / (sigma_predictions[i, :] + 1e-8) # Inverse variance weighting\n \n # Weighted smoothing spline\n try:\n from scipy.interpolate import UnivariateSpline\n \n wl_indices = np.arange(n_mu)\n # Higher smoothing for noisier spectra\n smoothing_factor = lambda_smooth * n_mu * np.mean(sigma_predictions[i, :])\n \n spline = UnivariateSpline(wl_indices, spectrum, w=weights, s=smoothing_factor)\n spectrum_smooth = spline(wl_indices)\n \n # Conservative blend with original (20% smoothed, 80% original)\n blend_factor = 0.2\n mu_smooth[i, :] = (1 - blend_factor) * spectrum + blend_factor * spectrum_smooth\n \n except Exception as e:\n # If smoothing fails (e.g., scipy not available), use simpler approach\n try:\n # Fallback: simple moving average with weights\n window_size = max(3, min(5, n_mu // 3))\n if window_size % 2 == 0:\n window_size += 1\n \n spectrum_smooth = savgol_filter(spectrum, window_length=window_size, polyorder=2)\n \n # Conservative blend\n blend_factor = 0.15 # Even more conservative for fallback\n mu_smooth[i, :] = (1 - blend_factor) * spectrum + blend_factor * spectrum_smooth\n \n except:\n # If all smoothing fails, keep original\n pass\n \n return mu_smooth\n\ndef detect_malformed_transits(preprocessed_signals, cfg):\n \"\"\"Detect incomplete/malformed transits that hurt performance\"\"\"\n quality_scores = []\n \n for single in preprocessed_signals:\n fgs = single[:, 0]\n \n # Basic quality checks\n quality = 1.0 # Start with perfect quality\n \n # 1. Check for sufficient data length\n if len(fgs) < 50:\n quality *= 0.3\n \n # 2. Check for excessive NaN/inf values\n nan_frac = np.sum(~np.isfinite(fgs)) / len(fgs)\n if nan_frac > 0.1:\n quality *= (1.0 - nan_frac)\n \n # 3. Check transit depth reasonableness\n try:\n # Use your existing phase detection\n smooth_fgs = savgol_filter(fgs, min(21, len(fgs)//3*2+1), 2)\n p1, p2 = _phase_detector_signal(smooth_fgs, cfg)\n \n if p2 > p1 + 5:\n transit_region = fgs[p1:p2]\n oot_region = np.concatenate([fgs[:p1], fgs[p2:]])\n \n if len(oot_region) > 0 and len(transit_region) > 0:\n depth = np.nanmean(oot_region) - np.nanmean(transit_region)\n depth_std = np.nanstd(oot_region)\n \n # Transit should be detectable (SNR > 2)\n if depth_std > 0:\n snr = depth / depth_std\n if snr < 2.0:\n quality *= 0.5\n elif snr > 15.0: # Too good to be true\n quality *= 0.7\n \n # Depth should be reasonable (0.001% to 5%)\n relative_depth = depth / (np.nanmean(oot_region) + 1e-9)\n if relative_depth < 1e-5 or relative_depth > 0.05:\n quality *= 0.6\n else:\n # No clear transit detected\n quality *= 0.4\n \n except:\n quality *= 0.5\n \n # 4. Check for systematic trends (detrending artifacts)\n try:\n # Polynomial detrending residuals should be small\n x = np.arange(len(fgs))\n if len(fgs) > 6:\n poly_fit = np.polyval(np.polyfit(x, fgs, 3), x)\n residuals = fgs - poly_fit\n trend_strength = np.std(poly_fit) / (np.std(residuals) + 1e-9)\n \n if trend_strength > 5.0: # Strong systematic trends\n quality *= 0.8\n except:\n pass\n \n # 5. Check for data completeness and regularity\n try:\n # Check for large gaps or irregular sampling\n if len(fgs) > 10:\n # Simple completeness check - are there big jumps in \"time\"?\n # Assume roughly regular sampling\n expected_std = np.std(fgs) * 0.1 # Expected noise level\n actual_std = np.std(fgs)\n \n # If much noisier than expected, reduce quality\n if actual_std > expected_std * 5:\n quality *= 0.9\n except:\n pass\n \n quality_scores.append(np.clip(quality, 0.1, 1.0))\n \n return np.array(quality_scores)\n\ndef find_ingress_egress(signal_1d, coarse_steps=(range(30,130,20)), refine_radius=20):\n # signal_1d: 1D time-series average across wavelengths (already smoothed)\n n = len(signal_1d)\n best = None\n def fit_piecewise(p1, p2):\n # fit quad on [0:p1], quad on [p2:n_half], connect with line between\n # but we only need error metric: use polyfit on each segment and RMSE\n seg1 = signal_1d[:p1]\n seg2 = signal_1d[p2:]\n if len(seg1) < 3 or len(seg2) < 3:\n return np.inf\n x1 = np.arange(len(seg1))\n x2 = np.arange(len(seg2))\n c1 = np.polyfit(x1, seg1, 2)\n c2 = np.polyfit(x2 + p2, seg2, 2)\n pred1 = np.polyval(c1, x1)\n pred2 = np.polyval(c2, x2 + p2)\n err = ((seg1 - pred1)**2).mean() + ((seg2 - pred2)**2).mean()\n return err\n\n # coarse search (assume ingress in first half)\n half = n//2\n best_err = np.inf\n for p1 in coarse_steps:\n for p2 in coarse_steps:\n if p1 < p2 and p2 < half:\n err = fit_piecewise(p1, p2)\n if err < best_err:\n best_err = err\n best = (p1, p2)\n\n # refine around best\n if best is None:\n return 0, n-1\n p1c, p2c = best\n rng1 = range(max(2, p1c - refine_radius), min(p1c + refine_radius, half))\n rng2 = range(max(rng1.stop+1, p2c - refine_radius), min(p2c + refine_radius, half))\n for p1 in rng1:\n for p2 in rng2:\n if p1 < p2:\n err = fit_piecewise(p1, p2)\n if err < best_err:\n best_err = err\n best = (p1, p2)\n return best\n\nfrom sklearn.decomposition import NMF\n\ndef nmf_denoise(spectrum, rank=3, max_iter=200):\n # spectrum: (n_ch,) or (n_samples, n_ch)\n single = False\n if spectrum.ndim == 1:\n spectrum = spectrum[None, :]\n single = True\n # ensure nonnegativity: shift by min if needed\n mn = spectrum.min()\n shift = 0.0\n if mn < 0:\n shift = -mn\n spectrum = spectrum + shift\n model = NMF(n_components=rank, init='nndsvda', max_iter=max_iter, tol=1e-4)\n W = model.fit_transform(spectrum)\n recon = model.inverse_transform(W)\n recon = recon - shift\n return recon[0] if single else recon\n\ndef compose_sigma(const_term, std_smoothed, gp_std=None, w0=1.0, a=1.0, b=1.0):\n # std_smoothed: (N, n_ch) or (N,) depending on how computed\n base = (w0 * const_term)**2\n s1 = (a * std_smoothed)**2\n if gp_std is not None:\n s2 = (b * gp_std)**2\n out = np.sqrt(base + s1 + s2)\n else:\n out = np.sqrt(base + s1)\n return out\n\ndef bootstrap_std_estimate(single_preprocessed, dip_estimator, n_boot=30, random_state=0):\n rng = np.random.RandomState(random_state)\n n = single_preprocessed.shape[0]\n dips = []\n for _ in range(n_boot):\n idx = rng.choice(n, size=n, replace=True)\n sample = single_preprocessed[idx]\n dips.append(dip_estimator(sample)) # returns vector (n_ch,)\n dips = np.vstack(dips)\n return np.std(dips, axis=0, ddof=1)\n\ndef compute_c_star_per_planet(y_true, mu_pred, sigma_shape):\n # y_true, mu_pred, sigma_shape: (N_planets, n_ch)\n eps = 1e-12\n R2_div_S2 = ((y_true - mu_pred)**2) / (np.clip(sigma_shape, eps, None)**2)\n c_star = np.sqrt(np.nanmean(R2_div_S2, axis=1))\n return c_star # shape (N,)\n\ndef _phase_detector_signal(signal: np.ndarray, cfg: Config):\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# Sigma estimators (baseline, from user's code)\n# ------------------------------\n\ndef estimate_sigma_fgs(preprocessed_data: np.ndarray, cfg: Config):\n \"\"\"Returns vector sigma_1 (for FGS1) length N_planets — soft multiplier to cfg.SIGMA.\"\"\"\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: np.ndarray, cfg: Config):\n \"\"\"Returns vector sigma_air length N_planets — soft multiplier to cfg.SIGMA for all AIRS channels.\"\"\"\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 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# Signal processing (from user, with minor generalization for dataset param)\n# ------------------------------\n\nclass SignalProcessor:\n def __init__(self, config: Config):\n self.cfg = config\n self.adc_info = pd.read_csv(f\"{self.cfg.DATA_PATH}/adc_info.csv\")\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, dataset):\n sensor_cfg = self.cfg.SENSOR_CONFIG[sensor]\n base = f\"{self.cfg.DATA_PATH}/{dataset}/{planet_id}\"\n\n signal = pd.read_parquet(f\"{base}/{sensor}_signal_0.parquet\").to_numpy()\n dark = pd.read_parquet(f\"{base}/{sensor}_calibration_0/dark.parquet\").to_numpy()\n dead = pd.read_parquet(f\"{base}/{sensor}_calibration_0/dead.parquet\").to_numpy()\n flat = pd.read_parquet(f\"{base}/{sensor}_calibration_0/flat.parquet\").to_numpy()\n linear_corr = pd.read_parquet(f\"{base}/{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 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 return signal\n\n def _detect_transit_mask(self, white_light):\n \"\"\"Detect approximate transit phase for phase-aware processing\"\"\"\n if len(white_light) < 10:\n return np.zeros(len(white_light), dtype=bool)\n \n # Simple transit detection: find the dimmest region\n smooth_wl = savgol_filter(white_light, window_length=min(7, len(white_light)//3*2+1), polyorder=2)\n median_val = np.nanmedian(smooth_wl)\n # Transit typically causes 0.1-2% dimming\n threshold = median_val - 2 * np.nanstd(smooth_wl)\n return smooth_wl < threshold\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 # Robust binning with sigma clipping to remove outliers\n binned = []\n for j in range(n_bins):\n chunk = cds_signal[j*binning : (j+1)*binning]\n # Apply sigma clipping to remove outliers, then take mean\n clipped = sigma_clip(chunk, sigma=2.5, axis=0, masked=False, maxiters=2)\n binned.append(np.nanmean(clipped, axis=0))\n binned = np.array(binned)\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 # PCA spectral denoising to preserve wavelength correlations\n n_bins, n_wl = binned.shape\n if n_wl > 5 and n_bins > 5 and not np.any(~np.isfinite(binned)):\n try:\n pca = PCA(n_components=min(8, n_wl//2), svd_solver='randomized', random_state=42)\n binned_pca = pca.fit_transform(binned)\n binned_reconstructed = pca.inverse_transform(binned_pca)\n # Conservative blend: 70% original, 30% denoised\n alpha = 0.7\n binned = alpha * binned + (1-alpha) * binned_reconstructed\n except:\n # If PCA fails, keep original\n pass\n \n # Phase-aware processing: different smoothing for transit vs out-of-transit\n white_light = np.nanmean(binned, axis=1) if binned.ndim > 1 else binned\n transit_mask = self._detect_transit_mask(white_light)\n \n if np.any(transit_mask) and np.sum(transit_mask) >= 3:\n # Apply gentle smoothing to in-transit data (lower SNR)\n in_transit_data = binned[transit_mask]\n n_transit = len(in_transit_data)\n if n_transit >= 3:\n window_len = min(5, n_transit//2*2+1) # Ensure odd\n try:\n smoothed = savgol_filter(in_transit_data, window_length=window_len, \n polyorder=min(2, window_len//2), axis=0)\n # Conservative blend: 80% original, 20% smoothed\n binned[transit_mask] = 0.8 * in_transit_data + 0.2 * smoothed\n except:\n pass\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 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 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 binned *= w[None, :]\n return binned\n\n def _process_planet_sensor(self, args):\n planet_id, sensor, dataset = args['planet_id'], args['sensor'], args['dataset']\n calibrated = self._calibrate_single_signal(planet_id, sensor, dataset)\n preprocessed = self._preprocess_calibrated_signal(calibrated, sensor)\n return preprocessed\n\n def process_all_data(self, dataset: str, planet_ids: np.ndarray):\n args_fgs1 = [dict(planet_id=pid, sensor=\"FGS1\", dataset=dataset) for pid in planet_ids]\n preprocessed_fgs1 = pqdm(args_fgs1, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n args_airs_ch0 = [dict(planet_id=pid, sensor=\"AIRS-CH0\", dataset=dataset) for pid in planet_ids]\n preprocessed_airs_ch0 = pqdm(args_airs_ch0, self._process_planet_sensor, n_jobs=self.cfg.N_JOBS)\n preprocessed_signal = np.concatenate([np.stack(preprocessed_fgs1), np.stack(preprocessed_airs_ch0)], axis=2)\n return preprocessed_signal\n\n# ------------------------------\n# Transit model (user's)\n# ------------------------------\n\nclass TransitModel:\n def __init__(self, config: 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 signal1 = signal[:min_index]\n signal2 = signal[min_index:]\n grad1 = np.gradient(signal1)\n grad1 /= max(grad1.max(), 1e-12)\n grad2 = np.gradient(signal2)\n grad2 /= max(grad2.max(), 1e-12)\n phase1 = int(np.argmin(grad1))\n phase2 = int(np.argmax(grad2)) + min_index\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 if phase1 - delta <= 0 or phase2 + delta >= len(signal) or phase2 - delta - (phase1 + delta) < 5:\n delta = 2\n y = np.concatenate([\n signal[: phase1 - delta],\n signal[phase1 + delta : phase2 - delta] * (1 + s),\n signal[phase2 + delta :]\n ])\n x = np.arange(len(y))\n coeffs = np.polyfit(x, y, deg=power)\n poly = np.poly1d(coeffs)\n error = np.abs(poly(x) - y).mean()\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 phase1, phase2 = self._phase_detector(signal_1d)\n phase1 = max(self.cfg.MODEL_OPTIMIZATION_DELTA, phase1)\n phase2 = min(len(signal_1d) - self.cfg.MODEL_OPTIMIZATION_DELTA - 1, phase2)\n result = minimize(\n fun=self._objective_function,\n x0=[0.0001],\n args=(signal_1d, phase1, phase2),\n method=\"Nelder-Mead\"\n )\n return result.x[0]\n\n def predict_all(self, preprocessed_signals):\n preds = [self.predict(single) for single in tqdm(preprocessed_signals)]\n return np.array(preds) * self.cfg.SCALE\n\n# ------------------------------\n# Feature extraction for ML layer\n# ------------------------------\n\ndef compute_ts_features(single: np.ndarray, cfg) -> dict:\n \"\"\"\n single: ndarray (n_bins, n_channels) - результат _preprocess_calibrated_signal для одного планета+сенсор-комбинации.\n Формат у тебя: column 0 == FGS, columns 1: == AIRS channels (может быть >=1).\n Возвращает словарь признаков (числа).\n \"\"\"\n # guard\n if single is None or single.size == 0:\n return {}\n\n # basic splits\n fgs = single[:, 0].astype(float).ravel()\n if single.shape[1] > 1:\n air_cols = single[:, 1:]\n air_mean_ts = np.nanmean(air_cols, axis=1).astype(float)\n else:\n air_mean_ts = np.array([], dtype=float)\n\n # helper safe funcs\n def safe_mean(a): \n return float(np.nanmean(a)) if a.size else 0.0\n def safe_std(a):\n return float(np.nanstd(a)) if a.size else 0.0\n def safe_mad(a):\n if a.size == 0: return 0.0\n m = np.nanmedian(a)\n return float(np.nanmedian(np.abs(a - m)))\n def safe_min(a):\n return float(np.nanmin(a)) if a.size else 0.0\n def safe_max(a):\n return float(np.nanmax(a)) if a.size else 0.0\n def safe_pct(a, q):\n if a.size == 0: return 0.0\n return float(np.nanpercentile(a, q))\n\n feats = {}\n\n # === simple statistics ===\n feats.update({\n \"fgs_mean\": safe_mean(fgs),\n \"fgs_std\": safe_std(fgs),\n \"fgs_mad\": safe_mad(fgs),\n \"fgs_min\": safe_min(fgs),\n \"fgs_max\": safe_max(fgs),\n \"fgs_ptp\": safe_max(fgs) - safe_min(fgs),\n \"fgs_skew\": float(skew(fgs)) if fgs.size > 2 else 0.0,\n \"fgs_kurt\": float(kurtosis(fgs)) if fgs.size > 4 else 0.0,\n \"fgs_p25\": safe_pct(fgs, 25),\n \"fgs_p50\": safe_pct(fgs, 50),\n \"fgs_p75\": safe_pct(fgs, 75),\n \"fgs_rms\": float(np.sqrt(np.nanmean((fgs - np.nanmean(fgs))**2))) if fgs.size else 0.0,\n })\n\n # === features from AIRS mean timeseries (if present) ===\n if air_mean_ts.size:\n feats.update({\n \"air_mean\": safe_mean(air_mean_ts),\n \"air_std\": safe_std(air_mean_ts),\n \"air_mad\": safe_mad(air_mean_ts),\n \"air_p50\": safe_pct(air_mean_ts, 50),\n \"air_p75\": safe_pct(air_mean_ts, 75),\n \"air_rms\": float(np.sqrt(np.nanmean((air_mean_ts - np.nanmean(air_mean_ts))**2))),\n \"air_skew\": float(skew(air_mean_ts)) if air_mean_ts.size > 2 else 0.0,\n \"air_kurt\": float(kurtosis(air_mean_ts)) if air_mean_ts.size > 4 else 0.0,\n })\n else:\n # put placeholders to keep feature matrix consistent\n feats.update({k: 0.0 for k in [\"air_mean\",\"air_std\",\"air_mad\",\"air_p50\",\"air_p75\",\"air_rms\",\"air_skew\",\"air_kurt\"]})\n\n # === ingress/egress & symmetry features ===\n # approximate transit phase detection like у тебя: min in search window\n try:\n # localized min on fgs -> pmin\n pmin = int(np.nanargmin(fgs)) if fgs.size else 0\n left = fgs[:pmin] if pmin>0 else fgs\n right = fgs[pmin+1:] if pmin+1= width else arr\n x = np.arange(len(seg))\n A = np.vstack([x, np.ones_like(x)]).T\n m, c = np.linalg.lstsq(A, seg, rcond=None)[0]\n return float(m)\n feats[\"fgs_ingress_slope\"] = slope_of_segment(left[::-1], width=8) # last part before min\n feats[\"fgs_egress_slope\"] = slope_of_segment(right, width=8) # first part after min\n feats[\"fgs_symmetry\"] = float(np.nanmean(left)) - float(np.nanmean(right)) if left.size and right.size else 0.0\n except Exception:\n feats.update({\"fgs_ingress_slope\":0.0,\"fgs_egress_slope\":0.0,\"fgs_symmetry\":0.0})\n\n # === autocorrelation / lag features ===\n def autocorr(a, lag=1):\n if a.size <= lag: return 0.0\n a = a - np.nanmean(a)\n denom = np.nansum((a)**2)\n if denom == 0: return 0.0\n return float(np.nansum(a[:-lag] * a[lag:]) / denom)\n feats[\"fgs_acf1\"] = autocorr(fgs, lag=1)\n feats[\"fgs_acf2\"] = autocorr(fgs, lag=2)\n feats[\"air_acf1\"] = autocorr(air_mean_ts, lag=1) if air_mean_ts.size else 0.0\n\n # === fft / spectral features ===\n def spectral_feats(ts):\n if ts.size < 4:\n return {\"spec_low_energy\":0.0, \"spec_high_energy\":0.0, \"spec_centroid\":0.0}\n x = ts - np.nanmean(ts)\n fft = np.abs(np.fft.rfft(x))\n freqs = np.fft.rfftfreq(ts.size)\n # energies:\n K = max(1, int(len(fft)*0.1))\n low_energy = float(np.sum(fft[:K]))\n high_energy = float(np.sum(fft[-K:]))\n # centroid\n denom = np.sum(fft) if np.sum(fft) > 0 else 1.0\n centroid = float(np.sum(freqs * fft) / denom)\n return {\"spec_low_energy\":low_energy, \"spec_high_energy\":high_energy, \"spec_centroid\":centroid}\n\n sf_fgs = spectral_feats(fgs)\n feats.update({f\"fgs_{k}\": v for k,v in sf_fgs.items()})\n sf_air = spectral_feats(air_mean_ts) if air_mean_ts.size else {\"spec_low_energy\":0.0,\"spec_high_energy\":0.0,\"spec_centroid\":0.0}\n feats.update({f\"air_{k}\": v for k,v in sf_air.items()})\n\n # === zero-crossings / number of peaks ===\n def zero_crossings(a):\n if a.size < 2: return 0\n s = np.sign(a - np.nanmedian(a))\n s[s==0] = 1\n return int(np.sum(np.abs(np.diff(s)))/2)\n feats[\"fgs_zero_cross\"] = zero_crossings(fgs)\n feats[\"air_zero_cross\"] = zero_crossings(air_mean_ts) if air_mean_ts.size else 0\n\n # === area/energy in 'inn' vs 'oot' using basic transit split ===\n try:\n # reuse pmin approximate split:\n inn = fgs[max(0, pmin-10): min(len(fgs), pmin+10)] if fgs.size else fgs\n oot = np.concatenate([fgs[:max(0,pmin-10)], fgs[min(len(fgs), pmin+10):]]) if fgs.size else fgs\n feats[\"fgs_inn_mean\"] = safe_mean(inn)\n feats[\"fgs_oot_mean\"] = safe_mean(oot)\n feats[\"fgs_depth_est\"] = float(safe_mean(oot) - safe_mean(inn))\n feats[\"fgs_inn_area\"] = float(np.nansum(inn))\n except Exception:\n feats.update({\"fgs_inn_mean\":0.0,\"fgs_oot_mean\":0.0,\"fgs_depth_est\":0.0,\"fgs_inn_area\":0.0})\n\n # keep only numeric scalars\n for k,v in list(feats.items()):\n if np.size(v) > 1 or (isinstance(v, (list,tuple))):\n feats[k] = float(np.nanmean(v)) if np.size(v) else 0.0\n\n return feats\n\n\ndef extract_features(preprocessed: np.ndarray, cfg: Config, star_df: pd.DataFrame) -> pd.DataFrame:\n \"\"\"\n New extract_features: вызывает compute_ts_features для каждого preprocessed[i]\n и соединяет с мета-признаками из star_df (например Rs, Ms, Ts...).\n Возвращает DataFrame с индексом == star_df.index (planet_id).\n \"\"\"\n feats = []\n for i, single in enumerate(preprocessed):\n f = compute_ts_features(single, cfg)\n # existing lightweight original features (keep them too)\n try:\n # original simple features (if you want to preserve old ones)\n fgs = single[:,0]\n air = single[:,1:] if single.shape[1] > 1 else np.empty((single.shape[0],0))\n air_white = np.nanmean(air, axis=1) if air.size else np.array([])\n # depth and width like old extractor\n p1, p2 = _phase_detector_signal(savgol_filter(np.nanmean(single[:,1:], axis=1) if single.shape[1]>1 else fgs, 20, 2), cfg)\n p1 = max(cfg.MODEL_OPTIMIZATION_DELTA, p1)\n p2 = min(len(fgs)-cfg.MODEL_OPTIMIZATION_DELTA-1, p2)\n inn = fgs[p1:p2] if p2>p1 else fgs\n oot = np.concatenate([fgs[:p1], fgs[p2:]]) if p2>p1 else fgs\n f.update({\n \"fgs_in_mean\": float(np.nanmean(inn)) if inn.size else float(np.nanmean(fgs)),\n \"fgs_in_std\": float(np.nanstd(inn)) if inn.size else float(np.nanstd(fgs)),\n \"fgs_oot_mean\": float(np.nanmean(oot)),\n \"fgs_oot_std\": float(np.nanstd(oot)),\n \"fgs_depth_analytic\": float(np.nanmean(oot) - np.nanmean(inn)) if inn.size else 0.0,\n \"fgs_width\": float(max(3, p2-p1)) if p2>p1 else 3.0,\n \"air_mean\": float(np.nanmean(air_white)) if air_white.size else 0.0,\n \"air_std\": float(np.nanstd(air_white)) if air_white.size else 0.0\n })\n except Exception:\n pass\n feats.append(f)\n\n X = pd.DataFrame(feats, index=star_df.index)\n # join with selected meta columns (numeric only)\n meta_cols = [c for c in [\"Rs\",\"Ms\",\"Ts\",\"Mp\",\"e\",\"P\",\"sma\",\"i\"] if c in star_df.columns]\n if meta_cols:\n meta = star_df[meta_cols].copy()\n X = X.join(meta)\n \n # optional: fill NaNs, small imputation\n X = X.fillna(0.0)\n return X\n\n\ndef add_multiscale_features(X: pd.DataFrame, preprocessed_signals: list, cfg) -> pd.DataFrame:\n \"\"\"Add features at multiple time scales using different smoothing windows\"\"\"\n X = X.copy()\n scales = [3, 7, 15, 31] # Different smoothing windows\n \n for i, single in enumerate(preprocessed_signals):\n if i >= len(X):\n break\n fgs = single[:, 0] # FGS signal\n \n for scale in scales:\n if len(fgs) > scale:\n try:\n smooth = savgol_filter(fgs, min(scale, len(fgs)//2*2+1), 2)\n detail = fgs - smooth\n \n X.loc[X.index[i], f'scale_{scale}_energy'] = np.var(detail)\n X.loc[X.index[i], f'scale_{scale}_trend'] = np.mean(np.gradient(smooth))\n except:\n # If smoothing fails, set to default values\n X.loc[X.index[i], f'scale_{scale}_energy'] = 0.0\n X.loc[X.index[i], f'scale_{scale}_trend'] = 0.0\n \n return X.fillna(0.0)\n\n\ndef add_transit_quality_features(X: pd.DataFrame, preprocessed_signals: list, cfg) -> pd.DataFrame:\n \"\"\"Add features that measure transit signal quality\"\"\"\n X = X.copy()\n \n for i, single in enumerate(preprocessed_signals):\n if i >= len(X):\n break\n fgs = single[:, 0] # FGS signal\n \n try:\n # Find transit phases using simple approach\n smooth_fgs = savgol_filter(fgs, min(20, len(fgs)//3*2+1), 2)\n \n # Simple transit detection\n baseline = np.median(smooth_fgs)\n threshold = baseline - 1.5 * np.std(smooth_fgs)\n in_transit = smooth_fgs < threshold\n \n if np.any(in_transit):\n # Find transit boundaries\n transit_indices = np.where(in_transit)[0]\n p1, p2 = transit_indices[0], transit_indices[-1]\n \n if p2 > p1 and p1 > 5 and len(fgs) - p2 > 5:\n transit_region = fgs[p1:p2]\n oot_left = fgs[:p1]\n oot_right = fgs[p2:]\n \n # Transit signal-to-noise ratio\n transit_mean = np.nanmean(transit_region)\n oot_mean = np.nanmean(np.concatenate([oot_left, oot_right]))\n oot_std = np.nanstd(np.concatenate([oot_left, oot_right]))\n \n if oot_std > 0:\n X.loc[X.index[i], 'transit_snr'] = (oot_mean - transit_mean) / oot_std\n \n # Transit consistency (inverse of scatter in transit)\n transit_std = np.nanstd(transit_region)\n if transit_std > 0:\n X.loc[X.index[i], 'transit_consistency'] = 1.0 / transit_std\n \n # V-shape quality (ingress vs egress symmetry)\n if len(fgs) > 10:\n gradients = np.gradient(fgs)\n ingress_grad = np.mean(gradients[max(0, p1-3):p1+3])\n egress_grad = np.mean(gradients[p2-3:min(len(gradients), p2+3)])\n X.loc[X.index[i], 'vshape_quality'] = -(ingress_grad - egress_grad) / 2.0\n \n except:\n # If any processing fails, use default values\n X.loc[X.index[i], 'transit_snr'] = 0.0\n X.loc[X.index[i], 'transit_consistency'] = 0.0\n X.loc[X.index[i], 'vshape_quality'] = 0.0\n \n return X.fillna(0.0)\n\n\ndef add_atmospheric_physics_features(X: pd.DataFrame, preprocessed_signals: list, cfg) -> pd.DataFrame:\n \"\"\"Add physics-based features using actual wavelength information from Ariel instruments\"\"\"\n X_new = X.copy()\n \n if not hasattr(cfg, 'wavelengths') or cfg.wavelengths is None:\n print(\"[WARN] No wavelength information available - skipping atmospheric physics features\")\n return X_new\n \n # AIRS-CH0 wavelengths (exclude FGS channel)\n wl_values = cfg.wavelengths # These are the actual wavelengths in micrometers\n \n for i, single in enumerate(preprocessed_signals):\n if i >= len(X) or single.shape[1] <= 1:\n # No AIRS channels - set default values\n for feat in ['rayleigh_slope', 'rayleigh_deviation', 'molecular_absorption_strength',\n 'water_2_7um', 'water_3_2um', 'co2_2_0um', 'co_2_3um', 'spectral_curvature',\n 'blue_slope', 'red_slope', 'transmission_range', 'spectral_roughness']:\n X_new.loc[X_new.index[i], feat] = 0.0\n continue\n \n try:\n airs_data = single[:, 1:] # AIRS spectroscopic channels\n n_wl = airs_data.shape[1]\n \n # Match wavelengths to AIRS channels (may need alignment)\n if len(wl_values) >= n_wl:\n wl_subset = wl_values[:n_wl] # Use first n_wl wavelengths\n else:\n # Interpolate if needed\n wl_subset = np.linspace(wl_values[0], wl_values[-1], n_wl)\n \n # Compute mean transmission spectrum across time\n mean_spectrum = np.nanmean(airs_data, axis=0)\n \n if len(mean_spectrum) >= 5 and len(wl_subset) >= 5:\n # 1. Rayleigh Scattering Analysis (λ^-4 dependence)\n try:\n # Fit power law in log-log space: spectrum ∝ λ^β\n log_wl = np.log(wl_subset)\n log_spec = np.log(np.abs(mean_spectrum) + 1e-12)\n \n # Remove any inf/nan values\n mask = np.isfinite(log_wl) & np.isfinite(log_spec)\n if np.sum(mask) >= 3:\n slope, intercept = np.polyfit(log_wl[mask], log_spec[mask], 1)\n \n X_new.loc[X_new.index[i], 'rayleigh_slope'] = slope\n # Theoretical Rayleigh scattering has slope ≈ -4\n X_new.loc[X_new.index[i], 'rayleigh_deviation'] = abs(slope + 4.0)\n else:\n X_new.loc[X_new.index[i], 'rayleigh_slope'] = 0.0\n X_new.loc[X_new.index[i], 'rayleigh_deviation'] = 0.0\n \n except:\n X_new.loc[X_new.index[i], 'rayleigh_slope'] = 0.0\n X_new.loc[X_new.index[i], 'rayleigh_deviation'] = 0.0\n \n # 2. Molecular Absorption Feature Detection\n try:\n # Smooth spectrum to find broad features\n if len(mean_spectrum) >= 5:\n smooth_spec = savgol_filter(mean_spectrum, min(5, len(mean_spectrum)//3*2+1), 2)\n \n # Find absorption features (local minima)\n peaks, properties = find_peaks(-smooth_spec, prominence=0.1*np.std(smooth_spec))\n \n X_new.loc[X_new.index[i], 'molecular_absorption_strength'] = len(peaks)\n \n # Spectral roughness (high-frequency variations)\n if len(mean_spectrum) > 1:\n roughness = np.std(np.diff(mean_spectrum))\n X_new.loc[X_new.index[i], 'spectral_roughness'] = roughness\n else:\n X_new.loc[X_new.index[i], 'spectral_roughness'] = 0.0\n else:\n X_new.loc[X_new.index[i], 'molecular_absorption_strength'] = 0.0\n X_new.loc[X_new.index[i], 'spectral_roughness'] = 0.0\n \n except:\n X_new.loc[X_new.index[i], 'molecular_absorption_strength'] = 0.0\n X_new.loc[X_new.index[i], 'spectral_roughness'] = 0.0\n \n # 3. Specific Molecular Band Analysis\n # Key atmospheric absorption bands in AIRS range\n molecular_bands = {\n 'water_2_7um': 2.7, # Water vapor\n 'water_3_2um': 3.2, # Water vapor\n 'co2_2_0um': 2.0, # CO2\n 'co_2_3um': 2.3 # CO (carbon monoxide)\n }\n \n for band_name, band_center in molecular_bands.items():\n try:\n if np.min(wl_subset) <= band_center <= np.max(wl_subset):\n # Find closest wavelength index\n idx = np.argmin(np.abs(wl_subset - band_center))\n \n if 2 <= idx < len(mean_spectrum) - 2:\n # Local absorption depth relative to continuum\n local_continuum = np.mean([mean_spectrum[idx-2], mean_spectrum[idx+2]])\n absorption_depth = local_continuum - mean_spectrum[idx]\n X_new.loc[X_new.index[i], band_name] = absorption_depth\n else:\n X_new.loc[X_new.index[i], band_name] = 0.0\n else:\n X_new.loc[X_new.index[i], band_name] = 0.0\n except:\n X_new.loc[X_new.index[i], band_name] = 0.0\n \n # 4. Spectral Shape Analysis\n try:\n # Quadratic fit for overall curvature\n wl_norm = (wl_subset - np.mean(wl_subset)) / (np.std(wl_subset) + 1e-9)\n curvature = np.polyfit(wl_norm, mean_spectrum, 2)[0]\n X_new.loc[X_new.index[i], 'spectral_curvature'] = curvature\n \n # Blue vs Red slopes (short vs long wavelength behavior)\n mid_idx = len(wl_subset) // 2\n if mid_idx >= 2:\n blue_slope = np.polyfit(wl_subset[:mid_idx], mean_spectrum[:mid_idx], 1)[0]\n red_slope = np.polyfit(wl_subset[mid_idx:], mean_spectrum[mid_idx:], 1)[0]\n \n X_new.loc[X_new.index[i], 'blue_slope'] = blue_slope\n X_new.loc[X_new.index[i], 'red_slope'] = red_slope\n else:\n X_new.loc[X_new.index[i], 'blue_slope'] = 0.0\n X_new.loc[X_new.index[i], 'red_slope'] = 0.0\n \n # Transmission spectrum dynamic range\n transmission_range = np.max(mean_spectrum) - np.min(mean_spectrum)\n X_new.loc[X_new.index[i], 'transmission_range'] = transmission_range\n \n except:\n X_new.loc[X_new.index[i], 'spectral_curvature'] = 0.0\n X_new.loc[X_new.index[i], 'blue_slope'] = 0.0\n X_new.loc[X_new.index[i], 'red_slope'] = 0.0\n X_new.loc[X_new.index[i], 'transmission_range'] = 0.0\n \n else:\n # Insufficient data - set all features to 0\n for feat in ['rayleigh_slope', 'rayleigh_deviation', 'molecular_absorption_strength',\n 'water_2_7um', 'water_3_2um', 'co2_2_0um', 'co_2_3um', 'spectral_curvature',\n 'blue_slope', 'red_slope', 'transmission_range', 'spectral_roughness']:\n X_new.loc[X_new.index[i], feat] = 0.0\n \n except Exception as e:\n # Complete fallback - set all features to 0\n for feat in ['rayleigh_slope', 'rayleigh_deviation', 'molecular_absorption_strength',\n 'water_2_7um', 'water_3_2um', 'co2_2_0um', 'co_2_3um', 'spectral_curvature',\n 'blue_slope', 'red_slope', 'transmission_range', 'spectral_roughness']:\n X_new.loc[X_new.index[i], feat] = 0.0\n \n return X_new.fillna(0.0)\n\n\nimport numpy as np\n\ndef add_interaction_features(X: pd.DataFrame, eps: float = 1e-9, depth_col: str = \"fgs_depth_analytic\") -> pd.DataFrame:\n \"\"\"\n Добавляет набор производных признаков (произведения/отношения/логи) + физические признаки.\n X: DataFrame, должен содержать базовые признаки (fgs_*, air_*, fgs_depth_analytic, fgs_width)\n и звездные параметры Rs, Ms, Ts, Mp, e, P, sma, i (в градусах).\n depth_col: имя колонки с глубиной транзита (если есть) — используется для оценки Rp и планетных величин.\n Возвращает новый DataFrame с добавленными колонками.\n \"\"\"\n X = X.copy()\n\n # --- helper getters (safety) ---\n def g(name):\n return X[name].values if name in X.columns else np.zeros(len(X))\n\n n = len(X)\n\n # stellar / planet inputs (with safe defaults)\n Rs = np.maximum(g(\"Rs\").astype(float), eps) # R★ в R_sun\n Ms = np.maximum(g(\"Ms\").astype(float), eps) # M★ в M_sun\n Ts = np.nan_to_num(g(\"Ts\").astype(float), nan=5778.0) # T★ в К\n Mp = np.nan_to_num(g(\"Mp\").astype(float), nan=0.0) # M_p в M_jup\n e = np.nan_to_num(g(\"e\").astype(float), nan=0.0) # эксцентриситет\n P = np.maximum(np.nan_to_num(g(\"P\").astype(float), nan=1.0), eps) # период в днях\n sma = np.maximum(np.nan_to_num(g(\"sma\").astype(float), nan=1.0), eps) # a/Rs (безразмер)\n incl_deg = np.nan_to_num(g(\"i\").astype(float), nan=90.0) # градусы\n\n # signal features (fallback = 0)\n fgs_depth = np.nan_to_num(g(\"fgs_depth_analytic\").astype(float) if \"fgs_depth_analytic\" in X.columns else g(\"fgs_depth\"), nan=0.0)\n fgs_width = np.nan_to_num(g(\"fgs_width\").astype(float) if \"fgs_width\" in X.columns else 0.0, nan=0.0)\n fgs_mean = np.nan_to_num(g(\"fgs_mean\").astype(float), nan=0.0)\n fgs_std = np.nan_to_num(g(\"fgs_std\").astype(float), nan=0.0)\n fgs_rms = np.nan_to_num(g(\"fgs_rms\").astype(float) if \"fgs_rms\" in X.columns else 0.0, nan=0.0)\n fgs_skew = np.nan_to_num(g(\"fgs_skew\").astype(float) if \"fgs_skew\" in X.columns else 0.0, nan=0.0)\n fgs_kurt = np.nan_to_num(g(\"fgs_kurt\").astype(float) if \"fgs_kurt\" in X.columns else 0.0, nan=0.0)\n\n air_mean = np.nan_to_num(g(\"air_mean\").astype(float) if \"air_mean\" in X.columns else 0.0, nan=0.0)\n air_std = np.nan_to_num(g(\"air_std\").astype(float) if \"air_std\" in X.columns else 0.0, nan=0.0)\n air_rms = np.nan_to_num(g(\"air_rms\").astype(float) if \"air_rms\" in X.columns else 0.0, nan=0.0)\n fgs_spec_centroid = np.nan_to_num(g(\"fgs_spec_centroid\").astype(float) if \"fgs_spec_centroid\" in X.columns else 0.0, nan=0.0)\n air_spec_centroid = np.nan_to_num(g(\"air_spec_centroid\").astype(float) if \"air_spec_centroid\" in X.columns else 0.0, nan=0.0)\n\n # trig\n incl_rad = np.deg2rad(incl_deg)\n sin_i = np.sin(incl_rad)\n cos_i = np.cos(incl_rad)\n\n # stellar density proxy\n rho = Ms / (Rs**3 + eps) # Ms / Rs^3 (solar units)\n\n # ========== Вставляем уже существующие interaction-признаки ==========\n X[\"depth_x_Rs\"] = fgs_depth * Rs\n X[\"depth_div_Rs\"] = fgs_depth / (Rs + eps)\n\n X[\"depth_x_sin_i\"] = fgs_depth * sin_i\n X[\"depth_x_cos_i\"] = fgs_depth * cos_i\n\n X[\"depth_x_Mp\"] = fgs_depth * Mp\n X[\"depth_x_log1pMp\"] = fgs_depth * np.log1p(Mp + eps)\n\n X[\"width_x_rho\"] = fgs_width * rho\n X[\"width_div_sma\"] = fgs_width / (sma + eps)\n X[\"width_x_inv_sma2\"] = fgs_width / (sma**2 + eps)\n\n X[\"stellar_rho\"] = rho\n X[\"stellar_rho_sqrt\"] = np.sqrt(np.maximum(rho, 0.0))\n X[\"stellar_rho_cbrt\"] = np.cbrt(np.maximum(rho, 0.0))\n\n X[\"width_x_1pluse\"] = fgs_width * (1.0 + e)\n X[\"depth_x_1pluse\"] = fgs_depth * (1.0 + e)\n\n X[\"logP\"] = np.log(P + eps)\n #X[\"logsma\"] = np.log(sma + eps)\n X[\"width_x_logP\"] = fgs_width * X[\"logP\"]\n X[\"depth_x_logsma\"] = fgs_depth * np.log(sma + eps)\n\n X[\"fgs_spec_x_Ts\"] = fgs_spec_centroid * Ts\n X[\"air_spec_x_Ts\"] = air_spec_centroid * Ts\n\n X[\"fgsmean_x_airmean\"] = fgs_mean * air_mean\n X[\"fgsrms_x_airrms\"] = fgs_rms * air_rms\n X[\"fgs_spec_x_air_spec\"] = fgs_spec_centroid * air_spec_centroid\n\n X[\"fgs_mean_x_std\"] = fgs_mean * fgs_std\n X[\"fgs_depth_x_rms\"] = fgs_depth * fgs_rms\n X[\"fgs_skew_x_kurt\"] = fgs_skew * fgs_kurt\n\n X[\"Rs_div_Ms\"] = Rs / (Ms + eps)\n X[\"Ms_div_Rs3\"] = Ms / (Rs**3 + eps)\n\n X[\"logRs\"] = np.log(Rs + eps)\n X[\"logMs\"] = np.log(Ms + eps)\n X[\"logMp\"] = np.log1p(Mp + eps)\n\n X[\"depth_x_Rs_sq\"] = fgs_depth * (Rs**2)\n X[\"width_x_sqrtMs\"] = fgs_width * np.sqrt(np.maximum(Ms, 0.0))\n\n # ========== Новые физически-мотивированные признаки ==========\n # константы (SI)\n G = 6.67430e-11 # m^3 kg^-1 s^-2\n M_sun = 1.98847e30 # kg\n R_sun = 6.957e8 # m\n M_jup = 1.89813e27 # kg\n sec_per_day = 86400.0\n T_sun = 5778.0\n\n # простые логи / тригонометрия\n X[\"sin_i\"] = sin_i\n X[\"cos_i\"] = cos_i\n X[\"log_P_days\"] = np.log(P + eps)\n X[\"log_sma\"] = np.log(sma + eps)\n\n # impact parameter (безразмерный)\n X[\"impact_b\"] = sma * cos_i\n\n # approximate transit duration (days) (упрощённая оценка)\n X[\"T_dur_days_approx\"] = P / (np.pi * sma + eps)\n\n # mean motion (1/day)\n X[\"mean_motion_1pd\"] = 2.0 * np.pi / (P + eps)\n\n # Kepler-based stellar mass estimate (переводим в SI)\n a_m = sma * Rs * R_sun # полусось в метрах\n P_s = P * sec_per_day # период в секундах\n with np.errstate(divide='ignore', invalid='ignore', over='ignore'):\n M_est_SI = (4.0 * np.pi**2) * (a_m**3) / (G * (P_s**2) + eps)\n M_est_Msun = M_est_SI / M_sun\n X[\"M_kepler_est_Msun\"] = np.nan_to_num(M_est_Msun, nan=0.0)\n X[\"M_kepler_ratio\"] = X[\"M_kepler_est_Msun\"] / (Ms + eps)\n\n # equilibrium temperature proxy and insolation proxy (без абсолютных констант)\n X[\"Teq_rel\"] = Ts * np.sqrt(1.0 / (2.0 * sma + eps))\n X[\"insolation_rel\"] = (Ts / T_sun)**4 * (Rs**2) / (sma**2 + eps)\n X[\"log_insolation\"] = np.log(X[\"insolation_rel\"] + 1e-12)\n\n # orbital velocity (m/s -> km/s)\n with np.errstate(divide='ignore', invalid='ignore'):\n v_orb = 2.0 * np.pi * a_m / (P_s + eps) # m/s\n X[\"v_orb_km_s\"] = np.nan_to_num(v_orb / 1000.0, nan=0.0)\n X[\"log_v_orb\"] = np.log(np.abs(X[\"v_orb_km_s\"]) + 1e-12)\n\n # Планетные признаки — если есть глубина транзита (depth_col или fgs_depth_analytic)\n depth_arr = None\n if depth_col is not None and depth_col in X.columns:\n depth_arr = np.nan_to_num(X[depth_col].astype(float).values, nan=0.0)\n elif \"fgs_depth_analytic\" in X.columns:\n depth_arr = np.nan_to_num(X[\"fgs_depth_analytic\"].astype(float).values, nan=0.0)\n else:\n depth_arr = None\n\n if depth_arr is not None:\n # Rp / Rs и Rp (в R_sun)\n Rp_over_Rs = np.sqrt(np.clip(depth_arr, 0.0, None))\n Rp_Rsun = Rp_over_Rs * Rs\n X[\"Rp_over_Rs\"] = Rp_over_Rs\n X[\"Rp_Rsun\"] = Rp_Rsun\n\n # абсолютные величины в SI\n Rp_m = Rp_Rsun * R_sun\n Mp_kg = Mp * M_jup\n\n # surface gravity planet (m/s^2) и плотность (kg/m^3)\n with np.errstate(divide='ignore', invalid='ignore'):\n g_p = (G * Mp_kg) / (Rp_m**2 + eps)\n vol = (4.0/3.0) * np.pi * (Rp_m**3 + eps)\n rho_p = Mp_kg / vol\n\n X[\"planet_g_m_s2\"] = np.nan_to_num(g_p, nan=0.0)\n X[\"planet_rho_kg_m3\"] = np.nan_to_num(rho_p, nan=0.0)\n\n # Hill radius (в радиусах Солнца -> затем в единицах R_sun, можно интерпретировать относительно Rs)\n with np.errstate(divide='ignore', invalid='ignore'):\n R_H_m = a_m * ((Mp_kg) / (3.0 * Ms * M_sun + eps))**(1.0/3.0)\n X[\"R_H_Rs\"] = np.nan_to_num(R_H_m / R_sun, nan=0.0)\n\n # detectability proxy\n X[\"depth_div_sma2\"] = (Rp_over_Rs**2) / (sma**2 + eps)\n else:\n # заполняем нулями, чтобы колонки были на месте\n X[\"Rp_over_Rs\"] = 0.0\n X[\"Rp_Rsun\"] = 0.0\n X[\"planet_g_m_s2\"] = 0.0\n X[\"planet_rho_kg_m3\"] = 0.0\n X[\"R_H_Rs\"] = 0.0\n X[\"depth_div_sma2\"] = 0.0\n\n # final cleaning: replace inf/nan and fill\n X = X.replace([np.inf, -np.inf], np.nan).fillna(0.0)\n\n return X\n\n\ndef add_spectral_correlation_features(X: pd.DataFrame, preprocessed_signals: list, cfg) -> pd.DataFrame:\n \"\"\"\n Add spectral correlation features that capture wavelength-to-wavelength correlations\n and transmission spectrum characteristics for atmospheric spectroscopy.\n \"\"\"\n X_new = X.copy()\n \n # Helper function to safely extract transit phases\n def safe_phase_extraction(fgs_signal, cfg):\n try:\n # Use your existing transit model phase detection\n model = TransitModel(cfg)\n p1, p2 = model._phase_detector(fgs_signal)\n return max(0, p1), min(len(fgs_signal), p2)\n except:\n # Fallback to simple detection\n n = len(fgs_signal)\n return n//3, 2*n//3\n \n for i, single in enumerate(preprocessed_signals):\n if i >= len(X) or single.shape[1] <= 1:\n # No AIRS channels available\n for feat in ['transmission_slope', 'transmission_intercept', 'transmission_curvature', \n 'spectral_roughness', 'spectral_coherence', 'first_channel_coherence',\n 'transit_spectral_depth', 'spectral_snr']:\n X_new.loc[X_new.index[i], feat] = 0.0\n continue\n \n fgs = single[:, 0]\n airs_data = single[:, 1:] # AIRS channels\n \n try:\n # 1. Transit phase detection\n p1, p2 = safe_phase_extraction(fgs, cfg)\n \n if p2 > p1 + 3 and airs_data.shape[1] >= 3:\n # 2. Compute transmission spectrum (relative depths per wavelength)\n in_transit = np.nanmean(airs_data[p1:p2], axis=0)\n \n # Out-of-transit: combine before and after transit\n oot_data = []\n if p1 > 2:\n oot_data.append(airs_data[:p1])\n if p2 < len(airs_data) - 2:\n oot_data.append(airs_data[p2:])\n \n if oot_data:\n out_transit = np.nanmean(np.concatenate(oot_data), axis=0)\n else:\n out_transit = np.nanmean(airs_data, axis=0)\n \n # Transmission spectrum: (out - in) / out (transit depth per wavelength)\n transmission = (out_transit - in_transit) / (out_transit + 1e-9)\n transmission = np.nan_to_num(transmission, 0.0)\n \n # 3. Spectral shape analysis\n if len(transmission) >= 3:\n wl_indices = np.arange(len(transmission))\n \n # Linear trend across wavelengths (atmospheric slope)\n try:\n slope, intercept = np.polyfit(wl_indices, transmission, 1)\n X_new.loc[X_new.index[i], 'transmission_slope'] = slope\n X_new.loc[X_new.index[i], 'transmission_intercept'] = intercept\n except:\n X_new.loc[X_new.index[i], 'transmission_slope'] = 0.0\n X_new.loc[X_new.index[i], 'transmission_intercept'] = 0.0\n \n # Curvature (2nd order polynomial - molecular features)\n if len(transmission) >= 5:\n try:\n curvature = np.polyfit(wl_indices, transmission, 2)[0]\n X_new.loc[X_new.index[i], 'transmission_curvature'] = curvature\n except:\n X_new.loc[X_new.index[i], 'transmission_curvature'] = 0.0\n else:\n X_new.loc[X_new.index[i], 'transmission_curvature'] = 0.0\n \n # Spectral roughness (high-frequency variations)\n if len(transmission) > 1:\n roughness = np.std(np.diff(transmission))\n X_new.loc[X_new.index[i], 'spectral_roughness'] = roughness\n else:\n X_new.loc[X_new.index[i], 'spectral_roughness'] = 0.0\n \n # Mean transit depth across all wavelengths\n mean_depth = np.nanmean(transmission)\n X_new.loc[X_new.index[i], 'transit_spectral_depth'] = mean_depth\n \n # Spectral SNR proxy\n if np.std(transmission) > 0:\n spectral_snr = mean_depth / (np.std(transmission) + 1e-9)\n X_new.loc[X_new.index[i], 'spectral_snr'] = spectral_snr\n else:\n X_new.loc[X_new.index[i], 'spectral_snr'] = 0.0\n \n # 4. Wavelength correlation analysis\n if airs_data.shape[1] >= 3:\n try:\n # Correlation matrix between wavelength channels\n corr_matrix = np.corrcoef(airs_data.T)\n \n # Mean correlation (spectral coherence)\n upper_tri = corr_matrix[np.triu_indices_from(corr_matrix, k=1)]\n mean_corr = np.nanmean(upper_tri)\n X_new.loc[X_new.index[i], 'spectral_coherence'] = mean_corr\n \n # Correlation with first channel (often cleanest reference)\n if corr_matrix.shape[0] > 1:\n first_ch_corr = np.nanmean(corr_matrix[0, 1:])\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = first_ch_corr\n else:\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = 0.0\n \n except:\n X_new.loc[X_new.index[i], 'spectral_coherence'] = 0.0\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = 0.0\n else:\n X_new.loc[X_new.index[i], 'spectral_coherence'] = 0.0\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = 0.0\n \n else:\n # Fallback values when transit detection fails\n for feat in ['transmission_slope', 'transmission_intercept', 'transmission_curvature', \n 'spectral_roughness', 'transit_spectral_depth', 'spectral_snr']:\n X_new.loc[X_new.index[i], feat] = 0.0\n \n # Still compute basic correlations if possible\n if airs_data.shape[1] >= 3:\n try:\n corr_matrix = np.corrcoef(airs_data.T)\n upper_tri = corr_matrix[np.triu_indices_from(corr_matrix, k=1)]\n mean_corr = np.nanmean(upper_tri)\n X_new.loc[X_new.index[i], 'spectral_coherence'] = mean_corr\n \n if corr_matrix.shape[0] > 1:\n first_ch_corr = np.nanmean(corr_matrix[0, 1:])\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = first_ch_corr\n else:\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = 0.0\n except:\n X_new.loc[X_new.index[i], 'spectral_coherence'] = 0.0\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = 0.0\n else:\n X_new.loc[X_new.index[i], 'spectral_coherence'] = 0.0\n X_new.loc[X_new.index[i], 'first_channel_coherence'] = 0.0\n \n except Exception as e:\n # Complete fallback - set all features to 0\n for feat in ['transmission_slope', 'transmission_intercept', 'transmission_curvature', \n 'spectral_roughness', 'spectral_coherence', 'first_channel_coherence',\n 'transit_spectral_depth', 'spectral_snr']:\n X_new.loc[X_new.index[i], feat] = 0.0\n \n return X_new.fillna(0.0)\n\n\n# ------------------------------\n# ML Trainer (safe blending around baseline)\n# ------------------------------\n\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.linear_model import Ridge\nfrom sklearn.decomposition import PCA\nfrom sklearn.model_selection import KFold\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nclass SafeMLCalibrator:\n \"\"\"\n Safe ML calibrator with feature scaling and robust handling of feature column changes.\n\n Usage:\n cal = SafeMLCalibrator(cfg, n_components=4)\n cal.fit(preproc_train, s_base_train, mu_train, X_train_df)\n mu_pred = cal.infer_mu(s_base_test, X_test_df, n_mu)\n sig_pred = cal.infer_sigma(X_test_df, base_sigma_scalar, n_mu)\n cal.feature_importance_summary(X_train_df)\n \"\"\"\n def __init__(self, cfg, n_components: int = 4,\n max_blend_scale: float = 1.0, max_blend_shape: float = 1.0, max_blend_sigma: float = 1.0):\n self.cfg = cfg\n self.K = n_components\n self.max_alpha = max_blend_scale\n self.max_beta = max_blend_shape\n self.max_gamma = max_blend_sigma\n\n # learned params\n self.a_ = None\n self.b_ = None\n self.alpha_ = 0.0\n\n self.pca_ = None\n self.ridge_shape_models_ = None\n self.beta_ = 0.0\n\n self.ridge_sigma_air_ = None\n self.ridge_sigma_fgs_ = None\n self.gamma_ = 0.55\n\n # feature scaler and names\n self.feature_scaler_ = None\n self.feature_names_ = None\n self.n_features_in_ = None\n\n # -------------------------\n # Fit\n # -------------------------\n def fit(self, preproc_train: np.ndarray, s_base_train: np.ndarray, mu_train: np.ndarray,\n star_feat_train: pd.DataFrame, random_state: int = 42, transit_quality: np.ndarray = None):\n \"\"\"\n Fit calibrator on training data, but choose blending params (alpha,beta) by minimizing\n Negative Log-Likelihood (NLL) on training set (OOF-style).\n \"\"\"\n # check inputs\n if not isinstance(star_feat_train, pd.DataFrame):\n raise ValueError(\"star_feat_train must be a pandas DataFrame with column names.\")\n self.feature_names_ = list(star_feat_train.columns)\n self.n_features_in_ = len(self.feature_names_)\n\n # --- 0) fit feature scaler\n self.feature_scaler_ = StandardScaler()\n X_all = star_feat_train[self.feature_names_].values.astype(float)\n self.feature_scaler_.fit(X_all)\n X_scaled = self.feature_scaler_.transform(X_all)\n\n N = s_base_train.shape[0]\n n_mu = mu_train.shape[1]\n \n # --- 0.5) Setup quality weights for training\n if transit_quality is not None:\n sample_weights = np.clip(transit_quality, 0.1, 1.0)\n sample_weights = sample_weights / np.mean(sample_weights) # Normalize\n print(f\"[INFO] Using quality weights: mean={np.mean(sample_weights):.3f}, \"\n f\"min={np.min(sample_weights):.3f}, max={np.max(sample_weights):.3f}\")\n else:\n sample_weights = np.ones(N)\n\n # --- 1) OOF linear calibration of scale (a + b * s) with quality weighting\n y_level = mu_train.mean(axis=1)\n kf = KFold(n_splits=5, shuffle=True, random_state=random_state)\n coefs = []\n for tr, va in kf.split(s_base_train):\n Xtr = s_base_train[tr].reshape(-1, 1)\n ytr = y_level[tr]\n \n # Apply quality weights for this fold\n weights_tr = transit_quality[tr] if transit_quality is not None else None\n \n mdl = Ridge(alpha=1e-6).fit(Xtr, ytr, sample_weight=weights_tr)\n coefs.append((float(mdl.intercept_), float(mdl.coef_[0])))\n self.a_ = float(np.mean([c[0] for c in coefs]))\n self.b_ = float(np.mean([c[1] for c in coefs]))\n\n # --- 2) Train sigma regressors on log(sigma) (they don't depend on alpha/beta)\n # compute residual-based sigma targets\n mu_base_train = np.tile(s_base_train.reshape(-1, 1), (1, n_mu))\n resid = mu_train - mu_base_train\n sigma_air_star = np.sqrt(np.mean(resid[:, 1:] ** 2, axis=1)) if n_mu > 1 else np.sqrt(np.mean(resid ** 2, axis=1))\n sigma_fgs_star = np.sqrt(resid[:, 0] ** 2) if n_mu > 0 else sigma_air_star\n\n # fit ridge regressors on X_scaled -> log(sigma)\n # clip targets to sensible bounds to avoid infinities\n sigma_air_clipped = np.clip(sigma_air_star, 1e-6, 0.1)\n sigma_fgs_clipped = np.clip(sigma_fgs_star, 1e-6, 0.1)\n self.ridge_sigma_air_ = Ridge(alpha=0.5).fit(X_scaled, np.log(sigma_air_clipped), \n sample_weight=transit_quality)\n self.ridge_sigma_fgs_ = Ridge(alpha=0.5).fit(X_scaled, np.log(sigma_fgs_clipped), \n sample_weight=transit_quality)\n\n # set gamma (keep conservative default or user-provided cap)\n self.gamma_ = min(0.55, self.max_gamma)\n\n # --- 3) PCA shape + Ridge on features (produce OOF normalized-shape predictions)\n level = mu_train.mean(axis=1, keepdims=True)\n norm = np.divide(mu_train, np.maximum(level, 1e-9))\n self.pca_ = PCA(n_components=self.K, random_state=random_state).fit(norm)\n Z = self.pca_.transform(norm)\n\n kf = KFold(n_splits=5, shuffle=True, random_state=random_state)\n oof_norm = np.zeros_like(norm)\n self.ridge_shape_models_ = []\n for tr, va in kf.split(X_scaled):\n # Apply quality weights for this fold\n weights_tr = transit_quality[tr] if transit_quality is not None else None\n \n mdl = Ridge(alpha=1.0).fit(X_scaled[tr], Z[tr], sample_weight=weights_tr)\n Z_hat = mdl.predict(X_scaled[va])\n # invert PCA to normalized mu-space\n oof_norm[va] = self.pca_.inverse_transform(Z_hat)\n self.ridge_shape_models_.append(mdl)\n\n # --- 4) Now choose (alpha, beta) by minimizing OOF Negative Log-Likelihood (NLL)\n # prepare sigma predictions on training set using ridge_sigma regressors\n sigma_air_ml_train = np.exp(self.ridge_sigma_air_.predict(X_scaled))\n sigma_fgs_ml_train = np.exp(self.ridge_sigma_fgs_.predict(X_scaled))\n # clip for numerical stability\n sigma_air_ml_train = np.clip(sigma_air_ml_train, 1e-6, 0.1)\n sigma_fgs_ml_train = np.clip(sigma_fgs_ml_train, 1e-6, 0.1)\n\n # combine with base sigma via gamma_ (geometric blend in log-space)\n base_sigma_scalar = self.cfg.SIGMA\n sigma_air_final_train = np.exp((1 - self.gamma_) * np.log(base_sigma_scalar) + self.gamma_ * np.log(sigma_air_ml_train))\n sigma_fgs_final_train = np.exp((1 - self.gamma_) * np.log(base_sigma_scalar) + self.gamma_ * np.log(sigma_fgs_ml_train))\n\n # Prepare constants for search\n flat = np.ones_like(oof_norm) # (N, n_mu)\n # ensure oof_norm has n_mu columns (pad or crop)\n if oof_norm.shape[1] != n_mu:\n if oof_norm.shape[1] < n_mu:\n pad = np.ones((oof_norm.shape[0], n_mu - oof_norm.shape[1]))\n oof_norm = np.concatenate([oof_norm, pad], axis=1)\n else:\n oof_norm = oof_norm[:, :n_mu]\n\n # grid of candidates (bounded by max_alpha/max_beta)\n alphas = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]\n betas = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]\n\n # restrict to allowed maxima\n alphas = [a for a in alphas if a <= self.max_alpha]\n betas = [b for b in betas if b <= self.max_beta]\n\n best_alpha, best_beta = 0.0, 0.0\n best_nll = np.inf\n\n # precompute log(2*pi) for efficiency\n log2pi = np.log(2.0 * np.pi)\n\n for alpha in alphas:\n # calibrated scale blend\n s_cal = self.a_ + self.b_ * s_base_train\n s_blend = (1 - alpha) * s_base_train + alpha * s_cal # shape (N,)\n\n for beta in betas:\n # shape blend (normalized)\n cand_shape = (1 - beta) * flat + beta * oof_norm # (N, n_mu)\n # mu prediction\n mu_pred = s_blend.reshape(-1, 1) * cand_shape # (N, n_mu)\n\n # sigma matrix for each mu-channel: first column = fgs, others = air\n sigmas = np.full((N, n_mu), base_sigma_scalar, dtype=float)\n sigmas[:, 0] = sigma_fgs_final_train\n if n_mu > 1:\n sigmas[:, 1:] = sigma_air_final_train.reshape(-1, 1)\n\n # numerical guard\n sigmas = np.clip(sigmas, 1e-8, 1.0)\n\n # compute NLL per entry: 0.5*(log(2*pi*sigma^2) + ((y - mu_pred)^2)/sigma^2)\n sq_err = (mu_train - mu_pred) ** 2\n denom = 2.0 * (sigmas ** 2)\n term = 0.5 * (np.log(denom) + (sq_err / denom) * 2.0) # rearranged to stable form\n # better compute robustly:\n # nll_matrix = 0.5 * (log2pi + 2*log(sigmas) + (sq_err / (sigmas**2)))\n nll_matrix = 0.5 * (log2pi + 2.0 * np.log(sigmas) + (sq_err / (sigmas ** 2)))\n\n # average NLL across all entries (stars × mu-channels)\n nll = np.nanmean(nll_matrix)\n\n if not np.isfinite(nll):\n continue\n if nll < best_nll:\n best_nll = nll\n best_alpha = alpha\n best_beta = beta\n\n # assign best found values (constrain to maxima)\n self.alpha_ = min(best_alpha, self.max_alpha)\n self.beta_ = min(best_beta, self.max_beta)\n\n # note: gamma_ left as previously set (could be tuned in future)\n # End of fit\n\n # -------------------------\n # Inference: mu\n # -------------------------\n def infer_mu(self, s_base: np.ndarray, star_feats: pd.DataFrame, n_mu: int) -> np.ndarray:\n \"\"\"\n Predict full mu matrix for new star_feats.\n - star_feats: DataFrame (N x F) — columns may differ; they will be reindexed to feature_names_ with fill_value=0.0\n - n_mu: number of mu channels to output\n \"\"\"\n if self.feature_scaler_ is None or self.feature_names_ is None:\n raise RuntimeError(\"SafeMLCalibrator not fitted. Call fit(...) before infer_mu.\")\n\n # align columns to training feature set (fill missing with zeros; ignore extra columns)\n star_feats_aligned = star_feats.reindex(columns=self.feature_names_, fill_value=0.0)\n X_query = self.feature_scaler_.transform(star_feats_aligned.values.astype(float))\n\n # scale blending\n s_cal = self.a_ + self.b_ * s_base\n s_blend = (1 - self.alpha_) * s_base + self.alpha_ * s_cal\n\n # shape prediction: average predictions of ridge_shape_models_\n Z_hat_list = [mdl.predict(X_query) for mdl in self.ridge_shape_models_]\n Z_hat = np.mean(Z_hat_list, axis=0)\n shape_ml = self.pca_.inverse_transform(Z_hat) # (N, k) -> (N, n_mu) maybe\n\n # safety: pad/crop to n_mu\n if shape_ml.shape[1] != n_mu:\n if shape_ml.shape[1] < n_mu:\n pad = np.ones((shape_ml.shape[0], n_mu - shape_ml.shape[1]))\n shape_ml = np.concatenate([shape_ml, pad], axis=1)\n else:\n shape_ml = shape_ml[:, :n_mu]\n\n shape_blend = (1 - self.beta_) * 1.0 + self.beta_ * shape_ml\n mu_ml = s_blend.reshape(-1, 1) * shape_blend\n mu_ml = np.clip(mu_ml, 0.0, None)\n return mu_ml\n\n # -------------------------\n # Inference: sigma\n # -------------------------\n def infer_sigma(self, star_feats: pd.DataFrame, base_sigma_scalar: float, n_mu: int):\n \"\"\"\n Predict sigma matrix for new star_feats.\n - base_sigma_scalar: fallback scalar (cfg.SIGMA)\n \"\"\"\n if self.feature_scaler_ is None or self.feature_names_ is None:\n raise RuntimeError(\"SafeMLCalibrator not fitted. Call fit(...) before infer_sigma.\")\n\n star_feats_aligned = star_feats.reindex(columns=self.feature_names_, fill_value=0.0)\n X_query = self.feature_scaler_.transform(star_feats_aligned.values.astype(float))\n\n sigma_air_ml = np.exp(self.ridge_sigma_air_.predict(X_query))\n sigma_fgs_ml = np.exp(self.ridge_sigma_fgs_.predict(X_query))\n\n sigma_air_final = np.exp((1 - self.gamma_) * np.log(base_sigma_scalar) + self.gamma_ * np.log(np.clip(sigma_air_ml, 1e-6, 0.1)))\n sigma_fgs_final = np.exp((1 - self.gamma_) * np.log(base_sigma_scalar) + self.gamma_ * np.log(np.clip(sigma_fgs_ml, 1e-6, 0.1)))\n\n sigmas = np.full((star_feats.shape[0], n_mu), base_sigma_scalar, dtype=float)\n sigmas[:, 0] = np.clip(sigma_fgs_final, 1e-6, 0.1)\n sigmas[:, 1:] = np.clip(sigma_air_final.reshape(-1, 1), 1e-6, 0.1)\n return sigmas\n\n # -------------------------\n # Feature importance / correlation summary\n # -------------------------\n def feature_importance_summary(self, star_feat_df: pd.DataFrame, top_n: int = 30,\n corr_top_k: int = 30, corr_threshold: float = 0.6,\n save_fig: bool = True, fig_name: str = \"feature_corr.png\"):\n \"\"\"\n Print and optionally save summary of feature importances and correlations.\n Returns (df_imp, corr).\n \"\"\"\n if self.feature_scaler_ is None:\n print(\"[WARN] feature_scaler_ not present — run fit() first.\")\n return None, None\n\n feat_names = list(star_feat_df.columns)\n n_feat = len(feat_names)\n X_scaled = self.feature_scaler_.transform(star_feat_df.reindex(columns=self.feature_names_, fill_value=0.0).values.astype(float))\n\n # shape importance from ridge_shape_models_\n if self.ridge_shape_models_ is None or len(self.ridge_shape_models_) == 0:\n shape_imp = np.zeros(n_feat)\n else:\n per_model = []\n for mdl in self.ridge_shape_models_:\n coef = getattr(mdl, \"coef_\", None)\n if coef is None:\n per_model.append(np.zeros(self.n_features_in_))\n continue\n coef = np.atleast_2d(coef)\n # mean abs across targets\n mean_abs = np.mean(np.abs(coef), axis=0)\n per_model.append(mean_abs)\n shape_imp = np.mean(np.vstack(per_model), axis=0)\n\n def _coef_abs_vector(mdl):\n if mdl is None:\n return np.zeros(self.n_features_in_)\n coef = getattr(mdl, \"coef_\", None)\n if coef is None:\n return np.zeros(self.n_features_in_)\n coef = np.asarray(coef)\n if coef.ndim == 2 and coef.shape[1] == self.n_features_in_:\n return np.mean(np.abs(coef), axis=0) if coef.shape[0] > 1 else np.abs(coef).ravel()\n elif coef.ndim == 1 and coef.shape[0] == self.n_features_in_:\n return np.abs(coef)\n else:\n c = np.abs(coef.ravel())\n if c.size < self.n_features_in_:\n c = np.pad(c, (0, self.n_features_in_ - c.size))\n return c[:self.n_features_in_]\n\n sigma_air_imp = _coef_abs_vector(getattr(self, \"ridge_sigma_air_\", None))\n sigma_fgs_imp = _coef_abs_vector(getattr(self, \"ridge_sigma_fgs_\", None))\n\n df_imp = pd.DataFrame({\n \"feature\": self.feature_names_,\n \"shape_imp\": shape_imp.astype(float),\n \"sigma_air_imp\": sigma_air_imp.astype(float),\n \"sigma_fgs_imp\": sigma_fgs_imp.astype(float),\n }).set_index(\"feature\")\n\n # normalize for comparability\n for c in df_imp.columns:\n mx = df_imp[c].max()\n df_imp[c] = df_imp[c] / (mx + 1e-12)\n df_imp[\"score_sum\"] = df_imp.sum(axis=1)\n df_imp = df_imp.sort_values(\"score_sum\", ascending=False)\n '''\n # print results\n print(\"\\n=== TOP feature importances (combined score) ===\")\n with pd.option_context('display.precision', 4, 'display.max_rows', top_n):\n print(df_imp.head(top_n))\n\n print(\"\\n=== TOP by shape importance ===\")\n print(df_imp.sort_values(\"shape_imp\", ascending=False).head(top_n).drop(columns=[\"score_sum\"]))\n print(\"\\n=== TOP by sigma_air importance ===\")\n print(df_imp.sort_values(\"sigma_air_imp\", ascending=False).head(top_n).drop(columns=[\"score_sum\"]))\n print(\"\\n=== TOP by sigma_fgs importance ===\")\n print(df_imp.sort_values(\"sigma_fgs_imp\", ascending=False).head(top_n).drop(columns=[\"score_sum\"]))\n '''\n # correlation matrix on provided DataFrame (original scale)\n try:\n corr = star_feat_df.reindex(columns=self.feature_names_, fill_value=0.0).corr()\n except Exception as e:\n print(\"[WARN] correlation calculation failed:\", e)\n corr = pd.DataFrame(np.zeros((self.n_features_in_, self.n_features_in_)), index=self.feature_names_, columns=self.feature_names_)\n\n # top correlated pairs\n pairs = []\n cols = corr.columns\n for i in range(len(cols)):\n for j in range(i+1, len(cols)):\n a = cols[i]; b = cols[j]\n val = corr.iloc[i, j]\n pairs.append((a, b, val, abs(val)))\n pairs_sorted = sorted(pairs, key=lambda x: x[3], reverse=True)\n '''\n print(f\"\\n=== Top {min(corr_top_k, len(pairs_sorted))} correlated feature pairs (abs corr desc) ===\")\n cnt = 0\n for a, b, val, aval in pairs_sorted[:corr_top_k]:\n if corr_threshold is None or aval >= corr_threshold:\n print(f\"{a:30s} <-> {b:30s} | corr = {val:.4f}\")\n cnt += 1\n if cnt == 0:\n print(f\"No pairs found with |corr| >= {corr_threshold}. Showing top {corr_top_k} anyway:\")\n for a, b, val, aval in pairs_sorted[:corr_top_k]:\n print(f\"{a:30s} <-> {b:30s} | corr = {val:.4f}\")\n '''\n # optional corr heatmap\n if save_fig:\n try:\n fig, ax = plt.subplots(figsize=(10, 10))\n im = ax.imshow(corr.values, cmap=\"RdBu_r\", vmin=-1, vmax=1)\n ax.set_xticks(np.arange(len(cols))); ax.set_yticks(np.arange(len(cols)))\n ax.set_xticklabels(cols, rotation=90, fontsize=8)\n ax.set_yticklabels(cols, fontsize=8)\n fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04)\n fig.tight_layout()\n fig.savefig(fig_name, dpi=200)\n plt.close(fig)\n print(f\"\\nCorrelation heatmap saved to {fig_name}\")\n except Exception as e:\n print(\"[WARN] failed to save correlation plot:\", e)\n\n return df_imp, corr\n\n\ndef drop_highly_correlated_by_importance(X: pd.DataFrame, df_imp: pd.DataFrame, corr_thresh: float = 0.95):\n \"\"\"\n X: DataFrame features\n df_imp: DataFrame returned by feature_importance_summary (index = feature names, contains 'score_sum')\n corr_thresh: threshold absolute correlation to consider a pair as duplicated\n returns: X_reduced, dropped_list\n \"\"\"\n Xc = X.copy()\n corr = Xc.corr().abs()\n upper = corr.where(np.triu(np.ones(corr.shape), k=1).astype(bool))\n to_drop = set()\n for col in upper.columns:\n high = upper[col][upper[col] > corr_thresh].index.tolist()\n for other in high:\n if other in to_drop or col in to_drop:\n continue\n imp_col = df_imp.loc[col, 'score_sum'] if col in df_imp.index else 0.0\n imp_other = df_imp.loc[other, 'score_sum'] if other in df_imp.index else 0.0\n # drop the less important\n if imp_col >= imp_other:\n to_drop.add(other)\n else:\n to_drop.add(col)\n dropped = [c for c in Xc.columns if c in to_drop]\n X_red = Xc.drop(columns=dropped)\n return X_red, dropped\n\n# ------------------------------\n# Submission helper\n# ------------------------------\n\nclass SubmissionGenerator:\n def __init__(self, config: Config):\n self.cfg = config\n self.sample_submission = pd.read_csv(f\"{self.cfg.DATA_PATH}/sample_submission.csv\", index_col=\"planet_id\")\n\n def create_from_full_mu(self, mu_full: np.ndarray, sigmas_full: np.ndarray, index: np.ndarray):\n n_mu = self.sample_submission.shape[1] // 2\n mu_full = np.asarray(mu_full, dtype=float)\n sigmas_full = np.asarray(sigmas_full, dtype=float)\n assert mu_full.shape[1] == n_mu, \"mu_full must have n_mu columns\"\n assert sigmas_full.shape == mu_full.shape, \"sigmas_full must match mu_full shape\"\n\n out = np.concatenate([mu_full, sigmas_full], axis=1)\n df = pd.DataFrame(out, columns=self.sample_submission.columns, index=index)\n\n # --- FIX: сбрасываем индекс в колонку planet_id ---\n df = df.reset_index()\n df.rename(columns={\"index\": \"planet_id\"}, inplace=True)\n\n # сохраняем в CSV без индекса\n df.to_csv(\"submission.csv\", index=False)\n return df\n\n# feature generator\n# --- Additional feature generator ------------------------------------------------\ndef augment_features_from_signals(X_train, X_test, pre_train, pre_test, cfg,\n savgol_windows=(11, 21, 31),\n fft_bins=8, spec_pca_components=6, lag_list=(1,2,3)):\n \"\"\"\n Add many feature variants derived from preprocessed signals.\n - pre_*: arrays shape (N, T, C) with channel0=FGS, channels1: AIRS\n - X_train/X_test: DataFrames (index aligned)\n Returns augmented (X_train_aug, X_test_aug), plus meta: ('spec_pca_cols', list)\n \"\"\"\n Xtr = X_train.copy()\n Xte = X_test.copy()\n Ntr, Ttr, Ctr = pre_train.shape\n Nte, Tte, Cte = pre_test.shape\n\n # 1) per-star white (AIRS mean) and FGS series stats with multiple savgol windows\n def compute_white(mat):\n if mat.shape[2] <= 1:\n return mat[:, :, 0] # fallback\n return np.nanmean(mat[:, :, 1:], axis=2) # (N, T)\n\n white_tr = compute_white(pre_train)\n white_te = compute_white(pre_test)\n fgs_tr = pre_train[:, :, 0]\n fgs_te = pre_test[:, :, 0]\n\n # helper to add vector to DataFrame\n def add_col(df, name, arr):\n df[name] = arr\n\n # 2) Savitzky variants: mean/std/slope for different windows\n from numpy.polynomial.polynomial import polyfit as np_polyfit\n def slope_array(mat):\n # mat: (N, T) -> slope per row (linear fit)\n N, T = mat.shape\n xs = np.arange(T)\n xs_mean = xs.mean()\n denom = ((xs - xs_mean)**2).sum()\n if denom == 0:\n return np.zeros(N)\n ys_mean = np.nanmean(mat, axis=1)\n cov = np.nansum((mat - ys_mean[:,None]) * (xs[None,:] - xs_mean), axis=1)\n return cov / denom\n\n for w in savgol_windows:\n for prefix, (mat_tr, mat_te) in [(\"white\", (white_tr, white_te)), (\"fgs\", (fgs_tr, fgs_te))]:\n # safe smoothing (skip if window too large)\n def smooth_safe(mat, w):\n out = np.zeros_like(mat)\n for i in range(mat.shape[0]):\n try:\n out[i] = savgol_filter(mat[i], w, 2, mode='nearest')\n except Exception:\n out[i] = mat[i]\n return out\n s_tr = smooth_safe(mat_tr, w)\n s_te = smooth_safe(mat_te, w)\n add_col(Xtr, f\"{prefix}_savgol_w{w}_mean\", np.nanmean(s_tr, axis=1))\n add_col(Xte, f\"{prefix}_savgol_w{w}_mean\", np.nanmean(s_te, axis=1))\n add_col(Xtr, f\"{prefix}_savgol_w{w}_std\", np.nanstd(s_tr, axis=1))\n add_col(Xte, f\"{prefix}_savgol_w{w}_std\", np.nanstd(s_te, axis=1))\n add_col(Xtr, f\"{prefix}_savgol_w{w}_slope\", slope_array(s_tr))\n add_col(Xte, f\"{prefix}_savgol_w{w}_slope\", slope_array(s_te))\n\n # 3) FFT bands on white and fgs (low/mid/high)\n def fft_bands(mat, nbins=fft_bins):\n N, T = mat.shape\n out = np.zeros((N, 3))\n for i in range(N):\n row = mat[i] - np.nanmean(mat[i])\n spec = np.abs(np.fft.rfft(row))\n if spec.sum() <= 0:\n out[i,:] = 0.0\n continue\n K = max(1, int(len(spec) / (nbins*1.0)))\n low = spec[:K].sum()\n mid = spec[K:K*4].sum() if K*4 < len(spec) else spec[K:].sum()\n high = spec[-K:].sum()\n out[i,0] = low\n out[i,1] = mid\n out[i,2] = high\n return out\n\n fft_tr = fft_bands(white_tr, nbins=fft_bins)\n fft_te = fft_bands(white_te, nbins=fft_bins)\n add_col(Xtr, \"white_fft_low\", fft_tr[:,0]); add_col(Xte, \"white_fft_low\", fft_te[:,0])\n add_col(Xtr, \"white_fft_mid\", fft_tr[:,1]); add_col(Xte, \"white_fft_mid\", fft_te[:,1])\n add_col(Xtr, \"white_fft_high\", fft_tr[:,2]); add_col(Xte, \"white_fft_high\", fft_te[:,2])\n\n fft_tr_fgs = fft_bands(fgs_tr, nbins=fft_bins)\n fft_te_fgs = fft_bands(fgs_te, nbins=fft_bins)\n add_col(Xtr, \"fgs_fft_low\", fft_tr_fgs[:,0]); add_col(Xte, \"fgs_fft_low\", fft_te_fgs[:,0])\n add_col(Xtr, \"fgs_fft_mid\", fft_tr_fgs[:,1]); add_col(Xte, \"fgs_fft_mid\", fft_te_fgs[:,1])\n add_col(Xtr, \"fgs_fft_high\", fft_tr_fgs[:,2]); add_col(Xte, \"fgs_fft_high\", fft_te_fgs[:,2])\n\n # 4) lag features: mean absolute diff at several lags for white and fgs\n for lag in lag_list:\n def lag_mean_abs(mat, lag):\n arr = np.abs(mat[:, lag:] - mat[:, :-lag])\n return np.nanmean(arr, axis=1)\n add_col(Xtr, f\"white_lag{lag}_mad\", lag_mean_abs(white_tr, lag))\n add_col(Xte, f\"white_lag{lag}_mad\", lag_mean_abs(white_te, lag))\n add_col(Xtr, f\"fgs_lag{lag}_mad\", lag_mean_abs(fgs_tr, lag))\n add_col(Xte, f\"fgs_lag{lag}_mad\", lag_mean_abs(fgs_te, lag))\n\n # 5) per-wavelength (AIRS) mean spectrum + PCA across wavelengths (train-fit -> transform both)\n if pre_train.shape[2] > 1:\n # mean spectrum per star across time: shape (N, n_wl)\n spec_tr = np.nanmean(pre_train[:, :, 1:], axis=1) # (N, WL)\n spec_te = np.nanmean(pre_test[:, :, 1:], axis=1)\n # safe fill nan\n spec_tr = np.nan_to_num(spec_tr, nan=0.0)\n spec_te = np.nan_to_num(spec_te, nan=0.0)\n # fit PCA on train\n try:\n pca_spec = PCA(n_components=min(spec_pca_components, spec_tr.shape[1], spec_tr.shape[0]), svd_solver='randomized', random_state=0)\n Z_tr = pca_spec.fit_transform(spec_tr)\n Z_te = pca_spec.transform(spec_te)\n for k in range(Z_tr.shape[1]):\n add_col(Xtr, f\"spec_pca_{k}\", Z_tr[:, k])\n add_col(Xte, f\"spec_pca_{k}\", Z_te[:, k])\n # also add spectral slope (linear fit across wavelengths)\n slopes_tr = np.polyfit(np.arange(spec_tr.shape[1]), spec_tr.T, 1)[0] if spec_tr.size else np.zeros(spec_tr.shape[0])\n # Note: np.polyfit above returns vectorized coeffs only if shapes align; do fallback:\n slopes_tr = np.array([np_polyfit(np.arange(spec_tr.shape[1]), spec_tr[i], 1)[1] if spec_tr.shape[1]>1 else 0.0 for i in range(spec_tr.shape[0])])\n slopes_te = np.array([np_polyfit(np.arange(spec_te.shape[1]), spec_te[i], 1)[1] if spec_te.shape[1]>1 else 0.0 for i in range(spec_te.shape[0])])\n add_col(Xtr, \"spec_slope\", slopes_tr)\n add_col(Xte, \"spec_slope\", slopes_te)\n except Exception:\n pass\n\n # 6) per-channel percentiles (for a few selected percentiles) but vectorized\n pct_list = (5,25,50,75,95)\n for ch in range(min(8, pre_train.shape[2])): # limit channels to first 8 to control feature explosion\n arr_tr = pre_train[:, :, ch]\n arr_te = pre_test[:, :, ch]\n for p in pct_list:\n add_col(Xtr, f\"ch{ch}_p{p}\", np.nanpercentile(arr_tr, p, axis=1))\n add_col(Xte, f\"ch{ch}_p{p}\", np.nanpercentile(arr_te, p, axis=1))\n\n # 7) counts: n_nan per star per signal\n add_col(Xtr, \"fgs_n_nan\", np.sum(~np.isfinite(fgs_tr), axis=1))\n add_col(Xte, \"fgs_n_nan\", np.sum(~np.isfinite(fgs_te), axis=1))\n add_col(Xtr, \"white_n_nan\", np.sum(~np.isfinite(white_tr), axis=1))\n add_col(Xte, \"white_n_nan\", np.sum(~np.isfinite(white_te), axis=1))\n\n # done\n return Xtr, Xte\n\n\n# ------------------------------\n# Main pipeline\n# ------------------------------\n\nif __name__ == '__main__':\n __t0 = time.perf_counter()\n cfg = Config() # This will now automatically load wavelengths.csv\n\n # --- IDs & metadata\n train_star_info = pd.read_csv(f\"{cfg.DATA_PATH}/{cfg.DATASET_TRAIN}_star_info.csv\", index_col='planet_id')\n test_star_info = pd.read_csv(f\"{cfg.DATA_PATH}/{cfg.DATASET_TEST}_star_info.csv\", index_col='planet_id')\n train_ids = train_star_info.index.astype(int).values\n test_ids = test_star_info.index.astype(int).values\n \n # --- Preprocess signals\n processor = SignalProcessor(cfg)\n print(\"[INFO] Processing TRAIN signals...\")\n pre_train = processor.process_all_data(cfg.DATASET_TRAIN, train_ids)\n print(\"[INFO] Processing TEST signals...\")\n pre_test = processor.process_all_data(cfg.DATASET_TEST, test_ids)\n \n # --- Baseline TransitModel predictions (s)\n model = TransitModel(cfg)\n print(\"[INFO] Predicting analytical s on TRAIN...\")\n s_train = model.predict_all(pre_train)\n print(\"[INFO] Predicting analytical s on TEST...\")\n s_test = model.predict_all(pre_test)\n \n # --- Load training targets and sample submission layout\n sample_sub = pd.read_csv(f\"{cfg.DATA_PATH}/sample_submission.csv\", index_col='planet_id')\n n_mu = sample_sub.shape[1] // 2\n mu_cols = list(sample_sub.columns[:n_mu])\n\n train_df = pd.read_csv(f\"{cfg.DATA_PATH}/train.csv\", index_col='planet_id')\n # best-effort alignment: use same mu columns as sample_sub\n missing = [c for c in mu_cols if c not in train_df.columns]\n if missing:\n print(f\"[WARN] Missing columns in train.csv: {missing}. ML shape may be disabled.\")\n mu_cols_present = [c for c in mu_cols if c in train_df.columns]\n if len(mu_cols_present) != n_mu:\n # fallback: use all columns that look like mu_* sorted\n mu_cols_present = sorted([c for c in train_df.columns if c.startswith('mu_')])[:n_mu]\n mu_train = train_df[mu_cols_present].reindex(train_ids).values.astype(float)\n\n\n # ------------------------------\n # Aggressive feature expansion + auto-prune\n # ------------------------------\n print(\"[INFO] Building base features...\")\n X_train = extract_features(pre_train, cfg, train_star_info)\n X_test = extract_features(pre_test, cfg, test_star_info)\n X_train = add_interaction_features(X_train)\n X_test = add_interaction_features(X_test)\n \n # --- Add Multi-Scale and Transit Quality Features ---\n print(\"[INFO] Adding multi-scale features...\")\n X_train = add_multiscale_features(X_train, pre_train, cfg)\n X_test = add_multiscale_features(X_test, pre_test, cfg)\n \n print(\"[INFO] Adding transit quality features...\")\n X_train = add_transit_quality_features(X_train, pre_train, cfg)\n X_test = add_transit_quality_features(X_test, pre_test, cfg)\n \n print(\"[INFO] Adding spectral correlation features...\")\n X_train = add_spectral_correlation_features(X_train, pre_train, cfg)\n X_test = add_spectral_correlation_features(X_test, pre_test, cfg)\n \n # --- Add Atmospheric Physics Features (using wavelengths.csv) ---\n print(\"[INFO] Adding atmospheric physics features...\")\n X_train = add_atmospheric_physics_features(X_train, pre_train, cfg)\n X_test = add_atmospheric_physics_features(X_test, pre_test, cfg)\n \n # --- Additional feature generator (user-defined above) ---\n print(\"[INFO] Augmenting features from signals (PCA, savgol variants, FFT bands, per-wl PCA...)\")\n X_train_aug, X_test_aug = augment_features_from_signals(X_train, X_test, pre_train, pre_test, cfg,\n savgol_windows=(11,21,31),\n fft_bins=8,\n spec_pca_components=6,\n lag_list=(1,2,3))\n \n X_train = X_train_aug\n X_test = X_test_aug\n \n # --- Enhanced Transit Quality Detection ---\n print(\"[INFO] Computing transit quality scores...\")\n transit_quality_train = detect_malformed_transits(pre_train, cfg)\n transit_quality_test = detect_malformed_transits(pre_test, cfg)\n \n X_train['transit_quality'] = transit_quality_train\n X_test['transit_quality'] = transit_quality_test\n \n # print(f\"[INFO] Transit quality stats - Train: mean={np.mean(transit_quality_train):.3f}, \"\n # f\"std={np.std(transit_quality_train):.3f}, min={np.min(transit_quality_train):.3f}\")\n # print(f\"[INFO] Low quality transits (< 0.5): {np.sum(transit_quality_train < 0.5)} / {len(transit_quality_train)}\")\n \n # quick sanitization: drop constant or near-constant features\n var = X_train.var(axis=0)\n const_feats = var[var <= 1e-12].index.tolist()\n if const_feats:\n print(f\"[INFO] Dropping {len(const_feats)} constant/near-constant features: {const_feats[:20]}\")\n X_train = X_train.drop(columns=const_feats)\n X_test = X_test.drop(columns=const_feats)\n \n # ------------------------------\n # Initial SafeML calibrator pass (feature importances)\n # ------------------------------\n print(\"[INFO] Running SafeMLCalibrator to compute feature importances (first pass)...\")\n cal_tmp = SafeMLCalibrator(cfg, n_components=min(4, max(2, n_mu//50)))\n cal_tmp.fit(pre_train, s_train, mu_train, X_train)\n \n df_imp, corr = cal_tmp.feature_importance_summary(X_train, top_n=40, corr_top_k=80, corr_threshold=0.5, save_fig=False)\n \n # ------------------------------\n # Drop low-importance + correlated features\n # ------------------------------\n imp_threshold = 1e-3\n low_imp = df_imp[df_imp['score_sum'] <= imp_threshold].index.tolist()\n print(f\"[INFO] Dropping {len(low_imp)} low-importance features (score_sum <= {imp_threshold}).\")\n \n X_train_red, dropped_corr = drop_highly_correlated_by_importance(X_train, df_imp, corr_thresh=0.95)\n dropped_total = set(low_imp) | set(dropped_corr)\n X_train_red = X_train_red.drop(columns=[c for c in dropped_total if c in X_train_red.columns], errors='ignore')\n X_test_red = X_test.drop(columns=[c for c in dropped_total if c in X_test.columns], errors='ignore')\n \n print(\"[INFO] Total dropped features count:\", len(dropped_total))\n try:\n with open(\"dropped_features.txt\", \"w\") as fh:\n for c in sorted(dropped_total):\n fh.write(c + \"\\n\")\n print(\"[INFO] Dropped feature names saved to dropped_features.txt\")\n except Exception as e:\n print(\"[WARN] Could not save dropped features file:\", e)\n \n X_train = X_train_red.copy()\n X_test = X_test_red.copy()\n \n X_train = X_train.replace([np.inf, -np.inf], np.nan).fillna(0.0)\n X_test = X_test.replace([np.inf, -np.inf], np.nan).fillna(0.0)\n \n print(f\"[INFO] Feature matrix sizes: X_train {X_train.shape}, X_test {X_test.shape}\")\n # -------------------------------------------------------------------------------\n\n\n # ------------------------------\n # Robust LightGBM OOF feature (mu_level)\n # ------------------------------\n from sklearn.model_selection import KFold\n from sklearn.metrics import mean_squared_error\n \n mu_level = mu_train.mean(axis=1)\n \n Xtr_mat = X_train.values.astype(np.float32, copy=False)\n Xte_mat = X_test.values.astype(np.float32, copy=False)\n \n oof_preds = np.zeros(Xtr_mat.shape[0], dtype=float)\n test_preds_folds = []\n \n n_splits = 5\n kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)\n \n use_lgb = True\n try:\n import lightgbm as lgb\n from lightgbm import LGBMRegressor\n use_lgb = True\n print(\"[INFO] lightgbm detected (will train without early_stopping to be robust).\")\n except Exception as e:\n print(\"[WARN] lightgbm import failed, falling back to sklearn. Error:\", e)\n from sklearn.ensemble import HistGradientBoostingRegressor\n \n if use_lgb:\n lgb_params = {\n \"n_estimators\": 500,\n \"learning_rate\": 0.01,\n \"num_leaves\": 31,\n \"max_depth\": 7,\n \"min_data_in_leaf\": 20,\n \"verbosity\": -1,\n \"random_state\": 42,\n \"n_jobs\": -1\n }\n else:\n hgb_params = {\n \"max_iter\": 300,\n \"learning_rate\": 0.05,\n \"max_depth\": 7,\n \"random_state\": 42\n }\n \n print(\"[INFO] Running OOF LightGBM/HGB for mu_level (no early stopping)...\")\n for fold, (tr_idx, va_idx) in enumerate(kf.split(Xtr_mat)):\n X_tr, X_va = Xtr_mat[tr_idx], Xtr_mat[va_idx]\n y_tr, y_va = mu_level[tr_idx], mu_level[va_idx]\n \n if use_lgb:\n try:\n mdl = LGBMRegressor(**lgb_params)\n mdl.fit(X_tr, y_tr)\n pred_va = mdl.predict(X_va)\n pred_test = mdl.predict(Xte_mat)\n except Exception as e:\n print(f\"[WARN] LGBMRegressor.fit failed on fold {fold}, error: {e}. Falling back to sklearn HGB.\")\n from sklearn.ensemble import HistGradientBoostingRegressor\n mdl2 = HistGradientBoostingRegressor(**hgb_params)\n mdl2.fit(X_tr, y_tr)\n pred_va = mdl2.predict(X_va)\n pred_test = mdl2.predict(Xte_mat)\n else:\n mdl = HistGradientBoostingRegressor(**hgb_params)\n mdl.fit(X_tr, y_tr)\n pred_va = mdl.predict(X_va)\n pred_test = mdl.predict(Xte_mat)\n \n oof_preds[va_idx] = pred_va\n test_preds_folds.append(pred_test)\n mse_fold = mean_squared_error(y_va, pred_va)\n print(f\"[OOF] fold {fold+1}/{n_splits} MSE: {mse_fold:.5e}\")\n \n test_pred = np.mean(np.vstack(test_preds_folds), axis=0) if len(test_preds_folds) else np.zeros(Xte_mat.shape[0], dtype=float)\n oof_mse = mean_squared_error(mu_level, oof_preds)\n print(f\"[OOF] Overall OOF MSE (mu_level): {oof_mse:.6e}\")\n \n X_train[\"lgb_oof_mu\"] = oof_preds.astype(np.float32)\n X_test[\"lgb_mu\"] = test_pred.astype(np.float32)\n \n print(f\"[INFO] Added 'lgb_oof_mu' and 'lgb_mu'. Shapes: {X_train.shape}, {X_test.shape}\")\n # ------------------------------\n\n # ------------------------------\n # Safe ML calibration (with c*-predictor for sigma scaling)\n # ------------------------------\n calibrator = None\n mu_test_final = None\n sigmas_test_final = None\n reg_c = None # will hold Ridge regressor for c* if trained\n try:\n print(\"[INFO] Fitting Safe ML Calibrator...\")\n calibrator = SafeMLCalibrator(cfg, n_components=min(4, max(2, n_mu//50)))\n # first fit on current features\n calibrator.fit(pre_train, s_train, mu_train, X_train)\n df_imp, corr = calibrator.feature_importance_summary(X_train, top_n=30, corr_top_k=50, corr_threshold=0.6, save_fig=True)\n\n # drop correlated by importance\n X_train_red, dropped = drop_highly_correlated_by_importance(X_train, df_imp, corr_thresh=0.95)\n X_test_red = X_test[X_train_red.columns]\n\n print(\"Dropped features (auto):\", dropped)\n\n # retrain calibrator on reduced feature set\n calibrator.fit(pre_train, s_train, mu_train, X_train_red)\n df_imp_new, corr_new = calibrator.feature_importance_summary(X_train_red, top_n=30, save_fig=True)\n\n # align X_train/X_test to reduced columns for subsequent modelling (including c* regressor)\n X_train = X_train_red.copy()\n X_test = X_test_red.copy()\n\n # ------------------------------\n # --- TRAIN c* regressor on train OOF (analytic optimal multiplier)\n # ------------------------------\n # helper: align per-wl sigma-form to mu air channels count\n def align_shape_to_mu(sigma_shape, n_mu_minus1):\n # sigma_shape: (N, m); we resample along axis=1 to length n_mu_minus1\n if sigma_shape.ndim != 2:\n return sigma_shape\n N, m = sigma_shape.shape\n if m == n_mu_minus1:\n return sigma_shape\n x_orig = np.linspace(0.0, 1.0, m)\n x_new = np.linspace(0.0, 1.0, n_mu_minus1)\n out = np.zeros((N, n_mu_minus1), dtype=float)\n for i in range(N):\n out[i] = np.interp(x_new, x_orig, sigma_shape[i])\n return out\n\n # obtain ML mu prediction on TRAIN (base predicted shape)\n mu_base_train = calibrator.infer_mu(s_train, X_train, n_mu) # (N_train, n_mu)\n\n # build per-wavelength empirical sigma shape from pre_train (std over time per wl for AIRS channels)\n if pre_train.ndim == 3 and pre_train.shape[2] > 1 and n_mu > 1:\n # pre_train: (N, n_bins, channels_total) where channels_total = 1 + n_air_ch\n sigma_shape_train = np.nanstd(pre_train[:, :, 1:], axis=1) # (N, n_air_ch)\n # align to mu air channels (n_mu-1)\n sigma_shape_train_al = align_shape_to_mu(sigma_shape_train, n_mu - 1)\n # compute c_star using user-provided compute_c_star_per_planet if available, else fallback\n try:\n c_star_train = compute_c_star_per_planet(mu_train[:, 1:], mu_base_train[:, 1:], sigma_shape_train_al)\n except Exception:\n # fallback analytic computation if function missing\n eps = 1e-12\n R2_div_S2 = ((mu_train[:,1:] - mu_base_train[:,1:])**2) / (np.clip(sigma_shape_train_al, eps, None)**2)\n c_star_train = np.sqrt(np.nanmean(R2_div_S2, axis=1))\n else:\n # if no AIRS channels or shapes mismatch -> fallback to scalar c* based on aggregated sigma\n agg_sigma_train = np.nanmean(np.nanstd(pre_train, axis=1), axis=1) # (N,)\n # avoid zero\n agg_sigma_train = np.clip(agg_sigma_train, 1e-9, None)\n mu_base_train = calibrator.infer_mu(s_train, X_train, n_mu)\n if n_mu > 1:\n resid_air = mu_train[:,1:] - mu_base_train[:,1:]\n # collapse resid to scalar per planet (mean squared across channels)\n sq = np.nanmean(resid_air**2, axis=1)\n else:\n sq = np.nanmean((mu_train - mu_base_train)**2, axis=1)\n c_star_train = np.sqrt(np.clip(sq, 0, None) / (agg_sigma_train**2 + 1e-12))\n\n # train small Ridge on log(c*)\n from sklearn.linear_model import Ridge\n # sanitize targets\n c_star_train = np.nan_to_num(c_star_train, nan=1.0, posinf=1.0, neginf=1.0)\n c_star_train = np.clip(c_star_train, 1e-6, 1e3)\n logc = np.log(c_star_train)\n reg_c = Ridge(alpha=1.0).fit(X_train.values, logc)\n print(\"[INFO] Trained c* Ridge regressor on train (log-space).\")\n\n # ------------------------------\n # infer μ for test\n mu_test_ml = calibrator.infer_mu(s_test, X_test, n_mu)\n\n # --- baseline sigmas (scalar per planet for AIRS and FGS)\n sigma_fgs_vec = estimate_sigma_fgs(pre_test, cfg) # shape (N_test,)\n sigma_air_vec = estimate_sigma_air(pre_test, cfg) # shape (N_test,)\n\n sigmas_base = np.full_like(mu_test_ml, cfg.SIGMA, dtype=float)\n sigmas_base[:,0] = np.clip(sigma_fgs_vec, 1e-6, 0.1)\n sigmas_base[:,1:] = np.clip(sigma_air_vec.reshape(-1,1), 1e-6, 0.1)\n\n # --- apply c* predicted scaling to AIRS sigmas (per-planet)\n if reg_c is not None and n_mu > 1:\n # compute per-planet empirical sigma_shape on test, align to mu channels\n if pre_test.ndim == 3 and pre_test.shape[2] > 1:\n sigma_shape_test = np.nanstd(pre_test[:, :, 1:], axis=1) # (N_test, n_air_ch)\n sigma_shape_test_al = align_shape_to_mu(sigma_shape_test, n_mu - 1)\n else:\n sigma_shape_test_al = None\n\n # predict c on X_test\n try:\n logc_test = reg_c.predict(X_test.values)\n c_pred_test = np.exp(logc_test)\n # clip to reasonable range to avoid pathological scaling\n c_pred_test = np.clip(c_pred_test, 0.3, 3.0)\n # apply multiplicative scaling to AIRS columns\n sigmas_base[:, 1:] = sigmas_base[:, 1:] * c_pred_test[:, None]\n print(\"[INFO] Applied predicted c* scaling to AIRS sigmas (clipped to [0.3,3.0]).\")\n except Exception as e:\n print(\"[WARN] c* regressor prediction failed on test: \", e)\n\n # refine sigmas via ML (geom blend)\n sigmas_ml = calibrator.infer_sigma(X_test, cfg.SIGMA, n_mu)\n\n # blend base and ml in log space (as previously implemented)\n log_final = 0.05*np.log(np.clip(sigmas_base,1e-12,0.1)) + 0.95*np.log(np.clip(sigmas_ml,1e-12,0.1))\n sigmas_test_final = np.exp(log_final)\n mu_test_final = mu_test_ml\n # Apply spectral smoothness constraint for physically realistic transmission spectra\n print(\"[INFO] Applying spectral smoothness constraint...\")\n mu_test_final = apply_spectral_smoothness_constraint(mu_test_final, sigmas_test_final, lambda_smooth=0.1)\n except Exception as e:\n print(f\"[WARN] ML calibration failed with: {e}\\nFalling back to baseline flat μ and baseline σ.\")\n # baseline μ (flat)\n mu_test_final = np.tile(s_test.reshape(-1,1), (1, n_mu))\n # baseline σ\n sigma_fgs_vec = estimate_sigma_fgs(pre_test, cfg)\n sigma_air_vec = estimate_sigma_air(pre_test, cfg)\n sigmas_test_final = np.full_like(mu_test_final, cfg.SIGMA, dtype=float)\n sigmas_test_final[:,0] = np.clip(sigma_fgs_vec, 1e-6, 0.1)\n sigmas_test_final[:,1:] = np.clip(sigma_air_vec.reshape(-1,1), 1e-6, 0.1)\n # Apply spectral smoothness constraint even in fallback case\n print(\"[INFO] Applying spectral smoothness constraint (fallback mode)...\")\n mu_test_final = apply_spectral_smoothness_constraint(mu_test_final, sigmas_test_final, lambda_smooth=0.05)\n # --- Create submission\n print(\"[INFO] Creating submission.csv ...\")\n subgen = SubmissionGenerator(cfg)\n submission = subgen.create_from_full_mu(mu_test_final, sigmas_test_final, index=test_ids)\n\n __t1 = time.perf_counter()\n print(f\"[TIMING] total runtime: {(__t1-__t0):.2f} s ({(__t1-__t0)/60:.2f} min)\")\n print(submission.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-18T00:14:14.175609Z","iopub.execute_input":"2025-09-18T00:14:14.175965Z","iopub.status.idle":"2025-09-18T00:15:03.850118Z","shell.execute_reply.started":"2025-09-18T00:14:14.175933Z","shell.execute_reply":"2025-09-18T00:15:03.846576Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport os\n\n# Пути – проверь, что sample_submission.csv есть в датасете\nsample_path = \"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\"\nsubmission_path = \"submission.csv\" # Файл, который ты создаёшь\n\n# Проверим наличие файлов\nfor p in [sample_path, submission_path]:\n if not os.path.exists(p):\n print(f\"Файл не найден: {p}\")\n else:\n print(f\"Файл найден: {p}\")\n\ntry:\n sample = pd.read_csv(sample_path)\n sub = pd.read_csv(submission_path)\nexcept Exception as e:\n print(\"Ошибка чтения CSV:\", e)\n raise\n\nprint(\"\\nКолонки sample_submission:\", sample.columns.tolist())\nprint(\"Колонки submission:\", sub.columns.tolist())\n\n# Проверка колонок\nif list(sample.columns) != list(sub.columns):\n print(\"\\n⚠️ Колонки не совпадают по именам или порядку!\")\n missing = set(sample.columns) - set(sub.columns)\n extra = set(sub.columns) - set(sample.columns)\n if missing:\n print(\" Не хватает колонок:\", missing)\n if extra:\n print(\" Лишние колонки:\", extra)\nelse:\n print(\"\\n✔️ Колонки полностью совпадают с образцом.\")\n\n# Проверка числа строк\nif len(sample) != len(sub):\n print(f\"⚠️ Количество строк не совпадает: sample={len(sample)}, submission={len(sub)}\")\nelse:\n print(f\"✔️ Количество строк совпадает: {len(sub)}\")\n\n# Проверка порядка id (берём первую колонку как id)\nid_col = sample.columns[0]\nif all(sample[id_col] == sub[id_col]):\n print(\"✔️ Порядок id совпадает с образцом.\")\nelse:\n print(\"⚠️ Порядок id отличается от sample_submission!\")\n print(\" Например, первые 5 id в sample:\", sample[id_col].head().tolist())\n print(\" Первые 5 id в submission:\", sub[id_col].head().tolist())\n","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}