3_clinical_finetuning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fine-tuning on clinical samples\n",
    "In this notebook, we'll demonstrate fine-tuning a CNN on clinical samples to illustrate the procedure described in the `Extension to clinical isolates` section and also in the `CNN architecture & training details` section in the Methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load clinical data and pre-trained CNN\n",
    "The first step is to load the clinical data that we're going to use for fine-tuning. This first set of data was collected in 2018."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import time\n",
    "t00 = time()\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10000, 1000) (10000,)\n"
     ]
    }
   ],
   "source": [
    "# Load clinical data\n",
    "X_fn = './data/X_2018clinical.npy'\n",
    "y_fn = './data/y_2018clinical.npy'\n",
    "X = np.load(X_fn)\n",
    "y = np.load(y_fn)\n",
    "print(X.shape, y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we load the ResNet model that has been pre-trained for the 8 class antibiotic treatment groupings task on the reference dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from resnet import ResNet\n",
    "import os\n",
    "import torch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CNN parameters\n",
    "layers = 6\n",
    "hidden_size = 100\n",
    "block_size = 2\n",
    "hidden_sizes = [hidden_size] * layers\n",
    "num_blocks = [block_size] * layers\n",
    "input_dim = 1000\n",
    "in_channels = 64\n",
    "n_classes = 8 # instead of 30, we use the 8 empiric groupings\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = '{}'.format(0)\n",
    "cuda = torch.cuda.is_available()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load trained weights for demo\n",
    "cnn = ResNet(hidden_sizes, num_blocks, input_dim=input_dim,\n",
    "                in_channels=in_channels, n_classes=n_classes)\n",
    "if cuda: cnn.cuda()\n",
    "cnn.load_state_dict(torch.load(\n",
    "    './clinical_pretrained_model.ckpt', map_location=lambda storage, loc: storage))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fine-tuning\n",
    "Now we can fine-tune the pre-trained CNN on our first set of clinical samples. In the experiments reported in the paper, we fine-tune within a cross-validation loop, but here we show just one fold for clarity. We also only train for one epoch here in the interest of time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train/val/test split\n",
    "For the clinical experiments, we need to be careful about splitting our spectra into training, validation, and test sets for each fold. In particular, we want to make sure that spectra from the same patient do not appear in both training and validation/test.\n",
    "\n",
    "Our 2018 clinical dataset contains 400 spectra from each of 25 patients (10k spectra total), with 5 patients coming from each of the 5 empiric treatment groups. For each treatment group, we randomly assign 3 patients to the training set, 1 patient to validation, and 1 patient to the test set. We then randomly sample 10 spectra from each patient for fine-tuning, validation, and testing -- this is to simulate a culture-free clinical setting where not many bacterial spectra may be available.\n",
    "\n",
    "The following code demonstrates how this procedure can be implemented for a single fold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each treatment group of 5 patients, patients are assigned to [train, train, train, val, test]\n",
    "patient_idxs = []\n",
    "for i in range(5):\n",
    "    x = list(range(5))\n",
    "    np.random.shuffle(x)\n",
    "    patient_idxs.append(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group 0 patients\n",
      " Tr: [3, 1, 2]\n",
      " Val: 4\n",
      " Te : 0\n",
      "Group 1 patients\n",
      " Tr: [3, 2, 1]\n",
      " Val: 4\n",
      " Te : 0\n",
      "Group 2 patients\n",
      " Tr: [1, 3, 2]\n",
      " Val: 4\n",
      " Te : 0\n",
      "Group 3 patients\n",
      " Tr: [2, 1, 3]\n",
      " Val: 0\n",
      " Te : 4\n",
      "Group 4 patients\n",
      " Tr: [2, 0, 1]\n",
      " Val: 4\n",
      " Te : 3\n"
     ]
    }
   ],
   "source": [
    "# Sample train/val/te spectra\n",
    "idx_tr, idx_val, idx_te = [], [], []\n",
    "for group_idx, patient_list in enumerate(patient_idxs):\n",
    "    print('Group {} patients'.format(group_idx))\n",
    "    print(' Tr: {}'.format(patient_list[:3]))\n",
    "    print(' Val: {}'.format(patient_list[3]))\n",
    "    print(' Te : {}'.format(patient_list[4]))\n",
    "    for j, patient in enumerate(patient_list):\n",
    "        start_idx = 5 * 400 * group_idx + 400 * patient\n",
    "        idx_range = list(range(start_idx, start_idx+400))\n",
    "        # Randomly sample 10 spectra for each patient\n",
    "        np.random.shuffle(idx_range)\n",
    "        idx_sample = idx_range[:10]\n",
    "        if j < 3:\n",
    "            idx_tr.extend(idx_sample)\n",
    "        elif j ==3:\n",
    "            idx_val.extend(idx_sample)\n",
    "        else:\n",
    "            idx_te.extend(idx_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fine-tune the pre-trained CNN\n",
    "Now that we've split our clinical data, we can fine-tune the pre-trained model on the clinical data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datasets import spectral_dataloader\n",
    "from training import run_epoch\n",
    "from torch import optim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting fine-tuning!\n",
      " Epoch 1: 0.00s\n",
      "  Train acc: 92.00\n",
      "  Val acc: 74.00\n",
      "  Test acc: 96.00\n",
      "Finished: 15.34s\n"
     ]
    }
   ],
   "source": [
    "# Fine-tune CNN\n",
    "epochs = 1 # Change this number to ~30 for full training\n",
    "batch_size = 10\n",
    "t0 = time()\n",
    "# Set up Adam optimizer\n",
    "optimizer = optim.Adam(cnn.parameters(), lr=1e-3, betas=(0.5, 0.999))\n",
    "# Set up dataloaders\n",
    "dl_tr = spectral_dataloader(X, y, idxs=idx_tr,\n",
    "    batch_size=batch_size, shuffle=True)\n",
    "dl_val = spectral_dataloader(X, y, idxs=idx_val,\n",
    "    batch_size=batch_size, shuffle=False)\n",
    "dl_te = spectral_dataloader(X, y, idxs=idx_te,\n",
    "    batch_size=batch_size, shuffle=False)\n",
    "# Fine-tune CNN for first fold\n",
    "best_val = 0\n",
    "no_improvement = 0\n",
    "max_no_improvement = 5\n",
    "print('Starting fine-tuning!')\n",
    "for epoch in range(epochs):\n",
    "    print(' Epoch {}: {:0.2f}s'.format(epoch+1, time()-t0))\n",
    "    # Train\n",
    "    acc_tr, loss_tr = run_epoch(epoch, cnn, dl_tr, cuda,\n",
    "        training=True, optimizer=optimizer)\n",
    "    print('  Train acc: {:0.2f}'.format(acc_tr))\n",
    "    # Val\n",
    "    acc_val, loss_val = run_epoch(epoch, cnn, dl_val, cuda,\n",
    "        training=False, optimizer=optimizer)\n",
    "    print('  Val acc: {:0.2f}'.format(acc_val))\n",
    "    # Test\n",
    "    acc_te, loss_te = run_epoch(epoch, cnn, dl_te, cuda,\n",
    "        training=False, optimizer=optimizer)\n",
    "    print('  Test acc: {:0.2f}'.format(acc_te))\n",
    "    # Check performance for early stopping\n",
    "    if acc_val > best_val or epoch == 0:\n",
    "        best_val = acc_val\n",
    "        no_improvement = 0\n",
    "    else:\n",
    "        no_improvement += 1\n",
    "    if no_improvement >= max_no_improvement:\n",
    "        print('Finished after {} epochs!'.format(epoch+1))\n",
    "        break\n",
    "print('Finished: {:0.2f}s'.format(time()-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions on treatment groups for test patients\n",
    "Now that we have a fine-tuned model, we demonstrate how it can be used for the analysis presented in Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting test indices\n",
    "idx_te = []\n",
    "for group_idx, patient_list in enumerate(patient_idxs):\n",
    "    start_idx = 5 * 400 * group_idx + 400 * patient_list[-1]\n",
    "    idx_te += list(range(start_idx, start_idx+400))\n",
    "dl_te = spectral_dataloader(X, y, idxs=idx_te,\n",
    "    batch_size=batch_size, shuffle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from training import get_predictions\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished: 46.68s\n"
     ]
    }
   ],
   "source": [
    "t0 = time()\n",
    "# Make predictions on the 5 test patients\n",
    "y_hat = get_predictions(cnn, dl_te, cuda)\n",
    "# Subset y according to which patients are chosen for the test set to generate y_te\n",
    "y_te = y[idx_te]\n",
    "print('Finished: {:0.2f}s'.format(time()-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate predictions of treatment groups\n",
    "We use 10 spectra per patient to make our treatment group predictions. Here we only make predictions for one set of 5 test patients; in the experiments reported in the paper, we fine-tune within a cross-validation loop to cycle all patients through the test set. To compute identification accuracies across all patients, repeat the fine-tuning process in a cross-validation loop so that each patient appears in the test set once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ran 10000 trials: 3.98s\n"
     ]
    }
   ],
   "source": [
    "t0 = time()\n",
    "\n",
    "# Experiment parameters\n",
    "trials = 10000\n",
    "n_spectra = 10\n",
    "\n",
    "preds2018 = np.zeros((5,trials))\n",
    "y_2018 = np.zeros((5, trials))\n",
    "\n",
    "for i in range(5):\n",
    "    y_2018[i] = y_te[400*i]\n",
    "    for k  in range(trials):\n",
    "        y_patient = y_hat[400*i:400*i+400]\n",
    "        idxs = np.random.permutation(np.arange(len(y_patient)))[:n_spectra]\n",
    "        preds2018[i,k] = stats.mode(y_patient[idxs])[0][0]\n",
    "\n",
    "print('Ran {} trials: {:0.2f}s'.format(trials, time()-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "from sklearn.metrics import confusion_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "from config import antibiotics, ab_order\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/inkjet/lib/python3.7/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in true_divide\n",
      "  \"\"\"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAFoCAYAAAAGr79JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXeclNX1h58vqDSxK7aIREWs0VhiDWDUqPEXbGiswRaxxxo0Kvbegg1Q7EbsBQuiIhpRUSQKKmIDlCIRkbqACOf3x73DDrOzO7s778zuzJ6Hz/uZmffeuee+M8t75p57iswMx3Ecx2msNGvoCTiO4zhOTbiichzHcRo1rqgcx3GcRo0rKsdxHKdR44rKcRzHadS4onIcx3EaNa6oHMdxnEaNKyrHcRynUeOKynEcx2nULNfQE3Acp/yQ1A5oD7TMbDOzt4o/I6eUcUXlOE5iSFoHeAjomq0ZMKB5USfllDyuqBzHSZLbgW2BfwCjgYUNOx2nHJAnpXUcJykkTQfON7N7G3ouTvngzhSO4ySJAd829CSc8sJNf47jJMlTwH7Aaw09kUIhaTfgSLI7i5iZ/aH4sypvXFE5jpMkjwD3SgIYBPyY2cHMRhd7Ukkh6VhgADAD+IKqe3Aq+qSaAL5H5eRE4a7zW6p3N/530SflNEokLUl7mXlzEWHFUbJef5LGEpxEjjEzdxQpEr6icmpE0sbAc0Ansv9aNMAVlZPi2IaeQIFpD5zuSqq4uKJycnEH0JZgk3d344SIq9QdqX6V+mDRJ5UAZvZAQ8+hwHwOrNXQk2hquOnPqRFJM4GeZjawoedSLkjqCDwPbEI1q9RCmMcktQQ2ILtiLNl9o2IiaT/gOmAvM/u+oefTVPAVlZOLCuCnhp5EmXEH0Ao4jCKsUiWtDvQHulFVMeadLULSvcAVZjY+Pq8JM7Pj6yurEXACsBLwpaRRVHUWMTM7uPjTKm98ReXUiKTrgDXN7Lgiy90CuAToDKwO/M7MRkm6DHjLzF4v5nySRNIs4G9m9liR5D0N7AncA4wFfs7sk4/JTtJ44AAz+1jSBKo6UWSIsl/XV1ZDU+7X11jxFZWTixHADZIep3p345eSFChpe+BNggvwq8ARac0tgZ5AySoqYB4wq4jy9gDOMrMBhRjczDqkPd+wEDIaC+V+fY0VX1E5NZLhbpyOUSB3Y0mvAcsDfwR+IawAto8rqu7A9ek3x1JD0jXAOmbWo0jyvgOON7MhxZDnOEnjKyonF3s1gMydgCPMbIGkTCU4DWjXAHNKklHAtdEkV90q9fkE5d1H8NosuKKSdCiwoZldn6XtfOAbM3uy0PNIEkkbAFPNbFF8XiNm5imkEsYVlVMjDbQXtARYXE3b6sD8Is6lEKT2pjoAB2RpT7oUxmVAP0mvAi8STKrLCkzOHf48Kq8vk0WxvaQUFTAe2Bl4H5hAzXtU4GVMEscVldMYGQX8hXBTzeRAwr5ZKZOtVlMh2QrYG1gfyJaHzoCkFFVH4KNq2kYDvROSU0yOA75Oe+77JUXGFZVTBUlDCNH34+LzmjAz+2PCU7geGCRpOULuOAO2l3Q0cDjhpluymNmbRRbZl7BK7UE1Xn8J0hxYsZq2toS9x5Ii3SPSzO5vwKk0WVxROdloRWW8TWuK/AvSzF6SdCJwI3BonEtfYDZwopm9Ucz5lAFbA0eZ2dNFkDUGOAR4NkvbocCnRZiDU2a4onKqYGa7pz3frYHmcK+kgYS9gXbAdGC4mc1riPnki6ShwClm9nl8XhNJl4qYTNgfKgZ3AA9KqgD6AZMIJseehADnks4FmCOgeQkh7OAD4BnPB5gcrqicRouZVVDa8VLppGeEaEbNq9SkS0XcDJwl6RUzK6TZDzN7WNLmwPlAegaKJcA1pZrDMI2uwMrAKoTQiR8JDj7LATMJ391ZwDhJXcxsWkNNtJzwOCqnRiSdDLQ3s15Z2q4BJphZvwLILcukrQ2BpJsJTijNgaFU9fozMzsnYZkdCKENawA/AK+a2YQkZTQEknYgeC2eDzxpZotjCEV3Qg7A7gSl9TTwgpmd0GCTLSNcUTk1Imk00N/Mbs/SdjJwkpltk7DMBknaWq7UELSdwj/PWiLpLeAJM7stS9uZQHcz203SWcC5ZrZe0SdZhrjpz8nFr4HPqmn7PLYnTVGTthYbSacSVqnnZ2m7nrBKvTMpeWbWLKmx6oKkNQnf4zKUeEDsdoS4tGx8Alwdn39EWE06CeCKysnFYmC1atpWpzClt3ckJG0ttcDQ2tKT4MWYjW+Ak4HEFFUxiSbbS4FTgVWr6VbKq7fZhETJ2fZOu8Z2gBbAnGJNqtxpkF9aTkkxihB/k43jgP8WQGaxk7YWmw6EeKZsjKMwq1Qk7SvpWkn3pFIBSdpJUpIpqU4CzgZuIfyIuToe3wFfEcpklDL/BnpJukzSlpLWjI+XE/atHo79fkuwODgJ4HtUTo1I2p+wXzSYqu7GexPKOwxKWGZRk7YWm1iM8kQzeyJLW3dggJmtlKC8FQk5BTsTVsjNgB1ikt+BwBQzOzshWf8FBhJi4BZRmUy4JfAG8LyZXZOErIZA0grA/YTMKek3TwGPAj1iTsC9gTlm9m7xZ1l+uKJyciLpFIJHU+vUKWAucL6ZVWfCykded+Ba4GOKk7S1qMQ4qoVmtm+WtleAlmbWOUF5txFKpfQgmKzmUqlAjgX+bma/SUjWHODPwDCCotrNzN6Lbd2B68qhXlN0we9MMIv/SKiRVt1erpMnvkfl5MTM7pT0MLAble7Gw81sds3vrDfFTtpabG4EXojpqTJXqXsQbvRJcjBwkZkNypKN/jtCefqkmA80NzOTNA3YEHgvts2m9DPfAxCVkiumIuGKyqkRSa3NrCIqpUQLJNZAsZO2FpWYIqoncAOVSWJFuJH3NLNsyXjzYTWq3y8RYeM/KT4nKCcIyYPPkvQfwurqfODLBGU1GHFfr7oYv7eKP6PyxhWVk4vJkh4A7jKzccUQ2ABJW4uOmfWX9G9gFypXqe8UKEXUtwS36mw5ErcjWeUxENg0Pr8SeCvKh5DJ4ZAEZRUdSesAD1H5Yyrl9bq0kCilvdpvlPgelVMjkvoARwMrEbIa3A4MMrNcQaRJyF6JUERxdcI+wHsFNDcWDUkrmtncIsq7HDiNkHn+NcLqZjvCTfUV4KZshQ4Tkt0e6BZlvWpmJe0JJ+kpoAtwDdXE+DWFH1rFxhWVkxNJrQnKqifwG8K+Rl/gHjP7oUAyLwQuZNlM7hXA5YW6qRYLSbMJ9Z/uLMYGvKQWBLNtF4LCXwP4HliLsMra18x+KfQ8ygFJ0wlORDUlp3USxuOonJzEPap+ZrYt8HtgOKEA3neSHpK0U5LyYomPKwml0w8mZFA/iPDr/5roqVbK3ENY3YyR9Iakg7M4OSRGzOK9F+HHxiuEVdUwQibzfZJUUpJOr6FtBUl3JSWrgTAqTZlOkfAVlVNnJO0FXATsTmVczn+AY81sfALjjwE+MLPjsrTdC2yXlDt1QxHjio4krFK3A6YQPADvNrPvG3Ju+RDzCj4DHGdms9LObwo8DmxkZtUVVmz0SOoLVCQVd+bUDldUTq2QtDLhF3hPQrLYUcBthEzS+xPinian17LKQ9YCoJuZvZKl7Y/As2ZWJYdcAnIbxJMrrkhPIWTebkbIvH2Hmb1dKJmFQtIBwABCyYu/mNkHko4h7G1OAg41s08aco75IGl34F5CfF91MX6jiz2vcscVlVMjkn5LuIn+hVBG/CngtsyI+6hAnjOzKjf5esicDpyTXgI8re1Y4HozWzNfOWljZnpyLdNMEbKLS9qXsCe3K8HZYTngXeCvZvZ1nmMvD/yDsIJrT1V3dDOzxDyAowPFY8A2wNuEz/Vh4ORYY6xkychEn3nzLMrfSlPE3dOdXIwEpgE3AX3NbGo1/cYTKpsmwevAlZJGmtnS0uWStiRkrn41ITkpbge2JdzMi5atXdJqhHyJPQnBze8DRxF+DPwJuB54gBBonQ/XEIr5vUwwyxX0+sxsYiwB8w4hgPl94Pgycdgo9f3RksRXVE6NSDoaGGhmxSplniq69w6wJkFxfA+sDWxNUJo7J1kqotieXJJ2JKxSDyWY+h4nrFI/yOi3NyE3Xl6rVEnfAvea2aX5jFMHeccBfYCvCSuriwnf42HlUDzRKT7u9efUiJk9VEwlFWWOJ5iNbiKYwTYBfiasMLYpQD2jYntyvUfwwrsG2MDMjslUUpEJBPNfvqwKFCW2JwaH30PIMv47M7uaENS8GvBfSSUd8Os0DL6icnIiaQ1C0GY2RwMzs38Uf1bJUWxPLkmHE6rEFsUUJull4BUzu7UIsmYTqj4/mnG+LXA3oQJuye7hRK/TmjAzO74ok2lCuKJyakTSH4BngTbVdEl881jSdmb2YQ3th5nZY9W110NeWXtySepI+A57AS8VUkFK2sTMqk3JJOkkM+tXKPmFRtIEqjpRrA6sSPB0nFkO2eEbG66onBqR9DEhI/YphGJ/P2f2MbPFCctcCJxrZrdlnG9BcIk/Pknl2BCeXJLWAg6k+lXqOQnK+ong6deCEPeWmb7JzGz1pOQ1RSTtQajK3N3MxjT0fMoN9/pzcrExcJCZjSqizLuAf0nqSgginiWpEzFgFPhbwvKK6skVA6afobK+VyYGJKaogOeoqoALRjQVn0VwS0/laXwduNXMqqxWywEzGyrpduBfBE9HJ0F8ReXUiKSRwA1JmtpqKbcbwRw3Kz72Ijg8lHTAKCzNvDGH4JY+ttjOKoUkxlANB9YBxgBT4/OtCAG/u5rZpIabYeGIq6rnSznzRmPFFZVTI5L2IdRN2rfYN5hYRfVDYAVCjFZXM5tfzDkUAknzgAPNbEhDzyVpYumSzsB+ZvZx2vmtCYlxh5nZUQ01v0Ii6TJCgPaGDT2XcsNNf06NmNlgSbsCX0r6FJhRtYv9MWm5krYAniDUMPoPsCchIe15SaxAovfWFWY2vgE8uT4luGsXjfh5XkJQIqsTXMdHxZvrW2b2ekKi9iLEpH2cftLMRku6BLguITkNQryGTFoQYvz2JfyocxLGFZVTI5JOA/4JLADakmw12OpkHkvIFvEVsL2ZjZN0CiGualdJhyaQ/LYrYT8Bwp5CTaaFpM0O/wRukfSumU1MeOwqSNqeEEc1g5DV44i05pYEE2RSiqotMLmatskE77hS5tIs5xYSYt56E2L9nIRx059TI5K+I2SJOM4KU302m8wlhJibM81sQdr5bQkOFWuY2arFmEuhkHQlcC6hdHu2Veofqr6r3rJeI+Rp/CNhhfoz4QfAKEndCbkTOyQkawyhwOWJWdoGADuY2dZJyHKaDr6icnKxGtC/WEoqcoSZDcw8aWb/jUly+xdxLokj6UxCAtp5BAVSaKW7E+EzXZCl7tU0oF2Csu4A7ozVmR+hMv3VUYTaYj0TlFUUJH1D2FP8OJr+7jGzKQ09r6aEKyonFx8QXNSTMg3lJJuSSmubQyg6mBeSNqjjnJJMsXQeMJAQD1YM55AlhPipbKxOiJNLBDPrK2lt4HwglS5JBNPxJWZ2d1KyisgGVMa69QYGE+qHOUXCFZWTizOARyWNB1618rEVT6Bue09JBvyuCgwoogfjKEKZlheztB0IjEhSmJldKulWQmXm1QhxVO9aWiHFEmMKsJekcQSl2yauGLNiZrOLNrMmgisqJxevEwJTXwYWxVxu6ZiZJWk6ApY6VJwGdCJ7IcN8FcdxFDEINoMRhES7xVqlXg8MkrQcwRxnwPYxM/7hwN5JCzSzmYS/mXLgEeByQokZA17L0b9kcxk2VlxRObl4hSLf0OMGfz9CLaZtCQG/LYA/E/Y8/p2vDDO7P98x8uAMYGAsv/FK0imoMjGzlySdCNxIKC0ioC8wGzjRzN7IZ3xJv6/jfApWLbkQmNkFkt4DtgSuIFQwLsug5caKe/05jQ5JIwirjYsJZT5SHmprAW8Bt5R4YtMGyb0nqTWh5MZawHRgeBJOMtFLM3UjEdX/sCn5CrjRBH5AZpyYU1h8ReU0RjYFLiI4ARjBMw4z+1906/4nYcVVb2oR5JtO0gG/Rc29l8JCGfhcZqv60LUAYzZKknLjd+qGKyonJ5I2JqxuMpOMXmFm3xRA5BLgFzMzST8A66W1TSVkHM+XXEG+6SSqVMysR5Lj1QZJKxNcw1OZKX4EhgH9EnBy6Ex02Y7elFPLKX9hJk0x6W5D46Y/p0Zi6p13CPn2XqcyyegeBLfm3cxsbMIy3yeUZn9I0mCCovi/+PhvQpXfTZOUWc5I2gh4A1if4O2Yip3aEPiOkEOx3j84JC0Gdjaz99Of5zvvxkhTTrrbkHgpeicX1xAcGDYys/3N7EQz259QbuN/wNUFkPki4Vc6hNxwexKyqM8gxOaUfJoaSZtKekTSFEkL4+NDkjYpgLg+BOvJTmb2azPbORb325lgVu2T5/izgDXj85r2qMqBawjX+Fsz28bM9jWzbQhOP82Baxt0dmWKr6icGpE0EzjZMkqLx7YjgDsKnc5I0u+Awwg3wBfNbGgCYy41UdUm+DfJgF9JWxF+lTcn5N5L/Srfk5DiaFcz+zRBefOAU8zsgSxtxxJWr/XOwRdXvVsD7xM8M98Efqqmu5nZwfWV1dBEU/T5ZnZflrbjgOvMbM2q73TywfeonFy0IPxizsZMipCk1sxGkHBQKjCesKJ4n9oF/ybpqXY1IUFrVzP7PnVSUjtgaGzvlqC8CoK5LxtTY3s+nAzcAmxO+Bw3JSRqzUap/zIu96S7jRJXVE4uPiUEx76Upe0EoFSLGB4HfJ32vJg30N2BnulKCsDMpkm6ghDjlCSPAT0IqX8y6UFI9FtvYib7A2Cpq/oB5bpHBXwJdAey1RI7NLY7CeOKysnFDYQUSm8THBlSSUaPJCQ7PTQJIRmJP8eTo+yGmW2Uj7x0M1gDBP+uQKjwm405sT1JhgPXSfoPQSmlvsPDCA4WvST9OdXZzJ7PQ1ZX4LM83t/YKbuku6WA71E5OZF0AnAVlRvmEExJFyR1k5d0H3B5LGR4PzlWOGZ2bBJyo+y2wIpmNjVL2zrAHDPLDMrNR94HwHdmdlCWtmeBdczsdwnKW5KjixEcBCChgFxJHYEuVLpvv2FmZbHakHQpIeluyuydSrp7lZld1VDzKmdcUTm1QlIzYAsqk4x+Zma5boAlQSyfvsjM/pql7T5gOTM7OkF5hxBWNu9TdZW6PXCwmT2boLzOuXtVYmZv5iGrGcF0eTyVyg9CbFx/MzulvmM3JiStQsjysSqln3S30eOKymlUSFqBkMz0GjMrRBaFbDInAWebWZW9mqhUbjazOpUFqYXMHgSnibXTTk8BepnZw0nKKiaSLgCuJFRPfpBKj8ajgb8Trs/LtTt1whWVUyOSegPrmdnfsrT1J5iwrkhY5izChnxeyVLrIG8h8EczG5alrQsw2MyqZHBPQK4I2eFTq9RxhSyjEvdVdqLSHPde0iUpJH0BPGdm52VpuwH4c6kHaxc4y4eTBQ/4dXJxOMFElY33SKCIYRbeJHjGFYv/EVyrs7ElVUvF54WkSySta4GxZjbczD6PKaPWiVVkE0XShYQV28sEJ4DBwBRJ5ycsqj0h4342XiGZ9FcNRszyMYYQ+NuJYNLsRAj0HS3p1w04vbLFFZWTiw2Ar6pp+4bC3Hh6Az0k9ZLUSdLKklZKPxKW9wpwUeZNJr7+J9nduvOhN8HbLhvrxvbEiCU+riS4VB9MiB87iHDd18Sg36SYAXSspq0TCSv9BqDQWT6cLLjpz6kRST8Cp9WQmeJOM1slYZkpJ41q/ziTLBUhaT3gQ2BlQnbxSQRFsichw8KOSeZvi9e3U7ZYI0l/AF4ws1YJyhsDfGBmx2VpuxfYzsx+k5CsuwnK8Egzeznt/H7Aw8ATZnZSErIagkJn+XCy43FUTi5GAGdKeizdyy96d51BMP8lzeUUMQDXzCZL2p5QFO+P8fiB4JHXOwklFYsLdkk7dYKkfTK6tSIk3/08X3kZbAKcW03bYyRrvr2A4A33QtxrTHk0rkwIHr8gQVkNQaGzfDhZcEXl5OIqwkbxR5LuIaSJWZ/gftyJAtQiMrNLkx6zFjInAUmawDLpSqVJzwhZPbLxDXB6wrLnsqx3YTrrUrVwY70xs+lR6fcgKObVgNGE1FAPmtmCpGQ1EAXN8uFkx01/Tk5i1oLbgF+lnZ4InG5mLxRA3r2EWlfjs7S1J6xyqpixGjOSWgAtCbFFMwhmxQ8zui0sxI1c0mOEVc4+6cluJW1JSI31tpkdkYCc5aKcr82sunx4JY2kwwkZ/SdSTZYP0lZVeWb5cCKuqJxaI2lzYA3gh6RrUGXIqWkPZzvg/Xz3qNKVYS2q/ZolWOE3KtupZvZzUmPmkNeBUFNsTcLqJnVz3Zpgxto5iezw0d1+IbC/mWXLhVfyNESWD8dNf04dMLPGkMOtHSFdTb50JQSlQu5qv0lX+J0IIGldYDcqY3HeNrMpScqK8sZL2gY4m2CO2yTKux64xcx+SEiOSZpIgTOIS/oTYX+0PXEv0czuKqTMNBI3dTu58RWVkxNJbYB9CDeGzMBXM7NrEpDRjcrSFj0IxRMzb6CtCEGW481s13xlNhSSmhOU5EksWz7kF0L6ob8nlZ4qOr2sD/xkZtUlwk0MSWcARwB/MLN5BRi/K8EzU4S/j9UIYTaXmdnlSctzGge+onJqRNIOBKWxRjVdjBD8mC+dgAPTxuwMLM7os5BQVuTsBOQ1JJcTajj1AwZSWRr+cOBUYDZwUUKyRHDQOABIfD8xCx2IJe8lvU7whEv/NWxmdk4e419I8IrcP64UVyakajpb0pUNkX8yKs+tgaFmNqbY8psCvqJyakTScGAl4ETC/kaVgnhmlqlQ8pVZ7R5VIYmFC7OtGjGztxKUMwV4pJo0QzcCR5jZugnK+5YQ+1NwRVWbPZx89m0kTQXOMrOBaec2AcYBmxY6Q3tMYCwzOzy+PhpIxVSlUnEl9rfiBDwzhZOLrYF/mtl7ZlZhZoszjySESBoqqVN8eSyVRQ0LTkxb9BohxdC7wBtpx7D4mCSrUH2aocGEmKMkuQ84OTo7FBQza5bjyNe5oB2hOnM6qdfFCLTdmWBhSHE+8CSwIcEF/59FmEOTw01/Ti5+pDjBt52pvEHfS7gh/FgEuQC3A9sC/6CaVWPCjCT8AMiWHf43VHVbz5fZcdxPJT1PdnPcbQnLLCSZq7ZimoXaAd/B0lplWwDHmdm3ku4k/O06CeOKysnFXYTg1EEFlvM/YCtCJgxR3JtPZ+A8MyvWTeYM4BlJs4HHzGxOLN74F+A0Kp1KkiJVVmNdwl5gJkaIk0sMSb9j2ezib5rZiISGv17ST+ni4uNNkmamnTczOzghmSnmE5x6IHhsziP88CA+b5uwPAdXVE5u5gFbSnqPYPKossoxszsTkDMI6B/3aAx4TVJ1ZkUzs9UTkLl0PCDvOKI68AahOmw/oJ+kBVTuiy0EhqVZ6ZK41g55vr/WRA/RgcB+LFs40SS9CPzFzPJJM/Qtwcy2Ycb5iVS9zkL82BkNHB33bk8kVC5OydmA8IPLSRh3pnBqpNCb42lyViB4wm0JHEcoR1FtfI8lW4q+L1BhZkXxJpR0P3W4iSZ5rYUmmr+OBy6mqkfj5cAAK+Eqv5L2IPxgWwH4GdjTzIbHtoHA8gVYxTV5XFE5NRLr79SImSXq+FBsrz9JuxP2FgbFI9uqcXQx5lJIJG1GMMetAdxjZt9L+hXwY56rnHQZ0wkVka/O0vZPgsdedaEOJYGkDYHtgP+a2Tdp508CPkrQxOlEXFE5TZ6MVWPmfwhR4qlwYg6+AcBRVO7/7WBmo6JzxWgzSyRuS9IcQnXm17O07Qk8bWb1ricWi0q+ZGYjc/TbDLjDzPaoryyn8eB7VE6jRFIrQoaK1ArgZDP7UtJBhBtrdcUc60ODmNYkbUVwbqhSe8rMHkxQ1CVAd+Asgvt7ehmRlwmm1qQCjN8kfGdVFBXwe+A/eY5/KXCBpNPMbEAN/VaK8ygYktYk+3dXzP3OJoErKicnsZbSqWS/qZqZbZqwvHYEh4NNCWVF1qPSm+pPwN5Az6TkWZYieIUkevi9QPAag7QkpmndklRUxxBSDPWJ6ZvSGU+yzhYXAU/F1E2PUZkA9y/xODC9QrOZza6HjCkEx5sdCUU9F+U/7doRY9EuJfx/WLWabiW7+m6seMCvUyOSuhACGdckuI+PB6YDvyb8/RRiH+l6QvDmtlFOuvfYGxTwl7KkTpJ2i95rheJSgqt4Z8K1HUgo+/EE8BWwfcLy1qH67+lnIMlrHUVQfBcCHxEU1UeEgokdgP8SqianjvpwNHArwevurZjct1icREjhdQvhu7s6Ht8Rvrvq6ow5eeArKicXvQmOBicDi4AL4t7G1gSz0ZMFkLkfcI6Zjc6yAkitsBIllhG/iuChBrADMErSk8CrZtYvQXH7E/Ijvhtff2dmo4Chkh4g/FpP8oY3lbA6zZZhY0uSdc0vRnXmX8zsHEkjgHuADyUdVqTURScBVwI3EipCPxP/P1xF+HyrK1Dp5IErKicXWxAKxaVuPs0heMFJupyw//FcwjJXJCikbLQkYUtALIY3AHiWoHzTldL7wKEZ5/LlV8CXZrZYUuaK5lGSNftBuK6LJP2Hyv0pk7QeYd/q0aQEWRGrM5vZ45I+AZ4mxN2db2a3FljsxoS/iSXxWCHOZYGkmwn/V5JI0uyk4aY/JxctCDFGSwiVadultX0BbFYAmV9SvXlvNyDpuli9gLvN7CCqpsD5nOSvcQaVe25TCD8GUqxM+MyT5FJgJsHs9i7hR8c9hM9xFmElWZLEGmk7EPb8bpL0SHTEKRTzgeYxyHcaywYez2bZ/x9OQviKysnFBCrNGZ8TVhepLNwHEP6zJs19wFWSJhM25AGWizWrTgPOTVjepkB1pSdmEmoeJcmHwOaEMvAvAL2jW/ciwq/xRPf9zGxmTGl0JrAvIeHvQoL56pZ860ZFl/F7zGxKfJ5jOnZFPvKyDDgHOEhSL4I5bgvCPmch+JxK5TQCOCuuVBcREtQWNHt7k8XM/PCj2gO4Gegfn3e/uPmqAAAgAElEQVQnmDs+J6SSWQxcVACZIpROWEK4AaQeFwP3FkDeD8Ch8XnzKO+38fVRwKSE5e0JnBSftwPGUGlK+grYrKG/9zpezxJgx7TnNR2Lk5JVTfsfCGmMfs5XVjXjnwLcEJ//Fpgb/y4XE5T//zX091GOhwf8OlWQtLxFl19JLYFWZvZTfH04IR2OAS+aWf8CzmMXwgqgHcHT8GUzyzcOJ5ucJwh52joTFOIiQuaBjwlxQV+bWY+k5abJF8Gj0oCxZvZLAmPWJcGumdnx+cosBrXJWhKzbTwFbGcFDtSW1J6QRNgITjef53iLUw9cUTlViGlwHgHuM7OPiiRzJYKH4ctmlq38RSqzwb7AJZZgmfOYxWAEYWX1NMH9eACwDbAJ4Yb3TfUj1EnWCsBhhODXdeLpyQSF+IQlFBMUb+gLqF3JEjOzeps3Y1zYpRThu4uKYUquz0nS8sC6ZjaxvrLSxmqwv00n0tBLOj8a30HIKpAyZ4wiBjcWWGYvwg27dQ19WgOTCK7rScvfChhCMBktAX6Jn8NWCcrYhpDlezFZTGKEkvFbJiRrJkFRPUa4gaocvrv4OVVr+iv16/Mj++Fef04VzOwPhODMywjeabcBUyQ9KmmvAontRtgLqzY5amzrBySWnVrSzpL+AqxFiG9qC6wPtDWzP5jZmITkrERIeNsWOI/gwNGKcIPrRCjauArwgqQkKtW2I6SgWiXK/U7S1ZI6JjB2JsX87gpepTgLDfK36aTR0JrSj8Z/EPZu7iO43y4mrAouAzokKGMWsHct+u0FzExA3koEV+3FacdECuTIAPydUNtrixr6bEnYnD8jYdm/IqQ2+jJe5zuErA4rldp3Rw5nigJ9d0X92/Sj6uErKicnZvamhZpIaxNqDU0g3vgkZUs+Wh9aEExVuVhIZZHBfLgA2BF4iODyfgvBDf2eBMbOxn7A/Wb2aXUdzOwTQrDvn5IUbGbfmdmVZrYJ0AUYB9xJcoHFxf7uir2xXuzrczLwOCqn1lgwb9wv6RlCLrdzCTe+JPieYA7LlQZnU5KJ3eoG3GRm56dOxJQ8AyWtbmZValLlyRZAbTwkh5J8KXoAJK1PcOLYneCGPymhoYv93T0rqbZOIjnrqdWCYl+fk4ErKqfWxP2pHoRA35aEmJ/7Ehr+LaCnpHvNLGsJ+pj37yRgWALyNiSkS0rnFcIeyHpkKZ6YJ6sScu7lYioJBhhLakFIenssIcZoGmEVeZ+ZjUtITLG/uy+pf0Lb+lDs63MyaWjbox+N+yDkNruSSm+1OQTltHvCcrYlxC89DbTL0t6OEBvzM7BNAvKq7HWQEeyb8PXVam8F+B0JBKoS0grdSUjXtAB4nOD916wA11a07662n2OpXp8f2Q9fUTlViF5nhxJ+he9CWGW8S8iM/ZiZzU1appn9V9LpwO3AnySNJOyFAbQnlL5YjlB/KKnYro0lpe89pIJDN5G0TNCtJVOKPlNeNjZJQA6EuLDZBAX1b4LCAtgyxBcvSz7X10DfXdEo9+srBTzg16mCpLkE1+mUmWiAmX1RJNm7EBwduhJctwEqCHs315jZu9W9t45ylpB9U14Z5xMpRV+DvKzyE5JHLWQmIi/KLPh3V5vMFIWiWH+bTlV8ReVk41VCFvGXrBqbfKEws3eA/4sVYteIp6dbyN6eJMUuP1/u8or53TUI5X59jRlfUTmO4ziNGo+jcuqEpJmSZro8l+fyii+vqeKKynEcx2nUuKJyHMdxGjWuqBzHcZxGjSsqx3EcZxkkrS/pX5LeljRXkknqUk3fIyR9LGmBpEmSro0FVzP7tZP0gKTpkuZJ+k90+c+JKyrHcRwnk40JlbznEuqyZUXSUYQiq8MJmU+uJtSvuz+jX8s4TmfgdEJarznA65K2zTUZj6NyHMdxMnnLzNYCkHQA8OfMDjG/4Q3A82Z2Sjz9hqRFQH9Jt5jZiHj+OEJi5u3MbFR8/5vAWIJy27emyfiKynEcx1mGWgYx70Qo/fNAxvlHCLkR04tIHgiMSSmpKGMh8Ciwl6S2NQnyFZXjOE4ZUpv4LjNbJQ8RW8bHTzLGrJD0dVp7qu8bWcYYTcixuRlQbVosV1RNmFYbHF7ntCQtf/WXer93/reX1fUtVFai+KIoKVRcnssrnryOVbMD56CO/+9m1XX8OrJ6fJyRpW1GWnuqb3X9yOhbBVdUjuM4ZUieq6U6iarl+ZqUbI0K2BWV4zhOiRDy4TYaUsVFV6dqodHVgPEZfbOtmlJFQrOttpbSqK7acRzHqZ5mWq7WRxH4ND6m70UhqTWwEcvuXX2a2S+yFaEg6+c1CXJF5TiOUyJIzWp9FIH3gO+BozPOHw4sT6iInOIZYCtJ26ROSFoh9n3NzGbXJMhNf47jOCVCturMBZR1SHy6Q3zsLGkNYJ6ZvWxmv0jqBdwv6XbgSYL33nXAk2b2XtpwAwiBwE9LuoBg6jsTWJdQTbxGXFE5juOUDEU1gj2R8frS+DgR2BDAzB6QtBj4B3AiMB3oC/ROf6OZLZC0ByFA+C6gJTAK2MvMPsw1EVdUjuM4JUIxnSnMrFbLNzN7GHi4Fv2ymQlrhSsqx3GcEqGRef0VjXpftaRnYwbclWro87qkGZJa1FdOUyNmK36tIWT3/OvevP3CVcz88kH639RzmbYuu27BR0Nv5Mdx9zN44EVssN4aS9tWWGE5+t5wEtM+HcD4kXdxxgn7JTKfmTPncOqpV7HNNofQtetxDBo0LJFxXZ7LK0V50Oi8/opGPldzL9CNsBF2T2ajpPZAV+DOmNPJqR1/I0fwW6GYOu0nruvzDHt23ppWLVdYen71VdsysN/ZnPKP/rz42ih6n9Odh+44g84HXALARWcdwsYd1mbTXc6g3ZorM3jgxYz9cjKvvvlxXvO5/PK+LL/8cgwf/hBjx37DSSddTqdOHdhkk/Z5jevyXF4pygNfUdWHlwiuicdW0/5XQASF5tQSM/vMzMY2hOznBn/AoCEjmfHT3GXOd9t3B8Z+MYmnXxzBwoWLuPKWp9hq8/Z03GhdAI48eHeu+dfTzJw1j3FfTeG+R4dydPff5zWXiooFDBnyDmeeeRRt2rRi++23YI89duS557KlC8sfl+fyGrO8FI3MPb1o1PtqzOwX4CFgF0kd09sUfCiPAT42s1GSWkm6RdJoSbMl/SBpmKSuGe/bOBboOl1SL0kTY9GudyTtmDkHSbtKeimaF+dLGifpsow+h0n6MBb1miVpkKTMALWHYzGv30h6U1JFLAB2emzfJ45REa9hz7T3HhPnvH2W+fWJMldMO7d/vPbZ0XQ6RtKpae3LmP4k7RnHP1jS9ZKmxjGHSNok1/eUBJt3XJ/RYycufV0xfyHfTJzG5h3XZ5WV27Du2qsxZuy3S9vHjJ3IZh3Xz0vmhAmTadasGR06rLf0XKdOHfjqq29reJfLc3nlKS+F6vCvnMhX7Q6Ijz0yzv+eEJmcWk21AtoC1xDqmvwV+IpQNKtLlnHPBnYFzgCOAlYCXkpPBS/pz8CbwBqx35+Am4D10vqcAAwEvgUOAU4BNgXeyVSuQOvYdyDBpPku0EfSlXHcmwhp6+cDz0laK75vIDAtjr0USW0IyvpBM5sbz50CDIpjnBDl3A3UxlZwI7AWoa7LyYR4hedjTZgqSJqZ66iFTADatG7J7DkVy5ybPaeCFdu0ZMXWYftxVlr7rDnzadumVW2Hz0pFxQLatm29zLm2bdswb978vMZ1eS6vFOWlaKorqrx23MxsnKR3gaMlXZRWw6QH8DOhLglmNoNwYwaWFtwaQlBmpwLDMoaeDvzZzCz2nwa8A+wDPBFXbLcB44Bd4uoOYCjQP03G1YTo6YPSxhoGfA1czLKukq2As83s5djvXWB/4Dygo5lNjOf/B4wkKhkz+1lSf+BcSeeY2U9xvKOAlYE74/tWBq4HhphZepGw2jpOjDKzHqkXkgz4N/Bb4INajlEv5lUsoO2KGf8pV2zF3HkLmFsRth9XWrEVPyxctPT5nDz/w7Zu3ZK5c5dVjnPnVtAmTwXo8lxeKcpL0axZeTlJ1JYk1O69wPrAXrB0JXEI8JyZLU1UKKm7pOGSfgR+IRTW6gJ0yjLmCynFEhkdH1Mrj82BDYB70pRUJlsAawKPpI9lZpMJJZG7ZvRfRFCeqX5zge+AT1NKKpLaP0pfBfUlKP30/bqTgaFp+027AW1i3/rwfMbrzM9kGcxslVxHbQV/9sUktt58g6WvW7dqwa/bt+OzLyYxc9Y8pk77ia02r5zGVpttwNgvJtX+yrKw4YbrsXjxEiZMmLL03Oefj2fjjTeo4V0uz+WVp7xKmtXhKB+SuJrHgHlUmv+6AyuS5kQh6XDgcWAC4Wa+EyEtx6uElUwmmZl4U16DLeNjyje6prthKivv91naplI1k+9sM1ucce5nqmb1/TljLpjZFOApoKcCuwG/Ae5Ie19t5lwTuT6TvGnevBktWixP8+bNlnn+/OAP2Lzjrzhg3x1p0WJ5Lvz7QXwy9lu++Dr8J33kqbfodfqBrLJyGzputC7HHrEHDz3xVl5zad26JXvttTN9+jxCRcUCPvzwM15/fQTdumX+vkgGl+fyGrO8FE3V9Jf31ZjZHEKOpwMkrUJQRJNIW50ARwBfmtmRZva8mY0ws5GEfav68EN8rGnHPnVjXztL2zpUvfHny23AJoSV5SmEzyB9FVSbOTcovc44kJlfPsh5p3bjiIN2Z+aXD9LrjAOZPmMOh/e8hUvPO5SpY+5hh2025ujT+ix93xU3P8k3E6cx7p0+DHn8Em7t+0LerukAvXufzIIFC9lll6M455wbuPTSkwvq+uvyXF5jlgdNV1FpWQtbPQeRfk9wbLgBOBe42swuSmt/Hvi1mW2Zdm5r4L/AeDPbOJ7bGPgSON3Mbk/ruxzBNHexmV0Z96jGA3OBbbKZ/+Ie1eQ4/s5p59cBvgGeMLNj4rmHgX3MbI2MMT4BvjezdC+/1FxuMrNzM/qPBCqA3wFXmNmVaW0rE5TXO2b2xxo+y7eBBSmZ0cPwVeD/zOyFtH6pz+romMKkztSnSm8+1KfCr+OUL3Wv8Lvelr1r/X928ieXlY3rXyJq18zeItw0Uzfu+zK6vAhsEd2195B0EjCYkNywPvIMOJ2wvzVc0pGSuko6TtLdsc9i4AJgJ0nPSPqTpCOANwhms8vrIzsHtwG7E+LH7s6Y8yzgfGBvSS/HPbs9JJ0i6doCzMVxnDKjqa6okrya+wg36LfM7OuMtv6EbLrdCIHCPQnefu/UV5iZDSI4RMwiZON9iaAIJqX1uY+QOeNXhD2ku4AvgF3N7Kv6yq6BJwlFwJ40s2lZ5nwXcADB5Hk/wTTYk3oqbMdxmhbNmjWv9VFOJGL6cwKSDiU4l+xmZsMbej65cNOf4zQkdTf9tf/N1bX+Pzvx4wvLxvTXNJ3yE0bSpoT6LNcBb5eCknIcp/QoN5NebXFFlQwDgB0JgcA9GnYqjuOUK66onHpjZrs19Bwcxyl/VGaBvLXFFZXjOE6JoCaaQqlpXrXjOE4JEkJImx6uqBzHcUoEN/05TY753z7a0FNwHKcOuDOF4ziO07hx05/jOI7TqGmaCypXVI7jOCVDs6apqVxROY7jlApNU0+5onIcxykVzPeoHMdxnEZN09RTrqgcx3FKhmZNU1MV1OIpaVtJD0qaKGmhpNmS3pN0Qax4W5exWkjqL+l7SUskDYvne0k6oI5j3S9pQl3e4ziO0+BItT/KiIKtqCT1BG4HxgBXA+OAlsDOwBnAWsBZdRjyZOBEQqHB0YSCiQC9gGfjUVuuAP5Vh/6O4zgNT/PyUkC1pSCKStIuwB3AC8AhZrYorXmwpBuAXeo47BbAHDPrl8e8WpjZwiwViB3HcRo/ZbZSqi2FMv1dACwBTs5QUgCY2VwzGwIgqZmksySNkbRA0vRoLlw71V+SAScAbSVZPHrE8ysDf007f398z6Xx9W8lDZI0Gxgc26qY/mL/DyT9JGmWpPdjxd70PofHMY/POH+rpJ8l7ZB27teSHpP0YzR7fiLpbxnv6xLHO1TStZKmRvPoa7EYIxn995P0pqQ5kuZJej1dpuM4ZY7qcOQjJmzbPCtpSrzXfBa3WVpk9NsrbufMl/Q/Sf0krZKf9KokvqKS1BzYAxhpZlNq8Zb7gMOAm4BhwHoE09wwSduZ2TyCufBi4PfAXvF9X8fzrwJvxfcA/JAx/tPA/cCtQPMa5rEBYRX4HeFz6Qr8W1JbMxsAYGaPSuoK3CbpfTMbI6kbcCZwlpl9ED+DtYHhBGV9HjA1XmM/SWuZ2ZUZsq8D/gMcT1C81wGDJG1mZovjmMcB9wADgRvjtZwNvCnpd2Y2Jn1ASTNruFbi9eTq4jhOY6IIzhSSOgHvELZr/g5MJ9zTrwI2B46J/boALxG2XS4C1iXcu7aUtLuZLUlqToUw/a0BtAYm5OoYTYTHAKea2Z1p5z8CRhGq5d5hZu9J+gFYbGbvpQ3xg6TFwA8Z59O528yuyjUXMzsuTX4z4PV4LacRKvimOBPYCXhC0sEERfu8md2a1ucsoB2wjZmNjudejg4kF0q63czSFcloMzsmTf4vwOOEqsHvSmpDUORPmtkRaf2GAF8AlwDdc12j4zglTnEsf38h+BMcnLZNMlRSe+BwScdHS9n1wCfAYSmlJGkqMIRwP3osqQk1dJzzfoRVx0BJy6UOwsVPBjonIOOZ2nSS1FXSK5KmAb8AiwgrnE7p/cxsPnAosD7wITCHquXnuxKUz+iM8w8CrQiKLp3nM16n3tc+Pu4MrAI8nPE5/UxYhVb5nMxslVxHtR+G4ziNEmverNZHHqS2a2ZlnJ8V2xZLWg/YAXgofeVkZq8S7t0H5zOBTAqhqKYDFVTeZGuiXZzDj4QPIP1Yn7CiyZepuTpI2olgQhRwKrAr4Uu4l/DLYhnM7HPgXaAFMMDMfsroshrwfQ1zWT3j/I8ZrxfGx5TsdvHxOap+TkeSzOfkOE5jpw57VJJm5jqqkfIQMAO4S1IHSSvFLY6/AjdFxbRl7PtJlvePSWtPhMRNf2a2WNJQ4I+S1jGzmhTFdMKKajcqtXg6c5KYUi36HBbl/5+ZpZQEklbI1lnSKcCewAfAPyQ9aWafpXX5EVg7y1vXSWuvC9Pj4ylRpuM4TZEieP2Z2bfxx/uzwDdpTVeb2cXxeerH9owsQ8wAfpvknAoVR3UNwax3p6TuZvZLemPcc9kFeJkQB7W2mdXKRJeFhQRzWj4Ywdy3dAkraS2gW2ZHSdsCNxNixHoBI4HHJe0QzYIAQwkKbOsM899RwHyguv206hgOzAY6mdlddXyv4zjlQh2cKepr3o97UYMIVqEDgZmE7YULJC1JU1ZQ/UIgUU+tgigqM3tH0mnAbcAISf0IHiQtgN8RgnYfN7Ozojv5g5JuA94GFhA8/7oCL5nZkznEfQJ0lvQnYBow3cwm1HHKLxIcIB6R1J+wGro4jtc21UlSW4KTw2fAuWa2MLqwjwD6EAKSAW4hLJNfknQJMIWwr3UQcGGGI0VOzGyupL8D90hajbDv9iPBJLg9sDDjj8dxnHKkOM4U1xLue9um/fgeprCau0TSACqtQpnbGBC2PrKttOpNwTJTmNldkkYQFMDFhJvqAsJN/l9AamVwHGGFcSLBFXIJYTPuTSqdCmribKAv8CRhT+cBqjo35Jrr6zHG6XxCkPJEgjt7O6B3Wte747l9UybC6KJ+BnC3pKFm9qiZ/S96NF4L3ACsSPDOO9HM7qnL3NLmeJ+kbwnu7gMIq8hphBVd3/qM6ThOiVGcgN9tgc/SlFSKkQSfgk5U7k1tSfDyS2crgnt7YshjaZo0/uU7TsNRZ62zcfeHa/1/9qsnjqqXVos+BpsBG5lZRdr5K4F/AlvHH+gfEBTXDmnu6X8AXgMON7OB9ZGfDc+e7jiOUyoUZ0XVh7C98IqkWwlu6V0IFqfX0pIL/IOwmno0bpmkAn5HAE8kOaGGjqNyHMdxaksRUiiZ2bOEDEALgTsJYTEHErL/HJDWbyiwP7AhYZ//5vi4byqjTlL4ispxHKdEsCLVozKz1wgmvFz9BhNzqBYSV1SO4zilQhPNnu6KynEcp1RomnrKFVXT5ouGnkBBabVB79ydEmT+t5cVVZ5T6nSs+1vyy+FXsriichzHKRV8ReU4juM0aorkTNHYcEXlOI5TKriichzHcRoz1jT1lCsqx3GcksGdKRzHcZxGjZv+HMdxnEZN01xQuaJyHMcpGZpoZoqS08+SekiytGO+pLGSrpW0cgHlTohFHuvyntRcN6xunNr0cRzHAYLpr7ZHGVFyiiqNY4CdCdl7XySkoB8iqVDXlMoeXBdeJMxxahFkFZyZM+dw6qlXsc02h9C163EMGjSspOT1/OvevP3CVcz88kH639RzmbYuu27BR0Nv5Mdx9zN44EVssN4aS9tWWGE5+t5wEtM+HcD4kXdxxgn75TWPFKX+ebq84soDMKnWRzlRyqa/MWb2UXz+uqQ1qVRew5MWZmb/rcd7fgB+KIasYnD55X1ZfvnlGD78IcaO/YaTTrqcTp06sMkm7UtC3tRpP3Fdn2fYs/PWtGq5wtLzq6/aloH9zuaUf/TnxddG0fuc7jx0xxl0PuASAC466xA27rA2m+5yBu3WXJnBAy9m7JeTefXNjxvV9bm88pYHwHLlpYBqSymvqDIZER/bA0jaTNITkqZLWihptKQj0t+QZnb7vaR+se+Pkp6WtG5G3yrmOEnrS7pH0iRJP0uaLOnRlAkym1mvNtRgHsw5z0JRUbGAIUPe4cwzj6JNm1Zsv/0W7LHHjjz33BslI++5wR8waMhIZvw0d5nz3fbdgbFfTOLpF0ewcOEirrzlKbbavD0dNwof7ZEH7841/3qambPmMe6rKdz36FCO7v77Rnd9Lq985S1Fqv1RRpSTovp1fPxB0tbA+4SCXqcD/0dQZI9IOj7Le+8D5gNHEKpWdgEerkmYpF8BIwmmx+uAfYBzgF+AFfO7lGqp9Twlzcx11EXwhAmTadasGR06rLf0XKdOHfjqq2/rfTGNRd7mHddn9NiJS19XzF/INxOnsXnH9Vll5Tasu/ZqjBlbKXfM2Ils1nH9vGSW8+fp8gpIE92jKmXTX3NJywFtgL2BnsBk4D/A88B0oIuZzYv9h0haA7hK0n1mtiRtrBfM7O+pF5JWBa6XtI6ZVbe/dDmwKrClmX2Zdn5gEhdXDfWZZyJUVCygbdvWy5xr27YN8+bNL3l5bVq3ZPqM2cucmz2nghXbtGTF1i0AmDWnYmnbrDnzadumVV4yy/nzdHkFpLz0T60pZUU1MuP1COCE+Lwr0AdYGJVZipcIpZQ7AZ+lnX8+Y6zR8bE91TtC7AO8lqGkCk2t52lmq+Qe7gurreDWrVsyd27FMufmzq2gTZ437MYgb17FAtqumHHTWbEVc+ctYG7FQgBWWrEVPyxctPT5nDxvSOX8ebq8wlGsCr+NjVI2/R0J7ABsC6xlZjuZ2SfA6gQFfDawKOPoH9+7RsZYP2a8XhgfW9Ygfw1gUr1nXz/qM89E2HDD9Vi8eAkTJkxZeu7zz8ez8cYblLy8z76YxNabV47bulULft2+HZ99MYmZs+YxddpPbLV55Qb5VpttwNgv8vvqy/nzdHkFpIma/kpZUX1mZiPN7KPoXZfiJ2AJMICgyLIdSXjV/QDkt1FRQrRu3ZK99tqZPn0eoaJiAR9++Bmvvz6Cbt26loy85s2b0aLF8jRv3myZ588P/oDNO/6KA/bdkRYtlufCvx/EJ2O/5Yuvw03okafeotfpB7LKym3ouNG6HHvEHjz0xFuN7vpcXvnKW0pz1f4oI0pZUWXFzCqAN4FtgI+iMss85iQgajCwp6SNExirJOjd+2QWLFjILrscxTnn3MCll55cUFfcpOX1OuNAZn75IOed2o0jDtqdmV8+SK8zDmT6jDkc3vMWLj3vUKaOuYcdttmYo0/rs/R9V9z8JN9MnMa4d/ow5PFLuLXvC3m7phfi+lxeecsDmqzXn8xqvU3RKJDUg+D9tm1aHFVmn60JThVjgH7At8AqhL2pHc3s4JrGktQFeAPoambD4rkJwDAz6xFf/wr4IL7lKiBldtwf+KeZTU4bv4OZTahmnLr0yTnPulH7PapSxEvRO42bjnXWJhte9kqt/89O6P3HstFWpexMUS1mNlrSdsAlBNfxNYAZwFjg8YRkfCdpR4L33z+B1YD/EVZzSazYHMdxlqXM9p5qS8mtqJwk8RVVkviKyqkbdV9Rtb/y1Vr/n5140V5lo9XKckXlOI5TlpSZk0RtcUXlOI5TKjRR058rKsdxnFLBFZXjOI7TqGmaesoVleM4TqnQVFMouaJyypZie+EV08vQPQybKGUWyFtbyi4zheM4TtlSxBRKkrpIGhLLAlVI+kzS3zL6HCHpY0kLYl2+ayUlnnvUV1SO4zglQrMiLS0k/ZWQL/Vu4GZCUu9OwAppfY4CHgLuAv4ObEZIsLAh8Jck5+OKynEcp0QohuUvpoe7C7jQzK5Pa3o9rU9z4AbgeTM7JZ5+Q9IioL+kW8xsBAnhpj/HcZwSoUg5aVNV0G+roc9OwNrAAxnnHyGsvg7OawYZuKJyHMcpESTV5ZiZ66hGzO8JeVEPkjRO0uK0/aeU6W/L+PhJ+htj9Yqv09oTwU1/juM4JUKR9qjWjcdtwMXAp8AewAXArwhFa1ePfWdkef+MtPZEcEXlOI5TIqgOisrMVqmnmGZAW+BwMxsYzw2T1Ao4V1J6HEZ1SXITTXjtpj/HcZwSoUh7VD/Gx1cyzr8cH3+b1ifbymk1sq+06k0iikpSD0mWdsyXNEXSq5LOkrRyEnJqkN9S0qWxkGDJI2mYpGENPQ/HcRoXzVT7Iw/GVHM+NeoSgjkQMvaiJLUGNiJj7ypfkl5RHQPsDOxN8OSUEnUAACAASURBVKv/ilD99lNJ2yYsK52WQG+gSwFlFJNT4tGomDlzDqeeehXbbHMIXbsex6BBw1xeDfT86968/cJVzPzyQfrf1HOZti67bsFHQ2/kx3H3M3jgRWyw3hpL21ZYYTn63nAS0z4dwPiRd3HGCfvlNY8Upf55NnV5ULQV1dPxMfMPbz+CSe8D4D3ge+DojD6HA8unjZEISe9RjckoD/+4pDuBt4DnJW1iZgsSlll2mNlnDT2HbFx+eV+WX345hg9/iLFjv+Gkky6nU6cObLJJe5eXhanTfuK6Ps+wZ+etadVyaZwkq6/aloH9zuaUf/TnxddG0fuc7jx0xxl0PuASAC466xA27rA2m+5yBu3WXJnBAy9m7JeTefXNjxvV9bm84sqD4sRRmdlgSS8Dd0hag0pnijOBvmY2McxFvYD7Jd0OPEllwO+TZvZeknMq+B6VmY0BrgTWBw4DkLSXpOclTY5mwi8k3SxppfT3Sro/ulFuK+k/se9USVdLWi722RD4Kb6ld5r58dK0cQ6RNDK+f5akFyRtXY2srSS9EVOGTJb097Q5fxDPfyJp77T3HhFl7pR5/fG65qRfm6T9ooxZkubF8U5Pa1/G9BdTmZikQ6OL6FRJsyW9JmnTOn4l9aKiYgFDhrzDmWceRZs2rdh++y3YY48dee65N1xeNTw3+AMGDRnJjJ/mLnO+2747MPaLSTz94ggWLlzElbc8xVabt6fjRusCcOTBu3PNv55m5qx5jPtqCvc9OpSju/++0V2fyyuevBTNmqvWR550B+4BzgdeAg4FLgJOS3UwswcIK6rOhP2si4C+BMtaohTLmWJQfEz9b9sIeBs4CdgHuB7oBryY5b0tgGeB52KffwO9gD6xfSqwV3w+gGB63JnwISOpB/AEMJnwYZ8MbAwMl7RZFlmPxf7dgP8At0i6CvgXcCtwEDAHeFbS2vF9T8R5LGOui14yPYCHzWx2PPe3eJ2LgL9FOf2A2vwMu47gNno84bPrCAyKUeLLkEcMRVYmTJhMs2bN6NBhvaXnOnXqwFdffVuXYVwesHnH9Rk9duLS1xXzF/LNxGls3nF9Vlm5DeuuvRpjxlbKHTN2Ipt1XD8vmeX8eTYFeSmKZPrDzOaZ2blmtp6ZrWBmG5vZtWa2JKPfw2a2lZm1iH3PN7P5+UmvSrHc01Pf3joAZtY31SBJwHBgHPCWpN+YWbqNoyVwtZn1i6+HxA27kyRdbWaTJI2MbZPSl5ySmgHXACOBA8zM4vk3CEFplxBsqumy/mFmg2K/4cD+wHlAJzP7Jp6fCnwEHAjcZWaLJPUFLpR0tplNj+MdAawK3BHf1xa4ERhqZnunyX2tdh8jo81s6a8VSb8AjwM7Au/Wcox6UVGxgLZtWy9zrm3bNsybl/jfZNnLa9O6JdNnzF7m3Ow5FazYpiUrtm4BwKw5FUvbZs2ZT9s2rfKSWc6fZ1OQl6KJJk8v2opqmY9X0tqS+kgaDywgrC7e+v/2zj1Oq6r6/+/PDMgwgCBaKCAXI0IUv1RqYpmAmWkZ4i2/onnpIorihUojFb4oUWqhlKmkyU+8J941gryghZKo5IUBvAEiiHKVYZiBcP3+2OcZzhyfYZ5nrs/Ms96+zuvM2WefvdY545wPe+119o5O90tz/b1pjguAb9Rgtx9hmo+7UiIFYGargNnAkET97cDMWL0ygsguSolUREm0j/eCbiHc59mxsnOB58wslQFzKOH7hJupHY8mjl9L4wcQvqGoacvGcHFxEaWlZVXKSkvLaFfHF2g+2ttcVk6H9omXXPu2lG4up7SsAoBd2++wu2v7tmyq4wuwJT/PfLCXorF6VLlGYwlVj2i/MurlzCaMV11HyBA8iBBSA0j+psvNbGOibHW0r+nr587R/sM051alub7UzLYlyraS+CbAzLZGPxbFylYTejcjJRVE41VfJepNRaRSu1bU4Hd1rE0cVyT9aCh69erG9u2fsnTpysqyRYveo0+fHju5yu2lY+GSFRzQf0e7xW3bsE/PLixcsoINGzezavV6BvTf8W+PAfv2oGRJbf+XCbTk55kP9lI0Unp6ztFYQnVstJ8DDCDk3v/czG40szlmNh/4pJpri9J8h9Ul2idf3ElS5/dMc26vDK7Plj8AvQnjbucRxPCh2PmPo33dBhyagOLiIo48chBTptxFWVk5L7+8kKeemsewYclOqdtLUVhYQJs2rSksLKjy86MzX6J/37057uiDadOmNWMvOp43Spaz5J3w0rtrxnNcdsFwOnVsR98vdOWsU4cy/a/P1WCt8e/P7TWevRTeo2ogJO0P/IoQQrufHVNrbE1U/clOmkmubXIKIUz3z+g41bNI9sYWAyuBEdFYWMqnLsC3gKczuIWMMbN/A/OAsYTEjamJHtpcgiCPTHN5zjNu3LmUl1dw6KGnMWbMtYwff26DpuI2d3uXjR7Ohrfu4OejhnHq8Yex4a07uGz0cNas28T/jpzM+J+fzKrXb+WggX04/fwplddd9fsHeHfZahbPncKs+6/k+psfr3NqekPcn9trXHvQqFl/OYViQze1byRk1t1OSEt8i5CksRch9/6HhNDZ98zsP5JaA4uAQkL23kZCQsMgQjbeWWY2LWp3GiFE+BGht/If4CjgEkISw6iYD+8QBGs0sAFYaWYrJZ0O3EEY3/kz0J6QRNEdOMjMFsdsHZccu5G0ANhgZoMT5QbcYGYXJcpHAHcC/wV6mtnKxPmfEsazZhMyE9cAXwL2MbOfR3WeBUjZjGbceAYYbmYPx9rqBbwXf2bZsaRe5+PKd3wpeic7+matJgf/9Z8Z/83++6RvtBi1qu8e1R2E7LOnCOnjfQi59fulMvmiHsb3CVl3fyasELmd6leErCCkcA8jiM1phDTtCxP1fkpIyniC8OX0TyN70wlro3QlfJQ2FXgX+HpKpOqZGZEfDyVFKvJnKuH+2wJ/IaTunwcsS9Z1HMeJk6+hv3rpUTUU1fVychlJxxPEarCZzWlqf3aO96jqE+9ROdmRfY/qkBmZ96hePKHl9Kh8mY96IpohoidheeYXc1+kHMdpbrS0bL5McaGqP24hfCf1Mg0whYjjOE7BZ+agyQ9yWqjM7Mym9iFTkskWjuM49U1LG3vKlJwWKsdxHGcHylOlcqFyHMdpJuSpTrlQOY7jNBdcqBzHqRONmTLemKnw4OnwuYILleM4jpPTtGqs2VlzDBcqx3GcZkKB8vMbfRcqx3GcZoJ/8Os4juPkNHka+XOhchzHaS546M9xHMfJaTz05ziO4+Q0rVyoHMdxnFxGeRr6q7exOUlnSrJo+1qa80WSNkTnp9WXXafx2LBhE6NGTWTgwBMZMuRsHnvsWbeXQ/ZGnvFt/vn4RDa8dQdTfzeyyrnBX9+PBU9fx9rF05h57+X06LZH5blddmnFzdeew+o3b+O9+Tcx+sfH1MmPFM39eeaaPQihv0y3lkRD9Kg2AWcB8xLlwwmr2lY0gE2nEZgw4WZat27Fv/41nZKSdznnnAn069ebL36xp9vLAXurVq/nt1Me4luHH0Dbol0qy3ffrQP33nIJ5106lSf+8QrjxpzE9BtHc/hxVwJw+cUn0qf3nnzp0NF0+VxHZt57BSVvfcDsOf/JqfvLd3uQv1l/DXHf9wOnSGqbKD+LsJR8eX0YkdRKUp6uztL4lJWVM2vWXC688DTatWvLgQfux9ChB/PII8+4vRyx98jMl3hs1nzWrS+tUj7s6IMoWbKCB5+YR0XFNq6ePIMB/XvS9wtdARhxwmFMuuFBNmzczOK3V3L7PU9z+knfzLn7y2d7KQpkGW8tiYYQqrsJPafhqQJJewNHANOSlSXtJul6ScskbY32V0lqHavTKwoZ/kzSBEnvA1uBvaPz+0i6T9JaSRWS3pD004SdwVEbp0q6UdIaSWWSnpTUJ41fx0iaI2mTpM2SnpJ0UKLOtCic+SVJMyWVSnpf0u8ktUnULYp8fyvycZWkP0rqkKhn0fM4T9LbkrZImifpy5E4T5C0QtJGSQ9J+lzGv5k6sHTpBxQUFNC7d7fKsn79evP228vdXo7b69+3O6+VLKs8LttSwbvLVtO/b3c6dWxH1z0783rJDruvlyxj377d62SzJT/PprCXopUy31oSDRH6W0foOZ1FEC2AM4DVwMx4RUntgeeBzsBEYBFwMHAl0As4PdH2xcBrwKjoeL2kPYF/AZ8CPwdWAT8AbpH0eTO7OtHGb4EXCKvwdgF+DTwraYCZrY/8Ohu4FbgXuA4oBC4B5kj6mpm9HmuvDfAIMJWwDP3hwOXARmBC1F4h8ATwVWASMB/oB1wFDJA0xMw+jbU5HHgPGAO0Bn4XPdNZwC7Ajwki/fvI7nASSNqQLEtitrimKpWUlZXToUNxlbIOHdqxefOWjNvIBrdXf7QrLmLNuk+qlH2yqYz27YpoXxz+PbVxU1nluY2bttChXTIgkh0t+Xk2hb0ULW3sKVMaKuvvduAxST3MbDlwJjDdzLYnFv4aTXhhf8XMXovKnpJUBlwv6Tdm9mas/ibge2a2PVUgaSxBcAbG2vibpI7AWEl/NLP4S/t9Mzs5dv1C4EVgJDBJUjuCMDxgZqfG6s0ClhBE9KRYe0XAL83soZj/XwVGEAkVQTiHAt81sydj9VYADwNHE4QsRSFwtJltiWy3I/RGe5nZETGf+gMXSio2szIakOLiIkpLq5ooLS2jXR1faG6v4e1tLiunQ/vES7V9W0o3l1NaFoaMd23flo8rtlX+vKmOL9yW/Dybwl6KlhbSy5SGGpv7O6EHdYakbwJfIE3YDzgGWAAsjMJarSS1Av4WnT88Uf/RuEhFDAFei4lUijsIIchDEuX3xQ/MbB6wNGZrENAJuDPh01bg2TQ+baeqyEDo9cVHVI8B1gKzEm3Ojq5Ptvl0SqQiSqJ90k4JIKBHohwz61TTlrxmZ/Tq1Y3t2z9l6dKVlWWLFr1Hnz6fMV0vuL36Y+GSFRzQf0e7xW3bsE/PLixcsoINGzezavV6BvTf8b/rgH17ULJkRZ1stuTn2RT2UuRr1l+DCFUkJtMJPamzgXlmVpKmahdCOGxbYkvFpPZI1F+Vpo3OwIdpylN1d0+Up6u7OlavS7R/JI1fI9L4VGpmWxNlFYSQYIouUfvJ9jYTek/JNtcljrfWUF6U5p7qleLiIo48chBTptxFWVk5L7+8kKeemsewYUPcXo7YKywsoE2b1hQWFlT5+dGZL9G/794cd/TBtGnTmrEXHc8bJctZ8k54yd414zkuu2A4nTq2o+8XunLWqUOZ/tfncu7+8tleioIstpZEQ37wezvwC8JY03nV1FkDlAI/qeb8ysRxun7vWmDPNOV7xc7HSVe3CzvEcU20Pw94qRq/smUNQQy/t5PzOc+4cecyduwNHHroaXTq1IHx489t0FRct5cdl40ezuUXn1h5fOrxh3H15AeYOHkG/ztyMpMnnMlfbhjFS6++zennT6msd9XvH2DKxB+xeO4UtpRv4/c3PVrn1HRo/s8z1+xB04X+JI0HxgH/MbOBiXNHEsbb/4cwPPMQcGliyKVu9s3q58YlnUkQpy+b2YKo7E+EQf/TU05Hg/wPm9mZkq4kJED0N7P3d9J2L0JywcVmdn3i3CTgUqqOUSFpBmHsp6uZbZA0GHgGeMHMDo3V+xphjGqsmU2KEjw+AKaZ2YU13PM04LhkGC31SzUzRcc/JIQ+DzSzV2po04AbzOyiWNlA4FXgLDObFis/k8Qzz44l+RnwbgH4Cr8tgb5ZB+jOf+GZjP9m/zhoSL0EACXtR0gA2wh8GBeq6L06mzDWfgvQlZCw9i5wWCJJrNY06BRKZlZdTyrFZOBE4HlJk4E3CFltvQjjOueb2bLqL69s4wzgyUj4VgInA8cTxCep6t0l3UcQjj0JWX8fADdHPpdKugi4VVJnwr8O1hJ6XQcCFWZ2RQ0+JbmbkGU4M7rPlwm9w72BowjCNDfLNh3HyTMaO6QnqQC4jZAFPYAwfh/nGsJ7+wcpUZK0ipChfBKJnIDa0qRz/ZnZJklfB34FnE9ICthM6D39nc+G7dK18ZGkQ4HfENLD2xOy835iZremueQy4BuEMbRiYA5wQSo1PWrzdknLCb292whJGasJ/6q4uRb3+V9JxxBS3EcQutAVwHLgqeh+HcdxdkoThP4uBroD3yZ8IlOJpG7AQcCYeM/JzGZL+gA4gXoSqnoL/eU6sdDfcDN7uIndyRE89Ndc8dBfSyD70N8l857O+G928iFHbKypzs6yfyXtA7wOjDCzhyU9C3RKhf4kHUX4NvYoM5uVuPZvQE8z65+pvzujpSWHOI7jtFgaK+tP4YPXPwN/38k/7FOZ0sls5FRZMuO61vgyH47jOM2EbL6PyvZbyQQ/IYzJZ9Ijqq6XV28Rm7wRKjN7lvBxrOM4TrOksKDho/WS9iAkSUwCNktKCV4roDA6LmdHDkG6nlNn0ve0aoWH/hzHcZoJjRT66w50JAjV+tj2dWD/6OfxQGp6u/3TtDGAkA1YL+RNj8pxHKe500hZf28TpqZLcj0hq/rHwHIzWyFpPjBC0vWx9PQjgG7Ag/XlkAuV4zhOM6Ex5vAzs1LCvKZVSK3IEA2jpLiU8M3UPZKmsuOD33nAX+vLJxeqvKZvUzvg1JIty+9pahecJiDXJps1s6clfQ/4P8Kk2ZsIs1T8Is0E4rXGhcpxHKeZ0LoJl/kws8HVlM8ksdZgfeNC5TiO00zItR5VY+FC5TiO00xwoXIcx3FymkIXKsdxHCeX8R6V4ziOk9M01cKJTY0LleM4TjOhdZ72qLKaaUPSlyXdIWmZpApJn0h6UdIvJXWM6li0wm2TIGkfSTMlrU/5Imlw9PPgpvIrGyQtjVYPdhzHqaRAmW8tiYx7VJJGAn8krE/ya2AxUAQMAkYDnycssjUIWFHvnmbO74GBhFV/P4p86dOE/tSG4cAnTe2E4zi5hYf+dkK0gu6NwOPAiWa2LXZ6pqRrgUMBzOzF+nIuWga5MGGvJvYD/m1mlatRSmpWQmVmrza1D47j5B75mvWXaejvl8CnwLnpRMPMSlMrPCZDf1HozaKw4ROSSiWtlXSTpPbxdqJ610u6UNLbwFbCjL1I+pykP0v6UNJWSW9JGiupMDo/WJIRek/HRm2ZpF7V3ZSkkZLeiMKYayTdE68vqa+kTZKmJ647Lmr7wljZKEnPS/o4uscFks6NxDZp97QoZLo5Cp/Ol3RK7HyV0J+kMyN735R0S+TrWkkPSupa3f05jtOy8NBfNURCMBSYb2Yr62DrQeAuwgy8XwOuBHoA303UOwn4ALgMKAXekVRMmCSxG3AFIex4FHA1sA9hNt9XCGHHh4AlhMkSAVYBvdLc13hgHHA78LOo7auAuZIGmtlHZrZE0jnAXZLmmNmtknoCfwEeMbMbYk3uA0wHlgLbgYOB6wiTNF4Rs/trgvDfS1jzpZQQqvyMj2m4HXgMOJXw7K4B7iT8fpL3t6GmxszyM4zgOM2VVnm6MFMmob89gGLCC7gu3Gtml0c/z5a0DfiNpEFm9kKsXhHwLTOrHKORdC5hpcljzezxqHhW1Fu5SNJ1ZrYIeFFSBbA+HoIMqyrvQNJuBCGbYWZnx8pfJgjexQQxwczuljQEmCLpFeBPhPGjs+JtmtmYWDsFwBygELhE0pVmZpL2iez+xcx+FLt8Vs2PD4DHzeyixH1cI2kvM1uVYRuO4zRTCvN0jKox9fneao4PT5Q/FRepiCHAxphIpbgj2g/O0pdBBEG8M15oZgsIySLJtVhGE9ZomQt8BTjFzNbHK0j6ShSKWwlsi7argN0IiSYARxKe+c1Z+pvi0cTxa9G+Z7KimXWqaaulD47jNBGNtHBizpHJ/awBykjzMsySDxPHq6N9chnjdD2Dzmmuj9dNtxTyzuhcjU+pNqu0Z2ZbgPuANsA/kwkj0bjW88DehDDiYcBBwMSoSttov0e0r21W5NrEcUW0L6ple47jNCPydYyqRqGK1hR5GjhQ0l51sLVn4rhLtE++fNP1bdemuR5gr9j5bEjVr67NKu1J+jJhnOklYIikMxPXDCOER08ws7vNbK6Zzeez9/JxtO+epb+O4zguVDUwiTDe8idJnxnXktRO0pE1tHFKNcdzMrD/NNAxWqArzmkEMXgmgzbivACUR9dXIukAYEBkL1W2K3A/8AbwDULCxB8l7Ru7NCVIW2PXFQGnJ+zOJiRajMzSX8dxHAplGW8tiYy+ozKzuZLOB/4AzJN0CyHzrg0hg28k4WU+eyfNnCLpv8BzhIy4ccCTiUSK6rgDOB+4U9IVwCLg24Skhz+b2eJM7iN2PxskTQSukvSXyPeuhDGllcDkWPWphDGm75jZ1iix4yDgPkkHm1k54b63AXdLugboAIyJyuJ235P0W2BslMl4P2FFzAOAXczsN9nch+M4+YVn/dWAmd0kaR5BHK4ghO7KgYXADcBNNTQxnDCjxcWEnsfthJd5Jra3RJl3k4BfEcaYlgFjCSnaWWNmV0v6GLgAGEFIE58FXGpmH0PlbBw/ICRPvBNdt1nSycA8wn2fY2Ylkk4iCN1DhPG3WwnjXbcm7P5K0lsE4b2bIGaLgWtrcx+O4+QPLS2klylq6G9pYt8r7WZmNX7b4zQqLSs+4DjNi6xlZ+aKv2X8N/ud7ke3GFnz2dMdx3GaCT7Xn+M4jpPT5OkQVcPft5mNNzN52M9xHKdu5Gt6uveoHMdxmgmtCzz05ziO4+QwLa2nlCkuVHnNkqZ2wHHS0rbHuEa1t2X5/zWqvUDfrK9woXIcx3FymnxNpnChchzHaSbIe1SO4zhOLuOhP8dxHCen8dCf4ziOk9MoT2emyFeBdhzHaXYoi63WNqQjJE2TtFhSmaQV0erlA9LUPVLSi5K2SPpI0i2S6n31cBcqx3GcZoKU+VYHRgI9CMsdHQ1cEh2/JOmQHb5oMPAk8D5wLGF18+8DT0iqV23x0J/jOE4zoZFyKUaZ2UdV7EqzgPeAnwMnRMXXEBaU/YGZfRrVW0VYLukk4L76csh7VPWIJMtw65VhvfGxtpcmzn0o6R+ShjbW/W3YsIlRoyYycOCJDBlyNo899qzbc3uNZm/kGd/mn49PZMNbdzD1d1UXyR789f1Y8PR1rF08jZn3Xk6PbntUnttll1bcfO05rH7zNt6bfxOjf3xMnfxI0djPE6BQmW+1JSlSUdkG4C2gO4CkboQFZKenRCqqNxv4gB1iVi94j6p+GZQ4/i3h8/PhifJVaeqmaEVY0bgboVsd5zngUsI/rHoDlwP/kDTMzB6rg98ZMWHCzbRu3Yp//Ws6JSXvcs45E+jXrzdf/GJPt+f2GtzeqtXr+e2Uh/jW4QfQtmiXyvLdd+vAvbdcwnmXTuWJf7zCuDEnMf3G0Rx+3JUAXH7xifTpvSdfOnQ0XT7XkZn3XkHJWx8we85/cur+MiGbkJ6kGicCN7OMxpMkfQ7YH7gnKto/2r+RpvrrsfP1gveo6hEzezG+AeuBimS5maUrS10znCBCF5jZvxMm1kf1XjCzu4HvEURrdEPfW1lZObNmzeXCC0+jXbu2HHjgfgwdejCPPPKM23N7jWLvkZkv8dis+axbX1qlfNjRB1GyZAUPPjGPioptXD15BgP696TvF7oCMOKEw5h0w4Ns2LiZxW+v5PZ7nub0k76Zc/eXCY2RTPEZm5KAqQS9uC4q3j3ar0tzybrY+XrBhSqHiJa4vwS4zcym1lTfzN4FPgYa7p9wEUuXfkBBQQG9e3erLOvXrzdvv73c7bm9JrXXv293XitZVnlctqWCd5etpn/f7nTq2I6ue3bm9ZIddl8vWca+fbvXyWZjP88U2QiVmXWqacvQ7LXAccBIMytJnKsuX75e8+g99JcjSOoP3Aa8BIzK8JqOQGdC7Dh5LoNu/+KM/SsrK6dDh+IqZR06tGPz5i0Zt5ENbs/tZUq74iLWrPukStknm8po366I9sVtANi4qazy3MZNW+jQrm2dbDb280zR2DNTSJoIjAEuNLNpsVNro326nlNn0ve0ao33qHIASbsCDwJbgBPMrKL6qmolqbWkLxLGsgrZETduMIqLiygtLatSVlpaRrs6/sG7PbdXVzaXldOhfUI02reldHM5pWXhT2nX9jvs7tq+LZvqKCiN/TxTNGboT9IEYCzwCzObkjj9ZrRPNxY1gPRjV7XGhSo3mAb0IaR5vr+Tet8HtgFbCWt0DAEmADcmK9Zjtx+AXr26sX37pyxdurKybNGi9+jTp0c2zbg9t1fvthYuWcEB/Xe0W9y2Dfv07MLCJSvYsHEzq1avZ0D/HdHxAfv2oGTJijrZbOznmaJAlvFWFySNA64ArjCza5PnzWwFMB8YEf9mStIRhESwB+vkQAIXqiZG0mWEBIpLzaymkdg5hJTQA4F+QGczG2dmDT6vSnFxEUceOYgpU+6irKycl19eyFNPzWPYsCFuz+01ir3CwgLatGlNYWFBlZ8fnfkS/fvuzXFHH0ybNq0Ze9HxvFGynCXvBBG5a8ZzXHbBcDp1bEffL3TlrFOHMv2vz+Xc/WVCY3zwK2kMMB54nJBVfEhs+3Ks6qXA/wD3RLNZnA5MB+YBf629B2l8aoR3XN4i6WFgoJn1qub8EcDfgQfM7JQa2loKLDCz4+rPwyVZ/fI3bNjE2LE3MHfuAjp16sCYMWdw7LGD688dt+f2ItItnPiri0/g8otPrFJ29eQHmDh5BkO+sT+TJ5xJj+6f46VX3+YnY25i+Yo1QPiOasrEHzH8mIPZUr6N39/0KFNurfrlR20WTqz78+ybtZws3fRYxn+zvTocWyu5kvQscHg1p5fF32eSvgP8H0GwNgEPE0KF62tju1qfXKgajp0JlaS9gVeAD4FDzGxzDW0tpYmFynEaizxZ4TdrIVlWmrlQ9WxfO6HKRTzrrwmQEdGkaAAAB/FJREFUtAswA9gNuBgYoPR99U/MbGFj+uY4Tu7SYpQnS1yomoauhLEmCDHd6pgDDG5wbxzHaRb4wolOvVNdmM7MlpLlP46qG+dyHCd/cKFyHMdxcpo81SkXKsdxnOZCvq7w60LlOI7TTPAeleM4jpPT1HHl3maLC5XjOE4zobCpHWgiXKgcx3GaCd6jcvKQvk3tgOOkZcvyBl8QoJmSn0rlQuU4jtNMkAuV4ziOk8vEVtTIK1yoHMdxmg3eo3Icx3FyGOXpEoIuVI7jOM0ED/05juM4OY6H/hzHcZwcJl+z/urcj5T0ZUl3SFomqULSJ5JelPRLSR3rw8mGRIFJkt6XtD1aSbepfTpTkknq1dS+OI6TOyiL/1oSdepRSRoJ/BF4Hfg1sBgoAgYBo4HPE1awzWWGAZcB44HZwJYm9SbwBOEZrmpqRxzHyR2k/JxEqdZCJelQ4EbgceBEM9sWOz1T0rXAoXX0rzHYL9r/wczWNaknEWb2MfBxU/vhOE6u0bJ6SplSl9DfL4FPgXMTIgWAmZWa2SwASaMkPS/pY0mlkhZIOleJFBZJQyU9K2mdpC2Slkq6N3Z+cBQSO1XSjZLWSCqT9KSkPkkfJI2U9EYUklwj6Z54OC0K810dHa6N2h4fnav8OdHmUknTYsepMN03Jd0S2Vkr6UFJXRPXtpH0O0kfRn7/U9LBO2mzV8Luw5K+K+nV6PksknT2Z34zjuO0SDz0lwUK/c+hwHwzW5nBJfsA04GlwHbgYOA6oCtwRdRmb+BJ4GHgGkIIbm/gu2na+y3wAvBDoAsh7PispAFmtj5qbzwwDrgd+BnQDbgKmCtpoJl9BAwHRgE/Ao4ESoEVmT+JKtwOPAacCvSI7uFOwnOK1zk58v9ZoD/wALBrhja+ErU7CfgI+Alwm6S3zey5eEVJG2pqzCw/F2FznOaLp6dnwx5AMUF4asTMxqR+jnpRcwgz1l8i6UoLb8yvAm2Ac8xsY+zyO9I0+b6ZnRxrcyHwIjASmCRpN+BSYIaZnR2r9zLwCmHc7Jdm9qqklDDNN7MaX+474XEzuyhmazfgGkl7mdkqSfsB/wtMNLPLo2qzJa0GMp2Bc3dgkJl9ENl4jiCEI4Dndnah4zjNn5bWU8qURpFnSV+JQmErgW3RdhWwGyHhAuBVYCvwgKSTk2GzBPfFD8xsHkE0D4+KBhGSOu5M1FtASPwYUqcbSs+jiePXon3PaJ/y7f5EvQeA/2Zo45WUSAGYWTmwJGaD2LlONW0Z2nQcJ0eQlPHWkqitUK0BykjzgkwSjbM8Twjj/Qw4DDgImBhVaQtgZu8A3wIqCCGyD6LxpR+mafbDNGWrCT0OgM47qbcqVq8+WZs4roj2RQmfVscrmdl/01ybqY2UnaI05Y7jtDBEYcZbS6JWQmVm24GngQMl7VVD9WGEMOEJZna3mc01s/nAZwZIzOx5M/seYczmEGAR8P8kHZGoumcaO13Y8SJfu5N6e5GZMFQQQpFJaityKZtd4oWSWtWhTcdx8gplsbUc6hL6m0QYZ/pT9LKtgqR2ko5khyBtjZ0rAk6vrmEz2x6F81JjWwckqvwgYetrQC/C2BeERIty4LREvQOAAQSRrYmlSbuShgLtM7g2HakxpJMT5SfiM4Q4jpMB+Rr6q/UL0szmSjof+AMwT9IthA9+2wBfIyQ23A9MJYxJ3S3pGqADQYCqpLRHHw8PIWT+LQfaAecSBC4pLN0l3QdMI/Safg18ANwc+bZB0kTgKkl/ifzoShgXWwlMzuAW7wQmSJpAEMD+wPnAxp1eVQ1m9qake4DLFP4veobwDdeYqM1Pa9Ou4zj5RMsSoEyp07/kzewmSfMIWXRXEMJa5cBC4AbgJjPbJOkkgkg8RBijuZUwVnRrrLkFwFFRvc8Dm6Ky75jZfxKmLwO+QUh5LyYIyQWp1PTIt6slfQxcQMiKKwVmAZdGH9TWREpUzySMrf2b0Bt6JINrq+Mswn3/GLiEcH8nAzOppQA6jpM/5OsyH2pO39JIGkzoiQw3s4eb2J16QWGGj38Bp5nZXY1svvn88h2n5ZF192jbpwsy/pttXTCwxXS/fGykEZF0FHAg4VuuCuB/CL3Dt4AHm9A1x3GaAQV5uh5Vft5107EJOAa4C/g78AvCmNzhZpYLk+E6jpPTFGSx1R5J7SVNkbQqmq5tvqTv16nRuvjTnEJ/Tr3jv3zHaTqyDs19am9m/DdboP1qHfqTNJswZdsvgPcIY/UjgGPN7Mnatltrf1yo8hr/5TtO01ELoVqYhVD1r5VQSTqGsNTQ8Wb2UFQmwsQNu5vZvrVpty546M9xHKeZ0EjfUQ0nZCFXZjhH87H+P6CfpP51u4vs8WQKx3GcZkI2UyNluIJCujk/9wcWmlny287X4uczdqQecKHKb1pM+qrj5Ad9s/mbre1qELsTJrtOsi52vlFxoXIcx2mB1HGFhJ2NhTX62LaPUTmO4zhx1pK+15RaAWJdmnMNiguV4ziOE+dNYN9okds4A6L9G43sjwuV4ziOU4WHgE7AsYnyHwKLzaxREynAx6gcx3GcqjxJmFP1Nkm7Ez74PYMwEfiwpnDIP/h1HMdxqiBpV8LySScSelcLgQlNNRm4C5XjOI6T0/gYleM4jpPTuFA5juM4OY0LleM4jpPTuFA5juM4OY0LleM4jpPTuFA5juM4OY0LleM4jpPTuFA5juM4Oc3/B0KBlWIOEGtWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_context(\"talk\", rc={\"font\":\"Helvetica\", \"font.size\":12})\n",
    "label = [antibiotics[i] for i in ab_order]\n",
    "cm = confusion_matrix(y_2018.flatten(), preds2018.flatten(), labels=ab_order)\n",
    "plt.figure(figsize=(5, 4))\n",
    "cm = 100 * cm / cm.sum(axis=1)[:,np.newaxis]\n",
    "ax = sns.heatmap(cm, annot=True, cmap='YlGnBu', fmt='0.0f',\n",
    "                 xticklabels=label, yticklabels=label)\n",
    "ax.xaxis.tick_top()\n",
    "plt.xticks(rotation=90) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fine-tuning on new clinical datasets\n",
    "In addition to fine-tuning on a first set of clinical data, we can continue to fine-tune the CNN on new clinical datasets. We use a dataset collected from new patients in 2019. As above, we show just one fold for clarity and train for only one epoch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2500, 1000) (2500,)\n"
     ]
    }
   ],
   "source": [
    "# Load 2019 clinical data\n",
    "X_fn = './data/X_2019clinical.npy'\n",
    "y_fn = './data/y_2019clinical.npy'\n",
    "X = np.load(X_fn)\n",
    "y = np.load(y_fn)\n",
    "print(X.shape, y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train/val/test split\n",
    "We split our patients into training, validation, and test sets for each fold, as above.\n",
    "\n",
    "Our 2019 clinical dataset contains 100 spectra from each of 25 patients (2500 spectra total), with 5 patients coming from each of the 5 empiric treatment groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each treatment group of 5 patients, patients are assigned to [train, train, train, val, test]\n",
    "patient_idxs = []\n",
    "for i in range(5):\n",
    "    x = list(range(5))\n",
    "    np.random.shuffle(x)\n",
    "    patient_idxs.append(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Group 0 patients\n",
      " Tr: [2, 3, 0]\n",
      " Val: 1\n",
      " Te : 4\n",
      "Group 1 patients\n",
      " Tr: [3, 2, 0]\n",
      " Val: 1\n",
      " Te : 4\n",
      "Group 2 patients\n",
      " Tr: [3, 4, 1]\n",
      " Val: 2\n",
      " Te : 0\n",
      "Group 3 patients\n",
      " Tr: [4, 1, 0]\n",
      " Val: 3\n",
      " Te : 2\n",
      "Group 4 patients\n",
      " Tr: [2, 3, 1]\n",
      " Val: 0\n",
      " Te : 4\n"
     ]
    }
   ],
   "source": [
    "# Sample train/val/te spectra\n",
    "idx_tr, idx_val, idx_te = [], [], []\n",
    "for group_idx, patient_list in enumerate(patient_idxs):\n",
    "    print('Group {} patients'.format(group_idx))\n",
    "    print(' Tr: {}'.format(patient_list[:3]))\n",
    "    print(' Val: {}'.format(patient_list[3]))\n",
    "    print(' Te : {}'.format(patient_list[4]))\n",
    "    for j, patient in enumerate(patient_list):\n",
    "        start_idx = 5 * 100 * group_idx + 100 * patient\n",
    "        idx_range = list(range(start_idx, start_idx+100))\n",
    "        # Randomly sample 10 spectra for each patient\n",
    "        np.random.shuffle(idx_range)\n",
    "        idx_sample = idx_range[:10]\n",
    "        if j < 3:\n",
    "            idx_tr.extend(idx_sample)\n",
    "        elif j ==3:\n",
    "            idx_val.extend(idx_sample)\n",
    "        else:\n",
    "            idx_te.extend(idx_sample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Continue fine-tuning the CNN\n",
    "We use the new 2019 data to continue fine-tining the CNN that has already been fine-tuned on the 2018 data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting fine-tuning!\n",
      " Epoch 1: 0.02s\n",
      "  Train acc: 93.33\n",
      "  Val acc: 84.00\n",
      "  Test acc: 84.00\n",
      "Finished: 11.85s\n"
     ]
    }
   ],
   "source": [
    "# Fine-tune CNN\n",
    "epochs = 1 # Change this number to ~30 for full training\n",
    "batch_size = 10\n",
    "t0 = time()\n",
    "# Set up Adam optimizer\n",
    "optimizer = optim.Adam(cnn.parameters(), lr=1e-3, betas=(0.5, 0.999))\n",
    "# Set up dataloaders\n",
    "dl_tr = spectral_dataloader(X, y, idxs=idx_tr,\n",
    "    batch_size=batch_size, shuffle=True)\n",
    "dl_val = spectral_dataloader(X, y, idxs=idx_val,\n",
    "    batch_size=batch_size, shuffle=False)\n",
    "dl_te = spectral_dataloader(X, y, idxs=idx_te,\n",
    "    batch_size=batch_size, shuffle=False)\n",
    "# Fine-tune CNN for first fold\n",
    "best_val = 0\n",
    "no_improvement = 0\n",
    "max_no_improvement = 5\n",
    "print('Starting fine-tuning!')\n",
    "for epoch in range(epochs):\n",
    "    print(' Epoch {}: {:0.2f}s'.format(epoch+1, time()-t0))\n",
    "    # Train\n",
    "    acc_tr, loss_tr = run_epoch(epoch, cnn, dl_tr, cuda,\n",
    "        training=True, optimizer=optimizer)\n",
    "    print('  Train acc: {:0.2f}'.format(acc_tr))\n",
    "    # Val\n",
    "    acc_val, loss_val = run_epoch(epoch, cnn, dl_val, cuda,\n",
    "        training=False, optimizer=optimizer)\n",
    "    print('  Val acc: {:0.2f}'.format(acc_val))\n",
    "    # Test\n",
    "    acc_te, loss_te = run_epoch(epoch, cnn, dl_te, cuda,\n",
    "        training=False, optimizer=optimizer)\n",
    "    print('  Test acc: {:0.2f}'.format(acc_te))\n",
    "    # Check performance for early stopping\n",
    "    if acc_val > best_val or epoch == 0:\n",
    "        best_val = acc_val\n",
    "        no_improvement = 0\n",
    "    else:\n",
    "        no_improvement += 1\n",
    "    if no_improvement >= max_no_improvement:\n",
    "        print('Finished after {} epochs!'.format(epoch+1))\n",
    "        break\n",
    "print('Finished: {:0.2f}s'.format(time()-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make predictions on treatment groups for test patients\n",
    "As above, we make predictions for one group of 5 test patients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting test indices\n",
    "idx_te = []\n",
    "for group_idx, patient_list in enumerate(patient_idxs):\n",
    "    start_idx = 5 * 100 * group_idx + 100 * patient_list[-1]\n",
    "    idx_te += list(range(start_idx, start_idx+100))\n",
    "dl_te = spectral_dataloader(X, y, idxs=idx_te,\n",
    "    batch_size=batch_size, shuffle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished: 11.86s\n"
     ]
    }
   ],
   "source": [
    "t0 = time()\n",
    "# Make predictions on the 5 test patients\n",
    "y_hat = get_predictions(cnn, dl_te, cuda)\n",
    "# Subset y according to which patients are chosen for the test set to generate y_te\n",
    "y_te = y[idx_te]\n",
    "print('Finished: {:0.2f}s'.format(time()-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate predictions of treatment groups\n",
    "We use 10 spectra per patient to make our treatment group predictions for 5 patients. To compute identification accuracies across all patients, repeat the fine-tuning process in a cross-validation loop so that each patient appears in the test set once."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ran 10000 trials: 3.73s\n"
     ]
    }
   ],
   "source": [
    "t0 = time()\n",
    "\n",
    "# Experiment parameters\n",
    "trials = 10000\n",
    "n_spectra = 10\n",
    "\n",
    "preds2019 = np.zeros((5,trials))\n",
    "y_2019 = np.zeros((5, trials))\n",
    "\n",
    "for i in range(5):\n",
    "    y_2019[i] = y_te[100*i]\n",
    "    for k  in range(trials):\n",
    "        y_patient = y_hat[100*i:100*i+100]\n",
    "        idxs = np.random.permutation(np.arange(len(y_patient)))[:n_spectra]\n",
    "        preds2019[i,k] = stats.mode(y_patient[idxs])[0][0]\n",
    "\n",
    "print('Ran {} trials: {:0.2f}s'.format(trials, time()-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/envs/inkjet/lib/python3.7/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in true_divide\n",
      "  \"\"\"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaoAAAFoCAYAAAAGr79JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXeclNX1h58vqDSxK7aIREWs0VhiDWDUqPEXbGiswRaxxxo0Kvbegg1Q7EbsBQuiIhpRUSQKKmIDlCIRkbqACOf3x73DDrOzO7s778zuzJ6Hz/uZmffeuee+M8t75p57iswMx3Ecx2msNGvoCTiO4zhOTbiichzHcRo1rqgcx3GcRo0rKsdxHKdR44rKcRzHadS4onIcx3EaNa6oHMdxnEaNKyrHcRynUeOKynEcx2nULNfQE3Acp/yQ1A5oD7TMbDOzt4o/I6eUcUXlOE5iSFoHeAjomq0ZMKB5USfllDyuqBzHSZLbgW2BfwCjgYUNOx2nHJAnpXUcJykkTQfON7N7G3ouTvngzhSO4ySJAd829CSc8sJNf47jJMlTwH7Aaw09kUIhaTfgSLI7i5iZ/aH4sypvXFE5jpMkjwD3SgIYBPyY2cHMRhd7Ukkh6VhgADAD+IKqe3Aq+qSaAL5H5eRE4a7zW6p3N/530SflNEokLUl7mXlzEWHFUbJef5LGEpxEjjEzdxQpEr6icmpE0sbAc0Ansv9aNMAVlZPi2IaeQIFpD5zuSqq4uKJycnEH0JZgk3d344SIq9QdqX6V+mDRJ5UAZvZAQ8+hwHwOrNXQk2hquOnPqRFJM4GeZjawoedSLkjqCDwPbEI1q9RCmMcktQQ2ILtiLNl9o2IiaT/gOmAvM/u+oefTVPAVlZOLCuCnhp5EmXEH0Ao4jCKsUiWtDvQHulFVMeadLULSvcAVZjY+Pq8JM7Pj6yurEXACsBLwpaRRVHUWMTM7uPjTKm98ReXUiKTrgDXN7Lgiy90CuAToDKwO/M7MRkm6DHjLzF4v5nySRNIs4G9m9liR5D0N7AncA4wFfs7sk4/JTtJ44AAz+1jSBKo6UWSIsl/XV1ZDU+7X11jxFZWTixHADZIep3p345eSFChpe+BNggvwq8ARac0tgZ5AySoqYB4wq4jy9gDOMrMBhRjczDqkPd+wEDIaC+V+fY0VX1E5NZLhbpyOUSB3Y0mvAcsDfwR+IawAto8rqu7A9ek3x1JD0jXAOmbWo0jyvgOON7MhxZDnOEnjKyonF3s1gMydgCPMbIGkTCU4DWjXAHNKklHAtdEkV90q9fkE5d1H8NosuKKSdCiwoZldn6XtfOAbM3uy0PNIEkkbAFPNbFF8XiNm5imkEsYVlVMjDbQXtARYXE3b6sD8Is6lEKT2pjoAB2RpT7oUxmVAP0mvAi8STKrLCkzOHf48Kq8vk0WxvaQUFTAe2Bl4H5hAzXtU4GVMEscVldMYGQX8hXBTzeRAwr5ZKZOtVlMh2QrYG1gfyJaHzoCkFFVH4KNq2kYDvROSU0yOA75Oe+77JUXGFZVTBUlDCNH34+LzmjAz+2PCU7geGCRpOULuOAO2l3Q0cDjhpluymNmbRRbZl7BK7UE1Xn8J0hxYsZq2toS9x5Ii3SPSzO5vwKk0WVxROdloRWW8TWuK/AvSzF6SdCJwI3BonEtfYDZwopm9Ucz5lAFbA0eZ2dNFkDUGOAR4NkvbocCnRZiDU2a4onKqYGa7pz3frYHmcK+kgYS9gXbAdGC4mc1riPnki6ShwClm9nl8XhNJl4qYTNgfKgZ3AA9KqgD6AZMIJseehADnks4FmCOgeQkh7OAD4BnPB5gcrqicRouZVVDa8VLppGeEaEbNq9SkS0XcDJwl6RUzK6TZDzN7WNLmwPlAegaKJcA1pZrDMI2uwMrAKoTQiR8JDj7LATMJ391ZwDhJXcxsWkNNtJzwOCqnRiSdDLQ3s15Z2q4BJphZvwLILcukrQ2BpJsJTijNgaFU9fozMzsnYZkdCKENawA/AK+a2YQkZTQEknYgeC2eDzxpZotjCEV3Qg7A7gSl9TTwgpmd0GCTLSNcUTk1Imk00N/Mbs/SdjJwkpltk7DMBknaWq7UELSdwj/PWiLpLeAJM7stS9uZQHcz203SWcC5ZrZe0SdZhrjpz8nFr4HPqmn7PLYnTVGTthYbSacSVqnnZ2m7nrBKvTMpeWbWLKmx6oKkNQnf4zKUeEDsdoS4tGx8Alwdn39EWE06CeCKysnFYmC1atpWpzClt3ckJG0ttcDQ2tKT4MWYjW+Ak4HEFFUxiSbbS4FTgVWr6VbKq7fZhETJ2fZOu8Z2gBbAnGJNqtxpkF9aTkkxihB/k43jgP8WQGaxk7YWmw6EeKZsjKMwq1Qk7SvpWkn3pFIBSdpJUpIpqU4CzgZuIfyIuToe3wFfEcpklDL/BnpJukzSlpLWjI+XE/atHo79fkuwODgJ4HtUTo1I2p+wXzSYqu7GexPKOwxKWGZRk7YWm1iM8kQzeyJLW3dggJmtlKC8FQk5BTsTVsjNgB1ikt+BwBQzOzshWf8FBhJi4BZRmUy4JfAG8LyZXZOErIZA0grA/YTMKek3TwGPAj1iTsC9gTlm9m7xZ1l+uKJyciLpFIJHU+vUKWAucL6ZVWfCykded+Ba4GOKk7S1qMQ4qoVmtm+WtleAlmbWOUF5txFKpfQgmKzmUqlAjgX+bma/SUjWHODPwDCCotrNzN6Lbd2B68qhXlN0we9MMIv/SKiRVt1erpMnvkfl5MTM7pT0MLAble7Gw81sds3vrDfFTtpabG4EXojpqTJXqXsQbvRJcjBwkZkNypKN/jtCefqkmA80NzOTNA3YEHgvts2m9DPfAxCVkiumIuGKyqkRSa3NrCIqpUQLJNZAsZO2FpWYIqoncAOVSWJFuJH3NLNsyXjzYTWq3y8RYeM/KT4nKCcIyYPPkvQfwurqfODLBGU1GHFfr7oYv7eKP6PyxhWVk4vJkh4A7jKzccUQ2ABJW4uOmfWX9G9gFypXqe8UKEXUtwS36mw5ErcjWeUxENg0Pr8SeCvKh5DJ4ZAEZRUdSesAD1H5Yyrl9bq0kCilvdpvlPgelVMjkvoARwMrEbIa3A4MMrNcQaRJyF6JUERxdcI+wHsFNDcWDUkrmtncIsq7HDiNkHn+NcLqZjvCTfUV4KZshQ4Tkt0e6BZlvWpmJe0JJ+kpoAtwDdXE+DWFH1rFxhWVkxNJrQnKqifwG8K+Rl/gHjP7oUAyLwQuZNlM7hXA5YW6qRYLSbMJ9Z/uLMYGvKQWBLNtF4LCXwP4HliLsMra18x+KfQ8ygFJ0wlORDUlp3USxuOonJzEPap+ZrYt8HtgOKEA3neSHpK0U5LyYomPKwml0w8mZFA/iPDr/5roqVbK3ENY3YyR9Iakg7M4OSRGzOK9F+HHxiuEVdUwQibzfZJUUpJOr6FtBUl3JSWrgTAqTZlOkfAVlVNnJO0FXATsTmVczn+AY81sfALjjwE+MLPjsrTdC2yXlDt1QxHjio4krFK3A6YQPADvNrPvG3Ju+RDzCj4DHGdms9LObwo8DmxkZtUVVmz0SOoLVCQVd+bUDldUTq2QtDLhF3hPQrLYUcBthEzS+xPinian17LKQ9YCoJuZvZKl7Y/As2ZWJYdcAnIbxJMrrkhPIWTebkbIvH2Hmb1dKJmFQtIBwABCyYu/mNkHko4h7G1OAg41s08aco75IGl34F5CfF91MX6jiz2vcscVlVMjkn5LuIn+hVBG/CngtsyI+6hAnjOzKjf5esicDpyTXgI8re1Y4HozWzNfOWljZnpyLdNMEbKLS9qXsCe3K8HZYTngXeCvZvZ1nmMvD/yDsIJrT1V3dDOzxDyAowPFY8A2wNuEz/Vh4ORYY6xkychEn3nzLMrfSlPE3dOdXIwEpgE3AX3NbGo1/cYTKpsmwevAlZJGmtnS0uWStiRkrn41ITkpbge2JdzMi5atXdJqhHyJPQnBze8DRxF+DPwJuB54gBBonQ/XEIr5vUwwyxX0+sxsYiwB8w4hgPl94Pgycdgo9f3RksRXVE6NSDoaGGhmxSplniq69w6wJkFxfA+sDWxNUJo7J1kqotieXJJ2JKxSDyWY+h4nrFI/yOi3NyE3Xl6rVEnfAvea2aX5jFMHeccBfYCvCSuriwnf42HlUDzRKT7u9efUiJk9VEwlFWWOJ5iNbiKYwTYBfiasMLYpQD2jYntyvUfwwrsG2MDMjslUUpEJBPNfvqwKFCW2JwaH30PIMv47M7uaENS8GvBfSSUd8Os0DL6icnIiaQ1C0GY2RwMzs38Uf1bJUWxPLkmHE6rEFsUUJull4BUzu7UIsmYTqj4/mnG+LXA3oQJuye7hRK/TmjAzO74ok2lCuKJyakTSH4BngTbVdEl881jSdmb2YQ3th5nZY9W110NeWXtySepI+A57AS8VUkFK2sTMqk3JJOkkM+tXKPmFRtIEqjpRrA6sSPB0nFkO2eEbG66onBqR9DEhI/YphGJ/P2f2MbPFCctcCJxrZrdlnG9BcIk/Pknl2BCeXJLWAg6k+lXqOQnK+ong6deCEPeWmb7JzGz1pOQ1RSTtQajK3N3MxjT0fMoN9/pzcrExcJCZjSqizLuAf0nqSgginiWpEzFgFPhbwvKK6skVA6afobK+VyYGJKaogOeoqoALRjQVn0VwS0/laXwduNXMqqxWywEzGyrpduBfBE9HJ0F8ReXUiKSRwA1JmtpqKbcbwRw3Kz72Ijg8lHTAKCzNvDGH4JY+ttjOKoUkxlANB9YBxgBT4/OtCAG/u5rZpIabYeGIq6rnSznzRmPFFZVTI5L2IdRN2rfYN5hYRfVDYAVCjFZXM5tfzDkUAknzgAPNbEhDzyVpYumSzsB+ZvZx2vmtCYlxh5nZUQ01v0Ii6TJCgPaGDT2XcsNNf06NmNlgSbsCX0r6FJhRtYv9MWm5krYAniDUMPoPsCchIe15SaxAovfWFWY2vgE8uT4luGsXjfh5XkJQIqsTXMdHxZvrW2b2ekKi9iLEpH2cftLMRku6BLguITkNQryGTFoQYvz2JfyocxLGFZVTI5JOA/4JLADakmw12OpkHkvIFvEVsL2ZjZN0CiGualdJhyaQ/LYrYT8Bwp5CTaaFpM0O/wRukfSumU1MeOwqSNqeEEc1g5DV44i05pYEE2RSiqotMLmatskE77hS5tIs5xYSYt56E2L9nIRx059TI5K+I2SJOM4KU302m8wlhJibM81sQdr5bQkOFWuY2arFmEuhkHQlcC6hdHu2Veofqr6r3rJeI+Rp/CNhhfoz4QfAKEndCbkTOyQkawyhwOWJWdoGADuY2dZJyHKaDr6icnKxGtC/WEoqcoSZDcw8aWb/jUly+xdxLokj6UxCAtp5BAVSaKW7E+EzXZCl7tU0oF2Csu4A7ozVmR+hMv3VUYTaYj0TlFUUJH1D2FP8OJr+7jGzKQ09r6aEKyonFx8QXNSTMg3lJJuSSmubQyg6mBeSNqjjnJJMsXQeMJAQD1YM55AlhPipbKxOiJNLBDPrK2lt4HwglS5JBNPxJWZ2d1KyisgGVMa69QYGE+qHOUXCFZWTizOARyWNB1618rEVT6Bue09JBvyuCgwoogfjKEKZlheztB0IjEhSmJldKulWQmXm1QhxVO9aWiHFEmMKsJekcQSl2yauGLNiZrOLNrMmgisqJxevEwJTXwYWxVxu6ZiZJWk6ApY6VJwGdCJ7IcN8FcdxFDEINoMRhES7xVqlXg8MkrQcwRxnwPYxM/7hwN5JCzSzmYS/mXLgEeByQokZA17L0b9kcxk2VlxRObl4hSLf0OMGfz9CLaZtCQG/LYA/E/Y8/p2vDDO7P98x8uAMYGAsv/FK0imoMjGzlySdCNxIKC0ioC8wGzjRzN7IZ3xJv6/jfApWLbkQmNkFkt4DtgSuIFQwLsug5caKe/05jQ5JIwirjYsJZT5SHmprAW8Bt5R4YtMGyb0nqTWh5MZawHRgeBJOMtFLM3UjEdX/sCn5CrjRBH5AZpyYU1h8ReU0RjYFLiI4ARjBMw4z+1906/4nYcVVb2oR5JtO0gG/Rc29l8JCGfhcZqv60LUAYzZKknLjd+qGKyonJ5I2JqxuMpOMXmFm3xRA5BLgFzMzST8A66W1TSVkHM+XXEG+6SSqVMysR5Lj1QZJKxNcw1OZKX4EhgH9EnBy6Ex02Y7elFPLKX9hJk0x6W5D46Y/p0Zi6p13CPn2XqcyyegeBLfm3cxsbMIy3yeUZn9I0mCCovi/+PhvQpXfTZOUWc5I2gh4A1if4O2Yip3aEPiOkEOx3j84JC0Gdjaz99Of5zvvxkhTTrrbkHgpeicX1xAcGDYys/3N7EQz259QbuN/wNUFkPki4Vc6hNxwexKyqM8gxOaUfJoaSZtKekTSFEkL4+NDkjYpgLg+BOvJTmb2azPbORb325lgVu2T5/izgDXj85r2qMqBawjX+Fsz28bM9jWzbQhOP82Baxt0dmWKr6icGpE0EzjZMkqLx7YjgDsKnc5I0u+Awwg3wBfNbGgCYy41UdUm+DfJgF9JWxF+lTcn5N5L/Srfk5DiaFcz+zRBefOAU8zsgSxtxxJWr/XOwRdXvVsD7xM8M98Efqqmu5nZwfWV1dBEU/T5ZnZflrbjgOvMbM2q73TywfeonFy0IPxizsZMipCk1sxGkHBQKjCesKJ4n9oF/ybpqXY1IUFrVzP7PnVSUjtgaGzvlqC8CoK5LxtTY3s+nAzcAmxO+Bw3JSRqzUap/zIu96S7jRJXVE4uPiUEx76Upe0EoFSLGB4HfJ32vJg30N2BnulKCsDMpkm6ghDjlCSPAT0IqX8y6UFI9FtvYib7A2Cpq/oB5bpHBXwJdAey1RI7NLY7CeOKysnFDYQUSm8THBlSSUaPJCQ7PTQJIRmJP8eTo+yGmW2Uj7x0M1gDBP+uQKjwm405sT1JhgPXSfoPQSmlvsPDCA4WvST9OdXZzJ7PQ1ZX4LM83t/YKbuku6WA71E5OZF0AnAVlRvmEExJFyR1k5d0H3B5LGR4PzlWOGZ2bBJyo+y2wIpmNjVL2zrAHDPLDMrNR94HwHdmdlCWtmeBdczsdwnKW5KjixEcBCChgFxJHYEuVLpvv2FmZbHakHQpIeluyuydSrp7lZld1VDzKmdcUTm1QlIzYAsqk4x+Zma5boAlQSyfvsjM/pql7T5gOTM7OkF5hxBWNu9TdZW6PXCwmT2boLzOuXtVYmZv5iGrGcF0eTyVyg9CbFx/MzulvmM3JiStQsjysSqln3S30eOKymlUSFqBkMz0GjMrRBaFbDInAWebWZW9mqhUbjazOpUFqYXMHgSnibXTTk8BepnZw0nKKiaSLgCuJFRPfpBKj8ajgb8Trs/LtTt1whWVUyOSegPrmdnfsrT1J5iwrkhY5izChnxeyVLrIG8h8EczG5alrQsw2MyqZHBPQK4I2eFTq9RxhSyjEvdVdqLSHPde0iUpJH0BPGdm52VpuwH4c6kHaxc4y4eTBQ/4dXJxOMFElY33SKCIYRbeJHjGFYv/EVyrs7ElVUvF54WkSySta4GxZjbczD6PKaPWiVVkE0XShYQV28sEJ4DBwBRJ5ycsqj0h4342XiGZ9FcNRszyMYYQ+NuJYNLsRAj0HS3p1w04vbLFFZWTiw2Ar6pp+4bC3Hh6Az0k9ZLUSdLKklZKPxKW9wpwUeZNJr7+J9nduvOhN8HbLhvrxvbEiCU+riS4VB9MiB87iHDd18Sg36SYAXSspq0TCSv9BqDQWT6cLLjpz6kRST8Cp9WQmeJOM1slYZkpJ41q/ziTLBUhaT3gQ2BlQnbxSQRFsichw8KOSeZvi9e3U7ZYI0l/AF4ws1YJyhsDfGBmx2VpuxfYzsx+k5CsuwnK8Egzeznt/H7Aw8ATZnZSErIagkJn+XCy43FUTi5GAGdKeizdyy96d51BMP8lzeUUMQDXzCZL2p5QFO+P8fiB4JHXOwklFYsLdkk7dYKkfTK6tSIk3/08X3kZbAKcW03bYyRrvr2A4A33QtxrTHk0rkwIHr8gQVkNQaGzfDhZcEXl5OIqwkbxR5LuIaSJWZ/gftyJAtQiMrNLkx6zFjInAUmawDLpSqVJzwhZPbLxDXB6wrLnsqx3YTrrUrVwY70xs+lR6fcgKObVgNGE1FAPmtmCpGQ1EAXN8uFkx01/Tk5i1oLbgF+lnZ4InG5mLxRA3r2EWlfjs7S1J6xyqpixGjOSWgAtCbFFMwhmxQ8zui0sxI1c0mOEVc4+6cluJW1JSI31tpkdkYCc5aKcr82sunx4JY2kwwkZ/SdSTZYP0lZVeWb5cCKuqJxaI2lzYA3gh6RrUGXIqWkPZzvg/Xz3qNKVYS2q/ZolWOE3KtupZvZzUmPmkNeBUFNsTcLqJnVz3Zpgxto5iezw0d1+IbC/mWXLhVfyNESWD8dNf04dMLPGkMOtHSFdTb50JQSlQu5qv0lX+J0IIGldYDcqY3HeNrMpScqK8sZL2gY4m2CO2yTKux64xcx+SEiOSZpIgTOIS/oTYX+0PXEv0czuKqTMNBI3dTu58RWVkxNJbYB9CDeGzMBXM7NrEpDRjcrSFj0IxRMzb6CtCEGW481s13xlNhSSmhOU5EksWz7kF0L6ob8nlZ4qOr2sD/xkZtUlwk0MSWcARwB/MLN5BRi/K8EzU4S/j9UIYTaXmdnlSctzGge+onJqRNIOBKWxRjVdjBD8mC+dgAPTxuwMLM7os5BQVuTsBOQ1JJcTajj1AwZSWRr+cOBUYDZwUUKyRHDQOABIfD8xCx2IJe8lvU7whEv/NWxmdk4e419I8IrcP64UVyakajpb0pUNkX8yKs+tgaFmNqbY8psCvqJyakTScGAl4ETC/kaVgnhmlqlQ8pVZ7R5VIYmFC7OtGjGztxKUMwV4pJo0QzcCR5jZugnK+5YQ+1NwRVWbPZx89m0kTQXOMrOBaec2AcYBmxY6Q3tMYCwzOzy+PhpIxVSlUnEl9rfiBDwzhZOLrYF/mtl7ZlZhZoszjySESBoqqVN8eSyVRQ0LTkxb9BohxdC7wBtpx7D4mCSrUH2aocGEmKMkuQ84OTo7FBQza5bjyNe5oB2hOnM6qdfFCLTdmWBhSHE+8CSwIcEF/59FmEOTw01/Ti5+pDjBt52pvEHfS7gh/FgEuQC3A9sC/6CaVWPCjCT8AMiWHf43VHVbz5fZcdxPJT1PdnPcbQnLLCSZq7ZimoXaAd/B0lplWwDHmdm3ku4k/O06CeOKysnFXYTg1EEFlvM/YCtCJgxR3JtPZ+A8MyvWTeYM4BlJs4HHzGxOLN74F+A0Kp1KkiJVVmNdwl5gJkaIk0sMSb9j2ezib5rZiISGv17ST+ni4uNNkmamnTczOzghmSnmE5x6IHhsziP88CA+b5uwPAdXVE5u5gFbSnqPYPKossoxszsTkDMI6B/3aAx4TVJ1ZkUzs9UTkLl0PCDvOKI68AahOmw/oJ+kBVTuiy0EhqVZ6ZK41g55vr/WRA/RgcB+LFs40SS9CPzFzPJJM/Qtwcy2Ycb5iVS9zkL82BkNHB33bk8kVC5OydmA8IPLSRh3pnBqpNCb42lyViB4wm0JHEcoR1FtfI8lW4q+L1BhZkXxJpR0P3W4iSZ5rYUmmr+OBy6mqkfj5cAAK+Eqv5L2IPxgWwH4GdjTzIbHtoHA8gVYxTV5XFE5NRLr79SImSXq+FBsrz9JuxP2FgbFI9uqcXQx5lJIJG1GMMetAdxjZt9L+hXwY56rnHQZ0wkVka/O0vZPgsdedaEOJYGkDYHtgP+a2Tdp508CPkrQxOlEXFE5TZ6MVWPmfwhR4qlwYg6+AcBRVO7/7WBmo6JzxWgzSyRuS9IcQnXm17O07Qk8bWb1ricWi0q+ZGYjc/TbDLjDzPaoryyn8eB7VE6jRFIrQoaK1ArgZDP7UtJBhBtrdcUc60ODmNYkbUVwbqhSe8rMHkxQ1CVAd+Asgvt7ehmRlwmm1qQCjN8kfGdVFBXwe+A/eY5/KXCBpNPMbEAN/VaK8ygYktYk+3dXzP3OJoErKicnsZbSqWS/qZqZbZqwvHYEh4NNCWVF1qPSm+pPwN5Az6TkWZYieIUkevi9QPAag7QkpmndklRUxxBSDPWJ6ZvSGU+yzhYXAU/F1E2PUZkA9y/xODC9QrOZza6HjCkEx5sdCUU9F+U/7doRY9EuJfx/WLWabiW7+m6seMCvUyOSuhACGdckuI+PB6YDvyb8/RRiH+l6QvDmtlFOuvfYGxTwl7KkTpJ2i95rheJSgqt4Z8K1HUgo+/EE8BWwfcLy1qH67+lnIMlrHUVQfBcCHxEU1UeEgokdgP8SqianjvpwNHArwevurZjct1icREjhdQvhu7s6Ht8Rvrvq6ow5eeArKicXvQmOBicDi4AL4t7G1gSz0ZMFkLkfcI6Zjc6yAkitsBIllhG/iuChBrADMErSk8CrZtYvQXH7E/Ijvhtff2dmo4Chkh4g/FpP8oY3lbA6zZZhY0uSdc0vRnXmX8zsHEkjgHuADyUdVqTURScBVwI3EipCPxP/P1xF+HyrK1Dp5IErKicXWxAKxaVuPs0heMFJupyw//FcwjJXJCikbLQkYUtALIY3AHiWoHzTldL7wKEZ5/LlV8CXZrZYUuaK5lGSNftBuK6LJP2Hyv0pk7QeYd/q0aQEWRGrM5vZ45I+AZ4mxN2db2a3FljsxoS/iSXxWCHOZYGkmwn/V5JI0uyk4aY/JxctCDFGSwiVadultX0BbFYAmV9SvXlvNyDpuli9gLvN7CCqpsD5nOSvcQaVe25TCD8GUqxM+MyT5FJgJsHs9i7hR8c9hM9xFmElWZLEGmk7EPb8bpL0SHTEKRTzgeYxyHcaywYez2bZ/x9OQviKysnFBCrNGZ8TVhepLNwHEP6zJs19wFWSJhM25AGWizWrTgPOTVjepkB1pSdmEmoeJcmHwOaEMvAvAL2jW/ciwq/xRPf9zGxmTGl0JrAvIeHvQoL56pZ860ZFl/F7zGxKfJ5jOnZFPvKyDDgHOEhSL4I5bgvCPmch+JxK5TQCOCuuVBcREtQWNHt7k8XM/PCj2gO4Gegfn3e/uPmqAAAgAElEQVQnmDs+J6SSWQxcVACZIpROWEK4AaQeFwP3FkDeD8Ch8XnzKO+38fVRwKSE5e0JnBSftwPGUGlK+grYrKG/9zpezxJgx7TnNR2Lk5JVTfsfCGmMfs5XVjXjnwLcEJ//Fpgb/y4XE5T//zX091GOhwf8OlWQtLxFl19JLYFWZvZTfH04IR2OAS+aWf8CzmMXwgqgHcHT8GUzyzcOJ5ucJwh52joTFOIiQuaBjwlxQV+bWY+k5abJF8Gj0oCxZvZLAmPWJcGumdnx+cosBrXJWhKzbTwFbGcFDtSW1J6QRNgITjef53iLUw9cUTlViGlwHgHuM7OPiiRzJYKH4ctmlq38RSqzwb7AJZZgmfOYxWAEYWX1NMH9eACwDbAJ4Yb3TfUj1EnWCsBhhODXdeLpyQSF+IQlFBMUb+gLqF3JEjOzeps3Y1zYpRThu4uKYUquz0nS8sC6ZjaxvrLSxmqwv00n0tBLOj8a30HIKpAyZ4wiBjcWWGYvwg27dQ19WgOTCK7rScvfChhCMBktAX6Jn8NWCcrYhpDlezFZTGKEkvFbJiRrJkFRPUa4gaocvrv4OVVr+iv16/Mj++Fef04VzOwPhODMywjeabcBUyQ9KmmvAontRtgLqzY5amzrBySWnVrSzpL+AqxFiG9qC6wPtDWzP5jZmITkrERIeNsWOI/gwNGKcIPrRCjauArwgqQkKtW2I6SgWiXK/U7S1ZI6JjB2JsX87gpepTgLDfK36aTR0JrSj8Z/EPZu7iO43y4mrAouAzokKGMWsHct+u0FzExA3koEV+3FacdECuTIAPydUNtrixr6bEnYnD8jYdm/IqQ2+jJe5zuErA4rldp3Rw5nigJ9d0X92/Sj6uErKicnZvamhZpIaxNqDU0g3vgkZUs+Wh9aEExVuVhIZZHBfLgA2BF4iODyfgvBDf2eBMbOxn7A/Wb2aXUdzOwTQrDvn5IUbGbfmdmVZrYJ0AUYB9xJcoHFxf7uir2xXuzrczLwOCqn1lgwb9wv6RlCLrdzCTe+JPieYA7LlQZnU5KJ3eoG3GRm56dOxJQ8AyWtbmZValLlyRZAbTwkh5J8KXoAJK1PcOLYneCGPymhoYv93T0rqbZOIjnrqdWCYl+fk4ErKqfWxP2pHoRA35aEmJ/7Ehr+LaCnpHvNLGsJ+pj37yRgWALyNiSkS0rnFcIeyHpkKZ6YJ6sScu7lYioJBhhLakFIenssIcZoGmEVeZ+ZjUtITLG/uy+pf0Lb+lDs63MyaWjbox+N+yDkNruSSm+1OQTltHvCcrYlxC89DbTL0t6OEBvzM7BNAvKq7HWQEeyb8PXVam8F+B0JBKoS0grdSUjXtAB4nOD916wA11a07662n2OpXp8f2Q9fUTlViF5nhxJ+he9CWGW8S8iM/ZiZzU1appn9V9LpwO3AnySNJOyFAbQnlL5YjlB/KKnYro0lpe89pIJDN5G0TNCtJVOKPlNeNjZJQA6EuLDZBAX1b4LCAtgyxBcvSz7X10DfXdEo9+srBTzg16mCpLkE1+mUmWiAmX1RJNm7EBwduhJctwEqCHs315jZu9W9t45ylpB9U14Z5xMpRV+DvKzyE5JHLWQmIi/KLPh3V5vMFIWiWH+bTlV8ReVk41VCFvGXrBqbfKEws3eA/4sVYteIp6dbyN6eJMUuP1/u8or53TUI5X59jRlfUTmO4ziNGo+jcuqEpJmSZro8l+fyii+vqeKKynEcx2nUuKJyHMdxGjWuqBzHcZxGjSsqx3EcZxkkrS/pX5LeljRXkknqUk3fIyR9LGmBpEmSro0FVzP7tZP0gKTpkuZJ+k90+c+JKyrHcRwnk40JlbznEuqyZUXSUYQiq8MJmU+uJtSvuz+jX8s4TmfgdEJarznA65K2zTUZj6NyHMdxMnnLzNYCkHQA8OfMDjG/4Q3A82Z2Sjz9hqRFQH9Jt5jZiHj+OEJi5u3MbFR8/5vAWIJy27emyfiKynEcx1mGWgYx70Qo/fNAxvlHCLkR04tIHgiMSSmpKGMh8Ciwl6S2NQnyFZXjOE4ZUpv4LjNbJQ8RW8bHTzLGrJD0dVp7qu8bWcYYTcixuRlQbVosV1RNmFYbHF7ntCQtf/WXer93/reX1fUtVFai+KIoKVRcnssrnryOVbMD56CO/+9m1XX8OrJ6fJyRpW1GWnuqb3X9yOhbBVdUjuM4ZUieq6U6iarl+ZqUbI0K2BWV4zhOiRDy4TYaUsVFV6dqodHVgPEZfbOtmlJFQrOttpbSqK7acRzHqZ5mWq7WRxH4ND6m70UhqTWwEcvuXX2a2S+yFaEg6+c1CXJF5TiOUyJIzWp9FIH3gO+BozPOHw4sT6iInOIZYCtJ26ROSFoh9n3NzGbXJMhNf47jOCVCturMBZR1SHy6Q3zsLGkNYJ6ZvWxmv0jqBdwv6XbgSYL33nXAk2b2XtpwAwiBwE9LuoBg6jsTWJdQTbxGXFE5juOUDEU1gj2R8frS+DgR2BDAzB6QtBj4B3AiMB3oC/ROf6OZLZC0ByFA+C6gJTAK2MvMPsw1EVdUjuM4JUIxnSnMrFbLNzN7GHi4Fv2ymQlrhSsqx3GcEqGRef0VjXpftaRnYwbclWro87qkGZJa1FdOUyNmK36tIWT3/OvevP3CVcz88kH639RzmbYuu27BR0Nv5Mdx9zN44EVssN4aS9tWWGE5+t5wEtM+HcD4kXdxxgn7JTKfmTPncOqpV7HNNofQtetxDBo0LJFxXZ7LK0V50Oi8/opGPldzL9CNsBF2T2ajpPZAV+DOmNPJqR1/I0fwW6GYOu0nruvzDHt23ppWLVdYen71VdsysN/ZnPKP/rz42ih6n9Odh+44g84HXALARWcdwsYd1mbTXc6g3ZorM3jgxYz9cjKvvvlxXvO5/PK+LL/8cgwf/hBjx37DSSddTqdOHdhkk/Z5jevyXF4pygNfUdWHlwiuicdW0/5XQASF5tQSM/vMzMY2hOznBn/AoCEjmfHT3GXOd9t3B8Z+MYmnXxzBwoWLuPKWp9hq8/Z03GhdAI48eHeu+dfTzJw1j3FfTeG+R4dydPff5zWXiooFDBnyDmeeeRRt2rRi++23YI89duS557KlC8sfl+fyGrO8FI3MPb1o1PtqzOwX4CFgF0kd09sUfCiPAT42s1GSWkm6RdJoSbMl/SBpmKSuGe/bOBboOl1SL0kTY9GudyTtmDkHSbtKeimaF+dLGifpsow+h0n6MBb1miVpkKTMALWHYzGv30h6U1JFLAB2emzfJ45REa9hz7T3HhPnvH2W+fWJMldMO7d/vPbZ0XQ6RtKpae3LmP4k7RnHP1jS9ZKmxjGHSNok1/eUBJt3XJ/RYycufV0xfyHfTJzG5h3XZ5WV27Du2qsxZuy3S9vHjJ3IZh3Xz0vmhAmTadasGR06rLf0XKdOHfjqq29reJfLc3nlKS+F6vCvnMhX7Q6Ijz0yzv+eEJmcWk21AtoC1xDqmvwV+IpQNKtLlnHPBnYFzgCOAlYCXkpPBS/pz8CbwBqx35+Am4D10vqcAAwEvgUOAU4BNgXeyVSuQOvYdyDBpPku0EfSlXHcmwhp6+cDz0laK75vIDAtjr0USW0IyvpBM5sbz50CDIpjnBDl3A3UxlZwI7AWoa7LyYR4hedjTZgqSJqZ66iFTADatG7J7DkVy5ybPaeCFdu0ZMXWYftxVlr7rDnzadumVW2Hz0pFxQLatm29zLm2bdswb978vMZ1eS6vFOWlaKorqrx23MxsnKR3gaMlXZRWw6QH8DOhLglmNoNwYwaWFtwaQlBmpwLDMoaeDvzZzCz2nwa8A+wDPBFXbLcB44Bd4uoOYCjQP03G1YTo6YPSxhoGfA1czLKukq2As83s5djvXWB/4Dygo5lNjOf/B4wkKhkz+1lSf+BcSeeY2U9xvKOAlYE74/tWBq4HhphZepGw2jpOjDKzHqkXkgz4N/Bb4INajlEv5lUsoO2KGf8pV2zF3HkLmFsRth9XWrEVPyxctPT5nDz/w7Zu3ZK5c5dVjnPnVtAmTwXo8lxeKcpL0axZeTlJ1JYk1O69wPrAXrB0JXEI8JyZLU1UKKm7pOGSfgR+IRTW6gJ0yjLmCynFEhkdH1Mrj82BDYB70pRUJlsAawKPpI9lZpMJJZG7ZvRfRFCeqX5zge+AT1NKKpLaP0pfBfUlKP30/bqTgaFp+027AW1i3/rwfMbrzM9kGcxslVxHbQV/9sUktt58g6WvW7dqwa/bt+OzLyYxc9Y8pk77ia02r5zGVpttwNgvJtX+yrKw4YbrsXjxEiZMmLL03Oefj2fjjTeo4V0uz+WVp7xKmtXhKB+SuJrHgHlUmv+6AyuS5kQh6XDgcWAC4Wa+EyEtx6uElUwmmZl4U16DLeNjyje6prthKivv91naplI1k+9sM1ucce5nqmb1/TljLpjZFOApoKcCuwG/Ae5Ie19t5lwTuT6TvGnevBktWixP8+bNlnn+/OAP2Lzjrzhg3x1p0WJ5Lvz7QXwy9lu++Dr8J33kqbfodfqBrLJyGzputC7HHrEHDz3xVl5zad26JXvttTN9+jxCRcUCPvzwM15/fQTdumX+vkgGl+fyGrO8FE3V9Jf31ZjZHEKOpwMkrUJQRJNIW50ARwBfmtmRZva8mY0ws5GEfav68EN8rGnHPnVjXztL2zpUvfHny23AJoSV5SmEzyB9FVSbOTcovc44kJlfPsh5p3bjiIN2Z+aXD9LrjAOZPmMOh/e8hUvPO5SpY+5hh2025ujT+ix93xU3P8k3E6cx7p0+DHn8Em7t+0LerukAvXufzIIFC9lll6M455wbuPTSkwvq+uvyXF5jlgdNV1FpWQtbPQeRfk9wbLgBOBe42swuSmt/Hvi1mW2Zdm5r4L/AeDPbOJ7bGPgSON3Mbk/ruxzBNHexmV0Z96jGA3OBbbKZ/+Ie1eQ4/s5p59cBvgGeMLNj4rmHgX3MbI2MMT4BvjezdC+/1FxuMrNzM/qPBCqA3wFXmNmVaW0rE5TXO2b2xxo+y7eBBSmZ0cPwVeD/zOyFtH6pz+romMKkztSnSm8+1KfCr+OUL3Wv8Lvelr1r/X928ieXlY3rXyJq18zeItw0Uzfu+zK6vAhsEd2195B0EjCYkNywPvIMOJ2wvzVc0pGSuko6TtLdsc9i4AJgJ0nPSPqTpCOANwhms8vrIzsHtwG7E+LH7s6Y8yzgfGBvSS/HPbs9JJ0i6doCzMVxnDKjqa6okrya+wg36LfM7OuMtv6EbLrdCIHCPQnefu/UV5iZDSI4RMwiZON9iaAIJqX1uY+QOeNXhD2ku4AvgF3N7Kv6yq6BJwlFwJ40s2lZ5nwXcADB5Hk/wTTYk3oqbMdxmhbNmjWv9VFOJGL6cwKSDiU4l+xmZsMbej65cNOf4zQkdTf9tf/N1bX+Pzvx4wvLxvTXNJ3yE0bSpoT6LNcBb5eCknIcp/QoN5NebXFFlQwDgB0JgcA9GnYqjuOUK66onHpjZrs19Bwcxyl/VGaBvLXFFZXjOE6JoCaaQqlpXrXjOE4JEkJImx6uqBzHcUoEN/05TY753z7a0FNwHKcOuDOF4ziO07hx05/jOI7TqGmaCypXVI7jOCVDs6apqVxROY7jlApNU0+5onIcxykVzPeoHMdxnEZN09RTrqgcx3FKhmZNU1MV1OIpaVtJD0qaKGmhpNmS3pN0Qax4W5exWkjqL+l7SUskDYvne0k6oI5j3S9pQl3e4ziO0+BItT/KiIKtqCT1BG4HxgBXA+OAlsDOwBnAWsBZdRjyZOBEQqHB0YSCiQC9gGfjUVuuAP5Vh/6O4zgNT/PyUkC1pSCKStIuwB3AC8AhZrYorXmwpBuAXeo47BbAHDPrl8e8WpjZwiwViB3HcRo/ZbZSqi2FMv1dACwBTs5QUgCY2VwzGwIgqZmksySNkbRA0vRoLlw71V+SAScAbSVZPHrE8ysDf007f398z6Xx9W8lDZI0Gxgc26qY/mL/DyT9JGmWpPdjxd70PofHMY/POH+rpJ8l7ZB27teSHpP0YzR7fiLpbxnv6xLHO1TStZKmRvPoa7EYIxn995P0pqQ5kuZJej1dpuM4ZY7qcOQjJmzbPCtpSrzXfBa3WVpk9NsrbufMl/Q/Sf0krZKf9KokvqKS1BzYAxhpZlNq8Zb7gMOAm4BhwHoE09wwSduZ2TyCufBi4PfAXvF9X8fzrwJvxfcA/JAx/tPA/cCtQPMa5rEBYRX4HeFz6Qr8W1JbMxsAYGaPSuoK3CbpfTMbI6kbcCZwlpl9ED+DtYHhBGV9HjA1XmM/SWuZ2ZUZsq8D/gMcT1C81wGDJG1mZovjmMcB9wADgRvjtZwNvCnpd2Y2Jn1ASTNruFbi9eTq4jhOY6IIzhSSOgHvELZr/g5MJ9zTrwI2B46J/boALxG2XS4C1iXcu7aUtLuZLUlqToUw/a0BtAYm5OoYTYTHAKea2Z1p5z8CRhGq5d5hZu9J+gFYbGbvpQ3xg6TFwA8Z59O528yuyjUXMzsuTX4z4PV4LacRKvimOBPYCXhC0sEERfu8md2a1ucsoB2wjZmNjudejg4kF0q63czSFcloMzsmTf4vwOOEqsHvSmpDUORPmtkRaf2GAF8AlwDdc12j4zglTnEsf38h+BMcnLZNMlRSe+BwScdHS9n1wCfAYSmlJGkqMIRwP3osqQk1dJzzfoRVx0BJy6UOwsVPBjonIOOZ2nSS1FXSK5KmAb8AiwgrnE7p/cxsPnAosD7wITCHquXnuxKUz+iM8w8CrQiKLp3nM16n3tc+Pu4MrAI8nPE5/UxYhVb5nMxslVxHtR+G4ziNEmverNZHHqS2a2ZlnJ8V2xZLWg/YAXgofeVkZq8S7t0H5zOBTAqhqKYDFVTeZGuiXZzDj4QPIP1Yn7CiyZepuTpI2olgQhRwKrAr4Uu4l/DLYhnM7HPgXaAFMMDMfsroshrwfQ1zWT3j/I8ZrxfGx5TsdvHxOap+TkeSzOfkOE5jpw57VJJm5jqqkfIQMAO4S1IHSSvFLY6/AjdFxbRl7PtJlvePSWtPhMRNf2a2WNJQ4I+S1jGzmhTFdMKKajcqtXg6c5KYUi36HBbl/5+ZpZQEklbI1lnSKcCewAfAPyQ9aWafpXX5EVg7y1vXSWuvC9Pj4ylRpuM4TZEieP2Z2bfxx/uzwDdpTVeb2cXxeerH9owsQ8wAfpvknAoVR3UNwax3p6TuZvZLemPcc9kFeJkQB7W2mdXKRJeFhQRzWj4Ywdy3dAkraS2gW2ZHSdsCNxNixHoBI4HHJe0QzYIAQwkKbOsM899RwHyguv206hgOzAY6mdlddXyv4zjlQh2cKepr3o97UYMIVqEDgZmE7YULJC1JU1ZQ/UIgUU+tgigqM3tH0mnAbcAISf0IHiQtgN8RgnYfN7Ozojv5g5JuA94GFhA8/7oCL5nZkznEfQJ0lvQnYBow3cwm1HHKLxIcIB6R1J+wGro4jtc21UlSW4KTw2fAuWa2MLqwjwD6EAKSAW4hLJNfknQJMIWwr3UQcGGGI0VOzGyupL8D90hajbDv9iPBJLg9sDDjj8dxnHKkOM4U1xLue9um/fgeprCau0TSACqtQpnbGBC2PrKttOpNwTJTmNldkkYQFMDFhJvqAsJN/l9AamVwHGGFcSLBFXIJYTPuTSqdCmribKAv8CRhT+cBqjo35Jrr6zHG6XxCkPJEgjt7O6B3Wte747l9UybC6KJ+BnC3pKFm9qiZ/S96NF4L3ACsSPDOO9HM7qnL3NLmeJ+kbwnu7gMIq8hphBVd3/qM6ThOiVGcgN9tgc/SlFSKkQSfgk5U7k1tSfDyS2crgnt7YshjaZo0/uU7TsNRZ62zcfeHa/1/9qsnjqqXVos+BpsBG5lZRdr5K4F/AlvHH+gfEBTXDmnu6X8AXgMON7OB9ZGfDc+e7jiOUyoUZ0XVh7C98IqkWwlu6V0IFqfX0pIL/IOwmno0bpmkAn5HAE8kOaGGjqNyHMdxaksRUiiZ2bOEDEALgTsJYTEHErL/HJDWbyiwP7AhYZ//5vi4byqjTlL4ispxHKdEsCLVozKz1wgmvFz9BhNzqBYSV1SO4zilQhPNnu6KynEcp1RomnrKFVXT5ouGnkBBabVB79ydEmT+t5cVVZ5T6nSs+1vyy+FXsriichzHKRV8ReU4juM0aorkTNHYcEXlOI5TKriichzHcRoz1jT1lCsqx3GcksGdKRzHcZxGjZv+HMdxnEZN01xQuaJyHMcpGZpoZoqS08+SekiytGO+pLGSrpW0cgHlTohFHuvyntRcN6xunNr0cRzHAYLpr7ZHGVFyiiqNY4CdCdl7XySkoB8iqVDXlMoeXBdeJMxxahFkFZyZM+dw6qlXsc02h9C163EMGjSspOT1/OvevP3CVcz88kH639RzmbYuu27BR0Nv5Mdx9zN44EVssN4aS9tWWGE5+t5wEtM+HcD4kXdxxgn75TWPFKX+ebq84soDMKnWRzlRyqa/MWb2UXz+uqQ1qVRew5MWZmb/rcd7fgB+KIasYnD55X1ZfvnlGD78IcaO/YaTTrqcTp06sMkm7UtC3tRpP3Fdn2fYs/PWtGq5wtLzq6/aloH9zuaUf/TnxddG0fuc7jx0xxl0PuASAC466xA27rA2m+5yBu3WXJnBAy9m7JeTefXNjxvV9bm88pYHwHLlpYBqSymvqDIZER/bA0jaTNITkqZLWihptKQj0t+QZnb7vaR+se+Pkp6WtG5G3yrmOEnrS7pH0iRJP0uaLOnRlAkym1mvNtRgHsw5z0JRUbGAIUPe4cwzj6JNm1Zsv/0W7LHHjjz33BslI++5wR8waMhIZvw0d5nz3fbdgbFfTOLpF0ewcOEirrzlKbbavD0dNwof7ZEH7841/3qambPmMe6rKdz36FCO7v77Rnd9Lq985S1Fqv1RRpSTovp1fPxB0tbA+4SCXqcD/0dQZI9IOj7Le+8D5gNHEKpWdgEerkmYpF8BIwmmx+uAfYBzgF+AFfO7lGqp9Twlzcx11EXwhAmTadasGR06rLf0XKdOHfjqq2/rfTGNRd7mHddn9NiJS19XzF/INxOnsXnH9Vll5Tasu/ZqjBlbKXfM2Ils1nH9vGSW8+fp8gpIE92jKmXTX3NJywFtgL2BnsBk4D/A88B0oIuZzYv9h0haA7hK0n1mtiRtrBfM7O+pF5JWBa6XtI6ZVbe/dDmwKrClmX2Zdn5gEhdXDfWZZyJUVCygbdvWy5xr27YN8+bNL3l5bVq3ZPqM2cucmz2nghXbtGTF1i0AmDWnYmnbrDnzadumVV4yy/nzdHkFpLz0T60pZUU1MuP1COCE+Lwr0AdYGJVZipcIpZQ7AZ+lnX8+Y6zR8bE91TtC7AO8lqGkCk2t52lmq+Qe7gurreDWrVsyd27FMufmzq2gTZ437MYgb17FAtqumHHTWbEVc+ctYG7FQgBWWrEVPyxctPT5nDxvSOX8ebq8wlGsCr+NjVI2/R0J7ABsC6xlZjuZ2SfA6gQFfDawKOPoH9+7RsZYP2a8XhgfW9Ygfw1gUr1nXz/qM89E2HDD9Vi8eAkTJkxZeu7zz8ez8cYblLy8z76YxNabV47bulULft2+HZ99MYmZs+YxddpPbLV55Qb5VpttwNgv8vvqy/nzdHkFpIma/kpZUX1mZiPN7KPoXZfiJ2AJMICgyLIdSXjV/QDkt1FRQrRu3ZK99tqZPn0eoaJiAR9++Bmvvz6Cbt26loy85s2b0aLF8jRv3myZ588P/oDNO/6KA/bdkRYtlufCvx/EJ2O/5Yuvw03okafeotfpB7LKym3ouNG6HHvEHjz0xFuN7vpcXvnKW0pz1f4oI0pZUWXFzCqAN4FtgI+iMss85iQgajCwp6SNExirJOjd+2QWLFjILrscxTnn3MCll55cUFfcpOX1OuNAZn75IOed2o0jDtqdmV8+SK8zDmT6jDkc3vMWLj3vUKaOuYcdttmYo0/rs/R9V9z8JN9MnMa4d/ow5PFLuLXvC3m7phfi+lxeecsDmqzXn8xqvU3RKJDUg+D9tm1aHFVmn60JThVjgH7At8AqhL2pHc3s4JrGktQFeAPoambD4rkJwDAz6xFf/wr4IL7lKiBldtwf+KeZTU4bv4OZTahmnLr0yTnPulH7PapSxEvRO42bjnXWJhte9kqt/89O6P3HstFWpexMUS1mNlrSdsAlBNfxNYAZwFjg8YRkfCdpR4L33z+B1YD/EVZzSazYHMdxlqXM9p5qS8mtqJwk8RVVkviKyqkbdV9Rtb/y1Vr/n5140V5lo9XKckXlOI5TlpSZk0RtcUXlOI5TKjRR058rKsdxnFLBFZXjOI7TqGmaesoVleM4TqnQVFMouaJyypZie+EV08vQPQybKGUWyFtbyi4zheM4TtlSxBRKkrpIGhLLAlVI+kzS3zL6HCHpY0kLYl2+ayUlnnvUV1SO4zglQrMiLS0k/ZWQL/Vu4GZCUu9OwAppfY4CHgLuAv4ObEZIsLAh8Jck5+OKynEcp0QohuUvpoe7C7jQzK5Pa3o9rU9z4AbgeTM7JZ5+Q9IioL+kW8xsBAnhpj/HcZwSoUg5aVNV0G+roc9OwNrAAxnnHyGsvg7OawYZuKJyHMcpESTV5ZiZ66hGzO8JeVEPkjRO0uK0/aeU6W/L+PhJ+htj9Yqv09oTwU1/juM4JUKR9qjWjcdtwMXAp8AewAXArwhFa1ePfWdkef+MtPZEcEXlOI5TIqgOisrMVqmnmGZAW+BwMxsYzw2T1Ao4V1J6HEZ1SXITTXjtpj/HcZwSoUh7VD/Gx1cyzr8cH3+b1ifbymk1sq+06k0iikpSD0mWdsyXNEXSq5LOkrRyEnJqkN9S0qWxkGDJI2mYpGENPQ/HcRoXzVT7Iw/GVHM+NeoSgjkQMvaiJLUGNiJj7ypfkl5RHQPsDOxN8OSUEnUAACAASURBVKv/ilD99lNJ2yYsK52WQG+gSwFlFJNT4tGomDlzDqeeehXbbHMIXbsex6BBw1xeDfT86968/cJVzPzyQfrf1HOZti67bsFHQ2/kx3H3M3jgRWyw3hpL21ZYYTn63nAS0z4dwPiRd3HGCfvlNY8Upf55NnV5ULQV1dPxMfMPbz+CSe8D4D3ge+DojD6HA8unjZEISe9RjckoD/+4pDuBt4DnJW1iZgsSlll2mNlnDT2HbFx+eV+WX345hg9/iLFjv+Gkky6nU6cObLJJe5eXhanTfuK6Ps+wZ+etadVyaZwkq6/aloH9zuaUf/TnxddG0fuc7jx0xxl0PuASAC466xA27rA2m+5yBu3WXJnBAy9m7JeTefXNjxvV9bm84sqD4sRRmdlgSS8Dd0hag0pnijOBvmY2McxFvYD7Jd0OPEllwO+TZvZeknMq+B6VmY0BrgTWBw4DkLSXpOclTY5mwi8k3SxppfT3Sro/ulFuK+k/se9USVdLWi722RD4Kb6ld5r58dK0cQ6RNDK+f5akFyRtXY2srSS9EVOGTJb097Q5fxDPfyJp77T3HhFl7pR5/fG65qRfm6T9ooxZkubF8U5Pa1/G9BdTmZikQ6OL6FRJsyW9JmnTOn4l9aKiYgFDhrzDmWceRZs2rdh++y3YY48dee65N1xeNTw3+AMGDRnJjJ/mLnO+2747MPaLSTz94ggWLlzElbc8xVabt6fjRusCcOTBu3PNv55m5qx5jPtqCvc9OpSju/++0V2fyyuevBTNmqvWR550B+4BzgdeAg4FLgJOS3UwswcIK6rOhP2si4C+BMtaohTLmWJQfEz9b9sIeBs4CdgHuB7oBryY5b0tgGeB52KffwO9gD6xfSqwV3w+gGB63JnwISOpB/AEMJnwYZ8MbAwMl7RZFlmPxf7dgP8At0i6CvgXcCtwEDAHeFbS2vF9T8R5LGOui14yPYCHzWx2PPe3eJ2LgL9FOf2A2vwMu47gNno84bPrCAyKUeLLkEcMRVYmTJhMs2bN6NBhvaXnOnXqwFdffVuXYVwesHnH9Rk9duLS1xXzF/LNxGls3nF9Vlm5DeuuvRpjxlbKHTN2Ipt1XD8vmeX8eTYFeSmKZPrDzOaZ2blmtp6ZrWBmG5vZtWa2JKPfw2a2lZm1iH3PN7P5+UmvSrHc01Pf3joAZtY31SBJwHBgHPCWpN+YWbqNoyVwtZn1i6+HxA27kyRdbWaTJI2MbZPSl5ySmgHXACOBA8zM4vk3CEFplxBsqumy/mFmg2K/4cD+wHlAJzP7Jp6fCnwEHAjcZWaLJPUFLpR0tplNj+MdAawK3BHf1xa4ERhqZnunyX2tdh8jo81s6a8VSb8AjwM7Au/Wcox6UVGxgLZtWy9zrm3bNsybl/jfZNnLa9O6JdNnzF7m3Ow5FazYpiUrtm4BwKw5FUvbZs2ZT9s2rfKSWc6fZ1OQl6KJJk8v2opqmY9X0tqS+kgaDywgrC7e+v/2zj1Oq6r6/+/PDMgwgCBaKCAXI0IUv1RqYpmAmWkZ4i2/onnpIorihUojFb4oUWqhlKmkyU+8J941gryghZKo5IUBvAEiiHKVYZiBcP3+2OcZzhyfYZ5nrs/Ms96+zuvM2WefvdY545wPe+119o5O90tz/b1pjguAb9Rgtx9hmo+7UiIFYGargNnAkET97cDMWL0ygsguSolUREm0j/eCbiHc59mxsnOB58wslQFzKOH7hJupHY8mjl9L4wcQvqGoacvGcHFxEaWlZVXKSkvLaFfHF2g+2ttcVk6H9omXXPu2lG4up7SsAoBd2++wu2v7tmyq4wuwJT/PfLCXorF6VLlGYwlVj2i/MurlzCaMV11HyBA8iBBSA0j+psvNbGOibHW0r+nr587R/sM051alub7UzLYlyraS+CbAzLZGPxbFylYTejcjJRVE41VfJepNRaRSu1bU4Hd1rE0cVyT9aCh69erG9u2fsnTpysqyRYveo0+fHju5yu2lY+GSFRzQf0e7xW3bsE/PLixcsoINGzezavV6BvTf8W+PAfv2oGRJbf+XCbTk55kP9lI0Unp6ztFYQnVstJ8DDCDk3v/czG40szlmNh/4pJpri9J8h9Ul2idf3ElS5/dMc26vDK7Plj8AvQnjbucRxPCh2PmPo33dBhyagOLiIo48chBTptxFWVk5L7+8kKeemsewYclOqdtLUVhYQJs2rSksLKjy86MzX6J/37057uiDadOmNWMvOp43Spaz5J3w0rtrxnNcdsFwOnVsR98vdOWsU4cy/a/P1WCt8e/P7TWevRTeo2ogJO0P/IoQQrufHVNrbE1U/clOmkmubXIKIUz3z+g41bNI9sYWAyuBEdFYWMqnLsC3gKczuIWMMbN/A/OAsYTEjamJHtpcgiCPTHN5zjNu3LmUl1dw6KGnMWbMtYwff26DpuI2d3uXjR7Ohrfu4OejhnHq8Yex4a07uGz0cNas28T/jpzM+J+fzKrXb+WggX04/fwplddd9fsHeHfZahbPncKs+6/k+psfr3NqekPcn9trXHvQqFl/OYViQze1byRk1t1OSEt8i5CksRch9/6HhNDZ98zsP5JaA4uAQkL23kZCQsMgQjbeWWY2LWp3GiFE+BGht/If4CjgEkISw6iYD+8QBGs0sAFYaWYrJZ0O3EEY3/kz0J6QRNEdOMjMFsdsHZccu5G0ANhgZoMT5QbcYGYXJcpHAHcC/wV6mtnKxPmfEsazZhMyE9cAXwL2MbOfR3WeBUjZjGbceAYYbmYPx9rqBbwXf2bZsaRe5+PKd3wpeic7+matJgf/9Z8Z/83++6RvtBi1qu8e1R2E7LOnCOnjfQi59fulMvmiHsb3CVl3fyasELmd6leErCCkcA8jiM1phDTtCxP1fkpIyniC8OX0TyN70wlro3QlfJQ2FXgX+HpKpOqZGZEfDyVFKvJnKuH+2wJ/IaTunwcsS9Z1HMeJk6+hv3rpUTUU1fVychlJxxPEarCZzWlqf3aO96jqE+9ROdmRfY/qkBmZ96hePKHl9Kh8mY96IpohoidheeYXc1+kHMdpbrS0bL5McaGqP24hfCf1Mg0whYjjOE7BZ+agyQ9yWqjM7Mym9iFTkskWjuM49U1LG3vKlJwWKsdxHGcHylOlcqFyHMdpJuSpTrlQOY7jNBdcqBzHqRONmTLemKnw4OnwuYILleM4jpPTtGqs2VlzDBcqx3GcZkKB8vMbfRcqx3GcZoJ/8Os4juPkNHka+XOhchzHaS546M9xHMfJaTz05ziO4+Q0rVyoHMdxnFxGeRr6q7exOUlnSrJo+1qa80WSNkTnp9WXXafx2LBhE6NGTWTgwBMZMuRsHnvsWbeXQ/ZGnvFt/vn4RDa8dQdTfzeyyrnBX9+PBU9fx9rF05h57+X06LZH5blddmnFzdeew+o3b+O9+Tcx+sfH1MmPFM39eeaaPQihv0y3lkRD9Kg2AWcB8xLlwwmr2lY0gE2nEZgw4WZat27Fv/41nZKSdznnnAn069ebL36xp9vLAXurVq/nt1Me4luHH0Dbol0qy3ffrQP33nIJ5106lSf+8QrjxpzE9BtHc/hxVwJw+cUn0qf3nnzp0NF0+VxHZt57BSVvfcDsOf/JqfvLd3uQv1l/DXHf9wOnSGqbKD+LsJR8eX0YkdRKUp6uztL4lJWVM2vWXC688DTatWvLgQfux9ChB/PII8+4vRyx98jMl3hs1nzWrS+tUj7s6IMoWbKCB5+YR0XFNq6ePIMB/XvS9wtdARhxwmFMuuFBNmzczOK3V3L7PU9z+knfzLn7y2d7KQpkGW8tiYYQqrsJPafhqQJJewNHANOSlSXtJul6ScskbY32V0lqHavTKwoZ/kzSBEnvA1uBvaPz+0i6T9JaSRWS3pD004SdwVEbp0q6UdIaSWWSnpTUJ41fx0iaI2mTpM2SnpJ0UKLOtCic+SVJMyWVSnpf0u8ktUnULYp8fyvycZWkP0rqkKhn0fM4T9LbkrZImifpy5E4T5C0QtJGSQ9J+lzGv5k6sHTpBxQUFNC7d7fKsn79evP228vdXo7b69+3O6+VLKs8LttSwbvLVtO/b3c6dWxH1z0783rJDruvlyxj377d62SzJT/PprCXopUy31oSDRH6W0foOZ1FEC2AM4DVwMx4RUntgeeBzsBEYBFwMHAl0As4PdH2xcBrwKjoeL2kPYF/AZ8CPwdWAT8AbpH0eTO7OtHGb4EXCKvwdgF+DTwraYCZrY/8Ohu4FbgXuA4oBC4B5kj6mpm9HmuvDfAIMJWwDP3hwOXARmBC1F4h8ATwVWASMB/oB1wFDJA0xMw+jbU5HHgPGAO0Bn4XPdNZwC7Ajwki/fvI7nASSNqQLEtitrimKpWUlZXToUNxlbIOHdqxefOWjNvIBrdXf7QrLmLNuk+qlH2yqYz27YpoXxz+PbVxU1nluY2bttChXTIgkh0t+Xk2hb0ULW3sKVMaKuvvduAxST3MbDlwJjDdzLYnFv4aTXhhf8XMXovKnpJUBlwv6Tdm9mas/ibge2a2PVUgaSxBcAbG2vibpI7AWEl/NLP4S/t9Mzs5dv1C4EVgJDBJUjuCMDxgZqfG6s0ClhBE9KRYe0XAL83soZj/XwVGEAkVQTiHAt81sydj9VYADwNHE4QsRSFwtJltiWy3I/RGe5nZETGf+gMXSio2szIakOLiIkpLq5ooLS2jXR1faG6v4e1tLiunQ/vES7V9W0o3l1NaFoaMd23flo8rtlX+vKmOL9yW/Dybwl6KlhbSy5SGGpv7O6EHdYakbwJfIE3YDzgGWAAsjMJarSS1Av4WnT88Uf/RuEhFDAFei4lUijsIIchDEuX3xQ/MbB6wNGZrENAJuDPh01bg2TQ+baeqyEDo9cVHVI8B1gKzEm3Ojq5Ptvl0SqQiSqJ90k4JIKBHohwz61TTlrxmZ/Tq1Y3t2z9l6dKVlWWLFr1Hnz6fMV0vuL36Y+GSFRzQf0e7xW3bsE/PLixcsoINGzezavV6BvTf8b/rgH17ULJkRZ1stuTn2RT2UuRr1l+DCFUkJtMJPamzgXlmVpKmahdCOGxbYkvFpPZI1F+Vpo3OwIdpylN1d0+Up6u7OlavS7R/JI1fI9L4VGpmWxNlFYSQYIouUfvJ9jYTek/JNtcljrfWUF6U5p7qleLiIo48chBTptxFWVk5L7+8kKeemsewYUPcXo7YKywsoE2b1hQWFlT5+dGZL9G/794cd/TBtGnTmrEXHc8bJctZ8k54yd414zkuu2A4nTq2o+8XunLWqUOZ/tfncu7+8tleioIstpZEQ37wezvwC8JY03nV1FkDlAI/qeb8ysRxun7vWmDPNOV7xc7HSVe3CzvEcU20Pw94qRq/smUNQQy/t5PzOc+4cecyduwNHHroaXTq1IHx489t0FRct5cdl40ezuUXn1h5fOrxh3H15AeYOHkG/ztyMpMnnMlfbhjFS6++zennT6msd9XvH2DKxB+xeO4UtpRv4/c3PVrn1HRo/s8z1+xB04X+JI0HxgH/MbOBiXNHEsbb/4cwPPMQcGliyKVu9s3q58YlnUkQpy+b2YKo7E+EQf/TU05Hg/wPm9mZkq4kJED0N7P3d9J2L0JywcVmdn3i3CTgUqqOUSFpBmHsp6uZbZA0GHgGeMHMDo3V+xphjGqsmU2KEjw+AKaZ2YU13PM04LhkGC31SzUzRcc/JIQ+DzSzV2po04AbzOyiWNlA4FXgLDObFis/k8Qzz44l+RnwbgH4Cr8tgb5ZB+jOf+GZjP9m/zhoSL0EACXtR0gA2wh8GBeq6L06mzDWfgvQlZCw9i5wWCJJrNY06BRKZlZdTyrFZOBE4HlJk4E3CFltvQjjOueb2bLqL69s4wzgyUj4VgInA8cTxCep6t0l3UcQjj0JWX8fADdHPpdKugi4VVJnwr8O1hJ6XQcCFWZ2RQ0+JbmbkGU4M7rPlwm9w72BowjCNDfLNh3HyTMaO6QnqQC4jZAFPYAwfh/nGsJ7+wcpUZK0ipChfBKJnIDa0qRz/ZnZJklfB34FnE9ICthM6D39nc+G7dK18ZGkQ4HfENLD2xOy835iZremueQy4BuEMbRiYA5wQSo1PWrzdknLCb292whJGasJ/6q4uRb3+V9JxxBS3EcQutAVwHLgqeh+HcdxdkoThP4uBroD3yZ8IlOJpG7AQcCYeM/JzGZL+gA4gXoSqnoL/eU6sdDfcDN7uIndyRE89Ndc8dBfSyD70N8l857O+G928iFHbKypzs6yfyXtA7wOjDCzhyU9C3RKhf4kHUX4NvYoM5uVuPZvQE8z65+pvzujpSWHOI7jtFgaK+tP4YPXPwN/38k/7FOZ0sls5FRZMuO61vgyH47jOM2EbL6PyvZbyQQ/IYzJZ9Ijqq6XV28Rm7wRKjN7lvBxrOM4TrOksKDho/WS9iAkSUwCNktKCV4roDA6LmdHDkG6nlNn0ve0aoWH/hzHcZoJjRT66w50JAjV+tj2dWD/6OfxQGp6u/3TtDGAkA1YL+RNj8pxHKe500hZf28TpqZLcj0hq/rHwHIzWyFpPjBC0vWx9PQjgG7Ag/XlkAuV4zhOM6Ex5vAzs1LCvKZVSK3IEA2jpLiU8M3UPZKmsuOD33nAX+vLJxeqvKZvUzvg1JIty+9pahecJiDXJps1s6clfQ/4P8Kk2ZsIs1T8Is0E4rXGhcpxHKeZ0LoJl/kws8HVlM8ksdZgfeNC5TiO00zItR5VY+FC5TiO00xwoXIcx3FymkIXKsdxHCeX8R6V4ziOk9M01cKJTY0LleM4TjOhdZ72qLKaaUPSlyXdIWmZpApJn0h6UdIvJXWM6li0wm2TIGkfSTMlrU/5Imlw9PPgpvIrGyQtjVYPdhzHqaRAmW8tiYx7VJJGAn8krE/ya2AxUAQMAkYDnycssjUIWFHvnmbO74GBhFV/P4p86dOE/tSG4cAnTe2E4zi5hYf+dkK0gu6NwOPAiWa2LXZ6pqRrgUMBzOzF+nIuWga5MGGvJvYD/m1mlatRSmpWQmVmrza1D47j5B75mvWXaejvl8CnwLnpRMPMSlMrPCZDf1HozaKw4ROSSiWtlXSTpPbxdqJ610u6UNLbwFbCjL1I+pykP0v6UNJWSW9JGiupMDo/WJIRek/HRm2ZpF7V3ZSkkZLeiMKYayTdE68vqa+kTZKmJ647Lmr7wljZKEnPS/o4uscFks6NxDZp97QoZLo5Cp/Ol3RK7HyV0J+kMyN735R0S+TrWkkPSupa3f05jtOy8NBfNURCMBSYb2Yr62DrQeAuwgy8XwOuBHoA303UOwn4ALgMKAXekVRMmCSxG3AFIex4FHA1sA9hNt9XCGHHh4AlhMkSAVYBvdLc13hgHHA78LOo7auAuZIGmtlHZrZE0jnAXZLmmNmtknoCfwEeMbMbYk3uA0wHlgLbgYOB6wiTNF4Rs/trgvDfS1jzpZQQqvyMj2m4HXgMOJXw7K4B7iT8fpL3t6GmxszyM4zgOM2VVnm6MFMmob89gGLCC7gu3Gtml0c/z5a0DfiNpEFm9kKsXhHwLTOrHKORdC5hpcljzezxqHhW1Fu5SNJ1ZrYIeFFSBbA+HoIMqyrvQNJuBCGbYWZnx8pfJgjexQQxwczuljQEmCLpFeBPhPGjs+JtmtmYWDsFwBygELhE0pVmZpL2iez+xcx+FLt8Vs2PD4DHzeyixH1cI2kvM1uVYRuO4zRTCvN0jKox9fneao4PT5Q/FRepiCHAxphIpbgj2g/O0pdBBEG8M15oZgsIySLJtVhGE9ZomQt8BTjFzNbHK0j6ShSKWwlsi7argN0IiSYARxKe+c1Z+pvi0cTxa9G+Z7KimXWqaaulD47jNBGNtHBizpHJ/awBykjzMsySDxPHq6N9chnjdD2Dzmmuj9dNtxTyzuhcjU+pNqu0Z2ZbgPuANsA/kwkj0bjW88DehDDiYcBBwMSoSttov0e0r21W5NrEcUW0L6ple47jNCPydYyqRqGK1hR5GjhQ0l51sLVn4rhLtE++fNP1bdemuR5gr9j5bEjVr67NKu1J+jJhnOklYIikMxPXDCOER08ws7vNbK6Zzeez9/JxtO+epb+O4zguVDUwiTDe8idJnxnXktRO0pE1tHFKNcdzMrD/NNAxWqArzmkEMXgmgzbivACUR9dXIukAYEBkL1W2K3A/8AbwDULCxB8l7Ru7NCVIW2PXFQGnJ+zOJiRajMzSX8dxHAplGW8tiYy+ozKzuZLOB/4AzJN0CyHzrg0hg28k4WU+eyfNnCLpv8BzhIy4ccCTiUSK6rgDOB+4U9IVwCLg24Skhz+b2eJM7iN2PxskTQSukvSXyPeuhDGllcDkWPWphDGm75jZ1iix4yDgPkkHm1k54b63AXdLugboAIyJyuJ235P0W2BslMl4P2FFzAOAXczsN9nch+M4+YVn/dWAmd0kaR5BHK4ghO7KgYXADcBNNTQxnDCjxcWEnsfthJd5Jra3RJl3k4BfEcaYlgFjCSnaWWNmV0v6GLgAGEFIE58FXGpmH0PlbBw/ICRPvBNdt1nSycA8wn2fY2Ylkk4iCN1DhPG3WwnjXbcm7P5K0lsE4b2bIGaLgWtrcx+O4+QPLS2klylq6G9pYt8r7WZmNX7b4zQqLSs+4DjNi6xlZ+aKv2X8N/ud7ke3GFnz2dMdx3GaCT7Xn+M4jpPT5OkQVcPft5mNNzN52M9xHKdu5Gt6uveoHMdxmgmtCzz05ziO4+QwLa2nlCkuVHnNkqZ2wHHS0rbHuEa1t2X5/zWqvUDfrK9woXIcx3FymnxNpnChchzHaSbIe1SO4zhOLuOhP8dxHCen8dCf4ziOk9MoT2emyFeBdhzHaXYoi63WNqQjJE2TtFhSmaQV0erlA9LUPVLSi5K2SPpI0i2S6n31cBcqx3GcZoKU+VYHRgI9CMsdHQ1cEh2/JOmQHb5oMPAk8D5wLGF18+8DT0iqV23x0J/jOE4zoZFyKUaZ2UdV7EqzgPeAnwMnRMXXEBaU/YGZfRrVW0VYLukk4L76csh7VPWIJMtw65VhvfGxtpcmzn0o6R+ShjbW/W3YsIlRoyYycOCJDBlyNo899qzbc3uNZm/kGd/mn49PZMNbdzD1d1UXyR789f1Y8PR1rF08jZn3Xk6PbntUnttll1bcfO05rH7zNt6bfxOjf3xMnfxI0djPE6BQmW+1JSlSUdkG4C2gO4CkboQFZKenRCqqNxv4gB1iVi94j6p+GZQ4/i3h8/PhifJVaeqmaEVY0bgboVsd5zngUsI/rHoDlwP/kDTMzB6rg98ZMWHCzbRu3Yp//Ws6JSXvcs45E+jXrzdf/GJPt+f2GtzeqtXr+e2Uh/jW4QfQtmiXyvLdd+vAvbdcwnmXTuWJf7zCuDEnMf3G0Rx+3JUAXH7xifTpvSdfOnQ0XT7XkZn3XkHJWx8we85/cur+MiGbkJ6kGicCN7OMxpMkfQ7YH7gnKto/2r+RpvrrsfP1gveo6hEzezG+AeuBimS5maUrS10znCBCF5jZvxMm1kf1XjCzu4HvEURrdEPfW1lZObNmzeXCC0+jXbu2HHjgfgwdejCPPPKM23N7jWLvkZkv8dis+axbX1qlfNjRB1GyZAUPPjGPioptXD15BgP696TvF7oCMOKEw5h0w4Ns2LiZxW+v5PZ7nub0k76Zc/eXCY2RTPEZm5KAqQS9uC4q3j3ar0tzybrY+XrBhSqHiJa4vwS4zcym1lTfzN4FPgYa7p9wEUuXfkBBQQG9e3erLOvXrzdvv73c7bm9JrXXv293XitZVnlctqWCd5etpn/f7nTq2I6ue3bm9ZIddl8vWca+fbvXyWZjP88U2QiVmXWqacvQ7LXAccBIMytJnKsuX75e8+g99JcjSOoP3Aa8BIzK8JqOQGdC7Dh5LoNu/+KM/SsrK6dDh+IqZR06tGPz5i0Zt5ENbs/tZUq74iLWrPukStknm8po366I9sVtANi4qazy3MZNW+jQrm2dbDb280zR2DNTSJoIjAEuNLNpsVNro326nlNn0ve0ao33qHIASbsCDwJbgBPMrKL6qmolqbWkLxLGsgrZETduMIqLiygtLatSVlpaRrs6/sG7PbdXVzaXldOhfUI02reldHM5pWXhT2nX9jvs7tq+LZvqKCiN/TxTNGboT9IEYCzwCzObkjj9ZrRPNxY1gPRjV7XGhSo3mAb0IaR5vr+Tet8HtgFbCWt0DAEmADcmK9Zjtx+AXr26sX37pyxdurKybNGi9+jTp0c2zbg9t1fvthYuWcEB/Xe0W9y2Dfv07MLCJSvYsHEzq1avZ0D/HdHxAfv2oGTJijrZbOznmaJAlvFWFySNA64ArjCza5PnzWwFMB8YEf9mStIRhESwB+vkQAIXqiZG0mWEBIpLzaymkdg5hJTQA4F+QGczG2dmDT6vSnFxEUceOYgpU+6irKycl19eyFNPzWPYsCFuz+01ir3CwgLatGlNYWFBlZ8fnfkS/fvuzXFHH0ybNq0Ze9HxvFGynCXvBBG5a8ZzXHbBcDp1bEffL3TlrFOHMv2vz+Xc/WVCY3zwK2kMMB54nJBVfEhs+3Ks6qXA/wD3RLNZnA5MB+YBf629B2l8aoR3XN4i6WFgoJn1qub8EcDfgQfM7JQa2loKLDCz4+rPwyVZ/fI3bNjE2LE3MHfuAjp16sCYMWdw7LGD688dt+f2ItItnPiri0/g8otPrFJ29eQHmDh5BkO+sT+TJ5xJj+6f46VX3+YnY25i+Yo1QPiOasrEHzH8mIPZUr6N39/0KFNurfrlR20WTqz78+ybtZws3fRYxn+zvTocWyu5kvQscHg1p5fF32eSvgP8H0GwNgEPE0KF62tju1qfXKgajp0JlaS9gVeAD4FDzGxzDW0tpYmFynEaizxZ4TdrIVlWmrlQ9WxfO6HKRTzrrwmQEdGkaAAAB/FJREFUtAswA9gNuBgYoPR99U/MbGFj+uY4Tu7SYpQnS1yomoauhLEmCDHd6pgDDG5wbxzHaRb4wolOvVNdmM7MlpLlP46qG+dyHCd/cKFyHMdxcpo81SkXKsdxnOZCvq7w60LlOI7TTPAeleM4jpPT1HHl3maLC5XjOE4zobCpHWgiXKgcx3GaCd6jcvKQvk3tgOOkZcvyBl8QoJmSn0rlQuU4jtNMkAuV4ziOk8vEVtTIK1yoHMdxmg3eo3Icx3FyGOXpEoIuVI7jOM0ED/05juM4OY6H/hzHcZwcJl+z/urcj5T0ZUl3SFomqULSJ5JelPRLSR3rw8mGRIFJkt6XtD1aSbepfTpTkknq1dS+OI6TOyiL/1oSdepRSRoJ/BF4Hfg1sBgoAgYBo4HPE1awzWWGAZcB44HZwJYm9SbwBOEZrmpqRxzHyR2k/JxEqdZCJelQ4EbgceBEM9sWOz1T0rXAoXX0rzHYL9r/wczWNaknEWb2MfBxU/vhOE6u0bJ6SplSl9DfL4FPgXMTIgWAmZWa2SwASaMkPS/pY0mlkhZIOleJFBZJQyU9K2mdpC2Slkq6N3Z+cBQSO1XSjZLWSCqT9KSkPkkfJI2U9EYUklwj6Z54OC0K810dHa6N2h4fnav8OdHmUknTYsepMN03Jd0S2Vkr6UFJXRPXtpH0O0kfRn7/U9LBO2mzV8Luw5K+K+nV6PksknT2Z34zjuO0SDz0lwUK/c+hwHwzW5nBJfsA04GlwHbgYOA6oCtwRdRmb+BJ4GHgGkIIbm/gu2na+y3wAvBDoAsh7PispAFmtj5qbzwwDrgd+BnQDbgKmCtpoJl9BAwHRgE/Ao4ESoEVmT+JKtwOPAacCvSI7uFOwnOK1zk58v9ZoD/wALBrhja+ErU7CfgI+Alwm6S3zey5eEVJG2pqzCw/F2FznOaLp6dnwx5AMUF4asTMxqR+jnpRcwgz1l8i6UoLb8yvAm2Ac8xsY+zyO9I0+b6ZnRxrcyHwIjASmCRpN+BSYIaZnR2r9zLwCmHc7Jdm9qqklDDNN7MaX+474XEzuyhmazfgGkl7mdkqSfsB/wtMNLPLo2qzJa0GMp2Bc3dgkJl9ENl4jiCEI4Dndnah4zjNn5bWU8qURpFnSV+JQmErgW3RdhWwGyHhAuBVYCvwgKSTk2GzBPfFD8xsHkE0D4+KBhGSOu5M1FtASPwYUqcbSs+jiePXon3PaJ/y7f5EvQeA/2Zo45WUSAGYWTmwJGaD2LlONW0Z2nQcJ0eQlPHWkqitUK0BykjzgkwSjbM8Twjj/Qw4DDgImBhVaQtgZu8A3wIqCCGyD6LxpR+mafbDNGWrCT0OgM47qbcqVq8+WZs4roj2RQmfVscrmdl/01ybqY2UnaI05Y7jtDBEYcZbS6JWQmVm24GngQMl7VVD9WGEMOEJZna3mc01s/nAZwZIzOx5M/seYczmEGAR8P8kHZGoumcaO13Y8SJfu5N6e5GZMFQQQpFJaityKZtd4oWSWtWhTcdx8gplsbUc6hL6m0QYZ/pT9LKtgqR2ko5khyBtjZ0rAk6vrmEz2x6F81JjWwckqvwgYetrQC/C2BeERIty4LREvQOAAQSRrYmlSbuShgLtM7g2HakxpJMT5SfiM4Q4jpMB+Rr6q/UL0szmSjof+AMwT9IthA9+2wBfIyQ23A9MJYxJ3S3pGqADQYCqpLRHHw8PIWT+LQfaAecSBC4pLN0l3QdMI/Safg18ANwc+bZB0kTgKkl/ifzoShgXWwlMzuAW7wQmSJpAEMD+wPnAxp1eVQ1m9qake4DLFP4veobwDdeYqM1Pa9Ou4zj5RMsSoEyp07/kzewmSfMIWXRXEMJa5cBC4AbgJjPbJOkkgkg8RBijuZUwVnRrrLkFwFFRvc8Dm6Ky75jZfxKmLwO+QUh5LyYIyQWp1PTIt6slfQxcQMiKKwVmAZdGH9TWREpUzySMrf2b0Bt6JINrq+Mswn3/GLiEcH8nAzOppQA6jpM/5OsyH2pO39JIGkzoiQw3s4eb2J16QWGGj38Bp5nZXY1svvn88h2n5ZF192jbpwsy/pttXTCwxXS/fGykEZF0FHAg4VuuCuB/CL3Dt4AHm9A1x3GaAQV5uh5Vft5107EJOAa4C/g78AvCmNzhZpYLk+E6jpPTFGSx1R5J7SVNkbQqmq5tvqTv16nRuvjTnEJ/Tr3jv3zHaTqyDs19am9m/DdboP1qHfqTNJswZdsvgPcIY/UjgGPN7Mnatltrf1yo8hr/5TtO01ELoVqYhVD1r5VQSTqGsNTQ8Wb2UFQmwsQNu5vZvrVpty546M9xHKeZ0EjfUQ0nZCFXZjhH87H+P6CfpP51u4vs8WQKx3GcZkI2UyNluIJCujk/9wcWmlny287X4uczdqQecKHKb1pM+qrj5Ad9s/mbre1qELsTJrtOsi52vlFxoXIcx2mB1HGFhJ2NhTX62LaPUTmO4zhx1pK+15RaAWJdmnMNiguV4ziOE+dNYN9okds4A6L9G43sjwuV4ziOU4WHgE7AsYnyHwKLzaxREynAx6gcx3GcqjxJmFP1Nkm7Ez74PYMwEfiwpnDIP/h1HMdxqiBpV8LySScSelcLgQlNNRm4C5XjOI6T0/gYleM4jpPTuFA5juM4OY0LleM4jpPTuFA5juM4OY0LleM4jpPTuFA5juM4OY0LleM4jpPTuFA5juM4Oc3/B0KBlWIOEGtWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " This demo was completed in: 96.60s\n"
     ]
    }
   ],
   "source": [
    "sns.set_context(\"talk\", rc={\"font\":\"Helvetica\", \"font.size\":12})\n",
    "label = [antibiotics[i] for i in ab_order]\n",
    "cm = confusion_matrix(y_2019.flatten(), preds2019.flatten(), labels=ab_order)\n",
    "plt.figure(figsize=(5, 4))\n",
    "cm = 100 * cm / cm.sum(axis=1)[:,np.newaxis]\n",
    "ax = sns.heatmap(cm, annot=True, cmap='YlGnBu', fmt='0.0f',\n",
    "                 xticklabels=label, yticklabels=label)\n",
    "ax.xaxis.tick_top()\n",
    "plt.xticks(rotation=90) \n",
    "plt.show()\n",
    "\n",
    "print('\\n This demo was completed in: {:0.2f}s'.format(time()-t00))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}