{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":101849,"databundleVersionId":13093295,"sourceType":"competition"}],"dockerImageVersionId":31040,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ariel Data Challenge 2025: Introductory model: training\nCredita : @AmbrosM - This notebook is based on ADC24 Intro training ⭐️⭐️⭐️⭐️⭐️\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.","metadata":{}},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:34:26.91205Z","iopub.execute_input":"2025-06-27T04:34:26.912385Z","iopub.status.idle":"2025-06-27T04:34:29.812498Z","shell.execute_reply.started":"2025-06-27T04:34:26.91236Z","shell.execute_reply":"2025-06-27T04:34:29.811498Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# A look at the data\n\nWe start by reading the metadata:","metadata":{}},{"cell_type":"code","source":"import os\nimport pandas as pd\nfrom glob import glob\nfrom collections import defaultdict\nimport numpy as np\n\n# --- Load CSVs ---\ntrain_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train.csv')\nwavelengths_df = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/wavelengths.csv')\ntrain_star_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train_star_info.csv')\n\n# --- 1. Basic Train Set Stats ---\nprint(\"🪐 Number of training planets:\", train_df.shape[0])\nprint(\"📈 Number of target labels (wavelengths):\", train_df.shape[1] - 1)\nprint(\"🔬 Length of wavelength grid:\", wavelengths_df.shape[0])\n\n# --- 2. Target Stats (per flux column) ---\ntarget_cols = [col for col in train_df.columns if col != 'planet_id']\nflux_summary = train_df[target_cols].agg(['min', 'max', 'mean', 'std']).T\nprint(\"\\n📊 Flux value summary (first 5 rows):\")\nprint(flux_summary.head())\n\n# --- 3. Unique Stars ---\nif 'planet_id' in train_star_info.columns:\n    num_stars = train_star_info.drop(columns='planet_id').drop_duplicates().shape[0]\nelse:\n    num_stars = train_star_info.drop_duplicates().shape[0]\nprint(\"\\n🌟 Number of unique stars in training:\", num_stars)\n\n# --- 4. Planets with Multiple Observations ---\nobs_counts = defaultdict(int)\ntrain_planets = os.listdir('/kaggle/input/ariel-data-challenge-2025/train')\n\nfor pid in train_planets:\n    air_obs = glob(f\"/kaggle/input/ariel-data-challenge-2025/train/{pid}/AIRS-CH0_signal_*.parquet\")\n    obs_counts[pid] = len(air_obs)\n\nmulti_obs = {pid: count for pid, count in obs_counts.items() if count > 1}\nprint(\"\\n🔁 Planets with multiple observations:\", len(multi_obs))\n\n# --- 5. Check Calibration File Coverage ---\nmissing_calibs = []\nexpected = {\"dark\", \"dead\", \"flat\", \"linear_corr\", \"read\"}\n\nfor pid in train_planets:\n    for band in [\"AIRS-CH0\", \"FGS1\"]:\n        calib_path = f\"/kaggle/input/ariel-data-challenge-2025/train/{pid}/{band}_calibration_0\"\n        calib_files = {os.path.splitext(f)[0] for f in os.listdir(calib_path)} if os.path.exists(calib_path) else set()\n        missing = expected - calib_files\n        if missing:\n            missing_calibs.append((pid, band, missing))\n\nprint(\"\\n🧪 Planets missing calibration files:\", len(missing_calibs))\nif missing_calibs:\n    print(\"   Example:\", missing_calibs[0])\n\n# --- 6. Optional: Distribution of Observations Per Planet ---\nobs_distribution = pd.Series(list(obs_counts.values())).value_counts().sort_index()\nprint(\"\\n🗂 Observation count distribution per planet (AIR-CH0):\")\nprint(obs_distribution)\n\n# --- 7. Planet-Star Uniqueness Check ---\nmerged = pd.merge(train_df[['planet_id']], train_star_info, on='planet_id', how='left')\nunique_links = merged[['planet_id'] + [col for col in train_star_info.columns if col != 'planet_id']].drop_duplicates()\nprint(\"\\n🔗 Unique planet-star mappings:\", unique_links.shape[0])\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T06:40:04.402946Z","iopub.execute_input":"2025-07-20T06:40:04.403289Z","iopub.status.idle":"2025-07-20T06:40:22.707164Z","shell.execute_reply.started":"2025-07-20T06:40:04.403258Z","shell.execute_reply":"2025-07-20T06:40:22.706346Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Key Dataset Insights (Ariel 2025)\n\n\n- There are 1100 training planets, each corresponding to a unique star.\n\n- Each planet has 283 target labels, representing spectral flux values.\n\n- The wavelength grid has 1 shared configuration across all targets.\n\n- All 1100 stars are unique — one per planet.\n\n- Each planet is mapped to exactly one star (1100 unique planet-star pairs).\n\n- There are no repeated observations per planet.\n\n- A total of 2200 calibration file sets are missing (e.g., both AIRS-CH0 and FGS1 missing for each planet).\n\n- Example: Planet 1253730513 is missing all five calibration files in AIRS-CH0.\n\n- Flux values show low variability with means around 0.0146 and standard deviations around 0.0105 for the first 5 wavelengths.\n\n- Observation count per planet (AIR-CH0) is 0 for all planets — signal files are not present in the current training directory.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/adc_info.csv')\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-2025/train.csv',\n                           index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/wavelengths.csv')\naxis_info = pd.read_parquet('/kaggle/input/ariel-data-challenge-2025/axis_info.parquet')\ntrain_labels","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T01:39:29.664089Z","iopub.execute_input":"2025-07-20T01:39:29.665097Z","iopub.status.idle":"2025-07-20T01:39:30.482002Z","shell.execute_reply.started":"2025-07-20T01:39:29.665052Z","shell.execute_reply":"2025-07-20T01:39:30.481061Z"}},"outputs":[],"execution_count":null},{"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":{}},{"cell_type":"code","source":"import pandas as pd\nplanet_id = 1010375142\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.parquet')\nf_signal","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T06:43:50.746398Z","iopub.execute_input":"2025-07-20T06:43:50.74672Z","iopub.status.idle":"2025-07-20T06:43:52.411796Z","shell.execute_reply.started":"2025-07-20T06:43:50.746689Z","shell.execute_reply":"2025-07-20T06:43:52.410963Z"}},"outputs":[],"execution_count":null},{"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":{}},{"cell_type":"code","source":"import seaborn as sns\nimport matplotlib.pyplot as plt\nimport pandas as pd\n_, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))\nsns.heatmap(f_signal.iloc[0].values.reshape(32, 32), ax=ax1, vmin=0, vmax=52000)\nax1.set_aspect('equal')\nsns.heatmap(f_signal.iloc[1].values.reshape(32, 32), ax=ax2, vmin=0, vmax=52000)\nax2.set_aspect('equal')\nplt.suptitle('A pair of FGS1 images')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T06:45:00.490924Z","iopub.execute_input":"2025-07-20T06:45:00.49204Z","iopub.status.idle":"2025-07-20T06:45:02.204262Z","shell.execute_reply.started":"2025-07-20T06:45:00.491955Z","shell.execute_reply":"2025-07-20T06:45:02.203236Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n_,((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, figsize=(12, 4))\nplanet_id =4290810553\t\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.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 [19000, 25000, 50000, 56000]:\n    ax3.axvline(time_step, color='gray')\n\nplanet_id =1873185\nf_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/FGS1_signal_0.parquet')\n\nmean_signal = f_signal.values.mean(axis=1)\n\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.set_ylim(110, 120) \nax2.legend()\nax4.plot(smooth_signal, color='c', label='smoothened signal')\nax4.legend()\nax4.set_xlabel('time step')\nfor time_step in [7500, 15000, 54000, 60000]:\n    ax4.axvline(time_step, color='gray')\n\n# plt.suptitle('FGS1 time series', y=0.96)\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T07:12:53.564299Z","iopub.execute_input":"2025-07-20T07:12:53.564607Z","iopub.status.idle":"2025-07-20T07:12:56.184144Z","shell.execute_reply.started":"2025-07-20T07:12:53.564585Z","shell.execute_reply":"2025-07-20T07:12:56.183071Z"}},"outputs":[],"execution_count":null},{"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":{}},{"cell_type":"code","source":"planet_id = 1240764363\na_signal = pd.read_parquet(f'/kaggle/input/ariel-data-challenge-2025/train/{planet_id}/AIRS-CH0_signal_0.parquet')\na_signal","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:34:39.77942Z","iopub.execute_input":"2025-06-27T04:34:39.779686Z","iopub.status.idle":"2025-06-27T04:34:42.339396Z","shell.execute_reply.started":"2025-06-27T04:34:39.779663Z","shell.execute_reply":"2025-06-27T04:34:42.338529Z"}},"outputs":[],"execution_count":null},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:34:42.340352Z","iopub.execute_input":"2025-06-27T04:34:42.340638Z","iopub.status.idle":"2025-06-27T04:34:42.764319Z","shell.execute_reply.started":"2025-06-27T04:34:42.340616Z","shell.execute_reply":"2025-06-27T04:34:42.763137Z"}},"outputs":[],"execution_count":null},{"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":{}},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T07:16:31.8623Z","iopub.execute_input":"2025-07-20T07:16:31.862625Z","iopub.status.idle":"2025-07-20T07:16:34.866069Z","shell.execute_reply.started":"2025-07-20T07:16:31.862592Z","shell.execute_reply":"2025-07-20T07:16:34.865148Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Reading the data\n\nWe now read the FGS1 data and the AIRS-CH0 data for all 1100 training planets. 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":{}},{"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-2025/{dataset}/{planet_id}/FGS1_signal_0.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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:34:43.250035Z","iopub.execute_input":"2025-06-27T04:34:43.250308Z","iopub.status.idle":"2025-06-27T04:34:43.256698Z","shell.execute_reply.started":"2025-06-27T04:34:43.250288Z","shell.execute_reply":"2025-06-27T04:34:43.255921Z"}},"outputs":[],"execution_count":null},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:34:43.25771Z","iopub.execute_input":"2025-06-27T04:34:43.25807Z","iopub.status.idle":"2025-06-27T04:50:57.603221Z","shell.execute_reply.started":"2025-06-27T04:34:43.258039Z","shell.execute_reply":"2025-06-27T04:50:57.60072Z"}},"outputs":[],"execution_count":null},{"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-2025/{dataset}/{planet_id}/AIRS-CH0_signal_0.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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:50:57.606828Z","iopub.execute_input":"2025-06-27T04:50:57.607792Z","iopub.status.idle":"2025-06-27T04:50:57.622305Z","shell.execute_reply.started":"2025-06-27T04:50:57.607738Z","shell.execute_reply":"2025-06-27T04:50:57.621352Z"}},"outputs":[],"execution_count":null},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T04:50:57.6236Z","iopub.execute_input":"2025-06-27T04:50:57.623844Z","iopub.status.idle":"2025-06-27T05:07:56.692023Z","shell.execute_reply.started":"2025-06-27T04:50:57.623825Z","shell.execute_reply":"2025-06-27T05:07:56.690695Z"}},"outputs":[],"execution_count":null},{"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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:07:56.693947Z","iopub.execute_input":"2025-06-27T05:07:56.694325Z","iopub.status.idle":"2025-06-27T05:07:57.150169Z","shell.execute_reply.started":"2025-06-27T05:07:56.694301Z","shell.execute_reply":"2025-06-27T05:07:57.149187Z"}},"outputs":[],"execution_count":null},{"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).","metadata":{}},{"cell_type":"code","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    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    df = pd.DataFrame({'a_relative_reduction': a_relative_reduction,\n                       'f_relative_reduction': f_relative_reduction})\n    \n    return df\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:07:57.151278Z","iopub.execute_input":"2025-06-27T05:07:57.151616Z","iopub.status.idle":"2025-06-27T05:07:57.159115Z","shell.execute_reply.started":"2025-06-27T05:07:57.151582Z","shell.execute_reply":"2025-06-27T05:07:57.158072Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import pandas as pd\ntrain_star_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/train_star_info.csv')\ntrain_star_info=train_star_info.drop(columns='planet_id')\ntrain_star_info","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-21T05:21:19.884413Z","iopub.execute_input":"2025-07-21T05:21:19.884612Z","iopub.status.idle":"2025-07-21T05:21:22.108552Z","shell.execute_reply.started":"2025-07-21T05:21:19.884594Z","shell.execute_reply":"2025-07-21T05:21:22.107673Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"exec(open('feature_engineering.py', 'r').read())\n\ntrain = feature_engineering(f_raw_train, a_raw_train)\nprint(train)\nresult=pd.concat([train_star_info,train], axis=1)\ntrain=result","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-07-20T09:09:09.546975Z","iopub.execute_input":"2025-07-20T09:09:09.547341Z","iopub.status.idle":"2025-07-20T09:09:09.630658Z","shell.execute_reply.started":"2025-07-20T09:09:09.547305Z","shell.execute_reply":"2025-07-20T09:09:09.629494Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from sklearn.linear_model import Ridge\nfrom sklearn.model_selection import GridSearchCV, cross_val_score\nfrom sklearn.metrics import mean_squared_error\nimport numpy as np\ndef find_optimal_ridge_alpha(X, y, alphas=None, cv=5, scoring='neg_mean_squared_error', random_state=42):\n    \"\"\"\n    通过交叉验证找到最优的岭回归正则化系数(alpha)\n    \n    参数:\n    - X: 特征矩阵 (numpy数组或pandas DataFrame)\n    - y: 目标变量 (numpy数组或pandas Series)\n    - alphas: 要测试的alpha值列表，默认使用对数空间的100个值\n    - cv: 交叉验证折数，默认5折\n    - scoring: 评估指标，默认负均方误差\n    - random_state: 随机种子，确保结果可复现\n    \n    返回:\n    - optimal_alpha: 最优的alpha值\n    - best_score: 对应的模型性能得分(取绝对值后的均方误差)\n    - model: 训练好的最优Ridge模型\n    \"\"\"\n    # 设置alpha搜索范围\n    if alphas is None:\n        alphas = np.logspace(-3, 3, 100)  # 从0.001到1000的对数空间\n    \n    # 创建Ridge模型和网格搜索对象\n    ridge = Ridge(random_state=random_state)\n    grid_search = GridSearchCV(\n        estimator=ridge,\n        param_grid={'alpha': alphas},\n        cv=cv,\n        scoring=scoring,\n        n_jobs=-1  # 使用所有CPU核心\n    )\n    \n    # 执行网格搜索\n    grid_search.fit(X, y)\n    \n    # 获取最优参数和结果\n    optimal_alpha = grid_search.best_params_['alpha']\n    best_score = -grid_search.best_score_  # 转换为正的均方误差\n    best_model = grid_search.best_estimator_\n    \n    # 交叉验证确认结果\n    cv_scores = cross_val_score(\n        best_model, X, y, cv=cv, scoring='neg_mean_squared_error'\n    )\n    cv_mse = -cv_scores.mean()\n    \n    print(f\"最优alpha: {optimal_alpha:.4f}\")\n    #print(f\"交叉验证MSE: {cv_mse:.4f}\")\n    \n    return optimal_alpha#, cv_mse, best_model","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"aa=find_optimal_ridge_alpha(result,train_labels)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"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":{}},{"cell_type":"code","source":"model = Ridge(alpha=aa)\n\nimport numpy as np \nfrom sklearn.model_selection import KFold  \nfrom sklearn.base import clone \nmse = np.zeros(5)\npredictions = np.zeros((len(result),283))\nkf = KFold(n_splits=5, shuffle=True, random_state=42)\nfold=0\nfor train_idx, val_idx in kf.split(result):\n    X_train, X_val = result.iloc[train_idx], result.iloc[val_idx]\n    y_train,y_val = train_labels.iloc[train_idx] ,train_labels.iloc[val_idx]\n    model_clone = clone(model)\n    model_clone.fit(X_train, y_train)\n    fold_predictions = model_clone.predict(X_val)\n    predictions[val_idx] = fold_predictions\n    fold_mse = mean_squared_error(y_val, fold_predictions)\n    mse[fold] = fold_mse\n    fold=fold+1\n#oof_pred = cross_val_predict(model, train, train_labels)\noof_pred=predictions\nprint(mse)\nprint(f\"# R2 score: {r2_score(train_labels,predictions):.3f}\")\nsigma_pred = mean_squared_error(train_labels,predictions, squared=False)\nprint(f\"# Root mean squared error: {sigma_pred:.6f}\")\n# R2 score: 0.971\n# Root mean squared error: 0.000293","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:07:57.238178Z","iopub.execute_input":"2025-06-27T05:07:57.238497Z","iopub.status.idle":"2025-06-27T05:07:57.40667Z","shell.execute_reply.started":"2025-06-27T05:07:57.238472Z","shell.execute_reply":"2025-06-27T05:07:57.405768Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"col = 1\nplt.scatter(predictions[:,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":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:07:57.407712Z","iopub.execute_input":"2025-06-27T05:07:57.408228Z","iopub.status.idle":"2025-06-27T05:07:57.645887Z","shell.execute_reply.started":"2025-06-27T05:07:57.408196Z","shell.execute_reply":"2025-06-27T05:07:57.644799Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"%%writefile competition_score.py\n# Adapted from https://www.kaggle.com/code/metric/ariel-gaussian-log-likelihood\n\n# Custom error for invalid submissions\nclass ParticipantVisibleError(Exception):\n    pass\n\ndef competition_score(\n    solution: pd.DataFrame,\n    submission: pd.DataFrame,\n    naive_mean: float,\n    naive_sigma: float,\n    sigma_true: float,\n    row_id_column_name='planet_id'\n) -> float:\n    '''\n    Computes a Gaussian Log Likelihood-based score.\n    '''\n    # Drop ID columns\n    solution = solution.drop(columns=[row_id_column_name], errors='ignore')\n    submission = submission.drop(columns=[row_id_column_name], errors='ignore')\n\n    # Validation checks\n    if submission.min().min() < 0:\n        raise ParticipantVisibleError('Negative values in the submission')\n\n    for col in submission.columns:\n        if not pd.api.types.is_numeric_dtype(submission[col]):\n            raise ParticipantVisibleError(f'Submission column {col} must be numeric: {col}')\n\n    n_wavelengths = len(solution.columns)\n    if len(submission.columns) != 2 * n_wavelengths:\n        raise ParticipantVisibleError('Submission must have 2x columns of the solution')\n\n    # Extract predictions and sigmas\n    y_pred = submission.iloc[:, :n_wavelengths].values\n    sigma_pred = np.clip(submission.iloc[:, n_wavelengths:].values, a_min=1e-15, a_max=None)\n    y_true = solution.values\n\n    # Compute log likelihoods\n    GLL_pred = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_pred, scale=sigma_pred))\n    GLL_true = np.sum(scipy.stats.norm.logpdf(y_true, loc=y_true, scale=sigma_true))\n    GLL_mean = np.sum(scipy.stats.norm.logpdf(y_true, loc=naive_mean, scale=naive_sigma))\n\n    score = (GLL_pred - GLL_mean) / (GLL_true - GLL_mean)\n    return float(np.clip(score, 0.0, 1.0))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:07:57.647152Z","iopub.execute_input":"2025-06-27T05:07:57.647486Z","iopub.status.idle":"2025-06-27T05:07:57.654686Z","shell.execute_reply.started":"2025-06-27T05:07:57.647453Z","shell.execute_reply":"2025-06-27T05:07:57.653547Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n\ndef postprocessing(pred_array, index, sigma_pred, column_names=None):\n    \"\"\"\n    Creates a submission DataFrame with mean predictions and uncertainties.\n\n    Parameters:\n    - pred_array: ndarray of shape (n_samples, 283)\n    - index: pandas.Index of length n_samples\n    - sigma_pred: float or ndarray of shape (n_samples, 283)\n    - column_names: list of wavelength column names (optional)\n\n    Returns:\n    - df: DataFrame of shape (n_samples, 566)\n    \"\"\"\n    n_samples, n_waves = pred_array.shape\n\n    if column_names is None:\n        column_names = [f\"wl_{i+1}\" for i in range(n_waves)]\n\n    if np.isscalar(sigma_pred):\n        sigma_pred = np.full_like(pred_array, sigma_pred)\n\n    # Safety check\n    assert sigma_pred.shape == pred_array.shape, \"sigma_pred must match shape of pred_array\"\n    assert len(index) == n_samples, \"Index length must match number of rows\"\n\n    df_mean = pd.DataFrame(pred_array.clip(0, None), index=index, columns=column_names)\n    df_sigma = pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i+1}\" for i in range(n_waves)])\n\n    return pd.concat([df_mean, df_sigma], axis=1)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:08:17.650894Z","iopub.execute_input":"2025-06-27T05:08:17.651264Z","iopub.status.idle":"2025-06-27T05:08:17.658815Z","shell.execute_reply.started":"2025-06-27T05:08:17.65124Z","shell.execute_reply":"2025-06-27T05:08:17.657772Z"}},"outputs":[],"execution_count":null},{"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_labels.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.123","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:09:52.304393Z","iopub.execute_input":"2025-06-27T05:09:52.304794Z","iopub.status.idle":"2025-06-27T05:09:52.433541Z","shell.execute_reply.started":"2025-06-27T05:09:52.304768Z","shell.execute_reply":"2025-06-27T05:09:52.432619Z"}},"outputs":[],"execution_count":null},{"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)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:09:56.560123Z","iopub.execute_input":"2025-06-27T05:09:56.560495Z","iopub.status.idle":"2025-06-27T05:09:56.574776Z","shell.execute_reply.started":"2025-06-27T05:09:56.560468Z","shell.execute_reply":"2025-06-27T05:09:56.573755Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# Submission","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport pickle\n\n# 1. Load required files\ntest_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/test_star_info.csv', index_col='planet_id')\nsample_submission = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/sample_submission.csv', index_col='planet_id')\nwavelengths = pd.read_csv('/kaggle/input/ariel-data-challenge-2025/wavelengths.csv')\n\n# 2. Load model and sigma\nwith open('model.pickle', 'rb') as f:\n    model = pickle.load(f)\n\nwith open('sigma_pred.pickle', 'rb') as f:\n    sigma_pred = pickle.load(f)\n\n# 3. Run your preprocessing + feature extraction on test set\n# These must be implemented in your own code — you can adapt from training logic\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_features = feature_engineering(f_raw_test, a_raw_test)\ntest_features=pd.concat([train_star_info,test_features], axis=1)\n# 4. Predict\ntest_pred = model.predict(test_features)\n\n# 5. Postprocessing\ndef postprocessing(pred_array, index, sigma_pred, column_names):\n    \"\"\"\n    Convert predictions and uncertainty into final submission DataFrame.\n    \"\"\"\n    if np.isscalar(sigma_pred):\n        sigma_array = np.full_like(pred_array, sigma_pred)\n    else:\n        sigma_array = sigma_pred\n    df_pred = pd.DataFrame(pred_array.clip(0, None), index=index, columns=column_names)\n    df_sigma = pd.DataFrame(sigma_array, index=index, columns=[f\"sigma_{i}\" for i in range(1, len(column_names)+1)])\n    return pd.concat([df_pred, df_sigma], axis=1)\n\nsubmission_df = postprocessing(\n    pred_array=test_pred,\n    index=sample_submission.index,\n    sigma_pred=sigma_pred,\n    column_names=wavelengths.columns\n)\n\n# 6. Save\nsubmission_df.to_csv('submission.csv')\n\n# 7. Preview\n!head submission.csv\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-06-27T05:09:59.668684Z","iopub.execute_input":"2025-06-27T05:09:59.669037Z","iopub.status.idle":"2025-06-27T05:10:01.883461Z","shell.execute_reply.started":"2025-06-27T05:09:59.669008Z","shell.execute_reply":"2025-06-27T05:10:01.882091Z"}},"outputs":[],"execution_count":null}]}