notebooks/model_training_xgboost.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "<h1>MassBalanceMachine XGBoost Model Training - Iceland Region Example</h1>\n",
    "<p style='text-align: justify;'>\n",
    "In this notebook, we will simulate the glacier surface mass balance for the Iceland region using a custom <a href='https://xgboost.readthedocs.io/en/stable/'>XGBoost</a> model. The XGBoost model is designed with a custom objective function that generates monthly predictions based on aggregated observational data. We will create an instance of <code>CustomXGBoostRegressor</code> and train it using this custom loss function on the stake data from Iceland, which we have prepared in earlier notebooks. If you haven't already, please review the <a href='https://github.com/ODINN-SciML/MassBalanceMachine/blob/main/notebooks/data_preprocessing.ipynb'>data preprocessing</a> and <a href='https://github.com/ODINN-SciML/MassBalanceMachine/blob/main/notebooks/data_preprocessing.ipynb'>data processing WGMS</a> notebooks for more details.</p>\n",
    "\n",
    "<p style='text-align: justify;'>\n",
    "The workflow includes several key steps:\n",
    "</p>\n",
    "\n",
    "<ol style=\"margin-left: 20px; padding-left: 0;\">\n",
    "    <li style=\"margin-bottom: 10px;\">\n",
    "        <p style='text-align: justify;'><strong>Data Loading and Preparation:</strong> A <code>Dataloader</code> object is created to handle the loading of data and the creation of a training and testing split. This object also manages the generation of data splits for cross-validation.</p>\n",
    "    </li>\n",
    "    <li style=\"margin-bottom: 10px;\">\n",
    "        <p style='text-align: justify;'><strong>Cross-Validation and Model Training:</strong> Using Scikit-learn's cross-validation techniques, we explore different hyperparameters and train the model on the prepared data splits. This approach ensures a robust evaluation and helps in selecting suitable parameters.</p>\n",
    "    </li>\n",
    "    <li style=\"margin-bottom: 10px;\">\n",
    "        <p style='text-align: justify;'><strong>Aggregated Predictions:</strong> After training, we will display the aggregated monthly predictions generated by the model to visualize and analyze the results.</p>\n",
    "    </li>\n",
    "    <li style=\"margin-bottom: 10px;\">\n",
    "        <p style='text-align: justify;'><strong>Model Evaluation:</strong> Finally, the model's performance is evaluated on the test set, providing insights into its predictive accuracy for glacier mass balance.</p>\n",
    "    </li>\n",
    "</ol>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import massbalancemachine as mbm\n",
    "import warnings\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "warnings.filterwarnings('ignore')\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set a random seed:\n",
    "data = pd.read_csv('./example_data/iceland/files/iceland_monthly_dataset.csv')\n",
    "print('Number of winter and annual samples:', len(data))\n",
    "display(data)\n",
    "\n",
    "cfg = mbm.Config()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {},
   "source": [
    "<h2>1. Create the Train and Test Dataset and the Data Splits for Cross Validation</h2>\n",
    "<p style='text-align: justify;'>\n",
    "First, we create a <code>DataLoader</code> object, which generates both training and testing datasets, as well as the data splits required for cross-validation. To conserve memory, the <code>set_train_test_split</code> method returns iterators containing indices for the training and testing datasets. These indices are then used to retrieve the corresponding data for training and testing. Next, the <code>get_cv_split</code> method provides a list indicating the number of folds needed for cross-validation.\n",
    "</p>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a new DataLoader object with the monthly stake data measurements.\n",
    "dataloader = mbm.DataLoader(cfg, data=data)\n",
    "# Create a training and testing iterators. The parameters are optional. The default value of test_size is 0.3.\n",
    "train_itr, test_itr = dataloader.set_train_test_split(test_size=0.3)\n",
    "\n",
    "# Get all indices of the training and testing dataset at once from the iterators. Once called, the iterators are empty.\n",
    "train_indices, test_indices = list(train_itr), list(test_itr)\n",
    "\n",
    "# Get the features and targets of the training data for the indices as defined above, that will be used during the cross validation.\n",
    "df_X_train = data.iloc[train_indices]\n",
    "y_train = df_X_train['POINT_BALANCE'].values\n",
    "\n",
    "# Get test set\n",
    "df_X_test = data.iloc[test_indices]\n",
    "y_test = df_X_test['POINT_BALANCE'].values\n",
    "\n",
    "# Create the cross validation splits based on the training dataset. The default value for the number of splits is 5.\n",
    "splits = dataloader.get_cv_split(n_splits=5)\n",
    "\n",
    "# Print size of train and test\n",
    "print(f\"Size of training set: {len(train_indices)}\")\n",
    "print(f\"Size of test set: {len(test_indices)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {},
   "source": [
    "<h2>2. Create a CustomXGBoostRegressor Model</h2>\n",
    "<p style='text-align: justify;'>\n",
    "Next, we define the parameter ranges for each XGBoost parameter. In the subsequent step, we use cross-validation to explore these parameter ranges and select the combination that yields the lowest loss. Additionally, we create a <code>CustomXGBoostRegressor</code> object.\n",
    "</p>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each of the XGBoost parameter, define the grid range\n",
    "parameters = {\n",
    "    'max_depth': [\n",
    "        3,\n",
    "        4,\n",
    "        5,\n",
    "        6,\n",
    "    ],\n",
    "    'learning_rate': [0.01, 0.1, 0.2, 0.3],\n",
    "    'n_estimators': [100, 200, 300],\n",
    "    'gamma': [0, 1]\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a CustomXGBoostRegressor instance\n",
    "params_init = {\"device\": \"cpu\"}\n",
    "custom_xgboost = mbm.models.CustomXGBoostRegressor(cfg, **params_init)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8",
   "metadata": {},
   "source": [
    "<h2>3. Train the CustomXGBoostRegressor Model</h2>\n",
    "\n",
    "<p style='text-align: justify; margin-bottom: 5px;'>\n",
    "In the following cell, we begin training our model using either <strong>GridSearchCV</strong> or <strong>RandomizedSearchCV</strong>:\n",
    "</p>\n",
    "\n",
    "<ul style=\"margin-left: 20px; padding-left: 0; margin-bottom: 5px;\">\n",
    "  <li style=\"margin-bottom: 10px;\">\n",
    "    <p style='text-align: justify;'><strong>GridSearchCV</strong> performs an exhaustive search across all possible parameter combinations to find the best set for optimal performance using cross-validation. While this method is thorough, it is often time-consuming and computationally expensive.</p>\n",
    "  </li>\n",
    "  <li style=\"margin-bottom: 0px;\">\n",
    "    <p style='text-align: justify;'><strong>RandomizedSearchCV</strong>, on the other hand, samples a fixed number of parameter combinations from the distribution, making it more efficient in terms of time and computational resources, especially with larger hyperparameter spaces. However, this approach may miss some of the best parameter combinations that aren't selected in the random sampling.</p>\n",
    "  </li>\n",
    "</ul>\n",
    "\n",
    "<p style='text-align: justify;'>\n",
    "You can choose either of the two training methods. Both methods will use all CPU cores by default. If you want to adjust the number of cores used, you can change the <code>num_jobs</code> parameter.\n",
    "</p>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GridSearch\n",
    "# custom_xgboost.gridsearch(parameters=parameters, splits=splits, features=df_X_train, targets=y_train, num_jobs=-1)\n",
    "\n",
    "# RandomisedSearch, with n_iter the number of parameter settings that are sampled. Trade-off between goodness of the solution\n",
    "# versus runtime.\n",
    "custom_xgboost.randomsearch(\n",
    "    parameters=parameters,\n",
    "    n_iter=20,\n",
    "    splits=splits,\n",
    "    features=df_X_train,\n",
    "    targets=y_train,\n",
    ")\n",
    "best_params = params = custom_xgboost.param_search.best_params_\n",
    "best_estimator = custom_xgboost.param_search.best_estimator_\n",
    "print(\"Best parameters:\\n\", best_params)\n",
    "print(\"Best score:\\n\", custom_xgboost.param_search.best_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(custom_xgboost.param_search.cv_results_['mean_train_score'])\n",
    "plt.plot(custom_xgboost.param_search.cv_results_['mean_test_score'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "<h4>3.1 Save the Trained Model</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "# custom_xgboost.save_model('model.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# mbm.models.CustomXGBoostRegressor.load_model('model.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14",
   "metadata": {},
   "source": [
    "<h3>3.1 Show the Predictions</h3>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15",
   "metadata": {},
   "outputs": [],
   "source": [
    "def predVSTruth(grouped_ids, mae, rmse, title):\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n",
    "    legend_xgb = \"\\n\".join(\n",
    "        (r\"$\\mathrm{MAE_{xgb}}=%.3f, \\mathrm{RMSE_{xgb}}=%.3f$ \" % (\n",
    "            mae,\n",
    "            rmse,\n",
    "        ), ))\n",
    "\n",
    "    marker_xgb = 'o'\n",
    "    sns.scatterplot(grouped_ids,\n",
    "                    x=\"target\",\n",
    "                    y=\"pred\",\n",
    "                    ax=ax,\n",
    "                    alpha=0.5,\n",
    "                    marker=marker_xgb)\n",
    "\n",
    "    ax.set_ylabel('Predicted PMB [m w.e.]', fontsize=20)\n",
    "    ax.set_xlabel('Observed PMB [m w.e.]', fontsize=20)\n",
    "\n",
    "    ax.text(0.03,\n",
    "            0.98,\n",
    "            legend_xgb,\n",
    "            transform=ax.transAxes,\n",
    "            verticalalignment=\"top\",\n",
    "            fontsize=20)\n",
    "    ax.legend([], [], frameon=False)\n",
    "    # diagonal line\n",
    "    pt = (0, 0)\n",
    "    ax.axline(pt, slope=1, color=\"grey\", linestyle=\"-\", linewidth=0.2)\n",
    "    ax.axvline(0, color=\"grey\", linestyle=\"-\", linewidth=0.2)\n",
    "    ax.axhline(0, color=\"grey\", linestyle=\"-\", linewidth=0.2)\n",
    "    ax.grid()\n",
    "    ax.set_title(title, fontsize=20)\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set to CPU for predictions:\n",
    "xgb = best_estimator.set_params(device='cpu')\n",
    "\n",
    "# Make predictions on test\n",
    "features_test, metadata_test = xgb._create_features_metadata(df_X_test)\n",
    "y_pred = xgb.predict(features_test)\n",
    "\n",
    "# Make predictions aggr to meas ID:\n",
    "y_pred_agg = xgb.aggrPredict(metadata_test, features_test)\n",
    "\n",
    "# Calculate scores\n",
    "score = xgb.score(df_X_test, y_test)  # negative\n",
    "mse, rmse, mae, pearson = xgb.evalMetrics(metadata_test, y_pred, y_test)\n",
    "\n",
    "# Aggregate predictions to annual or winter:\n",
    "df_pred = df_X_test.copy()\n",
    "df_pred['target'] = y_test\n",
    "grouped_ids = df_pred.groupby('ID').agg({'target': 'mean'})\n",
    "grouped_ids['pred'] = y_pred_agg\n",
    "\n",
    "predVSTruth(grouped_ids, mae, rmse, title='XGBoost on test')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "MBM",
   "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}