jet_notebooks/2.JetTaggingMLP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/makagan/SSI_Projects/blob/main/jet_notebooks/2.JetTaggingMLP.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Rf258xT0XIwV"
   },
   "source": [
    "# Training a Jet Tagging with **DNN** \n",
    "\n",
    "---\n",
    "In this notebook, we perform a Jet identification task using a multiclass classifier based on a \n",
    "Dense Neural Network (DNN), also called multi-layer perceptron (MLP). The problem consists on identifying a given jet as a quark, a gluon, a W, a Z, or a top,\n",
    "based on set of physics-motivated high-level features.\n",
    "\n",
    "For details on the physics problem, see https://arxiv.org/pdf/1804.06913.pdf \n",
    "\n",
    "For details on the dataset, see Notebook1\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4OMAZgtyXIwY"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "import glob\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2lbB-J3hXIwb"
   },
   "source": [
    "# Preparation of the training and validation samples\n",
    "\n",
    "---\n",
    "In order to import the dataset, we now\n",
    "- clone the dataset repository (to import the data in Colab)\n",
    "- load the h5 files in the data/ repository\n",
    "- extract the data we need: a target and jetImage \n",
    "\n",
    "To type shell commands, we start the command line with !\n",
    "\n",
    "**nb, if you are running locally and you have already downloaded the datasets you can skip the cell below and, if needed, change the paths later to point to the folder with your previous download of the datasets.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jWjxFaRPXIwb"
   },
   "outputs": [],
   "source": [
    "! curl https://cernbox.cern.ch/s/zZDKjltAcJW0RB7/download -o Data-MLtutorial.tar.gz\n",
    "! tar -xvzf Data-MLtutorial.tar.gz \n",
    "! ls Data-MLtutorial/JetDataset/\n",
    "! rm Data-MLtutorial.tar.gz "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "cCGhrKdwXIwc"
   },
   "outputs": [],
   "source": [
    "target = np.array([])\n",
    "features = np.array([])\n",
    "# we cannot load all data on Colab. So we just take a few files\n",
    "datafiles = ['Data-MLtutorial/JetDataset/jetImage_7_100p_30000_40000.h5',\n",
    "             'Data-MLtutorial/JetDataset/jetImage_7_100p_60000_70000.h5',\n",
    "             'Data-MLtutorial/JetDataset/jetImage_7_100p_50000_60000.h5',\n",
    "             'Data-MLtutorial/JetDataset/jetImage_7_100p_10000_20000.h5',\n",
    "             'Data-MLtutorial/JetDataset/jetImage_7_100p_0_10000.h5']\n",
    "# if you are running locallt, you can use the full dataset doing\n",
    "# for fileIN in glob.glob(\"tutorials/HiggsSchool/data/*h5\"):\n",
    "for fileIN in datafiles:\n",
    "    print(\"Appending %s\" %fileIN)\n",
    "    f = h5py.File(fileIN)\n",
    "    myFeatures = np.array(f.get(\"jets\")[:,[12, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 48, 52]])\n",
    "    mytarget = np.array(f.get('jets')[0:,-6:-1])\n",
    "    features = np.concatenate([features, myFeatures], axis=0) if features.size else myFeatures\n",
    "    target = np.concatenate([target, mytarget], axis=0) if target.size else mytarget\n",
    "    f.close()\n",
    "print(target.shape, features.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6a333RYPXIwe"
   },
   "source": [
    "The dataset consists of 50000 jets, each represented by 16 features\n",
    "\n",
    "---\n",
    "\n",
    "We now shuffle the data, splitting them into a training and a validation dataset with 2:1 ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ZBqFs1eBXIwf"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_val, y_train, y_val = train_test_split(features, target, test_size=0.33)\n",
    "print(X_train.shape, X_val.shape, y_train.shape, y_val.shape)\n",
    "del features, target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GkNz5UAhXIwg"
   },
   "source": [
    "# DNN model building"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "tTSDOiEHXIwh"
   },
   "outputs": [],
   "source": [
    "# keras imports\n",
    "from tensorflow.keras.models import Model\n",
    "from tensorflow.keras.layers import Dense, Input, Dropout, Flatten, Activation\n",
    "from tensorflow.keras.utils import plot_model\n",
    "from tensorflow.keras import backend as K\n",
    "from tensorflow.keras import metrics\n",
    "from tensorflow.keras.callbacks import EarlyStopping, ReduceLROnPlateau, TerminateOnNaN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rAl0DZTxXIwi"
   },
   "outputs": [],
   "source": [
    "input_shape = X_train.shape[1]\n",
    "dropoutRate = 0.25"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2l492G8BXIwj"
   },
   "outputs": [],
   "source": [
    "####\n",
    "inputArray = Input(shape=(input_shape,))\n",
    "#\n",
    "x = Dense(40, activation='relu')(inputArray)\n",
    "x = Dropout(dropoutRate)(x)\n",
    "#\n",
    "x = Dense(20)(x)\n",
    "x = Activation('relu')(x)\n",
    "x = Dropout(dropoutRate)(x)\n",
    "#\n",
    "x = Dense(10, activation='relu')(x)\n",
    "x = Dropout(dropoutRate)(x)\n",
    "#\n",
    "x = Dense(5, activation='relu')(x)\n",
    "#\n",
    "output = Dense(5, activation='softmax')(x)\n",
    "####\n",
    "model = Model(inputs=inputArray, outputs=output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xu8rRUkhXIwj"
   },
   "outputs": [],
   "source": [
    "model.compile(loss='categorical_crossentropy', optimizer='adam')\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2HfKWoOtXIwk"
   },
   "source": [
    "We now train the model with these settings:\n",
    "\n",
    "- the **batch size** is a hyperparameter of gradient descent that controls the number of training samples to work through before the model internal parameters are updated\n",
    "    - batch size = 1 results in fast computation but noisy training that is slow to converge\n",
    "    - batch size = dataset size results in slow computation but faster convergence)\n",
    "\n",
    "- the **number of epochs** controls the number of complete passes through the full training dataset -- at each epoch gradients are computed for each of the mini batches and model internal parameters are updated.\n",
    "\n",
    "- the **callbacks** are algorithms used to optimize the training (full list [here](https://keras.io/api/callbacks/)):\n",
    "    - *EarlyStopping*: stop training when a monitored metric (`monitor`) has stopped improving in the last N epochs (`patience`)\n",
    "    - *ReduceLROnPlateau*: reduce learning rate when a metric (`monitor`) has stopped improving in the last N epochs (`patience`)\n",
    "    - *TerminateOnNaN*: terminates training when a NaN loss is encountered"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "KzO-lyLEXIwk"
   },
   "outputs": [],
   "source": [
    "batch_size = 128\n",
    "n_epochs = 50\n",
    "\n",
    "# train \n",
    "history = model.fit(X_train, y_train, epochs=n_epochs, batch_size=batch_size, verbose = 2,\n",
    "                validation_data=(X_val, y_val),\n",
    "                callbacks = [\n",
    "                EarlyStopping(monitor='val_loss', patience=10, verbose=1),\n",
    "                ReduceLROnPlateau(monitor='val_loss', factor=0.1, patience=2, verbose=1),\n",
    "                TerminateOnNaN()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "044bCLqVXIwl"
   },
   "outputs": [],
   "source": [
    "# plot training history\n",
    "plt.plot(history.history['loss'])\n",
    "plt.plot(history.history['val_loss'])\n",
    "plt.yscale('log')\n",
    "plt.title('Training History')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['training', 'validation'], loc='upper right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "oESSmNLxXIwm"
   },
   "source": [
    "# Building the ROC Curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "a_AROD6SXIwm"
   },
   "outputs": [],
   "source": [
    "labels = ['gluon', 'quark', 'W', 'Z', 'top']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gjKT7EjUXIwn"
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "predict_val = model.predict(X_val)\n",
    "df = pd.DataFrame()\n",
    "fpr = {}\n",
    "tpr = {}\n",
    "auc1 = {}\n",
    "\n",
    "plt.figure()\n",
    "for i, label in enumerate(labels):\n",
    "        df[label] = y_val[:,i]\n",
    "        df[label + '_pred'] = predict_val[:,i]\n",
    "\n",
    "        fpr[label], tpr[label], threshold = roc_curve(df[label],df[label+'_pred'])\n",
    "\n",
    "        auc1[label] = auc(fpr[label], tpr[label])\n",
    "\n",
    "        plt.plot(tpr[label],fpr[label],label='%s tagger, auc = %.1f%%'%(label,auc1[label]*100.))\n",
    "plt.semilogy()\n",
    "plt.xlabel(\"sig. efficiency\")\n",
    "plt.ylabel(\"bkg. mistag rate\")\n",
    "plt.ylim(0.000001,1)\n",
    "plt.grid(True)\n",
    "plt.legend(loc='lower right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "lzbQ-d0RKVmV"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "Notebook2_JetID_DNN.ipynb",
   "provenance": []
  },
  "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.8.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}