Notebooks/linearprobes-PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c7b64b44",
   "metadata": {},
   "source": [
    "AI/ML directed high-throughput virtual screening of probes\n",
    "\n",
    "Siva Dasetty, Max Topel | Ferguson Lab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16e3398",
   "metadata": {},
   "source": [
    "Standard imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "153ee8bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import os.path\n",
    "import sys\n",
    "import importlib\n",
    "import yaml\n",
    "\n",
    "import subprocess\n",
    "\n",
    "from time import time\n",
    "from datetime import datetime\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "\n",
    "import rdkit\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "from rdkit import DataStructs\n",
    "from rdkit.Chem.Draw import IPythonConsole\n",
    "from rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions\n",
    "from rdkit.Chem import Draw\n",
    "from rdkit.Chem.Draw import rdMolDraw2D\n",
    "from rdkit.Chem import Descriptors\n",
    "from rdkit.Chem.Fragments import *\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "from matplotlib import rcParams\n",
    "import matplotlib as mpl\n",
    "from matplotlib  import cm\n",
    "from PIL import Image\n",
    "from matplotlib import colors\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.offsetbox import (TextArea, DrawingArea, OffsetImage,\n",
    "                                  AnnotationBbox)\n",
    "from matplotlib.cbook import get_sample_data\n",
    "\n",
    "import multiprocessing\n",
    "\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "sys.path.append(\"../Codes/\")\n",
    "\n",
    "\n",
    "import utils_activeLearning_Oliver as ut\n",
    "\n",
    "import utils_functions as utils_functions\n",
    "#import utils_figures as utils_figures\n",
    "\n",
    "importlib.reload(ut)\n",
    "importlib.reload(utils_functions)\n",
    "#importlib.reload(utils_figures)\n",
    "\n",
    "import selfies as sf\n",
    "\n",
    "from oliver_vae.VAE import VAE_encode, VAE_decode\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "\n",
    "import copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e68613e0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ea8432dd",
   "metadata": {},
   "source": [
    "Plot settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7d3700f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({'font.size': 20})\n",
    "\n",
    "# time stamp\n",
    "now = datetime.now()\n",
    "now = str(now)[:10]\n",
    "#figures_dir = './Figures/{}'.format(str(now)).replace(':','-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c52750a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c9191f68",
   "metadata": {},
   "source": [
    "Read all search molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e3b4a748",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read all potential linear probes\n",
    "\n",
    "all_potential_linear_probes = pd.read_csv('../Data/VAE_data/SMILES_linearProbes_All.csv')\n",
    "\n",
    "# convert to list\n",
    "all_potential_linear_probesList = all_potential_linear_probes.smiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6844121d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "dd235523",
   "metadata": {},
   "source": [
    "## get molecular weight and other properties using smiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7f255fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "molWts = []\n",
    "ccount = []\n",
    "fcount = []\n",
    "clcount = []\n",
    "brcount = []\n",
    "for i in range(all_potential_linear_probes.smiles.values.shape[0]):\n",
    "    molWts.append(Descriptors.ExactMolWt(Chem.MolFromSmiles(all_potential_linear_probes.smiles.values[i])))\n",
    "    ccount.append(all_potential_linear_probes.smiles.values[i].count('C') - all_potential_linear_probes.smiles.values[i].count('Cl'))\n",
    "    fcount.append(all_potential_linear_probes.smiles.values[i].count('F'))\n",
    "    clcount.append(all_potential_linear_probes.smiles.values[i].count('Cl'))\n",
    "    brcount.append(all_potential_linear_probes.smiles.values[i].count('Br'))\n",
    "    \n",
    "##smiles_latentspace\n",
    "molWtsArray = np.array(molWts)\n",
    "all_potential_linear_probes['molwt'] = molWtsArray\n",
    "all_potential_linear_probes['ccount'] = np.array(ccount)\n",
    "all_potential_linear_probes['fcount'] = np.array(fcount)\n",
    "all_potential_linear_probes['clcount'] = np.array(clcount)\n",
    "all_potential_linear_probes['brcount'] = np.array(brcount)\n",
    "\n",
    "\n",
    "# save updated dataframe [if running first time]\n",
    "#all_potential_linear_probes.to_csv('../Data/VAE_data/SMILES_linearProbes_All_Properties.csv')\n",
    "\n",
    "# read updated dataframe\n",
    "\n",
    "all_potential_linear_probes = pd.read_csv(\"../Data/VAE_data/SMILES_linearProbes_All_Properties.csv\")\n",
    "\n",
    "#all_potential_linear_probes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bceaae23",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "7f0890d5",
   "metadata": {},
   "source": [
    "Get SELFIE reprsentation of all molecules for passing to trained VAE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "33a41e78",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Representation: SMILES\n",
      "--> Translating SMILES to SELFIES...\n",
      "Finished translating SMILES to SELFIES.\n",
      "--> Creating one-hot encoding...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████| 3850/3850 [00:01<00:00, 2126.66it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finished creating one-hot encoding.\n",
      "Length of largest molecule:  130\n",
      "Encoding alphabets:  [')', '+', '[', 'C', ']', 'B', 'N', 'l', 'P', 'r', 'F', '(', ' ']\n",
      "selfies alphabet:  ['[epsilon]', '[C]', '[=C]', '[C@Hexpl]', '[N]', '[Ring1]', '[Branch1_1]', '[C@@Hexpl]', '[Branch1_3]', '[=O]', '[O]', '[Branch2_3]', '[Branch2_1]', '[=N]', '[c]', '[F]', '[Ring2]', '[Branch1_2]', '[o]', '[n]', '[S]', '[Branch2_2]', '[#N]', '[#C]', '[=S]', '[s]', '[-c]', '[Cl]', '[=c]', '[/C]', '[/c]', '[Br]', '[-n]', '[nHexpl]', '[=N+expl]', '[O-expl]', '[N+expl]', '[C@expl]', '[C@@expl]', '[/N]', '[\\\\S]', '[\\\\c]', '[.]', '[n+expl]', '[/O]', '[Expl-Ring1]', '[\\\\C]', '[P]', '[/S]', '[\\\\N]', '[\\\\O]', '[I]', '[Siexpl]', '[S@@expl]', '[S@expl]']\n",
      "Length of largest selfies:  137\n"
     ]
    }
   ],
   "source": [
    "# Encode molecules using Oliver's encoder\n",
    "\n",
    "print('Representation: SMILES')\n",
    "\n",
    "\n",
    "# all search probes\n",
    "selfies_list, selfies_alphabet, largest_selfies_len, encoding_list, encoding_alphabet, largest_molecule_len = ut.get_selfie_and_smiles_encodings_for_dataset(all_potential_linear_probes)\n",
    "\n",
    "# completed probes\n",
    "#selfies_list, selfies_alphabet, largest_selfies_len, encoding_list, encoding_alphabet, largest_molecule_len = ut.get_selfie_and_smiles_encodings_for_dataset(pfosfldataCopy)\n",
    "\n",
    "# unexplored linear probes\n",
    "#selfies_list, selfies_alphabet, largest_selfies_len, encoding_list, encoding_alphabet, largest_molecule_len = ut.get_selfie_and_smiles_encodings_for_dataset(unexplored_probes_list_df)\n",
    "\n",
    "# mobo probes\n",
    "#selfies_list, selfies_alphabet, largest_selfies_len, encoding_list, encoding_alphabet, largest_molecule_len = ut.get_selfie_and_smiles_encodings_for_dataset(mobo_candidates)\n",
    "\n",
    "#selfies_list, selfies_alphabet, largest_selfies_len, encoding_list, encoding_alphabet, largest_molecule_len = ut.get_selfie_and_smiles_encodings_for_dataset(mobo_candidates_batch2)\n",
    "\n",
    "print('--> Creating one-hot encoding...')\n",
    "# set largest selfies len according to Oliver data\n",
    "largest_selfies_len = 137\n",
    "onehotdata = ut.multiple_selfies_to_hot(selfies_list, largest_selfies_len, selfies_alphabet)\n",
    "        \n",
    "print('Finished creating one-hot encoding.')\n",
    "\n",
    "print(\"Length of largest molecule: \", largest_molecule_len)\n",
    "print(\"Encoding alphabets: \", encoding_alphabet)\n",
    "#print(\"Encoding list: \", encoding_list)\n",
    "\n",
    "#print(\"selfies list: \", selfies_list)\n",
    "print(\"selfies alphabet: \", selfies_alphabet)\n",
    "print(\"Length of largest selfies: \", largest_selfies_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a1c289f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "067265b1",
   "metadata": {},
   "source": [
    "Set device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f375f982",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "device(type='cuda')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "device"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3359364a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "f28693f9",
   "metadata": {},
   "source": [
    "Read trained VAE model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6d31ef3a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--> Acquiring data...\n",
      "-->cuda\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "VAE_decode(\n",
       "  (decode_RNN): GRU(100, 200, num_layers=3)\n",
       "  (decode_FC): Sequential(\n",
       "    (0): Linear(in_features=200, out_features=55, bias=True)\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load Oliver's model \n",
    "\n",
    "# Oliver VAE hyperparameters or settings\n",
    "\n",
    "if os.path.exists(\"OliverVAE-settings.yml\"):\n",
    "    settings = yaml.safe_load(open(\"OliverVAE-settings.yml\", \"r\"))\n",
    "else:\n",
    "    print(\"Expected a file settings.yml but didn't find it.\")\n",
    "\n",
    "print('--> Acquiring data...')\n",
    "\n",
    "len_max_molec = 137 #onehotdata.shape[1]\n",
    "len_alphabet = 55 #onehotdata.shape[2]\n",
    "len_max_molec1Hot = len_max_molec * len_alphabet\n",
    "\n",
    "encoder_parameter = settings['encoder']\n",
    "decoder_parameter = settings['decoder']\n",
    "training_parameters = settings['training']\n",
    "\n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")\n",
    "print('-->' + str(device))\n",
    "\n",
    "model_encode = VAE_encode(len_max_molec1Hot, **encoder_parameter).to(device)\n",
    "model_decode = VAE_decode(len_alphabet, **decoder_parameter).to(device)\n",
    "\n",
    "# load model\n",
    "model_dir='oliver_vae/Models/2021-08-05 12-14-37/'\n",
    "warm_path = model_dir #'./Models/2020-09-04 15-15-46'\n",
    "warm_epoch = 1999\n",
    "warm_model = os.path.join(warm_path, 'Epochs-' + str(warm_epoch))\n",
    "\n",
    "encoder_state = torch.load(os.path.join(warm_model, 'encode.tar'))\n",
    "decoder_state = torch.load(os.path.join(warm_model, 'decode.tar'))\n",
    "\n",
    "model_encode.load_state_dict(encoder_state['state_dict'])\n",
    "model_decode.load_state_dict(decoder_state['state_dict'])\n",
    "\n",
    "model_encode.to(device)\n",
    "model_decode.to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c57c7c00",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e5b393d3",
   "metadata": {},
   "source": [
    "Encode all search molecules using VAE and the SELFIE representation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9c6a7d59",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([3850, 7535])\n",
      "% correct:  0.0\n"
     ]
    }
   ],
   "source": [
    "# flatten one hot data\n",
    "inp_flat_one_hot =  torch.tensor(np.array(onehotdata), dtype=torch.float).clone().detach().to(device).flatten(start_dim=1)\n",
    "\n",
    "print(inp_flat_one_hot.shape)\n",
    "\n",
    "total_correct = 0\n",
    "all_correct_molecules = set()\n",
    "all_incorrect_molecules = set()\n",
    "\n",
    "encoded_list = [] #torch.zeros(inp_flat_one_hot.shape[0]).to(device)\n",
    "mean_encoded_list = []\n",
    "# encode and decode each probe \n",
    "for i in range(inp_flat_one_hot.shape[0]):\n",
    "\n",
    "    # encode and save to encoded list\n",
    "    fancy_latent_point, mus, log_vars = model_encode(inp_flat_one_hot[i].unsqueeze(0))\n",
    "    \n",
    "    # append encoded probe to encoded_list \n",
    "    #encoded_list.append(fancy_latent_point[0])\n",
    "    encoded_list.append(mus[0])\n",
    "    \n",
    "    mean_encoded_list.append(mus[0])\n",
    "    \n",
    "    \n",
    "    # decode the probe\n",
    "    fancy_latent_point = fancy_latent_point.unsqueeze(0)\n",
    "\n",
    "    # initialize hidden\n",
    "    hidden = model_decode.init_hidden()\n",
    "    #  runs over letters from molecules (len=size of largest molecule)\n",
    "    gathered_atoms = []\n",
    "    for _ in range(len_max_molec):\n",
    "        out_one_hot, hidden = model_decode(fancy_latent_point, hidden)\n",
    "\n",
    "        out_one_hot = out_one_hot.flatten().detach()\n",
    "        soft = nn.Softmax(0)\n",
    "        out_one_hot = soft(out_one_hot)\n",
    "        out_index = out_one_hot.argmax(0)\n",
    "        gathered_atoms.append(out_index.data.cpu().tolist())\n",
    "    #print(\"gathered atoms\")\n",
    "    #print(gathered_atoms)\n",
    "    molecule_pre = ''\n",
    "    for i in gathered_atoms:\n",
    "        molecule_pre += selfies_alphabet[i]\n",
    "    molecule = molecule_pre.replace(' ', '')\n",
    "    #print(\"molecule\", sf.decoder(molecule))\n",
    "    \n",
    "    molecule = sf.decoder(molecule)\n",
    "    if ut.is_correct_smiles(molecule):\n",
    "        #print(\"correct\")\n",
    "        total_correct += 1\n",
    "        all_correct_molecules.add(molecule)         \n",
    "    else:\n",
    "        all_incorrect_molecules.add(molecule)         \n",
    "        #print(\"incorrect molecule\", molecule)\n",
    "\n",
    "# Convert encoded list to Torch tensor\n",
    "encoded_torch = torch.stack(encoded_list)\n",
    "mean_encoded_torch = torch.stack(mean_encoded_list)\n",
    "\n",
    "encoded_array = encoded_torch.cpu().detach().numpy()\n",
    "mean_encoded_array = mean_encoded_torch.cpu().detach().numpy()    \n",
    "    \n",
    "print(\"% correct: \", 100*(total_correct)/1000)\n",
    "\n",
    "#onehotdata.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d2d358a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "3b4db147",
   "metadata": {},
   "source": [
    "Save/read encoded data of all search molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9011bf83",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('../Data/VAE_data/LatentVecs_AllSearchLinearProbes.npy', encoded_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bde9f255",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_search_linear_probes_array = np.load('../Data/VAE_data/LatentVecs_AllSearchLinearProbes.npy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f006a48",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ebbafa8c",
   "metadata": {},
   "source": [
    "Read latent space points of ZINC training set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "73c602b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Oliver's training data in latent space to perform latent space quality checks\n",
    "\n",
    "latentspace_extensive = np.load('../Data/VAE_data/LatentVecs_Extensive.npy')\n",
    "latentspace_primary = np.load('../Data/VAE_data/LatentVecs_Primary.npy')\n",
    "latentspace_ZINC = np.load('../Data/VAE_data/LatentVecs_ZINC.npy')\n",
    "\n",
    "# total latent space data\n",
    "latentspace_training = np.concatenate((latentspace_extensive, latentspace_primary, latentspace_ZINC), axis=0)\n",
    "\n",
    "# Load the corresponding SMILES labels\n",
    "\n",
    "smiles_extensive = pd.read_csv(\"../Data/VAE_data/Ext_Lib.csv\")\n",
    "smiles_primary = pd.read_csv(\"../Data/VAE_data/PS_Lib.csv\")\n",
    "smiles_ZINC = pd.read_csv(\"../Data/VAE_data/ZINC_Lib.csv\")\n",
    "\n",
    "#smiles_latentspace\n",
    "latentspace_smiles = pd.concat([smiles_extensive, smiles_primary, smiles_ZINC])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "589e2d07",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9a47847f",
   "metadata": {},
   "source": [
    "Space for calculating properties of molecules in ZINC set, if required"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f8b50b5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98fa52ad",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "2ea55340",
   "metadata": {},
   "source": [
    "Apply PCA over search set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "fffbd131",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[89.58982  91.56334  92.90315  94.01864  94.92091  95.50037  96.03267\n",
      " 96.5355   96.86376  97.1596   97.41349  97.65239  97.86564  98.05541\n",
      " 98.21808  98.37208  98.51624  98.64283  98.76724  98.87211  98.97498\n",
      " 99.06517  99.149445 99.222305 99.291756 99.358116 99.41777  99.47132\n",
      " 99.51957  99.56569  99.606476 99.642334 99.67528  99.70566  99.73209\n",
      " 99.75523  99.77686  99.79739  99.814415 99.83     99.84519  99.85845\n",
      " 99.87122  99.883316 99.89443  99.90424  99.913704 99.922714 99.93131\n",
      " 99.93866  99.945496 99.95228  99.95833  99.96304  99.96756  99.9719\n",
      " 99.97585  99.97932  99.9827   99.985725 99.98837  99.990875 99.993\n",
      " 99.994354 99.99535  99.996216 99.997025 99.99767  99.998276 99.998825\n",
      " 99.99927  99.99952  99.99974  99.999855 99.999954 99.99997  99.99997\n",
      " 99.99997  99.99997  99.99997  99.99997  99.99997  99.99997  99.99997\n",
      " 99.99997  99.99997  99.99997  99.99997  99.99997  99.99997  99.99997\n",
      " 99.99997  99.99997  99.99997  99.99997  99.99997  99.99997  99.99997\n",
      " 99.99997  99.99997 ]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtkAAAEQCAYAAABlQmh/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABgPElEQVR4nO3dd5xU1f3/8dd7d6UtAnaNJWA0oiZIUVGxYMVVFGOJmmiisfwSS2KJsXxVwBRbTFFjjMaIRsWosUTiigVRFEFBil0RURElUmWXPnx+f5w7MMzOlrs7s7Mz83k+HvOYnXPP3Htmy72fPfdzzpGZ4ZxzzjnnnMuesnw3wDnnnHPOuWLjQbZzzjnnnHNZ5kG2c84555xzWeZBtnPOOeecc1nmQbZzzjnnnHNZ5kG2c84555xzWeZBtnPOOeecc1nmQbZzzjnnnHNZFjvIlnSUpAclTZM0I6V8Z0m/krR1dpvonHPOOedcYaloakVJAkYAp0RFy4COKVUWAr8DBFyfpfY555xzzjlXcOL0ZJ8DnArcDWwM/D51o5l9CbwCHJm11jnnnHPOOVeA4gTZZwDTgLPMbDFgGep8CPTIRsOcc84555wrVHGC7J2AF8wsU3Cd9D9gs5Y1yTnnnHPOucIWJ8heDXRopM7WQE3zm+Occ84551zhixNkvwMMjAZA1iGpA3AQMCUbDXPOOeecc65QxQmy/wn0BP4oab33SSoH/gB8gzADiXPOOeeccyVLDadYp1QMgfR/gcOAL4AlwI7AY8BehAD7CTP7Xm6a6pxzzjnnXGFocpANIKkCuBI4F9gkZdMi4Bbg12a2OpsNdM45V1gkbQR0NrPP8t0W55zLl1hB9to3hbzsbxMC7cXAe2aWyHLbnHPOFQhJnYHhwA8Js0yZmVVE2/oDQ4ErzeyN/LXSOedaT7OC7LZu0003te7du+e7Gc65AjB58uR5ZuZTj7aApK7Ay8CuwFSgHbCzmZVH2zsBc4HbzeySfLWzqXJxDXn//fcB2GmnnbK6X+dcfjV0DYmzrPq3gAHAf81sfobtmwJHAC+b2czmNjYbunfvzqRJk/LZBOdcgZD0Sb7bUAT+jxBgn2Zm90oaClyd3GhmSyW9CBycrwbGkYtryMCBAwEYO3ZsVvfrnMuvhq4hcWYXuQy4Cfi6nu2LCUutt/leCuecc1l1LDDazO5toM4nhLUUnHOuJMQJsgcCz5nZqkwbo/JnCXNlO+ecKx3bANMbqVMDdI2zU0k7SPqbpGmSEpLGZqgjSVdI+kzSMkkvSeqdod4ukp6XtFTSHEnXRLNmOedcTsQJsrcGZjVS51PCVH7OOedKxxJg80bq9ADmxdzvroQ0xA+iRyaXAVcB1wNHEYL55yRtmawQzXbyHGDAEOAa4GLCQE3nnMuJJudkAyuBLo3U2ZBwEisIiUSC6upqpkyZQp8+faiqqqK83Ds2nHMupteBwZI2NLMl6RslbUUIlkfF3O+TZvZEtI9HgE3T9tuBEGRfa2a3RmWvEjqEziNMOQvwU6AjcKyZfQ08K6kLMEzSDVGZc21OapzSq1cvAKZOnUoikaC8vJzevXuvVyYJM8u4LdtlrXms1j5+v379shITxgmy3wKOlHRBppQRSe2AwYTl19u8RCLBoEGDmDhxIrW1tVRWVtK/f39Gjx7tgbZzzsXzZ6AaeErS2akbJO0M3Al0AG6Os1MzW9NIlX0InT8PpbynVtKTQBXrguwqQs54ajD9IKH3+wDgyTjtcsUvGdxOnjw5b4FfIpHgkksuYc6cOSxfvpzy8nLMjDVr1v1ZZCpraFu2y1rzWK15/M6dO2clJowTZN8H3AY8JOlnZvZlckN0W+52YFvghma3phVVV1czceJEampqAKipqWHixIlUV1czePDgPLfOOecKh5mNljQMGEbokFkFIGkesBEg4FIzG5/lQ/cEEsCHaeXvAiem1RuT1uZPJS2Ntq0XZE+ePJmwHMT6inHK22LWkkD5iSeeYMaMGdTW1q7dXzYCv7KyMsxsvd+lppYlEnWXI8lUFrd+S8pa81itefxsxYRxguw7gOMI+WyHSpoOfE7I1e4FdCLkvN3e7Na0oilTpqz3xwNQW1vL1KlTPch2zrmYzOwaSeOAnwN7ERYrM+Ap4I9mNqah9zfTRkBNhsXQFgKdJLUzs5VRvUUZ3r8w2uYKSEMpFMme5GwFyunHbUpZQ9syHaOpZa51ZSMmbHKQbWZrJB1BGCjyM8JJNGkR8CdgeBNu77UJffr0obKycm1PNkBlZeXa/2idc87FY2YvAC+09mEzlCnDtvrq1Snv16+fr7WQZ/X1QjclhaKx/TalzLlsxIRxerKT0/RdIelKwi22boQA+71CCa6Tqqqq6N+/P2PGjMHM1ubfVFVV5btpzjnnmmYhsKGk8rTe7G7A0pTxQwujsnRdydzD7VpRekDdUC/0mjVrGk2hKFaSkFQn/QQy93xn2pbtstY8VmsePzlOr6UxYawgOykKqAtigGN9ysvLGT16ND169KBLly5cd911PruIc841g6SDgR8CV5rZnAzbvwH8BrjXzMZm8dDvAeXADsD7KeU9o22p9XqmtWlboDKtnsuxpgbU9b03n1o78JNEIpGgQ4cOfOMb3+DGG2+kvLycadOmsXr1aioqKthtt90A1paVlZWxZs2ajNuyXdaax2rt4/ft2zcrMaGKcSDH7rvvbk291bfbbrux/fbb89hjj+W4Vc65tkjSZDPbPd/tKGSSHgd6mlnPBuq8C7xjZsc18xiPAJua2cCUsg7AXOBGM/tNVNaJMIXfHWZ2ZVR2OWE14m8mpxiU9EvCfNlbpk/hF+ca0lSlvKx6amD92GOPNSmgzpWmBrmdOnVixx135JhjjokVyGU78HvzzTfp3bu3dwK2YQ1dQ2L1ZEvaEfgFsCdhsEimn7iZ2bditzJPysvLWb16db6b4ZxzhawvYeB7Q14GDouz0yhgPiJ6uTXQRdLx0eunzGyppOuAqyQtJPRKX0RYaO2WlF3dThiQ+aik64HtCTOh/MHnyM6+1KB61apVPPDAA3zxxRcsX74868fKlEKR1NxAubGezCFDhjSpLG79uPt1bV+Tg2xJexNOoh2B1YTeg0zRad15j9qwioqKvN+Ccs65Arc5UCdNJM1cGl8VMtN+H04rS77uQeixvo4QVF9OmNFkEnComc1NvsHMFkYpLbcSputbBPyREGi7LGiN3uqmpFCk9iRnM1B2rjni9GRfC7QnrJz1DzMriu7fiooK78l2zrmWWUxYJ6Eh2wKxoi4zm0UjHTcWch5/Gz0aqvcOcFCc47vGJRIJRo0axUUXXbR2xo9saKgXOlMKhQfIri2KE2TvATxiZnfkqjH54OkizjnXYq8Bx0jaMnWhsqRo4OMxwCut3TCXfem91u+//36LgutMAXXcXmjn2qI4QfZK4NNcNSRfvCfbOeda7BbgSGCcpIsJS5ivkNQeOBy4CehMzGXVXduSjV7r9u3b07NnzyYH1M4VsjhB9nigT64aki8VFRWsXLky381wzrmCZWbPSPo1cBXwGGDRQMTkkuoCrjGzp/PYTNcM2eq1bt++PVtvvTV/+MMfGDx4sAfUriTECbKvAMZLOtXM/pmrBrU278l2zrmWM7Ohkl4Bzgf6ExZ/WQBMAG4xs2fz2DwXU0t6rb232rkgTpA9BBgDjJB0JjCZzCtlmZn9OgttaxWek+2cc9lhZs8Az+S7Ha5lEokEgwYN4pVXXokVXCdn/PDeaueCOEH2sJSv94semRhQMEG2T+HnnHPOrUsNefDBB5scYCd7rY899ljvrXYuTZwg+8BcNEDSScCvgG8TpoF6HrgsdWleSSLMgfozYFPgdeDnZja1pcf3dBHnnMseSZWEVJGMkZaZFd0A+kLXnNQQ77V2rnFNDrLN7MVsH1zS0cBI4C+EJW+3An4DjJK0u5kll3C6jDCg5hLWrej1nKTvZJouKg4Psp1zruUknQpcCuzcQDUj5krDLrfipIZ4r7Vz8eT7ZPcD4A0zOy9ZIOlr4AlgJ+BdSR0IQfa1ZnZrVOdVwkpf5wFXtqQBnpPtnHMtI+k04B9AAhgHfEbmFYFdGxEnNcR7rZ1rnnwH2RsQUkRSLYqek6t87QN0AR5KVjCzWklPAlW0MMj2nGznnGuxXwILgX3N7N18N8Y1LNl7PXHiRGpqauqt16FDB44//nhOPPFE77V2rhliBdmStiIEtYOArYF2GaqZmTV1v/8AHpf0I+BxYEtCusgL0RK4AD0JvSMfpr33XeDETDudPHkyIY27TsPqlHm6iHPOtdgOwD0eYLd9iUSC4cOHM27cuAbXiOjQoQMDBgxgxIgRHlw710xNDrIlbU1YOncL4G2gPfAJsALYPtrXVOr2TNfLzP4b3Wa8C7gnKh4PHJ1SbSOgxszSu5sXAp0ktTOzZq8m4+kizjnXYguA5q+r7VpFsge7oQDbU0Ocy56yGHWvJvQ0H25mu0Vld5tZT0KQPRroCBzb1B1KOhC4HfgzYfaSk4CNgcckpf5l1+2CXpdOUmdbv379MLM6j0y8J9s551psFDBQmW4hurxLzh7y4x//mJdffjljgN2hQwdOOeUUHn74YT744AOGDBniAbZzLRQnXWQQ8LSZPZe+wcxmSzoBeAsYDvy8ifu8CfiPmV2aLJA0lTCDyBDgUUKP9YaSytN6s7sBS81sVYzPUIfnZDvnXItdDrwC3C7pYjOrP9HXtaqm5F+3a9fOU0Ocy4E4PdlbEtJEkhKEnmsAopPqs4TguKl6ElJM1jKz94FlwLeiovcI863ukOG978U4Vkbek+2ccy32MLAUOBP4QtIbksZkeDyf53aWlNT864YC7EsvvZTRo0d7gO1clsXpyf6a9Qc6LiQMfky1GNgsxj4/AfqmFkjamRC8z4qKxkfHPoEwKBJJnYCjgDtiHCsjz8l2zpUSSWXAucAPCXNaVyYHq0vqA5wF/MnMPoix24EpX1cCveuplzlvz2VdU/KvO3fuTP/+/Rk6dKgH2M7lQJwg+xNg25TX04CDJHUys6XRifswYHaMfd4O/FHSHKCaMKjyakKA/RSAmS2XdB1wlaSFrFuMpgy4JcaxMvJ0EedcqZDUjnCuHUgYrLgE6JxS5WPgJ8BXwNCm7tfM4twVdTnW2AwiPjWfc60jTpD9PHC2pA2iPOh7gHuB8ZKeBfYFdgV+F2OfNwMrCcul/5QwR/bLwOVmVptS7zpCUH05sAkwCTjUzObGOFZGni7inCshlxAGmQ8j3Bm8mrCaLgBmtkjSS4QxOE0Osl3b0VgPtudfO9d64gTZdxFSRDYFvjCz+yT1A84HekV1HgR+29QdWpjy46/Ro7F6v42z76aqqKhgzZo1rFmzhrIy74xxzhW1HwKvmNk1AJIypW98TEjHcwWourqaiRMn1htgX3rppZ4e4lwraXKQbWYfAtenlV0o6XeEKfxmZaNnubUlTzSJRMKDbOdcsesB/LeROgsIU6k2i6RtCON12mfabmYvNXffrn7JZdL/8Ic/ZBzk2K5dO/bbbz8PsJ1rRS1eVt3MviLk7xWkiorwLUgkEmywwQZ5bo1zzuXUMsL0pw3ZjpC6F4ukw4A/EmZ+akjWIzxJxwDXADsBc4BbzOwPaXW2IqQzHgZ0Jawi/Hszuz/b7WltjU3T5z3YzuVHyXfdJoNsz8t2zpWAqcBh0QDIOiR1JeRjvxZnp5L6Exak6QbcSlgs7CXgTsJgdQFPEgLhrJI0gLCmwmuENJd/ANdLuiClThnwH+AA4FeEqWYnAPdJ+l6229Takiki6QG2JDp37uw92M7lSb092ZL+QZhu6Qozmxu9bgozszOy0rpW4EG2c66E3AncD9wvab3ztKRuwN3ARoSZn+K4grCs+h5mNkfS+cALZnZNtArkMOBi4P9a1vyMrgZeNrMzo9fPSNoIuFrSbWa2Evg2sDtwtJk9GdV7Pvrn4CTgsRy0q1UkEgkefPDBjD3YBx10EBdccIHPIOJcnjSULnIaIci+HpgbvW4KAwomyE6eeDzIds4VOzMbKekQ4HTgaMJgdiRNIswO1R74i5k9FXPXexNW752TUlYWHdOAoZKOIKwIfHzLPkUdvYG/pJU9A1wYtetFIJkLuDit3iJCL3tBSqaJvPjii3W2de7cmQsuuIDBgwfnoWXOOWg4yO4RPX+e9rqopOZkO+dcsTOzMySNA35BmBlKhEXB3gb+YGZ3N2O3XYFPU16vJCxKk+oV4AfN2HdjOkTHS7Uiet6ZEGS/BUwErpF0FqHj6FhgAFCVgza1iurqasaPH1+nk6hDhw7079+fqqqC/WjOFYV6g2wz+6Sh18XC00Wcc6XGzEYAIyR1JKSHLE5bmyCu/0X7SX39rbQ6GxBW8822GcAeaWV7Rs8bQ+hNl1QFPAEkV7JcBZxuZmPSdzh58mRClsv6Qqd8/iVnErn22mtZtmxZne3HH3+8z4PtXBvQ5NlFJM0Eqs3s3By2p9V5uohzrlRIKjOzNcnXZraMMONIS33A+kH1BKBK0rfN7ANJWwLHEWb0yLbbgb9GPdSPEALsi6NtCVg78PGfhMXMTiT8E3AEcJek+Wb2dA7alRPJFJEJEyZQW1v3/6LOnTtz4okneoDtXBsQZ3aRzaibz1bwPF3EOVdCZku6XtKuWd7v08ABkpLza/+Z0Gs9RdLrhBlGNgP+lOXjQphN5HbComYLCDONJGcxSa7dMBg4EjjGzB4ys7Fm9ivCgMcb0nfYr18/zKzOoy1IziSSHmAnZxLxNBHn2o44Qfbb1L39V/A8XcQ5V0I6EpZWny7pdUnnStokC/v9G7A/IQUDM3sFOIGweuR3gC+An5nZvVk41nrMLGFm5xGC+F7AFoSedFKeewJLo0XVUk2hgK5rjc0kMnLkSEaPHu292M61EXGC7JuBoyT1arRmAfEg2zlXQrYgTFn3NGFWjpuBzyX9W9LRkpoVnZnZ12Y20cyWpJQ9ZmbfMbOOZrazmd2RjQ/QQBsWmtmbZlYDnAOMN7P3os2fAJ0k7ZT2tn7ArFy2K1uSaSKPPPJInW2pM4l4gO1c2xFnxcfZwHPAK5L+BrwOfEmYsm89hbRsrudkO+dKRTRn9EPAQ5K2AE4Bfgx8DzgGmC/pfuBeM5uSt4bGIGkvYF/CQjtdgJMJC+rsm1LtKcLsJ49LuoawSvGRwPeBghhnlEwTWbFixXrlPpOIc21XnCB7LCGgFnARGYLrFAXzr7TnZDvnSpGZzQVuAm6S1IcQbJ9EmNrvfOJdH/JpFWEw4zBgDTAOGGBmbyYrmNkSSQcD1xI+cxfgI+CnQE572LNlypQpGQc6+kwizrVdcU6i19BwYF2QPF3EOVfqzGyKpBrC/NIX0Mi1IZptyoBDzOzj6HUTD2VZzYE2s8nUncIvU70ZhDzxgpKcru/tt9+uM/jSZxJxrm1rcpBtZsNy2I688SDbOVeqJHUl9F7/GOgfFS8BHm7krWWs3+mS/rreQ8ZtYylL5mFPnDhx7WDHsrIyzIzKykpPE3GujSuU24E54znZzrlSEs0ZfTghsD6KsJS6Ac8DI4BHzWx5Q/sws+4NvXbZkczDTp1NpF27dhx//PGceOKJVFVVeS+2c21YyQfZnpPtnCsVkn4P/BDYnNCr/AFwD/BPM5vdgv1uB6w0sy+z0lAHZM7DXrFiBTvttBODBw/OU6ucc00VK8hWWGf2eMLI7a0JPSDpzMwOzkLbWoWnizjnSshFhEXF7gTuMbNXs7TfjwnB+k+ytD8H9OnTh3bt2q03o0hlZSW9e/fOX6Occ00WZ1n19oRpkAYSekCSM40kWUp5wfB0EedcCfkB8JiZrWi0ZjyLgHlZ3mdJSyQS1NbWsnr1akL/Fp6H7VyBibMYzaXAgcBvCCtriTBl0jcIJ+7PgAeBdtltYm55uohzrlSY2YM5CLAhrKzYJwf7LUnJAY+nnnoqiUSCDTbYgB49enDffff5io7OFZA4QfYJwBtmNtTM5icLzexLM3sQOAgYTJj+qWB4uohzzrXYMGA/SWfmuyHFoLq6mgkTJrBq1SoAVq5cyf/+9z/Ky8s9wHaugMTJyf4WIY8vyYAN1r4wmynpv8BphMn+C4IH2c65YiVpDWGBll3M7IPodVNS+szM4lwfqggLlv1N0s+A18i8IrCZ2a9j7LckZRrwWFtby9SpU33Ao3MFJM5JdBWQOq3TEkLaSKpPgKNb2qjW5DnZzrki9hIh0F2a9jrbhqV83Yf6U0cM8CC7HsmFZ6ZMqbuivQ94dK7wxAmyZxNmFEn6ANg7rU4fYEFLG9WaPCfbOVeszGxgQ6+z6MAc7bdk+MIzzhWfOEH2K8AhKa8fB34j6S7gUcKsI4cAD2Srca3B00Wcc65lzOzFfLeh0PnCM84VnzgDHx8A3pfUPXr9J+B14HTgP8DFwEfAZXEaIKlC0mWSPpS0QtJsSX9MqyNJV0j6TNIySS9J6h3nOPXxINs5VyokzZT080bqnCtpZmu1yQWNLTzjAbZzhafJPdlmNpYwsCX5eqmkAcAQYAdgFvCkmS3N9P4G3A0cDAwH3gO2BXZJq3MZcBVwSVTnIuA5Sd9p6QpjnpPtnCsh3YFujdTpBnwz1w1x6+vTpw+dOnVaL9D2PGznCluLllU3s9XAv5v7fkmHAycBu5nZO/XU6UAIsq81s1ujslcJQf15wJXNPT54TrZzzqXpDKyM+yZJWxHOx8kVgTOtmRB31pKSUVVVxVZbbcWMGTOQ5HnYzhWBOCs+/gx40MwWZvH4PwHG1BdgR/YBugAPJQvMrFbSk4Rpo7ISZHtPtnOuGEnaLq2oW4YygHJgO+B4IFa6iKStCdP2bQG8DbQnzDa1AtiecK2ZSljS3aVJJBI8+eSTzJ07l29+85v85Cc/oW/fvp6H7VyBi5OT/RdgjqSHJB0pKc5769Mf+EDSrZK+lrRU0qOSvpFSpyeQAD5Me++70bY6Jk+ejKQ6j0w8XcQ5V+RmAR9HD4BfpLxOfcwAxlB3TYSmuBrYEjjczHaLyu42s56EIHs00BE4ttmfokglZxU56aSTWLJkCXPnzuWll17yANu5IhDntt0VwI8IvRzHAV9Jug+418ymN/P4WxIWr5lGSBvZELgBeEzSXmZmwEZAjZml53MsBDpJamdmsW9tJnm6iHOuyN1LmJ9ahHP4dEKvcroEMB943syeiXmMQcDTZvZc+gYzmy3pBOAtwtibBgdelprkrCIrVoTV7pcvX87EiROprq72hWecK3BxBj5eB1wnaXdCYHwiYQDihZKmASOAkWb2VYzjK3oMSS7VLukL4EXCMu3PJw9fz3szbuvXrx+TJk1qUgM8XcQ5V8zM7LTk15J+BDxmZtdk+TBbkpLSRwjYO6a0oUbSs4SB8h5kp/DVHZ0rXrFTPsxskpmdB3yD0Ks9CtiVMKXfbEmPx9jdQuDNZIAdeZkw6GaXlDobSkq/b9YNWGpmq+J+hlQeZDvnSoWZleUgwAb4mvUHOi5k/cXLIORjp68SXPL69OlDu3brjxH1WUWcKw7Nzqs2s1Vm9qiZDSEE3FdHm46KsZt36ykXsCb6+j3CgJwd0ur0jLa1iOdkO+dci31CmH41aRpwkKROANEYnsMIKwe7FFVVVXTu3Hnt2KHOnTv7rCLOFYkWTaWkMJrwUODHhNuAGxBuEzbVKGC4pE3NbF5Utn+0n2nR6/GEXpITgN9Ex+1ECObvaEn7o31RVlbmOdnOuZIhaQ/WTbXXPkMVM7MzYuzyeeBsSRtEdxfvIeSCj4/SRPYl3PH8XctaXnxWrFhBTU0NgwcPZs8996R3794+6NG5ItGsIFvSzoTA+hRgK0LP84eEk+q9MXZ1ByE/70lJvyMMfLweeM7MXgYws+WSrgOukrSQdYvRlAG3NKf96SoqKrwn2zlX9KKOkRGEc7dYNyAyyVLK4wTZdxFSRDYFvjCz+yT1A84HekV1HgR+25L210fSMcA1wE7AHOAWM/tDhnrfBa4F9iNcQ94FfmZmk3PRrqYYM2YMK1as4Oc//zmHHHJIvprhnMuBOPNkbwycTAiu+xFOxF8TTq4jzGx83IOb2deSDgJuJpyAVwJPABemVb2OcEK8HNgEmAQcamZz4x4zEw+ynXMl4jzgVEJnyM2Ec+mfCIMWBxIW/nqKcK5tMjP7kNBBklp2YdR5sj0wK1vn63TRysOPAv8AfkmYGvZ6SWvM7E8p9XoD4wjXmBOj4j1IGaDZmhKJBNXV1QwfPpwOHTowYMCAfDTDOZdDcXqyv4jqG/AcoTfkMTNb3pIGmNkM4IhG6hihByQnvSDl5eUeZDvnSsGPgfeTM45E6wcsMrMJwARJo4EJwLPA3S09WDTbVJwZp5rjauBlMzszev2MpI2AqyXdljLF6+3Ak2Z2Ssp7n85x2zJKzo09ceJEampqKC8v56ijjmL06NGeJuJcEYkz8PFj4P+A7cxskJmNbGmA3VZUVFR4TrZzrhTsRFhwJtXazhYzm0IYK3NOnJ1K+pekqiwtUhZXb0LHT6pnCGss7A0gaRdCD3dWUgxbKjk3dk1NDRCC7uTc2M654tHkE6KZ9TSz68xsTi4blA+eLuKcKxFi/aXNa4GN0+p8SD2r6TbgBEJw/rmkG6Pc59bSgZBqmGpF9Lxz9Nw/et5I0jRJqyV9JClj3nmcVYObo6G5sZ1zxSMfvQ5tjqeLOOdKxOesP3/1TMIYm1Q7EoLvOPYG/kaYK/tiYKqkSZLOl7RpcxvbRDMIudWp9oyek/9AbBk93wvcT5gV62ng75IaTFfMhT59+lBZWblemc+N7Vzx8SAbTxdxzpWM11g/qK4G9pR0laRdJZ1LmI51QpydmtlEMzuHMNvU9wmDJ3sBfyb0bj8q6RhJLZo2th63A0MknSVpI0mDCIE+rJtSNnmt+7uZ3WBmL5jZucALZBjk2a9fP8ysziNbqqqq6Ndv3Y/B58Z2rjh5kI2nizjnSsa/gXJJPaLXNxAWkhkOTCfkLC8izDISm5mtNLNHzOwoQo/5xYRp8o6Jjp2LdMN/EALtvwILCDONJFe1TM5osiB6fiHtvWNYt7pwqykvL+cXv/gFAKeffjojR470QY/OFaFc9CoUHA+ynXOlwMweBx5Peb1AUh/gLOBbwCzgXjP7IgvH+gr4o6Q/EaZlvZYwBWtWmVkCOE/SVcA2hEH6yZzyZI98U1YXblXPPfcclZWV3H777XWWVXfOFQcPsvGcbOdc6TKzxcDvs71fSTuxbtGyrVm3aFlOmNlCwoI4SDoHGG9m70Wbx0fbDgZGp7ztYNatLtyqnn76aQ488EAPsJ0rYh5k4znZzjmXDdH81CcRgus9WH/RsnvM7JUcHHMvwrLtU4EuhEXTBkVlQEhjkXQNcIOkRcDrwHHA/sAB2W5TQxKJBHfddRczZ87k0EMPJZFIeJqIc0XKg2w8XcQ5V5wk7d/c95rZSzGOM5gQWA8mzDCSXLTsHuDRHK+psIqwguMwQurHOGCAmb2ZWsnM/hTN431+VPd94HgzG5fDtq3HzBg0aBAvv/wyAP/85z+ZMWOG52M7V6TqDbIlbdfcnZrZp819bz54kO2cK1JjCQFvc8SJ+v4TPX9ACKzvNbPPm3ncWMxsMnWn8Kuv7h+AP+S2RfVbsGABs2bNYsWKMI330qVL1y5CM3jw4Hw1yzmXIw31ZM+ieSdna2S/bY7nZDvnitQ1ND/IjuMOYES0PLurR01NTb2L0HiQ7VzxaSgYvpe6J+cehBy2xYT8ty8Jk/z3BroCLxFGdhcUz8l2zhUjMxvWSsf5aWscp9B17tyZDh06sGzZsrVlvgiNc8Wr3iDbzE5LfR2NFH8V+CMw3My+TtnWhTDP6o+As3PS0hzydBHnnHO5tvHGG7Ns2TI++ugjJFFZWemL0DhXxOKkdVwHvGlmF6dviALuCyX1i+odm6X2tYqKigpWrlyZ72Y455wrYpLo0aMHq1at4qyzzqJ3795UVVX5oEfnilScIHt/wqpaDXkZ+H/Nb05+lJeXe7qIc67oSRrTxKpmZgfntDElaM2aNbzyyiucddZZXHnllflujnMux+IE2e0J+dcN2SqqV1A8XcQ5VyIGNrLdCHNbt8ZgyZKzZMkSli1bxoEHHpjvpjjnWkFZjLpTgJOiJXjriFJFTgTeyEbDWpMH2c65UmBmZZkewEbAYYQB7f8izHXtsmzRokVI4oADWnX9G+dcnsTpyR4OPA1MkHQ/YSaRucAWhBWzfkAI2odnu5G55lP4OedKWbS0+nOSDgXeAi4Gbshvq/Ln/fffZ+DAgVnbn5kxadIkli1bRvv27TnmmGOQlLX9O+fapib3ZJvZc4TlcmuA0wjL5I6Knn8clZ9kZs9nv5m55VP4OeccmNkC4CngzHy3pViYGdOnT2fp0qWYGStWrGD69OmYeUaOc8Uu1qIxZvaIpGpgCNCXMDf2YkKKyBNmVtvQ+9sqTxdxzrm1vgYaXPFX0o+au3Mzu7e5720tO+20E2PHjs3KvkaNGsXJJ5+89rWZsXr1an75y1/ywQcfcMghh7Dbbrtl5VjOudbX0F2p2CszRoH0A9GjKHiQ7ZxzIKkjcCTwv0aqjmD9wZFNGSyZrNPmg+xsmjJlSsZVHqdMmcL111/P5Zdfzrbbbsu5557LKaecwuabb56nljrnsq3Zy59L2gjobGafZbE9eeE52c65UtBAD3QFsC1hbM0OwO8b2dXpGcqOBY4CXgTGsm5F4AMJU8D+B3gsdqMLXJ8+faisrKSmpmZtWWVlJX379uWggw5i1KhRzJw5kyuvvJIrrriCAQMGcN5553HkkUfSrp2PP3WukMUKsiV1Jgxs/CGwGaFXoiLa1h8YClxpZgU1w4jnZDvnSsQIMvc4J+93rgHuAxqcxNnM7lnvzdIRwOHAEDN7Mq36cElDgIdofK2FolNVVUX//v15/vkwXKlz585rV3lcvHgxY8eOpaamZu1S62PGjOH111/HzPjBD37A2WefTd++fX2gpHMFqMlBtqSuhMVmdiVM8zQP2DmlypvAfsDJFNg0fp4u4pwrEZl6oCEE1wuBSWb2ZTP2+3/AYxkCbADM7AlJjwNXEWapKhnl5eWMHj2a9u3b07FjR0aOHLl2lcfDDz8842rDS5YsAeCuu+7ivvvuY7PNNuOcc87hRz/6EVtu2dhyFc65tiLOPNn/RwiwTzOzvsDDqRvNbCnhNmHBrRLmQbZzrhSY2T31PP5pZqOaGWAD7AbMaKTODKBXM/df0BYsWEAikWCLLbZg8ODBa5dR32STTfj2t79d7/sSiQRLly7lk08+YejQoWy99dY88cQTrdVs51wLxQmyjwVGNzIy/BNg65Y1qfX5surOOdciKwmBdkN2A1a1QlvanHfeeQcIudjpTjzxxCbnXvft29dXi3SugMQJsrcBpjdSp4YwrV+zSNpaUo0ki/K/k+WSdIWkzyQtk/SSpN7NPU4678l2zpWa6Ly6laTtMj1i7u554AhJ5ykteTg6zvlAFfBcttpfSJJBdqdOnepsO/roo9lggw0afH/Hjh0ZMGAAL730El26dMlJG51z2Rdn4OMSoLG5hXoQcrWb60ZCoJ7+7/5lhFy+S4D3gIsIq5N9pwW3N9fyINs5VyoknUA4p34XKK+n2tpB7U10GWEWkT8DF0h6mXUrAu9LuDYsiOqVnLfffpvy8nLat29fZ9t3v/td2rdvX2eav3QPPvggHTt2zFUTnXM5EOck+jowWNKGZrYkfaOkrYAjCKtAxiZpP8Lo9N8Rgu1keQfCiflaM7s1KnsVmAWcRyOj4JvCp/BzzpUCSecCNwOrCQPZP4++bhEz+0jSXsBtwCHA9mlVngXONbOZLT1WIXrnnXcy9mJDWMji6KOPZsSIEfW+P5FIMGjQIMaNG1fvfpxzbU+cIPvPQDXwlKSzUzdI2hm4E+hAOIHHIqkcuAW4BliUtnkfoAth+icgLIgj6UnC7ccWB9k+hZ9zrkRcSFhoZh8z+zibOzazGcBhkrYG+rBuReApZvZ5No9VaN55552M+dhJxx13HI8++ihff/11xu0rV67knXfe4YgjjuCZZ57x+bOdKxBNzsk2s9HAMGAA8BZwOYCkedHrfYDLzWx8M9rxU0KA/pcM23oCCeDDtPJ3o211TJ48GUl1HvWpqKjAzFizZk0zmu6ccwVja+DhbAfYqczs82imkvuj55wH2JKOkTRd0gpJH0u6qJH6f4rG/jS26E6LzZ8/n7lz5zbYA33QQQexYsWKta87duxYJzVk+fLlvPbaa3z/+9/3TiHnCkScgY+Y2TWEKfr+Q5hTNUHI3XsKOMTMbmzg7RlJ2gT4NXCRmWUaeb4RUGNm6WeVhUAnSS3+l76iInToe8qIc67IfQbUTQzOIkk9JX1P0qm5PE7K8QYAjwKvEVac/AdwvaQL6qm/C/ATIHO3cZa9++67QOaZRZI6derEHnvsAYQA+5ZbbuHxxx+vE2gvW7aMZ599ljPPPBOzxlaxd87lW6wgG8DMXjCz75nZVmbWzsw2M7OjzGxMM9vwW2CimT3V0GEzlKm+bf369cPM6jzqk5yz1INs51yRGwFUSdow2zuW1FvSJOBt4JHoWMltB0haKumobB8XuBp42czONLNnzOzXhPTDq+vphLmZkP64MAdtWU8ikeDBBx8EQspHQ9ehk08+mfLycu6//37OOOMMDjvsMO655546gfbSpUt56KGH+NWvfpXTtjvnWi52kJ1NknYl9CgMl9RNUjcgeU+tq6SOhBPhhlHedqpuwNJ6er9jSfZk+y0451yRu54wiP25KPDNSrAt6dvAWGAn1o3fSfUSYXaR47NxvDS9qTs14DOEu6B7p7XzeMJKxdfloB3rSQ5WvOOOOwD48MMPmT59er3XmTPOOIO33nqL733ve2vLTjjhBP785z9nDLRvu+02rr322tx9AOdcizUryJZUGc1p3dL5VXcENgBeJQTTC1mXlz2b0BvxHmGaqR3S3tsz2tZini7inCsFUdrdXwjn0zHAIkmJDI+4J8OhQDtgTzO7iBDIpx7XCOf5PVr+KeroQFgMJ1UywXnnZEHUaXMTcJmZNThfXtxxPZlUV1czceJEVq0K/UBr1qzh66+/pro6/f+PoH379vTsWXeY0VlnncXVV19dJ6d76dKl/OY3v+Fvf/tbrHY551pPnNlFiHLsLiXlxJVBnPlVXybMrZrq8OgYRwAzCatIfg2cAPwmakcnQu7dHU1te0M8yHbOlQJJQwipHOXAx8AcsjCFH2GszqNm9m4DdT4FDs3CsdLNoG7wvmf0vHFK2eXAF8B9OWhDHVOmTKkz9/WaNWuYOnUqgwcPjrWvyy67jPnz53PbbbexdOnSteVLly7lwgsvZKONNuL73/9+VtrtnMueJgfZkk4jDChJAOMIA2hadHI2s3mEW4ypx+kefTnOzGqisuuAqyQtZN1iNGWEnu4WS+Zke7qIc67IDQOWAkea2ctZ3G83wt3HhpQReruz7Xbgr5LOIvwDsSdwcbQtASCpB/BL4CBrwojBfv36MWnSpBY1qk+fPlRWVlJTU7O2rKysjN69ezdrfzfccAMLFizgwQcfXC/QXrZsGaeddhpdu3Zl0KBBLWqzcy674vRk/5KQzrFvI70VuXAd4QR9ObAJMAk41MzmZmPn3pPtnCsROwH3ZjnAhjD3dnpKX7pdCZ0z2fYPYDfgr4S7m0sJd0NvIaw6CeEaUg28F439gXBNaR+9XtyU4DuOqqoqevXqxfjx49emm3Tp0oWqqqpm7U8Sd955JwsWLGD06NEsW7Zs7bZly5Zx7LHH8uyzz7LPPvtk6yM451ooTk72DsAjuQ6wzWyEmSnZix2VmZn91sy2MbOOZrafmU3J1jE9yHbOlYh51M1fzoYxwFGSdsq0UdIehJSS0dk+sJklzOw8YDOgF2Ep9wnR5uTzTsCxrBv7sxDYlrBq8ELC/OFZVV5ezrBhwwA47bTT2GWXXejVq9faO6fNUVZWxkMPPcRee+1Fhw4d1tu2dOlSBg0axPTp01vSbOdcFsUJshcAy3PVkHzyKfyccyXi38ChkjbI8n6vJaQPviTpZ8A3IMwgFb1+ElgC5GzxFzNbaGZvRh005wDjzSw5OP5Mwvif1MdcwkrCBwJf5aJNn3zyCQBDhw5lk002iT14MpMNNtiAUaNGscsuu9RZ+bGmpoYDDjiAGTNmtPg4zrmWi5MuMgoYKEnZvq2Wbz6Fn3OuRFwJ9AcelnSBmc3Kxk7N7H1JxwEjgVujYgHTo+dFwLFm9mk2jpdK0l7AvsBUoAtwMjAoKku2r06CtaTlwGdmNjbbbUqaOXMmFRUVbLPNNlndb6dOnRgzZgx77rknM2fOXK+DaPHixey7775MnjyZrbfOege9cy6GOD3ZlxNWCrtdUucctScvPF3EOVci3iSkSRwFfCRpvqSZGR4fxd2xmT0N9CAMTH+IMHf1o8AlwA4tWLCsMauAE4HHCQvgdAIGmFne8yZmzpxJ9+7dW5QiUp+uXbsybtw4ttxyS8rK1l3KzYx58+YxYMAA5s+fn/XjOueaLk5P9sOEASVnAj+Q9CGhdyKdmdnBWWhbq/Eg2zlXIsoIaR2pPcqZchialddgZosIi9H8uTnvb+YxJ9OM+bfNrHv2W7O+mTNnsv322+ds/5tvvjmvvPIK/fr1Y/78+WtXlEwkEsyZM4f999+fCRMmsOGGWV/g0znXBHGC7IEpX1cSVtnKpOBSSTwn2zlXClojsHTrzJw5kxNOOCGnx9huu+14+eWX6d+/P4sXL15bvmrVKj766CMOOeQQXnzxxToDJZ1zudfkINvM8roEey55TrZzzrWcpM2B3QlLmmfMkTCze1u1UXmyePFi5s+fn9Oe7KSddtqJMWPGsP/++6+3AM6KFSuYPn06Rx99NE899dTaa51zrnX4XxyeLuKccy0RzVZyO/Aj6h/rI8KdzpIIsj/++GOAVgmyAfr27ctTTz3F4Ycfvt4c2suXL2fMmDF89NFH7LRTxhkWnXM54kE2HmQ750qDpB81tW7MHudfA6cDHwH3k4UVgQvdzJkzgdYLsgH2339/HnroIb7//e+vDbQ7derEww8/7AG2c3lQb5Ataf/oy9fMbHnK60aZ2Ustblkr8mXVnXMlYgSNj5tpTo/zD4APgD5mtqyxyqUgH0E2wODBg7nzzjs588wzqaiooLq6mn333bfxNzrnsq6hnuyxhBPtzoSTZ/J1U2R/vqIc8p5s51yJOL2e8m6EGTpOIixY89+Y+90cuM0D7HVmzpzJxhtvTNeuXVv92D/84Q/p0KEDO+64I7169Wr14zvngoaC7GsIQfW8tNdFx4Ns51wpMLN7Gtou6W5CgH1zzF1/SlgIxkVyPX1fY4477ri8Hds5F9QbZJvZsIZeFxOfws8558DMnpf0NKFT5aAYbx0BnCupq5ktbqxysUskEkyfPp2NN96YUaNGUVVVle8mOefyoGin5YvDp/Bzzrm1PiBMwxfHdcDLwHOSDpRUsr3aiUSCww47jC+++IK3336bk08+mUGDBq1dKMY5Vzo8yMbTRZxzLsUuxE8NXAUcD/QjLKe+UFIiw6PoT7LV1dVMmDBh7euamhomTpzIggUL8tiqurp3746kjI8tt9wy381zrijEmsJPkggn0kHA1kD7DNV8WXXnnCsgksqAbYGzgCqgOuYuxlGkY3bimjJlynrzVAPU1tZSU1PDJptskqdWZda1a1cuuOCCOuWdO3du/cY4V4SaHGRLag88RVhePTnFk1KqWEp5QfGcbOdcKZC0hobP0QLmA5fE2a+ZDWxBs4pKnz59aNeuHStWrFhbVllZ2SYD127dujFs2LB8N8O5ohWnJ/tS4EDCogM3A18Bw4A7CIH3dcB44NSstrAVeE62c65EvETmIHsNsBB4DbjbzL5q1VYVkaqqKrbZZhs++ugjJFFZWUn//v1ZtWpVvpvmnGtlcYLsE4A3zGwoQMgcATP7EnhQ0mvAVOAC4KastjLHPF3EOVcKvMc598rLyxk0aBAjRozg8ssvp3fv3lRVVXHwwW0vi3LFihXcd999fPrpp1RWVtKrVy/233//tXd3nXMtEyfI/hZwZ8prAzZY+8JspqT/AqfhQbZzzhUtSVcTrgF/MbMF0eumMDP7dQ6b1ibMmTOH7bffniuvvDLfTWnQl19+yamnrn/zuUePHtx9990ccMABeWqVc8UjTpC9Clie8noJsFlanU+Ao1vaqNbmy6o754qZpMHARsCDZpYxb0FSO+BEYIGZNbbi4zBCkP0vYEH0uimMkHJY1GbPns0222yT72Y06PTTT2e//fZj1113ZcMNN2TmzJnceuut3HHHHVRVVfHqq6+y22675buZzhW0OEH2bMKMIkkfAHun1elDOOEWFO/Jds4VK0nfBR4H/mpm/6yvnpmtlLQ7cI6k3mb2dgO7PTB6/jTttQM+++wzevfune9mNGjo0KHrvf7Od77D7bffTufOnbnpppsYNmwYjz32WJ5a51xxiBNkvwIckvL6ceA3ku4CHiUMfjwEeCBbjWstHmQ754rYT4CVwPAm1P01cGb0uLC+Smb2YkOvS9nKlSuZO3cu2267bb6b0iw//elPuemmm3jppZfy3RTnCl6cIPsBYFtJ3c1sFvAnYAhwOiEPW8AM4LLsNjH3PMh2zhWxgcBYM5vXWEUzmydpLPGWVHcp5syZA9Dm00Xqs/nmmwNhbm/nXMs0Ocg2s7HA2JTXSyUNIATaOwCzgCfNbGl2m5h7npPtnCti3wKej1H/XWDfHLWl6H322WdA4QbZr776KgDbb799nlviXOFr0bLqZrbazP5tZteb2b8KMcAG78l2zhW1doR0kaZaSebVfBskaStJf5E0Q9Ky1lxWXdIxkqZLWiHpY0kXZWjbjZKmSaqR9JmkeyR9I9ttmT17NtC2g+y333474zLvn3zyCeeddx4Ap5xySms3y7miE2tZ9WJVVhb+1/Ag2zlXhBYA28Wovx1h1ccmk7Q1YSGbLYC3CUH6J8AKYHvCtWYqsDjOfpt47AGEcUH/AH4J9Aeul7TGzP4UVesHfA/4OzAxaucwYLyk75hZTbbakwyy23JO9sMPP8x1113HgQceSI8ePdhwww356KOP+O9//8vy5cs54ogj+OUvf5nvZjpX8OoNsiX9qLk7NbN7m1JP0gmEFSL7AV2B94Hfm9nIlDoCLgd+BmwKvA783MymNrd9GdpBeXm5B9nOuWL0BnCopPZmtqKhipI6AIcSzrNxXA1sCQwys+ei5dvvNrNrJG1DWGOhO5CLFVmuBl42szOj189I2gi4WtJtZrYSeBnoaWZrT/KS3iBcc44D7slWY2bPnk2XLl3YcMMNs7XLrDvwwAN5//33mTJlCq+++iq1tbV069aNfffdl1NPPZVTTz117YJzzrnma6gnewSZl99tiKL3NCnIBi4CPiaMYp8HHAE8IGlTM7slqnMZcBVwCfBe9J7not6HL2O2r14VFRWek+2cK0aPEc6tvyX09DbkGkJnxr9jHmMQ8LSZPZe+wcxmRx0qbxFmOPl5zH03pjfwl7SyZwjXlb2BF81sUYZ2fSBpKbB5Nhvz2WeftelUEYADDjjAF5txrhU0FGSf3grHPyptxPuYKEfuIuCWqFflMuBaM7sVQNKrhEGW5wFZW06roqLCe7Kdc8XoXuBi4EJJWwDXmNmHqRUk7UjozPghYeBjvfNp12NL4KGU1wmgY/KFmdVIepYwUD7bQXYH6uacJ3vsdwYyTi8oqRfQCXgnfdvkyZMz9uSaNd7vNHv27DadKuKcaz31BtlmlrXbZw0cI9OUUlMIJ2KAfYAupJy8zaxW0pNAFR5kO+dcg8xslaRjCDOM/BD4gaTPCQuMGbBN9FBUdkxqWkUTfU0YYJm0kPUXL4OQj52+SnA2zAD2SCvbM3reONMbJJUBfwY+JPR6Z83s2bPp1atXNnfpnCtQLZpdJEf2YV3PQk9Cj8iHaXXejbZllOyFSH80pLy83NNFnHNFycw+IIx9uYvQ67sNsBchnWLbqOzvwO5mNqMZh/gk2k/SNOAgSZ1gbVB7GCGIz7bbgSGSzpK0kaRBhJ57CNePTK4lfPZTMy0z369fP8yszqMxK1eu5Msvv2zz6SLOudYRe3YRSZ0Jo7T7EAYrLib0Pj/W0hHakg4m9GL/JCraCKgxs/QT5UKgk6R20aCWFvOebOdcMTOz/wFnSTof2B3YitB7PQeYZGbLW7D754GzJW0QBa33ENJUxkdpIvsCuwK/a8lnqMc/gN2AvwJ3AEuBS4FbgLnplSWdQxjjc7KZTcxmQ7744gvMzNNFnHNAzCA7GrxyO9CNcHJOMuBPkv6fmT3SnIZI6k5YVfIJMxuRtu861RvYRr9+/Zg0aVKs43uQ7ZwrBVEw/XKWd3sXofNjU+ALM7tPUj/gfCCZO/EgYfBlVkWdMOdJuorQQ/8x6+50TkitK+k4QvD9KzP7V7bbUghzZDvnWk+T00UkHQqMBDYk9FCcTsiLPp0wSGZDYKSkQ+I2QtLGQDXwKZA6A/5CYENJ5Wlv6QYszXSbr7k8yHbOueYxsw+jRcm+SCm7kNBbvjewlZn9oIW95Y21YaGZvRndUT0HGG9m7yW3SxoI3A/cama/z/bxE4kETzzxBAAzZszw9EPnXKye7KsJI7b3M7M30rbdI+lW4KWoXp1pnOoT5eyNIgyaOdLMalM2vweUE5Ztfz+lvGe0LWs8J9s557LLzL4CvsrlMSTtRUhHmUoYKH8yYUrBfVPq7Aw8Trhu/Ct6T9JXZvZRS9qQSCQYNGgQ48aNA+DSSy/lscceY/To0ZSXp/cROedKRZyBj32Af2UIsAEws0mEWUD6NnWHkiqAh4EdgaooZzDVeMKo9RNS3tMJOIrQ85013pPtnHMFaRVwIiGIHkGYlm+AmU1PqdOfMIZoN+AV4NWUx1UtbUB1dTUTJ05k5cowRKi2tpaJEydSXZ3Vy5RzrsDE6cleAXzRSJ05rJuftCluIyyS8Atg47TehSlmtlzSdcBVkhaybjGaMkJeXdZ4kO2cc00j6R/NfKuZ2RnZbIuZTabuFH7pdUYQAvCcmDJlCrW1teuV1dbWMnXqVAYPHpyrwzrn2rg4QfY4Um6/1WMAIWWkqQ6Lnv+cYVsPwqIz1xGC6suBTYBJwKFmVmfUeEv4surOOddkpzXzfQZkNchuC/r06UNlZSU1Nesm2KqsrKR37975a5RzLu/iBNmXAq9GPcu/Ts2dllQJDAW+Q5jnuknMrHsT6hhhRHrWR6Wn8mXVnXOuyXrkuwFtSVVVFf379+eFF15gzZo1dO7cmf79+1NVVZXvpjnn8ihukD2dML/o2ZLeIMxBugUhD7sroRf70rSFX7J+ezAXPF3EOeeaxsw+yXcb2pLy8nJGjx7N5ptvTvfu3Rk+fDhVVVU+6NG5EhcnyD4t5etuwEEZ6hwQPVIVxO1BD7Kdc6VGUhWhA2XXqOgt4DozG52/VhWmsrIyamtrOfjggz0P2zkHxAuyi/r2oE/h55wrJZLOAv4GzADGAB2A/YCnJJ1pZnc3Y5/7EtZOSF0R+A1ghJllewGcNmXp0qWsWLGCTTfdNN9Ncc61EU0Osov99qD3ZDvnSswVwF/M7PxkgaSuhNUgrwBiBdmSbiEsAqO0Tb2B0yX9xcx+3qIWt2Hz588HYJNNNslzS5xzbUWcebKbRFKXbO+zNXiQ7ZwrRpKul9Q+w6ZtgUdTC8xsMfBMtC3OMc4HziUsaX464c5nx+j5J1H5uZLOjf0BCsS8efMAD7Kdc+vEWVZ9tKQG74NJ6k9YdavgeJDtnCtSZwBTJaXP/PQh8P8kdUwWSOoOHAt8EPMYPyWsk7C7md1jZp+Y2YroeQSwJ/Aloae7KHlPtnMuXZye7EOBaZIOzLRR0i8Js4tslY2GtTbPyXbOFaldgLeBlyT9KVo1F+BKwmq6cyRNkDSFEFxvG22LY3vg32a2KNNGM1sA/DuqV5Q8yHbOpYsTZJ9FGMjyjKThiubpk7SJpP8C1wOf0viCNW2S92Q754qRmf3PzI4HTooeb0o60Mz+TVgp8SlgA8JMUP8C+pnZf2IeZj6wspE6K4F5MfdbMJJBtg98dM4lNTnINrO7CLf83if0crwg6XhgGlBFODn3jZa4LTgeZDvnipmZPULo1Z4APCfpduBDM/uhmfUzs75mdqqZTWvG7h8Hjpa0QaaNktoBR0f1ilIyyN54443z3BLnXFsRa+Cjmb0D7A6MAPYnBNZdgTPN7AdmtiTrLWwlvqy6c67YmdkCM/shMAQ4Enhb0uFZ2PUVhOn6npO0T8qdTkkaADwHLIzqFaV58+bRtWtXKirizIzrnCtmzTkbbAR8K/o6OVWTZac5+ePLqjvnSoWZjZL0EvBHwrzY9wIX1JdT3QRTgXaEMTnjgNWS5gGbsu468wVhXE9aU+xbFIH58+d7PrZzbj2xerKj1cGmEnqx/w4cTui9+LukByRtmPUWthJPF3HOFStJVZJGSXozej7CzL42szOAQYSVet+RdEwzD1EGrCKMy/mUMNPIyug5WbaK0DGT+sj6NLL54kG2cy5dk3uyJd0IXAjUAieb2b+i8t2AewgDavaQdFIh5mV7kO2cK0aSTgRGAosIs4fsA1RJOsXMRprZs5K+A9wA/FvSw8B5ZtbkQYpm1j37LS8s8+fPZ7PNNst3M5xzbUicXoSLgSmEwY3/Shaa2XwzGwxcAmxHWC2s4PgUfs65InUFYQq/7ma2F9AdeBe4LFnBzGrN7FzgQKBvtN3FMG/ePO/Jds6tJ06QfQuwj5l9lGmjmd1EmL5vTjYa1tq8J9s5V6S2B6rN7GuA6PkpMsxZbWYvAb0Ig9ubrLGFylLq7R5nv4XE00Wcc+niTOH3CzNb1Uid14E+LW5VHniQ7ZwrUh8DAySlnu/3AWZlqmxmy83skpjHmCrpgIYqSLqIAr3T2ZiVK1eyZMkSD7Kdc+vJ+qCTZG9JofEg2zlXpK4H9iYMbHxA0juEIPuGLB5jY8L0fVcrbfoQSRtL+g/we+DzLB6zzViwYAHgqz0659bXYJAtaYykH6WV9Zf083rq/0LSzGw2sLV4TrZzrhiZ2f2EebFnArsBnwDfM7N/ZvEwexIGVQ4Fnpe0JYCkfQkzUg0GHqFA73Q2xld7dM5l0lhP9kDCIJlUhxPmVs2kG/DNFrUoT7wn2zlXrMzsSTM7wsx2NbMqM3siy/t/i7BQ2b2E68Y0SbcCYwhzZf/MzL6fqzudko6RNF3SCkkfR6kp6XUk6QpJn0laJuklSb2zcfx588JELN6T7ZxLVTRzlLaUB9nOOdd8ZrbMzE4nzDS1GfAzwiqPe5jZ33J13GhFyUeB14CjgH8A10u6IK3qZcBVhPSZo4AaQorLli1tQ7In24Ns51wqD7IjyXQRs4JfvNI55/JC0mGEIBtgCaEX+1eSKnN42KuBl83sTDN7xsx+TZgN62pJ7aJ2dSAE2dea2a1m9hxwAmG14vNa2gAPsp1zmXiQHamoCOvyrFmzJs8tcc65wiKpXNJ1hKkBOwE/AHYAngFOBSZnKzUjg97Ac2llzwAbEQZ8Qhjo2QV4KFnBzGqBJ4GqljbAg2znXCYeZEeSQbanjDjnXGzjgF8B0wgLlj1oZvPMrIrQg9wDeLW+QfMt1IGwhHuqFdHzztFzTyABfJhW791o23omT56MpDqP+syfP5+OHTvSqVOnZn0A51xx8iAbSCQSfPDBBwA8+eSTPsuIc87FsxdwK7C3mc1I3WBmNwD7A19S/6D5lpgB7JFWtmf0vHH0vBFQY2bpJ/eFQKdkWklz+WqPzrlMKppQZ2Daf/ADASRdBaT/az8wK61qRYlEgkGDBjFu3DgATjvtNG6//XZGjx5NeXl5nlvnnHMF4Tgze6y+jWY2UVIf4O85OPbtwF8lnUWYJnBP4OJoW2pQnWnAjTJt69evH5MmTWpyA3y1R+dcJk0KsskcPA9Pe22EE1ZORg5K2oUwmGVvYBHhZD08Q89ELNXV1UycOJGVK8PdxtraWl555RVOPfVUdthhB8rLy+nduzcA06dPp0+fPlRVVXkA7pxzkYYC7JQ6i4Djc3D4fxDm//4rcAewFLiUcL2YG9VZCGwoqTztmtENWNrYasaN8SDbOZdJY0F2eiCdF5I2IgxseYewqMK3gJsI6S5XtmTfU6ZMoba2dr2y5cuXM3LkyLWvkwF1IpGgQ4cOfOMb3+DGG2+kvLycqVOnkkgkKC8vp1+/fh6AO+dKlqTNgOMIudCVZnZmSnkP4E0zW5bNY0ZB83nR3dVtCMvIJ/OsJ0TP7wHlhMGY76e8vWe0rUXmz59Pr169Wrob51yRaTDINrM2EWQDPwU6AsdGixk8K6kLMEzSDS1Z4KBPnz5UVlZSU1NTb53UHO3ly5czc+ZMjj/+eCStNxtJZWUlO+ywA0OGDMHMkISZrdcbnhqU57qsNY9fSp8138cvpc/anOP7P7v5IekM4GbCQMTkXc0zo81bAK8CZwN35eL4ZraQ0GONpHOA8WaWDKDHA18Tpu37TVSnE2G+7DtactxEIsGcOXPo0qULo0aN8t8959xaTUkXaQuqgNFpwfSDhEUFDiBMw9S8HVdV0b9/fyZOnNhgoJ3OzOrMqV1bW8u0adOYNm1anfplZWV13lNeXo6ZrReoZ7usWI9V6scvpc8a9/idO3emf//+Pq6iFUk6lBCsTicsrT6I0DkChBUhJb0NHEOWg2xJewHJ5du7ACdHx9835fjLoykGr5K0kNB7fRHhbugtzT12IpHgsMMOY8mSJbz++uucfPLJ/rvnnFurUGYXqXNLz8w+JeTetWj6pfLyckaPHs3IkSM55ZRT6NChQy7az5o1a+oE5YlEok7gkO2yYj1WqR+/NY9VaMevqalh4sSJVFdXZ3yPy4lLgS+AA8zsP8D/MtSZDuySg2OvAk4EHgdGEObpHmBm09PqXQf8FrgcGEUIyA81s7k0U3JMT5L/7jnnUhVKkL0RYbBjuoXRthYpLy9n8ODBjBgxggEDBlBZuf7iZGVlZZSVFcq3yjlXW1vL1KlT892MUrI7MKqR1L3ZQIuXME9nZpPNbA8z62xmXczsSDN7M0M9M7Pfmtk2ZtbRzPYzsyktOfaUKVNYunTpemX+u+ecSyqUdBGof/qlOuVxp19KSvZqV1dX88Ybb7B69WoqKirYbbfdSCQSXHLJJcyZM4fly5c3eAvbOZdflZWVa/O1XatoB9Q2Uqcb60+pV/Ayjenx3z3nXFKhBNkLCSfodF3J3MPdbMle7cGDB9fZNmTIEKqrq5k6dSrf/e53AZg2bRqrV6+mrKyMxx9/nA8//LBOzwawtic8NShvjbJiPVapH7+UPmvc41dWVtK/f3+qqlq8WrZrullAv0bq9Gf9mT0KXuqYntraWv/dc86tp1CC7PdIy72WtC1QSQunX0rmaqfnS2eSKQAfMmTI2q+vuuqq9XrBy8rKWLNmzdrecFgXlLdGWX3HHzp0KACPPvpozo+V78+aWv+YY45Z+3PKx/HzcaxS/Fn37duXo446ioqKiib9XbuseAL4laQTzOzh9I2STgd6Af/X6i3LoeTdz4qKcCkdOXJkycwuEufaWUxK8XOX4meG7HxuFcI3TdLlwCXAN81sSVT2S+AaYMv0PMDdd9/dmpouUoq/PKX4maE0P3cpfmaI97klTTaz3XPdpmIWrWXwBrAt8G/CXcZDgV8A+wHHAh8B/cyssbSSvItzDYGm/b4NHDgQgLFjx7akaW2Gn1tK53OX4meGpn/uhq4hhdKTfTvwc+BRSdcD2wPDgD+0ZI5s55xzLWdmCyUdANxLmIs66eboeRzwg0IIsJ1zLluaHWRLqiJM27RrVPQWcJ2Zjc5Gw1JFJ/CDgVsJc2IvAv5ICLSdc87lWTSt6kBJvYC9gU2AxcAEM5uc18Y551weNCvIlnQW8DdgBjCGsMLXfsBTks40s7uz18TAzN4BDsr2fp1zzmVPND91+hzVzjlXcprbk30F8BczOz9ZIKkr8HK0LetBtnPOOeecc4WiwYGPUf7z1Wa2Iq18NWGlrBfSym8CzjWz3Cyb2ESSvgI+yWcbnHMF45tmtlm+G+HaDr+GOOdiqPca0liQPQ/4CjjDzManlL8LTANON7NlUVl34AVgiZn1yl7bnXPOOeecKyyNrRW+C/A28JKkP0nqFJVfSRhBPkfSBElTgA8I0zddmbPWOuecc845VwCaNE+2pOMJM3vUAmea2QuS+gIXExaJESEY/72ZTcthe51zzjnnnGvzmrwYjaSNgVuAk4A7gUuSC8M455xzzjnn1om94qOkwcBfAQPONrOnc9Ew55xzzjnnClVjOdl1mNkowgI0zxLmxR4hqVu2G+acc84551yharQnO1rZ8Vzgm4QpjW4zs6eibYcCdwDtgXPM7PGctraFJHUGDiDkkW9E6I1fBLwHvGhmNflrncslSRWAmVki321xzhUuv46UJr+GuOZobAq/E4GRhBPIB8C3ga7AKWY2MqpTCdwA/BR4GDjPzOblttnxSBIwHLgI6AQsBRYSBmx2BSqjspuAYRY3h6aASNoQOB3oS7g4TALuKbYLg6TNgfOAKsLFMDkzzlLCxfAp4FYz+yo/Lcy9UvlZgwc+Lvf8OrJOKZxb/BoSlMLPOikX15HGguxphJSSAWb2taQuwHggYWa7pdXdH/g7sFFbW9hB0nDCTCjDgQfN7LO07dsQBnQOBW4ys2Gt3sgckDQHOMLMpkavtwNeArYG3o+q7QTMAvYzsy/z0Mysk7Qb8BzhD+RJ4F3WXQy7Ef6AjopeHxItA13QSvhn7YGPaxV+HSmdc0spXkOgNH/WkOPriJnV+wCWADekld1AWHAmU/0OwI0N7TMfD+BzwiDNxuqdDXye7/Zm8XOvAfZMef0wMBvYNaXsO1HZXflubxY/94vAf4FODdTpFNUZm+/2+s+6RZ97OFADXAJsm2H7NsAvo3PZsHy31x+F+/DryNrXRX9uKcVrSKn+rKPPlLPrSGMHng68ApSllL0MvJnvb0rMb2AtcHAT6h0M1Oa7vVn83Ol/MIsI85yn1/sZMCff7c3i514KHNSEegcVy8+7hH/WJRn4+KP1H34dWfu66M8tpXgNKdWfdfR5cnYdaWx2keuBvYF3JD0g6R1gH0JvdiGZAPwqyh/PKNp2KfBqq7Wq9XUk3PZK9y6wSSu3JZfmEW5pNaYnMD/HbcmXUvlZdwM+akK9j6K6zjWXX0eCUji3+DUkKIWfNeTwOlLR0EYzu1/S14T/WnYjzC5yuZk9EecgbcB5hPyqTyWNJiSxLyLkW3Uj/KEMAlYQeiGKyfmSvoi+XkK47ZHuG4TvR7G4Hfh9tIDSw8CHFv0bGuVe7QCcAFwB/DZvrcy+UvxZJwOfCWZWm6lCiQQ+Lvf8OhKUwrmlVK8hUHo/a8jhdST2YjSFKprL+2fA4awbOQohuf09oBq43cwW5aN9uSBpFuECkOpxM7swrd4DwGZmdmhrtS3XJF1O+IPYEEgQ8q0sel1OOHlcZ2bX5a2RWVSqP2tJOxMCnw5Ao4GPmb2Xl4a6ouDXkbVK4dxSUtcQKOmfdc6uIyUTZLv6SToC+MzM3sx3W7JJUntCelOmi+F4M1uRr7blSzH+rEsx8HGurSnSc4tfQzIo0p91N3JwHfEg2yGpHbClmX2a77bkmqTuhAUFPsl3W7JJUn/gdTNbk++2OOdKT6lcR4r1GgJ+HcmF2Muqu8Ii6VxJH0laImmipFMzVOsDfNzabcsVSWdL2jSt7BeSviIMXJgpaa6kc/LTwpx4FfhC0q2S9s13Y/JJUoWk8ny3w7liUWrXkRK9hoBfR9bK1nXEg+wiJukk4BZCUv9wYA4wQtIjkjrmtXG59Vdg++QLSWcDfwSeB06MHi8At0j6QV5amBtTgZOBFyV9KukGSX3z3Kack7S5pGskvS5pCSFvbmUUELwuabikNrVAlnOFokSvI6V6DQG/jmT1OuLpIkVM0iRgjJn9KqXsYOB+wopNR5rZ/OgW0XgzK4reP0lrgL3M7LXo9XvABDM7La3eP4GeZrZH67cyu5KfGZhCWAb4JGAwYaWqGcADhFXq3q93JwVIJboym3OtpRSvI6V4DQG/jpCD64gH2UUs+m/sKDMbm1benZDEX074Q9qUIjk5QsYT5CpgsJmNTqt3JPCQmdU7722hSP/MUVlH4GjCifJwoB0wDRhpZjfmpaFZJulFwqj/E8xsaT11OhGm4ao0s4Gt2DznCl4pXkdK8RoCfh0hB9cRTxcpbosJJ771mNkswojpecB4oCj+C0/TQVKn6A9jPmEKpnQJ6k5XVDTMbJmZ/cvMvgdsAZwB/A/4TX5bllV7ADfVd2IEiLbdRHH+njuXa6V6HSn5awj4dSSpudcRD7KL22TgmEwbzGwhYcGEScDNrdim1vICYR7TJcDmwJ4Z6vQCPmvNRuWLmX1tZiPM7HDCQgLFwldmcy63SvU64teQNH4diX8daXDFR1fw7gMulLSxmS1I32hmyyQdTRjkURSTykdOz1D2RYayPYF/57gtreVF4OumVDSzYgo2S3llNudaQyleR0rxGgJ+Hcn6dcRzsp1zBU0luDKbc8657MnVdcSDbOdcwfOV2ZxzzrVELq4jHmQ754qepA7A5sW+Gp1zzrncaM51xAc+OudKwZEUyWp0zjnn8iL2dcSDbOecc84557LMZxdxzhUsSWOaWNWXVXfOOVdHLq8jHmQ75wrZ/sD7wDuN1OvQCm1xzjlXeHJ2HfEg2zlXyN4C3jezExuqJOl44F+t0yTnnHMFJGfXEc/Jds4VsonAXk2oZ4By3BbnnHOFJ2fXEZ/CzzlXsCR9C9jVzP7TSL2OhKmXPmmdljnnnCsEubyOeJDtnHPOOedclnm6iHPOOeecc1nmQbZzzjnnnHNZ5kG2c84559YjaaykosonlbSjpMckfSnJJC3Kd5tccfMg261HUk9Jt0h6S9JiSSslzZH0X0lnSPL5hvNM0mnRBeK0fLfFOVe/6O/UJH1S37lT0qyojk+pm0OSyoHHgSOAUcBw4Lp8tsllR1u+JvoftVtL0tXAUMI/XxOAe4AaYAtgIPB34GfA7nlqonPOFaLtgAvwoC6fegC7AHea2dn5bowrDR5kOwAkXUH4z/4z4AQzm5ihzmDg4tZum3POFbCFhPl1L5f0dzObl+8GlahvRM9z8toKV1I8XcQhqTswDFgFHJEpwAYws1HA4Wnv/b6kl6LUkmWS3pR0uaT2GY4zK3p0lvRHSZ9F75kq6ZioToWkKyR9KGm5pI8knZdhXwOj20PDJO0t6bmoDUskjZaUsbddUldJ10p6P9r/wqj+IY0co3eUMrNI0lJJL0rap55jVEg6R9IESV9H9adIOk9SWVrd7tExRkRfPyhpXtS2SdE/Nqn1xwJ3Ry/vTrkdbdHP0TnXtiwFfg10IdwpbFTquaee7bMkzUorW3vLXNKhksZJqpH0laS7JXWL6vWRNCo699VI+k9D5w5J7SX9RtLHklZE5+ShktrVU79ndD77LKo/V9IDknbKUHdE1ObtJZ0vaXp0TRjbxO9TP0n/lvS/6FifSLpN0lZp9Qx4MXo5NOWcOayJx9lT0r8kfR4d5wtJz0j6foa6fk2s+768XhMlbSjpKoU02K+j78lH0c+0X6Y2Z42Z+aPEH4QebANGxnzf76L3fQX8FbiRsDypAWOBDdLqzwI+B14F3gduBe4AlgAJ4GDg38Bs4M5o+9xofyem7WtgVF4NrAD+G7XnIWA1sAzYL+093YC3o/e9Rrh1+3fga2AN8P/qOcYowoXyeeD30TES0TF2SnvPBsDT0fveA24H/gRMi8r+mVa/e1T+AvA/QprOHwmpOsuj4xyYUv80Ql6hRc/DUh7d8v275A9/+GPdI/o7nR2dF2YAK4Fvp9WZFdWrSClLnnuG1bPfWcCstLLTovc8Gh3n0eh8NT7lnLwXUBudo34PjI62vQ2Upe1vbLTtCULv783ATdHnMOBJorU2Ut5zeHSuXBUd/wbggehcthjom1Z/RMq+FgH3R+fl3zbhezuYcO5fGR3jWuCZaH+fA91T6g5LOdbYlHPmwCYc5yzCNWUF8DDhOvN3YCowNq2uXxPb2DWRsELjK9H28cAfCL+XI4EvgPNyeg7I90nIH/l/RH8oBpwZ4z17R+/5FNgypbyCcMI04Iq098xi3Qm1fUr5flH5AuB1UoJFYHvCSXRK2r6Sf+yW/kcCDInKPyTlwgH8LSr/GykXB2BHwgVgBeufmFOPcVraMf5fVH5bWvmwqPwWoDylvBy4K9o2JKW8e8oxhqbta1BU/lRa+WmZ2uQPf/ijbT2iv9PZ0dfHR68fTauTPC9mK8heDRyQUl4GPJtyjv1h2vvqnJei8rFR+QfARinlHQhBoQGnppRvREiNmQfskravXQnje95IKx/BuqC4R4zva+foOAnqBo6XRvt8Jq28we9pPcfZhfAPwwLCioDp27dJ+dqviW3wmgh8N9r2WIZtZaT8bufi0SonGn+07QfwTvRLeHiM99wZvefsDNu+HZ38ZqaVJ08o38rwnpnRtoMybHuBcOFI/QNN/rGvd9JI2T422n5A9HoDQg/OEmDjDPV/HdW/OsMxXs5QfwPCyXdSSlkZ4cT/BSkXzJTt3Qi9Aw+llCVPKLNSP1/K9k+AeWll9Z5Q/OEPf7SdBylBdvQ62au8b0pZ8ryYrSD7nxnq/yja9lKGbQeQOaBJnkNPzfCeZPteSCn7RVR2bj1t/mO0fZeUshFR2S9ifl9/GL3vgQzbKoCPo+3bNfV7Ws9xbonec2ET6vo1sQ1eE1kXZNf5XWmNhw98dBBup0D4RWyqvtHzmPQNZvaBpNlAD0ndzGxRyuZFZvZRhv3NIYz+npxh2+eE/3q3jL5ONc7M1mR4z1jCxaMPIRevJ9AJeMXMFmSoPwa4MqqfblJ6gZmtkjSX0HuT9G1gE8JJ7kpJ6W+DcDtt5wzlU80skaH8M0IPiXOu8F1MCLRvkrSXRVFAltU5X7FusF9951eAberZ34sZysYRgrzU82XyPLVbPbnO346edyZ07KR6rZ5j16eh689qSS8RgrU+hJ7l5toreq5uYZv8mpi/a+I7hNSekyV9k5D+9DLhn4GVMfbTLB5kOwh/zD2p/ySbSdfo+Yt6tn9BmLaqKyHXLmlxPfVXA5hZpu2ro+cNMmybW8/+vkxrZ1PaC+E/63SL6nnPasKJLmmT6HlHGh7g1DnmMXyAsnNFwMxelfQIIXXk+8C/cnCYhs6hcc+vkOEca2YJSfOBzVOKk+e/sxppX6bz35cZyhrSkvN5HMn3pweymfg1sQ1eE6Pf1YOAqwl/d9dHm5ZIuge43Mxqmrq/uPzi7SD8VwdhkEVTJf/wt6xn+1Zp9XJli3rKk+1anPacy/Ym3/uYmamBR48WHMM5V9guI9xWv1b1zNBBuIUO9XeEda2nPBfqnGMVFnbZhDBALil5/tutkfPfPRmOEbdHv7WuP4ui562bUNeviXW1iWuimS00swvNbFtCwH8mYRDmeYQBqjnjQbaDMP3NKuA4Sbs0VDFlGqIp0fPADHV2IPSKf5x2WywX9k2fAiitXcl2vk8YDd1b0kYZ6h8YPb/Rgra8Rzgp7yWpvl6hbEjeQitvsJZzrs2JUgNuI6QCnF9PtYXR87bpG6Lza7ecNC6zAzKU7Uf4B2BKStmElG251tD1pwLYN3rZkvM5rPtMVS1sk18T28g10cxmmNldhN/rGsKg0JzxINthZrMII4DbAf9tYD7Nw1mXm/aP6PlKSZul1CknTOlTRhg5nGs7AuektXMI4Q9oBiF3kCj36n7Cbalr0up/C/g54R+Nfza3IWa2mjBQZivgZkkd0+tI2qqxf2SaYH70vF0L9+Ocy49rCMHH/5H5Vvl7hF7iIZLWpmRE55SbW6OBKa5KDcIUloe/Nnp5d0q9uwmfaaikPdN3IqlM0sAstelxwswbJ0vaK23bBYQZOJ4zs5bkY0Po5VxN+B7UOW9LSk2x9GtimrZwTZTUQ9KuGd6zEdCekBOeM56T7QAws99FPQBDgdcljScMbkguq74/4Y93UlR/vKQbgF8Bb0V5hrWE//i/Q0hBubEVmv40YRBRFWHezR2AYwnzaZ6RNgDkMkIvy3mS9iCM0N6UkBu5IWHao49b2J5fA7sBPwWOkjSGkM+3OeH7N4BwYU0f+BPHq4QeiAskbcy6HLxb6snfc861IWa2QNLvCPP1Ztq+StKfgauAKZIeI1yvDyWMoWnNVQvfBd6OzvGrCD1/3yLMw7w2ADOz+ZKOBx4DJkh6njAH8xpC8LM3IcWkQ0sbZGY1kn5CmLf6RUkPEwY49gMOI+Qf/78sHOcdSecQ5naeIukJwiC+TYDdCTNzHBjV9WtiZnm9JkbHfkzSZMKc5XOAzQi/xxuwLkc7N/IxpYk/2u6DMMr3luiX8WvCfJxfEHqwzyBlLs+o/kmEk8cSwh/x24Q/mA4Z9j2LtGmnUraNDb+OGbeNIOTsdU8pGxiVDSOcvJ+L2ruEsCDBHvXsqxvhj+pDwhygiwhzyB6Woe7aY9Szr4yfhzBby6mE+ccXRN/Dz6Pv0xXAtil1u0fHGBHn+0JY9OFVwj9Blv798Yc//JH/B2lT+KVta8+6qeaMtCnOovPIZcBH0TnkU0JQ3inTuYeGpzGr91xW3zmIdVO+tQd+E7V1BWFquaHp14K0/d0anWOXR+fl9wgB+TFpdeuc22N+f/cgBPVfpXyP/gp8I873oAnH2ZuwKMz/ouPMIQSzx2eo69fEuuV5uyYS0nR+R1iQ5svoM84mxDRV2fpbr++hqHHOFZTotuMLwHAzG5bXxjjnnHN55NfEtslzsp1zzjnnnMsyD7Kdc84555zLMg+ynXPOOeecyzLPyXbOOeeccy7LvCfbOeecc865LPMg2znnnHPOuSzzINs555xzzrks8yDbOeecc865LPMg2znnnHPOuSz7/1LAohZ8L/leAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Perform PCA on the total encoded data including probes\n",
    "\n",
    "pca = PCA() #n_components=2\n",
    "pca.fit(all_search_linear_probes_array) \n",
    "\n",
    "\n",
    "fig = plt.figure(figsize=[12,3])\n",
    "fig.tight_layout(pad=10.0)\n",
    "fig.subplots_adjust(bottom=0.0, hspace=.5, wspace=.5)\n",
    "\n",
    "ax = []\n",
    "for a in range(2):\n",
    "    ax.append('ax' + str(a)) \n",
    "    ax[a] = fig.add_subplot(1, 2, a+1)\n",
    "    ax[a].tick_params(axis=\"x\", labelsize=15, direction='in', length=6, width=2, colors='black', rotation=90)\n",
    "    ax[a].tick_params(axis=\"y\", labelsize=15, direction='in', length=6, width=2, colors='black')\n",
    "    ax[a].tick_params(which='minor', direction='in', length=3, width=2)  \n",
    "\n",
    "\n",
    "#print(\"Explained variance: \", 100*pca2.explained_variance_ratio_)\n",
    "\n",
    "ax[0].plot(100*pca.explained_variance_ratio_, color='black', marker='.', ms=10)\n",
    "ax[0].set_xlabel('Component', fontsize=20)\n",
    "ax[0].set_ylabel('% Explained variance', fontsize=20)\n",
    "\n",
    "ax[1].plot(np.cumsum(100*pca.explained_variance_ratio_), color='black', marker='.', ms=10)\n",
    "\n",
    "ax[1].set_xlabel('Number of components', fontsize=20)\n",
    "ax[1].set_ylabel('% Cumulative \\n explained variance', fontsize=20)\n",
    "ax[1].axhline(y=95, color='black')\n",
    "ax[1].axvline(x=4, color='black')\n",
    "\n",
    "ax[1].annotate('5', (5, 95),\n",
    "            xytext=(0.25, 0.25), textcoords='axes fraction',\n",
    "            arrowprops=dict(facecolor='black', shrink=0.05),\n",
    "            fontsize=20,\n",
    "            horizontalalignment='right', verticalalignment='bottom')\n",
    "\n",
    "print(np.cumsum(100*pca.explained_variance_ratio_))\n",
    "\n",
    "plt.savefig(figures_dir + '-PCAVariance_total-zoom.png', bbox_inches = \"tight\")\n",
    "\n",
    "#plt.savefig('../WeeklyAnalysis/Figures/PCAVariance_total-zoom.png', dpi=300, bbox_inches = \"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b16cabbd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d102896a",
   "metadata": {},
   "source": [
    "Transform latent points to PC space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "436f36b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "## search set\n",
    "\n",
    "projected_all_search_linear_probes_array = pca.transform(all_search_linear_probes_array) # search set\n",
    "projected_zinc_set = pca.transform(latentspace_training) # ZINC set\n",
    "\n",
    "\n",
    "## save the arrays (ONLY ONCE)\n",
    "\n",
    "#np.save('../Data/VAE_data/PCA_AllSearchLinearProbes.npy', projected_all_search_linear_probes_array)\n",
    "#np.save('../Data/VAE_data/PCA_ZINC_TrainingSet.npy', projected_zinc_set)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b26a8fe",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a75d9df1",
   "metadata": {},
   "outputs": [],
   "source": [
    "## read saved PCA arrays\n",
    "projected_all_search_linear_probes_array = np.load('../Data/VAE_data/PCA_AllSearchLinearProbes.npy')\n",
    "\n",
    "projected_zinc_set = np.load('../Data/VAE_data/PCA_ZINC_TrainingSet.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbe4aae0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d19f086e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "98b50743",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}