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": [
    "Please make sure you have enabled a GPU (Modify in the menu above) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eT7-MMpfrlHR"
   },
   "source": [
    "## Standard Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "QdmaHjz-xCm-"
   },
   "outputs": [],
   "source": [
    "#COLAB=False #if running on local anaconda installation https://docs.anaconda.com/anaconda/install/\n",
    "COLAB=True #if running on https://colab.research.google.com/notebooks/welcome.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
    "#from IPython import display\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": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "TvU6crowVYDO",
    "outputId": "4ff001e1-9a76-4894-f139-2c8429dba23f"
   },
   "outputs": [],
   "source": [
    "# Fix environment if necessary\n",
    "if COLAB:\n",
    "    # tensorflow already installed\n",
    "    pass\n",
    "else:    \n",
    "    #!pip install tensorflow # uncomment to install tensorflow\n",
    "    pass\n",
    "\n",
    "import tensorflow as tf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "IsChk6u-r7fk",
    "outputId": "7fe98ed4-589c-4483-b2d4-573194e5dcc5",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "if COLAB:\n",
    "    #### Reading files from Google Drive\n",
    "    # one need a google account to be identified\n",
    "    # select a google account, then cut and paste the long password in the pop up field\n",
    "    !pip install PyDrive\n",
    "    import os\n",
    "    from pydrive.auth import GoogleAuth\n",
    "    from pydrive.drive import GoogleDrive\n",
    "    from google.colab import auth\n",
    "    from oauth2client.client import GoogleCredentials\n",
    "    auth.authenticate_user()\n",
    "    gauth = GoogleAuth()\n",
    "    gauth.credentials = GoogleCredentials.get_application_default()\n",
    "    drive = GoogleDrive(gauth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dHSzcuTgVYDT"
   },
   "source": [
    "# Load events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kl-4W4Ifs8EV",
    "outputId": "54d366bc-8925-4c6a-e867-ef70016e016b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running .zshenv\n",
      "total 3104\n",
      "-rw-r--r--@ 1 rousseau  staff   40958 12 jan 11:22 Copy_of_Welcome_To_Colaboratory.ipynb\n",
      "-rw-r--r--+ 1 rousseau  staff    3881 13 jan 15:27 requirements.txt\n",
      "-rw-r--r--@ 1 rousseau  staff   63621 20 jan 10:54 SOS2021_practice_Basic_statistics_from_github.ipynb\n",
      "-rw-------@ 1 rousseau  staff   18810 24 jan 12:37 HEPML_HandsOn_install_check.ipynb\n",
      "-rw-r--r--+ 1 rousseau  staff    6475 24 jan 20:40 html_original.html\n",
      "-rw-r--r--@ 1 rousseau  staff  188022 24 jan 21:46 Notebook_separation_maury_permutation_importance.ipynb\n",
      "-rw-r--r--+ 1 rousseau  staff  143607 24 jan 23:25 XGBoost.model\n",
      "-rw-r--r--+ 1 rousseau  staff  361110 24 jan 23:25 LightGBM.model\n",
      "-rw-r--r--+ 1 rousseau  staff   16228 24 jan 23:25 Score_BDT_XGBoost_Hist.pdf\n",
      "-rw-r--r--+ 1 rousseau  staff   17086 24 jan 23:25 Score_BDT_LightGBM.pdf\n",
      "-rw-r--r--+ 1 rousseau  staff   48661 24 jan 23:25 ROC_comparing.pdf\n",
      "-rw-r--r--+ 1 rousseau  staff   15438 24 jan 23:25 Significance_comparing.pdf\n",
      "-rw-------@ 1 rousseau  staff  339507 24 jan 23:26 HEPML_HandsOn_BDT.ipynb\n",
      "-rw-r--r--@ 1 rousseau  staff  103055 24 jan 23:43 HEPML_HandsOn_NN_veryold_immersion_save.ipynb\n",
      "-rw-------@ 1 rousseau  staff   30321 24 jan 23:44 HEPML_HandsOn_NN.ipynb\n",
      "-rw-r--r--@ 1 rousseau  staff   89534 25 jan 00:05 HEPML_HandsOn_NN_veryold_immersion.ipynb\n",
      "running .zshenv\n",
      "total 70344\n",
      "-rwxr-xr-x+ 1 rousseau  staff  35905611 13 mai  2020 \u001b[31mdataWW_d1_600k.csv.gz\u001b[m\u001b[m\n",
      "Will take data from :  /Users/rousseau/Google Drive/GD_openData/dataWW_ATLAS_openData13TeV_filtered\n",
      "File loaded with  600000  events \n"
     ]
    }
   ],
   "source": [
    "if COLAB:\n",
    "\n",
    "    #attach dataset from google drive \n",
    "    download = drive.CreateFile({'id': '1nlXp7P-xq_jip4aPE0j0mnPhYnIOcBv4'})\n",
    "    download.GetContentFile(\"dataWW_d1_600k.csv.gz\")\n",
    "\n",
    "\n",
    "    datapath=\"\"\n",
    "\n",
    "\n",
    "    !ls -lrt\n",
    "else :\n",
    "    # make sure the file is available locally. \n",
    "    #Should be downloaded from https://drive.google.com/open?id=1nlXp7P-xq_jip4aPE0j0mnPhYnIOcBv4\n",
    "    !ls -lrt # what is in the local directory\n",
    "    datapath=\"/Users/rousseau/Google\\ Drive/GD_openData/dataWW_ATLAS_openData13TeV_filtered/\"\n",
    "\n",
    "    !ls -lrt {datapath} # what is in the data directory\n",
    "    datapath=os.path.abspath(datapath).replace(\"\\ \", \" \")  # try to normalise the path (annoyance with the space)\n",
    "    print (\"Will take data from : \",datapath)\n",
    "\n",
    "filename=os.path.join(datapath,\"dataWW_d1_600k.csv.gz\")\n",
    "#load data\n",
    "# data was created from ATLAS Open Data see doc\n",
    "# http://opendata.atlas.cern/release/2020/documentation/datasets/intro.html\n",
    "dfall = pd.read_csv(filename) \n",
    "\n",
    "#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 \")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7xYM7O_zVYDY"
   },
   "source": [
    "At this point, it should tell you \"File Loaded with XXX events\". If not, it could not access the datafile, no point going further !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LFPm8OHuVYDY"
   },
   "source": [
    "# Examine dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "xNeCAIuFVYDZ",
    "outputId": "b12cbd30-536e-429c-a1be-677ee0b937fb"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['eventNumber', 'label', 'met_et', 'met_phi', 'lep_n', 'lep_pt_0',\n",
       "       'lep_pt_1', 'lep_eta_0', 'lep_eta_1', 'lep_phi_0', 'lep_phi_1',\n",
       "       'lep_E_0', 'lep_E_1', 'lep_charge_0', 'lep_charge_1', 'lep_type_0',\n",
       "       'lep_type_1', 'jet_n', 'jet_pt_0', 'jet_pt_1', 'jet_eta_0', 'jet_eta_1',\n",
       "       'jet_phi_0', 'jet_phi_1', 'jet_E_0', 'jet_E_1', 'mcWeight', 'runNumber',\n",
       "       'channelNumber'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#dump list of feature\n",
    "dfall.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>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>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>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>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>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>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": [
       "   eventNumber  label   met_et  met_phi  lep_n  lep_pt_0  lep_pt_1  lep_eta_0  \\\n",
       "0       402756      1   25.609  0.42452      2    48.295    15.214    0.73991   \n",
       "1       101274      0  196.560  1.31140      2    69.459    21.081   -0.52666   \n",
       "2       468437      1   45.653 -2.76860      2    45.927    22.822   -1.61910   \n",
       "3       272337      1   49.415 -0.57805      2    45.929    14.263    1.55280   \n",
       "4       354546      1   71.988 -2.60390      2    62.029    21.453   -0.51082   \n",
       "\n",
       "   lep_eta_1  lep_phi_0  lep_phi_1   lep_E_0  lep_E_1  lep_charge_0  \\\n",
       "0    2.27420  -2.316400   -1.39410   62129.0  74721.0            -1   \n",
       "1    0.22380   0.023132   -0.67855   79317.0  21611.0            -1   \n",
       "2   -2.00770   1.901700   -0.11248  120480.0  86498.0            -1   \n",
       "3    0.12809   2.254700   -2.54810  113360.0  14380.0            -1   \n",
       "4   -0.38177   0.236200    1.40820   70300.0  23036.0            -1   \n",
       "\n",
       "   lep_charge_1  lep_type_0  lep_type_1  jet_n  jet_pt_0  jet_pt_1  jet_eta_0  \\\n",
       "0             1          13          11      2     26.32    20.064    -1.1350   \n",
       "1             1          13          13      2    192.12    36.217    -1.2643   \n",
       "2             1          11          11      0     -7.00    -7.000    -7.0000   \n",
       "3             1          13          13      0     -7.00    -7.000    -7.0000   \n",
       "4             1          11          13      0     -7.00    -7.000    -7.0000   \n",
       "\n",
       "   jet_eta_1  jet_phi_0  jet_phi_1   jet_E_0  jet_E_1  mcWeight  runNumber  \\\n",
       "0   -2.07540    -2.8092     2.0740   45397.0  81304.0  0.000002     284500   \n",
       "1   -0.69746    -2.2147    -1.9152  367670.0  45733.0  0.000288     284500   \n",
       "2   -7.00000    -7.0000    -7.0000      -7.0     -7.0  0.000029     284500   \n",
       "3   -7.00000    -7.0000    -7.0000      -7.0     -7.0  0.000029     284500   \n",
       "4   -7.00000    -7.0000    -7.0000      -7.0     -7.0  0.000002     284500   \n",
       "\n",
       "   channelNumber  \n",
       "0         345323  \n",
       "1         363492  \n",
       "2         345324  \n",
       "3         345324  \n",
       "4         345323  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#examine first few events\n",
    "dfall.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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>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>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>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>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>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>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>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>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>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": [
       "        eventNumber          label         met_et        met_phi  \\\n",
       "count  6.000000e+05  600000.000000  600000.000000  600000.000000   \n",
       "mean   7.029775e+05       0.666667      60.365356      -0.008556   \n",
       "std    4.568961e+05       0.471405      48.153210       1.812778   \n",
       "min    1.000000e+00       0.000000       0.044211      -3.141600   \n",
       "25%    3.498568e+05       0.000000      33.196000      -1.577600   \n",
       "50%    6.678815e+05       1.000000      50.082000      -0.018694   \n",
       "75%    1.001273e+06       1.000000      73.032250       1.559300   \n",
       "max    1.997995e+06       1.000000    7062.200000       3.141600   \n",
       "\n",
       "               lep_n       lep_pt_0       lep_pt_1      lep_eta_0  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean        2.005748      60.270328      28.133970       0.000269   \n",
       "std         0.076061      42.676858      20.352272       1.210614   \n",
       "min         2.000000      25.000000       7.000300      -2.700000   \n",
       "25%         2.000000      36.624750      15.316000      -0.899650   \n",
       "50%         2.000000      48.229500      23.644000       0.003413   \n",
       "75%         2.000000      68.979000      34.149000       0.901950   \n",
       "max         5.000000    7041.400000     744.900000       2.699800   \n",
       "\n",
       "           lep_eta_1      lep_phi_0      lep_phi_1       lep_E_0  \\\n",
       "count  600000.000000  600000.000000  600000.000000  6.000000e+05   \n",
       "mean       -0.000664       0.005219       0.008524  1.164008e+05   \n",
       "std         1.235890       1.811177       1.814205  1.083341e+05   \n",
       "min        -2.699900      -3.141600      -3.141600  2.501300e+04   \n",
       "25%        -0.946043      -1.560200      -1.562900  5.221675e+04   \n",
       "50%         0.002974       0.014583       0.018823  8.184450e+04   \n",
       "75%         0.942210       1.573000       1.577300  1.412400e+05   \n",
       "max         2.700000       3.141500       3.141600  7.041600e+06   \n",
       "\n",
       "            lep_E_1   lep_charge_0   lep_charge_1     lep_type_0  \\\n",
       "count  6.000000e+05  600000.000000  600000.000000  600000.000000   \n",
       "mean   5.630948e+04      -0.027027       0.027037      11.907273   \n",
       "std    5.748834e+04       0.999636       0.999635       0.995692   \n",
       "min    7.003800e+03      -1.000000      -1.000000      11.000000   \n",
       "25%    2.329100e+04      -1.000000      -1.000000      11.000000   \n",
       "50%    3.824900e+04      -1.000000       1.000000      11.000000   \n",
       "75%    6.730825e+04       1.000000       1.000000      13.000000   \n",
       "max    2.107900e+06       1.000000       1.000000      13.000000   \n",
       "\n",
       "          lep_type_1          jet_n       jet_pt_0       jet_pt_1  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean       11.958010       1.334185      59.606114      13.521164   \n",
       "std         0.999119       1.232632      79.586143      38.074439   \n",
       "min        11.000000       0.000000      -7.000000      -7.000000   \n",
       "25%        11.000000       0.000000      -7.000000      -7.000000   \n",
       "50%        11.000000       1.000000      38.744000      -7.000000   \n",
       "75%        13.000000       2.000000      84.291250      26.881000   \n",
       "max        13.000000      13.000000    1658.200000    1429.800000   \n",
       "\n",
       "           jet_eta_0      jet_eta_1      jet_phi_0      jet_phi_1  \\\n",
       "count  600000.000000  600000.000000  600000.000000  600000.000000   \n",
       "mean       -1.934521      -4.402337      -1.939645      -4.401941   \n",
       "std         3.349404       3.491304       3.489197       3.559421   \n",
       "min        -7.000000      -7.000000      -7.000000      -7.000000   \n",
       "25%        -7.000000      -7.000000      -7.000000      -7.000000   \n",
       "50%        -0.891010      -7.000000      -1.210600      -7.000000   \n",
       "75%         0.724705      -0.821005       0.967753      -1.088575   \n",
       "max         2.500000       2.500000       3.141600       3.141600   \n",
       "\n",
       "            jet_E_0       jet_E_1       mcWeight  runNumber  channelNumber  \n",
       "count  6.000000e+05  6.000000e+05  600000.000000   600000.0  600000.000000  \n",
       "mean   1.375410e+05  4.125468e+04       0.000218   284500.0  351379.666667  \n",
       "std    1.980717e+05  8.990068e+04       0.000736        0.0    8564.720183  \n",
       "min   -7.000000e+00 -7.000000e+00      -0.045557   284500.0  345323.000000  \n",
       "25%   -7.000000e+00 -7.000000e+00       0.000002   284500.0  345323.000000  \n",
       "50%    7.501700e+04 -7.000000e+00       0.000029   284500.0  345324.000000  \n",
       "75%    1.766200e+05  5.060500e+04       0.000197   284500.0  363492.000000  \n",
       "max    5.638800e+06  2.784300e+06       0.052420   284500.0  363492.000000  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#examine feature distribution\n",
    "dfall.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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",
    "\n",
    "label_nevents = (dfall[dfall.label==0].shape[0], dfall[dfall.label==1].shape[0] )\n",
    "print (\"total class number of events\",label_nevents)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rtI5u5GErlHq"
   },
   "source": [
    "## Event selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kaO2JM1hrlHr",
    "outputId": "bc59e98d-f6f1-42b7-e9f3-48052fe4ebc9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(600000, 29)\n",
      "(577357, 29)\n"
     ]
    }
   ],
   "source": [
    "print (dfall.shape)\n",
    "fulldata=dfall[(dfall.lep_n==2) & (dfall.mcWeight>0 )] # only keep events with exactly two leptons \n",
    "print (fulldata.shape)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iMQpKhDKrlH0",
    "outputId": "ccb36177-999c-4203-939e-4e04383c1954"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(577357, 27)"
      ]
     },
     "execution_count": 11,
     "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": {
    "id": "nviyIMgerlH3"
   },
   "source": [
    "\n",
    "# DO NOT MODIFY ANYTHING ABOVE\n",
    "... and always rerun from this cell whenever you change something below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# for simplicity of the exercise only keep some features\n",
    "# this is actually making a copy from fulldata\n",
    "data=pd.DataFrame(fulldata, columns=[\"met_et\",\"met_phi\",\"lep_pt_0\",\"lep_pt_1\",'lep_phi_0', 'lep_phi_1'])\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']\n",
    "print (data.shape)\n",
    "data.head()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Krg3pTRjVYDs"
   },
   "source": [
    "# Feature engineering\n",
    "To be switched on in a second iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "id": "DH5UqBF9VYDt"
   },
   "outputs": [],
   "source": [
    "if False: \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",
    "    #data[\"lep_deltaphi\"]=data.lep_phi_1-data.lep_phi_0\n",
    "\n",
    "\n",
    "    print (data.shape)\n",
    "    display(data.head())\n",
    "\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "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": "iVBORw0KGgoAAAANSUhEUgAAA3cAAAK7CAYAAABRbnZtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABlLklEQVR4nO39f7xld13ffb/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": [
    "\n",
    "\n",
    "plt.figure()\n",
    "\n",
    "#_, ax = plt.subplots()\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",
    "\n",
    "\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4kowHjX4rlIC"
   },
   "source": [
    "## Transformation of the features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Y9j5hdrmrlID",
    "outputId": "75c0925f-74c9-4e24-d169-8a884b510f4c"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(433017, 6)\n",
      "(433017,)\n",
      "(433017,)\n",
      "(144340, 6)\n",
      "(144340,)\n",
      "(144340,)\n",
      "class_weights_train: (119.66325440279999, 4.634016907200001)\n",
      "Test : total weight sig 6.1815253696\n",
      "Test : total weight bkg 160.1009499872\n",
      "Train : total weight sig 119.66325440279994\n",
      "Train : total weight bkg 119.66325440279999\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",
    "#reset index for dataseries, not needed for ndarray (X_train, X_test)\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",
    "print (X_train.shape)\n",
    "print (y_train.shape)\n",
    "print (weights_train.shape)\n",
    "print (X_test.shape)\n",
    "print (y_test.shape)\n",
    "print (weights_test.shape)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "#scale to mean 0 and variance 1\n",
    "scaler = StandardScaler()\n",
    "X_train = scaler.fit_transform(X_train)\n",
    "X_test = scaler.transform(X_test)  # applied the transformation calculated the line above\n",
    "\n",
    "\n",
    "class_weights_train = (weights_train[y_train == 0].sum(), weights_train[y_train == 1].sum())\n",
    "\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 (\"Test : total weight sig\", weights_test[y_test == 1].sum())\n",
    "print (\"Test : total weight bkg\", weights_test[y_test == 0].sum())\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",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yxybCOi-rlIM"
   },
   "source": [
    "# Testing NN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "q2Dw40gRVYD3"
   },
   "source": [
    "## Load significance function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "id": "iqiOrCfyVYD4"
   },
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "from math import log\n",
    "def amsasimov(s,b): # asimov (or Poisson) significance\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)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_TjIniZPWlVD"
   },
   "source": [
    "## NN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "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": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/10\n",
      "13532/13532 - 8s - loss: 0.4312\n",
      "Epoch 2/10\n",
      "13532/13532 - 7s - loss: 0.3870\n",
      "Epoch 3/10\n",
      "13532/13532 - 7s - loss: 0.3776\n",
      "Epoch 4/10\n",
      "13532/13532 - 7s - loss: 0.3731\n",
      "Epoch 5/10\n",
      "13532/13532 - 7s - loss: 0.3708\n",
      "Epoch 6/10\n",
      "13532/13532 - 8s - loss: 0.3693\n",
      "Epoch 7/10\n",
      "13532/13532 - 7s - loss: 0.3684\n",
      "Epoch 8/10\n",
      "13532/13532 - 7s - loss: 0.3676\n",
      "Epoch 9/10\n",
      "13532/13532 - 8s - loss: 0.3670\n",
      "Epoch 10/10\n",
      "13532/13532 - 8s - loss: 0.3662\n",
      "Training time: 75.30642509460449\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 hidden layer\n",
    "  tf.keras.layers.Dense(1,activation=\"sigmoid\") # output layer\n",
    "])\n",
    "model.compile(loss=\"binary_crossentropy\", optimizer=\"adam\")\n",
    "starting_time = time.time( )\n",
    "model.fit(X_train, y_train.values, epochs=10, verbose=2)#, class_weight=class_weights)\n",
    "training_time = time.time( ) - starting_time\n",
    "print(\"Training time:\",training_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "auc test: 0.8724323841062293\n",
      "auc train: 0.8721320565898362\n",
      "Z: 0.8159291774347175\n"
     ]
    }
   ],
   "source": [
    "y_pred_test = model.predict(X_test).ravel()\n",
    "y_pred_train = model.predict(X_train).ravel()\n",
    "auc_test = roc_auc_score(y_true=y_test, y_score=y_pred_test)\n",
    "auc_train = roc_auc_score(y_true=y_train.values, y_score=y_pred_train)\n",
    "print(\"auc test:\",auc_test)\n",
    "print (\"auc train:\",auc_train)\n",
    "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",
    "\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": "markdown",
   "metadata": {
    "id": "1JKKJbnWX7Ql"
   },
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "IQJuC2-1UkqP"
   },
   "source": [
    "*   improve NN AUC and significance by increasing the number of neurones, and layers, epochs, or by any other techniques (google) (beware of training time)\n",
    "*   draw NN score for signal and background, training and testing (see BDT notebook) \n",
    "*   enable feature permutation importance (see BDT notebook)\n",
    "*   draw plot comparing significance for NN vs BDT ( see BDT notebook), for the same features\n",
    "\n"
   ]
  },
  {
   "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": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}