{"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":12846694,"sourceType":"competition"}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import pandas as pd\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport scipy.stats\nfrom tqdm import tqdm\nimport pickle\nfrom sklearn.model_selection import KFold\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import mean_squared_error\nfrom catboost import CatBoostRegressor\nimport os\nfrom glob import glob\nfrom collections import defaultdict\nimport gc\nimport time # Import time module for timing\n\n# --- Configuration and Paths ---\nDATA_DIR = '/kaggle/input/ariel-data-challenge-2025'\nTRAIN_DIR = os.path.join(DATA_DIR, 'train')\nTEST_DIR = os.path.join(DATA_DIR, 'test')\nSAMPLE_SUBMISSION_PATH = os.path.join(DATA_DIR, 'sample_submission.csv')\n\n# --- Load Base Data ---\ntrain_df = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'))\nwavelengths_df = pd.read_csv(os.path.join(DATA_DIR, 'wavelengths.csv'))\ntrain_star_info = pd.read_csv(os.path.join(DATA_DIR, 'train_star_info.csv'))\ntrain_adc_info = pd.read_csv(os.path.join(DATA_DIR, 'adc_info.csv'))\ntrain_labels = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'), index_col='planet_id')\n\n# Extract target column names (flux_0 to flux_355)\ntarget_cols = [col for col in train_df.columns if col != 'planet_id']\nNUM_WAVELENGTHS = len(target_cols) # This should be 356\n\n# --- Define FGS1 and AIRS-Ch0 channel indices and weights ---\nFGS1_CHANNEL_IDX = 0 # Corresponds to flux_0\nAIRS_CH0_CHANNEL_INDICES = list(range(1, NUM_WAVELENGTHS)) # Corresponds to flux_1 to flux_355\n\n# Weights based on competition description\nWEIGHT_FGS1 = 2 * (0.2 / 1) # 0.4\nWEIGHT_AIRS_CH0_PER_POINT = 1.95 / 282 # ~0.0069\n\n# Ideal uncertainty values (from competition description)\nIDEAL_UNC_FGS1_PPM = 1 # 1 PPM\nIDEAL_UNC_AIRS_PPM = 10 # 10 PPM\n\nprint(\"Setup Complete. Data Loaded.\")\n\n# --- Custom GLL Metric Function (DEFINITIVELY CORRECTED) ---\ndef calculate_gll(y_true, mu_pred, sigma_pred, wavelengths_df_unused): # Renamed parameter for clarity\n    \"\"\"\n    Calculates the Gaussian Log-likelihood (GLL) score for the Ariel Data Challenge.\n    (MODIFIED to infer instrument from channel index, NOT wavelengths_df_for_gll)\n\n    Parameters:\n    y_true (pd.DataFrame or np.array): Ground truth flux values (shape: num_planets, num_wavelengths).\n    mu_pred (pd.DataFrame or np.array): Predicted mean flux values (shape: num_planets, num_wavelengths).\n    sigma_pred (pd.DataFrame or np.array): Predicted uncertainty values (shape: num_planets, num_wavelengths).\n    wavelengths_df_unused (pd.DataFrame): The wavelengths.csv dataframe (this parameter is now truly UNUSED\n                                          for instrument lookup to avoid KeyError).\n\n    Returns:\n    float: The total GLL value (L).\n    \"\"\"\n    y_true = np.asarray(y_true)\n    mu_pred = np.asarray(mu_pred)\n    sigma_pred = np.asarray(sigma_pred)\n\n    # Ensure all inputs have the same shape\n    if not (y_true.shape == mu_pred.shape == sigma_pred.shape):\n        raise ValueError(\"y_true, mu_pred, and sigma_pred must have the same shape.\")\n\n    # Ensure sigma_pred is strictly positive\n    sigma_pred = np.maximum(sigma_pred, 1e-9) # Prevent log(0) or division by zero\n\n    # Calculate log-likelihood for each point\n    log_likelihood_per_point = -0.5 * np.log(2 * np.pi * sigma_pred**2) - 0.5 * ((y_true - mu_pred) / sigma_pred)**2\n\n    total_gll = 0.0\n    \n    # Apply weights based on instrument, now determined by the column index (i)\n    # The problematic line 'instrument = wavelengths_df_for_gll.loc[i, 'instrument']' IS REMOVED HERE.\n    for i in range(y_true.shape[1]): # Iterate through each column (wavelength)\n        \n        # Determine instrument based on pre-defined global channel indices\n        if i == FGS1_CHANNEL_IDX: # Check if it's the FGS1 channel (flux_0)\n            weight = WEIGHT_FGS1\n        elif i in AIRS_CH0_CHANNEL_INDICES: # Check if it's one of the AIRS-CH0 channels (flux_1 to flux_355)\n            weight = WEIGHT_AIRS_CH0_PER_POINT\n        else:\n            # This case should ideally not be reached if NUM_WAVELENGTHS matches expected channels\n            print(f\"Warning: Unexpected wavelength index {i} found. Assigning weight 0.\")\n            weight = 0 \n\n        # Sum likelihoods for this wavelength across all planets, then apply weight\n        total_gll += np.sum(log_likelihood_per_point[:, i]) * weight\n\n    return total_gll\n\n# --- Function to calculate the final competition score ---\ndef calculate_competition_score(L_user, L_ideal, L_ref):\n    \"\"\"\n    Calculates the final competition score.\n\n    Parameters:\n    L_user (float): GLL value from user's submission.\n    L_ideal (float): GLL value for the ideal case.\n    L_ref (float): GLL value for the reference case.\n\n    Returns:\n    float: The competition score, clipped to [0, 1].\n    \"\"\"\n    score = (L_user - L_ref) / (L_ideal - L_ref)\n    return np.clip(score, 0, 1)\n\nprint(\"Custom GLL Metric function defined.\")\n\n# --- Feature Engineering Functions (now defined within the main script) ---\n\ndef f_read_and_preprocess(dataset_path, planet_ids, band_name=\"FGS1\"):\n    \"\"\"Read FGS1 files and extract comprehensive time series features.\"\"\"\n    extracted_features = []\n    for planet_id in tqdm(planet_ids, desc=f\"Processing {band_name} signals\"):\n        try:\n            signal_files = glob(os.path.join(dataset_path, str(planet_id), f\"{band_name}_signal_*.parquet\"))\n            \n            if not signal_files:\n                features_dict = {f'{band_name}_mean': np.nan, f'{band_name}_std': np.nan,\n                                 f'{band_name}_skew': np.nan, f'{band_name}_kurt': np.nan,\n                                 f'{band_name}_min': np.nan, f'{band_name}_max': np.nan}\n                extracted_features.append({'planet_id': planet_id, **features_dict})\n                continue\n            \n            f_signal = pl.read_parquet(signal_files[0])\n            \n            pixel_count = 32 * 32\n                \n            mean_signal = f_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy().flatten() / pixel_count\n            net_signal = mean_signal[1::2] - mean_signal[0::2]\n\n            features_dict = {\n                f'{band_name}_mean': np.mean(net_signal),\n                f'{band_name}_std': np.std(net_signal),\n                f'{band_name}_skew': scipy.stats.skew(net_signal),\n                f'{band_name}_kurt': scipy.stats.kurtosis(net_signal),\n                f'{band_name}_min': np.min(net_signal),\n                f'{band_name}_max': np.max(net_signal),\n            }\n            extracted_features.append({'planet_id': planet_id, **features_dict})\n        except Exception as e:\n            print(f\"Error processing planet_id {planet_id} for {band_name}: {e}\")\n            features_dict = {f'{band_name}_mean': np.nan, f'{band_name}_std': np.nan,\n                             f'{band_name}_skew': np.nan, f'{band_name}_kurt': np.nan,\n                             f'{band_name}_min': np.nan, f'{band_name}_max': np.nan}\n            extracted_features.append({'planet_id': planet_id, **features_dict})\n    return pd.DataFrame(extracted_features).set_index('planet_id')\n\ndef a_read_and_preprocess(dataset_path, planet_ids, band_name=\"AIRS-CH0\"):\n    \"\"\"Read AIRS-CH0 files and extract comprehensive time series features.\"\"\"\n    extracted_features = []\n    for planet_id in tqdm(planet_ids, desc=f\"Processing {band_name} signals\"):\n        try:\n            signal_files = glob(os.path.join(dataset_path, str(planet_id), f\"{band_name}_signal_*.parquet\"))\n            \n            if not signal_files:\n                features_dict = {f'{band_name}_mean': np.nan, f'{band_name}_std': np.nan,\n                                 f'{band_name}_skew': np.nan, f'{band_name}_kurt': np.nan,\n                                 f'{band_name}_min': np.nan, f'{band_name}_max': np.nan}\n                extracted_features.append({'planet_id': planet_id, **features_dict})\n                continue\n\n            a_signal = pl.read_parquet(signal_files[0])\n            pixel_count = 32 * 356\n            mean_signal = a_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy().flatten() / pixel_count\n            net_signal = mean_signal[1::2] - mean_signal[0::2]\n            \n            features_dict = {\n                f'{band_name}_mean': np.mean(net_signal),\n                f'{band_name}_std': np.std(net_signal),\n                f'{band_name}_skew': scipy.stats.skew(net_signal),\n                f'{band_name}_kurt': scipy.stats.kurtosis(net_signal),\n                f'{band_name}_min': np.min(net_signal),\n                f'{band_name}_max': np.max(net_signal),\n            }\n            extracted_features.append({'planet_id': planet_id, **features_dict})\n        except Exception as e:\n            print(f\"Error processing planet_id {planet_id} for {band_name}: {e}\")\n            features_dict = {f'{band_name}_mean': np.nan, f'{band_name}_std': np.nan,\n                             f'{band_name}_skew': np.nan, f'{band_name}_kurt': np.nan,\n                             f'{band_name}_min': np.nan, f'{band_name}_max': np.nan}\n            extracted_features.append({'planet_id': planet_id, **features_dict})\n            \n    return pd.DataFrame(extracted_features).set_index('planet_id')\n\n\n# --- 1. Visualization & Exploratory Data Analysis (EDA) ---\nprint(\"\\n--- 1. Visualization & Exploratory Data Analysis (EDA) ---\")\n\n# --- Basic Train Set Stats ---\nprint(\"🪐 Number of training planets:\", train_df.shape[0])\nprint(\"📈 Number of target labels (wavelengths):\", train_df.shape[1] - 1)\nprint(\"🔬 Length of wavelength grid:\", wavelengths_df.shape[0])\n\n# --- Target Stats (per flux column) ---\nflux_summary = train_df[target_cols].agg(['min', 'max', 'mean', 'std']).T\nprint(\"\\n📊 Flux value summary (first 5 rows):\")\nprint(flux_summary.head())\n\n# --- Unique Stars ---\nif 'planet_id' in train_star_info.columns:\n    num_stars = train_star_info.drop(columns='planet_id').drop_duplicates().shape[0]\nelse:\n    num_stars = train_star_info.drop_duplicates().shape[0]\nprint(\"\\n🌟 Number of unique stars in training:\", num_stars)\n\n# --- Planets with Multiple Observations ---\nobs_counts = defaultdict(int)\ntrain_planets_dirs = [d for d in os.listdir(TRAIN_DIR) if os.path.isdir(os.path.join(TRAIN_DIR, d))]\nfor pid_str in train_planets_dirs:\n    air_obs = glob(os.path.join(TRAIN_DIR, pid_str, \"AIRS-CH0_signal_*.parquet\"))\n    obs_counts[pid_str] = len(air_obs)\nmulti_obs = {pid: count for pid, count in obs_counts.items() if count > 1}\nprint(\"\\n🔁 Planets with multiple observations:\", len(multi_obs))\n\n# --- Check Calibration File Coverage ---\nmissing_calibs = []\nexpected_calibs = {\"dark\", \"dead\", \"flat\", \"linear_corr\", \"read\"}\nfor pid_str in train_planets_dirs:\n    for band in [\"AIRS-CH0\", \"FGS1\"]:\n        calib_path = os.path.join(TRAIN_DIR, pid_str, f\"{band}_calibration\")\n        calib_files = set()\n        if os.path.exists(calib_path):\n            calib_files = {os.path.splitext(f)[0] for f in os.listdir(calib_path)}\n        missing = expected_calibs - calib_files\n        if missing:\n            missing_calibs.append((pid_str, band, missing))\nprint(\"\\n🧪 Planets missing calibration files:\", len(missing_calibs))\nif missing_calibs:\n    print(\"   Example:\", missing_calibs[0])\n\n# --- Distribution of Observations Per Planet ---\nobs_distribution = pd.Series(list(obs_counts.values())).value_counts().sort_index()\nprint(\"\\n🗂 Observation count distribution per planet (AIR-CH0):\")\nprint(obs_distribution)\n\n# --- Planet-Star Uniqueness Check ---\nmerged_star_info = pd.merge(train_df[['planet_id']], train_star_info, on='planet_id', how='left')\nunique_links = merged_star_info[['planet_id'] + [col for col in train_star_info.columns if col != 'planet_id']].drop_duplicates()\nprint(\"\\n🔗 Unique planet-star mappings:\", unique_links.shape[0])\n\n# --- Visualizing FGS1 Images ---\nprint(\"\\n--- Visualizing FGS1 Images ---\")\nplanet_id_fgs = 1010375142 # Example ID\ntry:\n    f_signal_ex_pl = pl.read_parquet(os.path.join(TRAIN_DIR, str(planet_id_fgs), 'FGS1_signal_0.parquet'))\n    # Convert Polars DataFrame to Pandas DataFrame for .iloc access\n    f_signal_ex = f_signal_ex_pl.to_pandas()\n\n    _, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\n    sns.heatmap(f_signal_ex.iloc[0].values.reshape(32, 32), ax=ax1, vmin=0, vmax=52000, cmap='viridis')\n    ax1.set_aspect('equal')\n    ax1.set_title(f'FGS1 Image Frame 0 for Planet {planet_id_fgs}')\n    sns.heatmap(f_signal_ex.iloc[1].values.reshape(32, 32), ax=ax2, vmin=0, vmax=52000, cmap='viridis')\n    ax2.set_aspect('equal')\n    ax2.set_title(f'FGS1 Image Frame 1 for Planet {planet_id_fgs}')\n    plt.suptitle('A pair of FGS1 Images')\n    plt.show()\nexcept Exception as e:\n    print(f\"Could not load FGS1 example for visualization: {e}\")\n\n\n# --- Visualizing FGS1 Time Series ---\nprint(\"\\n--- Visualizing FGS1 Time Series ---\")\nplanet_id_strong_signal = 1048114509 # Planet with strong signal\nplanet_id_weak_signal = 1240764363 # Planet with weak signal\n\n_, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, figsize=(14, 8))\n\n# Strong signal planet\ntry:\n    f_signal_strong_pl = pl.read_parquet(os.path.join(TRAIN_DIR, str(planet_id_strong_signal), 'FGS1_signal_0.parquet'))\n    mean_signal_strong = f_signal_strong_pl.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy().flatten() / (32*32)\n    net_signal_strong = mean_signal_strong[1::2] - mean_signal_strong[0::2]\n    cum_signal_strong = net_signal_strong.cumsum()\n    window=800 # Define window for smoothing\n    smooth_signal_strong = (cum_signal_strong[window:] - cum_signal_strong[:-window]) / window\n\n    ax1.set_title(f'FGS1: Raw Signal (Planet {planet_id_strong_signal})')\n    ax1.plot(net_signal_strong, label='raw signal', alpha=0.7)\n    ax1.legend()\n    ax3.set_title(f'FGS1: Smoothed Signal (Planet {planet_id_strong_signal})')\n    ax3.plot(smooth_signal_strong, color='c', label='smoothened signal')\n    ax3.legend()\n    ax3.set_xlabel('time step')\n    for time_step in [20500, 23500, 44000, 47000]: # Example transit timings\n        ax3.axvline(time_step, color='gray', linestyle='--', alpha=0.6)\nexcept Exception as e:\n    print(f\"Could not load FGS1 strong signal example for visualization: {e}\")\n\n# Weak signal planet\ntry:\n    f_signal_weak_pl = pl.read_parquet(os.path.join(TRAIN_DIR, str(planet_id_weak_signal), 'FGS1_signal_0.parquet'))\n    mean_signal_weak = f_signal_weak_pl.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy().flatten() / (32*32)\n    net_signal_weak = mean_signal_weak[1::2] - mean_signal_weak[0::2]\n    cum_signal_weak = net_signal_weak.cumsum()\n    window=800 # Define window for smoothing\n    smooth_signal_weak = (cum_signal_weak[window:] - cum_signal_weak[:-window]) / window\n\n    ax2.set_title(f'FGS1: Raw Signal (Planet {planet_id_weak_signal})')\n    ax2.plot(net_signal_weak, label='raw signal', alpha=0.7)\n    ax2.legend()\n    ax4.set_title(f'FGS1: Smoothed Signal (Planet {planet_id_weak_signal})')\n    ax4.plot(smooth_signal_weak, color='c', label='smoothened signal')\n    ax4.legend()\n    ax4.set_xlabel('time step')\n    for time_step in [20500, 23500, 44000, 47000]: # Example transit timings\n        ax4.axvline(time_step, color='gray', linestyle='--', alpha=0.6)\nexcept Exception as e:\n    print(f\"Could not load FGS1 weak signal example for visualization: {e}\")\n\nplt.suptitle('FGS1 Time Series Analysis', y=1.02)\nplt.tight_layout(rect=[0, 0, 1, 0.98])\nplt.show()\n\n\n# --- Visualizing AIRS-CH0 Data ---\nprint(\"\\n--- Visualizing AIRS-CH0 Data ---\")\nplanet_id_airs = 1240764363 # Example ID\ntry:\n    a_signal_ex_pl = pl.read_parquet(os.path.join(TRAIN_DIR, str(planet_id_airs), 'AIRS-CH0_signal_0.parquet'))\n    a_signal_np = a_signal_ex_pl.to_numpy().reshape(a_signal_ex_pl.shape[0], 32, 356)\n\n    plt.figure(figsize=(12, 6))\n    sns.heatmap(a_signal_np[100], cmap='viridis') # Displaying an arbitrary time slice (e.g., 100th frame)\n    plt.ylabel('Spatial Dimension (Pixel Row)')\n    plt.xlabel('Wavelength Dimension (Pixel Column)')\n    plt.title(f'AIRS-CH0 Single Frame (Time Step 100) for Planet {planet_id_airs}')\n    plt.show()\n\n    mean_signal_airs = a_signal_ex_pl.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy().flatten() / (32*356)\n    net_signal_airs = mean_signal_airs[1::2] - mean_signal_airs[0::2]\n    cum_signal_airs = net_signal_airs.cumsum()\n    window_airs=80 # Smaller window for AIRS due to fewer total steps\n    smooth_signal_airs = (cum_signal_airs[window_airs:] - cum_signal_airs[:-window_airs]) / window_airs\n\n    _, (ax1_airs, ax2_airs) = plt.subplots(2, 1, sharex=True, figsize=(12, 6))\n    ax1_airs.plot(net_signal_airs, label='raw net signal', alpha=0.7)\n    ax1_airs.legend()\n    ax1_airs.set_title(f'AIRS-CH0: Raw Net Signal (Planet {planet_id_airs})')\n    ax2_airs.plot(smooth_signal_airs, color='c', label='smoothened net signal')\n    ax2_airs.legend()\n    ax2_airs.set_xlabel('Time Step (Paired Frames)')\n    ax2_airs.set_title(f'AIRS-CH0: Smoothed Net Signal (Planet {planet_id_airs})')\n\n    fgs1_total_frames_example = 135000\n    airs0_total_frames_example = 11250\n    scaling_factor = (airs0_total_frames_example / 2) / (fgs1_total_frames_example / 2) # Ratio of net steps\n    for time_step in [20500, 23500, 44000, 47000]:\n        ax2_airs.axvline(time_step * scaling_factor, color='gray', linestyle='--', alpha=0.6)\n    plt.suptitle('AIRS-CH0 Time Series Analysis', y=1.02)\n    plt.tight_layout(rect=[0, 0, 1, 0.98])\n    plt.show()\nexcept Exception as e:\n    print(f\"Could not load AIRS-CH0 example for visualization: {e}\")\n\nprint(\"\\nEDA and Visualization Complete.\")\n\n# --- 2. Feature Engineering & Preprocessing ---\nprint(\"\\n--- 2. Feature Engineering & Preprocessing ---\")\n\nprint(\"\\n--- Running FGS1 Preprocessing for Training Data ---\")\nstart_time = time.time()\nf_train_features = f_read_and_preprocess(TRAIN_DIR, train_labels.index.tolist(), band_name=\"FGS1\")\nend_time = time.time()\nprint(f\"FGS1 Preprocessing took: {end_time - start_time:.2f} seconds\")\n\nprint(\"\\n--- Running AIRS-CH0 Preprocessing for Training Data ---\")\nstart_time = time.time()\na_train_features = a_read_and_preprocess(TRAIN_DIR, train_labels.index.tolist(), band_name=\"AIRS-CH0\")\nend_time = time.time()\nprint(f\"AIRS-CH0 Preprocessing took: {end_time - start_time:.2f} seconds\")\n\n# Merge features and target labels\nX_train = f_train_features.merge(a_train_features, left_index=True, right_index=True, how='left')\ny_train = train_labels.copy()\n\n# --- Additional Feature Engineering (from star_info and adc_info) ---\n# Merge train_star_info using 'planet_id' (this is correct as star_info has planet_id)\nX_train = X_train.merge(train_star_info.set_index('planet_id'), left_index=True, right_index=True, how='left')\n\n# CORRECT WAY TO ADD ADC_INFO FEATURES:\n# Since train_adc_info contains global ADC settings and NO 'planet_id' column,\n# we add them as new columns to every row of X_train by broadcasting.\nif not train_adc_info.empty:\n    # Get the values from the first (and likely only) row of train_adc_info\n    # This assumes adc_info has only one row of global settings.\n    adc_features_values = train_adc_info.iloc[0].to_dict()\n\n    # Add these values as new columns to X_train\n    for col_name, value in adc_features_values.items():\n        X_train[col_name] = value\n    print(\"ADC info features added to X_train.\")\nelse:\n    print(\"Warning: train_adc_info is empty. ADC features will not be added to X_train.\")\n\n\n# Handle missing values (this should be done AFTER all merges/additions)\n# This will fill NaNs that might arise from missing signal files or star_info.\nX_train = X_train.fillna(0)\n\nprint(\"\\n--- Processed Training Features (X_train) Head ---\")\nprint(X_train.head())\nprint(f\"X_train shape: {X_train.shape}\")\nprint(f\"y_train shape: {y_train.shape}\")\n\nprint(\"\\nFeature Engineering and Preprocessing Complete.\")\n\n# --- 3. Model Training & Fitting ---\nprint(\"\\n--- 3. Model Training & Fitting ---\")\n\n# Define models for MEAN prediction\nridge_mean_model = Ridge(random_state=42)\ncat_mean_model = CatBoostRegressor(\n    iterations=500, learning_rate=0.05, depth=6,\n    loss_function='RMSE', eval_metric='RMSE',\n    random_seed=42, verbose=0, thread_count=-1\n)\n\n# Define models for LOG-UNCERTAINTY prediction\nridge_log_sigma_model = Ridge(random_state=42)\n\n# Pre-allocate OOF arrays for mean and sigma\nmu_oof_preds = np.zeros(y_train.shape)\nsigma_oof_preds = np.zeros(y_train.shape)\n\n# Cross-validation setup\nN_SPLITS = 5\nkf = KFold(n_splits=N_SPLITS, shuffle=True, random_state=42)\n\nprint(f\"\\nTraining {N_SPLITS}-Fold Cross-Validation Models (Mean & Log-Sigma)...\")\n\ncategorical_features_indices = np.where(X_train.dtypes == 'object')[0].tolist()\nfor col_idx in categorical_features_indices:\n    col_name = X_train.columns[col_idx]\n    X_train[col_name] = X_train[col_name].astype(str) # Ensure string type for CatBoost\n\nfor fold, (train_idx, val_idx) in enumerate(kf.split(X_train)):\n    print(f\"Fold {fold+1}/{N_SPLITS}\")\n    X_train_fold, X_val_fold = X_train.iloc[train_idx], X_train.iloc[val_idx]\n    y_train_fold, y_val_fold = y_train.iloc[train_idx], y_train.iloc[val_idx]\n\n    for i, col in tqdm(enumerate(target_cols), total=len(target_cols), desc=f\"Training for targets in Fold {fold+1}\"):\n        # --- Train Mean Models (mu_pred) ---\n        ridge_mean_model.fit(X_train_fold, y_train_fold[col])\n        \n        cat_mean_model.fit(X_train_fold, y_train_fold[col], cat_features=categorical_features_indices)\n        \n        # Ensemble mean predictions for this fold\n        fold_mu_pred = (ridge_mean_model.predict(X_val_fold) + cat_mean_model.predict(X_val_fold)) / 2\n        fold_mu_pred[fold_mu_pred < 0] = 0 # Ensure non-negative flux\n\n        mu_oof_preds[val_idx, i] = fold_mu_pred\n        \n        # --- Train Log-Sigma Models (sigma_pred) ---\n        cat_mean_preds_on_train = cat_mean_model.predict(X_train_fold)\n        residuals_sq = (y_train_fold[col].values - cat_mean_preds_on_train)**2\n        \n        target_log_var = np.log(residuals_sq + 1e-6) # Add a small epsilon to avoid log(0)\n\n        ridge_log_sigma_model.fit(X_train_fold, target_log_var)\n        \n        log_var_val_pred = ridge_log_sigma_model.predict(X_val_fold)\n        sigma_val_pred = np.sqrt(np.exp(log_var_val_pred))\n        \n        min_sigma_clamp = np.mean(y_train_fold[col].values) * (min(IDEAL_UNC_FGS1_PPM, IDEAL_UNC_AIRS_PPM) / 1e6)\n        sigma_val_pred = np.maximum(sigma_val_pred, min_sigma_clamp)\n\n        sigma_oof_preds[val_idx, i] = sigma_val_pred\n\n    gc.collect()\n\nprint(\"\\nModel Training & Fitting Complete.\")\n\n# --- 4. Computing GLL Score ---\nprint(\"\\n--- 4. Computing GLL Score ---\")\n\n# Ensure predictions are in the correct format for GLL calculation\ny_true_oof = y_train.values\nmu_oof = mu_oof_preds\nsigma_oof = sigma_oof_preds\n\n# Calculate User's GLL (L_user)\n# The calculate_gll function itself should now be using the channel index (i) for instrument type.\nL_user = calculate_gll(y_true_oof, mu_oof, sigma_oof, wavelengths_df)\nprint(f\"User's OOF GLL (L_user): {L_user:.6f}\")\n\n# --- Calculate L_ideal ---\ny_ideal_mu = y_true_oof.copy() # Mu perfectly matches y_true\n\nsigma_ideal = np.zeros_like(y_true_oof)\nfor i in range(NUM_WAVELENGTHS):\n    # This is the line that needed modification:\n    # Instead of looking up 'instrument' in wavelengths_df,\n    # we determine it based on the channel index 'i' (as per problem description).\n    if i == FGS1_CHANNEL_IDX: # Check if it's the FGS1 channel (flux_0)\n        sigma_ideal[:, i] = y_ideal_mu[:, i] * (IDEAL_UNC_FGS1_PPM / 1e6)\n    elif i in AIRS_CH0_CHANNEL_INDICES: # Check if it's one of the AIRS-CH0 channels (flux_1 to flux_355)\n        sigma_ideal[:, i] = y_ideal_mu[:, i] * (IDEAL_UNC_AIRS_PPM / 1e6)\n    # No 'else' needed here, as all columns should fall into one of these categories\n    # if NUM_WAVELENGTHS is correctly defined as 356.\n\nL_ideal = calculate_gll(y_true_oof, y_ideal_mu, sigma_ideal, wavelengths_df)\nprint(f\"Ideal GLL (L_ideal): {L_ideal:.6f}\")\n\n\n# --- Calculate L_ref ---\nmu_ref_base = y_train.mean(axis=0).values\nmu_ref = np.tile(mu_ref_base, (y_true_oof.shape[0], 1))\n\nsigma_ref_base = y_train.std(axis=0).values\nsigma_ref = np.tile(sigma_ref_base, (y_true_oof.shape[0], 1))\n\nL_ref = calculate_gll(y_true_oof, mu_ref, sigma_ref, wavelengths_df)\nprint(f\"Reference GLL (L_ref): {L_ref:.6f}\")\n\n\n# --- Calculate Final Competition Score ---\ncompetition_score = calculate_competition_score(L_user, L_ideal, L_ref)\nprint(f\"\\nFinal OOF Competition Score: {competition_score:.6f}\")\n\nprint(\"\\nGLL Score Computation Complete.\")\n\n# --- 5. Saving Submission File ---\nprint(\"\\n--- 5. Saving Submission File ---\")\n\n# --- Load Test Data for Prediction ---\nsample_submission_df = pd.read_csv(SAMPLE_SUBMISSION_PATH)\ntest_planet_ids = sample_submission_df['planet_id'].tolist()\n\nprint(\"\\n--- Generating FGS1 Features for Test Set ---\")\nstart_time = time.time()\nf_test_features = f_read_and_preprocess(TEST_DIR, test_planet_ids, band_name=\"FGS1\")\nend_time = time.time()\nprint(f\"FGS1 Test Feature Generation took: {end_time - start_time:.2f} seconds\")\n\n\nprint(\"\\n--- Generating AIRS-CH0 Features for Test Set ---\")\nstart_time = time.time()\na_test_features = a_read_and_preprocess(TEST_DIR, test_planet_ids, band_name=\"AIRS-CH0\")\nend_time = time.time()\nprint(f\"AIRS-CH0 Test Feature Generation took: {end_time - start_time:.2f} seconds\")\n\n# Merge features for test set\nX_test = f_test_features.merge(a_test_features, left_index=True, right_index=True, how='left')\n\n# --- Add adc_info features to X_test (same logic as X_train) ---\n# Assuming train_adc_info (loaded at the top of the script) contains the global ADC settings.\nif not train_adc_info.empty:\n    adc_features_values = train_adc_info.iloc[0].to_dict()\n    for col_name, value in adc_features_values.items():\n        X_test[col_name] = value\n    print(\"ADC info features added to X_test.\")\nelse:\n    print(\"Warning: train_adc_info is empty. ADC features will not be added to X_test.\")\n\n\n# --- Ensure Test Set Columns Match Training Set Columns ---\n# This re-creates a dummy X_train to get column names if X_train is not in memory.\ntry:\n    X_train_cols = X_train.columns\nexcept NameError:\n    print(\"X_train not found in memory. Re-generating dummy X_train to get column names for consistency.\")\n    \n    train_labels_dummy = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'), index_col='planet_id')\n    train_star_info_dummy = pd.read_csv(os.path.join(DATA_DIR, 'train_star_info.csv'))\n    train_adc_info_dummy = pd.read_csv(os.path.join(DATA_DIR, 'adc_info.csv')) # Load it again for dummy\n\n    f_train_features_dummy = f_read_and_preprocess(TRAIN_DIR, train_labels_dummy.index.tolist(), band_name=\"FGS1\")\n    a_train_features_dummy = a_read_and_preprocess(TRAIN_DIR, train_labels_dummy.index.tolist(), band_name=\"AIRS-CH0\")\n\n    X_train_dummy = f_train_features_dummy.merge(a_train_features_dummy, left_index=True, right_index=True, how='left')\n    X_train_dummy = X_train_dummy.merge(train_star_info_dummy.set_index('planet_id'), left_index=True, right_index=True, how='left')\n    \n    # CORRECTED: Add adc_info features to X_train_dummy using broadcasting\n    if not train_adc_info_dummy.empty:\n        adc_features_values_dummy = train_adc_info_dummy.iloc[0].to_dict()\n        for col_name, value in adc_features_values_dummy.items():\n            X_train_dummy[col_name] = value\n    \n    X_train_dummy = X_train_dummy.fillna(0)\n    X_train_cols = X_train_dummy.columns\n    del train_labels_dummy, train_star_info_dummy, train_adc_info_dummy, f_train_features_dummy, a_train_features_dummy, X_train_dummy\n    gc.collect()\n\nX_test = X_test.reindex(columns=X_train_cols, fill_value=0)\ncategorical_features_indices_test = np.where(X_test.dtypes == 'object')[0].tolist()\nfor col_idx in categorical_features_indices_test:\n    col_name = X_test.columns[col_idx]\n    X_test[col_name] = X_test[col_name].astype(str)\n\nprint(f\"X_test shape after feature engineering and alignment: {X_test.shape}\")\nprint(\"X_test head:\\n\", X_test.head())\n\n# --- Train Final Models on Full Training Data ---\nprint(\"\\nTraining final models on full training data (mean and log-sigma)...\")\nfinal_mean_models = {}\nfinal_log_sigma_models = {}\n\n# Re-define categorical_features_indices for the full X_train if it's not global\n# This ensures consistency for CatBoost when training on the full dataset.\n# It was defined in the cross-validation loop, but let's ensure it's accessible here.\ncategorical_features_indices_full_train = np.where(X_train.dtypes == 'object')[0].tolist()\n\n\nfor i, col in tqdm(enumerate(target_cols), total=len(target_cols), desc=\"Final Model Training\"):\n    # Train Mean Models (Ensemble of Ridge and CatBoost)\n    ridge_m = Ridge(random_state=42)\n    cat_m = CatBoostRegressor(\n        iterations=500, learning_rate=0.05, depth=6,\n        loss_function='RMSE', eval_metric='RMSE',\n        random_seed=42, verbose=0, thread_count=-1\n    )\n    ridge_m.fit(X_train, y_train[col])\n    cat_m.fit(X_train, y_train[col], cat_features=categorical_features_indices_full_train) \n    \n    final_mean_models[col] = (ridge_m, cat_m)\n\n    # Train Log-Sigma Model (Ridge)\n    log_sigma_m = Ridge(random_state=42)\n    cat_mean_preds_on_full_train = cat_m.predict(X_train)\n    residuals_sq_full_train = (y_train[col].values - cat_mean_preds_on_full_train)**2\n    target_log_var_full_train = np.log(residuals_sq_full_train + 1e-6)\n\n    log_sigma_m.fit(X_train, target_log_var_full_train)\n    final_log_sigma_models[col] = log_sigma_m\n\n# --- Generate Predictions on Test Set ---\nprint(\"\\nGenerating mean (mu) and uncertainty (sigma) predictions for the test set...\")\nmu_test_predictions = np.zeros((len(test_planet_ids), NUM_WAVELENGTHS))\nsigma_test_predictions = np.zeros((len(test_planet_ids), NUM_WAVELENGTHS))\n\nfor i, col in tqdm(enumerate(target_cols), total=len(target_cols), desc=\"Predicting for test targets\"):\n    ridge_m, cat_m = final_mean_models[col]\n    log_sigma_m = final_log_sigma_models[col]\n\n    # Predict mean\n    ridge_mu_pred = ridge_m.predict(X_test)\n    cat_mu_pred = cat_m.predict(X_test)\n    \n    ensemble_mu_pred = (ridge_mu_pred + cat_mu_pred) / 2\n    ensemble_mu_pred[ensemble_mu_pred < 0] = 0\n\n    mu_test_predictions[:, i] = ensemble_mu_pred\n\n    # Predict log-variance and convert to sigma\n    log_var_pred = log_sigma_m.predict(X_test)\n    sigma_pred = np.sqrt(np.exp(log_var_pred))\n    \n    min_sigma_clamp_test = np.mean(ensemble_mu_pred) * (min(IDEAL_UNC_FGS1_PPM, IDEAL_UNC_AIRS_PPM) / 1e6)\n    sigma_pred = np.maximum(sigma_pred, min_sigma_clamp_test)\n\n    sigma_test_predictions[:, i] = sigma_pred\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-14T15:41:01.08965Z","iopub.execute_input":"2025-07-14T15:41:01.09002Z","iopub.status.idle":"2025-07-14T17:07:01.462268Z","shell.execute_reply.started":"2025-07-14T15:41:01.089991Z","shell.execute_reply":"2025-07-14T17:07:01.46128Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- 5. Saving Submission File ---\nprint(\"\\n--- 5. Saving Submission File ---\")\n\n# --- Load Test Data for Prediction ---\nsample_submission_df = pd.read_csv(SAMPLE_SUBMISSION_PATH)\ntest_planet_ids = sample_submission_df['planet_id'].tolist()\n\nprint(\"\\n--- Generating FGS1 Features for Test Set ---\")\nstart_time = time.time()\nf_test_features = f_read_and_preprocess(TEST_DIR, test_planet_ids, band_name=\"FGS1\")\nend_time = time.time()\nprint(f\"FGS1 Test Feature Generation took: {end_time - start_time:.2f} seconds\")\n\n\nprint(\"\\n--- Generating AIRS-CH0 Features for Test Set ---\")\nstart_time = time.time()\na_test_features = a_read_and_preprocess(TEST_DIR, test_planet_ids, band_name=\"AIRS-CH0\")\nend_time = time.time()\nprint(f\"AIRS-CH0 Test Feature Generation took: {end_time - start_time:.2f} seconds\")\n\n# Merge features for test set\nX_test = f_test_features.merge(a_test_features, left_index=True, right_index=True, how='left')\n\n# --- Add adc_info features to X_test (same logic as X_train) ---\n# Assuming train_adc_info (loaded at the top of the script) contains the global ADC settings.\nif not train_adc_info.empty:\n    adc_features_values = train_adc_info.iloc[0].to_dict()\n    for col_name, value in adc_features_values.items():\n        X_test[col_name] = value\n    print(\"ADC info features added to X_test.\")\nelse:\n    print(\"Warning: train_adc_info is empty. ADC features will not be added to X_test.\")\n\n\n# --- Ensure Test Set Columns Match Training Set Columns ---\n# This re-creates a dummy X_train to get column names if X_train is not in memory.\ntry:\n    X_train_cols = X_train.columns\nexcept NameError:\n    print(\"X_train not found in memory. Re-generating dummy X_train to get column names for consistency.\")\n    \n    train_labels_dummy = pd.read_csv(os.path.join(DATA_DIR, 'train.csv'), index_col='planet_id')\n    train_star_info_dummy = pd.read_csv(os.path.join(DATA_DIR, 'train_star_info.csv'))\n    train_adc_info_dummy = pd.read_csv(os.path.join(DATA_DIR, 'adc_info.csv')) # Load it again for dummy\n\n    f_train_features_dummy = f_read_and_preprocess(TRAIN_DIR, train_labels_dummy.index.tolist(), band_name=\"FGS1\")\n    a_train_features_dummy = a_read_and_preprocess(TRAIN_DIR, train_labels_dummy.index.tolist(), band_name=\"AIRS-CH0\")\n\n    X_train_dummy = f_train_features_dummy.merge(a_train_features_dummy, left_index=True, right_index=True, how='left')\n    X_train_dummy = X_train_dummy.merge(train_star_info_dummy.set_index('planet_id'), left_index=True, right_index=True, how='left')\n    \n    # CORRECTED: Add adc_info features to X_train_dummy using broadcasting\n    if not train_adc_info_dummy.empty:\n        adc_features_values_dummy = train_adc_info_dummy.iloc[0].to_dict()\n        for col_name, value in adc_features_values_dummy.items():\n            X_train_dummy[col_name] = value\n    \n    X_train_dummy = X_train_dummy.fillna(0)\n    X_train_cols = X_train_dummy.columns\n    del train_labels_dummy, train_star_info_dummy, train_adc_info_dummy, f_train_features_dummy, a_train_features_dummy, X_train_dummy\n    gc.collect()\n\nX_test = X_test.reindex(columns=X_train_cols, fill_value=0)\ncategorical_features_indices_test = np.where(X_test.dtypes == 'object')[0].tolist()\nfor col_idx in categorical_features_indices_test:\n    col_name = X_test.columns[col_idx]\n    X_test[col_name] = X_test[col_name].astype(str)\n\nprint(f\"X_test shape after feature engineering and alignment: {X_test.shape}\")\nprint(\"X_test head:\\n\", X_test.head())\n\n# --- Train Final Models on Full Training Data ---\nprint(\"\\nTraining final models on full training data (mean and log-sigma)...\")\nfinal_mean_models = {}\nfinal_log_sigma_models = {}\n\n# Re-define categorical_features_indices for the full X_train if it's not global\n# This ensures consistency for CatBoost when training on the full dataset.\n# It was defined in the cross-validation loop, but let's ensure it's accessible here.\ncategorical_features_indices_full_train = np.where(X_train.dtypes == 'object')[0].tolist()\n\n\nfor i, col in tqdm(enumerate(target_cols), total=len(target_cols), desc=\"Final Model Training\"):\n    # Train Mean Models (Ensemble of Ridge and CatBoost)\n    ridge_m = Ridge(random_state=42)\n    cat_m = CatBoostRegressor(\n        iterations=500, learning_rate=0.05, depth=6,\n        loss_function='RMSE', eval_metric='RMSE',\n        random_seed=42, verbose=0, thread_count=-1\n    )\n    ridge_m.fit(X_train, y_train[col])\n    cat_m.fit(X_train, y_train[col], cat_features=categorical_features_indices_full_train) \n    \n    final_mean_models[col] = (ridge_m, cat_m)\n\n    # Train Log-Sigma Model (Ridge)\n    log_sigma_m = Ridge(random_state=42)\n    cat_mean_preds_on_full_train = cat_m.predict(X_train)\n    residuals_sq_full_train = (y_train[col].values - cat_mean_preds_on_full_train)**2\n    target_log_var_full_train = np.log(residuals_sq_full_train + 1e-6)\n\n    log_sigma_m.fit(X_train, target_log_var_full_train)\n    final_log_sigma_models[col] = log_sigma_m\n\n# --- Generate Predictions on Test Set ---\nprint(\"\\nGenerating mean (mu) and uncertainty (sigma) predictions for the test set...\")\nmu_test_predictions = np.zeros((len(test_planet_ids), NUM_WAVELENGTHS))\nsigma_test_predictions = np.zeros((len(test_planet_ids), NUM_WAVELENGTHS))\n\nfor i, col in tqdm(enumerate(target_cols), total=len(target_cols), desc=\"Predicting for test targets\"):\n    ridge_m, cat_m = final_mean_models[col]\n    log_sigma_m = final_log_sigma_models[col]\n\n    # Predict mean\n    ridge_mu_pred = ridge_m.predict(X_test)\n    cat_mu_pred = cat_m.predict(X_test)\n    \n    ensemble_mu_pred = (ridge_mu_pred + cat_mu_pred) / 2\n    ensemble_mu_pred[ensemble_mu_pred < 0] = 0\n\n    mu_test_predictions[:, i] = ensemble_mu_pred\n\n    # Predict log-variance and convert to sigma\n    log_var_pred = log_sigma_m.predict(X_test)\n    sigma_pred = np.sqrt(np.exp(log_var_pred))\n    \n    min_sigma_clamp_test = np.mean(ensemble_mu_pred) * (min(IDEAL_UNC_FGS1_PPM, IDEAL_UNC_AIRS_PPM) / 1e6)\n    sigma_pred = np.maximum(sigma_pred, min_sigma_clamp_test)\n\n    sigma_test_predictions[:, i] = sigma_pred\n\n\n# --- Create Submission DataFrame ---\nmu_cols = [f'flux_{i}' for i in range(NUM_WAVELENGTHS)]\nsigma_cols = [f'uncertainty_{i}' for i in range(NUM_WAVELENGTHS)]\n\nsubmission_mu_df = pd.DataFrame(mu_test_predictions, columns=mu_cols)\nsubmission_sigma_df = pd.DataFrame(sigma_test_predictions, columns=sigma_cols)\n\nsubmission_df = pd.concat([submission_mu_df, submission_sigma_df], axis=1)\nsubmission_df.insert(0, 'planet_id', test_planet_ids)\n\n# --- ADDED DEBUGGING PRINTS ---\nprint(\"\\n--- Final Submission DataFrame Info ---\")\nprint(f\"Submission DataFrame shape: {submission_df.shape}\")\nprint(f\"Submission DataFrame columns: {submission_df.columns.tolist()[:5]} ... {submission_df.columns.tolist()[-5:]}\")\nprint(\"\\nSubmission DataFrame Head (first 5 rows, first 10 columns):\")\nprint(submission_df.iloc[:, :10].head())\nprint(\"\\nSubmission DataFrame Head (first 5 rows, last 10 columns):\")\nprint(submission_df.iloc[:, -10:].head())\nprint(\"\\nChecking for any NaN values in submission_df:\")\nprint(submission_df.isnull().sum().sum()) # Should be 0\n\n# --- Save Submission File ---\nsubmission_file_name = 'submission.csv'\nsubmission_df.to_csv(submission_file_name, index=False)\n\nprint(f\"\\n--- Submission file '{submission_file_name}' created successfully! ---\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-14T17:24:46.888472Z","iopub.execute_input":"2025-07-14T17:24:46.889124Z","iopub.status.idle":"2025-07-14T17:33:00.504988Z","shell.execute_reply.started":"2025-07-14T17:24:46.88909Z","shell.execute_reply":"2025-07-14T17:33:00.503941Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# --- Verify the saved submission file ---\nprint(\"\\n--- Verifying the saved 'submission.csv' file ---\")\ntry:\n    verified_submission_df = pd.read_csv(submission_file_name)\n    print(\"Head of the saved submission.csv:\")\n    print(verified_submission_df.head())\n    print(f\"\\nShape of the saved submission.csv: {verified_submission_df.shape}\")\n    print(f\"Columns of the saved submission.csv (first 5 and last 5): {verified_submission_df.columns.tolist()[:5]} ... {verified_submission_df.columns.tolist()[-5:]}\")\nexcept Exception as e:\n    print(f\"Error reading the saved submission.csv: {e}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-14T17:37:53.4555Z","iopub.execute_input":"2025-07-14T17:37:53.456805Z","iopub.status.idle":"2025-07-14T17:37:53.488082Z","shell.execute_reply.started":"2025-07-14T17:37:53.456756Z","shell.execute_reply":"2025-07-14T17:37:53.486788Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"![Kaguya‑sama cover](https://c4.wallpaperflare.com/wallpaper/835/612/1004/kaguya-sama-love-is-war-kaguya-fan-art-anime-hd-wallpaper-preview.jpg)\n","metadata":{}},{"cell_type":"markdown","source":"**THANKS FOR VISITING MY NOTEBOOK FELLAS!!!!!**","metadata":{}}]}