{"metadata":{"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30746,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.13"},"papermill":{"default_parameters":{},"duration":1087.424982,"end_time":"2024-08-03T13:00:18.164826","environment_variables":{},"exception":null,"input_path":"__notebook__.ipynb","output_path":"__notebook__.ipynb","parameters":{},"start_time":"2024-08-03T12:42:10.739844","version":"2.5.0"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ariel Data Challenge 2024: Introductory model: training with more feature\n\nIn this notebook, we show how to train and cross-validate a model. At the end, we save the model so that it can be used for inference in the separate notebook. Origional version is from [ADC24 Intro training](https://www.kaggle.com/code/ambrosm/adc24-intro-training).\n","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","papermill":{"duration":0.009657,"end_time":"2024-08-03T12:42:14.074227","exception":false,"start_time":"2024-08-03T12:42:14.06457","status":"completed"},"tags":[]}},{"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\n\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error","metadata":{"_kg_hide-input":true,"papermill":{"duration":3.023083,"end_time":"2024-08-03T12:42:17.107326","exception":false,"start_time":"2024-08-03T12:42:14.084243","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T10:27:09.442466Z","iopub.execute_input":"2024-08-04T10:27:09.442997Z","iopub.status.idle":"2024-08-04T10:27:09.45567Z","shell.execute_reply.started":"2024-08-04T10:27:09.44295Z","shell.execute_reply":"2024-08-04T10:27:09.453632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# A look at the data\n\nWe start by reading the metadata:","metadata":{"papermill":{"duration":0.008428,"end_time":"2024-08-03T12:42:17.180297","exception":false,"start_time":"2024-08-03T12:42:17.171869","status":"completed"},"tags":[]}},{"cell_type":"code","source":"train_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_adc_info.csv',\n                           index_col='planet_id')\n# test_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n#                            index_col='planet_id')\ntrain_labels = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_labels.csv',\n                           index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2024/axis_info.parquet')\n","metadata":{"papermill":{"duration":0.188106,"end_time":"2024-08-03T12:42:17.37713","exception":false,"start_time":"2024-08-03T12:42:17.189024","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T09:53:37.065029Z","iopub.execute_input":"2024-08-04T09:53:37.06549Z","iopub.status.idle":"2024-08-04T09:53:37.365489Z","shell.execute_reply.started":"2024-08-04T09:53:37.065453Z","shell.execute_reply":"2024-08-04T09:53:37.363996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Some facts about the data:\n- We have 673 planets for training. These planets belong to two different stars.\n- There will be roughly 800 planets for testing (but the test data is hidden).\n- The competition is a multi-output regression task with 283 targets to predict.","metadata":{"papermill":{"duration":0.008486,"end_time":"2024-08-03T12:42:17.394634","exception":false,"start_time":"2024-08-03T12:42:17.386148","status":"completed"},"tags":[]}},{"cell_type":"markdown","source":"## The FGS1 data\n\nHaving read the metadata, we'll tackle the FGS1 data (Fine Guidance System). The FGS1 measurements consist of one file per planet (673 files for 673 planets for training). For now, we ignore the calibration files.\n\nEach file contains 135,000 rows of images taken at 0.1 second time steps. Each row is a 32\\*32 image at a single wavelength.\n\nWe read a sample file:","metadata":{"papermill":{"duration":0.008446,"end_time":"2024-08-03T12:42:17.411955","exception":false,"start_time":"2024-08-03T12:42:17.403509","status":"completed"},"tags":[]}},{"cell_type":"code","source":"planet_id = 14485303\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\nf_signal","metadata":{"papermill":{"duration":1.759363,"end_time":"2024-08-03T12:42:19.180097","exception":false,"start_time":"2024-08-03T12:42:17.420734","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T09:53:38.357806Z","iopub.execute_input":"2024-08-04T09:53:38.358191Z","iopub.status.idle":"2024-08-04T09:53:40.382055Z","shell.execute_reply.started":"2024-08-04T09:53:38.358163Z","shell.execute_reply":"2024-08-04T09:53:40.380822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Every row of the file corresponds to an image of a star. The images come in pairs, and the second image is lighter than the first one:","metadata":{"papermill":{"duration":0.009711,"end_time":"2024-08-03T12:42:19.199727","exception":false,"start_time":"2024-08-03T12:42:19.190016","status":"completed"},"tags":[]}},{"cell_type":"code","source":"_, (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":{"_kg_hide-input":true,"papermill":{"duration":0.559852,"end_time":"2024-08-03T12:42:19.769794","exception":false,"start_time":"2024-08-03T12:42:19.209942","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T09:54:46.155378Z","iopub.execute_input":"2024-08-04T09:54:46.155844Z","iopub.status.idle":"2024-08-04T09:54:46.979783Z","shell.execute_reply.started":"2024-08-04T09:54:46.155807Z","shell.execute_reply":"2024-08-04T09:54:46.97861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To see the time series, we first have to compute the difference between the even and the odd frames to get the net signal (67500 time steps). We then take the mean over all 1024 pixels. The net signal is very noisy, and we smoothen it by computing a moving average. The plot of the smoothened signal clearly shows that the signal intensity is reduced (i.e., the image gets darker) while the planet passes in front of the star (between time steps 23500 and 44000).\n\nThe left diagram shows a planet with a strong reduction of the signal intensity, the right diagram shows a planet with a weak reduction:","metadata":{"papermill":{"duration":0.011386,"end_time":"2024-08-03T12:42:19.791847","exception":false,"start_time":"2024-08-03T12:42:19.780461","status":"completed"},"tags":[]}},{"cell_type":"code","source":"_, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, figsize=(12, 4))\n\nplanet_id = 14485303\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n\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\n\nax1.set_title('FGS1: time series of planet with strong signal')\nax1.plot(net_signal, label='raw signal')\nax1.legend()\nax3.plot(smooth_signal, color='c', label='smoothened signal')\nax3.legend()\nax3.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax3.axvline(time_step, color='gray')\n\nplanet_id = 4249337798\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/FGS1_signal.parquet')\n\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\n\nax2.set_title('FGS1: time series of planet with weak signal')\nax2.plot(net_signal, label='raw signal')\nax2.legend()\nax4.plot(smooth_signal, color='c', label='smoothened signal')\nax4.legend()\nax4.set_xlabel('time step')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax4.axvline(time_step, color='gray')\n\n# plt.suptitle('FGS1 time series', y=0.96)\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-04T11:53:39.73513Z","iopub.execute_input":"2024-08-04T11:53:39.735573Z","iopub.status.idle":"2024-08-04T11:53:43.581762Z","shell.execute_reply.started":"2024-08-04T11:53:39.735538Z","shell.execute_reply":"2024-08-04T11:53:43.580456Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## The AIRS data\n\nAIRS is the other sensor of the satellite. It produces one file per planet as well. Each file contains 11,250 rows of images captured at constant time steps. Each 32 x 356 image has been flattened into 11392 columns.","metadata":{"papermill":{"duration":0.077595,"end_time":"2024-08-03T13:00:12.362301","exception":false,"start_time":"2024-08-03T13:00:12.284706","status":"completed"},"tags":[]}},{"cell_type":"code","source":"planet_id = 4249337798\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/train/{planet_id}/AIRS-CH0_signal.parquet')\na_signal","metadata":{"papermill":{"duration":2.313418,"end_time":"2024-08-03T13:00:14.755198","exception":false,"start_time":"2024-08-03T13:00:12.44178","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:55:20.139045Z","iopub.execute_input":"2024-08-04T11:55:20.140139Z","iopub.status.idle":"2024-08-04T11:55:22.915138Z","shell.execute_reply.started":"2024-08-04T11:55:20.140098Z","shell.execute_reply":"2024-08-04T11:55:22.913964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"a_signal = a_signal.values.reshape(11250, 32, 356)\n\nplt.figure(figsize=(10, 3))\nsns.heatmap(a_signal[1])\nplt.ylabel('spatial dimension')\nplt.xlabel('wavelength dimension')\nplt.show()","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.807236,"end_time":"2024-08-03T13:00:15.817499","exception":false,"start_time":"2024-08-03T13:00:15.010263","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:55:22.917083Z","iopub.execute_input":"2024-08-04T11:55:22.917443Z","iopub.status.idle":"2024-08-04T11:55:23.523336Z","shell.execute_reply.started":"2024-08-04T11:55:22.917413Z","shell.execute_reply":"2024-08-04T11:55:23.522136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The data again is a time series, and we can see how the star is obscured while the planet is passing in front of it.","metadata":{"papermill":{"duration":0.08019,"end_time":"2024-08-03T13:00:15.979865","exception":false,"start_time":"2024-08-03T13:00:15.899675","status":"completed"},"tags":[]}},{"cell_type":"code","source":"mean_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')\nax1.legend()\nax2.plot(smooth_signal, color='c', label='smoothened net signal')\nax2.legend()\nax2.set_xlabel('time')\nfor time_step in [20500, 23500, 44000, 47000]:\n    ax2.axvline(time_step * 11250 // 135000, color='gray')\nplt.suptitle('AIRS-CH0 time series', y=0.96)\nplt.show()\n","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.727242,"end_time":"2024-08-03T13:00:16.788841","exception":false,"start_time":"2024-08-03T13:00:16.061599","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:55:23.524891Z","iopub.execute_input":"2024-08-04T11:55:23.525274Z","iopub.status.idle":"2024-08-04T11:55:24.131163Z","shell.execute_reply.started":"2024-08-04T11:55:23.525234Z","shell.execute_reply":"2024-08-04T11:55:24.130046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Reading the data\n\nWe now read the FGS1 data and the AIRS-CH0 data for all 673 training planets. As the dataset doesn't fit into RAM completely, we keep only two one-dimensional time series for every planet. At the end we'll have\n1. A time series with 67500 steps per planet taken from the FGS1 data, and\n2. A time series with 5625 steps per planet taken from the AIRS-CH0 data.\n\nWe use the Jupyter `%%writefile` cell magic to save the function code to a file. This ensures that the inference notebook will process the test data in exactly the same way as this notebook processes the training data.","metadata":{"papermill":{"duration":0.011189,"end_time":"2024-08-03T12:42:20.928036","exception":false,"start_time":"2024-08-03T12:42:20.916847","status":"completed"},"tags":[]}},{"cell_type":"code","source":"%%writefile f_read_and_preprocess.py\n\ndef f_read_and_preprocess(dataset, adc_info, planet_ids):\n    \"\"\"Read the FGS1 files for all planet_ids and extract the time series.\n    \n    Parameters\n    dataset: 'train' or 'test'\n    adc_info: metadata dataframe, either train_adc_info or test_adc_info\n    planet_ids: list of planet ids\n    \n    Returns\n    dataframe with one row per planet_id and 67500 values per row\n    \n    \"\"\"\n    f_raw_train = np.full((len(planet_ids), 67500), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        f_signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/FGS1_signal.parquet')\n        mean_signal = f_signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / 1024 # mean over the 32*32 pixels\n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n        f_raw_train[i] = net_signal\n    return f_raw_train\n    ","metadata":{"papermill":{"duration":1067.258124,"end_time":"2024-08-03T13:00:08.197814","exception":false,"start_time":"2024-08-03T12:42:20.93969","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T12:27:23.283354Z","iopub.execute_input":"2024-08-04T12:27:23.283949Z","iopub.status.idle":"2024-08-04T12:27:23.294819Z","shell.execute_reply.started":"2024-08-04T12:27:23.283906Z","shell.execute_reply":"2024-08-04T12:27:23.29351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nexec(open('f_read_and_preprocess.py', 'r').read())\nf_raw_train = f_read_and_preprocess('train', train_adc_info, train_labels.index)\nwith open('f_raw_train.pickle', 'wb') as f:\n    pickle.dump(f_raw_train, f)\n","metadata":{"papermill":{"duration":1067.258124,"end_time":"2024-08-03T13:00:08.197814","exception":false,"start_time":"2024-08-03T12:42:20.93969","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T10:16:27.461497Z","iopub.execute_input":"2024-08-04T10:16:27.461961Z","iopub.status.idle":"2024-08-04T10:27:09.214198Z","shell.execute_reply.started":"2024-08-04T10:16:27.461928Z","shell.execute_reply":"2024-08-04T10:27:09.212147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile a_read_and_preprocess.py\ndef a_read_and_preprocess(dataset, adc_info, planet_ids):\n    \"\"\"Read the AIRS-CH0 files for all planet_ids and extract the time series.\n    \n    Parameters\n    dataset: 'train' or 'test'\n    adc_info: metadata dataframe, either train_adc_info or test_adc_info\n    planet_ids: list of planet ids\n    \n    Returns\n    dataframe with one row per planet_id and 5625 values per row\n    \n    \"\"\"\n    a_raw_train = np.full((len(planet_ids), 5625), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        signal = pl.read_parquet(f'/kaggle/input/ariel-data-challenge-2024/{dataset}/{planet_id}/AIRS-CH0_signal.parquet')\n        mean_signal = signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / (32*356) # mean over the 32*356 pixels\n        net_signal = mean_signal[1::2] - mean_signal[0::2]\n        a_raw_train[i] = net_signal\n    return a_raw_train\n    ","metadata":{"execution":{"iopub.status.busy":"2024-08-04T12:27:27.791719Z","iopub.execute_input":"2024-08-04T12:27:27.792857Z","iopub.status.idle":"2024-08-04T12:27:27.800511Z","shell.execute_reply.started":"2024-08-04T12:27:27.792811Z","shell.execute_reply":"2024-08-04T12:27:27.799003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nexec(open('a_read_and_preprocess.py', 'r').read())\na_raw_train = a_read_and_preprocess('train', train_adc_info, train_labels.index)\nwith open('a_raw_train.pickle', 'wb') as f:\n    pickle.dump(a_raw_train, f)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T10:49:43.17432Z","iopub.execute_input":"2024-08-04T10:49:43.1748Z","iopub.status.idle":"2024-08-04T11:10:07.5598Z","shell.execute_reply.started":"2024-08-04T10:49:43.174748Z","shell.execute_reply":"2024-08-04T11:10:07.557469Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As a plausibility check, we plot the means of all time series:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(6, 2))\nplt.plot(f_raw_train.mean(axis=0))\nfor time_step in [20500, 23500, 44000, 47000]:\n    plt.axvline(time_step, color='gray')\nplt.xlabel('time step')\nplt.title('FGS1: Overall mean')\nplt.show()\n\nplt.figure(figsize=(6, 2))\nplt.plot(a_raw_train.mean(axis=0))\nfor time_step in [20500, 23500, 44000, 47000]:\n    plt.axvline(time_step * 11250 // 135000, color='gray')\nplt.xlabel('time step')\nplt.title('AIRS-CH0: Overall mean')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:21:05.261032Z","iopub.execute_input":"2024-08-04T11:21:05.262209Z","iopub.status.idle":"2024-08-04T11:21:05.753651Z","shell.execute_reply.started":"2024-08-04T11:21:05.262167Z","shell.execute_reply":"2024-08-04T11:21:05.752301Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature engineering\n\nIn our feature engineering process, we aim to capture the key characteristics of how the brightness of stars changes when obscured by planets. Specifically, we want to quantify how much darker the images become during planetary transits. The time series data shows that planets reduce the brightness of stars by approximately 0.2% on average (e.g., from 228.2 to 227.6 or from 1371 to 1368).\n\nTo achieve this, we have extracted a variety of features:\n\nBasic Statistical Features: These features capture the overall distribution of the data, including:\n\nRelative Reduction: The difference in brightness between the obscured and unobscured regions, normalized by the unobscured brightness.\nSignal-to-Noise Ratio (SNR): The ratio of unobscured brightness to the overall standard deviation, indicating the quality of the signal.\nVariance, Skewness, Kurtosis: These statistical metrics describe the spread, symmetry, and peak characteristics of the brightness data.\nSubset Features: These features focus on specific time intervals, particularly the regions where the planet partially obscures the star:\n\nHalf Reduction: We calculate the relative reduction in brightness for the two half-obscured regions of the time series, providing a more detailed view of the transit effect.\nSliding Window Features: By applying a sliding window approach, we capture the local variations in the time series data:\n\nFor each window, we compute the mean, standard deviation, minimum, and maximum brightness values, which help us understand how brightness fluctuates over time.","metadata":{}},{"cell_type":"code","source":"%%writefile feature_engineering.py\n\nimport numpy as np\nimport pandas as pd\n\ndef feature_engineering(f_raw, a_raw, window_size=50, step_size=25):\n    \"\"\"Create a dataframe with combined features from the raw data, including sliding window and time-series statistics.\n    \n    Parameters:\n    f_raw: ndarray of shape (n_planets, 67500)\n    a_raw: ndarray of shape (n_planets, 5625)\n    window_size: int, size of the sliding window for time-series statistics\n    step_size: int, step size for the sliding window\n    \n    Return value:\n    df: DataFrame of shape (n_planets, several features)\n    \"\"\"\n    f_obscured = f_raw[:, 23500:44000].mean(axis=1)\n    f_unobscured = (f_raw[:, :20500].mean(axis=1) + f_raw[:, 47000:].mean(axis=1)) / 2\n    f_relative_reduction = (f_unobscured - f_obscured) / f_unobscured\n    f_std_dev = f_raw.std(axis=1)\n    f_signal_to_noise = f_unobscured / f_std_dev\n\n    a_obscured = a_raw[:, 1958:3666].mean(axis=1)\n    a_unobscured = (a_raw[:, :1708].mean(axis=1) + a_raw[:, 3916:].mean(axis=1)) / 2\n    a_relative_reduction = (a_unobscured - a_obscured) / a_unobscured\n    a_std_dev = a_raw.std(axis=1)\n    a_signal_to_noise = a_unobscured / a_std_dev\n\n    f_variance = f_raw.var(axis=1)\n    a_variance = a_raw.var(axis=1)\n    \n    f_skewness = pd.DataFrame(f_raw).skew(axis=1).values\n    a_skewness = pd.DataFrame(a_raw).skew(axis=1).values\n\n    f_kurtosis = pd.DataFrame(f_raw).kurtosis(axis=1).values\n    a_kurtosis = pd.DataFrame(a_raw).kurtosis(axis=1).values\n    \n    f_half_obscured1 = f_raw[:, 20500:23500].mean(axis=1)\n    f_half_obscured2 = f_raw[:, 44000:47000].mean(axis=1)\n    f_half_reduction1 = (f_unobscured - f_half_obscured1) / f_unobscured\n    f_half_reduction2 = (f_unobscured - f_half_obscured2) / f_unobscured\n\n    a_half_obscured1 = a_raw[:, 1708:1958].mean(axis=1)\n    a_half_obscured2 = a_raw[:, 3666:3916].mean(axis=1)\n    a_half_reduction1 = (a_unobscured - a_half_obscured1) / a_unobscured\n    a_half_reduction2 = (a_unobscured - a_half_obscured2) / a_unobscured\n\n    # Sliding window features\n    def sliding_window_features(data, window_size, step_size):\n        features = []\n        max_index = data.shape[1]\n        for start in range(0, max_index - window_size + 1, step_size):\n            end = start + window_size\n            window = data[:, start:end]\n            features.append([\n                np.mean(window, axis=1),\n                np.std(window, axis=1),\n                np.min(window, axis=1),\n                np.max(window, axis=1)\n            ])\n        if features:\n            return np.vstack(features).T  # Stack vertically and transpose to get the correct shape\n        else:\n            return np.empty((data.shape[0], 0))  # Return empty array with correct shape\n    \n    f_sliding_features = sliding_window_features(f_raw, window_size, step_size)\n    a_sliding_features = sliding_window_features(a_raw, window_size, step_size)\n\n\n    print(f'f_sliding_features.shape: {f_sliding_features.shape}')\n    print(f'a_sliding_features.shape: {a_sliding_features.shape}')\n\n\n    df = pd.DataFrame({\n        'f_relative_reduction': f_relative_reduction,\n        'f_signal_to_noise': f_signal_to_noise,\n        'f_variance': f_variance,\n        'f_skewness': f_skewness,\n        'f_kurtosis': f_kurtosis,\n        'a_relative_reduction': a_relative_reduction,\n        'a_signal_to_noise': a_signal_to_noise,\n        'a_variance': a_variance,\n        'a_skewness': a_skewness,\n        'a_kurtosis': a_kurtosis,\n        'f_half_reduction1': f_half_reduction1,\n        'f_half_reduction2': f_half_reduction2,\n        'a_half_reduction1': a_half_reduction1,\n        'a_half_reduction2': a_half_reduction2\n    })\n\n\n    if f_sliding_features.size > 0:\n        f_sliding_df = pd.DataFrame(f_sliding_features, columns=[f'f_slide_{i}' for i in range(f_sliding_features.shape[1])])\n        df = pd.concat([df, f_sliding_df], axis=1)\n\n    if a_sliding_features.size > 0:\n        a_sliding_df = pd.DataFrame(a_sliding_features, columns=[f'a_slide_{i}' for i in range(a_sliding_features.shape[1])])\n        df = pd.concat([df, a_sliding_df], axis=1)\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:49:32.362734Z","iopub.execute_input":"2024-08-04T11:49:32.363268Z","iopub.status.idle":"2024-08-04T11:49:32.371284Z","shell.execute_reply.started":"2024-08-04T11:49:32.363232Z","shell.execute_reply":"2024-08-04T11:49:32.370079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('feature_engineering.py', 'r').read())\n\ntrain = feature_engineering(f_raw_train, a_raw_train)","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:33:15.805608Z","iopub.execute_input":"2024-08-04T11:33:15.806053Z","iopub.status.idle":"2024-08-04T11:33:15.833728Z","shell.execute_reply.started":"2024-08-04T11:33:15.806018Z","shell.execute_reply":"2024-08-04T11:33:15.832575Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following scatterplot shows a strong correlation between the signal reduction when the planet is in front of the star and one of the targets we want to predict. The other targets have a similarly high correlation.","metadata":{"papermill":{"duration":0.073746,"end_time":"2024-08-03T13:00:08.349433","exception":false,"start_time":"2024-08-03T13:00:08.275687","status":"completed"},"tags":[]}},{"cell_type":"code","source":"\ncolor_array = np.array(plt.rcParams['axes.prop_cycle'].by_key()['color'])\nplt.scatter(train.a_relative_reduction, train_labels.wl_1, s=15, alpha=0.5,\n            c=color_array[train_adc_info.star])\nplt.xlabel('relative signal reduction when planet is in front')\nplt.ylabel('target')\nplt.title('Correlation between relative signal reduction and target')\nplt.gca().set_aspect('equal')\npoints = [plt.Line2D([0], [0], label=f'star {i}', marker='o', markersize=3,\n         markeredgecolor=color_array[i], markerfacecolor=color_array[i], linestyle='') for i in range(2)]\n\nplt.legend(handles=points)\nplt.show()","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.478644,"end_time":"2024-08-03T13:00:08.902093","exception":false,"start_time":"2024-08-03T13:00:08.423449","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:18.162145Z","iopub.execute_input":"2024-08-04T11:33:18.162554Z","iopub.status.idle":"2024-08-04T11:33:18.521622Z","shell.execute_reply.started":"2024-08-04T11:33:18.162523Z","shell.execute_reply":"2024-08-04T11:33:18.520319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The model and the cross-validation\n\nTo keep things simple, we predict the targets with ridge regression.\n\nWe are interested in three cross-validation metrics:\n1. The R2 score is above 0.9, which confirms the correlation we've seen in the scatterplot.\n2. The root mean squared error will be the predicted uncertainty.\n3. The competition metric gives an indication of the leaderboard score. Unfortunately the competition metric depends on the value of `sigma_true`, which I don't know.","metadata":{"papermill":{"duration":0.076791,"end_time":"2024-08-03T13:00:09.056736","exception":false,"start_time":"2024-08-03T13:00:08.979945","status":"completed"},"tags":[]}},{"cell_type":"code","source":"model = Ridge(alpha=1e-12)\n\noof_pred = cross_val_predict(model, train, train_labels)\n\nprint(f\"# R2 score: {r2_score(train_labels, oof_pred):.3f}\")\nsigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\nprint(f\"# Root mean squared error: {sigma_pred:.6f}\")\n# R2 score: 0.971\n# Root mean squared error: 0.000293","metadata":{"papermill":{"duration":0.700113,"end_time":"2024-08-03T13:00:09.837053","exception":false,"start_time":"2024-08-03T13:00:09.13694","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:20.382678Z","iopub.execute_input":"2024-08-04T11:33:20.383118Z","iopub.status.idle":"2024-08-04T11:33:20.508765Z","shell.execute_reply.started":"2024-08-04T11:33:20.383083Z","shell.execute_reply":"2024-08-04T11:33:20.507443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"col = 1\nplt.scatter(oof_pred[:,col], train_labels.iloc[:,col], s=15, c='lightgreen')\nplt.gca().set_aspect('equal')\nplt.xlabel('y_pred')\nplt.ylabel('y_true')\nplt.title('Comparing y_true and y_pred')\nplt.show()","metadata":{"papermill":{"duration":0.700113,"end_time":"2024-08-03T13:00:09.837053","exception":false,"start_time":"2024-08-03T13:00:09.13694","status":"completed"},"tags":[],"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-08-04T11:33:24.335234Z","iopub.execute_input":"2024-08-04T11:33:24.33565Z","iopub.status.idle":"2024-08-04T11:33:24.628829Z","shell.execute_reply.started":"2024-08-04T11:33:24.335622Z","shell.execute_reply":"2024-08-04T11:33:24.627451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile competition_score.py\n# Adapted from https://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\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    This is a Gaussian Log Likelihood based metric. For a submission, which contains the predicted mean (x_hat) and variance (x_hat_std),\n    we calculate the Gaussian Log-likelihood (GLL) value to the provided ground truth (x). We treat each pair of x_hat,\n    x_hat_std as a 1D gaussian, meaning there will be 283 1D gaussian distributions, hence 283 values for each test spectrum,\n    the GLL value for one spectrum is the sum of all of them.\n\n    Inputs:\n        - solution: Ground Truth spectra (from test set)\n            - shape: (nsamples, n_wavelengths)\n        - submission: Predicted spectra and errors (from participants)\n            - shape: (nsamples, n_wavelengths*2)\n        naive_mean: (float) mean from the train set.\n        naive_sigma: (float) standard deviation from the train set.\n        sigma_true: (float) essentially sets the scale of the outputs.\n    '''\n\n    del solution[row_id_column_name]\n    del submission[row_id_column_name]\n\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\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 a number')\n\n    n_wavelengths = len(solution.columns)\n    if len(submission.columns) != n_wavelengths*2:\n        raise ParticipantVisibleError('Wrong number of columns in the submission')\n\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    # Set a non-zero minimum sigma pred to prevent division by zero errors.\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=10**-15, a_max=None)\n    y_true = solution.values\n\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 * np.ones_like(y_true)))\n    GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean * np.ones_like(y_true), scale=naive_sigma * np.ones_like(y_true)))\n\n    submit_score = (GLL_pred - GLL_mean)/(GLL_true - GLL_mean)\n    return float(np.clip(submit_score, 0.0, 1.0))","metadata":{"_kg_hide-input":true,"papermill":{"duration":0.029243,"end_time":"2024-08-03T12:42:17.162838","exception":false,"start_time":"2024-08-03T12:42:17.133595","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:25.777887Z","iopub.execute_input":"2024-08-04T11:33:25.778327Z","iopub.status.idle":"2024-08-04T11:33:25.787298Z","shell.execute_reply.started":"2024-08-04T11:33:25.778293Z","shell.execute_reply":"2024-08-04T11:33:25.78591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile postprocessing.py\n\ndef postprocessing(pred_array, index, sigma_pred):\n    \"\"\"Create a submission dataframe from its components\n    \n    Parameters:\n    pred_array: ndarray of shape (n_samples, 283)\n    index: pandas.Index of length n_samples with name 'planet_id'\n    sigma_pred: float\n    \n    Return value:\n    df: DataFrame of shape (n_samples, 566) with planet_id as index\n    \"\"\"\n    return pd.concat([pd.DataFrame(pred_array.clip(0, None), index=index, columns=wavelengths.columns),\n                      pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i}\" for i in range(1, 284)])],\n                     axis=1)\n\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T11:33:27.453077Z","iopub.execute_input":"2024-08-04T11:33:27.453955Z","iopub.status.idle":"2024-08-04T11:33:27.46076Z","shell.execute_reply.started":"2024-08-04T11:33:27.453914Z","shell.execute_reply":"2024-08-04T11:33:27.459579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('competition_score.py', 'r').read())\nexec(open('postprocessing.py', 'r').read())\n\noof_df = postprocessing(oof_pred, train_adc_info.index, sigma_pred)\ndisplay(oof_df)\n\ngll_score = competition_score(train_labels.copy().reset_index(),\n                              oof_df.copy().reset_index(),\n                              naive_mean=train_labels.values.mean(),\n                              naive_sigma=train_labels.values.std(),\n                              sigma_true=0.000003)\nprint(f\"# Estimated competition score: {gll_score:.3f}\")\n# Estimated competition score: 0.259","metadata":{"papermill":{"duration":0.231627,"end_time":"2024-08-03T13:00:10.149224","exception":false,"start_time":"2024-08-03T13:00:09.917597","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T11:33:28.166304Z","iopub.execute_input":"2024-08-04T11:33:28.166706Z","iopub.status.idle":"2024-08-04T11:33:28.266475Z","shell.execute_reply.started":"2024-08-04T11:33:28.166678Z","shell.execute_reply":"2024-08-04T11:33:28.265186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Refitting and saving the model","metadata":{}},{"cell_type":"code","source":"# Refit the model to the full dataset\nmodel.fit(train, train_labels)\nwith open('model.pickle', 'wb') as f:\n    pickle.dump(model, f)\nwith open('sigma_pred.pickle', 'wb') as f:\n    pickle.dump(sigma_pred, f)\n","metadata":{"execution":{"iopub.status.busy":"2024-08-04T13:15:17.228909Z","iopub.execute_input":"2024-08-04T13:15:17.229421Z","iopub.status.idle":"2024-08-04T13:15:17.253336Z","shell.execute_reply.started":"2024-08-04T13:15:17.229387Z","shell.execute_reply":"2024-08-04T13:15:17.251897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submission","metadata":{"papermill":{"duration":0.077479,"end_time":"2024-08-03T13:00:10.305108","exception":false,"start_time":"2024-08-03T13:00:10.227629","status":"completed"},"tags":[]}},{"cell_type":"code","source":"# Load the data\ntest_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/test_adc_info.csv',\n                           index_col='planet_id')\nsample_submission = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/sample_submission.csv',\n                                index_col='planet_id')\nf_raw_test = f_read_and_preprocess('test', test_adc_info, sample_submission.index)\na_raw_test = a_read_and_preprocess('test', test_adc_info, sample_submission.index)\ntest = feature_engineering(f_raw_test, a_raw_test)\n\n# Load the model\nwith open('model.pickle', 'rb') as f:\n    model = pickle.load(f)\nwith open('sigma_pred.pickle', 'rb') as f:\n    sigma_pred = pickle.load(f)\n\n# Predict\ntest_pred = model.predict(test)\n\n# Package into submission file\nsub_df = postprocessing(test_pred, sample_submission.index, sigma_pred)\ndisplay(sub_df)\nsub_df.to_csv('submission.csv')\n#!head submission.csv","metadata":{"papermill":{"duration":1.822768,"end_time":"2024-08-03T13:00:12.204995","exception":false,"start_time":"2024-08-03T13:00:10.382227","status":"completed"},"tags":[],"execution":{"iopub.status.busy":"2024-08-04T13:15:48.046438Z","iopub.execute_input":"2024-08-04T13:15:48.046918Z","iopub.status.idle":"2024-08-04T13:15:51.037458Z","shell.execute_reply.started":"2024-08-04T13:15:48.046886Z","shell.execute_reply":"2024-08-04T13:15:51.036278Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}