bpc_interactive.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basis profile curve identification to understand electrical stimulation effects in human brain networks\n",
    "\n",
    "Method developed in Matlab by Kai J. Miller, available at Miller KJ, Mueller KR, Hermes D. Basis profile curve identification to understand electrical stimulation effects in human brain networks. doi: https://doi.org/10.1101/2021.01.24.428020\n",
    "\n",
    "This project was supported by the National Institute Of Mental Health of the National Institutes of Health under Award Number R01MH122258. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health\n",
    "\n",
    "This Jupyter Notebook was written by Tal Pal Attia, Harvey Huang and Dora Hermes (2021)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview\n",
    "\n",
    "This notebook is a walk through of the methods presented in [Miller et al. (2021)](https://www.biorxiv.org/content/10.1101/2021.01.24.428020v1).\n",
    "\n",
    "> **Abstract (Miller et al. 2021)**\n",
    "<br /> Brain networks can be explored by delivering brief pulses of electrical current in one area while measuring voltage responses in other areas. We propose a convergent paradigm to study brain dynamics, focusing on a single brain site to observe the average effect of stimulating each of many other brain sites. Viewed in this manner, visually-apparent motifs in the temporal response shape emerge from adjacent stimulation sites. This work constructs and illustrates a data-driven approach to determine characteristic spatiotemporal structure in these response shapes, summarized by a set of unique “basis profile curves” (BPCs). Each BPC may be mapped back to underlying anatomy in a natural way, quantifying projection strength from each stimulation site using simple metrics. Our technique is demonstrated for an array of implanted brain surface electrodes in a human patient. This framework enables straightforward interpretation of single-pulse brain stimulation data, and can be applied generically to explore the diverse milieu of interactions that comprise the connectome.\n",
    "\n",
    "This notebook will walk you through the following five steps to compute and visualize Basis Profile Curves (BPCs) in an interactive manner:\n",
    "1. Getting started with the Python environment and packages\n",
    "2. Loading and checking the BIDS data and metadata\n",
    "3. Preprocessing the data\n",
    "4. Calculating BPCs \n",
    "5. Visualizing BPCs\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Getting started with the Python environment and packages\n",
    "Before getting started, you need to have the necessary software installed and data downloaded. \n",
    "\n",
    "An iEEG dataset formatted according to the Brain Imaging Data Structure used this tutorial can be downloaded from OpenNeuro:\n",
    "https://openneuro.org/datasets/ds003708\n",
    "\n",
    "## 1.1 Setting up the environment\n",
    "This notebook requires the following Python packages to be installed: [bpc-requirements.txt](./bpc-requirements.txt) \n",
    "\n",
    "\n",
    "To easy setup an enviroment for this toturial use this file:  [bpc-dev-env.sh](./bpc-dev-env.sh) \\\n",
    "Run from terminal: <code>bash bpc-dev-env.sh</code> \\\n",
    "or from this notebook: <code>!bash bpc-dev-env.sh</code>\n",
    "\n",
    "You can run the following cell to setup your environment for this tutorial, after running, refresh the page and make sure the bpc kernel is selected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# creating the enviroment for this toturial and refreshing the notebook\n",
    "# !bash bpc-dev-env.sh"
   ]
  },
  {
   "attachments": {
    "image-2.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU4AAAEECAYAAABHkN1jAAAgAElEQVR4Ae2dafBVxZnwb74mmQ9TmeTDVKqcSlKZqiRvJiaZLJPKOm8qcZKIJnFXQhQcRBBBUBSRxSCbgCIiO8i+gyyyyqKIKIogCALucVxGJ7NoNpNJv/Vr3ufSt/9nvfeec8+59+mqc8/Sy+l+Tp/feXq5/VSMOpWASkAloBJIJYFKqtAaWCWgElAJqASMglMrgUpAJaASSCkBBWdKgWlwlYBKQCWg4NQ6oBJQCagEUkpAwZlSYBpcJaASUAkoOLUOqARUAiqBlBJQcKYUmAZXCagEVAIKTq0DKgGVgEogpQQUnCkFpsFVAioBlYCCU+uASkAloBJIKQEFZ0qBaXCVgEpAJaDg1DqgElAJqARSSkDBmVJgGlwloBJQCVTeeOMNo5vKQOuA1gGtA8nrQOUPf/iDcTcVXq3wXNnocW1dUXmoPDq1Dig4YzTuTq0YWm6FotaB8Dqg4FRw1rQ49GUJf1lUNiobqQMKTgWngtPrrpKXQ/cKyrA6oOBUcCo4FZxaB1LWAQWnglNfmpQvTZgWotc7R0NVcCo4FZwKTq0DKeuAglPBqS9NypdGNcvO0SzDnrWCU8Gp4FRwah1IWQeaDs6HH37YrF27tsu/kbi2ZMkS8+tf/7qLXxaT7g8cOBCYj7T3CvvixF3/r//6L7N69WpbXgn729/+1l57+umnayrqnj17zL59+8ypU6fMtm3bavwkLnviPfbYY6H+btgyH1POIDk899xz5oEHHogt/8mTJwPjl1kmmvdiablNB+fkyZNNr169auA4adIkc8EFF5gNGzbUXE8LsTThFy9e3CUfaeJL2HorLJC85JJLzLJly6ov+qFDh8yFF15obrjhhuo10r/sssvM/fffb9asWWN69+5d4+fe/6677jK/+tWvQv3dsGHH//mf/2m6d+9uXnjhhYbSCUu/GdenTZsWKIdVq1ZZmcbdI06OcfHVv1iQKuLzyBycEydOtNDcsmVLbtAEeq0GJw/7lltuMcOGDasCav78+fbFB55opIR56aWXLEwPHz5s3n33XQPYwipKM8D5H//xH/Z+x48fD71P2P3zuq7gVHDlVdfqvU9m4Hz99dfN+PHjLTQfeuihGmg+/vjj5rrrrrN+aFibNm2y/o888oi9/uCDD5qrrrrKLF++3Dbvx4wZY0aPHm0uuugic+mll5oFCxZU0wtLqwjgJP9onb///e8tpPr37281UMApTW60cMKgodI8dUFLU//KK6+0oBs0aJAZMmRIVeN87733zNSpU632iMYKlK+++mrb3Kcy7Nixw55zr8GDB5vnn3/e0ISV9Hr06GFWrlxp/ud//sdMmDChCnTySNh6K1Qz4iUFZ1AZub+rcVKXkNsdd9xhy4i2PX36dCvvZuRV0+hMyGcCTqA3btw4C0YqqTR72fPyAkC0J/oh582bZ8Pt37/fbN++3R5ffvnlZvbs2QYtbMaMGfYamithRIMlnai0igDOY8eOWejRd4kmCcRoIqOJIhdeurFjx1Zh6b7wO3futOGBG9rhrFmz7Lk01YEmwNy7d6+V07XXXmv9uScfE+4FNIgLNAj77//+74b+VPzWrVtntV36ngH3wYMHDX2IdCPceOONLQcngCOP7jZy5EibV+QWVka0dleOfLwoL3GPHj1qP06Ud86cOS0towK33MDNBJz0ZwLH66+/3u6PHDlShSfaImAEis8884zdeOnRegScaF4CW8BJeDRYrqENkT4aW1RaRQAnmiYvKVol2jQw4IWhr076MtH8eNG57r7wN998s33Z3RcMLR1wki4wWL9+ffXlR5ZcA5y33XabGT58uHnllVfshszw27Vrl/Gb6sAVP/L02muvmTfffNPwvNz75n2MxkmeBgwYULMhK+RJfqLK6MpRtP7//u//rpZp0aJF9kOSd7n0fuWGpfv8MgEn0EQjePHFFy0s+vTpY19KwCcaI/BzN7QcAeerr75aA06AISDFj3iPPvpoZFpFAKe84GiVU6ZMsfnlmoAMQAEIGahxX3ia3Wjj7sOSPk40R+I99dRTVf+3337bXgOcQBl/f1u6dGkXcNLknzt3rgFKhGe/efPmarru/fM6TtJUjyqjK0fAyYfZzTtaOmV1Yer663H7AC6rZ5kJON1RdTQtQMdLD/x4SV1/rj3xxBO2qSjgFEiyR+MMA2dUWkUBJ/2UwKhnz56G5rc8SLRPtErRQrnuvvB8SACuhGcvfZy/+93v7Iu/cePGqj9TeIAB4CRdmvJu3CeffNJqk77GCcTZCMtHiX5p0nnrrbdq4rtpZX2cBJxRZXTlCDjppqAPWfIt3RNyrnsFZdo6kDk4gR/9ecCTJjhTcjimmf3yyy/bJizndPSnBWdUWkUBJ/2bgIjtnXfeqb68dE1wDVDJQ3NfeKCIP32SjMAzK4Fz6eOk3xIgo7WeOHHCXHPNNdYfcKIx0qTlg0Rc5j5yDgw5Jx3kTb8rfX0MGNGspwuA5j/+aLCSr7z3ScAZVUZXjtLHed9991ltm35x5ObKPe/y6f3KD+pcwPlv//Zvpm/fvnZEnJecF5nRcYDJnmYsgA0CJ4MijPSKFuo21bkWllZRwAmM0Hj69etXAyLABaDcid7uC49WyQeHMGxopjQ5BZxAmAEP/ICiDB4hX+45c+bMaly0Xfc+jNATjzikI+eS1ooVK2rymveLDuSkD9i9N9q79HFGldGVI+AU2VE+NrTVVn4Y3DLpcTkh2nRwCuCS7BnFlUGfJOGjwjQzLfc+ra7YNDEZsPHzQbOfEXP8gSyj4kCBpriE5XpQXPF398RjcIg+T/d60Y/jyuj2cSILmT9b9HJp/ooN1JaC0wVUUY+LWoFlDif9dYyIo1XdfvvtpYJeHrJ1wZnH/fQexQZes56PgrOki3ygadIkBZY035nvWTZtsVmVOCodBs1a3fUQlT/1KydoE4GT5nSnblqxy1mx9bnpc8uyDgSCs1MhGVRuBiF0K5YMsnwhNG0FbpI60AWcQfDo5GsKzWJBM+x5JKnsGkah2Kw6YMHpVsakkGSKUSdsjNrqVgwZuPU06rhZL4emo6ANqwMVvwK64OwEMMaVUaFZDGiGPQe//sp5WIXX6wrDZtSBGnBSOYNAwvy+Tt0YvdateDIIAqlA09034yXRNBS2fh2w4HQrIeD0IYm5i07dWKZMt3xlwLSquM3/mLl12AUnx36l13MFYaN1oOJWOCqjQNMFJX9zjNpk+bJ23LPQr26tlUHQh8sFq0DUrcsKT4Vjo3CMil+RprkAk0Ul2Hbv3m031nDUTWVQtDog9VPqq9Rfqc/s3f56jov67zTN1xulezaVP//5z+ZPf/qT3d5//33z7LPPGnUqgTJJgDpL3ZV6TJ3+3//9X7v95S9/Me5WpnJpXosrgYpUNireH//4R7tMWXGzqzlTCXSVAEvrUXcVnl1lo1eykYAFp0CTNj0mGNSpBMokAeosdTcpPMtUNs1rMSVQcaFJhzqLA6tTCZRJAtRZGT134RnWZC9T2TSvxZRAhYrG15qKx+ikgrOYD0pzFS4B6ix1V+EZLiP1aa4EKliL3LdvnzUzi30gVk5XpxIokwSos9RdjLBRl6nTmJLGYCAbJkTYMEctGzaYmrkBb1bTj3OEIWwz750kraT5i8u/+p+WQHUCPPPimK+I5UR1KoEySYCpUswzZkoS05CY3sNq71gDxUQGG8BilXvZfvOb35hmbthzAk5xjjBox3k73u8k+cs7X2W9nwUnD5JJxphL5YusTiVQJgmgTWJ4jvpLPaY+Mxk+bMCIvk+3/1OmLjWy515ofnEuSZi4NOr1b+W9681zUePZfw6Jtok9FgVnUR9V/fm6//77ralhzA0n3YgT5ebPn584Le4Zl17UveL8qLPUXVpM1GUXmky3cyGZ1ZxOoJsETEnCxJW3Xv9W3rvePBc1XsXVNvlq0yekrr0kALjSurg4cf7+/dKG9+NHnVNnfY2TgaKgEXYFZ5Qk1S+pBCw4+VLzxabPh051de0lgXqgFRcnzt+XYNrwfvyoc+osdVe0TpQBGWGXSfFZa52qcUY9ofbzq0gznS82HeeMSGbhMN+7devWLJLWNGMkUA+04uLE+ftZShvejx91Tp2l7lKH/ea6gvOM5LSpfkYWjR5ZcNKpTqVj5DErcM6dO9ece+65ifI7cODAQkAWzeWcc86xI7aJMl7QQPVAKy5OnL8virTh/fhR59RZ6q4011vRz9lMjZOPwJYtW8yDDz6YeMOaZ5xTcMZJKLl/RUbTaepQ+ZgHl4VLA84rrrjC2grPIh9p0uQF/NKXvmSef/75NNEKF7YeaMXFifP3hZA2vB8/6pw6S92lDsvIet4DRM0E56BBg8zQoUPN+PHjE2984JmCFeUUnFHSSedX+dd//VfzhS98wfz4xz+2y8k9+uij6VJIGNoF5wMPPGArBC/TN7/5TfP973/frF692qYENIEV10aMGGGv0Xc1fPhwG/b88883jOjiWBWnd+/ehjxffPHFZtOmTTYOXQKEJyz3GjVqlA3PD023Cy64wLz00kuG9UM5njFjhtUsyQsVFs0Ff+KTFzTlVatWVdMo24EPLTQ0/6X0n7sfxy+z779s2bIuabKUmzg/vFxvxp68M1fTBaf0czJAVO/IOvbY+Xj6jmv4ua6Z4Lz66qvNqVOnzJIlS+wAF3vKtnHjRvPyyy/byf6HDx82bEz855rEcfPkHys4fYnUf15ZuXKlnYLUq1cv8/Of/9xCqP7kwmO64GRqCkC6++67zcmTJ82UKVPsOYCkyQGoJk6caOgXZRSUStG/f38b9uGHH7YABbRMQyEdvrb8e4TJzz179rT+N910k9m5c6edBgNcxfEiEef48eM2PY6BLpOD2Tjmfrx4VEr8gS9pl9UFQYt5h+vWrTNr1qyx8x39sgXFccP4/oCDNKdOnWpfaD5QrvPDu36NHrvgpA7JXE53ZD3t4BBg5KM6ePDgGngCTa7h58IzC3Bu2LDBQp89mjQT/amHAJB3g+322283Dz30kIKz0UqUMr5tqlPZ6E/5h3/4h9zACeyobDgqBYA6ceKEPXeb6jST8WPpMFmY9p577rGAE3C6U6gA55VXXlkVA5COA+fRo0er4dFiuR//RGn3pjoaDB+FIBcHujD/adOmmWPHjnVJMix8l4B1XBBw0jfYLHDy7Gkyu/B0oYkf5+KyAOeECRNsa4kPwb333ms1SzR7/t23bds2u3Xr1k3BKQ8hx30VnHy18gRnGMwouwtO0fqAmbvRrBZwokWKA5wzZ86U00QaJ/1h4tBSuA8wVnCKVLruw0DYKnDy98oocNJcB25sSedy+vAUTdOHJtLJApwAknqI9k6znBkDKBeUVczU8NGgGa9N9a51NMsrFZnDWVRw0hwBZORTHP8LRksUcMp19oDTXahkwYIF9pqEoQOd9Nym+gsvvCDe9quOP10ICs6qWLocFBmc1BW3j1OmJKUFJ4V24Yn2GQRNwmUBTmaXUI5f/OIXpnv37l02Wonz5s2zGqiCs0sVzfRCIcHZo0cP2+9Js4u+KQZpGOChSU9zHY0UrTIJOKlcgBAA80IxaOSDs2/fvvYrzqAQlZX7cV80UcLyUaECl9WFQS5pU53BB2z8uJp9WJpF0DibCU6eucAzDJqEyQKcMkqO4cQgY4k04amztJIUnPm+nS0B58KFC02fPn2qJXUHbLg4Z84cCyy0RxwaIYM2QIxtyJAhFmph4KQfSBzwo89T4jIgxbGrcTKSLv5oFe70I6CKn9v8l7TLsvchx6j6uHHjqv8159ifhubGobm4fPnymo+H648cli5dWpMmo/asUiTODy/Xm7Gnueo21ZsNTvIIPNnCXDPBySyP73znO+a73/1u4u28886reT5B+dRR9SCp1HetCk76h6h8VMKiOl4IV+tJk0/m+THq6zqa40ARuPJSoNG2o6sHWm4c+gcZgXed6+9eDztOGz4snaDreYAz6L7utWaC0023mccKzuZJs1TgbF6xT6fkgrPZaRcpvXqg5cZhpoE/Uu76Jylr2vBJ0pQwCk6RRPRewRktnzS+HQ1ONGzmHfpzDtMIsAxh64FWXJw4f18uacP78aPOFZxR0jnjp+A8I4tGjzoanI0Kryzx64FWXJw4f182acP78aPOFZxR0jnjp+A8I4tGjxScjUqwBPH5EwDgSrPFLTxcpIWMFZzJKqGCM5mckoRScCaRkoYptAQUnMkej4IzmZyShFJwJpGShim0BIoATuZUJjGGRphWzAlm1kiS/BX6QRcocxUWzWByM5O8t2/fbudQFih/mhWVQKwEmPdL3aUOU5ep00GmgrM0EQyUmPIW5whDWLS/PLek+YvLv/qfloBqnFoTSi+BImicpReiFiCVBBScqcSlgYsoAQVnEZ9Ke+dJwdnez7cjSqfg7IjHXKhCKjgL9TiyyYxOR3rfLgjM3yLZki4rl83T0FTbQQIKznZ4ijFlqGfyeVycOH8/S2nD+/GjzlXjjJKO+mUhgQrrWjLidvDgQbtMG6Y01LWXBOqBVlycOH9fgmnD+/GjzsX8C3WYukyd5v/1rIDFcoKsSYANH1a9YmO1LbYXX3xRN5VBXXXAghOzFFQ4VpzOyijZM888Y+3bYONGNuwLse6lumwlUA+04uLE+fslShvejx91Tp2l7lKHqcuAE2iyWjrQFFgKKDHEx8Y6o7qpDOqpAxVWvaGiATZAJtYmoypqPX4spsESbqzDycbCq5xjAuO1116rJ0kbh8nEjdg+56+D2HZpZ1cPtOLixPn78kwb3o8fdU6dpe5Sh6nL1GmgKVomwAwCpZif0P0rVVMcKotksqhQwfg6U+Gwa+KvuxhVYdP4Ac7LLrusJopAD796Hf+IAMDu4sNp0sLwGwsjt7MLghZrk+7Zs8dOGPfXKUUWQXFcGQX5M7kbI2JBf+0LCu+m18gxdZa662qaQJPmuAATIMgq6qyozsYHWzeVQT11oEIFA558pWnmrF+/vpE6HBo3CJwExoyvgAtTGUH20wnHCuRolkCSVdox4sZLEWT7HPvqch377LJghW+HndXd0XglzXbtNgiCFv+2mTx5st1mzZrV5bkFxXED+f4sMI3mPn36dLsS/I4dO9zgsSCuCZzyhDpL3XU1TenDpBkGNF1QYmKXTaym6v51lcXr6WTQBZxh5mJT1uUuwQEnwMOcAhtf/y1btlhoAboo++lUfuBGOOKJ+QtWhBcrmGL7HH/CYsSKeEuWLLHn2GYRUxtih33//v0W3HQb0EfGVJV2dD7kWNEdcxlomgCPY39lfT+OLxffn75EtHccAzJ+K8IP76fXyDnPvh5wKjzTwUI/MGfkVaGJS6WnuY5GhgGvLBwvEkDzN+wBASzygR8vgDwgsZ8OaPGbOHGibX5hnIr+LF52v6nO4sRbt261abLwAuUhLhAVcLp22DuxqY42NmnSpOpjvuuuu2wdqF6oo6nOR1BsPdEFAygBtLgswckzllF0GQyiX5NnzocUbdOFpHy82WMQTTeVQdo6UOElormOlkBTJ0twYnAN6LHxcjH/DqgBOtEcOXc3mtI4wmDdEj+uYU0R4PrgBKZopGiVhBNDbS44Xe2qE8HJYhfITxzNawy4uS4OdL4/xtrcuoM/dqzE+eHlejP23Je6Sx12+zbdJjofYwGm+5Jgalo3lUHaOlChn9DVOt3K34xKLWmE9XECQ+BFpQeKNL/Fif10XkA0ChwG1Zh+QlgGN3xw0hdKvyYaB+7tt9+2YV1wWo///9OJ4ESu7kyCO++802r5rlziQOf78+GTvmQGibBy6To/vOvX6DF1lhaTaJvUaYGmaJrAMujloIWim8ogbR2oapwyQJQ3OBkYGjFiRKT9dPofASUTnOkLRUvmnJfVt32OJnXuuedaAAPZO+64w4YFzNJUd19UwmOGmJeqXV0QtAAbH0w+MPRx+s6Nw0cnzq46gALADPCxxNvs2bNrknTTq/FowolonDIFifohzXRGTMmbaJt+8zwIpnpNNdC4OpAbOO+9994u05F4Z+hro1nNiDYVPsh+OrBkBBhYyjZq1KhqH5pr+5wXRkbfCUs8zhllDwIn/Z0yst6uRtuCoMXalVxnA3S+c+MksatO/Llz51bTpMnsOjc993ozjpOAU5vqCsM4GKbxr6BxuE11Rq5b7cLspzMKTOHoI41y9H2iesvgBPHoU+1UFwYtBs/YgpwbJ41ddbpV3D5kSdtNT641a0+ddZvq1Gm/qQ44pW/TfUHSNtE0vDbrqQPVUXWasvQjbtiwoVn1WdMpiATqgZYbh3rB4IvrXH/3ethx2vBh6QRdp86SR+pwUD+nTHCW2RphzXYBq+51lD2uDig4g97ENrtWD7Ti4sT5+yJMG96PH3Wu4FTQxYGu2f4VvtB8qdEosvznUFTFV79sJVAPtOLixPn7JUob3o8fde7+c0i0TvrLZXSduZxsonkyWCSbaKG6PzO5W2URLwsFZ9Qb2SZ+7b6QcRw46e9UcMbDQIGZXEYVvtB5LPLRJgzSYhRQAv4iH9RpdyI8My1kepL8k0gX+tDFPaQFUs++uqwcq8swT5K1MtWpBMokAeosdVdWSJJ/EclAEc12abrTfBeQskcb1U1lkLYOVBiNpG9TVoHPaj3OMr2ImtdySYA6K6u/U5ep07SiZEI80+1ksj9TlYCnbC5E9VgXNU5aBypUNNE2WUcxqxXgy/Uqam7LJAHqLHVXtE4fnmieNN0FoKKBsgekuqkM0tYBNdZWJkJoXgMloMbaAsWiFzOUgIIzQ+Fq0vlIQMGZj5z1LmckoOA8Iws9KqkEFJwlfXAlzraCs8QPT7N+WgIKTq0JeUtAwZm3xPV+TZeAgrPpItUEYySQGzjLbFedydKswMPyaknsErH6DpOyWeeS9T5Zzs51pMVf/rJ2rAjFknGsiE8+du3aZdpx6TwFZ9Y1SdP3JZAbOMXmUNnsqjONhXU9ZZ1Q1gGNcrzErO+JYboxY8aYHj162PjEY11RHP4swpylYz4aa5CyGj4LRQ8dOrRqSoRV8ZO4gQMHWpMlScLGhWlmWv69FJy+RPQ8awnkCs4y2lXHSiYrygM9DJIBvTDHivP4s6CFq5lihRH4il2frMFJXjFJArTJkzg0YWCKSeYkjjSaNa+3mWn5eVdw+hLR86wl0FJwUrii21UHdkCPSdSY+UBjDnOLFy+24PQXTQaiLLQh1jVdcKLR9u7d294D0LmaKE36JPbhCec6Weme/2z7DjtNlIe/mAF0V4OmGY+mzL9qAB3hRGMlPH4zZsyoGsJDi8XsMi5NWn6eGj1XcDYqQY2fVgK5grOsdtUFIpL/MCHTJB4wYECYd/W6gBP7PBzT/8i/WDA0J5BmMQqOk9iH9/tL0ZIBXpBjNXzS3blzp4U50BbHyu348XfFp59+2mramGSWpdrwo8uCv+eycYxNehwfhqRpyf2atVdwNkuSmk5SCeQKTl48fyuyXXXsIKGRATea6/3797e2kebPnx9o3A3/YcOGxcpewEkzGHtIgE+W9MJwHHaSWJ4/jX1496ZohWiqYY5nsHbt2kjYEddtXqNxE+/o0aPVZPlPONeAfBQ4/bSqCTTpQMHZJEFqMoklkCs40VDKZFeddR6BHGDD3AJaHF0LwIKlqHzHKDrN7SC3cOFCQ3o4ASej3aTlb4MGDbJ2e9LYh3fvyX1IM8j2D0DGD9vqaWAn4MSqqDieJWnRBZEmLYnfrL2Cs1mS1HSSSiBXcPqDQ2QSraaodtUBYb9+/aqyZPEIQBGmzaEh4k+/peswiAYsZaBFwLls2bIulj9pstOfmNY+vHs/QE8+6Hf0Hdoo9ydPCxYssKaRJQzmBYhHUx0XpHGSP3GM3BMeqKZJS+I3a6/gbJYkNZ2kEmg5OItsV33lypUWDNh1R9NiEAZQsMkIuStomvZo1TS/0ejQyNBMacKjrcoIt4ATP9JiWTTCStOXlX7o90xjH97NB8dTpkyxaQNP+lLRNIEb96NcuAcffNCe04eJZdHhw4fbcwEn2jPpEF80TkwxkxaDQkwx4iNCudOk5ee10XMFZ6MS1PhpJZAbOMtoV52mLn2cAkvgBzBlTqqrfYnggSMfA4nDHs0N8IgTcHLOpHTSJRx7NEIc2lwa+/CStuyZkrRo0SKrXUpeSB/AieNjQB+z+EvXgYBzzpw51o9+VwEnI+kSnsEylmrDpUlL7t+svYKzWZLUdJJKIDdwJs0Q4YpmVx0oYC88jWNqD/ATLTMuLum7cz8JzznaXaP24WmCv/POO6FZwI/R9ign4EQW7733Xmi5kqQVdZ96/BSc9UhN4zQigUKCs5ECadxsJOCCM5s71J+qgrN+2WnM+iSg4KxPbh0XC82XLooi/tddwdlx1bHlBVZwtvwRaAYalYCCs1EJavy0ElBwppWYhi+cBBSchXskbZ8hBWfbP+L2L6CCs/2fcdFKqOAs2hPR/KSWgIIztcg0QoMSUHA2KECN3noJKDhb/ww6LQcKzk574m1YXgVnGz7UghdJwVnwB6TZi5eAgjNeRhqiuRJQcDZXnppaCySg4GyB0Dv8lgrODq8A7VB8BWc7PMVylUHBWa7npbkNkICCM0AoeilTCSg4MxWvJp6HBBSceUhZ7+FKIDdwdoJddVaJX7NmTc2GcTRWFMrSsYAyNoR8J/bbWaGJfAWtCO/HKeO5grOMT63cec4NnLKGZTvbVecFZq1KDJhJOWWtzYceeiiTmiJG3bivv/SdrPuJ1Ur8ky5xl0lGM0xUwZmhcDXpQAnkCk7fdAZmdFmsF6jW61gbEijIgrpp08FsBwsPh7k0dtUFnC6gWB395ptvtosZh92jkessfAyoWS2ehYtdp+D8vdWyWdGJNU1Z35SNRZ5lc+WlxyqBpBJoKTjJZDvZVQ8CJ2XEVAUQ42WNsj+O3/jx483YsWNteLRVzNLrkIcAACAASURBVGqEOaBMGIyzzZ4925qx4B7iFJwKTqkLum+uBHIFp9glpy+QJuaWLVustkhfHC88mhP2eVg09+GHH7bwABxiFIxwxBMTD6wU/8gjj9g0gA5GyqTpmsQe+f79+y24uS92hfwV2EXUmL5Aq5X8y3V/L+Ak/5SR/Bw4cMACDcNvuChrkPhxH8pHGmI3CJs/QU7Kjhww8EZc7KGLU3AqOKUu6L65EsgVnLzY/tZOdtUFnH4Z6Y7AuBkuDpyEFYDT5CetEydOBD51zAjTDSAOwI8YMUJO7YcHuGofpzbVq5VCD5oigVzB2e521QWcrJYu9uMx6AYMx4wZYx9YHDh79+5dfbCMggNOMZ5W9TDGvP32210+QgJsNFCcapyqcbp1Ro+bJ4FcwekPDlEMtKR2sasu4HQHhygjg1+UExdlfzwKqjay80M69G9i2le2w4cPW5hKv6iCU8HpVBk9bKIEWg7OdrKrHgZOQAbkcFH2x5OCk/5gRtHpA/UdzXf5QCk4FZx+/dDz5kggN3B2gl11bK7TXH733Xdrng6T4LnOAE6U/fGFCxfa+Z8SOaypLppl0BQs7LRzLwaXBJwyuCZNeEm/XfZ8sOgeYR4rA2mUk6ludJcgQ52O1C5PujjlyA2caYpMxQ/6lwu2v9966y37QkSlx+BKM+yRyz3qsasuccP2rbA/HpaXsl9XcJb9CZYv/4UEZ/nEqDlupQQUnK2UfmfeW8HZmc+9rUqt4Gyrx1mKwig4S/GYNJNRElBwRklH/bKQgIIzC6lqmrlKQMGZq7j1ZsYYBadWg9JLQMFZ+kdYugIoOEv3yDTDvgSKBM7vfe975tOf/rTdPvOZz1T/auvnWc/LI4EdO3ZUnynPdvTo0apxlufxaU7DJJAXOJk/Cxhl+9WvftUlSx//+MdNpVKpbkyfK4O74YYbquWifMw5VndaAiwg5D7THj16KDi1cpRfAnmBc+/evTUv0I9+9KMuwisrOPlnmwuHJ554okvZOvWCgrNTn3ybl7tI4GQ5w4kTJ1Y3/r1UBqfgDH9KCs5w2ahPiSUg4GQlKv6RlfYvl8CN1aaC/q3miiWJxumGT3LMv9x+85vfJAkaGoa1C5LkPzQBY+xaCvVonFhgKPLHgfyllS/1B4sB4hScIgndt5UE/umf/sl8+MMfrjY1P/WpT5kBAwaYF198MfC/6vx3HeN13bp1M3/1V39VjQc4OGclK9dGFAuzsDQggz0uXAjLdbaNGzdamdL/JddoyvsrZRGIFxnLB9/4xjeq6ZEWfYsMPPA/e9+tWLGimi7pHzp0qLpItlsG0mS9gqTuwgsvtOm65eL4i1/8or2Ov+8oK4t6f+ITn6jmn0GTyy+/3LBeQ1rnl41z382cObOm/BgoFIecReY33nijhd7IkSNtGaRcH/vYx2ydCJIt6bBm7TXXXFPzjH/wgx+YDRs2WKsNkg5728d57bXX2q80CyTw/26+3upUAmWSgFup/WNWpnIX+UA7AgZ+uKBzmtw4zJkE+bvXsP2Ei+vj5EXkJXbj+scAyX8PsSDghhs4cGDNuevHMbaykjg/nn8OlMUBfD4qfhj/nLylsezqlw15+47BK/c+dImIc6/zcQPq7jX3GH9fQ6Yl4X583PAc80Fzr1lwcoFV2BWc8hh0XzYJuJU66Hjfvn1VY20sY+iHAVRon0FAQ0NJAk7WR8VFgRMT2f69w855kTG9Is6HS1g8uU6Z3OampOPvJXzYHpmIYyHysHD+9SgDiJKe7P2yNQJOPx9B56xCJo5ZD0Fhoq5VwUmggwcPZqpxltmuughZ98WUwAc/+EG7QDQmRjCMRzPTrfhUdLFy6TYvCTNr1qxqoejb8v2Zw0dcVubCdr2bLs1DrrOJkbwwcGJY72tf+1pNfM6xe4UmR9r+vd1msg8X8sH90ZaYJnXdddfVpI0/73Sck/z7WhUaL34C382bN3dJnwW6X3vtNWth9rbbbuviT76SOL9szQAnLYCjR4/a9W/9DyIfSXFBH9Jly5ZZFh45ciSwdVIDTtr3WTbVs7KrLgIo0n7+/PlmwoQJRcpSW+eFBZ3d9Tjpn3QBx4vDwA9rdf785z+3TTmaczQ7/QEh1jB147r9bUkGh8LACSDddNEo/f7PX//61zVhCA+YcD5ciA/YxAFuH7z0zSZ1caPqQNrN/5133tkl6X79+tWEgSlJnF+2RsHpmxun28LNOx9WHANz7nWO161bV5NlPhxBH+KKRPzsZz+bOThlZXLJGc0gHohfUPEv6z7OVntZy1XUfMuouruQMQMkWEVFU2KjX4sXhQ3IsHH8+uuv2wENmtpBfWPNAid9cvKusR88eHCgOP08oPHifLhcddVVXeL7/Y/SfdAlYMCFOHD6WhuzF3znd0UwUJXE+WVrFJws3O06FvV2ZS/dD3yU/Ou0DHw3ffr0mnA1Guff/M3f5A5OMpjUrjovBRYd0QjYRo0aVR19xILk8OHD7ZQKKgB+8jVnKgEPIshWOQKlf5fmHaYoWDmdUU0RHhoK6XK/888/36BJigvLDwMKhCctXgJJS+LpvvkSCAInH2VgiUbpDg6hQQDUoBF19yWS42aBE1BKmuzdwQ1XIowGu+FEqfDh4lozlfjMJHDjht1Dwrv7KHDSleCmK+Bx43OMrN1wHEsXhh/WPffL1ig4/cEfpmq5+ZL8Y4LcvY4MgpzfRVMDThJ48803u4zmBSVUzzUqACDB3nhau+oIn8yy8VVjygOFBFL4AT+0WfpUsCvOffr06WOziR0fIBZkq5y08Ovbt681eMYUB84xqUu6YXbeo/LDXEI+BnG22uuRocYJlkAacF5//fU1L4v74tDM9kdXmwVO/75hUMMaqpsnPvi4RuESLLkzV6PA6Q+gCHjOxD59hJLg5p1j6SP1w7rnftk4953fDeDKz7+nHxclxw0j+cekjXud1m+Q87t+LDiffPJJw9+rsh4cApxAyd+S2FVHMyQe8/LEAcnt27db2zr4EUacAJH/2wJOBEKzDIcmSngGEiSc9CPhD4CXL19uO7wJRwcxzTk2muCMLEblhzS0qW5FndtPEDgZOeVjSf8nzwOtye9n5KVhcIDrYve+V69eNS9Ts8BJ/XdfUjEX7QsJ89BuOFpMOB8uabUy/z7+eRQ4URTcPHHs9w2Tnt9H+/nPf96/TeC5X7Zbb721Szi/C6MZ4IQPbrnoIw5yixYtqglnwUknLwMZwCLrwaF67aqjUgMxgZ9bOLRD/NwmMXPIuPb4449bcIbZKhdwunF79uxp58BJuqTjbjTDo/JD3hSc7hPK/nj9+vU1g0N8VN0XAi2S5hvPxb3OSLTv/NHlKHBSF3wXNjhEX6V7b+YT+o7BHl/jxfQzzodL1uD0J9G7k/Upx9q1a/3s24+UW0b6XJM4Jre78egWcx1dLWiJbphmgDNIQ+bd9h3dOu69LTjlAn0rWYPTHxwigwiXCk0FAVBMCRFHE4EpBQI48icOLYMvwbPPPmvjuVojzWXSwgpklMldSdcFsoCz3vyQPwWnPKV89meffbZtjdAk45ny4km9Zt+9e3fbZPT7D+kTdF2QRhoFTtJmpJ3WCP3huDBwUq/9l59J3dKUBez+xHw0NvHPG5x0NfH+yHs1efLkGplSFt5NcX6zF9nw54Mkjg+f+7w4FnDzMQma7N8McJI3yunem+d36tQpm200bazzuv4cV8HJ39X4y1ErwElTiY5u6E+Fl4EdKiNQ5WuE8GhKDBs2zDa1gRqDOatWrbIVi3j9+/evTuKnfxPtFiDWC85684PEp02bZgBwWZYUS1K5ixxGKvbf/u3fdqnk+AFEAISWJGFlT//2XXfdZRillmvu3n1B6Spy/dzjJP8cAiRuHDnm74py7O5pMYnLGpxh5UcDxiE/QO7mj2Py7n8suP6Tn/xEsh67Z6K/n27cuftc/LD+DcP6OAkHJP34nFMmX/uXcFVwbtmypQodNLksXCN21ckPmiWARJOkiUQfkcxjA/pAVprU9Ju++uqrthhRtsr5YhInSOMkMporAJZ0gbz8lSwqP1R4GVmnmaEuWwlIhQ7aM8FdRtXRCsMgFRSXa25fJBqIP1dS4iUBJ1Lgv9QSJ2pP36zrsgbn0qVLA/Ml4CQv9PeHld8tC/Me6e9M43hn3TT8Y38uZbPASR79rgL/3v6HwYKTv6NRoaBylhpnGiHSrAnqfCYN8ur2Sbrpvvvuu9UpSu71Ro/rzU+j99X4ySTgV3TOpS/an47EoCCLNwTFYSoa/xpx/fw1N/mYok25YTgWcPpgCZrviHLig0DSo4+Vj7Lv/KZyUB8n0/UkHfYuXPz0gs7RvH1IuOAkDtO8mDUSpI3RfGe8RLoXgu4Rdg0lyJ9ORRm4D90lPtjdsrll5th3DAi7Ycin7xho9hdxIQ6DUn6fuQUnUCgaOP1C6blKIEoCdOgDIwY46WujtRE2j5PWBZojk6R3795t49FcS9syIDyKBn2AzBN0Wy1ReXX9aL3wUvIPn+PHj6fOg5tWM49hAs1nusuQY5gjDH28TA9EBs1wQI6BWVZg4rnwrPJ0MvuCVmNU2dXmUJ5PRe+ViQSAJhCj1YQSwIsfB05eyLxfykwKr4m2RAIKzpaIXW/aTAkoOJspTU0riQQUnEmkpGEKLQEFZ6EfT1tmTsHZlo+1swql4Oys512E0io4i/AUNA8NSUDB2ZD4NHIdElBw1iE0jVIsCSg4i/U8OiE3Cs5OeMptXkYFZ5s/4AIWT8FZwIeiWUonAQVnOnlp6MYloOBsXIaaQosloOBs8QPowNsrODvwobdbkRWc7fZEi18eBWfxn5HmMEYCCs4YAal30yVQwUQp/wst0iIfTS+lJtjWElBwtvXjLWThqlYuWVop69WRsGnE+ph33HGHXbOSBQ7EscZm0OrL4t/oHrshLD+nrv0koOBsv2da9BJVwfl3f/d3djUUKmEWTlZqZ4FfjKyxGjzrXLKiDc5fcDhNHpLYMWd5KFaaVtd+ElBwtt8zLXqJquBk7Tm0sqzAyert2P5xV6QR07tBK7WnEVwSUxUKzjQSLVdYBWe5nlc75LYGnCwWmhU4Me/JKs/uuoes54fGyWKzaJys4j506FCriWIqQ2yWYFUS8LlrHgJd1jEMsmNOepMmTbLmNtBqWcUbUxuicZIe5jXw4/rcuXPt4sisich96LLAYf2Tc7FBwhqOnLN0Gd0N2DzC9ALpsDq92/XQDpWjLGVQcJblSbVPPmvAiXmLrMApliEx1QvssPonJlkRJ+AEQFH2z11wilG1IDvmpA94KQswHDRokE0bcLJQKauDowEDUFZ+5hx4snI1x5s2bbJPGE2WPC1ZssSeA3kxOMeHAD8G1rgHABW/9qke5SiJgrMcz6mdclkDTuycZAVOhMbK3HPmzLFNdqDDxj1pvgPOOPvnQeAkXbepTlqk65ovfe+99+w1wMky/EDVNc3BvbmGw1gcUMWhAQ8ePNgu6c85BuFkyX7ASVhxYk5YznWfnwQUnPnJWu90WgK5gBMbQTR3MYEqjuPFixdboKGNAq84++dJwMnS+4Dz0KFDciu7p4kNOAE1JkFdx9L/xMHuCQBE65R0sFGDn8AXDRcHOPkIiDtw4IANJ+e6z08CCs78ZK13Oi2B3MAJfKQJ7Aqf6zSFo8CJdT3Cuf2j9E1KE9rVOIErYen/FId2CQwBJ7B2NVvC0JeKPw54Ep/pWfTJ4giPtUT6McUpOEUSrd8rOFv/DDotB7mAE6EymENzGCNIABCbMKJx0kcYBU6MRgGzNWvWWK2V5jbnAk7fjjlNaMCKdoimSL8p4QEnZks5xiohWi/zR8kX/aLiZIAKeOJIjzjTp0+XIKpxViXR+gMFZ+ufQaflIDdwAkoBEBBiA1hbt261Mo+yf04AwCbxGBGXwSH8fDvmNLNlQIg4hOVegBNHc5xzSW/EiBE1Fu2YloXfk08+acNjd57zY8eO2XN+3P5OzrWpXhVN7gcKztxF3vE3zA2cImm0PKb1oEWG2UeXsP4e64VAManjXsQJc/zN1G3+h4XT68WWgIKz2M+nHXOXOzjbUYhaptZKQMHZWvl34t0VnJ341NuszArONnugJSiOLitXgoekWYyWgIIzWj7q23wJKDibL1NNMWcJKDhzFrjezig4tRKUXgIKztI/wtIVQMFZukemGfYloOD0JaLnWUtAwZm1hDX9zCWg4MxcxHoDTwIKTk8gelo+CSg4y/fMyp5jBWfZn6Dm367oxRqq/KGBtVL50wP/VOMPEKxTwJ8cWDKQdQzYWEFLNhWfSqAeCSg465GaximUBFTjLNTj6IjMKDg74jG3dyF9cPK3XBZ3+d3vfmdXuxKtk7/4sonm6WufooXq/oxGrrI4Iwv3LVJwutLQ41JKQMF55uVW0OUjCwVnKVGhmUYCAom9e/eat956y7zzzjvWHItqnPnAQ+TfiXsFpzKoVBIIeknDwCkDRCxOTXOdASI2abL7zXa3Ca/HpwfSVA61g4lS/xScpcJGZ2dWKq2/F3C+/fbb1hgfI+vvvvtudWRdwakQbPYHIBCcmO3VTWVQtDrw2muvGX9jRX8Wpn7zzTeNgBPrqUxJkgEipiUBT6YlydSkIO3T1UT1+PRAmsqhVg4C4EBwdrZeo6UvkgR87VLOpQKzB5xvvPGGkbmcmICmn1O0TkbXZU6naJ803WW0XYCq+9Mflk6Wg3xQw/byIVFwFokSmpcaCQgk3b0LTKnEWEnFooAMELkT4YEnmqf0d4r2CUB9iApMdX/6o6JyOP3nCfdDIkBVcNa8qnpSJAm4wORYoCnApBIDyaVLl1qtk2lJmHrGBhU2oLAZdfDgQfP0009bc9GHDx82/vbMM88Y3VQGfh3Asq67HT161Mj27LPP6rJyRQKF5uWMBJJAE3BipZQmuro4CfyF+VtxgdrDP4dyqsbZHlWlrUqRBJo0n2hKogWoUwnkLQEFZ94S1/vFSiAInG7zXKBJHyVNLHUqgbwloODMW+J6v1gJuOB0+zVpmrvQZKCHPkt1KoG8JaDgzFvier9ICURBE3DSPEfTBJpMMzp06FBkeuqpEshCAoHgZOqGbiqDVtQBpg7JxlQiNiazMy+TfwQxuZ15mvwvnQnv+/bty+K90DRVApESUHDqR6JQH0mBJnuBpoBToMkUJKDJvE2mIKlTCeQtAQWngrOQ4AzTNoEm2iZTkPibJf9TV6cSyFsClf79+5uPfvSjplKpmLPOOsuMHDmyUC9SK5qLes/WdFP42qZonNJMp4nuapv8zZK/W6pTCeQtgQrAvO2228zixYvN97//fQtQBUdrwNHpcvfBSRPdb6a72iZ/s+TvlupUAnlLoDJo0CDb6c6XfNWqVQpObbq3pMXhQtPt3xRtk/5NX9tk9aY9e/bk/c7o/VQCpkIf0ZgxY8wll1xiPvzhD2cOzqeeesrMnDnTDB061EyZMsXwv0/RtPDbvHlz9VyuF2W/YsUKc+rUqZblj3/JrF+/vmX3z/I5pAUn2iZLzO3evTuT13jnzp1mzZo11W3r1q12ZL+em5FOHn8LZaoW07Nk0ZN68qpxkknANtXPO+88M3z4cDNq1KhMwQkoP/3pT5vLL7/c3HLLLeZ73/uePWeRBl7KadOmWb8sX9BG0j777LPNli1bWgauRYsWmW7durXs/o3ILi7uCy+8YF588UU7FSloYMjt32RQCHCyFueuXbuS1fSUoa644gpzzjnnmD59+tjtm9/8pvnSl75k5syZE5vSwIEDDaAVRzwWGsnSoZlffPHFhnySb+5JC1JdNhKo0L/J/DiaQXwZ6fOMq+T1+PMVBJrcw40/YMAAe52XQcEZ3bfazuCcNGmS+da3vmVOnDhhpyEFDQy5/Zs007MG591331196/gHE/UTMDERP8oBXRdaeYBz7ty55vzzz7fvFnmbOHGizSt/KFDXfAlUunfvbp5//nnz0EMPmc985jOZgbNXr16mR48eNdAEoK+88oqZPHmyXeWGisnDZ6QfyH7ta18z999/fzUOTWUGsMSP8KRB8+TCCy808+fPr2qxN998sx1YwP/VV181ffv2NWiMxL/nnnuqlQx/mr8/+tGPbLq//OUvzfHjx6v3xF82V+PknuQVDQR/ztGkydvPfvazqmb6xBNP2Ovbtm2z96AMaNvTp0+38iA8eWIit9wnLD/tDk5k4cJTBob4sMukd5mGBDh5rllqnC44efV4xkCQjz/aHfkT98ADDxg0TaBJGJ7piBEjrDfnfBjOPfdc64cWK013IEzXFX5AmXqK9o0jzfHjx5uxY8daP9JcvXq13LJmTxorV66sXkPj5b5xkK9G0INUErBNdbRMtjvuuMN88IMfrL7A8iI3Yw8EebhRaQFCXh7yQSWl75VzmmUsH8bx1KlTbb/orFmz7Dkv0P79++0xYKaflIpNWJrV9J3RvL300kvt+owADADiT15EE6ZC8r/nG2+80QKbwQg/rwJO1nrkmPyRPvnjfMKECXYNP2BK+pSBPjiOeSlmzJhh+0iBM9eWL19u80S+pQkelZ9OACdyAZ48bwGnTHwXjZNpSAJO+iKzcADQBSf3o1+eZjB/9eR58oETR3i0PprkQBCNjzLgBKSsE0r9JA1J+7777rNpobicPHnS3HTTTfacpjdKA3EJix9dXZzzIQlzQFdmyCxcuDAsmF5vUAIVmsg8YP6JwVdwx44dXYDhAyTtOXDhheCBRsWVphDNNMLxghCPyvjyyy9bIOLHP0aADn4MLgk4aeZJ+oBIKjLhKKP4AXCucc4Xvl+/flbrRvNG28QvaJAKON5+++3WH2hKelRwXiQGjkiDDY2EsAJONCMJDzgBtJwD+CT56RRwIgs+nlHgZGAIjTNLcAIpd6OFsX37dvvKAcbevXvbY/JBON4lXFBTfdOmTdaPH0DIbBYc9Ya6LE6gjLZJvQKydBPggCn3oZ6HOea1Al/CAVpWlVLXfAnk9s8hmhloZAIL2fNy0HSlqQo4ae6KH5oGLxEreXPMy0RFA2A0zfFzwUlaEhdwoZUCQMIBb/FzQQVg8fc3aYJLHPbcl3AMavESyf3Ilx+fc14sAaerwQJOXh5JW8JwHpWfTgFnz549bdM8CJx8NPmgAk66ebIEJxomH1w27uU6Po4CS7qI+PiKCwInK9GL48NNs57ykQZaqOtogdA6AZwCZ/wZNSc8H/c499JLL9mwunpUnKTq888NnFRCmuv0VQkw2DMdCsjwgKPACciIz5Qc4lGRieeCUzRV/AWcoo2604ikS4BwV199tW1ycywb/39m4EHOZQ84qcy8SNybdPAjb9LUlrCsE0kFd6EofoATbUDO3TBR+ekEcNIXTh1xB4ekqS7/T88LnNKcDnu1rrzySkNz+IILLrBjBBIuCJzuqLqAE00SELraKBoiSgaDS2nAyYd8wYIFkgW7R8nYuHFjzTU9aY4EcgMnX2hgg0ZJ3x8vB1BDewMkaIRR4ERbJSzNIV4ctFTS42stcAwCJ80b4vHlpolDk1k0R8BFfyjnNHGYWcDUKM7ptxSwyZ7rMh1J+lgBJOmSFyo6aTAgxDmTs10oSjpR4IzKTzuD89577zVAE0jyHH1w+pPfReOkbzALB/ziwAnwgBMb2qA4NEY+jNIXCRyDwEn40aNH2+Y4fZg006UfE806DTiZSkjfKu8DS+6hBXNfurjUNV8CuYETaKBV+k1RpiMJpBhpZpRfAMNLBIBoqqNZUkE5Z6P/kHO+znS6cy0InKRFXGnaE4emNRDEjzj0V0m6P/zhD0MnmbvgJB5p8qXnmC87GjHpsCdN0g8C51VXXWUHuaScbpio/LQzOGm2UnY+oOzZ3KZ6K8AJxKIcgAJOzApxHS0QrtPlgOPYbTKLxokfZaQ1Rhg2+jQxOIdDm2UEXlxUU513hdkokg7vhTuXVNLQfXMkkCs4BRQMQgEzXga5lmTPy8QXFC2S8GitNN+i4vLVR4tEQ+GYsGg3AM+NR9qEca/Vc0wTn7TqievGaVZ+3DSLfAwwZaPsbK0EZ5LXi/oLqOgGatQBRepzow6Yky+dv9moJKPjtwSceb/AaLn0edIMZjoTWuGSJUsahlve5Wjn+wk0i6JxRr82p+dYyr+K4sKqf/tJoCPASVcAzS6mHt16662BU43aGUplKFvZwMmo97hx46oT2dsPDVqiKAl0BDjLAI5Oz2PZwBn1Uqlf+0tAwelMQ+p0eLWy/ArO9odNO5VQwangLERfr4KznbDS/mVRcCo4FZzt/55rCZssgQrTPlg8gblkn/zkJ81f//VfF+JFamWzUe995l9UeclCNc4mv9maXKYSsOC8/vrr7ervLF/FqLM6lUDeEmDeoWz8FZGNvx++//779l85zE8E4nzomY/LXEXmA6t54LyflN4PCVSYsM2ScrNnz7aTybUiasVohQQEmuwVnK14AnrPNBKosLoM4OT/3PzXlb06lUDeElBw5i1xvV8jEqg8+OCDFpwXXXSRueaaa+yq5I0kqHFVAvVIoF5wPvLI3npup3FUAg1JoMLKPmica9eutYNErOzDYJE6lUCeEqgXnKxWpU4lkLcEKoxm/uAHPzCf+tSnDGYlWKmHNSE/9rGP2dWuOf/Qhz5kVwKaN2+evcYq2F//+tfNBz7wAbvaNauwsOo052eddZZdlYVr/M2Ra/xXnJV97rzzTvPRj37UrkyEP/8bx5+V4TlnI+xXv/pVe8xLwSIKjPRj2uLb3/62+fjHP27t9XDOKtrEZ3Ua/nv+5S9/2V6TtNw9K9UQljiEveGGG+w5pkIIt27dOms25HOf+5xdAJkVblgNCX/uRRziUybCf/e73zU//elP7bF7H45HTp5nvnHpbTV+V98yxVwy4E57bf6S1ebz5w0xY6cuMCvWrLfXrhg0yfz02nFm7qJVdum6njffbcOwX7JynZmzaJX5TOOiJAAAB5VJREFU5x4jzb/0Gm22bNlqCM85aa1bv8mMvnu+DT/j/pVm5ZoN5oLrxttrmzZvMV+54FYzbe6yan6GT5xruvUZWz0nz+s3brLhyMPCZWvNgmVrbBjyOWLSXBv2ov4TbLnunbvMrFi93vQbfq+9533zllt/8jN03Cx7PHX2UpsH8k36cXG5z7ipC2zZee43jZlpPn/eTWbkpLlm5v3LzS8GTzKf73ajubj/eLN+w0bz5Z8NMbMXrEhkdTLvl0rv1/4SqDBSyVqZwBPNk+0rX/mKNfyEnRUGj1ibkOtABZsmbJiO4BowZXUY6Svl/+CsLSgbUPvIRz5STfs73/mOHQnFn0U3SAM7OxKe9ToJI+fAnDCXXXaZhfoXvvCFalpcZ6EFtGZeTs5ZiVviuntWY2KdRMKw8TFg6S32Eg54/v3f/301DB8TFpTFn0Ez4sk5ZfriF79YjStpsJ+5fIf5v9dMrfEbctcq02f0EnvtyLPHzYU3zTZf6T7R3H7fA/bavNW77DnXDh85ZgZPXmnDXDxkTvX6pTfPNY88fsiG33fgsOGc8LKNnrnBnDhx0ux+7KC9Nnb2RrN0w8Pmx9fdZ547caKan/NvmGnmr9lVPZe8b979hOk2YHo1vX5jl5pv9bzbTFu6zYY98PQRWwa5H36zV562lUMaxL17wWYbdtaKHTadtdsejY1L/flK90lmwdrdNizP8/jx58yo+9bZ+//jFXeaniPvN/3HLTHXjl5kFq3daf6l71TzzJEjmS3Uiz0hMYXR/hjQEqaVgAUnWidLtbEYBhBlqTYqDecCT+yqsKQbS/KzYQObTUAat8fAGfCKC5fEHzMEmNpghfUk4d0wvJTEpZzudfeYtRNZeNa9lvfxTXevNteOWWbzcPCZY+bIsycC8/PM0efMgUNHzYmTpwL9nzz0rHnsqSNVv+eff8Hs2PtUaHjKeejIcXPsuZPVOH7Znztx0pAn/3qS86i4PBM2QLp+x36zdc+T9vjY8efsx/HCG2ebkdPWmscOHDJ79p1eo9VdPT1t5Y8KT0tH1sWMCqd+nSmBCqtOA07WPmR+HGsCijVBF56sVYn2CUAxWwFE3U2A2un7PY8fNl/tcVdDGzIcMmWt6Tt2hf1INZpe2eLzQR49a6P5dq97zOL1D5vNu580t05da7XSbQ8fsKZL0EL/8YoJ5v90uzGTN1fBmYlY2ybRijuxGK1T7Lu48MRUhWifAlAgKiAFpro1VwbTlu00kxds7yi58kGWj/ILL7xkJszbYn42eI4557rppufti826bY9ZjVS6Fmg9ZDU4BDhZNo51XDmm35/74bgnlgswV4FFAVZtx8yvWJRk0j6LZXOduHRfoWSoax8JVP7whz9YWyfyrwwfnmIgC+2TTSAqIKUpr5vKoBl1gI9y0IdZgCpdRHQJSJOevu0sHMADitg0Apj8JZlz3g/M+eI/bNgw60cY/PgTCQ5o0n9O3z02qbB3hUE3de0jAQtOtE5psvOXNiqHNNv5axvaJwDlL24YkZJNYKr70x8VlUNjcuCj7H6Y5ePstm6AKACV/lQGD7NwgBGbP+JQMIAjg0aAk2OuiWPWCLDE4cfApzjKsWLFCvOnP/1JLum+5BKoYOuECuDDU/o8pelO36dAVEAKTHVTGTSrDvBhdj/O8iECPC5EBZ70hTI1LgsHOP3BIbF3Djg5dh1hiYMs2LvG2dxwetweEqjIIgouPGm2y4CRGMxCAwWiLkiBadAGYHVTGSStA1KHJDwfZvk4C0iBKNooXQI056UvNCvzwMCP+bvi+P88GiWaI+DkWPo0CcNUNjRNrhHXBTr9/1g9dU0IS7q6L6cEKjQfwuApABWLgwJRmvL+Blh1Uxk0Ugf4KPsfaGCKFicAlaa8wJP5w1k44MfgznPPPWe7scTkL+CWPk76MlEw0HxZ5wGz07ghQ4bYPk26GAA+f7rgmrr2kUCFL6QLT7/pTt8nABWIUlHYgKluKoNm1QE+yrLJR1kgjEYqWij968AIgAGmXbt2ZfI2oj0OHDjQao9AlHPRbgEnUKW5jh8bYck3DtjzRwzxA5xoyOraRwIWnC48Xe2T5rs04aUPFJDKJkDVff4L/7abzOWD7H+UgSkABZ6ifQJPNE+0Tv4SnKVDkUDj5b/04tw+TvLEuxHkROkI8tNr5ZZARRaMFXiK9glABaKihQpIZU+F0U1l0Iw6IB9jgY18GACpD09ptjNgtGfPntzfQBecud9cb1gICVhwugvHugD1ISowdYEKVHVTGTRaB+RjHNTC8eGJBihaJ3Ml83bMzWSQSF3nSqBCE0TAKXvgKRvwjNtcoOrxaU1d5ZBODkHgBaKigdKXSh8iA0j0d6J10tf5yCOPdO7bqyVvmQQsOGUtRAGnuxeA6v7Mx0Rl0XxZuB9n+egAU7oBgKerdbp9nXv36kLGLaNHB9+4BpwC0CAt1IWpHp82JqZyaJ4c3I+RQBSAStOdPk/ROhkokua6grOD6dXCooeC04WoHCsomgcKlWWtLJOAU7ROaa7Tz6ngbCE9OvjWFSm7wFH3Z8zUqizyk4X7IRGIonnSXJe+TgWnvK26b7UE/h/wfku9N65IDgAAAABJRU5ErkJggg=="
    },
    "image.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAABSCAYAAADtsLBxAAALgUlEQVR4Ae1da1MUSRb1P/nZj0asuzPs7sT6HnV0HXQQ1nGEAQaQ1WE11kCnMVx11fCJT8TQcCRUdHUV0VghfELzfnYDAgLig6ruhrtxErIterqwGrParurbERVVlZl1O/PcczJvZkPlHEqiT0l5DX1VWE5/LiijP+adpS9zT9MX2aV0ouIxjY/8SsG+IsPxMwV7t1Ggp5AC/gLSu38ivSub9M5M0jt+IL1tI2mtGaS1bCCteT1pjd+S1vBX0ryracy7kibGOmxDtqamxrLta9eu0fz582lwcHDaMTY2Rhs3bhQHrl+/fk0rVqygHTt20Lt372jevHl07Ngxev/+PZWVldHcuXOpvLw8nHfo0CGamJig5uZmUfb69evk9XpFudbWVgoEApSdnS3uUVm/3y+uX716RePj45brH0tbLRuNKDgn4t7Vt5erGumrLWX0p/xzlDIlgN9nnaDqF10UGtj/QQC9EMBWCvi3UMCXR3p3LuldP5LesZn09k2kt/2NtNZ00pq/I61pHWmNa0lrWEOa9xsaq19BWmM60YR1R8cKeizEgAhA4Mjj4sWL1NvbS2vXrg3nbdq0iUZGRkR1Hj16RCkpKSIPZXBcvXpV5NXW1tKCBQvCz+3Zs0cQG+SGDXwXRLR9+3ZxjYd0XaeFCxeK+xcvXlhucixttWw0omBSiWDw9Xv6S2EZpfx0hr7IKaU//HiCFmQeo9pGPwX7PZMiEL3/30XvLwTQlUN6Z9aUAL6f3vs3pRp6/1VCAGN1yynYdyYCZrW3qonx9u1bwmggP5qm0blz50RaMBgkHCD9y5cvZRFxhmBCodC0NNwMDw8L0v8mYxYJqtsarQpJJQIAcPG/9fRlzqmwAH73wxGqbfRNikCGP7786eGP6P0R/qST1pJGWpMx/PlGhD9jdV8TBKA1ZhCNv4uGtbI0u4mBMAeh0fr16+n48eOUmppKGRkZyuofiyG724q6JJ0I0Oj/eX1Ucr6athyqpIKD16nFN0jBl8WG8Ae9/1T405lPwb5/U7D3AAV69lHAv5cC/j0U8JWQ3v0L6V27Se8sJr1zlxgBJoLDsfh4VmXjQQyMDJWVlWJEaGtrm1U9VTwUj7YmpQiiOSfY909D74/wB5Pf70nv2kah4UoKDd2g0KtrU0cFhQYraELzRzNle1o8iGF7Iyx+QTzayiKYckag5x+Tqz+mk1+s/nyY/CL8CQ1XWXSl2mLxIIbaGs/eWjzayiKY8o/u22Y++cXqj3c16e0/0/jIfdJackX8Pz7CIpg9va09ySKwhpOSUgF/sWHtPy1i7X9y8js++kh8V2joP5MiePtcyXfHaiQexIi1TnaVj0dbeSSY8l5o5BZprfjhC2v/qeG1/zGvXPr8mvS2IhofrSGttWBqFejDsqJdJIhmNx7EiPa9nyMtHm1lEYQ9O0GhgQuTAgj/8gsBTC59YvlTHlpzpq2/CIerZHIRD2KYfHXck+PRVhZBhFsnAv0UGvyVAr5/Tcb+9atJ835LWksOBfyHSYREE7/9gSjCjK238SCGrQ2IwXg82soiiMEhiVIUf6MDciTDgbba/WER2I0w2094BFgECe8irqDdCLAI7EaY7Sc8AiyChHcRV9BuBFgEdiPM9hMeARZBwruIK2g3AiwCuxFm+wmPAIsg4V3EFbQbARaB3Qiz/YRHYA7+2ZoPxiCZOTAH/1TNB2PgJA5AsCo/LALuBBzXCbIImLSOI63qUYZFwCJgEbgtHOrs7KSWlhZXHF1dXUlPUNW9fjR7rhsJIAC8vs8NB9oSzWmcpnbhgUWQwIJJNBEsWrTIlaJkEbAILBObRWBtIfWzL5FyOKQ2VDCGXk4RAV7bjpf7Gus+0zWPBDwSWCbLTCK4fPkyfc6JfF9fHx04cIDS09MJ9Vy6dCnl5eXRpUuXPtq+pBBBU1MT5ebmRj1OnTr1yZPo+vp6am9vj9kONp/YuXOn6XNOmhOAeFVVVR8l3Ew98mzzIABsELJkyRLKz8+n06dPTxPEwYMHZ6xXUogAGzpg5xQcd+7coeLi4vA9dlH51JWk8+fP0+3bt2O2A5KzCD4tfMNeCBDAsmXL6MGDB9PIjjdhY8MPCBSdnZnIkkIERpIDKI/HEybsyZMnCTudHD58mG7cuEFnzpyhhoaGcH5FRYV4FQls3Lx5kwoLC8WIcvToUXrz5g0hDEAaDuzWgnIgN8BHGsBHOaRDcGfPnhXp2MYIz7pJBLt27RJ7D4B0ICZwBfGAB3amwZZMK1euFKHK1q1bxSYdkphXrlwJP4v9C+7evWtKWvkMzo8fPxYkx/PGdHkNISAsWrduXdR8lEt6EUAQICtA7O7upn379tGzZ8/CIsDQCocgDyEVQiufzydIjvftYxIG5164cEGUGRoaEuUePnwo7o8cOUIQGkSAEQPfh9Cprq5OfK+bRIA4HB0JyA+S435gYIAQLkIYILd8txGuIRSQEGEU8rGbDZ7dvXu3uAfOksxm59LSUlEWRDYrg73RYB/7IkQrwyLweMI9OIhqJoKenh5Bbowk6NmxcR3SJLllOATBQBQy/MIv2BAP7qWI8AyO+/fvu2okACElydBeiACjnRTBkydPwvlPnz4VxMS8qKCgQIhGPosQBx0GsJNpZmd0IqtWrRJbQZmVgc8gArM5C4vA4xFDqiSmmQiQD4IXFRUJMu/fv184CenGOQFWI6JNwhESIB0jhfwut88JNmzYQJiUShEYly1xDWJWV1eLEQI7W5qReKZ04A072PXSrBy2iFq8ePG08MtYlkXg8QgnSWJCBOil5D3mCiC/7P2x6RxCIIQ5mGBHiuDWrVtiXoHJOA7EpAh/cIYIsIwobSOedVM4hJBQkgtbtS5fvlx0EFIEmGvJfJAW5EVeVlYWYT4h8zASYE6F0FOmmZ2lHfgtWhlsDpiWlkabN2+Omo9nWAQRIkCvgZ4dZAd5QVyIAM7CKIAeDCQ2hjJwPrYjBell+IM5BGygHByKawgH9gE8RIV0N4kAk97nz58LjNBWkLyjoyM8EmRmZooOBMQErmvWrBGLBcAPZTGPwhwC4SQENDo6akpcI+Exz8LzeA4jrcyD/7AFLMIyjDgyPfLMIvB4xKbRsneWYQuchJ5ejgQgOFaOkC5Xg+BwPIfJHtKxYoR7iAb3OEB0OAPpmEfs3bs3nAfbbhEBiIb2gow4cI/VNBBOjgSYLMt8rBbh5bjIx6pZSUlJOA87W967d8+UtJEkxqgjhYDfCnJyckTvL+uBjijyGeN90olAkn2mMwiPiV20MgAcowF69mj5Mg35cmlUpskz0vEd8t7sDEEanfW5r0Gqj9UBvTcWDBD+ybJSBMANO9D39/eH82QZnCEGjBLGNKvX8AsmwBgN0PujoyorK5txriBtswj4zyYsk86KCCSxjGejCIzpiXLNImAR2C4CLCRg4ooJb6IQ31gPFgGLwDIxZzsSGAmXiNcsAhYBi8Bt/2OMyaTZRNNp6Yk2MU7EXlxFnVw3ErjpH+2t/NmAChIkuw3XiSDZHcrtj/1Ps1kE/N4hy3MCtwqMRcAiYBG4bWLs1t6K2xV7mGMFM/x6j1+pVX741ez8anpHvZofAsCfa6j88CYdKtFkW45EgEXgSLdxpVUiwCJQiSbbciQCLAJHuo0rrRIBFoFKNNmWIxFgETjSbVxplQiwCFSiybYciQCLwJFu40qrRIBFoBJNtuVIBFgEjnQbV1olAiwClWiyLUciwCJwpNu40ioRYBGoRJNtORIBFoEj3caVVokAi0AlmmzLkQiwCBzpNq60SgRYBCrRZFuORIBF4Ei3caVVIsAiUIkm23IkAiwCR7qNK60SARaBSjTZliMRYBE40m1caZUIsAhUosm2HIkAi8CRbuNKq0SARaASTbblSARYBI50G1daJQIsApVosi1HIsAicKTbuNIqEfg/lPQ8cHvY8JYAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Check whether the environment is setup\n",
    "Verify the kernel on the top right states bpc -\n",
    "\n",
    "![image.png](attachment:image.png)\n",
    "\n",
    "You may have to refresh the browser window.\n",
    "\n",
    "Then, go to kernel, change kernel and select bpc - \n",
    "\n",
    "![image-2.png](attachment:image-2.png)\n",
    "\n",
    "\n",
    "\n",
    "## 1.2 Importing python packages\n",
    "\n",
    "We use several Python packages, and the following allow us to work with the BIDS dataset in this notebook:\n",
    " - [PyBIDS](https://bids-standard.github.io/pybids/#) ([Yarkoni et al., 2019](https://joss.theoj.org/papers/10.21105/joss.01294)), a Python library to centralize interactions with BIDS datasets. For more information on BIDS, see: [https://bids.neuroimaging.io](https://bids.neuroimaging.io).\n",
    " - [Pymef](https://pymef.readthedocs.io/en/latest/index.html#), a wrapper library for Multiscale Electrophysiology Format (MEF). For more information on MEF, see: [https://readthedocs.org/projects/pymef/downloads/pdf/latest/](https://readthedocs.org/projects/pymef/downloads/pdf/latest/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# packages importing\n",
    "import os\n",
    "\n",
    "# for scientific computing and data visualization \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from sklearn.decomposition import NMF\n",
    "import itertools\n",
    "from tqdm import tqdm\n",
    "\n",
    "# for handling neuroimaging data\n",
    "import bids\n",
    "bids.config.set_option('extension_initial_dot', True) # Disable extension_initial_dot branching\n",
    "from pymef.mef_session import MefSession\n",
    "from nilearn import plotting\n",
    "\n",
    "# for this tutorial\n",
    "import functions.helperFunctions as helper\n",
    "import functions.pyBPCs as pyBPCs\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Loading and checking the BIDS data and metadata\n",
    "\n",
    "An example CCEP dataset is available on OpenNeuro ([link](https://openneuro.org/datasets/ds003708)). This dataset is formatted according to the Brain Imaging Data Structure ([BIDS](https://bids.neuroimaging.io/)) and contains one subject to work with in this tutorial. \n",
    "\n",
    "This dataset includes an electrocorticography (ECoG) dataset with single pulse stimulation, and accompanying metadata, such as electrode positions, channel information, stimulation events etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change this path to the full hardcoded BIDS data path to be analyzed\n",
    "BIDS_dataset_path = '' # '/full/path/to/Basis_profile_curve/data'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Load BIDS metadata\n",
    "\n",
    "We use pyBIDS to initialize a BIDSLayout: this will index the files and metadata under the specified root folder (should be 1 subject, 1 session, 1 run)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the layout\n",
    "layout = bids.BIDSLayout(BIDS_dataset_path)\n",
    "\n",
    "# Print some basic information about the layout\n",
    "print(layout)\n",
    "print('Subjects in this BIDS layout:', layout.get_subjects())\n",
    "\n",
    "all_files = layout.get()\n",
    "df_layout = layout.to_df()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Select run to analyze\n",
    "\n",
    "User input regarding the subject, session and task to be analyzed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change this to make sure it fits your data\n",
    "bids_sub = '01' # The subject label\n",
    "bids_ses = 'ieeg01' # The session label\n",
    "bids_task = 'ccep' # The task name\n",
    "bids_run = '01' # The run name\n",
    "bids_space = 'MNI152NLin6Sym' # The electrodes space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can load filenames and metadata files needed for this analysis from the BIDSLayout:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve file names\n",
    "\n",
    "# pyBIDS does not yet support MEF3 data, we use the OS module.\n",
    "for (root, dirnames, _) in os.walk(BIDS_dataset_path):\n",
    "    for directory in dirnames:\n",
    "        if directory.endswith(\".mefd\"):\n",
    "             mefName = os.path.join(root, directory)\n",
    "print('We will load the following iEEG data:',mefName)\n",
    "                \n",
    "dataJson = layout.get(subject=bids_sub,session=bids_ses,run=bids_run, task=bids_task,suffix='ieeg',extension=\"json\")[0]\n",
    "# print('\\nSidecar JSON (*_ieeg.json):',dataJson)\n",
    "\n",
    "channels_tsv_name = layout.get(subject=bids_sub,session=bids_ses,run=bids_run, task=bids_task,suffix='channels',extension = \"tsv\")[0]\n",
    "# print('\\nChannels description (*_channels.tsv):',channels_tsv_name)\n",
    "\n",
    "events_tsv_name = layout.get(subject=bids_sub,session=bids_ses,run=bids_run, task=bids_task,suffix='events',extension = \"tsv\")[0]\n",
    "# print('\\nevents_tsv_name:',events_tsv_name)\n",
    "\n",
    "electrodes_tsv_name = layout.get(subject=bids_sub,session=bids_ses,suffix='electrodes',space = bids_space,extension = \"tsv\")[0]\n",
    "# print('\\nElectrode description (*_electrodes.tsv):',electrodes_tsv_name)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the pandas.DataFrame.head function to quickly test if our files have the expected data in them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_channels = channels_tsv_name.get_df() #  Get file contents as a pandas DataFrame (only works for TSV files)\n",
    "# np.shape(df_channels)\n",
    "df_channels.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_events = events_tsv_name.get_df() #  Get file contents as a pandas DataFrame (only works for TSV files)\n",
    "df_events.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_electrodes = electrodes_tsv_name.get_df() #  Get file contents as a pandas DataFrame (only works for TSV files)\n",
    "df_electrodes.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Render brain surface\n",
    "\n",
    "Render electrode labels in MNI space. This should open a separate browser window."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "marker_colors=[]\n",
    "\n",
    "xyz_list = df_electrodes[['x', 'y','z']].values.tolist()\n",
    "electrodeName_list = df_electrodes[['name']].name.values.tolist()\n",
    "\n",
    "# gray out bad channels\n",
    "for cords_name in electrodeName_list:\n",
    "    if df_channels[df_channels.name == cords_name].status.values[0]=='good':\n",
    "        marker_colors.append('red')\n",
    "    else:\n",
    "        marker_colors.append('gray')\n",
    "\n",
    "view = plotting.view_markers(marker_coords=xyz_list, marker_labels=electrodeName_list, marker_size=7.5,\n",
    "                             marker_color=marker_colors, title='Electrodes in MNI space rendered on a standard brain') # Insert a 3d plot of markers in a brain\n",
    "view.open_in_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Preprocessing data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Look at the rendering with electrode labels and define the electrode of interest (e.g LMS2) and time window to be analyzed (e.g. -1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define electrode of interest\n",
    "el_interest = 'LPS2'\n",
    "# epoch_size for visualization\n",
    "epoch_limits = [-1,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Extracting relevant information from metadata\n",
    "Because we have the data in BIDS, we can automatically pull information essential for analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get electrode of interest from channels.tsv \n",
    "el_interest_nr = df_channels.loc[df_channels['name'] == el_interest]\n",
    "\n",
    "# Sampling rate from sidecar JSON (*_ieeg.json).\n",
    "dict_dataJson = dataJson.get_dict()\n",
    "srate = dict_dataJson[\"SamplingFrequency\"]\n",
    "\n",
    "# Construct a time vector\n",
    "tt = pyBPCs.create_time_vector(epoch_limits,srate)\n",
    "\n",
    "# Number of trials with electrical stimulation, good and not stimulating electrode of interest\n",
    "df_include_trials = df_events.loc[(df_events[\"trial_type\"]=='electrical_stimulation') &\n",
    "                               (df_events[\"status\"]=='good') &\n",
    "                               (~df_events[\"electrical_stimulation_site\"].str.contains(el_interest))]\n",
    "df_include_trials.head()\n",
    "\n",
    "# Filter channels to select good channels, only ECOG and SEEG\n",
    "df_good_channels = df_channels[(df_channels[\"status\"]=='good') &\n",
    "                               ((df_channels[\"type\"]=='ECOG') | (df_channels[\"type\"]=='SEEG'))] \n",
    "# df_good_channels.head()\n",
    "\n",
    "# All Pymef operations are handled through pymef.mef_session.MefSession.\n",
    "session_path = mefName # path (absolute or relative) to the MEF3 session folder\n",
    "password = '' # Pass empty string/variable if not encrypted \n",
    "# Either for reading\n",
    "ms = MefSession(session_path, password)\n",
    "\n",
    "# To obtain information about time series channels we can use the Pymef function read_ts_channel_basic_info()\n",
    "# ms.read_ts_channel_basic_info()\n",
    "\n",
    "print('Sampling frequency is',srate,'Hz')"
   ]
  },
  {
   "attachments": {
    "image-4.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcQAAAGqCAYAAACVoTQ1AAAgAElEQVR4AeydB3gU17n3U2+Sm89xkhvnxnGq4/TmXPeCjbuNG3aMCy70KnovogoBogkQAkQRAgSii15EE0iiSyAQRRJFoogqegf7/Z7fIbPsSitpy+zu7O55n2fRMjtz5pz3zJn/efvXRJPmgOaA5oDmgOaA5oB8TfNAc0BzQHNAc0BzQHNANCDqh0BzQHNAc0BzQHMADmgJUT8HmgOaA5oDmgOaAxoQ9TOgOaA5oDmgOaA5cIcDWkLUT4LmgOaA5oDmgOaAlhD1M6A5oDmgOaA5oDlwhwNaQtRPguaA5oDmgOaA5oCWEPUzYHDg1q1bUlRUJIsWLZaJExOdfqZNmy7p6eskPz9fLl26bFyq/2oOaA54wIHz58/L5s1bZOrUZKfrLTFxksyfv0C2bcuW48ePy+3btz24i77EHQ5oCdEdbgXwXBbDmTNnpKCgUDZv3iwrVqTJ7NlzhEUTFxcvI0fGyahR8ZKYmCRz56Yq4Nq5c5ccPnxYXXf27FnhU1paKocOFcmWLVtl4cJFMmrUaOnYsZN8+unn8vHHtd36RES0kMGDh8jcufMkNzdXrly5EkAO6VtrDpjPAZ5p1hDPN5vBBQsWChvDsWMT1Jpj3fF9+vQUWbZsuQI4NozHj59Q6401d+LECdm9e4+sXr1GgV+/fv2lUaMmbq011uYXX9SVyMgeMmFCoqxdmy5FRcXy5Zdfmj/oMG5RA6JFJp8HG7Dat2+frFu3XoEdYNWrVx9p3ryFfPLJp5UuoM8/ryMtWrSUtm3bSevWbaV+/YaVnm8Pfm3btpcxYxJk8eIlkpGRqXak7EqNT1bWBlmzZq0sXbpM5syZK8nJ0yQ6ur/Uq9fA4R70sWPHzjJlylTZuXOn3Lx50yLc1d3QHHDOgRs3bsqRI0fVs75kyVKZNClJYmIGS4cOnco93/Zrhu887w0bNlbrjXXHOq1d+zOHNVH2Gv5fr159GTBgoKSkzFTratOmzba1xppDaly/PkPS0tLUpnXGjFkyfvwEdZ+y7TVo0EiGDBkqK1euVBtf56PUR13lgM8BMSoqSmrVqiUZGRmu9imkz7t69ZpSOSLhoZoEWAAwAM3+YWdhscC6dYuUQYMGy7hx42XGjJlqFwpA7dq1S6k4AdGKgOfKlasKXFm4SHJ8ADN2s8uXr1C73gsXLnrF73PnzisQBzBR/QwYECPNmjVXY6lbl4Ufo9SwejfrFZsrvDgnJ0datGDD9IlaZ40aNZIxY8bIxYvezWuFNwzSH06ePCmbNm2SWbNmy/DhI6Rr125OpTQ2ku3adZCoqH4yYkScTJ48RVJT5yvgAqwKCgqkpOS4XLp0yal0xsb2zJlSta4BK6RJ1ty8eanCui0uPuyV6vPGjRsKwNHwYN4YPXqMdOnSTT777Au15tq376C0Rhs2bNTPgAfPqk8BkYfmO9/5jnzve9+Tjz76yIPuhcYl6P9RtfTs2dtB0uMhBvB4qFkwmZmZsnfvXrXTM0MVAv8B1e7dIwUbIR9/0enTpxUYo8Zt3LipWqxNmjSTiRMnqTGaMT5/jcXK90lISJCvfe1r0q5dO+nTp48AiD/4wQ/koYceUqo6K/fdl31jk7h161a1tlgD9pvNiIiWSgpkA4eaMzs7W6lF2ayaQZmZWUpSzMvbbUZzLrUBUO7Zs0cBfvfuPdT92VQD7CtXrtLg6BIXfZypZtasWfL1r39d2rRpo0ARG1g4EA/n9u3b1U6tffuOtsXYqVMXpZLBFoFdwpdG8q+++kqGDYuVFi1aKVtGIPkO+OGww8snJmaQsoU0bdpcEhLGKVvm9evXA9m9oL63AYgXLlywjQMbMyDZsWNH27Fw+FJSUiJpaSuVRgV7GyDYoEFDGTgwRm04san7WnLOzy9Qzzf2/UDS1atXFdAnJU2WNm3aKoDs06evkip59/B+0FSeAz6VEGvXri1PPPGEbNu2TS3Q0aNHl+9BiBxB+tq+fYdyUqlTp55ajOzQ+vbtp2wBGNf9SewKcZQ5ePCgP2/r0r0AQGwkffpEKT7x8ho+fKTk5GzXC9UlDt49yRkgsiEDEBs0aHD3xBD9durUaaV9wXZtSIHNmkWojScSkz81EZcvXxY2epGRPX262fVkKg8ePCSAI2pc+IRaGA/W0lL/vpc86bs/r/EZIAIQP/nJT2TQoEHqobz//vulWrVq/hybz++FWgZJEK8vQy1oqClWr16t7Aw+74STG2BXBJSxf1idjhw5IhMmTLQ5MLRr116ys3Os3m3L9M8ZIHLsG9/4hvAMhiLxEkfb0KvXXRMEIMjzjjenP0HQnr+oYFl3mEisSryzsGUC2gAjZhs80wkB0eRDlenixYvlm9/8pvDCg1DfoD5FXA9WYue9b1++8vxC9Wd4WeJCHRc3Skk99qqrQI0zPn60tGzZKqjCIJAacczB3slCjY4eoJx1AsXDYLmvAYg1a9aUDz74QB577DG57777ZNmyZcEyhCr7eerUKcnKylIbTyQbng82nrzU8XouLNwfMBA0Os97DnBZv369ccjyf9EesRlFQ1O3bj3Fy3A3X/hMQvzwww/l6aeftj0UqC9Q48TFxdmOWf3LtWvXZMeOHcpTrHfvPjZPLkPlgAcZi9FK+vji4mLluIP6MVhp16486dq1u3rxDR06zLapCtbx+LLfBiDizR0TEyNNmzaVe++9V/7yl79IYWGhL2/tk7ZZS0ePHlUxe8TXIvmx3vggffE8EM9nNYkG0wghUsFIqHpx6sPDFk0X/PWnA56VeOYTQDx37pzyLn3ggQfkhRdesH2QGF966SUrjd+hLzi54BlGWILhqcVCvBMQ21MmTZqs4vTwoLQq4UhD8G6wEy9G5gKnIOYAgCQBwLFjJcE+NFP7bwCivWYCBybMFX//+98DLjm5MlhU/Ni8CYfAE9kAQOJqhw0brjQyaGYqCi9y5R6+PGf//v2qz1u3bvPlbXzeNg5HSIzwH3BE07RlyxYJJ6nRJ4CYkpKipMGePXsKC9b4EI8IKFoNULBJECtkbwfE4QOjsxXUMa6uBB5oVEnEGIYKIQmwQTHmhpjKoUNjVZaOUBmjN+NwBoi0V6dOHbUGkbasSGx4duzIVfF6BgAiDRJvu3HjpoB7RrvDM2yHjRs3sZwjjTtjsD8XG37fvlG2ELGGDRsp79RwkBp9AoioS5EOyzKQoFbsiLGxsfb8D9h3dqbsiAiKR/9PxgdsFagQgpEIHOblYthtg3EMFfWZZwm1Oy8fI64Mu+3p0+ERylMRXyoCxOeff15tPu0lx4ra8OdxgBCnjg4d7oQjkZSCEIVDhw5ZyvTgDk+wZY4ePdadS4LiXMLkVq1arRJ68H5EW4N3OHMYqmQ6ICJe33PPPdKhQwenPHvkkUfk4YcfVr+xWH/1q1+pTBtOT/bhwQ0bNqgYJQBk7NhxQsaVYCdSTzEeq6qWzOIvXoRk4sDFHXU2izRcqSwgspnr3bu3kg7r16+v2DJjxgwVrL9x48aAsonQI1KW8Yy2bNlaCGAPlEeomYzAqY7EG6FMaPVQaTN3ZJ+ymg3XLN6bDojz589Xi5HgYGcUGRmpfj9w4IDwwdHm8ccfd3aqT46xKGNjh6uJbd48QnDgsCrhBbZp40bZ62I8Fcm6eWB9GfDvCq82ZGbK+Ph6kjTmSUke/5hMiK8hU5NGmy55kyKLVFWMGVUbXsDhRgYg/vSnP1VamW9961vCp3HjxjYv486dO6t1NmHChICwB9Bbt26dTe1NXs5QmitSFLL2woHId8x6Q1oMRqetqubIdEAkZ2BeXl6FYjU7C343MkawUPwRinHgwEGVxgn1KA8wdilSm1mNUEcsXzJTZiS9LlvWVpdDu16THVkvyuwp1SUleVClfU5PT1cPa6DsRkimk8b3kNwNL8mXF95Tn5tna8rt8zXlYsnbMm3Ci7IzN9dUllONgBAYFmmnTp3VJsvUG1i8MWLelixZImSFmjNnjqxdu1YlibfvNlIj55Q1Ydif44vveGlTlYFE2cxPjx69VJygL+4VyDZR+1LxIpAEr0mAgnMSsdG+dIRBsifkDH8Fgv3JihMqZDogWokxgAvhB6QsYkGSuYW8odgOrUpj46MUCBqAwt/dW15R4HLt9LsSH/txhSpRbIeMM1C71TGjekhp8Zs2MKTvQ6L/Lm0i/mA7ljzhxUpB3ZN5QQKZOXOWcgJgkWKTCgVVnCe8sMI1eJmz4WTjyfNIJQgqOljd9gSokDCbZPr4E2CCcAVYCA/BKzNQUi+Of4Byj569JWbQUOke2VPZ2amc4ysiWb/hAU75OEpchQKFJCCePHlaBZm2atVGLUicMHhJ+jt9mrsPyILUGVKQcwf8DEBMX1pN/vn3HwpgyLFzR96SieP6Vtg04I+3nr+l34L8Alm39Dkb8NHXw3tel5/f/30p3PGq7fjlE+9I0sTBFfbfmx/IVYltipcwNeesvPHxZpxWvBbtABUYABI2JWw+qa6C6cTfkqkn/CH4n6oR/QcMkvETJ8m4CZMkun+MSr5f1XNErUSeObzS/U1o5ACmiZOmSMrMOf/5zJZJk6cKZd2q6rs3/UUyJBabuSZGlMQawU4hA4gsOrzXjPyYPKBUl8CFOxikBfo4fqQjoNworSl//P0PZO2SajZAAWhSk5+qsPbZsWPHlNcsfPDniyh1zgSHPtLPpg1+KxNG/V+544mjX/bZumFHbx9LRbFWTb7jAM8b9S+RkFhz1AckhMlqoVVVcYDY3dFjxtmByh1wGTU6Qb1Tqroe8McTc+9e30llzvqAJD58xKhy/QYch8WOVF7Zzq4z8xje30ZYFAn7/fneMXMctBX0gMiCJIchsTIsSHTbTApu3MFEOBhtWPm8A3hMHf+IfPbxb+XEgbfk6qk7EiJAc6boTVm7Zk2Fw0PtAy+o5m22g83tmzfk0rEiuVDkCDTz5/R16PvJA2/KL3/xE3n5hZ9J24i/OPw2I/GxCtW+FQ7KzR+wpRCzyIvaignO3RyOpU5HNUgRa3KJ8pzxIUsLEoIrKkZLDUZEqJLRslUbp6ACsEREtBIkyMro4sVLKl0iauLc3J2VnerSb1+6WKpt0KAhkjx9ptO+T5qcLL17V6xNcqkjLp6E9q1btzvZpXBwC1YKSkBkQfLCM9J7GeoZ3MoDpcf39gHI3pYt++1UiwDfsAH/lLdrPCD/+NtdlSnH+SyYP7fCW2KrSU1NVS8qHE7MCsM4uTVbdgzpIXuTYssB4oK5MQ6ghzMN6t2HHvyhHMi9qzKl71MTnvCL1E6iZ7L7I7XoDDcVPi4u/8AmMzFxkuInIIhJAlNEsPOWBAE9evZxCioAYrfuPVUMbFWMAlhbtWqtVIjUYvSG8pOmyYoaNSWjcQvZHj1Y9qfMlrN7ykufeMxPnjrdad8nTpos/fsP9KYbbl2LCnXIkGHqvYMqNRgp6ABx8+YttvRO7MaQDkOhzqIzCRHweKn6z2TdMkeVaVUSovEgwhteXGSz94bO7suXvPhxktWqueQO7y35U+OkNG+rfHX7tq3ZZUsWyvmjb9lAkT52bf8nGRdXXmU6Pv5T23W+/oJHHDzAtqrJMw7gET5ixEjFR3iJJEASiGAwRbgy4kOHiqRdu45OQQVAbN2mncqv6kpbSJJoq9BMGJ70rlxnf87Ny1ekePFyWVr9VZn/8KPqs/bjOnLZSdrC5cuXS1S/AU77znGqgviT0Ejh7c1zgiNVsFHQACJOIiNGxClGEwiLARuX+1AhXi7jRjo61PTr+Q+p9d5vBWnLkAz5O2vyM4InX1XEji0iooVycsD4zscVun3jhpRkbpCcvgMk/fMGsjUySoqXrJBT23Jk5/DBcnB+spzclik3L97tA679E+NfsPWztPht+dEPvyM7NzlKh3u2vCzr0v1nfEda7tcvWj035GXU5B4HeKkhYfOCw/sSdaDVvUXdG6EoYEfbNHbcxHLAMnrseFVVw50xs2mHX4QkuEqXjx2T/EnJktmkpZIK946bJFnN28qCR5+SPWMTHTaf9m3iGUve5cFDYmX6jNmq//wdMmyEUpeapR2yv2dV34ntZvw49wWbCj0oAJGMNkYBUOJ9Qinuxf7hWr5siUMMX6O6v5VZk5+wgQxgiOQ1OXG4/WWVfjey11BJnE9ldOnIUdnef4gsf+0d2dazn5zJ3VXZ6eV+27J5syxIecbW3yN737B9p+/nj74tExO6lLvO1wfYtRN/GgpJz33NK/v2CcLGDsuLbedO7+1i9m1b7TvxnKiAcUSZljJLfYYOG6HCGaqyHzobC3U92bhXRrdv3JSihUtlXd0msvbjurJ/5ly5cf6CugQAzopoJ+cLD1TWhPoNUCTbVqPGTVV/cXAZM2as6YkwquyI3QmGRoGMYMFElgdEVIm4FZMFP9izybvyYMxKGSvZ619UQIKXqb1keLzwDZmY0NwtLy6C9Nmt4Y3GpywhDRYvSZPMpq1lY6uOciRtjdy+ebPsaS7/f+fOXJk64X3ZveVluXXuTv9PH3pT5iY/KynJQwNm4126dJniA17HmirnAC9jPEV5bkaNig8pTUxlI8f0whohvdyAgTEqttXTUC027vDPmYR04VCR7IiJlfTPGsju+HFyYf/Bct366ssvxVXHGuNi/CcA70BIhUYfjL9o9MhX26NHT7feV8b1gfpraUAkGTi5KiMiWlq6CrW3k3f+4D65dPTuosjMWCdjhr8paxc/K/k5r0hW2vOSMOJFmZmS5LbdhvRmLEycH/gYxO6UHWna2x/IjoHD5MoJ19SpxvWV/eWFunXLVpmSNF4mJ8bL4kWLAv5SpU/kYiRGToNixbMHnyZNSlLPTGrq/IpP1L9UygHs90jX9iEIpXl7ZGPrjrL2k7pybM36kFM9l2UI4NygQUNV7BuTSjCQZQGR3RpSIWqMUHCaqehhKJg+UfKnxMnVk8fKnYJXH7lMyTjhqQMDNQUBRJJh8+GFh5om7c33lY3wShXu5OU6FcQH2K0THsCLCvtOsHok+3IKDDCcMWOmL29jibYvHT4iJzf5xq5MXCIp66DLx4/LpnZdlVr02NrQB0L7ycXblvXWpk07IX2m1cmSgEhJn5YtWykmeqqysDrj6R+qk21RnWTv5OFyZPV8uXBon+m7RnT5xGbuS12kPumf1pedQ0fKtdKzwcAi0/sICMbG3snaj106GBap6UyooEEKMLN5osQWG6dQoi9v3Vbr7cjKtcpjek2tz2VT+25yyweOeTi8kbllTtIUyekzQNbXbyolGcFlSzNz7kmfSfgT2plZs2abHhttZl8tB4i47eL8gKo0FOv6GZPH7nT1B5/Ipi5tZPuQHnI0fancvn7N+NmUvyT5bfL+v2VitRdkS9fe6nP97F3PUFNuEoSNIG0TJ8XLn0W6enXFSQ6CcHgedZkYXvgxcWKiR9db/aIjq9Jl2Utv2sIYAENfgD7P1qCYQdLztTdk2RvvytHV6VZnjV/6V1xcrMq18YyR5N2XKeW8GZClAJEHND5+jApsZVcRqnR6e64yqJfu3C0Yzy/5oKr5vtxcGfNWTVnTsJmc2eGet2io8r3suKi0wsaLRUrAebiqUHFcwwuX/K9mZzYqy/NA/f/W1WtKIlz4+DOy6v1P5NYV8ys0EAY2vXdfSaxWXfYkJMrt69cDNVxL3pdKR1FR/dR6wwOX3MNWI0sBIi8oXk4pKTOsxifT+nNi42ZZ+U4tuXLcd9nhN0yeKlOfeFqSe/QK2RecWRNSVFSkwgp47rp1i/Q4mNqs/vi7HZw+8AZEpRWq5gk22jiz5E+ermL7LpccN53Np48elTkf1pb4F1+RPB3vWiF/ed4SEsar9zwZxqzm4GYZQCS2kAroZGgP1Z36xeLDsuKNmsJfX9DN69dlbbceMuaZ5yRhWGzI8tFs3hGDZlTJICG8M1d5s+9plfaw6bAZCOWkBbvHTFDpz+D52X0FprP+yKbNMvO5FyX6k09FJ5N3jb2EQWGuwOGmomLyrrVk7lmWAET07nhlkYotVO2GxBStr99MSnftMXcG/9Pa+YL9suzjz6XfO+/LunRtt3CXycRNDR16Jw9jZGRPuXDhortNBN35e/bsVS8lkuGHKh1evkrW12sqhBmZTbS5K26MjH3+JYkfMMCydjGzx21We3ig4nxElRCrSIqWAEQK2rJLXbhwoVm8tlw7exImye6xE33Sr4LkmZL89HPSrUlTlbnfJzcJg0ZRrZGpn2exb9+okFY3YyskmwqhTaGUAtH+Mb125qyseP1d8YUj2fXzFyS9TiPp88rrMiOETTz2/PTFd/LIkhoQSTEnJ8cXt3CrzYADIoZWHBswspKCKBSJWL+VNT8SDPtm0q1r12Rzl14SU/1lad+ydchUrTaTR+62hbbCSIo+bNhwh8Bqd9uy8vlGSr9Vq1ZbuZte9W1HzDApXlp5ukJPblCSkSWpr78tbd59T+CjL7xVPelXsF6zf/9+BYoE8Qe6VFvAAZEaauzISRUVqrRrxGgpXmRu1vlrp8/ImlqfSuTbNVUmn1B1iAjUM4EakecyVNWJSIeYKKyQ5ssXc3zr8hWVmN7stgunz5Ipz70kdWt9qGLqzG4/XNs7fPiwipdGMDp16nTA2BBQQLxy5ap07txVvXgOHz4SMCb48sbsHslZaF8qydv7XThwSAiwnx49UD1EoWp39ZZP3lwPUGDXBhRDLU4RqZBxkbQhVKkgeYYUTDXPW53A/h2DhsuaFu2laaMmEhMzOKRV6oF4Lkggjz2RDD+BUuMHDBDxamvevIVamKTTClU6k5sn2X36mza80t17lafqgY2bVewYKhtNvuEALuLYEjH8nzjhuzAZ3/S+fKuYJ0aPHqPWHIC4aVPoJjpf+1EduXKy8ir35Tnk/Aib2i2de0ru4BGSNGmyCtO5fDl0Ss85H3VgjmZlbVDP54QJgUkQ4XdARLVHXTUWJLsBsumTszNUiQDdooXLTBke4Rqr3vtYLhwsUnYuNhShGqJiCsNMaITSYzie9OzZy+N8siZ0w+smiPElhR/rjrhD0rSF6rODujTt7X97zTMaAAxz+g5UZdHOnbvj70DRaU2+4wBZpHCyKSws9N1NKmjZr4CYm5urYg1ZmHiWmuVEg8t8xvr1snD+PElbscJSnpbUNHNW3qWC+ajwMOEagOHZPfuElzT2n4yMzArP1z+Yx4Hs7BwFJCtWpJnXqJ9aQvU0fPidaveUNSos3B/yTiBkgtrSzRyt0/Z+g1TuXzJKjR8/QZl4tBONbx9ePKDRGlLYwd+1b/0GiOvWrVeojzHfrOoVPJizZiTJtAnPSknBnWK0FKFdPvcZSYhvF9ACmcYjQ3klag56Q+f3H1Du4xeLilUzVCKgRqRemN5w1b1rR48eqzzhAmXbcK+3d85G5dunT5QCczag4ULFi1fIvgmuV6uviC84w+VExaifUTcTSL5y5aqKTtfHTeQATjZIifYl60xsvsKm/AKIlC8iTQ/ZQKhkYQbh9DB6ZDM5svd1WxHd2+dryoIZT8nFkrflysl3ZOyIN+TkSfPq/HnSbzLTeEN4k+JAY0iZjJuK2MGW3g5VuVUT+royPxcvXpSGDRtJSkpwlEViszRq1GgFhsuWmevhbPCroKBQVcbALknGG9LgWYEKp8+Wg/O82wDkjUqQTe272or0wkNMPGijgoGQso4cORoMXa2wj4Q/NWrU2K+Zo/wCiN2791Bof/CgebbC+fOmyJmiN21gSGX5JXOelocevEcuHX9HHQcgpyR2rJDh/viBwGBP6da16yq0wr5mGwlxsQPl55ufgsrTflZ2HS9mygnR57y8vMpOtfxv2DYIIg6GBNibN29RPB82LNYnfJ0/f4G0bNVGRsWPlcSkqRI7YpQ0a9ZCVq0KvASFhyk1Pz0lQqQymrR08Azv27efDBo0xNMm/XpdUVGxKp1HsvZgJooKE6OOdtFf5HNAJNCSlyHOM2bSrKnvOYAhgFjt6ftkzeJqDscXz3o2oAmbAURKzxjqTnd4gJt32ew2qEvr1KkbFAHjBLkb8XwTJkwMehUvWTV4lgkktjr16NFT9bWgwPyNExqfiIhWMi1llqTMnGP78H9AMpBxZMzLgVmpsn/mXI+miKT7a2p9Jjcu3pUEb9y4qTy609LMD/L3qJOVXIQjCv4FfPbty6/kzOD4afjwEeKrTZ0zDvgcEA3pwMyFifpqxbxnHIBv5uTH5fVX/leyVj4vJ/bXsP12LP8NWR2AXWtJ5gbZNTxe1WBb9NTzLqePMlzFqay9sU2nciDSu3cf6d9/oLO5tNyxKVOmqpcySRdCwd7JGEg+v3jxEsvx2r5DpaVnFd/pqy/4PnZsgiSMT7QBoT0oxg6PUxoB+/74+ztrZ1fsKLdve+vqVVn7cR05lbPD4VrqirIRIgm8lenYsWPKI7pZs4iQ8dzPzMxU3tHYw/1BPgdEXuD16tU3dWEeOnhINq+pbgM9pMPYgf+Uzz76rTz0u3uU/ZBjfG6UviszUlL8wUuHexCIn9Goha0gKV5qrlBmszZyLr9QxRqWzcGIxIUKAZd5q1NGRoZ6ibDD88VLOVDjx66B16aVaceOHYr3SUneO5Y4GyeqOHsQtP8+cdKUgG/YLhwqkswmrZx1vdJjuUNGqhJRZU9iQ0d5LCsT3pj4aCAZAoyhQvicsBnBXu0P8jkgskvt2LGzqWNh8pfMdpQQAb/33n5Atqa/4ACUh3a9JpkZGabe39XGbl66LGs/qScLH3va1UtkafVXZfEz1WX7gCFSsj7LwY6BDZaHgzAAK1NJSYkKZu/SpVvIxbpt2rRZWrZ0/2Xrz/kiHIfnxFdhIuR4TZoyzSkoxo9JkDFjxvpzuOXudfnoMcloWld2xQ2Xa2eqduLDca14WZpkNW/jdPNGCABlwaxMEycmqjkPtVAsNtOUBUxNne8X9vscENm1UHjVbJox+VMH4Gvf8lUzSOgAACAASURBVM/y2UcPyu1zNR2Oz02u5lcvpbLjvHLipCx8/GlBHVMVAaDzH35UfajsfWyNozHZyPtq9awpQ4YM+48T1cGqhhx0v1OmBrA5dcqcLCi+YEB6+p1C274KIF+/fr10697DKSB27NhF9uzxTYkzV3l1YN5s2Tmqn+yKj5YTWxzXUNk2bl27KjsG95LtgyOlcMb0sj+rZAxoZeLjR5f7zSoHDhw4oNbbyJFxVumSqf1ABUz4kD/Ip4BIoVVclQcOHGT6WLIy0yUn40Ub+HVo9SeZN/1J2/+RGE8felPmzLzzkNy+cV3O7NridAdoeufKNIg9AzVOVVS6M0+BIUZ9Z+ejAiMWyspejrt371GAEWxhIVXNjfE7tnAAMd3CNSfnzJmr+rhjR67RbVP/8vzhJNc9spfyMJ0+Y7aMmzBJOnfprjIoGTcLlKr8wqEC2TspVo5npsmhRdOlMnPF8Q2rpWD6WNk9fpDsTSwPeqdPn1a8xJnNisRcoIlBirp06bIVu+h1n7p1667C9vxhR/QpIFLfkJeHr+KgZk4fLxnLn3MAQcN2iDPNsMHNFXhgz7t56YIUTBstJRn+9xRD8jtfeKDKB6NowVLZ3DFSKOvkjKKj+0uzZs2d/WSJY7wAO3XqrJIG+OPhDdSgeQGRRcOsBBNmjuPixUuK/6Sb83VqNuzE+Ai0bdtOBg6MEfIT29PZvfmyf8Zcv29CKcZ9dM0C+er2LUEC5P/O6PbNm3KxuFCunjym3g0Xi8t7D+fl7VbvsBUrVjhrIuDHNm7c+J/+BV8WJVeZt2HDnTH6I0jfZ4BIoC6ZBnr16u3ThZmVmSljR7wt65ZWk8LtrygbYmL8MzIrJUGpO7766ku1U9w/a4LsTRou2dFd5MJBa6rysGVURryEO3XqUtkpAf1t/fo7jjRLl5qTuzWgg6nk5rt333lJkpXf36mlKumWcqZAvfT553UkN3dnZaf67bet3fsK6Quvnzvvt3tyI4DQVeIdcXJbhlPgJgaOTT0vZasRG9DIyJ5KOvT15ifQYyf0AjzBhu9L8gkg3rx5S6n2oqKi/RIvx4Nx5PARycnOlsKCgnL3/PLWTTm4cJpsH9RDNnZoJiTJDjZC4kJdipRoRcIbDE+8Fi1aBtRm6y/e4OnLi9JKXrRGn8gZbBUinm/ZSzVk4RPPyrYeUVKaFxj7IhlmKtNaAIrOaN68VDXPSIpWI6PIc6g50jjjMxtP0n6y2fNl0m+fACKd52XhzwwDzphof4yFyWI8NH+JCrwNNlAkBgqexsW5H19lzwdffMeOERXVT+3gqGkWDsQmrF+/aDUnvvLmdJePhAeQ49ZqRF5RnMUWPPKkbOnaWy6YmLGqsrEWFxXJ2FEdJHH00zIv+UlJHv+4JMR9LhuyXK9WQUJv1p3Vao4SWsMGmXUXLmRk6SJywVeFrX0CiOzGeIjWrk235Fyd3b1X0mvXl32JUyu0L/ii44QjZGVlCZ6K7tqfdi5ZJm3feVcIdndGpHkjywZ2m5NbsoXg5OIlK+Tg3AWyf8YcUQVTKZpa9jNtpjqHElVHVq6V41mbBOceSkxdLz3rEPZR9r7ERW7dulW5pKPOCLfEx2fOlKpUbuxaSUYcaAIQIyJaBrob5e6PKWDBo0+pmoIZjVvI+YLytrpyF3l5YO2aFbJ0zvNC+kbDr+DqqXfV9+2ZL0pKsmvJLXAI5F1GhRkrEB7mM2fOUs6KmFBw+gknIuMV88FGxRfkE0BEYuAFaWVPQxxt8uLGysp3PlDxfr5grtEmqpoxYxIkokUrieo3QPpERau8j9hZkTQqIyplkHE/9V+PSssPPpDZAwdJftI0yek7QHi5UE2D9HAra34k6xs0U045nE+mfnbmgCGgSG5HZx8kZs6hzbz4cbIjJlaVztnQor1KKr6ixnuqfTxfN7TsIDsGDlOgmjF5qnRpEaEeTgCBwp7hSEiHLNDo6AEBH76hMrWSXdNgCs8V3p6UL1tT63PZ3n+wXK3CZm5c6+7fvXv2yIp5z9qAEEA8caCG/M+P75GDu15Txw/ufFXmzbmbsIO+XTt9Us7k5cjB1DnCBhPCTlz/s8/kzM48OZy2RgqmzlBZcLJ795dNbTtLRqMIWVu7nqz+oLakvfmeSqhBQv+yH0q3rfnwC1n3RSMh+cbmDt3VGmad5k+ertbm8Q2blPOdfdo4+kBmLlJgDhgQo6RCnrf+/QdYBqTdnR9vzkcyJJQPHmDLN5t8Aoh0ErE2Jsb8cAuzGUB4Q1bztpL+WQOhbMztGzfNvoVMn54iPXv1FdzTjawe5H1s36FThVIV/difMluWv/q2LTZx4lPPSlqr9gq8jqVnKNXT7f8sXDM6DXAjxaKawIOQkkFkZhk6NFY6tmgtzf/9gbR/+13p8cobMqj6SzKp+suy4KXXZH6N92RrZJQUTJspp7Ztl1tXnXvJmtFHq7UBz3hpskAD7chC3CH9yM+3Xg5Le09PvDvzp6Qo0Njef4gQSO8usaG9VnpKbl2+WO7SlKmdHMAQQOza/g/SqulDDscnxN2NUaa9k1vSZe/k4ZI7oq9sbN9KAd30Z6vJtOdeko1tOqvN4L5JyULyb1Izlu7arfwRkIDxJK8svOPWlavKsejKiRMK9E5v36k24rnwoU9/mfx5Axn6/Csyqlp1SXz8GZnz/Iuy6LW3ZcKLr6i1xppr805NaVunnqSvWVvlRrocU0LogJEFq29f82MTfQaIEyYkqjJFwTIPlFdCYlz3RUO1ezswZ75Qh5CFYk/YPza26qjOJZPMhQOHhPyjt687r3lIxvY2bds7gKEBisnTZ0rHTl3KPdykbGPXuCNmmGyPHixIa7OefUESnqkmhSYa91F5btuWLZTvadOmnZLqeaHyIRgZ131c6snSgZcX4Ighf8uWrcqmgiYAwoPwVPYOOTA7VUhIntm0tayv11Q2te8m+yZOkVNbs0MaJEn6Tbxt69ZtfepRbf8cOvuOOo25CxYvX0Dy8PJVaj2h3WD9obIvKyGVHWvemNGyOyFGDqZOlltXy8fezU12TOtYvOd1efzRn8jlE3eq4Bgq1KK812TVyrthWLdvXJOipTMkb3R/yRkQLaf/k8zdzCoXqDhx1LmT0rKBmi/mDHtg/foNlcq7c+euKv3duLFjZVHKdNm5dJkcWrhUdo+ZIHjtohlaX7+ZqsiREzVQCqfPkpObtyoTR1leheL/eW+Rzxm+UdXFTPIZIK5evUZ12Cq6d3eYdunwEcmflKxUG6gjV3/wqQKl7D79lUoRlY+RUYa/K9/9UPan3FWz2N8LJ5NevaNskqEBhsbfjp26umQHaNOmrdT++BMpPX7CvnmPv+/atUvat++o5ojFGBMzWJXuKS4uFrL7m0HYNFHT4l2IOmnNh59L3sgxwu64st20Gff2dxtGIvO0tMDFg6F+ZyODnSXYiNhbMjPxgl/179qy/PV3ZH3dJsruiBS5c1icbB8wVDa27SjbojqoQPoTG1cLCTfsiQ1oWqpjWsfWzf8qTz52v/z2N/fJvuxXHaTEqZMTbZdfOLRPDqfNk5KMVXLl1CmhqgcvXWyz3hL9ovILpiTaRO1HOTE8gknGzkveXcKcgqMgm3dMKGs/+kKZNza27ih7xyfJqeztgjQeikR91bp160nz5i1M3YT6DBC3b7+TYJgXbLATO1YkQdSB6PmxH8z/12PKY64qBwH03D179akQEDt07Cxnz56rlEXozQEtXnZV2RwrbUhEJSpA0mNRsiMl2JVq4P6ga6dLlUSAjXNdncZKMkAtfDnIC5nCOxzJyBaCl2cgY8JwtCAUJNiJDRMbKl74SD/H1mUIdUHP5O6S/CmTZefIaJWF5vx+RzsSa2XJrKcdQA+JMCryz9K7258djl899Y7MTCmfrs3g3apVq9U6IWWiN0R2I+aFNYftD5OEt+u4ov5gasGxjo3ojgFDJaNhc+VbgKapeEmakEoyVIianPAUXwyzyGeAiFqNzh47VmJWXy3Tzo6BQwVPVVeIlyOqNHv7oSEdTp02Q6Vdqqqdo0ePKl527ux5UD794DNgwB1VA3ZBjPWBJHbhOCmkf95AVr5TS6mELpdYu8ROZfwiaJhnftGixZWd5tPf8DIN1ZyWrjJudplaqRtXvygvv/BzuVUmz/HaxdWksrJ0pEpkPvfscW2tO+sfgsGnn34ujRo1EV+l0nN2X/tjqKbP7NipnOZwAEp7833Z3m+Q2mDY23btrwmG72x+2GhgrjBrU+8zQMQhgweBToc74SZN3kd7UMR+2K5dR3Fl90l1CxYmIOYp4aZsxFQhIfpqh+pp/66dOascLfBmXVe3iRQvXl5hCjtP7+GP63r06KVefoHw9CR3sNW9u/0xB3NnjZeLJW/bpMFPav1SPv3oQbl59m4IBuEYI4Z8Xml38Bxm3aHS9ISIXaxXr4HSHFgpjhEHIEKysPHjDYtJAye92zedp7jzZOz+ugYfCObIrLRuPgFEXrZdunQNq6DRyh4AnE9wiW/broOqEtC1a6TyTFy9erVLwITEwaR7GsZiJHs22rAaGJblHerp3fHjZV2dRrK5Uw8pydxYzrmp7DVW+T/SBHzGs9jfRC5R7o19OJyJzUjcsI9toLhr08ty/cydGETDoWbW5JfkSBWxo40bN/XYTIGDEyp0tENkcbIqYYc8sXGz7Bwap0JCWG9HVqULx4OFYmNHKAfOy5eveN1lnwDigQMH1cJMS7vrweV1T0OgAaTloqJiBw9NV4aFMZ4XnSeZfwxnDzxJ+QQbnd2Xr+IiUamSVMCVMlqBHiM2PDPVOK6OBw0CL3FPHDRcvUewnIe0PDqurWSlOQbnH9z5msQOfEsOHqg82X5paalacx06dHR7yNiTcfawSsIGdwZAEQJqsRLuhWnIk5AYd+5nxrnnzp1X9VfNkBJ9AogLFtypcnEihAy4ZkxcZW1U5nWJGg5A3L/fvQwf2dl31Alcf+3aNfWprA9W/k2pVCclK3sjQd6ExFiVyEKEqsyMBerqGAFBnKQoFKvpLgf27t0r8+dOlkWpMbJg3iihGEDZDcP5A3vKeT1v27ZNrbnBg4fcbcyFb2iDCKBnQ7R9+3YXrrDmKTjnHF6xWja0aKccc8hkZWWVKtoz1hybEW/IJ4DIDtmKKaS8YZSvrz25eZtyOS9bFQD1JhONbYhdr6vEwiQ5QseOneXKFe9VCa7e19fnwQ9UOjgHEA96sciaXszY0FGZGbGavuZLcfFh9QLnRR6udOnIAbl2xjUvSp6jW1cuqdSNhxZOk5PbMh3Yht2fTejUqckOx6v6D1oxrlu+3Jrloqrqv7Pfr5WelR2DYpW9kaxXlW3enV3vj2MAYd269b2OwfUJILZr1yGsPd1uXbkihdNnKycRVx8G4rAWPv6MLHrqedk9ZrzNoQQvXRYYAfKuEjvgsWMT1DXnzlUe0uFqm1Y8Dzd8kgBg93Cl3qQ/x4BHL4kNyKrhD8I5i+ckGON+zeDP+YP5kp88Si4fPSS3blxzyTZ/KjtLCmeNl/wpcbI9prtcLrnrEW841BB64SoVFu5Xqrtp01Jcur+r7VrlPDbru0ePV5vRA7NSLScx4iuBAOANmQ6I7IiJmUNtGm50+/p12TsuSRY/86IK3Mf+5Q6R49AI+F/64uuC1Jievk696EaMiHO5KeyGqM9IwxYOdCpnh6pgQnxoVVlO/MkPArq7dYv0yy0TE5MUAPvlZha8yY7B/VSwPgBXtHRmuYB9Z11GSjyQOlVyR/SRDe2a2hL9G1oZNhgE57tCONEQWmFGEL8r9wvkOdjxd8WOkrS3/y2Hl62yDPiT/ABVdX5+gcfsMR0QT548qV7gVHQIJ7p8rETlQzUALe2t990efsG0WQoQqR9nGLPxGMQWRbkXV2jZsuXKmO9N7JQr97HaOahx4B+ZhUghx0st0LRnzx61FnhZ+prIG9ytW3df38ay7W+P6Sd5YwbI3kmxcu3MKZf7idRD0D8ZcKh+A2FvZ83NmDFTeMlWRTjg4EQzZMhQSzx3VfXXrN8poUfCfzbypK+0AhFn7U2mJtMB0fAw9UUmciswvLI+8BImxyBqT1KUuUsk/iVInfJNZMZ3p2YjalIWMdL5hg3hWXkCflOyaufQkSo7R+lOxywm7s6Ht+ejLWnYsLEsXrzE26aqvB7HKbKghCMpQGvXWZU8O7+/0CMW4DBCAm9y77qzmTp8+Iiy1ffq1UcCEXvq0WBNvojsXRRH2D16QsDDNci13LBhI4/j300HROrCoWrIyQleDytPnxdqEJKWDJsWycI9IeoQXrp0WTLnpMr8196StXPnqcwn5IYl5VPZWBsWL9IgCbjhO7GNmkRVAiF3KjkwyyZo9yd/iJHq16+/z29JsWJeyuFGaAbS3vq31ynJ2FAeLiqS6f/+WNZH9lZrjpfrhg0bVZhUWc9UpMjU1PlKRdepUxe3HN5CcY6YB8rNER5V+p8sXph8/O2ZeuDAAfUezPOwCILpgIgzAXpcM/PLBcsDlP5pfaGclDuSHWMjhdr69RkSHz9G7TaNBMAN3/9Axj9dTTrVeEtNMoCHBMgCJNt7r169lScjx/lgVNZ0lwMsUuo84pFK0eNAENIhOWhZF74kbIhIo/7yavXlWNxpmwwr23q6XzX+9u0vVZksKk9ERUUrZxhjHfV9+XWJe7a61K/1oW3dYR+kOj2bm06dOqt3HOe3bNlKTp2yhrrQHb756lwEgYxGLVTlEoQDqvb4k4j1JkOap/hjOiAyeF7WBAiziwoXulBUrDI9uDpedpxI0YMGDVYgx+LCbZiqEwsWLFAxTOQwLT19Wja16yqL6jaW7A0blfqNXJUsTsrSjBs3XqlKyZmoyTkHSB6OCpq8qUjU1KariFB74Uhx5PARt1RnFbVH/krm1tfZY4i34z6upAKsqK/BeHxLl15yvsB1NSn2PkIqIiJa2MAOr/hJk5Jk/fr1gqcoJZqK09fL4lfelIxxiUKdSRzVsE9FR/dXHvT8nwxS4fSOc/X5YCO6tXsf5Q+x5LmX5JYJGWRcvTfn4ZHP+9ET8gkgbt16J6g1nIKED8ya59JuCJUnahhS2/ECQ5qm1iA8o9isM+Ilzk5rRY2aUpIRvvZBZ7xx9djNy5clq3lr2dSxlWyJbCdn8rIdLi0uKpKE+PYyfcJTsml1dVk5/1kZH/eCTJ400it12PHjd2oUzpw50+F+vvhP9+49lHfxkRCoHuIqf6hA7woRp4kUTfYY1h0p1VB5HjtWcXFiqrPgE5DdO1oVAHblPvockYPzFik/CsPBkLJ5/iQSIxC7XVbN7UoffAKIvMCNEkN4v4WDGmdrtz6qRFRFTCc+jA1CnTr11IIkXgaPUHeytFO6hR0xBYMpi6PJPQ5QfHZHbC/ZN3mEFM4cJ+fyc1UDWRmrZNHMag7VEE4felMlhz5/9C0ZH1/LrXmy7xVrAQ9E1Ny+Jl76hNs0atRYdu8O/eeD9UCC6ooI3mdkZEqfPn3VmkN1TXwuGgBXX5ZIOwdTFwvmkPwpKSFd6LoiPnpynCoap7fnKgc3pEQEhrIELmRmZqo5IayMck5mxNFipmDT44mnvU8AkYHzMBo1wI4erXgXVpZJwfr/9fWbVth11FmNGzdRk9S1a3dVpR7+eErUO2NnDDCe3ed5zI0n979544Zs3bJVFi5YIJkZGUGVBedaaakUpIyT/bPwhruT9YcwofnTq9kqI5D8+drpd+V3v/2B5G58WR0nMXTs4IaesEtdgzqOBeqP8IsVK9LUvahsEup0OG2N7ElwnqoOWxK5e+E7NiUyznhT7uzWtevKi5Ji4FSo5/+aXOMA8dl471J43SDAEPMQVYAmTEySyVOnS8zgoWrzWOCGCtxoz/4vQgZz7kn4hc8AkQ4a9cTwPA11ohBnWcKRAuMualFqGaIW9QYI7dunHcIz2CFnNIxQMXioeHxFvGDmzBwuc5OrSX72K3Lu8FtyaNdrsmBGNZk2Ocqrl42v+uysXTxOT2Vnyq1rd+yI48f2lBuljpUQenT+k7SJeMgBJHdkvSjbc3KcNVnlMVzzcZTyh9MTKaxwvArGRO5VMrLMCXvGJkrx0vIFBPbvPyDt23dUzkwEypNb1iy6fvacSr5BmAGJICgajhSpyT0O8F7sFz2wXOH08ROTpGXL1i7Ff1Z2x7i4eKU2dXcT5FNAzMrKUju0/Hz3MrZUNlAr/sZLNqNxC4eukSWGRYnNAtWoWUDocJP//IfYO7wpV9b8SDIat5SDcxdK2Zyozq5zduyr27fk1nVHpxN2c0Nj6svZw285gMTKhc+qsjpIUKOGfVShDdTZfQJ9zJiPSaOfcxhTUd5r8ttf36vGev303fp5SI6pc+I97jYqIVK5uZOP1tObdezYSTlreXp9sFyX03eAnNp6d5PCnJJLlA1o3779BOnfl3RnQ9pVZWzJHTxCiHs1nqvK7ks4wq2rrjkcsvbI+oXKHaAYPHiocvyprH2r/4a6mrUwLWVWOUCkeHqfvv2UI5M34zh06JDCHndz0foUEAEFVBahnrXmxoWLKqiXCWRBrF2brlzg2aX7sxYa9z5fsF/2TZwqWc3byvq6TWRr976SPylZStZnqQ/JsFFhlKUjK1fLzpFRsnfycDlf6FhPb8G8RDl1sIYDcCyY8aT89c/32mrOAZaTE2PKNmvp//OySUl8zGFcU8c/Ko3r/0F+9MP/lrwtd1SmRg29+XM8j/PDbICUSF1MXxPPnbc5HX3dRzPa39yhu81kQOgDIRGER1AmzRVgMqMPtEHB3SNpaySnzwBZV6exZDRpqWxnRQuXSumu3UJybPv+kJx+4eNPqzV6YM58ISGHM0LaJ7SqX/QAmZKcogqMj02YoEKzCNMKVkJqa9a8hVMwBBCHDR8pJMf3lvAKZnPkzjvYp4DIQ4CRn1ifUCbSFm3p3FPFmiEJ4OFkhaoD8J9iuyxMstXzIRk2wbNUyibN2ap/fyKrP/hUslo3kV3x/aR4+Wy5fdMxZm72VEcwvHWupvzrn/8j27NedACTlIlPyYXzF4JqqpPGPOMwBsBvaP+/S6+ufy53fN7sYV6NDamFChiuOnR4erPly5erF0FFXsuetmu16whHYgNIyBExmFTZ8bb8jxljJFH/yc1bVaA674X02vXVWmPNsfZWf1BbFjzypEPe4jO5jptQ+kEd1IExQ8oBR/L0mdK4STMVHmJGf/3dRlUSYt+oaOUM5W2/du7cqQQywmNcJZ8CIp0YODBGRowY6Wp/gvK86+cvyKLPG0iHDp3UztzXqhozmXTz8hW1wy1auFiy+3dVCZJLMpYLqlPo7NmzsmrBsw7gMD3xMfnso9/Jzk0vSmnxXTXqgdxXJX3tWjO75/O2EuKbOYwtY0V1efjv9wnepYZkyN8T+2vI6lXl7VXudJBdLxqTwkLX4+bcad8414hJ9PV9jPsF6u+m9l1l6tBYZTM1Q6LwxzioRH/l5ElZ9GQ1WVr9NVX02pnGBtAglnv6jNnlABEpKqrfAL+kBPQVT3B4GTQkttzYkIQBezO8TdEAoTEg1ttV8jkg4gLbrJnvd8WuDtjs83hwVyxbLvHPVFdB8sGcz/D6ubPlpENSIW1Nf8EBHIbHPCwf/fsB+ftf75XLJ96x/Xa2+C1ZMH++2Sz2aXs7tm+Xjauq28YQ3esvUv+LB+Viyd1xAYjxw2p4nW0GgAIQCQz3JWGnxHbt6/v4cgxVtV1Sclym164jvRs3ETYAwUSoTAumzRTKxFVE1DBt2iyiHGAAhnxih8e5XauxonsF4viVK1eld+8+Et0/RpKmTBOk3hEj46VlyzayefMW07qE+aBBg4Yuh/75HBCPHz+uXgKe5pYzjTM+aAjdNEmVsQ2l1WnggzsEvkmy/S+fW16t+NYbD0juhpdsQAJo7NnysmzcsDHwnXazBzOmDVces4zBHuD5P5+ZSc8JGwNvCa9jANEfMYlki8IJw9525W3/rXI9mWMwSwys8bYcXrnGKt0ytR/MW7NmEQooDBC0/9urT5RyIDL1pn5ujPWAZN+zZy/BEYxQoVOnTpvai7S0O2FIrpbx8jkgMjpUiaQZCxXiYV23bp0KgG7RopUKgl7foHzYRaiMd+aUjx2Ar0u7P8sHNX8tt887emFOHPW8X7wozeYr87l08WyZnvimZK9/QamBj+W/ISvmPSvTEhtKQYF5sZ6owbAj+pp4PgFf7CihQqRUIwsJ4+LlWTArVQqSfZ8BKFD8I1C9a7fIcmpT4vbItONKaapA9d0q992+fbt6XrZtc8xMVVH//AKIxF/xELuK0hV11grHUZGSP5TxdO8eaTPiZ0W0E+xxoUjpa1ZK3ua7HpfNGj4ky+c5So3YDxfOnxP0w9+3d68sXrRIVq9aJeS9NJvatGnnF0Bk992gQSMFIGaPIRDtkWMUV320MYaTBGEO26Ndtw8Fot/e3JONGsDftl0HiR8zTiZOmiIxg4Yqyf/QoSJvmg6baylDyLvaUoCIkwkPMrkWAZRgJbz2yDt6Bwx72MCQ8ewYNFxObnVtFxKM408YHSX7d7yqJMXb5xwlw5MHasjYUe2CcVh+7zMaBaol+INSUmaoZ3XvXteqvvujT57cg7JnpF3DhX7z5s22JsgWk9mkle3/ofply5atEhs7XCjxZXaigVDlmTEuAxAp1u0K+UVCpCNG1ppg8QYryzx23ORlBQxx7+b/9lS8ZIVK7WR/LNS+r1wxX2ZMrilpqc/I9swXJX1JNZmRVEMWzJsaVEH5gZwXpBwqq/uDSCSPehYQNsNrzx99LnsPdvZ169ZT4VsAY1ki1hZJSpPmgDMOGCpTV7MV+Q0QSf2F4ZQcc97mqnM2cF8eQzIkyBMwRF3qTMq9fPy4SgDsy35YpW1yBZJ96JSPM4FYZbxm9QObD88Qie/9RRkZGeqeQ4d6F0Ppr/7a3wf7J2no2ERU5rKynAAAIABJREFUVJUiu1e0nN0T3BKw/Zj1d3M5QJQDz5Czd7azO/kNELk5RnEy/+NY4I9Ex84G7O4xYlnY0Rsvssp2o6ve/0RlpXD3Hvr88OCAURZt6dJlfhswzyu1+3h+PS2a6rfO2t2I/K8E2+NpWVk5q8PLVsqu4Z6n1LO7pf4aghzo2zdKIiJaujwyvwIivcrPL1CITT1ApEar0+zZc9TLhLIxlYEh49g9ZqLKI2r1Men+BYYDM2bMVM+Sr4sFlx0dz+2QIcPUvSm3Y3VCA8HGGUCsauN84+IlSXvzvSrXptXHrPtnPgcQZpAOo6KiXW7c74BIzzCMozrFhbqsLc7lnvvhRGIn6Wd8/GiXAjupz0YOUU2aA844ACASMH/VxcTOztrw9BjqWkqPYY+zspMNm2Ti0nCgcbWmIyncqDoBkWhfk+YAHAAQeX9PmzbdZYYEBBDp3ZIlS9WOdcyYBJc7688TyYlIDUMC72Gsq5TZrI1cLjnu6un6vDDiAM8S3oKBIp5papRSpJowBivS/Pnz1XvBneTVp7fvlK2RfeV8QaEUL02z4rB0nwLAARyyMBW4urGiiwEDRNQ4SF502Iqep7g3s0vFluEOFS9eLtv7h04SAnfGrs+tmAOoSXnWc3LuSDIVn+nbX8i4g5TqijrStz0p3/q5c+dVBhoc19wl6hNSaSIvboy7l+rzQ5ADONFQ8QUbtDsCTcAAkTmgo3i/8aJAnVSVjc5f83b06FH10li92v20ULdv3pS1n9WWXXHDpHRPnr+6rO9jYQ6QW5SAfDI2WeEZ37UrTwFP27bt3CqN40sWwxeSMPMSc8eMcmH/QaE495LnX1bVIza19U+Mpy95odv2ngOrVq1WuLJ8+Qq3GgsoINJTHn6SvAKKqanWKBPVu3dfiY4e4BYjjZMvFh+WLT3aSt6Y/rIncZjcuu5aIVDjev039DhACi6eb9KpWYVwbkN1ivOKO/XifNX/3NxcxSN31FtGX7Ddr/p3bQWIFMnWFN4coMACoTpt2rR1SzqEawEHRDpx7do1W2gDruGB3EWj0sIzqbi42KOnigz2+cmj5dDi6XI8M01Kd4du9hqPGBRmF5Hcvm7d+tK1azfLJS/IyclRmpB27ToEHBR79uwtVCbwlFh3OX0HyMLHn5Wvgjgblqfj19fd4QDYMXLkKJUZDU2Iu2QJQKTTqE+N4Hfsd4EiKlQnJk7y6vYXigrlTO4mr9rQFwc/B/CYpLIFtmjU8FYkbJvYFJs0aSqAdyCITWi9evUdUiF62o/jGzbJ9XPBVaTa07Hq68pzICMjU2kaRo8eW/5HF45YBhDpK5Iinniol6jl5m9J0XA4wLjvLd26csnbJvT1Qc4BI5eo1esS8hIh1zAeqL5IaF7VNLLmk5L8l72nqv7o34OTA6WlZ1VBYBzGzp0759EgLAWIjABQjIkZrECRqsrueAh5xAG7i6gh54mHm10Ttq8EF6fOTZJ5MyIkdcbHkjqjvsydNarKQGNbA/pLhRzAg4wchajYSPDQrl17iYzsKSNHxinnrJUrVwo2KaSyQBVsJhsNIEP+W1L/WZ1IIE2MIsW8SYjsLyLpMiaKqgLw/dUffR/nHDh79qzMnTtPoqL6Sdu27ZXzExo9qnEsWLBQNmzYKBTA5rxbt/xfwAE7OP3CPEFaSU/JcoDIQHiB8HJDUoyO7m+KKqUqBlGaivsdPHiwqlOr/J0K3hNGvShXTzlWXb9Z+q7MnPS0LF1sDeehKgdisRPQGGzcuEmweTFXVEAAcBISxgkpmnAS4XjZD4uExUK1AECU7EPr1q1XWZOIzTOTAGtyldKHPn2igiIbkzF+wAkeAlDuxAEa13vyl+oxw4eP9OTSSq+5ePGinDhx0uUclpU2FsY/ojEg9R9qdZ5p1tGYMWPV+5kC1GXXGv/n+SHcgUQQpL1kPbBBzM7OkZKS46bPyalTp5QXNxvQrVu3ejVblgRERsTLLzV1vmIunnA7duR6NdCqLuZlGRnZo6rTqvwdZ5wZk553KKh79vBb8tYb98vpQ2+q41vTX5D16TqAuEpm/ucEpLxly5bbgBDVHrGrZcEMMGJxsENkx0ryB8J5unTpJngO9+nTV32PiGihpDdjMZNbF9sx6fm4BmDwRLI8evSY2kHTLgWxg7GAK2Po0aOn4s/EiYly5YrvanwiFZJJBMnCLMJ7lnlu176jdO7STdq0ba9eyGieNLnGAd69rKG4uFFqfnieY2NHqGNlzVjw9fDhw0pjQ6gD64fwGdTgbAj52759ByW5GesNmzphSKwRzAqsVZ6Fsm1X1Vts9CTvrl+/oYqrRWvkLVkWEI2BkWaK8jUwE3WmL14yBw4cVO2vXZtu3Nbjv9Mmd3EAwy8vvCeRnf4ojeo+6HA8JfENtx8AjzsVpBcy19jfkPCYfx584ooAPm+J2EDKCaHy6dy5qxCTR6kkY9Gyy2UxI01WVYyVpPWABztUFjt9dHdxezseM6+3l3JbtWotnnjrudIfNiDt23d05VSXzuEl3qRpc0lMmiIpM+fYPv0HDpIRI+JcaiOcT+KZRfDguTfWAeBWVOSZx709L2kbpy02tv37D1D2aja2rBfjXkiVSJ9ZWRsq3ZCiQeQc1izX8gxVVA3Fvg+ufLc8IDIIdutxcfFq8Ng4CJg3MzE4tkp2qmaA7dypjtLhuSNvy5//9GNBSgQcjQ/1BHfn+c9W48rDYJVzMIijpjGAkLmZOjXZlPmpaozYflEXosoz1EQsOna0kyYlyerVq2Xz5i1K5TpnzlylqgUIOYcXycGDh6q6RdD8jvrJ2CRgwjDTzocqjpch82oWwf8xCRNsQGgPikiMVvX0NWv83rTD5pAcsgY4UcQaO7yvCR8RNjJIiqhjjftjDhk8eIgQwwv4kf96xYo0pcUhxpDzeD+gRTQTC4ICEI1JQVocOPBOkV6YQnjEpk2blSHXOMfdv7x8YSziu7d05vRpWbu4mg30AL+Ixr+TJx69T+rU/qUcy3/D9tuVk+/InFmzvb1lyFzPDnL//v1KYjNsga1btxFAB0N9IIiFhj2YZw6J0Visxl+AmrAKgJLzQpEuXbqs5gC1MsBPvOCiRYsEj2xvXkRjx45T/KTmoRnEZhbwtgdB+++Dh8TKokWLzbhVyLSBnRU1J3VqeabZALIR3LFjR8A0HCdPnlRlylq2vKMVNNaa8Zf3PlIrjnP032zSgKgB0exnyqP2NCB6xDafX6QB0ecsDtgNNCCWZ31QAaLRfYy47DANSYLdQ6NGTaR790gV3M8Ogg9qHrwKse8Q7I89iF1iWtpKwV64fv165T7M9Z7kLTX6Y/xVEuISRwnx3JF3ZOSgf0iT+o42RC0hily4cEGpPFCVYB9kHurVa6DUpWbYLYx5MePvxYuX5MyZUiWtIrHSdzNsmWb0zR9tYHNl3XTv3sMmKSMhM3c4TxhrjtzErDvs/XgX4tSEWgsvQ6QR6jFSjoe5RmXqTt7SysaJqjuiRasKJcRhsSNV2EBlbYT6bzyvJGJgrgy/DOYBJ6Q1a9b4xSThKo/pKw5yrDXjwzPoawpKQDSYwmJCVYXHIZUzmFgSKOMO3Lp1W2VsxeBa2ccw6uIY4S0h5cyZ+oJNLYrK9OCu1+Sl6j+TiyVvOxzftu4FIdQj3AjbEZsRFiUvVBYkqrhevfqo4/546MON52aPl7WCnRX7H/OIez02Vl6yhMRUtt74Da9x5p2YNrOItde2XXtJTJrqFBS7dotU3sNm3S9Y2sEBZc+evUogsFdDIkCwYSkuPhwsQ/FLP4MaEM3gEDkm8W4yi1KS+8mtczVt4FfjtZ/LZx/92uEYQJk88d9m3dJS7aCGYVNCsPyIESOVtEcWEuoA2hvNsVcQTE+wLztBTeHDAWJJAUQ0NmYSTkDNmreQqdNmOIDiiJHx0q9ffzNvZam24CObDCrDT5iQqDYqeE/fic2tq3gNvzln1Kh42bBhg1f2X0sN3uTOhD0gsmvCDdgsYvc8Pu4VuXn2DiiuWfysXDjmKB2mzX9OsrdtMeuWlmrHKLvCDhSHE1yriQMkQBe1NSqzI0eO+jUDkaUYpDsj6enr1EuaQG2zCSe7Dh07S/fIXtK7D1lVOihHETPSMZrdVzPau337S2U6QsvCekNSRwpH40I+T+L08vLyBFuwpqo5ENaASFApO6ekpMlVc8qNM/CUihv6vpQU3PUqRSo8d+QtGT2sumzdstmN1oLnVNShqKobNGjk04Du4OGI7qkzDmBXZN2dOXPG2c9eH0NNWFBQoGLqsPWGMmEu8oW0Hco8q2xsYQ2IpGnjYXK3iGRlDDV+w76ZOne2JCa0kqkT3pekcY1kZkpSpQGnxrXB9BfbDdlN2ImSVBd+kttQk+ZARRzAtb9hw0YV/ayPV8EBsgeRlQX1J+uNcBNfhCBU0Y2Q/DmsAZEAax4o6sJpqpwD7LpLSkpUPkI8BklmQGFnwzsUPqKuIXgWkNSkOVARB1ChE0SvqXIOIN0WFu4XqpGQsWno0Fi1xoxEEMTGDhs2POC1LCsfRXD9qgHx49rKrhVc02ZOb3FtZmdJUDOLj0weeKTh9ADokS0GL0KArmxgOhmDSKzNOdgNtbeaOXMSDq0AiDhYhSPhRY09EykPZzK0VKRLW7NmrXIyIkyMzQJmBzaZxgePbEJeSG1GGAubeDMya4XjHFQ25rAGRCN9FA9gqEs1AB6fjIwMZTMlLZJRNcJYdPZ/AUAWIZkhMM5jq8DVnjgmbKShzq/KFo3+zTsOYLPnWfNXRQ3veuv51ayRffvyFYCRKBuNCinJ7NeZ/Xd+q137U+UdSiqzlStXqTSBZHDS4Uiez4M7V4Y1IMIoks3yUJoRmO8O4/11LlIfu07AzYj7Y7yEmlABAjXMgAExygMUHlALD09ZpEc+qEo1aQ6YyQGC6NEwYEcsW7HEzPsEqi3WDbG2eFgbgEeYEbU7Bw6MUVoX4qZJEkKoCFKiUVVEr7dAzdqd+4Y9IMKG4cNHSP36DeTUKe+D8wM7nXfvjgct6kzD3sBi5LNs2TId93eXTfpbgDhALlQkIvLAhhKRJJvqKQAhmbTI2EO1EG/yvoYSf6w+Fg2IImp3RgwPSW6DfceKagX1JnGALEgCdVG5oL7Rak6rL8fw6l9OznaVUJrY1GAnwjz69x+ogLBDh44q+bSuwRh8s6oB8T9zVlp6VtnLgq3Kuf0jhzcaNj+kQozvjEmT5oCVOTBvXqpy2DKjuGsgxokfAlmYkAhRkeK5rjeegZgJc+6pAdGOj+jyKQWFNxfFLIOBqCeWmZklpKADCLEJVlXQNhjGpfsYHhzg+cXB64sv6qiQHexvwUAknMY5CNtgRERLWbFihWmJyoNh/KHaRw2IZWYW92d2e/Xq1Zfly5db2qmEEAlUvfS3b99+SjVaZjj6v5oDlucASSzat++gnuPo6AGW3oxevnxF5QqlKACJKHCMMatih+UnKgw6qAGxzCSzQwUUKQqLxEVCXLLCkyORwHQr7GALCwtV/1DRpKTMlMOHj5QZhf6v5kBwcQCNzMSJk5S0RcgPlTCWLFkq+fkFlsjDiRco1dt5H5C0Ho/sq1evBheTdW+r5IAGxEpYZHjCIYEZHzzj+vWLVo4r/larYhPEI5a+kLnCCuBcCfv0T5oDbnMAaQsVqrHejL9t2rQVKjiQENyfEhlrjGTk2OaRYqnFqil0OaABsYq5xY1627Zslb2FHWtCwjiVucVYqD179pJVq1b5NHAWbzUSIgPGzZtHqOz1VXRb/6w5ELQcwFFl06ZNKkaPeD6efSrS8Pyz7vCeBhx9DU7Y4nv27K3uSfiEDp0I2kfK5Y5rQHSZVY4nkq2F9Gbdu0eqBdOoUWNVCdzMDP54q5EWjSB61Ei8BHQSX8d50P8LHw4gGeKNOnFioi2RPLbzLVu2mqotOX78hCpXhsmEtIWkVtMUHhzQgGjCPBcVFansE+xe8TqbNWu2KWodEmjTJnaLgoJCE3qqm9AcCA0OEG9LhRXibVkjkZE9hHXoLdGG0ebIkaO0ndBbhgbZ9RoQTZyw/Px8iY7urxYodQGpTO2JnY/FTu5DFjpp1c6dO2diL3VTmgOhwwESXM+dO08lw0aiY92cOHHCowHm5uZKvXoNVFskuNcUfhzQgOiDOc/OzpaWLVsrQOvQoZNbKhdUrsRBAoaoSHWQrw8mSDcZchwgcT3JKABFtDTTp6e4nHWKTSvSJrl+kQ7NkDRDjsFhMiANiD6aaBxh8AQl0B9w69YtUtLSVqqKE85uScJjVK04DGAvpMSLJ9Kls7b1Mc2BcOEAYNa3b5Rac+Qnxt6Yl5fnNJ6Yzea+ffukT5++6nxykJLYXlP4ckADoo/nnvglEmp37dpdLTrAEelxwICBquI1aZ/4zahGMXToMDl06JCPe6Wb1xwIXQ4AdCStoKIECTZYc2xMe/ToqdKsUWmevKONGzdVvxHPS9Ubf4ZzhC73g3tkGhD9OH94pi5fvkKpdnr16qOAMDKyp7J7ELoR7InF/chKfSvNAZc4QDklqk0QukFBazQ1bEABRNSq/KY1MS6xMixO0oAYFtOsB6k5oDmgOaA5UBUHNCBWxSH9u+aA5oDmgOZAWHAgbAER12zCJDx10Q6Lp0MPUnNAc0BzIIw4EFaAiLE9MTFRfv/738vXvvY12+eVV14JoynXQ9UcCBwHKAb84IMPyn//93+r9fe9731P/va3v0nnzp0FG7s9NWvWTP7xj3/YH9LfNQd8yoGwAUQM5x9//LF84xvfkI4dOypXbGqa7d69W5YsWeJTJgdT4yQT+NOf/mTpEjzBxE/dV0cOLF68WAEhpdW2bdsm6enpEhMTI7/4xS/kJz/5iUrNZlzRtWtXef31143/BvXfZ555RpYuXRrUYwiHzocNII4cOVK+/e1vqxCIcJhYT8c4a9Ys9cLydyUPT/urrwsuDhiAWLbXJPR+4okn5M9//nNIJtFGI8XYNVmbA2EDiKhJGzRoUOVsEP7AwgQ877vvPiVN2scnff755zJ+/HiJioqS//3f/5V7771XWrRoYYth+vDDD6VXr14O98nKypJHHnlE9uzZo45TXfvhhx+Wr3/96/LHP/5RUlNTbefHxsZKjx49ZOLEifKjH/1I1qxZo36jX48//rjq1z//+U9ZvXq1anPXrl22a7kPff/mN78pv/71r2XSpEm238aNGyd86tWrpxYmUuB//dd/yWuvvWazo/I76iwWL6oq+qwDlW0s1F9M4EBFgEjTaGt49hYuXKju1Lt3b/n000/V9+joaHnvvfccenD79m15/vnnZdq0aeo4Afg1atSQ73znO/I///M/0rNnT1umJ1S1SGlbtmyRhx56SDp06KCu4fn+6KOP5Lvf/a488MADyqSCqnbKlCm2e9Huiy++KN/61rfk5z//uQwZMsTWLoH9rDlih19++WW1plhbxrolaxXriHFxX77Pnj3b1rb+Yi0OhAUg4jjDA7lgwYJKuU8APbaNQYMGKYcb1Dq/+93v5I033hAWH1S9enX52c9+Jk2aNJF169bJ0KFD1UJhwUKRkZEKyOyLh7Zq1UqpIVHbzpw5U+655x5BEqNfgBbgm5GRoa7nJQBIsqhGjRolR48eVeAHyEVERCg1EwD697//XY2JxQrt2LFDAXhSUpIqZDx//nzVD8YAMSY+vCheeuklAZR5ObFReOqpp9Q5BQUF0qVLF9VucnKyUiWTcUeT5oBZHKgMELkHoISqFKpbt64CEL5TDoo1vHPnTvUb/7BeATLWEcDE5rRPnz5y5MgRYXPI5m7gwIHqfFSzXP/kk09K//79ZfPmzUJ2KNZZtWrV1Fpev369PPfcc/LDH/5Qxo4dq64z2h02bJgyI7Ce2Kga7eKYR7v/93//p0AUIASUf/zjH6scxKRixCTDOaxtvuvUcLYptNyXsABEdmk8kNgsKiOAqE6dOg6nLFq0SF1r6P8BRHZ59lSzZk159tln1aGyu1wceX71q18pcAMQf/nLXyoJ0P567CSNGzdWh1jQ2Dnts9UAyGUdfwA+xmQA4jvvvCNNmza1b1aBNhIgZAAi11D42CAWPsdKSkrUIa0yNTij//qCA1UB4l/+8hdBQoPsAZH///Wvf5U2bdrYuoVkZ0iNaGkAInsaMWKE2iSyBg1AxHRiEGuCtWZfVxHVLRvUhIQEdZqzdrt3725r1wDEqVOnGs0q5yDW1Lx582zH+D9j12RtDoQFIO7fv1+99FeuXFnhbJw6dUqdgxeqPZE0mIcZ9QsEINaqVcv+FGnfvr2wkA1C3WgsanaiqHDYjbJzpS1ACgA0PuxSUbdAACJSqT3df//9Sk1jf8xY4AYgsiNFJWS0yV+AG9UpZA+I9u0Yu1ejHQ2I9tzR383mQFWAiPbFMDmUBUSACJUl2hpSIuKhaqhL0XLwvNs//6wp1htqUWO9sDk2CDD9wx/+YPzX9venP/2pDRCdtYtESbs45RmAWHaz/f3vf9/WBg1zvgZEG4st+yUsAJH0TahB8C6tiJDIeGhRFdoTpZg4zmKEXAFEdpeAGPdFtYpdEdq7d69qq2XLljJ69GiHj2E3ARDtwZXr6Lv9zpZjxgI3gAyVKvaWsu0aO1cNiGoK9D8B5kBlgIg6lLWGfRwqC4g86wawYIdDRWqkO8SE8O6775Z7/lkPqP3Lrhfaf/XVV1XIR1mW2ANiZe3ybtCAWJZ7wf3/sABEpghgYkdpOLaUnTZ2ndgGWrdu7fATdgUW4dy5c9VxVwCR3ev/+3//TzkH/OAHP7A5CVC7DUcWw/7gcKP//McZID799NOCWtaeJkyYoPplACKxXG3btrU/xeG7BkQHduj/BIgDFQEijmtIXjzHqDihsoDIMWyAH3zwgdSuXdtBfQoYltXc2A/RGSBiq8RRxj7+EU0Rm0tDZVpVuxoQ7bkc/N/DBhCxkaGaxKkEZ5WDBw+qUkyFhYXKoM5UYjgHyHBIASA559FHH1WL8ObNm2q2XQFETsQb9Te/+Y2SFNlJGoQNhJ0tEiESJHZFnFkuXryoTnEGiKgx8UhFwkXlw+4YKRKgNgARSRBnnZSUFDHuh9SLqhZyFRBxCqBdAFcXJjZmTf81iwMGIBpOZ4bTCd6iqP0xMRjkDBDx/kQdScwi68Yg1hP2QGIaMXNAAB1mCsgZIPJOQAWLfR4vVM5BzYqJwwBEo11D0qStY8eO2dp1FRDpL5oh1qN2VFNTYsl/wgYQ4T6ggxclkiAvfeNjb4zv16+fWhAsCn5///33bQuMNlwFRMM2h2eZPQG0hFXgzYrxHomRBW68CJwBItcDioRqAOg40ACA9M/+pYATAYDOrpf+IxEb7t+uAiI7dewmtE3fNGkOmMkBAxCNtcdfwAIPantJjXs6A0Q2e5zPOixLkydPVhtQgJH1xTpgTUDOAJHjxcXF8uabb6p3AmFNaIRQmdqHLNEuIVi0i1cr69Zo11VAxC5qjJkNtyZrciCsANGYAiQzFgIPszMpCJUn9j52r74idsjcHxds+uMuAYgsTPpqTwAaUi/SrSEp2v/uynekVhyRCPnQpDlgJgd41nmuCGfgY3YuYZ5dNCNsFMuuDVfGgQTHuiJsw57QEOGdzbrwVMJDWmVtGtom+/b1d2twICwB0Rqsd70XOA7k5OTYLuAlgq3FcNax/aC/aA5oDrjFgY0bN9qSagB0OKbhmY0mR1P4cUADYhDMOflFsSEStEyWGtSh2DZ1FpkgmDzdRUtzgBAnbPr/+te/lDkCO6bh5WrpjuvO+YQDGhB9wlbzG8WQT5aMOXPmqAQDegdrPo91i+HHAaRCskRho1+7dq1Hatbw41rojlgDYujOrR6Z5oDmgOaA5oAbHNCA6Aaz9KmaA5oDmgOaA6HLAQ2IoTu3emSaA5oDmgOaA25wQAOiG8zSp2oOaA5oDmgOhC4HNCCG7tzqkWkOaA5oDmgOuMEBDYhuMEufqjmgOaA5oDkQuhzQgBi6c6tHpjmgOaA5oDngBgc0ILrBLH2q5oDmgOaA5kDockADYujOrR6Z5oDmgOaA5oAbHNCA6Aaz9KmaA5oDmgOaA6HLAQ2IoTu3emSaA5oDmgOaA25wQAOiG8zSp2oOaA5oDmgOhC4HNCCG7tzqkWkOaA5oDmgOuMEBDYhuMEufqjmgOaA5oDkQuhzQgBi6c6tHpjmgOaA5oDngBgc0ILrBLH2q5oDmgOaA5kDockADYujOrR6Z5oDmgOaA5oAbHNCA6AazwuHUK1euyOnTZ+TGjZvhMFw9Rs0BzQHNARsHNCDaWBF+X27duiUbNm6S8eMnSu8+UdK6TXtp2aqt+uzfvz/8GKJHrDmgORDWHNCAGKbTf/ToMenTN9oGgAYQGn9ztu8IU87oYWsOaA6EKwc0IIbhzJ85c0Y6d+lmA8MBAwfJkqXLJDs7RwoL98uxYyVy7dq1MOSMHrLmgOZAOHNAA2IYzv6YseMUGLZp20G2bcsOQw7oIWsOaA5oDpTngAbE8jwJ6SOF+/fbJMP09HUhPVY9OM0BzQHNAXc4oAHRHW6FwLnxo8cqQBw4cLB8+eWXITAiPQTNAc0BzQFzOKAB0Rw+Wr6Vc+fOCR/DaUarSi0/ZbqDmgOaA37mgAZEPzM8ULfbuGmz8AEQO3bqIjdv6jjDQM2Fvq/mgOaANTmgAdGa82J6r1JmzBI+AOLoMQmmt68b1BzQHNAcCHYOaEAM9hl0sf+xw0cKHwBxwYJFLl6lT9Mc0BzQHAgfDmhADJO57h7ZS/gAiBs3bgqTUethag5oDmgOuM6BgABiUVGRNG7c2Onn/Pnzrvden+kSB/AmbdW6nfoAiPv25bt0nT5Jc0BzQHMgnDgQEEDctm2bfO1rX5Pnn39eatWq5fDBE1KTuRy4ePGSkgwBQz5Hjh419wa6Nc0BzQHNgRDgQEABcc2aNSHAQusP4eTJkw6ASDX4yTZxAAAgAElEQVQLTZoDmgOaA5oDjhzQgOjIj5D8H7lJDemQv1oKD8lp1oPSHNAc8JIDGhC9ZGAwXF58+LADIKJC1aQ5oDmgOaA54MiBgALihx9+KO3bt7d9+vXr59g7/T9TOKAB0RQ26kY0BzQHQpwDAQXEsk41zZo1C3F2B2Z4x48fd5AQS0vPBqYj+q5BwQEc3mbNmuVyX/Py8uTsWWs9U5Qvy852vZILZgTGYU8PPfSQREZG2h/S30OcAwEFRO1U45+nq6xTzZkzpf65sb5LUHLAXUDkfDzHrUQJCQnyl7/8xeUuEQaGx7s9jRs3TjZu3Gh/SH8PcQ5oQAzxCWZ4eJXaO9WcOnUqDEath+gpB+wB8eDBg7Jy5UopLS0VQKZ///42yYuYYSRJzh84cKD6bsQRI3ElJiZKVFSULFy4UL766ivVHXLocg1/J0+eLCtWrFDfOcY5ixYtkujoaElOTnYoUr106VI5cuSILFmyROLi4mxDy8jIUH0aM2aMTUpF0gPgfvGLX6h70X+opKRExo8fL71795bp06fLrVu31PENGzbIyy+/LE8++aQ635AUuc74zomcP2/ePNW/0aNHq/6oBkSETefy5cvlypUrqu8DBgyQzMxM42f9N0g4oAExSCbKm27ykrIHxJKS4940p68NcQ7YA2JKSoo8+OCD8uijjyqQeeONN+Tb3/62kpwAqI8//lgB4quvvqq+c4zn7U9/+pPUr19fRo4cKX/961+ldevWimsAIe03atRItQuwXL9+XR379NNPpUaNGtKuXTu577775M0337Rx+h//+Id88MEH6jjnQQDnH//4Rxk+fLh88skn8stf/lLY7M2ePVsee+wx+cEPfqD61LVrV3X+559/Lm3atFEAynW1a9dWx2NjY+V3v/udup7xcD30yCOPSK9evdT3ixcvKh48/PDD6r4fffSR3HvvvTbJeNOmTQqAuYZ2a9asqcZktKUa0f9YngOWBEQWDQ43PPCavOeABkTveRhOLdgDoiEBzp8/38aCl156SQxQ4iDn26tMO3bsqNavcUFubq584xvfUBKacT7Axzo3iDYAPEOS3LNnjwLedevuFLEGaH7961/LmTN3YmgPHTok3/zmN2X37t2qCbIxAVY9evRQ/69KZbp69Wr5+te/rgCUC5ypTO0BMSYmRn7+85/LhQsXjC7LF198oUCTA0aykaSkJNvvgOJbb71l+7/+Yn0OBBQQWQRlPzzIxcXF6vjf/vY363MwCHqoATEIJslCXWRNAoSQAYiGepFj3bp1UxKY0WXOtwfEX/3qVxIREaGu5Xo+SJWoQyHOR1VqTxwDpOwJp7suXbqoQ4AT0p1BvCcAKKN9/iK9Pvfcc+oUZ4AIACPxtmzZUkmi3NNwvKkKEOlLWQcb+ksbx44dswGiPcij2tXvMGPGguNvQADx8uXL6gFiEZX9GPattLQ02blzZ3Bw0eK91IBo8QmyWPd4yQMwEH/5vz2hRgSgDOJ3e0Dk/6+//nq5XMU5OTnqkrLnc9DZsffff18Mz3N7aY3zsQPef//95e6BLRMqC4gA1VNPPSXVqlVTdsAFCxY43LMqQATYUO/aE5Iv/d67d68NEO1/L9sH+9/0d2tywPFJt2Yfda+85IAGRC8ZGGaX85L3BhAJVxgxYkSFXHMGfhxbtmyZwzXYHrFBQmUBEYcdJER7ydX+4rJghETKPYzwkIKCArcA8e2335Y6derY30KmTp0q3/nOd5QjDRsC2rensn2w/01/tyYHHGfQmn3UvfKSAxoQvWRgmF3Oi90dQPzxj39sk56wAQJiDzzwgGAHNAi1okG0by9RcpxjgB42QqMNnGIMm2FZQDx9+rT88Ic/lO7du8vt27dV04AdXp5Qamqq3HPPPbb2Vq1ape6xa9cu1T7qV/t+oA5FCkSSNOyY9vekPeygBmjjPPT73/9eOnTooO6nAVGxIej/0YAY9FNY9QA0IFbNI33GXQ64C4h4kHLN9773PQGocHBp0aKFfOtb31KOMD/96U/l8ccft93AHoiMgxxr27atcO7PfvYz+dGPfiT2jjz24GRcQ8jFb37zG3Uuf7///e/b4gYvXbqkPEf/67/+S1C9Ikm+8sor6hzAGscf+meEZOzbt095jX73u98VQjigsvfEsYZ7YCPlvAYNGtgcgzQgGrMS3H81IAb3/LnUew2ILrFJn+QhB5CoiFfE89OQrmgKj0zsa8T/VUUGSBLegDqTUAxXCOmQex84cKDcNVevXlVSqpHMnr5xHqDtjIz+kuWmIkICpX9GmxWdp48HJwc0IAbnvLnVaw2IbrFLnxwADhiAGIBb61tqDtg4oAHRxorQ/aIBMXTnNlRGRpo1+6wwoTIuPY7g4oAGxOCaL496qwHRI7bpizQHNAfCjAMaEMNgwjUghsEk6yFqDmgOeM0BDYhes9D6DWhAtP4c6R5qDmgOBJ4DGhADPwc+74EGRJ+zWN9Ac0BzIAQ4oAExBCaxqiFoQKyKQ/p3TzhA7N2//vUvFZs3ePBgT5rQ12gOWIoDGhAtNR2+6YwGRN/wNdxbJZF2vXr1VKYY+6TW4c4XY/zENeqNgsGN4PirATE45smrXmpA9Ip9+uIKOEBmmvT09Ap+1YedZa/RXLE2BzQgWnt+TOmdBkRT2Bh2jZB/lDJOVImoVauWZGVlKR6QFo3qEATTv/fee7YSTWUZRHV5ahxyPfUTT5w4oU6hujxFgKk8QV1EI28qP3KPfv36qWwwFNp95plnbHlSjfsaNRKN+1EAeMmSJeq/ZKKpW7euPPHEE6ooMddAZL6hz8Q60heKAkPkSjXGyF8q7Bglp/idlG8UIKb8ExU8yIlqEMV/SfBNWSsKJFN+yugbYAhv4BH3NQoNG9fqv9bkgAbEKubl8OHDKt/hwoULVdmYadOmCZ/k5GRVYoaM9v+/vfMAj+o68/7mS7KbzWaT3U2ySTbFjuOW4o4NppreTDNgwHRsbNzA2GDcwKaaYNN7Eb2b3ptASEKogIQakmZURp3ee/H/e/5HnOHOaEbSSNP1nuc5uu3cc85979X5zXvK++pI/2fffvutw8h/cv5TVCTyH5INhrsi7UrWqPGCNQ4YMMBtedMockWeiWloC9KZfGg/UsuR2wULFig5U9YrV64E5c+Gh/YraStTgmclQOPVtCtKD/K7du1SPhDpkHfdunXKgDa7AtnY044p35d90HY/+X0QjHTXRLNnNKf217/+VYGQ75NeMahpaq8Wa9asUbZMmYbXeD/LXb9+vSqCoKNtUh0sFouqB0FKE26/+tWvFFB5TOe8BC4DvxnWl3Cm415+T+zSfOyxx0CHxzt27MCoUaPwxz/+UdVP50+bp88//7zy1chvkjZWafeUYcSIEXj88cdRt25dLFq0SOVL7xe8zqhtvFJW9M4hwf8lIEAs5x3RIHDvPm+gS5fueLVLd7Rq1RbNm7dWsUmTFmjYsElJbNQUDRs1Rf0GjVCnTn3UqWuMDVC3XgN17vnna0LFF2rhBR1r1sILTuKTTz6FRx97HI8+9hgefbQkPvLIo9Dx4UcehY5/efgRWONfHsFf/vKwig899DAeeugv1vjnP/8Fxvjgnx/Cgw/axgce/DN05LUHHnwIf3rgQfzpTw/ib39/Ak888TT+8cRT1vj3J56Cs/jkU8/gmWefxzPPPY9nn3vhfqxRE889Xwsv1q5vjXXqNEC9eg3RqFFzFZs2bYUO7TujW9ee6NmzLx5//K+q4SvntcnlKkpg4MCBSsvSniSY3VtvvWXjcomAoSZkH4qLixXk9u/fb38JY8eOxZNPPmnjtmnevHn49a9/rQxla/+L0dHR1nup0bVu3Vod0yM9DXafPXtWHVN7IzwJvC5duiiA6xtZD2MduW904TRy5Ej85je/Ae2n6kAtkVZzGHj/L37xC7XV1/ljkhofA3/k8bq2a0obqPTAQVAySJepEkNA/REglvO6+AvTkncC6eZ8pJnycTwjrySm5yLDVABzZhEy7WJWdhGyLSeQnXsCOXknYck/idz8UygoPI2i4jM4cfIsTp2+gLPnLuPc+Su4cPEqLly6iosqXsPFS9fUubMXruDsecbLOHXmIk6evoDik+dQdOIsCopOI6/gJHJyi5BlKURWTgHMWfkwZeUhIzMXaaYcHM/IRmp6FlKOZyIxKUPFhMR0MMYnHMeRoymIiU3E4ehjiIpOQNThBBw6HI+IQ0cRfugoDoTHIvRgDELDolXcuz8Ku0MPYde+SBV37D4Ixm07w1Tcsn0/tuwwxO37sXVHGLZsP6C2TLd9F+8Jx849Edi97xD27j+s4v6DsQiLOIpDh48hJi4ZcUeP42iCCckpFqRn8PlOoKDoDHr26qM8lJfz2uRyFSVQo0YNpc0bszlw4IACTGZmpjpthI0xHb3S0yOEo8CuR+0ySV8nkJhXampquQ6J2fVJLW3p0qXq9oYNG2LGjBlqnzBiNyyBpeOPf/xjhISEqOssgz07OrAuhK0xUNvVQKS2yrJ0XtxSQ2VkIBDpEcMYeMx0DAJEo2QCY1+AWM57IhAJsczsIpizCu/H7CJ1Liu7GBbLCeTmnrTGfAU/ArAkFhefwalT53GGEDx7CZcuXcPly9dx/fot3Lh5G3fu3lW/cPkr9/adu7h567aK12/exrUbN3H52g1cvnoD5y9dw5lzhOMFBccTp86j8MQZ5BedRn7hKVjyipGdW4RMwjH7PhzTzTkKjCVwNCM51YSklAwV4xNSwUg4xsYllQAyNhFRMYmIiIrHwcgjKh6IiMN+AjIsGvsOEGJR2HMPkDv3RmDnnnAFx+27CL37UUOTW6YhTPeEEoRR2HcgWkH3QHgcDkYShgmIiUtSMDyWmI3EpBz1QyQzpwh5Badx8tQF9OrdV4BYzjfrjstcTjFp0iSbrDiBhlDhOB2DMyBSi/vzn/9sc68+qF27NoYNG6YP1ZbjfMyL/hO1hmhMYA8eApXdptTg6IuRXZ8MhN/bb7+t8mA+OtILBwPL4DkdOD7Zo0cPfai27LrXQGQ35wMPPGDNR+cXGhqq0trXiycFiDbiDLgDAWI5r0wDMSf3BLIsxdbI49y8Uyrm559GfsEZa6QGd/LUeWskDE+fvqCAeOH8FVxWQLymgHjz5m3cun0Ht+9F7msgcnvtxi1cvX4TV6/fwJVr13H+0lWcv3gVp89dUlCk1qjBWFh8BrnUGvOKVdRQpNaYlpGt4vH0LDCmpmUiOYVaYxqOJR63xrijSYiOKwFi5OEERFJjjDqqoBhmgKLWGglHao2MBF5ZkTBkuhIYHsb+sBgFQsKQGmnskWQV446mgkBMM+UiIzMflryTKCg6q34M9O7TT4BYzjfrjsvsWqT2ZQwfffSR6trU3ajOgEiw/eAHP1ATVIz3c5/j4kbfiDzHMWJqYsyX0GG+xmAPHmqo9Es4fvx4G6BRc+MyEGfBHogEK7tMtVNh3scJORqIx44dU12/dDzsKNjXi2mMQIyPj1fPouXlKA85518SsP3y/KtuflGb2NhY5TB04sSJalII/wkYObON/1A6vvPOOzDGfv36qfEKNizGyEkKnLFXVuTsNILY3yN/pZf1HLxmfHa9z0bHKCvuv//++1ZZUqaUr5Y1JztMmDABfAccb7py5YpffBvBXonIyEg1S/Tpp59GmzZt1KQm4zPz/Wpt0Xie+xxX++yzz9SEE2qbnG3K2aUMR48eRbdu3cB8mzRpot6phkZUVJT6poz5EZL2WuXq1atVHva+C+nwl5NmnnjiCVU2vyMdWF/mbwycNMfvjf+XnBTHyO5iHeg7kdcJW9aX33xGRoa67KherKdxkhEnFz311FPKYbLOU7b+KwEBYjnvRoDoHMwCxHI+ngC/LECEciYsQAzwD9mF6gsQXRCWJBUJiASCWwKvv/66zUzV4H5aeTp7CQgQ7SUixyIBkUC1kAAnsXFclGsStUGAn/zkJ+BsWgnVUwICxOr53uWpRQLVXgJcwsGJObQywzFCjm1zIoyE6isBAaIH3z0nCly8eElNJigoLARjUXGxmib+/fffe7BkyVokIBIQCYgEXJWAANFViTlJT/jl5FgQdjAcq1avxbTpszB12kzMX7AIy1esxpq161XkNZ7j9QNhB22mfDvJ2ienaWIrOTkFmzZtwazZ8zB56nRMmTodk6dwOwMLQhZhz559VvuUPqmkg0L5Q+PcufPIyDAhJjYO4eGRCAsLV9uYmDikpKSioKBQWUVxcLucEgmIBKqxBPwOiCdPnlIN2NJlKzB9xmwFDsJj9ux5+G7dBhw5chTaYK+v3xsbX0Jwy9ZtCn4rVq7BgbBwpKSmobCoGMUnTjqNuXn52L1nH6bPmIWIiEN+0UBfuHARUVHRCFm4WMl946YtiDsaj7z8glLPkZWdg/DIQ5g7LwRr1q7DufPnffY6OBaUmZmFjZs2K2jPW7AQ6zdswp69oQgLj1D1DDsYgb17Q7F5yzYsW77y3o+Vhdh/4CBOnTrts7pLwSIBkYD/SMAvgEi7hAcPhmPhoqVYvWYdIg8dhsmchaLiE9aGOL+gECmpx7Fv/wEsWbZCNWpxR44qY8HeFifHHqJjYrFw0RLV8B45moCCwiJrXcsCof01wmbvvv1K40pMSva64Wpa+yeQly4tkSnBkZmVXeFn4Ts6Ep+AkIVL1A+Zm7duee118LsJDT2gZLdp81YkJCaV+0PEKP/snFwciopWGv2iRUsRHhGJ8xcueK3+UpB7JcBeDXqzsF+b6N5SJLdgloBPgcgF1ps2b8GsOfMQERnlElQsuXnYtXuv6r5jg64X9nryZVETiU84hhkzZ6uyCWljA1uVfWqM1LSonVFT82SgM1c+B4E+b/5CJfuqPgvByPcxc9Zcj2tc1MpXr/4Oc+YuQNThGJcg6Owd8fkPhkdi9pz5WLJ0OY7GJ+DateuefA2St5sloG2HEooMXKSv7Zy6uSjJLkgl4BMgsquR4znsbouOiasSVKiZsWts3vwQHD+e5pHXRBCmp2dg0eKl2LJ1O3IsuVWqs7NGmecTjiVi/oKFSEsrcTHjrgfKy8sH47ZtO7AgZLGCV4bJ7PbnSD2ervKn9k65uSvcuHETR4/GY9HiZarrnNq0sQehLJm6eo1yodbOsth1HxF5CEVFxZCJUO56m57Jxx6ItFJDd08SRAIVlYDXgXjz5k3VRcUxnvLG2VxpyCx5+Vi5ao3qSj156lRFn7/MdGwACaZZs+di9ZrvwHEzV+pU2bTUftkFye7AqjTC7NqNjY3DnLnzsXTZShWPJSZ5/BmobX23biOWLVuJ8+er1gV5+fIV7D8QpsY0d+zcDcqmsnKtzH18FvZelIxpz8LefaFW10Nlfjxy0esSsAei1ysgBQa8BLwORDbSh6NjPNaoUcPipIkVK1fj2LFEl8cYqdVwecSBAwfvjU1tgykzy2P1ddZIU/vZvmOXehaTyVzhLmF2Q6emHgcnxBCq1HRcGRN0Vp/KnOe7YLfsli3b1MzOisKdk6YSE5Pw3br1SksLjziEPDd2T1fmWXgPx3vZRcsx7A0bNimzXgHfArjxAWjmkF7laWeUtkR//vOfK+e6PK8D/SHSqe7vf/975UmCaXWgvdT58+erex599FH1Y7Bnz554+OGHldNgGgVPS7vfC0Tj2/RY8bvf/U7ZC6XtVBrw1l2m7C6l/VId6OiYvhjpmopOf2kXl4a82Sbxhzr3d+7cqepHR7+0iWo2m9GyZUv87Gc/U86EjX4adb6yDR4JeB2IbBQr2yC5cl9WtgXUKDiLk4DkOKPJbAZnsbLBpQ82Ts/Pzy9AcnIq9u8Pw6pVa9SY5PIVq3A4OtatGqwrdTemTc8wqy5CLnfg2Nau3XvUhB5qeknJKWrWLZcVbNi4WWmyXCKxddsOsOvSmI8v948lJmPJ0hVqZufmzVsRHR2rure5NCI9w6R+uPD9bOQzzJqrxvH47jzRpesuOcQnJGLmrDnYtn0HODFJQon/P4KEfgbpIYJeJOhv8MEHH7T+oEtKSlLwobzWrVunoEkgMdBTxG9/+1sFxezsbHWOhsB113v//v2Vv0NeoMunX/7yl8rSDK+zPELOCESjNwp64PjRj36k8ub9FotFQdnoWYP7dE/FyVpMT6CzPnR7xTrSSTH9LUoIXgkELRCNDd/xtHSE7g9TYFm8ZBnmzFugxi85lsalEpyKzxmGFVkuYczXm/vUGAkIgppaH7VHgo8TWTgZhNoYJ+Z4s06ulsXuxyNH49WYL2eFEuKM27bvUssfCBlHSzxcLcdb6flO2J0aH58QvC2EC0+muywJPR0IKnq4j4iI0KfUlp4vtHskrdERiEaNTt/AH9H5+fnKqgxNqzHQeS/z1b4QeW7Hjh1OgfjVV1/h73//u80QxLhx41R6XQ6ByDx06NKli9J49TE1zJ/+9Kf6ULZBKIFqAURvNZBSjvN1l8Esm9u3bwdh0+D6I2kg2neNs/tzwYIFKkO6TPrHP/6hXDN17NgRP/zhD5VneV4kEIcPH24tmJofTavxfrqJqlu3rhVgo0ePVt2y1sQGD/UasEYNke6maJrNGOimyV5D5DPoQM/3RkDr59PXZRt8EhAglrF4PpgbcXk298FbgFjSMGpg2C+B4nghu0fT09OV4+B9+/ZZW1I6+tUQIhAJMR02btyoHPTq7tPt27dbAUb/mA899JBOqraHDx9W1x0BkV4s6GfUGDjGKEA0SkT2/ROIxVVsrPzo/pCFizBy1CjXujL9qP6VAue9+ld6DDDAnl+AWNKQaiASZDqEhoaC3ZynTp0CzxNA2kM9He3+4Ac/cApETsB57LHHdFb49ttvrQCjdwrmZeyKpZNpnktOTlb3GDXEFStWqO5OQpmBi/g5XihAtIpXdgD4JxCDSGvjOFOHVzpi0qTJrkExwGWwfPlKtGr9crV4ZgFiSVuqgcgu0cGDB6tIDZDdngxFRUXgMTU1woozPjkxhnBjsNcQNfQGDBiADz/8UHWREmB6Eg5noP7iF78Ar7dq1QrsFuV1PYZpBCK7X7t27arS07E169iuXTsBopK8/NESECB6ATxcyN+qVWtMnTa9WgBixcpVqFOnLtLSM6rF8woQS5oTDcQLFy6oJQ0EYXh4uG5r1JZaIc/Pnj1bGYanNRmCkoH7urtT36T9FC5btkwtc+EYpDbNRsht2LABI0eOVFse8zpniTIwL+ZpDDExMVi/fj3y8vIwadIk/OY3v7Fe5r2cBKQDnycqKkofqmtMIyF4JSBA9AIQ2e1IizrNmjVHSMiioIbE1m3bUbNmLZ+tfaxUF28VvwEBYkkDqYHor82lcWyT76xGjRpKa/TX+kq9vC8BAWIVG0NXGmBqTE899TTYnejKfYGSdueu3XjkkUeU26VAqbM76ilALGm4/BmI7GblbNUPPvgAX3/9tZqxyjWGJpPJ+62ulOi3EhAgehGIbHxTj6ehbr16WL4iuKBIzfDpp5/xqBUid8DLE3kIEEvaN3Y3+muXIpeCsPtzypQp4JINTrIxrmH02xZaKuZVCQgQvQxENsgWSx7YrUiIeKKB9naedIr85JNPKfdc3i7bH8oTIHq1zZLCRAIek4AA0QdAZCNOs2vPP/+Ccn/lD416ZetAU3L/+McTOBR1OCjgXhk5CBA91j5JxiIBr0pAgOgjILLhpR3PGs8/j9D9BwISJjGxcXjiiSer3ZihPTQFiF5ts6QwkYDHJCBA9CEQ2bDGxh1BjRrPK6Pd9g2tPx/T9iu9A1BD9Od6eqNuAkSPtU+SsUjAqxIQIPoYiGyw6bmCUKTG5Y0GvKplsL7UDOmn0D6v48dzsHPnMRwMTy11zT5tsBwLEL3aZklhIgGPSUCA6AdAJBh279mr/LB5w4FvVUCUknpcaYYrV61WwCsoLMbwL7eic7dFqN9oBpq0mIN2ryzEBx+tFyB67N9WMhYJiAQ8IQEBop8AkZCibz0uXYiOifVLmMQdOYq///0fWLp0mbV+I77aijfeXINpU8Mxb+5hFWdMj0C7jgusaaoC4EC4VzRETzRNkqdIwPsSECD6ERDZ+HNs7tnnnlPrFf0JBpwAxKUV6zdstIJu7vyD6PLaUisINRC5rd94JgoKiq1p/elZ3F0XAaL3Gy4pUSTgCQkIEP0MiGysuT6RSzKOxif4BVC4RIQwXBCy0FqfDRtj0artfMycEekQiG06hCDuiMma3t0Q8qf8BIieaJokT5GA9yUgQPRDILKxpxm02rXr+NwmaFZ2Dp599jksXHjfBmt0TAaavzzPKQypIfZ9YxU2bT4iQPT+/7SUKBIQCVRSAv4JxADzh1dKW3FT/TlxpVJeI9xUPm2vEsrTps+wgo1aX8u28zB+fKhDzZAwnDv3MIYM3YIZs/Zb7yslo7J+iLip/i6VaayPi+WLhljJ1kduEwn4mQT8E4jGxqma769avQYNGjRAXn5B5eBSSfmxvDp162LK1GnWctPSLWjeei5GjdrjFIZ6HJHA/OSzTdZ7Kw2nStbfm+UFIhAv7NsHyxv9kdW+HfIHf4RL0TF+1jRJdUQC3peAADEAGtxvvp2o/CnmFxR6BTCEIf03Dh78obU8Ojp+6+1VGPrxlnJhSCjOnnUIfd9Ybr3fm4DydlmBBsTiSZNherEGCkeMwOllK5D/wWCYXnwOBV8Mx93bt73fCkmJIgE/kYAAMQCAyAZ+3LivlT/FzKxsj0KG+Tdu3ARvDXgbhUX3Z4kuWRaJnn1WVgiGWkvs2OX+uKO3IeXN8gIJiBf2hcJUtxYuxcbZNEFXM0zIfvVVWPr0we0LF2yuyYFIoLpIQIAYIEBkAz9z5iy0b9/BBlTubPgJwObNW2DIkKGloLtnTyJ69XUNiG1fWYicHO929bpTHhXNK5CAmNP9NZycPcdh+3b74iVkv9IBlt69RVN0KCE5GewSECAGEBDZQI8d9y5sO7MAACAASURBVDU6v/qq28cUcyy5aNOmLd4fOKgUDFluXn4R2r4S4pKGSIAeOJDsML+KwiYQ0gUKEK9n58BUrzbuXLrstF27WVyMrFYtUPztt07TyAWRQLBKQIAYYEAkIEZ8+RV69OzpNtBwfLBrt27o9/rr4L4zCPV5fTm+/nof/jk+VMWJ3x7A5ElhNnH6tHDQUg3jO++tw8LFEU7zc1ZOoJ0PFCCemDkLue++V25bdt2cCVO9WrgQGlpuWkkgEggmCQgQAxCIBMYbb/TH0KEfuwU21Apbv/xyuV2x80PC8fqbK6yxV9+leK3XEpvYqetCdOxSEjt0DkHrdvOUse9Ag5wr9Q0UILIr9OSChRVqv04tXARzgzq4UVhYofSSSCQQDBIQIAYoEDnjtGev3vj8i+FVguJHQ4aqSTSemqwTF2fCsE83oUefpZgzLwzH0yxVqq8roPJW2kAA4t3rN2Cq8zwux8dXqN36/u5d5A36AJZevXD31q0K3SOJRAKBLgEBYoACkY09uze7vfYaxv9zQqUgQ5jWq1/f7eORjkBkMudh1pwDaPPKfPQfsBLLV0bBZM6tVL0d5e/Lc4EAxKspqWppxd1r1yrcZt29cQPZHdqjcOTICt8jCUUCgSwBAWIAA5EQ4GSYtm3bYfjwEWWO/xmBQZC++977qFmzFjJMZq9CiWUfOnQc02bsw5sDVinNcfCQ9Zg4eQ9WrYlG6P5k0DScjpGHjmP/gRSncd++JNUlSx+MzuKOHQnKjBxNyTmL6zfGYs26GIeR9Vq8NLJUDFkUobTeFauOYNXaeBWXLo/DoiUx1jh3QRRmzY1UcebsCEycEoabN+94vc04u2GjmkHqasFX09Jhqvcizq7f4Oqtkl4kEHASECAGOBAJOi6XaN/hFYweM7ZCcPvwoyGoXbs2LJa8CqU3wtQT+wkJmdi4MQ5Tpu1V3avvDVoLHQd/tB6ffL7Jafxs+GZ8NWpbuXH8hF0oK34zcTemTt/nME6fuR+En30kJFesisLW7cnYvjNVxX2hGQg7mGmNUVE5iI3LVTHuSB5SU4tx585drzcUxRMnIn/YsEqVe37vPrV2kZNtJIgEglkCAsQgACIhVVBYhM6dX8WEb751CDlqZowjRnyJuvXqISvb8Vie2ZyHbyftQWGR89mmnoBiIOcZCF2meQMH4uT8kEq3ZcUTJyGrZXPcKCyodB5yo0jA3yUgQAwSIBIonGjTrFlzTJ4ytRQUR44aDcYnnnwSJnNmqeu8PzMrH126L1ZaTyADytt1DwQgcsH9hbDwSrdH39+5A0u/fshu3w53rlypdD5yo0jAnyUgQAwiIBIE7AZt0OAlG0e+ISGL8Oijj6mYknrcIQy58P61nkuwbEWUw+vehkwglRcIQOS6whv5+VVqi26dO4esVi1R+OVXVcpHbhYJ+KsEBIhBBkSChA59a734IjZs2IiJkybjoYcewuHoGBWdgebz4VswZ26YUxgmJWeD3i6c3V+dz/s7EAlC2i/lUoqqhusWCzKbNMSphYurmpXcLxLwOwn4JxBd9EdXqjGW+5GeYcYjjzyKP/zhj4hPOFYmyLZtj8f7g9beT+NAfu8NXIsJE3ffT1PWDwkH95d6R07uj40zITExy2E5qak5iInNcHjNJv8qlK/yKebs3UI1UadB45l4puYk9Ou/AsfTchyW7e9A5KSY7E4d3db4XDocrZZwnN+12215SkYiAX+QgH8C0UljadPoSRqHjbNRRq+/8QaWLF1WZrqs7Hy06TAf1ACN99rvuwTEKrybZq3nqKUP9uXz2Ft1YFmZmXlq1uuxY1lKNv3fWgnOeHVUL38HYvHUacrnoTsbnLObt8D8Uj1cijvizmwlL5GATyUgQKxC4+2ocfSnc6NHj8HKVasdNuK6njNmHVBLDfSxs60RRtu2J2DV2mjMW3AQ3XouUdplUlIJUC2WQowcvQ1xR0zqPCfpzJi531qHjZuPqOULupzwiFR8/c+d6jqXN9SsNxW9X1+ullnoNNyuXHUYhCXNwnEZRtThNLXlrFidbuXqw2qtYEysCcxr1+5jeOOtlejeeynWfBdtTcfZtgsXRaB77yXKSEBEpONxVZ0vt1u3xaNtxwXWPIzX/B2IuW++CdoxdXc4OW+B6oq9mnLc3VlLfiIBn0hAgFiNgcilFV17LIElt3zHw0Ygcj1e3UbTwbV7dAs14N3V6N2vRBMtLDyBZ1+crCC0c9cxsDu2RZu5mDx1r4IJu12ZlwYKF8MTdDzevz8ZDZrMxNivd2Dz5qPWNLx2+HCamvQz6MN16lpmVp6aETt7zv1xz9bt56uF92EHU1DnpenKofGBsBQQlDzesrUkz5mzDqgu0KPxZmzeclSVn52db1Oerp/e0srOW++scpjGn4HIcUNqchdC97u9gfn++++VV4ysFk0haxTdLl7J0AcSECBWYyASEKPGbHfYyGsQ6K09EKn56WvsVnym1iTr8RM1JoJQ0tdXr41G3YbT1XFZQGR6V7pMF4SEK+PhvI/AbNhsllo/Sa3zuRcnI9vgi/Gf3+xSEGTa5i/PRXz8/aUnNFi+bkOstb663nqbnpGrfgBQS9TnjFt/BqLyXFHzWdw6fcpjzUvh2HEKujQPJ0EkEMgSECBWYyAOGvwd9u5LctjIGxt87tsDkZNMdBp2WxKC+pj7xq5Mdm/yXFZWvpqY40xD5P2uAJGzXp+uOQlRUWnKCs0XI7aoOhCIjZrPstaH+dL8GjVV7j/1/EQ0aDzDJtLijK6/ccvu1bffXY0hH29weJ1p/RmINLmW1bqlR9sorlHMH/yh8o4hUPSoqCVzD0tAgFhNgchZlM1bz0FBYbHTht4IBleByBmh+v7tOxLwQt2pSnv7dvIe1ZWpr624Nzaoj10BIu/hRBdqf01bzcG+0BK4E4j1Gs2wls90s+aEqW5c7teqPw3G+umy7be5uUXgj4Z33l9Tppz8GYj5Qz9G0dfjPdyMQC3pKJ4yDZkN6+PCwcobAPB4RaUAkUAZEhAgVlMgsvuvLK3HHg6uApFdsdSuOE751tur8OnnmxSgOJ5Xq8E0pGdYQGMA1DQJQV3eK6+GYLaT9ZDDPiuxaarTckv4cXywVdu5VuPmPEeNlGUxDbtO6WWDE2l4PHDwdxg6bAM43snj/ILSPwoIw559l4FOkdPSLErjpdbLZ+I9xuivQLx7+zbMDevhSlJyGU2Aey+dXrVaLck4vXKVezOW3EQCXpCAANGucTM2dIG+X9Ys09FjdygXTBV9xq9Gb8PcBQcVCDgR5pPPSgDH+9kV+nKH+VZIEEajx2wHnQVTC/3gw3Xg8g6mteQW4Z331uDFBtPQoVMIZs4+gF73JuTw+tJlh1D7pel4tdsia366jnrSTbNWc2w0PAKVE310OgKxfuOZ+PiTjaoOTVvOxphxO6z2WVOP5yitr2WbuejcbRF69S29NIWzUvlM9pHPqsvRW38F4tmNm5T7Ji+0IzZFUEPk4v2Czz7HrfPnba7JgUjAnyUgQKymQORShPh4z7h+sh9D1ODwxJZjoE1azrbR8ghE+zFET5St8/RHIFI75Njh6cVLfdL+0DpOdscOyGzcEJeiY3xSBylUJOCqBASI1RSIjZvPdtj9pxv5qmy9BcQlyyKVNrlx0xEbrU2ACJyYPgOZTRvjzuXLrrYJbkt/9+ZNnJw3H+YGdVE4fARunTzhtrwlI5GAJyQgQKyGQDSZclVXYlWgV9a97L7k+GBZadxxbfqMUOVQ2D4vLpPgZB3785469jcNMXblCmTUehZbxozGsuWrfB7XzJiBhA5tkVb3Bex8ZwCWLVzs8zotDwnB8n9OwJqvvsTGj4di+7vvYc9br+NA9+6I7NIZ0Z3a4Uj71kho0RxJzRohuclLSGlUV8W0ei/Aaaxf05qO6ZOaNlQxoXkzlRfzjOvwMqI7d8Dhzq8gsuurCO/aBWHdXsP+Xj0Q+kYf7H2zH3YPeFPF7e+9i8SDYZ5o+yVPBxIQIFZDINKSC5cSeAoQ1S1ffwPihX2hoENgfwpcxH/2u/XIbPwSMps1xskFIbh1+rRHq0gN9Wp6Bs7v3Kk05vwPhyC7c2e1PMRU81mo+GINZDZtpLp3Lf1eR94Hg5U2W/TPCTgxbbqq5+kVq3B27Tqc3bJVRRo5cBZpN1anO7tpM06vXqPiyZBFOLlgocqTpvSKxv8ThWPGqbLyP/0M+R8NQd7AQch9+21Y3ugPS+/eKma/2hlXU8USkEc/FEPmAsRqCMTduxMx/MuSNXvVDV6eeF5/A6Lh/9vvdu9ev67AaOnVC6Y6L4BbmpW7FHkIdC/larh94SKuHk/DhX37cGrhIhSOHoPct95C1sutYCLsGr+E7K5dkT9kKIonT1FguxgRiWsms5rw4w4PIK7WWdL7rwQEiNUQiBs2xiqza56AQ3XMU4BYuQbu1pnTOL1ipQKYqd6LSmPLbNwAlh7dkTdoEApHfKk0KWpTSqP68kvkfzAYltffQFb7djDVq12i5dV8Vo2XWvr2ReFXI3F62XJciomVGa6Vey3V+i4BYjUE4vKVUZjlZK1fdQRaVZ9ZgFj1NpSzYq9nZauuSEKSWqOC4MhRKLwXeczJQqcWLQG7Iy8disL1zCzcuXq16hWQHEQCAPwTiG7wZ1elRi5Iyne2DnHBwnAsWhLpfAwxSJ6/0t+Ai89f3YF499ZtnJwzt2RMbOQonJwfglsnT0oDW0UJcIbwleSUKuYit7siAf8EYhBrbZVupCshE2dAnDMvDM5sd3qzfsFSVnUH4p0rV1XX5cm583B+715wQgonqgQTFKmletOjx8mQhaAXEU78keA9CQgQKwGaQGnInQGRroxWr/HesoRAkVdl6ylALAGiUZvhON/p5SusLRm7RK9mmHD7wgXrOe7QMDhB4wiet86fw9V0E3ivfeAMVXXt5i37S+r4RlEhbuTllbp269zZUvdxBuz1HAtuFBaWSq9PZLVo5rK2xvJpoMA+3L19R8ni1tkz9pesx3TAfCk2ToBolYh3dgSI1RCIU6fvw/oy3B1VFgzV9T4BogMg9uuLs1u2qFaMIOSsT8vrr6tZn+e271DnOeOUM0CzO3dCVsvmuHIsEVeSk5HVrIla7pDVqiUymzRCdocO0PC4e+MG8j/5RE2iye7aBVnNm+BKwjGV3+Wj8SC4qM0p7arOC2pfN6W0s0rNNee115D/4UcldbhxE3SgzIk8WW1aq/QEpDFktXkZplrPKmfIrC/BT83tek6ONZmlZ09lhODUwsXIHTBATQrKVhN/XgSXe2io0zgBl36oiUEtmoJ1chZ0Oc6uy3n3S0CAKEB0PpYYxLJxJ7wFiCVAPL1kKS6EhirvGpZ+/UB4MRBclw5Hq31qTIQWxx3Pbd0GpmOgpsjI5RCED2eL8pggIbzyBg5U6TjzlEs17ly9po7PbtiIzEYNcPvipRJQ1XoW2rD41bQMldeV+ASV1lS/Ni7FHVH7d2+VaJb0BMIuXoa7166B8LsUdVgdG//Ya4gKgLNnqyQ3i4rUEo8bBQVqzSHrf2779pJrJ06o5+V6RAauLzy7bv29a8Uw16/tUJNlAgGiEpNX/wgQg7jRd9ZlKhqirbeKqsJRgFgCxLz3B4LupmiqjWBkuFlcDPNLda0L1AkGrg2k1sjIpRMce7x98aJKr4BIh8aGSTmq67D28wqOpro1cWH/AWsjyXWEzP9ieIQVIEZzddQyT69Zq9JzWYalTx8b7x+sS/E331rrR02RE4Tsgz0Qz+/ejeyOr6hkp5etUFomD9TzNW1scztnzNJfJGVBzZL2ZfWCfRopINQdBQGiI6l49pwAUYAoGmIVvwEBom2X6aXYWGWNhrDionlCh9ZejFEvwmc3J9cWEgxXj6eXaIg1n8WdS/dtsF5NTVMguX3pktryHmNgd+v5PXvvA/HK/WUY7AbV2tmdy1dwctZsVVbhuK9VFly8Xzxxsk3dCCL7YA9EjgOyW5ZWZNj9eSG0BNIsi2Uaw+mly5H77rtKFgqIy5bblMelI46CANGRVDx7ToBYxcawqtqFJ+8XDdG9mqCzdyVAtAUimyy6fyLguEaQ3YLluYEqHDUGhaPHWoF4NSXV2vLRbFp2h/bqmAvyT8yYab3GLlg1nmexlAtEfRMNAphqP68Oszt1BN1klRfUGGdikk2y4omTlIsraqHa4g2BSCMDHB/VgcBnWsrCVKcmrqaUBq5Oa9wKEI3S8M6+AFGAKBpiFb8BAWJpIOa+976aGMNmjF2S2V1eBWdOnt+1ywogdnOyO/P8nj2gzU7aANVdppYePVQ6pXG1aIYL+0sMXF88GAHzS/XUBJbzO3aqSTl6DFADhMtAdNAaIu2a0r6rsvM6dZq6j2mYH2HHGbHs5qV3Do4l2gfaOaW9UV0PXufsWNOLz6nn0+k1EAlBPhdByB8HerYpl1NwkhC7Sc/v3oPTi5foW0tt9fOUuiAnPCYBAWIVG0NnWoM/nBcNUTREj7Uchow5QYYguHX2rPUsNbwLYQfVMWdtskuRxqyLJ3yDCwdK4EagEBiFo0YrCzVMrIF4PTsbJ6ZPR9HX/4RRW2QaLo8onjYdhWO/xsWIQ9YyWT7rwfrocG7bdgWuuzduKjNxnKxzcvZs27qmmxQsqaFy1ufda/e1O53P9dxcVX9a0dGB5bHL1TjbVAG8TWtcSUxE0bivcWLmTJuyeC/9Q/Ian40Ti5wF/TzOrst590tAgFgNgTh56l5s3BgnmqGb3n111xDd2SxpIBrHEN2Zv7vyIuSVl4qhH9tkqYFoc1IOAkYCAkQ3NYr+oBHa18GZhjhh4m7s2HlMgOimdy9AdF97FyhA5JpKLqGwB7cA0X3fgi9yEiC6qVG0h5E/HAsQpcvUF41KdSiTs0u/v3O3OjxqtXpGAaIAUTTFKn4DoiFWqzZTHjaIJSBArGJj6A+aoLM6iIYoGmIQt13yaCIBt0tAgChAFA2xit+AaIhub5ckQ5GATyTgn0B00R9dKQ1J7leQq7SGKPJz6UeCANEnbZcUKhJwuwT8E4hV/MVeCpDVNL9KA7Gayquy340A0e3tkmQoEvCJBASIQdz4CxBlDNEnrYoUKhIIUAkIEAWILnUPVlaLCub7REO0bf1u3ryDc+eu4tSpyzhz5gpu3LhvOcY2pRyJBPxLAgLEagjEseN3Yu++JAGhm959dQXipUvXEXYwE9NnheP9wRvQ4dVFePGl6XiixkQ8U2synq87Fc/VnqKOud+y7Xz0e3M1Ro7ejWUrjuBwdA4uXixtN9S/mkipTXWSgADRTY2iP2pAzrpMR47ejvCIVAGim959dQNifHw+Bn24Ec++OBm1X5qBD4duwpJlcYg8lI28vHOlNMJbt+7g5MlLOJZYgM1bk/HNpP3o238VXqg3FU8+PxFdeizFzNkROJ5WDHtv9dWpMZZn9b0EBIhuahQFiN4Zr/NHOVcXIOZYzmLg4A14+oVJ+OjjzTh0KBt3qmCt5fbtu0g4VqBg2KX7UqVJtmq3ADNmR8BkOuX71lFqUO0kIEAUIIqmWMVvwN+AWFx8Ebv3prmtMeOY4NQZB1GjzhR8PnwHsrPPuC1vY0bUIjdsSlTaZ60G09C+00JMnyVwNMpI9j0rAQFiFRtDf9RYdJ2ky9Q7Wqu/ATEmNleN4RUWXqhy63Hi5CV0fm0J2nQIQVr6iSrnV9EM2M26a08a+r+91tqtumFTUqnu2IrmJ+lEAhWRgABRgCgaYhW/AX8DIsfhevZdgS++3FGRNsBpmvT0E2j+8jwM+2yrTye/EOyLlsSgc7clatLOp19sU+OVVemudfrQcqFaS0CAWMXGUGtj/rgVDbF6aohs0aJjLHjqhUkoKrpYqQYuLe0EatafhklTS5z5VioTD9yUmlqM0eP2qLo1bTUXCxYexvkLMlPVA6KullkKEKsJEAuLToAxKzsfw7/cij17k5CTUyDaoRvev79piLol69ZzmdKs9HFFt5bcc2jcYg6mzQjH3bvfV/Q2r6bjko9NW5LwxoA1amnH4CEbsWdfOq5cvenVekhhwSUBAaIbGkR/1A5ZJ6OGWFBYDMZW7eap2XxcK/bZF5sFiG54//4KRK71e/u9dS61WOfPX0WLNvMwauwel+7zZWJ2qU6eGoY6jWYoOHJJyKYtycoogC/rJWUHngQEiG5oEAMBiLqOW7YeVUCs/dJ0ZJhyBYhueP/+CsT4hHy1GL6izRLHHt8duE7N7rx27VZFb/ObdDdv3saBMDM+H7EDdRvNVJNxXuu1HBMnH1BjjldFe/Sbd+WvFREguqFB1LDxt61RQzTWrV//FVi+Ikpg6KZ3769ALCy6gFoNple47flu/THUazwTuXnnKnyPvyZkVy/HQVesOqomF7XtuBBP15yEDp0XqUlCIYujlZWdnJyz4IxWCSIBSkCA6KZG0Qgcf9l3BsSU1ByBoRvfu78CkbZEaTKtIuHipeuo22gGQvebKpI8INPwGcMjstTaxgHvfqfgr83Mte+8CIOHbMKU6QexfmMioqMtyM45iytX3D8mSYME7OZNTCxEbFyuiimpReC4qDEQ1CxfbMEapeLZff8Eovjjqxqw7snPGRDLBbbI3yX5+xsQWZ/jaemIT0hR3ePcLy9+NnwLunafp9KlpqWVm76s/ALp/qjDCVi1OhQTJq7HBx8tRLceM9Cq7WQ0aDIJL9SZqORXo85UNGk5Fx27LEavfitBmBKe7Jod8/Ve/PPbUEycEnY/Tj6g9r+ddECZqeOsWFr26dlvBZq1nqc01edenIg6L01C/cYlsXaDSaqLl0MZzV+eq2BNs3YENg0i3LotBtI9i8KS3P0TiG789V5u4x/EZVUaiEEsE098D/4IxOSUVDCyQY07kqj29Tn77e49h1VjvGNnbJnp7O+rDsdJyccRF5uCgweTsGdPErZtT8S6DcewdHkcZs2NBKE3dvxejByzu1TkxCTGiVMOqNm+W7cl4kBYEuLjUxzKOTExFdHRKSrNoUMpOHas5B2yDrduBd6YrjcA5u4yBIhB3PgLEKvvOkTdUHAtYVFR2RZr2GBz+YIEkUB1l4AAUYDoUvegJzSsQM/T3zREY6PGxevHjxcbT9nsczyrfpNZoLk3CSKB6i4BAaIAUYBYxW/An4FIc2cRkdlO27l9oRlqaYa4XXIqIrlQjSQgQKxiY+jP2o10mUqXKSeArNuQ6LRJ6/PGKsyaE+n0ekUuREXlqDEyjqmtWhPvMW8YFalLMKThkpFjxwqwd186sjzkWSQY5OSJZxAgChBFQ6ziN+DPGuLwr3Zi2sxwh22HyXxaTbrJyz/v8HpFT77z/joM+3Qrtu1IxeRpYcpaTERkVkVv9/t0NHBA4+LeCDRY3qnrEjWblb4nn6s9BctXHvFG0VKGrEP0jgbhKy1SNETvvF9/BiK1vyHDNjts7CZM3I93Brpm2s1RRgQiF/XrMH1muE2Z7I6lkfFTpy/rJNbt6dOXUVBwoZTNVGpJNBBgvzaPN3LcMzf3HC44MerNtXt0ZmzvDYN+HWmn9epV2xmbtGDDBfo3nSzQX7w0Vi2bsFa6AjusN8uyrwNv5frQomLnRtdZFx2iDufgpWaz9aFsPSwB0RCrqB34CnYVKVeAKEDcvDUZXHRuH2iarfZLM5THevtrrh7bA3Hu/EP4YkSJ6ylCiPZUuX6PPhW5REGPV3Kf7qXob1F32743aL3SxugcmPV+ptZkrF2XYK0SDXi/1HQWaJKNhgTG/XOfusZF7FxismL1UXAiEdfzMc316yXr91JSitTkoR59VqBhs9nQpul27j6OVu0W4M131qJ1+wUKztbCALVIn1oardzQaTG1YM7cpVx1CFkUjT79V4GL65u0nAMCtGXb+aqcV7osxtmzV1VSyoI/Tl7tvhSvvFpiMUfLQudlv6Vt2edenGJ/Wo49JAGfAbHIkg2HMScThRnHy4wFaakoSEpwKeYfPYL8w9Eux7xDUcgN21+xuH8v8iIjUZh+HMXFxT7vimzfvgN69Ojp83pUBN6BnMafNUQ2uGzs337vO7Bx5TEn0tCCzf4DZrc0KwRilx5L8eHQTQpG304+YM2XcJo9NxK796ap2Lf/KixcEqP8K9ZvPBNMe+LEJWt6LmCntRhj6P36SgUgguatd74zXsL8kMNg+QwEYuSh+xOIxn8TqhbO89rcBVEKQryuIURwd+25zFq3hYuj8UK9aTb588BeQ8zMPK1ca9HaDDVcAprg5XnWgRqvDoejLUrW1HgpCz6blgWNkH8zab9OWmqbX3BeAdb4TKUSyQm3SkCAWA4kBYje0bIEiG79v7ZmJkDcq2QhQLR+ErJThgTcBsRL0THIatO6QtHUuhlMHdrD1KFD6djpFZh6di079uoGc/8+rsUB/WF+723X48D3YB4ytGLx449hHvg+TN27wvRyY5heaQfToPeQNXMa8sLDUFxY6DVtbfuOnejZsyeaNGmCrdu2e63cQAZbZevuzxoi//c5pkYti91+dJHEbsgly+PKaBZcu2TsMqXG9OyLk5GSWrL2kebH2DU6e94hazxyNE8VQE2KGiHrw65LBh5TIzMGmkmjKyu6eBoxcpfxEjZsSlRdozxJ7SzWsJ6SXi5oQYaBY3nMg12anFlL+6A0qdal+1JrvVhHgtM+2GuIvM5npqbLPL8cVVInrSFeNtg/zTCdUvVieZQF62OUBW2rOgqUY4s285UDZEfX5ZxnJOA2IFa0evzFWtmGJ9DuK8zORF7EQWQvWYTMLz6DqU93mPv3RdbECciLDK+wHApN6cjdsxPZC+bB/PmnMPXtAVObZihITXaYR25ePp577jlkZmYiNfU4nnuuBvLyC1QMNBkGQn39HYj83yQQEpMKVZfimTNXKvrvWqF0RiDyBh7TmwQDxwF37UkrMx+Oy9GZMQOBOH5CqDU9uxpf7rBAGQ7YsCkJBMMD2QAAIABJREFU9FrBSTU6jB67B4wMZQFRp+ePAzo/PhiepWZvsquXZZQVliyLU93BxjSc7NKj7wp067UM0TEWdUkD0Zx52pp06/YUNTbJE5TFGsN4qDWR3Q5hyHFGmV1qJxgvHAoQvTyphpDMWbMK5rf7w9SuJbKmT0FhdpYVbIWZGbDs2IKsSd/APKA/TK2bIKNFY2RQ6xz0HsxjRyFzziyYPh6CzOGfWu8zgmPc1+MxbtzX1s9n+IgRmDxlqorGdLLvnu7gQACi9WPwwI49EKl9DvpwgyqJ42UcOws7mInMrDNqUgonl3DW6aatycjKOgOO9VELZCAQOXll5eqjSE8/gREjdyoNkD+kCTM6L+bEFLp2Yhpa2WFeDGUBkWv6UlOLFbxq1psGam6XL99Qk3pYPpegHI3PR1JykcrL+IdjrnwGrrc8f/6a9VKbV0LU/RqoGojdey9XeR06lK3u274zVd1DWdAd145dqWqtJo/tPVmcO3dVTfIhbKmZ6lhcxqxUa4Vkp8oSECB6GYhGCBWaTciePxemnt1g6t4Fpk7tSrpZX+8D09CPYB4/VsEvc94clIpzZsLUojEKUpJsoJiWbkLDho1w8+Z9tzVXr15F/QYNVMyxiFNg4ztwx351ByI1PMJGh+ITF226ZGkp56vRu/DpF9vUbNIrV2+qiSdc9jH0ky2YNiMcZ86WaK0EIifhLFoai08+34aZsyNASOjA/QWLDiufhpOmhqnlEvoaDQMYJ7QQYAQxw5ZtKap8AtY4SeXs2SuYOSdS1YNdu/Hx+To765YaKevx8adbYbHc9xVJower18Zb02kgEqzMizNt6bDYGNhdTA8Zwz7bqtaHXrp8w3hZuYXic9hH43PZ3CAHbpWAANGHQDQ2xuZPhsI0ckRp8DmC4b1zpg8/QNY/v7YBYs9evbFhw8ZSH8mCBQvAOGTIUJv0xjrIfuU0xuoOxFIfWxVOOBpDrEJ2HruVWh6XkRi7bzUQjWOIHquAZOwRCfgnEKuhPz7Lti3IeK2LAqJ57uwKgTGjRzfkLF9mBdzBgxGoU6cu7t69P8aiv5obN26A8W9/+zuSklOs9ziEYDWUv40cXHx+AaL+yqq+JRC9ZRWmsrWl5R/6KKRzX2OwAvHy/d4Z43XZ938J+CcQ/URrs2kkPVynooICmFo1gXn6lArB0DxuDExdOqK4qMgKt7Zt2yEp+f6CYUef344dO/FG/zet93jzGYO1LAGioy+tcufY3cr1kv4c2CVaUFDa3N3NmyVrEfWYoj8/g9TNsQR8DkTL+nUwvT+w2jXQRXm5yF68GPlH463PnvX1GGS8/06FgGjq2M5mpurKVavx7rvvOn7LdmdfbtMGMbFx1nKDFVTeei4Bot0HJocigQCVgM+BmL1ooVqL6OnGK2vWDGSOGO4QApbNG2Dq7T2LLrl7dsHUoQ1M9WsjZ9Vqa53yk47B1KIRMsvpMjV9+QXMg+//iMgwmfHUU08r6zgV+Q5jYmLQpGlTFBWfsJbtafkHc/4CxIp8dZJGJOD/EvArIBblWJAXdRhF+fmwrFuDnLWrUVSQb220ea0wOxv5MYeRvXiRWuOnG9qinFzQqow+5jbvYDiYZ0FqCsyD3of5zbfUuUKT2ZquMDMTWV+Pg6l1K3WtICUVeeERKMww3U+TnVOSV35eybagALn7dyN7yWIUHLuv4bFM1jdn4zrkrFhus5zCWK/8xAQUpqepHwJGIDKNuX+/kqUVzibTcMlFs4bIT7hfbt9+/TBv3jyXvrbeffpg3vwF1mc01k/2XZtcI0B06dOTxCIBv5WAXwExL/IQTE3rw9SvJ8xvvaGWFZhee9U6TmZq2wqZg9+FuVM7mN97Bxm1ayDrm5JZloSfqXkjmwbeVPNZZYM0Z1mIgoipSQO1xMGya5c1Xe6BPTC1bQFTvVrqGrsxTX17InPUV9Y0WbNnqLqwm5N5Zn70Lkx9esD8el9VB8KbEFHjgD27wvwBrdsMVvUpOO548TzT01KPPRAt362G6fW+TrtNTR8NhvnzT6x1i4iMQs2aNXHnzh2XPrL8/AI8/te/gov4BYCuAdBeXgJElz49SSwS8FsJ+BUQaXybwGEXJhsdGvEm9CzbNqtjAsT8eh8UF5UYzrbs2oaMF59FQfwRpbk5AyLzyvxqhDK/Zt+Y8di+29aydaMCoDbQbe7cHllTJpbUgUD86n7Xa9a0STA1rIOivDxkTZ8K89APrYAxfzIE5s+GWY/ty3YERGqYpg4vOwSiefZMZHRoDa5f1Hm179ABCQn3vQG48qVNmjQJo0aPtual85Sta4AUILry1UlakYD/SsAvgVhkuW+5xdSqKbLnlXTtESA0X2ZssE2tm6vuybI0RKZ3BYjFxSdgatkEBG7+0Thk1KmBQnOGKldpnaF7rHUoyslSEM87dAimnq/B9Fo3q91TU9euMHXuaE1rrDf3HQGR59WaxFFfloJixoA3lGUbnc/ateswYMCASn9d1Crr1q2HrOwcp3XUZcnWOSQFiJX+BOVGkYBfScD/gdi6mS0QF4bYNN6ml1siZ+WKEg2xWUPba/e6TNmYuwTEEyeRNWE8zO8OQNa0KTB/9IE1XwXEA6HWY7qw4rm8qEMwdXsVmWNGIXfPDmssy2apMyDm7t6JjC4drUA0T50E0wfvw9SmJYosFlU2IcaJNIWFhVX6oNasWYM333zL+jwCPufgcyYbAWKVPkG5WSTgNxIIOCBmDv/M2nhzck1GrWdQkHAUBceTVPeqtguaH31IgYq+DBUQR42E+YP7YDM2btlLlsDU9mVrvryWfyxBaYamLq8gLyLMeo3wM2qpOcuXwPRSbTWZhl2p5iH3u0yNZTjadwZEdglzTaJp2EcwdWgLU/tWyJo6CQXJidZ6fP7FcMyYObPKHxJtRNauXQfhEZHWvB3VVc45B6UAscqfoWQgEvALCQQcEDkBxjzkI2SNH6cmymRNuG+6zNS7B0ydOyjPEuY3eqvrufv3qoY+Z81Ktcwhc/jnqhvU2MBTu+NYpfnjIcjdt9sKBlO/XjB3am895j0EorlzO6VxUus0Na6HnLWrVJpCU4Ya/zNzLeGkb2H+dBhyQ0vKN5an950CkeOa8+cgc8zokpm0dssjEpOS0aRJU5cn0jj74uLjE/BKx042z6nrKFvnINSyESA6+7LkvEggsCTgcyAWZmQgLzpGNcZFuSXLGozWV7jUoiD93vgdxxAXLQRnhnLJA8f3dKPELWeB5qxeqWKRJUfly2UaOo1lx1ZkL1yIwox06zl9Le9QOLIXhCA/4WjJNY4jdutsow0yreoyDd0Dy7ZNyFm2FAXpx23yUktGtmxU5eTu3qGWkOgy7Ld8buMSEPvrzo7btW+PvXv3ufVL69evH7Zt22HzLM7Kl/O2kBQguvVTlMxEAj6TgM+B6ErjSo2KQHTlnsqmzZ46EaYWTVCUa+sdQgHxXjdsZfOuyn07d+1G02bN3P7BpKam4oUXaspi/UqY6BMguv1zlAxFAj6RQEAB0fz227B8V7LmrypQKe/erLkz1dhd/tEjpeBratNaTaApL4/KXM8LD0NxYWGpMnVetCzTrHlzZGY69rJd1S9o4MBBslhfgFjVz0juFwkErAQCCogaDJ7esrtVr3X0dFk6//zYGJjatSqlkerr3M6ZMxeffPqpxz62y5cvK28ZhffWeRrLln3bblKjPERD9NgnKRmLBLwqAQFiJTQCY2Polv2iYphf72VjrNs+35TU43jiySdx7tx9B6We+FImT56MUaNksb69/Ms6FiB64kuUPEUC3peAfwLRRX90pRqrALs/a/Z0ZI4zQMhB/Tt27IQ1a9Z6/AvRPhMTjiU57botJW/7HxUO6l/uPcY8Aux+AaLHP0spQCTgFQn4JxCNjWOQ7xekJsP0cjNltNwZNPbs3YdGjRqDawa9EVavXo0uXbtVHohB/s7s35MA0RtfpZQhEvC8BASIPm68aQicy0HsG1l9zIk0LVq0hMlk8vzXcK8EgpdlckarrodsZQzRax+gFCQS8JEEBIg+BGLunp0wD3q3TOgsCFmIoUOHev3z4ExWQlFA6ByEWjaiIXr985QCRQIekYAA0UdALMrPg+mVtjZ+DXUDq7fpGSY8+eRTOH36jEdefnmZchnGSoMDY10v2dpCUoBY3pck10UCgSEBAaKPgJgdMh9ZX48pUwMb9MFgrFixwmdf0okTJ/D44+IzsbwfAAJEn32iUrBIwK0SECD6AohFRTB1fxXaELmjBjc27ghatWrttYk0zr6qiRMnYfSYsWWC21H9q9M5AaKzr0fOiwQCSwICRB8AkcbAM8eWrR1yVml0dLTPvyY29nXq1BWfiWV8JwJEn3+mUgGRgFskIEAso6HzlJZj6vUa8mOjnWpd69ZvQKdOnd3ygt2RyYIFC/DRkKFO6+spOQVKvgJEd3xlkodIwPcSECB6GYj5R2OVFw1njT3NplE7LCws8v3Xca8GN2/eRM2aNRETa+tdxNkzVLfzAkS/+VSlIiKBKklAgOhlIGZ9OwHZ82Y71bYmTPgGY8eOq9JL9cTNkZGH8Fr3Hk7rXd0gaHxeAaInvjjJUyTgfQkIEL0MRFO3TshPOuYQLEeOxuPpp5/GpcuXvf8lVKDErl27YV/ofod1NwKiuu0LECvw8UgSkUAASECA6E0gcnZpm6ZOgdK0WXPs3r3Hbz8b5TOxpvhMtAe+ANFvP1mpmEjAJQkIEL0IxIK0ZHBCjX2DyuONmzajXfv2Lr08XyTu3acP5s6b7/AZHD1XdTgnQPTFlyhligTcLwEBoheBmH84GuZ3B5SCCSfSNGnSFHl5+e5/w27OkYv1n3zySZjMmaWeozrAz9EzChDd/JFJdiIBH0lAgOhFIOZFRsI86P1SIJk0eQpGjhzpo0/A9WIXhIRg2CeflnoOR7CoDucEiK5/Q3KHSMAfJeCfQAwwf3ilGn0n9c8/Fg9z/742IElKTsEvf/lL5ZR38pQpYJw9ezbmzJmjYkhICBYvXmKNa9euhY5bt27Fnj17rDE8PAKMBw9G4PDhw4iLi7PGlJQUGGNWVhYKCgut8cyZszh7tiRevHgRV69es0Z7t1N37txRGm1ySqrNs1jl4OT5rdfL+xESYPcLEP2xaZM6iQRcl4B/ArG8BjNAryuD3m1b2EAkMSkZq9d8p4xo05D2ipWrELJwkTVOnTYNE775piRO+AYjR46yxs8//wKffPKpNQ4a9AF0fP2NN9C4SRM1a7VHj554pWMnm9imbTu0bNXaGhs2bISXGjZU8YWaNVW3KLtGGR955BE8/HBJ/Otf/4annnoaDz/yKD4e9onNs1QYeAH6/pw9nwDR9YZH7hAJ+KMEBIhebpxNndqhIM2JZuXmuqxZ+52CpbOGXM7beq2orDwEiP7YtEmdRAKuS0CA6GYIldeoZo0fh+z5zhfml3e/K9cFiO4BXnky92cg7t+/H8uWLcP169eRn5+PYcOGYdw4/zP84HrT5f93zJ07FydPnvT/ikoNrRIQIHoZiPlxMTC92sErXY0CxOAC4kMPPYT/+I//wH/+53/iV7/6FVq3bo2dO3da/5lXrlyJBx54wHr8zTff4OGHH0b37t1x/vx5/PGPf8Sbb76JRYsWWdN4cufSpUuezN7v8/6Xf/kXxMTE+H09pYL3JSBA9DIQqW2YB76DnE0bPA5FAWJwAZENbFRUlPrvJWwIwP/7v//D5MmT1bmcnByr/8y7d+/iv//7v7F37151bd++fQqinBDljZCRkYEf/ehH3ijKb8sQIPrtq3FaMQGiD4BYmJ4GU5cOKLJYPApFAWLwAdFe46CLMGqM165dA2cOr1q1Sv2zT5kyBWyQv/zyS7Drjprh//7v/6r93NxclSY0NBTvvfee6kK9fM9c4LFjx0B4HjhwwNq1euvWLcyYMQMDBgxQ2iVhy7Bjxw4UFBSoc++88w7WrVunztMYPMtl+SxbQ1ldBJCenq7KOHLkCMaOHatOc2bzhAkT8O6772LXrl06KdgdPXPmTHV+27Zt6jzrGBkZqbSvDz74QOVBDVgH3sNuYubFPI3dlnS4zR8T3377LQYOHIhDhw7p20AZMP2gQYPU8+sL2dnZ+OSTT/DRRx8hLS1Nn7bZ8llYlzFjxoBrdRn4/JQxNXXeazabrfdwf9SoUXj//fetcuNFzg5nXuzq5ruhfIzPxjR8x3wXEydOxNWrV1We/KFDrzS8x5dOxa0PGKA7AkQfAJFaomX9dzD17o7CjHSPQVGAGPxAZLvz+9//HuHh4Wo5DmcEM4wePVo1yGyI2cj36NFDaYjcz8zMVEt7mjRpomA1YsQINRuZ4OOSH84sbtCggYLZjRs3ULt2bQU4gq1du3YKDiyjU6dOqFGjBr766ivVOFMjXbJkCXgPyyUQWN7mzZtVnfSf1atXqzKef/55TJs2TY1v0psKIb59+3Y888wzmDp1qkpOaHXu3BmbNm2ygoOAfPzxx1X5BB/rwbwIQi4R6tChg7pn48aNasyUPwQ0jP7617+iUaNGmDRpEj7//HP827/9m7VbkzJ66623VDksj4FLlyjTNWvWKI2cXdB5eXn6UdSWkPyf//kfBSLWjWO1DHz+hg0bKqC//fbbSmPX15YuXareF8thnSg3hi+++EK9i969e6vyjM/G6/xh06dPHwXaDz/8EG3btlX3cVy4adOm4DPPnz9fnZM/rktAgOgjICoobt4Ic8e2yAvzjMFsAWL1ACLhwLFErk/VQKTGxQaZa08ZjNc4weZnP/uZFRK8/uCDDyIiIkJB8Mc//jHOnDmj7tuwYQPq1Kmj9vmH53/3u9+pY4KK2pQOBCPPMbBclu8osC4sQ2tSCxcuVOOcOu3u3bvxpz/9SR2+9NJLVk1VX6fWyfFU3f1LjZXHlAE1Rz6LvsZ7unbtqgDN/b/97W8gjHRo06aNgj2PCVUNJn29W7duNue0Fqivc0uIUyYXLlwwnlbPz7XCDKwjf7iw7jqwjtTqKUPCmIGaNd+n1sKNz8Zu6F/84hfqB4fOg3lSAyYk+eNBQtUk4PiLrVqeZd7NX3DlzdqrTtfzog6pSTbmwQORH+fcaXBlZCJADH4gUisi3NgNaYReWUBMSkpSjTUbUx0JKHZ5ssH+y1/+Yv0fZpfdT3/6U2s6pv9//+//qeuEH7sDdeC9LVu2VIflAdFYBrv52O2r6/LrX//aWgYhTe2rVatWIBAYjOXosqkVUuujltnezibw9OnTQfAxEIjsjtSBGteQIUPU4fr16/GTn/wEr776KoqLi9U5ep9hfXTdWE92VxoDNeLmzZurNDSooQ1Z8AeBsYv7ueees8qLmiF/vLRo0UKBWP+QIBC7dOlizF5pvHw2an98T7ou3HKcNj4+HkePHlX1pGbP9yuhchIQIPpQQ7RCrvgELNs2wzzoPZh6vIrM0aOQu3cXigoKqvTjQYAY/EDkhBp2MTJUFIgcE2NjbT82xTwIG61l8pgNcceOHVX+9n+qAkRjGRyfM2qa9uWcO3cO7NZ97LHH1CXWkd2DxsCuSWp+nEFLKBjD8OHD0bNnT3WqLCAywenTp8Huzaeeekqlf/HFF0GNtSIhLCwMf/jDH1QXKdM7A2JhYaECW2JiosqWEDUC0V7e+tlYD2q/zgLHE/kD5ec//7lVw3eWVs47loAA0R+AaKxDUREs2zfDPGQwTC83Q+bIEchPdOw/0QpU4/2GfQFi8AGRk0CogZhMJjWBgsswtBZSUSCyKeCYHSd16MDxQwZ7IB4/flxpTbrxZhpqpQxlAZFdgQSC1rTUDff+GOvJU9RwOFuWE3R00PW5cuWKOkWzgtRUGVhH7lMrZuAkFI5fEpzs0v2v//ovHDx4UF2zWCxKc6L2x1AWEPUEFcqYWjcDZcRuW05aYtDPrg7u/eE1rRVyzI9jewzOgJicnKyuEb7sEuW4rBGILJszhhmMz8ZJP9SW9eQlXtf10XVnty3L1bJRmcifCktAgGgASHmA8fb1otxcWNatgfmNPjAP+dBlMAoQgwuInNzy97//HU888YTqomNjzQZfB87O1F2WbDzZ+OvJJMZrTM+Zpkz7j3/8A7Vq1VJrGnmesNJ56Hx5jhoTu/yoMelu0sGDB4PjfzowHSelMLChZ9clx9Y4o9MY7OvCa/PmzVNjZ5ykw/roiSHMo27dukoLZtcnA4HIdOwG5bgfn5PdiTowf3Z1Mg3rzQk6Gli8R/+AYHp2Uernefnll1G/fn113/Lly1V2BA1nq1I7Zb34/KdOndJFqS3HCVkW60nNVU+6Yb2M3ZeEHuVF2bCrlte5lvSzzz6zdsOyPtRw+Q6Yp/2zcVYx68Dn45baLANnq7J83kNtWkLlJCBA9GMgGgFs2bpJaYyW7Vsr3I0qQAwuIFbuXzz47iIQ7aEdLE9JIHISkATfSECAGCBAJBy5RMPUtxcsG9dVCIoCRAGib5oVz5YqQPSsfKtz7gLEAAKigmJmJkwd2iAvPKxcKAoQBYjB2LhxXJNGAYIxcMyQhhEk+EYC/gnEAPOHZ+zaVPsern9+XCxMPbuqCQulyibg75VfaSB6uP4O62z8YRJg5euJDb75F5ZSRQIiAXdJwD+BaGwcZd+hJpg16RtkzZrm8JoGTqWBKDIvU65avnorQHRXcyT5iAR8KwEBYoA2/kV5uTC1bobCTLPTxluAKF2mvm1epHSRQGBJQIAYoECkdpI9fw6ypk4WIPr4HXpaQ6S1FlqecWegXdKq+EXkcgS9do6m4Gi+TJtL41pAWn/p1auXMmdGk2pcEsBlDb4MXHrB59b19GVdpGz/lIAA0ceNqe52q8y2KD9PWbZxdq9oiMGhIXLROtf4uTNUdaamNg3HOtF+6A9+8ANlNJzHNLPGdXZcHM8F9vTdSK8PevG4O5/DUV5ck6iNaBuv00jAD3/4Q2s9jdfcuc91kFx87yiUdc1RejnnXQkIEAMYiAShedD7yI91bANVgBgcQKRFE72w3F3NgzuByDppizTUlmlfMzY2VlWVM0IJdG8GLmbnbE1HgbL0dKABA2c/YMq65ul6Sf7lS0CAGOBAzFm9ElnjxznsNhUgBgcQaRFGaxw0hE2bmTQKTfNrtPBCDY2unOiVQk/Zpzkxdmmyi5Bpae+UroVoiJpBA5HmzmhM2ujdnvZL6VbJGNjNSIsutI5CU2Z0E0UTYQy0zMKyqAHSYgrPN2vWTJ1jvWhmjddpi5NgZ/60z1mvXj1VD+ZBU2+8l5ZcaIWFgfCiIW1abqFFF5p4Y2D3LOtPG6h8rldeecVq6ox50JYny6Xc7IOuJ8ujD0f6WKR1GlqhoWcPBnYBMx19E9JoOC0D0RwbzzNQpixHB1oDYnoGnqcZOeanz+l09tcoD75PY+DzpqamKrusdAA9cuRIZSGIz8Oucx1oAYfvjXV/7bXXrBoxjXzTaMELL7xQpS5xXU512woQAxyIhaYMmDu2ESD68D16egyRGo/20MB9wobmw7Zs2aK8QrDRJtgIKTbGdNDLQDARXkVFRQo4dJE0fvx4dU0DkQeEnDaVRmDSFqguTyUGlMFqeqggSAi+7t27W4GoPVuwG5UaIcsltOidgpChhsh9GhOno2Ha7mRXKu1u8nloSo150IYnn4XmyRi0BwvuE0L0G8jAujMtHeVSDhybpIcKBtoApc9ClqvtgaoL9/6wbqwny6MvRLp3ohsqOtWlXVjtKJnpaNaONmN5nYCj6ycGlm80Ts76Mz0Dy6R5PZqL09451AUH1/hDgl3NGvSU3W9+8xtln5QWa/iM2l/krFmz8Mtf/lL9cOEYLT1wsBuagT4ntUFwvm+anaN5OG1CTpcv2/IlIED0YUPqbOzP1fOmTu1QkJZaCoqiIQaHhmgPRHqz14GNr9bmqBWyYWYjzsB97Y+Px9Q4aPuTwQhEgpRdeQzfffedja9BdRJQ2ifv14EeGzQEjEAkUHie5xjYdfnAAw/o2xSsqCFSy2UkkD7++GOr/0Rddzb6hDufVaelXVQ6N2bdjR4teJ0Q1IHyctZlyrppIHKfz6HDn//8Z6Vt85jXtD1THhPw9EVI0JQFRKYtq1vU/hp/HHC8lYEar/ZpSCBS+zMGvju+A2rR1AC1XDguSbdUDNSk+/fvb2Mo3ZiH7JctAQFiEAAxc8TnyFm7WoDoo3fpbQ3R2NjbN7BGGHHfaMia+2z0GYxApOZGP4DshuOEGGoc9oFQ0+OCvGacVOMKEOlLkV2JhICO1M50HlpDo9ZEIOo0esvJMsa6sy720HUFiHwOHYw/PCg7o5ZM2fActWN3ApETj9ilzC5rAl9rx45smlIL5I8JdqPSibCWCbfaRyPHcilfvk9tDF0/n2zLl4AA0UeNqKtaYFnpc5YtQebY0QJEH71LfwaiHlNkU0BNUmtW9lBhFyjHGNltqCfIGJsPetowepqnJkdAMGiYES7laYgEK7Ut+6Dz0EBkVyi7brUbKGN6+7p7CohGDxrsniawGJYsWaK6NnWdCE4tC56z/5Gi0zm6xjFVdgXTcwl/rFADZSAQ7Q2YM1/Wg+PGHIMtK0RFRSmNtqw0cq20BASIPmpEywKcq9fyow/B/HpfAaKP3qU/A7FRo0ZISEhQSyDoXFZDzR4qcXFxqtvNmacFOrH905/+BHp6p79AdvX9+7//u2pRNMwqAkTmw4kw4eHhoKNijg1y3FLnoYHIjKnpENT0xUj/fiyXwb7u9kBkV+NXX33lcHkFwcV66vLK0hBZT4KFmjFhRGfMDOy2/dd//VdMmTJFdVvSRZURiOy25IQdTo6xD46u8Xk4XqjdUPEeApHy5RgixyLHjh2rtELKij8WCGemoVyovfIdM3DNJ7VoThbi+5LgmgQEiD5qRF2FXlnpiwoKYGrTTIDoo3fpaSDqho//2tw3+kAcM2aMzXgZu8/0uBgbaQKI3aDUKNg1qQO1X8LGAAAFCUlEQVS76tjY6kBNhesF9QQPfV5veZ2aCWFDGBJSQ4cOVTBjeSyXY5hssI11YF05NmYMHAcjqDlLtnfv3qDzX50H79eB3ZP0dk/tlJrtxIkT1SX7utuXweusJ40D2AddT10e66yDUc6UXUhIiJqsQqMChJYxEDisE2HImZ3MV09m4oxg+kXk+Kh9cHSNAOMYoBHOrAsnChGEnPREOdGHpQ68h9DlWCInWa1Zs0Zd4jNzAhDrZezi1vfJtmwJCBB91IiWBbjKXDO1bY7iwgIbKMqkmuCYVFP2v7Dzq2zUjWOIzlISdmx46TxXQokEKDvjGKKn5EKIcvnEsGHDbIogEJ1p6zYJ5cCtEhAgBgsQe3dHwfEkAaIP3qenNcTK/sdXFIiclfjb3/7WRvOsbJnBcp83gMjlHE8++aTqjr1y5YqN6ASINuLw2oEA0QcNaGU0wPLuMX8wCHnh4QJEH7xPfwUix8mMXYLOWhWOO3GZg4T7EqDs7CF1/6p79jhhiBOdjN3EOueTJ09aF9vrc7L1vAT8E4gB5g+vFKx8UP/M0V/CsmNrCRDvlV/pLlMf1N9GhgFWvr8C0fPNh5QgEgguCfgnEH3wK9+mQQ7A8ukf0fLdWtEQffDuBIjB1SjK01RfCQgQfdCAegK+WXPmInvJEgGiD96nALH6NqDy5MElAQGiDxpQTwAxZ8VS5R/RmHelu0yDRCZGWXhyX4AYXI2iPE31lYAAMUgaf8u6VciaNV00RB+8TwFi9W1A5cmDSwICRB80oJ7QViwb1yF76iQBog/epwAxuBpFeZrqKwEBog8aUAGidxbMe0LOjvIUIFbfBlSePLgkIEAUINpolY4afDlXNsAFiMHVKMrTVF8JCBAFiALEKn4DAsTq24DKkweXBASIVWwMfa49FZ9A5rAPYGrdDKZmDZE5fJgVcDLLtGzNzl3vToAYXI2iPE31lYAAMdCBeOIkcg/sQ0atZ5BR61nkRUUKEL38TgWI1bcBlScPLgkIEL3ceLpLK7HPx/z5MJiHfmiFIa+LhigaYnA1V/I0IgHPSkCAGCRALMrJQcHxFAGiD96naIiebaQkd5GAtyQgQPRBA2qv3XnqWDRE0RC91ZBIOSKBYJCAAFGAaKNVegrOwZyvaIjB0BTKM4gEAAGiAFGAWMVvQIAoTalIIDgk4J9ADDB/eKW0Hz+pf6W7TP2k/qXkWlFwebn+AsTgaAzlKUQC/gnEijZ8kq5M7a7SQBS5lilXe1ALEKUhFQkEhwR8AsTTZ8661ODYN0ByXLHJIgLEismpqt+TADE4GkN5CpGA14GoRX77zh3cunWravH2bbAxcne8c+cOgiFu3boVI0aMCIpn4ftw93tmfrduMVbtO/z+++/1Zy1bkYBIIIAl4DMgBrDMAqbq27Ztw5dffhkw9ZWKigREAiIBX0pAgOhL6Xu4bAGihwUs2YsERAJBJQEBYlC9TtuHESDaykOORAIiAZFAWRIQIJYlnQC/JkAM8Bco1RcJiAS8KgEBolfFLYWJBEQCIgGRgL9KQIDor29G6iUSEAmIBEQCXpWAANGr4pbCRAIiAZGASMBfJSBA9Nc3I/USCYgERAIiAa9KQIDoVXFLYSIBkYBIQCTgrxIQIPrrm5F6iQREAiIBkYBXJSBA9Kq4pTCRgEhAJCAS8FcJCBD99c1IvUQCIgGRgEjAqxIQIHpV3FKYSEAkIBIQCfirBP4/DYErhMu4Y8MAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Calculate convergent Matrix V\n",
    "\n",
    "**Start running the following cell because loading all the data and re-referencing to get CCEPs takes a few minutes.**\n",
    "\n",
    "In the meantime, we can discuss the convergent paradigm (A) used to calculate BPCs and contrast with the divergent paradigm (B) and the interpretation of the convergent paradign (E). \n",
    "\n",
    "![image-4.png](attachment:image-4.png)\n",
    "     \n",
    "    Note: Rereferencing is tailored to these CCEP data, using the custom function ccep_CAR64blocks() which applies common average referencing on CCEP data. This excludes channels with high variance from 500-2000 ms after stimulation onset or high variance from 10-100 ms after stimulation onset (channels with stimulation artifact or large evoked responses). The function also assumes that noise is shared between blocks of 64 channels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "We iterate through all trials with electrical stimulation, no artifacts and not stimulating the recorded electrode of interest.\n",
    "\"\"\"\n",
    "\n",
    "df_cceps = pd.DataFrame(columns = ['ccep_name_1','ccep_num_1','ccep_status_1','ccep_name_2','ccep_num_2','ccep_status_2'])\n",
    "V_pre = [] # single-trial stimulation-evoked voltage matrix\n",
    "\n",
    "for index, row in tqdm(df_include_trials.iloc[:].iterrows(), desc='Included trials', total=df_include_trials.shape[0]):\n",
    "    \n",
    "    df_data = pd.DataFrame(columns = ['channel','data']) # dataframe for mef3 data\n",
    "    samples = [row.onset*srate+element*srate for element in epoch_limits] # get start and end samples\n",
    "    # Read mef3 data\n",
    "    for channel in df_electrodes.name:\n",
    "        # pyMef does not apply unit conversion to microVolt, apply Natus conversion factor 0.2658 here after pyMef\n",
    "        data = 0.2658*ms.read_ts_channels_sample(channel, [samples]) # Returns 1D numpy array with data from sample X to sample Y\n",
    "        curr_channel = {'channel': channel, 'data': data,}\n",
    "        df_data = df_data.append(curr_channel,ignore_index = True) # append to mef3 data dataframe\n",
    "    \n",
    "    # Get stimulated channels\n",
    "    trial_ccep_names = row.electrical_stimulation_site.split(\"-\")\n",
    "    curr_ccep = {'ccep_name_1':trial_ccep_names[0],\n",
    "                 'ccep_num_1':df_channels.index[df_channels['name'] == trial_ccep_names[0]].tolist()[0],\n",
    "                 'ccep_status_1': df_channels['status'][df_channels.index[df_channels['name'] == trial_ccep_names[0]].tolist()[0]],\n",
    "                 'ccep_name_2':trial_ccep_names[1],\n",
    "                 'ccep_num_2':df_channels.index[df_channels['name'] == trial_ccep_names[1]].tolist()[0],\n",
    "                 'ccep_status_2': df_channels['status'][df_channels.index[df_channels['name'] == trial_ccep_names[1]].tolist()[0]],\n",
    "                }\n",
    "    df_cceps = df_cceps.append(curr_ccep,ignore_index = True)\n",
    "\n",
    "    # Delete current electrical stimulation site from dataFrame\n",
    "    curr_indexes = df_good_channels[(df_good_channels['name']==curr_ccep['ccep_name_1']) | (df_good_channels['name']==curr_ccep['ccep_name_2'])].index\n",
    "    df_incl_channels = df_good_channels[~df_good_channels.index.isin(curr_indexes)]\n",
    "    \n",
    "    # Run the custom function ccep_CAR64blocks() in functions/pyBPCs.py\n",
    "    df_data_CAR = pyBPCs.ccep_CAR64blocks(df_data,tt,df_incl_channels)\n",
    "    \n",
    "    # append to single-trial stimulation-evoked voltage matrix\n",
    "    V_pre.append((df_data_CAR[df_data_CAR['channel']==el_interest]['data'].values)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Plot convergent matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Calculate baseline subtraction and plot difference\n",
    "\"\"\"\n",
    "V_pre = np.array(V_pre)\n",
    "\n",
    "# calculate baseline\n",
    "baseline_V = np.mean(V_pre[:,np.logical_and(tt > -.5,tt < -0.05)],1)\n",
    "# make baseline into the same size\n",
    "baseline_Vrep = np.transpose(np.tile(baseline_V,(np.shape(tt)[0],1)))\n",
    "\n",
    "# subtract baseline from V_pre\n",
    "V_pre_baseSub = np.subtract(V_pre,baseline_Vrep)\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 10), dpi=80)\n",
    "\n",
    "fig.suptitle('Convergent Matrix (V)',fontsize= 15)\n",
    "vmin=-250\n",
    "vmax=250\n",
    "im = ax1.imshow(V_pre, aspect='auto',extent=helper.extents(tt) + helper.extents(df_include_trials.index),vmin=vmin, vmax=vmax)\n",
    "ax1.set_xlabel('Time from stimulation (s)')\n",
    "ax1.set_ylabel('Stimulation trial index (k)')\n",
    "ax1.set_title('Convergent Matrix (V)')\n",
    "im = ax2.imshow(V_pre_baseSub, aspect='auto',extent=helper.extents(tt) + helper.extents(df_include_trials.index),vmin=vmin, vmax=vmax)\n",
    "ax2.set_xlabel('Time from stimulation (s)')\n",
    "ax2.set_ylabel('Stimulation trial index (k)')\n",
    "ax2.set_title('Convergent Matrix (V) after baseline correction')\n",
    "\n",
    "fig.subplots_adjust(right=0.8)\n",
    "cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "fig.colorbar(im, cax=cbar_ax)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Calculate BPCs\n",
    "\n",
    "Identifying basis profile curves (BPCs) that group characteristic shapes in the convergent CCEPs. \n",
    "\n",
    "The input for the BPC calculation is the convergent matrix (V_pre, numpy array), the time vector (tt, numpy array), and the pandas dataframe pair_types that we will generate in this section.\n",
    "\n",
    "## 4.1 Select time-frame for BPC extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# select the epoch times to enter in the BPC analyses in seconds\n",
    "BPCs_epoch = [0.015,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Calculate the significance matrix\n",
    "\n",
    "To calculate the significance matrix, we project the unit-normalized stimulation trials into all other trials. We then calculate t-values across all subgroups of stimulation pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# V: single-trial stimulation-evoked voltage matrix (Dimensions of V are with K total stimulation events by T total timepoints)\n",
    "# tt_BPCs: BPCs time vector\n",
    "V,tt_BPCs = pyBPCs.bpcVoltage(V_pre_baseSub,tt,BPCs_epoch) # take only Voltage at BPC times\n",
    "\n",
    "# You can compare results to taking the V_pre matrix before baseline correction, just to test\n",
    "# V,tt_BPCs = pyBPCs.bpcVoltage(V_pre,tt,BPCs_epoch) # take only Voltage at BPC times\n",
    "\n",
    "# pair_types: structure with subgroups of stimulation pairs, with fields\n",
    "#    ccep_name_1, ccep_name_2: electrodes associated with the pair, cathode-anode ordering\n",
    "#    indices: indices of subgroup in V matrix (out of K total stims)\n",
    "\n",
    "pair_types = df_cceps.groupby(['ccep_name_1','ccep_num_1','ccep_name_2','ccep_num_2']).size().reset_index().rename(columns={0:'count'})\n",
    "pair_types = pair_types.sort_values(by=['ccep_num_1']) \n",
    "pair_types['indices'] = ''\n",
    "\n",
    "for index, row in pair_types.iterrows():\n",
    "    indices = df_cceps.index[(df_cceps['ccep_name_1'] == row['ccep_name_1']) & (df_cceps['ccep_name_2'] == row['ccep_name_2'])].tolist()\n",
    "    pair_types['indices'].loc[index] = indices\n",
    "\n",
    "pair_types = pair_types.reset_index(drop=True)\n",
    "# pair_types.head()\n",
    "\n",
    "### calculate\n",
    "V0 = V/(np.ones((V.shape[0], 1)) * (np.sum(V**2,axis=0)**0.5)) # normalize (L2 norm) each trial\n",
    "P = V0.T @ V # calculate internal projections\n",
    "\n",
    "tmat = pyBPCs.nativeNormalized(pair_types,P)\n",
    "\n",
    "# visualize the set of native normalized single stimulation cross-projections \n",
    "plt.figure(figsize=(10, 8), dpi=80)\n",
    "plt.imshow(tmat,aspect='auto') # Plot a matrix or an array as an image.\n",
    "plt.clim(0, 10)\n",
    "plt.ylabel('stimulation pair subgroup')\n",
    "plt.xlabel('stimulation pair subgroup')\n",
    "plt.title('Significance matrix '+r'$\\Xi$',fontsize= 15)\n",
    "plt.colorbar()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.3 Iteratively decrease inner components of non-negative matrix factorization\n",
    "\n",
    "Using Non-Negative Matrix Factorization (NMF) on the significance matrix to cluster sites that produce similar measured responses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Find two non-negative matrices (W, H) whose product approximates the factorize, \n",
    "non-negative and rescaled matrix (t0_scaled). \n",
    "At the end, the matrix H has size of number of clusters by stimulation pair sub-groups.\n",
    "\"\"\"\n",
    "\n",
    "#### Default settings and initialization\n",
    "\n",
    "# number of cluster dimensions\n",
    "cl_dim = 10\n",
    "    \n",
    "# number of iterations to re-run NNMF - because can get suboptimal factorization due to non-unique nature of NNMF (starts out with a random matrix)\n",
    "num_reruns = 20\n",
    "\n",
    "# convergence threshold: 1e-5\n",
    "conv_thresh = 0.00001\n",
    "\n",
    "t0 = tmat\n",
    "t0[t0 < 0] = 0 # factorize t-matrix, but first make non-negative\n",
    "\n",
    "# Globally rescale data to avoid potential overflow/underflow\n",
    "t0_scaled = t0/(np.max(t0))    \n",
    "\n",
    "nnmf_xcorr_score = 100; # start off-diagonal penalty score out with high value\n",
    "\n",
    "WH_struct = pd.DataFrame(columns=['W','H','nnmf_xcorr_score']) # saving structure\n",
    "\n",
    "#### Do the NMF\n",
    "\n",
    "# while nnmf_xcorr_score>.5:\n",
    "while nnmf_xcorr_score>1:\n",
    "  \n",
    "    cl_dim = cl_dim-1 # reduce number of inner dimensions\n",
    "    \n",
    "    if cl_dim == 0:\n",
    "        break\n",
    "\n",
    "    # multiple run-throughs of  nnmf\n",
    "    tmp_mat_W = []\n",
    "    tmp_mat_H = []\n",
    "    tmp_err = []\n",
    "\n",
    "    for k in range(num_reruns):\n",
    "\n",
    "        model = NMF(n_components=cl_dim, init='random', solver='mu',tol=conv_thresh,max_iter=10000)\n",
    "        W = model.fit_transform(t0_scaled)\n",
    "        tmp_mat_W.insert(k, W)\n",
    "    \n",
    "        H = model.components_\n",
    "        tmp_mat_H.insert(k, H)\n",
    "\n",
    "        rec_err = model.reconstruction_err_\n",
    "        tmp_err.insert(k, rec_err)\n",
    "\n",
    "    # select factorization with smallest error\n",
    "    k_ind = np.argmin(tmp_err)    \n",
    "    W_min = tmp_mat_W[k_ind]\n",
    "    H_min = tmp_mat_H[k_ind]\n",
    "\n",
    "    # Normalize so rows of H have unit norm\n",
    "    Hnorms = np.sqrt(np.sum(H_min.T**2,axis=0))\n",
    "    Hnorms_ = (np.reshape(np.tile(Hnorms.T, t0_scaled.shape[1]),(-1, Hnorms.shape[0]))).T\n",
    "    H_min = H_min/Hnorms_\n",
    "\n",
    "    # Score for this decomposition - off diagonal element weights\n",
    "    HHT = H_min @ H_min.T\n",
    "    HHT_triu = np.triu(HHT,1)\n",
    "    HHT_triu_reshaped = np.reshape(HHT_triu,(1,-1))\n",
    "    # nnmf_xcorr_score = np.max(HHT_triu_reshaped) # max of nonnegative matrix factorization component interdependencies\n",
    "    nnmf_xcorr_score = np.sum(HHT_triu_reshaped) # sum of off-diagonal nonnegative matrix factorization component interdependencies\n",
    "    print('Inner dimension: ', cl_dim, ', off diagonal score: ',nnmf_xcorr_score)\n",
    "\n",
    "    # save matrices and scores for plotting later\n",
    "    WH_struct = WH_struct.append({'W': W_min, 'H': H_min,'nnmf_xcorr_score':nnmf_xcorr_score}, ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.4 Winner take all\n",
    "\n",
    "Each stimulation pair subgroup is assigned to one cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "winner-take-all on columns of H and then threshold by 1/(2*sqrt(length(pair_types))) \n",
    " -- since all equal would be 1/sqrt(length(pair_types))\n",
    "\"\"\"\n",
    "\n",
    "H0 = 0*H_min\n",
    "\n",
    "for k in range(len(pair_types)):\n",
    "    k_ind = np.argmax(H_min,axis=0)[k]\n",
    "    H0[k_ind][k] = H_min[k_ind][k]\n",
    "\n",
    "H0_ = H0>1/(2*np.sqrt(len(pair_types)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.5 Identification of Basis Profile Curves \n",
    "\n",
    "BPCs are identified from the clustered groups (rows of H) using linear kernel PCA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "output will show the subgroup numbers clustered in each BPC\n",
    "\"\"\"\n",
    "\n",
    "B = []\n",
    "B_struct = pd.DataFrame(columns=['curve','pairs']) # saving structure\n",
    "\n",
    "for q in range(H0.shape[0]):\n",
    "    cl_pg = np.where(H0_[q]) # cluster pair groups\n",
    "    cl_inds = [] # cluster indices (e.g. concatenated trials from cluster pair groups)\n",
    "    for k in cl_pg:\n",
    "        for i in (pair_types['indices'][k]).values:\n",
    "            cl_inds.extend(i)\n",
    "\n",
    "    V_tmp = V[:,cl_inds]\n",
    "\n",
    "    E = pyBPCs.kpca(V_tmp)\n",
    "    B_tmp = E[:,0] # basis vector is 1st PC\n",
    "\n",
    "    if np.mean(B_tmp.T @ V_tmp)<0:\n",
    "        B_tmp = -B_tmp\n",
    "    \n",
    "    B_struct = B_struct.append({'curve': B_tmp, 'pairs':cl_pg[0]}, ignore_index=True)\n",
    "    \n",
    "# pairs not represented by any basis\n",
    "excluded_pairs = np.where(1-(np.sum(H0_,axis=0)))\n",
    "\n",
    "\"\"\"\n",
    "Calculate statistics for each basis curve, eliminating those where there is no significant representation\n",
    "\"\"\"\n",
    "\n",
    "B_struct = pyBPCs.curvesStatistics(B_struct,V,B,pair_types)\n",
    "\n",
    "\"\"\"\n",
    "Calculate projection weights\n",
    "\"\"\"\n",
    "\n",
    "B_struct = pyBPCs.projectionWeights(B_struct)\n",
    "\n",
    "B_struct.pairs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Visualize BPCs\n",
    "\n",
    "## 5.1 Plot Calculated BPCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = iter(cm.tab10(np.linspace(0, 1, 10)))\n",
    "\n",
    "plt.figure(figsize=(10, 8), dpi=80)\n",
    "for q in range(len(B_struct)): #  cycle through basis curves\n",
    "    plt.plot(tt[tt_BPCs[0]:tt_BPCs[1]],B_struct.curve[q],color=next(colors),label=q)\n",
    "    \n",
    "plt.xlabel('Time from stimulation (s)')\n",
    "plt.ylabel('Normalized weight of BPCs')\n",
    "plt.title('Calculated BPCs',fontsize=15)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.2 Spatial representation of the BPCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "marker_colors = [] # marker_color array\n",
    "marker_sizes = [] # marker_size array\n",
    "\n",
    "# all electrodes\n",
    "xyz = df_electrodes[['x', 'y','z']]\n",
    "xyz_list = xyz.values.tolist()\n",
    "electrodeName_list = df_electrodes[['name']].name.values.tolist()\n",
    "for cords in xyz_list:\n",
    "    marker_colors.append('white')\n",
    "    marker_sizes.append(7.5)\n",
    "\n",
    "# el_interest\n",
    "xyz_el_interest = df_electrodes.loc[df_electrodes.name == el_interest_nr.name.values[0]][['x', 'y','z']]\n",
    "xyz_list.append(xyz_el_interest.values[0].tolist())\n",
    "marker_colors.append('black')\n",
    "marker_sizes.append(20)\n",
    "\n",
    "pair_types['interpolated_locs'] = ''\n",
    "\n",
    "for index, row in pair_types.iterrows():\n",
    "    xyz_1 = xyz.iloc[pair_types.ccep_num_1[index]].values.tolist()\n",
    "    xyz_2 = xyz.iloc[pair_types.ccep_num_2[index]].values.tolist()\n",
    "    pair_types.interpolated_locs[index] = np.mean((xyz_1,xyz_2),axis=0).tolist()\n",
    "\n",
    "# non-significant stim pair sites    \n",
    "interpolated_locs_list = pair_types.iloc[excluded_pairs[0]].interpolated_locs.values.tolist()\n",
    "for cords in interpolated_locs_list:\n",
    "    marker_colors.append('gray')\n",
    "    marker_sizes.append(7.5)\n",
    "    \n",
    "xyz_all = [y for x in [xyz_list, interpolated_locs_list] for y in x]\n",
    "\n",
    "# plot BPCs, colored\n",
    "colors = iter(cm.tab10(np.linspace(0, 1, 10)))\n",
    "\n",
    "for q in range(B_struct.shape[0]):\n",
    "    # get electrodes for BPC\n",
    "    xyz_BPC_list = pair_types.interpolated_locs[B_struct.pairs[q]].tolist()\n",
    "    xyz_all = [y for x in [xyz_all, xyz_BPC_list] for y in x]\n",
    "    plotweights = B_struct.plotweights[q]\n",
    "    curr_color = next(colors)\n",
    "    for ind_cords in range(len(xyz_BPC_list)):\n",
    "        marker_colors.append(curr_color)\n",
    "        marker_sizes.append(B_struct.plotweights[q][ind_cords]*25)\n",
    "\n",
    "# plot all\n",
    "view = plotting.view_markers(marker_coords=xyz_all,marker_size=marker_sizes, # marker_labels=electrodeName_list,\n",
    "                             marker_color=marker_colors, title='Spatial representation of BPCs rendered on an MNI brain') # Insert a 3d plot of markers in a brain\n",
    "view.open_in_browser()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5.3 Optional parameters to change\n",
    "\n",
    "Once you have completed this, you can go back to Section 4.1 and change the time interval over which the BPCs are calculated (e.g. 0.2 - 1 sec) and look at the effects on the outputs.\n",
    "\n",
    "You can also select a different electrode in Section 3 to look at the various inputs into different regions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bpc",
   "language": "python",
   "name": "bpc"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}