fig_4_invivo_invitro_logos.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-01T22:16:27.476975Z",
     "start_time": "2021-04-01T22:16:26.574747Z"
    }
   },
   "outputs": [],
   "source": [
    "# Standard imports\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import logomaker\n",
    "import time\n",
    "import glob\n",
    "%matplotlib inline\n",
    "\n",
    "# Import marginalization function\n",
    "from utils import x_to_ct_df\n",
    "\n",
    "# For drawing rectangles\n",
    "from matplotlib.patches import Rectangle\n",
    "pink = np.array([247, 198, 219])/256\n",
    "purple = np.array([126, 66, 149])/256\n",
    "red = np.array([228, 33, 63])/256\n",
    "\n",
    "# Load figure style\n",
    "plt.style.use('fig_style.txt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-01T22:16:35.842366Z",
     "start_time": "2021-04-01T22:16:27.478730Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6-7: ..\n",
      "7-8: ..\n",
      "8-9: ..\n",
      "Execution time:6.65 sec.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-2-b3f3707e7b8d>:124: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.\n",
      "  fig.show()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEECAYAAADHzyg1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABpz0lEQVR4nO2dd5hU1dnAf+/OzO7O9srCLmXpnZUmYKEIikYRG7HF2CKWaIwm9kSJmsSYfDGWaGJsKJbYIaggRQSkSG8CLmVZ2u6yvc3OzM6c7497d9i+s31hz+955plz7z3lvXfeO++p7xGlFBqNRqPRAAS0twAajUaj6Thoo6DRaDQaH9ooaDQajcaHNgoajUaj8aGNgkaj0Wh8aKOg0Wg0Gh/aKGgaRESeEhElIrHNzMcuIk4RebGlZNN0DlpLB0XkZhG5o2WkPD0QvU5B0xAiEgMEK6WONTMfK9AVyFdKFbeIcJpOQWvpoIgcBdYrpa6oJW6AUsrbnPJORXRLQeMP3wOvA4jIUhHJEpFPRaRIRN6rHllEvhKRJWa4v4i4ReR+4BrgMDC0rjgiEi4i/xaRAhHJFZH/a7vb1HRgWkMH9wOJwOVmK6SniDhEZJWIbATeEZEEEflIREpEJENEftNmd9xOaKOgqRcRiQD6AFvMUyMBB/A0sBGYVUuyzUCKGX4a2A+8aKb1ANvriTMPmAlcDLwF3C8ik1vujjSnGq2og0+Z128EugORQDBQDvwceBT4ChgGTARWAH8Tkd4tdnMdEG0UNA1xBiDAFhHpCcQA/1ZKbcbQn321pNkExIvIJRgv7P1KKTfGC7lXKeWoLQ7GS3kp8LZSarUZB8DSWjenOSU4g9bRwe5m3OVKqaPmNYBHlVI/AHHmuReVUpuAHeb10/p/87S+OU2LUPGibKkU/l5EBONl3VRLms3m9xvAIqXUl+bxGZys7dUWRwFejKb9GcBvMV741S1xI5pTltbSwYGAGwgWkWAz7/JKaSvGE0aKyDjgdmClUmp/c2+oI6ONgqYhRgJFGM3vihdyK0ZzPoJaXkilVBqQC0QB9wGISDIQbaatNY5SKgf4FTAKWAtkAxcrpZwtfE+aU4tW0UHgGwyjkAoMN/PeU6FvSqktwB+BK4BlwAbg6pa8sY6Inn2k0Wg0Gh+6paDRaDQaH9ooaDQajcaHNgp+ICLjRGSNiKwWkefaWx5N50Xroqa10UbBPw4B5ymlzgG6iMjw9hZI02nRuqhpVaztLcCpgFIqo9KhG2Pxi0bT5mhd1LQ22ig0AhEZAcSbC1uqX5sNzDYPR+tZXZo6kBbJpA5d1HqoaQS16qKekuonpkOuz4GfVqut1RZX6eeqqYNmGwV/dVHroaYBatVFPabgB6ZnxXnAbxsyCBpNa6J1UdPaaKPgH7OAscCzIrJCRCa0t0CaTovWRU2r0mm6j0Qk1nSjUHEcpZTKb6WydLNdUxctMqbgV0FaDzX10+m7jz6qdvxqu0ih0Wg0HZjTfvaRiFyO4dBqiIi8bZ62YjjS0mg0Gk0lTnujACzHcIWbDvwbo8nkBvQgnUaj0VSjM40pJGG0GKIw+9KUUk+2Ulm6L7edcCxfi3P9VvA0vLWuBNqwzziPwMH9Wl+wSsW2WUFaDzX107nXKYjIBuDvgG/jb6XUt61U1mn5MuY/829KPvgfuMsbjmwJIOis0cS+8iTGXiitj3K6ODpkOjTi2dtGDCLh83+1olQ10EZB01GoVRc7Q/dRBYeUUu83NbGIJAILgSFAmFLKj3/G0wv37n2owmL/4+/60fiDbiOjIEGBRNx/C861W8BrtBS8pWW4t+8xIlgtBI2p5Coo0EbYNTPaRLaWQuuhprXpTC2FlUAo8APGto8opX7eiPTBgB34DJhW38t4utbQVHk55QePoNxu37msK+4Cl3Ec99azBMRFAyABFqzJSUhwULvIWoF7/yEyz78RgIDYKBI3fN6u8tDMloLWQ00L0ulbCjc0J7FSqgwoa6uukI6IWK3Y+idXPRcQQMXfjrV/MtZuXdpcrs6E1kNNa9OZ1iko4BGMPVePAOe3ZOYiMltENorIxpbMV6NpDFoPNc2lM7UU3gDuBl5WSnlE5FrgtZbKXCn1KuaCOBHRbXZNu6D1sOPi2vkjZSvWoco9BI0eTvC5Y9pbpFrpTEbBopTaU6nZ3ZlaSRqNph0p/XIFufc+CRYLEhRI0Qtzibj3JiLuvam9RatBZ/pjXC4i/wISReR5YEljEouITUSWAinAYhEZ1xpCajT1ofXw1MOxZDW5v34KPF5i/vowca/9GYICKXz+LQqff6u9xatBp2kpKKWeMrcuXAbsVUptb2R6NzCtVYTTaPxE6+GphWPFenLungPlHiIfvYuQGVMBiH3hcXLufNxnFDpSi6HTGAURsQN9gTBgjIiMUUq90c5iaTSaJqCcLrzFpX7FFUsAEhGGBLRtx0jZ6o3k3PE7cJdj6RZPybvzKXl3vu+6pXtXPOnHOpxh6DRGAVgMfAkcb29BOiN/mvMkjgL/Fr7ZI8N4dM7jPPirR8hMz/YrTULPOJ594c/NEVFziqBcbgqeexNvTp7vnDc3n7Jv1gEgoXbsF046mSAggOAJIwm57Pwm6SHQaF384713kDP7MXC5sV9yHjHPPYZYLDXiFr7yLoV//Q+Fz7+FpVsXQn/6E7/KaE06k1FwKKWeaW8hOiKlpaWsW7eOnJwcRISuXbty5plnEhgY2HBi60lFF2tNpa/AUVDMo7Nu80ueP330HwAy07NJLjrDrzRp6Vtrv1BlPn/Dc/u9Xi95eXkUFhbidDrxeDwEBAQQHBxMZGQk0dHRbea2Q1M7Emgj6uHbq5xzbdvtMwoBcdHE/PXhWtM2RQ+hCbpotWIbPRxKS/GWlJJz++/qjB84dgR4vUh4mF/5tzadySi8KSKfATs5uaL58fYVqXYcDgc7d+4kLS2N8PBwCgsLCQ8PZ9iwYXTv3r1F/5R2797NF198gdPp5IILLiA/P59ly5axadMmZsyYQZ8+fepNb79wEqUff0Xg6GEExEa3mFwthSUhDgm1o0ocWPv1rDNeeno6a9as4cCBA0RERNCnTx+ioqKwWCy4XC4yMzPZs2cPgYGBnHHGGZx//vmICEVFRezZs4djx46Rn5+PxWLBarWilMLpdOL1eomOjqZHjx4MHz6coKD2XeGtaRusPROJ+dNvqnRxOb78hqJ/zgMg6JwxRD5yp++aBFqx9khsczlro1WMgoicAbwEjATygauVUqsrXb8JeLNasjeVUre0hjwmDwJ/pZJDvLZAKUVGRgY5OTkopQgLC6Nnz55YqjUlvV4vBw4cYOvWraSmphISEsLQoUNJSkrC4XCwbds2duzYQUhICCNGjCAlJYXIyEhfeo/HQ2pqKjt37sTpdBIfH8/Ro0dJTExk+PDhdOvWrYoxKS0t5auvvmLnzp0AjB8/HqvVSlxcHMOHD2fHjh288847jBw5kgsuuIDg4GCUUrg27KDsu02Upx5EOV0ExMcQeu0MCA4i//fP4TmRCwEB2Ab2xn7eBAJTBrfNg64FpRSuLT9gv3QanvRjBI4cQvmho1h7JVWJd+jQId566y0ARo8ezcUXX1yr4S0sLOTQoUOEhYUhImzatIkvv/wSr9fLueeey7Rp0wgNDa2SJjs7m82bN/sM/MCBA1vtfjs1gTZfUPxp4bYyIoK1R7cq51ybd/rCARFhBA7u2yJlKaU4fvw47kruZxoiIiKC6OjaK3EtbhREJByj/z4VmAAkA9VHhD4Elprhy4EXaOQU0Sawv5kO8Z4DxgCblVL3NhTf6/Xyww8/sGbNGvLy8pgwYQLh4eEsXLgQt9vNmWeeyejRo7Hb7QCkpqayZ88eRIQpU6bQrVs3AgMDERFiYmJITEykvLyc9PR0cnNz2bhxI2effTY5OTls27aNXbt2ER0dTXx8PEOGDMFisRAbG8vevXv56KOPsFqtDB8+nBEjRlBaWsp7771HSUkJI0eO5Mwzz6wie7du3Rg/fjyrVq1iy5Yt7Nu3j+uvvx7rn17F8b9lAET94deEXn9pjcE7pRSujTtw7Uql/OARbCMGtUt3i2v3Pgr++DLONZsJueoiwm6ZRcFf/k3Rv94j7PqZhN/zcywxUQDk5Z3sm46Nja0ib4Ux9nq9eL1ejh49is1m45JLLmHjxo14vV6sVisTJ07EajVep/LycpYuXUpRUZHv2OPxsGfPnhYxCo3VxbbC4XCQk5NDSUkJHo+HoKAgrFar73l6PB7KysoQEUJDQ4mPj8flclFYWEhlH0379++nrKwMgOTk5CqG1maz0aVLF7zHsyhbuYHy1DSU04Wle1cifnMryuXCEhNN8duf4snOw5OVgzUxgeBJZ7ZrBcVfioqKKCgowF+fVRXPo1u3blXSHDt2jNdffx2AqKgo7rnnnirpAuoZdG+NlsLFQBfgKnPaZ42pn0qpUkxDISJXApnAJ60gS2USRGQThkM8Uwz/HOKJyCgMj5TnisgrIjJWKbWhvjQvvvgiZWVlDBo0iOHDhxMcHAzAWWedxfHjx9m6dSsrV64kJSWFs846i8jISLp06UJJSQkul4u8vDzy8vJwuVwAxMTEEBgYiNVqJTo6mrAwo//xwIEDOJ1O+vTpQ9++fQkPD8dqtRIQEEBISAhxcXGUlZWRkZGB2+0mKyuLmJgYBg8eTEFBAXa7nfz8fIKCgrBYLJSWllJSUoLb7SYuLo4xY8YQGxtLdHQ0ZYP74vh6FThdlB86iud4FpbEBN9Lr7xeyg8do2z5WpwbtmPt3QP7zLadPek5kUPB/71O6UdfERAdQfDUs7D27o57z35CLp6CY9G3lPz3C0o+WUz4ndcRfvNVjBgxguzsbNavX8+yZcsoLCykb9++REREkJycTI8ePSgpKSEzM5NDhw5RXl5OaWkpP/nJT1iwYAHZ2dm8++67jBkzhm7dumG32znvvPNwuVwUFRWxfv16tm3bhogwc+bMZt1fU3Sx4s86NzcXESEqKgq3201BQQFBQUHExMT49Ksy5eXl5OXlkZub62vlFhUVYbFYiI6OJjo62mcICwsLWbx4MdnZ2YgI55xzDgEBAVgsFgICAvB6vXg8HrxeL8ePH8dqtWK324mLiyM4OBiHw4Hb7cbj8ZCdnU12djYBAQH07t3b1+Vms9mw2+2oMidZl92BNycfCQ6i6/J5WLrG17hv5fHgXL0R954DeEvLmvXc24rw8HDCw8OrnPv3v/9NRoaxJ9gNN9xQo0vX6/Wyd+9esrOzKSkpITExkeDgYC644AKUUtjtdo4cOYLL5aKkpASLxUJiYiIxMTG1ytAaRqGX+f2MiPQHtgG3KaXSqkcUkaHAJOAppZSrFWSpTHMc4o3nZEtmKUYLqN4XcfTo0cTExBAdHU1JSQkZGRl4vV7sdjsjRoxg8ODB5ObmAkYNqmvXrnTt2hWAsrIy8vLyOHToEDk5OYDR3IuLi6sy2FmRX1BQUJVamtvt9v1RBwUF4Xa7fTU0r+lSOjAwkKCgIAoKCny13eovr9frJT8/n9LSUrp06UKfO64j7PqZODdsx/3jQYpe+5Athw6SXlJIUkQko+MTCbAHY+kWT8R9txA0ZlibthKUx0PxGx9TUlRM1i8up/jMYRSXGX82FosFEaH8thkEOpwErd5C3IbN9ImKIPzaGUybNo2JEyeSlpZGeno6GzZswOVy4fF4sFqtREREEBsby1VXXUXXrl2JiYnB6XRyww03kJ2dTUZGBmlpaWzcuBGPx4NSivDwcOx2OxEREcycOZNu3bo1fBMN0yhd3LhxI7m5uQQGBhISEoJSipKSEgCsVislJSUcP34cp9NJdHQ0w4cPJyMjg/T0dEpLS4mNjSUiIgK3243L5UIphc1mIz09nW3btmGz2UhKSqJ3795MmTIFh8OBy+WivLycI0eOcPjwYQCio6NJTk4mKCiI/v37Y7fbCQ8P5+jRo6SmplJYWEhoaCiRkZEMGDCAlJQU3G43JSUlvgF/j8dDREQEffv2peu37+PatBP3j2kUzf2U3aX5HCjMp9zrJcJuZ2T3ZKJdHsRmw5qchK1v3eNJrY3b7ebgwYOkW1wcG5mM2+kkrFsoMUuW0K1bN/r16+erNGZnZ5Oenk5RUZHvNxs7dixerxelFEVFRezcudNXWYyKiqJHjx5V3t+CggIKCwux2+0EBgbicDg4duwYVqsVr9dLYGAg5eV1e1xvDaOQY36vw3BAtxx4Q0TOMs8/rZR62gz/EijH2CaztXEDV1Jp5zXA353XooADZrgAGFo9gojMBmZXHC9durR6lGaRnp7ud1xfzb2WJuiaNWualG7Xrl01zi1evJh169b5jlNSUpg5cyZyJBuO7IYl//NdS0tL81v+tLQ05syZQ1paGsmxZ/id5ok//IF58+Zx8OBBRISrr76aXr16ERQU5Ls3r9dLcXExuw/s5sSJE/TIO0DK3k1+y9Zc5syZ09wsoqhHF6vr4cKFCxuV+eLFi5snXT2sX7++RfPLyclhx44d7NixA4/Hw8iRI+nSpQs7duxg37599OrVi2HDhjFo0CCCNq0CmqaHFeHG6OKcOXPYs2cP33//Penp6cTFxTFw4EBiYmIICgvCs/V78ld8zb59+0hPTycpKYnzzz+fHj16NPIpNJ26dLE1jMISjD9gF1CGMdPnECcVNRd8Yw8/A+YrpY62ghzVmY+x81qjVjKbFAARZjgCY/C8CtUdkVU00++8885ap3Z6PB52795Nfn4+ffv2Zfv27b4/2BtvvJHk5GRf3BMnTlBcbMytdrvdFBYW0qVLF3r06MGhQ4c4ePCg71xYWBgOh8NXowsKCqK4uJji4mLCw8MZOnRojcFQqPoHP2nSJKZMmVLvA3nllVdYt24dgYGBvPzyy9x9991s27aNyy+/nCeeeKJG/N/f92C9+VUmOTmZOXPmcHDrbVDkf5onn3ySRx55hK1bt3LgwAHy8vKwWCwcOHCA/Px8wOiGS05OZtCgQSQlJTFy5EgSEhJwOp1VjOE333zDpk2GsRg3bhznnnuu75qIEBwcTGZmJnv37iUvL8/XlWKxWHy1uvLyctxuN/n5+YSFhbXUIHO9ulhdD1vACPlFeXk5P/74I5mZmbjdbhISErDb7VgsFt9YQkXr8/Dhw1itVvr3709ISAg//PADeXl5JCQkEB8fT2BgIAEBAWRlZeFyufB6vRQVFZGbm0tISAgXXHAB7777Lrt37+acc87hoYceIi8vj/z8fFwuF2eddRbBwcGUlJSwY8cOjh49ym233ca5557bJD0EGq2L06dP5w9/+AMAY8eOZd26dXX24x86dIjs7GzDeAUFcfz4cV/XcVhYmG9ssXLFxuFwUFJSQnBwMH379uXVV1+lsLCQmJgYbr75Zl/3ssPh4NChQ3g8HsAwWCUlJXTv3r1O+VvcKCilDonIjcDTwD0YRuJRpVT1RWM/B8KBf7a0DHXQnJ3X1gK3YwyQTwPeaijBrFmzyM7OZuXKlVV+WKUUbreb4uJibDYbsbGxvvGECg4ePEjPnj19SpSfn+8btMzKyvLVuG688UZ27dpFWloabrebxMREwsLCiI+Px+PxcPToUd8YRVpaGqWlpWRlZXHJJZeQlZXFwYMHOXHiBB6Ph9DQUMaPH09AQABWq5VvvvnG1+fctWtXBg8e7GviLl261Ddw9fzzz3PrrbdisVi4+eabmTNnDn379uVnP/tZEx9187Db7UyYMIFu3brxm9/8hk8//bRGnJCQEH7zm99wyy23+KaPFhYWcvjwYfLy8ggICCA8PJwJEyb4BkV37NhBYWEhAQEBvlpfWloaO3bsIDc3l5EjRzJ8+HCioqJ8v7XL5SI/P5/Dhw+Tm5tLbm5ug1N8/aDRutgWWK1WhgwZwpAhQxqMO2rUqCrHdf1B9e1b9+ycGTNmsGPHDt544w3+85//cOmllzJo0CAiIyNxu90cOHCAzz//nP379zNp0qSW6rrzm0GDBnH22Wfz3XffsW/fPhYsWMD06dN9E0vA+MNetGgRzzzzDMePH+f//u//mDVrFj169Gh0i+Hee+8lMzOT7OxsduzYgdVq9RlXwGeQ4+Pj6dGjB/HxNcdgKmiVKanmn2+9f8BKqX/SBgZBRN7BaK10qTTQ3Kid15RSm0WkTERWAVuVUt83lMafl6MyI0eOJDExke3bt5OWlsauXbuIjIwkODgYm83m69OtGFTu0qULXbp08bUolFK4XC5ff67X66Vv376+P/lzzjnH92e1a9cuPv74YwAGDx7MVVddVWstxuPxsG/fviq16B9//JFZs2bh8XgYPXo0Q4YMYdWqVfTt29f3Etxyyy306NGDSZMm1cizLdixYwdTp07lxIkTAFgsFnr06EFRURE5OTmUlpby1FNPcezYMV577TVSU1N5/31DXfv168e1115b6/NwOp3s2LHD188dHh5ORESE77ioqAgRwWYzpke6XC4KCgp8Lb3Y2Nhm31tTdPF0o6ysjHHjxpGRkcHgwYNZv359lcHZ/Px83G43999/P+vWrcPj8ZCUlFRPji1PdHQ0q1atYtOmTSxfvpz58+fz8ssvV6ntezweevbsyVVXXcW0adMYOXJkk8sLCAigW7duLWL8OsPitbqXEjaCtpj6l5CQwPnnn9z7RynlG/S1WCw11jZURkQICgrya3HUoEGDuOyyyzhw4AC5ubl88cUXREREYLPZfIPRTqeT7OxsSktL6datGz169KCsrIxLLrmE/Px8unTpgt1u57HHHqsiQ1JSEkePHuXyyy9n7dq1bT4vv6ysjFmzZvkMwu23387TTz9NXFwcSim+/PJLrr32WoqKinj99de5/fbbGTt2LPfeey9HjhwhNzeXb7/9lrCwMN90Sq/X65u5YbfbSUhIICwsjGHDhjFs2DDfQrXKM2gAQkNDiY6OZsiQIf6tDveTjjQNtT2o+D0qqPxepKam8txzz+F0OgFYtGgRx44d44ILLmjV8ZLaEBGcTifbtm1j6dKlZGVlER8fT0REBIWFhZw4cYKwsDByc3M566yzGs6wjTjtjYJS6hCAiNyjlHrRDAvGhjsvtqdsDSEivil/LYnFYiElJYWUlJRGpbvssstITU0lPDyc5cuXM3RojfF2Dh06xJlnnklWVhaXXnopu3fvrndOdEvzwQcfsHfvXgCuueYaXnnlFV/tTES4+OKLefXVV7n22msBYx3C2LFjiYqKIioqqkllVowxVHSvaVqXkJAQ1qxZw5w5c1i4cCFDhw5lxowZDBkyhMjISCZOnEh+fj6HDh1CKUVQUBC9e/duczmfeOIJnnzSmMvy61//mvvvv79Kt9D+/ft58MEH+fTTT1m8eDFHjhypt1unrTjtjUIlLsc0AkopJSK+Y41//O53v+P8889n8uTJ9OnTh6+//podO3aQnZ2N2+0mKiqKIUOGsGTJElatWsXAgQN9BsEeGVbFl0x92CONQbKEnnF1+zSqRkLPOMAwChX84Q9/qHVK7NVXX82WLVtITU3liiuu8Ct/Tcdi4MCBvP/++7jdbvbt20daWpqvazAgIICoqCj69evHLbfcQu/evX2Vq6boITReF/fs2eMzCFdffTXPPfdcjXh9+/bl448/JjU11beYr0OglOoUH+AboL8Z7g+saMWy1OlKSUmJevDBB1VYWJjCGJtRoaGhKiQkxHccHBysHnvsMVVeXt7m8kVFRSlA9e7du83L9pO21Pm2uidNNb766ivf+/DEE09UufbOO++oAQMGqD59+vg+w4YNUx988EFbi1m73tR14XT7AAOBT4GVwMfAwFYsq1G/zKlCWVmZOuussxSgLBaL+t3vfqfS09OVUkp5vV71448/qiuvvNL3MjzzzDNtKp/D4fCVfd5559W4fvToUbV9+3bfZ9++fcrr9bapjEobhU5Bfn6+iouLU4Dq2bOnys/P911zuVxq7969Pj2Mjo5WgJo9e3Zbi1mr3nSa7iOl1F6gSX0FInIR8ByQrZQ6p0UFO4V44YUXfIvf/va3v/HrX//ad01E6N+/P//973955JFHWL16daPHLJpLXYOPAJs2beLMM8/0reiu4O677+bFF0+dXkSti6cGkZGRLFmyhOuvv54ffviBwYMHc8MNNzBq1Cji4+Nxu90cOXKEzZs3+9YgteXCtfoQpfxzvHSqIiLPK6XuNafwVblZpdREP/OIxvDVtMyfF1FE1On4XM866yzWrl2L1Wr1LaPvSCilfEv4hw4d6vMAC8Y01QkTJvhcPFQwZMiQWldrtyLN8vvRGF08XfXwVMLr9fLNN9/w9ddfs3v3bo4fP47D4fCtUUpOTmbkyJFMnjy51okbrUytunjaGwXwzTaarpRa1Mx8VndmozBgwABSU1OJjo72+W0Cw7PjhRdeWMUVR8XaiLfffrtN/R+NGjWKLVu2AHD8+HGfPynAt7oY4NZbb+V///sfAwYM8M1WaiNa5GH4o4unqx5qWoxadbFtNy1tJ8w34472luNUp3///oDhavrHH3/0nbfb7Vx11VVcddVVXHnllfTs2ZO0tDTmzZuHw+FoUxmnT5/uC//3v/+tcs1msxEfH098fLyePqrR1EGnaCkAiMh8wAZsBrxQc+c1EekKfFAtaYZS6hrzep21s2qOyEafjs914cKFzJgxAzDceLz//vu1Lqh75plneOSRRwAoLi6u1ddSa3H48GH69evn8xuzYMGCGn6cjhw5wowZM9i6dSuDBg1i9+7dbSYffrYUmqqLnUEPNS1G5+0+AhCRGj4XlFLfNjKPTt19BPD444/z1FNPAYZX1GuvvZaxY8disVg4duwYa9eu5eOPP+b48eNER0eTlZXVKgvw6uOll16qsqnI5MmTmTRpEuXl5WzevJklS5b4XAf/7Gc/45133mlL8XT3kaajULsu1jUt6XT7ACHADIx9FX4O/LwRacdg+K7PN7+DG4ivTmdWrFihZsyYoQIDA31TQCs+IqJGjBihHn/8cXX48OF2k/H5559XQUFBNeSr+AQGBqr77rtPFRcXt7VozdVjv3XxdNdDTbOpVW86U0thFfAllfZoVkrNbaWyVGd4ri6Xi9TUVHJycvB4PERFRfl2LOsIpKen89Zbb7FixQoyMjKwWCz06tWLSZMm8bOf/azNPWeatNmoe2fRQ02T6fTdRwuUUpe2UVn6ZdTUhTYKmo5CrbrYaRavAVYR+QzYyUnX2Y/Xn0Sj0Wg6F51iSqrJq8B64AIMVxedzg+9RqPRNERnMgq3AlsBp1JqKfDrdpVGo9FoOiCdySjYlbGiudw8brtlthqNRnOK0JmMwo8i8hAQIyL3AX47vBGR2SKyzvxc13oiajR1o/VQ0xZ0mtlHACJyCTAY2KuUWtCIdMlKqTQRsQHrlFKjG4ivZ31o6qLJLVSth5oWptPPPkIptRBY2IR0aWawnJPdTxpNm6L1UNMWdKbuo5bgDmB+ewuh6fRoPdS0GtooVEJEuorIimqfD8xr44CfAH+pI+1sEdkoIhsrTlX+iEiSiGwREafZ/JeGPiIyXkTWish3IvIPP9MMN9OsFpG3RCTAn3SV0t8vIt81In5vEckSkW9FZEkj0t0oIsvNdN39THORGf9bEckw99luKE2oiHxpplkgIsF+lmUTkf+aOvDXBuLW+G1F5B/mb/BCtfgN0tH0kFNAF5uqh2hdrEld/i/0p4oPmSRgNRDTjDyCgWhgBWD1M01XTN82wLvAcD/S2CqF3wTGNkLGIGAusLoRaZKBeU14nq838zdZD4T5Ee8K4HEz/Bgw08/8ZwGPmOEXgRR/f1tgFPAf89orjfkNOqIengq62BQ91LpY+0e3FPzjcSAB+NS01o3eckwpVaaUymtkmgylVJl56AY8fqRxVzp0AocbUeStGC9iY5kiIqvEmNXlD9MBi4gsE5EXRaSm/+16EJE+QKZSqtiP6PuBCt/dUUCOn8X0Abab4a3AWXVFrOW3HQ8sMcNLgQl+ltkQ7aKHZrpTQRcbq4egdbEG2ij4gVLqdqVUf6XUZPPTpjvHiMgIIF4p9YOf8S8VkZ0YfyB+KZ7Z1JyslFreSPGOAwOAKcA0U9aGSAAClVJTMbaWnNnIMq8APvMzbiowQUR2YXgYXeNnur1Ahbv1KRgvsb9EAYVmuKCRaeukvfUQOrQuNkUPQetiDbRR6OCISAzwEkbNyS+UUguUUsOAI8Alfia7AXivsfIppZxKqRKlVDnGzK5hfiQrAL41w8sxpgk3hhmAv1OKbwT+p5QaCnwB/MzPdP8D7CKyDKOWm9kI+QqAClexERhurk95OrIuNlEPQetiDbRR6MCIiBWYB/xWKZXhZ5qgSoeFgL+1yYHAnSKyCBgqIvc0lMAsL7zS4dkYTeSGWANU1OTOAA76KWPFjmQupZS/TW8BKjaUzgYi/UmklPIope4xa5AeYLG/MgJrgalmeBqwrhFpOyQdXRebqIegdbEG2ii0EeYMgqVACrDYnEXSELOAscCzZh+yP33TF1bMisBoGn/tj3xKqYeUUtOVUhcCu5RSL/qTDjhXRDaJyBrgqFJqvR9lbQUcIrIC4/4+9rMsMJr38xsR/z3gp2ZZ12MMkjaIOYtjhYgsB9YopY7WE7fKb4ux7WuZGHt4eJRSHcb5YhP1EDq+LjZaD82ytqJ1sWoe5qi0RqPRaDS6paDRaDSak2ijcBogIpNF5Ol2LP9CEbm4kWleMxfUJLWSTLeIyEERmVfpXKK5SGmNiExrjXI7O1oXa83/lNJFbRQ0zUYptUgp9UXlc2KsXq2PgUqpcyr3j/qRpjEsAM6vdu5h4PcYGy39rgXL0nQQtC42H20UTh/OFJGvzNpHjDngtExEVorIJyJiEZF+Yrgd+EZEHgUQkV+YC35WicioyhmKyBwReV9ElorI6+a5M8zBw/WV8rjJzCfZzPtj4CYRmWvG/abyS2bWJEeIyEKzZrlARBYA00XkITFcKSwXkZ5m/B0iMs/8vlZEvjAHFbvX9TCUUtnUdBo3HGOgrhgoEpGImik1LYDWxUqccrrozzJr/enYH2AysMQMX41RCxGMjYUAnsaoqdwK3GSeEyAOoxYjQAzwebV85wCPmeFXMFZH2jk5QeEb8/gm4BcYrgZ2ARaMWQ/LKsqqRebVlWRfZIa7AovN8DnAK2Y4E2P5/lnAZozKzHXA/Q08l2QquT4AVlYKzwN6tvdvd7p9tC6e+rqoWwqnD1vM761AP4zl9K+LMR3wKiAR+AijVvQucCHG8vkUjBfqU2pf7Vg9397Al2a+g4Eu1eJvU8a8ajcw1+xHfbqB5vhm8zuZk8v5N5rlAexThouFY8BupZTXDEdXzkRE3hVj6t4ZdZTjrRQ+bRaVdUC0Lp7CuqiNwulDSqXv/Rg+XX5USk0CPsGogbmVUvcDNwNPYizU2aBMtwnU7PesLd87gb+Y+e6jprdFL4AYPmTeV0r9DIjHmANeFxUvSFql8sZwcgFS5XnTlcNVylZKXW/ey9Y6ytkuIhNEJBSIUEoV1hFP0zy0Lp7CutipNtk5zXGLsQI0GLgSoyn9mIiMwVjqngpcKiJ3AyEYTdkTZp/oSoyVksuBp6rlO1iM5fXpSqm1IhIJvCQiPwCueuQJBxaYL2QhsKOhG1BKZZh9vmvMvG/0//arIsYuew8DfUXkE6XUlcCzwNsYz+aJpuataRCti5U41XRRL17T1ImIzMHob13a3rJoOjdaF9sO3X2k0Wg0Gh+6paDRaDQaH7qloNFoNBof2ihoNBqNxoc2ChqNRqPxoY2CRqPRaHxoo6DRaDQaH9ooaDQajcaHNgoajUaj8aGNgh+IyDgxNsNYLSLPtbc8ms6L1kVNa6ONgn8cAs5TSp0DdBGR4e0tkKbTonVR06poh3h+oJTKqHToxnDYpdG0OVoXNa2NNgqNQERGAPFKqR9quTYbmG0ejtbuQzR1UN29c9MyqUMXtR5qGkGtuqh9H/mJiMQAnwM/rVZbqy2u0s9VUwfNNgr+6qLWQ00D1KqLekzBD0TEirFl3m8bMggaTWuidVHT2mij4B+zMHZretbcYm9Cewuk6bRoXdS0Krr7qBXQzXZNPbTImIJfBWk91NRP5+4+EpHYasdR7SSKRqPRdFg6jVEAPqp2/Gq7SKHRaDQdmNN+SqqIXA5cAQwRkbfN01Ygov2k0mg0mo7JaW8UgOXAZiAd+DdGP5ob0DM3NBqNphqdZqBZRJIwWgxRmAMsSqknW6ksPcCnqQs90KzpKHTugWaMxT7ZwErgW/PjNyKSKCKbRaTMnCuu0bQ5Wg81rU1nUqpDSqn3m5E+F5gKfNZC8mg0TUHroaZV6UxGoYuIbAJ+ABSAUurn/iZWSpUBZSJt1vo/rWhMN4Z+xnWj9VDT2nQmo3BDa2ZezRGZxsR98Agnrr0Xb1aOfwmCAon87W2E3zqrdQU7TdF6qGkunWmguSfwKBAG3AjcrJR6rQn5rACmKaXK64mjB/hMlNNFyYdf4snK9p1z7z1I2dLvALD0TCRkxnm+axIcjP0nk7H17t7msrYRLeUldQVaDzXNo1Zd7EwthTeAu4GXlVIeEbkWaLRR0DQOCQok7IbLqpwr+fxrn1Gw9e5B5G9+0Q6SaTSa2uhMs48sSqk9lY4bde8iYhORpUAKsFhExrWodBqNH2g91LQ2namlsFxE/gUkisjzwJLGJFZKuYFprSKZRuMnWg+bTtnqjeQ/8TzewiLfOeUoQ7mNHjgJCkSCAn3XLHExxDz3O2yD+rS5rO1JpzEKSqmnzP1slwF7lVLb21smTcdAudy496eDy+VXfAkKwtq3J2LrNK/PaUHwOWPouuydKudyH36W0g+/BCDi/lsJv/3a9hCtQ9FptFpE7EBfjIHmMSIyRin1RjuLpekAlB84TOlni1HOk0bBtXEH7t37AbANG0DgyCG+axIcROhVF2Hrn9zWomo6GaVfrqD49Q/xZOYg9iCCJ40j4r6bCQgNabUyO41RABYDXwLH21sQjX/8ac6TOAqK/Yprjwzj0TmPA/Dgrx4hMz27gRQGCT3jePaFPxP16F1Vzuc//U+fUQiedGaDg+Eul4ucnByKiopwu90opbBarYSFhREXF0dwcLBf8mg6Jm2piwDewmLy5zxP6edLsA0bQNzbf6Xknc8pfuMjHEu/I+avDxM0dkTTbqYBOpNRcCilnmlvITRg7ZHoC1t6JdYZz1FQzKOzbvMrzz999B9fODM9m+SiM/xKl5a+1a94taGUYsuWLWzatIljx46RmJhI9+7dCQ0NBaCoqIh9+/aRn59PUlISF1xwAT179mxyeZr2oy110fn9NnLv/xOeY5lIWCjBU8bj+GolAXEx2FIG4962mxPX3EvYL35K5P23IEFBjb2demlxoyAiNwFvVjv9plLqlkpxBPgrxnqBEOB74Aal1JGWlqeyDCLyGbCTkyuaH2/F8gDweDwcOnSIgoIClFKEhITQq1cv7HZ7s/NWSpGamkppaalvxbDL5eLAgQMAWK1W+vXr54svIsTHx5OUlFQlD6fTidPppLy8HK/XS0BAAFarleDgYAIDA5u9wlgphXvPflwbtuM+cBhvdh7BF04Er0I5nOT86g8EhIZg6ZVE0KihBI4ZjgR0/IlxX3/9NevWrQPg6quvZtCgQTXieL1edu3aRV5eHjabra1F1NSC8nrxHD5O+ZEMrL2SCDp7NAQICqO7xhIfgzW5O5b4mDaXrfjd+eQ//g8kNIToZx4k6NwxVa6HXHEB3twC8n//d4r/81+cqzYS/+GLBIS1XHdSa7QUPgSWmuHLgReoOdNnBPAb4G0z/kLz+L5WkKeCBzEM0bGmJBaR54AxwGal1L0NxT9y5Ajbtm1j165dREREMHr0aEJCQti4cSOffPIJ/fr1Y/jw4QwYMACrtWk/g4gwYMCAKudOnDjBokWLAAgNDWXWrJorg/Pz81m1ahWHDx+mvLyclJQU4uPjsdvtWCwWvF4vJSUlZGRkkJOTQ3R0NOPHjyc8PLxKPtnZ2aSnp+N0OgkMDCQpKYmuXbtWieMtKCL71odxbd6FhNqJfuZBgiePq9InqjweXFt3U7pgKcVz9xLqdBF87tgmPZOWQDmdWHp0I+isUSCCJakrqrwcqfY7HTx4EACbzUb//v19591uN/Pnzyc/P993TkQ4ceIEV155ZbPla6wudnQq9M3hcOB0OvF4PIgIIoJSChHBZrNht9sJCwtr8vviOZFDwV9exbFkNQQEEH7TlQSOGkrw5PFIcBCqzInnaAalX3yD89vvsSQlEP3sw1gTu7TwHdeNbVBfwu+6HiwWPFnZlH6yiNKF3+DNKwDAPu1sLN3iCZoygcDxI7F0iUUFCMeOHaO8vM51jDUICwsjJqZ2o9fiRkEpVQqUAojIlUAm8Em1aMeBAow9DfaZ54poXfY31SGeiIwCwpRS54rIKyIyVim1ob40y5cvJzY2lunTp+P1enE6nRQVFTFw4EBSUlLIzs5m8+bNbNu2jQkTJpCcnOxLq5QiKyuLEydO4PV68Xq9hIeH0717d4LqaSoWFRWRmZlJ3759UUoRHBxMeno6Xbt2JTDw5FQ7h8NBXl4eeXl5xMXFER8fT2xsLMHBwQQEBODxeLDZbDgcDnJzc30vaoVs+/btY926dRw+fJhRo0aRkpLCjh07+Prrr0lISGDcuHEMGjQIi8WCJycf9w/GT2xJiKthEADEYiFwWH/K96fjycxGglu2Oewv5cezKHl3ASUf/A9rryTCb78Wb14hRf/5gKIX3yb0uhmEXn0JlrhoACZNmsQnn3yC2+1mwYIFTJkyhaioKGw2G1dddRVerxelFCtWrGD16tUcOXKk2UahKbrYFng8HkpKSqr4uCorK6O0tBSAwMBAwsLCqqSx2+3s3buXpUuXUlhYyODBg5kwYQJdunTBarUiIng8Hl833JEjR4iIiGDKlCmUlpb6dBKMis533xkLIkNDQ5k8eXKNsryLV1P66WIjznWXEnHvTTVvZHBfgs8ZQ9k5Y1AuNxIcWDNOKxI0ehhBo4dVOedYvArviVwA7BdNIvicMTXSdQ0Owuv1+o6PHz/OG28Y82iioqL45S9/WSW+xWKpU4ZWG1MQkaHAJOAppVT1uX4ngPcxau+/BfYAf2stWUwSKjnEA1CNcIg3npOtnaXABKDeF/Hqq6/2DTb++OOPLFu2DK/Xy4ABA5gxYwZKKV8tKCAggLy8PFJTU0lPT6esrAybzUb37t3p0aMHhYWF7Nq1i/Xr1+P1eomPjyc5OZnevXuzb98+0tPTKSwsJDw8nMjISM4++2wCAwMpLCzk4MGD7NixA7fbTXR0NH379iUhIYGZM2fi8XgoLi7G4XCwYcMGjh8/jlKKpKQkevXqRWRkJOecc46vhrZjxw6+/fZbSktLSUlJYfDgwdhsNrKyskhISOCCCy4gLy+PRYsWsXjxYsaPH8+ECRNIWDKXko++wvX9NrKuuAtrchKW+FijduZy483OpfzgESQijKCxI7Akda3v0bY45YePU/DXV3F8vZqg8WcQfsf1BMRG4S11QJCN8Lt+hjc7D8fS7yh88W3sF04i8oHbGDx4ML/61a/YunUraWlpvP7661gsFoKDg7FarTidTlwuF6GhoaSkpFQx/M2gUbqYn59fpUuwYjAcIDY2loiICJRSvjh2ux2n00lxcTFlZWUopSgvL8dmsxEUFITT6cTtdvu6wgIDAwkPD8dms3H8+HEKCgooLy8nMTERh8PB0aNHAYiJifFVNIqKinwty379+pGUlITb7fZ1Y3788cc+GUePHk1cXBzdunWjd+/e2O12AgICOHDgAEePHqWoqIi+ffsSERFBSkoKYLTcCgoKKCkpITc3F4vFQp8+feh2/aVYErvg+OpbXFt2ceKG32DtmUhATCRYLKgyJ96sHNypaQREhBE07gyCzxrVEr9ZswgIN8aosAQgITW7nb1eL/v27ePEiROUlpaSlJREcHAw06ZNw+v1EhISwokTJygrKyM/P5+AgAASExOJj4+vtbxW830kIi8DtwHJGMZnr3npHYyWw1fASxj7HCwC/q2UurtVhDHk6VX9nFLqkJ9pH8Voqi8SkWnAWdU36KnmiGz0E0880VyRT0vS09NZ8c03ZOfkYLNasdlsjD3zTEaOHElAtXGE/dt+4NUH/uhXvrP/+hh9U4xpo998/h2TYy/zK92KnM+ZPPMsDv6YSuqBA2Tn5hAREeH7owPjpSstLSUrKwtrgIXkXj0Zf/bZTRprmTNnTrMGaBrSRa2H/uHxeNiyZQsbNmwApSgrKyM6KoqzzzmHfv371/ht20oXz71kHLt37yYtLY2CggKio6MJCQkhMMCCEnCbLaecHENPhw0bxsCBAxt17xXUpYut0lIQkXDgZ8B8pdRRczOQijZRIXCGGS4zP14gqXo+LYwbuJJKO68B/u68VsDJPZ0jgPzqEZRSrwKvguGITEQYNGgQs2bN8v3ZFRQUkJGR4at9ZWVlYbVaGTFiBEFBQWRmZlJUVITFYiEsLAybzYaI4PV68Xg8KFN5i4qKsNls9OnTh7feeoucnByioqK4+eabiYgwxPR6veTl5eHxeCgvL+fQoUN4PB6GDRvGsWPH+OijjwCYPHkyEydO9L0ER44c4cSJE4AxaJ2bm0toaCgjR47E5XL5anAVLZ3//ve/vmb8zJkzCQ0N9fUHBwcHk5CQgIhw+PBhHnzwQT744AMmTpzIV4sWER8fz0MPPcRrr71Geno6f//735k6darvmf7+vgf9/HkgOTmZOXPmAHBw621+d0YmJyczYMAAnnzSUIVnn32WBx54wHd97dq1vhaa1+tl//79REdHc/XVV3Ps2DE2b96Mw+FgypQpDB48uEqft9vtpqioiPT0dHJycoiLi/P7fuqhXl2srocVz6QtcDgcFBQUUFpaiogQFBTk04UKKgyszWYjOjqaQ4cOsWbNGgoKChg5ciQDBgzwjW1lZ2dTWlqKw+EgMzOTwsJCunfvzvjx49m2bRsnTpxAKeVLY7VaKS8vJycnx/e+ZGdn4/V66d69O6NGjUIpxYIFC3jooYdITU3lrrvu4rHHHiM7O5v777+f995/n0mTJvHMM88wfvx4n9xtpYtZWVl8+umnAKxfv54zzzzTd/3DDz/k0KFDKKVwOBycOHGC3r17M3v2bLxeLzt37mTnzp2UlpYSFxfH8OHDCQkJ4fDhw/z4448UFxfTs2dPhg4dWmMssjKt1X30cyAc+CeA6cmxYuwA03fLS8DNGE7qtgCtPRNoPvB3oCkrmdcCt2MMik8D3moogdVq5fDhw2zatIkePXoQHh5OQEAAISEhOBwOXC4XGzZsoKysjLy8PGbOnNno7gWlFNHR0eTk5OBwOMjKyiIsLIyAgABf073iD3v79u1kZGSwfft2br75ZsaNG0dqaiqbN2/G6/WSkJDga5rb7XbKyspwOBx8//33BAUFER0dzfDhw4mNja0iw7nnnovD4cBqtTJw4MAqYxcATqeTZ599lmeeeQar1cqvfvUrYmNjfUapV69ePPzww8yfP59p06YxY8YM/v73v1eZNdXa/PSnP6WsrIyvvvqK9957j127dpGcnEx0dDR2ux2bzUZRUZGva6x79+5MnjyZCy64gKlTp1JSUuLr4ktPT8flcmG1WklMTKRfv3706NGDoUOHEhLSIjNEGq2LbYXdbm/0rLrhw4fTu3dvCgoKfO8FGAPzISEhhISE+PSzwtAEBwdz1lln1Zln7969az2/a9cu7rrrLlauXMn555/Pvffei4jw+eefA3DFFVdwySWX8O9//5sJEyZw+eWX8+KLL1aZrdfa/P73vyc8PJxvvvmGW265hQsvvJABAwYQERFBYGAg3bt3Jz8/n61bt7Jo0SKSk5O5/PLL8Xg8REdHc9lll9GlS5cqhnjAgAFMnToVt9vt657OycmhW7dutcrQKkZBKfVPTINQx3UvcI/5aSuavPOaUqpi+8NVwFal1PcNpXnooYfIzMwkLy/PN0PH6/X6ZlIEBQVx+eWXExUVVeOP1l9EhOuuu46MjAzS0tL44Ycf2LJlC0FBQVUG6ir6eIcMGULPnj2x2+1ceOGFXHjhhb7arMPh8E1JDQ4OJjw8nB49ejB+/Hhfja82zjnnnHpldLvdjBgxgrlz5zJp0iTS0tL47rvv2LdvH4WFhdhsNhITE3n66aexWq04HI5GbcjTEthsNq6++mrsdjvffvstu3fvJjU1lcjISGw2G16vl9zcXOx2OzNnzuSCCy5g2jTD/ZCIkJuby5IlSzh69Chjxozh7LPPZt26daxbt479+/czdepUEhISWkTWpuhiRycsLKzGIHRr0K9fP/71r39RXFzM1q1b2bhxI2lpaeTl5eF0OrHZbMTFxXHeeedxxx13MGXKlBb73fylR48ePPDAAwwbNoxNmzaxY8cO9uzZ43unK97XhIQEHnnkES666CL69PHPN5PNZqNv374NxjvtF6+JyDsY6xKau/Nao6b+Wa1WkpKSWr2WISJ069atTqvfEDabrc6paS1BWFgYM2fOZPny5UyZMoVdu3YBEB0dTUJCAnl5eWRmZgIQHx/Pv/71ryrTO9uCnTt3MnXqVLKysujXrx8ffvghZ5xxRhVDqJRi5cqVzJs3j3fffZdevXrRrVs3li5dyu7du0lISGDkyJFERkayc+dOwsLCGD16NMeOHeP999+nZ8+enH/++XTv3vx9Ik6HaajtgdVq5Y033uAf//gHAL/85S959NFH6du3L0FBQeTn57N69Wr+9Kc/8dJLLzFu3DjfOpS24rXXXuOuu+7C7XYza9Ys3n777RqGyeVysXDhQpYtW8Zbb73FU0891SLrnio47Y0C8Lv2FqCzs2DBAq644go8Hg9jxozhn//8J2PHjvX96W7cuJEbbriBPXv28NOf/pR9+/a11Ewdv9i2bRtZWVkAjB8/npEjR9aIIyJMmDCBjIwMCgsLiYyM5MSJE/To0YOUlBQGDBhQa2tKKcXBgwfJzMwkOzu7RYyCpmm89957/O1vxiTHv//979x3X9VlUTExMfTu3ZtrrrmG3bt3t8vq80OHDuF2uwEYNmwYXbrUXCNRMXMrJCSEwMBAXC6XNgqNoWKGkYjco5R60QwLxljGi+0pW2egvLycO++8E4/HQ0JCAkuXLiUyMrJKnDFjxrBixQoeffRRwsLC2rzJft111yEivPzyy3zwwQekp6czdepUunfvTmBgIMXFxaSmpvLZZ5+Rm5vLxRdfzPTp0/2a9SEi9OnTx+8mvqb1qJhAAdTosv3zn//M3Llz8Xq9uN1u0tLSCAkJ4dNPP2X69OltJuOTTz7J+eefz6effsrSpUv5+uuv6dmzJyEhIYgIDoeDjIwMioqKSElJ4bLLLvNNLmkpTnujUInLMY2AUkqJiO9Y03ps376dY8eMReSXXXZZDYNQMYspPj6e//znP1Vmq9gjw6r4kakPe+TJPumEnnF++zRK6BnnG5u57rrrKCsr89Xsi4uL8Xg8xMTEMHz4cO655x569uxZY+qs5tTg1ltv5cMPP2T9+vU88MAD2O12fvKTnxAaGsojjzzCI488gtfr5bvvvmPixImUlpZy+PBhoG11ceLEiUycOLHR99dSdKY9mr8BZiulUkWkP/AfpdTkVipL741rsnLlSiZNmgTAgw8+yF/+8hfftbfeeovbbrutyvJ8u93OP/7xD2bPPm33nm+RPZr9KkjrYQ08Hg+fffYZn3zyCRs2bCAjI4OoqCgiIyNxOBzk5OTQpUsXhg4dyoQJE5g9ezbR0dHtLXZrUasudiajMBD4MxAHZAGPKaX21p+qyWXpl9EkKyuLxMREPB4PI0eOZNOmTb6WQGZmJgsXLvS5QpgzZw65ublcd911vPvuu+0pdmuijUIHomLdhMvlwmazERYW1mwHkKcQndsoNAcRuQh4DshWStU/BxP9Mlbn1ltv9flhefnll7njjjtqvHhKKfr06UNaWpo2CvUlboQuaj3UNEDbrWjuSIjI80qpe8153VXeEKWUvx136zA2Sl/W0vJ1Bl544QVycnKYP38+d911F6+99ho/+clP6NmzJw6Hw+cbKi0tDYAzzjijXeXt4Ghd1LQqnaKlYM42mq6UWtTMfFbrlkLTUErx1Vdf8eabb7J69WoyMjJ814KCghg6dChTp07lmmuuYdSo9ndC1oq0SN+EP7qo9VDTAJ27+0hEPldKXdbMPOp8Eas7Iussz7WplJaW+nw4RUVFdaYZPa1qFLQeahpBpzcK8wEbsBnDAV+NnddEpCvwQbWkGUqpa8zruqWgaS5+GYWW0EWth5oG6JxjCpX4e0MRlFIZwOTWF0WjqR+ti5r2otO02TE2IokAegK9zI9fiMgY07PrMBFZKiLBrSSjRlMvWhc1rU1n6j5aBXxJpT2alVJzW6ks3WzX1IVep6DpKHT67qM8pdSf21sIjUaj6ch0JqNgFZHPgJ2cdJ3d2hv7aDQazSlFZxpTeBVYD1wArARO+c1JNBqNpqXpTEbhVmAr4FRKLQV+3a7SaDQaTQekMxkFu7miucIlZ6fxeqXRaDT+0pmMwo8i8hAQIyL3Abv8TSgis0Vknfm5rvVE1GjqRuuhpi3oNFNSAUTkEmAwsFcptaAR6ZKVUmkiYgPWKaVGNxBfTwXU1EWTW6haDzUtTKefkopSaiGwsAnp0sxgOSe7nzSaNkXroaYt6EzdRy3BHcD82i6YTfuNIrKx4lTlj4gkicgWEXGaNT1p6CMi40VkrYh8JyL/8DPNcDPNahF5S0QC/ElXKf39IvJdI+L3FpEsEflWRJY0It2NIrLcTNfdzzQXmfG/FZEMc0vVhtKEisiXZpoF5gpgf8qyich/RWSFiPy1gbg1flsR+Yf5G7xQLX5L0KZ6yCmgi03VQ7Qu1qRij1z9UQBdgRXVPh+Y18YB/wMsTcw7GIg287Q2Qp5gM/wuMNyPNLZK4TeBsY2QMQiYC6xuRJpkYF4jn0US8Hozf6v1QJgf8a4AHjfDjwEz/cx/FvCIGX4RSPH3twVGYWz3CvBKY36DjqiHp4IuNkUPtS7W/tEthUoopTKUUpOrfa4RkSTg/4AblVKeJuZdppTKa4I8ZeahG2iwbKWUu9KhEzjciCJvxXgRG8sUEVklxgC+P0wHLCKyTEReFBFLYwoTkT5AplKq2I/o+4FQMxwF5PhZTB9guxneCpxVV8RaftvxwBIzvBSY4GeZFfl1KD2sJFNH18XG6iFoXayBNgr+8TiQAHxqNuHsbVm4iIwA4pVSP/gZ/1IR2Ykhs1+KZzY1JyulljdSvOPAAGAKMM2UtSESgECl1FSgFJjZyDKvAD7zM24qMEFEdgFjgDV+ptsLTDLDUzBeYn+JAgrNcEEj09ZHu+ohdGhdbIoegtbFGmij4AdKqduVUv0r1docbVW2iMQAL2HUnPxCKbVAKTUMOAJc4meyG4D3GiufUsqplCpRSpVjDOIP8yNZAfCtGV6OMSOsMcwA/J09diPwP6XUUOAL4Gd+pvsfYBeRZRi13MxGyFeA4ZEX8zu/EWnrpD31EDq2LjZRD0HrYg20UejAiIgVmAf8Vhn+9f1JE1TpsBDw949jIHCniCwChorIPX6WF17p8GyMJnJDrAEqanJnAAf9lLFi8xmXUsrfprcAuWY4G4j0J5FSyqOUusesQXqAxf7KCKwFpprhaRj7Kp/SdHRdbKIegtbFGmij0EaYMwiWYmy6vlhExvmRbBYwFnjW7C7wp2/6wopZERhN46/9kU8p9ZBSarpS6kJgl1LqRX/SAeeKyCYRWQMcVUqt96OsrYBDRFZg3N/HfpYFRvO+1pk3dfAe8FOzrOsxBkkbxJzFsUJElgNrlFJH64lb5bfF2OGvTAx37R6lVIfxs9VEPYSOr4uN1kOzrK1oXayahzkqrdFoNBqNbiloNBqN5iTaKGg0Go3GhzYKpwEiMllEnm7H8i8UkYsbmeY1c5VlUivJdIuIHBSReZXOJYqxcnWNiExrjXI7O1oXa83/lNJFbRQ0zUYptUgp9UXlc2K4NKiPgUqpcyoPmvmRpjEsAM6vdu5h4PcYGy39rgXL0nQQtC42H20UTh/OFJGvzNpHjDkLYZmIrBSRT0TEIiL9xPBF842IPAogIr8wV4GuEpFRlTMUkTki8r6ILBWR181zZ5gzStZXyuMmM59kM++PgZtEZK4Z95vKL5lZkxwhIgvNmuUCEVkATBeRh8Twr7NcRHqa8XeIyDzz+1oR+cKcadK9roehlMqmptO44RizN4qBIhGJqJlS0wJoXazEKaeL/vje0J+O/QEmA0vM8NUYtRDB2FgI4GmMmsqtwE3mOQHiMGoxAsQAn1fLdw7wmBl+BWPJvJ2Ts9a+MY9vAn6B4X9mF2DBmAq3rKKsWmReXUn2RWa4K7DYDJ8DvGKGMzF8upwFbMaozFwH3N/Ac0mmkj8cYGWl8DygZ3v/dqfbR+viqa+LuqVw+rDF/N4K9MPwsfK6GHPErwISgY8wakXvAhdi+FRJwXihPqX2JfDV8+0NfGnmOxjoUi3+NmUstnEDc81+1KcbaI5vNr+TOenjZaNZHsA+ZfjdOQbsVkp5zXB05UxE5F0x5nOfUUc53krhFltprKmB1sVTWBe1UTh9SKn0vR/D0dePSqlJwCcYNTC3Uup+4GbgSYzVmxuU6TaBmv2eteV7J/AXM9991HTB6wUQw7HY+0qpnwHxGAuD6qLiBUmrVN4YTq5KrbyYpnK4StlKqevNe9laRznbRWSCiIQCEUqpwjriaZqH1sVTWBc71SY7pzluMdwCBANXYjSlHxORMRj+T1KBS0XkbiAEoyl7wuwTXYmxfH458FS1fAeL4XMlXSm1VkQigZdE5AfAVY884cAC84UsBHY0dANKqQyzz3eNmfeN/t9+VcTYZe9hoK+IfKKUuhJ4Fngb49k80dS8NQ2idbESp5ou6hXNmjoRkTkY/a1L21sWTedG62LbobuPNBqNRuNDtxQ0Go1G40O3FDQajUbjQxsFjUaj0fjQRkGj0Wg0PrRR0Gg0Go0PbRQ0Go1G40MbBY1Go9H40EahEYjIfSKyur3l0Gi0LmpaC20U/EREgoAz2lsOjUbroqY10UbBf24F5ra3EBoNWhc1rYg2Cn4gIjZgslJqeXvLouncaF3UtDbaS6p/3AC8V18EEZkNzDYPR2v3IZo6qO7eubHUq4taDzWNoFZd1L6P/EBE/oLRh6uAccDjSqkX64mv9HPV1EGzjEJjdFHroaYBtFFoCURktVLqnAbi6JdRUxfNbSmczKgBXdR6qGkAbRTaCv0yauqhxYxCgwVpPdTUT6262GnGFEQkVimVU+k4SimV344iaerBuXEHxW98hNdR1mBcQbD27k7kg7OR4KA2kE6jOX3pNEYBY6Pw8yodvwr8tJ1k0TSA44tvcCxa6X+CVRsI/dlMbH16tp5QGk0n4LTvPhKRy4ErMDYC/9o8bQVilFIXtlKZutneTJS7HNeOvahSh+9c/tP/pPzHgwBE3HsTgaOHGRdEsPZMxNqjW3uI2lh095Gmo9Bpu4+WA5uBdODfGA/CDWS0p1Ca+hGblaBRQ6ucC4gM94Vtg/sSfM6YthZLozntOe2NglKqACgQkZcxWgxRnLSQT7aXXBqNRtMROe2NQiU+B/4ObG9sQhEZBzwHeIENSqn7WlY0jcY/tC62PSWfLMK5bgt4vA1HDgggMGUwYTdc1upytRadySgcUkq939S0wHlKqTIReVdEhiuldrSkcBqNn2hdbGNcm3biXLcVvKZRUApPRrbvekBCHBJgdj4EBKBcLm0UThG6iMgm4AeM1aAopX7uT0KlVOXxBzfgaXnxNJqG0brY9kT/6bdVjpXXy9F+Jycydl38JgER4dWTnbJ0JqNwQ3MzEJERQLxS6ocWkEejaTJaFzWtRWfykqqAR4A/Akcwpqj6jYjEAC9huC2u7fpsEdkoIhubK6hGUx/16aLWQ01z6UxG4Q3gH0CiUsoDXOtvQhGxAvOA31ZrvvtQSr2qlBqjlNLzJDWtRkO6qPVQ01w6k1GwKKX2VDpuzL3PAsYCz4rIChGZ0LKiafzB1j/ZCAQEYO28K5e1LmpaldN+RXMFIvJ7IAmYDCwGMpVSf2qlsvRK0hbEW1RC2aoNuLbtwf1DKmKzYeneFVvfngSNPwPbwD7tLWJj0CuaT3GqDzQnbv3fqTrQ3GlXNAOglHpKRIYDy4C9SqlGr1fQtC3K6aTguTcpefszVJkTRLAN7IPYg3Dt+pGSE7kABI4aSvTTv8E26JQyDppTGasFys2JX9bT62+0M7UU7MB0IBrTQiql3milsnQNrZkodznZNz2Ac+0WAMJmX0v47KuxxEQZ15XCuXojub/5E97sPCQynG4r3qviCqODolsKpzDK66X80DFKP/wC59YfsPXrhf2CiViTk7AkJSABp1SPfOduKWB0GX0JHG9vQTQNU/r5Ep9BCJ4ynqiHb69yXUQIPncs8R+8QMEf/4mE2MF8If8050kcBcV+lWOPDOPROY8D8OCvHiEzPbuBFJDQM45nX/gznrwCvHkF4Ocfb0BYKJaEOL/iajoW7oNHKPr3e5R9vRpvfiGWpAQsiQmUHzhC3u/+D8/h4wTERmO/aBIRd12PpWs80DRd9FcP4aQutiSdySg4lFLPtLcQnR2Xy8WmTZvIyMjA6/USFRXFqFGjiI6OrhLP+f1WXzhkxtQ687P16UHc61V/VkdBMY/Ous0vef700X984cz0bJKLzmgwTVq6IZvneBaerNwqRiHnzt+Dyw1A9F8fJsBs2SBCQKhdG4VTkLJVG8i+7VFwubEN6kvc23/DNrQ/Iicr2s6tP5D38F8pmfc5ro07SPjydaBpuuivHsJJXWxJWtwoiEgQ8B9gBmABvgFuUEoVVov3AHCfGed14HdKKT+cizSZN0XkM2AnJ1c0P96K5WkqoZRi165dfP3115SVlTF16lRiYmL44osvWLt2LRMmTODcc88lMDDQSBBgOZnWU/+iXdfu/YjNiq1fr9a8hRoEDukPQ6qeE6sFZRqFoAmjsCZ2aVOZNC1PwZ//ddLQP/MAgcMG1IgTdMYQEha8ivvHg75WwqlKa7QULsRYPfwnDD8t/wZuBp6viCAi5wDPAo+ZceZhOKr7oBXkqeBB4K/AsaYkFpHngDHAZqXUvS0pWFtTWlrK0aNH8VT6s83JyaGgoACAqKgoYmJifNcsFgvdu3fHbrc3qbzMzEy++uorDh06REREBEOHDqWsrIxjx44xdOhQ9u7dy3fffcfWrVuZOnUqKSkp2C84h9KPvjTk/d8yQi6/oNb+2pL/fkHeI38FIGHRm9gG9G6SjM1FKYVr226CJ4/Dk1NAQKid8gPpWBJiEYul4Qwawemki6cCqtLufwHRkXXHc7nxZOUggTYscdF1xmsLPJnZuFPTUO5yAiLDCRzaHwkK9CttaxiFA4ALY9XwIfNcUbU4483vz4D9GEbhYlrXKOxvqkM8ERkFhCmlzhWRV0RkrFJqQ31pFi5ciNPpJDk5mbCwMKxWq6+56fF4cDgcFBcXExwcTJ8+fYiKimq0XEopioqKqvy5N0RgYCChoaH079+/yvlFixaxYYNxSxMnTmTQoEE10qanp7Nv3z7y8vIIDw8nISGBoKAgAgICEBG8Xi9Op5OCggKUUnTt2pWkpCTWrl2L3W7n/PPPJzExkfLyct+zUEqRnJxMYWEhP/74Iz/88ANhYWH0PW8CIT/9CaUffolz5QZyfzmH8DuuxTaoL1gteI5mUvLpYopeeRcAW8pgrL17NPoZNhf33gOULlhG6f+W4c0vJGTGVEJmnk/Z6o3k3PMHJCiQkJ9Mxn7JeQSOHFKly6EpNEUXlVK43W6cTifl5eV4TcduIkJAQAA2m43AwEBsNluzZDsVKCkp4ciRIxw+fJjjx43hxa5du+J0Ojlx4gR2u53ExES6d+9OYmIiwcHBhN97I3n3G7PX8596iZi/PVxjCmp5+jFy7vw97t37sXTvSreVrflXVjtKKZzrtlL85sc4v9tEyJUXEjhyCIXPvYHnWCah18wg9PpLsXarv/XaGkYhFVgEvIzh3ncVMLdanIraegoQZoZjW0GWyiRUcogHoPx1iIdhxJaY4aXABKDeF9Hj8eD1enG5XBQXF+P1ehERRASLxeIzCm63m7KyMkpKSsjJycHhcBAQEEBERITPkIgISik8Hg8ulwuPx4Pdbic2NpawsDDfS+7xeDh+/DgLFy7E6/USFhbGlVdeSUhIiO/PyGKxkJ+fz969eyksLCQoKIgePXqQnJxMdHQ0IkJsbCxHjx7F6XTidDoJDg4mKSmJwsJC8vPzyc3NJSIiwidXRflKKZRSlJeXU15ejtvtxuPxcPToUbKzsykoKKBfv37ExsZisVgoKSnhxIkTvj+qimdW8ayi//wAgSmDKXp5Ho7FK3Esrrk9p4TaCf35FUT8+ibE1rZDZMUffEHZt+uQoEAiHrjN2PnNZkMEQq66EPul0/Dm5OFcvYniuZ8SuGMv4Tde0dxiG6WLb7/9NocPHyY2Npazzz6buLg47HY7FouF8vJyioqKOHLkCOnp6YSGhjJp0iROnDhBVlYWTqeTxMREoqOjsVgsPl2s+K0KCwtxuVyEh4eTlJTUpJtxOByUlpZWObd9+3ZKSkoAGDJkCJGRJ2vnFouF0NBQ1q5dS3p6Ok6nk5SUFLp06UJgYCABZmuy4t4yMjIoLy+nb9++WK1WNm3aRFlZGYmJiUyePBmr1UpAQADFxcUkJydTXl5OdnY2O3fuJC0tjXHjxhF62QVY4mMp/MeblC39juPjryJw5BAsXWJRHi+eI8dx7dgLItgvnEjkw3c06Vk0FVVeTulnSyh+82PKDx/DfuEkIn59EwGx0aAUYTdegfvHg7g276To3+9hP/8cwn5xdY1NrCpo8SmpInIH8ArwKHAQeB94wjwGeBr4CzAfuAgowRhXmKuUarWnKSI1OpyVUodqi1tL2kcxmuqLRGQacJZS6slqcWYDs83D0U888URzRdZUoBTR+SUc2/Mju9MPERUUzLkDh1AeE0l2XDieat0z+7f9wKsP/NGvrGf/9TH6phgDA998/h2TYy9rMM2KnM+ZctnZhBc5iCwsxepwscdRyKGCXPJKSnC6XVg9XoIsVmLCwugVFkW/iGhUkI2b3ny5WU2FhnRR62HrYne4kMMZbNm9i+LSUiKD7YzvNxBrfCw5sWG4q1VMmqKL/uohGLo4ecZ4Cn48QHZONoWiiImLI8RqI6bUhQDllgDyQoNwOp14i4qJCLARn9yTx5//e5tNSa0YLHYAFZ1x0YC5oS655vdTwBxgFIZzr7daQZbKuIEradrOawVAhBmOAPKrR1BKvQq8Csb8cBGhT58+XH311b7B04KCAo4ePYpSirKyMg4eNPYbHjVqFF26dKGwsBC3201oaCjBwcG+mpnL5cLr9fpaC263G4vFQkxMDLt37+bQoUMUFBSQlJREWFiYr0sHqnZVWSwW+vTpQ0REBEuWLOHgwYPEx8czZswYXw0yICDAV8tKT08nJyeHbt26MX78eA4ePMjBgwcpKiqiV69ehIeHU1RUhNfrpby8nLCwMNxuNyUlJdhsNvr27UtAQACvvfYapaWlREVFccMNN/haJJUpLS1l//79vtZE165dq1xfvHgxv37h/3zHh6NDWP7B8lq7Y35/34N+/qyQnJzMnDlzADi49baaHZ31pCkrK+OPf/wjr7zyCjk5OUydOpXpV8ykW7duBAYGcvz4cRYvXsw7X3xCTEwMDz7ov1z1UK8uVtfDlJQUjh07RnR0NAMGDCAyMpLAwEBKS0spLi6mrKyM3NxcMjMz6dKlCxMnTsRms5Gbm0tpaSkigs1mIzIykqCgIMrKyigoKMDtdqOUIiQkhOjo6Cq1+UbdTEEBWVlZlJaWEhoaSlhYmK9VUlJSgsfjMbpFnE5KS0ux2Wz079+fXbt2ceLECZRS9O3bl+DgYKxWK+Xl5eTm5vrS5eTkoJQiMTGRPn36cODAAUpKSvB6vYSEhJCdnU1RkfGjx8XFERUVhcvlwuVyYbfbGThwIGFhYT55165dy6WXXkp2djaXXnop7y5YwJKMQ8yfP5/x48fXuL+m6KK/eliR5vdPPMH8+fPZvHkzOTk5JCcnExQUREhIiK/3INK8p1ybFWdQEL3HjKwzz9YwCu8A52H84Vsx9kj+P6XUkYoIIhIKvIvhduIA8HOl1LpWkKUy82nizmvAWuB24ENgGn4YsIsvvpi8vDzWrFlDWFgYNpvN90dd0T0SGhpKYGAgIkJYWFgV5atMSEhIneWkpKSQkpJS5VxFNw7gMyyVOXLkCCdOnMDpdBIaGkpcXBxxcXEEBgZy7NgxMjIyfH3QO3fuZPfu3djtduLi4nx/DBXdPBVjIfn5+TgcDrxeL6WlpZSWlpKfn895553HT37yE9LS0sjPz2fVqlWEh4djs9l83U/l5eUUFxdTWFhIZGQk5eXlVeR1Op3cc889AAwbNoydO3eyYsUK5s2bxw03NNsjepO54447mDvX6Bn9xz/+wb331hzzvffee1mzZg3z589v8kB9NRqli5dffrkvXNHVV6EfFWMKtelIeHjtiwDDw8OJj2+52TWRkZF1GpS4uLqn744dO7bOa8nJyXVeO+OMM6ocr1q1itTUVESEAQMGMGLEiDrTfvTRR9xwww243W7eeustbrzxRl566SXuueceJk+ezBtvvMF1111XZ/rWwuv18uGHH/LVV1/hcDh4+eWXmTJlCpGRkdhsNpRSlJaWkpaWxmeffcbx48fr/U9pcaOglHIA1zQQpwRoa58ETd55TSm1WUTKRGQVsFUp9X1DacaMaT8nlbW95JXp3r07d999N4WFhWRnZ5OTk+ObjVRhuAICAkhISODaa6/1/RFUzEJqDC6Xi+DgYPr16+c7l5aWxqpVqwDo1q0bkyZNqiJ7dYV97rnnSE1NxWKx8PHHH3P//ffz5Zdf8tvf/pYZM2Y0aZC+Jdi6dasvPG3atDrjjRkzBqvV2iJ/pk3RxQoqxrM0Jzn33HM599xz642jlOLZZ5/l4YcfBmDgwIGsXbuWtWvXAvhaINdffz179uxhzpw5vveoLQgMDOS///0vAMXFxeTl5VFaWkpmZiZer5eAgAAsFguJiYk88sgjREdH1yvfab94TUTewViX0OSd18y4p93Uv4iICCIiIhqO2AwCAwPp27dvlXNOp9MXDg0NZeDAgXWmP3r0KE8//TQAv/jFLxg4cCB/+ctfWLRoEVlZWfz+97/nxRdfbB3hG+Avf/kLs2bNoqioiDvvvJMXX3yRESNG+Ayy1+vl+++/55e//CWbN2+mX79+pKamNrvc01EXOzL33HMP//znPwkNDWXu3LmMHj26yvWHH36YDRs2cNNNN/HUU0+xf/9+3n333XaRtb4eB3857Y0C8Lv2FkBTlR49ejBx4kQ8Hg89etQ/jfSBBx6gpKQEu93O448baw2HDRvGzTffzOuvv87LL7/MzTffzKhRo9pC9CpMnz6d1NRU3nrrLb7++msuvPBCnE4ncXFxeL1esrOzCQsLY8SIEbzwwgvcdNNNbS6jpvnMnj2b8847j9GjRxMZGcmXX37J5s2bfTPnYmNjGT16NGvXrmXfvn31dl+dCpz2RqFihpGI3KOUetEMC3A30D5VzE5OdHQ0U6ZMaTDeypUref99o8cvODi4yvhBbq4xX8Hr9XLXXXexZs2aNm2yV5CQkMBDDz3EQw89BOAbFAUICgpq9roETfszYsQIBg4cyB/+8Aeee+45ysrKao0XEhLCH//4R664otnTjtuV094oVOJyTCOglFIi4jvWdDzKy8t9g8sRERE1jEhUVBQ2m40NGzawfv16Xn/9dW67zfAxY48Mq+LTqD7skSeb2gk94/zyJZPQs+4BUBEhODjYr7I1pwZer5errrqKhQsXAjBr1izuu+8+hgwZgtfrZc2aNfzmN79h79693HfffSQmJvLTn/4UaJou+quHFXFbnIqZCKf7B8MHU38z3B9Y0YplKU3zeOWVVxTG2I/629/+VmucgoIC1aVLFwWo2NhYVVBQ0MZSNom21Pm2uqfTmrlz5/p08fbbb1der7dGnIyMDNWnTx8FqBdeeKEdpGwStepNZ2op3AH8RUTigCyMaX2aDkpBQQEiwrhx43wthupERETwr3/9i2uuuQar1UpBQUGrD5xrOh/vvPMOYLQC//CHP9TaJZiQkMDOnTvJzMw85ccUOs0mO83BXCV6i3n4glLqvQbiK/1cm4/L5cJms7F9+3bmzp3L2rVrOXz4MCUlJURGRjJ06FAuv/xyrr76akJDQ9tlTKEJNHmQQeth+9CjRw+OHDlCUlISR44cqXJt9erVVc5ZrVbGjRvX4ASKDkKtunjaGwUReV4pda85r7vKzSqlJvqZR7JSKk1EbMA6pdToBuLrl7EFcLlc3H333fznP0afrIj4Vq8eOHDA5zNn8ODBLF68+JR+Ef1KqPWwXejduzdpaWl06dKFjIwMX0shPT2dXr1qumvv1q0bx441yRlzW1OrLp4SVavmYBoEAf6olJpY+dOIPNLMYLn50bQBs2fP9hmEG264gaNHj5KamsqOHTvIzs7md78zZhvv3r27zi6m0wmth+1DxSrnrKwstmzZ4jvfs2dPtm7dyoIFC5g/fz6zZs0C8HlfPVU57Y0CGKNtGGMKzeUODHcZNRCR2SKyUUQ2tkA5nZ5Nmzb5XEhMnz6duXPn0q1bN991u93OU089xaOPPorVaq2xQO40R+thG3LHHSf/Oh5//HHcbrfvOCUlhRkzZnDppZcyYEDNzXdORU777qMKRGQ+YAM2YzrtU9V2XhORrtTc0yFDKXWNiIzDWAh3mVKq3g0MdLO9+Tz22GP86U+GD/svv/ySiy66qM64yvTjc4rQoKBaDzse9957Ly+88AJgGIKbb76ZCRMmEBAQQFpaGsuXL+eDDz4gLy+PLl26kJmZ2c4S+0XnHFOoQEQmVT+nlPrWz7RJwH+BS5VSuX7E1y9jM7n++ut57z1jHHX//v306XPSVdZ3333HunXrKqZdIiIMGTKkXsPRgWjOmILWw3ZCKcWbb77JnDlzOHz4cK1xQkJCuOqqq3jkkUdq3aSqA9JmrrM7KhuAqVR1ne0vjwMJwKdmjfQiZTj+07QSlReAVd+E5cYbb2T//v010mzfvp3hw4e3umztiNbDdkJEuOWWW7jxxhvZsGGDz82FUorY2FiGDh3KhAkTWsoTbrvSmYzCYuBLmrBHs1JKr2loYyq7OP7iiy8YNmyY73jdunXs2bMHj8fDxo0b+e1vfwucdH1xuqL1sP2xWCyMHz++1r0TThc6U/fRAqXUpW1Ulm62N5OsrCz69+9PYWEhSUlJbNiwocpAcwVLly7l/PPPB2DFihVV3HB3UNps8EProaYBOn33kVVEPgN2ctJ19uP1J9G0F126dOHVV1/l2muv5ejRo4waNYrbbruNiy66iPDwcDIyMvjuu++YN2+eL019m7JoNBr/6EwthcuAQRiO8R4DgpVSC1upLF1DayGWLVvGPffcw+7du+uMM2HCBB588EEuu+yythOs6eiWgqaj0OlnH/0P+CfwqFJqoogsVUrVvV1W88rSL2MLopRi/fr1rFmzhiNHjlBWVkZkZCQDBgxg8uTJ9O7du71FbAzaKGg6Cp3eKCxVSk0TkeVKqfNEZJlSamorlaVfRk1daKOg6Sh0TjcXlfhRRB4CYkTkPmBXYzMQkfki8nTLi6bRNA6ti5rWotMYBaXUXRiG4F1gv1LqV41JLyIjgFN/ErLmlEfroqY16UyzjzAHlps6uPwr4GVgTMtJpNE0Ca2LmlajUxmFpiIig4ATQH49cWYDsyufamWxNJ2QhnRR66GmuXSa7iN/EJGuIrKi2ucD4H7g+frSKqVeVUqNMT81XkQRSRSRzSJSJiJ+GWMRGScia0RktYg852eaYWaaVSLypjTSU5yI3CciqxsRP1lEMs1n9XUj0v1cRJaZ6ZL8THNhpd/luDnNuKE0ISLyhZlmvogE+VmWVUQ+EJFvROTZBuLW+G1F5DnzN6hXb+rJs0m62Bp6aKbr0LrYVD0002pdrExd+3TqT5W9bhcDi4Dvgf3ApCbkEQxEAysAq59pumKspwBjLGS4H2lslcJvAmMbIWMQMBdY3Yg0ycC8Rj6LJOD1Zv4m64EwP+JdATxuhh8DZvqZ/yzgETP8IpDi728LjAL+Y157pTG/QWvrYlP08FTQxaboodbF2j+6peAHSqnpSqkLgQeB95Wf3lWr5VGmlMprZJoMpVSZeegG6nWVbKZxVzp0ArW7dKydWzFexMYyxayJ3Odn/OmAxaydvSgilsYUJiJ9gEylVLEf0fcDoWY4Csjxs5g+wHYzvBU4q66Itfy244ElZngpMMHPMhukubrYFD00050KuthYPQStizXQRqERKKVWKKV+19blmrNN4pVSP/gZ/1IR2YnhUdMvxRNji8fJSqnljRTvODAAmAJMM2VtiAQgUBnrREqBmY0s8wrgMz/jpgITRGQXxsDsGj/T7QUqHClNwXiJ/SUKKDTDBY1M6xdaF2vQFD0ErYs10EahgyMiMcBLGDUnv1BKLVBKDQOOAJf4mewGoN6N4Osoy6mUKlFKlWPM7BrWUBoM5ayo4S4HBjey2BnAAj/j3gj8Tyk1FPgC+Jmf6f4H2EVkGUYttzG7phQAEWY4gnomKJxKdGRdbKIegtbFGmij0IExB4rmAb9VSmX4maby4FUh4K+//YHAnSKyCBgqIn5teiwi4ZUOz8ZoIjfEGqCiJncGcNBPGSt2JXMppfxtegtQ4VM7G4j0J5FSyqOUusesQXow+vL9ZS3G3h0A04B1jUjbIenouthEPQStizXQRqGNEBGbiCwFUoDFYmyr2BCzgLHAs+aMBX/6pi8UkW9F5FuMprFfMzGUUg9V6q/epZR60Z90wLkisklE1gBHlVLr/ShrK+AQkRUY9/exn2WB0byvdX/iOngP+KlZ1vUYg6QNIiJJ5jNfDqxRSh2tJ26V3xZj29cyEVkFeJRS3zdC3laliXoIHV8XG62HZllb0bpYNQ9zVFqj0Wg0Gt1S0Gg0Gs1JtFE4DRCRydKOztHMhTwXNzLNa+ZCKL8WCzVBpltE5KCIzKt0LlFElouxoKpV3KZ3drQu1pr/KaWL2ihomo1SapFS6ovK50SkId0aqJQ6p3L/qB9pGsMC4Pxq5x4Gfg9cALT5dE5N66N1sfloo3D6cKaIfGXWPmLMAadlIrJSRD4REYuI9BORtWIsl38UQER+YS74WSUioypnKCJzROR9EVkqIq+b584wBw/XV8rjJjOfZDPvj4GbRGSuGfebyi+ZWZMcISILzZrlAhFZAEwXkYdE5DvzPnqa8XeIyDzz+1oxXAVsEpHudT0MpVQ2UF7t9HCMgbpioEhEImqm1LQAWhcrccrpoj/LrPWnY3+AycASM3w1Ri1EALt57mmMmsqtwE3mOQHiMGoxAsQAn1fLdw7wmBl+BWN1pJ2TExS+MY9vAn6B4WpgF2DBmPWwrKKsWmReXUn2RWa4K7DYDJ8DvGKGMzGW758FbMaozFwH3N/Ac0mmkusDYGWl8DygZ3v/dqfbR+viqa+LuqVw+rDF/N4K9MNYTv+6GNMBrwISgY8wakXvAhdiLJ9PwXihPqX21Y7V8+0NfGnmOxjoUi3+NmXMq3YDc81+1KcbaI5vNr+TObmcf6NZHsA+ZbhYOAbsVkp5zXB05UxE5F0xpu6dUUc53krh02ZRWQdE6+IprIvaKJw+pFT63o/h0+VHpdQk4BOMGphbKXU/cDPwJMZCnQ1KqclKqcnU7PesLd87gb+Y+e6jpmtmL4AYPmTeV0r9DIjHmANeFxUvSFql8sZwcgFS5XnTlcNVylZKXW/ey9Y6ytkuIhNEJBSIUEoV1hFP0zy0Lp7Cuqj3Uzh9cIuxAjQYuBKjKf2YiIzBWOqeClwqIncDIRhN2RNmn+hKjJWSy4GnquU7WIzl9elKqbUiEgm8JCI/AK565AkHFpgvZCGwo6EbUEplmH2+a8y8b/T/9qsiIpdgdF30FZFPlFJXAs8Cb2M8myeamremQbQuVuJU00W9eE1TJyIyB6O/dWl7y6Lp3GhdbDt095FGo9FofOiWgkaj0Wh86JaCRqPRaHxoo6DRaDQaH9ooaDQajcaHNgoajUaj8aGNgkaj0Wh8aKOg0Wg0Gh//D924UTVXw4KhAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 381.6x252 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "### Draw figure\n",
    "\n",
    "# Set replicate\n",
    "rep = 'rep2'\n",
    "primer = 'NpN'\n",
    "\n",
    "# Set figure font properties\n",
    "from matplotlib.font_manager import FontProperties\n",
    "figtext_font = FontProperties()\n",
    "figtext_font.set_weight('semibold')\n",
    "figtext_font.set_size('9')\n",
    "\n",
    "# Rows index positions\n",
    "poss = [6,7,8]\n",
    "poss_labels = [f'{pos}-{pos+1}' for pos in poss]\n",
    "nrows = len(poss)\n",
    "\n",
    "# Concentrations to use\n",
    "conc = 160\n",
    "ncols = 2\n",
    "\n",
    "# Set figsize\n",
    "figsize=[5.3, 3.5]\n",
    "\n",
    "# Create figure\n",
    "fig = plt.figure(constrained_layout=True, figsize=figsize)\n",
    "gs = fig.add_gridspec(nrows=nrows, ncols=ncols+1,\n",
    "                      width_ratios=[.3,1,1],\n",
    "                      left=.05, right=.99, bottom=.05, top=.90,\n",
    "                      wspace=.2, hspace=.2)\n",
    "\n",
    "# Draw logos within figure\n",
    "t = time.time()\n",
    "\n",
    "# Create dict of logo csv files\n",
    "logo_files_dict = {}\n",
    "for i, pos in enumerate(poss):\n",
    "    logo_files_dict[(i,0)] = f'csv_logos/invivo_p{pos}_{primer}_{rep}_logo.csv'\n",
    "    logo_files_dict[(i,1)] = f'csv_logos/invitro_p{pos}_{conc}uM_{rep}_logo.csv'\n",
    "\n",
    "col_titles = ['in vivo', 'in vitro']\n",
    "    \n",
    "# Iterate over rows\n",
    "for i, pos in enumerate(poss):\n",
    "\n",
    "    # Print position for feedback\n",
    "    pos_label = poss_labels[i]\n",
    "    print(f'{pos_label}: ', end='')\n",
    "    \n",
    "    # Write position\n",
    "    ax = fig.add_subplot(gs[i,0])\n",
    "    ax.set_ylim([0,1])\n",
    "    ax.set_xlim([0,1])\n",
    "    ax.text(x=.5, y=.5, s=pos_label, ha='center', va='center', \n",
    "            fontproperties=figtext_font)\n",
    "    ax.axis('off')\n",
    "    \n",
    "    # Iterate over columns\n",
    "    for j in [0,1]:\n",
    "    \n",
    "        # Provided feedback on progress\n",
    "        print('.', end='')\n",
    "\n",
    "        # Load logo\n",
    "        logo_file = logo_files_dict[i,j]\n",
    "        logo_df = pd.read_csv(logo_file, index_col=0)\n",
    "\n",
    "        # Draw logo\n",
    "        ax = fig.add_subplot(gs[i,j+1])\n",
    "        logo = logomaker.Logo(df=logo_df,\n",
    "                               ax=ax,\n",
    "                               font_name='Arial Rounded MT Bold',\n",
    "                               show_spines=False,\n",
    "                               flip_below=False)\n",
    "        ax.spines['left'].set_visible(True)\n",
    "\n",
    "        # Color logo\n",
    "        logo.style_glyphs(color=red)\n",
    "        logo.style_glyphs_below(color='k')\n",
    "        for c in 'ACGT':\n",
    "            for k in range(1,11):\n",
    "                if (k != pos-1) and (k != pos+2):\n",
    "                    logo.style_single_glyph(c=c, p=k, color='gray')\n",
    "                  \n",
    "        # Style plot\n",
    "        if i==0:\n",
    "            ax.set_title(col_titles[j], style='italic', fontproperties=figtext_font)\n",
    "        if i in [0,1]:\n",
    "            yticks=[-2,-1,0,1,2]\n",
    "            ylim=[-2,2]\n",
    "            rect_height = 1\n",
    "        elif i in [2]:\n",
    "            yticks=[-4,-2,0,2,4]\n",
    "            ylim=[-4,4]\n",
    "            rect_height = 2  \n",
    "        ax.set_yticks(yticks)\n",
    "        ax.set_ylim(ylim)\n",
    "        ax.set_xticks(range(1,11))\n",
    "        if j==0:\n",
    "            ax.set_ylabel('enrichment', labelpad=.5)\n",
    "        ax.set_xlabel('base pairs from -10', labelpad=.5)\n",
    "        ax.tick_params(axis='x', length=0, pad=0)\n",
    "        \n",
    "        # Gray out positions with zero entries\n",
    "        for l in logo_df.index:\n",
    "            if l not in [pos, pos+1]:\n",
    "                if (logo_df==0).all(axis=1)[l]:\n",
    "                    patch = Rectangle((l-.5, ylim[0]), 1, height, \n",
    "                                      facecolor='lightgray', edgecolor='lightgray')\n",
    "                    ax.add_patch(patch)\n",
    "        \n",
    "        # Create rectangular patches\n",
    "        rect_kwargs = {'linewidth':.5, 'edgecolor':'k'}\n",
    "        ax.add_patch(Rectangle((pos-.4,-.5*rect_height), .8, rect_height, facecolor=pink, **rect_kwargs))\n",
    "        ax.add_patch(Rectangle((pos+.6,-.5*rect_height), .8, rect_height, facecolor=purple, **rect_kwargs))\n",
    "\n",
    "    # Feedback\n",
    "    print('')\n",
    "    \n",
    "# Show execution time\n",
    "print(f'Execution time:{time.time()-t:.2f} sec.')\n",
    "\n",
    "# Tighten layout and save\n",
    "fig.show()\n",
    "fig.savefig(f'figures/fig_4_{rep}_logos.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.0"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}