{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"sourceType":"competition"},{"sourceId":12696895,"sourceType":"datasetVersion","datasetId":7978785},{"sourceId":12757112,"sourceType":"datasetVersion","datasetId":8064654},{"sourceId":13041871,"sourceType":"datasetVersion","datasetId":8024087},{"sourceId":13162856,"sourceType":"datasetVersion","datasetId":8164246}],"dockerImageVersionId":31089,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ariel Data Challenge 2025 Final Submission\n\nNotebook for my final, and best scoring submission to this years Ariel Data Challenge. This competition focused on revealing the atmospheric composition of 1,000 exoplanets with simulated data from the ESA's ariel mission. This competition used Gaussian Log-Likelihood as the evaluation metric, requiring both accurate predictions and proper uncertainty quantification. \n\n## My Approach\n\nOne of my biggest challenges in this competition was understanding and working with the large amount of complex data. Working within the Kaggle notebook and working with small samples of the training data increased my productivity immensely as I didn't need to download or save any of the data, and Kaggle provides 60Gb of disk space and 30Gb of Ram (better than my laptop). While the output format is straightforward (a single prediction vector per planet), another challenge lied in accurately modeling the complex relationships between spectral channels while quantifying prediction uncertainty. \n\nI started with an EDA, then worked on creating submissions that would not throw an error (I think I misunderstood the sample file at first), I then implemented a Random Forest Regressor to start with. I immediately noticed that the model was overfitting to the data, and I had already figured that it wouldn't alone suffice. I looked around and got some inspiration from the 1st place team in the 2024 challenge (https://www.kaggle.com/code/cnumber/neurips-ariel-data-challenge-2024-final-submission). I noticed that they implemented an ensemble method, including using Gaussian Process Regression. I played around with GPR for a long time, and finally settled on my own ensemble with GPR, XGBoost, and Ridge regression. Tuning the models, the mix of the ensemble, the features all took a long time to work on, and eventually got me to the score of 0.311, which I am incredibly happy with. If I had more time to work on this, or if I could start fresh, I would experiment with different or at least more models in the ensemble, better uncertanty quantification, and more robust testing methods to hopefully see larger increases in scores when submitting. ","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport os\nimport pickle\nimport scipy.stats as stats\nfrom scipy.stats import skew, kurtosis\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport itertools\nimport glob \nfrom astropy.stats import sigma_clip\nfrom tqdm import tqdm\nfrom joblib import Parallel, delayed, dump, load\nfrom scipy.optimize import minimize\n\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.linear_model import Ridge, LinearRegression\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import mean_squared_error, r2_score\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.multioutput import MultiOutputRegressor\nfrom sklearn.utils import resample\nimport xgboost as xgb\nfrom xgboost import XGBRegressor\nfrom sklearn.gaussian_process import GaussianProcessRegressor\nfrom sklearn.gaussian_process.kernels import RBF, WhiteKernel, Matern, ConstantKernel as C","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"ariel_path = \"/kaggle/input/ariel-data-challenge-2025/\"\nall_files = os.listdir(ariel_path)\npath_out = \"/kaggle/tmp/light_data_raw/\"\noutput_dir = \"/kaggle/tmp/light_data_raw/\"\n\n#training files\ntrain_folder = \"/kaggle/input/ariel-data-challenge-2025/train\"\nsample_path = \"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\"\ntrain_df = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/train.csv\")\ntrain_star_info = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/train_star_info.csv\")\n\n#test files\ntest_folder = \"/kaggle/input/ariel-data-challenge-2025/test/\"\ntest_file = os.listdir(test_folder)[0]\ntest_df = pd.read_csv(sample_path)\n\ntest_star_info = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/test_star_info.csv\")\n\n## other files\naxis_info_df = pd.read_parquet(\"/kaggle/input/ariel-data-challenge-2025/axis_info.parquet\")\nwavelengths_df = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/wavelengths.csv\")\nadc_info_df = pd.read_csv(\"/kaggle/input/ariel-data-challenge-2025/adc_info.csv\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:28.400661Z","iopub.execute_input":"2025-09-21T14:43:28.401082Z","iopub.status.idle":"2025-09-21T14:43:28.839374Z","shell.execute_reply.started":"2025-09-21T14:43:28.401048Z","shell.execute_reply":"2025-09-21T14:43:28.838334Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Data Preprocessing","metadata":{}},{"cell_type":"code","source":"if not os.path.exists(path_out):\n    os.makedirs(path_out)\n    print(f\"Directory {path_out} created.\")\nelse:\n    print(f\"Directory {path_out} already exists.\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:28.840215Z","iopub.execute_input":"2025-09-21T14:43:28.840492Z","iopub.status.idle":"2025-09-21T14:43:28.847441Z","shell.execute_reply.started":"2025-09-21T14:43:28.840469Z","shell.execute_reply":"2025-09-21T14:43:28.846471Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"CHUNKS_SIZE=1","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:28.848276Z","iopub.execute_input":"2025-09-21T14:43:28.848574Z","iopub.status.idle":"2025-09-21T14:43:28.870581Z","shell.execute_reply.started":"2025-09-21T14:43:28.848543Z","shell.execute_reply":"2025-09-21T14:43:28.869533Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"#Functions for preprocessing\n# restore the original dynamic range of the data as outlined in the competition\ndef ADC_convert(signal, gain=0.4369, offset=-1000):\n    signal = signal.astype(np.float64)\n    signal /= gain\n    signal += offset\n    return signal\n\ndef mask_hot_dead(signal, dead, dark):\n    hot = sigma_clip(\n        dark, sigma=5, maxiters=5\n    ).mask\n    hot = np.tile(hot, (signal.shape[0], 1, 1))\n    dead = np.tile(dead, (signal.shape[0], 1, 1))\n    signal = np.ma.masked_where(dead, signal)\n    signal = np.ma.masked_where(hot, signal)\n\n    return signal\n\n#linearity correction\ndef apply_linear_corr(linear_corr, clean_signal):\n    linear_corr = np.flip(linear_corr, axis=0)\n\n    for x, y in itertools.product(\n        range(clean_signal.shape[1]), range(clean_signal.shape[2])\n    ):\n        poli = np.poly1d(linear_corr[:, x, y])\n        clean_signal[: x, y] = poli(clean_signal[:, x, y])\n\n    return clean_signal\n\n#dark current subtraction\ndef clean_dark(signal, dead, dark, dt):\n    dark = np.ma.masked_where(dead, dark)\n    dark = np.tile(dark, (signal.shape[0], 1, 1))\n\n    signal -= dark*dt[:, np.newaxis, np.newaxis]\n\n    return signal\n\n#correlated double sampling\ndef get_cds(signal):\n    cds = signal[:,1::2,:,:] - signal[:,::2,:,:]\n\n    return cds\n\n#time binning (observations binned together by frequency)\ndef bin_obs(cds_signal, binning):\n    cds_transposed = cds_signal.transpose(0,1,3,2)\n    cds_binned = np.zeros((cds_transposed.shape[0], \n                           cds_transposed.shape[1]//binning, \n                           cds_transposed.shape[2], \n                           cds_transposed.shape[3]))\n\n    for i in range(cds_transposed.shape[1]//binning):\n        cds_binned[:,i,:,:] = np.sum(cds_transposed[:,i*binning:(i+1)*binning,:,:],\n                                    axis=1)\n\n    return cds_binned\n\n#flat field correction\ndef correct_flat_field(flat, dead, signal):\n    flat = flat.transpose(1,0)\n    dead = dead.transpose(1,0)\n    flat = np.ma.masked_where(dead, flat)\n    flat = np.tile(flat, (signal.shape[0], 1, 1))\n    signal = signal / flat\n\n    return signal\n\n#get index for training data\ndef get_index(files, CHUNKS_SIZE):\n    index = []\n\n    for file in files:\n        file_name = file.split('/')[-1]\n\n        if file_name.split('_')[0] == 'AIRS-CH0' and file_name.split('_')[1] == 'signal' and file_name.split('_')[2] == '0.parquet':\n            file_index = os.path.basename(os.path.dirname(file))\n            index.append(int(file_index))\n\n    index = np.array(index)\n    index = np.sort(index)\n\n    index = np.array_split(index, len(index)//CHUNKS_SIZE)\n\n    return index    ","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:28.871808Z","iopub.execute_input":"2025-09-21T14:43:28.87215Z","iopub.status.idle":"2025-09-21T14:43:28.999999Z","shell.execute_reply.started":"2025-09-21T14:43:28.872116Z","shell.execute_reply":"2025-09-21T14:43:28.999254Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Save the train data so we can retreive it quickly instead of running this everytime.","metadata":{}},{"cell_type":"code","source":"# files = glob.glob(os.path.join(ariel_path + 'train/', '*/*'))\n# index = get_index(files, CHUNKS_SIZE)\n\n# axis_info = pd.read_parquet(os.path.join(ariel_path, 'axis_info.parquet'))\n# DO_MASK = True\n# DO_THE_NL_CORR = False\n# DO_DARK = True\n# DO_FLAT = True\n# TIME_BINNING = True\n\n# cut_inf, cut_sup = 39, 321\n# l = cut_sup - cut_inf\n\n# for n, index_chunk in enumerate(tqdm(index)):\n#     AIRS_CH0_clean = np.zeros((CHUNKS_SIZE, 11250, 32, l))\n#     FGS1_clean = np.zeros((CHUNKS_SIZE, 135000, 32, 32))\n\n#     for i in range(CHUNKS_SIZE):\n#         df = pd.read_parquet(os.path.join(ariel_path, f'train/{index_chunk[i]}/AIRS-CH0_signal_0.parquet'))\n#         signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 356))\n#         signal = ADC_convert(signal,)\n\n#         dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n#         dt_airs[1::2] += 0.1\n#         chopped_signal = signal[:, :, cut_inf:cut_sup]\n        \n#         del signal, df\n\n#         flat = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/AIRS-CH0_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n#         dark = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/AIRS-CH0_calibration_0/dark.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n#         dead_airs = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/AIRS-CH0_calibration_0/dead.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n#         linear_corr = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/AIRS-CH0_calibration_0/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n\n#         if DO_MASK:\n#             chopped_signal = mask_hot_dead(chopped_signal, dead_airs, dark)\n#             AIRS_CH0_clean[i] = chopped_signal\n#         else:\n#             AIRS_CH0_clean[i] = chopped_signal\n            \n#         if DO_THE_NL_CORR: \n#             linear_corr_signal = apply_linear_corr(linear_corr,AIRS_CH0_clean[i])\n#             AIRS_CH0_clean[i,:, :, :] = linear_corr_signal\n#         del linear_corr\n\n#         if DO_DARK: \n#             cleaned_signal = clean_dark(AIRS_CH0_clean[i], dead_airs, dark, dt_airs)\n#             AIRS_CH0_clean[i] = cleaned_signal\n#         else: \n#             pass\n#         del dark\n\n#         df = pd.read_parquet(os.path.join(ariel_path, f'train/{index_chunk[i]}/FGS1_signal_0.parquet'))\n#         fgs_signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 32))\n\n#         fgs_signal = ADC_convert(fgs_signal, )\n#         dt_fgs1 = np.ones(len(fgs_signal))*0.1\n#         dt_fgs1[1::2] += 0.1\n#         chopped_FGS1 = fgs_signal\n\n#         del fgs_signal, df\n\n# # cleaning fgs1\n#         flat = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/FGS1_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n#         dark = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/FGS1_calibration_0/dark.parquet')).values.astype(np.float64).reshape((32, 32))\n#         dead_fgs1 = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/FGS1_calibration_0/dead.parquet')).values.astype(np.float64).reshape((32, 32))\n#         linear_corr = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/FGS1_calibration_0/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 32))\n        \n#         if DO_MASK:\n#             chopped_FGS1 = mask_hot_dead(chopped_FGS1, dead_fgs1, dark)\n#             FGS1_clean[i] = chopped_FGS1\n#         else:\n#             FGS1_clean[i] = chopped_FGS1\n\n#         if DO_THE_NL_CORR: \n#             linear_corr_signal = apply_linear_corr(linear_corr,FGS1_clean[i])\n#             FGS1_clean[i,:, :, :] = linear_corr_signal\n#         del linear_corr\n        \n#         if DO_DARK: \n#             cleaned_signal = clean_dark(FGS1_clean[i], dead_fgs1, dark,dt_fgs1)\n#             FGS1_clean[i] = cleaned_signal\n#         else: \n#             pass\n#         del dark\n        \n#     AIRS_cds = get_cds(AIRS_CH0_clean)\n#     FGS1_cds = get_cds(FGS1_clean)\n    \n#     del AIRS_CH0_clean, FGS1_clean\n    \n#     if TIME_BINNING:\n#         AIRS_cds_binned = bin_obs(AIRS_cds,binning=30)\n#         FGS1_cds_binned = bin_obs(FGS1_cds,binning=30*12)\n#     else:\n#         AIRS_cds = AIRS_cds.transpose(0,1,3,2) ## this is important to make it consistent for flat fielding, but you can always change it\n#         AIRS_cds_binned = AIRS_cds\n#         FGS1_cds = FGS1_cds.transpose(0,1,3,2)\n#         FGS1_cds_binned = FGS1_cds\n    \n#     del AIRS_cds, FGS1_cds\n    \n#     for i in range (CHUNKS_SIZE):\n#         flat_airs = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/AIRS-CH0_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n#         flat_fgs = pd.read_parquet(os.path.join(ariel_path,f'train/{index_chunk[i]}/FGS1_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n#         if DO_FLAT:\n#             corrected_AIRS_cds_binned = correct_flat_field(flat_airs,dead_airs, AIRS_cds_binned[i])\n#             AIRS_cds_binned[i] = corrected_AIRS_cds_binned\n#             corrected_FGS1_cds_binned = correct_flat_field(flat_fgs,dead_fgs1, FGS1_cds_binned[i])\n#             FGS1_cds_binned[i] = corrected_FGS1_cds_binned\n#         else:\n#             pass\n\n#     np.save(os.path.join(path_out, 'AIRS_clean_train_{}.npy'.format(n)), AIRS_cds_binned)\n#     np.save(os.path.join(path_out, 'FGS1_train_{}.npy'.format(n)), FGS1_cds_binned)\n#     del AIRS_cds_binned\n#     del FGS1_cds_binned","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:29.000991Z","iopub.execute_input":"2025-09-21T14:43:29.0013Z","iopub.status.idle":"2025-09-21T14:43:29.021872Z","shell.execute_reply.started":"2025-09-21T14:43:29.001269Z","shell.execute_reply":"2025-09-21T14:43:29.020821Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Copy the same code but we need to run it each time we submit for the hidden test set.","metadata":{}},{"cell_type":"code","source":"#test files\nfiles = glob.glob(os.path.join(ariel_path + 'test/', '*/*'))\nindex = get_index(files, CHUNKS_SIZE)\n\naxis_info = pd.read_parquet(os.path.join(ariel_path, 'axis_info.parquet'))\nDO_MASK = True\nDO_THE_NL_CORR = False\nDO_DARK = True\nDO_FLAT = True\nTIME_BINNING = True\n\ncut_inf, cut_sup = 39, 321\nl = cut_sup - cut_inf\n\nfor n, index_chunk in enumerate(tqdm(index)):\n    AIRS_CH0_clean = np.zeros((CHUNKS_SIZE, 11250, 32, l))\n    FGS1_clean = np.zeros((CHUNKS_SIZE, 135000, 32, 32))\n\n    for i in range(CHUNKS_SIZE):\n        df = pd.read_parquet(os.path.join(ariel_path, f'test/{index_chunk[i]}/AIRS-CH0_signal_0.parquet'))\n        signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 356))\n        signal = ADC_convert(signal,)\n\n        dt_airs = axis_info['AIRS-CH0-integration_time'].dropna().values\n        dt_airs[1::2] += 0.1\n        chopped_signal = signal[:, :, cut_inf:cut_sup]\n        \n        del signal, df\n\n        flat = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/AIRS-CH0_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        dark = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/AIRS-CH0_calibration_0/dark.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        dead_airs = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/AIRS-CH0_calibration_0/dead.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        linear_corr = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/AIRS-CH0_calibration_0/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 356))[:, :, cut_inf:cut_sup]\n\n        if DO_MASK:\n            chopped_signal = mask_hot_dead(chopped_signal, dead_airs, dark)\n            AIRS_CH0_clean[i] = chopped_signal\n        else:\n            AIRS_CH0_clean[i] = chopped_signal\n            \n        if DO_THE_NL_CORR: \n            linear_corr_signal = apply_linear_corr(linear_corr,AIRS_CH0_clean[i])\n            AIRS_CH0_clean[i,:, :, :] = linear_corr_signal\n        del linear_corr\n\n        if DO_DARK: \n            cleaned_signal = clean_dark(AIRS_CH0_clean[i], dead_airs, dark, dt_airs)\n            AIRS_CH0_clean[i] = cleaned_signal\n        else: \n            pass\n        del dark\n\n        df = pd.read_parquet(os.path.join(ariel_path, f'test/{index_chunk[i]}/FGS1_signal_0.parquet'))\n        fgs_signal = df.values.astype(np.float64).reshape((df.shape[0], 32, 32))\n\n        fgs_signal = ADC_convert(fgs_signal, )\n        dt_fgs1 = np.ones(len(fgs_signal))*0.1\n        dt_fgs1[1::2] += 0.1\n        chopped_FGS1 = fgs_signal\n\n        del fgs_signal, df\n\n# cleaning fgs1\n        flat = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/FGS1_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n        dark = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/FGS1_calibration_0/dark.parquet')).values.astype(np.float64).reshape((32, 32))\n        dead_fgs1 = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/FGS1_calibration_0/dead.parquet')).values.astype(np.float64).reshape((32, 32))\n        linear_corr = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/FGS1_calibration_0/linear_corr.parquet')).values.astype(np.float64).reshape((6, 32, 32))\n        \n        if DO_MASK:\n            chopped_FGS1 = mask_hot_dead(chopped_FGS1, dead_fgs1, dark)\n            FGS1_clean[i] = chopped_FGS1\n        else:\n            FGS1_clean[i] = chopped_FGS1\n\n        if DO_THE_NL_CORR: \n            linear_corr_signal = apply_linear_corr(linear_corr,FGS1_clean[i])\n            FGS1_clean[i,:, :, :] = linear_corr_signal\n        del linear_corr\n        \n        if DO_DARK: \n            cleaned_signal = clean_dark(FGS1_clean[i], dead_fgs1, dark,dt_fgs1)\n            FGS1_clean[i] = cleaned_signal\n        else: \n            pass\n        del dark\n        \n    AIRS_cds = get_cds(AIRS_CH0_clean)\n    FGS1_cds = get_cds(FGS1_clean)\n    \n    del AIRS_CH0_clean, FGS1_clean\n    \n    if TIME_BINNING:\n        AIRS_cds_binned = bin_obs(AIRS_cds,binning=30)\n        FGS1_cds_binned = bin_obs(FGS1_cds,binning=30*12)\n    else:\n        AIRS_cds = AIRS_cds.transpose(0,1,3,2) ## this is important to make it consistent for flat fielding, but you can always change it\n        AIRS_cds_binned = AIRS_cds\n        FGS1_cds = FGS1_cds.transpose(0,1,3,2)\n        FGS1_cds_binned = FGS1_cds\n    \n    del AIRS_cds, FGS1_cds\n    \n    for i in range (CHUNKS_SIZE):\n        flat_airs = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/AIRS-CH0_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 356))[:, cut_inf:cut_sup]\n        flat_fgs = pd.read_parquet(os.path.join(ariel_path,f'test/{index_chunk[i]}/FGS1_calibration_0/flat.parquet')).values.astype(np.float64).reshape((32, 32))\n        if DO_FLAT:\n            corrected_AIRS_cds_binned = correct_flat_field(flat_airs,dead_airs, AIRS_cds_binned[i])\n            AIRS_cds_binned[i] = corrected_AIRS_cds_binned\n            corrected_FGS1_cds_binned = correct_flat_field(flat_fgs,dead_fgs1, FGS1_cds_binned[i])\n            FGS1_cds_binned[i] = corrected_FGS1_cds_binned\n        else:\n            pass\n\n    #test files\n    np.save(os.path.join(\"/kaggle/working/\", 'AIRS_clean_test_{}.npy'.format(n)), AIRS_cds_binned)\n    np.save(os.path.join(\"/kaggle/working/\", 'FGS1_test_{}.npy'.format(n)), FGS1_cds_binned)\n    del AIRS_cds_binned\n    del FGS1_cds_binned","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:29.024847Z","iopub.execute_input":"2025-09-21T14:43:29.025123Z","iopub.status.idle":"2025-09-21T14:43:45.412761Z","shell.execute_reply.started":"2025-09-21T14:43:29.025101Z","shell.execute_reply":"2025-09-21T14:43:45.411951Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def load_data(file, chunk_size, nb_files):\n    data0 = np.load(file + '_0.npy')\n    data_all = np.zeros((nb_files*chunk_size,\n                        data0.shape[1], data0.shape[2], \n                        data0.shape[3]))\n    data_all[:chunk_size] = data0\n    \n    for i in range (1, nb_files):\n        data_all[i*chunk_size:(i+1)*chunk_size] = np.load(file + '_{}.npy'.format(i))\n    \n    return data_all\n\n#train files\n# data_train = load_data(path_out + 'AIRS_clean_train', CHUNKS_SIZE, len(index))\n# data_train_FGS = load_data(path_out + 'FGS1_train', CHUNKS_SIZE, len(index))\n\n#test files\ndata_test = load_data(\"/kaggle/working/\" + 'AIRS_clean_test', CHUNKS_SIZE, len(index))\ndata_test_FGS = load_data(\"/kaggle/working/\" + 'FGS1_test', CHUNKS_SIZE, len(index))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:45.413515Z","iopub.execute_input":"2025-09-21T14:43:45.413789Z","iopub.status.idle":"2025-09-21T14:43:45.430465Z","shell.execute_reply.started":"2025-09-21T14:43:45.413768Z","shell.execute_reply":"2025-09-21T14:43:45.429468Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# np.save('./' + 'data_train.npy', data_train)\n# np.save('./' + 'data_train_FGS.npy', data_train_FGS)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:45.431376Z","iopub.execute_input":"2025-09-21T14:43:45.431691Z","iopub.status.idle":"2025-09-21T14:43:45.452544Z","shell.execute_reply.started":"2025-09-21T14:43:45.431669Z","shell.execute_reply":"2025-09-21T14:43:45.451704Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# data_train = np.load(\"/kaggle/input/adc-preprocessed-data/data_train (1).npy\")\n# data_train_FGS = np.load(\"/kaggle/input/adc-preprocessed-data/data_train_FGS.npy\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T15:27:53.697953Z","iopub.execute_input":"2025-09-21T15:27:53.6983Z","iopub.status.idle":"2025-09-21T15:31:19.47612Z","shell.execute_reply.started":"2025-09-21T15:27:53.698273Z","shell.execute_reply":"2025-09-21T15:31:19.474955Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"planet_test_ids = [int(planet_id[0]) for planet_id in index]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T15:32:08.070714Z","iopub.execute_input":"2025-09-21T15:32:08.071136Z","iopub.status.idle":"2025-09-21T15:32:29.963298Z","shell.execute_reply.started":"2025-09-21T15:32:08.071095Z","shell.execute_reply":"2025-09-21T15:32:29.962487Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Feature extraction and train test split","metadata":{}},{"cell_type":"code","source":"def extract_physical_features(signal: np.ndarray, dt: float = 1.0):\n    \"\"\"\n    extract physically motivated features from a detector time series cube (T,H,W).\n    \"\"\"\n    T, H, W = signal.shape\n    features = {}\n\n    # total flux\n    F_tot = signal.sum(axis=(1, 2))\n    features[\"flux_mean\"] = F_tot.mean()\n    features[\"flux_std\"] = F_tot.std()\n    features[\"flux_skew\"] = skew(F_tot)\n    features[\"flux_kurt\"] = kurtosis(F_tot)\n\n    # centroid\n    x_coords = np.arange(W)\n    y_coords = np.arange(H)\n    X, Y = np.meshgrid(x_coords, y_coords)\n\n    x_c = (signal * X[None, :, :]).sum(axis=(1, 2)) / F_tot\n    y_c = (signal * Y[None, :, :]).sum(axis=(1, 2)) / F_tot\n    features[\"centroid_x_std\"] = x_c.std()\n    features[\"centroid_y_std\"] = y_c.std()\n\n    # centroid scatter skew/kurt\n    features[\"centroid_x_skew\"] = skew(x_c)\n    features[\"centroid_y_skew\"] = skew(y_c)\n\n    # background mean (outer rows/cols)\n    outer_rows = np.r_[0, 1, H-2, H-1]\n    outer_cols = np.r_[0, 1, W-2, W-1]\n    background = np.concatenate([\n        signal[:, outer_rows, :].reshape(T, -1),\n        signal[:, :, outer_cols].reshape(T, -1)\n    ], axis=1)\n    features[\"background_mean\"] = background.mean()\n    features[\"background_std\"] = background.std()\n\n    # spatial RMS\n    spatial_rms = np.sqrt(np.mean(\n        (signal - signal.mean(axis=(1, 2), keepdims=True)) ** 2,\n        axis=(1, 2)\n    ))\n    features[\"spatial_rms_mean\"] = spatial_rms.mean()\n    features[\"spatial_rms_std\"] = spatial_rms.std()\n\n    # flux derivative\n    flux_deriv = np.diff(F_tot) / dt\n    features[\"flux_deriv_mean\"] = flux_deriv.mean()\n    features[\"flux_deriv_std\"] = flux_deriv.std()\n\n    return features\n\n\ndef extract_frequency_features(signal: np.ndarray):\n    \"\"\"\n    extract frequency domain features using FFT of total flux.\n    \"\"\"\n    features = {}\n    F_tot = signal.sum(axis=(1, 2))\n    fft_vals = np.fft.fft(F_tot)\n    fft_mag = np.abs(fft_vals[:len(fft_vals)//2])  # keep positive freqs\n    freqs = np.fft.fftfreq(len(F_tot))[:len(fft_vals)//2]\n\n    # dominant frequency\n    dom_idx = np.argmax(fft_mag)\n    features[\"dom_freq\"] = freqs[dom_idx]\n    features[\"dom_amp\"] = fft_mag[dom_idx]\n\n    # spectral centroid\n    power = fft_mag ** 2\n    if power.sum() > 0:\n        features[\"spec_centroid\"] = (freqs * power).sum() / power.sum()\n        features[\"spec_spread\"] = np.sqrt(((freqs - features[\"spec_centroid\"]) ** 2 * power).sum() / power.sum())\n    else:\n        features[\"spec_centroid\"] = 0\n        features[\"spec_spread\"] = 0\n\n    return features\n\n\n\n#all features utilizing above functions \n\ndef extract_all_features(data_airs, data_fgs, planet_ids, star_info, spectra_df):\n    all_features = []\n    all_labels = []\n\n    for i, planet_id in enumerate(planet_ids):\n        airs_signal = data_airs[i]\n        fgs_signal = data_fgs[i]\n\n        # physics-inspired + FFT features\n        airs_feats = {**extract_physical_features(airs_signal),\n                      **extract_frequency_features(airs_signal)}\n        fgs_feats = {**extract_physical_features(fgs_signal),\n                     **extract_frequency_features(fgs_signal)}\n\n        # metadata\n        meta = star_info[star_info[\"planet_id\"] == int(planet_id)].iloc[0]\n        label = spectra_df[spectra_df[\"planet_id\"] == int(planet_id)].drop(columns=[\"planet_id\"]).iloc[0].to_numpy()\n\n        combined = {\n            \"planet_id\": planet_id,\n            \"Ts\": meta[\"Ts\"],\n            \"Mp\": meta[\"Mp\"],\n            \"P\": meta[\"P\"],\n            **{f\"AIRS_{k}\": v for k, v in airs_feats.items()},\n            **{f\"FGS_{k}\": v for k, v in fgs_feats.items()}\n        }\n\n        all_features.append(combined)\n        all_labels.append(label)\n\n    return pd.DataFrame(all_features), np.array(all_labels)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:45.536957Z","iopub.execute_input":"2025-09-21T14:43:45.53718Z","iopub.status.idle":"2025-09-21T14:43:45.554781Z","shell.execute_reply.started":"2025-09-21T14:43:45.537162Z","shell.execute_reply":"2025-09-21T14:43:45.553749Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Extract features, train test split, and scale (scalar saved for test data).","metadata":{}},{"cell_type":"code","source":"# X_df, y = extract_all_features(data_train, data_train_FGS, planet_ids, train_star_info, train_df)\n# scaler = StandardScaler()\n# X_numeric = X_df.drop(columns=['planet_id'])\n# scaler = StandardScaler()\n# X_scaled = scaler.fit_transform(X_numeric)\n# X_train, X_val, y_train, y_val = train_test_split(X_scaled, y, test_size=0.3, random_state=1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T15:32:33.962079Z","iopub.execute_input":"2025-09-21T15:32:33.962357Z","iopub.status.idle":"2025-09-21T15:33:12.324208Z","shell.execute_reply.started":"2025-09-21T15:32:33.962336Z","shell.execute_reply":"2025-09-21T15:33:12.323314Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Build the models ","metadata":{}},{"cell_type":"code","source":"# def build_custom_kernel(y):\n#     kernel1 = C(y.max() - y.min(), (1e-9, 1e3)) * RBF(10, (1, 1e5))\n#     kernel2 = C(y.max() - y.min(), (1e-9, 1e3)) * Matern(length_scale=10, length_scale_bounds=(1, 1e5), nu=1.5)\n#     return kernel1 + kernel2\n\n# # Fit GPR\n# def train_gpr(X_train, y_train, noise_level=1e-2):\n#     kernel = build_custom_kernel(y_train)\n#     gpr = GaussianProcessRegressor(kernel=kernel, alpha=noise_level**2, n_restarts_optimizer=10, normalize_y=True)\n#     gpr.fit(X_train, y_train)\n#     return gpr","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:45.576258Z","iopub.execute_input":"2025-09-21T14:43:45.576583Z","iopub.status.idle":"2025-09-21T14:43:45.595236Z","shell.execute_reply.started":"2025-09-21T14:43:45.576554Z","shell.execute_reply":"2025-09-21T14:43:45.594353Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ridge = Ridge(alpha=1e-12)\n# xgb = XGBRegressor(n_estimators=200, max_depth=5, learning_rate=0.1, random_state=42)\n\n# ridge.fit(X_train, y_train)\n# xgb.fit(X_train, y_train)\n# gpr = train_gpr(X_train, y_train)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:45.596257Z","iopub.execute_input":"2025-09-21T14:43:45.596536Z","iopub.status.idle":"2025-09-21T14:43:45.614797Z","shell.execute_reply.started":"2025-09-21T14:43:45.596513Z","shell.execute_reply":"2025-09-21T14:43:45.613517Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Load the models and the saved scaler, extract features for test data, and normalize the data","metadata":{}},{"cell_type":"code","source":"ridge = load(\"/kaggle/input/ensemble-models/ridge.joblib\")\nxgb = load(\"/kaggle/input/ensemble-models/xgb.joblib\")\ngpr = load(\"/kaggle/input/ensemble-models/gpr.joblib\")\nscaler = load(\"/kaggle/input/ensemble-models/scaler.joblib\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:45.615714Z","iopub.execute_input":"2025-09-21T14:43:45.616488Z","iopub.status.idle":"2025-09-21T14:43:47.304897Z","shell.execute_reply.started":"2025-09-21T14:43:45.616454Z","shell.execute_reply":"2025-09-21T14:43:47.303867Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ## test data for competition submission\nX_test_df, y_test_df = extract_all_features(data_test, data_test_FGS, planet_test_ids, test_star_info, test_df)\nX_test_df = X_test_df.drop(columns=[\"planet_id\"])\nX_test_scaled = scaler.transform(X_test_df)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:47.305722Z","iopub.execute_input":"2025-09-21T14:43:47.30596Z","iopub.status.idle":"2025-09-21T14:43:47.379339Z","shell.execute_reply.started":"2025-09-21T14:43:47.30594Z","shell.execute_reply":"2025-09-21T14:43:47.37848Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Predict the mean and get the std. Std calculation could have been further improved if I had additional time. Other methods such as bootstrapping, quantile regression, and residual based methods all performed worse, or I simply needed more time explore their implementations. ","metadata":{}},{"cell_type":"code","source":"# Predict means\nridge_mean = ridge.predict(X_test_scaled)\nxgb_mean = xgb.predict(X_test_scaled)\ngpr_mean, gpr_std = gpr.predict(X_test_scaled, return_std=True)\n\n#weighted mean\nensemble_mean = (gpr_mean*2 + xgb_mean*4 + ridge_mean) / 7\n\n# weighted uncertainty (conservative aggregation)\nridge_std = np.std(ridge_mean) * np.ones_like(ridge_mean)\nxgb_std = np.std(xgb_mean) * np.ones_like(xgb_mean)\n\nensemble_std = (gpr_std*2 + xgb_std*4 + ridge_std*2) / 8","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T16:23:00.93871Z","iopub.execute_input":"2025-09-21T16:23:00.938985Z","iopub.status.idle":"2025-09-21T16:23:01.069782Z","shell.execute_reply.started":"2025-09-21T16:23:00.938966Z","shell.execute_reply":"2025-09-21T16:23:01.066356Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"X_train = load(\"/kaggle/input/ensemble-models/X_train.joblib\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:47.415261Z","iopub.execute_input":"2025-09-21T14:43:47.415576Z","iopub.status.idle":"2025-09-21T14:43:47.428042Z","shell.execute_reply.started":"2025-09-21T14:43:47.415552Z","shell.execute_reply":"2025-09-21T14:43:47.42727Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Create and score submit sumbmission","metadata":{}},{"cell_type":"markdown","source":"Training ensemble mean and std was saved for use in the following functions.","metadata":{}},{"cell_type":"code","source":"ensemble_train_mean = load(\"/kaggle/input/ensemble-models/ensemble_train_mean2 (1).joblib\")\nensemble_train_std = load(\"/kaggle/input/ensemble-models/ensemble_train_std2.joblib\")\ny_val = load(\"/kaggle/input/ensemble-models/y_val.joblib\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T14:43:47.48733Z","iopub.execute_input":"2025-09-21T14:43:47.487761Z","iopub.status.idle":"2025-09-21T14:43:47.57159Z","shell.execute_reply.started":"2025-09-21T14:43:47.487734Z","shell.execute_reply":"2025-09-21T14:43:47.570808Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def gll_score(y, mu, sigma):\n    # gaussian log-likelihood per sample\n    return -0.5 * (np.log(2*np.pi) + 2*np.log(sigma) + ((y-mu)**2)/(sigma**2))\n\ndef gll_mean(y, mu, sigma):\n    return np.mean(gll_score(y, mu, sigma))\n\n# grid search on lambda, epsilon\ndef find_best_scale(y_val, mu_val, sigma_val, lambdas=None, epsilons=None):\n    lambdas = np.concatenate([\n        np.linspace(0.01, 0.1, 20),   # very conservative\n        np.linspace(0.1, 1.0, 40),    # conservative  \n        np.linspace(1.0, 5.0, 40),    # normal\n        np.linspace(5.0, 50.0, 20)    # aggressive\n    ])\n    epsilons = [1e-8, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]\n    \n    best_grid = (-1e9, None, None)\n    \n    for lam in lambdas:\n        for eps in epsilons:\n            sigma_cal = np.maximum(eps, lam * sigma_val)\n            score = gll_mean(y_val, mu_val, sigma_cal)\n            if score > best_grid[0]:\n                best_grid = (score, (lam, eps))\n\n    return best_grid\n\n# usage\nbest_score, (best_lambda, best_eps) = find_best_scale(y_val, ensemble_train_mean, ensemble_train_std)\nsigma_test_cal = np.maximum(best_eps, best_lambda * ensemble_std)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-21T16:33:51.971288Z","iopub.execute_input":"2025-09-21T16:33:51.971621Z","iopub.status.idle":"2025-09-21T16:33:52.757957Z","shell.execute_reply.started":"2025-09-21T16:33:51.971595Z","shell.execute_reply":"2025-09-21T16:33:52.756991Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def create_submission(predictions, uncertainties, planet_ids, sample_path=\"/kaggle/input/ariel-data-challenge-2025/sample_submission.csv\", output_path=\"submission.csv\"):\n    sample_columns = pd.read_csv(sample_path, index_col=\"planet_id\").columns\n\n    submission_df = pd.DataFrame(\n        np.concatenate([predictions, uncertainties], axis=1),\n        columns=sample_columns,\n        index=[int(pid) for pid in planet_ids]\n    )\n\n    submission_df.index.name = \"planet_id\"\n    submission_df.reset_index(inplace=True)\n\n    submission_df.to_csv(output_path, index=False)\n    print(f\"Submission saved to {output_path}\")\n    return submission_df\n\ncreate_submission(ensemble_mean, sigma_test_cal, planet_test_ids)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:16:42.169473Z","iopub.execute_input":"2025-09-24T21:16:42.170123Z","iopub.status.idle":"2025-09-24T21:16:42.18254Z","shell.execute_reply.started":"2025-09-24T21:16:42.170084Z","shell.execute_reply":"2025-09-24T21:16:42.181151Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Score the submission. This function wasn't incredibly accurate to the acutal competition score that would be recieved on attempts, but it at least gave a ballpark idea if the GLL score increased between attempts, and helped me avoid scores of 0.000. I would attribute the difference mainly to the difference in the training data vs. the hidden test set, my calibration in the functions above, and the fact that my method was sensitive to outliers, possibly too specific to outliers in the training data, and thus gave a better local score. ","metadata":{}},{"cell_type":"code","source":"def competition_score_fixed(solution, submission, naive_mean=0.0025517145902829823, \n                           naive_sigma=0.0017261973793536417, sigma_true=1e-5):\n    \"\"\"\n    naive, sigma true from last years winners\n    \"\"\"    \n    n_wavelengths = 283\n    \n    # Remove planet_id columns if present\n    solution_clean = solution.copy()\n    print(\"Solution shape:\", solution.shape)\n    submission_clean = submission.copy()\n    print(\"Submission shape:\", submission.shape)\n    \n    if 'planet_id' in solution_clean.columns:\n        solution_clean = solution_clean.drop('planet_id', axis=1)\n    if 'planet_id' in submission_clean.columns:\n        submission_clean = submission_clean.drop('planet_id', axis=1)\n    \n    # Extract predictions and uncertainties correctly\n    y_pred = submission_clean.iloc[:, :n_wavelengths].values\n    sigma_pred = submission_clean.iloc[:, n_wavelengths:].values  # FIX: Get uncertainties, not predictions\n    \n    # ensure non-zero sigma\n    sigma_pred = np.clip(sigma_pred, a_min=1e-15, a_max=None)\n    \n    # ground truth\n    #y_true = solution_clean.values\n    y_true = solution_clean.values\n    \n    # calculate GLLs\n    GLL_pred = np.sum(stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n    GLL_true = np.sum(stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true))\n    GLL_mean = np.sum(stats.norm.logpdf(y_true, loc=naive_mean, scale=naive_sigma))\n    \n    # normalize the score\n    submit_score = (GLL_pred - GLL_mean) / (GLL_true - GLL_mean)\n    \n    print(f\"GLL_pred: {GLL_pred:.2f}\")\n    print(f\"GLL_true: {GLL_true:.2f}\")\n    print(f\"GLL_mean: {GLL_mean:.2f}\")\n    print(f\"Raw score: {submit_score:.6f}\")\n    print(GLL_pred > GLL_mean)\n    \n    return float(np.clip(submit_score, 0.0, 1.0))\n\n\n#solution = train_df[:330]\nsolution = test_df.iloc[:, :284]\nsubmission = pd.read_csv(\"/kaggle/working/submission.csv\")\nscore = competition_score_fixed(solution, submission)\nprint(f\"Final Score: {score:.6f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-09-24T21:14:09.489928Z","iopub.execute_input":"2025-09-24T21:14:09.490206Z","iopub.status.idle":"2025-09-24T21:14:09.595095Z","shell.execute_reply.started":"2025-09-24T21:14:09.490176Z","shell.execute_reply":"2025-09-24T21:14:09.593774Z"}},"outputs":[],"execution_count":null}]}