HEPML_HandsOn_NN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "O3V3tNoJVYDF"
   },
   "source": [
    "# Hands on : introduction to NN on HEP dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pfxlugZJVlR0"
   },
   "source": [
    "### Many thanks to David Rousseau, Yann Coadou, and Aishik Gosh. This tutorial is based on their work!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eT7-MMpfrlHR"
   },
   "source": [
    "## Import Packages\n",
    "\n",
    "[This should look familiar!]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "FO_W4IvbVYDL"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import uproot as ur\n",
    "\n",
    "import tensorflow as tf\n",
    "os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'\n",
    "\n",
    "from IPython.display import display, HTML\n",
    "%matplotlib inline\n",
    "import time\n",
    "pd.set_option('display.max_columns', None) # to see all columns of df.head()\n",
    "np.random.seed(31415) # set the random seed for the reproducibility"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.4.1\n"
     ]
    }
   ],
   "source": [
    "print(tf.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Num GPUs Available:  0\n"
     ]
    }
   ],
   "source": [
    "print(\"Num GPUs Available: \", len(tf.config.list_physical_devices('GPU')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dHSzcuTgVYDT",
    "tags": []
   },
   "source": [
    "# Load events"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "data was created from ATLAS Open Data see doc\n",
    "\n",
    "http://opendata.atlas.cern/release/2020/documentation/datasets/intro.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kl-4W4Ifs8EV",
    "outputId": "54d366bc-8925-4c6a-e867-ef70016e016b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'tree_event;1': 'TTree'}\n"
     ]
    }
   ],
   "source": [
    "filename=(\"dataWW_d1.root\")\n",
    "file = ur.open(filename)\n",
    "print(file.classnames())\n",
    "#print(file.keys())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['index', 'eventNumber', 'label', 'met_et', 'met_phi', 'lep_n', 'lep_pt_0', 'lep_pt_1', 'lep_eta_0', 'lep_eta_1', 'lep_phi_0', 'lep_phi_1', 'lep_E_0', 'lep_E_1', 'lep_charge_0', 'lep_charge_1', 'lep_type_0', 'lep_type_1', 'jet_n', 'jet_pt_0', 'jet_pt_1', 'jet_eta_0', 'jet_eta_1', 'jet_phi_0', 'jet_phi_1', 'jet_E_0', 'jet_E_1', 'mcWeight', 'runNumber', 'channelNumber']\n"
     ]
    }
   ],
   "source": [
    "tree = file[\"tree_event\"]\n",
    "print(tree.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'uproot.models.TTree.Model_TTree_v20'>\n",
      "name                 | typename                 | interpretation                \n",
      "---------------------+--------------------------+-------------------------------\n",
      "index                | int64_t                  | AsDtype('>i8')\n",
      "eventNumber          | int64_t                  | AsDtype('>i8')\n",
      "label                | int64_t                  | AsDtype('>i8')\n",
      "met_et               | double                   | AsDtype('>f8')\n",
      "met_phi              | double                   | AsDtype('>f8')\n",
      "lep_n                | int64_t                  | AsDtype('>i8')\n",
      "lep_pt_0             | double                   | AsDtype('>f8')\n",
      "lep_pt_1             | double                   | AsDtype('>f8')\n",
      "lep_eta_0            | double                   | AsDtype('>f8')\n",
      "lep_eta_1            | double                   | AsDtype('>f8')\n",
      "lep_phi_0            | double                   | AsDtype('>f8')\n",
      "lep_phi_1            | double                   | AsDtype('>f8')\n",
      "lep_E_0              | double                   | AsDtype('>f8')\n",
      "lep_E_1              | double                   | AsDtype('>f8')\n",
      "lep_charge_0         | int64_t                  | AsDtype('>i8')\n",
      "lep_charge_1         | int64_t                  | AsDtype('>i8')\n",
      "lep_type_0           | int64_t                  | AsDtype('>i8')\n",
      "lep_type_1           | int64_t                  | AsDtype('>i8')\n",
      "jet_n                | int64_t                  | AsDtype('>i8')\n",
      "jet_pt_0             | double                   | AsDtype('>f8')\n",
      "jet_pt_1             | double                   | AsDtype('>f8')\n",
      "jet_eta_0            | double                   | AsDtype('>f8')\n",
      "jet_eta_1            | double                   | AsDtype('>f8')\n",
      "jet_phi_0            | double                   | AsDtype('>f8')\n",
      "jet_phi_1            | double                   | AsDtype('>f8')\n",
      "jet_E_0              | double                   | AsDtype('>f8')\n",
      "jet_E_1              | double                   | AsDtype('>f8')\n",
      "mcWeight             | double                   | AsDtype('>f8')\n",
      "runNumber            | int64_t                  | AsDtype('>i8')\n",
      "channelNumber        | int64_t                  | AsDtype('>i8')\n"
     ]
    }
   ],
   "source": [
    "print(type(tree))\n",
    "tree.show()\n",
    "dfall = tree.arrays(library=\"pd\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[uproot](https://indico.cern.ch/event/686641/contributions/2894906/attachments/1606247/2548596/pivarski-uproot.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#csv_file = \"/global/cfs/cdirs/ntrain5/atlas-ml-training/dnn_bdt_tutorial/dataWW_d1.csv\"\n",
    "#dfall = pd.read_csv(csv_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File loaded with  600000  events \n"
     ]
    }
   ],
   "source": [
    "#shuffle the events [already done but just to be safe]\n",
    "dfall = dfall.sample(frac=1).reset_index(drop=True)\n",
    "print (\"File loaded with \",dfall.shape[0], \" events \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7xYM7O_zVYDY"
   },
   "source": [
    "#### At this point, you should see \"File Loaded with [XXX] events\". If not, the data file could not be accessed. No point going further!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LFPm8OHuVYDY"
   },
   "source": [
    "# Examine Pandas Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "xNeCAIuFVYDZ",
    "outputId": "b12cbd30-536e-429c-a1be-677ee0b937fb"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['index', 'eventNumber', 'label', 'met_et', 'met_phi', 'lep_n',\n",
       "       'lep_pt_0', 'lep_pt_1', 'lep_eta_0', 'lep_eta_1', 'lep_phi_0',\n",
       "       'lep_phi_1', 'lep_E_0', 'lep_E_1', 'lep_charge_0', 'lep_charge_1',\n",
       "       'lep_type_0', 'lep_type_1', 'jet_n', 'jet_pt_0', 'jet_pt_1',\n",
       "       'jet_eta_0', 'jet_eta_1', 'jet_phi_0', 'jet_phi_1', 'jet_E_0',\n",
       "       'jet_E_1', 'mcWeight', 'runNumber', 'channelNumber'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#dump list of feature\n",
    "dfall.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 226
    },
    "id": "j9l7wkkorlHe",
    "outputId": "7d824317-2287-44b0-e153-112ecf7ddc13"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>eventNumber</th>\n",
       "      <th>label</th>\n",
       "      <th>met_et</th>\n",
       "      <th>met_phi</th>\n",
       "      <th>lep_n</th>\n",
       "      <th>lep_pt_0</th>\n",
       "      <th>lep_pt_1</th>\n",
       "      <th>lep_eta_0</th>\n",
       "      <th>lep_eta_1</th>\n",
       "      <th>lep_phi_0</th>\n",
       "      <th>lep_phi_1</th>\n",
       "      <th>lep_E_0</th>\n",
       "      <th>lep_E_1</th>\n",
       "      <th>lep_charge_0</th>\n",
       "      <th>lep_charge_1</th>\n",
       "      <th>lep_type_0</th>\n",
       "      <th>lep_type_1</th>\n",
       "      <th>jet_n</th>\n",
       "      <th>jet_pt_0</th>\n",
       "      <th>jet_pt_1</th>\n",
       "      <th>jet_eta_0</th>\n",
       "      <th>jet_eta_1</th>\n",
       "      <th>jet_phi_0</th>\n",
       "      <th>jet_phi_1</th>\n",
       "      <th>jet_E_0</th>\n",
       "      <th>jet_E_1</th>\n",
       "      <th>mcWeight</th>\n",
       "      <th>runNumber</th>\n",
       "      <th>channelNumber</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>543448</td>\n",
       "      <td>402756</td>\n",
       "      <td>1</td>\n",
       "      <td>25.609</td>\n",
       "      <td>0.42452</td>\n",
       "      <td>2</td>\n",
       "      <td>48.295</td>\n",
       "      <td>15.214</td>\n",
       "      <td>0.73991</td>\n",
       "      <td>2.27420</td>\n",
       "      <td>-2.316400</td>\n",
       "      <td>-1.39410</td>\n",
       "      <td>62129.0</td>\n",
       "      <td>74721.0</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>11</td>\n",
       "      <td>2</td>\n",
       "      <td>26.32</td>\n",
       "      <td>20.064</td>\n",
       "      <td>-1.1350</td>\n",
       "      <td>-2.07540</td>\n",
       "      <td>-2.8092</td>\n",
       "      <td>2.0740</td>\n",
       "      <td>45397.0</td>\n",
       "      <td>81304.0</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>284500</td>\n",
       "      <td>345323</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>580260</td>\n",
       "      <td>101274</td>\n",
       "      <td>0</td>\n",
       "      <td>196.560</td>\n",
       "      <td>1.31140</td>\n",
       "      <td>2</td>\n",
       "      <td>69.459</td>\n",
       "      <td>21.081</td>\n",
       "      <td>-0.52666</td>\n",
       "      <td>0.22380</td>\n",
       "      <td>0.023132</td>\n",
       "      <td>-0.67855</td>\n",
       "      <td>79317.0</td>\n",
       "      <td>21611.0</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>13</td>\n",
       "      <td>2</td>\n",
       "      <td>192.12</td>\n",
       "      <td>36.217</td>\n",
       "      <td>-1.2643</td>\n",
       "      <td>-0.69746</td>\n",
       "      <td>-2.2147</td>\n",
       "      <td>-1.9152</td>\n",
       "      <td>367670.0</td>\n",
       "      <td>45733.0</td>\n",
       "      <td>0.000288</td>\n",
       "      <td>284500</td>\n",
       "      <td>363492</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>112856</td>\n",
       "      <td>468437</td>\n",
       "      <td>1</td>\n",
       "      <td>45.653</td>\n",
       "      <td>-2.76860</td>\n",
       "      <td>2</td>\n",
       "      <td>45.927</td>\n",
       "      <td>22.822</td>\n",
       "      <td>-1.61910</td>\n",
       "      <td>-2.00770</td>\n",
       "      <td>1.901700</td>\n",
       "      <td>-0.11248</td>\n",
       "      <td>120480.0</td>\n",
       "      <td>86498.0</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "      <td>0</td>\n",
       "      <td>-7.00</td>\n",
       "      <td>-7.000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.00000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>0.000029</td>\n",
       "      <td>284500</td>\n",
       "      <td>345324</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>121430</td>\n",
       "      <td>272337</td>\n",
       "      <td>1</td>\n",
       "      <td>49.415</td>\n",
       "      <td>-0.57805</td>\n",
       "      <td>2</td>\n",
       "      <td>45.929</td>\n",
       "      <td>14.263</td>\n",
       "      <td>1.55280</td>\n",
       "      <td>0.12809</td>\n",
       "      <td>2.254700</td>\n",
       "      <td>-2.54810</td>\n",
       "      <td>113360.0</td>\n",
       "      <td>14380.0</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "      <td>13</td>\n",
       "      <td>13</td>\n",
       "      <td>0</td>\n",
       "      <td>-7.00</td>\n",
       "      <td>-7.000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.00000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>0.000029</td>\n",
       "      <td>284500</td>\n",
       "      <td>345324</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>55912</td>\n",
       "      <td>354546</td>\n",
       "      <td>1</td>\n",
       "      <td>71.988</td>\n",
       "      <td>-2.60390</td>\n",
       "      <td>2</td>\n",
       "      <td>62.029</td>\n",
       "      <td>21.453</td>\n",
       "      <td>-0.51082</td>\n",
       "      <td>-0.38177</td>\n",
       "      <td>0.236200</td>\n",
       "      <td>1.40820</td>\n",
       "      <td>70300.0</td>\n",
       "      <td>23036.0</td>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "      <td>11</td>\n",
       "      <td>13</td>\n",
       "      <td>0</td>\n",
       "      <td>-7.00</td>\n",
       "      <td>-7.000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.00000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.0000</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>284500</td>\n",
       "      <td>345323</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    index  eventNumber  label   met_et  met_phi  lep_n  lep_pt_0  lep_pt_1  \\\n",
       "0  543448       402756      1   25.609  0.42452      2    48.295    15.214   \n",
       "1  580260       101274      0  196.560  1.31140      2    69.459    21.081   \n",
       "2  112856       468437      1   45.653 -2.76860      2    45.927    22.822   \n",
       "3  121430       272337      1   49.415 -0.57805      2    45.929    14.263   \n",
       "4   55912       354546      1   71.988 -2.60390      2    62.029    21.453   \n",
       "\n",
       "   lep_eta_0  lep_eta_1  lep_phi_0  lep_phi_1   lep_E_0  lep_E_1  \\\n",
       "0    0.73991    2.27420  -2.316400   -1.39410   62129.0  74721.0   \n",
       "1   -0.52666    0.22380   0.023132   -0.67855   79317.0  21611.0   \n",
       "2   -1.61910   -2.00770   1.901700   -0.11248  120480.0  86498.0   \n",
       "3    1.55280    0.12809   2.254700   -2.54810  113360.0  14380.0   \n",
       "4   -0.51082   -0.38177   0.236200    1.40820   70300.0  23036.0   \n",
       "\n",
       "   lep_charge_0  lep_charge_1  lep_type_0  lep_type_1  jet_n  jet_pt_0  \\\n",
       "0            -1             1          13          11      2     26.32   \n",
       "1            -1             1          13          13      2    192.12   \n",
       "2            -1             1          11          11      0     -7.00   \n",
       "3            -1             1          13          13      0     -7.00   \n",
       "4            -1             1          11          13      0     -7.00   \n",
       "\n",
       "   jet_pt_1  jet_eta_0  jet_eta_1  jet_phi_0  jet_phi_1   jet_E_0  jet_E_1  \\\n",
       "0    20.064    -1.1350   -2.07540    -2.8092     2.0740   45397.0  81304.0   \n",
       "1    36.217    -1.2643   -0.69746    -2.2147    -1.9152  367670.0  45733.0   \n",
       "2    -7.000    -7.0000   -7.00000    -7.0000    -7.0000      -7.0     -7.0   \n",
       "3    -7.000    -7.0000   -7.00000    -7.0000    -7.0000      -7.0     -7.0   \n",
       "4    -7.000    -7.0000   -7.00000    -7.0000    -7.0000      -7.0     -7.0   \n",
       "\n",
       "   mcWeight  runNumber  channelNumber  \n",
       "0  0.000002     284500         345323  \n",
       "1  0.000288     284500         363492  \n",
       "2  0.000029     284500         345324  \n",
       "3  0.000029     284500         345324  \n",
       "4  0.000002     284500         345323  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#examine first few events\n",
    "dfall.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 320
    },
    "id": "Oz-lWJhgrlHg",
    "outputId": "7ad3adfa-9819-47cb-c321-6d0d884c96cf"
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>eventNumber</th>\n",
       "      <th>label</th>\n",
       "      <th>met_et</th>\n",
       "      <th>met_phi</th>\n",
       "      <th>lep_n</th>\n",
       "      <th>lep_pt_0</th>\n",
       "      <th>lep_pt_1</th>\n",
       "      <th>lep_eta_0</th>\n",
       "      <th>lep_eta_1</th>\n",
       "      <th>lep_phi_0</th>\n",
       "      <th>lep_phi_1</th>\n",
       "      <th>lep_E_0</th>\n",
       "      <th>lep_E_1</th>\n",
       "      <th>lep_charge_0</th>\n",
       "      <th>lep_charge_1</th>\n",
       "      <th>lep_type_0</th>\n",
       "      <th>lep_type_1</th>\n",
       "      <th>jet_n</th>\n",
       "      <th>jet_pt_0</th>\n",
       "      <th>jet_pt_1</th>\n",
       "      <th>jet_eta_0</th>\n",
       "      <th>jet_eta_1</th>\n",
       "      <th>jet_phi_0</th>\n",
       "      <th>jet_phi_1</th>\n",
       "      <th>jet_E_0</th>\n",
       "      <th>jet_E_1</th>\n",
       "      <th>mcWeight</th>\n",
       "      <th>runNumber</th>\n",
       "      <th>channelNumber</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>600000.000000</td>\n",
       "      <td>6.000000e+05</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>6.000000e+05</td>\n",
       "      <td>6.000000e+05</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>6.000000e+05</td>\n",
       "      <td>6.000000e+05</td>\n",
       "      <td>600000.000000</td>\n",
       "      <td>600000.0</td>\n",
       "      <td>600000.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>299999.500000</td>\n",
       "      <td>7.029775e+05</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>60.365356</td>\n",
       "      <td>-0.008556</td>\n",
       "      <td>2.005748</td>\n",
       "      <td>60.270328</td>\n",
       "      <td>28.133970</td>\n",
       "      <td>0.000269</td>\n",
       "      <td>-0.000664</td>\n",
       "      <td>0.005219</td>\n",
       "      <td>0.008524</td>\n",
       "      <td>1.164008e+05</td>\n",
       "      <td>5.630948e+04</td>\n",
       "      <td>-0.027027</td>\n",
       "      <td>0.027037</td>\n",
       "      <td>11.907273</td>\n",
       "      <td>11.958010</td>\n",
       "      <td>1.334185</td>\n",
       "      <td>59.606114</td>\n",
       "      <td>13.521164</td>\n",
       "      <td>-1.934521</td>\n",
       "      <td>-4.402337</td>\n",
       "      <td>-1.939645</td>\n",
       "      <td>-4.401941</td>\n",
       "      <td>1.375410e+05</td>\n",
       "      <td>4.125468e+04</td>\n",
       "      <td>0.000218</td>\n",
       "      <td>284500.0</td>\n",
       "      <td>351379.666667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>173205.225094</td>\n",
       "      <td>4.568961e+05</td>\n",
       "      <td>0.471405</td>\n",
       "      <td>48.153210</td>\n",
       "      <td>1.812778</td>\n",
       "      <td>0.076061</td>\n",
       "      <td>42.676858</td>\n",
       "      <td>20.352272</td>\n",
       "      <td>1.210614</td>\n",
       "      <td>1.235890</td>\n",
       "      <td>1.811177</td>\n",
       "      <td>1.814205</td>\n",
       "      <td>1.083341e+05</td>\n",
       "      <td>5.748834e+04</td>\n",
       "      <td>0.999636</td>\n",
       "      <td>0.999635</td>\n",
       "      <td>0.995692</td>\n",
       "      <td>0.999119</td>\n",
       "      <td>1.232632</td>\n",
       "      <td>79.586143</td>\n",
       "      <td>38.074439</td>\n",
       "      <td>3.349404</td>\n",
       "      <td>3.491304</td>\n",
       "      <td>3.489197</td>\n",
       "      <td>3.559421</td>\n",
       "      <td>1.980717e+05</td>\n",
       "      <td>8.990068e+04</td>\n",
       "      <td>0.000736</td>\n",
       "      <td>0.0</td>\n",
       "      <td>8564.720183</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.000000e+00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.044211</td>\n",
       "      <td>-3.141600</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>25.000000</td>\n",
       "      <td>7.000300</td>\n",
       "      <td>-2.700000</td>\n",
       "      <td>-2.699900</td>\n",
       "      <td>-3.141600</td>\n",
       "      <td>-3.141600</td>\n",
       "      <td>2.501300e+04</td>\n",
       "      <td>7.003800e+03</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000e+00</td>\n",
       "      <td>-7.000000e+00</td>\n",
       "      <td>-0.045557</td>\n",
       "      <td>284500.0</td>\n",
       "      <td>345323.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>149999.750000</td>\n",
       "      <td>3.498568e+05</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>33.196000</td>\n",
       "      <td>-1.577600</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>36.624750</td>\n",
       "      <td>15.316000</td>\n",
       "      <td>-0.899650</td>\n",
       "      <td>-0.946043</td>\n",
       "      <td>-1.560200</td>\n",
       "      <td>-1.562900</td>\n",
       "      <td>5.221675e+04</td>\n",
       "      <td>2.329100e+04</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-7.000000e+00</td>\n",
       "      <td>-7.000000e+00</td>\n",
       "      <td>0.000002</td>\n",
       "      <td>284500.0</td>\n",
       "      <td>345323.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>299999.500000</td>\n",
       "      <td>6.678815e+05</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>50.082000</td>\n",
       "      <td>-0.018694</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>48.229500</td>\n",
       "      <td>23.644000</td>\n",
       "      <td>0.003413</td>\n",
       "      <td>0.002974</td>\n",
       "      <td>0.014583</td>\n",
       "      <td>0.018823</td>\n",
       "      <td>8.184450e+04</td>\n",
       "      <td>3.824900e+04</td>\n",
       "      <td>-1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>11.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>38.744000</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-0.891010</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>-1.210600</td>\n",
       "      <td>-7.000000</td>\n",
       "      <td>7.501700e+04</td>\n",
       "      <td>-7.000000e+00</td>\n",
       "      <td>0.000029</td>\n",
       "      <td>284500.0</td>\n",
       "      <td>345324.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>449999.250000</td>\n",
       "      <td>1.001273e+06</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>73.032250</td>\n",
       "      <td>1.559300</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>68.979000</td>\n",
       "      <td>34.149000</td>\n",
       "      <td>0.901950</td>\n",
       "      <td>0.942210</td>\n",
       "      <td>1.573000</td>\n",
       "      <td>1.577300</td>\n",
       "      <td>1.412400e+05</td>\n",
       "      <td>6.730825e+04</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>13.000000</td>\n",
       "      <td>13.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>84.291250</td>\n",
       "      <td>26.881000</td>\n",
       "      <td>0.724705</td>\n",
       "      <td>-0.821005</td>\n",
       "      <td>0.967753</td>\n",
       "      <td>-1.088575</td>\n",
       "      <td>1.766200e+05</td>\n",
       "      <td>5.060500e+04</td>\n",
       "      <td>0.000197</td>\n",
       "      <td>284500.0</td>\n",
       "      <td>363492.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>599999.000000</td>\n",
       "      <td>1.997995e+06</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>7062.200000</td>\n",
       "      <td>3.141600</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>7041.400000</td>\n",
       "      <td>744.900000</td>\n",
       "      <td>2.699800</td>\n",
       "      <td>2.700000</td>\n",
       "      <td>3.141500</td>\n",
       "      <td>3.141600</td>\n",
       "      <td>7.041600e+06</td>\n",
       "      <td>2.107900e+06</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>13.000000</td>\n",
       "      <td>13.000000</td>\n",
       "      <td>13.000000</td>\n",
       "      <td>1658.200000</td>\n",
       "      <td>1429.800000</td>\n",
       "      <td>2.500000</td>\n",
       "      <td>2.500000</td>\n",
       "      <td>3.141600</td>\n",
       "      <td>3.141600</td>\n",
       "      <td>5.638800e+06</td>\n",
       "      <td>2.784300e+06</td>\n",
       "      <td>0.052420</td>\n",
       "      <td>284500.0</td>\n",
       "      <td>363492.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               index   eventNumber          label         met_et  \\\n",
       "count  600000.000000  6.000000e+05  600000.000000  600000.000000   \n",
       "mean   299999.500000  7.029775e+05       0.666667      60.365356   \n",
       "std    173205.225094  4.568961e+05       0.471405      48.153210   \n",
       "min         0.000000  1.000000e+00       0.000000       0.044211   \n",
       "25%    149999.750000  3.498568e+05       0.000000      33.196000   \n",
       "50%    299999.500000  6.678815e+05       1.000000      50.082000   \n",
       "75%    449999.250000  1.001273e+06       1.000000      73.032250   \n",
       "max    599999.000000  1.997995e+06       1.000000    7062.200000   \n",
       "\n",
       "             met_phi          lep_n       lep_pt_0       lep_pt_1  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean       -0.008556       2.005748      60.270328      28.133970   \n",
       "std         1.812778       0.076061      42.676858      20.352272   \n",
       "min        -3.141600       2.000000      25.000000       7.000300   \n",
       "25%        -1.577600       2.000000      36.624750      15.316000   \n",
       "50%        -0.018694       2.000000      48.229500      23.644000   \n",
       "75%         1.559300       2.000000      68.979000      34.149000   \n",
       "max         3.141600       5.000000    7041.400000     744.900000   \n",
       "\n",
       "           lep_eta_0      lep_eta_1      lep_phi_0      lep_phi_1  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean        0.000269      -0.000664       0.005219       0.008524   \n",
       "std         1.210614       1.235890       1.811177       1.814205   \n",
       "min        -2.700000      -2.699900      -3.141600      -3.141600   \n",
       "25%        -0.899650      -0.946043      -1.560200      -1.562900   \n",
       "50%         0.003413       0.002974       0.014583       0.018823   \n",
       "75%         0.901950       0.942210       1.573000       1.577300   \n",
       "max         2.699800       2.700000       3.141500       3.141600   \n",
       "\n",
       "            lep_E_0       lep_E_1   lep_charge_0   lep_charge_1  \\\n",
       "count  6.000000e+05  6.000000e+05  600000.000000  600000.000000   \n",
       "mean   1.164008e+05  5.630948e+04      -0.027027       0.027037   \n",
       "std    1.083341e+05  5.748834e+04       0.999636       0.999635   \n",
       "min    2.501300e+04  7.003800e+03      -1.000000      -1.000000   \n",
       "25%    5.221675e+04  2.329100e+04      -1.000000      -1.000000   \n",
       "50%    8.184450e+04  3.824900e+04      -1.000000       1.000000   \n",
       "75%    1.412400e+05  6.730825e+04       1.000000       1.000000   \n",
       "max    7.041600e+06  2.107900e+06       1.000000       1.000000   \n",
       "\n",
       "          lep_type_0     lep_type_1          jet_n       jet_pt_0  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean       11.907273      11.958010       1.334185      59.606114   \n",
       "std         0.995692       0.999119       1.232632      79.586143   \n",
       "min        11.000000      11.000000       0.000000      -7.000000   \n",
       "25%        11.000000      11.000000       0.000000      -7.000000   \n",
       "50%        11.000000      11.000000       1.000000      38.744000   \n",
       "75%        13.000000      13.000000       2.000000      84.291250   \n",
       "max        13.000000      13.000000      13.000000    1658.200000   \n",
       "\n",
       "            jet_pt_1      jet_eta_0      jet_eta_1      jet_phi_0  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean       13.521164      -1.934521      -4.402337      -1.939645   \n",
       "std        38.074439       3.349404       3.491304       3.489197   \n",
       "min        -7.000000      -7.000000      -7.000000      -7.000000   \n",
       "25%        -7.000000      -7.000000      -7.000000      -7.000000   \n",
       "50%        -7.000000      -0.891010      -7.000000      -1.210600   \n",
       "75%        26.881000       0.724705      -0.821005       0.967753   \n",
       "max      1429.800000       2.500000       2.500000       3.141600   \n",
       "\n",
       "           jet_phi_1       jet_E_0       jet_E_1       mcWeight  runNumber  \\\n",
       "count  600000.000000  6.000000e+05  6.000000e+05  600000.000000   600000.0   \n",
       "mean       -4.401941  1.375410e+05  4.125468e+04       0.000218   284500.0   \n",
       "std         3.559421  1.980717e+05  8.990068e+04       0.000736        0.0   \n",
       "min        -7.000000 -7.000000e+00 -7.000000e+00      -0.045557   284500.0   \n",
       "25%        -7.000000 -7.000000e+00 -7.000000e+00       0.000002   284500.0   \n",
       "50%        -7.000000  7.501700e+04 -7.000000e+00       0.000029   284500.0   \n",
       "75%        -1.088575  1.766200e+05  5.060500e+04       0.000197   284500.0   \n",
       "max         3.141600  5.638800e+06  2.784300e+06       0.052420   284500.0   \n",
       "\n",
       "       channelNumber  \n",
       "count  600000.000000  \n",
       "mean   351379.666667  \n",
       "std      8564.720183  \n",
       "min    345323.000000  \n",
       "25%    345323.000000  \n",
       "50%    345324.000000  \n",
       "75%    363492.000000  \n",
       "max    363492.000000  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#take a look at feature distribution\n",
    "dfall.describe()\n",
    "#dfall.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "lQsalTmorlHj",
    "outputId": "58b3c0c9-93c6-4757-cf69-66278940e60f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "total label weights (124.6622944642, 6.2000533976)\n",
      "total class number of events (200000, 400000)\n"
     ]
    }
   ],
   "source": [
    "label_weights = (dfall[dfall.label==0].mcWeight.sum(), dfall[dfall.label==1].mcWeight.sum() ) \n",
    "print(\"total label weights\",label_weights)\n",
    "\n",
    "label_nevents = (dfall[dfall.label==0].shape[0], dfall[dfall.label==1].shape[0] )\n",
    "print (\"total class number of events\",label_nevents)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rtI5u5GErlHq"
   },
   "source": [
    "## Event Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook essentially tries to classify events containing a Higgs Boson.\n",
    "\n",
    "The simulation includes top-quark-pair production, single-top production, production of weak bosons in association with jets (W+jets, Z+jets), production of a pair of bosons (diboson WW, WZ, ZZ) and __SM Higgs__ production.\n",
    "\n",
    "We will only keep events with exactly two leptons dfall.lep_n==2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kaO2JM1hrlHr",
    "outputId": "bc59e98d-f6f1-42b7-e9f3-48052fe4ebc9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(600000, 30)\n",
      "(577357, 30)\n"
     ]
    }
   ],
   "source": [
    "print (dfall.shape)\n",
    "\n",
    "# Also only keep events with positive weight. This is in principle wrong. \n",
    "#Many Data Science tools break given a negative weight.\n",
    "\n",
    "fulldata=dfall[(dfall.lep_n==2) & (dfall.mcWeight>0 )] # only keep events with exactly two leptons \n",
    "print (fulldata.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iMQpKhDKrlH0",
    "outputId": "ccb36177-999c-4203-939e-4e04383c1954"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(577357, 28)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#hide label and weights in separate vectors\n",
    "#they are not real features\n",
    "\n",
    "#WARNING : there should be no selection nor shuffling later on !\n",
    "target = fulldata[\"label\"]\n",
    "del fulldata[\"label\"]\n",
    "\n",
    "#hide weight in separate vector\n",
    "weights = fulldata[\"mcWeight\"]\n",
    "del fulldata[\"mcWeight\"]\n",
    "fulldata.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nviyIMgerlH3"
   },
   "source": [
    "\n",
    "# Try not to change the cells above $\\uparrow$\n",
    "...and return to this cell (or rerun the whole notebook) after changing things below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For simplicity, we'll only use some features on the first pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 225
    },
    "id": "6e0Hlpv6rlH4",
    "outputId": "8705b9a0-9dcc-4eb6-ab49-b3fd0b8ea1a6",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(577357, 6)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>met_et</th>\n",
       "      <th>met_phi</th>\n",
       "      <th>lep_pt_0</th>\n",
       "      <th>lep_pt_1</th>\n",
       "      <th>lep_phi_0</th>\n",
       "      <th>lep_phi_1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>25.609</td>\n",
       "      <td>0.42452</td>\n",
       "      <td>48.295</td>\n",
       "      <td>15.214</td>\n",
       "      <td>-2.316400</td>\n",
       "      <td>-1.39410</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>196.560</td>\n",
       "      <td>1.31140</td>\n",
       "      <td>69.459</td>\n",
       "      <td>21.081</td>\n",
       "      <td>0.023132</td>\n",
       "      <td>-0.67855</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>45.653</td>\n",
       "      <td>-2.76860</td>\n",
       "      <td>45.927</td>\n",
       "      <td>22.822</td>\n",
       "      <td>1.901700</td>\n",
       "      <td>-0.11248</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>49.415</td>\n",
       "      <td>-0.57805</td>\n",
       "      <td>45.929</td>\n",
       "      <td>14.263</td>\n",
       "      <td>2.254700</td>\n",
       "      <td>-2.54810</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>71.988</td>\n",
       "      <td>-2.60390</td>\n",
       "      <td>62.029</td>\n",
       "      <td>21.453</td>\n",
       "      <td>0.236200</td>\n",
       "      <td>1.40820</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    met_et  met_phi  lep_pt_0  lep_pt_1  lep_phi_0  lep_phi_1\n",
       "0   25.609  0.42452    48.295    15.214  -2.316400   -1.39410\n",
       "1  196.560  1.31140    69.459    21.081   0.023132   -0.67855\n",
       "2   45.653 -2.76860    45.927    22.822   1.901700   -0.11248\n",
       "3   49.415 -0.57805    45.929    14.263   2.254700   -2.54810\n",
       "4   71.988 -2.60390    62.029    21.453   0.236200    1.40820"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data=pd.DataFrame(fulldata, columns=[\"met_et\",\"met_phi\",\"lep_pt_0\",\"lep_pt_1\",'lep_phi_0', 'lep_phi_1'])\n",
    "print (data.shape)\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Krg3pTRjVYDs",
    "tags": []
   },
   "source": [
    "### Feature engineering (Two variations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. See if using more features improves model performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "more_features = False\n",
    "if (more_features):\n",
    "    data=pd.DataFrame(fulldata, columns=[\"met_et\",\"met_phi\",\"lep_pt_0\",\"lep_pt_1\",'lep_eta_0', 'lep_eta_1', 'lep_phi_0', 'lep_phi_1','jet_n','jet_pt_0',\n",
    "       'jet_pt_1', 'jet_eta_0', 'jet_eta_1', 'jet_phi_0', 'jet_phi_1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Engineer our own feature, $\\Delta\\varphi_l$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "id": "DH5UqBF9VYDt"
   },
   "outputs": [],
   "source": [
    "use_deltaphi = False\n",
    "if use_deltaphi: \n",
    "    data[\"lep_deltaphi\"]=np.abs(np.mod(data.lep_phi_1-data.lep_phi_0+3*np.pi,2*np.pi)-np.pi)\n",
    "\n",
    "    print (data.shape)\n",
    "    display(data.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 734
    },
    "id": "o2mf1bLVrlH7",
    "outputId": "b4484ca0-0983-4477-9af7-26bf496aa618"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3cAAAK7CAYAAABRbnZtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABlLklEQVR4nO39f7xld13ffb/eTn6gJLlDbnAMk9SEdhycWsR0msRSvY4imomWAavXlZlCMNh7TE0EWqkd5LKOF/W+qSJKau7MNUAsKYQU+aFTHBsourU+boOTxBAYwugYgRkyEBUJHNISJvncf+w1YefknH32OXufs/de+/V8PPZj9lrr+1378znnzPmez1rftVaqCkmSJEnSdPu6cQcgSZIkSRqexZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmS1EjyY0n+qM/2303ysvWMSRrUaeMOQJIkSRpWkg7w9qp6y1p+TlVtX8v9S8PwzJ00Jkk6Sf7FCPbzn5L8+1HEJEmSpOllcSdJkqSxSfLJJP8myb1JvpzkrUk2NtMfv5Tkvyd5WtP28iT/vyRfSPKRJHPN+l8Evgv49STzSX59mc+sJK9Icn+Sv07yy0m+bkGbNyT52yR/mWR7z/qRHJyV1oLFnbRCYxqEnp3kg0k+n+RIkv+9Wb8b+OfAzzT7+a9rmbskSWvknwEvAL4F+KfA7wI/Czyd7t+rr0iyCfgd4N8D5wGvBt6T5BlV9VrgfwDXV9VZVXX9AJ/5YmAbcAmwA3h5z7bLgCPN5/8S8NYkGTpLaY1Z3Emrs26DUJKnAh8EbgW+EdgJ/H+T/P2q2g+8A/ilZj//dG3SlSRpTf3HqvpcVX2G7vj44ar606r6CvA+4DuAlwAHq+pgVT1WVR8E7gSuXOVn/oeq+nxVfRr4Nbrj6ymfqqo3V9WjwNuA84GNq/wcad1Y3Emrs56D0A8Bn6yq36iqk1V1N/Ae4EdGl44kSWP1uZ73/3OR5bOAbwZ+tJkN84UkXwD+Cd3CazWO9bz/FPDMnuXPnnpTVQ83b89a5edI68a7ZUqrs5JBqPds2unA76/ws74ZuKwZxE45DfjPK9yPJEnT7Bjwn6vq/7XE9lrh/i4EDjfv/w7wwGoDkyaFxZ20dkY1CB0D/qCqXjDkfiRJmmZvBw4l+QHgv9M9YHo5cLSqjtM90PqsFezv3yT5MN0Dsq8E3jjieKV157RMae28HfinSX4gyYYkT0kyl+SCZvugg9D7gW9J8tIkpzevf5TkW1e4H0mSplZVHaN745OfBf6K7sHPf8PX/p59E/AjzR0ubxhgl78N3AXcQ/ca+beOOmZpvaXKg/7SSiT5JPAvquq/N8tvp3vUcG+z/C+Aq6rq+5JcRvcuW/8AeBT4E+BfVtWnk3wn3Yu0n0H3DN8r+nzmFrpHFC+lO4h9BPjXVXVPks3AbwIXAZ2qetHIk5YkqUWSFLC5qo6OOxZplCzuJEmSNFMs7tRWXnMnSZKkVknyXXQfU/QkVeVdL9VanrmTJoCDkCRJkoZlcSdJkiRJLTBV0zKf/vSn10UXXbTq/l/+8pd56lOfOrqApsCs5Wy+7Wa+0+Ouu+7666p6xrjj0OCGHWOXM80/z0tpW05tywfMaRq0LR9Y+5z6jbFTVdxddNFF3Hnnnavu3+l0mJubG11AU2DWcjbfdjPf6ZHkU+OOQSsz7Bi7nGn+eV5K23JqWz5gTtOgbfnA2ufUb4z1OXeSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJE2gJFckOZLkaJI9i2x/dpI/TvKVJK9esO3cJO9O8okk9yX5zvWLXJI0LlP1nLthnTgBe/eOZl+j2o8kSQsl2QDcCLwAOA4cSnKgqj7e0+zzwCuAFy2yizcB/62qfiTJGcA3rHHI6jGtfyNs2TLa2Kf16yBNs5kq7iRJmhKXAker6n6AJLcBO4DHi7uqehB4MMkP9nZMcg7w3cCPNe0eAR5Zn7AnwHpWFEtUQ3Od9QthKZ25veMOQdIYWNxJkjR5NgHHepaPA5cN2PdZwF8Bv5Hk24G7gFdW1ZcXNkyyG9gNsHHjRjqdzjAxM/9nJ5bc9uh5Z/L+/e8cav8DeeaWtf+MxqOnn8n7F/u8XesWwpK2sfKv9aOnn8m2Z47ue9RZh2/3cubPPnvon+uhY+jz/2I1VvN/6axvOX+kMYzS/Pz82L9HozbOnCzuJEmaPFlkXQ3Y9zTgEuCnqurDSd4E7AF+7kk7rNoP7AfYtm1bzc3NrS7aRqfPWbP5XVs469YjQ+1/0rQtp7blA3DWT2xh7s47xxrDqP/GX833aa6zc7RBjFCn02HY3z0DW6cz+50tfX7u1jiGgYq7JFfQnb+/AXhLVb1+wfY0268EHgZ+rKruXq5vkp8CrgdOAr9TVT8zdEaSJE2/48CFPcsXAA+soO/xqvpws/xuusWdNHPmvzT64moqTcAFkHtZPIZRX+vZz3pNmZ5/5tI/d3Nr/NnLFncDXtS9HdjcvC4DbgIu69c3yffQvX7gOVX1lSTfOMrEJEmaYoeAzUkuBj4DXMWAk/2q6rNJjiXZUlVHgOfTc62eJI3DXGfvouvnn7llyW1auUHO3C17UXezfEtVFXBHcwvm84GL+vT9l8Drq+or8PiF4ZIkzbyqOpnkeuB2ujNfbq6qw0mubbbvS/JNwJ3AOcBjSV4FbK2qLwI/BbyjuVPm/cA148hD0mTw7OXsGKS4G+Si7sXabFqm77cA35XkF4H/Bby6qg4t/PBRXux95pnzbNmy+v69puU/SRsvUu3HfNvNfDVLquogcHDBun097z9Ld7rmYn3vAbatZXySpMkzSHE3yEXdS7Xp1/c04GnA5cA/At6V5FnN2b+vNR7hxd7vfGeHI0dW37/Xzsm9LvUJ1vUi1Qlgvu1mvpIkSUsbpLgb5KLupdqc0afvceC9TTH3J0keA55O9/bNkiRJkqQV+LoB2jx+UXczd/8q4MCCNgeAq9N1OfBQVZ1Ypu9vAd8LkORb6BaCfz1sQpIkSZI0i5Y9czfIRd10rwm4EjhK91EI1/Tr2+z6ZuDmJB8DHgFetnBKpiRJkiRpMAM9526Ai7oLuG7Qvs36R4CXrCRYSZIkSdLiBpmWKUmSJEmacBZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRNoCRXJDmS5GiSPYtsf3aSP07ylSSvXmT7hiR/muT96xOxJGncLO4kSZowSTYANwLbga3AziRbFzT7PPAK4A1L7OaVwH1rFqQkaeJY3EmSNHkuBY5W1f1V9QhwG7Cjt0FVPVhVh4CvLuyc5ALgB4G3rEewkqTJcNq4A5AkSU+yCTjWs3wcuGwF/X8N+Bng7H6NkuwGdgNs3LiRTqezoiAXmt+1Zcltj553Zt/t06htObUtHzCnadC2fKB/TsP+nl2OxZ0kSZMni6yrgTomPwQ8WFV3JZnr17aq9gP7AbZt21Zzc32bL6uzd++S2+Z3beGsW48Mtf9J07ac2pYPmNM0aFs+0D+nuc7ONf1sp2VKkjR5jgMX9ixfADwwYN/nAS9M8km60zm/N8nbRxueJGkSWdxJkjR5DgGbk1yc5AzgKuDAIB2r6jVVdUFVXdT0+72qesnahSpJmhROy5QkacJU1ckk1wO3AxuAm6vqcJJrm+37knwTcCdwDvBYklcBW6vqi+OKW5I0XhZ3kiRNoKo6CBxcsG5fz/vP0p2u2W8fHaCzBuFJkiaQ0zIlSZIkqQUs7iRJkiSpBSzuJEmSJKkFBiruklyR5EiSo0n2LLI9SW5ott+b5JLl+ibZm+QzSe5pXleOJiVJkiRJmj3LFndJNgA3AtuBrcDOJFsXNNsObG5eu4GbBuz7q1X13OZ1EEmSJEnSqgxy5u5S4GhV3V9Vj9B9IOqOBW12ALdU1x3AuUnOH7CvJEmSJGlIgzwKYRNwrGf5OHDZAG02DdD3+iRX031Oz09X1d8u/PAku+meDWTjxo10Op0BQl7cmWfOs2XL6vv3GiKMdTU/Pz/U12zamG+7ma8kSdLSBinussi6GrBNv743Aa9rll8H/Arw8ic1rtoP7AfYtm1bzc3NDRDy4t75zg5Hjqy+f6+dO0eymzXX6XQY5ms2bcy33cxXkiRpaYMUd8eBC3uWLwAeGLDNGUv1rarPnVqZ5M3A+weOWpIkSZL0BINcc3cI2Jzk4iRnAFcBBxa0OQBc3dw183Lgoao60a9vc03eKS8GPjZkLpIkSZI0s5Y9c1dVJ5NcD9wObABurqrDSa5ttu8DDgJXAkeBh4Fr+vVtdv1LSZ5Ld1rmJ4GfGGFekiRJkjRTBpmWSfOYgoML1u3reV/AdYP2bda/dEWRSpIkSZKWNNBDzCVJ0vpKckWSI0mOJtmzyPZnJ/njJF9J8uqe9Rcm+f0k9yU5nOSV6xu5JGlcBjpzJ0mS1k+SDcCNwAvo3rTsUJIDVfXxnmafB14BvGhB95N0Hy90d5KzgbuSfHBBX0lSC3nmTpKkyXMpcLSq7q+qR4DbgB29Darqwao6BHx1wfoTVXV38/5LwH10nzsrSWo5iztJkibPJuBYz/JxVlGgJbkI+A7gw6MJS5I0yZyWKUnS5Mki62pFO0jOAt4DvKqqvrhEm93AboCNGzfS6XRWGOYTze/asuS2R887s+/2adS2nNqWD5jTNGhbPtA/p2F/zy7H4k6SpMlzHLiwZ/kC4IFBOyc5nW5h946qeu9S7apqP7AfYNu2bTU3N7eqYE/p7N275Lb5XVs469YjQ+1/0rQtp7blA+Y0DdqWD/TPaa6zc00/22mZkiRNnkPA5iQXJzkDuAo4MEjHJAHeCtxXVW9cwxglSRPGM3eSJE2YqjqZ5HrgdmADcHNVHU5ybbN9X5JvAu4EzgEeS/IqYCvwHOClwEeT3NPs8meb585KklrM4k6SpAnUFGMHF6zb1/P+s3Snay70Ryx+zZ4kqeWclilJkiRJLWBxJ0mSJEktYHEnSZIkSS1gcSdJkiRJLWBxJ0mSJEktYHEnSZIkSS1gcSdJkiRJLWBxJ0mSJEktYHEnSZIkSS0wUHGX5IokR5IcTbJnke1JckOz/d4kl6yg76uTVJKnD5eKJEmSJM2uZYu7JBuAG4HtwFZgZ5KtC5ptBzY3r93ATYP0TXIh8ALg00NnIkmSJEkzbJAzd5cCR6vq/qp6BLgN2LGgzQ7gluq6Azg3yfkD9P1V4GeAGjYRSZIkSZplpw3QZhNwrGf5OHDZAG029eub5IXAZ6rqI0mW/PAku+meDWTjxo10Op0BQl7cmWfOs2XL6vv3GiKMdTU/Pz/U12zamG+7ma8kSdLSBinuFqu8Fp5pW6rNouuTfAPwWuD7l/vwqtoP7AfYtm1bzc3NLddlSe98Z4cjR1bfv9fOnSPZzZrrdDoM8zWbNubbbuarWZLkCuBNwAbgLVX1+gXbnw38BnAJ8NqqesOgfSVJ7TTItMzjwIU9yxcADwzYZqn1fxe4GPhIkk826+9O8k0rCV6SpDYa8Hr3zwOvAN6wir6SpBYapLg7BGxOcnGSM4CrgAML2hwArm7umnk58FBVnViqb1V9tKq+saouqqqL6BaBl1TVZ0eVmCRJU2zZ692r6sGqOgR8daV9JUnttOy0zKo6meR64Ha60zturqrDSa5ttu8DDgJXAkeBh4Fr+vVdk0wkSWqPQa53H7rvKK9rB5jftWXJbY+ed2bf7dOobTm1LR8wp2nQtnygf05rfS39INfcUVUH6RZwvev29bwv4LpB+y7S5qJB4pAkaUYMcr370H1HeV07QGfv3iW3ze/awlm3Hhlq/5OmbTm1LR8wp2nQtnygf05znbW9ccdADzGXJEnrapDr3deiryRpilncSZI0eQa53n0t+kqSpthA0zIlSdL6GeR69+YO03cC5wCPJXkVsLWqvuj17pI0myzuJEmaQANc7/5ZulMuB+orSWo/p2VKkiRJUgvM1Jm7s790grnO3tHsbJjd9LmbmCRJkiSthmfuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZpASa5IciTJ0SR7FtmeJDc02+9NcknPtn+V5HCSjyV5Z5KnrG/0kqRxsLiTJGnCJNkA3AhsB7YCO5NsXdBsO7C5ee0Gbmr6bgJeAWyrqm8DNgBXrVPokqQxsriTJGnyXAocrar7q+oR4DZgx4I2O4BbqusO4Nwk5zfbTgO+PslpwDcAD6xX4JKk8TltkEZJrgDeRPfo31uq6vULtqfZfiXwMPBjVXV3v75JXkd3YHoMeLDp4+AjSRJsAo71LB8HLhugzaaqujPJG4BPA/8T+EBVfWCxD0mym+5ZPzZu3Ein0xkq6PldW5bc9uh5Z/bdPo3allPb8gFzmgZtywf65zTs79nlLFvc9UwNeQHdgeNQkgNV9fGeZr1TQy6jOzXksmX6/nJV/VzzGa8A/h1w7cgykyRpemWRdTVImyRPo3vw9GLgC8BvJnlJVb39SY2r9gP7AbZt21Zzc3PDxExn794lt83v2sJZtx4Zav+Tpm05tS0fMKdp0LZ8oH9Oc52da/rZg0zLHGZqyJJ9q+qLPf2fypMHLUmSZtVx4MKe5Qt48tTKpdp8H/CXVfVXVfVV4L3AP17DWCVJE2KQaZmrnhqyXN8kvwhcDTwEfM9iHz7KKSOjPO3bOXuYzp2RxDCI+fn5NT/9O0nMt93MVzPkELA5ycXAZ+jeEGXXgjYHgOuT3EZ3bH2oqk4k+TRweZJvoDst8/nAnesXuiRpXAYp7lY9NWS5vlX1WuC1SV4DXA/8/JMaj3DKyPv3v3Nkp32Hmrmyc21Px/bqdDoMO81mmphvu5mvZkVVnUxyPXA73WvWb66qw0mubbbvAw7Svdb9KN3r3a9ptn04ybuBu4GTwJ/SjKOSpHYbpLgbZmrIGQP0BbgV+B0WKe4kSZpFVXWQbgHXu25fz/sCrlui78/jmCpJM2eQa+4enxqS5Ay6U0MOLGhzALi6eaDq5TRTQ/r1TbK5p/8LgU8MmYskSZIkzaxlz9wNOTVk0b7Nrl+fZAvdRyF8Cu+UKUmSJEmrNtBz7oacGvKkvs36f7aiSCVJkiRJSxpkWqYkSZIkacJZ3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILWNxJkjSBklyR5EiSo0n2LLI9SW5ott+b5JKebecmeXeSTyS5L8l3rm/0kqRxsLiTJGnCJNkA3AhsB7YCO5NsXdBsO7C5ee0GburZ9ibgv1XVs4FvB+5b86AlSWNncSdJ0uS5FDhaVfdX1SPAbcCOBW12ALdU1x3AuUnOT3IO8N3AWwGq6pGq+sI6xi5JGhOLO0mSJs8m4FjP8vFm3SBtngX8FfAbSf40yVuSPHUtg5UkTYbTxh2AJEl6kiyyrgZscxpwCfBTVfXhJG8C9gA/96QPSXbTndLJxo0b6XQ6w8TM/K4tS2579Lwz+26fRm3LqW35gDlNg7blA/1zGvb37HIs7iRJmjzHgQt7li8AHhiwTQHHq+rDzfp30y3unqSq9gP7AbZt21Zzc3NDBd3Zu3fJbfO7tnDWrUeG2v+kaVtObcsHzGkatC0f6J/TXGfnmn620zIlSZo8h4DNSS5OcgZwFXBgQZsDwNXNXTMvBx6qqhNV9VngWJJTh42fD3x83SKXJI2NZ+4kSZowVXUyyfXA7cAG4OaqOpzk2mb7PuAgcCVwFHgYuKZnFz8FvKMpDO9fsE2S1FIWd5IkTaCqOki3gOtdt6/nfQHXLdH3HmDbWsYnSZo8A03LHPJBqov2TfLLzcNV703yviTnjiQjSZIkSZpByxZ3wzxIdZm+HwS+raqeA/wZ8Jqhs5EkSZKkGTXImbtVP0i1X9+q+kBVnWz630H3Ll+SJEmSpFUYpLgb5kGqg/QFeDnwuwPEIkmSJElaxCA3VBnmQarL9k3yWuAk8I5FP3yED1gd5UMSO2cP07kzkhgGMT8/v+YPS5wk5ttu5itJkrS0QYq7YR6keka/vkleBvwQ8Pzmrl9PMsoHrL5//ztH9pDEoZ7zunNtH17Yq9PpMOxDaaeJ+bab+UqSJC1tkGmZq36Qar++Sa4A/i3wwqp6eET5SJIkSdJMWvbM3TAPUl2qb7PrXwfOBD6YBOCOqrp2lMlJkiRJ0qwY6CHmQz5I9Ul9m/V/b0WRSpIkSZKWNNBDzCVJkiRJk83iTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpKkCZTkiiRHkhxNsmeR7UlyQ7P93iSXLNi+IcmfJnn/+kUtSRoniztJkiZMkg3AjcB2YCuwM8nWBc22A5ub127gpgXbXwnct8ahSpImiMWdJEmT51LgaFXdX1WPALcBOxa02QHcUl13AOcmOR8gyQXADwJvWc+gJUnjNdBz7iRJ0rraBBzrWT4OXDZAm03ACeDXgJ8Bzu73IUl20z3rx8aNG+l0OsPEzPyuLUtue/S8M/tun0Zty6lt+YA5TYO25QP9cxr29+xyLO4kSZo8WWRdDdImyQ8BD1bVXUnm+n1IVe0H9gNs27at5ub6Nl9WZ+/eJbfN79rCWbceGWr/k6ZtObUtHzCnadC2fKB/TnOdnWv62U7LlCRp8hwHLuxZvgB4YMA2zwNemOSTdKdzfm+St69dqJKkSWFxJ0nS5DkEbE5ycZIzgKuAAwvaHACubu6aeTnwUFWdqKrXVNUFVXVR0+/3quol6xq9JGksnJYpSdKEqaqTSa4Hbgc2ADdX1eEk1zbb9wEHgSuBo8DDwDXjileSNBks7iRJmkBVdZBuAde7bl/P+wKuW2YfHaCzBuFJkiaQ0zIlSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFBiruklyR5EiSo0n2LLI9SW5ott+b5JLl+ib50SSHkzyWZNto0pEkSZKk2bRscZdkA3AjsB3YCuxMsnVBs+3A5ua1G7hpgL4fA34Y+MPh05AkSZKk2TbImbtLgaNVdX9VPQLcBuxY0GYHcEt13QGcm+T8fn2r6r6qOjKyTCRJkiRphp02QJtNwLGe5ePAZQO02TRg376S7KZ7NpCNGzfS6XRW0v0JHj3vTOZ3bVl1/16ds4fp3BlJDIOYn58f6ms2bcy33cxXkiRpaYMUd1lkXQ3YZpC+fVXVfmA/wLZt22pubm4l3Z/g/fvfyVm3juZk4RBhwM6dI4lhEJ1Oh2G+ZtPGfNvNfCVJkpY2SHF3HLiwZ/kC4IEB25wxQF9JkiRJ0pAGuebuELA5ycVJzgCuAg4saHMAuLq5a+blwENVdWLAvpIkaYHV3qk6yYVJfj/Jfc1dqV+5/tFLksZh2eKuqk4C1wO3A/cB76qqw0muTXJt0+wgcD9wFHgz8JP9+gIkeXGS48B3Ar+T5PaRZiZJ0pQa5k7VwEngp6vqW4HLgesW6StJaqFBpmVSVQfpFnC96/b1vC/gukH7NuvfB7xvJcFKkjQjHr/bNECSU3eb/nhPm8fvVA3ckeTcJOc3M2dOAFTVl5LcR/cGZx9HktRqAxV3kiRpXQ1zp+oTp1YkuQj4DuDDi33IKO9IDfS9I/Uo71g9KdqWU9vyAXOaBm3LB/rntNZ3wba4kyRp8gxzp+ruxuQs4D3Aq6rqi4t9yCjvSA3Q2bt3yW3zu7aM7I7Vk6JtObUtHzCnadC2fKB/TnOdtb1r/iA3VJEkSetrmDtVk+R0uoXdO6rqvWsYpyRpgljcSZI0eVZ9p+okAd4K3FdVb1zfsCVJ4+S0TEmSJkxVnUxy6m7TG4CbT92putm+j+7Nyq6ke6fqh4Frmu7PA14KfDTJPc26n21ucCZJajGLO0mSJtBq71RdVX/E4tfjSZJazmmZkiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AKnjTuAadXpDNF37xOX9+5drJUkSZIkDW6gM3dJrkhyJMnRJHsW2Z4kNzTb701yyXJ9k5yX5INJ/rz592mjSUmSpOm3FmOvJKndli3ukmwAbgS2A1uBnUm2Lmi2HdjcvHYDNw3Qdw/woaraDHyoWZYkaeat4dgrSWqxQc7cXQocrar7q+oR4DZgx4I2O4BbqusO4Nwk5y/Tdwfwtub924AXDZeKJEmtsVZjrySpxQa55m4TcKxn+Thw2QBtNi3Td2NVnQCoqhNJvnEFcbfKqK6589o9SWqNtRp7JUktNkhxl0XW1YBtBunb/8OT3XSnmwDMJzmykv4LPB346yH6j8Yf/MKa7PYXFt/tZOS8fsy33cx3enzzuAOYcusy9o54jO3vD6b653lxbcupbfmAOU2DtuUD/XPKSOqAJcfYQYq748CFPcsXAA8M2OaMPn0/l+T85qzd+cCDi314Ve0H9g8Q57KS3FlV20axr2kxazmbb7uZr2bIWo29TzDKMXY5bfx5bltObcsHzGkatC0fGG9Og1xzdwjYnOTiJGcAVwEHFrQ5AFzd3LnrcuChZsplv74HgJc1718G/PaQuUiS1BZrNfZKklps2TN3VXUyyfXA7cAG4OaqOpzk2mb7PuAgcCVwFHgYuKZf32bXrwfeleTHgU8DPzrSzCRJmlJrOPZKklpsoIeYV9VBuoNI77p9Pe8LuG7Qvs36vwGev5JgR2Bdpp5MmFnL2XzbzXw1M9Zi7B2zNv48ty2ntuUD5jQN2pYPjDGndMcGSZIkSdI0G+SaO0mSJEnShJuZ4i7JFUmOJDmaZM+44xmVJJ9M8tEk9yS5s1l3XpIPJvnz5t+n9bR/TfM1OJLkB8YX+WCS3JzkwSQf61m34vyS/MPm63Q0yQ1JFrtV+Ngtke/eJJ9pvsf3JLmyZ9u053thkt9Pcl+Sw0le2axv5fe4T76t/R5LvZK8Lsm9zc/5B5I8c9wxDSPJLyf5RJPT+5KcO+6YhpXkR5vfT48lmdo7GLbx777F/kaYZkuNidMqyVOS/EmSjzT5rM2zz5ZTVa1/0b2g/C+AZ9G9RfRHgK3jjmtEuX0SePqCdb8E7Gne7wH+Q/N+a5P7mcDFzddkw7hzWCa/7wYuAT42TH7AnwDfSff5T78LbB93bivIdy/w6kXatiHf84FLmvdnA3/W5NXK73GffFv7Pfblq/cFnNPz/hXAvnHHNGQ+3w+c1rz/D6d+V03zC/hWYAvQAbaNO55V5tDKv/sW+xthml9LjYnjjmuIfAKc1bw/HfgwcPl6xzErZ+4uBY5W1f1V9QhwG7BjzDGtpR3A25r3bwNe1LP+tqr6SlX9Jd07rF26/uENrqr+EPj8gtUryi/d5yieU1V/XN3/cbf09JkoS+S7lDbke6Kq7m7efwm4D9hES7/HffJdylTnKy1UVV/sWXwqSzxcfVpU1Qeq6mSzeAfdZwpOtaq6r6rW7mH266OVf/et8G+EibeKMXGiVdd8s3h681r333GzUtxtAo71LB9nin94FijgA0nuSrK7Wbexus86ovn3G5v1bfk6rDS/Tc37heunyfXNtJ+be6YotirfJBcB30H3SFfrv8cL8oUZ+B5LAEl+Mckx4J8D/27c8YzQy+meRdf4teXvnZmxyJg4lZJsSHIP8CDwwapa93xmpbhb7FqUqT5a2ON5VXUJsB24Lsl392nb5q8DLJ3ftOd9E/B3gecCJ4Bfada3Jt8kZwHvAV614Mj+k5ousm7qcl4k39Z/jzU7kvz3JB9b5LUDoKpeW1UXAu8Arh9vtMtbLp+mzWuBk3RzmniD5DTl/B05RVbwN8DEq6pHq+q5dM/iX5rk29Y7hoGec9cCx4ELe5YvAB4YUywjVVUPNP8+mOR9dKcifC7J+VV1opm+9WDTvC1fh5Xmd5wnTpWZqryr6nOn3id5M/D+ZrEV+SY5ne4v9XdU1Xub1a39Hi+Wb9u/x5otVfV9Aza9Ffgd4OfXMJyhLZdPkpcBPwQ8v5kmPfFW8D2aVm35e6f1lvgbYOpV1ReSdIArgHW9Ac6snLk7BGxOcnGSM4CrgANjjmloSZ6a5OxT7+le2P0xurm9rGn2MuC3m/cHgKuSnJnkYmAz3ZsyTJsV5ddM6/tSksubOwpe3dNn4jXFzSkv5mu/JKY+3ya+twL3VdUbeza18nu8VL5t/h5LvZJs7ll8IfCJccUyCkmuAP4t8MKqenjc8ehxrfy7r236/A0wlZI849Qdc5N8PfB9jOF33Eycuauqk0muB26newelm6vq8JjDGoWNwPu6/zc4Dbi1qv5bkkPAu5L8OPBp4EcBqupwkncBH6c7feS6qnp0PKEPJsk7gTng6UmO0z3C+3pWnt+/BP4T8PV0r4mYyOsilsh3Lslz6U4p+STwE9COfIHnAS8FPtrMUQf4Wdr7PV4q350t/h5LvV6fZAvwGPAp4NoxxzOsX6d7N9sPNmPxHVU11TkleTHwH4FnAL+T5J6qmvhHJ/Vq6999i/2NUFVvHW9UQ1l0TKyqg+MLaSjnA29LsoHuCbR3VdX7l+kzcpmSGQSSJEmSpD5mZVqmJEmSJLWaxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnbSMJJ9M8n3jjkOSpLZxjJVGy+JOmmArGfSSPDfJXUkebv597hqHJ0nS1FrhGLs/yZEkjyX5sTUOTVo1izupBZKcAfw28HbgacDbgN9u1kuSpOF8BPhJ4O5xByL1Y3EnDSjJ1yXZk+QvkvxNknclOa/ZdlGSSrI7yQNJTiT56QH2uTfJu5P8lyRfSnJ3km9vtv1n4O8A/zXJfJKf6bOrOeA04Neq6itVdQMQ4HuHzVuSpLU24WMsVXVjVX0I+F8jSFdaMxZ30uBeAbwI+N+AZwJ/C9y4oM33AJuB7wf2DDjdYwfwm8B5wK3AbyU5vapeCnwa+KdVdVZV/VKfffx94N6qqp519zbrJUmadJM8xkpTw+JOGtxPAK+tquNV9RVgL/AjSU7rafMLVfXlqvoo8BvAzgH2e1dVvbuqvgq8EXgKcPkKYzsLeGjBuoeAs1e4H0mSxmGSx1hpapy2fBNJjW8G3pfksZ51jwIbe5aP9bz/FPAPBtjv432q6rEkx+ketVyJeeCcBevOAb60wv1IkjQOkzzGSlPDM3fS4I4B26vq3J7XU6rqMz1tLux5/3eABwbY7+N9knwdcEFPv1q0x5MdBp6TJD3rntOslyRp0k3yGCtNDYs7aXD7gF9M8s0ASZ6RZMeCNj+X5BuS/H3gGuC/DLDff5jkh5upJ68CvgLc0Wz7HPCsAfbRoXuE8xVJzkxyfbP+9wboK0nSuE3yGEuSM5I8he7Nyk5P8pSmWJQmij+U0uDeBBwAPpDkS3QHh8sWtPkD4CjwIeANVfWBAfb728D/Qffi8ZcCP9xcGwDw/wH+zyRfSPLqpXZQVY/QvRD9auALwMuBFzXrJUmadBM7xjY+APxP4B8D+5v33z3A50vrKk+8uZ6k1UhyEfCXwOlVdXIF/fYCf6+qXrJGoUmSNNUcY6XBeeZOkiRJklrA4k5aY0l+t3lA6sLXz65wP/98if140xRJ0kxyjJWeyGmZkiRJktQCnrmTJEmSpBaYqoeYP/3pT6+LLrpoRX2+/OUv89SnPnVtAppws5q7ec+eWc19kvO+6667/rqqnjHuODS4WRljjXl9GPP6mLaYpy1emMyY+42xU1XcXXTRRdx5550r6tPpdJibm1ubgCbcrOZu3rNnVnOf5LyTfGrcMWhlZmWMNeb1YczrY9pinrZ4YTJj7jfGOi1TkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklpgqp5zN6y9eydzX5IktcKwg6ODqyQNxTN3kiSNUZIrkhxJcjTJnkW2J8kNzfZ7k1zSrL8wye8nuS/J4SSv7OmzN8lnktzTvK5cz5wkSeMxU2fuJEmaJEk2ADcCLwCOA4eSHKiqj/c02w5sbl6XATc1/54Efrqq7k5yNnBXkg/29P3VqnrDeuUiSRo/z9xJkjQ+lwJHq+r+qnoEuA3YsaDNDuCW6roDODfJ+VV1oqruBqiqLwH3AZvWM3hJ0mTxzJ0kSeOzCTjWs3yc7lm55dpsAk6cWpHkIuA7gA/3tLs+ydXAnXTP8P3twg9PshvYDbBx40Y6nc6Kgp+fn39iny1bVtT/SVb4+avxpJingDGvD2Nee9MWL0xfzBZ3kiSNTxZZVytpk+Qs4D3Aq6rqi83qm4DXNe1eB/wK8PIn7aRqP7AfYNu2bTU3N7ei4DudDk/oM+wNUXbuHK7/AJ4U8xQw5vVhzGtv2uKF6YvZaZmSJI3PceDCnuULgAcGbZPkdLqF3Tuq6r2nGlTV56rq0ap6DHgz3emfkqSWs7iTJGl8DgGbk1yc5AzgKuDAgjYHgKubu2ZeDjxUVSeSBHgrcF9VvbG3Q5LzexZfDHxs7VKQJE0Kp2VKkjQmVXUyyfXA7cAG4OaqOpzk2mb7PuAgcCVwFHgYuKbp/jzgpcBHk9zTrPvZqjoI/FKS59KdlvlJ4CfWJSFJ0lhZ3EmSNEZNMXZwwbp9Pe8LuG6Rfn/E4tfjUVUvHXGYkqQpMNC0zNU+YHW5vkl+qtl2OMkvDZ+OJEmSJM2mZc/cDfOA1X59k3wP3Wf3PKeqvpLkG0eZmCRJkiTNkkHO3K36AavL9P2XwOur6isAVfXgCPKRJEmSpJk0yDV3wzxgtV/fbwG+K8kvAv8LeHVVHVr44aN8wOqwz1btNQ3PMpy2hy6OinnPnlnNfVbzliRJixukuBvmAav9+p4GPA24HPhHwLuSPKu5cPxrjUf4gNVhn63aax2eszq0aXvo4qiY9+yZ1dxnNW9JkrS4QYq7YR6wekafvseB9zbF3J8keQx4OvBXA0cvSZIkSQIGu+Zu1Q9YXabvbwHfC5DkW+gWgn89bEKSJEmSNIuWPXM3zANWl+rb7Ppm4OYkHwMeAV62cEqmJEmSJGkwAz3EfLUPWF2qb7P+EeAlKwlWkiRJkrS4gR5iLkmSJEmabBZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkjVGSK5IcSXI0yZ5FtifJDc32e5Nc0qy/MMnvJ7kvyeEkr+zpc16SDyb58+bfp61nTpKk8bC4kyRpTJJsAG4EtgNbgZ1Jti5oth3Y3Lx2Azc1608CP11V3wpcDlzX03cP8KGq2gx8qFmWJLXcQMXdao8q9uubZG+SzyS5p3ldOZqUJEmaGpcCR6vq/qp6BLgN2LGgzQ7gluq6Azg3yflVdaKq7gaoqi8B9wGbevq8rXn/NuBFa5yHJGkCnLZcg56jii8AjgOHkhyoqo/3NOs9qngZ3aOKlw3Q91er6g0jy0aSpOmyCTjWs3yc7ji6XJtNwIlTK5JcBHwH8OFm1caqOgFQVSeSfONiH55kN92zgWzcuJFOp7Oi4Ofn55/YZ8uWFfV/khV+/mo8KeYpYMzrw5jX3rTFC9MX87LFHT1HFQGSnDqq2FvcPX5UEbgjyblJzgcuGqCvJEmzKousq5W0SXIW8B7gVVX1xZV8eFXtB/YDbNu2rebm5lbSnU6nwxP67N27ov5PsnPncP0H8KSYp4Axrw9jXnvTFi9MX8yDTMtc6ojhIG2W63t9M43zZi/2liTNoOPAhT3LFwAPDNomyel0C7t3VNV7e9p8rjnISvPvgyOOW5I0gQY5czfMUcV+fW8CXtcsvw74FeDlT/rwEU4ZGXa2SK9pODs7baeRR8W8Z8+s5j6rebfMIWBzkouBzwBXAbsWtDlA92DobXSnbD7UTLUM8Fbgvqp64yJ9Xga8vvn3t9cwB0nShBikuBvmqOIZS/Wtqs+dWpnkzcD7F/vwUU4ZGXa2SK91mDkytGk7jTwq5j17ZjX3Wc27TarqZJLrgduBDcDNVXU4ybXN9n3AQeBK4CjwMHBN0/15wEuBjya5p1n3s1V1kG5R964kPw58GvjRdUpJkjRGgxR3wxxV/Kul+p6601fT/8XAx4bORpKkKdMUYwcXrNvX876A6xbp90csPkOGqvob4PmjjVSSNOmWLe6GOaq4VN9m17+U5Ll0p2V+EviJEeYlSZIkSTNlkDN3qz6quFTfZv1LVxSpJEmSJGlJAz3EXJIkSZI02SzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUGKu6SXJHkSJKjSfYssj1Jbmi235vkkhX0fXWSSvL04VKRJEmSpNm1bHGXZANwI7Ad2ArsTLJ1QbPtwObmtRu4aZC+SS4EXgB8euhMJEmSJGmGDXLm7lLgaFXdX1WPALcBOxa02QHcUl13AOcmOX+Avr8K/AxQwyYiSZIkSbPstAHabAKO9SwfBy4boM2mfn2TvBD4TFV9JMmSH55kN92zgWzcuJFOpzNAyF8zPz//eJ8tW1bUta8VhjEWvbnPEvOePbOa+6zmLUmSFjdIcbdY5bXwTNtSbRZdn+QbgNcC37/ch1fVfmA/wLZt22pubm65Lk/Q6XQ41Wfv3hV17WvnztHta6305j5LzHv2zGrus5q3JEla3CDTMo8DF/YsXwA8MGCbpdb/XeBi4CNJPtmsvzvJN60keEmSpt2QNy27OcmDST62oM/eJJ9Jck/zunI9cpEkjdcgxd0hYHOSi5OcAVwFHFjQ5gBwdTMAXQ48VFUnlupbVR+tqm+sqouq6iK6ReAlVfXZUSUmSdKkG+amZY3/BFyxxO5/taqe27wOjjRwSdJEWnZaZlWdTHI9cDuwAbi5qg4nubbZvg84CFwJHAUeBq7p13dNMpEkafo8fuMxgCSnbjz28Z42j9+0DLgjyblJzq+qE1X1h0kuWveoJUkTaZBr7miO+B1csG5fz/sCrhu07yJtLhokDkmSWmaYm5adWGbf1ye5GrgT+Omq+tuFDUZ50zJg+DuXrcMNgqbxRkTGvD6Mee1NW7wwfTEPVNy1xVxn7+h2NsyuRnlnF0nSNBvmpmX93AS8rmn3OuBXgJc/aScjvGkZMPz4tg53K5vGGxEZ8/ow5rU3bfHC9MU8yDV3kiRpbQxz07IlVdXnqurRqnoMeDPd6Z+SpJazuJMkaXyGuWnZkpKc37P4YuBjS7WVJLXHTE3LlCRpkgxz0zKAJO8E5oCnJzkO/HxVvRX4pSTPpTst85PAT6xXTpKk8bG4kyRpjIa8admiF6lV1UtHGaMkaTo4LVOSJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawIeYS5Kkkeh0huy/d/gY9o5gH5I0rTxzJ0mSJEktYHEnSZIkSS1gcSdJkiRJLTBQcZfkiiRHkhxNsmeR7UlyQ7P93iSXLNc3yeuatvck+UCSZ44mJUmSJEmaPcsWd0k2ADcC24GtwM4kWxc02w5sbl67gZsG6PvLVfWcqnou8H7g3w2djSRJkiTNqEHO3F0KHK2q+6vqEeA2YMeCNjuAW6rrDuDcJOf361tVX+zp/1SghsxFkiRJkmbWII9C2AQc61k+Dlw2QJtNy/VN8ovA1cBDwPcMHLUkSZIk6QkGKe6yyLqFZ9mWatO3b1W9FnhtktcA1wM//6QPT3bTnerJxo0b6azwITrz8/OP95nftWVFffvpnD1M586owuirN/dZYt6zZ1Zzn9W8JUnS4gYp7o4DF/YsXwA8MGCbMwboC3Ar8DssUtxV1X5gP8C2bdtqbm5ugJC/ptPpcKpPZ4RPNl1hGE+0c+eowuirN/dZYt6zZ1Zzn9W8JUnS4ga55u4QsDnJxUnOAK4CDixocwC4urlr5uXAQ1V1ol/fJJt7+r8Q+MSQuUiSJEnSzFr2zF1VnUxyPXA7sAG4uaoOJ7m22b4POAhcCRwFHgau6de32fXrk2wBHgM+BVw70swkSZIkaYYMMi2TqjpIt4DrXbev530B1w3at1n/z1YUqSRJkiRpSQM9xFySJEmSNNks7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJGqMkVyQ5kuRokj2LbE+SG5rt9ya5pGfbzUkeTPKxBX3OS/LBJH/e/Pu09chFkjReFneSJI1Jkg3AjcB2YCuwM8nWBc22A5ub127gpp5t/wm4YpFd7wE+VFWbgQ81y5KklrO4kyRpfC4FjlbV/VX1CHAbsGNBmx3ALdV1B3BukvMBquoPgc8vst8dwNua928DXrQWwUuSJstAz7mTJElrYhNwrGf5OHDZAG02ASf67HdjVZ0AqKoTSb5xsUZJdtM9G8jGjRvpdDorCn5+fv4JfeZ3bVlR/4W2nL2yz1/MciksjHkaGPP6MOa1N23xwvTFbHEnSdL4ZJF1tYo2q1JV+4H9ANu2bau5ubkV9e90OvT26ezdO1Q8d87tHKo/wM5ldrEw5mlgzOvDmNfetMUL0xez0zIlSRqf48CFPcsXAA+sos1Cnzs1dbP598Eh45QkTQGLO0mSxucQsDnJxUnOAK4CDixocwC4urlr5uXAQ6emXPZxAHhZ8/5lwG+PMmhJ0mSyuJMkaUyq6iRwPXA7cB/wrqo6nOTaJNc2zQ4C9wNHgTcDP3mqf5J3An8MbElyPMmPN5teD7wgyZ8DL2iWJUkt5zV3kiSNUVUdpFvA9a7b1/O+gOuW6LvoFWZV9TfA80cYpiRpCnjmTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklpgoOIuyRVJjiQ5mmTPItuT5IZm+71JLlmub5JfTvKJpv37kpw7kowkSZIkaQYtW9wl2QDcCGwHtgI7k2xd0Gw7sLl57QZuGqDvB4Fvq6rnAH8GvGbobCRJkiRpRg1y5u5S4GhV3V9VjwC3ATsWtNkB3FJddwDnJjm/X9+q+kDz8FaAO4ALRpCPJEmSJM2kQR5ivgk41rN8HLhsgDabBuwL8HLgvyz24Ul20z0byMaNG+l0OgOE/DXz8/OP95nftWVFffvpnD1M586owuirN/dZYt6zZ1Zzn9W8JUnS4gYp7rLIuhqwzbJ9k7wWOAm8Y7EPr6r9wH6Abdu21dzc3DLhPlGn0+FUn87evSvq288Kw3iinTtHFUZfvbnPEvOePbOa+6zmLUmSFjdIcXccuLBn+QLggQHbnNGvb5KXAT8EPL+qFhaMkiRJkqQBDXLN3SFgc5KLk5wBXAUcWNDmAHB1c9fMy4GHqupEv75JrgD+LfDCqnp4RPlIkiRJ0kxa9sxdVZ1Mcj1wO7ABuLmqDie5ttm+DzgIXAkcBR4GrunXt9n1rwNnAh9MAnBHVV07yuQkSZIkaVYMMi2TqjpIt4DrXbev530B1w3at1n/91YUqSRJkiRpSQM9xFySJEmSNNks7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJkiSpBSzuJEmSJKkFLO4kSZIkqQUs7iRJGqMkVyQ5kuRokj2LbE+SG5rt9ya5ZLm+SfYm+UySe5rXleuVjyRpfCzuJEkakyQbgBuB7cBWYGeSrQuabQc2N6/dwE0D9v3Vqnpu8zq4tplIkiaBxZ0kSeNzKXC0qu6vqkeA24AdC9rsAG6prjuAc5OcP2BfSdIMsbiTJGl8NgHHepaPN+sGabNc3+ubaZw3J3na6EKWJE2q08YdgCRJMyyLrKsB2/TrexPwumb5dcCvAC9/0ocnu+lO9WTjxo10Op2Bgj5lfn7+CX3md21ZUf+Ftpy9ss9fzHIpLIx5Ghjz+jDmtTdt8cL0xWxxJ0nS+BwHLuxZvgB4YMA2ZyzVt6o+d2plkjcD71/sw6tqP7AfYNu2bTU3N7ei4DudDr19Onv3rqj/QnfO7RyqP8DOZXaxMOZpYMzrw5jX3rTFC9MXs9MyJUkan0PA5iQXJzkDuAo4sKDNAeDq5q6ZlwMPVdWJfn2ba/JOeTHwsbVORJI0fgMVd2t0m+YfTXI4yWNJto0mHUmSpkdVnQSuB24H7gPeVVWHk1yb5Nqm2UHgfuAo8GbgJ/v1bfr8UpKPJrkX+B7gX61XTpKk8Vl2WmbPrZZfQHdqyKEkB6rq4z3Nem/TfBnduf6XLdP3Y8APA//3CPORJGmqNI8pOLhg3b6e9wVcN2jfZv1LRxymJGkKDHLmbk1u01xV91XVkZFlIkmSJEkzbJDibi1v0yxJkiRJGoFB7pa5VrdpHsgob9M87C2ae3XOHqZzZ1Rh9DVtt24dFfOePbOa+6zmLUmSFjdIcbcmt2ke1Chv0zzsLZp7DXVH1OXu0zwi03br1lEx79kzq7nPat6SJGlxg0zLXJPbNEuSJEmSRmfZM3dVdTLJqVstbwBuPnWb5mb7Prp36rqS7m2aHwau6dcXIMmLgf8IPAP4nST3VNUPjDpBSZIkSZoFg0zLXKvbNL8PeN9KgpUkSZIkLW6gh5hLkiRJkiabxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1wGnjDkCSJAlgrrN3BHsZxT4kaTp55k6SJEmSWsAzd6vU6QzRd+8Tl/fuXayVJElaqeXG1C1b1n7cdVyXNC4DnblLckWSI0mOJtmzyPYkuaHZfm+SS5brm+S8JB9M8ufNv08bTUqSJE0Px1hJ0qgsW9wl2QDcCGwHtgI7k2xd0Gw7sLl57QZuGqDvHuBDVbUZ+FCzLEnSzHCMlSSN0iBn7i4FjlbV/VX1CHAbsGNBmx3ALdV1B3BukvOX6bsDeFvz/m3Ai4ZLRZKkqeMYK0kamUGuudsEHOtZPg5cNkCbTcv03VhVJwCq6kSSb1zsw5PspnukEmA+yZEBYu71dOCvV9hnbf3BLzxh8Rd+YYl2w5u83NeHec+eWc19kvP+5nEHMCUcY0ftD5YdVNc85jUY1yfv67w8Y14f0xbztMULkxnzkmPsIMVdFllXA7YZpG9fVbUf2L+SPr2S3FlV21bbf5rNau7mPXtmNfdZzbtlHGPXmTGvD2NeH9MW87TFC9MX8yDTMo8DF/YsXwA8MGCbfn0/10wrofn3wcHDliSpFRxjJUkjM0hxdwjYnOTiJGcAVwEHFrQ5AFzd3NHrcuChZjpIv74HgJc1718G/PaQuUiSNG0cYyVJI7PstMyqOpnkeuB2YANwc1UdTnJts30fcBC4EjgKPAxc069vs+vXA+9K8uPAp4EfHWlmX7Pq6SYtMKu5m/fsmdXcZzXv1nCMHQtjXh/GvD6mLeZpixemLOZUrWh6viRJkiRpAg30EHNJkiRJ0mSzuJMkSZKkFmh1cZfkiiRHkhxNsmfc8Yxakk8m+WiSe5Lc2aw7L8kHk/x58+/Tetq/pvlaHEnyA+OLfGWS3JzkwSQf61m34jyT/MPm63U0yQ1JFruN+ERZIve9ST7TfN/vSXJlz7ZW5J7kwiS/n+S+JIeTvLJZ3+rve5+8W/891/SZ1DF2VGPGOsY7st936xjzU5L8SZKPNDH/wqTH3BPHhiR/muT90xBzpvBvvSTnJnl3kk80P9ffOckxJ9nSM77ek+SLSV41yTH3VVWtfNG9uPwvgGcBZwAfAbaOO64R5/hJ4OkL1v0SsKd5vwf4D837rc3X4Ezg4uZrs2HcOQyY53cDlwAfGyZP4E+A76T7bKjfBbaPO7dV5r4XePUibVuTO3A+cEnz/mzgz5r8Wv1975N367/nvqbrNclj7KjGjHWMd2S/79Yx5gBnNe9PBz4MXD7JMffE/q+BW4H3T/rPRhPHJ5myv/WAtwH/onl/BnDupMfcE/sG4LN0HxI+FTEvfLX5zN2lwNGqur+qHgFuA3aMOab1sIPufyqaf1/Us/62qvpKVf0l3buuXbr+4a1cVf0h8PkFq1eUZ7rPeTqnqv64uv8zb+npM7GWyH0prcm9qk5U1d3N+y8B9wGbaPn3vU/eS2lF3ppKEzvGjmLMWI84TxnV77t1jrmqar5ZPL15FRMcM0CSC4AfBN7Ss3qiY17CxMac5By6B1jeClBVj1TVFyY55gWeD/xFVX2K6Yn5Cdpc3G0CjvUsH6f/H0nTqIAPJLkrye5m3cbqPv+I5t9vbNa37eux0jw3Ne8Xrp9W1ye5t5l+dGqaQCtzT3IR8B10jwzPzPd9Qd4wQ99zTYVpG1OmYmwc8vfdumqmN94DPAh8sKomPmbg14CfAR7rWTfpMU/b33rPAv4K+I1m+utbkjyVyY6511XAO5v30xLzE7S5uFvs+pK2PffheVV1CbAduC7Jd/dpOwtfD1g6zzblfxPwd4HnAieAX2nWty73JGcB7wFeVVVf7Nd0kXVTm/siec/M91xToy0/YxOTxwh+362rqnq0qp4LXEB3xsC39Wk+9piT/BDwYFXdNWiXRdaN42dj2v7WO43utOibquo7gC/TndK4lEmIGYAkZwAvBH5zuaaLrJuY339tLu6OAxf2LF8APDCmWNZEVT3Q/Psg8D66p4Q/10zJovn3waZ5274eK83zePN+4fqpU1WfawbVx4A387WpAK3KPcnpdP/QeUdVvbdZ3frv+2J5z8r3XFNl2saUiR4bR/T7biyaKXcd4AomO+bnAS9M8km604i/N8nbmeyYp/FvvePA8eZMLsC76RZ7kxzzKduBu6vqc83yNMT8JG0u7g4Bm5Nc3FTiVwEHxhzTyCR5apKzT70Hvh/4GN0cX9Y0exnw2837A8BVSc5McjGwme4NF6bVivJsTqd/KcnlSQJc3dNnqpz6RdN4Md3vO7Qo9ybOtwL3VdUbeza1+vu+VN6z8D3X1Jm2MXZix8ZR/b5br3gBkjwjybnN+68Hvg/4xCTHXFWvqaoLquoiuj+vv1dVL5nkmKfxb72q+ixwLMmWZtXzgY8zwTH32MnXpmTCdMT8ZIPeeWUaX8CVdO869RfAa8cdz4hzexbdO/V8BDh8Kj/g/wl8CPjz5t/zevq8tvlaHGGK7pxH9z/aCeCrdI+W/Phq8gS20f2l+BfArwMZd26rzP0/Ax8F7qX7C+b8tuUO/BO6UxzuBe5pXle2/fveJ+/Wf899Td9rUsfYUY0Z6xjvyH7frWPMzwH+tIn5Y8C/a9ZPbMwL4p/ja3fLnNiYmdK/9eheQnBn8/PxW8DTpiDmbwD+Bvh/9Kyb6JiXeqUJUJIkSZI0xdo8LVOSJEmSZobFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnSRJkiS1gMWdJEmSJLWAxZ0kSZIktYDFnbSIJJ9M8n2T9tlJvivJkfWOSZKkUXGMldaOxZ00Rarqf1TVluXaJTkvyfuSfDnJp5LsWo/4JEmaVisYY69PcmeSryT5T+sQmjSw08YdgKQ1cSPwCLAReC7wO0k+UlWHxxqVJEnT7wHg3wM/AHz9mGORnsAzd1IfSb4uyZ4kf5Hkb5K8K8l5zbaLklSS3UkeSHIiyU8PsM+9Sd6d5L8k+VKSu5N8+4Jmz01yb5KHmnZPafrOJTm+zP6fCvwz4Oeqar6q/gg4ALx0VV8ESZLWwDSOsQBV9d6q+i3gb1aRtrSmLO6k/l4BvAj434BnAn9L96xYr+8BNgPfD+wZ8DqCHcBvAucBtwK/leT0nu3/O3AFcDHwHODHVhDztwCPVtWf9az7CPD3V7APSZLW2jSOsdJEs7iT+vsJ4LVVdbyqvgLsBX4kSe+U5l+oqi9X1UeB3wB2DrDfu6rq3VX1VeCNwFOAy3u231BVD1TV54H/Sndq5aDOAh5asO4h4OwV7EOSpLU2jWOsNNG85k7q75uB9yV5rGfdo3SvZTvlWM/7TwH/YID9Pt6nqh5rpoE8s2f7Z3veP7xg23LmgXMWrDsH+NIK9iFJ0lqbxjFWmmieuZP6OwZsr6pze15PqarP9LS5sOf936F7ofVyHu+T5OuACwbsN4g/A05Lsrln3bcD3kxFkjRJpnGMlSaaxZ3U3z7gF5N8M0CSZyTZsaDNzyX5hiR/H7gG+C8D7PcfJvnhZurJq4CvAHeMIuCq+jLwXuD/SvLUJM+je/3Bfx7F/iVJGpGpG2ObOE9rbsKyAdiQ5CkLppJKY2NxJ/X3Jrp3mvxAki/RHRwuW9DmD4CjwIeAN1TVBwbY728D/wfdi8dfCvxwc23AqPwk3dszPwi8E/iXPgZBkjRhpnWM/T+B/wnsAV7SvP8/R7h/adVSVeOOQZpKSS4C/hI4vapOrqDfXuDvVdVL1ig0SZKmmmOstDqeuZMkSZKkFnB+sLQGkvwu8F2LbPp/j2j/fwf4+BKbt1bVp0fxOZIkTRrHWGlpTsuUJEmSpBZwWqYkSZIktYDFnSRJkiS1wFRdc/f0pz+9LrroojXb/5e//GWe+tSnrtn+x6FtObUtHzCnadC2fGDtc7rrrrv+uqqesWYfoJFzjF25tuXUtnzAnKZB2/KB8Y6xU1XcXXTRRdx5551rtv9Op8Pc3Nya7X8c2pZT2/IBc5oGbcsH1j6nJJ9as51rTTjGrlzbcmpbPmBO06Bt+cB4x1inZUqSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSCwxV3CW5IsmRJEeT7Flk+7OT/HGSryR59YJt5yZ5d5JPJLkvyXcOE4skSZIkzbJV3y0zyQbgRuAFwHHgUJIDVfXxnmafB14BvGiRXbwJ+G9V9SNJzgC+YbWxSJIkSdKsG+ZRCJcCR6vqfoAktwE7gMeLu6p6EHgwyQ/2dkxyDvDdwI817R4BHhkiFmlqnTgBe/eON4a9jDiALVvGn9QorWc+bfq6SZImwiQMLZMQwywYprjbBBzrWT4OXDZg32cBfwX8RpJvB+4CXllVXx4inmUt90O1Xn+/+cPdWM0XYsTfpJEXNauwZcu4I5A0iZJcQXeWywbgLVX1+gXbnw38BnAJ8NqqekPPtnOBtwDfBhTw8qr643UKXZI0JsMUd1lkXa3gcy8BfqqqPpzkTcAe4Oee9CHJbmA3wMaNG+l0OquLluX/iD7zzHm2bFn9/gc1RAorNj8/v+jX7MSJ9YthKeevoqqZP/NMOiOshrbQGdm+Vmu9fu766TDaCnM136f5L400hFU56+zF14/6566vdfoFsdTvBk0GL32QJK3GMMXdceDCnuULgAdW0Pd4VX24WX433eLuSapqP7AfYNu2bTXM0947c3v7bp/ftYXzbr1z1fsf1Bo+sP5JOlu2MHfnk3OahDNW53dWEcOuLXDrkdHFwOj2tVrr9XO3rlbxfTprjUIZiRH/3PUz19m5Lp/T6XQY5vep1txUXvrQb2KFs2M0FhNw7cOo/+Zazf+ludX8zTVye8cdwPr9LPT7Jq1xDMMUd4eAzUkuBj4DXAXsGqRjVX02ybEkW6rqCPB8egYsSZpVyx2EGpX5XVvoLDHATMYfATNvXS59GOXsGOg/Q2aWZsdMhFVM0Zk/80w673zn6GI4//zR7WuV1nXmxRJGPUtoNf+X5p85/mtAlvq/sp7/j9br6/Do6Wfy/iU+66w1znXVxV1VnUxyPXA73esBbq6qw0mubbbvS/JNwJ3AOcBjSV4FbK2qLwI/BbyjmS5yP3DNcKlIktQa63Lpwyhnx8ByZ+46HDky3P4HsXN9Tn4DE34GfBVnBzpbtjB3ZISzFEa5r1V6/zPXb+bFUkY9S2h6Z/ws8XXYtQXWKZ/1miU0v2sLZy3xc7fWM3SGOXNHVR0EDi5Yt6/n/WfpTtdcrO89wLZhPn9aretRxWcu8Xlz6xeDJGnF1uXSh1Hrd9Z3/plb1ums8Hp8Rn/rdQZ+1Jb8m2GVJrXuldpsqOJO08tpV5I00bz0YZXWs7DqN71ZksbB4k6SpAnjpQ+SpNWwuJMkaQJ56YOm3UTca2ag891Se3zduAOQJEmSJA3P4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklpgqOIuyRVJjiQ5mmTPItufneSPk3wlyasX2b4hyZ8mef8wcUiSJEnSrFt1cZdkA3AjsB3YCuxMsnVBs88DrwDesMRuXgnct9oYJEmSJEldw5y5uxQ4WlX3V9UjwG3Ajt4GVfVgVR0Cvrqwc5ILgB8E3jJEDJIktZKzYyRJK3XaEH03Acd6lo8Dl62g/68BPwOc3a9Rkt3AboCNGzfS6XRWFGSv+V1b+m5/9Lwzl20zbdqWU9vyAXOaBm3LB/rnNMzvWY1Gz+yYF9AdXw8lOVBVH+9pdmp2zIuW2M2p2THnrGGokqQJMkxxl0XW1UAdkx8CHqyqu5LM9WtbVfuB/QDbtm2rubm+zfvq7N3bd/v8ri2cdeuRVe9/ErUtp7blA+Y0DdqWD/TPaa6zc52j0SIenx0DkOTU7JjHi7uqehB4MMkPLuzcMzvmF4F/vS4RS5LGbpji7jhwYc/yBcADA/Z9HvDCJFcCTwHOSfL2qnrJEPFIktQWUzc7BvrPkJm1M+DTqG35gDlNg7blA+OdHTNMcXcI2JzkYuAzwFXArkE6VtVrgNcANGfuXm1hJ0nS46Zudgz0nyEza2fAp1Hb8gFzmgZtywfGOztm1cVdVZ1Mcj1wO7ABuLmqDie5ttm+L8k3AXfSne//WJJXAVur6ovDhy5JUms5O0aStGLDnLmjqg4CBxes29fz/rN0B6R+++gAnWHikCSpZZwdI0lasaGKO0mSNHrOjpEkrYbFnSRJE8jZMZKklRrmIeaSJEmSpAlhcSdJkiRJLWBxJ0mSJEktYHEnSZIkSS1gcSdJkiRJLWBxJ0mSJEktYHEnSZIkSS1gcSdJkiRJLWBxJ0mSJEktYHEnSZIkSS0wVHGX5IokR5IcTbJnke3PTvLHSb6S5NU96y9M8vtJ7ktyOMkrh4lDkiRJkmbdaavtmGQDcCPwAuA4cCjJgar6eE+zzwOvAF60oPtJ4Ker6u4kZwN3Jfnggr6SJEmSpAENc+buUuBoVd1fVY8AtwE7ehtU1YNVdQj46oL1J6rq7ub9l4D7gE1DxCJJUqs4O0aStFKrPnNHtxg71rN8HLhspTtJchHwHcCHl9i+G9gNsHHjRjqdzko/4nHzu7b03f7oeWcu22batC2ntuUD5jQN2pYP9M9pmN+zGg1nx0iSVmOY4i6LrKsV7SA5C3gP8Kqq+uJibapqP7AfYNu2bTU3N7fCML+ms3dv3+3zu7Zw1q1HVr3/SdS2nNqWD5jTNGhbPtA/p7nOznWORot4fHYMQJJTs2MeL9Cq6kHgwSQ/2Nuxqk4AJ5r3X0pyanaMxZ0ktdwwxd1x4MKe5QuABwbtnOR0uoXdO6rqvUPEIUlS20zd7BjoP0Nm1s6AT6O25QPmNA3alg+Md3bMMMXdIWBzkouBzwBXAbsG6ZgkwFuB+6rqjUPEIElSG03d7BjoP0Nm1s6AT6O25QPmNA3alg+Md3bMqou7qjqZ5HrgdmADcHNVHU5ybbN9X5JvAu4EzgEeS/IqYCvwHOClwEeT3NPs8mer6uCqM5EkqT2cHSNJWrFhztzRFGMHF6zb1/P+s3QHpIX+iMWPSkqSJGfHSJJWYajiTpIkjZ6zYyRJq2FxJ0nSBHJ2jCRppYZ5iLkkSZIkaUJY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSCwxV3CW5IsmRJEeT7Flk+7OT/HGSryR59Ur6SpIkSZIGt+riLskG4EZgO7AV2Jlk64JmnwdeAbxhFX0lSZpZHkCVJK3UMGfuLgWOVtX9VfUIcBuwo7dBVT1YVYeAr660ryRJs8oDqJKk1ThtiL6bgGM9y8eBy0bdN8luYDfAxo0b6XQ6Kw70lPldW/puf/S8M5dtM23allPb8gFzmgZtywf65zTM71mNzOMHQQGSnDoI+vFTDarqQeDBJD+40r6SpHYaprjLIutq1H2raj+wH2Dbtm01Nzc34Ec8WWfv3r7b53dt4axbj6x6/5OobTm1LR8wp2nQtnygf05znZ3rHI0WMXUHUKH/QdRZO0gyjdqWD5jTNGhbPjDeA6jDFHfHgQt7li8AHliHvpIktd3UHUCF/gdRZ+0gyTRqWz5gTtOgbfnAeA+gDnPN3SFgc5KLk5wBXAUcWIe+kiS1nQdQJUkrtuozd1V1Msn1wO3ABuDmqjqc5Npm+74k3wTcCZwDPJbkVcDWqvriYn2HzEWSpLZ4/CAo8Bm6B0F3rUNfSdIUG2ZaJlV1EDi4YN2+nvefpXvEcKC+kiTJA6iSpNUZqriTJElrwwOokqSVGuaaO0mSJEnShLC4kyRJkqQWsLiTJEmSpBawuJMkSZKkFrC4kyRJkqQWsLiTJEmSpBawuJMkSZKkFrC4kyRJkqQWsLiTJEmSpBawuJMkSZKkFrC4kyRJkqQWGKq4S3JFkiNJjibZs8j2JLmh2X5vkkt6tv2rJIeTfCzJO5M8ZZhYJEmSJGmWrbq4S7IBuBHYDmwFdibZuqDZdmBz89oN3NT03QS8AthWVd8GbACuWm0skiS1jQdQJUkrNcyZu0uBo1V1f1U9AtwG7FjQZgdwS3XdAZyb5Pxm22nA1yc5DfgG4IEhYpEkqTU8gCpJWo3Thui7CTjWs3wcuGyANpuq6s4kbwA+DfxP4ANV9YHFPiTJbrqDFhs3bqTT6aw64PldW/puf/S8M5dtM23allPb8gFzmgZtywf65zTM71mNzOMHUAGSnDqA+vGeNo8fQAXuSLLYAdSv4gFUSZoZwxR3WWRdDdImydPoDkoXA18AfjPJS6rq7U9qXLUf2A+wbdu2mpubW3XAnb17+26f37WFs249sur9T6K25dS2fMCcpkHb8oH+Oc11dq5zNFrE1B1Ahf4HUWftIMk0als+YE7ToG35wHgPoA5T3B0HLuxZvoAnHxlcqs33AX9ZVX8FkOS9wD8GnlTcSZI0g6buACr0P4g6awdJplHb8gFzmgZtywfGewB1mGvuDgGbk1yc5Ay68/kPLGhzALi6uej7cuChqjpB92ji5Um+IUmA5wP3DRGLJEltMpIDqFX1VeDUAVRJUsuturirqpPA9cDtdAuzd1XV4STXJrm2aXYQuB84CrwZ+Mmm74eBdwN3Ax9t4ti/2lgkSWoZD6BKklZsmGmZVNVBugVc77p9Pe8LuG6Jvj8P/Pwwny9JUhtV1ckkpw6gbgBuPnUAtdm+j+74eyXdA6gPA9c02z6c5NQB1JPAn+IBVEmaCUMVd5IkaW14AFWStFLDXHMnSZIkSZoQFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AIWd5IkSZLUAhZ3kiRJktQCFneSJEmS1AJDFXdJrkhyJMnRJHsW2Z4kNzTb701ySc+2c5O8O8knktyX5DuHiUWSJEmSZtmqi7skG4Abge3AVmBnkq0Lmm0HNjev3cBNPdveBPy3qno28O3AfauNRZKktvEAqiRppYY5c3cpcLSq7q+qR4DbgB0L2uwAbqmuO4Bzk5yf5Bzgu4G3AlTVI1X1hSFikSSpNTyAKklajWGKu03AsZ7l4826Qdo8C/gr4DeS/GmStyR56hCxSJLUJh5AlSSt2GlD9M0i62rANqcBlwA/VVUfTvImYA/wc0/6kGQ33SOSbNy4kU6ns+qA53dt6bv90fPOXLbNtGlbTm3LB8xpGrQtH+if0zC/ZzUyix0cvWyANpuAk3ztAOq3A3cBr6yqLy/8kFGOsdB/nJ21/0fTqG35gDlNg7blA+MdY4cp7o4DF/YsXwA8MGCbAo5X1Yeb9e+mW9w9SVXtB/YDbNu2rebm5lYdcGfv3r7b53dt4axbj6x6/5OobTm1LR8wp2nQtnygf05znZ3rHI0WsS4HUEc5xkL/cXbW/h9No7blA+Y0DdqWD4x3jB1mWuYhYHOSi5OcAVwFHFjQ5gBwdXPR9+XAQ1V1oqo+CxxLcqqkfT7w8SFikSSpTYY5gHqcJx9AvQRJUuut+sxdVZ1Mcj1wO7ABuLmqDie5ttm+DzgIXAkcBR4GrunZxU8B72gKw/sXbJMkaZY9fgAV+AzdA6i7FrQ5AFyf5Da6UzYfqqoTAEmOJdlSVUfwAKokzYxhpmVSVQfpFnC96/b1vC/guiX63gNsG+bzJUlqIw+gSpJWY6jiTpIkrQ0PoEqSVmqYa+4kSZIkSRPC4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJawOJOkiRJklrA4k6SJEmSWsDiTpIkSZJaYKjiLskVSY4kOZpkzyLbk+SGZvu9SS5ZsH1Dkj9N8v5h4pAkqW0cYyVJK7Xq4i7JBuBGYDuwFdiZZOuCZtuBzc1rN3DTgu2vBO5bbQySJLWRY6wkaTWGOXN3KXC0qu6vqkeA24AdC9rsAG6prjuAc5OcD5DkAuAHgbcMEYMkSW3kGCtJWrHThui7CTjWs3wcuGyANpuAE8CvAT8DnN3vQ5LspntEko0bN9LpdFYd8PyuLX23P3remcu2mTZty6lt+YA5TYO25QP9cxrm96xGZl3GWElSuwxT3GWRdTVImyQ/BDxYVXclmev3IVW1H9gPsG3btpqb69u8r87evX23z+/awlm3Hln1/idR23JqWz5gTtOgbflA/5zmOjvXORotYl3G2FEeQIX+B1Fn7SDJNGpbPmBO06Bt+cB4D6AOU9wdBy7sWb4AeGDANj8CvDDJlcBTgHOSvL2qXjJEPJIktcW6jLGjPIAK/Q+iztpBkmnUtnzAnKZB2/KB8R5AHeaau0PA5iQXJzkDuAo4sKDNAeDq5o5elwMPVdWJqnpNVV1QVRc1/X7Pwk6SpMc5xkqSVmzVZ+6q6mSS64HbgQ3AzVV1OMm1zfZ9wEHgSuAo8DBwzfAhS5LUbo6xkqTVGGZaJlV1kO7g0rtuX8/7Aq5bZh8doDNMHJIktY1jrCRppYZ6iLkkSZIkaTJY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILWNxJkiRJUgtY3EmSJElSC1jcSZIkSVILDFXcJbkiyZEkR5PsWWR7ktzQbL83ySXN+guT/H6S+5IcTvLKYeKQJKltHGMlSSu16uIuyQbgRmA7sBXYmWTrgmbbgc3NazdwU7P+JPDTVfWtwOXAdYv0lSRpJjnGSpJWY5gzd5cCR6vq/qp6BLgN2LGgzQ7gluq6Azg3yflVdaKq7gaoqi8B9wGbhohFkqQ2cYyVJK3YaUP03QQc61k+Dlw2QJtNwIlTK5JcBHwH8OHFPiTJbrpHJNm4cSOdTmfVAc/v2tJ3+6Pnnblsm2nTtpzalg+Y0zRoWz7QP6dhfs9qZNZljJUktcswxV0WWVcraZPkLOA9wKuq6ouLfUhV7Qf2A2zbtq3m5uZWFSxAZ+/evtvnd23hrFuPrHr/k6htObUtHzCnadC2fKB/TnOdnescjRaxLmPsKA+gQv+DqLN2kGQatS0fMKdp0LZ8YLwHUIcp7o4DF/YsXwA8MGibJKfTHXTeUVXvHSIOSZLaZl3G2FEeQIX+B1Fn7SDJNGpbPmBO06Bt+cB4D6AOc83dIWBzkouTnAFcBRxY0OYAcHVzR6/LgYeq6kSSAG8F7quqNw4RgyRJbeQYK0lasVWfuauqk0muB24HNgA3V9XhJNc22/cBB4ErgaPAw8A1TffnAS8FPprknmbdz1bVwdXGI0lSWzjGSpJWY5hpmTQDxcEF6/b1vC/gukX6/RGLXysgSZJwjJUkrdxQDzGXJEmSJE0GiztJkiRJagGLO0mSJElqAYs7SZIkSWoBiztJkiRJagGLO0mSJElqAYs7SZIkSWoBiztJkiRJagGLO0mSJElqAYs7SZIkSWoBiztJkiRJagGLO0mSJElqAYs7SZIkSWqBoYq7JFckOZLkaJI9i2xPkhua7fcmuWTQvpIkzTLHWEnSSq26uEuyAbgR2A5sBXYm2bqg2XZgc/PaDdy0gr6SJM0kx1hJ0moMc+buUuBoVd1fVY8AtwE7FrTZAdxSXXcA5yY5f8C+kiTNKsdYSdKKnTZE303AsZ7l48BlA7TZNGBfAJLspntEEmA+yZEhYu7vD3g68Ndrtv9xaFtObcsHzGkatC0f6J9TfmEUn/DNo9jJDHOMnQZty6lt+YA5TYO25QNjHWOHKe6yyLoasM0gfbsrq/YD+1cW2uokubOqtq3HZ62XtuXUtnzAnKZB2/KBdubUMo6xU6BtObUtHzCnadC2fGC8OQ1T3B0HLuxZvgB4YMA2ZwzQV5KkWeUYK0lasWGuuTsEbE5ycZIzgKuAAwvaHACubu7odTnwUFWdGLCvJEmzyjFWkrRiqz5zV1Unk1wP3A5sAG6uqsNJrm227wMOAlcCR4GHgWv69R0qk9FYl6kp66xtObUtHzCnadC2fKCdObWGY+zUaFtObcsHzGkatC0fGGNOqVp0Gr4kSZIkaYoM9RBzSZIkSdJksLiTJEmSpBawuFsgyeuS3JvkniQfSPLMccc0jCS/nOQTTU7vS3LuuGMaVpIfTXI4yWNJpvbWuUmuSHIkydEke8YdzygkuTnJg0k+Nu5YRiHJhUl+P8l9zc/cK8cd0zCSPCXJnyT5SJPPSB62Iw3KMXbyOcZOLsfYyTYpY6zX3C2Q5Jyq+mLz/hXA1qq6dsxhrVqS7wd+r7nA/j8AVNW/HXNYQ0nyrcBjwP8NvLqq7hxzSCuWZAPwZ8AL6N7O/BCws6o+PtbAhpTku4F54Jaq+rZxxzOsJOcD51fV3UnOBu4CXjSt36ckAZ5aVfNJTgf+CHhlVd0x5tA0IxxjJ59j7ORyjJ1skzLGeuZugVODTuOpLPHg12lRVR+oqpPN4h10n3c01arqvqo6Mu44hnQpcLSq7q+qR4DbgB1jjmloVfWHwOfHHceoVNWJqrq7ef8l4D5g03ijWr3qmm8WT29eU/07TtPFMXbyOcZOLsfYyTYpY6zF3SKS/GKSY8A/B/7duOMZoZcDvzvuIAR0f3kd61k+zhT/QpsFSS4CvgP48JhDGUqSDUnuAR4EPlhVU52Ppo9jrNaBY+yUcYwdnZks7pL89yQfW+S1A6CqXltVFwLvAK4fb7TLWy6fps1rgZN0c5p4g+Q05bLIuqk+gt1mSc4C3gO8asGZh6lTVY9W1XPpnmG4NMnUT+3RZHGMnXyOsZokjrGjteqHmE+zqvq+AZveCvwO8PNrGM7QlssnycuAHwKeX1NykeUKvkfT6jhwYc/yBcADY4pFfTTz5t8DvKOq3jvueEalqr6QpANcAbTi4nxNBsfYyecYq0nhGDt6M3nmrp8km3sWXwh8YlyxjEKSK4B/C7ywqh4edzx63CFgc5KLk5wBXAUcGHNMWqC5OPqtwH1V9cZxxzOsJM84dTe/JF8PfB9T/jtO08UxVuvEMXYKOMauURxTcpBp3SR5D7CF7p2iPgVcW1WfGW9Uq5fkKHAm8DfNqjum+c5kAEleDPxH4BnAF4B7quoHxhrUKiS5Evg1YANwc1X94ngjGl6SdwJzwNOBzwE/X1VvHWtQQ0jyT4D/AXyU7u8EgJ+tqoPji2r1kjwHeBvdn7mvA95VVf/XeKPSLHGMnXyOsZPLMXayTcoYa3EnSZIkSS3gtExJkiRJagGLO0mSJElqAYs7SZIkSWoBiztJkiRJagGLO0mSJElqAYs7SZIkSWoBiztJkiRJaoH/P4YC3QPijwhzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1080x864 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "\n",
    "ax=data[target==0].hist(weights=weights[target==0],figsize=(15,12),color='b',alpha=0.5,density=True)\n",
    "ax=ax.flatten()[:data.shape[1]] # to avoid error if holes in the grid of plots (like if 7 or 8 features)\n",
    "data[target==1].hist(weights=weights[target==1],figsize=(15,12),color='r',alpha=0.5,density=True,ax=ax)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4kowHjX4rlIC"
   },
   "source": [
    "# Transform the features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Split the Data into Training and Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Y9j5hdrmrlID",
    "outputId": "75c0925f-74c9-4e24-d169-8a884b510f4c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Xtrain Shape:  (433017, 6)\n",
      "ytrain Shape:  (433017,)\n",
      "Training Weights:  (433017,) \n",
      "\n",
      "Xtest Shape:  (144340, 6)\n",
      "ytest Shape:  (144340,)\n",
      "Test Weights:  (144340,)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "train_size = 0.75 # fraction of sample used for training\n",
    "\n",
    "X_train, X_test, y_train, y_test, weights_train, weights_test = \\\n",
    "    train_test_split(data, target, weights, train_size=train_size)\n",
    "\n",
    "y_train, y_test, weights_train, weights_test = \\\n",
    "    y_train.reset_index(drop=True),y_test.reset_index(drop=True), \\\n",
    "    weights_train.reset_index(drop=True), weights_test.reset_index(drop=True)\n",
    "\n",
    "print (\"Xtrain Shape: \",X_train.shape)\n",
    "print (\"ytrain Shape: \",y_train.shape)\n",
    "print (\"Training Weights: \",weights_train.shape,\"\\n\")\n",
    "print (\"Xtest Shape: \",X_test.shape)\n",
    "print (\"ytest Shape: \",y_test.shape)\n",
    "print (\"Test Weights: \",weights_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Doing an extra data split. Test _and_ Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_test, X_val, y_test, y_val, weights_test, weights_val, = \\\n",
    "    train_test_split(X_test, y_test, weights_test, train_size=0.5, shuffle=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- __Training Dataset:__ The sample of data used to fit the model.\n",
    "- __Validation Dataset:__ The sample used to provide an unbiased evaluation of a model fit on the training dataset while tuning  hyperparameters.\n",
    "- __Test Dataset:__ The sample of data used to provide an unbiased evaluation of a final model fit on the training dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standardize the Data\n",
    "\n",
    "**Scale to Mean of 0 and Variance of 1.0:**   $\\ \\ \\ \\ (x-\\mu)/\\sigma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "X_train = scaler.fit_transform(X_train)\n",
    "X_val = scaler.transform(X_val)\n",
    "X_test = scaler.transform(X_test) #applies the transformation calculated the line above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adjust the Test and Train Signal/Background Weights\n",
    "Train on equal amount of Signal and Background, Test on 'natural' ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class_weights_train: (119.66325440279999, 4.634016907200001)\n",
      "Train : total weight sig 119.66325440279995\n",
      "Train : total weight bkg 119.66325440279999\n",
      "Test : total weight sig 3.0960555168000004\n",
      "Test : total weight bkg 80.03703343919999\n"
     ]
    }
   ],
   "source": [
    "class_weights_train = (weights_train[y_train == 0].sum(), weights_train[y_train == 1].sum())\n",
    "print (\"class_weights_train:\",class_weights_train)\n",
    "for i in range(len(class_weights_train)):\n",
    "    weights_train[y_train == i] *= max(class_weights_train)/ class_weights_train[i] #equalize number of background and signal event\n",
    "    weights_test[y_test == i] *= 1/(1-train_size) #increase test weight to compensate for sampling\n",
    "    \n",
    "print (\"Train : total weight sig\", weights_train[y_train == 1].sum())\n",
    "print (\"Train : total weight bkg\", weights_train[y_train == 0].sum())\n",
    "print (\"Test : total weight sig\", weights_test[y_test == 1].sum())\n",
    "print (\"Test : total weight bkg\", weights_test[y_test == 0].sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yxybCOi-rlIM"
   },
   "source": [
    "# Trying NN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_TjIniZPWlVD"
   },
   "source": [
    "## NN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(119.66325440279999, 4.634016907200001)\n",
      "0         0.000752\n",
      "1         0.000752\n",
      "2         0.000197\n",
      "3         0.001966\n",
      "4         0.000223\n",
      "            ...   \n",
      "433012    0.000052\n",
      "433013    0.000752\n",
      "433014    0.001966\n",
      "433015    0.000264\n",
      "433016    0.000752\n",
      "Name: mcWeight, Length: 433017, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "#Quickly take a look at weights\n",
    "print(class_weights_train)\n",
    "print(weights_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/10\n",
      "13532/13532 [==============================] - 20s 1ms/step - loss: 2.6647e-04\n",
      "Epoch 2/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3702e-04\n",
      "Epoch 3/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3429e-04\n",
      "Epoch 4/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3384e-04\n",
      "Epoch 5/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3278e-04\n",
      "Epoch 6/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3336e-04\n",
      "Epoch 7/10\n",
      "13532/13532 [==============================] - 20s 1ms/step - loss: 2.3168e-04\n",
      "Epoch 8/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3028e-04\n",
      "Epoch 9/10\n",
      "13532/13532 [==============================] - 20s 1ms/step - loss: 2.2908e-04\n",
      "Epoch 10/10\n",
      "13532/13532 [==============================] - 19s 1ms/step - loss: 2.3111e-04\n",
      "Training time: 193.77659010887146\n"
     ]
    }
   ],
   "source": [
    "model = tf.keras.models.Sequential([\n",
    "    tf.keras.layers.Flatten(input_shape=(X_train.shape[1],)), # input layer\n",
    "    tf.keras.layers.Dense(128, activation='relu'), # 1st hiddden layer\n",
    "    tf.keras.layers.Dense(128, activation='relu'), # 2nd hiddden layer\n",
    "    tf.keras.layers.Dense(1,activation=\"sigmoid\") # output layer\n",
    "])\n",
    "\n",
    "model.compile(loss=\"binary_crossentropy\", optimizer=\"adam\")\n",
    "starting_time = time.time( )\n",
    "the_fit = model.fit(X_train, y_train.values, epochs=10,\n",
    "                     sample_weight=weights_train)\n",
    "#not using validation dataset here; keras does not support val weights\n",
    "#Solution: Define your own loss function that takes train and val loss! \n",
    "\n",
    "training_time = time.time( ) - starting_time\n",
    "print(\"Training time:\",training_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x155501a7beb0>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAD4CAYAAADVTSCGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAz3ElEQVR4nO3deXxV9bnv8c+TeR42JDEkhIRJmVQgCXi0VlEQ6wAOFdRWrPYidTj36D1H0WrxWNsDnus9DqUq9vRUa4+ICkKrFRFrT0sdCMQyKkOYQpAEEjKS+bl/7JW4E0KyAwkrw/N+vfZrZ631++317C3mm/Vbv72WqCrGGGNMdwlwuwBjjDF9mwWNMcaYbmVBY4wxpltZ0BhjjOlWFjTGGGO6VZDbBfQ0AwcO1PT0dLfLMMaYXmXDhg1HVDWhrW0WNK2kp6eTk5PjdhnGGNOriMi+k22zoTNjjDHdyoLGGGNMt7KgMcYY060saIwxxnQrCxpjjDHdyoLGGGNMt7Lpzcb0A2VlZRQWFlJXV+d2KaYXCg4OJjExkZiYmFPqb0HTRY5V1fLK3/Zx2ahExqbEul2OMc3Kyso4fPgwKSkphIeHIyJul2R6EVXl+PHjHDx4EOCUwsaGzrpIYIDw7NodrNl22O1SjGmhsLCQlJQUIiIiLGRMp4kIERERpKSkUFhYeEqvYUHTRaLDghmVHMP6vcVul2JMC3V1dYSHh7tdhunlwsPDT3no1YKmC2VneNi4v4Ta+ka3SzGmBTuSMafrdP4NWdB0oex0D9V1jWwpKHW7FGOM6TEsaLpQZroHgPV7bPjMGGOaWNB0oYToUIYOjLTzNMZ0sWXLlvGb3/ymS1/z448/RkTYsmVLp/r95je/QUSoqKjo0npOZu/evYgIf/jDH87I/rqDBU0Xy87wsH5vCY2N6nYpxvQZ3RE0EyZM4JNPPmHYsGGd6nfVVVfxySefEBER0aX19GUWNF0sK91D6fE6dhSWu12KMf1OXV0dDQ0NfrWNiYlh8uTJnZ6Rl5CQwOTJkwkIsF+f/rJPqotlZ9h5GmO60u23387bb7/Nn//8Z0QEEeHxxx8H4JJLLuHGG29kyZIlDBs2jLCwMAoKCvjyyy+ZPXs2gwcPJiIigjFjxvDMM8/Q2PjNjNC2hs5EhGeffZZHHnmEhIQEEhMTueeee6ipqWlu03rorGloa9myZdx1113ExsaSmprKggULWuwP4M0332TEiBGEh4dz6aWXkpubi4h0+mitoaGBxx9/nLS0NEJDQxkzZgz//d//3aLN1q1bmT59Oh6Ph8jISEaNGsXixYubt//1r3/lW9/6FjExMcTExHD++efz5ptvdqoOf/l1ZQARmQ48CwQCv1LVha22i7P9O0AVcLuqbmyvr4h4gDeAdGAvcJOqlojIVGAhEALUAv+iqh85fT4GkoHjzq6nqWqhiIQCrwITgaPALFXd6/SZAzzqtH9SVV/x87M5Janx4STHhvHZnmK+f0F6d+7KmH7hscceY//+/Rw7doxf/vKXAKSmpjZvX7duHbt372bRokVEREQQGxvLjh07OPvss7n11luJjo7miy++YMGCBRw/fpyHH3643f09/fTTTJkyhddee41Nmzbx8MMPM2TIEB588MF2+z344IPccMMNvPXWW6xdu5YnnniCMWPGcNNNNwGQk5PD7NmzufHGG3n++efZvn07s2bNOqXP5Cc/+QlPPfUUCxYsICsri7fffptbb70VEeHmm28G4Nprr+Wcc87htddeIzQ0lK+++oqysjLAe7WIq6++mhkzZvCTn/wEVWXz5s0cO3bslOrpSIdBIyKBwGJgKpAPrBeRVaq6zafZlcAI5zEJeAGY1EHf+cBaVV0oIvOd5YeAI8A1qlogImOB1UCKz75uVdXW91q+EyhR1eEiMhtYBMxywmwBkAkosMHZf4nfn1AniQhZ6R4+23MUVbXvL5ge6V9/v5VtBWWu7Hv0oBgWXDPG7/bDhg3D4/HQ2NjI5MmTT9h+7NgxcnNzOeuss5rXXXbZZVx22WWA9xIqF110EVVVVbz88ssdBk16enrzEcYVV1zBunXrWL58eYdBc/HFF/P0008DMHXqVN5//32WL1/eHDSLFi1i1KhRLF26FBFh+vTp1NXV8dBDD/n9WQAUFxfzzDPP8Oijj/Loo48215mfn8/jjz/OzTffzJEjR8jLy+Odd95h3LhxzZ9Jkx07dlBaWsovfvELoqOjAZg2bVqn6ugMf4bOsoFdqpqnqrXAUmBGqzYzgFfV61MgTkSSO+g7A2g6ungFmAmgqrmqWuCs3wqEOUcs7fF9rbeAy5yjrCuANapa7ITLGmC6H+/5tGRleDhcVsOB4uMdNzbGnJaJEye2CBmA6upqFixYwPDhwwkNDSU4OJgf//jH7Nmzh/r6+nZfr/Uv3NGjR5Ofn99hHR31W79+Pddcc02LPz6vvfbaDl+3tS1btlBVVcV3v/vdFutnzZrFjh07KCwsxOPxMHjwYObNm8cbb7xxwqVjhg0bRlRUFLfccgsrV67stiOZJv4MnaUAB3yW8/EetXTUJqWDvkmqeghAVQ+JSGIb+74ByFXVGp91/yUiDcDbeIfC1Hc/qlovIqXAgHbqakFE5gJzAdLS0tooo3MmOedpPttzlLQBNjPF9DydOaLo6ZKSkk5Y99BDD/GrX/2KBQsWMGHCBOLi4li5ciVPPvkk1dXVREVFnfT14uLiWiyHhIRQXV3dYR0d9fv6669JSEho0ab1sj8OHToEnPi+m5ZLSkpITEzkgw8+4Mc//jF33HEHx48f58ILL+S5555j/PjxxMfH88EHH/Cv//qv3HTTTTQ2NjJt2jSef/55hg4d2umaOuLPEU1bYz+t5+6erI0/fdveqcgYvENgd/msvlVVxwHfch7f74r9q+oSVc1U1cxT+Q/f2vCEKOIigu37NMacAW0NT7/55pvcd999PPjgg1x++eVkZmYSFOTuxerPOussioqKWqxrveyP5ORkgBOOUg4f9l7Q1+Px/qF7zjnn8Pbbb3Ps2DE+/PBDqqurueqqq5onKFxwwQW8//77HDt2jOXLl7Njxw5uueWWTtfjD3+CJh8Y7LOcChT42aa9voed4TWc5+ZPTURSgRXAbaq6u2m9qh50nsuB/8Y7NNdi/yISBMQCxX7W3uUCAoTMId7v0xhjTp+/RxVNjh8/TmjoNyPuDQ0NLF26tDtK81tWVha///3v8Q7CeK1atarTrzN27FgiIiJOmCG2bNkyRo4cecJRUnBwMFOmTOGBBx7g0KFDJwyThYeHc80113DHHXewbds2uoM/Eb8eGCEiGcBBYDbQOvZWAfeKyFK8Q2OlznBYUTt9VwFz8M4wmwOsBBCROOBd4GFVXde0AydA4lT1iIgEA1cDH7Z6rU+AG4GPVFVFZDXwcxGJd9pNA9o/E9hFJmV4+HD7YQrLq0mMDjsTuzSmzzrnnHNYuXIl77zzDqmpqQwaNIhBgwadtP3UqVNZvHgxw4cPx+PxsHjx4hZTlN3w0EMPMWnSJGbPns0PfvADtm/fzssvvwzQqe/keDwe/umf/oknn3ySoKAgMjMzWb58Oe+99x6vv/46AJs2beKf//mfmTVrFkOHDqWkpIRFixZx3nnn4fF4ePfdd/n1r3/NzJkzSUtL4+DBg7z00ktMmTKlW957h0HjnPO4F+/sr0Dg16q6VUTmOdtfBN7DO7V5F97pzT9or6/z0guBZSJyJ7AfaDqzdS8wHHhMRB5z1k0DKoHVTsgE4g2Zl53t/wn8VkR24T2Sme3sv1hEfoo3LAGeUNUzMp6V1fx9mhKuOjf5TOzSmD7r7rvvJjc3lzvuuIOSkhIWLFjQ/F2atjz//PPMmzePe+65h/DwcObMmcN1113H3Llzz1zRrWRmZvL666/zyCOPsHLlSjIzM3nhhReYOnVqp28m9sQTTxAUFMQLL7zA4cOHGT58OK+99hqzZ88GvMN0SUlJ/OxnP6OgoIC4uDguvfRSFi1aBMDw4cMRER555BEKCwtJSEjg6quv5uc//3mXv28A8T2MM5CZmak5Oa1nT3deXUMj5z7+AbOyBvP4tX3nxKvpfbZv386oUaPcLsO04bXXXuP73/8+eXl5ZGRkuF1Oh9r7tyQiG1Q1s61tdivnbhIcGMCEIXF8blcIMMY4fvSjHzF16lTi4+PZuHEjTz75JFdddVWvCJnTYUHTjbLTB/DM2h2UHq8jNjzY7XKMMS47evQod999N0ePHmXAgAHMmjWLp556yu2yup0FTTfKyohHFTbuK+HSc9r6mpAxpj9ZtmyZ2yW4wi6q2Y3GD44nOFD43L5PY4zpxyxoulF4SCDjUmLtPI1xnU36MafrdP4NWdB0s6wMD5vyj1Fd5989MozpasHBwRw/btfdM6fn+PHjBAef2rlmC5pulp3uoa5B+eLAMbdLMf1UYmIiBw8epKqqyo5sTKepKlVVVRw8eJDExFM712yTAbpZ5hAPIvD5nmImDx3gdjmmH2r6MmBBQQF1dXUuV2N6o+DgYJKSkjr9xdImFjTdLDYimLOTou0Cm8ZVTXdRNMYNNnR2BmRneNi4r4T6hsaOGxtjTB9jQXMGZKV7qKxtYNshd+5oaIwxbrKgOQOynQts2jRnY0x/ZEFzBiTFhDFkQIQFjTGmX7KgOUOy0j2s31ts00uNMf2OBc0Zkp3hoaSqjl2FFW6XYowxZ5QFzRmSne6cp7FpzsaYfsavoBGR6SLylYjsEpH5bWwXEXnO2b5JRCZ01FdEPCKyRkR2Os/xzvqpIrJBRDY7zyfcW1REVonIFp/l/xCRL5zHDhE55rOtwWdb52/Q3UWGDIggITqU9XaexhjTz3T4hU0RCQQWA1OBfGC9iKxS1W0+za4ERjiPScALwKQO+s4H1qrqQieA5gMPAUeAa1S1QETG4r0NdIpPPdcDLcafVPV+n+33AeN9Nh9X1fP9+TC6k4iQneGxCQHGmH7HnyOabGCXquapai2wFJjRqs0M4FX1+hSIE5HkDvrOAF5xfn4FmAmgqrmqWuCs3wqEiUgogIhEAQ8AT7ZT783A6368rzMuO91DQWk1+SVVbpdijDFnjD9BkwIc8FnOx+cIo4M27fVNUtVDAM5zW1druwHIVdUaZ/mnwNNAm7+pRWQIkAF85LM6TERyRORTEZl5kn5znTY5RUVFbTXpElnOeRq7HI0xpj/xJ2ikjXWt5+ierI0/fdveqcgYYBFwl7N8PjBcVVe002028Jaq+l6TP01VM4FbgGdEZNgJBakuUdVMVc1MSEjwp7xTcvZZ0USHBdnwmTGmX/EnaPKBwT7LqUCBn23a63vYGV7DeS5saiQiqcAK4DZV3e2svgCYKCJ7gb8CI0Xk41Z1zKbVsFnTMJyq5gEf0/L8zRkVGCBkpdt5GmNM/+JP0KwHRohIhoiE4P1l3nr21irgNmf22WSg1BkOa6/vKmCO8/McYCWAiMQB7wIPq+q6ph2o6guqOkhV04GLgB2qeknTdhE5G4gHPvFZF+9zfmcgcCHgO4nhjMtK97C7qJKjFTUdNzbGmD6gw6BR1XrgXryzv7YDy1R1q4jME5F5TrP3gDxgF/AycHd7fZ0+C4GpIrIT76y0hc76e4HhwGM+05L9udvOzcBSbfnV+1FAjoj8HfgTsLDVbLkzLjsjHoD1e0vcLMMYY84YsUuitJSZmak5OTnd9vq19Y2Me3w1t04awk+uGd1t+zHGmDNJRDY458NPYFcGOMNCggIYnxZnM8+MMf2GBY0LstM9bC0opaKm3u1SjDGm21nQuCA7YwCNChv22XkaY0zfZ0HjgvFpcQQGiF33zBjTL1jQuCAyNIixg2LsSs7GmH7BgsYl2RkevjhwjJr6ho4bG2NML2ZB45KsdA+19Y1syi91uxRjjOlWFjQuabrApl2OxhjT11nQuCQ+MoQRiVH2fRpjTJ9nQeOi7AwPG/aW0NBoV2cwxvRdFjQuys7wUF5Tz/ZDZW6XYowx3caCxkV2IzRjTH9gQeOiQXHhpMaH24QAY0yfZkHjsux0D+v3FmNX0TbG9FUWNC7LyvBwpKKWPUcq3S7FGGO6hQWNy7Iz7Ps0xpi+za+gEZHpIvKViOwSkfltbBcRec7ZvklEJnTUV0Q8IrJGRHY6z/HO+qkiskFENjvPU9rY3yoR2eKzfLuIFPnckfOHPtvmOPvYKSJzWr+W24YOjGRgVIhd98wY02d1GDQiEggsBq4ERgM3i0jrW0NeCYxwHnOBF/zoOx9Yq6ojgLXOMsAR4BpVHQfMAX7bqp7rgYo2Sn1DVc93Hr9y2nqABcAkIBtY0BRoPYWIkDnEYzPPjDF9lj9HNNnALlXNU9VaYCkwo1WbGcCr6vUpECciyR30nQG84vz8CjATQFVzVbXAWb8VCBORUAARiQIeAJ708/1dAaxR1WJVLQHWANP97HvGZGV4OFB8nEOlx90uxRhjupw/QZMCHPBZznfW+dOmvb5JqnoIwHlObGPfNwC5qlrjLP8UeBqoaqutM2z3logM7kTtrptk52mMMX2YP0EjbaxrPRf3ZG386dv2TkXGAIuAu5zl84Hhqrqijea/B9JV9VzgQ745UvJr/yIyV0RyRCSnqKjIn/K61KjkGKJCg2z4zBjTJ/kTNPnAYJ/lVKDAzzbt9T3sDK/hPBc2NRKRVGAFcJuq7nZWXwBMFJG9wF+BkSLyMYCqHvU56nkZmNiJ2lHVJaqaqaqZCQkJbXwE3SswQJgwJJ71e+zWzsaYvsefoFkPjBCRDBEJAWYDq1q1WQXc5sw+mwyUOsNh7fVdhfdkP87zSgARiQPeBR5W1XVNO1DVF1R1kKqmAxcBO1T1EqdPsk8t1wLbnZ9XA9NEJN6ZBDDNWdfjTMrw8NXhckoqa90uxRhjulRQRw1UtV5E7sX7CzoQ+LWqbhWRec72F4H3gO8Au/CeP/lBe32dl14ILBORO4H9wHed9fcCw4HHROQxZ900VW0+4mnDP4rItUA9UAzc7uy/WER+ijfwAJ5Q1R45PtV03bOcfSVMHZ3kcjXGGNN1xC590lJmZqbm5OSc8f1W1zVw7uMfcPuF6TzynVFnfP/GGHM6RGSDqma2tc2uDNBDhAUHcv7gOD6zmWfGmD7GgqYHycqIZ+vBUqpq690uxRhjuowFTQ+Sle6hvlHJ3X/M7VKMMabLWND0IBOHxBMg2PCZMaZPsaDpQaLDghk9KIb1FjTGmD7EgqaHyUr3kHughNr6RrdLMcaYLmFB08Nkp3uormtkS0Gp26UYY0yXsKDpYbLsApvGmD7GgqaHGRgVytCESDtPY4zpMyxoeqDsdO+N0Bob7aoNxpjez4KmB8rO8FBWXc9Xh8vdLsUYY06bBU0P1HSBTbs/jTGmL7Cg6YFS48NJjg2zCQHGmD7BgqYHEhGyMzx8vqcYu7q2Maa3s6DpobLSPRSW17C/uMrtUowx5rRY0PRQ2fZ9GmNMH2FB00MNT4giLiLYgsYY0+v5FTQiMl1EvhKRXSIyv43tIiLPOds3iciEjvqKiEdE1ojITuc53lk/VUQ2iMhm53lKG/tbJSJbfJYfEJFtzr7XisgQn20NIvKF81jl/0fjroAAIcv5Po0xxvRmHQaNiAQCi4ErgdHAzSIyulWzK4ERzmMu8IIffecDa1V1BLDWWQY4AlyjquOAOcBvW9VzPVDRav+5QKaqngu8BTzls+24qp7vPK7t6P32JNnpHvYeraKwvNrtUowx5pT5c0STDexS1TxVrQWWAjNatZkBvKpenwJxIpLcQd8ZwCvOz68AMwFUNVdVC5z1W4EwEQkFEJEo4AHgSd+dq+qfVLXprPmnQKof76vHa7ru2fo9JS5XYowxp86foEkBDvgs5zvr/GnTXt8kVT0E4DwntrHvG4BcVa1xln8KPA20NxXrTuCPPsthIpIjIp+KyMy2OojIXKdNTlFRUTsvfWaNGRRDREggn+856nYpxhhzyoL8aCNtrGv95Y6TtfGnb9s7FRkDLAKmOcvnA8NV9X4RST9Jn+8BmcC3fVanqWqBiAwFPhKRzaq6u0VBqkuAJQCZmZk95osrwYEBTEiL5/O9dkRjjOm9/DmiyQcG+yynAgV+tmmv72FneA3nubCpkYikAiuA23xC4QJgoojsBf4KjBSRj336XA78GLjW5wiIpmE4Vc0DPgbG+/Gee4ysdA9ffl1G6fE6t0sxxphT4k/QrAdGiEiGiIQAs4HWs7dWAbc5s88mA6XOcFh7fVfhPdmP87wSQETigHeBh1V1XdMOVPUFVR2kqunARcAOVb3E6TMeeAlvyPgGVrzP+Z2BwIXANj/ec4+RneFBFTbss9lnxpjeqcOgUdV64F5gNbAdWKaqW0VknojMc5q9B+QBu4CXgbvb6+v0WQhMFZGdwFRnGaf9cOAxn2nJbZ2/8fXvQBTwZqtpzKOAHBH5O/AnYKGq9qqgGZ8WR3Cg8LlNCDDG9FJi19JqKTMzU3Nyctwuo4Xrf7kOEeHtH/2D26UYY0ybRGSDqma2tc2uDNALZGcMYFP+MarrGtwuxRhjOs2CphfIzoinrkHJ3X/M7VKMMabTLGh6gYlDPIjYjdCMMb2TBU0vEBsezNlJ0RY0xpheyYKml5iU4WHDvhLqGxrdLsUYYzrFgqaXyMrwUFXbwNaCMrdLMcaYTrGg6SWy050LbNrwmTGml7Gg6SUSY8JIHxDBZ3YjNGNML2NB04tkpXvI2VtMY6N9ydYY03tY0PQiWRkeSqrq2F3U+r5vxhjTc1nQ9CKTnBuh2fCZMaY3saDpRdI8ESRGh9qEAGNMr2JB04uICFkZHj7fU4xdDNUY01tY0PQy2ekeDpVWk19y3O1SjDHGLxY0vUx2hn2fxhjTu1jQ9DJnJ0UTExZkQWOM6TX8ChoRmS4iX4nILhGZ38Z2EZHnnO2bRGRCR31FxCMia0Rkp/Mc76yfKiIbRGSz8zyljf2tEpEtPsuhIvKGs4/PRCTdZ9scZx87RWRO69fqbQIChMx073kaY4zpDToMGhEJBBYDVwKjgZtFZHSrZlcCI5zHXOAFP/rOB9aq6ghgrbMMcAS4RlXHAXOA37aq53qg9RdJ7gRKVHU48B/AIqetB1gATAKygQVNgdabZWd42F1UyZGKGrdLMcaYDvlzRJMN7FLVPFWtBZYCM1q1mQG8ql6fAnEiktxB3xnAK87PrwAzAVQ1V1ULnPVbgTARCQUQkSjgAeDJNvbf9FpvAZeJiABXAGtUtVhVS4A1wHQ/3nOPluVc9yzHhs+MMb2AP0GTAhzwWc531vnTpr2+Sap6CMB5Tmxj3zcAuara9Kf7T4GngaqT7V9V64FSYICftSMic0UkR0RyioqK2iijZxmXEktYcACf7ylxuxRjjOmQP0Ejbaxr/SWOk7Xxp2/bOxUZg3cI7C5n+XxguKqu6ESNfu1fVZeoaqaqZiYkJPhTnqtCggIYPziez/cedbsUY4zpkD9Bkw8M9llOBQr8bNNe38PO8BrOc2FTIxFJBVYAt6nqbmf1BcBEEdkL/BUYKSIft96/iAQBsUCxn7X3SlkZHrYVlFFeXed2KcYY0y5/gmY9MEJEMkQkBJgNrGrVZhVwmzP7bDJQ6gyHtdd3Fd6T/TjPKwFEJA54F3hYVdc17UBVX1DVQaqaDlwE7FDVS9p4rRuBj9T71fnVwDQRiXcmAUxz1vV62ekeGhU27j/mdinGGNOuDoPGOedxL95f0NuBZaq6VUTmicg8p9l7QB6wC3gZuLu9vk6fhcBUEdkJTHWWcdoPBx4TkS+cR1vnb3z9JzBARHbhnSww39l/Md7zOuudxxPOul5vfFocgQHC53ts+MwY07OJXTOrpczMTM3JyXG7DL/MWLyO0MAAls27wO1SjDH9nIhsUNXMtrbZlQF6sez0eL7IP0Z1XYPbpRhjzElZ0PRiWekeausb2ZRf6nYpxhhzUhY0vVjTFzftumfGmJ7MgqYXi48MYWRSlF33zBjTo1nQ9HJZ6R427CuhodEmdRhjeiYLml4uO8NDRU092w+VuV2KMca0yYKml2u6EZoNnxljeioLml4uOTac1PhwmxBgjOmxLGj6gOwM743Q7Mu3xpieyIKmD8hO93C0spa8I5Vul2KMMSewoOkDspzzNOvtPI0xpgeyoOkDhg6MZGBUiE0IMMb0SBY0fYCIkJXu4XObEGCM6YEsaPqIrHQP+SXHOVR63O1SjDGmBQuaPsK+T2OM6aksaPqIUckxRIUGWdAYY3ocv4JGRKaLyFcisktE5rexXUTkOWf7JhGZ0FFfEfGIyBoR2ek8xzvrp4rIBhHZ7DxP8enzvoj8XUS2isiLIhLorP8Pn7tx7hCRYz59Gny2tb4FdZ8RGCBMHBJvX9w0xvQ4HQaN88t8MXAlMBq4WURGt2p2JTDCecwFXvCj73xgraqOANY6ywBHgGtUdRwwB/itz35uUtXzgLFAAvBdAFW9X1XPV9XzgeeB5T59jjdtU9VrO3q/vVl2hocdhysoqax1uxRjjGnmzxFNNrBLVfNUtRZYCsxo1WYG8Kp6fQrEiUhyB31nAK84P78CzARQ1VxVLXDWbwXCRCTU2dZ05cggIARo66vwNwOv+/G++pym8zR2VGOM6Un8CZoU4IDPcr6zzp827fVNUtVDAM5zYhv7vgHIVdWaphUishooBMqBt3wbi8gQIAP4yGd1mIjkiMinIjKzrTcoInOdNjlFRUVtNekVzk2NJSQowILGGNOj+BM00sa61kcSJ2vjT9+2dyoyBlgE3NWis+oVQDIQCkxp1W028JaqNvisS1PVTOAW4BkRGXZCQapLVDVTVTMTEhL8Ka9HCg0K5PzUOD7fW+J2KcYY08yfoMkHBvsspwIFfrZpr+9hZ3gN57mwqZGIpAIrgNtUdXfrglS1GljFiUN4s2k1bNY0DKeqecDHwPi232bfkJURz5aDpVTW1LtdijHGAP4FzXpghIhkiEgI3l/mrWdvrQJuc2afTQZKneGw9vquwnuyH+d5JYCIxAHvAg+r6rqmHYhIlE8wBQHfAb702X42EA984rMuvun8jogMBC4Etvnxnnutb41IoKFRmbl4HWu3H7YrOhtjXNdh0KhqPXAvsBrYDixT1a0iMk9E5jnN3gPygF3Ay8Dd7fV1+iwEporITmCqs4zTfjjwmM+05EQgElglIpuAv+M9AnrRp9SbgaXa8jfrKCBHRP4O/AlYqKp9OmgmDx3Ai9+bQH2jcucrOcxa8ilfHDjmdlnGmH5M7C/eljIzMzUnJ8ftMk5bXUMjSz/fz7Nrd3KkoparxiXzL1ecTfrASLdLM8b0QSKywTkffuI2C5qW+krQNKmoqWfJ/+Tx8v/kUdfQyPcmD+G+KcMZEBXqdmnGmD7EgqYT+lrQNCksq+aZtTt5Y/0BwoMDmfftodxxUQYRIUFul2aM6QMsaDqhrwZNk12FFTz1/pd8sO0widGhPDB1JDdOTCUo0C57Z4w5de0Fjf126WeGJ0ax5LZM3pp3Aanx4cxfvpkrn/0LH26zGWrGmO5hQdNPZaZ7ePtH/8CL35tAQ6Pyw1dzmPXSp+Tuty97GmO6lgVNPyYiTB+bzOr7L+anM8eSd6SC6375N+753Ub2Hql0uzxjTB9h52ha6evnaNpTUVPPy/+Tx8t/yaO2vpFbJ6Vx32UjGGgz1IwxHbDJAJ3Qn4OmSWF5Nc9+uJOl6w8QFhTAvG8P485v2Qw1Y8zJWdB0ggXNN3YXeWeord7qnaF2/9SRfNdmqBlj2mCzzswpGZYQxUvf985QG+yJ4OHlm5n+7F9YYzPUjDGdYEFjOpSZ7uGteRfw4vcm0tio/C+boWaM6QQLGuMX7wy1s1h9/8U8OXMseUcque6Xf+Pu321gj81QM8a0w87RtGLnaPxTWVPPy3/JY8n/eGeo3TIpjX+0GWrG9Fs2GaATLGg6p7C8mufW7uT1z70z1O769jB+aDPUjOl3LGg6wYLm1OwuquDf3/+K97d+TUJ0KPdfPpKbMm2GmjH9hQVNJ1jQnJ4N+4r5t/e+JGdfCRkDI7lxYiozx6eQEhfudmnGmG502tObRWS6iHwlIrtEZH4b20VEnnO2bxKRCR31FRGPiKwRkZ3Oc7yzfqqIbBCRzc7zFJ8+74vI30Vkq4i8KCKBzvrbRaTI546cP/TpM8fZx04Rabp1tOkmE4d4eHPeBbz0/YkMiAzh31d/xUWLPmL2kk9Ytv4AZdV1bpdojDnDOjyicX6Z78B7u+V8YD1ws+8tkUXkO8B9wHeAScCzqjqpvb4i8hRQrKoLnQCKV9WHRGQ8cFhVC0RkLLBaVVOc/cSoapmICPAW8KaqLhWR24FMVb23Ve0eIAfIBBTYAExU1ZPOy7Ujmq61/2gV73xxkBW5B9lzpJLQoAAuH53E9eNTuHhkAsE2tGZMn9DeEY0/Z2yzgV2qmue82FJgBrDNp80M4FX1ptanIhInIslAejt9ZwCXOP1fAT4GHlLVXJ/X3QqEiUioqtaoaplP3SF4w6M9VwBrVLXY2f8aYDrwuh/v23SBtAER/ONlI7hvynC+OHCMd3IP8vtNh3h30yE8kSFcc24y101I5bzUWLx/Pxhj+hp/giYFOOCznI/3qKWjNikd9E1S1UMAqnpIRBLb2PcNQK6q1jStEJHVeMPvj3iPaprbisjFeI+g7lfVA+3U1YKIzAXmAqSlpbVRhjldIsL4tHjGp8Xz6NWj+fNXRazIPcjr6w/wyif7GDowkuvGpzBzfAqDPRFul2uM6UL+BE1bf2a2PpI4WRt/+ra9U5ExwCJgWovOqleISBjwO2AKsAb4PfC6qtaIyDy8R0hT/N2/qi4BloB36Myf+sypCw70Dp9dPjqJsuo6/rj5EMs3HuTpNTt4es0OstLjuW58KleNSyY2Itjtco0xp8mfAfJ8YLDPcipQ4Geb9voedobXcJ4LmxqJSCqwArhNVXe3LkhVq4FVeIffUNWjPkc9LwMTO1G7cVFMWDCzstJ4464LWDd/Cv9yxdmUVNXxyIrNZP3sQ+b9dgOrt35NbX2j26UaY06RP0c064ERIpIBHARmA7e0arMKuNc5BzMJKHWGw4ra6bsKmAMsdJ5XAohIHPAu8LCqrmvagYhEAdHO6wbhnXjwF2dbctMwHHAtsN35eTXw86YZbXiPjh724z0bF6TEhXPPpcO5+5JhbDlYxvLcfH7/9wLe3/o1cRHBXDUumesnpDAhLd7O5xjTi/j1PRpnVtkzQCDwa1X9mTNEhaq+6MwC+wXeE+1VwA9UNedkfZ31A4BlQBqwH/iuqhaLyKN4w2CnTwnT8A6D/QEIdV7rI7znYupF5N/wBkw9UAz8SFW/dPZzB/CI8zo/U9X/au+92qyznqW+oZG/7DrCio0H+WDb11TXNTJkQAQzz0/huvEppA+MdLtEYwz2hc1OsaDpuSpq6nl/y9esyM3nb7uPogrj0+K4fnwKV507CE9kiNslGtNvWdB0ggVN73Co9DirvihgRe5Bvvy6nKAA4ZKzE7l+QgpTzkkkLDjQ7RKN6VcsaDrBgqb32VZQxorcfFZ+UUBheQ3RYUFcNS6Z68ankJXuISDAzucY090saDrBgqb3amhU/rbbez7n/a1fU1XbQEpcONPHnsW5qbGMS4klfUCkBY8x3cCCphMsaPqGqtp6Pth6mOW5B/k072jz9Oio0CBGD4phXIo3eMamxDJ0oIWPMafLgqYTLGj6nrqGRnYermDLwVI2O4/th8qoccInMiSQ0YNiGOuEz7iUWIYmRBFo4WOM3yxoOsGCpn+ob2hkZ6E3fJoCaNuhMqrrvOETERLI6GRv+DQF0LCESLu/jjEnYUHTCRY0/Vd9QyO7iyqbg2fLwVK2FpRxvK4BgLDgAEYne4fdxjjhMyIxysLHGCxoOsWCxvhqaFTyiirYUlDK5vwyJ3xKqaz1hk9oUACjnPAZm+I9AhqZFG23PzD9jgVNJ1jQmI40Nip5RyrZWlDK5nzv0c/WgjIqauoBCAkKYNRZ0c1Dbk3hExJk4WP6LguaTrCgMaeisVHZe7SSLQXeo57N+aVsKSilvNoJn8AAhiZEMiwximEJUQxLiGRYQhQZAyOJDPXnkoPG9Gyne+MzY0wHAgKEoQlRDE2I4trzBgHe8NlfXOUddjtY2jzz7Y+bD9Ho8/fdoNgwhiVGMXSgbxBFkRQTahcPNX2CBY0x3SQgQEgfGEn6wEiuPndQ8/qa+gb2Ha1id2EFu4sq2F1USV5RBW9vPNg8/AbeaddDnaOfoU74DEuMJH1ApF1ix/QqFjTGnGGhQYGMTIpmZFJ0i/WqSmF5TXP4NAXR+r0lvPPFN7dREoHU+PDmI5+hzjDcsIQoBkaF2FGQ6XEsaIzpIUSEpJgwkmLC+IdhA1tsO17bQN6RCvKKKlsE0ad5R5u/+wMQExbU4uhn6MAohidGkuaJtMkIxjUWNMb0AuEhgYwZFMuYQbEt1jc2KofKqpuPfpqC6K+7inh7Y35zu8AAYYgngqEJkQxPjGZkUhQjk6IZnhhlw3Cm21nQGNOLBQQIKXHhpMSFc/HIhBbbyqvr2HPEOQIqrCTvSAW7Civ4844i6hq8sxFEIM0TwQgnfM4+K5oRidEMTbDzQKbr+BU0IjIdeBbvnS1/paoLW20XZ/t38N5h83ZV3dheXxHxAG8A6cBe4CZVLRGRqXhv7xwC1AL/oqofOX3eB5Kduv8C3KOqDSLyAPBDvHfYLALuUNV9Tp8GYLNT6n5VvbYzH5AxvVV0WDDnpsZxbmpci/V1DY3sO1rJjsMV7Dhczk7n+eOvCql3psMFCKQPiGSEc+QzIskbREMHRtkQnOm0Dr9HIyKBwA5gKpAPrAduVtVtPm2+A9yHN2gmAc+q6qT2+orIU0Cxqi4UkflAvKo+JCLjgcOqWiAiY4HVqpri7CdGVcucYHsLeFNVl4rIpcBnqlolIj8CLlHVWU6fClWN8vcDse/RmP6qtr6RPUcqnfAp9wZRYTn7jlbR4ARQYICQPiCiOXzOdgIofWCkXQ2hnzvd79FkA7tUNc95saXADGCbT5sZwKvqTa1PRSRORJLxHq2crO8M4BKn/yvAx8BDqprr87pbgTARCVXVGlUt86k7BFAAVf2TT59Pge/58b6MMT5CggI4+6xozj6r5Wy4mvoG8ooqWxz9fPl1Oe9v/Zqmv1ODA4WMgZHeIx9nGG5EUjTpAyLsWnDGr6BJAQ74LOfjPWrpqE1KB32TVPUQgKoeEpHENvZ9A5CrqjVNK0RkNd7w+yPeo5rW7nS2NQkTkRy8w2oLVfWd1h1EZC4wFyAtLa2NlzSm/woNCmRUcgyjkmNarK+ua2BXYQU7C71HPzsPl7M5v5T3Nh9qDqCmKyKMTPomfEYmRZPmibDbMPQj/gRNW/8aWo+3nayNP33b3qnIGGARMK1FZ9UrRCQM+B0wBVjj0+d7QCbwbZ8uac4w3FDgIxHZrKq7W73mEmAJeIfO/KnPmP4uLDiw+TYKvo7XegNox+FydhR6j4I27i9h1d+/+S5QSGAAA6NC8ESFEB8RwoDIEOIjvc+eyFA8kcEtnmPDgy2YejF/giYfGOyznAoU+NkmpJ2+h0Uk2TmaSQYKmxqJSCqwAritdSgAqGq1iKzCO/y2xulzOfBj4Nu+R0CqWuA854nIx8B44ITXNMZ0jfCQQMalxjIutWUAVdbUs6uwgq8Ol7O7qIIj5bWUVNVytLKWfUerKK6sbXFlBF8BAnERIcRHBDMgMpT4VkHkG1RNzzZrrufwJ2jWAyNEJAM4CMwGbmnVZhVwr3MOZhJQ6gRIUTt9VwFz8M4wmwOsBBCROOBd4GFVXde0AxGJAqKd1w3CO/HgL8628cBLwHRV9Q2seKBKVWtEZCBwIfCUX5+MMaZLRYYGcd7gOM4bHHfSNjX1DZRU1lFcWet9VNVSXFHzzc/O+j1HKtmw7xglVbXNExVaCw8OxBMZctJHfEQIA6NCSIwOIzEm1IKpG3UYNKpaLyL3AqvxTlH+tapuFZF5zvYXgffw/uLfhXd68w/a6+u89EJgmYjcCewHvuusvxcYDjwmIo8566bhHYZbJSKhzmt9BLzobP93IAp407n8RtM05lHASyLSCATgPUfjO4nBGNODhAYFclZsIGfFhvnVvrFRKa+u52hljffoqOKbo6Tiim/CqaSylt1FFZRU1jbfS6i1uIhgkpzQ8V6hwfucGP3NzwnRoTa77hTYbQJasenNxvRt1XUNzaF0tLKWwrJqCstrOFxW7TxqmtfVt3G01HQU1BxETaEUHdYcUAOiQvvdOSW7TYAxxjjCggNJjg0nOTa83XaNjUpxVW2LAPINosPl1WwpKONIRQ2t/14PEEiIPvGIKCkm1BtMzrr4iBAC+kEgWdAYY0wbAgKEgVGhDIwKPeEac77qGxo5UuETSOVOEDmhlF9Sxcb9JRRX1p7QNzhQSIwOY2BUCNFhwUSHBTmP4JbPoUFtbu8t55UsaIwx5jQEBQZwVmxYh+eVauobKCqv+eaIyAmlw6XVHKmspby6jq/LqimvrqO8up6qk5xL8hUSGEB0WBBRTQEU2jKgYnx+jjrJ+vDgwG6/tYQFjTHGnAGhQYGkxkeQGh/hV/v6hkYqaxooc4KnKYDKa+qoqK6nrLq+5XrneX9xFeXV9ZRV11FRU3/CsF5rgQHiDaLQIM4fHMcvbpnQBe+2JQsaY4zpgYICA4iNCCA2IviUX6OxUamsbQokJ4xq2g6o8up6kv2c7ddZFjTGGNNHBQSIM0R26mHVJXW4undjjDF9ngWNMcaYbmVBY4wxpltZ0BhjjOlWFjTGGGO6lQWNMcaYbmVBY4wxpltZ0BhjjOlWdpuAVpybte07jZcYCBzponJ6O/ssWrLPoyX7PL7RFz6LIaqa0NYGC5ouJiI5J7snQ39jn0VL9nm0ZJ/HN/r6Z2FDZ8YYY7qVBY0xxphuZUHT9Za4XUAPYp9FS/Z5tGSfxzf69Gdh52iMMcZ0KzuiMcYY060saIwxxnQrC5ouIiLTReQrEdklIvPdrsdNIjJYRP4kIttFZKuI/G+3a3KbiASKSK6I/MHtWtwmInEi8paIfOn8G7nA7ZrcJCL3O/+fbBGR10Wke25z6SILmi4gIoHAYuBKYDRws4iMdrcqV9UD/0dVRwGTgXv6+ecB8L+B7W4X0UM8C7yvqucA59GPPxcRSQH+EchU1bFAIDDb3aq6ngVN18gGdqlqnqrWAkuBGS7X5BpVPaSqG52fy/H+Iklxtyr3iEgqcBXwK7drcZuIxAAXA/8JoKq1qnrM1aLcFwSEi0gQEAEUuFxPl7Og6RopwAGf5Xz68S9WXyKSDowHPnO5FDc9AzwINLpcR08wFCgC/ssZSvyViES6XZRbVPUg8H+B/cAhoFRVP3C3qq5nQdM1pI11/X7euIhEAW8D/6SqZW7X4wYRuRooVNUNbtfSQwQBE4AXVHU8UAn023OaIhKPd/QjAxgERIrI99ytqutZ0HSNfGCwz3IqffDwtzNEJBhvyPxOVZe7XY+LLgSuFZG9eIdUp4jIa+6W5Kp8IF9Vm45w38IbPP3V5cAeVS1S1TpgOfAPLtfU5SxousZ6YISIZIhICN6Teatcrsk1IiJ4x+C3q+r/c7seN6nqw6qaqqrpeP9dfKSqfe4vVn+p6tfAARE521l1GbDNxZLcth+YLCIRzv83l9EHJ0cEuV1AX6Cq9SJyL7Aa76yRX6vqVpfLctOFwPeBzSLyhbPuEVV9z72STA9yH/A754+yPOAHLtfjGlX9TETeAjbina2ZSx+8HI1dgsYYY0y3sqEzY4wx3cqCxhhjTLeyoDHGGNOtLGiMMcZ0KwsaY4wx3cqCxhhjTLeyoDHGGNOt/j8/glRHNsk0vwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(the_fit.history['loss'],label=\"training loss\")\n",
    "plt.legend(fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Use the model to make predicions!\n",
    "Evaluate the model based on predictions made with X_test $\\rightarrow$ y_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred_test = model.predict(X_test).ravel()\n",
    "y_pred_train = model.predict(X_train).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "### ROC curves and Area Under the Curve (AUC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "id": "SfB3mkkxW4CP"
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_auc_score # for binary classification if x > 0.5 -> 1 else -> 0\n",
    "from sklearn.utils import class_weight # to set class_weight=\"balanced\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Receiver Operating Characteristic')"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAArU0lEQVR4nO3debwWZf3/8debTVAQUtAQRFxwz/UImqmYpqj1JTM3TFMrslzyZ5qmfc3MzLL6qrlHZmZqhku4YptiqSGaC26EKwgqiIKIiMDn98c1p3N7e5b7wJkz5z73+/l43I+Zuee6Zz5z3+fMZ+aametSRGBmZrWrS9EBmJlZsZwIzMxqnBOBmVmNcyIwM6txTgRmZjXOicDMrMY5EVirSHpK0sii4+goJJ0uaVxB675a0jlFrLutSTpM0j0r+Fn/Ta4kJ4IqJuklSe9JWijptWzH0DvPdUbEFhFxb57rqCdpFUk/lvRKtp3/kXSKJLXH+huJZ6SkmaXvRcS5EfHVnNYnSSdImirpXUkzJf1R0ifyWN+KknSWpGtXZhkR8fuI2KuCdX0k+bXn32Rn5URQ/T4XEb2BbYBtge8WG07rSerWxKw/AnsA+wJ9gMOBscCFOcQgSR3t/+FC4FvACcAawMbArcB+bb2iZn6D3BW5bstEhF9V+gJeAvYsmf4pcEfJ9I7AA8DbwOPAyJJ5awC/AWYBbwG3lsz7LPBY9rkHgK3K1wmsA7wHrFEyb1tgLtA9mz4aeCZb/kRgvZKyARwL/Ad4sZFt2wNYDKxb9v4IYBmwUTZ9L/BjYDIwH/hTWUzNfQf3Aj8C/plty0bAUVnM7wAvAF/Pyq6WlVkOLMxe6wBnAddmZYZm2/Vl4JXsuzijZH29gN9m38czwHeAmU38tsOy7RzezO9/NXAJcEcW77+ADUvmXwjMABYAjwC7lMw7CxgPXJvN/yowHHgw+65mAxcDPUo+swXwZ2Ae8DpwOjAKWAJ8kH0nj2dl+wK/zpbzKnAO0DWbd2T2nf9ftqxzsvf+kc1XNu+N7Dd9AtiSdBDwQba+hcBt5f8HQNcsruez7+QRyv6G/Grkb6noAPxaiR/vw/8Ag4EngQuz6UHAm6Sj6S7AZ7LpAdn8O4A/AB8DugO7Ze9vl/0Djsj+qb6crWeVRtb5N+BrJfGcD1yejX8emA5sBnQDvgc8UFI2sp3KGkCvRrbtPOC+Jrb7ZRp20PdmO5otSTvrm2jYMbf0HdxL2mFvkcXYnXS0vWG2M9oNWARsl5UfSdmOm8YTwa9IO/2tgfeBzUq3KfvOB5N2cE0lgmOAl1v4/a8m7UiHZ/H/HrihZP6XgDWzed8GXgN6lsT9QfY7dcni3Z6UOLtl2/IMcGJWvg9pp/5toGc2PaL8OyhZ963AFdlvshYpUdf/ZkcCS4Hjs3X14sOJYG/SDrxf9jtsBgws2eZzmvk/OIX0f7BJ9tmtgTWL/l/t6K/CA/BrJX689A+wkHTkE8BfgX7ZvFOB35WVn0jasQ8kHdl+rJFlXgb8sOy952hIFKX/dF8F/paNi3T0uWs2fRfwlZJldCHtVNfLpgP4dDPbNq50p1Y27yGyI23Szvy8knmbk44Yuzb3HZR89uwWvuNbgW9l4yOpLBEMLpk/GTgkG38B2Ltk3lfLl1cy7wzgoRZiuxoYVzK9L/BsM+XfArYuiXtSC8s/EbglGz8U+HcT5f77HWTTa5MSYK+S9w4F/p6NHwm8UraMI2lIBJ8GppGSUpdGtrm5RPAcMHpl/7dq7dXR6kSt9T4fEX1IO6lNgf7Z++sBB0p6u/4FfIqUBNYF5kXEW40sbz3g22WfW5dUDVJuPLCTpHWAXUk7wftLlnNhyTLmkZLFoJLPz2hmu+ZmsTZmYDa/seW8TDqy70/z30GjMUjaR9JDkuZl5fel4Tut1Gsl44uA+gv465Str7ntf5Omt7+SdSHp25KekTQ/25a+fHhbyrd9Y0m3ZzceLADOLSm/Lqm6pRLrkX6D2SXf+xWkM4NG110qIv5Gqpa6BHhd0pWSVq9w3a2J0zJOBJ1ERNxHOlr6WfbWDNLRcL+S12oRcV42bw1J/RpZ1AzgR2WfWzUirm9knW8D9wAHAWOA6yM7LMuW8/Wy5fSKiAdKF9HMJv0FGCFp3dI3JQ0n/bP/reTt0jJDSFUec1v4Dj4Sg6RVSFVLPwPWjoh+wJ2kBNZSvJWYTaoSaizucn8FBkuqW5EVSdqFdEZ0EOnMrx+pvr30jqvy7bkMeBYYFhGrk+ra68vPIFWZNaZ8OTNIZwT9S7731SNii2Y+8+EFRlwUEduTqu02JlX5tPi5FuK0JjgRdC4XAJ+RtA3pIuDnJO0tqaukntntj4MjYjap6uZSSR+T1F3SrtkyfgUcI2lEdifNapL2k9SniXVeBxwBHJCN17sc+K6kLQAk9ZV0YKUbEhF/Ie0Mb5K0RbYNO5LqwS+LiP+UFP+SpM0lrQqcDYyPiGXNfQdNrLYHsAowB1gqaR+g9JbG14E1JfWtdDvK3Ej6Tj4maRBwXFMFs+27FLg+i7lHFv8hkk6rYF19SPXwc4Buks4EWjqq7kO6cLxQ0qbAN0rm3Q58XNKJ2W29fSSNyOa9Dgytv+sq+/u6B/i5pNUldZG0oaTdKogbSTtkf3/dgXdJNw0sK1nXBs18fBzwQ0nDsr/frSStWcl6a5kTQScSEXOAa4D/jYgZwGjSUd0c0pHSKTT85oeTjpyfJV0cPjFbxhTga6RT87dIF3yPbGa1E0h3uLweEY+XxHIL8BPghqyaYSqwTys36QDg78DdpGsh15LuRDm+rNzvSGdDr5EuZJ6QxdDSd/AhEfFO9tkbSds+Jtu++vnPAtcDL2RVHo1VlzXnbGAm8CLpjGc86ci5KSfQUEXyNqnKY3/gtgrWNZGU7KeRqssW03xVFMDJpG1+h3RA8If6Gdl38xngc6Tv+T/A7tnsP2bDNyU9mo0fQUqsT5O+y/FUVtUFKWH9Kvvcy6Rqsvoz3V8Dm2ff/62NfPYXpN/vHlJS+zXpYrQ1Qw1n8mbVR9K9pAuVhTzduzIkfYN0IbmiI2WzvPiMwKydSBooaeesqmQT0q2YtxQdl5mf6DNrPz1Id8+sT6rquYF0HcCsUK4aMjOrca4aMjOrcVVXNdS/f/8YOnRo0WGYmVWVRx55ZG5EDGhsXtUlgqFDhzJlypSiwzAzqyqSXm5qnquGzMxqnBOBmVmNcyIwM6txTgRmZjXOicDMrMbllggkXSXpDUlTm5gvSRdJmi7pCUnb5RWLmZk1Lc8zgqtJ/Zk2ZR9Sq5XDSH2RXpZjLGZm1oTcniOIiEmShjZTZDRwTdaRyUOS+kkamLVlblZTli2D996DxYth0SJYsACWLIEPPkjDd99N87t0SdNLlsArr8DHPgZLl6ZyH3yQxqdNg3XWgQhYvrzhtWxZ0+PLlqXykIZNjbc0vy0+19T8Z5+FwYOhVysblV6RVnTa4zMrso7NN4fLL2/951pS5ANlg/hw++gzs/c+kggkjSWdNTBkyJB2Cc6sNZYtg1mz4I034NVXYd68tOPq0wcWLkyvBQvS68030/ynnoI11kjjVpmZM4uOoFhLluSz3CITgRp5r9EcGRFXAlcC1NXVuZU8azcRaec9bx5Mn56Gr72WXi++mHbmU6eCtGJHePVJoP7za6wBPXumBNKrF/ToAd27p9esWbD22unVo0d6/ec/sM020K1bKlM/fPll2GQT6No1nUXUv0qny+d16ZLiqI+nqfGW5rfF55qav3Ah9F2B/uHU2N6mA3ymteX7NNVP4EoqMhHM5MN9tg4GZhUUi9WoJUvSDv3552H27DT+yivpqH7mzDR8992WlxMBa66ZqmTWXRf690878nnz0o66d29YffX0WmON9OrTB1ZdNb3Xs+eK7XjM2kKRiWACcJykG4ARwHxfH7A8zZwJDz6YjuKnT4c774S33mr5c716wVprpdfrr8P//A8MGABDh8IGG6Qd/zrrpCNxs2qUWyKQdD0wEugvaSbwfaA7QERcDtwJ7EvqE3cRcFResVjtWb4cpkyBK65IF1+nTk2vxqy9djoq79MH9tsP1lsPhgxJ7w8ZkqoifLRunVmedw0d2sL8AI7Na/1WW2bPhj/9CZ55Bm6/HV544aNlVlkFdtsNttsONtoI1l8ftt023XljVsuqrhlqM4C33047/IceSgmgsbtJevdO1Tef+Qzsuy/stBOstlp7R2rW8TkRWNV49FG4+GKYNCnt+N9/v2Fez55QVwe77ALbbw8jRqR7zs2sZU4E1qH95z/w05/CbbelC7WlttkGdtgB9t8fdt89JQMzaz0nAutwnn4azjoL7rgjXeit17Nnupg7dCicdlq6RdPMVp4TgXUIs2fDTTel1333ffjhrM9/Pt2yefDB6b57M2tbTgRWmCVL4He/gx/9KD3IVa97d/jc5+CAA+Cgg9LTsmaWH/+LWbtbtiw1nHXeeR++22evvWDMmHT071s6zdqPE4G1mwi46660s58/P703cGC6HnDYYb6106woTgSWuzfegHvugZNPbrjzZ/XV4ZhjUrWQq37MiuV/QcvNLbfAtdfCzTc3vNerF5x0Epx+ui/8mnUUTgTW5mbNgv/3/+DGGxveGzEi3ev/7W/7tk+zjsaJwNrMvHlw4onpFtBFi1KVz4EHwve/n9rGN7OOyYnA2sT48XDssel6AMDHPw5//jNsuWWxcZlZy/LsvN5qwDPPwMiR6cj/jTdSU8533ZU6dHESMKsOTgS2QhYuTE/6brllehK4R490RvDSSzBqVOr20Myqg6uGrNWuuQZOOaWhGqhvX3jkEdhww2LjMrMV4+M2q9iSJXDkkfDlL6ck0Ls3/OUvqW8AJwGz6uUzAqvI3LmpZ68ZM9L04YfDVVf5YTCzzsD/xtaiyZPTcwD1JkxIjcKZWefgqiFr1h13wKc+lcY33BCmTXMSMOtsnAisSeeeC5/9LHzwQboe8MQTMGxY0VGZWVtzIrBGXXwxnHFGGt9++9RYnNsGMuucnAjsI666Co4/Po1///vw8MNOAmadmROBfcgZZ8BXvpLGDz889RUgFRqSmeXMdw0ZkDqNOeSQhhZDjzsOLryw2JjMrH04ERiQmoeoTwLf/Cb88pfFxmNm7ceJoMYtXw5nnw2XXZamr7gCxo4tNiYza19OBDVsyRLYeWeYMiVNX3qpk4BZLXIiqFELFsCuu8Ljj6fpceMaLhKbWW1xIqhBzz+fmoqePh369YNbb4Xddis6KjMrihNBjXn/fdhoozTeuzfcf787kDGrdbk+RyBplKTnJE2XdFoj8/tKuk3S45KeknRUnvHUunfegbq6NL7GGjB1qpOAmeV4RiCpK3AJ8BlgJvCwpAkR8XRJsWOBpyPic5IGAM9J+n1ELMkrrlo1dy4MHAhLl6bexG68EdZbr+iozKwjyPOMYDgwPSJeyHbsNwCjy8oE0EeSgN7APGBpjjHVpLlzYZNNUhKA1Iz0HnsUG5OZdRx5JoJBwIyS6ZnZe6UuBjYDZgFPAt+KiOXlC5I0VtIUSVPmzJmTV7ydUgTsuSfMm5e6lHz0Udh776KjMrOOJM9E0FgLNVE2vTfwGLAOsA1wsaTVP/KhiCsjoi4i6gYMGNDWcXZay5bB8OENt4j+85+w7bbFxmRmHU+eiWAmsG7J9GDSkX+po4CbI5kOvAhsmmNMNeWssxoeFjv/fNhii0LDMbMOKs9E8DAwTNL6knoAhwATysq8AuwBIGltYBPghRxjqhkTJ8I556Txa66Bk08uNh4z67hyu2soIpZKOg6YCHQFroqIpyQdk82/HPghcLWkJ0lVSadGxNy8YqoV//pXemAMUmNyX/pSsfGYWcemiPJq+46trq4uptTXd9hHjB+f+hFYvBh22AEefBC6di06KjMrmqRHIqKusXnumKYTueEGOPDAlAR23x3+8Q8nATNrmRNBJ/H+++lMAGDMGLj77vTgmJlZS9zWUCdxyinpgbGNN4bf/ha6+Zc1swr5jKATmD8/9SUAcOqpTgJm1jpOBJ3Al7+cHh7bcUc48siiozGzauNEUOV++Uv405/S+HnnQRf/ombWSt5tVLGpU+GEE9L4CSe4cxkzWzFOBFXqzTdTY3KQnhe44IJCwzGzKuZEUKUOPRRefx3WXx/+/GdQY038mZlVwImgCt1yS9r5A/zhD6l5aTOzFVVxIpC0Wp6BWGVmzYIvfCGNH3JIqhYyM1sZLSYCSZ+U9DTwTDa9taRLc4/MGlV/e+iQIfDrXxcaipl1EpWcEfwfqQOZNwEi4nFg1zyDssZdd11DldDtt8OqqxYbj5l1DhVVDUXEjLK3luUQizVjwYKG5qR/8hP4xCeKjcfMOo9KGiOYIemTQGQdzJxAVk1k7eess1L/w5tuCieeWHQ0ZtaZVHJGcAxwLKnj+ZmkvoW/mWNMVubtt+GKK9L4mWe6VVEza1uVnBFsEhGHlb4haWfgn/mEZKUiYPRoWLQIPvWpdKeQmVlbquSM4JcVvmc5uOACmDQpjf/8535wzMzaXpNnBJJ2Aj4JDJB0Usms1Ul9EFvO3n4bTsq++TPPhOHDCw3HzDqp5qqGegC9szJ9St5fAHwxz6AsOfnkNNxgg3Sx2MwsD00mgoi4D7hP0tUR8XI7xmTAlCnpgbEuXeCSS1wlZGb5qeRi8SJJ5wNbAD3r34yIT+cWVY1btixdGIbU2cyoUcXGY2adWyUXi38PPAusD/wAeAl4OMeYat4VV6TO6Pv3h3Hjio7GzDq7ShLBmhHxa+CDiLgvIo4Gdsw5rpr11ltw2mlp/PjjYbPNio3HzDq/SqqGPsiGsyXtB8wCBucXUm078EB45x1Ye2343veKjsbMakElieAcSX2Bb5OeH1gdODHPoGrV66/DP/6RxusvFJuZ5a3FRBARt2ej84Hd4b9PFlsb++lP07WB9daDffctOhozqxXNPVDWFTiI1MbQ3RExVdJngdOBXsC27RNibZg/Hy6+OI2fc45vFzWz9tPcGcGvgXWBycBFkl4GdgJOi4hb2yG2mhEB++wDS5ZAXR2MGVN0RGZWS5pLBHXAVhGxXFJPYC6wUUS81j6h1Y6TToIHH0zjl1/uawNm1r6a2+UsiYjlABGxGJjW2iQgaZSk5yRNl3RaE2VGSnpM0lOS7mvN8juD115LDctBen5g++0LDcfMalBzZwSbSnoiGxewYTYtICJiq+YWnF1juAT4DKkfg4clTYiIp0vK9AMuBUZFxCuS1lrxTalOv/hFGo4eDWPHFhuLmdWm5hLByj7KNByYHhEvAEi6ARgNPF1SZgxwc0S8AhARb6zkOqvK4sVw/vlp/Oiji43FzGpXc43OrWxDc4OA0r6OZwIjyspsDHSXdC+phdMLI+Ka8gVJGguMBRgyZMhKhtVxfOc7abjaarDXXsXGYma1K8/Lko3dABll092A7YH9gL2B/5W08Uc+FHFlRNRFRN2AAQPaPtICLFkC116bxn/2M+jZs/nyZmZ5qeTJ4hU1k3T7ab3BpOYpysvMjYh3gXclTQK2BqblGFeHcOyxqV2hwYNdLWRmxarojEBSL0mbtHLZDwPDJK0vqQdwCDChrMyfgF0kdZO0Kqnq6JlWrqfqzJjR0KroRRe5M3ozK1aLiUDS54DHgLuz6W0kle/QPyIilgLHARNJO/cbI+IpScdIOiYr80y23CdID66Ni4ipK7gtVeOoo9Jwyy1h//2LjcXMTBHl1fZlBaRHgE8D90bEttl7T7R0+2he6urqYsqUKUWsuk3ceCMcfHAaf/DB1PGMmVneJD0SEXWNzaukamhpRMxv45hqVn3fw0cf7SRgZh1DJReLp0oaA3SVNAw4AXgg37A6p/vug2eege7dU8NyZmYdQSVnBMeT+it+H7iO1Bz1iTnG1ClFwFe+ksa/8x0YOLDYeMzM6lVyRrBJRJwBnJF3MJ3ZzTfD88+nBOCex8ysI6nkjOAXkp6V9ENJW+QeUScUAaefnsbHjPHDY2bWsbSYCCJid2AkMAe4UtKTknxM2wp33QXTpqWmJP73f4uOxszswyp6oCwiXouIi4BjSM8UnJlnUJ3ND36Qhl/8IvTtW2wsZmblKnmgbDNJZ0maClxMumNocO6RdRLz58PkyWm8PiGYmXUklVws/g1wPbBXRJS3FWQt+O1v03DQoNQpvZlZR9NiIogIP/a0gl57Db71rTR+7rnFxmJm1pQmE4GkGyPiIElP8uHmoyvqocxS89KQzgQOP7zYWMzMmtLcGUF2LMtn2yOQzmbBAvj5z9P4BReAGuudwcysA2jyYnFEzM5GvxkRL5e+gG+2T3jV66CD0nDLLVN/xGZmHVUlt49+ppH39mnrQDqT5cth4sQ0/rWv+WzAzDq25q4RfIN05L+BpCdKZvUB/pl3YNXsuuvScPXV4bjjio3FzKwlzV0juA64C/gxcFrJ++9ExLxco6pyV1+dhkcdBV3y7BXazKwNNJcIIiJeknRs+QxJazgZNO6VV+Cvf4Vu3RraFzIz68haOiP4LPAI6fbR0pruADbIMa6qddttaThyJKy1VqGhmJlVpMlEEBGfzYbrt1841a++WuiwwwoNw8ysYpW0NbSzpNWy8S9J+oWkIfmHVn3mzIEpU1K10Be+UHQ0ZmaVqeRS5mXAIklbA98BXgZ+l2tUVWrChDTcZZd0x5CZWTWotPP6AEYDF0bEhaRbSK3Md7+bhgcfXGwcZmatUUnro+9I+i5wOLCLpK5A93zDqj4vvZSqhgAOOKDQUMzMWqWSM4KDSR3XHx0RrwGDgPNzjaoKfe1rabjTTtC/f7GxmJm1RiVdVb4G/B7oK+mzwOKIuCb3yKrI0qXp2QGAs88uNhYzs9aq5K6hg4DJwIHAQcC/JH0x78Cqybnnpg7qe/aEPfYoOhozs9ap5BrBGcAOEfEGgKQBwF+A8XkGVk2uvDIN99vPDcyZWfWp5BpBl/okkHmzws/VhL//HV59NbUpNG5c0dGYmbVeJWcEd0uaSOq3GNLF4zvzC6m6nHNOGn7ta9CvX6GhmJmtkEr6LD5F0heAT5HaG7oyIm7JPbIqsHgx/O1vafzoo4uNxcxsRTXXH8Ew4GfAhsCTwMkR8Wp7BVYN/vjHNFx3Xdhhh2JjMTNbUc3V9V8F3A4cQGqB9JetXbikUZKekzRd0mnNlNtB0rJquxupvnP60aN9kdjMqldzVUN9IuJX2fhzkh5tzYKzJ5AvIXV1ORN4WNKEiHi6kXI/ASa2ZvkdwRNZv21f+lKxcZiZrYzmEkFPSdvS0A9Br9LpiGgpMQwHpkfECwCSbiC1V/R0WbnjgZuAqqpcWby4YXyrrYqLw8xsZTWXCGYDvyiZfq1kOoBPt7DsQcCMkumZwIjSApIGAftny2oyEUgaC4wFGDKkY7SAfc89abjZZtCrV7GxmJmtjOY6ptl9JZfdWK15lE1fAJwaEcvUTCV7RFwJXAlQV1dXvoxC/PjHaeizATOrdpU8R7CiZgLrlkwPBmaVlakDbsiSQH9gX0lLI+LWHONaaTNnwkMPpfH6pqfNzKpVnongYWCYpPWBV4FDgDGlBUq7wZR0NXB7R08CANdnj9Z9/vOw9daFhmJmttJySwQRsVTScaS7gboCV0XEU5KOyeZfnte681bfL/Ho0YWGYWbWJlpMBEr1NocBG0TE2Vl/xR+PiMktfTYi7qSsOYqmEkBEHFlRxAV7/nl4+un03MDeexcdjZnZyquk8bhLgZ2AQ7Ppd0jPB9SksWPTcPRoGDiw2FjMzNpCJVVDIyJiO0n/BoiItyT1yDmuDimioW2hUaOKjcXMrK1UckbwQfb0b8B/+yNYnmtUHdQDD6Rh375w1FHFxmJm1lYqSQQXAbcAa0n6EfAP4Nxco+qg6s8GDjoIetTkOZGZdUaVNEP9e0mPAHuQHhL7fEQ8k3tkHdC116ahLxKbWWdSyV1DQ4BFwG2l70XEK3kG1tHMnQvTpqXxHXcsNhYzs7ZUycXiO0jXBwT0BNYHngO2yDGuDuf229Nw991h0KBiYzEza0uVVA19onRa0nbA13OLqIM6+eQ03HffYuMwM2trre6EPmt+uqqajF5ZEfDmm2nc1UJm1tlUco3gpJLJLsB2wJzcIuqAHi3peWHnnYuLw8wsD5VcI+hTMr6UdM3gpnzC6ZjOOy8NjzjCXVKaWefTbCLIHiTrHRGntFM8HdL48Wm4557FxmFmlocmrxFI6hYRy0hVQTVr3Lg07NIFxoxpvqyZWTVq7oxgMikJPCZpAvBH4N36mRFxc86xdQiTJqXhF74AXbsWG4uZWR4quUawBvAmqV/h+ucJAqiJRPDww2l40knNlzMzq1bNJYK1sjuGptKQAOp1iH6D28Ozz6bhFjX1+JyZ1ZLmEkFXoDeVdULfKd1/fxquuSb06dN8WTOzatVcIpgdEWe3WyQd0H33peEmm/i2UTPrvJp7srjmd3133JGGX/1qsXGYmeWpuUSwR7tF0QEtWAAPPZRuG91//6KjMTPLT5OJICLmtWcgHc0//pGG660H/foVGoqZWa5a3ehcrfjXv9Jwyy2LjcPMLG9OBE24LeuG57DDio3DzCxvTgSNmD0b/v3vNL7rrsXGYmaWNyeCRtyUta36iU/AwIHFxmJmljcngkZMmJCGvlvIzGqBE0GZd99taGjuqKOKjcXMrD04EZT5zW/g/fdhnXVg6NCiozEzy58TQZnrrkvDM88sNg4zs/biRFBi8WJ48ME0fsABxcZiZtZeck0EkkZJek7SdEmnNTL/MElPZK8HJG2dZzwtqU8CAwdC//5FRmJm1n5ySwRZf8eXAPsAmwOHStq8rNiLwG4RsRXwQ+DKvOKpxIwZafj++0VGYWbWvvI8IxgOTI+IFyJiCXADMLq0QEQ8EBFvZZMPAYNzjKdFU6em4dFHFxmFmVn7yjMRDAJmlEzPzN5ryleAuxqbIWmspCmSpsyZM6cNQ/yw889Pw5Ejc1uFmVmHk2ciqLhnM0m7kxLBqY3Nj4grI6IuIuoGDBjQhiE2ePnlhvHdd89lFWZmHVIlndevqJnAuiXTg4FZ5YUkbQWMA/aJiDdzjKdZkyen4Sc/CauuWlQUZmbtL88zgoeBYZLWl9QDOASYUFpA0hDgZuDwiJiWYywtuuyyNHS1kJnVmtzOCCJiqaTjgIlAV+CqiHhK0jHZ/MuBM4E1gUuVOgVeGhF1ecXUnFdfTcPdditi7WZmxVFEo9X2HVZdXV1MmTKlTZc5bx6suWYaf+cd6N27TRdvZlY4SY80daDtJ4uBu+9Ow113dRIws9rjRAC8lT3JoMbuczIz6+ScCGjon3j06ObLmZl1Rk4EwLPPpuGmmxYbh5lZEZwIaLhjaMiQYuMwMytCzSeChQth1izo1g023rjoaMzM2l/NJ4L770/DYcOge/diYzEzK0LNJ4L6W0f33rvYOMzMilLziWDixDQcNarYOMzMilLTieD55+G556BPH7cxZGa1q6YTQf31gbXWglVWKTYWM7Oi1HQiuPnmNPT1ATOrZTWbCObPh9tuS+NHHFFsLGZmRarZRPDYY2m42mowfHihoZiZFapmE8Hf/56GBx/sxubMrLbVbCL4wQ/ScJ99io3DzKxoNZkIFixoGHdH9WZW62oyEYwfn4Zrr93QM5mZWa2qyUQwZ04a7rVXsXGYmXUENZkI7r03DX23kJlZjSaC+obmttyy2DjMzDqCmksEixc3jI8YUVwcZmYdRc0lgjvuSMOttoJevYqNxcysI6i5RHDWWWl48MGFhmFm1mHUXCKYPTsNhw4tNAwzsw6jphLB8uWwaFEa33PPYmMxM+soaioRvPACvPcerLNO6oPAzMxqLBFMm5aGG21UbBxmZh1JTSWCyZPTcKutio3DzKwjqalEUP9E8c47FxqGmVmHUjOJYMkSuO++NO6O6s3MGuSaCCSNkvScpOmSTmtkviRdlM1/QtJ2ecVSnwQ23hg+/vG81mJmVn1ySwSSugKXAPsAmwOHStq8rNg+wLDsNRa4LK947rwzDQcNymsNZmbVKc8zguHA9Ih4ISKWADcAo8vKjAauieQhoJ+kgTnGxLBheS7dzKz65JkIBgEzSqZnZu+1tgySxkqaImnKnPrOBFqpWzdYdVXYddcV+riZWaeVZyJorEv4WIEyRMSVEVEXEXUDBgxYoWDOPx/efRfGjFmhj5uZdVp5JoKZwLol04OBWStQpk2psdRjZlbD8kwEDwPDJK0vqQdwCDChrMwE4Ijs7qEdgfkRMTvHmMzMrEy3vBYcEUslHQdMBLoCV0XEU5KOyeZfDtwJ7AtMBxYBR+UVj5mZNS63RAAQEXeSdval711eMh7AsXnGYGZmzauZJ4vNzKxxTgRmZjXOicDMrMY5EZiZ1Til67XVQ9Ic4OUV/Hh/YG4bhlMNvM21wdtcG1Zmm9eLiEafyK26RLAyJE2JiLqi42hP3uba4G2uDXlts6uGzMxqnBOBmVmNq7VEcGXRARTA21wbvM21IZdtrqlrBGZm9lG1dkZgZmZlnAjMzGpcp0wEkkZJek7SdEmnNTJfki7K5j8habsi4mxLFWzzYdm2PiHpAUlbFxFnW2ppm0vK7SBpmaQvtmd8eahkmyWNlPSYpKck3dfeMba1Cv62+0q6TdLj2TZXdSvGkq6S9IakqU3Mb/v9V0R0qhepyevngQ2AHsDjwOZlZfYF7iL1kLYj8K+i426Hbf4k8LFsfJ9a2OaScn8jtYL7xaLjboffuR/wNDAkm16r6LjbYZtPB36SjQ8A5gE9io59JbZ5V2A7YGoT89t8/9UZzwiGA9Mj4oWIWALcAIwuKzMauCaSh4B+kga2d6BtqMVtjogHIuKtbPIhUm9w1ayS3xngeOAm4I32DC4nlWzzGODmiHgFICKqfbsr2eYA+kgS0JuUCJa2b5htJyImkbahKW2+/+qMiWAQMKNkemb2XmvLVJPWbs9XSEcU1azFbZY0CNgfuJzOoZLfeWPgY5LulfSIpCPaLbp8VLLNFwObkbq5fRL4VkQsb5/wCtHm+69cO6YpSGO9EpffI1tJmWpS8fZI2p2UCD6Va0T5q2SbLwBOjYhl6hydVVeyzd2A7YE9gF7Ag5IeiohpeQeXk0q2eW/gMeDTwIbAnyXdHxELco6tKG2+/+qMiWAmsG7J9GDSkUJry1STirZH0lbAOGCfiHiznWLLSyXbXAfckCWB/sC+kpZGxK3tEmHbq/Rve25EvAu8K2kSsDVQrYmgkm0+CjgvUgX6dEkvApsCk9snxHbX5vuvzlg19DAwTNL6knoAhwATyspMAI7Irr7vCMyPiNntHWgbanGbJQ0BbgYOr+Kjw1ItbnNErB8RQyNiKDAe+GYVJwGo7G/7T8AukrpJWhUYATzTznG2pUq2+RXSGRCS1gY2AV5o1yjbV5vvvzrdGUFELJV0HDCRdMfBVRHxlKRjsvmXk+4g2ReYDiwiHVFUrQq3+UxgTeDS7Ah5aVRxy40VbnOnUsk2R8Qzku4GngCWA+MiotHbEKtBhb/zD4GrJT1JqjY5NSKqtnlqSdcDI4H+kmYC3we6Q377LzcxYWZW4zpj1ZCZmbWCE4GZWY1zIjAzq3FOBGZmNc6JwMysxjkRWIeUtRb6WMlraDNlF7bB+q6W9GK2rkcl7bQCyxgnafNs/PSyeQ+sbIzZcuq/l6lZi5v9Wii/jaR922Ld1nn59lHrkCQtjIjebV22mWVcDdweEeMl7QX8LCK2WonlrXRMLS1X0m+BaRHxo2bKHwnURcRxbR2LdR4+I7CqIKm3pL9mR+tPSvpIS6OSBkqaVHLEvEv2/l6SHsw++0dJLe2gJwEbZZ89KVvWVEknZu+tJumOrP37qZIOzt6/V1KdpPOAXlkcv8/mLcyGfyg9Qs/ORA6Q1FXS+ZIeVmpj/usVfC0PkjU2Jmm4Uj8T/86Gm2RP4p4NHJzFcnAW+1XZev7d2PdoNajotrf98quxF7CM1JDYY8AtpKfgV8/m9Sc9VVl/RrswG34bOCMb7wr0ycpOAlbL3j8VOLOR9V1N1l8BcCDwL1LjbU8Cq5GaN34K2BY4APhVyWf7ZsN7SUff/42ppEx9jPsDv83Ge5BakewFjAW+l72/CjAFWL+ROBeWbN8fgVHZ9OpAt2x8T+CmbPxI4OKSz58LfCkb70dqg2i1on9vv4p9dbomJqzTeC8itqmfkNQdOFfSrqSmEwYBawOvlXzmYeCqrOytEfGYpN2AzYF/Zk1r9CAdSTfmfEnfA+aQWmjdA7glUgNuSLoZ2AW4G/iZpJ+QqpPub8V23QVcJGkVYBQwKSLey6qjtlJDL2p9gWHAi2Wf7yXpMWAo8Ajw55Lyv5U0jNQSZfcm1r8X8D+STs6mewJDqO72iGwlORFYtTiM1PvU9hHxgaSXSDux/4qISVmi2A/4naTzgbeAP0fEoRWs45SIGF8/IWnPxgpFxDRJ25Pae/mxpHsi4uxKNiIiFku6l9R08sHA9fWrA46PiIktLOK9iNhGUl/gduBY4CJSezt/j4j9swvr9zbxeQEHRMRzlcRrtcHXCKxa9AXeyJLA7sB65QUkrZeV+RXwa1J3fw8BO0uqr/NfVdLGFa5zEvD57DOrkap17pe0DrAoIq4Ffpatp9wH2ZlJY24gNRS2C6kxNbLhN+o/I2njbJ2Nioj5wAnAydln+gKvZrOPLCn6DqmKrN5E4Hhlp0eStm1qHVY7nAisWvweqJM0hXR28GwjZUYCj0n6N6ke/8KImEPaMV4v6QlSYti0khVGxKOkaweTSdcMxkXEv4FPAJOzKpozgHMa+fiVwBP1F4vL3EPql/YvkbpfhNRPxNPAo0qdll9BC2fsWSyPk5pm/inp7OSfpOsH9f4ObF5/sZh05tA9i21qNm01zrePmpnVOJ8RmJnVOCcCM7Ma50RgZlbjnAjMzGqcE4GZWY1zIjAzq3FOBGZmNe7/A1azVo3VzQkQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.metrics import roc_curve\n",
    "fpr,tpr,_ = roc_curve(y_true=y_test, y_score=y_pred_test,sample_weight=weights_test)\n",
    "plt.plot(fpr, tpr, color='blue',lw=2)\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('Receiver Operating Characteristic')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "auc test: 0.8746071344141351\n",
      "auc train: 0.8763695135819273\n"
     ]
    }
   ],
   "source": [
    "auc_test = roc_auc_score(y_true=y_test, y_score=y_pred_test,sample_weight=weights_test)\n",
    "auc_train = roc_auc_score(y_true=y_train.values, y_score=y_pred_train,sample_weight=weights_train)\n",
    "print(\"auc test:\",auc_test)\n",
    "print (\"auc train:\",auc_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AIsMSGl-Kwql",
    "tags": []
   },
   "source": [
    "## Significance Function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathrm{med}[Z_0|1] = \\sqrt{q_{0,A}} = \\sqrt{2+((s+b)\\ln(1+s/b)-s)}$\n",
    "\n",
    "**asimov significance [arXiv:1007.1727](https://arxiv.org/pdf/1007.1727.pdf) [Eq. 97]**\n",
    "\n",
    "Likelihood-based statistical test for significance. Need to esimate your sensitivity to MC. Running a toy MC thousands of times, should converge to 'truth'. Asimov is representative of number of sigmas in the gaus case.\n",
    "\n",
    "Essentially: For an observed number of signal events $s$, what is the significance $Z_0$ with which we would reject the $s = 0$ hypothesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "id": "qubY3CMNKwql"
   },
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "from math import log\n",
    "def amsasimov(s,b):\n",
    "        if b<=0 or s<=0:\n",
    "            return 0\n",
    "        try:\n",
    "            return sqrt(2*((s+b)*log(1+float(s)/b)-s))\n",
    "        except ValueError:\n",
    "            print(1+float(s)/b)\n",
    "            print (2*((s+b)*log(1+float(s)/b)-s))\n",
    "        #return s/sqrt(s+b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "id": "iqiOrCfyVYD4"
   },
   "outputs": [],
   "source": [
    "#from extra_functions import amsasimov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z: 0.5916129085239864\n"
     ]
    }
   ],
   "source": [
    "int_pred_test_sig = [weights_test[(y_test ==1) & (y_pred_test > th_cut)].sum() for th_cut in np.linspace(0,1,num=50)]\n",
    "int_pred_test_bkg = [weights_test[(y_test ==0) & (y_pred_test > th_cut)].sum() for th_cut in np.linspace(0,1,num=50)]\n",
    "vamsasimov = [amsasimov(sumsig,sumbkg) for (sumsig,sumbkg) in zip(int_pred_test_sig,int_pred_test_bkg)]\n",
    "Z = max(vamsasimov)\n",
    "print(\"Z:\",Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(np.linspace(0,1,num=50),vamsasimov, label='Significance (Z = {})'.format(np.round(Z,decimals=2)))\n",
    "\n",
    "\n",
    "plt.title(\"NN Significance\")\n",
    "plt.xlabel(\"Threshold\")\n",
    "plt.ylabel(\"Significance\")\n",
    "plt.legend()\n",
    "plt.savefig(\"Significance_xgb.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting NN Score for Signal and Background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_train_test(y_pred_train, y_train, y_pred_test, y_test, \n",
    "                   xlabel=\"NN Score\", title=\"NN\", \n",
    "                   weights_train=weights_train.values, weights_test=weights_test.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What does overtraining look like?\n",
    "\n",
    "Recipe:\n",
    "1. Add More layers\n",
    "2. Add more nodes per layer\n",
    "3. Train on less data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Crazy Example\n",
    "N = len(X_train)\n",
    "n = int(N/1000)\n",
    "print(\"Using\",n,\"/\",N, \"events\")\n",
    "\n",
    "X_small = X_train[:n]\n",
    "y_small = y_train[:n]\n",
    "weights_small = weights_train[:n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ot_model = tf.keras.models.Sequential([\n",
    "    tf.keras.layers.Flatten(input_shape=(X_small.shape[1],)), # input layer\n",
    "    tf.keras.layers.Dense(256, activation='relu'), # 1st hiddden layer\n",
    "    tf.keras.layers.Dense(256, activation='relu'),\n",
    "    tf.keras.layers.Dense(256, activation='relu'),\n",
    "    tf.keras.layers.Dense(1,activation=\"sigmoid\")\n",
    "])\n",
    "ot_model.compile(loss=\"binary_crossentropy\", optimizer=\"adam\")\n",
    "starting_time = time.time( )\n",
    "\n",
    "the_overfit = ot_model.fit(X_small, y_small.values, epochs=25 ,validation_data=(X_val, y_val))\n",
    "#Here, we are not using any sample weights. \n",
    "#This is just an example so we omits weights to show val loss\n",
    "\n",
    "training_time = time.time( ) - starting_time\n",
    "print(\"Training time:\",training_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ot_y_pred_test = ot_model.predict(X_test).ravel()\n",
    "ot_y_pred_train = ot_model.predict(X_small).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from extra_functions import compare_train_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, sharex=False, figsize=(15,5))\n",
    "\n",
    "axes[0].plot(the_overfit.history['loss'],label=\"training loss\")\n",
    "axes[0].plot(the_overfit.history['val_loss'],label=\"validation loss\")\n",
    "axes[0].legend(fontsize=15)\n",
    "\n",
    "axes[1].set_title(\"Overtrained NN\")\n",
    "compare_train_test(ot_y_pred_train, y_small, ot_y_pred_test, y_test, \n",
    "                   xlabel=\"NN Score\", title=\"Overtrained NN\", \n",
    "                   weights_train=weights_small.values, weights_test=weights_test.values)\n",
    "axes[1].legend(loc=\"upper center\",fontsize = 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1JKKJbnWX7Ql"
   },
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IQJuC2-1UkqP"
   },
   "source": [
    "1.   Improve NN AUC and significance by increasing the number of neurons, and layers, epochs, or by any other techniques (google) \n",
    "        - (beware of training time)\n",
    "        - Explore!\n",
    "        - Loss: [BCE, MSE,]\n",
    "        - Activations: [relu, leakyrelu, selu, tanh]\n",
    "2.   Draw NN score for signal and background, training and testing (see overtraining example, or BDT Notebook)\n",
    "        - feel free to look into *extra_functions.py*\n",
    "3.   Add Features more features and engineer aditional Features        \n",
    "___\n",
    "4.   Draw plot looking at significance for NN ( see BDT notebook), for the same features as BDT\n",
    "5.   Enable feature permutation importance (see BDT notebook)\n",
    "6.   Calculate the permutation importance with respect to significance and accuracy. See how Asimov significance changes as features shuffle.\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Some help for BDT comparisons:\n",
    "y_pred_xgb = np.load(\"y_pred_xgb.npy\")\n",
    "y_pred_train_xgb = np.load(\"./y_pred_train_xgb.npy\")\n",
    "weights_train_xgb = np.load(\"./weights_train_xgb.npy\")\n",
    "y_test_xgb = np.load(\"./y_test_xgb.npy\")\n",
    "weights_test_xgb = np.load(\"./weights_test_xgb.npy\")\n",
    "\n",
    "from xgboost import XGBClassifier\n",
    "import xgboost\n",
    "xgb = xgboost.XGBClassifier({'nthread': 4})  # init model\n",
    "xgb.load_model('xgb.model')  # load data\n",
    "\n",
    "#The model and predictions from the BDT notebook are now in memory!\n",
    "#plt.bar(data.columns.values, xgb.feature_importances_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extra: Unweighted Score Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare_train_test(y_pred_train, y_train, y_pred_test, y_test, \n",
    "                   xlabel=\"NN Score\", title=\"NN\", \n",
    "                   weights_train=weights_train.values,\n",
    "                   weights_test=weights_test.values,density=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "name": "HEPML_HandsOn_NN_veryold.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "ml_tutorial",
   "language": "python",
   "name": "ml_tutorial"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}