Notebook_7_run_catboost_lightgbm_xgboost_regressions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce6051bb",
   "metadata": {},
   "source": [
    "# Creating Boosted Tree regression models from GHG emission model input-output data\n",
    "\n",
    "### The analysis is performed via introspection of the fitted CatBoost, XGBoost and LightGBM boosted random-forest regression models using various explainable AI techniques in DALEX package and the `shap` package\n",
    "\n",
    "- Author: Tomasz Janus\n",
    "- E-mail: tomasz.k.janus@gmail.com, tomasz.janus@manchester.ac.uk\n",
    "- Mui Ne, 22/10/2023\n",
    "- Modified on 28/06/2024"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbd11ac6",
   "metadata": {},
   "source": [
    "### The notebooks proceeds in the following steps:\n",
    "  1. Load required libraries and the input/output data\n",
    "  2. Visualise relationships in the input data\n",
    "  3. Prune / clean the input dataset\n",
    "  4. Fit the catboost model using pre-set hyperparameter values and a fixed train/validation data-split to serve as a baseline quick check of what we can expect from boosted tree models\n",
    "  5. Fit the catboost, lightgbm and xgboost models using hyperparameter tuning and KFOLD cross-validation\n",
    "  6. Save the fitted models to files\n",
    "  7. Explore the model structure using DALEX (that provides interface to SHAP and LIME) and additionally using the `shap` package\n",
    "  \n",
    "## NOTE:\n",
    "#### Creating composite explanation figures requires prior calculation of emission intensity predictions that are calculated in `Notebook_9b_process_additional_information_from_water_resource_models`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6bb10f0",
   "metadata": {},
   "source": [
    "## 1a. Import the required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fce804f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Protocol, Dict, Protocol, List, Any, Literal, Tuple\n",
    "import os\n",
    "import sys\n",
    "import pathlib\n",
    "import re\n",
    "from dataclasses import dataclass\n",
    "from collections import defaultdict\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.pylab as plt2\n",
    "import seaborn as sns # For plotting data\n",
    "\n",
    "# Load tree-based regression models\n",
    "import catboost as cb\n",
    "import xgboost as xgb\n",
    "import lightgbm as lgbm\n",
    "from catboost import CatBoostRegressor\n",
    "from xgboost import XGBRegressor\n",
    "from lightgbm import LGBMRegressor\n",
    "\n",
    "# Load scikit-learn's classes for model and feature selection, validation and data transformation\n",
    "from sklearn.model_selection import train_test_split, cross_val_score, StratifiedKFold, KFold\n",
    "from sklearn.inspection import permutation_importance\n",
    "from sklearn.feature_selection import VarianceThreshold # Feature selecto\n",
    "from sklearn.pipeline import Pipeline, make_pipeline\n",
    "#from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score\n",
    "from sklearn.compose import make_column_transformer, ColumnTransformer\n",
    "from sklearn.feature_selection import SelectKBest, f_regression, mutual_info_regression\n",
    "\n",
    "# Enable the output data from scikit-learn's Pipeline to be in Pandas, rather than numpy ndarray format\n",
    "from sklearn import set_config\n",
    "set_config(transform_output=\"pandas\")\n",
    "\n",
    "# Import the hyperparameter optimal tuning tool\n",
    "import lib.hypertune as hypertune\n",
    "\n",
    "# Import library for automatic data profiling\n",
    "# Although, it seems to be running very slowly on our input data hence we might resort to manual \n",
    "# data exploration\n",
    "from ydata_profiling import ProfileReport\n",
    "\n",
    "# Import explainers\n",
    "import shap\n",
    "from lime import lime_tabular\n",
    "import dalex as dx\n",
    "\n",
    "# Model loading/saving\n",
    "import joblib\n",
    "\n",
    "# Import from local library folder\n",
    "from lib.pipeline import DataFrameOneHotEncoder\n",
    "from lib.pipeline import (\n",
    "    ProfileDropper, ReduceResAreaFractions, ReplaceTempProfileWithMean,\n",
    "    ColumnDropper, StaticColRemover)\n",
    "\n",
    "import importlib\n",
    "from lib.hypertune import HyperTuner, hypertune_model\n",
    "from lib.utils import (\n",
    "    save_model, load_model, plot_gini_feature_importances, plot_permutation_feature_importances,\n",
    "    plot_shap_feature_importances, model_check, plot_scores)\n",
    "from lib.utils import (\n",
    "    calculate_gini_feature_importances, calculate_permutation_feature_importances,\n",
    "    calculate_shap_feature_importances)\n",
    "from lib.utils import model_feature_importances\n",
    "import pickle\n",
    "from math import isclose\n",
    "\n",
    "import svgutils.transform as sg \n",
    "\n",
    "plt.style.use('ggplot')\n",
    "pd.set_option(\"display.max_columns\", 200)\n",
    "plt.rcParams['figure.figsize'] = (6,6)\n",
    "\n",
    "# -----------------------------------------\n",
    "# Internal functions\n",
    "# -----------------------------------------\n",
    "def filter_corr_matrix(corr_df: pd.DataFrame, threshold: float = 0.80) -> pd.DataFrame:\n",
    "    \"\"\"Removes rows and columns that do not feature any correlation coefficient larger equal threshold\"\"\"\n",
    "    selected_columns = []\n",
    "    for col in corr_df.columns:\n",
    "        series = corr_df[col].drop(col, axis=0)\n",
    "        if (series.abs() >= threshold).any():\n",
    "            selected_columns.append(col)\n",
    "\n",
    "    return corr_df.loc[selected_columns,selected_columns]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10d27829",
   "metadata": {},
   "source": [
    "### Execution options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e41dc749",
   "metadata": {},
   "outputs": [],
   "source": [
    "# EXECUTION OPTIONS\n",
    "fitted_models = set(['xgboost', 'catboost', 'lightgbm'])\n",
    "prune_additional_features: bool = True # Remove some of the potentially redundant features\n",
    "\n",
    "rerun = False # Runs model fitting even if the model in saved form already exists\n",
    "override = False # Saves the model (after fitting) even if saved model already exists\n",
    "recalculate_breakdown_interactions = False"
   ]
  },
  {
   "cell_type": "raw",
   "id": "704e24ad",
   "metadata": {},
   "source": [
    "TODO: Check autofeat package"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "268a7a32",
   "metadata": {},
   "source": [
    "## 1b. Input data load, check and merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfc90980",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "data_full = pd.read_csv(\n",
    "    pathlib.Path(\"outputs/reemission/combined/combined_outputs.csv\"))\n",
    "elev_data = pd.read_csv(\n",
    "    pathlib.Path(\"config/elev.csv\"))\n",
    "# Add elevation data to the full emission estimates dataset\n",
    "merged_df = pd.merge(\n",
    "    data_full, elev_data[['ifc_name', 'fsl_masl']], left_on='Name', right_on='ifc_name', how='inner')\n",
    "merged_df.drop(columns=['ifc_name'], inplace=True)\n",
    "# Perform prefiltering of data outside the pipeline - required for sorting data with regards to:\n",
    "#    soil types, landuse intensities and treatment factors (levels of ww treatment in the catchment)\n",
    "merged_df[['Soil', 'Landuse intensity', 'Treatment']] = \\\n",
    "            merged_df['Scenario'].str.split('_', expand=True)\n",
    "# Create a version with CH4 emissions without degassing - required to test a 'what-if' scenario if\n",
    "# reservoirs are operated with shallow water intakes\n",
    "merged_df['ch4_net_nodegassing'] = merged_df['ch4_net'] - merged_df['ch4_degassing']\n",
    "# Make sure that we haven't removed any rows during merging\n",
    "assert len(merged_df) == len(data_full)\n",
    "# Create a separate dataframe for plotting purposes\n",
    "df_plot = merged_df.copy()\n",
    "df_plot['Name_Code'] = pd.factorize(df_plot['Name'])[0]\n",
    "df_plot['Scenario_Code'] = pd.factorize(df_plot['Scenario'])[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e396ea85",
   "metadata": {},
   "source": [
    "## 1c. Filter data using pipelines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f3c08b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify columns to be dropped from data\n",
    "corr_vars = ['trophic_status'] # It's probably correlated with other data as it's calculated from \n",
    "                              # other data within reemission\n",
    "aux_vars = ['id', 'type', 'gasses_0', 'gasses_1', 'gasses_2']\n",
    "interm_result_vars = [\n",
    "    'co2_diffusion', 'co2_diffusion_nonanthro', 'co2_preimp', 'co2_minus_nonanthro', 'co2_net', \n",
    "    'co2_total_lifetime', 'ch4_diffusion', 'ch4_ebullition', 'ch4_degassing', 'ch4_preimp', 'ch4_net',\n",
    "    'ch4_total_lifetime', 'n2o_methodA', 'n2o_methodB', 'n2o_mean', 'n2o_total_lifetime',\n",
    "    'co2_total_per_year', 'ch4_total_per_year', 'n2o_total_per_year', 'nitrogen_downstream_conc']\n",
    "marginal_vars = [\n",
    "    'catch_riv_length', 'res_water_intake_depth', \n",
    "    'surface_density', 'bottom_density'] \n",
    "# comment: res water intake depth does not play part in regression (deep intake only)\n",
    "duplicated_vars = ['mean_radiance_lat'] # duplicated with res_mean_radiance\n",
    "additional_vars_prune = ['coordinates_0', 'coordinates_1', 'bottom_temperature', \n",
    "                         'surface_temperature', 'fsl_masl', 'thermocline_depth',\n",
    "                         'nitrogen_load', 'phosphorus_load', 'res_max_depth',\n",
    "                         'littoral_area_frac', 'retention_coeff']\n",
    "\n",
    "columns_to_drop = corr_vars + aux_vars + interm_result_vars + marginal_vars + duplicated_vars + \\\n",
    "    additional_vars_prune\n",
    "# Add additional columns that should not play part in the model as they're either internally calculated\n",
    "# or highly correlated with other variables\n",
    "if prune_additional_features:\n",
    "    additional_columns_to_drop = [\n",
    "            'reservoir_tn', 'reservoir_tp', 'inflow_p_conc', 'inflow_n_conc',\n",
    "            'global_radiance', 'catch_precip', 'ch4_net_nodegassing']\n",
    "else:\n",
    "    additional_columns_to_drop = []\n",
    "# reservoir_tn and reservoir_tp are internal variables, inflow_p_conc and inflow_n_conc are correlated with \n",
    "# nitrogen and phosphorus loads, global radiance is correlated with res_mean_radiance, catch_precip is highly\n",
    "# correlated with catch_runoff\n",
    "columns_to_drop.extend(additional_columns_to_drop)\n",
    "\n",
    "if 'bottom_temperature' in columns_to_drop or 'surface_temperature' in columns_to_drop:\n",
    "    plot_pairplot1 = False\n",
    "if 'population' in columns_to_drop or 'phosphorus_load' in columns_to_drop:\n",
    "    plot_pairplot2 = False\n",
    "plot_pairplot3 = True\n",
    "if 'res_max_depth' in columns_to_drop or 'res_mean_depth' in columns_to_drop:\n",
    "    plot_pairplot4 = False\n",
    "plot_pairplot5 = True\n",
    "    \n",
    "print(f\"Columns selected from dropping from the dataset: \\n{' ** '.join(columns_to_drop)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e36e9c78",
   "metadata": {},
   "source": [
    "### Create a preprocessing pipeline\n",
    "\n",
    "#### Think of other columns to drop, e.g. `nitrogen_load`, `phosphorus_load`\n",
    "#### Perhaps think of adding other preprocessing steps such as removing low variance columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08a1bb9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "pipe_1_options = [\n",
    "    ('prof_dropper', ProfileDropper()), # drop emission profile outputs from data\n",
    "    ('red_res_area_fractions', ReduceResAreaFractions()),\n",
    "    ('mean_temp', ReplaceTempProfileWithMean()),\n",
    "    ('col_drop_1', ColumnDropper(columns_to_drop + ['Scenario', 'Name'])),\n",
    "    ('stat_col_rem', VarianceThreshold())\n",
    "    #('stat_col_rem', StaticColRemover())\n",
    "]\n",
    "pipe_1 = Pipeline(pipe_1_options)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27500e13",
   "metadata": {},
   "source": [
    "### Create different data options for different types of analyses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b37e10",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Divide raw data into mineral soils and organic soils\n",
    "merged_df_min = merged_df[merged_df['Soil']=='MIN'].drop(\"Soil\", axis=1)\n",
    "merged_df_org = merged_df[merged_df['Soil']=='ORG'].drop(\"Soil\", axis=1)\n",
    "# Drop other options, i.e treatment and landuse intensity\n",
    "merged_df_min_prim_low = merged_df_min.query(\n",
    "    \"`Landuse intensity` == 'LOW' & Treatment == 'PRIM'\").drop([\"Landuse intensity\", \"Treatment\"], axis=1)\n",
    "merged_df_org_prim_low = merged_df_org.query(\n",
    "    \"`Landuse intensity` == 'LOW' & Treatment == 'PRIM'\").drop([\"Landuse intensity\", \"Treatment\"], axis=1)\n",
    "\n",
    "# Create datasets for CO2 regression and CH4 regression tasks for scenario with mineral soil,\n",
    "# primary treatment and low landuse intensity\n",
    "X_co2_min = merged_df_min_prim_low.drop(columns=['co2_net'])\n",
    "y_co2_min = merged_df_min_prim_low['co2_net']\n",
    "X_ch4_min = merged_df_min_prim_low.drop(columns='ch4_net')\n",
    "y_ch4_min = merged_df_min_prim_low['ch4_net']\n",
    "# Same for organic soil, primary treatment and low landuse intensity\n",
    "X_co2_org = merged_df_org_prim_low.drop(columns=['co2_net'])\n",
    "y_co2_org = merged_df_org_prim_low['co2_net']\n",
    "X_ch4_org = merged_df_org_prim_low.drop(columns='ch4_net')\n",
    "y_ch4_org = merged_df_org_prim_low['ch4_net']\n",
    "\n",
    "# Perform data splitting - use 90% train and 10% test\n",
    "_co2_data_random_seed = 42\n",
    "_ch4_data_random_seed = 42\n",
    "\n",
    "X_co2_train, X_co2_test, y_co2_train, y_co2_test = \\\n",
    "    train_test_split(X_co2_min, y_co2_min, train_size=0.9, test_size=0.1, random_state=_co2_data_random_seed)\n",
    "X_ch4_train, X_ch4_test, y_ch4_train, y_ch4_test = \\\n",
    "    train_test_split(X_ch4_min, y_ch4_min, train_size=0.9, test_size=0.1, random_state=_ch4_data_random_seed)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f98e0144",
   "metadata": {},
   "source": [
    "### Apply pipelines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f092ab7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fit_pipeline = pipe_1.fit(X_co2_train, y_co2_train)\n",
    "X_co2_train = fit_pipeline.transform(X_co2_train)\n",
    "X_co2_test = fit_pipeline.transform(X_co2_test)\n",
    "X_ch4_train = pipe_1.fit_transform(X_ch4_train, y_ch4_train)\n",
    "X_ch4_test = pipe_1.transform(X_ch4_test)\n",
    "# Merge test and train data - for model introspection we would like to look into the model trained\n",
    "# on all data. Train/test split is done to check how the model trained on train data generalizes to\n",
    "# test data in order to sense if it's under or overfitting\n",
    "X_co2_train_test = pd.concat([X_co2_train, X_co2_test])\n",
    "y_co2_train_test = pd.concat([y_co2_train, y_co2_test])\n",
    "X_ch4_train_test = pd.concat([X_ch4_train, X_ch4_test])\n",
    "y_ch4_train_test = pd.concat([y_ch4_train, y_ch4_test])\n",
    "### Rename some columns to improve the understanding of variables\n",
    "col_name_map = {\n",
    "    \"catch_area_fractions_0\": \"catchment bare soil fraction\",\n",
    "    \"catch_area_fractions_1\": \"catchment snow and ice fraction\",\n",
    "    \"catch_area_fractions_2\": \"catchment urban area fraction\",\n",
    "    \"catch_area_fractions_3\": \"catchment water area fraction\",\n",
    "    \"catch_area_fractions_4\": \"catchment wetland area fraction\",\n",
    "    \"catch_area_fractions_5\": \"catchment crop area fraction\",\n",
    "    \"catch_area_fractions_6\": \"catchment shrub area fraction\",\n",
    "    \"catch_area_fractions_7\": \"catchment forest area fraction\",\n",
    "    \"catch_area_fractions_8\": \"catchment unknown area fraction\",\n",
    "    \"res_area_fractions_red_0\": \"reservoir bare soil fraction\",\n",
    "    \"res_area_fractions_red_1\": \"reservoir snow and ice fraction\",\n",
    "    \"res_area_fractions_red_2\": \"reservoir urban area fraction\",\n",
    "    \"res_area_fractions_red_3\": \"reservoir water area fraction\",\n",
    "    \"res_area_fractions_red_4\": \"reservoir wetland area fraction\",\n",
    "    \"res_area_fractions_red_5\": \"reservoir crop area fraction\",\n",
    "    \"res_area_fractions_red_6\": \"reservoir shrub area fraction\",\n",
    "    \"res_area_fractions_red_7\": \"reservoir forest area fraction\",\n",
    "    \"res_area_fractions_red_8\": \"reservoir unknown area fraction\",\n",
    "    \"coordinates_0\": \"latitude\",\n",
    "    \"coordinates_1\": \"longitude\",\n",
    "    'catch_runoff': 'catchment runoff',\n",
    "    'catch_area': 'catchment area',\n",
    "    'catch_population': 'population',\n",
    "    'catch_slope': 'catchment slope',\n",
    "    'catch_etransp': 'evapotranspiration',\n",
    "    'catch_soil_wetness': 'catchment soil wetness',\n",
    "    'catch_mean_olsen': 'catchment mean olsen',\n",
    "    'res_volume': 'reservoir volume',\n",
    "    'res_area': 'reservoir area',\n",
    "    'res_max_depth': 'max depth',\n",
    "    'res_mean_depth': 'mean depth',\n",
    "    'res_soil_carbon': 'reservoir soil carbon',\n",
    "    'res_mean_radiance': 'reservoir mean radiance',\n",
    "    'res_mean_radiance_may_sept': 'reservoir mean radiance may-sept',\n",
    "    'res_mean_radiance_nov_mar': 'reservoir mean radiance nov-mar',\n",
    "    'res_mean_monthly_windspeed': 'reservoir mean monthly windspeed',\n",
    "    'retention_coeff': 'retention coefficient',\n",
    "    'littoral_area_frac': 'littoral area fraction',\n",
    "    'bottom_temperature': 'bottom temp',\n",
    "    'surface_temperature': 'surface temp',\n",
    "    'thermocline_depth': 'thermocline depth',\n",
    "    'nitrogen_load': 'N load',\n",
    "    'phosphorus_load': 'P load',\n",
    "    'fsl_masl': 'fsl',\n",
    "    'ave_temp': 'air temperature'\n",
    "}\n",
    "for data_frame in [X_co2_train, X_co2_test, X_ch4_train, X_ch4_test, X_co2_train_test, X_ch4_train_test]:\n",
    "    data_frame.rename(columns=col_name_map, inplace=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b8a6d97",
   "metadata": {},
   "source": [
    "## 2. Exploratory data analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17588462",
   "metadata": {},
   "source": [
    "#### Manual feature selection / analysis of features data\n",
    "\n",
    "doing EDA, it can also be used for checking multi co-linearity in data\n",
    "\n",
    "* Information gain\n",
    "* Correlation with target\n",
    "* Pairwise correlation\n",
    "* Variance threshold \n",
    "* ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ad96990",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check distribution of outputs\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2)\n",
    "fig.set_figwidth(6)\n",
    "fig.set_figheight(3)\n",
    "#fig.suptitle('Distributions of outputs variables')\n",
    "# CO2 emissions\n",
    "y_co2_train.plot(kind=\"hist\", bins=10, ax=ax1, color='red')\n",
    "ax1.title.set_text(\"CO$_2$ net unit emissions\")\n",
    "ax1.set_xlabel(\"Emission, gCO$_{2,eq}$/m$^2$/year\", fontsize = 11)\n",
    "ax1.set_ylabel(\"Frequency\", fontsize = 11)\n",
    "# CH4 emissions\n",
    "y_ch4_train.plot(kind=\"hist\", bins=10, ax=ax2, color='green')\n",
    "ax2.title.set_text(\"CH$_4$ net unit emission\")\n",
    "ax2.set_xlabel(\"Emission, gCO$_{2,eq}$/m$^2$/year\", fontsize = 11)\n",
    "ax2.set_ylabel(\"Frequency\", fontsize = 11)\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/data_exploration/co2_ch4_distributions.png'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46d1db84",
   "metadata": {},
   "source": [
    "### Find outliers - it's done via visual inspection  of histograms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fafaad0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "fig.set_figwidth(7)\n",
    "fig.set_figheight(3)\n",
    "ax.set_ylim([0,20])\n",
    "X_co2_train['catchment area'].hist(bins=75)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2965cb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find reservoirs with very large catchment areas indicated in the plot above\n",
    "large_c_area_indices = X_co2_train[['catchment area']].query('`catchment area` > 150000').index\n",
    "merged_df_min.loc[large_c_area_indices, :][['Name', 'co2_net', 'ch4_net', 'res_area', 'catch_area']]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7de1f06",
   "metadata": {},
   "source": [
    "## Find correlations among variables in the input data space and between inputs and the target variable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7884231",
   "metadata": {},
   "source": [
    "### A) For $CO_2$ emissions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61423137",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Xy_co2_train = pd.concat([X_co2_train, y_co2_train], axis=1).rename(columns={\"co2_net\": \"Net CO2 emissions\"})\n",
    "corr_matrix = Xy_co2_train.corr()\n",
    "mask_matrix = np.triu(corr_matrix)\n",
    "plt.figure(figsize=(28, 16))\n",
    "sns.set(font_scale=1.1)\n",
    "heatmap = sns.heatmap(corr_matrix, vmin=-1, vmax=1, annot=True, cmap=\"YlGnBu\", mask=mask_matrix)\n",
    "heatmap.set_title('Correlation Matrix Heatmap - all features', fontdict={'fontsize':22}, pad=14);\n",
    "heatmap.tick_params(axis='both', which='major', labelsize=24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc7fa70e",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_threshold = 0.80 # Value below which correlations are not displayed (e.g. in correlation plots)\n",
    "\n",
    "corr_matrix_no_large_c_area = Xy_co2_train[~Xy_co2_train.index.isin(large_c_area_indices)].corr()\n",
    "without_outliers: bool = True\n",
    "if without_outliers:\n",
    "    corr_matrix_filt_co2 = filter_corr_matrix(corr_matrix_no_large_c_area, threshold=corr_threshold)\n",
    "else:\n",
    "    corr_matrix_filt_co2 = filter_corr_matrix(corr_matrix, threshold=corr_threshold)\n",
    "fig_corr = plt.figure(figsize=(14, 8))\n",
    "sns.set(font_scale=1.1)\n",
    "mask_matrix = np.triu(corr_matrix_filt_co2)\n",
    "heatmap = sns.heatmap(corr_matrix_filt_co2, vmin=-1, vmax=1, annot=True, cmap=\"YlGnBu\", mask=mask_matrix)\n",
    "heatmap.set_title(\n",
    "    f'Feature Correlation Matrix for correlation coefficients > {corr_threshold}', pad=14, fontsize = 20)\n",
    "heatmap.tick_params(axis='both', which='major', labelsize=18)\n",
    "heatmap.tick_params(axis='x', rotation=90)\n",
    "fig_corr.savefig(\n",
    "    pathlib.Path('figures/data_exploration/feature_correlation_matrix.png'), \n",
    "    transparent = True, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5315d1",
   "metadata": {},
   "source": [
    "### A) For $CH_4$ emissions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23f854f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xy_ch4_train = pd.concat([X_ch4_train, y_ch4_train], axis=1).rename(\n",
    "    columns={\"ch4_net\": \"Net CH4 emissions\"})\n",
    "corr_matrix = Xy_ch4_train.corr()\n",
    "mask_matrix = np.triu(corr_matrix)\n",
    "plt.figure(figsize=(28, 16))\n",
    "sns.set(font_scale=1.1)\n",
    "heatmap = sns.heatmap(corr_matrix, vmin=-1, vmax=1, annot=True, cmap=\"YlGnBu\", mask=mask_matrix)\n",
    "heatmap.set_title('Correlation Matrix Heatmap - all features', fontdict={'fontsize':22}, pad=14);\n",
    "heatmap.tick_params(axis='both', which='major', labelsize=24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89dca5c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "corr_threshold = 0.80 # Value below which correlations are not displayed (e.g. in correlation plots)\n",
    "\n",
    "corr_matrix_no_large_c_area = Xy_ch4_train[~Xy_ch4_train.index.isin(large_c_area_indices)].corr()\n",
    "without_outliers: bool = True\n",
    "if without_outliers:\n",
    "    corr_matrix_filt_ch4 = filter_corr_matrix(corr_matrix_no_large_c_area, threshold=corr_threshold)\n",
    "else:\n",
    "    corr_matrix_filt_ch4 = filter_corr_matrix(corr_matrix, threshold=corr_threshold)\n",
    "    \n",
    "plt.figure(figsize=(14, 8))\n",
    "sns.set(font_scale=1.1)\n",
    "mask_matrix = np.triu(corr_matrix_filt_ch4)\n",
    "heatmap = sns.heatmap(corr_matrix_filt_ch4, vmin=-1, vmax=1, annot=True, cmap=\"YlGnBu\", mask=mask_matrix)\n",
    "heatmap.set_title(\n",
    "    f'Feature Correlation Matrix for correlation coefficients > {corr_threshold}', pad=14, fontsize = 20)\n",
    "heatmap.tick_params(axis='both', which='major', labelsize=18)\n",
    "heatmap.tick_params(axis='x', rotation=90)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19c44fdc",
   "metadata": {},
   "source": [
    "We should remove latitude and longitude form the list of features since they have not been used in the calculations. We might want to eliminate either P load or population as they're significantly correlated. Perhaps bottom temp and surface temp can be removed and only air temperature be used as an input variable since it correlates with both bottom and surface temperatures. Mean depth and max depth can be replaced with a single variable depth. **See plots below**:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c179defa",
   "metadata": {},
   "source": [
    "### Plot pairplots of selected features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e18bc5db",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_pairplot1:\n",
    "    X_co2_train_binned = X_co2_train.copy()\n",
    "    X_co2_train_binned['max_depth_bins'] = pd.cut(X_co2_train_binned['max depth'], bins=5)\n",
    "    pairplot = sns.pairplot(\n",
    "        X_co2_train_binned, vars=['air temperature', 'bottom temp', 'surface temp'], hue='max_depth_bins',\n",
    "        markers=\"+\",\n",
    "        kind='reg',\n",
    "        corner=True,\n",
    "        diag_kws= {'color': 'orange'})\n",
    "    fig.show()\n",
    "    pairplot.figure.savefig(\n",
    "        pathlib.Path('figures/data_exploration/temp_corr_pairplots.png'), \n",
    "        transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5668ef09",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note we plot the correlation for filtered data where we removed 'extremely' large catchments\n",
    "if plot_pairplot2:\n",
    "    catchment_area_threshold = 50_000\n",
    "    fig = plt.figure(figsize=(3, 3), dpi= 100, facecolor='w', edgecolor='k')\n",
    "    p_pop_pairplot = sns.pairplot(\n",
    "        X_co2_train.query(f\"`catchment area` < {catchment_area_threshold}\"), vars=['P load', 'population'],\n",
    "        markers=\"+\",\n",
    "        kind='reg',\n",
    "        diag_kind=\"hist\",\n",
    "        plot_kws={'line_kws':{'color':'black'},\n",
    "               'scatter_kws': {'alpha': 0.9,\n",
    "                               'color': 'green'}},\n",
    "        corner=True,\n",
    "        diag_kws= {'color': 'orange'})\n",
    "    p_pop_pairplot.fig.suptitle(\n",
    "        f\"Relationship between population and P load for catchments < {catchment_area_threshold} km2\",\n",
    "        y=1.04, fontsize = 12) \n",
    "    fig.show()\n",
    "    p_pop_pairplot.figure.savefig(\n",
    "        pathlib.Path('figures/data_exploration/pop_phosphorus_corr_pairplots.png'), \n",
    "        transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfd61e42",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_pairplot3:\n",
    "    fig = plt.figure(figsize=(2, 2), dpi= 100, facecolor='w', edgecolor='k')\n",
    "    urb_area_fractions_pairplot = sns.pairplot(\n",
    "        X_co2_train.query(\n",
    "            '`catchment urban area fraction` < 0.5'), \n",
    "        vars=['catchment urban area fraction', 'reservoir urban area fraction'],\n",
    "        markers=\"+\",\n",
    "        kind='reg',\n",
    "        corner=True,\n",
    "        diag_kws= {'color': 'orange'})\n",
    "    fig.show()\n",
    "    urb_area_fractions_pairplot.figure.savefig(\n",
    "        pathlib.Path('figures/data_exploration/res_catch_urban_fractions_corr_pairplots.png'), \n",
    "        transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a113aca",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_pairplot4:\n",
    "    fig = plt.figure(figsize=(3, 3), dpi= 100, facecolor='w', edgecolor='k')\n",
    "    min_max_depth_pairplot = sns.pairplot(\n",
    "        X_co2_train.query('`mean depth` < 20'), vars=['max depth', 'mean depth'],\n",
    "        markers=\"+\",\n",
    "        kind='reg',\n",
    "        corner=True,\n",
    "        diag_kws= {'color': 'orange'})\n",
    "    fig.show()\n",
    "    min_max_depth_pairplot.figure.savefig(\n",
    "        pathlib.Path('figures/data_exploration/min_max_depth_corr_pairplots.png'), \n",
    "        transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0dbef73e",
   "metadata": {},
   "outputs": [],
   "source": [
    "if plot_pairplot5:\n",
    "    fig = plt.figure(figsize=(3, 3), dpi= 100, facecolor='w', edgecolor='k')\n",
    "    area_volume_pairplot = sns.pairplot(\n",
    "        X_co2_train.query('`reservoir area` < 100'), vars=['reservoir volume', 'reservoir area'],\n",
    "        markers=\"+\",\n",
    "        kind='reg',\n",
    "        corner=True,\n",
    "        diag_kws= {'color': 'orange'})\n",
    "    fig.show()\n",
    "    area_volume_pairplot.figure.savefig(\n",
    "        pathlib.Path('figures/data_exploration/area_volume_corr_pairplots.png'), \n",
    "        transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a30aa50",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_co2_train['air temperature'].hist(bins = 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32a7603c",
   "metadata": {},
   "source": [
    "### Observations:\n",
    "1. Unit CO$_2$ net emissions are normally distributed whilst unit CH$_4$ net emissions have an L shaped distribution (strong left skew). \n",
    "2. Dataset contains two outliers for reservoirs with VERY LARGE catchment areas (Ywathit and Mong Tong reservoirs). It is not clear if this is an error or it is physically possible to have such large catchment areas for reservoirs, e.g. large flat surfaces, etc.\n",
    "3. There are a few correlations between features, namely:\n",
    "  - Population vs P. load, although the correlation does not hold well for very small catchments, which are predominant in the dataset. \n",
    "  - Catchment urban area fraction vs. reservoir urban area fraction seem correlated but it's an arficact of sparase histogram of data. In reality, these two variables are not correlated for small fracton values as shown in one of the pair-plots above\n",
    "  - Similar to above correlation between mean depth and max depth is high for the entire dataset but as the distributions of data are left skewed, it falls down for eg. small depths. HOWEVER, we decided to drop max_depth as less significant for predictions and only use mean_depth as a measure of the depth of the reservoir.\n",
    "  - Correlations between air, bottom and surface temperatures are also high. The distributions are right skewed (J-shape). We decided to only use air temperature as a feature as we assume that reservoir temperature is positivekly correlated with air temperature and negatively correlated with depth\n",
    "4. Reservoir volume vs reservoir area - although the correlation seems high, it does not hold for small reservoir volumes / areas. Therefore, we keep both the volume and the surface area in the feature dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b9f96d7",
   "metadata": {},
   "source": [
    "## Feature Scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcee59ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_FEAT = 7\n",
    "best_features_f_reg_co2 = SelectKBest(score_func = f_regression, k=N_FEAT)\n",
    "best_features_f_info_co2 = SelectKBest(score_func = mutual_info_regression, k=N_FEAT)\n",
    "best_features_f_reg_ch4 = SelectKBest(score_func = f_regression, k=N_FEAT)\n",
    "best_features_f_info_ch4 = SelectKBest(score_func = mutual_info_regression, k=N_FEAT)\n",
    "# ============== CO2 ==============\n",
    "fit_f_reg_co2 = best_features_f_reg_co2.fit(X_co2_train, y_co2_train)\n",
    "# X_new_co2_f_reg = fit_f_reg_co2.transform(X_co2_train)\n",
    "fit_f_info_co2 = best_features_f_info_co2.fit(X_co2_train, y_co2_train)\n",
    "# X_new_co2_f_info = fit_f_info_co2.transform(X_co2_train)\n",
    "# ============== CH4 ==============\n",
    "fit_f_reg_ch4 = best_features_f_reg_ch4.fit(X_ch4_train, y_ch4_train)\n",
    "fit_f_info_ch4 = best_features_f_info_ch4.fit(X_ch4_train, y_ch4_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c96a3e74",
   "metadata": {},
   "outputs": [],
   "source": [
    "#importlib.reload(lib.utils)\n",
    "#from lib.utils import plot_scores\n",
    "# PLOT CO2 REGRESSION FEATURE SCORES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6be8fb54",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_features = 10\n",
    "fig, axs = plt.subplots(1, 2, figsize=(10,6))\n",
    "fig.suptitle(\"CO$_2$ regression feature scores\")\n",
    "for ix, ax in enumerate(axs.flat):\n",
    "    if ix == 0:\n",
    "        plot_scores(\n",
    "            fit_f_reg_co2, X_co2_train, n_features, title = \"F value\", ax=ax, tick_fontsize = 11,\n",
    "            title_fontsize = 13)\n",
    "    if ix == 1:\n",
    "        plot_scores(\n",
    "            fit_f_info_co2, X_co2_train, n_features, title = \"mutual information\", ax=ax, tick_fontsize = 11,\n",
    "            title_fontsize = 13)\n",
    "plt.tight_layout()\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/data_exploration/co2_regression_feature_scores.png'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac518c71",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 2, figsize=(10,6))\n",
    "fig.suptitle(\"CH$_4$ regression feature scores\")\n",
    "for ix, ax in enumerate(axs.flat):\n",
    "    if ix == 0:\n",
    "        plot_scores(\n",
    "            fit_f_reg_ch4, X_ch4_train, n_features, title = \"F value\", ax=ax, tick_fontsize = 11,\n",
    "            title_fontsize = 13)\n",
    "    if ix == 1:\n",
    "        plot_scores(\n",
    "            fit_f_info_ch4, X_ch4_train, n_features, title = \"mutual information\", ax=ax, tick_fontsize = 11,\n",
    "            title_fontsize = 13)\n",
    "plt.tight_layout()\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/data_exploration/ch4_regression_feature_scores.png'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b77fef5",
   "metadata": {},
   "source": [
    "#### SEEMS TO BE CRASHING THE KERNEL - DO NOT RUN\n",
    "```\n",
    "profile = ProfileReport(X_min_train_co2, title=\"Profiling Report\")\n",
    "profile.to_file(\"profile_report.html\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15c06265",
   "metadata": {},
   "source": [
    "# Run CATBOOST, XGBOOST AND LIGHTGBM REGRESSIONS\n",
    "\n",
    "Currently runs two regressions per regression model - one for CO2 emissions and one for CH4 emissions.\n",
    "Both regressions do not include any scenario options (Soil Type, Treatment Factor or Landuse intensity)\n",
    "Instead the regressions assume the regressions are done for data assumed constant for all reservoirs\n",
    "in Myanmar, i.e. Mineral Soil, Primary Treatment, Low Landuse Intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3408fc57",
   "metadata": {},
   "source": [
    "### Make fast and dirty fitting first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "584b1e2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "simu_type = \"local\"\n",
    "cat_features = [] # Removed as we're only fitting the model to mya data with preset-values of those\n",
    "#                   [\"Landuse intensity\", \"Treatment\"]\n",
    "if simu_type == \"local\":\n",
    "    co2_model = CatBoostRegressor(loss_function = 'RMSE', task_type=\"CPU\", iterations=5000)\n",
    "    ch4_model = CatBoostRegressor(loss_function = 'RMSE', task_type=\"CPU\", iterations=5000)\n",
    "elif simu_type == \"colab\":\n",
    "    co2_model = CatBoostRegressor(loss_function = 'RMSE', task_type=\"GPU\" )\n",
    "    ch4_model = CatBoostRegressor(loss_function = 'RMSE', task_type=\"GPU\" )\n",
    "co2_model_catboost_quick_path = pathlib.Path(\"bin/regression_models/co2_model_catboost_quick.cbm\")\n",
    "ch4_model_catboost_quick_path = pathlib.Path(\"bin/regression_models/ch4_model_catboost_quick.cbm\")\n",
    "\n",
    "saved_model_path = pathlib.Path(\"saved_models\")\n",
    "\n",
    "if rerun or not os.path.isfile(co2_model_catboost_quick_path):\n",
    "    co2_model.fit(X_co2_train, y_co2_train, cat_features = cat_features, silent=True)\n",
    "    if override:\n",
    "        if not saved_model_path.exists():\n",
    "            saved_model_path.mkdir()\n",
    "        co2_model.save_model(\n",
    "            saved_model_path / \"co2_model_catboost_quick.cbm\", format=\"cbm\")\n",
    "else:\n",
    "    co2_model.load_model(co2_model_catboost_quick_path) \n",
    "        \n",
    "if rerun or not os.path.isfile(ch4_model_catboost_quick_path):\n",
    "    ch4_model.fit(X_ch4_train, y_ch4_train, cat_features = cat_features, silent=True)\n",
    "    if override:\n",
    "        if not saved_model_path.exists():\n",
    "            saved_model_path.mkdir()\n",
    "        ch4_model.save_model(\n",
    "            saved_model_path / \"ch4_model_catboost_quick.cbm\", format=\"cbm\")\n",
    "else:\n",
    "    ch4_model.load_model(ch4_model_catboost_quick_path) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7194a238",
   "metadata": {},
   "source": [
    "#### Make quick check of the quality of the models"
   ]
  },
  {
   "cell_type": "raw",
   "id": "1ccfa7e1",
   "metadata": {},
   "source": [
    "model_check(\n",
    "    model=co2_model, \n",
    "    X_train = X_co2_train, X_test = X_co2_test, \n",
    "    y_train = y_co2_train, y_test = y_co2_test)"
   ]
  },
  {
   "cell_type": "raw",
   "id": "0832d267",
   "metadata": {},
   "source": [
    "model_check(\n",
    "    model=ch4_model, \n",
    "    X_train = X_ch4_train, X_test = X_ch4_test, \n",
    "    y_train = y_ch4_train, y_test = y_ch4_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fd59a33",
   "metadata": {},
   "source": [
    "## Boosted tree model fitting with hyperparameter tuning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3371486d",
   "metadata": {},
   "source": [
    "### Tune the XGBoost models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8071054a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CO2 regression\n",
    "# Change override to True to retune the model\n",
    "model_co2_xgboost = hypertune_model(\n",
    "    X_co2_train, y_co2_train, num_evals = 2_000, hypertuner=HyperTuner.XGBOOST,\n",
    "    file=os.path.join('bin', 'regression_models', 'xgboost_co2.pkl'), override = override)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ac2a62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrain on full data set (for model explainability analysis)\n",
    "model_co2_xgboost.fit(X_co2_train_test, y_co2_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8840688",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"XGBOOST CO2 MODEL REGRESSION STATISTICS\")\n",
    "print(\"---------------------------------------\")\n",
    "model_check(model=model_co2_xgboost, \n",
    "    X_train = X_co2_train, X_test = X_co2_test, \n",
    "    y_train = y_co2_train, y_test = y_co2_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57e149a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CH4 regression\n",
    "# Change override to True to retune the model\n",
    "model_ch4_xgboost = hypertune_model(\n",
    "    X_ch4_train, y_ch4_train, num_evals = 2_000, hypertuner=HyperTuner.XGBOOST,\n",
    "    file=os.path.join('bin', 'regression_models', 'xgboost_ch4.pkl'), override = override)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48a40095",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrain on full data set (for model explainability analysis)\n",
    "model_ch4_xgboost.fit(X_ch4_train_test, y_ch4_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e1ccce4",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"XGBOOST CH4 MODEL REGRESSION STATISTICS\")\n",
    "print(\"---------------------------------------\")\n",
    "model_check(model=model_ch4_xgboost, \n",
    "    X_train = X_ch4_train, X_test = X_ch4_test, \n",
    "    y_train = y_ch4_train, y_test = y_ch4_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1acfb56f",
   "metadata": {},
   "source": [
    "### Tune the LightGBM models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "536db1ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CO2 regression\n",
    "model_co2_lightgbm = hypertune_model(\n",
    "    X_co2_train, y_co2_train, num_evals = 1_000, hypertuner=HyperTuner.LIGHTGBM,\n",
    "    file=os.path.join('bin', 'regression_models', 'lightgbm_co2.pkl'), override = override)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7b828d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove warnings in the LightGBM CO2 regression model - ONLY WORKS FOR PRE-TRAINED MODELS\n",
    "# IF TRAINING NEW MODELS - COMMENT OUT AND SEE THE WARNINGS FIRST BEFORE TURNING SOME CONFLICTING\n",
    "# REGRESSION PARAMETERS OFF\n",
    "import lightgbm as lgb\n",
    "\n",
    "model_co2_lightgbm.min_child_samples = None\n",
    "model_co2_lightgbm.min_split_gain=None\n",
    "model_co2_lightgbm.subsample=None\n",
    "model_co2_lightgbm.boosting_type=None\n",
    "model_co2_lightgbm.colsample_bytree=None\n",
    "model_co2_lightgbm.reg_alpha = None\n",
    "model_co2_lightgbm.reg_lambda = None\n",
    "model_co2_lightgbm.params={'verbose': -1, 'verbose_eval' : -1}\n",
    "model_co2_lightgbm.free_raw_data=False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "410768a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrain on full data set (for model explainability analysis)\n",
    "#model_co2_lightgbm.predict_raw_score = False\n",
    "model_co2_lightgbm.metric = {'rmse'}\n",
    "model_co2_lightgbm.fit(X_co2_train_test, y_co2_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfaea8cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"LIGHTGBM CO2 MODEL REGRESSION STATISTICS\")\n",
    "print(\"---------------------------------------\")\n",
    "model_check(model=model_co2_lightgbm, \n",
    "    X_train = X_co2_train, X_test = X_co2_test, \n",
    "    y_train = y_co2_train, y_test = y_co2_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5310c56",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CH4 regression\n",
    "model_ch4_lightgbm = hypertune_model(\n",
    "    X_ch4_train, y_ch4_train, num_evals = 1_000, hypertuner=HyperTuner.LIGHTGBM,\n",
    "    file=os.path.join('bin', 'regression_models', 'lightgbm_ch4.pkl'), override = override)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a98f30fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_ch4_lightgbm.min_child_samples=None\n",
    "model_ch4_lightgbm.colsample_bytree=None\n",
    "model_ch4_lightgbm.boosting_type=None\n",
    "model_ch4_lightgbm.min_split_gain=None\n",
    "model_ch4_lightgbm.reg_alpha=None\n",
    "model_ch4_lightgbm.subsample=None\n",
    "model_ch4_lightgbm.reg_lambda=None\n",
    "model_ch4_lightgbm.data_sample_strategy='goss'\n",
    "model_ch4_lightgbm.params={'verbose': -1, 'verbose_eval' : -1}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0618dfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrain on full data set (for model explainability analysis)\n",
    "model_ch4_lightgbm.fit(X_ch4_train_test, y_ch4_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "004602d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"LIGHTGBM CH4 MODEL REGRESSION STATISTICS\")\n",
    "print(\"---------------------------------------\")\n",
    "model_check(model=model_ch4_lightgbm, \n",
    "    X_train = X_ch4_train, X_test = X_ch4_test, \n",
    "    y_train = y_ch4_train, y_test = y_ch4_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f2b01fd",
   "metadata": {},
   "source": [
    "### Tune the CATBoost models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fa0fe7a",
   "metadata": {},
   "source": [
    "#### Encountered problems with Catboost errors and therefore the models have not been fitted\n",
    "* It is possible that the errors are caused by certain combinations of parameters during hypertuning. It may be possible to rectify this problem by removing some hyperparameters in hypertune or by reducing ranges of some hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34999308",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CO2 regression\n",
    "model_co2_catboost = hypertune_model(\n",
    "    X_co2_train, y_co2_train, num_evals = 40, hypertuner=HyperTuner.CATBOOST,\n",
    "    file=os.path.join('bin', 'regression_models', 'catboost_co2.pkl'), override = override)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f70d2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_co2_catboost.params = {\n",
    "    'silent': True, 'verbose': False, 'logging_level': 'Silent',\n",
    "    'metric_period':100}\n",
    "model_co2_catboost.metric_period = 10000\n",
    "model_co2_catboost.logging_level = 'Silent'\n",
    "model_co2_catboost.verbose = False\n",
    "model_co2_catboost.silent = True\n",
    "# None of the above f**** work!!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b7b6d75",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrain on full data set (for model explainability analysis)\n",
    "model_co2_catboost.fit(X_co2_train_test, y_co2_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "285440b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"CATBOOST CO2 MODEL REGRESSION STATISTICS\")\n",
    "print(\"---------------------------------------\")\n",
    "model_check(model=model_co2_catboost, \n",
    "    X_train = X_co2_train, X_test = X_co2_test, \n",
    "    y_train = y_co2_train, y_test = y_co2_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46f2fbd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CH4 regression\n",
    "model_ch4_catboost = hypertune_model(\n",
    "    X_ch4_train, y_ch4_train, num_evals = 40, hypertuner=HyperTuner.CATBOOST,\n",
    "    file=os.path.join('bin', 'regression_models', 'catboost_ch4.pkl'), override = override)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b620305b",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_ch4_catboost.verbose = -1\n",
    "model_ch4_catboost.logging_level = 'Silent'\n",
    "# Retrain on full data set (for model explainability analysis)\n",
    "model_ch4_catboost.fit(X_ch4_train_test, y_ch4_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11649778",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"CATBOOST CH4 MODEL REGRESSION STATISTICS\")\n",
    "print(\"---------------------------------------\")\n",
    "model_check(model=model_ch4_catboost, \n",
    "    X_train = X_ch4_train, X_test = X_ch4_test, \n",
    "    y_train = y_ch4_train, y_test = y_ch4_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fface8e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot figure for the visual abstract\n",
    "fig, axs = plt.subplots(1,1, figsize=(6,3))\n",
    "plot_shap_feature_importances(\n",
    "                model_co2_xgboost, X_co2_train_test,\n",
    "                max_vars = 15,\n",
    "                title=\" \",\n",
    "                plot_type = 'bar', ax=axs)\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/model_explanation/model_shaps_for_graphical_abstract.png'), \n",
    "    dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04bfbfdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_feat_importances(\n",
    "        model, X_train, X_test, y_test, title: str = \"Feature importances\",\n",
    "        file_name: str| None = None, **kwargs) -> None:\n",
    "    fig, axs = plt.subplots(2, 2, figsize=(10,7))\n",
    "    fig.suptitle(title)\n",
    "    for ix, ax in enumerate(axs.flat):\n",
    "        if ix == 0:\n",
    "            plot_gini_feature_importances(\n",
    "                model, X_train, 15, \n",
    "                'GINI-based Feature Importances', ax = ax)\n",
    "        if ix == 1:\n",
    "            # Computed on test data\n",
    "            plot_permutation_feature_importances(\n",
    "                model, X_test, y_test, max_vars = 15,\n",
    "                n_repeats = 7,\n",
    "                title='Permutation-based Feature Importances', ax = ax)\n",
    "        if ix == 2:\n",
    "            plot_shap_feature_importances(\n",
    "                model, X_test, \n",
    "                max_vars = 15,\n",
    "                title='Mean SHAP values',\n",
    "                plot_type = 'bar', ax=ax)\n",
    "\n",
    "    fig.delaxes(axs[1,1])\n",
    "    plt.tight_layout()\n",
    "    if file_name:\n",
    "        fig.savefig(file_name, dpi = 300, bbox_inches='tight', **kwargs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d2187a4",
   "metadata": {},
   "source": [
    "# Feature Importances for the CO$_2$ regression models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a2be7c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Doubling arguments (see plot_feat_importances) is weird, I know, but there must have been a reason for \n",
    "# it, that I've already forgotten (TJ)\n",
    "plot_feat_importances(\n",
    "    model_co2_xgboost, X_co2_train_test, X_co2_train_test, y_co2_train_test, \n",
    "    title = \"Feature importances - XGBoost model - CO$_2$ emissions\",\n",
    "    file_name = pathlib.Path('figures/model_explanation/feature_importances_xgboost_co2.png'),\n",
    "    transparent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f7ff7cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_feat_importances(\n",
    "    model_co2_lightgbm, X_co2_train_test, X_co2_train_test, y_co2_train_test, \n",
    "    title = \"Feature importances - LightGBM model - CO$_2$ emissions\",\n",
    "    file_name = pathlib.Path('figures/model_explanation/feature_importances_lightgbm_co2.png'),\n",
    "    transparent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e4128c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_feat_importances(\n",
    "    model_co2_catboost, X_co2_train_test, X_co2_train_test, y_co2_train_test, \n",
    "    title = \"Feature importances - CatBoost model - CO$_2$ emissions\",\n",
    "    file_name = pathlib.Path('figures/model_explanation/feature_importances_catboost_co2.png'),\n",
    "    transparent=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46258368",
   "metadata": {},
   "source": [
    "# Feature Importances for the CH$_4$ regression models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab707cc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_feat_importances(\n",
    "    model_ch4_xgboost, X_ch4_train_test, X_ch4_train_test, y_ch4_train_test, \n",
    "    title = \"Feature importances - XGBoost model - CH4 emissions\",\n",
    "    file_name = pathlib.Path('figures/model_explanation/feature_importances_xgboost_ch4.png'),\n",
    "    transparent=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53fb65bf",
   "metadata": {},
   "source": [
    "### Temporary code to export feature importances for the next notebook - move into a separate function / class and remove this temporary code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "564faffc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# --------------------- XGBOOST CH4 model ----------------------\n",
    "feature_importance = permutation_importance(\n",
    "        model_ch4_xgboost, X_ch4_train_test, y_ch4_train_test, n_repeats = 5, random_state = 42)\n",
    "num_features = len(feature_importance.importances_mean)\n",
    "sorted_idx = np.argsort(feature_importance.importances_mean)[::-1][:num_features]\n",
    "importances_df = pd.DataFrame(\n",
    "    data=feature_importance.importances_mean[sorted_idx]).T\n",
    "importances_df.columns = X_ch4_train_test.columns[sorted_idx]\n",
    "importances_df.to_csv(pathlib.Path(\"intermediate/ave_feature_importances_xgbost_ch4.csv\"))\n",
    "# ------------------- End of temporary code --------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0aa085ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_feat_importances(\n",
    "    model_ch4_lightgbm, X_ch4_train_test, X_ch4_train_test, y_ch4_train_test, \n",
    "    title = \"Feature importances - LightGBM model - CH4 emissions\",\n",
    "    file_name = pathlib.Path('figures/model_explanation/feature_importances_lightgbm_ch4.png'),\n",
    "    transparent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6139a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_feat_importances(\n",
    "    model_ch4_catboost, X_ch4_train_test, X_ch4_train_test, y_ch4_train_test, \n",
    "    title = \"Feature importances - CATBoost model - CH$_4$ emissions\",\n",
    "    file_name = pathlib.Path('figures/model_explanation/feature_importances_catboost_ch4.png'),\n",
    "    transparent=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf553eb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Save all feature importances to a file\n",
    "shap_values_folder = pathlib.Path('intermediate/shap_values')\n",
    "if not shap_values_folder.exists():\n",
    "    shap_values_folder.mkdir()\n",
    "\n",
    "output_folder = shap_values_folder/'model_avg_feat_importances'\n",
    "if not output_folder.exists():\n",
    "    output_folder.mkdir()\n",
    "\n",
    "model_feat_container = {}\n",
    "feat_imp_type: str = 'permutation' # shap, gini\n",
    "\n",
    "model_data_maps = {\n",
    "    ('xgboost', 'co2') : (model_co2_xgboost, X_co2_train_test, y_ch4_train_test),\n",
    "    ('xgboost', 'ch4') : (model_ch4_xgboost, X_ch4_train_test, y_ch4_train_test),\n",
    "    ('lightgbm', 'co2') : (model_co2_lightgbm, X_co2_train_test, y_ch4_train_test),\n",
    "    ('lightgbm', 'ch4') : (model_ch4_lightgbm, X_ch4_train_test, y_ch4_train_test),\n",
    "    ('catboost', 'co2') : (model_co2_catboost, X_co2_train_test, y_ch4_train_test),\n",
    "    ('catboost', 'ch4') : (model_ch4_catboost, X_ch4_train_test, y_ch4_train_test)}\n",
    "\n",
    "for key, pars in model_data_maps.items():\n",
    "    feats, cols = model_feature_importances(pars[0], pars[1], pars[2], feature_type=feat_imp_type)\n",
    "    model_feat_container[key] = (feats, cols)\n",
    "\n",
    "with open(output_folder / 'model_feats.pkl', 'wb') as handle:\n",
    "    pickle.dump(model_feat_container, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14af4018",
   "metadata": {},
   "source": [
    "## Publication permutation feature importances with removed colinearities between features\n",
    "## Figure for CO$_2$ and CH$_4$ XGBoost models\n",
    "* We played a bit with removal of multicollinear fetures before creating model feature importance plots using clustering. However, the methodology in scipy's documentation implemented here, removes whole clusters, but we'd rather choose one of the features in the cluster and remove the remaining correlated ones because we'd still want one of the correlated features to be included in the feature space. The code needs more fine-tuning. In the meantime, we removed correlated features manually in the beginning of the scripts and used cluster separation threshold on 0.0 which basically means that no features / clusters of features are being removed.\n",
    "* Instead, we've removed littoral area fraction and used mean_depth as a proxy for littoral area fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6953252",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.cluster import hierarchy\n",
    "from scipy.spatial.distance import squareform\n",
    "from scipy.stats import spearmanr\n",
    "\n",
    "def plot_feature_correlation_clusters(data) -> np.ndarray:\n",
    "    \"\"\" \"\"\"\n",
    "    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))\n",
    "    corr = spearmanr(data).correlation\n",
    "\n",
    "    # Ensure the correlation matrix is symmetric\n",
    "    corr = (corr + corr.T) / 2\n",
    "    np.fill_diagonal(corr, 1)\n",
    "\n",
    "    # We convert the correlation matrix to a distance matrix before performing\n",
    "    # hierarchical clustering using Ward's linkage.\n",
    "    distance_matrix = 1 - np.abs(corr)\n",
    "    dist_linkage = hierarchy.ward(squareform(distance_matrix))\n",
    "    dendro = hierarchy.dendrogram(\n",
    "        dist_linkage, labels=data.columns.to_list(), ax=ax1, leaf_rotation=90\n",
    "    )\n",
    "    dendro_idx = np.arange(0, len(dendro[\"ivl\"]))\n",
    "\n",
    "    ax2.imshow(corr[dendro[\"leaves\"], :][:, dendro[\"leaves\"]])\n",
    "    ax2.set_xticks(dendro_idx)\n",
    "    ax2.set_yticks(dendro_idx)\n",
    "    ax2.set_xticklabels(dendro[\"ivl\"], rotation=\"vertical\")\n",
    "    ax2.set_yticklabels(dendro[\"ivl\"])\n",
    "    _ = fig.tight_layout()\n",
    "    return dist_linkage\n",
    "dist_linkage = plot_feature_correlation_clusters(X_co2_train_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9ffc8e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_threshold = 0.0\n",
    "\n",
    "# Reduce input features based on the results of clustering\n",
    "cluster_ids = hierarchy.fcluster(dist_linkage, cluster_threshold, criterion=\"distance\")\n",
    "cluster_id_to_feature_ids = defaultdict(list)\n",
    "for idx, cluster_id in enumerate(cluster_ids):\n",
    "    cluster_id_to_feature_ids[cluster_id].append(idx)\n",
    "selected_features = [v[0] for v in cluster_id_to_feature_ids.values()]\n",
    "selected_features_names = X_co2_train_test.columns[selected_features]\n",
    "\n",
    "# Get subsets of features \n",
    "X_co2_train_sel = X_co2_train[selected_features_names]\n",
    "X_co2_test_sel = X_co2_test[selected_features_names]\n",
    "X_co2_train_test_sel = X_co2_train_test[selected_features_names]\n",
    "\n",
    "X_ch4_train_sel = X_ch4_train[selected_features_names]\n",
    "X_ch4_test_sel = X_ch4_test[selected_features_names]\n",
    "X_ch4_train_test_sel = X_ch4_train_test[selected_features_names]\n",
    "\n",
    "# columns to drop\n",
    "columns_to_drop = []\n",
    "X_co2_train_sel = X_co2_train_sel.drop(columns=columns_to_drop)\n",
    "X_co2_test_sel = X_co2_test_sel.drop(columns=columns_to_drop)\n",
    "X_co2_train_test_sel = X_co2_train_test_sel.drop(columns=columns_to_drop)\n",
    "X_ch4_train_sel = X_ch4_train_sel.drop(columns=columns_to_drop)\n",
    "X_ch4_test_sel = X_ch4_test_sel.drop(columns=columns_to_drop)\n",
    "X_ch4_train_test_sel = X_ch4_train_test_sel.drop(columns=columns_to_drop)\n",
    "\n",
    "model_co2_xgboost.fit(X_co2_train_sel, y_co2_train)\n",
    "model_ch4_xgboost.fit(X_ch4_train_sel, y_ch4_train)\n",
    "mm1 = XGBRegressor()\n",
    "mm2 = XGBRegressor()\n",
    "mm1.fit(X_co2_train_sel, y_co2_train)\n",
    "mm2.fit(X_ch4_train_sel, y_ch4_train)\n",
    "print(\n",
    "    \"Baseline accuracy on test data with features removed - CO2 model:\"\n",
    "    f\" {mm1.score(X_co2_test_sel, y_co2_test):.2}\"\n",
    ")\n",
    "print(\n",
    "    \"Baseline accuracy on test data with features removed - CH4 model:\"\n",
    "    f\" {mm2.score(X_ch4_test_sel, y_ch4_test):.2}\"\n",
    ")\n",
    "print(\"selected_model\")\n",
    "print(\n",
    "    \"Baseline accuracy on test data with features removed - CO2 model:\"\n",
    "    f\" {model_co2_lightgbm.score(X_co2_test_sel, y_co2_test):.2}\"\n",
    ")\n",
    "print(\n",
    "    \"Baseline accuracy on test data with features removed - CH4 model:\"\n",
    "    f\" {model_ch4_lightgbm.score(X_ch4_test_sel, y_ch4_test):.2}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b8ab5d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make a figure with two subplots - left subplot shows CO2 model permutation feature importances\n",
    "# whilst the right subplot show the permutation feature importances for the CH4 model\n",
    "# Save the figure to svg, pdf and png files.\n",
    "import random\n",
    "plot_boxplot = True\n",
    "num_features: int = 8\n",
    "random_state = 100 #random. randint(1,1_000)\n",
    "model_co2 = model_co2_lightgbm\n",
    "model_ch4 = model_ch4_lightgbm\n",
    "\n",
    "sns.set_theme(style=\"white\")\n",
    "sns.set_context(\"paper\", rc={\"grid.linewidth\": 0.00})\n",
    "\n",
    "result_co2 = permutation_importance(\n",
    "    model_co2, X_co2_train_test_sel, y_co2_train_test, n_repeats=50, \n",
    "    random_state=random_state, n_jobs=4, \n",
    "    scoring='neg_root_mean_squared_error'\n",
    ")\n",
    "result_ch4 = permutation_importance(\n",
    "    model_ch4, X_ch4_train_test_sel, y_ch4_train_test, n_repeats=50, \n",
    "    random_state=random_state, n_jobs=4, \n",
    "    scoring='neg_root_mean_squared_error'\n",
    ")\n",
    "\n",
    "sorted_importances_idx_co2 = result_co2.importances_mean.argsort()[::-1][:num_features]#[::-1]\n",
    "sorted_importances_idx_ch4 = result_ch4.importances_mean.argsort()[::-1][:num_features]#[::-1]\n",
    "importances_co2 = pd.DataFrame(\n",
    "    result_co2.importances[sorted_importances_idx_co2].T,\n",
    "    columns=X_co2_train_test_sel.columns[sorted_importances_idx_co2],\n",
    ")\n",
    "importances_ch4 = pd.DataFrame(\n",
    "    result_ch4.importances[sorted_importances_idx_ch4].T,\n",
    "    columns=X_ch4_train_test_sel.columns[sorted_importances_idx_ch4],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b96b209",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2,figsize=(12,3))\n",
    "rename_labels: bool = False\n",
    "\n",
    "if plot_boxplot:\n",
    "    boxplot_co2 = sns.boxplot(\n",
    "        data = importances_co2, orient = \"h\", hue=None, width=0.30,\n",
    "        linewidth=0.75, ax=ax1,\n",
    "        showcaps=True,\n",
    "        flierprops={\"marker\": \"o\", \"markersize\": 2},\n",
    "        whiskerprops={\"color\": \"k\"},\n",
    "        boxprops={\"facecolor\": (.3, .6, .8, 1.0), \"edgecolor\": 'k'},\n",
    "        medianprops={\"color\": \"k\", \"linewidth\": 0.7}) \n",
    "    boxplot_ch4 = sns.boxplot(\n",
    "        data = importances_ch4, orient = \"h\", hue=None, width=0.30,\n",
    "        linewidth=0.75, ax=ax2,\n",
    "        showcaps=True,\n",
    "        flierprops={\"marker\": \"o\", \"markersize\": 2},\n",
    "        whiskerprops={\"color\": \"k\"},\n",
    "        boxprops={\"facecolor\": (.3, .6, .8, 1.0), \"edgecolor\": 'k'},\n",
    "        medianprops={\"color\": \"k\", \"linewidth\": 0.7}) \n",
    "\n",
    "barplot_co2 = sns.barplot(\n",
    "    data=importances_co2, orient=\"h\", ax=ax1,\n",
    "    facecolor = 'gainsboro', alpha=0.75,\n",
    "    linewidth=0.5, edgecolor='k', estimator = np.median,\n",
    "    width = 0.65, ci = 0, zorder=0) \n",
    "\n",
    "barplot_ch4 = sns.barplot(\n",
    "    data=importances_ch4, orient=\"h\", ax=ax2,\n",
    "    facecolor = 'gainsboro', alpha=0.75,\n",
    "    linewidth=0.5, edgecolor='k', estimator = np.median,\n",
    "    width = 0.65, ci = 0, zorder=0) \n",
    "\n",
    "if rename_labels:\n",
    "    new_labels_left = [\"fFr\", \"evap\", \"T_air\", \"pop\", \"slope\", \"hmean\", \"runoff\", \"fSr\"]\n",
    "    new_labels_right = [\"hmean\", \"Ac\", \"T_air\", \"pop\", \"V\", \"evap\", \"slope\", \"fWc\"]\n",
    "    ax1.set_yticklabels(new_labels_left, fontsize=12)\n",
    "    ax2.set_yticklabels(new_labels_right, fontsize=12)\n",
    "\n",
    "\n",
    "ax1.set_xlabel(\"RMSE loss after permutations, gCO$_{2e}$ m$^{-2}$ yr$^{-1}$\", fontsize=15)\n",
    "ax1.spines['left'].set_linewidth(0.5)\n",
    "ax1.spines['bottom'].set_linewidth(0.5)\n",
    "\n",
    "ax2.set_xlabel(\"RMSE loss after permutations, gCO$_{2e}$ m$^{-2}$ yr$^{-1}$\", fontsize=15)\n",
    "ax2.spines['left'].set_linewidth(0.5)\n",
    "ax2.spines['bottom'].set_linewidth(0.5)\n",
    "\n",
    "ax1.tick_params(axis='x', which='major', labelsize=13)\n",
    "ax1.tick_params(axis='y', which='major', labelsize=13)\n",
    "ax2.tick_params(axis='x', which='major', labelsize=13)\n",
    "ax2.tick_params(axis='y', which='major', labelsize=13)\n",
    "\n",
    "sns.despine(offset=10, trim=False)\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7315ef37",
   "metadata": {},
   "outputs": [],
   "source": [
    "# INFORMATION: We use lightgbm models for the creation of feature importance plots\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/model_explanation/permutation_feat_importance_pub.png'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/model_explanation/permutation_feat_importance_pub.svg'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')\n",
    "fig.savefig(\n",
    "    pathlib.Path('figures/model_explanation/permutation_feat_importance_pub.pdf'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "976f1e44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot one feature-importance plot for graphical abstract\n",
    "fig2, ax1 = plt.subplots(1,1,figsize=(11,4))\n",
    "boxplot_co2 = sns.boxplot(\n",
    "    data = importances_co2, orient = \"h\", hue=None, width=0.30,\n",
    "    linewidth=0.75, ax=ax1,\n",
    "    showcaps=True,\n",
    "    flierprops={\"marker\": \"o\", \"markersize\": 2},\n",
    "    whiskerprops={\"color\": \"k\"},\n",
    "    boxprops={\"facecolor\": (.3, .6, .8, 1.0), \"edgecolor\": 'k'},\n",
    "    medianprops={\"color\": \"k\", \"linewidth\": 0.7}) \n",
    "barplot_co2 = sns.barplot(\n",
    "    data=importances_co2, orient=\"h\", ax=ax1,\n",
    "    facecolor = 'gainsboro', alpha=0.75,\n",
    "    linewidth=0.5, edgecolor='k', estimator = np.median,\n",
    "    width = 0.65, ci = 0, zorder=0) \n",
    "\n",
    "ax1.set_xlabel(\"RMSE loss after permutations, gCO$_{2e}$ m$^{-2}$ yr$^{-1}$\", fontsize=18)\n",
    "ax1.spines['left'].set_linewidth(0.5)\n",
    "ax1.spines['bottom'].set_linewidth(0.5)\n",
    "\n",
    "ax1.tick_params(axis='x', which='major', labelsize=16)\n",
    "ax1.tick_params(axis='y', which='major', labelsize=16)\n",
    "\n",
    "sns.despine(offset=10, trim=False)\n",
    "fig2.tight_layout()\n",
    "fig2.savefig(\n",
    "    pathlib.Path('figures/model_explanation/permutation_feat_importance_abstract.svg'), \n",
    "    transparent = False, dpi = 300, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5329a6c",
   "metadata": {},
   "source": [
    "# Model and predictions explanation with DALEX"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44d38841",
   "metadata": {},
   "source": [
    "### Create DALEX explainers for all 6 models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34c75f6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_co2_xgboost = dx.Explainer(\n",
    "    model_co2_xgboost, X_co2_train_test, y_co2_train_test, \n",
    "    label='xgboost model co2 emissions') # Uses dalex model explainer\n",
    "exp_co2_lightgbm = dx.Explainer(\n",
    "    model_co2_lightgbm, X_co2_train_test, y_co2_train_test, \n",
    "    label='lightgbm model co2 emissions') # Uses dalex model explainer\n",
    "exp_co2_catboost = dx.Explainer(\n",
    "    model_co2_catboost, X_co2_train_test, y_co2_train_test, \n",
    "    label='catboost model co2 emissions') # Uses dalex model explainer\n",
    "\n",
    "exp_ch4_xgboost = dx.Explainer(\n",
    "    model_ch4_xgboost, X_ch4_train_test, y_ch4_train_test, \n",
    "    label='xgboost model ch4 emissions') # Uses dalex model explainer\n",
    "exp_ch4_lightgbm = dx.Explainer(\n",
    "    model_ch4_lightgbm, X_ch4_train_test, y_ch4_train_test, \n",
    "    label='lightgbm model ch4 emissions') # Uses dalex model explainer\n",
    "exp_ch4_catboost = dx.Explainer(\n",
    "    model_ch4_catboost, X_ch4_train_test, y_ch4_train_test, \n",
    "    label='catboost model ch4 emissions') # Uses dalex model explainer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9502bb6",
   "metadata": {},
   "source": [
    "## Model diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a12c4c91",
   "metadata": {},
   "outputs": [],
   "source": [
    "variable = \"y\"\n",
    "yvariable = \"residuals\"\n",
    "exp_co2_xgboost.model_diagnostics().plot(variable=variable, yvariable=yvariable)\n",
    "exp_co2_lightgbm.model_diagnostics().plot(variable=variable, yvariable=yvariable)\n",
    "exp_co2_catboost.model_diagnostics().plot(variable=variable, yvariable=yvariable)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a5fef5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "variable = \"y\"\n",
    "yvariable = \"residuals\"\n",
    "exp_ch4_xgboost.model_diagnostics().plot(variable=variable, yvariable=yvariable)\n",
    "exp_ch4_lightgbm.model_diagnostics().plot(variable=variable, yvariable=yvariable)\n",
    "exp_ch4_catboost.model_diagnostics().plot(variable=variable, yvariable=yvariable)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc9aad98",
   "metadata": {},
   "source": [
    "## Model performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfc1f333",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_co2_xgboost.model_performance()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fea7a26",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_co2_lightgbm.model_performance()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7554f44a",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_co2_catboost.model_performance()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e115c56",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ch4_xgboost.model_performance()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b79c579",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ch4_lightgbm.model_performance()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d15af21",
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_ch4_catboost.model_performance()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2d5fdbc",
   "metadata": {},
   "source": [
    "## Model explanations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cb56eda",
   "metadata": {},
   "source": [
    "### Variable importances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5a1c9b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "importances_type = 'permutational'\n",
    "no_permutations = 30\n",
    "n_processes = 4\n",
    "# Other attributes for the `model_parts` method\n",
    "# types: permutational, variable_importance, feature_importance, ratio, difference, shap_wrapper, \n",
    "# shap_explainer = 'TreeExplainer'\n",
    "# e.g. co2_lightgbm_shp_vals = exp_co2_lightgbm.model_parts(type='shap_wrapper', shap_explainer='TreeExplainer')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fc4fc3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "co2_xgboost_importances = exp_co2_xgboost.model_parts(\n",
    "    type = importances_type , keep_distributions = True, label=\"XGBoost CO2 emissions\", B=no_permutations,\n",
    "    processes=n_processes)\n",
    "co2_lightgbm_importances = exp_co2_lightgbm.model_parts(\n",
    "    type = importances_type, keep_distributions = True, label=\"LightGBM CO2 emissions\", B=no_permutations,\n",
    "    processes=n_processes)\n",
    "co2_catboost_importances = exp_co2_catboost.model_parts(\n",
    "    type = importances_type, keep_distributions = True, label=\"CATBoost CO2 emissions\", B=no_permutations,\n",
    "    processes=n_processes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faf9f327",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Experimental (comment it out)\n",
    "mp = exp_co2_xgboost.model_parts(type='shap_wrapper', shap_explainer_type=\"TreeExplainer\")\n",
    "mp.plot(plot_type='bar', axis_color='k', color='k', show=True, max_display=6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc6a5d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "co2_xgboost_importances.plot(\n",
    "    title=None, max_vars=10, bar_width=20, vertical_spacing=0, split='model', digits=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b629840c",
   "metadata": {},
   "outputs": [],
   "source": [
    "co2_xgboost_importances.plot([co2_lightgbm_importances, co2_catboost_importances], max_vars=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04cfb658",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ch4_xgboost_importances = exp_ch4_xgboost.model_parts(\n",
    "    type = importances_type , keep_distributions = True, label=\"XGBoost CH4 emissions\", B=no_permutations)\n",
    "ch4_lightgbm_importances = exp_ch4_lightgbm.model_parts(\n",
    "    type = importances_type, keep_distributions = True, label=\"LightGBM CH4 emissions\", B=no_permutations)\n",
    "ch4_catboost_importances = exp_ch4_catboost.model_parts(\n",
    "    type = importances_type, keep_distributions = True, label=\"CATBoost CH4 emissions\", B=no_permutations)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "754abd89",
   "metadata": {},
   "outputs": [],
   "source": [
    "ch4_xgboost_importances.plot([ch4_lightgbm_importances, ch4_catboost_importances], max_vars=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cb5a343",
   "metadata": {},
   "source": [
    "### Partial, local and accumulated dependence profiles for XGBoost, LightGBM and CATBoost models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f319f11f",
   "metadata": {},
   "outputs": [],
   "source": [
    "## TODO: add groups and make grouped profiles with keyword : value pair of groups = 'cat var col name'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2d43e02",
   "metadata": {},
   "source": [
    "### CO2 emissions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "909239d2",
   "metadata": {},
   "source": [
    "### 1. CO2 XGBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "580aaea1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot for visual abstract\n",
    "selected_variables_co2 = ['reservoir forest area fraction', 'evapotranspiration']\n",
    "# Partial dependence profiles\n",
    "pd_co2_xgboost = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    N=50, label = 'Partial dependence')\n",
    "# Local dependence profiles\n",
    "ld_co2_xgboost = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence')\n",
    "# Accumulated dependence profiles\n",
    "ad_co2_xgboost = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence')\n",
    "agg_plot = pd_co2_xgboost.plot([ld_co2_xgboost, ad_co2_xgboost], show=False, y_title=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "155c0518",
   "metadata": {},
   "outputs": [],
   "source": [
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08aa96f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_plot2 = copy.deepcopy(agg_plot)\n",
    "tick_font_size = 16\n",
    "label_font_size = 18\n",
    "agg_plot.update_layout(\n",
    "    font=dict(color='black'),\n",
    "    yaxis_title=\"prediction\",\n",
    "    legend_title=\"\"\n",
    ")\n",
    "\n",
    "agg_plot2.update_xaxes(\n",
    "    linewidth=2,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=1)\n",
    "agg_plot2.update_yaxes(\n",
    "    linewidth=2,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=1)\n",
    "agg_plot2.update_xaxes(\n",
    "    linewidth=2,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=2)\n",
    "agg_plot2.update_yaxes(\n",
    "    linewidth=2,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=2)\n",
    "agg_plot2.update_traces(\n",
    "    line_width=2, \n",
    "    opacity=1\n",
    ")\n",
    "agg_plot2['data'][0]['line']['color']=\"#799ed9\"\n",
    "agg_plot2['data'][1]['line']['color']='#89b38a'\n",
    "agg_plot2['data'][2]['line']['color']='#c7644c'\n",
    "agg_plot2['data'][3]['line']['color']=\"#799ed9\"\n",
    "agg_plot2['data'][4]['line']['color']='#89b38a'\n",
    "agg_plot2['data'][5]['line']['color']='#c7644c'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6437c79",
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_plot2.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "455b2ef8",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_variables_co2 = [\n",
    "    'reservoir forest area fraction', 'evapotranspiration', \n",
    "    'population', 'catchment runoff', 'air temperature', 'catchment slope']\n",
    "# Partial dependence profiles\n",
    "pd_co2_xgboost = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    N=50, label = 'Partial dependence XGBoost CO2 emissions')\n",
    "# Local dependence profiles\n",
    "ld_co2_xgboost = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence XGBoost CO2 emissions')\n",
    "# Accumulated dependence profiles\n",
    "ad_co2_xgboost = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence XGBoost CO2 emissions')\n",
    "pd_co2_xgboost.plot([ld_co2_xgboost, ad_co2_xgboost])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d9b3831",
   "metadata": {},
   "source": [
    "### 2. CO2 LightGBM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "424d36dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partial dependence profiles\n",
    "pd_co2_lightgbm = exp_co2_lightgbm.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    N=50, label = 'Partial dependence LightGBM CO2 emissions')\n",
    "# Local dependence profiles\n",
    "ld_co2_lightgbm = exp_co2_lightgbm.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence LightGBM CO2 emissions')\n",
    "# Accumulated dependence profiles\n",
    "ad_co2_lightgbm = exp_co2_lightgbm.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence LightGBM CO2 emissions')\n",
    "pd_co2_lightgbm.plot([ld_co2_lightgbm, ad_co2_lightgbm])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fb78e45",
   "metadata": {},
   "source": [
    "### 3. CO2 CATBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "402cf4fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partial dependence profiles\n",
    "pd_co2_catboost = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    N=50, label = 'Partial dependence CATBoost CO2 emissions')\n",
    "# Local dependence profiles\n",
    "ld_co2_catboost = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence CATBoost CO2 emissions')\n",
    "# Accumulated dependence profiles\n",
    "ad_co2_catboost = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence CATGBoost CO2 emissions')\n",
    "pd_co2_catboost.plot([ld_co2_catboost, ad_co2_catboost])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d0c490c",
   "metadata": {},
   "source": [
    "### 4. CH4 XGBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6af4c2ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_variables_ch4 = [\n",
    "   'mean depth', 'population', 'catchment area',\n",
    "    'air temperature', 'reservoir shrub area fraction']\n",
    "#  'littoral area fraction', 'N load', 'retention coefficient', \n",
    "\n",
    "# Partial dependence profiles\n",
    "pd_ch4_xgboost = exp_ch4_xgboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    N=50, label = 'Partial dependence XGBoost CH4 emissions')\n",
    "# Local dependence profiles\n",
    "ld_ch4_xgboost = exp_ch4_xgboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence XGBoost CH4 emissions')\n",
    "# Accumulated dependence profiles\n",
    "ad_ch4_xgboost = exp_ch4_xgboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence XGBoost CH4 emissions')\n",
    "pd_ch4_xgboost.plot([ld_ch4_xgboost, ad_ch4_xgboost])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aebf973",
   "metadata": {},
   "source": [
    "### 5. CH4 LightGBM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69bd7511",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partial dependence profiles\n",
    "pd_ch4_lightgbm = exp_ch4_lightgbm.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    N=50, label = 'Partial dependence LightGBM CH4 emissions')\n",
    "# Local dependence profiles\n",
    "ld_ch4_lightgbm = exp_ch4_lightgbm.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence LightGBM CH4 emissions')\n",
    "# Accumulated dependence profiles\n",
    "ad_ch4_lightgbm = exp_ch4_lightgbm.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence LightGBM CH4 emissions')\n",
    "pd_ch4_lightgbm.plot([ld_ch4_lightgbm, ad_ch4_lightgbm])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4ba3ed",
   "metadata": {},
   "source": [
    "### 6. CH4 CATBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d718f746",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Partial dependence profiles\n",
    "pd_ch4_catboost = exp_ch4_catboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    N=50, label = 'Partial dependence CATBoost CH4 emissions')\n",
    "# Local dependence profiles\n",
    "ld_ch4_catboost = exp_ch4_catboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence CATBoost CH4 emissions')\n",
    "# Accumulated dependence profiles\n",
    "ad_ch4_catboost = exp_ch4_catboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence CATBoost CH4 emissions')\n",
    "pd_ch4_catboost.plot([ld_ch4_catboost, ad_ch4_catboost])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ace83522",
   "metadata": {},
   "source": [
    "## Accumulated profiles for all models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd8f5009",
   "metadata": {},
   "source": [
    "### CO2 accumulated dependence profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80a16519",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a plot for visual abstract\n",
    "selected_variables_co2 = ['reservoir forest area fraction', 'evapotranspiration']\n",
    "# Partial dependence profiles\n",
    "pd_co2_xgboost = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    N=50, label = 'Partial dependence')\n",
    "# Local dependence profiles\n",
    "ld_co2_xgboost = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='conditional',\n",
    "    N=50, label = 'Local dependence')\n",
    "# Accumulated dependence profiles\n",
    "ad_co2_xgboost = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2,\n",
    "    type='accumulated',\n",
    "    N=50, label = 'Accumulated dependence')\n",
    "agg_plot = pd_co2_xgboost.plot([ld_co2_xgboost, ad_co2_xgboost], show=False, y_title=\"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40587852",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_variables_co2_vis_abstract = ['reservoir forest area fraction', 'air temperature']\n",
    "prof_lightgbm_co2 = exp_co2_lightgbm.model_profile(\n",
    "    variables=selected_variables_co2_vis_abstract,\n",
    "    type='accumulated',\n",
    "    label=\"LightGBM\",\n",
    "    N=50)\n",
    "\n",
    "prof_xgboost_co2 = exp_co2_xgboost.model_profile(\n",
    "    variables=selected_variables_co2_vis_abstract,\n",
    "    type='accumulated',\n",
    "    label=\"XGBoost\",\n",
    "    N=50)\n",
    "\n",
    "prof_catboost_co2 = exp_co2_catboost.model_profile(\n",
    "    variables=selected_variables_co2_vis_abstract,\n",
    "    type='accumulated',\n",
    "    label=\"CATBoost\",\n",
    "    N=50)\n",
    "\n",
    "plot_visabstract = prof_lightgbm_co2.plot(\n",
    "    [prof_xgboost_co2, prof_catboost_co2], show=False, y_title=\"\",\n",
    "    title=\"Aggregated Profiles - CO2 Models\") # type = 'accumulated', type = 'conditional', geom='profiles'\n",
    "plot_visabstract.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6cdeb9a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_plot2 = copy.deepcopy(plot_visabstract)\n",
    "tick_font_size = 16\n",
    "label_font_size = 18\n",
    "agg_plot2.update_layout(\n",
    "    width=800,\n",
    "    height=400,\n",
    "    font=dict(color='black'),\n",
    "    yaxis_title=\"prediction\",\n",
    "    legend_title=\"\"\n",
    ")\n",
    "line_width = 3\n",
    "\n",
    "agg_plot2.update_xaxes(\n",
    "    linewidth=line_width,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=1)\n",
    "agg_plot2.update_yaxes(\n",
    "    linewidth=line_width,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=1)\n",
    "agg_plot2.update_xaxes(\n",
    "    linewidth=line_width,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=2)\n",
    "agg_plot2.update_yaxes(\n",
    "    linewidth=line_width,\n",
    "    showgrid=False,\n",
    "    tickfont=dict(size=tick_font_size), title_font=dict(size=label_font_size), color='black', row=1, col=2)\n",
    "agg_plot2.update_traces(\n",
    "    line_width=line_width, \n",
    "    opacity=0.85\n",
    ")\n",
    "agg_plot2['data'][0]['line']['color']=\"#799ed9\"\n",
    "agg_plot2['data'][1]['line']['color']='#89b38a'\n",
    "agg_plot2['data'][2]['line']['color']='#c7644c'\n",
    "agg_plot2['data'][3]['line']['color']=\"#799ed9\"\n",
    "agg_plot2['data'][4]['line']['color']='#89b38a'\n",
    "agg_plot2['data'][5]['line']['color']='#c7644c'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d28210c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_plot2.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08393e3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_plot2.write_image(pathlib.Path(\"figures/model_explanation/dependency_profiles_visabstract.svg\"))\n",
    "agg_plot2.write_image(pathlib.Path(\"figures/model_explanation/dependency_profiles_visabstract.png\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34d74b21",
   "metadata": {},
   "source": [
    "### CH4 accumulated dependence profiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ed2fcba",
   "metadata": {},
   "outputs": [],
   "source": [
    "prof_lightgbm_ch4 = exp_ch4_lightgbm.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='accumulated',\n",
    "    N=50)\n",
    "prof_xgboost_ch4 = exp_ch4_xgboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='accumulated',\n",
    "    N=50)\n",
    "prof_catboost_ch4 = exp_ch4_catboost.model_profile(\n",
    "    variables=selected_variables_ch4,\n",
    "    type='accumulated',\n",
    "    N=50)\n",
    "\n",
    "prof_lightgbm_ch4.plot(\n",
    "    [prof_xgboost_ch4, prof_catboost_ch4]) # type = 'accumulated', type = 'conditional', geom='profiles'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6850ffb",
   "metadata": {},
   "source": [
    "# Instance explanations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5dcb16a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_index_by_name(name: str, df_full: pd.DataFrame = merged_df_min_prim_low) -> pd.Int64Index | None:\n",
    "    \"\"\"Uses full dataset with Name column to obtain an index of a row containing the input data for the\n",
    "    reservoir which can be used to select data in the train/test dataset, e.g. for inspecting variable\n",
    "    importance for each reservoir\"\"\"\n",
    "    ix = merged_df_min_prim_low[df_full['Name']==name].index\n",
    "    if not ix.empty:\n",
    "        return ix\n",
    "    else:\n",
    "        print(f\"Reservoir with name {name} not found\")\n",
    "        return None\n",
    "    \n",
    "def loc_index_to_iloc(loc_index: pd.Int64Index, data: pd.DataFrame = X_co2_train_test) -> int:\n",
    "    \"\"\" \"\"\"\n",
    "    loc_index_int = int(np.mean(loc_index))\n",
    "    return data.index.get_loc(loc_index_int)\n",
    "\n",
    "def reservoir_names(df_full: pd.DataFrame = merged_df_min_prim_low) -> List[str]:\n",
    "    return list(df_full['Name'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ecb64ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(reservoir_names())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e541ce60",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"We have {len(reservoir_names())} reservoirs ...\")\n",
    "print(\" -- \".join(reservoir_names()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49adcefd",
   "metadata": {},
   "source": [
    "### PICK RESERVOIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb18025a",
   "metadata": {},
   "outputs": [],
   "source": [
    "reservoir_name = 'Zawgyi II'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c6f4be4",
   "metadata": {},
   "source": [
    "## Feature Importances on an instance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529fc368",
   "metadata": {},
   "source": [
    "## Breakdown interactions for all reservoirs - may take some time to compute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94ea5112",
   "metadata": {},
   "outputs": [],
   "source": [
    "# May take a while to run - uses the dalex interface to calculate SHAP values\n",
    "from typing import List\n",
    "import pickle\n",
    "from ipywidgets import IntProgress\n",
    "from IPython.display import display\n",
    "\n",
    "def run_and_save_breakdown_interactions_via_dalex(\n",
    "        reservoir_list: List[str], \n",
    "        input_data: pd.DataFrame = X_co2_train_test,\n",
    "        input_data_trimmed: pd.DataFrame = X_co2_train_test_sel,\n",
    "        output_path: str = \"outputs/model_explanations\",\n",
    "        B: int = 50, \n",
    "        interaction_preference: int = 1,\n",
    "        random_state: int = 42) -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    def to_dataframe(shaps, reservoir_name) -> pd.DataFrame:\n",
    "        \"\"\" \"\"\"\n",
    "        df = \\\n",
    "            shaps.result[['contribution', 'variable_name']]\\\n",
    "            .groupby('variable_name').mean().T  \n",
    "        df['reservoir name'] = reservoir_name\n",
    "        return df\n",
    "    \n",
    "    shp_conversion_config = {\n",
    "        'breakdown_lightbm_co2': (exp_co2_lightgbm, 'CO2 emissions for '),\n",
    "        'breakdown_xgboost_co2': (exp_co2_xgboost, 'CO2 emissions for '),\n",
    "        'breakdown_catboost_co2': (exp_co2_catboost, 'CO2 emissions for '),\n",
    "        'breakdown_lightgbm_ch4': (exp_ch4_lightgbm, 'CH4 emissions for '),\n",
    "        'breakdown_catboost_ch4': (exp_ch4_catboost,'CH4 emissions for '),\n",
    "        'breakdown_xgboost_ch4': (exp_ch4_xgboost, 'CH4 emissions for ')\n",
    "    }\n",
    "    num_iter = len(reservoir_list) * len(shp_conversion_config)\n",
    "    f = IntProgress(min=0, max=num_iter) # instantiate the bar\n",
    "    display(f)\n",
    "    \n",
    "    print(\"Calculating breakdown interactions using DALEX\")\n",
    "    print(\"Note that the pre-calculated interaction values with DALEX can be found in `model_explanations_precalculated`\")\n",
    "\n",
    "    for identifier, parameters in shp_conversion_config.items():\n",
    "        if identifier == 'breakdown_xgboost_co2' or identifier == 'breakdown_xgboost_ch4':\n",
    "            input_data = input_data_trimmed\n",
    "        else:\n",
    "            input_data = input_data\n",
    "        print(f\"Calculating breakdown interaction values for {identifier}...\")\n",
    "        # Initialise empty containers for data\n",
    "        shaps_dict = dict()\n",
    "        shaps_df = pd.DataFrame()\n",
    "        for reservoir_name in reservoir_list:\n",
    "            num_row = loc_index_to_iloc(find_index_by_name(name=reservoir_name), input_data)\n",
    "            input_reservoir = input_data.iloc[[num_row]]\n",
    "            #print(f\"Processing SHAP values for reservoir {reservoir_name}\")\n",
    "            shaps = parameters[0].predict_parts(\n",
    "                input_reservoir, 'break_down_interactions', \n",
    "                interaction_preference = interaction_preference, \n",
    "                label = f'{parameters[1]}{reservoir_name}',\n",
    "                B=B)\n",
    "            shaps_dict[reservoir_name] = shaps\n",
    "            # Add to a dataframe of shaps\n",
    "            shap_df = to_dataframe(shaps, reservoir_name)\n",
    "            shaps_df = pd.concat([shaps_df, shap_df])\n",
    "            f.value += 1\n",
    "        # Sanitise the dataframe\n",
    "        shaps_df.set_index('reservoir name', drop=True, inplace=True)\n",
    "        # Save the results\n",
    "        # Binary file with pickle\n",
    "        subfolder = f'interaction_preference_{interaction_preference}'\n",
    "        full_pickle_path = os.path.join(output_path, subfolder)\n",
    "        if not os.path.exists(full_pickle_path):\n",
    "            # Create the folder\n",
    "            os.makedirs(full_pickle_path)\n",
    "        pickle_file_path = os.path.join(full_pickle_path, identifier+'_dalex.pkl')\n",
    "        with open(pickle_file_path, 'wb') as fp:\n",
    "            pickle.dump(shaps_dict, fp)\n",
    "        # csv file with pandas\n",
    "        shaps_df.to_csv(os.path.join(full_pickle_path, identifier+'_dalex.csv'))\n",
    "        # xlsx file with pandas\n",
    "        shaps_df.to_excel(os.path.join(full_pickle_path, identifier+'_dalex.xlsx'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28821bbf",
   "metadata": {},
   "outputs": [],
   "source": [
    "recalculate_breakdown_interactions = False\n",
    "# BE CAREFUL WHEN RERUNNING BREAKDOWN INTERACTIONS - TAKES A VERY LONG TIME\n",
    "# WE ARE RUNNING BREAKDOWN INTERACTIONS FOR 4 DIFFERENT LEVELS OF INTERACTION PREFERENCE AND FOR 3 DIFFERENT MODELS\n",
    "\n",
    "if recalculate_breakdown_interactions == True:\n",
    "    run_and_save_breakdown_interactions_via_dalex(\n",
    "        reservoir_names(), interaction_preference = 0)\n",
    "    run_and_save_breakdown_interactions_via_dalex(reservoir_names(), interaction_preference = 1)\n",
    "    run_and_save_breakdown_interactions_via_dalex(reservoir_names(), interaction_preference = 2)\n",
    "    run_and_save_breakdown_interactions_via_dalex(reservoir_names(), interaction_preference = 3)\n",
    "else:\n",
    "    print(\"Breakdown interactions have not been recalculated\")\n",
    "    print(\"You can find pre-calculated values in `bin/model_explanations_precalculated/breakdown_interactions`\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1357b47",
   "metadata": {},
   "source": [
    "### Breakdown interactions (individual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27311082",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_predch4 = exp_ch4_lightgbm.predict_profile(input_reservoir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0d1baa2",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_predch4.result.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16a756bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1ch4 = output_predch4.plot(\n",
    "    variables = [\n",
    "        'Tair',\n",
    "        'hmean'], show=False, size=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8094dc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "p1ch4.update_layout(\n",
    "    xaxis=dict(\n",
    "        tickfont=dict(color='rgba(0, 0, 0, 0.8)', size=16),  # Set x-axis tick font color to black\n",
    "        title=dict(font=dict(size=18, color='black'))\n",
    "    ),\n",
    "    yaxis=dict(\n",
    "        tickfont=dict(color='rgba(0, 0, 0, 0.8)', size=16),  # Set y-axis tick font color to black\n",
    "        title=dict(font=dict(size=18, color='black'))\n",
    "    ),\n",
    "    xaxis2=dict(\n",
    "        tickfont=dict(color='rgba(0, 0, 0, 0.8)', size=16)  # Set x-axis tick font color to black\n",
    "    ),\n",
    "    yaxis2=dict(\n",
    "        tickfont=dict(color='rgba(0, 0, 0, 0.8)', size=16)  # Set y-axis tick font color to black\n",
    "    ),\n",
    "    title=dict(\n",
    "        text='Ceteris Paribus Plots for CH4 Aerial Emission Prediction',  # Main plot title\n",
    "        font=dict(size=18, color='black')  # Increase font size and change color to black for the main title\n",
    "    ),\n",
    "    font=dict(color='black'),\n",
    "    width=1000,\n",
    "    height=450\n",
    ")\n",
    "p1ch4.update_traces(\n",
    "    line=dict(width=3, color='rgba(0, 0, 0, 0.6)'),\n",
    "    opacity=0.70)\n",
    "p1ch4.show()\n",
    "p1ch4.write_image(\"figures/model_explanation/zagyi2cp4ch4.svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "009e2b73",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_pred = exp_co2_lightgbm.predict_profile(input_reservoir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb35ed92",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_pred.result.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37b684ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_pred.plot(variables = ['evapotranspiration', 'reservoir forest area fraction', 'mean depth', 'reservoir soil carbon'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "249f4ed7",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_row = loc_index_to_iloc(find_index_by_name(name=reservoir_name), X_co2_train_test)\n",
    "input_reservoir = X_co2_train_test.iloc[[num_row]]\n",
    "output_true = y_co2_train_test.iloc[num_row]\n",
    "output_pred = exp_co2_lightgbm.predict(input_reservoir)\n",
    "# Calculate, explain and plot the prediction using DALEX\n",
    "explanation_sample = exp_co2_lightgbm.predict_parts(\n",
    "    input_reservoir, type='break_down_interactions', interaction_preference = 1, \n",
    "    label=f'CO2 emissions for {reservoir_name}', B=25) \n",
    "# type=\"shap_wrapper\", type='break_down' keep_distributions = True does not have any effect\n",
    "explanation_sample.plot(max_vars=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e30b4af",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "num_row = loc_index_to_iloc(find_index_by_name(name=reservoir_name), X_ch4_train_test)\n",
    "input_reservoir = X_ch4_train_test.iloc[[num_row]]\n",
    "output_true = y_ch4_train_test.iloc[num_row]\n",
    "output_pred = exp_ch4_lightgbm.predict(input_reservoir)\n",
    "# Calculate, explain and plot the prediction using DALEX\n",
    "explanation_sample = exp_ch4_lightgbm.predict_parts(\n",
    "    input_reservoir, type='break_down_interactions', interaction_preference = 2, \n",
    "    label=f'CH4 emissions for {reservoir_name}') \n",
    "# type=\"shap_wrapper\", type='break_down' keep_distributions = True does not have any effect\n",
    "explanation_sample.plot(max_vars=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74241940",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "0ed77052",
   "metadata": {},
   "source": [
    "### SHAP values for individual reservoirs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d789c175",
   "metadata": {},
   "outputs": [],
   "source": [
    "# May take a while to run - uses the dalex interface to calculate SHAP values\n",
    "from typing import List\n",
    "import pickle\n",
    "from ipywidgets import IntProgress\n",
    "from IPython.display import display\n",
    "\n",
    "def run_and_save_shaps_via_dalex(\n",
    "        reservoir_list: List[str], \n",
    "        input_data: pd.DataFrame = X_co2_train_test,\n",
    "        output_path: str = \"outputs/model_explanations/dalex\",\n",
    "        B: int = 50, \n",
    "        random_state: int = 42) -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    def to_dataframe(shaps, reservoir_name) -> pd.DataFrame:\n",
    "        \"\"\" \"\"\"\n",
    "        df = \\\n",
    "            shaps.result[['contribution', 'variable_name']]\\\n",
    "            .groupby('variable_name').mean().T  \n",
    "        df['reservoir name'] = reservoir_name\n",
    "        return df\n",
    "    \n",
    "    shp_conversion_config = {\n",
    "        'shap_xgboost_co2': (exp_co2_xgboost, 'CO2 emissions for '),\n",
    "        'shap_lightbm_co2': (exp_co2_lightgbm, 'CO2 emissions for '),\n",
    "        'shap_catboost_co2': (exp_co2_catboost, 'CO2 emissions for '),\n",
    "        'shap_xgboost_ch4': (exp_ch4_xgboost, 'CH4 emissions for '),\n",
    "        'shap_lightgbm_ch4': (exp_ch4_lightgbm, 'CH4 emissions for '),\n",
    "        'shap_catboost_ch4': (exp_ch4_catboost,'CH4 emissions for ')\n",
    "    }\n",
    "    num_iter = len(reservoir_list) * len(shp_conversion_config)\n",
    "    f = IntProgress(min=0, max=num_iter) # instantiate the bar\n",
    "    display(f)\n",
    "    \n",
    "    print(\"Calculating instance-level SHAP values using DALEX\")\n",
    "    print(\"Note that the pre-calculated SHAP values with DALEX can be found in `model_explanations_precalculated`\")\n",
    "\n",
    "    for identifier, parameters in shp_conversion_config.items():\n",
    "        print(f\"Calculating SHAP values for {identifier}...\")\n",
    "        # Initialise empty containers for data\n",
    "        shaps_dict = dict()\n",
    "        shaps_df = pd.DataFrame()\n",
    "        for reservoir_name in reservoir_list:\n",
    "            num_row = loc_index_to_iloc(find_index_by_name(name=reservoir_name), input_data)\n",
    "            input_reservoir = input_data.iloc[[num_row]]\n",
    "            #print(f\"Processing SHAP values for reservoir {reservoir_name}\")\n",
    "            shaps = parameters[0].predict_parts(\n",
    "                input_reservoir, type='shap', \n",
    "                shap_explainer_type=\"TreeExplainer\",\n",
    "                keep_distributions=True,\n",
    "                processes=4,\n",
    "                label = f'{parameters[1]}{reservoir_name}',\n",
    "                B=B)\n",
    "            shaps_dict[reservoir_name] = shaps\n",
    "            # Add to a dataframe of shaps\n",
    "            shap_df = to_dataframe(shaps, reservoir_name)\n",
    "            shaps_df = pd.concat([shaps_df, shap_df])\n",
    "            f.value += 1\n",
    "        # Sanitise the dataframe\n",
    "        shaps_df.set_index('reservoir name', drop=True, inplace=True)\n",
    "        # Save the results\n",
    "        # Binary file with pickle\n",
    "        if not os.path.exists(output_path):\n",
    "            # Create the folder\n",
    "            os.makedirs(output_path)\n",
    "        pickle_file_path = os.path.join(output_path, identifier+'_dalex.pkl')\n",
    "        with open(pickle_file_path, 'wb') as fp:\n",
    "            pickle.dump(shaps_dict, fp)\n",
    "        # csv file with pandas\n",
    "        shaps_df.to_csv(os.path.join(output_path, identifier+'_dalex.csv'))\n",
    "        # xlsx file with pandas\n",
    "        shaps_df.to_excel(os.path.join(output_path, identifier+'_dalex.xlsx'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b2e9756",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Already pre-saved. Run only if you want to rerun all shaps again\n",
    "recalculate_shaps = False\n",
    "if recalculate_shaps:\n",
    "    run_and_save_shaps_via_dalex(reservoir_names())\n",
    "else:\n",
    "    print(\"SHAPS values have not been recalculated\")\n",
    "    print(\"You can find pre-calculated values in `bin/model_explanations_precalculated/dalex`\")    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a636e848",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_reservoir"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "679d5efd",
   "metadata": {},
   "source": [
    "## Make plots for visual abstract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89f9adb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_shaps(\n",
    "        reservoir_name: str, model: str = 'lightgbm', \n",
    "        file_location: str = \"figures/model_explanation/\") -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    if model == \"lightgbm\":\n",
    "        co2_explainer = exp_co2_lightgbm\n",
    "        ch4_explainer = exp_ch4_lightgbm\n",
    "    elif model == \"xgboost\":\n",
    "        co2_explainer = exp_co2_xgboost\n",
    "        ch4_explainer = exp_ch4_xgboost\n",
    "    elif model == \"catboost\":\n",
    "        co2_explainer = exp_co2_catboost\n",
    "        ch4_explainer = exp_ch4_catboost\n",
    "    else:\n",
    "        raise ValueError(f\"Model {model} not recognized.\")\n",
    "    res_location = loc_index_to_iloc(find_index_by_name(name=reservoir_name))\n",
    "    input_reservoir = X_co2_train_test.iloc[[res_location]]\n",
    "    explanation_sample_shap_co2 = co2_explainer.predict_parts(\n",
    "        input_reservoir, type='shap',\n",
    "        keep_distributions=True,\n",
    "        label = f'CO2 emissions for {reservoir_name}',\n",
    "        B=50,\n",
    "        processes=4,\n",
    "        random_state = 42)\n",
    "    exp_plot1 = explanation_sample_shap_co2.plot(\n",
    "    max_vars=6, title=\"\", bar_width=15, vertical_spacing = 0.05,\n",
    "    vcolors=(\"#799ed9\", '#89b38a', '#c7644c'), show=False) \n",
    "    # shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True\n",
    "    exp_plot1.update_layout(\n",
    "        xaxis=dict(\n",
    "            showgrid=False,  # Remove x-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set x-axis tick font color to black\n",
    "        ),\n",
    "        yaxis=dict(\n",
    "            showgrid=False,  # Remove y-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set y-axis tick font color to black\n",
    "        ),\n",
    "        font=dict(color='black')  # Set general font color to black\n",
    "    )\n",
    "    exp_plot1.show()\n",
    "    explanation_sample_shap_ch4 = ch4_explainer.predict_parts(\n",
    "        input_reservoir, type='shap', \n",
    "        keep_distributions=True,\n",
    "        processes=4,\n",
    "        label = f'CH4 Emission Intensity for {reservoir_name}',\n",
    "        B=50, random_state = 42)    \n",
    "    exp_plot2 = explanation_sample_shap_ch4.plot(\n",
    "        max_vars=6, title=\"\", bar_width=15, vertical_spacing = 0.05,\n",
    "        vcolors=(\"#799ed9\", '#89b38a', '#c7644c'), show=False) \n",
    "    # shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True\n",
    "    exp_plot2.update_layout(\n",
    "        xaxis=dict(\n",
    "            showgrid=False,  # Remove x-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set x-axis tick font color to black\n",
    "        ),\n",
    "        yaxis=dict(\n",
    "            showgrid=False,  # Remove y-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set y-axis tick font color to black\n",
    "        ),\n",
    "        font=dict(color='black')  # Set general font color to black\n",
    "    )\n",
    "    exp_plot2.show()\n",
    "    file_co2_svg = \"shap_\" + reservoir_name + \"_co2.svg\"\n",
    "    file_co2_png = \"shap_\" + reservoir_name + \"_co2.png\"\n",
    "    file_ch4_svg = \"shap_\" + reservoir_name + \"_ch4.svg\"\n",
    "    file_ch4_png = \"shap_\" + reservoir_name + \"_ch4.png\"    \n",
    "\n",
    "    exp_plot1.write_image(pathlib.Path(file_location) / file_co2_svg)\n",
    "    exp_plot1.write_image(pathlib.Path(file_location) / file_co2_png)\n",
    "    exp_plot2.write_image(pathlib.Path(file_location) / file_ch4_svg)\n",
    "    exp_plot2.write_image(pathlib.Path(file_location) / file_ch4_png)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5677921f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Rename columns for presentation purposes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "076d51cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "col_rename = {'catchment runoff': 'runoff',\n",
    " 'catchment area': 'Ac',\n",
    " 'population': 'pop',\n",
    " 'catchment bare soil fraction': 'fBS_c',\n",
    " 'catchment snow and ice fraction': 'fSI_c',\n",
    " 'catchment urban area fraction': 'fU_c',\n",
    " 'catchment water area fraction': 'fW_c',\n",
    " 'catchment wetland area fraction': 'fWt_c',\n",
    " 'catchment crop area fraction': 'fC_c',\n",
    " 'catchment shrub area fraction': 'fS_c',\n",
    " 'catchment forest area fraction': 'fF_c',\n",
    " 'catchment slope': 'slope',\n",
    " 'evapotranspiration': 'ET',\n",
    " 'catchment soil wetness': 'SWet',\n",
    " 'catchment mean olsen': 'OlsenP',\n",
    " 'reservoir volume': 'V',\n",
    " 'reservoir area': 'Ar',\n",
    " 'mean depth': 'hmean',\n",
    " 'reservoir soil carbon': 'soilC',\n",
    " 'reservoir mean radiance': 'Le',\n",
    " 'reservoir mean radiance may-sept': 'Le_May-Sept',\n",
    " 'reservoir mean radiance nov-mar': 'Le_Nov-Mar',\n",
    " 'reservoir mean monthly windspeed': 'vspeed',\n",
    " 'reservoir urban area fraction': 'fU_r',\n",
    " 'reservoir water area fraction': 'fW_r',\n",
    " 'reservoir wetland area fraction': 'fWt_r',\n",
    " 'reservoir crop area fraction': 'fC_r',\n",
    " 'reservoir shrub area fraction': 'fS_r',\n",
    " 'reservoir forest area fraction': 'fF_r',\n",
    " 'air temperature': 'Tair'}\n",
    "\n",
    "# Rename columns of the data and of the explained for visualisation purposes\n",
    "X_co2_train_test_renamed = X_co2_train_test.rename(\n",
    "    columns = col_rename)\n",
    "X_ch4_train_test_renamed = X_ch4_train_test.rename(\n",
    "    columns = col_rename)\n",
    "\n",
    "X_co2_train_test_renamed['V'] = X_co2_train_test_renamed['V'] / 1e6\n",
    "X_ch4_train_test_renamed['V'] = X_ch4_train_test_renamed['V'] / 1e6\n",
    "\n",
    "model_co2_xgboost.fit(X_co2_train_test_renamed, y_co2_train_test)\n",
    "model_ch4_xgboost.fit(X_ch4_train_test_renamed, y_ch4_train_test)\n",
    "model_co2_lightgbm.fit(X_co2_train_test_renamed, y_co2_train_test)\n",
    "model_ch4_lightgbm.fit(X_ch4_train_test_renamed, y_ch4_train_test)\n",
    "model_co2_catboost.fit(X_co2_train_test_renamed, y_co2_train_test)\n",
    "model_ch4_catboost.fit(X_ch4_train_test_renamed, y_ch4_train_test)\n",
    "# LightGBM explainers\n",
    "exp_co2_lightgbm = dx.Explainer(\n",
    "    model_co2_lightgbm, X_co2_train_test_renamed, y_co2_train_test, \n",
    "    label='CO2 net emission - LightGBM')\n",
    "exp_ch4_lightgbm = dx.Explainer(\n",
    "    model_ch4_lightgbm, X_ch4_train_test_renamed, y_ch4_train_test, \n",
    "    label='CH4 net emission - LightGBM')\n",
    "# CATBoost explainers\n",
    "exp_co2_catboost = dx.Explainer(\n",
    "    model_co2_catboost, X_co2_train_test_renamed, y_co2_train_test, \n",
    "    label='CO2 net emission - CATBoost')\n",
    "exp_ch4_catboost = dx.Explainer(\n",
    "    model_ch4_catboost, X_ch4_train_test_renamed, y_ch4_train_test, \n",
    "    label='CH4 net emission - CATBoost')\n",
    "# XGBoost explainers\n",
    "exp_co2_xgboost = dx.Explainer(\n",
    "    model_co2_xgboost, X_co2_train_test_renamed, y_co2_train_test, \n",
    "    label='CO2 net emission - XGBoost')\n",
    "exp_ch4_xgboost = dx.Explainer(\n",
    "    model_ch4_xgboost, X_ch4_train_test_renamed, y_ch4_train_test, \n",
    "    label='CH4 net emission - XGBoost')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d566eca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_breakdowns(\n",
    "        reservoir_name: str, model: str = 'lightgbm', \n",
    "        file_location: str = \"figures/model_explanation/\", interaction_preference: int = 1,\n",
    "        max_vars = 6,\n",
    "        rounding_digits: int = 3,\n",
    "        print_titles: bool = False,\n",
    "        input_data=X_co2_train_test_renamed, save_to_fig: bool = True) -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    if model == \"lightgbm\":\n",
    "        co2_explainer = exp_co2_lightgbm\n",
    "        ch4_explainer = exp_ch4_lightgbm\n",
    "    elif model == \"xgboost\":\n",
    "        co2_explainer = exp_co2_xgboost\n",
    "        ch4_explainer = exp_ch4_xgboost\n",
    "    elif model == \"catboost\":\n",
    "        co2_explainer = exp_co2_catboost\n",
    "        ch4_explainer = exp_ch4_catboost\n",
    "    else:\n",
    "        raise ValueError(f\"Model {model} not recognized.\")\n",
    "    #title_1 = f'Unit CO2 emission in gCO2e/m2/year - {reservoir_name}'\n",
    "    #title_2 = f'Unit CH4 Emission in gCO2e/m2/year - {reservoir_name}'\n",
    "    if print_titles:\n",
    "        title_1 = 'Unit carbon dioxide emission'\n",
    "        title_2 = \"Unit methane emission\"\n",
    "    else:\n",
    "        title_1, title_2 = \" \", \" \"\n",
    "    res_location = loc_index_to_iloc(find_index_by_name(name=reservoir_name))\n",
    "    input_reservoir = input_data.iloc[[res_location]]\n",
    "    explanation_sample_shap_co2 = co2_explainer.predict_parts(\n",
    "        input_reservoir, type='break_down_interactions', \n",
    "        interaction_preference = interaction_preference,\n",
    "        keep_distributions=True,\n",
    "        label = title_1,\n",
    "        B=50,\n",
    "        processes=4,\n",
    "        random_state = 42)\n",
    "    exp_plot1 = explanation_sample_shap_co2.plot(\n",
    "        max_vars=6, title=title_1, bar_width=15, vertical_spacing = 0.05,\n",
    "        digits = rounding_digits,\n",
    "        vcolors=(\"#2471a3\", '#89b38a', '#c7644c'), show=False) \n",
    "    # shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True\n",
    "    exp_plot1.update_layout(\n",
    "        paper_bgcolor='rgba(0,0,0,0)',\n",
    "        plot_bgcolor='rgba(0,0,0,0)',\n",
    "        xaxis=dict(\n",
    "            showgrid=False,  # Remove x-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set x-axis tick font color to black\n",
    "        ),\n",
    "        yaxis=dict(\n",
    "            showgrid=False,  # Remove y-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set y-axis tick font color to black\n",
    "        ),\n",
    "        font=dict(color='black')  # Set general font color to black\n",
    "    )\n",
    "    exp_plot1.update_traces(opacity=0.90)\n",
    "    exp_plot1.data[0].connector.line.color = 'black'\n",
    "    for shape in exp_plot1.layout.shapes:\n",
    "        if shape.type == 'line':\n",
    "            shape.line.color = '#424345'  # Set line color to black\n",
    "            shape.line.width = 2\n",
    "    exp_plot1.show()\n",
    "    explanation_sample_shap_ch4 = ch4_explainer.predict_parts(\n",
    "        input_reservoir, type='break_down_interactions', \n",
    "        interaction_preference = interaction_preference,\n",
    "        keep_distributions=True,\n",
    "        processes=4,\n",
    "        label = title_2,\n",
    "        B=50, random_state = 42)    \n",
    "    exp_plot2 = explanation_sample_shap_ch4.plot(\n",
    "        digits = rounding_digits,\n",
    "        max_vars=max_vars, title=title_2, bar_width=15, vertical_spacing = 0.05,\n",
    "        vcolors=(\"#2471a3\", '#89b38a', '#c7644c'), show=False) \n",
    "    #vcolors=(\"#799ed9\", '#89b38a', '#c7644c')\n",
    "    # shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True\n",
    "    exp_plot2.update_layout(\n",
    "        paper_bgcolor='rgba(0,0,0,0)',\n",
    "        plot_bgcolor='rgba(0,0,0,0)',\n",
    "        xaxis=dict(\n",
    "            showgrid=False,  # Remove x-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set x-axis tick font color to black\n",
    "        ),\n",
    "        yaxis=dict(\n",
    "            showgrid=False,  # Remove y-axis grid lines\n",
    "            tickfont=dict(color='black')  # Set y-axis tick font color to black\n",
    "        ),\n",
    "        font=dict(color='black')  # Set general font color to black\n",
    "    )\n",
    "    exp_plot2.update_traces(opacity=0.90)\n",
    "    exp_plot2.data[0].connector.line.color = 'black'\n",
    "    for shape in exp_plot2.layout.shapes:\n",
    "        if shape.type == 'line':\n",
    "            shape.line.color = '#424345'  # Set line color to black\n",
    "            shape.line.width = 2\n",
    "    exp_plot2.show()\n",
    "    file_co2_svg = reservoir_name + \"_breakdown_interactions\" + \"_co2.svg\"\n",
    "    file_co2_png = reservoir_name + \"_breakdown_interactions\" + \"_co2.png\"\n",
    "    file_ch4_svg = reservoir_name + \"_breakdown_interactions\" + \"_ch4.svg\"\n",
    "    file_ch4_png = reservoir_name + \"_breakdown_interactions\" + \"_ch4.png\"    \n",
    "\n",
    "    if not save_to_fig:\n",
    "        return\n",
    "    exp_plot1.write_image(pathlib.Path(file_location) / file_co2_svg)\n",
    "    exp_plot1.write_image(pathlib.Path(file_location) / file_co2_png)\n",
    "    exp_plot2.write_image(pathlib.Path(file_location) / file_ch4_svg)\n",
    "    exp_plot2.write_image(pathlib.Path(file_location) / file_ch4_png)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58f0e858",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\" -- \".join(sorted(list(set(merged_df['Name'])))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b31ca567",
   "metadata": {},
   "outputs": [],
   "source": [
    "reservoir_names = [\n",
    "    \"Thapanzeik\", \"Sedawgyi\", \"Zawgyi II\", \"Belin\", \"Laza\", \n",
    "    \"Mone Chaung\", \"Yeywa (upper)\", \"Mone Chaung\",\n",
    "    \"Kyee Ohn Kyee Wa\", \"Hawkham (upper)\", \"Myitsone\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da3b6714",
   "metadata": {},
   "outputs": [],
   "source": [
    "for reservoir_name in reservoir_names:\n",
    "    plot_breakdowns(\n",
    "        reservoir_name=reservoir_name,\n",
    "        interaction_preference = 1,\n",
    "        model = 'lightgbm')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68220c28",
   "metadata": {},
   "source": [
    "### This section requires all svg files to be generated including em intensity explanations that are calculated further on. Move this code to a new notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "266807be",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pprint\n",
    "# Collect plots for merging into composite figures\n",
    "directory = pathlib.Path(\"figures/model_explanation/\")\n",
    "# Initialize an empty dictionary to store the mappings\n",
    "reservoir_files = {}\n",
    "\n",
    "# Define the pattern for matching file names\n",
    "pattern = r'([^_]+)_breakdown_.*\\.svg'\n",
    "pattern = r'(.+)_breakdown_.*\\.svg'\n",
    "#pattern = r'([^_]+) breakdown.*\\.svg'\n",
    "pattern = r'(.+?)_breakdown_.*\\.svg'\n",
    "\n",
    "# Iterate over the files in the directory\n",
    "for filename in os.listdir(directory):\n",
    "    filepath = os.path.join(directory, filename)\n",
    "    if os.path.isfile(filepath):\n",
    "        # Check if the file matches the pattern\n",
    "        match = re.match(pattern, filename)\n",
    "        if match:\n",
    "            reservoir_name = match.group(1)\n",
    "            if reservoir_name == 'Thaphanseik':\n",
    "                reservoir_name = 'Thapanzeik'\n",
    "            # Add the file to the dictionary\n",
    "            if reservoir_name in reservoir_files:\n",
    "                reservoir_files[reservoir_name].append(filename)\n",
    "            else:\n",
    "                reservoir_files[reservoir_name] = [filename]\n",
    "\n",
    "# Print the dictionary\n",
    "pprint.pprint(reservoir_files)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "069ec793",
   "metadata": {},
   "outputs": [],
   "source": [
    "def combine_breakdown_plots(\n",
    "        reservoir_name: str, input_folder: str, input_files: List[str], output_file: pathlib.Path,\n",
    "        plot_offsets: Tuple[float, float, float], text_offsets: Tuple[float, float, float]) -> None:\n",
    "    \"\"\" \"\"\"\n",
    "    # TODO: Move all configs to here and later expose it to the caller\n",
    "    y_margin: int = 40\n",
    "    x_margin: int = 0\n",
    "    title_font_size: int = 19\n",
    "    subtitle_font_size: int = 16\n",
    "    #create new SVG figure\n",
    "    fig = sg.SVGFigure(\"10cm\", \"10cm\")\n",
    "    # Only accept three files (figures)\n",
    "    try:\n",
    "        assert len(input_files) == 3\n",
    "    except AssertionError:\n",
    "        raise ValueError(f\"Only three subplots supported. Entered {len(input_files)} files.\")\n",
    "    file_paths = [pathlib.Path(input_folder) / input_file for input_file in input_files]\n",
    "    # Load the fiure svg giles\n",
    "    fig1 = sg.fromfile(file_paths[0])\n",
    "    fig2 = sg.fromfile(file_paths[1])\n",
    "    fig3 = sg.fromfile(file_paths[2])\n",
    "    # get the plot objects\n",
    "    plot1 = fig1.getroot()\n",
    "    plot2 = fig2.getroot()\n",
    "    plot3 = fig3.getroot()\n",
    "    plot1.moveto(0, plot_offsets[0] + y_margin)\n",
    "    plot2.moveto(0, plot_offsets[1] + y_margin)\n",
    "    plot3.moveto(0, plot_offsets[2] + y_margin)\n",
    "    # add text labels\n",
    "    title = sg.TextElement(25,30, reservoir_name, size=title_font_size, weight=\"bold\")\n",
    "    txt1 = sg.TextElement(\n",
    "        25+x_margin, text_offsets[0] + y_margin, \n",
    "        \"Emission intensity\", size=subtitle_font_size, weight=\"bold\")\n",
    "    txt2 = sg.TextElement(\n",
    "        25+x_margin, text_offsets[1] + y_margin, \n",
    "        \"Unit carbon dioxide emissions\", size=subtitle_font_size, weight=\"bold\")\n",
    "    txt3 = sg.TextElement(\n",
    "        25+x_margin, text_offsets[2] + y_margin, \n",
    "        \"Unit methane emissions\", size=subtitle_font_size, weight=\"bold\")\n",
    "    # append plots and labels to figure\n",
    "    fig.append([plot1, plot2, plot3])\n",
    "    fig.append([title, txt1, txt2, txt3])\n",
    "    # save generated SVG files\n",
    "    fig.save(output_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b3a02c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "res_file_config = {\n",
    "    'Belin': {\n",
    "        \"files\" :\n",
    "        [\n",
    "             'Belin_breakdown_interactions_em_intensity.svg',\n",
    "             'Belin_breakdown_interactions_co2.svg',\n",
    "             'Belin_breakdown_interactions_ch4.svg'\n",
    "        ],\n",
    "        \"plot_offsets\": [-40, 160, 420],\n",
    "        \"text_offsets\": [30, 220, 480]\n",
    "    },\n",
    "    'Hawkham (Upper)': {\n",
    "        \"files\" : [\n",
    "            'Hawkham (upper)_breakdown_interactions_em_intensity.svg',\n",
    "            'Hawkham (upper)_breakdown_interactions_co2.svg',\n",
    "            'Hawkham (upper)_breakdown_interactions_ch4.svg'\n",
    "         ],\n",
    "        \"plot_offsets\": [-40, 160, 420],\n",
    "        \"text_offsets\": [30, 220, 480]\n",
    "    },\n",
    "     'Kyee Ohn Kyee Wa': {\n",
    "        \"files\" : [\n",
    "             'Kyee Ohn Kyee Wa_breakdown_interactions_em_intensity.svg',\n",
    "             'Kyee Ohn Kyee Wa_breakdown_interactions_co2.svg',\n",
    "             'Kyee Ohn Kyee Wa_breakdown_interactions_ch4.svg'\n",
    "        ],\n",
    "        \"plot_offsets\": [-40, 180, 440],\n",
    "        \"text_offsets\": [30, 240, 500]\n",
    "     },\n",
    "     'Laza': {\n",
    "        \"files\" : [\n",
    "            'Laza_breakdown_interactions_em_intensity.svg',\n",
    "            'Laza_breakdown_interactions_co2.svg',\n",
    "            'Laza_breakdown_interactions_ch4.svg'\n",
    "         ],\n",
    "        \"plot_offsets\": [-40, 165, 425],\n",
    "        \"text_offsets\": [30, 225, 485]\n",
    "     },\n",
    "     'Mone Chaung': {\n",
    "        \"files\" : [\n",
    "             'Mone Chaung_breakdown_interactions_em_intensity.svg',\n",
    "             'Mone Chaung_breakdown_interactions_co2.svg',\n",
    "             'Mone Chaung_breakdown_interactions_ch4.svg'\n",
    "        ],\n",
    "        \"plot_offsets\": [-40, 180, 440],\n",
    "        \"text_offsets\": [30, 240, 500]\n",
    "     },\n",
    "     'Sedawgyi': {\n",
    "        \"files\" : [\n",
    "            'Sedawgyi_breakdown_interactions_em_intensity.svg',\n",
    "            'Sedawgyi_breakdown_interactions_co2.svg',\n",
    "            'Sedawgyi_breakdown_interactions_ch4.svg'\n",
    "         ],\n",
    "        \"plot_offsets\": [-40, 165, 425],\n",
    "        \"text_offsets\": [30, 225, 485]\n",
    "     },\n",
    "     'Thapanseik': {\n",
    "        \"files\" : [\n",
    "            'Thaphanseik_breakdown_interactions_em_intensity.svg',\n",
    "            'Thapanzeik_breakdown_interactions_co2.svg',\n",
    "            'Thapanzeik_breakdown_interactions_ch4.svg'\n",
    "         ],\n",
    "        \"plot_offsets\": [-30, 150, 410],\n",
    "        \"text_offsets\": [30, 210, 470]\n",
    "     },\n",
    "     'Yeywa (Upper)': {\n",
    "        \"files\" : [\n",
    "            'Yeywa (upper)_breakdown_interactions_em_intensity.svg',\n",
    "            'Yeywa (upper)_breakdown_interactions_co2.svg',\n",
    "            'Yeywa (upper)_breakdown_interactions_ch4.svg'\n",
    "         ],\n",
    "        \"plot_offsets\": [-40, 180, 440],\n",
    "        \"text_offsets\": [30, 240, 500]\n",
    "     },\n",
    "     'Zawgyi II': {\n",
    "        \"files\" : [\n",
    "            'Zawgyi II_breakdown_interactions_em_intensity.svg',\n",
    "            'Zawgyi II_breakdown_interactions_co2.svg',\n",
    "            'Zawgyi II_breakdown_interactions_ch4.svg'\n",
    "         ],\n",
    "        \"plot_offsets\": [-25, 170, 430],\n",
    "        \"text_offsets\": [30, 230, 490]\n",
    "     }}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00e6ccb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_folder=\"figures/model_explanation\"\n",
    "\n",
    "for reservoir_name, config_data in res_file_config.items():\n",
    "    file_list = config_data['files']\n",
    "    plot_offsets = tuple(config_data[\"plot_offsets\"])\n",
    "    text_offsets = tuple(config_data[\"text_offsets\"])\n",
    "    file_name = reservoir_name + \"_breakdowns.svg\" \n",
    "    output_folder = pathlib.Path(input_folder, \"combined_breakdowns\")\n",
    "    output_folder.mkdir(parents=True, exist_ok=True)\n",
    "    output_file = output_folder / file_name\n",
    "    combine_breakdown_plots(\n",
    "        reservoir_name = reservoir_name,\n",
    "        input_folder = input_folder,\n",
    "        input_files = file_list,\n",
    "        output_file = output_file,\n",
    "        plot_offsets = plot_offsets,\n",
    "        text_offsets = text_offsets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3f47f29",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot_shaps(reservoir_name = \"Yeywa\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da196fad",
   "metadata": {},
   "source": [
    "## We haven't used any of the code or figures beyond this point..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef14094c",
   "metadata": {},
   "outputs": [],
   "source": [
    "reservoir_name = \"Lemro 2\"\n",
    "res_location = loc_index_to_iloc(find_index_by_name(name=reservoir_name))\n",
    "input_reservoir = X_co2_train_test.iloc[[res_location]]\n",
    "explanation_sample_shap = exp_co2_xgboost.predict_parts(\n",
    "    input_reservoir, type='shap',\n",
    "    keep_distributions=True,\n",
    "    processes=4,\n",
    "    label = f'CO2 Emission Intensity for {reservoir_name}',\n",
    "    B=50, random_state = 42)\n",
    "explanation_sample_shap.plot(max_vars=8) \n",
    "# shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87f2b227",
   "metadata": {},
   "outputs": [],
   "source": [
    "explanation_sample_shap = exp_ch4_xgboost.predict_parts(\n",
    "    input_reservoir, type='shap', \n",
    "    processes=4,\n",
    "    keep_distributions=True,\n",
    "    label = f'CH4 Emission Intensity for {reservoir_name}',\n",
    "    B=50, random_state = 42)\n",
    "explanation_sample_shap.plot(max_vars=8) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e6c7d15",
   "metadata": {},
   "source": [
    "## Back to the rest of the script..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83c2e23e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# How can the box plots be added using the Python version of DALEX?\n",
    "explanation_sample_shap = exp_co2_xgboost.predict_parts(\n",
    "    input_reservoir, type='shap',\n",
    "    processes=4,\n",
    "    keep_distributions=True,\n",
    "    label = f'CO2 emissions for {reservoir_name}',\n",
    "    B=50, random_state = 42)\n",
    "explanation_sample_shap.plot(max_vars=5) \n",
    "# shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8026fe79",
   "metadata": {},
   "outputs": [],
   "source": [
    "# How can the box plots be added using the Python version of DALEX?\n",
    "explanation_sample_shap = exp_co2_catboost.predict_parts(\n",
    "    input_reservoir, type='shap', \n",
    "    processes=4,\n",
    "    keep_distributions=True,\n",
    "    label = f'CO2 emissions for {reservoir_name}',\n",
    "    B=25)\n",
    "explanation_sample_shap.plot(max_vars=5) \n",
    "# shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7301eb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# How can the box plots be added using the Python version of DALEX?\n",
    "explanation_sample_shap = exp_ch4_lightgbm.predict_parts(\n",
    "    input_reservoir, type='shap', \n",
    "    processes=4,\n",
    "    keep_distributions=True,\n",
    "    label = f'CH4 emissions for {reservoir_name}',\n",
    "    B=50)\n",
    "explanation_sample_shap.plot(max_vars=5) \n",
    "# shap_explainer_type=\"TreeExplainer\" type=\"shap_wrapper\", type=\"break_dowo\", keep_distributions = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd216911",
   "metadata": {},
   "source": [
    "## Ceteris Paribus plots on an instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "342b002b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp_lightgbm_co2 = exp_co2_lightgbm.predict_profile(\n",
    "    input_reservoir, \n",
    "    variables=selected_variables_co2)\n",
    "cp_catboost_co2 = exp_co2_catboost.predict_profile(\n",
    "    input_reservoir, \n",
    "    variables=selected_variables_co2)\n",
    "cp_xgboost_co2 = exp_co2_xgboost.predict_profile(\n",
    "    input_reservoir, \n",
    "    variables=selected_variables_co2)\n",
    "cp_catboost_co2.plot([cp_xgboost_co2, cp_lightgbm_co2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "809a075f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cp_lightgbm_ch4 = exp_ch4_lightgbm.predict_profile(\n",
    "    input_reservoir, \n",
    "    variables=selected_variables_ch4)\n",
    "cp_catboost_ch4 = exp_ch4_catboost.predict_profile(\n",
    "    input_reservoir, \n",
    "    variables=selected_variables_co2)\n",
    "cp_xgboost_ch4 = exp_ch4_xgboost.predict_profile(\n",
    "    input_reservoir, \n",
    "    variables=selected_variables_ch4)\n",
    "cp_catboost_ch4.plot([cp_xgboost_ch4, cp_lightgbm_ch4])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfe65c56",
   "metadata": {},
   "source": [
    "# Model explanation with SHAP using the `SHAP` package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44ccf27a",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.plots.initjs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da385df0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CO2 explainers\n",
    "explainer_co2_xgboost = shap.TreeExplainer(model_co2_xgboost)\n",
    "explainer_co2_lightgbm = shap.TreeExplainer(model_co2_lightgbm)\n",
    "explainer_co2_catboost = shap.TreeExplainer(model_co2_catboost)\n",
    "# CH4 explainers\n",
    "explainer_ch4_xgboost = shap.TreeExplainer(model_ch4_xgboost)\n",
    "explainer_ch4_lightgbm = shap.TreeExplainer(model_ch4_lightgbm)\n",
    "explainer_ch4_catboost = shap.TreeExplainer(model_ch4_catboost)\n",
    "\n",
    "# SHAP VALUES - CO2\n",
    "shaps_co2_xgboost = explainer_co2_xgboost(X_co2_train_test, y_co2_train_test, check_additivity = True)\n",
    "shaps_co2_lightgbm = explainer_co2_lightgbm(X_co2_train_test, y_co2_train_test, check_additivity = True)\n",
    "shaps_co2_catboost = explainer_co2_catboost(X_co2_train_test, y_co2_train_test, check_additivity = True)\n",
    "# SHAP VALUES - CH4\n",
    "shaps_ch4_xgboost = explainer_ch4_xgboost(X_ch4_train_test, check_additivity = True)\n",
    "shaps_ch4_lightgbm = explainer_ch4_lightgbm(X_ch4_train_test, check_additivity = True)\n",
    "shaps_ch4_catboost = explainer_ch4_catboost(X_ch4_train_test, check_additivity = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2b068f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def validate_shaps(shaps, model, raw_score: bool = True):\n",
    "    \"\"\" \"\"\"\n",
    "    model_prediction = model.predict(shaps.data, raw_score=raw_score)\n",
    "    \n",
    "    # Test that mean model_prediction is equal to base_value\n",
    "    mean_model_prediction = np.mean(model_prediction)\n",
    "    assert isclose(mean_model_prediction, shaps.base_values[0], abs_tol=1e-6)\n",
    "    \n",
    "    # convert base vals vector to matrix\n",
    "    base_vals_matrix = shaps.base_values.repeat(shaps.data.shape[1]).reshape(shaps.data.shape)\n",
    "    \n",
    "    # Test that shaps add up to raw prediction\n",
    "    shap_predictions = np.sum(shaps.values, axis=1) + mean_model_prediction    \n",
    "    np.testing.assert_array_almost_equal(shap_predictions, model_prediction, decimal=6)\n",
    "    \n",
    "    # Test that shaps converted to real values match prediction\n",
    "    model_prediction_actual = model.predict(shaps.data, raw_score = False)\n",
    "    \n",
    "    # Find prediction from shaps\n",
    "    y_shap_actual = np.prod(np.exp(shaps.values), axis=1) * np.mean(model_prediction_actual, axis=0)\n",
    "    np.testing.assert_array_almost_equal(\n",
    "        y_shap_actual, model_prediction_actual\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "252333a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# validate_shaps(shaps_co2_lightgbm, model_co2_lightgbm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac34ede4",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.plots.waterfall(shaps_co2_lightgbm[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e24026fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check shap values from dalex and lightgbm shap output\n",
    "ix = find_index_by_name(reservoir_name)\n",
    "iloc_ix = loc_index_to_iloc(ix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5db2a9f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_model_input(data: pd.DataFrame, iloc_ix: int) -> pd.DataFrame:\n",
    "    \"\"\" \"\"\"\n",
    "    return data.iloc[iloc_ix].to_frame().transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5657410e",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_data = get_model_input(X_co2_train_test, iloc_ix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf51bd6c",
   "metadata": {},
   "outputs": [],
   "source": [
    "reservoir_name = 'Baingda Dam'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "679d0a19",
   "metadata": {},
   "outputs": [],
   "source": [
    "iloc_ix = loc_index_to_iloc(loc_index = find_index_by_name(reservoir_name), data = X_co2_train_test)\n",
    "\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(3, 1, sharex=False, sharey=False)\n",
    "fig.suptitle(f'SHAP values for CO$_2$ regression models for {reservoir_name} reservoir')\n",
    "plt.sca(ax1)\n",
    "shap.plots.waterfall(shaps_co2_xgboost[iloc_ix])\n",
    "ax1.title.set_text(\"XGBoost Regression Model\")\n",
    "plt.sca(ax2)\n",
    "shap.plots.waterfall(shaps_co2_lightgbm[iloc_ix])\n",
    "ax2.title.set_text(\"LightGBM Regression Model\")\n",
    "plt.sca(ax3)\n",
    "shap.plots.waterfall(shaps_co2_catboost[iloc_ix])\n",
    "ax3.title.set_text(\"CATBoost Regression Model\")\n",
    "fig.subplots_adjust(hspace=0)\n",
    "fig.set_figheight(10)\n",
    "fig.set_figwidth(14)\n",
    "fig.tight_layout()\n",
    "fig.savefig(pathlib.Path('figures/model_explanation/shap_values_per_reservoir_co2.png'))\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60b405bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "iloc_ix = loc_index_to_iloc(loc_index = find_index_by_name(reservoir_name), data = X_ch4_train_test)\n",
    "#fig, axes = plt.subplots(3)\n",
    "fig = plt.figure()\n",
    "fig.suptitle(f'SHAP values for CH$_4$ regression model for reservoir {reservoir_name}')\n",
    "plt.subplot(311)\n",
    "#fig.add_subplot(311)\n",
    "shap.plots.waterfall(shaps_ch4_xgboost[iloc_ix])\n",
    "plt.subplot(312)\n",
    "#fig.add_subplot(312)\n",
    "shap.plots.waterfall(shaps_ch4_lightgbm[iloc_ix])\n",
    "#fig.add_subplot(313)\n",
    "plt.subplot(313)\n",
    "shap.plots.waterfall(shaps_ch4_catboost[iloc_ix])\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "970690c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "explainer_co2_lightgbm.expected_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3630469",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "shap.plots.force(shaps_co2_lightgbm[iloc_ix, :])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23c6b3bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.plots.force(shaps_co2_catboost[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b9cbf64",
   "metadata": {},
   "source": [
    "### Plot beeswarm plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1462be1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "shap.plots.beeswarm(shaps_co2_xgboost, 15, axis_color='black', color=plt.get_cmap(\"viridis\"))\n",
    "fig.tight_layout()\n",
    "fig.savefig(pathlib.Path('figures/model_explanation/shap_beeswarm_co2_xgboost.png'),dpi=700)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e887120",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.summary_plot(shaps_co2_xgboost, X_co2_train_test )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7afef2c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.plots.bar(shaps_co2_catboost)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107e573b",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.plots.heatmap(shaps_co2_catboost[:], max_display = 42, plot_width=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47a8369a",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.plots.scatter(\n",
    "    shaps_ch4_catboost[:, \"retention coefficient\"])#,\n",
    "    #color=shaps_co2_catboost[:, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "075103fe",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7533677c",
   "metadata": {},
   "source": [
    "## Create dataframes with shap values and save them to files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80fe44a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_shp_to_dataframe(\n",
    "        shp_data: np.ndarray, \n",
    "        train_data: pd.DataFrame, \n",
    "        full_data: pd.DataFrame,\n",
    "        model = None,\n",
    "        relative: bool = False) -> pd.DataFrame:\n",
    "    \"\"\"Takes a numpy ndarray of shap values, information about column names and indices in\n",
    "    train_data, information about which reservoir matches which index in full data, and model (optional)\n",
    "    for predicting y_hat if the returned shap values should be in percentage terms relative to\n",
    "    the prediction\"\"\"\n",
    "    if relative:\n",
    "        # Calculate predictions\n",
    "        try:\n",
    "            y_hat = model.predict(train_data)\n",
    "        except AttributeError:\n",
    "            raise AttributeError(\"Model not provided or does not contain the predict method\")\n",
    "        # Get shap values in percentage\n",
    "        shap_data_scaled = (shp_data.T / y_hat * 100).T\n",
    "        shp_data = shap_data_scaled\n",
    "    \n",
    "    shaps_df = pd.DataFrame(\n",
    "        data=shp_data, index=X_co2_train_test.index, columns=X_co2_train_test.columns)\n",
    "    shaps_with_names = pd.concat(\n",
    "        [shaps_df, merged_df_min_prim_low['Name']], axis=1).set_index('Name', drop=True)\n",
    "\n",
    "    return shaps_with_names"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e9e41e5",
   "metadata": {},
   "source": [
    "## Convert all shap matrices to dataframes and store them in files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f87921f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "shp_conversion_config = {\n",
    "    'shap_xgboost_co2': (\n",
    "        shaps_co2_xgboost, X_co2_train_test, merged_df_min_prim_low, model_co2_xgboost),\n",
    "    'shap_lightbm_co2': (\n",
    "        shaps_co2_lightgbm, X_co2_train_test, merged_df_min_prim_low, model_co2_lightgbm),\n",
    "    'shap_catboost_co2': (\n",
    "        shaps_co2_catboost, X_co2_train_test, merged_df_min_prim_low, model_co2_catboost),\n",
    "    'shap_xgboost_ch4': (\n",
    "        shaps_ch4_xgboost, X_ch4_train_test, merged_df_min_prim_low, model_ch4_xgboost),\n",
    "    'shap_lightgbm_ch4': (\n",
    "        shaps_ch4_lightgbm, X_ch4_train_test, merged_df_min_prim_low, model_ch4_lightgbm),\n",
    "    'shap_catboost_ch4': (\n",
    "        shaps_ch4_catboost, X_ch4_train_test, merged_df_min_prim_low, model_ch4_catboost)\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a38ee73",
   "metadata": {},
   "source": [
    "## Save shap values calculated in the shap package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f4b288d",
   "metadata": {},
   "outputs": [],
   "source": [
    "for identifier, data in shp_conversion_config.items():\n",
    "    # Save absolute shap values\n",
    "    shp_df_absolute = convert_shp_to_dataframe(\n",
    "        shp_data = data[0].values,\n",
    "        train_data = data[1],\n",
    "        full_data = data[2],\n",
    "        model=data[3],\n",
    "        relative=False)\n",
    "    output_dir = os.path.join('outputs', 'model_explanations', 'shap')\n",
    "    if not os.path.exists(output_dir):\n",
    "        # Create the folder\n",
    "        os.makedirs(output_dir)\n",
    "    shp_df_absolute.to_csv(os.path.join(output_dir, identifier + '_absolute.csv'))\n",
    "    shp_df_absolute.to_excel(os.path.join(output_dir, identifier + '_absolute.xlsx'))\n",
    "    # Save percentage shap values\n",
    "    shp_df_relative = convert_shp_to_dataframe(\n",
    "        shp_data = data[0].values,\n",
    "        train_data = data[1],\n",
    "        full_data = data[2],\n",
    "        model=data[3],\n",
    "        relative=True)\n",
    "    shp_df_relative.to_csv(os.path.join(output_dir,identifier + '_relative.csv'))\n",
    "    shp_df_relative.to_excel(os.path.join(output_dir,identifier + '_relative.xlsx'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c4d40b",
   "metadata": {},
   "source": [
    "## Save shap values calculated in DALEX"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8401f95",
   "metadata": {},
   "source": [
    "Left for lated, if required..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a1e2921",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5f1545d7",
   "metadata": {},
   "source": [
    "# The END"
   ]
  }
 ],
 "metadata": {
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}