lesson_1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Science at IGFAE 2024\n",
    "## Lesson 1\n",
    "\n",
    "## Pietro Vischia (Universidad de Oviedo and ICTEA), pietro.vischia@cern.ch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment and run this if you are running on Colab (remove only the \"#\", keep the \"!\").\n",
    "# You can run it anyway, but it will do nothing if you have already installed all dependencies\n",
    "# (and it will take some time to tell you it is not gonna do anything)\n",
    "\n",
    "\n",
    "#from google.colab import drive\n",
    "#drive.mount('/content/drive')\n",
    "#%cd \"/content/drive/MyDrive/\"\n",
    "#! git clone https://github.com/vischia/data_science_school_igfae2024.git\n",
    "#%cd machine_learning_tutorial\n",
    "#!pwd\n",
    "#!ls\n",
    "##!pip install livelossplot shap\n",
    "#%pip install shap torchinfo livelossplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "import os\n",
    "import re\n",
    "import math\n",
    "import socket\n",
    "import json\n",
    "import pickle\n",
    "import gzip\n",
    "import copy\n",
    "import array\n",
    "import numpy as np\n",
    "import numpy.lib.recfunctions as recfunc\n",
    "\n",
    "from scipy.optimize import newton\n",
    "from scipy.stats import norm\n",
    "\n",
    "import uproot\n",
    "\n",
    "import datetime\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.utils import shuffle\n",
    "from sklearn.metrics import roc_curve, auc, accuracy_score\n",
    "from sklearn.tree import export_graphviz\n",
    "from sklearn.inspection import permutation_importance\n",
    "try:\n",
    "    # See #1137: this allows compatibility for scikit-learn >= 0.24\n",
    "    from sklearn.utils import safe_indexing\n",
    "except ImportError:\n",
    "    from sklearn.utils import _safe_indexing\n",
    "    \n",
    "from livelossplot import PlotLossesKeras\n",
    "\n",
    "from keras.losses import mean_squared_error, binary_crossentropy\n",
    "from keras.models import Sequential, Model\n",
    "from keras.layers import Input, Dense, Dropout, Concatenate\n",
    "from keras.layers import Lambda, Activation\n",
    "from keras.optimizers import SGD\n",
    "from tensorflow.keras.optimizers.legacy import Adam # for macos\n",
    "# from keras.optimizers import Adam # for non-macos\n",
    "from keras.regularizers import l2\n",
    "from tensorflow.keras.layers import BatchNormalization\n",
    "from tensorflow.keras.layers import LeakyReLU\n",
    "\n",
    "import keras\n",
    "from keras import backend as K\n",
    "import tensorflow as tf\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The data set\n",
    "\n",
    "We will use simulated events corresponding to three physics processes.\n",
    "\n",
    "- ttH production\n",
    "- ttW production\n",
    "- Drell-Yan ($pp\\to Z/\\gamma^*$+jets) production\n",
    "\n",
    "We will select the multilepton final state, which is a challenging final state with a rich structure and nontrivial background separation.\n",
    "\n",
    "<img src=\"figs/2lss.png\" alt=\"ttH multilepton 2lss\" style=\"width:40%;\"/>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "INPUT_FOLDER = './'\n",
    "HAVE_GPU = True\n",
    "# IMPORTANT: Uncomment this if you haven't installed the data yet\n",
    "#!mkdir data/; cd data/; wget https://www.hep.uniovi.es/vischia/cmsdas2024/ft_tth_multilep_igfae2024.tar.gz; tar xzvf ft_tth_multilep_igfae2024.tar.gz; mv igfae2024/* .; rmdir igfae2024; rm ft_tth_multilep_igfae2024.tar.gz; cd -;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "sig = uproot.open('data/signal.root')['Friends'].arrays(library=\"pd\")\n",
    "bk1 = uproot.open('data/background_1.root')['Friends'].arrays(library=\"pd\")\n",
    "bk2 = uproot.open('data/background_2.root')['Friends'].arrays(library=\"pd\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data inspection\n",
    "\n",
    "The first thing you need to do when building a machine learning model is to forget about the model, and **just look at the data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sig.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look in more detail at which features we have"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sig.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting histograms of some observables using matplotlib\n",
    "\n",
    "(see also examples on [matplotlib](https://matplotlib.org/3.5.3/api/_as_gen/matplotlib.pyplot.html) website)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "matplotlib.rcParams['figure.figsize'] = (8, 6)\n",
    "matplotlib.rcParams['axes.labelsize'] = 14\n",
    "\n",
    "\n",
    "plt.hist(sig[\"Hreco_Lep0_pt\"], bins=100)\n",
    "plt.xlabel(\"Lepton0 pT\")\n",
    "plt.ylabel(\"Events\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot one vs the other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(sig[\"Hreco_Lep0_pt\"], sig[\"Hreco_Lep1_pt\"])\n",
    "plt.xlabel(\"Lepton0 pT\")\n",
    "plt.ylabel(\"Lepton1 pT\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and another one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(sig[\"Hreco_Lep2_pt\"], bins=100)\n",
    "plt.xlabel(\"Lepton0 pT\")\n",
    "plt.ylabel(\"Events\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What is going on there???**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also check that \"event tag\" variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(sig[\"Hreco_evt_tag\"])\n",
    "plt.xlabel(\"Event tag\")\n",
    "plt.ylabel(\"Events\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What could we do?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = sig.drop([\"Hreco_Lep2_pt\", \"Hreco_Lep2_eta\", \"Hreco_Lep2_phi\", \"Hreco_Lep2_mass\", \"Hreco_evt_tag\"], axis=1 )\n",
    "signal.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to drop two more features, that are actually generation-level features (we will use them for regression tomorrow)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = signal.drop([\"Hreco_HTXS_Higgs_pt\", \"Hreco_HTXS_Higgs_y\"], axis=1 )\n",
    "signal.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also do the same for the backgrounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg1 = bk1.drop([\"Hreco_Lep2_pt\", \"Hreco_Lep2_eta\", \"Hreco_Lep2_phi\", \"Hreco_Lep2_mass\", \"Hreco_evt_tag\",\"Hreco_HTXS_Higgs_pt\", \"Hreco_HTXS_Higgs_y\"], axis=1 )\n",
    "bkg2 = bk2.drop([\"Hreco_Lep2_pt\", \"Hreco_Lep2_eta\", \"Hreco_Lep2_phi\", \"Hreco_Lep2_mass\", \"Hreco_evt_tag\",\"Hreco_HTXS_Higgs_pt\", \"Hreco_HTXS_Higgs_y\"], axis=1 )\n",
    "print(bkg2.columns)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Label events: signal = 1; background = 0\n",
    "Now let's put it all together, assigning some class labels. Let's start with signal and merge the backgrounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal['label'] = 1\n",
    "print(bkg1.shape, bkg2.shape)\n",
    "bkg = pd.concat([bkg1, bkg2])\n",
    "bkg['label'] = 0\n",
    "\n",
    "print(f\"bkg1 shape {bkg1.shape}\")\n",
    "print(f\"bkg2 shape {bkg2.shape}\")\n",
    "print(f\"bkg1+bkg2 shape {bkg.shape}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create a compound data set of signal & background together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.concat([signal,bkg])\n",
    "data.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This data set is still ordered, ie. all signal events are before the background events. ML training requires a shuffled data set instead!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.sample(frac=1).reset_index(drop=True)\n",
    "data.head(10)\n",
    "\n",
    "print(bkg.columns)\n",
    "print(data.isna().any())\n",
    "\n",
    "X = data.drop([\"label\"], axis=1)\n",
    "y = data[\"label\"]\n",
    "\n",
    "\n",
    "print(f\"data shape {data.shape}\")\n",
    "print(f\"input feature shape {X.shape}\")\n",
    "print(f\"label (=target) shape {y.shape}\")\n",
    "\n",
    "print(\"check for NAN numbers\")\n",
    "print(data.isna().any())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split the data set into training and test set\n",
    "\n",
    "When we train a machine learning algorithm, we are trying to solve an interpolation problem (*find the function of the input features that provides the best approximation of the true function*) by also requiring that the solution generalizes sufficiently well (*the interpolating function must also predict correctly the value of the true function for new, unseen data*).\n",
    "\n",
    "\n",
    "When we have a labelled dataset, we will therefore split it into: a *training set*, which we will use to train the machine learning algorithm; a *test set*, which we will use to evaluate the performance of the algorithm for various realizations of the algorithm (e.g. tuning hyperparameters); and an *application set*, which are the data we are really interested in studying in the end.\n",
    "\n",
    "For many applications, when the amount of hyperparameters tuning is moderate, application set and test set can be collapsed into a single set (usually called *test set*). This is what we will do in this tutorial.\n",
    "\n",
    "![Blah](figs/trainingNetwork.png)\n",
    "\n",
    "(Image: P. Vischia, [doi:10.5281/zenodo.6373442](https://doi.org/10.5281/zenodo.6373442))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the data set that we have loaded, we don't need to do any splitting, because the data set is already split into a training and test dataset.\n",
    "\n",
    "If our data set needed to be split, we would have to split it. Libraries like `scikit-learn` provide some ready-made functions that just do that:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(X, y, test_size=0.33, random_state=42)\n",
    "\n",
    "print(f\"We have {len(X_train)} training samples with {sum(y_train)} signal and {sum(1-y_train)} background events\")\n",
    "print(f\"We have {len(X_test)} testing samples with {sum(y_test)} signal and {sum(1-y_test)} background events\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect the data set\n",
    "\n",
    "The first and most important thing to do, when developing a ML application, is to inspect the data set.\n",
    "\n",
    "Looking at the data is crucial, because it may reveal the presence of malformed or missing data (e.g. when one feature has been filled with meaningless values for some or all of the data points) and may hint at whether it is necessary to preprocess the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize dataset\n",
    "\n",
    "import seaborn as sns\n",
    "sns.set()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols_to_plot = [\n",
    "    'Hreco_Lep1_pt',\n",
    "    'Hreco_HadTop_pt',\n",
    "    'Hreco_All5_Jets_pt',\n",
    "    'Hreco_More5_Jets_pt',\n",
    "    'Hreco_Jets_plus_Lep_pt',\n",
    "    'label'\n",
    "]\n",
    "pp=sns.pairplot(data=data.sample(1000)[cols_to_plot], hue='label', diag_kws={'bw_method': 0.2})\n",
    "pp.map_lower(sns.kdeplot, levels=4, color=\".2\") # Contours\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercise: what happens if you omit `.sample(100)` from the command above?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#pp=sns.pairplot(data=train_set[cols_to_plot], hue='label', diag_kws={'bw': 0.2})\n",
    "#pp.map_lower(sns.kdeplot, levels=4, color=\".2\")\n",
    "#plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Correlation matrices\n",
    "\n",
    "For classification problems, another important thing to take a look at is the correlation matrix between all the variables, in events belonging to each class separately.\n",
    "\n",
    "Looking at the correlation between features can highlight features that are strongly correlated with other ones, leading sometimes to dropping them (\"*if they are almost fully correlated, then including both does not add new information*\").\n",
    "\n",
    "We look at the correlation for each class because we are very interested in pairs of features that have different correlation in one class or the other (in our example, signal or backgroun)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def corrmatrix(corr, label):\n",
    "    ax = sns.heatmap(\n",
    "        corr, \n",
    "        vmin=-1., vmax=1., center=0.,\n",
    "        cmap=sns.diverging_palette(20., 220., n=200, as_cmap=True),\n",
    "        square=True\n",
    "    )\n",
    "    ax.set_xticklabels(\n",
    "        ax.get_xticklabels(),\n",
    "        rotation=45,\n",
    "        horizontalalignment='right'\n",
    "    );\n",
    "\n",
    "    ax.set_title('Correlation matrix for %s events' % label)\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "corrmatrix(X_train[y_train==1.].corr(), 'signal')\n",
    "corrmatrix(X_train[y_train==0.].corr(), 'background')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we have plotted is the Pearson correlation coefficient, which leads to an important limitation of this diagnostic tool.\n",
    "\n",
    "The Pearson correlation coefficient captures only **linear** correlation between variables, and is blind to many nonlinear correlations that there may be. Don't trust the above matrices blindly.\n",
    "\n",
    "\n",
    "![Figure from BDN2010](figs/corrcov.png)\n",
    "\n",
    "(figure from C. Delaere slides at the 2010 BND school)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preprocess the data\n",
    "\n",
    "Gradient minimization techniques like to mix numbers that are of the same order of magnitude and range.\n",
    "\n",
    "In this data set, features are already scaled such that the range falls in the neighbourhood of unity. Sometimes this happens naturally, but in this case several variables come directly from particle physics and represent momenta of particles produced in high-energy interaction: when expressed in GeV, these variables will most certainly **not** be in a range close to unity. Therefore, we can infer that these few variables have already been preprocessed.\n",
    "\n",
    "Common choices of preprocessing are minmax scaling or normalization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Minmax\n",
    "\n",
    "Compress the range linearly:\n",
    "\n",
    "$$X_{scaled} = \\frac{X-X_{min}}{X_{max}-X_{min}}$$\n",
    "\n",
    "A drawback is that this results in an artificially smaller variance (the range is compressed linearly), which can deform the effect of outliers.\n",
    "\n",
    "###### Standardization\n",
    "\n",
    "Compress the range and the shape:\n",
    "\n",
    "$$X_{normalized} = \\frac{X - \\mu}{\\sigma}$$\n",
    "\n",
    "where $\\mu$ is the mean of the feature values and $\\sigma$ is the standard deviation.\n",
    "\n",
    "###### Which one?\n",
    "\n",
    "Typically one would use minmax scaling when your features are remarkably nongaussian and your ML algorithm of choice doesn't require Gaussian inputs. The price is that it affects outliers.\n",
    "Typically one would use normalization when the features are approximately Gaussian or when your ML algorithm of choice requires Gaussian inputs. However, it also results in numbers close to 1 (minimization algorithms and gradient descend love numbers that are not too large or too small), so it can be used for any algorithm: the good news is that it doesn't affect outliers.\n",
    "\n",
    "For now, let's not apply any scaler.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build a simple tree-based classifier\n",
    "\n",
    "The act of selecting regions of a data set by \"cutting\" (imposing thresholds) on some of the features is very natural for the particle physicist, so let's start by training a tree-based classifier.\n",
    "\n",
    "Decision trees are precisely that:\n",
    "\n",
    "<img src=\"figs/bdt_en_edit.png\" alt=\"bdtexample\" style=\"width:80%;\"/>\n",
    "\n",
    "(image from [r2d3.us](http://www.r2d3.us/visual-intro-to-machine-learning-part-1/))\n",
    "\n",
    "Decision trees depend however on the random order of imposing cuts on the data set. To reduce the dependence on the starting point and the ordering of the selection process, we will use an *ensemble* of decision trees, which will cut at random the data set, and we will *pool* the classification answers from all the trees via e.g. a majority vote.\n",
    "\n",
    "Not all the trees resulting from the random cuts will make sense, so we will use a procedure called *boosting*, which consists in weighting each random tree based on its own performance. The weights will be used to decide how to generate the next trees.\n",
    "\n",
    "\n",
    "<img src=\"figs/boosting.png\" alt=\"boosting\" style=\"width:80%;\"/>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import AdaBoostClassifier, GradientBoostingClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "\n",
    "learning_rate = 1.0\n",
    "\n",
    "bdt_ada = AdaBoostClassifier(estimator=DecisionTreeClassifier(max_depth=3, criterion='log_loss'), n_estimators=100, learning_rate=learning_rate, random_state=42)\n",
    "bdt_grad = GradientBoostingClassifier(n_estimators=100, learning_rate=learning_rate, max_depth=3, random_state=42, verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If it takes too much to run, try downsampling:\n",
    "Ntrain=10000\n",
    "Ntest=2000\n",
    "X_train = X_train[:Ntrain]\n",
    "y_train = y_train[:Ntrain]\n",
    "X_test = X_test[:Ntest]\n",
    "y_test = y_test[:Ntest]\n",
    "fitted_bdt_ada=bdt_ada.fit(X_train, y_train)\n",
    "fitted_bdt_grad=bdt_grad.fit(X_train, y_train)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate the performance of the classifier\n",
    "\n",
    "A simple performance estimate for the classifier is the mean accuracy on a certain data set.\n",
    "Let's print it out for the training set and for the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Adaptive boost: train accuracy', fitted_bdt_ada.score(X_train, y_train),', test accuracy', fitted_bdt_ada.score(X_test, y_test))\n",
    "print('Gradient boost: train accuracy', fitted_bdt_grad.score(X_train, y_train),', test accuracy', fitted_bdt_grad.score(X_test, y_test))\n",
    "\n",
    "\n",
    "def get_loss_vs_iteration(x_features, y_labels, classifier):\n",
    "    return np.array(list(map(\n",
    "        lambda score: sklearn.metrics.log_loss(y_labels,score), classifier.staged_predict_proba(x_features)\n",
    "    )))\n",
    "\n",
    "train_loss_bdt_ada = get_loss_vs_iteration(X_train, y_train, fitted_bdt_ada)\n",
    "test_loss_bdt_ada = get_loss_vs_iteration(X_test, y_test, fitted_bdt_ada)\n",
    "\n",
    "\n",
    "train_loss_bdt_grad = get_loss_vs_iteration(X_train, y_train, fitted_bdt_grad)\n",
    "test_loss_bdt_grad = get_loss_vs_iteration(X_test, y_test, fitted_bdt_grad)\n",
    "\n",
    "\n",
    "plt.figure() \n",
    "plt.plot(np.arange(len(train_loss_bdt_ada)),train_loss_bdt_ada,label=\"Ada boost (train)\",color='royalblue',linestyle='-')\n",
    "plt.plot(np.arange(len(test_loss_bdt_ada)),test_loss_bdt_ada,label=\"Ada boost (test)\",color='royalblue',linestyle='--')\n",
    "\n",
    "plt.plot(np.arange(len(train_loss_bdt_grad)),train_loss_bdt_grad,label=\"Gradient boost (train)\",color='orange',linestyle='-')\n",
    "plt.plot(np.arange(len(test_loss_bdt_grad)),test_loss_bdt_grad,label=\"Gradient boost (test)\",color='orange',linestyle='--')\n",
    "         \n",
    "plt.ylabel(\"log(loss)\")\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.close()                                                                                                    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can observe two things:\n",
    "\n",
    "For each dataset, the performance of gradient-based boosting is better than of adaptive boosting, which is normally expected.\n",
    "\n",
    "For each algorithm, the performance on the test set is remarkably lower than that on the training set. Furthermore, the performance on the training set is very high. All of this is an indication that our algorithm is able to separate the two classes (signal and background) very effectively, but generalizes somehow poorly to data not seen during the training.\n",
    "\n",
    "To study the generalization properties of machine learning algorithms we will switch to neural networks, where we can control in an intuitive way the complexity of the network.\n",
    "\n",
    "For the moment, however, let's look at some additional ways of exploring the algorithm's performance.\n",
    "\n",
    "What happens if you change the learning rate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ROC curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_rocs(scores_and_names, y):\n",
    "    pack=[] \n",
    "    for s, n in scores_and_names: \n",
    "        fpr, tpr, thresholds = roc_curve(y.ravel(), s)\n",
    "        pack.append([n, fpr,tpr,thresholds])\n",
    "\n",
    "    plt.figure()\n",
    "    lw=2\n",
    "    for n, fpr, tpr, thresholds in pack:\n",
    "        plt.plot(fpr, tpr, lw=lw, label=\"%s (AUC = %0.2f)\" % (n, auc(fpr, tpr))) \n",
    "\n",
    "    plt.plot([0, 1], [0, 1], color=\"navy\", lw=lw, linestyle=\"--\")\n",
    "    plt.xlim([0.0, 1.0])\n",
    "    plt.ylim([0.0, 1.05])\n",
    "    plt.xlabel(\"False Positive Rate\")\n",
    "    plt.ylabel(\"True Positive Rate\")\n",
    "    plt.title(\"Receiver Operating Characteristic curve\")\n",
    "    plt.legend(loc=\"lower right\")\n",
    "    plt.show()\n",
    "\n",
    "    \n",
    "y_score = fitted_bdt_ada.decision_function(X_test)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now plot the ROC curve for our two classifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_rocs([ [fitted_bdt_ada.decision_function(X_test), 'AdaBoost'],\n",
    "            [fitted_bdt_grad.decision_function(X_test), 'GradBoost']],\n",
    "          y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspection of the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plots the model structure and it correlation with a dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn.tree\n",
    "\n",
    "plt.figure(figsize=[25,18])\n",
    "sklearn.tree.plot_tree(fitted_bdt_grad.estimators_[42, 0],feature_names=data.columns, fontsize=14)\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###Alternative method:# Pick one of the trees (maybe modify to pick the tree with largest weight or something like that)\n",
    "#sub_tree_42 = fitted_bdt_grad.estimators_[42, 0]\n",
    "## Visualization\n",
    "## Install graphviz: https://www.graphviz.org/download/\n",
    "#from pydotplus import graph_from_dot_data\n",
    "#from IPython.display import Image\n",
    "#dot_data = export_graphviz(\n",
    "#    sub_tree_42,\n",
    "#    out_file=None, filled=True, rounded=True,\n",
    "#    special_characters=True,\n",
    "#    proportion=False, impurity=False, # enable them if you want\n",
    "#)\n",
    "#graph = graph_from_dot_data(dot_data)\n",
    "#Image(graph.create_png())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Explainability\n",
    "\n",
    "A good question to pose yourself is whether the features you have chosen for your data are meaningful variables (i.e. if they are actually relevant to your classifier). Another good question is which features drive the prediction for a given event or set of events.\n",
    "\n",
    "All these questions can be answered by using different concepts:\n",
    "\n",
    "- **Permutation importance**: the decrease in a model score when a single feature value is randomly shuffled (scikit-learn docs) (akin to impacts for profile likelihood fits)\n",
    "Shapley values: based on game theory (see other contribution)\n",
    "Correlation-based: e.g. parallel coordinates in TMVA: look where each variable is mapped\n",
    "to/correlated with\n",
    "\n",
    "- **Perturbational approach**: perturbing the value of a feature and looking at the change of the prediction gives hints on how important the variable is for the method. This is at the basis of LIME (`pip install lime`). You can read more [here](https://www.oreilly.com/content/introduction-to-local-interpretable-model-agnostic-explanations-lime/)\n",
    "\n",
    "- **Game-theoretical approach**: consider the prediction task as a game in game theory, and the features as players who bet via their values. The payout, as difference of prediction with respect to the true value, estimates how much a feature pushes the prediction away from the truth.\n",
    "\n",
    "- **Visual approach**: parallel coordinates, which were implemented in ROOT TMVA, let you handily select a range in the prediction, and have a visual representation of which ranges of each feature is mapped into that region of the prediction.\n",
    "\n",
    "![Parallel coordinates (reference in the figure)](figs/parcoord.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Permutation importance\n",
    "\n",
    "The idea is: randomly shuffle one single feature value, then check how much does the prediction change. If the prediction decreases by a lot, then the value of the feature is crucial to the prediction.\n",
    "\n",
    "Note: the importance is always **relative to a specific model**, it has no absolute validity. A feature that is deemed low-importance for a badly designed model may be deemed high-importance for a good model, and viceversa. Permutation importance scores don't \"talk\" across different models.\n",
    "\n",
    "Also, if the model has a performance which is near-chance, then it is not strongly predictive, so the answers one may get from permutation importance scores are not really reliable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.inspection import permutation_importance\n",
    "\n",
    "scoring = ['r2', 'neg_mean_absolute_percentage_error', 'neg_mean_squared_error']\n",
    "\n",
    "rs = permutation_importance(\n",
    "    fitted_bdt_grad, X_test, y_test, n_repeats=30, random_state=0, scoring=scoring)\n",
    "\n",
    "for metric in rs:\n",
    "    print(f\"{metric}\")\n",
    "    r = rs[metric]\n",
    "    for i in r.importances_mean.argsort()[::-1]:\n",
    "        if r.importances_mean[i] - 2 * r.importances_std[i] > 0:\n",
    "            print(f\"    {X_train.columns[i]:<8}\" # +1 to skip the label in the naming\n",
    "                  f\"{r.importances_mean[i]:.3f}\"\n",
    "                  f\" +/- {r.importances_std[i]:.3f}\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Besides looking at the more important variables, you may also look at the less important, to **prune** them.\n",
    "\n",
    "Pruning consists in dropping the least important variables and retraining your machine learning algorithm.\n",
    "The idea behind it is that the variables dropped don't influence the prediction anyway, and retraining without them should give more or less the same performance but with a simpler model. Why would we do that? Well, for example inference may be time-sensitive, and simpler models are computationally **faster** to evaluate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Shapley values\n",
    "(based on [this blog post](https://www.analyticsvidhya.com/blog/2019/11/shapley-value-machine-learning-interpretability-game-theory/))\n",
    "\n",
    "Shapley values are a construct based on game theory.\n",
    "The main idea behind Shapley values is to consider the prediction task for a single event as game played by the feature values of that event. The features collaborate together to play the game by betting. The value of the feature is the amount each feature bets on the prediction task. The **Shapley Value** for each feature is the payout of the game, and consists in the correct weight such that the sum of all Shapley values for the features is the difference between the predictions and the average value of the model. In other words, the Shapley value represents how much each variable pushes the prediction far from the expected value.\n",
    "\n",
    "More concretely, the Shapley value for a feature A is computed as follows:\n",
    "\n",
    "- Get all subsets of features that do not contain A\n",
    "- Compute the effect of adding A to each of these subsets \n",
    "- Aggregate all the contributions (i.e. compute the marginal contribution of the feature over all the subsets)\n",
    "\n",
    "In principle we should retrain the model for each of these subsets, but instead we (the `shap` package, actually) will just compute predictions by replacing the value of the feature with its own average value.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shap\n",
    "shap.initjs()\n",
    "explainer = shap.TreeExplainer(fitted_bdt_grad)\n",
    "shap_values = explainer.shap_values(X_train)\n",
    "i = 500\n",
    "shap.force_plot(explainer.expected_value, shap_values[i], features=X_train.iloc[i], feature_names=X_train.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Values in blue represent features that push the prediction towards negative values, values in red represent features that push the prediction towards positive values, *for the event number 4776*.\n",
    "\n",
    "We naturally want a summary of Shapley values over all observations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap.summary_plot(shap_values, features=X_train, feature_names=X_train.columns, use_log_scale=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EXERCISE: Variable preprocessing\n",
    "\n",
    "Apply one of the scalers we have seen in class, or any other one from [this sklearn documentation page](https://scikit-learn.org/stable/auto_examples/preprocessing/plot_all_scaling.html).\n",
    "\n",
    "Once you have imported the scaler (see cell below), you can just apply the transform method like so:\n",
    "\n",
    "`std_scaled_X_train = StandardScaler().fit_transform(X)`. Make sure you apply this to both `X_train` and `X_test`\n",
    "\n",
    "Then create BDTs what use the transformed data instead of the original data, train them, and compare their performance (for instance using the ROC curve seen above) with the BDTs that use the transformed data.\n",
    "\n",
    "What is the result?\n",
    "\n",
    "Now create a BDT that uses the original data passed through a PCA tranformation (again, the `fit_transform`method works well), and compare it with a BDT that uses the standarditez data passed through a PCA transformation. What can you say about performance?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import (\n",
    "    MaxAbsScaler, # maxAbs\n",
    "    MinMaxScaler, # MinMax\n",
    "    Normalizer, # Normalization (equal integral)\n",
    "    StandardScaler# standard scaling\n",
    ")\n",
    "from sklearn.decomposition import PCA"
   ]
  },
  {
   "cell_type": "markdown",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}