{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"isSourceIdPinned":false,"sourceType":"competition"}],"dockerImageVersionId":31192,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2025-11-14T19:45:04.931803Z","iopub.execute_input":"2025-11-14T19:45:04.932119Z","iopub.status.idle":"2025-11-14T19:45:05.400555Z","shell.execute_reply.started":"2025-11-14T19:45:04.932083Z","shell.execute_reply":"2025-11-14T19:45:05.399569Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport lightgbm as lgb\nfrom sklearn.model_selection import KFold\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import mean_squared_error, r2_score\nfrom tqdm import tqdm\nimport time\nimport os\nimport matplotlib.pyplot as plt\nimport glob\nfrom scipy import signal, stats\nfrom scipy.stats import skew, kurtosis\nfrom scipy.optimize import minimize_scalar\nfrom scipy.stats import norm\nimport warnings\nimport sys\nimport joblib\nimport pickle\n\nwarnings.filterwarnings('ignore')\n\nclass Config:\n    \"\"\"\n    Central configuration class. All paths, model parameters,\n    and feature engineering settings are defined here.\n    \"\"\"\n    BASE_PATH = '/kaggle/input/ariel-data-challenge-2025/'\n    TRAIN_CSV = f\"{BASE_PATH}/train.csv\"\n    TRAIN_STAR_INFO_CSV = f\"{BASE_PATH}/train_star_info.csv\"\n    TEST_STAR_INFO_CSV = f\"{BASE_PATH}/test_star_info.csv\"\n    ADC_INFO_CSV = f\"{BASE_PATH}/adc_info.csv\"\n    TRAIN_DATA_PATH = f\"{BASE_PATH}/train/\"\n    TEST_DATA_PATH = f\"{BASE_PATH}/test/\"\n\n    SUBMISSION_FILE = 'submission.csv'\n    MODEL_DIR = 'trained_models'\n    IMPORTANCE_FILE = 'feature_importance.csv'\n\n    QUICK_RUN = False\n    QUICK_RUN_SIZE = 200\n\n    ROLLING_WINDOWS = [50, 100, 500, 1000, 2000, 5000]\n    DETREND_WINDOW = 5000\n    AUTOCORR_LAGS = [10, 50, 100, 500, 1000]\n    GRADIENT_WINDOWS = [10, 50, 100]\n    N_SEGMENTS = 4\n\n    N_SPLITS = 5\n    LGBM_PARAMS = {\n        \"n_estimators\": 500,\n        \"learning_rate\": 0.05,\n        \"num_leaves\": 127,\n        \"max_depth\": 8,\n        \"min_child_samples\": 20,\n        \"subsample\": 0.8,\n        \"colsample_bytree\": 0.8,\n        \"reg_alpha\": 0.1,\n        \"reg_lambda\": 0.1,\n        \"random_state\": 42,\n        \"verbose\": -1\n    }\n    \n    FSG_SIGMA_TRUE = 1e-6\n    AIRS_SIGMA_TRUE = 1e-5\n    FGS_WEIGHT = 0.4\n    SIGMA_OPTIMIZE_BOUNDS = (0.25, 4.0)\n\ndef load_initial_data(cfg):\n    \"\"\"Loads the initial CSV files.\"\"\"\n    train_df = pd.read_csv(cfg.TRAIN_CSV)\n    print(\"Train\")\n    display(train_df.head(2))\n\n    train_star_info_df = pd.read_csv(cfg.TRAIN_STAR_INFO_CSV)\n    print(\"Train Star Info\")\n    display(train_star_info_df.head(2))\n\n    test_star_info_df = pd.read_csv(cfg.TEST_STAR_INFO_CSV)\n    print(\"Test\")\n    display(test_star_info_df)\n\n    if len(test_star_info_df) > 1:\n        cfg.QUICK_RUN = False\n\n    if cfg.QUICK_RUN:\n        train_star_info_df = train_star_info_df.head(cfg.QUICK_RUN_SIZE)\n        print(f\"Only loading {cfg.QUICK_RUN_SIZE} records - this run is not a full training!\")\n    \n    return train_df, train_star_info_df, test_star_info_df\n\ndef load_demo_data(cfg, train_star_info_df):\n    \"\"\"Loads a single planet's data for demonstration.\"\"\"\n    planet_id = int(train_star_info_df.iloc[0]['planet_id'])\n    file_path = f\"{cfg.TRAIN_DATA_PATH}{planet_id}/FGS1_signal_0.parquet\"\n    df = pd.read_parquet(file_path)\n    signal_data = df.values.reshape(135000, 32, 32)\n    print(f\"Demo data loaded for Planet ID: {planet_id}\")\n    print(f\"Signal data shape: {signal_data.shape}\")\n    print(\"=\"*60)\n    return signal_data, planet_id\n\ndef apply_adc_correction(signal_array, instrument, adc_info_path):\n    \"\"\"Applies ADC correction to raw signal data using known gain and offset.\"\"\"\n    adc_df = pd.read_csv(adc_info_path)\n    gain_col = f\"{instrument}_adc_gain\"\n    offset_col = f\"{instrument}_adc_offset\"\n\n    if gain_col not in adc_df.columns or offset_col not in adc_df.columns:\n        raise ValueError(f\"Missing columns: {gain_col} or {offset_col} in {adc_info_path}\")\n    \n    gain = adc_df.at[0, gain_col]\n    offset = adc_df.at[0, offset_col]\n\n    calibrated_signal = signal_array.astype(np.float32) * gain + offset\n    return calibrated_signal\n\ndef visualize_transit_data(signal_data, total_flux, planet_id):    \n    base_frames = [0, len(signal_data)//4, len(signal_data)//2, 3*len(signal_data)//4, -1]\n\n    brightest_idx = np.argmax(total_flux)\n    darkest_idx = np.argmin(total_flux)\n\n    extra_frames = []\n    for idx in [brightest_idx, darkest_idx]:\n        if idx not in base_frames and (idx != -1 and idx != len(signal_data)-1):\n            extra_frames.append(idx)\n    frames = base_frames + extra_frames\n\n    plt.style.use('default')\n    fig = plt.figure(figsize=(18, 10))\n    \n    for i, idx in enumerate(frames):\n        ax = plt.subplot(2, len(frames), i + 1)\n        frame = signal_data[idx]\n        \n        vmin, vmax = np.percentile(frame, [2, 98])\n        im = ax.imshow(frame, cmap='hot', aspect='equal', vmin=vmin, vmax=vmax)\n        \n        time_min = idx * 0.1 / 60\n        if idx == brightest_idx:\n            title = f'Brightest\\nT={time_min:.1f} min'\n        elif idx == darkest_idx:\n            title = f'Darkest\\nT={time_min:.1f} min'\n        else:\n            title = f'T={time_min:.1f} min'\n        ax.set_title(title, fontsize=11)\n        ax.set_xticks([])\n        ax.set_yticks([])\n    \n    ax_light = plt.subplot(2, 1, 2)\n    time_hours = np.arange(len(total_flux)) * 0.1 / 3600\n    sample = slice(None, None, max(1, len(total_flux)//2000))\n\n    ax_light.plot(time_hours[sample], total_flux[sample], \n                  color='lightsteelblue', alpha=0.4, linewidth=0.5, label='Raw flux')\n\n    window = 500\n    moving_avg = pd.Series(total_flux).rolling(window, center=True).mean()\n    ax_light.plot(time_hours[sample], moving_avg.iloc[sample], \n                  color='darkblue', linewidth=3, label=f'{window}-frame average')\n\n    colors = ['red', 'orange', 'green', 'purple', 'brown', 'lime', 'black']\n    for i, idx in enumerate(frames):\n        time_point = idx * 0.1 / 3600\n        ax_light.axvline(time_point, color=colors[i % len(colors)], alpha=0.8, linewidth=1.5, linestyle='--')\n\n    ax_light.set_xlabel('Time (hours)', fontsize=12)\n    ax_light.set_ylabel('Total Flux (counts)', fontsize=12)\n    ax_light.set_title(f'Transit Light Curve - Planet {planet_id}', fontsize=14, pad=15)\n    ax_light.grid(True, alpha=0.3, linestyle='-', linewidth=0.5)\n    ax_light.legend(loc='upper right', framealpha=0.9)\n\n    smooth_flux = moving_avg.dropna()\n    flux_min, flux_max = smooth_flux.min(), smooth_flux.max()\n    flux_range = flux_max - flux_min\n    margin = flux_range * 0.1\n    ax_light.set_ylim(flux_min - margin, flux_max + margin)\n\n    plt.tight_layout(pad=2.0)\n    plt.show()\n\n    duration = time_hours[-1]\n    transit_depth = np.mean(total_flux) - np.min(total_flux)\n\n    print(f\"🌟 Planet {planet_id} Transit Observation\")\n    print(f\"    Duration: {duration:.2f} hours ({len(signal_data):,} frames)\")\n    print(f\"    Brightness: {np.min(total_flux):,.0f} → {np.max(total_flux):,.0f} counts\")\n    print(f\"    Brightest Frame: {brightest_idx} | Darkest Frame: {darkest_idx}\")\n    print(f\"    Transit depth: {transit_depth:,.0f} counts ({transit_depth/np.mean(total_flux)*100:.3f}%)\")\n\ndef extract_global_flux_features(total_flux):\n    features = {}\n    features['global_flux_mean'] = np.mean(total_flux)\n    features['global_flux_std'] = np.std(total_flux)\n    features['global_flux_min'] = np.min(total_flux)\n    features['global_flux_max'] = np.max(total_flux)\n    features['global_flux_range'] = features['global_flux_max'] - features['global_flux_min']\n    features['global_flux_skew'] = skew(total_flux)\n    features['global_flux_kurtosis'] = kurtosis(total_flux)\n    features['global_flux_cv'] = features['global_flux_std'] / features['global_flux_mean']\n    for p in [1, 5, 10, 25, 50, 75, 90, 95, 99]:\n        features[f'global_flux_p{p}'] = np.percentile(total_flux, p)\n    features['global_flux_depth'] = features['global_flux_mean'] - features['global_flux_min']\n    features['global_flux_depth_ratio'] = features['global_flux_depth'] / features['global_flux_mean']\n    return features\n\ndef extract_rolling_statistics_features(total_flux, window_sizes):\n    features = {}\n    for window in window_sizes:\n        if window < len(total_flux):\n            rolling_mean = pd.Series(total_flux).rolling(window=window, center=True).mean().dropna()\n            rolling_std = pd.Series(total_flux).rolling(window=window, center=True).std().dropna()\n            rolling_min = pd.Series(total_flux).rolling(window=window, center=True).min().dropna()\n            rolling_max = pd.Series(total_flux).rolling(window=window, center=True).max().dropna()\n            if len(rolling_mean) > 0:\n                features[f'rolling{window}_mean_min'] = rolling_mean.min()\n                features[f'rolling{window}_mean_max'] = rolling_mean.max()\n                features[f'rolling{window}_mean_std'] = rolling_mean.std()\n                features[f'rolling{window}_mean_range'] = rolling_mean.max() - rolling_mean.min()\n                features[f'rolling{window}_std_mean'] = rolling_std.mean()\n                features[f'rolling{window}_std_max'] = rolling_std.max()\n                features[f'rolling{window}_deepest_dip'] = rolling_min.min()\n                features[f'rolling{window}_highest_peak'] = rolling_max.max()\n                features[f'rolling{window}_volatility'] = rolling_std.std()\n    return features\n\ndef extract_transit_detection_features(total_flux, detrend_window):\n    features = {}\n    baseline = pd.Series(total_flux).rolling(window=detrend_window, center=True).median().fillna(method='bfill').fillna(method='ffill')\n    detrended = total_flux - baseline\n    features['detrended_min'] = np.min(detrended)\n    features['detrended_std'] = np.std(detrended)\n    features['detrended_skew'] = skew(detrended)\n    features['detrended_neg_excursions'] = np.sum(detrended < -2 * np.std(detrended))\n    features['detrended_deep_excursions'] = np.sum(detrended < -3 * np.std(detrended))\n    threshold = np.mean(total_flux) - 1.0 * np.std(total_flux)\n    below_threshold = total_flux < threshold\n    if np.any(below_threshold):\n        diff = np.diff(np.concatenate(([False], below_threshold, [False])).astype(int))\n        starts = np.where(diff == 1)[0]\n        ends = np.where(diff == -1)[0]\n        durations = ends - starts\n        features['longest_dip_duration'] = np.max(durations) if len(durations) > 0 else 0\n        features['num_dip_periods'] = len(durations)\n        features['total_dip_time'] = np.sum(durations)\n        features['avg_dip_duration'] = np.mean(durations) if len(durations) > 0 else 0\n        deepest_idx = np.argmin(total_flux)\n        features['deepest_time_fraction'] = deepest_idx / len(total_flux)\n        features['deepest_in_first_half'] = float(deepest_idx < len(total_flux) / 2)\n        features['deepest_in_middle_third'] = float(len(total_flux) / 3 < deepest_idx < 2 * len(total_flux) / 3)\n        transit_flux = total_flux[below_threshold]\n        features['transit_depth_mean'] = np.mean(transit_flux)\n        features['transit_depth_std'] = np.std(transit_flux)\n        features['transit_assymetry'] = skew(transit_flux)\n        features['transit_flatness'] = kurtosis(transit_flux)\n    else:\n        features.update({\n            'longest_dip_duration': 0, 'num_dip_periods': 0, 'total_dip_time': 0, 'avg_dip_duration': 0,\n            'deepest_time_fraction': 0.5, 'deepest_in_first_half': 0, 'deepest_in_middle_third': 0,\n            'transit_depth_mean': np.mean(total_flux), 'transit_depth_std': 0, \n            'transit_assymetry': 0, 'transit_flatness': 0\n        })\n    first_quarter = total_flux[:len(total_flux)//4]\n    last_quarter = total_flux[-len(total_flux)//4:]\n    middle_half = total_flux[len(total_flux)//4:-len(total_flux)//4]\n    features['first_quarter_mean'] = np.mean(first_quarter)\n    features['last_quarter_mean'] = np.mean(last_quarter)\n    features['middle_half_mean'] = np.mean(middle_half)\n    features['middle_vs_edges'] = features['middle_half_mean'] - (features['first_quarter_mean'] + features['last_quarter_mean']) / 2\n    return features\n\ndef extract_frequency_features(total_flux, lag_list):\n    features = {}\n    fft_flux = np.fft.fft(total_flux - np.mean(total_flux))\n    fft_power = np.abs(fft_flux)\n    fft_freqs = np.fft.fftfreq(len(total_flux))\n    features['fft_peak_power'] = np.max(fft_power[1:len(fft_power)//2])\n    features['fft_total_power'] = np.sum(fft_power[1:len(fft_power)//2])\n    features['fft_mean_power'] = np.mean(fft_power[1:len(fft_power)//2])\n    features['fft_std_power'] = np.std(fft_power[1:len(fft_power)//2])\n    low_freq_mask = np.abs(fft_freqs) < 0.01\n    features['fft_low_freq_power'] = np.sum(fft_power[low_freq_mask])\n    features['fft_low_freq_ratio'] = features['fft_low_freq_power'] / features['fft_total_power']\n    power_spectrum = fft_power[1:len(fft_power)//2]\n    freqs = np.abs(fft_freqs[1:len(fft_freqs)//2])\n    if np.sum(power_spectrum) > 0:\n        features['spectral_centroid'] = np.sum(freqs * power_spectrum) / np.sum(power_spectrum)\n        features['spectral_bandwidth'] = np.sqrt(np.sum(((freqs - features['spectral_centroid'])**2) * power_spectrum) / np.sum(power_spectrum))\n    else:\n        features['spectral_centroid'] = 0\n        features['spectral_bandwidth'] = 0\n    autocorr = np.correlate(total_flux - np.mean(total_flux), total_flux - np.mean(total_flux), mode='full')\n    autocorr = autocorr[autocorr.size // 2:]\n    autocorr = autocorr / autocorr[0]\n    for lag in lag_list:\n        if lag < len(autocorr):\n            features[f'autocorr_lag{lag}'] = autocorr[lag]\n    peaks, _ = signal.find_peaks(autocorr[1:1000], height=0.1)\n    features['autocorr_num_peaks'] = len(peaks)\n    features['autocorr_first_peak'] = peaks[0] if len(peaks) > 0 else 0\n    features['autocorr_strongest_peak'] = np.max(autocorr[peaks]) if len(peaks) > 0 else 0\n    return features\n\ndef extract_spatial_features(signal_data):\n    features = {}\n    key_frames = [0, len(signal_data) // 4, len(signal_data) // 2, 3 * len(signal_data) // 4, -1]\n    centroids_x, centroids_y, concentrations = [], [], []\n    for i, frame_idx in enumerate(key_frames):\n        frame = signal_data[frame_idx]\n        y_indices, x_indices = np.indices(frame.shape)\n        total_spatial_flux = np.sum(frame)\n        if total_spatial_flux > 0:\n            centroid_x = np.sum(x_indices * frame) / total_spatial_flux\n            centroid_y = np.sum(y_indices * frame) / total_spatial_flux\n        else:\n            centroid_x = centroid_y = 16\n        centroids_x.append(centroid_x)\n        centroids_y.append(centroid_y)\n        center_region = frame[12:20, 12:20]\n        concentration = np.sum(center_region) / total_spatial_flux if total_spatial_flux > 0 else 0\n        concentrations.append(concentration)\n        features[f'frame{i}_spatial_mean'] = np.mean(frame)\n        features[f'frame{i}_spatial_std'] = np.std(frame)\n        features[f'frame{i}_centroid_x'] = centroid_x\n        features[f'frame{i}_centroid_y'] = centroid_y\n        features[f'frame{i}_concentration'] = concentration\n    features['centroid_x_range'] = np.max(centroids_x) - np.min(centroids_x)\n    features['centroid_y_range'] = np.max(centroids_y) - np.min(centroids_y)\n    features['centroid_total_movement'] = np.sum(np.sqrt(np.diff(centroids_x)**2 + np.diff(centroids_y)**2))\n    features['concentration_range'] = np.max(concentrations) - np.min(concentrations)\n    features['concentration_std'] = np.std(concentrations)\n    return features\n\ndef extract_gradient_features(total_flux, window_list, n_segments):\n    features = {}\n    total_flux = np.array(total_flux, dtype=np.float64)\n    flux_diff1 = np.diff(total_flux)\n    flux_diff2 = np.diff(flux_diff1)\n    features['flux_diff1_mean'] = np.mean(flux_diff1)\n    features['flux_diff1_std'] = np.std(flux_diff1)\n    features['flux_diff1_min'] = np.min(flux_diff1)\n    features['flux_diff1_max'] = np.max(flux_diff1)\n    features['flux_diff1_range'] = features['flux_diff1_max'] - features['flux_diff1_min']\n    features['flux_diff1_skew'] = skew(flux_diff1)\n    features['flux_diff1_kurtosis'] = kurtosis(flux_diff1)\n    features['flux_diff2_mean'] = np.mean(flux_diff2)\n    features['flux_diff2_std'] = np.std(flux_diff2)\n    features['flux_diff2_extremes'] = np.sum(np.abs(flux_diff2) > 3 * np.std(flux_diff2))\n    for window in window_list:\n        if window < len(flux_diff1):\n            rolling_diff_std = pd.Series(flux_diff1).rolling(window=window).std()\n            features[f'diff_volatility_w{window}_max'] = rolling_diff_std.max()\n            features[f'diff_volatility_w{window}_mean'] = rolling_diff_std.mean()\n    segments = n_segments\n    segment_size = len(total_flux) // segments\n    segment_means = []\n    for i in range(segments):\n        start_idx = i * segment_size\n        end_idx = (i + 1) * segment_size if i < segments - 1 else len(total_flux)\n        segment_flux = total_flux[start_idx:end_idx]\n        if len(segment_flux) > 1:\n            trend = np.polyfit(np.arange(len(segment_flux)), segment_flux, 1)[0]\n            features[f'segment{i}_trend'] = trend\n            features[f'segment{i}_mean'] = np.mean(segment_flux)\n            features[f'segment{i}_std'] = np.std(segment_flux)\n            features[f'segment{i}_range'] = np.max(segment_flux) - np.min(segment_flux)\n            segment_means.append(features[f'segment{i}_mean'])\n    features['segment_mean_range'] = np.max(segment_means) - np.min(segment_means)\n    features['segment_mean_std'] = np.std(segment_means)\n    global_mean = np.mean(total_flux)\n    global_std = np.std(total_flux)\n    transit_segments = sum(1 for mean in segment_means if mean < global_mean - global_std)\n    features['num_transit_segments'] = transit_segments\n    features['transit_segment_fraction'] = transit_segments / segments\n    return features\n\ndef extract_enhanced_transit_features(signal_data, cfg, verbose=True):\n    n_frames = signal_data.shape[0]\n    signal_data = apply_adc_correction(signal_data, instrument='FGS1', adc_info_path=cfg.ADC_INFO_CSV)\n    total_flux = np.sum(signal_data, axis=(1, 2))\n    if verbose:\n        print(f\"Extracting enhanced transit features from {n_frames} frames...\")\n    features = {}\n    features.update(extract_global_flux_features(total_flux))\n    features.update(extract_rolling_statistics_features(total_flux, cfg.ROLLING_WINDOWS))\n    features.update(extract_transit_detection_features(total_flux, cfg.DETREND_WINDOW))\n    features.update(extract_frequency_features(total_flux, cfg.AUTOCORR_LAGS))\n    features.update(extract_spatial_features(signal_data))\n    features.update(extract_gradient_features(total_flux, cfg.GRADIENT_WINDOWS, cfg.N_SEGMENTS))\n    if verbose:\n        print(f\"Generated {len(features)} enhanced transit features\")\n    return features\n\ndef prepare_data_with_enhanced_fgs1(train_df, star_info_df, cfg, dataset_type='train'):\n    \"\"\"Enhanced FGS1 feature extractor\"\"\"\n    \n    if dataset_type == 'train':\n        data_path = cfg.TRAIN_DATA_PATH\n    else:\n        data_path = cfg.TEST_DATA_PATH\n\n    print(f\"Extracting enhanced FGS1 features for {len(star_info_df)} planets from {data_path}...\")\n    fgs1_features = []\n    \n    star_info_df['planet_id'] = star_info_df['planet_id'].astype(str)\n    if train_df is not None:\n        train_df['planet_id'] = train_df['planet_id'].astype(str)\n\n    for i, row in star_info_df.iterrows():\n        base_id = int(float(row['planet_id']))\n        print(f\"\\rProcessing planet {i+1}/{len(star_info_df)} (ID: {base_id})\", end='', flush=True)\n        try:\n            signal_paths = sorted(glob.glob(f\"{data_path}{base_id}/FGS1_signal_*.parquet\"))\n            for j, path in enumerate(signal_paths):\n                df = pd.read_parquet(path)\n                signal = df.values.reshape(135000, 32, 32)\n                features = extract_enhanced_transit_features(signal, cfg, verbose=False)\n                signal_id = f\"{base_id}_{j}\"\n                features['planet_id'] = signal_id\n                fgs1_features.append(features)\n        except Exception as e:\n            print(f\"\\n❌ Failed to process planet {base_id}: {e}\")\n            continue\n\n    print(\"\\n✅ Feature extraction complete.\")\n    \n    features_df = pd.DataFrame(fgs1_features)\n    features_df['planet_id'] = features_df['planet_id'].astype(str)\n    features_df = features_df.set_index('planet_id')\n    print(f\"→ Extracted features for {len(features_df)} entries with {features_df.shape[1]} columns\")\n    \n    expanded_meta = []\n    star_info_df['planet_id'] = star_info_df['planet_id'].apply(lambda x: str(int(float(x))))\n    \n    for pid in features_df.index:\n        base_id = str(int(float(pid.split(\"_\")[0])))\n        row = star_info_df[star_info_df['planet_id'] == base_id].copy()\n        if not row.empty:\n            row['planet_id'] = pid\n            expanded_meta.append(row)\n\n    star_info_df_expanded = pd.concat(expanded_meta, ignore_index=True)\n    full_df = star_info_df_expanded.set_index('planet_id').join(features_df, how='left')\n    X = full_df.select_dtypes(include=[np.number]).fillna(0).astype(np.float32)\n\n    if train_df is not None:\n        targets_df = train_df.set_index('planet_id')\n        extended_targets = []\n        for pid in X.index:\n            base_id = pid.split(\"_\")[0]\n            if base_id in targets_df.index:\n                y_row = targets_df.loc[base_id].copy()\n                y_row.name = pid\n                extended_targets.append(y_row)\n        y = pd.DataFrame(extended_targets).astype(np.float32)\n        print(f\"✅ Final shapes — X: {X.shape}, y: {y.shape}\")\n        return X, y\n    else:\n        print(f\"✅ Final shape — X_test: {X.shape}\")\n        return X\n\ndef train_validate_multioutput_cv(X, y, cfg):\n    print(f\"Training with CV: {X.shape[1]} features → {y.shape[1]} wavelengths\")\n    \n    kf = KFold(n_splits=cfg.N_SPLITS, shuffle=True, random_state=42)\n\n    all_models = [[] for _ in range(y.shape[1])]\n    all_fold_metrics = []\n\n    oof_preds = np.zeros_like(y.values, dtype=float)\n    oof_sigmas = np.zeros_like(y.values, dtype=float)\n    residuals_all = [[] for _ in range(y.shape[1])]\n\n    for fold, (train_idx, val_idx) in enumerate(kf.split(X)):\n        print(f\"\\n🔁 Fold {fold+1}/{cfg.N_SPLITS}\")\n        X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]\n        y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]\n\n        for i in tqdm(range(y.shape[1]), desc=f\"Training wavelength models for Fold {fold+1}\"):\n            model = lgb.LGBMRegressor(**cfg.LGBM_PARAMS)\n\n            y_single = y_train.iloc[:, i]\n            model.fit(X_train, y_single)\n            all_models[i].append(model)\n\n            pred_val = model.predict(X_val)\n            oof_preds[val_idx, i] = pred_val\n\n            pred_train = model.predict(X_train)\n            residuals = y_single - pred_train\n            residuals_all[i].extend(residuals.tolist())\n\n            sigma = np.std(residuals)\n            oof_sigmas[val_idx, i] = max(sigma, 1e-6)\n\n    target_uncertainties = [max(np.std(res), 1e-6) for res in residuals_all]\n    r2_scores = [r2_score(y.values[:, i], oof_preds[:, i]) for i in range(y.shape[1])]\n    rmses = [mean_squared_error(y.values[:, i], oof_preds[:, i], squared=False) for i in range(y.shape[1])]\n\n    print(\"\\n📊 CV Performance Summary:\")\n    print(f\"→ Mean R²:    {np.mean(r2_scores):.4f}\")\n    print(f\"→ Mean RMSE: {np.mean(rmses):.4f}\")\n\n    return {\n        'models': all_models,\n        'feature_columns': X.columns.tolist(),\n        'target_columns': y.columns.tolist(),\n        'oof_predictions': oof_preds,\n        'oof_uncertainties': oof_sigmas,\n        'y_true': y.values,\n        'cv_metrics': {\n            'r2_per_target': r2_scores,\n            'rmse_per_target': rmses,\n            'mean_r2': np.mean(r2_scores),\n            'mean_rmse': np.mean(rmses)\n        },\n        'target_uncertainties': target_uncertainties\n    }\n\ndef save_trained_models(model_dict, save_dir):\n    os.makedirs(save_dir, exist_ok=True)\n    for target_idx, model_list in enumerate(model_dict['models']):\n        for fold_idx, model in enumerate(model_list):\n            model_path = os.path.join(save_dir, f'model_target{target_idx}_fold{fold_idx}.pkl')\n            joblib.dump(model, model_path)\n    metadata = {\n        'feature_columns': model_dict['feature_columns'],\n        'target_columns': model_dict['target_columns'],\n        'cv_metrics': model_dict['cv_metrics'],\n        'target_uncertainties': model_dict['target_uncertainties']\n    }\n    with open(os.path.join(save_dir, 'trained_model_info.pkl'), 'wb') as f:\n        pickle.dump(metadata, f)\n    print(f\"✅ All models and metadata saved to: {save_dir}\")\n\ndef analyze_feature_importance(trained_models, X, cfg, save_file=True):\n    print(f\"\\n🔍 Feature Importance Analysis:\")\n    all_models_nested = trained_models['models']\n    all_models_flat = [m for models_per_target in all_models_nested for m in models_per_target]\n    \n    feature_importance = pd.DataFrame({\n        'feature': trained_models['feature_columns'],\n        'importance': np.mean([model.feature_importances_ for model in all_models_flat], axis=0)\n    }).sort_values('importance', ascending=False)\n    \n    print(\"Most important features:\")\n    for i, row in feature_importance.head(25).iterrows():\n        print(f\"  {i+1:2d}. {row['feature']:<35} {row['importance']:.4f}\")\n    \n    categories = {}\n    for _, row in feature_importance.iterrows():\n        category = row['feature'].split('_')[0] if '_' in row['feature'] else 'other'\n        categories.setdefault(category, []).append(row['importance'])\n    \n    print(f\"\\nFeature Category Performance:\")\n    for category, importances in sorted(categories.items(), key=lambda x: np.sum(x[1]), reverse=True):\n        total_contrib = np.sum(importances) / feature_importance['importance'].sum() * 100\n        print(f\"  {category:<15}: {len(importances):>3d} features, {total_contrib:>5.1f}% contribution\")\n    \n    zero_features = (feature_importance['importance'] == 0).sum()\n    print(f\"\\nQuick Stats:\")\n    print(f\"  Zero importance features: {zero_features}\")\n    print(f\"  Top feature: {feature_importance.iloc[0]['feature']}\")\n    \n    if save_file:\n        feature_importance.to_csv(cfg.IMPORTANCE_FILE, index=False)\n        print(f\"\\n💾 Saved feature importance to: {cfg.IMPORTANCE_FILE}\")\n    \n    return feature_importance\n\ndef predict_with_uncertainty(models, X, fixed_uncertainty=None, target_uncertainties=None):\n    predictions = []\n    uncertainties = []\n    is_cv = isinstance(models[0], list)\n    n_targets = len(models)\n\n    for i in range(n_targets):\n        model_group = models[i] if is_cv else [models[i]]\n        preds = []\n        for model in model_group:\n            preds.append(model.predict(X))\n        \n        pred_avg = np.mean(preds, axis=0)\n        predictions.append(pred_avg)\n\n        if fixed_uncertainty is not None:\n            unc = fixed_uncertainty\n        elif target_uncertainties is not None:\n            unc = target_uncertainties[i]\n        else:\n            unc = 0.01  # fallback\n\n        unc_array = np.full_like(pred_avg, max(unc, 1e-6))\n        uncertainties.append(unc_array)\n\n    y_pred = np.column_stack(predictions)\n    sigma_pred = np.column_stack(uncertainties)\n    \n    pred_df = pd.DataFrame(y_pred).apply(pd.to_numeric, errors='coerce').fillna(0).clip(lower=0)\n    sigma_df = pd.DataFrame(sigma_pred).apply(pd.to_numeric, errors='coerce').fillna(1e-6).clip(lower=1e-15)\n    \n    return pred_df.values, sigma_df.values\n\ndef score(solution, submission, row_id_column_name, naive_mean, naive_sigma,\n          fsg_sigma_true, airs_sigma_true, fgs_weight):\n    y_true = solution.drop(columns=[row_id_column_name])\n    y_pred = submission[[col for col in submission.columns if not col.endswith('_std') and col != row_id_column_name]]\n    sigma = submission[[col for col in submission.columns if col.endswith('_std')]]\n    log_likelihoods = -0.5 * np.log(2 * np.pi * sigma.values ** 2) - ((y_true.values - y_pred.values) ** 2) / (2 * sigma.values ** 2)\n    gll = np.mean(log_likelihoods)\n    return np.clip((gll + 10) / 10, 0, 1)\n\ndef evaluate_with_gll_metric(y_true, y_pred, sigma_pred, naive_mean, naive_sigma, cfg):\n    if not isinstance(y_true, pd.DataFrame):\n        y_true = pd.DataFrame(y_true)\n    y_true = y_true.reset_index(drop=True)\n    n_samples, n_waves = y_pred.shape\n    sigma_pred = np.clip(sigma_pred, 1e-15, None)\n    solution_df = y_true.copy()\n    solution_df['row_id'] = np.arange(n_samples)\n    submission_df = pd.DataFrame()\n    for i in range(n_waves):\n        submission_df[f'wavelength_{i}'] = y_pred[:, i]\n    for i in range(n_waves):\n        submission_df[f'wavelength_{i}_std'] = sigma_pred[:, i]\n    submission_df['row_id'] = np.arange(n_samples)\n    submission_df = submission_df.clip(lower=1e-15)\n    \n    gll_score = score(\n        solution=solution_df, submission=submission_df,\n        row_id_column_name='row_id', naive_mean=naive_mean, naive_sigma=naive_sigma,\n        fsg_sigma_true=cfg.FSG_SIGMA_TRUE, airs_sigma_true=cfg.AIRS_SIGMA_TRUE, fgs_weight=cfg.FGS_WEIGHT\n    )\n    return gll_score\n\ndef fast_gll_score_numpy(y_true, y_pred, sigma_pred, naive_mean, naive_sigma, cfg):\n    sigma_pred = np.clip(sigma_pred, 1e-15, None)\n    n_samples, n_waves = sigma_pred.shape\n    sigma_true = np.append([cfg.FSG_SIGMA_TRUE], np.full(n_waves - 1, cfg.AIRS_SIGMA_TRUE))\n    sigma_true = np.tile(sigma_true, (n_samples, 1))\n    weights = np.append([cfg.FGS_WEIGHT], np.ones(n_waves - 1))\n    weights = np.tile(weights, (n_samples, 1))\n    gll_pred = norm.logpdf(y_true, loc=y_pred, scale=sigma_pred)\n    gll_true = norm.logpdf(y_true, loc=y_true, scale=sigma_true)\n    gll_naive = norm.logpdf(y_true, loc=naive_mean, scale=naive_sigma)\n    ind_scores = (gll_pred - gll_naive) / (gll_true - gll_naive + 1e-9)\n    final_score = np.average(ind_scores, weights=weights)\n    return float(np.clip(final_score, 0.0, 1.0))\n\ndef optimize_sigma_per_wavelength(y_true, y_pred, sigma_pred, naive_mean, naive_sigma, cfg):\n    n_waves = sigma_pred.shape[1]\n    best_scales = np.ones(n_waves)\n    print(f\"⚡ Optimizing {n_waves} sigma scalers with fast NumPy GLL metric...\\n\")\n    for i in tqdm(range(n_waves), desc=\"Optimizing wavelengths\", unit=\"λ\"):\n        def objective(scale):\n            sigma_scaled = sigma_pred.copy()\n            sigma_scaled[:, i] *= scale\n            return -fast_gll_score_numpy(y_true, y_pred, sigma_scaled, naive_mean, naive_sigma, cfg)\n        \n        result = minimize_scalar(objective, bounds=cfg.SIGMA_OPTIMIZE_BOUNDS, method='bounded')\n        best_scales[i] = result.x\n    print(\"✅ Optimization complete.\")\n    return best_scales\n\ndef create_submission(trained_models, test_star_info_df, cfg, sigma_scalers=None):\n    X_test = prepare_data_with_enhanced_fgs1(\n        train_df=None,\n        star_info_df=test_star_info_df,\n        cfg=cfg,\n        dataset_type='test'\n    )\n    \n    X_test_aligned = X_test.reindex(columns=trained_models['feature_columns'], fill_value=0)\n\n    y_pred, sigma_pred = predict_with_uncertainty(\n        trained_models['models'],\n        X_test_aligned,\n        target_uncertainties=trained_models['target_uncertainties']\n    )\n\n    if sigma_scalers is not None:\n        sigma_pred = sigma_pred * sigma_scalers\n\n    base_ids = [idx.split(\"_\")[0] for idx in X_test_aligned.index]\n    wl_cols = trained_models['target_columns']\n    sigma_cols = [f\"sigma_{i+1}\" for i in range(y_pred.shape[1])]\n\n    pred_df = pd.DataFrame(y_pred, index=base_ids, columns=wl_cols)\n    sigma_df = pd.DataFrame(sigma_pred, index=base_ids, columns=sigma_cols)\n\n    pred_mean = pred_df.groupby(pred_df.index).mean()\n    sigma_mean = sigma_df.groupby(sigma_df.index).mean()\n\n    submission_df = pd.concat([pred_mean, sigma_mean], axis=1).reset_index()\n    submission_df = submission_df.rename(columns={'index': 'planet_id'})\n    submission_df['planet_id'] = submission_df['planet_id'].astype(int)\n\n    submission_df.to_csv(cfg.SUBMISSION_FILE, index=False, float_format='%.5f')\n    print(f\"✅ Submission saved to: {cfg.SUBMISSION_FILE}\")\n    \n    return submission_df\n\nmain_start_time = time.time()\n\ncfg = Config()\n\ntrain_df, train_star_info_df, test_star_info_df = load_initial_data(cfg)\n\nif cfg.QUICK_RUN:\n    print(\"--- Running Demo Cell ---\")\n    demo_start = time.time()\n    demo_signal, demo_pid = load_demo_data(cfg, train_star_info_df)\n    demo_flux = np.sum(apply_adc_correction(demo_signal, 'FGS1', cfg.ADC_INFO_CSV), axis=(1, 2))\n    \n    visualize_transit_data(demo_signal, demo_flux, demo_pid)\n    \n    all_features = extract_enhanced_transit_features(demo_signal, cfg, verbose=True)\n    print(f\"--- Demo Cell Finished ({time.time() - demo_start:.2f}s) ---\")\n\nprint(\"--- Preparing Training Data ---\")\nprep_start = time.time()\nX, y = prepare_data_with_enhanced_fgs1(train_df, train_star_info_df, cfg, dataset_type='train')\nprint(f\"--- Training Data Prepared ({time.time() - prep_start:.2f}s) ---\")\n\nprint(\"--- Training Model ---\")\ntrain_start = time.time()\ntrained_models_cv = train_validate_multioutput_cv(X, y, cfg)\nprint(f\"--- Model Training Finished ({time.time() - train_start:.2f}s) ---\")\n\nprint(\"--- Saving & Analyzing Model ---\")\nsave_trained_models(trained_models_cv, save_dir=cfg.MODEL_DIR)\nfeature_importance_df = analyze_feature_importance(trained_models_cv, X, cfg)\n\nprint(\"--- Running Dummy GLL Score ---\")\ny_true_oof = trained_models_cv['y_true']\ny_pred_oof = trained_models_cv['oof_predictions']\nsigma_pred_oof = trained_models_cv['oof_uncertainties']\nnaive_mean = y_true_oof.mean()\nnaive_sigma = y_true_oof.std()\n\nbest_scales = optimize_sigma_per_wavelength(\n    y_true_oof, y_pred_oof, sigma_pred_oof,\n    naive_mean, naive_sigma, cfg\n)\n\nfinal_score = fast_gll_score_numpy(\n    y_true_oof, y_pred_oof, sigma_pred_oof * best_scales,\n    naive_mean, naive_sigma, cfg\n)\nprint(f\"🎯 Calibrated OOF GLL Score: {final_score:.6f}\")\n\nprint(\"--- Creating Submission ---\")\nsub_start = time.time()\nsubmission = create_submission(\n    trained_models_cv,\n    test_star_info_df,\n    cfg,\n    sigma_scalers=best_scales\n)\nprint(f\"--- Submission Created ({time.time() - sub_start:.2f}s) ---\")\n\nprint(\"\\n--- 🚀 Pipeline Finished ---\")\ndisplay(submission.head())\nprint(f\"Total runtime: {(time.time() - main_start_time) / 60:.2f} minutes\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T19:45:05.402919Z","iopub.execute_input":"2025-11-14T19:45:05.403405Z","iopub.status.idle":"2025-11-14T21:16:36.232101Z","shell.execute_reply.started":"2025-11-14T19:45:05.403377Z","shell.execute_reply":"2025-11-14T21:16:36.230054Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ... (Continuing from the end of the original script's main execution block)\n\n# --- 6. Save & Analyze Model ---\nprint(\"--- Saving & Analyzing Model ---\")\nsave_trained_models(trained_models_cv, save_dir=cfg.MODEL_DIR)\nfeature_importance_df = analyze_feature_importance(trained_models_cv, X, cfg)\n\n## ===================================================================\n## Data Visualization and Analysis\n## ===================================================================\n\n# Define the visualization functions here (or ensure they are defined globally)\ndef plot_feature_distributions(X, feature_importance_df, n_features=9):\n    \"\"\"Plots the distribution of the top N features.\"\"\"\n    plt.style.use('seaborn-v0_8-whitegrid')\n    \n    # Select the top N features\n    top_features = feature_importance_df.head(n_features)['feature'].tolist()\n    \n    # Calculate the grid size\n    n_cols = 3\n    n_rows = (n_features + n_cols - 1) // n_cols\n\n    fig, axes = plt.subplots(n_rows, n_cols, figsize=(16, 4 * n_rows))\n    axes = axes.flatten()\n\n    fig.suptitle(f'Distribution of Top {n_features} Features', fontsize=16, y=1.02)\n    \n    for i, feature in enumerate(top_features):\n        ax = axes[i]\n        data = X[feature].dropna()\n        \n        # Plot histogram\n        data.hist(ax=ax, bins=50, color='skyblue', edgecolor='black', alpha=0.7, density=True)\n        \n        # Fit a normal distribution for comparison\n        mu, std = norm.fit(data)\n        xmin, xmax = ax.get_xlim()\n        x = np.linspace(xmin, xmax, 100)\n        p = norm.pdf(x, mu, std)\n        ax.plot(x, p, 'r', linewidth=2, label=f'Norm Fit (μ={mu:.2f}, σ={std:.2f})')\n        \n        ax.set_title(feature, fontsize=10)\n        ax.set_xlabel('Value', fontsize=8)\n        ax.set_ylabel('Density', fontsize=8)\n        ax.tick_params(axis='both', which='major', labelsize=7)\n        ax.legend(fontsize=7, loc='upper right')\n\n    # Hide unused subplots\n    for i in range(n_features, len(axes)):\n        fig.delaxes(axes[i])\n        \n    plt.tight_layout()\n    plt.show()\n\ndef plot_oof_predictions(y_true, y_pred, target_columns, n_targets=4):\n    \"\"\"Plots OOF predictions vs True values for the first N target wavelengths.\"\"\"\n    plt.style.use('seaborn-v0_8-whitegrid')\n    \n    # Select the first N targets\n    n_targets = min(n_targets, y_true.shape[1])\n    \n    fig, axes = plt.subplots(1, n_targets, figsize=(18, 5))\n    \n    fig.suptitle(f'OOF Predictions vs. True Values (First {n_targets} Wavelengths)', fontsize=16, y=1.05)\n    \n    for i in range(n_targets):\n        ax = axes[i]\n        true_values = y_true[:, i]\n        pred_values = y_pred[:, i]\n        \n        # Scatter plot\n        ax.scatter(true_values, pred_values, alpha=0.5, s=5, color='darkblue')\n        \n        # Add the ideal y=x line\n        min_val = min(true_values.min(), pred_values.min())\n        max_val = max(true_values.max(), pred_values.max())\n        ax.plot([min_val, max_val], [min_val, max_val], 'r--', linewidth=2, label='Ideal Fit (y=x)')\n        \n        # Calculate R2 for the plot title\n        r2 = r2_score(true_values, pred_values)\n        \n        ax.set_title(f'{target_columns[i]} (R²: {r2:.4f})', fontsize=12)\n        ax.set_xlabel('True Value', fontsize=10)\n        ax.set_ylabel('OOF Prediction', fontsize=10)\n        ax.tick_params(axis='both', which='major', labelsize=8)\n        ax.set_aspect('equal', adjustable='box')\n        ax.legend(fontsize=8)\n        \n    plt.tight_layout()\n    plt.show()\n\nprint(\"\\n--- 🖼️ Visualization Stage ---\")\n\n# 1. Plot feature distributions (X and feature_importance_df are now defined)\nplot_feature_distributions(X, feature_importance_df, n_features=9)\n\n# 2. Plot OOF predictions vs True values (trained_models_cv is now defined)\nplot_oof_predictions(\n    y_true=trained_models_cv['y_true'],\n    y_pred=trained_models_cv['oof_predictions'],\n    target_columns=trained_models_cv['target_columns'],\n    n_targets=4 # Plotting the first 4 wavelengths\n)\n\n# --- 7. Dummy Score (for demonstration) ---\nprint(\"--- Running Dummy GLL Score ---\")\n# ... (rest of the script)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-11-14T21:16:36.234033Z","iopub.execute_input":"2025-11-14T21:16:36.23442Z","iopub.status.idle":"2025-11-14T21:16:50.026012Z","shell.execute_reply.started":"2025-11-14T21:16:36.234373Z","shell.execute_reply":"2025-11-14T21:16:50.024804Z"}},"outputs":[],"execution_count":null}]}