drugs_classification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import sys \n",
    "from atol import Atol\n",
    "from heat_kernel_homology import GeneralFiltrationGraphHomology\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import pandas as pd\n",
    "from pysmiles import read_smiles\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "from networkx.convert_matrix import to_numpy_matrix\n",
    "from utils import blockPrinting, plot_diagram\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data load"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the ClinTox dataset which contains the data of drugs approved by the FDA and those that have failed clinical trials for toxicity reasons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('./data/clintox/clintox.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>smiles</th>\n",
       "      <th>FDA_APPROVED</th>\n",
       "      <th>CT_TOX</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>*C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[C@@H]1([C@@H]([C@@H]([C@H]([C@@H]([C@@H]1Cl)C...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[C@H]([C@@H]([C@@H](C(=O)[O-])O)O)([C@H](C(=O)...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[H]/[NH+]=C(/C1=CC(=O)/C(=C\\C=c2ccc(=C([NH3+])...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[H]/[NH+]=C(\\N)/c1ccc(cc1)OCCCCCOc2ccc(cc2)/C(...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              smiles  FDA_APPROVED  CT_TOX\n",
       "0            *C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC             1       0\n",
       "1  [C@@H]1([C@@H]([C@@H]([C@H]([C@@H]([C@@H]1Cl)C...             1       0\n",
       "2  [C@H]([C@@H]([C@@H](C(=O)[O-])O)O)([C@H](C(=O)...             1       0\n",
       "3  [H]/[NH+]=C(/C1=CC(=O)/C(=C\\C=c2ccc(=C([NH3+])...             1       0\n",
       "4  [H]/[NH+]=C(\\N)/c1ccc(cc1)OCCCCCOc2ccc(cc2)/C(...             1       0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Smiles are a standard way of encoding chemical information into a string. For further information and for a dictionary of the characters used visit the following [link](https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Smiles converted into Graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are going to analyze the problem in two different ways: in the first case we threat the **drugs molecules** as **graphs** (each atom is a node and the edges corresponds to chemical bounds), considering only the structure of the molecules (the information about the node elements is not considered in this part); in the second example we also add chemical information as *composition* and *stereoisomerism*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "@blockPrinting\n",
    "def convert_smiles(smiles):\n",
    "    return [to_numpy_matrix(read_smiles(smile, explicit_hydrogen=False, \n",
    "                                        zero_order_bonds=True, \n",
    "                                        reinterpret_aromatic=True)) \n",
    "            for smile in smiles]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['smiles_graph'] = convert_smiles(df.smiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>smiles</th>\n",
       "      <th>FDA_APPROVED</th>\n",
       "      <th>CT_TOX</th>\n",
       "      <th>smiles_graph</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>*C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[C@@H]1([C@@H]([C@@H]([C@H]([C@@H]([C@@H]1Cl)C...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1.]]]], ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[C@H]([C@@H]([C@@H](C(=O)[O-])O)O)([C@H](C(=O)...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[H]/[NH+]=C(/C1=CC(=O)/C(=C\\C=c2ccc(=C([NH3+])...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[H]/[NH+]=C(\\N)/c1ccc(cc1)OCCCCCOc2ccc(cc2)/C(...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              smiles  FDA_APPROVED  CT_TOX  \\\n",
       "0            *C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC             1       0   \n",
       "1  [C@@H]1([C@@H]([C@@H]([C@H]([C@@H]([C@@H]1Cl)C...             1       0   \n",
       "2  [C@H]([C@@H]([C@@H](C(=O)[O-])O)O)([C@H](C(=O)...             1       0   \n",
       "3  [H]/[NH+]=C(/C1=CC(=O)/C(=C\\C=c2ccc(=C([NH3+])...             1       0   \n",
       "4  [H]/[NH+]=C(\\N)/c1ccc(cc1)OCCCCCOc2ccc(cc2)/C(...             1       0   \n",
       "\n",
       "                                        smiles_graph  \n",
       "0  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....  \n",
       "1  [[[[[0. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1.]]]], ...  \n",
       "2  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1....  \n",
       "3  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....  \n",
       "4  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we show an example of conversion from **SMILES** to **graph**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I don't quite know how to handle stereo yet...\n",
      "*C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfXRU1d0v8O+ZF2aiZEAlJkjyGCTAQDWBBCFaIeHtUhMeCxI0aArPdd2ijX2z15elsVevmt7i9dbaPqAt1drHCIEVaSwQUFFCoBqQAEEgQwiSNtEEk2iYBDPDvJz7R5pAIG9zzplzzsx8P2u57CrOZhMg37P32b/9E0RRFEFERBQhDFpPgIiISE0MPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigMPiIiiigmrScQqVo73SipaoSj2Qmnywub1QR7nA3L0+Jx3UiL1tMjIgpbAvvxqau6oR1ry+uwp7YFAOD2+nt/zGoyQASQOTkG+RlJSEkYrdEsiYjCF4NPRUWV9Sgsc8Dl9WGwr7ogAFaTEQVZduSlJ6o2PyKiSMCtTpV0h14Nujz+If9bUQS6PD4UltUAAMOPiEhBPNyiguqGdhSWOYYVepfq8vhRWObA0cb2IM2MiCjyMPhUsLa8Di6vT9JnXV4f1pXXKTwjIqLIxeALstZON/bUtgz6Tm8wogjsPtmCtk63shMjIopQDL4gK6lqlD2GAKDkkPxxiIiIwRd0jmZnn5IFKVxePxxNHQrNiIgosjH4gszp8io0jkeRcYiIIh2DL8hsVmUqRmxWsyLjEBFFOgZfkNnjbLCY5H2ZrSYD7GOjFZoREVFkY/AFWU5avOwxRAA5qfLHISIiBl/QjRlpQcakGAiCtM8LAjB3cgwvriYiUgiDTwUPZybBajJK+qzVZER+ZpLCMyIiilwMPhWkJIxGQZYdUebAvtxRZgMKsuxIjmeXBiIipfCSapX0XDTN7gxERNpiWyKVHW1sx7ryOuw+2QIB3cXpPQz/eg+YdP1I/DJrKmZPitFmkkREYYzBp5GK2hY8v/0ETrd0QhS7T272YENaIqLgYfBpgA1piYi0w3d8KmNDWiIibfFUp4rYkJaISHsMPhWxIS0RkfYYfCphQ1oiIn1g8KmEDWmJiPSBwacSNqQlItIHBp9K2JCWiEgfGHwqYUNaIiJ9YPCphA1piYj0gcGnEjakJSLSBwafStiQlohIHxh8KmJDWiIi7TH4VMSGtERE2uMl1SpjQ1oiIm2xLZFGBmtI29OPb+7kGORnJnGlR0SkIAafxto63Sg51Ij/s+4vuC1jPuKuHQX72GjkpMbzIAsRURAw+HQiNjYW1dXViIuL03oqRERhjYdbdKKrqwtRUVFaT4OIKOzxcIuGWjvdKKlqhKPZiahFv8Avy+ow9YbRWJ7GbU4iomDhVqcGqhvasba8DntqWwCgT9eGnoMtmZNjkJ+RhJQEHmwhIlISg09lRZX1LGUgItIQtzpV1B16NejyDN2XTxSBLo8PhWU1AMDwIyJSCA+3qKS6oR2FZY5hhd6lujx+FJY5cLSxPUgzIyKKLAw+lawtr4PL65P0WZfXh3XldQrPiIgoMjH4VNDa6cae2pZB3+kNRhSB3Sdb0NbpVnZiREQRiMGngpKqRtljCABKDskfh4go0jH4VOBodvYpWZDC5fXD0dSh0IyIiCIXg08FTpdXoXE8ioxDRBTJGHwqsFmVqRqxWc2KjENEFMkYfCqwx9lgMcn7UltNBtjHRis0IyKiyMXgU0FOWrzsMUQAOanyxyEiinQMPhWMGWlBxqQYCIK0zwtCd1NaXlxNRCQfg08lD2cmwWoySvqs1WREfmaSwjMiIopMDD6VpCSMRkGWHVHmwL7kUWYDCrLsSI5nlwYiIiXwkmoV9Vw0ze4MRETaYVsiDRxtbMfa3XXYebQRFssIXPBd/C3o6cc3d3IM8jOTuNIjIlIYg08jBw4cwMof5uOR32+Co6kDTpcHNqsZ9rHRyEllB3YiomDhVqdGSktLsTRrIR6cM0HrqRARRRQebtHIX//6VyxdulTraRARRRwGnwYcDgecTidmzJih9VSIiCIOg08DpaWlWLJkCQwGfvmJiNTGd3xB1trpRklVIxzNTjhdXtisJuyqOovn/uP7Wk+NiCgi8VRnkFQ3tGNteR321LYAQJ9+fKLHDWtUFDInxyA/IwkpCSxZICJSC4MvCIoq61mkTkSkU9zqVFh36NWgyzN0x3VRBLo8PhSW1QAAw4+ISAU8XaGg6oZ2FJY5hhV6l+ry+FFY5sDRxvYgzYyIiHow+BS0trwOLq9P0mddXh/WldcpPCMiIrocg08hrZ1u7KltGfSd3mBEEdh9sgVtnW5lJ0ZERH0w+BRSUtUoewwBQMkh+eMQEdHAGHwKcTQ7+5QsSOHy+uFo6lBoRkRE1B8Gn0KcLq9C43gUGYeIiPrH4FOIzapMZYjNalZkHCIi6h+DTyH2OBssJnlfTqvJAPvYaIVmRERE/WHwKSQnLV72GCKAnFT54xAR0cAYfAoZM9KCjEkxEARpnxcEYO7kGHZeJyIKMgafgh7OTILVZJT0WavJiPzMJIVnREREl2PwKSglYTQKsuyIMgf2ZY0yG1CQZUdyPLs0EBEFG4NPYXnpiSjImgKr2QD4B6/rEwQgymxEQdYUXlBNRKQSdmcIgrz0RPhb6/HMpo8xInE6BHQXp/ewmgwQ0f1OLz8ziSs9IiIVsR9fkCxbtgzz5s1D7qr/gZJDjXA0dcDp8sBmNcM+Nho5qfE8yEJEpAEGXxCcPn0as2bNwj/+8Q9cffXVWk+HiIguwXd8QfC73/0OP/zhDxl6REQ6xBWfwtrb23HTTTfhs88+w7hx47SeDhERXYYrPoWtX78eWVlZDD0iIp3iik9BHo8HN910E0pLS5GWlqb1dIiIqB8sZ5CgtdONkqpGOJqdcLq8sFlNsMfZYPznp5gwYQJDj4hIx7jiC0B1QzvWltdhT20LAPRpPGs1GeByu5E8xojnV8xGSgJr84iI9IjBN0xFlfUoLHPA5fVhsK+YIHTfu1mQZedtLER0hYF2jJansbZXLQy+YegOvRp0eQa/guxS3fdv8ioyIuo21I6RCCBzcgzyM5K4YxRkDL4hVDe0I3d9Jbo8voA/G2U2YtPqdF5JRhThuGOkLyxnGMLa8jq4vIGHHgC4vD6sK69TeEZEFEou7hgNHnoAIIpAl8eHwrIaFFXWqzK/SMTgG0Rrpxt7aluG/MM6EFEEdp9sQVunW9mJEVFIqG5oR2GZI6DXJADQ5fGjsMyBo43tQZpZZGPwDaKkqlH2GAKAkkPyxyGi0MMdI31i8A3C0ezs8wJaCpfXD0dTh0IzIqJQwR0j/WIB+yCcLq9C43gUGYeIQoeSO0YPzpkgf0IDiMTyCgbfIGxWZb48NqtZkXGIKHTofcdo8PKKZry8qzZsyysYfIOwx9lgMTXL+sNrNRlgHxut4KyIKBQotWN0Lgg7RkOVV7j+9T3v/RNnUVHbGnblFXzHN4ictHjZY4gAclLlj0NEoUWpHaP3tv4VCxcuxGOPPYYNGzbgxIkT8PmkHZgBWF4BMPgGNWakBRmTYiAI0j4vCMDcyTFhu09ORAPr3jGS9y3WajLg5/+xHL/4xS9wzTXXoLS0FN///vdhs9kwa9YsPPTQQ3jttdewf/9+fPvtt0OOx/KKbry5ZQi8uYWIpGjtdOO7az6S9arEYjLg4yfmXfHw3NHRgerqahw5cgSHDx/GkSNHUFNTg8TEREybNg3Tp0/v/feYMWN6P7f6rYP4oOaspJOmggAsmhqL1/JmSP716AWDbxh4VycRSaFm0Fy4cAEOh6M3CHv+PXLkSEyfPh2TU2ag1J8GryhxCwsDB3GoYfANU1FlPZ7ffgKuCz4IhoG3L3jXHhH10HrHSBRF1NfX48iRI3jrYBMOe8dBNEh/92g1GfDIwklBLa9QA9/xDVNeeiLsjTuRaDoHi8kA62V791aTARaTAYumxmLT6nSGHhEhJWE0CrLsiDIH9q22e8fILvs1iSAIGD9+PJYuXYp/S75NVugB4XMhB8sZhqm6uhr7d2zGqVOFuCCMQMmhRjiaOuB0eWCzmmEfG42c1PAt+CQiaXoegrXuzsALOS5i8A1TQUEBnnrqKURHd9fkhfpSn4jUk5eeiOT40VhXXofdJ1sg4GKtHHCxH9/cyTHIz0wKyoE4XshxEYNvGPbt24djx47hnXfe0XoqRBSikuNH47W8GWjrdGuyY8QLOS7i4ZYhiKKI2bNnY/Xq1Vi5cqXW0yEikiSY5RWhhodbhlBWVob29nbcf//9Wk+FiEgyuRdyQBSRMXFMyIcewOAblN/vx5NPPonCwkIYjUatp0NEJMvDmUmwmqR9LxNELw4V/RoNDQ0Kz0p9mr3jC4VWGBs3bsTVV1+Nu+66S+upEBHJ1lNeIeVCjifvTEHTtSdx66234vXXX0d2dnYQZxpcqr/jG7wVRvfJJjVbYQwUwEuSY3F7WjL+/Oc/IyMjI+jzICJSy1DdGXr0V16xb98+rFixAvfddx9eeOEFmM2hd8pT1eCT88VW2lAB7PF6EdV+Bm8/tTLselERER1tbJdcXtHa2oof/OAH6OjowMaNG5GQkNDvz6HXnT3Vgk9P910OO4ABWM28foyIwpfU8gq/348XX3wRv/3tb/HGG28gKyur98f0trN3OVWCT+v76i6lpwAmIgp1e/fuxX333Yf7778fzz//PDZVfaGbnb0Bf241gk8vrTD0FMBEROGipaUFK1euRNPVE9A1JQtu7/C/2WuxsAh6OUNrpxt7alskhR7Q3QF498kWtHW6Zc9lbXkdXF5pnYtdXh/WldfJngMRUbiJiYnBr/7wNjonLgoo9ABtmtwGPfhKqhpljyEAKDkkbxw9BTARUbh5dc/n8AnSIkXthUXQ6/gczU5ZV+QA3aeN/vLuLtS+uw5RUVGS/tlUdVb2r6UngHlBNRHRRUouLNQ47Rn04FOqFUaU7VqMvXosurq6cP78ebS2tqKrq2vY/1gzH8RV38mUNYdw6UVFRKQkJXf21FhYBD34lGqFkWyfiEfvXS758w/85VN85PhK9jzCoRcVEZGSlNrZU2thEfR3fN2tMOT9NEq0wmAvKiKi4Ai1JrdBD76ctHjZY4gAclLljaOXACYiCjehtrAIevDJbYUhCN3X5sh94amXACYiCjehtrBQpS2RnFYYVpMR+ZlJsueglwAmIgo3obawUCX4elphRJkD++m6K/rtit2WoocAJiIKN6G2sFCtEW1eeiIKsqYgymwc8osjCN1XhCl9jY1eApiIKNyE0sLC+Oyzzz6r1k+WHD8acyaOwTfnL6Dhmy6YDQK8/osVj1aTAUaDgAVTrseLy5KxcGpcUOYwOsqMTz7/Gr4hqi0FBCeAiYjCTdwoK0ZHmfDJ5219vq8PpeeuzmB8vx+I6o1oe0hthaGUoXpRXfB6cb33K/zxZ8u40iMiGiY99V0d8OfWKvj0YqAA/u4NJsyZOR1nzpzBqFGjtJ4mEVHIkNPkVg0RH3yDWbFiBWbNmoWf//znWk+FiCjkXLqwqDldj5Yv/4nV92SrtrM3EAbfID755BPk5eWhtrYWRqO0l7ZERATs378f+fn5qKqq0noq6p3qDEXp6em49tprsWPHDq2nQkQU0iZMmIDTp09DD2utoF9SHcoEQcDPfvYzvPLKK1i8eLHW0yEiClmiZSRGpGQh/639cItG2Kwm2ONsWJ6m/rYntzqH4Ha7kZiYiA8//BBTp07VejpERCGluqEda8vrsKe2BW63CzCO6P2xnoMumZNjkJ+RhJQEdQ66MPiG4dlnn8XZs2fx6quvaj0VIqKQodfSBgbfMDQ1NWHq1Kk4+JkDH5zuhKPZCafLq+lSnYhIz7pDrwZdnuH36espZg92+DH4hqG6oR3//aVinLtqHEwmU5+Gi1ot1YmI9Kq6oR256yvR5fEF/NkosxGbVqcHtb6PpzqHUFRZj9z1lfjmqgT4YLiiy7DL64fb68f7J84id30liirrtZkoEZFOrC2vg8sbeOgBgMvrw7ryOoVn1BeDbxAXl+o+DLUsFkWgy+NDYVkNw4+IIlZrpxt7alsGfac3GFEEdp9sQVunW9mJXYLBN4DqhnYUljkC2p8GgC6PH4VlDhxtbA/SzIiI9KukqlH2GAKAkkPyxxkIg28Ael+qExHpkaPZecUroUC5vH44mjoUmtGVGHz9CIWlOhGRHjldXoXG8SgyTn8YfP0IhaU6EZEe2azKXAhms5oVGac/DL5+hMJSnYhIj+xxNlhM8qLFajLAPjZaoRldicHXj1BYqhMR6VFOWrzsMUQAOanyxxkIg68fobBUJyLSozEjLciYFANBkPZ5QehuUhvM27AYfP0IhaU6EZFePZyZBKtJWg9Tq8mI/MwkhWfUF4OvH6GwVCci0quUhNEoyLIjyhxYxHTf1WkP6nVlAIOvX6GwVCci0rO89EQ8MjcRoteNob6VCkL3HZ1qXFANMPgGpPelOhGR3p149w9I76zEou/EwmIywHrZKySryQCLyYBFU2OxaXW6KqEHsDvDoKS01RD8Hjx7VzJWffemIM6MiEjfjh07hnnz5uH48eOIiYlBW6cbJYca4WjqgNPlgc1qhn1sNHJS2YFddwJppGgxGWCr+wCptm/x6quvQpC6V0pEFMJEUcT8+fNx991348c//rHW07kCg28Yjja2Y115HXafbIGA7uL0Hj39+OZOjkF+ZhISbQbMmTMH9957L5588knN5kxEpJXNmzfjhRdewKFDh2AyKVMepiQGXwCGu1T/4osvcPvtt6OwsBB5eXkazpiISF3nz5/HlClTUFRUhDlz5mg9nX4x+ILk+PHjmDdvHjZs2ID58+drPR0iIlU8/fTT+Pzzz7FhwwatpzIgBl8QlZeX45577sGHH36IW265RevpEBHJ1trpRklVIxzNTjhdXtisJtjjbFieFo9vmhswa9YsHD16FOPGjdN6qgNi8AVZcXExHn/8cXz88ceIj2dBOxGFpuqGdqwtr8Oe2hYA6HORf89ZB0tbHRbc4Mdvnv6ZRrMcHv29dQwzubm5aGhowJ133ol9+/Zh1KhRvT822JMTi9+JSC+GOt3ec+DPHZ2IsgsmpFbWq1aTJwVXfCoQRRE//elPcfz4cezcuRM1Z78d8skpc3IM8jOSkJIQ3Kt7iIgGI6WeufvqMXVuYZGCwacSn8+HnJwcOK9PwRexs+D2+oesC7SajCjIsuv2Dw8Rhbfqhnbkrq9El8cX8GejzEZsWp0e9Hs3peCVZSoxGo3490fW4PSoFLg8g4ceAIgi0OXxobCsBkWV9arMkYjoUmvL6+DyBh56AODy+rCuvE7hGSmDwaeS6oZ2/N9dnwPGEQF9rsvjR2GZA0cb24M0MyKiK7V2urGntmXIh/SBiCKw+2QL2jrdyk5MAQw+lYTrkxMRhaeSqkbZYwgASg7JH0dpDD4VhPOTExGFJ0ezs8/BOylcXj8cTR0KzUg5DD4VhPOTExGFJ6fLq9A4HkXGURKDTwXh/OREROHJZlWmzNtmNSsyjpIYfCoI5ycnIgpP9jgbLCZ5EWE1GWAfG63QjJTD4FNBOD85EVF4ykmTf8WiCCAnVX9XNTL4VBDOT05EFJ7GjLQgY1IMpPbTFoTuPqV6vH6RwaeCcH5yIqLw9XBmEkYYpCWf1WREfmaSwjNSBoNPBXKfnCD6MevfRuryyYmIwte5M0fhrHgTZiGwWqzuuzrturyuDGDwqebhzCRYTUZJnzUKInb+v0fw0ksvwePhARciCr53330Xy5Ytw1v/azWeuetmRJmNQz68C0L3HZ16vqAaYPCpJiVhNAqy7IgyB/YljzIb8L+/n4y/b92IXbt2Yfr06aioqAjSLImIgDfeeAMPPfQQduzYgQULFiAvPRGbVqdj0dRYWEwGWC87s2A1GWAxGbBoaiw2rU7XdegB7M6guqH6WvXorzuDKIp455138Mgjj2DevHl48cUXERsbq87EiSjsiaKINWvW4A9/+APee+89TJo06Yr/pq3TjZJDjXA0dcDp8sBmNcM+Nho5qaHTR5TBp4Gjje1YV16H3SdbIOBiE0fgYj++uZNjkJ+Z1O8eeUdHB5577jm8+eabeO6557B69WoYjdK2UYmIAMDv9+PRRx/FBx98gJ07d2LcuHFaTyloGHwakvvkdOzYMeTn56Orqwvr1q3DrbfeqsKsiSjceDwePPDAAzhz5gy2bt2Ka665RuspBRWDL8SJooi33noLTzzxBJYuXYrCwsKw/0NLRMo5f/48cnJyYDKZsGnTJlx11VVaTynoGHxh4ptvvsHTTz+NLVu24Ne//jVWrlwJYYgjWK2dbpRUNcLR7ITT5YXNaoI9zoblaaGzV09E0rW1tSE7Oxt2ux3r16+H2RwZt0Mx+MLMwYMH8aMf/QhRUVFYt24dbr755iv+m+qGdqwtr8Oe2hYA6HOBds87xszJMcjPSEJKgj7rcIjoSoE8zDY0NGDRokVYvHgx1qxZM+SDcjhh8IUhn8+HP/7xj3jmmWewatUqPPPMMxg5ciQAeadKiUifAn2Yrampwfe+9z385Cc/waOPPqrRrLXD4AtjX331FR5//HF8+OGHePnll9E1bgZ+taMGXZ7ht0jqvoFB38WoRJEs0IfZ+6ZYsfaRXKxZswarVq1Sb6I6wuCLABUVFXjwyV/B9d2HIBoD38OPMhuxaXW6bq8fIopU3aEX2MOs6HEjd/IIrPnh4iDOTN94c0sEmDNnDu548HmIRmntkVxeH9aV1yk8KyKSo7qhHYVljoBCDwAEswV/azDhaGN7kGamfwy+CNDa6UbFqVYA0l5eiyKw+2QL2jrdyk6MiCRbW14Hl9cn6bOR/jDL4IsAJVWNsscQAJQckj8OEcnX2unGntqWQd/pDSbSH2YZfBHA0ezsc8pLCpfXD0dTh0IzIiI5+DArj7SXPhRSnC6vQuOwJRKRHuj5YTYULsZg8EUAm1WZ32abNTJudSDSOz0+zA5eS9iMl3fV6uZiDAZfBLDH2WAxNct6QjSIXozoaoHf74fBwB1yIi3p7WF2qFrCng407584i4raVs0vxuB3sAiQkxYvewxBMGDb73+J8ePH46mnnsKJEycUmBkRSRFr8cIIeVudoseNvds24YUXXkBFRQVcLpekcS7WEg5eQA90H6rp8vhQWFaDosp6ST+fEhh8EWDMSAsyJsVA6lV8ggAs/M5YHK+qxNatW+H1erFw4UKkpqbiN7/5DZqampSdMBFdoampCa+88gpuu+02/OYn90Lu3SMWqxWPLLkN586dw2OPPYYxY8Zg9uzZKCgowM6dO+F0OoccQ2otYZfHj8Iyh2a1hLy5JUJUN7Qjd30lujyB1/30d3OLz+dDeXk5ioqKUFpaipkzZyIvLw9Lly7tvReUiORpa2vDli1bsHHjRhw+fBh33XUXVqxYgfnz5+Ph4mp8UHNWUkmDIACLpsbitbwZvf9fZ2cnKisrUVFRgYqKChw8eBB2ux2zZ8/uvgTjjjsQExPTZ5zVbx1UdA5qYfBFECnXGw3nrs5vv/0WW7duRVFREfbu3Yvs7Gzk5eVh4cKFMJnkv4sIhVNiRErp6OjAu+++i+LiYuzduxeLFi1Cbm4usrKyYLVae/87pR9mL+d2u3Hw4EHs3bsXFRUV+Pjjj3HDDTf0BuEtM27H8iKHrLMDFpMBHz8xT/W/xwy+CBPs7gwtLS3YvHkzioqKcObMGeTm5iIvLw9paWkBtz1h+ySKFF1dXSgrK0NxcTHef/99zJkzB7m5ubjrrrsQHR094OeC9TDbH5/Ph6NHj/auCD9pvxojUpcCphEBjXMpq8mARxZOwoNzJkgeQwoGXwQ62tiOdeV12H2yBQIunrgCLgbK3MkxyM9MknUx9alTp/D222+jqKgIJpMJeXl5uP/++zF+/PghP8v2SRTuPB4PPvjgAxQXF2Pr1q1IS0tDbm4u7r77blx77bXDHkervys/33QYpUe+lD3O0mnj8PK902SPEwgGXwRr63Sj5FAjHE0dcLo8sFnNsI+NRk6qsluIoihi//79KCoqwqZNmzB58mTk5eXhnnvu6fcvuJpPsUSXCva2us/nw969e7Fx40Zs2bIFEydOxIoVK7B8+XLExcVJHleth9lLPfCXT/GR4yvZ48y3X4/XV92qwIyGj8FHqrpw4QLee+89FBUVYefOnZg3bx7y8vKQnZ0Nq9Ua9PcWRP0J5ra6KIo4cOAAiouLsXnzZlx//fXIzc3Fvffei8TERAV/Feo9zAJc8RFJcu7cOWzZsgVFRUU4fPgwli1bhrapy3C4Zeh6oP5oeUqMQlcwtgpFUcRnn32G4uJiFBcXw2w2Y8WKFcjNzYXdblf2F6CR1/acxsu7amUdbuE7PopoDQ0NeP3tzXjz6wmAhGa5PbQ6JUahSelt9VOnTqG4uBgbN27E+fPnkZubixUrViAlJSXgw11619rpxnfXfBSSpzpZwE66kJCQgLjblsBikfcXIJJvnKfAKFV83dDQgJdeegkzZszA7Nmz0dLSgj/96U84c+YM1qxZg2nTpoVd6AHKXIwxd3KMJg+pvKuTdEPPN85T+JHVyNXjw/98/X24dv0eJ06cwN133401a9YgIyNDkdrVUPFwZhL2nmqV9E7eajIiPzMpCLMaWuT8DpHu6fHGeQpPshu5AqjrsuJXv3gCy7L/G0aMkF7LFspSEkajIMsucbvYrtlBNG51km4odeP8547j2Lp1K86ePavIeBR+lGjkajGb4bxuSsSGXo+89EQUZE1BlNk45LanIHSfvta69IgrPtINJdonmQUR1wjf4j//cyMOHDiAUaNGYebMmZg1axZmzpyJtLQ0XHXVVQrOevh49Zp+cFtdWXnpiUiOH616LaFUPNVJuqH0KTFRFHHq1CkcOHAA+/fvx4EDB3Ds2DFMnDixNwhnzZqFKVOmwGg0Kvgr6YtXr+lPKBdf652atYRSMfhIV4J92zdd+B4AAAmMSURBVLvb7UZ1dXVvEO7fvx/Nzc1IS0vrszKMj5ffwxDg1Wt69dONVfjb0WbZ42hRfE3yMfhIV7S4ueXrr7/Gp59+2huE+/fvx4gRIzBz5szeMJwxYwZsNltA4/LqNX05e/YsysrKsH37duxttSBq5nJZNaNaFV+TfAw+0h2tA0MURdTX1/fZIj1y5AhuvPHGPqvCW265BWZz/984efWa9vx+Pw4fPoxt27Zh+/btqK2txcKFC5GdnY1ZGQuw5I3PQrL4muRj8JEu6W2L0OPx4NixY322SOvr6zFt2rTeIJw5cybGjx8PQRBCtkFnqOvo6MCuXbuwbds2lJWVYdSoUcjOzkZ2djbuuOOOPicw+XsUuRh8pFta3DgfCKfTiYMHD/bZIvV4PEi9LQN131kFn4xqIa4mhq+urg7bt2/Htm3bUFlZidtuu6037JKSBi6Q5qo8cjH4SPdC4ZRYj8bGRvz63SrsaBDgE6SfFA2V90dalGhcuHAB+/bt6w07p9PZG3QLFiwYtHHr5bTeVidtsI6PdO+6kRbdB0CP+Ph4mMa0wNcor12L3mvEBi/RaMbLu2oVLdE4e/YsduzYgW3btmHXrl2YNGkSFi9ejA0bNmD69OkwGKStrnvCS0/b6hR8XPERKUypGrFx+Bqr7X5MnDgRkyZNwnXXXaeLy47VeP/aczClZ1VXW1uLBQsWYPHixbjzzjsRGxsr7xdxGb1vq5OyGHxEClOqQecEYxvG1r+PU6dOoba2FgAwadKk3iC89N+BllpIFcytwZ6DKdu3b8f27dths9mQnZ2NxYsXX3EwJVhCaVudpGPwESksGA06RVFEW1sbamtre4Ow53+fOnUK0dHR/YbihAkTEBUVpcivKxiHQXoOpmzfvh2ffPIJ0tPTe9/XTZw4UZF5E12OwUekMLUbdIqiiC+//LLfUDxz5gxiY2P7DcXExMQB6xD7o8Tx/9/dk9x7MGX79u04d+4csrKykJ2djYULFwZ0MIVIKgYfURDopUbM6/Xin//8Z58w7Pn3F198gRtvvLE3CC8Nxfj4+D4HRpQIc8Hvxbn/+gmSEuJ6V3WpqamSD6YQScXgIwqCUKgRc7vd+Pzzz/tdKX7zzTeYMGFCbxC2jElBxTc2BNisvA+T4MdDtyfg0cW825K0xXIGoiAIhQadFosFU6ZMwZQpU674sc7OTtTV1fUG4UcNX8NzlbwDNF7RgC/OyxqCSBEMPqIgCeUasZEjR2LatGmYNq17dXZaoRINp8sjewwiuRh8REEUag06B2KzKvOtwmaV3g2BSCkMPqIgS44fjdfyZoR0jZg9zgaLqVl2iYZ9LE9tkvZ4uIWIhqR2iQZRMPEcMRENacxICzImxUDqjWmC0L2dy9AjPWDwEdGwPJyZBKtJWscJq8mI/MyBWwQRqYnBR0TD0lOiEWUO7NuGmiUaRMPBwy1ENGyhXKJB1IOHW4goYD1tfN479iWMBgFe8eIqMFRKNChyMfiISBJRFJGQNAU/evEvaPVaQq5EgyIXtzqJSJK6ujrA3Ymn7p6piwa5RMPFwy1EJMlHH32EefPmMfQo5DD4iEiSnuAjCjV8x0dEARNFEbGxsfj0009x4403aj0dooBwxUdEATt+/Diio6MZehSSGHxEFDBuc1Io46lOIhpSa6cbJVWNcDQ74XR5cdjhx6ybF6Gt082yBQo5fMdHRAOqbmjH2vI67KltAYA+3RksRgEQBGROjkF+RhJSElioTqGBwUdE/SqqrOfVZBSWuNVJRFfoDr0adHmG7r8nikCXx4fCshoAYPiR7vFwCxH1Ud3QjsIyx7BC71JdHj8Kyxw42tgepJkRKYPBR0R9rC2vg8vrk/RZl9eHdeV1Cs+ISFkMPiLq1drpxp7alkHf6Q1GFIHdJ1vQ1ulWdmJECmLwEVGvkqpG2WMIAEoOyR+HKFgYfETUy9Hs7FOyIIXL64ejqUOhGREpj8FHRL2cLq9C43gUGYcoGFjOQBRCLr9BxWY1wR5nw/I0ZRq/2qzKfEuwWc2KjEMUDAw+ohAw2A0qVlMzXt5VK+sGFafTiZ07d+Lw7pMQo78DwSQ9RK0mA+xjoyV/nijYeHMLkc4F6waVpqYm/O1vf0NpaSn+/ve/44477sDCf78br35xAy74pH9bsJgM+PiJebzDk3SLKz4iHVP6BpWTJ0+itLQUpaWlcDgcyMrKwgMPPIDNmzcjOrp7lVbz1kF8UHNWUkmDIABzJ8cw9EjXuOIj0qnqhnbkrq9ElyfwYvIosxGbVqfj5hts+PTTT3vDzul0YsmSJViyZAkyMjIwYsSIoPy8yfG8sJr0i8FHpFOr5ay8AFx/4Us0bHwG11xzTW/YpaWlwWAY+jB3ICvNHlFmAwqypvCuTtI9bnUS6ZDsG1QAtI6IQ+nOXZiZPDXgz/eE13DeLYp+P6JGmBh6FDJYx0ekQ0rcoGI2mXC4Xfq7trz0RGxanY5FU2NhMRlgNfX9dmE1GWAxGRDnPYtbvtrF0KOQwRUfkQ7p5QaV5PjReC1vBto63Sg51AhHUwecLg9sVjPsY6ORkxqPEeIFJCcno6xsPrKysgAEv96QSA4GH5EO6e0GletGWvDgnAkD/KgFr7/+OlatWoWN73+M/zp4Nmj1hkRKYPAR6VCo3aAyb948JC/7MVa+eQii0dTvO0HXv0Lw/RNnUVHbyo7tpBm+4yPSIXucDRaTvL+eat6gUlRZj9PRyfAb+g+9S11ab1hUWa/K/IguxeAj0qGctHjZY4gAclLljzOUno7trgDfSbJjO2mFwUekQ2NGWpAxKQaCIO3zat6gwo7tFGoYfEQ69XBmEqwmo6TPWk1G5GcmKTyjK7FjO4UiBh+RTqUkjEZBlh1R5sD+mnbfoGJX5dowdmynUMRTnUQ6FsgNKoF2Z1CCXuoNiQLB4CPSubz0RCTHj8a68jrsPtkCAehzkMRqMkBE9zu9/MwkVS+I1lu9IdFwMPiIQsBwblDR4kaUUKs3JAIYfEQhZfAbVNTXXW/YLGu7kx3bSW083EJEkoVSvSFRDwYfEUkWSvWGRD0YfEQkSyjUGxJdisFHRLKEQr0h0aV4uIWIZNN7vSHRpQRRlHrZEBFRX0cb23VZb0h0KQYfESlOb/WGRJdi8BERUUTh4RYiIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIoooDD4iIooo/x/cgFaDqouxqgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "smile = df.smiles.iloc[0]\n",
    "graph = read_smiles(smile, explicit_hydrogen=False, zero_order_bonds=True, \n",
    "                    reinterpret_aromatic=True)\n",
    "print(smile)\n",
    "nx.draw(graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From Graph to diagrams with HKS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the graphs we compute persistence diagrams. The time parameter in the HKS algorithm can be consider as a kind of **scale** parameter!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "hkh = GeneralFiltrationGraphHomology(time=.1, n_jobs=-1)\n",
    "diagrams_hkh_01 = hkh.fit_transform(df.smiles_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "hkh = GeneralFiltrationGraphHomology(time=10, n_jobs=-1)\n",
    "diagrams_hkh_10 = hkh.fit_transform(df.smiles_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "hkh = GeneralFiltrationGraphHomology(time=1, n_jobs=-1)\n",
    "diagrams_hkh_1 = hkh.fit_transform(df.smiles_graph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagrams_hkh_10[:, :, -1] = diagrams_hkh_10[:, :, -1] + 4\n",
    "diagrams_hkh = np.concatenate([diagrams_hkh_01, diagrams_hkh_10], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the HKS agorithm gives as output 4 subdiagrams that for seek of simplicity we encoded in the giotto-framework in the same way as different homology dimensions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1484"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.07547169811320754"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.CT_TOX.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset is highly unbalanced, with 92.5% of the molcules that are safe drugs and only 7.5% of them labeled as toxic. Thus, the accuracy is not a good metric for this kind of tasks, and for this reason we chose to use the AUC-ROC metric as measure of the performance of our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXhTdfbH8fehUhEFEXAZ1qIikgJlKUtFsAoKDv5cBnQYRkeGlpZNZXMAFUSWQRAEZa+tLIIssgoUChTKWpayQxFlEIcCKntZLW2/vz9SnFoClCbpbZLzep485Cb3Jp9cbU/PvcmJGGNQSimlcipkdQCllFIFkxYIpZRSDmmBUEop5ZAWCKWUUg5pgVBKKeXQHVYHcKXSpUubgIAAq2MopZTH2LZt20ljzP2O7vOqAhEQEEBSUpLVMZRSymOIyE83uk8PMSmllHJIC4RSSimHtEAopZRySAuEUkoph7RAKKWUckgLhFJKKYe0QACJiTBkCCT2WgDNmkFUlNWRlFLKcm4tECLSXEQOiMhBEent4P5QETknIjuzLv2ybi8vIqtFZL+I7BORd9yVMTERmjSBvu9n0mTYcyQuT4XISC0SSimPsH79eoYNG+aWx3ZbgRARP2As8DxgA/4mIjYHq64zxtTMugzIui0d6GGMqQo0ADrfYFunJSRAWhpkmEKkUZgEQu13zJ3rjqdTSimXOH/+PF26dKFRo0ZMnDiRixcvuvw53NlB1AMOGmMOGWPSgJnAS7nZ0Bhz3BizPev6eWA/UNYdIUNDwd8f/Apl4s9VQkmw39GypTueTimlnBYXF0e1atUYN24c77zzDrt27eLuu+92+fO4c9RGWeBItuUUoL6D9UJEZBdwDOhpjNmX/U4RCQBqAZvdETIkBOLjISGhEKFnlxOyszi0nAgREe54OqWUyrNTp07RvXt3pk6dStWqVdmwYQMhISFuez53FghxcFvO7zfdDlQ0xlwQkT8DC4DKvz+AyD3AXKCrMSbV4ZOIRAARABUqVMhT0JAQ+wVezroopVTBYYxh7ty5dO7cmdOnT/PBBx/wwQcfcOedd7r1ed15iCkFKJ9tuRz2LuF3xphUY8yFrOuxQGERKQ0gIoWxF4fpxph5N3oSY0yUMSbYGBN8//0OBxIqpZTHOn78OC1btuTVV1+lfPnyJCUlMXDgQLcXB3BvgdgKVBaRSiLiD7QGvs2+gog8JCKSdb1eVp5TWbfFAPuNMZ+6MaNSShVIxhgmTZqEzWZj6dKlDB06lE2bNhEUFJRvGdx2iMkYky4iXYA4wA/40hizT0Q6ZN0/AWgFdBSRdOAy0NoYY0TkSeANYI+I7Mx6yPeyugyllPJqP/74IxEREaxcuZJGjRoRHR3NY489lu85xJicpwU8V3BwsNHvg1BKeaqMjAzGjBnDe++9h5+fH0OHDiUyMpJChdx3sEdEthljgh3d51VfGKSUUp4qOTmZ8PBwEhMTef7555k4cSLly5e/9YZupKM2lFLKQlevXmXQoEHUqlWL77//nmnTprFkyRLLiwNogfhd4pFEhqwbQuKRRKujuJw3v7Y/iIrSWVrKo2zbto3g4GD69u3LK6+8QnJyMn//+9/Jeu+O5fQQE/ZfoE2mNiEtIw1/P3/i/xFPSHn3ffgkP3nza/uDqCj7DC2A5cvt/+qHHVUBdfnyZfr378/w4cN58MEHWbBgAS+9lKtBE/lKOwgg4XACaRlpZJgM0jLSSDicYHUkl/Hm1/YHOWdn6SwtVUCtXbuWoKAghg0bRlhYGMnJyQWyOIAWCABCA0Lx9/PHT/zw9/MnNCDU6kgu482v7Q9yzs7SWVqqgElNTaVTp0489dRTZGRkEB8fT1RUFCVKlLA62g3p21yzJB5JJOFwAqEBoV53CMabX9sfREXZO4eWLfXwkipQYmNjiYyM5NixY3Tt2pUBAwa4ZbheXtzsba5aIJRSyk1OnjxJ165dmT59OjabjS+//JL69R3NLLXOzQqEHmJSSikXM8Ywa9YsbDYbs2bN4sMPP2T79u0Frjjcir6LSSmlXOjYsWN07NiRb7/9luDgYOLj46levbrVsfJEOwillHIBYwzR0dHYbDZWrFjB8OHDSUxM9NjiANpBKKWU0/7zn/8QERHBqlWrCA0N5YsvvuDRRx+1OpbTtINQSqk8ysjI4NNPP6V69eokJSUxceJE4uPjvaI4gHYQSimVJ3v37iUsLIwtW7bwwgsvMH78eMqVK2d1LJfSDkIppW5DWloaH330EbVr1+bQoUN8/fXXfPvtt15XHEA7CKWUyrWtW7fSrl079u7dS5s2bRg1ahTe/FXH2kEopdQtXLp0iZ49e9KgQQPOnDnDokWLmD59ulcXB9AOQimlbmr16tWEh4dz6NAhIiMjGTp0KPfee6/VsfKFdhBKKeXAuXPniIyM5JlnnkFEWL16NRMmTPCZ4gBaIJRS6jqLFi3CZrMRHR1Nz5492b17N6GhoVbHyndaIJRSKsuJEydo06YNL774IqVKlWLTpk188sknFC1a1OpoltACoZTyecYYvv76a6pWrcqcOXMYMGAASUlJ1K1b1+poltKT1Eopn5aSkkLHjh1ZvHgx9evXJyYmhsDAQKtjFQjaQSilfFJmZiYTJ07EZrOxatUqRo4cyYYNG7Q4ZKMdhFLK5/zwww+0b9+eNWvW0KRJE6Kionj44YetjlXgaAehlPIZ6enpDB8+nBo1arBz506io6NZsWKFFocb0A5CKeUTdu/eTVhYGElJSbz00kuMGzeOMmXKWB2rQNMOQinl1X777Tf69etHnTp1+Omnn5g1axbz58/X4pAL2kEopbzWpk2bCAsLIzk5mTfeeIORI0dSqlQpq2N5DO0glFJe5+LFi3Tr1o0nnniC8+fPExsby9SpU7U43CbtIJRSXiU+Pp727dvz448/0qlTJ4YMGULx4sWtjuWRtINQSnmFs2fPEh4eTtOmTbnjjjtYs2YNY8eO1eLgBC0QSimPt3DhQmw2G5MnT6ZXr17s2rWLxo0bWx3L4+khJqWUx/rll194++23mT17NkFBQSxatIg6depYHctraAehlPI4xhi++uorbDYbCxYsYPDgwWzdulWLg4tpB6GU8ij//e9/6dChA0uXLiUkJISYmBiqVq1qdSyv5NYOQkSai8gBETkoIr0d3B8qIudEZGfWpV9ut1VK+ZbMzEzGjRtHYGAga9eu5fPPP2fdunVaHNzIbR2EiPgBY4FngRRgq4h8a4xJzrHqOmPMC3ncVinlA77//nvCw8NZt24dzz77LFFRUQQEBFgdy+u5s4OoBxw0xhwyxqQBM4GX8mFbpZSXSE9PZ+jQodSoUYM9e/YwadIk4uLitDjkE3cWiLLAkWzLKVm35RQiIrtEZKmIXBvEntttEZEIEUkSkaQTJ064IrdSqgDYuXMn9evXp3fv3rRo0YLk5GTatm2LiFgdzWe4s0A4+q9ocixvByoaY4KA0cCC29jWfqMxUcaYYGNM8P3335+3pImJ0LGj/ZKYmLfH8ABRI1+nWddSRI183eooSt3QlStXeP/99wkODubo0aPMmTOHuXPn8qc//cnqaD7Hne9iSgHKZ1suBxzLvoIxJjXb9VgRGScipXOzrcskJkJoKKSl2ZcnTYLVqyEkxC1PZ5Woka8TeW46lIDl56bDSIjoNs3qWEr9wcaNGwkLC+O7777jzTff5NNPP6VkyZJWx/JZ7uwgtgKVRaSSiPgDrYFvs68gIg9JVr8oIvWy8pzKzbYuk5AAV6/+bzktzX6bl5n701L7FcmxrFQBcOHCBd5++22efPJJLl26xLJly5g8ebIWB4u5rUAYY9KBLkAcsB+YbYzZJyIdRKRD1mqtgL0isgv4HGht7Bxu65agoaFQuPD/lv397bd5mZYVn7dfMTmWlbLY8uXLqVatGmPGjKFz587s3buXZs2aWR1LAWKMw0P7Hik4ONgkJSXd/oaJiTB1qv36P/7hdYeXroka+Tpzf1pKy4rP6+ElZbnTp0/To0cPJk+eTJUqVYiJiaFhw4ZWx/I5IrLNGBPs8D4tEEqp/DZ37lw6d+7MyZMn6dWrF3379qVIkSJWx/JJNysQOmpDKZVvfv75Z7p06cLcuXOpVasWy5Yto2bNmlbHUjegw/qUUm5njGHy5MnYbDYWL17MkCFD2Lx5sxaHAk47CKWUWx0+fJjIyEiWL1/Ok08+SXR0NFWqVLE6lsoF7SCUUm6RmZnJ6NGjqVatGhs3bmTMmDGsWbNGi4MH0Q5CKeVy3333HeHh4WzYsIHmzZszYcIEKlasaHUsdZu0g1BKuczVq1f597//TVBQEPv372fq1KnExsZqcfBQ2kEopVxi+/bthIWFsXPnTl599VVGjx7Ngw8+aHUs5QTtIJRSTrl8+TJ9+vShXr16/Pzzz8ybN4/Zs2drcfAC2kEopfJs/fr1hIWF8f3339OuXTuGDx/OfffdZ3Us5SLaQSilbtv58+fp0qULjRo1Ii0tjRUrVhATE6PFwctogVBK3ZalS5cSGBjIuHHj6Nq1K3v27KFp06ZWx1JuoIeYlFK5curUKbp168ZXX31F1apV2bBhAyFeOthS2WkHoZS6KWMM33zzDTabjRkzZtC3b1927NihxcEHaAehlLqh48eP06lTJxYsWECdOnVYvnw5QUFBVsdS+UQ7CKXUdYwxfPnll1StWpVly5YxbNgwNm3apMXBx2gHoZT6g0OHDhEZGcnKlStp3LgxX3zxBY899pjVsZQFtINQSgGQkZHBqFGjqF69Ops3b2b8+PGsXr1ai4MP0w5CKUVycjJhYWFs2rSJP//5z0yYMIHy5ctbHUtZTDsIpXxYWloaAwcOpFatWvzwww9MmzaNxYsXa3FQgHYQSvmspKQkwsLC2L17N61bt+azzz7jgQcesDqWKkC0g1DKx1y6dIl//etf1K9fn5MnT7Jw4UJmzJihxUFdRzsIpXzImjVrCA8P5+DBg7Rv355hw4ZRokQJq2OpAko7CKV8QGpqKh07diQ0NJTMzEzi4+OJiorS4qBuSguEUl5uyZIlBAYGEhUVRffu3dmzZw/PPPOM1bGUB9ACoZSXOnnyJK+//jovvPAC9957Lxs3bmTEiBEULVrU6mjKQ2iBUMrLGGOYOXMmVatWZfbs2Xz44Yds376d+vXrWx1NeRg9Sa2UFzl69CidOnXi22+/pW7dusTExFC9enWrYykPpR2EUl7AGMMXX3yBzWZjxYoVjBgxgsTERC0OyinaQSjl4f7zn//Qvn17Vq9ezdNPP80XX3zBI488YnUs5QW0g1DKQ2VkZPDpp59SvXp1tm3bRlRUFPHx8VoclMtoB6GUB9q7dy9hYWFs2bKF//u//2P8+PGULVvW6ljKy2gHoZQHSUtLo3///tSuXZtDhw4xY8YMFi5cqMVBuYV2EEp5iC1bttCuXTv27dtHmzZt+OyzzyhdurTVsZQX0w5CqQLu0qVL9OjRg5CQEM6ePcuiRYuYPn26FgfldrkqECLyFxH5QUTOiUiqiJwXkVR3h1PK161evZrq1avz6aefEhERQXJyMi+88ILVsZSPyG0HMQx40RhzrzGmuDGmmDGm+K02EpHmInJARA6KSO+brFdXRDJEpFW227qJyD4R2SsiM0SkSC6zKuXxzp07R0REBM888wyFChUiISGB8ePHU7z4LX/slHKZ3BaIX4wx+2/ngUXEDxgLPA/YgL+JiO0G6w0F4rLdVhZ4Gwg2xlQD/IDWt/P8SnmqRYsWYbPZiImJ4d1332XXrl089dRTVsdSPuimJ6lF5C9ZV5NEZBawAPjt2v3GmHk32bwecNAYcyjrsWYCLwHJOdZ7C5gL1HWQ7S4RuQoUBY7d/KUo5dl+/fVX3nnnHWbOnEn16tVZuHAhwcHBVsdSPuxW72L6v2zXLwHPZVs2wM0KRFngSLblFOAP08KyOoVXgGfIViCMMUdFZDjwX+AysNwYs9zRk4hIBBABUKFChVu8HKUKHmMMX3/9Ne+88w6pqakMGDCAXr164e/vb3U05eNuWiCMMf8EEJGGxpgN2e8TkYa3eGxx9JA5lkcBvYwxGSL/W11E7sPebVQCzgLfiMjrxphpDjJGAVEAwcHBOR9fqQLtyJEjdOzYkSVLltCgQQOio6MJDAy0OpZSQO4/BzEaqJ2L27JLAcpnWy7H9YeJgoGZWcWhNPBnEUkHCgM/GmNOAIjIPOAJ4LoCoZQnyszMJCoqin/9619kZGQwatQounTpgp+fn9XRlPrdrc5BhGD/xXy/iHTPdldx7CeOb2YrUFlEKgFHsZ9kbpN9BWNMpWzPNRlYbIxZICL1gQYiUhT7IaYmQFKuXpFSBdwPP/xA+/btWbNmDU2aNCEqKoqHH37Y6lhKXedW72LyB+7BXkiKZbukAq1ush3GmHSgC/Z3J+0HZhtj9olIBxHpcIttNwNzgO3AnqycUbd8NUoVYOnp6XzyySfUqFGDnTt3EhMTw4oVK7Q4qAJLjLn1YXsRqWiM+Skf8jglODjYJCVpo6EKnl27dhEWFsa2bdt46aWXGDduHGXKlLE6llKIyDZjjMO3y+X2HMQlEfkECAR+/8CaMUa/+Vypm/jtt98YNGgQH3/8MSVLlmT27Nm0atWK7G/KUKqgyu0H5aYD32F/V9FHwGHs5xiUUjeQmJhIrVq1GDRoEG3atCE5OZlXX31Vi4PyGLktEKWMMTHAVWPMGmNMO6CBG3Mp5bEuXrxI165dadiwIRcuXCA2NpYpU6ZQqlQpq6MpdVtye4jpata/x0WkBfa3q5ZzTySlPNfKlStp3749hw8fpnPnzgwZMoRixYpZHUupPMltgRgkIvcCPbB//qE40M1tqZTyMGfOnKFnz558+eWXVK5cmbVr19KoUSOrYynllFwVCGPM4qyr54Cn3RdHKc8zf/58OnXqxIkTJ+jduzf9+vXjrrvusjqWUk7L7fdBPCYi8SKyN2u5hoh84N5oShVsv/zyC6+99hp/+ctfeOihh9iyZQtDhgzR4qC8Rm5PUn8B9CHrXIQxZjc6flv5KGMMU6dOpWrVqixcuJDBgwezZcsWate+2eQZpTxPbs9BFDXGbMnx9rx0N+RRqkD773//S2RkJMuWLeOJJ54gJiaGxx9/3OpYSrlFbjuIkyLyCFnTWLO++e2421IpVcBkZmYyduxYAgMDWbduHZ9//jnr1q3T4qC8Wm47iM7YZyE9LiJHgR+Bv7stlVIFyIEDBwgPD2f9+vU8++yzREVFERAQYHUspdzuVtNcs09wjQVWY+86LgItgU/dF00pa129epURI0bQv39/7rrrLiZNmsSbb76pn4RWPuNWHcS1T/hUwf6NbwuxfxHQG8BaN+ZSylI7duwgLCyMHTt20LJlS8aMGcNDDz1kdSyl8tWtvlHuIwARWQ7UNsacz1ruD3zj9nRK5bMrV64wcOBAhg4dSunSpZkzZw4tW7a0OpZSlsjtOYgKQFq25TQgwOVplLLQhg0bCAsL48CBA7Rt25YRI0ZQsmRJq2MpZZncFoivgC0iMh/7O5leAaa4LZVS+ejChQu89957jBkzhgoVKhAXF8dzzz1ndSylLJfbURuDRWQpcG24zD+NMTvcF0up/BEXF0dERARHjhyhS5cu/Pvf/+aee+6xOpZSBUJuOwiMMduxfwWoUh7v9OnTdO/enSlTpvD444+zbt06GjZsaHUspQqU3H5QTimvMXfuXGw2G9OmTeP9999nx44dWhyUciDXHYRSnu748eN06dKFefPmUatWLZYtW0bNmjWtjqVUgaUdhPJ6xhgmT56MzWZjyZIlfPzxx2zZskWLg1K3oB2E8mqHDx8mIiKCFStW8OSTTxIdHU2VKlWsjqWUR9AOQnmljIwMPv/8c6pVq0ZiYiJjx45lzZo1WhyUug3aQSivs3//fsLDw9m4cSPNmzdn4sSJVKhQwepYSnkc7SB8SOLSKIYMakbi0iiro7jF1atXGTx4MDVr1uS7775j6tSpxMbGanFQKo+0g/ARiUujaLIhkjQ/8N+wnHgg5PkIq2O5zPbt22nXrh27du3itdde4/PPP+fBBx+0OpZSHk07CB+RsG0uaX6QUQjSCtmXvcHly5fp3bs39erV45dffmH+/PnMmjVLi4NSLqAFwkeE1mmJfwb4ZYB/pn3Z061du5agoCCGDh1K27ZtSU5O5uWXX7Y6llJeQw8x+YiQ5yOIx945hNZp6dGHl1JTU+nTpw/jxo0jICCAFStW0LRpU6tjKeV1xBhjdQaXCQ4ONklJSVbHUG60dOlSIiMjSUlJ4Z133mHQoEHcfffdVsdSymOJyDZjTLCj+7SDUB7h1KlTdOvWja+++gqbzcbGjRtp0KCB1bGU8mp6DkIVaMYYZs+eTdWqVZkxYwZ9+/Zl+/btWhyUygfaQagC69ixY3Tq1ImFCxdSp04dVq5cSY0aNayOpZTP0A5CFTjGGGJiYrDZbMTFxTFs2DA2bdqkxUGpfKYdhCpQDh06RPv27Vm1ahWNGzcmOjqaypUrWx1LKZ+kHYQqEDIyMhg1ahTVq1dn69atTJgwgdWrV2txUMpCbi0QItJcRA6IyEER6X2T9eqKSIaItMp2WwkRmSMi34nIfhEJcWdWZZ19+/bRsGFDunXrxtNPP01ycjKRkZEUKqR/vyhlJbf9BIqIHzAWeB6wAX8TEdsN1hsKxOW46zNgmTHmcSAI2O+urMoaaWlpDBw4kFq1anHw4EGmT5/OokWLKFeunNXRlFK49xxEPeCgMeYQgIjMBF4CknOs9xYwF6h77QYRKQ40BtoCGGPSgDQ3ZlX5bOvWrYSFhbFnzx5at27N559/zv333291LKVUNu7s4csCR7Itp2Td9jsRKQu8AkzIse3DwAlgkojsEJFoEXH4cVkRiRCRJBFJOnHihOvSK7e4dOkS7777Lg0aNODUqVMsXLiQGTNmaHFQqgByZ4EQB7flnOsxCuhljMnIcfsdQG1gvDGmFnARcHgOwxgTZYwJNsYE6y+Zgi0hIYGgoCCGDx9OeHg4ycnJvPjii1bHUkrdgDsPMaUA5bMtlwOO5VgnGJgpIgClgT+LSDqwCUgxxmzOWm8ONygQquA7d+4cvXr1YuLEiTzyyCOsWrWKp59+2upYSqlbcGeB2ApUFpFKwFGgNdAm+wrGmErXrovIZGCxMWZB1vIREalijDkANOH6cxfKAyxZsoTIyEiOHz9Ojx49GDBgAEWLFrU6llJcvXqVlJQUrly5YnWUfFGkSBHKlStH4cKFc72N2wqEMSZdRLpgf3eSH/ClMWafiHTIuj/neYec3gKmi4g/cAj4p7uyKtc7ceIEXbt25euvv6ZatWrMmzePevXqWR1Lqd+lpKRQrFgxAgICyDqK4bWMMZw6dYqUlBQqVap06w2yuPWT1MaYWCA2x20OC4Mxpm2O5Z3YD0EpD2KMYebMmbz99tucO3eO/v3706dPH/z9/a2OptQfXLlyxSeKA4CIUKpUKW73jTw6akO5TEpKCh07dmTx4sXUq1ePmJgYqlWrZnUspW7IF4rDNXl5rfpRVeW0zMxMoqKiCAwMJD4+nhEjRrBx40YtDkp5OC0QyikHDx6kSZMmREZGUqdOHfbs2UP37t3x8/OzOppSHqd///4MHz6cfv36sXLlytvadsqUKVSuXJnKlSszZcoUl+TRQ0y+IDEREhIgNBRCXDPS6tpwvb59+1K4cGG++OILwsLCfKplVz7IDT9LjgwYMOC21j99+jQfffQRSUlJiAh16tThxRdf5L777nMqh3YQ3i4xEZo0gb597f8mJjr9kHv27CEkJISePXvStGlTkpOTCQ8P1+KgvJsbfpYABg8eTJUqVWjatCkHDhwAoG3btsyZMweAgIAA3nvvPUJCQggODmb79u00a9aMRx55hAkT7O/5iYuL49lnn6VkyZLcd999PPvssyxbtszpbFogvF1CAqSlQUaG/d+EhDw/1G+//caHH35I7dq1OXz4MDNnzmThwoWULVv21hsr5elc+LN0zbZt25g5cyY7duxg3rx5bN261eF65cuXJzExkUaNGv1ePDZt2kS/fv0AOHr0KOXL/+9zyeXKlePo0aNO59NDTN4uNBT8/e3/Q/v725fzYPPmzYSFhbFv3z7+/ve/M2rUKEqXLu3SqEoVaC76Wcpu3bp1vPLKK79/ePRGo2eu3V69enUuXLhAsWLFKFasGEWKFOHs2bMYk3OKkWveoaUdhLcLCYH4eBg40P7vbR43vXjxIt27dyckJIRz586xePFipk2bpsVB+R4nf5ZuJDe/yO+8804AChUq9Pv1a8vp6emUK1eOI0f+Nxs1JSWFMmXKOJ1NC4QvCAmBPn1u+3/oVatWUaNGDUaOHEmHDh3Yt28fLVq0cFNIpTxAHn+WbqRx48bMnz+fy5cvc/78eRYtWpSnx2nWrBnLly/nzJkznDlzhuXLl9OsWTOn8+khJnWds2fP8u677xIdHc2jjz5KQkICTz31lNWxlPI6tWvX5q9//Ss1a9akYsWKNGrUKE+PU7JkSfr27Uvduvav1enXrx8lS5Z0Op84OnblqYKDg01SUpLVMTzawoUL6dixI7/88gs9e/akf//+3HXXXVbHUsrl9u/fT9WqVa2Oka8cvWYR2WaMcTjWSA8xKQB+/fVXWrduzcsvv0zp0qXZvHkzQ4cO1eKglA/TAuHjjDFMmzaNqlWrMn/+fAYOHEhSUhLBwTonUSlfp+cgfNiRI0fo0KEDsbGxNGjQgJiYGGw2m9WxlFIFhHYQPigzM5Px48cTGBhIQkICo0aNYv369VoclFJ/oB2Ej/n+++9p3749a9eupWnTpkRFRd3WF4gopXyHdhA+Ij09nWHDhhEUFMSuXbuIiYlh+fLlWhyUUjekHYQP2LVrF+3atWP79u28/PLLjB071iWfslRKuVb//v255557SE1NpXHjxjRt2jTX2zZv3pxNmzbx5JNPsnjxYpfk0QLhxX777TcGDRrExx9/TMmSJfnmm29o2bKlTl1VKo/yadr3bY/7Bnj33Xe5dOkSEydOdFkOPcTkpRITE6lVqxaDBg2iTZs2JCcn06pVKy0OSncawv0AAA39SURBVOWRm6Z9u2TcN0CTJk0oVqyYa0Jl0QLhZS5cuEDXrl1p2LAhFy9eZOnSpUyZMoVSpUpZHU0pj+aGad8uG/ftLnqIyYusWLGCiIgIDh8+TOfOnRkyZIjL/6JQyle5Ydq3y8Z9lyhRwvkwDmiB8AJnzpyhR48eTJo0iccee4y1a9fmeeiXUsqxa9O+XX0OwhXjvt1FDzF5uPnz52Oz2Zg6dSq9e/dm165dWhyUchMXT/t22bhvd9EOwkP9/PPPvPXWW8yZM4eaNWuyZMkSateubXUspdRtcNW4b4BGjRrx3XffceHCBcqVK0dMTIzT3wmh4749jDGGr776iq5du3Lp0iU+/PBDevbsSeHCha2OppRH0XHfdjcb960dhAf56aefiIyMJC4ujieeeIKYmBgef/xxq2MppbyUnoPwAJmZmYwZM4bAwEDWr1/P6NGjWbdunRYHpZRbaQdRwB04cICwsDA2bNjAc889x8SJEwkICLA6llLKB2gHUUBdvXqVIUOGEBQURHJyMpMnT2bZsmVaHJRS+UY7iAJox44dhIWFsWPHDlq1asXo0aN56KGHrI6llPIx2kEUIFeuXOG9996jbt26HDt2jLlz5/LNN99ocVBKWUI7iAJi/fr1hIeHc+DAAf75z38yYsQI7rvvPqtjKaXyUV7Hfe/cuZOOHTuSmpqKn58f77//Pn/961+dzqMFwmLnz5+nT58+jB07looVKxIXF8dzzz1ndSyllAOJRxJJOJxAaEAoIeXdN+/7dsd9Fy1alKlTp1K5cmWOHTtGnTp1aNasmdMzmvQQk4Xi4uKoVq0a48aN4+2332bv3r1aHJQqoBKPJNJkahP6ru5Lk6lNSDzimnnfrhj3/dhjj1G5cmUAypQpwwMPPMCJEyeczubWAiEizUXkgIgcFJHeN1mvrohkiEirHLf7icgOEXHN1yMVEKdPn+bNN9+kefPmFC1alHXr1vHZZ59xzz33WB1NKXUDCYcTSMtII8NkkJaRRsLhBKcf0x3jvrds2UJaWhqPPPKI0/ncdohJRPyAscCzQAqwVUS+NcYkO1hvKBDn4GHeAfYDxd2VM7/NmTOHzp07c/r0ad5//30++OADihQpYnUspdQthAaE4u/nT1pGGv5+/oQGhDr9mK4e9338+HHeeOMNpkyZQqFCzv/9785zEPWAg8aYQwAiMhN4CUjOsd5bwFygbvYbRaQc0AIYDHR3Y858cfz4cTp37sz8+fOpXbs2cXFx1KxZ0+pYSqlcCikfQvw/4l1+DsJV475TU1Np0aIFgwYNokGDBi7J5s5DTGWBI9mWU7Ju+52IlAVeASZwvVHAv4BMdwXMD8YYJk2ahM1mIzY2lo8//pjNmzdrcVDKA4WUD6FPoz4uKw6uGvedlpbGK6+8wj/+8Q9effVVl2QD93YQjspiztGxo4BexpiM7FVURF4AfjXGbBOR0Js+iUgEEAFQoUIFpwK72o8//khERAQrV66kUaNGREdH89hjj1kdSylVQLhq3Pfs2bNZu3Ytp06dYvLkyQBMnjzZ6T9E3TbuW0RCgP7GmGZZy30AjDFDsq3zI/8rJKWBS9h/2dcH3gDSgSLYz0HMM8a8frPnLCjjvjMyMhg7dix9+vShUKFCDB06lA4dOrjkmKBSyjV03LedVeO+twKVRaQScBRoDbTJvoIxplK2kJOBxcaYBcACoE/W7aFAz1sVh4Ji//79hIWFkZiYyPPPP8+ECRMKXGejlFK54bY/aY0x6UAX7O9O2g/MNsbsE5EOItLBXc9rlatXrzJ48GBq1qzJgQMH+Oqrr1iyZIkWB6WUx3LrJ6mNMbFAbI7bHJ2QxhjT9ga3JwAJLo7mUtu2baNdu3bs3r2b1157jdGjR/PAAw9YHUsppZyiB8WdcPnyZXr16kW9evU4ceIE8+fPZ9asWVoclFJeQWcx5dHatWsJDw/nhx9+ICwsjOHDhzs990QppQoS7SBuU2pqKp06deKpp54iPT2dlStXEh0drcVBKeV1tIO4DbGxsXTo0IGUlBS6devGwIEDufvuu62OpZTyEnkd9/3TTz/xl7/8hYyMDK5evcpbb71Fhw7OvxdIC0QunDx5km7dujFt2jRsNhsbN2502UfZlVKeI/HcORLOniW0RAlC7r3Xbc9zu+O+//SnP7Fx40buvPNOLly4QLVq1XjxxRcpU6aMUzn0ENNNGGOYNWsWNpuNmTNn0q9fP7Zv367FQSkflHjuHE127aLvjz/SZNcuEs+dc8njumLct7+//+8zmn777TcyM10zoUg7iBs4duwYHTt25NtvvyU4OJiVK1dSo0YNq2MppSyScPYsaZmZZABpmZkknD3rdBeRfdx3eno6tWvXpk6dOtetd23cd7du3Wjbti0bNmzgypUrBAYG/n4o6ciRI7Ro0YKDBw/yySefON09gHYQ1zHGEB0djc1mY/ny5XzyySckJiZqcVDKx4WWKIF/oUL4Af6FChHqgjemZB/3Xbx48VyN+65fvz7FihXj/vvv/33cN9iLyO7duzl48CBTpkzhl19+cTqfFohsDh06RNOmTWnfvj01a9Zkz5499OzZkzvu0EZLKV8Xcu+9xAcFMbBSJeKDglx2DsJV476vKVOmDIGBgaxbt87pbFogsA/XGzlyJNWqVWPr1q1MmDCBVatW8eijj1odTSlVgITcey99KlZ0WXFw1bjvlJQULl++DMCZM2fYsGEDVapUcTqfz/9pfObMGZ5//nk2b95MixYtmDBhAuXKlbM6llLKB7hq3Pf+/fvp0aMHIoIxhp49e1K9enWn87lt3LcV8jLu2xjD66+/TosWLfjb3/6Wq3ZPKeX5dNy3nVXjvj2CiDB9+nSrYyilVIGj5yCUUko5pAVCKeWzvOkQ+63k5bVqgVBK+aQiRYpw6tQpnygSxhhOnTpFkSJFbms7nz8HoZTyTeXKlSMlJYUTJ05YHSVfFClS5LbfoakFQinlkwoXLkylSpWsjlGg6SEmpZRSDmmBUEop5ZAWCKWUUg551SepReQE8JPVOW6iNHDS6hAFjO6T6+k+uZ7ukz9y5f6oaIy539EdXlUgCjoRSbrRR9p9le6T6+k+uZ7ukz/Kr/2hh5iUUko5pAVCKaWUQ1og8leU1QEKIN0n19N9cj3dJ3+UL/tDz0EopZRySDsIpZRSDmmBUEop5ZAWCBcQkeYickBEDopI75usV1dEMkSkVY7b/URkh4gsdn/a/OHMPhGREiIyR0S+E5H9IhKSP6ndy8l90k1E9onIXhGZISK3N5azgLrVPhGRUBE5JyI7sy79crutp8rrPhGR8iKyOutnZp+IvON0GGOMXpy4AH7Af4CHAX9gF2C7wXqrgFigVY77ugNfA4utfj0FYZ8AU4DwrOv+QAmrX5OV+wQoC/wI3JW1PBtoa/Vryo99AoQ6+rnI7f70tIuT++RPQO2s68WA753dJ9pBOK8ecNAYc8gYkwbMBF5ysN5bwFzg1+w3ikg5oAUQ7e6g+SjP+0REigONgRgAY0yaMeas+yO7nVP/n2CfvHyXiNwBFAWOuTNsPsntPnH1tgVZnl+XMea4MWZ71vXzwH7sf1zkmRYI55UFjmRbTiHHfxQRKQu8AkxwsP0o4F9AprsCWsCZffIwcAKYlHXYLVpE7nZn2HyS531ijDkKDAf+CxwHzhljlrs1bf645T7JEiIiu0RkqYgE3ua2nsaZffI7EQkAagGbnQmjBcJ54uC2nO8dHgX0MsZk/GFDkReAX40x29wVziJ53ifY/1KuDYw3xtQCLgLecHzZmf9P7sP+V2QloAxwt4i87paU+Ss3+2Q79llBQcBoYMFtbOuJnNkn9gcQuQd7F9rVGJPqTBj9wiDnpQDlsy2X4/r2PxiYKSJgH7L1ZxFJB+oDL4rIn4EiQHERmWaM8fQffmf2ySYgxRhz7S+fOXhHgXBmnxQGfjTGnAAQkXnAE8A0d4d2s1vuk+y/4IwxsSIyTkRK52ZbD5XnfWKMOSkihbEXh+nGmHlOp7H6pIynX7AX2UPY/7q7dlIp8CbrTybHSWpzkxNPnnhxdp8A64AqWdf7A59Y/Zqs3CfY/5DYh/3cg2A/if+W1a8pP/YJ8BD/+0BvPeyH2eR296enXJzcJwJMBUa5Ko92EE4yxqSLSBcgDvs7EL40xuwTkQ5Z9zs67+DVXLBP3gKmi4g/9h+Wf7o1cD5wZp8YYzaLyBzshxbSgR14weiJXO6TVkDHrE7qMtDa2H8zOtzWkhfiQs7sExF5EngD2CMiO7Me8j1jTGxe8+ioDaWUUg7pSWqllFIOaYFQSinlkBYIpZRSDmmBUEop5ZAWCKWUUg5pgVAqj0QkQET2Org9WkRsN9imq4gUzbZ8wZ0ZlXKGFgilXMwYE26MSc55u4j4AV2xf+BNqQJPC4RSzrlDRKaIyO6s77AoKiIJIhIM9g5BRAaIyGbgfeyzlFaLyOprDyAig7MGr20SkQcteh1KXUcLhFLOqQJEGWNqAKlApxz33w3sNcbUN8YMwD5X52ljzNPZ7t9k7IPX1gLt8ym3UrekBUIp5xwxxmzIuj4NeDLH/RnYh6fdSBpw7ZsEtwEBLk2nlBO0QCjlnJyzanIuXzHXjzTP7qr537ybDHTCsipAtEAo5ZwK2b4z+2/A+lusfx7710EqVeBpgVDKOfuBN0VkN1ASGH+L9aOApdlPUitVUOk0V6WUUg5pB6GUUsohLRBKKaUc0gKhlFLKIS0QSimlHNICoZRSyiEtEEoppRzSAqGUUsqh/we9DIjpDQWcpAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_diagram(diagrams_hkh_1[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The picture above represents an exemple of diagram used for feeding ATOL. The HKS algorithm gives as output 4 subdiagrams. Note that in the plot we used the terminology dim0, dim1, dim2, dim3 only in order to distinguish the four subplots. They are not to be considered as Homology dimensions in the real TDA meaning!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_chemical_info(smiles):\n",
    "    stereo = np.zeros(3)\n",
    "    for char in smiles:\n",
    "        if char in '/':\n",
    "            stereo[0] +=1\n",
    "        if char in '@':\n",
    "            stereo[2] += 1\n",
    "        if char in '\\\\':\n",
    "            stereo[1] += 1\n",
    "    chemical_el = 'B C N O P S F Cl Br I * b c n o s p'.split()\n",
    "    chemicals = np.zeros(len(chemical_el))\n",
    "    smiles = '|'+smiles+'|'\n",
    "    for i, el in enumerate(chemical_el):\n",
    "        if el != 'C':\n",
    "            chemicals[i] += len(smiles.split(el)) - 1\n",
    "        else:\n",
    "            chemicals[i] += len(''.join(smiles.split('Cl')).split(el)) - 1\n",
    "    return np.concatenate([chemicals, stereo])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['chemistry'] = df.smiles.apply(add_chemical_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>smiles</th>\n",
       "      <th>FDA_APPROVED</th>\n",
       "      <th>CT_TOX</th>\n",
       "      <th>smiles_graph</th>\n",
       "      <th>chemistry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>*C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....</td>\n",
       "      <td>[0.0, 14.0, 2.0, 7.0, 0.0, 0.0, 0.0, 0.0, 0.0,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[C@@H]1([C@@H]([C@@H]([C@H]([C@@H]([C@@H]1Cl)C...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1.]]]], ...</td>\n",
       "      <td>[0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[C@H]([C@@H]([C@@H](C(=O)[O-])O)O)([C@H](C(=O)...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1....</td>\n",
       "      <td>[0.0, 6.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[H]/[NH+]=C(/C1=CC(=O)/C(=C\\C=c2ccc(=C([NH3+])...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....</td>\n",
       "      <td>[0.0, 10.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[H]/[NH+]=C(\\N)/c1ccc(cc1)OCCCCCOc2ccc(cc2)/C(...</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>[[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....</td>\n",
       "      <td>[0.0, 7.0, 4.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              smiles  FDA_APPROVED  CT_TOX  \\\n",
       "0            *C(=O)[C@H](CCCCNC(=O)OCCOC)NC(=O)OCCOC             1       0   \n",
       "1  [C@@H]1([C@@H]([C@@H]([C@H]([C@@H]([C@@H]1Cl)C...             1       0   \n",
       "2  [C@H]([C@@H]([C@@H](C(=O)[O-])O)O)([C@H](C(=O)...             1       0   \n",
       "3  [H]/[NH+]=C(/C1=CC(=O)/C(=C\\C=c2ccc(=C([NH3+])...             1       0   \n",
       "4  [H]/[NH+]=C(\\N)/c1ccc(cc1)OCCCCCOc2ccc(cc2)/C(...             1       0   \n",
       "\n",
       "                                        smiles_graph  \\\n",
       "0  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....   \n",
       "1  [[[[[0. 1. 0. 0. 0. 1. 0. 0. 0. 0. 0. 1.]]]], ...   \n",
       "2  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1....   \n",
       "3  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....   \n",
       "4  [[[[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0....   \n",
       "\n",
       "                                           chemistry  \n",
       "0  [0.0, 14.0, 2.0, 7.0, 0.0, 0.0, 0.0, 0.0, 0.0,...  \n",
       "1  [0.0, 6.0, 0.0, 0.0, 0.0, 0.0, 0.0, 6.0, 0.0, ...  \n",
       "2  [0.0, 6.0, 0.0, 8.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...  \n",
       "3  [0.0, 10.0, 4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0,...  \n",
       "4  [0.0, 7.0, 4.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, ...  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## We split the dataset in train and test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "diag_train, diag_test, ext_feat_train, ext_feat_test, y_train, y_test = train_test_split(diagrams_hkh, df.chemistry, \n",
    "                                                                                         df.CT_TOX, test_size=0.2,\n",
    "                                                                                         random_state=52, \n",
    "                                                                                        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Validation for Simple Atol usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_validate(parameters, X, y, k_fold=5, **atol_params):\n",
    "    set_up = []\n",
    "    for parameter in parameters:\n",
    "        pipeline_model = Pipeline([('atol_layer', Atol(parameter, **atol_params)), \n",
    "                                   ('model', RandomForestClassifier(n_estimators=100))])\n",
    "        set_up.append((parameter, cross_val_score(pipeline_model, X, y, scoring='roc_auc', cv=k_fold)))\n",
    "    return set_up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We get a 5-fold cross validation accuracy estimation for the simpler model using ATOL! Note in this exemple only notions about the structure of the molecules are used! **No chemical info**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_res_hkh = cross_validate(np.arange(1, 25).astype(int), diag_train, y_train, \n",
    "                            5, padding=0, ignore_hom_dim=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEHCAYAAAC9TnFRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXjb53Xg++8BQIALuIkgJZLaN8uSvMtbYstbEitL4zRLEzuJk3Q6jjvxTGfa27npbZNpp5NpWrdzM504dR3HbTbXk9s4sZt6bCfeHW+SbFmWrJVauInivu/AuX8APwqiQBA7IPJ8nkePSfAH4BVM4eB933POK6qKMcYYMx9XvgdgjDHm/GABwxhjTEIsYBhjjEmIBQxjjDEJsYBhjDEmIZ58DyCTAoGArl69Ot/DMMaY88bu3bu7VbU2kWsXVMBYvXo1u3btyvcwjDHmvCEiJxO91pakjDHGJMQChjHGmIRYwDDGGJMQCxjGGGMSYgHDGGNMQixgGGOMSYgFDGOMMQmxgGEWtP3tA7T1j+V7GMYsCBYwzIL25R/u5t4nD+Z7GMYsCBYwzII1PhWkrX+MUwPj+R6KMQuCBQyzYLX1j6EKXcMT+R6KMQuCBQyzYDX3jgLQPWQBw5hMsIBhFqyWSMAYHJ9mYjqY59EYc/6zgGEWrOae0Zmve4Yn8zgSYxYGCxhmwXKWpAC6bR/DmLRZwDALVnPvKDVlXsAChjGZYAHDLEiqSkvvKJetrAage8iWpIxJlwUMsyD1jkwyMhnk8lVVgKXWGpMJFjDMguTsX2ysK6fM67YlKWMywAKGWZCcgLGyppTach9dVothTNosYJgFyanBWFFdSsDvsxmGMRmQ1YAhIjtE5JCIHBWRr85xzY0iskdE9ovIC1G3nxCRdyI/25XNcZqFp6V3jNpyHyVedyRg2Ka3MenyZOuBRcQN3Ae8H2gFdorI46r6btQ1VcB3gB2q2iwidbMe5iZV7c7WGM3C1dw7ysolpQAEyr28dtxmGMakK5szjKuAo6p6TFUngUeA22ZdcwfwqKo2A6hqZxbHYxaRswKG30f/6BRTwVCeR2XM+S2bAaMRaIn6vjVyW7SNQLWIPC8iu0XkzqifKfB05Pa75noSEblLRHaJyK6urq6MDd6cvyanQ5waGGNFVMAAaw9iTLqytiQFSIzbNMbzXwHcApQAr4rIa6p6GHivqrZHlql+KSIHVfXFcx5Q9QHgAYBt27bNfnyzCLX3jxFSzpphQLjae1llcT6HZsx5LZszjFZgRdT3y4H2GNc8qaojkb2KF4FLAFS1PfLfTuBnhJe4jJnXTEptJGDUlofbg1jxnjHpyWbA2AlsEJE1IuIFPgM8Puuax4DrRcQjIqXA1cABESkTkXIAESkDPgDsy+JYzQJyTsDwh2cVVothTHqytiSlqtMicg/wFOAGHlLV/SJyd+Tn96vqARF5EtgLhIAHVXWfiKwFfiYizhgfVtUnszVWs7C09I7i9bioKw8vRQXKrQGhMZmQzT0MVPUJ4IlZt90/6/t7gXtn3XaMyNLUYtcxMM49D7/Jdz53OXXltv6eiObeUVZUl+ByhbfRSr0eSr1ua0BoTJqs0rvAvXGil10n+9jbMpDvoZw3olNqHVbtbUz6LGAUuLa+McCWUxKlqjT3xAoYXnsNjUmTBYwC19Yf3sDtGbHllEQMjE0xNDE9U4PhsBmGMemzgFHgWiMzDMvwSYyTIXVOwCi3flLGpMsCRoGzJankzE6pddT6ffSNTlp7EGPSYAGjgKnqzAzDAkZiWnrDr1esGYZq+CQ+Y0xqLGAUsL7RKcamgoD1QUpUc+8oNWVe/L6zM8Zr/ZFqb1vaMyZlFjAKWGtfeHllWUWxzTAS1NI7es7sAs7uJ2WMSY0FjALm7F9csqKSPmvPnZBYNRgQHTBspmZMqixgFLDWmYBRBdj6+3ymgyHa+sdiB4xym2EYky4LGAWsrX+M8mIPawNlgK2/z+fUwDjBkMYMGGVeN8VFLrrtNTQmZRYwClhr3yiNVSVnDgCyGUZcc9VgAIiIFe8ZkyYLGAWstW+M5dWlZ9bf7dNxXDM1GDXnBgyAWiveMyYtFjAKWFvfGMurS2z9PUHNvaMUuYVlFbG7+gb8PlvWMyYNFjAKlNMTaXl1CWVeNz6PywLGPJp7R1leXYrbFet0YOsnZUy6LGAUKKcGo7GqZGb93Yr34purBsNR6/fSOzrJtKUnG5MSCxgFyqnBWF4dfgMMlPvsTOp5OAcnzWWmPcioBV5jUmEBo0A5NRiNkTfAWr/XNmzjGBibon90KmZKreNM8oC9jsakwgJGgWrrH6PU66a6tAiAmjJbf4+nZY4utdGsPYgx6bGAUaCcGgyR8AZuoNxL78gkoZDmeWSFydnzibeHEYg0ILSAYUxqLGAUqLb+cEqtI+D3EQwp/WNTeRxV4ZqvBgPCdRhgAcOYVFnAKFCtfWMz+xdgyynzae4dpaq0iIriojmv8fs8+Dwuq8UwJkUWMArQ8MQ0/aNTMxlSgFV7z6O5N3bTwWhn2oPYprcxqbCAUYCclNrGqugZRuQAIJthxDRfDYYjfLa3vYbGpCKrAUNEdojIIRE5KiJfneOaG0Vkj4jsF5EXZv3MLSJvicgvsjnOQtPWH16Pn72HAXbyXizBkNLaF/scjNlq/V5bkjImRVkLGCLiBu4DPghsBm4Xkc2zrqkCvgN8VFW3AJ+a9TC/BxzI1hgL1ewaDIDKkiI8LrFPxzF0DI4zFYzd1nw2W5IyJnXZnGFcBRxV1WOqOgk8Atw265o7gEdVtRlAVTudH4jIcuDDwINZHGNBausbw+dxURuZVQC4XEKN32sBI4bmnvlrMBwBv4/ekQmClp5sTNKyGTAagZao71sjt0XbCFSLyPMisltE7oz62beA/wwsusY/rX1jZ9VgOMLFe/bpeLZEivYcAb+XkEKftQcxJmnZDBixWobO/ljnAa4gPJO4FfiaiGwUkY8Anaq6e94nEblLRHaJyK6urq60B10IWvvPTql12IZtbM29o7hdQn1l7Lbm0WrLw9fk+3VUVV4+0s0f/+wdOgbG8zoWYxLlyeJjtwIror5fDrTHuKZbVUeAERF5EbgEuBz4qIh8CCgGKkTkR6r6udlPoqoPAA8AbNu2bUGsM7T1jbJ589Jzbg/4vTR1DudhRIWtuXeUhqpiPO75P//MZJsNTbBpWbZHdq6J6SCP72nney8f52DHEABLK4r5D7dsyP1gjElSNmcYO4ENIrJGRLzAZ4DHZ13zGHC9iHhEpBS4Gjigqn+kqstVdXXkfs/GChYL0fhUkO7hybNqMBy1/nDHWtUFERczprk3sQwpIG+HUfWNTPLtZ49w3V8+xx/+815U4a8+eTGb6yt48fDCmBmbhS9rMwxVnRaRe4CnADfwkKruF5G7Iz+/X1UPiMiTwF7CexUPquq+bI3pfNAaowbDEfD7mJwOMTQxHbeiebFp6R3lA1vOnZHFkuuOtU1dwzz08nF++mYr41MhbthYy+/81hquWx9ARDjZM8L9LxxjcHzK/p+agpfNJSlU9QngiVm33T/r+3uBe+M8xvPA81kYXkFq63fOwTg3YNQ4zfOGJuzNJWJkYpqekcmEivYAKoo9eN3ZPb1QVXntWC8PvnSMZw524vW4+Phljfz2dWvYuLT8rGu3b6jlvueaeOVoDzu25mGNzJgkZDVgmOTNnLQXa9N7pp/UJGtrczqsgtXSl3iGFDjtQbxZqZifnA7xr++08+BLx9nfPkhNmZffu2UDn7921cz/u9kuX1VNmdfNi0e6LGCYgmcBo8C09Y1R5Bbqys/N+DlT7W2ZUo5kajAc4WyzzC5JNXUN89nvvk7H4Djr6/x88+MX8bHLGikucse9X5HbxbXrArx4uAtVPSeV2phCYr2kCkxr3xj1lSW4Xee+cQTK7TyH2ZqTqMFwBPy+jDdx/D/vnKJjcJx/+OKV/PI/beczV62cN1g4btgYoLVvjBOR4FcoXjrSxRceesPOQDczLGAUmNnnYERbUupFBLqseG9GS+8o5cUeKksS39Op9We+nuV49yjLKoq5aVNd0rOE7RvD64uFli31o9dO8sLhLk72FlYgM/ljAaPAOCftxeJxu6gutfYg0ZyU2mTepAPlXnoyfHrh8e5h1gTKUrrvqpoyVi4pLaiAMTEd5KUj3QActdofE2EBo4BMTAfpHJqIWYPhCPi9tocRJZkaDIdzemEm24Mc7x5hTW1qAQNg+8YArx7rYXK6MJZ/3jjey+hkELCAYc6wgFFATvWPoxo7Q8ph3VbPCIWUlr75D06aLTrbLBP6RyfpG51ibYozDAin145OBtl9si8jY0rXMwc68Xlc1l3AnMUCRgGJV4PhCGRh/f181Tk0weR0KOEaDEemj7s93j0CwOqa1APGtetq8LiEF4/kf1lKVXn2YCfvWVfDpmUVHO2ygGHCLGAUkJkajDn2MCBcvGfHtIalkiEFUJvhbDMnYKSzJFVeXMTlK6t5qQACRlPXCM29o9x84VLW1/k52jls7WgMYAGjoLT1jc3bdTXg9zEyGWQssr68mKUaMJwZRqZO3jvePYLbJayIs/eUiO0bA+xrG8z7DPLZg6cBuHlTHevr/IxOBmm3jroGCxgFpbVvjGUV8buu1mZ4OeV81tw7igg0xJmRxVJZUkSRWzK2h3G8e4Tl1SV4Pen9c7p+Qzi99uVIdlK+PHuwk03LymmsKmF9nR+wjW8TZgGjgMx1DkY0K947o6V3lIbK5N+ow+1BMrcXdLx7JOWU2mhbGyupLi3K6z7GwNgUO0/0cfOmOgALGOYsFjAKSFvfGMvn+bSc6Qyf81lz7ygrliQ3u3AE/L6MLEmpasYChtslXLehlpeOdOdtz+ClI10EQzoTMGrKvFSVFlnAMIAFjIIxHQzRMTgeN0MKoMaWpGa0pFCD4Qhk6Hz0rqEJRieDGQkYANs3BOgamuDAqaGMPF6ynj3QSVVpEZetrAbCs7H1tX5LrTWABYyCcWpgnGBI512SqikLL0kt9uK9sclwkWPqASMzS1LHnAypDAUMZx8jH9lSwZDy3KFObrqg7qxeZuvr/JZaawALGAXDqcForIr/Blhc5Ka82LPol6ScFORkazAcgXIfPcPptwc5nuGAsayymAuWludlH2NPSz99o1PcFFmOcqyv89M7MrnoP6QYCxgFwzlpb74lKThzVOtilmpKrSPg9zEdUgbGptIax4nuEbweFw2Vqe2lxLJ9Y4Cdx/sYnZzO2GMm4tmDp3G7hBs2nH3Yim18G4cFjALRFgkY9VVz12A4rHgvEwEjM9lmx7pHWF1TiitGO/pUXb+hlslgiNeP9WbsMRPx7MEurlhVTWXp2Z1/ZwKGLUstehYwCkRr3yhLK3z4PPOfoWDtQcIBo8zrZklkTydZTj1LujO1TGVIRbtqzRJ8HldOl6Xa+8c4cGqQW2YtRwE0VJZQUuS2GYaxgFEo2vrH4rYEiRbw++gZWdx7GC29o6xIsq15tNry9NOTgyGluWeUNQF/yo8RS3GRm6vX1uS03fmzBzsBZtJpo7lcwrq6MgsYxgJGoWjtG4vb1jxawO+jf3SKqUV8Eloqbc2jzdSzpLG0194/xmQwxJpAei1BYtm+IUBT18hMMkS2PXewkxVLzlR2z2aptQYsYBSEYEg5NTB/lbfDqfbuWaSZUqqadsCoLCnC45K0lqTOpNRmdoYBZ07heykHs4zxqSC/burmlk1L55yxbVhaTvvAOMMT2dmIt+aG5wcLGAWgc2icqaAmlCEFUFO2uIv3uoYnGJ8KsbIm9YDhcknayQPHI5vAmd7DANhQ52dZRXFO9jFebephfCp0TjpttHW14aCYjVnGdDDE9nuf4y/+z4GMP7bJLAsYBcDJkEp0DyPT7bnPNy296dVgONJNHjjRM4rf55nJuMokEeH6DQFePtLNdJaXHp85eJpSr5ur1yyZ85psptYe6RympXeMv3/hGI+80Zzxx492qGOIx/a0ZfU5FrI5A4aIfE5EPh/j9n8rIndkd1iLy5kajMT3MGDx9pNyUmrTbSee7umFxyIZUqluvM9n+8ZaBsen2ds2kJXHh/BS0HMHu3jv+gDFRXNn6K2qKcXjkqyk1r7d0g/A1sYK/uTn+3jtWE/GnwPgaOcQn37gVX7vkT2cHrR27amIN8P4A+DnMW5/JPKzeYnIDhE5JCJHReSrc1xzo4jsEZH9IvJC5LZiEXlDRN6O3P5niTzf+epMlXfiWVKweGcYzT2JFznGk+4M43j3cFaWoxzXrQ8gQlazpQ6dHqKtfyxmOm20IreL1YHsZEq93dpPRbGHH//ONayqKeV3f7Sb5p7RjD5HW/8Yn//eGwQjlf3PRbLCTHLiBQy3qp7TAS1yW1GM688iIm7gPuCDwGbgdhHZPOuaKuA7wEdVdQvwqciPJoCbVfUS4FJgh4hck8DfJ2mqSufQeF7bHrT2jRLweynxzl+DAVDqdVNc5Fq0xXvNvaMsqyiO+4k4EYFyLz3DkyltuE5MB2nrG2N1FgNGdZmXi5dXZTVgPHMg/MYZb//Cka1Mqbea+7lkRRWVJUU8+IUrCSn8m+/vZGg8vSp8R+/IJHd+73WGx6d55K5raKwq4RkLGCmJFzCKROScfw0iUg4ksmh7FXBUVY+p6iThmclts665A3hUVZsBVLUz8l9VVec3syjyJytpFKrw3m8+y3dfOp6Nh09Ia1/iNRiQ+fMczjctfellSDlq/T4mgyEGx5LP/GnpHSWksDaLAQPC6bV7WvrTbmEyl+cOdrK1sYKlFfN3GNiw1M/J3lEmpzO3pzI6Oc3h00NctqIKCCcQ/N1nL+d49wj/4Z/empkRpGpkYpov/eNOWvvGePAL29jSUMktF9bx8pFuxqfs1MpkxQsY3wP+WURWOzdEvn4k8rP5NAItUd+3Rm6LthGoFpHnRWS3iNwZ9VxuEdkDdAK/VNXXYz2JiNwlIrtEZFdXV/KfxFwuYVllMe05ynePpS2JGgzHYi7ec4r20uUU76WSWnusK7NNB+eyfWMtIYVXjmb+FL6+kUnebO7j5k1LE7p+fZ2fYEg50TOSsTHsaxskpHBJJGAAvGd9gD/96BaeO9TFN9PInJqcDnH3j3azr22Ab99xOVevrQHCxYljU8Gs7ZUsZHMGDFX9a+Ax4AUR6RGRHuAF4Beqem8Cjx1rJ3D2xwUPcAXwYeBW4GsisjHy/EFVvRRYDlwlIlvnGOcDqrpNVbfV1tbGumRe9ZUlnBrIT8BQ1XCVd5Lr8Zk6AOh8Mz4VpGNwPCMzjHTO9nbeNLO5JAVw6Yoqyn2erKTXvnC4i5DGru6OxUmtPXI6c8tSzob3xcurzrr9c9es4s5rV/Hdl47zk10tse4aVzCk/P5P9vDSkW7+4uMX8f7NZ4LiNWtrKClyz1S3m8TFTatV1ftVdRWwClitqqtU9e8SfOxWYEXU98uB9hjXPKmqI6raDbwIXDJrDP3A88COBJ83aY1VJbT35ydromt4gonpUNIbuLXl3kWZJdXWP4YqrKxJvztsOskDx7tHqCnzUlky73ZeWorcLq5dV8OLhzN/Ct8zBzsJ+L1c3FiZ0PXrav2IZDa1dk9rP41VJTOzvWhf/8hmrlsf4I9/9g47TyTeiFFV+bN/2c8v9p7ijz64id/atuKsnxcXubluQ4BnDnRawWCS4gYMEdkqIt8n/Ib9rIh8X0QuSvCxdwIbRGSNiHiBzwCPz7rmMeB6EfGISClwNXBARGojG+KISAnwPuBgwn+rJNVXFtMxOJ72emkqkq3BcNSU+egdmcjLmPMp3S610dLpWHusK/NNB+eyfWMtbf1jM5XlmTAdDPHCoU5uvKAu4U67JV43jVUlGU2t3dPcz6Urq2L+zON2cd8dl7OiupQv/3D3TP3NfP7nM0f4wasnuWv7Wr58w7qY19yyqY62/jEOZ3C2tBjEq8O4DfgZ4WWo3wZ+J/L1o5GfxaWq08A9wFPAAeAnqrpfRO4Wkbsj1xwAngT2Am8AD6rqPqAeeE5E9hIOPL9U1V+k/teMr76qhGBI87LE46TUJr+H4SWk0D+6uGYZmSraA6gu9eJ2ScozjFwFjBsibUIymS21+2Qfg+PT86bTzra+zp+xGUbX0ARt/WNcujx2wACoLC3iwS9sYzoY4ne+v2ve1iQ/fPUE3/rVET55xXL+6IOb5rzOyQp75uDplMa+WMWbYfxX4P2q+pCq7lXVt1X1IeD9kZ/NS1WfUNWNqrpOVb8Rue1+Vb0/6pp7VXWzqm5V1W9Fbturqpep6sWR2xN6vlQ1Rs6gyFWjt2hO0V7SexgZ6LZ6PmruGaW4yDXTnjwdLpewpMxL91Byr+HIxDSdQxNZ379wrFhSyuqaUl46krmN72cPdlLkFq7bEEjqfhvq/BzrGs7IzHZva3j/InrDO5a1tX7u++zlHO0a5j8+Mnfm1C/2tvP1x/fzvgvr+ObHL4pbULm0opitjRVWj5GkuGm1qnpi9o2R27K7cJtj9ZHT0vKx8d3WN0ZVaRF+nyep+y3W4j2n6WCmqqtrU0hPdo5lzXZKbbTtG2t5tamHienMpII+e7CTq9Ysobw4uX/K6+v8TEyHZpZS0/F2Sz8uCVd4z+f6DbV8/SOb+dWBTu596tA5P3/pSBf/6X/vYduqar59x+V43PN3Pbp501J2n+yjb5FmG6Yi3qs6JSIrZ98oIquA3J4dmWUNkf2DU3nY+G7tG016/wIsYGRKoDz1gLGmNocBY0MtY1NBdp/oS/uxWnpHOdI5nHA6bbQzp++dU9ObtD2tA2xcWk6pN7EPS3deu4rPXr2S+19o4qe7W888Tks/X/7hbtbV+nnwC1cmXNB5y6Y6QhrOFjOJiRcw/gvwKxH5oohcFNkA/xLwNPD13AwvNyqKPZR53XlZkmrrH0upxYWzYbuYUmtVNWM1GI6AP/lssxORgLFqSe4CxjXravC4hBczsCwV77Ck+ayvLQfST61VVd5u6eeyOTa8YxER/vSjW7h2bQ1/9Og77D7Zy9HOYb70D29Q4/fyg9++KqmstYsaKwn4fVb1nYR4dRg/J9yq42bgH4EfADcBvxX52YIhItRX5b4WQ1UjVd7JvwFWlhRR5JZFVbzXOzLJyGQw7aaD0Woj9SzJpFce7x6hobI44VYumeD3ebhiVXVGNr6fOdjJ2kBZSpv2laVFBPy+tDe+T/SMMjA2xSVxNrxjKXK7+M5nL6e+qpgv/3A3d37vddwuFz/87aupS6BaPZrLJdy8qZYXDnUu6sPIkjFfHcbbqnqnql6hqpdHvn47siy1oDRUlXBqILdLUv2jU4xOBlOaYYgINWW+RdVPKpMptY6A0x5kPPFV1mPdIzldjnJs31jLu6cG05pVjkxM81pTT0K9o+ayvq4s7dRap2Bvvg3vWKrLvHzvC9uYmAoxND7N93/7ypQTEG7etJTB8Wl2n0x/qW8xmK8O41oR+aSI1EW+v1hEHgZezsnocqghD+1BUs2QcgTKvRndw3jgxSYeeLEpY4+XaTMBI42Dk2YLpHC2yImeEVbX5D5gOOm1Lx9NfZbx66PdTAZDSafTRnNSa9MpetvT0k9JkZsNcxwJO/8YyvnZV97LY/e8ly0NiRUexnLdhgBFbrGq7wTFq8O4F3gI+ATwryLyX4BfAq8DG3IzvNypryyhe3gyY1koiWjrD78Bptqmu6YsvfMcZnvwpeP89ycO8qPXTmbsMTOpJUPnYERL9mzvvpFJ+kenclaDEW1zfQU1ZV5ePJz6PsazBzsp93nYtnruw5Lms6GunKHx6bRmOnta+rmosTKhbKa5rK/zs7Y2veNx/T4P16ytyUrAOHBqMKfvJ7kQ7//Wh4HLVPV24APAV4HrVPV/quqCO32kIVKL0ZHDZamZg5NS2MOASAPCDM0w+kYm6RyaoNzn4euP7ePZAixoaukdo7bcl9G9g2QPo3KqrdfmYUnK5QrXTbx0pItQCnUQqsqzBzu5fmMArye9N2pIvUXI5HSId9sH56zwzrWbN9VxtHOYkxlsqnjg1CAf+tuX+F/PHM3YYxaCeL81Y05gUNU+4JCqHsnNsHLPSa3NZaZUa98Y5T4PFSXJ1WA4ApF+Upnoh3OwI5wm+VefvJjNDRXc8/BbvNOavZPeUpHplFo407E20SWpmZTaQHqfbFO1fUMt3cOTHOgYTPq++9sH6RyaSCmdNpoTMI6kGDAOdgwyGQwlveGdLU62WCZnGfc9dxRVeGRny4LaUI8XMNaJyOPOH2D1rO8XlHzUYrT2hbvUplqEVpvChu1cDp8OB4zLV1Xz0BeupLrUy29/fyetfZk9+Swd2QgY1aVeXJJ4wDjRPYLbJWmf9peq6yOV2X/55CGe3HeKwSQOGXrmQCcicOMFqXV1dtSV+yj3eVKeYZzZ8E597yGTVtWUsb7On7GA0dQ1zL++c4pLVlTRPTzBL98tvNl6quJ9tJ3dL+pvsjmQfKuvDC9J5TK1NtUaDEd08V66XVMPdgxRWVJEXbkPEeEfvnQln/i7V/jSP+zkn3/3PVnvyjqfyekQpwbGMlqDAeB2CUvKEi/eO949wsolpRSlsfaejrqKYr58w1p+/FozLx7uwu0SLl9ZxfYNtdxwQS1bGyrnbCb47KFOLlleNfN7kyoRYV0aPaX2tAwQ8HtTKljNlls21fHQr48zPDGddNeF2f7u+SZ8Hhff/fwV/OZ3XuHh15v50EX1GRppfsWrw3hh9h9gKOrrBaW4yM2SMi9tOZ1hpFbl7ahxuq1mILX2UMcgFywrn5ntbFxazt9/7gpO9Ixw9w93Z/SUtdTGN0RIYV0W9g4Cfm/CG7jHcth0cC5/9MELeevr7+cnX76Wu29Yy/hUiL/55WE++u1fs+0bv+L3HnmLR99sPevv1DU0wdst/SkV68Wyvs6fcmrtnpY+Ll1RlbH2Lplw86Y6poLKy2meO9LSO8rP3mrj9qtWUldRzKevXMHLR7tnij3Pd8l+THowK6MoEA1VxTmbYQyMTTE0Pp10l9pozifFdIv3VJXDp4fZtKz8rNvfsz7ANz9+Ma8e6+GrP92b17MDXmkKZwZdGzk1LXX3czgAACAASURBVJNqy310JbDpraqc6M5PSu1sRW4XV61Zwh/euol/+ffXsetP3se3Pn0pN2ys5eUj3fz+T97mym/8ig//7Uv85ZMHefDlY0Bq1d2xbKjz0zU0kfTRsYPjUzR1jRTM/oXjilXVVBR7Zs44T9Xfv9iES+Cu7WsB+PSVK3C7hH/a2ZyJYeZdsnOvwvlIkAX1lSUZzZSIpy3NGgzIXD+ptv4xhiemuWBWwAD4xBXLae0b4//91WGWLynl99+/Ma3nStUrTT2sr/MnXc2biIDfN3PkajynBycYmwrmpWhvPgG/j49d1sjHLmskFFLePTXIC4e7eOFwF9998RjTIWVphY8tDfM3+ktEdKbUFauqE76fk0iRSsFeNnncLm64oI7nDnUSCmnCZ4REOz04zk92tvLJK1bMNDRdWlHMLZvq+OddrfzB+y9IKzutECQbMP4sK6MoEA2VxbzWlJtzfs+cg5F6wFhS5kUk/SWpQ5EMqQuWnhswAP7DLetp6Rvlb585wvLqknNOMMu2qWCInSd6+eQVy7Py+OF+UuH2IPGWSY51h5dgctmlNhUul7C1sZKtjZV85ab1DI1P8UpTD8sqijO2DOQEjKYkA8aemSNZC2PDO9otm+r4l7fbeadtIKWA9sCLxwiq8ruzDm367DWrePrd0zy1v4PfuKQhU8PNi3iFe7eKyCejb1PVn4vIZ0Xk/dkfWu41VJUwNDGdVOZJqpzso3T2MNwuYUmpN6HllHiclNqNMWYYEN7k/IuPX8R16wP8P4++w8sZPJchEXtb+xmdDGZlOQrCn84npkPzHs5zojv8/yxX52BkSnlxEbduWZbRT/XLq0vxelwc6Uyua+3bLf2sCZRRVerN2Fgy5YaNtbiElJoR9gxP8PDrzdx2ScM5nQiuXx9gxZISfvx6YRbEJiPe/OjPCJ+wN9szJHiA0vmmPoeptW19Y5RENtrTkYnivUMdQzRWlVAR52yEIreL73zuctbX+fndH+3mYAp1AKl65Wh41ndNlgJGbYKHUR3vHsbncVGfhWWx843bJawNlCWVKaWq7Gnp59ICW45yVJd5uWJVdUpFqw/9+jjj00H+3U3nHgnrcgmfuXIlrx3rpSmDx9vmQ7yAUaqq56QMqGoHcH59xEpQQyS1tj0HG9/p1mA4MtFP6lDHUMz9i9kqiot46ItXUupz86V/2JmzqvhXmnrYXF9BdZrBdS6J7gU5x7Kmsr69ECWbKdUxOE7n0ASXFOBylOPmTUvZ1zaY1O/2wNgUP3jlJB/aWs/6utj/jj61bTkel/BPr2d+83twfIqhHKyKQPyAUSwi5+xxiEgRUDgJ1BmUy+K9dGswHAF/ev2kJqdDNHUNJxQwIPwaPfTFKxkcm+JL/7hz3mWcdI1PBdnd3Md71mVndgGJ95M6ViAZUoViQ105rX1jjE8l1i8pnQ61uXLLheEssucOJb4s9YNXTjA0MR1zduGoKy/mA1uW8s9vtib8eiXqL//PQXZ86yVGJ7N/rl28gPEo8F0RmfkXEvn6/sjPFpy6ch8uISdda9OtwXAEUjhiNNrx7hGmQzrnhncsWxoque+zl3P49BBf+fGbWW198GZzH5PTId6zPosBI9KxtivO6zgdDNHSO1qQGVL5sr7OjyoJL7PsaRmgyC1cWJ+ZTK1s2FDnp7GqJOH02pGJab736+Pcsqlu3q65d1y1iv7RKZ7c15GJoQLhJIKH32jm1i3LEj65MB3xAsafAKeBkyKyW0TeBE4AXZGfLTget4tlFcVZX5IamZimb3QqrRoMR43fy+hkMOVPF85eRKIzDMeNF9Tx3z62lRcOd/Fn/7I/pedOxKtNPbhdwpVpdFedz5LS+bPN2vrHmApq3ov2CkmyTQj3tPSxub4i4SNU80FEuOXCOn59tDuhmcCPXz9J/+gUX7l5/bzXvmddDatqSnk4Q8tS08EQf/yzd6gr9/H7H8hNunu8Su9pVf0qsAL4IvAFYKWqflVVc7Nglgf1VSVZX5JyUmrTqcFwzBTvpbgsdahjCI9LWJdCm+jbr1rJndeu4sevN9M5lJ3X7JWmHi5qrKQ8zoZ8ujxu17zZZjNdai1gzFgdKMUliQWMYEh5pzW1dNVcu3lTHWNTQV49Fj/FfnwqyHdfOs5719dw+cr5U4tdLuH2q1byxolejpxO/0z0H752kv3tg3z9I1vSbmeSqHhptR8XkY8DHyR8/sV6YJuIJPdR9DxTX5n9GcZM0V4GlqRqIwEj3nJKPIdPD7G2tizlgqI7rl6JKllpsDY8Mc3bLf1Z3b9wzLe057R2ON9SarPJ53GzqiaxTKmmrmFGJoMFV+EdyzVraygpcvPcPOm1P9nVQtfQBPfclPjxQJ+8YjlFbuHhN9KbZZweHOdvnj7M9o21fOiiZWk9VjLivUv8xqw/HwX+L2CviNycg7HlRWPkqNZUzhtIlFODsSKDM4xUi/cOdgxxwbLU15QvWFrOqppSntqf+YCx80Qv0yHlPesCGX/s2ebLNjvePUJ5sYeaLGVqna/W1SbWhHDPebDh7SgucnPdhgDPHOicsx3O5HSIv3/hGNtWVXPN2sSXSwN+Hzu21vPT3eltfv/5L95lMhjiv350S057csVbkvpSjD+3ATcCf5GzEeZYfWUxk9OhtPszxdPaP4bX7Uq7ayhENSBMYUlqeGKa1r4xLlia+tkOIsKOLct45Wh30n2F5vNaUw9FbkmqkjhVtfPMMI53j7A2UFZQDfMKwfo6Pyd6RpieJ/Hh7ZZ+yos9582S3i2b6mjrH+Pw6djB8OdvtdHWP8ZXbl6f9O/EHVetZHB8ml/sPZXS2F460sUv9p7iKzeuz/mMN+l1CFU9CSS0oCwiO0TkkIgcFZGvznHNjSKyR0T2i8gLkdtWiMhzInIgcvvvJTvOVM0U72VxWcqpwchEPr8TMFIp3ptpCZLGDAPg1q3LmA5pxk/pe6Wph8tWVmf0hL25BPw+uofmDrrHu0dsOSqGDXV+poLKyd7456bsaennkuVV500Ny02RJo3PxPidng6G+M7zR9naWMGNG5M/W+SatUtYW1vGwylUfo9PBfnaz/exJlDG3TeuTfr+6Uo6YIjIJmDedycRcQP3Ed4D2QzcLiKbZ11TBXwH+KiqbgE+FfnRNPAHqnohcA3wldn3zRZnX6E9ixvfbX1jGTsLwOdxU1HsSSm11gkYs7vUJuvS5VXUlft4al/mAsbA6BT72gdysn8BECj3MTYVZCRGXcn4VJC2/jHLkIohkUyp8akgBzuGCubApEQsrShma2MFz8ZIr/3Xd05xomeUe27akNKMU0S446qVvNncn3THhPtfaOJEzyh/fttWfJ7cZ5vF2/T+l+gT9iJ/Xgb+Ffj9BB77KuCoqh5T1UngEc49lOkO4FFVbQZQ1c7If0+p6puRr4eAA0Bjsn+5VDgHKWWzFqO1LzNFe45AeWrFe4dPD1HmdacdvFwu4dYty3j+cCdjk5kpSnrteA+q5GT/As7sBcU6F6O5dxRVLGDEsC6BgLG/fYBgSM+LDe9oN29aypvNffRGLU+HQsp9zx1l41I/H9ic+lG3n7h8OV6PK6kU2+PdI3zn+SZ+45IGrtuQm38Xs8WbYfw14VP2nD9/DdwNfAn4XAKP3Qi0RH3fyrlv+huBahF5PlLrcefsBxGR1cBlwOuxnkRE7hKRXSKyq6srvcNPINwB1udxZW1JanwqSPfwREZPGwv4fSllSR3sGGTD0vKMLBPs2LqM8akQL6Z5AI3j1aYeiotcOes7FJjZCzr3dTxzjrcFjNn8Pg/1lcVxA8aelnBL80LtITWXWzbVEVJ44fCZWcbT757m8OlhvnLT+rT+3VSXefnQ1mX87M22hGqoVJWvP7YPn9vF1z58YcrPm66ETtwDBoCPAL8g3JTwQAKPHevVnJ1y4AGuAD4M3Ap8TURmKlBExA/8FPiPqhpz7qaqD6jqNlXdVlub3lnFkeekoaqE9iz1SZppa74kkwEj+X5SqsqhjqG0l6McV61ZQmVJEU9lqIr11aYerly9JGfnB8TrJ3XcUmrjWj/Pca17WvppqCzOylkm2XRRYyUBv2+m6ls1PLtYVVPKhzNw5OodV69iaGKaX7w9/+b3L/ae4qUj3fzBBzbm9XWMtyS1UUS+LiIHgG8Tni2Iqt6kqt9O4LFbCRf9OZYD7TGueVJVR1S1G3gRuCTy/EWEg8WPVTWnrUjqK4uztiR1pgYjc2dThzvWJrck1TU0Qd/oVNIV3nMpcrt434VL+dWB02m3CukamuDQ6SGuzdH+BZzpWBureO941wgBvy9uN9/FbF2tn6au4TlT0d9u6T8v0mlnc7mEmzfV8uLhLqaCIV443MU7bQP8uxvX4cnAme5Xrq5mfZ2fH89TkzE0PsWf/+JdtjZW8PlrV6f9vOmI97c+CNwC/IaqXqeq/wtIZoF6J7BBRNaIiBf4DPD4rGseA64XEY+IlAJXAwckvJP0PeCAqv6PJJ4zIxqyWO3d2pf+wUmzBfw+Bsamkjp3++BMhlTm6jB3bF3G4Pg0r81TITsf5/652r8AZtrMx6pncVJqTWwblvoZnQxyavDcfzO9I5M0946elwEDwvsYg+PT7DrRx7efPUpDZTG/eVlmDvJyNr/fbulnX9vAnNf9zdOH6Rqe4Bsfuwh3nrPM4gWMTwAdwHMi8l0RuYUkjmhV1WngHuApwktYP1HV/SJyt4jcHbnmAPAksBd4A3hQVfcB7wU+D9wcSbndIyIfSuHvl5KGymI6h8az0lSvrX8Uj0tYmsFp5ZmzvRNfljp8Ov4pe6m4fkOAUq877eZqrzT1UO7zsDVDx4kmosjtYklZ7KW94z0jrA5kbka40KyvnXvj++3WSMHeebbh7bhuQwCv28VfP32IXSf7+PIN6zK6TPqJy5fj87jmrPze1zbAD149wWevXlkQQTfeHsbPVPXTwCbgeeA/AUtF5O9E5AOJPLiqPqGqG1V1nap+I3Lb/ap6f9Q196rqZlXdqqrfitz2sqqKql6sqpdG/jyRxt8zKfVVJYQ0XH6faa19Y9RXFWf0k8LMhm2cOoLZDnYMEfD7qMlA8aCjuMjNjRfU8vS7p9OqlH+1qZur1y7JyLQ/GbH2gobGp+gammBNIPXixoUuXmrt2y39iMBFBXwGRjx+n4er1y5h98k+An4fn74ys8cTV5YW8ZGLG3jsrbZzjgoIhpQ//vk+lpR5+cNbN2X0eVM177/IyP7Cj1X1I4T3IfYAMYvwFoqZczGysPGdyRoMh/Om353EDCOTG97Rbt2yjK6hCd5q6Uvp/u39Y5zoGc3a6XrxxDpbxDmW1TKk5lbj91FdWsTRGMe17mnpZ2Ndec6a42XDzZEivru2r8lKp907rl7JyGSQx/ecvcX7T28083ZLP3/84QupLCmM/bOkPsKpaq+q/r2qLtheUhB18l4WNr7DNRiZXd6oTbKfVDCkHD6d2Cl7ybppUx1Fbkl5WerVptzvXzgCft85dRjHeyylNhGxMqVUNbLhfX7OLhyfuGI5f/D+jXz+mtVZefzLV1axaVk5D79xpvK7a2iCv3ryINeureFjl+akBC0huZ3znyfqs1TtPTkd4vTQeMZnGM4BQIkW7zX3jjIxHcpKwKgoLuK96wM8tf/0nI3b4nmlqYfq0qKszH7mE6tj7fGuEURgVY3tYcQTK2C09I7RNzpVEGvv6agoLuLf37Ihay1qRIQ7rl7JvrZB9kb2fP7iiQOMTQX5849tLaj+ZRYwYvD7PFQUezJevHdqYAzVzGZIAZR6PZR63QnXYhyKtCPI1pvyji3LaO4d5cCp5Hr+qyqvNnVz7bqavPQcCpSfexjV8e5hGipLCvrQn0Kwvq6cvtGps3qa7TnPN7xz6WOXNVJS5Obh15t5tamHR99q467ta2f2hwqFBYw5NFSVZHxJaqYGI8MBA8JNCBMNGAc7hhAJn8mcDe/bvBSXwJP7k1uWau4dpX1gnGvzsBwF0a3iz8zUjveM2nJUAmJtfL/d0k9xkSsrM9mFpqK4iN+4pJ7H9rTzxz9/h+XVJUmds5ErFjDmEA4YmV2ScmowVmR4DwOSK9471DHEqiWlWZtiB/w+tq1ewtNJBoxXIvsX1+ZhwxvOPYxKVTneNWwBIwEzASPqfO89Lf1sbaikKMfZbuerO65exdhUkGNdI/z5bVtz0qU5WfZ/cg71lcUZX5I62TuC2yUsq8x8af98J8ZFO9SRnQ3vaDu2LONgx9BMW41EvNLUQ125j3W1+XmDdqq9ndexd2SSwfFpCxgJaKgsptTrnplhTAVD7Gs7P45kLRSXLK/k6jVLuO3Shpn26oXGAsYcGqpK6Budylj3VYD97YNsqPNn5RNXogFjfCrIiZ6RtM/AmM8HtoQ7eT6V4CzD2b94z7qavG3yze4nZU0HEyciZ52+d6hjiInpkAWMJIgIj9x1Dd/69KX5HsqcLGDMoaEqklqbwVnG/vZBtjRkJ8Ww1u+ld2SS4DwFc0c7hwlp9ja8HcurS7mosTLhgHGkc5ju4cm8pNM6amYVQFrASE50ppRT4X2pbXgnRUQKKitqNgsYc6ivjBTvZWgfo3NwnK6hCbZkqd1Fjd9HSKFvNP4+htNDamMGW4LM5dYtS3mruZ+OBAognfqLXDYcnK3I7aKqtOisGYbHJRnPaluo1tf5OTUwzvDENHua+1lS5mVFBrsym/yzgDGHhkqnFiMzM4z97eFU1q2N2ZlhxGvPHe1QxyBej4vVOagr2LF1GQBPvzv/LOOVpm6WV5ewYkl+6x2ii/eOd4+wsqY05y1KzlfOxndT5zBvt/ZzyfLKgv60bJJn/xLmsLTSh0jmlqT2t4e7UV5Yn51P9on2kzrYMcSGOn9O3gTX15WzrrZs3mWpYEh57Vhvzo5jjSe6n9Tx7hHW1NhyVKKcgPF2az9HOodt/2IBsoAxB5/HTcDvy9iS1L62QdYEyijP0pkKgfLEZhjZagkyl1u3LOO1Y730jcwdyA6cGmRgbCqv+xcOJ3kgFFJO9IzY/kUSVi0ppcgt/OytNlSxgLEAWcCIo6GyOHMzjFMDbM5iu+5A2fwBo390ktODExltaT6fHVuXEQwpzxzsnPOaV5q6gfzuXzicBoQdg+OMT4VYk6cU3/ORx+1idU0ZbzXbhvdCZQEjjkxVew+MTtHSO5a1DW+AihIPXrcrbj+pbByaNJ+LGitpqCyO24zw1aYe1tWWZfSMkFTVlvsYnpjmYKR9ii1JJcdZllpVU0p15FAqs3BYwIijvjJc7Z1KE71o+0+F9y+2ZimlFsLpePO1BzkUCRibslyDEU1E+MCWZbx4pIuRiXMPu58KhnjjeG9BzC7gTLX3G8fD7dlthpEcJ2BY/6iFyQJGHA1VxYxNBRkYm0rrcd6NZEhlc4YB8xfvHTo9RGVJEUsrMndoUiJ2bF3G5HT4TOTZ9rYOMDIZLIj9CzjT+XfniV5KitwsLc//rOd8MhMwbP9iQbKAEUdDhtqc72sboL6yOKOn28US68S4aE5LkFynOl65egk1Zd6Yy1KvRvYv8nFgUixOevLe1n5W1ZTmpWvu+eyqNUvYUOfnpgtq8z0UkwUWMOKoz9BBSuEK7+wvA9XEaUCoqhzuGMrphrfD7RLed+FSnj3YycT02a1WXmnq4cL6CpYUyHq3EzCmgspaW45KWn1lCb/8/RtYW1tYbblNZljAiKNx5qjW1APG2GSQpq5hNmdx/8LhdKyNtefS1j/G0MR03lpN79i6jOGJ6ZmOtBDua7X7ZF9B1F84nPYgYC1BjJnNAkYcAb+PIrfQnsbZ3gc6BgkpbM3BDCPg9zIZDDE4du7m8pkN7/wEjPesr8Hv8/BU1LLUW839TEyH8tbOPBafx01Fcfj86dWWIWXMWSxgxOFyCUsritNaktrfFs6Q2pKlliDRnPbcXTH2MQ6djvSQylPA8Hnc3LSpjl++e3qmQeKrTd24BK5auyQvY5qL8zrakpQxZ7OAMY+GqpK0qr33tw9SXVpEQxbOwJitJk7x3qGOIRqrSqjIUqV5InZsWUbPyCS7TvQC4f2Li5ZX5XVMsTj7GGsCtg5vTDQLGPNIt9rbaWmei8wkJyU01sb3oY4hNi7N7xvgjRfU4vW4eHJ/ByMT0+xp6S+o/QtHbbmPimIP1aWFFciMybesBgwR2SEih0TkqIh8dY5rbhSRPSKyX0ReiLr9IRHpFJF92RzjfOqrSugYGJ/3nIlYpoIhDnUM5SRDCubuWDsVDNHUNZz1Q5PmU+bzsH1DgKf3n2bniV6mQ1qQAePL29fxzU9cbJ1WjZklawFDRNzAfcAHgc3A7SKyedY1VcB3gI+q6hbgU1E//kdgR7bGl6iGqhKmQ5rw8afRjpweZjIYysn+BUB1qReXnBswjnWNMBXUvG14R7t1yzLa+sf47kvHKHIL21YV1v4FwEXLK/nQRfX5HoYxBSebM4yrgKOqekxVJ4FHgNtmXXMH8KiqNgOo6kyHOlV9EejN4vgS0pBGLca+SEvzXM0w3C5hSdm51d7Ohne+Umqjve/Cpbhdwq+P9nDZiuqCPOjeGBNbNgNGI9AS9X1r5LZoG4FqEXleRHaLyJ1ZHE9K0qn2frd9kDKvO6cN7MLV3mfvYRzqGMTtkoLI+qku83L1mvCsolD6RxljEpPNgBFrAXj2RoAHuAL4MHAr8DUR2ZjUk4jcJSK7RGRXV9e5vYrS5Zy8l0rx3v72AS6sr8hpe4lY/aQOdQyxNlCGz1MYn+adk/jeu74w+kcZYxLjyeJjtwIror5fDrTHuKZbVUeAERF5EbgEOJzok6jqA8ADANu2bUuvrWwMFSUeSr1u2pJckgqFlHfbB/nUthXzX5xBAb+Xk80jZ912sGOISwuoGdxnrlzJsopirlxdne+hGGOSkM0Zxk5gg4isEREv8Bng8VnXPAZcLyIeESkFrgYOZHFMSRORlGoxTvSMMDIZzOqhSbEE/L6zjmkdnpimtW+sIDa8HV6Piw9sWWZZSMacZ7IWMFR1GrgHeIpwEPiJqu4XkbtF5O7INQeAJ4G9wBvAg6q6D0BE/gl4FbhARFpF5N9ka6zzqa8sTnpJal+OWprPFij3MTYVnDl74vDMhnd+U2qNMee/bC5JoapPAE/Muu3+Wd/fC9wb4763Z3NsyWioLOHAqaGk7rO/fQCv28WGutx+sq8pO1O8V+bzzPSQykeXWmPMwmKV3gloqCqhe3jinNbc8bzbPsjGZX68nty+xIFZ/aQOdQxR6nWzvLokp+Mwxiw8FjASUF8VrsU4PZBY8Z6qsq9tIKtHss6ldla198GOQTYuLbeDgIwxabOAkQAntTbRTKlTA+P0jU7lfP8Czm4Poqoc6hgqqA1vY8z5ywJGAhoiM4xEN773RVqa5+LQpNmWRO1hdA1P0Dc6VRAV3saY858FjATUzxTvJZZau799EJfAhfW5f6P2elxUlhTRPTxhG97GmIyygJGAEq+bJWXehJek9rcPsLbWT6k3q0locwq3B4kKGDbDMMZkgAWMBNVXFnMq4YAxmJMjWefiFO8d7Bgi4PdRE9nXMMaYdOTnI/B5qL6yhJbe0Xmv6xme4NTAOFvysH/hCJT7ONA+yPh00Da8jTEZYzOMBDVWJXby3v48VXhHC5R56Rqa4PDpIVuOMsZkjAWMBNVXlTA0Ps3Q+FTc684EjDzOMPw+hiamGZ8K2Ya3MSZjLGAkqL7SSa2Nnym1r32A5dUlVObxPGin2htsw9sYkzkWMBLUOHOQUvxlqXfbB/NS4R3NKd4TgY02wzDGZIgFjATVV81fizE0PsXx7pG87l8A1PjDxXurlpTaEajGmIyxgJGgpeU+XBJ/huF0tN3amN8ZhtNPymYXxphMsoCRII/bxdKK4rhne+9vD7cEyfcMo7bcR5Fb8h64jDELi9VhJGG+g5T2tQ0S8PuoqyjO4ajOVVzk5qe/+x7W1frzOg5jzMJiM4wkNFSVxF2S2t8+wNbGwjjZ7uLlVZT57POAMSZzLGAkoaGqhFMD46jqOT8bnwpypHM478tRxhiTLRYwklBfWczEdIjekclzfnb49BDBkOY9pdYYY7LFAkYSGmZqMc7d+C6ECm9jjMkmCxhJcE7ei9VTal/bAOXFHlYssbOzjTELkwWMJDhne8fa+N7fPsiWhgpE7OxsY8zCZAEjCTVlXrwe1znV3tPBEAdODdpylDFmQbOAkQQRoaGy+JwZxrHuESamQwWTUmuMMdlgASNJsWoxzlR42wzDGLNwZTVgiMgOETkkIkdF5KtzXHOjiOwRkf0i8kIy982H+sqSc5ak9rUN4vO4WBsoy9OojDEm+7JWCiwibuA+4P1AK7BTRB5X1XejrqkCvgPsUNVmEalL9L750lBVzOnBcaaDITzucLzd3z7AhfUVM98bY8xClM13uKuAo6p6TFUngUeA22ZdcwfwqKo2A6hqZxL3zYuGqhJCCqeHJgBQ1ZkMKWOMWciyGTAagZao71sjt0XbCFSLyPMisltE7kzivgCIyF0isktEdnV1dWVo6HObOXkvso/R0jvG0Pi0dYY1xix42exOF6sgYXYTJg9wBXALUAK8KiKvJXjf8I2qDwAPAGzbti3mNZk0U+0d2cfYVyAtzY0xJtuyGTBagRVR3y8H2mNc062qI8CIiLwIXJLgffPCmWE4mVL72wdwu8QOKzLGLHjZXJLaCWwQkTUi4gU+Azw+65rHgOtFxCMipcDVwIEE75sX5cVFlBd7Zpak9rcPsqHOT3GRHYVqjFnYsjbDUNVpEbkHeApwAw+p6n4RuTvy8/tV9YCIPAnsBULAg6q6DyDWfbM11mQ1VJbQFmlAuK9tkBs21uZ5RMYYk31ZPWFHVZ8Anph12/2zvr8XuDeR+xaKhqrwyXudg+N0D09YhbcxZlGwwoEU1EcOUtpnFd7GmEXEAkYK3AarIAAABtFJREFUGqtK6B2ZZPfJPgA2W4aUMWYRsICRAidT6lfvdrImUIbfzs42xiwCFjBSUB85SOnQ6SGbXRhjFg0LGClorDpzqp6d4W2MWSwsYKRgaaVv5mur8DbGLBYWMFLg87gJ+MNBwwKGMWaxsN3aFDVWFVPkFmr8vvkvNsaYBcACRoo+d80qxqeC+R6GMcbkjAWMFH1q24r5LzLGmAXE9jCMMcYkxAKGMcaYhFjAMMYYkxALGMYYYxJiAcMYY0xCLGAYY4xJiAUMY4wxCbGAYYwxJiGiqvkeQ8aISBdwEggA3XkeTiGw1yHMXocwex3C7HUIc16HVapam8gdFlTAcIjILlXdlu9x5Ju9DmH2OoTZ6xBmr0NYKq+DLUkZY4xJiAUMY4wxCVmoAeOBfA+gQNjrEGavQ5i9DmH2OoQl/TosyD0MY4wxmbdQZxjGGGMyzAKGMcaYhCyogCEiO0TkkIgcFZGv5ns8+SIiJ0TkHRHZIyK78j2eXBKRh0SkU0T2Rd22RER+KSJHIv+tzucYc2GO1+FPRaQt8nuxR0Q+lM8x5oKIrBCR50TkgIjsF5Hfi9y+qH4n4rwOSf1OLJg9DBFxA4eB9wOtwE7gdlV9N68DywMROQFsU9VFV5wkItuBYeAHqro1cttfAb2q+s3IB4lqVf2/8znObJvjdfhTYFhV/zqfY8slEakH6lX1TREpB3YDHwO+yCL6nYjzOvwWSfxOLKQZxlXAUVU9pqqTwCPAbXkek8kxVX0R6J11823A9yNff5/wP5QFbY7XYdFR1VOq+mbk6yHgANDIIvudiPM6JGUhBYxGoCXq+1ZSeEEWCAWeFpHdInJXvgdTAJaq6ikI/8MB6vI8nny6R0T2RpasFvQyzGwishq4DHidRfw7Met1gCR+JxZSwJAYty2M9bbkvVdVLwc+CHwlsjxhzN8B64BLgVPA3+R3OLkjIn7gp8B/VNXBfI8nX2K8Dkn9TiykgNEKrIj6fjnQnqex5JWqtkf+2wn8jPBy3WJ2OrKG66zlduZ5PHmhqqdVNaiqIeC7LJLfCxEpIvwm+WNVfTRy86L7nYj1OiT7O7GQAsZOYIOIrBERL/AZ4PE8jynnRKQssqmFiJQBHwD2xb/Xgvc48IXI118AHsvjWPLGeYOM+E0Wwe+FiAjwPeCAqv6PqB8tqt+JuV6HZH8nFkyWFEAkJexbgBt4SFW/kech5ZyIrCU8qwDwAA8vptdBRP4JuJFw6+bTwH8Bfg78BFgJNAOfUtUFvSE8x+twI+GlBwVOAF921vEXKhG5DngJeAcIRW7+fwiv3y+a34k4r8PtJPE7saAChjHGmOxZSEtSxhhjssgChjHGmIRYwDDGGJMQCxjGGGMSYgHDGGNMQixgGGOMSYgFDLOoiMhvioiKyCYRuSiqrXOviByPfP2ryLVbRORZETkcaYP9tUgBFCLyRRH5dn7/NsbklgUMs9jcDrwMfEZV31HVS1X1UsKVv38Y+f59IlISue2bqroRuAR4D/DvcjlYEfHk8vmMicd+Gc2iEWm89l7gJsLB4E/jXH4H8GtVfRpAVUdF5B7geeC+JJ93GPj7yPP2EQ5WXSLyb4G7AC9wFPh85Hn+kXBr8suAN0XkfxPuYFACjAFfUtVDIvJFwm253cBWwo3jvMDngQngQwu5etnkns0wzGLyMeBJVT0M9IrI5XGu3UL4kJkZqtoE+EWkIsnnLQPejHQQfoFwmw6AR1X1SlW9hPD5BP8m6j4bgfep6h8AB4HtqnoZ8HXgv0ddt5VwcLsK+AYwGrnuVeDOJMdpTFw2wzCLye2EP6lD+ICt24E357hWmLs9frL9dELA/458/SPA6Zi6VUT+G1AF+IGnou7z/6lqMPJ1JfB9EdkQee6iqOueixyIMyQiA8C/RG5/B7g4yXEaE5cFDLMoiEgNcDPhN2klvIyjIvKfNXZDtf3A9lmPsZbwcZZDkb3vVDnP94/Ax1T17cjy0o1R14xEff3nhAPDb0YOv3k+6mcTUV+Hor4PYf++TYbZkpRZLD5J+HzrVaq6WlVXAMeB6+a4/sfAdSLyPoDIJvjfAn+VwnO7Is8P4eWjlyNflwOnIucUfDbO/SuBtsjXX0zh+Y3JCAsYZrG4nTNt3x0/JfwGfg5VHSN87vOfiMghwks8O4HoVNovikhr1J/lczz3CLBFRHYTnuX818jtXyPcZvuXhPcp5vJXwF+IyK8Jz4yMyQtrb25MlonIsKr68z0OY9JlMwxjjDEJsU0xYzJERF4HfLNu/rzNLsxCYUtSxhhjEmJLUsYYYxJiAcMYY0xCLGAYY4xJiAUMY4wxCfn/AanDVMEq2/2mAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cv_res_mean_hkh = [(x, y.mean()) for x, y in cv_res_hkh]\n",
    "plt.plot([x[0] for x in cv_res_mean_hkh], [x[1] for x in cv_res_mean_hkh])\n",
    "plt.xlabel('ATOL_param')\n",
    "plt.ylabel('AUC-ROC')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best ATOL parameter according to the 5-fold cross validation is n_centers=23"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_parameter = np.arange(1, 25).astype(int)[np.argmax([x[1] for x in cv_res_mean_hkh])]\n",
    "best_model = Pipeline([('atol_layer', Atol(best_parameter, padding=0, ignore_hom_dim=False)), \n",
    "                       ('model', RandomForestClassifier(n_estimators=100))])\n",
    "best_model.fit(diag_train, y_train)\n",
    "test_accuracy = roc_auc_score(y_test, best_model.predict_proba(diag_test)[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7243498817966902"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Features for more complex Atol usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.base import BaseEstimator, TransformerMixin, ClusterMixin\n",
    "from atol import lapl_feats\n",
    "\n",
    "class AtolX(BaseEstimator, ClusterMixin, TransformerMixin):\n",
    "    def __init__(self, n_centers=5, cluster_model=None, method=lapl_feats, aggreg=np.sum, order=None,\n",
    "                 ignore_hom_dim=False, padding=None, external=True, n_jobs=None):\n",
    "        self.atol = Atol(n_centers, cluster_model, method, aggreg, order, \n",
    "                         ignore_hom_dim, padding, n_jobs)\n",
    "        self.external = external\n",
    "    \n",
    "    def fit(self, X, y=None):\n",
    "        if self.external:\n",
    "            ex = np.array([x[1] for x in X])\n",
    "            X = np.array([x[0] for x in X])\n",
    "        self.atol.fit(X, y)\n",
    "        return self\n",
    "    \n",
    "    def transform(self, X, y=None):\n",
    "        if self.external:\n",
    "            ex = np.array([x[1] for x in X])\n",
    "            X = np.array([x[0] for x in X])\n",
    "        Xt = self.atol.transform(X, y)\n",
    "        return np.concatenate([ex, Xt], axis=1)\n",
    "    \n",
    "    def get_params(self, **params):\n",
    "        atol_params = self.atol.get_params(**params)\n",
    "        atol_params['external'] = self.external\n",
    "        return atol_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_validate_with_ext(parameters, X, y, k_fold=5, **atol_params):\n",
    "    set_up = []\n",
    "    for parameter in parameters:\n",
    "        pipeline_model = Pipeline([('atol_layer', AtolX(n_centers=parameter, **atol_params)),\n",
    "                                   ('model', GradientBoostingClassifier(n_estimators=100))])\n",
    "        set_up.append((parameter, cross_val_score(pipeline_model, X, y, scoring='roc_auc', cv=k_fold)))\n",
    "    return set_up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "external_diags_train = [(diag, ex) for diag, ex in zip(diag_train, ext_feat_train)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "external_diags_test = [(diag, ex) for diag, ex in zip(diag_test, ext_feat_test)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cv_res_hkh = cross_validate_with_ext(np.arange(1, 25).astype(int), external_diags_train, y_train, \n",
    "                                     5, padding=0, ignore_hom_dim=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEHCAYAAABr66s0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXycZ3no/d812vfdWm3Lqxw7XuM4NlkICUkTAgkE6ElKCISmOel7oECXl5TTlgPvaZsW2k85bd5SllBoaSC0SQgQGmgSYhLsxI53x5Is7xpZi61tJI3Wuc8f8zzyeDySZns0M5rr+/n4E+mZZ0aPJqO55r6v+74uMcaglFJKxcKV6AtQSimV+jSYKKWUipkGE6WUUjHTYKKUUipmGkyUUkrFTIOJUkqpmGUm+gLmQ2VlpWlsbEz0ZSilVEp56623LhhjqsI5Ny2CSWNjI3v37k30ZSilVEoRkTPhnqvTXEoppWKmwUQppVTMNJgopZSKmQYTpZRSMdNgopRSKmaOBhMRuUNEWkSkTUQeC3F7mYg8KyKHRORNEbk64LbPishRETkiIk+JSK51/MPWcZ+IbHXy+pVSSoXHsWAiIhnAE8CdwFrgfhFZG3Ta54EDxpgNwIPAV6371gO/B2w1xlwNZAD3Wfc5AtwL7HTq2pVSSkXGyZHJNqDNGHPSGDMOfB+4J+ictcBLAMaYZqBRRKqt2zKBPBHJBPKBDuu8Y8aYFgevOybNnYO8/4nXGRydSPSlKKXUvHEymNQD5wK+b7eOBTqIf5SBiGwDlgINxhg38BXgLHAeGDDG/NzBa42b19sucuBcP8e7PIm+FKWUmjdOBhMJcSy4rePjQJmIHAA+BewHJkWkDP8oZhlQBxSIyAMR/XCRR0Rkr4js7enpifzqo9TR7wXg/MDovP1MpZRKNCeDSTuwOOD7BqypKpsxZtAY85AxZhP+nEkVcAp4N3DKGNNjjJkAngHeEckPN8Z83Riz1RiztaoqrNIyceHu8weTTg0mSqk04mQw2QOsEpFlIpKNP4H+fOAJIlJq3QbwMLDTGDOIf3pru4jki4gAtwLHHLzWuHH3azBRSqUfx4KJMWYS+CTwIv5A8LQx5qiIPCoij1qnXQUcFZFm/Ku+Pm3d9w3g34F9wGHrOr8OICIfEJF2YAfwUxF50anfIRrT01yDGkyUUunD0arBxpgXgBeCjn0t4OtdwKoZ7vsF4Ashjj8LPBvfK40P7/gUF4fHAejSkYlSKo3oDvg46hjwj0qyM12agFdKpRUNJnFkJ9831JfQNTiKzxe8eE0ppRYmDSZxZCffr1laxqTPTE95KaXUQqfBJI46+r1kuISNi0sBXdGllEofGkziyN3npaY4l4ayPADOWzkUpZRa6DSYxFF7v5e60lxqinMB6NLlwUqpNKHBJI46+r3Ul+ZRUZhDpkt0RZdSKm1oMImTKZ+hc2CUutI8MlxCdXEunToyUUqlCQ0mcdLtGWXSZ6i38iXVxTmagFdKpQ0NJnFi7zGpL/UHk9qSPA0mSqm0ocEkTuw9JnYwsae5jNGNi0qphU+DSZzYwaRuemSSy8j4FIOjk4m8LKWUmhcaTOLE3eelND+Lghx/7cyaEv/yYJ3qUkqlAw0mceK2lgXbpoOJruhSSqUBDSZx0hEcTIrtkYnugldKLXwaTOLAGIO7zzudLwF/Ah6gc2AsUZellFLzRoNJHAx6Jxken5quyQX+niaVhdl0DurIRCm18GkwiYP2/hGAy6a5wJ830ZIqSql0oMEkDuwNi3XBwaQ4V1dzKaXSggaTOOiwNyyWXTky0dVcSql0oMEkDtz9XnIyXVQUZF92vLYkj/6RCUYnphJ0ZUopNT80mMRBR/8o9aV5iMhlxy+t6NLRiVJqYdNgEgft/d4rprjAX1IF0CS8UmrB02ASB+4+L3UlVwaTS7vgdXmwUmph02ASo9GJKS4MjYUcmdToxkWlVJpwNJiIyB0i0iIibSLyWIjby0TkWRE5JCJvisjVAbd9VkSOisgREXlKRHKt4+Ui8gsROW79t8zJ32Eu9hRW8B4TgIKcTIpyM7WkilJqwXMsmIhIBvAEcCewFrhfRNYGnfZ54IAxZgPwIPBV6771wO8BW40xVwMZwH3WfR4DXjLGrAJesr5PmJn2mNhqdXmwUioNODky2Qa0GWNOGmPGge8D9wSdsxZ/QMAY0ww0iki1dVsmkCcimUA+0GEdvwf4jvX1d4D3O/crzM3eY9IQYpoLrCZZmoBXSi1wTgaTeuBcwPft1rFAB4F7AURkG7AUaDDGuIGvAGeB88CAMebn1n2qjTHnAaz/LnLsNwhDe78XkUvLgIPVakkVpVQacDKYSIhjwT1sHwfKROQA8ClgPzBp5UHuAZYBdUCBiDwQ0Q8XeURE9orI3p6ensivPkwd/V6qi3LJzgz9VNYU59IzNMbElM+xa1BKqURzMpi0A4sDvm/g0lQVAMaYQWPMQ8aYTfhzJlXAKeDdwCljTI8xZgJ4BniHdbcuEakFsP7bHeqHG2O+bozZaozZWlVVFc/f6zLuvtB7TGw1JXkYAz0eXdGllFq4nAwme4BVIrJMRLLxJ9CfDzxBREqt2wAeBnYaYwbxT29tF5F88W8rvxU4Zp33PPAx6+uPAT9y8HeYk7vfO2PyHXTjolIqPTgWTIwxk8AngRfxB4KnjTFHReRREXnUOu0q4KiINONf9fVp675vAP8O7AMOW9f5des+jwO3ichx4Dbr+4Tw+QznB7whlwXb7FxKl67oUkotYJlOPrgx5gXghaBjXwv4ehewaob7fgH4QojjF/GPVBLOnwsxs05z6chEKZUOdAd8DNqtPSb1paFXcgGU5meRk+nSkYlSakHTYBKD6T4mpfkzniMi2nFRKbXgaTCJgbvf3v0+88gE7I6LWlJFKbVwaTCJQUe/l+LcTIpys2Y9T0uqKKUWOg0mMfDvMZl5istWXZJL18AYPl/wnk2llFoYNJjEwN3vnTX5bqstzmV8ykfvyPg8XJVSSs0/DSYx8AeTmZcF26abZGkSXim1QGkwidLg6ASe0clZ95jYaqwujBpMlFILlQaTKHX0z97HJFDtdPteDSZKqcj1DY/z72+10+1J3vcQR3fAL2Tu6Q2LcweTysIcMlyiIxOlVFT+ZfcZ/vYXrbgEblhVxb2b67l9XTX52cnzFp48V5Ji3P3hB5MMl7CoKEc3LiqlotLa5aG6OIcPX7OYZ/e7+cwPDpCfncEd62p4/+Z6rl9ZSYYrVNeP+aPBJErufi/ZGS4qC3PCOr+mJFdLqiilonKiZ5i1tcX84W808fu3rWbP6V6eO+DmJ4fO88x+N4uKcrh7Yx0f2FLP2tpi/MXW55cGkyi5+7zUlebiCvPTQE1xLq1dHoevSim10Ph8hpM9Q1y/ogIAl0u4bnkF1y2v4AvvW8fLzd08u9/Nd3ad5puvnWJ1dSEf2NzAPZvqwsrpxosGkyjN1cckWE1JLjtbnev4qJRamNz9XsYmfaxcVHjFbblZGbxnfS3vWV9L3/A4Pzl8nmf3tfNX/9nMX7/YzPZlFTx25xo2Li51/Do1mESpo9/LTavC7+BYW5LL8PgUntGJOcuvKKWUra1nCIAVIYJJoLKCbD66fSkf3b6UMxeHeXa/m+f2u8nNypiPy9RgEo3xSR/dnrGw9pjY7CZZnQOjGkyUUmE70e0PJiurZg8mgZZWFPCZd6/m07eumrf8ie4ziULnwCjGhLfHxFZrbVzUFV1KqUic6BmivCCbsoLsuU8OMp+JeA0mUWjvHwGgIZKcSbFuXFRKRa6te4gVVQWJvow5aTCJgr1hMZKRyaJi/xJi3biolIrEiZ7hkMn3ZKPBJAod/f6AUBtGxWBbblYGFQXZOjJRSoWtd3ic3uFxVkSQL0kUDSZRcPePsKgoh5zMyFZJVBfn6shEKRW2E2Gu5EoGGkyi0NE/GtVmoFrtBa+UikA0K7kSRYNJFNz93oiWBdu0pIpSKhJt3UPkZLrCqgGYaBpMIuTzGdz93ohWctlqinPpHR5ndGLKgStTSi00J3qGWF5VGHbZpkTSYBKhi8PjjE/6oprmsjsu6uhEKRWOtp6hlFjJBRpMIhZJ6flgtdpxUSkVptGJKdr7vCmxxwQcDiYicoeItIhIm4g8FuL2MhF5VkQOicibInK1dbxJRA4E/BsUkc9Yt20UkV0iclhEfiwixU7+DsGi2WNiqymx9proyEQpNYeTPcMYg45MRCQDeAK4E1gL3C8ia4NO+zxwwBizAXgQ+CqAMabFGLPJGLMJuAYYAZ617vNN4DFjzHrr2B859TuEYrfrjS4BryMTpVR4ppcFp8BKLnB2ZLINaDPGnDTGjAPfB+4JOmct8BKAMaYZaBSR6qBzbgVOGGPOWN83ATutr38BfNCJi5+Ju99LUU4mJXmRF2sszMmkMCdTlwfP4iPf3M333zyb6MtQKuHauocQgWWVOs1VD5wL+L7dOhboIHAvgIhsA5YCDUHn3Ac8FfD9EeBu6+sPA4vjdL1hibSPSbCaEt24OJMezxivt13kqT3n5j5ZqQXuRM8Qi8vy562EfKycDCah1rKZoO8fB8pE5ADwKWA/MDn9ACLZ+APHDwPu8wngf4jIW0ARMB7yh4s8IiJ7RWRvT0/8mlK5+6LbY2KrLcnVnMkMmjsHATjU3s/FobEEX41SiXWiZzhlku/gbDBp5/JRQwPQEXiCMWbQGPOQlRt5EKgCTgWcciewzxjTFXCfZmPM7caYa/CPWE6E+uHGmK8bY7YaY7ZWVYXfxGou7n5vTBuItKTKzFo6/W2NjYFfHb+Q4KtRKnGmrFa9qZJ8B2eDyR5glYgss0YY9wHPB54gIqXWbQAPAzuNMYMBp9zP5VNciMgi678u4E+Arzl0/VcYGptkwDsR0zRXbUku3Z5RJqd8cbyyheHYeQ+VhTmUF2Tzy5buRF+OUgnTYbXqTZXkOzgYTIwxk8AngReBY8DTxpijIvKoiDxqnXYVcFREmvGPQj5t319E8oHbgGeCHvp+EWkFmvGPdL7t1O8QLJaVXLaaklx8Bi4MhZydS2stXYNcVVvETasq2Xn8Aj5f8KyoUumhza7JlUIjE0fb9hpjXgBeCDr2tYCvdwGrZrjvCFAR4vhXsZYQzzd7j0l9BKXng9lNss4PeKd3xCuYnPLR2jXEx9/RyNraYp470MEh9wCbFpcm+tKUmneptiwYdAd8RC7tfs+P+jHsAKJ5k8udvjjC+KSPpuoiblpdhQg61RWh7795lnueeF1rvy0Abd1DVETZqjdRNJhEwN3vJStDWFSUE/VjaPve0OyVXGtqiygvyGZDQym/bInfKrxEGhiZYGoepuxeau7m4Ll+vvmrk47/LOWsEz1DKTUqAQ0mEeno909NxVLBs7wgm+wMV1xGJn/w9EG+98aZuU9MAS2dHjJcMj1HfPPqKg6299M7nNq5pbHJKW768it8d9dpx3/W8S7/argnXjmhI98oHHEPJE2erq17KCUaYgWaMZiIyAMi8tEQx39HRH7L2ctKTu6+2JYFA4iIf+NijCOTjn4v/7GvnR/ubY/pcZLFsfMellcWTHevvLmpyloinNqjk/P9owx4J3jrTJ+jP2d0YoozvSPcu6WeKWP4q/9sdvTnLTTHuzy89+9f45n97kRfCr3D4/SNTKTUHhOYfWTyB8BzIY5/37ot7cS6+91WUxx7x8VXrHzC2x2DjE2m/hx5c+cgTTVF099vaCilLD8r5ae67DxbqzVqcEpb9xDGwK1rqnnkxuU8u9/teABbSPaf7Qfg9bbE729KpVa9gWYLJhnGmCv+AqxjkRemSnETUz66BkejaooVLB4lVV4+5g8m41M+jrgH5zg7uXlGJ2jv83JV7aUC0Bku4abVVexs7UmaqYdotPeNAP4KsOOTzu0tsoPV6upCfvfmFVQX5/ClHx9N6eduPh12DwCw68RFjEnsc9aWQq16A80WTLJE5IpxlogUAamzxCBOOgdG8ZnY9pjY7JIq0b5oveNTvNZ2gbvW1wKw/2xqfwK13wjXBIxMwD/VdXF4fPoPPRXZy8knfYaTF4Yc+zmtXUNkZQiNlQUU5GTy2J1rONg+wH/sWxjToE6zX2Odg6OcuTiS0Gs50T1EblZqtOoNNFsw+Rbw7yLSaB+wvv6+dVtasTcsxmOaq7o4l/FJH30jE1Hdf9fJC4xN+rhv22LqS/PYl+LBpNkqo9IUFExuWmUvEU7dqa72Pi/2eg27XIwTjnd5WF5ZSFaG/0/6no31bF5Syl/9Zwue0eheZ+licsrHsfOD3NzkL7u0++TFhF5PW88QyytTo1VvoBmDiTHmK8CPgFdF5KKIXAReBX5ijPnyfF1gsoilw2Kw2pJLGxej8dKxbgqyM9i2rJwtS8um53tTVfN5D0U5mVc8txWFOWyoL+GXram736S938uGhlIyXeJo3qSly8Oq6kvTIi6X8IX3rePC0BhPvBKyfJ2yHO8eYmzSx/s31VNZmJPwYHKiJ/VWcsEcS4ONMV8zxizFXxq+0Riz1Bjzj/Nzackllg6LwWLpBW+M4eXmbm5cVUVOZgabF5dyfmA06sCUDFo6PTTVFCFy5SexdzYt4sC5fvpSdImwu8/LssoCllUW0NLpzDTX8Ngk7X1eVldfPrLbtLiUD25p4MnXTnH6wrAjP3shONzun+Ja31DC9uXl7DqZuLyJ3ao31fIlMEcwEZGrReQ7wC+Bl0XkOyKyfl6uLMl0DHipLMyOS2+BmumRSeTB5Nh5D+cHRrnlqkUAbFlaBsC+M6k5OjHGcKxzkDW1RSFvf5e1RHhnCi4Rnpzy0Tk4Sn1pHk01RbR0ObNQwk7YBgcTgM/d0URWhvDnLxxz5GcvBIfdAxTmZLKsooAdKyroGhzjdILyJnar3hWLUmtZMMy+z+Qe/G1xX8XfQ+Rh6+tnrNvSSnsc9pjYqgpzcAl0RRFMXm72V+N/V5M/mKytLSYn05WyeZOOgVE8o5M01RSHvN1eIvxqCuZNOgdHmfIZGsryaKou4lyvl+GxybnvGKHAlVzBFhXn8slbVvGLt7tSfs+OUw67B1hbV4zLJWxf7i8HmKiprrae1CvwaJttZPIl4DZjzJPGmEPGmIPGmCfxV/L90vxcXvKI1x4TgMwMF1VFOVGNTF5q7mbj4lKqrJIu2Zku1teXpOyKrharjMpVNaFHJvYS4VdTcIlwe9+lKtOrrd/veHf8p7pauzxkZ7pYWhH60+wnbmhkaUU+X/rx29r6IIidfF9fXwLA8soCqooSlzc50T2ES6Bxhv+XyWzWpcHGmNPBB61jabXPxBhDR4xNsYLVlORFvAv+wtAYB871c+uaRZcd37yklCPu1Ny8eOy89al6hmACl5YIH+lIrSXCdp6toSyfJmsKyg6e8dTa5a/jlDHD6p+czAz+53uu4nj3EN9742zcf34qs5PvdjAR8Y9OErXf5ETPEIvLU6dVb6DZgsmEiCwJPigiSwlorZsOeofHGZ3wxWWPia02io6Lv2zpwRi4JSiYbFlSxviUj6Mdqbd5saXTQ31pHsW5M38+SdUlwvbIpLYklyXl+eRmuRxJwh/v8tAUYoor0G1rq7lhZSV/+4vWlF3M4AR7f8nVVjAB2LG8gm7PGKcSsGihrTv1CjzaZgsmXwD+S0Q+LiLrrWT8Q8DPgT+bn8tLDh39/jf9eE1zQXS74F9u7qK6OId1dZfnFy4l4VNvqqu5098Qazb2EuFXUqwkvbt/hEVFOeRmZeByCauri+K+PHhwdIKOgVFWhUi+BxIR/vS9axkam+Rvf9Ea12tIZYfbByjIzmB55aVppe3LywHYfbJ3Xq9lymc4dSG1+r4Hmm2fyXPAh4FbgH8Gvgu8C/hN67a04e73r+yI7zRXLp6xSYbCTMiOT/rY2XqBW9ZUX7GEtro4l/rSPPafS60VXWOTU5zsGb5is2IoqbhEuL3Pe9lodnV1ES1xDibHu2ZeyRWsqaaIB65bwvfeODNd8j/dHXYPsK6+5LINgssqC1iUgLyJu8/fqjcVk+8w9z6Tg8aYB40x1xhjtlhfH7SmutJG+/TcdxynuSJskrXndC9DY5NX5Etsm5eUsj/FRiYnuoeZ9BnWzLCSK9DNKbhE2B2UZ1tTU0SPZyyuZfXtsvNNYQQTgM/etprivCy++PzbCa9BlWjByXebnTfZPc/7TVKxu2KgufaZ7BCRD4nIIuv7DSLyb8Br83J1SaKjf5T87AxK8uK37qC6OLJg8tKxbnIyXVy/sjLk7ZuXlNExMJpSfSzsfRfBNblC2ZhiS4R9Pv+ijYayS105V08n4eM3OmntGiIvKyPsDzql+dn8/m2r2XXyIi8e7YrbdaSi4OR7oB0r/HmTk/OYN7H3Cy24YCIiXwaeBD4I/FREvgD8AniDGfq2L1Tu/hHqS/NC7tCOViQlVYwxvNTcxTtWVJCXHXqVx5Yl/l7pqbTfpPm8h+wMF8sq554jznAJN65KnSXC3Z4xJqbMZdNc9nRePPMmrV0eVi6KrI7Tb21bwurqQv78hbfTusVvqOS7LRH7TU70pF6r3kCzjUzuAjYbY+4HbgceA24wxnzVGJM6H3/jIJ57TGz2yCSckionLwxz5uIIt1xVPeM56+pKyM50pdR+k2Od/jfCzIzwGn6m0hJhO88WOGJYVJRDSV7WdGHLeGgNqskVjswMF1943zrO9Xr51mun4nYtqeaI+8rku62xIp/q4px5TcKnYnfFQLP9FXvtoGGM6QNajDHH5+eykktH/2hclwUD5GZlUJafFdbGRbt3SfCS4ED25sV9KVT0sWWWMiqh3LTaX9U1FZYIT+fZAj6EiAhNNfFb0dU/Mk63ZyzsfEmg61dWcvvaap54pS2qGnELwaH2K5PvtkTkTVKx73ug2YLJChF53v4HNAZ9nxZGxifpHR53pLdATUleWH/ILzV3saamaM5r2Ly4lMPuAUebMMVL3/A4XYNjYeVLbJWFOWxoKOGXKbBEOHD3e6Cm6iJaOz1xeYNqjWAlVyh/ctdaJqfSs8XvTMn3QDuWV9DjGeNEj/N5E7tVb6qu5ILZg8k9wN8E/Av+Pi3Ye0wcCSbFc5dUGfBOsOd036yjEtuWpWWMT/o4mgLTQPZUTzgruQLdvLqKA+f66R9J7iXC7X1eyguyyc/OvOz46poiPGOTMbdthku5l0inuWxLKvJ5+MZlPLPPPb0qLF3Mlny3zWfe5FLyPTX3mMDs+0xeDf4HeAK+TgvTfUziPM0FVkmVOd5Udrb2MOUz3HpVGMFkiX/zYir0N7H3OUQyzQX+/SY+AzuPJ75X92yClwXb7JFYPFZ0He/yUJCdEdMHnU/csAyXwHMH3DFfTyqZLfluW1qRT01x7rwEk1RfFgxzLA0O4ZuRnCwid4hIi4i0ichjIW4vE5FnReSQiLwpIldbx5tE5EDAv0ER+Yx12yYR2W0d3ysi2yL8HSJi11dyYmRSW5LLxeHxWWtqvdzcTVl+FpsWl835eDUludSV5KbEiq6WTg/lBdlUFeZEdL9Ni0spzc9K+qkud99IyOW6qxdZwSQOI4HWriFWVYfuAxOuysIcrl9ZyY8Pnk+rfSezJd9t/rxJObtP9jr+3LSlaKveQJEGk7BftSKSATwB3AmsBe4XkbVBp30eOGCM2QA8CHwVwBjTYozZZIzZBFwDjOAvhw/w18AXrdv+zPreMR39XjJcwqKiyN70wlFjrejqHhwLefuUz/BKSzfvalo0YxG/YJuXpEbnxWOdHtbM0BBrNhku4aZVVexM4iXCxpgZRyYl+VnUFOfSGoeRSWuXJ2TZ+Ui9b2MdZ3tHONie/NOj8RJq53so25dXcGFobHrk4JQTKdqqN1CkweSLEZy7DWgzxpw0xozj7x0f3AdlLfASgDGmGX+SP3j9663ACWPMGet7A9gT7SVARwTXFDF3v5ea4tywl69GYq4mWfvP9tE/MjHdCCscm5eU4u73JvUKHZ/P0Gp1V4zGzU1VXBgaT9rClhetwqAzbSRsqimKeXnwxaExLg6PR518D/Qb62rIznDx/AFH/5SSxuSUj7c7Zk++23as8OdNdjm8RLiteyilk+8w+6bF3xCRDwUeM8Y8JyIfEZHbwnjseuBcwPft1rFAB4F7rZ+3DX974Iagc+4Dngr4/jPAl0XkHPAV4I/DuJaozfQJMx6mS6rM8Mb/UnM3mdZmvXDZRR+Teb/J2d4RvBNTXBVh8t12aYlwck51XVrJlR/y9qaaItp6hmLqLRLrSq5AJXlZ3NxUxU8OdTCVpKO9eAon+W5bUp5PbYmzeRPv+BTufm9K50tg9pHJF/F3Vgz2EuE1xwo1Xgt+pT4OlInIAeBTwH4CytuLSDZwN/DDgPv8LvBZY8xi4LPAt0L+cJFHrJzK3p6e6PcluIOK9cVT9XR9rtC74F8+1s21jeURlXFZV1dMdoYrqfeb2Mn3aEcm00uEW5Nzv4l7jlpuq6uLGJ/0caY3+tawl7orxh5MwD/V1e0Z441TiWkKNZ/CSb7b7P0mbzi43+TUBX+r3gU7MgHyjTFX/LUaYzqBcNavtQOLA75vIGhKyhgzaIx5yMp/PAhUAYFbcu8E9hljAosIfQx4xvr6h/in065gjPm6MWarMWZrVVX4n+wDBfbwdkJRTiYF2Rl0DlyZMznXO0JLlyesVVyBcjIzWFdfnNTl6Js7PYjE9kZ48+oqaxow+ZYIt/dZVaZnmuayfu9Y8iatXR6KcjOpLo5PLu/dV1WTn53Bjw+ej8vjJbNwku+Bti8v58LQ+PTy3XizW/WmYt/3QLMFk1wRyQw+KCJZQDjvrnuAVSKyzBph3AdcttlRREqt28DfY36nMSZwIvx+Lp/iAn9Aeqf19S2AY7vyuz1jTPlM3Eup2ETE39dk8MqRid27I5z9JcG2LClL6s2Lzec9NFYUzFhnLBz2EuFfJeESYXe/l+LczBkbfq2qLkSEmPImx7uGaIpxJVegvOwMbltbzc+OnE/a1028HHYPsK5u7uS7bcdyf3FVp6a6UrlVb6DZgskzwDdEZPo3tL7+GpdGBjMyxkwCnwReBI4BTxtjjorIoyLyqHXaVcBREWnGP4Pd+10AACAASURBVAr5dMDPysffbz74Z/0O8DcichD4C+CRua4lWk7uMbHVlOSGTMC/dKybZZUFLI9iHnXLkjLGJv07fJNRS5cnop3voVxaIpx8U13+qdHQ+RLwl9JprCiIuqyKMYbWbs+cDbEidffGOvpHJnitLfme03iZTr43zD3FZVtcnkddSa5jdbraUrhVb6ArRh4B/gT438AZETmDPweyGH+O4k/DeXBjzAvAC0HHvhbw9S5mqEBsjBkBKkIcfw3/cmHHObnHxFZTnMeuE5d/uh4em2TXiYt8dEd0bWO2LL1UQXjj4tKYrzGeRsYnOX1xmPdvCl6LEZngKsLJtKSyvc/L4vKZgwnA6urCqPea9HjG6B+ZiMuy4EA3rqqiJC+L5w90cMuamYuKprJIku82O2/yamsPxpi4Vg8H/8hkZYon32H2HfCTxpjH8AeQj+PPVSwxxjxmjJmYp+tLKHtkUlea69jPqCnJocuaTrO93naB8SnfjI2w5lJbkkdNcW5SJuGPdw1hTPTJ90A3r67iwtAYbyfRCMzeYzJXf5GmmmJOXxiOqgR8PFdyBcrOdHHn1TX84u0uvOMLszR9JMn3QNuXV3BxOP55kymf4eSF4ZSuFmybbWnwvSJyL/7pp1XASmCriMT3FZzE3P2h6yvFU01JHlM+w4WhS0n4l5u7KcrJZGtjedSPu2VpaVIuD7ZXcs3V9z0cybhEeMA7wdDY5NzBpLoInyGqN6d4r+QKdPfGOobHp3i5OXme03iKNPluu7TfJL55E3efl/FJX0rX5LLNljN5X9C/u4E/BA6JyC3zcG0J96lbVvLkx6919GfUBnVc9PkMLzd3c9PqKrIzo98ouWVJGe19Xro9ybV58dh5D/nZGSyeJacQrqqiHNbXlyRV3iTcFs9NNf5PotHkTVq7PJTlZ1FZGP8mStctr6CqKIfnDy7MWl2RJt9tDWV51JfmxT0J39bj//+f6suCYfZprodC/LsHuBn4y3m7wgSqLcljk8M5h+Bd8Ec7Bun2jEW1iivQZqvo474zsU91vXHyIh//9puMjE/OffIcWjo9rK4uiluO4+amKvad7WNgJDlmXqc3LJbOHiyXVhSQneGKKm/ib4gVv5VcgTJcwl3ra3mlpYfB0eR4TuPFLjsf6RQX+PMm1zlQp+tEt7+8fapvWITIy6lglTWJXzP0NFcTtHHxpeYuRPxvkrG4ut6/eTHWqa4pn+FPf3SEX7b08Oz+2D6tGmNo7hyMeSVXoJubqvxLhJNkBZKdZ5trZJKV4WLFosKIqwcbY6aXBTvl7k11jE/6+PkC6xHf1jPE6ISPDRGs5Aq0fXkFvcPjHI9j3uREzxCVhdmU5qdmq95AEQcTEVkDhK5MqCJWnp9NdoaLTqvY48vN3WxeXEpFhNV0g+VkZrC2rjjmoo/P7nfT2jVEcW4mT752Kqbiij2eMfpGJuIaTDYtLqMkL4v/ejs53vjcfV7yszMozZ/781ZTdWHEGxc7B0fxjE3GfSVXoM2LS2koy+P5gwurVteh9uiS77YdVn+TXSfiN9XV1j0U1fL/ZDRbAv7HgZ0VrX+vAT8Ffn/+LnFhc7mERcU5dA546R4c5VD7ALfO0us9EluWlHHI3c9ElDWgRiem+Nuft7ChoYQv3rOOEz3D7Dwe/QjgmPXG2RRlTa5QMlzCPZvqeO5AB99/82zcHjda7X0j1JfmhTUFtbqmiI6B0Yimk+yRTLz3mAQSEd63sY7X2y5wcWjhfG6MNvluW1yeH/e8yYme1C/waJttZPIVLu+s+BXgUeAh4AHnLy191Jbk0jk4GtOu91C2LC1ldCL6zYv/uvsMHQOjPHbHGu5aX8eiohyefP101NfTYjfEiuPIBPztZ9+5uoo/fvZwwj9Nh7Ms2GY/D5GMTo47tCw42N0b65jyGV44vHDKq0SbfA+0fXkFb5zqjUv7g4tD/pH6QsiXQJidFoEB4L3AT/AXgDw2T9eXFqqLc+kcGOWlY93UleTG7c12y3QSPvK8yYB3gn94pY0bV1XyjpWVZGe6eHDHUna29kTd4rX5vIea4lzKCuI7P5yd6eJrD1zDtY3l/P4PDiR0yqs9gsKgdkCIJAnf2uWhsjCH8jg/h8HW1BSxalHhgqnVFUvyPdD25eVxy5vYveUX/MhERFaLyJ+JyDHgH/CXkxdjzLuMMf8wb1eYBmpLcukYGOW1tgvcctWiuK3SqS3Jpbo4h/3nIs+b/NOrJ+gfmeBzd6yZPnb/tiXkZLr49q9PR3U9zTH0MJlLXnYG3/rYVtbWFfP//Ns+ft02/zW7PKMTDHgnaAhz2XN9aR6FOZkRjUzi1RBrLiLC3RvrePN0Lx39oatapxI7+b6+IbYp1u3TeZPYX18Loe97oNmmuZrxN6Z6nzHmBmPM3wMLc1tsgtWU5DE+6WNkfIpb41jGQkTYsqQs4ja+XYOjPPn6Ke7eWHfZJ7mKwhw+sLmeZ/a10zccWbXeiSkfbd1DEfd8j0RRbhbfeWgbyyoKePi7e3lrnisnT9dyC7P8joiwurow7IKPPp/hePeQ41NctvdtrAPgJ4dSPxF/2Eq+r6+Pban/4vJ8Gsry4lKn60TPEHlZGdSVpG6r3kCzBZMPAp3AKyLyDRG5lQja9qrw2e17c7Nc0ztt42XLkjLO9Xrp8YSfSP27/zrOlM/wh7c3XXHbQ9cvY3TCx1N7Ikt2n7owzPiUL+75kmBlBdn8y29vY1FRDg99+02OdsxfK9q5+piE0lRTRGuXJ6y9C+5+LyPjU/MWTBorC9jQUJLwPFQ8HI4x+R7Inze5GHPexL+SqyCp6srFYracybPGmP8GrAF+ib8RVbWI/KOI3D5P15cW7L0mN6ysjHvl0MCij+E40TPE03vP8ZHrlrKk4srpmqaaIm5YWcl3f30molVi9qfvNXFcyTWTRcW5/OvD11GYk8mD33rTsT4UwS51WAw/mKyuLqJvZIKeMFZNHe+2y6jM3xz73RvrOOIe5KTDPdCdFo/ku2378gr6RiZo7Y6t9fKJnqEFk3yHMPaZGGOGjTHfM8a8F3+DqwPAY45fWRpprMgnO8PFe9bXxv2x19WVkJUhYe83+cqLLeRmuvjkLStnPOcTNzTSOTjKz450hn0dzecHyXTJvP3xNJTl868PX4cIPPDNNzgXQ1fDcLn7vWRnuqgsCH+P0KVGWXO/WbdY5zi5LDjYezfUIUJKJ+LjlXy3bV/ur5m3O4b9Jnar3oWSfIcINy0aY3qNMf9kjEmL2lzzpaIwh1//8S18YHNsZdlDyc3KYG1dSVgjk/1n+/jZkU5+56blVM6yafLm1YtYVlnAk6+dmvGcYC2dHlZUFcZUbyxSy6sK+Zffvo6R8Uk+8s036Bp0tk6Zu89LQ2leRJ9+7QUJdgHM2Rzv8lBdnBNRG+dY1ZTkcm1jOc8fdDvWttZp8Uq+2xrK8llcnhdT0ceTF/zVs9NqZKLmR2VhjiO1lgC2LCnlUPvsmxeNMTz+s2YqCrJ5+Mblsz6eyyU8dH0jB871hz195uRKrtlcVVvMdz6xjYtDYzzwzTfojXDhQCTa+0YibqRWUZhDZWF2WAUfW7s985YvCXT3xjpO9AwnVan/SFxKvsdnZAKwfVls+00W2rJg0GCSFrYsKWN0wkfz+ZnfsH7Z2sMbp3r5vVtXUZgzd8n9D25poMgqsTKXwdEJ3P1eR1dyzWbzkjK++bFrOds7woNPvuFYAcNINiwGWl1dREvX7NNcUz5/Ta5EBJP3rK8lwyUpO9Vl73xfVhm/N+7tyyvoH5mIusFZm92qtzL26tnJQoNJGti8xJ+E338u9CjC5zP81c+aWVKez/3bloT1mAU5mdy/bQk/O9I55z4EuwTIVfOQfJ/JjhUV/OMDW2g+7+ET394TlwrIgbzjU1wYGo+qK2dTTRHHuzyzfso91zvC2KRvXpPvtvKCbG5YWcmPD3ak5FTXISv5nhHHVVPX2XmTKKe6TvQMsaQ8n5zM1G7VG0iDSRqoL81jUVHOjDvhf3TQTXOnhz+4fXVEOY0HdyzFGMN3d52Z9bxma3okEdNcgW5ZU81X79vMvrN9/Pd/eYuxyfhtm7pULTjyT5pN1UWMjE9NrwYLxcmGWOG4e2Md7n5vxHuWEi3eyXfbdN4kyiT8ie6FtZILNJikhUubF69c0TU2OcXf/LyVdXXFvG9DXUSP21CWzx1X1/DUm2dn/aTf3OmhODeT2hLn2h+H664NtTx+7wZ+dfwCn/3Bgbh90m7v868WizRnAv6CjzB7WRW7fMd8ruQKdPu6anIyXTx/ILX2nMQ7+R5oR5R1uuxWvQspXwIaTNLGlqWlnO0duaw9MMD3dp+lvc/L5+5YE9Ua/E9cv4wB7wTP7Ju510lzp4c1NcWOLTCI1G9eu5hPvmslLxzunB5RxCrcPiah2KON2ZLwLZ2e6fIriVCUm8Utaxbx08PnmYyyCnUiOJF8t21fXsGAdyLsCga29r4Rq1WvBhOVguzOi4H7TTyj/mKO16+s4MZVlVE97jVLy9jQUMK3Xw/d68QYQ0unJ2HJ95ncub4GgL2n4zNt4+7zkukSFhVFPvoqzMmkoSxv1jclf3fFxL75vG9jHReGxqMqJWKM4Yh7YN5L2h9xD5Af5+S77TqrTtcTr7TxowNuDp7rD6vj5wlrA+iKBTYySczHHDXv1teXkOkS9p3t47a1/vpf39h5kt7hcT53x5qoRw0iwieuX8ZnfnCAncd7uLnp8vL57X1ehsYmE54vCbamppiinEzePN3L++Owv6e9z0ttaW7USd6m6qIZCz5OTvk42TPMO1fH1n0zVresWURhTibPH3RzQ5gfPobGJnluv5t/3X2G5k4Pd22o5Ynf2uLwlV5yyD3A1XFOvtvqS/O4cVUlPz18np8GlOovzc+isaKAxop8llYU0Fhp/beigLL8rIBWvQujwKNNg0mayM3KYF1d8XQSvtszyjd+dYq7NtSyoSG24nfvWV/LX7xwjCdfP31FMGmZxzIqkchwCVuWlrH3dOwF+8BaFjxH3/fZrK4p4tXWHsYnfVcsgjh9cYTxKV/C8iW23KwMbl9bzX8e6eT/e//Vs65EertjkO+9cYbn9rsZHp9ibW0xV9UWs38ei2/ayfff2rbUsZ/xL799HaMTU5ztHeH0hWHOXBzh9MVhTl8cZs/pPn50sIPAtFxxbiYulyyYVr2BNJikkc1LyvjBnnNMTvn4+5famJjyhSzmGKnsTBcfe0cjX36xheNdnsve9Oyd3ck2MgG4trGMr/y8h/6R8Zj/sNv7RrhxVfQjhzU1RUz6DKcuDF/xXB3vmv+aXDN536Y6ntnvZmfrhekRrm10YooXDp/nX3efYd/ZfnIyXbx3Qx0PbF/CpsWlfONXJ/mLF5q5ODQWc1vqcDiZfA+Um5XB6uqikCvtxianONfr5czFYU4FBJtrG8sdvaZE0GCSRrYsLeOff32a/zzayVNvnuW+bYtZFocqquDvdfJ/XjrOk6+f5i/vXT99vLnTw+LyxCWOZ7PV+oN+60xfTK2Sxyan6PaMRZV8twU2ygoOJq1dQ4gkx27pG1ZWUpafxfMHO6aDyakLw/zbG2f44Vvt9I9MsLyygD+56yo+dE3DZUHaXp572D1wxQjWCU4m38OVk5nBykWFSfH/zmmO/oWLyB3AV4EM4JvGmMeDbi8DngRWAKPAJ4wxR0SkCfhBwKnLgT8zxvydiPwAsD9OlwL9xphNTv4eC8Xmxf7prD9+5jBZGS5+79ZVcXvs8oJs7t3i73Xy//5G03Q3RXslVzLatLiUrAzhzdO9MQWT8/2jGBN+H5NQllcVkOESf95k4+W3tXZ7WFyWT3524gNyVoaLO9fX8uw+N88f7ODpPed4re0CmS7h9nXVPHDdUnasqAiZg7ODyZF5CiZOJt/VlRxbzSUiGcATwJ3AWuB+EVkbdNrngQPGmA3Ag/gDD8aYFmPMJitIXAOMAM9at/23gNv+A3jGqd9hoWkoy6OqKAfP6CQP37gsqpVHs3no+mWMTfr4tzf9vU5GJ6Y4dWHY8R4m0crNymB9fUnMK7pi2bBoy8nMYFllQci9Jq2d89NdMVx3b6zDOzHF7z21n5M9Q/zBbav59WO38P9/5BresbJyxsUcxblZLKss4FD7/PSYOexg8l1dycmlwduANmPMSWPMOPB94J6gc9YCLwEYY5qBRhEJ/oh4K3DCGHPZNmvxv2J/E3jKiYtfiESE65aVU16QzSM3zV7MMRqrq4u4cVUl3911erqz4pTPJO3IBODaxnIOtfczOhH9bvhommKF0lRTNL1gwTY+6ePUheGEJ98DbWss5/PvWcM3H9zKrz53C5+6dRWLisP7YHJ1fQlH3M4Hk8kpH287sPNdzczJYFKPv2+8rd06FuggcC+AiGwDluLvmRLoPkIHjBuBLmPM8bhcbZr43++/mh9/6gaKcp0pY/6J65fRNTjGC4fPT++bSMbku+3axnImpgwHz4XX7yWU9r4RXHKpyVm0mqqLONs7clk1gdMXh5n0mem+J8nA5RIeuWkF715bHfGn/g31JXQMjDq+32S+ku/qEieDSahXWfCutseBMhE5AHwK2A9M/yWJSDZwN/DDEI91P7OMSkTkERHZKyJ7e3p6Ir32Bas0Pzumuf25vHN1FcutXifN5wfJyXTRGKJjY7K4Zql/M+feGJastvd7qSnOJSsjtj8nOwl/PKCCsD1SSfSGxXgJTMI7KRmS7+nGyWDSDiwO+L4BuKywjzFm0BjzkJX/eBCoAgJrmt8J7DPGdAXeT0Qy8Y9oApP0lzHGfN0Ys9UYs7WqKrGbvdKJ3evkYPsAzx/sYFV1IZkxvsk6qawgm1WLCtkTw36T9j5vVDW5gtm5pcCpruNdHlyycJoorav3jxQOO5w30eT7/HPyr3wPsEpEllkjjPuA5wNPEJFS6zaAh4GdxpjADjwzjT7eDTQbY9oduG4Vo3u3NFCcm0m3Zyyp8yW2rY3lvHW6j6koGx25+7wxJd9ti8vzyc1yXZaEb+0aYmlFAblZC6NUuZ2Ed3xk4h5gXV2xJt/nkWPBxBgzCXwSeBE4BjxtjDkqIo+KyKPWaVcBR0WkGf8o5NP2/UUkH7iN0Ku1ZsqjqCRg9zoBknYlV6Bty8rwjE1ekfwOx+SUj87B0bhMHWa4hFWLii4r+OjvrriwPl2vdzgJbyff19fHVtlBRcbRhevGmBeAF4KOfS3g611AyM0OxpgRoGKG2z4ev6tUTvj49Y3sPnkx4fWkwrF1qX/z4t4zvayti2wk1Tk4ypTPxLySy9ZUU8TOVn+Ob3RiitMXhrlrfW1cHjtZrK8v4fmDHVwYGqPSgZ3wmnxPjOSdzFYprbYkjx998oakWtI6k4ayPGqKc3nzVOR5E3tZcDxyJuBf0dXtGaNveJyTPcP4TOJ6mDjF6ST8IU2+J4QGE5X2RIRrl5Wz53RvxM2y2vti37AYKLBR1vFua2n1AgsmdhL+iENJ+DdO9lKWn8VyTb7PKw0mSuEv+tg1ODZr69xQ7N3v8eoi2RTQKKul00OmS+JWPy1ZFOdmsdzBJPzukxe5bllFVM3eVPQ0mCjF5XmTSLT3jVBVlBO31VbVxTmU5GXR3OmhtWuIxsqCK0rSLwRO7YQ/1zuCu9/L9uULrypvslt4r1KlotBUU0RRbiZvnops86K73xu35Dv4p9zsRlnHF+BKLtt6ayd8cBvpWO0+eRGAHSui6xyqoqfBRCn8y3KviaJZVnufN+4VBVbXFHLs/CBne0dC9shYCJxKwu8+2Uu5tRFVzS8NJkpZrm0s53j3EH3D42Gd7/MZzvePxi35bmuqKWZ4fApjWMDBxJkkvD9fUq75kgTQYKKUZWuEdbp6hsYYn/LFbVmwLXD11kKd5ipyIAl/KV8ScnuacpgGE6UsGxeXkp3hCnuqq71vBIi99HwwO4BkZ7hYWrGwVnIFurq+JK7BxM6XaDBJDA0mSllyszJY31ASdtHH6T0mcc6ZlOZnU12cw/KqgpgrESez9fUlnI9jEl7zJYm1cF+pSkVha2MZh90DYTXLao/z7vdAv33DMj66Y2ncHzeZrG+IbxJe8yWJpcFEqQDbrGZZB8JoluXu91JekO1Ib/ZHblrBR65b2MFkXV38kvCaL0k8DSZKBZhulhXGVJcTy4LTiZ2EPxSHkYnmSxJPg4lSAUrzs1ldXcie03Ov6HL3jcQ9+Z5u1jfEZyf8rpMXNV+SYBpMlApybWM5+87M3izLGIO7X0cmsYpHEt4YwxsnezVfkmAaTJQKcm1jOZ6xSZo7B2c85+LwOKMTPh2ZxCgeO+Hb+7y4+73sWKFTXImkwUSpIFsb7bzJzFNdl/qYxHf3e7qxk/Cx9ITfpfmSpKDBRKkg9aV51Jbk8uYsSfjpZcE6zRWTotwsllfFthN+t+ZLkoIGE6WCiAjXNpazd5ZmWe5+/+53J/aYpJtYesLb+ZLty8sR0XxJImkwUSqEuZpltfd5KcrNpCQva56vbOGxk/A9nsiT8Ha+RKe4Ek+DiVIhbG30N1eaqS+8u88b92rB6cpOwkczOtF8SfLQYKJUCE3V/mZZM3Ve1A2L8bOurhiR6FZ0ab4keWgwUSoEl0vYurQs5OZFe4+JLguOj6LcLJZFUY5e8yXJRYOJUjPY2lhOW/cQvUHNsga9kwyNTWowiaP19SURLw/WfEly0WCi1Ay2LfPnTYLrdJ1zqI9JOltfX0LnYGRJeM2XJBdHg4mI3CEiLSLSJiKPhbi9TESeFZFDIvKmiFxtHW8SkQMB/wZF5DMB9/uU9bhHReSvnfwdVPpaX1/ib5YV1HnR3W/vMdEEfLysjyIJv/vkRSo0X5I04l872yIiGcATwG1AO7BHRJ43xrwdcNrngQPGmA+IyBrr/FuNMS3ApoDHcQPPWt+/C7gH2GCMGRORRU79Diq95WZlsCFEsywn+5ikq3X1JdNJ+HetmftP2hjD7hMX2b68QvMlScLJkck2oM0Yc9IYMw58H38QCLQWeAnAGNMMNIpIddA5twInjDFnrO9/F3jcGDNm3a/bqV9Aqa2N5RxuH8A7fqlZlrvPS352BmX5usckXgpzMllWWcChMPMm53q9dAyMsn15ucNXpsLlZDCpB84FfN9uHQt0ELgXQES2AUuBhqBz7gOeCvh+NXCjiLwhIq+KyLVxvWqlAmxbVsak7/JmWe19I9SX5ukn4jjbEMFOeO1fknycDCah/tKCa1M8DpSJyAHgU8B+YHL6AUSygbuBHwbcJxMoA7YDfwQ8LSH+qkXkERHZKyJ7e3p6YvpFVPq6ZsmVSXhdFuyMqyNIwtv5kpWaL0kaTgaTdmBxwPcNQEfgCcaYQWPMQ8aYTcCDQBVwKuCUO4F9xpiuoMd9xvi9CfiAyuAfboz5ujFmqzFma1VVVXx+I5V2SvKzaKouuqzoo7vfq/kSB4SbhDfGsPuk5kuSjZPBZA+wSkSWWSOM+4DnA08QkVLrNoCHgZ3GmMAmEvdz+RQXwHPALdb9VwPZwAUHrl8pAK5dVsa+M31MTvkYGpukf2RCS6k4wE7Cz5U30XxJcnIsmBhjJoFPAi8Cx4CnjTFHReRREXnUOu0q4KiINOMfhXzavr+I5ONfCfZM0EM/CSwXkSP4k/ofMzOVdlUqDq5tLGd4fIrmTs+lPiZaSiXuCnMyWR7GTnjNlyQnx5YGAxhjXgBeCDr2tYCvdwGrZrjvCHDFq8VaGfZAfK9UqZnZRR/3nO5lSbl/RKLTXM5YX1/C7pMz95EBzZckK90Br9Qc6kvzqCvJZe/pvukNi5qAd4adhO/2jIa8XfMlyUuDiVJhuHZZOXtO99Le5yU700VlQU6iL2lBmisJr/mS5KXBRKkwbG0sp9szxq4TF2kozcPl0k/FTpjeCd8+GPJ2O1+yY4XmS5KNBhOlwnBtYxngL/eh+RLnzJWE333yIpWF2ayo0nxJstFgolQYVi8qojjXv15F8yXOmqknvJ0vuU7zJUlJg4lSYXC5ZHpVly4LdtZMSfhL+RKd4kpGGkyUCtNWa6pLp7mctaGhFLgyCb/rpH9v8g5NviclDSZKhemWNYvIzXJNrzhSzpjuCR+UhN99slfzJUnM0U2LSi0ka2qKOfalO3S+3mEFIZLwmi9JfjoyUSoC+kY2P9bXl3DYfans/9neEc5rviSpaTBRSiWd9Q2ldA2OTSfhp/eXaL4kaWkwUUolneCd8P58SY7mS5KYBhOlVNIJTMJfqsdVrtOMSUyDiVIq6RTkZLKiqpDD7n7Nl6QIDSZKqaTkT8IPaP+SFKHBRCmVlK6uL6FrcIwfHzxv5UsKEn1JahYaTJRSSclOwr/WdkHzJSlAg4lSKinZSXjQKa5UoMFEKZWU7CQ8aDBJBVpORSmVtK5tLMM7PqX5khSgwUQplbT+511r+ey7JzVfkgI0mCilklZhTiaFOfo2lQo0Z6KUUipmGkyUUkrFTIOJUkqpmGkwUUopFTMNJkoppWKmwUQppVTMNJgopZSKmRhjEn0NjhORHuAMUAlcSPDlJAN9Hvz0efDT5+ESfS787OdhqTGmKpw7pEUwsYnIXmPM1kRfR6Lp8+Cnz4OfPg+X6HPhF83zoNNcSimlYqbBRCmlVMzSLZh8PdEXkCT0efDT58FPn4dL9Lnwi/h5SKuciVJKKWek28hEKaWUAzSYKKWUillaBBMRuUNEWkSkTUQeS/T1JJKInBaRwyJyQET2Jvp65ouIPCki3SJyJOBYuYj8QkSOW/8tS+Q1zocZnof/JSJu6zVxQETek8hrnA8islhEXhGRYyJyVEQ+bR1Pq9fELM9DxK+JBZ8zEZEMoBW4DWgH9gD3G2PeTuiFJYiInAa2GmPSamOWiNwEDAHfNcZcbR37a6DXGPO49SGjzBjzuURep9NmeB7+FzBkjPlKIq9tPolILVBrjNknHOxQQAAABDxJREFUIkXAW8D7gY+TRq+JWZ6H3yTC10Q6jEy2AW3GmJPGmHHg+8A9Cb4mNc+MMTuB3qDD9wDfsb7+Dv4/ogVthuch7Rhjzhtj9llfe4BjQD1p9pqY5XmIWDoEk3rgXMD37UT5ZC0QBvi5iLwlIo8k+mISrNoYcx78f1TAogRfTyJ9UkQOWdNgC3pqJ5iINAKbgTdI49dE0PMAEb4m0iGYSIhjC3tub3bXG2O2AHcC/8Oa9lDp7R+BFcAm4DzwN4m9nPkjIoXAfwCfMcYMJvp6EiXE8xDxayIdgkk7sDjg+wagI0HXknDGmA7rv93As/inAdNVlzVnbM8ddyf4ehLCGNNljJkyxviAb5AmrwkRycL/Bvo9Y8wz1uG0e02Eeh6ieU2kQzDZA6wSkWUikg3cBzyf4GtKCBEpsJJsiEgBcDtwZPZ7LWjPAx+zvv4Y8KMEXkvC2G+elg+QBq8JERHgW8AxY8zfBtyUVq+JmZ6HaF4TC341F4C1rO3vgAzgSWPMnyf4khJCRJbjH40AZAL/li7PhYg8BdyMv7R2F/AF4DngaWAJcBb4sDFmQSenZ3gebsY/nWGA08B/t/MGC5WI3AD8CjgM+KzDn8efL0ib18Qsz8P9RPiaSItgopRSylnpMM2llFLKYRpMlFJKxUyDiVJKqZhpMFFKKRUzDSZKJSER+biI1CX6OpQKlwYTpZLTx4GIgomIZDpzKUrNTYOJUnMQkUarRPc3rDLdPxeRvBnOXSki/yUiB0Vkn4issI7/kYjssWodfXG2xxWRDwFbge9Z5b/zROQaEXnVqqn2YsAu7V+KyF+IyKvAp0XkwyJyxPr5O+fpKVJKg4lSYVoFPGGMWQf0Ax+c4bzvWedtBN4BnBeR2637b8O/EeyagJpoVzyuMebfgb3AR4wxm4BJ4O+BDxljrgGeBAI3m5YaY95pjPkb4M+A37B+/t3x+uWVmosOi5UKzyljzAHr67eAxuATrFI19caYZwGMMaPW8dvxl67Zb51aiD+InA3ncYEm4GrgF/7qF2TgL75n+0HA168D/ywiTwPPoNQ80WCiVHjGAr6eAkJNc4WqUG0f/0tjzD9ddtBf8jvcxz1qjNkxw+MP218YYx4VkeuAu4ADIrLJGHNxhvspFTc6zaVUnFilu9tF5P0AIpIjIvnAi8AnrDLfiEi9iMzVJ8MDFFlftwBVIrLDun+WiKwLdScRWWGMecMY82fABS6vmK2UY3RkolR8fRT4JxH5EjCBv1Dgz0XkKmCXNU01BDyAfyQyk38GviYiXmAH8CHg/4hICf6/278Djoa435dFZBX+0cxLwMG4/FZKzUELPSqllIqZTnMppZSKmU5zKRUFEXkCuD7o8FeNMd9OxPUolWg6zaWUUipmOs2llFIqZhpMlFJKxUyDiVJKqZhpMFFKKRUzDSZKKaVipsFEKaVUzP4vBus3u6uRSnkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cv_res_mean_hkh = [(x, y.mean()) for x, y in cv_res_hkh]\n",
    "plt.plot([x[0] for x in cv_res_mean_hkh], [x[1] for x in cv_res_mean_hkh])\n",
    "plt.xlabel('n_centers')\n",
    "plt.ylabel('AUC-ROC')\n",
    "plt.savefig('cv_output.png')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_parameter = np.arange(1, 25).astype(int)[np.argmax([x[1] for x in cv_res_mean_hkh])]\n",
    "best_model = Pipeline([('atol_layer', AtolX(best_parameter, padding=0, ignore_hom_dim=False)), \n",
    "                       ('model', GradientBoostingClassifier(n_estimators=100))])\n",
    "best_model.fit(external_diags_train, y_train)\n",
    "test_accuracy = roc_auc_score(y_test, best_model.predict_proba(external_diags_test)[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9773049645390071"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## only using classical features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_tot = np.concatenate([x.reshape(1, -1) for x in df.chemistry], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1484, 20)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_tot.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = GradientBoostingClassifier(n_estimators=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_out = cross_val_score(model, X_tot, df.CT_TOX, scoring='roc_auc', cv=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9581783560774356"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(cv_out)"
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}