13_plot_ziegenhain_Fig4_Fig5.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle\n",
    "\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "from readcount_tools import broken_zeta, add_largedot_legend\n",
    "#import plotting constants\n",
    "from readcount_tools import PAPER_CONTEXT, PAGEWIDTH_IN, POINTSIZE_HIGHLIGHT,POINTSIZE_SMALL,SMALL_SIZE,MEDIUM_SIZE,BIGGER_SIZE,LINEWIDTH,POINTSIZE,TICKLENGTH,LEGEND_FONTSIZE,SPINEWIDTH,LETTER_LOC_X,LETTER_LOC_Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Load preproc results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "basepath = 'data/reads_per_umi_tables/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "tcmalloc: large alloc 1073741824 bytes == 0x7a8d6000 @ \n",
      "tcmalloc: large alloc 2147483648 bytes == 0xe08de000 @ \n",
      "tcmalloc: large alloc 4294967296 bytes == 0x1b86a2000 @ \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 59.4 s, sys: 30.9 s, total: 1min 30s\n",
      "Wall time: 1min 35s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#1min30sec\n",
    "with open(f'{basepath}reads_per_umi_dfs_hd1_preproc.pickle','rb') as f:\n",
    "    dfs = pickle.load(f) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Prepare Fig3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "exps = [df.index[0] for df in dfs]\n",
    "protocols = [e.replace('_Fibroblast_PE','P').replace('_Fibroblast','S') for e in exps]\n",
    "runs = [p[-1] for p in protocols]\n",
    "runs = [r.replace('P','PE') for r in runs]\n",
    "runs = [r.replace('S','SE') for r in runs]\n",
    "protocols = [p[:-1] for p in protocols]\n",
    "colormap = {'CELseq2':'tab:blue', \n",
    "            'DropSeq':'tab:orange', \n",
    "            'MARSseq':'tab:green', \n",
    "            'SCRBseq':'tab:red', \n",
    "            'Smartseq3':'tab:brown',\n",
    "            '':'tab:olive',}\n",
    "prot2print = dict(CELseq2='CEL-seq2',\n",
    "                    DropSeq='Drop-seq',\n",
    "                    MARSseq='MARS-seq',\n",
    "                    SCRBseq='SCRB-seq',\n",
    "                    Smartseq3='Smart-seq3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Figure 3: Paper version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAADpCAYAAACk03fwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9OUlEQVR4nO3de1xUdf748ddwE7xmKKhg4rW84KV0s9yMNru4imGSlKZQNF8x+blqtZGoM5iEl7VcTXMdXcmypETJSxd1N7K0Mg0Vb6UpraCGkikXEZDz+2OcURSGizNzZob38/HwIQznzLzhMzPv+VzO+6NRFEVBCCGEEJVyUzsAIYQQwpFJohRCCCEskEQphBBCWCCJUgghhLBAEqUQQghhgSRKIYQQwgIPtQMQdXfw4EESExNp0aIFgYGBtGvXji+++AJvb28CAgIYMGAAOp2OTp060axZM9588030ej0///wzTZo04Y477iA+Pl7tX6NeSU5OJi4ujjNnzqgdilWUlJTw0ksvoSgKJSUlbNmyhUGDBuHh4cFDDz2Ev78/Bw4cIDY29pYf65VXXuHixYv88ccfLF++nCZNmljhNxCWJCcns3btWjp27Eh+fj7Tpk2jQ4cOACiKgkajUTlC+6g2UaalpTF16lQuXrzIkCFD+Ne//mWPuEQNbNmyhTFjxjB48GDy8vKYMmUK7777LmB8A9u5cydPPfUUsbGxPP300+bzpk6dSo8ePXj22Wcr3F96ejr//ve/8ff3Z8qUKRw4cIDNmzdz6dIlRowYwR133IFer6dLly7s2rWLzz//3K6/r6u4dOkSjz32GPv37yctLY17771X7ZDqzGAw8Ne//pXBgwcDMGrUKBYsWEDjxo0B43OqMmfPnmXChAkEBQUxdOhQ2rVrx/z581EUhY4dOxITE0N0dDRBQUF89tlnbN++nXnz5gHw1ltvsXfvXh544AG7/I41dejQIcaMGUNOTg6RkZHMmTNH7ZCsIiYmhqFDh5KXl8fkyZP55ZdfCA0N5eGHH+af//wn3t7ehIaG0qtXL0aPHs3IkSP56aefePvtt3FzMw5aOvt7S7VDrxqNhitXrqAoCsuWLePcuXP2iEvUQHR0NN988w3R0dEsX76c4OBg88+8vLwASE1N5fHHH6dHjx7mn82ZM4dnn322wm0AJ0+epF27dmi1Wlq3bs3ChQu57bbbaN26Nbt27cJgMDBnzhxmzJhBaWmpfX5JF1RSUsLatWt5+OGHzW/+zurgwYP069fP/L2XlxeTJk0iJiaGL774osrzLly4gLu7O2FhYTzwwAMsWbIEHx8ffH19yczMZOvWrTz66KMkJibSunVr83lnzpxh9+7d3H///Tb9vepi5syZ3HHHHZw5c4Zp06apHY7V+fr6UlpaSrdu3YiLiyMlJYXXX3+d5cuX8/777wPQtWtX/va3v9G5c2e+//5787nO/t5SbY9y2bJlPPjggzzzzDP85S9/oaCggBYtWtgjNlGNpk2bkpiYCMDgwYNp2bKl+WclJSUAjBgxgvHjxxMZGcnly5cBePXVV+nRowfjx4/n7NmzJCYmEhAQwCuvvMLRo0dZsGAB4eHhlJeXM23aNDw8jE+Tl156CTB+eKovQy624ip/w+7du7Nnzx4ee+wxwPi8e/vtt6vsUS5YsICsrCySkpJYuHAh69atY+vWrZSXlzNmzBh69uwJwMaNG83nmP5OOTk5xMXFsXjxYtzd3e3w29WeqQdl+t+V5OXl4eXlRbNmzYCKQ6+m/8vKygAoLS2lrKyMSZMmucR7S7WJcsSIEUyZMoWzZ8/aIx5RC2lpaXzxxRd4eHjQo0cP+vTpQ3R0NN7e3gQGBnLfffcB4O7uTkREBMuWLQPgjTfe4Pbbb6ekpARfX18WLFgAwNq1a/nuu+/4448/CAwMZOLEibzwwgvcfvvt9O3blxdeeIG4uDg6duxofoKL2vPy8mL48OFkZmayfv16tcO5JVqtlilTprBx40auXLnC6dOnmTRpEh4eHtx77720b9+elJQUDhw4QLNmzczDkZmZmaxYsYLi4mIGDRpEv379mDp1Kq1bt6ZJkya8+uqrvPDCC/z000+cOnUKgGHDhtGlSxfi4uIYP348vXr1UvNXv8m0adMYM2YMrVq1YuzYscydO1ftkKxi6dKlbN26lfz8fPR6PYsXLwZg3LhxTJ8+nYYNG/LMM88AcPToUaZOncrp06d5+eWXzcPjzv7eopFar6IuwsPDWbt2rdphiHogKiqqQi9VOKasrCzefvtt/vGPf9zS/Tjie4skSiGEEMIC1xtIF0IIIaxIEqUQQghhgSRKIYQQwgJJlEIIIYQFTpMok5OT1Q5BCNVs376d2bNns2rVqmqPldeKc5G2dXyOd8FKFbKystQOQQibOX78OImJiVy4cIG1a9dSWFjIiy++iJeXFyEhIWRnZxMXF1ejsmjyWnEs0rbOz2kSpRCurEOHDqxYsYLw8HAA1q1bR3h4OKGhoURERHDPPfdYPD89Pd1cBUfeTB2LtK3zk0QphAPKzs421+51d3enf//+zJ49m4CAgEqPDwkJISQkBAC9Xm+nKEVdSNs6H9US5fbt29m5cydt2rRh7NixaoUhhEMKDAwkOzub3r17U15ezsCBAxk4cKDaYQkrkLZ1PjZbzHP8+HGio6PNww2FhYVERkai1WpZvXo13377LXFxcZw+fdpWIQhRY1cKCijKyOBKQUGNjp87dy6rV6+22uPn5eURExNDRkYGSUlJPPnkk6SmpjJ+/HhCQ0Ot9jj1SX5xKXt+PU9+cc12o7B2m5rYo21/Swyu/iBXUnwRTu4y/l8Dt9q2NutR3uq4PMjYvLCPCxs3clqnRykqQtOwIa0T9DSz8AZWWFhIixYtrLpRgK+vL0uXLq1w28qVK612/wBhi3eQNmGAVe/TUaVl5DB1fSZFJVdo6OXOG8ODCetT+dAm2KZNTezRtv7xmVa9P4e2/yPYOAlKC8GzEYQugJ4jqzzcGm1rt8tDsrOzadu2LVBxXL5NmzZVnhMSEoJer0ev1xMUFGSnSEV9cqWgwJwkAZSiIs7o9BZ7luvXrycnJ4dDhw5x9OhRe4V6y+pLkswvLjUnSYCikivEr8+02LN01jatd4ovXkuSYPx/02SLPUtrtK3dEqVpXB4wj8vHxcUxZswYe4UgxE0uHz1qTpIm5UVFlBw7VuU5R44cYfr06SxcuJCUlBRbh2g1BoNB7RDs4uffCsxJ0qSw5ApHc6v+8OOsbVrvnD1yLUmalBTA2Z+qPMUabWuzRClzLsIZNOjcGU3DhhVuc2vYEK9Onao8Z9asWQB4e3s71U72Wq1W7RDsoot/Yxp6VdzYuZGXO539qt6my1nbtN5peZdxuPV6Xo2h5Z1VnmKNtrXZHKU9xuWrYjAY6s2bgrg17o0b0zpBbx5+dWvYkFYJetxl70On1cTbkzeGB5uHXxt5uZM4PJgm3p5qhyZulXdT45ykafjVqzEMfct4uw255HWUg7ZtA0mUooaahYbS+KGHKDl2DK9OnepVknTVBT5hfQJ4uKsfR3ML6OzXuH4mScPDoP2P2lFYX8+R0OVx43BryzttniTBiWq91kZ7mWMQtVTq4cE5X19KPVzys2OVhngdUjsEm/HkCi01BXhypfqDXZCBp9UOwWaK8eIkrSjGyy6PV7/eFYSoxP79+9m4cSOlpaV4enoSGhpKz549qzxer9fTtGlTDhw4wMKFC2nsxD1QV52iqG2bgmu1K1TSti7Sw1SjbV2yRylETRUXF5tfdAClpaVs2rSJ4uJii+cpikJhYSENb1gI5OxcYWVsXdsUXLdd4boepuFhdQO5BWq1rSRKUa+dPXvW/KIzKSkpqfbi5L/97W+MHDmSvXv32jA6+9tc0g24ljCdMXHWtU3BddsVrvUwnTlhqtW2MvQq6rWWLVvi6elZ4cXn5eVFy5YtLZ63cOFCfvnlF15//XVbh2hXpoU9pjdVrVZrTpamr3Nzc/Hz86vwP4Cfn59DDOXWtU3Bddv1etcnTK2TDceq1bYaRVGUWp+lAlOFHiGs7fo5Dy8vL4YOHVrtnIcjs8ZrpS6rYRMTEwGIj4+/pce2BldrUxNrvw8aDAaG5S7E38/PaRKmGm0rPUpR7/Xs2ZMuXbpw9uxZWrZsibe3t9ohqc6UJMMW76jwvSXx8fEYDAYSExNVT5bSpjVj7F0aRwqGJQYbEyY4dNJUo21ljlIIjFU72rZtK2+oN0ibMIC0CQPMCdP0f1W0Wi1+fn4kJCTYIzyLpE1rTqvV4h+fiYGnjXOYiVXX4HYE9m5b6VEKIapl6lGmTRhAtxmf08W/CT//lm/+H6C41HS9YhugDSvjNgEa8330bnubSxY3cCWm+cvExFxeSAzGnzzjD+JPqRiV+iRRCiFq5dDMx2t8rF6vJyAgAK1WS1DcZhtGJazJOHQeT2JiIlFl7xGQ2Ab8ujr0kKwtydCrEMJm9Ho9OTk5JCQkkDV7iCRLJxMfH8+nbV4ikQmU5OyDhOZOeVnJrar3ifJERITaIQjh0kyrNBMSEnDTQLcZn6sbkKgVrVZLfHw87wYkkcAkynP2GBNmPVLvE6XUhRXC9nQ6HYqiEN147017RQrnoNVq0el0rAiYR7bih6JvBvpmaodlF/U+UQoh7EOv11+9UFwhKG5ztStohWPSarUE6n/iDc+/k02repEw632idMYSXUI4K71eT5T3blpqCth78g+Zs3Ri8fHxBOp/YqZmCuVojAnTRecv632idISSW0LUJ3q9niENDvOcz24ASZZOTqfT4ab/gwQmo+TsdsneZb1PlEII+9Pr9SiKgrbJXkAW+LgCvV6PRn+By3i63HCsJEohhCpMc5YtNQWywMeFNNCfM/YuwZgwXYAkSiGEakzDsKYFPsI1mHqXOVcX+zh7wpREKYRQlWmBjyRL12Na7AM49XCsJEoLpBiBcDRJSUmsWbNG7TCsTpKl67atTqe71rvkasJ08KLrN5JEaYEUIxC2cvz4caKjowkPDwegsLCQyMhItFotq1evrvScr776iuDgYHuGaVfXkmX1u5Q4MmnbylXoXZYWOtVwrCRKIVTQoUMHVqxYYf5+3bp1hIeHYzAY2LBhA2lpaURFRREVFUXK1Q9se/bsYdeuXXz33XdqhW1zAQEBRHn/wN6T5502WUrbVk2n0/GG599JYDKA01x7KbuHWGAwGOQ6S2EX2dnZ5h6Fu7s7YWFhhIWFVThmypQpZGVlVfpmmp6eTnp6OgBZWVk2jtZ2tFqteS9LU7J09q25pG0rMm3qrdeDjrfM115q9BfUDcwCSZQWSJIU9hIYGEh2dja9e/emvLy8yuOCgoIICgq66faQkBBCQkKAa0XInZVOp0PR60ku7svek3+oHc4tk7atnPF30XNZ3wIvSlH0zcj1vAP/+Ey1Q7uJDL0KoYK8vDxiYmLIyMggKSmJJ598ktTUVMaPH09oaKja4amu4uKeTWqHUyvStrWzKiCJHFoB4Ff6Py7rW6gc0c00iqIoagdRE3q93qU+TQlhK670Wun7WgrnlEY09PKo1YbRrsqV2vZGer0eHW8BkEMrPgt42WFG9aRHWQtyuYgQ9rU7yfiaKyop437depWjEbbk6elp7lkGcIacnByH2bRCEmUtyOUiQthf1uyhAJy67EmnuA0qRyNsxbQbSQmegHGhT05OjkP0oCVR1oKjfLoRor7p3bY5AGW4SbJ0cQ305zilaQ1gHoo1rYRWiyTKWnCU8XIh6pu0CQPQXP26DDf6viajO64sQHfE/PULfIiiKKr2LCVRCiGcwomrQ7AA55RG/GnqRypGI2xNE9AXgDb8Zr4tMTFRlVgkUQohnIan5toi/dzyhtKzdGXa/6ABNFxr89LSUlVCUTVRbt68mYSEBN5//301w6iz7x96SO0QhKhXjiaF0tDrWp2Uc0ojWTvgyjTGFGWaqwR1ii5YLVHWpRDwkCFDePnll8nOzrZWGHZ175dfqh2CEPVOcekVMM9YQuIvrVUbkhM21uZuc0u/wIeqhWG1RFmXQsDl5eXMmzePmJgYa4WhGvlUK4R9HE8actNthvxe9H0tRRKmq9H+x/xlAGfMydLevUqb1XqtSSHguXPncuHCBb7//nsee+yxm+7DmYoBy4pYIewna/aQq/tWXpu/+kPxobS0VDYzcGEtyVPlcW02R2kqBAxUWQg4Li6Ot956q9IkCcZiwKaSTZUVC3Ykpqo90rMUwj7cNHD9EGwZbiQX9yUnJ0e1mIQN6C+YW9mLUl7jbePNduxVWq1HmZeXR3x8vLkQ8MSJE4mNjWXz5s31ohCwqWqPfJIVwj68Pd0pKrly0+3JxX2Nezhh3N9SXpOuwZQsG1CKPuDrCsOytma1ROnr68vSpUsr3LZy5Upr3b3TOBERYU6aMgQkhO108W/Cz7/lU1RSdtPPkouN1+A9d2oPIK9Fp6fyXpWyH6WVXV8PVqvVygtUCBu5fkPnhIQEKtsISVGMQ3Senp7yWnRm+mbc2LoKGmYyCU9PT/z8/MjNzcXPz88mbSyJ0oZORESgld6lEDYTtnjH1c2d+8JNb6XXKQby4V+vpRAYGFghyQoncN08pYkG0Nvp4SVR2tCNvUuQhCmENUnCqycMD0PO7sp/FtDX5vOVkiht7Po5S5DhWCGEqDU7LtypjCRKGzMlyesTpiRLIazPYDDU6dIQ08pYg8FAbm4u8fHxJCYmEh8fb4MohTOSRGknNyZMGYoVwrpu9XV0/fmSJMX1LBYcKC4utlcc9Ub7lBRzcQK41rsUQgjhmCwmypiYGF577TV+/fVXe8VTL1SVLCVhCiGE47GYKJOTkxk3bhwTJkxg+PDh7N5dxaojUWvXD8WCMVlK71IIIRyPxUT58ccfM2PGDAYNGsQ777zDP//5T3vFVW9U1ruUHRCEEMJxWEyU7u7uJCcnM2nSJFq1asW8efPsFVe9cmOyjI+Pl56lEEI4CIuJ8uTJk7i5GQ9ZsmQJrVq1sktQ9ZEpWV4/FGswGKR3KYQQKrN4ecj//vc/89eyoMf2TPOWR+6+h7t+3FPhEhKQnUnqu/Xr1/Prr79y8eJFZsyYoXY4woqkbR1btUOv06dPZ8aMGVXuKSms764f99zUuwTZ69KVHD9+nOjoaMLDwwEoLCwkMjISrVbL6tWrKz3H3d2d06dP07x5c3uGKmpJ2tb1WEyUc+fOZdSoUTzzzDMyP2ln7VNSaJ+SwpG77wEkWbqaDh06sGLFCvP369atIzw8HIPBwIYNG0hLSyMqKoqoqChSro40nDp1ijlz5nDhgrpbDgnLpG1dj8Wh13fffZcvvvgCRVHQaDR88MEH9opLXGXqXQJoU1LM11vKMKxryc7OJjg4GDD2LsLCwggLC6twjI+PD/PmzaNRo0Y3nZ+enk56ejoAWVlZNo5W1Ia0rfOzmCiPHDkiydEBXL/Qx7RtV0JCAjqdTuXIhLUEBgaSnZ1N7969q5zmiIyMrPL8kJAQQkJCAOP+i8JxSNs6P4uJsrS0lB9++MH8Kadbt252CUrczJQsTQt9dDqdFG52Ynl5ecTHx5ORkUFSUhITJ04kNjaWzZs3Exoaekv3nZ2dXekbalZWFkFBQbd039bmTDEFBQURFRVV7fnStkbOFFO1batYoNfrlYSEBPP/atLpdKo+viM5PnKkcrjP3YqiKMqsWbOUWbNmqRyRcAaO+BqSmKzDEWN2pZgsLuYJCwujTZs26HQ67rvvvlvM5cJa2qek0KBzJ05ERBAfH4+fn58s8hFCCBuxmCgXLVrEiRMnAPj000/tEpComeuvuRy0bRuAFCcQFpnmuRyJxGQdjhizK8VkMVE2bdqUBg0aAFBUVFSnBxC20z4lhbt+3MPlo8cYtG0bfn5+JCQkqB2WcFCu9MZlS44YU3UcMWZXisliouzfvz/btm3j8ccf589//nOdHkDY3l0/7gHggXeWotPpZGWcEEJYkUZRFEXtIGpCr9dLAqjGiYgIijMP4B3cg1XdutGmTRu53lIIIW6RxctDnnnmGTQaDQUFBfzvf/9j7969dgpL1IVp3vJwt+5EHD1GyhPDpDiBEELcIouJ8sMPPzR/vWDBAlvHIqyk66GDnIiIYOSHa1g35lnpjYsKtm/fzs6dO2nTpg1jx45VOxwANm/ezO7du+nYsSPPPvus2uGYJSUl0b59e55++mm1Q6kRaduaq03bWkyUS5YsAYyFB3bv3m2d6IRdmHqXT3brTpm7OwkajVTyqSeOHz9OYmIiFy5cYO3atRQWFvLiiy/i5eVFSEgI2dnZxMXFMWfOHIeJafTo0YSEhLBo0SKHiSkwMJDg4GAKCgrsFlN1pG2tE1Nt29biYp7u3bvTo0cP+vfvz/Lly63yCwj76nroIB5XrvDUh2ukV1lPVFeUW6PROFxM5eXlzJs3j5iYGIeJac+ePezatYvvvvvObjFVR9rWOjHVtm0t9ig3bdpk/sOnpqYCxh1FhHPpeuggABF3dUWP1Iusb24syt2/f39mz55NQECAw8Q0d+5cLly4wPfff89jjz3mEDFNmTKFrKwsh0qUN5K2rVtMtW1bi4nyypUrjBo1Co1Gw/vvv8+kSZNuOWChHu9ePRm5JkWSZT1zY1HugQMHMnDgQIeKKS4uTtV4KosJjDVAHa1e6fWkbesWE9SubS0OvRYVFdG3b1/uueceLl26RLt27W45YKGe9ikp+PTqScSaFL4ZMEDK3rmovLw8YmJizEW5n3zySVJTUxk/fvwtF+WWmNTliDHXh5gsXkf50Ucf8dFHH+Hm5sZTTz3FU089dUvB3wpZuWk9JyIiKN63n1IPd9aFh8vfVQghLLA49DpgwAB8fHwIDQ2VayhdiPl6y7u6ylCsEEJUw+LQ6/Tp0807a//73/+2RzzCjroeOUy3I4cZuSaF/T16SFF1IYSohMVEedttt9GsWTN7xSJU4tOrJx5lVwh7733pWQohxA0sJsouXbrw+eefM3LkSO644w6bBJCUlMSaNWtsct+iZtqnpNDtyGHcgJFXF/rILiRCCGFU5Ryloij4+/uzdetWFEWhcePGFu+oJtUZbvTVV185XOWL+qzrkcMc7tYd37zfGf7xx+gVhYCAAKkVK4So16rsUWo0GjIzM7l48SJubm7V7kdZXSWEtLQ0oqKiiIqKIuXqYhJHrHxR33U9dJCuRw7jWXaFQVu3kpOTo3ZIQgihKourXk+cOMHUqVMBY+KszYKeGyshhIWFERYWVuGY6qojpKenmxcTZWVl1fixxa3z7tUT3337eerqqliN1IoVQtRTVSbKefPm8dhjj9W5an5llRAqY6k6QkhIiHlHallkYl/tU1LM11tGrEkh5ekIEhISJFkKIeqdKodev/76a2bNmsXIkSPN/yxxxOoM4ta0T0mh65HDAESsSUFRFPnAIoSodyxW5tm8eTNDhgyxZzxVkso86joREcGlffsB+OjpCGkLIUS9YfHyEEdJkkJ97VNScGvYEOBqgYJglSMSQgj7sJgohbjeXT/uwadXTwA8ysrY36OHyhEJIYTtWUyUW7du5Y033gCQogACuFac4HyLFniUXeHQXV3VDkkIIWzKYqJMTU2lrKwMgN27d9slIOEcBnzzNedbtACMxdWFEMJVWUyUPj4+AJSXl5OXl2eXgITzGPDN1/zuezsKkiyFEK6r2sU8u3fv5oknnmDUqFH2ikk4ka/++tcKyfJERITaITmN/v37y+phIZyAxUSZnZ3NBx98wMaNG3nkkUfsFZNwIvHx8fzn0Ucp83BHAYozD6gdkkP75Zdf6Ny5M/feey9nzpxRO5xbkpycTO/evQEoLS2lXbt2bNq0iR9//JF27dpx6dIlAKKiohg3bhwjRoxgy5YtgPF5Exsby/jx4zly5Ihav4JVHTx4kH79+tGuXTuGDh1abdlPR3Pw4EFGjRrFxIkTmTt3bp3u48SJE7z88su3HEtaWhoxMTEMGTLEIUqcWixh16pVK6ZMmYKiKISFhcnlIqJSOp0O/dXLcSPWGCv6mDaHFhW98847+Pv7k56eTqdOndQO55bdeeed7Ny5k99++4377rsPgFWrVjFz5kw+/vhjxo4dC8D8+fPJz89n7ty5PProo+zevZvPPvsMN7eKn9XPnj3LhAkTCAoKYujQobRr14758+ejKAodO3YkJiaG6OhogoKC+Oyzz9i+fXu1GzbYy+uvv87Ro0e57bbb2LFjB1u3buWJJ55QO6wa27JlC2PGjGHw4MEA9O7dm8GDB3P69Gl69erFd999R0xMDHfeeSeLFy8mLy+Pxx9/nLCwMLp3786YMWMoKSlh586dvP3228TGxgJQUlJCZGQkgYGBDBgwgIceegidToeiKDRp0oRZs2bxf//3f7Rs2ZKdO3eyaNEic8nTjIwMtm3bRv/+/dX801juUT7++OPMmzeP4OBgpkyZYq+YhBMKCAgAMPYq9+3ncLfu6gbkwDQaDRqNRu0wrCI8PJzU1FS2bNnCo48+CsC5c+cYPXo0GzduNB/3yiuv8MgjjzBu3DgAJk+eTExMDFqttkLh/QsXLphrQz/wwAMsWbIEHx8ffH19yczMZOvWrTz66KMkJibSunVr+/6y1VAUhcGDB5OVlcX27dvN5TedRXR0NN988w3R0dH861//ws/Pj6SkJJo3b87DDz/MW2+9RVpaGh4eHly+fBl/f39Wr14NQJs2bYiLi2Ps2LHcf//95iQJcPnyZQoLCxk8eDBDhw7lgw8+4NKlSzRv3pzjx4+zb98+AgICSExMpFu3bubz5s+fz+TJk3n88cft/re4kcVEGRUVxauvvkq3bt04fPiwvWISTsi0FddHT0egAJSXy3xlJWJiYjh9+jT9+vWzWAPZWZgW/Pn7++Pm5kZeXh5nzpwhNjaWn3/+mZ9//hkw1o5evXq1eeegxx9/nGXLlhEdHc3y5ctZsGABkyZNIiAggIULF5KZmcnMmTMpLy9n9OjR6PX6CrsTAQ73YSM+Pp5Dhw7h7+9PREQEJSUlaodUK02bNiUxMZEVK1awYcMGc0+9QYMGNG3aFC8vLy5fvsx7773HsGHDmDp1Kvn5+QA0a9YMoMIIwapVq5g0aRIXL17kvffeIy8vj9jYWMrLyxkyZAh6vZ4PPvigwjnXt+lLL73ERx99xJtvvmmPX98ii0Ov77zzjvmFIER1TGUGyzzc8Sy7QnHmARmGvUGnTp04duyY2mFY1dy5c9FoNKxatYrJkyeTk5ODj48P+/fvr5DcevXqxfz58/ntt9+YOnUqjRo14ty5c7zyyiv06dMHgMzMTFasWEFxcTGDBg2iX79+TJ06ldatW9OkSRNeffVVXnjhBX766SdOnTql1q9cqZ49e7Jv3z61w6iztLQ0vvjiCzw8POjWrRsnTpyo9Lj777+fpUuXsmPHDry8vCr8rFWrVhw7dow333yTKVOmMHbsWE6fPs20adNwd3ene/fuPPvss8TGxvL1119TUlLCP//5T95++22mTp3KgQPGNQ7Lly9n3759XLhwwSH2w62y1uvf//5389ZWiqKg0Wj46KOP7BlbBVLr1XkkJCTw1Idr0IC5qLoQ1hYVFcXbb7/tMHOU4tbp9XrCw8Pp4WBVvywWRd+xYwcDBgwAYPv27QwcONBugd1IEqXzMLVTxBpjT1KSpRDCmVU5R3no0CE+/PBDDh06xKFDh1i5cqU94xJO7PqFPQBH7r5HvWCEEOIWVZkoP/74Y3766SfWrl3LunXriJCFGaKGrl/YA6A42fVkQghxvSoX8+h0OsrKysjIyKCoqMjhVpgJxxYQEEBOTg55vr745uVx+K6uMgQrhHBKFi8PiY2NZebMmXz11VfMmzfPXjEJF2DqVW57ZBDeV7fmknqwdbd9+3Zmz57NqlWrqj02OTnZ9gEJq5G2dXwWE2Xjxo0JDg5mxowZqi7kEc7JNFdZcPCQ+TaZr6zc8ePHiY6OJjw8HIDCwkIiIyPRarWsXr2ab7/9lri4OE6fPl3tfZlWqwvHIG3r/Cwmyj59+tC9e3eeeOIJjh8/bq+YhIsw9SoVFDQNGwLQoLPzl22zhQ4dOlS45nDdunWEh4djMBjYsGGDTH04MWlb52ex4MDo0aMr/C9EbQUEBLAuPJxBW7fiW1RE8b79aofkFLKzswkODgbA3d2d/v37M3v2bHMv/Ubp6emkp6cD0utwdNK2zqfKRPnUU0/d9ElHzYIDwjlptVr0ej3bHnnEfF2lLOypXmBgINnZ2fTu3Zvy8nIGDhxocfojJCTEXFtUrjd2bLZo2x92j6Bf31QbRCvAQqL8+OOPKSgo4PDhw/Tr14/s7Gx7xiVckGkFLCCl7W6Ql5dHfHw8GRkZJCUlMXHiRGJjY9m8eTOhoaFqh+eUCkoKOPbHMTrd1onGXupV77FH2zprkiwvLqP0tyI8/Rvi5m1xgFNVFiObPHkyXl5e9OvXj6SkJBYvXmyvuIQLMVVV2vbIICIPH6Z4334Zgr2Br68vS5curXCbFPmou03HNzHz25lcKruEj4cPM+6bwdAOQy2eM3fuXAICAqw+1WSPtjUYDA5RE7U2ijJyOb/+KEpJORovN5oP70zDPn5VHm+r9qkJi4t5mjRpgr+/PwDe3t52CUi4Lk9PzwoJUnYXsa9BKwepHYJdFJQUmJMkwKWyS7z+7esUlBRUeU5hYSEtWrTg7Nmz9grTqp4o6ad2CLVSXlxmTpIASkk559cfo7y4rNLj1W4fi4myRYsW7Ny5k5dffhl3d3d7xSRckEajobS0lG9euravqfQq7Wvbc9vUDsEujv1xzJwkTYrKivjlwi9VnrN+/XpycnI4dOgQR48etXWIVuc3obfaIdRK6W9F5iRpopRcoTS38ipearePxUQ5depUFixYwPPPP8+DDz5or5iEC9LpdACcOnVKFvIIm+p0Wyd8PCpuD9jQoyEdm3Ws8pwjR44wffp0Fi5caN4z09lt/vQBtUOokqd/QzReFdOPxssdT7+GlR6vdvtUmSjnzp3LqFGj+PHHH5k9e7YsSxZWoSgKiYmJ5u9l+FVdrjgc29irMTPum2FOlg09GjL9vukWF/TMmjULME4xTZs2zS5x2lq/E4sqfH9At0WlSG7m5u1B8+GdzclS4+VO8+GdqlzQo3b7VLmY5/Dhw3zwwQf8+c9/Ztu2bTJHKW6Zqf5raWmpeWeRy0ddaxNjZ+Oqw7FDOwwlJDCEXy78QsdmHVVd9aqWG4djeyQ8ChgTpulrNTXs44d319spzS3C08+xV71W2aM8d+4cn376KS1atOC///0vn376qT3jEi7o+h3pNW5uaJCdRRyNK/UwG3s1plfLXvUySVpyfcJUm5u3Bw3uaOrQSRIsJMrw8HDOnj3L8OHDOXv2rNOuBhOOQ6fTmYtYlLpdK2Zx5O57ZAjWQZh6mK6QMIuLizl58iTFxcU1Ol6v1/Pmm2/y/PPPU1BQ9QpZV9Ej4VEO6LaQu3ivajHUpo3UbJ8q03hkZKQ94xD1QGJiIopiHHTd+MwzhHz2ObefO4dSVCTFBxzM9QnTGYdn9+/fz8aNGyktLcXT05PQ0FB69uxZ7XmKolBYWEjDhpUvKnFmlVXv6ZHwKAaDgScW23/lbF3aSK32sbjqVQhrio+PR6PRmHuVhyLHmn8mPUrHtO25bU7XuywuLja/AQOUlpayadOmGvVa/va3vzFy5Ej27t1r4yjtr6rqPVqtlk+8fiB38V679S7r2kZqtY8kSmFXOp0ORVEoLS2la/K75kU9ck2l43K2ZHn27FnzG7BJSUlJjaaPFi5cyH//+1+CgoJsFJ1jMiVLvwm97ZIs69pGarWPY8+gCpdkWv16OCqSAfPfNN7o5ibF0h3YjXOXjjwc27JlSzw9PSu8EXt5edGyZUuL59X3YvJardZYCm+ClpwZOwiYOcBmj1WXNlKzfaRHKewuNzcXgJycHM63aIECaLy98e5V/RySUJczLPbx9vYmNDQUT09PwPgGPHToULnErQZMyTJg5gCb9iydrY1U7VEePHiQDRs2cPfdd/PYY4+pGYqwIz8/P3JycgDMi3uUoiKKMw/IriJOwJQs+6zsgy++Dtm77NmzJ126dOHs2bO0bNnSYd+A1VDdllz26lk6UxtZrUd5/PhxoqOjCQ8PB4xFbCMjI9FqtaxevbrSc9avX0+jRo2sFYJwEqZdDjQaDdseGWS8nhKg3Fj7URb2OIeM5zLY9tw2+qzsQ8iyEAwGg9ohVeDt7U3btm0d+g1YDbXZksvWPUsPj1KaNj2Lh0dp9QeryGqJskOHDqxYscL8/bp16wgPD8dgMLBhwwbS0tKIiooiKirKXKsvNzeX6Ohovv32W2uFIZxAYmIier0eRVF4cu1aUp6+lhhlUY/zyXgug4CAABZ6LqT3yt5qhyOswNSrBCj9rdAmj3HmzCd8s+N+du95im923M+ZM59YPH79+vVs27YNg8HAjz/+aJOYqmKzodfs7GyCg4MBcHd3JywsjLCwsArHRERE8NZbb+Hr61vpfaSnp5Oeng4gtWZdSHx8PHB1cl6vx7Nbd3BzM/coLx89JkOwTmb1EOOo0ejNowlONr7un8552tzWwvmYh2BnWn8ItqwsnyM/TePKFWNlritXijjy03RatPgLHh5NKj1n+PDhvPzyy/j5+fHmm2/y17/+lTNnznD33XezceNGIiIiSElJoVGjRowePZrVq1fTunVrhg8fzuzZs7nzzjuJjIys0zWYNkuUgYGBZGdn07t3b8rLyys9ZsCAAQwYUPUfPyQkhJCQEEBWpLmaxMTEayveRj7FyDXGpHjF0xOP68raScJ0LqaECRCcHMya5DX4K/4OOY8pasaULHMX77VaUYLCwqPmJGly5UohhYXHaNasT6XnlJaWcvnyZQoLC7ntttsYNWoUL7/8snmNy6+//mreQ/nEiRP07NmTY8eOUV5eTs+ePcnNzTWviagtqw295uXlERMTQ0ZGBklJSTz55JOkpqYyfvx4QkNDrfUwwkX4+VXcydz09PW4mjxNxdIlSTqvzKhMMqMy+U3zG8HJwQ43h1nf/bB7RLXHmNYTgHWHYBs16oy7e8Wenbt7Ixo16lTlOQsXLmTKlCn06dOHTz75hKVLl9KiRQsA3NzcGDhwIEVFRZSXl9OpUyf++OMPLl++zMmTJ8nPz+f8+fPk5eXVLWDFSeh0OrVDEMIpOOprpcfKHkqPlT3UDsOpqdG2y5YtM3+dPf0bq93v6dNpypfpPZRt/+mgfJkerJw+nVbjc+Pj460WR01IwQGhmuuHX59cuxaPsivmn2nc3PAO7iE9SheSGZVJn5V9zHOYmVGZKkckasI8V6nV4ulvvasUWrV6ghYt/kJh4TEaNepU5dxkZUz7U9qLJEqhmgoLPfR6Dt/V9dr35eUU79vP4W7d6XrooP2DEzaR8VwGYJy/DE4OxhNPfoyy7wpGUTe2uLbSw6NJlXOSjkQq8wiH0fXI4Zv/SZJ0SZlRmUwsnUgppQQnBzN682i1Q6qXajJPCRXnKusj6VEKVRgMBk6dOoWiKAzauhWA2/N+v+k4n6tl7WQI1vVotVq0aOm9sjf7z+4nODlYhmPtrDbFB+DaCtjs174mMOkBG0XleCRRClVU+IQql/7Ua3uf20tCQgJr260lODmYni17VrjMRDiG64sQeAU2serlIo5Ohl6FEKrT6XTm3uT+s/tlKNaBGQwG/Cb0piQ7X+1Q7EYSpXAYR+6+h8N3da3035G77+Fwt+5qh6i6pKQk1qxZo3YYNnN9srx39b0qR2NfarVtTecpoeJIUGDSA+TM2GGLkByODL0KVZiWmxsMBvNOIgyrWJjCtAUPGFfIGgwGuuIajh8/TmJiIhcuXGDt2rUUFhby4osv4uXlRUhICKNH39yj+uqrrwgODqagoECFiO1nYulEFnsspqi0iHtX38v3o79XO6Racba2re08pWmbPAClpPKqa65GEqVQhemTaW1W07nSyjvTJgKm3XZMmwiEhoYSERFBo0aNSEtLA2Dw4MFERESwZ88e/vjjDy5evMjTTz+tYvS2pdVqwQALPRdSVFrE6M2jnWrO0tXbNj4+nsTEROLj4wmc/UC9WNgjiVIIB1CTTQSmTJlCVlYW33333U3nu9oGAlqtlhx9DqlBqew/59w7yrh622o8XX8GTxKlcBimAugVCg9UouuRw3aKyH5qsokAQFBQEEFBQTfd7oobCBh3l4HUoFSC3w0mM9I5Lx1xhratbjPnG13fqwyYOYDsuK8JnO26vUpJlMJhmK6VdMVEeKO8vDzi4+PNmwhMnDiR2NhYNm/eLJsIXOf6ZOksQ7DO2La1nae8icY6cTgqSZRCqMDX15elS5dWuG3lypUqReMcnGUI1lnb9lZ6lYFJD7h0r1ISpVCNwWAgNzf32r6U1dBojB9bdTpdhfuAioWbTbffeJtwTnq9ni//9SW/N/jdKVfBOou69CrLysoqfO+qRQgkUQrVWCOBXX8flX0tSdI1PHT6IVKDUikqK6r+YFFnte1VtmnTpuIK2LivbRidelx/uZIQwunp9Xpuv3w7KBD8brDa4bis2vYqtVpthV6lxsuN7NdcL1lKohQO40REhMXqPFX9AzjcrTsnIiIq3Fd1jyWcy0OnHzJ+oVDvqvbY0w+7R9SqWo9OpzOvxg2YOQAUGwWmIhl6FQ7jVnYIuXE7ruruS3YjcT6enp6MyBohQ7A2Vpe5So1GY14PoPFyc7mFPdKjFEI4hQobfSvQ691e6gXj4urSqzSVojRt6py7eK8tQlOFJEohhNMICAgwf11OuQzB2ki/vqn065tK+lc1nw/WaDQkJiYC4NW2CSUnXWd3EUmUQgincvvl241fKMgQrI2FPJhZ456lTqczX+rlN6E3aHCZVbCSKIUQTkOr1V5b1AMuuXDE0RQW/lzjZBkQEGBe2OMV2MSGUdmXJEohhNMx9yoxXi4i85W2E/JgZo0X+Fxf8MOVepWSKIUQTiUgIICHTj/EiKwRjMgaQWZkJvsi96kdlrhKr9ebF/aYtt9y9g2eJVEKIZyKVqvlaL+jpLZPJbW9cWeR4HeDGb355g2RhXWYhl5rsxJWr9djMBjwatsEpaTcqXuWkiiFEE6ny+4ujDgxwvxvYslEp9hZxFmZVsD265tao2Rpmqc8depUhdqvztqzlIIDQgin83Pfn2/aTWThuwsBaOjRkOKyYnq06CHJ04oaNepS6x6lXq8nISEB3Wwd2XFfo5RUvR+nI5MepRDC6XT+oXOFHqXpX2ZkJi8UvcC+yH2EnApRO0yXYlrQU9NeJRirKSmKgsFgQONlTDfOOAQrPUohhNMxDe1VRnaOsZ0bk2V1q2Hj4+PNi3tyy9rQgiZo0JAd9zVebZs4zZZc0qMUQghRY+lfBfOf/3YGara4x/Sh5hOvH7jQoNh8e8nJfLJf+5qcGTscvtyd9CiFEELUWMiDmbU+R6PRoCgKazU7eZ6/4O7lbpyvvFowwtF7ltKjFEIIYVM6nQ5PT0+GXe7Lv73/S3nJFZTryio5+mpY6VEKIYSoNdOwa2Hhz1y5UkzTpj25eHEvD//lF/PPr5/D7N3nczbsKmXY5b4AaNAYk+XVlbCmDZ81nm7mHUgchSRKIYQQtVbdQp4bfz7kr19zKscAx6GMK3jgXvGEqx1MpaScnBk78PRvhN+E3uQu3qv60KwkSiGEEHazp8Npc4k7sLyCGRxj/lLVOcr169ezYMECZs6cqWYYQggh7ECr1VZIknCt1J0js1qP8vjx4yQmJnLhwgXWrl1LYWEhL774Il5eXoSEhDB69M11GN3d3Tl9+jSBgYHWCkMIIYQD02g06HQ6TsZtv3bjL5i/16C56RyvtsYtu9TqXVotUXbo0IEVK1YQHh4OwLp16wgPDyc0NJSIiAgaNWpEWloaAIMHDyYiIoJTp04xZ84cZs2aVel9pqenk56eDkBWVpa1QhVCCKESnU4HQNvZA1WOpOZsNkeZnZ1NcHAwYOw5hoWFERYWVuEYHx8f5s2bR6NGjSq9j5CQEEJCQoDqx7GFcHXr16/n119/5eLFi8yYMUPtcIQVSds6NpslysDAQLKzs+nduzfl5ZUXwo2MjLTVwwvh0GSqwnVJ29aMpc5PQEAAp06dok2bNmi1WgwGg6olCa2WKPPy8oiPjycjI4OkpCQmTpxIbGwsmzdvJjQ01FoPI4RLsMVUhXAM0rY1U5tRQrXr9lotUfr6+rJ06dIKt61cudJady+ES7vVqQqZz3dc0rbOT66jFMIB3OpUhcznOy5pW+cntV6FUEFeXh4xMTHmqYonn3yS1NRUxo8fL1MVTk7a1vVoFEVRqj9MfS+88EKNJ7qzsrIICgqybUC1JDHVnCPGVZuYgoKCiIqKsmk8llT1WnH2v6u9VBWT2u0K0ra3qs5tq7ggnU6ndgg3kZhqzhHjcsSYassRfweJyTocMWZXikmGXoUQQggLXDJRmia+HYnEVHOOGJcjxlRbjvg7SEzW4Ygxu1JMTjNHKYQQQqjBJXuUQgghhLVIohRCCCEskEQphBBCWODSlXm2b9/Ozp07adOmDWPHjlU7HMAY05IlS1izZo3aoVSwefNmdu/eTceOHXn22WfVDgeAH374gS+//BI/Pz/Vr1+7XlJSEu3bt+fpp59WO5Q6ccTXhSM+/8D52lratuZq07ZOnSirq9KfnZ1NXFwcc+bMcZiYRo8ezc6dO+0WT23iCgkJYdGiRQ4VU2pqKn5+fg4TU2BgIMHBwRQUFNgtptpy1teFoz3/HLGtpW2tE1Nt29aph15NVfpNTFX6DQYDGzZsQKO5eadstWNSS3VxlZeXM2/ePGJiYhwmJoDZs2fb9Y2qupj27NnDrl27+O677+wWU2054+vCEZ9/jtjW0rbWiam2bevUifJG2dnZtG3bFjBW6e/fvz+zZ8+mTZs2DhNTZmYmO3bsYMuWLarFVFlcc+fO5cKFC3z//fcOE9OGDRvQ6/V4eno6TExTpkzh+eefp3///qrFVFvO8LpwxOefM7S1tG3dYqpt2zr10OuNbqzSP3DgQAYOHOhQMQUHB7Nx40ZVY6osrri4OLVDuimmYcOGMWzYMIeKCYx1IR2thqUlzvC6cMTnHzh+W0vb1i0mqGXbWrOOnr2dO3dOGTdunNKhQwfljTfeUAoKCpSoqCglJiZGef/99yUmB49LYrINR/wdJCbrcMSY60NMUplHCCGEsMCl5iiFEEIIa5NEKYQQQlggiVIIIYSwQBKlEEIIYYEkSiGEEMICSZRCCCGEBZIohRBCCAskUdZScnIyQ4cOJTIysk6Fh8PDw60Sh6m6BIBer+fAgQMAfP755yQnJ5OcnEzv3r0BKC0tpV27dmzatInk5GQ2bdpU4b6Sk5PZsWOHVeKyltr8nU6fPs306dNtGE395IrPdWEkbVs7LlXCzl5iYmIYOnSoeXuWxMREzp07R35+Pm+99RaffPIJGRkZ5Ofns3jxYr799luWLVtG586dATh8+DAJCQnccccdjB07lh49egCQnp7OP/7xDx5++GFyc3NJSkrinXfe4eeff+aPP/7g9ddfZ9q0aQQFBREcHMyIESMsxnnnnXeyc+dOfvvtN+67774qj9u2bRtRUVH8+uuvTJs2DT8/P4YPH46HhwdLliwB4MUXX6S4uLhCfBMnTmTOnDksWLCAN998k/vuu8/8OHq9nvz8fDw8PLjrrrt47rnnCA8PZ+3atSxdupS77rqL9PR0zp8/z/nz5+nZsyd5eXlkZ2fz3nvvcfnyZXQ6HT/99BNTp06lWbNmzJ8/H0VR6NixI2FhYYwePZrQ0FCio6PJycmhpKQELy+vW25fcY2rPdfFNdK2NSc9yjowGAzcf//9DB06lMOHD7N9+3aaNWuGp6cnhw8fBsDT05OcnBwyMjJ45513WL58ubl6fm5uLs2bN2fUqFHmJ5dJ//79mTx5MkVFRZw8eZJVq1bRrFkzbrvtNn788UcAtFpttU8uMH7qS01NZcuWLTz66KOVHpOfn2/+evHixcyYMYP58+fz5z//mUWLFrFs2TKWLVvG4sWLb4oPoKSkhPPnz7Nr166bnsQjR45kzpw5FgvAR0REsGjRInbu3ElSUhI+Pj6cP3+e0tJSpk6dyqJFi1i8eDFLlizBx8cHX19fMjMzAejWrRtxcXG0bNmSDh06kJGRUe3fRNSOKz3XRUXStjUnPco60Gq1/OUvf2HcuHH06dOH7t27o9frzT+fNWsWGzZsICEhgaKiItzc3PDw8KBBgwYAPPjgg+ZtYPbt20dkZKT53LKyMsA4zAAQEBBQ4b7XrVtHs2bNKsTTvHlzzp49CxifvLfffju///47Pj4+APj7++PmVvlnovz8fJo2bQqAoigVjlMUxbxtj6nS4Y3xabVannnmGUaNGnXTfTdq1KjCuab7LiwsNB/TtGlTGjRoYI7By8uLy5cvA6DRaMyPX15ezpgxY+jZsycAWVlZFf4OzZs35+LFi5X+jqLuXOm5LiqStq05SZR11LBhQ/70pz9x/Phx3NzcmDJlCpcuXWLq1Km0bt2auXPnsmvXLh588EHGjRvH3//+d1q1agXAl19+ycaNG/n999955JFHKtzvrl27ePXVV/Hx8aFt27b86U9/4v/9v/+Hoig8//zzlcbyzDPPMHnyZD755BMuXLjA0qVL+fDDDwGYO3cuGo2GVatWVXquv78/eXl5gHF4Va/X07p1a4YNG0ZsbKz50+OLL75ISUlJhfhat25N69atKSoqIiIiotq/WUBAAP/4xz/45ptvuOeeeywe6+npyaxZszh69CivvfYazZo1M/9tmzRpUuFFCfDLL78wdOjQamMQtecqz3VxM2nbmpGi6A4kPT2dAwcOEBsba9fHffXVV83zgJZUFl9cXBxt27ZlwoQJtg6zSoqiEBkZKW+QTkSt57qwPVdsW0mUgpMnT5Kfn0+3bt3UDqVOcnNzOXnyZLW9VCGEqAtJlEIIIYQFMusthBBCWCCJUgghhLBAEqUQQghhgSRKIYQQwgJJlEIIIYQF/x+EBcxdcpfcOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 450x225 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ncols = 2\n",
    "figheight_in = PAGEWIDTH_IN/ncols\n",
    "\n",
    "path='figures/main_observed_amplification.pdf'\n",
    "with mpl.rc_context(PAPER_CONTEXT):\n",
    "    \n",
    "    fig = plt.figure(figsize=(PAGEWIDTH_IN,figheight_in), constrained_layout=True)\n",
    "    spec = fig.add_gridspec(2, 4)\n",
    "\n",
    "    ax_big = fig.add_subplot(spec[:2, :2])\n",
    "    ax00 = fig.add_subplot(spec[0, 2])\n",
    "    ax01 = fig.add_subplot(spec[0, 3])\n",
    "    ax10 = fig.add_subplot(spec[1, 2])\n",
    "    ax11 = fig.add_subplot(spec[1, 3])\n",
    "\n",
    "    axmap = {'CELseq2':ax00, 'DropSeq':ax01, 'MARSseq':ax10, 'SCRBseq':ax_big, 'Smartseq3':ax11}\n",
    "\n",
    "    for df,exp,protocol,run in zip(dfs,exps,protocols,runs):  \n",
    "        \n",
    "        if protocol=='Smartseq3xpress':\n",
    "            protocol = 'Smartseq3'\n",
    "            run = 'Xpress'\n",
    "            color = 'tab:olive'\n",
    "            zorder=4\n",
    "        elif run=='SE':\n",
    "            color = 'tab:pink'\n",
    "            zorder=4\n",
    "        elif run=='A':\n",
    "            color = colormap[protocol]\n",
    "            zorder=4\n",
    "        else:\n",
    "            color = 'tab:gray'\n",
    "            zorder=3\n",
    "            \n",
    "        ax = axmap[protocol]\n",
    "        if run == 'SE':\n",
    "            runlabel = 'A'\n",
    "        elif run == 'PE':\n",
    "            runlabel = 'B'\n",
    "        else:\n",
    "            runlabel = run\n",
    "        \n",
    "        reads_per_umi,abs_frequency = np.unique(df['N'],return_counts=True)\n",
    "        rel_frequency = abs_frequency/sum(abs_frequency)\n",
    "        ax.scatter(reads_per_umi,rel_frequency,\n",
    "                   label=runlabel,c=color,\n",
    "                   linewidth=0,s=POINTSIZE,rasterized=True,zorder=zorder)\n",
    "\n",
    "        ax.set_yscale('log')\n",
    "        ax.set_xscale('log')\n",
    "        if protocol == 'SCRBseq':\n",
    "            ax.set_ylabel('Relative frequency')\n",
    "            ax.set_xlabel('Reads per UMI (copy number)')\n",
    "            plt.sca(ax)\n",
    "            plt.xticks([10**0, 10**1, 10**2, 10**3, 10**4])\n",
    "            plt.yticks([10**-8,10**-6,10**-4,10**-2,10**0]) \n",
    "        elif protocol in ['MARSseq','Smartseq3']:\n",
    "            ax.set_xlabel('Reads per UMI')\n",
    "            plt.sca(ax)\n",
    "            plt.xticks([10**0,10**2,10**4])\n",
    "            plt.yticks([10**-8,10**-4,10**0])\n",
    "        else:\n",
    "            plt.sca(ax)\n",
    "            plt.xticks([10**0,10**2,10**4])\n",
    "            plt.yticks([10**-8,10**-4,10**0])\n",
    "\n",
    "        ax.minorticks_off()\n",
    "        \n",
    "\n",
    "        add_largedot_legend(ax,loc='upper right',kwargs=dict(fontsize=LEGEND_FONTSIZE,frameon=False))\n",
    "        ax.set_title(prot2print[protocol])\n",
    "        \n",
    "        ax.set_ylim(10**-8.3,10**0.3)\n",
    "        ax.set_xlim(10**-0.3,10**4.3)\n",
    "        \n",
    "        for axis in ['bottom','left']:\n",
    "            ax.spines[axis].set_linewidth(SPINEWIDTH)\n",
    "    \n",
    "    sns.despine()   \n",
    "        \n",
    "\n",
    "   \n",
    "    letters=['a','b','c','d','e']\n",
    "    x_corrections=[0.1,0,0,0,0]\n",
    "    y_corrections=[0.085,0.14,0.14,0.14,0.14]\n",
    "        \n",
    "    for ax,letter,x_corr,y_corr in zip([ax_big,ax00,ax01,ax10,ax11],letters,x_corrections,y_corrections):\n",
    "        ax.text(LETTER_LOC_X+x_corr,LETTER_LOC_Y+y_corr,letter,transform=ax.transAxes,fontweight='bold')\n",
    "    \n",
    "    plt.savefig(path, dpi=300, format=None, bbox_inches = 'tight', pad_inches = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_observed_stats(dfs,quick=True):\n",
    "    if quick:\n",
    "        print('           \\tmax\\tmean\\tff\\tnum\\talpha',end='\\n')\n",
    "    else:\n",
    "        print('           \\tn\\tmax\\tmean\\tff\\tnum\\talpha',end='\\n')\n",
    "    stat_rows = []\n",
    "    for df in dfs:\n",
    "        \n",
    "        row = dict()\n",
    "        \n",
    "        exp_name = df.index.values[0]\n",
    "        if 'Smartseq3' in exp_name:\n",
    "            exp_name = exp_name.replace('Smartseq3_Fibroblast','Smrtsq3')\n",
    "            exp_name = exp_name.replace('Smartseq3xpress','Smrtsq3xpr')\n",
    "    \n",
    "        print(exp_name,end=':\\t')\n",
    "        if not quick:\n",
    "            n_cells = len(np.unique(df['RG']))\n",
    "            print(n_cells,end='\\t',sep='')\n",
    "\n",
    "        data = df['N'].values\n",
    "        max_=np.max(data)\n",
    "        print(int(max_),end='\\t',sep='')\n",
    "        mean=np.mean(data)\n",
    "        var=np.var(data)\n",
    "        print(np.round(mean,decimals=1),end='\\t',sep='')\n",
    "        print(np.round(var/mean,decimals=1),end='\\t',sep='')\n",
    "        print(len(data),end='\\t',sep='')\n",
    "        print(np.round(mean + (var/mean),decimals=1),end='\\t',sep='')\n",
    "        print('',end='\\n')\n",
    "        \n",
    "        row['exp'] = exp_name\n",
    "        row['mean'] = mean\n",
    "        row['ff']=var/mean\n",
    "        row['max']=max_\n",
    "        row['alpha']=row['mean']+row['ff']\n",
    "        row['n_UMIs']=len(data)\n",
    "        row['n_cells']=n_cells\n",
    "        stat_rows.append(row)\n",
    "    return pd.DataFrame(stat_rows)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "           \tn\tmax\tmean\tff\tnum\talpha\n",
      "CELseq2A:\t34\t3476\t27.2\t107.5\t2140365\t134.8\t\n",
      "CELseq2B:\t37\t11092\t24.4\t199.8\t4303956\t224.2\t\n",
      "DropSeqA:\t42\t2463\t29.6\t284.2\t2506244\t313.8\t\n",
      "DropSeqB:\t34\t3718\t31.1\t414.1\t1272895\t445.2\t\n",
      "MARSseqA:\t29\t1624\t24.1\t132.1\t1342232\t156.2\t\n",
      "MARSseqB:\t36\t1719\t20.9\t107.0\t1903673\t128.0\t\n",
      "SCRBseqA:\t39\t4840\t9.7\t218.8\t9429371\t228.5\t\n",
      "SCRBseqB:\t45\t3276\t9.4\t159.3\t6800371\t168.7\t\n",
      "Smrtsq3:\t145\t1216\t5.4\t8.0\t21549849\t13.4\t\n",
      "Smrtsq3_PE:\t319\t989\t3.8\t4.5\t48073893\t8.3\t\n",
      "Smrtsq3xprA:\t96\t173\t1.3\t0.3\t6940889\t1.6\t\n",
      "CPU times: user 3min 48s, sys: 21.6 s, total: 4min 9s\n",
      "Wall time: 4min 9s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "observed_amp_stats_df = get_observed_stats(dfs,quick=False)\n",
    "observed_amp_stats_df.set_index('exp',inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "856"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observed_amp_stats_df.n_cells.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "106263738"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "observed_amp_stats_df.n_UMIs.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Prepare Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_fig4_name = 'DropSeqA'\n",
    "protocol_fig4='DropSeq'\n",
    "exp_fig4_id = np.where([np.array(exps)==exp_fig4_name])[1][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t\tE_Z \t\t\tFF_Z\n",
      "observed\t29.633605107882552\t284.17879402825855\n",
      "simulated\t29.07421144948377\t284.12615547072755\n"
     ]
    }
   ],
   "source": [
    "obsdata = dfs[exp_fig4_id]['N']\n",
    "\n",
    "a1_sim= 1.4\n",
    "a2_sim= 4.5\n",
    "b_sim=500\n",
    "\n",
    "a1_grid= np.linspace(1.0,5.0, num=20)\n",
    "a2_grid= np.linspace(3.0, 10.0, num=20)\n",
    "breakpoint_grid= [b_sim]\n",
    "\n",
    "n_sim = len(obsdata)\n",
    "simdata,p = broken_zeta(a1=a1_sim,\n",
    "                      a2=a2_sim,\n",
    "                      breakpoint=b_sim,\n",
    "                      size=n_sim,return_p=True)\n",
    "\n",
    "print(f'\\t\\tE_Z \\t\\t\\tFF_Z')\n",
    "print('observed',np.mean(obsdata),np.var(obsdata)/np.mean(obsdata),sep='\\t')\n",
    "print('simulated',np.mean(simdata),np.var(simdata)/np.mean(simdata),sep='\\t')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 4s, sys: 1.15 s, total: 1min 5s\n",
      "Wall time: 1min 5s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "means = np.zeros((len(a1_grid),len(a2_grid),len(breakpoint_grid)))\n",
    "maxs = np.zeros((len(a1_grid),len(a2_grid),len(breakpoint_grid)))\n",
    "FFs = np.zeros((len(a1_grid),len(a2_grid),len(breakpoint_grid)))\n",
    "alphas = np.zeros((len(a1_grid),len(a2_grid),len(breakpoint_grid)))\n",
    "\n",
    "for i_a1,a1 in enumerate(a1_grid):\n",
    "    for i_a2,a2 in enumerate(a2_grid):\n",
    "        for i_breakpoint,breakpoint in enumerate(breakpoint_grid):\n",
    "            \n",
    "            x = broken_zeta(a1=a1,\n",
    "                      a2=a2,\n",
    "                      breakpoint=breakpoint,\n",
    "                      size=n_sim,z_max=1000000)\n",
    "            \n",
    "            means[i_a1,i_a2,i_breakpoint] = np.mean(x)\n",
    "            var  = np.var(x)\n",
    "            maxs[i_a1,i_a2,i_breakpoint] = np.max(x)            \n",
    "            FFs[i_a1,i_a2,i_breakpoint] = var / means[i_a1,i_a2,i_breakpoint]\n",
    "            alphas[i_a1,i_a2,i_breakpoint] = means[i_a1,i_a2,i_breakpoint] + FFs[i_a1,i_a2,i_breakpoint]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "164318.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(maxs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_inset(ax):\n",
    "\n",
    "    ax.plot([10**0, 10**2],[10**0,10**-4],lw=LINEWIDTH,c='tab:red')\n",
    "    ax.text(10**1,10**-5.5,r'slope $a_1$',ha='center',rotation=334,c='tab:red',fontweight='bold',fontsize=LEGEND_FONTSIZE)\n",
    "\n",
    "    ax.plot([10**2, 10**3],[10**-4,10**-12],lw=LINEWIDTH,c='tab:green')\n",
    "    ax.text(10**2.8,10**-10.5,r'slope $a_2$',ha='center',rotation=298,c='tab:green',fontweight='bold',fontsize=LEGEND_FONTSIZE)\n",
    "\n",
    "    ax.vlines(10**2,10**-12,10**-1.2,linestyles=':',colors='tab:gray',lw=LINEWIDTH)\n",
    "    ax.text(10**2.2,10**-0.8,r'breakpoint $b$',ha='center',c='tab:gray',fontweight='bold',fontsize=LEGEND_FONTSIZE)\n",
    "\n",
    "    ax.set_ylabel(r'$p(z)$') \n",
    "    ax.set_yscale('log')\n",
    "    ax.set_xscale('log')\n",
    "\n",
    "    ax.set_xlabel('z')\n",
    "    \n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    \n",
    "\n",
    "    for axis in ['bottom','left']:\n",
    "        ax.spines[axis].set_linewidth(SPINEWIDTH)\n",
    "\n",
    "    sns.despine(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcQAAADXCAYAAACEaiGeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABR50lEQVR4nO3deVxUVf/A8c8wzLCpgKLmwiKm5oJLpunjhruoGBmGS6m54l5mZVmKmWmalbn7aJLmjua+Pb+ELC33StS0RFTcWBSQHWbu74+RcQYBGQXmznDerxcv78zc5cscnO+cc8+ikCRJQhAEQRDKOBtzByAIgiAIciASoiAIgiAgEqIgCIIgACIhCoIgCAIgEqIgCIIgACIhCoIgCAIgEqIgyFpoaCjPPfecucMQhDJBJERBkLn09HS6d+9OtWrVOH78uLnDKZLmzZsTHBxMcHAwJ0+eNHc4xerChQs0b96c5557jg8++MDc4RSJYXlkZGQ89ljQsS3sxR07dvDRRx+RnJxMr169WLFiRWnFJQjCQ1lZWYSFhTFmzBjmz59PWFiYuUMq1I0bN2jRogXLly83dygl4tNPP8XDw4PTp0/z4MEDc4fzRHnLw9rL51kUWkNUKBRoNBokSWLlypXEx8eXVlyCIBhQKBQoFApzh1Ekp0+f5uLFiwQHBxMSEkJkZCTffPMN3bt355dffjF3eMXCxsbG6F85MyyP8PDwxx5bY/k8rUJriCtXrqRDhw4MGDCATp06kZKSgpubW2nFJggCoFarefXVVzl37hw//vijucN5otOnT/PNN9/QrFkz/XP37t2jfPnytGvXzoyRFY+PP/6YN998k+eee47Bgwczb948c4dUqLzl8cknn1h1+TwLRWFzmX733XdMnjyZjh07smPHDq5evYqXl1cphicIgqXp2bMn1atXx9bWlu7du1OvXj127tzJhx9+aO7QyqTXXnuNzZs3Y2trm+/jCxcuiPJ5qNCEKAiC8CwyMjJo3749AwYMwM3NjTfffNPcIQkGRPkYs9iEGBoaytChQ80dhiAIgmAl5H9HuADR0dHmDkEQBEGwIhabEAVBEAShOImEKAiCIAiIhCgIgiAIwBPGIQoWLCMZ4v6Gyi+AfQVzRyMUUVRUFLNnzyYpKYmwsDA2bNhAeHg4mZmZLFu2DCcnJ/2+U6dOJS0tDUdHR+bOnVvoeV999VWqVq0KQGZmpv759PR0/XZaWprRMYavGU7vlZ2drd/WaDT5bmu1Wv123n57BfXje5r+fUU5plWrVqxatcrkc5eG/v37U61atceeN/y98v6Ohu9zVlaWftuwXFNTU/XbKSkp+m3DMs5b3oZlXNC24TUM/w4MtwFycnLy3Tb8uxg+fLjsykUkRGv01xbY/TZkp4LKCfy/gcavF7h7SEgIzs7OZGZm0rRpU3r06FFiof3xxx8kJCTQuXNnk46bNGkSr776Kr6+vixYsIB33323hCI0L29vb1avXk1gYCAAP/74I1u3bmXPnj1s375d3y3++vXrZGdn8+233/Lee+9x48YN3N3dCzxv1apVCQ4OBjCabuzevXv67bi4OKNjDF9LSkrSbxt+2Bp+IBtuF5Qc8z4u7IP/Sc8/6bVcNWvWfOI+5lKtWjUmT5782POG75HhewnG77NhsktOTtZvx8bG6rfv3Lmj3zacbSwhIcHovPfv39dvJyYm5rtteA3DhGr45Sm/x/nFLsdysb6EqMmGywegbg9QqswdTenLSH6UDEH37553dO9HITXFCRMmYGtry4cffoi9vT27d++mb9++7NixAxsbGyZMmMDUqVPp2bMnd+7c0f8nPnbsGD/99BPNmjWjfPny/Pzzz6SmpjJt2jSmT5+Ol5cXiYmJhISEALBnzx4+/vhjJkyYQKtWrYiNjcXf359ff/2Vtm3bMnXqVLp3745Wq2X06NGArkdxly5d+OOPP/D19aVcuXLcuXOnTKwCkTtdm6enJ+fOndM/f/PmTX0C9PDwICYm5rGEGBERQUREBAC3bt0qnYAFwYJZ3z3Eywdg8xtw+aC5IzGPuL8fJcNcWSkQd8mk03Tv3p2YmBj69evH66+/zpEjR3BxcWHgwIHcunWLEydOsHjxYry9vXFwcCA5OZldu3bh4eFBuXLlOHv2LG3atGHAgAHGoTz8hli+fHkGDRpExYoVkSRJ/y24Tp06jBw5khs3buiP2bZtG5cuXeL06dMkJyfj4eHB1atXn+LNsVzXr183+kZdo0YNYmJiAN1kzfl92/b19SUkJISQkBCqV69earEKzy537tonzV9ruF9RjxEKZn01xLo9IGg91O1u7kjMo/ILumZSw6SoLgeV6xV62KJFi8jMzMTX1xfQTVrctm1bvvnmG30Ncf/+/Sxfvhw3NzdatmxJy5Yt+euvv3BwcODy5cv6ml7FihV56aWX+Oijj7h586bRdezs7PTnB2jatCnr1q3j2rVrdOjQAaVSabS/VqslOTmZmTNncuPGDTZt2oQkSbz44ovP+EbJU0JCAtOmTePs2bPMmTOHgIAAxowZQ3p6OkuWLOG3337j7NmzjB07FpVKxeTJk7Gzsyu0uVQQhKKxvoSoVEH93uaOwnzsK+juGeY2m6rLQe+vC20uzW3OzM/8+fP1256envr7ULkaN25M48aN9Y9zEyrAwoULycnJ4bPPPtM/17NnT3766Sf9c02aNKFJkyb613Ofz/3XxsaGmTNnAuDu7s6oUaNYsGCBvoOItalUqdJjy/IMHDhQv926dWtat24NwJw5c0o1NkGwdtaXEIFx48bh6+tLv379zB2KeTR+XVdTjrukqxkWUy9Tw8RWVLa2tkYJ13CG/adlrR1qBEEwL+u7hwgMGzaMqVOnMn78eKNuwmWKfQVwbyGGXAiCIBSRVSbE5s2bc+bMGe7cucO6devMHY4gCIJgAawyIQI4OzuzdetWhg0bxv79+9m6dau5Q5KtkJAQvvrqK4YNG2Y0rslQaGgo//77b7Fcy3CgrqHo6Gj69+/PuHHjuHz5Mk2aNOH69ev8/vvv+Pn5ERERwdixY/nmm28eG5slFJ3omSgI+bPKe4i5cv+jV6lShfHjx/Pzzz/z5ZdfYm9vb+7QZEeSJFJTU3F0dCQwMJABAwYQGRmJi4sLrq6ugK7r/6ZNm2jXrp1+vOGsWbPo06cPnTt3pn79+vTurevQtG7dOn766Sf69u3L8ePHcXNzo2bNmkRFRbF27Vratm3LwYMHuXTpEosXL9bH0aVLFxo3bszp06fp0qULO3bsIDMz06hXqYuLy2O9UQVBEJ6V1dYQDTVv3pzTp09z584d/WBvwdikSZN4/fXX+eOPP6hbty6vvfYaycnJTJo0iT///BOAL774gsmTJxuNN4yNjaVZs2ZMmTKF48eP68/XokULGjVqROvWrbly5Qqurq7ExcXh7e3N4MGD9VM9xcbGGs1+8X//93/s3buXnj174uDgQGpqKvb29voaTN++fcU6mIIglIgykRBBV6vYunUrX3/9NcnJyezevdvcIcnKt99+y+HDh/Hy8tKPEXR2dmbhwoX6YRGfffYZn376KR07diQmJoaKFSvi5uaGUqlEoVAYTaM1atQonJ2diY2NpW7dujx48ICGDRtSp04dFi1aRGRkJPb29kiSZDRNVZcuXZg5cybOzs6Abr7DwYMHl+I7IQhCWaWQnmZGXRnInYHjaVy6dImePXvSs2dPvvzyS/1gcUGwVsHBwfoxpIb3iQ3nKzWc/zLva5Y8l2m9evWe+rOipL3zzjv6aRAN7+Mavkd577kbvs+GZWHY0nL37l39tmG5Gm6X1lymhj39DWOfMWOG7MqlzNQQDdWrV4/Tp09z69Yt2rdvX2AnD0EQhPycPHmSefPmERoaau5QBAPPWi5lMiGCrgk1LCyMpUuXYmtrS2RkpLlDEgRBJqKiohg+fLh+1ZHU1FSGDBnCyJEjWb9+PS1atDCqQQulo6TLxap7mT6JQqGgefPmpKSkEBAQQI8ePViwYIFFNqGGhITopzgDOHXqFAAvvfSS/rncJoqVK1cyatSoUo9RMC8xtEI33V3u+n63b9/WT3ru6+trNO1g3mW4tm/fTmBgIP7+/gQFBTFo0CDmzp1r1EP6aRQ05KWslZVcykV2CfHIkSMsXbqUTZs2ldo1y5Urx6lTpxg+fDj/+c9/2L17t8WtDlDQPdX87rHklwwzMjJYtWoV48ePN+m6kZGRfPzxx3z55Zc8//zzRq9t27aNiIgIFi1aZNI585PfuUJDQ4mKiqJp06b07dvXaP/Y2Fh27tzJyJEjCz2v4Tnc3d3Zt28fkiQREhLCyZMn9Y9feOEF2rdvb3F/F4KxtLQ0pk+fDsDs2bOLfA8rJiYGHx8fAJRKJbt27eLMmTP5Lu4rmE4u5VKiCTHv6t+pqamMHTsWtVqNr68vgwYNeuyY9u3bc+zYsZIMK1+5Tag//PADFStWJDk5mQoVys60Z3v37qVz587ExcWxbds2fv/9d8aPH6//F+DPP/8kPDxcf8yECRNo1KgRAQEBj53v3LlzODo66nuLGjK8xogRIzh27BiTJ0/G1tY232tcuHAh33NVqFABOzu7fKfnO3XqFNHR0YwYMYLZs2dz586dx86rVCqNzrF7924+/fRTFi1aRFxcnNHj1q1bs2HDBqZMmWLyeyvIh1arfao+AzVr1iQmJoamTZui1Wrp06cPffr0KYEIyya5lEuJJsQnVW+dnJzYsWMHAH5+fgQFBRV6PsMFT6Ojo4s9XoVCoV+R/JVXXqFu3bplphdqVFQU3bt358yZM6jVanJycnjppZee+t7q4cOHkSSJEydOcOPGDapXr64fTH/x4kX9Ndq2bcuJEyewtS34T9HwXIYLA+fWCqdOnUpgYCA2Njb6a/z555988sknbNu2zWiF+LwMz6FWqwvcz9HR0ainZVljLU14Rf3gzbsM18SJExk/fjx79+7F39+/FCJ9pLD3vqDXLK285FIupdpkmrd6GxAQ8Fjt4ty5cxw9epRDhw7RrVs3o9cM25OLvbuuJlu3uHDdHqBUsXnzZoYPH06bNm3YvHkztWvXLt7ryYy3tzc3btzg2LFjPP/88/o1Dv/8808yMjKwt7d/bKkm0C1Oe+jQIf79919mzpzJ4cOHadiwIZMmTQIgPj6eAwcO0KVLF2rVqgVgdI0vv/wSNzc3srOzUalU+V7D8FzPPfcc//vf/2jYsCGRkZGcOnUKlUpFaGio0TViYmJYvnw5Wq1W35Sb97wAhw4d0p+jd+/efPrpp0iSROXKlY0eOzk55VvbLU3r16/n5s2bBAcHc/z4cbp27WrWeCyRVqvVTwpRmPyW4VqzZk1JhVXmyaVcSjUh5q3e5sfHx8c8g+YvH4DNb+gWF67fW9+EunjxYm7duoW3t7fFfesyRa9evVi1ahVTp04F0Nfqv/7660KPc3d3Z8OGDfrHeT+kP/vsM2JjY6lSpYr+ubzXKKrc5adyr1G9enX9l6a811iyZEmRztmtWzejL14tW7Y02s59vGnTJqN1Cc3hypUrTJkyhfnz55OcnCwS4lPQaDRGY+EEeZBLuZTosIuEhASCg4P11du+ffuybds2xowZU+rNDk9Ut4cuGdbtrn9KoVAwYcIE2rVrx9y5c5kwYYLVLidlb29vcoeaojJMVCWlpK/Rv39/s3eoadGiBY6OjkyfPr1U3tOCWPKE4LkfvHL48M3rWd9XUS7PrkRriBbV7KBUQf3eBb48ZsyYMtWEKsjLrFmzOHbsGNu2baNZs2Z88MEH5g7JImk0miI1zQmlSy7lUmYH5psqtwl1yJAhrFixwtzhCGVMQkICrVq1Ytq0aVy6dMnc4VgsjUZDZmam1bb0WCq5lIvsxiHKWW4TKsAff/xBaGgoX3zxRZnohSqYl6urKxqNhipVqlCxYsUSuYYlNrWZKicnRz8AXC5Mfd+tsZzkUi6ihlgUmmy4uFv370NeXl7cuHGDNm3acOXKFTMGJ5QFM2bMIDg4mIkTJ5q9t6sly8nJkUVNRDAml3IRCbEocnugXj6of8qwCbVfv34F9poVhOJSvXp1Vq9ezTvvvGPuUCyWXJrmBGNyKReREIsinx6o8KgJ9bfffkOj0fDVV1+ZvUAFIT+GPRAtuTfis8rOziY9Pf2x5YrkoKByKazMrKUs5VIuIiEWRW4PVKXK+PmHTal2tjZkZGTw66+/iiZUoUQsXLhQPy/rrFmzzByN5crJyZHFB69gTC7lIhLiszBoSi1fvjzbtm1j8ODBtG7dmps3b5o7OsGKXLlyBXd3d4BCp6ITCpednU1aWprR4raC+cmlXERCfBZ5mlIVCgUTJ07k5MmT1KhRg19++UU0oQrFQqFQkJ6eTmRkJLdu3TJ3OBYrKyuL1NRUo5XmBfOTS7mIhPgs8jalPmxC9axZHUmSWL58OW3atCEqKsq8cQoW791330WSJNatW8fnn39ebOe1lntQRSWXprlcRb1vWJR7wIbP29jY6H8sgVzKRYxDLE4G86Eq6vfmhx9+YNGiRbRq1YqwsDDat29v7ggFC+Xh4cHcuXPNHYbFy+28IciLXMpFJMTilLcJVZvDxK61aN1yJ7Wer0NcXJx+/T1BMMWCBQs4deoUOTk5tGzZkvfee8/cIVmk3JoIYDG1p7JALuUi/iKKU94m1Ic1xhbOCbi5ubFkyRLRhCo8Fa1Wy8aNG9m6dSt37941dzgWK3dGFDnMilKYwpqyTR1CYwnN4nIpF5EQS5JhjVGTzYzXX2TwG4No1aqVfmFkQSiK8+fPs2HDBtauXcvt27fZt2+fuUOySHLpzSgYk0u5iCbTkmS4gsbF3Si2vMnEft/TuvpUbkkaJEkiJycHlUpV+HmEMq9jx476FcW7detGfHy8mSOyTHKZM1MwJpdyETXE0pJbWwRaXJjJKy/YsX//flq3bi2aUIUncnV1Zc+ePezduxcXFxcGDx5s8jlMnanGEpraTCWXpjlDpvQkLep5LI1cykUkxNKSW1t8oZcuMdbuhJ9XDm8OGkirVq3Ytm2buSMUZGzPnj1s2bKFzZs3c+DAAXOHY7HkMkWYnqTF7ur/Hi0coMnWPdaaf23A0iSXchEJsbTlJsYrP6HY8iaTetRl766dHNr8X6Qc+a3iLchDeno6165d49q1a2YfvGzJsrOzZTEAPJdtUjQVD03A7vrPANhdj6DioQnYP3xcVsilXMQ9RHMx6HDT4vIBWjT4jdjtHzHku0iWLF2Kt7e3uSMUZCQkJITFixcDMH369GI7r2HzmiU2tZkqd0YUudA4e3G/+xKyPDqgUCjI8uzI/e5LyHBvp9/nacqloGPkWsZyKRdRQzQXwyEadXtAm7epfH4lPZrWoFWLF9ketsXcEQoysW/fPi5duoSvry++vr5cvnzZ3CFZrKysLB48eCCz+WClJzy2fnIpF1FDlAOlCjp9jKJmCyZpc2idsJmR0z6gffVs3F5+/fFVNoQyJS4uDtB9u9dqtWJA+TOQS00klzIpGteD33G/+xIya3XB7noErgfHI3VbRIZXF3OHV2rkUi7if5ZcGHS6aTl5E2fDFlDp4BhChnYh6p9L5o5OMKO4uDiGDBlCRkYGYWFhnD9/3uRziPUQdeQyZ2YujUstEnssIcvTV99kmthjCZkevgWWmeE8pYY/ps59KidyKReREOXmYWK0aeAPbd7GNe44rVq8yLatm80dmWAmMTExAJw8eZI9e/aY/UPDkpXkAPAjR47Qv39/0w+UpPwfl6Eep3IpF5EQ5UqpQtH5Eya9+wF7+yl4b/Ik/rf600fds4UyIy4ujoMHD1K1alUA/QB9wXRF7d4fFRXF8OHDCQwMBCA1NZUhQ4YwcuRI1q9fn+8x7du3p2nTpibFY5t4FZeD47G7HgGA3bXwh49/1vc4dbh+xKRzWiK5lIu4hyhnD+8ttqjZgrPDkyi/fww/f/UP7j5t8O46XNxbLCOmTZvGoUOHmDp1KgDdu3d/qvPkNpnl/tutW7fHOjE4OjqyaNGiZ4hW3j777DP94O9bt24REhICoO+wlMvb25vVq1frP3i3b99OYGAg/v7+BAUF4eTkpJ9+0c/Pj6CgoKeKR+PqTZJfP7I9fbGxsSG7VmeS/JZRu8dYEhOTdDvNfAtnZ2d9U3lBPYMtucewXMpFJES5e9iE6qzJhgeXufztPPqFrGf5lMP0nb5eJMUyoEGDBjRo0ED/+JVXXimW8z548IC9e/ei0OZQJfEsf2VWZ9iIUcVybrnK/VJhqpiYGHx8fABQKpUEBAQQEBBgtM+5c+c4evQohw4dolu3bkU7scKGLO+ujx4rVWR5dyUxMYk7N29gdy2c1JrtqOlR66nithRyKReREC3Fw9riyOrNaPLTLoJmhBJ9+z6TJ43TzX4jEqPwlKoknqX5v19x322YuUORrZo1axITE0PTpk3RarX57uPj48Pu3buL7Zq5zaeart8W2zmtTXGXi0iIlkSpgoYBtHyhF2dqupF0ZBnJPxwnof5Qag1aIJKiUCDDXoZ5m9biXF/kdJ13uZlVA3i0Hp2lNbsVp4SEBKZNm8bZs2eZM2cOEydOZPz48ezduxd/f/9iu46NjQ1KpTLf13KbT7NqtNXva/gvYHSs4fOG29ZUjiVdLiIhWiKlCtdXZuNa7z/s376RIWMWs+LKLV6dvlEkRaHIKlSoQM+ePY2ec3JyMlM08lKpUiWWL19u9NyaNWtK7fouLi5UqVbD6LFQ8uUiEqKlelhb9AP2XN9L0Lfb+TWyHQs+nSqaUIUiOXToEPfu3aPqvRO8+M9XHK48nOtOjc0dlgBGsxGJXsWlRwy7sHQv9KLl2z9wZtVk2tj+BVuHkBTaH86FiSEaQpHEuuiaTG84NjR3KIJgViIhWrqHNUXXgDn0nb6BG26+1Jv0Iz9+Nhi2DoUssTK4gNEMJ3mf1/08PptN/vtZz/0oOVAoFCiVyif/oMXh2k8o0Ro9bzhTTUHPW9qsNeYkEqK1eJgY3UdtZNfnQ5kcbsOkr8PI/K+fqC0KhaqceIYXLy+gZprpU8IJpUN9LRyXA+NQX4swdyhWTSREa6N2pOXE7zhz7hIq56pob51FChsmaotCgeJcXuRM3SnEiCZT2dLNcbqULE9fc4di1URCtFKuVWvw5Z6/UbeZQLtNdvz4448Q2kvUFsuw/CaH1jehUfAE4KJ5reQolUpUKhUqlQpbW1v9T+5zKpUKtVqN2sEJXuiF2sEJlQ2Ui4lAZYPRfoZNpobnKuoE4IJIiNZN7Yiy52wWrA5j8k8Sk1YfI3PzMPhppkiKgl7l+2dodvlLaqZFmjsUoQhyl4iyu/6zuUOxOiIhlgEvt+3AmfP/clOqyp93tXBsEWx+UzShCgDEub7I2XpTiHFsZO5QhCLI9PDVrZ/o0cHcoVgdMQ6xjHCtWoOw367AlsF8vupH6l/cxaup/jBkN6gdzR2eUIiLFy+ycOFC4uPj6dy5M2PGjAEgMjKSOXPmAPDhhx/SqFHBCS23qRSMZzextbUFbLlX5T/Y3Lnz2DH5bQvFx9bWFgcHBwCjqcc0Gk2+23r1e2MP5Egp+qfs7e3123Z2dkbXyCXKsXAiIZYlShW8tpoud+0I+mwzP1/7lXlJrVA37Amdp4vEKFP169dn+fLlaLVaBg8erE+ICxcuZMmSJSgUCt5//31WrFhhdFxERAQREREA3Lx5s7TDFgSLY3KT6bx587h3715JxCKUBrUjLd8O5cz5f4lOc+CrA5fh+DL4TnS4kbNdu3bRq1cvo6nWkpKScHFxwdnZ+bFlnEC3dE5ISAghISHUqFHjsdcFQTBmcg2xa9eufPLJJ9jY2DBixAiaNGlSEnEJJcy1ag1+/D2anNA+nP3zDNf+Pk7AnTO6pBj4nagtykyfPn3o06cPvXr1YuDAgQA4OzuTlJSEQqGgfPnyhR6fnp5e4AKqRWE4oN/Z2fmpz2MOMTEx5g6hQDdv3mTmzJlmuXaVKlUKfVzS5FguCkmSJFMP+umnn1i5ciUVK1ZEpVLx7bfFtzzJyZMnCQ8Pp0qVKgwdOrTA/XK/+QrPSJPNqa1f0W98CAHPa/miqx3qmi/CW/tEUpSJiIgItm/fTmZmJo0bN+b3339n3bp1REZGMm/ePADef//9Qu8hCoLwZCYnxAEDBtC9e3cGDhyIWq1m7dq1DB48ON99o6KimD17NklJSYSFhZGamsrYsWNRq9X4+voyaNCgfI+bOnUqL7zwgkiIpej+3Zu81b0JNplJbA9yBJda0PFDaNRXTBQuCEKZYHJC3Lt3L7169QJ09zX69OnzxGMCAwMJCwtj3bp1uLi44O/vT1BQEAMGDGDHjh0A+Pn5ERQUpD9m8eLFjB8/3ug8hp0EoqOjCQ0NNSV04QmkzFRurhhA1bgIfrmuoVMtW3g5GLp9JpKiIAhWz6R7iHv37mXjxo0oFAokSWLDhg1FSoi5YmJi8PHxAXRdvwMCAggICDDaZ9euXZw5c4Zq1ao9dryvry++vr4AonZYAhR2TtQct42rB5cxbMm7vFo3hy80y1Dfuwr9QkUTqiAIVs2kXqbx8fE4ODgQHx/P/fv3+eijj0y6WM2aNfU3Ug3H3Bjq06cPISEhjB492qRzC8VEqaJWz4mcOfcPUQ9UtFuTStqFA/BNYzi7XvRCFQTBapmUEGvXrk1AQABubm64uro+sZdQQkICwcHBnD17ljlz5tC3b1+2bdvGmDFj8Pf3f6bAhZJVsYYXO07GMLuPF44qBddv34WdY2HjADHDjSAIVsmkJtP79+9z//59ACRJeuKsB5UqVWL58uVGz61Zs8bEEAVzUTi60OWbc2TvnUa3UYvoUVvJvK6HUH/nB8P2iyZUQRCsikk1RBsbG9zc3Ix+BCundkT16tccO32Rq8k2tP0uleuXzsDXPqIJVRAEq2JSDTE+Pl5fKyxKDVGwHhU96rDj1G2WDaqLg206mQ/isNs5FrQaaJ7/sBtBEARLYlJCrF27NsnJyeSO1BAJsWxROLowdsM/8L0/A5b8RrVyNszVTECt1cCLb4ihGYIgWDSTmkzv379PXFwcCQkJJCQkEB8fX1JxCXLl4Awjf2LJwoVcua+l3ZoUotdPhG+bQ4r4exAEwXKZlBD9/f158cUXOXbsGL/99hsvvvhiScUlyJlSRcVOY9lx9F/6N1Rz4qYGkq7Bl3Vg9zuiF6ogCBbJ5Mm958yZw6JFi5AkiUmTJj3ThMGCZVNU8uSdLVGwpCWhJ+L5446WeZrVqP/9Pwj+VVebFARBsBAmL/9Ur149nJ2dcXFxoU6dOiURk2BJnKvBlEv06dqRq4la2q1JJfpaNHzhBXf+Nnd0giAIRWZSQuzXrx+nT5/G398ff39/zp49W1JxWQQxl+pDakcqBu9kxw+rCGpoz/TwTEALy18WTaiCIFgMk5pMt27dWlJxWKTo6GhzhyAfShWKl4Ywea4t0o4x3EjSsvxUFjM0q1Hfi4KX3oIXeomeqIIgyJbJTaZz587l5ZdfplOnTjRu3LgkYrIomgcPuD5iJCm//MJTLC1pfRq/jqL/Bpyeb8u52IdNqGcPw9YhsH+qGMgvCIJsmZwQ79y5Q48ePTh8+DBDhgwpiZgsio2jI86vBnB37hdc6z+AlF+PPlNiDA8PJzw8vESPO3v2bL7N3XnPkZKSwunTp00LRKmC+r2pOGwTOye2IKihii5rU8nSSHBqFZz63rTzCYIglBKTe5m6ubmh1WqZNWsWf/31V0nEZFEUSiXOvXpRoUcPkvcf4O7nn6N0dsZt/Dic/vOfp5q84J9//iEyMpKcnBzs7Ozw9PTEy8uLX375BaVSyWuvvUZMTAxnzpwhLS2NN954A4C0tDS+++47lEolGRkZZGVl0bZtW1QqFT///DM2NjYEBgaSmJgI6BLgpUuXyMzMpGfPnvz1119otVpcXFxo1qwZt27dIjExkZ07d9K5c2fKlStX9F/CwRlF8M9MfvE7hu98D5UNrP0zi/6ayahvnIB63aDBK6IJVRAE2TA5IX788ceArpYxadKkYg/IUimUSpx796KCXw+S9+3n7mezUbq6Unn8OBxbtzYpMXp5eVGjRg3Onj1LpUqV8PPzY/369Wi1WiRJ4sqVK2i1WnJyckhJSeH27duArkz+8zAJZ2dn68+hUCjo1KkTV69e5fLly0bX8vb2xtbWlps3b+qbwJs1awbA3bt3ad++PRcvXiQzM9O0hAi6ZNd6NM4+r5HxVXPCLtxiycksNgduxCtyMxxbBIN3ieEZRdCqVSt69OhhFeuAvv3223z++ec4Oj6aHP706dOcPHmS4OBgM0Zm3dLS0pgxYwbJyck4Ojri7+9Pp06dzB2WrJicEOfPn8+JEydQKBS0aNGC9957ryTislgKpRJn/95U6OlH8r593Pl0FspKlXSJsVWrIiXG6OhoLl26RE5ODjVq1ACgadOm/PLLLzg5OVGlShUuXboEQPny5fVrS7788stcv36dpKQkHBwcuHTpEm3btsXW1pbDhw+jUCjo168fkZGR+mupVLoamiRJVK5cmYiICDw8PPTT9J06dQpJkqhUqdLTvynl3LCfepGd1QL4OuwIL69K5chQR+rxB6ztAyP+T9QU83HlyhV69OhBxYoVuXv3rrnD0QsNDSU8PBwHBweqVatGdnY2kZGRzJw5k82bN3P37l1GjBhB69at+fHHH7lx4wYNGjTgxIkT+nVO4+PjWblyJVevXgVg/fr1fP311xaREM+fP8/QoUOJjY3Fx8eHLVu2GCV3uVq2bBl+fn54eHgwcuRIhg8fzjfffMP+/fv5+OOPadeunblDND/JRFOmTNFvv/vuu6YeXmxmzJhhtmubEoM2J0dK3LlT+rdbd+nqoEFSym+/SVqttkTjOnz4sHT48OESvcZTycmSpFB/6dwYJ0kzvbx0boyTlPlxeUla9LIkpSWaOzrZeffdd6U2bdpI2dnZkqenpyz+5iVJktasWSOtX79ekiRJ6tSpkyRJkjR79mzpxIkT0jvvvCN9+OGH0pgxY/T7jxs3TgoMDJRycnKk3bt3SytWrJAkSZKioqKkN998U0pJSZEkSZKGDx8u3blzp5R/G9MFBQVJzs7Okqenp+Ti4iLt2LHD3CEVyaBBg6SsrCxpy5Yt0tKlSyVJkqSff/5ZWrVqlZkjkw+Taoh79+7l9u3brFu3DtB1sBEKp1Aqce7Thwo9e5K8dy93ZoRgW7kybuPH49Tq5RK5ZseOHUvkvM9MqYIBm2iU2hViI1nwWxYX4jRsDjyP1+ru0HGqGJqRh0KhkOUk+hUqVACgcuXKAKjVaj766CNWr16NJEnMmDED0LU8JCYmYmNjQ05ODsnJybi4uBAdHU1ISAjLli3DyckJAFdXV5KTk6latap5fqkikiQJPz8/Nm7cyLlz5/Dw8DB3SEXi6+vLRx99BMDNmzfp0KEDR48e5cMPPzRzZPJhUi/T+Ph4unbtiiRJSJJE165dSyouq6OwtcX5lVfw3rsH58DXuD1jOtfeHEzq8RPFep2DBw8W6/mKndoRRvwPmg7huz72BDVU8fKqVH46cU43NOPCLnNHKBvBwcHcvn2bFi1a6JvF5axjx4588cUXfP/9o57ECxcuZMCAAbz//vt88skn1KlTh+joaHr27EnFihWZM2cO9+7dA3Qf0paQXKZNm8aFCxeoWrUqQUFBZGVlmTukIhkxYgTz589n/vz5fPfddwwdOhR7e3t9BUcAhSSZNkbgypUrLF++HIDRo0fz/PPPl0hgTxISEmL2DgbPEoOUk0PS7j3EL1uG6rnndL1SW7Z85piOHTvGf/7zn2c+T6k4txW2jeD3mBwqOiioXt4GtVKBOvgo1Ghk7uiEEiBJEmPGjNF/huRKSUlhypQpjz0vCKXJ5HGIM2fOZOLEiUycOJFZs2aVRExlgsLWFpdXA6i9by/OAQHcnvYx14YMJe3kyWc6r8UkQ4AGAfDaalrVrkzdSkpWnMrSLSc1vzX8tlwM4rdCCoWCQYMGkZZmPJ3frVu3RAc9wexMToh2dna4u7vj7u6Og4NDScRUpihsbXHp+6ouMfbpw60PP+La0LdIO3Xqqc735ZdfFnOEJUipAp9AmHAK1FWZ3Fqtb0LduXAyfNdTzINqhdq1a/dYr8y6detSu3ZtM0UkCDomJ0R/f38CAwPp168ffn5+JRGT1ZI0mgJfU6hUuLzWl9r79+Hs35tbUz/k2ltvkXbmjEnXGDVq1LOGWfrKucEH51HU7szk1nbsCHLgdooEN0+gXRcI53eI2qIgCCXOpIQoSRLR0dGEhYWxdetWXnnllZKKy2rkxMVxf/MWcuLiUCiV5MTFFbq/LjG+pkuMvXpx6733uT5sGGlniraySEJCQnGEXfqUKui/Htwa0drdluCX1ByPyaFtyCGu/fdN+EtMLC8IQskyadiFQqHg1KlTbNy4EWdn3ewiPXv2LJHArEX27dukHf+dBwcPUr5rF9LPn6dc+/ZU6NZNP+dpft3qFSoVLoGBOPfpQ+LOndyaMgW1lxduE8bj+HAmmfxERERQq1atEvt9SpTaEUb9D7YOg3/207KGksAGtrRclcp/Y0fTZ1lH3fqLViwqKorZs2eTlJREWFgYGzZsIDw8nMzMTKMhCgBTp04lLS0NR0dH5s6dW+h5X375ZVxcXACMeqwa9qnL25PV8LHhtsagpaOgfQraHzDqlZmZmanfTk9Pz/d5w/2zs41bCgzPnZOTk+/1hw8fzqpVq5CjgQMHUr16daDgsshbLoa/s+H7YfieGd6jTU1N1W8X9B7nPVdBZWwYY0HbeRX02ssvvyy7cjEpIc6fP599+/bRsmVLi+lqbE6a5GSSdu9BYWePQ1MvNImJqKpURVmhAtq0NGwe3kcpNDGq1bj264fLK6+QuGMHt96dgrpWLdzGj8s3Mb711lsl+0uVNLUj9F8H619HEXWYya3taF1TydyjGfT6qgHKyResOil6e3uzevVqAgMDAfjxxx/ZunUre/bsYfv27bz55psAXL9+nezsbL799lvee+89bty4gbu7e4HndXFxoXv37oDxh11BCSXvY8Ntww/OgvYx3M77wWv4AZ2UlKTfvn//fr7PG+6ftzOO4bkLSig1a9ZErqpXe46PX2tCpmdHtAql/vmC3ksoOPElJyfrt+Pj4/XbsbGx+u3cIS559wfjZGn4+W54/aIkyryPC9qWY7mY1GR64cIF4uPjOXfuHEOGDBGrXTyBskIFXAJfA0ki4+LfaDMykbKzSd5/gGtvDubWR9PIvnu3SIOvFWo1rq+/Tu0D+ynfrRs3332X6yNGkv7HH0b77dplBeP4lCpoNlD/sLW7LTv7O5KerWFQ+1qk3ys7E0Lk/l14enoSExOjf/7mzZv6BOjh4WH0Wq6IiAj90CDDD0JBPmwTr+JycDx21yPMHYqAiTXE+Ph49u3bx927d9m3bx8gmkyfxL5eParPnUP2rVuknTyJJiUFbUYGlUaPImH1am5/NA1V9Wq4BPXHoVFD/XGSJBVcYwx6HZdXA0jc/iMx70zG7vnnqTx+HA5NmuibXyxegwCwdYDkO7D/XQCcVBCfmsPmUQ0YOnwUdJ6uq1GWAdevXzf6Rl2jRg19Erxx4wYBAQGPHePr64uvry8Av/32m/7vycbm0fdgw2/5ef/eDB+bui0UTY5LLRJ7LCHTw9foefFemodJNcTAwEDi4uJ49dVXiYuLI+4JHUQEkB5+4KiqV8f5lVdwefVVtGmpxH65AO2DFCSNBoW9AwkrlpN17Rqah0szFanG2D+I2gcPUL5zJ2LefofoUSOx01hJ7enhuoq8NAQ82wK692RiSzVLf0uC48sgbJjV9T5NSEggODiYs2fPMmfOHAICAhgzZgxbt26lb9++/PbbbyxduhQPDw9UKhWTJ09GqVQW2lwqyJxYWFw2TKohiiZS0ykefhvPrfHZODpSffZs0s6cIXHLVjKvRoFWC7a2JKxareuNqlJR9eNpqIowp6ONWo1r//449+3Lv2GhbD7wO9Ubt6aSwzOsTiEnShW89BZc+xWAHs/bsuRkNvfSJSpe3g9/bITmg80cZPGpVKnSY7O1DBxo0HzcujWtW7cGYM6cOaUam1D8dE2mq0nssYR0z87mDqfMM3kcovB0cmt8klaLpNWirlULu7p1cX09CJV7TZx79ULSaMi+GUPGpUtkRV/L9zwF9diyUaupO3AUdp3tWHthbYn9HmbR4BV4dQXYVUJpo2DfIEcqOjysQf88z7yxWQgbGxuUSiVKpVJ/z1qhUGBjY1Pgj+F+Bf0Iz6agJlND4j0vPSIhljKFjQ0KGxtsXV2pNOwtXF7rS6WhQynn64vjSy+hzcjEwccHhc2jP3xNSgrZD1cWedJ/iC7luxB2OYz7GfcL3c+iKFXQpD/01NWINFqJ1qtTSUjTQnIMnF1vdU2nQhkimkxlQyREM8m9t5g7e03uFG5eWzZTecJ47Js0ASAzKoqMyPPEL1/OvbXryLlfeKK7ev4qPbx68P357wvdzyI16gu1fFHaKKhT0YY1f2QDEuwcC6esrFYslAmil6m8iISYR1xcHFeuXHlsAHBxy723qFDqxh5JkoSk1WLr6oraywsbtRopO5uUw4fJib2LJikJxxYvkXbiJOnnIgs878CBAxnhM4Kwf6yslggP11PcCLYOjGuhYtmpLLS53673T4Z7180bn4wV1DxaWPNpQc2kBe2T93qiie/JitJkmpd4z0uOSIgPrVixgsGDBzNt2jRWrFjB0KFDGTt2LFevXi2V6ysUCn2S1D+nUuE6eDBqDw/svGujrFQJTWIi9vVfKPA8YWFhVCtXjW6e3azvXiLohll0n0PLGkq6etsSn2bQ3PRfmS6MLAiFEU2msmFSL1Nr1qJFC3799VcWLlyIg4MD9+/fR61Wc+vWLbPFJEkSNmo1Dk2bYu/jQ/pff2Hj6IjCtuBiq1u3LgAjfEbw+p7XGdJgCC72LqUUcSl58Q0U8ZdZrliKRmvwYZIeD7vfge6zy8z4RMGyiV6m8iJqiA+9+OKL+sVL09PT+frrr3FycqJOnTpmi0nfM1WSUCiVODZrhrN/70KPady4MQDVy1Wnq2dX66wlKlXQ7VO0ryzDZ1ka0YkGcz2e/k70PBUsxtM0mQolRyREA7Vq1SIkJIQxY8YYzZ9obqbcEwgJCdFvj/AZwZbLW0jMSHxsv7g03aQKhU3KK2tKFTbNBtKjcVVWnMozr27r8eaJScYM7/3lDr9QKpWF3jM09V6hqfuIe12Q+w6I90UeREI0MGrUKLy8vAgJCeHAgQPmDuepGCbEGuVq0MWji1Et8XbKbU7dOcUXJ78g5kEMCoXCcpMiEDzhHVafzSYjx+B3CBshhmEIJeqbb77hq6++4oyJ65XmpUyKxuXAONTXIoonsDLu5MmTzJs3j9DQ0Kc6XiTEhzIyMvTTX3l5eXH+/Hn985bkr7/+Mno8svFItlzeQlKmbvWAg9EHORB9AI1Ww8Hog4Blz5tY138Sb/drT1KGQUKMDodTVjjsRCg1UVFRDB8+XL/qSGpqKkOGDGHkyJGsX7+eihUrkp2d/dhKFKbSOHuR2GMpWZ6+xRC19XtSubRo0eKZJrIXnWoe+vDDD1EoFDRo0AAXFxeio6OJiopi0KBBtGnTxtzhFdnly5f19xHhUS3x+/PfM/HFibSv2Z73j7yPr7svKdkpfHvmW7p4dqFBpQZmjPoZKFV8tPp/JM/0BAyauQ9O1a2YITrXAI+aTMG4mdxwom+bPL2cDR9byyTec+bM0X/JvX37tn4yfMOJ0OHxZbi2b99OYGAg/v7+BAUFsXnzZgA+++wzWrZs+YxRPd5CU9T32JLLwlBxlcugQYOYO3cuixcvfqo4REJ86Ouvv+bevXvs27ePnJwcfH19mTJlirnDMlnuH4qhET4j6L+3P4MbDMbbxZve3r2JTY+lskNlYtNiiU2LtdyECEg2trT8LoONvTQ0q/ZwTTlttm7y76B1uk44goBu/cBPPvkEgM8//9zoFkNhYmJi8PHxAUCpVLJ//35OnDjB888//0zx6JpMV5PYYykZXmW3l2lxlcuuXbs4c+YM1ao93ZqpsmsyLa62+acxatQojh49yq5du3jw4EGpX784bNiw4bHnapavSWePzqy7sA6AV+u8ir3Snti0WGo516Llc7pvuFpJ+9ixlkChUDB49CSWns3zbfnyfrhgBetDCsVGq9WSk5NjclNnzZo19cttabVa/Pz8mDFjBoMGDXqmeESTqU5xlUufPn0ICQlh9OjRTxVHiSbEJ7X35qe42uafRr169Vi2bBkbNmwgLCys1K9fHJo3b57v8yN8RrD58mYSMxJxtnOml3cvNJKGmJQY3g5/GwAbhey+HxXZ8FHBhP2tJTEjT/NTzEnzBCQzNjY22NraYmtra9R7tKAep0X9sbTeo0X94M27DFffvn3Ztm0bY8aMwd/fv/gCslGS5d210FYMS3lvn4VcyqVEm0yf1N7r5OTEjh07APDz8yMoKIjBg3VL+RRP27xpDhw4QMWKFWnSpAmOjpZ578nLyyvf593Lu9PJvRPr/17PuKbjqOJYhZN3TvJWo7c4evMof8T+QdMqTQtcmFjuqlatysJvviHr6mzAYMq6uxd1PU6toNl0/fr13Lx5k+DgYI4fP07Xrl2LfrCkxTP9HNftLbdpvDhoNJoiTcuY3zJca9asKf6AtBrUUYfI8uwICmXxn99CyKVcSrVKEBMTo+/JqVQqCQgIIDQ0lNDQUIKCggDYv38/M2fOpFatWo8dHxERQUhICCEhIURHRxd7fAcOHKB+/focPXqUmzdvWuT6jwsWLCjwtZGNR7Lp703cz7hPeXV5etbqyYaLG4hNi+Xw9cOAZd+kH9z6OdIeJDya3xQgOgL+3mu2mIrTlStXGD9+PF9//bXJw4Kcc+LomhCKR8aFEorOMmg0GrKyssjKynryzqVADLvQkUu5lGqnmtz23qZNm6LV5n+/ys/PDz8/v3xfM+xxVNSbrqaoXLkyPXv2pGfPnsV+7tLy0UcfFfiae3l3fN192XBxA+OajcPZzhmP8h762mLY5TAC6wZabC2Ruj14dZvE/I5aungb/Gknmm/6veLUokULHB0dmT59OvPmmTYbT7KqCv/n9h9iHBpga/Bfz7DHqebhyiu5itLLtCjbcpL7wSsXGpdaJPktI9uEe4hyfW+fhVzKpURriKXWDi/onTp1qtDXR/mMYtOlTSRnJVNeXZ7UnFQ0koY+z/chKikKsOD/cEoVwYNfZ8nJPP+xDn9qnniK0axZs/j2228ZMWIES5cu5YMPPjD5HAqD7v0KSYNX+jlsJE0hR1gfudREjFnuxBjFRS7lUqIJMbe998qVK3z44Yc4OTmxZs0ali1b9sy9s4T8PWkycvcK7rSv0Z71F9bjV8uPxZ0W85zjc5y8c5KK9hX1S0ZZ6uw1g/p05OfoHG4kGVSDVPaQEm++oIpBQkICrVq1Ytq0aVy6dMnk451z4ugcvwb3dF2TqWfGBbre+x7PrL+LO1RZk8sHby5l4lWc948VTaYyKRfL7VYo5KtPnz5P3Gd0k9Fs+HsD9zLuoVAomNR8Ek0rN+Vs7FnmHJ9DtjbbYmuJ5cqVY1OgI+XUBvFn3If1j4/PtCSurq5oNBqqVKlCxYoVTT4+xe45Ip4byZ0KTVAqlcQ4+fB/bsO4Yd8g3x6neXudGj5vaT1LDeXk5JCRkSGbGahym0zzDrso6vta2DyxolxMJxKilSlKjyuPCh60r9mezZd0s21UUFfgxoMbXEm8QviNcA5ctcx5XAF4oTfdBoznn3tasjQGtdysNN2PhZoxYwbBwcFMnDgRZ2dnk4+XUHDDqQmSQU9GRRlsqtNoNGRmZpKZmWnuUAyUvXLISy7lIhKilTGc5qgwoxqPYuPFjfo5Tju4d8DNwY1uXt1oVqVZCUZYwpQq8GzF1P/LYPtFgzFNCZd0M9dY8KTf1atXZ/Xq1bzzzjvPfC6P9PN0jl+DR2bZajLNycmRxQdvLtFkqiOXchFTt1mZSpUqFWk/zwqetKvZjo1/byS4STBta7TFzcENraSlernqJRxlCXuhN2MH+vPtD7vp38hg/OHl/bohGA0DzBaaueQOzAfd/eFb5ZsQrhzBTdu62D6sNebt+V3QPKeG20qlQY3TAnqZZmdnk5Ymn5aCgppMC2PqeyvXsjAkl3IRNUQrs3LlyiLvO6rxKDZc3MCDLN00dS9UfIEGlRpY9Iw1AChVvPLJeq4k2XDubt5elPL/cMjPwoULGTlyJKDrcfqsJIVS34SqkDR4lpEep7kfvHL48H1ENJnKpVws/JNPyMuUCck9K3jStkZbNlx8NP9pSlYKv9/+najEKK4mXQUss8epSqVi72A3alc0+BNXqOH6bxZ5L/HKlSv6SS2Ke55dj4wLdE0ILRM9TuXywZtLNJnqyKVcREK0MseOHTNp/1GNR7H+4npSslIASMlO4WrSVf6M+5Pvz3/P3dS7FtHk8pjLB2haMZXwqzkkZz5M6FIWHF8GEXN09xIv7raYe4oKhYL09HQiIyOfOLTGVNftG/C/SkO5pn6hWM8rRzk5OaSnp5Oenm7uUICnazK1RnIpF5EQrYyptQcvZy/a1GjDhr91tcTnnJ7Do7wHv978Fc8KntxJu1MSYZa8uj2gxUi+/zObdX/mSXrVmsLlA7D5Dbh80Czhmerdd99FkiTWrVvH559/bvLxSqUSlUpV4I+tUonS1haVUsHz2RdRKRX6ycBtbW0LHI5RlK7+choGkJ2dLYsPXmO6L2wKhQKFNgf7q//TLV9WhsilXERCtDLdu3c3+Zi8tcS6rnWp41qHoHpBRMZHEp9ugYPalSroMYexEyaz9FT2o2ZfpQNE/Qxe7SFoPdQ1/f0yBw8PD+bOncsXX3yBh4dHsZ7bPS2SjrGr8ci4iGfGBbrdX4tHxsVivYZcyKVpLlfeJlO7a+G4HByP/fWfzRtYKZNLuYiEaGVWrFhh8jG1nGvRunprNv69EYDKjpXxcfNh95Xd1HOth5uDW3GHWTqUKjo4XAIkfr3+sMOIJh3Ofg/f+0PtThazCsaCBQsYMGAA/fr1Y/78+cV67huOjThcZQTX7etzzb4BB12HcN2+frFeQy7k0jSXK2+TaaZnRxJ7LCHDo4N5AytlcikXMezCyjztHLGjGo/irQNvMbD+QJxUTrSq1gqpmoStjWX/iSgSr3JgkCM1KuRpqrvzJ/w8D7qGmCUuU2m1WjZu1H1hMaXjVC6lUolarX7seV3NWcUddXOUGRkoJA22St0QDaXBxAYFDbUoygTgcpLbNAdQvnx5M0ej6/OsQHrUnJzP+ybX97I4yaVcRA3RyuT3oVcU3s7etKrWSl9LVNooLT4ZAvDmbtwrurL+r2xuPzAYZ+dUFVw84fwOi+hYc/78eTZs2MDatWu5ffs2+/btK5Hr5A7Yz53z1NrIZYqwXMrEq1TYPxb1tXAA1NfCcTkwrsw1mcqlXERCtDIbNmx48k4FGN14NOsurCM1O7UYIzKzu3+CJoljNzSsPG2Q+FLvwt63YesQi+hY07FjR/1q4t26dSM+vmTu6153aKib49TBOhcSlssHby6NSy2Sey7TLRAMZHl2JLHH0jLZZCqHchEJ0cpMnDjxqY/1dvHm5Wov62uJVqFuD+g6h7Et1Kw8k0W2Js+YynLVISVO9jVFV1dX9uzZw969e3FxcWHw4MEmHa9UKrGzs3vsR61W639UKhW2antuO7+Irdoeta0NdXL+Rm1rU6QepwX1PpUTuXTe0LNRkl27OwpbtcF7BgqePGl3UXv2WgK5lItIiFYmPDz8mY4PbhzMugvrSMuWyQfGs1Kq4PR/8amqpLarDeHReWZjSbllETXFPXv2sGXLFjZv3syBA6Uz+XruElHW1ONULjWRgqivheO8fyx2osnULKzgJpFQnLxdvGlboy1nY8/SpkYbc4dTPIbsg5Vd2DvwNuXtchOiEtCA2wvQejzYV5D1EIz09HSuXbsGQGpq6TRpX7NvwKGKQ3U9TjPksX7gs8rKyiIlJcXcYRQoy7MjSX7LyKzR1tyhlCq5lItIiFamY8eOz3yOz9p8ZjFNLUVSzg3av035/e+z+EQWXbyVvJA7kiT+b9g9HpyqQdo9cHSFF3rJbjhGSEgIixcvBmD69OkmH29ra4udnR1g3DPUcFo+jca49qzNycI220bXlJrzqENSQb1MLWGi7+zs7FL7QlEU+TdrSgX2OM095knbhV1PjuRSLqLJVHiMXP/TPLXLB2D/FEDLnRQtS0/mc68w9bZsm0737dvHpUuX8PX1xdfXl8uXL5fKdd3TIukU9x3uaedL5XqlISsri9TUVFl8+OanrDaZyqVcRA3xGcTExBASEmLuMEqUl5cXQ4cONXcYz6ZuD/BbAPvfZ3RzNU1XpPJ5ZzvKqQ0Sf7nq0OEDXQ1RZk2ncXFxgO6LilarNaqVlaQbjo04XHk4NxwbQro81g98VnKpiRRE32Ras525QylVcikXkRCfwapVq8wdglAUShW8PAJeHoE70CmmH7+0G4afn5+5IyuSuLg4pkyZwsqVK9m1axcNGzY0+Ry5vUyBAhNq3lVNtFotd+1boAaj3rm5wz8A/RqLec9bWCuDOVsgcmsicqRQKMBWTXbtbigMmq/z60H6pO2iloVcyKVcREIUrF9WGvz0KaQlwrmNbKovoTx+CI7nt7MtVKoDXm3AqzU0eMXs9xNjYmIAOHnyJHv27HmmoTVPSyFp8Mi4wHV7yx6fmDtF2LM6cuQIS5cuZdOmTcUQlQFNNurow6S768Yh2l0LJ7tGW7CR1z3t4lZc5ZIfU8pK3EMUrN/Pc3XLPp3LnYVHwdxfM/k9JiefnXMg4SKcXgXbhsvifmJcXBwHDx6katWqgHENrbTkrpnokWHZM9gUdVWFqKgohg8fTmBgIKDr2TtkyBBGjhzJ+vXrad++PU2bNi32+NTRh6mwfyx21yMMJvo+UuzXkZviKpf8mFJWooYoWL8OUyEnB9LuwzndTD72trDoRBataub9L/CwhlirDXj+Rxb3E6dNm8ahQ4eYOnUq8HQrmqhUKhwcHIBHzWndunUjOTnZaL/y5cuzbds2wLjX6Z0KTQm3HcEdx0bYpj760CqoybSgbXObPXu2fqzbrVu39H0Acjss5fL29mb16tX6D97t27cTGBiIv78/QUFBDBo0qFjiUSgURu9PjncXHvRaQa1uY0hMTNQ9OXMoLi4unD9/Xn+M4fFP2i7s2nJRXOXi5OTEjh07APDz8yMoKMikOERCFKyf2hH8Hq4h+NoyAIYmJjKzVi3ujvlbX/OSqwYNGtCgwaOmyldeeaVYzpucnMyRI49qHxkZGXTr1i3ffSWFkhtOTYrluuaU+6XCVDExMfj4+AC6+7Hnzp3j6NGjHDp0qMD37KkoVWTX7kZiYiKxt2+ivhZOao221HD3Kr5ryFBxlUtAQAABAQFG+5hSVvL56iYUasuWLXz88cc0b97c6ENMeDouLi4MGTKE06dPmzsUwQLUrFlTfy9Xq9Xi4+PD7t27izcZ5pE7BKOsTfRtirzlkh9TykrUEC3E66+/jru7O5mZmbRv397c4ViFb775xtwhlBqVSoWjoyNg3MyZ24xqyN7eHjD+gDFsPs3OfjSOMzPz0XCMp+lxKlcJCQlMmzaNs2fPMmfOHCZOnMj48ePZu3fvUy+xlp+8TaaGcmp1JrnncrIezlqTu5+pS3HJqcn6WZV0uYiEaCFu3rzJ4sWL+f77780dilWZMWMGLVq0oHfv3uYOpdRVqFCBFi1aGD1XlLXoFJIGj/TzXHcwffiHpahUqRLLly83em7NmjWldn0XFxfcqlYzeiyUfLmIhGghJkyYQNWqVQkJCWHQoEHUr2+dK5qXtjp16rBo0aIymRCPHj1qNJlyUbu9e6Sfp0tCKP9X6S3O4VlS4ZVp//77r37bsEYulCyREC3E9u3bzR2CZdNkw4UdcPn/4Nyj8UiBORKTf03h8oRy1K2kfPw4lRMMP6wbioFClvOclrbrDg35v0pvcd2hATww/2BqQSguIiEKZcPlA7BtxGNP29sqeO8/aqLuS9StlM9x2amwuiPkLocVtB7qW15tUqVS4eTkpN/OVdB9v7wM7yfm5OQQ5/QyDkBm9qN7i4b3Ewua6Du/x2VZYfcQc2m1Wt2A/WvhZHl2LNLwlrJwP7EkiIQolA11e8Brq+HyT/qxiLnea6Ob0kySpMc/rPU1xL8BhSzGJQplT26P0yS/ZWRUL1vznJYmkRCFskGpAp9A3c/DsYiGpk6dSp06dRg+fHj+xz/3QgkHKAgFy530O8vTF7I1T9xfeDoiIQoC0KFDB6ZNm8awYcOssknPzs6OChUqALqJlHMZdqoxbOaEgpvXDKeOy85Mxz0tkhuOjUhPV+ufN2yWLayZrqw34dnY2OibrQ2bpR8bU6dSIdXriQrQkGHwdP7N30VtCheMiXdKENBNh5aUlMSJEyfMHYpFcU+LpGPsampa0ZqJQtklEqIgoPsWPXPmTFJSUswdikW54diIw1VGEONovWMShbJDNJkKwkNvvPEGWq2WzMxM/dqBcnDx4kUWLlxIfHw8nTt3ZsyYMQBERkYyZ84cAD788EMaNWpU4Dns7Oz0g+4Nx7UZ/p5qtdromIKaPQ3XTczOziaxQhvKAWkZj5piDde2y9sUa41N0k/L1tb2iTMDFfb+Ge5nWJaG24bNp3Jdp1IuREIUBANTp07Fzc2N999/39yh6NWvX5/ly5ej1WoZPHiwPiEuXLiQJUuWoFAoeP/991mxYoXRcREREURERAAQHR1dylELguURTaaCYKBfv34sX77c6Ju3HOzatYtevXrRs2dP/XNJSUm4uLjg7OzMgwcPHjvG19eXkJAQQkJC8PLyKsVoBcEyiRqiIBho0aIFlSpV4uDBg0bJx9z69OlDnz596NWrFwMHDgTA2dmZpKQkFArFE+cgjYmJYejQoWZJjHnn4SyteTmjo6Px8vLSr4YgR+Ysl7wTu+c30XtJkHO5KCTDGwIWJPebryAUt4MHD1KuXDnatGlj7lAAXdPn9u3byczMpHHjxvz++++sW7eOyMhI5s2bB8D7779f6D1EKHv/Zyzl97WUOIuLnH9f2dUQ9+/fz59//klKSgqfffaZucMRyqDu3buTkZFBUlISzs7O5g7nsVXDx40bB0CjRo1Yu3atmaISBOtTogkxKiqK2bNnk5SURFhYGKmpqYwdOxa1Wo2vry+DBg167Bg/Pz8yMjLw9vYuydAEoVCzZs0iOztbXwOzBoZJtSywlN/XUuIsLnL+fUu0U423tzerV6/WP96+fTuBgYH897//ZdeuXezYsYOhQ4cydOhQNm/erN/vr7/+okmTJiUZmiAU6q233mLNmjVFXhLJEsj5g6gkWMrvaylxFhc5/76l2mQaExODj48PoBtbExAQQEBAgNE+qamp+ln58zLsRn7y5MkC26Fzb9paEkuLuSzEO27cOL744osi7+/l5cXQoUNNC0wQBNko1YRYs2ZNYmJiaNq06eNz9T3k5OTElClT8n0t772Ugsj5pm1BLC1mEa8gCNamRJtMExISCA4O5uzZs8yZM4e+ffuybds2xowZg7+/f0leWhDKvKioKIYPH05gYKDR85GRkQwaNIhBgwYRGRlppuiKV0REBO3atSM4OFjfigQQHh7OkCFDGDRoELdu3TJfgAZEucizXACQrFB4eLi5QzCZpcUs4rUcr732mtHjESNGSPfv35cSExOlUaNGmSmq4hURESH16NFDGjJkiPTPP//on3/99dcljUYjnTt3Tvr000/NGOHjRLnIr1xkN+yiOMj5pm1BLC1mEa/lyp3hBsh3hhtL1K5dOzp06MDdu3eZPHky69evB3TzrtrY2ODp6SnLgeCGRLmYn1UmREEQCmbKDDeWInfycVdXVzIzM42e12q1XL9+nZo1a5orvCIR5WJ+FjtTjSAIhUtISGDatGn873//Y8SIEVy4cOGpZrixBNu3b+fgwYMkJiYyZswYVq9ezbp16zh8+DA//PAD2dnZfPHFF1SvXt3coYpykWm5gEiIgiAIggBYeZPpkSNHOHbsGNWrV2fw4MHmDqdIjhw5wtKlS9m0aZO5QymSvXv3curUKWrXrs0bb7xh7nCe6OTJk4SHh1OlShUxZlAQBCMWnRCfNDVcTEwMU6dONWlwdUkrynR2x44dM3eYekWJ19fXl0WLFpk7VKBo8W7bto0qVaqYO1RBEGTGotdDfNLUcHJcAfpJMcvNk+LVarXMnz+f4OBgM0b5SFHe37lz55KSkmKuEAVBkCmLToh5xcTE4O7uDuimhmvVqhVz586VzQ3b/OSN+dy5cxw9epRDhw6ZObL85Y133rx5JCUlcfz4cTNHlr+88e7atYuQkBBUKpWZIxMEQW4susk0r7xTw7Vv35727dubO6xC5Y3Zx8eH3bt3mzusAuWNd+rUqeYOqVB5481daFcQBCEvi64hWuLUcJYWs4hXEISyQgy7EARBEAQsvIYoCIIgCMVFJERBEARBQCREQRAEQQBEQhQEQRAEQCREQRAEQQCsbByiIAgl5/z588yePRs3Nzdq1qxJlSpVcHNzo3fv3uYOrUwT5VJ8REIUBKFIDh06xJtvvomfnx8AoaGhAPz777+EhIRgb2+Pv78/9+/fJzw8nAYNGqBUKpkyZQrXrl1jwYIFSJJE7dq1efvtt/XnDQkJ4f79+1SqVInp06eb4TezbKJcio9oMpWx0NBQevfuzZAhQ55qgvLAwMBiiUOr1eq3Q0JCiIyMBODAgQOEhoYSGhpK06ZNAcjOzsbT05M9e/YQGhrKnj17iiUGwfyGDx/Or7/+yvDhw1mxYoX++eXLlzNr1ixWrVrFDz/8AEDXrl354IMPOHXqFFqtlqVLl+Lg4EClSpU4d+6c/tibN2+Sk5ODi4sLR48e5ciRI/Tv37/UfzdLJsql+IgaoswFBwfTu3dv/R/j7NmziY+P58GDB3z99dfs3LmTs2fP8uDBA5YsWcJvv/3GypUrqVOnDgAXL15k5syZeHh4MHjwYP2ioxEREXz55Zd07tyZ2NhY5syZw7Jly7h8+TKJiYnMmjWLjz/+GC8vL3x8fHjttdcKjbNevXocO3aMu3fv0rp165J9UwSzqFChArNnzwagV69e+r8JSZL0E+nn/puTkwPoviCB7kvVm2++SePGjY3O+cknn7Bw4ULi4uK4ceMG7du3l9VqL5ZAlEvxEQlR5v773//y+eefM3bsWC5evMiRI0do3bo1GRkZXLx4EQCVSsXNmzc5e/Ysy5YtY82aNSQlJTFu3DhiY2NxdXVl4MCBj63A3apVK9555x0mTZrEjRs3WLt2Ld27dwfgzJkzAIwcOZIaNWo8Mc7AwEC2bdtGWloa3bp1K+Z3QZCDHTt2cPDgQWxtbWnQoAE2NroGptGjR/PJJ5/g6OjIgAEDSEpK4tChQ/z111+0bNkSGxsbxo8fz0cffUS1atUoX748M2bMAKBhw4Z8+eWXJCQk0KxZM3P+ehZLlEvxEQlR5kaOHEmnTp0YPXo0zZo1o2HDhoSEhOhf/+yzz9i1axczZ84kLS0NGxsbbG1tsbOzA6BDhw76JZH+/PNPhgwZoj8277fFGjVqGJ17+/btODs7G8Xj6upKXFwcALGxsVSsWJF79+7h4OAAQNWqVfX/IQXrEhAQQEBAQL6vrVu3Tr8dGhrKwIEDjTp1eHp6sn79+seOe/fdd40eG672Ir5YFY0ol+IjEqIFcHR0pGXLlkRFRWFjY8PkyZNJT0/Xf7ObN28eJ06coEOHDowePZr333+f5557DoDw8HB2797NvXv36Nq1q9F5T5w4wQcffICDgwPu7u60bNmSCRMmIEkSw4YNyzeWAQMG8M4777Bz506SkpJYvnw5GzduBGDevHkoFArWrl1bsm+IIGtDhw596mPlvtqLJRPl8mRicu8yKiIigsjISMaPH2/uUARBEGRBJERBEARBQAy7EARBEARAJERBEARBAERCFARBEARAJERBEARBAERCFARBEARAJERBEARBAOD/AYdPNb2G58ZVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 450x225 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ncols = 4\n",
    "figheight_in = PAGEWIDTH_IN/ncols * 2\n",
    "path='figures/main_dropseq_broken_zeta.pdf'\n",
    "with mpl.rc_context(PAPER_CONTEXT):\n",
    "    \n",
    "    fig = plt.figure(figsize=(PAGEWIDTH_IN,figheight_in))\n",
    "    spec = fig.add_gridspec(2, 4)\n",
    "\n",
    "    ax_big = fig.add_subplot(spec[:2, :2])\n",
    "    ax10 = fig.add_subplot(spec[1, 2])\n",
    "    ax11 = fig.add_subplot(spec[1, 3],sharey=ax10)\n",
    "    ax00 = fig.add_subplot(spec[0, 2],sharex=ax10)\n",
    "    ax01 = fig.add_subplot(spec[0, 3],sharey=ax00,sharex=ax11)\n",
    "    \n",
    "    ax_inset = ax_big.inset_axes([0.2,0.15,0.35,0.35])    \n",
    "\n",
    "    \n",
    "    reads_per_umi,abs_frequency = np.unique(obsdata,return_counts=True)\n",
    "    rel_frequency = abs_frequency/sum(abs_frequency)\n",
    "    ax_big.scatter(reads_per_umi,rel_frequency,\n",
    "                   label='Drop-seq (run A)',c=colormap[protocol_fig4],\n",
    "                   linewidth=0,s=POINTSIZE_HIGHLIGHT,rasterized=True)\n",
    "    ax_big.set_yscale('log')\n",
    "    ax_big.set_xscale('log')\n",
    "    ax_big.set_ylabel('Probability')\n",
    "    ax_big.set_xlabel('Reads per UMI')\n",
    "    ax_big.minorticks_off()\n",
    "    ax_big.set_ylim([10**-7,1])\n",
    "    sns.despine(ax=ax_big)\n",
    "\n",
    "    obs_max = np.max(obsdata)\n",
    "    x=np.arange(len(p))+1\n",
    "    pmf_label='Broken zeta PMF\\n' fr'($a_1={a1_sim}$, $a_2={a2_sim}$, $b={b_sim}$)'\n",
    "    ax_big.plot(x[:obs_max],p[:obs_max],'--',color='k',linewidth=LINEWIDTH,label=pmf_label)\n",
    "    add_largedot_legend(ax_big,loc='best',kwargs=dict(fontsize=LEGEND_FONTSIZE,frameon=False))\n",
    "    plot_inset(ax_inset)\n",
    "\n",
    "    \n",
    "    \n",
    "    ### HEATMAPS ###\n",
    "    all_data = np.array([means,FFs,maxs,alphas])\n",
    "    shared_map_data = np.array([means,FFs,alphas])\n",
    "    logmin_shared = np.log10(np.min(shared_map_data))\n",
    "    logmax_shared = np.log10(np.max(shared_map_data))\n",
    "    \n",
    "    for i,(ax,data,stat,title) in enumerate(zip([ax00,ax01,ax10,ax11],\n",
    "                                           [means,FFs,maxs,alphas],\n",
    "                                           ['mean','ff','max','alpha'],\n",
    "                                           [r'$E_Z$','$FF_Z$',r'$\\max(z)$',r'$\\alpha_Z$'])):\n",
    "        \n",
    "        data=data[:,:,0].T ##transpose to make a1 (first dimension in array) the x-axis dimension\n",
    "        tolerance_level=0.1\n",
    "        obs_stat_lo_tol = observed_amp_stats_df.loc['DropSeqA',stat] * (1-tolerance_level)\n",
    "        obs_stat_hi_tol = observed_amp_stats_df.loc['DropSeqA',stat] * (1+tolerance_level)\n",
    "        \n",
    "        sim_matching_obs_idx = (data >= obs_stat_lo_tol) & (data <= obs_stat_hi_tol)\n",
    "        y_match,x_match = np.where(sim_matching_obs_idx)\n",
    "        \n",
    "        sim_a1_idx = np.argmin(np.abs(a1_grid - a1_sim))\n",
    "        sim_a2_idx = np.argmin(np.abs(a2_grid - a2_sim))\n",
    "        \n",
    "        \n",
    "        if i in [2,3]:\n",
    "            ax.set_xlabel(r'Slope $a_1$',labelpad=10.25)\n",
    "            ax.set_xticks(ticks=np.arange(len(a1_grid))[::len(a1_grid)-1])\n",
    "            ax.set_xticklabels(labels=a1_grid[::len(a1_grid)-1])\n",
    "\n",
    "        if i in [0,2]:\n",
    "            ax.set_ylabel(r'Slope $a_2$',labelpad=-15)\n",
    "            ax.set_yticks(ticks=np.arange(len(a2_grid))[::len(a2_grid)-1])\n",
    "            ax.set_yticklabels(labels=np.round(a2_grid[::len(a2_grid)-1],decimals=1))\n",
    "            \n",
    "        ax.set_title(title)\n",
    "        if stat in ['mean','ff','alpha']:\n",
    "            im = ax.imshow(np.log10(data),origin='lower',cmap='Greys_r',vmin=logmin_shared,vmax=logmax_shared)\n",
    "            ticks=[logmin_shared,logmax_shared]\n",
    "        else:\n",
    "            im = ax.imshow(np.log10(data),origin='lower',cmap='Greys_r')\n",
    "            ticks=[np.log10(np.min(data)),np.log10(np.max(data))]\n",
    "        ax.scatter(x_match,y_match,s=POINTSIZE_HIGHLIGHT,c='tab:orange',linewidth=0)\n",
    "        ax.scatter(sim_a1_idx,sim_a2_idx,marker='s',facecolor='none',edgecolor='k',linewidth=LINEWIDTH,s=20)\n",
    "        \n",
    "        cax = ax.inset_axes([1.05,0.25,0.05,0.5])\n",
    "        cbar = plt.colorbar(im,cax=cax,ticks=ticks)\n",
    "        cbar_y =cbar.ax.get_yaxis()\n",
    "        cbar_ticklocs = cbar_y.get_ticklocs()\n",
    "        cbar_y.set_ticks(cbar_ticklocs)\n",
    "        cbar_y.set_ticklabels([fr'$10^{{{np.round(v,decimals=0).astype(int)}}}$' for v in cbar_ticklocs])\n",
    "        cbar.outline.set_linewidth(SPINEWIDTH)\n",
    "    \n",
    "    #hide y ticklabels in 2nd clm\n",
    "    plt.setp(ax01.get_yticklabels(), visible=False)\n",
    "    ax01.tick_params(axis='y',length=0)\n",
    "    plt.setp(ax11.get_yticklabels(), visible=False)\n",
    "    ax11.tick_params(axis='y',length=0)\n",
    "    \n",
    "    #hide x ticklabels in 1st row\n",
    "    plt.setp(ax00.get_xticklabels(), visible=False)\n",
    "    ax00.tick_params(axis='x',length=0)\n",
    "    plt.setp(ax01.get_xticklabels(), visible=False)\n",
    "    ax01.tick_params(axis='x',length=0)        \n",
    "    \n",
    "    plt.tight_layout()\n",
    "    \n",
    "    #panel letters\n",
    "    letters=['a','b','c','d','e']\n",
    "    x_corrections=[0.065,0,0,0,0]\n",
    "    y_corrections=[0.065,0.175,0.175,0.14,0.14]\n",
    "        \n",
    "    for i,(ax,letter,x_corr,y_corr) in enumerate(zip([ax_big,ax00,ax01,ax10,ax11,cax],letters,x_corrections,y_corrections)):\n",
    "        ax.text(LETTER_LOC_X+x_corr,LETTER_LOC_Y+y_corr,letter,transform=ax.transAxes,fontweight='bold')\n",
    "        if i == 0:\n",
    "            axes = ['bottom','left'],\n",
    "        else:\n",
    "            axes = ['bottom','left', 'top', 'right']\n",
    "        for axis in axes:\n",
    "            ax.spines[axis].set_linewidth(SPINEWIDTH)\n",
    "    \n",
    "    \n",
    "    fig.savefig(path, dpi=1000, format=None, bbox_inches = 'tight', pad_inches = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}