{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "53e5fdf4-68ef-4084-a400-ff984db89cd3",
   "metadata": {},
   "source": [
    "# Multiplicity test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e520279-ba5f-4872-a38d-27487ed26c56",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-10 14:39:13.224762: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  SSE4.1 SSE4.2 AVX AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "from creme import creme\n",
    "from creme import utils\n",
    "from creme import custom_model\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import shuffle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3083d60-240d-4c05-afc2-02bded899c81",
   "metadata": {},
   "source": [
    "Load Enformer and example sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5e3979fd-a872-41fa-afc4-109ee830b141",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-06-10 14:39:17.638837: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  SSE4.1 SSE4.2 AVX AVX2 FMA\n",
      "To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.\n",
      "2024-06-10 14:39:18.081015: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1613] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 78954 MB memory:  -> device: 0, name: NVIDIA A100 80GB PCIe, pci bus id: 0000:85:00.0, compute capability: 8.0\n"
     ]
    }
   ],
   "source": [
    "data_dir = '../../../data'\n",
    "track_index = [5111]\n",
    "model = custom_model.Enformer(track_index=track_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "5b84bb66-1c29-448f-b78f-3ec690e39d22",
   "metadata": {},
   "outputs": [],
   "source": [
    "fasta_path = f'{data_dir}/GRCh38.primary_assembly.genome.fa'\n",
    "seq_parser = utils.SequenceParser(fasta_path)\n",
    "\n",
    "gene = 'KIT_chr4_54657941_+'\n",
    "\n",
    "gene_name, chrom, start, strand = gene.split('_')\n",
    "wt_seq = seq_parser.extract_seq_centered(chrom, int(start), strand, model.seq_length)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "61c69d44-b67f-4f71-b10b-5dc6c7b867f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TSS bin indeces\n",
    "bins = [447, 448]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "502ab1fc-bcb5-444d-ab7b-8351db4f7463",
   "metadata": {},
   "outputs": [],
   "source": [
    "wt = model.predict(wt_seq)[0,:,0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "34adb807-c057-4993-aa1a-dc2bf161d288",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: >"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABj0AAADFCAYAAAAPFjDeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxaklEQVR4nO3de5QU5Z3/8U/1fe4DA8OAghLjigbvF5jVxD3KSvxxcjSyiUlIlqg5ObsZXZBsIia/ePllDV5+uaxKMDEuZn/RmOhZ4upGE4IR1giIIBuMBiVRQZkZrjM909PXquf3x6TK7qZn6IGBrhner3PmwFTVVH+r+6mnnurv8zxlGWOMAAAAAAAAAAAARrhApQMAAAAAAAAAAAAYDiQ9AAAAAAAAAADAqEDSAwAAAAAAAAAAjAokPQAAAAAAAAAAwKhA0gMAAAAAAAAAAIwKJD0AAAAAAAAAAMCoQNIDAAAAAAAAAACMCr5LehhjFI/HZYypdCgAAAAAAAAAAGAE8V3So6enRw0NDerp6al0KAAAAAAAAAAAYATxXdIDAAAAAAAAAADgUJD0AAAAAAAAAAAAowJJDwAAAAAAAAAAMCqQ9AAAAAAAAAAAAKMCSQ8AAAAAAAAAADAqkPQAAAAAAAAAAACjAkkPAAAAAL6RSCRkWVbBTyKRqHRYAAAAAEYIkh4AAAAAAAAAAGBUIOkBAAAAAAAAAABGBZIeAAAAAAAAAABgVCDpAQAAAAAAAAAARgWSHgAAAAAAAAAAYFQg6QEAAADA17pSXZUOAQAAAMAIEap0AAAAAADgCofDuvXWW2Xbtp5/+3n15fq0O7Vbx+m4SocGAAAAYAQg6QEAAADANyKRiG677TZJ0uLfLNb+5H6Fwty2AAAAACgP01sBAAAA8KWMnZEkOcapcCQAAAAARgq6TAEAAADwDcdx9Prrr8txHLVva1dOOeXOzVU6LAAAAAAjBEkPAAAAAL6RTCY1ffr0gmULP7qwMsEAAAAAGHGY3goAAACArzG9FQAAAIBykfQAAAAA4GskPQAAAACUi6QHAAAAAF8j6QEAAACgXCQ9AAAAAPgaSQ8AAAAA5TqspMedd94py7K0cOFCb1kqlVJbW5uamppUW1uruXPnqrOz83DjBAAAAHCMIukBAAAAoFyHnPTYsGGDfvCDH+iMM84oWH7jjTfqqaee0uOPP67Vq1dr586duuqqqw47UAAAAADHJpIeAAAAAMoVOpQ/6u3t1bx58/Tggw/qX/7lX7zl3d3deuihh/Too4/qkksukSQtX75cp556qtatW6eZM2cOT9QAAAAARqVwOKx//ud/VmdPpza3b5YVtBQIMisvAAAAgPIc0t1DW1ub5syZo1mzZhUs37hxo7LZbMHyadOmacqUKVq7dm3JfaXTacXj8YIfAAAAAMemSCSie+65RwtvXajWL7Rq5jUzFQwHKx0WAAAAgBFiyCM9HnvsMW3atEkbNmw4YF1HR4cikYgaGxsLlk+YMEEdHR0l97dkyRLdfvvtQw0DAAAAwCiWsTPe/5neCgAAAEC5hjTSY8eOHVqwYIEeeeQRxWKxYQng5ptvVnd3t/ezY8eOYdkvAAAAgJHHcRy9/fbbeuvPbyneHldPZ49ydq7SYQEAAAAYIYY00mPjxo3atWuXzjnnHG+Zbdtas2aN7r//fv3qV79SJpNRV1dXwWiPzs5OtbS0lNxnNBpVNBo9tOgBAAAAjCrJZFJTp04tWPbJ//lkhaIBAAAAMNIMKelx6aWXasuWLQXLrrnmGk2bNk033XSTJk+erHA4rFWrVmnu3LmSpK1bt2r79u1qbW0dvqgBAAAAHDOY3goAAABAuYaU9Kirq9P06dMLltXU1Kipqclbft1112nRokUaO3as6uvrdcMNN6i1tVUzZ84cvqgBAAAAHDNIegAAAAAo15AfZH4w3/3udxUIBDR37lyl02nNnj1b3//+94f7ZQAAAAAcI0h6AAAAACiXZYwxlQ4iXzweV0NDg7q7u1VfX1/pcAAAAAAcRYlEQrW1tQXLfrrxp/rUOZ+qUEQAAAAARpJApQMAAAAAgMEw0gMAAABAuUh6AAAAAPA1kh4AAAAAyjXsz/QAAAAAgEMVCoX0pS99SZt3btb+1H5ZAUtW0Kp0WAAAAABGCJIeAAAAAHwjGo1q6dKluud392jbvm2SpGA4WOGoAAAAAIwUTG8FAAAAwHcydsb7P9NbAQAAACgXIz0AAAAA+IYxRnv27NG+PfvUl+iTZVmyHbvSYQEAAAAYIUh6AAAAAPCNvr4+NTc3Fyz78NoPVygaAAAAACMN01sBAAAA8DWmtwIAAABQLpIeAAAAAHyNpAcAAACAcpH0AAAAAOBrJD0AAAAAlIukBwAAAABfM8ZUOgQAAAAAIwRJDwAAAAC+Zhu70iEAAAAAGCFIegAAAADwNaa3AgAAAFCuUKUDAAAAAABXKBTSJz/zSW3auUmWLFmB/h8AAAAAKAdJDwAAAAC+EY1GdffSu/Wt//6WtywQZoA6AAAAgPJw9wAAAADAVzJ2puB3prcCAAAAUC5GegAAAADwDWOMuuJdyiT7Ex+WZclxSHoAAAAAKA9JDwAAAAC+0dfXpws+cEHBsjt+c0eFogEAAAAw0jC9FQAAAABfY3orAAAAAOUi6QEAAADA14wxlQ4BAAAAwAhB0gMAAACAr9nGrnQIAAAAAEYIkh4AAAAAfI3prQAAAACUi6QHAAAAAF9jeisAAAAA5SLpAQAAAMDXmN4KAAAAQLlClQ4AAAAAAFzBYFAz/naGOno6VBepUyKXoKsWAAAAgLKR9AAAAADgG7FYTDf83xu05p01Om38aXpt92sKhoOVDgsAAADACEGfKQAAAAC+4k5nFQ6GC34HAAAAgIMh6QEAAADAV2znL0mPQH/SwzFOJcMBAAAAMIIwvRUAAAAA30gkErr2nGslST/UDyVJC59cWMGIAAAAAIwkjPQAAAAA4GuOGOkBAAAAoDwkPQAAAAD4GtNbAQAAACgXSQ8AAAAAvmaMqXQIAAAAAEYIkh4AAAAAfM02dqVDAAAAADBCkPQAAAAA4GtMbwUAAACgXCQ9AAAAAPga01sBAAAAKFeo0gEAAAAAgCsYDOqUmacokUloUt0ktSfa5ViM9AAAAABQHpIeAAAAAHwjFovps3d9Vju6d+ias6/R8leWS+of7WFZVoWjAwAAAOB3TG8FAAAAwFdsp//B5aHA+320jJjiCgAAAMDBkfQAAAAA4Cu26U96hANhbxkPMwcAAABQDqa3AgAAAOAbiURCd/6vO+U4jv49+O/KOTl97v99jqQHAAAAgLKQ9AAAAADgK9lUVpJky/aWkfQAAAAAUA6mtwIAAADgeyQ9AAAAAJSDpAcAAAAA3yPpAQAAAKAcJD0AAAAA+B5JDwAAAADlIOkBAAAAwPdIegAAAAAoB0kPAAAAAL5H0gMAAABAOUKVDgAAAAAAXJZlaeL0iTKO0dQxU/VO/B1ZlkXSAwAAAEBZSHoAAAAA8I1oLKqPfetjkqTvzP6Ovrbqa0rlUiQ9AAAAAJSF6a0AAAAA+IZtbO//wUBQAav/loWkBwAAAIByDCnpsWTJEp1//vmqq6tTc3OzrrzySm3durVgm1Qqpba2NjU1Nam2tlZz585VZ2fnsAYNAAAAYHSynbykh0XSAwAAAMDQDCnpsXr1arW1tWndunVauXKlstmsLrvsMiUSCW+bG2+8UU899ZQef/xxrV69Wjt37tRVV1017IEDAAAAGH3ivXH9+2f/XT+e92NNnDBRy65epmwqS9IDAAAAQFmG9EyPZ599tuD3hx9+WM3Nzdq4caM+8pGPqLu7Ww899JAeffRRXXLJJZKk5cuX69RTT9W6des0c+bM4YscAAAAwKhjO7ZS8ZQkKa20t5ykBwAAAIByHNYzPbq7uyVJY8eOlSRt3LhR2WxWs2bN8raZNm2apkyZorVr15bcRzqdVjweL/gBAAAAcGzKf6ZHPpIeAAAAAMpxyEkPx3G0cOFCXXjhhZo+fbokqaOjQ5FIRI2NjQXbTpgwQR0dHSX3s2TJEjU0NHg/kydPPtSQAAAAAIxw+c/0yEfSAwAAAEA5Djnp0dbWpldffVWPPfbYYQVw8803q7u72/vZsWPHYe0PAAAAwMg1UHKDpAcAAACAcgzpmR6u66+/Xk8//bTWrFmj448/3lve0tKiTCajrq6ugtEenZ2damlpKbmvaDSqaDR6KGEAAAAAGGUY6QEAAADgcAxppIcxRtdff71WrFih5557TlOnTi1Yf+655yocDmvVqlXesq1bt2r79u1qbW0dnogBAAAAjFo80wMAAADA4RjSSI+2tjY9+uijevLJJ1VXV+c9p6OhoUFVVVVqaGjQddddp0WLFmns2LGqr6/XDTfcoNbWVs2cOfOIHAAAAACA0cPIaPwHxyuggI6rP04diQ5ZlkXSAwAAAEBZhpT0WLZsmSTpb/7mbwqWL1++XJ///OclSd/97ncVCAQ0d+5cpdNpzZ49W9///veHJVgAAAAAo1s4FtbHv/NxNdc065uXfFN3rLlD27u3k/QAAAAAUJYhJT2MMQfdJhaLaenSpVq6dOkhBwUAAADg2OQ+0yMYCEqSAlb/jLwkPQAAAACUY0jP9AAAAACAI8l9pkfQIukBAAAAYOiGNNIDAAAAAI6k3t5ePfqFRxW0gvpR9EfqyfTo4/d9nKQHAAAAgLKQ9AAAAADgG7ax1burV5LUrW5J/dPskvQAAAAAUA6mtwIAAADgG+4zPYqR9AAAAABQDpIeAAAAAHzDfaZHMZIeAAAAAMpB0gMAAACAbzDSAwAAAMDhIOkBAAAAwDcY6QEAAADgcJD0AAAAAOAbjPQAAAAAcDhClQ4AAAAAAFyOHI2ZPEbRYFRjqsaoK9Uly7JIegAAAAAoC0kPAAAAAL4Rjob1iaWf0AXHXaDrzrlOD258UC/vfJmkBwAAAICyML0VAAAAAN9wn+kRDAQlSQGr/5aFpAcAAACAcpD0AAAAAOAb7jM9ghZJDwAAAABDx/RWAAAAAHwj2ZfU422P66nQU1oSW6LudLc+etdHSXoAAAAAKAtJDwAAAAC+kXNy2r9jvySpXe2SJGMMSQ8AAAAAZWF6KwAAAAC+4U5vVYykBwAAAIBykPQAAAAA4Bvug8yLkfQAAAAAUA6SHgAAAAB8g5EeAAAAAA4HSQ8AAAAAvsFIDwAAAACHg6QHAAAAAN8YKLlB0gMAAABAOUKVDgAAAAAAXLaxVdtcq+pgtWLhmJK5pCzLIukBAAAAoCyM9AAAAADgG6FoSJ/50Wf009/9VO+8844eXPWgQtEQSQ8AAAAAZSHpAQAAAMA33Gd6BK2gJClg9d+ykPQAAAAAUA6SHgAAAAB8w3b+kvQIkPQAAAAAMHQ80wMAAACAbySTSa1YtELPx55XbaRWPZkeffjWD5P0AAAAAFAWkh4AAAAAfCNn57R7227t1m5v2UXmIpIeAAAAAMrC9FYAAAAAfGOg5AZJDwAAAADlIOkBAAAAwDfcZ3oUI+kBAAAAoBwkPQAAAAD4hm1IegAAAAA4dCQ9AAAAAPgGIz0AAAAAHA6SHgAAAAB8g5EeAAAAAA5HqNIBAAAAAIDLdmzF6mOKhWIKWkEvCULSAwAAAEA5GOkBAAAAwDeCsaD+/id/r03bNmnPnj365eZfKhwLk/QAAAAAUBaSHgAAAAB8w32mRzAQlCQFrP5bFpIeAAAAAMpB0gMAAACAb7jTWQUtkh4AAAAAho5negAAAADwjXQyraf+91PaXLdZ0VBUiVxCZy06i6QHAAAAgLKQ9AAAAADgGzknp/ZX29Wudm/ZmeZMkh4AAAAAysL0VgAAAAB8w32mRzFjzFGOBAAAAMBIRNIDAAAAgG+4z/QodzkAAAAA5CPpAQAAAMAXHOMMOKKD6a0AAAAAlIOkBwAAAABfGGhqK4nprQAAAACUh6QHAAAAAF8YbDQH01sBAAAAKEeo0gEAAAAAgPR+YiMUDSkcCMuyLBn1j/BgeisAAAAA5WCkBwAAAABfsB1b4VhY1z5+rRKJhBKJhLbs2KJwLMz0VgAASPrvd/5b//bKvynn5CodCgD4FkkPAAAAAL7gjvQIWAFZluX9P38dAADHKmOMVvxxhda/u16v7X6t0uEAgG+R9AAAAADgC+6DzIOBoLfMTXoMNL3Vn/b9SVs6tyiejh/5AAEAqKCeTI8SmYQk6b34exWOBgD8i2d6AAAAAPAF29jKZXL69V2/1rZ7t8myLGWcjE784okysQOnt9rZs1P3vHiPN/XV3532d/rbk/72aIcNAMBRkZ/oeK+HpAcADISkBwAAAABfsB1bxjF6e8Pbeltve8uv+cI1Jae32rhzo4wxCgfDytpZPbPtGV184sWKBCNHMWoAAI6O/EQHIz0AYGBMbwUAAADAFwZ7bkep6a02tW+SJH3m9M+oqbpJiUxC699df8TiAwCgkvITHR29HTzMHAAGQNIDAAAAgC+4z/QoxZ3CytXZ26mdPTsVsAI6q+UsXTL1EknSqrdWHbAtAACjwc6end7/HeOoo7ejgtEAgH+R9AAA4AjLOTm9uONFHrILAAcx2EiP4nXuKI9Tx5+q6nC1LppykWKhmNp72vX7zt8f0TgBADjajDFe0qM+Wi+JKa4AYCBHLOmxdOlSnXjiiYrFYpoxY4ZeeumlI/VSAAD42hOvPaEfb/6x7l1/L0PQAWAQg430yJ/eqi/bp7XvrpUknd1ytiQpForp4hMvliT97A8/UzqXPoKRAgBwdO3p26OMnVEoENKZLWdK4mHmAA7PaB4dfUQeZP6zn/1MixYt0gMPPKAZM2boe9/7nmbPnq2tW7equbn5SLzkqNSX7dO78XcVC8U0vnq8qsJVR+R1bMfWnr49aow1KhqKHpHXONJsx5ZjHIUCIb3d9bb+uOePOrnpZJ005iRZljUsr/Fu/F09uuVR2Y6tT5/+aZ3YeOKw7HckcowjS9awvbfGGO1K7NKO+A6lc2mdPuF0r+cKSss5OQWt4LB9BhicYxwFrEPrJ7Bt3zY9//bzkqQd3Tv0X2/8l66YdsUwRgeX7djqy/apJlJzyJ9XJbl1oSQ11zRzfuOIcW9u3DJ2OHXccEvbAycqHOMolUups7dTj2x5RJ29naqJ1Oicied428w5eY5e3vmy9vbt1ROvPaGPnfIx1UXqvPV92T519HYoFAhpcsNk3xz30WCMUc7JKRwMD7pdOpdWe2+7xlWPU22k9ihFd/SN1LaUMUaOceQYR7axZTu2oqGoQoGQbMdWMpdUTbjmkI7LLSMZO+Pt82hJ5VKS+pOXh6u4jstnO7aMzFE9tqNtsOM/2owxB8ThJrCPVP2bzqWVc3KqidQckf3nv45tbFWHq4/o6+RzExwT6yZqSsOU/mWDjPQwxqg73a36aL1vr3e2Y+u9nve0K7FLvZlehQNhNcQaNLVx6hH/DAdijFFftk97+vZod99u7e3bq4ZYg6Y3Tx/yddG9pu5L7pMkTRs37aiWmSOpK9WlTe2btCuxy/se56yWswYsa45xtLNnpyLBiMZXj69oHZWxM9rbt1c1kRrVRmqH5fxIZBJKZBMaVz3O258xRslcUrFQzFuWc3Lq7O1UY6yxYmV8b99evfTeS1r/3np1pbo099S5umjKRb64bgwnyxyBlM6MGTN0/vnn6/7775ckOY6jyZMn64YbbtDixYsLtk2n00qn37+5icfjmjx5sr7+y68rVlPY4BksVKPCdeUcVsAKyLKsw/rydigxlSvn5JTMJtWV6ipYXhupVUOsQQEr0B+7LPVmetWX7VM4GFYkGFE4EC5oqEaCERljZJv+pICbHHAbGqFASF2pLuWcnAJWQBNqJyhoBUsehzHG+7/7926j290uf7mRKfjXbZhXh6u9C0X+++dud7DXzTk5WbIUDoYVDoRlG1u7ErtkO7ZqIjVKZBLePpqqm2Sp/7ONBCOyLMs7fmOM9166ZcF77RKf686enQUNtEl1k4ZUbtw48qXttHrSPbKNraAVVDAQLPlvwAoc8P4Uc4yjjJ2RkTlgHwEr4P1tqfd1oM+h1LbJbFLxdFzRUFRNVU0KBoLe8RW/H+4x558L+fty497Tt6egN6ZlWWqINqgn06NYKFbQSCt+Hwd6zYGWF/+dJUtGRqlcyrvZOtjxFP99ObH0ZfvUk+7xPp/6aP2gicz8vy2OpTvVra5Ul6rCVWqqahrSBbr4s8hfPlj5cuPI/13qr69sxy6Is3i74s9koPO8VIxSf0O4N9OrnJM7aHLWrSNsx/a+CAhYAVWFqw4oA6FASLaxvbLn1gNu/RqwAnKMo/2p/UrlUmqMNao2Uvv+Fw1/qU9zTu6AG/X813bXT26YrB3dO2RZlsZXj++vn/PqSUkF7+FA/2adrJLZZH/ysUQ5zH+vS60r/myKPyN32UDXDsc4igQjBe9pqc+sJ9Mj27FVF61TNDh4Qn0o18t0Lq1ENqFYKKbqcLUCVkAZO6N4Ou7V/+FgWA3RBvVl+yT1D/93v9woVdZLvX5xPTiQ/PPzYGXf3W/aTst2bAUDQYUCIe9Lt3g6ru5UtySpIdaghmhDyRgGOlcHOr8H2sYt58FAUJas/nPH2Mo5uYJzO/9aGbAC3nXJva67ifD84zlUxdcdY/rr55yTU1W4SkEreNDrlKSC8zj/vC7VNvHejwHOp6HEXbBsCE3soZwDpfY7lL/P2Bl1p7oVsAJqiPWfJ8lsUnXROtVF6g64lpVzHIfS5s1/n/O/wN2X3KdMMqPln1xesP01P79G4Vjhl/W1kVrd2Hqjjq8/vmD5q7te1X3r7ztoDNFQVDXhmpJtwfw4B6on848j/3zIrzMty1IoEFI4EPbaTIMZ7NzN2lklc0lFghFVh6tLXvNioZh3ncw/T3JOTnv79npfBNaES99kO8bRvuQ+79xormk+aJKkOO5DvQc6mEM5NwfSm+lVPB1XJBjR2KqxB21LlXtdLZa1s0rlUl45KFVHlmoLuXWx7eTVyX/5jPPrrfzXrg5Xqy/bJ2P6v9CvjdQW1H+lyq57vU/baaVzaaVyqYL9h4Nh1YRrFA1FC+IrbkMMdi9Vqs52tzEyytgZJbNJry01pmqMwoGwMnZGkkrer7n7K3Uf6t5LuyNrLctS0AqqLlqncCCsPX17ZGTUGGtUOBAuuC6494jlLAtYAdVGahUJRrxt8ttKtrEVDoS9tkrx53ywa9lA22XsjDJ2RrWRWq8ecH/ca6R7rtdH6xUNRQu+Pxis3BV/vpK8a3BxG9G7p857j9zjzl8v9de17pe8GTtT0F6LhWKKBqMKB8PK2ln1ZHpkyVI0FFU02J94yzk57z132ywDtbPSubR29+2WMUZ10bqCL6iLt89vxxZ/xvl/U+oakMqltDe5V8YY1URqVBV6/75uoPum/P0NFNNA37m4sSVzSSUyCc08fqY+csJHdPfv7pZlWWqqalIkGDng9ff07VEql1I4GNaY2Bj1ZnrlGEc1kZoDEn/l3NeWcqjbum287nS3sna25LYTaiYUfBeS334rdZ5K/fd44WBYoUBIGTujvmyfqkJVioViSuVSso3t1cf570F+vbI/tV/JbLLkMcRCMVWF+/cXCoQGrDOk/vpof2p/wXkWsALedaf43nOw93Uwg7UNy203DnYdH2gfe/r2HLCuNlKraChaUNbdbXoyPd69d3203vvCf6DXdu9HJXnXz2AgWFAG3HaX+/9S3/W57637N7Zjqzvd7b1uwAqoPlqvWChWsj7I/7dUjFL/Z+3ee4YCIdVF6xSwAoqn48raWVmWpapQlWoiNepKdXllvqm6acBzcaB2acHvg9y/FJ9v7veFASugzt7OA7YfWzW2ZJvvUBNCw9luy3fhlAs16wOzyothuJMemUxG1dXVeuKJJ3TllVd6y+fPn6+uri49+eSTBdvfdtttuv322w/Yz+cf+7wi1ZHhDG1EGlc9zvti/EhyGxOjQSgQ0injTtEbe98oefE8HGdPPFuhQEgb3tswrPtFf6P3+PrjZYzR211vVzocYFg1xhp169/cqhWvr9Cad9ZUOhz4WCgQkmVZw379AkaSiXUTdcUpV+jsif3TVhlj9O2139abe9+UJNVEajS1car+7rS/08S6iSX38Zs//0YvbH9BHb0dB9yQj6kao3Qu7d2c4kDFHYkAAP4x/6z5On/S+fo/q/+PN0p4JKsKV2lS3STVR+uVtbPam9yr9p72SoelhliDxlWPU1NVk3b27NS78XcPaT/10XqNrxmvvmyfL45rOJ009iSdPPZk2cbWiztePGjbIRaKecn8SouFYkrb6WGd3slN4B5MNBSt6DSslmXplKZTNOP4GerN9OrJPz7pi8+kHLM/OFtXnXpVWdsOe9Jj586dOu644/Tiiy+qtbXVW/7Vr35Vq1ev1vr16wu2H2ikx7pt61Rb15+VL6dnTfF2B6wboFdDfqZ4OIbxDEcmKxgIqjpcraaqJi/zmcqltDux28vMuxnKmnCNqsPV3uiOrJNV0AoqEox4vUDcHhHFPTSl/p4WddE6ja8er65Ulzp6Ow48phI9E8r9V3o/K+j2ZO3L9nkV4UC9IAbqhRGwAv3Z9L/0dnNPyuaaZkWCEe1L7tO46nHejdp7Pe95WVO3QnHfCymvZ2HR/NGlylxNpEaT6yfLsiy9G3+35AOJh3o6hYNh1UfrFbSCXu+t4n/dHjOletEXx+z27ijeT8F2Q+yVXLx9NBRVY6xRyWyyoCdgqZ5bA71u/mtblqWxVWPVXNPsfS57+/aqN9Or2kit0nZa8XR80Ox/vsF6JQy2zu0ZWarn+mAZ/fzlg8XojlgJBoLK2lnF0/GCkSXFSvUkcN/r2kitxlePV0+mR/uT+4fco7Lcz32g9yD/OEOBkNfbIr+XXf7f5P/dUEfISPJ60gWsgLpT3V7Pv1Lc+i6/J4jt2Erb6QN6RrrTWuT3iC3usSb1JyyqwlXal9ynZDbp1SFunRYKhPp7hfylp6Tb4z1gBbxRKhPrJnpTtu1K7FI8HVc6ly4YcZj/Pg30b/4oC3cEVzk9cdx9F29T3KMtv4fTQNcOy7K8XpmDfQ510ToFraDi6fiwNqDc3oKpXKqgN019tF51kTpVh6u1L7lP8XRc1eFqGRlvRN1AZb14WfE67/eDjFYarHdmvlgo5pVNtzevMUaRYMSbOvGd7ndKNoQPNhphoLZMqXPOLedu7/T8Xsj553Z+L9L8Xpxu71m37s4fJXKoSvVCjoViClpBJXNJb1SJexwD9Vou7nma3zO3VFvFVfzZDqVdN5R25MF6Yx3K35cbUygQUkO0QUamf9RgqEq1kVp1pbrUm+ktK47Djb/4+pZfxzTGGjWmakzJv3OMo6yd9UbvliOdSyvrZL24Q4GQ12u9o7fDuz7k1/sFsZbocVvw+1/W54+uLR4RlXNyyjrZQZ9Xkm+g8zwSjCgainp1sNsuzh8ZnLEzXvsiv5y7PUvdOrK4DZL/WmOqxmhs1Vj1pHuGPFf8YHXq4RjOW1Yjo6pQlcZWjVUylxxSW2qoI+3cnuzGGK8c5NdjruL3y62HvXr5L20bt2zlt0UCVkB92T51p7tVF6lTTaRG3alu9WX7SrZrSvVmjwajioaiXq/7SDCitJ327tvctldxb/SBesEOdo9XfA2OBCOKhWJqjDXKNrbae9rlGMc7T/N77+a/d94ocOv9HtJuG6wqVOXdF7n3jfF0XBk7o+aaZgUDQe1L7is5Wrl4NoiBltnm/ZHI+du415lgIOj1/C117RrofTnYdu60Yz3pHiVzSa88uDEZY9RU3aRwIOz1oM///uBgbfDi8uGOtsgvO/nHmt9LvXhUjrsumUuqL9vnfd7uKCR3pL07TZQ7Ikfqr7vTdv80VaFASJasgl7a+fG4jPpHOU2qm6RoMOrV8aXaaAMdT/G5OdA1IBQIqaW2RZFgRLv7ditrZ0u2A0v97i4rFVNxOSu1LBaKqaW2RZbV/57sT+7X/tR+r/2Vf50aUzVGE2omaG9yr7pSXaqL1Hl1hlteio81P8ZSy0spp44utZ9YKObd2xaXx65Ul9p72gc8H4vbdO6IENuxlXWyXnuhKlylZDaptJ322pRuezV/RJj0/rnYEG1QU/X7o2dciUxCvZler+xmnWxBHKXameOqxxVM3b23b683yqBUHV3qvRro/T3Yd6WD3d+Uc+9zsNdojDWqqbrJ+z1jZwpmSMm/JkhSVahKLbUtyjk57YjvOOCeodTruyMPsnbWu9fIv992P3u3HnTLR/51w32v87cbWzVWddE6OcZRPB33Rl8M9H3mQPcP7j1RwAp4I672Jfd5393WRevUEG1QKpdSb6ZXiWxCtZFaTaiZ4E25WuqesdT5MtA5Wc5IH/f7Qkv99/IT6yaqMdbobR9Pxw9IoBZf24diqN9RDYX7/WE5Kp70KBaPx9XQ0KDu7m7V1zOnPwAAAAAAAAAAKM+wP8lo3LhxCgaD6uwsnB+ss7NTLS0tw/1yAAAAAAAAAAAAko5A0iMSiejcc8/VqlWrvGWO42jVqlUFIz8AAAAAAAAAAACGU+jgmwzdokWLNH/+fJ133nm64IIL9L3vfU+JRELXXHPNkXg5AAAAAAAAAACAI5P0uPrqq7V7927dcsst6ujo0FlnnaVnn31WEyZMOOjfuo8YiccPfEg0AAAAAAAAAAA4NtXV1RU8VL6UYX+Q+eH685//rJNOOqnSYQAAAAAAAAAAAB/p7u5WfX39oNsckZEeh2Ps2LGSpO3bt6uhoaHC0QAHF4/HNXnyZO3YseOgJxzgB5RZjDSUWYw0lFmMNJRZjDSUWYw0lFmMNJRZ+FldXd1Bt/Fd0iMQ6H+2ekNDAycVRpT6+nrKLEYUyixGGsosRhrKLEYayixGGsosRhrKLEYayixGqkClAwAAAAAAAAAAABgOJD0AAAAAAAAAAMCo4LukRzQa1a233qpoNFrpUICyUGYx0lBmMdJQZjHSUGYx0lBmMdJQZjHSUGYx0lBmMdJZxhhT6SAAAAAAAAAAAAAOl+9GegAAAAAAAAAAABwKkh4AAAAAAAAAAGBUIOkBAAAAAAAAAABGBZIeAAAAAAAAAABgVCDpAQAAAAAAAAAARgXfJT2WLl2qE088UbFYTDNmzNBLL71U6ZBwjFqzZo0+9rGPadKkSbIsS7/4xS8K1htjdMstt2jixImqqqrSrFmz9OabbxZss2/fPs2bN0/19fVqbGzUddddp97e3qN4FDhWLFmyROeff77q6urU3NysK6+8Ulu3bi3YJpVKqa2tTU1NTaqtrdXcuXPV2dlZsM327ds1Z84cVVdXq7m5WV/5yleUy+WO5qHgGLFs2TKdccYZqq+vV319vVpbW/XMM8946ymv8Ls777xTlmVp4cKF3jLKLfzktttuk2VZBT/Tpk3z1lNe4UfvvfeePvvZz6qpqUlVVVU6/fTT9fLLL3vruQeDn5x44okH1LOWZamtrU0S9Sz8x7ZtfeMb39DUqVNVVVWlk046Sd/85jdljPG2oZ7FaOGrpMfPfvYzLVq0SLfeeqs2bdqkM888U7Nnz9auXbsqHRqOQYlEQmeeeaaWLl1acv3dd9+te++9Vw888IDWr1+vmpoazZ49W6lUyttm3rx5+sMf/qCVK1fq6aef1po1a/TFL37xaB0CjiGrV69WW1ub1q1bp5UrVyqbzeqyyy5TIpHwtrnxxhv11FNP6fHHH9fq1au1c+dOXXXVVd5627Y1Z84cZTIZvfjii/rxj3+shx9+WLfcckslDgmj3PHHH68777xTGzdu1Msvv6xLLrlEV1xxhf7whz9IorzC3zZs2KAf/OAHOuOMMwqWU27hNx/60IfU3t7u/bzwwgveOsor/Gb//v268MILFQ6H9cwzz+i1117Tt7/9bY0ZM8bbhnsw+MmGDRsK6tiVK1dKkj7xiU9Iop6F/9x1111atmyZ7r//fr3++uu66667dPfdd+u+++7ztqGexahhfOSCCy4wbW1t3u+2bZtJkyaZJUuWVDAqwBhJZsWKFd7vjuOYlpYWc88993jLurq6TDQaNT/96U+NMca89tprRpLZsGGDt80zzzxjLMsy77333lGLHcemXbt2GUlm9erVxpj+8hkOh83jjz/ubfP6668bSWbt2rXGGGN++ctfmkAgYDo6Orxtli1bZurr6006nT66B4Bj0pgxY8yPfvQjyit8raenx5x88slm5cqV5uKLLzYLFiwwxlDPwn9uvfVWc+aZZ5ZcR3mFH910003moosuGnA992DwuwULFpiTTjrJOI5DPQtfmjNnjrn22msLll111VVm3rx5xhjqWYwuvhnpkclktHHjRs2aNctbFggENGvWLK1du7aCkQEHeuutt9TR0VFQXhsaGjRjxgyvvK5du1aNjY0677zzvG1mzZqlQCCg9evXH/WYcWzp7u6WJI0dO1aStHHjRmWz2YIyO23aNE2ZMqWgzJ5++umaMGGCt83s2bMVj8e93vfAkWDbth577DElEgm1trZSXuFrbW1tmjNnTkH5lKhn4U9vvvmmJk2apA984AOaN2+etm/fLonyCn/6z//8T5133nn6xCc+oebmZp199tl68MEHvfXcg8HPMpmMfvKTn+jaa6+VZVnUs/Clv/7rv9aqVav0xhtvSJL+53/+Ry+88IIuv/xySdSzGF1ClQ7AtWfPHtm2XVDZS9KECRP0xz/+sUJRAaV1dHRIUsny6q7r6OhQc3NzwfpQKKSxY8d62wBHguM4WrhwoS688EJNnz5dUn95jEQiamxsLNi2uMyWKtPuOmC4bdmyRa2trUqlUqqtrdWKFSt02mmnafPmzZRX+NJjjz2mTZs2acOGDQeso56F38yYMUMPP/ywTjnlFLW3t+v222/Xhz/8Yb366quUV/jSn//8Zy1btkyLFi3S1772NW3YsEH/9E//pEgkovnz53MPBl/7xS9+oa6uLn3+85+XRLsA/rR48WLF43FNmzZNwWBQtm3rjjvu0Lx58yTxXRdGF98kPQAAw6OtrU2vvvpqwbzdgB+dcsop2rx5s7q7u/XEE09o/vz5Wr16daXDAkrasWOHFixYoJUrVyoWi1U6HOCg3F6bknTGGWdoxowZOuGEE/Tzn/9cVVVVFYwMKM1xHJ133nn61re+JUk6++yz9eqrr+qBBx7Q/PnzKxwdMLiHHnpIl19+uSZNmlTpUIAB/fznP9cjjzyiRx99VB/60Ie0efNmLVy4UJMmTaKexajjm+mtxo0bp2AwqM7OzoLlnZ2damlpqVBUQGlumRysvLa0tGjXrl0F63O5nPbt20eZxhFz/fXX6+mnn9Zvf/tbHX/88d7ylpYWZTIZdXV1FWxfXGZLlWl3HTDcIpGIPvjBD+rcc8/VkiVLdOaZZ+pf//VfKa/wpY0bN2rXrl0655xzFAqFFAqFtHr1at17770KhUKaMGEC5Ra+1tjYqL/6q7/Stm3bqGfhSxMnTtRpp51WsOzUU0/1pmXjHgx+9c477+g3v/mNvvCFL3jLqGfhR1/5yle0ePFifepTn9Lpp5+uz33uc7rxxhu1ZMkSSdSzGF18k/SIRCI699xztWrVKm+Z4zhatWqVWltbKxgZcKCpU6eqpaWloLzG43GtX7/eK6+tra3q6urSxo0bvW2ee+45OY6jGTNmHPWYMboZY3T99ddrxYoVeu655zR16tSC9eeee67C4XBBmd26dau2b99eUGa3bNlS0IBZuXKl6uvrD7gBBY4Ex3GUTqcpr/ClSy+9VFu2bNHmzZu9n/POO0/z5s3z/k+5hZ/19vbqT3/6kyZOnEg9C1+68MILtXXr1oJlb7zxhk444QRJ3IPBv5YvX67m5mbNmTPHW0Y9Cz/q6+tTIFD4VXAwGJTjOJKoZzHKVPpJ6vkee+wxE41GzcMPP2xee+0188UvftE0Njaajo6OSoeGY1BPT4955ZVXzCuvvGIkme985zvmlVdeMe+8844xxpg777zTNDY2mieffNL8/ve/N1dccYWZOnWqSSaT3j4++tGPmrPPPtusX7/evPDCC+bkk082n/70pyt1SBjF/vEf/9E0NDSY559/3rS3t3s/fX193jb/8A//YKZMmWKee+458/LLL5vW1lbT2trqrc/lcmb69OnmsssuM5s3bzbPPvusGT9+vLn55psrcUgY5RYvXmxWr15t3nrrLfP73//eLF682FiWZX79618bYyivGBkuvvhis2DBAu93yi385Mtf/rJ5/vnnzVtvvWV+97vfmVmzZplx48aZXbt2GWMor/Cfl156yYRCIXPHHXeYN9980zzyyCOmurra/OQnP/G24R4MfmPbtpkyZYq56aabDlhHPQu/mT9/vjnuuOPM008/bd566y3zH//xH2bcuHHmq1/9qrcN9SxGC18lPYwx5r777jNTpkwxkUjEXHDBBWbdunWVDgnHqN/+9rdG0gE/8+fPN8YY4ziO+cY3vmEmTJhgotGoufTSS83WrVsL9rF3717z6U9/2tTW1pr6+npzzTXXmJ6engocDUa7UmVVklm+fLm3TTKZNF/60pfMmDFjTHV1tfn4xz9u2tvbC/bz9ttvm8svv9xUVVWZcePGmS9/+csmm80e5aPBseDaa681J5xwgolEImb8+PHm0ksv9RIexlBeMTIUJz0ot/CTq6++2kycONFEIhFz3HHHmauvvtps27bNW095hR899dRTZvr06SYajZpp06aZH/7whwXruQeD3/zqV78ykg4oh8ZQz8J/4vG4WbBggZkyZYqJxWLmAx/4gPn6179u0um0tw31LEYLyxhjKjLEBAAAAAAAAAAAYBj55pkeAAAAAAAAAAAAh4OkBwAAAAAAAAAAGBVIegAAAAAAAAAAgFGBpAcAAAAAAAAAABgVSHoAAAAAAAAAAIBRgaQHAAAAAAAAAAAYFUh6AAAAAAAAAACAUYGkBwAAAAAAAAAAGBVIegAAAAAAAAAAgFGBpAcAAAAAAAAAABgVSHoAAAAAAAAAAIBR4f8D6b8c4kz+k6cAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 2000x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "utils.plot_track([wt], color='green', zoom=[0, 896], marks=bins)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "727193db-efab-48b9-9fe8-7fd4ed36b89a",
   "metadata": {},
   "source": [
    "**Multiplicity test**\n",
    "\n",
    "\n",
    "In this example we will test the saturation effect of (5Kb) CRE tile in the context of an enhancing context sequence of ZNF93 gene. \n",
    "\n",
    "To run the test we need:\n",
    "- a loaded model\n",
    "- onehot encoded sequence (WT) of the sequence\n",
    "- list of start and end coordiantes of the TSS tile\n",
    "- list of start and end coordiantes to place the CRE in for the control case (we will use the original position)\n",
    "- onehot encoded sequence of the CRE tile\n",
    "- list of lists of start and end coordiantes to embed/test CRE effect\n",
    "- num_shuffle - number of shuffled\n",
    "- number of iterations to run or copies to insert\n",
    "- np.argmax or np.argmin to optimize for enhancers or silencers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "fa851a1f-f266-4142-821b-333d08aa684b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enhancing tile at position 105804 - 110804\n",
      "TSS tile at center position 95804 - 100804\n"
     ]
    }
   ],
   "source": [
    "\n",
    "perturb_window = 5000\n",
    "N_shuffles = 10\n",
    "num_copies = 3\n",
    "tss_tile, cre_tiles = utils.set_tile_range(model.seq_length, perturb_window)\n",
    "enhancing_cre = [105804, 110804]\n",
    "print(f'Enhancing tile at position {enhancing_cre[0]} - {enhancing_cre[1]}')\n",
    "print(f'TSS tile at center position {tss_tile[0]} - {tss_tile[1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0193ee8-b6d3-4e89-ad06-7448ae91ef62",
   "metadata": {},
   "source": [
    "Multiplicity effect of an enhancing tile of ZNF93 gene TSS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "52086ff5-48a3-45b2-99ad-818cd054e529",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████████████████████████████████████████| 3/3 [08:41<00:00, 173.86s/it]\n"
     ]
    }
   ],
   "source": [
    "model.bin_index = [447, 448]\n",
    "cre_seq = wt_seq[enhancing_cre[0]:enhancing_cre[1]].copy()\n",
    "result_summary = creme.multiplicity_test(model, wt_seq, tss_tile, enhancing_cre, cre_seq,\n",
    "                                         cre_tiles, N_shuffles, num_copies, np.argmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "4fe7f59a-ab6a-4e25-934c-0b08b25f2798",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['only_tss_pred', 'tss_and_cre_pred', 'best_tss_signal', 'selected_tile_order', 'all_mutants'])"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_summary.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "2f2ac3f0-90d1-4c21-8f1a-51683cdc1d51",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TSS signal with CRE copies at the optimal positions: [3.631044265255332, 19.56945643424988, 32.70301297903061]\n",
      "TSS signal in shuffled context: 0.12562379240989685\n",
      "TSS signal in shuffled context + CRE at original position: 3.631044387817383\n",
      "Tile order: [[105804, 110804], [100804, 105804], [120804, 125804]]\n"
     ]
    }
   ],
   "source": [
    "print(f\"TSS signal with CRE copies at the optimal positions: {result_summary['best_tss_signal']}\")\n",
    "print(f\"TSS signal in shuffled context: {result_summary['only_tss_pred']}\")\n",
    "print(f\"TSS signal in shuffled context + CRE at original position: {result_summary['tss_and_cre_pred']}\")\n",
    "print(f\"Tile order: {result_summary['selected_tile_order']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a3bd00-43bd-4d16-8d35-e60ab8df038d",
   "metadata": {},
   "source": [
    "Now, let's recreate the final sequence at each iteration and plot the entire track (instead of only looking at the TSS signal)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "ce80c893-e966-4543-bb37-521538567a9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.bin_index = None\n",
    "background_seqs_wth_TSS = shuffle.dinuc_shuffle(wt_seq, 10)\n",
    "background_seqs_wth_TSS[:, tss_tile[0]:tss_tile[1], :] = wt_seq[tss_tile[0]:tss_tile[1], :]\n",
    "just_tss_pred = model.predict(background_seqs_wth_TSS)[0,:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "b6f13879-7303-48b4-b4a6-55f02d5a2169",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[105804, 110804]\n",
      "[100804, 105804]\n",
      "[120804, 125804]\n"
     ]
    }
   ],
   "source": [
    "CREs_and_TSS = background_seqs_wth_TSS.copy()\n",
    "preds_per_iter = []\n",
    "for tile_position in result_summary['selected_tile_order']:\n",
    "    print(tile_position)\n",
    "    CREs_and_TSS[:, tile_position[0]:tile_position[1],:] = cre_seq.copy()\n",
    "    preds_per_iter.append(np.array(model.predict(CREs_and_TSS)[:,:,0]).mean(axis=0))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "9be2fd38-9d9c-47d0-be26-0ef0b7de1ddc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fd7404818b0>"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABkoAAADFCAYAAAAbg64/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZjElEQVR4nO3deXxU9b3/8fc5s2WdhARIQAKigOwYxQrSq1apiL1erTy0Un5Vqbd6LW5V6y11XyoutK5IfbRV660WV7zWWxe0KhYBBQkVAUFEASFAgGSyzXrO749ZkkkmIQnZeT0fOY+Z+Z7v+Z7vmeU758wn3+/XsG3bFgAAAAAAAAAAwGHI7OoKAAAAAAAAAAAAdBUCJQAAAAAAAAAA4LBFoAQAAAAAAAAAABy2CJQAAAAAAAAAAIDDFoESAAAAAAAAAABw2CJQAgAAAAAAAAAADlsESgAAAAAAAAAAwGGr2wVKbNuWz+eTbdtdXRUAAAAAAAAAANDLdbtASWVlpXJyclRZWdnVVQEAAAAAAAAAAL1ctwuUAAAAAAAAAAAAdBYCJQAAAAAAAAAA4LBFoAQAAAAAAAAAABy2CJQAAAAAAAAAAIDDFoESAAAAAAAAAABw2HJ2dQUAAAAAdBzbtvXm12+qMLNQxf2Lu7o6AAAAOExFIhGFQqGurga6MbfbLdPsmr4dBEoAAACAXmxH1Q59XPqxHIZDw3OHK8ud1dVVAgAAwGHEtm2VlpaqvLy8q6uCbs40TQ0dOlRut7vT902gBAAAAOjFfAGfJCliR/Tpnk918qCTu7hGAAAAOJzEgyT9+/dXRkaGDMPo6iqhG7IsSzt37tSuXbs0ePDgTn+fECgBAAAAerGKYEXi/urdq/XdI74r02CqQgAAAHS8SCSSCJLk5+d3dXXQzfXr1087d+5UOByWy+Xq1H1zhQQAAAD0Yr6gL+n+pgOburA2AAAAOJzE5yTJyMjo4pqgJ4gPuRWJRDp93wRKAAAAgF4sPvRWpjNTkvRJ6SddWR0AAAAchhhuCy3Rle8TAiUAAABALxbvUfLdQd+VIUNfVXylstqyLq4VAAAAAHQfBEoAAACAXizeo2Rw9mAN7zNckrRq96qurBIAAAAAdCsESgAAAIBeKmJFVBWqkiTleHI0sWCiJGntnrUKRoJdWTUAAAAA6DYIlAAAAAC9VGWwUrZsOQyHMpwZGpY7TH08feSP+LWubF1XVw8AAADoti655BIZhqF77703Kf3VV19tl7k0gsGg7r//fk2YMEEZGRnq27evpkyZoqeeekqhUCipDoZhyOVyaejQobrxxhvl9/uTyornabgsWrTokOt5uHB2dQUAAAAAdIzKYKUkyev2Ji7mJhZM1JJtS/RJ6Scq7l/MxJoAAABAE9LS0nTffffp8ssvV58+fdqt3GAwqGnTpmnt2rW66667NGXKFHm9Xq1YsULz589XcXGxjj32WEnSmWeemQierF69WhdffLEMw9B9992XVOZTTz2lM888MyktNze33erc2xEoAQAAAHqpimCFpGigJO7Y/sfqve3vqbSmVDuqdqgou6irqgcAAIDDkG3bqqmp6fT9ZmRktPqfhKZOnaovv/xS8+bN0/33399kvpdfflm33nqrvvzySw0YMEBXXXWVrr/++ibzP/TQQ1q6dKlWrVql4uLiRPpRRx2l888/X8Fg3TC5Ho9HhYWFkqSioiJNnTpVS5YsaRQoyc3NTeRD67Vq6K2FCxdq/Pjx8nq98nq9mjx5st54443Eer/frzlz5ig/P19ZWVmaMWOGdu/e3e6VBgAAAHBw8Yncczw5ibQMV4bG9h0rSVpVyqTuAAAA6Fw1NTXKysrq9KUtwRmHw6F77rlHjz76qHbs2JEyz+rVq3XBBRfowgsv1Geffabbb79dt9xyi55++ukmy3322Wc1derUpCBJnMvlUmZmZsrt1q1bp48++khut7vVx4LmtSpQMmjQIN17771avXq1Vq1apdNOO03nnHOOPv/8c0nSL37xC/3tb3/Tiy++qA8++EA7d+7Ueeed1yEVBwAAANA8XzAaKMl2ZyelTyyMTur++b7PVR2q7vR6AQAAAD3FD3/4Qx177LG67bbbUq7/3e9+p9NPP1233HKLRowYoUsuuURXXnmlHnjggSbL3Lx5s0aOHNmi/b/++uvKyspSWlqaxo0bpz179uiXv/xlo3wzZ85sFBzatm1byw4SrRt66+yzz056/Jvf/EYLFy7UihUrNGjQIP3pT3/Sc889p9NOO01SdFy0UaNGacWKFZo0aVL71RoAAADAQcUDJfV7lEjSEVlHaGDmQO2s3qmSPSWacsSUrqgeAAAADkMZGRmqqqrqkv221X333afTTjtNN9xwQ6N1GzZs0DnnnJOUNmXKFD300EOKRCJyOByNtrFtu8X7/t73vqeFCxequrpaDz74oJxOp2bMmNEo34MPPqipU6cmpQ0cOLDF+znctXmOkkgkohdffFHV1dWaPHmyVq9erVAolPRijBw5UoMHD9by5cubDJQEAgEFAoHEY5/P19YqAQAAAKgnHiipP0dJ3AmFJ+h/t/yvVu1epckDJ8s0WtXZHAAAAGgTwzCaHFqquzr55JM1bdo0zZ07V5dccskhlzdixAht3LixRXkzMzM1bNgwSdKTTz6pCRMm6E9/+pMuvfTSpHyFhYWJfGi9Vl8NffbZZ8rKypLH49F//dd/afHixRo9erRKS0vldruVm5ublL+goEClpaVNljdv3jzl5OQklqIiJpMEAAAA2kNFoPFk7nFj+o5RujNd5YFybT6wubOrBgAAAPQo9957r/72t79p+fLlSemjRo3SsmXLktKWLVumESNGpOxNIkk//vGP9c4772jNmjWN1oVCIVVXpx4e1zRN/frXv9bNN9+s2traNh4JUml1oOSYY45RSUmJVq5cqSuuuEIXX3yx1q9f3+YKzJ07VxUVFYll+/btbS4LAAAAQFTYCifmH/F6GgdKXKZLxf2jk0eu2s2k7gAAAEBzxo0bp1mzZumRRx5JSr/++uv17rvv6q677tKmTZv05z//WY899ljKYbrirr32Wk2ZMkWnn366FixYoLVr1+qrr77SCy+8oEmTJmnz5qb/ken888+Xw+HQggULktLLy8tVWlqatDQVcEFjrQ6UuN1uDRs2TMcff7zmzZunCRMm6OGHH1ZhYaGCwaDKy8uT8u/evVuFhYVNlufxeOT1epMWAAAAAIemKlglW7YchkMZztTjMR9fcLwkaUv5Fu337+/M6gEAAAA9zp133inLspLSjjvuOL3wwgtatGiRxo4dq1tvvVV33nlns0N0eTweLVmyRDfeeKOeeOIJTZo0SSeccIIeeeQRXX311Ro7dmyT2zqdTl155ZW6//77kwIhs2fP1oABA5KWRx999JCP+XBh2K2ZOSaF0047TYMHD9bDDz+sfv366a9//WtiMpkvvvhCI0eObHaOkoZ8Pp9ycnJUUVFB0AQAAABoo22+bXrq86fUx9NHVx93dZP5nt3wrL4s/1InDTxJ3x/y/U6sIQAAAHo7v9+vrVu3aujQoUpLS+vq6qCb68r3S6smc587d66mT5+uwYMHq7KyUs8995zef/99vfXWW8rJydGll16q6667Tnl5efJ6vbrqqqs0efLkFgdJAAAAALSP+ETuOZ6cZvNNLJioL8u/1Jo9a3Rq0alyma7OqB4AAAAAdButCpTs2bNHF110kXbt2qWcnByNHz9eb731lr7//eh/nj344IMyTVMzZsxQIBDQtGnT9Pjjj3dIxQEAAAA0rbmJ3Osb3me4cj25Kg+U6/Oyz3Vs/2M7oXYAAAAA0H20KlDypz/9qdn1aWlpWrBgQaOJZAAAAAB0rniPkoMFSkzD1PEFx+vdbe9q1e5VBEoAAAAAHHZaPZk7AAAAgO7PF4gFSjwHn/evuH+xHIZD31Z9q2+rvu3oqgEAAABAt0KgBAAAAOiFWtqjRJIyXZkakz9GkrSqdFWH1gsAAAAAuhsCJQAAAEAv1NLJ3OMmFk6UJK0rW6eaUE2H1QsAAAAAuhsCJQAAAEAvE7bCqgpVSWpZjxJJGpQ1SIUZhQrbYZXsKenA2gEAAABA90KgBAAAAOhlKoOVkiSn4VS6M71F2xiGoRMKT5Akrdq9SrZtd1j9AAAAAKA7IVACAAAA9DKJ+Uk8XhmG0eLtxvYdqzRHmg4EDmhL+ZaOqh4AAACAHswwDL366qtdXY12RaAEAAAA6GV8gZZP5F6f2+HWhP4TJEmf7P6k3esFAAAA9BRLly7V2WefrYEDB7YqMBAMBnX//fdrwoQJysjIUN++fTVlyhQ99dRTCoVCkqRLLrlEhmHIMAy5XC4NHTpUN954o/x+f1JZ8TwNl0WLFrX34bbKrl27NH369C6tQ3tzdnUFAAAAALSvRI+SVgZKJGliwUSt3LVSmw9sVrm/XLlpue1cOwAAAKD7q66u1oQJE/TTn/5U5513Xou2CQaDmjZtmtauXau77rpLU6ZMkdfr1YoVKzR//nwVFxfr2GOPlSSdeeaZieDJ6tWrdfHFF8swDN13331JZT711FM688wzk9Jyc3Pb4xDbrLCwsEv33xHoUQIAAAD0MvWH3mqtvul9dVTOUbJla/Xu1e1dNQAAAKBHmD59uu6++2798Ic/bPE2Dz30kJYuXap3331Xc+bM0bHHHqujjjpKP/7xj7Vy5UoNHz48kdfj8aiwsFBFRUU699xzNXXqVC1ZsqRRmbm5uSosLExa0tLSmqxDeXm5Lr/8chUUFCgtLU1jx47V66+/nlj/8ssva8yYMfJ4PDryyCP129/+Nmn7I488UnfddZdmzpypzMxMHXHEEVqwYEFSnoY9bLZv364LLrhAubm5ysvL0znnnKOvv/46sf7999/Xd77zHWVmZio3N1dTpkzRN99809KntVMQKAEAAAB6mYpAhaS29SiRlJjU/dM9nypkhdqtXgAAAIBt2wqEI52+2Lbd4cf27LPPaurUqSouLm60zuVyKTMzM+V269at00cffSS3231I+7csS9OnT9eyZcv0l7/8RevXr9e9994rh8MhSVq9erUuuOACXXjhhfrss890++2365ZbbtHTTz+dVM4DDzygCRMmaM2aNfrVr36la665JmUQR5JCoZCmTZum7Oxsffjhh1q2bJmysrJ05plnKhgMKhwO69xzz9Upp5yif/3rX1q+fLkuu+yyVs2l2BkYegsAAADoZSqDlZLaHigZ0WeEvG6vfEGfNuzboPH9xrdn9QAAAHAYC0Ys3f7a+k7f7+3/MVoep6ND97F582adeuqpLcr7+uuvKysrS+FwWIFAQKZp6rHHHmuUb+bMmYlAR9z69es1ePDgRnnfeecdffzxx9qwYYNGjBghSTrqqKMS63/3u9/p9NNP1y233CJJGjFihNavX68HHnhAl1xySSLflClT9Ktf/SqRZ9myZXrwwQf1/e9/v9E+n3/+eVmWpT/+8Y+J4MdTTz2l3Nxcvf/++5o4caIqKir07//+7zr66KMlSaNGjWrRc9SZ6FECAAAA9DLxHiU5npw2bW8apo4vOF6S9Ekpk7oDAAAALdGaXivf+973VFJSopUrV+riiy/W7NmzNWPGjEb5HnzwQZWUlCQtAwcOTFlmSUmJBg0alAiSNLRhwwZNmTIlKW3KlCnavHmzIpFIIm3y5MlJeSZPnqwNGzakLHPt2rX68ssvlZ2draysLGVlZSkvL09+v19btmxRXl6eLrnkEk2bNk1nn322Hn74Ye3atavZ56Yr0KMEAAAA6EXCVljV4WpJbe9RIknH9T9OS3cs1Y6qHSqtLlVhZu+bsBEAAACdz+0wdft/jO6S/Xa0ESNGaOPGjS3Km5mZqWHDhkmSnnzySU2YMEF/+tOfdOmllyblKywsTOQ7mPT09NZVuB1UVVXp+OOP17PPPttoXb9+/SRFe5hcffXVevPNN/X888/r5ptv1pIlSzRp0qTOrm6T6FECAAAA9CLxYbechlPpzrZfKGW5szQqL9olnl4lAAAAaC+GYcjjdHT60hlzYvz4xz/WO++8ozVr1jRaFwqFVF1dnXI70zT161//WjfffLNqa2vbvP/x48drx44d2rRpU8r1o0aN0rJly5LSli1bphEjRiQN77VixYqkPCtWrGhyuKzjjjtOmzdvVv/+/TVs2LCkJSenrod7cXGx5s6dq48++khjx47Vc88919bD7BAESgAAAIBexBf0SZK8Hu8hXwzGJ3X/195/qTbc9gs2AAAAoKepqqpKDHUlSVu3blVJSYm2bdvW5DbXXnutpkyZotNPP10LFizQ2rVr9dVXX+mFF17QpEmTtHnz5ia3Pf/88+VwOLRgwYKk9PLycpWWliYtTQVcTjnlFJ188smaMWOGlixZoq1bt+qNN97Qm2++KUm6/vrr9e677+quu+7Spk2b9Oc//1mPPfaYbrjhhqRyli1bpvvvv1+bNm3SggUL9OKLL+qaa65Juc9Zs2apb9++Ouecc/Thhx9q69atev/993X11Vdrx44d2rp1q+bOnavly5frm2++0dtvv63Nmzd3u3lKCJQAAAAAvUh8fpJDGXYrrii7SP0z+itsh/XF/i8OuTwAAACgp1i1apWKi4tVXFwsSbruuutUXFysW2+9tcltPB6PlixZohtvvFFPPPGEJk2apBNOOEGPPPKIrr76ao0dO7bJbZ1Op6688krdf//9SYGQ2bNna8CAAUnLo48+2mQ5L7/8sk444QTNnDlTo0eP1o033piYf+S4447TCy+8oEWLFmns2LG69dZbdeeddyZN5C5FAyrx47/77rv1u9/9TtOmTUu5v4yMDC1dulSDBw/Weeedp1GjRunSSy+V3++X1+tVRkaGNm7cqBkzZmjEiBG67LLLNGfOHF1++eVNHkNXMOzWzDDTCXw+n3JyclRRUSGv99Av7gAAAIDDyT+//afe3fauxvcdrx8O/+Ehl/fG1jf0cenHOmngSfr+kO+3Qw0BAABwuPD7/dq6dauGDh2qtLS0rq4OWuDII4/Utddeq2uvvbbT992V7xd6lAAAAAC9iC9QN/RWe+ib3leSVFZb1i7lAQAAAEB3Q6AEAAAA6EXic5TkuHMOkrNl8tPyJUn7ave1S3kAAAAA0N04u7oCAAAAANpPPFCS7c6OJoSD0sdPSDmDpDGtH4orPz0aKDngP6CIFZHDdLRbXQEAAAB0L19//XVXV6FL0KMEAAAA6EXiQ2/leGI9SvZ/Je37Uvrqfal8W6vL87q9cpkuWbJ0IHCgHWsKAAAAAN0DgRIAAACglwhZIVWHqyVFAxySpKrddRk2vdXqMg3DYPgtAAAAAL0agRIAAACgl6gMVkqSnIZT6c70aGL9QMnudVLFjlaXGx9+i0AJAAAAgN6IQAkAAADQS9QfdsswjGhiPFDiyozetqFXSd/0vpKkMn/ZIdcRAAAAALobAiUAAABALxGfyD0x7JYkVZZGb8f8UJIhlf5L8u1sVbkMvQUAAACgNyNQAgAAAPQSFYEKSZLXEwuUBKulYFX0/oDx0oAJ0fub325VuQy9BQAAAKA3I1ACAAAA9BLxOUoaTeSelis5PdLwM6KPd5bU9TRpgXigpDpcrdpwbTvVFgAAAEBP8/XXX8swDJWUlHR1VdoVgRIAAACgl6gIxnqUJAIle6K32YXR25wjpMLxkmxp85IWl+txeJTtzpYkldUyTwkAAAB6v3nz5umEE05Qdna2+vfvr3PPPVdffPHFQbfz+Xy66aabNHLkSKWlpamwsFBTp07VK6+8Itu2JUmnnnqqDMOQYRhKS0vTiBEjNG/evMR6qS4gkWpZsWJFhx33wRQVFWnXrl0aO3Zsl9WhIxAoAQAAAHqJ+GTuiaG34j1KsvrXZRoxLXr77eq6QEoLME8JAAAADicffPCB5syZoxUrVmjJkiUKhUI644wzVF1d3eQ25eXlOumkk/TMM89o7ty5+vTTT7V06VL96Ec/0o033qiKiopE3p/97GfatWuXvvjiC82dO1e33nqrfv/73zcq85133tGuXbuSluOPP75DjrklHA6HCgsL5XQ6u6wOHYFACQAAANBLNJrMvTIeKCmsy5QzSCoYq9b2KknMU+InUAIAAIDe780339Qll1yiMWPGaMKECXr66ae1bds2rV69usltfv3rX+vrr7/WypUrdfHFF2v06NEaMWKEfvazn6mkpERZWVmJvBkZGSosLNSQIUM0e/ZsjR8/XkuWND4/z8/PV2FhYdLicrmarMOOHTs0c+ZM5eXlKTMzUxMnTtTKlSsT6xcuXKijjz5abrdbxxxzjP7nf/4naXvDMLRw4UJNnz5d6enpOuqoo/TSSy8l1qcaemvdunWaPn26srKyVFBQoJ/85CcqK6vrif7SSy9p3LhxSk9PV35+vqZOndpswKkrECgBAAAAeoGQFVJNuEZS/aG3YvOQZBUkZ47PVfLtKqm6ZUNp9U3vK4keJQAAADhEti2FA52/1BvWqi3ivUHy8vJSrrcsS4sWLdKsWbM0cODARuuzsrJS9sKwbVsffvihNm7cKLfbfUh1rKqq0imnnKJvv/1Wr732mtauXasbb7xRlmVJkhYvXqxrrrlG119/vdatW6fLL79cs2fP1nvvvZdUzi233KIZM2Zo7dq1mjVrli688EJt2LAh5T7Ly8t12mmnqbi4WKtWrdKbb76p3bt364ILLpAk7dq1SzNnztRPf/pTbdiwQe+//77OO++8pGHGuoNW9Y+ZN2+eXnnlFW3cuFHp6ek66aSTdN999+mYY45J5PH7/br++uu1aNEiBQIBTZs2TY8//rgKCgqaKRkAAADAoYgPu+UyXUp3pkuRkFSzP7qy/tBbktRniNR/tLRnfbRXybEzD1o+Q28BAACgXUSC0hs3dv5+p98vOT1t2tSyLF177bWaMmVKk3NzlJWV6cCBAxo5cmSLynz88cf1xz/+UcFgUKFQSGlpabr66qsb5TvppJNkmsn9HaqqqlKW+dxzz2nv3r365JNPEgGdYcOGJdbPnz9fl1xyiX7+859Lkq677jqtWLFC8+fP1/e+971EvvPPP1//+Z//KUm66667tGTJEj366KN6/PHHG+3zscceU3Fxse65555E2pNPPqmioiJt2rRJVVVVCofDOu+88zRkyBBJ0rhx41r0HHWmVvUoacm4bL/4xS/0t7/9TS+++KI++OAD7dy5U+edd167VxwAAABAnfrDbhmGEZt/xJZcGZInu/EG8V4lOz6uC6g0I9GjxL9Plm21V7UBAACAbm/OnDlat26dFi1a1GSe1vaQmDVrlkpKSrRs2TJNnz5dN910k0466aRG+Z5//nmVlJQkLU0pKSlRcXFxk71eNmzYoClTpiSlTZkypVFvkcmTJzd63FSPkrVr1+q9995TVlZWYokHi7Zs2aIJEybo9NNP17hx43T++efrD3/4gw4cONDkMXSVVvUoefPNN5MeP/300+rfv79Wr16tk08+WRUVFfrTn/6k5557Tqeddpok6amnntKoUaO0YsUKTZo0qf1qDgAAACCh0fwkiYncCyTDaLxB3lCp7zFS2RfSl+9I4y9otvwcT44chkMRO6KKQIX6pPVpz+oDAADgcOFwR3t3dMV+2+DKK6/U66+/rqVLl2rQoEFN5uvXr59yc3O1cePGFpWbk5OT6O3xwgsvaNiwYZo0aZKmTp2alK+oqCipV0hz0tPTW5SvPVVVVenss8/Wfffd12jdgAED5HA4tGTJEn300Ud6++239eijj+qmm27SypUrNXTo0E6vb1MOaY6ShuOyrV69WqFQKOnFHDlypAYPHqzly5enLCMQCMjn8yUtAAAAAFonPvRWjicnmlA/UNKUEdOit9tWSLXN/1eXaZiJ4bfKals2rwkAAADQiGFEh8Dq7CXVPw81w7ZtXXnllVq8eLH+8Y9/HPRHfdM0deGFF+rZZ5/Vzp07G62PD0GVSlZWlq655hrdcMMNhzR3x/jx41VSUqL9+1P3GB81apSWLVuWlLZs2TKNHj06KW3FihWNHo8aNSplmccdd5w+//xzHXnkkRo2bFjSkpmZKSk6QfyUKVN0xx13aM2aNXK73Vq8eHFbD7NDtDlQkmpcttLSUrndbuXm5iblLSgoUGlpacpy5s2bp5ycnMRSVFTU1ioBAAAAh614j5Jsd2yYrUSgpH8TW0jKP1rKHybZEenLdw+6j/z02DwlfuYpAQAAQO82Z84c/eUvf9Fzzz2n7OxslZaWqrS0VLW1tU1u85vf/EZFRUU68cQT9cwzz2j9+vXavHmznnzySRUXFzc5t4gkXX755dq0aZNefvnlpPR9+/Yl9h1f/H5/yjJmzpypwsJCnXvuuVq2bJm++uorvfzyy4lODL/85S/19NNPa+HChdq8ebN+97vf6ZVXXtENN9yQVM6LL76oJ598Ups2bdJtt92mjz/+WFdeeWWTz9P+/fs1c+ZMffLJJ9qyZYveeustzZ49W5FIRCtXrtQ999yjVatWadu2bXrllVe0d+/eJgMvXaXNgZKWjMvWEnPnzlVFRUVi2b59+yGVBwAAAByO4oGSuh4le6K32YXNbzjizOjttuVSbXmzWeOBEnqUAAAAoLdbuHChKioqdOqpp2rAgAGJ5fnnn29ym7y8PK1YsUL/7//9P919990qLi7Wv/3bv+mvf/2rHnjgAeXk5DS77UUXXaTbb79dllU3J+DUqVOT9j9gwAC9+uqrKctwu916++231b9/f5111lkaN26c7r33XjkcDknSueeeq4cffljz58/XmDFj9MQTT+ipp57SqaeemlTOHXfcoUWLFmn8+PF65pln9Ne//rVRr5O4gQMHatmyZYpEIjrjjDM0btw4XXvttcrNzZVpmvJ6vVq6dKnOOussjRgxQjfffLN++9vfavr06U0+F12hVXOUxDU1LlthYaGCwaDKy8uTepXs3r1bhYWpL9A8Ho88Hk9bqgEAAAAgJj70ltftlSyrXo+SgwRK8odJeUdJ+7+StvxDGnte01ljQ2/tq6VHCQAAAHq3tg6BlZOTo3nz5mnevHlN5nn//fdTpv/+979P3D/yyCPbVIchQ4bopZdeanL9FVdcoSuuuKLZMgYOHKi333475bpU9Ro+fLheeeWVlPlHjRrVaO7z7qhVPUoONi7b8ccfL5fLpXffreu2/8UXX2jbtm2aPHly+9QYAAAAQCMVwej8gV63V6rdL1lhyXRK6QeZdN0w6nqVfPOR5G96zkCG3gIAAADQG7WqR8mcOXP03HPP6X//938T47JJ0ShZenq6cnJydOmll+q6665TXl6evF6vrrrqKk2ePFmTJk3qkAMAAAAADnchK6TacHSsZK/HK+37Kroiq0AyW/C/UX1HSH2OlA58He1VMubc1NnS+0qSKoOVCkQC8jjoGQ4AAACg52tVoGThwoWS1GjMsqeeekqXXHKJJOnBBx+UaZqaMWOGAoGApk2bpscff7xdKgsAAACgsfiwW27TrTRHmlQZ/YemZidyr88wpOHTpI+fkL5ZJg07XfJkN8qW7kxXpjNT1eFq7avdp4FZA9vrEAAAAAB0A20dcqyna1WgpCVPUlpamhYsWKAFCxa0uVIAAAAAWi4x7JbHK8Mw6iZyP9j8JPX1HyXlFEkV26Ut70mj/yNltvz0fFVXEigBAAAA0Hu0ao4SAAAAAN1PZbBSUmx+EkmqamWPEil5rpKv/ykFqlJmY54SAAAAAL0NgRIAAACgh6sI1JvI3bbr9SgpaF1BBWMk7yApEpC2fpAyS9+06DwlZbVlba4vAAAAAHQnBEoAAACAHs4XjM5R4nV7pUClFKqRZLSuR4kU61VyRvT+1qVSsKZRlkSPklp6lAAAAADoHQiUAAAAAD1cfDJ3r8db15skI09yuFpfWOF4KXugFPan7FVSf+itw3WiRwAAAAC9C4ESAAAAoIdL6lFStTua2JqJ3Our36vkqw+kUG3S6j6ePjJlKmSFEvsFAAAAgJ6MQAkAAADQw8UDFjmenLqJ3LNbOT9JfQOOjQZawrXS1g+TVjlMh/qk9ZHE8FsAAADA4ebrr7+WYRgqKSnp6qq0KwIlAAAAQA8WioRUG472+sh2Z0uV8R4lhxAoSepV8r4U8ietrj/8FgAAANAbLVy4UOPHj5fX65XX69XkyZP1xhtvHHQ7n8+nm266SSNHjlRaWpoKCws1depUvfLKK4mha0899VQZhiHDMJSWlqYRI0Zo3rx5SUPbxgMSqZYVK1Z02HEfTFFRkXbt2qWxY8d2WR06grOrKwAAAACg7SqCFZIkt+lWmiOt3tBbrZzIvaEBxVLmm1L1Humbf0rDpiZW9U3vq00HNqmstuzQ9gEAAAB0U4MGDdK9996r4cOHy7Zt/fnPf9Y555yjNWvWaMyYMSm3KS8v13e/+11VVFTo7rvv1gknnCCn06kPPvhAN954o0477TTl5uZKkn72s5/pzjvvVCAQ0D/+8Q9ddtllys3N1RVXXJFU5jvvvNNof/n5+R1yzC3hcDhUWNjGYX67MXqUAAAAAD1Y/WG3jHBA8pdHV7R1jpI405SGx3qVbHlPCgcSq/LTYj1KGHoLAAAAvdTZZ5+ts846S8OHD9eIESP0m9/8RllZWc325vj1r3+tr7/+WitXrtTFF1+s0aNHa8SIEfrZz36mkpISZWVlJfJmZGSosLBQQ4YM0ezZszV+/HgtWbKkUZn5+fkqLCxMWlwuV5N12LFjh2bOnKm8vDxlZmZq4sSJWrlyZWL9woULdfTRR8vtduuYY47R//zP/yRtbxiGFi5cqOnTpys9PV1HHXWUXnrppcT6VENvrVu3TtOnT1dWVpYKCgr0k5/8RGVldf9U9dJLL2ncuHFKT09Xfn6+pk6dqurq6iaPoSsQKAEAAAB6MF8gGijJdmdHe39IkidbcmcceuFHHC9l9JWCVdLOkkQyQ28BAACgrWzbVjAS7PSl/rBWrRWJRLRo0SJVV1dr8uTJKfNYlqVFixZp1qxZGjhwYKP1WVlZcjobD/Bk27Y+/PBDbdy4UW63u811lKSqqiqdcsop+vbbb/Xaa69p7dq1uvHGG2VZliRp8eLFuuaaa3T99ddr3bp1uvzyyzV79my99957SeXccsstmjFjhtauXatZs2bpwgsv1IYNG1Lus7y8XKeddpqKi4u1atUqvfnmm9q9e7cuuOACSdKuXbs0c+ZM/fSnP9WGDRv0/vvv67zzzjuk16MjMPQWAAAA0IMlTeTeHvOT1Gea0sBi6csl0v4t0uATJUWH3pKkikCFQlZILrPp/2gDAAAA6gtZIc37eF6n73fud+bK7WhdIOKzzz7T5MmT5ff7lZWVpcWLF2v06NEp85aVlenAgQMaOXJki8p+/PHH9cc//lHBYFChUEhpaWm6+uqrG+U76aSTZJrJ/R2qqqpSlvncc89p7969+uSTT5SXlydJGjZsWGL9/Pnzdckll+jnP/+5JOm6667TihUrNH/+fH3ve99L5Dv//PP1n//5n5Kku+66S0uWLNGjjz6qxx9/vNE+H3vsMRUXF+uee+5JpD355JMqKirSpk2bVFVVpXA4rPPOO09DhgyRJI0bN65Fz1FnokcJAAAA0IPFAyVet7fe/CTtFCiRpPyjo7f7v0okZTgzlOZIky2b4bcAAADQax1zzDEqKSnRypUrdcUVV+jiiy/W+vXrU+ZtbQ+JWbNmqaSkRMuWLdP06dN100036aSTTmqU7/nnn1dJSUnS0pSSkhIVFxcngiQNbdiwQVOmTElKmzJlSqPeIg17zUyePLnJHiVr167Ve++9p6ysrMQSDxZt2bJFEyZM0Omnn65x48bp/PPP1x/+8AcdOHCgyWPoKvQoAQAAAHqwikB0MvdooGRLNLE9AyV9jpRkSNV7Jb9PSvPKMAzlp+fr26pvtc+/T4WZvW8yRwAAAHQMl+nS3O/M7ZL9tpbb7U70yDj++OP1ySef6OGHH9YTTzzRKG+/fv2Um5urjRs3tqjsnJycRNkvvPCChg0bpkmTJmnq1KlJ+YqKipJ6hTQnPT29RfnaU1VVlc4++2zdd999jdYNGDBADodDS5Ys0UcffaS3335bjz76qG666SatXLlSQ4cO7fT6NoUeJQAAAEAPVhmslNSgR0l2OwYuXOmSNzbGcr1eJUzoDgAAgLYwDENuh7vTF8MwDrnulmUpEAikXGeapi688EI9++yz2rlzZ6P18SGoUsnKytI111yjG2644ZDm7hg/frxKSkq0f//+lOtHjRqlZcuWJaUtW7as0XBiDSesX7FihUaNGpWyzOOOO06ff/65jjzySA0bNixpyczMlBR9zadMmaI77rhDa9askdvt1uLFi9t6mB2CQAkAAADQgyV6lLgyo70+JCmrf/vuJO+o6G29QEl8nhICJQAAAOiN5s6dq6VLl+rrr7/WZ599prlz5+r999/XrFmzmtzmN7/5jYqKinTiiSfqmWee0fr167V582Y9+eSTKi4ubnJuEUm6/PLLtWnTJr388stJ6fv27VNpaWnS4vf7U5Yxc+ZMFRYW6txzz9WyZcv01Vdf6eWXX9by5cslSb/85S/19NNPa+HChdq8ebN+97vf6ZVXXtENN9yQVM6LL76oJ598Ups2bdJtt92mjz/+WFdeeWXKfc6ZM0f79+/XzJkz9cknn2jLli166623NHv2bEUiEa1cuVL33HOPVq1apW3btumVV17R3r17mwy8dBWG3gIAAAB6qGAkKH8kepHkDYck25IcHiktt313lDdU+vrD5B4l6bEeJX4CJQAAAOh99uzZo4suuki7du1STk6Oxo8fr7feekvf//73m9wmLy9PK1as0L333qu7775b33zzjfr06aNx48bpgQceUE5OTrPbXnTRRbr99tt13nnnJdIbDsUlSX/961914YUXNkp3u916++23df311+uss85SOBzW6NGjtWDBAknSueeeq4cffljz58/XNddco6FDh+qpp57SqaeemlTOHXfcoUWLFunnP/+5BgwYoL/+9a9NTmI/cOBALVu2TP/93/+tM844Q4FAQEOGDNGZZ54p0zTl9Xq1dOlSPfTQQ/L5fBoyZIh++9vfavr06U0+F13BsA+lL08H8Pl8ysnJUUVFhbxeb1dXBwAAAOi2ymrLtKBkgTwOj35VNF1a9aSUUySdfMPBN26N2gPSO7dLMqQz75VcadpTs0cL1y6Ux+HRf5/w3+0ylAEAAAB6F7/fr61bt2ro0KFKS0vr6uqgBQzD0OLFi3Xuued2+r678v3C0FsAAABAD+UL+CTF5yfZE01sz4nc49L7SOl5kmyp/BtJUp+0PjJkKBAJqDpU3f77BAAAAIBOQqAEAAAA6KF8wXqBksrSaGJ7TuReX4N5SlymS7meXEnRni0AAAAA0FMxRwkAAADQQyUCJR6vVLYpmtjeE7nH5R0lfbuq0TwlBwIHtM+/T0fmHNkx+wUAAADQabrZTB2dhh4lAAAAQA9VEaiQJHld2fWG3urgHiUHvpasiKS6Cd3pUQIAAACgJyNQAgAAAPRQ8R4lOYZDigQkw5Qy+3bMzrILJVeGFAlKFTskSflp0UDJvtp9HbNPAAAA9AqHay8FtE5Xvk8IlAAAAAA9VDxQkh0KRhMy+0mmo2N2ZhiN5imJ9yjZ5ydQAgAAgMZcLpckqaampotrgp4gGIxe1zgcHXRN0wzmKAEAAAB6KF8g1qMkGLvwzCro2B3mDZV2r4sGSo7+nvqmR3uvlPvLFbbCcppcXgAAAKCOw+FQbm6u9uyJDhObkZEhwzC6uFbojizL0t69e5WRkSGns/OvK7iSAQAAAHqgYCQof8QvSfL6q6KJHR4oqdejxLaV7cqW23QraAV1wH9A/TL6dez+AQAA0OMUFkbn0IsHS4CmmKapwYMHd0kwjUAJAAAA0APFh93yODzy1MSGvmoiULLu2wr1yXTriNz0Q9tpzmDJdErBKql6r4ys/spPz9eu6l3a599HoAQAAACNGIahAQMGqH///gqFQl1dHXRjbrdbptk1s4UQKAEAAAB6oPiwW163VyrfGU3Mbhwo2bavRs+u3CaP09Qvvj9COemutu/U4ZRyh0j7t0R7lWT1V35aLFDChO4AAABohsPh6JK5J4CWYDJ3AAAAoAeK9yjJcaRLgcpoYooeJet2VkiSAmFLf1u789B33GBC9/g8JWW1ZYdeNgAAAAB0AQIlAAAAQA9UEYgGQLyWFU1Iy5WcnqQ8tm3r81igRJI+3+nT+p2+Q9txg0BJfnq+JGmfnx4lAAAAAHomAiUAAABADxTvUeKNhKMJKXqT7PYFtL86JJfD0OSjowGN19bulD8UafuO+xwpyZCq90p+Hz1KAAAAAPR4BEoAAACAHigRKAkFogkp5ieJ9yYZ3j9L08cWKj/TrYrakJas3932HbszJO/A6P39XykvLU+SVBuuVU2opu3lAgAAAEAXaXWgZOnSpTr77LM1cOBAGYahV199NWm9bdu69dZbNWDAAKWnp2vq1KnavHlze9UXAAAAgOpN5h6IBSeyChvliQ+zNXqgVy6HqXOLowGO5V/t0/b9hxDUqDf8ltvhjk4oLzGhOwAAAIAeqdWBkurqak2YMEELFixIuf7+++/XI488ot///vdauXKlMjMzNW3aNPn9/kOuLAAAAICoRI8Sf2wOkgZDb+2vDmpnhV+GIY0sjAYyhvXPVnFRrmxbWrzmW0Usu207zxsa20lsnpI05ikBAAAA0HM5W7vB9OnTNX369JTrbNvWQw89pJtvvlnnnHOOJOmZZ55RQUGBXn31VV144YWNtgkEAgoEAonHPt8hTi4JAAAA9HKBSED+iF+yLHn9VZJhSln9k/Js2BU9rx6an6lMT91p/1njB+iL3ZXaVeHXsi/LdPKIfq2vQLxHie9bKRxQ3/S+2urbyjwlAAAAAHqkdp2jZOvWrSotLdXUqVMTaTk5OTrxxBO1fPnylNvMmzdPOTk5iaWoqKg9qwQAAAD0OvFht9IsSx7DlFwZkic7KU98fpIxA71J6Vkep84aFx2m650Nu7W/Otj6CqT3kdLzJNuSDnyj/PRYjxKG3gIAAADQA7VroKS0tFSSVFCQ3O2/oKAgsa6huXPnqqKiIrFs3769PasEAAAA9DqJYbfiCVn9JcNIrK8KhPX1vugcJKMGeNXQcYP76Oh+mQpFbP1vybey7TYMwZWYp2QLQ28BAAAA6NHaNVDSFh6PR16vN2kBAAAA0LREoCRiRRMaTOS+cZdPti0dkZumPpnuRtsbhqFzjj1CTtPQpt1V+teOitZXot6E7n3T+0bv+vfLsq3WlwUAAAAAXahdAyWFhdELtN27dyel7969O7EOAAAAwKGJD73ljYSiCQ3mJ1kfm59kzMCcJsvol+3RaSOj273+r52qCYZbV4l4oOTAN8pxZclpOBWxIzrgP9C6cgAAAACgi7VroGTo0KEqLCzUu+++m0jz+XxauXKlJk+e3J67AgAAAA5biR4lQX80Ibvun5L8oYg2766SJI0e2Hxv7X8b3lf9sz2qCkT0xmeph8ptUnZhdG6USEBG5U7lpedJYvgtAAAAAD1PqwMlVVVVKikpUUlJiaToBO4lJSXatm2bDMPQtddeq7vvvluvvfaaPvvsM1100UUaOHCgzj333HauOgAAAHB48gV9km0rJ1AdTciqmyNw8+4qhS1bfbPc6p/tabYcp8PUD4uPkCSt+uaAtpZVt7wShiH1GRq9v/+runlKmNAdAAAAQA/T6kDJqlWrVFxcrOLiYknSddddp+LiYt16662SpBtvvFFXXXWVLrvsMp1wwgmqqqrSm2++qbS0tPatOQAAAHCY8gV8UjgorwzJdEqx3hyStH5XdL6R0QO8MupN8N6UI/tm6sSh0e0Xf7pD4Ugr5hjJbzxPCT1KAAAAAPQ0ztZucOqpp8q27SbXG4ahO++8U3feeechVQwAAABAahXBCilUI6/hjPYmMaP//xSOWNpYWimp+flJGpo2plDrd/m0tyqo97/Yq6mjCw6+kVQ3T8m+LcoffFz0Lj1KAAAAAPQw7TpHCQAAAICOFYgEFIgEYoESV9JE7l+VVcsfspSd5lRRXnqLy0x3O/Tv4wdIkj7YtFd7Kv0t2zBncLRHS7BKfW2HJAIlAAAAAHoeAiUAAABAD+ILRCdyT4uE5TbMpPlJ1u+MrmvpsFv1jTsiRyMLsxW2bL265ttme5EnOJxS7mBJUn5NuSSpMlQZDeQAAAAAQA9BoAQAAADoQSoC0TlIvJFwNCEWKLFtWxt2xQIlA72tLtcwDP3HhIFyOwxtLavRqm8OtGzD2PBbaRU7lOXKkkSvEgAAAAA9C4ESAAAAoAepDFVKtuQNB6MJsUDJ9v218vnD8jhNHdU3s01l98l06/ujCyVJb3xWqkp/6OAb5dVN6J6fli9JKqsta9P+AQAAAKArECgBAAAAepCKQIUUCclr2ZKMxBwln++M9jQZWZgtp6Ptp/knHZ2vI3LTVBuK6P/+tevgG/QZGq1H9R7lu6IBGnqUAAAAAOhJCJQAAAAAPYgv6JPCtcoxnFJGnuRwybZtrY8NuzVmYM4hlW+ahn543CAZhrR2R4U27a5sfgN3hpQdnQg+PxSdm2Sfn0AJAAAAgJ6DQAkAAADQg/gCPilUK6/hkrKiw2TtqQyorCoop2loeEHWIe/jiNx0TTm6ryTp1TXfKhCONL9BbPit/NpoUIWhtwAAAAD0JARKAAAAgB7EF/RJoRp5TWdi2K31O6O9SYb1z1Kay9Eu+5k6ur9yM1w6UBPSPzbsaT5zLFDStzrak2Rf7T7Ztt0u9QAAAACAjkagBAAAAOhBooGSWI+S7GiPkrpht7ztth+P06Fzjh0oSfrnl2XaWV7bdOa8oZKkPpVlcthS2A5H6wkAAAAAPQCBEgAAAKCH8If9CkQC0R4lRrRHSXlNUDsO1MowpJED2i9QIkkjC70ad0SOLFtavOZbWVYTvUQy8qS0XJmy1cc2JDGhOwAAAICeg0AJAAAA0EP4gj7JiigtEpHbMKWsgsSwW0fmZyjL42z3ff77hAFKc5nacaBWK75qJviRf3T0JhyWxDwlAAAAAHoOAiUAAABAD5GYyN10Su4syZ2ZGHZr9ICcDtmnN82lM8dEh/h6e/1uVdSEUmeMz1MSqJEk7fPTowQAAABAz0CgBAAAAOgh6uYncUrZhaoJhrW1rFqSNLod5ydp6DtD8zQkP0OBsKXX1n6beqL2WKAkv7ZCsm16lAAAAADoMQiUAAAAAD1EPFCSY7ikrAJt2OWTZUsDctKUl+nusP0ahqEfFh8hhymt31Wpz3emmKg9e4DkTFdfW1KwmjlKAAAAAPQYBEoAAACAHiIaKKmbyD0+P8mYDuxNElfgTdPJw/tJkv62dqf8oUhyBsOQ8oYq33RLgUpVBCsUijQxTBcAAAAAdCMESgAAAIAeoiJQIYVr5TVcCqb30+Y9VZI6dtit+r43sr/6Zrnl84f11ueljTPkHaUMw6n0UK0k5ikBAAAA0DMQKAEAAAB6CF+gXAr55TWd+rI2S6GIrbxMlwq9aZ2yf5fD1LnFR0iSVm7dr237apIzxOcpCfolWwy/BQAAAKBHIFACAAAA9AC2bctXvUeSLa8zU+vKDEnS6AE5Mgyj0+pxdL8sHT+kj2xbemXNDkWsehO75w6WTKf6WrYU9tOjBAAAAECPQKAEAAAA6AH8Eb+CgQpJUnbWQG3cHR12qzPmJ2norHGFynQ7tNsX0NLNe+tWOFxSTpHyDbcU8NGjBAAAAECPQKAEAAAA6AEqg5VSqFbpcqjcyFNtKKIsj0OD8zI6vS4Zbqd+MH6AJOm9jXtUVhWoW5l3VHRCd79PZbVlnV43AAAAAGgtAiUAAABAD1ARqJBCtfKaTn3lz5IkjRrglWl23rBb9R1blKth/aPzpLy65lvZdmwIrryj1Nd0S4FK7fPvq0sHAAAAgG6KQAkAAADQA/iCvmigRC5tqEyXJI0ZmNNl9TEMQ+ceO1Auh6Ete6u1Znt5dEXeUPUxXDLCtQoEq1QVquqyOgIAAABASxAoAQAAAHoAX6xHiccyVBrJkcdp6uh+mV1ap/wsj04b2V+S9Pd/7VJ1ICy5M+X0HqFcwxXtVcI8JQAAAAC6OQIlAAAAQA/gqy6V7IgUMFXl7KNjCrPldHT96fy/De+nQm+aqoMR/f2zXdFE5ikBAAAA0IN0/ZUVAAAAgIPyVUaDEP5wtmzDodEDvF1coyiHaeiHxUfIMKRPt5Xryz1VsXlKPFLAp31+epQAAAAA6N4IlAAAAAA9QEVVqSKWrRqrjxymdExhdldXKWFwfoZOHJonSfrfkm8VyjlS+YZLClarrLq0i2sHAAAAAM0jUAIAAAB0c7Ztq9K/T/5QRBFnoY7ul6U0l6Orq5Vk2phCedOdKqsK6r1vQuqb3leSrX3lX3d11QAAAACgWQRKAAAAgG7OH/ErGKhUMGwp5CjUmIE5XV2lRtJcDp09fqAk6YPNZXJnHi1JKq/cobAV7sqqAQAAAECznF1dAQAAAADN8wV8igRrZEYM1Wb116gB3WfYrfrGDPRq9IBsrd9VqU/25cstU0G/Twf8B9Qvo19XV6/DWbal/f792lOzR2W1ZXIYDuV6cpXjyVGuJ1eZrkwZhtHV1QQAAADQAIESAAAAoJvzVe9RMBBQhu1Qn36DlJ3m6uoqpWQYhv5jwhHasneT1tf0U1rIVNCoVFnNnl4VKIkHRPbW7NXe2r3aW7NXe2r3aF/tPkXsSJPbOQ2ncjw56pPWJxE8qR9IyXJlEUgBAAAAugCBEgAAAKCb81V8rUDYUobh1TGDunfAISfDpTNGF+j1tWG5fE5Zzlrt27dJ6jumq6vWarZt60DgQDQQEuslsqdmj/bW7q0LiNiSrLAUCUlWSC7bUj9HhvqaHskwVG6HVW4FVWmFFHY4tS/s177afVKKeIjDcCSCJl63Vx6nRx7TI4/DI7fDLbfDnbjfKM10E2QBAAAA2qjDAiULFizQAw88oNLSUk2YMEGPPvqovvOd73TU7gAAAIBeq+zANoXClgwjV2MGepvMF4lEVFZWpt27d2vPnj3avXt34n5GRoYGDRqkoqIiFRUVadCgQcrKyuqQ+k46Kl9rtper+kCBKv1fat/+TR2yn9YIWSHVhmvlD/uTbmvCNaoN16o2VCt/qFq1wSrVhqpVG6pWVbBS4bBfigQlKyQ7HJQdCcmOhOSMhJVnmeprS/m2U3lyK09uZdsO2XatLFsyDck0DZmGIVu2asyIKhWWz4io0umUz2HIZxoqlyWfIoqYTu13uLTfdEmmU0oEPoxYYMVocL9unSFTLodLHodbbodHGa4sZaXlKsvtVZY7S5muTGW5spTlzlKWK/rYafJ/cwAAAIAkGbZt2+1d6PPPP6+LLrpIv//973XiiSfqoYce0osvvqgvvvhC/fv3b3Zbn8+nnJwc7d6zS16vV/ELANN0RG/rXSwk7hvxPMxNj97BikQUCgcVtsIKhYIKRYIKhcOKREIKh6OPTdOUw3DK6XDJ4XTJYUYvjp1OtxwOp5ymU06HW6bD7NH/XWhZlsLhsCKRiMLhcPS4HQ45nU45HI4efWwAgO7Ltm1ZtiVLlmzLki1blhVW4tQ5dmvLarRd3X0rlieRonA4oECoJvqjfLBGgVCtaoO1qg1WqcZfI3+wVrXBagXCAflDtQqE/QqGA9pXtU3VtfvVPzJSR2d/NxEAabiUlZXJspLr1Jzc3NykwEnD+4MGDVJGRkabnsNdFbV65tX52hZ4V32y81SYO0KmYUqGIVOmTCN6jmIapkzDkCFDhuGQkXhsxtKM6HlRJKxwJKSIFVY4ElbEii5hKyI7dhuxworYliw7LMuyFLEjClkh+e2gIlZYsq3o62JZkiwZtiXZlgzblmxbtmzZtmTLVuxPpi1lRRzKDpvyhh3KDkfvZ0SidY6/tgHbpRrbo+rYUmO75TRsZTrCyjRDsSUgjxGJHp9pRI/VNOQ0TRmmFHTaqnXYqnFa8jtshWxLYcNW0LYUMmwFZSkoS2FZCspWyLAUitY6Xo3E+81Q7DLJMCXDIZkO2aZThumQYTolh1NpjjRlujKV6cpSlscrb1quctJzlen2ynQ45DScMk2HHIZDpsMhh2HKYTpj56EuORyxx4ZDTtMhhyOaN7r32GfEsmMVq6unbdvRz45tK/GRse1Yuh19TYxomixLdmwbxV4v27ISZSqW37Ys2batsB1ROBxROP6eiIQVsSIKWdG0iBVROGJF3zu2JcuKKGxZilgRWbYtt9sjjytNHle6PG6PPK50paVlKM2TLo8rQ263R25nmlxOj5zONLkcHjkczuh7GwAOY3bi3KjuzCf+Pd4RrNj3vmWFFImEZFvRcwQrEpJlR2TEzyXM6LlF9H7s3MJ0xs43zMQ6mWZSmqQu+b3Btm2FIrZqQxEFQhHVxpdgRP6QJX8oIn8ooppgRP5wND0QtlQbjCgYsWLf/0b0nzVit4ZhyIg9dpix31Dr5YmvczlNeZym0lwOeWL3PU6H0lzRW48reX2ayyGn2XGv8eHItm1ZVkRWJKSIFYzdhmTF3+vx97gVkmVbMk2nHKZbDqdbDsMp0+GSw+GWw+GSaUbvx39Tx8F1SKDkxBNP1AknnKDHHntMUvSHzqKiIl111VX61a9+lZQ3EAgoEAgkHvt8PhUVFelHd4+XO61tL6Qdu2Axki4TVP8qNaH+R9luIk+qfIZd735sn42ahZY+sy3J1yhPikaonfZnGPWyJHZjJz9ujzawQT2MesmG3dLDafq1bW6fdZdurdSubb8dLc+U7PgbyoimteeHsu4lNBLPa3R3TTx3Tey8uc9HYk9G3WfOaCZzo4t5u/6pVOISuuVPROxL2TDi+44n171gKT6hrVD/3SnJTl1W8gnhQdgN34UtqV/j94dR/zNqNKhnE0U2fFqNJlY2lS/Vy5JYZx/sSJpfW/8bqdnPaTNtdas+P214Wxiq95ltUIyt1M9T0neI0Ti9TZ/5Q2go6n+HHbz8Fj5JjRqJRt+wTZfXpu/Blm/aYkbbSmx4RHYT6akk7a09v2OaqlR7lNUJmmhm26v06E38a9Bo8LlucK5jJLZp30rFf4hv7Yvz7QeWPvj7umbzGIah/Px8FRQUJJZ+/fqppqZG27dv144dO7R9+3b5fL4W7TMvL0+FhYWxCzdLkUikxbfDJo7V2LP8DY6y3S8/WsUVMeSyJJclOevdd0Ukl2XE0iWFTSnsVDjgVFXEraqIS5WxpSriki/kUGXIUGXQUFVICkVs2ZGwZIVlR8KyrYgkQ4bTLcPhlOF0yXC45HQ4lOUxleUMK9sRUpYjpEwzHL11RNMyzbDSHWGZsSsb04jdypZhRK9z4umSJduQbNNWxIzeWg5bIdNW0CkFHLaCDingkIIOO3Hb/OcsxUqjqYctb+3sxkkp83VvdY1D3fe4ETtNa1nj29Tz0Oias8kSWsJu9mHrpXjFW3gic9C3Woua2FRnW60vu70v5xo/GS3XsuuT5p9cO+lO47y23Shnk7Wpu2skpac8642f8ycCnbG19R83WNe4Go2vR+0UDY1tG0nZjdg1c3Sp990dSzeMePsWvU38vmHYSScYdafzRlJdE69L/UYyfu3W8Fyx/mMj9t1e77CMBtdvDX/zqP+0NLwma+6KOrGumd9UGjXySU93w29lI5E91fd10lNVv94teAsnfp+rd53Y8GWTjLr1iZfXiJ2b2UpeGX9dO/b7ouH7Pum9Xf9k0U6Rp9H2qXbQ3JNnJH8O23Ko9SpjqOF3fvJ7OVVVjIZZE8Wmbn+N+D9FKMXnuK6olIkde97fOodel9a+WIbqfge0kz8o7cmWZMdaN9tQ9P+8jLoLoUSmRK0aXMHbTR9ZE/VN/ops38/rM//9abuWV1+797UOBoNavXq15s6dm0gzTVNTp07V8uXLG+WfN2+e7rjjjkbpViSiSNPzIPYMLX1zd6NGoVMdrsfdnPonlKlWJ84TYj/b1DtJbNXTaTRxv8dKPoHrCZfY3VaveD8AQOs093ut3XyuFqnftJqWIYctmZYky5RhGTIsQ4pIikRv7YghOywZlinbMqSIIdsyFap1qHJPtc4444ykIEhBQYH69++fFBRxOg9+mu/z+ZICJw3vb9++XdXV1dq/f7/279/fpmNf++FqDcg9Vjn5DhmGnfjByjBjpzCGZBtG3brED12x5y1+8WhJtmXEOxjEHtuJNFnREyLbjl70xW8lU1bYlGk5ZIedkmWq1nDINhyy5ZAczugQV6ZDMqLBDNN0yOF0ye1yyuUw5cyUHEZ0GK0sU8o1o/+B6TBjPWJMU2aK+4ZhKBwOq7a2Vn6/X36/P3a/OnY/oJpAUGXBsHaEIvIHwwqEwgpGbAVCEYUtO1FmXdnxXihmvXUOORxmUnr8P0kdhiW3actl2rHb6GO3w5LbLbkzbLnSJKdHcqTZMj2W5JYMp5W4fo6+cZU450z84GjU+4El9muXLTX6kSH5Qrv+j2zJPwLaibQUJyN2wzyNy44+Nup+U2v041vsRyzbjp5LxzY2Yx2wjHhHLDO6sWFG7xux4JMMO/Y82LKMaF+yhu1Dm1qJZn4La9dz2nb88adRvQ5S9kGPoyPPP3tq2a3VVF06qo7d6dibc5Dr60ba4zq5Qz/IDRzkt/b2KjLF76ipH7dlJ/UCNC0urpnMSe29lAia1d1vZx3920rsR6DoP9jEFlmxXpa26npWRpeGvXka17VhJZMDkXVBmuhtNF4TP2mLngwkpcXRYaF9NPfejp+71L+VEudg8SBxw3/qbFyQeL2a0e6BkrKyMkUiERUUFCSlFxQUaOPGjY3yz507V9ddd13icbxHyYxRc5Ue6+Yf/6A3jKfbdv3u24lU1U9ocn2DaJgdfzdaSo7cNnqX2vVS6s7Yk0qrP+RC/RPKRp13bDWMyRlG43x1wzY0861kx+vVIE/DshLPR4qyjNhwB4o2qIaZuJe4jXffj9434/8bkGiA7XrPb4MaJtWnYW2Tm+qGr1n90uxGkc3Gr3HDclINP9Hwkq1ph9KF0FDzXfCj3fQdMk2HTJlymC45TEcszSnTdMoZX286m62LbVmK2OFoFz1Z0WCjHZEVW2zLig7UYNmy7HBsowavmR29TXp14r09Yl/AlmyZSc9Z8n8U2PHeHY2ejLr3icMRHT7LNB0yTFMuhynD4ZDD4ZTDiN53mqZMhxl9jkxDztjwYrYtRSLhxJBc0e62liIRS1YkXHffisTSQopYduy/WyNJz2Gji+z4SUaq5zcpT7M5Gj1K9ao115mv4efYlh3/BMowoi2hKUO2YUYHJTHrPqdS3bCERuyzqoZNWjP7qVfB5Eol/hUrdRvauPjGa+PHbKheO5piNym3T6q31bKTaNuKncylFu2KHa1R8r4lNfHZNerd2oaZ8rU14ieWiccNt26qbTDU8CPe1Pskeaif+POazDroe7WJ1zDVCxNbYSt+fI2/v5qqa7xedtI2Dba1Gn9nJH8VN/XOaunlVONXKsXLEntopsgU//ZpvEnTrXLqy8zUzXjTv2409xXUcFXiP/YabJS6iJb/ohL/MTtF7dTUf4ylLrp136eJz1s7d+WPv48Nw5QZH9op/t0bG5LBjJ0TyTCj30uGKdNwxtpbs975S/3PesveGQ6nS2meNLndLnk8bnlcTrndbjmdTrlcrsRt/AfvjuT1ejVmzBiNGZN6knXbtlVRUaHt27dr7969ieEv23IbHyLXMIw23dYfdrP+wlAPbWfbdmJo01SLHR8Kq8FiWZYs246dU0ZvI1Yk6XE0UBM9bzHMWJCnyTQzkWYYhhymGbsuiZ3dGIYMM/7jjJFYJ6lBuiQZcjoccsTKd8Tmh3E66t6DLRU/1wyHwwqFQgqGQgoGQ9FgVjCoQMCvQKBWgUB0CD1/oDpp+4bf/8k9nZWUFj+/a5hWP2/yuuTHSfdTrat31dagconvdaned3niJvmaOppHSXkbXos1KiPpoi75OsOOP7aT99/wtu4aMmlH0WBV8u7q7S/FdW8T+2vs4O2KVS9X/LdHo9592fVf5walt7DZaqp6tm3FPgdmYrgcxYcxNGP3E21u8ufLjP2qEr1GtmVFosMRWrGh7SwrukSHsovms61IYkjHSMSKBrbN+Pdp7JjM+EEbMmTVBcBj+zPMWC+A+PWmbUWD3fWuOeuGP4y9S2JtTt2wfFbsPM0RnffJNiXFr3dMGYoO+xfNE1tnR7/jo0HURAUU/6HZiA2xGa1DdJjN6K1VVx/ZshWJlRsrP/H5jO4nfm4RvVKO1cmMR5jr2p7odV2q34UUHQoy6XWv+7Q0Cg/bdZ/OeHr8Mx//lBn1PjeNWwMlUo16OeqaLaPBNnbDHPWOqt6IDinSbNV97utuY+8BxX9jUPR90yBP9Md4R/QczHBEn1vDIUMuyXTIMJySzNitQ4oN39nkuUHivVV/+EYrtsf4a27JNOrlieU3EnWK1tOIDyFZ73iMRDl1u4y/Kg0qkrp+sXWmKbnM6D9pNL4+TH7c+HrAUnTYsNg5ruLXrfHfE8zY5yX6/rVj7UI8CJL41Np23WdD8Wvw5LZbtqWQZStk24pYtkJWfKydaD6zQc1S1TZxHLad+C2poTaf6TW4Lm+XMlu1Zerrx0a5DFOGnFLifR4dLjWa5oimGc7GbwZFX4tG73c7+pufYUdk22EZRnSIWssOy7YjshWW7LBsK5yizOTzhaTfz5LahXiuFHWS0aBORuyvcUvRtKY/I1bK33jbT5fP3ufxeOTxeBqlT592dmyOEgAAAACdzTAM5ebmKjc3t6urgg5gGIacTmeLeh8djkzTlNvtltvt7uqqAAAAoBO0+4xzffv2lcPh0O7du5PSd+/ercLCwvbeHQAAAAAAAAAAQJu1e6DE7Xbr+OOP17vvvptIsyxL7777riZPntzeuwMAAAAAAAAAAGizDulnfd111+niiy/WxIkT9Z3vfEcPPfSQqqurNXv27I7YHQAAAAAAAAAAQJt0SKDkRz/6kfbu3atbb71VpaWlOvbYY/Xmm282muA9lfiEaj6fryOqBgAAAAAAAAAAeqDs7OzGE9m3A8OORya6ia+++kpHH310V1cDAAAAAAAAAAB0I3v27FG/fv3avdwO6VFyKPLy8iRJ27ZtU05OThfXBgAOnc/nU1FRkbZv3y6v19vV1QGAQ0a7BqC3oV0D0NvQrgHobeLtmtvt7pDyu12gxDSj88vn5OTQkAPoVbxeL+0agF6Fdg1Ab0O7BqC3oV0D0Nt0xLBbkmR2SKkAAAAAAAAAAAA9AIESAAAAAAAAAABw2Op2gRKPx6PbbrtNHo+nq6sCAO2Cdg1Ab0O7BqC3oV0D0NvQrgHobTq6XTNs27Y7pGQAAAAAAAAAAIBurtv1KAEAAAAAAAAAAOgsBEoAAAAAAAAAAMBhi0AJAAAAAAAAAAA4bBEoAQAAAAAAAAAAhy0CJQAAAAAAAAAA4LDVJYGSe++9V4Zh6Nprr02k+f1+zZkzR/n5+crKytKMGTO0e/fupO22bdumH/zgB8rIyFD//v31y1/+UuFwuJNrDwCNNWzX9u/fr6uuukrHHHOM0tPTNXjwYF199dWqqKhI2o52DUB3lep8Lc62bU2fPl2GYejVV19NWke7BqC7aqpdW758uU477TRlZmbK6/Xq5JNPVm1tbWL9/v37NWvWLHm9XuXm5urSSy9VVVVVJ9ceABpL1a6VlpbqJz/5iQoLC5WZmanjjjtOL7/8ctJ2tGsAuovbb79dhmEkLSNHjkys78yYgfOQj6aVPvnkEz3xxBMaP358UvovfvEL/d///Z9efPFF5eTk6Morr9R5552nZcuWSZIikYh+8IMfqLCwUB999JF27dqliy66SC6XS/fcc09nHwYAJKRq13bu3KmdO3dq/vz5Gj16tL755hv913/9l3bu3KmXXnpJEu0agO6rqfO1uIceekiGYTRKp10D0F011a4tX75cZ555pubOnatHH31UTqdTa9eulWnW/U/hrFmztGvXLi1ZskShUEizZ8/WZZddpueee66zDwMAEppq1y666CKVl5frtddeU9++ffXcc8/pggsu0KpVq1RcXCyJdg1A9zJmzBi98847icdOZ13IolNjBnYnqqystIcPH24vWbLEPuWUU+xrrrnGtm3bLi8vt10ul/3iiy8m8m7YsMGWZC9fvty2bdv++9//bpumaZeWlibyLFy40PZ6vXYgEOjMwwCAhKbatVReeOEF2+1226FQyLZt2jUA3dPB2rU1a9bYRxxxhL1r1y5bkr148eLEOto1AN1Rc+3aiSeeaN98881Nbrt+/Xpbkv3JJ58k0t544w3bMAz722+/7chqA0CTmmvXMjMz7WeeeSYpf15env2HP/zBtm3aNQDdy2233WZPmDAh5brOjhl06tBbc+bM0Q9+8ANNnTo1KX316tUKhUJJ6SNHjtTgwYO1fPlySdH/9Bk3bpwKCgoSeaZNmyafz6fPP/+8cw4AABpoql1LpaKiQl6vNxEZp10D0B01167V1NToxz/+sRYsWKDCwsJG62nXAHRHTbVre/bs0cqVK9W/f3+ddNJJKigo0CmnnKJ//vOfiTzLly9Xbm6uJk6cmEibOnWqTNPUypUrO+0YAKC+5s7XTjrpJD3//PPav3+/LMvSokWL5Pf7deqpp0qiXQPQ/WzevFkDBw7UUUcdpVmzZmnbtm2SOj9m0GlDby1atEiffvqpPvnkk0brSktL5Xa7lZubm5ReUFCg0tLSRJ76BxxfH18HAJ2tuXatobKyMt1111267LLLEmm0awC6m4O1a7/4xS900kkn6Zxzzkm5nnYNQHfTXLv21VdfSYqOjT1//nwde+yxeuaZZ3T66adr3bp1Gj58uEpLS9W/f/+k7ZxOp/Ly8mjXAHSJg52vvfDCC/rRj36k/Px8OZ1OZWRkaPHixRo2bJgk0a4B6FZOPPFEPf300zrmmGO0a9cu3XHHHfq3f/s3rVu3rtNjBp0SKNm+fbuuueYaLVmyRGlpaZ2xSwDoUK1p13w+n37wgx9o9OjRuv322zunggDQSgdr11577TX94x//0Jo1a7qgdgDQegdr1yzLkiRdfvnlmj17tiSpuLhY7777rp588knNmzevU+sLAAfTkuvQW265ReXl5XrnnXfUt29fvfrqq7rgggv04Ycfaty4cZ1cYwBo3vTp0xP3x48frxNPPFFDhgzRCy+8oPT09E6tS6cMvbV69Wrt2bNHxx13nJxOp5xOpz744AM98sgjcjqdKigoUDAYVHl5edJ2u3fvTgzrUFhY2GhG+/jjVEM/AEBHOli7FolEJEmVlZU688wzlZ2drcWLF8vlciXKoF0D0J0crF1bsmSJtmzZotzc3MR6SZoxY0ZiKAfaNQDdSUuuQyVp9OjRSduNGjUqMeRDYWGh9uzZk7Q+HA5r//79tGsAOt3B2rUtW7boscce05NPPqnTTz9dEyZM0G233aaJEydqwYIFkmjXAHRvubm5GjFihL788ksVFhZ2asygUwIlp59+uj777DOVlJQklokTJ2rWrFmJ+y6XS++++25imy+++ELbtm3T5MmTJUmTJ0/WZ599ltSYL1myRF6vt9GJLQB0tIO1aw6HQz6fT2eccYbcbrdee+21Rv/xQ7sGoDs5WLt200036V//+lfSekl68MEH9dRTT0miXQPQvRysXTvqqKM0cOBAffHFF0nbbdq0SUOGDJEUbdfKy8u1evXqxPp//OMfsixLJ554YqceDwAcrF2rqamRJJlm8s99Docj0YuOdg1Ad1ZVVaUtW7ZowIABOv744zs1ZtApQ29lZ2dr7NixSWmZmZnKz89PpF966aW67rrrlJeXJ6/Xq6uuukqTJ0/WpEmTJElnnHGGRo8erZ/85Ce6//77VVpaqptvvllz5syRx+PpjMMAgISDtWvxIElNTY3+8pe/yOfzyefzSZL69esnh8NBuwagW2nJ+Vqq/8gZPHiwhg4dKonzNQDdS0vatV/+8pe67bbbNGHCBB177LH685//rI0bN+qll16SFO1dcuaZZ+pnP/uZfv/73ysUCunKK6/UhRdeqIEDB3b6MQE4vB2sXQuFQho2bJguv/xyzZ8/X/n5+Xr11Ve1ZMkSvf7665Jo1wB0LzfccIPOPvtsDRkyRDt37tRtt90mh8OhmTNnKicnp1NjBp02mfvBPPjggzJNUzNmzFAgENC0adP0+OOPJ9Y7HA69/vrruuKKKzR58mRlZmbq4osv1p133tmFtQaA1D799FOtXLlSkhKT5sVt3bpVRx55JO0agF6Hdg1AT3PttdfK7/frF7/4hfbv368JEyZoyZIlOvrooxN5nn32WV155ZU6/fTTE9esjzzySBfWGgBSc7lc+vvf/65f/epXOvvss1VVVaVhw4bpz3/+s84666xEPto1AN3Fjh07NHPmTO3bt0/9+vXTd7/7Xa1YsUL9+vWT1LkxA8O2bbvdjgwAAAAAAAAAAKAH6ZQ5SgAAAAAAAAAAALojAiUAAAAAAAAAAOCwRaAEAAAAAAAAAAActgiUAAAAAAAAAACAwxaBEgAAAAAAAAAAcNgiUAIAAAAAAAAAAA5bBEoAAAAAAAAAAMBhi0AJAAAAAAAAAAA4bBEoAQAAAAAAAAAAhy0CJQAAAAAAAAAA4LBFoAQAAAAAAAAAABy2/j/KLWNp/Jd+IQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 2000x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax=utils.plot_track([just_tss_pred], label='No CRE', color='k', alpha=1)\n",
    "for i, p in enumerate(preds_per_iter):\n",
    "    utils.plot_track([p], zoom=[400, 500], ax=ax, label=f\"{i+1} CRE copies\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f25d81dc-793a-4692-a6e8-15050ab459f1",
   "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}