04_join_annotate_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Joining together various data types into one big DataFrame\n",
    "This notebook will accomplish many things to put together the data into a big DataFrame for downstream analyses/figures:\n",
    "* Show that the two Rho libraries are calibrated\n",
    "* Join together the two libraries with Rho and Polylinker\n",
    "* Compute fold changes between WT and MUT sequences in the Rho assay\n",
    "* Annotate sequences for their activity as WT, MUT, and Polylinker\n",
    "* Annotate sequences for binding patterns with CRX, NRL, and MEF2D ChIP data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import itertools\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from pybedtools import BedTool\n",
    "from IPython.display import display\n",
    "\n",
    "sys.path.insert(0, \"utils\")\n",
    "from utils import fasta_seq_parse_manip, modeling, plot_utils, sequence_annotation_processing\n",
    "\n",
    "data_dir = os.path.join(\"Data\")\n",
    "figures_dir = os.path.join(\"Figures\")\n",
    "all_seqs = fasta_seq_parse_manip.read_fasta(os.path.join(data_dir, \"library1And2.fasta\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_utils.set_manuscript_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show reproducibility of raw data for all experiments in one figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAJBCAYAAABBBGGtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXt4XFXV8H8rTVrStGmotEC0JAVS+31FKVARolwEgm0FxIxcCiKavIAi8vGKAZWKBQqoeauVV6Eg5OUiBbQXbulbykWuCljQIiVDI5dYKUiBJmnTlt7W98c5ZziZTOaSzGRmzqzf8+wnc/bZ55w1k3XOWXvttdcWVcUwDMMwDKPQKMq2AIZhGIZhGNnAjCDDMAzDMAoSM4IMwzAMwyhIzAgyDMMwDKMgMSPIMAzDMIyCxIwgwzAMwzAKEjOC4iAic0REfeUdEXlQRD49xHLsLyI3ishLIrJTRB4fyutHyTJMRC4VkadE5H23rBCRz2RLJmNgmH73K88VIvJ3EekWkY0islJETsumTEbq5JB+nyIi94vIWyKySUReEJFZQylDlDznuzq9QUQ2u7p+vohItmTKJsXZFiAP6AKmu5+rgSuBh0Xk/6jqB0MkwxRgJvAsUDJE1+yPUuAHwP8A1wIKXAA8LSK1qvpCNoUzUsb0uy/lwK3AK8BO4KvA3SKyU1UXZVMwI2VyQb+/B7wB/CfwHo6uLxSRPVT1v4dIBj+7A0uBl4DNwLHAr4GRwH9lQZ6sIpYssX9EZA5wgaru4as7DPgzcKaqLhwiOYpUdZf7eRGwh6oePRTXjiHLMKBcVTf46oYDa4A/quo3syGXkTqm38kjIs8A76vqSdmWxUiOHNLvPVT1vai6hcDhqjpxKGRIhIjcCXxKVYfUS5YL2HBY6qxy/07wV4rIWBG5SUT+LSJbReRPIvLZqDYqIt8TkV+JyAci0iki/+0aEf3ivSBSQUQmutf7UlT9MNctPNfd/oSI/F5E3hWRLSLymohcFUeWnX4DyK3bBqwGKlOV08g5Clq/4/A+EPd7GHlBNvT7vRjVfyXO81JEykSkR0S+E2PfX0Tkd+7nChG5WUTWuXL/U0R+G0+efihY/bbhsNTZx/37hlchIiOAR4AKoAl4F/g28IiI1KjqO77jL8Zx+5+JMwxwNbDVPS5tqOobIvI8cCrQ6tt1FLAncLe7fTvOENe5QCewLzA5lWu53/9gwIYK8h/TbxcRKQZGAV8CjgdOT4vwRjbJFf0+HMd7HhNV7RGRB3H0+zc+WfcFpgFXuFW/AGpxhtrewTHujkxGAFe/dwOOAL4O/DjF7xAMVNVKPwWYgzOGW+yW/YCHcaz4Eb52jcA2oMZXVwy8BjT76hQIA0W+ustwxmXHJinTIuDxJNv+J86D3y/rjcDLvu1NwImD/J2uBD4EPpnt/5mVlP5vpt/9n/sw9/sosB04L9v/Lysp/w9zTr/dY44FdgHfSNDuKzgxaZW+uh8CHwAl7vbLwHcH8Nvs5dNvBa7K9v8rW8WGwxLzMZyH4HbgH8BBQL2qfuhrcxzwAvCGiBS7FjbAEzhWu5/7tLf7fwlOT/WADMj+e5wgz+kQsfzrgXt8bf4GXCsi3xCRffqeIj7ucMRlwKWq+urgRTaGGNPv2Pwd+AxQhxM0+mvJ4oweY8DklH6LSDWw0D3PrQma/y+OEX+Kr+40YKmqbne3/wY0iTO7a1IyMri8h6PfX8DxKjWJyCUpHB8YzAhKTBeOshwGnIczbrpQRPy/3R7u/u1R5ZtEjT3juFpjbe+dXrFBVd8Cnsa5ccDpgezBR0MFuPtWAr8EOkTkbyJybDLnF2da/D3AAlWdnzbBjaHE9Dv2uXtUdaWqPqKq/wncAfwsrV/AGApyRr9FZCyOYdOBM5wWF1XdCtyHq98i8kngQHrr9wXAvcDlwKsi0i4iCYdtVXWHq9+Pq+oc4BrgChEZmejYoGExQYnZoaor3c/PicgWnDiDU/iox/kBzoP22zGO/zBqe3w/22+nQdZY3AP8VERKcW6mv6pqu7fTfZF8w30oHIrjQr5fRPZR1ff7O6nb62gFHgUuzJDsRuYx/U6OF4Fvikixqu5In/hGhskJ/XaNiwdxjLATVHVzcuJzD/CA68U8DVgPPObtVNVOnOfvheLkP7oEuFNEXlLVV5K8Bjj6vRtOsPY/Ujgu7zFPUOr8Dmcm1KW+ukeB/YF/uta1v/w96vgvR/VC6oEtOGO7meAPOO7ar7jl7liNVHWXqj6L4xodCVT1d0IR2Rt4CGfMfJaq7ky30EbWKHj97ofPAf8yAyjvGXL9dofX/gDUANNVNdqbFI8VOHFvp+IYQYv6e96q6ks4AdpFpDi5BUe/PwTWpXhc3mOeoBRRVRWRa3Cs7WNV9VGcnsW3gMdF5L+A13HGog8F3lHVX/pOMRr4gzuNcQpORP5vNE7iLrcXMdPd/DhQLiJfdbeXxetVqOq74mTg/S+c2Q+/9513DI4xczvOTIUROLMf3gHa+pGlFMeluzuOK/bT8lGi0Q9V9a/9yWLkPqbfUgW04BhTr+HMDvsKzsywWJ4CI4/Ihn4D1+Po9/8DPiYiH/Pt+2tUfFK0vNtFZAlOwsW9gfP9+0XkaZzEhy/jBDifA/QAz/d3ThH5C3Ab8CpOctI6nGf5vBQ8VMEh25HZuVxwZxfEqB+G81B9yFc3BvgVsBZnpsG/cILmPudrozjK/GtgA8549W/wzVToR45qekfy+0t1Et/jP9y2f46qHwH8Fudm2IwTLPcgTtKsgcjyZrb/Z1aSL6bfMc81Bif+5w2cqc/v4Aw/zMz2/8tKaiWH9PvNQer3cW7bt/DNTHP3NeME8W/E8Rj9ETgiwfl+635/7574M/A13OTJhVYsY/QQIiKKM53x19mWxTDSjem3EWRMv4OJxQQZhmEYhlGQmBFkGIZhGEZBYsNhhmEYhmEUJBnxBInIBSKyUkQ+FJFbffXV7iJ0m3zlx779I0SkRUS6xVkE8XuZkM8wDMOIjT2/jUIiU1Pk1wFzgS/i5PCIpkJj59uYg5NLoQpnbZM/isgrqro8Q3IahmEYvbHnt1EwZMQIUtUlACIyDfhECoeejbOo3AZgg5uL4RtA3Jtojz320Orq6oEJa+QFL7zwwnuqOi7bcmQD0+/gk0v6PdTP76KiIi0tjWVrGUFh8+bNqqo5GYOcrWSJHe50w4eBJlV9T0R2x0kGtcrXbhVwcqKTVVdXs3LlykTNjDxGRDqyLYOHiFyA83D/FHCXqn7Dt+9YnNwh+wDP4bwUOtx9I4AbgK/i5Oj4uar+ItH1TL+DTy7pdxKk9fldWlpKT09PZiQ1cgJ3uZKcZKgtM2/l2irgEJzsm3e6+0a5f7t87bvcNn0QkXPdceuV69evz5C4hhETb7igxV8pInvgJFj7MTAWZz0i/4rmc/houOALwCUiMn0I5DWMdJCR5/eOHbYSiZE9htQIUtVN6qzHskNV/42Tqvt4ERkNbHKblfsOKcfJhBnrXDep6jRVnTZuXGIv8vPPP8/hhx/OkUceyaxZs9i+ffsgv41RqKjqElW9F4hegLMeWK2qf1BnBeg5wIEi4q3jczZwlapuUNU2nMyt30iHTKbfRqbJ1PO7uDj+gITptpFJsj1G583PL3LHkd8GDvTtPxBnsbtBM2HCBB577DGefPJJqqurue+++9Jx2rwhHA7T2NhIOBzOtihBZgq+4QBV7cFZf2pKnOGCKbFOlKqns9D1O9sU6P01JM9v020jk2QkJshdNbcYZ42WYSKyG7ADx4XaCbTjLMB5HfC4qnou1NuB2SKyEtgTZzG4b6ZDpr333jvyefjw4RQVZdv+G1qam5tpaXFGb2655ZYsSxNYRgHRFos3JJDScIGq3gTcBDBt2rSEybwKXb+zTZDur1x7fptuG5kkU9o0G9gC/ABnYbYtbt2+ODMFNuKsevshMMt33E9wes4dwBNAc7qnV3Z0dLBixQpOPPHElI/94IMP+MpXvkJZWRlVVVUsXLiw37ZtbW0cc8wxjBkzhv3335+lS5f22j9q1KheZdiwYXz3u99NWaZkaWpqoqGhgaampoxdw2ATvYcD4KMhgZSGCwZKoep3tgnY/ZWTz2/TbSMjZHsF13SUQw45RBOxY8cO7erq0iOOOELD4XDC9rE4/fTT9dRTT9WNGzfqU089peXl5fryyy/3abd9+3atqanRefPm6Y4dO/TRRx/VkSNH6quvvhrzvBs3btSysjJ94oknBiRXIQCs1BzQNX/BCY6+1bd9LvCMb7sMZxbYZHd7HVDn238lcHei65h+B59c1O+hKiNHjoz725hu5z9Aj+aArsUqWRcgHSXWS+Lmm2/W4447ThsaGrSiokJ//vOf64wZM/SRRx6J97/ql02bNmlJSUmvm+FrX/uaXnrppX3a/v3vf9eysjLdtWtXpK6urk5nz54d89y33nqrTpw4sVd7oze59JLAGSrYDbgWuMP9XAyMwxniCrl1PwOe9R33U5we8u7AZJwYiumJrmf6HXxySb+HukQbQabbwSOXjaDADq6uWrWKZ599li9/+cu8//77jB8/nueee46rrrqKo48+mnvu+Wjm8gknnEBFRUXMcsIJJwCwZs0aiouLmTRpUuS4Aw88kNWrk4v7U1VefvnlmPtuu+02vv71ryMig/jGxhASc7hAVdfjGEBXAxuAzwKn+45L23CB6bcRVEy3jSEl21ZYOkqsnvKRRx6pc+bM6ccuTZ0nn3xS99xzz151N910kx511FF92m7btk0nTpyoP/vZz3Tbtm360EMPaUlJiR5//PF92r755ptaVFSkr7/+etpkDSIUcE/Z9Ds/aWtr04aGBm1ra0vYtpD1O9oTZLodPDBP0NDz0ksvccopp6TtfKNGjaK7u7tXXXd3N6NH953cU1JSwr333ktrayt77bUX8+bN49RTT+UTn+ibgf6OO+7g85//PBMnTkybrEbwMf3OfbwZY83NzdkWJa8w3TaGkkAaQR0dHWzfvp3JkycnbgzMmDGjT8S/V2bMmAHApEmT2LFjB+3t7ZHjVq1axZQpMdO88OlPf5onnniC999/n4ceeojXX3+dQw89tE+722+/nbPPPnsA39IoVEy/84OAzRgbEky3jSEn266odJTo4YL77rtPDzvssEQeupQ57bTT9PTTT9dNmzbp008/3e8MA1XVVatW6ZYtW7Snp0ebm5u1urpat27d2qvNM888oyNHjtTu7u60yxo0KODhAtPv4FPI+u0fDjPdDibYcNjQsmrVKqZOnZr2815//fVs2bKF8ePHM2vWLG644YZevYkZM2ZwzTXXAI6rdO+992b8+PE8+uijPPzww4wYMaLX+W677Tbq6+tjumUNoz9Mv42gYrptDDXiGGn5zbRp09RW2Q42IvKCqk7LthzZwPQ7+BSyfpeVlamtIh9sRGSzqpZlW45YBNITZBiGYRiGkQgzggzDMAzDKEjMCDIMwzAMoyAxI8gwDMMwjILEjCDDMAzDMAoSM4IMwzAMwyhIzAgyDMMwDKMgMSPIMAzDMIyCxIwgwzAMwzAKEjOCDMMwDMPICUSkWERahup6ZgQZhhFowuEwjY2NhMPhbItiGEZihgFnD9XFiofqQoZhGNmgubmZlhanY3nLLbdkWRrDMETksTi7hw2ZIJgRZBhGwGlqaur11zCMrPNZ4Frg7Rj7SoDPD5UgZgQZhpGThMNhmpubaWpqYvLkyQM+z+TJk80DZBi5xd+AsKouit4hIiOA64dKkIzEBInIBSKyUkQ+FJFbo/YdKyJhEdksIn8UkSrfvhEi0iIi3SLyjoh8LxPyFToWI2HkA94wVnNzc7ZFKSjs+W0MAfOBD/rZtx345lAJkilP0DpgLvBFoNSrFJE9gCXAfwAPAFcB9wCHuU3mADVAFbAX8EcReUVVl2dIzoLEYiSMfMCGsbKGPb+NjKKqf4izbxdw21DJkhFPkKouUdV7gfejdtUDq1X1D6q6FeemOVBEPF/32cBVqrpBVduA3wLfyISMhUxTUxMNDQ32cjFyGm8YK9WhMPN0Dg57fhuFxFBPkZ8CrPI2VLUHeA2YIiK7A3v797ufpwyphAXAQF8uhjEUDNaIsWG0jGHPb6Nf8rXzMdSB0aOA9VF1XcBod5+3Hb2vDyJyLnAuwD777JNeKQ3DyBqDHa61YbSMkZHn9/Dhw9MrpZEV8jXMYqiNoE1AeVRdObDR3edtb43a1wdVvQm4CWDatGmadkkNw8gKgzVibDZYxsjI87usrMye3wEgXzsfQz0ctho40NsQkTJgP5xx5g04OQMO9LU/0D2moMlXN6NhDAQbrs1Z7Plt9MtA7lsRWSoiJ4tISQZFi0umpsgXi8huOJkfh4nIbiJSDCwFDhCRkLv/cuAlVfXe7rcDs0VkdzfY7hzg1kzImE9YjINhpA/rVMTHnt/GEPIUjh69IyI3iEhtqidw9e3rIvJD9+/YVI7PlCdoNrAF+AHwNffzbFVdD4SAq4ENOFkjT/cd9xOcQLsO4Amg2aZX2mwuw0gn1qlIiD2/85R8M/BV9ReqejBwJNAJ3CUi7SJyuYjsl+h4ETkcR+e+BXwaOA/4h1ufFKKa/8Ox06ZN05UrV2ZbDCODiMgLqjot23JkA9Pv9JKuTNTppJD1u6ysTHt6erItRiBobGykpaWFhoaGnIqLE5HNqlqWRLsjgF8DB+DEmf0FuFhVV/XT/jngl6p6t6/uNOD7qvqZZGSzVeQNwxhyBtpjTUdP12KOjKCSj6MGIvJJEblKRF7DCZa/B6gG9gSWAffGOXwS8PuoukXA/sle34wgwzCGnIEOSfV3XL4NAxhGJsg3A19EVgLPAGOBM1T1/6jqNaq6VlW3quovEpyind5DsgCn4AyRJYUtoJpH5KIb3zAGQvR02mR1u79puPmao8QwChUREeBu4DpV3dZfO1WdGOc0FwEPisiFOLFo1ThLt5yQrBxmBOUR9qA3gkJ0Lp9kdbu/HED5mqPEMAoVVVURuQJI5O2Jd44/uQHUXwIqcda0W6aq/S3O2gczgvIIe9AbQaWpqYmuri46OzsJh8MpezotQaJh5CV/xYnrGfA4tpuj6ncDPd6MoDzCHvRGUJk8eTJjxoyhpaWFiooK03PDKAweB5aLyK3AWiAyXV1VWxIdLCITcVI2TOWjpVu845NaT8sCo7OMBXQGDxF5XES2isgmt7zq23eGiHSISI+I3JtqYq8g0J/O5+PMFsMwBsXngDeAo3ByUp3llq8lefxCYBdwse9YrySFeYKyjMX5BJYLVPVmf4WITAFuxBm/fhFnOuj19J3dEGj603nzdBpGYaGqXxjkKaYAn1PVXQM9gXmCMkx0rzd623q/BcWZwAOq+qSqbgJ+DNSLSMyVtoNKU1MToVAoEv9jGEbyBHX0QByKvJLkYU8CBw3muuYJyjDRvd7o7f56vzYdPu+5VkR+CrwKXKaqj+P0Wv7kNVDV10RkG05g4AtZkTJLrFy5ko6ODkSERYsWmb4bRpIEafRARD6OkyH6SKAiavewfo650rf5Jk5M0VLgHX87Vb08GRnMCMow0TO6kp3hFSRFL0AuBV4BtuEMdT0gIl7gXldU2y6gjydIRM4FzgXYZ5+k4vvyhubmZjo6OgDwlu0xfTeM5AjYLOEFwGbgWJz15o4E5uBkiu6PCVHbDwIlMeqTQ1XzvhxyyCGaa7S1tWlDQ4O2tbVl5figAazUHNC1gRRgOfBd4D7gkqh9G4FD4h2fi/qdCp4ut7a2Rv6GQiGtr6+P6Heh63s+6/dgy8iRI9PxE+YVhabvQI/Gfja+D5S5nzvdv2OBcKz2mSjmCcoQg+3ZWpBooFBAgNXAgV6liOwLjADWZEmuIcG7F5566ina29sBWLRoUZalMozsYZ7PCDuBHe7nThEZB3QDHx/oCUXkU8DlqnpKMu0tMDpDDDTgOahBb4WCiFSIyBdFZDcRKRaRM3FcvMuBO4ETReQIESkDrgSWqOrGbMqcKTxdPvTQQ6mpqeHiiy/u954Y6FpihpFvhMNhOjs7CYVCQRnSGgzPATPdzw/hLJ66BFgZ7yARGekuuvqAiPxCRMpFZF83NujPwLvJCmBGUI5hL4O8pwSYC6wH3sMZBjtZVdeo6mrgWzjG0Ls4sUDnZ0vQTLFs2TKqq6s59thjaWlp4ZJLLqG9vZ3nn38+0uuNNvRtlqQRJOJ1Zpubm1myZAljxoyxSQBOPp8n3M8XAY8BLwNnJDjuN8CJOLGXxwGL3fOsBqpV9TtJSzBU426ZLLkYM9HQ0KCANjQ0ROqSGQcutLHiZKGAYyZyUb/jUVNTozhDgFpWVhb5W1tbG4kFir43Cp1C1u8gxgTFev57FOIznhgxQTizv24DRkTvS1SAdcB49/MncBImHpHqeVTVjKBMEUvR490YRnwK+SWRi/odj9bW1ojxU1dX18soAjQUChXcSyARhazfQTSCBmPoBNFIimUEOdW8DZTE2hevAN3xtlMpNhyWRvwuUC+w2e/uNJe/ERTiJQGdOXMmK1eupKGhgeuuu4758+dTVVVFbW0toVCIhoaGLEtvGJkl1vM/WQosJOKXwBUiUpLiccUi8gUROUZEjgHwb3t1STFQ6ymXSrZ7ym1tbVpfX6+VlZWR3m/QLPlsQwH3lLOt37GI9mp62/X19X103xv+qq+vj3msUdj6HURPkOrAwx8KzBO0FtgObHU//9Mrsdr7jnsTZ82x/srr8Y73F5siP0jC4TAnnXRSZOovwJo1a3j44YeBgp/+aAQML7NzKBQCIBQK0djYGNleu3YtS5YsobOzk6uvvprm5mY2bnQmv4kIELhkb4YRE39qiPvvvz+mVyjWVPkCS4+S7EKpvVDV6nQJYEbQIAiHw0yfPp2Ojg5KS0vZsmULlZWVTJo0iWnTpkUe8suWLeOiiy5i/vz5zJw5M8FZDSN38R7aXV1dqCrnnHMO69atY+3atUyY8FHCVhHhsssuY8mSJdTW1lJVVUV3d3evoWLDCDJNTU2R3FjNzc29dD66M1GoHQJVfSJxq+QQkftV9aSBCDHkBXgcx/21yS2v+vadAXQAPcC9wNhE58vWcIHn1ge0srJS6+vrI67/qqqqSEZcLzC0pqYmK3IGAQp4uCCXhsO8od+qqqpewc7etjcc1traGhke9v5iQ2D9Usj6na/DYYOZ7VtoQ8L0Pxw2HCdfWrv7zm8HrgJ2i9U+XgE+SPUY1ewOh12gqjf7K0RkCnAj8CXgReAm4Hqc9ZdyjqamJrq6uiKLQW7cuJE1a9ZQWVlJR0dHZIHI6upqtm3bxvz587MtsmH0S7xFTMPhMLNnz/YeNnR0dFBWVsaBBx5IWVkZANOmTaOhoYHFixfT0tLCunXrAJgyZQqHH344qlqwPd4gISKPA4fxUabft1T1k+6+M4BrgT2Ah4EGVf0gG3JmGs/T2dnZyeLFi2O26c/raUPCEW4APglciOP8qAJ+hJMxekhmUOTacNiZwAOq+iSAiPwYaBOR0ZqjWXXHjBnD9ddfz+LFi3n00Ufp6OigqqqKUCgUsTQffvhhGhoabCjMyGn6S+XvH/YFJw6osrKSdevWUVZWxoQJE2hpaaGmpoaWlhYWL15MfX099fX1bNy4kdGjRzN37lxLDBcs8r4TO1i8GDfvbyrYkHCEk4H9VLXT3X5FRJ4D/kECI0hEoleJ381fp6pXkgTZNIKuFZGfAq8Cl6nq48AU4E9eA1V9TUS2AZOAF7IiZRz8L42mpibWrl0LwDe/+U3uvPNO5s+fz7777ktFRYVZ/EbO01/v1L/qe2VlJa+99lrEy/P000/z+c9/nqqqKtrb2/nUpz7VKw2EN2mgoqLCHvrBJ+86sYOhoaGBl156KemUD/E8rQXMO8BIoNNXV4qTPygR0danxKhLzEDG0AZbgM/iLBkwAjgbZyXt/YBHgW9FtX0LODrGOc7FWV9k5T777JPqEOWAaGtr01AopHV1dVpfXx9ZFbutrU3r6ur6JIezGKD0QQHHTGQ7JsjT76qqKq2tre0VC+SVWGkhvLiHmpqaXveK0Zd80m+cmE5vWZhnvOczcB9waVTbTcAh8c6XrzFBfv1ua2tLGCNUaHFAfvDFBAHH+MoPgJeAc4AZ7nt9VbQeJVPIp5ggVX3Ot3mbiMzCWURtE1Ae1bwcx0iKPsdNOO5Wpk2bphkSFfjIgu/q6uoz9utt/+1vf4v8vfXWWyOzwQwjX/H0/tBDD+W5556ju7s7bvumpqbIjJeWlha6u7upr6+PTJW3VbMDw6U4azZtwxnqekBEpgKjgK6otl04Hd5eiMi5OC88hg8fnlFhM4V/9tf06dMBIh5TiwOKS6wHwI+its8DfpbieVP3AkFuJEsE/hcnMOoa4E5f/b44N9roeMdnuqfsTwRXW1urIqKA1tbWRmbL7LbbbpGZMEb6IY96yukuQ+0J8nq0oVBIAS0vL494fHbffXcdO3ZsTE+Q3wNK1GywICaASyf5rN/AcpyFgu8DLonat5GAeYL8uuyf/es9/03H+0I/s8PSWXCGYlM/LtOCxRC0AvgisBtOTNKZOFPjJuHEBHUDRwBlwO+AuxOdM1MvCf/wVygU0tbW1l4vhMrKyj6LRba2tmZElkInn18Sgy2ZNoKiH+redPfdd99dy8rK9Dvf+Y6WlJRE9NzrBHiluLg4MkxWWVmptbW1WllZqXV1dfZCSJJ81u/BdmJzzQjqz2CP7hx4hv+cOXO0uLg4kirC6MtQGEEDLdm4YcYBf3F7CJ3As0Cdb/8ZOGmze9yeRdbyBPnzAHkxEYCOGDEi8rD3tr02RmbI55fEYMtQeTr9D3h/GTt2bC8jyCtTp06N1Hv3hj9/UCgUyqjcQSJf9DsTndhcM4K8+yF6oV8vB5wX/+bpuncPlJeX64IFC8zjGQN6xwS1+T73WiqDJJfNSGcZ8pggVV0PfCbO/oXAwqGTKDbhcJjOzs7IVODHHnuMnTt3RvZffvnl/OhHzjCmv94w8g0vRiEUCnHOOef02f/BB33TvFRWVvLuu++yfft2ysrK+OEPf8jzzz/PoYceysUXX0xPT4/3kDOCRQkwF5gM7ATCwMmqugamMPKlAAAgAElEQVRARL4F3Al8DHgE+GaW5EyZ6CzOnZ2dtLS00NnZiYjw9NNPR9recsstLFu2jFmzZtHd3U1NTQ3333+/xb4lh/8hM6BlM9JJruUJygn864FNnTqVdevW9TJ0PvzwQy655JJIoOiOHU7OsPLy6Jhuw8h9vJwljY2NkanvidiyZQsbNmwAoKenhxUrVlBRUcGKFSvo6emhpqaGq6++OpNiG1kgXzqxA8G/JMyYMWNobGykoqKCtWvXRtaCBBg92onznjlzJnfddVdkEszkyZMt+DkJVPVp3+c+y2aIyDDgJ0DaltSIR9FQXCTfmD17Nu3t7YwYMSIy6yuafffdl7FjxwLO7Ia6ujrmzp07lGIaxqAJh8M0NjYSDofZZ599kj7OM4AqKyupr69HRGhpaUFEaGho6HfBSMPIVQ499FDKy8t5++23aWlp4fzzzyccDvPkk09G2lRVVdHY2Bi5ZxYvXkx7e3tklrDXoTDdHxTFwGVDdrWhGnfLZElHzIQX9LZgwQItLS2NmQslutTV1WlZWVlk5piROciTmIlMlEzEBHn6PnXqVAUisxtTKf6ZMNHB1RYXkRqFrN/Zigny579qbW2NTHIpLi7W8ePH96vv/hg60/XkIIXAaJz8gbuSbT/YYsNhLp4r9Pbbb48Mb8XDWy+pp6cHGFjqdMMYaqJzXnl6u3Xr1oTHFhUVsWvXrsj24YcfzuTJk/tkwm1sbLS4CCPnaW5ujgxzebFwRUVF7Nixg40bndR0I0aMYNeuXWzfvp0pU6b0GfKy5S8yxpAFFJoRhPNiCIfDkRsgGTzjxxsKsKEwIx/wjP36+npqampob29P+thhw4ax33770d7eTk1NTUTno4NBLS7CyHXC4TBr166lvLycTZs2sX79erZv3x7pFGzbto3S0lK2bNlCVVUVHR0dkZhPzxCyJTAGjogcE2f3kGbPNCMI5yH+pz/9KXFDF69HPHr06H5XDzaMXMRvoLz++uuEQqGkvEAABxxwAAsXLuzz8I82eqx3bOQ6fi8QEPFwusMx7Ny5ky1btlBeXh6Z+eg36m0W2KBJ9KP9c0ikwIygiBeopKSE7du3J2zvtbPZL0a+ED1c1dTURGNjI88++2yv4a1E/Otf/wL6PvTN6DHygXA4zGWXXcbGjRvp6elh7NixMdM/gBPeUFxcTHd3Nw8//DCLFi3qtd+8nYNDVSdmWwaPgjeCUvUCDR8+nLPOOsvcoEbe4PVaOzs7qaio4KGHHuKtt95K+nhvWGD9+vWcdNJJNvPLyCu8TkD0VPf+KC8vp7u7O9Ip9rxD/nM1NTWZ4R8QCt4ICoVC3HHHHUl5gYYNG8a8efM477zzhkAyw0gPoVCI5cuX09rayocffpjSsZWVlVx++eVce+21bN++nfb2dpqbm+0FYOQNXicg2YVae3p6qKurA5ycQH6Pvw2DBY+CNYLC4TCzZ8/mscceS8oAqqqqYvny5dYDNvKOW265JekkiNEcdthhnHfeeZx33nm9esGGkQ+Ew2FefPFFwAl2ToadO3dSXl7eawgsOpu03QPBoWCNoMsuu4wlS5Yk1Xby5MksXbrUDCAj77jxxhtZunRpSsdUVlYyZcqUPr1gi/0x8olly5ZRX1+fkvfTG/r1VgPwMA9QcClII+iKK65I2gACqK2tNQPIyDuWLVvGt7/97V4xDclw2GGH2axHI2/xvDZLly5NyQAaN24cU6dOjcQNNTY29ppMAOYBCiIFZwSFw2HmzJmTVFuvR2yKb+Qb4XCYU089NSUDaOrUqRx88MGm70beMhDvDxBZABUcr4+3eCoQWQbDPEDBpKCMoHA4zGc+0+/af72oq6tjxYoVGZbIMDJDXV1dJKFnMlRWVnLXXXeZx9PIW8LhMCeffHJSMZ7FxcXs2LGDgw46iIMOOqjXbN9bbrmFcDhMRUWFdQgKgIIxgsLhMEcddRSbNm1K2La0tJTrrrtuCKQyjPThDQOUlpZGcvokw5gxY3j00UfNADLyiuhA/WOPPTYpAwjg17/+dSQBYiy9N89P4VAwRtDs2bN59913E7YrKyvj97//vb0QjLwjlWB/P8cddxzQOwbCMHIdL1j5wQcfTOrZ7lFbWxuZ8WgYRdkWYChYtmxZ0i+H2tpaZs6cmWGJDCP9PPbYYym1Lyoqora2lrlz50ZeKM3NzRmSzjDSQzgc5vjjj+eBBx4ASNoAGjNmDPX19ebhMXoReE/QjTfeyLe+9a2E7SorK1m3bl1kkTzDyBfC4TCNjY10dnamdNzJJ58cmQVms1+MfOHCCy9MKvOzR1FREcceeyzXXXedeTmNPgTeCErGACopKeG3v/0tixcvtpeAkXdceOGFSS/9UlVVxfXXX99H1y0GwsgHbrzxxqQNoOLiYsaPH8+UKVPMADL6JdBG0Pjx4xO28WKAZs6cacNgRt6xbNmypF8KZWVlkaznputGvnHFFVcknd4E4Atf+AITJkyIDPOakW/EIrBG0NixY9mwYUPcNnPmzOEnP/nJEElkGOll2bJlfOlLX0qqbV1dnfWGjbzlzDPPZOHChUm1HT58ONu2baO8vNyGeY2E5GRgtIiMFZGlItIjIh0ickYqx5955plmABk5yWB120+yBlBtbS0rVqwwA8jIS5YtW5a0AVRXV8eqVatoaGhg7ty5kWFe032jP3LVE/QbYBuwJzAVaBWRVaq6OpmDE90wU6dONQPIyBaD0m0PEUmqTX19PXPnzh2QoIaRbT71qU/x8ssvJ9W2pqYm4u20oS8jWXLOCBKRMiAEHKCqm4CnReR+4CzgB0kcH3f/1KlTueuuu9IhqmGkxGB1O8Vr8eCDD1rsjzGkiMhY4BbgeOA94IeqmpwbJ4ri4mJ27tyZsF1tbS177bUXV199tXl8jJTJOSMImATsUNU1vrpVwFGDPXFra6u9FIxskhbdTmToH3TQQSxcuNBeCEY2GDJPZ1FREatXrzY9NwZFLsYEjQK6o+q6gNH+ChE5V0RWisjK9evXJzypGUBGDpCUbkPq+u3R2trKiy++aC8GY8jxeTp/rKqbVPVpwPN0pnKehG0qKirMADLSQi4aQZuA6IyF5cBGf4Wq3qSq01R12rhx4/o9mYigqmYAGblAUroNyeu3x2c+8xnTcyPb9OfpnJLOi8yZM4cNGzaYAWSkhVwcDlsDFItIjaq2u3UHAkm5U1U1Y4IZxiAZlG57mI4bOUpKnk7gXHCmtCeDF89pxo+RTnLOCFLVHhFZAlwpIv+BM678ZaA2u5IZxuAw3TYCTkqeTuAmgLKysrhWvRn9RibJxeEwgPOBUuBd4C7g26kG1hlGjmK6bQSViKfTVzcgT6e/GEYmkSAomYisBzp8VXvgTM80kiMffq8qVU0cHBNAYug35Mf/bDAE+fvF+m6B0G8RuRtQwPN0LgNq4xn6IrIL2BJVXQzsyJScASTXf69SVc1Jp0sgjKBoRGSlqk7Lthz5gv1e+UfQ/2dB/n4B/25jgRagDngf+MFA8gQF+TfKBPZ7DZyciwkyDMMw8hNV/QA4OdtyGEay5KR7yjAMwzAMI9ME1Qi6KdsC5Bn2e+UfQf+fBfn7Bfm7pQv7jVLDfq8BEsiYIMMwDMMwjEQE1RNkGIZhGIYRFzOCDMMwDMMoSAJlBInIWBFZKiI9ItIhImdkW6ZsIiIXuItwfigit0btO1ZEwiKyWUT+KCJVvn0jRKRFRLpF5B0R+d6QC2/EJGg6LiKPi8hWEdnklld9+85wv2OPiNzrTr/OWex+GzxB0+/BYjqVeQJlBAG/AbYBewJnAjeISFoX78sz1gFzcfJ2RBCRPYAlwI+BscBK4B5fkzlADVAFfAG4RESmD4G8RmKCqOMXqOoot3wSwP1ON+KsQL4nsBm4PosyJoPdb4MniPo9GEynMkxgAqNFpAzYABzgrWIsIncAb6nqD7IqXJYRkbnAJ1T1G+72ucA3VLXW3S7DyWB7kKqGRWSdu3+Fu/8qoEZVT8/KFzCAYOq4iDwO/E5Vb46qvwaoVtUz3O39gDbgY6raZy2qXMLut4ERRP1OF6ZTmSNInqBJwA7v5nFZBRRyL6I/puD8NoCzsCfwGjBFRHYH9vbvx37HXCGoOn6tiLwnIs+IyNFuXbSOvobjIZiUBfkGi91vyRFU/c4EplNpIkhG0CigO6quCxidBVlynVE4v40f77ca5duO3mdklyDq+KXAvsDHcXKdPOB6feLpaL5h91tyBFG/M4XpVJoIkhG0CSiPqisHctp1niXi/VabfNvR+4zsEjgdV9XnVHWjqn6oqrcBzwAzCdZ3tfstOYL0P880plNpIkhG0BqgWERqfHUHAv2uXlzArMb5bYDIePJ+wGpV3QC87d+P/Y65QiHouAJCXx3dFxiB8xvkG3a/JUch6He6MJ1KE4Exgtwx0SXAlSJSJiKfA74M3JFdybKHiBSLyG7AMGCYiOwmIsXAUuAAEQm5+y8HXlLVsHvo7cBsEdldRCYD5wC3ZuErGD6CpuMiUiEiX/T0UkTOBI4ElgN3AieKyBHuA/5KYEkuB0Xb/TY4gqbf6cB0aghQ1cAUnKmC9wI9wD+BM7ItU5Z/jzk4PWt/mePuOw4IA1uAx3Fm4njHjcCZktkN/Bv4Xra/i5XI/yYwOg6MA/6C46bvBJ4F6nz7z3C/Yw9wHzA22zIn+D52vw3+NwyMfptO5UcJzBR5wzAMwzCMVAjMcJhhGIZhGEYqmBFkGIZhGEZBYkaQYRiGYRgFiRlBhmEYhmEUJGYEGYZhGIZRkJgRZBiGYRhGQWJGkGEYhmEYBYkZQYZhGIZhFCRmBBmGYRiGUZCYEWQYhmEYRkFiRpBhGIZhGAWJGUGGYRiGYRQkZgQZhmEYhlGQmBFkGIZhGEZBYkaQYRiGYRgFiRlBhmEYhmEUJGYEGYZhGIZRkJgRZBiGYRhGQVKcbQHSwR577KHV1dXZFsPIIC+88MJ7qjou23JkA9Pv4FPI+l1UVKSlpaXZFsPIIJs3b1ZVzUmnSyCMoOrqalauXJltMYwMIiId2ZYhW5h+B59C1u/S0lJ6enqyLYaRQURkS7Zl6I+ctMwMwzAMwzAyTcEYQc8//zyHH344Rx55JLNmzWL79u3ZFskw0obptxFUTLeNTFIwRtCECRN47LHHePLJJ6murua+++7LtkhGFOFwmMbGRsLhcLZFyTtMv42gYrptZJKCMYL23ntvvOC74cOHU1RUMF89b2hubqalpYXm5uZsizJgRGRTVNkpIv/t7qsWEY3a/+N0XNf02wgqpttGJhlybcrWS8Kjo6ODFStWcOKJJ6Z87AcffMBXvvIVysrKqKqqYuHChf22bWtr45hjjmHMmDHsv//+LF26tNf+N998k5kzZ7L77ruz1157ccEFF7Bjx46UZQoSTU1NNDQ00NTUlG1RBoyqjvIKsBewBfhDVLMKX7ur0nl9028jk2Tz+Z0rup1ov5FnqGrCAlQADcBtwBPA80ArcBVQm8w5+jnvKGATcKS7XQ0oUJzKeQ455BBNxI4dO7Srq0uPOOIIDYfDCdvH4vTTT9dTTz1VN27cqE899ZSWl5fryy+/3Kfd9u3btaamRufNm6c7duzQRx99VEeOHKmvvvpqpM2MGTP07LPP1i1btujbb7+tBxxwgP7qV78akFyFALBSB6hn2SrA2cDrgKjpt+l3HPJUv9Py/B45cmTc3yaXdDsZ3Tf6AvRoDuhsrJJIySuBm3F6s68BdwHzgLnA9cCTQA/wCnBayhfP4Evi5ptv1uOOO04bGhq0oqJCf/7zn+uMGTP0kUceSe2/57Jp0yYtKSnppexf+9rX9NJLL+3T9u9//7uWlZXprl27InV1dXU6e/bsyPbkyZO1tbU1sv39739fzz333AHJVgjk6UviMWCOb9vT77eAfwH/A+yR6Dym38EnT/U7Lc/vaCMol3U7Gd03+pLLRlCi4bC/ARuAaaq6n6rOUtWLVXW2qp6vqkcCe7geoe+JyPcTnC+as4Hb3R/JT4eI/EtE/kdE9kjxnACsWrWKZ599li9/+cu8//77jB8/nueee46rrrqKo48+mnvuuSfS9oQTTqCioiJmOeGEEwBYs2YNxcXFTJo0KXLcgQceyOrVq5OSR1V5+eWXI9sXXXQRd999N5s3b+att97if//3f5k+ffpAvqqRg4hIFXAUjvfU4z3gM0AVcAgwGrizn+PPFZGVIrJy/fr1ffabfhs5QEae37mu26nuN3KceBYSMC4ViyqV9jgvgp3ARF/dKGAaThLHPYFFwEP9HH8usBJYuc8++/SxPI888kidM2dOPOM0JZ588kndc889e9XddNNNetRRR/Vpu23bNp04caL+7Gc/023btulDDz2kJSUlevzxx0favPLKK3rwwQfrsGHDFNCzzz67V+/C6A151lMGZgNPJGizF07PeXS8drE8QabfwSIP9Tttz+/hw4f3+i1yWbeT0X2jL+SrJ0hV+3ZB09f+LOBpVX3Dd/wmVV2pqjtU9d/ABcDxIjI6xrVuUtVpqjpt3Li+2eZfeuklTjnllFTEj8uoUaPo7u7uVdfd3c3o0X1Eo6SkhHvvvZfW1lb22msv5s2bx6mnnsonPvEJAHbt2sX06dOpr6+np6eH9957jw0bNnDppZemTV4j63yd3l6gWHg96JQnKJh+G1kmbc/v4uLeCxfksm4n2m/kIclaS8CpwPG+7ctx4hoeAvZO1foC1gANCdrsifOiGBOvXXRP+c0339SysjLduXNnv5apn+nTp2tZWVnMMn36dFX9aFx5zZo1kePOOuusmOPKsTj88MN1wYIFqqq6fv16BbSzszOyf+nSpTplypSkzlWIkEc9ZaAWJ1ZudFT9Z4FP4hg9HwPuAf6Y6Hym38Enn/Rb0/z89scE5bpuD2S/oTntCUpF6V/xjCDgYGArcAlO8OfClC6a4ZfEfffdp4cddliK/6bEnHbaaXr66afrpk2b9Omnn+53hoGq6qpVq3TLli3a09Ojzc3NWl1drVu3bo3snzhxol577bW6fft23bBhg5588sk6a9astMscFPLpJQHcCNwRo34W8Iar+28DtwN7JTqf6XfwyTP9Tuvz228E5YNuJ9pv9CWXjaBU3PBVwKvu568A96rqz4HvAcemcB5wAuqWqOrGqPp9geXARuBl4EOcF0dKrFq1iqlTp6Z6WEKuv/56tmzZwvjx45k1axY33HADU6ZMieyfMWMG11xzDQB33HEHe++9N+PHj+fRRx/l4YcfZsSIEZG2S5YsYfny5YwbN47999+fkpISfvnLX6Zd5lwk6JmhVfU8VT0rRv1dqjpRVctUdW9V/bqqvpPq+U2/jSyTsed3Puh2ov1GfuFNbUzcUOR94ChVfVlE/gS0qOrNIjIRWK2qIzMpaDymTZumtsp2/tDY2EhLSwsNDQ3ccsstSR0jIi+o6rQMi5aTmH4Hn0LW77KyMrVV5IONiGxW1bJsyxGL4sRNIjwFzBORp3FmAHzVrZ8ErE23YEZw8TJC53NmaMMwDCP/ScUIugC4Acf4+ZaqrnPrZ+AERxtGUkyePDlpD5BhGIZhZIqkjSBV/RfQZ9EWVb0orRIZhmEYhmEMAbYcr2EYhmEYBUlcT5CI7OKjhG5xUdVhaZHIMAzDMAxjCEg0HHYqHxlBewJXAkuBP7t1hwMnAz/JiHSGYRiGYRgZIq4RpKqLvM8icj/wQ1X9ra9Ji4g8j2MIXZ8ZEQ3DMAzDMNJPKjFBxwB/jFH/R+DotEhjGIYRYIKeKNQw8o1UjKD3+Cg3kJ+vAikttGoYhlGINDc309LSQnNzc7ZFMYy0kq8Gfip5gi4H/kdEvsBHMUGHAccBjekWzDAMI2hYolAjqHgGPpBXeeBSyRN0u4i8ClwInORWtwGfU9XnMiGcYRhGkLBEoUZQSdXAF5HDcCZXrVbVFVH7fqCqP023jLFIxROEa+ycmSFZDMMwDMPIQ1Ix8EXkLOBXOMtxNYnIX4HTVHWT2+RHwJAYQSknSxSRShGZKiIH+0smhDMMw8hH8jU+wjCGiB8C01X1y8B+ODHHfxSRCne/DJUgSRtBInKQiKzGWSz1RWClr/wlM+IZhmHkHxYAbeQ7GTbkP66qzwOo6hZVPRt4HHhSRMaTZJLmdJDKcNhNOAbQOcA6hlBIwzCMfMICoI18J8OBzv8WkRpVbfcqVLVJRDYDTwMl6b5gf6QyHPZ/gQtV9U+q+qaqdvhLpgQ0MkMyVr659A1jYHjxEZMnT862KIYxIJqammhoaMiUIX8fcEZ0par+BPgfYEQmLhqLVDxBfwf2AtZkSBZjCFi2bBkXXXQR1dXVPPzww0BvKz8cDtPc3ExTU1OkJ9DZ2UlFRQVNTU32UDcMwygAMjmTUVX7taxU9Vrg2oxcOAapGEE/An4uIrNxDKLt/p2q+kE6BTMyw0UXXUR7ezvbtm2LaeV7hk9XVxcA9fX1iEhe5n8wDMMwjHikYgQ94v5dQe94IHG3bRX5PGD+/PlcdNFFzJ8/n5kzZ/bZ7xlFnZ2dLFmyJGIojRkzxuIbDMMwjECRihH0hXRdVEQex8k2vcOtektVP+nuOwPHFbYH8DDQYF6m9DFz5syYxo+H5wINh8O9hsDMA5Q8pt9GkDH9NoJE0oHRqvpEvDKAa1+gqqPc4t1AU4AbgbOAPYHN2Or0aSXZgGgvLshigAaM6bcRZEy/jUCQUsZocJIlAvsAw/31qvpkGuQ5E3jAO5eI/BhoE5HRqroxDecveJKZ9piva8DkAabfRpAx/TZSQkSWArcBraq6PVH7fs6xO3Ai8HHgLeDBVLyPqSRLrHTdoP8CnsFJbPRHX0mVa0XkPRF5RkSOduumAKu8Bqr6GrANmDSA8xcUyU5nT2baY4anRhYKpt9GkDH9NtLBUziLs78jIjeISG0qB4vI4cBrwLeATwPnAf9w65MilTxB84GdOPmCNgNHAKfgLKI6PYXzAFwK7Itjud0EPCAi+wGjgK6otl3A6OgTiMi5IrJSRFauX78+xcsHj1gZav2GkfcZSJi/xHKcDBrT7xwmHA7z1a9+lVAoZDmwBkZa9XvHjh3Ru40UyOd8bqr6C1U9GDgS6ATuEpF2Ebnc1alEzAfOV9VaVZ2lqp8Dvg1cl4oQSRXg38A093M3MMn9/CXg2WTP08+5lwPfxUmgdEnUvo3AIfGOP+SQQ7TQaWtr04aGBm1ra4vU1dfXK6D19fXa0NCggDY0NPQ5prW1tc+xuQawUgehY9kspt+5hXcvRN8P2aSQ9XvkyJHp+yELgOhnfaxne64B9GhyunQEjjdxJ44B/QhwYJz2G4CiqLphwIZkrqeqKcUEleIscgbwATAeJ3HiKzhuqMGgOFPtVwMHepUisi9O5khL0JiAWDO4RCTyN1Yaf8979NRTT9He7mQvtxigjGD6nUM0NTXR1dWFqkbuB5sMMChMv4eQ6JhNT587OzsJh8N5p78i8kngazgZpLcBdwAnAOuB84F7gYn9HN4OnA4s9NWdgjNElhzJWkvA8zirvuIK9TugCvgvoD2F81QAXwR2wwnMPhPowRk3noLjZToCKHOvcXeic1pPOTbRPYb+ts0TlNZesel3HpLt3nQh67d5glIjltc/2/qbCPrxBOEswP4e8Bvgs/20eSNWvbuvFscp8yxwD/Ccu13b3zF9zpF0Q0fZv+F+Phh4F8dltRk4JYXzjMNZdX4jzhjgs0Cdb/8ZwD/dG+s+YGyic9pLIjFtbW1aU1OT0zdKPPLoJWH6nYfEerEMJYWs32YEDZ5s628iYhlBON7D7wPDo/elUoDdcTxJl7h/E+qcvyQ9HKaqd/o+vygi1cBk4J+q+l5/x8U4z3rgM3H2L6S3a8tIA83NzbS3t1NTU2OzvjKI6Xf28Q9tAUkNc1lC0OQw/TbShaqqiFwB/GKQ59mA43UcEKnMDou+8GYc1+emgZ7DSA/e7IBly5b1O0sgFApRU1PD/Pnz827M2DBSwYuZmD17NieddFKfWZOGkQ+kMuvL0/mTTjop32aJ/ZVBpFAQkYkislBEXhGRf/pLsudI2hMkItcAr6rqbeJE3K4AjgW6RGS6qj6X+lcwBkI4HOayyy5DRJg7dy6zZ89m8eLFLF++nHXr1tHV1cWiRYsibZubm+nq6qK9vZ3FixfHXTbDMPKZcDhMZ2cnoVAIVY3r/bRgaCOXSSVpbVNTU2SCS3Nzcz55NR8HlovIrcBafOuSqmpLEscvxAmCvhgnNCdlUpkddiZwmvt5BjAVZ/2YM4Gfksa1xYz+CYfDnHTSSZHZXGPGjPHGRSP4t70bqaqqilAoZENhRqBpbm7utfCvf/07D8/48RYJBpsVaeQesWb0xiIcDjN79myqq6v59Kc/nW/P+M8BbwBHRdUrkIwRNAX4nKruGrAEKQQfbQU+4X7+NfAb9/P+QOdgApsGW4IaOBpvFkBZWZnW1dVpW1tb3Fle+R4Q7UGeBI5mogRVv1MhUeBnonsgVl6VUCiUM8GkhazfFhjdG7++JhPw7M97VV9fP4SSJg9J5glKtQAPkiAPVcJzpHCxt1yLC5y8D/Xu58lAVya+YLIlqC+JWNMeW1tbtby8PCWjJtdnDiRDIb8kgqrfqZBoCrC3v7y8XFtbW+Mem4v3QyHrtxlBvfHrq/fZS3gbS2fb2tq0qqoqYtjnIskYQTizxYq8Eqfdlb7ya5x8QjdF1V+Z6HpeSWU4bDGwUETWAGOBh9z6qcA/UjiPkSSx3KGLFy+mu7ubqqoq1q5dSygUorGxkcWLFxMKhVi8eHEf97/NfDHynURDA01NTSxatIju7m5mzZrFc889F7kHoo+1+8HIZUKhEE899RShUIh9992Xzs5OXnjhBTo6Oujq6mLMmDG9nvGTJ09m+fLlvWZE5gsi8nEcQ+ZInBxUfob1c9iEqO0HgZIY9cmRrLWEEz90MfAr4NI1X90AACAASURBVCBf/X8C/5HseTJRCqmn3NraqjU1NVpXVxdxgVZWVioQiGGv/qCAe8qFpN+DYSBe0lyhkPXbPEG9ifZcets1NTWRpZDyUL/7S5b4AE6Sw6k4y2QcCCwFzonVPhMl6zdAOkrQXhLR7nrP8FmwYEHE0Kmvr9fS0lIFdMSIERoKhfIi8/NAKeSXRND0OxmSGbJqa2vTuro6raqqigyB5eJQVzIUsn4XuhGUTGb/UCik9fX1efuMj2MEvQ+UuZ873b9jgXCs9skU4FPAH5Jun8KJD45XBipwOkrQXhKhUCgyvtvW1hbp3ZaVlSmgVVVVkReA5w3Kt55BqhTySyJo+p0MycTx+ANCa2pqsiVqWihk/S50I8h73nvPdT+e3uerB8gjjhH0LjDC/fwmTkbyEcDGWO19x40ErnI9Sb8AyoF9XS/SJtyJW8mUVJIlrsRJl77SV/7iK0aacP/JdHd3M336dLq7uxk2bBjDhjlDpJMmTWLy5Mlcd911hEIh6urqIovnGUY+0V9CuKamJkKhUESvvVQP/qSHTU1N1NXVUVVVxcUXX5x0YjnDGGqi9dy/7T3vOzo6+iT19PReRCJpHwLGc4CXuO4hnKGxJTj2RTx+A5yIs4D7cTgxy0/gLOJbrarfSVqCZK0lnMVS/WV/nNVa/wbMSPY8mShB6CnHmhbp9/T4S2VlZV4tnpcOKOCechD0uz/8uuvpvTfs6+m/t88bEog1FJDv90Ah63cheIL6i/Px0jqMGzdOS0tLdcGCBb2Oy9fh3Wjo3xNUgbvWF1AKzAZ+Buwdq73vuHXAePfzJ4BdwBHxjun3XAM5KEqY44FnBnuewZR8fEn0l7ukpqYmUue5QHHjfvyGkP9hH5QbJR6F/JLIR/1OFr/u+qe5e8Z+TU1NJN4nnqGT7/dAIet3EI2gRHE+XpynF+NDQIZ1+yOWEYQz++s23OGwVArQHW87pXMN9EDfxWv6s/KGquTjSyL6gd7a2hqJ+amrq4v0ELz8D3V1dVpXV6fjxo3TysrKPrlQgk4hvyTyUb8Hguftqa2t1aqqqognyEtqmK9BoclQyPodRCMokcHun8nb2tqqlZWVgX6ux/EEvQ2UxNoXr+AskfEF4Bi3dEdtH5P0uVK46Nio8jHgAGAR8GKqXyKdJR9fEv4ZXw0NDVpbW9vH61NZWRkZGojuMeSr23+gFPJLIh/1O1Wih4Dr6+t1wYIFWl5eHrk3/B6jeBmh85FC1u8gGkGxdDI6vYln+HgpToL8TI9jBF0CXJOqIYQTRP1GnPJ60udK4aK7gJ1RZRfQARyWyhdId8m3l4S/J+C5/vsrnnfI7xnylssoJAr5JZFv+j0QPAPfS/vgTX33Xhae/sfKhZXv8UCqha3fQTSCYj2rPX3ebbfdtKqqKmL8RIdBBJE4RtBaYDvOslxrgX96JVb7TJRUMkZHL5C6Cydd9T9UdUcK5ylYold0Lykpobu7m+LiYnbs+OgnHDFiBB9++CEAPT09kVWwm5ub6ejoYPjw4dn6CoaRdsLhMF1dXVRWVrJu3ToARo8ezaRJk+jo6KC6upquri7OP/98Ojo6+qwKn+xCk4YxVFx00UV0dHTQ0dHBhRdeyIQJE/j4xz9OR0cHW7dupaOjg7q6OkpKSpg0aRLXXXddryz/BcTX0nUiEblfVU9K+cChsrYyWfKlp+z1WKuqqiK9guLiYh01alSkR1BcXBxxl5aWlvbqSQRlMdSBQAH3lPNFv5PFP9PLG/Ly7gvvr3+WpJdHhQD3mAtZv4PkCfLPcPS8msOHD1cgEuvm1/FCgSGIGwY+GNBxKV5kT5zFyRYBfwCuAPbM9JdLVHLxJRFryrs/0Hnq1KkqIr2GvoqLi3XOnDlaX18faRdt7AQh/mEgFPJLIhf1ezD4Uz+UlJRE/i5YsKBPBmjV2EMLQaOQ9TvfjSC/Ue/N6I1e2sjfqe0vzUOQ6c8IAoa7NkU70OP+vQrYLVb7eCXjRhDwOWAjzmKpd7jlHzhR2YcP5OLpKrn4kvCv/ut5b8rKynTs2LEKaFFRUZ8bxB/8HG/V4EKkkF8SuajfydBfR8DrIUffA1OnTo0YRf6pwkGI+UlEIet3PhtB/hXcAR03blxku7S0VMeOHdvL0x90Pe6POEbQLcDTwAzg/7p/nwJaYrWPOvbyqLLZv53o+Mh5km4If8ZZrr7IV1fk1v0phfOMcL94h2tURZItAtWuomzylR8nOmeuvSS8nkFdXZ2OHz8+buCzv3izwOrr6yNLZhgO+fKSKAT9ThZ/7itvhpdnAAF97g3PMyoivZLGFYL3s5D1O1+NoGgDyCuVlZWRCS3RHd1Cfa7HMYLeByqi6saShFcH+ElU2eLfTnR85DxJN3Qu8MkY9ZOBLSmcpwyY494wRcAJ7s1U7buJipM9n+bgSyI66Vt/xf9C8D43NDQURM83VfLoJRF4/U4WfwxbcXFxRNerqqoiC/76Y+P8hlCh6X4h63e+GkH+lCXRyWzr6uoihpC3r66uLtsiZ404RtBqoDKq7uPA6ljt45VkDKdYJZXZYV3ARODVqPqJQGeyJ1HVHvcm8nhQRN4ADgFeSEGenCUUCrFo0SK6u7sZO3YsH3zwQZ82tbW1lJWVsXr1atatW8fnP/95JkyYEHPWi5E/FIJ+J0M4HOayyy6jurqazs5O1q9fT1FREWPGjGHYsGHce++9vPbaa/zwhz9k3rx5XHzxxcybN4/29vY+s7+M3KGQ9dub3RsKhVi8eDGHHnooy5cvZ+vWrYDjUNi2bRulpaWEQiFGjx4dWfNr8eLFptMuInKMb/MOYLmI/DfwL2AC8B3g9oGcekACpWBlzQfeAs7EMXwm4kxve+v/s3fn8VGV9+LHP1+SsCRhVTYXgrRgvFihiksREbWhglYto1bA1ja2tLWWH2pzbYsV6lK1EaX2ugtXaxW9NbgVioILYt0u9gqiGUlVIm4VEQgEkAS+vz/OmXAyTDLnJDOZ7ft+vc4rM2d58szMd2a+85znPA9wU1syMLfc/jhjBJSy95fER+4T8t/A/vHKSOUvZW9fh+jh//Pz8zUvL6/V1qDI6NC52EQaBBnySzl6yfT4bk1rp6lamvcueokeOdc7TUYuyeX4TveWoOirFCMtO9GtP94lchYg11o0W4KnJYjWBzmMLL4HO/SU+2TQY1SDnQ7rDPwR+JK9gyXuBG4GOrfpn0MBsAy4071fDIwC8t031yPAUy0cOw13NvtBgwYl7tUKKHqQt9GjR2tZWdk+nT5bWiZNmpSyumeSTPySyIb4bk30aVtvX7iuXbu2GPPFxcXaqVMn7dOnT9OI6d4fD7n4xZHL8d25c+dEP50J5b1YxU8fz8h3Qa5dBt8aUjy1VmtLWwK/EPiauxS2+R8755MfAhbTwpDZwAA3sLq3VlYqfykvWrQobt+flpbS0lJ7k/iUaV8S2RLfrWlp1vdYS2FhYdMv58iVj9EJTy50gG5JLsd3prQE3XHHHc1iOnqIk8g679QvxpHOSVCQPkEAqOp24M2gx3mJiOBcYdAfmKiqDS39O/dvp/b8v0SKPi+8fv166urqApVRUFBAQ0MDu3fvTlItTSplcny3JBL3FRUVlJaWNt0/5phjuPjii2lsbKSmpqbF47dv304oFKJnz54t9nsrLS1l3rx5SX0cpv2yMb69vLE9Z84c5s6dSygU4qyzzmq2X15eHqpKcXExW7ZsQUSYNWsW3/3ud5veK2ZfIlKtqoe5t9ezN06aUdVBHVIhv9kS0BW4HHga57LI1d4lSOYF3AG8AhRHrT8WOBTnTbMf8DDwXLzyOvKXcvSVX62dF4619O7du8U5kEzLyKBfypkc3y2J9IcIhUKq6r/PT2Q57LDDcnKQOL9yOb7TqSXIe9m73y4N2KmvuGjeJ2iM5/aJLS3ajhgNsgRpCboN+A7OSNEv0UL2Fo+IlAA/welb9KnzowLcdXtwZpTthzMI41Jgclv+T7KEQiEWLFjQ1PoTmeOrNQcddBAbNmygR48e3HvvvUycOJEhQ4bYr4UslOnx3RL3A4uVK1dy55138uyzz/o6rqSkhNraWhobG1m4cCG9evWy1p4Mlq3x7Z3Xsba2FoA9e/a0eky/fv3YunUrO3bsYNSoUbk691dgqvqi5/by6O0ikocz1s8+25JVIb/Z/xfANzsqOwuyJPuXcnV1ddNw/iNHjgzc9yfS6oO1/LQZGfRLOdFLKlqCvP1zFi1apP369Wu60tE75k9rS6TPW/TVk2ZfuRzfqWwJisRnpGWztLTUV2zn5+c3Gwnd4rp1BOgThDMg526/+7d3CdIStB1nqvuc4f11sHTpUoCmXwmtifT5AejWrRuXXXYZS5cuRVWt5cekvXA4zBlnnEFNTQ3PPPMMn3/+OfX19U3bGxsbfZXzH//xH836+UycODEp9TUmqMhn+/r161m6dCl9+vRpWt+a3r17s2vXLubMmdPU8mMtm0nRtjF/2iBIEvQH4FIR+amb2WW9yspK5s+fH7iZ8+ijj+all16iqKiI+vp6XnvtNR555JEk1dKYxJo+fTo1NTXk5eX5Svq9SktLGT58OKrKtddem6QaGtM+kc/2oqIigJgD2kbr0aMHL730kp326hgdlmO0mgSJyBNRq8YCp4rI20CzKwJU9YwE1y3lQqEQK1asYO3atb6Pifzy9V5BZq0/Jh15r/h67733mDFjBlOnTmXZsmUA7N69m06dOsXtGxGRn5/PjBkz+MlPfpLMahvTZosXL2bGjBmMHz+evLy8Zi2crenWrRsLFiywBChBokaNjta5wyoCrfcJwhnx09fSUefvYi2J7DPh7cMQmfE9yJLL88MkEzncZyJZfYIiVzqWlZU1zd7eliU/P7/pasnI7O/WVyKYXI7vZPcJ8sait3+mn6VTp05aVlZmcdxORPUJwseo0dpB8ddqS5Cq/tBvMpUtIs2kf/7zn333ffDq3r17EmplTOJFWiiffPLJpj5sbXHMMccwc+ZMZsyYwdy5c4G97yOwPhMmtSKxuHnzZnbs2BHo2FNOOYWnn346STXLXap6SKrrEOG7T5CIPAtMUtXNUet7AI+pamvNWxlh8eLFPPnkk4D/zp8AnTt3ZteuXZSUlFg/CJMRIhOcigi7du1qczmjR49m3rx5lJaWNuv4HEmw7FSwSbVjjjmGRx55hDfeeIMPP/zQ93H9+vXjlltuSWLNTDoI0jF6HLHP1XUFTkhIbVIk0jdiyZIlbNiwIdCxZWVlzJgxo6nvj50zNuksHA5zxRVX8PLLL/Pxxx8D4BnrxZfi4mJUlTlz5rTY/8dGfzbpIBwOM336dHbt2uV7ZP/Ro0czcOBArrnmGvs8zwFxkyAROdJz9wgR8XajzwO+hTNrcMbyNt0HdfDBBzNx4kS7/NdkhJkzZ7Jw4cJm69xz9L54W36MSTfeaY3mzp3L8uXLfbd0durUidtuu8069ucYPy1BK9nbUSzWydEdwC8SWamOFgqFAiVBeXl53Hrrrbz22mvW3G8yQuTL4f3332/T8ZF5kWbNmpXgmhmTOFdccQVVVVVUVVWxZcuWQMd+5zvfsQQoB/lJgg7BGbjoPeAYwHu+aBfwmapm7Eyg4XCYUCgU6JiTTz6Zn/zkJ/aGMWkvelC4tujduzebNm3igw8+SHDtjEmccDjMypUrAQIlQIsWLbKhTHJY3CRIVSOjpWXUTMCtiXwxvPvuuyxfHmx6kr59+1pnOZMRvCM/5+cH6f6319ChQ5k7d659SZi0Fg6HOe644wIlPwUFBfzpT3+y7gw5Lt5giWPUM9lZnH2LgUNU9c2E1CyJpk+fHvhX8c9//nN27NhhnZ9NxqisrGxKgIJc7di1a9emUwqReLcvCZPOJk+eHCgB+vrXv84///nPJNbIZIp4Pw/vEZGPgHuARaq6T/d6ETkCOB/4PvBLIO2ToJdeesn3vgUFBTz22GP2JWAySjgcbvqQDzre1QknnGC/jk3GOP3003njjTd879+lSxcefPDBJNbIZJJ4p7gOBxbiTGv/hYi8IyLPicjfReQVEdkEvAocCJysqn9Jcn3bJRwOM378eN9DpY8cOZLVq1fbl4HJKOFwmFNOOSXQFwNAz549KSsrs9O9Ju0tXryYAw88kIKCAhYtWuTrmP79+1NSUsLChQutNd80iTdidCNwK3CriIwCxgAlQDfgdaASeE5V488+lwbGjBnDxo0bfe07adIkqqqqklwjYxJr8eLFfPvb3/Y931fEokWLLNk3GSMUCrFz507f+9vnuWmJ796SqroS53L5jDRs2DBfCdCUKVPo2rWrdQI1GWfx4sWcdtppgY8bPXq0JUAmI9x555387Gc/8z22VUFBAePGjbOR/E2L2nbJSIY5/fTTqampibtf//79eeCBBzqgRsYkVtAEKC8vj2OPPbZpZFxj0l04HOanP/2pr31FhL/97W+W3Ju4sj4Jmjp1qq9zxkOGDPF9btmYdBIOhwMlQDbqs8kk4XCYsrIy3/N+TZo0iWuvvdbi2/iSlmP/iEgfEXlUROpFpFZEprSlnHA47OsqgNNOO413333X3jQm4yxevJjDDjvM1775+fksWrSIf/zjHxbrJmOMGTPGdwJUVlZGVVWVxbfxLS2TIJzO2LuA/sBU4HYRGR60kHhfDiNHjqS6upq//e1vbaqkMUElKsEHp3+E3xag0aNH09DQYKcHTFIlMr7d8nz15Tz66KMpLy+3KxtNYGl3OkxEioAQcLiqbgNeFJEngO8BvwpQTqvb77jjDpv2wqSCN8EfCSwSkVWq+lbQgvz0j8jLy+PMM8+0jqGmoyQsvuN9hkfYlY2mPXy3BInIuSIy3nP/ShH5UESeEpGBCazTMKBRVdd61q0CArcEtURVLQEyHc6T4P9WVbe5o7FHEvygZcXdp6ysjDVr1tjpAdMhOjq+Z8+ejapaAmTaJcjpsNmRGyJyJPAb4BagAJiTwDoVA9EjU28BuntXiMg0EVkpIis3bNiAH/n5+b4vrTQmCZKe4INz6qu6upqnn37akh/TkRIS3/ESoPz8fKqrq5k1a1bwGhoTJUgSVAK8497+DvCYqv4BuBQ4JYF12gb0iFrXA9jqXaGqd6nqKFUd1bdv37iFTpkyhYaGhsTV0pjgfCX40LYkH5w57qzjs0mRNsV3kGld7rjjDhoaGiy+TcIE6RO0k73BfAow370dM8jbYS2QLyJDVTUyuM8IINA5ZWvxMWnIV4IPTpIP3AUwatSouME8e/Zs+2VsUq1N8V1UVBQ3vo8++mhee+21RNTRmGaCJEErgDki8iIwCjjbXT8MWJ+oCqlqvYgsBK4SkR/hdK47ExidqP9hTIokJMEHS/JNWrIfsCbjBDkddjFOr/+zgZ+q6sfu+gnAUwmu10U485N9BiwAftaWqwuMSSeqWo8zIfFVIlIkIsfjJPj3p7ZmxrSfxbfJRJINWbeIbABqY2zaH/i8g6uT7jL1OSlR1fidv9KciPTBOZVcBmwEfqWqrY7oafEdWCY+L7kc33uAHVGr8wH/nYVyR6Y+L91UNS3HJfSdBInIe8DRqroxan0v4J+qOiQJ9WsXEVmpqqNSXY90Ys9J9rDXMjZ7XjKfvYax2fOSeEEys8FAXoz1XYCDElIbY4wxxpgOErdjtIhM8tw9TUS2eO7n4Vwp9n6iK2aMMcYYk0x+rg57xP2rwLyobQ3AOuCyBNYpke5KdQXSkD0n2cNey9jsecl89hrGZs9LggXpE/Q+Tp+gTOtwaIwxxhizj6y4OswYY4wxJqhAs8iLyLE4fYD6EdWpWlWnJ7BexhhjjDFJFWQW+V8CLwM/wBnF+Wue5fBkVK6tRKSPiDwqIvUiUisiU1Jdp44mIs+LyE4R2eYu73i2TXGfl3oRecwd28NkCItvh8V49rHYdlhsd5wgl8j/P2C6qg5T1XGqepJnOTlZFWyjW3FGt+4PTAVuF5GEztSdIS5W1WJ3ORTAfR7uBL6H8/xsB25LYR1NcBbfe1mMZxeL7b0stjtAkCSoB7A4WRVJFBEpAkLAb1V1m6q+CDyBEzTG+WB5UlVfUNVtwG+BSSKSyElwTZJYfPtiMZ6BLLZ9sdhOsCBJ0ALg1GRVJIGGAY2qutazbhWQi78mrhORz0XkHyIyzl03HOf5AEBV38X55TUsBfUzwVl8N2cxnj0stpuz2O4AQTpGrwd+506KtxpnjKAmqnpTIivWDsVAXdS6LUCuZcqXA2/jvEHOA54UkZE4z8+WqH1z8fnJVBbfe1mMZxeL7b0stjtIkCToR8A2YLS7eCmQLknQNpxTd149gK0pqEvKqOqrnrv3ichkYCL2/GQ6e/1cFuNZx143l8V2x/GdBKnqIcmsSAKtBfJFZKiq1rjrRgBvpbBO6UABwXkeRkRWisgQnPnf1rZwnEkvFt8tsxjPbBbbLbPYTpKsHCxRRB7CCZof4VzOvxgYrao58WYSkV7AscByoBH4Ls5w618HCnCGOjgN+CfOlQb5qnpeamprgsr1+AaL8WxlsW2x3dF8twSJyC2tbU+zwRIvAuYDnwEbgZ/l0psI541yDVAK7AbCwFmRDoci8lPgAWA/YBnwwxTV07RNrsc3WIxnK4tti+0OFWTusOeiVhXgvEh5wP+l4VhBxhhjjDEtCtIn6KTodSLSFWdm+RWJrJQxxhhjTLK1u0+QO4LlElU9ODFVMsYYY4xJviCDJbZkf5yxC4wxxhhjMkaQjtGXRq8CBuIM453202kYY4wxxngF6Rj9ftSqPcAG4FngOlW1wZqMMcYYkzGycpwgY4wxxph42tQnSESK3Rl/jTHGGGMyUqAkSER+LiIf4EzYVicitSJyUXKqZowxxhiTPEE6Rv8G+DVwI/Ciu/oE4HoR6aGq1yehfsYYY4wxSRGkY/QHwOWquiBq/VTg96pakoT6GWOMMcYkRZDTYf2A/42x/jWgf2KqY4wxxhjTMYIkQWuBKTHWTwHeSUx1jDHGGGM6hu8+QcBs4H9EZCzwD3fd8cCJwDkJrpcxxhhjTFIFGidIRI4CLgEOc1dVA3NU9f+SUDdjjDHGmKSxwRKNMcYYk5N89wkSkXNE5MwY688UkbMTWy1jjDHGmOQK0jF6NrAzxvp6d5sxWUFELhaRlSLypYjcG7XtFBEJi8h2EXlOREo827qIyHwRqRORT6MnHW7tWGM6isW3MXsFSYKGEPsqsH+524zJFh8D1wDzvStFZH9gIfBboA+wEnjYs8tsYChQApwE/KeInOrzWGM6isW3Ma4gV4dtwnkDrItaPwxI6Qzy+++/vw4ePDiVVTBJ9vrrr3+uqn074n+p6kIAERkFHOTZNAl4S1X/6m6fDXwuIqWqGgYuAH6gqpuATSJyN/ADYImPY1tk8Z39cjm+O3XqpN26dUvgIzTpZvv27aqqbZqrNNmCJEGPAzeLyCRVXQsgIocCNwGPJaNyfg0ePJiVK1emsgomyUSkNtV1AIYDqyJ3VLVeRN4FhovIv4GB3u3u7bPiHQu0+iVh8Z39cjm+u3XrRn19fWIegUlLIrIj1XVoSZDM7HKciVPfFpH1IrIeeAuoAyqSUTlj0kwxznvAawvQ3d1G1PbItnjH7kNEprn9NlZu2LChXZU2xqeUxHdjY2O7Km1Me/hOglS1TlWPByYAt7jLqcDxqlqXpPolzGuvvcY3vvENxo4dy+TJk2loaEh1lUzm2Qb0iFrXA+d08DbP/eht8Y7dh6repaqjVHVU377xz5JYfJsESEl85+e3fkLCYtskU+BzdKq6VFUr3WWZZshAQwcffDDPPvssL7zwAoMHD+bxxx9PdZXaJRwOc+GFFxIOt9rSbBLrLWBE5I6IFAFfwekLsQn4xLvdvf1WvGMTUbFMiG+L2bSXlvGdCbGdqew92YYkKFMNHDiQSOe7zp0706lTZj/0yspK5s+fT2VlZaqrknVEJF9EugJ5QJ6IdBWRfOBR4HARCbnbrwRWezp+/hm4QkR6i0gp8GPgXndbvGPbJRPi22I2PWRafGdCbGcqe08CqprwBadp1LvsBv7kbhsMaNT233qO7YJz6WYd8Clwabz/d9RRR6lf69at0+OOO0537drl+5iIjRs36llnnaWFhYU6aNAgfeCBB1rc9+2339aTTjpJe/TooV/5yld04cKFbS4rlurqai0vL9fq6urAjyMTASs1CbEaa8G5FFijltnutm/idPTcATwPDPYc543df0fHbmvHtrZkYnxPnTpVBwwYoN27d9ehQ4fq3XffnXMxG0Qux3dhYaGv5yhdYjti7dq12qVLF506dWrg+qSLjnpPAvXaQfEddOmIN1yxm+iMde9HkqD8Fva/DlgB9MaZo+xT4NTW/oefL4nGxkbdsmWLnnDCCRoOh+PuH8t5552n5557rm7dulVXrFihPXr00DVr1uyzX0NDgw4dOlTnzJmjjY2N+swzz2hhYaG+8847gcsyjo78kki3JRPje82aNbpz505VdT5o+/fvrytXrmxTvXJBLsd3vCQo3WI7oqysTMeMGZPRSVBHyfUk6ALgPfbOUxYvCfoYGO+5fzXwUGv/I9aXxD333KPf/OY3tby8XHv16qV/+MMfdMKECbps2TJ/r1qUbdu2aUFBQbM3w/nnn6+XX375Pvu++eabWlRUpHv27GlaV1ZWpldccUXgsowjl78kMi2+o4XDYR0wYIA+/PDDbapbLsjl+I5OgjIhthcsWKDnnHOOzpo1y5IgH9I5CQp0ctU9d3y2iFwuIr3cdV8RkT6tHHYB8Gf3ifCqFZEPReS/3dFGEZHexB6LYniQegKsWrWKV155hTPPPJONGzfSr18/Xn31Va6++mrGjRvHww/vHcz09NNPp1evXjGX008/HYC1a9eSn5/PsGHDmo4bMWIEb73lr9+fqrJmzZqElGVMOsd3xEUXXURhYSGlpaUMHDiQiRMnCzDJ2gAAIABJREFUJuCRm2yX7rFdV1fHlVdeyU033ZSgR2xSym+2BHwVeB/4HGgEhrjrbwTuaeGYEpz+QId41hUDo3AGauwPPAI85W47GKeVqKtn/zJgXYyyp+EMzb5y0KBB+2SeY8eO1dmzZwfMV1v2wgsvaP/+/Zutu+uuu/TEE0/cZ99du3bpIYccojfccIPu2rVLn3rqKS0oKNDx48cHLss4yOFfyrFagtI5vr0aGxt1xYoVevXVV7epL0euyOX4jm4JSvfYnj59ul5//fWqqtYS5BNZ0hI0F1jqJi7e0R+fwJlHJpbvAS+q6vuepGubqq5U1UZV/TdwMTBeRLoTfywKPOW0Oo7K6tWrOeecc/w/ujiKi4upq2s+HFJdXR3du+87FlhBQQGPPfYYixYtYsCAAcyZM4dzzz2Xgw46KHBZxsSSzvHtlZeXx5gxY/jwww+5/fbbE1Zfk73SObbfeOMNli1bxiWXXJKw+pnUCpIEjQZuVNXdUes/AA5o4ZjvA/fFKTdymqyTxh+Lwpfa2loaGhooLS31tf+ECRMoLi6OuUyYMAGAYcOG0djYSE1NTdNxq1atYvjw2GfqjjjiCJYvX87GjRt56qmneO+99zjmmGPaVJYxXuke37E0Njby7rvvBniUJhele2w///zzrFu3jkGDBjFgwABuvPFGqqqqOPLII9v5yE3K+G0yAr4Ahru3t7L3dNhY4NMY+48G6oHuUeuPBQ7FScD2w5lp+DnP9uuB5ThXh5XiJEWBrg57/PHH9bjjjvPfVufTd7/7XT3vvPN027Zt+uKLL7Z6RdeqVat0x44dWl9fr5WVlTp48OCmq2WClmU0p08XZFp8//vf/9YFCxbo1q1btbGxUZcsWaKFhYX6+OOPJ7zO2SKX49t7OizdY7u+vl4/+eSTpuWyyy7TUCikn332WcLrnE3IktNhTwOXevMnEekB/A5YFGP/C4CFqhp9KmsIzqzDW4E1wJfAZM/2WcC7QK2bDFWq6pIA9WTVqlWMHDkyyCG+3HbbbezYsYN+/foxefJkbr/99ma/JiZMmMDvf/97AO6//34GDhxIv379eOaZZ1i6dCldunTxXZYxLUn3+BYRbr/9dg466CB69+7NL3/5S+bOncsZZ5yR8Dqb7JLusV1YWMiAAQOaluLiYrp27YqfqW1Meopcth5/R5EDgOfcu0OA/8PpLP1vnDGAUjbL46hRo9Rm2c5uIvK6qo5KdT1SweI7++VyfBcVFanNIp/dRGS7qhaluh6xtD5znYeqfiwiI3FabY7EOZ11F/CAqu5o9WBjjDHGmDTjOwkCcJOd+e5ijDHGGJOxWk2CROT7fgtS1T+3vzrGGGOMMR0jXkvQrVH3OwMFwB73fiegAadzsyVBxhhjjEk6EcnDOSs1TVW/bGs5rV4dpqrdIwtwHrAaOAHo6i4nAG8AU9paAWOMMcaYINQZs3A8extl2iTIJfI3AtNV9R/qjPbcqKr/AGYAc9pTCWOMMcaYgG4GficiBW0tIEjH6ME4gx9G2w4MamsFjDHGGGPa4BfAAOBSEdnA3hkoUFVfeUmQJOhV4BYRmaqqHwGIyIE4mdgrAcoxxhhjjGmv89tbQJAk6ELgMWCdiHzkrjsQeAc4q70VMcYYY4zxS1WXt7eMIIMlvisiRwBlOHN6AVQDy9TvsNPGGGOMMQkgIl2AK3EGcd5PVXuKyHhgmKr+l58ygnSMxp0L7WlVvcVdlloCZIwxxphEEJF8EfE7IPPNwOHAVPb2B3oL+Jnf/xcoCRKR00TkBRH5XEQ2iMhyEZkYpAxjjDHGmBbk4UzA7sd3gCmq+jLupfJun+UD/f4z36fDRORHwG3AA8B97uoTgEdF5GeqalNpGGOMMaZVIvJsK5vzAhS1i6g8RkT6Ahv9FhCkY/TlwKVR59nmicjrwK+w+cSMMcYYE9+xwHXAJzG2FQBjfJbzV+A+EbkEQEQGAnOBh/xWJEgSNAhYEmP933EGUjTGGGOMiecNIKyqj0RvcDs73+aznN8ANwBvAoVADXA3cJXfigTpE/QBzpVh0cYDtQHKMcYYY0zumgt80cK2BuCHfgpR1V2qeomqFgP9ge6qegmeQRPjCTptxh9F5G4R+aG73IPTO9tagozJEeFwmAsvvJBwOJzqqhiTMex9s5eq/lVVY/YLUtU9qnpfrG3RRGSu57gNqqoi0hV4wm9dgowTdKeIfAZcBkxyV1cD56rq437LMcZktsrKSubPd7oAzps3L8W1MSYz2PsmKY4Qkd+p6iwAESkEngQ+9FtAkD5BqOqjwKOBqmiMySoVFRXN/hpj4rP3TVKcCSwTkS3AnTh9lN8BpvktwPfpMBE5UURObGH92BjrnxeRnSKyzV3e8WybIiK1IlIvIo+JSB/Ptj4i8qi7rVZEpvitozEm+UpLS5k3bx6lpaXxdzbGAPa+SQZV3QpMwOlD9AawSlV/HGQQ5yAtQTcTu8d1D2A2cFSMbRer6j3eFSIyHCdjOw34J3AXTk/w89xdbsW59r8/MBJYJCKrVPWtAHU1xhhjTJYRkVh5yGs4OcWmyHZVvdJPeUGSoEOBVTHWr3G3+TUVeFJVXwAQkd8C1SLSHWfExxBwuKpuA14UkSeA7+GMRWSMMcaYLCAij+IMvrxIVRt8HnZwC+v/3sq2FgVJgnYAA4H3o9YfiNNyE8t1InI9zjm6mar6PDAceCmygzsx6y5gGE4S1Kiqaz1lrAL2OQ1njDHGmIy2AmcC1Hki8j/A/ar6UmsHqKqvy+f9CnKJ/FPADSLSO7LC7ctznbst2uXAEJwk6S7gSRH5ClAMbInadwvQ3d1W18K2ZkRkmoisFJGVGzZsCPAwjDHGGJNqqnqTqh4JjAU2AwtEpEZErnTzhbhEpKeIHCMiJ3sXv3UI0hL0S+AFYJ2IrHbXHQF8Bnw3emdVfdVz9z4RmQxMBLbh9CPy6gFsxWkJamlbdPl34SRXjBo1ymayN8YYYzKQ2+f31yKyGPgvYBZwmYj8L3CZqsbqioOI/ACnH/E2YLu3SJxGmLh8twSp6ifACJxkaLW7XAaMUNWP/RQBCM409yMiK0VkCNAFWOsu+SIy1HPcCPeYrGIDZxljjMl1InKoiFwtIu/iNGw8DAzGuThqMfBYK4dfC5ytqv1V9RDP4isBAp8tQSJSAPwF+I2q3u1j/144E6QtBxpxWorGAv8PZ3K0l0XkBJyrw64CFrqXuiEiC4Gr3FnrR+KMAzDa7wPKFDZwljHGmFwmIitxEp6HgSlRZ5AAbhKRX7RSRD7wdHvq4CsJUtUGERkP/NpnuQXANUApsBsIA2dFOjyLyE+BB4D9gGU0nyfkIpwZ6T8DNgI/y8bL423gLGOMMblKRARntvdbVLWli6tQ1UNaKeYG4AoRuVpV97SpHn7HFBKReUC1qqbdPGGjRo3SlStXproaJolE5HVVHZXqeqSCxXf2y+X4Lioq0vr6+lRXwySRiGxX1aIY6+txJj1tWwIjsh4YgHOF+kbvNlUd5KeMIB2jP8DJuE4AVgLNolZVbwpQljEZTUSeB47DOd0L8JGqHupum4Jz1eT+wFKgXFW/cLf1AeYB44HPgV+r6oMdW3tjWmaxbTrQ/+EMj9PWzrHnt7cCQZKgHwCbcK4IOyJqmwKWBJlcYyOim2xlsW06wvPAEhG5F1iPk0sAoKrz4x2sqsvbW4Egs8i3dl7OGOOwEdFNtrLYNol2PM4AzNEDIitO3+B9iMhMVb3WvR1rCg2ngCRMm2GMac5GRDfZymLbJJ2qntSGww7y3A48TUa0QEmQiAwDzgYGAZ2921S1vL2VMSaDXA68jdP8fx7OiOgjaX1E9N0EGBEdmAYwaJCv/n3GJEpSYxuax3fnzp1j7WJyjHu1mETut9RZWlV/5rnd7ik0fCdBInIaUIXTkeko4H+Br+AMdLiivRUxJpPYiOgmWyU7tt3/0RTfRUVFFt85SkQOxBkheizQK2pzXgvH+BoIUVXf87NfkJagq4Dfqep1IrIV51zvx8D9wMsByjEmG/kZEX0P7ojoqlrj7pKVI6KbrGKxbZLlDpzpLk7BGVx5LDAbZ6TolvyLvTHZEqWFJCpakAlUD8UZ1RGgAShU1Z04ydGMAOUYk9FEpJeIfEtEuopIvohMxXnzLsEZBPTbInKCiBThGRFdVeuByIjoRSJyPM6I6Pen6rEY42WxbTrYaJxhFt4A1J0j7EKcKbliUtVOqprn/m1p8ZUAQbAkaCvQ1b39CfBV93Y+0DvmEVnK5v3KeZER0TfgjIfyC9wR0d3LgSMjon+G0yfiIs+xFwHd3G0LyNIR0U3Gstg2HWk3e8ej2iwifXHGIDzQz8EicqaI+E54YglyOuxVYAxOh7lFwBwRGQF8hxw7HTZz5kwWLlzI5s2bqaqqSnV1TAdT1Q3A0a1sfxCIOUicO7DcWUmqmjHtYrFtOtirOP3NHgWewjnbtANnQGY/rgLuEZGHgftjzD0WV5Ak6FKcqwPAOWfXHWdciLXutpzhdGLf+9cYY4wxgX2PvWekZuCcBusOzPVzsKqOcBtjzgeq3Gk47gf+oqrr/JTh+3SYqr6nqqvd29tV9WeqeoSqnq2qH/gtJxtcc801lJeXc8011zSts1NkxhhjjD/uaaw/4k7Bpao7VPUaVb1cVT/xW46qrlLVCpwxg34OnAO8KyIviMhUEWk1zwk8WKKInAz8h3v3bVV9NmgZ2aiyspL5850BLufNm5fi2hhjjDHpS1V3i8h4nCsL20VEvoLTGnS+W96VOPOdXoxzxmpSS8cGGSfoEJxxgo7AuTQe4AAReRMI+b0mPxOFw2EqKyupqKigtLQ0ZsJTUVHR7K8xxhhjWnUz8DsRmaWqDUEPFpGf45xSG4rTn+h7qvqKZ3sVTkf9FgW5OmwezhViQ1R1kDtN/RBgM3BPq0dmuEjSU1lZCTiJTnl5ebOEp7S0lHnz5lFaWpqqahrji/fUrZ3GNSa+cDhMKBTi7LPPtvdKYv0CqAC2ish6Efkgsvg8fgIwBzhAVS/yJkDgdN2hlVYgCHY67BvAcd7+P6r6gYhcQpZfHRbdyhNJeCJfIJEWImMygbclE2jWqhnd6mmMcd4zCxcuBGD16tU88cQT9v5IjPPbc7Cqnu5jn6db2x4kCfoAZwyIaF2B9QHKyTiRpCea9QMymSjWqdvIbYtpY/ZVUVHB5s2bef3116mpqaGystLeHwmgqsuDHiMi9+OMCB2v7O/7KS9IEnQZcIuITMeZN0yBY3AuZWtxdMdsZv2ATCaKTuq9rZqhUAiwmDYGmvcHraqqanbftJ+IdAauACYDB+D0N34IuNadkSKWfyWyDq0mQe4cYd6MqyvwD/b25u6EM+LjA+w7eV7Wa6mFyJhMEflQ37JlS9PAnxbTxjiiW0a9n/l26jghbseZkms6UAuUAL/BGTG6PNYBqvq7RFYgXkvQxW0pVES6ALcB3wT6AO8Cv1bVv4vIYOB93LEBXDeo6tWeY28HzsaZWO0PqnpTW+rRFhbYJleEw2HOOOMMampqmDRp0j6d/Y3JdS219nvfO2A/HNrhLOArqrrZvf+2iLyK09oTMwmKJiLjgO/jJE4f4Ywc/ZzfCrSaBKnqfX4LilHueuBEnL5EE4H/EZGvefbppaqNMY6djXO5WwkwAHhORN5W1SVtrEsg1ifC5IrKykpqamoYOnQo1157rSX9xkRprT9o5L1jPxza5VOgEOcq84huOPOTxiUiPwJ+j3OF+qvAIGCBiPxWVe/2U0bgwRL9cGcUnu1Z9TcReR84Cng9zuEXAD9Q1U3AJhG5G/gBzizGSWf9fEw2aa1l0xvrlgAZ4/9MgL132s4dcDnifmCJiPwJ+JC9oz7/2Wdx/wmUubPPR8p/GGdMQ19JEKqa9AXoD+wESoHBOP2MPnIf9H8D+7v79Xa39fccezbwZmvlH3XUUWqyG7BSOyBW03FpT3yXl5croOXl5TG3V1dXa3l5uVZXV7f5f5j2y+X4LiwsTMRTGFis2I/3fjFtA9Sr+3rjdIeJt7ynPmIH2AgURK3rAmz0c7yqJqclyEtECnA6Tt+nqmERKcaZpfgNYD/gVnf7t9g7QesWTxFbcCZUiy53GjANYNCgQQmvt/UNMtkgXsumnf41ucpG/k8NVT0kgcW9CNwkIper6nYRKQKuA14KUqGkLThXjz0ELCYqW/PsMwCn9ac7e1uC+nm2h0hBS1Brvwjs13PHI4d/KSezpdNiOT3kcnynU0uQSQ48LUGJXICBwAtAA/Bv9+9ynBGkfZXhe9oMETnAd2bl7C84U230x5lbrKV5QSKX4HdSpx/QJ8AIz/YRwFtB/ncixJoaIyJ6Gg1j0llrU2PYdC8mV0S/Dyz2U0NEqj23m02VEXTaDFX9RFXHAocA3wYOUdUTVfXjOIc2CXI67EMR+RfwfGSJ849uBw4DvqmqOyIrReRYnJ7gNTgtP7e4ZUVOgf0ZuEJEVuIkUD8GfhigngnR2hhA1mRqMkHklO7mzZubhvy3U14mF9kl7Wnlx57bbZo2Q0QKcQZZPBz4J3Cdqn7YlrKCTKA6FLgB5/K163GSorUicpeITI6qYAnwE2Ak8KmIbHOXqTiTri7BmYx1DfAlzmiREbNwxhWqxWnWqtQOujzeL/sFYTJBpMVSRJg0aRJbtmyxyR9NTmrpknabQLjjqeqLntvLoxecfj6nxCnmVpyWnzDOxVM3tqdCbT0XVwrMB3YBu5Nxvs/vkug+E5FzxYsWLbJzxmmCHO4z0db49vZ5CIVCCmgoFGpTWSa5cjm+E9EnKF7/npa229VgHYMAfYJwru5qNafA6TYz0L19MPC+3/L3KStAxTrhzBV2OfB3oA7nUrb/Bi5oawUSsSQiCfK+SSJvjKFDh9obJE3k8pdE0PiO9YE/adIkBbSkpMSS+jSUy/GdiCQoOpmJ9R7wu84kXhuSoD1x9qmLuv+F3/L3KStAxeqAz9zWn+8DJW39p4legnxJRAd95H7kl7I3Ebrjjjt06NChumjRIt/lm+TI5S+JoEmQN4n3xrkl9ekrl+M70S1BLcW6tfqkThJagrYDJwEnu0td1P2T/f6/IB2jVwOjgGPdCtSLyDZV3RigjJSLHhsicr+srIyhQ4cSCoWa+vxceOGF1NTUUFVVxcSJE1Ncc2PiC4fDbN68mZKSEmpqaqisrGzqv/bEE0/YDNgmK5WWllJRUdE0GXCs/j92QUv6iBo1OlpnH0VEGmQiNkbdV5z+x3H5ToJUdYyIdANGA+OAGcD97hVjz6nq//NbVipFvxEifzdv3rxPwmNvGpNpKisrWbhwIUVFRYwePZrNmzcTDocpLS1t9YpHYzJd9A/auXPnAhAKhRARrrnmGov/9BHvhWj1EnlVHZyoigQaMVqdS92fEZE1wNvAacC5wHAgI5Ig7xdBOBxm5syZiAgXXnghvXr1apbw2JeGyTQVFRU88sgj1NXVsWbNGl566SV69eplcWyyXuSzO9ISVFVVBdA0PETPnj3tfZAmNLGjRgMgIk+o6hlBj/OdBInIuTgtQCcBw3Bmf30B+AXOuEEZZ+bMmU1vEHDeJMZkIu80LwsWLGDGjBlcdtllvPbaa01fDjYVjMkF5eXl9OzZsynu169fz9q1awmFQimumUmyMW05KEhL0FyccXvm4gxu+E5b/mE6cQa1hpKSElTV5lAyGSuS0K9fv57u3btzxBFHcOKJJ/KTn/ykaR+bJ8xkg8WLFzNjxgzmzp3b1HWhtcEQDz74YJYuXWp9O01MQfoEBZo2I51F3kRTp05l9erVzJ07lyFDhtCrVy9CoRAXXnih/Vo2GSWS0K9du5ba2tqm9Y888kjTbevjZrLBRRddRG1tLRdddBHr1q0DWh4MESzus5WIXBm1qqt3nape5aecICNGe//5ABEZ5F3aUk5HiR4VdMaMGdTU1HDDDTdQU1PDj3/sjOI9b948qqqqmD9/PjNnzrSRRE3GKC8vp6SkhAMPPJC+ffsCUFdX12wfG+ncZINhw4Y1+xsOh1m/fj0lJSXMnTvX4jt3SJzFnwDX7vcE7gN2ALujF7/lJGOJN45KZHyIUCjUNP7PAQccoJ06dVKcS+maxlSJNW6QST1yeByVePG9aNEi7dGjR1Msl5SUKKCTJk3y/fya1Mrl+PYzTlD0OEDesd4in9UtfV7b+ECpR5JmkfcutHHAxCB9gm7EmdH9LGAhUA4ciHNV2GUByukwkY6gkQ5xmzdvZv78+axYsQKAPXv2kJeXR//+/ZuNqTJv3jzC4XCzznXGpBNvbE+ePJm6ujo6d+7M7t276dSpEwcccABbt25tasm0DtEmk0X6s61YsaLp0vfly5dzxhlnkJeXB0B+fj5Llixh8eLFzfr+2OmwnOG/9ccrQJb1IXCCe7sO+Kp7ezKwNNlZXmtLS7+UvS1AoVBIy8rK9IADDlBACwsLFdADDzxQJ02apKFQKO7Q6TbEeuqQw7+UY8V3JLYjrT75+fnapUuXpl/EeH4Zx/olbLGcXnI5vv22BEVGgY78zcvL2yfeI9tNeqFjWoKebMtxQVqCeuHM7A6wBdgP+BfwMnBP0OSrI3gvkVy6dCkARUVFAGzfvh2Ajz/+mIULF1JeXh73V7JdXWPSRUVFBevXr+fFF50JmRsbGzn44IN5//33AejWrRtjxoyJOWIuWCyb9Odt7ayqqmLu3LlUVVURCoWaOkdHFBUVMXToUDZt2tTUUmRyi6p+u60H+s2yVgHj3NtPAzfjND9dCqxPdpbX2tJan4nq6uqm1p/I0qNHDz3ttNNURPSggw7SkpISX/OD2a/n1CGHfynHiu/q6motKipqFteRPm75+flx+0BYLKeXXI7vllqCols7u3XrpqNHj9bRo0drt27dmsW+ze+Y3uiAlqC2Lv53hEuA6e7tk3HmD2vA6Rh9cSofRPSXRHV1tZaVlWnfvn21oKCg2ZtFRJq9saDlDnUmfeTyl0SsJMjbGTR6GTlypCU4GSaX47ulJCiSqI8cOTJmnHs/28vKyizm01g6J0FBxgm62XP7WREpxZlQtUZV3/RbTkeorKxsOv0VTVXp0aMHBx54ILW1tRxwwAF84xvfsE5zJqN88sknLW4bMmSIneIyGS8yKerXvva1fbaJCOPGjWv6nF+7dm3TbYt9E0SbxgkCUNUPVHWhqr4pIgcnslJtFRkP6JhjjqGgoKDF/erq6igqKqK8vJxnnnmGRx55xK6aMRlj8eLFvPzyy/us79OnD2VlZYiIjW9lMlbkc3zx4sWMHTuWxsbGffYREWbMmMGkSZMIhULcdtttlJeXN5sixsZ5M760pxkJGADcCuxIZXNW5HRB5Bxy3759WzxVEFlKSkqs6TSDkMOnC7ynw+64444WTw2UlZU1nea107uZJZfjO/p0WORzvKWrv/BxCszGBkovpPHpsLgtQSLSS0QeEJENIvKxiEwXxyzgPeBYnDGDUq6iooJJkyaxYcOGuPvW1tZSWVnZAbUyJjHC4TA//elPY27bs2cPS5cupba2NubUAcZkilAoRJcuXdi9e3er+61du5b58+fH/ByvqKho1jJkTEv8nA77PTAWZ7ToL3CuCnsCOBGYoKqjVHVBIislIn1E5FERqReRWhGZ4ue45cuXN5sVvjX2RWFSoa2xDc4kqS352te+1nRq4IknnrDTuyYl2hPfEdOmTePLL79sdZ+RI0fucwrMy6aIMX756Rh9GvBDVV0mIrfhjA30rqrOSGK9bgV2Af2BkcAiEVmlqm+1dlBLv5IjCgoKmDlzJh988IGNnmtSpU2xDbSY4Ee+CCyeTRpoc3xHfPTRRy1u69evH5999hlHHnkkEydOtFnhTbv5SYIOAN4GUNX3RGQncHeyKiQiRUAIOFxVtwEvisgTwPeAX7W13FAoxDXXXGNfFCZlkhHbX//61+1qGJMWkvXZ7aWqdprLJJSfJKgTznhAEbtxxghKlmFAo6qu9axbhXP6rU1GjhzJI4880u6KGdNOCY3t/v378+CDDyakYsYkQMI/uyNEBFVl5MiRlvSbhPKTBAnwFxGJnKTtCtwtIs0SIVU9I0F1KsaZm8xrC9C9WaVEpgHTAAYNGtRiYYcddhgLFiS0y5IxbeUrtsFffH/66acJrp4x7dKm+O7cuXOrhXbq1Iknn3ySqqoqawEyCecnCbov6v5fklERj21Aj6h1PYCt3hWqehdwF8CoUaM0VkGzZ89m1qxZyaijMW3hK7YhfnxXV1cno37GtEeb4tud/iWmwsJC/vrXv1r/H5M0cZMgVf1hR1TEYy2QLyJDVbXGXTcCiNuxzhmOwJi01ebYBotvk/baFd8RFuemI7V5xOhkUdV6YCFwlYgUicjxwJnA/amtmTHtY7FtspnFt8lEko5Zt4j0AeYDZcBG4Feq2mIPUBHZANS6d/cHPk96JTNTJj83JaraN9WVaK+gse0ek4vxnWuPM5fjew+wI2p1PrDvfBnZJZceYzdVTbtGF0jTJKg9RGSlqo5KdT3SkT03mS9XXkN7nLktF54Xe4zpIS0zM2OMMcaYZLMkyBhjjDE5KRuToLtSXYE0Zs9N5suV19AeZ27LhefFHmMayLo+QcYYY4wxfmRjS5AxxhhjTFyWBBljjDEmJ2VNEiQifUTkURGpF5FaEZmS6jp1FBF5XkR2isg2d3nHs22K+3zUi8hj7jgekW05+5xlmmx5rbI1VkXkYhFZKSJfisi9UdtOEZGwiGwXkedEpMSzrYuIzBeROhH5VEQu9XtsNkr319mPbIzxbI7vrEmCgFuBXUB/YCpwu4gMT22VOtTFqlrsLocCuI//TuB7OM/LduA2zzG5/pxlkmx6rbIxVj8GrsEZKLCJiOyPM4ryb4E+wErgYc8us4GhQAlwEvCfInKqz2OzUbq/zn5lW4xnb3yrasYvQBFO8AzzrLsfuD7VdetAE+EVAAAgAElEQVSgx/888KMY638PPOi5/xX3eeqe689ZJi3Z9Fple6zifFHc67k/DXgp6rXcAZS69z8Gxnu2Xw085OfYbFsy6XWO8ziyNsazMb6zpSVoGNCoqms961YB6ZJFd4TrRORzEfmHiIxz1w3HeR4AUNV3cd9o2HOWSbLttcqlWI1+XPXAu8BwEekNDPRup/njavHYJNc5VTL5dY6WKzGe8fEddxb5DFEM1EWt24KTYeeCy4G3cd5Q5wFPishInOdlS9S+kedlN7n9nGWSbIrvXIvVYmBD1LpI3Ys996O3xTs2G2VLnOdSjGd8fGdLErQN6BG1rgewNQV16XCq+qrn7n0iMhmYSOvPy55Wtpn0kjXxnYOx2trj2ua5vzNqW7xjs1FWPN4ci/GMj+9sOR22FsgXkaGedSOAt1JUn1RTQHAe/4jIShEZAnTBeb7sOcsc2fxaZXusRj+uIpy+IG+p6ibgE+92mj+uFo9Ncp1TJZNf59Zkc4xnfnynuqNVAjtsPQQswOlcdTxOs9rwVNerAx53L+BbQFeclr2pQD3OOebhOE2sJ7jPy19wO6Xl8nOWiUs2vFbZHKvu4+kKXIfToTXyGPu6dQ25624AXvEcdz2wHOgNlOJ8aZzqbmv12Gxc0v119lH/rIzxbI7vlD+5CXyR+gCPuQH3ATAl1XXqoMfdF/hfnCbEzcArQJln+xT3+agHHgf65PpzlolLNrxW2RyrOJcCa9Qy2932TSCMc+XL88Bgz3FdcC47rgP+DVwaVW6Lx2bjku6vs4/6Z2WMZ3N829xhxhhjjMlJ2dInyBhjjDEmEEuCjDHGGJOTLAkyxhhjTE6yJMgYY4wxOcmSIGOMMcbkJEuCjDHGGJOTLAkyxhhjTE6yJMgYY4wxOcmSIGOMMcbkJEuCjDHGGJOTLAkyxhhjTE6yJMgYY4wxOcmSIGOMMcbkJEuCjDHGGJOTLAkyxhhjTE6yJMgYY4wxOcmSIGOMMcbkJEuCjDHGGJOTLAkyxhhjTE7KT0ahIrItalU34DZV/YWIDAbeB+o9229Q1avdY7sAtwNnA9uBP6jqTa39v/33318HDx6cmMqbtPT6669/rqp9U12PVLD4zn7pEt8d/dkN0KlTJ+3WrVsCam/S1fbt21VV07LRJSlJkKoWR26LSDHwKfDXqN16qWpjjMNnA0OBEmAA8JyIvK2qS1r6f4MHD2blypXtrrdJXyJSm+o6pIrFd/ZLl/ju6M9ugG7dulFfX9/aLibDiciOVNehJR2RmYWAz4AVPve/ALhaVTepajVwN/CD9lbitdde4xvf+AZjx45l8uTJNDQ0tLdIk8NE5DwRqRaRehF5V0ROcNefIiJhEdkuIs+JSInnmC4iMl9E6kTkUxG5NFH1sfg2SWCf3SbrdUQSdAHwZ1XVqPW1IvKhiPy3iOwPICK9gYHAKs9+q4Dh7a3EwQcfzLPPPssLL7zA4MGDefzxx9tbZMYKh8OEQiHOPvtswuFwqquTcUSkDLgB+CHQHRgLvOfG8ULgt0AfYCXwsOfQ2ez9pXwS8J8icmoi6mTxbZLAPrtN1kvK6bAI91fwicCFntWfA0cDbwD7AbcCDwDfAiJNsVs8+2/B+aKJLnsaMA1g0KBBcesycODAptudO3emU6e0PD3ZISorK1m4cCEAPXv2ZN68eSmuUcb5HXCVqr7i3v8ImmLyLVX9q3t/NvC5iJSqahjnS+UHqroJ2CQikV/KrZ4u8MPi2yRSMj+73fKbPr87d+7cal0stk0yJTuavge8qKrvR1ao6jZVXamqjar6b+BiYLyIdAcinfJ6eMroAWyNLlhV71LVUao6qm9f//0Ja2trefrpp/n2t78d+MF88cUXfOc736GoqIiSkhIefPDBFvetrq7m5JNPpmfPnnz1q1/l0UcfbbZ93LhxdO3aleLiYoqLizn00EMD16etKioqmDRpEqFQiIqKijaXEw6HufDCC3OqNUlE8oBRQF8R+Zf7i/i/RKQbzq/epl/CqloPvAsMT+YvZa90ie+Impoaunbtyvnnnx+4PialkvbZ7ZbV9Pmdn+/vt3i6xPa6deuYOHEivXv3ZsCAAVx88cU0NsbqImUygqombQHWAuVx9ukPKNDTvf8xUObZfhXwUGtlHHXUURpPY2OjbtmyRU844QQNh8Nx94/lvPPO03PPPVe3bt2qK1as0B49euiaNWv22a+hoUGHDh2qc+bM0cbGRn3mmWe0sLBQ33nnnaZ9TjzxRL377rvbVI90UV5eroCWl5cn/X8BKzWJsep3AQ5w43UlTlKzP/AP4FpgHnB91P7/wGntOdg9rqtnWxmwroX/M839HysHDRoU9/lJt/iOKCsr0zFjxujUqVPbVKdckS7xHVk66rNbVSksLGz1uUm32J4wYYJecMEFumPHDv3kk0/08MMP1z/+8Y9tqleuAOo1DeI61pLMN9FonEspu0etPxY4FKcVaj+cPhPPebZfDywHegOlwCfAqa39r1hJ0D333KPf/OY3tby8XHv16qV/+MMfdMKECbps2TL/r5zHtm3btKCgoNmb4fzzz9fLL798n33ffPNNLSoq0j179jStKysr0yuuuKLpfjYkQdXV1VpeXq7V1dVJ/1/p8iXhxqUCF3jWhYD/A/6Iczmxd/833e2R4/pFHfdmvP+ZifGtqrpgwQI955xzdNasWZYExZEu8a0d/Nmtum8SlO6xXVpaqosWLWq6/8tf/lKnTZvWprrlinROgpJ5OuwCYKGqRjeHDsHpA7EVWAN8CUz2bJ+Fcwqh1n1DVWqcSyxjWbVqFa+88gpnnnkmGzdupF+/frz66qtcffXVjBs3jocf3ttf9fTTT6dXr14xl9NPPx2AtWvXkp+fz7Bhw5qOGzFiBG+99Zav+qgqa9asabbu17/+Nfvvvz/HH388zz//fNCHmHKlpaXMmzeP0tLSVFelw6jTn+dDnISmabX79y1gRGSliBQBX8HpJ7QJ50thhOe4Ee4xgaV7fNfV1XHllVdy001xh4kx6cc+uz2iY3vGjBk89NBDbN++nY8++oi///3vnHpqQq5vMKmQ6iwsEUusX8pjx47V2bNnx8hJ2+aFF17Q/v37N1t311136YknnrjPvrt27dJDDjlEb7jhBt21a5c+9dRTWlBQoOPHj2/a55VXXtG6ujrduXOn3nvvvVpcXKz/+te/ElbfbEN6/VK+CvhfoB/Or94VwNVAX5zOoCGgK84VZK94jmvTL+VMjO/p06fr9ddfr6pqLUE+pFN8d/QS3RKU7rH99ttv65FHHql5eXkK6AUXXNCs5cjsixxtCUqp1atXc8455ySsvOLiYurq6pqtq6uro3v3fS9+KCgo4LHHHmPRokUMGDCAOXPmcO6553LQQQc17XPsscfSvXt3unTpwgUXXMDxxx/P4sWLE1Zfk1RX4yRBa4FqnFNh16rqBpwE6FpgE87pg/M8xyXklzKkd3y/8cYbLFu2jEsuuSRh9TO5I51je8+ePZx66qlMmjSJ+vp6Pv/8czZt2sTll1+esPqajpWVSVBtbS0NDQ2+T9NMmDCh6Sqt6GXChAkADBs2jMbGRmpqapqOW7VqFcOHx76454gjjmD58uVs3LiRp556ivfee49jjjmmxTqISKS1oMPl4lVe7aGqDap6kar2UtUBqjpdVXe625apaqmqdlPVcaq6znPcl6parqo9VLW/+phSIJZ0j+/nn3+edevWMWjQIAYMGMCNN95IVVUVRx55ZFserskh6R7bX3zxBR988AEXX3wxXbp0Yb/99uOHP/yh/YDNZKluikrEEn264PHHH9fjjjtOE+273/2unnfeebpt2zZ98cUXW7zCQFV11apVumPHDq2vr9fKykodPHiw7ty5U1VVN23apEuWLNEdO3ZoQ0OD/uUvf2nx6pqO0JFXebUVOXy6INPiu76+Xj/55JOm5bLLLtNQKKSfffZZwuucLXI5vr2nw9I9tlVVDznkEL3uuuu0oaFBN23apGeddZZOnjw54XXOJtjpsI61atUqRo4cmfByb7vtNnbs2EG/fv2YPHkyt99+e7NfExMmTOD3v/89APfffz8DBw6kX79+PPPMMyxdupQuXboA0NDQwBVXXEHfvn3Zf//9+dOf/sRjjz3WrONeR6qoqKC8vLxdYwaZjpPu8V1YWMiAAQOaluLiYrp27UqQ8bxMbkr32AZYuHAhS5YsoW/fvnz1q1+loKCAm2++OeF1Nh1DnCQts40aNUptgsnsJiKvq+qoVNcjFSy+Uy8cDlNZWUlFRUVSrobM5fguKipSm0A1u4nIdlUtSnU9YknqtBnGGJMNKisrmT9/PoBNM2NMFsnK02HGGNMe0RcLJOqUsV2EYEx6sSQoy9iHrDHtF2n5qaysBBI3MGh0ucaY1LLTYVnGmu2Nab9Ii0+iLxZIVrnG5Bp3Muv5wDRV/bKt5VgSlGXsQ9aY9ou0/GRKucbkGlXdLSLjgT3tKcdOh2WZXJzPyxhjTE66GfidiBS0tQBLgowxWc/6yhmTlX4BVABbRWS9iHwQWfwWYKfDjDFZz/rKGZOVzm9vAZYEGWOynvWVMyb7qOry9pZhSZAxJutZh2Rjso+IdAGuBCYD+6lqT7ez9DBV/S8/ZVifIGOM8bD+Q8Ykn4gcJyKXuElL9LZf+SzmZuBwYCoQmQPsLeBnfuthSZAxxnjYgIbGJJeIfA9YDIwD7hWRRSJS7NnlNz6L+g4wRVVfxr1UXlU/Ag70Wxc7HWaMMR7Wf8iYpPs1cKqqviYi3YA7gOdEpExVNwPis5xdROUxItIX2Oi3ItYSZIwxHjbWljFJd6CqvgagqjtU9QLgeeAFEenH3lNb8fwVuE9EDgEQkYHAfwEP+a2IJUH/v70zD5OrKPf/5zuGbANJCAQ00QyLidEgiSSiBpRoEgWUqBmuSyJeGQUUkIvguBFEMFyXuXiRKyhcEwVZrwwI/IgooEEQEREBhTSMIMMSlqBJSAYMJnl/f1T15EzT3dMz03u/n+epp/ucWs7b1e85562qt6ocx3Ecx8lKiXzknpE0JXnCzNqBq4HbgEIXP/wK8Dfgz8A4oAtYA5xRqCBuBDmO4ziOk5US+chdAyzOPGlmpwE/AkYUUoiZvWRmnzOzHYHdgZ3M7HMU3pNUOiNI0ipJ/5S0KYYHE3GLJXVL6pH0M0njE3HjJV0d47olvayiHMdxnNLgz24nSXt7O21tbXl95FKpFIcffjitra0F9RiZWbuZnZ4j7htmVpBtIunsRL61ZmaSRgLXFpIfSt8TdLyZ7RjD6wAkTQfOB44gWG4vAOcl8pxLcHbanTDt7fsxj+NUFZKmxJfFxYlz/pJw6gF/djtAYT5yHR0ddHZ2ctVVV5V7VuW+knqNKUmjgeuBZwotoBKzw5YA15nZbwAknQqslrQTYYpbK7CPmW0CbpN0LeGmK3TdAMcpF+cCf0gfJF4S7wXuBi4gvCQ+kkiffknMBK6XdK+Z3V9OoR1nkPiz28lKe3s7GzZswMzKPavy/cBNkjYQnr0/Bx4Eji60gFL3BH1D0nOSfitpbjw3Hbg3ncDMHia8GKbGsMXMHkqUcW/M4zhVg6SPAOuBmxOne18S8UVwKrBI0k6SmgkviVPNbJOZ3Ubosj2i3LI7TgH4s9spmGnTpnHllVfS2dlZ1lmVZrYROAQ4ErgHuNfMjjKzyvsEAV8E9iIsWnQBcJ2kvYEdgQ0ZaTcAO8W453PE9UHS0ZLuknTX2rVriy17w+Or5uZG0hjC7IOTMqL8JeHUAyV9dkPf5/eWLVuKKbszCGrpeS/pjGQATgTuJOjausT5gijZcJiZ/T5xeKGkjwKHApuAMRnJxwAbCV2queIyy7+AcIMye/bsgq0+pzB81+28fB1YbmZPSH3W9Mr3ktjKAF8SxC7dyZMnF0FkxymMUj+74zV6n9/Nzc3+/K4wlXreS7oauBC43sz+VWC21+Q4//M8cTkpp0+QEVaBvB+YkT4paS/CdLiHCDfSMElTzKwrJpkR8zhlxFfNzY6kmcB84E1ZokvyknAj36kw/uyucyr4vL+VsAHqckn/B/zEzG7Pl8HMjiymACUZDpM0TtJ7JI2UNEzSEuAdwA3AJcBhkt4e/STOAK4ys41m1gNcBZwhqVnSAQTHp5+UQs5SU0tdjJn4qrk5mQvsATwm6Wng80CrpLvJ/5J4iPiSSJTlLwmnqvBnd2NSjOf9YN53ZvYdM9uPoGPrgcskdUn6ahyC7RdJYyXtL+ldyTAQIYoegAmEWTMb4w+7A1iQiF8MPAb0EBZNGp+IGw/8LMY9RtgcLe/1Zs2aZdVIW1ubAdbW1lZpUWoe4C4rga4ONACjgVcmwn8BV0adn04Y8no70AxcDFyeyHs5cFmMO4AwHDa9v2tWq347xaOK9Lusz24zY/To0aWoUqfM5HvfAT1WmP69neAruTU+H28CZuRJ/4mob88QVo5Oh0cKuZ6ZlWY4zMzWAm/OE38pcGmOuH8AHyiFXOXGh5TqDzN7gbA+CgCSNgH/jDq/VtKnCS3mXQg3cLLr9lhgBfAsYYO/z5hPj3eqCH92O4NlsO87Sa8DPkYwsF8i9B6+D1hLeGb+DNgzR/YzgcPN7OeDEhpQtKZqmtmzZ9tdd91VaTGcEiLpj2Y2u9JyVALX7/qnkfW7ubnZenp6Ki2GU0IkvWBmzVnO30VwL7gCuMj6OuWn0/zNzLIaQZKeASaa2dbByuZ7hzmO4ziOU1YUptZeTjBijstmAAHkMoAi3wKWShq0LeNGkOM4JaE/R8lanjjgOM7QiL5CpwNDWSjqc8BSYKOkx5Kh0AIqsW2G4zgNQH9rj/haVI7T8PyJsJDsYFtCHxuqAG4EVYhUKkVHRwft7e0+Dd2pS/pzlKyHiQN+HzvOkFgF3CDpx8DjhDWpADCzFf1lNrNbhiqAG0EVwlvBTr2TXntksPG1gN/HjjMkDiBMaT8o47wRZtK+DEmnmNmZ8XvO7THM7KuFCOBGUIWoh1aw4zQ6fh87zuAxs3cOIturE98HvE1GJj5F3qkJGnkKset3/dPI+u1T5OufXFPkM9KIsD0LAGa2reSC4bPDSorPfnEcx3Gc7EiaJOlqSX8nzBL7VyLkyrNXIaFQGXw4rIS4v4DjOI7j5OQHhBX45wG3EPYQ+xqwMk+ev7J9U99cGPCKQgRwI6iEuL+A4ziO4+RkDjDZzHokmZndK+mTwO3A/2bLYGZFHcHy4bASkrkzrw+POY2O3wOO4yTYyvbFEtdLmkDYEHVSIZklvV9SQT0+uXAjqIykh8c6OjoqLYrjVAS/BxynL4U0DOq48fB74ND4/ReEPcSuAgqdCXIG8LSk70l6y2AE8OGwMuLDY06j4/eA4/SlEN/RZJr29vZ6WqDzCLZ3xpwInAzsBJxdSGYzmyFpBmHl6E5JPYRd6C82s0cLksDMaj7MmjXLnPoGuMuqQNcqEWpVv1evXm1tbW22evXqSotS9TSyfo8ePboYVVizFHKfJNO0tbUZYG1tbWWUcmgAPZbxvxMcly8ERmTGDSYQHKXnA/cShtl+AywBmvLl854gx3FKgs+OdJz+KWTl9GSaeulNNbOtkt4NDHk9IEl7E3qDPhbL+yrwGHA80AosypXXjSDHcUpCvTysHaeaqIftZhL8N3C6pNPMLOfaQLmQdBxhSG0KwZ/oCDO7IxHfCTybrww3ghzHKQl19rB2HKf4fBZ4JXCSpLX03UB1cgH5DwHOAq41s82ZkWb2gqScvUDgRpDjOI7jOJXhY0PJbGbvKyDNL/PFuxFURFKpVD157TuO4zhOyTCzWwaaR9JPSPQY5Sn744WUV5J1giSNkLRcUrekjZLukXRIjNtDkknalAinZuRdIel5SU9LOqkUMpYCXwPFcRzHSVPH6/sUBUnDJZ0hqUtST/z8uqSRebL9FXi4gFAYxZialmWqWjNh/489CIbW+4CN8XgPghU3LEfebwC3AjsDrweeBg7Od71yTCEe6DRGp7hQJVOIgRHAcqA76vQ9wCGJ+HlAirAfzq+Bloy8K4Dno16fVMg1a22KvN8HA6cW9Dvx7N6UCKdm5B2wftfzFPlanM5eCsgyRT6cZjlwG8G35w3x81ZgRbb0pQjlvLnuI0xV688IWgO8O3H8deDyfGWX4yWRTZn9YV8+quglkc/A3xXYAPwbMBLoAO5I5B2wgW9l0u9i4g/+gVMj+l30BqxZfRtBme+IRn1n5DGC/g6Myzg3HvhHtvQ5ypgbje9fxM93FprXrExGELA78E9gWuJGehJ4AvgRsGtMt3OM2z2R93Dgz/nKr1RPUL6HfaMqe6molpdEtpAw8I8Gbk+cbwZeBKbF4wEb+FYm/S4mrvsDp0b0u+gNWLP6NoIyadQGQh4j6H5gYsa5ScD92dJnyf8pwhT4/wSOAc6MxvdRheQ3K8NiiZJ2AC4BLjSzlKQdgTcTull3Ac6N8e8BdozZNiSK2EBYRjuz3KMJLx0mTy5kJt3QyDbdt729nQ0bNrB+/XpSqVQfZ2hfKK4xkLQ7MJVwM3+GsFopABZ2Rn4YmC7pGeBVyfj4/QNlFLco9DcBwKfG1w8Z+p2mW5IBNwLtZvacpJ2pE/0uJb52Fkh6V+LwJ8ANkv6H0CnyGuA44KICi/sCsMDMevVO0hVAJzl2oc+kpEaQpCbCj3yJsHIjZraJ7ZujPSPpeOApSTsRxpgBxhB6jtLfN2aWbWYXABcAzJ49u19P8VIwbdo0xo4dy4oVKxg3blyfB78re/2Tw8Bfm5EsbcQXbODHsstq5A8EN/Abg1I1YGPZvfo9fPjwkshfjXgDAQh+QJl8JeP4GOBbBZS1C/BAxrkHCUNqBVEyI0iSCD92d+BQy70aZNqAaTKzdZKeAmYQWhnE7/dnzVlhUqkUGzZsYNGiRS8zdlzZ65tsBj7BiB+TkTRtxBds4EN1GPm5enzcwK9/StmAjWX16ndzc3NF9NupDGa2ZxGLuw34jqQvWlgYsZngm3Z7oQWUZIp85PsE57jDzOzF9ElJb5H0OklNknYBzgFWmVm6BXERsFTSzpKmAUcBPy6hnIOmo6ODzs5Oxo0b5+sCNRAZBn5rwsC/n2C0p9M1A3sTxrfXAU8l46liAx9yL/mQNvBd5+uTPPqdSZ8GLDWm305d8GmCnm2ILgfr4/ExhRZQqnWCWqIQM4GnE+sBLQH2Am4gtBD+AmwGPprIfhphjn83cAvQYWY3lELOodLe3k5bW5u3iBuPrAY+cDWwj6TWuM7FV4H7zCy9SEjNGPiQX799/ZO6pu4bsAPF9b14SFqd+P64pMeyhULKMrOnzOwdwJ7AYcCeZnaQma0pWKBCPairOZRr9ozPeqkcVMnsGaCF0AL+J33XS1kS4+cT1gl6EVgF7JHIm1xH5RlqeJ0gXzKiuNSCfhMaq38Degi9PhcBr0zkHZR+18LssEad1VUsSMwOAw5MfD8oV7D8ejqaMCPsWsKSDiPypc8XfNuMAeAOoY6ZdQPKE38TYSmIbHGbgbYYappMv6BUKsXChQvp6uoC/P6oVfrTb+CyPHnrRr8zSet5a2srn/zkJ31rpCFgZrclvr9s2wxJryCMCOXbUuNcYDbwc8IyOrsQNmMdMKX0CapZcnV9tra2MmXKFFpbWweUz3FqlVw6nekX1NHRQVdXF1OmTPHhYafuSOv78uXLWbFiBaecckqlRapnhgH9VfDBhDWpvkBYZbrfjVRz4UZQFnI5hHZ2dtLV1UVnZ2fWfEuXLmXFihUsXbq0HGI6TslJ3wsLFy7Ma9yn/YeuvfZabyE7dUvwGd/+6ZSM/iq42cyeAjCzx4Gxg72QD4dlId2S3X///Zk6dSonn3wyd955Z28PUK6Wbhyr7P10nFolPT2+tbWVW2+9la6uLjo6OnqHuTKnz/uSEE49k9b3trY2xo4d672dpae/l+gwSe9ku7GUeYyZ/aqwK1WBM95QQ6kcR6dMmWKAxXUsrLW1NW96dwwtHVSJ42glQiUc/5OOoAPdMsYZOI2s37XqGO3P+8IhY9sM4F15wsHAVsujM8CjBCf9XOGRfPmTwXuC8nD22Wdz4oknMmHCBG6//fZ05efEW8NOLZN0/E86PufaMib56Ti1Sn/bwEB2ffeJMkOivwrLO0XezPYomiSFWkvVHIrVUs5l2bvFX3lo4JZyMXuC8uny6tWrbdGiRdba2mrXX3/9gHXe75PB08j6XcmeoNWrV/f2+A+0V9P1vXDIsYFqMQNw7aDylVqwcoRivSS8i796aeSXRDGNoP669dPxhb4Ycg2hOQOjkfW7kkZQUt8HY/g7hVEmI+gfg8nX8MNhya7QwXbxF9Kd6jjVQH/d+sn1UDo7O/u9F3INoTlONZJ8VgOsX7+e1tZWli1b9jJd9md6g1Bq66wcYSgt5XRLoKWlxVpbWwfVCvAWcOmhgVvKxXaMvv7663tbvmbZu/UL7er3IYHi0Mj6Xc6eoOSzOvO57b2apYMS9AQRtiVKhheSx4WW09DrBKV3gZ84cSLd3d10dnYye/ZsVq5c+bJ0+RZB9D3EnFrixBNPpKuriw996EMcfvjhPPLIIy9Lk2utrEx8M1WnGsn1zE4/q1tbW0mlUowaNYpUKkUqleqjy/5MrwnUTyiMYltnlQiDbSmnrf3ddtvNCOsS9I4Pe6uguqCBW8rFdoxesGCBjRo1qlffW1pafPpvhWlk/S5FT1D6mb1o0aKsepyOTwd/tpcW3CeoukiPC++///7ceuutPPfcc71xzc3NnHzyyX32QXJfB6de6Ojo4MYbb+xzburUqcybN6+PfvtyD06tsnLlSm6++WbmzJnDH//4R7q7u4G+09jb29t5/PHHuf/++5k+fbo/2+uDQS3j3ZBG0NKlS+ns7OTCCy9k69atjB0bVtwePnw4Z511FmeddZbvg+TULMnVntPOzY888gjHHlBej+IAABMwSURBVHsskyZNYtiwYWzZsgWAiRMnMmbMGHcAdWqetN7ffPPNdHd3s27dOp5//vne/R6TG59OmzaNX/7yl5UW2Skut/WfJAul7qIqRyh0uCDdxT9nzpw+XaFNTU2938eMGdNnSMzMHZ+rARp4uGCgw2GZ09xbWlr6DH9NmDDBRo0aZXPmzLFFixb1robuw1+Vo5H1u1jDYWldnjlzpo0aNcrGjx9vCxYscJeGKoAyDIcNNjRUT9App5zCVVddxYQJE/qc37ZtGwBNTU29LYfkRpA+HObUEu3t7axfv56nn36aJ554onc4AGDEiBGsXbuWtrY2li9fTiqVYty4caxfv95Xv3VqllQqxR133AHA6tWr2bx5My+++CK333474M9wJw+VtsKKEfprKadXwh0/fnyfHqDMMGfOHG8NVyk0cEt5MI7RCxYsyKrju+22m7W0tPROj0+T6QjtjtHlpZH1e7A9QWkdTS/5kOu57r2clQfvCaosJ5xwwsucQZPssMMOLFy4kGXLlrlfhFOzrFy5kiOOOIL169cjZfcRfPbZZwE49thjueGGG3r1PdMR2vdFcqqZVCrVO3nliiuuoKenpzduxIgRzJo1i+bmZnbaaSfMzHXZyUldG0GpVIoTTjiBm2++OW+6uXPncuWVV5ZJKscpPueffz6f/vSn+00nCTOju7ubhQsX9hn2TeLDB061kjSARo0a1ccAAti8eXMfoz495Ou67GSjrhdLTE8HTvv8ZKO5uZlzzjmnjFI5TnFJpVIcd9xxedNIYtSoUZgZLS0ttLS00NXVlXMxRF8E0ak2UqkUra2tHHDAAb3Ll7z44osvS9fc3Exra2vvseuyk4+67QlKpVL89Kc/zZtmzJgxXHbZZX5zODXNCSecwNatW/OmMTMOPfRQzAxJtLW1FbQ3mONUA6lUinnz5rFmzZq86dI9Q52dnRx66KFlks6pZaqyJ0jSeElXS+qR1C1p8UDyp1Ip5syZw8aNG3OmmThxIpdddpnfKE5ZGapuZ7Jy5cq8/m5pJk6cyLJlyxg3bhydnZ10dnZ669ipelauXMmkSZOYMWNGvwYQwIEHHujbXTgDolp7gs4FXgJ2B2YC10u618zuLyTz/PnzWbduXc74pqYm1qxZ460FpxIMSbeTLFmyhEsvvbSgtG9961t790QC9/VxSoOk8cBy4N3Ac8CXzawwJc3CkUce2evMn42mpia2bdvGbrvtxoEHHsiZZ57phr0zIKrOCJLUDLQC+5jZJuA2SdcCRwBfKqSMJ598Mm/8vvvuy3777ecvAqesFEO3E2UVlG7BggWMGTOGZcuWAb4dhlNyimbkv/GNb8xrAAGcd9553Hnnnb7iuTNoqs4IAqYCW8zsocS5e4GDhlrw2LFjmTdvnrcWnEpRMt1OMnz4cLZt28bSpUs57bTTilm04+SkWEb+yJEj2bx5c3/X4thjj+WYY47hmGOOGYrYToNTjUbQjsDzGec2ADslT0g6GjgaYPLkyQUVPH/+fJ8K71SSgnQbBqffaQ477DDXc6cSFMXI788A+sEPfuCGj1M0qtEI2gSMyTg3Bujj5WxmFwAXAMyePdvyFbh48WJGjhzpw19OpSlIt2Fg+g1hgbh77rmHjo4O13OnUgzKyB8+fHhBhTc1NXHddde5H6dTVKrRCHoIGCZpipl1xXMzgILHlMMq3Y5TdQxZt9Pk0nH393EqyKCM/Obm5n4f2AsWLOCcc85xNwan6FSdEWRmPZKuAs6Q9CmCc937gTmVlcxxhobrtlPnFMXI90asU05UjQoXp1muABYAfwe+lG+apaS1QHfG6V0JUzSd/NRKPbWY2YRKCzFUBqrbMU82/Yba+e/S1Jq8UD6Z60W/LydsXJo28lcCc/LNDpO0Dchc+nkYsKVUctYRtVJPo8ysOtclrEYjqBhIusvMZldajmrH66l2qbX/rtbkhdqUuZIMxsjPUY7XewF4PQ2dqhsOcxzHcWoTM/sH8IFKy+E4hVKV3VOO4ziO4zilpp6NoAsqLUCN4PVUu9Taf1dr8kJtylwPeL0XhtfTEKlbnyDHcRzHcZx81HNPkOM4juM4Tk7cCHIcx3EcpyGpOyNI0nhJV0vqkdQtaXGlZaoEklZJ+qekTTE8mIhbHOumR9LP4rTWdJzXX5VTbf+RpBGSlkdZNkq6R9IhMW4PSZbQw02STs3Iu0LS85KelnRSGeX2e6TK8LoNuG6Wj7ozgoBzgZeA3YElwPclTa+sSBXjeDPbMYbXAcS6OJ+ws/PuwAvAeYk8Xn/VT7X9R8OAxwkbZY4FlgL/J2mPRJpxCV38euL814ApQAvwTuALkg4uh9ARv0eqC6/b7bhuloG6coyW1AysA/ZJ72Qs6SfAk2b2pYoKV2YkrQIuNrMfZpz/T2APM1scj/cGVgO7ANvw+qtqakXHJd0HnA78EfgbsIOZvWxlW0lrgE+Y2S/j8deBKWb2kTLIuAq/R6qGWtHtcuC6WT7qrSdoKrAlrQCRe4FGtYS/Iek5Sb+VNDeem06oEwDM7GFCy2EqXn+1QNX/R5J2J8iZ3CqhW9ITkn4kadeYbmfgVST0kfL/Fr9Hqgev2764bpaBejOCdgSezzi3AdipArJUmi8CewGTCGtJXBdbDTsS6iRJuo68/qqfqv6PJO0AXAJcaGYpwr5bbyYMd80iyHlJTL5j/EzqYzl/i98j1YXX7XZcN8tEvW2bsQkYk3FuDLCxArJUFDP7feLwQkkfBQ4lfx1tyxPnVAdVq+OSmoCfEFqmxwOY2SbgrpjkGUnHA09J2onwWyDI/8/E97L8Fr9Hqo6q1e1y47pZPuqtJ+ghYJikKYlzM+jbLd+oGCBCXcxIn5S0FzCCUHdef9VPVf5HkgQsJzhktprZv3IkTTshNpnZOuApEvpIZX+L3yOVxes2N66bpcLM6ioAlwOXAc3AAYTuwOmVlqvMdTAOeA8wktDbtwToIYwZTyd0mb491tHFwOVef7UTqvE/An4A3AHsmHH+LcDrCA2uXYArgF8n4r8J3ALsDEwjGEUHl0Fev0eqMHjdum6Wvb4rLUAJFGg88LOoNI8BiystUwXqYALwB0I36Pr4clqQiF8c66YHuAYY7/VXO6Ha/iOCv48RhrQ2JcIS4KOE2WE90cC5CHhlIu8IYEV8sD8DnFQmmf0eqcLgdeu6We5QV1PkHcdxHMdxCqXefIIcx3Ecx3EKwo0gxykBkj4haVP/KauHuDz/JyotRzUhaQdJD0p6Rz/pfizp/5VLrnIj6XBJBQ0bSNpN0lpJry61XI4zVNwIcuqK+DKyGLZIekzS9+PCfE6DEP//w4tQ1NHAGjP7TSw3vRfa7CKUXZeY2bME36/TKy2L4/SHG0FOPXITYSXiPYBPAYfRd3+dmkXS8ErL0CjEaf8nEKb+OwPjR8CS5OaejlONuBHk1CObzexpM3vCwp5UVwDvTiaQdJKk++Juy09K+qGkcYn4pyR9JHF8m8IO6cPi8Wtjj0DeLn9Jh0l6SGFH6F/HdT3ScXtLukZh9/QeSXdLel9G/kclfU1hp/X1xNWWJU2UdImkv0t6QWHn9ncm8h0j6a+SXoqfR2WU+1pt36n6wczrxjSTJF0uaV0M12esQZLt946NPW9PxbJXS/pwIn6RpD9L2izpcUmnRGMj+Xs/n1HmKknfy0izVNL5CrvPPyGpPRkfv/40/kePxvOvifX9j1hnqeR/nIVZhM1dk8Ncf4uff4hlr8qQ9T+iPq1T2CJkdCJuhKSzJT0T6+YOSQcm4ufGMndNnOvT8xSH586RtCZRh99MpP+YpD9EXX1W0k8lTcpyjXmSfh/r4S5J+2X8jo8r7EL+gsIw3+4Z8Xnr0sz+AqwBFuWpX8epOG4EOXVNNDoOBjIX79sGnEhYd2MxsD/wP4n4W4C5sYzRhK0fNgPpYZC5wMNm9kSey48ATgOOBN4GvAK4KvHS3xH4ObCAsKhZZ4yfllHOSUAqXvsrChtN3kLo6foA8EbgjMRv/iDwPeBsYB/gu8B5kg6L8U3A1YT7/21AG2E39xGJMkYDvyZMez8opnsKuCn5Yk8Sf9fKmP5I4A1R9pdi/Czgp8BVUeYvAV8mri49QD4H/BnYD/gW8G1Jb4txb46fRxF6BNPH5wGjCbvVTyf8/+vzXOPthP84mWb/+HlwLHtRRvp9gPnAh4EPAv+RiP92PN8GvCnKf4OkV/X/c3s5IZb7EYKB9mHgwUT8cILOzQDeB+xKWDcmk28Q6n8/4O/AJWm9lPQW4MeE7RpmAteR0K9IIXV5J0EXHKd6qfQcfQ8eihkID+8thHVqXiSsX2PA5/rJdzDByGmKx58GHozf5xN2av4x8OV47mLgh3nK+0S87gGJcy3AVmB+nnx3AEsTx48C12WkOYqwhsiuOcr4LbAiS73cFr+/O8oxORF/YJT3E/G4DeiCsIxGPPcKwgvzQzmuu4BgXL4+R/wlwK8yzn0NeCLj934+I80q4HsZaS7LSNOVUW8GHJ6R5j7gtAHo0tnALRnn9ohlz85Sv48Dr0ic+1/gpvi9mWAMfjyjPh8GlsXjubHsXXNdDzgHuDn5v/TzG6bF/K/OuMZ7EmkOyEhzKXBjRjk/BGwgdQl8B7i10Pr24KESwXuCnHrkN4QWbLp3ZyXh5dGLpHdJujEOpWwk9E4MB14Zk6wCpsZW+lxCr8iq+B1CC3dVP3JsI7SGATCzbsIQwRuiDM2Svi3pgTh8sonQ2zM5o5y7Mo7fBNxnZs/luO7rCYZQktvS143xT5rZY4n430d508wC9gQ2Kswa20RYfXZnYO8c130T8JSZrR6gXJMkZe551B/3ZRyvAXbrJ893gaWSfidpWeyZyscotu9pVggPmNnWHDLtDexA4vfHtL9j+/9SCD8m6PZDks6V9N7YsweApP3iMFV31Ou07mTqVLL+1sTPtKyvj3IlyTwupC5fJNSh41QtbgQ59cgLZvZXM/uzmZ1A6LY/NR0pqQW4ntC782+EF35bjB4OYGEH9KcJ3f1z2W4EHSDp9cCr6d8Igu17ZWXjv+L1TyUYVTMJRlOm83NPAdcphIGsjNoE3BNlSoapwPlFkidJWrZthD2SkuyQJX3m8KbRz/PMzJYTDLsfEX7H7ZK+lifLcwSjr1AGLFMiHWw3QpO/v89vN7O7Cb1DX45lXwjcKKkpDpP+AngBOIIwDHhwzJqpU0lZe/dzK0DWtByF1OV4YG2hZTpOJXAjyGkETge+KGliPJ5NeCl8zsx+Z2YPAROz5LsFeG9Mv8rMHiW8GL9A//5AEO6vtA8JkibH66R7Sg4ELjKzTjO7D3iC3L0sSf4E7Jt0oM1gNWGII8mBwAOJ+EmSXpOI35++z4O7gdcCz0WDMhn+kUeuV0UjcSByPWFm6Z2u1xJ8bQCQNJIwpDNQ/kUYbuqDBWf5C8zsQ8BXCVPgc/En4HXJnhaif1O2svvh4Zi39/dLegXB1yr9v6QNhqSP0MzMgsxso5ldaWafIejnuwj/1TSCD9BXzOw30ZDvr3csG6uBt2acyzwupC73IeiR41QtbgQ5dY+ZrSK8aJbGU10E3T9R0p6SPkpw7MxkFfAh4K9mtjZx7mMU1gu0BThb0tskzSS02u8nTOGHsOvzB+MQxhsJfkYjCyj3UuBZ4BpJb5e0l6SF2j47rAM4QtJxkqZI+ixhH69vx/ibCI7WF0maGR2K/zvKm+YSwl5e10g6KNbTOySdpdwzxG4mDKt1SnpPzLNA0gdi/FnAQQqz3aZKWgKcnJAL4FeEqdVzJU0n7Cs2rIA6yeRRYJ6kVyquESXpu5IOjvU1k9BL8kCeMn5N+D/2TZx7ljDM8x5Ju0saW4gwZtYDfB/4lqRDo6H4fcKsq/TyDX8l+BWl6+fdbNdZ4m84SdJHJb1e0msJTv3PEwzoxwh+bcfH3/he4OuFyJfBOcB8SV+O+nMUwRk7KUfeuozO87OAGwZxfccpG24EOY3CWcAnJbXEXpf/IMxceoCwltDns+RZRXgBr+rnXC42A2cSFo77PeF+W2Rm6eGHkwgv1VsJs8TuiN/zEl+oBxFefNcBfyH0dlmM/xnwWcIMqgfibz3WzK6L8dsIL7WmKNdFwLIob/oaLwDvAB4hzOhKEYy4nYF1OeTaBhxC8Hu5mNCj8F22DzHeTRj+a40yfzOG7yWK+QbBELoG+CXBZ+hP/dVJFk4mDGU+nsjfRPARewC4kWDk/XuuAszs7wRfsSWJc1sIM7Q+RfCluWYAMn2RsFzDjwhDjfsCB5vZU7HsfxFmfe0F3Ev4T7+SUcZGoJ0wbHo3oafoEDN7IRrq/06YMfgAYZbYSQOQL/0b7wA+CXyG4Du0iODAnqS/unw/8JiZ9avPjlNJfANVx3GcHMTeqF8DrzWz5ystT60g6U7gbDO7tNKyOE4+vCfIcRwnB2Z2P6GXcM9Ky1IrSNoNuJLs6xM5TlXhPUGO4ziO4zQk3hPkOI7jOE5D4kaQ4ziO4zgNiRtBjuM4juM0JG4EOY7jOI7TkLgR5DiO4zhOQ+JGkOM4juM4DYkbQY7jOI7jNCT/H9fIZU+UlLV5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 16 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# File names of the raw barcode counts\n",
    "raw_data_files = [os.path.join(data_dir, dirname, filename) for dirname, filename in itertools.product([\"Rhodopsin\", \"Polylinker\"], [\"library1RawBarcodeCounts.txt\", \"library2RawBarcodeCounts.txt\"])]\n",
    "raw_data_names = [\"Library 1\\n+Rho\", \"Library 2\\n+Rho\", \"Library 1\\n+Polylinker\", \"Library 2\\n+Polylinker\"]\n",
    "comparison_columns = [\"Rep 1 vs 2\", \"Rep 1 vs 3\", \"Rep 2 vs 3\"]\n",
    "fig, ax_list = plt.subplots(nrows=4, ncols=3, figsize=(8, 8))\n",
    "\n",
    "# Read in each dataset\n",
    "for row, filename in enumerate(raw_data_files):\n",
    "    row_df = pd.read_csv(filename, sep=\"\\t\")\n",
    "    # Get all 3 pairs of combinations and plot them\n",
    "    for col, (x, y) in enumerate(itertools.combinations([\"RNA1\", \"RNA2\", \"RNA3\"], 2)):\n",
    "        rsquared = stats.pearsonr(row_df[x], row_df[y])[0] ** 2\n",
    "        ax = ax_list[row, col]\n",
    "        ax.scatter(row_df[x] / 1000, row_df[y] / 1000, color=\"k\")\n",
    "        ax.text(0.02, 0.98, fr\"$r^2$={rsquared:.2f}\", transform=ax.transAxes, ha=\"left\", va=\"top\")\n",
    "        max_value = max(ax.get_xlim()[1], ax.get_ylim()[1])\n",
    "        ax.set_xlim(right=max_value)\n",
    "        ax.set_ylim(top=max_value)\n",
    "        \n",
    "# Add \"axis\" labels\n",
    "fig.text(0.5, 0.025, \"Raw barcode counts (thousands)\", ha=\"center\", va=\"top\", fontsize=14)\n",
    "fig.text(0.025, 0.5, \"Raw barcode counts (thousands)\", rotation=90, ha=\"right\", va=\"center\", fontsize=14)\n",
    "\n",
    "# Add column labels at the top\n",
    "for col, text in enumerate(comparison_columns):\n",
    "    ax_list[0, col].set_title(text)\n",
    "    \n",
    "# Add row labels on the right\n",
    "for row, text in enumerate(raw_data_names):\n",
    "    twinax = ax_list[row, 2].twinx()\n",
    "    twinax.set_ylabel(text)\n",
    "    twinax.set_yticks([])\n",
    "\n",
    "plot_utils.save_fig(fig, os.path.join(figures_dir, \"supplementalFigure1\"), timestamp=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Rho data, assess calibration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scrambled sequences from L1 and L2 are drawn from the same distribution, KS test p = 0.087, D = 0.14\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARQAAAEUCAYAAADqcMl5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcFNW99/HPV3CZCwyIchGCOBowGlTIFROJD5FEk7gGvPK8NAKSGCRRyXJNzI3cGIkajD4uMU9cYnCLRjSJCybqc73uSxa3SHRcMLKJK0RZhkVAfs8fVYNN2z3TPVPTPQ3f9+vVL6rrnDp9qpj+9alTdU4pIjAzy8JW1a6AmW0+HFDMLDMOKGaWGQcUM8uMA4qZZcYBxcwy44BiZplxQDGzzDigmFlmula7ApWy4447RkNDQ7WrYVaTnnrqqSUR0ae1fFtMQGloaODJJ5+sdjXMapKkBaXk8ymPmWXGAcXMMuOAYmaZcUAxs8w4oJhZZhxQzCwzDihmlhkHFDPLjAOKmWVmi7lTthY1NDQwY8YMDj74YKZPn87cuXOZMWMG8+fPZ9ddd2XdunV07er/wo5w26JZbd72qAGjM6xJbXELpUZMnTqVGTNmVLsaRZ1xxhnsvffedO3alWnTplW7OlYl/nnLc9XNj3Vo+V875oAOLb+Q9evXl9WSKTc/wKBBgzj//PO54ooryq2ebUbcQqkR06ZNY/z48Zusu/rqq+nfvz/9+vXjggsu2CTv2LFjGT9+PPX19Vx77bU8/vjjjBgxgl69etGvXz+mTJnC2rVrN24jiUsvvZTBgwczePBgTjnlFL773e9u8nlf+tKXuPjiiwvWb+LEiRx66KH06NEjw722WuOAUsMeeOABXn75Ze655x7OO+887r333o1ps2bNYuzYsSxdupRx48bRpUsXLr74YpYsWcKf//xn7rvvPi677LJNyrv99tv561//yvPPP8/EiROZOXMmGzZsAGDJkiXce++9HHfccRXdR6stDig17Mwzz6Rbt27svffefPWrX2XmzJkb00aMGMGYMWPYaqutqKurY99992X//fena9euNDQ08PWvf52HHnpok/JOP/10evfuTV1dHZ/85Cfp2bMn9913HwA33XQTo0aNom/fvhXdR6stDig1bOedd964vMsuu/D6668XTAOYM2cORxxxBDvttBP19fVMnTqVJUuWFC0PktOYG264AYAbbriBCRMmZL0LtplxQKlhr7766sblhQsX0r9//43vJW2S96STTmKPPfbg5ZdfZvny5UyfPp3851rnbzN+/HhmzZrF7NmzeeGFFxgzZkwH7IVtThxQatjZZ5/NqlWraGxs5JprruGYY44pmnfFihXU19fTvXt3XnzxRS6//PJWyx8wYAD77bcfEyZM4Oijj6aurq5o3nXr1rFmzRo2bNjA+vXrWbNmDe+//36b9stqlwNKDTvwwAMZNGgQBx10EN/73vf4whe+UDTvBRdcwI033kiPHj048cQTWww+uSZOnMizzz7b6unOiSeeSF1dHTNnzuQnP/kJdXV1XH/99WXtj9U+5Td7N1fDhw8Pzylbvocffpjx48ezYMGCD50Sbc58p+ymJD0VEcNby+cWihW1bt06LrnkEiZNmrRFBRNrOwcUK+iFF16gV69evPHGG3znO9+pdnWsRvjWeytozz33ZOXKldWuhtWYirVQJPWWdJuklZIWSCp4y6Wk0yQ9J2mFpHmSTstLny9ptaSm9HVPZfbAzFpTyRbKpcBaoC8wDLhT0uyIaMzLJ+B44O/AR4F7JL0aETfl5DkyIu7FzDqVirRQJHUDjgbOiIimiHgUuAP40LXIiDg/Ip6OiPUR8RIwC6j8EF0zK1ulTnl2B9ZHxJycdbOBIS1tpOTSwkggvxXzG0mLJd0jaWgL20+W9KSkJxcvXtzWuptZiSoVULoDy/PWLQNaG+s+jaSO1+SsGwc0ALsADwD/LalXoY0j4sqIGB4Rw/v0afU5z2bWTpUKKE1Afd66emBFsQ0kTSHpSzk8It5rXh8Rj0XE6ohYFRHnAktJWjFmVmWVCihzgK6SBuesG8qHT2UAkHQC8APgoIhY1ErZQdKRu9lpaGjYOMfJ9OnTmTRpEgDz589HEuvXr69m9cw+pCJXeSJipaRbgbMkTSK5yjMa+HR+XknjgOnAZyNibl7aQGBn4AmSYPhNYEcgs3kb23PLdSnaelv21KlTM65Jdt5++22+/e1v89BDD7Fy5Ur22msvLrroIj71qU9Vu2pWYZW8U/ZkoA54G5gJnBQRjZJGSmrKyXcOsAPwRM69Js0TlfYALgfeBV4DDgEOjYh/VmwvalC5LZly8zc1NbHffvvx1FNP8c477zBx4kQOP/xwmpqaWt/YNisVCygR8U5EjImIbhExMCJuTNc/EhHdc/LtGhFbR0T3nNc30rTGiNgnLWOHiDgoIraIEX+deU7Z3XbbjVNPPZV+/frRpUsXJk+ezNq1a3nppZcyPgrW2XksTw3rrHPKPvPMM6xdu5ZBgwZlu8PW6Tmg1LDOOKfs8uXLmTBhAmeeeSY9e/bMfqetU3NAqWGdbU7Z1atXc+SRR7L//vtz+umnt2mfrLY5oNSwzjSn7HvvvceYMWMYMGAAv/zlL9uzW1bDPH1BDTv77LP51a9+xbx587jmmms2tiYKKTSnbGt3D5c6p+y6desYO3YsdXV1XHfddWy1Vba/U549rXY4oOSppT/A5jllN2zYUNKcspMnT+b888/nE5/4BMcccwz3339/q58xceJEJkyYwCWXXFI0z5/+9Cf++Mc/UldXR69eH4yCuPvuuxk50jcxb0k8p6y1qDPMKVuNFopbRZsqdU5Zt1CsqCznlG3PQ+h7e/KKmuFOWSvIc8paW7iFYgV5TllrC7dQzCwzDihmlhkHFDPLjAOKmWXGnbJmBTzduLDN2x41IMOK1Bi3UMwsMw4oZpYZBxQzy4wDipllxgHFzDLjgGJmmXFAMbPMOKCYWWYcUMwsMw4oZpYZ33pvm7WOfla1bcotFDPLjAOKmWWmYgFFUm9Jt0laKWmBpIIPyZV0mqTnJK2QNE/SaXnpDZIekLRK0ouSDq7MHphZayrZh3IpsBboCwwD7pQ0OyIa8/IJOB74O/BR4B5Jr0bETWn6TODPwGHp6/eSBkfE4krshFVee6YSsMqqSAtFUjfgaOCMiGiKiEeBO4APPSw3Is6PiKcjYn1EvATMAg5Iy9kd+DfgzIhYHRG3AM+mZZtZlVXqlGd3YH1EzMlZNxsY0tJGSh4GMxJobsUMAeZGxIpyyjGzyqhUQOkOLM9btwzo0cp200jqeE1OOctKLUfSZElPSnpy8WKfEZl1tEoFlCagPm9dPbCiQF4AJE0h6Us5PCLea0s5EXFlRAyPiOGtPRjczNqvUgFlDtBV0uCcdUP54FRmE5JOAH4AHBQRi3KSGoHdJOW2SIqWY2aVVZGAEhErgVuBsyR1k3QAMBq4Pj+vpHHAdODzETE3r5w5wDPAmZK2k3QUsA9wS0fvg5m1rpI3tp0M1AFvk1z6PSkiGiWNlNSUk+8cYAfgCUlN6euKnPRjgeHAu8BPgbG+ZGzWOVTsPpSIeAcYU2D9IySdrc3vd22lnPnAqIyrZ2YZ8K33ZpYZBxQzy4wDipllxgHFzDLjgGJmmXFAMbPMOKCYWWYcUMwsMw4oZpYZBxQzy4wDipllxgHFzDLjgGJmmXFAMbPMOKCYWWYcUMwsMw4oZpYZBxQzy0zJAUXSaEmVfHSpmdWYclooZwFvSPqFpE91VIXMrHaVHFAiYihwMLAauEXSS5J+KKmhg+pmZjWmrD6UiJgdEacBOwOnAP8beEXSw5LGSXKfjNkWrOw+EUkfBcanrw3Aj4CFwBTgaODfs6ygbR7m93q62lWwCig5oEg6BZgADAZuBiZExF9y0m8heYiXmW2hymmhHApcCNyR8/DyjSJilSS3Tsy2YOX0eTwYEb/LDyaSTm1ejoh7MquZmdWccgLKj4qs/2EWFTGz2tfqKY+kzzXnlfRZQDnJuwErOqJiZlZ7SulDuSr9d1vg6pz1AbwJfDPrSplZbWr1lCcido2IXYHfNC+nr90i4tMRcUcpHySpt6TbJK2UtEDScUXyfVbSA5KWSZpfIH2+pNWSmtKX+23MOoly7pQ9vp2fdSmwFugLjAMulzSkQL6VJC2h01oo68iI6J6+vtDOeplZRlo85ZH0QkTsmS6/SnKa8yERMbCVcrqR3PS2V0Q0AY9KuoPkvpYf5JX1OPC4pINL3gsz6xRa60M5MWd5fDs+Z3dgfUTMyVk3GziwjeX9Jr3N/2/AaRExux11M7OMtBhQIuLRnOWH2vE53YHleeuWAT3aUNY44GmSq03fBv5b0h4RsTQ/o6TJwGSAgQNbbESZWQbKmQ/lVEnD0uX9JS2UNE/SiBI2bwLq89bV04ZLzhHxWESsjohVEXEusBQYWSTvlRExPCKG9+nTp9yPMrMylXNj238A89Llc4GLgHOAn5Ww7RyS+1gG56wbCjSW8fnFBJveG2NmVVLOWJ6eEbFMUg+SYHBwRLwv6cLWNoyIlZJuBc6SNAkYBowGPp2fN+0b2QbYOnmr7YANEbFW0kCSqROeIAmG3wR2BB4rYz/MrIOU00J5VdKngWOBh9NgUg+8X+L2JwN1JCOSZwInRUSjpJGSmnLyfYZkEqe7gIHpcvO9Jj2Ay4F3gdeAQ4BDI+KfZeyHmXWQcloopwG/J7mX5Oh03RHA46VsHBHvAGMKrH+EpNO2+f2DFDmFiYhGYJ8y6mxmFVRyQImIu4D+eat/l77MzMqbsU1ST+Bj5LQoUvdnViMzq1nlzNj2FZLb55uAVTlJQTLq2My2cOW0UH4CjI2IuzuqMmZW28oJKF354GqLmRVx26JZbdruqAGjM65J5ZVz2fg84Id+VIaZFVNOC+U/gJ2A70va5L6P1kYbm9mWoZyA0p7Rxma2BSjnPpT2jDY2sy1AOZeNtyWZ+f7LwA4R0VPSF4DdI+IXHVVB61yuurmNw6Z6ZVsP65zK6WC9GNiLZD6S5pnbGoGTsq6UmdWmcvpQjgIGpSOHNwBExGuSPtIxVTOzWlNOC2UteQFIUh/AI33NDCgvoPwOuE7SrgCS+gG/AG7qiIqZWe0pJ6BMBeYCz5J0sb0MvAH8uAPqZWY1qJw+lEHAS8B0oAtwe0Q82yG1MrOaVMqzjUXyONKJwCLgdeAjwJmSrgdOiIiCz+sx2xI93biwTdsdNSDjilRBKac8k4FRwP4RsUtEjEhvtR9BMtv81zuwfmZWQ0oJKBOAb0XEE7kr0/ffSdPNzEoKKB8Hit12/1CabmZWUkDpEhEFH8iVrvd0BmYGlHaVZ2tJn6X4w7TKmpfWzDZfpQSDt4GrW0k3M2s9oEREQwXqYWabAfd/mFlmHFDMLDMOKGaWGQcUM8uMA4qZZaZiAUVSb0m3SVopaYGk44rk+6ykByQtkzS/QHpDmr5K0ouSDu7wyptZSSrZQrmUZNa3viTz0l4uaUiBfCtJ7ns5rUg5M4G/ATsA/wX8Pp05zsyqrCIBRVI34GjgjIhoiohHgTsoMLAwIh6PiOtJJnPKL2d34N+AMyNidUTcQjLh09EdugNmVpJK3Ta/O7A+IubkrJsNHFhmOUOAuXlji2an6z9E0mSS6RcYONAPN8zC/F5PV7sK1olV6pSnO7A8b90yoEcbyllWajkRcWVEDI+I4X36+KzIrKNVKqA0AfV56+qBgqOYK1COmXWASgWUOUBXSYNz1g0leVBYORqB3STltkjaUo6ZdYCKBJSIWAncCpwlqZukA4DRwPX5eSVtJWk7YOvkrbaTtE1azhzgGZL5bLeTdBSwD3BLJfbDzFpWycvGJwN1JNMdzAROiohGSSMlNeXk+wywGrgLGJgu35OTfiwwHHgX+CkwNiIWV6D+ZtaKik2OFBHvAGMKrH+EpLO1+f2DFJ/MiYiYTzJptpl1Mr713swy44BiZplxQDGzzDigmFlmHFDMLDMOKGaWGQcUM8uMA4qZZcYBxcwy44BiZplxQDGzzDigmFlmHFDMLDMVG21sncdVNz/W9o17ZVcP21R7/l++dswBGdak7dxCMbPMOKCYWWYcUMwsMw4oZpYZBxQzy4wDipllxgHFzDLjgGJmmfGNbWadRPseRO8b28xsM+OAYmaZcUAxs8w4oJhZZhxQzCwzDihmlpmKBRRJvSXdJmmlpAWSjiuST5LOk/TP9HWeJOWkR1pGU/qaUal9MLOWVfI+lEuBtUBfYBhwp6TZEdGYl28yMAYYCgTwP8A84IqcPEMj4h8dX2UzK0dFWiiSugFHA2dERFNEPArcAUwokH0icGFELIqI14ALga9Uop5m1j6VOuXZHVgfEXNy1s0GhhTIOyRNaynfw5LelHSrpIZiHyppsqQnJT25ePHittXczEpWqYDSHViet24Z0KNI3mV5+brn9KMcCDQAewCvA3+UVPDULSKujIjhETG8T58+7ai+mZWiUgGlCajPW1cPrCghbz3QFBEBEBEPR8TaiFgKfBvYFdgz+yqbWbkq1Sk7B+gqaXBEvJyuGwrkd8iSrhsKPN5KvmYBqIX0zVa7Zq836wAVCSgRsVLSrcBZkiaRXOUZDXy6QPZfA6dKuoskWHwX+L8AkoYAWwPPAnXAOcBrwAsdvhObkfaNarXO6LZFs9q03VEDRmdaj0re2HYySRB4G5gJnBQRjZJGSmrKyfdL4A8kQeM54M50HSSXnG8m6Y+ZS9KXckRErKvIHphZiyp2H0pEvENyf0n++kdIOmKb3wfw/fSVn/d+4GMdWE0zawffem9mmXFAMbPMOKCYWWYcUMwsM56kusraerkvsWNm9TDLglsoZpYZBxQzy4wDipllxgHFzDLjgGJmmfFVnip7unFhm7dt8FUe62TcQjGzzDigmFlmHFDMLDMOKGaWGQcUM8uMA4qZZcYBxcwy4/tQcrRn5G/Wk/2WwpNNW2fjFoqZZcYBxcwy44BiZplxQDGzzDigmFlmHFDMLDO+bGy2GWjrNBhHDci2Hm6hmFlmHFDMLDMVCyiSeku6TdJKSQskHVcknySdJ+mf6es8ScpJHybpKUmr0n+HVWofzKxllWyhXAqsBfoC44DLJQ0pkG8yMAYYCuwDHAl8HUDSNsAs4AZge+A6YFa63syqrCIBRVI34GjgjIhoiohHgTuACQWyTwQujIhFEfEacCHwlTRtFElH8s8i4r2I+Dkg4HMdvAtmVoJKtVB2B9ZHxJycdbOBQi2UIWlaoXxDgL9HROSk/71IOWZWYZW6bNwdWJ63bhnQo0jeZXn5uqf9KPlpLZWDpMkkp1AATZJeKrPepdgRWNIB5VaS96FzqPg+nMO3Ss26SymZKhVQmoD6vHX1wIoS8tYDTRERksoph4i4EriyTTUukaQnI2J4R35GR/M+dA6bwz5U6pRnDtBV0uCcdUOBxgJ5G9O0QvkagX1yr/qQdNwWKsfMKqwiASUiVgK3AmdJ6ibpAGA0cH2B7L8GTpX0EUn9ge8C16ZpDwLvA9+StK2kKen6+zuy/mZWmkpeNj4ZqAPeBmYCJ0VEo6SR6alMs18CfwCeBZ4D7kzXERFrSS4pHw8sBU4AxqTrq6VDT6kqxPvQOdT8PmjTCyZmZm3nW+/NLDMOKGaWGQeUMkl6UNIaSU3pq+i9La2NS6qGtDP7qnQ81QpJz0g6tIX8X5H0fs7+NkkaVcEq59Ylk/Fg1VLOse9Mx70cng+lbaZExIwS8uWOSwrgf4B5wBUdWLfWdAVeBQ4EFgKHAb+VtHdEzC+yzZ8j4n9VqH4tyR0PNgy4U9LsiMi/baAzHnco/9h3luNeMrdQOlZL45KqIiJWRsS0iJgfERsi4o8kX7Z9q1mv1mQ4HqxqavXYl8MBpW3OlbRE0mOtNENbGpfUKUjqSzLWqqWbAz+R7u8cSWdIqkbLNqvxYJ1GCce+Mxz3snT6CnZC/wk8T9L0Phb4g6RhEfFKgbxFxyVFJ7heL2lr4DfAdRHxYpFsDwN7AQtIvpQ3A+uBcytSyQ9kMh6sMxx3KOnYd5bjXha3UHKkHa5R5PUoQET8NSJWpNMnXAc8RnIuXEjRcUnV3Ic031YkdyqvBaYUKy8i5kbEvLSJ/ixwFjC2o+rfgkzGg3VQ3cpSyrHvRMe9LG6h5IiIUW3ZjGROlkKaxyU9nr4vNn4pM6XsQ3rF4yqSzs3DImJdOR9B8f3tSBvHg0XEy+m61saDVey4l6odx75ax708EeFXiS+gF/BFYDuSYDwOWAnsXiT/N4AXgI8A/Un+qL/RCfbjCuAvQPcS8h4K9E2X9yAZDnFmlep9E8mwjW7AASSnMkNq5biXc+w703Eva/+qXYFaegF9gCdImtlL0z+Mz+ekjyRpWje/F3A+8E76Op90uEMV92EXkl+7NSSnBs2vcWn6wPT9wPT9BcBbaeCcS9L03rpKde8N3J7WZSFwXK0c99aOfWc+7uW8PJbHzDLjTlkzy4wDipllxgHFzDLjgGJmmXFAMbPMOKCYWWYcUNpI0nxJB1e7Hp1FS8dD0ihJi9pYbpu3bcNnTZVUyrQUpPOT7NbRdUo/a5qkGyrxWe3lgFJlkk6T9Fw64c48Sae1kv9rkl5M878l6S5JBR90tqVIJ1SaK+n5Mrb5UKCKiOkRMamU7SOie0TMTcu6VtI5ZdZ5m3QkcXdtOmnXEkm3SupXTnmdhQNK9YlkFv/tgUOAKZKOLZhROhCYDnw5InoAe5KMQi3vA2tgGHyZPgP8K7CbpP2qXZkSfQZ4JiKan/gwJSK6A4NIRktfULWatYMDSgbSqf1+Jun19PUzSdvmpH9f0htp2qR05O8ggIg4PyKejoj1EfESMItknEoh+5HM4vW3dNt3IuK6iFiRfk6dpAvTKQaXSXo0XdeQfubXJC0kfY6RpN9JejPN+7CkjXOGpL+6l0m6O/3lfEzSTum+vZu2kj6RXz9Jz6fp10jarsjx6i/pFkmL01bZt3LS6tLPfjdtcZQSICamx+2udDn3s3qndXk9LfN2JZM13Q301wfTK/bPPbVI93tKXlmzJf17uhySBil53O044PtpOX9IW5235G37c0mX5Kw6LK3vJiJiKcnwgmF5SdtI+nXaMm2UtPEJg5L2TFs5S9O0L5VwzDpGte/9r9UXMB84OF0+i2Rcz7+SjPf5E3B2mnYI8CbJnBb/AtxAMp5jUIEyBfyNIgPZSMasrAZ+TBJ0ts1Lv5TkYWgfAboAnwa2BRrSz/w1ycC6ujT/CSTziWwL/IzkF7O5rGtJnrO7L8lgyPtJZhc7Pi37HOCBvOPxHLAzyZibx4Bz0rRRwKJ0eSvgKeBHwDbAbiRjVb6Ypv8UeCQtY+e0zEUt/D/8C8k8KYeRzOi2BNgmJ/1Oklbc9sDWwIH5dcrJOw24IV0+HngsJ+3jJOO3tk3fb/w/TI/VOTl5+5GMwemVvu9K8jyqfXPyvAh8LF1+EJiULu8A3AvMyqvXmnQfu5DMifKXNG1r4B/A1PR4fo5krNnHqvK9qPYXs1ZfbBpQXiEZit6c9kVgfrp8NXBuTtogigeUH5PMLrZtC597KMmD0JaSDCa7KP0j24ok2AwtsE1D+pm7tVBurzRPz/T9tcCvctK/CbyQ835vYGne8fhGzvvDgFfS5Y1fXuBTwMK8zz4duCZdngsckpM2Of+Ln7fteGBx+qXdjmQE8lFpWj9gA7B9ge1G5ZfLpgGlRxoUdknf/wS4Oidv0YCSrrsbODFdPgJ4Pifto8A/ct4/CKxK6x7AM6SDBHPqdW/O+48Dq9PlkSQ/WFvlpM8EplXje+FTnmz0J5lZq9mCdF1z2qs5abnLG6XN6+OBwyPivWIfFBF3R8SRJL/go0nmSp0E7EjyhSo0c9yHPltSF0k/lfSKpOUkAYG0nGZv5SyvLvC+e7Hy2fQY5NqF5FRjafOL5Ne1b5qef7wW5BeQZyLw20hOGdcAt/DBac/OwDsR8W4rZXxIJKeRd5LMygfwZZIZ1kp1HUmwI/0397G7h5EEnFzfioieJM/q3h4YkJf+Zs7yKmC7tC+sP/BqRGzISV9A0kqtOAeUbLxO8kVpNjBdB/AGm/5x7Jy/saQTgB8AB0VESZdII5nJ6z6SU5G9SJr6a0h+/YpulrN8HElAOhjoSdKKgfZN4pO7b7nHINerwLyI6JXz6hERzbPevVGgnIIkDSBp4o9P+4LeJJnV7DBJO6af1VtSrwKblzLMfibwZUkjSIL1A0XyFSrrdmAfSXuRtFByg1HB/hOASGZnOwe4VCrp0R+vAzsrmQWu2UDgtRK2zZwDSjZmAj+U1Cf9Q/4RSV8JwG+Br6YdZ/8CnJG7oaRxJFduPh/pZchiJI2WdKyk7ZX4JMkjGf6S/kJdDVyUdjB2kTQit3M4Tw/gPeCfJP0Q09u055s6RdIASb2B/6LwFajHgRWS/jPtgO0iaS99cHXmt8Dp6T4OIDnVKmYCyUxuHyPpxBxGMunzIpIrYW+QtAQuS8vbWtJn0m3fAnaQ1LOF8u8i+aE4C7g5rxWQ6y2SvqCN0tbS74EbgccjYiFA+jfwSYoHJ0haN32BUjpX/0rSYvl+un+jgCNJJqOqOAeUbJwDPAn8neQh70+n64iIu4Gfk/wB/YOk8xaSL3PztjsAT+RccSj2/Jh3gROBl0k6Im8A/k9ENP/6fS/9/CdIJhY6j+L/x78maRq/RjLp9l+K5CvHjcA9JP0gr5Aeg1wR8T7JL/Ywkk7eJcAMklYSJP1IC9K0e9j0VCHfROCyiHgz90UyK1rzac8EYB1JJ+jbwHfSerxI8kMwNz31+tDpWXrqeStJK+7GFupxFfDxtJzbc9ZfR9LXlLsPnyO5UremWGERsRa4hLwfnxbyHknSt7YEuAw4PopPOt6hPMFShUnak+TKxbYRsb7a9bGOI2kgSSDbKSKWp+suA56LiMuqWrkO4hZKBUg6Ssm9KtuTtBr+4GCyeUv7NE4FbmoOJqlngNuqU6uO5xZKBUj6f8AI4H3gIeDk9PzeNkPpjXNvkZy6HRIRBa/sbY4cUMwsMz7lMbPMOKCYWWYcUMwsMw4oZpaaQt7NAAAAEUlEQVQZBxQzy4wDipll5v8Dxi8os0YOkuIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rho_dir = os.path.join(data_dir, \"Rhodopsin\")\n",
    "\n",
    "# Read in data\n",
    "library1_df = pd.read_csv(os.path.join(rho_dir, \"library1TotalExpressionSummary.txt\"), sep=\"\\t\", index_col=0)\n",
    "library1_df[\"library\"] = 1\n",
    "library2_df = pd.read_csv(os.path.join(rho_dir, \"library2TotalExpressionSummary.txt\"), sep=\"\\t\", index_col=0)\n",
    "library2_df[\"library\"] = 2\n",
    "\n",
    "# Get scrambled sequences from each library with RNA barcodes measured\n",
    "scrambled_library1_df = library1_df[library1_df.index.str.contains(\"scrambled\") & (library1_df[\"expression\"] > 0)]\n",
    "scrambled_library2_df = library2_df[library2_df.index.str.contains(\"scrambled\") & (library2_df[\"expression\"] > 0)]\n",
    "\n",
    "# Compare distributions of log2 expression\n",
    "scrambled_library1_expr = np.log2(scrambled_library1_df[\"expression\"])\n",
    "scrambled_library2_expr = np.log2(scrambled_library2_df[\"expression\"])\n",
    "ks_stat, pval = stats.ks_2samp(scrambled_library1_expr, scrambled_library2_expr)\n",
    "print(f\"Scrambled sequences from L1 and L2 are drawn from the same distribution, KS test p = {pval:.3f}, D = {ks_stat:.2f}\")\n",
    "\n",
    "# Show the two histograms\n",
    "fig, ax = plt.subplots()\n",
    "ax.hist([scrambled_library2_expr, scrambled_library1_expr], bins=\"auto\", histtype=\"stepfilled\", density=True, label=[\"library 2\", \"library 1\"], color=plot_utils.set_color([0.75, 0.25]), alpha=0.5)\n",
    "ax.set_xlabel(\"log2 Scrambled Activity/Rho\")\n",
    "ax.set_ylabel(\"Density\")\n",
    "ax.legend(loc=\"upper left\", frameon=False)\n",
    "plot_utils.save_fig(fig, os.path.join(figures_dir, \"supplementalFigure2\"), timestamp=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Join two Rho libraries together, annotate sequences as strong/weak enhancer, inactive, silencer, or ambiguous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cutoff to call something a strong enhancer: activity is above 2.84\n"
     ]
    }
   ],
   "source": [
    "color_mapping = {\n",
    "    \"Strong enhancer\": \"#1f78b4\",\n",
    "    \"Weak enhancer\": \"#a6cee3\",\n",
    "    \"Inactive\": \"#33a02c\",\n",
    "    \"Silencer\": \"#e31a1c\",\n",
    "    np.nan: \"grey\"\n",
    "}\n",
    "\n",
    "# Join the libraries and add a pseudocount to take log2\n",
    "rho_df = library1_df.append(library2_df)\n",
    "rho_pseudocount = 1e-3\n",
    "rho_df[\"expression_log2\"] = np.log2(rho_df[\"expression\"] + rho_pseudocount)\n",
    "\n",
    "# Define cutoff for a strong enhancer based on scrambled sequences\n",
    "scrambled_mask = rho_df.index.str.contains(\"scrambled\")\n",
    "scrambled_df = rho_df[scrambled_mask]\n",
    "scrambled_df = scrambled_df[scrambled_df[\"expression\"].notna()]\n",
    "strong_cutoff = scrambled_df[\"expression_log2\"].quantile(0.95)\n",
    "print(f\"Cutoff to call something a strong enhancer: activity is above {strong_cutoff:.2f}\")\n",
    "\n",
    "# Drop scrambled sequences\n",
    "rho_df = rho_df[~scrambled_mask]\n",
    "\n",
    "# Helper function to label and color a sequence\n",
    "def label_color_sequence(row, alpha=0.05, strong_cutoff=strong_cutoff, inactive_cutoff=1, color_mapping=color_mapping):\n",
    "    expr_log2 = row[\"expression_log2\"]\n",
    "    qval = row[\"expression_qvalue\"]\n",
    "    # Inactive\n",
    "    if (np.abs(expr_log2) <= inactive_cutoff) & (qval >= alpha):\n",
    "        group = \"Inactive\"\n",
    "    # Silencer\n",
    "    elif (expr_log2 < -inactive_cutoff) & ((qval < alpha) | (row[\"expression\"] == 0)):\n",
    "        group = \"Silencer\"\n",
    "    # Enhancer\n",
    "    elif (expr_log2 > inactive_cutoff) & (qval < alpha):\n",
    "        # Strong\n",
    "        if expr_log2 > strong_cutoff:\n",
    "            group = \"Strong enhancer\"\n",
    "        # Weak\n",
    "        else:\n",
    "            group = \"Weak enhancer\"\n",
    "    # Ambiguous\n",
    "    else:\n",
    "        group = np.nan\n",
    "    \n",
    "    color = color_mapping[group]\n",
    "    return pd.Series({\"group_name\": group, \"plot_color\": color})\n",
    "\n",
    "# Annotate both WT and MUT sequences\n",
    "rho_df = rho_df.join(rho_df.apply(label_color_sequence, axis=1))\n",
    "rho_df[\"group_name\"] = sequence_annotation_processing.to_categorical(rho_df[\"group_name\"])\n",
    "\n",
    "# Save to file\n",
    "sequence_annotation_processing.save_df(rho_df, os.path.join(data_dir, \"rhoActivityJoinedAnnotated.txt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in Polylinker data, annotate for autonomous activity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "polylinker_dir = os.path.join(data_dir, \"Polylinker\")\n",
    "\n",
    "# Read in data\n",
    "library1_poly_df = pd.read_csv(os.path.join(polylinker_dir, \"library1TotalExpressionSummary.txt\"), sep=\"\\t\", index_col=0)\n",
    "library2_poly_df = pd.read_csv(os.path.join(polylinker_dir, \"library2TotalExpressionSummary.txt\"), sep=\"\\t\", index_col=0)\n",
    "\n",
    "# Join two dataframes, compute log2\n",
    "poly_df = library1_poly_df.append(library2_poly_df)\n",
    "poly_pseudocount = 1e-2\n",
    "poly_df[\"expression_log2\"] = np.log2(poly_df[\"expression\"] + poly_pseudocount)\n",
    "\n",
    "# Annotate for autonomous activity\n",
    "poly_df[\"autonomous_activity\"] = (poly_df[\"expression_log2\"] > 0)\n",
    "\n",
    "# Save to file\n",
    "sequence_annotation_processing.save_df(poly_df, os.path.join(data_dir, \"polyActivityJoinedAnnotated.txt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute WT vs. MUT activity in the Rho assay, join with WT Polylinker activity\n",
    "Only compute effect sizes when the MUT is CRX sites (i.e. not DNA shape mutants)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ryan/miniconda/envs/bclab/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in greater\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/home/ryan/miniconda/envs/bclab/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:879: RuntimeWarning: invalid value encountered in less\n",
      "  return (self.a < x) & (x < self.b)\n",
      "/home/ryan/miniconda/envs/bclab/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:1821: RuntimeWarning: invalid value encountered in less_equal\n",
      "  cond2 = cond0 & (x <= self.a)\n"
     ]
    }
   ],
   "source": [
    "wt_mask = rho_df.index.str.contains(\"_WT$\")\n",
    "mut_mask = rho_df.index.str.contains(\"_MUT-allCrxSites$\")\n",
    "\n",
    "# Add variant info as a column, then trim it off the index\n",
    "rho_df_no_variant_df = rho_df.copy()\n",
    "rho_df_no_variant_df[\"variant\"] = rho_df_no_variant_df.index.str.split(\"_\").str[2:].str.join(\"_\")\n",
    "rho_df_no_variant_df = sequence_annotation_processing.remove_mutations_from_seq_id(rho_df_no_variant_df)\n",
    "\n",
    "# Separate out WT and MUT, then join them together on the same row\n",
    "wt_df = rho_df_no_variant_df[wt_mask]\n",
    "mut_df = rho_df_no_variant_df[mut_mask]\n",
    "wt_vs_mut_rho_df = wt_df.join(mut_df, lsuffix=\"_WT\", rsuffix=\"_MUT\")\n",
    "wt_vs_mut_rho_df[\"wt_vs_mut_log2\"] = wt_vs_mut_rho_df[\"expression_log2_WT\"] - wt_vs_mut_rho_df[\"expression_log2_MUT\"]\n",
    "\n",
    "# Compute parameters for lognormal distribution to do stats\n",
    "wt_cov = wt_vs_mut_rho_df[\"expression_std_WT\"] / wt_vs_mut_rho_df[\"expression_WT\"]\n",
    "wt_log_mean = np.log(wt_vs_mut_rho_df[\"expression_WT\"] / np.sqrt(wt_cov**2 + 1))\n",
    "wt_log_std = np.sqrt(np.log(wt_cov**2 + 1))\n",
    "mut_cov = wt_vs_mut_rho_df[\"expression_std_MUT\"] / wt_vs_mut_rho_df[\"expression_MUT\"]\n",
    "mut_log_mean = np.log(wt_vs_mut_rho_df[\"expression_MUT\"] / np.sqrt(mut_cov**2 + 1))\n",
    "mut_log_std = np.sqrt(np.log(mut_cov**2 + 1))\n",
    "\n",
    "# Do t-tests and FDR\n",
    "wt_vs_mut_rho_df[\"wt_vs_mut_pvalue\"] = stats.ttest_ind_from_stats(wt_log_mean, wt_log_std, wt_vs_mut_rho_df[\"expression_reps_WT\"], mut_log_mean, mut_log_std, wt_vs_mut_rho_df[\"expression_reps_MUT\"], equal_var=False)[1]\n",
    "wt_vs_mut_rho_df[\"wt_vs_mut_qvalue\"] = modeling.fdr(wt_vs_mut_rho_df[\"wt_vs_mut_pvalue\"])\n",
    "\n",
    "# Pull out WT polylinker measurements\n",
    "poly_wt_df = poly_df.copy()\n",
    "poly_wt_df = poly_wt_df[poly_wt_df.index.str.contains(\"WT\")]\n",
    "\n",
    "# Drop the variant ID\n",
    "poly_wt_df = poly_wt_df.rename(index=lambda x: x[:-3], columns={\"expression\": \"expression_POLY\", \"expression_SEM\": \"expression_SEM_POLY\", \"expression_log2\": \"expression_log2_POLY\"})\n",
    "\n",
    "# Join with Rho\n",
    "activity_df = wt_vs_mut_rho_df.join(poly_wt_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annotate for binding patterns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get info on CRX binding from the seq ID strings\n",
    "activity_df[\"crx_bound\"] = activity_df.index.str.contains(\"_C...$\")\n",
    "\n",
    "# Read in BED files\n",
    "library_bed = BedTool(os.path.join(data_dir, \"library1And2.bed\"))\n",
    "nrl_chip_bed = BedTool(os.path.join(\"Data\", \"Downloaded\", \"ChIP\", \"nrlPeaksMm10.bed\"))\n",
    "mef2d_chip_bed = BedTool(os.path.join(\"Data\", \"Downloaded\", \"ChIP\", \"mef2dPeaksMm10.bed\"))\n",
    "\n",
    "# Get binding patterns for NRL and MEF2D\n",
    "library_nrl_bound_df = library_bed.intersect(nrl_chip_bed, wa=True).to_dataframe()\n",
    "activity_df[\"nrl_bound\"] = activity_df.index.isin(library_nrl_bound_df[\"name\"])\n",
    "\n",
    "library_mef2d_bound_df = library_bed.intersect(mef2d_chip_bed, wa=True).to_dataframe()\n",
    "activity_df[\"mef2d_bound\"] = activity_df.index.isin(library_mef2d_bound_df[\"name\"])\n",
    "\n",
    "# Helper function to \"reverse one hot encode\" binding patterns\n",
    "def annotate_binding(row):\n",
    "    crx, nrl, mef2d = row[[\"crx_bound\", \"nrl_bound\", \"mef2d_bound\"]]\n",
    "    if crx:\n",
    "        if nrl:\n",
    "            if mef2d:\n",
    "                return \"All three\"\n",
    "            else:\n",
    "                return \"CRX+NRL\"\n",
    "        elif mef2d:\n",
    "            return \"CRX+MEF2D\"\n",
    "        else:\n",
    "            return \"CRX only\"\n",
    "    elif nrl:\n",
    "        if mef2d:\n",
    "            return \"NRL+MEF2D\"\n",
    "        else:\n",
    "            return \"NRL only\"\n",
    "    elif mef2d:\n",
    "        return \"MEF2D only\"\n",
    "    else:\n",
    "        return \"No binding\"\n",
    "\n",
    "activity_df[\"binding_group\"] = activity_df.apply(annotate_binding, axis=1)\n",
    "\n",
    "# Write to file\n",
    "sequence_annotation_processing.save_df(activity_df, os.path.join(data_dir, \"wildtypeMutantPolylinkerActivityAnnotated.txt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write different groups to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_dir = os.path.join(data_dir, \"ActivityBins\")\n",
    "\n",
    "masks_variant_names = [\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Strong\"), \"variant_WT\", \"strongEnhancer.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Weak\"), \"variant_WT\", \"weakEnhancer.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Inactive\"), \"variant_WT\", \"inactive.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Silencer\"), \"variant_WT\", \"silencer.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Strong\") & activity_df[\"group_name_MUT\"].str.contains(\"Strong\"), \"variant_MUT\", \"strongToStrong.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Strong\") & activity_df[\"group_name_MUT\"].str.contains(\"Weak\"), \"variant_MUT\", \"strongToWeak.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Strong\") & activity_df[\"autonomous_activity\"], \"variant_WT\", \"strongAutonomousEnhancer.fasta\"),\n",
    "    (activity_df[\"group_name_WT\"].str.contains(\"Strong\") & ~(activity_df[\"autonomous_activity\"]), \"variant_WT\", \"strongNonautonomousEnhancer.fasta\")\n",
    "]\n",
    "\n",
    "for mask, variants, name in masks_variant_names:\n",
    "    mask = mask & mask.notna()\n",
    "    ids = activity_df.loc[mask, variants]\n",
    "    fasta_seq_parse_manip.write_fasta(all_seqs[ids.index + \"_\" + ids], os.path.join(out_dir, name))"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}