{"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 (ipykernel)","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.8.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:  training\n\nThe original version of the Train code is at:  [ADC24 Intro training](https://www.kaggle.com/code/ambrosm/adc24-intro-training).\nThe original version of the infer code is at: [ADC24 Intro inference](https://www.kaggle.com/code/ambrosm/adc24-intro-inference).\n\nI have added some new features, which resulted in a slight performance improvement(LB:0.388), \n\nThis is the training code。","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_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('../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('../input/ariel-data-challenge-2024/train_labels.csv',\n                           index_col='planet_id')\nwavelengths = pd.read_csv('../input/ariel-data-challenge-2024/wavelengths.csv')\naxis_info = pd.read_parquet('../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_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.\n\n我们有673个行星用于训练。这些行星属于两颗不同的恒星。\n测试数据大约有800个行星（但测试数据是隐藏的）。\n这场比赛是一个多输出回归任务，需要预测283个目标。","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:\n\n阅读了元数据之后，我们将处理FGS1数据（精细导星系统）。FGS1的测量数据由每个行星一个文件组成（训练数据中为673个行星，共有673个文件）。目前，我们先忽略校准文件。\n\n每个文件包含135,000行图像数据，每0.1秒拍摄一次。每一行是单一波长下的32x32像素图像。","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'../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_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:\n\n文件的每一行对应一颗恒星的图像。这些图像成对出现，第二幅图像比第一幅亮一些：","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_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:\n\n要观察时间序列，我们首先需要计算偶数帧和奇数帧之间的差异，以获得净信号（67500个时间步）。然后我们对所有1024个像素取平均值。净信号非常嘈杂，我们通过计算移动平均值来平滑它。平滑后的信号图清楚地显示，在行星经过恒星前方时（在时间步23500到44000之间），信号强度会降低（即，图像变得更暗）。\n\n左侧的图表显示了一个信号强度大幅减少的行星，右侧的图表显示了一个信号强度减少较弱的行星：","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'../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'../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,"tags":[]},"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.\n\nAIRS 是卫星上的另一个传感器。它同样为每个行星生成一个文件。每个文件包含11,250行图像，这些图像以恒定的时间间隔被捕捉。每个32x356像素的图像被展平为11,392列。\n\n这段描述说明AIRS传感器的数据结构与之前提到的FGS1传感器的数据结构有所不同。AIRS的数据包括：\n\n每个行星一个文件，与FGS1相同。\n\n文件中包含11,250行图像，这意味着时间序列的长度比FGS1的135,000行要短。\n\n每张图像原始尺寸为32x356像素，但已经被转换成了一维数据，即11,392列，这可能是通过将二维图像展平为一维数组来实现的。\n\n这种数据结构可能意味着在处理AIRS数据时，需要考虑图像的一维表示，以及如何处理和分析这些数据以提取有用的特征或信息。\n\n","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'../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_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_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.\n数据再次表现为时间序列，我们可以观察到当行星从恒星前面经过时，恒星是如何被遮挡的。\n\n在这种情况下，时间序列数据通常用于分析行星凌星现象，即行星在恒星前面经过时，恒星的亮度会发生变化。\n\n这种变化可以通过时间序列分析来观察，通常表现为恒星亮度的周期性下降。这种亮度下降的幅度和持续时间可以提供有关行星大小、轨道和大气成分等信息。\n\n通过分析这些数据，科学家可以更好地了解行星的物理特性和其与恒星的相互作用。","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_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.\n\n我们现在读取所有673个训练行星的FGS1数据和AIRS-CH0数据。由于数据集不能完全装入RAM，我们为每个行星只保留两个一维时间序列。最后我们将得到：\n\n来自FGS1数据的每个行星有67500个步骤的时间序列，\n\n来自AIRS-CH0数据的每个行星有5625个步骤的时间序列。\n","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'../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_count":null,"outputs":[]},{"cell_type":"code","source":"#%%time\n#exec(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_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'../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":{"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#%%time\n#exec(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":{"tags":[]},"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":{"_kg_hide-input":true,"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature engineering\n\nWe want to know how much darker the images get when the planet obscures the star. The time series diagrams above show that the planets reduce the brightness of the stars (on average) by 0.2 % (from 228.2 to 227.6 or from 1371 to 1368).\n\n我们想知道当行星遮挡恒星时，图像会变暗多少。上面的时间序列图显示，行星（平均）使恒星的亮度降低了0.2%（从228.2降至227.6，或从1371降至1368）。","metadata":{}},{"cell_type":"markdown","source":"#%%writefile feature_engineering.py\n\ndef feature_engineering(f_raw, a_raw):\n    \"\"\"Create a dataframe with two features from the raw data.\n    \n    Parameters:\n    f_raw: ndarray of shape (n_planets, 67500)\n    a_raw: ndarray of shape (n_planets, 5625)\n    \n    Return value:\n    df: DataFrame of shape (n_planets, 2)\n    \"\"\"\n    obscured = f_raw[:, 23500:44000].mean(axis=1)\n    unobscured = (f_raw[:, :20500].mean(axis=1) + f_raw[:, 47000:].mean(axis=1)) / 2\n    f_relative_reduction = (unobscured - obscured) / unobscured\n    \n    half_obscured1 = f_raw[:, 20500:23500].mean(axis=1)\n    half_obscured2 = f_raw[:, 44000:47000].mean(axis=1)\n    f_half_reduction1 = (unobscured - half_obscured1) / unobscured\n    f_half_reduction2 = (unobscured - half_obscured2) / unobscured\n\n    \n    obscured = a_raw[:, 1958:3666].mean(axis=1)\n    unobscured = (a_raw[:, :1708].mean(axis=1) + a_raw[:, 3916:].mean(axis=1)) / 2\n    a_relative_reduction = (unobscured - obscured) / unobscured\n    \n    half_obscured1 = a_raw[:, 1708:1958].mean(axis=1)\n    half_obscured2 = a_raw[:, 3666:3916].mean(axis=1)\n    a_half_reduction1 = (unobscured - half_obscured1) / unobscured\n    a_half_reduction2 = (unobscured - half_obscured2) / unobscured\n\n\n    df = pd.DataFrame({'a_relative_reduction': a_relative_reduction,\n                       'f_relative_reduction': f_relative_reduction,\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    return df\n","metadata":{"tags":[]}},{"cell_type":"code","source":"\"\"\"\nobscured：取数组中第 23,500 到 44,000 列的平均值，代表行星遮挡星光时的亮度。\nunobscured：取数组中前 20,500 列和从 47,000 列之后的平均值，然后取平均，代表行星未遮挡星光时的亮度。\nhalf_obscured1,half_obscured2： 取数组中前 20,500 列到从23,500 列，44,000列到47,000 列的平均值，然后取平均，代表行星半遮挡星光时的亮度。\nf_relative_reduction：计算遮挡和未遮挡状态下亮度的相对减少量。\n\"\"\"\n\ndef feature_engineering(f_raw, a_raw):\n    \"\"\"Create a dataframe with two features from the raw data.\n    \n    Parameters:\n    f_raw: ndarray of shape (n_planets, 67500)\n    a_raw: ndarray of shape (n_planets, 5625)\n    \n    Return value:\n    df: DataFrame of shape (n_planets, 2)\n    \"\"\"\n    obscured = f_raw[:, 23500:44000].mean(axis=1)\n    unobscured = (f_raw[:, :20500].mean(axis=1) + f_raw[:, 47000:].mean(axis=1)) / 2\n    unobscured1 = f_raw[:, :20500].mean(axis=1)\n    unobscured2 = f_raw[:, 47000:].mean(axis=1)     \n    f_relative_reduction = (unobscured - obscured) / unobscured\n\n    half_obscured1 = f_raw[:, 20500:23500].mean(axis=1)\n    half_obscured2 = f_raw[:, 44000:47000].mean(axis=1)\n    f_relative_reduction_half1 = (unobscured - half_obscured1) / unobscured1\n    f_relative_reduction_half2 = (unobscured - half_obscured2) / unobscured2\n    \n    obscured = a_raw[:, 1958:3666].mean(axis=1)\n    unobscured = (a_raw[:, :1708].mean(axis=1) + a_raw[:, 3916:].mean(axis=1)) / 2\n    unobscured1 = a_raw[:, :1708].mean(axis=1)\n    unobscured2 = a_raw[:, 3916:].mean(axis=1)    \n    a_relative_reduction = (unobscured - obscured) / unobscured\n\n    half_obscured1 = f_raw[:, 1708:1958].mean(axis=1)\n    half_obscured2 = f_raw[:, 3666:3916].mean(axis=1)\n    a_relative_reduction_half1 = (unobscured - half_obscured1) / unobscured1\n    a_relative_reduction_half2 = (unobscured - half_obscured2) / unobscured2\n\n    df = pd.DataFrame({'a_relative_reduction': a_relative_reduction,\n                       'a_relative_reduction_half1':a_relative_reduction_half1,\n                       'a_relative_reduction_half2':a_relative_reduction_half2,\n                       'f_relative_reduction': f_relative_reduction,\n                       'f_relative_reduction_half1':f_relative_reduction_half1,\n                       'f_relative_reduction_half2':f_relative_reduction_half2,    \n                        \n                      })\n    \n    return df","metadata":{},"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":{"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train","metadata":{},"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.\n\n下面的散点图显示了当行星位于恒星前方时信号减少与我们想要预测的目标之一之间存在很强的相关性。其他目标也具有类似的高相关性。","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_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.\n\n为了保持简洁，我们使用岭回归来预测目标。\n\n我们关注三个交叉验证指标：\n\n1.R2 分数高于0.9，这证实了我们在散点图中看到的相关性。\n\n2.均方根误差将作为预测的不确定性。\n\n3.竞赛指标给出了排行榜分数的指示。不幸的是，竞赛指标取决于 sigma_true 的值，而我不知道这个值。","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_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":{"_kg_hide-input":true,"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_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_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":{"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#exec(open('competition_score.py', 'r').read())\n#exec(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_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":{"tags":[]},"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('../input/ariel-data-challenge-2024/test_adc_info.csv',\n                           index_col='planet_id')\nsample_submission = pd.read_csv('../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_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}