notebooks/advanced_01_hyperopt_on_bbbp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f623d875-4b23-4fa4-9115-83fd5be390a7",
   "metadata": {},
   "source": [
    "# Hyperparameter search on the BBBP data set\n",
    "\n",
    "This is a more advanced example of how to use MolPipeline in combination with Sklearn to make molecular machine learning experiments. In this notebook we show how hyperparameter optimization can be performed using MolPipeline and Sklearn's [GridSearchCV](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.GridSearchCV.html#sklearn.model_selection.GridSearchCV). We will train multiple Sklearn estimators (like RandomForest) with Morgan fingerprints and scaffold splits and select the best performing model.\n",
    "\n",
    "We selected the BBBP data sets from [MoleculeNet](https://moleculenet.org/) because it is a widely used data sets for benchmarking molecular property prediction methods and it contains some problematic molecules that are erroneous which highlights the error handling capabilities of MolPipeline. See [the blog post](http://practicalcheminformatics.blogspot.com/2023/08/we-need-better-benchmarks-for-machine.html) by Patrick Walters for details on the erroneous molecules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "da554c7b-2dc4-4efd-8774-a8a3a9628293",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import roc_auc_score\n",
    "from sklearn.model_selection import GridSearchCV, GroupShuffleSplit\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "\n",
    "from molpipeline import ErrorFilter, FilterReinserter, Pipeline, PostPredictionWrapper\n",
    "from molpipeline.any2mol import AutoToMol\n",
    "from molpipeline.estimators import MurckoScaffoldClustering\n",
    "from molpipeline.mol2any import MolToMorganFP\n",
    "from molpipeline.mol2mol import (\n",
    "    ElementFilter,\n",
    "    MetalDisconnector,\n",
    "    SaltRemover,\n",
    "    StereoRemover,\n",
    "    TautomerCanonicalizer,\n",
    "    Uncharger,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ae49e047-cf6c-457a-867a-e9613268a95b",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_JOBS = 16\n",
    "RANDOM_STATE = 123456"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab18103b-3ede-43d4-a909-446058b91e82",
   "metadata": {},
   "source": [
    "Let's start by getting the data set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "39da2bef-0aa7-4078-ae9e-bd66bafc964a",
   "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>name</th>\n",
       "      <th>p_np</th>\n",
       "      <th>smiles</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Propanolol</td>\n",
       "      <td>1</td>\n",
       "      <td>[Cl].CC(C)NCC(O)COc1cccc2ccccc12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Terbutylchlorambucil</td>\n",
       "      <td>1</td>\n",
       "      <td>C(=O)(OC(C)(C)C)CCCc1ccc(cc1)N(CCCl)CCCl</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>40730</td>\n",
       "      <td>1</td>\n",
       "      <td>c12c3c(N4CCN(C)CC4)c(F)cc1c(c(C(O)=O)cn2C(C)CO...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>24</td>\n",
       "      <td>1</td>\n",
       "      <td>C1CCN(CC1)Cc1cccc(c1)OCCCNC(=O)C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>cloxacillin</td>\n",
       "      <td>1</td>\n",
       "      <td>Cc1onc(c2ccccc2Cl)c1C(=O)N[C@H]3[C@H]4SC(C)(C)...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2045</th>\n",
       "      <td>licostinel</td>\n",
       "      <td>1</td>\n",
       "      <td>C1=C(Cl)C(=C(C2=C1NC(=O)C(N2)=O)[N+](=O)[O-])Cl</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2046</th>\n",
       "      <td>ademetionine(adenosyl-methionine)</td>\n",
       "      <td>1</td>\n",
       "      <td>[C@H]3([N]2C1=C(C(=NC=N1)N)N=C2)[C@@H]([C@@H](...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2047</th>\n",
       "      <td>mesocarb</td>\n",
       "      <td>1</td>\n",
       "      <td>[O+]1=N[N](C=C1[N-]C(NC2=CC=CC=C2)=O)C(CC3=CC=...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2048</th>\n",
       "      <td>tofisoline</td>\n",
       "      <td>1</td>\n",
       "      <td>C1=C(OC)C(=CC2=C1C(=[N+](C(=C2CC)C)[NH-])C3=CC...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2049</th>\n",
       "      <td>azidamfenicol</td>\n",
       "      <td>1</td>\n",
       "      <td>[N+](=NCC(=O)N[C@@H]([C@H](O)C1=CC=C([N+]([O-]...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2050 rows × 3 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                   name  p_np  \\\n",
       "0                            Propanolol     1   \n",
       "1                  Terbutylchlorambucil     1   \n",
       "2                                 40730     1   \n",
       "3                                    24     1   \n",
       "4                           cloxacillin     1   \n",
       "...                                 ...   ...   \n",
       "2045                         licostinel     1   \n",
       "2046  ademetionine(adenosyl-methionine)     1   \n",
       "2047                           mesocarb     1   \n",
       "2048                         tofisoline     1   \n",
       "2049                      azidamfenicol     1   \n",
       "\n",
       "                                                 smiles  \n",
       "0                      [Cl].CC(C)NCC(O)COc1cccc2ccccc12  \n",
       "1              C(=O)(OC(C)(C)C)CCCc1ccc(cc1)N(CCCl)CCCl  \n",
       "2     c12c3c(N4CCN(C)CC4)c(F)cc1c(c(C(O)=O)cn2C(C)CO...  \n",
       "3                      C1CCN(CC1)Cc1cccc(c1)OCCCNC(=O)C  \n",
       "4     Cc1onc(c2ccccc2Cl)c1C(=O)N[C@H]3[C@H]4SC(C)(C)...  \n",
       "...                                                 ...  \n",
       "2045    C1=C(Cl)C(=C(C2=C1NC(=O)C(N2)=O)[N+](=O)[O-])Cl  \n",
       "2046  [C@H]3([N]2C1=C(C(=NC=N1)N)N=C2)[C@@H]([C@@H](...  \n",
       "2047  [O+]1=N[N](C=C1[N-]C(NC2=CC=CC=C2)=O)C(CC3=CC=...  \n",
       "2048  C1=C(OC)C(=CC2=C1C(=[N+](C(=C2CC)C)[NH-])C3=CC...  \n",
       "2049  [N+](=NCC(=O)N[C@@H]([C@H](O)C1=CC=C([N+]([O-]...  \n",
       "\n",
       "[2050 rows x 3 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv(\n",
    "    \"https://deepchemdata.s3-us-west-1.amazonaws.com/datasets/BBBP.csv\",\n",
    "    usecols=lambda col: col != \"num\",\n",
    ")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3009a4ae-f519-4b39-9787-879486bacd05",
   "metadata": {},
   "source": [
    "Make a function to setup some pipeline elements for standardization that we will use later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6d6f5fd3-4ba6-42ce-a742-4fecffae1ed2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_standardization_elements():\n",
    "    \"\"\"Make a new list of pipeline elements for molecular standardization.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    list[tuple[str, molpipeline.mol2mol.Mol2Mol]]\n",
    "        List of tuples with the name and the pipeline element.\n",
    "\n",
    "    \"\"\"\n",
    "    return [\n",
    "        (\n",
    "            \"element_filter\",\n",
    "            ElementFilter(\n",
    "                allowed_element_numbers=[\n",
    "                    1,\n",
    "                    3,\n",
    "                    5,\n",
    "                    6,\n",
    "                    7,\n",
    "                    8,\n",
    "                    9,\n",
    "                    11,\n",
    "                    12,\n",
    "                    14,\n",
    "                    15,\n",
    "                    16,\n",
    "                    17,\n",
    "                    19,\n",
    "                    20,\n",
    "                    34,\n",
    "                    35,\n",
    "                    53,\n",
    "                ]\n",
    "            ),\n",
    "        ),\n",
    "        (\"metal_disconnector\", MetalDisconnector()),\n",
    "        (\"salt_remover\", SaltRemover()),\n",
    "        (\"uncharge\", Uncharger()),\n",
    "        (\"canonical_tautomer\", TautomerCanonicalizer()),\n",
    "        (\"stereo_remover\", StereoRemover()),\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3abb637-ac6b-4e7f-bd87-f676060d430a",
   "metadata": {},
   "source": [
    "### Make a clustering pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31cfbec2-9384-4d87-9861-f20e4941c9b0",
   "metadata": {},
   "source": [
    "We first setup a pipeline that will do the scaffold clustering for us using [Murcko scaffolds](https://www.rdkit.org/docs/GettingStartedInPython.html#murcko-decomposition).\n",
    "\n",
    "In this example we also use Molpipelines error handling to detect and filter out all invalid molecules before clustering. After clustering we will re-add fill-values(np.nan in this case) at the corresponding positions of the invalid molecules. This keeps the resulting list of cluster labels (also called 'groups') in the same shape as the input molecule list. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "13a5feab-d4ac-4896-91a1-8055342d5614",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0        993.0\n",
       "1       1001.0\n",
       "2        671.0\n",
       "3        828.0\n",
       "4        258.0\n",
       "         ...  \n",
       "2045     616.0\n",
       "2046    1028.0\n",
       "2047     324.0\n",
       "2048     715.0\n",
       "2049    1001.0\n",
       "Name: groups, Length: 2050, dtype: float64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# create an error filter and reinserter.\n",
    "error_filter = ErrorFilter()\n",
    "# the Reinserter is wrapped into a PostPredictionWrapper because we execute\n",
    "# the re-insertion step at the end of the pipeline, even after the clustering.\n",
    "error_reinserter = PostPredictionWrapper(\n",
    "    FilterReinserter.from_error_filter(error_filter, np.nan)\n",
    ")\n",
    "\n",
    "# set up the clustering pipeline\n",
    "pipeline_clustering = Pipeline(\n",
    "    [(\"auto2mol\", AutoToMol())]\n",
    "    + make_standardization_elements()\n",
    "    + [\n",
    "        (\"error_filter\", error_filter),  # filter out invalid samples before clustering\n",
    "        (\n",
    "            \"scaffoldClustering\",\n",
    "            MurckoScaffoldClustering(linear_molecules_strategy=\"own_cluster\"),\n",
    "        ),\n",
    "        (\n",
    "            \"error_reinserter\",\n",
    "            error_reinserter,\n",
    "        ),  # re-add fill values to invalid samples to keep the input shape\n",
    "    ],\n",
    "    n_jobs=N_JOBS,\n",
    ")\n",
    "\n",
    "# cluster the smiles by their murcko scaffold\n",
    "df[\"groups\"] = pipeline_clustering.fit_predict(df[\"smiles\"])\n",
    "\n",
    "# the 'groups' are the cluster labels\n",
    "df[\"groups\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e4b0659-957c-4d14-9e29-5096857fae1b",
   "metadata": {},
   "source": [
    "We can have a look at the clustering results by plotting the distribution of cluster sizes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4f8fa876-f26c-4728-9fe6-8ed9176c9812",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of cluster: 1033\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHHCAYAAAC2rPKaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABKMUlEQVR4nO3deVxU9f7H8feAsiqgiQuK4JIZppBr5p4amdfUssWscMnqiivdSusWZpZLZbZQVvemXtvN9LaJW5ZXM7dEJcx9S8UNAUVTYb6/P3owP0dQGTgTOL6ej8c8Yr7nO+d8ztdR3p3zPefYjDFGAAAAHsirtAsAAABwF4IOAADwWAQdAADgsQg6AADAYxF0AACAxyLoAAAAj0XQAQAAHougAwAAPBZBBwAAeCyCDlAKfvjhB9lsNv3www/F+vzYsWNls9mc2iIjI9W/f/+SF3cZu3fvls1m04wZMxxt/fv3V4UKFdy+7Xw2m01jx479y7ZXHGvWrNHNN9+swMBA2Ww2paSkFHtdhY351aB///6KjIws7TJwhSPo4IqzadMm9enTRxEREfLz81PNmjXVtWtXvfnmm5Zv66WXXtK8efMsX29Z8d1335XZwFCWa7ucc+fO6e6771ZGRoZee+01zZo1SxEREaVd1kWlpaVp7Nix2r17d2mXAliuXGkXALjip59+UqdOnVS7dm0NHjxY1atX1759+/Tzzz/r9ddf17Bhwyzd3ksvvaQ+ffqoV69elq7XHbZs2SIvL9f+3+W7775TUlKSS4EiIiJCp0+fVvny5V2s0DWXqu306dMqV67s/vO1Y8cO7dmzR++//74efvjh0i7nstLS0vT888+rY8eOZeoIyvvvvy+73V7aZeAKV3b/pQAK8eKLLyo4OFhr1qxRSEiI07LDhw9bsg1jjP744w/5+/tbsr6/iq+vr1vXn5ubK7vdLh8fH/n5+bl1W5dT2tu/nPzv4oXf0atNTk6OAgMDi/15d4dpXB04dYUryo4dO9SoUaNCf4FUrVrV6X1ubq5eeOEF1atXT76+voqMjNTTTz+tM2fOOPWLjIzU3/72Ny1YsEDNmzeXv7+/3n33XdlsNuXk5GjmzJmy2Wyy2WxOc2D279+vgQMHqlq1avL19VWjRo30wQcfFKjr999/V69evRQYGKiqVatq1KhRBWq4lOXLl6tFixby8/NTvXr19O677xba78I5OufOndPzzz+va6+9Vn5+frrmmmvUtm1bLVq0SNKf8x+SkpIkybF/+fN+8ueEvPLKK5o6dapjDNPS0i45X2Tnzp2KjY1VYGCgwsLCNG7cOBljHMsvNjfpwnVeqrb8tguP9Kxfv17dunVTUFCQKlSooM6dO+vnn3926jNjxgzZbDatWLFCCQkJCg0NVWBgoHr37q0jR44U/gdwge+//17t2rVTYGCgQkJC1LNnT23evNmxvH///urQoYMk6e6775bNZlPHjh0vuc7MzEyNGjVKkZGR8vX1Va1atfTQQw/p6NGjF/1Mx44dC11vYfNaPv30UzVr1kwVK1ZUUFCQGjdurNdff90xJnfffbckqVOnTo6xPv/PaP78+Y59rlixorp3765ff/21wHYrVKigHTt26Pbbb1fFihXVr1+/i9Z/4sQJjRw50rHPVatWVdeuXfXLL79cdF86duzo9H04/3X+9zEzM1MjR45UeHi4fH19Vb9+fU2aNKnA0aFLjQs8B0d0cEWJiIjQypUrlZqaqhtuuOGSfR9++GHNnDlTffr00eOPP65Vq1ZpwoQJ2rx5s+bOnevUd8uWLerbt68effRRDR48WNddd51mzZqlhx9+WC1bttQjjzwiSapXr54k6dChQ7rppptks9k0dOhQhYaGav78+Ro0aJCys7M1cuRISX+eYuncubP27t2r4cOHKywsTLNmzdL3339fpP3dtGmTbr31VoWGhmrs2LHKzc1VYmKiqlWrdtnPjh07VhMmTHDsQ3Z2ttauXatffvlFXbt21aOPPqoDBw5o0aJFmjVrVqHrmD59uv744w898sgj8vX1VeXKlS96KiEvL0+33XabbrrpJk2ePFnJyclKTExUbm6uxo0bV6T9zVeU2s7366+/ql27dgoKCtKTTz6p8uXL691331XHjh31448/qlWrVk79hw0bpkqVKikxMVG7d+/W1KlTNXToUH322WeX3M7ixYvVrVs31a1bV2PHjtXp06f15ptvqk2bNvrll18UGRmpRx99VDVr1tRLL72k4cOHq0WLFpf88zp58qTatWunzZs3a+DAgWratKmOHj2qr776Sr///ruqVKlStEG7iEWLFqlv377q3LmzJk2aJEnavHmzVqxYoREjRqh9+/YaPny43njjDT399NO6/vrrJcnx31mzZikuLk6xsbGaNGmSTp06pXfeeUdt27bV+vXrnYJIbm6uYmNj1bZtW73yyisKCAi4aF2PPfaYvvjiCw0dOlRRUVE6duyYli9frs2bN6tp06aFfuaZZ54pcCrwww8/1IIFCxz/o3Pq1Cl16NBB+/fv16OPPqratWvrp59+0pgxY3Tw4EFNnTq1SOMCD2KAK8jChQuNt7e38fb2Nq1btzZPPvmkWbBggTl79qxTv5SUFCPJPPzww07t//jHP4wk8/333zvaIiIijCSTnJxcYHuBgYEmLi6uQPugQYNMjRo1zNGjR53a77vvPhMcHGxOnTpljDFm6tSpRpL5/PPPHX1ycnJM/fr1jSSzdOnSS+5vr169jJ+fn9mzZ4+jLS0tzXh7e5sL//pGREQ41RodHW26d+9+yfXHx8cXWI8xxuzatctIMkFBQebw4cOFLps+fbqjLS4uzkgyw4YNc7TZ7XbTvXt34+PjY44cOWKMMWbp0qWF7ndh67xYbcYYI8kkJiY63vfq1cv4+PiYHTt2ONoOHDhgKlasaNq3b+9omz59upFkunTpYux2u6N91KhRxtvb22RmZha6vXwxMTGmatWq5tixY462DRs2GC8vL/PQQw852vL3c/bs2ZdcnzHGPPfcc0aS+fLLLwssy6+xsPHp0KGD6dChQ4HPxMXFmYiICMf7ESNGmKCgIJObm3vRGmbPnl3on8uJEydMSEiIGTx4sFN7enq6CQ4OdmrP/w6MHj36Env7/4KDg018fPwl+1y4LxdasWKFKV++vBk4cKCj7YUXXjCBgYFm69atTn1Hjx5tvL29zd69e40xRRsXeAZOXeGK0rVrV61cuVJ33HGHNmzYoMmTJys2NlY1a9bUV1995ej33XffSZISEhKcPv/4449Lkr799lun9jp16ig2NrZINRhjNGfOHPXo0UPGGB09etTxio2NVVZWluPw+3fffacaNWqoT58+js8HBAQ4jhBdSl5enhYsWKBevXqpdu3ajvbrr7++SLWGhITo119/1bZt24q0X4W56667FBoaWuT+Q4cOdfycf7Tr7NmzWrx4cbFruJy8vDwtXLhQvXr1Ut26dR3tNWrU0P3336/ly5crOzvb6TOPPPKI06mwdu3aKS8vT3v27Lnodg4ePKiUlBT1799flStXdrQ3adJEXbt2dXznXDVnzhxFR0erd+/eBZZdeAuB4ggJCVFOTo7jlKUrFi1apMzMTPXt29fpe+7t7a1WrVpp6dKlBT7z97//vch1rVq1SgcOHHC5LklKT09Xnz59FBMTo7ffftvRPnv2bLVr106VKlVyqrlLly7Ky8vTsmXLHNsv7rjgykLQwRWnRYsW+vLLL3X8+HGtXr1aY8aM0YkTJ9SnTx+lpaVJkvbs2SMvLy/Vr1/f6bPVq1dXSEhIgV9oderUKfL2jxw5oszMTL333nsKDQ11eg0YMEDS/09G3bNnj+rXr1/gF9Z1111XpO2cPn1a1157bYFlRfn8uHHjlJmZqQYNGqhx48Z64okntHHjxqLsooMr4+Ll5eUUNCSpQYMGkuTWy5aPHDmiU6dOFTom119/vex2u/bt2+fUfn5wlKRKlSpJko4fP37R7eR/Zy62naNHjyonJ8fl+nfs2HHZ07AlMWTIEDVo0EDdunVTrVq1NHDgQCUnJxfps/kh+ZZbbinwXV+4cGGBCwDKlSunWrVqFWndkydPVmpqqsLDw9WyZUuNHTtWO3fuLNJnc3Nzdc899ygvL09ffvml00T8bdu2KTk5uUC9Xbp0kfT/fzdLMi64sjBHB1csHx8ftWjRQi1atFCDBg00YMAAzZ49W4mJiY4+Rf0/YleusMqfo/LAAw8oLi6u0D5NmjQp8vrcpX379tqxY4f++9//auHChfrXv/6l1157TdOmTSvyJc9WX3l2sT+PvLw8S7dzOd7e3oW2m/MmTpd1Nput0HovHMuqVasqJSVFCxYs0Pz58zV//nxNnz5dDz30kGbOnHnJbeR/12fNmqXq1asXWH7hJf6+vr5FvsXBPffco3bt2mnu3LlauHChXn75ZU2aNElffvmlunXrdsnPPvHEE1q5cqUWL15cIFjZ7XZ17dpVTz75ZKGfzQ/fJRkXXFkIOvAIzZs3l/Tn6QXpz0nLdrtd27Ztc0yqlP6cRJyZmVnkm7cV9os5NDRUFStWVF5enuP/Ei8mIiJCqampMsY4rWvLli2X3XZoaKj8/f0LPfVUlM9LUuXKlTVgwAANGDBAJ0+eVPv27TV27FhH0LHi1Eg+u92unTt3On6RSNLWrVslyTFhNf/ISWZmptNnCztlVNTaQkNDFRAQUOiY/Pbbb/Ly8lJ4eHiR1nUp+d+Zi22nSpUqxbqUul69ekpNTXX5c5UqVSr0CEhhY+nj46MePXqoR48estvtGjJkiN599109++yzhR5xPL826c9QcLnvenHUqFFDQ4YM0ZAhQ3T48GE1bdpUL7744iWDzqeffqqpU6dq6tSpjqvbLqz55MmTRar3cuMCz8CpK1xRli5dWuj/xebPj8g/rXD77bdLkuMKi3xTpkyRJHXv3r1I2wsMDCzwS9nb21t33XWX5syZU+gvqPMvU7799tt14MABffHFF462U6dO6b333rvstr29vRUbG6t58+Zp7969jvbNmzdrwYIFl/38sWPHnN5XqFBB9evXd7q0Pf8X84X7WFxvvfWW42djjN566y2VL19enTt3lvRnWPD29nbMk8h3/hwLV2vz9vbWrbfeqv/+979Op8gOHTqkjz/+WG3btlVQUFAx9+j/1ahRQzExMZo5c6ZTTampqVq4cKHjO+equ+66Sxs2bChwJaB06SNM9erV02+//eb0fduwYYNWrFjh1O/C74GXl5fjiGP+d+FiYx0bG6ugoCC99NJLOnfuXIEainpJ/oXy8vKUlZXl1Fa1alWFhYVd8tYLqampevjhh/XAAw9c9Mqoe+65RytXriz070hmZqZyc3MlFW1c4Bk4ooMryrBhw3Tq1Cn17t1bDRs21NmzZ/XTTz/ps88+U2RkpGOOTHR0tOLi4vTee+8pMzNTHTp00OrVqzVz5kz16tVLnTp1KtL2mjVrpsWLF2vKlCkKCwtTnTp11KpVK02cOFFLly5Vq1atNHjwYEVFRSkjI0O//PKLFi9erIyMDEnS4MGD9dZbb+mhhx7SunXrVKNGDc2aNeuSl92e7/nnn1dycrLatWunIUOGKDc3V2+++aYaNWp02fk2UVFR6tixo5o1a6bKlStr7dq1jst5z98/SRo+fLhiY2Pl7e2t++67r0i1XcjPz0/JycmKi4tTq1atNH/+fH377bd6+umnHROag4ODdffdd+vNN9+UzWZTvXr19M033xR6s0dXahs/frwWLVqktm3basiQISpXrpzeffddnTlzRpMnTy7W/hTm5ZdfVrdu3dS6dWsNGjTIcXl5cHBwsR9X8cQTT+iLL77Q3XffrYEDB6pZs2bKyMjQV199pWnTpik6OrrQzw0cOFBTpkxRbGysBg0apMOHD2vatGlq1KiR0+Trhx9+WBkZGbrllltUq1Yt7dmzR2+++aZiYmIcRztjYmLk7e2tSZMmKSsrS76+vrrllltUtWpVvfPOO3rwwQfVtGlT3XfffQoNDdXevXv17bffqk2bNk7htqhOnDihWrVqqU+fPoqOjlaFChW0ePFirVmzRq+++upFP5f/97t9+/b68MMPnZbdfPPNqlu3rp544gl99dVX+tvf/qb+/furWbNmysnJ0aZNm/TFF19o9+7dqlKlSpHGBR6i9C74Alw3f/58M3DgQNOwYUNToUIF4+PjY+rXr2+GDRtmDh065NT33Llz5vnnnzd16tQx5cuXN+Hh4WbMmDHmjz/+cOoXERFx0cuwf/vtN9O+fXvj7+9vJDldvn3o0CETHx9vwsPDTfny5U316tVN586dzXvvvee0jj179pg77rjDBAQEmCpVqpgRI0aY5OTkIl1ebowxP/74o2nWrJnx8fExdevWNdOmTTOJiYmXvbx8/PjxpmXLliYkJMT4+/ubhg0bmhdffNHpUvzc3FwzbNgwExoaamw2m2Od+Zczv/zyywXqudjl5YGBgWbHjh3m1ltvNQEBAaZatWomMTHR5OXlOX3+yJEj5q677jIBAQGmUqVK5tFHHzWpqakF1nmx2owpeHm5Mcb88ssvJjY21lSoUMEEBASYTp06mZ9++smpT/7l5WvWrHFqv9hl74VZvHixadOmjfH39zdBQUGmR48eJi0trdD1FeXycmOMOXbsmBk6dKipWbOm8fHxMbVq1TJxcXGO2xcUNubGGPPhhx+aunXrGh8fHxMTE2MWLFhQ4JLsL774wtx6662matWqxsfHx9SuXds8+uij5uDBg07rev/9903dunUdty44fyyWLl1qYmNjTXBwsPHz8zP16tUz/fv3N2vXrnX0yf8OFMWZM2fME088YaKjo03FihVNYGCgiY6ONm+//bZTvwv3Jf9WEIW9zh+bEydOmDFjxpj69esbHx8fU6VKFXPzzTebV155xfH9L+q44MpnM+YKmn0HAADgAuboAAAAj0XQAQAAHougAwAAPBZBBwAAeCyCDgAA8FgEHQAA4LGu+hsG2u12HThwQBUrVrT0dvgAAMB9jDE6ceKEwsLCLvmMtas+6Bw4cMCS5+AAAIC/3r59+wo83PV8V33QqVixoqQ/B8qK5+EAAAD3y87OVnh4uOP3+MVc9UEn/3RVUFAQQQcAgCvM5aadMBkZAAB4LIIOAADwWAQdAADgsQg6AADAYxF0AACAxyLoAAAAj0XQAQAAHougAwAAPBZBBwAAeCyCDgAA8FgEHQAA4LEIOgAAwGMRdAAAgMci6AAAAI9F0AEAAB6LoAMAADwWQQcAAHgsgg4AAPBYBB0AAOCxCDoAAMBjEXQAAIDHIugAAACPRdABAAAei6ADAAA8FkEHAAB4LIIOAADwWAQdAADgsQg6AADAYxF0AACAxyLoAAAAj3XFB519+/apY8eOioqKUpMmTTR79uzSLgkAAJQR5Uq7gJIqV66cpk6dqpiYGKWnp6tZs2a6/fbbFRgYWNqlAQCAUnbFB50aNWqoRo0akqTq1aurSpUqysjIIOgAAIDSP3W1bNky9ejRQ2FhYbLZbJo3b16BPklJSYqMjJSfn59atWql1atXF7qudevWKS8vT+Hh4W6uGgAAXAlKPejk5OQoOjpaSUlJhS7/7LPPlJCQoMTERP3yyy+Kjo5WbGysDh8+7NQvIyNDDz30kN57772/omwAAHAFsBljTGkXkc9ms2nu3Lnq1auXo61Vq1Zq0aKF3nrrLUmS3W5XeHi4hg0bptGjR0uSzpw5o65du2rw4MF68MEHL7mNM2fO6MyZM4732dnZCg8PV1ZWloKCgqzfKQAAYLns7GwFBwdf9vd3qR/RuZSzZ89q3bp16tKli6PNy8tLXbp00cqVKyVJxhj1799ft9xyy2VDjiRNmDBBwcHBjhenuQAA8FxlOugcPXpUeXl5qlatmlN7tWrVlJ6eLklasWKFPvvsM82bN08xMTGKiYnRpk2bLrrOMWPGKCsry/Hat2+fW/cBAICrVeTob0u7hCv/qqu2bdvKbrcXub+vr698fX3dWBEAACgryvQRnSpVqsjb21uHDh1yaj906JCqV69eSlUBAIArRZkOOj4+PmrWrJmWLFniaLPb7VqyZIlat25dipUBAIArQamfujp58qS2b9/ueL9r1y6lpKSocuXKql27thISEhQXF6fmzZurZcuWmjp1qnJycjRgwIBSrBoAAFwJSj3orF27Vp06dXK8T0hIkCTFxcVpxowZuvfee3XkyBE999xzSk9PV0xMjJKTkwtMUAYAALhQmbqPzl8pKSlJSUlJysvL09atW7mPDgAAFosc/a12T+zulnV7xH103Ck+Pl5paWlas2ZNaZcCAADc5KoNOgAAwPMRdAAAgMci6AAAAI9F0AEAAB6LoAMAADwWQQcAAHisqzboJCUlKSoqSi1atCjtUgAAgJtctUGH++gAAOD5rtqgAwAAPB9BBwAAeCyCDgAA8FgEHQAA4LEIOgAAwGMRdAAAgMci6AAAAI911QYdbhgIAIDnu2qDDjcMBADA8121QQcAAHg+gg4AAPBYBB0AAOCxCDoAAMBjEXQAAIDHIugAAACPRdABAAAei6ADAAA8FkEHAAB4rKs26PAICAAAPN9VG3R4BAQAAJ7vqg06AADA8xF0AACAx3Ip6Jw7d06dO3fWtm3b3FUPAACAZVwKOuXLl9fGjRvdVQsAAIClXD519cADD+jf//63O2oBAACwVDlXP5Cbm6sPPvhAixcvVrNmzRQYGOi0fMqUKZYVBwAAUBIuB53U1FQ1bdpUkrR161anZTabzZqqAAAALOBy0Fm6dKk76gAAALBcsS8v3759uxYsWKDTp09LkowxlhUFAABgBZeDzrFjx9S5c2c1aNBAt99+uw4ePChJGjRokB5//HHLCwQAACgul4POqFGjVL58ee3du1cBAQGO9nvvvVfJycmWFudOPOsKAADP5/IcnYULF2rBggWqVauWU/u1116rPXv2WFaYu8XHxys+Pl7Z2dkKDg4u7XIAAIAbuHxEJycnx+lITr6MjAz5+vpaUhQAAIAVXA467dq103/+8x/He5vNJrvdrsmTJ6tTp06WFgcAAFASLp+6mjx5sjp37qy1a9fq7NmzevLJJ/Xrr78qIyNDK1ascEeNAAAAxeLyEZ0bbrhBW7duVdu2bdWzZ0/l5OTozjvv1Pr161WvXj131AgAAFAsLh/R2bt3r8LDw/XMM88Uuqx27dqWFAYAAFBSLh/RqVOnjo4cOVKg/dixY6pTp44lRQEAAFjB5aBjjCn0mVYnT56Un5+fJUUBAABYocinrhISEiT9eZXVs88+63SJeV5enlatWqWYmBjLCwQAACiuIged9evXS/rziM6mTZvk4+PjWObj46Po6Gj94x//sL5CAACAYipy0Ml/avmAAQP0+uuvKygoyG1FAQAAWMHlOTrTp093CjnZ2dmaN2+efvvtN0sLAwAAKCmXg84999yjt956S5J0+vRpNW/eXPfcc48aN26sOXPmWF4gAABAcbkcdJYtW6Z27dpJkubOnStjjDIzM/XGG29o/PjxlhcIAABQXC4HnaysLFWuXFmSlJycrLvuuksBAQHq3r27tm3bZnmB7pKUlKSoqCi1aNGitEsBAABu4nLQCQ8P18qVK5WTk6Pk5GTdeuutkqTjx49fUffRiY+PV1pamtasWVPapQAAADdx+REQI0eOVL9+/VShQgVFRESoY8eOkv48pdW4cWOr6wMAACg2l4POkCFD1LJlS+3bt09du3aVl9efB4Xq1q3LHB0AAFCmuBx0JKl58+Zq3ry5U1v37t0tKQgAAMAqLgedgQMHXnL5Bx98UOxiAAAArORy0Dl+/LjT+3Pnzik1NVWZmZm65ZZbLCsMAACgpFwOOnPnzi3QZrfb9fe//1316tWzpCgAAAAruHx5eaEr8fJSQkKCXnvtNStWBwAAYAlLgo4k7dixQ7m5uVatDgAAXMEiR39b2iVIKsapq4SEBKf3xhgdPHhQ3377reLi4iwrDAAAoKRcDjrr1693eu/l5aXQ0FC9+uqrl70iCwAA4K/kctBZunSpO+oAAACwnGVzdAAAAMqaIh3RufHGG2Wz2Yq0wl9++aVEBQEAAFilSEGnV69ebi4DAADAekUKOomJie6uAwAAwHIuz9FZs2aNVq1aVaB91apVWrt2rSVFAQAAWMHloBMfH699+/YVaN+/f7/i4+MtKQoAAMAKLgedtLQ0NW3atED7jTfeqLS0NEuKAgAAsILLQcfX11eHDh0q0H7w4EGVK+fybXlKTVJSkqKiotSiRYvSLgUAALiJy0Hn1ltv1ZgxY5SVleVoy8zM1NNPP62uXbtaWpw7xcfHKy0tTWvWrCntUgAAgJu4fAjmlVdeUfv27RUREaEbb7xRkpSSkqJq1app1qxZlhcIAABQXC4HnZo1a2rjxo366KOPtGHDBvn7+2vAgAHq27evypcv744aAQAAiqVYk2oCAwP1yCOPWF0LAACApXjWFQAA8FgEHQAA4LEIOgAAwGMRdAAAgMcq9h3+zp49q8OHD8tutzu1165du8RFAQAAWMHloLNt2zYNHDhQP/30k1O7MUY2m015eXmWFQcAAFASLged/v37q1y5cvrmm29Uo0YN2Ww2d9QFAABQYi4HnZSUFK1bt04NGzZ0Rz0AAACWcXkyclRUlI4ePeqOWgAAACzlctCZNGmSnnzySf3www86duyYsrOznV4AAABlhcunrrp06SJJ6ty5s1M7k5EBAIAkRY7+trRLcHA56CxdutQddQAAAFjO5aDToUMHd9QBAABguSIFnY0bN+qGG26Ql5eXNm7ceMm+TZo0saQwAACAkipS0ImJiVF6erqqVq2qmJgY2Ww2GWMK9GOODgAAKEuKFHR27dql0NBQx88AAABXgiIFnYiIiEJ/BgAAKMt4ejkAAPBYBB0AAOCxCDoAAMBjEXQAAIDHcjno7Nu3T7///rvj/erVqzVy5Ei99957lhYGAABQUi4Hnfvvv9/xGIj09HR17dpVq1ev1jPPPKNx48ZZXiAAAEBxuRx0UlNT1bJlS0nS559/rhtuuEE//fSTPvroI82YMcPq+twmKSlJUVFRatGiRWmXAgAA3MTloHPu3Dn5+vpKkhYvXqw77rhDktSwYUMdPHjQ2urcKD4+XmlpaVqzZk1plwIAANzE5aDTqFEjTZs2Tf/73/+0aNEi3XbbbZKkAwcO6JprrrG8QAAAgOJyOehMmjRJ7777rjp27Ki+ffsqOjpakvTVV185TmkBAACUBUV6BMT5OnbsqKNHjyo7O1uVKlVytD/yyCMKCAiwtDgAAICSKNZ9dIwxWrdund59912dOHFCkuTj40PQAQAAZYrLR3T27Nmj2267TXv37tWZM2fUtWtXVaxYUZMmTdKZM2c0bdo0d9QJAADgMpeP6IwYMULNmzfX8ePH5e/v72jv3bu3lixZYmlxAAAAJeHyEZ3//e9/+umnn+Tj4+PUHhkZqf3791tWGAAAQEm5fETHbrcrLy+vQPvvv/+uihUrWlIUAACAFVwOOrfeequmTp3qeG+z2XTy5EklJibq9ttvt7I2AACAEnH51NWrr76q2NhYRUVF6Y8//tD999+vbdu2qUqVKvrkk0/cUSMAALhCRI7+trRLcOJy0KlVq5Y2bNigTz/9VBs3btTJkyc1aNAg9evXz2lyMgAAQGlzOehIUrly5fTAAw9YXQsAAIClihR0vvrqqyKvMP8hnwAAAKWtSEGnV69eRVqZzWYr9IosAADg+cra/BypiEHHbre7uw4AAADLFetZVwAAAFcClycjjxs37pLLn3vuuWIXAwAAYCWXg87cuXOd3p87d067du1SuXLlVK9ePYIOAAAoM1wOOuvXry/Qlp2drf79+6t3796WFAUAAGAFS+boBAUF6fnnn9ezzz5rxeoAAAAsYdlk5KysLGVlZVm1OgAAgBJz+dTVG2+84fTeGKODBw9q1qxZ6tatm2WFAQAAlJTLQee1115zeu/l5aXQ0FDFxcVpzJgxlhUGAABQUi4HnV27drmjDgAAAMu5PEcnKytLGRkZBdozMjKUnZ1tSVEAAABWcDno3Hffffr0008LtH/++ee67777LCkKAADACi4HnVWrVqlTp04F2jt27KhVq1ZZUhQAAIAVXA46Z86cUW5uboH2c+fO6fTp05YUBQAAYAWXg07Lli313nvvFWifNm2amjVrZklRAAAAVnD5qqvx48erS5cu2rBhgzp37ixJWrJkidasWaOFCxdaXiAAAEBxuXxEp02bNlq5cqVq1aqlzz//XF9//bXq16+vjRs3ql27du6oEQAAoFhcPqIjSTExMfr444+trgUAAMBSxXrW1Y4dO/TPf/5T999/vw4fPixJmj9/vn799VdLiwMAAFeGyNHflnYJhXI56Pz4449q3LixVq1apTlz5ujkyZOSpA0bNigxMdHyAgEAAIrL5aAzevRojR8/XosWLZKPj4+j/ZZbbtHPP/9saXEAAAAl4XLQ2bRpk3r37l2gvWrVqjp69KglRQEAAFjB5aATEhKigwcPFmhfv369atasaUlRAAAAVijWs66eeuoppaeny2azyW63a8WKFfrHP/6hhx56yB01AgAAFIvLQeell15Sw4YNFR4erpMnTyoqKkrt27fXzTffrH/+85/uqBEAAKBYXA46Pj4+ev/997Vjxw598803+vDDD/Xbb79p1qxZ8vb2dkeNl9W7d29VqlRJffr0KZXtAwCAsqlYNwyUpNq1a6t27dpW1lJsI0aM0MCBAzVz5szSLgUAAJQhRQo6CQkJRV7hlClTil1McXXs2FE//PDDX75dAABQthUp6Kxfv75IK7PZbC4XsGzZMr388stat26dDh48qLlz56pXr15OfZKSkvTyyy8rPT1d0dHRevPNN9WyZUuXtwUAAK4uRQo6S5cudVsBOTk5io6O1sCBA3XnnXcWWP7ZZ58pISFB06ZNU6tWrTR16lTFxsZqy5Ytqlq1qtvqAgAAV75iz9GRpN9//12SVKtWrWKvo1u3burWrdtFl0+ZMkWDBw/WgAEDJEnTpk3Tt99+qw8++ECjR492eXtnzpzRmTNnHO+zs7NdLxoAAFwRXL7qym63a9y4cQoODlZERIQiIiIUEhKiF154QXa73dLizp49q3Xr1qlLly7/X7CXl7p06aKVK1cWa50TJkxQcHCw4xUeHm5VuQAAoIxx+YjOM888o3//+9+aOHGi2rRpI0lavny5xo4dqz/++EMvvviiZcUdPXpUeXl5qlatmlN7tWrV9Ntvvzned+nSRRs2bFBOTo5q1aql2bNnq3Xr1oWuc8yYMU6Tq7Ozswk7AAB4KJeDzsyZM/Wvf/1Ld9xxh6OtSZMmqlmzpoYMGWJp0CmqxYsXF7mvr6+vfH193VgNAAAoK1w+dZWRkaGGDRsWaG/YsKEyMjIsKSpflSpV5O3trUOHDjm1Hzp0SNWrV7d0WwAAwPO4HHSio6P11ltvFWh/6623FB0dbUlR+Xx8fNSsWTMtWbLE0Wa327VkyZKLnpoCAADI5/Kpq8mTJ6t79+5avHixI2ysXLlS+/bt03fffedyASdPntT27dsd73ft2qWUlBRVrlxZtWvXVkJCguLi4tS8eXO1bNlSU6dOVU5OjuMqLAAAgItxOeh06NBBW7duVVJSkmNC8J133qkhQ4YoLCzM5QLWrl2rTp06Od7nTxSOi4vTjBkzdO+99+rIkSN67rnnlJ6erpiYGCUnJxeYoAwAAHAhmzHGlHYRpSEpKUlJSUnKy8vT1q1blZWVpaCgoNIuCwCAK0bk6G8lSbsndnf8fKHdE7u7ZdvZ2dkKDg6+7O/vYt0w8I8//tDGjRt1+PDhAvfOOf9qrLIsPj5e8fHxjoECAACex+Wgk5ycrIceekhHjx4tsMxmsykvL8+SwgAAAErK5auuhg0bprvvvlsHDx6U3W53ehFyAABAWeJy0Dl06JASEhKYDAwAAMo8l4NOnz599MMPP7ihFAAAAGu5PEfnrbfe0t13363//e9/aty4scqXL++0fPjw4ZYVBwAAUBIuB51PPvlECxculJ+fn3744QfZbDbHMpvNRtABAABlRrGeXv78889r9OjR8vJy+cxXmXH+fXQAAEDRnH/vnCuBy0nl7Nmzuvfee6/okCP9eR+dtLQ0rVmzprRLAQAAbuJyWomLi9Nnn33mjloAAAAs5fKpq7y8PE2ePFkLFixQkyZNCkxGnjJlimXFAQAAlITLQWfTpk268cYbJUmpqalOy86fmAwAAFDaXA46S5cudUcdAAAAlruyZxQDAABcAkEHAAB4LIIOAADwWFdt0ElKSlJUVJRatGhR2qUAAFCm5d8kMP+/V5IiBZ2mTZvq+PHjkqRx48bp1KlTbi3qr8ANAwEA8HxFCjqbN29WTk6OJOn555/XyZMn3VoUAACAFYp0eXlMTIwGDBigtm3byhijV155RRUqVCi073PPPWdpgQAAAMVVpKAzY8YMJSYm6ptvvpHNZtP8+fNVrlzBj9psNoIOAAAoM4oUdK677jp9+umnkiQvLy8tWbJEVatWdWthAAAAJeXynZHtdrs76gAAALCcy0FHknbs2KGpU6dq8+bNkqSoqCiNGDFC9erVs7Q4AACAknD5PjoLFixQVFSUVq9erSZNmqhJkyZatWqVGjVqpEWLFrmjRgAAgGJx+YjO6NGjNWrUKE2cOLFA+1NPPaWuXbtaVhwAAEBJuHxEZ/PmzRo0aFCB9oEDByotLc2SogAAAKzgctAJDQ1VSkpKgfaUlJQr6kosHgEBAMDlXYmPfTify6euBg8erEceeUQ7d+7UzTffLElasWKFJk2apISEBMsLdJf4+HjFx8crOztbwcHBpV0OAABwA5eDzrPPPquKFSvq1Vdf1ZgxYyRJYWFhGjt2rIYPH255gQAAAMXlctCx2WwaNWqURo0apRMnTkiSKlasaHlhAAAAJVWs++jkI+AAAICyzOXJyAAAAFcKgg4AAPBYBB0AAOCxXAo6586dU+fOnbVt2zZ31QMAAGAZl4JO+fLltXHjRnfVAgAAYCmXT1098MAD+ve//+2OWgAAACzl8uXlubm5+uCDD7R48WI1a9ZMgYGBTsunTJliWXEAAAAl4XLQSU1NVdOmTSVJW7dudVpms9msqeovkJSUpKSkJOXl5ZV2KQAAlLrCnmm1e2L3UqjEWi4HnaVLl7qjjr8cz7oCAMDzFfvy8u3bt2vBggU6ffq0JMkYY1lRAAAAVnA56Bw7dkydO3dWgwYNdPvtt+vgwYOSpEGDBunxxx+3vEAAAIDicjnojBo1SuXLl9fevXsVEBDgaL/33nuVnJxsaXEAAAAl4fIcnYULF2rBggWqVauWU/u1116rPXv2WFYYAABASbl8RCcnJ8fpSE6+jIwM+fr6WlIUAACAFVwOOu3atdN//vMfx3ubzSa73a7JkyerU6dOlhYHAABQEi6fupo8ebI6d+6stWvX6uzZs3ryySf166+/KiMjQytWrHBHjQAAAMXi8hGdG264QVu3blXbtm3Vs2dP5eTk6M4779T69etVr149d9QIAABQLC4f0ZGk4OBgPfPMM1bXAgAAYKliBZ3jx4/r3//+tzZv3ixJioqK0oABA1S5cmVLiwMAACgJl09dLVu2TJGRkXrjjTd0/PhxHT9+XG+88Ybq1KmjZcuWuaNGAACAYnH5iE58fLzuvfdevfPOO/L29pYk5eXlaciQIYqPj9emTZssLxIAAKA4XD6is337dj3++OOOkCNJ3t7eSkhI0Pbt2y0tDgAAoCRcPqLTtGlTbd68Wdddd51T++bNmxUdHW1ZYe6WlJSkpKQk5eXllXYpAAC4ReTobyVJuyd2d/x8oUst8wRFCjobN250/Dx8+HCNGDFC27dv10033SRJ+vnnn5WUlKSJEye6p0o3iI+PV3x8vLKzsxUcHFza5QAAADcoUtCJiYmRzWaTMcbR9uSTTxbod//99+vee++1rjoAAIASKFLQ2bVrl7vrAAAAsFyRgk5ERIS76wAAALBcsW4YeODAAS1fvlyHDx+W3W53WjZ8+HBLCgMAACgpl4POjBkz9Oijj8rHx0fXXHONbDabY5nNZiPoAACAMsPloPPss8/queee05gxY+Tl5fJteAAAAP4yLieVU6dO6b777iPkAACAMs/ltDJo0CDNnj3bHbUAAABYyuVTVxMmTNDf/vY3JScnq3HjxipfvrzT8ilTplhWHAAAQEkUK+gsWLDA8QiICycjAwAAlBUuB51XX31VH3zwgfr37++GcgAAAKzj8hwdX19ftWnTxh21AAAAWMrloDNixAi9+eab7qgFAADAUi6fulq9erW+//57ffPNN2rUqFGBychffvmlZcUBAACUhMtBJyQkRHfeeac7agEAALCUy0Fn+vTp7qgDAADActzeGAAAeCyXj+jUqVPnkvfL2blzZ4kKAgAAsIrLQWfkyJFO78+dO6f169crOTlZTzzxhFV1uV1SUpKSkpKUl5dX2qUAAKDI0d8W2r57YnfHsvN/vtwy/MnloDNixIhC25OSkrR27doSF/RXiY+PV3x8vLKzsxUcHFza5QAAADewbI5Ot27dNGfOHKtWBwAAUGKWBZ0vvvhClStXtmp1AAAAJebyqasbb7zRaTKyMUbp6ek6cuSI3n77bUuLAwAAKAmXg06vXr2c3nt5eSk0NFQdO3ZUw4YNraoLAACgxFwOOomJie6oAwAAwHLcMBAAAHisIh/R8fLyuuSNAiXJZrMpNze3xEUBAABYochBZ+7cuRddtnLlSr3xxhuy2+2WFAUAAGCFIgednj17FmjbsmWLRo8era+//lr9+vXTuHHjLC0OAACgJIo1R+fAgQMaPHiwGjdurNzcXKWkpGjmzJmKiIiwuj4AAIBicynoZGVl6amnnlL9+vX166+/asmSJfr66691ww03uKs+AACAYivyqavJkydr0qRJql69uj755JNCT2UBAACUJUUOOqNHj5a/v7/q16+vmTNnaubMmYX2+/LLLy0rDgAAoCSKHHQeeuihy15eDgAAUJYUOejMmDHDjWUAAABYjzsjAwAAj0XQAQAAHougAwAAPBZBBwAAeCyCDgAA8FgEHQAA4LEIOgAAwGMRdAAAgMci6AAAAI9F0AEAAB6LoAMAADxWkZ915WmSkpKUlJSkvLy80i4FAK56kaO/lSTtntjd8fOFirrswn6XWlaW1g/3uGqP6MTHxystLU1r1qwp7VIAAICbXLVBBwAAeD6CDgAA8FgEHQAA4LEIOgAAwGMRdAAAgMci6AAAAI9F0AEAAB6LoAMAADwWQQcAAHgsgg4AAPBYBB0AAOCxCDoAAMBjEXQAAIDHIugAAACPRdABAAAei6ADAAA8FkEHAAB4LIIOAADwWAQdAADgsQg6AADAYxF0AACAxyLoAAAAj0XQAQAAHougAwAAPBZBBwAAeCyCDgAA8FgEHQAA4LEIOgAAwGMRdAAAgMci6AAAAI9F0AEAAB6LoAMAADwWQQcAAHgsgg4AAPBYBB0AAOCxCDoAAMBjEXQAAIDHIugAAACPRdABAAAeyyOCzjfffKPrrrtO1157rf71r3+VdjkAAKCMKFfaBZRUbm6uEhIStHTpUgUHB6tZs2bq3bu3rrnmmtIuDQAAlLIr/ojO6tWr1ahRI9WsWVMVKlRQt27dtHDhwtIuCwAAlAGlHnSWLVumHj16KCwsTDabTfPmzSvQJykpSZGRkfLz81OrVq20evVqx7IDBw6oZs2ajvc1a9bU/v37/4rSAQBAGVfqQScnJ0fR0dFKSkoqdPlnn32mhIQEJSYm6pdfflF0dLRiY2N1+PDhv7hSAABwpSn1oNOtWzeNHz9evXv3LnT5lClTNHjwYA0YMEBRUVGaNm2aAgIC9MEHH0iSwsLCnI7g7N+/X2FhYRfd3pkzZ5Sdne30AgAAnqlMT0Y+e/as1q1bpzFjxjjavLy81KVLF61cuVKS1LJlS6Wmpmr//v0KDg7W/Pnz9eyzz150nRMmTNDzzz/v9trzRY7+ttD23RO7F2nZhf0utaw462fbV9e2XVl/aW67OOtn256xbcBqpX5E51KOHj2qvLw8VatWzam9WrVqSk9PlySVK1dOr776qjp16qSYmBg9/vjjl7ziasyYMcrKynK89u3b59Z9AAAApadMH9EpqjvuuEN33HFHkfr6+vrK19fXzRUBAICyoEwf0alSpYq8vb116NAhp/ZDhw6pevXqpVQVAAC4UpTpoOPj46NmzZppyZIljja73a4lS5aodevWpVgZAAC4EpT6qauTJ09q+/btjve7du1SSkqKKleurNq1ayshIUFxcXFq3ry5WrZsqalTpyonJ0cDBgwoxaoBAMCVoNSDztq1a9WpUyfH+4SEBElSXFycZsyYoXvvvVdHjhzRc889p/T0dMXExCg5ObnABGUAAIALlXrQ6dixo4wxl+wzdOhQDR061NLtJiUlKSkpSXl5eZauFwAAlB1leo6OO8XHxystLU1r1qwp7VIAAICbXLVBBwAAeD6CDgAA8FgEHQAA4LEIOgAAwGMRdAAAgMci6AAAAI911QadpKQkRUVFqUWLFqVdCgAAcJNSv2FgaYmPj1d8fLyysrIUEhKi7Oxst2zHfuZUoe3Z2dlFWnZhv0stK8762fbVtW1X1l+a2y7O+tk22+a7XHa37Q75673cTYdt5nI9PNzvv/+u8PDw0i4DAAAUw759+1SrVq2LLr/qg47dbteBAwdUsWJF2Ww2y9abnZ2t8PBw7du3T0FBQZatF39ifN2HsXUfxta9GF/3KYtja4zRiRMnFBYWJi+vi8/EuWpPXeXz8vK6ZBIsqaCgoDLzpfBEjK/7MLbuw9i6F+PrPmVtbIODgy/b56qdjAwAADwfQQcAAHgsgo6b+Pr6KjExUb6+vqVdikdifN2HsXUfxta9GF/3uZLH9qqfjAwAADwXR3QAAIDHIugAAACPRdABAAAei6ADAAA8FkHHTZKSkhQZGSk/Pz+1atVKq1evLu2SyrQJEyaoRYsWqlixoqpWrapevXppy5YtTn3++OMPxcfH65prrlGFChV011136dChQ0599u7dq+7duysgIEBVq1bVE088odzc3L9yV8q8iRMnymazaeTIkY42xrZk9u/frwceeEDXXHON/P391bhxY61du9ax3Bij5557TjVq1JC/v7+6dOmibdu2Oa0jIyND/fr1U1BQkEJCQjRo0CCdPHnyr96VMiUvL0/PPvus6tSpI39/f9WrV08vvPCC07ONGNuiW7ZsmXr06KGwsDDZbDbNmzfPablVY7lx40a1a9dOfn5+Cg8P1+TJk929a5dmYLlPP/3U+Pj4mA8++MD8+uuvZvDgwSYkJMQcOnSotEsrs2JjY8306dNNamqqSUlJMbfffrupXbu2OXnypKPPY489ZsLDw82SJUvM2rVrzU033WRuvvlmx/Lc3Fxzww03mC5dupj169eb7777zlSpUsWMGTOmNHapTFq9erWJjIw0TZo0MSNGjHC0M7bFl5GRYSIiIkz//v3NqlWrzM6dO82CBQvM9u3bHX0mTpxogoODzbx588yGDRvMHXfcYerUqWNOnz7t6HPbbbeZ6Oho8/PPP5v//e9/pn79+qZv376lsUtlxosvvmiuueYa880335hdu3aZ2bNnmwoVKpjXX3/d0YexLbrvvvvOPPPMM+bLL780kszcuXOdllsxlllZWaZatWqmX79+JjU11XzyySfG39/fvPvuu3/VbhZA0HGDli1bmvj4eMf7vLw8ExYWZiZMmFCKVV1ZDh8+bCSZH3/80RhjTGZmpilfvryZPXu2o8/mzZuNJLNy5UpjzJ9/ib28vEx6erqjzzvvvGOCgoLMmTNn/todKINOnDhhrr32WrNo0SLToUMHR9BhbEvmqaeeMm3btr3ocrvdbqpXr25efvllR1tmZqbx9fU1n3zyiTHGmLS0NCPJrFmzxtFn/vz5xmazmf3797uv+DKue/fuZuDAgU5td955p+nXr58xhrEtiQuDjlVj+fbbb5tKlSo5/bvw1FNPmeuuu87Ne3RxnLqy2NmzZ7Vu3Tp16dLF0ebl5aUuXbpo5cqVpVjZlSUrK0uSVLlyZUnSunXrdO7cOadxbdiwoWrXru0Y15UrV6px48aqVq2ao09sbKyys7P166+//oXVl03x8fHq3r270xhKjG1JffXVV2revLnuvvtuVa1aVTfeeKPef/99x/Jdu3YpPT3daXyDg4PVqlUrp/ENCQlR8+bNHX26dOkiLy8vrVq16q/bmTLm5ptv1pIlS7R161ZJ0oYNG7R8+XJ169ZNEmNrJavGcuXKlWrfvr18fHwcfWJjY7VlyxYdP378L9obZ1f9Qz2tdvToUeXl5Tn9QpCkatWq6bfffiulqq4sdrtdI0eOVJs2bXTDDTdIktLT0+Xj46OQkBCnvtWqVVN6erqjT2Hjnr/savbpp5/ql19+0Zo1awosY2xLZufOnXrnnXeUkJCgp59+WmvWrNHw4cPl4+OjuLg4x/gUNn7nj2/VqlWdlpcrV06VK1e+qsd39OjRys7OVsOGDeXt7a28vDy9+OKL6tevnyQxthayaizT09NVp06dAuvIX1apUiW31H8pBB2UOfHx8UpNTdXy5ctLuxSPsG/fPo0YMUKLFi2Sn59faZfjcex2u5o3b66XXnpJknTjjTcqNTVV06ZNU1xcXClXd2X7/PPP9dFHH+njjz9Wo0aNlJKSopEjRyosLIyxRZFx6spiVapUkbe3d4ErVg4dOqTq1auXUlVXjqFDh+qbb77R0qVLVatWLUd79erVdfbsWWVmZjr1P39cq1evXui45y+7Wq1bt06HDx9W06ZNVa5cOZUrV04//vij3njjDZUrV07VqlVjbEugRo0aioqKcmq7/vrrtXfvXkn/Pz6X+jehevXqOnz4sNPy3NxcZWRkXNXj+8QTT2j06NG677771LhxYz344IMaNWqUJkyYIImxtZJVY1kW/60g6FjMx8dHzZo105IlSxxtdrtdS5YsUevWrUuxsrLNGKOhQ4dq7ty5+v777wsc+mzWrJnKly/vNK5btmzR3r17HePaunVrbdq0yekv4qJFixQUFFTgF9HVpHPnztq0aZNSUlIcr+bNm6tfv36Onxnb4mvTpk2BWyFs3bpVERERkqQ6deqoevXqTuObnZ2tVatWOY1vZmam1q1b5+jz/fffy263q1WrVn/BXpRNp06dkpeX868pb29v2e12SYytlaway9atW2vZsmU6d+6co8+iRYt03XXXlcppK0lcXu4On376qfH19TUzZswwaWlp5pFHHjEhISFOV6zA2d///ncTHBxsfvjhB3Pw4EHH69SpU44+jz32mKldu7b5/vvvzdq1a03r1q1N69atHcvzL4G+9dZbTUpKiklOTjahoaFcAl2I86+6MoaxLYnVq1ebcuXKmRdffNFs27bNfPTRRyYgIMB8+OGHjj4TJ040ISEh5r///a/ZuHGj6dmzZ6GX7d54441m1apVZvny5ebaa6+9Ki+BPl9cXJypWbOm4/LyL7/80lSpUsU8+eSTjj6MbdGdOHHCrF+/3qxfv95IMlOmTDHr1683e/bsMcZYM5aZmZmmWrVq5sEHHzSpqanm008/NQEBAVxe7onefPNNU7t2bePj42Natmxpfv7559IuqUyTVOhr+vTpjj6nT582Q4YMMZUqVTIBAQGmd+/e5uDBg07r2b17t+nWrZvx9/c3VapUMY8//rg5d+7cX7w3Zd+FQYexLZmvv/7a3HDDDcbX19c0bNjQvPfee07L7Xa7efbZZ021atWMr6+v6dy5s9myZYtTn2PHjpm+ffuaChUqmKCgIDNgwABz4sSJv3I3ypzs7GwzYsQIU7t2bePn52fq1q1rnnnmGadLlxnbolu6dGmh/87GxcUZY6wbyw0bNpi2bdsaX19fU7NmTTNx4sS/ahcLZTPmvFtMAgAAeBDm6AAAAI9F0AEAAB6LoAMAADwWQQcAAHgsgg4AAPBYBB0AAOCxCDoAAMBjEXSAq5TNZtO8efNKu4wywdWxGDt2rGJiYtxWT2EiIyM1derUv3SbgCcg6AAeKD09XcOGDVPdunXl6+ur8PBw9ejRw+k5Nlb64YcfZLPZCjwYFNZZs2aNHnnkkdIuA7jilCvtAgBYa/fu3WrTpo1CQkL08ssvq3Hjxjp37pwWLFig+Ph4/fbbb6Vd4kUZY5SXl6dy5S7/T9O5c+dUvnz5v6CqsiE0NLS0SwCuSBzRATzMkCFDZLPZtHr1at11111q0KCBGjVqpISEBP3888+FfqawIzIpKSmy2WzavXu3JGnPnj3q0aOHKlWqpMDAQDVq1Ejfffeddu/erU6dOkmSKlWqJJvNpv79+0uS7Ha7JkyYoDp16sjf31/R0dH64osvCmx3/vz5atasmXx9fbV8+fIC9e3evVs2m02fffaZOnToID8/P3300Uc6duyY+vbtq5o1ayogIECNGzfWJ5984vTZjh07avjw4XryySdVuXJlVa9eXWPHjr3kGCYmJqpGjRrauHHjJfu9++67Cg8PV0BAgO655x5lZWVJkpYtW6by5csrPT3dqf/IkSPVrl27QtdljNHYsWNVu3Zt+fr6KiwsTMOHD3csP//U1YwZM2Sz2Qq8zt+vf/3rX7r++uvl5+enhg0b6u23377kvgAeq1SftAXAUseOHTM2m8289NJLl+0rycydO9cY8/8P+zt+/Lhjef4Tjnft2mWMMaZ79+6ma9euZuPGjWbHjh3m66+/Nj/++KPJzc01c+bMMZLMli1bzMGDB01mZqYxxpjx48ebhg0bmuTkZLNjxw4zffp04+vra3744Qen7TZp0sQsXLjQbN++3Rw7dqxArbt27TKSTGRkpJkzZ47ZuXOnOXDggPn999/Nyy+/bNavX2927Nhh3njjDePt7W1WrVrl+GyHDh1MUFCQGTt2rNm6dauZOXOmsdlsZuHChQXGwm63m6FDh5rIyEizbdu2i45dYmKiCQwMNLfccotZv369+fHHH039+vXN/fff7+jToEEDM3nyZMf7s2fPmipVqpgPPvig0HXOnj3bBAUFme+++87s2bPHrFq1yunhoBEREea1114zxhhz6tQpc/DgQcfrk08+MeXKlXPs04cffmhq1KjhGKs5c+aYypUrmxkzZlx0nwBPRdABPMiqVauMJPPll19etq+rQadx48Zm7Nixha6rsM//8ccfJiAgwPz0009OfQcNGmT69u3r9Ll58+Zdstb8oDN16tTL7lf37t3N448/7njfoUMH07ZtW6c+LVq0ME899ZTjvSQze/Zsc//995vrr7/e/P7775fcRmJiovH29nbqN3/+fOPl5eV46vukSZPM9ddf71g+Z84cU6FCBXPy5MlC1/nqq6+aBg0amLNnzxa6/Pygc77t27ebypUrO4WqevXqmY8//tip3wsvvGBat259yf0CPBGnrgAPYoxx27qHDx+u8ePHq02bNkpMTLzsaZ3t27fr1KlT6tq1qypUqOB4/ec//9GOHTuc+jZv3rxINVzYLy8vTy+88IIaN26sypUrq0KFClqwYIH27t3r1K9JkyZO72vUqKHDhw87tY0aNUqrVq3SsmXLVLNmzcvWUrt2bad+rVu3lt1u15YtWyRJ/fv31/bt2x2nC2fMmKF77rlHgYGBha7v7rvv1unTp1W3bl0NHjxYc+fOVW5u7iVryMrK0t/+9jd1795dTzzxhCQpJydHO3bs0KBBg5zGffz48QXGHbgaEHQAD3LttdfKZrO5POHYy+vPfwrOD0rnzp1z6vPwww9r586devDBB7Vp0yY1b95cb7755kXXefLkSUnSt99+q5SUFMcrLS3NaZ6OpIv+8r/Qhf1efvllvf7663rqqae0dOlSpaSkKDY2VmfPnnXqd+GkZZvNJrvd7tTWtWtX7d+/XwsWLChSLZdTtWpV9ejRQ9OnT9ehQ4c0f/58DRw48KL9w8PDtWXLFr399tvy9/fXkCFD1L59+wJ/Dvny8vJ07733KigoSO+9956jPX/c33//fadxT01NvegcLcCTcdUV4EEqV66s2NhYJSUlafjw4QWCQWZmpkJCQgp8Lv+KnoMHD6pSpUqS/pyMfKHw8HA99thjeuyxxzRmzBi9//77GjZsmHx8fCT9+cs3X1RUlHx9fbV371516NDBoj10tmLFCvXs2VMPPPCApD8nP2/dulVRUVEur+uOO+5Qjx49dP/998vb21v33XffJfvv3btXBw4cUFhYmCTp559/lpeXl6677jpHn4cfflh9+/ZVrVq1VK9ePbVp0+aS6/T391ePHj3Uo0cPxcfHq2HDhtq0aZOaNm1aoO+oUaO0adMmrV27Vn5+fo72atWqKSwsTDt37lS/fv1cGQLAIxF0AA+TlJSkNm3aqGXLlho3bpyaNGmi3NxcLVq0SO+88442b95c4DP169dXeHi4xo4dqxdffFFbt27Vq6++6tRn5MiR6tatmxo0aKDjx49r6dKluv766yVJERERstls+uabb3T77bfL399fFStW1D/+8Q+NGjVKdrtdbdu2VVZWllasWKGgoCDFxcWVeF+vvfZaffHFF/rpp59UqVIlTZkyRYcOHSpW0JGk3r17a9asWXrwwQdVrlw59enT56J9/fz8FBcXp1deeUXZ2dkaPny47rnnHlWvXt3RJzY2VkFBQRo/frzGjRt3yW3PmDFDeXl5atWqlQICAvThhx/K399fERERBfpOnz5db7/9tubOnSubzea4uiv/NNXzzz+v4cOHKzg4WLfddpvOnDmjtWvX6vjx40pISCjW2ABXrNKeJATAegcOHDDx8fEmIiLC+Pj4mJo1a5o77rjDLF261NFH501GNsaY5cuXm8aNGxs/Pz/Trl07M3v2bKfJyEOHDjX16tUzvr6+JjQ01Dz44IPm6NGjjs+PGzfOVK9e3dhsNhMXF2eMMcZut5upU6ea6667zpQvX96Ehoaa2NhY8+OPPxpjCp/EXJj8ycjr1693aj927Jjp2bOnqVChgqlatar55z//aR566CHTs2dPR58OHTqYESNGOH2uZ8+ejhoLG4vPPvvM+Pn5mTlz5hRaT2JioomOjjZvv/22CQsLM35+fqZPnz4mIyOjQN9nn33WeHt7mwMHDlxyH+fOnWtatWplgoKCTGBgoLnpppvM4sWLHcvPn4wcFxdnJBV4JSYmOvp/9NFHJiYmxvj4+JhKlSqZ9u3bF2mSOuBpbMa4cfYiAFzlBg0apCNHjuirr74q7VKAqxKnrgDADbKysrRp0yZ9/PHHhBygFBF0AMANevbsqdWrV+uxxx5T165dS7sc4KrFqSsAAOCxuI8OAADwWAQdAADgsQg6AADAYxF0AACAxyLoAAAAj0XQAQAAHougAwAAPBZBBwAAeCyCDgAA8Fj/B1zqIWM2Vz4DAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "values, counts = np.unique(df[\"groups\"], return_counts=True)\n",
    "print(\"Number of cluster:\", len(values))\n",
    "plt.bar(range(len(counts)), sorted(counts))\n",
    "plt.xlabel(\"Cluster rank by size\")\n",
    "plt.ylabel(\"Number of molecules in cluster\")\n",
    "plt.yscale(\"log\")\n",
    "_ = plt.title(\"Sorted distribution of cluster sizes\")\n",
    "# plt.savefig(\"cluster_sizes.pdf\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e092303b-a165-4967-8049-6d340335afb3",
   "metadata": {},
   "source": [
    "### Analyzing failed molecules\n",
    "\n",
    "For the BBBP data sets some molecules fail the processing by the clustering pipeline.\n",
    "Analysing the molecules that failed is an optional step but can give interesting insights."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "642c53ff-3e61-42b4-ac3a-fb2e4f1f1883",
   "metadata": {},
   "source": [
    "Since we replaced all failed molecules with np.nan we can count them with numpy's isnan function "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7e53e730-328e-42d4-a719-c054b83bad4a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df[\"groups\"].isna().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2de6b36-5ef9-4c0f-b718-bbe4c6a86598",
   "metadata": {},
   "source": [
    "Using the nan values in the `df[\"groups\"]` array we can map back to the original molecules in the BBBP data set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "52784f83-0de3-49f6-94a1-81cd459b9258",
   "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",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>59</th>\n",
       "      <td>O=N([O-])C1=C(CN=C1NCCSCc2ncccc2)Cc3ccccc3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>61</th>\n",
       "      <td>c1(nc(NC(N)=[NH2])sc1)CSCCNC(=[NH]C#N)NC</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>391</th>\n",
       "      <td>Cc1nc(sc1)\\[NH]=C(\\N)N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>614</th>\n",
       "      <td>s1cc(CSCCN\\C(NC)=[NH]\\C#N)nc1\\[NH]=C(\\N)N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>642</th>\n",
       "      <td>c1c(c(ncc1)CSCCN\\C(=[NH]\\C#N)NCC)Br</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>645</th>\n",
       "      <td>n1c(csc1\\[NH]=C(\\N)N)c1ccccc1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>646</th>\n",
       "      <td>n1c(csc1\\[NH]=C(\\N)N)c1cccc(c1)N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>647</th>\n",
       "      <td>n1c(csc1\\[NH]=C(\\N)N)c1cccc(c1)NC(C)=O</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>648</th>\n",
       "      <td>n1c(csc1\\[NH]=C(\\N)N)c1cccc(c1)N\\C(NC)=[NH]\\C#N</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>649</th>\n",
       "      <td>s1cc(nc1\\[NH]=C(\\N)N)C</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>685</th>\n",
       "      <td>c1(cc(N\\C(=[NH]\\c2cccc(c2)CC)C)ccc1)CC</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              smiles\n",
       "59        O=N([O-])C1=C(CN=C1NCCSCc2ncccc2)Cc3ccccc3\n",
       "61          c1(nc(NC(N)=[NH2])sc1)CSCCNC(=[NH]C#N)NC\n",
       "391                           Cc1nc(sc1)\\[NH]=C(\\N)N\n",
       "614        s1cc(CSCCN\\C(NC)=[NH]\\C#N)nc1\\[NH]=C(\\N)N\n",
       "642              c1c(c(ncc1)CSCCN\\C(=[NH]\\C#N)NCC)Br\n",
       "645                    n1c(csc1\\[NH]=C(\\N)N)c1ccccc1\n",
       "646                 n1c(csc1\\[NH]=C(\\N)N)c1cccc(c1)N\n",
       "647           n1c(csc1\\[NH]=C(\\N)N)c1cccc(c1)NC(C)=O\n",
       "648  n1c(csc1\\[NH]=C(\\N)N)c1cccc(c1)N\\C(NC)=[NH]\\C#N\n",
       "649                           s1cc(nc1\\[NH]=C(\\N)N)C\n",
       "685           c1(cc(N\\C(=[NH]\\c2cccc(c2)CC)C)ccc1)CC"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.loc[df[\"groups\"].isna(), [\"smiles\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7916de94-248d-4a9f-815f-39ef0126574f",
   "metadata": {},
   "source": [
    "These molecules contain uncharged tetravalent nitrogen atoms which is not correct. Therefore, these SMILES represent invalid structures and RDKit can't read them. See [Pat Walters great blog post](http://practicalcheminformatics.blogspot.com/2023/08/we-need-better-benchmarks-for-machine.html) for more details. Molpipeline can automatically handle these cases through it's error handling. Without the error handling the execution of the pipeline would fail and we would have to make manually interventions.   "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67cafdb7-4f3c-419a-9727-c24f8f7c91ff",
   "metadata": {},
   "source": [
    "### Group clusters into train/valid/test split"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d73a32b5-1d9e-478f-b097-c1541eb54fe6",
   "metadata": {},
   "source": [
    "We now will group the computed clusters into train/valid/test splits that we can use with Sklearn's GridSearchCV. The split ratio we choose is 80/10/10."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6fc46a83-2c7f-41c0-ad08-b0155c33766e",
   "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>name</th>\n",
       "      <th>p_np</th>\n",
       "      <th>smiles</th>\n",
       "      <th>groups</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Propanolol</td>\n",
       "      <td>1</td>\n",
       "      <td>[Cl].CC(C)NCC(O)COc1cccc2ccccc12</td>\n",
       "      <td>993.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Terbutylchlorambucil</td>\n",
       "      <td>1</td>\n",
       "      <td>C(=O)(OC(C)(C)C)CCCc1ccc(cc1)N(CCCl)CCCl</td>\n",
       "      <td>1001.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>40730</td>\n",
       "      <td>1</td>\n",
       "      <td>c12c3c(N4CCN(C)CC4)c(F)cc1c(c(C(O)=O)cn2C(C)CO...</td>\n",
       "      <td>671.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>24</td>\n",
       "      <td>1</td>\n",
       "      <td>C1CCN(CC1)Cc1cccc(c1)OCCCNC(=O)C</td>\n",
       "      <td>828.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>cefoperazone</td>\n",
       "      <td>1</td>\n",
       "      <td>CCN1CCN(C(=O)N[C@@H](C(=O)N[C@H]2[C@H]3SCC(=C(...</td>\n",
       "      <td>546.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2045</th>\n",
       "      <td>licostinel</td>\n",
       "      <td>1</td>\n",
       "      <td>C1=C(Cl)C(=C(C2=C1NC(=O)C(N2)=O)[N+](=O)[O-])Cl</td>\n",
       "      <td>616.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2046</th>\n",
       "      <td>ademetionine(adenosyl-methionine)</td>\n",
       "      <td>1</td>\n",
       "      <td>[C@H]3([N]2C1=C(C(=NC=N1)N)N=C2)[C@@H]([C@@H](...</td>\n",
       "      <td>1028.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2047</th>\n",
       "      <td>mesocarb</td>\n",
       "      <td>1</td>\n",
       "      <td>[O+]1=N[N](C=C1[N-]C(NC2=CC=CC=C2)=O)C(CC3=CC=...</td>\n",
       "      <td>324.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2048</th>\n",
       "      <td>tofisoline</td>\n",
       "      <td>1</td>\n",
       "      <td>C1=C(OC)C(=CC2=C1C(=[N+](C(=C2CC)C)[NH-])C3=CC...</td>\n",
       "      <td>715.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2049</th>\n",
       "      <td>azidamfenicol</td>\n",
       "      <td>1</td>\n",
       "      <td>[N+](=NCC(=O)N[C@@H]([C@H](O)C1=CC=C([N+]([O-]...</td>\n",
       "      <td>1001.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1866 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                   name  p_np  \\\n",
       "0                            Propanolol     1   \n",
       "1                  Terbutylchlorambucil     1   \n",
       "2                                 40730     1   \n",
       "3                                    24     1   \n",
       "5                          cefoperazone     1   \n",
       "...                                 ...   ...   \n",
       "2045                         licostinel     1   \n",
       "2046  ademetionine(adenosyl-methionine)     1   \n",
       "2047                           mesocarb     1   \n",
       "2048                         tofisoline     1   \n",
       "2049                      azidamfenicol     1   \n",
       "\n",
       "                                                 smiles  groups  \n",
       "0                      [Cl].CC(C)NCC(O)COc1cccc2ccccc12   993.0  \n",
       "1              C(=O)(OC(C)(C)C)CCCc1ccc(cc1)N(CCCl)CCCl  1001.0  \n",
       "2     c12c3c(N4CCN(C)CC4)c(F)cc1c(c(C(O)=O)cn2C(C)CO...   671.0  \n",
       "3                      C1CCN(CC1)Cc1cccc(c1)OCCCNC(=O)C   828.0  \n",
       "5     CCN1CCN(C(=O)N[C@@H](C(=O)N[C@H]2[C@H]3SCC(=C(...   546.0  \n",
       "...                                                 ...     ...  \n",
       "2045    C1=C(Cl)C(=C(C2=C1NC(=O)C(N2)=O)[N+](=O)[O-])Cl   616.0  \n",
       "2046  [C@H]3([N]2C1=C(C(=NC=N1)N)N=C2)[C@@H]([C@@H](...  1028.0  \n",
       "2047  [O+]1=N[N](C=C1[N-]C(NC2=CC=CC=C2)=O)C(CC3=CC=...   324.0  \n",
       "2048  C1=C(OC)C(=CC2=C1C(=[N+](C(=C2CC)C)[NH-])C3=CC...   715.0  \n",
       "2049  [N+](=NCC(=O)N[C@@H]([C@H](O)C1=CC=C([N+]([O-]...  1001.0  \n",
       "\n",
       "[1866 rows x 4 columns]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "TEST_SET_SIZE = 0.1  # 10%\n",
    "VALID_SET_SIZE = 0.1  # 10%\n",
    "\n",
    "# create a mask to remove invalid samples from all data\n",
    "df_has_grouping = df.query(\"groups.notna()\")\n",
    "\n",
    "# we first do a 90/10 split to create the hold-out test set\n",
    "grouper_test = GroupShuffleSplit(n_splits=1, test_size=0.10, random_state=RANDOM_STATE)\n",
    "\n",
    "# test_indices are the indices of the hold-out test set.\n",
    "tmp_indices, test_indices = next(\n",
    "    grouper_test.split(\n",
    "        df_has_grouping[\"smiles\"],\n",
    "        df_has_grouping[\"p_np\"],\n",
    "        groups=df_has_grouping[\"groups\"],\n",
    "    )\n",
    ")\n",
    "\n",
    "# calculate what 10% of the full data set is on the 90% chunk\n",
    "propotion_valid_set = df[\"smiles\"].shape[0] * VALID_SET_SIZE / len(tmp_indices)\n",
    "\n",
    "# then, we split the 90% of the data again to create a train and a valid set.\n",
    "grouper_valid = GroupShuffleSplit(\n",
    "    n_splits=1, test_size=propotion_valid_set, random_state=RANDOM_STATE\n",
    ")\n",
    "\n",
    "# the 10% as hold-out test set\n",
    "df_test = df_has_grouping.iloc[test_indices, :]\n",
    "\n",
    "# the 90% containing both train and valid set (will be used later with the grouper_valid)\n",
    "df_train_and_valid = df_has_grouping.iloc[tmp_indices, :]\n",
    "df_train_and_valid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c18106f8-85ac-4512-addf-853db01c8c84",
   "metadata": {},
   "source": [
    "### Hyperparameter search and model selection\n",
    "\n",
    "Here we will use Sklearn's GridSearchCV with Molpipeline to determine the best hyperparameter configuration on BBBP using the train and valid data sets. Afterwards we can determine the best models' performance on the hold-out test set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb5c49cf-3e01-432e-aa5c-5bbda44aa964",
   "metadata": {},
   "source": [
    "Before starting the actual hyperparameter search we pre-compute the standardized molecules. While this pre-computation is not strictly necessary this gives a speed up of about 5-6x for the whole hyperparameter search. The slowest step in the standardization is finding a canonical tautomer. By pre-computation the canonical tautomer is determined once and not for each hyperparameter combination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b55094d2-d07d-43b6-87f9-026bc2f87caf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 463 ms, sys: 123 ms, total: 586 ms\n",
      "Wall time: 7.58 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# Let's pre-compute the molecular standardization\n",
    "pipeline_standardization = Pipeline(\n",
    "    [(\"auto2mol\", AutoToMol()), *make_standardization_elements()],\n",
    "    n_jobs=N_JOBS,\n",
    ")\n",
    "\n",
    "stand_mols = pipeline_standardization.transform(df_train_and_valid[\"smiles\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2bab5bd1-5cf8-4bb3-a01d-912b0b9b1dc1",
   "metadata": {},
   "source": [
    "Let's do the hyperparameter search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "91e40fb7-79ca-4487-81cb-887bb2121557",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/joblib/externals/loky/process_executor.py:752: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
      "  warnings.warn(\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n",
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 27s, sys: 3.24 s, total: 1min 30s\n",
      "Wall time: 4min 49s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "# setup the estimators and their hyperparameters we want to consider in the grid search\n",
    "estimator_dict = {\n",
    "    \"RF\": {\n",
    "        \"model\": RandomForestClassifier(),\n",
    "        \"hyperparameters\": {\n",
    "            \"n_estimators\": [200, 1000],\n",
    "            \"max_depth\": [8, 32, 128, None],\n",
    "            \"min_samples_split\": [2, 5, 7],\n",
    "            \"min_samples_leaf\": [1, 3, 5],\n",
    "            \"random_state\": [RANDOM_STATE],\n",
    "        },\n",
    "    },\n",
    "    \"KNN\": {\n",
    "        \"model\": KNeighborsClassifier(),\n",
    "        \"hyperparameters\": {\n",
    "            \"n_neighbors\": [3, 5],\n",
    "            \"metric\": [\"jaccard\"],\n",
    "            \"weights\": [\"uniform\", \"distance\"],\n",
    "        },\n",
    "    },\n",
    "    \"Logit\": {\n",
    "        \"model\": LogisticRegression(),\n",
    "        \"hyperparameters\": {\n",
    "            \"C\": [0.001, 0.1, 1, 10, 100, 1000],\n",
    "            \"max_iter\": [1000],\n",
    "            \"random_state\": [RANDOM_STATE],\n",
    "        },\n",
    "    },\n",
    "}\n",
    "\n",
    "# Molpipline lets you also natively include hyperparameters of the descriptor\n",
    "morgan_grid_spaces = {\"morgan__n_bits\": [1024, 2048], \"morgan__radius\": [1, 2]}\n",
    "\n",
    "# we loop over all estimators and do a grid search separately\n",
    "cv_dict = {}\n",
    "for model_name, model_dict in estimator_dict.items():\n",
    "\n",
    "    pipeline_estimator = Pipeline(\n",
    "        [\n",
    "            (\"auto2mol\", AutoToMol()),\n",
    "            (\"morgan\", MolToMorganFP(return_as=\"dense\")),\n",
    "            (\"estimator\", model_dict[\"model\"]),\n",
    "        ],\n",
    "        n_jobs=N_JOBS,\n",
    "    )\n",
    "\n",
    "    pipeline_param_grid = {\n",
    "        f\"estimator__{k}\": v for k, v in model_dict[\"hyperparameters\"].items()\n",
    "    }\n",
    "    pipeline_param_grid.update(morgan_grid_spaces)\n",
    "\n",
    "    cv = GridSearchCV(\n",
    "        estimator=pipeline_estimator,\n",
    "        param_grid=pipeline_param_grid,\n",
    "        cv=grouper_valid,\n",
    "        scoring=\"roc_auc\",\n",
    "        n_jobs=N_JOBS,\n",
    "    )\n",
    "\n",
    "    cv.fit(\n",
    "        X=stand_mols,\n",
    "        y=df_train_and_valid[\"p_np\"],\n",
    "        groups=df_train_and_valid[\"groups\"],\n",
    "    )\n",
    "\n",
    "    cv_dict[model_name] = cv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "146f0dbf-4f50-4b34-97e5-b1fb2eea8090",
   "metadata": {},
   "source": [
    "Let's have a look at the best estimators selected with the grid search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "64ec2757-fe0a-4a5d-b589-99231fd8c713",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Pipeline(n_jobs=16,\n",
       "          steps=[('auto2mol', AutoToMol()),\n",
       "                 ('morgan', MolToMorganFP(radius=1, return_as=dense)),\n",
       "                 ('estimator',\n",
       "                  RandomForestClassifier(max_depth=32, min_samples_leaf=3,\n",
       "                                         n_estimators=200,\n",
       "                                         random_state=123456))]),\n",
       " Pipeline(n_jobs=16,\n",
       "          steps=[('auto2mol', AutoToMol()),\n",
       "                 ('morgan', MolToMorganFP(radius=1, return_as=dense)),\n",
       "                 ('estimator',\n",
       "                  KNeighborsClassifier(metric='jaccard', weights='distance'))]),\n",
       " Pipeline(n_jobs=16,\n",
       "          steps=[('auto2mol', AutoToMol()),\n",
       "                 ('morgan', MolToMorganFP(return_as=dense)),\n",
       "                 ('estimator',\n",
       "                  LogisticRegression(C=0.001, max_iter=1000,\n",
       "                                     random_state=123456))])]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_estimators = [cv.best_estimator_ for cv in cv_dict.values()]\n",
    "best_estimators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10dff86b-caf3-4c71-93ad-d1a944b86bf8",
   "metadata": {},
   "source": [
    "The three estimators achieved the following ROC AUCs on the validation set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "5338c16f-e970-45a7-b82b-5992379bf22c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.9608421052631579, 0.9332631578947368, 0.9372631578947368]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_valid_scores = [cv.best_score_ for cv in cv_dict.values()]\n",
    "best_valid_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15b72672-2999-4ad7-8967-ccc66e786b97",
   "metadata": {},
   "source": [
    "### Determine performance on the hold-out test set\n",
    "\n",
    "After selecting the best predictors on the valid set we want to see an estimate of their performance on the hold-out test set  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c0387335-ea6d-4541-83e3-8acfb8a891e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pre-compute standardized molecules\n",
    "pipeline_standardization = Pipeline(\n",
    "    [(\"auto2mol\", AutoToMol()), *make_standardization_elements()],\n",
    "    n_jobs=N_JOBS,\n",
    ")\n",
    "test_smiles_standardized = pipeline_standardization.transform(df_test[\"smiles\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "fe30a622-de78-4b71-90af-d98ee24d368d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RF valid_roc_auc: 0.96 test_roc_auc: 0.95\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/siegjo/miniforge3/envs/ml24_1/lib/python3.11/site-packages/sklearn/metrics/pairwise.py:2350: DataConversionWarning: Data was converted to boolean for metric jaccard\n",
      "  warnings.warn(msg, DataConversionWarning)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "KNN valid_roc_auc: 0.93 test_roc_auc: 0.91\n",
      "Logit valid_roc_auc: 0.94 test_roc_auc: 0.89\n"
     ]
    }
   ],
   "source": [
    "test_roc_aucs = {}\n",
    "for model_name, cv in cv_dict.items():\n",
    "    best_estimator = cv.best_estimator_\n",
    "    test_predictions = best_estimator.predict_proba(test_smiles_standardized)\n",
    "    test_roc_aucs[model_name] = roc_auc_score(df_test[\"p_np\"], test_predictions[:, 1])\n",
    "    print(\n",
    "        model_name,\n",
    "        f\"valid_roc_auc: {cv.best_score_:.2f}\",\n",
    "        f\"test_roc_auc: {test_roc_aucs[model_name]:.2f}\",\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b004220-01e9-4f72-a16c-1a0544274e2e",
   "metadata": {},
   "source": [
    "As shown above the best model on the validation and test set is RandomForestClassifier. However, the other two simpler models are also good predictors with ROC AUCs of around 0.90. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55fb7c57-61a0-413b-9dc8-d326ff029ec7",
   "metadata": {},
   "source": [
    "### Plot the performances on the validation set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "095b0288-40c2-4808-9658-14b2785d82dc",
   "metadata": {},
   "source": [
    "Let's extract the result on the validation set from all hyperparameter combinations we tried in a DataFrame for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "169cfe81-5c83-45c1-a5db-a4ecde295bdc",
   "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>mean_test_score</th>\n",
       "      <th>estimator</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.954316</td>\n",
       "      <td>RandomForest</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.935368</td>\n",
       "      <td>RandomForest</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.955158</td>\n",
       "      <td>RandomForest</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.936421</td>\n",
       "      <td>RandomForest</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.953684</td>\n",
       "      <td>RandomForest</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>323</th>\n",
       "      <td>0.842947</td>\n",
       "      <td>LogisticRegression</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>324</th>\n",
       "      <td>0.805474</td>\n",
       "      <td>LogisticRegression</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>325</th>\n",
       "      <td>0.780211</td>\n",
       "      <td>LogisticRegression</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>326</th>\n",
       "      <td>0.832842</td>\n",
       "      <td>LogisticRegression</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>327</th>\n",
       "      <td>0.846947</td>\n",
       "      <td>LogisticRegression</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>328 rows × 2 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     mean_test_score           estimator\n",
       "0           0.954316        RandomForest\n",
       "1           0.935368        RandomForest\n",
       "2           0.955158        RandomForest\n",
       "3           0.936421        RandomForest\n",
       "4           0.953684        RandomForest\n",
       "..               ...                 ...\n",
       "323         0.842947  LogisticRegression\n",
       "324         0.805474  LogisticRegression\n",
       "325         0.780211  LogisticRegression\n",
       "326         0.832842  LogisticRegression\n",
       "327         0.846947  LogisticRegression\n",
       "\n",
       "[328 rows x 2 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dict_data = {\"mean_test_score\": [], \"estimator\": []}\n",
    "for model_name, cv in cv_dict.items():\n",
    "    dict_data[\"mean_test_score\"].extend(cv.cv_results_[\"mean_test_score\"])\n",
    "    dict_data[\"estimator\"].extend([model_name] * len(cv.cv_results_[\"mean_test_score\"]))\n",
    "\n",
    "rename_dict = {\n",
    "    \"RF\": \"RandomForest\",\n",
    "    \"KNN\": \"K-NearestNeigbor\",\n",
    "    \"Logit\": \"LogisticRegression\",\n",
    "}\n",
    "\n",
    "df_cv = pd.DataFrame(dict_data)\n",
    "df_cv[\"estimator\"] = df_cv[\"estimator\"].replace(rename_dict)\n",
    "df_cv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9cdc2f7-b81c-449b-a5dd-ef5ca4a4380b",
   "metadata": {},
   "source": [
    "Now, let's plot the performances of all tried hyperparameter configurations of the three models. The red X's are the performances of the best estimators on the hold-out test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6938bc29-a2f2-45a2-a5b5-f5fab7f94f69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxYAAAN5CAYAAACVHK0oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAADAzklEQVR4nOzdd3hUVf7H8c9MkplJDymkEQIEBCkCgnQFEUFQVJoUVxERyyq7C64dC64u7v4UdV1dlLXispa1C2JBVKpIkSLSIdQEEkjvmfv7AzMyZpJMMkkm5f16njzKuefe+51kJpnP3HvOMRmGYQgAAAAAPGD2dgEAAAAAGj+CBQAAAACPESwAAAAAeIxgAQAAAMBjBAsAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAADxGsAAAAADgMYIF8Is2bdrIZDI5fVmtVrVu3VoTJ07UypUr67WejIwM3X777UpMTJTFYpHJZNKQIUPqtQY0DwcPHpTJZFKbNm28XUq9eu2112QymXTDDTc4tXvy/Sj7PXLw4MFaqbEqFT0GVOyRRx6RyWTSI4884u1SgCaHYAH8xsCBAzV16lRNnTpVI0eOlN1u1zvvvKPBgwdr/vz59VbHzTffrBdeeEFms1ljx47V1KlTddlll9Xb+dF0DBkyRCaTSd988423S0E1NNfA1xB98803DeLDnYZSR2V43jZvvt4uAGhobrrpJqdP/woKCnTLLbfojTfe0N13360rrrhC55xzTp3WUFxcrA8++EA2m01btmxRSEhInZ4PwK/i4+P1888/y8/Pz9ulVGnMmDHq16+fQkNDvV0KAHDFAqiKzWbT888/r8DAQJWWlur999+v83MeP35cJSUlio6OJlQA9czPz0+dOnVSUlKSt0upUmhoqDp16qTY2FhvlwIABAvAHUFBQerYsaMklbt3evfu3brllluUlJQkm82m0NBQXXTRRXrzzTddHuvs21JWrlyp0aNHKyoqSmaz2XG/dGJioiQpOTnZaczH2beylJSUaMGCBRowYIBCQ0Nls9nUoUMH/eEPf9DRo0ddnrvsOJL06quvqn///goNDXXcE372JWy73a5//OMfOu+88xQQEKDY2FjdeuutOnXqlCSpsLBQf/nLX9SpUyf5+/srLi5Of/zjH5Wbm1vuvNnZ2Vq4cKHGjh2rDh06KDAwUIGBgerWrZseeOABZWRkuKz37PvVV6xYoeHDh6tFixby9/fX+eefrzfeeKPCn5lhGHr//fd1xRVXKCYmRhaLRTExMRo0aJD+9re/KT8/v9w+Gzdu1LXXXqvWrVvLarUqPDxcI0aM0NKlSys8T2Xy8vL0xBNP6Pzzz1dwcLACAgLUpUsXzZkzR6dPny7X/+zvv2EYeumll9SrVy8FBgYqNDRUw4cP19q1a90+f9ltE99++60k6eKLL3Z6Pr322mvl9qnJefPz8/XUU0+pX79+CgsLk81mU8eOHXX33XcrPT3d7XpffPFFmUymSm/5S09Pl9VqlcVi0cmTJx3tX331lWbOnKkePXooMjJSVqtVrVq10sSJE/XDDz+4XYNU9a0cO3bs0IQJExQZGSl/f3917dpVTz75pEpLSys85o4dO/Twww9r4MCBio+Pl8ViUUREhIYNG6Z33nmnXP8bbrhBbdu2lVT+90DZa1iqeozF+vXrdc011yguLk4Wi0UtW7bU6NGj9eWXX7rsf8MNNzieGwcOHNB1112nmJgYWa1WJSUlac6cOSosLKzwcbpSXFysN998U9dee606deqkkJAQ+fv7q2PHjvrDH/6gY8eOudzv7N+VP/74o8aOHev42Xbu3FlPPfWUDMNwuW9+fr4eeeQRdejQQVarVbGxsZo6daoOHTpUrdrL6rj44oslSd9++63Tz8HVc2T58uUaO3asYmNjHd/zMWPGVPga2rNnj2688Ua1bdtWVqtVQUFBSkxM1OWXX65XX321xnW4Yrfb9dJLL2ngwIEKCwuTn5+fWrZsqe7du2vmzJkuxwaVlJTo3//+t4YMGaLw8HBZrVa1bdtWt912mw4fPuzU193nLZowA4BhGIaRmJhoSDJeffVVl9vbt29vSDL+8Ic/ONreeecdw2azGZKMTp06GWPGjDGGDh1qBAYGGpKMadOmlTvO4MGDDUnG73//e8NsNhudO3c2Jk2aZAwfPtxYvHixMXXqVGPcuHGGJCMwMNCYOnWq4+vnn382DMMwCgoKjGHDhhmSDJvNZowcOdKYOHGikZCQYEgyIiMjjY0bN5Y7tyRDknHHHXcYZrPZGDRokDF58mSjb9++xsGDB40DBw4YkozExERj8uTJhr+/v3HZZZcZV199tdGyZUtDktGzZ08jJyfHGDRokBESEmJceeWVxhVXXGGEhoYakoyRI0eWO+/KlSsNSUZUVJQxaNAgY+LEicbw4cONiIgIQ5LRvn17Iy0trcKfyYMPPmiYTCajV69exqRJk4x+/fo5HsvTTz9dbr+ioiJj7NixhiTDbDYb/fr1MyZPnmxceumlRnx8vCHJOHDggNM+zzzzjGE2mw1JRo8ePYzx48cbgwYNMiwWiyHJmDt3rsvnRUXS09ONHj16GJIc36dx48YZkZGRhiSjbdu25Wo4+/s/depUw8/Pzxg6dKhxzTXXGOecc44hybBarca6devcquHnn382pk6dakRHRxuSjBEjRjg9n1auXOnxeY8ePWp069bNkGSEh4cbw4YNM8aMGeP42bVp08Y4ePCgW/VmZGQY/v7+htlsNo4cOeKyzz/+8Q9DkjF27Fin9qSkJMNisRg9e/Y0rrzySmPs2LFG586dDUmGr6+v8b///a/csV599VVDkjF16lSn9rO/H7+1cuVKx+u7Xbt2xqRJk4xhw4YZfn5+xrhx4xyP+7c/2+nTpzt+T4wYMcKYOHGi0b9/f8dzbtasWU79Fy5cWOHvgbPrregxGIZhvPTSS47j9+zZ05g8ebIxYMAAx2vnkUceKbfP1KlTDUnGH//4RyMkJMRITEw0rrnmGmPYsGGGv7+/Icm4+uqry+1XmcOHDxuSjNDQUKNfv37GhAkTjFGjRhlxcXGO3wt79uwpt1/Z78p7773XsFgsxrnnnmtMmjTJGDx4sOHj4+Oo87dyc3MdvyMCAwONK664wpgwYYIRHR1tREREGNdff70hyXj44Yfdqn/evHnGiBEjDElGdHS008/hzjvvdOp75513On7v9OnTx5gwYYLRt29fw2QyGT4+PsYrr7zi1H/btm1GSEiIIcno2LGjMXbsWGPChAlG//79jaCgIKN79+41qqMi06ZNc/zdGDZsmDF58mRjxIgRRocOHQxJxgcffODUPysryxgyZIghyQgKCjIGDx5sjB8/3ujYsaMhyYiIiDA2bdrk6O/u8xZNF8EC+EVlwWLLli2OP9Blfxi2bt1qWK1Ww2azGe+9955T/4MHDzrebL3++utO28r+WEoynn/+eZe1VPbGxjAM45577jEkGUlJSU5vYIqKihxvYNq2bWsUFhY67Vd23pCQEGPt2rUVnrfs2Ge/IUxLS3P88enWrZvRp08fpzCwf/9+o0WLFoYkY9WqVU7HPXz4sPHVV18ZpaWlTu25ubmOP/K///3vy9VT9jPx8/MzPvnkE6dtZW+oQkNDjby8PKdts2fPdryp/fHHH5222e1246uvvjIyMjIcbcuWLTNMJpMRGRlpfPvtt079t27darRq1cqQZHzzzTflaqzIxIkTDUlG3759nb5P2dnZxsiRIw1JxoABA5z2Ofv7n5iYaOzatcuxraSkxLjxxhsNScbw4cPdrsMwfn3OrVixwuX2mp7XbrcbAwcONCQZ06dPN7KyshzbiouLHW+yLr74Yrdrvfbaaw1Jxrx581xu79mzpyGp3PPhgw8+ME6dOlWu/wcffGD4+voaERER5Z4n1Q0W+fn5jvD+pz/9ySgpKXFs27JliyM0ugoW33zzjbFv375y9e3cudPx/Pr+++/dqsOdx7B161bD19fXMJlMxhtvvOG0benSpY7A/MUXXzhtKwsWkowHHnjA6TFu27bNEarWrFlTYU2/lZWVZXz00Uflfh8VFRUZ9913nyHJGDVqVLn9zv5duWDBAqdty5cvd7xZP3z4sNO2P//5z44Qd/ToUUd7bm6ucdVVVzmO6W6wMAzDWLFihSHJGDx4cIV9XnrpJccHJVu2bHHa9u233xrBwcGGxWIxdu/e7Wgve6P/2GOPlTteXl5eud9H7tRRkeTkZEOS0apVK+P48ePltu/YscNITk52apsyZYohybjiiiuM1NRUp21PP/20Icno0KGD0/PEnectmi6CBfALV8EiIyPDWLJkiZGUlGRIMuLi4oycnBzDMH594/jkk0+6PN769esNSUavXr2c2sv+WA4dOrTCWir7xZyfn28EBQUZkoyPP/643Pbc3FzHJ9T/+c9/nLaV/UF99NFHKz2vJGPJkiXlts+fP9+QZJhMJmPbtm3lts+cObPan+7n5uYavr6+RlRUVLltZT+T2bNnu9y3U6dOhiTju+++c7SlpqY63jRt2LDBrRr69u1rSHL5qbZhnLkyJckYN26cW8dLTk42zGazYTKZyr3BMAzDOHLkiONK1+rVqx3tZ3//Xf1sjx8/7rh6UFRU5FYthlG9YFGd83722WeOKzzFxcXl9istLTW6du1qSHL5fHFl+fLljk9vf+vHH380JBkxMTEuz1eRyZMnu3xOVzdYvPnmm4YkIyEhweX3v+yNlqtgUZkXX3zRkGTcddddbtXhzmMo+4Dht1d2ytxxxx2GJOPSSy91ai8LFr169TLsdnu5/W699dZKf4fURFxcnGE2m52CqWH8+ryt6DFcdtllhiSn4JSXl2cEBwcbkozPPvus3D7Hjx93vPZqM1iUlpY6rsBU9Hvn73//uyHJ6erCqFGjDElOn/p7Ukdlyv4mXXnllW7137Fjh2EymYy4uLhyP5syZfWfHfQJFs0bs0IBvzFt2jRNmzatXHtSUpLee+89BQYGym6367PPPpMkTZw40eVxevfuraCgIG3evFkFBQWy2WxO28ePH1+j+jZs2KCcnByFh4dr9OjR5bYHBARo0qRJevbZZ7VixQpNmTKlXJ+qzu3r66vhw4eXa+/QoYMkqXXr1uratWuF2yu6Z3rNmjVauXKlDh06pLy8PMf90WX3y58+fVotWrQot5+rxylJ5557rnbu3Ok0pmTFihUqKipSr1691KtXr0ofpySlpaVp/fr18vf3r/A8ZVM7rlmzpsrjSdJ3330nu92u888/X+edd1657fHx8RoxYoQ++ugjrVixQgMGDHDa7uvr63KcQUxMjFq0aKHTp08rPT1dMTExbtXjruqed8mSJZKkcePGyde3/J8Ts9msiy66SNu3b9eaNWtcPmd+6+KLL1abNm20a9curV27Vv3793dsK7vf/Prrr3d5vmPHjmnJkiXauXOnMjMzVVJSIkn66aefJEm7du3SqFGjqqyhImVjnK655hqXM0ZNnTpVs2bNqnD/nJwcffbZZ9q8ebPS0tJUVFQk6cxkDWX11ZayWisaezF9+nT985//1MqVK1VaWiofHx+n7VdccYXLe+LPPfdcSapwHFdltmzZouXLl+vAgQPKzc2V3W6XdOYefrvdrr1796pnz57l9qvs9b9s2TKnWjZt2qTs7GxFRkZW+FwePny4Pv7442rXX5nNmzfr2LFjSkpKqvD3jqvfI3369NHSpUt12223ae7cuRo8eHC5vxW1pVOnTgoODtbSpUv1+OOPa8qUKY7xEK4sXbpUhmFo5MiRCg4OdtlnyJAhWrp0qdasWaMrrriiTupG40KwAH5j4MCBat++vSQ5Bt7169dPl112mePNTHp6urKysiRJCQkJVR4zPT1d8fHxTm01neO77I9oZX8QymazqeiPf1Xnjo2NdfnGLSgoSNKZYOFK2R+fgoICp/YTJ05o3LhxWrVqVaXnzcrKchksKjpf2YxZZ58vOTlZ0pk/ou44cOCADMNQfn6+rFZrpX3PHixcGU9/RrGxsRVOdRoSEqLTp0+X+x7Xhuqed//+/ZKkBx98UA8++GClx3b3e1c2EPmRRx5xTDAgnRkA/J///EeSXAb/uXPn6vHHH1dxcXGFxy57zdbUkSNHJFX8c23RooVCQ0OVmZlZbtsnn3yiadOmVTqY3dP6zlbVc7Ds+VdQUKD09HS1bNnSaXt1XnNVyc3N1XXXXacPPvig0n4VPf7q1FL2M6rsd1xlr8uaKnst7Nu3r8pByme/Fu666y6tWrVKX331lS677DL5+fmpe/fuuuiiizRp0iRdcMEFtVZjcHCwXn31VU2bNk1z5szRnDlzFBsb6/j7NmXKFMfv+LMf08svv6yXX37Z7ceE5o1gAfzGb9excKXskzbpzKeUVXH1htXf37/atdWWqs5tNlc+YVxV23/rpptu0qpVq9S/f3/NnTtX3bt3V4sWLRxvYuPi4nT8+PEKZ3ip7vmqo+xnGRQUpHHjxtXZeaqjLh9vbZ637Hs3aNCgKqdm7dKli9vHveGGGzR37ly98847evbZZ+Xv769PPvlEaWlp6tevX7nQ+P777+uRRx5RUFCQ/vnPf2ro0KGKi4uTv7+/TCaT7r//fs2bN6/C51ddO3r0qCZOnKj8/Hzdfffduvbaa9WmTRsFBQXJbDbriy++0IgRI7xWnyu1+Ry877779MEHH6hTp0564okndMEFFygyMlIWi0WSNGDAAK1du9Yrr//aUvZaiImJ0YgRIyrtGxkZ6fj/gIAAffnll/rhhx+0bNkyrVmzRmvWrNGGDRs0f/58/f73v9fzzz9fa3WOGzdOw4YN08cff6yVK1dq9erV+uCDD/TBBx/ooYce0pdffqlu3bo5PaYePXqoe/fulR63b9++tVYjGjeCBVADZdNM5ufn68knn3T6Q1HXyq58HDhwoMI+ZZ80/fYqiTfk5uZq6dKlMpvNWrp0qcLCwsptT0lJqbXzlX26uXPnTrf6l11xMplMeuWVV2rlTUzZ973s5+BKQ/oZ1VTZ9+6qq67Sn//851o7bmJiooYOHarly5fr/fff17XXXuuYGvfGG28s179sutbHH39cN998c7nte/bsqZW6yn5WrqbklKSMjIwKr1bk5+drzJgx+tvf/lZn9Z0tPj5e+/bt0/79+13eglb2/LPZbAoPD6/185+t7Ofz9ttvu7w1sDYff1U/o6q21VTZayEiIsLlNM5VueCCCxxXJ0pKSvThhx/q+uuv1wsvvKDx48c7ppmtDaGhobruuut03XXXSZIOHz6smTNn6qOPPtIdd9zhmJ667DENHDhQ//znP2vt/GjaGv7HAEAD5OPjo0svvVSSXM5BX5fKxm6cOnXK5X3C+fn5euuttySpVv8Y1VRmZqZKS0sVEhJSLlRI0ptvvlmrn9QOHTpUFotFGzdu1KZNm6rsHxcXp/POO0/Z2dlatmxZrdRw0UUXyWw268cff9SWLVvKbT9+/LjjXPXxMyr7ZLhszEFtGTlypCTp3XffrfVP28sCxGuvvabU1FR99tln8vf3dzmmqWxtlbL1X8524sSJCtdsqK7BgwdLOvOad3XLVUXrqlRWn2EYWrx4scv9PPm5ld3PX9Gb3FdeeUWSdOGFF7q87bE2Vfb4P//8c6WlpdXauXr16qWgoCClpaXpiy++KLc9NTXVZXtVqvpZlF2F2bFjh2NMT035+vpq/PjxjisfP/74o9t11ERCQoLmzp1b7lxlr++PP/64Wre+1dXvGzQOBAughh5++GFZLBbdddddev31151ujyqzffv2Wl+p22az6fbbb5ck3XnnnY4xBdKZ+9D/+Mc/KiUlRW3btq3xAPHaFB0drRYtWigjI0OLFi1y2rZu3Trdd999tXq+li1b6rbbbpMkTZgwQdu3b3fabhiGvv76a6dPlh977DFJZ+7d/+STT8od0zAMff/9926/IWndurUmTJggwzB0yy23ON1Xn5ubq5tvvlkFBQUaMGBAuYHbdaFVq1aS5PEbnt+66qqrdMEFF2j9+vWaNm2ay/usT58+rQULFlT7TcbYsWMVFhamr7/+Wo8//rhKSko0btw4lyvRlw0ofumllxwDoqUzoXbq1KkuryLUxPjx4xUfH69Dhw7pvvvuc3rNb9++3fE8qqi+//3vf46B2pJUWlqqhx56qMJJAaKiomSxWJSSkuJ4c+6uP/7xj/L19dWHH35YbrHOL774Qi+++KIk1eqVpoqUPf7nnnvOqX3Xrl269dZba/Vc/v7+jqtWs2bNcvp+5+fn67bbbnO5OGZVyl5De/bscRkq/fz89PDDD8swDI0ZM8bleLLS0lJ9/fXXWrdunaPthRdecDloPyUlRRs2bJDkHMiqqqMymzdv1ttvv+3y8Zf93jv7XD179tS4ceN0+PBhjR071uWVntzcXP3nP/9Ramqqo82T5y2aAG9MRQU0RFUtkOfKO++8YwQEBDjmBh8+fLhx7bXXGiNHjnTMTT9x4kSnfaqa+tMwqp6ur6CgwLjkkksMSYa/v78xatQoY+LEiUbr1q0N/bJokaspD/XLVJg1PW9VUx1WNPXl2dNw9u3b15g8ebIxcOBAw2QyGdddd12Fi4pV1F6mbGrM3/7MCgsLjSuvvNKQzixU1b9/f2PKlCnG8OHDK1wg79lnnzV8fX0d89BffvnlxpQpU4xLL73UsTjgPffc4/ob50JaWprRvXt3Q7+stXH11Vcb48ePN6Kiogyp6gXyKlLV98SVTz/91JBkWCwW44orrjBuvPFGY/r06Y6pbj0579GjRx0LAQYGBhoDBgwwJk2aZIwdO9bo0aOHYyGz/Px8t+stUza1adnX119/7bLf/v37jbCwMEOSER8fb4wbN8648sorjdDQUCM2NtaxDsdvpxetyQJ533zzjeM1n5SUZEyaNMm49NJLDT8/P2Ps2LEuv0/FxcVGr169DOnMImOXX365cc011xiJiYmGn5+fY10aV6+r8ePHO6a4nTx5sjF9+nRj+vTpVT4GwzgzjW3Z+jvnn3++MWXKFMfrTqp8gbyKfg9Wdr6KvPfee45zduvWzZg0aZIxdOhQx0KMZYv2/fZ3YlW/Kx9++GGXP9ecnByjT58+ju/36NGjjQkTJhgxMTE1WiCvTO/evR1TIV977bXG9OnTy/1OuOuuuxzP1y5duhhXXXWVMWnSJGPIkCGO5+i//vUvR/+y3xFt27Y1Ro8ebVx77bXG8OHDHYsRDh06tNzUyu7U4coHH3zg+JsxcOBAY9KkSU6L3VkslnJT9GZlZTn+1lgsFuOCCy4wrrnmGmPChAnGBRdc4Jjau2zx1jJVPW/RdBEsgF/UJFgYxpk3IbNmzTK6du1qBAYGGjabzUhMTDSGDBliPPHEE8bevXud+tdGsDCMM29WXnjhBaNfv36OhZeSkpKMmTNnVrhqsbeChWEYxocffmgMGDDACAsLM4KCgozevXsbL7zwgmG322s9WBjGmcXbFi9e7Fjh28/Pz4iJiTEuvPBC4//+7/9cvtHdtm2bcfPNNxsdOnQwbDabERAQYLRr184YMWKE8Y9//MNpsS135ObmGvPmzTN69OhhBAQEGDabzTj33HON+++/3+VibnUVLAzjzIq4559/vuNN8dnfN0/PW1BQYCxYsMC4+OKLjYiICMPX19do2bKl0aNHD+P22283Pv/882rVWqZs3n3pzGKHrtZVKHPgwAHj2muvNVq3bm1YrVYjMTHRuPXWW42UlJQK34DWJFgYxpnnydixY43w8HDDarUa5557rjFv3jyjuLi4wu9Tdna2cf/99xsdO3Y0bDab0bJlS+Pqq682NmzYUOnrKj093bjllluM1q1bG35+fuVew1W90V+3bp0xfvx4IyYmxrFQ4OWXX15uYbwydREsDMMwvvvuO+OSSy4xIiMjjYCAAKNr167G448/bhQWFlb4O7GmwcIwzrz2HnzwQceK7NHR0ca1115rHDhwoNL9KpOcnGxMmTLFiI2NdXwI4eo5snr1auPaa681EhMTDavVagQHBxvnnHOOcfXVVxv//ve/nV77n376qXHbbbcZPXv2NKKiogyLxWK0atXKGDJkiPH666+7XC/F3Tp+6/jx48YTTzxhjBo1ymjbtq0REBBghISEGJ07dzZuv/12Y+fOnS73Ky0tNRYvXmyMGjXKiI6ONvz8/IyIiAija9euxrRp04wPPvigXJ1VPW/RdJkMowFNQwEAAACgUWKMBQAAAACPESwAAAAAeIxgAQAAAMBjBAsAAAAAHiNYAAAAAPBY3S632YzY7XYdO3ZMwcHBMplM3i4HAAAA8JhhGMrOzlZcXJzM5sqvSRAsasmxY8eUkJDg7TIAAACAWnf48GHH6u8VIVjUkuDgYElnvukhISFergYAAADwXFZWlhISEhzvdStDsKglZbc/hYSEECwAAADQpLhzqz+DtwEAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAAIDHCBYAAAAAPEawAAAAAOAxggUAAAAAjxEsAAAAAHiMYAEAAADAYwQLAAAAAB4jWAAAAADwGMECAAAAgMcIFgAAAAA8RrAAAAAA4DGCBQAAAACPESwAAAAAeIxgAQAAAMBjBAsAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAADxGsAAAAADgMYIFAAAAAI8RLAAAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAAIDHCBZoVux2o1xbqYu2ivY1DPf6VnZcV+3u1gAAANBQ+Xq7AKA+7E7N1uNLftbKPScVaPHVuF6tNKxTlJ74fJe2H81SRKBF1/ZL1B+Gtpevj3PeTk7P1WNLftbXO0/I6mvWVT3idd+oTgqx+bk819Jtx/XMV7u1OzVHCeH+unVwkq7tm6ivdqRq/pe7teN4luJCbbrpwnYKtProhW/2KTk9T+1bBukPl3TQld3j6uNbAgAAUKtMRnU+gkWFsrKyFBoaqszMTIWEhHi7HJwlLadQw5/+Tqdyi6rse8OANnrkyi6Of2cXFOvS+d8pJavAqV+/duF66+b+5fb/emeqpr++Qb99VU0f1Eavrj4ody5MvHhdL43oElN1RwAAgDpWnfe4XLGA1xQUFCg5ObnOz/P21tNuhQpJWvx9ska3MSnI6iNJ+nRnZrlQIUnr9p/Sx6u2qGOUzan96c+OlAsVkvTm2mS3QoUkPbPsJ7XxzXSvczUkJibKZrNV3REAAKAGuGJRS7hiUX27du3SjBkz6vw86UkjlRt9ntv9Y358WZa8NEnS6TYXKzuuj8t+Ebs/UWDaDqe2o71/r1JLcM2LlWQuzlOrH57z6BiuLFy4UB07dqz14wIAgKarOu9xCRa1hGBRffV1xeLDHRn61/dpbvW1+Jj07uS2svmdGWfx6c8Zem6d631fvDpBbVpYndoe+OKYNhzNK9c3wM+kvGL3Xmo9Yv31t8vi3epbHVyxAAAA1cWtUGgUbDZbvXyC/vvEYn26Z6UOn8qvsu8NA9qoe9dzHf/+5OBOSeWDhb+vSSP6lb8Kcq81SpMXrlNx6a8hwmSS/nBJRz391W4Vltid+psknR03fM0m3XNFd3VMiqiyVgAAgIaEYIEmL9jmp3dvGaB/fL1H3+w8oUCrrwpLSnXoN0FjaMeWuv/yzk5ta/aluzxmfomhoxn5ig/zd2rv3SZc/53RT/9csVc/H89Sm4hA3TK4nYZ2ila/pAj98+u92nY0Q63DAzR9UDvZ/Mx66bv92ncyRx1jQnT7kCT1bUeoAAAAjQ+3QtUSboVqXAqKS/XW+kNavvOEAi2+Gnt+vIa7mIlp2qvrtWLXyXLtZpO0+cHhCg1wPeUsAABAU8CtUEAVbH4+umFgW90wsG2l/SZekOAyWAxMDCJUAAAAnIVgAVTisq6xuvPSc/T8N3tVUHxmfITt9AHdcs0Q/Xvlfn2984QCLD4ae34rjeoWW61jf7zlmD7cfFSFJaUadm60pvRtLauvT108DAAAgDrHrVC1hFuhmrbMvGL9dCxTeaeOa+5dMxV69YPanuq8vsWtg5N078hObh1v7ic/6dXVB53aBiRF6M3pfWU2m2qrbAAAAI9U5z2uuZ5qAhq10AA/DWgfqYRQi/LCzykXKiTp3yv3KyWzfPtvHT6Vp9fXHCzXvmZfupbvPFEb5QIAANQ7ggXgBrvd0Le7T2rZ7izlh7d32afEbmjTodNVHmtj8ukKV+H+4eApT8oEAADwGsZYAFU4mpGvqa+s194TOWcaorpU2Dc6xFrhtjItgyvuU9k2AACAhowrFkAV7n1v66+hohJd4kLUKzG8yn79kyLUoWVQufZgq6/G9Kz9FbcBAADqA8ECqMSp3CKt2lt+5W1J8vf7dQanQe0j9e+pvd06Zond0EvX91K/dr+GkI7RwXrtxgsUEcQVCwAA0DhxKxSaleT0XJXYDSVFlb9icLbMvGIdy8xXgJ+PKpo37ZzoIM0c2kFWX5OSWgYrNtTfdUedGaPxzxV79dJ3+5RTWCrpzG1Tcy4/V5d2jlZiRGCNHxMAAEBDQLBAs7D3RLbufGeLthzJlHQmFPx9fHf1SAhz6ldcatejn+zQ2xsOq6jErmCrr2JCbS5ne9pyJFM3vbHB8e8LO0TqqWu6q2WwrVzfv32+Uy9+u9+pLTWrUI8t+VkRQRaCBQAAaPS4FQpNXnGpXVNf+cERKiRpd2qOpr6yXlkFxU59n/lqtxatS1ZRyZnF8LILS5SSWaBAS9UL163ck6bfv7mpXHt+UaneXJtc4X6/DRwAAACNEcECTd7XO0/oaEZ+ufbM/GJ99ONRx78Nw9B/vj/k8hhtowL12NVdddW5oarw3ihJG5JPa2dKllPbyexC5RaVVrjPntRszfzvZr2y6oCyfxN0AAAAGguCBZq8tJzCCre9uvqg7L8sKlFcaigjz/Ub+9O5xfpdv0TFBPtJpspXxj6Z7Xy+mFCbwgMtFfYvNaRPthzTo5/u0FXPr9ap3KJKjw8AANAQESzQ5PVtW/EUsPtP5uqDzUe0YtcJXfPiWlUUGU5mF+jBD7frve1VL2A3++0t+r/Pd6qwpFTHMvJ11/+2KKegxK1a95/M1cKV3BoFAAAaH5NhVHJfB9yWlZWl0NBQZWZmKiQkxNvlNGspmQX6OSVLieEBavfL7E/Xv/y9vtvjetrYunJ5txhtOZKpI6edb8MySfLzMamo1PVLr2t8iD6deWE9VAgAAFC56rzHZVYoNBl2u6GHP/5Ji9cfUukvtzedEx2kf0zqqduGJNV7sFi6LUWuooMhVRgqJCnMv+LbpgAAABoqboVCk/HG2oNatC7ZESqkM7M/XfbsSv1w4JTaRtbvlK41vRR4zQUJtVoHAABAfSBYoMl4Z8ORCrfN/2qPxp0fL3Pl4669ymyS/jC0va7sHuftUgAAAKqNW6HQZPx2TYrfeuqL3TW+ilAf7IbUrVWYt8sAAACoEa5YoMmID/OvdHtDDhVlvt55wtslAAAA1AjBAk1GZWtFNBZB1qpX+AYAAGiICBZoMhIj6ndwdl2IDLJ6uwQAAIAaIVigyehXyUJ4FfFtYKO5V9bzlLgAAAC1hWCBJiPQVv25CO4b1Umtwysfm1GRiEA/je0Zp8SImu3vyum8olo7FgAAQH0iWKDJ6BYfqgBL9cYomEwmzbm8c7lpaM0mVbnuRXpusX7XL1HJ6fmV9quOge0jau1YAAAA9YlggSZj06HTyi8qrdY+6TmFuu/9bbL/Zsoofz8f9WsbruAqroIEWH1rdW2MZdtTtf7Aqdo7IAAAQD0hWKDJeHnlgWpPKfv8in1Kzy1/+1FuUan++8NhZReUVLhvVJBF57QM1qD2kS63d4sPrWY10qFTebrxtR+UnlNY7X0BAAC8iQXymoi8vDwVFTXv+/P3n8yu1/OdzCnSRX9frqMZ5UOAzdesTlFWFRYHaPeJPEmSSdLlXaJUUGLXV7vSKzxuTmGJ/rt2n353AStwWywWBQQEeLsMAADgBoJFE5CXl6ePl36uzNwCb5fiVf7F/pL86vWcR1yECkkqKLHr3R9PSDLk72PILCnAx1DeicOy+BiKsfopo9isArvri4arf9wpv2Nb6q7wRiI00KYrR40gXAAA0AgQLJqAoqIiZeYWyC+ilaz+zfcN2CWhxdqz/pSKfztgwqtMyi89Mwgjt1T6+qR7g8sjIiMV1Krxr8vhicL8PGWmH1FRURHBAgCARoBg0YRY/QPkHxjs7TK8pm2gNPOiAH21+7T2p+Upp8ju7ZJq7EBmqYY2459lmWJvFwAAANxGsECT0irMphv6xEqSjmUW6ulvDqm0IV3AcNPhCm6xAgAAaKiYFQpNVlyoVcM7Nc51IcID6nesCAAAgKcIFmjShnZooX5tQmp1rYn6cGG7MG+XAAAAUC3cCoUmzcds0pVdonRRuzDtPpmnj7enlVsMryEwSTIkhdp8dWnHcPVsxfgKAADQuDS4KxbPP/+82rRpI5vNpr59+2r9+vUV9i0uLtajjz6qpKQk2Ww2de/eXcuWLXPq88gjj8hkMjl9derUyalPQUGBbr/9dkVERCgoKEjjxo1TampqnTw+1J/swhK9vv645izdp79/fUjLd59WXIhVoTb3ZmZqH2lTRV39/ap+6Vh83L9MYkjyMUlJkf7qGR/k9n4AAAANRYMKFm+//bZmz56thx9+WJs2bVL37t01YsQInThxwmX/OXPm6MUXX9Rzzz2nHTt26NZbb9WYMWO0efNmp35dunTR8ePHHV+rVq1y2j5r1ix98sknevfdd/Xtt9/q2LFjGjt2bJ09TtSPf689pi3HchxXKLILS3Uks1CZBaVqE27TXUMTNLBtqNqE29QuwqYwfx9ZfU2KCvLTpR1baG9agQpKXR87v7jqGaem94vViE7hCrK6F2RKDWnTkWy9/kOKuw8RAACgwWhQwWL+/PmaMWOGpk2bps6dO2vBggUKCAjQK6+84rL/okWLdP/992vUqFFq166dbrvtNo0aNUpPPfWUUz9fX1/FxMQ4viIjIx3bMjMz9fLLL2v+/PkaOnSoevXqpVdffVVr1qzRunXr6vTxou4cSM+vdGalg6cKtOtEviID/RQfalXLIIu6xATpii6RmjU4oVZmZTqeVaSCYru6RFdvPYrdJ/OUksWsUAAAoHFpMGMsioqKtHHjRt13332ONrPZrGHDhmnt2rUu9yksLJTNZnNq8/f3L3dFYs+ePYqLi5PNZlP//v01b948tW7dWpK0ceNGFRcXa9iwYY7+nTp1UuvWrbV27Vr169evwnMXFv765i8rK6t6Dxh16nR+SZV9lvyU5nIq2o+2pcmvGrcxVeTDbWk13vd0foliQqwe1wAAAFBfGswVi7S0NJWWlio6OtqpPTo6Wikprm8NGTFihObPn689e/bIbrfryy+/1Pvvv6/jx487+vTt21evvfaali1bpn/96186cOCALrzwQmVnZ0uSUlJSZLFYFBYW5vZ5JWnevHkKDQ11fCUkJNTwkaMuJIRZVVU0qGh9ixK74datTnXFxyTFhxIqAABA49JggkVNPPvss+rQoYM6deoki8WiO+64Q9OmTZPZ/OvDGjlypCZMmKDzzjtPI0aM0NKlS5WRkaF33nnHo3Pfd999yszMdHwdPnzY04eDWhQVZFFCi8b55vzCpDCF2BrMxUQAAAC3NJh3L5GRkfLx8Sk3G1NqaqpiYmJc7hMVFaUPP/xQBQUFSk9PV1xcnO699161a9euwvOEhYXpnHPO0d69eyVJMTExKioqUkZGhtNVi8rOK0lWq1VWa+N849oclJTadTK7YY9TCLL6KDrITwEWH53KK5HV16w+rUPUJzHE26UBAABUW4O5YmGxWNSrVy8tX77c0Wa327V8+XL179+/0n1tNpvi4+NVUlKi9957T1dddVWFfXNycrRv3z7FxsZKknr16iU/Pz+n8+7atUuHDh2q8rxouE7nl6iyYRY239p/6vv/csxAi3vHHpwUptsvTNC0vnG68+LWuuPCVoQKAADQaDWYKxaSNHv2bE2dOlW9e/dWnz599Mwzzyg3N1fTpk2TJF1//fWKj4/XvHnzJEnff/+9jh49qh49eujo0aN65JFHZLfbdffddzuO+ec//1mjR49WYmKijh07pocfflg+Pj6aPHmyJCk0NFTTp0/X7NmzFR4erpCQEM2cOVP9+/evcOA2Gr4Qm6/8zFJFQyUSw23KL7br0OmCWjmfj9mkey5JlNksBVh8tGLPaS37Ob3CcRyStPZgpi45J7xWzg8AAOBtDSpYTJw4USdPntRDDz2klJQU9ejRQ8uWLXMM6D506JDT+ImCggLNmTNH+/fvV1BQkEaNGqVFixY53dJ05MgRTZ48Wenp6YqKitKgQYO0bt06RUVFOfo8/fTTMpvNGjdunAoLCzVixAi98MIL9fa4UfusvmYNbBemb/ZmlNtmknTJOS3UPjJAB0/l63ReiXzN0nf7MnUgPV+usojpl6+KhnRf0DpYIf6/vpwuOSdcfVqHaG9avt7alKpiF8t95xRWsEgGAABAI2QyDKOSz1ThrqysLIWGhiozM1MhIfV7O0tGRobe/ugzBbU6R/6BwfV67obMbhj6/Od0fbMvQ8W/XDqICPDVVd2i1DXW9erWOYUlysgr0U8puVp9MFM5haWKD7Vq5LkRat3Cqsz8UuUWlWjZzlM6eKpA/n5m9W8TqpHnRsjH7HoeqpfXHdNPKbnl2jtHB+im/vG194CbmPzcbOUc2a2JV40sN2sbAACoH9V5j9ugrlgAtclsMmlk50hddm6ESkrtkskkP5/Kxz8EWX0VZPVVqxY2De8UrlK7Id+z9gmy+kqy6pyWgSoptcvHbJLJVPnEtiPPjdD+9HynKWxtvmaNPDeykr0AAAAaF4IFmjyTySQ/X58a7edbyUJ5vlWElDJxoVbdOaS1Vh3IUEp2kaKDLBrULkwRgX7VrgkAAKChIlgA9SA80E9Xdo2quiMAAEAj1WCmmwUAAADQeBEsABfshqGsghKVupjNqSql9jP72g1DJWf9/9ky8oqVllvkOH5Bcalyi5glCgAANF7cCgX8xur9Gfpy9yllFZQq0GLWiFg/XRpvUX7L8iux+59IUXFAoEqCgmUYhpbvPq1v951WbpFdFh+TDEnFpYbC/H01vGO4Qmy+emtzqmOqWR+TFBHop5M5xTIktQ23aVz3looLZVV3AADQuBAs0GzlFZXqp5Rc2Q1DnWMCFWz11cbDWXpv60lHH3N2tn738kNqXZKjd5/4t/ZawxUdYlH7yAD5px7Xhbdfp8IWEVr99L+1PKVES39Od+xbdNbqeBn5JXrnxxMySTr72kWpIZ3IKXb8+8CpAi1YfVT3XZoof7/qDzgHAADwFoIFmqWtx3K0eGOK482/j9mk8d2j9N2+DKd+gYX5isjLVGRGikbNnqZJk+fpeEiUzi0+rdfeuEdBacdkGNJPe1K07Gjl085KzqGiIjlFpdp0OFsD24VV/4EBAAB4CWMs0OzkFZU6hQrpzLiId388odTsIqe+KSGRmjR5npLDYpSYkaK3/nufzj/ysxa8cpei044pOSxGo8b8RS8mS4UltbfW5On8klo7FgAAQH3gigWanIz8En2xK127UvNk8zOrb2KIBrULk9lkUkFxqf7zm1BRpqJx2sdDojRp8jy99d/7lJiRovf/c5ckKTks5swVjKDaX+iudQtbrR8TAACgLnHFAk1KflGpnvvusNYdzNLp/BIdzyrSh9vS9MHWk7IbhhasOaafU/OqfdzjIVGadfmdTm2zLr9Tx0Nqf22KNuE2dY0NrPXjAgAA1CWuWDQR+Xm5yj9+WDZb835Duj6lVKfzy0/buvZgpkKNHB06XfmUri0s0umi8u2xWSf19JKnnNqeXvKUY8yFJ7qEm5RZJJXYpXPCzOodXaq0o4c8OmZTUFCQK+XlersMAADgJoJFE7Fn5w5t27zB22V4XVr7y6WWXcu12w3pm/U/Si27Vbhv2MEVKvYNkFr1dWqPzTrpuA0qOSxGsy6/U08vecox5sLTcHHsh88UnLpFFkkHf/nCGd169vZ2CQAAwE0mwzBqb8RpM5aVlaXQ0FBlZmYqJCSkXs+dkZGh1/77Pyk4stlfsVh1rFSrjrm+KnFJK7OWH7G73BZulUIsJh3Mdn45xGSl6e3/3usIFWUh4rdhY+LkJ5QSUrOxFpPO8VWbEO5K/K2CglwpO003TB6vsLAwb5cDAECzVJ33uFyxaCL8AwIVFJsg/8Bgb5fiVZdElGjjyWTlFzsHiPPigjSiZ7RWpxxQQUn5cHGqUDpVWD5j51r9lR4QKklOVybOHtCdHhCqXKt/jepNCLPqgnMTZDJVPVVtc5Ofm62c0nxvlwEAANxEsECTEurvq98PbKWPtp/Q3rQCWXxM6psYqss7Ryi7sFSGWytJ/CrbGqip1zyqwML8clckjodEaeLkJ5Rr9Ve2tXpXinzNJvWID9KVXaOUV2zXkdMFMpmkhDCb/C0sjAcAABofggWaFLthaNORbB06XSjpzMJ3ZpO0YM1RHTxVUKNjZlsDKwwONbn9KcDPrEdHtVNRiaFFG46Xm6XqvLgg/a5XtHx9uD0KAAA0HrxzQZOyfPdprdh72rFORX6xXd/uy6hxqKgLARazzCaT3tqc6nLq263HcvTBtpNeqAwAAKDmCBZoUlbtz/B2CVUanNRC2QUl2nYsp8I+6w9lqajU9UBzAACAhohggSbDMAzlFFa+TkVVooMtuqprhNqG21ST8dQ+VewTFeSrge3CfhnvUbFSu1RQTLAAAACNB2Ms0GSYTCa1CbfpQDVve+oRF6ghHcLlYzYpPtQqSRrcPlzpucV6akWyCkrcG/AdZPHRPcNaa+lP6VqbnOWyz6XnREiSWgZb5O9nLjd7VZkwf18FWxnEDQAAGg+uWKBJubxzpHzNzpcNXF1F8DFJNl+zzCbpx2O5entzqvKKnK927DyR63aoMEm6qlukAi2+mtAzWpec06Jcn9YtrIoI9NPek3n6Yme6SuwVH3tMtyimoAUAAI0KVyzQpLSL9NefBido5f4MncwpVlyIRYOSwpR8qkCbj2TLbpyZdemn4zn6+cSvA6ePZxXp3+uO6c8Xt1ZUkEWSdDDdvSsfZpM0tU+susUGOdou7xyphDCbfjiUpcISu+LDrNpxPEfPrTxS6XHiQiwac15LtY2o2boYAAAA3kKwQJMTF2rVxJ7RTm0tgyxqGWzRlqPZOpyR7xQqyhSXGlp3MEuju56ZQjbU372Xh92QUrOK1C3Wuf28uCCdFxcku2Hor18e1Km8kkqPk9jCppkXJbh1TgAAgIaGYIFm4fOd6fp856kq+2XkFzv+v1+bEK3cn6Hi0qpvh8rIrzg07DmZV2WoqOoYAAAADR1jLNDkncwp0hduhApJahXq55iNKTLQoun94hQdfObWqMqGPCSG2ySdWaDvt7M55RW5N7tT2TEAAAAaI65YoMn7OTW30qldz/bJjtP6ZMdpRQX6alTnSHWPD9Y9lyTqdF6xbL5mfbHrlL7dl+G0j9XXpNgQP32+M12r9mcot8iu6GCLLu8coa6xQWof6S8fs0mllQzWtvqaNeyc8Jo/SAAAAC8jWKDJs/hUfGEuwM+kvOLyb/hP5pbo9R9SNCyzUK3CbOrUMkAWX7Ou7BqpbcdznG5tKiwx9PzKoyo865ap1Owivbb+uH4/qJVyi0oVH2rRodOFTucItfkqyOqj+FCrLu7QwnFlBAAAoDEiWKDJOppZqG3HclRqGPIzS79dMsLiYyrX9ltf7T4tSfL3M+v6C2Lkaza7HC9R6GIcht2QXlx9xOkcZpPUKSpAHaMDlF9sl92QusUFESoAAECjR7BAk/TZz+n6ctev4ypMOhMkis4KAH4+JuW6Of4hv9iuhWuPqXN0YLXq+G1wsRvSvlNnZqUqq+SLXad0acdwjTw3olrHBgAAaEgYvI0mpdRu6KNtJ51ChSQZOjOdrNX31xHY7oaKMnZD2p6S63GNhSVGuTEfX+46pdfWH9Op3GKX+wAAADR0BAs0KW9uSCk3uLqMoTNv6huqrcdy9ex3h5l2FgAANEoECzQZRzMLteVYTr2cK9hiVnyoVdHBFrUM9qu142YXlmrl/oxaOx4AAEB9IVigyTh8uqDezpVTZFf/NqG655JElZbW7rGP1OPjAAAAqC0M3kaTER5Ye1cOqmJIem/LCQVZzCo1avf2qhYB9fc4AAAAagvBogkpzM/zdgle1crfUGywr45n188YBUPSaz+k1OoxfUzSBXF+ys/NrtXjNkbN/fkMAEBjQ7BoAiwWi0IDbcpMP6LmPqfQlS1NWpgdJLtMVfQ0pCr71K9Q31JdFlugkMws5WR6u5qGITTQJouFNT4AAGgMCBZNQEBAgK4cNUJFRUXeLsXrXlt3RPZ9yZX28TNLr/+uu2579ydl5tfyAIka8vcza+kdfRVo5SV5NovFooCAAG+XAQAA3MC7mCYiICCAN2CS3tzwQ5V9iu3SjLd+Um5RwwgVkjT3qq6Kj470dhkAAAA1xqxQaDJ+Pp6ltBz3rto0pFAxqH2krumd4O0yAAAAPMIVCzQZxzPzvV1CjZzfOszbJQAAmqhTBae0+OfF2nxis1oGtNTEjhPVo2UPb5eFJopggSajW3yY/HxMKi5tuKtr/1aYv6+m9E30dhkAgCYoPT9d1y69Vkdzjjralh5Yqr9d+Ddd1vYyL1aGpopbodBkRAVbNX1QO2+XUS0zh3ZQTKjN22UAAJqgN39+0ylUSJLdsGv+xvkqtTecW4LRdHDFAk3KvSM76dzYYL274Ygy8oq043iW7HV4ASPA4qM8D8ZrvLbmoIZ0jNKmQxmKDfXXgKQImc0NaxpcAEDdKigoUHJy5TMa1sSqg6tcth/PPa5V21YpxhZT6+dsLBITE2Wz8cFebTMZRi0vG9xMZWVlKTQ0VJmZmQoJCfF2OfjF+5uO6K53t6gu7o4KsflqQFKklv3k2SJ5Jp1ZVUOSOrQM0ms39lF8mL/H9QEAGoddu3ZpxowZtX7cnKE5Km7rYoWrUil0cajMRc33xpWFCxeqY8eO3i6jUajOe1yCRS0hWDRcxzPzNXPxZm1IPl07B7SX6KY+0frDFb1033vbtGTb8do57i8u7BCpRdP71uoxAQANV11dsdiWuU2P/vyoDDm/1RsSOUR3tL/D4+MnJyfrscce05w5c5SY2LjGC3LFwn3VeY/LrVBo8mJD/fXC785Xn8eXe3ysVqF+Klq9WBOmP6TtRzNdhoqzr0DUxKq9aUrLKVRkkNWDowAAGgubzVYnn553VEf5tPDRs5ueVXpBunxNvrqs7WV6sN+DCvCrvbWvEhMT+fQfkggWaCbe3XDErX4PjOqkl1cfVEpmQbltZpO0cExr3fLZYUnSlztSXR7D00uAhiEVl9o9PAoAANKYDmN0RdIVOpx1WOG2cIXZwrxdEpowggWahdV706rsM6FXK53ILlSqi1AhSRZfs0y/+XddOK9VqGJDGWMBAKgdfmY/tQtrXLMmonFqvqN20KyEB1oq3d41PkSxoTYtXHmgwisOV3aPk8n0a7S4qnu8TC4mcHJ3Tqe4UJt6JbZwagsL8NPjV3dz8wgAANRMZmGm3t75thZuXaif0n/ydjloIrhigSZp9d40fbr1mErthkZ2jdWUPq316daKB1nPG3Oepr/+Q4Xbeye20AOjOivl8H5J0r5Thfph92GdFx+q7ceyVPrLnLZWX7NK7Ibj3xVpFxWo16edmf3p650ntPHQacWG2nRV93iFBvjV4BEDAOCe9cfXa+bXM5VXkidJ+sfmf2j8OeP1cP+HvVwZGjuCBZqcJz7bqQXf7nP8+50NR/S7fq018YIEvf3D4XL920QEqGt8iNJyCl0eL9TfT/+7bYAkKUVSblQX3fHxYaf1MUL9/ZSZX6zCkorHRozpGaeB7aMUHWLVwKRIx3oVwzpHa1jn6Bo8UgAAqqfUXqr7V93vCBVl/rf7fxrSaogGJwz2UmVoCrgVCk3KgbRcvfjdvnLtb647pN/1ba2re8Y5tYfYfPXUNd1lMpl0QZtwl8cc1D7S8f8f7chQevtR5Rbdy8x3MU/4WbrEhejh0V00vlcrXdghikXwAABesS1tm1LzXE8+8kXyF/VcDZoarligSVm156QqWpnl/z7fpTem99W0AW21el+awgMsGnVerEJsZ249uvuyTrr23+tUUPzrVYcQm6/+OKyDJOml7/bphe/TJJP7eXxIxyhd2zdRQzu1lA9hAgDgZeZK/ob5mHzqsRI0RQQLNCkh/hWPT9iQfFpFJXZ1TwhT94Swctt7JbbQpzMH6dXVB3UgLVedYkI0bWAbJYQHqLjUrhe/3V/tei7vFqtLuc0JANBAdI3sqvigeB3NOVpu22VtL/NCRWhKCBZoUi7tHC0fs8nl4Om8olKl5RQqLqziqVzbtwzW42PKz8p0KrdI6blF1aol2Oarkd1iq7UPAAB1yWwy6+8X/V13LL9DpwtPS5JMMumGrjdoQNwAL1eHxo5ggSYlwOKrC9tH6pvdJ8ttCw+01Hg16zP7WpSWUz5c+JikhPAAHcsoUNEvC9vFhdr0zKSeCrLyEgMalP3fShtflXLTpDYXSn1mSAGux1cBTdV5Uefp8/Gfa8WhFcoqylL/uP5KDEn0dlloAnjXgybnnpGdtO5AutNYCUm6+aJ2NV7Uzs/HrFsHJ+mxJT87tZtNUqkhHUz/dXaNpKhAfXzHIAUSKoCGZcOr0qd/+vXfB1dK296VbvpK8g/zVlWAV/j7+mtUu1HeLgNNDLNCock5NzZEi2f004UdIhVo8VH7lkF6fExX3To4SYfS87Qx+bTyi0qrfdybLmynPw6Ikl/eSQX4mdS3bbh8XQzI3ncyV59uPVbl8ZLTc7Ux+bQKiqtfC4BqKi6Qlj9avj19z5krGAAAj/GRKpqk81u30KLpfR3/PpVbpOtfWa/vfrlFKsTmq7sv66SB7SP1w8FTahls1YUdoqqcuWlUx1B98OQrWrhwoVJN4br+lfUu+322PUWS1CE6WOe3dl5d+1Rukf741mat3JMm6cwaGHdf1lHX9uUyNOpfQUGBkpOTvV1GnbOe2qk2+adcbsvZ8aWORtXfJ7eJiYmy2Wz1dj4AqC8ECzQLf353iyNUSFJWQYnmfLhdJkllw7zbRQbqtWl91DoiwK1jtgypeLzGN7tO6ptdZ843sH2EXrqut+PWqDvf+dERKqQza2A88MF2JUUFqV+7iOo9MMBDycnJmjFjhrfLqHORlmL9b+CZ2xd/65sNP+vvb9bf92DhwoXq2LFjvZ0PAOoLwQJNXmpWgVbsOuFy29lzR+1Py9Vd/9uit2/p79ZxO8WEqE+bcK0/6PpT0DKr96br6S93a84VnXU8M9/lwHJJevuHwwQL1LvExEQtXLjQ22VUS3Jysh577DHNmTNHiYnuX+nLXXm3go9849RmmHzU9Yb5WhjRuZarrFh1agaAxoRggSbvdF5RhYvm/db3B04pJbNAkvT8ir1atTdNof5+mnRBgib1aV2u/wu/O193vbtF3+w+szCf2aRyq3JL0oc/HtOcKzorI6+4wlpOVXM6W6A22Gy2RvvpeWJiYvVqT3xD+vgP0s+fSEapFNJKpuF/UZuuY+quSABoRggWaPLaRwUpOsSq1KxCt/qn5xTq5kUbdTQj39H24+EMHUjP1dh2zvMdRAZZ9eq0PjqRXaDcwlJd/fxqZeYXlztm8S/T0HZoGaSWwVadyC5fy6D2kdV5WACqyxYqXfO6lJsu5Z+SwttJZlYaBoDawqxQaPJ8fcx68IrOVQ7MlqT2LYO0cm+aU6go8+rqg8oscD2DU8tgm9pGBla4yvaILtGV1tI5NkST+5a/IgKgDgRGSJEdCBUAUMu4YoFm4Yrz4tQmIlBv/XBI3+46qcOnywcHSRrfq5W2Hc10ua2oxK7kjMpvV7p7REdtSj6t/Wm5jrZ2kYH684hfb9cY3T1ObSMD9d/1h5SeU6R+7cJ1zQUJCrDwcgQAAI0X72TQbHSND9Vj8d204Nt9euKznS77DD4nyuWtTJJkMknRQZW/ZFqG2PTZny7UZ9tStDs1W+dEB2tktxhZfZ0/Ge0aH6rHx3Sr2QMBAABogAgWaHYm9GqlF7/dp9N5zgFiUPtInRsboiCrr95Yc1C5v1lEb1TXWEUH+VV5fKuvj67uGV+rNQMAADR0jLFAsxMRZNV/buqn/r9M7Wr1Neua3q30wu/OlyQlhAdo0U191SvxzMJ2ARYfXd8/UU9O6O61mgEAABo6rligWeocF6L/3txP+UWl8vUxyc/HOWOf37qF3rttgPKKSmTxMcvXhwwOAABQGYIFmjV/S+WzwjCgGmgESoulXZ9JmYel+N5S677erggAmiXeNQEAGq/TydKiq6VT+39tO+cy6ZpFkq/Fa2UBQHPE/R0AgMZryZ3OoUKSdi+Tvl/gnXoAoBkjWAAAGiVzcY609yvXG3d8WK+1AAAIFgCApsiwe7sCAGh2CBYAgEbJ7hckJQ11vbHzVfVbDACAYAEAaMQuf1IKbe3cljRU6nubd+oBgGaMWaEAAI1XeDtp5gZp56dSxqEz0822vdDbVQFAs0SwAAA0br5Wqes4b1cBAM0et0IBAAAA8BjBAgAAAIDHuBUKANB45aZJ29+T8jPODNpOuMDbFQFAs0WwAAA0SgEp66X/3SMV555p+OavUs/rpCufk0wm7xYHAM0Qt0IBABodH5OhmHVzfw0VZTYvknYv805RANDMESwAAI3OuSF58ss/6Xrjz5/UbzEAAEkECwBAI2Q3Ktlo9qm3OgAAvyJYAAAanZ+zAlQU1Mr1xq7j67cYAIAkggUAoBEyZNLxAY9JgVG/NprM0oV3Su0Ge68wAGjGmBUKANDwnPhZSv1JiuwgxXZ32aUgorP0p+3S7s9+mW72YqlFm3otEwDwqwZ3xeL5559XmzZtZLPZ1LdvX61fv77CvsXFxXr00UeVlJQkm82m7t27a9ky59lA5s2bpwsuuEDBwcFq2bKlrr76au3atcupz5AhQ2QymZy+br311jp5fPCOU7lFKigubbbnBxqNkkLp7d9JL/ST3psuvXiRtGiMVJjjur+fTeoyRuo9jVABAF7WoILF22+/rdmzZ+vhhx/Wpk2b1L17d40YMUInTpxw2X/OnDl68cUX9dxzz2nHjh269dZbNWbMGG3evNnR59tvv9Xtt9+udevW6csvv1RxcbGGDx+u3FznKQpnzJih48ePO77+/ve/1+ljRf1YvTdNI59dqfP/8qV6PPqF7nt/m/KL6u8N/tc7U3Xp/G91/l++1Pl/+VKPfPyTCksIGECFVj5VflanfV9Lyx/1Tj0AALc1qGAxf/58zZgxQ9OmTVPnzp21YMECBQQE6JVXXnHZf9GiRbr//vs1atQotWvXTrfddptGjRqlp556ytFn2bJluuGGG9SlSxd1795dr732mg4dOqSNGzc6HSsgIEAxMTGOr5CQkDp9rKh7e09ka9prP+jn41mSpIJiu/67/pDufm9rvZx/+9FM3fzGRu05ceaT1ryiUr225qAe/uinejk/0Chtect1+9YK2gEADUaDGWNRVFSkjRs36r777nO0mc1mDRs2TGvXrnW5T2FhoWw2m1Obv7+/Vq1aVeF5MjMzJUnh4eFO7f/5z3/05ptvKiYmRqNHj9aDDz6ogICACo9TWFiowsJCx7+zsrIqfnBwqaCgQMnJyXV2/BfWnVRRib1c+5KtxzSpo0WRgdV/+pfV607d/1iVqhIXc2K+t/GwxnXwVYi1fqfETExMLPd6ARqc4nzX7UV59VsHAKDaGkywSEtLU2lpqaKjo53ao6OjtXPnTpf7jBgxQvPnz9dFF12kpKQkLV++XO+//75KS13famK32/WnP/1JAwcOVNeuXR3tU6ZMUWJiouLi4rR161bdc8892rVrl95///0K6503b57mzp1bg0eKMsnJyZoxY0adHf9kxzFSxDnl2u2GNHvOo7LmHK/xsR977LEq+6R2niiFtSnXXmyXZt41R5b8tBqfvyYWLlyojh071us5gWrrOFLa9Hr5dnuxtP/bupnxyV4qbXtX2vGxZDZLXcedGbcBAKiWBhMsauLZZ5/VjBkz1KlTJ5lMJiUlJWnatGkV3jp1++23a/v27eWuaNx8882O/+/WrZtiY2N1ySWXaN++fUpKSnJ5rPvuu0+zZ892/DsrK0sJCQm18Kiaj8TERC1cuLDOjv/fLaf02qZT5dr9fU1a8PdHFOBXt3cCvrIxXW9vPV2uPcRq1sL5j8viW793IiYmJtbr+ZqrvLw8FRUVebuMRsvU63aFbP6PTEZJuW3F3zyp3PDuysk5c3thTk6OMjIyPD5nwJLfy7Jnya8NP3+iwl3LlT+06g8QmgOLxVLpFXzgbIWlhbL6WN3un5qbqpe2vqQ1x9Yo2BKssR3GamLHiTKZTHVYJepKgwkWkZGR8vHxUWpqqlN7amqqYmJiXO4TFRWlDz/8UAUFBUpPT1dcXJzuvfdetWvXrlzfO+64Q59++qm+++47tWpVwaJKv+jbt68kae/evRUGC6vVKqvV/RcOyrPZbHX6CfofWxXpy/2rdDTD+daKOy45Rz27tq+z85a5M7ZAKw6s0onsQqf2WcM7qVuXtnV+ftS/vLw8Lf/sIxXnZnq7lEbLbC/W1S5ChSQVHdumrz9+S7l5+erbs6t2/PCdkn/6waPzReQf0ODjS8q1W7cu0renWirH0tKj4zcFfoGhumTkVYQLVOqTfZ9owZYFOpR9SLGBsbqx642a1GlSpftkFmbq+s+u17HcY462x79/XMlZybqnzz11XTLqQIMJFhaLRb169dLy5ct19dVXSzpz69Ly5ct1xx13VLqvzWZTfHy8iouL9d577+maa65xbDMMQzNnztQHH3ygb775Rm3bVv2G7scff5QkxcbG1vjxwPtaBFr0we8H6MXv9mv13jSFB1o0uU9rje4eVy/nbxli04e3D9SL3+7T9wdOKSrYqt/1S9SILq6DMhq/oqIiFedmqmtLHwX688FDTZWkRsk3/2S5dp/QOPVNDJAUoKHnRtTKuYIPHq5wW5+Aw8qNb1Mr52mscvMLtf1EpoqKiggWqNC69HV6cs+Tjn8fzz2ux79/XL5mX40/Z3yF+32w5wOnUFHmrV1vaXq36Yr0j6yTelF3GkywkKTZs2dr6tSp6t27t/r06aNnnnlGubm5mjZtmiTp+uuvV3x8vObNmydJ+v7773X06FH16NFDR48e1SOPPCK73a67777bcczbb79dixcv1kcffaTg4GClpKRIkkJDQ+Xv7699+/Zp8eLFGjVqlCIiIrR161bNmjVLF110kc4777z6/yagVrUMsenBKzp77fxxYf6ae1XXqjuiSQn0tyo0yN/bZTRaJR1Hy/dH51taDZOP7B2vqPXvq29QeIXbrCER8uXnKImB86jcx8c/dtn+2k+vVRosfj71s8v2EnuJ9mbsJVg0Qg0qWEycOFEnT57UQw89pJSUFPXo0UPLli1zDOg+dOiQzOZf70svKCjQnDlztH//fgUFBWnUqFFatGiRwsLCHH3+9a9/STqzCN7ZXn31Vd1www2yWCz66quvHCEmISFB48aN05w5c+r88QIAyitJGKB8H4ss+5bJnHtCpaGJKupwhewRHWr/XPEXyNj5nkzFzm+e7bYWKo12veI3AGcpBSku2w9nV3xFUJISgl2PTTXJpFZBld+2joapQQUL6cxYiIpuffrmm2+c/j148GDt2LGj0uMZRvnpPs+WkJCgb7/9tlo1AgDqVmlcb+XH9a77E/kFKL/vn2T78RWZc868OSoNaa3C86dL5gb3JxJokNoGttWWzC3l2juFd6p0v/HnjNfinxcruzjbqf3SxEvVKphg0RjxWxMA0KzZW7RT3sWPyZR9TDKZZQQxDgqojnHx4/RT9k8qsf868YLZZNbvu/++0v1iAmO0cMRCPbXhKf2Q8oMCfAN0VfurNKvXrLouGXWEYAEAgCQjuH4mdgCams4hnfXqiFf18raXtfv0brUJbaNpXaepX2y/KvftEtFFr4x4RUWlRfI1+8psqt+p2FG7CBYAAADwSI+WPfTcJc/VeH+Lj6UWq4G3EAsBAAAAeIxgAQAAgGrZkb5D/9j7D2VdkaWFBxbqUNYhb5eEBoBgAQAAALetPbZWv1v6O32X9p1Ko0v1eernmvTpJO3L2Oft0uBlBAsAAAC47R+b/qFie7FTW3Zxtl7c8qKXKkJDQbAAAACAW4rtxdqevt3lts0nN9dzNWhoCBYAAABwi5/ZT+G2cJfbWga0rOdq0NAQLAAAAOC2SZ0muWyf0mlKPVeChoZ1LAAAAOC2W867RQUlBVq8Y7EK7AUK8Q3RLT1u0eXtLvd2afAyrlgAAADAbWaTWbN6zdLLvV5WyDshevH8F3V9l+u9XRYaAIIFAAAAqs3qY5VPto/8zH7eLgUNBMECAAAAgMcIFgAAAAA8RrAAAAAA4DGCBQAAAACPESwAAAAAeIxgAQAAAMBjBAsAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAADxGsAAAAADgMYIFAAAAAI8RLAAAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAAIDHCBYAAAAAPEawAAAAAOAxggUAAAAAjxEsAAAAAHiMYAEAAADAYwQLAAAAAB4jWAAAAADwGMECAAAAgMcIFgAAAAA8RrAAAAAA4DGCBQAAAACPESwAAAAAeIxgAQAAAMBjBAsAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAADxGsAAAAADgMYIFAAAAAI8RLAAAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAAIDHCBYAAAAAPEawAAAAAOAxggUAAAAAjxEsAAAAAHiMYAEAAADAYwQLAAAAAB7z9XYBAIBmwDDke3i1fI+slam0SCXR3VXcbpjka/N2ZQCAWkKwAADUOcv2xbIcXOH4t0/GAfmmblH+wHskM3+KAKAp4FYoAECdMuWlye/gN+XafTIOyPf4xvovCABQJwgWAIA65ZNxQCYZLreZT++v52oAAHWFYAEAqFN2W3iF2wz/ircBABoXggUAoE7Zw5NUGtamXLvhF6DiVgPcO0hOnkyp6S43mVLTpZw8DyoEUFsyCzO1PHm51hxbo1J7qbfLQT1jxBwAoM4V9Jkp69Y35ZO6RSbDrtKwtirsNkWyBle9c06e/Gc+LtOpLOW/+LCMmEjHJlNKmvxvmSsjPET5zz0gBQXU4aMAUJl3dr2j//vh/1RQWiBJig2M1T+G/kOdwjt5uTLUF4IFAKDOGdZQFVxwu1ScJ5O9WIY11O19Tbn5Mp3KkvloqvxvmesIF2Whwnw0VfZf+hkEC8Ardp7aqcfWPSbjrPFUx3OP608r/qSlY5fKbOImmeaAnzIAoP74BVQrVEiSER2h/Bcflj0+2hEuzFt2/Roq4qPPhI3oiDoqGkBVPt33qVOoKHM056g2pjL7W3NBsAAANHhGTKRTuAiY/qBzqDjr9igA9S+/JL9G29C0ECwAAI2CEROpgkfvcGorePQOQgXQAAxOGOyyPcgvSL2je9dzNfAWggUAoFEwpaTJ9tA/ndpsD/1TppQ0L1UEoMyF8Rfq8naXO7X5mHx0f9/7FeDH2KfmgsHbAIAGz2mgdny0Ch69Q7aH/lluQDcA7zCZTJo3aJ6ubHelvjv6nQL9AnVFuyvUNrStt0tDPSJYAAAaNFNqevmB2r+MuShr979lrvJfeoQB3IAXmUwmDYgfoAHxbq5PgyaHW6EAAA2aEegvIzyk3EDtswd0G+EhMgL9vVwpADRvXLEAADRsQQHKf+6BM+tU/OaKhBETeeZKRaA/i+MBgJcRLAAADV9QQIWL33H7EwA0DNwKBQAAAMBjBAsAAAAAHiNYAAAAAPAYYywAAA2CqTBTvodWy5yXptKwRJXE95N8rd4uCwDgJoIFAMDrzJmH5L/2SZmK8yRJfock+/6vlDfgbska7OXqAADu4FYoAIDXWbf/1xEqyphzjsuy9zMvVQQAqC6CBQDAu0oK5HNqj8tNvie21XMxAICaIlgAALzL5CPD7PrOXMPXVs/FAABqimABAPAuHz+VxPdxuak4YUA9FwMAqCmCBQDA6wq7TFJJZGfHvw2TWcWJQ1SSOMR7RQEAqoVZoQAA3ucXoIL+s2XOOiJTXprsIQkyAiK8XRUAoBoIFgBQi3Lz8nXgWJ6C/Fl/oWYskuKk9FIp/YS3i2n2cvILlZtn8nYZABqJBhcsnn/+ef3f//2fUlJS1L17dz333HPq06eCe2+LizVv3jy9/vrrOnr0qDp27Ki//e1vuuyyy6p1zIKCAt1555166623VFhYqBEjRuiFF15QdHR0nT5WAE3P1p/3aOHm7d4uA6g1fXt21WhvFwGgUWhQweLtt9/W7NmztWDBAvXt21fPPPOMRowYoV27dqlly5bl+s+ZM0dvvvmmFi5cqE6dOunzzz/XmDFjtGbNGvXs2dPtY86aNUtLlizRu+++q9DQUN1xxx0aO3asVq9eXa+PH0Djd965HXRVv/ZcsUCTkJNfqAOZXLEA4B6TYRiGt4so07dvX11wwQX65z//KUmy2+1KSEjQzJkzde+995brHxcXpwceeEC33367o23cuHHy9/fXm2++6dYxMzMzFRUVpcWLF2v8+PGSpJ07d+rcc8/V2rVr1a9fP7dqz8rKUmhoqDIzMxUSEuLR9wFA45SRkaGvP35LfRMDFBrk7+1yUMZeIt+j6+Vzcrvk66/ihIGyt2jn7aoahcycfH2fnKehV05SWFiYt8tBA7Nr1y7NmDFDCxcuVMeOHb1dDupIdd7jNphZoYqKirRx40YNGzbM0WY2mzVs2DCtXbvW5T6FhYWy2ZznOPf399eqVavcPubGjRtVXFzs1KdTp05q3bp1hectO3dWVpbTFwCggbGXyrb+Odl+fEV+R9fLL/lb+a+aJ9+DK7xdGQA0OQ0mWKSlpam0tLTcuIbo6GilpKS43GfEiBGaP3++9uzZI7vdri+//FLvv/++jh8/7vYxU1JSZLFYyn0SU9l5JWnevHkKDQ11fCUkJFT3IQMA6pjv8Y3yPfmTU5tJhqw7/ieVFHipKgBomhpMsKiJZ599Vh06dFCnTp1ksVh0xx13aNq0aTKb6/5h3XfffcrMzHR8HT58uM7PCQCoHp+TO1y2m0oL5XNqbz1XAwBNW4MJFpGRkfLx8VFqaqpTe2pqqmJiYlzuExUVpQ8//FC5ublKTk7Wzp07FRQUpHbt2rl9zJiYGBUVFSkjI8Pt80qS1WpVSEiI0xcAoGExLIE12gYAqL4GEywsFot69eql5cuXO9rsdruWL1+u/v37V7qvzWZTfHy8SkpK9N577+mqq65y+5i9evWSn5+fU59du3bp0KFDVZ4XANCwlSQMkGEq/6euNKSV7GFtvVARADRdDWq62dmzZ2vq1Knq3bu3+vTpo2eeeUa5ubmaNm2aJOn6669XfHy85s2bJ0n6/vvvdfToUfXo0UNHjx7VI488IrvdrrvvvtvtY4aGhmr69OmaPXu2wsPDFRISopkzZ6p///5uzwgFAGiY7MHxKuw5XdZt/5WpOEeSVBraWgW9bvNyZQDQ9DSoYDFx4kSdPHlSDz30kFJSUtSjRw8tW7bMMfj60KFDTuMnCgoKNGfOHO3fv19BQUEaNWqUFi1a5DQQu6pjStLTTz8ts9mscePGOS2QBwBo/Eri+6ok5nyZMw5Ifv6yhzDZBgDUhQa1jkVjxjoWAFjHAk0N61igMqxj0Tw0ynUsAAAAADReBAsAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAADxGsAAAAADgMYIFAAAAAI81qAXyAAAAPJWXl6eioiJvl9Hk5eTkOP6bkZHh3WKaOIvFooCAAG+XUSWCBQAAaDLy8vL08bKPlZWf5e1Smry83Dx169VNKzeu1MadG71dTpMW4h+iKy+7ssGHC4IFAKDZMBVmSaWFMgKivF0K6khRUZGy8rNkibXIGmD1djlNWrCCFX1etLfLaPIK8wqVdTxLRUVFBAsAALzNVJgp65bX5ZO6TSYZsgfFqrDbtSqN7OTt0lBHrAFWBQQ17DdhgLuK1Dhu7WPwNgCgybP98IJ8U7fKJEOSZM45Ltv6f8iUl+7lygCg6SBYAACaNHNGsnxO7yvXbiotkt/hVV6oCACaJoIFAKBJMxWcrmRbZj1WAgBNG8ECANCk2Vu0lWF2PaSwNDypnqsBgKaLYAEAaNIMa6iKky4r114a1kYlcX28UBEANE3MCgUAaPKKOl2t0pAE+R1ZLVNxgUqiu6m4zcWSj5+3SwOAJoNgAQBoFkrjeqk0rpe3ywCavFOFp2QxWxTkF+TtUlDPCBYAAABwi292rnxz81UQE1lu24l9m/TxqRU6aD4tk0zqENJBY9qMUbBfsBcqhTcwxgIA0OiYsw7Lb9/n8j20UirO83Y5QLPgm52r/rc+qkHTHpAt5aTTtuJD+zXy5if02OMbFJRXKkOGdmft1pt73/RStfAGrlgAABoVy7bFshz82vFvY9tilUZ2kj0kXsUJg2QExXixOqDp8s3Nl/VUhgKPpGrQtDla9epjKoiJki3lpAbeNFeRJ4pkGH4KLLArJ8BHknQ076iO5B5Rq8BWXq4e9YErFgCARsMndatTqJAkk71Yvie2ybJ3mQK+fUQ+KT96pzigiSuIidSqVx9TbqtoR7ho8eNODZo2R5HHM3U4yk/T7m2r1HDnSRGyirO8VDHqG8ECANBo+B7fWOl2k71E1m3/kQx7PVUENC8FMVFO4eKi6+5V4JFUnYprcSZURFic+vuYfJQQmOClalHfCBYAgCbFXHBa5qyj3i4DaLIKYqK0cd4sp7Yt8+6UpXX5BScHRg9k8HYzwhgLAECjURJ3gfwOr66yn+Fnq4dqgObJlnJSve572qmtzwPP6ff/fljfhu3XrsxdsvhY1DOip7q26OqlKuENBAsAqGW5+YXeLqHpCkhSaKuLFXRkRYVdCsPaK8MeJOXk12NhTRPPZfyWLeWkBk2bo8AjqcptFa2N82ap131PK/BIqi6+aa78Xn1MgzoO8naZ8BKCBQDUEovFIr/AUG0/kSmJKVDrSm7JYKVkBmpIq2LF2I8rNu9n+RrFkqRT1lb6PmSC8pP5/tcWv8BQWSyWqjuiybOlpDmFirJZoVa9+pij/cxsUY+7XOcCTR/BAgBqSUBAgC4ZeZWKioq8XUqTtm/fPi1c/IGG/u5pBSclKbcgUz4pm2UERMjcspv6e7vAJsZisSggIMDbZaABKAn0V2F4mCQ5QoUkp3BRGB6mkkB/L1YJbyJYAEAtCggI4E1YHQsKClKoX4laH3pPYYdOSy3PlXr+TgoI93ZpQJNWEhyotQsecrny9plw8bhKAv1VEhzopQrhbQQLAECj4peVrNf77lL4lp9+bVz3L+nGz6QWbbxWF9AclAQHVhgcuP0JTDcLAGhUorb8U+GWEufG7GPSinneKQgAIIlgAQBoZAKPr3O9Ye+X9VsIAMAJwQIA0KjY/Sq4f9saUr+FAACcECwAAI1KZtvRrjf0vLZ+CwEAOCFYAAAalfRuM/RVapgMmc40mMxSj2ulgbO8WxgANHPMCgUAaFQMH4se/SlR7W56RO1C7FLUOVJYa2+XBQDNHsECANCo+OSnqVNwnkqsYVKH871dDgDgF9wKBQBoHIoLpPdvUdJHo/XSBXvU/sPLpe/+z9tVAQB+QbAAADQOXz4kbX1LJqNUkmQuyZO+fkza+o6XCwMASAQLAEBjUFIkbX7T9baNr9VrKQAA1wgWAICGryRfKs51vS33ZP3WAgBwiWABAGj4bKFSbA/X29oOrtdSAACuESwAAI3D8MckX5tzW0i8NOhPXikHAOCMYAEAaBzaXijd8p1On3ONVqeFKK3bzdLN30qhrbxdGQBABAsAQGMS1VEnev1Z921tq/SuN0lBUd6uCADwC4IFAAAAAI8RLAAAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAAIDHCBYAAAAAPEawAAAAAOAxggUAAAAAjxEsAAAAAHiMYAEAAADAYwQLAAAAAB7z9XYBAABIklK2S9vekUoKpY4jpXZDvF0RAKAaCBYAAO9bv1Baepck48y/v18g9bpBGv2sN6sCAFQDt0IBALwrN136/AE5QkWZja9JyWu9UREAoAYIFgAA79q/QiotdL1t97L6rQUAUGMECwCAd1kCK9kWVH91AAA8QrAAAHhX0iVSUHT5drOv1G18/dcDAKgRggUAwLt8LdLE/0jBcb+2WYKlq/8lhbf1Xl0AgGphVigAgPclXCD9aZt0cKVUUiC1uVCychsUADQmBAsAQMPg4yslXeztKgAANcStUAAAAAA8RrAAAAAA4DGCBQAAAACPESwAAAAAeIxgAQAAAMBj1QoWubm5ysnJqbRPTk6OcnNzPSoKAAAATcfR3KPak7VHhaWF3i4FdcjtYLFr1y61aNFCf/3rXyvtN2/ePIWHh2vfvn0eFwcAAIDG63Thab3w8wv6185/6fU9r+tvW/+m70987+2yUEfcDhYLFixQVFSUHn744Ur7Pfjgg4qKitKCBQs8Lg4AAACN19sH3taxvGOOfxfZi/TJ4U90KOeQF6tCXXE7WHz11VcaP368rFZrpf1sNpsmTJigL774wuPiAAAA0DidyD+hI7lHXG7blL6pnqtBfXA7WBw4cEBdunRxq++5557LrVAAgMbBXirlnZLsdm9XAjQpBaUFNdqGxsvX3Y52u11ms3s5xGw2y84vaABAQ7fmn9Kaf0g5qVJIK+miO6XeN3q7KqBJiAuIU6BvoHJLyk/q0yGkgxcqQl1z+4pFbGysduzY4VbfHTt2KDY2tsZFAQBQ59YvlL544EyokKSsI9Kns6Rt//NuXUAT4Wv21eUJl8v8m7ebbYPbqnt4dy9Vhbrk9hWLCy+8UIsXL9ajjz6qoKCgCvvl5ORo8eLFGjlyZK0UCACoOwUFBUpOTvZ2GdVSVq+ndbf97hlZXLTnr3hShyzdPDp2ZRITE2Wz2ers+EBDcl74eYr2j9am9E3KK8lTUnCSuoV3k4/Jx9uloQ64HSzuuOMOLVq0SGPGjNHbb7+t8PDwcn1Onz6ta665RmlpaZo5c2atFgoAqH3JycmaMWOGt8uokccee8yj/VdcfEQylW8vTNldp9+ThQsXqmPHjnV2fKChifaP1shWfODcHLgdLHr37q2HHnpIc+fOVdu2bTV27Fh1795dISEhys7O1ubNm/Xhhx8qKytLc+fO1fnnn1+XdQMAakFiYqIWLlzo7TK8ouiLG+Wfvr1cu29iHy2c+mydnTcxMbHOjg0A3uR2sJCkhx9+WAkJCXrggQf0+uuvS5JMJpMMw5AkxcTE6Omnn9a0adNqv1IAQK2z2WzN99Nzn0elxddI9pKz2qwKGvWoOiY00+8JAHigWsFCkm688UZdd911Wr16tbZv366srCyFhISoa9euGjhwoPz8/OqiTgAAalf7S6QblkhrnpPS9kgtz5UG/UmK6+ntygCgUap2sJAkPz8/DRkyREOGDKnlcgAAqEet+535AgB4zO3pZgEAAACgIm5fsWjXrl2F20wmk/z9/dW2bVuNGTNGU6dOlY8P04gBAAAAzYXbVyzsdrsMw3D5ZbfbdeLECS1ZskQzZszQ8OHDVVxcXJd1AwAAAGhA3L5icfDgwSr7pKWl6amnntLf/vY3PfPMM7rrrrs8qQ0AAABAI1GrYywiIyM1b948XXbZZXrrrbdq89AAAAAAGrA6Gbx9ySWXaPfu3TXa9/nnn1ebNm1ks9nUt29frV+/vtL+zzzzjDp27Ch/f38lJCRo1qxZKigocGxv06aNTCZTua/bb7/d0WfIkCHltt966601qh8AAABojmo03WxV/P39VVJSUnXH33j77bc1e/ZsLViwQH379tUzzzyjESNGaNeuXWrZsmW5/osXL9a9996rV155RQMGDNDu3bt1ww03yGQyaf78+ZKkH374QaWlpY59tm/frksvvVQTJkxwOtaMGTP06KOPOv4dEBBQ7foBAAAaixP5J7TmxBqdLDipKFuUBrQcoJb+5d9vAe6qk2CxefNmtWrVqtr7zZ8/XzNmzHCs3L1gwQItWbJEr7zyiu69995y/desWaOBAwdqypQpks5cnZg8ebK+//57R5+oqCinfZ544gklJSVp8ODBTu0BAQGKiYmpds0AAACNzdHco3p598sqshdJkpJzkrX11FbdeM6NahVY/fdwgFQHt0J9+eWXeuONNzRy5Mhq7VdUVKSNGzdq2LBhvxZnNmvYsGFau3aty30GDBigjRs3Om6X2r9/v5YuXapRo0ZVeI4333xTN954o0wmk9O2//znP4qMjFTXrl113333KS8vr9J6CwsLlZWV5fQFAADQGHx17CtHqChTZC/S8mPLvVQRmgK3r1jceOONlW7Py8vTzz//rO3btysmJkb3339/tQpJS0tTaWmpoqOjndqjo6O1c+dOl/tMmTJFaWlpGjRokAzDUElJiW699dYKz/3hhx8qIyNDN9xwQ7njJCYmKi4uTlu3btU999yjXbt26f3336+w3nnz5mnu3LnVeowAAAANQXJOcrXaAXe4HSxee+21KvsEBQVp0qRJeuKJJ+rltqJvvvlGf/3rX/XCCy+ob9++2rt3r/74xz/qL3/5ix588MFy/V9++WWNHDlScXFxTu0333yz4/+7deum2NhYXXLJJdq3b5+SkpJcnvu+++7T7NmzHf/OyspSQkJCLT0yAACAuhPsF6z0wnSX7UBNuR0sDhw4UOl2f39/RUZGymyu2d1VkZGR8vHxUWpqqlN7ampqhSHlwQcf1HXXXaebbrpJ0plQkJubq5tvvlkPPPCAUy3Jycn66quvKr0KUaZv376SpL1791YYLKxWq6xWq1uPDQAAoCHp17Kflhxe4rIdqCm3g0ViYmJd1iGLxaJevXpp+fLluvrqqyWdWe17+fLluuOOO1zuk5eXVy7I+Pj4SJIMw3Bqf/XVV9WyZUtdfvnlVdby448/SpJiY2Or+SgAAAAavn5R/ZRXkqc1qWtUaC+U1WzVgOgB6hdFsEDN1fqsUKmpqXr99df12muvaceOHdXad/bs2Zo6dap69+6tPn366JlnnlFubq5jlqjrr79e8fHxmjdvniRp9OjRmj9/vnr27Om4FerBBx/U6NGjHQFDOhNQXn31VU2dOlW+vs4Ped++fVq8eLFGjRqliIgIbd26VbNmzdJFF12k8847z8PvBgAAQMNjMpl0SdwlujDmQmUVZSnEEiKL2eLtstDI1UqwsNvt+vTTT/Xyyy/rs88+U0lJiYKDq3+P3sSJE3Xy5Ek99NBDSklJUY8ePbRs2TLHgO5Dhw45XaGYM2eOTCaT5syZo6NHjyoqKkqjR4/W448/7nTcr776SocOHXI5AN1iseirr75yhJiEhASNGzdOc+bMqXb9AAAAjYnFbFGkLdLbZaCJMBm/vWeoGnbt2qVXXnlFixYtUmpqqlq0aKErr7xS48aN06WXXiqLpfkk36ysLIWGhiozM1MhISHeLgcAgGYpIyND7yx5R8FJwQoIYrFbNH55OXnK3petay6/RmFhYfV+/uq8x632FYvc3Fy9/fbbeuWVV7R27Vr5+vpq4MCBSk1N1UsvvaSxY8fWuHAAAAAAjZPbUzitXr1a06dPV2xsrG666Sbl5+fr6aef1rFjx/TSSy+VGywNAAAAoPlw+4rFhRdeqOjoaN18882aOnWqunXr5tiWkZFRF7UBAAAAaCSqtehEfn6+MjMzlZWVVVf1AAAAAGiE3A4WO3bs0E033aRPPvlEF110kdq3b6+//OUvSk5m6XcAAACguXM7WHTq1ElPPvmkjhw5onfffVcdO3bUo48+qqSkJE2YMEEmk4lxFgAAAEAzVa1boSTJ19dXY8eO1ZIlS3To0CE9+uijys7OlmEYuu666zRmzBi9+eabyszMrIt6AQAAADRA1Q4WZ4uNjdX999+vPXv2aMWKFRo/fry++OILXX/99Y5F7QAAAAA0fR4Fi7MNHjxYb7zxho4fP64XXnhB5513Xm0dGgAAAEADV2vBokxISIhuvfVWrV+/vrYPDQAAAKCBqvVgAQAAAKD5IVgAAAAA8BjBAgAAAIDHCBYAAAAAPEawAAAAAOCxagWLnTt3aufOnR73AQAAANC0uB0sNmzYoC5duujzzz+vtN/nn3+url27auvWrR4XBwAAAKBxcDtYLFy4UG3bttXMmTMr7Tdz5ky1a9dOCxYs8Lg4AAAAAI2D28Him2++0bhx42Q2V76L2WzWuHHjtGLFCo+LAwAAANA4uB0sDh8+rHPOOcetvu3bt1dycnKNiwIAAADQuLgdLMxms0pKStzqW1JSUuWVDQAAAABNh9vv/hMSEvTjjz+61ffHH39Uq1ataloTAAAAgEbG7WBx8cUX67///a9SU1Mr7Zeamqr//ve/uuSSSzwuDgAAAEDj4HawmDVrlvLz8zVixAjt2rXLZZ/du3dr5MiRKigo0J/+9KfaqhEAAABAA+frbscOHTrohRde0C233KIuXbpo0KBB6tGjh0JCQpSdna3Nmzdr9erVMgxDCxcuVIcOHeqybgAAAAANiNvBQpKmT5+uxMRE3XPPPfruu+/03XffOW3v2bOn/va3v2nYsGG1WiQAAACAhq1awUKShg0bpo0bN+rgwYPavn27srKyFBISoq5du6pNmzZ1UCIAAACAhq7awaJMmzZtCBIAAAAAJNUwWKSnp2vp0qXatm2bMjMzFRoaqm7dumnkyJGKjIys7RoBAAAANHDVChaGYWju3Ll68sknlZ+fL8MwHNtMJpNsNpvuvvtuPfTQQzKZTLVeLAAAAICGqVrBYtq0aXrjjTfUunVr/e53v1OvXr0UEhKirKwsbdy4UW+++aYeffRRHThwQK+99lodlQwAAACgoXE7WHz00Ud64403NHXqVC1YsEBWq9Vp+5gxY/Tggw/qtttu0+uvv66xY8fqyiuvrPWCAQAAADQ8bi+Q99JLL6lr1656+eWXy4WKMlarVf/+97/VtWtXvfjii7VWJAAAAICGze1gsWHDBk2ZMkVmc+W7mM1mTZ48WRs2bPC4OAAAAACNg9vBIiMjQzExMW71jYmJUWZmZo2LAgAAANC4uB0sIiMjdeDAAbf6Hjx4UBERETUuCgAAAEDj4nawGDBggBYtWqSCgoJK+xUUFOiNN97QgAEDPC4OAAAAQOPgdrCYOXOmDh48qKuvvlrp6eku+5w6dUpjxoxRcnKyZs6cWWtFAgAAAGjY3J5u9qKLLtK9996rJ554Qu3atdPVV1+tnj17KjQ0VJmZmdq0aZM++ugjZWdn66677tJFF11Ul3UDAAAAaECqtUDeX//6VyUlJWnOnDlatGiRFi1aJJPJ5FiBOzo6Wk8++aRmzJhRJ8UCAAAAaJiqFSwkafr06br++uu1evVqbd++XVlZWQoJCVHXrl01YMAAWSyWuqgTAAAAQANW7WAhSX5+fhoyZIiGDBlSy+UAAAAAaIzcHrxdHatXr9Yll1xSF4cGAAAA0ABVO1ikp6dr/fr12rt3b7lt69at0/Dhw3XRRRfpu+++q5UCAQAAADR8bgeL0tJS3XrrrYqOjlb//v3VsWNHDRgwQCdOnFBWVpamTJmigQMHasWKFZoyZYq2bdtWl3UDAAAAaEDcHmPx3HPP6aWXXlKrVq3Ur18/7d27V+vWrdPtt9+uI0eOaP369bruuuv04IMPKikpqS5rBgAAANDAuB0sFi1apG7dumnt2rUKCAiQJN1+++3617/+pYiICK1atUr9+/evs0IBAAAANFxu3wq1e/duXX/99Y5QIUm33XabJOmee+4hVAAAAADNmNvBIjc3VzExMU5tZf/u1q1b7VYFAAAAoFGp1qxQJpPJ5b/9/PxqryIAAAAAjU61FshbunSpUlJSHP/Oy8uTyWTSu+++qx9//NGpr8lk0qxZs2qlSAAAAAANm8kwDMOdjmZz9Za8MJlMKi0trVFRjVFWVpZCQ0OVmZmpkJAQb5cDAECzlJGRoXeWvKPgpGAFBAVUvQPQwOXl5Cl7X7auufwahYWF1fv5q/Me1+0rFitWrPC4MAAAAABNk9vBYvDgwXVZBwAAAIBGrHr3NwEAAACACwQLAAAAAB4jWAAAAADwGMECAAAAgMcIFgAAAAA8RrAAAAAA4DGCBQAAAACPESwAAAAAeIxgAQAAAMBjBAsAAAAAHvP1dgEAAABoHIrsRdqUtkn7svYpwDdAvSN7KyEowdtloYEgWAAAAKBKRfYivbzrZR3NO+po25S+SVcnXq1ekb28WBkaCm6FAgAAQJU2p212ChWSZMjQsiPLVGwv9lJVaEgIFgAAAKjSvux9LtvzS/N1PO94PVeDhohgAQAAgCoF+gbWaBuaD4IFAAAAqtQ7srdMMpVrbxvcVhG2CC9UhIaGYAEAAIAqxQfGa2ybsU5XJ9oFt9M1ba/xYlVoSJgVCgAAAG7pGdFT3Vp0U0p+igJ8AxRuDfd2SWhACBYAAABwm6/ZV60CW3m7DDRA3AoFAAAAwGNcsQAAAGhGDmQf0I6MHTLLrPPCz1N8YLy3S0ITQbAAAABoJpYeXqo1J9Y4/r36xGpdGn+pBscM9mJVaCq4FQoAAKAZOJZ3zClUlFl+dLkyijLqvyA0OQQLAACAZmBX5i6X7XbZtSdzTz1Xg6aIYAEAANAMWM3WCrdZfCz1WAmaKoIFAABAM9AtvJt8TeWH1/r7+Ovc0HO9UBGaGoIFAABAMxDsF6xr2l0jfx9/R1uQb5CmJE3higVqBbNCAQAANBOdwzqrw3kdtD97v8wms9oFt5OPycfbZaGJIFgAAAA0I35mP3UM7ejtMtAEcSsUAAAAAI8RLAAAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAAIDHGlyweP7559WmTRvZbDb17dtX69evr7T/M888o44dO8rf318JCQmaNWuWCgoKHNsfeeQRmUwmp69OnTo5HaOgoEC33367IiIiFBQUpHHjxik1NbVOHh8AAADQFDWoYPH2229r9uzZevjhh7Vp0yZ1795dI0aM0IkTJ1z2X7x4se699149/PDD+vnnn/Xyyy/r7bff1v333+/Ur0uXLjp+/Ljja9WqVU7bZ82apU8++UTvvvuuvv32Wx07dkxjx46ts8cJAAAANDUNaoG8+fPna8aMGZo2bZokacGCBVqyZIleeeUV3XvvveX6r1mzRgMHDtSUKVMkSW3atNHkyZP1/fffO/Xz9fVVTEyMy3NmZmbq5Zdf1uLFizV06FBJ0quvvqpzzz1X69atU79+/WrzIQIAAABNUoO5YlFUVKSNGzdq2LBhjjaz2axhw4Zp7dq1LvcZMGCANm7c6Lhdav/+/Vq6dKlGjRrl1G/Pnj2Ki4tTu3btdO211+rQoUOObRs3blRxcbHTeTt16qTWrVtXeF5JKiwsVFZWltMXAAAA0Fw1mCsWaWlpKi0tVXR0tFN7dHS0du7c6XKfKVOmKC0tTYMGDZJhGCopKdGtt97qdCtU37599dprr6ljx446fvy45s6dqwsvvFDbt29XcHCwUlJSZLFYFBYWVu68KSkpFdY7b948zZ07t+YPGAAAAGhCGswVi5r45ptv9Ne//lUvvPCCNm3apPfff19LlizRX/7yF0efkSNHasKECTrvvPM0YsQILV26VBkZGXrnnXc8Ovd9992nzMxMx9fhw4c9fTgAAABAo9VgrlhERkbKx8en3GxMqampFY6PePDBB3XdddfppptukiR169ZNubm5uvnmm/XAAw/IbC6fm8LCwnTOOedo7969kqSYmBgVFRUpIyPD6apFZeeVJKvVKqvVWt2HCQAAADRJDeaKhcViUa9evbR8+XJHm91u1/Lly9W/f3+X++Tl5ZULDz4+PpIkwzBc7pOTk6N9+/YpNjZWktSrVy/5+fk5nXfXrl06dOhQhecFAAAA4KzBXLGQpNmzZ2vq1Knq3bu3+vTpo2eeeUa5ubmOWaKuv/56xcfHa968eZKk0aNHa/78+erZs6f69u2rvXv36sEHH9To0aMdAePPf/6zRo8ercTERB07dkwPP/ywfHx8NHnyZElSaGiopk+frtmzZys8PFwhISGaOXOm+vfvz4xQAAAAgJsaVLCYOHGiTp48qYceekgpKSnq0aOHli1b5hjQfejQIacrFHPmzJHJZNKcOXN09OhRRUVFafTo0Xr88ccdfY4cOaLJkycrPT1dUVFRGjRokNatW6eoqChHn6efflpms1njxo1TYWGhRowYoRdeeKH+HjgAAADQyJmMiu4ZQrVkZWUpNDRUmZmZCgkJ8XY5AAA0SxkZGXpnyTsKTgpWQFCAt8sBPJaXk6fsfdm65vJrys1iWh+q8x63wYyxAAAAANB4NahboQAAANA8ZRRl6IeTP+hU4SnFBsSqd2RvBfhy1akxIVgAAADAq47mHtUru19Rob1QkrTt9DZ9f/J7zeg4Q2GWMO8WB7dxKxQAAAC86rMjnzlCRZnMokx9e/xbL1WEmiBYAAAAwGtK7CU6mHPQ5ba9WXvrtxh4hGABAAAArzGbzLKarS63+fv613M18ATBAgAAAF5jNpl1fuT5Lrf1iuxVz9XAEwQLAAAAeNXw+OHq1qKbTDJJknxMPhrYcqD6RPbxcmWoDmaFAgAAgFf5mf00sd1EjSgaoVOFpxRti1agX6C3y0I1ESwAAADQIIRZwphethHjVigAAAAAHiNYAAAAAPAYt0IBAIAmJS83TwVHCmQLsHm7FMBjBXkFsufavV2GWwgWAACgSdnz8x5t27jN22UAtaZbr27eLsEtBAsAANCkdDi3gzoO7MgVCzQJBXkFsp/iigUAAEC9CwgMUHCrYAUEBXi7FMBjeTl5yi7M9nYZbmHwNgAAAACPccUCAAAAte5I7hFtStuk/NJ8tQ9pr+7h3eVr5q1nU8ZPFwAAALVqQ9oGfZT8kQwZkqRtp7dpc/pm3dDhBsJFE8atUAAAAKg1RfYiLTuyzBEqyhzMOaitp7d6qSrUB4IFAAAAas3R3KMqKC1wuW1v1t56rgb1iWABAACAWuPv61/htgAfZupqyggWAAAAqDUx/jFqFdCqXLtJJp0feb4XKkJ9IVgAAACgVk1KmqRWgb+GC38ff41tM1ZxAXFerAp1jWH5AAAAqFVhljDd2ulWpeanKr80X/EB8fIz+3m7LNQxggUAAADqRLR/tLdLQD3iVigAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAqBG7YVexvdjbZaCBYPA2AAAAqqXUKNXyY8v1w8kfHLM+DY8frqSQJG+XBi/iigUAAACq5bPDn+m7lO+UX5ovSTqad1SL9i5SSl6KlyuDNxEsAAAA4LaC0gJtSNtQrr3EKNHaE2u9UBEaCoIFAAAA3JZVlKUSo8TltlNFp+q5GjQkBAsAAAC4rYW1hfx9/F1uiwuIq+dq0JAQLAAAAOA2P7OfLoq5qFx7gG+A+rfs74WK0FAwKxQAAACq5cKYCxXsF6z1J9cruzhbbYLbaEjMEIVZwrxdGryIYAEAAIBq6xHRQz0ieni7DDQgBAsAAABUS15JnjalbVJqQapa2lrq/MjzFegb6O2y4GUECwAAALjtVOEpLdy1UNnF2Y621amrdVPHmxRpi/RiZfA2Bm8DAADAbV8d/copVEhSTkmOvjz6pZcqQkNBsAAAAIDb9mTtqVY7mg+CBQAAANxm9bFWqx3NB8ECAAAAbjs/4vxqtaP5IFgAAADAbYNjB6tHeA+ZZJIkmWTSeeHn6eLYi71cGbyNWaEAAADgNh+Tj8a3Ha9L4i7RyYKTirRFKtwa7u2y0AAQLAAAAFBtLawt1MLawttloAHhVigAAAAAHiNYAAAAAPAYwQIAAACAxwgWAAAAADxGsAAAAADgMYIFAAAAAI8RLAAAAAB4jGABAAAAwGMECwAAAAAeI1gAAAAA8BjBAgAAoJnLKsrSvqx9Ol142tuloBHz9XYBAAAA8A67Ydenhz7VhrQNsssuk0zqFt5NYxPHytfM20RUD1csAAAAmqk1qWu0Pm297LJLkgwZ2npqq5YfW+7lytAYESwAAACaqY3pG6vVDlSGYAEAANBM5Zfku2wvKCmo50rQFBAsAAAAmqn2Ie1dtncI7VDPlaApIFgAAAA0U5fEXaJgv2CntgCfAA2PH+6litCYMdwfAAD8f3v3HR5Fub9//N6EbHpCS0JPICC9hiK9RRAQD0WliCQgRQQLeEBQukcjipAjihQR/KoIUgQURUHBAkgTxIKhSRHpJQkJEJI8vz/4ZQ/LJhDYhA3h/bquvTTPPDPzmWV2du+dZ2ZxlyrkWUhPVXlKW09t1bELx1TUs6jqBdVzCBtAdhAsAAAA7mI+BXzUrFgzV5eBfIChUAAAAACcRrAAAAAA4DSCBQAAAACnESwAAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGkECwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAyKfSTbpOXDihhJQEV5eCu0ABVxcAAACAnBcXH6fPD32usylnJUnlA8qra1hX+Xv4u7gy5FecsQAAAMhnTl08pfn75ttChSTtTdirj/d97MKqkN8RLAAAAPKZbae2Kc2kObQfSjqko8lHXVAR7gYECwAAgHzmfOr5rKddznoa4AyCBQAAQD4T5heWabuHxUOlfEvd3mJw18hzweLtt99WWFiYvLy81KBBA23evPm6/WNjY1WxYkV5e3urdOnSGjp0qC5evGibHhMTo3r16snf31/BwcHq1KmT4uLi7JbRokULWSwWu8cTTzyRK9sHAACQ22oUrpFpgGhRooW8C3i7oCLcDfLUXaEWLlyoYcOGacaMGWrQoIFiY2PVtm1bxcXFKTg42KH//PnzNXLkSL333ntq1KiRdu/erejoaFksFk2ZMkWS9N1332nw4MGqV6+eUlNT9cILL6hNmzb6448/5Ovra1tW//79NXHiRNvfPj4+ub/BAAAAucDDzUN9K/TVllNbtDt+t7zcvRRRNEL3BN7j6tKQj+WpYDFlyhT1799fffr0kSTNmDFDK1eu1HvvvaeRI0c69N+wYYMaN26snj17SpLCwsLUo0cPbdq0ydZn1apVdvPMmzdPwcHB2rZtm5o1a2Zr9/HxUbFixXJjswAAAG47q7tVjUMaq3FIY1eXgrtEnhkKlZKSom3btikyMtLW5ubmpsjISG3cuDHTeRo1aqRt27bZhkvt379fX3zxhdq3b5/leuLj4yVJhQsXtmv/6KOPVLRoUVWrVk2jRo1ScnLydeu9dOmSEhIS7B4AAADA3SrPnLE4deqU0tLSFBISYtceEhKiP//8M9N5evbsqVOnTqlJkyYyxig1NVVPPPGEXnjhhUz7p6en69lnn1Xjxo1VrVo1u+WEhoaqRIkS2rlzp55//nnFxcVp6dKlWdYbExOjCRMm3MKWAgAAAPlPngkWt2LdunV65ZVXNH36dDVo0EB79+7VM888o5deekljxoxx6D948GD99ttv+vHHH+3aBwwYYPv/6tWrq3jx4mrdurX27dun8PDwTNc9atQoDRs2zPZ3QkKCSpcunUNbBgAAANxZ8kywKFq0qNzd3XX8+HG79uPHj2d57cOYMWP02GOPqV+/fpKuhIKkpCQNGDBAL774otzc/jfSa8iQIfr888/1/fffq1Sp699mrUGDBpKkvXv3ZhksPD095enpme3tAwAAAPKzPHONhdVqVUREhL755htbW3p6ur755hs1bNgw03mSk5PtwoMkubu7S5KMMbb/DhkyRJ9++qm+/fZblS1b9oa17NixQ5JUvHjxW9kUAAAA4K6TZ85YSNKwYcMUFRWlunXrqn79+oqNjVVSUpLtLlG9e/dWyZIlFRMTI0nq2LGjpkyZotq1a9uGQo0ZM0YdO3a0BYzBgwdr/vz5Wr58ufz9/XXs2DFJUmBgoLy9vbVv3z7Nnz9f7du3V5EiRbRz504NHTpUzZo1U40aNVzzRAAAAAB3mDwVLLp166aTJ09q7NixOnbsmGrVqqVVq1bZLug+dOiQ3RmK0aNHy2KxaPTo0Tpy5IiCgoLUsWNHvfzyy7Y+77zzjqQrP4J3tblz5yo6OlpWq1Vr1qyxhZjSpUura9euGj16dO5vMAAAQD53Of2yfjz+o34785vSla6qBauqabGm8nRnSHl+YzEZY4bglISEBAUGBio+Pl4BAQGuLgcAgLvSuXPn9MnKT+Qf7i8fP37sNi/4vz3/p90Ju+3ayviWUb+K/eRmyTOj8vOs5PPJStyXqEc6PKKCBQve9vXfzGdc/jUBAACQKw6dP+QQKiTpUNIh7UnY44KKkJsIFgAAAMgV/yT/c0vTcGciWAAAACBXFPIslOW0wp6Fb2MluB0IFgAAAMgVFQIqKMQrxKG9oLWgqhas6oKKkJsIFgAAAMgVbhY3Rd8TrWqFqslNbrLIokqBlfT4PY+rgFueujkpcgD/ogAAAMg1/h7+6l6uu1LTUyWJQJGP8S8LAACAXEegyP8YCgUAAADAaQQLAAAAAE4jWAAAAABwGsECAAAAgNO4igYAAOQ7l5IvuboEIEfcSfsywQIAAOQbVqtVAd4BSjiaoBSluLqcfC05KVl7du1RhcoV5OPr4+py8rUA7wBZrVZXl3FDBAsAAJBv+Pj46MH7H1RKCqEit+3bt09LP1yqvr36Kjw83NXl5GtWq1U+Pnk/vBEsAABAvuLj43NHfAi70/n5+dn+W7BgQdcWgzyBi7cBAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGkECwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4jWAAAAABwGsECAAAAgNMIFgAAAACcRrAAAAAA4LQCri4AAAAAd48dJ3bo5xM/K9gnWJFlIuVVwMvVJSGHECwAAACQ61LTUzX8u+Fac2iNrW2q91TNajNL4QXDXVgZcgpDoQAAAJDrlu5ZahcqJOnEhRMat2GciypCTiNYAAAAINd9ffDrTNt/OfmLjicdv83VIDcQLAAAAOBSFovF1SUgBxAsAAAAkOvahrXNtL12cG0F+wTf5mqQG7h4GwAAADlqf/x+rdy/UhdSL6h5qeZqULyBupTvok1HN+mrA1/Z+oX4hGh8o/GuKxQ5imABAACAHPPpnk81fuN4pZt0SdIHf3ygLhW6aEKjCZrcfLKiqkRp+4ntCvYJVqsyrWR1t7q4YuQUggUAAAByRGJKomI2x9hCRYale5aqfdn2alC8gaoHVVf1oOouqhC5iWssAAAAkCM2H92sC6kXMp227vC621oLbj+CBQAAAHLE9X5F27uA922sBK5AsAAAAECOqF+8fqZ3eHKzuKlDuQ4uqAi3E8ECAAAAOcLDzUOxLWIV7P2/cOFdwFtj7x2r8ILhLqwMtwMXbwMAACDHVA+qrlUPrdKmo5t0MfWi6hevrwBrgKvLwm1AsAAAAECO8nDzUJOSTVxdBm4zhkIBAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4jWAAAAABwGsECAAAAgNMIFgAAAACcRrAAAAAA4DSCBQAAAACnESwAAAAAOK2AqwsAAABA3vT939/rwz8+1NGko6oRVEOPV39c5QLLubos5FEECwAAADhYvne5Rq8fbfv7QMIBrT28VvPbz1dYYJjrCkOexVAoAAAA2Ek36Xp7x9sO7YkpiZr3+7zbXxDuCAQLAAAA2Dl94bSOJh3NdNrvp3+/zdXgTkGwAAAAgJ1Az0D5efhlOq2Eb4nbXA3uFAQLAAAA2LG6W9W9UneHdjeLm3pV6eWCinAn4OJtAAAAOBhSa4jcLe5aELdA8ZfiVTawrJ6q/ZTqFavn6tKQRxEsAAAA4MDdzV1Dag/REzWf0IXUC/K3+ru6JORxDIUCAABAlgq4FSBUIFsIFgAAAACcRrAAAAAA4DSCBQAAAACnESwAAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGl5Lli8/fbbCgsLk5eXlxo0aKDNmzdft39sbKwqVqwob29vlS5dWkOHDtXFixdvapkXL17U4MGDVaRIEfn5+alr1646fvx4jm8bAAAAkF/lqWCxcOFCDRs2TOPGjdPPP/+smjVrqm3btjpx4kSm/efPn6+RI0dq3Lhx2rVrl+bMmaOFCxfqhRdeuKllDh06VJ999pkWLVqk7777Tv/884+6dOmS69sLAAAA5Bd5KlhMmTJF/fv3V58+fVSlShXNmDFDPj4+eu+99zLtv2HDBjVu3Fg9e/ZUWFiY2rRpox49etidkbjRMuPj4zVnzhxNmTJFrVq1UkREhObOnasNGzbop59+yrLWS5cuKSEhwe4BAAAA3K3yTLBISUnRtm3bFBkZaWtzc3NTZGSkNm7cmOk8jRo10rZt22xBYv/+/friiy/Uvn37bC9z27Ztunz5sl2fSpUqqUyZMlmuV5JiYmIUGBhoe5QuXfrWNx4AAAC4w+WZYHHq1CmlpaUpJCTErj0kJETHjh3LdJ6ePXtq4sSJatKkiTw8PBQeHq4WLVrYhkJlZ5nHjh2T1WpVwYIFs71eSRo1apTi4+Ntj8OHD9/sJgMAAAD5Rp4JFrdi3bp1euWVVzR9+nT9/PPPWrp0qVauXKmXXnop19ft6empgIAAuwcAAABwtyrg6gIyFC1aVO7u7g53Yzp+/LiKFSuW6TxjxozRY489pn79+kmSqlevrqSkJA0YMEAvvvhitpZZrFgxpaSk6Ny5c3ZnLa63XgAAAAD28swZC6vVqoiICH3zzTe2tvT0dH3zzTdq2LBhpvMkJyfLzc1+E9zd3SVJxphsLTMiIkIeHh52feLi4nTo0KEs1wsAAADAXp45YyFJw4YNU1RUlOrWrav69esrNjZWSUlJ6tOnjySpd+/eKlmypGJiYiRJHTt21JQpU1S7dm01aNBAe/fu1ZgxY9SxY0dbwLjRMgMDA/X4449r2LBhKly4sAICAvTUU0+pYcOGuvfee13zRAAAAAB3mDwVLLp166aTJ09q7NixOnbsmGrVqqVVq1bZLr4+dOiQ3RmK0aNHy2KxaPTo0Tpy5IiCgoLUsWNHvfzyy9lepiRNnTpVbm5u6tq1qy5duqS2bdtq+vTpt2/DAQAAgDucxRhjXF1EfpCQkKDAwEDFx8dzITcAAMj34uLi1L9/f82ePVsVK1Z0dTnIJTfzGTfPXGMBAAAA4M6Vp4ZCAQAA4M4VdyZOBxIOqELBCipXsJyry8FtRrAAAACAU5IvJ+vf3/1bPxz5wdZ2X+h9erXpq7K6W11YGW4nhkIBAADAKdO2T7MLFZK0+uBqvfvruy6qCK5AsAAAAIBTPtv/WabtK/atuM2VwJUIFgAAAHDKxdSLmbZfSL1wmyuBKxEsAAAA4JRmpZpl2t6idIvbWwhcimABAAAApwyNGKoQnxC7ttL+pTW41mAXVQRX4K5QAAAAcEpp/9Ja9q9l+mz/Z/or/i9VKFRBHcp2kI+Hj6tLw21EsAAAAIDT/Kx+6lGph6vLgAsxFAoAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGkECwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4jWAAAAABwGsECAAAAgNMIFgAAAACcRrAAAAAA4DSCBQAAAACnESwAAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGkECwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4jWAAAAABwGsECAAAAgNMIFgAAAACcRrAAAAAA4DSCBQAAAACnESwAAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABAAAAwGkECwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnFbA1QUAAAAg95y6cEpzfp2jjf9sVIBngLpU6KJO5TvdlnUfPX9Uc36boy3HtqiwV2F1q9hN95e9/7asG7cfwQIAACCfSkhJUO8ve+tw4mFb2/YT23Ug/oCejXg2V9d96sIp9fqil05cOCFJ2h+/X1uPb9XRpKPqU61Prq4brsFQKAAAgHzq0z2f2oWKDB/88YHOXjybq+v++M+PbaHiarN3zlby5eRcXTdcgzMWAAAALnTx4kUdPHgwV5a98a+NmbanpKfo253fqlpgtVtedkbNWdW++eDmTNsTLyfq+1+/V1nfsre8bmeFhobKy8vLZevPryzGGOPqIvKDhIQEBQYGKj4+XgEBAa4uBwAA3CHi4uLUv3//XFn2hYgLuljrouMEIwUsCpB7onuurFeSkholKaVyiuOENCnw40C5XXLdwJnZs2erYsWKLlv/neRmPuMSLHIIwQIAANyK3DxjceLiCQ3bOUwX0+3DRf1C9TWi4ohcWWeGQ8mHNOLXEUo1qXbtrYJa6cnwJ3N13TfCGYvsI1i4AMECAADkRdtPbNfkLZO189ROebl7qUO5DhpRb4R8PHxyfd0b/tmgKVunKO5snHwK+Khzhc4aFjFMVndrrq8bOYNg4QIECwAAkJclpiTK093TJR/qE1IS5O3uLQ93j9u+bjjnZj7jcvE2AADAXcDf6u+ydQdY+dL1bsDtZgEAAAA4jWABAAAAwGkECwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4jWAAAAABwWp4MFm+//bbCwsLk5eWlBg0aaPPmzVn2bdGihSwWi8OjQ4cOtj6ZTbdYLHr99ddtfcLCwhymv/rqq7m6nQAAAEB+UcDVBVxr4cKFGjZsmGbMmKEGDRooNjZWbdu2VVxcnIKDgx36L126VCkpKba/T58+rZo1a+rhhx+2tR09etRuni+//FKPP/64unbtatc+ceJE9e/f3/a3v79/Tm0WAAAAkK/luWAxZcoU9e/fX3369JEkzZgxQytXrtR7772nkSNHOvQvXLiw3d8LFiyQj4+PXbAoVqyYXZ/ly5erZcuWKleunF27v7+/Q9+sXLp0SZcuXbL9nZCQkK35AAAAgPwoTw2FSklJ0bZt2xQZGWlrc3NzU2RkpDZu3JitZcyZM0fdu3eXr69vptOPHz+ulStX6vHHH3eY9uqrr6pIkSKqXbu2Xn/9daWmpma5npiYGAUGBtoepUuXzlZ9AAAAQH6Up85YnDp1SmlpaQoJCbFrDwkJ0Z9//nnD+Tdv3qzffvtNc+bMybLP+++/L39/f3Xp0sWu/emnn1adOnVUuHBhbdiwQaNGjdLRo0c1ZcqUTJczatQoDRs2zPZ3QkIC4QIAAAB3rTwVLJw1Z84cVa9eXfXr18+yz3vvvadHH31UXl5edu1Xh4QaNWrIarVq4MCBiomJkaenp8NyPD09M20HAAAA7kZ5aihU0aJF5e7uruPHj9u1Hz9+/IbXPiQlJWnBggWZDnHK8MMPPyguLk79+vW7YS0NGjRQamqqDhw4kK3aAQAAgLtZngoWVqtVERER+uabb2xt6enp+uabb9SwYcPrzrto0SJdunRJvXr1yrLPnDlzFBERoZo1a96wlh07dsjNzS3TO1EBAAAAsJfnhkINGzZMUVFRqlu3rurXr6/Y2FglJSXZ7hLVu3dvlSxZUjExMXbzzZkzR506dVKRIkUyXW5CQoIWLVqkN954w2Haxo0btWnTJrVs2VL+/v7auHGjhg4dql69eqlQoUI5v5EAAABAPpPngkW3bt108uRJjR07VseOHVOtWrW0atUq2wXdhw4dkpub/YmWuLg4/fjjj/r666+zXO6CBQtkjFGPHj0cpnl6emrBggUaP368Ll26pLJly2ro0KF2110AAAAAyJrFGGNcXUR+kJCQoMDAQMXHxysgIMDV5QAAAABOu5nPuHnqGgsAAAAAdyaCBQAAAACnESwAAAAAOI1gAQAAAMBpBAsAAAAATiNYAAAAAHAawQIAAACA0/LcD+TdqTJ+DiQhIcHFlQAAAAA5I+OzbXZ++o5gkUMSExMlSaVLl3ZxJQAAAEDOSkxMVGBg4HX78MvbOSQ9PV3//POP/P39ZbFYXF0OcklCQoJKly6tw4cP8wvrgIvwOgTyBl6LdwdjjBITE1WiRAm5uV3/KgrOWOQQNzc3lSpVytVl4DYJCAjgIAq4GK9DIG/gtZj/3ehMRQYu3gYAAADgNIIFAAAAAKcRLICb4OnpqXHjxsnT09PVpQB3LV6HQN7AaxHX4uJtAAAAAE7jjAUAAAAApxEsAAAAADiNYAEAAADAaQQL5GvR0dEKCwtzdRkA8pEWLVqoRYsWtzxvtWrVcrYg3PGc2acyExYWpujo6Bxb3t1k3bp1slgsWrdunatLuSMRLJBj5s2bJ4vFYnsUKFBAJUuWVHR0tI4cOeLq8m676Ohou+fj6seqVatcXZ6Df/75R+PHj9eOHTtcXQr+v4zX1NatW+3a4+PjVb9+fXl5eWW6L2XM5+XllelrL798uE1OTtb48eMz/QAwfvx4WSwWhYSEKDk52WF6WFiYHnjggdtQJe40Wb3u8poNGzZo/PjxOnfuXI4s78CBA3bvU25ubipcuLDatWunjRs35sg6kP/xy9vIcRMnTlTZsmV18eJF/fTTT5o3b55+/PFH/fbbb/Ly8nJ1ebeVp6en3n33XYf2mjVruqCa6/vnn380YcIEhYWFqVatWq4uB1lISEhQmzZttHPnTn366ae6//77s+x76dIlvfrqq5o2bdptrPD2SU5O1oQJEyQpy297T5w4oXfeeUfPPfdcjq3366+/zrFlAdKt7VMbNmzQhAkTFB0drYIFC9pNi4uLk5vbrX133KNHD7Vv315paWnavXu3pk+frpYtW2rLli2qXr36LS3zTtKsWTNduHBBVqvV1aXckQgWyHHt2rVT3bp1JUn9+vVT0aJFNWnSJK1YsUKPPPKIi6u7vQoUKKBevXrlyrKTk5Pl4+OTK8tG3pSYmKi2bdtqx44dWrp0qdq1a3fd/rVq1dLs2bM1atQolShR4jZVmbWkpCT5+vre1nXWqlVLr7/+up588kl5e3vnyDLz+geO1NRUpaen5/k68T85/W/lzO9K1KlTx+59q2nTpmrXrp3eeecdTZ8+PSfKyzZXHDPc3Nzuui9BcxJDoZDrmjZtKknat2+fJCklJUVjx45VRESEAgMD5evrq6ZNm2rt2rV282Wclp08ebJmzZql8PBweXp6ql69etqyZYvDepYtW6Zq1arJy8tL1apV06effpppPUlJSXruuedUunRpeXp6qmLFipo8ebKu/UkXi8WiIUOGaNGiRapSpYq8vb3VsGFD/frrr5KkmTNnqnz58vLy8lKLFi104MCBW3p+pk+frqpVq8rT01MlSpTQ4MGDHU5tZwxd2bZtm5o1ayYfHx+98MILkq58Kz1u3DiVL19enp6eKl26tEaMGKFLly7ZLWP16tVq0qSJChYsKD8/P1WsWNG2jHXr1qlevXqSpD59+thOhc+bN++Wtgk57/z587r//vv1888/a8mSJerQocMN53nhhReUlpamV199NVvr+PDDDxURESFvb28VLlxY3bt31+HDh+36/PDDD3r44YdVpkwZ2/42dOhQXbhwwa5fdHS0/Pz8tG/fPrVv317+/v569NFHJUnp6emKjY1V1apV5eXlpZCQEA0cOFBnz561W8bWrVvVtm1bFS1aVN7e3ipbtqz69u0r6crxISgoSJI0YcIE2z47fvx4u2WMHTtWx48f1zvvvHPD7c9uXZmNhz948KAefPBB+fr6Kjg4WEOHDtVXX32V5Vjtbdu2qVGjRrbtmjFjhkOfEydO6PHHH1dISIi8vLxUs2ZNvf/++3Z9rj5OxsbG2o6Tf/zxxw23F7du+/btateunQICAuTn56fWrVvrp59+cui3c+dONW/eXN7e3ipVqpT+85//aO7cubJYLHbvGZntU9OmTVPVqlXl4+OjQoUKqW7dupo/f76kK0P9hg8fLkkqW7asbf/PWGZm11icO3dOQ4cOVVhYmDw9PVWqVCn17t1bp06duu62XvsefvXynn32Wdt7afny5TVp0iSlp6fb9Tt9+rQee+wxBQQEqGDBgoqKitIvv/zi8B6T28eMDAsWLFBERIT8/f0VEBCg6tWr67///a9telbXWCxatMh2fCxatKh69erlMNQ0YxuOHDmiTp06yc/PT0FBQfr3v/+ttLS06z7P+QVnLJDrMg50hQoVknRlKMe7776rHj16qH///kpMTNScOXPUtm1bbd682WEYzvz585WYmKiBAwfKYrHotddeU5cuXbR//355eHhIunIauWvXrqpSpYpiYmJ0+vRp9enTR6VKlbJbljFGDz74oNauXavHH39ctWrV0ldffaXhw4fryJEjmjp1ql3/H374QStWrNDgwYMlSTExMXrggQc0YsQITZ8+XU8++aTOnj2r1157TX379tW3337rsP3XHrQ9PDwUGBgo6cqbw4QJExQZGalBgwYpLi5O77zzjrZs2aL169fbtk+6cnBu166dunfvrl69eikkJETp6el68MEH9eOPP2rAgAGqXLmyfv31V02dOlW7d+/WsmXLJEm///67HnjgAdWoUUMTJ06Up6en9u7dq/Xr10uSKleurIkTJ2rs2LEaMGCA7Y2kUaNG2fo3Ru5KSkpSu3bttGXLFi1evDjb1waULVtWvXv31uzZszVy5MjrnrV4+eWXNWbMGD3yyCPq16+fTp48qWnTpqlZs2bavn27bajFokWLlJycrEGDBqlIkSLavHmzpk2bpr///luLFi2yW2Zqaqratm2rJk2aaPLkybYzbAMHDtS8efPUp08fPf300/rrr7/01ltvafv27bb9/sSJE2rTpo2CgoI0cuRIFSxYUAcOHNDSpUslSUFBQXrnnXc0aNAgde7cWV26dJEk1ahRw66Gpk2bqlWrVnrttdc0aNCg6561yE5dmUlKSlKrVq109OhRPfPMMypWrJjmz5/v8GVJhrNnz6p9+/Z65JFH1KNHD33yyScaNGiQrFar7UPQhQsX1KJFC+3du1dDhgxR2bJltWjRIkVHR+vcuXN65pln7JY5d+5cXbx4UQMGDJCnp6cKFy6c5XbCOb///ruaNm2qgIAAjRgxQh4eHpo5c6ZatGih7777Tg0aNJAkHTlyRC1btpTFYtGoUaPk6+urd999N1tnE2bPnq2nn35aDz30kJ555hldvHhRO3fu1KZNm9SzZ0916dJFu3fv1scff6ypU6eqaNGikmQL29c6f/68mjZtql27dqlv376qU6eOTp06pRUrVujvv/+2zZ+Za9/DpStnzJs3b64jR45o4MCBKlOmjDZs2KBRo0bp6NGjio2NlXQlEHTs2FGbN2/WoEGDVKlSJS1fvlxRUVGZris3jxnSlS/YevToodatW2vSpEmSpF27dmn9+vUOr6mrZay7Xr16iomJ0fHjx/Xf//5X69evtzs+SlJaWpratm2rBg0aaPLkyVqzZo3eeOMNhYeHa9CgQVmuI98wQA6ZO3eukWTWrFljTp48aQ4fPmwWL15sgoKCjKenpzl8+LAxxpjU1FRz6dIlu3nPnj1rQkJCTN++fW1tf/31l5FkihQpYs6cOWNrX758uZFkPvvsM1tbrVq1TPHixc25c+dsbV9//bWRZEJDQ21ty5YtM5LMf/7zH7v1P/TQQ8ZisZi9e/fa2iQZT09P89dff9naZs6caSSZYsWKmYSEBFv7qFGjjCS7vlFRUUaSw6N58+bGGGNOnDhhrFaradOmjUlLS7PN99ZbbxlJ5r333rO1NW/e3EgyM2bMsKv7gw8+MG5ubuaHH36wa58xY4aRZNavX2+MMWbq1KlGkjl58qTJypYtW4wkM3fu3Cz74PbKeE2FhoYaDw8Ps2zZspuab8uWLWbfvn2mQIEC5umnn7ZNb968ualatart7wMHDhh3d3fz8ssv2y3n119/NQUKFLBrT05OdlhfTEyMsVgs5uDBg7a2jP1/5MiRdn1/+OEHI8l89NFHdu2rVq2ya//0009t25CVkydPGklm3LhxDtPGjRtn2+e/++47I8lMmTLFNj00NNR06NDhpusy5srzl/E6NsaYN954w0iy+/e5cOGCqVSpkpFk1q5dazevJPPGG2/Y2i5dumRq1aplgoODTUpKijHGmNjYWCPJfPjhh7Z+KSkppmHDhsbPz892/Mk4TgYEBJgTJ05k+Vwh+65+/WSmU6dOxmq1mn379tna/vnnH+Pv72+aNWtma3vqqaeMxWIx27dvt7WdPn3aFC5c2OH94tp96l//+pfdazQzr7/+usNyMoSGhpqoqCjb32PHjjWSzNKlSx36pqenG2P+ty9NmDDBnDx50hw7dsz88MMPpl69ekaSWbRokW2el156yfj6+prdu3fbLWvkyJHG3d3dHDp0yBhjzJIlS4wkExsba+uTlpZmWrVq5fB+czuOGc8884wJCAgwqampWfZZu3at3es2JSXFBAcHm2rVqpkLFy7Y+n3++edGkhk7dqzDNkycONFumbVr1zYRERFZrjM/YSgUclxkZKSCgoJUunRpPfTQQ/L19dWKFStsZw/c3d1t40nT09N15swZpaamqm7duvr5558dltetWze7b0oyvk3fv3+/JOno0aPasWOHoqKibGcCJOm+++5TlSpV7Jb1xRdfyN3dXU8//bRd+3PPPSdjjL788ku79tatW9vdrjbjm6iuXbvK39/foT2jpgxeXl5avXq13eONN96QJK1Zs0YpKSl69tln7S6y69+/vwICArRy5Uq7ZXl6eqpPnz52bYsWLVLlypVVqVIlnTp1yvZo1aqVJNm+Mc34NmX58uUOp6mR9x0/flxeXl4qXbr0Tc9brlw5PfbYY5o1a5aOHj2aaZ+lS5cqPT1djzzyiN1+VKxYMVWoUMHum/erv/FPSkrSqVOn1KhRIxljtH37dodlX/sN3aJFixQYGKj77rvPbl0RERHy8/Nz2Gc///xzXb58+aa3+2rNmjVTy5Yt9dprrzkM2brZujKzatUqlSxZUg8++KCtzcvLS/3798+0f4ECBTRw4EDb31arVQMHDtSJEye0bds2SVeOVcWKFVOPHj1s/Tw8PPT000/r/Pnz+u677+yW2bVr1yy/rUbOSUtL09dff61OnTqpXLlytvbixYurZ8+e+vHHH5WQkCDpyn7RsGFDu7PwhQsXtg3vuZ6CBQvq77//znTY761YsmSJatasqc6dOztMs1gsdn+PGzdOQUFBKlasmO0sxxtvvKGHHnrI1mfRokVq2rSpChUqZPd6iYyMVFpamr7//ntJV54DDw8Pu9eCm5ubbRRAZnLzmFGwYEElJSVp9erV13m27G3dulUnTpzQk08+aXftRYcOHVSpUiWH92pJeuKJJ+z+btq0qcPng/yKYIEc9/bbb2v16tVavHix2rdvr1OnTjmc+n3//fdVo0YNeXl5qUiRIgoKCtLKlSsVHx/vsLwyZcrY/Z0RMjLGVh48eFCSVKFCBYd5K1asaPf3wYMHVaJECbtQIF0ZCnT1srJad0ZwufYDXkb7teM93d3dFRkZafeIiIiwW9e1NVqtVpUrV86hlpIlSzpc4Ldnzx79/vvvCgoKsnvcc889kq6M0ZauhLPGjRurX79+CgkJUffu3fXJJ58QMu4QM2fOlNVq1f3336+4uDhJVz7gHDt2zO6RkpKS6fyjR49Wampqltda7NmzR8YYVahQwWFf2rVrl20/kqRDhw4pOjpahQsXto0fbt68uSQ5vH4LFCjgMBxxz549io+PV3BwsMO6zp8/b1tX8+bN1bVrV02YMEFFixbVv/71L82dO9fh2qHsGj9+vI4dO5bptQw3U1dmDh48qPDwcIcPaOXLl8+0f4kSJRwuSM14zWYMOzl48KAqVKjgcGefrI5VZcuWzbI+5JyTJ08qOTnZ4bgtXfm3SU9Pt12XdPDgwUz3gaz2i6s9//zz8vPzU/369VWhQgUNHjzYNnT1Vuzbty/bt5geMGCAVq9erc8++8x2/dS11wfs2bNHq1atcnitREZGSvrfe8/BgwdVvHhxhxuNZPUc5PYx48knn9Q999yjdu3aqVSpUurbt+8Nb/+e1Xu1JFWqVMnhtejl5eUQ8gsVKuTw+SC/4hoL5Lj69evb7grVqVMnNWnSRD179lRcXJz8/Pz04YcfKjo6Wp06ddLw4cMVHBwsd3d3xcTEOFwcJl35cJ4Zc83F1rkhq3W7oqbMxoanp6erevXqmjJlSqbzZAQgb29vff/991q7dq1WrlypVatWaeHChWrVqpW+/vrrLLcHeUOVKlX0xRdfqHXr1rrvvvu0fv16paWlOXyYXLt2baa3XS1Xrpx69eqlWbNmaeTIkQ7T09PTZbFY9OWXX2a6L/j5+Um6Embuu+8+nTlzRs8//7wqVaokX19fHTlyRNHR0Q5B1dPT0+GDcXp6uoKDg/XRRx9luq0Zb8gWi0WLFy/WTz/9pM8++0xfffWV+vbtqzfeeEM//fSTrabsatasmVq0aKHXXnvN4dvEm6krr8qpO14hb6hcubLi4uL0+eefa9WqVVqyZImmT5+usWPH2m6xnFsqVKhgCwgPPPCA3N3dNXLkSLVs2dL23p6enq777rtPI0aMyHQZGUH5ZuX2MSM4OFg7duzQV199pS+//FJffvml5s6dq969ezvcGOFW3e3vpwQL5KqMwNCyZUu99dZbGjlypBYvXqxy5cpp6dKldt/wjRs37pbWERoaKunKtxrXyvh29+q+a9asUWJiot1Ziz///NNuWbdDxrri4uLsTqmnpKTor7/+sh3Yryc8PFy//PKLWrdu7fBt6bXc3NzUunVrtW7dWlOmTNErr7yiF198UWvXrlVkZOQN54dr1a9fX8uWLVOHDh1033332YbWXe16v48yevRoffjhh7YLFq8WHh4uY4zKli173Q8Ev/76q3bv3q33339fvXv3trXfzLCC8PBwrVmzRo0bN87Wh+F7771X9957r15++WXNnz9fjz76qBYsWKB+/frd9D47fvx4tWjRQjNnznS6rquFhobqjz/+kDHGrqa9e/dm2v+ff/5xuI3m7t27Jck29DI0NFQ7d+5Uenq63QctVxyr8D9BQUHy8fFxeG+RrvzbuLm52b7QCQ0NzXQfyGq/uJavr6+6deumbt26KSUlRV26dNHLL7+sUaNGycvL66b2//DwcP3222/Z7n+1F198UbNnz9bo0aNt3+6Hh4fr/PnzN3yfCg0N1dq1ax1uj57d5yBjXTl1zJCujAro2LGjOnbsqPT0dD355JOaOXOmxowZk+mZlKvfqzOGGWeIi4vjtXgNhkIh17Vo0UL169dXbGysLl68aEvzV3+7v2nTplv+Zc/ixYurVq1aev/99+2GYqxevdrhlosZP/rz1ltv2bVPnTpVFovlhr8LkJMiIyNltVr15ptv2j0Xc+bMUXx8fLZuJ/rII4/oyJEjmj17tsO0CxcuKCkpSZJ05swZh+kZ434zThNnfMjJqV9xRc5r3bq1Pv74Y+3du1edOnVS/fr17YbZXX0t0rXCw8PVq1cvzZw5U8eOHbOb1qVLF7m7u2vChAkOZ92MMTp9+rQkZfraNcbY3arxRh555BGlpaXppZdecpiWmppq2//Onj3rUMu1+2zGB5Xs7rPNmzdXixYtNGnSJF28ePGW6spM27ZtdeTIEa1YscLWdvHixUxflxnLuzrcpKSkaObMmQoKCrINlWzfvr2OHTumhQsX2s03bdo0+fn52Yaf4fZyd3dXmzZttHz5crvbxR4/flzz589XkyZNFBAQIOnKfrFx40bt2LHD1u/MmTNZfvN+tYzXXAar1aoqVarIGGO7fuBmjtldu3bVL7/8kult2G90pr1gwYIaOHCgvvrqK9u2PPLII9q4caO++uorh/7nzp1TamqqpCvPweXLl+1eC+np6Xr77bdvWHOGnDxmXPu8urm52e4kl9Uwy7p16yo4OFgzZsyw6/Pll19q165d2XqvvptwxgK3xfDhw/Xwww9r3rx5euCBB7R06VJ17txZHTp00F9//aUZM2aoSpUqOn/+/C0tPyYmRh06dFCTJk3Ut29fnTlzxnYP8KuX2bFjR7Vs2VIvvviiDhw4oJo1a+rrr7/W8uXL9eyzzyo8PDynNvmGgoKCNGrUKE2YMEH333+/HnzwQcXFxWn69OmqV69etn5Y77HHHtMnn3yiJ554QmvXrlXjxo2VlpamP//8U5988om++uor1a1bVxMnTtT333+vDh06KDQ0VCdOnND06dNVqlQpNWnSRNKVD54FCxbUjBkz5O/vL19fXzVo0ICx23lM586dNXv2bPXt21cPPvigVq1ale0fc3rxxRf1wQcfKC4uTlWrVrW1h4eH6z//+Y9GjRqlAwcOqFOnTvL399dff/2lTz/9VAMGDNC///1vVapUSeHh4fr3v/+tI0eOKCAgQEuWLLmpscPNmzfXwIEDFRMTox07dqhNmzby8PDQnj17tGjRIv33v//VQw89pPfff1/Tp09X586dFR4ersTERM2ePVsBAQFq3769pCvDf6pUqaKFCxfqnnvuUeHChVWtWrXrjiUfN26cWrZsect1ZWbgwIF666231KNHDz3zzDMqXry4PvroI9u/y7XfLJcoUUKTJk3SgQMHdM8992jhwoXasWOHZs2aZbul7YABAzRz5kxFR0dr27ZtCgsL0+LFi7V+/XrFxsY6XCeGnPfee+9lOv5+/Pjxtt8FevLJJ1WgQAHNnDlTly5d0muvvWbrN2LECH344Ye677779NRTT9luN1umTBmdOXPmumcc2rRpo2LFiqlx48YKCQnRrl279NZbb6lDhw62f/uMEPriiy+qe/fu8vDwUMeOHTP9Qbnhw4dr8eLFevjhh9W3b19FRETozJkzWrFihWbMmHHds52S9Mwzzyg2NlavvvqqFixYoOHDh2vFihV64IEHFB0drYiICCUlJenXX3/V4sWLdeDAARUtWtT2Bchzzz2nvXv3qlKlSlqxYoXty67snHXJyWNGv379dObMGbVq1UqlSpXSwYMHNW3aNNWqVct2/dK1PDw8NGnSJPXp00fNmzdXjx49bLebDQsL09ChQ2+4DXeV234fKuRb17tFX1pamgkPDzfh4eEmNTXVvPLKKyY0NNR4enqa2rVrm88//9xERUXZ3Ro249Z3r7/+usPylMktJpcsWWIqV65sPD09TZUqVczSpUsdlmmMMYmJiWbo0KGmRIkSxsPDw1SoUMG8/vrrtlvuXb2OwYMH27VlVVPG7emuvh1fVFSU8fX1vd5TZoy5cnvZSpUqGQ8PDxMSEmIGDRpkzp49a9fn2tuDXi0lJcVMmjTJVK1a1Xh6eppChQqZiIgIM2HCBBMfH2+MMeabb74x//rXv0yJEiWM1Wo1JUqUMD169HC4VeDy5ctNlSpVTIECBbj1bB5wvdfU5MmTjSTzwAMPmMuXL2d7vozbIWa2Py1ZssQ0adLE+Pr6Gl9fX1OpUiUzePBgExcXZ+vzxx9/mMjISOPn52eKFi1q+vfvb3755ZdMbx15vf1/1qxZJiIiwnh7ext/f39TvXp1M2LECPPPP/8YY4z5+eefTY8ePUyZMmWMp6enCQ4ONg888IDZunWr3XI2bNhgIiIijNVqtTsuXH272Wtl3PL16tvNZreujPmvvjWoMcbs37/fdOjQwXh7e5ugoCDz3HPP2W61+dNPP9nNW7VqVbN161bTsGFD4+XlZUJDQ81bb73lUMvx48dNnz59TNGiRY3VajXVq1d3eE1e7ziJW5Px+snqcfjwYfPzzz+btm3bGj8/P+Pj42NatmxpNmzY4LCs7du3m6ZNmxpPT09TqlQpExMTY958800jyRw7dszW79p9aubMmaZZs2amSJEixtPT04SHh5vhw4fbjukZXnrpJVOyZEnj5uZmd+vZa283a8yVW90OGTLElCxZ0litVlOqVCkTFRVlTp06ZYy58b4UHR1t3N3dbbdlT0xMNKNGjTLly5c3VqvVFC1a1DRq1MhMnjzZdttkY67cFrpnz57G39/fBAYGmujoaLN+/XojySxYsMDW73YcMxYvXmzatGljgoODjdVqNWXKlDEDBw40R48etfW59nazGRYuXGhq165tPD09TeHChc2jjz5q/v77b7s+WW1DxvHobmAx5jZcAQsAwF0mNjZWQ4cO1d9//62SJUu6uhzkEc8++6xmzpyp8+fP37UX+i5btkydO3fWjz/+qMaNG7u6HOQgggUAAE66cOGC3YWlFy9eVO3atZWWlma7MBt3n2v3i9OnT+uee+5RnTp1buqmB3eya5+DtLQ0tWnTRlu3btWxY8e4o1k+wzUWAAA4qUuXLipTpoxq1aql+Ph4ffjhh/rzzz+zdaEu8q+GDRuqRYsWqly5so4fP645c+YoISFBY8aMcXVpt81TTz2lCxcuqGHDhrp06ZKWLl2qDRs26JVXXiFU5EMECwAAnNS2bVu9++67+uijj5SWlqYqVapowYIF6tatm6tLgwu1b99eixcv1qxZs2SxWFSnTh3NmTNHzZo1c3Vpt02rVq30xhtv6PPPP9fFixdVvnx5TZs2TUOGDHF1acgFDIUCAAAA4DR+xwIAAACA0wgWAAAAAJxGsAAAAADgNIIFAAAAAKcRLAAAAAA4jWABALgtoqOjZbFYXF0GACCXECwA4C62bt06WSyWLB8FCtzczx3NmzdPsbGxuVNsDhk/fryWLVvm6jIAIN/hdywA4C62bt06tWzZUj169FD79u0dpru5ualnz57ZXl6LFi104MABHThwwGHa5cuXlZaWJi8vL2dKdprFYlFUVJTmzZvn0joAIL/hl7cBAKpTp4569eqVq+vw8PCQh4dHrq7D1RITE+Xv7+/qMgDAJRgKBQDIlv/7v/9T/fr1VbBgQfn6+qpcuXJ69NFHdfLkSUlSWFiYvvvuOx08eNBuONW6deskZX6NRUbb6dOnFR0draJFi8rf31+dOnXSsWPHJEmzZs1S5cqV5eXlpUqVKmn58uUOtU2fPl1t2rRRyZIlZbVaVbx4cfXq1cvuzMmBAwds63///fftarzau+++qzp16sjb21uBgYFq06aNfvzxR4d1WiwWRUdH65tvvlGTJk3k5+enjh07SpLOnDmjoUOHKjw8XF5eXipSpIgiIiL0+uuv39qTDwB3AM5YAACUnJysU6dOObRbrVYFBATogw8+UFRUlJo2baqJEyfK29tbhw8f1hdffKETJ04oKChIsbGxGjVqlE6dOqWpU6fallG5cuUbrv/+++9XqVKlNHHiRO3du1dvvvmmOnfurC5dumjWrFl6/PHH5eXlpTfffFMPPfSQdu/erbJly9rmnzx5su699149/fTTKly4sH777Te9++67+vbbb/Xrr7+qSJEiCgoK0gcffKDHHntMTZs21YABAxzqeP755/Xaa6+pfv36euWVV5SYmKhZs2apZcuWWr58ucNwsa1bt2rJkiXq37+/oqKibO0PP/ywvv/+ez3xxBOqUaOGLly4oF27dmndunUaPnx4tv5NAOCOYwAAd621a9caSVk+OnToYIwxpnPnzsbf399cvnz5ustr3ry5CQ0NzXRaVFSUufZtJ6PtySeftGsfOnSokWRKly5t4uPjbe2//PKLkWRGjhxp1//8+fMO61uzZo2RZCZNmmTXLslERUU59P/zzz+NxWIxjRs3NpcuXbK1HzlyxAQGBprQ0FCTmppqtxxJZvXq1XbLOXfunJFkBg0alOnzAAD5FWcsAAAaMGCAHn74YYf2oKAgSVJgYKCSk5O1cuVKPfjggzl+29hnn33W7u+mTZtq6tSp6t27twICAmztNWrUUEBAgPbs2WPX39fXV5KUnp6uxMREXb58WTVr1lRgYKA2bdqUrRqWL18uY4xGjBghq9Vqay9RooT69Omj2NhYbd++XXXr1rVNq1mzpiIjI+2W4+3tLU9PT23atEkHDhxQWFhYttYPAHc6ggUAQBUqVHD4gHy1F154Qd9//706deqkIkWKqHnz5mrXrp26deuWIxcrlytXzu7vQoUKSZLdcKerp50+fdqu7dtvv9XEiRO1adMmXbx40W7a2bNns1XDX3/9JUmqWrWqw7SMtv3799sFi3vuucehr9VqVWxsrJ555hmVLVtWVapUUatWrdSpUye1bt06W7UAwJ2Ii7cBADdUoUIF/fHHH1q5cqWioqJ08OBB9e/fX5UqVdK+ffucXr67u/tNtZur7pS+ZcsWtWnTRseOHdOrr76q5cuX6+uvv9bq1atVpEgRpaenO11fVnx8fDJtf+KJJ3TgwAHNnj1bderU0eLFixUZGanu3bvnWi0A4GqcsQAAZIunp6fat29vu4D5iy++UIcOHTRlyhS9/fbbkuSSX9aeP3++0tLS9OWXX9qd4UhKSsr22Qrpf2dNfv/9d4WHh9tN++OPP+z6ZEfx4sXVr18/9evXT2lpaXrsscf08ccf67nnnlO9evWyvRwAuFNwxgIAcEOZ3TGqTp06kq7cWjWDn5+fzp49a3dGIbdlnNW4dp2vvPJKpmcr/Pz87GrOkHHtyOuvv67Lly/b2o8ePaq5c+cqNDRUtWvXvmE9ycnJSk5OdqixRo0akpTpugEgP+CMBQBAP//8sz788MNMp3Xq1Elt2rRRwYIF1bRpU5UuXVrnzp3TvHnzZLFY9Nhjj9n63nvvvfr88881ZMgQNWrUSO7u7mrVqpWCg4NzrfbOnTtr6tSpat++vQYMGCCr1arVq1dr586dKlq0qEP/e++9V2vWrNGkSZNUpkwZWSwWde/eXRUrVtTw4cP12muvqVmzZurWrZvtdrPnz5/XRx99lOXQrKvt3r1bzZs3V+fOnVWtWjUVKlRIu3bt0jvvvKOyZcuqadOmufE0AIDLESwAAPr444/18ccfZzptz549GjRokD755BPNnDlTZ86cUZEiRVS7dm1NmzZNLVu2tPUdOnSo9u/fr8WLF2vGjBlKT0/X2rVrczVYNG7cWEuWLNFLL72kMWPGyNvbW5GRkfruu+/UrFkzh/7Tp0/X4MGD9fLLLysxMVGSbNc+TJo0SeXLl9f06dM1cuRIWa1WNWjQQPPnz892IChdurT69u2rtWvXatmyZbp06ZJKliyp/v376/nnn8/yugwAuNNZzO08Xw0AAAAgX+IaCwAAAABOI1gAAAAAcBrBAgAAAIDTCBYAAAAAnEawAAAAAOA0ggUAAAAApxEsAAAAADiNYAEAAADAaQQLAAAAAE4jWAAAAABwGsECAAAAgNMIFgAAAACc9v8A87UcwkvDXqgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x900 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(8, 9))\n",
    "\n",
    "# plot the validation set performance\n",
    "ax = sns.boxplot(\n",
    "    data=df_cv,\n",
    "    x=\"estimator\",\n",
    "    y=\"mean_test_score\",\n",
    "    hue=\"estimator\",\n",
    "    boxprops={\"alpha\": 0.4},\n",
    "    showmeans=False,\n",
    ")\n",
    "sns.stripplot(\n",
    "    data=df_cv,\n",
    "    x=\"estimator\",\n",
    "    y=\"mean_test_score\",\n",
    "    hue=\"estimator\",\n",
    "    ax=ax,\n",
    ")\n",
    "\n",
    "# plot the test set performance of the best models\n",
    "for model_name, metric_value in test_roc_aucs.items():\n",
    "    plt.scatter(\n",
    "        rename_dict[model_name], metric_value, marker=\"x\", color=\"red\", zorder=10\n",
    "    )\n",
    "\n",
    "plt.title(\"Performance on the validation and test set\", fontsize=16)\n",
    "plt.ylabel(\"ROC AUC\", fontsize=13)\n",
    "plt.xlabel(\"Estimators\", fontsize=13)\n",
    "plt.xticks(fontsize=12)\n",
    "plt.tight_layout()\n",
    "# plt.savefig(\"hyperopt_performance.pdf\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bd3b7c5-0da0-4752-9d3d-756b7c053813",
   "metadata": {},
   "source": [
    "The results illustrate the ROC AUC distribution of the models' hyperparameter configurations during the optimization on the validation set and after the parameter selection on the test set. The models vary in their performance. Overall, the prediction performance on the BBBP data set is very good, even with KNN and logistic regression."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}