Week4/MakeROCfigure/MakeROCfigure.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ROC-curves\n",
    "\n",
    "Python script for illustrating ROC curves with Gaussian examples.\n",
    "\n",
    "For more information, see e.g. the Wikipedia page on __[ROC curves](https://en.wikipedia.org/wiki/Receiver_operating_characteristic)__. \n",
    "\n",
    "\n",
    "### Authors: \n",
    "- Christian Michelsen (Niels Bohr Institute)\n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 07-12-2022 (latest update)\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then set up variables for the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "r = np.random             # Random generator\n",
    "r.seed(42)                # Set a random seed (but a fixed one)\n",
    "\n",
    "save_plots = True        # For now, don't save plots (once you trust your code, switch on)\n",
    "N_entries  = 1_000_000    # Notice that Python (3) allows us to write undercores in numeric values without it changing the value. \n",
    "                          # Useful as a 1000 separater, in phone numbers or similar cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variables in the calculation of ROC-curves:\n",
    "\n",
    "For the computation of the ROC-curve, we rely heavily on the the Wikipedia page for __[ROC curves](https://en.wikipedia.org/wiki/Receiver_operating_characteristic)__ and especially the following figure:\n",
    "\n",
    "<img src=\"ROC_curves.png\" height=\"500\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, <span style=\"color:red\"> red is signal </span> and <span style=\"color:blue\">blue is background </span>. We have the following variables:\n",
    "\n",
    "- True Positives (__TP__, signal, power):  <span style=\"color:red\"> strong red in figure </span>\n",
    "- True Negatives (__TN__, background): <span style=\"color:blue\"> strong blue in the figure </span>\n",
    "- False Positives (__FP__, type 1 error):  <span style=\"color:IndianRed\"> light red in figure </span>\n",
    "- False Negatives (__FN__, type 2 error): <span style=\"color:MediumSlateBlue\"> light blue in figure </span>\n",
    "\n",
    "The __TPR__ is the ratio between TP and all true positives (TP, FN): $\\text{TPR}=\\frac{\\text{TP}}{\\text{TP} + \\text{FN}}$.\n",
    "\n",
    "The __FPR__ is the ratio between FP and all true negatives (FP, TN): $\\text{FPR}=\\frac{\\text{FP}}{\\text{FP} + \\text{TN}}$.\n",
    "\n",
    "The idea is that we move the cut (the vertical line, also known as threshold) from left to right and along the way measure the TPR and FPR as a function of the cut. Below we have made a function that, given two histograms, calculates the TPR and FPR as described: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "# Calculate ROC curve from two histograms (hist1 is signal, hist2 is background):\n",
    "def calc_ROC(hist1, hist2) :\n",
    "\n",
    "    # First we extract the entries (y values) and the edges of the histograms:\n",
    "    # Note how the \"_\" is simply used for the rest of what e.g. \"hist1\" returns (not really of our interest)\n",
    "    y_sig, x_sig_edges, _ = hist1 \n",
    "    y_bkg, x_bkg_edges, _ = hist2\n",
    "    \n",
    "    # Check that the two histograms have the same x edges:\n",
    "    if np.array_equal(x_sig_edges, x_bkg_edges) :\n",
    "        \n",
    "        # Extract the center positions (x values) of the bins (both signal or background works - equal binning)\n",
    "        x_centers = 0.5*(x_sig_edges[1:] + x_sig_edges[:-1])\n",
    "        \n",
    "        # Calculate the integral (sum) of the signal and background:\n",
    "        integral_sig = y_sig.sum()\n",
    "        integral_bkg = y_bkg.sum()\n",
    "    \n",
    "        # Initialize empty arrays for the True Positive Rate (TPR) and the False Positive Rate (FPR):\n",
    "        TPR = np.zeros_like(y_sig) # True positive rate (sensitivity)\n",
    "        FPR = np.zeros_like(y_sig) # False positive rate ()\n",
    "        \n",
    "        # Loop over all bins (x_centers) of the histograms and calculate TN, FP, FN, TP, FPR, and TPR for each bin:\n",
    "        for i, x in enumerate(x_centers): \n",
    "            \n",
    "            # The cut mask\n",
    "            cut = (x_centers < x)\n",
    "            \n",
    "            # True positive\n",
    "            TP = np.sum(y_sig[~cut]) / integral_sig    # True positives\n",
    "            FN = np.sum(y_sig[cut]) / integral_sig     # False negatives\n",
    "            TPR[i] = TP / (TP + FN)                    # True positive rate\n",
    "            \n",
    "            # True negative\n",
    "            TN = np.sum(y_bkg[cut]) / integral_bkg      # True negatives (background)\n",
    "            FP = np.sum(y_bkg[~cut]) / integral_bkg     # False positives\n",
    "            FPR[i] = FP / (FP + TN)                     # False positive rate            \n",
    "            \n",
    "        return FPR, TPR\n",
    "    \n",
    "    else:\n",
    "        AssertionError(\"Signal and Background histograms have different bins and/or ranges\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian distributions (signal & background) of varying separation:\n",
    "\n",
    "Now we generate the data that we are going to use. We want to generate data for signal and a background (both Gaussian) that are shifted unit Gaussians (i.e. mean different from 0) such that their separation is $N_\\sigma = $ `N_sigma`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_sigma = np.array([0.0, 0.5, 1.0, 1.5, 2.0, 3.0])      # We choose six signal-background separations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the above, define the means and widths of the signal and background Gaussians: We choose the widths to be unit, and the means to move symmetrically around $5.0$ (arbitrary choice):\n",
    "$$N_\\sigma = \\frac{\\mu_s - \\mu_b}{\\sqrt{ {\\sigma_s}^2+{\\sigma_b}^2 }} =  \\frac{\\Delta_\\mu}{\\sqrt 2}$$\n",
    "\n",
    "Which in code looks like:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_sig = 5 + N_sigma / np.sqrt(2.0)\n",
    "mean_bkg = 5 - N_sigma / np.sqrt(2.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having defined the means, we can now generate the random Gaussianly distributed numbers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "random_sig = r.normal(loc=mean_sig, scale=1, size=(N_entries, len(N_sigma)))\n",
    "random_bkg = r.normal(loc=mean_bkg, scale=1, size=(N_entries, len(N_sigma)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the data:\n",
    "\n",
    "Before we plot the data, we set up the figure consisting of 6 smaller plots (one for each value in `N_sigma`) and a large plot to show the ROC-curve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "from matplotlib.gridspec import GridSpec\n",
    "fig = plt.figure(figsize=(16*1.2,6*1.2))\n",
    "\n",
    "gs = GridSpec(2,5) # 2 rows, 3 columns\n",
    "ax = [fig.add_subplot(gs[0,0]), # First row,   first column\n",
    "      fig.add_subplot(gs[1,0]), # Second row,  first column\n",
    "      fig.add_subplot(gs[0,1]), # First row,   second column\n",
    "      fig.add_subplot(gs[1,1]), # Second row,  second column\n",
    "      fig.add_subplot(gs[0,2]), # First row,  third column\n",
    "      fig.add_subplot(gs[1,2])] # Second row,  third column\n",
    "ax2 = fig.add_subplot(gs[:,3:]) # Span both rows and columns from 3 and onwards"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we loop over all our data (`random_sig` and `random_bkg`) and plot it as histograms in the 6 small subplots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "N_bins = 200\n",
    "xmin, xmax = 0, 10\n",
    "\n",
    "hist_sig = []\n",
    "hist_bkg = []\n",
    "\n",
    "y_sig = np.zeros((N_bins, len(N_sigma)))\n",
    "y_bkg = np.zeros_like(y_sig)\n",
    "\n",
    "for i in range(len(N_sigma)):\n",
    "    hist_sig.append(ax[i].hist(random_sig[:, i], bins=N_bins, range=(xmin, xmax), histtype='step', lw=1, label='Signal', color='red'))\n",
    "    hist_bkg.append(ax[i].hist(random_bkg[:, i], bins=N_bins, range=(xmin, xmax), histtype='step', lw=1, label='Background', color='blue'))\n",
    "    ax[i].set_xlabel('Separating variable')\n",
    "    ax[i].set_ylabel('Frequency')\n",
    "    ax[i].text(0.75, 0.9, r'{:1.1f}$\\sigma$'.format(N_sigma[i]), size=16,\n",
    "              verticalalignment='center', transform=ax[i].transAxes)\n",
    "\n",
    "ax[0].legend(loc='upper left')\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having plotted the input distributions, we can now focus on the actual ROC-curves. We loop over all the histograms and calculate the FPR and TPR and plot it along with small text annotations describing the number of $\\sigma$ of separation: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Locations of the text placement on plot:\n",
    "text_sigma_loc = [(0.5, 0.45), (0.36, 0.59), (0.23, 0.72), (0.15, 0.82), (0.1, 0.9), (0.02, 0.96)]\n",
    "\n",
    "# Dictionary to save the FPRs and TPRs:\n",
    "d = {}\n",
    "\n",
    "for i in range(len(N_sigma)) :\n",
    "    \n",
    "    # Number of sigmas written as string ('r' to make it a raw string literal, i.e. allows LaTeX-style):\n",
    "    text_sigma = fr'{N_sigma[i]:1.1f}$\\sigma$'\n",
    "    \n",
    "    # Calculate FPR and TPR by calc_ROC:\n",
    "    FPR, TPR = calc_ROC(hist_sig[i], hist_bkg[i])    \n",
    "    d[int(N_sigma[i]*10)] = {'FPR': FPR, 'TPR': TPR}\n",
    "    \n",
    "    # Plot the roc-curve:\n",
    "    ax2.plot(FPR, TPR, '-', label=text_sigma)\n",
    "    \n",
    "    # Add text:\n",
    "    ax2.text(text_sigma_loc[i][0], text_sigma_loc[i][1], text_sigma, size=12, verticalalignment='center', transform=ax2.transAxes)\n",
    "    \n",
    "\n",
    "# Plot diagonal line:\n",
    "ax2.plot([0, 1], [0, 1], '--', label='Random Guess')\n",
    "\n",
    "# Set figure options:\n",
    "ax2.set(xlabel='False Positive Rate (Background Efficiency)', ylabel='True Positive Rate (Signal Efficiency)', xlim=(0, 1), ylim=(0, 1))\n",
    "ax2.legend()\n",
    "\n",
    "fig.tight_layout()\n",
    "fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the plots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "if save_plots:\n",
    "    fig.savefig('ROCcurves_GaussianSeparations.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple alternative:\n",
    "\n",
    "Scikit-learn is a package for Machine Learning (ML), and Scikit-plot is package for plotting in relation to ML. Therefore, it has a build-in ROC curve plotting tool, which is shown below. It requires two arrays, one with the true labels (which we often denote \"y\" or \"y_true\"), that is the target of the learning (i.e. typically two cases: 0 and 1), and one with the predictions of an (ML-) algorithm (\"y_hat\" or \"y_predict\"), which are floating numbers attempted - for the two cases - to be as separated from each other, as possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import scikitplot as skplt\n",
    "# \n",
    "# y_true = # ground truth labels\n",
    "# y_probas = # predicted probabilities generated by sklearn classifier\n",
    "# skplt.metrics.plot_roc_curve(y_true, y_probas)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Below we show how to save files using Pickle. This is saved as a binary file and then read in again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "filename = f'ROC_curve_FPR_TPR_N_{N_entries}.pkl'\n",
    "\n",
    "with open(filename, 'wb') as outfile:\n",
    "    pickle.dump(d, outfile)\n",
    "\n",
    "with open(filename, 'rb') as infile:\n",
    "    d_new = pickle.load(infile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# \"Question\"\n",
    "- Make sure that you understand what the ROC represents, and how it is calculated!"
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}