capper_rocs_paired.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import heapq\n",
    "import logging\n",
    "import sys\n",
    "import time\n",
    "import math\n",
    "from pairedCapper import PairedReads\n",
    "import capper\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Got  55640  in  103.59122705459595  seconds\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "# Parse the 1kg reads\n",
    "\n",
    "reads_1kg_small = PairedReads(\"minimizers_1kg_paired/minimizers_correct\",\n",
    "                              \"minimizers_1kg_paired/minimizers_incorrect\",\n",
    "                              \"minimizers_1kg_paired/minimizers_mixed\",\n",
    "                              max_correct_reads=40000)\n",
    "\n",
    "print(\"Got \", len(reads_1kg_small.reads), \" in \", time.time() - start_time, \" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Got  72794  in  152.98134016990662  seconds\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "# Parse the hgsvc reads\n",
    "\n",
    "reads_hgsvc_small = PairedReads(\"minimizers_hgsvc_paired/minimizers_correct\",\n",
    "                                \"minimizers_hgsvc_paired/minimizers_incorrect\",\n",
    "                                \"minimizers_hgsvc_paired/minimizers_mixed\",\n",
    "                                max_correct_reads=40000)\n",
    "\n",
    "print(\"Got \", len(reads_hgsvc_small.reads), \" in \", time.time() - start_time, \" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def current_cap(r):\n",
    "    return r.capped_map_q  # round(min(r.uncapped_map_q, f(r.xian_cap), r.vg_computed_cap, 60))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5QcVZn/8feTmUyCBCGBkAkJ+cEK+I0oiYzooiYqGMIXT1BkISEqIhgEVBbXrwQQxQy4sniWRRchIBkQkCCRPaBBQwBhjARhshkC4UfIRljGdIawIUsEQk9mnu8fVZ3pmXT3dGe6u7qrP69z+kx31e3qp2tm6ql769a95u6IiEhtGhJ1ACIiEh0lARGRGqYkICJSw5QERERqmJKAiEgNq486gEIccMABPmnSpKjDEBGpKqtXr37N3UdnWldVSWDSpEm0tbVFHYaISFUxs5ezrVNzkIhIDVMSEBGpYUoCIiI1TElARKSGKQmIiNQwJQERkQrU2AhmYPsksDNnYPtsxixYXkxKAiIiFaizM3wyoxkmrITpzX2XF0leScDMZpnZC2a2wcwWZFh/jZm1h4/1ZrYtXD7VzFaZ2TozW2tmp6W9Z7KZ/dnMXjSzu8ysoXhfS0QkBkYkYGoLDOmBaS0wYnPRP2LAJGBmdcB1wAnAFGCumU1JL+PuF7r7VHefCvwUuCdc9RbwJXd/HzAL+Dcz2y9cdxVwjbsfCrwOnFWMLyQiEhszmsF6gufWvas2UEz51ASOBja4+0Z3TwJLgJNylJ8L3Ang7uvd/cXw+SbgVWC0mRnwKWBp+J5bgc/u2VcQEYmhVC2gPhm8rk+WpDaQTxIYB7yS9rojXLYbM5sITAYezrDuaKAB+C9gf2Cbu+/MY5vzzazNzNq2bNmSR7giIjGQXgtIKUFtIJ8kYBmWZZuTcg6w1N27+2zAbCxwG3Cmu/cUsk13v9Hdm9y9afTojOMfiYjETv3EP/XWAnYtTFI/+bHifk4eZTqAg9Nejwc2ZSk7Bzg/fYGZvRtYBnzX3R8PF78G7Gdm9WFtINc2RURqy86ddD2yE5LvgXXroKF0/WbyqQk8CRwa9uZpIDjQ39e/kJkdDowEVqUtawD+A/iFu9+dWu7B7PZ/AE4JF50B3LunX0JEJFYuvRSefRauvrqkCQDySALhmfrXgeXAc8Cv3H2dmS00s9lpRecCS8IDfMqpwHTgy2ldSKeG6y4CvmVmGwiuEdxchO8jIlLdOjrgJz+Bc86Bz5a+v4z1PWZXtqamJtd8AiISe2vXwmGHwfDhRdmcma1296ZM63THsIhIJdixA37zm+D5Bz5QtAQwECUBEZFK8M1vwuzZ8PTTZf1YJQERkajdeivcdBNcfDG8//1l/WglARGRKD39NJx7LnziE7BwYdk/XklARCQq77wDp5wC++4Ld94J9fnculVc5f9EEREJDBsGl18O48YVf6KAPCkJiIhEYetWGDUK5s6NNAw1B4mIlNuqVTBxIixfHnUkSgIiImX1xhtw6qlw4IHw4Q9HHY2ag0REyuqRR4KhIZYvh/32G7B4qakmICJSTg8+GPQC+tjHoo4EUBIQESkf96Am8PGPw7veFXU0gJKAiEhpPfFEcA3g9dfBLLgf4Ne/jjqqXZQERERKYdUqOOGE4OLvQw/1jgn0vvfByJHRxpZGSUBEpJh27ICZM+GYY6CtDX70I3jpJZg+PerIMlLvIBGRwXKHF16A9743GAJ6/Hj48Y/ha1+DvfeOOrqclARERPaUe9DbZ+HCoPln/Xo45BBYvDjqyPKm5iARkUK5w+9+FzT5zJwZNPf827/B2LG7ijQ2BteB+z8iGiIoK9UEREQKlUjASScFB/0bboAvfzkYDC5NZ2fmt2ZbHhUlARGRgfT0wL33Bn38r70WDjoIHn4Yjj4aGhqijm5Q1BwkIpJNTw/cfTdMnQonnxw0Ab3+erDuYx/LLwGMSMCXZ8CIzaWNdQ8pCYiIZLJuXTDV46mnQlcX3H47PPts4X38ZzTDhJUwvbk0cQ6SkoCISMrOncFFXoCDD4b994clS+CZZ2DevMJn/hqRgKktMKQHprVUZG1ASUBEpKsLWlqCfv4nnhg0A7373dDaCqedBnV1BW9yzBiCWoD1BAusG6Y3B8sriJKAiNSuZBJuugkOOwy+8pVgrt8f/jDoyzlIa15MMPzvW6A+GSyoT7LXMS20b6is2oCSgIjUrnvugfnzgwlefvvbYJiHk04qShJobm2mx3v6LOv2bpofraxrA+oiKiK14+234ec/D4ZxPussOOUUOOAAOPbYohz4063qWEWyO9lnWbI7yWMdjxX1cwZLSUBE4u+tt2DRIviXf4HNm4OD/1lnBRd6jzuuJB+55pw1Jdlusak5SETi7e67YfJk+Na3YMqU4Iavu++OOqqKoZqAiMTPG29Ad3fQp3///YObvS67rGKmdKwkqgmISHxs2wbNzTBpEvzgB8GyT30qmNRdCSAj1QREpPpt3RqM6XPttfC//wuzZ8MXvhB1VFVBSUBEqt+3vx3c7HXyyfDd78K0aVFHVDXUHCQi1aezE77zHXjqqeD1d78La9cGE7grARRENQERqR6JBFx9dTCG/zvvBNM4HnlkMJuX7BElARGpDpdcAtdcE4zzM29e8Prww6OOquqpOUhEKtemTcFUjgBDhsDpp8Pzz8OttyoBFImSgIhUnr/8JRjTZ9Ik+P3vg2VXXAE33wzveU+kocVNXknAzGaZ2QtmtsHMFmRYf42ZtYeP9Wa2LW3d781sm5n9tt97bjGzv6S9b+rgv46IVLUNG4LRPA89NDjb/+pX4Ygjoo4q1ga8JmBmdcB1wKeBDuBJM7vP3Z9NlXH3C9PKfwNIvzx/NfAu4JwMm/9/7r50D2MXkTjp7g7G8enshPPPD3r/jBsXdVSxl09N4Ghgg7tvdPcksAQ4KUf5ucCdqRfu/hCwfVBRikg8PfssfOMbwbj+dXVw222wcWNw05cSQFnkkwTGAa+kve4Il+3GzCYCk4GH8/z8K81sbdicNCzLNuebWZuZtW3ZsiXPzYpIRVu7Npi794gjgpu82tuD5R//OIwdG21sNSafJJBpkG3PUnYOsNTdu/PY7sXAe4EPAaOAizIVcvcb3b3J3ZtGjx6dx2ZFpGJ1dAR39R55ZHDB95JLgjl9jz466shqVj73CXQAB6e9Hg9sylJ2DnB+Ph/s7onw6Ttm1gJ8O5/3iUgV6+oKJnD/3vfgggtg1KioI6p5+SSBJ4FDzWwy8FeCA/3p/QuZ2eHASGBVPh9sZmPdPWFmBnwWeCbvqEWkOk2eHLT3S8UYMAm4+04z+zqwHKgDFrv7OjNbCLS5+31h0bnAEnfv01RkZn8kaPYZYWYdwFnuvhy4w8xGEzQ3tQNfK9q3EpHK4g533QV77RVM5ThiRNQRScj6HbMrWlNTk7e1tUUdhojk66mn4MorYeXKYNyfadPg0Udhn32ijqymmNlqd2/KtE5jB4lIaWzaBCecEHT//PSng3sAzjgjmNdXKoZ+GyJSXFu2BM0/XV3BNYBFi3TXbwXT2EEiUjxvvQXvex+sWAETJwbNQEoAFU1JQESKZ+PGoCawLRw+zDLdZiSVRElARIrjzTeD8X7q6uAzn4k6GsmTrgmIyOB1d0NTE7zwAlx/fdAUJFVBSUBEBs8dtm+HZcuCHkFSNZQERGRwurpg6NBgXCCpOromICJ77vnnYcKE4AYwqUpKAiJVprEx6HTT/9HYGEEwV14Jb7wBU6ZE8OFSDEoCIlWms7Ow5SWzYQP88pdw7rmgYd6rlpKAiBTOPZgLYOhQ+LZGga9mSgIiVSKxPUHDOTNgxOaoQ4FVq+Duu+Gf/imidigpFiUBkSrQ2AgHzW2mq3ElTG+OOhw45hh4/HG44oqoI5FBUhIQqUCNjWD7JLAzZ2CNT9H5mY/A1BYY0gPTWqKrDWzeDA8+GDQHffjDGhYiBpQERCpIYnuC+vM/SOeZ74J/nAAT/wgnz4Pxf4a6ZFDIujPWBsaMKUOAN9wAM2cG8wJLLCgJiFSAxPYEM26ZweT5F9M9eg0MexvqdoI5HLgumH9vSE9QuD6ZsTawudSVg64uuPFGmDUrGCJaYkFJQCRiie0JjrrxKFpfbuWdw34RHPBTj2yy1AZK6t57g9nBzjuvvJ8rJaUkIBKxQ+YvILE9EbywDNO9ZkoG9Uk4+LFdL8vSFHTddcHAcBobKFaUBEQikGr+eWrzU+w49I7eA/1A11l3NsAT58HlDovWMGZMcI225E1B27bB+vXBjWF1dSX+MCknDSAnEoHm1mZW/vdKTr37VBjSnf8b02oAQ4aU4eCfst9+wcXgrq4yfaCUi5KASJkltidoaW+hx3tY/z/rBz77B+gazpifv9170L+hlBH2k0wGGWfo0OAhsaIkIFJmza3N9HhP7kKJqbBoza6XY8aU8ay/v5YWaG6GtjbdHRxDSgIiZdSeaGfR6kW9SSBTLSAtAXiG68Rl5Q4/+xkceGCZrj5LuSkJiJRR0z9/gZ5RPX0P/jsb4D/Phvuv61O2Io65f/oTrF0LN92ku4NjSr2DREosvSdQ98hndz/779fdMyWy5p90P/sZ7LsvzJ0bdSRSIqoJiJRYqifQvHvmQc9QGJLMevZfUV59FZYuDW4O23vvqKORElFNQKSE0nsCrduyLjjrh6xDP6RURFPQAQfAb34DF1wQdSRSQkoCIiWUsydQv6Ef3HsfFdEUNGQIHH+8xgmKOSUBkRJJ1QKS3cnMBbJcC6gIDzwA3/kObN8edSRSYromIFIiGWsBldwTKN2//is88wz88IdRRyIlppqASIms6li1ey2gPsnU2Y/1afqpmOaflA0bYPlyOOccqNd5YtzpNyxSImvOWTNwoUp0/fXBwf/ss6OORMpANQER6fXWW8EwESefDGPHRh2NlIFqAiLS64034MQTYf78qCORMlESEJFejY1w221RRyFlpOYgEQmsWAH//u+aM6DGKAmISOCXv4Qf/xi6C5jkRqpeXknAzGaZ2QtmtsHMFmRYf42ZtYeP9Wa2LW3d781sm5n9tt97JpvZn83sRTO7y8waBv91RGSP3XJL0F91+PCoI5EyGjAJmFkdcB1wAjAFmGtmU9LLuPuF7j7V3acCPwXuSVt9NfDFDJu+CrjG3Q8FXgfO2rOvICKD9tZbwc9Ro6KNQ8oun5rA0cAGd9/o7klgCXBSjvJzgTtTL9z9IaDPvedmZsCngKXholuBzxYQt4gUS08PfPKTwfNLL402Fim7fHoHjQNeSXvdAXw4U0EzmwhMBh4eYJv7A9vcfWfaNsdl2eZ8YD7AhAkT8ghXRAqSTML3vgd1dTBzZtTRSJnlkwQyTSeUbdK7OcBSdx/oylLe23T3G4EbAZqamqKebE8kXtxh2LDg3gCpSfk0B3UAB6e9Hg9sylJ2DmlNQTm8BuxnZqkklGubIlIqt9wCn/ucRgutYfkkgSeBQ8PePA0EB/r7+hcys8OBkcCqgTbo7g78ATglXHQGcG++QYtUksbGYPrd/o/GxqgjG0B7e3Bn8Natmjmshg2YBMJ2+68Dy4HngF+5+zozW2hms9OKzgWWhAf4Xczsj8DdwLFm1mFmx4erLgK+ZWYbCK4R3Dz4ryNSXontCTpPmJFxhrDOzggCytdbb8HppwfzB992WzCBjNSkvIaNcPf7gfv7Lftev9eXZ3nvx7Ms30jQ80ik6rQn2vngTz+Br58J71sZzBBWyfMF93fJJfDcc8FdwhMnRh2NREjpXySLbM08ZjCt+Qt4w//CEUthSE/O+YIrUmcnzJ4Nxx0XdSQSMSUBkX4S2xM0nDODzjfDg/qIBHw5rclnTDscuC7s4xa2fvabL7ji3XILXH551FFIBVASEAFGH5LAzpyB7bOZg+Y209W4svegPqMZJqS9PvkLvW9MdXauT1ZHbWDnTnj11aBb6LRpUUcjFUBJQGpSYnuCGbfMYPPfgoP2a1PCA/2xC2BqS28Tz5in+r6e/GBaLaCfaqgN/P73MG4cPP541JFIhVASkFhJP7hnbNPfJ2jqmTz/Ylr/spKxpzVj+yR6D/QfuB0snBzeuuHkeX1f/8Op2T+8PgkHP7brZcVNHg9w883B+EBHHRV1JFIhNKmMVLXGxn5dMU9shqNWcsiZzbzdmaG3zoxmuhr/CI1/6j27H/pm74F+SHdwsIfgoJ5+1l+fhLpk5lpA13C48m0gOPhX1MTxKZ2d8Nvfwj/+IwwdGnU0UiFUE5Cq0b8JB/olgBG9Z/RvH56hfX7Xeu890NtOOPL24AAPmQ/w6bob4Inz4HKHy50x1zv+fceveBv3YBSGikwAENwPsHMnnKUBe6WXkoBUnGxdMw+a29zbhBMu62NGc9+mm/7t8zOa0w7+4bL6rt5lmfT/jPokU2c/VvkH/ExaWuCYY+C97406Eqkgag6SitP5ZgK+PAeW3gV/C8deSDvLZ1oLtF7Wuy59feqMPtVbJ1Vu1/oMUydmOvtPTIVFayq3aWdPLFsWDBEhkkY1ASmrTE06u+nfJTO1bMCz/J6+y9LLZVrfJ7Cpu5p4/PuO37Cm+s70BzJpEnzwg1FHIRVGSUDKqrm1mZX/vZLmR3fvStnYSN+eOql+99nO8tPb/Mev6l2fkt5bJ9N66D34L1oDVGiPnsHavh0+/3lYvTrqSKQCWb/x3ipaU1OTt7W1RR2GFCixPcGcX8/hJ7N+wkdu/gg7du5gr/q92HjBRhpH9DbpmAEnngfTbg4O2Dsb4D/PBvPeZSmpdVnG64lVM85g3XwznH02/OlPwTUBqTlmttrdmzKtU01ASq65tZnWl1YydeE8duwImmTe3tHN2NOa+w63nO2Mf0Jr7rN82HWhtiov2Jba4sXBxeC///uoI5EKpAvDMmipM/27Trmrz5l9al1Le0vQHt+/z/20FjpbLwPC92Rr1395Blz/DEOGBNPh9hfLJpxiee45eOwxuPrqDN2pRFQTkCLI1s7f2Bh060yd/e+m/wXeHO367tDdvfsZv876B7B4MdTXw5e+FHUkUqFUE5CCpZ/5uzst7S30eA8t7S1cNuOyXbWBzjf7Ne9k6HOf3qSTujib0Q3F/Q41Y/JkOP98OPDAqCORCqUkIAVLP/N3nB4PzvS7vZvmR5u57sTwYm2m5p1MF3TDA/yYMZln41JzzyCcd17UEUiFU3OQFCTVxt/jPSxes5iWNS0ku4Mz/WR3kpb2lt57AAbqttnP5s1q7imaHTvgqqugK8PNcSJpVBOQgjS3Nu8680927z6YWp/aQK7mnZDO8kvkvvtgwQI46CD44hejjkYqmJKA5C1VC0id+ffQs2tirZRkd5LHOjKf6aerottTqtOtt8I++8DcuVFHIhVOSUDyll4LSGmoa+DsaWf3XgdIozb+iJxxBtx/P1xxRdAzSCQH/YVI3lZ1rNpVC0jJdeavtvwIPPMM/OIXwfOvfjXaWKQqKAlI3tacM3Abv0TsiCPg5Zdh//1h772jjkaqgJKASNxMmBB1BFJF1EVUJC7++lf4yEdg1aqoI5EqoiQgEhf33gt//jPst1/UkUgVURIQiYOenqBb6GGHafpIKYiuCYjEwR//CE88AYsWabRQKYhqAiJx8OCDMGQInHZa1JFIlVESEImDD30ILroI9t036kikyqg5SEqmsTH7HcO6kazIZs8OHiIFUk1ACpbYnmDGLTN6RwvNIlMCyLVc9tArrwQ3iInsASUBKdiCBxfQ+nIrCx5cQGNjcB2y/6OxceDtSJFcey0cfngwfLRIgZQEJO8z+1TZO56+A4Db195O55uZ36Oz/TJ66CE45hgYPjzqSKQKKQlI1jmCM1nw4AK6vRsI5g7g2AWlDk9y+Z//gfZ2OPbYqCORKqUkUOMS2xMsXrM4mCmsfXHO2sDoQxL8ov2OvguPvB1G6CpvZB55JPj5qU9FGoZULyWBGpCruae5tZmunmAKwmR3kuZHm7O287/2gQVg3X03YKoNROqhh4LJYz70oagjkSqVVxIws1lm9oKZbTCz3f7jzewaM2sPH+vNbFvaujPM7MXwcUba8kfCbabed2BxvlJtKKQdP1tzT3otANhVG8jWzs9hy3abThILl2eQbfIYTSpTRJddFowZpMljZA8NmATMrA64DjgBmALMNbMp6WXc/UJ3n+ruU4GfAveE7x0FfB/4MHA08H0zG5n21nmp97n7q0X5RjUivYdOLukTw/eZBJ6+tYCUZHcSpme5NvDGwVmWj8+4WBPHl8HYsfDJT0YdhVSxfGoCRwMb3H2juyeBJcBJOcrPBe4Mnx8PrHD3re7+OrACmDWYgOOgkLP4bO9P76GTazvpU0KmJoFPuXF5627TRfZ4D0x8NPPGFq2By333Rx4TyksJrFgBn/88bN0adSRSxfJJAuOAV9Jed4TLdmNmE4HJwMN5vrclbAq6zKx2Rr3K9yw+1/vTe+hk207/ieGT3ck+tYHujdNhZ0PfN+1sgJdn7FFc6dTkUwZr18I990B398BlRbLIJwlkOjh7lrJzgKXunvqrzPXeee7+fuDj4eOLGT/cbL6ZtZlZ25YtW/IIt7IVchY/0PtTsm0n08TwfWoD41dBfd85g6lPwsGZ5wzOZswYNflE4pFH4O/+DkaPjjoSqWL5JIEOIL0xeDywKUvZOfQ2BeV8r7v/Nfy5HfglQbPTbtz9Rndvcvem0TH4Y8/3LD6f96dk286AE8PvQfOO2vgrRDIJf/gDzJwZdSRS5fLpUvAkcKiZTQb+SnCgP71/ITM7HBgJpM9ttxz4YdrF4JnAxWZWD+zn7q+Z2VDgM8CDe/41qkO2s/gfHfcjGkfkN87Cshcz98RZtn735ZoYPsY2bYI334Sjjoo6EqlyA9YE3H0n8HWCA/pzwK/cfZ2ZLTSz9GEL5wJL3N3T3rsVaCZIJE8CC8Nlw4DlZrYWaCdILjcV6TvtZrAXYoulkLP4bA7eN3MPnfH77t5Dp9jj+qidv4LU18NHP6pfigyapR2zK15TU5O3tbUV/L7zlp3HotWL+NpRX+O6E68rQWT5GXXVKF7f8fpuy0cOH8nWi4rfwyPXpXb3gdeLSDyY2Wp3b8q0LvZ3DOfqJy9SlZJJeO21qKOQmIh9EsjVT16kKi1dGrTpPfFE1JFIDMQ6CQzUT77chljm3V1ndWWORKpWVxdceCF84AMwbVrU0UgMxDoJDNhPvswKuahbDhrbpwo98AC8+ir84AcwdGjU0UgMxHrUqQH7yZdZpXXZVP/+KnT77TBqFBx/fNSRSEzEOgmsOWcNie0JJl87mXe632F4/XD+csFf8u6TL1JR3nkHli2DefOgoWHg8iJ5iHVzEGQeL1+kKg0bBi++CJdeGnUkEiOxTgLZxstXN1GpWmPGwPhoriFJPMU6CWQbL1+1Aak6mzfDccfBHtwsKZJLrJNA68uZx8t/9OUs4+WLVKq77gqmknzXu6KORGIm1klg+sTpNNT1vYDWUNfAjImDHy9fpKzuuCO4L2DKlIHLihQg1kmg0rqIlpvuA4iJ9evhySeDXkEiRRb7LqK1TPcBxMQddwSj/c2dG3UkEkOxrgmIxMJ73gPnnw8HHRR1JBJDsa4JiMTCF78YPERKQDUBkUr25JPwxhtRRyExVhNJoFJmFiuGurrMs4XVaSDS+Pnb3+DUU2HWrKgjkRiriSSw4MEFtL7cWvCk7pWop6ew5VLFFi6El16CBdX/dyuVK/ZJIH1y99vX3h6L2oDUgO5uuP9+OPJImD174PIieyj2SSB9cvdCJ3UXicwdd8C6daoFSMnFOgmk1wJSVBuQqrBtG8yYEVwTECmhWCeB9FpAimoDUhW++U34wx9gSKz/RaUCxPovbNmLyzIvX595uUhFMYs6AqkBsU4C2WYQGzNCg+eIiEDMk4BGERURyS3WSaDWRxEVERlIrJPA/affz7C6YX2WDa8fzu/m/S6iiEREKkusk4CmlxQRyS3WSUDTS4qI5BbrJDB94nSGWN+vOMSG6MKwiEgo1klANQERkdxinQSmT5zOkH5fcQiqCYiIpMQ6CbS+3EoP/WoCqCYgIpIS6yRw1NijMi5vOqipzJGIiFSmWCcBjR0kIpJbrJOAxg4SEckt1klAXURFRHKLdRJQF1ERkdzySgJmNsvMXjCzDWa224wsZnaNmbWHj/Vmti1t3Rlm9mL4OCNt+VFm9nS4zZ+YFX/wdNUERERyqx+ogJnVAdcBnwY6gCfN7D53fzZVxt0vTCv/DWBa+HwU8H2gCXBgdfje14HrgfnA48D9wCygqCO7qSYgIpJbPjWBo4EN7r7R3ZPAEuCkHOXnAneGz48HVrj71vDAvwKYZWZjgXe7+yp3d+AXwGf3+FtkofkERERyyycJjANeSXvdES7bjZlNBCYDDw/w3nHh83y2Od/M2sysbcuWLXmE20vzCYiI5JZPEsjUVu9Zys4Blrrvmt0923vz3qa73+juTe7eNHr06AGDTbfmnDWc23TurtpAQ10D5zWdx5pz1hS0nWJpbAymje3/aMzck1VEpOTySQIdwMFpr8cDm7KUnUNvU1Cu93aEz/PZ5h5LbE/Q0t6yqzaQ7E7S0t7C5r9tLvZH5aWzs7DlIiKllk8SeBI41Mwmm1kDwYH+vv6FzOxwYCSwKm3xcmCmmY00s5HATGC5uyeA7Wb2kbBX0JeAewf5XXbT3Nq824Xhbu+u6kllxmS5zy3bchGRXAbsHeTuO83s6wQH9DpgsbuvM7OFQJu7pxLCXGBJeKE39d6tZtZMkEgAFrr71vD5ucAtwF4EvYKKPudjHK8JbI6mEiMiMWVpx+yK19TU5G1tbXmXT2xPcMhPDmHHzh27lu1VvxcbL9iYdUiJUsp1J0QV/RpEpMqY2Wp3zzhyZqzvGI5jc5CISDHFOgnEsTlIRKSYYp0E7j/9fobVDeuzbHj9cH43r+iXH0REqlKsk0BzazNdPV19liW7k2oOEhEJxToJaOwgEZHcYp0ENHaQiEhusU4ClXZhWDd6iUilGfBmsWoW1RhB2ehGLxGpNLGuCYiISG5KAiIiNUxJQESkhikJiIjUMCUBEZEaVlWjiJrZFuDlEmz6AOC1ErhOfAsAAARnSURBVGx3sCoxrkqMCRRXISoxJlBchSg0ponunnFqxqpKAqViZm3ZhlmNUiXGVYkxgeIqRCXGBIqrEMWMSc1BIiI1TElARKSGKQkEbow6gCwqMa5KjAkUVyEqMSZQXIUoWky6JiAiUsNUExARqWFKAiIiNaymk4CZfdvM3MwOyLK+28zaw8d9FRTXGWb2Yvg4o8SxNJvZ2nAfPGBmB2UpV9Z9VUBcZdtX4eddbWbPh7H9h5ntl6XcS2b2dBh/W4XENMvMXjCzDWa2oJQxhZ/3D2a2zsx6zCxrd8dy7qsC4yrb/jKzUWa2Ivw7XmFmI7OUK/z/0N1r8gEcDCwnuPnsgCxl/lZpcQGjgI3hz5Hh85EljOfdac+/CdxQCfsqn7jKva/Cz5wJ1IfPrwKuylLupWx/d1HEBNQB/wUcAjQATwFTShzX/wEOBx4BmnKUK9u+yjeucu8v4F+ABeHzBTn+rgr+P6zlmsA1wHeASrsyPlBcxwMr3H2ru78OrABmlSoYd38j7eXeOeIqqzzjKuu+CuN6wN13hi8fB8aX8vPykWdMRwMb3H2juyeBJcBJJY7rOXd/oZSfsSfyjKvc++sk4Nbw+a3AZ4u14ZpMAmY2G/iruz81QNHhZtZmZo+bWdF2+iDjGge8kva6I1xWyriuNLNXgHnA97IUK+u+yjOusu+rfr4C/C7LOgceMLPVZja/AmKKel/lEtW+yqXc+2uMuycAwp8HZilX8P9hbGcWM7MHgcYMqy4FLiGoIg9kgrtvMrNDgIfN7Gl3/6+I47IMywZ1dp4rJne/190vBS41s4uBrwPfz1C2rPsqz7iKvq/yiSsscymwE7gjy2Y+Gu6vA4EVZva8u7dGGFNk+yoPRd1XRYqrrP+HBWym4P/D2CYBdz8u03Izez8wGXjKzCCoGv+nmR3t7n0mgHT3TeHPjWb2CDCNoB0wyrg6gE+kvR5P0HZZ9Jgy+CWwjAxJoJz7qoC4ir6v8okrvAD9GeBYDxtqM2wjtb9eNbP/IGhe2OMDWxFi6iC4HpUyHti0p/HkG1ee2yjqvipSXEXfX7liMrNOMxvr7gkzGwu8mmUbBf8f1lxzkLs/7e4Huvskd59E8Mv8YP8EYGYjzWxY+PwA4KPAs1HHRXDReGYY30iCmsPyUsVlZoemvZwNPJ+hTFn3Vb5xUeZ9FcY1C7gImO3ub2Ups7eZ7ZN6Hsb1TJQxAU8Ch5rZZDNrAOYAZesRl02591UByr2/7gNSvdvOAHarrezx/2GprmZXy4O0ngdAE/Dz8PkxwNMEV/2fBs6qhLjC118BNoSPM0scx68J/unWAr8BxlXCvsonrnLvq/DzNhC0FbeHjxvC5QcB94fPDwn31VPAOoImiEhjCl//X2A9wZljSWMKP+9zBCc77wCdwPKo91W+cZV7fwH7Aw8BL4Y/R4XLB/1/qGEjRERqWM01B4mISC8lARGRGqYkICJSw5QERERqmJKAiEgNUxIQEalhSgIiIjXs/wPliFfp4nLTNwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5xcdX3/8ddnd5MQ2CAEkl3IXY0IFdiYVSgU1sotFAwI2gajJhFESWls+/CylFr7YK0WC0WFtHKNRECiUCQRMFwsCTQR2ZgLBgSScEnIZFkCYuAX2Nvn98c5szk7e2Z3NtmdM5f38/GYx8z5zpmznxnC9/P9fs8536+5OyIiUn4qkg5ARESSoQQgIlKmlABERMqUEoCISJlSAhARKVNVSQcwEIceeqhPnjw56TBERIrKmjVrXnP3MZnlRZUAJk+eTHNzc9JhiIgUFTN7Ka5cQ0AiImVKCUBEpEwpAYiIlCklABGRMqUEICJSppQAREQKUG0tmIGNSmHzGrBROzALygeLEoCISAFqaQlfNDTBxMfh5Kae5YNACUBEpFBVp6BuEVR0wbRFUL1jUA+vBCAiUqgamsC6gtfW2d0LGCxKACIihSjd+q9qC7ar2ga9F6AEICJSSLq64Kaberb+0wa5F6AEICJSKHbvhgsugC9+kapJ/7en9Z9W1UbVlFWD9ueKajI4EZGS1dIC55wDv/0tfO97tH/1q8F1oENICUBEJGlPPw1nnRUkgbvvhk9+Mi9/VglARCRp7e0wbBisXAn19Xn7s0oAIiJJaWuD4cPh2GODXkBVfqtknQQWEUlCVxecfjpcdlmwnefKH5QARESScf31sGIFHHFEYiEoAYiI5NvLL8PXvw6nnQZz5iQWhhKAiEg+ucOXvxwMAV1//ZBf6tkXJQARkXx6/vlg6Oc734EpUxINRVcBiYjk0wc+EFzxM3580pGoByAikjePPRYMAU2aBJWVSUejBCAikhfLlsHJJ8MddyQdSTclABGRofbmm3DJJXD00fDpTycdTTedAxARGWrf+AakUnDPPcGdvwVCPQARkaG0YkVwuec//AN85CNJR9ODEoCIyFB65x34i7+AK65IOpJeNAQkIjKUzjgjmPMnwRu+slEPQERkKDQ3w3e/G0z1XICVPygBiIgMvvZ2uPBCuO46ePvtpKPJSkNAIiKD7aabYMMG+MUv4KCDko4mq5x6AGY2w8yeNbNNZtYY8/5cM2s1s3Xh46LIe52R8qWR8ilm9oSZPW9mS8yscK6NEhHZW52dcPXVcNxxMHNm0tH0qd8EYGaVwELgTOAo4AIzOypm1yXuXhc+boqU746UR3+NK4Fr3H0q8AZw4d5/DRGRAnHvvbB5M3ztawU79p+WSw/go8Amd9/i7m3AncA5+/JHzcyAjwN3hUW3AufuyzFFRArCiSfCxz4G5xZ+lZZLAhgHbI1sbwvLMp1vZhvM7C4zmxAp38/Mms3sN2aW/kUOAf7o7h39HBMzuzj8fHNra2sO4YqI5FkqBY2NwfBPTQ38+tcFMdlbf3JJAHF9GM/YXgZMdvdjgIcJWvRpE929HvgM8H0ze1+OxwwK3W9w93p3rx8zZkwO4YqI5NEvfwnHHAM//CGsWxeUFfjQT1ouCWAbEG3Rjwe2R3dw953u/m64eSMwPfLe9vB5C/AoMA14DTjIzNJXIfU6pohIQdu9Gy69FD7xCRg3DtasgenT+/9cAcklATwJTA2v2hkOzAKWRncws8MimzOBZ8Lyg81sRPj6UOBE4Gl3d+B/gU+Fn5kD3LsvX0REJK9mzYKFC4M5fp54Ao48MumIBqzf+wDcvcPMLgWWA5XALe6+0cyuAJrdfSmwwMxmAh3A68Dc8ONHAtebWRdBsvl3d386fO8bwJ1m9m1gLXDzIH4vEZHB5w4dHTBsGFx+OcyfH0z1UKQsaIwXh/r6em9ubk46DBEpR6++CvPmwfvfDz/4QdLRDIiZrQnPxfagqSBERPrzq18FJ3ofeQSmTk06mkGjBCAiks077wRj/GeeCWPGwJNPBid+S4QSgIhINlu3wo03wt/9Hfz2t8GSjiVECUBEJModli8PnqdOpXb/N7Frf4jtPxIzuh+1tUkHuu+UAERE0l57LZjCYcaMIAkALa3xd/S2tOQzsKGhBCAiAvDQQ8EQz69+BddcE6zilVadgrkNUL0jufiGgBKAiMi3vx1U+KNHB2P9f//3UBGpHhuaYOLjcHJTcjEOASUAEZGjj4ZLLgmu8jn22J7vVaegbhFUdMG0RSXVC1ACEJHy4w433BBM4AZwzjnwX/8F++/fe9+GJrCu4LV1llQvQAlARMrLzp1w/vnwpS/tudoni9SuVNDqr2oLCqraunsBNTV5incIKQGISPn49a+DIZ5f/hKuugqWLetz6uamlU0MH9HVo2z4fp3Mv7OJHSUwEqRF4UWkPLz8cjBx2/veB0uXwoc/3O9HVm9bTVtnW4+yts42Vm1bNVRR5pUSgIiUttdfD67umTgR7r4bTjkFDjggp4+u/dLaIQ4uWRoCEpHS5A633AKTJ8PDDwdlM2fmXPmXAyUAESk9b7wBf/M3cOGFUF8PH/xg0hEVJCUAESktK1cGJ3rvuQf+/d+DO3zHj086qoKkcwAiUlo2bIARI2DVKvjIR5KOpqCpByAixW/z5u7J2/jbv4V161T550AJQESKlzssXgx1dcGNXe3twXX9OtGbEyUAESlOf/wjfOYzMGcOTJsWjP0PG5Z0VEVF5wBEpPi8/npwI9e2bcFMno2NUBk/b79kpwQgIsVn9Gj43OfgrLPg+OOTjqZoaQhIRIpDa2twF+/69cF2U5Mq/32kBCAiha+rC2bPDi7t3Lo16WhKhoaARKSwtbfD/PnBDV3XXw9nn510RCVDCUBECtfmzfDXfw2/+x380z/BF7+YdEQlRQlARArXpEnwiU8Ea/R+7nNJR1NylABEpHC4wyOPwE9+Ap/6VFD5/+u/Jh1VyVICEJHC4B7c1PWTn8Dw4XDwwcF4fx8rdsm+0VVAIlIYbrghqPy//nV48034/vdV+Q8x9QBEJHk7d8JXvwozZsB3vwsVapvmQ06/spnNMLNnzWyTmTXGvD/XzFrNbF34uCjj/QPN7BUzuy5S9mh4zPRnxu771xGRovLCC9DRAYccEqza9bOfqfLPo35/aTOrBBYCZwJHAReY2VExuy5x97rwcVPGe03AipjPzI585tWBBi8iRezHPw5W6rrllmD7uONg1KhEQyo3uaTajwKb3H2Lu7cBdwLn5PoHzGw6UAM8uHchikhJ6ewMhnvmzYOTTgrm85FE5JIAxgHRe6+3hWWZzjezDWZ2l5lNADCzCuBq4GtZjr0oHP75pln82R4zu9jMms2subW1NYdwRaRgvfMOnHceXH01XHopPPAAjIurTiQfckkAcRWzZ2wvAya7+zHAw8CtYfl84H53j5u8Y7a7Hw2cFD5i7/Jw9xvcvd7d68eMGZNDuCJSsB59NFi569prg4fm709ULlcBbQMmRLbHA9ujO7j7zsjmjcCV4es/B04ys/lANTDczN5y90Z3fyX87C4zu4NgqGnx3n0NESkKM2YEk7mpMVcQcukBPAlMNbMpZjYcmAUsje5gZodFNmcCzwC4+2x3n+juk4GvAovdvdHMqszs0PCzw4Czgd/v87cRkcK0eTP8z/8Er1X5F4x+ewDu3mFmlwLLgUrgFnffaGZXAM3uvhRYYGYzgQ7gdWBuP4cdASwPK/9KgmGjG/f+a4hIwWpvD5ZufO45+Mu/DO7wlYJg7pnD+YWrvr7em5ubkw5DRAbi8svhO9+Bn/88mN9H8s7M1rh7fWa57rgQKRG1tcHMCZmP2toEg3rxxeDO3i98QZV/AVICECkBtbXQ0hL/XrbyvFiwIHi+7LIEg5BslABEili61d/ydgrmNkD1jqRD6umLX4Rbb4X3vz/pSCSGJoMTSUhfrfYBqU7BxdODyv/kJrh/4SAcdJB84hNJRyB9UA9AJAGDUvlXh63+M74Co1JQ4TBtUWH0Au6+O5jb/803k45E+qAEIJJHqV0pGn7cQAvrew/ZpCv0mvU9n6t39H6vegc0NMHEx+BDP99zv751Br2ApN18M6xYocndCpwSgEie1NbC4Rc0sfKFx+G82TDx8Z6VdUNTUJZ+L7pP5nunNELdoqDVH1XV1qsXUFOTpy+Y1tICDz4YXPuvqZ0Lmv7riORBaleKlrOPDyvtLhi7MXhOV9bVqd7vde9zS+/3jrktaO1D79m6Ir2AmhrYke8RoSVLghk/P/vZPP9hGSidBBYZIj3G+c9qgvonwDPaXOnK2hysK/5AlW29yyo69ySATFVtVE1ZRXtS93jefjvU1cFRccuGSCFRAhAZIt2Vf3UK6m4JWurpSj7daq9qC1r46dfR99LPFZHEkPlejLraOtZ+a+2+hL73OjrgxBNV+RcJJQCRIZDalYK5s+CuJcH4fWV79p3jWvh7oa62jrVfSqjiT6uqgv/8z2RjkJwpAYgMgaaVTXtO1n5oSc9WfKa+3utHQVT6ae7w2GNBD6CyMuloJAdKACKDLLUrxaJ1i/acrPWMyrBjOPzuou4bthI5UTsUfvMbaGiAn/xEJ4CLhK4CEhkk6WkZDr+giXfeCVv1FZ17xvbTqtqom7kK96DRXBKVP8Btt8F++8HMmUlHIjlSD0BkL/W6m7c6BRd+EmrX9zyh2z4SfrAF3gqm5aypgbWlUumntbcHl3/OnAkHHph0NJIj9QBE9lKvqRwammD8E71P6lon8+9sKr0Wf9Ty5bBzp4Z+iowSgMjeqk71nKqhblHQ4s88qVvVxqptqxIJMW9+8Qs45BA444ykI5EB0BCQyN5KT8+QeSNXxklegLXFs/De3vnRj+Af/xGGD086EhkA9QBE9kJqV2TqhvRUDelx/5j5eEpeVZVu/ipCSgAie6FpZdOeFn9lW+y4f3Q+npI2bx58//tJRyF7QQlAZIC6r/NPt/grumLH/dOXepbkSd+07dth8eLgBLAUHZ0DEBmgppVNdHnPCn945XAumnYRC88qoNW48uHOO6GrC2bPTjoS2QvqAYgM0Optq2nr7Dnk09ZZBlf6xLntNqivhw9+MOlIZC+oByAyQAUz907Snn4a1q6Fa65JOhLZS+oBiMjeMYNZs4KHFCX1AERk7xx5JPz0p0lHIftAPQARGbhdu+CZZ6BtcNYykGQoAYjIwD34YHDj1/r1SUci+0AJQEQG7r77gudJk5KNQ/aJEoCIDMzNN8OiRXDJJTB2bNLRyD7IKQGY2Qwze9bMNplZY8z7c82s1czWhY+LMt4/0MxeMbPrImXTzeyp8Jg/NLM+lrkWkYKweTNcFP7vfdllycYi+6zfBGBmlcBC4EzgKOACM4ub9WmJu9eFj5sy3msCVmSU/TdwMTA1fMwYaPAikifPPQdbtsCUKcFdv6+8AhMmJB2V7KNcegAfBTa5+xZ3bwPuBM7J9Q+Y2XSgBngwUnYYcKC7r3Z3BxYD5w4ochHJj64umDED7rkHKiqCu38PPzzpqGQQ5JIAxgFbI9vbwrJM55vZBjO7y8wmAJhZBXA18LWYY27L4ZiY2cVm1mxmza2trTmEKyKDasUKeOEFqK5OOhIZZLkkgLix+czlLZYBk939GOBh4NawfD5wv7tvzdg/l2MGhe43uHu9u9ePGTMmh3BFZFBdfz285z3w+c8nHYkMslzuBN4GRAf7xgPbozu4e3Qu2BuBK8PXfw6cZGbzgWpguJm9BfwgPE7WY4oUil6Lv4dqakp8qmeABx4IFntfsABGjkw6GhlkufQAngSmmtkUMxsOzAKWRncIx/TTZgLPALj7bHef6O6Tga8Ci9290d1TwC4zOz68+ufzwL37/nVEBkdtbTDVjRm0vB1Z+zciLimUnMWLYepU+O53k45EhkC/CcDdO4BLgeUEFfvP3H2jmV1hZjPD3RaY2UYzWw8sAObm8LcvAW4CNgGbgQf2In6RQdOj0o9W7tG1f8vN7bfDI4/A/vsnHYkMAQsuwikO9fX13tzcnHQYUqJsVAo+NQvuWgJ48PqBH8JFx8Owd6B9JPxgC7xV2/2ZIvrfZ2C2bw/m+zniiKQjkUFgZmvcvT6zXHcCS1mLtvp7tPTTr8+bvWft38g6vyXv2mvhQx+C119POhIZQkoAUtKiFXz0UVmZMb5fsx7qFgVr+067Zc/rsRv3rP1b1QbTFvU6F1By3INpnk87DUaPTjoaGUJKAFIy4ir7lhaguvdJ3K6usPzi6TDxsZ4t/cq24BEn0guoqRna75OY1avhpZfggguSjkSGmBKAFJ1srfpsV+tkPYl7aiOMSkGF92zpV3QFD+h9x0pVG3UzV+FewpeA3nEH7LcfnKub80udEoAkIlslXvGeFDavARu1o0d57Z7zrgOr6KtTkaGdyPBNdQqOvj3+lsRMHcOZXz8f/5bj3/LSXhPYHZYuhbPPhlGjko5GhpgSgAy5rEMzMfyk+NZ6v5dlZqvoG5riT+Ke2ggVnXs+31ciqGpj1bZVOX3XomcWLPLyH/+RdCSSB0oAMiiytejHvDdFy5kxrfWYcfmslXhEatcAKvr08TJP4tasj2/9dwyH386Hf3Vq/tu7W/wl3+rPdPDBMHly0lFIHigByD5L7cpSyQOvHZVl/D2uFZ+ttR7RtHIAFf0pl+3ZN806gxO+0dZ/WlUbVVNKfHy/L+++C6efDg8/nHQkkidKALLPmlZmqeT7Gn/PLM9WiUeSSmpXikXrBlDRf+C+PfumVbWx37jNsUM+dbV1tF9XRi39TA88AA89BJ0xyVFKUi6TwYn00GNytOoUfGURDAsr85Xf3HOnbFyL/v6F8eXm8ZV4+jMEiabLc6/o68aPZ+2XNIV4zn76UxgzBk45JelIJE+UAGRAguGecLqEt2qzV/LZWvRrLo4vf/29sZU4E4KTrzU1sHrbato6VdEPibfegmXLYN48qFK1UC70X1oGpMdwz8p/jq/MV36zZ2JIS4+/Z5QP36+Ti85qYOFZv4//oz9Kvyjj4Zmhdu+9sHu3bv4qM0oA0q/uIZ/M4Z5hb2cfthm/Ouv4+zsdPcvbOsvoMstCNXo0nH8+nHBC0pFIHikBSNYFTw6dkuKof5lFy9tLgJjhnmxj7zNXlddlk6XgzDODh5QVJYAy12tMP+K1o5p4/OU+hnv8bbgqVR7TI5ey5cuhvh4OOSTpSCTPdBlomevvEs4u7+r7UsvI50p2crRSlkrBjBnwsY8lHYkkQAmgjHVfVx935+0AhnvcKd+bp4pdY2PwfM01ycYhidAQUBnrcV19f5dwxgz3AKzVkE/xeuqpYM3fs86CU09NOhpJgHoAZSrd+u++rj565222Szgzhok05FPkVqwInq+9Ntk4JDFKAGUq6121fVzCGR3u0ZBPCTjxxOCuX038VrY0BFSmst5Vq0s4y0NHB0ybponfypwSQJlSJV/G/vjHYMrnf/5naCqTRe4lloaARMrNddcFz9OmJRuHJE4JQKScvPxysNrXKafAeeclHY0kTAlApJwsWQJ/+hNceWXSkUgBUAIQKScvvQSTJsH06UlHIgVACUCknFx3HTz9dNJRSIFQAihx2RZrr63t/7NSovbfP+kIpEAoAZSQuMq+5e0UzO29YHvc9M9S4p54Aj75Sdi8OelIpEAoARShbK362Mq+Ictsn1J+Hn8cfvELGDUq6UikQOSUAMxshpk9a2abzKwx5v25ZtZqZuvCx0Vh+SQzWxOWbTSzL0c+82h4zPRnxg7e1yptPRZkj1b4mZV9elK3uNk+pfz87ncwfjyM1f9qEug3AZhZJbAQOBM4CrjAzI6K2XWJu9eFj5vCshRwgrvXAccBjWZ2eOQzsyOfeXXfvkrpydbSB4LK/eLpMPGxoMKPq+zjFmyX8rVmja7+kR5y6QF8FNjk7lvcvQ24Ezgnl4O7e5u7vxtujsjx70koa0sf4NRGGJWCCu+9YIt1wimN8Qu2qxdQnnbtgueegw9/OOlIpIDkUiGPA7ZGtreFZZnON7MNZnaXmU1IF5rZBDPbEB7jSnffHvnMonD455tm3W3bspXZ4u8WN7Rz9O2Q3sc64Njbelb2x9zW55TOmsq5zLz6Khx3HBx/fNKRSAHJJQHEVcyZy4AsAya7+zHAw8Ct3Tu6bw3L3w/MMbN01TPb3Y8GTgofn4v942YXm1mzmTW3trbmEG7xim3xxw3tnNoIFZ17PljVHlTuURWdfU7prKmcy8z73gerV8PppycdiRSQXGYD3QZMiGyPB6KteNx9Z2TzRqDXfebuvt3MNhJU9ne5+yth+S4zu4NgqGlxzOduAG4AqK+vL4/1p6ItfvPeQzvR1n9azHZdbZ1m/ZSAe0a3UiS3BPAkMNXMpgCvALOAz0R3MLPD3D0Vbs4EngnLxwM73X23mR0MnAj8p5lVAQe5+2tmNgw4m6DnINUpmHZL0OL/8E3gFb2Hdio6Yz+qCl9itbXBQQfB/Plw1VVJRyMFpN8E4O4dZnYpsByoBG5x941mdgXQ7O5LgQVmNhPoAF4H5oYfPxK42sycoI16lbs/ZWYHAMvDyr+SoPK/cZC/W9GorQ2Hf9JX9lSG580r24IEEFXRGTsop8pfsrrvPti9W9f/Sy/mXjyjKvX19d7c3Jx0GIOishK6Ms7Rcu4cOHZx/FmXCFX2MiDnnQcPPQQ7d8Lw4UlHIwkwszXuXp9ZrhXB8qi7pR9Vsw7mfQzuvDt+bN+BJ+fj9y3MT5BSWjo7YdkyOOMMVf7Si67Lz6PYq3zO+yyMeBP++tPxY/sGTHkkn2FKKXnllWD9349/POlIpACpB5CE9FU+ZyyAsRuDSn7kG9kvuO3SfybZS4ccAv/3f3DEEUlHIgVINUu+Ra/r/7Of97+/AaM1e6PspQMOgBNOSDoKKVAaAsq3HvPzELmbt/euVa11+Lcc//bufEUnpaSzEz796eAcgEgM9QDyKd36T1/XH3fC99U/w//r9/mOTErRiy/CXXfBX/1V0pFIgVIPYAh0z+kzKoXNa8BG7Qhuwoy2/uMY7DdOwz0ySJ57Lnj+wAeSjUMKlnoAQ6D7ap/olA73L4Txq3vPzwOQqsN/pOv6ZZAtXx486wSwZKEewFCpTkFdOKXDtFuCSz5vvx/a9+u5X/tIuP2BZGKU0vbkkzBpUnAlkEgM9QCGSkMTVLYHryvbek/slmadjJzRRLDmjsgg+uAH4bTTNAmcZKUEMBSirX/Y0wt4/X2xUzQfceqq/Mcope/mm5OOQAqcEsBQiLb+0yrb4KUGXeEjIgVD5wCGwsSVe1r/aRVdMGlFMvFI+fnVr2DcOHjqqaQjkQKmBDAERraeDB0ZE291DGdka0MyAUn52bwZtm+HQw9NOhIpYEoAQ+CIU2Mu99RYv+TTCy/AiBFa/Fn6pHMAQ0Bz9UviXngBJk+GCrXxJDv96xApRS++CFOmJB2FFDj1APZR7CIvBD3vHTvyH48IECwAM2FC0lFIgVMPYC+kdqVo+HEDO97aEVv5Q5AUovuJ5NV3vgOXXJJ0FFLglAAGID3J2+Gfa2Tliys57LONfe7ftLKJx19+nKYVTXmKUARYtw6uvz7pKKQIKAEMQEsLwV2+x4Zr9x57WzDHT5zqFIvWLaLLu1i0bpF6AZI/ixfDNdfAbq0jIX1TAhioUxrBwrV7rTPYjtPQRJcHN4N1eqd6AZI///u/0NoKI0cmHYkUOCWAgYi2/iF7LyBc+KWtM7gXoK2zTb0AyZ+XXoLx45OOQoqAEsBARFv/aXG9gJiFX9QLkLzYuRPeeAPmzEk6EikCSgAD8YH7ei/jaGF5VMzCL22dbazapjuBZYg9/3zwPHVqsnFIUdB9AAPxpwlwwM6Y8vG4Rwt0J7AkZNo0WL8+WAhGpB/qAWTRva5v5MHt90NXZc8duyq1opcUjhEj4Jhj4D3vSToSKQJKAFl03+BVsw4aD4KxG3I/ByAiUgQ0BNSf8z4LI96E8z8Do7bHnwM4YmkSkYmI7BP1APpSsw7Gbgwq+bEbgzV943hlfLmISAFTAujLeZ/tuV2V5c7KtzTnuogUn5wSgJnNMLNnzWyTmfUa8DazuWbWambrwsdFYfkkM1sTlm00sy9HPjPdzJ4Kj/lDM8scXElWtPUPwXPVu9AxrOd+WulLRIpUv+cAzKwSWAicBmwDnjSzpe7+dMauS9z90oyyFHCCu79rZtXA78PPbgf+G7gY+A1wPzADKJzLaTJb/2lVGYu9a6UvESlSuZwE/iiwyd23AJjZncA5QGYC6MXdo3dDjSDscZjZYcCB7r463F4MnEshJYDRm+JP+LaPwL/9ThIRiYgMqlyGgMYBWyPb28KyTOeb2QYzu8vMuleiMLMJZrYhPMaVYet/XHic/o6JmV1sZs1m1tza2ppDuINk3RdiF3Zn7YX5i0FEZAjlkgDixuYzL4dZBkx292OAh4Fbu3d03xqWvx+YY2Y1OR4z/fkb3L3e3evHjBmTQ7iDJGY6B6raYIKGe0SkNOSSALYB0bXlxgPbozu4+053fzfcvBGYnnmQsOW/ETgpPGZ0usJex0zc7fdD+349y9pH6q5fESkZuSSAJ4GpZjbFzIYDs4Aedz6FY/ppM4FnwvLxZjYyfH0wcCLwrLungF1mdnx49c/ngXv3+dsMUPd0D6NS2LwGbNQOzILyuBk9sU44WTN6ikhp6DcBuHsHcCmwnKBi/5m7bzSzK8xsZrjbgvAyz/XAAmBuWH4k8ERYvgK4yt2fCt+7BLgJ2ARsJoETwN3TPTQ0wcTHuyv3lhY0BCQiJc/cs9zdWoDq6+u9ubl50I5nRrB4y1feC8PeCYZ4frAF3qrtWZ4Wvu+7agctBhGRoWZma9y9PrNcdwJHh3oiQzwjZ8QPAY2coSEgESkN5Z0AwqUbu4d6qtpg2iKo3sERp8YPAemmLxEpFeU9G2gfJ3rXfkmLuohIaSvpHkDcoi7dV/mATvSKSFkr6QTQfZVPdQrmNkD1jp7l16+Fq7bvud6/fSRclQrKRURKXEkngG4Zl3n2ei/mJLCISKkr/QSQPtFb0dV9grfXezEngUVESl3pJ4C+Wvi621dEymOEpCIAAARjSURBVFhpJ4B+WvhVk+NPAldN0UlgESl9JX0Z6MgZTezOejPXQtqv08leESlfJd0D0M1cIiLZlXQPQDdziYhkV9I9ABERyU4JQESkTCkBiIiUKSUAEZEypQQgIlKmimpFMDNrBV7K0587FHgtT39rMBVr3FC8sRdr3FC8sSvugZnk7mMyC4sqAeSTmTXHLaFW6Io1bije2Is1bije2BX34NAQkIhImVICEBEpU0oA2d2QdAB7qVjjhuKNvVjjhuKNXXEPAp0DEBEpU+oBiIiUKSUAEZEypQQQMrMmM9tgZuvM7EEzOzzLfp3hPuvMbGm+44yJJ9e455jZ8+FjTr7jjGNm/2Fmfwjjv8fMDsqy34tm9lT4HZvzHWdMPLnGPcPMnjWzTWbWmO8445jZp81so5l1mVnWyxEL8DfPNe6C+s3NbLSZPRT+f/eQmR2cZb9k6hV31yM4D3Jg5PUC4EdZ9nsr6VgHGjcwGtgSPh8cvj64AGI/HagKX18JXJllvxeBQ5OOdyBxA5XAZuC9wHBgPXBUAcR+JHAE8ChQ38d+hfab9xt3If7mwPeAxvB1Yx//xhOpV9QDCLn7nyKbBwBFcXY8x7jPAB5y99fd/Q3gIWBGPuLri7s/6O4d4eZvgPFJxpOrHOP+KLDJ3be4extwJ3BOvmLMxt2fcfdnk45joHKMuxB/83OAW8PXtwLnJhhLL0oAEWb2b2a2FZgN/EuW3fYzs2Yz+42ZFcR/zBziHgdsjWxvC8sKyReAB7K858CDZrbGzC7OY0y5yBZ3MfzmfSnk3zybQvzNa9w9BRA+j82yXyL1SkmvCJbJzB4GamPeutzd73X3y4HLzewy4FLgWzH7TnT37Wb2XuDXZvaUu28ewrAHI26L+Wxeejj9xR7ucznQAdye5TAnhr/5WOAhM/uDu68cmogDgxB3Qf/mOSjI37y/Q8SUDflv3lfcAzhM3usVKLME4O6n5rjrHcB9xCQAd98ePm8xs0eBaQTjjkNmEOLeBnwssj2eYCx1yPUXe3hC+mzgFA8HQ2OOkf7NXzWzewi6+kNaGQ1C3NuACZHt8cD2wYswuwH8e+nrGAX3m+cgkd+8r7jNrMXMDnP3lJkdBrya5Rh5r1dAQ0DdzGxqZHMm8IeYfQ42sxHh60OBE4Gn8xNhvFziBpYDp4fxH0xwEnN5PuLri5nNAL4BzHT3/5dlnwPMbFT6NUHsv89flLEx9Rs38CQw1cymmNlwYBaQ+FVjuSjE3zxHhfibLwXSV93NAXr1ZBKtV5I8Q15ID+Bugn/kG4BlwLiwvB64KXx9AvAUwdUFTwEXFkPc4fYXgE3hY17ScYcxbSIYs10XPn4Ulh8O3B++fm/4e68HNhIMBxR83OH2XwHPEbTkEo87jOmTBC3ld4EWYHmR/Ob9xl2IvzlwCPAI8Hz4PDosL4h6RVNBiIiUKQ0BiYiUKSUAEZEypQQgIlKmlABERMqUEoCISJlSAhARKVNKACIiZer/AzGCAVnfv8FSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Took  2.1038260459899902  seconds\n"
     ]
    }
   ],
   "source": [
    "## Make the small plots\n",
    "def add_log(x, y):\n",
    "    if (x > y):\n",
    "        return x + math.log(1.0 + math.exp(y - x))\n",
    "    else:\n",
    "        return y + math.log(1.0 + math.exp(x - y))\n",
    "\n",
    "def calc_map_q(r):\n",
    "    if len(r.alignment_scores) == 0:\n",
    "        return 0.0\n",
    "   \n",
    "    c = 1.38333\n",
    "    \n",
    "    def alignment_score(a):\n",
    "        return a.score1 * c + a.score2 * c + a.log_likelihood\n",
    "    \n",
    "    def modified_alignment_score(a):\n",
    "        #m = a.equivalent_or_better_clusters # Equivalent clusters\n",
    "        #n = a.equivalent_clusters_kept # Equivalent clusters kept\n",
    "        #assert m >= n\n",
    "        assert a.multiplicity >= 0\n",
    "        k = a.multiplicity #+ (m - n)/n\n",
    "        return alignment_score(a) + (math.log(k) if k > 1 else 0.0)\n",
    "    \n",
    "    log_sum_exp = modified_alignment_score(r.alignment_scores[0]) \n",
    "    for a in r.alignment_scores[1:]: # Add from smallest to largest\n",
    "        log_sum_exp = add_log(log_sum_exp, modified_alignment_score(a))\n",
    "    \n",
    "    max_score = alignment_score(r.alignment_scores[-1])\n",
    "    \n",
    "    if max_score == log_sum_exp:\n",
    "        return capper.p_inf\n",
    "\n",
    "    return -10 * math.log10(1.0 - math.exp(max_score - log_sum_exp))\n",
    "\n",
    "def proposed_cap(r):    \n",
    "    map_q = calc_map_q(r)\n",
    "    escape_bonus = 1.0 if map_q < capper.p_inf else 2.0\n",
    "    return int(min(r.xian_cap, (r.faster_cap_precomputed + \n",
    "                                r.pair.faster_cap_precomputed) * escape_bonus, map_q, 120)/2.0)\n",
    "\n",
    "def make_rocs(reads):\n",
    "    # Make ROC curves\n",
    "    roc_unmodified = reads.get_roc()\n",
    "    #print(\"Roc unmodified\", roc_unmodified)\n",
    "    roc_adam_modified = reads.get_roc(map_q_fn=current_cap)\n",
    "    #print(\"Roc adam modified \", roc_adam_modified)\n",
    "    start_time = time.time()\n",
    "    roc_new_sum_modified = reads.get_roc(map_q_fn=proposed_cap)\n",
    "    #print(\"Roc mode modified (time:{}) \".format(time.time()-start_time), roc_new_sum_modified)\n",
    "    #Reads.plot_rocs([roc_unmodified, roc_adam_modified, roc_new_sum_modified])\n",
    "    \n",
    "    PairedReads.plot_rocs([roc_unmodified, roc_adam_modified, roc_new_sum_modified ])\n",
    "\n",
    "start_time = time.time()\n",
    "make_rocs(reads_1kg_small)\n",
    "make_rocs(reads_hgsvc_small)\n",
    "print(\" Took \", time.time() - start_time, \" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start_time = time.time()\n",
    "# Parse the hgsvc reads\n",
    "reads_hgsvc = PairedReads(\"minimizers_hgsvc_paired/minimizers_correct\",\n",
    "                          \"minimizers_hgsvc_paired/minimizers_incorrect\",\n",
    "                          \"minimizers_hgsvc_paired/minimizers_mixed\")\n",
    "print(\"Got \", len(reads_hgsvc.reads), \" in \", time.time() - start_time, \" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Got  2000000  in  1953.9742600917816  seconds\n"
     ]
    }
   ],
   "source": [
    "start_time = time.time()\n",
    "# Parse the 1kg reads\n",
    "reads_1kg = PairedReads(\"minimizers_1kg_paired/minimizers_correct\",\n",
    "                        \"minimizers_1kg_paired/minimizers_incorrect\",\n",
    "                        \"minimizers_1kg_paired/minimizers_mixed\")\n",
    "print(\"Got \", len(reads_1kg.reads), \" in \", time.time() - start_time, \" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de5jVZbn/8ffNwAgpHodmTE5SupPUC3Oi01bcaEVZHilRa3ew1Iusrq2WUFZuJvRn9ktzSx5SMLa6wZ9SYUIKikxuqByaAQQECUXBgcYDBhoOs+b+/fH9rmHNOsysNbNmfdda83ld17pmzfc091ow6577eZ7v85i7IyIikmhA1AGIiEjxUXIQEZEUSg4iIpJCyUFERFIoOYiISIqBUQeQi6qqKh89enTUYYiIlJRVq1a96u7DcjmnpJLD6NGjaWhoiDoMEZGSYmZbcz1HzUoiIpJCyUFERFIoOYiISAolBxERSaHkICIiKZQcRESKUE0NmIENbca+OgEbugOzYHshKDmIiBShnTvDJxPqYOTTcGpd5+19TMlBRKRYHdQM4+bAgHY4aQ4ctKNgP1rJQUSkmLS1waJFwfMJdWDtwXOLdVQPhaDkICJSLN54Az70ITjzzP1Vw8DWYN/A1oJWD0oOIiJRaw0TwKGHwsknwwMPdK4a4gpYPSg5iIhEpa0NbrsNRo+GF18MhifdfTdceCEDR6/cXzXEDWxl4NErChJaSU28JyJSNlasgG9+E5qa4PTTob1zlbDvtsaIAguochARKSR3+PrX4eMfh5YWePBBWLIExoyJOrJOlBxERArBPfhqBgceCNdcA889B5//fLCtyKhZSUQkz2pqYOdbzTB5Cjw0H/bUAEb14a3seK0SbrmlKBNCIlUOIiJ5tnMnKXc2A+x8vTJ4UuSJAZQcRETyL8I7m/NFyUFEJN8ivLM5X5QcRETyoaUFLrmE5vtuj/TO5nzJKjmY2SQz22hmm81sWpr9o8zsCTNbY2ZPmdnwhH03mtmz4eOChO1Hm9mfw2vON7PK/LwkEZECamuD22+HY4+FuXOpe3FupHc250u3ycHMKoBZwKeBscCFZjY26bCfAXPd/URgBnBDeO6ZwAeBccCHgavN7ODwnBuBm939fcAbwCW9fzkiIgV0221QXQ1Tp8IHPwhr1rBy2N60dzYzYgXV1dGE2RPZDGUdD2x29y0AZjYPOBtYn3DMWODK8Pky4LcJ2+vdvQ1oM7M1wCQz+3/AROCi8LhfA9cBt/f8pYiI9KE9e2DxYnj4YfjpT2HkyCAxnHkmXHABfOYzYEbjcV3c2XxH4cLtrWyalY4CXk74flu4LdFq4Lzw+bnAUDM7Itw+yczeZWZVwL8BI4AjgF1h0sh0TQDM7FIzazCzhpaWlmxek4hIfrz1FsydC2efDVVV8IUvwLJl8Pzzwf7Pfz7Yf+aZJTE8NRf56pC+GphgZo3ABGA7EHP3x4FFwArgf4CVQCyXC7v7Xe5e6+61w4YNy1O4IiIZ7NgB68OGkT174CtfgcZGuPxyWL4cXnklmAupzGXTrLSd4K/9uOHhtg7u/gph5WBmBwHnu/uucN9MYGa47wFgE/AacKiZDQyrh5RriogUzNatsGBB8Pjf/4UzzoDHHw+ajdatg/e/v+wqg+5kUzk8AxwTji6qBKYACxMPMLMqM4tfazowO9xeETYvYWYnAicCj7u7E/RNTA7P+TLwu96+GBGRnH3ta8GU2VdeCbt3w3XXwc9/vn//ccf1u8QAWVQO7t5mZlcAjwEVwGx3X2dmM4AGd18InAbcYGYO1APfDE8fBPzRgjf2H8AXE/oZrgHmmdlPgEbgnvy9LBGRJO7B9NgLFsCjjwZNREOHBlXCccfBuefC+94XdZRFwzw+U2AJqK2t9YaGhqjDEJFS8uKLwZDTBQvghRdgwACYMAHuuqvfJAMzW+XutbmcozukRaS8tLXBE08EVQIEncq33hpUB3ffHXQ4P/lkv0kMPaUpu0Wk9O3dC0uXBtXB734Hr78ejDKaMwc+8AF49VU4+OBuLyP7KTmISGlqa4OB4UfY+PGwdi0ccgh87nNw3nnwqU8F+8yUGHpAyUFESscbb8DChUGF0NQEW7ZARQVce22QACZOhEpN05YPSg4iUvz++EeoqwvuTm5rg+HDg+rgrbeCpPCFL0QdYdlRchCR4hO/Ke2MM+CEE+Cdd4KRRlddBeefD7W1/fLeg0JSchCR4rBxYzCp3YIFsGpVsO2mm4LkcPrpsGmTEkIBKTmISDTcg1FEw4ZBayt86EPBHcof/jDceGPQbBQfbqqkUHBKDiJSOO3t8Kc/7Z/H6F3vgmefDTqRH3wQjj8+6E+QyCk5iEhh3H03/OhH0NwMgwYF/QnnnRckjAEDYNKkqCOUBEoOIpJ/77wDS5YE1cGPfhRMbHfIIfDRjwYJ4bOfDb6XoqXkICL5sXcv/P73Qafyo48G/QcHHwyTJwfJ4fOfDx5SEpQcRKTn3ngDWlrg2GODZHDBBXD44cHX88/XTWklTMlBRHKzY0cwf9GCBcEEdhMmBPMaDRsGDQ3B0NOB+mgpdfoXFJHsXXEF/PKXwTDU9743WCBn8uT9+086KbrYJK+UHET6qebdzUx5eArzJ8+n5qCa1AM2bgyqg0cegcWLgw7kj3wEqqqCJqPjj9f9B2VM6zmIlIHm3c1MuHcCO/bs6Hi+esfqjm3DxjRjX52ADd2BDQ2ev+dL06l/4WmOvKAOM6ipAXbuDEYXHX98sG7y978PsRi88krwg774xWAZzRNOUGIoc6ocREpE4l/67s6Uh6dw66Rb+fYfvs3Rhx7N0y89Td3yOhzn6Zee5uIFF7Ph1Q3ULa/j1bEOI5+GU+vAHEb+EUb+Lwxoh5PmQP0P2bmzBl57DWbOhFNPhV/8Ilg6c8SIqF+6REDLhIoUqeRmn6mPTuXOVXdy+cmX4zh3rrqT46qOY33LegbYAGIeY3DFYAD2xvbuv9C+YBuD9nZ+7oABbZXw16/Doll4WywYgVRVVdDXKn1Ly4SKlLDEpiGAuvq6jmqgeXczc5rm0O7tzG6czZzG4Pm6lnU4TsxjALTGWmltb+184YrW4NHx/J3gebxVaGBrUD0ctCNYG0GJQVByECkamZLBnKY5TH9iOu3eDiQlgKTCv532juM6DGgPHh3P07QWWCxochIJqc9BpAgkJ4O39r3V8SHf1t7GfWvu66gO2mnfnxTy1Sc8sBVGrMjTxaQcqHIQKZDkZqNEdfV1KcmgNRZUB/va93Ukhr4JbBxc51T/trHvfoaUnKwqBzObBPwCqADudvf/k7R/FDAbGAa8DnzR3beF+34KnEmQiJYA33F3N7OngCOBf4aX+aS7/73Xr0ikSCU2G806c1bH9njVkJgMeqx5HNyZ3Yd8p7Eod/T8R0p56jY5mFkFMAv4BLANeMbMFrr7+oTDfgbMdfdfm9lE4AbgS2b2MeDjwInhcU8DE4Cnwu8vdncNP5Kyl9xs9MMJP+y48Syxauj6Itl/8Henujovl5Eylk2z0nhgs7tvcfdWYB5wdtIxY4Enw+fLEvY7MBioBA4ABgE7exu0SKlJTAAxj1G3fH/n78ptKzuqhkTjaoLmno5HDxODH3Y4vuE53Ol47Eht2RLpJJtmpaOAlxO+3wZ8OOmY1cB5BE1P5wJDzewId19pZsuAZoKus9vcfUPCeXPMLAY8DPzE09x0YWaXApcCjBw5MrtXJVJEkpuNWmOtnaqHxssyf+jb5b372dVD3oTf/ja421kkB/nqkL4amGBmjQTNRtuBmJm9DzgOGE6QZCaa2SnhORe7+wnAKeHjS+ku7O53uXutu9cOGzYsT+GKFE66ZqPk6qG3EqsCb3d815tBhfD2IcHdziI5yiY5bAcS758fHm7r4O6vuPt57n4S8INw2y6CKuJP7r7H3fcAi4GPhvu3h193Aw8QNF+JlJSuRiDFpWs2ao21smJbfoaOpvQfXHkl1NYGU2GI9FA2yeEZ4BgzO9rMKoEpwMLEA8ysyszi15pOMHIJ4CWCimKgmQ0iqCo2hN9XhecOAj4LPNv7lyOSf90NQY2PQMqk8bJG/Mee8mi8rJGammD+uuRHTZpJUpOl7T/4y1/gllvgE58IFt0R6aFuk4O7twFXAI8BG4AH3X2dmc0ws7PCw04DNprZJqAamBlufwj4G7CWoF9itbs/QtA5/ZiZrQGaCCqRX+XtVYnkUaYEkDwCqavqIZOdGYZnxLdnGlWUcbRRXV2QFG68UbOmSq9o4j2RLjTvbmbMrWPY27aXIQOHsOU7WzqGoE59dCr3NN5Da6yVyopKvn7S1zvdv5CNrj6/c/7VbGyED34wSBDXXpvjyVLONPGeSJ5lGoKaaQRST6qHvHnwwWBBnm99K7oYpGwoOYhk0FUCKMQIpJxdfz00NQUJQqSXlBxEMugqAfT1CKSc7d4dtFGNHh3Nz5eyo1lZRTLoKgF0deNaLqqr03dK5zS9xfr18OEPw/z58JnP5CUuESUHkQzylQC6kpdpLGbODHqvx+tWIckfNSuJlLJNm2DePPjmN7WCm+SVkoNIKbv+ejjgALjqqqgjkTKj5CBSqlpa4IEH4PLL4d3vjjoaKTPqcxApVcOGBTe+qTlJ+oCSg0gpam+HAQPgAx+IOhIpU2pWEilF3/oWfOlLPZhjQyQ7Sg4ipWbbNrj7bjjwQE2uJ31GyUGk1Pz0p0Gz0rRpUUciZUzJQaSUNDfDXXfBl7+sqTKkTyk5iJSSn/8c2tpg+vSoI5Eyp9FKIqXkqqvg5JPhve+NOhIpc6ocREpJTQ1MmRJ1FNIPKDmIlIJXX4VPfhL++teoI5F+QslBpBTccgssXRrMoyRSAEoOIsXujTfg1lth8mTdES0Fo+QgUuxuvTVY6e3aa6OORPoRJQeRYvbmm0GT0rnnwoknRh2N9CMayipSzCor4brrYMKEqCORfiarysHMJpnZRjPbbGYp9+yb2Sgze8LM1pjZU2Y2PGHfT81snZltMLNbzYLJYMzsZDNbG16zY7tIPjXvbmbCvRPYsScf63HuV1MTTGuU6VFTk6cfNGQIfOc7MG5cni4okp1uk4OZVQCzgE8DY4ELzWxs0mE/A+a6+4nADOCG8NyPAR8HTgSOBz4ExP8Euh34BnBM+JjU2xcj/UvyB3+6RFBXX8fTLz1N3fK6vP3cmhrYubPrY7rbn5Xbb4c77sjDhURyl03lMB7Y7O5b3L0VmAecnXTMWODJ8PmyhP0ODAYqgQOAQcBOMzsSONjd/+TuDswFzunVK5F+J/mDP/n75t3NzGmaQ7u3M6dpTq+qh8RKIS8f/N1pa4Pvfx/+8IcC/DCRVNkkh6OAlxO+3xZuS7QaOC98fi4w1MyOcPeVBMmiOXw85u4bwvO3dXNNAMzsUjNrMLOGlpaWLMKVcpOuIkj+4F+9Y3VKIqirr6Pd2wGIeaxX1UNBEkKce9AJvWuX7oaWyORrtNLVwAQzayRoNtoOxMzsfcBxwHCCD/+JZnZKLhd297vcvdbda4cNG5ancKXYdNU3kK5pKPmD/+IFF3f6ftqSacxpmkNrrBWA1lhrj6qHeMVQUD/8IXz3u3D44eqIlshkkxy2AyMSvh8ebuvg7q+4+3nufhLwg3DbLoIq4k/uvsfd9wCLgY+G5w/v6ppSfnJNAPFzkiuC+LbED/51Les6fX/f2vs6kkVcLtVDPCkUtGIAaG2FRYvgggtg+3Y48sgCByASyCY5PAMcY2ZHm1klMAVYmHiAmVWZWfxa04HZ4fOXCCqKgWY2iKCq2ODuzcA/zOwj4Silfwd+l4fXI0UslwSQeE5y01DitkxiHutIFnGtsVZWbFuR9vjk0Ue9TQrV1T08sbISVqyAe+6BwYN7F4RIL3SbHNy9DbgCeAzYADzo7uvMbIaZnRUedhqw0cw2AdXAzHD7Q8DfgLUE/RKr3f2RcN9U4G5gc3jM4ry8IikK6UYS5ZIAEs9Jbhqq31qf8sGfzriacfiPvdOj8bLGtMfmo0Korg66C9xhR0/6vv/4x2CqjMGDgyVARSJkXkILlNfW1npDQ0PUYUgGzbubmfLwFOZPns+M5TO4c9WdXH7y5cw6cxZTH53KPY330BprpbKikq+f9HVmnTmL5t3NjLl1DHvb9nZcZ8jAIWz5zhZmLJ/RcU5c4rk9lc1Q1FxUV/cwGSTavTtY2e300+HBB/MRlkgHM1vl7rW5nKPpM6TX4lXC9Cem8/RLT3d0BiePJErXOZyuiShePazctjKnpqFs5TMx9LhKSPbLX8Lrr8PVV+fhYiK9p8pBUiRWADUHdX+r79RHp3JHwx0MsAHEPEaFVVAxoKKjSjjm8GN4/vXn01YAK7atoGlHU8o1x9WMy9gE1Fv5Gn2Ul4oBYM8eOPpoqK2FxWpdlfzrSeWguZUkxbSl06jfWs+0pdO495x7uzw23i/gODGPAcFf/rFY8Lw11sr6lvU4nf8IiVcAfZUAEuWzGSlvCSHRTTcFi/n8+Md5vrBIz6lZqZ/KNKy0eXcz96+9H4D71tzX7X0BdfV1xNpjXR4zqGIQU2unZt05nG+9SQyJncx5a0JK9re/BbOufuQjfXBxkZ5RcuinMg0rnbZ0WqcKYNrSlHkWO8Srhn3t+7r8WfnoJyiUxETQZ8kg2X33wfz5BfhBItlTcuhH4tVCuqkm4vvjVUNcV9VDd/cbJA4lLVSV0Bs9vjehp9zh5XBmmkGDCvzDRbqm5NCPxKuF5Kkm4tVDYtUQ11X1kG40EexPCqWQEApeJSRasiQYvrpsWYF/sEj3NFqpn0h3P0Fc/L6C4395PK/987WU/VVDqmj5XulOetjV6KRI//ufeSasWgVbt8IBB0QYiJQ73ecgGXXVBBSvHkYcMiLt/uGHDE+7vVRkai4qeDNSoo0bgzmUpk5VYpCipKGs/UDyNBTJCjmsNAoFby7Kxq23BvMoXXZZ1JGIpKXKoUx0N+NpctVQWVHZaXhpuSaGorR3LzzwAFx0UcTli0hmSg4lKNflMPtqGgrpocGD4dln4T//M+pIRDJSh3QJmvro1E6T2iV2Nsc7l7OZ9kJE+gd1SPcD6aa+zudymNLHFiyAiROLtCNEZD8lhxKTnAjytRymFMgtt8ALL4CWvJUip+RQQtItftPb5TClgFatChb0+fa3oaIi6mhEuqShrCUk09oH8RlQ49TZXKR+8Qs46CD42teijkSkW0oOJaSr6So0FLXI7dgB8+bB5ZfDIYdEHY1It5QcSogSQAk78EC48Ub43OeijkQkK0oOIoUwdCj8x39EHYVI1tQhLdLXfv97uPdeiHW9KJJIMVFykKJXUxPMrJrpUVPM9/u5w/XXQ10dDNCvm5QO/W+VotfdMp/5Wh+6TyxeDCtXwre+1fXc4SJFJqvkYGaTzGyjmW02s5SVX8xslJk9YWZrzOwpMxsebv83M2tKeOw1s3PCffea2QsJ+8bl96VJKUpXJZS0+fODuZS++c2oIxHJSbcd0mZWAcwCPgFsA54xs4Xuvj7hsJ8Bc93912Y2EbgB+JK7LwPGhdc5HNgMPJ5w3nfd/aH8vBQpRTU1Rf6Xf280NcHcufCNb2gZUCk52VQO44HN7r7F3VuBecDZSceMBZ4Mny9Lsx9gMrDY3d/uabBS2tJVBWWbGADGjIFPfhKuvTbqSERylk1yOAp4OeH7beG2RKuB88Ln5wJDzeyIpGOmAP+TtG1m2BR1s5mlXQ7LzC41swYza2hpKd2lKvuLTJ3HZV0hZHLwwfDYYzByZNSRiOQsXx3SVwMTzKwRmABsBzrG7ZnZkcAJwGMJ50wH3g98CDgcuCbdhd39LnevdffaYZqsrCj0JAH0ZWIouvVympvhU58K1mwQKVHZ3AS3HUhcXHh4uK2Du79CWDmY2UHA+e6+K+GQLwC/cfd9Cec0h0/fMbM5BAlGSkAUCSBR0S9B8qtfweOPa21oKWnZVA7PAMeY2dFmVknQPLQw8QAzqzKz+LWmA7OTrnEhSU1KYTWBmRlwDqA/s7LU1ZKgfSWxWohS0VUJyfbtgzvuCCqHY46JOhqRHus2Obh7G3AFQZPQBuBBd19nZjPM7KzwsNOAjWa2CagGZsbPN7PRBJXH8qRL329ma4G1QBXwk169kn5k2tJp1G+tZ9rSlFHFfabQ/QXV1UGFkPwo+jVyVq0KmpW++tWoIxHpFS0TGrHm3c1MeXgK8yfPz2ppz+bdzYy4eQQxj1FhFWy7cltBlgTNR8VQXZ05yVRXl8AHfzYeeQTOOitYt+Ff/zXqaEQALRNakurq63j6paezXpxn2tJpxDzo6495rM+rh3hzUrYyNfvEP/zTVQMlURFk68ADYdIkqM3p91Ck6KhyiFDz7mbG3DqGvW17GTJwCFu+s6XLKiCxaojr6+oh18RQNh/yImVElUOJSV4PurvqIbFqiCtE9dCdsvvrv6feegtefz3qKETyQskhIunWg57TNKfLEUiPPv9o+u2b0m8vhKIfPVRICxZAVRU891zUkYj0mpJDRDKtB91V9TDikBFptw8/ZHheY8uWqoUkS5bAEUfAscdGHYlIr2kluIikWw+6NdbKim0rMp6jZUKLmDssXQqnn651G6Qs6H9xH+rqZrVFFy2iwio6bauwChZfvLhQ4WWlq9FHkmD9+uD+hk98IupIRPJCyaEPdXWzWrF2LifLNPxUzUlJliwJvp5xRrRxiOSJkkMfad7dzP1r7wfgvjX3pVQPj2x6JO15CzcuTLtditw558Ds2TBqVNSRiOSFkkMf6e5mtQGW/q1PbmqSEjF6tKbMkLKi5JCDbCe8S6wa4pKrh2IbeSS9sGED/Pd/B/c5iJQJjVbKQWIfwr3n3NvlcZn6E+LnaeRRGXngAbj+evjc56KORCRvVDlkqbs+hETqT+hnli6F8ePh0EOjjkQkb/pFcsjH+geFnvBOSsSbb8Jf/gITJ0YdiUhe9YvkkOvMp8my6UNIVEqdzV0t+SlZmDMH2tvh7LOjjkQkr8o+OTTvbmZ242zavZ3ZTbN7VD3kek9CKXU2R73kZ8lbuzZoUho/PupIRPKq7Duk6+rr2NceLF3dGmulbnkds86cldM1cp3wTp3N/cg998A//hF1FCJ5V9aVQ2LVAPS4eiilSkAK6J//DL4efHC0cYj0gbKuHBKrhrieVA/lVgnU1KjZqNeWL4dLLw3ezCVLoLIy6ohE8qqsK4f6rfUp02K3ezvLty6PKKK+kWunshJDL730UjD7KsBttykxSFkq6+Rw6qhTqazo/ItbWVHJhFETIoqodzIlAXUqF1hzM8RicPPNcMIJUUcj0ifKulmpJ2smFAs1/RSxt9+OOgKRPlfWlcOiixYxeODgTtuGDBxSdGsmpKPEUMRaW4PV3jR8VcpYVsnBzCaZ2UYz22xmKYP7zWyUmT1hZmvM7CkzGx5u/zcza0p47DWzc8J9R5vZn8NrzjezvDfc9mQpTgloMZ8ufOpTwWR7VVVRRyLSZ7pNDmZWAcwCPg2MBS40s7FJh/0MmOvuJwIzgBsA3H2Zu49z93HAROBt4PHwnBuBm939fcAbwCV5eD2dlHKzUl/qanU3LebThd274c47oa1NS4FK2cumz2E8sNndtwCY2TzgbGB9wjFjgSvD58uA36a5zmRgsbu/bWZGkCwuCvf9GrgOuD3XF9CVchuCmqtMSUAf/D1UVwc33QS1tXDyyVFHI9Knsvnz5yjg5YTvt4XbEq0GzgufnwsMNbMjko6ZAvxP+PwIYJe7t3VxTQDM7FIzazCzhpaWlizC7X8S/+LXUp59ZP36YHTSJZcoMUi/kK/a+Gpggpk1AhOA7UDHZERmdiRwAvBYrhd297vcvdbda4cNG5ancKPXmwnvlAQKzB2uuAKGDoUbbog6GpGCyKZZaTuQOH/E8HBbB3d/hbByMLODgPPdfVfCIV8AfuPu8duVXwMONbOBYfWQcs1SlmkYanX1/g9y3ZtQQubPh2XL4PbboYz+QBHpSjaVwzPAMeHookqC5qFOq9aYWZVZxzzV04HZSde4kP1NSri7E/RNTA43fRn4Xe7hF6d8fPB31WksBXbUUXDRRfCNb0QdiUjBdJscwr/sryBoEtoAPOju68xshpmdFR52GrDRzDYB1cDM+PlmNpqg8kies+Ia4Eoz20zQB3FPr15JmdmxQ/0IReOUU+D++6Gi+NbjEOkrFvwRXxpqa2u9oaEh5/Oadzcz5eEpzJ88n5qD+n4VG7PM++JvdzbHSMTWrw+m5L7uuqC/QaREmdkqd6/N5Zx+MVh72tJp1G+t19Kekr14J/ScOfDOO1FHI1JwZZ8cEpf47Gppz0JTn0KRe+ihoBP6+ut1J7T0S2WfHBKX+Oxqac98yuaDX30KRe4PfwgW8VEntPRTZZ0cEquGuEJUD/rgL3EvvQQLFsC//Is6oaXfKuvkkFg1xBWqepASdtRRwSpvDzwQdSQikSnr9Rwe2fRI2u0LNy5Mu10ECKqFG2+MOgqRSJV15TDA0r+8ClNTgaTxzjvwsY/BokVRRyISubJODpnuaag+SEOCJI1582DlSk3HLUKZJ4dTR53KgKSXOIABOa8h3ZtJ8qSELF4Mhx4KZ5wRdSQikSvr5FC/tZ52Oq8E1047y7cmz+TRNU2S1w9s3AgPPwxf+QoMLOuuOJGslHVyOPnI9PPu174np7vIpT+46ip417tg+vSoIxEpCmWdHB59/tH02zel3y792KWXBlNzv/vdUUciUhTKun4eccgIXvvnaynbhx8yPIJopKiddVb3x4j0I2VdOSy6aFHKsNUKq2DxxYsjikiK0urV8NRT0N7e7aEi/UVZJwfdIS1Z+a//gsmTu55HXaSfKevkoD4HyUpTE5x0kpKDSIKyTg66CU66tW8frF0L48ZFHYlIUSnr5HDqqFNTptAYYLnfBNcV3QhX4p57Dlpbg8pBRDqUdXKo31pPuyfdBOe53wTXFd0IV+KamoKvqhxEOinr5JCv6TOkjE2eDH/+Mxx7bNSRiBSVsr7PIV/TZ0gZGzIExo+POgqRolPWlYOmz5Au7d0L3/se/OUvUUciUnTKOjnkayhrpjWhpcRdfz3cdBM0N0cdiSx6DfgAAAmBSURBVEjRySo5mNkkM9toZpvNLOUOMjMbZWZPmNkaM3vKzIYn7BtpZo+b2QYzW29mo8Pt95rZC2bWFD7y3iM44pARabfnOn2G1n4uUzt2wBFHwNlnRx2JSNHpNjmYWQUwC/g0MBa40MzGJh32M2Cuu58IzABuSNg3F7jJ3Y8DxgN/T9j3XXcfFz6aevE60lp00SIOqDig07bBAwdr+gwJmGl6bpEMsqkcxgOb3X2Lu7cC84DkP7XGAk+Gz5fF94dJZKC7LwFw9z3u/nZeIs9CXX0d+9r3ddrWGmulbnldoUKQYjZggOZTEskgm+RwFPBywvfbwm2JVgPnhc/PBYaa2RHAscAuM1tgZo1mdlNYicTNDJuibjazA0jDzC41swYza2hpacnqRcUV4j4HKWGjR+v+BpEM8tUhfTUwwcwagQnAdiBGMFT2lHD/h4AxwFfCc6YD7w+3Hw5ck+7C7n6Xu9e6e+2wYcNyCurUUadSWVHZaVtlRaXuc+jvYjF480245hp4/PGooxEpStkkh+1AYs/u8HBbB3d/xd3Pc/eTgB+E23YRVBlNYZNUG/Bb4IPh/mYPvAPMIWi+yquV21bSGmvttK011sqKbSvy/aOklCxbBpddppEGIl3IJjk8AxxjZkebWSUwBViYeICZVZl1TGI0HZidcO6hZhb/k38isD4858jwqwHnAM/25oWks+iiRQweOLjTtiEDh6hDur/bsydY9W379u6PFemnuk0O4V/8VwCPARuAB919nZnNMLP48lmnARvNbBNQDcwMz40RNCk9YWZrAQN+FZ5zf7htLVAF/CRvrypUV1+X0ucQ85g6pPu7TZuCr2PGRBuHSBHLahyfuy8CFiVt+1HC84eAhzKcuwQ4Mc32iTlF2gP5bFaqrk4/yZ5ukCtBa9bAiBFw2GFRRyJStMp6kHfjZY15u5aap8vImjVwwglRRyFS1Mo6OYikcIeRI+GUU6KORKSoKTlI/2IGv/991FGIFL2ynnhPRER6pqyTQ01N8Idi8kNLe/Zj114bLAnqHnUkIkWtrJNDpiU8tbRnP/bXvwaJwSzqSESKWlknB5EUGqkkkhUlB+k/3n47uCv6uOOijkSk6Ck5SP8RiwVfD0g7AbCIJFBykP5j0CD493+HsclrVYlIMt3nIP3H4MHw619HHYVISVDlICIiKZQcREQkhZKDiIikUHIQEZEUSg4iIpJCyUFERFIoOYiISAolBxERSaHkICIiKco6OVRX57ZdREQCZT19xo4dUUcgIlKayrpyEBGRnskqOZjZJDPbaGabzWxamv2jzOwJM1tjZk+Z2fCEfSPN7HEz22Bm681sdLj9aDP7c3jN+WZWma8XJSIivdNtcjCzCmAW8GlgLHChmSXPefwzYK67nwjMAG5I2DcXuMndjwPGA38Pt98I3Ozu7wPeAC7pzQsREZH8yaZyGA9sdvct7t4KzAPOTjpmLPBk+HxZfH+YRAa6+xIAd9/j7m+bmQETgYfCc34NnNOrVyIiInmTTXI4Cng54ftt4bZEq4HzwufnAkPN7AjgWGCXmS0ws0YzuymsRI4Adrl7WxfXBMDMLjWzBjNraGlpye5ViYhIr+RrtNLVwG1m9hWgHtgOxMLrnwKcBLwEzAe+Avwu2wu7+13AXQBm1mJmW/MUc1+pAl6NOohuKMb8UIz5UQoxQmnEmSnGUbleKJvksB0YkfD98HBbB3d/hbByMLODgPPdfZeZbQOa3H1LuO+3wEeA2cChZjYwrB5SrpmOuw/LIt5ImVmDu9dGHUdXFGN+KMb8KIUYoTTizGeM2TQrPQMcE44uqgSmAAuTAqoys/i1phN8+MfPPdTM4h/qE4H17u4EfROTw+1fJodqQkRE+la3ySH8y/4K4DFgA/Cgu68zsxlmdlZ42GnARjPbBFQDM8NzYwRNTk+Y2VrAgF+F51wDXGlmmwn6IO7J26sSEZFeyarPwd0XAYuStv0o4flD7B95lHzuEuDENNu3EIyEKjd3RR1AFhRjfijG/CiFGKE04sxbjBa08IiIiOyn6TNERCSFkoOIiKRQcuglM7vOzLabWVP4+EyG4140s7XhMQ1FGmOXc2gVgpldZWZuZlUZ9scSXsfCdMf0tSxi/LKZPR8+vlzg2OrCOc6awjnN3pPhuMjexxxijOx9DH/+TWb2XBjrb8zs0AzHRfm7nW2Muf9uu7sevXgA1wFXZ3Hci0BVscYIVAB/A8YAlQR3vY8tcJwjCEbFbc30XgF7Iv737jJG4HBgS/j1sPD5YQWM7+CE598G7ii29zGbGKN+H8MYPkkw/Q8Ec8HdmOG4KH+3u42xp7/bqhwkLps5tPrazcD3gGIeJdFdjJ8Clrj76+7+BrAEmFSo4Nz9HwnfHkgRvpdZxhjp+wjg7o/7/il+/kRws25RyTLGHv1uKznkxxVhWTfbzA7LcIwDj5vZKjO7tJDBhbqLMZs5tPqMmZ0NbHf31d0cOjica+tPZlbQyRqzjDHS9xHAzGaa2cvAxcCPMhwW2fsIWcUY+fuY5GvA4gz7ov7djssUY4/ey7JeCS5fzGwpUJNm1w+A24E6gv8gdcD/JfhHSvav7r7dzN4NLDGz59y9vshi7FPdxPh9ghK5O6PC93EM8KSZrXX3vxVZjH2qqxjd/Xfu/gPgB2Y2neAG1h+nOTay9zGHGPtcd3GGx/wAaAPuz3CZyH63c4gxZ0oOWXD3M7I5zsx+Bfw+wzW2h1//bma/ISj18vYfKA8xdjuHVm9litHMTgCOBlabWfxn/9XMxrt7p8VeE97HLWb2FMGkjnn7UMtDjNsJZgyIGw48la/4uooxjfsJbl5N+eCN6n3MIcY+fx+h+zgtmEz0s8DpHjbgp7lGpL/bWcTYs9/tKDpRyukBHJnw/D+AeWmOORAYmvB8BTCpyGIcSNDpdzT7O60+ENF7+iLpO3sPAw4In1cBz1PgTvMsYjwceCGM9bDw+eEFjOuYhOffAh4qtvcxyxgjfR/DGCYB64FhXRwT9e92NjH26He7YG90uT6A/wbWAmsIJiQ8Mtz+HmBR+HxM+A+yGlhHUA4WVYzh958BNhH8BVnQGJPi7fjgBWqBu8PnHwtfx+rw6yXFFmP4/deAzeHjqwWO62Hg2fDf+hHgqGJ7H7OJMer3Mfz5mwna6pvCxx3h9mL63e42xvD7nH+3NX2GiIik0GglERFJoeQgIiIplBxERCSFkoOIiKRQchARkRRKDiIikkLJQUREUvx/MVo5NiEkX5wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dfZhT9Zn/8ffNwCgVVBCcUUEc6lPpamGd+tP6q2O3u1ttrY/tSrUPPq221K3d6uVC7a69nLW21ba/PmDVtszKjlVa1K2iroqKoGAvQQYQREVaZCBD0UqBKoZJ7t8fJxkymSSTGZKcPHxe15VrknNOkjsBbu65z/d8v+buiIhI9RoSdgAiIlJcSvQiIlVOiV5EpMop0YuIVDklehGRKjc07ADSjRkzxo844oiwwxARqSjLli17093HZtpXdon+iCOOYOnSpWGHISJSUcxsQ7Z9at2IiFQ5JXoRkSqnRC8iUuWU6EVEqpwSvYhIlVOiFxEJUWMjmIGNjGCXtGAjuzALtheKEr2ISIi2bEncaWmFw5+FU1t7by8AJXoRkbBEo8HPERGY3AZD4jClDUZ0FfRtlOhFRIopFoP162HnzuDxokVwxhlw5JHwvvcF21paweLBfYv1VPWFokQvIlJIGzbANdfApz8Nxx4Lw4fD+98PTz0V7I/Fgr7MCSfAjBl7qvmhiep+aLTgVX3ZTYEgIlJxNm+Gv/wFPvABePdduO02OOoo+OAH4ZxzgvuTJwfHnnYavPjinucun7anmk/qqepnFiQ8JXoRkcF67TX4/vdh9mw49VR44gk45hj4619hSH4Nk6FHLKE7Wc33bIwytGlxwcJUohcRGaiODvjOd2DuXKivh8sug2uvDfaZBbc87f7Z8iIFuYcSvYhIPtyD25AhMH8+PPYYTJ8OV18NDQ1hR5eTTsaKiOQSj8MDD8BJJ8HddwfbvvIVeOONoKov8yQPSvQiIplFo9DWBpMmwXnnwZtv7hkOud9+cMAB4cY3AGrdiEhNa2zMfBVqQ/1OuqKXwoc+BPfeC+efD0MrM2WqoheRmtaT5EdE4OKWnvHrW6Kj4dFHYflyuOCCik3yoEQvIhJIm2sGgNNPH9AImnKlRC8itWn3brj00uB+keeaCZsSvYjUlg2JNbSHDQtOsELR55oJmxK9iFQ/9+Cq1Y9/PJiOYOPGYPvvfleSuWbCpkQvItUrFguuXv3wh+Ef/xHWroWbb4ZRo4L9Zgw/vTXjXDPDT6+eqj6vRG9mp5vZK2a2zsymZ9g/wcyeNLOVZrbAzMal7Pu+ma02s5fN7CdmVXBmQ0TKm3vwc/NmmDoVtm+HX/wimC74mmtgxIieQ4/5+yV7qvmkoVGO+fvCzTUTtn7HC5lZHcEUav8AdAIvmNmD7r4m5bBbgdnufpeZ/R1wM/AFM/sIcApwfOK4Z4EWYEHhPoKISMKrr8JPfwqRSFDJjx8Pzz4bVPR1dRmfsvzK4s81E7Z8KvoTgXXuvt7do8C9wNlpx0wCEpMt83TKfgf2BeqBfYBhQAEXyBKRmucOjz8On/pUMHPknXfCyJHQ3R3sP+mkrEm+VuST6A8DNqY87kxsS7UCOC9x/1xgpJkd5O5LCBJ/JHF7zN1fTn8DM7vCzJaa2dKtW7cO9DOISC274w74xCdg2TL49reDOWja2ir6AqdCK9Q3cS3wMzO7GFgIbAJiZnYk8AEg2bN/wsw+6u6LUp/s7ncCdwI0Nzd7gWISkWq0cSPMnAknnhjMQfNP/xTMQXPBBbDPPmFHV5byqeg3AeNTHo9LbOvh7pvd/Tx3nwJcn9i2jaC6f97dd7r7TuBR4OSCRC4itcMdFi8OknlTE9xyy55VmkaPhi9+UUk+h3wS/QvAUWbWZGb1wFTgwdQDzGyMmSVfawYwK3H/DaDFzIaa2TCCE7F9WjciIjl96UtwyilBL/5f/zUYPfOf/xl2VBWj30Tv7t3AVcBjBEn6N+6+2sxuNLOzEoedBrxiZq8CDcBNie1zgdeBVQR9/BXu/lBhP4KIVJ2tW4Px7jt2BI/POSdo12zcGFTzEyaEG1+FMffyaok3Nzf70qVLww5DRMKwahX8+MfQ3g7vvQf33Rf04aVfZrbM3Zsz7dNpaREJ3zvvwKc/DU89BcOHw8UXw9e+Fiz6IXtNUyCISDi2bw/WXYVg1ExDA3zve9DZCbffriRfQKroRaS0Xn89uHp11qygPbN5Mxx0EPz612FHVrVU0YtIaaxdC2efHcweedttwf3nnguSvBSVKnoRKZ5334U//xkOOwzq6+H3v4dvfQu+8hU45JCwo6sZSvQiUnibNwdV++23w8knw0MPwcSJQf9dUxOUnL5xESmcF1+EH/wAfvObYC74s86Cr399z34l+VDoWxeRvbN7NwwZEswQ+fDDQfV+1VXwL/8SVPESOp2MFZHBeest+O53g2R+//3BtquvDtozP/qRknwZUUUvIgOzejX85Cfw3/8dnGz9+Mf3nFjdf/9wY5OMlOhFJH/ucP75sGEDfP7zwdWrxx0XdlTSDyV6Eclu5064666gen/ySdhvP7j77mBSsTFjwo5O8qQevYgA0NgIZonbyAh2SQt2yE4arzo/OCASCX6ecIKSfIVRohepEZEdEVr+q4UVXSs46ZcncfKvTqZrZxdjJwZJfctfu2BEBC5ugY/PgMOfhVNb2UIjPP88HHlk2B9BBkmtG5EqFNkRYep9U5nzmTm4O1Pvm0rTgU08+8azXHT/RazeuhqA1mdaeXOS9yR1zOHwRXD4czAkDlPaYOG/A43hfiDZK5qPXqQKpCb2xhGNTHt4Gncsu4Mvn/BlHOf2pbczxIYQ81jvJ+7eBzAYtgt27xtsG7YLPNhMdz28eDn+8MwSfyIZqFzz0at1I1IFWhe28uwbz9L6TCuRHRHaOtqIe5xZy2fRtrwNx/smeYCh70FdNLhfF4W694L7ltwfhSltdO3sKsnnkOJQohepAMn+eqaEm5rY2zramPHkDOIeByAai/Je7L3sL2wELRoSP4dk+A3fYrQ+01qATyFhUaIXqQCpFXumfcnE3h3vpn1lO9FYUKXHiePsZXt2aJTFnYv37jUkVEr0ImVgIBV76jHJfcnEvju+m1g8Q4tmkCY3TsZvcJZfubxgrymlp1E3ImUgtWKf+amZffYlK/aYx3odk7qvh5FbZDLckT1xl9n4DCkAjboRCVlkR4SJP5nIru5dDB86nPVXr6dxRGOffUmpx0y5YwodXR0ZXjR3Ms+moQG6dN61IuUadaOKXiRkA63Y390V45ALWuGRmTQ0LMe7gHvuwS783KBjKLN6TwpMPXqRAsrVa892fGqPPRqL9urDL+lc0rOvx9AojA9Ojm7ZArS3BxOMiWShil6kgHL12rMdn16xp1b1yZOglqvv/s1vQksLPD24mBsaBvc8qRyq6EUKJNfomGwyVezRWJTbHlqMWTDRWL+eeQbmzcuasBsagtZMtpt68tVPFb1IgeTqtWfTX8W+ZUseb9zUBChhS3aq6EUKoL9eu0iYlOhFCiBXr70QcrVlRPqjRC9SANl67YWaOqCrC/ynP8Mx/Oxz8F3vqb8ueVOPXqQf6VMAZ1KSKQIOPxwuuABmz4b6+uK/n1QNVfQi/cg1oVh/ei3Pl3JLH02TvTXj8OKLwYOzzoJ771WSlwFTohdJk3rR02CGTKbKNmomfXtXV5ahj7OfCNZo/d3vBvlpRNS6EekjtYJ3fMBDJgvGHW68EcaNg9NPL817SlVSRS+SItPqTKENmVywAJ57DqZPh332Kc17SlVSohdJkTpMMhqLEo33HklTyCGT/QfTCoccApddVpr3k6qlRC+SkH7RU5x4n7HxhRwymTuYCKxYAdddB/vuW/z3k6qmHr1IQqaLnurr6rl8yuWD7ss3NGQ+IdvvhU6HHAJ//CMM1T9R2Xt5VfRmdrqZvWJm68xseob9E8zsSTNbaWYLzGxcYvvHzKwj5bbLzM4p9IcQKYRiXPSUdTRNrjb/W29BLAYjR8Lw4YN+b5GkfleYMrM64FXgH4BO4AXgc+6+JuWY3wLz3P0uM/s74BJ3/0La64wG1gHj3P2dbO+nFaak5p15JmzfHsxKmXN+YpE9cq0wlU9FfyKwzt3Xu3sUuBc4O+2YScBTiftPZ9gP8Bng0VxJXqTmLVsGDz8cDKdUkpcCySfRHwZsTHncmdiWagVwXuL+ucBIMzso7ZipwD2Z3sDMrjCzpWa2dOvWrXmEJFKlWlvhwAPhqqvCjkSqSKFG3VwLtJjZcqAF2ATEkjvN7BDgOOCxTE929zvdvdndm8eOHVugkEQqzIoVwRWwX/867L9/2NFIFcnnlP4mYHzK43GJbT3cfTOJit7MRgDnu/u2lEP+CXjA3XfvXbgiVeyXvwwS/Ne+FnYkUmXyqehfAI4ysyYzqydowTyYeoCZjTGz5GvNAGalvcbnyNK2EZGEH/0IFi2CUaPCjkSqTL+J3t27gasI2i4vA79x99VmdqOZnZU47DTgFTN7FWgAbko+38yOIPiN4JmCRi5STXbvDsbMH3982JFIFep3eGWpaXil1Jy1a+FjH4N77oHTTgs7GqlQezu8UkSK6Tvfgb/8BSZNCjsSqVJK9CJhev11+PWv4StfgYMPDjsaqVJK9CJhuvlmGDYMrr027EikiinRi4Rl82a46y74538OJjETKRJNjScSlkMPhfnz4cgjw45EqpwSvUiYWlrCjkBqgFo3ImH4t3+Da64J5iwWKTIlepFSi0Tgxz8OhlRqhkopASV6kVK79Vbo7oYZM8KORGqEEr1UpciOCC3/1ULXzlxLOfWvri4outNvdXWDfMFXX4XbboMLL4T3v3+vYhPJlxK9VI3U5N66sJVn33iW1mdaB/VajY1BQo/HM+/Ptr1fl18O++0XjJ8XKRGNupGKF9kRYep9U2k6sIln33iW6U9MZ86aOcQ9TltHG//e8u80jmjM67UaGzMv5l0wv/pVMH7+sPS1e0SKRxW9VLzWha0s2rCI9pXtxD1O+6rgJ0DMY3lV9ckKvihJPhqFxx8Pfg046igNqZSSU6KXipPaoonsiNDW0YbjxDxY1CzmMaKxKADRWJS2jrY+vfpkYk/eilbF/+//BpOVfeITatdIaJTopSJk67+3LmwlFo/lfG5qVV/Uyj2VO1x3HZxxBtTXwyOPwDe/WeQ3FclMiV7KUvqomWRyn/7EdNo62oh7nFnLZ9G2vI3d8dwrVEZjUW57aHFBE/yQ/v7l3Hor3HILXHklLFsWJHyNmZeQ6GSslKXUqv1bp36rJ7m3r2qnbkgwtjEai0Ku3BmZDHcsL2hcDQ3Qlc+IzbPOgrfeCto1SvASMiV6CVVyxMycz8zpGRmT7LsnR838dfdfe51cjcWCVk2cOGSaQSDMBL9tGxxwABxzDHz3uwWNQWSw1LqR0ER2RDjhzhNYtGFRr5ExrQtbexJ7d7yb9pXtPSdX09XX1TOteRp82/fcCpjkGxqCdnteSX7XLvjoR+Hqqwv2/iKFoEQvJZXae58+fzqRnREc7xkZk6zmk4l9d3x3z2iaTJL990JJJvbkLa8En3T99fDSS0E/XqSMqHUjRZfankk9qXr3qrt7jkmOjHG8p5rP/aKFbc/k3ZrJ5qmn4Ic/hGnTlOil7JiX2TSpzc3NvnTp0rDDkEHI1G8HmPbwNO5YdgdfOO4LzFkzh13duxjCkKDHnmL40OFMHDWR1VtX93ntyY2TWX7lnsS+t+c39zqxp9q2DY4/HoYPhxdfDKY4ECkxM1vm7s2Z9qmil4JJHSkz81Mzgd4nVlNHzKQneQiq+pYJLbw07aWCxVTQhJ7N2rXB1a9z5yrJS1lSj14KIn2kTOr499QRM9lOqkLQb1/cuaffnn71avLW2M+0NYPusQ/WSSfBH/4AJ55YgjcTGTgleslbrql/0xN66zOtfU6s9ie9PZPt4qaiX9Warz/8AW66Cd5+O2jbiJQptW4kb5laM0CfhJ6cXyZ1/Hs26ck9Xw0NmRN+Q8OAX2rwfvrT4HbJJTBqVAnfWGRgVNFLRunVe7bWDPSu5pNiHuPhVx/OWM1PbpyM3+D4DT6oJA9BSya1RVPSVg3Ajh3BlMOf/SwcemiJ3lRkcFTRS0bp1Xum1kyyql/SuaRPQo/Goow7YBxbr9ta8thL4q67YPt2XRwlFUHDK6WPyI4IE38ykV3duxg+dDhLLlvCSb86iV3du3qOGT50OOuvXp/3gh6DkWsIZah/beNxOPZYGD0ann8+xEBE9sg1vFKtG+kjvXq/6P6LMrZmBrtMX76y9dtL2ofPZNs2OPpo+PrXQw5EJD9q3UgvmU6srtm6Bk+bPSx9KGQxlKzfPlCjR8O8eWFHIZI3JXrpJdOJ1WF1w7h8yuW9RtrUrM5O2L0bmprCjkQkb2rd1Lj00TXZTqwWu3qvGDfdBH/zN8GoG5EKoYq+BmWaZCw5imawwx1rwttvw+zZcMEFMHJk2NGI5E0VfQ3KtCxfpgW0Jc0vfwnvvKMhlVJxlOhrTPokY+lj4yWL7m742c+gpQU+9KGwoxEZECX6GpNtkrHktAWq6rNYvhwiEVXzUpHySvRmdrqZvWJm68xseob9E8zsSTNbaWYLzGxcyr7DzexxM3vZzNaY2RGFC18Gor9JxlTV5/DhD8MbbwSLfotUmH4TvZnVATOBM4BJwOfMbFLaYbcCs939eOBG4OaUfbOBW9z9A8CJwJ8KEbgMXKahk6k0uiaLaOI/xsZGqKsLNxaRQchn1M2JwDp3Xw9gZvcCZwNrUo6ZBHwjcf9p4H8Sx04Chrr7EwDuvrNAccsgZBo6CYOfQbJmXHEFbN0aXCS1t0tbiYQgn0R/GLAx5XEn8H/SjlkBnAf8GDgXGGlmBwFHA9vM7H6gCZgPTHfvvdqzmV0BXAFw+OGHD+JjSD6UzAdh82a45x64/HIlealYhToZey3QYmbLgRZgExAj+I/ko4n9HwYmAhenP9nd73T3ZndvHjt2bIFCknKQbZWogawYFapLLw2C1ElYqWD5JPpNwPiUx+MS23q4+2Z3P8/dpwDXJ7ZtI6j+O9x9vbt3E7R0/rYgkUtFyGc1qLJZMSrdggXw2GNw883BJGYiFSqfRP8CcJSZNZlZPTAVeDD1ADMbY2bJ15oBzEp57oFmlizT/47evX2pMukVfEWbOxfe976gRy9SwfpN9IlK/CrgMeBl4DfuvtrMbjSz5Fiz04BXzOxVoAG4KfHcGEHb5kkzWwUY8IuCfwopuWwtmbKtzgfjhz+E+fNhv/3CjkRkr2jhERmUQlfrZfXX8N134aWXgrHzIhVCC4/IoOQ6kVrVZs+Gk06C114LOxKRgtDsldJLY2Pp2y+hrxiV7s474bjj4Mgjw45EpCCU6KWXYib5smrPZLNyJbz4YjCBWdX/6iK1Qq0b6VHM8exlV7Vn88ADQYL/7GfDjkSkYJToa0D6KlLZFKOab2gIKvmyXf813bx5cMopcPDBYUciUjBq3VSx5EpSTQc29VpFKqkQ/fiKaMcMxBNPVND/SiL5UUVfRdIr99aFrSzasIj2le0ZV5Ha2yRfMe2YgTjwQDj22LCjECkoJfoqkrr+a3LueceJJeaQ25v55pMtmNRb1RW+X/0q/PrXYUchUnBK9FWiI9LBHcvu6KncZzw5g1i81yShRGNRbnu+DRvZNeABJVWX1NP96U/w85/DunVhRyJScEr0FS7Zrpl639SeRUW64920r2xnd3x3hmfE4NSBVfVV2aJJ98gjwa8pn/502JGIFJxOxla41oWtLNywsNe2zAk+YWgUxve/ilTVnWTtz7x5cOihMHly2JGIFJwSfQWL7Igwa/ms/g9MmNw4mY4v97/4SE1U8Kmi0WA64gsv1EVSUpWU6CtY68LWrAt9Q+YlAu3L2V+v5qr4pC1bggnMzjkn7EhEikKJvgJFdkQ4d865dHR14PTNzvV19Vw+5fJeY+Ylh/Hj4amnwo5CpGiU6CtI6gVQv9/0+6zHRWNRFndm7sM3NGQeP19z7Zokd9i2DUaNCjsSkaJRoq8AqQl+0YZFPPfGc1mP/eDYD/LStJey7q/6YZID9fLLwUyVc+fCueeGHY1IUSjRV4DkFa7PvfFcrwug0tXX1dMyoaXE0VW4hx6CeBxOPDHsSESKRuPoy1zyQqhcCT4pV8tGspg3D6ZMgcMOCzsSkaJRRV/GIjsinDLrlJ4LoTLRide98NZbsHgxfOtbYUciUlRK9GUo2ZPfv35/3ul+J+exquL3QrJtc+aZYUciUlRK9GUo2ZPPNHQyVaZx8jIAJ58Md90FzRnXUxapGurRhyjTgiCps05mM7lxMn6D90ny2RbzLubKURXtmGPgi1/U1bBS9ZToQ5Q6rXDqtvRZJ5Pq6+qZ1jwtYxWfaxGRUi/2XTHuvx82bAg7CpGiU6IPSfq0wl07u3qq+WyTkuXqxyuZD9Dzz8P558Nvfxt2JCJFpx59SD7/wOd7RtMkFwRxPOcIm9SefCGWAaxp994L9fVw5ZVhRyJSdEr0IeiIdLB66+qex9FYlLaONiaOmphxkrJMJ12V5PdCLAZz5sAnPwkjR4YdjUjRKdGXUHLY5JadfbN0zGO0TGjJOX2BFMgzzwRzQVx0UdiRiJSEEn0J5Ro2me94+MG0bGp2wrJsOjpg//2Dil6kBuhkbInkGjaZazRNusEkeU1kluYb34BNm+B97ws7EpGSUEVfIq0LW7OeaC3G1a01u4hIvkaMCDsCkZJRoi+BZDWfeqJ1+NDhrL96PY0j8r+aSRc+FcCaNfDVr8IPfgB/+7dhRyNSEmrdlECmaj45pDIfySteNdKmAFatggULYNiwsCMRKRkl+hJY0rmkz7DJgbRrlOALaO3a4H/No44KOxKRklHrpgRKPfGYRtnksHYtNDXBvvuGHYlIyaiiD0ExJh9raAhOwLprlE1Oa9fCsceGHYVISamiL6H+xsDvTYtGyT1PTU2allhqjhJ9CRWr165WzQDcf3/YEYiUnFo3FSBbIk+2a1TNi0gueSV6MzvdzF4xs3VmNj3D/glm9qSZrTSzBWY2LmVfzMw6ErcHCxl8ucu0sMhgdHXt6b+n3pTgB2jmTDj6aNi+PexIREqq30RvZnXATOAMYBLwOTOblHbYrcBsdz8euBG4OWXfu+4+OXE7q0Bxl7Vkgp/x5Iw+C4tIiNauhc5OzVgpNSefiv5EYJ27r3f3KHAvcHbaMZOApxL3n86wv6ZMnz+dhRsW0r6yvdfCIv1Rr73I1q8PTsZq6UCpMfkk+sOAjSmPOxPbUq0AzkvcPxcYaWYHJR7va2ZLzex5Mzsn0xuY2RWJY5Zu3bp1AOGXn8iOCHevuhsIrn5N/mx9plW99rC5a/y81KRCnYy9Fmgxs+VAC7AJSC58OsHdm4ELgf9nZu9Pf7K73+nuze7ePHbs2AKFVDqpvfjp86f3JPik5MIiHeu61GsPS3c3LFsG7+/z10+k6uWT6DcB41Mej0ts6+Hum939PHefAlyf2LYt8XNT4ud6YAEwZe/DDl9qck8u8j39iek91Xy6gcxtI0Xw5pvwxS8GUxSL1Jh8xtG/ABxlZk0ECX4qQXXew8zGAH929zgwA5iV2D4KeMfd30sccwrw/QLGX3LJVaKaDmzqSe5z1swh7nHaV7YTI5bxecWYilgGoLERbrkl7ChEQtFvonf3bjO7CngMqANmuftqM7sRWOruDwKnATebmQMLga8mnv4B4A4zixP89vBdd19ThM9RMslVop5747kgua9qp25IHUDWJD9m+Bi2XlfZ5x4qWiQCjz4KZ54JBx8cdjQiJZfXlbHu/gjwSNq2/0i5PxeYm+F5i4Hj9jLGspG6SlTqidZYLHOCT3pz/bheAz206lOJLV8Ol10GS5Yo0UtN0pWxWWS62Kl1YSuxeO6knlwW0G9w+Hbidkfv2Ss17XCJdXYGPw89NNw4REKiRJ9Fciz89PnBhcDJan53fHfO56kXX4Zefx3q6+Gw9FHBIrVBiT6D1LHw7Svbe0bWZFvzFWBy42Qafh5U8B1fXq5rcsrJunXBsMq6urAjEQmFZq/MIHUsfMxjTJ8/nRVbVvRZJQqCBJ9cWMS+XNIwJV/r1sGRR4YdhUholOjTpFbzSe0r2+n8RueAFvKWMjJ/PrzzTthRiIRGrZs0ma5sTVb1UqHGjoUJE8KOQiQ0SvRpHnr1oYzb71ryYEGW/ANNXlZSL78MN96o8axS02o60WcaQrnt7Sxfie85kTeY4ZGa2yYkS5bADTeodSM1raYTfXKOmomXtPZU6749S7m+U2V4RXrtNRg2DA4/POxIREJTE4k+vXJvbAQbGeHnS9qIe5x3j2mDEYkye/MJ4Gkv4MBmLShdkdatC+agH6pxB1K7aiLRJyv35OyRW7YALa1giXHxFoNTEzNLHv0wpI+Bt8R2qTxvvw1jxoQdhUioqj7RR3ZEmLV8FnGPM6tjVlDVj4jA5DYYmhgXPzQKUxJV/fbxmV9o+7jM21PkWlhEQjJkCMSzX+gmUguq/vfZ1oWtPdMWRGPRoKpv8T3VfFKyqk+blyaTbIlbJ1nL0Lx5WjpQal5VVvTJnvyKrhU91TxA3OPctmQWHL5wTzWfNDQK47PPUZNc7k+jZipMfX1wMlakhlVlRZ/syV90/0V9JyGri8KGFvj5S/2+jqYTrnD/8z9w993wi1/AgQeGHY1IaKquok/OMhn3OGu2ruk7EdmQOEx4JuvzNd69imzYAHPnqkcvNa/qEn3qLJOGMcTSPmJ3fVDRS/VLzlbZz8IwItWuqlo3yWo+OctknHjfMfH99OKlinj6H75IbaqqRJ9pzvj6unoun3I5t505s9/naxhkldmwAfbZB0aPDjsSkVBVVetmSeeSPnPGJ1d8yjXGXT35KjVsGHzkI1pwRGqeeZn9etvc3OxLly4NOwwRkYpiZsvcPeNcLVVV0YuISF9K9FKdXn8dPvQhePrpsCMRCV3NJPrGRnqmIk697e0iIlKmXn4ZVq6EffcNOxKR0NVMos+2WMhgFhGRCvDqq8HPo48ONw6RMlAziV5qzCuvwEEHBTeRGleViT5Tm0ZqzNtvB4uCi0h1XBQJfXMAAAWqSURBVDDV2KgWjKSZMgUOOCDsKETKQlUkeiV56WPGjLAjECkbVdm6ERGRPWo+0Wt+GxGpdlXRuhmIMpvxQUSk6Gqqolf1LiK1qOorelXwIlLraqqiFxGpRUr0IiJVToleRKTK5ZXozex0M3vFzNaZ2fQM+yeY2ZNmttLMFpjZuLT9+5tZp5n9rFCBp8q1epSISK3rN9GbWR0wEzgDmAR8zswmpR12KzDb3Y8HbgRuTtvfCizc+3Az6+rasxxg6k1LA4qI5FfRnwisc/f17h4F7gXOTjtmEvBU4v7TqfvN7ASgAXh878MVEZGByifRHwZsTHncmdiWagVwXuL+ucBIMzvIzIYAPwCuzfUGZnaFmS01s6Vbt27NL3IREclLoU7GXgu0mNlyoAXYBMSAacAj7t6Z68nufqe7N7t781hNLSsiUlD5XDC1CRif8nhcYlsPd99MoqI3sxHA+e6+zcxOBj5qZtOAEUC9me109z4ndEVEpDjySfQvAEeZWRNBgp8KXJh6gJmNAf7s7nFgBjALwN0vSjnmYqBZSV5EpLT6TfTu3m1mVwGPAXXALHdfbWY3Akvd/UHgNOBmM3OC0TVfHWxAy5Yte9PMNgz2+TmMAd4swuuWQqXGXqlxQ+XGrrhLr1xin5Bth3mNTAZjZkvdvTnsOAajUmOv1LihcmNX3KVXCbHrylgRkSqnRC8iUuVqKdHfGXYAe6FSY6/UuKFyY1fcpVf2sddMj15EpFbVUkUvIlKTlOhFRKpc1SZ6M/u2mW0ys47E7ZNZjss5BXOYzOwaM/PEBWmZ9sdSPt+DpY4vmzzi/pKZvZa4fanU8WViZq2JabY7zOxxMzs0y3Fl9Z0PIO6y+s7N7BYzW5uI/QEzOzDLcX80s1WJz7e01HFmMoDYyye3uHtV3oBvA9f2c0wd8DowEagnmJxtUtixJ2IbT3CR2gZgTJZjdoYd50DjBkYD6xM/RyXujyqDuPdPuf814PZK+M7zibscv3PgH4GhifvfA76X5bg/Zvv7X86xl1tuqdqKPk/5TMEclh8B1wGVdra8v7g/ATzh7n9297eBJ4DTSxVcNu6+PeXhflTI955n3GX3nbv74+7enXj4PMEcWhUhz9jLKrdUe6K/KvHr1SwzG5Vhfz5TMJecmZ0NbHL3Ff0cum9ieufnzeycUsSWS55xl+V3DmBmN5nZRuAi4D+yHFZW3znkFXfZfucJlwKPZtnnwONmtszMrihhTPnKFntZfef5TGpWtsxsPtCYYdf1wM8JVrbyxM8fEPyhlIV+Yv8mwa+H/Zng7pvMbCLwlJmtcvfXCxlnugLFHYpcsbv779z9euB6M5sBXAXckOHYsvrOBxB3yfUXd+KY64Fu4O4sL/N/E9/3wcATZrbW3Yu2Wl1SgWIvGxWd6N397/M5zsx+AczLsKvfKZiLJVvsZnYc0ASsMLNkTC+a2Ynu3mtxRHfflPi53swWAFMI+oLlHPcmgknwksYBC4oSbJp8/74Q/MN9hAwJs5y+8wyyxR3Kd95f3IkZbc8EPu6JxnaG10h+338yswcIWiJFT/QFiD203JJR2Cc2inUDDkm5/6/AvRmOGUpwYqqJPSdMPhh27Gkx/pHMJzVHAfsk7o8BXqNMTiT3E/do4A+J+Ecl7o8ug3iPSrn/L8DcSvjO84y77L5zgnMEa4CxOY7ZDxiZcn8xcHoZ/F3JJ/ayyi2hfmFF/sP4b2AVsBJ4MJn4gUMJVr1KHvdJ4FWCquz6sOPO8Dl6EibQDPwycf8jic+3IvHzsrBjzSfuxONLgXWJ2yVhx5qI6T7gpcTfl4eAwyrhO88n7nL8zhNxbAQ6ErfbE9t7/n0SjFhZkbitLpd/n/nEnnhcNrlFUyCIiFS5ah91IyJS85ToRUSqnBK9iEiVU6IXEalySvQiIlVOiV5EpMop0YuIVLn/D1UGgkWW2avaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Took  130.8554928302765  seconds\n"
     ]
    }
   ],
   "source": [
    "## Make the big plots\n",
    "\n",
    "start_time = time.time()\n",
    "make_rocs(reads_1kg)\n",
    "make_rocs(reads_hgsvc)\n",
    "print(\" Took \", time.time() - start_time, \" seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hgsvc\n",
      "Total wrong 17 Wrong stages single end: {'cluster': 4, 'none': 7, 'extend': 4, 'winner': 2}\n",
      "1kg\n",
      "Total wrong 3 Wrong stages single end: {'cluster': 1, 'none': 2}\n"
     ]
    }
   ],
   "source": [
    "#  Print some of the funky reads\n",
    "def print_stages(reads, map_q):\n",
    "    stages = {}\n",
    "    total_wrong = 0\n",
    "    for i, read in enumerate(reads.reads):\n",
    "        if not read.correct and proposed_cap(read) >= map_q: # and read.stage == \"cluster\":\n",
    "            #print(\"Read {} {}\".format(i, read))\n",
    "            if read.stage not in stages:\n",
    "                stages[read.stage] = 0\n",
    "            stages[read.stage] += 1\n",
    "            total_wrong += 1\n",
    "    print(\"Total wrong\", total_wrong, \"Wrong stages single end:\", stages)\n",
    "    \n",
    "print(\"hgsvc\")\n",
    "print_stages(reads_hgsvc_small, 60)\n",
    "\n",
    "print(\"1kg\")\n",
    "print_stages(reads_1kg_small, 60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}