{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import pickle\n",
    "import math\n",
    "\n",
    "from xgboost import XGBClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import accuracy_score, roc_auc_score\n",
    "from architectures.data_processing import *\n",
    "from architectures.utils.toolbox import *\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IMPORTING THE DATASET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "l_index=1\n",
    "mother_ID=[\"Ds\",\"Dplus\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bkg data amounts to 13821 while signal MC amounts to 5881 Ds and 9639 Dplus samples\n"
     ]
    }
   ],
   "source": [
    "MC_Dplus_sig_dict, MC_Ds_sig_dict, data_bkg_dict = load_datasets(l_index)\n",
    "\n",
    "m_plus=MC_Dplus_sig_dict[\"Dplus_ConsD_M\"].shape[0]\n",
    "m_s=MC_Ds_sig_dict[\"Ds_ConsD_M\"].shape[0]\n",
    "n=data_bkg_dict[\"Ds_ConsD_M\"].shape[0]\n",
    "\n",
    "print('Bkg data amounts to {0} while signal MC amounts to {1} Ds and {2} Dplus samples'.format(n,m_s,m_plus))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Normalising the Chi2 vertex fits to the NDoF\n",
    "\n",
    "MC_Ds_sig_dict[\"Ds_ENDVERTEX_CHI2\"]=MC_Ds_sig_dict[\"Ds_ENDVERTEX_CHI2\"]/MC_Ds_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "MC_Ds_sig_dict[\"Ds_OWNPV_CHI2\"]=MC_Ds_sig_dict[\"Ds_OWNPV_CHI2\"]/MC_Ds_sig_dict[\"Ds_OWNPV_NDOF\"]\n",
    "MC_Ds_sig_dict[\"Ds_IPCHI2_OWNPV\"]=MC_Ds_sig_dict[\"Ds_IPCHI2_OWNPV\"]/MC_Ds_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "\n",
    "del MC_Ds_sig_dict[\"Ds_OWNPV_NDOF\"]\n",
    "del MC_Ds_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "#del MC_sig_dict[\"Ds_M\"]\n",
    "\n",
    "MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_CHI2\"]=MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_CHI2\"]/MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_NDOF\"]\n",
    "MC_Dplus_sig_dict[\"Dplus_OWNPV_CHI2\"]=MC_Dplus_sig_dict[\"Dplus_OWNPV_CHI2\"]/MC_Dplus_sig_dict[\"Dplus_OWNPV_NDOF\"]\n",
    "MC_Dplus_sig_dict[\"Dplus_IPCHI2_OWNPV\"]=MC_Dplus_sig_dict[\"Dplus_IPCHI2_OWNPV\"]/MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_NDOF\"]\n",
    "\n",
    "del MC_Dplus_sig_dict[\"Dplus_OWNPV_NDOF\"]\n",
    "del MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_NDOF\"]\n",
    "\n",
    "data_bkg_dict[\"Ds_ENDVERTEX_CHI2\"]=data_bkg_dict[\"Ds_ENDVERTEX_CHI2\"]/data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "data_bkg_dict[\"Ds_OWNPV_CHI2\"]=data_bkg_dict[\"Ds_OWNPV_CHI2\"]/data_bkg_dict[\"Ds_OWNPV_NDOF\"]\n",
    "data_bkg_dict[\"Ds_IPCHI2_OWNPV\"]=data_bkg_dict[\"Ds_IPCHI2_OWNPV\"]/data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "\n",
    "del data_bkg_dict[\"Ds_OWNPV_NDOF\"]\n",
    "del data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
    "\n",
    "data_bkg_dict[\"phi_ENDVERTEX_CHI2\"]=data_bkg_dict[\"phi_ENDVERTEX_CHI2\"]/data_bkg_dict[\"phi_ENDVERTEX_NDOF\"]\n",
    "#data_bkg_dict[\"phi_OWNPV_CHI2\"]=data_bkg_dict[\"phi_OWNPV_CHI2\"]/data_bkg_dict[\"phi_OWNPV_NDOF\"]\n",
    "data_bkg_dict[\"phi_IPCHI2_OWNPV\"]=data_bkg_dict[\"phi_IPCHI2_OWNPV\"]/data_bkg_dict[\"phi_ENDVERTEX_NDOF\"]\n",
    "\n",
    "#del data_bkg_dict[\"phi_OWNPV_NDOF\"]\n",
    "del data_bkg_dict[\"phi_ENDVERTEX_NDOF\"]\n",
    "del MC_Ds_sig_dict[\"phi_M\"]\n",
    "del MC_Dplus_sig_dict[\"phi_M\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def return_branches(mother_index=None):\n",
    "\n",
    "    branches_needed = [\n",
    "                    #________________________________________\n",
    "                    #D\n",
    "                    #________________________________________\n",
    "                    #D Geometric variables, pT and FD\n",
    "        \n",
    "                    mother_ID[mother_index]+\"_ENDVERTEX_CHI2\",\n",
    "                    mother_ID[mother_index]+\"_IPCHI2_OWNPV\",\n",
    "\n",
    "                    mother_ID[mother_index]+\"_OWNPV_CHI2\",\n",
    "                    mother_ID[mother_index]+\"_IP_OWNPV\",\n",
    "                    mother_ID[mother_index]+\"_DIRA_OWNPV\",\n",
    "        \n",
    "                    mother_ID[mother_index]+\"_PT\",\n",
    "                    mother_ID[mother_index]+\"_FD_OWNPV\",\n",
    "                    \n",
    "                    #________________________________________\n",
    "                    #PHI\n",
    "                    #________________________________________\n",
    "                    #phi geometric variables, pT and FD\n",
    "        \n",
    "                    \"phi_ENDVERTEX_CHI2\",\n",
    "                    \"phi_IPCHI2_OWNPV\",\n",
    "                    \n",
    "                    \"phi_PT\",\n",
    "        \n",
    "                    #phi Reconstructed mass\n",
    "        \n",
    "                    #\"phi_M\",\n",
    "        \n",
    "                    #________________________________________\n",
    "                    #PION\n",
    "                    #________________________________________\n",
    "                    #pi Geometric variables and pT\n",
    "                    'pi_PT',\n",
    "        \n",
    "                    #________________________________________\n",
    "                    #LEPTONS\n",
    "                    #________________________________________\n",
    "                    #leptons Geometric variables and pT\n",
    "        \n",
    "                    l_flv[l_index]+\"_plus_PT\",\n",
    "                    l_flv[l_index]+\"_minus_PT\",\n",
    "        \n",
    "        \n",
    "                    #D Reconstructed mass\n",
    "                    mother_ID[mother_index]+\"_ConsD_M\",\n",
    "                  ] \n",
    "    return branches_needed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Number of input features\n",
    "\n",
    "dim=len(return_branches(mother_index=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert data dictionaries to arrays for NN\n",
    "\n",
    "MC_Ds_sig = extract_array(MC_Ds_sig_dict, return_branches(mother_index=0), dim, m_s)\n",
    "MC_Dplus_sig = extract_array(MC_Dplus_sig_dict, return_branches(mother_index=1), dim, m_plus)\n",
    "data_bkg = extract_array(data_bkg_dict, return_branches(mother_index=0), dim, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA64AAAHVCAYAAADxfKZZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHhJJREFUeJzt3XGwpfVd3/HPN6yJVaOQsIkIpEsV0kmdqnGHxDq1QYQQ6oRMK0rGMWuk3WklVkkdA9qWmTgOcbRi0lEsFZTMZEJoTCfbFhu3cW2mMxJZMCYhJLCNFlYwrCWibaop+u0f99lws3v33t177t7zu+e+XjN37jm/57n3/C48e3bf9/ec51R3BwAAAEb1nHlPAAAAAFYjXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhrZj3hNYzdlnn927du2a9zQAAAA4De6///4/7u6da+03dLju2rUrBw8enPc0AAAAOA2q6n+ezH5OFQYAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhrZj3hMAYIs5cPPxY5fcuPnzAAC2DSuuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQ9sx7wkAMLgDN897BgDANmfFFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIa2ZrhW1R1V9WRVfXyFbT9aVV1VZ0/3q6reUVWHquqjVfXyZfvuqapHpo89G/tjAAAAsKhOZsX1V5NccexgVZ2f5LIkjy4bfk2SC6ePvUlunfZ9QZKbkrwiycVJbqqqs2aZOAAAANvDmuHa3R9K8tQKm25J8mNJetnYVUne2UvuTXJmVZ2T5NVJ9nf3U9392ST7s0IMAwAAwLHW9RrXqnptkj/s7t87ZtO5SR5bdv/wNHai8ZW+996qOlhVB48cObKe6QEAALBATjlcq+rLkvxEkn+10uYVxnqV8eMHu2/r7t3dvXvnzp2nOj0AAAAWzHpWXL82yQVJfq+q/iDJeUkeqKqvztJK6vnL9j0vyeOrjAMAAMCqTjlcu/tj3f2i7t7V3buyFKUv7+4/SrIvyRumqwu/MsnT3f1Ekg8kubyqzpouynT5NAYAAACrOpm3w3l3kt9O8tKqOlxV166y+z1JPp3kUJJ/l+QHk6S7n0ryk0numz7eOo0BAADAqnastUN3v36N7buW3e4k151gvzuS3HGK8wMAAGCbW9dVhQEAAGCzCFcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABjajnlPAIAFcODm48cuuXHz5wEALCQrrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADG3NcK2qO6rqyar6+LKxn6mqT1bVR6vqP1TVmcu23VhVh6rqU1X16mXjV0xjh6rqho3/UQAAAFhEJ7Pi+qtJrjhmbH+Sr+/uv53k4SQ3JklVvSzJNUn+1vQ1v1hVZ1TVGUl+IclrkrwsyeunfQEAAGBVa4Zrd38oyVPHjP1Gdz8z3b03yXnT7auS3NXdf9Hdv5/kUJKLp49D3f3p7v58krumfQEAAGBVG/Ea1x9I8uvT7XOTPLZs2+Fp7ETjAAAAsKqZwrWqfiLJM0nedXRohd16lfGVvufeqjpYVQePHDkyy/QAAABYAOsO16rak+Q7k3xvdx+N0MNJzl+223lJHl9l/DjdfVt37+7u3Tt37lzv9AAAAFgQ6wrXqroiyVuSvLa7P7ds074k11TV86rqgiQXJvmdJPclubCqLqiq52bpAk77Zps6AAAA28GOtXaoqncneVWSs6vqcJKbsnQV4ecl2V9VSXJvd/+T7n6wqu5O8oksnUJ8XXf/5fR93pTkA0nOSHJHdz94Gn4eAAAAFsya4drdr19h+PZV9v+pJD+1wvg9Se45pdkBAACw7W3EVYUBAADgtBGuAAAADG3NU4UB2EYO3DzvGQAAHMeKKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwtDXDtaruqKonq+rjy8ZeUFX7q+qR6fNZ03hV1Tuq6lBVfbSqXr7sa/ZM+z9SVXtOz48DAADAojmZFddfTXLFMWM3JPlgd1+Y5IPT/SR5TZILp4+9SW5NlkI3yU1JXpHk4iQ3HY1dAAAAWM2a4drdH0ry1DHDVyW5c7p9Z5LXLRt/Zy+5N8mZVXVOklcn2d/dT3X3Z5Psz/ExDAAAAMdZ72tcX9zdTyTJ9PlF0/i5SR5btt/haexE48epqr1VdbCqDh45cmSd0wMAAGBRbPTFmWqFsV5l/PjB7tu6e3d37965c+eGTg4AAICtZ73h+pnpFOBMn5+cxg8nOX/ZfucleXyVcQAAAFjVesN1X5KjVwbek+T9y8bfMF1d+JVJnp5OJf5Aksur6qzpokyXT2MAAACwqh1r7VBV707yqiRnV9XhLF0d+G1J7q6qa5M8muTqafd7klyZ5FCSzyV5Y5J091NV9ZNJ7pv2e2t3H3vBJwAAADjOmuHa3a8/waZLV9i3k1x3gu9zR5I7Tml2AAAAbHsbfXEmAAAA2FDCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIa2Y94TAGBBHbj5i+9fcuN85gEAbHnCFQDYFm7Z//BxY9dfdtEcZgLAqXKqMAAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADG3HvCcAADAvt+x/+Lix6y+7aA4zAWA1VlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGgzhWtVXV9VD1bVx6vq3VX1pVV1QVV9uKoeqar3VNVzp32fN90/NG3ftRE/AAAAAItt3e/jWlXnJvlnSV7W3f+3qu5Ock2SK5Pc0t13VdUvJbk2ya3T589299dV1TVJfjrJ98z8EwAAHGOl92cFYOua9VThHUn+WlXtSPJlSZ5I8u1J3jttvzPJ66bbV033M22/tKpqxscHAABgwa07XLv7D5P8bJJHsxSsTye5P8mfdPcz026Hk5w73T43yWPT1z4z7f/CY79vVe2tqoNVdfDIkSPrnR4AAAALYt3hWlVnZWkV9YIkX5Pky5O8ZoVd++iXrLLt2YHu27p7d3fv3rlz53qnBwAAwIKY5VTh70jy+919pLv/X5L3Jfk7Sc6cTh1OkvOSPD7dPpzk/CSZtn9VkqdmeHwAAAC2gVnC9dEkr6yqL5teq3ppkk8kOZDku6Z99iR5/3R733Q/0/bf7O7jVlwBAABguVle4/rhLF1k6YEkH5u+121J3pLkzVV1KEuvYb19+pLbk7xwGn9zkhtmmDcAAADbxLrfDidJuvumJDcdM/zpJBevsO+fJ7l6lscDAABg+5n17XAAAADgtBKuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQdsx7AgDMyYGb5z0DAICTYsUVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGJpwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAICh7Zj3BAAAZnXL/ofnPQUATiMrrgAAAAxNuAIAADC0mcK1qs6sqvdW1Ser6qGq+paqekFV7a+qR6bPZ037VlW9o6oOVdVHq+rlG/MjAAAAsMhmXXF9e5L/0t1/M8k3JHkoyQ1JPtjdFyb54HQ/SV6T5MLpY2+SW2d8bAAAALaBdYdrVX1lkm9LcnuSdPfnu/tPklyV5M5ptzuTvG66fVWSd/aSe5OcWVXnrHvmAAAAbAuzrLj+jSRHkvxKVf1uVf1yVX15khd39xNJMn1+0bT/uUkeW/b1h6cxAAAAOKFZwnVHkpcnubW7vynJ/8mzpwWvpFYY6+N2qtpbVQer6uCRI0dmmB4AAACLYJb3cT2c5HB3f3i6/94shetnquqc7n5iOhX4yWX7n7/s689L8vix37S7b0tyW5Ls3r37uLAFADidjn1P2Osvu2hOMwHgqHWvuHb3HyV5rKpeOg1dmuQTSfYl2TON7Uny/un2viRvmK4u/MokTx89pRgAAABOZJYV1yT5oSTvqqrnJvl0kjdmKYbvrqprkzya5Opp33uSXJnkUJLPTfsCAADAqmYK1+7+SJLdK2y6dIV9O8l1szweAAAA28+s7+MKAAAAp5VwBQAAYGjCFQAAgKEJVwAAAIYmXAEAABiacAUAAGBowhUAAIChCVcAAACGJlwBAAAYmnAFAABgaMIVAACAoQlXAAAAhiZcAQAAGNqOeU8AgG3iwM3Hj11y4+bPAwDYcqy4AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEPbMe8JALAJDtw87xkAAKybFVcAAACGJlwBAAAYmnAFAABgaMIVAACAoc0crlV1RlX9blX9p+n+BVX14ap6pKreU1XPncafN90/NG3fNetjAwAAsPg24qrCP5zkoSRfOd3/6SS3dPddVfVLSa5Ncuv0+bPd/XVVdc203/dswOMDAJw2t+x/+Lix6y+7aA4zAdi+Zlpxrarzkvz9JL883a8k357kvdMudyZ53XT7qul+pu2XTvsDAADACc264vrzSX4syfOn+y9M8ifd/cx0/3CSc6fb5yZ5LEm6+5mqenra/4+Xf8Oq2ptkb5K85CUvmXF6AMCiWWkFFIDFtu4V16r6ziRPdvf9y4dX2LVPYtuzA923dffu7t69c+fO9U4PAACABTHLiuu3JnltVV2Z5Euz9BrXn09yZlXtmFZdz0vy+LT/4STnJzlcVTuSfFWSp2Z4fAAAALaBda+4dveN3X1ed+9Kck2S3+zu701yIMl3TbvtSfL+6fa+6X6m7b/Z3cetuAIAAMByp+N9XN+S5M1VdShLr2G9fRq/PckLp/E3J7nhNDw2AAAAC2Yj3g4n3f1bSX5ruv3pJBevsM+fJ7l6Ix4PAACA7eN0rLgCAADAhhGuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0IQrAAAAQxOuAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQhCsAAABD2zHvCQCwjR24+fixS27c/HkAAEOz4goAAMDQhCsAAABDE64AAAAMTbgCAAAwNOEKAADA0FxVGGARrXS1XgCALcqKKwAAAEOz4goM6Rc/8otfdP8Hv/EH5zQTAADmTbgC6yIsAQDYLE4VBgAAYGhWXIG5sGILAMDJEq7Aik41LE93iApdYCS37H/4uLHrL7toDjMB2B6EKzCEY8MUAACO8hpXAAAAhmbFFbYgp80C28VKp+QCsP0IVyCJU3UBABiXcIUFdDIrskIVAICtQrgCp8VopzOPNh8AAE6ecAUWwqxv33MyXwMAwHwIVxjAdlgN3OhTk53qDACwfXg7HAAAAIZmxRVYSFZkAQAWh3AFTsp2CMHtcMo2AMBWJFxhAawVldshOmFR+YUKAMwQrlV1fpJ3JvnqJH+V5LbufntVvSDJe5LsSvIHSb67uz9bVZXk7UmuTPK5JN/f3Q/MNn0A2N6ELQDbwSwrrs8k+efd/UBVPT/J/VW1P8n3J/lgd7+tqm5IckOStyR5TZILp49XJLl1+gyswYopcNSpnmEhZAFYBOsO1+5+IskT0+0/q6qHkpyb5Kokr5p2uzPJb2UpXK9K8s7u7iT3VtWZVXXO9H0AYFsSmgCwtg15jWtV7UryTUk+nOTFR2O0u5+oqhdNu52b5LFlX3Z4GvuicK2qvUn2JslLXvKSjZgebLhZ/6HpNanAifjzDwDHmzlcq+orkvxakh/p7j9deinryruuMNbHDXTfluS2JNm9e/dx22Er8g9R4CjPBwBw6p4zyxdX1ZdkKVrf1d3vm4Y/U1XnTNvPSfLkNH44yfnLvvy8JI/P8vgAAAAsvnWH63SV4NuTPNTdP7ds074ke6bbe5K8f9n4G2rJK5M87fWtAAAArGWWU4W/Ncn3JflYVX1kGvvxJG9LcndVXZvk0SRXT9vuydJb4RzK0tvhvHGGx4ahORWQTXXg5nnPAADgtJrlqsL/PSu/bjVJLl1h/05y3XofDwCAMbgaNrDZNuSqwgCwYVZaQb7kxs2fBwAwDOEKG8CpwcCJeH5gO7ACC5xuwhUAgE23VuyKYWA54QoAwMITwrC1CVcAADbUeiLRafXAatb9Pq4AAACwGYQrAAAAQ3OqMAAAq5r1NF6nAQOzEq4AAAzPxZVgexOuALCBrCzB5pg1ZIUwbC3CFQCAheOXSLBYhCsAMIRb9j887ykAMCjhCifBb20BYGz+robFJlwBAGANXhML8yVcAQA2wLGnOl9/2UVzmsnsrF4CoxGuAABwjFONdyuycHoJVwDGd+DmL75/yY3zmQcAMBfCFWArOTbgANgQs66wAqeXcAUAgA3m1GHYWM+Z9wQAAABgNVZcYQVO/wEA5s2qLTxLuEKEKgAAjMypwgAAAAzNiisAwGlwy/6Hjxu7/rKL5jATRuDsLpiNcGVb8pcHwPytFHYAsBLhysJZKUpdzAAAALYu4QoAAFvAWlcZdhViFplwBQCALchLn9hOXFUYAACAoVlxZVvwG0m2rAM3z3sGY1rpv8slN27+PAC2EKcSs5VZcQUAAGBoVlzZ8qymAgDAYhOuAMBp5z1bAZiFU4UBAAAYmhVXAIBNstLK8/WXXTSHmQBsLcIVAAAW0FrXAXGVYbYS4cpwXGyJbcHb3AAwmLX+DSZsmSfhCsBi8N6uALCwhCtzZ4UVOG2OjVkhu+FcLRi2L6cas5mEKwAAsCaLDcyTcGXTedJj2/F6VgC2Aa+R5XQSrgBsH14Hy4C8RQ7A2oQrwEazwsoC85rWzXHsf2chC2x3wpXTzqnBLAyrdYvJ/1cAGJ5wBZiF1VUWmNVV4HRaz+KG18luX8KVDWeFlYUhSrcvb6PDYE72lwhOKWZk/o3ILDY9XKvqiiRvT3JGkl/u7rdt9hyYjScdYNvZ4qcTWznd3rxeFlgEmxquVXVGkl9IclmSw0nuq6p93f2JzZwHAMAiOplfUpzMPg/86f/Kt3ztCzdiSrChNnoBxanHW8dmr7henORQd386SarqriRXJRGu67TSH961/gBaMQXYAFt8FRaA453qv5NnDd9jH09In1h19+Y9WNV3Jbmiu//RdP/7kryiu9+0bJ+9SfZOd1+a5FObNsHNd3aSP573JGADOJZZBI5jFoVjmUXhWN4e/np371xrp81eca0Vxr6onLv7tiS3bc505quqDnb37nnPA2blWGYROI5ZFI5lFoVjmeWes8mPdzjJ+cvun5fk8U2eAwAAAFvIZofrfUkurKoLquq5Sa5Jsm+T5wAAAMAWsqmnCnf3M1X1piQfyNLb4dzR3Q9u5hwGsy1OiWZbcCyzCBzHLArHMovCscwXbOrFmQAAAOBUbfapwgAAAHBKhCsAAABDE64brKruqKonq+rjy8beU1UfmT7+oKo+smzbjVV1qKo+VVWvXjZ+xTR2qKpu2Oyfg+3tBMfxN1bVvdNxfLCqLp7Gq6reMR2rH62qly/7mj1V9cj0sWcePwvb2wmO5W+oqt+uqo9V1X+sqq9cts1zMsOpqvOr6kBVPVRVD1bVD0/jL6iq/dNz7P6qOmsa97zMkFY5lq+e7v9VVe0+5ms8L5PEa1w3XFV9W5L/neSd3f31K2z/10me7u63VtXLkrw7ycVJvibJf01y0bTrw0kuy9JbCN2X5PXd/YlN+BFgxeO4qn4jyS3d/etVdWWSH+vuV023fyjJlUlekeTt3f2KqnpBkoNJdmfp/ZrvT/LN3f3ZOfxIbFMnOJbvS/Kj3f3fquoHklzQ3f/SczKjqqpzkpzT3Q9U1fOz9Hz6uiTfn+Sp7n7b9A/3s7r7LZ6XGdUqx3In+ask/zZLz88Hp/09L/MFVlw3WHd/KMlTK22rqkry3Vn6A5gkVyW5q7v/ort/P8mhLP3BvDjJoe7+dHd/Psld076wKU5wHHeSoytTX5Vn34P5qixFQXf3vUnOnP5ienWS/d391PSPov1Jrjj9s4dnneBYfmmSD0239yf5h9Ntz8kMqbuf6O4Hptt/luShJOdm6Ti8c9rtziwFQOJ5mUGd6Fju7oe6+1MrfInnZb5AuG6uv5vkM939yHT/3CSPLdt+eBo70TjM048k+ZmqeizJzya5cRp3HLPVfDzJa6fbVyc5f7rtWGZ4VbUryTcl+XCSF3f3E8lSECR50bSbY5nhHXMsn4hjmS8Qrpvr9Xl2tTVJaoV9epVxmKd/muT67j4/yfVJbp/GHcdsNT+Q5Lqquj/J85N8fhp3LDO0qvqKJL+W5Ee6+09X23WFMccyw3Assx7CdZNU1Y4k/yDJe5YNH86zv+lPkvOydPrlicZhnvYked90+99n6TSdxHHMFtPdn+zuy7v7m7P0y8T/MW1yLDOsqvqSLP1D/13dffS5+DPTKcBHXzv45DTuWGZYJziWT8SxzBcI183zHUk+2d2Hl43tS3JNVT2vqi5IcmGS38nSC8wvrKoLquq5Sa6Z9oV5ejzJ35tuf3uSo6e870vyhukqlq/M0sXHnkjygSSXV9VZ05UuL5/GYK6q6kXT5+ck+RdJfmna5DmZIU3XyLg9yUPd/XPLNu3L0i8VM31+/7Jxz8sMZ5Vj+UQ8L/MFO+Y9gUVTVe9O8qokZ1fV4SQ3dfftWfoDtfw04XT3g1V1d5JPJHkmyXXd/ZfT93lTlv4yOSPJHd394Ob9FGx3Kx3HSf5xkrdPZw/8eZK90+73ZOnKlYeSfC7JG5Oku5+qqp/M0l8uSfLW7l7xwmVwupzgWP6Kqrpu2uV9SX4l8ZzM0L41yfcl+Vg9+5Z6P57kbUnurqprkzyapddsJ56XGdeJjuXnJfk3SXYm+c9V9ZHufrXnZZbzdjgAAAAMzanCAAAADE24AgAAMDThCgAAwNCEKwAAAEMTrgAAAAxNuAIAADA04QoAAMDQ/j//WVdQ4r67kAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(MC_Ds_sig[:,dim-1],alpha=0.5,bins=30);\n",
    "plt.hist(MC_Dplus_sig[:,dim-1],alpha=0.5,bins=30);\n",
    "plt.hist(data_bkg[:,dim-1],alpha=0.5,bins=200);\n",
    "fig=plt.gcf()\n",
    "fig.set_size_inches(16,8)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Cut on Ds Mass\n",
    "\n",
    "#data_bkg_cut=data_bkg[np.where(data_bkg[:,dim-1]>1915)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Add 0/1 label for bkg/sig\n",
    "\n",
    "MC_Dplus_sig_labelled=add_labels(MC_Dplus_sig,signal=True)\n",
    "MC_Ds_sig_labelled=add_labels(MC_Ds_sig,signal=True)\n",
    "data_bkg_labelled=add_labels(data_bkg,signal=False)\n",
    "\n",
    "#Merge MC sig and data bkg, shuffle it\n",
    "\n",
    "data=np.concatenate((data_bkg_labelled,MC_Ds_sig_labelled), axis =0)\n",
    "data=np.concatenate((data,MC_Dplus_sig_labelled), axis =0)\n",
    "np.random.seed(1)\n",
    "np.random.shuffle(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get train size\n",
    "train_size=data.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Strip away the label column and convert it to a one-hot encoding\n",
    "\n",
    "X=data[:,0:dim]\n",
    "Y_labels=data[:,dim].astype(int)\n",
    "Y_labels=Y_labels.reshape(train_size,1)\n",
    "#Y_labels_hot = to_one_hot(Y_labels)\n",
    "#Y_labels=Y_labels_hot\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Divide the dataset k \"equi populated\" sets\n",
    "\n",
    "k=5 #number of subsets\n",
    "i=0 #number of subset that is train set\n",
    "\n",
    "X_train, Y_train, X_test, Y_test, X_dict, Y_dict = k_subsets(i, k, X, Y_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.528952660100201"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.float((m_s+m_plus)/(train_size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#np.argmax(Y_train, axis=1).sum()/Y_train.shape[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_mean=X_train.mean(axis=0)\n",
    "X_train_1=X_train-X_mean\n",
    "X_std=X_train_1.std(axis=0)\n",
    "X_train_2=X_train_1/X_std\n",
    "X_train_2.std(axis=0)\n",
    "\n",
    "\n",
    "X_mean=X_test.mean(axis=0)\n",
    "X_test_1=X_test-X_mean\n",
    "X_std=X_test_1.std(axis=0)\n",
    "X_test_2=X_test_1/X_std\n",
    "X_test_2.std(axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# SETTING UP THE NETWORK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,\n",
       "       colsample_bytree=1, gamma=0, learning_rate=0.1, max_delta_step=0,\n",
       "       max_depth=3, min_child_weight=1, missing=None, n_estimators=100,\n",
       "       n_jobs=1, nthread=None, objective='binary:logistic', random_state=0,\n",
       "       reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,\n",
       "       silent=True, subsample=1)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# fit model to training data\n",
    "model = XGBClassifier()\n",
    "\n",
    "y=Y_train.reshape(Y_train.shape[0])\n",
    "\n",
    "model.fit(X_train_2, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "output=model.predict_proba(X_test_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA8EAAAIgCAYAAAClJ2uTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmcjfX///HnyxiGyWCYhUGUtaRisqRF9owlpMS3bKVF9S31Gfqo0DcU+ejTHmWpRFKWIltJm0K/TwuhFLKv0WIf798f55rzmeXMGMxC1+N+u53bnPO+3td1vc9p8p7neb+v92XOOQEAAAAA4AeFCroBAAAAAADkF0IwAAAAAMA3CMEAAAAAAN8gBAMAAAAAfIMQDAAAAADwDUIwAAAAAMA3CMEAgDxlZj3N7Gsz+8PMfjOz/5jZv7Ko+4OZ3ek9d2kex81sq5m9ZWZV8qndG8zsqVw+ZqyZDTGzyhnKm3jvs3Zunu8EbelgZqvN7IiZbciv8+ZUmv/2jTKU1/bKm6Qp+9greyjEcXab2ZC8bzEA4GxBCAYA5BkvlLwiab6kTpJukTRLUvsQdatIqiVpbpri0ZIaSWos6UFJdSXNMbPCedvyPBMrabCkyhnK/58C7/Pn/GiEmYVJek3St5KaSuqYH+c9RQ+fRN37zax4nrUEAPC3cLb+EQEAODvcLell59w/05S9Z2ZDQ9RNkrTSObcxTdkG59yX3vOlZrZP0hxJ1SX9kCctLgDOud8lfXnCirmnnKQoSW865z7LqpKZFXPOHcy/ZmXysaQ2Znapc+4/J6i7VIEvSfpKejqvGwYAOHsxEgwAyEulJG3PWOiccyHqJikQcLPzh/czPLtKZtbHzFaZ2UFvOuwSM7swzfYIMxtpZpvM7LCZfWtmbU5wbpnZFd6xDpjZHjMbZ2YlMtQ518ymeOc9YGbfmVk3bwr09161xanTfb19Mk2HNrPiZvaMmW03s0NmttzMWmY418dmNt07/joz+93MPjCzCtm8h56SNnkvZ3nnHeJtc2bW38yeNrNdadorM7vbzH7yPq91ZnZ/huMO8d5zAzNb4X32n5lZFW8a+Ewz+9Obgt30RJ+1510FvuwYlIO6WyVNkPSgmRXN4fEBAD5ECAYA5KX/J+keM+thZmWyquRNYW2izCG4kJkVNrNwM6suaaiknyStzOZYV0l6SdIbkq6V1FvSF5JKpqk2XVJPScMltZO0XNJsM7skm+M2lvShAqH+ekn3SWqjQPBKrROrwIjkZQpM324n6VVJFSVtk9Tdq9pPgenP6a53zWCcpF6ShikwXXmTAlPBr8hQr4ECI+4PKDAKWlfS2GyOO0eBqeny2thIgSnrqf6hwEjxzZLu9d7XbZKelTTbe09vSxptZgMzHLu4d+4xkm6SVEnS65KmSPrMO+8WSW/ncNqyU+C/USczuyAH9Z+UFKfA5wYAQEhMhwYA5KV+kmZKmijJmdlqSe9IesqbApyqmaSDCoTVtP7tPVJtltTGOZeSzTnrS/rOOTciTdns1Cdm1kyBUecmzrklXvECL2QPktQli+M+IekL59yNaY61RdKHZlbbObdS0v0KhO16zrltXrUP09T/znv6Q5pp3pmYWS0FQmQv59wkr2y+pO8kPSKpVZrqUZKSnHO/efXiJY3Jaiqzc26XmaVOLV4boh3bM7zHQpKGSJronHvAK15gZiUlPWRmTzvnDnnlxSTdm/q5mll5Sc9LGuyce8or2yxplaSrJX2Q1WeQxlQFvvx4SIFgniXn3AYzmyxpgJm94pw7loPjAwB8hpFgAECecc59p8BiV+0lvSDJFAhxK8zsnDRVkyTNDxFuRykwqnqZV+c7SXPNLCGb034j6VIzG2NmV5lZkQzbmyswmvu5N8pc2Fto60NJiaEO6I1aNpI0LcM+n0k6KqmeV7WppHlpAvCpukyBz+rt1ALn3HHvdcaR4OWpAdiTeq10dp9RdjKOxleQVD5tWzxvKRDAL0pTdkTSp2ler/N+fhSiLEft834nnpB0k5mdn4NdhiswAt39RBUBAP5ECAYA5Cnn3GHn3HvOubudcxdIulVSNUl90lRro9DXA//qnFvhPeYqMJ02QoER16zOt0iB6bBXKbCw0m4ze8HMIr0qZSXFKxBe0z6GKDBtOZTSksIUCPJp9zmswPXJqfuVUWDa8+kqJ+lP59yBDOU7JBXPcM3rvgx1jng/I07x3DtCtCVUeerr6DRlf3hhPWNbgm10zp1K+15T4JrfASeq6Jz7UYHp7g95o9gAAKTDdGgAQL5yzr1qZiMl1ZQkM6ujwKjgvBzse9jMflFgdDm7epMkTTKzGAWC8xhJv0saKGmvAtelXncSzd6nwPWpQ5T+Fk6ptno/9+i/ofF0bJN0jpkVzxCE4yQdcM4dzoVzZCXjomWpoT42Q3mc93NvHrZFUiA4m9koSU8psFjWiQxTYEbA9XnaMADAWYlvSAEAecZbKCpjWYwC182mjiQmSfrKObc7B8eLkHS+/ru6cbacc7uccy8rMEU3dWGlDxUYCf4zzShz8JHFcf5S4BZGNULt45zbmubYrcwsLtRxlPNR2uUKhNFgiDMz815neUujPLJZgZCf8VrpGxT4YuH7THvkjXGSfpOUfKKK3jT89yT9U4Fp5QAABDESDADIS9+b2SxJCyTtlHSuAisSH5A0yauT3a2RKptZQ+95jKS7FAjQr2Z1Qgvcgzha3lRoSZcqsAhT6krGCyXNl7TQzJ5UYJGmKEmXSIpwzj2UxaGTFVgE67gC023/UODa0yRJg7xpuGMk3SLpUzMbpkBYryUp0jk3UtKvCiwA1sPM9ks6Gip4O+dWm9kUSc+ZWZQC19HepsDo+Z1Zvfe84Jw77t1C6WUz26PA53e1145/plkUK6/bccjM/qXACtA5MUzSV3nYJADAWYoQDADIS49J6iDpGQWC6XYFVoC+0Tm33syiJTVU4BY/oTzgPaTAVOPvJbV0zi3P5pzLFbhmuKukEpI2KjCN+d9S4B7FZtZJgVHC+xQIsnsVmD77bFYHdc595t1+aagCt/0J8449T96otrfycmNJIyU9LamoArd0GuFtP+TdbmiwpCUKXE+c1UjlbQoEvkcUuN/y95LaOufyeyRYzrlx3nXI90n6XwVGhx9wzo3J56a8oMB1wdEnquicW2ZmCyW1yPNWAQDOKuZcxkt/AADIH2bWTdJI51yFgm4LAADwB0IwAAAAAMA3WBgLAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4Rr6GYDMbb2Y7zWxlhvJ7zGytma0ys5Fpyh8ys3XetlZpylt7ZevMbGB+vgcAAAAAwNnLnHP5dzKzqyT9Kek151xtr+waSYMkJTnnDptZrHNup5ldIGmKpPqSyktaJKm6d6gfJbWQtFnSckk3Oed+yLc3AgAAAAA4KxXOz5M55z4xs8oZiu+U9IRz7rBXZ6dX3kHSVK98vZmtUyAQS9I659wvkmRmU726hGAAAAAAQLbOhGuCq0u60sy+MrMlZnaZV54gaVOaepu9sqzKMzGzvma2wnv0zYO2AwAAAADOIvk6EpyFwpJKS2oo6TJJ08zsPEkWoq5T6OAeck63c26spLGSVLZsWZeYmPhyrrQYAOB7X3/99W7nXExBt+NsVrZsWVe5cuWCbgYA4G8ip33zmRCCN0t61wUuTl5mZscllfXKK6apV0HSVu95VuVZqly5slasWJE7LQYA+J6ZbSzoNpzt6JsBALkpp33zmTAdeqakppJkZtUlFZG0W9JsSV3NrKiZVZFUTdIyBRbCqmZmVcysiKSuXl0AAAAAALKVryPBZjZFUhNJZc1ss6TBksZLGu/dNumIpB7eqPAqM5umwIJXxyT1c86leMe5W9J8SWGSxjvnVuXn+wAAAAAAnJ3ye3Xom7LY9D9Z1B8maViI8rmS5uZi0wAAAAAAPnAmTIcGAAAAACBfEIIBAAAAAL5xJqwODeBv4Pfff9fOnTt19OjRgm4KcNrCw8MVGxurqKiogm4KgHxGfwacmXKzbyYEAzhtv//+u3bs2KGEhAQVK1ZMZqFu8w2cHZxzOnjwoLZs2SJJBGHAR+jPgDNTbvfNTIcGcNp27typhIQEFS9enD8YcNYzMxUvXlwJCQnauXNnQTcHQD6iPwPOTLndNxOCAZy2o0ePqlixYgXdDCBXFStWjOmQgM/QnwFnttzqmwnBAHIF35jj74bfacCf+H8fOHPl1v+fhGAAAAAAgG8QggFA0pAhQ2RmqlatWsjtVatWlZlpyJAhmba9++67atq0qUqVKqWiRYuqevXqevjhh7V79+7Tak/ZsmVPef/c0KRJE11//fUF2oZTMXbsWM2cObOgmwEABWbixImqV6+eSpQoodKlS+vSSy9V//798/y8H3/8scxMK1euzNXjLliwQE8//XSuHrNs2bIh+/S0tm3bpjZt2qhkyZIyM3388ce52obTlZf9XX72pUeOHNGQIUP0zTff5Mv5JEIwAARFRERo/fr1WrFiRbry5cuXa+PGjYqIiMi0zwMPPKAuXbrovPPO0+uvv64FCxbo/vvv13vvvafbbrvtlNty6623av78+ae8v58RggH42YgRI3TrrbeqVatWevfdd/Xaa6+pQ4cOmj17dkE37ZTlRQjOiWHDhunbb7/VlClTtHTpUtWtWzff25Cdv1MIHjp0aL6GYG6RBCDPVB44p0DOu+GJpFPaLzIyUnXr1tXUqVOVmJgYLJ86daqaNm2qr7/+Ol399957T//617/06quvqnfv3sHyq6++Wn379tWCBQtO7Q1IqlChgipUqHDK+wMAcs9D735fIOcd0emik97nueee0+23367hw4cHy9q1a6fBgwfnZtN8Yc2aNWrQoIHatGmTZZ2UlBSlpKSoSJEi+dgynC5GggEgja5du2ratGlyzkkK3Jdu2rRp6tq1a6a6Y8aMUd26ddMF4FRhYWG69tprszzPvn37dOutt6p8+fKKiIhQpUqV0o0ch5oO/d133+nyyy9XRESELrzwQs2dO1eJiYnq2bNnsE7Pnj2VmJiohQsXqk6dOoqMjNQVV1yhVatWpTvW6NGjddlll6lkyZKKi4tTu3bttG7duhx9RhnNmjVLiYmJioiIUHx8vJKTk4MrNy5evFhmlun8v/32m4oUKaJXX301WPbZZ5/p6quvVvHixVWmTBnddttt+uOPP4LbJ06cKDPT999/rxYtWigyMlI1a9bUu+++G6zTpEkTff3115o0aZLMTGamiRMnSpJmz56tevXqKTIyUqVLl1aDBg20ZMmSU3rPfmFmpcxsupmtMbPVZtbIzKLNbKGZ/eT9LO3VNTN7xszWmdl3ZnZmDZkAPrFv3z7Fx8dnKj/RgkIn6pckaeXKlUpKSlKJEiVUokQJdenSRdu3b8/2uMePH9cTTzyhqlWrBi8ZmjRpUqZ6M2bMUP369VWsWDGVKVNGbdq00caNGzVkyBCNHj1aGzduDP67nrbfO1HfIUmffPKJLr74YkVERKhevXr64osvsm2zFPi8PvzwQ82YMUNmpsqVK0v6bz87c+ZMXXjhhYqIiNBXX30lSfrmm2/UrFkzFS9eXKVLl1b37t21Y8eO4DE3bNggM9PUqVPVq1cvRUVFqUKFCnrjjTckSSNHjlT58uUVExOjAQMG6Pjx41m2L7v+TpJeeeUVXXjhhSpatKjOPfdcjRw5Mt3+q1atUuvWrRUdHa3IyEjVqlVLzz//fI6OndGIESNUtWpVRUREKC4uTq1bt073e7F3717dfvvtiouLU0REhC6//PLgZyZJJUqUkCT16tUreL4NGzZk/R8nFxCCASCNTp06aceOHfrss88kSZ9++ql27dqljh07pqt39OhRffHFF2rduvUpnad///767LPPNGbMGM2fP1/Dhw/P9g+UAwcOqFWrVjp48KCmTJmihx9+WPfff79+/fXXTHV//fVX/eMf/9CgQYM0ZcoU7dy5UzfccEMw2EvS5s2bdffdd2vWrFkaN26cUlJS1LhxY+3fv/+k3se0adPUqVMn1a9fX7Nnz9bgwYM1duxYPfTQQ5ICo+LlypXTtGnT0u03Y8YMSQp+rp9//rmaNWum+Ph4TZ8+XU8//bTmzp2rXr16ZTpnt27d1L59e82YMUPVqlVT165dtXnzZknSCy+8oJo1a6pNmzZaunSpli5dqqSkJP3888+6/vrr1bRpU7333nuaPHmy2rZtq717957U+/Whf0ua55yrKeliSaslDZT0oXOumqQPvdeSdK2kat6jr6QX87+5AOrWratnn31WkyZN0p49e3K834n6pXXr1qlx48Y6dOiQXn/9dU2cOFGrVq1Su3bt0vUvGd1zzz16/PHH1bdvX82ZM0cdO3ZU79699f777wfrvP766+rUqZPOP/98TZs2TRMmTFD16tW1a9cu3XrrrerWrZvi4+OD/64/8sgjknLWd2zdulXXXnutoqOjNX36dN1+++3q3r27Dhw4kO3nsXTpUl166aW65pprtHTp0mC/JQXCbHJysh566CHNnTtXVapU0a5du9SkSRMdOHBAb775pp599lktWbJELVq00JEjR9Ide8CAASpXrpzeeecdXXnllerRo4ceeOABLVu2TOPHj9d9992nkSNHZuo708qqv5OkUaNG6c4779R1112n999/X3feeaceeeQRPffcc8H927dvr7CwML3xxhuaPXu27rnnnuCXB9kdO6PXXntNw4cPV//+/TV//ny9+OKLqlq1qv766y9J0uHDh9W8eXMtXLhQo0aN0syZMxUTE6PmzZsHg/JHH30kSXr44YeD5ytXrly2/31OF9OhASCNUqVKqXXr1po6daquvPJKTZ06Va1bt1apUqXS1duzZ48OHz6sSpUqndJ5li1bpn79+unGG28Mlv3P//xPlvUnTJigPXv2aMWKFUpISJAknX/++WrQoEGmunv37tXnn38eXOTr+PHj6tixo9auXauaNWtKCoxip0pJSVGLFi0UGxurWbNm6ZZbbsnRe3DO6R//+IduueUWvfDCC8HyokWLql+/fnrooYdUpkwZdenSRW+99ZaGDh0arPPWW2+pZcuWio6OliQNHDhQl19+ud56661gnYSEBDVr1kwrV65U7dq1g+X3339/cPS9Xr16iouL0/vvv6877rhDF1xwgSIjIxUTE6OGDRsG91myZIlKlCihUaNGBcuym94GycyiJF0lqackOeeOSDpiZh0kNfGqTZL0saQBkjpIes0F/hr+0htFLuec25bPTQd87fnnn9d1112nnj17ysxUq1Ytde7cWQ8++KCioqKy3O9E/dLQoUMVHx+vDz74IDj1t06dOqpZs6bmzp0bMiStW7dOL774oiZMmKAePXpIkpo3b65t27Zp6NChatu2rY4fP66BAweqY8eOmjJlSnDf9u3bB5+XK1dORYsWTffvupSzvuPpp59WRESE5syZo+LFi0sKXP6UXZ8rSQ0bNlRUVJSio6MznXfPnj1atGiRLrnkknRtkaT58+cHP+fq1aurQYMGeuedd3TTTTcF6zZt2jQ4Xb1BgwaaPn26Zs+erTVr1igsLEytW7fWrFmzNGPGjJAz0SRl2d/9/vvvGjp0qB5++OHgFPgWLVrowIEDevzxx3XnnXfqt99+0y+//KKZM2fqoosCU+6bNWt2wmOHsmzZMrVs2VJ33XVXsKxTp07B52+88YZWrlypVatWBf8uad68uWrUqKHRo0dr1KhRuuyyyyQF/q450flyCyPBAJBB165dNX36dB0+fFjTp0/PsgOSTv1+dZdccolGjRqlF154QT/++OMJ6y9fvlz16tULBmBJql+/vuLi4jLVrVy5crpVri+44AJJCo6WStKXX36pFi1aqEyZMipcuLCKFy+uP//8M0dtSfXjjz/q119/1Q033KBjx44FH02bNtWhQ4eCq4PeeOONWrt2rb799ltJ0u7du/XRRx8F/9A6cOCAli5dmuk4V1xxhcLDwzNdi92yZcvg8zJlyig2Njbdewvloosu0v79+9WjRw8tWLAg+A01snWepF2SJpjZf8zsFTOLlBSXGmy9n7Fe/QRJm9Lsv9krA5CP6tSpo9WrV2v27Nm666675JzT//3f/ykxMVF//vlnlvudqF9atGiROnbsqEKFCgX/na5SpYoqV66caUHJVB9++KEKFSqkjh07pvv3vVmzZvrmm2+UkpKitWvXauvWrSFn/mQnp33HsmXL1KJFi2AAltKHtFORkJCQLgCnnqdly5bpvmioX7++KleuHJxdlipt4IyKilJMTIyuvvpqhYWFBcurVq2qLVu2nHTbli5dqr/++ktdunTJ1Dfv2LFDmzdvVnR0tCpWrKg77rhDb731lnbu3HnS50l1ySWXaO7cuRo8eLCWLVumlJSUdNsXLVqkevXqqUqVKsG2SIGZYln93uQHQjDyVeWBczI9gDNN+/bt9eeff2rQoEH666+/1K5du0x1ypQpo6JFi4acjpwTzz33nK677jo99thjqlGjhqpVq6apU6dmWX/79u2KiYnJVB6qLOOodeo39ocOHZIUmC7dsmVLOef08ssv6/PPP9fy5csVGxsbrJMTqbeAatOmjcLDw4OPKlWqSJI2bQrkoUaNGqlSpUrBb+rfeecdFS5cWNddd52kwPXBKSkpuuuuu9Idp2jRojp69GjwONm9vxO1u0aNGpo1a5Z++eUXtWnTRmXLllW3bt20a9euHL9fHyosqa6kF51zl0r6S/+d+hxKqG+EMs2RNLO+ZrbCzFbw+eOscvxY5scZqmjRomrXrp2ee+45/fDDD3rllVf0008/pVuHIaMT9Uu7d+/Wk08+me7f6fDwcP3yyy+Z/p1Ou09KSopKliyZbp+ePXvq2LFj2rZtW3DK9slOf81p37F9+3bFxsam27dYsWI655xzTup8aYX6Anrbtm0hy+Pi4jJdehOqHzuVvi2U1L75wgsvTPe5XHPNNZICfXOhQoW0YMECxcfHq3fv3oqPj9eVV16p//znPyd9vt69e2v48OGaNm2aGjRooLi4OD3yyCPBMLx79259+eWXmX5vJkyYkOXvTX5gOjQAZBAZGam2bdtqzJgx6tKliyIjIzPVCQ8PV+PGjTV//nw9/vjjJ32OUqVK6ZlnntEzzzyj7777TiNHjlT37t1Vp06d4MhtWvHx8Vq7dm2m8lMJEfPmzdOBAwc0a9as4Hs7duzYSV8fmzqVeezYsbr00kszbU8Nw2amG264QW+99ZaGDx+ut956S9dee21wIYxSpUoF78Ecaopy+fLlT6pdWUlKSlJSUpL279+vOXPm6L777tM999yT7ZcPPrdZ0mbnXOrqJdMVCME7Uqc5m1k5STvT1K+YZv8KkrZmPKhzbqyksZKUmJiY9YWEAHJNnz59lJycrDVr1mRZ50T9UnR0tDp27Khbb701075Z3dc+OjpahQsX1ueff65ChTKPvcXGxgavQ9227eSunMhp3xEfH59ppPPgwYPZjoqfSKhZYOXKlQs5orpjxw7Vq1fvlM91slL75vfffz9kKK9Ro4YkqWbNmnrnnXd09OhRffrppxowYICSkpK0efPmkP+tslKoUCHdf//9uv/++7Vp0yZNnjxZgwYNUkJCgu644w5FR0crMTFRL76YeZmIokWLnuK7PH2EYAAI4c4779Thw4d1xx13ZFnnvvvuU/v27TVp0qTgtU6pjh8/rgULFuRo4aw6depo1KhRmjx5stasWRMyBF922WWaPHmytmzZEpwSvWzZsnSrTubUwYMHVahQIRUu/N8uYNq0acEpSjlVo0YNJSQkaMOGDSe8J3LXrl311FNP6f3339eSJUvSXfcVGRmphg0bau3atXr00UdP7s2EcKJvz0uWLKlu3bppyZIlWrp06Wmf7+/KObfdzDaZWQ3n3FpJzST94D16SHrC+znL22W2pLvNbKqkBpL2cz0wkP927tyZaeRz165d2r9/f8hQFEqofin1Ott69erl+FKgpk2bKiUlRfv371eLFi1C1kntSyZNmhRy5pUU+t/1nPYdl112mcaPH68DBw4Ep0SnvatAbmnQoIFefPFF/fHHH8EveZcvX64NGzboiiuuyPXzSaE/l0aNGqlYsWLaunVrlotZpRUeHq6mTZuqf//+6tatm/bt26fo6OhTGomuWLGiBg4cqAkTJuiHH36QFJj6vWDBAlWqVCnT72Xa9yHplEa+TxUhGABCaNKkiZo0aZJtnXbt2ql///7q06ePPv/8c3Xo0EHnnHOO1qxZo5deekmVK1fOMgRfccUV6tixo2rXri0z07hx4xQZGan69euHrN+rVy89/vjjatu2rQYPHqyDBw9q8ODBiomJOalvbKX//lHSq1cv9enTR6tWrdJTTz2VaSrWiRQqVEijR4/WzTffrN9//13XXnutihQpElxsY/r06cE/OOrVq6eqVauqb9++KlasmNq2bZvuWCNHjlSzZs1UqFAhXX/99SpRooR+/fVXzZkzR8OGDVP16tVz3K6aNWtq/vz5mj9/vsqUKaMqVapo+vTpWrp0qVq3bq3y5cvrp59+0ttvv53jRcB87B5Jk82siKRfJPVS4FKqaWbWR9Kvkrp4dedKaiNpnaQDXl0A+eyiiy5Shw4d1LJlS8XGxmrjxo166qmnVLx48Uxf2KZ1on5pyJAhql+/vpKSktR6QALEAAAgAElEQVS7d2+VLVtWW7Zs0cKFC9WzZ8+QfWaNGjV0xx13qGvXrkpOTlZiYqIOHTqkVatW6ccff9Qrr7yiQoUKBUedu3fvrptuuklmpo8++kg33XSTEhMTVbNmTe3YsUMTJ05U7dq1VbZsWVWuXDlHfcd9992n559/Xm3btlX//v21detWjRgxQsWKFcvVz71///568cUX1apVKw0YMEB//vmnBg4cqIsuukidO3fO1XOlCtXflSlTRkOGDNH//u//auPGjbrqqqt0/Phx/fjjj1q8eLFmzJih7777Tg8++KBuvPFGnXfeefrtt9/05JNP6uKLLw6OJGd17Ixuv/324OJhJUuW1OLFi/XTTz/pySeflCTdcssteumll9SkSRM9+OCDOu+887Rnzx4tW7ZM8fHxuv/++1WkSBFVqVJF06ZNU+3atRUREaE6derk6b2XCcEAcBpGjx6tyy+/XM8995y6deumgwcPqnLlymrfvr0efPDBLPdr1KiRJk6cqA0bNigsLEyXXnqpPvjgA1WoUCFk/eLFi2vevHm68847deONNwY7/+Tk5GxX+wzloosu0oQJEzR06FDNmDFDF198sd5+++10K4Lm1I033qioqCgNHz5c48ePV1hYmM477zy1bds2U+d14403atiwYeratWu6BUqkwB9fn3zyiQYPHqybb75ZKSkpOvfcc9W6descj1ykevjhh4MLdv3++++aMGGC6tSpo9mzZ6t///7au3evypUrp9tuu02PPfbYSb9nP3HOfSMpMcSmZiHqOkn98rxRALL16KOPatasWbr33nu1d+9excfHB1dQTr1MJZQT9UvVq1fXl19+qYcfflh9+/bVwYMHgysxV61aNcvjPv/886pevbrGjRunRx99VFFRUbrgggvUp0+fYJ1u3bopIiJCw4YN0/XXXx8c5U1d9+KGG27Q4sWLlZycrF27dqlHjx6aOHFijvqOhIQEzZ07V/fee686d+6sWrVq6Y033lCHDh1y4+MOiomJ0eLFi/XAAw/opptuUpEiRdSmTRuNGTMmz8JcqP6uZ8+eSk5OVvny5TVmzBiNHj1aERERql69erCfj4+PV1xcnIYNG6atW7eqVKlSuuaaa4LBNbtjZ9SoUSONGzdOL7/8sg4dOqSqVatq3LhxwXU/IiIitHjxYj366KMaPHiwduzYodjYWNWvXz/dCuAvvfSSHnzwQTVv3lyHDx/W+vXrg/dmzguW3X29/k4SExNdQa5AhoBQC2FteOLEUzVwZlu9erVq1apV0M3wnfXr16t69eoaO3bsSa+qiZzJ7nfbzL52zoUKiMgh+macabLtz46EuK9skeKZywDkqdzomxkJBoCzxIgRI1S+fHmde+65+vXXXzVixAjFxMTk2TQrAACAvyNCMACcJcxMQ4cO1datW1W0aFFdeeWVeuqpp056OjQAAICfEYIB4CwxcOBADRyY3W1aAQAAcCInt6QoAAAAAABnMUIwAAAAAMA3CMEAcoVfVpqHf/A7DQDA3xMhGMBpCw8P18GDBwu6GUCuOnjwoMLDwwu6GQAAIJcRggGcttjYWG3ZskUHDhxg9AxnPeecDhw4oC1btig2NragmwMAAHIZq0MDOG2pt+jZunWrjh49WsCtAU5feHi44uLiuP0UAAB/Q4RgALkiKiqKwAAAAIAzHtOhAQAAgL8BMzvh4+OPP86Vc/3www8aMmSI/vzzz1w53onMnTtXzz33XI7rb9iwQeecc442bdoULPvll1/UrVs3VaxYUREREapUqZI6duyopUuXnlbb4uPj033GsbGxateunVatWpWuXp8+fdSvX7/TOhdyByPBAAAAQHa2/kc6eihzeXhE3p63/KUnVT1tmDt48KCaNm2qhx9+WElJScHyCy64IFea9sMPP2jo0KG64447dM455+TKMbMzd+5cLVq0SHfffXeO6g8dOlRdunRRxYoVJUm7du1SgwYNVKVKFY0aNUpxcXFav369Zs6cqa+++kqNGjU6rfb17NlTt99+u5xz2rJlix5//HG1atVKq1evVokSJSRJycnJuuSSSzRgwABVqlTptM6H00MIBgAAAP4GGjZsGHyeOkJ7/vnnpyv3gz179ujNN9/UokWLgmVTp07Vvn37NH/+fJUuXVqSdM0116h3797ZLuo5b948XXfddTp0KMSXIGkkJCSk+5zPO+881atXT8uXL1fTpk0lSTVq1FC9evX08ssva9iwYafzFnGamA4NAAAA+Mz69evVpUsXlSpVSpGRkUpKStLPP/8c3O6c02OPPabzzjtPERERio+PV5s2bbRnzx7NmzdPXbp0kSSVK1dOZqaaNWtmea5vv/1WLVq0UOnSpXXOOefowgsv1Lhx49LVmT59uurWrauIiAiVL19egwYNUkpKiiRp4MCBev7557V27drglOM77rgjy/NNmTJF0dHRuuKKK4Jl+/btU7FixVSyZMlM9c0sZx/aSUgd/c24YGjnzp31+uuv5/r5cHIYCQYAAAB8ZOfOnWrcuLESEhL0yiuvqEiRIho2bJhatmyp1atXq0iRIho3bpxGjx6tkSNHqlatWtq1a5cWLVqkgwcPqlGjRho+fLj++c9/as6cOYqOjlaxYsVCnuv48eNKSkpSYmKi3nzzTRUpUkSrV6/W/v37g3Vee+019erVS3fffbeeeOIJrV27Vv/85z9lZnr88cfVr18//fzzz1q+fLmmTp0qSYqLi8vy/X344Ydq2LBhunBbt25d/fHHH+rZs6f69++viy++OFfDr3NOx44dkyRt2bJFAwcOVGxsrBo3bpyu3uWXX67+/ftr7dq1qlGjRq6dHyeHEAwAAAD4yKhRo3T8+HEtWrQoODLaqFEjValSRa+//rr69OmjZcuWqW3btrr99tuD+3Xu3Dn4vFq1apIC4TI+Pj7Lc23dulVbtmzR4sWLg/s0a9YsuD0lJUUDBgxQ37599e9//1uS1LJlS4WFhSk5OVnJycmqWLGi4uLiFBERkaOp3V9//bV69OiRriwpKUl33XWXXnjhBb3++uuKiopSy5Yt1a9fPzVp0iRYzzkXHIGWAiFeUjDgSoGR47CwsHTHHz58uIYPHx58XaZMGc2cOTPT9dIXX3yxJGnZsmWE4ALEdGgAAADARxYtWqTWrVsrMjJSx44d07Fjx1S6dGldfPHFWrFihSTpkksu0cyZM/XYY49pxYoVwTB4suLi4hQfH6/bbrtNb7/9tnbt2pVu+8qVK7V9+3Z16dIl2JZjx46padOm+uuvv7R69eqTPueOHTtUtmzZTOWpU6qffPJJXXnllZozZ46aNm2qiRMnBuu8/PLLCg8PDz6SkpJ0+PDhdGWRkZGZjt27d28tX75cy5cv17x589SqVSt16NAhU/sjIiIUGRmp7du3n/T7Qu4hBAMAAACpjhzI/Ai1MvRZbPfu3Zo0aVK6YBceHq4vvvgieEuhO++8U4MHD9bkyZN12WWXKT4+XkOHDj3pMBweHq6FCxeqVKlS6tGjh+Lj49WkSRN9//33wbZIgdHhtG2pVauWJKW7xVFOpKSk6OjRoypatGjI7dWrV1dycrLef/99rV+/XhdeeKEeeuih4PbOnTsHw+zy5cv1zDPPqEiRIunKvvjii0zHLVeunBITE5WYmKhWrVrp9ddfV0xMTMgFsIoWLXrChbaQt5gODQAAAPhIdHS0GjZsqAEDBmTaljo9Ou105I0bN+q1117T4MGDde6556pnz54ndb7atWtr5syZOnLkiJYsWaLk5GS1a9dOGzZsUHR0tCRp0qRJIW/fdP7555/UucLCwhQVFaV9+/adsG5cXJxuueUWJScna//+/SpZsqRiYmIUExMTrLN7926ZmRITE0+qHYUKFVKNGjUyjQQ757R///7g+0bBIAQDAAAAPtKsWTN98MEHqlOnjooUKXLC+ueee64eeeQRvfLKK/rhhx8kKbjfyYxoFilSRC1atNC9996r3r1766+//tJFF12kmJgYbdy4Ubfccku2++b0XDVq1ND69evTle3cuVOxsbGZ6v7000+KjIwMOcX5dBw/flyrV6/OFOw3b96slJQUVa9ePVfPh5NDCAYAAAB8JDk5WVOnTlWzZs3Ur18/lStXTtu3b9fHH3+s5s2bq3PnzurVq5cSEhJUv359RUVFacGCBdq0aZOuueYaSQreEumFF15Q586dg7c+ymjZsmUaPHiwbrjhBlWpUkW7d+/W6NGj1aBBg2DwHDVqlG677Tbt3btXLVu2VOHChfXzzz9rxowZmjt3rsLCwlSzZk1t2rRJkydPVo0aNRQbG6tKlSqFfH+NGzfWJ598kq5s7Nixmjlzpm6++WbVqVNHR44c0bx58/TKK6/ogQceUOHCpxeLtmzZoi+//FKStHfvXk2aNEk//fRTcLGvVCtWrFBYWJjv7t18piEEAwAAAD4SHx+vr776SoMGDdK9996r33//XeXKldNVV12l2rVrSwrcymf8+PF6/vnndeTIEVWrVk0TJ07UtddeKylwbe3w4cP14osvavTo0apWrZrWrFmT6VwJCQkqXbq0HnvsMW3btk2lS5dW8+bN9cQTTwTr9OjRQ9HR0RoxYoRefvllFS5cWFWrVlW7du1UqFBgCaPu3bvr008/1X333afdu3fr9ttv10svvRTy/XXq1En//ve/tX379uDK1R06dNCWLVv00ksvadOmTQoPD1fVqlX18ssvq0+fPqf9mU6cODG4wFapUqVUq1YtzZo1K/h5pZo3b55atGgRvI8wCoY55wq6DfkiMTHRpa52h4JTeeCcTGUbnkgqgJYAwOkxs6+dcyd3kRjSoW/GmWb16tWqdf65Od+hSPG8awxOS82aNdWvXz/dc889Bd2UoKNHjyohIUEvvPCCrr/++oJuzllr9erVwYXTMspp38zq0AAAAAD+VgYNGqRnn332lG/tlBcmT56smJgYderUqaCb4nv5GoLNbLyZ7TSzlSG2PWhmzszKeq/NzJ4xs3Vm9p2Z1U1Tt4eZ/eQ9emQ8FgAAAAD/6t69u3r16qVt27YVdFOCwsLCNG7cuOAUbxSc/L4meKKk5yS9lrbQzCpKaiHp1zTF10qq5j0aSHpRUgMzi5Y0WFKiJCfpazOb7Zz7Lc9bDwAAAOCMV6hQoXT3/z0T3HzzzQXdBHjy9WsI59wnkvaG2DRGUrICoTZVB0mvuYAvJZUys3KSWkla6Jzb6wXfhZJa53HTAQAAAAB/AwU+Fm9m7SVtcc59m2FTgqRNaV5v9sqyKgcAAAAAIFsFeoskMysuaZCklqE2hyhz2ZSHOn5fSX0lZXkfMQAAACCVc05mof7cBFDQcuvORgU9Eny+pCqSvjWzDZIqSPp/ZhavwAhvxTR1K0jamk15Js65sc65ROdcYkxMTB40HwAAAH8X4eHhOnjocEE3A0AWDh48qPDw8NM+ToGGYOfc9865WOdcZedcZQUCbl3n3HZJsyXd4q0S3VDSfufcNknzJbU0s9JmVlqBUeT5BfUeAAAA8PcQGxurLVu36sDBQ7k24gTg9DnndODAAW3ZskWxsbGnfbx8nQ5tZlMkNZFU1sw2SxrsnHs1i+pzJbWRtE7SAUm9JMk5t9fM/k/Scq/eY865UIttAQAAADkWFRUl7TmirVs26WhObi8bViTP2wQgIDw8XHFxcYH/T09TvoZg59xNJ9heOc1zJ6lfFvXGSxqfq40DAACA70UVNUUVzUkCllS+Vt42BkCeKOhrggEAAAAAyDeEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPhG4YJuAP6+Kg+cU9BNAAAAAIB0GAkGAAAAAPgGIRgAAAAA4BuEYAAAAACAb3BNMApcqGuHNzyRVAAtAQAAAPB3x0gwAAAAAMA3CMEAAAAAAN8gBAMAAAAAfIMQDAAAAADwDUIwAAAAAMA3CMEAAAAAAN8gBAMAAAAAfIMQDAAAAADwDUIwAAAAAMA3CMEAAAAAAN8gBAMAAAAAfIMQDAAAAADwDUIwAAAAAMA3CMEAACAkM9tgZt+b2TdmtsIrizazhWb2k/eztFduZvaMma0zs+/MrG7Bth4AgNAIwQAAIDvXOOcucc4leq8HSvrQOVdN0ofea0m6VlI179FX0ov53lIAAHKAEAwAAE5GB0mTvOeTJF2Xpvw1F/ClpFJmVq4gGggAQHYIwQAAICtO0gIz+9rM+nplcc65bZLk/Yz1yhMkbUqz72avDACAM0rhgm4AAAA4YzV2zm01s1hJC81sTTZ1LUSZy1QpEKb7SlKlSpVyp5UAAJwERoIBAEBIzrmt3s+dkmZIqi9pR+o0Z+/nTq/6ZkkV0+xeQdLWEMcc65xLdM4lxsTE5GXzAQAIiRAMAAAyMbNIMyuR+lxSS0krJc2W1MOr1kPSLO/5bEm3eKtEN5S0P3XaNAAAZxKmQwMAgFDiJM0wMynw98Kbzrl5ZrZc0jQz6yPpV0ldvPpzJbWRtE7SAUm98r/JAACcGCEYAABk4pz7RdLFIcr3SGoWotxJ6pcPTQMA4LQwHRoAAAAA4BuEYAAAAACAbzAdGgAAAMjg2Y/WZSq7p2nVAmgJgNzGSDAAAAAAwDcIwQAAAAAA3yAEAwAAAAB8gxAMAAAAAPANQjAAAAAAwDcIwQAAAAAA3yAEAwAAAAB8I19DsJmNN7OdZrYyTdkoM1tjZt+Z2QwzK5Vm20Nmts7M1ppZqzTlrb2ydWY2MD/fAwAAAADg7JXfI8ETJbXOULZQUm3nXB1JP0p6SJLM7AJJXSVd6O3zgpmFmVmYpOclXSvpAkk3eXUBAAAAAMhWvoZg59wnkvZmKFvgnDvmvfxSUgXveQdJU51zh51z6yWtk1Tfe6xzzv3inDsiaapXFwAAAACAbJ1p1wT3lvSB9zxB0qY02zZ7ZVmVAwAAAACQrTMmBJvZIEnHJE1OLQpRzWVTHuqYfc1shZmt2LVrV+40FAAAAABw1jojQrCZ9ZDUVlJ351xqoN0sqWKaahUkbc2mPBPn3FjnXKJzLjEmJib3Gw4AAAAAOKsUeAg2s9aSBkhq75w7kGbTbEldzayomVWRVE3SMknLJVUzsypmVkSBxbNm53e7AQAAAABnn8L5eTIzmyKpiaSyZrZZ0mAFVoMuKmmhmUnSl865O5xzq8xsmqQfFJgm3c85l+Id525J8yWFSRrvnFuVn+8DAAAAAHB2ytcQ7Jy7KUTxq9nUHyZpWIjyuZLm5mLTAAAAAAA+UODToQEAAAAAyC+EYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuEYAAAAACAbxCCAQAAAAC+QQgGAAAAAPgGIRgAAAAA4BuFC7oBAAAAwNng2Y/WpXu9tXhhjeh0UQG1BsCpYiQYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAWTKzMDP7j5m9772uYmZfmdlPZvaWmRXxyot6r9d52ysXZLsBAMgKIRgAAGTnfyWtTvP6SUljnHPVJP0mqY9X3kfSb865qpLGePUAADjjEIIBAEBIZlZBUpKkV7zXJqmppOlelUmSrvOed/Bey9vezKsPAMAZhRAMAACy8rSkZEnHvddlJO1zzh3zXm+WlOA9T5C0SZK87fu9+gAAnFEKF3QDgFAqD5wTsnzDE0n53BIA8Cczaytpp3PuazNrklocoqrLwba0x+0rqa8kVapUKRdaCgDAycnXkWAzG29mO81sZZqyaDNb6C2wsdDMSnvlZmbPeAtsfGdmddPs08Or/5OZ9cjP9wAAgE80ltTezDZImqrANOinJZUys9Qv0StI2uo93yypoiR520tK2pvxoM65sc65ROdcYkxMTN6+AwAAQsjv6dATJbXOUDZQ0ofeAhsfeq8l6VpJ1bxHX0kvSoHQLGmwpAaS6ksanBqcAQBA7nDOPeScq+Ccqyypq6SPnHPdJS2WdL1XrYekWd7z2d5reds/cs5lGgkGAKCg5WsIds59oszfCqddSCPjAhuvuYAvFfjmuZykVpIWOuf2Oud+k7RQmYM1AADIGwMk9TezdQpc8/uqV/6qpDJeeX/990ttAADOKGfCNcFxzrltkuSc22ZmsV55cIENT+riG1mVAwCAPOCc+1jSx97zXxSYiZWxziFJXfK1YQAAnIIzeXXorBbYyNHCG1Jg8Q0zW2FmK3bt2pWrjQMAAAAAnH3OhBC8w5vmLO/nTq88uMCGJ3XxjazKM2HxDQAAAABAWmdCCE67kEbGBTZu8VaJbihpvzdter6klmZW2lsQq6VXBgAAAABAtvL1mmAzmyKpiaSyZrZZgVWen5A0zcz6SPpV/72eaK6kNpLWSTogqZckOef2mtn/SVru1XvMOZfpFgwAAAAAAGSUryHYOXdTFpuahajrJPXL4jjjJY3PxaYBAAAAAHzgTJgODQAAAABAviAEAwAAAAB8gxAMAAAAAPANQjAAAAAAwDcIwQAAAAAA3yAEAwAAAAB8I19vkYS/p8oD5xR0EwAAAAAgRxgJBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvkEIBgAAAAD4BiEYAAAAAOAbhGAAAAAAgG8QggEAAAAAvlG4oBsAAAAAwGe2/idzWflL878d8KVTHgk2s9JmdomZFc3NBgEAgNxHvw0AQECOQrCZDTWzJ9K8birpV0lfS/rZzC7Mo/YBAICTRL8NAEDWcjoS3F3SmjSvR0v6TFJjSWsljcjldgEAgFNHvw0AQBZyGoLLS/pFksysoqSLJQ12zn0p6V+SGuZN8wAAwCmg3wYAIAs5DcF/SCrpPW8q6Tfn3DLv9SFJxXO7YQAA4JTRbwMAkIWcrg69RNJAMzsu6UFJs9Jsqy5pU243DAAAnDL6bQAAspDTkeD7JR2WNFXSPkmD0my7RdInudwuAABw6ui3AQDIQo5Ggp1zWxSYThVKK0kHc61FAADgtNBvAwCQtZzeIukjM6uZxeZ4SfNzr0kAAOB00G8DAJC1nE6HbiIpKottUZKuypXWAACA3NBE9NsAAISU0xAsSS5jgZkVUWC61fZcaxEAAMgN9NsAAISQ5TXBZjZY0qPeSyfpSzPLqvqoXG4XAAA4CfTbAADkTHYLY82VtFuSSXpG0mhJGzLUOSJpjXPu0zxpHQD4XOWBc0KWb3giKZ9bgrMA/TYAADmQZQh2zi2XtFySzOwPSXOcc7vzq2EAACDn6LcBAMiZnN4iaVJeNwQAgP/f3r3Hy1mVhx7/PRLuIOESkJBgQPBuj2CkWLVasBbEGrRc4g2kWKwfoNpqa+T0qD2nHqG1XqrnYGkRwaNcVDxEQC0SrG2PUCIgiKgEiJAmQJSbGhGQ5/yx1iaTyczes/ee2TOz5/f9fOYzM+t9Z+Z5LzNrnvdda73qDuttSQNl7fX9jkDaREdJcERsCbwDeB2wANimeZ7M3L27oUmSpKmw3pYkqb2OkmDgo8DbgEuBqyh9iiRJ0mCy3pYkqY1Ok+CjgWWZ+Xe9DEaSJHWF9bYkSW10ep3gAG7sZSCSJKlrrLclSWqj0yT4H4HX9zIQSZLUNdbbkiS10Wlz6HuAN0bEVcAVwANN0zMzz+xqZJIkaaqstyVJaqPTJPhj9X5v4GUtpidgZSpJ0mCw3pYkqY1OrxPcabNpSZLUZ9bbkiS1ZyUpSZIkSRoZnTaHJiJ2B94FLAYWAq/NzJsj4h3Af2Tmt3sUoySNhEXLLut3CJpFrLclSWqtozPBEXEQcCvwB8Bq4GnA1nXynpRKVpIkDQDrbUmS2uu0OfRHgauApwNvo1x/cMx/AAd1OS5JkjR11tuSJLXRaXPoA4Elmfl4RETTtJ8Cu3c3LEmSNA3W25IktdHpmeAHgXltpu1LuR6hJEkaDNbbkiS10WkSfAnwVxGxb0NZRsRuwLuBi6cbSET8aUTcHBHfi4jzI2KbiNgnIq6JiFsj4sKI2KrOu3V9vqpOXzTdz5ckaRbpeb0tSdKw6jQJXgY8BHwf+FYt+xTwQ+CXwPumE0RE7AX8CbA4M58LbAEsBc4APpqZ+wP3AyfWl5wI3J+Z+1H6PZ0xnc+XJGmW6Wm9LUnSMOsoCc7M+4GDgZOBHwPfAO6gVLIvzsyfdSGWOcC2ETEH2A5YBxwCfLFOPxc4sj5eUp9Tpx/aos+TJEkjaYbqbUmShlLH1wnOzEeAs+utqzLzPyPiw8CdlCPU/wx8B3ggMx+rs60B9qqP9wLuqq99LCIeBHYFftLt2CRJGka9rLclSRpmnV4n+F8i4u0R0W6QjWmJiJ0pZ3f3AeYD2wOHt5g1x14yzrTG9z0pIlZGxMr169d3K1xJkgZar+ttSZKGWad9gtcDHwbWRsQVEfGHNXHtllcAd2Tm+sx8lDJgx28Bc2vzaIAFwNr6eA2wEKBO3wm4r/lNM/OszFycmYvnzfN/gCRpZPS63pYkaWh12if4KMo1BY8Dfg78L+DuiLg0It4cETtOM447gYMjYrvat/dQymAeVwFH1XmOp4x2CbC8PqdOX5GZm50JliRpFM1AvS1J0tDq9EwwmfmLzDw/M19LqVj/qE76R+Du6QSRmddQBri6DripxnUW8B7gzyJiFaXP78HeHFwAACAASURBVFi/prOBXWv5n1EG+pAkSVUv621JkoZZxwNjNcrMn0XEbZSRJh8CdptuIJn5fuD9TcW3Awe1mPdh4OjpfqYkSaOgF/W2JEnDquMzwQARcVBE/F1E3Em57uDLgI8D+/ciOEmSNHXW25Ikba6jM8ERcTpwDPBU4FbgHOCCzLylh7FJkqQp6Ea9HRHbUBLnrSn/F76Yme+PiH2AC4BdKN2Y3pyZj0TE1sB5wAuAnwLHZubq7i2VJEnd0Wlz6GOAiygV6A09jEeSJE1fN+rtXwGHZObPI2JL4N8i4quUsTg+mpkXRMSngBOBM+v9/Zm5X0QsBc4Ajp32kkiS1GUdJcGZuW+vA5EkSd3RjXq7XnXh5/XplvWWwCHAG2r5ucAHKEnwkvoYymCXn4yI8OoNkqRB03Gf4IjYOiLeHhFnR8Q/R8T+tfzYiHhW70KUJEmT1Y16OyK2iIgbgHuBK4DbgAcy87E6yxpgr/p4L+AugDr9QcqVHSRJGigdJcER8XTgR8CHgEWU6/iOXWPwpcB7exGcJEmavG7V25n568x8PrCAcrWGVsnz2JneGGdaY2wnRcTKiFi5fv36TsKQJKmrOj0T/PfAnZSK9PfYtKL7F+Al3Q1LkiRNQ1fr7cx8APgmcDAwNyLGulMtANbWx2uAhQB1+k7AfS3e66zMXJyZi+fNmzeZMCRJ6opOk+CXAh+qlWDzUd17gD27GpUkSZqOadfbETEvIubWx9sCrwBuAa4CjqqzHQ9cUh8vr8+p01fYH1iSNIg6HR36YWDbNtP2Ah7oTjiSJKkLulFv7wmcGxFbUA6aX5SZl0bE94ELIuKvgeuBs+v8ZwOfjYhVlDPAS6ezAJIk9UqnSfAVwGkR8Q02jhSZ9ZqApwKX9yI4SZI0JdOutzPzRuCAFuW3U/oHN5c/DBw9naAlSZoJnSbBfw78O7CKUrEm8D7gOcBWwOt6Ep0kSZoK621JktroqE9wZt4F/BfgU5RBNm6jNJP6AvCCzLy7VwFKkqTJsd6WJKm9Ts8Ek5n3A/+t3iRJ0gCz3pYkqbVOR4eWJEmSJGnomQRLkiRJkkaGSbAkSZIkaWSYBEuSJEmSRoZJsCRJkiRpZHQ8OnREHAcsz8wHehiPJI2ERcsu63cImuWstyVJam0yZ4LPAfYGiOJ9EfGU3oQlSZKmyXpbkqQW2p4JjojLgO/W241AAFknPwl4P3ApcHePY5QkSROw3pY0LD6xYtVmZacesl8fItGoGq859BXAAcCrgWdSKtJPRsRVwLVsWrlKkqT+st6WJKkDbZPgzPzY2OOI2Br4JXAd8AzgzZSK9LMR8TXgG5n5tR7HKkmS2rDeliSpM237BEfEqRHxkojYMTN/VYvPyczXUyrUAM4HdgA+2ftQJUlSO9bbkiR1Zrzm0K8B/hLYLSJWU44gL42IbYGb6jxfzczrehuiJEnqgPW2JEkdaHsmODN/NzP3ABYAp1COIL8C+CpwH6VyfXtEHFqbXUmSpD6x3pYkqTMTXiIpM9dl5lfr07dm5i7AYkrluhD4DHB/zyKUJEkds96WJGl8k7lOcKNb6v1pmbkQeEGX4pEkSd1nvS1JUjVen+BNZGZjwpzAj4Ff1Wm3tHyR1GWLll22Wdnq04/oQySSNNistyVJaq3jJLhRZj4O7NPlWCRJUg9Yb0uStNGUkmBJUn/ZKkKSJGlqptonWJIkSZKkoWMSLEmSJEkaGSbBkiRJkqSRYRIsSZIkSRoZDowlSZKkkffei28CYP6GVX2ORFKveSZYkiRJkjQyTIIlSZIkSSPDJFiSJEmSNDJMgiVJkiRJI8MkWJIkSZI0MkyCJUmSJEkjw0skSZIkSeqZT6zwslMaLJ4JliRJkiSNDJNgSZIkSdLIMAmWJEmSJI2MgUmCI2JuRHwxIn4QEbdExIsiYpeIuCIibq33O9d5IyL+PiJWRcSNEXFgv+OXJEmSJA2+gUmCgY8DX8vMZwL/BbgFWAZcmZn7A1fW5wCHA/vX20nAmTMfriRJkiRp2AxEEhwRTwZ+GzgbIDMfycwHgCXAuXW2c4Ej6+MlwHlZXA3MjYg9ZzhsSZIkSdKQGYgkGNgXWA+cExHXR8Q/RcT2wB6ZuQ6g3u9e598LuKvh9Wtq2SYi4qSIWBkRK9evX9/bJZAkSZIkDbxBSYLnAAcCZ2bmAcAv2Nj0uZVoUZabFWSelZmLM3PxvHnzuhOpJEmSJGloDUoSvAZYk5nX1OdfpCTF94w1c6739zbMv7Dh9QuAtTMUqyRJkiRpSA1EEpyZdwN3RcQzatGhwPeB5cDxtex44JL6eDlwXB0l+mDgwbFm05IkSZIktTOn3wE0OBX4XERsBdwOnEBJ0i+KiBOBO4Gj67yXA68CVgEb6rySJEmSJI1rYJLgzLwBWNxi0qEt5k3g5J4HJUmSJEmaVQaiObQkSZIkSTNhYM4EazgsWnZZv0OQJEmSpCnzTLAkSZIkaWSYBEuSJEmSRoZJsCRJkiRpZNgnWJIkSZqC+Rt+AGsfayo8oD/BSOqYZ4IlSZIkSSPDJFiSJEmSNDJsDi1JkqTRtPb6Jx7O37Cqj4FImkmeCZYkSZIkjQyTYEmSJEnSyDAJliRJkiSNDJNgSZIkSdLIMAmWJEmSJI0Mk2BJkiRJ0sgwCZYkSZuJiIURcVVE3BIRN0fEO2r5LhFxRUTcWu93ruUREX8fEasi4saIOLC/SyBJUmsmwZIkqZXHgHdl5rOAg4GTI+LZwDLgyszcH7iyPgc4HNi/3k4Czpz5kCVJmphJsCRJ2kxmrsvM6+rjnwG3AHsBS4Bz62znAkfWx0uA87K4GpgbEXvOcNiSJE1oTr8DkCR1x6Jll21Wtvr0I/oQiWabiFgEHABcA+yRmeugJMoRsXudbS/groaXrall62YuUkmSJuaZYEmS1FZE7AB8CXhnZj403qwtyrLF+50UESsjYuX69eu7FaYkSR0zCZYkSS1FxJaUBPhzmXlxLb5nrJlzvb+3lq8BFja8fAGwtvk9M/OszFycmYvnzZvXu+AlSWrDJFiSJG0mIgI4G7glMz/SMGk5cHx9fDxwSUP5cXWU6IOBB8eaTUuSNEjsEyxJPdaqr640BF4MvBm4KSJuqGWnAacDF0XEicCdwNF12uXAq4BVwAbghJkNV5KkzpgES5KkzWTmv9G6ny/AoS3mT+DkngYlSVIX2BxakiRJkjQyTIIlSZIkSSPDJFiSJEmSNDJMgiVJkiRJI8MkWJIkSZI0MkyCJUmSJEkjwyRYkiRJkjQyTIIlSZIkSSPDJFiSJEmSNDJMgiVJkiRJI2NOvwOQJEmSJNZe37p8/gEzG4dmPc8ES5IkSZJGhkmwJEmSJGlkmARLkiRJkkaGSbAkSZIkaWSYBEuSJEmSRoZJsCRJkiRpZJgES5IkSZJGhkmwJEmSJGlkmARLkiRJkkaGSbAkSZIkaWQMVBIcEVtExPURcWl9vk9EXBMRt0bEhRGxVS3fuj5fVacv6mfckiRJkqThMKffATR5B3AL8OT6/Azgo5l5QUR8CjgROLPe35+Z+0XE0jrfsf0IWJIkSdL0fGLFqk2en3rIfn2KRKNgYM4ER8QC4Ajgn+rzAA4BvlhnORc4sj5eUp9Tpx9a55ckSZIkqa1BOhP8MeAvgB3r812BBzLzsfp8DbBXfbwXcBdAZj4WEQ/W+X8yc+FKkiRJ2sTa6/sdgTShgUiCI+LVwL2Z+Z2IePlYcYtZs4Npje97EnASwN57792FSDWIFi27bLOy1acf0YdIJEmSJA26gUiCgRcDr4mIVwHbUPoEfwyYGxFz6tngBcDaOv8aYCGwJiLmADsB9zW/aWaeBZwFsHjx4s2SZEmSJGk6NuvL+qYD+hSJpE4NRJ/gzHxvZi7IzEXAUmBFZr4RuAo4qs52PHBJfby8PqdOX5GZJrmSJEmSpHENypngdt4DXBARfw1cD5xdy88GPhsRqyhngJf2KT5JkiR1ol1f0fmeOZU0swYuCc7MbwLfrI9vBw5qMc/DwNEzGpgkSZIkaegNRHNoSZIkSZJmwsCdCZYkSdIAsjmzpFnCM8GSJEmSpJHhmWBJkiR1V7uzxpI0AEyC1daiZZf1OwRJkiRJ6iqbQ0uSJEmSRoZJsCRJkiRpZJgES5IkSZJGhkmwJEmSJGlkmARLkiRJkkaGo0NLkiRpU17iSNIsZhIsSbNYu0udrT79iBmORJIkaTCYBEtSl3htbUmSpMFnEixJkjTbtGvOPP+AmfssSRpQDowlSZIkSRoZJsGSJEmSpJFhc2hJkqRRZnNmSSPGJFiSJEmaLVod1OhFX3BpiNkcWpIkSZI0MkyCJUmSJEkjw+bQkiRJ6h+b70qaYZ4JliRJkiSNDM8ES5IkjQpHgpYkzwRLkiRJkkaHSbAkSZIkaWSYBEuSJEmSRoZJsCRJkiRpZDgwliRJkkbKey++CYD5G1b1ORJJ/eCZYEmSJEnSyPBMsCRJ0jDzskeSNCmeCZYkSZIkjQzPBEuSpJYi4tPAq4F7M/O5tWwX4EJgEbAaOCYz74+IAD4OvArYALwlM6/rR9yS+ucTK+xnrcHnmWBJktTOZ4DDmsqWAVdm5v7AlfU5wOHA/vV2EnDmDMUoSdKkeCZYkiS1lJnfiohFTcVLgJfXx+cC3wTeU8vPy8wEro6IuRGxZ2aum5lopRFkf3BpSjwTLEmSJmOPscS23u9ey/cC7mqYb00t20REnBQRKyNi5fr163serCRJzTwTLAAWLbus3yFIkoZbtCjLzQoyzwLOAli8ePFm0yVJ6jXPBEuSpMm4JyL2BKj399byNcDChvkWAGtnODZJkiZkEixJkiZjOXB8fXw8cElD+XFRHAw8aH9gSdIgsjm0JElqKSLOpwyCtVtErAHeD5wOXBQRJwJ3AkfX2S+nXB5pFeUSSSfMeMCSJHXAJFiSJLWUma9vM+nQFvMmcHJvI5KjAUvS9JkEa1ZqN9DX6tOPmOFIJEmSNBRaHWSaf8DMx6GeMwmWJEmSZjOTu6lz3c1KDowlSZIkSRoZJsGSJEmSpJFhEixJkiRJGhkD0Sc4IhYC5wFPAR4HzsrMj0fELsCFwCJgNXBMZt4fEQF8nHIphg3AWzLzun7ELkmSJEkdj95un+K+G5QzwY8B78rMZwEHAydHxLOBZcCVmbk/cGV9DnA4sH+9nQScOfMhS5IkSZKGzUAkwZm5buxMbmb+DLgF2AtYApxbZzsXOLI+XgKcl8XVwNyI2HOGw5YkSZIkDZmBaA7dKCIWAQcA1wB7ZOY6KIlyROxeZ9sLuKvhZWtq2bqZi1SSJElDo6Gp6vwNq/oYiKR+G4gzwWMiYgfgS8A7M/Oh8WZtUZYt3u+kiFgZESvXr1/frTAlSZIkSUNqYJLgiNiSkgB/LjMvrsX3jDVzrvf31vI1wMKGly8A1ja/Z2aelZmLM3PxvHnzehe8JEmSJGkoDERz6Dra89nALZn5kYZJy4HjgdPr/SUN5adExAXAbwIPjjWbliRNbNGyyzYrW336EX2IRJKkIdPpKNAaWAORBAMvBt4M3BQRN9Sy0yjJ70URcSJwJ3B0nXY55fJIqyiXSDphZsOVJEmSJA2jgUiCM/PfaN3PF+DQFvMncHJPg5KkcbQ6kypJkqTBNxBJsCRJkiSNhHbNqecfMLNxjLCBGRhLkiRJkqReMwmWJEmSJI0Mm0NLkiRJ3WJTV01Vq33H/aYnPBMsSZIkSRoZJsGSJEmSpJFhc+gR42VdJEmSJI0yk2BJkiRJA+UTK1Y98XjtdiVl+dDrntevcDTL2BxakiRJkjQyTIIlSZIkSSPDJFiSJEmSNDLsEyxJkqTB4vVSJfWQSbAkSdIgapUIShotnf4OeJBoUmwOLUmSJEkaGZ4J1khpdZ3k1acf0YdIJEnSbNR4aZ8xpx6yXx8ikdSOSbAkSZKkgTV/ww/Kg7WPNRTa/FdTZ3NoSZIkSdLI8EywJE2gVTN6SdLweO/FNzF/w+bNlCWNJpNgSZIkaZA5UrjUVTaHliRJkiSNDM8Ez2I24ZQ0GY6eLkmSRoFJsCRJUj/NQFPX5sv2eMkeSaPMJFgjz7NfkiRJ0ugwCZYkSVJLnkGWNBuZBEuSJKlrmhNnMHmelRyxWkPMJFiSJGkWaZWEdjJPJ4lqt15nUix1WbuDEvMPmNk4hoRJsCRJkjpKniVpNjAJliRJ0pSZPA++2XImvnE51m43hw+97nl9jEbDzCRYktSWo6dLkjTLtGo6PWLNpk2CJUmSJE2arQA0rEyCJalBqzOfktQ10xxRd7aMvDxblkMaeI7i3ZJJsCRJ0hDzbJwauT9IEzMJliRJkobQjCW8nk3ULGMSPEvYhFOSJKlqStrmb+jv2dGSrG6M4Ymm3yM2GNFAMbEfaSbBkkaSB44kSVI79luf3Z7U7wAkSZIkSZopngmWJElS303pzJtNWiVNgUnwAGvXXHP16UfMcCSjp9W6d71LkiQNhvkbfgBrH2sqnHof614NMmaz6sFkEjyE7MvYHx6UkCRNimcpJWkgmQRLkiRpaHld3Jnjuu6O5vU4EGeGWx20m8Wjl5sES5pVbCnRe7aKkKTpGUuC1m638a/4h173vH6F09J7L77picf9vsSU1G0mwQPCP+7Dy/7DUuF3QVJPzcLm5fM3/GDjk7H+reOcffNM7OBzGw0Hk2CpB0wGJEmSNFUDMaBWuwNPs6CZtElwl5j0SJKk6fAM0uZcJ31QE59haALdvH+s3W5OT5uV96ov79Dt57Og/7BJcB/Y9FnqDr9Lg88DhFIxEGd1JEmASXBP+Qdd6g6/S5Jmo6E7+6OZMQv7Ps+Ubl87eCJ+h4fXUCfBEXEY8HFgC+CfMvP0mfhc/5BrKjwjJhWT+Q31OzJ8+lU3952Ji6ZpIC+bM4Q2X499CqQHupV0d7JvTbr1ypA1kR7aJDgitgD+F/C7wBrg2ohYnpnf729kUudGJRnwwJGmqtN9Z5i/H7OJdXPh2SFpcNgVYXP+Rg1xEgwcBKzKzNsBIuICYAkwUhWtJKk9W2DMuNlVN3t2V30004nK2HWBh2FArOkyCeyzAThrPMxJ8F7AXQ3P1wC/2adYpJ7zz7zUnq0NBkb/6uYB+FMlDYtPrFjF2u2GOQ3QIJjwLPsAH0gc5r0/WpTlJjNEnAScVJ/+PCJ+OM777Qb8pEuxDYLZtDwuSxtxRrfeacrcNoPJZRlHF783T+3aO80e3a6bu2U2fScGieu1+1ynveF67YI/2bxos/XaYp6Z1lHdPMxJ8BpgYcPzBcDaxhky8yzgrE7eLCJWZubi7oXXX7NpeVyWwTWblsdlGUyzaVlGRFfr5m5xP+oN12v3uU57w/XaG8O8Xp/U7wCm4Vpg/4jYJyK2ApYCy/sckyRJo8y6WZI08Ib2THBmPhYRpwBfp1yG4dOZeXOfw5IkaWRZN0uShsHQJsEAmXk5cHmX3m5Gm2bNgNm0PC7L4JpNy+OyDKbZtCwjoct1c7e4H/WG67X7XKe94XrtjaFdr5GZE88lSZIkSdIsMMx9giVJkiRJmpRZmQRHxKcj4t6I+F6b6TtFxFci4rsRcXNEnFDLnx8R365lN0bEsQ2v+UxE3BERN9Tb8wd9eeq0XzfEvLyhfJ+IuCYibo2IC+sAJgO7LBHxOw3LcUNEPBwRR9Zpfdk2HSzLzhHx5bov/UdEPLdh2mER8cOIWBURyxrKB3W7tFyWiFgYEVdFxC11e72j4TUfiIj/bNgur5qJZZnO8tRpqyPiphrzyobyXSLiirptroiInQd5WSLiGU3fmYci4p11Wl+2zXj7S8M8ERF/X78bN0bEgQ3Tjq/r/9aIOL6h/AV1m62qr211mR7Nch18V/68YZ//XpT6cZeZjnPYTLXe1vimU0+ptenWMWqtw/X6zCg51K8i4t39iHPSMnPW3YDfBg4Evtdm+mnAGfXxPOA+YCvg6cD+tXw+sA6YW59/BjhqmJanPv95m9dcBCytjz8FvH3Ql6Vhnl1q+Xb93DYdLMvfAu+vj58JXFkfbwHcBuxb97vvAs8e8O3Sbln2BA6sj3cEftSwLB8A3j3T22U6y1OfrwZ2a/GavwGW1cfLxvbTQV6Whnm2AO4GntrPbTPe/tIwz6uAr1KuN3swcE0t3wW4vd7vXB/vXKf9B/Ci+pqvAof3Y7/z1t/bRN+Vpnl/H1jR75iH4dbBb9CE9ba3Ka3XCX/bvW22zqZcx3ib9nrdHXgh8MF+/feb7G1WngnOzG9RfoTbzgLsWM8W7FDnfSwzf5SZt9b3WAvcS/lB76upLk+7met8hwBfrEXnAkd2J9rxdWlZjgK+mpkbehNlZzpYlmcDV9Z5fwAsiog9gIOAVZl5e2Y+AlwALBnw7dJyWTJzXWZeV8t/BtwC7NXreCcyjW0zniWUbQJDsG2a5jkUuC0zf9ybKDvT4f6yBDgvi6uBuRGxJ/B7wBWZeV9m3g9cARxWpz05M7+dpSY+jxnaNhosHXxXGr0eOL+H4cwa3f4PoqJH9dRIm2YdozY6Wa+ZeW9mXgs82ocQp2RWJsEd+CTwLGAtcBPwjsx8vHGGiDiIcpbutobiD9amEx+NiK1nLNqJjbc820TEyoi4OmrzYWBX4IHMHKuk1jAAiUs14bahXHey+c/LIG6b7wKvgyf2p6cCCyjr+q6G+cbW/yBvl3bL8oSIWAQcAFzTUHxK3S6fnqnmwx0ab3kS+OeI+E5EnNTwmj0ycx2UCoFy1HMQTLhtaP2d6eu2abO/QPvvx3jla1qUSy1FxHbAYcCX+h3LLNFJva3J6+S3XW1MoY5RB8ZZr0NnVJPg3wNuoDR5fj7wyYh48tjEekTos8AJDT/k76U0R3khpTnee2Y04vGNtzx7Z+Zi4A3AxyLiaZQmIM0GZZjwTrbN8yjXoBwzqNvmdGDniLgBOBW4nnJ0vN36H+Tt0m5ZAIiIHSh/KN+ZmQ/V4jOBp1G24zrg72Y04vGNtzwvzswDgcOBkyPit/sUY6cm2jZbAa8BvtDwmr5umzb7yxOTW7xkvO/HIH9vNJh+H/j3zOz0rLHGN269rSkb97dd7U2xjtEEJlivQ2dUk+ATgItrU4hVwB2UJIr6w30Z8Je1mQTwRFOAzMxfAedQmrQOirbLU5t1k5m3A9+kHL35CaX5x9h1ohdQjuAOgrbLUh0DfDkzn2huMajbJjMfyswTMvP5wHGUpvV3UI46LmyYdWz9D+x2GWdZiIgtKT+Kn8vMixtec09m/roeSPpHBmS7wPjL0/CduRf4MhvjvmesyVS9v3fGA29hvGWpDgeuy8x7Gl7Tt23Tbn9p0O77MV75ghblUjutWkZo6iaqtzUFHfy2q4Vp1DEaRwfrdeiMahJ8J6WPHLV/xTOA2+sZky9T+go0njWh4c9vUPqbtRzNr0/aLc/OY02DI2I34MXA92u/uasofWsBjgcumfGoW2u5LA3TN+vHNajbJiLmxsbRnd8KfKseObsW2D/KSNBbUf6QLR/k7dJuWeo6Pxu4JTM/0vSaxj42r2VAtguMuzzbR8SOdZ7tgVeyMe7llG0CQ7BtGmZp+52pZmzbjLe/NFgOHFdH8DwYeLA2P/868Mr6u7YzZdt8vU77WUQcXN//OAZk22jwRMROwMtwH+mmieptTUEHv+1qMs06Rm10uF6HTpT/3bNLRJwPvBzYDbgHeD+wJUBmfioi5lNGFN6T0izi9Mz8PxHxJsqZxJsb3u4tmXlDRKygHIULSrOfP87Mnw/48vwW8A/A45QDHh/LzLPre+5LGZBpF0oTmzfVM6kDuSz1tYuAfwcWNvY36te26WBZXkQZpOfXwPeBE7MM6EOUS9J8jDJq76cz84O1fFC3S8tliYiXAP9K6Qc2tk1Oy8zLI+KzlKZxSRlx+W0zVdFMY3n2pRwIA5gDfL5h2+xKGb17b8qfvqNnojnlNPez7Sh9n/bNzAcb3rMv26bd/kJZp2PLE5Q+hocBGyjdUlbW1/9hnR/gg5l5Ti1fTPnd2JYy6uepORsrN41rou9KnectwGGZubQ/UQ6f6dTbam86v+1qbbp1jFrrcL0+BVgJPLnO83PKCNIDe+BmVibBkiRJkiS1MqrNoSVJkiRJI8gkWJIkSZI0MkyCJUmSJEkjwyRYkiRJkjQyTIIlSZIkSSPDJFhdExEfiIhsuN0dEZdGxG/0O7bpiIit6rI9f4Y/d/f6uYu6+J4fjojVE8zTuB0fj4j7I+LaiPhgHQK/pyJidf3s/9pi2ksbYlvU61g6ERHvioirWpS/LCIuiYh7I+LRev+ViDi8XqKh0/e/NCJuGmf6J+s22joijo6IH0bEFlNdHkmzi3Vz1z/XunnzadbNm0+3bh5wJsHqtgeBF9XbO4GnA1dExC59jWp6tqJcv29GK1pg9/q5i2b4c2HjdvwtYClwMfBm4KaIeMEMfP7Pgde3KF9apw2EiNgBeA9welP5O4GrKNd3PBU4FDiFEvulwCGT+JjzgedGxHNafP4WwFHAxfV60l+iXKfzzZNeGEmzmXVz91g3b866edPPsW4eAibB6rbHMvPqersAOI5SYRzW57hmRERs2+8YuqRxO349Mz8E/AawDrhwBo5mXgo8OyKeO1bQUKks7/FnT8brgV8B/zxWEBEHAh8G/ntmvi4zL8zMb2XmRZn5euAlwE8m8RmXABsofzKa/Q6wB6UyJjMfB86jVO6SNMa6eXawbu6MdbMmZBKsXvtuvV/YWBgRu0TEP0TEPRHxcET8v4j4zaZ5toiI90bEjyLiVxGxJiI+0zTPKRFxa52+KiL+tGn6ByLiJxFxQERcHREbIuL6iHhp03yviYjvRMQvavOVayLiZXXyz+r9OY3NfeotI+KNEXFeRDwAfKW+X0bEKa1iaSp7akScX2PcEBE3RsQbanOisWY2V4197iTX39yI+HxdpnWtmjBNRmY+APwF8DTgdzt5NrnmzQAACe9JREFUTUQ8JyK+FhH31ThuiYiTO3jpfwL/xqaVyyHADrSoaGuzp2sj4sG6Tr4SEfs1zfOSiPjXiHio3m6IiKMbpo+3D7RzPOVIbzaUnQrcC/x1qxdk5rcz87uNZRHx1oi4ue7HP46Iv2iYf+wI9bEt3m4pcA/lyPaYLwEHNv5JkaQm1s1NsTSVWTe3Zt28cX7r5iFnEqxe27ve3zFWEBFbA9+g/Fj/OXAksB74Rmzar+UfgL8CLgJeDbwL2L7hff4I+ATlh/f3gS8AfxcRy5pi2A44t77fH1CODn45Irar7/M04IvAivo+b6T8sI01ExtrHvPXbGxOtq7h/T9MqYyPBv5nR2ulfO7uwLeBFwLvrp99NuVPyboaB8DJDZ87mfV3DnA4penbScAraX3EcjKuAh4DDu5w/uWUZkdvAl5D2V47dvja89k03tdT/sj8osW8C4BPAkuAPwK2AP49InYCiIgnU7bp7ZR94Cjgs8DcOn2ifWAzEbE98JvA/2ua9NvAisx8rJOFjIg/B84E/i9lPz8T+B9Nf9TOB/aPhuZuEbEl8Frgosz89Vh5Zt4C3E9p5iVJrVg3t2HdPCHr5o2sm4dZZnrz1pUb8AFKU5I59fY04ArgemDrhvlOBB4B9m8omwPcBvxtff5MIIE/afNZT6IckTynqfx/U/rMbNMQUwKHNMzz/Fp2WH1+FPDTcZZrhzr/W5rKF9XyL7d4TQKntFo/Dc8/RKk09mzzuc+t7/PypvJO1t9z6muPbVqO+4DVnWzHcaavA87sYH/YrcbwvEnuR6spf17mAY9S/ohsRak8jqRURgksavP6LYBtKX9+jqtli+trdmzzmnH3gTav+a36ns9pKv8l8KGmsmj4XswBnlTLn0zpi/T+pvn/O3A3sEV9vnVd/r9tmGdsPbyoRWzfBD43meXx5s3b7Lxh3dw4zbrZurmxzLp5hG+eCVa37Ur5cXwUWAUcALwuy8AAY14BfAe4IyLmRMScWv4vlB9EKP0pAD7T5nMWAPMpR5gbXUj58XpeQ9mjlB+eMd9veA8oTZt2iohzI+KV9SjiZFw2yfnHHAJ8LTPXTTjnpjpZfy+s9080T8rSdOeKKcbaqNPRE+8D7gI+FRHH1qPrHcvM9ZSjv0sp/dYC+GrLgCIOjogrIuKnlKPhGyh/LJ5eZ7mNUqF9PiKWRMTcpreYyj4wdmS/VR+ibHr+B2z8XjwK/E0tfxHlDMoXxrZl3Z4rKP2JFgDU78+XgWMinhi98ljgx8DVLT7/Jw3xSZJ1c+esm8dh3WzdPFuYBKvbHqT8yB8MvI1ylPDzEdG4r+1Wpz/adDuBjf2TdgV+kZkPtfmcPev9PU3lY88bm8o8lGVQAgAy85H6cJv6/IeUpjr7ApcDP6n9deZNuLStY+jUrmzadKtTnay/pwA/y8xfNr323qmFWkTENpS4J1zmus5fSTlq+mng7trv54BJfOQFwDHAG4D/2/SHbSymvSmDXwRln3sxZR+8l43b+P4ay5aUJnzrI+KyiNi3Tp/KPrBNvW+OaS0b/8SNubLG9EI23ea71fub2XRbjvUjauyvdz6lCeOL6nZYApyf9fByk181xCdJ1s2ds26emHXzRtbNQ2rOxLNIk/JYZq6sj6+JiF9SRsQ7mnIkGMpRyJXA21u8fuxH66fA9hHx5DaV7diPVfMRzD0aPqNjmXkZcFntp3IE8DFKH5lO+um0+6HbqqmsuQ/LT9n4h2EyOll/dwM7RsS2TZXtpI74tvA7lN+Nb3cyc2b+APiD2kfmpcAZlPW8oPHPzzguBj5F2X+OaDPPYZS+ZUsy8xcA9YjtJus7M78NHBZllNBXAB8BPk/tQzWFfWBsH5sLPNBQ/i3glRGxRdb+QLWiX1lje6TFe7ya1n9eftjweEWdZyllv9mROvJkC3OZ5HdA0qxm3VxYN2PdbN0s8Eyweu//UI6kvaeh7EpgP+DOzFzZdBsbdXFFvT+uzfuuoRzVO7qp/BjgITaO3jgpmflgZn6e0rzl2bV4k6PTHVoDPGvsST3a3nz9uSuB34uIPWit3ed2sv6urfevaYhhBzocObKV2kzpDEpTum9M5rWZ+WhmrqBUbntSB73o4HUP1s/80jifuS3wOKWp1ZhjaHOQLzN/mZlfoRwBf3aL6a32gVbGKsF9mso/QfnDd9o4rx3zbUo/pfkttuXKzBwb/ZRaaX+Bss+/AbglM29s876LgB918PmSRpN1M9bN1s1tWTePAM8Eq6cyMyPifwKfi4hDM/NKytHnPwa+GREfpowKuCtwEHB3Zn40M38YEWdRRpTcnXIEby5wVGYuzczHI+IDwD/UviZXAC+jHIE9LTMf7jTGiHgbpf/H1yiV9/6UH7Pz6jI8EhF3UPp8fA94GGj3Azfmy8DJEXF9Xb63UvpDNfoo5Y/Ev0bEByl9dJ4FbJ+ZfwPcSfkRPj4iHgQerUfyO1l/N0fEcuDMOvriOspolRs6XC1zImJslMkdgRdQ1u12lEFLft32lVVE/AZlEI0La4w7U/5wfTczOz4Smpnvm2CWFZQBN86JiLMpA4+8m4YjwBFxBPCHlFEe7wT2ojTPWlGnj7sPtInrjohYR1k3VzWUXxcR7wY+EhHPr8u/DtiJcsT9KZQ+UGTmA3U//nhEPJWynz+J0l/qdzLztU0fez5wCmXkyZbrpfaZeibw39rFLmm0WTdbN2PdbN086nIARufyNjtutBm5kPIj+CPg6w1lOwEfp1Quj1COzl4MvLjpdadRfqTH5jmn6b1PoRz9fKTO96cdxvTECJGUH9jLKD+wD1MuGXEGm46a+UpK5fpwfe0iNo5A+eoW778D5dIP91GaP/1lq1iAp1J+iO+nVILfBZY2TH9jXXePlK/rpNbfzpR+O7+gNNV5H6XiW93Bdsx6e5xSYa0EPgg8ZRL7w+6USx3cXtfb3dS+MxO8bjXw4XGmbzYCJeUPy22UPyZXUy6P8MT7AM+gXGbhLkqztDWUply7dLoPtInlk8CVbaa9nDL4yXpKX6J7KX2algLRNO+bKAOq/LLuC9cAf9biPaPGlsB+bT73tZTRN7fv92+CN2/e+n9rVffUcutm62brZuvmkb1F3TCSpEmqA4lcCyzIzLv7HQ9ARJxPGbjmrf2ORZKkmWbdrE6YBEvSNETEZcD1mfmXAxDLQkp/qN/IzFX9jkeSpH6wbtZEHBhL0pRExJMar5/XfOt3fDPoXZRmVYNgAfDHVrKSNJqsm59g3axxeSZY0pRExGeA48eZZZ/MXD0z0UiSJOtmqTMmwZKmJCIWsfGC8q3cmJmPjDNdkiR1kXWz1BmTYEmSJEnSyLBPsCRJkiRpZJgES5IkSZJGhkmwJEmSJGlkmARLkiRJkkaGSbAkSZIkaWT8f8Svns9UedUxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x576 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.suptitle('S/B selection from NN', fontsize='15')\n",
    "Ds_mass_MC =[MC_Ds_sig_dict[\"Ds_ConsD_M\"][i][0]/1000 for i in range(m_s)]\n",
    "Dplus_mass_MC =[MC_Dplus_sig_dict[\"Dplus_ConsD_M\"][i][0]/1000 for i in range(m_plus)]\n",
    "NN_selected = X_dict[i][np.argmax(output,1).astype(np.bool)]\n",
    "Ds_mass_sel_NN = NN_selected[:,dim-1]/1000\n",
    "Ds_mass_train_NN =X_dict[i][:,dim-1:dim]/1000\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.hist(Ds_mass_MC+Dplus_mass_MC,bins=70, label='MC signal events');\n",
    "plt.legend(fontsize='15')\n",
    "plt.ylabel('# events', fontsize=15)\n",
    "plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
    "plt.subplot(1,2,2)\n",
    "\n",
    "plt.hist(Ds_mass_sel_NN,alpha=0.6,bins=70, label='S selected from test set');\n",
    "plt.hist(Ds_mass_train_NN,alpha=0.2,bins=70, label='Test set (S+B)');\n",
    "plt.legend(fontsize='15')\n",
    "plt.ylabel('# events', fontsize=15)\n",
    "plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
    "fig=plt.gcf();\n",
    "fig.set_size_inches(16,8)\n",
    "plt.savefig('/home/hep/davide/D_s_NN.png', format='png', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f6790669c18>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAF9JJREFUeJzt3Xtw1OW9x/HPl5iIFDxaLtUpxagDDjSGlIaAAmqMArVgEZSLXLz1ADoUhYrGOtPjeBktbcVjByeHQWkVjkY5tTBFR9Aj3gpygkQEcxBELvHYSrFyYAyYkOf8sXFPbmR/Ifvb3Wf3/ZphyGZ/+e332d188uzze57fz5xzAgD4o1OyCwAAtA/BDQCeIbgBwDMENwB4huAGAM8Q3ADgGYIbADxDcAOAZwhuAPDMKWHstEePHi43NzeMXQNAWtq8efPfnXM9g2wbSnDn5uaqoqIijF0DQFoys71Bt2WoBAA8Q3ADgGcIbgDwTChj3AASp7a2VtXV1Tp69GiyS0EAnTt3Vu/evZWdnX3S+yC4Ac9VV1erW7duys3NlZkluxy0wTmngwcPqrq6Wueee+5J74ehEsBzR48eVffu3QltD5iZunfv3uFPRwQ3kAYIbX/E47UiuAHAM4xxA2lm0bqP4rq/eVf2a/P+gwcPqqSkRJL017/+VVlZWerZM7IAcNOmTcrJyYlrPR110003qbS0VH379tXChQtVWloqSTp+/Lguu+wyvfXWW0muMLb0Cu7XH256u/ie5NQBZJDu3bursrJSknTfffepa9euuvPOO5ts45yTc06dOiX/Q/6yZcskSXV1dXrkkUeiwZ2VleVFaEsMlQAIya5du5SXl6fZs2dr0KBB2r9/v84444zo/c8995x++tOfSpL+9re/afz48SosLFRRUZE2btzYYn9Lly7VNddco1GjRumCCy7Qgw8+GL1v4cKFysvLU15enn73u99Jkg4fPqwf/ehHGjhwoPLy8rRy5UpJ0vDhw1VZWanS0lIdPnxYBQUFmjFjhurq6qL1TZgwQWvXro3uf9q0aVq1apXq6uo0f/58FRUVKT8/X0uXLpUkffrppxo+fLgKCgqUl5env/zlL3F+NptKrx43gJTy4YcfatmyZSorK1NdXd0Jt5s7d67uuusuDR06VHv27NGYMWO0bdu2Fttt2rRJ27ZtU05OjgYPHqwxY8bo66+/1ooVK7Rp0yYdP35cRUVFuvTSS1VVVaXc3Fy9/PLLkqRDhw412dcjjzyipUuXRj8tNK5v8uTJKi8v18iRI3X06FG98cYbevLJJ7VkyRL16tVLmzZt0rFjxzR06FCNHDlSzz77rMaOHau7775bx48fV01NTTyevhMiuAGE5vzzz9fgwYNjbvfqq69qx44d0dv/+Mc/VFNTo9NOO63JdqNGjdKZZ54pSRo3bpzefvttHTt2TBMmTFCXLl2afL+4uFilpaUqLS3V2LFjNWzYsMB1//jHP9bPf/5z1dbWas2aNbr88st16qmnau3ataqqqtJzzz0nKfLHYOfOnRo8eLBmzZqlo0ePaty4cRo4cGDgxzoZBDeA0HzrW9+Kft2pUyc556K3G89lds4FOpDZfCqdmTXZZ2P9+/dXRUWFXnrpJS1YsEBjxozRL37xi0B1d+nSRcOGDdO6detUXl6um266KVrnE088ET0Y29j69eu1Zs0aTZ06Vffcc4+mTp0a6LFOBmPcABKiU6dOOvPMM7Vz507V19frxRdfjN53xRVXaPHixdHb3wxfNLd27Vp9+eWX+uqrr7Rq1SoNGzZMl1xyiV588UXV1NToyJEjWrVqlUaMGKFPP/1UXbt21fTp0zV//ny99957TfZ1yimRfuuJhnAmT56sJ598Uhs2bNAVV1whKdLjf+KJJ6I/s2PHDtXU1Gjv3r0666yzNHPmTN14443asmXLyT9RAdDjBtJMrOl7yfSrX/1Ko0ePVp8+fTRgwAAdO3ZMkrR48WLdeuutWrZsmerq6lRcXNwkyL8xfPhwXX/99fr44481ffp0FRQUSJKmTJkSHZK59dZbdeGFF+qll15SaWmpOnXqpJycHJWVlbXY3y233KL8/HwVFhbqqaeeanLf6NGjdcMNN+i6666Lnldk1qxZ2rdvX/Rxe/XqpVWrVum1117To48+quzsbHXt2lXLly+P35PWCjvRx4yOKCwsdEm5kALTAZGBqqqq1L9//2SXEbqlS5dq27Zteuyxx5JdSoe19pqZ2WbnXGGQn/erx00wA4BnwQ0gY30z5xscnAQA7xDcAOAZv4dKmo95A0AGCNzjNrMsM9tiZn8OsyAAQNva0+O+XVKVpNNDqgVAPMT7k2iA2VtZWVm68MILo7f/9Kc/KTc3t9Vt2zoXSaJVVFTo6aef1uOPP67169crJydHF198sSSprKxMXbp00YwZM5JcZUuBgtvMekv6saSHJM0PtaJ4avwGZuogEJrTTjvthKsdU1lhYaEKCyNTp9evX6+uXbtGg3v27NnJLK1NQYdKHpN0l6T6EGsBkEb27NmjESNGaNCgQRo0aFCrpzrdvn27ioqKVFBQoPz8fO3cuVOStHz58uj3Z82apePHj7f42dzcXN19990qKipSUVGRdu3aJUnau3evSkpKlJ+fr5KSEu3bt0+S9MILLygvL08DBw7UJZdcIikS1mPGjNGePXtUVlamRYsWqaCgQG+99Zbuu+8+/eY3v1FVVZWKioqatCs/P1+StHnzZl166aX64Q9/qFGjRumzzz6TJD3++OMaMGCA8vPzNXny5Dg+qxExg9vMxkj63Dm3OcZ2M82swswqDhw4ELcCAaS+mpoaFRQUqKCgQNdcc42kyHLwdevW6b333lN5ebnmzp3b4ufKysp0++23q7KyUhUVFerdu7eqqqpUXl6ud955R5WVlcrKytKKFStafdzTTz9dmzZt0pw5c3THHXdIkubMmaMZM2Zo69atmjp1avRx77//fr3yyit6//33tXr16ib7yc3N1ezZszVv3jxVVlZqxIgR0fv69++vr7/+Wrt375YklZeXa+LEiaqtrdXPfvYzrVy5Ups3b9bNN9+se++9V1LklLFbtmzR1q1bW11q31FBhkqGSbrazK6S1FnS6Wa23Dk3rfFGzrklkpZIkSXvca8UQMpqbaiktrZWc+bMiYbvRx+1vKTaRRddpIceekjV1dUaP368+vbtq9dee02bN2+OnnukpqZGvXr1avVxp0yZEv1/3rx5kqQNGzboj3/8oyRp+vTpuuuuuyRJw4YN04033qiJEydq/Pjx7WrfxIkT9fzzz6u0tFTl5eUqLy/Xjh07tG3bNl155ZWSIpc+O/vssyVJ+fn5mjp1qsaNG6dx48a167GCiBnczrl7JN0jSWZ2maQ7m4c2ADS3aNEifec739H777+v+vp6de7cucU2119/vYYMGaI1a9Zo1KhRWrp0qZxzuuGGG/Tww7EPsjY+zeuJrp7+zffLysr07rvvas2aNSooKGjXmPykSZN03XXXafz48TIz9e3bVx988IG+//3va8OGDS22X7Nmjd58802tXr1aDzzwgLZv3x49G2E8sAAHQCgOHTqks88+W506ddIzzzzT6jj17t27dd5552nu3Lm6+uqrtXXrVpWUlGjlypX6/PPPJUlffPGF9u7d2+pjlJeXR/+/6KKLJEkXX3xx9EIHK1as0PDhwyVJH3/8sYYMGaL7779fPXr00P79+5vsq1u3bjp8+HCrj3P++ecrKytLDzzwgCZNmiRJuuCCC3TgwIFocNfW1mr79u2qr6/X/v37VVxcrIULF+rLL7/UkSNH2vXcxdKuPwHOufWS1se1AgDxlSIzqG677TZNmDBBL7zwgoqLi5tcVOEb5eXlWr58ubKzs3XWWWfpl7/8pb797W/rwQcf1MiRI1VfX6/s7GwtXrxY55xzToufP3bsmIYMGaL6+no9++yzkiIHBm+++Wb9+te/Vs+ePaMXB16wYIF27twp55xKSko0cOBAvfHGG9F9jR07Vtdee61WrVoVvW5lY5MmTdKCBQv0ySefSJJycnK0cuVKzZ07V4cOHVJdXZ3uuOMO9evXT9OmTdOhQ4fknNO8efOaXGszHvw6rWtH5qemyJsZiLdMOa1rc7m5uaqoqFCPHj2SXUq7dfS0rgyVAIBn/D5XCYCMtWfPnmSXkDT0uIE0EMaQJ8IRj9eK4AY817lzZx08eJDw9oBzTgcPHmx1amR7MFQCeK53796qrq4WK5b90LlzZ/Xu3btD+yC4Ac9lZ2fr3HPPTXYZSCCGSgDAMwQ3AHiG4AYAz2TuGHfzVZisrATgCXrcAOCZzOlxc0V4AGmCHjcAeIbgBgDPENwA4BmCGwA8Q3ADgGcIbgDwDMENAJ4huAHAMwQ3AHiG4AYAzxDcAOCZzDlXCYC0tWjdR01uz7uyX5IqSQx63ADgGYIbADxDcAOAZwhuAPAMwQ0AniG4AcAzBDcAeCb153FzrUggI2Xa3Oz2oMcNAJ5J/R43ALSieY88kxDcADKaj0MyDJUAgGcIbgDwDMENAJ6JOcZtZp0lvSnp1IbtVzrn/iXswgAgGXwY8w5ycPKYpMudc0fMLFvS22b2snNuY8i1AUAofJ+REjO4nXNO0pGGm9kN/1yYRQFAR8QzmFOxBx5ojNvMssysUtLnktY5595tZZuZZlZhZhUHDhyId50AgAaB5nE7545LKjCzMyS9aGZ5zrltzbZZImmJJBUWFtIjB9Auvg9fJFK7ZpU4576UtF7S6FCqAQDEFDO4zaxnQ09bZnaapCsk/XfYhQEAWhdkqORsSX8wsyxFgv5559yfwy0LAHAiQWaVbJX0gwTUAgAIgJWTAOAZghsAPMNpXQGEJhUXr6QDetwA4Bl63ACSpj2Lblig8//ocQOAZ+hxA0gYes3xQY8bADxDcAOAZwhuAPAMwQ0AniG4AcAzBDcAeIbgBgDPENwA4BmCGwA8Q3ADgGcIbgDwDMENAJ7hJFMn8vrDTW8X35OcOgCgGXrcAOAZghsAPENwA4BnCG4A8AzBDQCeIbgBwDNMBwyq8fRApgYCSCJ63ADgGXrcAOKGq7gnBj1uAPAMwQ0AniG4AcAzBDcAeIaDk99ofjZAAEhRBHc8cApYAAnEUAkAeIbgBgDPENwA4JmYwW1m3zOz182sysy2m9ntiSgMANC6IAcn6yT93Dn3npl1k7TZzNY55z4MuTYAQCtiBrdz7jNJnzV8fdjMqiR9VxLBHRRnFgQQR+2aDmhmuZJ+IOndMIrxBnO+ASRR4OA2s66S/kPSHc65/23l/pmSZkpSnz594lYggNTS/AyA867sl6RKMlegWSVmlq1IaK9wzv2xtW2cc0ucc4XOucKePXvGs0YAQCNBZpWYpCclVTnnHg2/JABAW4L0uIdJmi7pcjOrbPh3Vch1AQBOIMiskrclWQJqAQAEwMpJAPAMZwcEgHZoPKsmWTNq6HEDgGcIbgDwDMENAJ4huAHAMxycBNAhzZfAI3z0uAHAMwQ3AHiG4AYAzxDcAOAZghsAPENwA4BnmA4YBi5tBiBEBDeAJrg0WepjqAQAPENwA4BnCG4A8AzBDQCe4eAkgDZxEqnUQ48bADxDcAOAZwhuAPAMY9xAGoq1iIZFNn4juFNN8+Xyxfckpw4AKYuhEgDwDD3uRKNHjSRgSl96occNAJ4huAHAMwQ3AHiG4AYAz3Bw0jeND25yYBPISPS4AcAzBDcAeIahEgDM8/YMPW4A8Aw97nTGKk0gLRHcydY8XAEgBoZKAMAzMXvcZvaUpDGSPnfO5YVfEoCTwQHGzBGkx/17SaNDrgMAEFDM4HbOvSnpiwTUAgAIIG4HJ81spqSZktSnT5947RbACTA0krnidnDSObfEOVfonCvs2bNnvHYLAGiGWSUA4BmCGwA8EzO4zexZSRskXWBm1WZ2S/hlAQBOJObBSefclEQUggzGOcaBdkm9Je8sAQcktZw1Mu/KfkmqJL6G7lvS5PbGPjOTVIm/Ui+4ASQFgeoPghsIWeOec7r0mpFcBDdODuPSCZfMBTf0xluXrOEsghuIs44ErK+rIQn2xGIeNwB4hh53qmOWTVqLdw+7cc83zF5v8x52mPum994SwY3WcdmztBdm+CJcBDeQQO3tYdP7RGsY4wYAz9DjziRtjZfHcygk1jAL4/ZJwdBI+iC4fZaB49AsZkmcjgQ9QzzhIriRkeK5cMLXudfwF8GdThiCSAn0NltimCa+CG5kjkZ/2IbuO9h2oMb6I5gBw1KJQqi3H8F9Aht2H2xy+6LzuiepkgRJx956szY1f03bEvP1b/ZHIFHozUMiuIF2a88fACAMBDeCSdUeOWcpTHttLePP1E8gBDe80uQX9fUTD1+kkvaGC2O+iIXgBgJgeASphOBGUrU4CFjc9P7mc6SHtmdfKXpAOZ496rb2Rc89fWVUcPvyi53u6L0CHZNRwY0k6cDYc1u9xnj+AWhPzx7JE+tTRKYcrOTsgADgGXrcGcSLoaIEzgzJlN4Z/l+irhAUNoLbM/EM38b7SpUQj+fwB2PpSFcEdwg6Eq7t/dlUCN8WZ9rjXQWEil8xoAHT5zKbT0NnBHcKaOsjfaI+7sd6nLZ68y0CL0WGXYB05X1wtxU48Rw68OHAng81Ar5KpQOb3gd3snSkJ5zuB83SvX1AshHcCZAOQUZvHj5K11MCsAAHADyT1j3udOjpdkSY7c/05xbpL5V75N4FNws0AGS6lA9uwhUAmmKMGwA8Q3ADgGcCDZWY2WhJ/yopS9JS59wjoVYFACks2cvjY/a4zSxL0mJJP5I0QNIUMxsQdmEAgNYFGSopkrTLObfbOfe1pOck/STcsgAAJxIkuL8raX+j29UN3wMAJEGQMW5r5XuuxUZmMyV9M9BzxMx2tKOOHpL+3o7t0wXtziyZ2m4p7dv+W0nS/JZ3tKfd5wR9tCDBXS3pe41u95b0P803cs4tkXRSS43MrMI5V3gyP+sz2p1ZMrXdUua2Pax2Bxkq+S9Jfc3sXDPLkTRZ0up4FwIACCZmj9s5V2dmcyS9osh0wKecc9tDrwwA0KpA87idcy9JeinEOlL3bC7hot2ZJVPbLWVu20NptznX4jgjACCFseQdADyTsOA2s9FmtsPMdplZaSv3n2pm5Q33v2tmuYmqLWwB2j7fzD40s61m9pqZBZ4WlMpitbvRdteamTOztJh1EKTdZjax4TXfbmb/nugawxDgfd7HzF43sy0N7/WrklFnvJnZU2b2uZltO8H9ZmaPNzwvW81sUIcf1DkX+j9FDmp+LOk8STmS3pc0oNk2t0kqa/h6sqTyRNSWIm0vltSl4etb06HtQdrdsF03SW9K2iipMNl1J+j17itpi6QzG273SnbdCWr3Ekm3Nnw9QNKeZNcdp7ZfImmQpG0nuP8qSS8rsiZmqKR3O/qYiepxB1k2/xNJf2j4eqWkEjNrbfGPb2K23Tn3unPuq4abGxWZK++7oKdKeEDSQklHE1lciIK0+58lLXbO/UOSnHOfJ7jGMARpt5N0esPX/6RW1oP4yDn3pqQv2tjkJ5KedhEbJZ1hZmd35DETFdxBls1Ht3HO1Uk6JCkdrkjb3lMG3KLIX2ffxWy3mf1A0vecc39OZGEhC/J695PUz8zeMbONDWff9F2Qdt8naZqZVSsyS+1niSkt6eJ+2pBEXQEnyLL5QEvrPRS4XWY2TVKhpEtDrSgx2my3mXWStEjSjYkqKEGCvN6nKDJccpkin67eMrM859yXIdcWpiDtniLp986535rZRZKeaWh3ffjlJVXcsy1RPe4gy+aj25jZKYp8lGrr44cvAp0ywMyukHSvpKudc8cSVFuYYrW7m6Q8SevNbI8iY3+r0+AAZdD3+irnXK1z7hNJOxQJcp8Fafctkp6XJOfcBkmdFTmXR7oLlAHtkajgDrJsfrWkGxq+vlbSf7qGkX3PxWx7w5DBvykS2ukw3inFaLdz7pBzrodzLtc5l6vI2P7VzrmK5JQbN0He639S5IC0zKyHIkMnuxNaZfwFafc+SSWSZGb9FQnuAwmtMjlWS5rRMLtkqKRDzrnPOrTHBB55vUrSR4oceb634Xv3K/LLKkVexBck7ZK0SdJ5yT5anMC2vyrpb5IqG/6tTnbNiWh3s23XKw1mlQR8vU3So5I+lPSBpMnJrjlB7R4g6R1FZpxUShqZ7Jrj1O5nJX0mqVaR3vUtkmZLmt3o9V7c8Lx8EI/3OSsnAcAzrJwEAM8Q3ADgGYIbADxDcAOAZwhuAPAMwQ0AniG4AcAzBDcAeOb/AA8fpFoRigoNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "true_positives=output[:,1][np.where(Y_test==1)[0]]\n",
    "false_positives=output[:,1][np.where(Y_test==0)[0]]\n",
    "plt.hist(true_positives,alpha=0.5,bins=80,density=True,label=\"True positives\");\n",
    "plt.hist(false_positives,alpha=0.5,bins=80,density=True, label=\"False positives\");\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "pAUC 0.8704610495805598\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAHVCAYAAADLvzPyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmczfXix/HXx9iyryF7ISQpE4UyxtjDRbKErlRuJYXW23ZvdbXThqsQUciS3Mk2lpnpp81UchEhxWTft7HMzOf3x5ebNMwZc2Y+Z3k/H4/zOOfMfJvzdpqZ93w+3+/38zXWWkRERCRw5HEdQERERP5I5SwiIhJgVM4iIiIBRuUsIiISYFTOIiIiAUblLCIiEmBUziIiIgFG5SwiIhJgVM4iIiIBJq+rFy5TpoytVq2aq5cXERHJVd9+++0ea21ZX7Z1Vs7VqlUjKSnJ1cuLiIjkKmPMr75uq2ltERGRAKNyFhERCTAqZxERkQCjchYREQkwKmcREZEAo3IWEREJMCpnERGRAKNyFhERCTAqZxERkQCjchYREQkwKmcREZEAo3IWEREJMCpnERGRAKNyFhERCTCZlrMxZoIxZpcxZvV5Pm+MMW8ZYzYaY1YZY67zf0wREZHw4cvIeSLQ9gKfbwfUPH27BxiT/VgiIiLhK9NyttYmAvsusEln4APr+QooYYyp4K+AIiIiuSktDY4ehf373WXI64evURHYetbz5NMf2+6Hry0iIgLAiRNw6BAcPpzx7ehROH7cu6Wk/P743OcX+tzx43DqlPd6JUqksX9/hJN/qz/K2WTwMZvhhsbcgzf1TZUqVfzw0iIiEqjS0uDIkT8W6IXKNbPbmdL0Rb58ULCgd7vkkj8/Llnyz58783zx4lhWrEikY8e/AE1y7P25EH+UczJQ+aznlYBtGW1orX0XeBcgMjIywwIXERE3rPVGkZmVpK8Fe+yYb6+bJw8ULfr7rVgx7758+T9+PLNb4cK/l21ENga8AwdeTULCPvr1c1PM4J9yngsMMsZMAxoDB621mtIWEclF1nqFuHevd9u37/f7gwd9K9gjR7zRri8KFfpzoVasmLUyPXO75BIwGc3B5qKjR48yZswYhg4dStWqVenXr5/TPJmWszFmKhAFlDHGJAPPAvkArLX/BuYB7YGNwDGgf06FFREJBydO/Llkz76d+7F9+7zbhaZ98+X7cymWKAFVqmS9TIsUyd7INNAcPnyYDh06sHz5cpo2bcqNN97oOlLm5Wyt7ZXJ5y1wv98SiYiEiLQ0OHAgayW7d693YNP5FCwIpUtDqVLefd263v3ZHzv7VqoUFC8OBQrk3r87mBw4cIB27dqxYsUKPvroo4AoZvDPtLaISEiz1ivMzEr13I/t3+/9txnJk8crzjOFWqkSXHPN+Qv2zONChXL33x7K9u3bR5s2bfjhhx+YMWMGXbp0cR3pf1TOIhKWrIU9eyA5+Y+3XbsyLt6TJ8//tYoU+WOZVq+eeckWL+4VtLjz448/snHjRmbPns0tt9ziOs4fqJxFJOSkp8Pu3b8X7tatfy7h5GRv3+7ZIiKgbNnfC7RmzQsX7Jnn+fO7+XfKxTlx4gQFChSgadOmbN68mRIlSriO9CcqZxEJKmlp3uj2fKW7dSv89tufD47Kl887mrhSJWjUCLp29R6ffStXLrQOdJI/27ZtGzExMTzyyCP0798/IIsZVM4iEkDS0mDHjvOXbnIybNsGqal//O8KFPi9YJs2/f1x5cq/Py5bVtPI4W7r1q1ER0ezY8cOatSo4TrOBamcRSRXWOvtu/31199vW7f+sYi3b//zebaXXPJ7wUZF/Xm0W6kSlCnj/jxZCWybN28mOjqaffv2sWjRooA5Kvt8VM4i4hepqd508pYtfyzgM8+3bPnzilGFC/8+uo2J+WPhnvl4yZIqXsmegwcP0rx5c44cOcKSJUuIjIx0HSlTKmcR8cnRo38s2nML+Lff/jzqLVsWqlaFq66C9u29BS+qVv39puKV3FC8eHGGDh1KVFQUDRo0cB3HJypnEfnfaUVnl+65Bbx37x//m7x5vZFtlSrQvPnvhXumgKtU8aakRVxZs2YNx44d4/rrr+ehhx5yHSdLVM4iYWbvXli9Gtas8e7PPN53zlXbCxf+vXCvv/7P5XvZZTqyWQLXDz/8QExMDJdeeimrVq0iIsi+WVXOIiHq0CFYu/b3Aj5Twjt2/L5N8eJQrx7ceivUqQPVqv1ewKVKacpZglNSUhKtW7emcOHCzJkzJ+iKGVTOIkEvJQV+/PGPBbx6tTctfUahQt5+37ZtvTI+c7vsMhWwhJavvvqKNm3aUKpUKZYuXUr16tVdR7ooKmeRIJGe7pXwf//7xynpTZt+X785f36oXRuaNfu9gK+6yhsR6xxfCQejRo2ibNmyLF26lCpVqriOc9FUziIBbNs2iIuDRYu8+927vY9HRHhLSzZoAH36eAVcrx7UqOEdqCUSbqy1GGMYN24cBw4coFy5cq4jZYt+jEUCSEoKfP65V8YLF3ojY4BLL4U2baBlS7juOrjySl0CUOSMRYsW8fTTT/PZZ59RpkyZoC9mUDmLOGWtN029aJF3S0z0LsaQPz/cdBP07euV8tVXa1paJCOfffYZXbt2pU6dOtjzXZ8zCKmcRXLZzp1/nKo+c/T0VVfBffdB69Zw8826bq9IZubMmcNtt91G/fr1WbRoEaVKlXIdyW9UziK54NAhmDEDPvjAGx2Dtx50q1ZeGbdq5V0xSUR8ExsbS/fu3YmMjGT+/PkBe3Wpi6VyFskhaWneyPiDD+CTT+D4cahVC557Djp08A7m0lS1yMW57rrr6N27N2+//TbFihVzHcfvVM4ifrZ6NUyaBB9+6F1lqWRJuPNO6NfPu46wzisWuXjLli3jpptu4rLLLmPSpEmu4+QY/d0u4ge7dsEbb3hHUl99tfe4USOYNcsr6FGjoHFjFbNIdrz77rtER0fzxhtvuI6S4zRyFrlIhw/DnDkwbZp32lNaGkRGwltvQc+e3hWZRMQ/Ro0axaBBg2jfvj2DBg1yHSfHqZxFsiAlBebNg6lT4bPPvP3IVarAww9709Z167pOKBJ6RowYwbBhw+jcuTPTp0+nQBic5K9yFsnEqVPegV3Tpnkj5cOHoVw5uPtu6NULbrhB09UiOeW3337jmWeeoXv37nz44Yfky5fPdaRcoXIWyUBamnfK07RpMHOmdznFkiWhRw9vyjoqSpdLFMkNFStW5Msvv6ROnTrkDaO1acPnXyrigxMn4P334aWX4NdfvWsad+7sjZBbt/ZW7hKRnGWt5amnnqJChQoMGjSIq6++2nWkXKejtUXw9iW//TZccQXce693KcVp07yjsD/8EG65RcUskhustTz88MMMHz6cNWvWhNSSnFmhkbOEtaNHYexYePVVbxnNm27yzlGOjtZ+ZJHcZq1l8ODBvPPOOzzwwAO8+eabmDD9QVQ5S1g6fBjGjIHXXvMuwxgd7Y2Umzd3nUwkPFlruffeexk7dizDhg3j1VdfDdtiBpWzhJndu72R8siR3kFebdrA009D06auk4mEN2MMV155JU888QT/+te/wrqYQeUsYcBa+Pprb5Wujz+Gkye9ta2fftpbtUtE3ElNTeWnn36ibt26DBkyxHWcgKEDwiRkHTsG48dDw4Zw443w6afeuclr1kBsrIpZxLVTp07Ru3dvbrjhBrZv3+46TkDRyFlCzsaNMHq0d0rUgQPedZJHj4Y+faBoUdfpRATgxIkT9OjRg08//ZTXX3+dChUquI4UUFTOEjJ++gleeME79SlPHujWDe67zzsCO8x3X4kElOPHj9OtWzfmzZvH22+/HRZrZWeVylmC3tmlXKAADBkCw4aB/hAXCUxvvfUW8+fPZ+zYsdxzzz2u4wQklbMErYxK+ZFHvHWvRSRwDRkyhIYNG9KyZUvXUQKWDgiToLNhA9xxB9Sp4617PWQIbN7snbOsYhYJTIcOHaJ///7s3LmTfPnyqZgzoXKWoJGW5k1X164NM2aolEWCxYEDB2jdujVTpkwhKSnJdZygoGltCQrHj8Ptt8Ps2TBwIPzznypkkWCwb98+WrduzapVq5g5cyYdOnRwHSkoqJwl4B044F0ZKjHRW9nroYdcJxIRX+zevZtWrVqxbt065syZQ/v27V1HChoqZwloycnQrh2sXw9Tp3rXUhaR4GCtJW/evPznP/+hVatWruMEFZWzBKy1a6FtW2/kvGCBd3EKEQl8O3fupFSpUlx66aV888035Mmjw5uySu+YBKTly6FZMzh1ypvOVjGLBIctW7bQtGlTBg4cCKBivkh61yTgzJ0LMTFQpgx88QU0aOA6kYj4YvPmzTRv3pw9e/b8r5zl4qicJaC89x506QL163uj5+rVXScSEV9s2LCBm2++mUOHDrFkyRIa68oy2aJyloBgrXd61D33ePuZly6FsmVdpxIRX6SlpdG5c2eOHz/O0qVLadiwoetIQU8HhIlzqalw//3w7rvw17969/nyuU4lIr6KiIhgwoQJFC1alKuuusp1nJCgkbM4lZICt97qFfKTT8KECSpmkWCxcuVK3nnnHQBuuOEGFbMfaeQszuzbBx07wpdfwttvg64aJxI8kpKSaN26NUWKFKFfv34UK1bMdaSQonIWJ7Zs8fYtb9oEH3/sjZ5FJDh8+eWXtG3bllKlSrFs2TIVcw5QOUuuW7fOO1Xq8GFYuBCiolwnEhFfJSYm0qFDB8qXL8/SpUupXLmy60ghSeUsuWr1ajhzpbjPP/dOmRKR4LF+/XoqVarEkiVLuOyyy1zHCVnGWuvkhSMjI60uHRZevv8eWrWCAgVgyRLv0o8iEhwOHTr0v+nr48ePU7BgQceJgo8x5ltrbaQv2+pobckV33zjLcFZuLC3HKeKWSR4xMbGUq1aNb766isAFXMuUDlLjlu+3NvHXKqUV8xXXOE6kYj46pNPPqFr165cccUV1KpVy3WcsKFylhy1bBm0aQMVKkBCAlSt6jqRiPjq448/pnv37jRs2JDFixdTqlQp15HChspZcszChdC+vVfICQlQqZLrRCLiq+XLl9OrVy+aNGnCokWLKF68uOtIYUXlLDli4ULo1AmuvBLi46F8edeJRCQrbrjhBoYPH878+fMpWrSo6zhhR+UsfpeYCH/5C9StqwtYiASbyZMnk5ycTEREBI899hiFCxd2HSksqZzFr1asgFtugWrVYNEi7yAwEQkOb731Fv369ePll192HSXsqZzFb1av9pbkLFMGFi/WiFkkmLz66qs8+OCDdOnShddff911nLCncha/2LDBO12qYEFvgZGKFV0nEhFf/etf/+LRRx+lR48eTJ8+nfz587uOFPZUzpJtW7Z4xZyW5o2Yq1d3nUhEfJWSksLHH39M3759mTJlCvl0zdaAoLW1JVt27PCK+eBB75zmOnVcJxIRX1hrSU1N5ZJLLiEhIYGiRYsSERHhOpacppGzXLQ9e6B1a9i2DebNg2uvdZ1IRHxhrWXYsGF069aNU6dOUaJECRVzgFE5y0XZt8+7iMWGDfDpp9CkietEIuKL9PR0HnjgAUaOHEn16tXJm1cTqIFI/1ckyw4e9JbkXLvWK+Yzl4AUkcCWnp7OwIEDGTduHA8//DCvvPIKxhjXsSQDGjlLlhw+DO3awQ8/wKxZ3qlTIhIchgwZwrhx43jyySdVzAFOI2fx2dGj0KGDd/nHGTO8xUZEJHj069ePihUr8uijj7qOIpnQyFl8cuwYdOzoXf7xo4+gSxfXiUTEF6dOnWLmzJkANGzYUMUcJFTOkqnjx70yjo+HSZPgtttcJxIRX5w4cYJbb72V7t27k5SU5DqOZIFP5WyMaWuMWW+M2WiMeTyDz1cxxiwzxnxvjFlljGnv/6jiwokT0K2bt072+PHQp4/rRCLii5SUFLp06cLcuXMZNWoUkZGRriNJFmRazsaYCGAU0A6oC/QyxtQ9Z7OngI+ttdcCPYHR/g4quS81FXr29M5hHjsW+vd3nUhEfHHs2DE6derEggULeO+997jvvvtcR5Is8mXk3AjYaK392Vp7EpgGdD5nGwsUO/24OLDNfxHFlZEjYc4cePNNuOce12lExFfx8fHEx8fz/vvvc9ddd7mOIxfBl6O1KwJbz3qeDDQ+Z5t/AIuMMQ8AhYGYjL6QMeYe4B6AKlWqZDWr5KING+CZZ7zrMj/wgOs0IuILay3GGNq3b8/69eu5/PLLXUeSi+TLyDmjE+HsOc97AROttZWA9sBkY8yfvra19l1rbaS1NrKsricYsNLT4e67oUABGDUKdCqkSODbv38/UVFRxMXFAaiYg5wvI+dkoPJZzyvx52nrAUBbAGvtl8aYgkAZYJc/Qkrueu89SEiAcePgsstcpxGRzOzdu5dWrVqxevVqUlJSXMcRP/Bl5LwCqGmMqW6MyY93wNfcc7bZArQEMMbUAQoCu/0ZVHJHcjI88ghER8Odd7pOIyKZ2bVrFy1atGDt2rV8+umndOrUyXUk8YNMR87W2lRjzCBgIRABTLDWrjHGPAckWWvnAsOA94wxQ/CmvP9qrT136lsCnLXwt795R2m/956ms0UC3Zmp7F9++YXY2FhiYjI83EeCkE/Ld1pr5wHzzvnYM2c9Xgs09W80yW3TpsFnn8GIEaDdVSKBr3jx4sTExNCtWzeaN2/uOo74kXE1wI2MjLRasSZw7NkDdep4pfzFF6BLu4oErl9//RVrLdWqVXMdRbLAGPOttdan1WB04QsB4KGHvEtBjh+vYhYJZD///DMtWrSgdOnSfPvtt7qyVIhSOQuffQYffgjPPgv16rlOIyLn89NPPxEdHU1KSgqffPKJijmEqZzD3N693upfV10FTzzhOo2InM/atWtp2bIlaWlpLFu2jPr167uOJDlI5RzGrIWBA2H3bvjPf7xFR0QkMD366KNYa4mPj6du3XMvbyChRuUcxiZOhFmz4KWX4LrrXKcRkQuZPHkye/bsoWbNmq6jSC7Q9ZzD1KZNMHgwNG8ODz/sOo2IZOSbb76hR48eHD9+nJIlS6qYw4hGzmEoNdW7LnPevDB5so7OFglEX3zxBW3btqVMmTLs3buXihUruo4kuUgj5zD0wgvw1Vfw739D5cqZby8iuSshIYHWrVtTvnx5EhMTVcxhSOUcZr78Ep5/Hvr2hR49XKcRkXMtXbqUdu3aUaVKFRISEqhUqZLrSOKAyjmMHD7sTWdXqQLvvOM6jYhkpEyZMjRu3Jj4+HgqVKjgOo44on3OYWTwYPjlF0hMhGLFXKcRkbP9+OOP1K5dm/r167N06VItMBLmNHIOEzNmeKdO/f3v0FSXKBEJKLNmzaJ+/fqMGzcOQMUsKudwkJzsLTbSqBE880zm24tI7pk2bRo9evSgUaNG3Hbbba7jSIBQOYc4a+Guu+DkSZgyBfLlc51IRM6YPHkyt99+O02bNmXBggUUL17cdSQJENrnHOI+/BAWLoS33gKtXyASOH799VfuvPNOoqKimDt3LoULF3YdSQKIyjmE7d7tXQryhhvgvvtcpxGRs1WtWpV58+bRrFkzLrnkEtdxJMBoWjuEDRkChw7BuHFaBUwkULz99tvExsYC0KpVKxWzZEjlHKLmz/emtJ94wrscpIi49+qrrzJ48GCmTp3qOooEOJVzCDpyBP72N6hTxzt1SkTce+GFF3j00Ufp2bMnkyZNch1HApz2OYegp56CrVvh//5P12gWcc1ay7PPPsvzzz9P3759ef/994nQfibJhEbOIebrr70js++7D5o0cZ1GRAD27dvHgAEDVMziM2OtdfLCkZGRNikpyclrh6qTJ6FhQzhwANas0RKdIi5Za9m1axflypUjPT0dgDx5NB4KZ8aYb621kb5sq++UEPLKK7B6NYwZo2IWcSk9PZ3777+fhg0bsnv3bvLkyaNilizRd0uIWLfOuxRkjx5wyy2u04iEr7S0NO655x7GjBlD7969KVOmjOtIEoRUziEgPR3uvhsKF4Y333SdRiR8paam0r9/f8aPH8/TTz/Nyy+/rItYyEXR0doh4L33vCOz338fypVznUYkfL300ktMnjyZ559/nqeeesp1HAliKucg99tv8Oij0LIl3HGH6zQi4W3w4MFUrVqVvn37uo4iQU7T2kHulVcgJQXGjgXNnonkvhMnTvDMM89w9OhRihUrpmIWv1A5B7GDB2HCBOjZE664wnUakfCTkpJC586def7551m8eLHrOBJCNK0dxMaN85bqfOgh10lEws/Ro0fp1KkTy5YtY9y4cXTu3Nl1JAkhKucglZrqrQR2881w3XWu04iEl8OHD9OhQweWL1/OpEmTNJUtfqdyDlKffAJbtujUKREXdu7cyc8//8xHH31Ejx49XMeREKRyDlJvvAGXXw4dO7pOIhI+jh49SqFChahRowY//fQThQoVch1JQpQOCAtC33wDX3wBDz4IWkNfJHfs2bOHZs2a/e/8ZRWz5CSVcxAaOdJbO7t/f9dJRMLDrl27iI6OZt26ddx0002u40gYUDkHma1bYcYMuOsuKFrUdRqR0Ld9+3aioqLYuHEjsbGxtG3b1nUkCQPa5xxk3nkHrIUHHnCdRCT0nTp1ipiYGLZs2cL8+fNp3ry560gSJlTOQeTIEXj3XejaFapVc51GJPTly5ePf/7zn1x22WU0adLEdRwJIyrnIDJpEhw4AEOGuE4iEto2bdrEmjVr6NSpE7feeqvrOBKGVM5BIj3dO6e5USO48UbXaURC1/r164mOjiY9PZ2WLVtSuHBh15EkDKmcg8S8ebBhA3z0kS5wIZJT1q5dS3R0NNZaFi9erGIWZ3S0dpAYORIqVQLNsInkjFWrVhEVFYUxhvj4eK6++mrXkSSMqZyDwA8/wNKlMGgQ5MvnOo1IaJo9ezb58+cnISGBOnXquI4jYc5Ya528cGRkpE1KSnLy2sGmf3/4+GNIToaSJV2nEQktqamp5M2bF2ste/bsoWzZsq4jSYgyxnxrrY30ZVuNnAPcjh3efua//lXFLOJvy5cvp27duqxfvx5jjIpZAobKOcCNGQMnT3rraIuI/8THx9OmTRuMMRQpUsR1HJE/UDkHsEOHYNQo6NABatVynUYkdCxevJj27dtTtWpV4uPjqVixoutIIn+gcg5gr70Ge/fCs8+6TiISOr744gtuueUWatSowbJly6hQoYLrSCJ/onIOUDt2wOuvw223wfXXu04jEjquueYa7rzzTpYtW8all17qOo5IhlTOAeq557x9zS+84DqJSGiIi4vj8OHDFC5cmNGjR1O6dGnXkUTOS+UcgDZsgPfeg7vvhpo1XacRCX5Tp06lXbt2PP30066jiPhE5RyAnnoK8ueHZ55xnUQk+E2aNIk+ffrQrFkzXtBUlAQJlXOAWbHCW3Bk2DAoX951GpHgNm7cOPr37090dDTz5s3TKVMSNLRCWACxFmJiYNUq2LQJihVznUgkeB05coQ6depw9dVXM3v2bAoWLOg6koS5rKwQpqtSBZC4OG8N7TfeUDGLZIe1liJFivD5559ToUIFChQo4DqSSJZoWjtApKfDY49BtWrwt7+5TiMSvF566SUefPBBrLVUq1ZNxSxBSeUcIKZNg5UrvVOn9LtEJOustTz33HM88cQT7Nmzh7S0NNeRRC6ayjkAnDzpHaF9zTXQq5frNCLBx1rLU089xbPPPssdd9zB5MmTyZtXe+0keOm7NwCMHQubN8P8+ZBHfy6JZNmTTz7Jiy++yN13382///1v8ugHSYKcvoMdO3TIWw2sRQto08Z1GpHgdOONN/Lggw+qmCVk6LvYsddfhz174OWXwRjXaUSCR3p6Ol9//TUAHTt25I033lAxS8jQd7JDO3d65dy9uy5uIZIVaWlp3HXXXTRp0oRVq1a5jiPid9rn7NDzz8Px4/Cvf7lOIhI8UlNT6d+/P1OmTOHZZ5/l6quvdh1JxO9Uzo5s3OgdCKaLW4j47tSpU/Tt25fp06fzwgsv8OSTT7qOJJIjVM6O6OIWIlk3c+ZMpk+fzquvvsrDDz/sOo5IjlE5O/DddzB9ulfQFSq4TiMSPHr27EmVKlVo2rSp6ygiOUoHhDnw1ltQtCjoD3+RzKWkpNC7d29Wr16NMUbFLGFB5ZzLDh2CGTOgZ08oXtx1GpHAdvToUTp06MC0adP44YcfXMcRyTWa1s5l06bBsWMwYIDrJCKB7fDhw3To0IHly5fzwQcfcPvtt7uOJJJrfBo5G2PaGmPWG2M2GmMeP882txlj1hpj1hhjPvJvzNAxfjxcdRU0auQ6iUjgOnjwIG3atOGLL77go48+ok+fPq4jieSqTEfOxpgIYBTQCkgGVhhj5lpr1561TU3gCaCptXa/MebSnAoczP77X/jmGxgxQquBiVxI/vz5KVmyJDNmzKBLly6u44jkOl+mtRsBG621PwMYY6YBnYG1Z21zNzDKWrsfwFq7y99BQ8H48ZAvH/Tt6zqJSGDas2cPERERlCxZktjYWIz+ipUw5cu0dkVg61nPk09/7Gy1gFrGmOXGmK+MMW0z+kLGmHuMMUnGmKTdu3dfXOIgdeIETJ4Mf/kLlCnjOo1I4Nm5cydRUVF06dIFa62KWcKaLyPnjH5CbAZfpyYQBVQCPjfG1LPWHvjDf2Ttu8C7AJGRked+jZD26aewb58OBBPJyLZt22jZsiVbtmzhrbfeUjFL2PNl5JwMVD7reSVgWwbbfGqtPWWt3QysxytrOW38eKhcGWJiXCcRCSxbt26lefPmJCcns2DBAqKjo11HEnHOl3JeAdQ0xlQ3xuQHegJzz9lmDtACwBhTBm+a+2d/Bg1mv/4KcXHQvz9ERLhOIxJY7rjjDnbt2sWiRYu46aabXMcRCQiZTmtba1ONMYOAhUAEMMFau8YY8xyQZK2de/pzrY0xa4E04BFr7d6cDB5MJk707vv3dxpDJCCNGzeOffv2ERkZ6TqKSMAw1rrZ9RsZGWmTkpKcvHZuSk+Hyy/3rjwVF+c6jUhgWLduHRMnTmT48OHkyaOFCiU8GGO+tdb69Feofipy2JIl3rS2DgQT8axevZqoqCg+x9TqAAAgAElEQVTef/99fvvtN9dxRAKSyjmHjR8PJUt6p1CJhLsffviBFi1akCdPHhISEqhcuXLm/5FIGFI556C9e+GTT6BPHyhY0HUaEbe+/fZbWrRoQcGCBUlISKB27dquI4kELJVzDvrwQzh5UlPaIgB79+6lXLlyJCYmUrOmzrQUuRAdEJZDrIUGDbzlOkP4nymSqT179lDm9LJ4qamp5M2ri+FJeNIBYQHg229h1Sq46y7XSUTciY+P5/LLL2fWrFkAKmYRH6mcc8j48XDJJdCrl+skIm7ExcXRvn17KleuTNOmTV3HEQkqKucccOwYfPQR3HorFC/uOo1I7ps3bx4dO3akVq1axMfHU758edeRRIKKyjkHzJwJhw7pQDAJTxs2bOAvf/kL9erVY+nSpZQtW9Z1JJGgo3LOAePHQ40acPPNrpOI5L6aNWsyZswYFi9eTKlSpVzHEQlKKmc/27ABEhPhzjtBV72TcDJt2jS+++47AAYMGECJEiUcJxIJXipnP5swAfLkgTvucJ1EJPdMnDiR3r17M3z4cNdRREKCytmPUlNh0iRo3x4uu8x1GpHc8e6779K/f39iYmL44IMPXMcRCQkqZz+aPx+2b9eBYBI+3nnnHQYOHEj79u2ZO3cuhQoVch1JJCSonP1o/HgoVw46dHCdRCTnpaens3DhQv7yl7/wySefUFALyIv4jZbr8ZMjR2DePHjgAW/JTpFQlpKSwiWXXMKMGTOIiIggn77pRfxKI2c/SUiAU6egXTvXSURyjrWWf/zjHzRp0oSDBw9SsGBBFbNIDlA5+0lcnHdZyGbNXCcRyRnWWp588kn++c9/0qBBA4oUKeI6kkjI0rS2n8TFwU036brNEpqstTz88MOMGDGCgQMHMnr0aPLk0d/2IjlFP11+8NtvsHYttGrlOolIzhg+fDgjRozggQceYMyYMSpmkRymkbMfLFni3cfEuM0hklP69etHREQEjz32GEZL34nkOP356wdxcVC2LFxzjeskIv6TlpbGuHHjSEtLo3Llyjz++OMqZpFconLOJmth8WJo2dJbtlMkFKSmptKvXz/uvvtuYmNjXccRCTua1s6m1athxw7tb5bQcerUKW6//XZmzJjB8OHD6dy5s+tIImFH5ZxNcXHevcpZQsGJEyfo2bMnc+bM4fXXX2fo0KGuI4mEJZVzNi1eDFdeCZUru04ikn1r165l0aJFvP322wwaNMh1HJGwpXLOhhMnvJXB7rzTdRKR7ElLSyMiIoJrr72WjRs3UqFCBdeRRMKaDmHKhi+/hGPHdAqVBLcjR47QqlUrxo4dC6BiFgkAKudsiIuDiAiIinKdROTiHDp0iLZt25KQkEDRokVdxxGR0zStnQ1xcdC4MRQv7jqJSNYdOHCAtm3b8u233zJt2jS6d+/uOpKInKaR80Xavx+SknSUtgSnEydOEBMTw3fffcfMmTNVzCIBRuV8kZYu9RYgUTlLMCpQoAC9evVizpw5Oo9ZJABpWvsixcVB0aLQqJHrJCK+27FjB1u3buX6669n2LBhruOIyHmonC9SXJx3IJiuMy/B4rfffiM6OpqjR4+yceNGCur6piIBS9PaF+Hnn72bprQlWGzZsoXmzZuzfft2pk+frmIWCXAaOV+ExYu9e5WzBIPNmzfTokULDhw4QFxcHI0bN3YdSUQyoXK+CHFxUKmSt2ynSKB77bXXOHz4MEuWLKFhw4au44iIDzStnUVpabBkiTdq1qVtJRiMHDmSr776SsUsEkRUzln03XfeOc6a0pZAtnr1alq2bMnu3bvJnz8/NWvWdB1JRLJA09pZdOYSkS1bus0hcj4rV64kJiaGAgUKsH//fsqWLes6kohkkUbOWbR4MVxzDVx6qeskIn+WlJREdHQ0hQoVIiEhgVq1armOJCIXQeWcBceOwfLlmtKWwJSUlETLli0pXrw4iYmJ1KhRw3UkEblIKucsSEyEkydVzhKYKlWqxE033URiYiLVqlVzHUdEskHlnAVxcVCgANx0k+skIr9buXIlqamplC9fntjYWCpXruw6kohkk8o5C+LioFkzuOQS10lEPAsXLuTGG2/k2WefdR1FRPxI5eyjHTvgv//VlLYEjtjYWDp16kTt2rUZMmSI6zgi4kcqZx8tWeLdx8S4zSEC8Mknn9C1a1fq16/PkiVLKFOmjOtIIuJHKmcfxcVB6dJw7bWuk0i4O3jwIAMGDKBhw4YsXryYUqVKuY4kIn6mRUh8YK1Xzi1bQh79OSOOFS9enEWLFnHllVdStGhR13FEJAeoanzw44+wbZv2N4tbEyZM4M033wQgMjJSxSwSwlTOPtAlIsW1f//73wwYMID58+eTlpbmOo6I5DCVsw/i4qBmTaha1XUSCUdvvfUW9957Lx06dGDOnDlERES4jiQiOUzlnIlTpyA+XqNmceO1117jwQcfpEuXLsyePZuCBQu6jiQiuUDlnIlvv4UjRyA62nUSCUcFChSgR48eTJ8+nfz587uOIyK5ROWciaQk7/6GG9zmkPBhreWXX34B4IEHHmDq1Knky5fPbSgRyVUq50wkJUH58nDZZa6TSDiw1vLEE09Qr1491q9fD4AxxnEqEcltKudMrFgBkZGg34+S06y1DBs2jJdffpm+fftSs2ZN15FExBGV8wUcOeKd43z99a6TSKhLT0/ngQceYOTIkQwePJjRo0eTRyveiIQt/fRfwPffe6uDRUa6TiKhbuLEiYwaNYqHH36YN954Q1PZImFOy3dewIoV3n3Dhm5zSOjr168fhQoVokePHipmEdHI+UKSkqByZShXznUSCUWpqak8+uijbN++nbx589KzZ08Vs4gAKucLSkrS/mbJGSdPnqRnz568+uqrfPbZZ67jiEiAUTmfx4EDsGGD9jeL/504cYLu3bsza9YsRowYwV133eU6kogEGO1zPo9vv/XuVc7iTykpKXTr1o358+czatQo7rvvPteRRCQAqZzP48zKYCpn8adjx46RnJzMe++9pxGziJyXyvk8kpLgiiugZEnXSSQUHDlyhPz581O6dGmSkpK0TraIXJD2OZ/HmZXBRLLr4MGDtGnThj59+mCtVTGLSKZUzhnYvRt+/VVHakv27d+/n9atW/PNN99w22236VQpEfGJprUzoIPBxB/27t1Lq1atWL16NbNmzaJTp06uI4lIkFA5Z2DFCu9CF9de6zqJBCtrLV27dmXt2rV8+umntGvXznUkEQkiPk1rG2PaGmPWG2M2GmMev8B2txpjrDEmqMecSUlw5ZVQrJjrJBKsjDG88sorxMbGqphFJMsyLWdjTAQwCmgH1AV6GWPqZrBdUWAw8LW/Q+Y2rQwmF+u3337j3XffBaBx48bExMQ4TiQiwciXkXMjYKO19mdr7UlgGtA5g+2eB14BjvsxX67bts27aX+zZNWWLVto3rw5Dz/8MNu3b3cdR0SCmC/lXBHYetbz5NMf+x9jzLVAZWtt7IW+kDHmHmNMkjEmaffu3VkOmxu0+IhcjJ9//pmbb76ZPXv2EBcXR4UKFVxHEpEg5ks5Z3Tuh/3fJ43JA4wEhmX2hay171prI621kWXLlvU9ZS5KSoKICGjQwHUSCRYbNmygefPmHD58mKVLl9K4cWPXkUQkyPlSzslA5bOeVwK2nfW8KFAPiDfG/ALcAMwN1oPCVqyAq66CQoVcJ5Fg8dVXX3Hy5EmWLVvGdddd5zqOiIQAX8p5BVDTGFPdGJMf6AnMPfNJa+1Ba20Za201a2014Cugk7U2KUcS5yBrvZGzprTFFydOnACgb9++/PTTT9SvX99xIhEJFZmWs7U2FRgELAR+BD621q4xxjxnjAmpVRW2bIE9e3SktmTu+++/p2bNmiQkJABQvHhxx4lEJJT4tAiJtXYeMO+cjz1znm2jsh/LjRUrvHuNnOVCvvnmG9q0aUOxYsWoVKmS6zgiEoK0tvZZkpIgXz64+mrXSSRQffHFF8TExFCyZEkSExO54oorXEcSkRCkcj5LUhLUrw8FCrhOIoHoxx9/pHXr1pQvX57ExESqVq3qOpKIhCiV82np6VoZTC7syiuvZMiQISQkJGg6W0RylMr5tE2b4OBB7W+WP1u8eDFbtmwhT548PP/881pgRERynMr5NK0MJhn5z3/+Q4cOHRg6dKjrKCISRlTOpyUlQcGC3gIkIgCzZs2ia9euXHPNNbz33nuu44hIGFE5n7ZihXf95ry6wrUAU6dOpUePHjRq1Ii4uDhKlizpOpKIhBGVM5CWBt99pylt8aSlpTFy5EiaNm3KggULtMCIiOQ6jROB9evh6FEdqS2Qnp5OREQE8+fPp2DBghQuXNh1JBEJQxo5o5XBxDN69Gg6derEiRMnKF26tIpZRJxROeMdDFakCNSq5TqJuPLGG29w//33kyePfiRExD39JsIr5+uu867jLOHnlVdeYciQIXTr1o2ZM2dSQEvEiYhjYV/Op07BypXa3xyuXn/9dR577DF69uzJtGnTyJ8/v+tIIiIq5zVr4Phx7W8OVy1atOD+++9nypQp5NV5dCISIMK+nLUyWPix1hIXFwfAddddxzvvvEOE9mmISAAJ+3JeuRKKFQNd+S88WGsZMmQIrVu3ZsGCBa7jiIhkKOzn8datg9q1wRjXSSSnpaenM2jQIMaMGcODDz5ImzZtXEcSEclQ2I+c16/3yllCW3p6OgMHDmTMmDE8+uijjBw5EqO/yEQkQIV1OR8+DMnJcOWVrpNITvviiy8YP348Tz/9NC+99JKKWUQCWlhPa//0k3evkXPoa9asGd999x0NGjRwHUVEJFNhPXJev96718g5NJ08eZLevXv/78AvFbOIBIuwLud16yBPHqhRw3US8bcTJ05w6623MnXqVDZs2OA6johIloT1tPb69XD55aDVGkNLSkoKXbt2ZcGCBYwePZp7773XdSQRkSwJ63Jet05T2qHm+PHjdOzYkaVLlzJu3DgGDBjgOpKISJaF7bR2erp3QJgOBgst+fPnp1atWkyaNEnFLCJBK2xHzlu2eGtqa+QcGg4ePMj+/fupVq0ao0ePdh1HRCRbwrac163z7jVyDn779++nTZs27N+/nzVr1ujKUiIS9MK2nHUaVWjYs2cPrVq1Yu3atcycOVPFLCIhIWzLed06KFkSypZ1nUQu1s6dO4mJiWHjxo18+umntG3b1nUkERG/CNtyPrOmtlZxDF6PPvoomzZtIjY2lpYtW7qOIyLiN2F7tLZOowp+b775JsuWLVMxi0jICctyPnQItm/XwWDB6JdffuHOO+8kJSWFEiVK0LhxY9eRRET8LizLWQeDBadNmzbRvHlzPvnkE37++WfXcUREckxYlrNOowo+69ev5+abb+bo0aMsXbqUq666ynUkEZEcE5YHhK1fD3nzwhVXuE4ivli7di3R0dFYa1m2bBlXX32160giIjkqbEfOl18O+fK5TiK+SE9Pp1y5csTHx6uYRSQshO3IWVPagS85OZmKFStSr149vv/+e/LkCcu/JUUkDIXdb7u0NO+CFzoYLLB9/fXX1KtXjxEjRgComEUkrITdb7xffoGTJzVyDmTLly+nVatWlC5dmltvvdV1HBGRXBd25azTqAJbfHw8bdq0oUKFCiQmJlK1alXXkUREcl3YlbNOowpc+/bto1OnTlStWpWEhAQqVqzoOpKIiBNhd0DY+vVQpgyULu06iZyrVKlSfPjhh9xwww2U1RVJRCSMhV05a03twDN37lwAOnXqRMeOHR2nERFxLyyntTWlHThmzZpFt27deO2117DWuo4jIhIQwqqc9++HXbs0cg4UU6dOpUePHjRq1IjY2FiMrt8pIgKEWTmfOVJbI2f3Jk2aRJ8+fWjWrBkLFy6kWLFiriOJiASMsCxnjZzdW7lyJdHR0cybN48iRYq4jiMiElDC6oCwdeu89bSrV3edJHwdPHiQ4sWLM2LECE6ePEmBAgVcRxIRCThhN3KuUUMXvHBl5MiR1K1bly1btmCMUTGLiJxHWJWzTqNy5+WXX2bo0KE0adKEChUquI4jIhLQwqacU1Nh40YdDObC888/z+OPP06vXr2YOnUq+TR1ISJyQWFTzps3w6lTGjnntvHjx/PMM89wxx13MHnyZPLmDavDHERELkrY/KbUmtpu3Hbbbezbt49hw4bpso8iIj4Km9+WOo0q91hrefvttzly5AhFixblkUceUTGLiGRB2PzGXLcOLr0USpZ0nSS0paenc9999zF48GCmTJniOo6ISFAKm2nt9es1pZ3T0tLSuPvuu3n//fd5/PHHGThwoOtIIiJBKaxGzprSzjmpqan89a9/5f333+fZZ59l+PDhWitbROQihcXI+fBh2LPHW4BEcsb27dtZvHgx//rXv/j73//uOo6ISFALi3Levt27v+wytzlC0alTp8ibNy+VK1dmzZo1lCpVynUkEZGgFxbT2jt2ePfly7vNEWqOHz9O165deeSRRwBUzCIifhJW5axVI/0nJSWFzp07ExsbS82aNV3HEREJKWE1ra2Rs38cPXqUjh07Eh8fz4QJE+jfv7/rSCIiISUsynnHDu9KVJp1zT5rLZ07dyYhIYEPPviAPn36uI4kIhJywmJae/t2b9SsM3uyzxjDfffdx9SpU1XMIiI5JGxGzprSzp59+/bxzTff0LZtW7p27eo6johISAuLct6+HapVc50ieO3Zs4dWrVqxYcMGNm/eTNmyZV1HEhEJaWExra2R88XbuXMnUVFRrFu3jtmzZ6uYRURyQciPnFNTYfdulfPF2LZtGy1btmTLli189tlnREdHu44kIhIWQr6cd+0Ca3WO88WYOnUqycnJLFiwgJtuusl1HBGRsBHy09paHSzrrLUADB06lFWrVqmYRURymU/lbIxpa4xZb4zZaIx5PIPPDzXGrDXGrDLGLDHGVPV/1ItzZgESjZx9s3HjRho1asSPP/6IMYbq1au7jiQiEnYyLWdjTAQwCmgH1AV6GWPqnrPZ90CktbY+MBN4xd9BL5ZGzr5bt24dN998M7/88gsnTpxwHUdEJGz5MnJuBGy01v5srT0JTAM6n72BtXaZtfbY6adfAZX8G/PinSnncuXc5gh0q1evJioqirS0NJYtW0aDBg1cRxIRCVu+lHNFYOtZz5NPf+x8BgDzM/qEMeYeY0ySMSZp9+7dvqfMhu3boWRJKFgwV14uKK1bt44WLVqQJ08eEhISqFevnutIIiJhzZdyzmjRS5vhhsb0ASKBVzP6vLX2XWttpLU2MrfOl9U5zpmrUqUKbdq0ISEhgdq1a7uOIyIS9nw5lSoZqHzW80rAtnM3MsbEAE8Cza21AbPDcvt2HQx2Pt9//z2XX345xYsXZ8qUKa7jiIjIab6MnFcANY0x1Y0x+YGewNyzNzDGXAuMBTpZa3f5P+bF08g5Y//3f//HzTffzP333+86ioiInCPTcrbWpgKDgIXAj8DH1to1xpjnjDGdTm/2KlAEmGGMWWmMmXueL5errFU5ZyQ+Pp42bdpQsWJFXn75ZddxRETkHD6tEGatnQfMO+djz5z1OMbPufzi8GE4dkzT2meLi4ujc+fOVK9enSVLllBef7mIiASckF4hTOc4/9GpU6e47777qFmzJvHx8SpmEZEAFdJra2t1sD/Kly8fCxYsoESJEpQuXdp1HBEROQ+NnMPAjBkzGDx4MNZarrjiChWziEiAC4tyDueR84cffkjPnj35/vvvSUlJcR1HRER8ENLlvH075MvnrRAWjiZOnEjfvn1p3rw58+fPp1ChQq4jiYiID0K6nM+cRmUyWuMsxI0bN47+/fsTExNDbGwsRYoUcR1JRER8FNLlHM6rg1WoUIEuXbowd+5cjZhFRIJMSJdzOC5A8uOPPwLQoUMHZs+eTUFd8UNEJOiEfDmH08j5xRdfpF69eiQmJrqOIiIi2RCy5ZyaCrt3h8fI2VrLP//5T/7+97/Tq1cvmjRp4jqSiIhkQ8guQrJrl7e2dqiXs7WWJ598khdffJG//vWvjBs3joiICNexREQkG0J25Bwuq4PFxcXx4osvMnDgQMaPH69iFhEJASE7cg6X1cFatWrFf/7zHzp06IAJx3PGRERCUMiOnEN5dbD09HQeeeQRfvjhB4wx3HLLLSpmEZEQErIj5zPT2uXKuc3hb2lpadx1111MnDiRkiVLcs0117iOJCIifhay5bxjh7dsZ4ECrpP4T2pqKnfccQcfffQR//jHP3jiiSdcRxIRkRwQsuUcaquDnTp1ittvv50ZM2YwfPhwFbOISAgL2XIOtdXB0tPTOXz4MK+//jpDhw51HUdERHJQSJfzjTe6TpF9x48fJyUlhZIlSxIbG6tTpUREwkBIHq1trTetHewj52PHjtGpUyfatm1LamqqillEJEyE5Mj58GFISQnucj5y5AgdO3YkISGBCRMmkDdvSP6vEhGRDITkb/ydO737YD2N6tChQ7Rv354vv/ySKVOm0Lt3b9eRREQkF4VkOe/a5d0HazkPHDiQr7/+mmnTptG9e3fXcUREJJeF5D7nMyPnSy91m+NivfTSS8yZM0fFLCISpkKynINx5Lx7926ee+450tPTqVq1Kh06dHAdSUREHAnJcj4zci5Txm0OX+3YsYOoqChefPFF1q5d6zqOiIg4FrL7nEuXhnz5XCfJ3G+//UZ0dDTJycnMmzePevXquY4kIiKOhWQ579wZHPubt2zZQnR0NLt27WLhwoU0a9bMdSQREQkAIVnOu3YFRzn//PPPHD16lEWLFnHDDTe4jiMiIgEiJMt5505o0MB1ivM7cuQIRYoUISoqik2bNlGoUCHXkUREJICE5AFhgTxyXrduHbVr12by5MkAKmYREfmTkCvnEyfgwIHAPI1q9erVNG/enNTUVK677jrXcUREJECFXDnv3u3dB9rIeeXKlURFRZE3b14SEhK46qqrXEcSEZEAFXLlHIjrau/cuZPo6GgKFSpEQkICV155petIIiISwEKunM+sDhZII+dy5crx3HPPkZiYSI0aNVzHERGRABdyR2sH0tKdn3/+OQULFuT6669n0KBBruOIiEiQCLmRc6Bc9GLJkiW0bduWhx56CGut2zAiIhJUQq6cd+2CQoWgSBF3GRYuXMgtt9zC5ZdfzuzZszHGuAsjIiJBJ+TK2fXSnbGxsXTq1InatWuzbNkyygXC/LqIiASVkCvnXbvc7m+eNGkS9evXZ8mSJZQJlstiiYhIQAm5A8J27oQqVXL/dVNTU8mbNy9Tpkzh+PHjFC9ePPdDiIhISNDI2Q+mTJlCo0aN2Lt3LwUKFFAxi4hItoRUOaen5/662hMmTKBfv36UKFGCggUL5t4Li4hIyAqpct6/H9LScq+cx44dy4ABA2jVqhWxsbEULlw4d15YRERCWkiVc24u3fnBBx/wt7/9jQ4dOvDpp5/q6lIiIuI3IVXOubl0Z8uWLRk8eDCzZ8/WdLaIiPhVSJVzboyc58yZQ1paGhUrVuTNN98kf/78OfdiIiISlkKqnHNy5Gyt5R//+AddunRhwoQJ/n8BERGR00LqPOedOyFPHihd2r9f11rL3//+d1566SX69+/PnXfe6d8XEBEROUtIlfOuXVC2rFfQ/mKtZdiwYYwcOZK//e1vjBo1ijz+fAEREZFzhFTL7Nzp//3NGzduZOzYsQwePJjRo0ermEVEJMeF3MjZX/ubrbUYY6hZsyYrV66kRo0aurqUiIjkipAaBvpr5JyWlsadd97J6NGjAahZs6aKWUREck1IlbM/Rs6pqan07duXiRMnsnfvXv8EExERyYKQmdY+etS7ZWfkfPLkSXr37s2sWbN46aWXeOyxx/wXUERExEchU87ZPcc5PT2d7t27M3fuXEaMGMGQIUP8F05ERCQLQqacs7s6WJ48eWjevDmtW7fm/vvv918wERGRLAqZcr7YkfOxY8f46aefaNCgAUOHDvV/MBERkSwKmQPCLqacjxw5Qvv27WnRogX79+/PmWAiIiJZFDIj5zPT2r6W86FDh2jXrh1ff/01kydPpmTJkjkXTkREJAtCppx37YJixcCXqzfu37+ftm3b8t133zF9+nS6deuW8wFFRER8FFLl7Ouo+fXXX+f7779n1qxZdOrUKWeDiYiIZFHIlPP+/VCqlG/bPvvss3Tq1IlGjRrlbCgREZGLEDIHhO3fDxfabbx9+3a6du3Kzp07yZcvn4pZREQCVsiMnA8cgOrVM/5ccnIy0dHRbNu2jU2bNlHO35euEhER8aOQGTkfOAAlSvz547/++ivNmzdnx44dLFy4kCZNmuR+OBERkSwIiZGztRmX8+bNm4mKiuLQoUMsXrxYU9kiIhIUQmLknJICJ0/+eZ9zoUKFqFy5MkuWLFExi4hI0AiJkfOBA979mZHz5s2bqVSpEuXKlePzzz/XtZhFRCSohMTI+exy/u9//0vjxo156KGHAFTMIiISdHwqZ2NMW2PMemPMRmPM4xl8voAxZvrpz39tjKnm76AXcqac9+zZSIsWLcifPz8PPvhgbkYQERHxm0zL2RgTAYwC2gF1gV7GmLrnbDYA2G+trQGMBF72d9ALOVPOjz8+kMKFC5OQkECtWrVyM4KIiIjf+DJybgRstNb+bK09CUwDOp+zTWdg0unHM4GWJhfnk3ftOgV409qJiYlcccUVufXSIiIifudLOVcEtp71PPn0xzLcxlqbChwESp/7hYwx9xhjkowxSbt37764xBkoWTIfdeoc4bPPJlO1alW/fV0REREXfCnnjEbA9iK2wVr7rrU20lobWbZsWV/y+aRzZ1i7tgj161/mt68pIiLiii/lnAxUPut5JWDb+bYxxuQFigP7/BFQREQk3PhSziuAmsaY6saY/EBPYO4528wF7jj9+FZgqbX2TyNnERERyVymi5BYa1ONMYOAhUeDQ2wAAAUISURBVEAEMMFau8YY8xyQZK2dC4wHJhtjNuKNmHvmZGgREZFQ5tMKYdbaecC8cz72zFmPjwPd/RtNREQkPIXECmEiIiKhROUsIiISYFTOIiIiAUblLCIiEmBUziIiIgFG5SwiIhJgVM4iIiIBRuUsIiISYFTOIiIiAUblLCIiEmBUziIiIgFG5SwiIhJgVM4iIiIBRuUsIiISYIy11s0LG7Mb+NWPX7IMsMePXy9c6X3MPr2H2af3MPv0Hmafv9/Dqtbasr5s6Kyc/c0Yk2StjXSdI9jpfcw+vYfZp/cw+/QeZp/L91DT2iIiIgFG5SwiIhJgQqmc33UdIETofcw+vYfZp/cw+/QeZp+z9zBk9jmLiIiEilAaOYuIiIQElbOIiEiACbpyNsa0NcasN8ZsNMY8nsHnCxhjpp/+/NfGmGq5nzKw+fAeDjXGrDXGrDLGLDHGVHWRM5Bl9h6etd2txhhrjNEpLRnw5X00xtx2+vtxjTHmo9zOGOh8+HmuYoxZZoz5/vTPdHsXOQOVMWaCMWaXMWb1eT5vjDFvnX5/VxljrsuVYNbaoLkBEcAm4HIgP/ADUPecbe4D/n36cU9guuvcgXTz8T1sARQ6/fhevYdZfw9Pb1cUSAS+AiJd5w60m4/fizWB74GSp59f6jp3IN18fA/fBe49/bgu8Ivr3IF0A24GrgNWn+fz7YH5gAFuAL7OjVzBNnJuBGy01v5srT0JTAM6n7NNZ2DS6cczgZbGGJOLGQNdpu+htXaZtfbY6adfAZVyOWOg8+X7EOB54BXgeG6GCyK+vI93A6OstfsBrLW7cjljoPPlPbRAsdOPiwPbcjFfwLPWJgL7LrBJZ+AD6/kKKGGMqZDTuYKtnCsCW896nnz6YxluY61NBQ4CpXMlXXDw5T082wC8vxrld5m+h8aYa4HK1trY3AwWZHz5XqwF1DL/3979hNgUxmEc/z66ZMFultRYWKhZWrCiSLKYlQU1MbJlIVlZkK1kSyJlodhwF2o2UhJltpSa0DRlIWU2Sv48Fu8NzWjmVXPvPWfm+azurVPn13PvOb/O731vV3om6YWkgwOrrh1qMrwITEiaAx4BpwdT2qrxv/fMFdHp9wlW2L+egBf+FqzmmLWsOh9JE8BOYE9fK2qfJTOUtA64CkwOqqCWqvkudiij7b2UCc5TSWO2P/e5traoyfAocNv2FUm7gTu9DH/2v7xVYSg9pW1PznPA1r/eb2HxiOb3MZI6lDHOUiOLtaYmQyTtB84D47a/Dqi2tlguw83AGPBE0nvKOlU3m8IWqb2eH9r+Zvsd8IbSrKOoyfAkcA/A9nNgI+UPHaJO1T1zpbWtOb8EtkvaJmkDZcNXd8ExXeB47/Vh4LF7q/oBVGTYG8lepzTmrPEttmSGtudtj9getT1KWbcftz09nHIbq+Z6fkDZoIikEcqY++1Aq2y2mgxngX0AknZQmvPHgVbZbl3gWG/X9i5g3vaHfp+0VWNt298lnQKmKLsUb9l+JekSMG27C9ykjG1mKE/MR4ZXcfNUZngZ2ATc7+2lm7U9PrSiG6Yyw1hGZY5TwAFJr4EfwDnbn4ZXdbNUZngWuCHpDGUcO5kHlj8k3aUsm4z01uUvAOsBbF+jrNMfAmaAL8CJgdSVzygiIqJZ2jbWjoiIWPXSnCMiIhomzTkiIqJh0pwjIiIaJs05IiKiYdKcIyIiGibNOSIiomF+AZ1JhTLwF3v4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "tp=output[:,1][np.where(Y_test==1)[0]]\n",
    "fp=output[:,1][np.where(Y_test==0)[0]]\n",
    "\n",
    "def sig_eff(tp,threshold):\n",
    "    \n",
    "    sig_eps = np.where(tp>threshold)[0].shape[0]/tp.shape[0]\n",
    "    return sig_eps\n",
    "\n",
    "def bkg_eff(fp,threshold):\n",
    "    \n",
    "    bkg_eps = np.where(fp>threshold)[0].shape[0]/fp.shape[0]\n",
    "    return bkg_eps\n",
    "\n",
    "threshold_range=np.linspace(0.0,1.,num=30)\n",
    "sig_eps_vals=[sig_eff(tp,threshold_range[i]) for i in range(len(threshold_range))]\n",
    "bkg_eps_vals=[bkg_eff(fp,threshold_range[i]) for i in range(len(threshold_range))]\n",
    "\n",
    "\n",
    "plt.plot(threshold_range,threshold_range, 'black', linestyle='dashed')\n",
    "plt.plot(bkg_eps_vals,sig_eps_vals,'b',label=\"NN ROC Curve\")\n",
    "fig=plt.gcf()\n",
    "fig.set_size_inches(8,8)\n",
    "pAUC=roc_auc_score(Y_test,output[:,1])\n",
    "print(\"pAUC {0}\".format(pAUC))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39.04553343189967"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S=X_test[np.where(output[:,1]>0.6)].shape[0]\n",
    "B=X_test[np.where(output[:,1]<0.6)].shape[0]\n",
    "S/(np.sqrt(S+B))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}