{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"nvidiaTeslaT4","dataSources":[{"sourceId":101849,"databundleVersionId":12846694,"sourceType":"competition"}],"dockerImageVersionId":31040,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# Import libraries\nimport pandas as pd\nimport polars as pl\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nfrom tqdm import tqdm\nimport pickle\nimport os\nfrom glob import glob\nfrom collections import defaultdict","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:32.556559Z","iopub.execute_input":"2025-07-01T09:48:32.556863Z","iopub.status.idle":"2025-07-01T09:48:34.151511Z","shell.execute_reply.started":"2025-07-01T09:48:32.556844Z","shell.execute_reply":"2025-07-01T09:48:34.150894Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load main CSV files\ntrain_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train.csv')\nwavelengths_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/wavelengths.csv')\ntrain_star_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train_star_info.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.152924Z","iopub.execute_input":"2025-07-01T09:48:34.153557Z","iopub.status.idle":"2025-07-01T09:48:34.363079Z","shell.execute_reply.started":"2025-07-01T09:48:34.153512Z","shell.execute_reply":"2025-07-01T09:48:34.362549Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load main CSV files\ntrain_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train.csv')\nwavelengths_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/wavelengths.csv')\ntrain_star_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train_star_info.csv')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.3637Z","iopub.execute_input":"2025-07-01T09:48:34.363897Z","iopub.status.idle":"2025-07-01T09:48:34.368537Z","shell.execute_reply.started":"2025-07-01T09:48:34.363874Z","shell.execute_reply":"2025-07-01T09:48:34.36791Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Summary stats of flux values (excluding planet_id)\ntarget_cols = [col for col in train_df.columns if col != 'planet_id']\nflux_summary = train_df[target_cols].agg(['min', 'max', 'mean', 'std']).T\nprint(\"Flux value summary (first 5 wavelengths):\")\nprint(flux_summary.head())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.369061Z","iopub.execute_input":"2025-07-01T09:48:34.369241Z","iopub.status.idle":"2025-07-01T09:48:34.567111Z","shell.execute_reply.started":"2025-07-01T09:48:34.369227Z","shell.execute_reply":"2025-07-01T09:48:34.566588Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Count unique stars in star info\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(\"Number of unique stars:\", num_stars)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.568699Z","iopub.execute_input":"2025-07-01T09:48:34.568903Z","iopub.status.idle":"2025-07-01T09:48:34.579763Z","shell.execute_reply.started":"2025-07-01T09:48:34.568888Z","shell.execute_reply":"2025-07-01T09:48:34.578821Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Count AIRS observations per planet by checking parquet files in train folder\nobs_counts = defaultdict(int)\ntrain_planets = os.listdir('/kaggle/input/ariel-data-challenge-2025/train')\n\nfor pid in train_planets:\n    air_obs = glob(f\"train/{pid}/AIRS-CH0_signal_*.parquet\")\n    obs_counts[pid] = len(air_obs)\n\nmulti_obs = {pid: count for pid, count in obs_counts.items() if count > 1}\nprint(\"Planets with multiple AIRS observations:\", len(multi_obs))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.580626Z","iopub.execute_input":"2025-07-01T09:48:34.58091Z","iopub.status.idle":"2025-07-01T09:48:34.615461Z","shell.execute_reply.started":"2025-07-01T09:48:34.580879Z","shell.execute_reply":"2025-07-01T09:48:34.614938Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check for missing calibration files per planet and band\nmissing_calibs = []\nexpected = {\"dark\", \"dead\", \"flat\", \"linear_corr\", \"read\"}\n\nfor pid in train_planets:\n    for band in [\"AIRS-CH0\", \"FGS1\"]:\n        calib_path = f\"train/{pid}/{band}_calibration\"\n        if os.path.exists(calib_path):\n            calib_files = {os.path.splitext(f)[0] for f in os.listdir(calib_path)}\n        else:\n            calib_files = set()\n        missing = expected - calib_files\n        if missing:\n            missing_calibs.append((pid, band, missing))\n\nprint(\"Planets missing calibration files:\", len(missing_calibs))\nif missing_calibs:\n    print(\"Example missing calib:\", missing_calibs[0])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.616091Z","iopub.execute_input":"2025-07-01T09:48:34.616317Z","iopub.status.idle":"2025-07-01T09:48:34.633733Z","shell.execute_reply.started":"2025-07-01T09:48:34.6163Z","shell.execute_reply":"2025-07-01T09:48:34.633065Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Distribution of AIRS observations per planet\nobs_distribution = pd.Series(list(obs_counts.values())).value_counts().sort_index()\nprint(\"Distribution of AIRS observations per planet:\")\nprint(obs_distribution)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.634342Z","iopub.execute_input":"2025-07-01T09:48:34.634604Z","iopub.status.idle":"2025-07-01T09:48:34.645654Z","shell.execute_reply.started":"2025-07-01T09:48:34.634582Z","shell.execute_reply":"2025-07-01T09:48:34.645021Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Check unique planet-star mappings by merging train and star info\nmerged = pd.merge(train_df[['planet_id']], train_star_info, on='planet_id', how='left')\nunique_links = merged.drop_duplicates()\nprint(\"Unique planet-star mappings:\", unique_links.shape[0])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.646421Z","iopub.execute_input":"2025-07-01T09:48:34.646672Z","iopub.status.idle":"2025-07-01T09:48:34.665902Z","shell.execute_reply.started":"2025-07-01T09:48:34.646656Z","shell.execute_reply":"2025-07-01T09:48:34.665089Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Load additional metadata files for reference\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/adc_info.csv')\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train.csv', index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2025/axis_info.parquet')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:34.66694Z","iopub.execute_input":"2025-07-01T09:48:34.667404Z","iopub.status.idle":"2025-07-01T09:48:35.003196Z","shell.execute_reply.started":"2025-07-01T09:48:34.667336Z","shell.execute_reply":"2025-07-01T09:48:35.002218Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Plot a pair of FGS1 images for one planet\nplanet_id = 1010375142\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.parquet')\n\n_, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\nsns.heatmap(f_signal.iloc[0].values.reshape(32, 32), ax=ax1, vmin=0, vmax=52000)\nax1.set_aspect('equal')\nsns.heatmap(f_signal.iloc[1].values.reshape(32, 32), ax=ax2, vmin=0, vmax=52000)\nax2.set_aspect('equal')\nplt.suptitle('A pair of FGS1 images')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:35.004227Z","iopub.execute_input":"2025-07-01T09:48:35.004514Z","iopub.status.idle":"2025-07-01T09:48:37.01413Z","shell.execute_reply.started":"2025-07-01T09:48:35.004485Z","shell.execute_reply":"2025-07-01T09:48:37.013416Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Plot FGS1 time series for planets with strong and weak signals\n_, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, figsize=(12, 4))\n\n# Strong signal planet\nplanet_id = 1048114509\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.parquet')\nmean_signal = f_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow = 800\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\nax1.set_title('FGS1: strong signal planet')\nax1.plot(net_signal)\nax3.plot(smooth_signal, color='c')\nax3.set_xlabel('time step')\n\n# Weak signal planet\nplanet_id = 1240764363\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.parquet')\nmean_signal = f_signal.values.mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\nax2.set_title('FGS1: weak signal planet')\nax2.plot(net_signal)\nax4.plot(smooth_signal, color='c')\nax4.set_xlabel('time step')\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:37.014796Z","iopub.execute_input":"2025-07-01T09:48:37.015028Z","iopub.status.idle":"2025-07-01T09:48:40.093778Z","shell.execute_reply.started":"2025-07-01T09:48:37.015011Z","shell.execute_reply":"2025-07-01T09:48:40.092978Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Plot AIRS-CH0 signal heatmap for one planet\nplanet_id = 1240764363\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/AIRS-CH0_signal_0.parquet')\na_signal = a_signal.values.reshape(11250, 32, 356)  # reshape to (time, spatial, wavelength)\n\nplt.figure(figsize=(10, 3))\nsns.heatmap(a_signal[1])\nplt.ylabel('spatial dimension')\nplt.xlabel('wavelength dimension')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:40.094631Z","iopub.execute_input":"2025-07-01T09:48:40.094885Z","iopub.status.idle":"2025-07-01T09:48:41.861314Z","shell.execute_reply.started":"2025-07-01T09:48:40.094861Z","shell.execute_reply":"2025-07-01T09:48:41.860329Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Plot AIRS-CH0 net signal and smoothened net signal for the same planet\nmean_signal = a_signal.mean(axis=2).mean(axis=1)\nnet_signal = mean_signal[1::2] - mean_signal[0::2]\ncum_signal = net_signal.cumsum()\nwindow = 80\nsmooth_signal = (cum_signal[window:] - cum_signal[:-window]) / window\n\n_, (ax1, ax2) = plt.subplots(2, 1, sharex=True)\nax1.plot(net_signal, label='raw net signal')\nax2.plot(smooth_signal, color='c', label='smoothened net signal')\nax2.set_xlabel('time')\nplt.suptitle('AIRS-CH0 time series', y=0.96)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:41.864189Z","iopub.execute_input":"2025-07-01T09:48:41.864399Z","iopub.status.idle":"2025-07-01T09:48:42.21544Z","shell.execute_reply.started":"2025-07-01T09:48:41.864383Z","shell.execute_reply":"2025-07-01T09:48:42.21478Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Function to read and preprocess FGS1 signals for multiple planets\ndef f_read_and_preprocess(dataset, adc_info, planet_ids):\n    \"\"\"\n    Read FGS1 parquet signals, compute net signals for all planets.\n\n    Returns:\n        numpy array with shape (n_planets, 67500)\n    \"\"\"\n    f_raw = np.full((len(planet_ids), 67500), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(enumerate(planet_ids)):\n        f_signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/{dataset}/{planet_id}/FGS1_signal_0.parquet')\n        mean_signal = f_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / 1024\n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n        f_raw[i] = net_signal\n    return f_raw","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:42.216141Z","iopub.execute_input":"2025-07-01T09:48:42.216677Z","iopub.status.idle":"2025-07-01T09:48:42.221713Z","shell.execute_reply.started":"2025-07-01T09:48:42.21666Z","shell.execute_reply":"2025-07-01T09:48:42.220926Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example: Load and preprocess train FGS1 signals\nf_raw_train = f_read_and_preprocess('train', train_adc_info, train_labels.index)\n\n# Save processed data for later reuse\nwith open('f_raw_train.pickle', 'wb') as f:\n    pickle.dump(f_raw_train, f)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T09:48:42.222648Z","iopub.execute_input":"2025-07-01T09:48:42.222945Z","iopub.status.idle":"2025-07-01T10:10:07.652819Z","shell.execute_reply.started":"2025-07-01T09:48:42.222919Z","shell.execute_reply":"2025-07-01T10:10:07.652229Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def a_read_and_preprocess(dataset, adc_info, planet_ids):\n    \"\"\"\n    Reads AIRS-CH0 parquet files and extracts net signal time series.\n    Returns array of shape (num_planets, 5625) with processed signals.\n    \"\"\"\n    # Preallocate array with correct shape\n    a_raw_train = np.full((len(planet_ids), 5625), np.nan, dtype=np.float32)\n    \n    for i, planet_id in tqdm(enumerate(planet_ids), total=len(planet_ids)):\n        # Read parquet file for current planet\n        a_signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/{dataset}/{planet_id}/AIRS-CH0_signal_0.parquet')\n        \n        # Sum horizontally across spatial dimensions and convert to float\n        mean_signal = a_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / (32*356)\n        \n        # Calculate net signal by differencing consecutive pairs (length 5625)\n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n        \n        # Assign net_signal to preallocated array row\n        a_raw_train[i] = net_signal\n    \n    return a_raw_train\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:10:07.6535Z","iopub.execute_input":"2025-07-01T10:10:07.653773Z","iopub.status.idle":"2025-07-01T10:10:07.659208Z","shell.execute_reply.started":"2025-07-01T10:10:07.653748Z","shell.execute_reply":"2025-07-01T10:10:07.658402Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Example: Load and preprocess train AIRS-CH0 signals\na_raw_train = a_read_and_preprocess('train', train_adc_info, train_labels.index)\n\n# Save processed AIRS data for later reuse\nwith open('a_raw_train.pickle', 'wb') as f:\n    pickle.dump(a_raw_train, f)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:10:07.659992Z","iopub.execute_input":"2025-07-01T10:10:07.66018Z","iopub.status.idle":"2025-07-01T10:24:14.961009Z","shell.execute_reply.started":"2025-07-01T10:10:07.660166Z","shell.execute_reply":"2025-07-01T10:24:14.959986Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Plot example AIRS and FGS1 net signals for a planet\nplanet_idx = 0  # first planet in train_labels\nplt.figure(figsize=(10, 4))\n\nplt.plot(a_raw_train[planet_idx], label='AIRS-CH0 net signal')\nplt.plot(f_raw_train[planet_idx], label='FGS1 net signal')\nplt.legend()\nplt.title('Example planet net signals: AIRS-CH0 vs FGS1')\nplt.xlabel('Time step')\nplt.ylabel('Signal intensity')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:14.961852Z","iopub.execute_input":"2025-07-01T10:24:14.962592Z","iopub.status.idle":"2025-07-01T10:24:15.322683Z","shell.execute_reply.started":"2025-07-01T10:24:14.962568Z","shell.execute_reply":"2025-07-01T10:24:15.322086Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"correlations = []\nfor i in range(len(train_labels)):\n    # Find minimum length of the two signals\n    length = min(len(a_raw_train[i]), len(f_raw_train[i]))\n    \n    # Trim both signals to the same length\n    a_signal_trimmed = a_raw_train[i][:length]\n    f_signal_trimmed = f_raw_train[i][:length]\n    \n    # Calculate correlation coefficient\n    corr = np.corrcoef(a_signal_trimmed, f_signal_trimmed)[0, 1]\n    correlations.append(corr)\n\n# Plot the distribution of correlations\nplt.hist(correlations, bins=30)\nplt.title('Distribution of correlations between AIRS-CH0 and FGS1 signals')\nplt.xlabel('Correlation coefficient')\nplt.ylabel('Count of planets')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:15.323518Z","iopub.execute_input":"2025-07-01T10:24:15.323801Z","iopub.status.idle":"2025-07-01T10:24:15.652486Z","shell.execute_reply.started":"2025-07-01T10:24:15.323778Z","shell.execute_reply":"2025-07-01T10:24:15.651589Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Print some statistics on the correlations\nprint(\"Mean correlation:\", np.mean(correlations))\nprint(\"Median correlation:\", np.median(correlations))\nprint(\"Min correlation:\", np.min(correlations))\nprint(\"Max correlation:\", np.max(correlations))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:15.653114Z","iopub.execute_input":"2025-07-01T10:24:15.653344Z","iopub.status.idle":"2025-07-01T10:24:15.662041Z","shell.execute_reply.started":"2025-07-01T10:24:15.653319Z","shell.execute_reply":"2025-07-01T10:24:15.661179Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Feature engineering function to extract key features from raw signals\ndef feature_engineering(f_raw, a_raw):\n    \"\"\"\n    Extract features capturing relative signal reductions during transit.\n\n    Parameters:\n    - f_raw: ndarray, FGS1 net signals (n_planets, 67500)\n    - a_raw: ndarray, AIRS-CH0 net signals (n_planets, 33750)\n\n    Returns:\n    - DataFrame with two columns: relative reductions for AIRS and FGS1\n    \"\"\"\n    # For FGS1: mean signal during obscured period\n    obscured_f = f_raw[:, 23500:44000].mean(axis=1)\n    # Mean signal during unobscured periods (before and after transit)\n    unobscured_f = (f_raw[:, :20500].mean(axis=1) + f_raw[:, 47000:].mean(axis=1)) / 2\n    # Calculate relative reduction (transit depth proxy)\n    f_relative_reduction = (unobscured_f - obscured_f) / unobscured_f\n\n    # For AIRS-CH0: same logic, adjusted indices for smaller time series\n    obscured_a = a_raw[:, 1958:3666].mean(axis=1)\n    unobscured_a = (a_raw[:, :1708].mean(axis=1) + a_raw[:, 3916:].mean(axis=1)) / 2\n    a_relative_reduction = (unobscured_a - obscured_a) / unobscured_a\n\n    # Build dataframe with features\n    df_features = pd.DataFrame({\n        'a_relative_reduction': a_relative_reduction,\n        'f_relative_reduction': f_relative_reduction\n    })\n\n    return df_features","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:15.662808Z","iopub.execute_input":"2025-07-01T10:24:15.663014Z","iopub.status.idle":"2025-07-01T10:24:15.67422Z","shell.execute_reply.started":"2025-07-01T10:24:15.662969Z","shell.execute_reply":"2025-07-01T10:24:15.673579Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Generate training features from preprocessed signals\ntrain_features = feature_engineering(f_raw_train, a_raw_train)\nprint(train_features.head())\n\n# Initialize Ridge regression model with a very small regularization parameter\nfrom sklearn.linear_model import Ridge\nmodel = Ridge(alpha=1e-12)\n\n# Cross-validate with out-of-fold predictions on training data\nfrom sklearn.model_selection import cross_val_predict\noof_predictions = cross_val_predict(model, train_features, train_labels)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:15.675063Z","iopub.execute_input":"2025-07-01T10:24:15.675457Z","iopub.status.idle":"2025-07-01T10:24:16.09501Z","shell.execute_reply.started":"2025-07-01T10:24:15.67543Z","shell.execute_reply":"2025-07-01T10:24:16.094284Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Evaluate model performance\nfrom sklearn.metrics import r2_score, mean_squared_error\n\nr2 = r2_score(train_labels, oof_predictions)\nrmse = mean_squared_error(train_labels, oof_predictions, squared=False)\nprint(f\"R2 score: {r2:.3f}\")\nprint(f\"Root Mean Squared Error (RMSE): {rmse:.6f}\")\n\n\n# Scatter plot for predictions vs true values for one target (wavelength index 1)\nimport matplotlib.pyplot as plt\n\ncol = 1  # target column index to plot\nplt.scatter(oof_predictions[:, col], train_labels.iloc[:, col], s=15, c='lightgreen')\nplt.gca().set_aspect('equal')\nplt.xlabel('Predicted flux')\nplt.ylabel('True flux')\nplt.title('Predicted vs True Flux at Wavelength Index 1')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:16.095896Z","iopub.execute_input":"2025-07-01T10:24:16.096209Z","iopub.status.idle":"2025-07-01T10:24:16.310853Z","shell.execute_reply.started":"2025-07-01T10:24:16.096176Z","shell.execute_reply":"2025-07-01T10:24:16.310212Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# competition_score.py\n# This function computes a Gaussian Log-Likelihood based competition score.\n\nimport scipy.stats\nimport numpy as np\nimport pandas as pd\n\n# Custom error to show user-friendly message for invalid submissions\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef competition_score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    naive_mean: float,\n    naive_sigma: float,\n    sigma_true: float,\n    row_id_column_name='planet_id'\n) -> float:\n    \"\"\"\n    Calculate the competition score using Gaussian Log Likelihood.\n    \n    Parameters:\n    - solution: true values DataFrame\n    - submission: predicted values DataFrame (means + uncertainties)\n    - naive_mean: mean of training targets (baseline)\n    - naive_sigma: std of training targets (baseline)\n    - sigma_true: assumed true measurement noise (small value)\n    - row_id_column_name: column to drop for matching\n    \n    Returns:\n    - score: float between 0 and 1 (higher is better)\n    \"\"\"\n    # Remove ID column if exists\n    solution = solution.drop(columns=[row_id_column_name], errors='ignore')\n    submission = submission.drop(columns=[row_id_column_name], errors='ignore')\n\n    # Basic checks\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission are not allowed.')\n\n    for col in submission.columns:\n        if not pd.api.types.is_numeric_dtype(submission[col]):\n            raise ParticipantVisibleError(f'Submission column {col} must be numeric.')\n\n    n_wavelengths = len(solution.columns)\n    # Submission must have means + sigmas (2x columns)\n    if len(submission.columns) != 2 * n_wavelengths:\n        raise ParticipantVisibleError('Submission must have twice the number of columns as solution.')\n\n    # Extract predictions and uncertainties\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=1e-15, a_max=None)\n    y_true = solution.values\n\n    # Calculate Gaussian Log Likelihoods\n    GLL_pred = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n    GLL_true = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true))\n    GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean, scale=naive_sigma))\n\n    # Normalize score to [0,1]\n    score = (GLL_pred - GLL_mean) / (GLL_true - GLL_mean)\n    return float(np.clip(score, 0.0, 1.0))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:16.311618Z","iopub.execute_input":"2025-07-01T10:24:16.311915Z","iopub.status.idle":"2025-07-01T10:24:16.319714Z","shell.execute_reply.started":"2025-07-01T10:24:16.311898Z","shell.execute_reply":"2025-07-01T10:24:16.318906Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Postprocessing function to prepare submission DataFrame with uncertainties\n\ndef postprocessing(pred_array, index, sigma_pred, column_names=None):\n    \"\"\"\n    Prepare submission DataFrame combining predicted means and uncertainties.\n\n    Parameters:\n    - pred_array: ndarray of predictions (n_samples, n_wavelengths)\n    - index: DataFrame index (planet IDs)\n    - sigma_pred: float or ndarray of uncertainties (same shape as pred_array)\n    - column_names: list of wavelength column names (optional)\n\n    Returns:\n    - DataFrame concatenating means and uncertainties side by side\n    \"\"\"\n    n_samples, n_waves = pred_array.shape\n\n    if column_names is None:\n        column_names = [f\"wl_{i+1}\" for i in range(n_waves)]\n\n    # If sigma is a scalar, expand it to match pred_array shape\n    if np.isscalar(sigma_pred):\n        sigma_pred = np.full_like(pred_array, sigma_pred)\n\n    assert sigma_pred.shape == pred_array.shape, \"Shape of sigma_pred must match pred_array\"\n\n    df_mean = pd.DataFrame(pred_array.clip(0, None), index=index, columns=column_names)\n    df_sigma = pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i+1}\" for i in range(n_waves)])\n\n    return pd.concat([df_mean, df_sigma], axis=1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:16.321281Z","iopub.execute_input":"2025-07-01T10:24:16.321714Z","iopub.status.idle":"2025-07-01T10:24:16.336849Z","shell.execute_reply.started":"2025-07-01T10:24:16.321686Z","shell.execute_reply":"2025-07-01T10:24:16.336134Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Define sigma_pred (example small uncertainty value; adjust as needed)\nsigma_pred = 0.00001  \n\n# Use the postprocessing function on out-of-fold predictions\noof_df = postprocessing(oof_predictions, train_labels.index, sigma_pred)\n\n# Display the first few rows of the postprocessed predictions\ndisplay(oof_df.head())\n\n# Compute and print the competition score on the training data with OOF predictions\ngll_score = competition_score(\n    solution=train_labels.reset_index(),\n    submission=oof_df.reset_index(),\n    naive_mean=train_labels.values.mean(),\n    naive_sigma=train_labels.values.std(),\n    sigma_true=0.000003  # assumed small noise level in ground truth\n)\n\nprint(f\"Estimated competition score: {gll_score:.3f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-01T10:24:16.33741Z","iopub.execute_input":"2025-07-01T10:24:16.337644Z","iopub.status.idle":"2025-07-01T10:24:16.441015Z","shell.execute_reply.started":"2025-07-01T10:24:16.337619Z","shell.execute_reply":"2025-07-01T10:24:16.440246Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true},"outputs":[],"execution_count":null}]}