{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing predicted occupancies, information content, and related metrics\n",
    "This notebook takes all sequences in the libraries and computes the predicted occupancy of 8 TFs. Using those predicted occupancies, we then calculate the total occupancy, TF diversity, and information content (entropy) of each sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from IPython.display import display\n",
    "\n",
    "sys.path.insert(0, \"utils\")\n",
    "from utils import fasta_seq_parse_manip, modeling, plot_utils, predicted_occupancy, sequence_annotation_processing\n",
    "\n",
    "data_dir = os.path.join(\"Data\")\n",
    "figures_dir = os.path.join(\"Figures\")\n",
    "all_seqs = fasta_seq_parse_manip.read_fasta(os.path.join(data_dir, \"library1And2.fasta\"))\n",
    "# Drop scrambled sequences\n",
    "all_seqs = all_seqs[~(all_seqs.index.str.contains(\"scr\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_utils.set_manuscript_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute predicted occupancy of all TFs on all sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in PWMs\n",
    "pwms = predicted_occupancy.read_pwm_files(os.path.join(\"Data\", \"Downloaded\", \"Pwm\", \"photoreceptorAndEnrichedMotifs.meme\"))\n",
    "ewms = pwms.apply(predicted_occupancy.ewm_from_letter_prob).apply(predicted_occupancy.ewm_to_dict)\n",
    "mu = 9\n",
    "\n",
    "# Do predicted occupancy scans\n",
    "occupancy_df = predicted_occupancy.all_seq_total_occupancy(all_seqs, ewms, mu, convert_ewm=False)\n",
    "occupancy_df = occupancy_df.rename(columns=lambda x: x.split(\"_\")[0])\n",
    "# Save to file\n",
    "sequence_annotation_processing.save_df(occupancy_df, os.path.join(data_dir, \"predictedOccupancies.txt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute total occupancy, diversity, information content"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "entropy_df = occupancy_df.apply(predicted_occupancy.boltzmann_entropy, axis=1)\n",
    "sequence_annotation_processing.save_df(entropy_df, os.path.join(data_dir, \"entropyDiversityOccupancy.txt\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot predicted occupancy vs. relative $K_D$ for various values of $\\mu$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAEUCAYAAABqA/5qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXl8ldW1v599huRknicgE1OYA0mQSURARKyzVetUe3sr9LZWb+3tpdpq1VZb21vv716HVlBrtdpqL85UHFAmGSQJYwhhykjmeU7OsH9/vDkBMUBC3vO+gezHTz7hvNlnrx3fnHXW+e611xJSShQKhULhOyxmL0ChUCgudJSjVSgUCh+jHK1CoVD4GOVoFQqFwscoR6tQKBQ+RjlahUKh8DHK0SoUCoWPUY5WoVAofIxytAqFQuFjbGYvQC+io6NlSkqK2csY1uTk5ACQmZlp8koUZyMnJ6dWShlj9jqGC+JCOYKblZUls7OzzV7GsEYIAcCF8jd1ISOEyJFSZpm9juGCkg4UCoXCx1ww0oHCfFQkq1D0jYpoFQqFwscoR6tQKBQ+xjBHK4S4RwiRLYToEkK8fJaxPxZCVAohmoUQLwkh/A1apmIQZGZmqowDhaIPjNRoy4FfA0uBgNMNEkIsBX4GLOp5ztvAoz3XFEOY3Nxcs5egUAxJDHO0Usq3AIQQWcCoMwy9C3hRSpnXM/5XwGsoR6u4gDndRqLkTBuMp/+ZRVgHuSKFngzFrIPJwLsnPd4DxAkhoqSUdSatSTFE8EgPHe5W2lzNtLmaaHe30OVup8vTidPTSZenk25PJ92eLro9nTg9XbilC49045ZuPNKFW7pxS9eJ67jxSDcSD1JKvP9x0r9PfuyR8mtjNYfo/TKX5vYYnp79vNnLUJzEUHS0wUDTSY+9/w4BvuJohRDLgeUASUlJhixO4XuanfVUdRZT311FfVcVDU7te7OzjnZ3KxKP2Us0lDNlzfX1oy6X22drUZwbQ9HRtgKhJz32/rvl1IFSylXAKtBOhvl+aQq9cXq6KGkr4Fjbfo53HKWio5BWV+MZnxNgDSbIGkqQLZRAWygOayB+FscpX/74WRzYLf5YhQ2rsGLxfsfac82GxMLhplr21ldysLGWgqZqyttbQJ5wYhJx0mNBsM2faEcQIXYHwXYHIXYHIXZ/gv0chNj8CbT542+x4W+14We1nvi3xYrDZsdusWITFixCIBCU1TdxsLyG/OPVHDheTUVjM0jhNf617yEOf8ICHIQGOAgLdBAW4Eeow0GIw0Ggv50FMZG63yfF4BiKjjYPSAfe7HmcDlQp2eDCoc3VRF7TDvKatlHUdgCXdH7l5w5LIAkBqUT5JxDhF0ukXzwRfrGE2aMJsoViFYP7s213dbO+/DAfluaztaqIFmfXV34eYA0kLSyGlJAoUkIiSA6OZERgKHEBIcQ4gnDY7IOy3+1ys/VwMZ/lH2VD/jHqWtu/8nO71UZKdARjYiMZHRPJyMgwEsJCiAsLJi4shEC/wdlXGI9hjlYIYeuxZwWsQggH4JJSuk4Z+grwshDiNbSsg18ALxu1TsW5c/fdd5/2Zx7p4UjrbrbXfsihll29H/8FggRHKmOCp5IYlMYIRyoRfnG9dRP0ZG9dOX85vJOPygrocJ9w7snBEcyLS2V61EimRiYwJiQKq0X/zMfCmnrW7NzPO7kHaGjr6L0eFxpMRspIZiSPYEZyAuPio7Fb1WbWhYRhRWWEEI8Avzzl8qPAS8ABYJKUsqRn7P3ASrQ0sDXA96WUXZwBVVRmaOKRbnY3bGRD9RrquisAsGBlbEg6U8LmMCF0JkG20LPMcu5IKdlQcYQ/5W8ju7a09/qMqJFclTSZy0aOY1RQuM/sAxRU1PDc+u18mnek99q4uCgunzqeRZPGkBYf7ZM3ljOhisoYi6repfAZBc25fFT5KlWdxQCE2aOZFXUFWZGLCbKF+d5+YzW/3v0JW6uKAAix+/OtMTO4fUwmicG+da4A5Q3N/NeHm/ho32EA/G1Wrpo+kW9eNJWpo3wTtfcX5WiNZShqtIrzFG892rRpY3ivfBV5TdsBCLfHcFn8rUwLn4/VgPzOTpeT/9q3gb8c3olHSsL8HPxg4jy+NWYGwXbfHzL0eCSvfJHLM59spcPpws9m5eaLpvG9BVnEhAb73L5i6KEcrUI3srK0AOnX+++i3d2Mn8XBoribmR11JXaLnyFrONhYxX3b3uFIcy1WIbhr3EzunTyfcP/THkbUlZrmVh74x0dsO1ICwNKp4/jPbywgPizEEPvd3S78/NTLeqih7ohCF06WoNrdzYwOnsoNo35IhF+sYWv4Z2k+/7njfTrcTsaERPGH2dcyNTLBMPv7Siu559V3qW1pJzIogMduXMLCiWMMs19Z1cQDD77JzTfNYtkV0wyzqzg7ytEqBo1bunm37E+9jxfHfYtLY7+JRRhXHG7VwW08ueczAK5PmcqvM5cNOg1rIHyad4SVb3xIp9PFzNGj+P0tywyVCY4ereJnD/6DuvpW3n43h8uXTMFqVcX5hgrK0SoGhVu6eKP4KfKat/deWxR3s2H2pZT8v/2beObAFgTwwPTL+O74iwzdaPpwbwEr3/gQt0dyY9YUHrpukaHpWfkHy/nPlW/Q1t7F9PQkHnvkBuVkhxjK0SrOGY90s6b0afKat+OwBpmyhmcPfMEzB7ZgEYLfXXQ116dMNdT+R/sO9TrZ5Qsv4t4lcw118keOVLHyAc3JXjI/jQd/drXSaIcg6m1Pcc58UP4Sexo342dx8C+pDxtu/42ju/jv/RsRwP+bfZ3hTja36Dg/e3Mdbo9kxcJZhjvZyqomVj7wBq2tXcybN45fPHiNcrJDFOVoFefEjrqP2FH3ITZh586UBxkVOM5Q+1urivhFzocAPJp5Bd9ImmSo/eLaRn706nt0u9x8a3Y6P1oyx1An29HRzS8e/j8aGtvJzEjhoQevxWZTp8mGKurtTzFgitry+eD4CwBcN+r7jA6eAoBRB0Yq25u5b9vbeKRkxYQ53D7W2K4OXU4XP379AxrbO7kkLZUHrrrUUCcrpeS3v/uAY8dqGDUqkl8+dJ2KZIc46u4oBkSHq5U3Sv6ABzfzoq9mRsTC3p8Z0cbG5fHwo61vU9/Vzry4VH4y9VKf2zyVJ9dupKCihqSocH7/rWXYDN54ev+D3WzecoigIH8e/9U3CQ52GGpfMXCUdKAYEO8dX0Wzs57EwPEsTfi24fZXH9xGbl0Z8QEh/Pfsa31S/OVMrM87whs79mK3Wnnqtm8Q7DC2nV1xcS1/fH49AD++bymJo1RJxPMB5WgV/WZv4xb2Nm3Bz+LgpsT7vnacdvny5Sxfvtxn9gsaq/mfvE0A/O6iq4lyGJvp0NTRyWPvak7uJ8suZuII4w5jALjdHn77+7V0dbm4fMkUFi00VpdWnDvK0Sr6RYe7jbXlLwGwLOE7RPl//cTV6tWrWb16tU/suz0eVn75AU6Ph9vGZDAvPtUnds7EH/65mdqWdmYkj+D2OTMMt//e+7soKKggOjqEH/1wieH2FeeOcrSKfvFp5eu0uhpJDpxAVuRlhtv/R+Ee9jVUEB8Qwsr0RYbbzyk6zprs/ditVh67YQkWi7GVt2prW3jxzxsB+NEPLyMoyFjJQjE4lKNVnJXyjkJ21H2EBQvXjFxu6NFagKbuDv6wbwOgnfwyogLXyXg8kic/0Jzc9xZkMTrWeF101QsbaG/vZs7ssVw8b7zh9hWDQzlaxVn5qOJVJB5mR19JfECK4fafzttCfVc7F8Uk8Y3EiYbbX7vnIHnHq4gNDeK7C2Yabv/IkSrWf5aH3W7lnh9cZmodW8W5oRyt4owcbd3HkdbdOCyBLIy9yXD75W1NvHYkBwE8NGOJ4U6m0+ni/320BYD7Lp9nSr+u1S9uQEq49poMEhJ8X7BcoT/K0SpOi5SSjyv+CsD82OsItBlTU/Vknj3wBd0eN99ImsSkiHjD7a/ZuY/KplbSEmK4Zobxu/y5u4rYmV1IUKA/t986x3D7Cn1QBxYUp+VgSzZlHYcJtoUzN/qqs47PyMjQ1X5JawP/V7gHixDcN/kSXefuD11OFy9s3AnADxfPNnwDDOCvr20F4OabLyIsLNBw+wp9UI5WcVo2Vb8FwPyY6/CznP30kbeVjV78MX8rLunh+pSpjA6N0nXu/vB2Th7VzW2Mj482tIC3l/15ZezeU0JQkD/XX2fsMWOFvijpQNEnRW35lLQXEGANZmak8TmbtZ2tvFO0DwH8cNI8w+073e7eaPb7i2aZEs2+/vdtAFx3bSbBQeqY7fmMcrSKPtlU/TYAs6OW4W81pt/Wybx6OIduj5vFI8eTGmJ8NLs+7wgVjS2kxkSwZLKxlckAjhVWs337Ufz9bdx4vWpWe76jHK3ia9R0llHQko1d+DEn+sp+P08IoUtWQIfLyWtHNBnie2mzBj3fufDq1l0A3D5nhinR7NvvaL//FUunER6utNnzHeVoFV9jR91HAKRHLCDIFma4/beL9tHQ3UF65AiyohMNt7+vtJLdxRWEOPy5JsP4vN2Wlk7Wf3YAgOuvVdrshYBytIqv0O3pZFfD5wDMirrCcPtSSv52NBeA74yfaUpy/ms90eyNM6cQ5G9Mm/STWffxXjo7nWRmpJCUZLxsotAf5WgVX2FPw2Y6Pe0kBaYxIsD4wi376is40FhFhF8AS0dNMNx+Y3snH+0/jBBw6+x0w+17PJL33tMc/XXX6psupzAP5WgVX2FH3TrAnGgW4G/HNCdzQ+o0/K3GZx+u3Z1Pt8vNnLHJjIo0XjbZvbuY4+UNxMaGMnvWWMPtK3yDcrSKXio6CqnoLCTAGszkMONPIbU4u/igJA+Ab402vgyhlJI12Zr9G7MmG24fYN3H+wC4YulU1TL8AkLdSUUvuxo2ADAtfD52i/Ha5D9LDtDucnJRTJIpBxTyy6spqKghPNDBoknGH1Boa+ti85YCAC5fYmxHX4VvUSfDFAC4pZvdjVr3ghkRl57THM8///yg1vBu8X4AbkydNqh5zpW3eqLZq6ZPxM9m/Etj0+aDdHW5mDY1kRGqeMwFhWF/TUKISOBF4HKgFnhASvl6H+P8gf8BrgfswBfA96WUx41a63DkSMtu2lxNxPiPZFTAuWmDg2ljU97ezI6aEvytNlM2wZxuN+v2HQLgepNlg6WXq2j2QsNI6eBZoBuIA24H/iiE6Osv+j5gDjANGAE0AE8btcjhilc2mBFhbOtsL+8Xa9Hk4hHjCDG4sDfA9iOlNLR1MDo2krT4aMPtV1Y1sW9fGQ6HnQWXpBluX+FbDHG0Qogg4EbgISllq5RyC/AecGcfw1OBj6SUVVLKTuANwJwQY5jQ7enkYLN2rj89/NyrZK1atYpVq1ad03O9ssE1yebc6n/uOQjAldPSTHmj2bhJsz9n9lgCA1WbmgsNoyLa8YBLSnnopGt76NuBvgjME0KMEEIEokW/H/Y1qRBiuRAiWwiRXVNTo/uihwsFzTk4ZTdJgWmE+8Wc8zwrVqxgxYoVA37eoaYaCpqqCfNzsCDe+JSmTqeL9QeOAnBlujnR5IaNmqNdcInxsonC9xjlaIOB5lOuNQF9VZI+DJQCx3ueMxF4rK9JpZSrpJRZUsqsmJhzdxDDnf1NWs3TKWFzTbG/rjQfgKUj0/CzWs8yWn82FRTS1tXN5JGxJEdHGG6/oqKRgoIKHA47sy4abbh9he8xytG2AqGnXAsFWvoY+yzgD0QBQcBbnCaiVQyebk8nBc1aAZPJYbNNWcNHZVpK09JEc6K5j/cdBuCKaeZEs17ZYO6ccfj7G98qR+F7jHK0hwCbEOLkenPpQF4fY6cDL0sp66WUXWgbYRcJIYzfoRgGHGrOxSm7SQwcNyjZ4FwpaqnnYFM1wXZ/5sSmGG6/2+ViU0EhAJdNNuckllc2uHSBkg0uVAxxtFLKNrTI9DEhRJAQYh5wLfBqH8N3At8WQoQJIezAD4ByKWWtEWsdbuxv0opLmyUbeKPZxSPGmXLk9stjZbR1dTM+PpqkKONzV6urmzl0uBKHw87MLONrSyiMwcj0rh8AAUA18Dfg36SUeUKI+UKI1pPG/QfQiabV1gBXouXUKnTG5XFyqEWrlDXJJNng4zItmls6ypyP7Z/1bIKZcRIMYNv2IwBkZaYq2eACxrAQQkpZD1zXx/XNaJtl3sd1aJkGCh9T1HaALk8Hsf6JRPrFGW6/qqOF3fXlOKw2Lok33tF5PLLX0ZolG2zdpunDc+eoAjIXMuoI7jDmYHM2ABNCZ+oyn5RyQOM3lGvR3MVxqQTYjI/m9pdVUtPSRkJ4CBMSjNen29q62LW7GItFqEpdFziqqMwwRUrJwRavozWnJ9VnFZqjvXSEOU7m84PHAFg0cYwphxR2ZhficnmYPGmkaldzgaMc7TCluquUhu4qAq2hJAYa33ywy+1ia5W2239pgjmOdnNBEQCXTDBnE8orG8yZY/z/f4WxKEc7TPHKBmmhmViEPocEMjMzyczsX4+rL2tKaHc5mRgeS0LgqSnWvqemuZX88mocdhszU0cZbt/jkWTnaG80c2aZsxGnMA6l0Q5TDrdonQzSQvRr/pebm9vvsZ/36LNmRbNbDhUDMGt0Iv52418GR45W0djYTmxsqOoLNgxQEe0wpMvdQUl7AQILY4LNKcn3eY8+u2iEOR+bNx/SosmL01JMse+NZrMyU3XTh5vrW3j23pfY8c/+v+EpjEE52mFIYVsebuliZMAYAm19lZvwLSWtDZS0NhDm5yA9coTh9l1uD1sPlwAwfwg4Wr3Y/dl+3nnmQ978/bu6zanQB+VohyFHW/cAMDbE+C6vAFurigCYE5uC1WL8n+CeknJaOrtIjYkgMdL402AdHd3s31+GxSLImJGs27y5n+wFIOMyczpUKE6PcrTDkCMt2gtybLA5jvaLnmyDeXEpptjffrQUgLlj9XNyA2HP3hJcLg9p4xMIDQ3Qbd7c9VqHhswlytEONZSjHWY0dddS3VWKn8VBYuB4w+17pGRbT0Q7N86ctKrtRzXZYPbYJFPsZ2d7ZYMU3easOFZFxbEqQiKCGJepSi0ONVTWwTDjaKsWzY4OmoLNou9prLvvvvusY/Ibq2jo7mBkYBjJwcbXfm3r6mZvSSUWIchKHWm4fYDsnCIAsnQsIpPTIxtMXzQFqwk1fRVnRjnaYcZhH+qz/Wljc7JsYMZprNyiclweD1NHxRMa4DDcflV1EyWldQQG+jFxgn4bgbmfavc14zJz5CDFmVHSwTDCIz0cbelxtMHTTVmDdyNsbrw5ssGOHtlg1phEU+zn9ESzM2YkY7PpE3m63W52f6b1XMu4THXQHYooRzuMqOwsos3dTJg9mmh//dOqcnJyyMnJOe3Pu9wudtZojs6MIt9wYiNs9lhzHG3uLu2gRFaGfm80R3YV0dLQRnxqLCPGxOs2r0I/lHQwjDjWqu1Kjw1O98nH9qwsrTjN6ap45daW0el2MTE8lmhHkO72z0ZjWwcHK6qxW63MSDZen5VSsmeP9kaTnq7fRtzejQcAmH6pahY9VFER7TDiWKvWOSg12JwXpFefNSvb4MvCMqSE6ckJOEw4dlt2vIG6+lbCwwNJ1vHY7d6N2n2dtkA52qGKcrTDBI90U9ymRT6pQea8ILdXax+b55qVP3ukJ61rjDlpXXv2avanTU3U7ROF2+1m32ati/C0BZN0mVOhP8rRDhMqO4vp9LQTYY81pQlju6ubffUVWlpVtDn66JfHevRZkzbCvLLBdB1lg8K9JbQ1tROfEkNcsvH3VdE/lKMdJhT2yAYpJskGu+qO45IeJofHE2z3N9x+dXMrhTUNBPrZmTzK+LY9Ukr27NUcffo0/Rztng3afZ2qotkhjXK0w4TCth59NsicF6Q32+CiGHM+tucUHQdgRvII7CYk9JeXN1Jb20JYWADJydG6zbt3kyYHpSt9dkijHO0wwCM9FJmsz35Zo0VzM2NNyl8t1BxtpkmnwU7WZy0WffRZj8fDvh5Hq/TZoY1K7xoGVHeW0uFuJdQeSYQPu91mZ2f3eb3L7WJ3nebozNJnvRFtZoo5jna3N61LR9mgcF8JLQ1txCRGEZ8Sq9u8Cv1RjnYYcEI2mOzTY6+na2Ozr76CLreL8WExRPgb34SwuaOTw1W12K1Wpo4yPqFfSslerz7rg/zZ9Et9e18Vg0dJB8OAoh5Hm2KSbGC2PruruBwpYcqoOFPa1lRUNlFd00xoiIPUFP0yA7z67LRLlGww1FGO9gJHSkmhQfrs8uXLWb58+deuf9njaGea5GizC82VDfb26LNTddRnpZRKnz2PUI72Aqe26zhtriaCbeE+qW9wMqtXr2b16tVfuebyeMitLQNgZoxJ9QVM1mfzDmj2p0zRr9vu8cMVNNW2EBEXpuobnAcMyNEKIYKFEKOEEMG+WpBCX4rbCwBIDppgio6X31hFq6ub5OAI4gKM70/W6XSx/3gVQmhHb83gwIFyACZP1M/R523V7uukuWlKnz0POKujFUJMEUI8LYQ4BjQBJUCTEOKoEOIZIYSqyzaEKe1xtEmBaabYN1uf3VtagcvtYXx8jCn1Z1vbOikqrsFmszB+vH6RZ/62QwBMmm18lwzFwDmjoxVC/B14HagA7gCiAb+e73cCx4HXesYphiAlbZqjTTTJ0eb0yAampXWZrM/m51cgJYwbG4+fn34bcQe29zjauebcV8XAOFtE+1cp5TQp5RNSyq1SygYppavn+1Yp5W+klNOAv57NkBAiUgjxthCiTQhRLIS47QxjM4QQm4QQrUKIKiHEfQP9xRTQ4W6jpqsMq7AxIsD4PlJSSnbVaY52RrQ5js6rz5rVtuZAvmZ/0iT99PG2pjaK9pdis1sZr/qDnRec8S1WSvlBfybp57hngW4gDpgOrBVC7JFS5p08SAgRDawDfgz8H1oErd8uwjCirP0QEsmIgNHYLX6G269ob6aqo5UwPwepIfqVBewvbo+HPaWVAGQk+3Yj8HTk5WmOdvIk/f6ED355BCklY2ek4ucw/r4qBs4ZHa0Q4rv9mURK+dJZ5gkCbgSmSClbgS1CiPfQ5IefnTL8fuAjKeVrPY+7gPz+rEPxVbyygVH6bEZGxlce5/acBpseNRKLCRs2R6rqaOvqZmREKDGhxu/fejyS/IM9G2E6RrQHvPrsHCUbnC+cTTS686R/C2AeUAmUAolo0ekXwBkdLTAecEkpD510bQ+woI+xs4F9QoitwFhgB/BDKWXJqQOFEMuB5QBJSeZstgxlStq1/91G6bOntrHZ1aPPzogyqb5ASQUA6UnmZBsUl9TS1tZFbEwoMTGhus17YFtPxsEctRF2vnBGjVZKudD7BewDfiqlTJRSzpVSJgI/7bl+NoKB5lOuNQF95fuMAu4C7gOSgELgb6dZ3yopZZaUMismRtXiPBmP9FDa42iTAs15Qeb26LMZ0eYoP7uKtWhyukmO1ps/q6c+6/F4yN9+GICJytGeNwxkG9SbdXAyzwC1wL1neW4rcOpbeijQ0sfYDuBtKeVOACHEo0CtECJMStk0gPUOa2q6SunytBNmjybMT7+yfP2l0+XkQEMVFiFIjzRHHzU7oj3Q62j1i+hLDx6nramdmFFRxCYaf18V58ZADixUAteccu1qoLofzz0E2IQQ4066lg7k9TF2L3Byd7++O/0pzkhJm/HRrBCiN3l+f0MlLulhfGiMKYW+G9o6KK5rxGG3kZZgzqcdb0Q7WUdHm7dVu68qmj2/GIijvRf4ixBiqxDiDSHENuAvwI/O9kQpZRvwFvCYECJICDEPuBZ4tY/hfwauF0JMF0LYgYeALSqaHRgl7ebmz+aanNbljWanjIozpdB3c3MHpaX1+PnZGDtGv9KU+T367GS1EXZe0W9HK6X8BEgF/gjk9nwfLaX8uJ9T/AAIQIuA/wb8m5QyTwgxXwjRepKdz4AHgbU9Y8cCp825VfSN19EmBU0wxf6u2p6OBlHm6LOmywb5mj6cNj4eu10/R5+3TUW05yMDOqoipayj7yi0P8+tB67r4/pmtM2yk6/9Ec2RK86BdlcLtV3HsQk/Ehwphts/+aBChkkR7a4S70aYSfmzPtBnm+tbKD14HLu/nbEzUnSbV+F7+u1ohRCpwONohw1OdYwqt2oI4c02GBkwGpvFbrj9srYmajrbiPALICU40nD7LreH/T0HFUzfCJuon6M/uOMIAOOzRmP3M/6+Ks6dgUS0rwNHgZ8A7b5ZjkIPzNZnTz52a0ZlqcNVtXQ4XSRGhhEVbHxHB7fbw8ECTbrQcyPsgLdilyokc94xEEc7GZgnpfT4ajEKfTihz5q0EWayPru7WHNyZuXPFhbW0NHRTUJ8GJGR+p1IU4Vkzl8G4mg3ATOAnLMNVJiHR7opa9cS2o0ujfj8888DsLbO7BNhmj6bbpY+21NIZvJk/d5o3G43B3do91WdCDv/GIijLQLWCSHeRsup7UVK+bCei1KcO1WdJXR7OomwxxJijzDU9vLly2l3dfNfb/0XFiGYZtJBhd09GQdmFfr2FpLRU58t2l9KR2sn8SkxRMYbe18Vg2cgjjYI+ACwo9U5UAxBeusbmCQb7KuvwC0lk8LjCLIbX1mqrrWd0vomAvzsjIsz5+TUidKIOuqz25RscD7Tb0crpfwXXy5EoQ8lbQcBczoqrFq1io0VR2FCqOmFZKaMisNmNb4lXkNDG+XljTgcdsaMjtVt3vwefXai2gg7LxlIetdpKwxLKY/psxzFYDlRSMZ4R7tixQoARv/916YVktldYm4hGW80OyEtAauOjr63R5iP9NmcnJxYm832AjAF1bR1oEghRJPb7f6zx+P5Y2ZmZvepAwYiHRxBqztwcr6Otw6B8WccFV+jzdVEXXcFduFHfECyqWuZYZKj9Ua05h1U8Naf1S+ib6xpovxIJY5Af0ZP8819tdlsL8THx0+MiYlpsFgsqr7IAJBS0t3dbS8vL/9Rc3NzBlr1wa8wkCO4Fimltee7BRgBrOKrNWsVJuLVZ0cGjsUq9OtPNVAi/QNJCgo33K7T7WZ/WRUwBA4q6OhovWURx88cg83Rz3mYAAAgAElEQVTus/s6JSYmplk52YEjhMDf39+ZnJzcBFzc15hz/oggpawE/h34zbnOodCXEx0VzKlv4CUjypyDCgUVNXQ6XSRHhRMRFGC4fZfLTcEhLaLWM+Og96CCbwvJWJSTHRw9///6/HQ/WC0mDTD+6I2iT060Fjd3w8Rs2cCsaPbo0Wq6ulyMGhVJWJh+L4vegwoqf/a8pd+OVgixuaczrfcrG63NzFO+W17/ycnJ6a2HeurXqlWresetWrXqtONOjcIyMzNPO2758uX9si2E+EqLl+XLl592XGZm5lfsn2nOvn6nu8f+mifS32ZS+CxTficv/zZpnm6/00Du0x3zMjjw2/v57S3LTLlPEyaM4PNPH+DVl1fo+rf31IbH+FT+H3Ovmanb76QwloEIPi+c8rgN2COlPKzjehQKheKCQ0h5YcgyWVlZMjs72+xlmMa22n/yQfkLTA9fwE1J9xlvv6qIOza8xpSIeN69/F8Nt1/T0salT6wi0M/O9l/+AKvF+AylW+94jqqqZl5Y9V1Gp+qTQ/vOMx/y7L0vseTbC/jPl+/RZU4AIUSOlDLL+3jPnj1F6enptboZOE944oknYl5//fXoQ4cOBVx99dX1a9asKRrMfHv27IlOT09POfX6QKQDuxDiUSHEMSFEZ8/3R4UQqrH8EKDU7EIyvfUNzNVnpybGm+Jka2tbqKpqJijQn+Qk/U6kneh4q06E+YKRI0c6V65cWXHTTTf59E1mIH+RvwMuA76P1u/r+8Ai4EkfrEsxQIxuLX4q3o4KZhX63t1bSMbkgwoT9D2okL9NbYQBrFy5Mv62227rrXtdU1NjtdlsGe3t7YMSnO+6667GO++8szEqKso1+FWenoFotDcB6T1dFgAKhBC5wB7gx7qvTNFvWpwNNHRX4WdxEOcwvga7lJLddccpe+A5fh72D67ZtdvwNZw4qGBWa3H9DyrUVTRQWVRDYEgAyTpWAusPSyw3ZZ591OD5xPOPflUDzMvLC7zkkkt6u2Zv3749MDU1tSswMPAr2ufChQvHZmdn91mbMisrq/Xzzz8/MrgVnxsDcbSne+dQW5gm460/OypwHFZh/CG9otZ6Gro76C4sZz/lhtvvdrnJ8x5USDTL0erf8dZbSGbCrLFYTWgwOZQ4ePBgwI9//OMq7+Pc3NyAiRMnfq0BgVmO9GwMxNH+A3hfCPEoUAIkA78A3vTFwhT9x8z6BgA5tWWm2PVSUFFDl8tNakwE4SYcVOjudnH4sFY5dKKOBxW8HW/NKCTT30jTCDo7O0Vpaan/RRdd1OG9tnfv3sD09PTzptPLQBztf6I51mfRjt+Wo3Wz/bUP1qUYAN4TYYkmHVTw6rNmYfZBhcNHqnA63aSkRBMc7NBt3hMHFYb3RtiuXbscsbGx3SEhIR4Aj8fD9u3bQ26//fb6U8decskl484kHWzatMmUdNSBlEnsBh7u+VIMEVweJ8c7tE9LZkW03owDs9htdmvx3kaM+skGzm4nh7K1ongTZ4/Tbd7zkdzc3ID6+np7Xl6ef3JycveDDz6YUF5e7jd27NiuU8cO1JE6nU6cTqdwu93C7XaL9vZ2Ybfbpd2ub/PLAW2PCiEWCSFWCyHW9nxfrOtqFAOmorMIl3QS7T+SQFuI4fZbujs53FSD3YSUKi/mV+zy6rP62T+yqwhnl5OkiSMJidCv79j5yL59+wLnz5/ftHjx4rTRo0dPDQkJ8cTFxTkfeeSRQb+zrly5ckRQUFDGc889F//uu+9GBgUFZaxcuVL3P6SB1KP9CbAS+DOwC0gCXhdC/E5K+Qe9F6boH2bXN9hdX44EJofHU2CC/ZrmVsobmwny92NMrPGtzaWUva1rJk/SLzNAdbw9wYEDBwK++93v1qxbt6637vXjjz9eeabn9Jennnqq/KmnnvL5Du5ANNr7gUVSyv3eC0KIV4FPAOVoTcLbUcG8/FlNNsiIHkXU3Xcbbt8rG0wz6aBCVVUTdfWthIY4SEzUz9HnbdXu6+R55lZiGwoUFBQETJ06tdPsdQyGgRa3PDV14hgnin8rDEZKSXG79oJMDjLnBZlb19NaPHoUPz+pgIpRmL0Rtt/biHGyfqUhpZTkfaFFtJPnDe+NsJqaGmt9fb1typQpX9NjzyfOGAIIISzeL+AR4EUhxDghRIAQYjxa4e9fGrBORR80OmtodtbjsAYR42/80VePlOzyOlqTeoTtNvugQo+jnaLjgYLKomrqKxsJjQph1HhzdOehQkxMjNvpdOb6+/uf1wHd2SJaFyciVu/b9a2nXLuNr1f2UhjAiULfaViE8R+bDzfV0OrsIiEwlITA0N6SfKeWEfQV3S43ece1HPZppkW0mnSi60GFrd6Ot+NVScMLhLM52lRDVqE4J7yygVkdFbzRbEZPIZmsLK0YlFEV4Q5W1NDtcjM6JpKwAP3yV/tLW1sXhUU1WK0WJqTp5+jzvujRZ4d5/uyFxBkdrZSyWC9DQohI4EXgcqAWeEBK+foZxvuh1VEIkVKaUxJqiOPdCDNNn+3ZCJthViGZYnMLyeQfLMfjkUyYEI+/v355l3nbvPqs2gi7UDibRvuUECL+LGPihRD96bLwLNANxAG3A38UQkw+w/ifAjX9mHdY0uXuoLKzGAtWRgWak9DujWgzTS6NaJaj9ebP6qnPtjW3U7SvBJvdyvis0brNqzCXs0kHBcCXQoh8YGPP4xYgBBgPXIrWN+yMx3CFEEHAjcAUKWUrsEUI8R5aB92f9TE+FbgDLaVs9QB+n2FDafshJB5GBIzFz+JvuP2GrnaOtdThb7UxITzOcPtg/kbY/v3667MHdxzG45GMnzka/wDj76vCN5xNOnheCPEScC2wDLgOCAcagL3An4D3pZRnq+U4HnBJKQ+ddG0PsOA0458GHgQ6TvNzAIQQy4HlAElJxpcHNJMTaV3m6Hi7e6LZaZEJ+JlQWaqqqZXKphaC/f0YExtluH2320P+Qf1LI/amdQ3z+rMXGmfNo5VSOoUQa4BcoKQfTrUvgoHmU641oUXGX0EIcT1glVK+LYS49CxrW4WWYkZWVtZ5nf4xULz6bFLgRFPseyt2mdVRYVdxj6NPSsBiMX5nvqi4lvb2buLjw4iO1u/oszqocGHSrwMLUkophNhHH46xn7QCoadcC0WTIXrpkRh+B1x5jnaGBR7p7i2NaPpBBZPyZ3OKtGgyM8Uc+75I63K73eRv12qiTJqrMg4uJAZyMmwXmgRw8BzsHAJsQohxJ3XNTQfyThk3DkgBNvfkD/oBYUKISmC2lLLoHGxfcFR1ltDl6SDCHkuo3fjz/d1ud690kBl9IqI1sjlmTlGPfZMcrS8OKhTuK6GjtZP41FiiEiJ0m1dhPgNxtBuAdUKIl4FSTjp6K6V86UxPlFK2CSHeAh4TQnwPmI6m+849Zeh+IPGkx3OBZ4AMVAZCL8Ve2cCkaHZ/QwVdbhdjQ6OJcgT1XjfqoEJLZxeHKmuwWS1MTTxjUozPyPPhQYXJKpo1jGuvvTZ169atIR0dHdbo6GjnvffeW3n//ffr3qhxII52HlDI1zewJHBGR9vDD3rGVQN1wL9JKfOEEPOBD6WUwT36b29VHiFEPeCRUupSqedCoaT3oII5L8idNSUAZEUnnmWkb9hdXI6UMGVkHA77QMt1DJ66ulYqKpsICPAjNTVGt3l79VnlaA3joYceqpg8eXJRQECA3LVrl2PJkiVpM2fObJ8/f76u3RvOem5TCBEohHgCTWfdBFwhpVx40tei/hiSUtZLKa+TUgZJKZO8hxWklJullH0W3JRSblCHFb6O9+itWfrszppSAC6K+Wqmx/Lly1m+fLnP7XtlgwyzZIMDWjQ7aeIIXTveeksjqo2wr+OrLrhZWVmdAQEBEsBisUghhCwoKNA9r64/4cCzQBbwIVoubCTwI70Xougfzc56GpzV+FsCTel46/Z4yK7VHO3MmK9GtKtXaynPq3xcxctsfXbvvh7ZYLJ+9mvK6kzreNsXi5b81hAd6LNPfmZ6F9w77rgjac2aNVGdnZ2WiRMntt90001NA/09zkZ/HO0VQIaUskII8TRaVKscrUkUt+UDWn8wiwkdbw811dDi7GJEYCgjgsIMt9/ldLGvVCskMz3ZnMpWe/dq0kn6NP3e6PZuPADAlPkThn3H277wZRfcv/71ryUvv/xyyfr164PXr18f4nA4dE8V7Y+jDZJSVgBIKUuFEMa/uhS9HGvT6q6nBk8yxf7OWs3JnCobGMX+41U43W7GxUURHmh8IZmWlk6OHqvGbrcySceOt3s3agk40y4506l04+hvpGkERnTBtdlsLF26tPUvf/lL5O9///uYX/ziF9V6zQ39c7Q2IcRCTpRJPPUxUsrP9FyU4vQUtmovyNSgKabY9+qzWTHmbITlmqzP7ttfipQwIS1B10IyezdpEe20Bea8gQ5ljOyC63a7xdGjR03RaKv5alZB3SmPJaCqXxhAq7ORmq4y7MKfkQFjDLcvpezNOJhpUkSb23NQwSxHu2dPj2yQrt/vX1fRQNmhChxB/ozLUJVJT8VXXXCPHz9uW7t2bcgtt9zSFBQU5Hn33XdD33vvvcgXXnjh2NmfPTD6cwQ3RW+jinOjsE2LepKDJmCz6NsOuT8UtzZQ09lGpH8gY0JMqC/g8bCr2NwTYXv2ahG9T/TZiydgMyFdbahzchdcj8fD3XffXe3tgvvWW28Vneu8QghWr14d+5Of/CRZSilGjBjR9etf/7r09ttvN2UzTDFEKOzVZ83R8bK9skF0Yp+V/zMyMnxqP7+8mpbOLhIjw0gIN761emtbJ0eOVmG1Wi5ofXao4asuuCNGjHDt3LnTkObNytGeR5zQZ815QW6v0erAn5rW5cXbysZX7Djak787xhx9eP/+43g8kkmTRhAQ4KfbvEqfPTMXQhdc4xtNKc6JVlcj1V2l2IUfIwPGGm5fSsm2qiIA5sSmGG4fTjjaWaPNcbQn0rr0s99Q3URJ/nEcgf6q0HcfXChdcFVEe55Q2KpFPUkm6bOFLfVUdrQQ6R9IWnis4fa7Xe7ejAOzIlqvPjvdB/rspLnjsfsZf1+HOt4uuGavY7CoiPY8obDN3LSurVWFgBbNWk7TmVUI4bOurfvKKulwuhgTG0lMSNDZn6Az7e1dFByqwGIRup4IU/rs8EA52vOEwlZtI2y0SRthW6uLAJgbl2KK/V7ZYIxJByXyNH12/Ph4AgP1S7P0RrRKn72wUY72PKDVaa4+6/Z42F6tbYSZ52g1fXSWSbJBdo4W0WfOSNFtztryeorySnEE+TNhlvH3VWEcytGeBxxp3QNASvBkU/TZA41VNHV3MioojKRg4wtSd3Q72V1SgRCQlWpOwRWvo83K1O9AQe4newGYvnCK0mcvcJSjPQ843LIbgHHB002xv7Un22CuSdkGu4rLcbk9TEyINaW+QW1tC0VFtTgcdibpWOg75xPtDTRzSbpucyqGJsrRDnE80sPh1h5HG2KOo93Wo8/OiTPneGivPjvWHNkgJ7cIgOnTk7Db9ams5fF4eiPazMun6TKnYuiiHO0Qp7KzmDZXE6H2SGL8jf/Y3OV29dY3mBuXbLh9gG1HevRZk/JnfSEbHNtTTGNNMzGJUYwab065R4VxqDzaIc6Rll0AjAue4bPUqTOxu+44nW4X48NiiHb0WRSpl+eff153+/Wt7Rwor8LPZiXTBH3W45G9EW2WjgVfsj8+IRuYcV8VxqIc7RCnV581STbYVKkdL5/XD9nAF21svjhcjJSQlTKSQBM2jI4eq6KxsZ3YmFASE/XrOOzVZ7MuV/rscEBJB0OYbk8nxe0HEQjGBJuj420o1wrWL0wwJ/1oy6EiAC5OM0cfzs7W7GdlpugWeXa0dZK35SBCCGYsnqrLnIpzZ9WqVRGjR4+eHBAQMCMxMXHKunXrzvzR7RxQEe0Q5ljrftzSxaiAcQTajK9WVdHezMGmagJt9n4V+vb2CtMrsvV4JF8c1vJ3549P0WXOgdKrz2bp5+j3bcrH2e0ibeYYQqOMv6+KE7z99tuhjzzyyKhXX3312KWXXtpWUlLik49NKqIdwhzpkQ3Ghpjz8XJTxVFAS+vyt579PXnFihWsWLFCN/t5x6toaOtgRHgoqTEm5O92dLM/rwwhYIaOBxWyP9Luq0rr6j++6oL72GOPjfjpT39asXjx4jar1UpqaqozNTXVOfgVfxUV0Q5hDnk3wkJmmGJ/Q4+jXWCSbLC5RzaYn6bfx/aBkJ1TiNPpZuKEEYSFBugyp5SSHWu1cpIzrzBHd+8Pkx/4b0O64Ob95semdcF1uVzs378/sKamxpaUlDSlq6vLsnTp0sbnnnuuNDg4WNcGjcrRDlFqOo9T111BgDWYxMDxhtvvdrv5oqeQzKUJxrfNAdiQr23EXWySbLBtu/aanDd3nG5zluSXUX60itCoECbOMf6+nq/4ogtuWVmZ3eVyiffeey9i8+bNBX5+fvIb3/jG2AceeGDE008/fVyvtYNytEOW/OYvAUgLycRqQlvxnTUltLm6GR8WY0pb8aqmVvKOV+Gw25gz1vhCMm63h+07tIh+zmz9Ivpt72sB3KyrMoZ0W/H+RppG4KsuuEFBQR6AFStWVCcnJzsB7r333sonn3xyBKAc7XDgYPNOACaGzjTF/sfHtQ4fl400J+r6PL9HHx6XTIAJaV0HCypobGwnIT6MlJRo3ebd9n42AHOuNue+no/4qgtuTEyMOy4uznmyLOUriUo52iFIq6uRkvYCrMJmij7rkZJPehztEpMc7WcHNEe7aKI5ssXWbdrrcc6ccbq9+Bqqm8jfdgi7n40sdey23/iqCy7At771rdo//elPsddff32Tn5+f/N///d+4yy+/vFG/1WuorIMhyMHmHCSSMcHT8LfqswkzEPbVl1PV0UpCYChTIxIMt9/a2cWOY6VYhGDBBHPyZ7du1V6vc3WUDXaszUVKyfRFUwgINv6+nq+c3AV39OjRU0NCQjzeLriDnfvJJ5+smD59etuECROmTJw4ccrUqVPbf/Ob31Tose6TURHtECS/aQcAE0yTDQ4BWjQ7kGhOSn02arccKsLl9pCRMoLI4EBd5hwIRcW1FJfUERLiYJqO/cE2r9kGwJxrlGwwEHzVBRfA399f/vWvfy0BSvSY73QYFtEKISKFEG8LIdqEEMVCiNtOM+6nQoj9QogWIUShEOKnRq1xKNDpbuNw624Ewjx9tswrG6SZY3+/Fk0unmROWtnGTQcBuHjeeGw2fTasWhpayf1kLxaL4OIbZuky53DhQuiCa2RE+yzQDcQB04G1Qog9Usq8U8YJ4NvAXmAM8LEQolRK+XcD12oa+c1f4pYuUoOmEGrX72x9fylorOZYSx3hfgFcFGP8bn9bVzcbD2ppZUunmqMPex3tgksm6Dbn1nd34nK6mbF4KhGxxmdxnK9cKF1wDYlohRBBwI3AQ1LKVinlFuA94M5Tx0opfyelzJVSuqSUBcC7wDwj1jkU2Nv4BQBTw835ld8v0d73rhg1AZtlYH8emZmZZGYOLs99Q/4xOp0upicnkBBu/PHUouJaiopqCQlxkDFDv7KQG/+hyQaXfHOObnMOB7xdcP39/XU9QGA0RkW04wGXlPLQSdf2AAvO9CShCYTzgT7r7wkhlgPLAZKSzGnapyftrhaOtOxBYGFymPEfL6WUfFCiNQu8OnngTSBzcwffFXrdPu1PZNlUc2QL38sGF+kyp+L8wiiNNhhoPuVaE3C2kOURtDX+ua8fSilXSSmzpJRZMTExg16k2Rxo/hIPbkYHTyHYFm64/d115ZS2NRIXEMzMaOOLbDd3dLK5oAghYOlU/U5j9RcpJZ99rr3R6CkbbF6zA7fLTfrCKYTHKNlgOGJURNsKhJ5yLRRo6WMsAEKIe9C02vlSyvNan+kvuxs2AjA1zBzZ4L0SraX5VYmTsQ5QNtCDT/OO4HS7mTl6FDGhuleqOysHD1ZQWlpPREQQmRkpus37ySsbAFh023zd5lScXxj1ajoE2IQQJ4cp6cCpG2EACCG+C/wMWCylLDNgfaZT31VJYdt+7MKPqeFzDbfv9LhZW5IPwFXJkwy3D/B2jvbncM2MiabY/+iTfQBctmgSVqs+L43yo5Xs33IQR6A/l3xzti5zKs4/DHG0Uso24C3gMSFEkBBiHnAt8OqpY4UQtwNPAEuklMdO/fmFyq6GDQBMCpuNwxpkuP3Py49Q19XG2NBoUw4pFNc2kFtUToDdZkq2QXe3i883aG80l1+uXzHuT17RPqVcfOMsAkPUIYXhipGfD38ABADVwN+Af5NS5gkh5gshWk8a92sgCtgphGjt+fqTges0HI/0kNvwOQCZEYtMWcObx7QaqTePnm5KScJ3cjRt9PKp4wny9zPc/rbtR2hp6WTsmFjGjI7VZU6Px8Onr2qOdsm3L9VlTsX5iWF5tFLKeuC6Pq5vRtss8z4258yliRS25dHorCHcHk1q8BTD7Ve2N7Ox8ih2i4Xrks/d/t13331Oz3N7PLybqznaG7IGnu2gBx+u01p/6xnN7tmQR2VRDTGJUUxfaM7vpRgaqCO4Q4Av6z4CYEbEIizC+E2oNUV78UjJFaMmEOU4d9nC28pmoGwuKKKquZXEyDAyU0aes/1z5fjxBr7ceQw/PxtLLtPvje6957T7uuxfF2MxYXNRMXRQd99kmrprOdC0HQsWZkYtMdy+y+Ph70e1Tg43pZrTWuW1bZpsccusaabIFu+9r+X/Lrx0om6dFGrK6tj67k6sNitX3n2ZLnMqzl+UozWZL+s/xoOHSWGzCbNHGW7/k+MFlLc3kxIcycXxowc1V05ODjk5A6sXfay6nq2Hi3HYbVyfZbxs0tnpZN3HWrbBdddk6Dbv2uc/weP2MP/GWUQlGN/vTNE/cnNzHbNnzx4fEhIyPSkpacorr7zikwR25WhNxOVxsrP+EwDmRF9pyhr+fEjr5PCd8TOxDDKazMrKIisra0DPeb0nmr1q+gTCAx2Dsn8ufPb5AVpaOklLSyAtTZ9si+4uJ/984VMArvnBFbrMqdAfp9PJDTfcMPaKK65oamho2P3MM88Ur1ixInXv3r3+ettSjtZE9jRuos3VRLwjmeRA43NH99aVk1NbRojdnxtSjC9E3dTRyTs9m2C3zzWhwLlH8sY/tJKU11+rXy/CT1/ZSENVE6PTk5lysX4nzIYzvuiCu3v3bkd1dbX94YcfrrLZbFxzzTUtGRkZrS+++KLuHy3VZphJeKSbTdVvA3BxzLWmaJMvFGhO5pbRMwiyG59S9drW3XR0O5kzNonx8fq1i+kvW744RGlpPXFxoSxaqM8bndvt5o3fvwvALf95nSn3VQ/GvPG4IV1wj97yc9O64PaFlJIDBw7onvCsHK1J5DVtp7a7nAh7LNPCjT+aebiphn+WHsDPYuWu8cbXvW3r6ubVL7RNqBULzSmg87c3tgNw8zdn6VZAZvP/baf8SCUJo+NYcJOq1KUXvuiCO23atK7IyEjXww8/HPfQQw9Vr127NmTnzp0hs2bNOm1pgHNFOVoTkFKyoXoNAPNjrzely+0zB7Yg0TINRgSeWobC9/x9+x6aO7rISBlBVqrxKV05uUUUFFQQHh7Isiv0kU08Hg9/+632KeWm/7gGq07O2wz6G2kaga+64Pr7+8s1a9Ycueeee5KeeeaZhClTprRdeeWVDf7+/p7Br/qrKEdrAnlN26jsLCLEFkFGxELD7R9prmVtyQHsFgvfn2R8AZvWzi5e3qy9jlcsnGX4x2uPR/LCi9qJrW/eOBOHQ58uuxvf3MaxPcVEjYhg6Xcu1WVOhe+64ALMmjWrY+fOnQXexzNmzJhw22231er9OyhHazAuj5OPK18DYGHczdgtxmujf9i7oSeanW5KNPvipmzq2zqYnpTAvHH6FdfuL5s2H+TQ4UqiIoO54bqBZUmcDme3k5cf+hsA337kFvwcxt/XCxVfdsHdsWNHwNSpUzvdbrf4/e9/H1NdXW2/55576vRbvYbKOjCYnfWfUNddQbT/SLIiFxtuf1tVER8fLyDQZueeyRfrOnd2djbZ2dlnHFPV1MorWzRt9j+uvMTwaNblcvPSnzcB8O075+kWzf5z9XrKj1aROGGkimZ1xpddcF966aWo+Pj49Li4uPQNGzaEfvTRR4cCAgJ07+agIloDaXe18FnVmwAsjb8DqzD2f7/b4+GJ3Vp+54oJc4kL0LdVTH/a2PzPx1/Q6XSxZMpYZiSP0NV+f1jzdjZlxxsYNTJCN222qbaZVx7R7ut3H7/1vNZmhyK+7IL7/PPPlz3//PM+L8WqIloDWVfxCu3uZlKDJjMx1PiWJn8/tosDjVUkBIbyvTTjd/q/PFbKu7kHsFut/HipvtF0f6iubuYvr2wB4J4fLtEt0+CFlX+lua6FGYunMu861apGby6ELrjK0RpEYWseOQ3rsQob145cYfhH5vK2Jp7c8xkAP59+GQ6bPh+ZT2b58uUsX768z591OV08+vZ6bdzCi0iONv5Y6jPPfUpnp5NL5qdx0czBHTf2sm9zPuv+/Dl2Pxv3Pvu98zZvdqhyoXTBVdKBAXR7Onmn7I8ALIi5gRjHKEPtSyn5efY/aXN1s3RUGssSfXMKbfXq1UDfVbz++Nl2imobGB0byfcW6LMBNRA+XZ/Hli8OERDgxw9/oE+Rl47WDv7rX58D4JaV1zFqvPFSyIWOtwuu2esYLCqiNYC15S9R211OrH8il8TeYLj9147ksKnyGGF+Dh7NMP7s/fajJbywcScWIXj0+svwsxn7/l5Z2cj//O/HAPzw3xYTE62PNv3cv79M+ZFKRk9L5tYHrtdlTsWFiXK0PmZf41ay6z/FJuzcknS/4elc++oreLxnA+yxzGXEBBjb9LC+tZ2fvbEOKWHFwovIMLjerAhmOK0AABV5SURBVNPp5onffkBbexcXzxuv2wbYxje3su6lz/Bz2HngtftUOpfijChH60MqOgp5q+wZAJYlfIf4AGNzRuu72rl361t0e9zcPiaDq5KMbbrY7XJz/+trqWlpIzNlJN9fZGxzQiklTz/7CfvzyoiODuEn9y/TRUM9sruQ//quJhnc/bs7SZlsfGt2xfmFcrQ+osXZwKtFv6Hb00l6+CXMijL2I3uny8mKzf+gpK2RyRHx/HyGsUXFpZT86t317CwsIyYkiN99axk2nTrL9pd33s3lg7W78fOz8dgjN+hS1Lu+soGHr32SzvYulty1gGt/qMogKs6OcrQ+oN3VwsuFv6LJWUti4HiuH/UDQ3ejnR439+94l9y6MhICQ1k9/2b8rcbpolJK/vfjrbyVnYfDbuOZb19DfJi+Obtn49P1eTzznFbr9z/uX8YEHWrNNte1sPLyX1FTWsfE2eP49z8uV1kGin6hsg50psPdxp8LH6Oys4hovxHckfwzQ3VZp8fNv297h4/KCgix+/Pi/Ft0P5hwOjIytA4Fz63fzqoNX2K1CH53yzKmjIo3xL6XDRsP8tvffYCUcPe/XspliwffGLG5voUHlj1O0f5SEieM5NF3VipdVtFvlKPVkcbuWv5S+Cuqu0qJ9Ivnu2MeJdjuk84YfdLm7Oa+bW/zecURQuz+vLzgVtLC9Wmd3R++3LmTJz/YyHPrt2MRgidvWcbiyWMNsw/w7vu5/O/THyMl3HH7XG791uB14eqSGh5Y9jgl+cdJGB3H7z55iIjYMB1WqxguKEerEyVtBbxe/HtaXPXE+I/iO6kPGdoDrLytieVb3iS/sZowPwcvL7iVaZHG5XW2dHbxwJvr+Dz/GHarld/cvJRl09IMs+9yuVn1wgb+b81OAP7lO/O547a5g543b2sBv7r5D9SVN5AyOZHH//kg0SON7+2mOL9RjnaQeKSHrbXv81HFX/HgJiVoEnck/4wAm3FpVB+W5vPgzn/S7OwkNSSS1fNvITUk0jD7+0or+Y+/raWsoZnQAH+evvMaslKNO5RRWdnIE7/9gP15ZVitFn5831KuXDa4jr4ej4c1/72WFx94DbfLzbQFk3jkrZ8SEmFsepziwkA52kFQ3VnK22XPUdKulbOcG30VS+PvxGbR/3hrn/Y7Wnh896d8UKL13VqYMJY/zL6GMD/dO3H0SUe3k2c/3cYrX+Ti9kgO/PZ+ALIe/oEh9t1uD++9v4vVL26gs9NJdHQID//iWqZMHpyTLz5QylPLn+fAVu2+3vjjq/jeb2/HZlcvlwuNJ554Iub111+PPnToUMDVV19dv2bNmiLvzzo6OsRdd92VtGXLltCmpiZbYmJi12OPPVZ28803Nw/UjvrLOQeaumv5rPof5Navx4OHYFs41478PpPCjCko0tLdycuHd/JCwQ5anV04rDZWpi/izrFZhuyCu9we3tt1gGc+2UZVcysWIbjr4gxW+tyyhpSSLV8c4sWXNlFSqpUOXXDJBO770eWEhwee87w1ZXW8+ug/+OjPn+HxSCITIrjvj3cz9xrjW/0ojGHkyJHOlStXVqxbty60s7PzK1lYTqdTjBo1yvn5558XjB07tvvNN98M+5d/+Zcx6enpeWlpad0DsaMc7QCo6ixmR91H5NSvxyWdCCzMjLycpfF3GCIVVLY388ax3bx8aCfNTq2Y0aIRY3l4xlISg32/6dbS2cWanft5fdtujjdob+qTRsTy8HWLmZoY73NH293tYsPGg7z9TjYFh7QqeQnxYaxYvohL5p+7Hly4v4T3nl3Hx3/ZQHenE4vVwlUrLuO7T9ympIIhwsqVK+NLS0v9Xn/99RLQis0kJCSkNzc37zq1QeNAuOuuuxoBdu7cGXj8+PGvpJGEhoZ6nnrqqXLv41tvvbXpl7/8Zde2bdsClaPVmVZXI/lNO9nVsIHi9vze61PD5rI47lZiHL49Utrq7GJz5THeKtrLhoqjeKT2NzUrNpl7J89ndqxvT5s53W62Hynlw70FfLL/MO3dTgCSosK557I5LJuWhsXiuyja45HkHyxn46aDfLo+j8ZGrU1UeHggd94+j6u+MR27feDlDhuqGvninZ189rfN7Nt04r4uuHkO33nsW8O+QMzP995gSBfcx6e9NaS64J6J0tJSW3FxsSM9PX3AJRuVoz2FLncHpe2HKG7L51jbforb8pFo99LP4mB6xAJmR11BnMM3Dq7L7eJAQxW5dWVsqTzG9upiuj1uAOwWC0tHpnHnuCxm+cjBejyS4roGdh4rY9uRErYfLfn/7Z15kBzVfcc/35lZacUuu0bRgQToAoSEhI5CHDGnMYYyiQ1EieVQ4DgHVnBcKbATg1NxYpNUEieRSFLlHEAMBSlwbFBsKBIcF5j7CtgSNzrQgS1ptRLa2Wt2do5f/uie3dFodrXH9Owxv8/Wq+n3+vV7v19PzXdfv+5+P9pT/SvUnbvoFD534WouOWNRJAJrZrS2drDl9T1s3rKH/3t1JwcP9gclPe3UWVz96bP5+GVnDis6Qm9PL++8vI0tP3mLnz35Bm+/8B75fPC9Tmus5/IbLuHTX7zSX6cdp0QRBXc4pNNprVu3btHatWsPrV69evwKraTpwL8DVwAHga+Z2QNl6gn4G+D3wqK7gdvMrKLhJTL5XpKZgxxK7+NAzwe0pPewP7Wblp7d5OkPghlXglMbV7Cs+XzOar6AqfHK3GjK5vMcSHXwfschtrcfYnt7K++1tfLm4X19wgog4OwZJ3PlyUu4ev5yZtQ3VKT/XD5PS3snHxxK8sGHbexo+ZC397bwzt5WutJHXhUtmjWdq1acwSdXnMGCmZVZR9bMOHy4i337k+zb18b7O1vZsaOF7dtbONx25O9n1qwmLr7oDC69eAlLl84ddB66pztNy64D7N95gJ9v3cf2zTvZsXkXe975Bbls/3mtm5LgnE+s4MJfO5+L1p5HQ9PI53YnI0MdaVaDqKLgDpVcLsfatWsX1tXV5e+55549I2mjmiPabwO9wGxgFfCYpC1m9lZJvS8A1wArAQN+DOwE/nU0nf/s8FO8lXyJZOYQyd5WunLlbxzGiHHStFNZ0HAm8xuWcmrjCurjR/8IzYys5enJZUnnMqRzOXpyGXpyWdp7e2jP9JDsLaQU7b09HEp3s7+7nf2pDg70dPZNA5RyetMMVs84mTUzTuGSOacOKK5mRiaXJ5XJ0J3O0N3bG372bydTPRzuSnGos5sPu7r5sDPFwY4u9rZ1kMnlyrY7u6mRFfNO5KOnzef80+ZxyvTmAcUtm83R05MhleoX5y2v7yGVytDT00uyPUWyrZu2ZDfJZIpksptDhzrZ35Iknc6WbbOhYSpnLpnD0sUnsnjRTGaf0EC6O03P3kO8sPUXdLZ1kWxtJ3mwI/xsp+1AkpbdrRxuSZZtUxILz5rHykuXsepjy1n1sWU0NFfmn5YTLVFGwT0W+XyedevWLWhtbU088cQT26ZOnTqiAV9VhFZSA7AWWG5mncBzkh4BbgBuK6n+W8AGM/t5eOwG4EZGKbT//MwjzFq4qy+fz4tUbx1d6akkU9NoT00j2VNPsvs4shYH7cPYB3oiGFYqUH0EFidYJWI0V84GykAsJeLdMWIpEUuJREeMA7lOHte7PK53+Qv9GBOYhBXsCBOjvHSP9xqJdJ66tFHXk2dqd54pnXliuRRbaeU93uA7galhUsH0vkSRAJ+xJFiT9ZavHHWhUhblcsTSGWLpDOrugfYu8skuUl09/PQRGMlqz4m6OLPmz+TEhbOYu2g2i1Yu4LTVC1mw/BSmNdSPoEVnrIkyCm4mkyGTySiXyymXy6m7u1t1dXVWVxdMS11//fXztm3bVv/MM89sbWxsHPFVdbVGtIuBrJltLSrbAlxSpu6ycF9xvbIvq0v6AsEImHnz5g1qwLb9x7OtbiHdmSl0ZaeQztb1CccRHAehhBwbA/KFT6HCdk4oJ5Tt3yYHygplhHoFWSE7sv880DtcLcgbyhOk3NHbsSzEskYsY0Xbgcgqf2RThkgTH94/EDPI5SGbY+70FZDNQ1tH8JnPo0wWerMQfiqTCfKpNMrl+/wuRsDUaVOY1lhPfcNU6hvqqe/bnkpD83F8ZEYTzTObaZ5xPM0zm2ie2cSseTP4pbknEI97cMTJRHEU3Hw+z4033nigEAV306ZNu0bT9q233jr3jjvu6FtxqKGhYfott9yyb+PGjXu3bt065cEHH5w5ZcoUmzNnTt8bMBs2bNh90003HTWaHoxqCW0jUHqtngTKrXbSGO4rrtcoSaXztGZ2J3AnwJo1awZVx2umX8Cr7+xBCFkwGJSJmIK8+v4gZkKCGAq2CesiYkCCOPEwLxXrksJBZlBS2KdCPgFKKBDzcCRYeHCvcGmucFdcQfuJWIyYRFwx4tIRKYb6bkhJ/W322VM02iy9cSXpiHqFqpKIx0U8FiMWE4l4jFg8Rjwm4vFYkMLtRDx2ZLtF/ZVONSTq4sTr4sQTcRJ1cRJ1CeJ1wXahLF6X6Mv7qlhOgSij4G7cuHFv8SNcxSxevLjXzCoyV10toe0EmkrKmoCOIdRtAjpHezPs5t+4YjSHO0OgECtsoACNjjMSPAru0NkKJCSdXlS2Eii9EUZYtnII9Zxxxvr161m/fv1Ym+FMIiZLFNyqCK2ZdQGbgNslNUi6ALgauL9M9fuAL0s6SdJc4CvAvdWw03Gc8UUhCu5I7/aPF6oZYeGLwDTgAPAgcJOZvSXpIkmdRfX+DXgUeAN4E3gsLHMcx5mQVO05WjP7kOD52NLyZwlugBXyBnw1TI7jOBMejxnmOI4TMS60juMA5PP5vD9TNwrC81f2dUsXWsdxAN5sbW1tdrEdPmZGOp2u271790eA58rVUYXXahkzJCWBwut3zRz50kMhP4NgQZtKUNrHaOoOtL9c+UC+DbY9Hv2ulM+l+bH2u1rfdWl+uH7PN7OZhcxrr702K5FI3A0sxwdgwyUvKZnL5e7N5/P/cvbZZx+9Vq2ZTYoE3FluuzgPvBpFf6OtO9D+cuUD+XaM7XHnd6V8Hm9+V+u7rpbfniqTJtN/rkcH2C6Xr3R/o6070P5y5YP5Ntg5qBSV8rtSPpfmx9rvan3Xpfmo/HYqwKSZOhgKkl41szVjbUe1cb9ri1r1ezwzmUa0Q+HOsTZgjHC/a4ta9XvcUlMjWsdxnLGg1ka0juM4VceF1nEcJ2JcaEuQNFvSC5KelvSkpDnHPmriI+lcSS9KekbSg5KGHmJ2giKpWdIrkjolLR9re6JG0rckPSvp/lr4fscTLrRHcxC40MwuIViy8XfH2J5q8QFwmZldDOwiWMZystMN/Arw0FgbEjWSVgInmdlFwLvAr4+xSTWFC20JZpYzs0IYq+OpkUXHzWyfmRXCOfdydCivSYeZZcysdaztqBIfBf433H4cuGAMbak5JrTQSvqSpFclpSXdW7JvuqT/ktQlabek64bR7ipJLwNfYmTBWCMlKr/D4+cDVzDOHoCP0ueJxCjOwwn0x+1LAtOrZLJDFdejjYi9wF8CVxIsKl7MtwlGZrOBVcBjkrZYsNj4icB3y7T3WTPbb2abgfMkfQb4GvD7kXkwMiLxW1ITQdSLz5tZJjrzR0QkPkdpcESM6DwAbfTH4msGhhXF1RkdE1pozWwTgKQ1wMmFckkNwFpguZl1As9JegS4Abgt/IFdWq5NSVPMrLAoRJJgHm9cEZHfCQJB+qaZvRetB8MnCp8nIiM9D8ALwJcJ7jtcCTxfZdNrmgk9dTAIi4GsmW0tKtsCLBvCsavCO+8/AW4G/i4KAyNiNH7/JnAe8HVJT0laF4WBETAan5H03wRTJXdJ+nzlzasag56H8CqtRdKzYdnD1TexdpnQI9pBaKR/PqpAkuDm1qCY2SvAxVEYVQVG4/f9lA+WOd4Zsc8AZnZVxS0aG455Hszsj6tqkdPHZB3RdtI/H1WgCegYA1uqSS36XYs+l8PPwzhmsgrtViAh6fSispVM/ke1atHvWvS5HH4exjETWmglJSTVA3EgLqleUsLMuoBNwO2SGiRdQPAA/kS8ND6KWvS7Fn0uh5+HCcpYrzw+mgR8A7CS9I1w33TgB0AXsAe4bqztdb/dZz8PtZl8mUTHcZyImdBTB47jOBMBF1rHcZyIcaF1HMeJGBdax3GciHGhdRzHiRgXWsdxnIhxoXUcx4kYF1rHcZyIcaF1BkTSLkmXj+L4tyRdWkGTHGdC4kI7yQnFMhVGet0v6V5JjRH1c4Qom9kyM3uqwv2cIMnCkDso4Paw/1WV7MtxKoULbW3wKTNrJAhvspogPM9EZRVw2Mx2h1EFHgIuA861YHFrxxl3uNDWEBaEdfkRgVgBIGmupIcltUraKekPBzpe0m2SdkjqkPS2pGvD8vuBecCj4cj5q2H5LkmXS7pV0kMlbf2jpH8arg2h7ZslzQOeI1jc+jIzOzCik+I4VcCFtoaQdDLwSWB7mI8RRLvdApwEfBy4WdKVAzSxA7iIILjfN4H/kDTHzG4gWC3qU2bWaGZ/W3Lcd4GrJB0f9hsHPgM8MAIbVgN1wEvA/Wb2O9Yf481xxiUutLXBDyR1AB8AB4A/D8vPAWaa2e1m1mtm7wN3AZ8t14iZfd/M9ppZ3sz+E9gGnHuszs1sN0HY9mvDosuAbjN7abg2EIxolwGvm9nG0p2Svifp+TDu2Y8kLT2WfY4TNS60tcE1ZnY8QTTYJcCMsHw+MFdSWyEBf0IQrvooJH1O0uaiusuL2joWDxAEgAS4LswPywZJU4GlBNFel0q6uUw/pwEXm9mlYTvfGaJ9jhMZkzU4o1MGM3ta0r3A3wPXEIxwd5rZ6YMeCIR3+e8iuLR/0cxykjYDKjR/jCa+D2wIpy+uBX45LB+yDQTCngOeDdt4WtLrZvZkaOMUIGdmudDf1yTNCCMQZIfQvuNEgo9oa49/AD4haSXwCtAR3qyaJikuabmkc8oc10Agpq0Akn6bQPgKtACLBurUzFqBp4B7CIT1nXDXcGxYDbxpZlkz+ynwB8D3JC0I9y8hiJ1VzHEuss5Y40JbY4SCdx/wZ+HI71cJ5j13AgeBuwludpUe9zawAXiRQFTPAp4vqvLXwJ+Gl/9/NED3DwCX0z9twHBsCOtsLjr2vrCtH4aPei2nKBihpLPw4ITOOMBD2TiTBkl/BbxsZj+UNAt4GPh6pV+acJzh4kLrTBokPULwiFg3wVzut8zsf8bWKsdxoXUcx4kcn6N1HMeJGBdax3GciHGhdRzHiRgXWsdxnIhxoXUcx4kYF1rHcZyIcaF1HMeJGBdax3GciHGhdRzHiZj/B/b2yg5qQLULAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "relative_kd = np.logspace(-3, 0, 200)\n",
    "mu_values = [1, 3, 6, 9, 12]\n",
    "mu_colors = plot_utils.set_color(np.arange(len(mu_values)) / len(mu_values))\n",
    "\n",
    "width, height = mpl.rcParams[\"figure.figsize\"]\n",
    "fig, ax = plt.subplots(figsize=(width * 1.25, height))\n",
    "for mu, color in zip(mu_values, mu_colors):\n",
    "    occupancy_probability = 1.0 / (1 + np.exp(-2.5 * np.log(relative_kd) - mu))\n",
    "    ax.semilogx(relative_kd, occupancy_probability, color=color, label=rf\"$\\mu = {mu}$\")\n",
    "    \n",
    "ax.set_xlabel(r\"Relative $K_D$\")\n",
    "ax.set_ylabel(\"Pr(bound)\")\n",
    "ax.legend(loc=(1.02, 0))\n",
    "\n",
    "# Show where the 50% cutoff is for mu of 9\n",
    "ax.axhline(0.5, color=\"k\", linestyle=\"--\")\n",
    "ax.axvline(np.exp(9 / -2.5), color=\"k\", linestyle=\"--\")\n",
    "\n",
    "plot_utils.save_fig(fig, os.path.join(figures_dir, \"predictedOccupancyCurveVsMu\"), timestamp=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}