{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7219000f-db8e-4038-813b-9931173ffb99",
   "metadata": {},
   "source": [
    "# Sufficiency test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0e520279-ba5f-4872-a38d-27487ed26c56",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3083d60-240d-4c05-afc2-02bded899c81",
   "metadata": {},
   "source": [
    "Load Enformer and example sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "id": "5e3979fd-a872-41fa-afc4-109ee830b141",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = '../../../data'\n",
    "track_index = [5111]\n",
    "model = custom_model.Enformer(track_index=track_index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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 = 'GATA2_chr3_128487916_-'\n",
    "\n",
    "gene_name, chrom, start, strand = gene.split('_')\n",
    "wt_seq = seq_parser.extract_seq_centered(chrom, int(start), '-', model.seq_length)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "61c69d44-b67f-4f71-b10b-5dc6c7b867f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TSS bin indeces\n",
    "bins = [447, 448]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "502ab1fc-bcb5-444d-ab7b-8351db4f7463",
   "metadata": {},
   "outputs": [],
   "source": [
    "wt = model.predict(wt_seq)[0,:,0]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "34adb807-c057-4993-aa1a-dc2bf161d288",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: >"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABkYAAADFCAYAAAABv06xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA84UlEQVR4nO3de3yU9Z33//c1xxxnQoAcOATwLIoUQTHVtlapiLQ/rba3dllLq62/uuh66O1pV63auljd3VqtxXbbaveultbeq64+6gGx4loRAaUKKIKgoOQAhGSSSTKn63v/Ea/LTJhAAoE55PV8POYBc13XXPOdyTdX5ro+8/l8LGOMEQAAAAAAAAAAwDDgyfYAAAAAAAAAAAAADhUCIwAAAAAAAAAAYNggMAIAAAAAAAAAAIYNAiMAAAAAAAAAAGDYIDACAAAAAAAAAACGDQIjAAAAAAAAAABg2CAwAgAAAAAAAAAAho28DIwYYxSJRGSMyfZQAAAAAAAAAABAHsnLwEh7e7vC4bDa29uzPRQAAAAAAAAAAJBH8jIwAgAAAAAAAAAAsD8IjAAAAAAAAAAAgGGDwAgAAAAAAAAAABg2CIwAAAAAAAAAAIBhg8AIAAAAAAAAAAAYNgiMAAAAAAAAAACAYYPACAAAAIC8Eo1GZVlW2i0ajWZ7WAAAAADyBIERAAAAAAAAAAAwbBAYAQAAAAAAAAAAwwaBEQAAAAAAAAAAMGwQGAEAAAAAAAAAAMMGgREAAAAAAAAAADBsEBgBAAAAAAAAAADDhi/bAwAAAACAwfD7/TrnO+doe9t2HVd1nI4YdYT8fn+2hwUAAAAgTxAYAQAAAJBXAoGAzvnuOXqr6S2dd8x5mnPknGwPCQAAAEAeoZQWAAAAgLyTMilJUjwVz/JIAAAAAOQbAiMAAAAA8opt2/r4/Y+164Nd2vjORq1bt062bWd7WAAAAADyhGWMMdkexGBFIhGFw2G1tbUpFAplezgAAAAADqFoNKqysrK0ZR0dHSotLc3SiAAAAADkEzJGAAAAAAAAAADAsEFgBAAAAAAAAAAADBsERgAAAAAAAAAAwLBBYAQAAAAAAAAAAAwbgwqMLFq0SCeccIJCoZBCoZDq6+v1zDPPuOu7u7u1YMECjRw5UmVlZbrgggvU1NSUto+tW7dq7ty5KikpUVVVla677jolk8mheTUAAAAAAAAAAAB7MajAyLhx43TXXXdp9erVWrVqlc444wyde+65WrdunSTpmmuu0VNPPaXHHntMy5Yt0/bt23X++ee7j0+lUpo7d67i8bheffVV/fa3v9XDDz+sW2+9dWhfFQAAAAAAAAAAQAaWMcYcyA4qKyt1zz336Gtf+5pGjx6tRx99VF/72tckSe+++66OPfZYLV++XKeccoqeeeYZffnLX9b27dtVXV0tSXrwwQd1ww03aMeOHQoEAgN6zkgkonA4rLa2NoVCoQMZPgAAAIA8E4/Hdfo3T1d7d7uKfEU6fdLpuvPOOwd8PgEAAABgeNvvHiOpVEqLFy9WNBpVfX29Vq9erUQioVmzZrnbHHPMMaqrq9Py5cslScuXL9eUKVPcoIgkzZ49W5FIxM06ySQWiykSiaTdAAAAAAxPgUBAs/7/WfrspZ/V5777Od1zzz0ERQAAAAAM2KADI2+//bbKysoUDAb1ve99T48//rgmT56sxsZGBQIBVVRUpG1fXV2txsZGSVJjY2NaUMRZ76zrz8KFCxUOh93b+PHjBztsAAAAAAXENrYkKZ6KZ3kkAAAAAPLNoAMjRx99tNasWaMVK1bo8ssv1/z587V+/fqDMTbXTTfdpLa2Nve2bdu2g/p8AAAAAHKXbdtqaWhRpCGiXR/t0uYtm2XbdraHBQAAACBP+Ab7gEAgoCOOOEKSNH36dK1cuVI//elPdeGFFyoej6u1tTUta6SpqUk1NTWSpJqaGr3++utp+2tqanLX9ScYDCoYDA52qAAAAAAKUFdXlxb9/SL3/mItVkdHh0pLS7M4KgAAAAD5Yr97jDhs21YsFtP06dPl9/u1dOlSd92GDRu0detW1dfXS5Lq6+v19ttvq7m52d1myZIlCoVCmjx58oEOBQAAAAAAAAAAYK8GlTFy0003ac6cOaqrq1N7e7seffRRvfTSS3ruuecUDod16aWX6tprr1VlZaVCoZCuvPJK1dfX65RTTpEknXXWWZo8ebIuvvhi3X333WpsbNTNN9+sBQsWkBECAAAAAAAAAAAOukEFRpqbm/XNb35TDQ0NCofDOuGEE/Tcc8/pS1/6kiTpJz/5iTwejy644ALFYjHNnj1bP//5z93He71ePf3007r88stVX1+v0tJSzZ8/X3fcccfQvioAAAAAAAAAAIAMLGOMyfYgBisSiSgcDqutrU2hUCjbwwEAAABwCEWjUZWVlaUto8cIAAAAgIE64B4jAAAAAAAAAAAA+YLACAAAAAAAAAAAGDYG1WMEAAAAALLN5/PpuHOOUyqVkmVZOq76OPl8nNoAAAAAGBjOHgAAAADklWAwqNMuP01Ou8T5n5mvYDCY5VEBAAAAyBeU0gIAAACQV4wxblBEkmLJWBZHAwAAACDfkDECAAAAIK+k7JS62rrc4Ehzc7PMRCPLsrI8MgAAAAD5wDK9v2qVJyKRiMLhsNra2hQKhbI9HAAAAACHUGukVSPCI9KWdXR0qLS0NEsjAgAAAJBPKKUFAAAAIK/Yxs72EAAAAADkMQIjAAAAAPJKyqSyPQQAAAAAeYzACAAAAIC8kofVgAEAAADkEAIjAAAAAPIKpbQAAAAAHAgCIwAAAADyCqW0AAAAABwIAiMAAAAA8gqltAAAAAAcCF+2BwAAAAAAg2F5LR11xlGybVuWLJUVlcnn49QGAAAAwMBw9gAAAAAgr/gDfp1+9enu/dryWgWDwewNCAAAAEBeoZQWAAAAgLzSt5RWPBXP0kgAAAAA5CMyRgAAAADklaSdVKI74QZIoiYqY4wsy8ryyAAAAADkA8vkYefCSCSicDistrY2hUKhbA8HAAAAwCH0XsN7OnrM0WnLOjo6VFpamqURAQAAAMgnlNICAAAAkFdsY++xLA+/7wUAAAAgSwiMAAAAAMgrmQIjCTuRhZEAAAAAyEcERgAAAADklUyBkVgyloWRAAAAAMhHBEYAAAAA5JVMgZF4Kp6FkQAAAADIRwRGAAAAAOSVTIGRpJ3MwkgAAAAA5CMCIwAAAADySqbASKZlAAAAAJCJL9sDAAAAAIDBsDyWJn12koq8RUraSSVNUsYy2R4WAAAAgDxBYAQAAABAXvEH/frSjV/ShIoJ6oh3aFfnLvkCnNoAAAAAGBhKaQEAAAA5jN4Ze3LKZlmy5LW8kqSUSWVzSAAAAADyCIERAAAAIEf9detfddUzV2ld87psDyWnOIERj+WRz9OTKUIACQAAAMBAkW8OAAAA5KhNLZuUtJPavHuzjqs6LtvDyRnRaFS//P9+mbbsO5u/I43K0oAAAAAA5JVBZYwsXLhQJ510ksrLy1VVVaXzzjtPGzZsSNumu7tbCxYs0MiRI1VWVqYLLrhATU1Nadts3bpVc+fOVUlJiaqqqnTdddcpmeQbXgAAAEBvThZEPBXP8khyi5Mx0hultAAAAAAM1KACI8uWLdOCBQv02muvacmSJUokEjrrrLMUjUbdba655ho99dRTeuyxx7Rs2TJt375d559/vrs+lUpp7ty5isfjevXVV/Xb3/5WDz/8sG699dahe1UAAABAASAwklnK3jMIkilYAgAAAACZDKqU1rPPPpt2/+GHH1ZVVZVWr16tz3/+82pra9Ovf/1rPfroozrjjDMkSQ899JCOPfZYvfbaazrllFP0/PPPa/369XrhhRdUXV2tz3zmM/rhD3+oG264QbfddpsCgcDQvToAAAAgjxVqYOTdne8qGo9q+pjp+/V4I7PHMjJGAAAAAAzUATVfb2trkyRVVlZKklavXq1EIqFZs2a52xxzzDGqq6vT8uXLJUnLly/XlClTVF1d7W4ze/ZsRSIRrVuXualkLBZTJBJJuwEAAACFrlADI79c/Uv9xxv/oWg8uu+NM8iUMULzdQAAAAADtd+BEdu2dfXVV+vUU0/V8ccfL0lqbGxUIBBQRUVF2rbV1dVqbGx0t+kdFHHWO+syWbhwocLhsHsbP378/g4bAAAAyBvOxf5YKpblkQwdY4yi8aiMMepKdu3XPjL2GMkQLAEAAACATPY7MLJgwQKtXbtWixcvHsrxZHTTTTepra3NvW3btu2gPycAAACQbYWYMdI7qLG/wQxKaQEAAAA4EIPqMeK44oor9PTTT+vll1/WuHHj3OU1NTWKx+NqbW1NyxppampSTU2Nu83rr7+etr+mpiZ3XSbBYFDBYHB/hgoAAADkrYSdkFRYgZHeJa/2u/yVJY2fMV7l/nJZshRJRCRriAYIAAAAoOANKjBijNGVV16pxx9/XC+99JImTZqUtn769Ony+/1aunSpLrjgAknShg0btHXrVtXX10uS6uvrdeedd6q5uVlVVVWSpCVLligUCmny5MlD8ZoAAACAguBkVBRSYKR3Zsf+BkZ8AZ/m3DpHJ409SR7LoxUfrZDHf0DtEwEAAAAMI4MKjCxYsECPPvqonnzySZWXl7s9QcLhsIqLixUOh3XppZfq2muvVWVlpUKhkK688krV19frlFNOkSSdddZZmjx5si6++GLdfffdamxs1M0336wFCxaQFQIAAAD0UoiltHqXz9rf8ldOOS6P5ZHX8h7QvgAAAAAMP4MKjCxatEiSdPrpp6ctf+ihh/Stb31LkvSTn/xEHo9HF1xwgWKxmGbPnq2f//zn7rZer1dPP/20Lr/8ctXX16u0tFTz58/XHXfccWCvBAAAACgwTimtWLJwmq/3zhLZ3x4jaYERj/eA9gUAAABg+Bl0Ka19KSoq0gMPPKAHHnig320mTJigP//5z4N5agAAAGDYoZRWZp3RTv3m67/RQ3pIlmXJNrbOePWMoRoiAAAAgAK3X83XAQAAABx8hVhKKy1j5ABKaSVj6UEVSmkBAAAAGCg6FAIAAAA5yimlZRu7YEpF9Q6M7G/GiFNKq7dCeX8AAAAAHHwERgAAAIAc1ftif6FkjaQ1Xz/AHiNp+yVjBAAAAMAAERgBAAAAcpBt7LQAQCxVGA3Yh6LHSKbAyP7uCwAAAMDwQ2AEAAAAyEF9L/QXSsbIUPUYGcgyAAAAAMiEwAgAAACQgwo1MNK7fNaQ9hihlBYAAACAAfJlewAAAAAA9lSwgRFz4D1GZEm1x9eqIlghn8ennV07yRgBAAAAMGAERgAAAIAcVKiBkaEopeUL+vSVf/mKvnzUlxUuCuuRtx6RN+gdqiECAAAAKHCU0gIAAABy0HAIjBxoKS3LsuS1vAe0LwAAAADDD4ERAAAAIAcVamBkKHqMOJkmXssrr8e7x34BAAAAYG8opQUAAADkoIINjAxBj5Huzm7959//pxZ7FsuyLMVSMd3yX7cM1RABAAAAFDgCIwAAAEAOKtTAyFD0GEmZlLoj3epWt7uM5usAAAAABopSWgAAAEAOKtTAyFCU0soUBKGUFgAAAICBIjACAAAA5KC+QYNYMpalkQyttIyR/QxmGGP2WLa/2ScAAAAAhh8CIwAAAEAOKtiMkd49Rg6glNYey8gYAQAAADBABEYAAACAHFSogZHer2soS2nt774AAAAADD8ERgAAAIAcVKiBkYPWY4RSWgAAAAAGyJftAQAAAADYU6EGRoakx4iMRh8xWiOKRijgDagp2iRbewZLAAAAACATAiMAAABADirUwEjvzI79zRjxBX366r9/VZdMu0Tjw+N1+0u3yxMgGR4AAADAwHD2AAAAAOSggg2M2AfefN0ppeWxPPJaXkn0GAEAAAAwcARGAAAAgBzkXOgv9hdLKpzAyFCU0koLjHi8B7QvAAAAAMMPpbQAAACAHOQEEEr8JepKdBVkYGR/szy6u7r16Hce1X/7/ls+r0+tsVZd9POLhmqIAAAAAAocgREAAAAgB/UOjOzSLsVSsSyPaGgMRY8R27bV0dyhDnV8ul87JWOMLMs64DECAAAAKGyU0gIAAABykBM0KPWXSiqcUlpD2WNkoMsBAAAAoDcCIwAAAEAO6p0xIhVOYGQoe4zsbd8AAAAA0B8CIwAAAEAOKtTm60NSSqufwMj+ZqAAAAAAGF4IjAAAAAA5qG8prZSd2u8Mi1xyMEtpFcL7AwAAAODgIzACAAAA5KC+pbSkwsgaGUwprc27N+v+FferqaMpbTkZIwAAAAAOhC/bAwAAAACwJyeAUOQrkmVZMsYonoq7pbXyVe/AyL5Kaf1161+1tnmtJo2YpC8f9WV3uW1sjRg/QpVFlQr4AmqMNsqyLHqMAAAAABiQQWeMvPzyy/rKV76iMWPGyLIsPfHEE2nrjTG69dZbVVtbq+LiYs2aNUsbN25M26alpUXz5s1TKBRSRUWFLr30UnV0dBzQCwEAAAAKiXOR3+fxKeANSCqMjJHB9BjpTnZLkmLJWNpyX5FPX3/g63rh9Re0fv16ffMX35Qv6KOUFgAAAIABGXRgJBqNaurUqXrggQcyrr/77rt133336cEHH9SKFStUWlqq2bNnq7u7291m3rx5WrdunZYsWaKnn35aL7/8si677LL9fxUAAABAgXGCBn6vv7ACI4PoMeK83r6v29mHx+o5nfFa3gHtDwAAAACk/SilNWfOHM2ZMyfjOmOM7r33Xt18880699xzJUn/+Z//qerqaj3xxBO66KKL9M477+jZZ5/VypUrNWPGDEnS/fffr3POOUf/+q//qjFjxhzAywEAAAAKQyKVkNSTMRL0BtWudsVSsX08KvcNpseImzHS53U7PUacwIjP4xvQ/gAAAABAGuLm61u2bFFjY6NmzZrlLguHw5o5c6aWL18uSVq+fLkqKircoIgkzZo1Sx6PRytWrMi431gspkgkknYDAAAACpmT/eC1vG5fka5EVzaHNCR6Z3XYxu63kbr0aaZI31Ja8e64HlvwmL540hc1efJk/ea7v1EyliRjBAAAAMCADGnz9cbGRklSdXV12vLq6mp3XWNjo6qqqtIH4fOpsrLS3aavhQsX6vbbbx/KoQIAAAA5rXcprWJfT2CkM9GZzSENib59RVJ2Sh5v5u9rOZkimUpp7d62W7u1211mjKH5OgAAAIABGdKMkYPlpptuUltbm3vbtm1btocEAAAAHFS9S2mV+EskSV3JAsgY6VPuam9ZHk6mSN/ASH9ZJpTSAgAAADAQQxoYqampkSQ1NTWlLW9qanLX1dTUqLm5OW19MplUS0uLu01fwWBQoVAo7QYAAAAUst6ltJzASKFmjPTHyRjpr8fIvvYNAAAAAJkMaWBk0qRJqqmp0dKlS91lkUhEK1asUH19vSSpvr5era2tWr16tbvNiy++KNu2NXPmzKEcDgAAAJC30kpp+QunlFbfDJG9BTMGnTFCjxEAAAAAAzDoHiMdHR3atGmTe3/Lli1as2aNKisrVVdXp6uvvlo/+tGPdOSRR2rSpEm65ZZbNGbMGJ133nmSpGOPPVZnn322vvvd7+rBBx9UIpHQFVdcoYsuukhjxowZshcGAAAA5LOMpbQKoPl630BIf4ER29juOkppAQAAABhKgw6MrFq1Sl/84hfd+9dee60kaf78+Xr44Yd1/fXXKxqN6rLLLlNra6tOO+00PfvssyoqKnIf88gjj+iKK67QmWeeKY/HowsuuED33XffELwcAAAAoDA42Q8+j69gmq/bxpYxJm1Zf1keTrZI3/87+8mEjBEAAAAAAzHowMjpp5++x8lMb5Zl6Y477tAdd9zR7zaVlZV69NFHB/vUAAAAwLDhZEsUUvP13hkdPo9PSTvZb8ZI774ivTNGjDGyLEtlVWWqCFbIY3nUHm+XZVlkjAAAAAAYkCHtMQIAAABgaGQqpZXvGSO9MzqCvmDPsn6CGb2zROKpuPvlrJRJyRf06e9+9Xfa8P4Gffjhh7rjyTvkC/povg4AAABgQAiMAAAAADnGGOOWi/J5fG7z9XzvMdI7cBHwBiT1X/6qb18R537v7HWP1XM647W8e90XAAAAAPRGYAQAAADIMb0DCAWVMfJJdojH8riBkYGU0pI+DYz0Dn44ARGvx5u2fwAAAADYm0H3GAEAAABwcPUNjDjN1/O9x4jzurwerxvU6DcwkswcGLGNrWQsqaduekorb1spS5Zault0xu1nHNSMEWOMjIybpQIAAAAgfxEYAQAAAHJM72CB1/K6GSOJVEKJVEJ+rz9bQzsgTuDCa3n3meXRN2PEuW+MkTFGOzbt0A7tcNcbYw5ajxFjjO559R51Jjp18+dvls/DaRQAAACQz/i6EwAAAJBjnAv8Po9PlmWpyFcky7Ik5XfWSO/X5QQX+svy6C9jZG9ZIQerlFZ7vF3vt7yvhvYGfRz5+KA8BwAAAIBDh8AIAAAAkGN6l5ySJMuyPi2nlccN2J3ARe/AyGB7jDhN6TPu/yCV0mpob3D//1Hko4PyHAAAAAAOHQIjAAAAQI7pnVnhKPb3BEbyuQG7W0qrV4+Rfktp7aXHSH8OVimtxo5G9//bItsOynMAAAAAOHQIjAAAAAA5xrnA7/d82kvEyRjJ58CImwnTu8dIP1keTiDE4QRK9poxcpBKafUOjJAxAgAAAOQ/AiMAAABAjulbSkuS24C90HqMDLSUlnM/G6W0+gZGjDEH5XkAAAAAHBq+fW8CAAAA4FAq2FJa9p6ltPoNjOyjlFZxqNgNFiXsRNr+h1pDx6c9RroSXWrpatHIkpEH5bkAAAAAHHxkjAAAAAA5JlMpLTdjJJ+br5tPm6+7pbT66zGyl+br/iK/Lv/j5dq5c6d27typxSsWy1/kPyg9RmLJmHZ37ZYkjSgeIYlyWgAAAEC+IzACAAAA5JhMGSNOYCSfM0Z69xhxXlt/5a+cjBGP5Um772SMOMsl7XNfB6Ip2iRJKg+W6+iRR0uiATsAAACQ7wiMAAAAADkmU48Rp/l6PvcYcbJDBtNjpCxQJmnPUlq9AyP7yj45EA3tPWW0aspqND48XhIZIwAAAEC+o8cIAAAAkGOcrBAnGCIVRsaIk9HRu8dIv6W0PskQKQ+WKxKLpAVGkrGkFt+8WK+GXpVlWWqPt+vE/33iQckYcRqv15TVaEz5mLRlAAAAAPITgREAAAAgx0RiEUlSKBhylznN1/O5x0jvEmFOlkd/GSNOIKQ8UC7p0wwS29gyxmjrW1u1VVvd7aeZaQPKGHm76W1tbNmoc48+Ny0jpz+9AyMjinp6jLR2t+7zcQAAAAByF4ERAAAAIMc4gZFwUdhdVggZI4PqMfJJIMQJDvUtpbW3/e/NH9f9Uc3RZo0PjddJY0/a5/ZOYKS2rFYVRRWSeoJTsWRMQV9wn48HAAAAkHvoMQIAAADkmEwZI05gpBB6jPQupdVvj5FepbR6399bYGRfpbRsY2tn505J0t+a/rbP8drGVnO0WVJPxkiRr8gNhrTF2vb5eAAAAAC5icAIAAAAkGOci+5ppbQ+6TeSzxkjTuCid/P1fnuMfJIx4pTSGoqMkd1du93Hr2tet8/SWzs7dyppJ+X3+lVZXCnLshQO9mTxtHUTGAEAAADyFYERAAAAIMe4pbSCn5bSKguUSZLaY+17DQ7ksoH2GDHGuIGQvqW09hbM2FegY1fXLvf/nYlObWrZtNftnTJa1aXVsixLktxyWvQZAQAAAPIXgREAAAAgxzjZCL0zRsJFYXk9XtnGztuL8m4prX30GEnYCRljJPUqpfVJBomR6X//+yil5ZTRcrzV9NZet3f7i5TXusucvi/5+jMAAAAAQGAEAAAAyCnxVFzdyW5J6c3XPZZHI4tHStrzAn++SMsY+aTHSKYsD6efiPRppkzfjBF/kV/FxcUqKSlRcUlxv/vqbVdnT8aI877212fEGCPb2G5gpKasxl1HxggAAACQ/wiMAAAAADnEKaPl9/oV9AbT1o0sKYzAiNfzacZIplJaTnZIwBtQka9IUnqPEX+RXz964Ufq7OxUNBrV37b+Tf4i/4AzRurH1cuyLO2I7lB7rD1tm2g8quuXXK8HXn9ADe0NkjIHRmi+DgAAAOQvAiMAAABADundX8Tpa+EYVTJK0qeZD/nGCVx4Le9ee4w4QZCgL6iANyDp0ywSp5SWx/r0VCbo6wkgReNRtwRXJk6PkXGhcRpdMlqStC2yLW2b93a9p0gsorXNa/VB6weS0gMjTt8XMkYAAACA/EVgBAAAAMghmfqLOPK9lJZT6iqtlFaGLA+nlFjQ+2lgJJ6Kyxjj7qN3YKS2rFYey6OOeIdaulr6fX7nfRtZMlJ14TpJ0ra29MDIh20fuv+3jS3LslRdWu0uczNGuskYAQAAAPKVL9sDAAAAAPApN2OkV38Rh5Mxkq+Bkd49Rtzm6xn6gjgZIwFvIK2cWMJOyDa2kvGkfvW/f6Xny5+XZVnyeDw6+ZqTtb1ru7a0bnFLjvV9bifLY1TJKI0LjdOq7av2yBj5sPXDtPsji0fK7/W793s3XzfG7JHVAwAAACD3ERgBAAAAcogTGMmUMZLvgRG3lNY+eox0JjolpZfSknrKadnGlrGN3nn1Hb2jd9x1X7v1az2Bkd1bNGPMjD32ubtrt4wx8nv9Kg+UZ8wYMca4GSMVRRVq7W5NK6PlLJd6gjfdyW4V+4sH/T4AAAAAyC5KaQEAAAA5xGnqvbfASFusLWNAIdc52SG9e4xkKqW14qMVkqS6cJ0sy3IzNuKpuGxjZ9z3xBETJUmbd2/OuN7pLzKyeKQsy9L48HhJUlO0ye1f0tLVomg8Kq/Hq0tPvFSVxZU6ZdwpafsJeAMq8ZdIos8IAAAAkK8IjAAAAAA5ZG8ZI2WBMgW8ARlj8rIBe8JOSErvMdI3wNMcbdbfmv4mSTpj0hmSlNZnJFMgRZImVkyUJG1t25q2T9vYev3j17V6+2pJcstshYIhhYIhGWP0cfvHkj7tLzK2fKyOGnmUFs5aqJPGnrTHcznltJwgFgAAAID8QiktAAAA5LyOeIf+58P/0artq1RVWqVvTv1mwZYwcpp6h4N79hixLEujSkZpe/t27ezcqeqy6j22OZiMMdrYslE1ZTUKBUPa1LJJ63esV1VplcaWj1Vtea1bIquvNxve1LrmdZKk8mC5yoPlkqSG9gZ90PqBJlZMVNJO6vn3n5cxRlOqp7hlrIp8RYrGo1rdsNrNJulrdMlolQZKFY1H9VHkI02smKiUndJv3vyNVm1f5W7nZN1IPRkpa5vXalvbNh024jC3v8iEigl7fR8qiirU0N4w5Bkjxhgt+3CZNuzcoLZYm+rH1etzEz43pM8BAAAAIMuBkQceeED33HOPGhsbNXXqVN1///06+eSTszkkAACAvbKNrTca3lBduE5VpVXZHs6wEIlFdNcrd7kZEh9FPlJDR4O+N+N7e/R/2B8pO6XmaLPa4+2aEJ6goO/TZt9diS7t7t4tSaotqz0kjbb3ljEiyQ2MOKWhEqmE3tn5jjbs3KDqsmp9dvxn+w1O7G+zcGOM1jav1ZMbntS2tm0q9hdr5tiZWvbhMhlj3O08lkdlgTKVBcpUF65Tka9I7+16Ty1dLepOdkuSTqw9UVOqpshjeXRi7Yl6o+EN/XL1LzWqZJQ27trolsqaddgsd7+fHf9ZPbXhKT214SlJn/b56M2yLE2smKh1zev05LtPalTJKG1p3aJtbdvk9XhVXVqttlibptVMcx8zPjxea5vX6tVtrypcFNb6HeslSRPCew+MOEGrD1o/0MyxMzO+p4lUQpZl9fuzcKTslLZFtimRSui595/T201vu+s2796s0aWjdcyoY7S7a7c2tmzUlt1b9FHkIxkZlQXKdPrE01UeKNczm55RTVmNzj7i7H0+Z19JO6lEKpEx2JiyU1q3Y53eanpLrd2tKvIVqdhXrBJ/iarLqjUuNE61ZbVpDeptYyueissYoyJfUb/vTyQWUbgoPOjxAgAAAAcqa59A//CHP+jaa6/Vgw8+qJkzZ+ree+/V7NmztWHDBlVVcZEBAAAcPMYY2cZ2exz0FUvGtLFloxKphEoDpbLUc3GzLFCm36/9vdY1r1ORr0gLTl6go0YedcjGnbST6kx0KuANqMhXJGOMYqmYovGojIxK/CUq8hXJY3mUSCXUmeiUx/Ko2F/sXnjsTnZrZ+dOtXS1aFfnLu3q2uX+m0glVD++XqeOP1VST9mjSCyit5reUmNHo2rKajS6ZLSKfEUKF4VVUVShRCoh29iqKKpwAwpJO6l4Ku6OJWWnZFmWPFZ6FdeUnVIsFVN3sluJVEJBX1DFvmIFvAFtbduqt5vfVqm/VCu3r9Suzl0aVTJKX5z0RS15f4ka2ht0+0u363MTPqfaslq3SXdrd6t2du7U6JLRCvqCWtu8VvFUXBPCEzSieIQkaXv7dkXjUVUWV6qho0Hrd6xXItVT4snr8aouXKeKogo1dTRpe/t2d7xHjjxSX5jwBVUUVag8WK5Sf6n7mBJ/iYwxao+3q6mjSW2xNvfn1Pfm8/jUEe9wL5zviO5Qc7TZ/RntKzDilIJa27xWrd2teumDlxSNR931f974Zzdg1JXoUneyW13JLnUlupSwEyr1l6rYX6xEKqESf4lqymo0pnyMQsGQtrdvV1usTZYslQXKNKJ4hCKxiNbvWK+mjiZJPQGIrkSXXvrgJUnS5NGTlbAT+jjysToTnYrEIorEImnvnfO4z0/4vC46/iJ3Llw89WJ92PZhzxz8JPBV7C/WtJppOnrk0e5j5x45Vyk7pT9v/LN8Hp+++5nv6if6yR7vzZGVR2pd8zo3wCFJfq9f35vxPR1fdfwe2x898mg9s/EZfdD6gRatXOQud8py9cfJ1PnLlr/oraa3VF1arTHlY3TYiMO0pnGN1jSuUTwVl8fyaGLFRBX7ixWNR3V81fE6re40bW3b6ga31jSuUXusPW285xx5jj6KfKTV21frF6t+odJAqXZEd2Qcy5sNb8qyLDdAtaZxjc6YdIbGh8bLY3nUEe/Qjs4dbnkxj+VRsa9YR1QeobXNa/X8+89rR+cOGWN0Yu2JOufIc1RTViO/16+G9gb96o1f6aPIR3t9PyzLUrGvWAk7oaSdTAuWST1zdnxovBJ2Qru7dqu1u1WdiU5JcoM7k0dPVk1ZjQLegGxjK5aKyTa2PJbHPaa1dLUoZVIKeAMKent+5wPegIK+oPwe/yEJXAIAkM/cfm/9nAcBg9Gd7NZ7u95TR7xDtrE1sWKixpaPzfiZrD3Wrq5kl0aVjNrjvDBbLNP3U+shMnPmTJ100kn62c9+JkmybVvjx4/XlVdeqRtvvDFt21gsplgs5t6PRCIaP368bnv2NhWVFh3ScaPwGGXlVwDDjCVO1AfLyCiRSiiWiskYI6/H614c6X3b3/f2QH73s/Sns+e5+4y791gG8pr29X4N9THReb6+H4ycb4Q74zcyGV+LMSZtTAPZfl/3bWMrmojKGKPSQGnPt5llybIs999dnbsG1Nja6/EqHAzLNrZsY8vIuP93bl7L6wYIknZSKZPq+ddOuRfUpZ6TlKSddC8sWrLci34ey6OuZJfbIFqSgr6gG5RIe88tS36PX/FUPG1ZRVGFYsmYezHyYPBYHlmW5Z5wWZYlr+VNuyAb8Abk9XgVS8b6fY97X+R1lPhLdONpN6q6rFqt3a36P3/7P1rbvHbIxu4EZTKVRioNlO51vM74UiaV9jM6ED6PT/eefW/at/AdL255UX9Y+4e0ZSOKR2jy6Mla17zuoDUEL/IV6fMTPq9Zh83SXz74i17Z+oq+dNiXdNbhZ7k/s7ZYmzriHWrtbtX7Le8rnorriMoj3MBLpoyEbW3b9OSGJ3XYiMN00piTNKpkVL8XuN/Z8Y7Kg+Ua4R2hsrKytHUdHR3yBX3667a/uoGG6rJqHVl5pBtMyuSD1g/0Px/+j97f/b5GFo/UsaOP1ZmTztzrRXYnu+OFzS+oK9E1kLdvr0oDpSr1lyoUDOl/Hfe/NKFighKphO565S43KGFZlurCdTqi8gjVhesU8Aa0YecGvfzhy7KNrSnVU7R59+a0INmB8Hl8SpmUe6w8acxJGh8er+5kt7qT3WqPtauho0EfRT7a7+f0WJ49jmH7yzneOEFY5zNC2mcGy8r4OcJZd7A/rxVC4CbfP9Me7J/B/r4/AxlXf58/M31OOtDnP9Q/5/4+fw7m9Qz2M2ymfQ/FPvrdNo9fo9MfLOANyO/xu32/3M/zfT6X955PvefSoTjOYu+6kl3uF12qy6pV7MtOWdpD/fdwoL9nmX4/DuSxgxnPQB+fS6/FGKOdnTv3OEcK+oIq8hUp6A0q6AvKGKPORKdaulok9RxLRhSPcM9zpQP7++k81jkWnX3E2ZpWO20fj/zk8dkIjMTjcZWUlOhPf/qTzjvvPHf5/Pnz1draqieffDJt+9tuu0233377Hvv51uJvKVASONjDBQAAw9DIkpEKB8PqTHTKyCieiqu1u1WjS0br29O+refff15vNrx5yMeVKWjg8/jksTxpwZD+tpV6LsKOKhmlyuJKjSweqZElIzWyeKTa4+16dtOz7jfTvR6vAt6Ajqw8UoeNOEyNHY1q7W5Vd7Jbrd2tblaEpAMKCPg8Pvk8PsVTcfck2+/1a0rVFMVSMbV0tegbx39DR486Ou1x63es1+rtq92gUSwVU3mgXKNKRqmxo1HRRFSTR09WWaBMW9u2KhqPKmVSqi6tVigYUktXi0LBkKbWTNXY8rGSehp/f9z+sdq62xQuCuuokUepLFCm3V279fz7z+uD1g/UEe9Qe7w94wVxy7JUWVypyuJKJVIJ9yKyc+u9nc/jc7Ntastq3cBaZXGlplZP1dSaqRnfr454h/7v+v+rtlibvJZXp4w7RdNqp7nfql+/Y726k90yMir2FavYX+z+6/f41RHvUHeyW36v3724vb19uyKxiGrLat0gQiQW0e6u3SoPlqu2rFbTaqepyMeXknrrTnZry+4taulq0fu739fm3ZtVF67TGZPOUE1ZjTriHdrUssnNmlry/hJtb9+uqtIqHTbiMI0oHqHDRxyuyaMnZ/zmZlt3m1776DWNDY3V4SMOzxhY2tW5S52JTo0Pj1dbd5te3PKiNrVsUlO058JHsa9Yo0tHK+jtyegyMmrpatHWtq0q8Zfoy0d9WTPGzFB7rF1PbnhS7+58N+33efLoyfrWZ77lNpzvywmIdSe73d/l3hfOupJd2t6+XR9HPlaxv1gVRRXurchXpDca3tDybcv1UeSjtKCec6LrHBM8lkcjikfI7/ErloopnorvM2AJAACAg2t06WhVl1YrZVLavHtzv+eFzvmPk6l/sPzdlL/TFyZ+YUDbZiUwsn37do0dO1avvvqq6uvr3eXXX3+9li1bphUr0hsq9pcx8sp7r6isPP2bYsD+KIRvbyF3DeQwa2T49kwGPo/PvQjX91v4zm2o5cPxoL8sjL5zyFnfdw4eqvnWX6ZH70yf3t/wcMbb3/3+tne22df2zn2P5VGpv1Rej1eRWMS9GO+M0Ta2ygJlqimr2eM9Ttkp95vFxhhtb9+uhJ3Y41vJvb+ZbBtbXYkuGRn5PD55LW/Pvx6vWxpLknsx0bmw6ARj4qm4UnZKJf4SlQZKVewrdoM0AW9ApYFSt4SMs79YMqayQJn7+9Meb9fOzp0q8hVpZPHItB4ae/zcTM/z+r3+faY49+5X4Vz4N8a4ZW26k91K2An3gmw8FVfCTihlp3q+RfRJ+SunhJTz3F3JLpX4S9wLq7nKNrYsWYqn4trVtUsey6NRJaP22t8jluq5kFviL8mZFHIcOs4cyIUAU3eyW36Pf4+AjDFGXcmeEmyWerLNDtXfxpSdco+pAW9gj79jmcbh9DOJJWNKmdSnGXyfZAhmyuRz1vVdPxxkI+t1uGToH+g3//f1+P6+hZ9pfX/b7e1nMZCMlMEazOfNwR5nBvoeDHYf/W7bz74L/TU65WSlTz/HORc1ne0zfZ7pm/XtLDsUOLfun9/rV21ZrSSpoaPhoF+gziRbfxMGOicG+jt2oL+fBzKeXHot4WBYVaVV7vZJO6mWrhb3S2tOkCTgDWhcaJyCvqB2RHcoEou41UEc+zM3+h5XLMvS2PKxe80W7y0vutwFg0EFg3uewE+pnqJQKHPtZQAAgH0pCwzuCxa9LyBalqWxobFDPaQBCfqCbo+D3nweX09fij4fm0LBUL/9KvqyLGuvgZO+2zqc/hm99S0DVarSAT33QJ8/25wLAUFfUGPKx+xze8uycuKCOLInl+ZAf+OwLEsl/hK3zN+h5PV4M2bO7O3k3GN5Mh5/AABA/5zee8BQ8nl8qirde+/w6rLqjOey2ZCVr6mNGjVKXq9XTU1NacubmppUU1OTjSEBAAAAAAAAAIBhICuBkUAgoOnTp2vp0qXuMtu2tXTp0rTSWgAAAAAAAAAAAEMpa6W0rr32Ws2fP18zZszQySefrHvvvVfRaFTf/va3szUkAAAAAAAAAABQ4LIWGLnwwgu1Y8cO3XrrrWpsbNRnPvMZPfvss6qu3neNMaexSiQSOdjDBAAAAAAAAAAAeaK8vHyfTeQt07d9ex7YvHmzDj/88GwPAwAAAAAAAAAA5JC2tjaFQqG9bpO1jJEDUVlZKUnaunWrwuFwlkcD7FskEtH48eO1bdu2ff5SArmAOYt8w5xFvmHOIt8wZ5FvmLPIN8xZ5BvmLHJZeXn5PrfJy8CIx9PTMz4cDvOLh7wSCoWYs8grzFnkG+Ys8g1zFvmGOYt8w5xFvmHOIt8wZ5GvPNkeAAAAAAAAAAAAwKFCYAQAAAAAAAAAAAwbeRkYCQaD+sEPfqBgMJjtoQADwpxFvmHOIt8wZ5FvmLPIN8xZ5BvmLPINcxb5hjmLfGcZY0y2BwEAAAAAAAAAAHAo5GXGCAAAAAAAAAAAwP4gMAIAAAAAAAAAAIYNAiMAAAAAAAAAAGDYIDACAAAAAAAAAACGDQIjAAAAAAAAAABg2MjLwMgDDzygiRMnqqioSDNnztTrr7+e7SFhGHr55Zf1la98RWPGjJFlWXriiSfS1htjdOutt6q2tlbFxcWaNWuWNm7cmLZNS0uL5s2bp1AopIqKCl166aXq6Og4hK8Cw8nChQt10kknqby8XFVVVTrvvPO0YcOGtG26u7u1YMECjRw5UmVlZbrgggvU1NSUts3WrVs1d+5clZSUqKqqStddd52SyeShfCkYJhYtWqQTTjhBoVBIoVBI9fX1euaZZ9z1zFfkurvuukuWZenqq692lzFvkUtuu+02WZaVdjvmmGPc9cxX5KKPP/5Yf//3f6+RI0equLhYU6ZM0apVq9z1nIchl0ycOHGP46xlWVqwYIEkjrPIPalUSrfccosmTZqk4uJiHX744frhD38oY4y7DcdZFIq8C4z84Q9/0LXXXqsf/OAHeuONNzR16lTNnj1bzc3N2R4ahploNKqpU6fqgQceyLj+7rvv1n333acHH3xQK1asUGlpqWbPnq3u7m53m3nz5mndunVasmSJnn76ab388su67LLLDtVLwDCzbNkyLViwQK+99pqWLFmiRCKhs846S9Fo1N3mmmuu0VNPPaXHHntMy5Yt0/bt23X++ee761OplObOnat4PK5XX31Vv/3tb/Xwww/r1ltvzcZLQoEbN26c7rrrLq1evVqrVq3SGWecoXPPPVfr1q2TxHxFblu5cqV+8Ytf6IQTTkhbzrxFrjnuuOPU0NDg3l555RV3HfMVuWb37t069dRT5ff79cwzz2j9+vX6t3/7N40YMcLdhvMw5JKVK1emHWOXLFkiSfr6178uieMscs+Pf/xjLVq0SD/72c/0zjvv6Mc//rHuvvtu3X///e42HGdRMEyeOfnkk82CBQvc+6lUyowZM8YsXLgwi6PCcCfJPP744+5927ZNTU2Nueeee9xlra2tJhgMmt///vfGGGPWr19vJJmVK1e62zzzzDPGsizz8ccfH7KxY/hqbm42ksyyZcuMMT1z1O/3m8cee8zd5p133jGSzPLly40xxvz5z382Ho/HNDY2utssWrTIhEIhE4vFDu0LwLA0YsQI86tf/Yr5ipzW3t5ujjzySLNkyRLzhS98wVx11VXGGI6zyD0/+MEPzNSpUzOuY74iF91www3mtNNO63c952HIdVdddZU5/PDDjW3bHGeRk+bOnWsuueSStGXnn3++mTdvnjGG4ywKS15ljMTjca1evVqzZs1yl3k8Hs2aNUvLly/P4siAdFu2bFFjY2PaXA2Hw5o5c6Y7V5cvX66KigrNmDHD3WbWrFnyeDxasWLFIR8zhp+2tjZJUmVlpSRp9erVSiQSafP2mGOOUV1dXdq8nTJliqqrq91tZs+erUgk4n6LHzgYUqmUFi9erGg0qvr6euYrctqCBQs0d+7ctPkpcZxFbtq4caPGjBmjww47TPPmzdPWrVslMV+Rm/77v/9bM2bM0Ne//nVVVVVp2rRp+o//+A93PedhyGXxeFy/+93vdMkll8iyLI6zyEmf/exntXTpUr333nuSpL/97W965ZVXNGfOHEkcZ1FYfNkewGDs3LlTqVQq7Q+CJFVXV+vdd9/N0qiAPTU2NkpSxrnqrGtsbFRVVVXaep/Pp8rKSncb4GCxbVtXX321Tj31VB1//PGSeuZkIBBQRUVF2rZ9522mee2sA4ba22+/rfr6enV3d6usrEyPP/64Jk+erDVr1jBfkZMWL16sN954QytXrtxjHcdZ5JqZM2fq4Ycf1tFHH62Ghgbdfvvt+tznPqe1a9cyX5GTNm/erEWLFunaa6/VP/3TP2nlypX6x3/8RwUCAc2fP5/zMOS0J554Qq2trfrWt74lic8FyE033nijIpGIjjnmGHm9XqVSKd15552aN2+eJK53obDkVWAEADA0FixYoLVr16bVEQdy0dFHH601a9aora1Nf/rTnzR//nwtW7Ys28MCMtq2bZuuuuoqLVmyREVFRdkeDrBPzrc/JemEE07QzJkzNWHCBP3xj39UcXFxFkcGZGbbtmbMmKF/+Zd/kSRNmzZNa9eu1YMPPqj58+dneXTA3v3617/WnDlzNGbMmGwPBejXH//4Rz3yyCN69NFHddxxx2nNmjW6+uqrNWbMGI6zKDh5VUpr1KhR8nq9ampqSlve1NSkmpqaLI0K2JMzH/c2V2tqatTc3Jy2PplMqqWlhfmMg+qKK67Q008/rb/85S8aN26cu7ympkbxeFytra1p2/edt5nmtbMOGGqBQEBHHHGEpk+froULF2rq1Kn66U9/ynxFTlq9erWam5t14oknyufzyefzadmyZbrvvvvk8/lUXV3NvEVOq6io0FFHHaVNmzZxnEVOqq2t1eTJk9OWHXvssW4JOM7DkKs+/PBDvfDCC/rOd77jLuM4i1x03XXX6cYbb9RFF12kKVOm6OKLL9Y111yjhQsXSuI4i8KSV4GRQCCg6dOna+nSpe4y27a1dOlS1dfXZ3FkQLpJkyappqYmba5GIhGtWLHCnav19fVqbW3V6tWr3W1efPFF2batmTNnHvIxo/AZY3TFFVfo8ccf14svvqhJkyalrZ8+fbr8fn/avN2wYYO2bt2aNm/ffvvttA85S5YsUSgU2uMkFTgYbNtWLBZjviInnXnmmXr77be1Zs0a9zZjxgzNmzfP/T/zFrmso6ND77//vmpraznOIiedeuqp2rBhQ9qy9957TxMmTJDEeRhy10MPPaSqqirNnTvXXcZxFrmos7NTHk/65WKv1yvbtiVxnEWByXb398FavHixCQaD5uGHHzbr1683l112mamoqDCNjY3ZHhqGmfb2dvPmm2+aN99800gy//7v/27efPNN8+GHHxpjjLnrrrtMRUWFefLJJ81bb71lzj33XDNp0iTT1dXl7uPss88206ZNMytWrDCvvPKKOfLII803vvGNbL0kFLjLL7/chMNh89JLL5mGhgb31tnZ6W7zve99z9TV1ZkXX3zRrFq1ytTX15v6+np3fTKZNMcff7w566yzzJo1a8yzzz5rRo8ebW666aZsvCQUuBtvvNEsW7bMbNmyxbz11lvmxhtvNJZlmeeff94Yw3xFfvjCF75grrrqKvc+8xa55Pvf/7556aWXzJYtW8xf//pXM2vWLDNq1CjT3NxsjGG+Ive8/vrrxufzmTvvvNNs3LjRPPLII6akpMT87ne/c7fhPAy5JpVKmbq6OnPDDTfssY7jLHLN/PnzzdixY83TTz9ttmzZYv7rv/7LjBo1ylx//fXuNhxnUSjyLjBijDH333+/qaurM4FAwJx88snmtddey/aQMAz95S9/MZL2uM2fP98YY4xt2+aWW24x1dXVJhgMmjPPPNNs2LAhbR+7du0y3/jGN0xZWZkJhULm29/+tmlvb8/Cq8FwkGm+SjIPPfSQu01XV5f5h3/4BzNixAhTUlJivvrVr5qGhoa0/XzwwQdmzpw5pri42IwaNcp8//vfN4lE4hC/GgwHl1xyiZkwYYIJBAJm9OjR5swzz3SDIsYwX5Ef+gZGmLfIJRdeeKGpra01gUDAjB071lx44YVm06ZN7nrmK3LRU089ZY4//ngTDAbNMcccY375y1+mrec8DLnmueeeM5L2mIfGcJxF7olEIuaqq64ydXV1pqioyBx22GHmn//5n00sFnO34TiLQmEZY0xWUlUAAAAAAAAAAAAOsbzqMQIAAAAAAAAAAHAgCIwAAAAAAAAAAIBhg8AIAAAAAAAAAAAYNgiMAAAAAAAAAACAYYPACAAAAAAAAAAAGDYIjAAAAAAAAAAAgGGDwAgAAAAAAAAAABg2CIwAAAAAAAAAAIBhg8AIAAAAAAAAAAAYNgiMAAAAAAAAAACAYYPACAAAAAAAAAAAGDb+HzdwMBTqtt39AAAAAElFTkSuQmCC",
      "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": [
    "**Sufficiency test**\n",
    "\n",
    "\n",
    "In this example we will test the necessity of (5Kb) tiles in the context of an enhancing context sequence of GATA2 gene. \n",
    "\n",
    "To run the test we need:\n",
    "- a loaded model\n",
    "- onehot encoded sequence (WT) of the sequence\n",
    "- TSS tile coordinates (or any other fixed tile coordinate that is always embedded in the shuffled backgrounds)\n",
    "- a list of tile coordinates to test\n",
    "- num_shuffle - number of shuffles to perform\n",
    "- tile sequence (onehot) which can be extracted based on coordinates (if tile_seq is None)\n",
    "- optionally, we can set mean=False to not average the shuffle results\n",
    "- optionally, we can set return_seqs to True to return the shuffled sequeneces for future use\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "fa851a1f-f266-4142-821b-333d08aa684b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enhancing tile at position 100804 - 105804\n",
      "TSS tile at center position 95804 - 100804\n"
     ]
    }
   ],
   "source": [
    "perturb_window = 5000\n",
    "N_shuffles = 10\n",
    "tss_tile, cre_tiles = utils.set_tile_range(model.seq_length, perturb_window)\n",
    "cre_tiles = cre_tiles[19:20]\n",
    "print(f'Enhancing tile at position {cre_tiles[0][0]} - {cre_tiles[0][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": [
    "Sufficiency of an enhancing tile of GATA2 gene TSS:\n",
    "\n",
    "\n",
    "As we can see, the GATA2 TSS activity is low on its own (control case of just TSS) compared to the high activity of the WT. \n",
    "However, when we embed the enhancing CRE together with the TSS the activity goes back up (not fully but to a much higher level than \n",
    "the control case)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "52086ff5-48a3-45b2-99ad-818cd054e529",
   "metadata": {},
   "outputs": [],
   "source": [
    "_, pred_mut, pred_control, control_sequences = creme.sufficiency_test(model, wt_seq, tss_tile, cre_tiles, N_shuffles, mean=False, return_seqs=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "0a03fe86-99c1-4f06-badc-4eebdfa9864a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABlMAAADFCAYAAADAOAQdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByd0lEQVR4nO3deZxcVZ03/s9dau3q6uq9ek1n63Q6IQlZyMKSBMIaETQqOsgmwqiAA47OTx4ZQZ0RX+gMDMqIzygEH2QYUEEWZQiBTiAkIQmELJ109vS+pffab937+6Nyb1d1V+9L9fJ5J/dV1Xc9tZ2qe773nK+gaZoGIiIiIiIiIiIiIiIiiktMdAGIiIiIiIiIiIiIiIgmMgZTiIiIiIiIiIiIiIiI+sFgChERERERERERERERUT8YTCEiIiIiIiIiIiIiIuoHgylERERERERERERERET9YDCFiIiIiIiIiIiIiIioHwymEBERERERERERERER9WNSBlM0TUNHRwc0TUt0UYiIiIiIiIiIiIiIaIqblMGUzs5OpKSkoLOzM9FFISIiIiIiIiIiIiKiKW5SBlOIiIiIiIiIiIiIiIjGC4MpRERERERERERERERE/WAwhYiIiIiIiIiIiIiIqB8MphAREREREREREREREfWDwRQiIiIiIiIiIiIiIqJ+yIkuABERERERERERERHRWAuHwwiFQokuBg2TyWSCJEkJOz6DKURERERE05zH44HD4QAAdHV1ISkpaVy2JSIiIiIaD5qmob6+Hm1tbYkuCo2Qy+WC2+2GIAjjfmwGU4iIiIiIiIiIiIhoytIDKVlZWbDb7QlpiKeR0TQNXq8XjY2NAICcnJxxLwODKUREREREREREREQ0JYXDYSOQkp6enuji0AjYbDYAQGNjI7KyssZ9yC8moCciIiIiIiIiIiKiKUnPkWK32xNcEhoN+uuYiNw3DKYQERERERERERER0ZTGob2mhkS+jgymEBERERERERERERER9YPBFCIiIiIiIiIiIiIion4wAT0RERER0TSnaRrWrl1r3B+vbYmIiIiIiCYL9kwhIiIiIprmHA4HysrKUFZWBofDMW7bEhERERFRfIIg9Ds98sgjAIBXXnkFq1atQkpKCpKTk7FgwQLcf//9xn7C4TB+9rOfoaSkBDabDWlpaVi5ciV++9vfjqh87733Hq677jqkp6fDbrejtLQU//iP/4iamhoAQFlZWUx5MzMzcd111+HgwYMx+7n99tvjPr5rrrlmROUbCwymEBERERERERERERFNIHV1dcb0xBNPwOl0xsz77ne/i61bt+Kmm27Cpk2b8NFHH2Hfvn3413/9V4RCIWM/P/rRj/D444/jJz/5CcrLy/Hee+/h7rvvRltbW5/HXrduHTZv3tzn8t/85jfYsGED3G43/vSnP6G8vBxPP/002tvb8W//9m8x61ZUVKCurg7/+7//i0AggI0bNyIYDMasc80118Q8trq6Ovz3f//3sJ63scRhvoiIiIiIpjlFUfDmm28CADZu3AhZHvxpwki2JSIiIiJKBE3TEAwHB15xDJglMwRBGHA9t9tt3E9JSYEgCDHzAOD111/HxRdfjO9973vGvOLiYtx4443G36+99hq+9a1v4Ytf/KIxb/HixcMuf3V1Nb797W/j29/+Nh5//HFjflFRES677LJeQZqsrCy4XC643W7cf//9+OxnP4ujR49i0aJFxjoWi6XXY5uIeKZDRERERDTNtbS0GCdcDQ0NyMrKGpdtiYiIiIgSIRgO4tt/+3ZCjv3ktU/CIltGZV9utxsvvPACDh06hIULF/a5zrvvvotvfetbyMzMHPExX375ZQSDQfzTP/1T3OUulyvu/Pb2drz44osAALPZPOJyJAKH+SIiIiIimuYCgUCii0BEREREREN03333YcWKFbjgggtQVFSEL3/5y3jmmWdift//+7//O5qamuB2u7Fo0SJ84xvfwN/+9rdhH/P48eNwOp3IyckZ1Pr5+flwOBxwuVx44YUX8NnPfhYlJSUx67zxxhtwOBwx009/+tNhl3GssGcKEREREdE0pyeJJCIiIiKaDsySGU9e+2TCjj1akpKS8Oabb+LkyZN47733sGvXLvzjP/4j/uM//gM7d+40EsMfOnQI+/btw44dO7B9+3Zcf/31uP32240k9D/96U9jghc+nw+7du3Cvffea8wrLy9HYWEhNE0b1DBluvfffx92ux27du3CT3/6Uzz99NO91lm/fj1+/etfx8xLS0sb6tMx5hhMISIiIiKa5nw+n3G/s7OTQ3URERER0ZQmCMKoDbU1EcyePRuzZ8/G17/+dfzgBz9AcXEx/ud//gd33HEHAEAURaxYsQIrVqzA/fffj+effx633HILfvCDH2DmzJn4xje+gS996UvG/m6++WZs2rQJn//85415ubm5ACI5Wdrb21FXVzeo3ikzZ86Ey+XCvHnz0NjYiJtuugnbt2+PWScpKQlz5swZjadiTHGYLyIiIiKiac7r9Rr3q6urE1gSIiIiIiIaiaKiItjtdng8nj7XKS0tBQBjnbS0NMyZM8eYbDYbsrKyYubJcqRfxhe+8AWYzWY89thjcffdMwF9tHvuuQeHDh3CK6+8MsxHl1jsmUJERERENM35/X7jfl1dHcLhMCRJSmCJiIiIiIhoII888gi8Xi+uu+46zJgxA21tbXjyyScRCoVw5ZVXAogEPy6++GKsWbMGbrcbp0+fxoMPPoji4uJeuUsGo6CgAI8//jjuvfdedHR04NZbb0VRURGqq6vx+9//Hg6HA//2b/8Wd1u73Y677roLDz/8MG688UZjuLBAIID6+vqYdWVZRkZGxpDLN5bYM4WIiIiIaBrTNC0mmBIKhXqdyBARERER0cSzdu1anDp1CrfeeitKSkpw7bXXor6+Hm+//TbmzZsHALj66qvx+uuv4/rrr0dxcTFuu+02lJSU4O233zZ6mwzVt771Lbz99tuoqanB5z73OZSUlODrX/86nE4nvvvd7/a77b333osjR47g5ZdfNua99dZbyMnJiZkuueSSYZVtLAmapmmJLsRQdXR0ICUlBe3t7XA6nYkuDhERERHRpOXz+fDaa6/hoYcegtVqxUMPPYTCwkKsXr16UNt3dXVh/fr1AID33nsPDodjLItLRERERDQkfr8fp0+fxsyZM2G1WhNdHBqhRL6e7JlCRERERDSNeb1e2O12PP3009i9ezeSkpLQ3Nwck0elPw6HA3v27MGePXsYSCEiIiIioimLwRQiIiIiomnM5/MBiIxfbLfbjXGJq6qqElksIiIiIiKiCYUJ6ImIiIiIpjGv1wtFUXD06FGEw2EUFhaiqakJVVVVKC4uNpJC9kVRFOzatQsAsGrVqmGPu0xERERERDSRDalnyq9//WssWrQITqcTTqcTq1evxt/+9jdjud/vxz333IP09HQ4HA5s2rQJDQ0NMfuorKzExo0bYbfbkZWVhe9973tQFGV0Hg0REREREQ2Jz+dDV1cXvva1r+HSSy+F2WyG2WyGz+dDU1PTgNu3tLTg0ksvxaWXXoqWlpZxKDEREREREdH4G1IwJT8/Hz/72c+wb98+7N27F5dffjluuOEGHD58GADwwAMP4PXXX8fLL7+Mbdu2oba2Fp///OeN7cPhMDZu3IhgMIgPP/wQzz33HDZv3owf/vCHo/uoiIiIiIhoUHrmRhFFEfn5+QAiF0IRERERERHREIMp119/Pa677jrMnTsXxcXF+Nd//Vc4HA7s2rUL7e3t+N3vfod///d/x+WXX45ly5bh2WefxYcffmh0+3/77bdRXl6O559/HkuWLMG1116Ln/zkJ3jqqacQDAbH5AESEREREVHf4iWaLywsBADU19cjEAiMd5GIiIiIiIgmnGEnoA+Hw3jxxRfh8XiwevVq7Nu3D6FQCBs2bDDWKSkpQWFhIXbu3AkA2LlzJy644AJkZ2cb61x99dXo6OgwerfEEwgE0NHRETMREREREdHIaJpmJKCPlpycjNTUVGiahurq6gSUjIiIiIiIaGIZcjDl4MGDcDgcsFgs+MY3voFXXnkFpaWlqK+vh9lshsvlilk/Ozsb9fX1ACJXtkUHUvTl+rK+PProo0hJSTGmgoKCoRabiIiIiIh68Pv90DQNotj7tED/zV1ZWQlN08a7aERERERERBPKkIMp8+bNw/79+7F7925885vfxG233Yby8vKxKJvhwQcfRHt7uzFVVVWN6fGIiIiIiKYDfYgvq9Xaa1leXh4kSUJXVxdaW1vHu2hEREREREQTypCDKWazGXPmzMGyZcvw6KOPYvHixfiP//gPuN1uBINBtLW1xazf0NAAt9sNAHC73WhoaOi1XF/WF4vFAqfTGTMREREREdHI6MEUm83Wa5ksy8jNzQXARPREREREROPt6aefRnJyMhRFMeZ1dXXBZDJh3bp1MeuWlZVBEIQBp7KysvF9EFPMsHOm6FRVRSAQwLJly2AymbB161ZjWUVFBSorK7F69WoAwOrVq3Hw4EE0NjYa62zZsgVOpxOlpaUjLQoREREREQ2Bni8lJSUFCxcuxMKFC2N6qeiJ6Gtra2NO4qJZrda42xIRERER0fCtX78eXV1d2Lt3rzHv/fffh9vtxu7du+H3+4357733HtxuN+rq6ozpS1/6Eq655pqYeWvWrEnEQ5ky5KGs/OCDD+Laa69FYWEhOjs78cILL6CsrAz/+7//i5SUFNx55534zne+g7S0NDidTtx3331YvXo1Vq1aBQC46qqrUFpailtuuQWPPfYY6uvr8dBDD+Gee+6BxWIZkwdIRERERETx6T1TsrKycPDgwV7LU1NT4XA40NXVhZqaGsyYMaPXOk6nM+62REREREQ0fPPmzUNOTg7KysqM9vWysjLccMMNePfdd7Fr1y6jh0pZWRmuvvrqmNGfbDYbAoFAvyNC0dAMqWdKY2Mjbr31VsybNw9XXHEF9uzZg//93//FlVdeCQB4/PHH8ZnPfAabNm3CZZddBrfbjT//+c/G9pIk4Y033oAkSVi9ejW++tWv4tZbb8WPf/zj0X1UREREREQ0oP6G+QIAQRCM3inMW0hEREREU4WmaVAUJSGTpmmDLuf69evx3nvvGX+/9957WLduHdauXWvM9/l82L17N9avXz/qzxPFGlLPlN/97nf9LrdarXjqqafw1FNP9bnOjBkz8Ne//nUohyUiIiIiojGgD/NltVpx5MgRAJEr4ESx+5qr/Px8HDlyBK2trejs7ERycnLMPlRVRUVFRdxtiYiIiIgmonA4jL/97W8JOfa1114LWR5cs/z69etx//33Q1EU+Hw+fPLJJ1i7di1CoRCefvppAMDOnTsRCAQYTBkHPNMhIiIiIpqGNE0zgilerxelpaUoLS1Fc3NzzHoWi8UYGiBeIvrm5uY+tyUiIiIiouFbt24dPB4P9uzZg/fffx/FxcXIzMzE2rVrjbwpZWVlmDVrltGjnMbOkHqmEBERERHR1ODz+aBpGkRRHDB/YWFhIerq6lBdXY358+ez9wkRERERTWqSJOHaa69N2LEHa86cOcjPz8d7772H1tZWrF27FgCQm5uLgoICfPjhh3jvvfdw+eWXj1VxKQqDKURERERE05DeK8Vms0EQhH7XzczMhNVqhd/vR319PXJzc8ejiEREREREY0IQhEEPtZVo69evR1lZGVpbW/G9733PmH/ZZZfhb3/7Gz766CN885vfTGAJpw9eUkZERERENA3pyeftdvuA6wqCgIKCAgDxh/oiIiIiIqKxsX79enzwwQfYv3+/0TMFANauXYvf/OY3CAaDzJcyThhMISIiIiKahoYSTAFgjMHc3NxsbEtERERERGNr/fr18Pl8mDNnDrKzs435a9euRWdnJ+bNm4ecnJwElnD6mBx9mYiIiIiIaFTpARGbzTao9e12OzIyMtDc3IyqqirMmzdvLItHREREREQAioqKoGlar/kzZsyIO1+3efPmMSzV9MSeKURERERE05CeM2WwPVOA7t4pVVVV/Z64ERERERERTTXsmUJERERENA1FD/MlSRJmz54NALBarX1u43a7YTKZ4PP50NTUhKysLFit1kFtS0RERERENJkxmEJERERENM2oqgq/3w8gEkyxWCw4ceLEgNtJkoT8/HycPn0aVVVVyMrKgtPpHNS2REREREREkxmH+SIiIiIimmb8fj80TYMkSTCbzUPaVh/qq76+HsFgcCyKR0RERERENOEwmEJERERENM1EJ58XBAGqqqKxsRGNjY1QVbXfbZ1OJ1wuF1RVRXV19ZC2JSIiIiIimqwYTCEiIiIimmai86UAQHNzM7Kzs5GdnY3m5uYBt9d7p1RWVqKpqWlI2xIREREREU1GDKYQEREREU0z0T1ThiM3NxeSJKGzsxPt7e2jWTQiIiIiIqIJicEUIiIiIqJpxufzAejumTJUJpMJOTk5AICamppRKxcREREREdFExWAKEREREdE003OYr+GITkRPREREREQ01TGYQkREREQ0zeg9U4Y7zBcApKWlISkpCYqijFaxiIiIiIioh9tvvx033nhjootBYDCFiIiIiGhaUVUVfr8fwMh6pgiCYPROISIiIiKisfEf//Ef2Lx586jsq6ioCE888USfy8vKyiAIQr9TWVkZwuEwfvazn6GkpAQ2mw1paWlYuXIlfvvb3xr7ampqwje/+U0UFhbCYrHA7Xbj6quvxo4dO0blsSSCnOgCEBERERHR+PH5fNA0DZIkwWw2j2hfBQUF2LVr1yiVjIiIiIiIekpJSRm3Y61ZswZ1dXXG3//wD/+Ajo4OPPvss8a8tLQ0/OhHP8JvfvMb/OpXv8Ly5cvR0dGBvXv3orW11Vhv06ZNCAaDeO655zBr1iw0NDRg69atOHfu3Lg9ntHGYAoRERER0TQSnS9FEAQAQIfSgZTMFFhl65ACLBaLBTk5OUhPT4csyyMOzhARERERUazbb78dbW1tePXVV1FUVIT7778f999/v7F8yZIluPHGG/HII49A0zT86Ec/wjPPPIOGhgakp6fjC1/4Ap588kmsW7cOZ8+exQMPPIAHHngAAKBpWsyxzGYz3G638bfNZkMgEIiZBwCvvfYavvWtb+GLX/yiMW/x4sXG/ba2Nrz//vsoKyvD2rVrAQAzZszARRddNGrPSyIwmEJERERENI3Ey5dyqOMQbvrdTUixpsDmGFoelYULF+LZZ5+FxWKB0+kc1bISEREREY0FTdMQDoYTcmzJLBkXNY22P/3pT3j88cfx4osvYsGCBaivr8enn34KAPjzn/+MxYsX4+6778Zdd901ouO43W68++67+Na3voXMzMxeyx0OBxwOB1599VWsWrUKFotlRMebKBhMISIiIiKaRqJ7puhqOmoAAO3+drx7+l1cO/faQe8vKysLZrMZgUAA7e3tSE1NHd0CExERERGNsnAwjL99+28JOfa1T14L2TI2zfKVlZVwu93YsGEDTCYTCgsLjd4gaWlpkCQJycnJvXqaDNW///u/4wtf+ALcbjcWLFiANWvW4IYbbsC110bOI2RZxubNm3HXXXfh6aefxtKlS7F27Vp8+ctfxqJFi0b8OBOFCeiJiIiIiKaReMGU6o5q4/5bJ95CV7Br0PsTBAEulwsA0N7ePjqFJCIiIiKiIfviF78In8+HWbNm4a677sIrr7wCRVFG/TilpaU4dOgQdu3aha997WtobGzE9ddfj69//evGOps2bUJtbS1ee+01XHPNNSgrK8PSpUuxefPmUS/PeGHPFCIiIiKiaaTnMF/ekBfVddV4/tbnAQBf/f1X8dfjf8WXFnxpUPtrbGzEqlWrAADvvvsuioqKRr/QRERERESjSDJLuPbJwffGHu1jD5coir3ynIRCIeN+QUEBKioq8M4772DLli341re+hZ///OfYtm0bTCbTsI/bV1lWrFiBFStW4P7778fzzz+PW265BT/4wQ8wc+ZMAIDVasWVV16JK6+8Ev/8z/+Mr3/963j44Ydx++23j2pZxgt7phARERERTSM9e6boQ3xFKztThmZv85D3zZ4pRERERDQZCIIA2SInZBpJvpTMzEzU1dUZf3d0dOD06dMx69hsNlx//fV48sknUVZWhp07d+LgwYMAIgnmw+GxyRVTWloKAPB4PP2u09/yiY49U4iIiIiIpglVVeH3+wF090yJHuILAIrTi1EdrsZfjv4Fdy69c0j77+rqQjgchiQN/2o7IiIiIiKK7/LLL8fmzZtx/fXXw+Vy4Yc//GHMb+/NmzcjHA5j5cqVsNvteP7552Gz2TBjxgwAQFFREbZv344vf/nLsFgsyMjIGFY5vvCFL+Diiy/GmjVr4Ha7cfr0aTz44IMoLi5GSUkJzp07hy9+8Yv42te+hkWLFiE5ORl79+7FY489hhtuuGFUnotEYM8UGpS9e/di27ZtYzLGHhERERGND32IL1mWYTabAfQOpmws3ggA+KjmI1S1Vw1p/5qmobOzcxRKSkREREREPT344INYu3YtPvOZz2Djxo248cYbMXv2bGO5y+XCf/3Xf+Hiiy/GokWL8M477+D1119Heno6AODHP/4xzpw5g9mzZyMzM3PY5bj66qvx+uuv4/rrr0dxcTFuu+02lJSU4O2334Ysy3A4HFi5ciUef/xxXHbZZVi4cCH++Z//GXfddRd+9atfjfh5SBRB6znI2iTQ0dGBlJQUtLe3w+l0Jro4U57f78eWLVsAAIsXL0ZhYWGCS0REREREw9HU1IRdu3YhOTkZ69atAwA8+v6jKD9TbuRMaWhowOvVr+Ojmo8wP3M+7l91f7/7bGxsRHZ2NgDg97//PS699FLmTSEiIiKiCcPv9+P06dOYOXMmrFZrooszZF/5ylcgSRKef/75RBdlQkjk6zmknimPPvooVqxYgeTkZGRlZeHGG29ERUVFzDp+vx/33HMP0tPT4XA4sGnTJjQ0NMSsU1lZiY0bN8JutyMrKwvf+9732ONhAjt37pxxv6pqaFcnEhEREdHE0TNfiqqpqO2s7bXeDSU3QBIlHGk6giNNR4Z0DOZNISIiIiIaOUVRUF5ejp07d2LBggWJLg5hiMGUbdu24Z577sGuXbuwZcsWhEIhXHXVVTFJYx544AG8/vrrePnll7Ft2zbU1tbi85//vLE8HA5j48aNCAaD+PDDD/Hcc89h8+bN+OEPfzh6j4pGVXQwpaWlZVInCSIiIiKazvRhvvR8KU2eJgTDQZgkU8x6GfYMrCtaBwD405E/YSid2dva2kalrERERERE09mhQ4ewfPlyLFiwAN/4xjcSXRzCEBPQv/XWWzF/b968GVlZWdi3bx8uu+wytLe343e/+x1eeOEFXH755QCAZ599FvPnz8euXbuwatUqvP322ygvL8c777yD7OxsLFmyBD/5yU/w//1//x8eeeQRY+xmmjhaWloAACaTCaFQCFVVVSgpKUlwqYiIiKa+UCiEiooKzJgxA8nJyYkuDk0BPXum1HTWAAByXbnGmMn67/Hr5l6HHZU7UNVehb21e7Eib0XcfZrNZmRmZkLTNMiyjM7OTiahJyIiIiIaoSVLlhi/32liGFECer0Lf1paGgBg3759CIVC2LBhg7FOSUkJCgsLsXPnTgDAzp07ccEFFxjjKgORhDUdHR04fPhw3OMEAgF0dHTETDQ+AoGAkURUD6BUV1cP6epEIiIiGp7Tp0/j9OnTOHr0aKKLQlNEz2CKnny+OK8YjY2NaGxshMvlAgA4zA5cPedqAMCrR1+FosYfltflchnbpqWlMQk9ERERERFNScMOpqiqivvvvx8XX3wxFi5cCACor6+H2Ww2TsB02dnZqK+vN9aJDqToy/Vl8Tz66KNISUkxpoKCguEWm4ZIH+LL6XSioKAAJpMJPp8Pzc3NCS4ZERHR1Nfa2mrc8kIGGg16MEUf5ksPpuQ78+Ouf8XMK5BiTUGztxnbzmzrd9+CIBjnARzqi4iIiIgmGp5TTQ2JfB2HHUy55557cOjQIbz44oujWZ64HnzwQbS3txsTk6CPHz2Ykp6eDkmSkJeXByDSO4WIiIjGjqZpRoN0IBCA3+9PbIFo0guHwwgEAgB690zpK5hikS24vvh6AMCbx9+EL+Tr9xgpKSkAmISeiIiIiCYOkymSH5BDZk0N+uuov67jaUg5U3T33nsv3njjDWzfvh35+d0nXm63G8FgEG1tbTG9UxoaGuB2u411Pvroo5j9NTQ0GMvisVgssFgswykqjVB0MAUACgoKcObMGdTV1WHhwoUJedMSERFNB16vF8Fg0Pi7ra3N6E1ANBx68nlZliO9jUM+nPNGfuuZfCYIggAg8ts8KyvL2O7iwoux5dQWNHQ14O2Tb+OGkhti9tvY2Gj0ND9w4AAABlOIiIiIaOKQJMkYmhaIXFik//alyUPTNHi9XmNo4kTkaBxSMEXTNNx333145ZVXUFZWhpkzZ8YsX7ZsGUwmE7Zu3YpNmzYBACoqKlBZWYnVq1cDAFavXo1//dd/RWNjo3GStmXLFjidTpSWlo7GY6JREgwGjfGu9bw4KSkpSE5ORmdnJ+rq6lBYWJjIIhIREU1ZPYdJam1tRU5OTmIKQ1NCdL4UQRBQ21kLAEi1pSLJnNTndqIg4nMln8PTe5/GllNbsK5oHVKsKXHXdTqdAICOjg4moSciIiKiCUO/iF8PqNDk5XK5+uyUMdaGFEy555578MILL+Avf/kLkpOTjRwnKSkpsNlsSElJwZ133onvfOc7SEtLg9PpxH333YfVq1dj1apVAICrrroKpaWluOWWW/DYY4+hvr4eDz30EO655x72Pplg9F4pycnJxmsjCALy8/Nx5MgRVFVVMZhCREQ0RvRgitlsNnr+Eo2E3jNF7+FU1REZOjcvOW/AbZe4l2B22mycbDmJ14+9jq8u+mrc9axWq/Ge7ezs7JVLkYiIiIgoEQRBQE5ODrKyshAKhRJdHBomk8mU0Au2hhRM+fWvfw0AWLduXcz8Z599FrfffjsA4PHHH4coiti0aRMCgQCuvvpq/Od//qexriRJeOONN/DNb34Tq1evRlJSEm677Tb8+Mc/HtkjoVHX0tICoHuIL11+fj6OHj2KlpYWeDweJCX1fSUjERERDY8ePCksLMSJEyfQ3t4OTdPYHZ2GLbpnCgDUdNQA6DtfSjRBEPD5+Z/Hz3f8HDsqd2DDrA1wO3pfDaYnoW9sbER7ezuDKUREREQ0oUiSxN7TNGxDSkCvaVrcSQ+kAJGr0Z566imjof3Pf/5zr243M2bMwF//+ld4vV40NTXhF7/4BWR5WOlbaAz1zJeis1qtyMzMBABUVVWNe7mIiIimOlVVjZwT+fn5kGUZiqIYw28SDUfPYMpAyed7mpM2B4vdi6FqKl49+mqf6+lJ6NmbioiIiIiIppIhBVNo+giFQujo6ADQO5gCRBLRA0B1dTU0TRvXshEREU11nZ2dCIfDMJlMcDgcbJymURE9zJemaajpHHzPFN3nSj4HQRDwSd0nONlyMu46+vuVSeiJiIiIiGgqYTCF4mppaYGmaXA4HHFz2WRnZ8NkMsHn8xk9WIiIiGh06EETl8sFQRCQmpoaM59oOKJ7pjR7mxFQApBFGdmO7EHvIyc5B2sK1gAA/nzkz3EvqtGDKXoSeiIiIiIioqmAwRSKq68hvnSSJCEvL5KslEN9ERERja7W1lYAMPJN6LcMptBwhcNhBAIBAJFgij7EV25yLkRBhCzLSElJQUpKyoDD73523mdhkkw40XICBxoO9NrWZrPBbDZD0zQOTUdERERERFMGgykUlx5MSUtL63Od/PzIkBB1dXVQFGVcykVERDQd6EETvUeKHkzhlf40XHqvFFmWIctyr3wpaWlpaGtrQ1tbW7+//wDAZXXhiplXAABeOfoKXKmumG0FQeBQX0RERERENOUwmEK9KIpinPj21TMFiDTsOBwOhMNh1NbWjlfxiIiIpjRFUdDV1QWgO4hitVphtVqhaRobp2lY9HwpdrsdgiAY+VLynHnD2t/Vc65GkjkJdZ112Fm1s9dy9qYiIiIiIqKphsEU6kXPl2K322Gz2fpcTxAEIxE9h/oiIiIaHW1tbdA0DTabzchbJggCG6dpRKLzpQDo1TNlqOwmO66bex0A4LWK1xAMB2OWs2cKERERERFNNQymUC8D5UuJlp+fD0EQ0NLSAo/HM9ZFIyIimvKik89H0//W86kQDUV0MMWv+NHkaQIA5CVHeqY0NjZCFEWIoojGxsZB7XNd0Tqk29NRW18Lq8kas60eTOns7ISqqqP9cIiIiIiIiMYdgynUy1CCKVarFZmZmQDYO4WIiGg09MyXotP/Zs8UGg59mC+bzYbazsjwrCnWFCRbko11NE2DpmmD3qcsyrhh3g1xt9WT0Kuqio6OjtF4CERERERERAnFYArFUBTFaKQZTDAFgDHUV3V19ZBOwImIiKi3vnqm6Ff6e71eBAKBcS4VTXbRPVNGOsRXtIvyLkJucm6v+UxCT0REREREUw2DKRSjtbXVGKddH1N7INnZ2TCZTPD5fEavFiIiIho6v98Pn88X0xCtM5lMcDgcANg7hYYuOgF9TUck+fxoBFMEQcDGuRuNv7sCXcZ9PSDIYAoREREREU0FDKZQjKEM8aWTJAm5uZErEjnUFxER0fDpQZLk5GTIstxrOZPQ03AoimL0ZrLZbKjqiPxe0/OljNS8jHnG/YpzFcZ9PSDI9ysREREREU0FDKZQjOEEU4Duob7q6uqgKMqol4uIiGg66GuILx3zptBw6L1STCYTZFke1Z4pQKR3ii5eMIVJ6ImIiIiIaCpgMIUM4XB4yPlSdC6XCw6HA+FwGLW1tWNQOiIioqmvtbUVQN/BlOieKcxTRoMVnS+lxdcCv+KHJEpwO9yjfqxj544Z700moSciIiIioqmEwRQytLa2QlVVWK3WQedL0QmCEJOInoiIiIZG0zTjoga9B0pPTqcToigiGAwaDeREA9F7pthsNiP5fI4jB5IoGevIsgy73Q673R53iLn+6NvKFhmekMc4BpPQExERERHRVMJgChmih/iKHq5hsPLz8yEIAs6dOwePxzPaxSMiIprSPB4PFEWBJElITk6Ou44oisxDQUMW3TNFD3T0HOIrLS0NHo8HHo8HaWlpQ9q/vu3j2x+H1WlFeVO5sYzBFCIiIiIimioYTCHDcPOl6KxWKzIzMwGwdwoREdFQ6cGRlJSUfi9qYBJ6GqrBBFNGw4LMBQAQE0zh+5WIiIiIiKYKBlMIQCRfij5O+3CDKUCkdwoAVFVVcSx3IiKiIRgoX4pOX66vTzSQeMN8jUUwpTSzFABwouUEAkoAAJPQExERERHR1MFgCgGIXC2oqiosFguSkpKGvR+32w2TyQSfz2f0dCEiIqKBDZQvRacvb29vZ+M0DYreM0UyS2jyNgEA8px5Mes0NzdDlmXIsozm5uYh7V/fNs+VB1vQBkVVcOzcMQCRAI7JZGISeiIiIiIimvQYTCEAI8+XopMkCbm5uQAivVOIiIhoYOFw2GhoHqhnit1uNxqnOzs7x6F0NJkpioJgMAgA6FA7oGkanBYnnBZnzHqqqiIcDiMcDg85SBe97dy0uQC6h/oSBMF4TzNvChERERERTWYMphCAkedLiVZQUAAAqKurg6IoI94fERHRVNfR0WH0ELXZbP2uG904zaG+aCB6rxSz2YxaTy2A3r1SRlNJRgkA4HDTYWMek9ATEREREdFUwGAKQVXVUcmXonO5XHA4HAiHw6itrR3x/oiIiKY6fYgvl8s1qB6iTOpNgxWdL6WmswbA2ORL0c1Omw1RENHQ1YBz3sjFOgymEBERERHRVMBgCqGtrQ3hcBhmsxkOh2PE+xMEweidUl1dPeL9ERERTXXRwZTB0POmMJhCA9F7ptjt9jFNPq+zm+yYmToTQPdQX/r7Wu+BRURERERENBkxmEKjli8lWn5+PgRBwLlz5+DxeEZln0RERFOV3kN0sMEUfb2uri4OqUn90oMpNpsNNR1j3zMFAEozSwF0B1Oik9Azzw8REREREU1WDKYQWlpaAIzOEF86q9WKzMxMAOydQkRE1J9gMGhceDDYYIrFYoHdboemaeydQv3Sh/kKS2F4Q16Iggi3wz2mx9SDKUeaj0DV1Jg8P3y/EhERERHRZMVgyjSnadqYBFOASO8UAKiqqoKmaaO6byIioqlCzyORlJQEs9k86O3YOE2DofdMaVcj7zO3ww1ZlHutJ4oiLBYLLBYLRHFopwg9ty1yFcFussMX8uFM2xkAzJtCRERERESTH4Mp01x7ezsURYHJZEJycvKo7tvtdsNkMsHn8xlDiREREVGsoeZL0TGYMjaCwSAOHTo0ZYaj0oMpLaHIxTMFKQVx18vIyIDf74ff70dGRsaQjtFzW1EQMT9zPgDgcONhAAymEBERERHR5MdgyjTX3NwMYHTzpegkSUJubi6ASO8UIiIi6m2o+VJ0+vr69jQ6KioqcPr0aRw4cCDRRRmxUCiEUCgEAGgINAAA8pLzxuXYPfOmMAk9ERERERFNdkMOpmzfvh3XX389cnNzIQgCXn311Zjlmqbhhz/8IXJycmCz2bBhwwYcP348Zp2WlhbcfPPNcDqdcLlcuPPOO9HV1TWiB0LDE518fiwUFESufqyrq2OCXCIioh6ic56kpqYOaduUlBQIgmD0CKCRC4VCRq63lpaWSf/7VM+XYjabUeMZn+TzugWZCwAAp9tOwxvyMgk9ERERERFNekMOpng8HixevBhPPfVU3OWPPfYYnnzySTz99NPYvXs3kpKScPXVV8ec5N988804fPgwtmzZgjfeeAPbt2/H3XffPfxHQcMSnS8lLS1tTI7hcrngcDgQDodRW1s7JscgIiKarHw+HwKBAARBgNPpHNK2siwbQ3RyqK/RUV1dHXPxR2VlZQJLM3L6EF8WqwUNnkjPlL6CKc3NzbBarbBarUbP5cGKt22qLRU5yTnQNA1Hm49CEARjqC++X4mIiIiIaDIacjDl2muvxb/8y7/gc5/7XK9lmqbhiSeewEMPPYQbbrgBixYtwu9//3vU1tYaPViOHDmCt956C7/97W+xcuVKXHLJJfjlL3+JF198kY3t40zPlyLLsnFyO9oEQTB6p5w4cYLDOhAREUXRG5WdTickSRry9hzqa/RomoYzZ84AALKysgBEgiuT+beLHkwJikFomgaH2QGnJX7QTlVVBAIBBAKBIT/mvrbVh/rS86bo71fmTSEiIiIiosloVHOmnD59GvX19diwYYMxLyUlBStXrsTOnTsBADt37oTL5cLy5cuNdTZs2ABRFLF79+64+w0EAujo6IiZaOT0XiljkS8lWlFRESwWCzweD06dOjVmxyEiIppshpt8Xsck9KOnubkZXV1dkGUZS5cuhdVqRSAQQENDQ6KLNmz6MF9dWmS4snxn/pj+5uspOm+KpmlMQk9ERERERJPaqAZT6uvrAQDZ2dkx87Ozs41l9fX1xtV+OlmWkZaWZqzT06OPPoqUlBRj0ns60MiMdb4UnSzLKC2NnEwfP36c47oTERGdN1C+lM5AJ144+AKqO6rjLte3a2trg6ZpY1LG6ULvlVJQUACTyWT83pzMQ33pPVPaw5HgxXjlS9EVpxdDFmW0+FrQ4GkwgilMQk9ERERERJPRqAZTxsqDDz6I9vZ2Y6qqqkp0kSY9TdPGLZgCAHl5eUhNTYWiKDhy5MiYH4+IiGiii04+31fPlFePvoptZ7bhuf3PxQ2WJCcnQ5IkKIoy6ZOlJ5LX6zV6oBQVFQGAEUxpamoyghKTjd4zpSnUBADIc+aN6/HNkhlz0+cCiAz1ZbfbmYSeiIiIiIgmrVENprjdbgDoNRxCQ0ODscztdqOxsTFmuaIoaGlpMdbpyWKxwOl0xkw0Mp2dnQiFQmOaLyWaIAhYuHAhgMj44xzbnYiIprvOzk6Ew2HIsgyHw9F7eaATu6p3AQAq2ytx7NyxXusIgsChvkbB2bNnoWkaMjMzjdciKSkJGRkZ0DRt0l7IoweBGgL9J58fS9FDfUUnoedQX0RERERENNmMajBl5syZcLvd2Lp1qzGvo6MDu3fvxurVqwEAq1evRltbG/bt22es8+6770JVVaxcuXI0i0P90HulpKWljdvY2S6XC4WFhQCAQ4cOcTgSIiKa1qJ7pcT7Ln6/8n0oqmL8veXUlrj7YTBlZMLhsDGUl94rRaf/bqmqqpp0v1tCoRBCoRCC4SC61C6IgogcR864l0MPphw7dwyKqhjBFL5fiYiIiIhoshlyMKWrqwv79+/H/v37AUSSzu/fvx+VlZUQBAH3338//uVf/gWvvfYaDh48iFtvvRW5ubm48cYbAQDz58/HNddcg7vuugsfffQRduzYgXvvvRdf/vKXkZubO5qPjfoRHUwZTyUlJZBlGW1tbZP2Kk8iIqLR0N8QX4qq4L3T7wEArp93PQRBwMGGg6jrrOu1LoMpI1NbW4tgMAibzdYr75/b7YbZbIbP50NTU1OCSjg8eq+UAAKACGQ7smGSTH2uL4oiJEmCJEkQxaGdIvS3bV5yHpwWJ4LhIE60nDDer+yZQkREREREk82Qgyl79+7FhRdeiAsvvBAA8J3vfAcXXnghfvjDHwIA/umf/gn33Xcf7r77bqxYsQJdXV146623YLVajX384Q9/QElJCa644gpcd911uOSSS/B//+//HaWHRAMZ73wp0SwWC+bNmwcAOHr0KEKh0Lgen4iIaKLQh7yMF0zZW7sXHYEOuKwuXDPnGizOXgwAeOfUO73WjW6cDofDY1beqUjTNCPxfFFRUa8eQpIkIT8/MjTWZEtEr+dL8SFyO9AQXxkZGVAUBYqiICMjY0jH6m9bQRBihvpiEnoiIiIiIpqshhxMWbduHTRN6zVt3rwZQOSE6cc//jHq6+vh9/vxzjvvoLi4OGYfaWlpeOGFF9DZ2Yn29nY888wzcccKp7HR1dWFYDAISZL6THg7loqKipCcnIxAIIBjx3qP/05ERDTVKYpiJOBOTU2NWaZpmhE0WVe0DrIo46rZVwEAdlXvQkegI2Z9m80Gi8UCTdPQ0RG7jPrX1taGtrY2iKJoJJzvSR/qq76+HoFAYDyLNyJ6z5QurQtAYvKl6BZkLQDAJPRERERERDS5jWrOFJoc9F4pqampQx7GYTSIoogFCyIn1adPn+aJNBERTTvt7e3QNA1WqzWm9y4AnGg5gar2KpgkEy6bcRkAYFbqLMxMnQlFVVB2pixmfSahHz69V0pubi4sFkvcdZKTk5GamjrpEtHrwZQ2pQ1AYoMp8zPmAwCqO6rRGexkEnoiIiIiIpqUGEyZhhI1xFe0zMxMuN1uaJqGw4cPT7qkrkRERCPRX76Urae3AgBW5a9CkjkJQCRgovdOKTtThmA4GLMNgylDFwgEUFtbCwCYOXNmv+vqvVMqKysnzW8Wn88HVVNxTon87stLzut3/ZaWFiQlJSEpKQktLS1DOtZA2yZbklGYEnkOo4f64vuViIiIiIgmEwZTpplE5kvpacGCBRBFEU1NTWhoaEhoWYiIiMZTX8GUZm8z9tfvBwBcMfOKmGVL3EuQYc+AJ+jBzqqdMcv0/eh5WGhglZWVUFUVLpdrwGFPc3NzIcsyPB7PkAMNieL1euEL+RCWwkgyJ8FldfW7vqIo8Hq98Hq9UBRlSMcazLb6UF/lTeVMQk9ERERERJMSgynTjMfjQSAQgCiKCcmXEs1ut2P27NkAgMOHDzNpLhERTRt6MKVnvpR3T78LTdOwIGsBcpJzYpaJgogNszYAiCSiV7Xu5N36d7rH40EwGNtrhXrTNA1nz54FMHCvFACQZRl5eZGeHZMhEb2mafB6vfCEPFBlFXnJeRAEIaFlik5C73Q6ATAJPdF4aWtrw8cff4zq6upJ07uOiIiIaCJiMGWaic6XIklSgksDzJkzB1arFV6vF6dOnUp0cYiIiMZcIBCA1+uFIAjGcEcA4Ff82FG5A0DvXim6NQVrYDfZ0ehpxIGGA8Z8s9mMpKTIkGAcOmlg9fX18Pl8MJvNyMnJGXgDdA/1VVdXh1AoNJbFG7FQKARFUdAV7IIqqwnNl6KblToLFtmCzkAnzinnmISeaByoqoqjR4/igw8+QE1NDT755BPs2LGD3xNEREREw8RgyjQzUYb40smyjNLSyJWKx48fh8/nS3CJiIiIxpbeiOVwOGAymYz5Oyp3wK/4kZOcY1zF35NFtmBt0VoAwJaTW2KW6b1c2Eg2MD3x/IwZMwZ9cUlKSgqcTifC4TCqq6vHsHQjp/+e8qpeQEhs8nmdLMqYlz4PAHCk+QiT0BONsfb2drz//vs4fvw4NE1DZmYmZFlGa2sr3n//fXzyySfw+/2JLiYRERHRpMJgyjQykfKlRMvNzUVaWhrC4TCOHDmS6OIQERGNqXj5UlRNxbun3wUQ6ZXS35BM64rWQRIlnGg5gdOtp435TEI/OJ2dnWhuboYgCJgxY8agtxMEYdIkovd6vQCADrUDwMQIpgDdeVMONx5mMIVojKiqimPHjuH9999HR0cHLBYLli9fjlWrVmH9+vUoKCgAAFRXV+O9997D8ePHOdwyERER0SAxmDKNeL1e+P1+iKLYa4z2RBIEAQsXLoQgCKipqTECPkRERFNRvGDKp/WfotnbjCRzElbmr+x3e5fVhZV5kXW2nOrunRIdTJnIDf2JpvdKyc7Ohs1mG9K2eXl5EEURHR0dEzoI4PP5EAwH4Rf8EAShV/6dRNF7XJ1sPQm7ww6AwT+i0dTZ2YkPPvgAFRUV0DQNOTk5WLt2rTGcodVqxZIlS3DppZciNTUViqLg6NGjKCsrQ11dHb87iIiIiAbAYMo0ogcpXC7XhMiXEi0lJcW42vPQoUP8IU9ERFOSpmlobW0FEBtM2Xp6KwDgshmXwSyZB9yPnoj+47qP0extBgA4nU6IoohAIMBhM/ugKIoxRNdgEs/3FJ1jZSInoo9OPp+VlDWo9xQQucBluInqB7Ntpj0TGfYMhNUwmtQmAExCTzQaNE3D8ePHsX37drS3t8NsNmPp0qVYtmwZLBZLr/VdLhcuvvhiLF261MhfuXfvXuzcuRMdHR0JeAREREREkwODKdPIRBziK1pJSQlMJhM6OjomdAMFERHRcHm9XoRCIYiiCKfTCQA423YWx88dhyiIWFe0blD7yXPmYUHWAmiahq2nIoEYSZKMffJq//iqq6uhKAocDsewfw/pQ4PV1NRAUZTRLN6o8Xq98AQ9UCUVBc6CQW2TlZUFVVWhqiqysrKGdLzBbisIgtE75UTnCciyDFVV0dXVNaTjEVG3rq4u7NixA0ePHoWqqsjOzsa6deuQl5fXb4BTEATk5eVh/fr1KC4uhiRJOHfuHLZv344DBw4gEAiM46MgIiIimhwYTJlGJnowxWw2Y968SGLSo0ePIhgMJrhEREREo0sPcqSkpEAUIz/D9F4pK/JWwGV1DXpfV866EgCwo2oHvKFIjgzmTembpmk4fTqSY6aoqGjYPTDS0tKQlJQERVFQV1c3mkUcNdE9U/KceYkuTgw9mFLeVG7kTeH7lWjoNE3DqVOnsH37drS2tsJkMmHJkiVYsWJF3N4ofZFlGfPmzcO6deuQm5sLTdNw9uxZvPfeezh16hR7jhERERFFYTBlmvB6vfD5fBAEYULlS+mpqKgIycnJCAaDqKioSHRxiIiIRlXPIb7a/G3YW7sXQCTx/FCUZJQg35mPgBLA9rPbY/fLxulezp07h66uLsiybCRgHo6eiegnGk3T4PP5Ij1TZHXCJJ/XlWSUQBRENHoaIdoipyITOf8M0UTk8Xjw4Ycf4vDhwwiHw8jKysLatWtRUFAw7ECx3W7HsmXLsGbNGqSkpCAUCuHw4cPYtm0bGhsbR/kREBEREU1ODKZME7W1tQAijSyyLCe4NH3Tk9EDwNmzZzlmLxERTSk9k89vO7MNYTWMOWlzMMM1Y0j7EgQBV86O9E559/S7UFTFuGCCSeh703ul5Ofnj/i3kN5g2dLSgs7OztEo3qgJhUIIKSF4Q94hBVNaWlrgcrngcrnQ0tIypGMOZVubyYZZqbMAAM1qJN8Pg39Eg6P3sNu2bRtaWlogyzIWLVqEiy66CDabbVSOkZ6ejksvvRSLFy+GxWJBV1cXdu/ejT179kzYoQ2JiIiIxguDKdNAfX09jh49CiDSgDDRZWRkGF3MmYyeiIimClVVjSvwU1NTEQqHsO3sNgDAFbOG1itFtzx3OVxWF9r97dhTswdJSUmQZRnhcHjCNfInks/nQ0NDA4BIL9iRslgsyM7OBjDxeqd4vd5IIEVUYTPbkGodXI9kRVHQ3t6O9vb2ITeYDnVbfaivSn/kuWMSeqKBeb1e7Nq1C4cOHUI4HEZGRgbWrl2LGTNmDLs3Sl/0Hnjr16/H7NmzIYoi6uvrsWfPHoTD4VE9FhEREdFkwmDKFNfa2oqPP/4YmqahsLDQSJo60c2fP99IgjhRxyMnIiIaCr3B2GQywW63Y1f1LniCHqTb07HEvWRY+5RFGZfPvBwAsOXUFgDdvV70IcUo0ttV0zRkZGQgOTl5VPapD/VVXV09oQIBPYf4Gu1G1tGwIGsBAKCivYJJ6BOkvr4eZWVl+Pjjj9HQ0DCh3sMUS89hsm3bNjQ3N0OSJFxwwQVYtWoV7Hb7mB7bZDKhtLQUa9asgSzLaG5uxr59+/h+ISIiommLwZQpzOPxGFcPZWVl4YILLpiQJ9Tx2O12zJkzBwBQXl7OK6CIiGjS6znEl554/vKZl0MUhv+T7NIZl8IiW1DTUYMjzUdihvoiIBwO4+zZswBGp1eKLisrC1arFcFgEPX19aO235GKTj4/0fKl6ApTCpFkToI/7IdiivRk4ft1fGiahhMnTmDv3r3o7OxETU0NPvroI2zZsgUHDhzAuXPn2Ct8AvH5fPjoo49w4MABKIqCtLQ0rF27FkVFReN6XpeamoqLLroIkiShoaHBuFiPiIiIaLphMGWKCgQC2L17NwKBAFJSUrBs2TKI4uR6uWfPng273Q6fz4cTJ04kujhEREQjojcWp6am4kjzEdR11sEiW3BxwcUj2q/dZMclhZcAAN4++TaT0PdQV1eHYDAIm80Gt9s9avudqInovV6v0TMlLzkv0cWJSxREzM+YDwBoQSTHCpPQjz1VVfHpp5/iyJEjRq/1WbNmwWKxIBgM4uzZs/jwww+xdetWlJeXo729nQ3mCaJpGqqqqozk76IoYsGCBVizZg2SkpISUqb09HSsWLECoiiirq4O+/fv5/uDiIiIpp3J1bpOgxIOh7Fnzx54PB7Y7XasXLlyQied74skSSgtjYypfeLECdTW1ia4RERERMMX3TNl66lIr5SLCy6GzTTypMFXzLwCgiDgSNMReCUvAKCzs5PJgtGdeH4oeQW8Ie+gGgn1RPRNTU3wer0jKudo6TnM10Sl502pCdYAYDBlrAUCAezcuRNVVVUQBAEXXHABFi9ejAULFuDKK6/E6tWrUVhYCJPJBJ/Ph5MnT2L79u0oKyvDsWPH4PF4Ev0Qpo1AIIA9e/Zg//79CIVCcLlcWLt2LWbNmpXwUQYyMzOxbNkyCIKA6upqHDx4kAEVIiIimlYYTJliNE3Dxx9/jNbWVphMJqxcuRIWiyXRxRo2t9uNnJwcqKqKffv2oaKigj/YiYho0gmFQkZOCL/kx6HGQxAEwch3MlLp9nQsy1kGANhesx02mw2apk37Buq2tja0tbVBFEWjF0l/uoJdePaTZ/HAWw/giV1PoMXX0u/6drsdGRkZACZO75TWjlYE1SA0WUNucm6ii9MnPZhSHaiGoiro6Ojgb7wx0tnZiQ8++AAtLS3G+UH0kHeCICAjIwOLFy/GlVdeieXLlyM3NxeSJKGrqwsVFRV499138f777+PUqVPw+/2JezBTmKZpqKmpQVlZGRoaGiCKIubPn49LLrkEDodjxPvuDHRC1Uae68TtdmPp0qUQBAFnz55FeXk5P7tEREQ0bUy+7grUJ03TcPjwYdTX10MURVx00UUj/uGdaIIgYNmyZSgvL8epU6dw7NgxdHZ2YsmSJZOytw0REU1P+nA5drsdH9R+AABYlL0ImUmZo3aMK2dfib21e/FRzUeYlTYLPp8PbW1tSE9PH7VjTDZnzpwBAOTm5vZ7cYmmadhdsxsvHX4JnmDkCvyjzUfxo7If4SsXfAUr81b2eUV4YWEhmpqaUFVVhXnz5iX0ynFN09Dc0QwASE1OhUUe3AU1mqrh9NbTkCBBgwZfiw/IGsuSAqm2VOQk56Cuow4dwQ6kiWno7OyE0+kc2wNPM42Njdi3bx8URUFSUtKA5weSJCEnJwc5OTlQFAX19fWoqalBU1OTEZwsLy9Heno6CgoKkJeXl/DeElNBIBDAwYMHUVdXBwBISUnBhRdeiOTk5BHvu7azFi8dfglHmo7AZXXh0hmX4tLCS5FiTRn2PnNzcxEOh7F//36cOnUKkiShpKRkxGUlIiIimujYGj2FnDp1CqdPn4YgCLjwwguRlpaW6CKNCkEQsGDBAjidThw4cAB1dXXwer1YsWIFbLaRD41CREQ01vQhvmwOG3ZW7wQAbJi1YVSPUeQqwtz0uTh+7jjO+M8gFanTOm9KIBBATU1kCKn+Es83e5vxhwN/QHlTOQAgz5mHzxR/Bm+ffBunW0/j2U+exf76/bj5gpuRbOndsOl2u2E2m+H3+9HY2Ijs7OwxeTyDEQwG0enrBADkpQ0uX4qvxYdPnvkE546fwyt3vwIAOPTLQ/Bu8GLutXMhWwc+XcjKyhrWlekLMhegrrMOHehAGtLQ3t7OYMoo0TQNp0+fNnoNpKenY/ny5TCbzYPehyzLyM/PR35+PgKBAOrq6lBTU4OWlhY0NzejubkZJ0+eRGlpKTIzRy8wPFGpqoqqqir4fD6kpqYiLS0NJpNpxPutq6vDwYMHEQgEIAgCiouLMWfOnBHnu+wKduH1itex7ew24/PZ5m/D6xWv481jb2JpzlKsK1qHOWlzhhUQKygoQDgcxsGDB3H8+HFIkoS5c+eOqMxEREREEx2DKVNETU0NyssjjQClpaXIzZ24wzoMV0FBAZKSkrB37160t7fj/fffx/Lly6dM0IiIiKYuPahRGahEKBxCQUoB5qaNfqPTlbOuxPFzx/Fp26e4TLpsWgdTqqqqoKoqXC4XXC5Xr+WqpuKdU+/gtYrXEAqHIIsyPlP8GVw1+ypIooQl7iV468RbeL3idXxS9wlOtJzALYtuwWL34pj9iKKI/Px8nDp1CpWVlQkNpvh8PnhCHqiSigJXwYDr1+6txYHnDyDkC0G2yCj+TDEaDzWiuaIZJ946gaodVZh3wzwUXlwIQRz93gelmaV459Q7qFPqUGQpwokTJyDLMtxuN3s7jICqqjh06BDOnj0LINJ76oILLhhR47zFYkFRURGKiorg9XpRU1ODkydPoqOjA7t27UJWVhZKS0tHpSfFRKNpGqqrq1FRUQGfz2fMFwQBTqcT6enpSE9PR1pa2pCCVcFgEIcOHTKCvk6nE0uWLEFKyvB7jABAWA2j7EwZ3jj2BryhSC6nC3MuxA3zbkBleyW2nd2Gky0nsbd2L/bW7kWeMw/ritZhZd7KQfdm0xUVFSEcDqO8vBxHjx6FJEmYNWvWiMpPRERENJEJ2iQc4LSjowMpKSm8eu28c+fOYdeuXVBVFbNmzcKCBQsSXaQx5fV6sWfPHnR0dEAURSxatAgFBQM3GBARESXKli1b4PP78KH2IVqFVtxx4R1Ylb9q1I+jaRoeLnsYDR0NWOpdilxHLq666qpJnT9tODRNw9atW+Hz+bBkyZJevxMq2yvx/z79f6hsj+Q5KU4vxlcXfRXZjt6BkKr2KjzzyTOo7awFAKwpWIMvLfgSbKbu3rGdnZ0oKyuDIAjYsGEDrFbrGD66vtXW1uK3b/4W7WjHrdff2ivwo1P8Cg7+90FU76oGAKTOTMWFX7sQSVlJ0DQNDQcaUP7HcngaI0OeOfOcKP1iKTLnj27vg2A4iAfeegBqQMUGeQPMQqQhOjk5GcXFxcjJyWFQZYiCwSD27t2Lc+fOQRAElJaWYubMmWPyPAaDQRw/fhynT5+GpmkQBAGFhYWYN2/elKhzNE1DfX09Kioq0NkZ6fFltVqRkZGB1tZWeDyeXtv0DK709Tw0NDTgwIED8Pv9EAQBc+bMQXFx8Yh7oxxqPISXDr+Ehq4GAEC+Mx83LbwJxenFMetVtVeh7EwZdtfsRigcijw22YrVBauxrmgd3A73kI577NgxVFRUAAAWLVqEGTNmjOhxEBEREU1UDKZMcp2dndixYwdCoRBycnKwbNmyaXHSqSgK9u/fb4wrPHv2bMyfP39aPHYiIppc/H4/tmzZgiZvE3Zbd8Npc+LRDY9CFsemg/D2s9vxhwN/QFZzFi7KuAgXXXRRQntLJEJ9fT327NkDs9mMDRs2QJIkAJHG+9crXsc7p96Bqqmwm+z4QukXsKZgTb+/IRRVwWsVr+Htk29D0zSk2dJw+5LbMS9jnrHOjh070NLSgpKSkoQNdXPs+DE88/Yz8Nv9eOimh5Bhz+i1TuupVnz8u4/hbfZCEATMvW4u5m6ci47ODhQXRxpcjx07BqfDiTPbzuDYG8cQ8kYaW7MXZaP0C6VwZMfm3Ghra4vZNl5PoL48sesJHGk6gk3FmzADM3D69GmEQpHjORwOzJ07l3k5BqmrqwsfffQRPB4PZFnG0qVLx+Wz7/F4cOTIEeN3uSzLmDt3LmbOnGl89iab5uZmHDlyxOjdZzabMWfOHBQVFRmPye/349y5czh37hxaWlqMgEu05ORkpKWlGQEWSZJw+PBhVFVVAYi8xy+88MIhfWbiqeusw8vlL+Nw4+HIcS3JuGHeDbi48GKIQt8BGm/Ii51VO1F2pgyNnkZjfklGCdYVrcNi9+J+t9dpmoajR4/ixIkTEAQBS5YsQX5+/ogeExEREdFExGDKJOb3+/HBBx/A5/MhLS0Nq1atmrQnLMOhaRqOHTuGY8eOAYiM171s2TImpiciogmlrq4Oe/fuxf7W/ahOrcZn530WG4s3jtnxQuEQvv/O9xGuC2OxbTFWL149LRIDq6qKtrY2nDt3DlVVVfB4PJgzZw7mz58PADjSdATPH3gezd5Igvbluctx08Kb4LQM/rfkiZYTePaTZ419XDHrCnyu5HMwSSZUVVVh//79SEpKwvr16xPS+L/9o+14Y88bUFNV/PwrP48pg6ZqOPbmMRz/63FoqgZ7uh0Xfu1CpM2JDJcane+loaEBWVmRDPRBTxDH3jiGM2VnoKkaBFFA0boiFH+mGOYkc7/bDsaWk1vwx/I/YkHWAnx75bcRCoVw+vRpnDp1ygiqJCUlGUGVkV65P1U1NTVh3759CIVCsNvtuOiiiwYccksf3m603qvnzp1DeXl5d44omw3z589Hbm7upAmGtbW14ciRI2hujnzGZVnGzJkzMXv27AHzowQCAbS0tBgBlo6Ojl7rSJKEcDgMQRAwa9YszJs3b0Tnb56gB68fex3bzmyDqqmQRAlXzLwC1829Lqb33EA0TcOR5iMoO1OGAw0HjBwrqbZUXDbjMlxSeMmAdaWmaTh8+LCRw3Pp0qVTcuhpIiIimt4YTJmkFEXBhx9+iPb2djgcDlx88cVDGqN3KqmtrcX+/fsRDoeRnJyMFStWICkpKdHFIiIiAgAcOXIEHx/+GHu69iCUGcLPNvwsbiLz0fR6xet4e9/byPJk4coLrsSqVaM/pFiiKYqClpYWY2ptbYWqqsZyURSxfv16qLKKlw+/jF3VuwBEGgf/7oK/w6LsRcM6bkAJ4I/lf8T2s9sBAG6HG3dceAfyHfnYsmULFEXB6tWrkZHRu1fIWHvp7Zew9/hepBal4sGNDxrzvc1efPy7j9F6qhUAkL8yHwu/shAmW3fj8EABka6GLpT/sRwNByLDB5nsJhR/phhFa4vQ3NI87GBKTUcNfrztxzBJJjx+9eMwSZEyKYpiBFWCwSAAwG63Y+7cucjPz2dQJcqZM2dw6NChSK+ptDQsX748ZngpTdPQ5G1CdUc1qtqrIrcdVWj1tcIqW+F2uJHtyEZWUlbkflLk/lDzZ+jHqqmpwdGjR438Ii6XCwsWLJjQeQ47OztRUVFh9K4RRREzZszA3Llzhz1kWTAYNHqt6MEVTdOQlJSEJUuWjOj5CKthbD+7Ha9VvGbkRVnsXowvlH4BWUmD//zFc857Du9Xvo8PKj9AZyDS20YSJSzOXoyFWQtRmlmKVFtq3G01TcOBAwdQWVkJQRCwYsWKadczkoiIiKY2BlMmIVVV8dFHH6GpqQkWiwWXXHIJ7HZ7oouVUG1tbdizZw/8fj/MZjOWLVuWkEYMIiIiTdMQCATg9/sRCARQUVGBXSd3ocpahRWlK3Dr4lsHva9wKAxPgwfWVKvRC2AwOgOd+D9v/h/Y6mxYmrcU61evR1JSEux2+6RthNav+u7ZMBnNYrEgLS0NaWlpyMzMRHl7OV46/BK6gl0QBAHri9bjhpIbYJVHntPkUOMh/P7T36Pd3w5REHHd3OtQ4C9AVWUV8vLycOGFF4751fiqqkaSzns88Hq9eHPXmzh77ixKLyzF7WtujzRs767Bwf8+CMWvQLbKWHTzIuRdlNdrX4PtXdJ0pAnlL5ejoyZy1X1SVhKyL8/GwssXDrhtPJqm4fvvfB9t/jbcv+p+zM+cH7NcURScPXsWJ0+eRCAQABDp8TB37lwUFBRMyvezpmlQFMWoIwRBgCAIEEURoijG3O/5t76uvh+9JwAA5Ofno2RBCeo8dTGBk5rOGgSUwJDL6bK6kO3IRnZSdsxthj1jwKGfwuEwTp06hRMnTkBRFABATk4O5s+fP6EuevL5fKioqEB1dbWR9yU/Px/FxcWjfn4VCoXg8XiQnJw8ot4ohxsP4+Xyl1HXGQn85Dnz8KUFX0JJxuj2QFRUBftq96HsTBlOtZ6KWZabnIvSzFKUZpaiOL3YCIICkfflJ598gpqaGoiiiIsuugiZmaObb4mIiIgoURhMmWQ0TcOnn36KqqoqSJKENWvWjHiM3anC7/djz549aGtrgyAIWLhwIYqKihJdLCKaIjRNQzgchqIoCIVCCIfDkGUZZrMZJpNpzBpNNU1DKBRCIBAwJv0qbb2BLd5tf8sGup0sw7GMN1VVjSCJ3gga7zYYDMY08geUAD6q/Qid7k48dOVDyHP2bsgGIq+1p9GDttNtaD3dirbTbWivaoemRvZlT7cjZUYKUgq7J0ty31dM/79P/x8O7DiANEsa5qbNhUW2QBAE2O12JCUlxUwOhwM2m23cX3tN02ImVVWN+4qioLW11QigdHV19drebrcbiZ7T0tJgt9vRFepCk6cJbxx7A+VN5QAijY23LLoFM1Nnjmr5PUEP/vvQf2NPzR4AQIGlAHO75sJuskOSJFitVlgsFlit1j7vy3Lfwyzpn3+v12sETKInn88X81471HgIrf5WXH3l1bgs/zIc+MMB1O6tBQCkz03HkjuWwJ4ev4F4KEN1aaqGyh2VqPhLBQKdAXR5u/B3z/8dAKD6TDXyZsR/j/fluf3P4cOqD2GWzMiwZyDDnoE0WxrS7enIsGcg3ZYOl8WF5trmmKCK1WrFnDlzUFhYOCGGuu0ZJOlZV0T/HQ6Hh3UMVVOhQoUGDWEtDG/AC0/IAyVNQZO5CU2+pl5BRgAwSSbkJuci35mPAmcB8p35cDvc6Ax2oqGrAY2eRtR31aPB04CGrgZ0BXt/3nSiICIzKROZ9kyk2lLhsrqQao3cuqwupNpSYZMj9YkeUK6srISmaRBFETNnzsTcuXMHHDZrLAUCAZw4cQJnzpwxerS53W6UlJQMODzaWNC0yOsZCocQUkMIhoPGff02oASw/ex2HGo8BABwmB24oeQGXFJ4yaDymoxEVXsV9tfvx+GmwzjTdibmPSaLMuamz8WCzAVYkLUAOY4caJqGffv2ob6+HpIkYeXKlUhPTx/TMhJNFpqmIRgMIhgMGr/rA4EAFEXp9btInwAYv5H0ffRcbrFYjN91drsddrudv+mJiMZAQoMpTz31FH7+85+jvr4eixcvxi9/+UtcdNFFA243nYMpx44dQ0VFBbtN9yEcDuPAgQOorq4GABQVFWHBggUJuWpRP6FWVRXhcLjf276WCYIAWZaNyWQyxf1bbwAlolh646yiKEYgJPp+OBxGKBQyAiT68uj70X/39ZUpCAJMJhPMZnOvqa/5AIwTKL1xra9pvL+qBwrG6PXNUG6j70dfYS1JUq+rsPu6OlsURaPOjDeFw2HjNY+uX/Wp50nnYO8DkRNYPX9DPKqmIqyFEVbPT1oYokkEZKDeV4/T/tOYOXsmHlj9gLFNsCuItjNtaD3VGgmenGkzEn1HM9lMCPniH9uWaosEVs4HWVwzXLA4IwGW+q56/OTNn8DsMUNURMhhGXbJDqtshd1kh022wWaKTCbRBFEU4wZagMj3a/TnZrCT/rz3vI2+P1iCIMCWZIPZYQZsQMgcQqfaiRZfC1p8LWj1t6LV1wpFVYxtZFHGZ4o/g6tmXwVJHFxju6qoUAIKZIsMUR7c74c9NXvwwsEX4A16kXwuGTPkGbDKVpgkE0yiybiVxd553XoGXcxmMwKBgBFA0a/qB4BgOAi/4kcgHIjcKgH4w374RT88mgchKQTFouCuWXeh+U/N8LX6IIgC5n12HuZcPQeC2PdvheHkPVH8Ck68dQKf/uVTfHnzlwEAL3z1BaRnp8OeaUdSVlLslJkEk713A/qxc8fw5O4nEQr3/RkDALNkRro1HU6/E3K7DJNqglW2ItmejIL8AqM+ESBAFM7XIecbmUVBhCAKEBE7X6+f9PojrJ1/r6tKd92ihntN0XVLSAkhGAgaQRIN5+uaAW4FSYAgC93HDCtQVMU4fjjcfaywFqlXetIEDd50LxR79/skxZqCfGd+TOAk25E9pAZ3T9CDRk+jEVyJvh3odQIiwZvoAItDcyBYH4TSpcAiW5BkTUKaKw1m2QxZliFJEiRJirkfPfWcLwhCr7qkZz3fV73j9Xpx5swZ47OVkZGBkpISpKamQtVUdAW70BnoRGewM+a+fhtQAsZrCKDXfQDG6xx9X18nOkASHTQZbH0oCiIun3k5NhZvhN00/qMTeIIeHG0+isNNh1HeVI5WX2vMcpfVhdLMUpSklyBQGUB7SztkWUZBQUHMeUy8yWQyxbzG01X0e2E6Pw9jLfo3oz6N5De3pmnwBXzw+Dzw+Dzw+r3w+r2RQHogMgWCAYSCIYS1cOS3oxo2fkOqmhpTX6ia2v29cX6+qkWCv/GWmSQTLJIFFtli3CYnJSMlOQUZKRlwJbvgcDiMYMtYXISgX3wW7/d4IBRAUAlGpnAw8t2pRG5D4RCUsBKpD6PuK2ElZgqrYeNW/67s9f3a47tWhQpoMJ6/6OWiKMJsMsNqscJqtsJuscNmtsFutcNhc8BhdSDZmoxkWzKSrcmwm+wwib0votM0DcFwEJ6gB52+Tnj8Hnj8HngDkfeAL+iDz+9DIBSAPxh5HwCA1WKFzWJDkjUJSdYkOOwOOK1OOJOccCW5kGxNjun5R9ND9HmofsHHdPpe1OsR/QLWgdpTo2+j65zFixePWRkTFkz5n//5H9x66614+umnsXLlSjzxxBN4+eWXUVFRMeDJmx5Mqa+vh9PpHPANFW95vJP5wd7qT9lgG5F6zuvraoO+Jv3D4/f7cebMGQDAokWLMGPGjH4f93SlaRpOnjyJo0ePQtM0ZGRkoKCgIGadeK9Nz/v633pQpGcj7EB/D/eKw+HQG3J7npDEey9Gb9Pf/oDYSnwotz3n9TV/KPsa6Bg9l+mNIz0bcwfzd7xj9Lzf17K+egaMx/2B6q+B6raRiC7HYKZ4Q5X0rPN6Nor0VT/2bOSN/iyO9KQonuggpx5kGS690V1RFaMhXlGVuPPDQhiqpCIshBESIw1ZAgQImhBzC+38/Dj3AUDQhO5bfd75xkd9n5H/Qu9lQ70/BvuKOWGKOrHs1VjZzzIgaj9x5gGIu43emBkWw1BEBYqoIIggQmIIqqhClSKTJmnQRA2IrmbDwB25dyC7K9vodeJp8vR6T0gmKRIUmelC6sxUuGa6YEuzQfEpaK9qR/vZdrRXRqauhvhXjltTrEZw5ZB2CAeCB9AcbEYY4ci7QRUgaRLEsAhRFSGFJJhUE+xyd6Al+tb47TLAyX2vk9fz8/TnVn9e9RNa4+8eyzXt/HMnAKpJRcAcQJfUhVa0wqtGcgNAA6BGJiEsAOHIcyyokcfnlJ3Is+XhirwrkCKmQPEpCPlCUHwKFH+c+37FWEdVuvOuyFYZ5iQzzI7IZEoy9f7bYYY5yQyv5MXLp17G4dbDgAaIYRFCWIi5lVQJVlhhgQUmzQQZMkyiCWbJDFmUYZIiAZdQOAS/4u8OnGh+eAUvFEmBKqsxU/R7TVAFzDg8AwvOLoAAAUlZSVh651K4ilxx3yvRRpJE/uzRsyiaXwQgEkxx2B19rmt2mGOCK0nZ54MsaSZ0qB1oCbTgnO8cWvwtaPY245z3HM75zqHd3x67Iw0wd5lh6bBADI/sohkBgvFeHA2aGKkDVPl8nXD+viZpRh2hSmpsHdHnzmLv63W3RbLAIlqQZEtCfnI+cm25yLXlwm11wy7aEQ6FEQ6GoYZUhEM9boPhyH1FhSiJEOUekylyK5mkXssEWUCn0onmYDNag61oC7ahPdiOtmBbZAq0wRPsXbfp5Zf9MmytNohK92smizIkQYrcihJk4fzt+b+jl+n3RUGMqXOA3kGLmOXn6x19nqIqUC0qxAwRXrM3EjQJdsIb8kZ6A0bXLyoAJfL5Qvj8fO38fK17ipmn9piH7nkAABHdr794vi7UnxIBkKXz9YEsG/dNkglp9jSsL1yPNFta5DGpWmTSIrcx86Lm68sESYAoicO7lcVe84BI4L68qRyHmw7j2LljMcE2AQLyOvLgUl0wS2bje10PaOr3jQBo1DL93MYkmyLvU2EQvzUhQBAF4zeDHkSVRAmyFAnIxdyKkcCNKIrGbc/7AOL+xgyFIg29wVAQoVAo0jAcCkYuwlEifyuKEvkuVCPfifpkfEfGuY35/tRgPHajnKIEURKNsve8lcXux2d8JtTY/fc6nqoa6+mfEeNCITlyYZDFZDFuzSYzrGYrLGZLv+eg0ecg4XAYQSWIQCiAUDgUuVVCfTaqq5pqBKf1oLJe7rAW7n5OVTXmVtM0hMIhhJW+L/JQw2qvcyMV3dsPXC2PzveF8V1x/ntCE/Uf5uePEf0d0d+888SwCFGJ/K4TwyJ6FlMURFgkC8ySGRbZApvVhmRHMlzJLtistkjAQlGMwL6ixN4az6N+YUF0Q6Z2/nnV4k+j+R2bSJoYuRBCD/6Gw2GElUjdIKhCr+d8NIhypD40mU2wmC1G4MdsMkOWur8jZEmOfF4lE0yyybi1yBaYpMhn2SybYTFZIItyd/0bVR8DgwvgxmtL1YODxkUgcS5CCavd5+Z6Pa3f7/l3dFnirhv5I2659ccWvS/9OyG67o+u63tO+rKebRbhcDhyMWTQHxM08wV8kaBZwAd/0I9gKIhAMBAJoCohqOHIjwERUc911L+Y8/aoxxD9vShLMmRT5HvRZIpMZpMZZlN3PW0xR+poi9nS/T45Xy/r32mDaa8Z7HthMO8V/TmL7pEXCAYiweaAH/5gJNgcDHQ/Z0a9Ev29id7fo3Hrm/Pr/p87/8+Iy9+XhAVTVq5ciRUrVuBXv/oVgEijWUFBAe677z58//vfj1lXvzJX19HRgYKCAnz/M9+H1TTyMa8nHRuARAw1PIJ3yqh/eQ5md0EAnYNct99DDWIH/a0i9HMr9Phb6zFP33d/01iaGr95Yo3hYxpOdSpg5F9Qo8F4n0+213wI5dWE7hOUmBOSnvcFLeY2erkxL3q/mmY0bkU3pOg/qI2gRdRyQRNi9qkJkUaUePeNhtKeb5XReq167EcPtsRd3scxo4M1fe631xPXe71e+9H6mB998ijEmRe1nrEroe9lcen763ncOO+B6PIbDYCCFNMY6FAcKEwq7LWtI9sREzhx5jshSoNrGFb85wMsld1Blq76rvh1kQb4w374Qj74FF/kNuSDL+yDX/Mb77PoH+76Y4p+/0V/jqI/N5FDaL3vR/99fruez2nM+hpi1otpvDzfqGnSTLAI3VddWiUrzLIZVskKi2wxGuwSQgOaA81oU9siV52rQQS1IIJqEIqm9P7e1+9qQkygU0N3o4rx3J2vO2ySDVbJCqt4fpKsxvNhFrofe+HFhVhw0wLIlt69YeIZSTAletuaszVIEpLgafT0mgIdg8/bIYiRRlFRijSGatAQ0ALwq374lMj71qf44A174ff6IxexRAXqhl1H9qxLhNiTTv1kHIidL8qRcgpSbK+Y6MZdURDjNiBLiAQJ4v0TtchyUf+nRaaYRvIJRBAFQET3ex9BBLRAJCio+o0pGAx2Bxn070f9flSdHzM/6n5kYc+Dx/lu6Lns/HxVUqGKaqQM4e76xahnBJMRwNADnnovs+hebjF1ZdSxes6LaczB+cYa/VboMUEc+DtqghAEwQiuiLIITdDQHmrHucA5NPmb0BHqAARAUqTYizj0Oq/na9/zN8h4PY4eDXVAj0akqGAcMLzf/FOZ0eB3/nmMuYAUWmzdPIH1+q04HOd/KxrfB1JUw6wsGoEuWZa7g8RRQeTohtPonpT6d4kIsdd3jP6dFAgH4Av54Ff88Ia8RqNuIBCAElKMC04EVRi/z1ofv9ljvj/170YxqmFfn6/XkVG3kiBFgqTnb2MCAVGf35iG6R4N6/pzFlbD3UE8pbs3TEwQTu1u1B2M6CCu8dpLImRRNgLkshT5bRYKh4zAbEgJ9bqyfkwN4i3Q8/yt50W3wOT4bPenr2BMr3n6xatDeC9MJH29lvGW93keNcB7Jl7bVvR3QczFdYMUUy9Ht8tEn3f20c7z2IuPDfo4Q5WQYEowGITdbscf//hH3Hjjjcb82267DW1tbfjLX/4Ss/4jjzyCH/3oR732892V34VF7nus8EGLd6LfQ88XcCj1Rbwvqrhf1INo2FFFFapp8n1wE0VQBUhBaUSvV8wHU/873oc2alnPRqQxpTfSovcJyZgajcfW87M31sfrqyF4sD8mB1neuA3L8cqgxf/C6bdxua/1NKF3PdW7YP3Wc6Ou5+PoMW9AcRq04h6mR+N2f5/RiSimAT7qZKrnfePK3PON9ED38xPvylt9fswt4qwXZ5ux2pex3TD3pYtpjBS6h+2JvnrHOBntuUw/0Ucf8+Lc1xu2JEEynv94t6LQfwOYOckcEzhJnZkad8ijkVACCjqqO4zeK+1n2+Fr9UELa1AV1bgqKpqqqZFeECE/vIo3JuASDAdj1jWuDI56/vp7PYC+n9dejZBxlpslc/eQFeeHrYg3XJfRmC0JkExS5NYswWQzQbbKkG0yZKsc+dsmx8yPt45kkaD4FQS7ggh5Qgh2BRH0RN0//3fP5dG9WnrSNK17eB89F4J6/spmNdh9qyowi5GrRq2y1QgYWWUrzJJ5wPrM7DBj0c2LkLM0p/8VE0DxK/A0eeBt8qKroQveJq8RaPG3+0f3YEZMZRD1GLrfez3f3xP5+6MnQYy8/0VTpFeJZO6+HzNPFo1bTdWMXipqSI3cRk3GsjjLR3QaqUXqHkVTjKtW9fvxemv2XE/V1N71OBDzuvV8TaPnyaIc6REmyTCL5pjAiUk0xbzuev1i9M6Rzu9TjJ0goNe8mOn8Npp2/juuZ2+SqPnR84x1zs+PezwhTjmE3n9rqgY1rEa+E/q5VZU4y4YRuAuGg2jztaE90B7pTdCjB6kKNe59vdeB8Tx0f6Bj3kORG814/XtdBBFnm7i/rXsG76LXRffvycGcFxrffyJi65Oo70z9t0W85THBW/19onU/HzE9TKB1v0+inlf97/4az3o1lunfy/rsqGMbwdvz78GY56uv57Q/QvfxBroquufvtuiy9lyu3xrDOvZsjD9/X2/kNnr86PP0Rvkh1Ps9n2NRiDT2T7TvDlVTEQxHhtgKhAPwh/zG1fN6rr+YXgr6/eiArxgnyNHj+dT/1p/PeNNEe24GdP5zp4ZVKKqCoBI0euuIohjTe1CWIkOwj8Zj1FTNGBbN6L11PuhjDEXas2eber4uVePUCVr3ORg0jE4ZB9O2oH9s+7nQaaD2lZjloxHz7FHXA8OrxzRo3QG98wE0/Vb/XOi9IY1AVPS5cNS/nn9Hn0cbr2lUL6Do3wgx35f6b4U+LlJJpLgXrerB5/O/WURRNG6NuqZHENkIwEb/jd7rfu35r43ZY0lIMKW2thZ5eXn48MMPsXr1amP+P/3TP2Hbtm3YvXt3zPp99Ux5+423jbG8x1vPL81xP36irrjUjz/Sxy+M4mMY5m7G9DUcxK5H4/En6n04mt3+RstwyzKmj2Ew74PBvoZjVMwRPX7jfGbi1If9nTTGWx69Tq9t4+QWiG6sGex++5ovQYJFjjQMD/Z1mEifOcMoFmkqPz6T3QR7RuITgWqaFhNYUUOq0XBmNJDqyxUVwUAkmCKK3VdF6oZz1VLMNn1tHv2ZlmIbMUXp/BBEPYYm6i8XyHjRNA3hYBjBriDCgXD8BtNBNqIO1CgrSqKxjn5f78FhspsGnetlItEb7Y3hicJazHBFavh8Q2JYjZlvrNffKc0AZzsjPh0aSUxB02IbxKUejeN9NNwb0/kAomSWxvVzoL9X9cb26EZ64/751yfeOvECu0MrwMjKPlCdMtHql4nAaFCPDrKc/y4ZKDijhRPbitNzOC192CF9yChokaE89Xn6umE1Mmyz3jAsiZLRYKNfRKE3nusNP1OdpmndQ9BqscPT6kOdhZRQdwOzHBleSB/KUr8AJdG/h4iAUfj+n4T0YfFigi09GvB7BWh7NPLrF0/pjdb6Zzq6N1V0UG4sxO3lEH0BX48hP2MeP6K+D84/F9EB/5jh7M5PmqpFhjg0WSOTbI0MQzkB6/3oIdf0+jr6YkXj8Z9/7DG35wNyYbX7uxEAYp/mOJ+bPj5KoihGhomUI98HshDpladfSKq/T8aCe7F7TPYLAIPr959gFosFFkvvHigrL1057RLQExEREQ2XIAgQZGFSNrZPdIIgQLbIgx5Wi2JJJgmSafST4dLY0Ic8k0QJYG7caUEP9kGKXBBCRERENB0l5Ew6IyMDkiShoaEhZn5DQwPc7rGLHBEREREREREREREREQ1VQoIpZrMZy5Ytw9atW415qqpi69atMcN+ERERERERERERERERJVrCxiH4zne+g9tuuw3Lly/HRRddhCeeeAIejwd33HFHoopERERERERERERERETUS8KCKTfddBOamprwwx/+EPX19ViyZAneeustZGdnD7itnjyoo6NjrItJRERERERERERERESTRHJy8pgkuBc0PTIxiZw6dQqzZ89OdDGIiIiIiIiIiIiIiGgCaWxsRGZm5qjvN2E9U0YiLS0NAFBZWYmUlJQEl4aIaOQ6OjpQUFCAqqoqOJ3ORBeHiGjEWK8R0VTDeo2IphrWa0Q01ej1mtlsHpP9T8pgiiiKAICUlBRW9kQ0pTidTtZrRDSlsF4joqmG9RoRTTWs14hoqhmLIb4AQByTvRIREREREREREREREU0RDKYQERERERERERERERH1Y1IGUywWCx5++GFYLJZEF4WIaFSwXiOiqYb1GhFNNazXiGiqYb1GRFPNWNdrgqZp2pjsmYiIiIiIiIiIiIiIaAqYlD1TiIiIiIiIiIiIiIiIxguDKURERERERERERERERP1gMIWIiIiIiIiIiIiIiKgfDKYQERERERERERERERH1g8EUIiIiIiIiIiIiIiKifkzYYMrPfvYzCIKA+++/35jn9/txzz33ID09HQ6HA5s2bUJDQ0PMdpWVldi4cSPsdjuysrLwve99D4qijHPpiYh661mvtbS04L777sO8efNgs9lQWFiIb3/722hvb4/ZjvUaEU1U8X6v6TRNw7XXXgtBEPDqq6/GLGO9RkQTVV/12s6dO3H55ZcjKSkJTqcTl112GXw+n7G8paUFN998M5xOJ1wuF+688050dXWNc+mJiHqLV6/V19fjlltugdvtRlJSEpYuXYo//elPMduxXiOiieKRRx6BIAgxU0lJibF8PGMG8ogfzRjYs2cPfvOb32DRokUx8x944AG8+eabePnll5GSkoJ7770Xn//857Fjxw4AQDgcxsaNG+F2u/Hhhx+irq4Ot956K0wmE376058m4qEQEQGIX6/V1taitrYWv/jFL1BaWoqzZ8/iG9/4Bmpra/HHP/4RAOs1Ipq4+vq9pnviiScgCEKv+azXiGii6qte27lzJ6655ho8+OCD+OUvfwlZlvHpp59CFLuvTbz55ptRV1eHLVu2IBQK4Y477sDdd9+NF154YbwfBhGRoa967dZbb0VbWxtee+01ZGRk4IUXXsCXvvQl7N27FxdeeCEA1mtENLEsWLAA77zzjvG3LHeHNcY1ZqBNMJ2dndrcuXO1LVu2aGvXrtX+4R/+QdM0TWtra9NMJpP28ssvG+seOXJEA6Dt3LlT0zRN++tf/6qJoqjV19cb6/z617/WnE6nFggExvVxEBHp+qrX4nnppZc0s9mshUIhTdNYrxHRxDRQvfbJJ59oeXl5Wl1dnQZAe+WVV4xlrNeIaCLqr15buXKl9tBDD/W5bXl5uQZA27NnjzHvb3/7myYIglZTUzOWxSYi6lN/9VpSUpL2+9//Pmb9tLQ07b/+6780TWO9RkQTy8MPP6wtXrw47rLxjhlMuGG+7rnnHmzcuBEbNmyImb9v3z6EQqGY+SUlJSgsLMTOnTsBRK4YuuCCC5CdnW2sc/XVV6OjowOHDx8enwdARNRDX/VaPO3t7XA6nUaEnfUaEU1E/dVrXq8Xf/d3f4ennnoKbre713LWa0Q0EfVVrzU2NmL37t3IysrCmjVrkJ2djbVr1+KDDz4w1tm5cydcLheWL19uzNuwYQNEUcTu3bvH7TEQEUXr7/famjVr8D//8z9oaWmBqqp48cUX4ff7sW7dOgCs14ho4jl+/Dhyc3Mxa9Ys3HzzzaisrAQw/jGDCTXM14svvoiPP/4Ye/bs6bWsvr4eZrMZLpcrZn52djbq6+uNdaKfFH25voyIaLz1V6/11NzcjJ/85Ce4++67jXms14hoohmoXnvggQewZs0a3HDDDXGXs14joommv3rt1KlTACJjdf/iF7/AkiVL8Pvf/x5XXHEFDh06hLlz56K+vh5ZWVkx28myjLS0NNZrRJQQA/1ee+mll3DTTTchPT0dsizDbrfjlVdewZw5cwCA9RoRTSgrV67E5s2bMW/ePNTV1eFHP/oRLr30Uhw6dGjcYwYTJphSVVWFf/iHf8CWLVtgtVoTXRwiohEbSr3W0dGBjRs3orS0FI888sj4FJCIaIgGqtdee+01vPvuu/jkk08SUDoioqEbqF5TVRUA8Pd///e44447AAAXXnghtm7dimeeeQaPPvrouJaXiGgggzkP/ed//me0tbXhnXfeQUZGBl599VV86Utfwvvvv48LLrhgnEtMRNS/a6+91ri/aNEirFy5EjNmzMBLL70Em802rmWZMMN87du3D42NjVi6dClkWYYsy9i2bRuefPJJyLKM7OxsBINBtLW1xWzX0NBgDCHhdrvR0NDQa7m+jIhoPA1Ur4XDYQBAZ2cnrrnmGiQnJ+OVV16ByWQy9sF6jYgmkoHqtS1btuDkyZNwuVzGcgDYtGmTMWwE6zUimkgGcx4KAKWlpTHbzZ8/3xhewu12o7GxMWa5oihoaWlhvUZE426geu3kyZP41a9+hWeeeQZXXHEFFi9ejIcffhjLly/HU089BYD1GhFNbC6XC8XFxThx4gTcbve4xgwmTDDliiuuwMGDB7F//35jWr58OW6++WbjvslkwtatW41tKioqUFlZidWrVwMAVq9ejYMHD8ZU+Fu2bIHT6ez145eIaKwNVK9JkoSOjg5cddVVMJvNeO2113pdOcR6jYgmkoHqtR/84Ac4cOBAzHIAePzxx/Hss88CYL1GRBPLQPXarFmzkJubi4qKipjtjh07hhkzZgCI1GttbW3Yt2+fsfzdd9+FqqpYuXLluD4eIqKB6jWv1wsAEMXYJkFJkozeeKzXiGgi6+rqwsmTJ5GTk4Nly5aNa8xgwgzzlZycjIULF8bMS0pKQnp6ujH/zjvvxHe+8x2kpaXB6XTivvvuw+rVq7Fq1SoAwFVXXYXS0lLccssteOyxx1BfX4+HHnoI99xzDywWy7g/JiKa3gaq1/RAitfrxfPPP4+Ojg50dHQAADIzMyFJEus1IppQBvN7Ld6VPYWFhZg5cyYA/l4joollMPXa9773PTz88MNYvHgxlixZgueeew5Hjx7FH//4RwCRXirXXHMN7rrrLjz99NMIhUK499578eUvfxm5ubnj/piIaHobqF4LhUKYM2cO/v7v/x6/+MUvkJ6ejldffRVbtmzBG2+8AYD1GhFNLN/97ndx/fXXY8aMGaitrcXDDz8MSZLwla98BSkpKeMaM5gwwZTBePzxxyGKIjZt2oRAIICrr74a//mf/2kslyQJb7zxBr75zW9i9erVSEpKwm233YYf//jHCSw1EVF8H3/8MXbv3g0ARqI/3enTp1FUVMR6jYimHNZrRDTZ3H///fD7/XjggQfQ0tKCxYsXY8uWLZg9e7axzh/+8Afce++9uOKKK4xz1ieffDKBpSYiis9kMuGvf/0rvv/97+P6669HV1cX5syZg+eeew7XXXedsR7rNSKaKKqrq/GVr3wF586dQ2ZmJi655BLs2rULmZmZAMY3ZiBomqaN2iMjIiIiIiIiIiIiIiKaYiZMzhQiIiIiIiIiIiIiIqKJiMEUIiIiIiIiIiIiIiKifjCYQkRERERERERERERE1A8GU4iIiIiIiIiIiIiIiPrBYAoREREREREREREREVE/GEwhIiIiIiIiIiIiIiLqB4MpRERERERERERERERE/WAwhYiIiIiIiIiIiIiIqB8MphAREREREREREREREfWDwRQiIiIiIiIiIiIiIqJ+MJhCRERERERERERERETUj/8fB02Cti1ALqAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 2000x200 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax=utils.plot_track(pred_mut.mean(axis=1)[:,:,0], color='green', label='TSS+CRE')\n",
    "utils.plot_track([wt], ax=ax, zoom=[400, 500], marks=[447,448], color='grey', label='WT')\n",
    "utils.plot_track(pred_control.mean(axis=1)[:,:,0], color='purple', label='just TSS', ax=ax, zoom=[400, 500], marks=[447,448])\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b54ad9c7-8331-41dc-bd60-59fcacde34b1",
   "metadata": {},
   "source": [
    "# Fine-tile search"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9a4c8ca-49f1-434b-b4b6-e88f21c1585e",
   "metadata": {},
   "source": [
    "In this example we will use the same enhancing CRE of GATA2 from above. We will perform a fine-tile search to identify the smallest subset of tiles that explain the CRE effect of the this particular sufficient enhancer.  \n",
    "The search output depends on these factors which need to be defined by the user depending on the specific aims:\n",
    "1. scale: The test can run multiple stages of sequence pruning/randomization at different resolutions. This ensures that we can run a broad sweep at first to eliminate unnecessary sequences and then run a finer brush search for shorter sequences.\n",
    "2. threshold: fraction explained/recovered below which we will stop the search or search stage. For instance, 0.9 means that once the current version of pruned sequence only recovers 80% of the original activity of the entire CRE tile the search will stop.\n",
    "3. batch size: sometimes we will want to remove large scale windows, one at a time. At other times, we will want to prune smaller-scale windows but we will still want to speed up the process. This parameter determines how many sub-tiles we prune at a time. Note, setting window size to 500bp and batch size to 1 is not equivalent to setting window size to 50bp with batch size of 10 as these 10 sub-tiles can be discontinuous (occur in different parts of the sequence).\n",
    "\n",
    "Note, we need to supply these as lists, e.g. scales=[500, 50], thresholds=[0.9, 0.7], N_batches=[1, 10] will run a fine-tile seach with 2 stages: (i) at 500bp, 1 tile at a time until fraction explained drops below 0.9, (ii) at 50bp, 10 tiles at a time until fraction explained drops below 0.7.\n",
    " \n",
    "Once we set these custom parameters we can run the search by using the following:  \n",
    "\n",
    "- a loaded model\n",
    "- onehot encoded sequence (WT) of the sequence\n",
    "- control sequences. These should contain the TSS embedded already.\n",
    "- TSS activity when the whole CRE is embedded (with TSS) in background sequences\n",
    "- tile start coordinate\n",
    "- tile end coordinate\n",
    "- scales (window sizes)\n",
    "- thresholds (to stop the search)\n",
    "- fraction of step size to move the sub-tile \n",
    "- batch sizes\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "8a0e16fa-611b-4f17-b3d2-8b3630b5989f",
   "metadata": {},
   "outputs": [],
   "source": [
    "track_index = [5111]\n",
    "bins = [447, 448]\n",
    "model = custom_model.Enformer(track_index=track_index, bin_index=bins)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "12d31d33-9680-45d5-835a-f42e27e169dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean CRE+TSS activity across bins 447 and 448 is 193.565185546875\n"
     ]
    }
   ],
   "source": [
    "mut = pred_mut[:,:,bins,:].mean()\n",
    "print(f'Mean CRE+TSS activity across bins 447 and 448 is {mut}')\n",
    "scales = [500, 50]\n",
    "thresholds = [0.9, 0.7]\n",
    "frac=1\n",
    "N_batches = [1, 10]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ff583a5-d5f9-4cf9-b052-c64169781aa0",
   "metadata": {},
   "source": [
    "**What smallest set of sub-tiles recovers 80% of the activity of the GATA2 enhancer?**\n",
    "\n",
    "\n",
    "GATA2 TSS activity is low on its own (control case of just TSS embedded in shuffled sequences) compared to the high activity of the WT. \n",
    "However, when we embed an enhancing CRE (identified previously) together with the TSS the activity goes back up (not fully but to a much higher level than the control case). But do we need the whole CRE for this effect?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "fbcd07e1-dd86-4e5b-9079-ca35e7796387",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tile size = 500, threshold = 0.9\n",
      "Starting score: 1.0\n",
      "10\n",
      "Starting optimization...\n",
      "score = 1.0\n",
      "Number of tiles to test: 10\n",
      "Number of tiles at the end of iteration: 10, score = 1.0085970163345337, bps = 4500.0\n",
      "score = 1.0085970163345337\n",
      "Number of tiles to test: 9\n",
      "Number of tiles at the end of iteration: 9, score = 1.011488437652588, bps = 4000.0\n",
      "score = 1.011488437652588\n",
      "Number of tiles to test: 8\n",
      "Number of tiles at the end of iteration: 8, score = 1.007453203201294, bps = 3500.0\n",
      "score = 1.007453203201294\n",
      "Number of tiles to test: 7\n",
      "Number of tiles at the end of iteration: 7, score = 0.9873055219650269, bps = 3000.0\n",
      "score = 0.9873055219650269\n",
      "Number of tiles to test: 6\n",
      "Number of tiles at the end of iteration: 6, score = 0.9530054926872253, bps = 2500.0\n",
      "score = 0.9530054926872253\n",
      "Number of tiles to test: 5\n",
      "Number of tiles at the end of iteration: 5, score = 0.9063980579376221, bps = 2000.0\n",
      "score = 0.9063980579376221\n",
      "Number of tiles to test: 4\n",
      "Number of tiles at the end of iteration: 4, score = 0.8719760775566101, bps = 1500.0\n",
      "Tile size = 50, threshold = 0.7\n",
      "Starting score: 0.9063980579376221\n",
      "40\n",
      "Starting optimization...\n",
      "score = 0.9063980579376221\n",
      "Number of tiles to test: 40\n",
      "Number of tiles at the end of iteration: 40, score = 0.9034410119056702, bps = 1500.0\n",
      "score = 0.9034410119056702\n",
      "Number of tiles to test: 30\n",
      "Number of tiles at the end of iteration: 30, score = 0.8665218353271484, bps = 1000.0\n",
      "score = 0.8665218353271484\n",
      "Number of tiles to test: 20\n",
      "Number of tiles at the end of iteration: 20, score = 0.7589319348335266, bps = 500.0\n",
      "score = 0.7589319348335266\n",
      "Number of tiles to test: 10\n",
      "Number of tiles at the end of iteration: 10, score = 0.09304965287446976, bps = 0.0\n",
      "CPU times: user 39.3 s, sys: 8.79 s, total: 48.1 s\n",
      "Wall time: 7min 45s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "optimization_results = creme.prune_sequence(model, \n",
    "                                          wt_seq, \n",
    "                                          control_sequences, \n",
    "                                          mut, \n",
    "                                          cre_tiles[0][0], \n",
    "                                          cre_tiles[0][1],\n",
    "                                          scales, \n",
    "                                          thresholds, \n",
    "                                          frac,\n",
    "                                          N_batches)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7732f682-1604-4de8-a7ef-7f5792d691ab",
   "metadata": {},
   "source": [
    "Breakdown of results:\n",
    "The run above went through several iterations of pruning and testing if the threshold of acceptable fraction explained has been crossed. \n",
    "At iteration 20 (the second to last iteration) we get a score which is still above the threshold. The last iteration outputs a score lower than the threshold. To now re-create the final sequence of 500bp sub-tiles we can use the 'insert_coords' values in the results as shown below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "id": "a0d2eef1-d481-43cc-a75f-7e40cbaecca1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.75893193"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "final_sequences = control_sequences.copy()\n",
    "for embed_start, embed_end in optimization_results[50]['insert_coords']:\n",
    "    final_sequences[:, embed_start: embed_end, :] = wt_seq[embed_start: embed_end, :].copy()\n",
    "final_prediction = model.predict(final_sequences)\n",
    "final_prediction.mean() / mut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "b6ae897f-0c1a-4caa-8a14-73a40255bb15",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bbf3e1e-9e4c-4efc-a6b8-58483b33a371",
   "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
}