{"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":"gpu","dataSources":[{"sourceId":70367,"databundleVersionId":9188054,"sourceType":"competition"}],"dockerImageVersionId":30747,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Import The Libraries","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","metadata":{"execution":{"iopub.status.busy":"2024-08-20T09:33:36.685498Z","iopub.execute_input":"2024-08-20T09:33:36.68594Z","iopub.status.idle":"2024-08-20T09:33:40.26494Z","shell.execute_reply.started":"2024-08-20T09:33:36.685914Z","shell.execute_reply":"2024-08-20T09:33:40.263972Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from 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-20T09:33:40.266101Z","iopub.execute_input":"2024-08-20T09:33:40.266541Z","iopub.status.idle":"2024-08-20T09:33:40.271885Z","shell.execute_reply.started":"2024-08-20T09:33:40.266514Z","shell.execute_reply":"2024-08-20T09:33:40.270915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load & Read The Data","metadata":{}},{"cell_type":"code","source":"PATH = \"/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-20T09:33:43.054803Z","iopub.execute_input":"2024-08-20T09:33:43.055169Z","iopub.status.idle":"2024-08-20T09:33:43.194218Z","shell.execute_reply.started":"2024-08-20T09:33:43.055139Z","shell.execute_reply":"2024-08-20T09:33:43.193423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile process.py\n# Function to load signal data\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).to_numpy().flatten()\n    # Compute the mean of the signal and then the net signal (odd - even indices)\n    mean_signal = signal.reshape(-1, img_size).mean(axis=1)\n    net_signal = mean_signal[1::2] - mean_signal[0::2]\n    return net_signal\n\n# Function to read and preprocess signal data for all planet IDs\ndef read_and_preprocess(dataset, planet_ids, instrument=\"AIRS-CH0\"):\n    img_size = 1024 if instrument == \"FGS1\" else 32 * 356\n    signal_length = 67500 if instrument == 'FGS1' else 5625\n    raw_data = np.empty((len(planet_ids), signal_length), dtype=np.float32)\n    \n    # Load the signal for each planet and store it in the array\n    for i, planet_id in tqdm(enumerate(planet_ids), total=len(planet_ids)):\n        raw_data[i] = load_signal_data(planet_id, dataset, instrument, img_size)\n    \n    return raw_data\n\n# Function to perform feature engineering\ndef feature_engineering(f_raw, a_raw, adc_info):\n    # Calculate basic statistics like mean, std, variance, skewness, and kurtosis\n    def calculate_stats(data):\n        return {\n            'mean': data.mean(axis=1),\n            'std': data.std(axis=1),\n            'variance': data.var(axis=1),\n            'skewness': pd.DataFrame(data).skew(axis=1).values,\n            'kurtosis': pd.DataFrame(data).kurtosis(axis=1).values\n        }\n    \n    # Define obscured and unobscured ranges\n    def calculate_reduction(raw_data, obscured_range, unobscured_range):\n        obscured_mean = raw_data[:, obscured_range].mean(axis=1)\n        unobscured_mean = (raw_data[:, :unobscured_range].mean(axis=1) + raw_data[:, -unobscured_range:].mean(axis=1)) / 2\n        return (unobscured_mean - obscured_mean) / unobscured_mean, unobscured_mean / raw_data.std(axis=1)\n    \n    # Calculate the reductions and statistics for both f and a signals\n    f_reduction, f_signal_to_noise = calculate_reduction(f_raw, slice(23500, 44000), 20500)\n    a_reduction, a_signal_to_noise = calculate_reduction(a_raw, slice(1958, 3666), 1708)\n    \n    f_stats = calculate_stats(f_raw)\n    a_stats = calculate_stats(a_raw)\n    \n    # Combine all the features into a single DataFrame\n    features = pd.DataFrame({\n        'f_relative_reduction': f_reduction,\n        'f_signal_to_noise': f_signal_to_noise,\n        'f_variance': f_stats['variance'],\n        'f_skewness': f_stats['skewness'],\n        'f_kurtosis': f_stats['kurtosis'],\n        'a_relative_reduction': a_reduction,\n        'a_signal_to_noise': a_signal_to_noise,\n        'a_variance': a_stats['variance'],\n        'a_skewness': a_stats['skewness'],\n        'a_kurtosis': a_stats['kurtosis']\n    })\n\n    # Append adc_info to the features\n    features = pd.concat([features, adc_info.reset_index(drop=True).iloc[:, 1:6]], axis=1)\n    \n    return features\n","metadata":{"execution":{"iopub.status.busy":"2024-08-20T09:33:43.828011Z","iopub.execute_input":"2024-08-20T09:33:43.828384Z","iopub.status.idle":"2024-08-20T09:33:43.836932Z","shell.execute_reply.started":"2024-08-20T09:33:43.828355Z","shell.execute_reply":"2024-08-20T09:33:43.835902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile -a process.py\ndef postprocessing(pred_array, index, sigma_pred):\n    # Clip predictions to a minimum of 0 and create DataFrames for predictions and sigma\n    pred_df = pd.DataFrame(pred_array.clip(0, None), index=index, columns=wavelengths.columns)\n    sigma_df = pd.DataFrame(sigma_pred, index=index, columns=[f\"sigma_{i}\" for i in range(1, 284)])\n\n    # Concatenate predictions and sigma DataFrames\n    submission_df = pd.concat([pred_df, sigma_df], axis=1)\n    return submission_df","metadata":{"execution":{"iopub.status.busy":"2024-08-20T09:33:44.53483Z","iopub.execute_input":"2024-08-20T09:33:44.535197Z","iopub.status.idle":"2024-08-20T09:33:44.541354Z","shell.execute_reply.started":"2024-08-20T09:33:44.535169Z","shell.execute_reply":"2024-08-20T09:33:44.540358Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile -a process.py\nclass ParticipantVisibleError(Exception):\n    pass","metadata":{"execution":{"iopub.status.busy":"2024-08-20T09:33:45.387099Z","iopub.execute_input":"2024-08-20T09:33:45.387565Z","iopub.status.idle":"2024-08-20T09:33:45.394891Z","shell.execute_reply.started":"2024-08-20T09:33:45.387533Z","shell.execute_reply":"2024-08-20T09:33:45.39334Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%writefile -a process.py\ndef competition_score(solution, submission, naive_mean, naive_sigma, sigma_true, row_id_column_name='planet_id'):\n    # Remove the row ID column from the solution and submission DataFrames\n    solution = solution.drop(columns=[row_id_column_name])\n    submission = submission.drop(columns=[row_id_column_name])\n\n    # Check for invalid values and data types in the submission\n    if (submission < 0).any().any():\n        raise ParticipantVisibleError('Negative values in the submission')\n    if not all(pd.api.types.is_numeric_dtype(submission[col]) for col in submission.columns):\n        raise ParticipantVisibleError('All submission columns must be numeric')\n\n    # Validate the number of columns in the submission\n    n_wavelengths = solution.shape[1]\n    if submission.shape[1] != n_wavelengths * 2:\n        raise ParticipantVisibleError('Incorrect number of columns in the submission')\n\n    # Extract predicted values and sigma from the submission\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    # Calculate Gaussian 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    # Compute the submission score\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-20T09:33:46.468251Z","iopub.execute_input":"2024-08-20T09:33:46.46918Z","iopub.status.idle":"2024-08-20T09:33:46.475569Z","shell.execute_reply.started":"2024-08-20T09:33:46.469146Z","shell.execute_reply":"2024-08-20T09:33:46.474597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"exec(open('process.py', 'r').read())","metadata":{"execution":{"iopub.status.busy":"2024-08-20T09:33:49.366192Z","iopub.execute_input":"2024-08-20T09:33:49.366891Z","iopub.status.idle":"2024-08-20T09:33:49.372968Z","shell.execute_reply.started":"2024-08-20T09:33:49.366857Z","shell.execute_reply":"2024-08-20T09:33:49.371872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-20T09:33:50.243564Z","iopub.execute_input":"2024-08-20T09:33:50.243932Z","iopub.status.idle":"2024-08-20T09:45:28.06828Z","shell.execute_reply.started":"2024-08-20T09:33:50.243903Z","shell.execute_reply":"2024-08-20T09:45:28.067227Z"},"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-20T09:45:28.070718Z","iopub.execute_input":"2024-08-20T09:45:28.071285Z","iopub.status.idle":"2024-08-20T10:16:05.564069Z","shell.execute_reply.started":"2024-08-20T09:45:28.071256Z","shell.execute_reply":"2024-08-20T10:16:05.563094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ntrain = feature_engineering(f_raw_train, a_raw_train, train_adc_info)","metadata":{"execution":{"iopub.status.busy":"2024-08-20T10:16:39.938845Z","iopub.execute_input":"2024-08-20T10:16:39.939704Z","iopub.status.idle":"2024-08-20T10:16:42.454541Z","shell.execute_reply.started":"2024-08-20T10:16:39.939668Z","shell.execute_reply":"2024-08-20T10:16:42.453582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{"execution":{"iopub.status.busy":"2024-08-20T10:16:48.26338Z","iopub.execute_input":"2024-08-20T10:16:48.264259Z","iopub.status.idle":"2024-08-20T10:16:48.288145Z","shell.execute_reply.started":"2024-08-20T10:16:48.264226Z","shell.execute_reply":"2024-08-20T10:16:48.28718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = train.iloc[:,:-1]","metadata":{"execution":{"iopub.status.busy":"2024-08-20T10:16:54.763638Z","iopub.execute_input":"2024-08-20T10:16:54.76436Z","iopub.status.idle":"2024-08-20T10:16:54.769087Z","shell.execute_reply.started":"2024-08-20T10:16:54.76433Z","shell.execute_reply":"2024-08-20T10:16:54.767989Z"},"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-20T10:17:04.025132Z","iopub.execute_input":"2024-08-20T10:17:04.025883Z","iopub.status.idle":"2024-08-20T10:17:04.453966Z","shell.execute_reply.started":"2024-08-20T10:17:04.025848Z","shell.execute_reply":"2024-08-20T10:17:04.453098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Build & Train the model","metadata":{}},{"cell_type":"code","source":"from xgboost import XGBRegressor\n\n# Instantiate the XGBoost Regressor model\nmodel = XGBRegressor(\n    n_estimators=100,\n    learning_rate=0.1,\n    max_depth=6,\n    subsample=0.8,\n    random_state=42\n)\n\n# Generate out-of-fold predictions using cross-validation\noof_pred = cross_val_predict(model, train, train_labels, cv=5)\n\n# Calculate and print the R2 score\nprint(f\"# R2 score: {r2_score(train_labels, oof_pred):.4f}\")\n\n# Calculate and print the Root Mean Squared Error (RMSE)\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-20T10:19:55.395558Z","iopub.execute_input":"2024-08-20T10:19:55.396038Z","iopub.status.idle":"2024-08-20T10:20:41.600002Z","shell.execute_reply.started":"2024-08-20T10:19:55.395997Z","shell.execute_reply":"2024-08-20T10:20:41.598808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"col = 1\nplt.scatter(oof_pred[:,col], train_labels.iloc[:,col], s=15, c='lightgreen')\nplt.gca().set_aspect('equal')\nplt.xlabel('y_pred')\nplt.ylabel('y_true')\nplt.title('Comparing y_true and y_pred')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-08-20T10:21:21.078677Z","iopub.execute_input":"2024-08-20T10:21:21.079347Z","iopub.status.idle":"2024-08-20T10:21:21.35456Z","shell.execute_reply.started":"2024-08-20T10:21:21.079313Z","shell.execute_reply":"2024-08-20T10:21:21.353627Z"},"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-20T10:21:31.816251Z","iopub.execute_input":"2024-08-20T10:21:31.816646Z","iopub.status.idle":"2024-08-20T10:21:31.897828Z","shell.execute_reply.started":"2024-08-20T10:21:31.816614Z","shell.execute_reply":"2024-08-20T10:21:31.896862Z"},"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-20T10:21:49.77327Z","iopub.execute_input":"2024-08-20T10:21:49.773981Z","iopub.status.idle":"2024-08-20T10:21:59.615372Z","shell.execute_reply.started":"2024-08-20T10:21:49.773946Z","shell.execute_reply":"2024-08-20T10:21:59.614597Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Testing & Release the output","metadata":{}},{"cell_type":"code","source":"test_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, sigma_pred, 0.001), (1, 283)))\ndisplay(sub_df)\nsub_df.to_csv('submission.csv')","metadata":{"execution":{"iopub.status.busy":"2024-08-20T10:22:36.181125Z","iopub.execute_input":"2024-08-20T10:22:36.182021Z","iopub.status.idle":"2024-08-20T10:22:40.83906Z","shell.execute_reply.started":"2024-08-20T10:22:36.181977Z","shell.execute_reply":"2024-08-20T10:22:40.838324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}