{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"},{"sourceId":191155259,"sourceType":"kernelVersion"}],"dockerImageVersionId":30747,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Initialization\n\nThis competition seems requires strong scientific background and I had lot of confusion during EDA process. Therefore, I just build a simple starter for future coding.\n\n## Load Library","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","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:10:51.106998Z","iopub.execute_input":"2024-08-20T07:10:51.107698Z","iopub.status.idle":"2024-08-20T07:10:58.029091Z","shell.execute_reply.started":"2024-08-20T07:10:51.10766Z","shell.execute_reply":"2024-08-20T07:10:58.02778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load Meta-Data\nPATH = \"/kaggle/input/ariel-data-challenge-2024\"\ntrain_adc_info = pd.read_csv('/kaggle/input/ariel-data-challenge-2024/train_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(f'{PATH}/wavelengths.csv')","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:10:58.031333Z","iopub.execute_input":"2024-08-20T07:10:58.031987Z","iopub.status.idle":"2024-08-20T07:10:58.202278Z","shell.execute_reply.started":"2024-08-20T07:10:58.031945Z","shell.execute_reply":"2024-08-20T07:10:58.200999Z"},"trusted":true},"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 = \"/kaggle/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=15):\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":{"iopub.status.busy":"2024-08-20T07:10:58.209538Z","iopub.execute_input":"2024-08-20T07:10:58.209992Z","iopub.status.idle":"2024-08-20T07:10:58.222331Z","shell.execute_reply.started":"2024-08-20T07:10:58.209953Z","shell.execute_reply":"2024-08-20T07:10:58.221176Z"},"trusted":true},"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":{"iopub.status.busy":"2024-08-20T07:10:58.22363Z","iopub.execute_input":"2024-08-20T07:10:58.223947Z","iopub.status.idle":"2024-08-20T07:10:58.245472Z","shell.execute_reply.started":"2024-08-20T07:10:58.223919Z","shell.execute_reply":"2024-08-20T07:10:58.244004Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('utils.py', 'r').read())","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:10:58.247162Z","iopub.execute_input":"2024-08-20T07:10:58.247654Z","iopub.status.idle":"2024-08-20T07:10:58.264358Z","shell.execute_reply.started":"2024-08-20T07:10:58.247612Z","shell.execute_reply":"2024-08-20T07:10:58.263048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load Data","metadata":{}},{"cell_type":"code","source":"%%time\nif os.path.exists(\"/kaggle/input/adc24-intro-training/f_raw_train.pickle\"):\n    f_raw_train = np.load('/kaggle/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":{"iopub.status.busy":"2024-08-20T07:10:58.265909Z","iopub.execute_input":"2024-08-20T07:10:58.266368Z","iopub.status.idle":"2024-08-20T07:11:00.747942Z","shell.execute_reply.started":"2024-08-20T07:10:58.266326Z","shell.execute_reply":"2024-08-20T07:11:00.746569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nif os.path.exists(\"/kaggle/input/adc24-intro-training/a_raw_train.pickle\"):\n    a_raw_train = np.load('/kaggle/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":{"iopub.status.busy":"2024-08-20T07:11:00.749788Z","iopub.execute_input":"2024-08-20T07:11:00.750256Z","iopub.status.idle":"2024-08-20T07:11:00.911727Z","shell.execute_reply.started":"2024-08-20T07:11:00.750212Z","shell.execute_reply":"2024-08-20T07:11:00.910682Z"},"trusted":true},"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":{"iopub.status.busy":"2024-08-20T07:11:00.913192Z","iopub.execute_input":"2024-08-20T07:11:00.913539Z","iopub.status.idle":"2024-08-20T07:11:04.749914Z","shell.execute_reply.started":"2024-08-20T07:11:00.913508Z","shell.execute_reply":"2024-08-20T07:11:04.748788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:11:04.751095Z","iopub.execute_input":"2024-08-20T07:11:04.751429Z","iopub.status.idle":"2024-08-20T07:11:04.787574Z","shell.execute_reply.started":"2024-08-20T07:11:04.751402Z","shell.execute_reply":"2024-08-20T07:11:04.786433Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train.iloc[:,:-1]","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:11:04.790557Z","iopub.execute_input":"2024-08-20T07:11:04.790898Z","iopub.status.idle":"2024-08-20T07:11:04.81565Z","shell.execute_reply.started":"2024-08-20T07:11:04.790869Z","shell.execute_reply":"2024-08-20T07:11:04.814338Z"},"trusted":true},"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":{"iopub.status.busy":"2024-08-20T07:11:04.81734Z","iopub.execute_input":"2024-08-20T07:11:04.817682Z","iopub.status.idle":"2024-08-20T07:11:05.37462Z","shell.execute_reply.started":"2024-08-20T07:11:04.817643Z","shell.execute_reply":"2024-08-20T07:11:05.373434Z"},"trusted":true},"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":{"iopub.status.busy":"2024-08-20T07:11:05.37603Z","iopub.execute_input":"2024-08-20T07:11:05.376482Z","iopub.status.idle":"2024-08-20T07:11:05.721582Z","shell.execute_reply.started":"2024-08-20T07:11:05.376434Z","shell.execute_reply":"2024-08-20T07:11:05.720394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model\n## Rigde Model","metadata":{}},{"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):.4f}\")\nsigma_pred = mean_squared_error(train_labels, oof_pred, squared=False)\nprint(f\"# Root mean squared error: {sigma_pred:.7f}\")","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:11:05.722982Z","iopub.execute_input":"2024-08-20T07:11:05.723359Z","iopub.status.idle":"2024-08-20T07:11:12.235242Z","shell.execute_reply.started":"2024-08-20T07:11:05.723326Z","shell.execute_reply":"2024-08-20T07:11:12.233952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"oof_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:.4f}\")","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:11:12.236591Z","iopub.execute_input":"2024-08-20T07:11:12.237587Z","iopub.status.idle":"2024-08-20T07:11:12.361663Z","shell.execute_reply.started":"2024-08-20T07:11:12.237547Z","shell.execute_reply":"2024-08-20T07:11:12.360405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model.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":{"execution":{"iopub.status.busy":"2024-08-20T07:11:12.363107Z","iopub.execute_input":"2024-08-20T07:11:12.363517Z","iopub.status.idle":"2024-08-20T07:11:13.295769Z","shell.execute_reply.started":"2024-08-20T07:11:12.363486Z","shell.execute_reply":"2024-08-20T07:11:13.294528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Inference","metadata":{}},{"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 = read_and_preprocess('test', sample_submission.index, 'FGS1')\na_raw_test = read_and_preprocess('test', sample_submission.index)\ntest = feature_engineering(f_raw_test, a_raw_test, test_adc_info)\ntest = test.iloc[: , :-1]\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 = sub_df = postprocessing(test_pred,\n                        test_adc_info.index,\n                        sigma_pred=np.tile(np.where(test_adc_info[['star']] <= 1, 0.0001555, 0.00085), (1, 283)))\ndisplay(sub_df)\nsub_df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-08-20T07:11:13.29767Z","iopub.execute_input":"2024-08-20T07:11:13.298518Z","iopub.status.idle":"2024-08-20T07:11:16.817251Z","shell.execute_reply.started":"2024-08-20T07:11:13.298481Z","shell.execute_reply":"2024-08-20T07:11:16.815795Z"},"trusted":true},"execution_count":null,"outputs":[]}]}