nbs/GT-Antibiotics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data from the paper: [Discovery of a structural class of antibiotics with explainable deep learning](https://www.nature.com/articles/s41586-023-06887-8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/pawel/anaconda3/envs/gt/lib/python3.7/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Numpy version: 1.21.6\n",
      "Rdkit version: 2022.09.5\n",
      "Torch version: 1.13.1\n",
      "TorchMetrics version: 0.11.4\n"
     ]
    }
   ],
   "source": [
    "# Standard\n",
    "import logging\n",
    "import importlib\n",
    "\n",
    "# Third party\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import rdkit\n",
    "from rdkit import RDLogger\n",
    "import torch\n",
    "from torch import nn\n",
    "import matplotlib.pyplot as plt\n",
    "from torch.optim.lr_scheduler import ReduceLROnPlateau\n",
    "from torch_geometric.loader import DataLoader\n",
    "import torchmetrics\n",
    "from torchmetrics import AUROC\n",
    "\n",
    "# GT-PyG\n",
    "import gt_pyg\n",
    "from gt_pyg.data.utils import (\n",
    "    get_tensor_data,\n",
    "    get_node_dim,\n",
    "    get_edge_dim,\n",
    "    clean_df,\n",
    ")\n",
    "from gt_pyg.nn.model import GraphTransformerNet\n",
    "\n",
    "\n",
    "# Turn off majority of RDKit warnings\n",
    "lg = RDLogger.logger()\n",
    "lg.setLevel(RDLogger.CRITICAL)\n",
    "\n",
    "\n",
    "# Set a random seed for a reproducibility purposes\n",
    "torch.manual_seed(1928374650)\n",
    "\n",
    "# Setup the logger\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)\n",
    "\n",
    "# Log the used versions of RDkit and torch\n",
    "print(f'Numpy version: {np.__version__}')\n",
    "print(f'Rdkit version: {rdkit.__version__}')\n",
    "print(f'Torch version: {torch.__version__}')\n",
    "print(f'TorchMetrics version: {torchmetrics.__version__}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The training dataset\n",
    "\n",
    "**Note**: 39,312 compounds for the antibacterial activity againts _S. aureus_ (Gram-positive pathogen)\n",
    "\n",
    "The training set was downloaded from the [original repo](https://github.com/felixjwong/antibioticsai) on 12/31/2023.\n",
    "\n",
    "Reference: _Wong et al., Nature, 2023_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>smiles</th>\n",
       "      <th>activity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Nc1nnc(o1)-c1ccc(o1)[N+](=O)[O-]</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O[C@H]1COC[C@@H]2O[C@H](CC[C@H]2N(C1)C(=O)Nc1c...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>CC(C)C[C@@H](N)C(=O)N[C@@H]1[C@H](O)c2ccc(c(c2...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[O-][N+](=O)c1ccc(o1)/C=N/N1CC(=O)NC1=O</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Cn1cnc(c1)CCNC(=O)C[C@@H]1CC[C@@H]2[C@H](COC[C...</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                              smiles  activity\n",
       "0                   Nc1nnc(o1)-c1ccc(o1)[N+](=O)[O-]         1\n",
       "1  O[C@H]1COC[C@@H]2O[C@H](CC[C@H]2N(C1)C(=O)Nc1c...         1\n",
       "2  CC(C)C[C@@H](N)C(=O)N[C@@H]1[C@H](O)c2ccc(c(c2...         1\n",
       "3            [O-][N+](=O)c1ccc(o1)/C=N/N1CC(=O)NC1=O         1\n",
       "4  Cn1cnc(c1)CCNC(=O)C[C@@H]1CC[C@@H]2[C@H](COC[C...         1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset_filename = './data/train.antibioticsai.csv'\n",
    "df = pd.read_csv(dataset_filename, sep=',')\n",
    "df = df.rename(columns={'SMILES': 'smiles', 'ACTIVITY': 'activity'})\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The paper seems to be employing either a random or stratified split for evaluating the model. In the exercise below, we opt for the random split, as when N approaches infinity, the statistics for both random and stratified splits converge. We generate training, validation, and test sets, maintaining the same proportions as outlined in the original paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_random_split(df, frac=[0.8, 0.1, 0.1], min_num_atoms=0):\n",
    "    \n",
    "    # Assert that the sum of fractions is equal to 1\n",
    "    assert sum(frac) == 1, \"The sum of fractions must be equal to 1\"\n",
    "\n",
    "    # Shuffle the DataFrame rows to get a random split\n",
    "    shuffled_df = df.sample(frac=1, random_state=42)\n",
    "\n",
    "    # Calculate the number of rows for each set\n",
    "    total_rows = len(shuffled_df)\n",
    "    train_rows = int(frac[0] * total_rows)\n",
    "    val_rows = int(frac[1] * total_rows)\n",
    "\n",
    "    # Split the DataFrame\n",
    "    train_set = shuffled_df.iloc[:train_rows]\n",
    "    val_set = shuffled_df.iloc[train_rows:train_rows + val_rows]\n",
    "    test_set = shuffled_df.iloc[train_rows + val_rows:]\n",
    "\n",
    "    train_set = clean_df(train_set, min_num_atoms=min_num_atoms, x_label='smiles', y_label='activity')\n",
    "    val_set = clean_df(val_set, min_num_atoms=min_num_atoms, x_label='smiles', y_label='activity')\n",
    "    test_set = clean_df(test_set, min_num_atoms=min_num_atoms, x_label='smiles', y_label='activity')\n",
    "    \n",
    "    return (train_set, val_set, test_set)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the minimum number of atoms in a molecule to 4, as we utilize 6 eigenvectors for spectral positional encoding. It is important to note that this exclusion affects compounds in both the training and validation sets, but not in the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:root:Removed 7 compounds that did not meet the size criteria.\n",
      "INFO:root:Removed 1 compounds that did not meet the size criteria.\n",
      "INFO:root:Removed 0 compounds that did not meet the size criteria.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of training examples: 31442\n",
      "Number of validation examples: 3930\n",
      "Number of test examples: 3932\n"
     ]
    }
   ],
   "source": [
    "PE_DIM = 6\n",
    "(tr, va, te) = get_random_split(df, frac=[0.8, 0.1, 0.1], min_num_atoms=4)\n",
    "tr_dataset = get_tensor_data(tr.smiles.to_list(), tr.activity.to_list(), pe_dim=PE_DIM)\n",
    "va_dataset = get_tensor_data(va.smiles.to_list(), va.activity.to_list(), pe_dim=PE_DIM)\n",
    "te_dataset = get_tensor_data(te.smiles.to_list(), te.activity.to_list(), pe_dim=PE_DIM)\n",
    "NODE_DIM = get_node_dim()\n",
    "EDGE_DIM = get_edge_dim()\n",
    "\n",
    "print(f'Number of training examples: {len(tr_dataset)}')\n",
    "print(f'Number of validation examples: {len(va_dataset)}')\n",
    "print(f'Number of test examples: {len(te_dataset)}')\n",
    "\n",
    "train_loader = DataLoader(tr_dataset, batch_size=256)\n",
    "val_loader = DataLoader(va_dataset, batch_size=1024)\n",
    "test_loader = DataLoader(te_dataset, batch_size=1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Removed 8 compounds.\n"
     ]
    }
   ],
   "source": [
    "original_number_of_cmpds = len(df)\n",
    "num_tr = len(tr)\n",
    "num_ve = len(va)\n",
    "num_te = len(te)\n",
    "print(f'Removed {original_number_of_cmpds - num_tr - num_ve - num_te} compounds.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train and eval the GT model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Auxiliary functions\n",
    "\n",
    "Due to the significant imbalance in the training set, we apply an upweighting technique to the positive class during the calculation of the loss value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import precision_recall_curve, auc\n",
    "\n",
    "from torch.nn.functional import one_hot\n",
    "\n",
    "def calculate_class_ratio(labels):\n",
    "    pos_class = sum(labels)\n",
    "    neg_class = len(labels) - pos_class\n",
    "    \n",
    "    ratio = max(1, neg_class) / max(1, pos_class)\n",
    "    return torch.tensor(ratio)\n",
    "\n",
    "\n",
    "def train(epoch):\n",
    "    model.train()\n",
    "    auroc = AUROC(task=\"binary\")\n",
    "\n",
    "    for data in train_loader:\n",
    "        pos_weight = calculate_class_ratio(data.y)\n",
    "        loss_func = nn.BCEWithLogitsLoss(pos_weight=pos_weight)\n",
    "        \n",
    "        data = data.to(device)\n",
    "        optimizer.zero_grad()\n",
    "        \n",
    "        # randomly flip sign of eigenvectors\n",
    "        batch_pe = data.pe * (2 * torch.randint(low=0, high=2, size=(1, PE_DIM)).float() - 1.0)\n",
    "        (out, _) = model(data.x, data.edge_index, data.edge_attr, batch_pe, data.batch, zero_var=True)\n",
    "        loss = loss_func(out.squeeze(), data.y)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        \n",
    "        auroc.update(out.squeeze(), data.y)\n",
    "\n",
    "    return auroc.compute()\n",
    "\n",
    "\n",
    "@torch.no_grad()\n",
    "def test(loader):\n",
    "    model.eval()\n",
    "    auroc = AUROC(task=\"binary\")\n",
    "    \n",
    "    for data in loader:\n",
    "        data = data.to(device)\n",
    "        # randomly flip sign of eigenvectors\n",
    "        batch_pe = data.pe * (2 * torch.randint(low=0, high=2, size=(1, PE_DIM)).float() - 1.0)\n",
    "        (out,_) = model(data.x, data.edge_index, data.edge_attr, batch_pe, data.batch)\n",
    "        \n",
    "        auroc.update(out.squeeze(), data.y)\n",
    "        \n",
    "    return auroc.compute()\n",
    "\n",
    "@torch.no_grad()\n",
    "def calc_auprc(loader):\n",
    "    (scores, labels) = score(loader)\n",
    "    precision, recall, thresholds = precision_recall_curve(labels, scores)\n",
    "    return (auc(recall, precision), scores, labels)\n",
    "\n",
    "@torch.no_grad()\n",
    "def score(loader):\n",
    "    model.eval()\n",
    "\n",
    "    all_logits = []\n",
    "    all_labels = []\n",
    "\n",
    "    with torch.no_grad():\n",
    "        for data in loader:\n",
    "            data = data.to(device)\n",
    "            batch_pe = data.pe * (2 * torch.randint(low=0, high=2, size=(1, PE_DIM)).float() - 1.0)\n",
    "            logits, _ = model(data.x, data.edge_index, data.edge_attr, batch_pe, data.batch)\n",
    "            all_logits.append(logits.cpu().numpy())\n",
    "            all_labels.append(data.y.cpu().numpy())\n",
    "\n",
    "    all_logits = np.concatenate(all_logits)\n",
    "    all_labels = np.concatenate(all_labels)\n",
    "\n",
    "    return (all_logits, all_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GraphTransformer architecture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GraphTransformerNet(\n",
      "  (node_emb): Linear(in_features=76, out_features=128, bias=False)\n",
      "  (edge_emb): Linear(in_features=10, out_features=128, bias=False)\n",
      "  (pe_emb): Linear(in_features=6, out_features=128, bias=False)\n",
      "  (gt_layers): ModuleList(\n",
      "    (0): GTConv(128, 128, heads=8, aggrs: sum,mean, qkv_bias: False, gate: True)\n",
      "    (1): GTConv(128, 128, heads=8, aggrs: sum,mean, qkv_bias: False, gate: True)\n",
      "    (2): GTConv(128, 128, heads=8, aggrs: sum,mean, qkv_bias: False, gate: True)\n",
      "    (3): GTConv(128, 128, heads=8, aggrs: sum,mean, qkv_bias: False, gate: True)\n",
      "  )\n",
      "  (global_pool): MultiAggregation([\n",
      "    SumAggregation(),\n",
      "    MeanAggregation(),\n",
      "    MaxAggregation(),\n",
      "    StdAggregation(),\n",
      "  ], mode=cat)\n",
      "  (mu_mlp): MLP(\n",
      "    (mlp): Sequential(\n",
      "      (0): Linear(in_features=512, out_features=128, bias=True)\n",
      "      (1): GELU(approximate='none')\n",
      "      (2): Linear(in_features=128, out_features=1, bias=True)\n",
      "    )\n",
      "  )\n",
      "  (log_var_mlp): MLP(\n",
      "    (mlp): Sequential(\n",
      "      (0): Linear(in_features=512, out_features=128, bias=True)\n",
      "      (1): GELU(approximate='none')\n",
      "      (2): Linear(in_features=128, out_features=1, bias=True)\n",
      "    )\n",
      "  )\n",
      ")\n",
      "Number of params: 1004 k\n",
      "Epoch: 01, AUPRC= Train: 0.1064 (AUROC= 0.8207), Val: 0.0933, Test: 0.1451\n",
      "Epoch: 02, AUPRC= Train: 0.1535 (AUROC= 0.8880), Val: 0.1434, Test: 0.2012\n",
      "Epoch: 03, AUPRC= Train: 0.2014 (AUROC= 0.9065), Val: 0.1513, Test: 0.2662\n",
      "Epoch: 04, AUPRC= Train: 0.2087 (AUROC= 0.9179), Val: 0.1744, Test: 0.2612\n",
      "Epoch: 05, AUPRC= Train: 0.2192 (AUROC= 0.9285), Val: 0.1707, Test: 0.2380\n",
      "Epoch: 06, AUPRC= Train: 0.1823 (AUROC= 0.9325), Val: 0.1762, Test: 0.2098\n",
      "Epoch: 07, AUPRC= Train: 0.2071 (AUROC= 0.9319), Val: 0.1303, Test: 0.2219\n",
      "Epoch: 08, AUPRC= Train: 0.1276 (AUROC= 0.9281), Val: 0.0886, Test: 0.1295\n",
      "Epoch: 09, AUPRC= Train: 0.2660 (AUROC= 0.9219), Val: 0.1750, Test: 0.2378\n",
      "Epoch: 10, AUPRC= Train: 0.1987 (AUROC= 0.9264), Val: 0.1597, Test: 0.2278\n",
      "Epoch: 11, AUPRC= Train: 0.2818 (AUROC= 0.9349), Val: 0.1949, Test: 0.2839\n",
      "Epoch: 12, AUPRC= Train: 0.2205 (AUROC= 0.9444), Val: 0.1715, Test: 0.2481\n",
      "Epoch: 13, AUPRC= Train: 0.2781 (AUROC= 0.9423), Val: 0.2068, Test: 0.2402\n",
      "Epoch: 14, AUPRC= Train: 0.2765 (AUROC= 0.9450), Val: 0.1755, Test: 0.2773\n",
      "Epoch: 15, AUPRC= Train: 0.2817 (AUROC= 0.9447), Val: 0.1857, Test: 0.2849\n",
      "Epoch: 16, AUPRC= Train: 0.2219 (AUROC= 0.9416), Val: 0.1363, Test: 0.2146\n",
      "Epoch: 17, AUPRC= Train: 0.2434 (AUROC= 0.9342), Val: 0.1588, Test: 0.2183\n",
      "Epoch: 18, AUPRC= Train: 0.2584 (AUROC= 0.9506), Val: 0.1694, Test: 0.2962\n",
      "Epoch: 19, AUPRC= Train: 0.1744 (AUROC= 0.9375), Val: 0.1069, Test: 0.1342\n",
      "Epoch: 20, AUPRC= Train: 0.2817 (AUROC= 0.9443), Val: 0.1667, Test: 0.2730\n",
      "Epoch: 21, AUPRC= Train: 0.2912 (AUROC= 0.9579), Val: 0.1644, Test: 0.2317\n",
      "Epoch: 22, AUPRC= Train: 0.3187 (AUROC= 0.9640), Val: 0.1949, Test: 0.3076\n",
      "Epoch: 23, AUPRC= Train: 0.3342 (AUROC= 0.9653), Val: 0.1970, Test: 0.2960\n",
      "Epoch: 24, AUPRC= Train: 0.3857 (AUROC= 0.9689), Val: 0.1951, Test: 0.2992\n",
      "Epoch: 25, AUPRC= Train: 0.3233 (AUROC= 0.9705), Val: 0.1817, Test: 0.2257\n",
      "Epoch: 26, AUPRC= Train: 0.3553 (AUROC= 0.9676), Val: 0.2093, Test: 0.2785\n",
      "Epoch: 27, AUPRC= Train: 0.3770 (AUROC= 0.9745), Val: 0.2090, Test: 0.3162\n",
      "Epoch: 28, AUPRC= Train: 0.4431 (AUROC= 0.9779), Val: 0.2158, Test: 0.3020\n",
      "Epoch: 29, AUPRC= Train: 0.4565 (AUROC= 0.9791), Val: 0.2377, Test: 0.3450\n",
      "Epoch: 30, AUPRC= Train: 0.4445 (AUROC= 0.9804), Val: 0.2161, Test: 0.3724\n",
      "Epoch: 31, AUPRC= Train: 0.4751 (AUROC= 0.9800), Val: 0.2374, Test: 0.3347\n",
      "Epoch: 32, AUPRC= Train: 0.4668 (AUROC= 0.9815), Val: 0.2278, Test: 0.3482\n",
      "Epoch: 33, AUPRC= Train: 0.4984 (AUROC= 0.9825), Val: 0.2405, Test: 0.3682\n",
      "Epoch: 34, AUPRC= Train: 0.4517 (AUROC= 0.9814), Val: 0.2378, Test: 0.3646\n",
      "Epoch: 35, AUPRC= Train: 0.4992 (AUROC= 0.9816), Val: 0.2492, Test: 0.3270\n",
      "Epoch: 36, AUPRC= Train: 0.4898 (AUROC= 0.9826), Val: 0.2615, Test: 0.3614\n",
      "Epoch: 37, AUPRC= Train: 0.4845 (AUROC= 0.9846), Val: 0.2226, Test: 0.3573\n",
      "Epoch: 38, AUPRC= Train: 0.5175 (AUROC= 0.9838), Val: 0.3051, Test: 0.3373\n",
      "Epoch: 39, AUPRC= Train: 0.4443 (AUROC= 0.9811), Val: 0.2174, Test: 0.3472\n",
      "Epoch: 40, AUPRC= Train: 0.5048 (AUROC= 0.9846), Val: 0.2537, Test: 0.3353\n",
      "Epoch: 41, AUPRC= Train: 0.5428 (AUROC= 0.9863), Val: 0.2773, Test: 0.3577\n",
      "Epoch: 42, AUPRC= Train: 0.5099 (AUROC= 0.9869), Val: 0.2485, Test: 0.3360\n",
      "Epoch: 43, AUPRC= Train: 0.4830 (AUROC= 0.9831), Val: 0.2712, Test: 0.2888\n",
      "Epoch: 44, AUPRC= Train: 0.5017 (AUROC= 0.9858), Val: 0.2175, Test: 0.2774\n",
      "Epoch: 45, AUPRC= Train: 0.5504 (AUROC= 0.9873), Val: 0.2487, Test: 0.3358\n",
      "Epoch: 46, AUPRC= Train: 0.5482 (AUROC= 0.9890), Val: 0.2463, Test: 0.3122\n",
      "Epoch: 47, AUPRC= Train: 0.5689 (AUROC= 0.9903), Val: 0.2547, Test: 0.3226\n",
      "Epoch: 48, AUPRC= Train: 0.5797 (AUROC= 0.9891), Val: 0.2694, Test: 0.3473\n",
      "Epoch: 49, AUPRC= Train: 0.5498 (AUROC= 0.9894), Val: 0.2608, Test: 0.3497\n",
      "Epoch: 50, AUPRC= Train: 0.5815 (AUROC= 0.9901), Val: 0.2650, Test: 0.3455\n",
      "\n",
      "Model's performance on the test set\n",
      "===================================\n",
      "AUROC=0.3372854818392364\n",
      "Epoch=38\n"
     ]
    }
   ],
   "source": [
    "device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
    "model = GraphTransformerNet(node_dim_in=NODE_DIM,\n",
    "                            edge_dim_in=EDGE_DIM,\n",
    "                            pe_in_dim=PE_DIM,\n",
    "                            num_gt_layers=4, \n",
    "                            hidden_dim=128,\n",
    "                            num_heads=8,\n",
    "                            norm='bn',\n",
    "                            gt_aggregators=['sum', \"mean\"],\n",
    "                            aggregators=['sum','mean','max', 'std'],\n",
    "                            dropout=0.2,\n",
    "                            act='gelu',\n",
    "                            gate=True).to(device)\n",
    "\n",
    "if int(torch.__version__.split('.')[0]) >= 2:\n",
    "    model = torch_geometric.compile(model) \n",
    "\n",
    "optimizer = torch.optim.AdamW(model.parameters(), lr=0.001)\n",
    "scheduler = ReduceLROnPlateau(optimizer, mode='max', factor=0.5, patience=5,\n",
    "                              min_lr=0.00001)\n",
    "\n",
    "print(model)\n",
    "print(f\"Number of params: {model.num_parameters()//1000} k\")\n",
    "\n",
    "\n",
    "best_epoch = 0\n",
    "best_validation_loss = -np.inf\n",
    "test_set_auroc = -np.inf\n",
    "\n",
    "best_labels, best_scores = None, None\n",
    "for epoch in range(1, 51):\n",
    "    tr_auroc = train(epoch)\n",
    "    #va_auroc = test(val_loader)\n",
    "    #te_auroc = test(test_loader)\n",
    "    \n",
    "    (tr_auprc, _, _) = calc_auprc(train_loader)\n",
    "    (va_auprc, _, _) = calc_auprc(val_loader)\n",
    "    (te_auprc, scores, labels) = calc_auprc(test_loader)\n",
    "    scheduler.step(va_auprc)\n",
    "    #print(f'Epoch: {epoch:02d}, AUROC: Train: {tr_auroc:.4f}, Val: {va_auroc:.4f}, '\n",
    "    #      f'Test: {te_auroc:.4f}')\n",
    "    print(f'Epoch: {epoch:02d}, AUPRC= Train: {tr_auprc:.4f} (AUROC= {tr_auroc:.4f}), Val: {va_auprc:.4f}, '\n",
    "          f'Test: {te_auprc:.4f}')\n",
    "        \n",
    "    va_loss = va_auprc\n",
    "    te_loss = te_auprc\n",
    "    \n",
    "    if va_loss > best_validation_loss:\n",
    "        best_epoch = epoch\n",
    "        best_validation_loss = va_loss\n",
    "        test_set_auroc = te_loss\n",
    "        \n",
    "        best_labels = labels\n",
    "        best_scores = scores\n",
    "        \n",
    "print(\"\\nModel's performance on the test set\\n\"\n",
    "        \"===================================\\n\"\n",
    "        f'AUROC={test_set_auroc}\\n'\n",
    "        f'Epoch={best_epoch}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correction: AUROC=0.3372854818392364 --> AUPRC=0.3372854818392364 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Area under Precision-Recall curve (AUPRC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the original paper, the AUPRC for the ensemble of 20 optimized MPNNs with RDKit features is approximately ~0.35-0.37; however, an exact comparison is not feasible due to the lack of specific information regarding the training, validation, and test sets. Additionally, the best RF model performs at AUPRC ~0.31, and the best ensemble of MPNNs without RDKit features at AUPRC ~0.27. As we can see in our first attempt, we already achieved performance comparable to the optimized ensemble model used in the paper, and better than RF, with AUPRC -- statistical significance tests are left as an exercise for the reader ;-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArMAAAIjCAYAAAAQgZNYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB9i0lEQVR4nO3dd3wT9f8H8FeSZnQPukuh7DLKKkM2QksZoigKsgQUHIAidaICKgoiiKAylKnfnwqCC6UCZa8is8hebSmFTqC7adb9/igNhKY7TXrwej4efZBcPnf3vnw6Xlw+9zmJIAgCiIiIiIhESGrrAoiIiIiIqophloiIiIhEi2GWiIiIiESLYZaIiIiIRIthloiIiIhEi2GWiIiIiESLYZaIiIiIRIthloiIiIhEi2GWiIiIiESLYZaIHhrjxo1DUFBQpdbZvXs3JBIJdu/eXSM1iV3v3r3Ru3dv4/OEhARIJBKsXbvWZjUR0cOFYZaIaszatWshkUiMXyqVCk2bNsWUKVOQmppq6/JqveJgWPwllUrh4eGBAQMGICYmxtblWURqairefPNNBAcHw8HBAY6OjggNDcUnn3yCzMxMW5dHRCJgZ+sCiOjB9/HHH6NBgwZQq9XYv38/li1bhqioKJw+fRoODg5Wq2PFihUwGAyVWqdnz54oKCiAQqGooarKN2LECAwcOBB6vR4XL17E0qVL8eijj+LIkSMICQmxWV3VdeTIEQwcOBC5ubkYPXo0QkNDAQBHjx7FZ599hr1792Lbtm02rpKIajuGWSKqcQMGDECHDh0AABMmTECdOnWwcOFC/PnnnxgxYoTZdfLy8uDo6GjROuRyeaXXkUqlUKlUFq2jstq3b4/Ro0cbn/fo0QMDBgzAsmXLsHTpUhtWVnWZmZl48sknIZPJcOLECQQHB5u8/umnn2LFihUW2VdNfC8RUe3BYQZEZHV9+vQBAMTHxwMoGsvq5OSEK1euYODAgXB2dsaoUaMAAAaDAYsWLULLli2hUqng4+ODl156Cbdv3y6x3X/++Qe9evWCs7MzXFxc0LFjR/z000/G182NmV23bh1CQ0ON64SEhGDx4sXG10sbM7thwwaEhobC3t4enp6eGD16NK5fv27Spvi4rl+/jiFDhsDJyQleXl548803odfrq/z+9ejRAwBw5coVk+WZmZl4/fXXERgYCKVSicaNG2PevHklzkYbDAYsXrwYISEhUKlU8PLyQv/+/XH06FFjmzVr1qBPnz7w9vaGUqlEixYtsGzZsirXfL9vv/0W169fx8KFC0sEWQDw8fHBBx98YHwukUjw4YcflmgXFBSEcePGGZ8XD23Zs2cPJk2aBG9vb9StWxcbN240LjdXi0QiwenTp43Lzp8/j6effhoeHh5QqVTo0KEDNm3aVL2DJqIawTOzRGR1xSGsTp06xmU6nQ4RERHo3r07FixYYBx+8NJLL2Ht2rUYP348XnvtNcTHx+Obb77BiRMncODAAePZ1rVr1+L5559Hy5YtMX36dLi5ueHEiRPYsmULRo4cabaO6OhojBgxAn379sW8efMAAOfOncOBAwcwderUUusvrqdjx46YO3cuUlNTsXjxYhw4cAAnTpyAm5ubsa1er0dERAQ6d+6MBQsWYPv27fjiiy/QqFEjvPLKK1V6/xISEgAA7u7uxmX5+fno1asXrl+/jpdeegn16tXDwYMHMX36dCQnJ2PRokXGti+88ALWrl2LAQMGYMKECdDpdNi3bx8OHTpkPIO+bNkytGzZEo8//jjs7Ozw119/YdKkSTAYDJg8eXKV6r7Xpk2bYG9vj6effrra2zJn0qRJ8PLywsyZM5GXl4dBgwbByckJv/zyC3r16mXSdv369WjZsiVatWoFADhz5gy6deuGgIAAvPvuu3B0dMQvv/yCIUOG4Ndff8WTTz5ZIzUTURUJREQ1ZM2aNQIAYfv27UJ6erpw7do1Yd26dUKdOnUEe3t7ISkpSRAEQRg7dqwAQHj33XdN1t+3b58AQPjxxx9Nlm/ZssVkeWZmpuDs7Cx07txZKCgoMGlrMBiMj8eOHSvUr1/f+Hzq1KmCi4uLoNPpSj2GXbt2CQCEXbt2CYIgCBqNRvD29hZatWplsq+///5bACDMnDnTZH8AhI8//thkm+3atRNCQ0NL3Wex+Ph4AYDw0UcfCenp6UJKSoqwb98+oWPHjgIAYcOGDca2s2fPFhwdHYWLFy+abOPdd98VZDKZkJiYKAiCIOzcuVMAILz22msl9nfve5Wfn1/i9YiICKFhw4Ymy3r16iX06tWrRM1r1qwp89jc3d2FNm3alNnmXgCEWbNmlVhev359YezYscbnxd9z3bt3L9GvI0aMELy9vU2WJycnC1Kp1KSP+vbtK4SEhAhqtdq4zGAwCF27dhWaNGlS4ZqJyDo4zICIalxYWBi8vLwQGBiIZ599Fk5OTvj9998REBBg0u7+M5UbNmyAq6srwsPDkZGRYfwKDQ2Fk5MTdu3aBaDoDGtOTg7efffdEuNbJRJJqXW5ubkhLy8P0dHRFT6Wo0ePIi0tDZMmTTLZ16BBgxAcHIzNmzeXWOfll182ed6jRw/ExcVVeJ+zZs2Cl5cXfH190aNHD5w7dw5ffPGFyVnNDRs2oEePHnB3dzd5r8LCwqDX67F3714AwK+//gqJRIJZs2aV2M+975W9vb3xcVZWFjIyMtCrVy/ExcUhKyurwrWXJjs7G87OztXeTmkmTpwImUxmsmz48OFIS0szGTKyceNGGAwGDB8+HABw69Yt7Ny5E8OGDUNOTo7xfbx58yYiIiJw6dKlEsNJiMi2OMyAiGrckiVL0LRpU9jZ2cHHxwfNmjWDVGr6f2k7OzvUrVvXZNmlS5eQlZUFb29vs9tNS0sDcHfYQvHHxBU1adIk/PLLLxgwYAACAgLQr18/DBs2DP379y91natXrwIAmjVrVuK14OBg7N+/32RZ8ZjUe7m7u5uM+U1PTzcZQ+vk5AQnJyfj8xdffBHPPPMM1Go1du7cia+++qrEmNtLly7hv//+K7GvYve+V/7+/vDw8Cj1GAHgwIEDmDVrFmJiYpCfn2/yWlZWFlxdXctcvzwuLi7Iycmp1jbK0qBBgxLL+vfvD1dXV6xfvx59+/YFUDTEoG3btmjatCkA4PLlyxAEATNmzMCMGTPMbjstLa3Ef8SIyHYYZomoxnXq1Mk4FrM0SqWyRMA1GAzw9vbGjz/+aHad0oJbRXl7eyM2NhZbt27FP//8g3/++Qdr1qzBc889h++//75a2y52/9lBczp27GgMyUDRmdh7L3Zq0qQJwsLCAACPPfYYZDIZ3n33XTz66KPG99VgMCA8PBxvv/222X0Uh7WKuHLlCvr27Yvg4GAsXLgQgYGBUCgUiIqKwpdfflnp6c3MCQ4ORmxsLDQaTbWmPSvtQrp7zywXUyqVGDJkCH7//XcsXboUqampOHDgAObMmWNsU3xsb775JiIiIsxuu3HjxlWul4gsj2GWiGqtRo0aYfv27ejWrZvZcHJvOwA4ffp0pYOGQqHA4MGDMXjwYBgMBkyaNAnffvstZsyYYXZb9evXBwBcuHDBOCtDsQsXLhhfr4wff/wRBQUFxucNGzYss/3777+PFStW4IMPPsCWLVsAFL0Hubm5xtBbmkaNGmHr1q24detWqWdn//rrLxQWFmLTpk2oV6+ecXnxsA5LGDx4MGJiYvDrr7+WOj3bvdzd3UvcREGj0SA5OblS+x0+fDi+//577NixA+fOnYMgCMYhBsDd914ul5f7XhJR7cAxs0RUaw0bNgx6vR6zZ88u8ZpOpzOGm379+sHZ2Rlz586FWq02aScIQqnbv3nzpslzqVSK1q1bAwAKCwvNrtOhQwd4e3tj+fLlJm3++ecfnDt3DoMGDarQsd2rW7duCAsLM36VF2bd3Nzw0ksvYevWrYiNjQVQ9F7FxMRg69atJdpnZmZCp9MBAIYOHQpBEPDRRx+VaFf8XhWfTb73vcvKysKaNWsqfWylefnll+Hn54c33ngDFy9eLPF6WloaPvnkE+PzRo0aGcf9Fvvuu+8qPcVZWFgYPDw8sH79eqxfvx6dOnUyGZLg7e2N3r1749tvvzUblNPT0yu1PyKqeTwzS0S1Vq9evfDSSy9h7ty5iI2NRb9+/SCXy3Hp0iVs2LABixcvxtNPPw0XFxd8+eWXmDBhAjp27IiRI0fC3d0dJ0+eRH5+fqlDBiZMmIBbt26hT58+qFu3Lq5evYqvv/4abdu2RfPmzc2uI5fLMW/ePIwfPx69evXCiBEjjFNzBQUFYdq0aTX5lhhNnToVixYtwmeffYZ169bhrbfewqZNm/DYY49h3LhxCA0NRV5eHk6dOoWNGzciISEBnp6eePTRRzFmzBh89dVXuHTpEvr37w+DwYB9+/bh0UcfxZQpU9CvXz/jGeuXXnoJubm5WLFiBby9vSt9JrQ07u7u+P333zFw4EC0bdvW5A5gx48fx88//4wuXboY20+YMAEvv/wyhg4divDwcJw8eRJbt26Fp6dnpfYrl8vx1FNPYd26dcjLy8OCBQtKtFmyZAm6d++OkJAQTJw4EQ0bNkRqaipiYmKQlJSEkydPVu/giciybDmVAhE92IqnSTpy5EiZ7caOHSs4OjqW+vp3330nhIaGCvb29oKzs7MQEhIivP3228KNGzdM2m3atEno2rWrYG9vL7i4uAidOnUSfv75Z5P93Ds118aNG4V+/foJ3t7egkKhEOrVqye89NJLQnJysrHN/VNzFVu/fr3Qrl07QalUCh4eHsKoUaOMU42Vd1yzZs0SKvLrt3iaq/nz55t9fdy4cYJMJhMuX74sCIIg5OTkCNOnTxcaN24sKBQKwdPTU+jatauwYMECQaPRGNfT6XTC/PnzheDgYEGhUAheXl7CgAEDhGPHjpm8l61btxZUKpUQFBQkzJs3T1i9erUAQIiPjze2q+rUXMVu3LghTJs2TWjatKmgUqkEBwcHITQ0VPj000+FrKwsYzu9Xi+88847gqenp+Dg4CBEREQIly9fLnVqrrK+56KjowUAgkQiEa5du2a2zZUrV4TnnntO8PX1FeRyuRAQECA89thjwsaNGyt0XERkPRJBKOMzOCIiIiKiWoxjZomIiIhItBhmiYiIiEi0GGaJiIiISLQYZomIiIhItBhmiYiIiEi0GGaJiIiISLQeupsmGAwG3LhxA87OzpBIJLYuh4iIiIjuIwgCcnJy4O/vD6m07HOvD12YvXHjBgIDA21dBhERERGV49q1a6hbt26ZbR66MOvs7Ayg6M1xcXGp8f1ptVps27bNeBtOEh/2ofixD8WPfShu7D/xs3YfZmdnIzAw0JjbyvLQhdnioQUuLi5WC7MODg5wcXHhD7BIsQ/Fj30ofuxDcWP/iZ+t+rAiQ0J5ARgRERERiRbDLBERERGJFsMsEREREYnWQzdmloiIqLoEQYBOp4Ner7d1KaKg1WphZ2cHtVrN90ykaqIP5XI5ZDJZtbfDMEtERFQJGo0GycnJyM/Pt3UpoiEIAnx9fXHt2jXO8S5SNdGHEokEdevWhZOTU7W2wzBLRERUQQaDAfHx8ZDJZPD394dCoWA4qwCDwYDc3Fw4OTmVOwE+1U6W7kNBEJCeno6kpCQ0adKkWmdoGWaJiIgqSKPRwGAwIDAwEA4ODrYuRzQMBgM0Gg1UKhXDrEjVRB96eXkhISEBWq22WmGW31FERESVxEBGVH2W+lSDP41EREREJFoMs0REREQkWgyzRERERCRaDLNEREQPiZiYGMhkMgwaNKjEa7t374ZEIkFmZmaJ14KCgrBo0SKTZbt27cLAgQNRp04dODg4oEWLFnjjjTdw/fr1Uve/du1aSCQSSCQSSKVS1K1bF+PHj0daWpqxTfHrEokELi4u6NixI/78888qH3NF3bp1C6NGjYKLiwvc3NzwwgsvIDc3t8x1XnrpJTRq1Aj29vbw8vLCE088gfPnz5tte/PmTdStW7fU95iqjmGWiIjoIbFq1Sq8+uqr2Lt3L27cuFHl7Xz77bcICwuDr68vfv31V5w9exbLly9HVlYWvvjiizLXdXFxQXJyMpKSkrBixQr8888/GDNmjEmbNWvWIDk5GUePHkW3bt3w9NNP49SpU1WutyJGjRqFM2fOIDo6Gn///Tf27t2LF198scx1QkNDsWbNGpw7dw5bt26FIAjo16+f2ZsKvPDCC2jdunVNlf9Q49RcREREVSQIAgq0trmjlb1cVqmrwXNzc7F+/XocPXoUKSkpWLt2Ld57771K7zcpKQmvvfYaXnvtNXz55ZfG5UFBQejZs2e5Zx0lEgl8fX0BAP7+/njttdcwY8YMFBQUwN7eHgDg5uYGX19f+Pr6Yvbs2Vi8eDF27dqFkJCQStdbEefOncOWLVtw5MgRdOjQAQDw9ddfY+DAgViwYAH8/f3Nrndv2A0KCsInn3yCNm3aICEhAY0aNTK+tmzZMmRmZmLmzJn4559/auQYHmY2DbN79+7F/PnzcezYMSQnJ+P333/HkCFDylxn9+7diIyMxJkzZxAYGIgPPvgA48aNs0q9RERE9yrQ6tFi5lab7PvsxxFwUFT8z/gvv/yC4OBgNGvWDKNHj8brr7+O6dOnV3p6pA0bNkCj0eDtt982+7qbm1ultmdvbw+DwQCdTlfiNZ1Oh1WrVgEAFApFmdtp2bIlrl69WurrPXr0KDVIxsTEwM3NzRhkASAsLAxSqRT//vsvnnzyyXKPIy8vD2vWrEGDBg0QGBhoXH727Fl8/PHH+PfffxEXF1fudqjybBpm8/Ly0KZNGzz//PN46qmnym0fHx+PQYMG4eWXX8aPP/6IHTt2YMKECfDz80NERIQVKiYiIhKnVatWYfTo0QCA/v37IysrC3v27EHv3r0rtZ1Lly7BxcUFfn5+1a7p0qVLWL58OTp06ABnZ2fj8hEjRkAmk6GgoAAGgwFBQUEYNmxYmduKioqCVqst9fXis77mpKSkwNvb22SZnZ0dPDw8kJKSUuZ+ly5dirfffht5eXlo1qwZoqOjjcG7sLAQI0aMwPz581GvXj2G2Rpi0zA7YMAADBgwoMLtly9fjgYNGhjH4zRv3hz79+/Hl19+WWvD7NnkbJy8KYHsTCrs7Mq/u4W3iwrtAt14e0QiIhGwl8tw9mPb/P2xl1f8jkkXLlzA4cOH8fvvvwMoCmrDhw/HqlWrKh1mBUGo0N8oJycn4+NRo0Zh3rx5AICsrCw4OTnBYDBArVaje/fuWLlypcm6X375JcLCwhAXF4dp06bhq6++goeHR5n7q1+/fqWOw1JGjRqF8PBwJCcnY8GCBRg2bBgOHDgAlUqF6dOno3nz5sb/RFDNENWY2ZiYGISFhZksi4iIwOuvv17qOoWFhSgsLDQ+z87OBgBotdoy/wdnKesOX8PPF2VYffFkhdf5eUJHdKjvXoNVUWUUf59Y4/uFagb7UPxqSx9qtVoIggCDwQCDwQAAUNnZ5lpqQRAgCEKF2q5cuRI6nc5k7KcgCFAqlfjqq6/g6upqDJ+3b9+Gi4uLyfqZmZlwdnaGwWBAkyZNkJWVhevXr5d5dvb48ePGx/eedXV2dsbRo0chlUrh5+dnPGNa/H4CgLe3Nxo2bIiGDRti1apVeOyxx3D69OkSZ0/vFRISUuYwg+7duyMqKsrsa97e3khLSzOpQafT4datW/D29jZZfj9nZ2c4OzujUaNG6NSpE+rUqYNff/0VI0aMwM6dO3Hq1Cls3LgRAIz95enpiffeew8ffvhhqdutbYprL/7+twSDwQBBEMzezrYyP+uiCrMpKSnw8fExWebj44Ps7GyTgeP3mjt3Lj766KMSy7dt22aV+2rnp0vQwLliv+iS8wG1XoItew4hzbNiv6DIeqKjo21dAlUT+1D8bN2HdnZ28PX1RW5uLjQajU1rqSidTocffvgBn3zyCR599FGT10aPHo01a9bg+eefh4+PD6RSKfbv3w9397snVBISEpCVlYWAgABkZ2ejX79+UCgU+PTTTzFnzpwS+8vKyoKrq6vZ4KlWqyGRSIyvlXZiqaCgwHjyKTg4GG3atMGHH36Izz77rNTj/Pnnn82Ouy2mUqmM27xfSEgIMjMzsXfvXrRt2xYAsHPnThgMBrRo0aLU9e5XWFgIQRCQlZWF7OxsrFmzBgUFBcbXT5w4gSlTpiAqKgoNGjSo8HZrk5ycHIttS6PRoKCgAHv37i3Rd/n5+RXejqjCbFVMnz4dkZGRxufZ2dkIDAxEv379SvzPsyaEa7WIjo5GeHg45HJ5mW3HrD6CQ/G30a5dOwwM8a3x2qhitJXoQ6qd2IfiV1v6UK1W49q1a3BycoJKpbJZHZXxxx9/IDMzE5MmTYKrq6vJa08//TR+/vlnvP7663BxccELL7yAmTNnwtnZGSEhIbh27RqmT5+ORx55BOHh4ZBIJGjRogUWLlyIV199FWq1GmPGjEFQUBCSkpLwv//9D05OTliwYIHJfgRBQE5ODlQqlXH+2LLY29ubtImMjMTQoUPx/vvvIyAgwOw6rVq1quI7BHTs2BERERGIjIzE0qVLodVq8e6772L48OFo1qwZAOD69esIDw/H2rVr0alTJ8TFxeGXX35BeHg4vLy8kJSUhHnz5sHe3h5PPfUUXFxc0KZNG5P9qNVqAECHDh0qfaGcrRX3obOzs8WGQqrVatjb26Nnz54lfp4qE/RFFWZ9fX2Rmppqsiw1NRUuLi6lDuxWKpVQKpUllsvlcqv+QqzI/iSSojO4MpmMf3BrIWt/z5DlsQ/Fz9Z9qNfrjRP+S6XimKp9zZo1CAsLMznbWuzpp5/G/Pnzcfr0abRu3RpfffUVPvvsM0yfPh1Xr16Fr68vwsPD8emnn5p8DDx58mQ0a9YMCxYswNChQ1FQUICgoCA89thjiIyMLPHe3P+xdHnv3f3v78CBA9GgQQPMnTsXS5curcrbUK6ffvoJU6ZMQXh4OKRSKYYOHYqvvvrKWIder8eFCxegVqshlUrh4OCA/fv3Y/Hixbh9+zZ8fHzQs2dPHDx40Dj1mLnjMnd8YlDch8Xf/5YglUohkUjM/lxX5udcVGG2S5cuJca7REdHo0uXLjaqiIiIqHb766+/Sn2tU6dOJuNuVSoVPvzwwwqN5QwLCytxHUt5xo0bh+eff77MNubGAUskEpw7d65S+6osDw8P/PTTT6W+HhQUZFKbv79/qWNwS9O7d+8Kj3OmirPpfwtyc3MRGxuL2NhYAEVTb8XGxiIxMRFA0RCB5557ztj+5ZdfRlxcHN5++22cP38eS5cuxS+//IJp06bZonwiIiIisjGbhtmjR4+iXbt2aNeuHYCiMTHt2rXDzJkzAQDJycnGYAsADRo0wObNmxEdHY02bdrgiy++wMqVK2vttFxEREREVLNsOsygvNPta9euNbvOiRMnarAqIiIiIhILcY0+JiIiIiK6B8MsEREREYkWwywRERERiRbDLBERERGJFsMsEREREYkWwywRERERiRbDLBEREVnc7t27IZFIkJmZafb1hIQESCQS442T7m+/du1auLm5lbufVatWoV+/fpYpmkSJYZaIiIisLjAwEMnJyWjVqlWVt6FWqzFjxgzMmjWrxGtJSUlQKBRmt39/kL5X79698frrr5ssO3HiBJ555hn4+PhApVKhSZMmmDhxIi5evFhqbcXhvPjLx8cHQ4cORVxcnLFNUFCQ8XUHBweEhIRg5cqVFX8DqkitVmPy5MmoU6cOnJycMHToUKSmppa5jiAImDNnDgICAmBvb4+wsDBcunTJpM2nn36Krl27wsHBoUL/EbEUhlkiIiKyOplMBl9fX9jZVf3+TRs3boSLiwu6detW4rW1a9di2LBhyM7Oxr///lvlffz999945JFHUFhYiB9//BHnzp3D//3f/8HV1RUzZswod/0LFy7gxo0b2LBhA86cOYPBgwdDr9cbX//444+RnJyM06dPY/To0Zg4cSL++eefKtdbEdOmTcNff/2FDRs2YM+ePbhx4waeeuqpMteZP38+vv32WyxduhT//vsvHB0dERERAbVabWyj0WjwzDPP4JVXXqnR+u/HMEtERFRdmrzSv7TqSrQtqFjbSurduzdee+01vP322/Dw8ICvry8+/PBDkzaJiYl44okn4OTkBBcXFwwbNqzcs3VXr17F4MGD4e7uDkdHR7Rs2RJRUVFm2+bn52PgwIHo1q0bMjMzyzw7WlHr1q3D4MGDSywXBAFr1qzBmDFjMHLkSKxatapK28/Pz8f48eMxcOBAbNq0CWFhYWjQoAE6d+6MBQsW4Ntvvy13G97e3vDz80PPnj0xc+ZMnD17FpcvXza+7uzsDF9fXzRs2BDvvPMOPDw8EB0dXaV6KyIrKwurVq3CwoUL0adPH4SGhmLNmjU4ePAgDh06ZHYdQRCwePFivPnmm3jiiSfQunVr/PDDD7hx4wb++OMPY7uPPvoI06ZNQ0hISI3Vb45Nb2dLRET0QJjjX/prTfoBozbcfT6/MaDNN9+2fndg/Oa7zxeFAPk3S7b7MKvSJX7//feIjIzEv//+i5iYGIwbNw7dunVDeHg4DAaDMcju2bMHOp0OkydPxvDhw7F79+5Stzl58mRoNBrs3bsXjo6OOHv2LJycnEq0y8zMxJNPPglXV1dER0fDwcGh1LG0lbF//36MGTOmxPJdu3YhPz8fYWFhCAgIQNeuXfHll1/C0dGxUtvfunUrMjIy8Pbbb5t9vbIfpdvb2wMoOoN5P4PBgN9//x23b9+GQqEoczsDBgzAvn37Sn29fv36OHPmjNnXjh07Bq1Wi7CwMOOy4OBg1KtXDzExMXjkkUdKrBMfH4+UlBT07t3buMzV1RWdO3dGTEwMnn322TLrrWkMs0RERA+B1q1bG8eWNmnSBN988w127NiB8PBw7NixA6dOnUJ8fDwCAwMBAD/88ANatmyJI0eOoGPHjma3mZiYiKFDhxrPxDVs2LBEm5SUFAwfPhxBQUFYv349VCqVRY4nMzMTWVlZ8Pcv+R+JVatW4dlnn4VMJkOrVq3QsGFDbNiwAePGjavUPorHhAYHB1e73uTkZCxYsAABAQFo1qyZcfk777yDDz74AIWFhdDpdPDw8MCECRPK3NbKlStRUFBQ6utyubzU11JSUqBQKEoEcR8fH6SkpJS6DgB4eXlVeB1rYpglIiKqrvdulP6aRGb6/K3L5tsBgOS+0X+vn6p6Tfdp3bq1yXM/Pz+kpaUBAM6dO4fAwEBjkAWAFi1awM3NDefOnUPHjh3RsmVLXL16FQDQo0cP/PPPP3jttdfwyiuvYNu2bQgLC8PQoUNL7Cc8PBwdO3bEihUryj3jWBnFYe7+cJyZmYnffvsN+/fvNy4bPXo0Vq1aVekwKwhChdqZe2+K1a1bF4IgID8/H23atMGvv/5q8j689dZbGDduHJKTk/HWW29h0qRJaNy4cZn7CwgIqNRxPOgYZomIiKpLUYmPr2uqbTnuP1snkUhgMBgqvH5UVBS0Wi2Aux+XT5gwAREREdi8eTO2bduGuXPn4osvvsCrr75qXG/QoEH49ddfcf78eXTp0sUCR1KkTp06kEgkuH37tsnyn376CWq1Gp07dzYuEwQBBoMBFy9eRNOmTeHi4gKgaPzo/TIzM+Hq6goAaNq0KQCUW7u596bYvn374OLiAm9vbzg7O5dY19PTE40bN0bjxo2xYcMGhISEoEOHDmjRokWp+6vOMANfX19oNBpkZmaanJ1NTU2Fr69vqesAQHp6uvE9KV6nbdu2pdZhLbwAjIiI6CHXvHlzXLt2DdeuXTMuO3v2LDIzM42hqn79+sbQde+ZwcDAQLz88sv47bff8MYbb2DFihUm2/7ss8/w3HPPYciQITh79qzFalYoFGjRokWJba5atQpvvPEGYmNjjV8nT55Ejx49sHr1agCAh4cHPD09cezYMZN1s7OzcfnyZWNg69evHzw9PfH555+braF43G9p7w0ANGjQAI0aNTIbZO8XGBiI4cOHY/r06WW2W7lypcnx3f9V2kV4ABAaGgq5XI4dO3YYl124cAGJiYmlBvYGDRrA19cXe/bsMS4rniXCkv9BqSqemSUiInrIhYWFISQkBKNGjcKiRYug0+kwadIk9OrVCx06dCh1vddffx0DBgxA06ZNcfv2bezatQvNmzcv0W7+/PkoKChAWFgYdu/ebZExqAAQERGB/fv3G+eFjY2NxfHjx/Hjjz+W2MeIESPw8ccf45NPPoGdnR0iIyMxZ84c+Pj44JFHHsHNmzcxe/ZseHl5GaepcnR0xMqVK/HMM8/g8ccfx2uvvYbGjRsjIyMDv/zyCxITE7Fu3TqLHEuxqVOnolWrVjh69Gip7311hhm4urrihRdeQGRkJDw8PODi4oJXX30VXbp0Mbn4Kzg4GHPnzsWTTz4JiUSCqVOn4rPPPkOrVq3QqFEjzJgxA/7+/hgyZIhxncTERNy6dQuJiYnQ6/XGmSoaN25s9sJAS2GYJSIieshJJBL8+eefePXVV9GzZ09IpVL0798fX3/9dZnr6fV6TJ48GUlJSXBxcUH//v3x5Zdfmm07Z84cyGQy9OnTB7t377bI+NkXXngBHTp0QFZWFlxdXbFq1Sq0aNHCbFh+8sknMWXKFERFReHxxx/H22+/DScnJ8ybNw9XrlyBh4cHunXrhl27dpkMFXjiiSdw8OBBzJ07FyNHjkR2djYCAwPRp08ffPLJJ9U+hvu1aNEC/fr1w8yZM8s8w1odX375JaRSKYYOHYrCwkJERERg6dKlJm0uXLhgMgzjrbfewq1bt/Dyyy8jMzMT3bt3x5YtW0zGLM+cORPff/+98Xm7du0AFM0uce9MCJYmESo6uvkBkZ2dDVdXV2RlZRnHzNQkrVaLqKgoDBw4sMyrCwFgxHeHEBN3E1+PaIfBbcqY5oWsqjJ9SLUT+1D8aksfqtVqxMfHo0GDBha7Kv9hYDAYkJ2dDRcXF0illh3h+Mwzz6B9+/blfjRP1VMTfVjWz1Nl8hrHzBIREZFozZ8/v0Y/wqbaj2GWiIiIRCsoKMhk9gR6+DDMEhEREZFoMcwSERERkWgxzBIREVXSQ3btNFGNsNTPEcMsERFRBRXPpJCfn2/jSojET6PRAABkMlk5LcvGeWaJiIgqSCaTwc3NDWlpaQAABwcHSCQSG1dV+xkMBmg0GqjVaotPzUXWYek+NBgMSE9Ph4ODA+zsqhdHGWaJiIgqofg+9cWBlsonCAIKCgpgb2/P8C9SNdGHUqkU9erVq/b2GGaJiIgqQSKRwM/PD97e3tBqtbYuRxS0Wi327t2Lnj178sYlIlUTfahQKCxylpdhloiIqApkMlm1x/o9LGQyGXQ6HVQqFcOsSNXmPuTAFSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLZuH2SVLliAoKAgqlQqdO3fG4cOHy2y/aNEiNGvWDPb29ggMDMS0adOgVqutVC0RERER1SY2DbPr169HZGQkZs2ahePHj6NNmzaIiIhAWlqa2fY//fQT3n33XcyaNQvnzp3DqlWrsH79erz33ntWrpyIiIiIagObhtmFCxdi4sSJGD9+PFq0aIHly5fDwcEBq1evNtv+4MGD6NatG0aOHImgoCD069cPI0aMKPdsLhERERE9mOxstWONRoNjx45h+vTpxmVSqRRhYWGIiYkxu07Xrl3xf//3fzh8+DA6deqEuLg4REVFYcyYMaXup7CwEIWFhcbn2dnZAACtVgutVmuhoyld8T4qsi9BMAAA9Hq9VWqjiqlMH1LtxD4UP/ahuLH/xM/afViZ/dgszGZkZECv18PHx8dkuY+PD86fP292nZEjRyIjIwPdu3eHIAjQ6XR4+eWXyxxmMHfuXHz00Ucllm/btg0ODg7VO4hKiI6OLrfNzZtSAFKcOHECkmtCzRdFlVKRPqTajX0ofuxDcWP/iZ+1+jA/P7/CbW0WZqti9+7dmDNnDpYuXYrOnTvj8uXLmDp1KmbPno0ZM2aYXWf69OmIjIw0Ps/OzkZgYCD69esHFxeXGq9Zq9UiOjoa4eHhkMvlZbb9OeUILmXfRrt27TAwxLfGa6OKqUwfUu3EPhQ/9qG4sf/Ez9p9WPxJekXYLMx6enpCJpMhNTXVZHlqaip8fc0HuRkzZmDMmDGYMGECACAkJAR5eXl48cUX8f7770MqLTkEWKlUQqlUllgul8ut+gNVkf1JJEX1y2Qy/rDXQtb+niHLYx+KH/tQ3Nh/4metPqzMPmx2AZhCoUBoaCh27NhhXGYwGLBjxw506dLF7Dr5+fklAqtMJgMACAI/liciIiJ62Nh0mEFkZCTGjh2LDh06oFOnTli0aBHy8vIwfvx4AMBzzz2HgIAAzJ07FwAwePBgLFy4EO3atTMOM5gxYwYGDx5sDLVERERE9PCwaZgdPnw40tPTMXPmTKSkpKBt27bYsmWL8aKwxMREkzOxH3zwASQSCT744ANcv34dXl5eGDx4MD799FNbHQIRERER2ZDNLwCbMmUKpkyZYva13bt3mzy3s7PDrFmzMGvWLCtURkRERES1nc1vZ0tEREREVFUMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWgyzRERERCRaDLNEREREJFoMs0REREQkWjYPs0uWLEFQUBBUKhU6d+6Mw4cPl9k+MzMTkydPhp+fH5RKJZo2bYqoqCgrVUtEREREtYmdLXe+fv16REZGYvny5ejcuTMWLVqEiIgIXLhwAd7e3iXaazQahIeHw9vbGxs3bkRAQACuXr0KNzc36xdPRERERDZn0zC7cOFCTJw4EePHjwcALF++HJs3b8bq1avx7rvvlmi/evVq3Lp1CwcPHoRcLgcABAUFWbNkIiIiIqpFbBZmNRoNjh07hunTpxuXSaVShIWFISYmxuw6mzZtQpcuXTB58mT8+eef8PLywsiRI/HOO+9AJpOZXaewsBCFhYXG59nZ2QAArVYLrVZrwSMyr3gfFdmXIBgAAHq93iq1UcVUpg+pdmIfih/7UNzYf+Jn7T6szH5sFmYzMjKg1+vh4+NjstzHxwfnz583u05cXBx27tyJUaNGISoqCpcvX8akSZOg1Woxa9Yss+vMnTsXH330UYnl27Ztg4ODQ/UPpIKio6PLbXPzphSAFCdOnIDkmlDzRVGlVKQPqXZjH4of+1Dc2H/iZ60+zM/Pr3Bbmw4zqCyDwQBvb2989913kMlkCA0NxfXr1zF//vxSw+z06dMRGRlpfJ6dnY3AwED069cPLi4uNV6zVqtFdHQ0wsPDjUMjSvNzyhFcyr6Ndu3aYWCIb43XRhVTmT6k2ol9KH7sQ3Fj/4mftfuw+JP0irBZmPX09IRMJkNqaqrJ8tTUVPj6mg9yfn5+kMvlJkMKmjdvjpSUFGg0GigUihLrKJVKKJXKEsvlcrlVf6Aqsj+JpGhyCZlMxh/2Wsja3zNkeexD8WMfihv7T/ys1YeV2YfNpuZSKBQIDQ3Fjh07jMsMBgN27NiBLl26mF2nW7duuHz5MgwGg3HZxYsX4efnZzbIEhEREdGDzabzzEZGRmLFihX4/vvvce7cObzyyivIy8szzm7w3HPPmVwg9sorr+DWrVuYOnUqLl68iM2bN2POnDmYPHmyrQ6BiIiIiGzIpmNmhw8fjvT0dMycORMpKSlo27YttmzZYrwoLDExEVLp3bwdGBiIrVu3Ytq0aWjdujUCAgIwdepUvPPOO7Y6BCIiIiKyIZtfADZlyhRMmTLF7Gu7d+8usaxLly44dOhQDVdFRERERGJg89vZEhERERFVFcMsEREREYkWwywRERERiRbDLBERERGJFsMsGWWrtXhu9WGs2Btn61KIiIiIKoRhlox+P34dey+m48d/r9q6FCIiIqIKqdLUXHq9HmvXrsWOHTuQlpZmckcuANi5c6dFiiPr+iP2OgBAsHEdRERERBVVpTA7depUrF27FoMGDUKrVq0gkUgsXRdZWUJGHk4kZtq6jBqj0xsgkUggk/J7lYiI6EFSpTC7bt06/PLLLxg4cKCl6yEb+TP2hq1LqDHJWQUYueJfKGRS/DO1B6QMtERERA+MKoVZhUKBxo0bW7oWshFBEPDnnSEGD5octRbj1xxBfEYeACBPo4OzSm7jqoiIiMhSqnQB2BtvvIHFixdDEDi68kFw6noW4u6EvQeJVm/ApB+P43xKjq1LISIiohpSpTOz+/fvx65du/DPP/+gZcuWkMtNz3T99ttvFimOrOP3E0VnZQPc7HE9s8DG1ViGIAh4//dT2HcpA0o7KQp1hvJXIiIiItGp0plZNzc3PPnkk+jVqxc8PT3h6upq8kXiodMb8NfJZADAkHb+Nq7Gcr7ZeRm/HE2CVAIsGt7W1uUQERFRDanSmdk1a9ZYug6ykYNXbiIjtxAejgr0aOKFJbuu2LqkavvteBK+iL4IAPj4iVZ4NNjbxhURERFRTalSmC2Wnp6OCxcuAACaNWsGLy8vixRF1vPHnSEGg0L8IJeJ/yr/g5cz8M6v/wEAXurVEKMfqQ+1Vm/jqoiIiKimVGmYQV5eHp5//nn4+fmhZ8+e6NmzJ/z9/fHCCy8gPz/f0jVSDSnQ6LH1TAoAYEi7ABtXU30XU3Pw0v8dg1Yv4LHWfngnItjWJREREVENq1KYjYyMxJ49e/DXX38hMzMTmZmZ+PPPP7Fnzx688cYblq6Rakj0uVTkafSo5+GA9vXcbF1OtaRmqzFu9WHkqHXoGOSOBc+04XyyRERED4EqDTP49ddfsXHjRvTu3du4bODAgbC3t8ewYcOwbNkyS9VHNejPO0MMnmjrL+q7uOUV6vD82iO4kaVGQy9HrHiuA1Ryma3LIiIiIiuo0pnZ/Px8+Pj4lFju7e3NYQYicStPgz0X0wEAT7QV7xADnd6AyT8dx5kb2ajjqMDacZ3g5qCwdVlERERkJVUKs126dMGsWbOgVquNywoKCvDRRx+hS5cuFiuOas7m/25AZxAQEuCKxt5Oti6nSgRBwIw/z2D3hXSo5FKsGtcR9eo42LosIiIisqIqDTNYvHgxIiIiULduXbRp0wYAcPLkSahUKmzdutWiBVLN+CP2BoCiIQZitXxPHH4+nAiJBPjq2XZoG+hm65KIiIjIyqoUZlu1aoVLly7hxx9/xPnz5wEAI0aMwKhRo2Bvb2/RAsnyEm/m49jV25BKgMfbiDPM7r+Ugc+3Fn3vfTi4Jfq19LVxRURERGQLVZ5n1sHBARMnTrRkLWQlf8YWXfjVrbEnvF1UNq6m8lKy1Ji67gQEAXi2YyDGdg2ydUlERERkIxUOs5s2bcKAAQMgl8uxadOmMts+/vjj1S6MaoYgCPgjtngWA/Fd+KXVGzDlp+O4madBCz8XfPh4S1uXRERERDZU4TA7ZMgQpKSkwNvbG0OGDCm1nUQigV7POy7VVmduZONKeh6UdlJEtCw5I0Vt9/mW8zh69TacVXZYNro9p+AiIiJ6yFU4zBoMBrOPSVx+vzO3bFgLHzir5DaupnK2nE7Bin3xAID5T7dB/TqONq6IiIiIbK1KU3OZk5mZaalNUQ3RGwT8dbJoFoMnRTbE4OrNPLy14SQAYGKPBujfihd8ERERURXD7Lx587B+/Xrj82eeeQYeHh4ICAjAyZMnLVYcWVbMlZtIyymEm4McPZt62bqcClNr9Xjl/44jp1CHDvXd8Xb/YFuXRERERLVElcLs8uXLERgYCACIjo7G9u3bsWXLFgwYMABvvfWWRQskyym+8GtQiB8UdhY7KV/jPvrrDM4mF93h65uR7SGXiad2IiIiqllVmporJSXFGGb//vtvDBs2DP369UNQUBA6d+5s0QLJMtRaPbacTgEADGknniEGvx5Lws+Hr0EiARY/2w6+ruKbSoyIiIhqTpVOcbm7u+PatWsAgC1btiAsLAxA0bRPnMmgdtp+LhW5hTrUdbdHaD13W5dTIedTsvH+H6cAAK/3bYruTTxtXBERERHVNlU6M/vUU09h5MiRaNKkCW7evIkBAwYAAE6cOIHGjRtbtECyjD9O3L19rVQqsXE15cst1GHS/x2HWmtAz6ZeeLUPv6+IiIiopCqF2S+//BJBQUG4du0aPv/8czg5OQEAkpOTMWnSJIsWSNV3O0+D3RfSAABDRDCLgSAIeOfX/xCXkQc/VxUWDW8rigBORERE1lelMCuXy/Hmm2+WWD5t2rRqF0SWt/lUMnQGAS39XdDEx9nW5ZTrh5ir2PxfMuykEnwzsj08HBW2LomIiIhqKd7OVqROXstE4q18DG7jX27bP+/MYmDJs7I6vQF/xN5A10Z14O9mb7Htnki8jU82nwUATB/YHKH1rTu+VxAEbDp5Aw09nRBS17XMtmqtHn+cuI6ujTxRr46DlSokIiKie/F2tiL1yv8dw40sNdoGuiHQo/QglZajxpGE25BIUKHgW1Hf7o3D/K0X8FT7ACwc1tYi2yzQ6PH6+lho9QIGtPLF892CLLLdihIEAR/9dRZrDyagsbcTtkf2KrVtbqEOE74/gkNxt9C/pS+Wjwm1YqVERERUjLezFaGsfC1uZKkBFIWqshyKuwUAaOHnYrFprdRaPdYcSAAAZBdoLbJNAFiw7QKu3syHn6sK855uDYnEeuNkBUHAZ1vOY+3BBABAjrr047qVp8G4NYfxX1IWgPL7gIiIiGpOlcbMkm1dycitcNuYKxkAgC4N61hs/7+fuI6M3EKLbQ8AjibcwuoD8QCAOU+FwEUlt+j2y7No+yV8uyeu3HYpWWqMWfUvLqVVvA+IiIio5lRpntnXXnsNX331VYnl33zzDV5//fXq1kTliEvPq3DbmCs3AQBdGlkmzBoMAlbsLT/0VYZaq8fbG/+DIABPh9bFo828Lbr98izZdRmLd1wCADwTWrfUdldv5ePp5QdxKS0Xfq4qvMbpwoiIiGyuSmH2119/Rbdu3Uos79q1KzZu3FjtoqhscekVOyuYnFWAhJv5kEkl6NTAwyL7jj6XiriMiofpilgYfRFxGXnwcVFixqAWFt12eVbtj8f8rRcAAO8OCMbYrkFm293IA0asOIyk2wUIquOADS93QSNvJytWSkREROZUKczevHkTrq4lr/R2cXFBRkZGtYuislX0zGzxWdlWAa5wttDH9t/dOSvr62KZ8bfHE29j5b6ibc55MgSuDtYbXvC/Q1cx+++imRNeD2uCl3s1MtvuxLVMfH1GhvRcDYJ9nfHLy11Q152zFxAREdUGVQqzjRs3xpYtW0os/+eff9CwYcNqF0Vli6vgmFnjEAMLjZc9mnALx67ehkImxajO9aq9PbVWj7c2nIRBAJ5qF4C+zX0sUGXF/HLkGmb8cRoA8ErvRpjat4nZdvsvZWDc2mPI10vQLtAV61/sAm9nywR5IiIiqr4qXQAWGRmJKVOmID09HX369AEA7NixA1988QUWLVpkyfroPnqDgISM/Aq1PWjh8bLf3jkr+2S7AHi7KKu9vcU7LuFKeh68nJWYOdh6wwv+jL2Od377DwDwfLcGeDuimdmZE7acTsFrP5+ARm9AsKsBa8eFWvXMMREREZWvSmH2+eefR2FhIT799FPMnj0bABAUFIRly5bhueees2iBZCrpdj40+vKnRrt2Kx/XMwtgJ5WggwVuPHAlPRfbz6UCACb2bIhjV29Va3snr2Xi2z1XAACfDmkFNwfr3OXrn9MpmP7bKQgCMKpzPcx4rLnZIHszV4NJPx6DQQAiWnijn/MNOCg4+QcREVFtU6VhBgDwyiuvICkpCampqcjOzkZcXByDrBVUdrxs20A3OCqrH8JW7ouDIABhzX3QuJoXPhXq9HhrY9Hwgsfb+KNfS99q11dR7/76H/QGAU+H1sXsJ1qVOpetziDAIADDOtTFomGtYVflnxQiIiKqSVX+E63T6bB9+3b89ttvEAQBAHDjxg3k5nL+zZp0pYIzGcTEWW6IQVqOGr8eK7ol7ku9qj8m+usdl3ExNReeTgp8+HjLam+vMooD9LyhrSGVlgyyctndH4kJ3Rtg3tDWsJMxyRIREdVWVTpld/XqVfTv3x+JiYkoLCxEeHg4nJ2dMW/ePBQWFmL58uWWrpPuqMi0WIIgWPTir+8PJkCjN6B9PbdqD1k4fT0Ly+4ML5j9RCt4ONb88IJ7T772b+mLhcPaQGYmyAJAE28nTOzRAEGejhjZqZ5V70JGRERElVelMDt16lR06NABJ0+eRJ06d8PSk08+iYkTJ1qsOCqpInPMxmfkISVbDYVMivbVDJ95hTr8L+YqAODFno2qFe40OgPe3HASeoOAQa39MCDEr1q1VZTSTobXw5ogu0CHdwcEl3mmVSqV4H0rz3VLREREVVelMLtv3z4cPHgQCoXpWbWgoCBcv37dIoWReRUZM1s8xKBdPTeo5LJq7W/dkWvIVuvQwNMR4S2qN3XWkl2XcT4lBx6OCnxs5eEFr4c1ter+iIiIyDqqNBjQYDBAr9eXWJ6UlARnZ+dqF0Xm5ai1SMspBACo5KV3naVuYavVG7B6fzwAYGKPhqV+NF8RZ25kYcmuywCAj59oiTpO1Z/ai4iIiKhKYbZfv34m88lKJBLk5uZi1qxZGDhwoKVqo/sUn5X1dFKWekcvQRBwKK5o2qyujTyrtb+oU8m4nlkATycFnmofUOXtaPUGvLXhP+gMAvq39MUgKw0vIKqNsgq0xotmiYio+qoUZhcsWIADBw6gRYsWUKvVGDlypHGIwbx58yxdI91RfOevhl6Opba5nJaLjNxCqORStAksecvhihIEAd/uKbpJwriuQdUarrBiXxzOJmfDzUGO2UNKnw6LSIwy8zVYfyQRl9NySm1TqNPjz9jreHrZQbT5aJvxBiRERFR9VRozGxgYiJMnT2L9+vU4efIkcnNz8cILL2DUqFGwt7e3dI10R/GZ2UZejogvZVaD4rt+dajvAaVd1QPo/ssZOJucDQeFDKMfqV/l7aTlqLFkZ9HwghmDWsDLmcML6MFwMTUHaw4k4PcTSVBrDXikoQfWvdjFpE3S7Xz89G8i1h+5hpt5GuPy88nZ1i6XiOiBVekwq9VqERwcjL///hujRo3CqFGjaqIuMqM4zDb0dAKQZraNpcbLfnfnzNGwDoHVujvXF1svIk+jR5tANzzZrupDFYiqQ6s34Ei6BMuWxEAikWDTlG4mcwpXlMEgYOf5NKw5GI8Dl2+avJZbqDO22XspHf936Cp2nE9D8YgCXxcV6rrb4+jV29U+HiIiuqvSYVYul0OtVtdELVSO4hsmlDbMwGAQcCi+6A/sI9WYXzYlS42rN/Mhk0rwQvcGVd7O6etZ+OXYNQDAzMdamL1JAVFNUmv12HD0Gr7dcwVJmTIARUMBUrLUCPRwqPB2ctRabDiahO9jEnD1Zj4AQCoB+rXwRTNfZyzecQn5Gj2+23sF/3coEYm38o3rdm/sidGP1EdYc2+sPZjAMEtEZGFVGmYwefJkzJs3DytXroSdHe9Xbw0Gg4CEm8XDDMzfTvZ8Sg4y87VwUMjQum7Vx8sW6gwAgEEhfpX6g38vQRDwyeazEARgcBt/hFZzvluiyshRa/F/hxKxan88MnKLZgBxkgvI10lgqMS1V/EZefj+YAI2HL2GPE3RDC4uKjs826kexjxSH4EeDth9oehTkrj0PMyJOg8AcFbZ4enQuhj9SP1Sf16JiMgyqpREjxw5gh07dmDbtm0ICQmBo6PpmcLffvvNIsXRXTeyCqDWGiCXSVDX3fy45OL5ZTs18KjSR6j3e7Fn1W9du+1sKg7F3YLSTop3+jerdi0Pk6s383Ao7iaeaBtQ7XmCrUGt1UNpJ60VF/Zl5BZizYF4/BBzFTnqoo/9A9zsMaF7fTilncaMEwqotYYytyEIAmLibmLVvnjsvHB3mEBjbyeM6xqEp9oHwEFx91fnvTOLtPBzwXNd6uPxtv4mbYiIqOZU6betm5sbhg4daulaqAzF42Xr13Es9Q5WMVcyAFjmFrbdG3uiVUDVzu4W6vSYE3UOADChRwPUda/a2d2HUdSpZLy14STyNHqo5DI80bb2jjPOVmuxePslfH8wAc91CcLMwba7c9r1zAKs2BuHdUcSjWG1sbcTXunVCI+39QcMekRFnS5zG1q9AVGnkrFiXxxOX797gdajzbwwvlsD9GjiaTawtwt0w/LR7eHtokK7QDeLhPpstRZ/nriO2GtZmBbehD9DRERlqFSYNRgMmD9/Pi5evAiNRoM+ffrgww8/5AwGVmAcL+tpfrys3iDg3/ii+WWre/EXUL2zsj8cvIqrN/Ph5azEK70bV7uWh4FWb8C8f85j5Z2bVABA9p0zi9Z2KikLC7ZdQJu6rojsV/KsusEgYOOxJHy+9Twycouu0P8vKdPKVRa5ejMPS3Zdxm/Hr0N3Z/xAm7qumPRoY4Q39zGO09YaSt7kpViOWot1h69hzYF43Mgquh5AJZfi6dC6eL5bAzQsZ5iAVCpB/1bVnztZEAQcu3obPx++hs2nbpiG8t6Nqr19IqIHVaXC7KeffooPP/wQYWFhsLe3x1dffYX09HSsXr26puqjO4wzGZTyh/XMjSzkqHVwVtmhpX/VzqgGujtAIZOibaAbejSp2g0XbuYW4qudlwAAb/VrBiclP2otT1q2GlN+OoHDCUX/GXFW2Rk/Irekg1cyoDcI6NHEy+zruYU6fLHtAr4/mACDAJxMyiwRZk8k3saHm87gZFIWAMBBIUO+pvSgWFMSMvLw9c7L+CP2OvR3QmzXRnUwqXdjdGtcp0JnR29kFmDNgXisO3wNOXdmIvB0UuC5LkEY/Uh9eDhWfRaPyridp8FvJ65j3eFEXErLNS63k0qgMwjQG8oeFgEAWfla3MwrLDd4ExE9iCqVNH744QcsXboUL730EgBg+/btGDRoEFauXAmptPpjNKl05d0woXhKrs4NPKp821lvFxUOvdcXDgpZlT8q/XL7ReSodWjh54KhoXWrtI2HyaG4m5jy0wlk5BbCSWmHBc+0xp+xN/DP6RSL7eNWngazNp3BXydvwE4qQeysfiX+k7H1TApm/XkGKdl3Zyq59yZVaTlqfL7lAjYeSwIAOCntMLVvE/i6qvDqzycqVEdCRh42n0rGk+0C4O9WtU9z4tJz8c2dEFt8IVevpl54rW+TSl1k+MEfp7H/coYxCDf2dsLEHg2sNk5ZQNF/LtYdvoYtp1Og0RcFVpVcisda+2NEp0BsOJqEdUeulboNvUHA/ssZ2HD0GradTYVGZ8Afk7uhbaBbjddPRFSbVCrMJiYmmtyuNiwsDBKJBDdu3EDdugwuNeneGyaYU3zxV5dq3sK2OmejLqbm4Kd/EwEAMx5rUeVQ/TAQBAHf7Y3D51svQG8Q0MzHGctGt0dDLyf8GXvDYvvZcjoZH/xx2jgcQGcQUKDRG8Ps9cwCzPrzDLafSwUA1PNwwPhuQfjor7MAAI3OgO8PJmDxjkvGeVSfDq2Lt/s3g7ezClsqELpv5Wnw1Y5L+PHfq9DqBaTnFOLDx1tW6jgup+Xim52XsOnkDWOIfbSZF6aGNa1SeNtzMR1A0fjyF3s2RK+mXladOm7TyRsm/dzS3wXPdqqHJ9r6w+XOBWXF/3G4X1x6LjYeS8Jvx6+b/OcDAFKyCgCGWSJ6yFQqzOp0OqhUKpNlcrkcWq3WokWRqXyNDsl3xvIV3TDBlFZvwOHi8bIWuPirKgRBwOy/z8IgABEtfSwybvdBla3W4q0NJ7H1TFGAfLJdAD59spVFr36/mVuIWZvO4O//kgEAzXyccSH17u1WdXoD1h5MwMLoi8jX6CGXSfBSz0aY0qcxkm4XAAAKtHr0X7zX+B+p1nVd8eHjLdG+XsXOgKq1eqw5kICluy4bP8YvXl5Rl9Ny8NWOy/jrvxvGM8Vhzb3xWt8maF3XrcLbKdbQ0xEXUnPxWGs/TOzRsMoXOVaVwq7oEyxBKDq7/Xhbf4zoWA8h5Uyll1uow+b/bmDD0SSTeWrdHOQY0jYAey+mI66UuwISET3oKvXXUxAEjBs3Dkrl3VuSqtVqvPzyyybTc3FqLssqDhMejgq4mzlz+l9SFvI1erg7yBHs62zt8gAAuy+kY9+lDMhlEkwf0NwmNYjB+ZRsvPJ/xxGfkQeFTIqZg1tgVOd6Fp3WKupUMmb8cRo38zSQSSWY1LsopDb7YAuAoou1vth2EWfv3FK1Y5A7Pn0yBE19TL93NDoD4tLzUMdRgXf6B+Pp0LoVOntpMAj4I/Y6Fmy9YLygqoWfCwLc7RF9NrVCx3A5LQeLtl/C5lPJxhAb3sIHU/s2qVYA/XlCR0BiB1cHefmNa8CAVn64nJaLVv6uGNTaD44VGFO+8VgSluy6goI7/wmQSoqGVjzTIRB9m3tDaSfD08sOMswS0UOrUmF27NixJZaNHj3aYsWQecV/pEqbyeBQXPF42To2ucuWVm/AJ5uLPpYe1zUIQaXU+bD7/UQSpv92CmqtAQFu9lgyqr1Fxzdm5BZi1p9nsPlU0dnYYF9nzH+6TYmzfi98fxQA4Govx3sDg/FMaKDJ9427gxxSCSCRSDCuaxBe69sErvYVC3/7L2VgTtQ5Y1D2d1XhjX7N8GS7ACzbc6XcMJt4Mx+Ltl80GRMb0dIHr/VtUuULG+/loLCDXG6bIAsAXs5KfPxEq0qtk3DnjmMNvRzxTGggnmofAB8XVTlrERE9PCoVZtesWVNTdVAZ4sq5jW3xxV+2+mj/p38TcSU9Dx6OCkzp08QmNdRmOr0Bn0adw5oDCQCAHk08sfjZdha9Wv7v/25g5p9ncOues7Gv9mli/Fj7fk+2C8D7g5rD00lZ4rU6TkpETe0BR4Vdhe8Al5FbiLGrDxvHojor7TDp0cYY3y2oQhdUJWcV4Oudl/HLkWvGKbYiWvrg9bCmaO7nUqEaHjR9gn1w4PJNdGtcB0+HBqJ9PcvMYUtE9KDhvEkiUNa0XBqdAUevWm5+2crKKtDiy+0XAQDTwptW+AzewyKrQItXfz6BvXdC3mt9m2Bq3yYWuzguI7cQM/88jahTRRdiBfs6Y8Ezbcx+FB/R0geJtwrw/sDm6F7O1GvBvpULkAk385FwMx92UglGP1Ifr/VtUqGwnpFbiGW7r+B/h65Cc+c2yr2aeuHNfs3KHUf6oAtv4YPwFj62LoOIqNZjmBWBsm6YEHstE2qtAZ5OCjTxtv4ck0cSii5GaerjhBEdA62+/9osPiMPL3x/BHHpebCXy/Dl8DYWmVy/2I5zqXh743+4maeBnVSCSY82xpRHG5d6NvbbMR0stu9i907xNTDEF29HBFdomElWvhbf7buCNQcSjPPUdmrggTf7NUOnBh4Wr5OIiB5cDLO1nCAIiM8o/cxs8RCDRxpWbKL4mvL+oBal3mb3YXTwcgZe+fE4sgq08HdVYcXYDhYZ8wkABRo9Po06i/87VDQNWjMfZ3wxzPzZ2Jr2SEMPfD60NZr4OKFdBWc5iIm7iR6f7zTe4axNXVe80a9ZqbeLJSIiKgvDbC2Xkq1GvkYPmVSCembGL/4bb9vxskDRfJ+9mpq/q9TD6H+HruLDTWegNwhoV88N344JhbezZS7YOX09C1PXncCVO0NPXujeAG9FNLPKRP/m2MmkGFbJM/JX71zQ1MzHGW/0a4rwFj4MsUREVGUMs7Vc8XjZeh4OZj8+vp1fNMevreaXlUkleH8Qp+ICimZ1mP33WfwQcxVA0UVWc58KsUjQ1BuKbrKwMPoCtHoBPi5KLHimTam3pq2Nii8ma+DpiNfDmmBwa3+bzL5BREQPFobZWq54JoPS7vwFAD4uSjSw8nRY7eu5w81BjgndG6Cxt23mtq1NsvK1mPTTMRy4fBMSCfBWRDO80quRRc44Xs8sQOT6WPx758YY/Vv6Yu5TIWbnHK7NHm/jjzZ1XRHgZs8hKUREZDEMs7XclTJmMijWxQbjZZv4OOPEjHB+PIyioSBDlh5AfEYeHBQyLBreFv1a+lpk23/GXscHf5xGjloHB4UMHw5uiWc61BXt+16/DucgJiIiy2KYreXKu2ECAHRtVPY0SzVFrIHK0i6nFZ09D3Czx8qxHSw2L+q3e64Yby3brp4bFg1vyzBIRER0n1rxWd+SJUsQFBQElUqFzp074/DhwxVab926dZBIJBgyZEjNFmhDd2+YUMaZWRte/EVFOtR3x59Tull0gv+k2wWQSoCpfZtgw0tdGGSJiIjMsHmYXb9+PSIjIzFr1iwcP34cbdq0QUREBNLS0spcLyEhAW+++SZ69OhhpUqtT63V43pm0Zm50u7+FeBmX+G7NJFltQ10Q4CbPUZ2rocfJ3Y2ezetqnBWFX1gUs/DARte7opp4U05xpSIiKgUNh9msHDhQkycOBHjx48HACxfvhybN2/G6tWr8e6775pdR6/XY9SoUfjoo4+wb98+ZGZmWrFi64nPyIMgAC4qO9Qp5WIfnpW1nfp1HHHg3T4W3+4b/ZohtL47BrX2N7kpAREREZVk07+UGo0Gx44dw/Tp043LpFIpwsLCEBMTU+p6H3/8Mby9vfHCCy9g3759Ze6jsLAQhYWFxufZ2dkAAK1WC61WW80jKF/xPiqyL0Eoup2nXq+HVqvFpZQsAEVTGel0uvsbAwA61XezynE8zCrTh5bgYS/DU239AAjsWwuxdh9am3Dn94FOp39gj/FB78MHHftP/Kzdh5XZj03DbEZGBvR6PXx8TO8/7uPjg/Pnz5tdZ//+/Vi1ahViY2MrtI+5c+fio48+KrF827ZtcHCw3sfz0dHR5ba5eVMKQIoTJ05Ack3A1iQJABkUhbcRFRVl0raOTIpcmQSaxFhEJcfWSM1kqiJ9SLXbg9qHt27LAEhw/Phx6K8Kti6nRj2offiwYP+Jn7X6MD8/v8JtRfUZZk5ODsaMGYMVK1bA07NiV/BPnz4dkZGRxufZ2dkIDAxEv3794OJiuYt1SqPVahEdHY3w8HDI5fIy2/6ccgSXsm+jXbt2GBjii50bTwHXktGzTVMM7NXQpG2/CAO0egH2Ctvc+elhUpk+pNrpQe/DH64fRnxOJtq3b4+Ilj7lryBCD3ofPujYf+Jn7T4s/iS9ImwaZj09PSGTyZCammqyPDU1Fb6+JefpvHLlChISEjB48GDjMoOh6KN5Ozs7XLhwAY0aNTJZR6lUQqkseWGOXC636g9URfYnkRRd5COTySCXy5Fw57afTXxdSqwrlwP2NVMqlcLa3zNkeQ9qHxZPk2dnJ3sgj+9eD2ofPizYf+JnrT6szD5seom0QqFAaGgoduzYYVxmMBiwY8cOdOnSpUT74OBgnDp1CrGxscavxx9/HI8++ihiY2MRGFi5e8TXZoIgGG9lW9a0XEREREQPM5sPM4iMjMTYsWPRoUMHdOrUCYsWLUJeXp5xdoPnnnsOAQEBmDt3LlQqFVq1amWyvpubGwCUWC526bmFyCnUQSoB6tfh1FtERERE5tg8zA4fPhzp6emYOXMmUlJS0LZtW2zZssV4UVhiYiKk0odvjs0raUVnZeu6O0Bpx3GxRERERObYPMwCwJQpUzBlyhSzr+3evbvMddeuXWv5gmqBuIziO3/xrk9EREREpXn4TnmKhHG8rCfHyxIRERGVplacmaWS4tJ5ZpaIbEOt1WPPxXRsO5OKRt6OmNS7sa1LIiIqFcNsLRWXUXRmthFnMiAiC4rPyIPeIKCxt+nvFrVWj70X07H5VDK2n01FnkYPAJBJJQyzRFSrMczWQhqdAdduFc0x24hnZomompJu5+Pv/5KxKfYGziZnQyaV4PB7feGotMO+SxnY/N8NbD+XhtzCu7fN9nRSIiO3EHrDg31HMSISP4bZWujqrXwYBMBJaQcv55I3fCAiKk9ajhpR/yXjr/+ScezqbZPX9AYB0345ieNXb5sEWD9XFQaG+GFgiB8CPezR6dMd92+WiKjWYZithe4dL1t8Zx8iovJkFWix/kgiNp28gZgrN1F8UlUiATo38MDjbQIwN+occgp12HsxHUBRgB3Qyg+DWvuhXaAbpNKi3zkZuYW2OgwiokphmK2FrhhnMuAQAyKquHd+PWXyvG2gGwa38cdjrf3g46ICAJxIvI1/42+hb3NvPNbaD+0C3Y0BlohIjBhma6GEDN7Glogqzl5x98Yqwb7OGNzGH4Nb+6OembsHzn+mjTVLIyKqcQyztVCBtugqYk7LRUQV8f6g5th3MQO9mnmhqY+zrcshIrIqhtlajDdMIKKKCPZ1QbCvi63LICKyCd4BrJaSSIAGHDNLREREVCaG2VrK39XeZBwcEREREZXEMFtLcbwsERERUfkYZmsp3saWiIiIqHwMs7UUz8wSERERlY9htpbiTAZERERE5WOYraV4ZpaIiIiofAyztZC9XAbfO7eeJCIiIqLSMczWQg08HXmvdCIiIqIKYJithRp5c7wsERERUUUwzNZCDXnnLyJ6AKm1euy5mI6luy8j8Wa+rcshogeEna0LoLsclUXd0dzP2caVEBFZRlqOGrvOp2HHuTTsv5yBfI0eABCXnocFz7SxcXVE9CBgmK1F3h/UHL2beSG8ha+tSyEiqhJBEHDmRjZ2nEvDzvOpOJmUZfK6QiaFRm9ArlpnowqJ6EHDMFuLNPB0RAMOMSAikSnQ6LH/cgZ2nk/FzvNpSM0uNHm9TV1X9An2Qd/m3jhxLRMz/jhd7jZv5Wmw71I6km4XYGSnenBS8KJYIjKPYZaIiCrtemYBdp5Pw85zqTh45SYKdQbjaw4KGXo08UTfYB/0DvaCt/PdqQZjr2Wa3Z5Ob8DJpCzsuZCGPRfT8d/1LAhC0WtSiQQTutWrycMhIhFjmCUionLpDQJOJmVix7lU7DiXhvMpOSavB7jZI6y5N/o090HnBh5QyWXlbjM1W409F9Ox52I69l/KQFaB1uR1e7kMBVo98go5JIGISscwS0REZXrjl5PYfSENN/M0xmVSCRBa3904fKCJtxMkkooPBdh5Pg1bzqSYLHNR2aFHUy/0uvO1bPcVrD2YUOo2BEHAxdRc7LuUjv2XM5BXqMPy0aGo46Ss9DESkXgxzBIRUZl+PZ4EAHBW2aFXUy/0be6NXk294eGoqPS2nFVFf3Y0egMkEqB1XTdjeG1T1xV2srJnjEzLKcS/CanYfykD+y9nIC3HdHzu/ssZeKJtgPl1s9WIibuJg5dvwsXeDu8NbF6pAE5EtRPDLBERlVDHUYHwFj6Iz8hD76Ze6NvcBx2C3CEvJ2yWJ6KlL2Y81gKeTgr0aOJVoUB8ITUHc/65gC2xMiTH7DF5TSWXonODOriclovrmQXGcbYAkFWgxaG4m4i5chMHLmfgUlquybpDQ+si2NelWsdDRLbHMEtERCVIJBKseK6DxberksvwQvcGlVon+mxqcVWQSICQAFd0b+yJ7k08EVrfHUo7GUav/BfXMwvwX1IWzqfk4OCVDJy+ngXDPeFWIgFa+rvgUmouCnUG6PSC2f0RkbgwzBIRUa3Urp4bfogB/Fzt0a2RBxyyEzHp6b7wdi19CsPVB+JNnjf0ckS3Rp7o2qgOHmlYB+6OCjwyZwdSstU1XT4RWQnDLBER1UpPtA1A/1a+UMik0Ol0iIq6CncH88MSWgW4Yv/lDPi5qtC1kSe6Na6Dro084euqMtueiB4cDLNERFRrKe3Kn+ILAN7p3wwv92oIV3s5L+oiesgwzBIRkehJJBK4lXLWlogebNW7LJWIiIiIyIYYZomIiIhItBhmiYiIiEi0GGaJiIiISLQYZomIiIhItBhmiYiIzFBr9TiacAvXMwtsXQoRlYFTcxEREQHIKtDi+NXbOJxwC4fjb+FUUhY0egO8nJU4NL0vZFLOX0tUGzHMEhHRQyktR42//8vDkfhbOJxwG+dTsiEIJdul5xRCqzdAJq3YDRyIyLoYZomI6KH0/NqjJZYF1XFAxyAPdGzggRZ+Lnjs6/02qIyIKoNhloiIHir+biqkZKshkQDBvi7oFOSOTg3qoGOQO7xdVMZ2uYW6MrcjCAISb+Xj2NXbOHr1Ns7cyMYTbfzxfPcGNX0IRHQPhlkiInqofDumAy6l5aClvytc7eUVXq9Qp8fp61k4dvX2na9MZOQWmrS5lVfIMEtkZQyzRET0UPFyVsLLWVmpdZ797hDO3siGRm8wWa6QSdEqwAV+rvbYfCrZ7JhbIqpZDLNERERmyCQSSCWAQQBir2UCAOo4KhBa3x2h9d3RIcgdLf1doZLLEHstE5tPJdu2YKKHFMMsERGRGfYKGT5+ohXOp2SjXWBRgK1fxwESCafoIqpNGGaJiIhKMfqR+rYugYjKwTuAEREREZFoMcwSERERkWhxmAEREZEV3MwtxMmkTMRey0JCRh7GdKmPjkEeti6LSPQYZomIiCysQKPHmRtZiL2WidhrmTiZlIlrtwpM22j1DLNEFsAwS0REZCEZuYUYuHgfLqTmQG8oOelsIy9HOCrt8F9SFrT3zVlLRFXDMEtERFRNKnnRJShqrQFnk7MBFN2coW2gm/ErpK4rXFRybDyWhDc3nLRluUQPFIZZIiKiamrm44y3Ipohu0BbFF7rucHXRcU5aYmsgGGWiIiomiQSCSY/2tii21Rr9biQkoOM3EJ0a+wJlVxm0e0TPSgYZomIiGxMrdXjfEoOTl3PwumkLJy6noWLqTnQ3Rl3+/7A5pjYs6GNqySqnRhmiYiIbODqzXy8s/G/EsH1XlIJYBCA1Gy1DSokEgeGWSIiIisqHkUbn5GH+Iw843IPRwVaBbgiJMAFIQGuaBXgiv8duopv98TZplAikWCYJSIisqKeTb3wSEMPKO1kxtAaUtcV/q4lLxiTgBeQEZWHYZaIiMiKvJyVWPdiF1uXQfTAYJglIiISMbVWj4upOTh7IxsavQHPdqwHhZ3U1mURWQ3DLBERkUhk5Bbi7I1snEvOxtnkbJy9kY0r6bm499qxOo5KDGrtZ7siiayMYZaIiKiW++u/G9h08gbScgrNvu7hqIBWb0COWocctbba+xMEAclZalxIycHF1Bw09XHGo8He1d4uUU1gmCUiIqqlHBVFN0pIzS4KsRIJ0KCOI5r7u6CF350vfxd4Oysx8Yej2H4uzex2BEFA0u0CnE/JwYWUbJxLyYGLSo5PhrRCrlqHC6lFy8/fCa8XUnKQrdYZ11fIpDg5qx/sFbxxA9U+DLNERES11IjO9SCg6Mxrcz8XBPs6w1FZ9p/uHLUORxJu4XxyUTgtCrA5yC3UlWj78+HEUrdjJ5UgyNMRl9NyodEboNEZGGapVmKYJSIiqqU8nZR4rW+TSq3zadQ5s8vlMgkaeTmhuZ8Lfj9x3eS1ADd7NPN1LvryKfq3oZcjZBIJGr//T5XrJ7IGhlkiIqIHQFAdR+Njf1cVgv1c0MzXGcG+zgj2dUFDL0fIZUWzHLwe1gSH42+hoZcjmvg4w0UlN7tNnd5gldqJqoNhloiI6AHwzoBgDA2tC39Xe7g6mA+nxerXcUT9e8IvkZgxzBIRET0A5DIpmvu52LoMIqvjrMpEREREJFoMs0REREQkWgyzRERERCRatSLMLlmyBEFBQVCpVOjcuTMOHz5catsVK1agR48ecHd3h7u7O8LCwspsT0REREQPLpuH2fXr1yMyMhKzZs3C8ePH0aZNG0RERCAtzfxdTHbv3o0RI0Zg165diImJQWBgIPr164fr16+bbU9EREREDy6bh9mFCxdi4sSJGD9+PFq0aIHly5fDwcEBq1evNtv+xx9/xKRJk9C2bVsEBwdj5cqVMBgM2LFjh5UrJyIiIiJbs+nUXBqNBseOHcP06dONy6RSKcLCwhATE1OhbeTn50Or1cLDw8Ps64WFhSgsLDQ+z87OBgBotVpotdpqVF8xxfuwxr6oZrAPxY99KH7sQ9u496YJWp0WVX372X/iZ+0+rMx+bBpmMzIyoNfr4ePjY7Lcx8cH58+fr9A23nnnHfj7+yMsLMzs63PnzsVHH31UYvm2bdvg4OBQ+aKrKDo62mr7oprBPhQ/9qH4sQ+tSy8AxVEhOjoaDtVMDew/8bNWH+bn51e4rahvmvDZZ59h3bp12L17N1Qqldk206dPR2RkpPF5dna2cZyti0vNTy6t1WoRHR2N8PBwyOVl35GFaif2ofixD8WPfWgbOr0BkYe2AwDCw8Phal+19579J37W7sPiT9IrwqZh1tPTEzKZDKmpqSbLU1NT4evrW+a6CxYswGeffYbt27ejdevWpbZTKpVQKpUllsvlcqv+QFl7f2R57EPxYx+KH/vQuiTSu8MM5HbVf+/Zf+JnrT6szD5segGYQqFAaGioycVbxRdzdenSpdT1Pv/8c8yePRtbtmxBhw4drFEqEREREdVCNh9mEBkZibFjx6JDhw7o1KkTFi1ahLy8PIwfPx4A8NxzzyEgIABz584FAMybNw8zZ87ETz/9hKCgIKSkpAAAnJyc4OTkZLPjICIiIiLrs3mYHT58ONLT0zFz5kykpKSgbdu22LJli/GisMTEREild08gL1u2DBqNBk8//bTJdmbNmoUPP/zQmqUTERERkY3ZPMwCwJQpUzBlyhSzr+3evdvkeUJCQs0XRERERESiYPObJhARERERVRXDLBERERGJFsMsEREREYkWwywRERERiRbDLBERERGJFsMsEREREYkWwywRERERiRbDLBERERGJVq24aQIRERGJi8EgIC2nEIm38nHtVj4Sb+UjOasAYc190K+lr63Lo4cIwywRERGVa83BeGTkFuLarQJcu5WPpMwCaHSGEu0OXL7JMEtWxTBLRERE5Vq0/VKJZTKpBAFu9gj0sIejwg7bzqai0EzAJapJDLNERERklp1Mikm9G+FQ3E3UdXdAoIc96nk4INDdAYEeDvBzVcFOVnT5zfmUbGw7m2rjiulhxDBLREREpXq7f7CtSyAqE2czICIiIiLRYpglIiIiItFimCUiIiIi0WKYJSIiIiLRYpglIiIiItFimCUiIiIi0WKYJSIiIiLRYpglIiIiItFimCUiIiIi0WKYJSIiIiLR4u1siYiIyOI0OgPSctRIzS7E9Vu52J0swemtF5Geq0FKthpZBTq80rsRHm/jb+tSSeQYZomIiMhibuUVInR2NG7mae57RQYkJJgs+fnfRIZZqjaGWSIiIqq2Oo5KyKQS6A2CMcgqZFJ4uyjh7ayEIe8W2jYNgp+bA5JuF+B/h65CgGDjqulBwDBLRERE1eblrMRfU7ojNVsNHxcVfFyU8HBUQCKRQKvVIioqCgMHBkMul+Pv/27gf4euVmi7Gp0Bap0eLip5DR8BiRXDLBEREVlEC38XtPB3qXD7vEI9Yq7cRFqOGuk5hcavtJxC47Lb+VpIJMD8p9vg6dC6NVg9iRXDLBEREdnEqetZGLHiULntBAGIvXabYZbMYpglIiIiq2ob6AZvZyUKtHp4Oyvh5ayEt7Pqzr9Fz4uXbTh6DSv3x9u6ZKrFGGaJiIjIquq6O+Df9/pCIpGU29ZRyahCZeNNE4iIiMjqKhJkiSqCYZaIiIiIRIvn7omIiEi01Fo9buVpcCtPg5t5GshlEnRuUAcyKc/8PiwYZomIiKjW+79DiVBrDcbQeiuvELdyNcjT6Eu0XT66Pfq38rNBlWQLDLNERERUa6nkMuPjjceSzLaxk0rg4ahAXqEOeRo9UrLU1iqPagGGWSIiIqq1hoYGICO3EIIA1HFSoI6jAh6OCtRxUsDDseguYy4qO0gkEkz+6Tg2/5ds65LJyhhmiYiIqNbydlZhxmMtqr0dvUFAVoEWt/M1yMzX4HaeFnWcFGhXz90CVZItMcwSERHRA+XDv85iz8V03M7XFgXXfC2y1VoIQsm22yN7orG3s/WLJIthmCUiIqIHgpPibqzZdSHdbBtnpR3cHOVIzS6ERmdAanYhw6zIMcwSERHRA+HVvo0R4G4PpZ0U7g4KuDnI4e6ogJu9HG53nstlRVPs91+0F+dTcmxcMVkCwywRERE9EOq6O+C1vk1sXQZZGe8ARkRERESixTBLRERERKLFMEtEREREosUwS0RERESixTBLRERERKLF2QyIiIjooXUk4RY0OgOy1VrkFuqQq9YhR61DbqGuaNmdx8lZaozqXM/YJrdQh5x7HhetV7SNJ9oGYPaQVgCK7jyWp9Ehr7hN4d3HhToDujX2hJez0mxtgiAgX6NHXqEOAgAfF5UV3xnxYJglIiKih9ai7Zcq3PaTzecq1O5/h65i65kU5BbqkK/Rl9u+Xwsf5GnuBuO8wqIAm6vRmdy1bMZjLTC8YyDy7gTivEL9nfa6O4G5aD2VXIqhoXXhoHg4Yt7DcZRERERE93iqfQBW7ouHg0IGJ5UdnJXyO//awVllV/RYJYeT0g6b/0uGUi6Fi+puGyelaRtnlR0EARi96l8AQFpOocn+7KQSOKvs4Hhn3Xtv2LDtbGqFap7991nM/vtshdrayaQY0aleBd8NcWOYJSIioofOiz0b4cWejSrUdvQj9Su83YPv9sG1W/lwLA7FyqIAq7STQiKRGNsZDAI2Hk/C7TyNMeAW/1v0uChkOyntcDj+Fp5fewSGO2dppRLAUVHU3lEpg5PSDg53np9PyUbS7QJkFWgr9X6IGcMsERERkYX4u9nD382+3HZSqQTDOgRWaJu9m3nj+IxwaPUCnJR2UMlNg/G93vjlJJJuJ1WqZrFjmCUiIiKq5dwcFLYuodbi1FxEREREJFo8M0tERET0ABIEARq9AWqNAQVaPQq0eqiL/9UU/evppESbQDdbl1otDLNERERED5jP/jmPz7ecN140VpafJz6CRl6Od8KuwRh6734ZkKvW4HiyBD0LdXCXy2v+ACqBYZaIiIjoAdHM18n4+N4gK5NK4CCXQaWQwV5e9HUhtWh6sBErDlVw6zJMztPA3an8C9ysiWGWiIiI6AHxYs9G6N/SDwCgUkihuhNc5bKSl0nN/eccvt0TB6Bouq/itiq5DCq51OS5QibB7YxUKOxq3+VWDLNERERED5B6dRwq1G76gOZ4rU8TyGVSyGWSUqf7AgCtVouoqCj41sJb6jLMEhERET2kHJXij4K171wxEREREVEFMcwSERERkWgxzBIRERGRaDHMEhEREZFoMcwSERERkWgxzBIRERGRaDHMEhEREZFoMcwSERERkWgxzBIRERGRaDHMEhEREZFoMcwSERERkWiJ/4a8IuBccA2S838DMpn5Bg17AyqXosfpF4q+ShPUHXDwKHqccRlIO1t62/pdAUfPose34oGUU6W3DewEOPsWPc5MBG7Elt42IBRwDSh6nHUduH6s9Lb+bQG3ekWPc1KBa/+W3ta3FeDRsOhxXgZw9WDpbb1bAJ6Nix4X3Abi95Xe1qtZ0RcAqLOBuN2lt63TCPBpWfRYkwdc3gGJXg+/zGOQnDeY9qF7fcCvTdFjXSFwcWvp23WtCwS0L3qs1wEXokpv6+IP1O1w9/nZTaW3dfIG6j1y9/n5zYBBb76tQx0gqNvd5xe3FtVtjsoVaNjr7vNL2wFtvvm2SiegUZ+7z6/sBApzzbeVOwBNwu4+j98LFGSabytTAM36332ecADIv2m+rVQGBA+6+zzxXyA31fi0RB+2ePxu26SjQPYN89sFgGYDAdmdX5XXjwNZSaW3bdIPkKuKHif/B9xOKL1t476AwrHoceoZ4OaV0tvydwQkej0cCtPvLqhFvyNKxd8RRZROQL0ed5/Xwt8RJfB3RJF7f0fcToDUoCm9rQ0xzFpB4K0DsPu1jF9Mkw7d/SY8/Ruw57PS207Yefeb8MJmIHpm6W3H/g00uPML5PJ2IOrN0tuO/OXuH6r4fcCfk0pv+/QawPWposdJh4EN40pv+8RSoN2oosfJJ4FfxpTeduACoNOdMJt2ruy2YR8B3V8venwrruy2vd4FHp1e9Dj7Rtltu74K9Puk6HFeBvDLGNgB6AQA8fe17fAC8NjCoseFOWVvt81I4MllRY/1mrLbtngCGPbD3edltW0cDozeePf5rxNK/4NSvzswfvPd539MAvIzzLf1bwe8uPvu87+nAVmJ5tt6NgOmHL77fMt0IP28+bau9YBp9wSm6FnAjePm2zrUAd6Ou/t81xzg6n7zbe3sgQ9S7j7ftwC4tO3uy7ivDz/Mutv24NfA2T/MbxcA3rtx9w/V4RXAyZ9Kb/vm5bt/qI5/DxxZWXrbqf/d/UN1ch1w8KvS2/J3BOwA+NR97u6CWvQ7olT8HVHEsxnw0oG7z2vh74gS+DuiyD2/I6RXdkKuV5Xe1oYYZq0gX+EJQ93OkEok5hvI7e8+dq0LBD5ivh1w9xsbAJz9y25b/I0NAE4+5bR1u/vY0avstg517j629yi7raPXPW3dym7r5HNPPS5lt3Xxv/tY4VR2W9e6dx/LVWW3dat/97GdEgh8BAZBwO3bt+Du7mHahx4N7j6Wysrebp2Gdx9LpGW39Wxq+rystsVnk4rV7Vj6mRSfFqbPA0IBdZb5tp5NTJ/7tzV9z+9VfFatmG9r0++nezl539e2VdHZFXNUrqbPfVoABp35tnb3bcOrWdEZtjtK7UOg6FjLeo8l94zGqtOw7Laye36lugeV3dZOefexW72y2/J3BAyCALX8nu+JWvQ7olT8HVFEBL8jysTfEQAAwckbgiSv9LY2JBEEQbB1EUuWLMH8+fORkpKCNm3a4Ouvv0anTp1Kbb9hwwbMmDEDCQkJaNKkCebNm4eBAwdWaF/Z2dlwdXVFVlYWXFxcyl+hmrRaLaKiojBw4EDI5fIa3x9ZHvtQ/NiH4sc+FDf2n/hZuw8rk9dsfgHY+vXrERkZiVmzZuH48eNo06YNIiIikJaWZrb9wYMHMWLECLzwwgs4ceIEhgwZgiFDhuD06dNWrpyIiIiIbM3mYXbhwoWYOHEixo8fjxYtWmD58uVwcHDA6tWrzbZfvHgx+vfvj7feegvNmzfH7Nmz0b59e3zzzTdWrpyIiIiIbM2mY2Y1Gg2OHTuG6dOnG5dJpVKEhYUhJibG7DoxMTGIjIw0WRYREYE//vjDbPvCwkIUFt4dH5SdXTRGRqvVQqvVVvMIyle8D2vsi2oG+1D82Ifixz4UN/af+Fm7DyuzH5uG2YyMDOj1evj4+Jgs9/Hxwfnz5q90TElJMds+JSXFbPu5c+fio48+KrF827ZtcHBwqGLllRcdHW21fVHNYB+KH/tQ/NiH4sb+Ez9r9WF+fimzbpjxwM9mMH36dJMzudnZ2QgMDES/fv2sdgFYdHQ0wsPDOehdpNiH4sc+FD/2obix/8TP2n1Y/El6Rdg0zHp6ekImkyE11XTi4tTUVPj6+ppdx9fXt1LtlUollEplieVyudyqP1DW3h9ZHvtQ/NiH4sc+FDf2n/hZqw8rsw+bXgCmUCgQGhqKHTvu3kHFYDBgx44d6NKli9l1unTpYtIeKDrlXVp7IiIiInpw2XyYQWRkJMaOHYsOHTqgU6dOWLRoEfLy8jB+/HgAwHPPPYeAgADMnTsXADB16lT06tULX3zxBQYNGoR169bh6NGj+O6772x5GERERERkAzYPs8OHD0d6ejpmzpyJlJQUtG3bFlu2bDFe5JWYmAip9O4J5K5du+Knn37CBx98gPfeew9NmjTBH3/8gVatWtnqEIiIiIjIRmweZgFgypQpmDJlitnXdu/eXWLZM888g2eeeaaGqyIiIiKi2s7mN00gIiIiIqoqhlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEq1acdMEaxIEAQCQnZ1tlf1ptVrk5+cjOzsbcrncKvsky2Ifih/7UPzYh+LG/hM/a/dhcU4rzm1leejCbE5ODgAgMDDQxpUQERERUVlycnLg6upaZhuJUJHI+wAxGAy4ceMGnJ2dIZFIanx/2dnZCAwMxLVr1+Di4lLj+yPLYx+KH/tQ/NiH4sb+Ez9r96EgCMjJyYG/vz+k0rJHxT50Z2alUinq1q1r9f26uLjwB1jk2Ifixz4UP/ahuLH/xM+afVjeGdlivACMiIiIiESLYZaIiIiIRIthtoYplUrMmjULSqXS1qVQFbEPxY99KH7sQ3Fj/4lfbe7Dh+4CMCIiIiJ6cPDMLBERERGJFsMsEREREYkWwywRERERiRbDLBERERGJFsOsBSxZsgRBQUFQqVTo3LkzDh8+XGb7DRs2IDg4GCqVCiEhIYiKirJSpVSayvThihUr0KNHD7i7u8Pd3R1hYWHl9jnVvMr+HBZbt24dJBIJhgwZUrMFUrkq24eZmZmYPHky/Pz8oFQq0bRpU/4+taHK9t+iRYvQrFkz2NvbIzAwENOmTYNarbZStXS/vXv3YvDgwfD394dEIsEff/xR7jq7d+9G+/btoVQq0bhxY6xdu7bG6zRLoGpZt26doFAohNWrVwtnzpwRJk6cKLi5uQmpqalm2x84cECQyWTC559/Lpw9e1b44IMPBLlcLpw6dcrKlVOxyvbhyJEjhSVLlggnTpwQzp07J4wbN05wdXUVkpKSrFw5FatsHxaLj48XAgIChB49eghPPPGEdYolsyrbh4WFhUKHDh2EgQMHCvv37xfi4+OF3bt3C7GxsVaunASh8v33448/CkqlUvjxxx+F+Ph4YevWrYKfn58wbdo0K1dOxaKiooT3339f+O233wQAwu+//15m+7i4OMHBwUGIjIwUzp49K3z99deCTCYTtmzZYp2C78EwW02dOnUSJk+ebHyu1+sFf39/Ye7cuWbbDxs2TBg0aJDJss6dOwsvvfRSjdZJpatsH95Pp9MJzs7Owvfff19TJVI5qtKHOp1O6Nq1q7By5Uph7NixDLM2Vtk+XLZsmdCwYUNBo9FYq0QqQ2X7b/LkyUKfPn1MlkVGRgrdunWr0TqpYioSZt9++22hZcuWJsuGDx8uRERE1GBl5nGYQTVoNBocO3YMYWFhxmVSqRRhYWGIiYkxu05MTIxJewCIiIgotT3VrKr04f3y8/Oh1Wrh4eFRU2VSGarahx9//DG8vb3xwgsvWKNMKkNV+nDTpk3o0qULJk+eDB8fH7Rq1Qpz5syBXq+3Vtl0R1X6r2vXrjh27JhxKEJcXByioqIwcOBAq9RM1Veb8oyd1ff4AMnIyIBer4ePj4/Jch8fH5w/f97sOikpKWbbp6Sk1FidVLqq9OH93nnnHfj7+5f4oSbrqEof7t+/H6tWrUJsbKwVKqTyVKUP4+LisHPnTowaNQpRUVG4fPkyJk2aBK1Wi1mzZlmjbLqjKv03cuRIZGRkoHv37hAEATqdDi+//DLee+89a5RMFlBansnOzkZBQQHs7e2tVgvPzBJVw2effYZ169bh999/h0qlsnU5VAE5OTkYM2YMVqxYAU9PT1uXQ1VkMBjg7e2N7777DqGhoRg+fDjef/99LF++3NalUQXs3r0bc+bMwdKlS3H8+HH89ttv2Lx5M2bPnm3r0kiEeGa2Gjw9PSGTyZCammqyPDU1Fb6+vmbX8fX1rVR7qllV6cNiCxYswGeffYbt27ejdevWNVkmlaGyfXjlyhUkJCRg8ODBxmUGgwEAYGdnhwsXLqBRo0Y1WzSZqMrPoZ+fH+RyOWQymXFZ8+bNkZKSAo1GA4VCUaM1011V6b8ZM2ZgzJgxmDBhAgAgJCQEeXl5ePHFF/H+++9DKuW5ttqutDzj4uJi1bOyAM/MVotCoUBoaCh27NhhXGYwGLBjxw506dLF7DpdunQxaQ8A0dHRpbanmlWVPgSAzz//HLNnz8aWLVvQoUMHa5RKpahsHwYHB+PUqVOIjY01fj3++ON49NFHERsbi8DAQGuWT6jaz2G3bt1w+fJl439EAODixYvw8/NjkLWyqvRffn5+icBa/B8TQRBqrliymFqVZ6x+ydkDZt26dYJSqRTWrl0rnD17VnjxxRcFNzc3ISUlRRAEQRgzZozw7rvvGtsfOHBAsLOzExYsWCCcO3dOmDVrFqfmsrHK9uFnn30mKBQKYePGjUJycrLxKycnx1aH8NCrbB/ej7MZ2F5l+zAxMVFwdnYWpkyZIly4cEH4+++/BW9vb+GTTz6x1SE81Crbf7NmzRKcnZ2Fn3/+WYiLixO2bdsmNGrUSBg2bJitDuGhl5OTI5w4cUI4ceKEAEBYuHChcOLECeHq1auCIAjCu+++K4wZM8bYvnhqrrfeeks4d+6csGTJEk7NJWZff/21UK9ePUGhUAidOnUSDh06ZHytV69ewtixY03a//LLL0LTpk0FhUIhtGzZUti8ebOVK6b7VaYP69evLwAo8TVr1izrF05Glf05vBfDbO1Q2T48ePCg0LlzZ0GpVAoNGzYUPv30U0Gn01m5aipWmf7TarXChx9+KDRq1EhQqVRCYGCgMGnSJOH27dvWL5wEQRCEXbt2mf3bVtxvY8eOFXr16lVinbZt2woKhUJo2LChsGbNGqvXLQiCIBEEns8nIiIiInHimFkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIiIiIhIthlkiIiIiEi2GWSIiIiISLYZZIqKHmEQiwR9//AEASEhIgEQiQWxsrE1rIiKqDIZZIiIbGTduHCQSCSQSCeRyORo0aIC3334barXa1qUREYmGna0LICJ6mPXv3x9r1qyBVqvFsWPHMHbsWEgkEsybN8/WpRERiQLPzBIR2ZBSqYSvry8CAwMxZMgQhIWFITo6GgBgMBgwd+5cNGjQAPb29mjTpg02btxosv6ZM2fw2GOPwcXFBc7OzujRoweuXLkCADhy5AjCw8Ph6ekJV1dX9OrVC8ePH7f6MRIR1SSGWSKiWuL06dM4ePAgFAoFAGDu3Ln44YcfsHz5cpw5cwbTpk3D6NGjsWfPHgDA9evX0bNnTyiVSuzcuRPHjh3D888/D51OBwDIycnB2LFjsX//fhw6dAhNmjTBwIEDkZOTY7NjJCKyNA4zICKyob///htOTk7Q6XQoLCyEVCrFN998g8LCQsyZMwfbt29Hly5dAAANGzbE/v378e2336JXr15YsmQJXF1dsW7dOsjlcgBA06ZNjdvu06ePyb6+++47uLm5Yc+ePXjsscesd5BERDWIYZaIyIYeffRRLFu2DHl5efjyyy9hZ2eHoUOH4syZM8jPz0d4eLhJe41Gg3bt2gEAYmNj0aNHD2OQvV9qaio++OAD7N69G2lpadDr9cjPz0diYmKNHxcRkbUwzBIR2ZCjoyMaN24MAFi9ejXatGmDVatWoVWrVgCAzZs3IyAgwGQdpVIJALC3ty9z22PHjsXNmzexePFi1K9fH0qlEl26dIFGo6mBIyEisg2GWSKiWkIqleK9995DZGQkLl68CKVSicTERPTq1cts+9atW+P777+HVqs1e3b2wIEDWLp0KQYOHAgAuHbtGjIyMmr0GIiIrI0XgBER1SLPPPMMZDIZvv32W7z55puYNm0avv/+e1y5cgXHjx/H119/je+//x4AMGXKFGRnZ+PZZ5/F0aNHcenSJfzvf//DhQsXAABNmjTB//73P5w7dw7//vsvRo0aVe7ZXCIiseGZWSKiWsTOzg5TpkzB559/jvj4eHh5eWHu3LmIi4uDm5sb2rdvj/feew8AUKdOHezcuRNvvfUWevXqBZlMhrZt26Jbt24AgFWrVuHFF19E+/btERgYiDlz5uDNN9+05eEREVmcRBAEwdZFEBERERFVBYcZEBEREZFoMcwSERERkWgxzBIRERGRaDHMEhEREZFoMcwSERERkWgxzBIRERGRaDHMEhEREZFoMcwSERERkWgxzBIRERGRaDHMEhEREZFoMcwSERERkWj9P29Z3PAugDI/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate precision and recall values at different probability thresholds\n",
    "precision, recall, thresholds = precision_recall_curve(best_labels, best_scores)\n",
    "\n",
    "# Calculate the area under the precision-recall curve\n",
    "auc_pr = auc(recall, precision)\n",
    "\n",
    "# Add precision-recall curve for a random model\n",
    "random_predictions = np.random.rand(len(best_labels))\n",
    "random_precision, random_recall, _ = precision_recall_curve(best_labels, random_predictions)\n",
    "random_auc_pr = auc(random_recall, random_precision)\n",
    "\n",
    "positive_ratio = np.sum(best_labels) / len(best_labels)\n",
    "random_skill_recall = np.linspace(0, 1, len(best_labels))\n",
    "random_skill_precision = np.full_like(random_skill_recall, positive_ratio)\n",
    "random_skill_auc_pr = auc(random_skill_recall, random_skill_precision)\n",
    "\n",
    "\n",
    "# Plot the precision-recall curve\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(recall, precision, label=f'AUC-PR = {auc_pr:.2f}')\n",
    "plt.plot(random_skill_recall, random_skill_precision, linestyle='--', label=f'no-skill (AUC-PR = {random_skill_auc_pr:.2f}')\n",
    "\n",
    "plt.xlabel('Recall')\n",
    "plt.ylabel('Precision')\n",
    "plt.title('Precision-Recall Curve')\n",
    "plt.legend(loc='best')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}