{"metadata":{"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":191155259,"sourceType":"kernelVersion"}],"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":true},"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.9.19"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"I experimented with a variety of regression algorithms to compare their performance differences and find the model that best suits this competition. The results of the experiment are as follows:\n# Ridge and Linear Regression is the best\nI used MultiOutputRegressor to wrap a single regressor to make it suitable for multi-target regression tasks, but the efficiency is somewhat low. It takes about 7 hours to compute the results on my 13900K.\n\n## Ridge\n% R2 score: 0.9945\n\n% Root mean squared error: 0.0001241\n\n% Estimated competition score: 0.379\n## Linear Regression\n% R2 score: 0.995\n\n% Root mean squared error: 0.000124\n\n% Estimated competition score: 0.379\n## k-Nearest Neighbors\n% R2 score: 0.938\n\n% Root mean squared error: 0.000428\n\n% Estimated competition score: 0.203\n## Decision Tree\n% R2 score: 0.943\n\n% Root mean squared error: 0.000412\n\n% Estimated competition score: 0.209\n## Random Forest\n% R2 score: 0.976\n\n% Root mean squared error: 0.000267\n\n% Estimated competition score: 0.272\n## SGD\n% R2 score: -341251262544492403707458665757358839300096.000\n\n% Root mean squared error: 1007245294093768960.000000\n\n% Estimated competition score: 0.000\n## XGBoost\n% R2 score: 0.976\n\n% Root mean squared error: 0.000266\n\n% Estimated competition score: 0.273\n## AdaBoost\n% R2 score: 0.974\n\n% Root mean squared error: 0.000277\n\n% Estimated competition score: 0.266\n## ExtreTrees\n% R2 score: 0.980\n\n% Root mean squared error: 0.000242\n\n% Estimated competition score: 0.286\n\n","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport seaborn as sns\nimport pandas as pd\nimport polars as pl\nimport numpy as np\nimport torch\n# import torch.nn as nn\n# import torch.optim as optim\n# from torch.utils.data import DataLoader, Dataset\nfrom tqdm import tqdm\nimport pickle\nimport time\nimport os\nimport pickle\nimport seaborn as sns\nimport scipy.stats\nfrom sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\n\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load Meta-Data\nPATH = \"../input/ariel-data-challenge-2024\"\ntrain_adc_info = pd.read_csv('../input/ariel-data-challenge-2024/train_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(f'{PATH}/wavelengths.csv')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pre-Processing\n## Load Functions","metadata":{}},{"cell_type":"code","source":"%%writefile utils.py\nimport pandas as pd\nimport polars as pl\nimport numpy as np\nfrom tqdm import tqdm\nimport pickle\nPATH = \"../input/ariel-data-challenge-2024\"\ndef load_signal_data(planet_id, dataset, instrument, img_size):\n    file_path = f'{PATH}/{dataset}/{planet_id}/{instrument}_signal.parquet'\n    signal = pl.read_parquet(file_path)\n    mean_signal = signal.cast(pl.Int32).sum_horizontal().cast(pl.Float32).to_numpy() / img_size # mean over the 32*32 pixels\n    net_signal = mean_signal[1::2] - mean_signal[0::2]\n    return net_signal\n\ndef read_and_preprocess(dataset, planet_ids, instrument = \"AIRS-CH0\"):\n    \"\"\"Read the files for all planet_ids and extract the time series.\n    Parameters\n    dataset: 'train' or 'test'\n    planet_ids: list of planet ids\n    instrument: the instrument of observation, 'AIRS-CH0' or 'FGS1', default to 'AIRS-CH0'\n    Returns\n    dataframe with one row per planet_id and 67500 values per row for FGS1 and 5624 for AIRS-CH0\n    \"\"\"\n    img_size = 1024 if instrument == \"FGS1\" else 32*356\n    column_num = 67500 if instrument == 'FGS1' else 5625\n    raw_train = np.full((len(planet_ids), column_num), np.nan, dtype=np.float32)\n    for i, planet_id in tqdm(list(enumerate(planet_ids))):\n        raw_train[i] = load_signal_data(planet_id, dataset, instrument, img_size)\n    return raw_train\n\ndef feature_engineering(f_raw, a_raw, adc_info, 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    df = pd.concat([df, adc_info.reset_index().iloc[:, 1:6]], axis=1)\n    \n    return df","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile -a utils.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\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":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('utils.py', 'r').read())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load Data","metadata":{}},{"cell_type":"code","source":"%%time\nif os.path.exists(\"../input/adc24-intro-training/f_raw_train.pickle\"):\n    f_raw_train = np.load('../input/adc24-intro-training/f_raw_train.pickle', allow_pickle=True)\nelse:\n    f_raw_train = read_and_preprocess('train', train_labels.index, 'FGS1')\n    with open('f_raw_train.pickle', 'wb') as f:\n        pickle.dump(f_raw_train, f)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif os.path.exists(\"../input/adc24-intro-training/a_raw_train.pickle\"):\n    a_raw_train = np.load('../input/adc24-intro-training/a_raw_train.pickle', allow_pickle=True)\nelse:\n    a_raw_train = read_and_preprocess('train', train_labels.index)\n    with open('a_raw_train.pickle', 'wb') as f:\n        pickle.dump(a_raw_train, f)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Feature Engineering","metadata":{}},{"cell_type":"code","source":"%%time\ntrain = feature_engineering(f_raw_train, a_raw_train, train_adc_info)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train.iloc[:,:-1]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Data Plot","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_count":null,"outputs":[]},{"cell_type":"code","source":"color_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')\n# plt.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":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import cross_val_predict\nfrom sklearn.linear_model import Ridge\nfrom sklearn.metrics import r2_score, mean_squared_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.model_selection import KFold\n\nfrom sklearn.linear_model import LinearRegression, SGDRegressor, Ridge\nfrom sklearn.neighbors import KNeighborsRegressor\nfrom sklearn.tree import DecisionTreeRegressor\nfrom sklearn.svm import LinearSVR\nfrom sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, ExtraTreesRegressor\nfrom xgboost import XGBRegressor\n\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import cross_val_score, RepeatedKFold\nfrom sklearn.multioutput import MultiOutputRegressor, RegressorChain","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Candidate models\nmodels = {\n    'Linear Regression'  : LinearRegression(),\n    'k-Nearest Neighbors': KNeighborsRegressor(),\n    'Decision Tree'      : DecisionTreeRegressor(random_state = 42),\n    #'Linear SVM'         : LinearSVR(max_iter = 1000),\n    'Random Forest'      : RandomForestRegressor(random_state = 42),\n    'SGD'                : SGDRegressor(random_state = 42),\n    #'Ridge'              : Ridge(alpha=1e-12),\n    'XGBoost'            : XGBRegressor(random_state = 42),\n    'AdaBoost'           : AdaBoostRegressor(random_state = 42),\n    'ExtreTrees'         : ExtraTreesRegressor(random_state = 42)\n}","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# ","metadata":{}},{"cell_type":"code","source":"print('Ridge')\nmodel = Ridge(alpha=1e-12)\noof_pred = cross_val_predict(model, train, train_labels)\nprint(f\"# R2 score: {r2_score(train_labels, oof_pred):.4f}\")\nsigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\nprint(f\"# Root mean squared error: {sigma_pred:.7f}\")\noof_df = postprocessing(oof_pred, train_adc_info.index, sigma_pred)\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# Performance of candidate models\nfor name, model in zip(models.keys(), models.values()):\n    print(name)\n    regressor = MultiOutputRegressor(model, n_jobs = -1)\n    oof_pred = cross_val_predict(regressor, train, train_labels)\n    print(f\"# R2 score: {r2_score(train_labels, oof_pred):.3f}\")\n    sigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\n    print(f\"# Root mean squared error: {sigma_pred:.6f}\")\n    oof_df = postprocessing(oof_pred, train_adc_info.index, sigma_pred)\n    gll_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)\n    print(f\"# Estimated competition score: {gll_score:.3f}\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}