Channel_Capacity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fab8de68",
   "metadata": {},
   "source": [
    "# Importing the libraries required for this code\n",
    "This code is written in Jupyter Notebook (6.4.11).\n",
    "This is an opensource code, can be modified and distributed for any purpose.\n",
    "The following is the library packages that we used in the code:\n",
    "Numpy (1.22.1),\n",
    "Pandas (1.4.0),\n",
    "SciPy (1.7.3),\n",
    "Matplotlib (3.5.1),\n",
    "Seaborn (0.11.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "060ebd37",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import itertools\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from scipy.optimize import minimize\n",
    "from scipy import stats\n",
    "\n",
    "\n",
    "\n",
    "#from decimal import *\n",
    "#from mpl_toolkits.axes_grid1.inset_locator import inset_axes\n",
    "#from matplotlib.collections import PathCollection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbaa3f9a",
   "metadata": {},
   "source": [
    "# Importing data (dectin-2 U937 cells stimulated with FurFurMan)\n",
    "\n",
    "The imported input_output consist of input (Doses, first column) and output (GFP, second column) information (See SI section 1).\n",
    "The total measured input and output is 124488."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "7f27b666",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Data:  [[ 0. 11.]\n",
      " [ 0.  2.]\n",
      " [ 0.  3.]\n",
      " ...\n",
      " [ 7. 10.]\n",
      " [ 7.  5.]\n",
      " [ 7. 17.]]\n",
      "Shape: (124488, 2)\n"
     ]
    }
   ],
   "source": [
    "input_output = np.genfromtxt(\"dectin-2 FFM.csv\", delimiter = \",\")\n",
    "\n",
    "print(\"Data: \", input_output)\n",
    "print(\"Shape:\", input_output.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e65fbe2e",
   "metadata": {},
   "source": [
    "# Plotting joint distribution of the input (doses) and out (GFP) information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "049420a8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entropy of input = 3.0\n",
      "Entropy of output = 1.53\n",
      "Joint Entropy = 4.23\n",
      "Mutual Information = 0.3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcUAAAG4CAYAAAA0ffybAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABVA0lEQVR4nO3dd3xUVfrH8c8zk4ReBUuCK6goiiL6i6hrXV17wbrq2lBWcHXtvSxGXey9rIqKIGvv2EVdd9fdtQRF7IodcBUboLRk5vn9cW9gjEEySWbOTPJ9v17zytw7d+Z+MynPnHPPPdfcHREREYFE6AAiIiKFQkVRREQkpqIoIiISU1EUERGJqSiKiIjESkIHaCm9evXyvn37ho4hIq3Q5MmTv3b33qFzSO61mqLYt29fqqurQ8cQkVbIzD4NnUHyQ92nIiIiMRVFERGRmIqiiIhITEVRREQkpqIoIiISU1EUAfr06YOZ5e3Wp0+f0N+yiDSg1ZySUV+fPn2YMWNG3vZXUVHB9OnT87a/ptL70rAZM2ZQVVWVt/3lc19Npd8VaYtabVEspH9yhfTPRe9L4f/Tzef7ot+VhhXL74q0vFZbFAtJIf1zKSR6XxqWz/dF70nDiuV9kZanY4oiIiIxFUURkRwxsx/qLQ8zs2ub+FpDzGxKfHvdzPZYynb9zOwlM5tmZnebWVlT9tdWqSiKiBSHN4FKdx8M7ADcaGYNHQK7CLjC3VcHvgOG5y9i8VNRFBEJwMweNrOD4/sjzez2X9re3ee5e2282B7wBl7TgK2B++JV44HdWyx0G6CBNiIiudPBzKZkLPcEJsb3RwD/NrOPgROBjQHM7AgAd7+h/ouZ2UbAWGAV4KCMIllnOeD7jPXTgYqW+VbaBhVFEZHcmR93dwLRMUWgEsDdvzSzUcDfgT3c/dt4/c+KYR13fwkYaGZrAePN7Al3X5DD/G2Ouk9FRMJZF/gGKM/mSe7+DvADsE69h74Bumcca+wD5O8Ez1ZARVFEJAAzGwLsCKwPnGRm/Zaxfb+6YmdmqwADgE8yt3F3J2p57h2vOgR4uGWTt24qiiIieWZm7YCbgMPcfSbRMcWxFjmi7rhiPZsBr8fHKB8EjnT3r+PXe9zM6lqbpwInmNk0omOMt+T422lVdExRRCRH3L1zveVxwLh4cb2M9RNZMgCnwWOK7j4BmLCUx3bKuP8RMKSpmds6tRRFRERiKooiIiKxnBZFM9vBzN6Lpxs6rYHHtzCzV82s1sz2rvfYIWb2QXw7JJc5RUREIIdF0cySwHVEo6vWBvY3s7XrbfYZMAy4o95zewJnAxsR9Y2fbWY9cpVVREQEcttSHAJMc/eP3H0RcBcwNHMDd//E3acC6XrP3R6Y5O7fuvt3wCSiuf5ERERyJpdFsQL4PGM5m+mGGvVcMxthZtVmVj1r1qwmBxUREYEiH2jj7mPcvdLdK3v37h06joiIFLlcFsUZwMoZy9lMN9Sc54qIiDRJLoviK0D/eGqiMmA/lpycuixPAduZWY94gM128ToREZGcyVlRjC9d8ieiYvYOcI+7v2Vm55rZbgBmtqGZTQf2Ibpg5lvxc78FziMqrK8A59bNIC8iIpIrOZ3mzd0fBx6vt25Uxv1XiLpGG3ruWKLrhomIiORFUQ+0ERERaUkqiiIiIjEVRRERkZiKooiISExFUUQkR8zsh3rLw8zs2ia+1lgz+8rM3vyFbbYys9lmNiW+jYrXtzezl83sdTN7y8zOaUqGtkAXGRYRKQ7jgGuB25ax3b/cfZd66xYCW7v7D2ZWCrxgZk+4+4s5yFnU1FIUEckzM+tiZh/HBQoz65q53BB3/yfQpPO1PVLXai2Nb96U12rtVBRFRHKnQ0ZX5hTgXAB3nws8D+wcb7cf8ACwnpnd3Mx9bhJ3kz5hZgPrVppZMs7wFdFViF5q5n5aJRVFEZHcme/ug+tuwKiMx24GDo3vHwrc6u7V7v6HZuzvVWAVd18PuAZ4qO4Bd0/FGfoAQ8xsnWbsp9VSURQRCcDd/w30NbOtgKS7L3UATRavOaeumzSeUazUzHrV2+Z74O/oGrUNUlEUEQnnNuAO4NaWeDEzW9HMLL4/hOh//Ddm1tvMusfrOwDbAu+2xD5bGxVFEZFwbgd6AHcCmFnl0o4pmtmdwH+BNc1supkNj9cfYWZHxJvtDbxpZq8DVwP7ubsDKwF/N7OpRBdZmOTuj+byGytWOiVDRCRH3L1zveVxRKdW1NkMuC/u0sTdq4EGjym6+/5LWX9Dxv1riU7bqL/NVGD9rMK3UY1qKZrZ2g2s26qlw4iItBVmdg1wIdFl8qRANLaleI+ZTQAuBtrHXyuBTXIVTESkNXP3o0NnkJ9r7DHFjYCVgf8Q9UfPBDbNVSgREZEQGlsUa4D5QAeiluLH7p7OWSoREZEAGlsUXyEqihsCmwP7m9m9OUslIiISQGOPKQ6PR0UBfAEMNbODcpRJREQkiF8simbWHjgCWN3M/g+4xd1rAdx9Qh7yiYiI5M2yuk/HE40yfQPYEbgs54lEREQCWVb36druvi6Amd0CvJz7SCIiImEsq6VYU3enrttURESktVpWS3E9M5sT3zeia4PNie+7u3fNaToREZE8+sWi6O7JfAUREREJTVfJEBERieW0KJrZDmb2nplNM7PTGni8nZndHT/+kpn1jdf3NbP5ZjYlvt3wsxcXERFpYTm7dJSZJYHriC5mOR14xcwmuvvbGZsNB75z99XNbD/gImDf+LEP3X1wrvKJiIjUl8uW4hBgmrt/5O6LgLuAofW2GUp0LiTAfcA2dVeNFhEpdmb2Q73lYWb2s+sdtsB+usYXHm7wtc3sEjN718ymmtmDZtY9Xj8ko0fudTPbo6WzFZtcFsUK4POM5enxuga3iU/5mA0sFz/Wz8xeM7N/mNnmDe3AzEaYWbWZVc+aNatl04uIFI/zgH/+wuOTgHXcfRDwPnB6vP5NoDLuldsBuNHM2vTF5wt1oM0XwK/cfX3gBOAOM/vZ6R/uPsbdK929snfv3nkPKSLSVGb2sJkdHN8faWa3m9lqZvZqxjb9M5eX8jr/B6wAPL20bdz96YxzzV8E+sTr52Wsbw9407+j1iGXnwhmEF2DsU6feF1D20yPP510A75xdwcWArj7ZDP7EFgDqEZEpHh0MLMpGcs9gYnx/RHAv83sY+BEYGN3/9bMZpvZYHefAhwK3GpmuxG16EZlvriZJYim3zwQ+G0jMx0G3J3xGhsBY4FVgIPa+kQtuWwpvgL0N7N+ZlYG7MeSX4Y6E4FD4vt7A8+5u5tZ73igDma2KtAf+CiHWUVEcmG+uw+uuwGLi5q7fxkv/x040d2/jR+6GTg0/h+4L3CHu0+sXxBjRwKPu/v0xoQxszOBWuD2jBwvuftAoksDnh5fCKLNyllL0d1rzexPwFNAEhjr7m+Z2blAtbtPBG4BJpjZNOBbosIJsAVwrpnVAGngiIxfGBGR1mJd4BugPGPd/cDZwHPAZHf/5heevwmwuZkdCXQGyszsB3dv6BS4YcAuwDZxb9xPuPs78cCgdWjDvXI5PaDq7o8Dj9dbl/lJaQGwTwPPu5/oF0NEpFUysyFEVx9aH/iHmT3t7h+7+wIzewq4nui0taVy9wMyXm8YURdrQwVxB+AUYEt3n5exvh/wedyIWQUYAHzS7G+uiBXqQBsRkVbLzNoBNwGHuftMomOKYzNOSbudqJfs6Xj73eJetmz2cbOZVcaL1wJdgEn1JkTZDHg9Pu75IHCku3/djG+t6LXpobciIrnk7p3rLY8DxsWL62Wsn8hPx1xsBtzq7qmlPN7QvjJfG3f/Q8b91ZfynAmALhifQUVRRKSAmNmDwGrA1qGztEUqiiIiBcTd2/ysMiHpmKKIiEhMRVFERCSmoigiIhJTURQREYmpKIqIiMRUFEVERGIqiiIiIjEVRRERkZiKooiISExFUUREJKaiKCIiElNRFBERiakoioiIxFQURURyxMzczP6WsVxiZrPM7NF4eZiZXduI1znEzD6Ib4csZZvzzGxqfBHhp82svOW+k7ZDRVFEJHd+BNYxsw7x8rbAjGxewMx6AmcDGwFDgLPNrEcDm17i7oPcfTDwKDCqyanbMBVFEZHcehzYOb6/P3Bnls/fHpjk7t+6+3fAJGCH+hu5+5yMxU6ANyFrm6eiKCKSW3cB+5lZe2AQ8FJDG5nZbmZ2bgMPVQCfZyxPj9c19Bqjzexz4ADUUmwSFUURkRxy96lAX6JW4uO/sN1Ed29WIXP3M919ZeB24E/Nea22SkVRRCT3JgKXkn3XKUTHIFfOWO7Dso9L3g7s1YR9tXkqiiIiuTcWOMfd32jCc58CtjOzHvEAm+3idT9hZv0zFocC7zYpaRunoigikmPuPt3dr/6lbZZ2TNHdvwXOA16Jb+fG6zCzm82sMt70QjN708ymEhXOY1v0m2gjSkIHEBFprdy9cwPrngeej++PA8bF9ycSdbM29DpjiVqb9df/IeO+uktbQE5bima2g5m9Z2bTzOy0Bh5vZ2Z3x4+/ZGZ9Mx47PV7/npltn8ucIiIikMOiaGZJ4DpgR2BtYH8zW7veZsOB79x9deAK4KL4uWsD+wEDic7H+Wv8eiIiIjmTy5biEGCau3/k7ouIztUZWm+bocD4+P59wDZmZvH6u9x9obt/DEyLX09ERCRnzD03kx6Y2d7ADnV93mZ2ELCRu/8pY5s3422mx8sfEk1lVAW86O5/i9ffAjzh7vfV28cIYES8uCbwXk6+GRFp61Zx996hQ0juFfVAG3cfA4wJnUNERFqHXHafNuaE08XbmFkJ0A34ppHPFRERaVG5LIqvAP3NrJ+ZlRENnKk/3HgiUHcZlL2B5zzqz51INFdgOzPrB/QHXs5hVhERkdx1n7p7rZn9iWjmhSQw1t3fik9OrY7PybkFmGBm04BviQon8Xb3AG8DtcBR7p7KVVYRERHI4UAbERGRYqNp3kRERGIqiiIiIjEVRRERkZiKooiISExFUUREJKaiKCIiElNRFBERiRX13KeZOnbs6N27dw8dQ0RaoS+++OLrbCcE79Wrl/ft2zdHiaSxJk+enNXPrtUUxe7duzPzo89CxwDgnItGc/apZ4aOASjL0hRallGnFEaWcy8urPelULJYh9JPs31O3759qa6uzkUcyYKZZfWzU/epiEgOzJw5M3QEaQIVRRGRHPjiiy9CR5AmUFEUERGJqSiKiIjEVBRFRERiKooiIiIxFUUREZGYiqKIiEhMRVFERCSmoigiIhJTURQREYmpKIqIiMRUFEVERGIqiiIiIjEVRZFWYPNttgwdQaRVUFEUaQX+/tSzoSOItAoqiiKtQElJq7leuEhQ+ktqoh1225mvv/4ar1vhxPejOw89+HC05NE698WbgXv8dclzFj8er/N43ZLtoucseRzcf7qOjOdkvs5fx99Num77eH36Z8+PH/MGXrvedktCLQn9k/U/26buhZ2qy69vcP32v67kyYmPtcSPpk3acbedeULvn0izqSg20dfffM3LL7xI2h1PO2mHdNpJu0df4/uedlLupFLRulR6yddUOp1xP1pfW+9Wk07/bF1t2qnxuseX3GrT6Yz70TY1aWdRKs2itLMwlWZRylmUjpbr1tdtU5NyalJpaut9rUk5tRmPp1JpSNfGtxSkU3gq9ZPlxfdTNZCqwVM1UFuzeHnx+toaqF3E19+8G/pHWtTumnBH6AgirYK6T0VagW7duoWOINIqqCiKiIjEVBRFRERieSmKZjbWzL4yszfrrT/azN41s7fM7OKM9aeb2TQze8/Mts9HRhGRlrTSSiuFjiBNkK+BNuOAa4Hb6laY2W+AocB67r7QzJaP168N7AcMBMqBZ8xsDXdP5SmriEizlZeXh44gTZCXlqK7/xP4tt7qPwIXuvvCeJuv4vVDgbvcfaG7fwxMA4bkI6eIiLRtIU/JWAPY3MxGAwuAk9z9FaACeDFju+nxup8xsxHACIhG351z0ejcJs7w2eef02f1VfK2v5CS8a3dMrabN28eHTt2bNI+vu/YMa8/PyDv+/sl515cOFkK6X0ppCzSNoQsiiVAT2BjYEPgHjNbNZsXcPcxwBiA8vJyP/vUM1s85NI8MnEi1f9+qcHHzrloNPnM8kvymaVy042W+p7kO8uyFFqWUaecEToGAOdefD6jTimM9+XciwvnZ1RVVRU6guRJyNGn04EHPPIy0SQrvYAZwMoZ2/WJ14mIiORUyKL4EPAbADNbAygDvgYmAvuZWTsz6wf0B14OFVJEpClmzpwZOoI0QV66T83sTmAroJeZTQfOBsYCY+PTNBYBh7i7A2+Z2T3A20AtcFQhjjwdcdgfQkcoOHpPRJb44osvQkeQJshLUXT3/Zfy0IFL2X40UNBH2EcMPzx0hIKj90REip1mtBEREYmpKIqIiMRUFJfhyaefYs1BA1l94AAuvOTinz2+cOFC9j3w96w+cAAbbf5rPvn0k/yHzLNlvSf/fOFfbLDJhpR0bs99D9wfIKGISNOoKP6CVCrFUccdwxMPP8Lbr03lznvv4u133v7JNreMG0uPHt2Z9ta7HH/0sZx6ZmGcb5YrjXlPfrXyyowbcwu/33e/QClFRJpGRfEXvPzKy6y+2mqs2m9VysrK2G+ffXn40Ud+ss3Djz7CIQccBMDee+7Fs88/h9ddkb4Vasx70neVvgxadxCJhH69RKS46L/WL5gxcyYr9+mzeLlPRQUzZsxoYJtoroGSkhK6de3GvHnz8poznxrznoiIFCsVRRERkZiK4i+oKC/n8+nTFy9PnzGDioqKBrb5HIDa2lpmz5nd5Emxi0Fj3hMRkWKlovgLNqzckA+mTePjTz5m0aJF3HXv3ey28y4/2Wa3nXdh/O0TALjvgfvZesvfYGYh4uZFY94TEZFipaL4C0pKSrj2iqvYftedWWvwuvxur30YuPZARp1bxcR4cMnwYYfxzTffsvrAAVx+9ZVc+JeCnoin2RrznrxS/Qp9VuvLvQ/cz8ijj2TgBuuFDS0i0kghLx1VFHbaYUd22mHHn6w7d1TV4vvt27fn3jvuynOqsJb1nmxYuSHTP/wkv6FERFqAWooiIiIxFUUREZGYiqKIiEhMRVFERCSmoigiIhLT6NMil+zUjnXXWYfa2hRrrTmA8TePbdWTB4iI5JJaikWuQ4cOTHlpMm9OnkJZWSk33HRj6EgiIkVLRbEV2XzTzZj20Ye8Uv0KgzZcnwULFrBo0SIGbrAeb771JgcPH8ZDEx9evP0Bww7i4UcmBkws0nqttNJKoSNIE6gothK1tbU88fRTrDtwHTas3JDddt6Vs6pGMWnSJA7c//esM3Adhg87jHETbgNg9uzZ/OfFF9l5x50CJxdpncrLy0NHkCZQUSxy8+fPZ/BG/0flphvxqz4rM3zYYQCMOuMsJj33DDNnzuSUE04CYMvNt+CDD6cxa9Ys7rznLvbafQ9KSnRYWUSkjrWWC+KWl5f7yJEjQ8fIu9GjR3PmmWf+bP3cuXMZO3YsyWSSESNGUFZWBsALL7xAMpnkzTffZOjQoSy//PL5jixSdKqqqia7e2U2z6msrPTq6upcRZJGMrOsfnZ5K4pmNhbYBfjK3dep99iJwKVAb3f/2qLLTFwF7ATMA4a5+6u/9Prl5eU+86PPchM+S+dcNJqzT/15ocqFzr2688PX3/9s/W5778F+e/+OWyaMY63+a3LtlVcD8OWXXzJk81+z4gor8NK//pOXjHXy+b4sS6Fl+fPJZ4SOAcB5l5zPqFMKI8u5F59fMD8j61Cqoliksi2K+ew+HQfsUH+lma0MbAdkVrQdgf7xbQRwfR7ytRq33T6B0pJSfr/f/my22Wa8Mrma557/OwArrLACaw0YwKEHHxI4pbSkLbf9TegIIq1C3oqiu/8T+LaBh64ATgEym6xDgds88iLQ3cw0lKsBDbUSDz7gIO6/6x4AEokEL/3rP2y9VfRPc968eXwwbRr7/26/fMaUHHv2iUmhI0g9U6dOpU+fPqFjSJaCDrQxs6HADHd/vd5DFcDnGcvT43X1nz/CzKrNrHrevHk5TNo6PPPcs6w1eF2O/uORdOvWLXQcaUEaMFV4ampqmDFjRugYkqVgf0lm1hE4g6jrtEncfQwwBqJjii0UrdX67dbb8On7H4aOISJSsEJ+vFwN6Ae8Ho2roQ/wqpkNAWYAK2ds2ydeJyIikjPBuk/d/Q13X97d+7p7X6Iu0g3c/X/AROBgi2wMzHb3L0JlFRGRtiFvRdHM7gT+C6xpZtPNbPgvbP448BEwDbgJODIPEUVEpI3LW/epu++/jMf7Ztx34KhcZxIREcmkad5ERERiKooiIiIxFUWRVuDLL78MHUGkVVBRFGkFRh79x9ARRFoFFUWRVuChex4IHUGkVVBRFBERiakoioiIxFQURUREYiqKIiIiMRVFERGRmIqiiIhITEVRREQkpqIoIiISU1EUERGJqSiKiIjEVBRFpMWkUil22WNo6BgFobS0lIqKitAxJEsqiiLSYpLJJH+7dXzoGAVh0KBBTJ8+PXQMyZJFF7kvfuXl5T5y5MjQMUSkFaqqqprs7pXZPKeystKrq6tzFUkaycyy+tm1qqI486PPQscA4JyLRnP2qWeGjgEoy9Kcc9FoRp1yRugYAJx78fmcdVJhZPnLpefz55MLI8t5l5xfML8v1qFURbFIZVsU1X0qIiISU1EUEcmBqVOn0qdPn9AxJEsqiiIiOVBTU8OMGTNCx5As5aUomtlYM/vKzN7MWHeJmb1rZlPN7EEz657x2OlmNs3M3jOz7fORUUREJF8txXHADvXWTQLWcfdBwPvA6QBmtjawHzAwfs5fzSyZp5wiItKG5aUouvs/gW/rrXva3WvjxReBus73ocBd7r7Q3T8GpgFD8pFTRETatkI5pngY8ER8vwL4POOx6fG6nzGzEWZWbWbV8+bNy3FEERFp7YIXRTM7E6gFbs/2ue4+xt0r3b2yY8eOLR9ORETalJKQOzezYcAuwDa+ZBaBGcDKGZv1ideJiIjkVLCWopntAJwC7ObumX2fE4H9zKydmfUD+gMvh8goIiJtS15aimZ2J7AV0MvMpgNnE402bQdMMjOAF939CHd/y8zuAd4m6lY9yt1T+cgpIiJtW16Korvv38DqW35h+9HA6NwlEmldJtwxgYN+f1DoGCJFL/hAGxFpvh9/1OhrkZagoijSChxxuC6bJtISVBRFRERiKooiIiIxFUUREZFYVqNPzaw3cDjQN/O57n5Yy8YSERHJv2xPyXgY+BfwDKBzB0VEpFXJtih2dPdTc5JERIpeKpXi4OHDuH3chNBRRJok22OKj5rZTjlJIiJFL5lMcv65fwkdQ6TJsi2KxxIVxgVmNje+zclFMBEpTqv8apXQEUSaLKvuU3fvkqsgIiIioWU996mZ7QZsES8+7+6PtmwkEZHiV1payvLLLx86hmQpq+5TM7uQqAv17fh2rJldkItgIiLFbNCgQUyfPj10DMlSti3FnYDB7p4GMLPxwGtEl4ESEREparbkgveN2NhsKrCVu38bL/ck6kIdlKN8jVZeXu4jR2pSZGmbampqKC0tDR2j1aqqqprs7pXZPKeystKrq6tzFUkaycyy+tll21K8AHjNzP4OGNGxxdOyfI2cOfvUM0NHAOCci0YrSwMKLctZJ50ROgYAf7n0fM44oXl/RtvuvB2THnu62VnOv/xCzjyxMDp+Rl92AaNOKYzfl6qqqtARJE+yHX16p5k9D2wYrzrV3f/X4qlEJCstURClZU2dOhUzo6KiQscWi0ijBtqY2YD46wbASsD0+FYerxMRkQw1NTVUVVUxY8aM0FEkC41tKZ4AjAAua+AxB7ZusUQiRWKr7bbm6UefpKysLHQUEWkhjSqK7j4i/vqb3MYRKR7PP/1c6Agi0sKyPU9xHzPrEt8/y8weMLP1cxNNREQkv7Kd+/TP7j7XzDYDfgvcAtzQmCea2Vgz+8rM3sxY19PMJpnZB/HXHvF6M7OrzWyamU3VcUsREcmHbIti3TUUdwbGuPtjQGMPqIwDdqi37jTgWXfvDzzLktM7dgT6x7cRwPVZ5hSRQJ546onQEUSaLNuiOMPMbgT2BR43s3aNfQ13/yfwbb3VQ4Hx8f3xwO4Z62/zyItAdzNbKcusIpJn6XSaKVNfDx1DpMmyLYq/A54Ctnf374GewMnN2P8K7v5FfP9/wArx/Qrg84ztpsfrRKSAuTsvvvRi6BgiTZbtjDYrAY+5+0Iz2woYBNzWEkHc3c2s8XPOAWY2gqh7lW7durVEDBFphmQyybibbg0dQ6TJsm0p3g+kzGx1YAywMnBHM/b/ZV23aPz1q3j9jPi16/SJ1/2Eu49x90p3r+zYsWMzYohIS+nRo0foCCJNlm1RTLt7LbAncI27n0zUemyqicAh8f1DgIcz1h8cj0LdGJid0c0qIvUcMOyg0BFEWoVsu09rzGx/4GBg13hdo6bmN7M7ga2AXmY2HTgbuBC4x8yGA58SHbMEeJzoMlXTgHnAoVnmFGlTDj/sD6EjiLQK2RbFQ4EjgNHu/rGZ9QMmNOaJ7r7/Uh7apoFtHTgqy2wibdb4v41nqy22DB1DpOhl1X3q7m8DpwKvxssfu/tFuQgmIo3XqVOn0BEASKVSHHPisaFjiDRZttO87QpMAZ6Mlweb2cQc5BKRLAwfNjx0BCAafXrowTraIcUr24E2VcAQ4HsAd58CrNqiiUQka+uvNzh0hMUKKYtItrItijXuPrveunRLhREREQkp26L4lpn9HkiaWX8zuwb4Tw5yiUgWFixYEDqCSKuQbVE8GhgILCQ6aX82cFwLZxKRLO26126hI4i0Co0+JcPMkkRTvP0GODN3kUQkW5Meezp0BJFWodEtRXdPAWkz0ySjIiLSKmV78v4PwBtmNgn4sW6lux/ToqlEJCubb7Ml/3r2H6FjSIbS0lKqqqqoqNAFfopJtkXxgfgmIgVEBbHwDBo0iOrq6tAxJEtZFUV3H29mZcAAwIH33H1RTpKJiIjkmUXTjDZyY7OdgBuBDwED+gEj3f2J3MRrvPLych85cmToGCJtWjqdZvz48Rx6aOua1aaqqmqyu1dm85zKykpXSzE8M8vqZ5dtUXwX2MXdp8XLqxGNSB2QddIWVl5e7jM/+ix0DADOuWg0Z59aGAN0laVh51w0mjNOOC10DADOv/xCTj72lGa9xpNPP8kO2+3Q7CyXXHUxpx3fvPfl0ccfZZeddml2lguvuJA/n3xGs1+nJSQ7lakoFqlsi2K25ynOrSuIsY+AuVm+hoi0sNemvBo6AhC1FN94643QMUSaLNuiWG1mj5vZMDM7BHgEeMXM9jSzPXOQT0Qa4fRTCqNF5e68/MrLoWMAMG/ePPb+/e+WvWGOTJ06FTPDzOjTp0+wHJKdbItie+BLYEuiCwbPAjoQXXC4+f0lItIkr1QXRiFKJpOMvXFs6BgAdOzYkZuuuzHY/mtqaqiqqqKqqooZM2YEyyHZyXb06c+OnptZmUagioR1z/33smHlkNAxAOjRo0foCIsVUhYpDtleT/F5M+ubsbwh8EpLhxIpBoeNHM6iRYXxefCSCy4JHUGkVcj25P0LgCfN7GqgAtgJaF1jr0UaadQZf6asrCx0DBFpQdl2nz5lZkcAk4CvgfXd/X85SSZS4Pqu0jd0BBFpYdl2n/4ZuAbYAqgCnjeznXOQS0REJO+yHX26HDDE3f/r7jcC26PrKYpIhllfzwodQaTJsiqK7n4cgJmtGS9/6u7bNieAmR1vZm+Z2ZtmdqeZtTezfmb2kplNM7O74/lWRaTApVIpjjzmyNAxRJos2+7TXYEpwJPx8mAzm9jUnZtZBXAMUOnu6wBJYD/gIuAKd18d+A4Y3tR9iEj+JJNJ7r3j3tAxRJos2+7TKmAI8D2Au08BVm1mhhKgg5mVAB2BL4Ctgfvix8cDuzdzHyKSJwsWLAgdQaTJsi2KNe4+u966dFN37u4zgEuBz4iK4WxgMvC9u9fGm00nOv3jZ8xshJlVm1n1vHnzmhpDRFpIKpVi9312Dx1DpMmyLYpvmdnvgaSZ9Teza4D/NHXnZtYDGEp0CapyoBPQ6Kn+3X2Mu1e6e2XHjh2bGkNEWkgymeTJR54MHUOkybItikcDA4GFwB1ELbvjmrH/3wIfu/ssd68BHgA2BbrH3akAfQBNHCgiIjmX7cn784Az41tL+AzY2Mw6AvOBbYBq4O/A3sBdwCHAwy20PxERkaXKtqXYotz9JaIBNa8Cb8R5xgCnAieY2TSicyNvCRZSZCnuvOcu0ukmH1KXHKutreWue+8OHUOKTLZzn7Y4dz8bOLve6o+IRrmKFKwffpiLu4eOIUvh7sz94YfQMaTIZHueYq9cBREpNocfdjjJZDJ0DAC226lZc2i0SqWlpRx+qE5xluw0qiia2a5mNgt4w8ymm9mvc5xLRLLw2MOPh46w2FOTngodQaTJGttSHA1s7u4rAXsRXUJKRApEaWlp6AgApNNpXp3yaugYIk3W2KJY6+7vwuLBMV1yF0lEilUikeD0k08PHQOAefPmsc8B+4aOIUWmsQNtljezE5a27O6Xt2wsEclG9eRXqPy/DUPHAODl6pcZUhl+nFzHjh0Zc+0NwfZfWlpKVVUVABUVDU7KJQWosS3Fm4hah3W3+ssiElChnHqQTqe59/7CmRC8R48ewfY9aNAg3B13Z/r06cFySHYa1VJ093NyHUREmm7WrMK4hqG789Wsr0LHEGkya8x5Vmb2tLtvF98/3d0LbqBNeXm5jxw5MnQMaUO++eYblltuudAxAPjuu++Ctooyffzxx/Tr1y90jBZVVVU12d0rs3lOZWWlV1dX5yqSNJKZZfWza2xRfM3d14/vv+ruGzQjY06Ul5f7zI8+Cx0DgHMuGs3Zp7bUTHjNU2hZ/nzyGaFjAHDeJedz8rGnNOs1/nj0EVx12dWUlTXvGtiXXHUxxx15UrNeo6Vc+ddLOeFPJzf5+dFFho/gxutuanaWy6+9hDNPLIxBO6Vd2qkoFqlsi2Jjjylq2g6Req6/5oZmF8SWcuIpx4WOsFiXLl1DRxBpssaOPl3VzCYClnF/MXffrcWTiRS4WV/Ponev3qFjAHDQgcNCRwCiS0ftsuMuoWMUhKlTp2JmP1lXUVGhQTcFrrFFcWjG/UtzEUSk2Bx7wjGMu3l8QbQWBw8aHDoCEHWf3nDT9Wy15W9CRwmupqZm8SkZdeovS+Fp7OjTf+Q6iEixueO2O0NHKDjJZJK7/nZP6BgiTdbYuU+HmtlRGcsvmdlH8W3v3MUTkcZYsGBB6AgirUJjB9qcAmQeR2wHbAhsBfyxhTOJSJb23n+P0BFEWoXGHlMsc/fPM5ZfcPdvgG/MrFMOcolIFnbW4BaRFtHYluJPzgp29z9lLBbG8DuRNmz53suHjiDSKjS2KL5kZofXX2lmI4GXWzaSiGTrhjHXh44ARKNPt91pm9AxRJqssUXxeOBQM/u7mV0W354HhgHH5SibSEG7+967SafToWMA8IfDRnDv/fcEz5NMJpn0+LNBM4g0R2NPyfgK+LWZbQ0MjFc/5u7P5SyZSIGbO3cOjZkmMR/mzJlNMpksmDwixaqxA20AiIugCqEI8ODEhzj4wENIJpOho3D4cE2GL9ISGtt9KlIQJj37DLW1taFjAPDgPQ8WxGw2hSKdTvPMc5N4+pmnQkcRabLgRdHMupvZfWb2rpm9Y2abmFlPM5tkZh/EXwvjmjgS3NQ3prJo0aLQMQBUEOupra3l9alTeG3Ka6GjiDRZ8KIIXAU86e4DgPWAd4DTgGfdvT/wbLwswonHnUDHjh1Dx5AGlJWVcdzRJ1A9+ZXQUUSaLGhRNLNuwBbALQDuvsjdvyeagHx8vNl4YPcQ+USKxeRXC+O6fclkkj+OPGrZG4oUqNAtxX7ALOBWM3vNzG6OZ8hZwd2/iLf5H7BCQ082sxFmVm1m1fPmzctTZJHCc/d9d4WOAETHFZ+e9GToGCJNFroolgAbANe7+/rAj9TrKvVojHmD48zdfYy7V7p7pbrUpC27+PzCuKJbIpHgwtEXh44BwLx58zj8SI3KleyELorTgenu/lK8fB9RkfzSzFYCiL9+FSifFJgPP/owdAQpEh07duSMUzQcQbITtCi6+/+Az81szXjVNsDbRFfkOCRedwjwcIB4UoAuv/pK5syZEzqGFIl+ffuFjiBFJnRLEeBo4HYzmwoMBs4HLgS2NbMPgN/GyyJ07dIFMwsdA4BTzjilYE4PKRTpdJqTTjsxdAyRJstqRptccPcpQGUDD2lWYfmZ3+21D507dw4dA4Adt99R5yrWk0gkOHD/A0PHEGmyQmgpijTa+oPXL5iW4k23jFFLsQGD11s/dASRJgveUhQpVuuus27BFOhCsHDhQq685gpqa2tJe5o/nz4qdKSgSktLqaqq+sm6ioqKMGGk0VQURZpooyEbU1KiP6E6paWlbDRkIxYsWKBuZWDQoEFUVxfGpArSeNZaLjVTXl7uI0fqnCSRkNLpNLfddhvDhg0LHaVFVVVVTXb3hsY+LFVlZaWrKIZnZln97FrVx9yzTz0zdAQAzrlotLI04JyLRnPGiaeHjgHA+ZddwJ9GntCs17h57I0MO3h4s1uL1954OcMOOrpZrzH21hs57NDmfygcN+EaDjvkmCY/P51Ok6ADhx3S/Cxjx1/N8Ued1OzXaQn1u0Gl9dJAG5Em6tWrN4lEYfwJ9erVO3QEAMyMXr16hY4h0mSF8RctUoR2323PgimKu+26Z+gILFgwn333341bx41h6B7bho7DvHnz2GWPnYLtf+rUqZjZMm99+vQJllF+rlV1n4q0VQ89fC+7D90naIb27Ttwz12PBM2QqaysjEMPPizY/mtqahrV7aqu2cJSGB9zRaRZ5syZHTpCwXF35szVlICSHRVFkSbaa9+hBXPy/sEH/SF0BBYsmM8BB+3JvvvvFjoKELXUHn2scFquUhzUfSrSRHdOuLdgzsd7/h/PstWWYWdGbN++A7feUhjXdYToKhl3Trg7dAwpMmopijRRoRREgCmvTw4dAYjek0J6XwopixQHFUWRVuD111/l8JEHUlNTEyzDggXzGXHEQRw6fP9gGTLNmzePgw/T5OSSHXWfirQCV1x2PRBNtRZK+/YduPjCq4Ptv76OHTty5WWFk0eKg4qiSCvQvXuP0BGAwslRp2ePnqEjSJFR96kUlaOOPYo5czTMvlClUimOPU5zEEvxUktRispJx59E165dQ8fgo48/BKJpzfr1XTVwGvjkk4/oWwA5kskkJ554RtAM6XSazz7/jHQ6DUTnK6626mpBM0nxUFGUotKvb7/QEQC4+NILKCsro6SkhHPOHk2Xzl2C5rny6ou48vIbg2ao86uVVwm6/6+++orrrr+G2tpaamprqK2p5YbrxgTNJMVDRVGkCW647ubQEX6iS5fwredCseKKK3LJhZcxf/58/nLheYw+5/zQkaSI6JiiSBN88803fPPNN6FjAPDtt9+w/Xa78O23hZOnEHTo0IF99vxd6BhSZNRSFGmCo4//I6Ulpdx0w61BTxBPp9OcctoxuDslJSVce/UtQU/LSKVSnHzqMdxy0+3BMmQavN7g0BGkyKilKEXlkssv5ccffwwdgztuu4d11xmEmQXNkUgkuHnM7Qxce11uvP62oAVx4cKFXPfXKxg4cF0uu1xdllKcCqKlaGZJoBqY4e67mFk/4C5gOWAycJC7F8bMyxLUeusOon379qFjALBh5RBKSgriT4jBg/8vdARKS0vZYIMNWbhwYdDiLNIchdJSPBZ4J2P5IuAKd18d+A4YHiSVFJzrb7qhIFqKAFtu8ZvgLcU6N950begIJBIJNtl4M266+Tq22Pw3oeOINEnwomhmfYCdgZvjZQO2Bu6LNxkP7B4knBScB+95oCDOUwS45dYx1NbWho4BwN13TAwdAYgK4/bb7Rw6BgCLFi3i5ltvCh1DikzwoghcCZwCpOPl5YDv3b3uv810oKKhJ5rZCDOrNrPqefPm5TyoSKbllutFIlEIf0KFZbnllgsdAYgK9HI9CyOLFI+gf9Fmtgvwlbs36bo37j7G3SvdvbJjx44tnE7kl908tnBainvsvX3oCEA0GvbW8YXROlu0aBG3jCus80ml8IUeJbApsJuZ7QS0B7oCVwHdzawkbi32AWYEzCjSoEcfejJ0hMUevO+p0BGAaJq3B+8rjPelrKyMQw8+LHQMKTJBi6K7nw6cDmBmWwEnufsBZnYvsDfRCNRDgIdDZRSR4uTuzJ4zO9j+S0tLqaqqWuZ2FRUNHh2SQAr1gMipwAlmNo3oGOMtgfO0afsddABffPFF6BgFZdc9dmSvfYeyaJHOFKqzYMF8DjhoT/bdf7fQUfj2u2858NDf89jjjwbLMGjQINx9mbfp06cHyyg/Z+4eOkOLKC8v95EjdcmaXKitrS2Y8/EKxbvvvouZscYaaxTEaRkffPABZsZqq60WNE9tbS3Tpk1jwIABwTJkZgFa5He3qqpqsrtXZvOcyspKr66ubva+pXnMLKufXav6T3fmiaeHjgDA6Msu4MyTCiTLpRdwRoG8L+dfdgF/GnlC6BgAXHvj5eyw6x+a/PzPZ15GSUkp2+44jJJmnqj+5CM3s90uzTsV97MbLyeZTLLNjsOaVQSefvQWNtv+0CY/P51O8/mt1zTrNeq88NSt7LTb4c1+nZbQmG5QaR0KtftUpKAdPvJE3n37dZzC6Gl55+3XOfQPxwZv0ScSCQ4afmzQDCLN0apaiiL5dMaoSyktDTcZeJ2335rCvr//A2+/NYVEIsGaA9bNexdqOp3mww/eIp2OjpMBDFh7vbxmaMjbb7/O2oFyTJ06NSc/h4qKCh2HzCEVRZEm6tatR+gIpNNpnnriQTydJllSQjKZZJW+q9OhQ37P250//0eef+ZRUqkU6VSKVDoVvCguWDCfZ55+JFhRrKmpyUm3q7pyc0vdp1JUTjzlOObOnRM6RsFIJBIcf9I5zJnzPceecDZ/OvasvBdEgE6dunDw8OOZO+d75sz5jqNPPDfvGepr374Dxxx3VugYUmTUUpSicsyfjiuYq8xP//wT+qzcN3QMAHbbY//QEWjXvj0HHnYMX36huTakeKmlKEVllV/1DR1hsXFjr6GmpjDOU5z40F2hIwCw/AoVPPnoPaFjiDSZiqJIE1x60Zl07twldAzS6TSXXzKKLp27ctVlVVxz5XnMnx9mcvxFCxdw03Xn07lAWvILFszn6ivOCx1Dioy6T0WaYMutdsSM4KNPE4kEO+2yNz/MnUvnLl0wM9q37xAkS1m79myz/e7MDTi1Wqb27Tuwyaa6rqNkRy1FkSa4646buPeecQXRfTpgrUHcdcdNDFhrUJDTMTKtuvraPHjvuGD7z7Rg/nweuH9C6BhSZFQURZrgsqvGs+aa62CEn+INYMBa63LrzVcGv5RVMpnkvIsL43JN7Tt04IKLbgwdQ4qMuk9FmmjNtdYlWQBzwi5atJC11l4PMyOV0jy1mRYtWkhZWbvQMaSIqKUo0kR333ELtbU1QTOk02lOOPog7phwI3ffcQtnnDKCefN+DJqpUCyYP5+TT2jenLLS9ugjpRSV8RNuZd999qd9+/aho3D1X+8IHYFEIsG1N4Y/BWLBgvlMeuJ+PJ0mlU6BR9l22eNAkslkkEzJkiRbbLVdkH1L8VJRlKLSs0dPEgl1cBSaRCJB12498HR68XUCzSzoz8osURBT8UlxUVGUorLrLkNDRwBg+CG7UlpaxjXX3xn8tIxCUFbWji233plUKsXxR+zN1Tc9GDTP7NnfcdbpR1FbW8Nvt901aBYpLiqKGZ557hkGrTuI5XsvHzqKFLhbxj8SOsJiTz/5IO7ELbMkW229AyUlzbvGY1Mlk8ngBRGiydqvKYDubSk+KooZfvjhR2pqwg6cEMmGuzN3zuyoKCaMhCVIp9OhY4kULRXFDLvvVhhdc9KwvfcdSklpKTddPzbYpODuzpEj9wF3SkrLKC0p5bwL/kqnTp2D5DEz9vrdMI4csTd/HXNfkAyZUqkUpxzzey677u5gGb78ciZXXFpFTc0iampqqKlZxPVj7g2WR4qLimKBeuE//2b99QbTqVOn0FEKxh0Ton9sZWXhjuGZGVdfdwcv/fefmBmWSASbVi3TVdcVRldhMplk7/1HBM2wwgrlnH3OFbz22su4O0M22ixoHikuGsZXoN5+523m6LqBP1FWVha0INYpLS3jnXde5/333+KD996kZlH4qd4KZbBPOp1m2ntvho5B2p0P3n+bD95/q2DeGykOQVuKZrYycBuwAuDAGHe/ysx6AncDfYFPgN+5+3ehcoYwYvjhoSPIL3j37amUlJTwlwuvD/ZPd+HCBVSddTS1tbW4OyUlJZSUlHLGqEvp2DFMD4O78/67U4PsO1OnTp3Zd7/DqDr72NBRpMiE7j6tBU5091fNrAsw2cwmAcOAZ939QjM7DTgNODVgTikAr015lUQiwaB11ws66TXAqHOu5P333gzaCmnXrj2nnXkx72e0zNZYc51gBRGi7tOTzro02P7rpNNpPvlkGnvtdXDoKFJkghZFd/8C+CK+P9fM3gEqgKHAVvFm44HnUVFs8x58+H4SlmDVVVejS+BrGXbr3oP//PvvDN5go6CFsUvXbrzwr2fwdJpkSQnr/9/GwbLUKYQT5uf9+APPTHqEdDrNhjqmKFkI3VJczMz6AusDLwErxAUT4H9E3attxsmnn0JtbS3HHHU0/fr2Cx2nYJx79mhOPOU4cA8dBYBjTxgVOgIQnahO3H0a0sIFC7j60rNIpWpJpVKcee41wbKkPc3cObOjKecCKS0tpaqqqsVft6KiosVfU5YoiKJoZp2B+4Hj3H1OZteYu7uZNfhf0MxGACMAunXrlo+oeXHUEUcyY8YMKsr1y5/p408+Yo/d9wp2OkZ90z//hD4r9w2aIZFIsOvQ/RYvhzppH6Bd+/YccOgxfPW/GSy/Ytjf3a5du7P9DruT9nDnbA4aNIjq6upg+5emCT761MxKiQri7e7+QLz6SzNbKX58JeCrhp7r7mPcvdLdKzt27JifwHnQd5W+PDXpaWbMnBE6SkG54qpLuf+Be5lbIKNyb73l6uAXGU6n0zz68N088tCdPPHovcGv2rHCihU8+ejdlFf8KmiOOXO+5+mnHuapJx8KmkOKj3nAriiLmoTjgW/d/biM9ZcA32QMtOnp7qf80muVl5f7yJEjc5o3n1544QXWXXfdgmgBP//882y44YbBz5l86qmnFk8yvckmmwTJ4+488cQTuDvJZJJEIsGWW25Ju3Zhr9n32GOPkUwm2XbbbYNdlQKiIv3kk0+y0047Bcvwww8/8NJLL5FOp39y23nnnZv8mlVVVZPdvTKb51RWVrpaiuGZWVY/u9Ddp5sCBwFvmNmUeN0ZwIXAPWY2HPgU+F1jXuzcJ9/NRcasjdphQLOz+LwSnv33dCzxxbI3XlaWp95vXpYfO/CPf32OWfM6FkZtvwZ/+bh3k5+f7lwJGLTryH//V4E14woMZ/Wbxdy1mjaDUf92/Vk47wcA2nfqysJ+A1jUjNGwXd55mE23H9bk5wN06Lk6ABts2LxBJf9+ahyDt85+xOac2UvOmNqpe39W7b82AF2bMehmynO38exyW2f9PO9eyw+JAaTm/0D0od8p6dCZZ5cb0OQsUNWM50oxCT369AVgaf9NtslnlkJjHbuHjrCYdQo/mhDA3/sXnkiQ3OT3zSqIzbXiagN58KJjSCRL2eW4C4OfHgLwwF1jKSkpYdD6G+X9uOK3X3/JhBsvobY2mlItVVuLe5pEsoRTz72OZJ4HAFmyhPY9V+TDe68knarFUzUMHHFBXjMATJ06tSB+N9qSiooKpk+f3qzXCN1SFGmU1JTHoGcfLJEg/d6/SKyxGdYu3HHk5futRTJZwiuPjOf/djqA0nZhp3r7y2Vjg+27Z68VGHH8OTx2/3jSniadSuGeZv/Djs97loWzv2bW5OgUlU7lq0bXdkyn+PSxsZBI8KvtD8IS+elerqmpycnoU1m6lni/VRSlKFjvVcGiaY+MBJSEnbprxdUGghlmRjLwNGKLFi1cfL+sLMyxzbKydvRfaz3AqVm0iGSgUbAl7TvSuc8aOI7X1ADRBY8TyehfXb4KohQvFUUpComKtUJH+InqR28jWVJKIllCnwHrU9YhzCCkdDrNmScMI51OU1JSyujLbw1yWkaypIT1KjcllUpRdcLBnHfV7XnPAJBs15Hua2xAzY9zeP/2C0jX1uCpWgYdfWWQPFJ8VBQLlH/9KXRbESsNO6qxkKTf/SduRmL1jbFkuPPxAPoP2YZEIsmgbfciEbD1kUgkuOTaO3no3nGYJZo9GKo5amtrmPTI3Wy8xfbBMgDUzpvL16//g55rbxwNtAl4Ar8UHxXFQpUsBR2k/6l2da2x8O9Lhy7dsUQCK4AsAF279sASFnRgh1mCLl27Bz1hHsCSSUo6dAFP41AwMyBJcQh+8n4hSX9cjf9YGBfj8G8+hZqFy94wD9Lv/hNfND90DBL9/g+m/RdSYU9QB3jtiTuZ8tTd3Df6SBbN/zF0HB66bxwT759AOmCrKJlMstk2u7DFb3cLlgGiLtTe629Fz3V+zf9eeJjl/69ND2SXLKmlmMF69oGy8BeMBbDlVoHAXaf+9afRJ+3effG5s6Lu3IADXNIfT4bVN4la0QG988ITrL/j/tFFhgtgoA3A7nsPY+vtm3beZUtJpVL85/nHcffghTFds4hv3/wPK266G19VP8Pylb8NmkeKh4piplQNBO76WSxVE77bp7ZuCjOLunJD51n4Y9ylHC6HuzN/7vcAJBJJLJHA0+F/Z+bMCd/D4Z7mh7lzSKfCH8NzT1M7b250TLFQ/qalKKgoAukP/hPdsQR89wX07ot1XT5oJp/9JdapZ9DWoq3YH3/3H4BBIgHdlg+Wp/bJK6OfT7KE9Mx3SGy8L9YhzBR4H7z8LLiTKCklWVLKWpvtCIRt1e++z7Cg+4doMvIddz8gdAwWfvcVHz5wDV5bQ7o2Onm/fIs9Q8eSIqGiCNhqG/10uRDOZeq1SmF05a7QP2ohmkFpuDw2cBswi9uIBmX5P3Hf3fn0jZeo3PmgOFTUffrFh29iGCv1Xzf4SfwC7XoszxoHns7saVOjFWopShZUFCmQIliPLZgLHbpCMuyPyOZ9j+NRK61Tj3B5vvkMLBFP72bQu1/ejy2mU7X8b9ob8YwtTjJZAgmLulEtQa+VV1dRLBQO82ZMWzyjzXLr/Dp0IikSKoqx9CevRuczmWG9+mJdegXN4/PnYN3CX1vZKtbGP3wJW3kQlueWq6dqSf3jFiz+pO+WhGQSS5SQrr4fEkmwBInVN8F6rZL7PO58+fG7pOPLMyWSJXEXaglbHXQiHbv1zHmG+tLpNGefcjgeX2T4rNHXBbum4sKFC7jivONJ1daSTqf488X5n3pu/lef8/kzd+CpFJ6qjU7eT9fy5vWnYMkkaw8/Dwv8QVMKm347IDoNo0cFdZOIuRnULsRKwhwnqsvjpR0K4yy4Xn2DHEu0ZAnJX0fHqPzb6YvPwXMM674SlkhGrdg8deuWlJax/na/i/ZZl5EEmNG+S5jjm4lEgt33PXTxcqiCmE6nmf7JNHbc48CgA7La9+5Dv6F/ZGmDsVQQZVn0GwLw3YzouIPFXXMA3VaALoEGT3w3A9JprH3n4N2nADbnK+jYDUrb53/f7TsDkP58Km4JsGjEp/Uox9p1zOuHhlRtDR+99gKpuKWYLCnBEkkskaT3Kv3p0KV7HtNE0uk0L/7rmXiatxLWC3CVDIAF83/kheceI5WqJZ1KkU6nWK+yeZexagozo7RT17zvV1qP8P9xC4DXnYoRF0VbbmWs83J523/6g//GtTgBCYtyWAL/4v2o1Rq3kKzr8lj3lfKXa9qL0b7NYPqbsNKAxUUqHzxVS/q/d4LHQ/zjokgySfr1JyCRiK6Wkadu5mRJKVsedAJPXHtWdEyxpJREMkkiWYIHHMwxd+5sPJ0mWVLC1ZecxR+PHUWHjvmdizWddn78YQ6pVC21tTVBT1NZOPtrPn9qAp5O46na6OYpBhwyKlgmKR5ttij6gh+WLHSviO9E58B5shRfEM1SYiWlOT9h3VZed0muhT8uzrH4/MDSdlExyNFJ4u4enQNYX+9+9bZLR++LsXhEqiVz9/5YsoTEejtG+/7x24xWoeEde0TdqZ2652TfS81kxq9/90fmfP0/ALr2WhEzo0OXMNecTCQSHPKH4/nqy5kALL9Ced4LIkDnLl3Z84CRAMz6cia9VyjPe4Y6pZ2603uDrYN+UJHi1XaL4vQ341ZQYklrqK77dPH8kQade0DPlXOaJbP15V+8G69cksWW+xXWOYddQrUL8ZlvR4XO0/EQdmPx+Ym2pPW6+P0xovudl4OefXIWrW7AU3rKo0u6T5NJEuvtjHXokrP9/pLuK/Th33f/lWRJKduOOCvYZZLqrFTxK8aNuYxkSQknnH5hsBwrlv+KVCrFnWOv5PizLg+Wo+bH75n16nN4OrV4wI2nU3Tvv36wTFI82mxRTKy+MelPpyz5555ZFCFjfe5O13BP459OWXIeIAaJkiVFmiiTf/8FPjtqmVinHi3ehWql7bFVhyxeTn9cXfdI1J1LRmEErOfKWB5aaJ6qJT35oahIl3Vc0n2aSJB+9/lo5OmqQ7CuvXOepb6djzk/7/v8JV26dudPJ54TNMOihQv5282X0iXAsdVM7br1YvXfHU9q0UI+e3ws/Xb/Y9A8UlzaTFH0hibX7rbiklahARkNJDwB5nhJO7xmEZiDR91nTe0udPd4MmtfMjiue71uJqu7km7cfVqXzROA46WZeaJu1uZkWpwt8/3p0UDLr+6YJ+AlJYu3t0QyZyP6LFmC9VnnpxGIT+C3BBikk2Uk0un4/MX8qZvqDQgywKa+bXfaK3QEytq1Y4vf7sYPc+eEjgJAumYB3dcaQs2PsyntFGZ0sBSfNlEUvWYh/smr0YCNui7CulZhxsCWJV2pmf9g4+N77oDjJWXYyoOaFmTBXHzmuz/tpqy7LW4pJuMMyZ9mscXVOuMbi17HkyXYKk3vGvKaBdH7k87MA9HAn0QD709dKxbo3CuafScHPFWLv/fPxdfDyzxPkURi8XmKrL5JzjIszZN/HUUiWUqipITtR44KdpHhOmuu1cTfyRbWb/W1ufDMIxhcuWnQHDU/zuHDe69afJ7iwBEXBM0jxaNVF0VfNA//+jPAo1MK3DEc96h15WT+g69rlWUeU2TxeJcltchJz3wnaiGtuEbjsyyYi3/3RTRLzeKWYhrzzDw0fNyuLo/XrctoaZLGEk0/puU1C/CvPoJ4HtHo/Ulji1tkSzveGj+/diH+xfvR6REtdGzRF80j/d4L4Gms58qLB0yYxUUa+0mR9q8+xmd9jPXul5eT+AH2OPVqXrz/JjYcOiz4McVCYmastuY6y94wh2p+nM3//vs4nSpWi0agplN8+thYVtn5sKC5pDi06qJIsjSamsxZPKzf44LidQUnU906t592X9Y9weJKVNdiykZJu3ikZPxa9fNAvdfMKIp1GTLz1h2HtLh12UjRtFfp+ER4j0756Nhtybe5JFVGV25mDoi6Uf2ncZLtopZm5ontTZ0+L1kWXTorbqFbuqbuBaMGe0ZrdfE+zOIPHPmz4urrkCiA80gLSSKRYM2BYQe0JMva03nlNaLDDES/zSF+TqWlpVRVVeV9v21ZRUXFsjdahlb9F23JUqzbClH36YcvZXRX1g1syTgvsK4rLj5W9bPrL7sDaUiWkui7QfZZSsoWn0/nNQvxj16OsqQb6j5NLCl2P+tCJa47DqXtSWScztEo382A72bGxTijGzedWvLepNNLCmFdnsRSulCXvEH4V+nouCwpwLD+TZtv0pIlWPma0au6k3rqSha3DpNJLFka5UqUYOv8FuvR/D+Epui73iZB9lvo1h+yRdD9J0rb0blPf97/2/mLr5Ix6Jir8p5j0KBBVFdXL3tDKSituijWsdJ22ICl/6H6rI+jk43jf/bWpXdOT1K30nbYmpvn7PV/cd89+2TVzelfTlvcOrOefbA8T/dmZpTscPzi5fT7/8ZWGxIVRpGlKO3UlYEjw52eIsXLPPSFY1uImc0CPm2Bl+oFfN0Cr9MSlKVhytKwQslSKDmg5bKs4u5ZnftjZnOB91pg39I0dT/7rH52raYothQzq3b3ytA5QFmWRlkaVihZCiUHhM1SSO9DW9TU9z+/J3eJiIgUMBVFERGRmIriz40JHSCDsjRMWRpWKFkKJQeEzVJI70Nb1KT3X8cURUREYmopioiIxFQURUREYiqK0ihm2c5r17qZWVn8Ve9LG2dmO5jZe2Y2zcxOa+DxdmZ2d/z4S2bWN0DMVqkR7/0VZjYlvr1vZt8v6zVVFGNmtqaZbWJmpWY5vIhi4/MUQobNzOwgAHf3QikAoXOY2TbADWbWzXVQfqlC/5ziDDn9Hxf/nV4H7AisDexvZmvX22w48J27rw5cAVyUy0xtRWPee3c/3t0Hu/tg4BrggWW9bpuY5m1ZzGxP4HxgRnyrNrNx7p73C8OZ2Rru/r67p8ws6R7PHJ7fDAmgI3BjtGid3P2GuDAmvO6yFfnLs0GcZ5G7vxyyEJnZ9kSj2j4GegKzQ7wnDeQqhAybEr0nSXd/KNTPycxWATq7+1t174mZWY7yDAGmuftH8X7uAoYCb2dsMxSoiu/fB1ybwzxtSWPe+0z7A2cv60XbfEvRzEqBfYHh7r4N8DCwMnCqmeX1sgtmtgswxczuAKgrjPnMEO837e4/AOOBW4Bfm9nxdY/lM0v8ntwCjABOMrOR+dx/A1nOAw4ApgIXQ/7fk4w8m5jZfnUZct0qWkaWnYAbgHWBs8xs/4zH8tZiNLO9gKeB68zsTjPb3cw657CnowL4PGN5eryuwW3cvRaYDSyXgyxtTWPee2DxB6V+wHPLetE2XxRjXYH+8f0HgUeBUuD3+fqDNrNOwJ+A44BFZvY3CFcYY7VEHxDGA0PM7HIzu8AiOf/dMbP1iVrww9z9YOBeYECu97uULAngYOA0d3+BqDh2NLPfBsqzHfAMsLeZHQ6LC2PeuyzNbB3gL8CR7n4+cCfR73C/OFdeut7jv6GDgQPcfSvgRWBL4IC6wpjrDFKw9gPua0zPW5sviu5eA1wO7Glmm8ef+l8ApgCb5THHj8BhwB3ASUD7zMKYrxz1PAz8z92fBaqBI4CuHslH66gD8Fd3fz1efg3Y1MxWznPrYx1gHeBEd38uLpDzgHeB9eNt8l2M+gGjgbFAZUZh9Hy1GOt9zyPc/V9m1hs4laj35UYzu74uVx4ipYHuwGrxPq8CXib6ILVlnLml35sZRB8c6/SJ1zW4jZmVAN2Ab1o4R1vUmPe+zn5EH9aWzeMrv7flG9CeqJU2BtgiY/1zwOBAmZYD7gf+Fi9vAAzIc4Zy4FbgcOADYBTwCDAyx/tdI+N+7/hrkui44iNEhRmgfx7egx2BN4g+IDwBVGQ8thnRLPy/DvD7URrfOsd/8Ddk/lyARJ5+Pwwoi5dLgIOIWmoAXYD3gZ1ynMOIjmMC7E00+GKDjMdOBu7O0b5LgI+IPqSUAa8DA+ttcxRwQ3x/P+CefP++tMZbY977eLsBwCfEk9Us66aBNoC7LzCz24ku33u6mQ0AFgIrAF8EyvRNfPzsEjN7l6go/CbPGWaa2efAn4Gj3P0RM/sNMC1X+4yP291jZhPdfT93nxUPIkmZ2QLiwWHxqNh9zOwQd/8uR1m2Aq4CDnT3l83sQWAtYIaZlbj7C2Z2MbCNRTPyL8pFjow8XT0e/OVRDwdAjZk9TvTPf0sz24OoUM4FHsphlh2IPiRNAb40s0vd/Uczu8/d58eDxOaa2RPAjznMMRTYPb5/M/ASMBDYLR7MMpnob+hZM1vV40EZLcXda83sT8BTRH+jY939LTM7F6h294lEx8QnmNk04FuiwijN1Mj3HqL3+y6PK+SyaJq3DBade7YpMBJYAFzl7q8FznQ8UXfUtu7+RoD9rwwsH/9zyekox/iY0P1Ew6Z/DZS4+4HxY0mi7v47iAYqDAYOdveljTRriTxrASu6+9/NbEXgVaLuuC+JjlfdRjSibU/g9+6+IIdZ9iQaOTcCeMXrjaqMB4VtCFwArAlslqvfFzPbFrgMOJboePyuRMcTF8WPJzw6vrkPcBqwt7t/nIMc6wF3AycAqwBHEn2Amw1sTtS19k+iY+NnAZu7+7ctnUNaFxXFBsT/gD1X//yzyNEDuIfoWNbUwFnyMoTczMqBOURd2jcAC+oKY/z4Q8AawB7unrcLuJrZmUR/L38xs2HATsAR7v6tmS3v7l/lcN99gb8B84mORV0KTK7/8zCzk4Hjgd/m+MPCaOA5d382HkzzBNEHma+Be9398/gY5x+Bg9z9rRzl2B74o7vvHi/vRHQY5EaiDy0bEX2ImAtcHPoDrhQHFcUCZ2btc9kCKWRmthzRcd757n6gmfUHDiU6zpqzf/qNzPY4UOXuL+dhX78C+rn7P8xsFFAJnAtMibuQ6lqLJwPPuvurOcqxDtGxm1lx4etCNCL4n0THXTcgap0dDmwDfJiLFmJGnhWAS4iOIb4St053itf9wd3/a2btiD7g5rRrW1oPFUUpaGbWi+if3Kbxqs3d/cs8Z/hJK9mic+HOAHZ29//lcL9ruPv78f1u7j47vv9nohOXz3X3V8xsPV8yQjdXWXYkmonlY6LjuofHx5wX79uik/cPdfc/5DDHRkS9CD+6e3V8/KgzcDXweXzs+RigL1EPi/7BSVba/CkZUtjc/WuiE+W7AnvluyDGGaJhjNEclsOJWmmH5Lgg1k3kcGecYXZ8zBt3P4/o2ObxZnYh0SCO5XOYZSuiAUd/cPehwCKik/QhaiHWWRWoMLMuuThFJS7MfyOaPGGUmV3u7qOIRiUfQ3QcGqJTM9qpIEpTqKUoBa3AjquWAtsSdQvm7HjmMgYctXP3hfH954mOr26fy0FYyxhw9F+iqcuGE504f3AujiHGx/lvBx5z9wnxwKKngTfc/fC49bwm0IuoC/f3uW49S+ukoigFry0eV23EgKM1iEZeDsvnP/8GBhztQDSv59nA+TkuzqcCM919Qsa6/wD/dveT4w9Q6wAfu/v0XOVorcwsRdTyLwHeIeoNmRc2Vf6pKIoUuAYGHA0m6k5+O+5eDpntSaLTMVr0/L+M1888rnog0SkeO7n7Z/G6XkTvzVmhB18VOzP7wd07x/dvJxrhfHngWHmnY4oiBc7dvyE6d7bGookc7iW6OkBeC2L944TxgKPeRKeK5GJ/dcdV7wJw978RzU3873hEbt0x50VEg22k5fwLWN3MzjWz4+pWmtloMzs2XKzcU1EUKQIZA466A3u6+8wAGZY24KjFZ32yn06QvyBjwNGfgXHAI2Y2Mu7OHQTMaukMbZVF87PWTW84luhYcd28sfsRDXZqtTTNm0gRiI+X7QRsF2Jmo3rSRNMf7pmrAUceTRl3GNFx1QeILup8p7vv7+5/NrNXgRWJ5r38XS7Ph2xDOpjZlPj+v4Bb3H2RmX1j0RVrVgBei3suWi0dUxQpEm1xwFGdjOOqi9x9fzMbCPzg7p8GjtZqZB5TrLd+X6JR0CsC49398byHyyMVRREpChkTOfyaaALorTTKtOX8QlEsI+pKLSW6Mk2oS9nlhbpPRaQouPvXZjaV6HjXtiqI+RF3of4d+L61F0RQURSRIlFgx1VbnYZaibB4gM3GwD75TRSGRp+KSFHw6LqZu4ae2agtMbO1ia6f+qy7fxA6Tz7omKKIiEhMLUUREZGYiqKIiEhMRVFERCSmoigiIhJTUZRmMbOUmU3JuPXN4rlbmdnsjOc+k+W+q8zMzWz1jHXHxesqs3mtX9jH+mZ2SzNfY2MzuymL7ceZ2d711vWOr0ghIjmk8xSluea7++BsnxRPOgzwL3ffpRnPf4NokuK/xMv7AC15kdszMl67qXYEmlXQ3H2WmX1hZpu6+7+bmUdElkItRWlxZvZJPCUXZlYZXyG+rmU3wcz+DUxYynP7mtmbGcsnmVlVfP95M7vSzKqBusvXPAQMjR9fDZgNfJ3x/OvNrNrM3jKzc+plPMfMXjWzN8xsQANZugCD6i7iG+c/KePxN+taxmb2ZzN7z8xeMLM7M7cDtgGeMbNhZvaQmU2K9/8nMzvBzF4zsxfNrOcy3tqHgAOWsY2INIOKojRXh4zuzwcbsf3awG/dff94efOM55/ZiOeXuXulu18WL88BPjezdYhajHfX2/5Md68kurzQlmY2KOOxr919A+B64CR+rhJ4s4H1P2FmGwJ7AesRtQorMx7rBdS4++x41TrAnsCGwGhgnruvD/yX+BI9v6Aa2HxZeUSk6dR9Ks2VbffpRHfPvCjtT7pPG3FMsn7RA7iLqCBuT9QqOzTjsd+Z2Qii3/WViIpy3YwoD8RfJxMVqvpWonHX6dsUeDi+gsUCM3sk47HtgKczlv/u7nOBuWY2G6jb9g2iwv1LvgLKG5FHRJpILUXJhVqW/G61r/fYj1k8t7HPfxQ4CPjM3efUrTSzfkQtwG3cfRDwWL3XWxh/TdHwB8T5Dew/8+rzpUv5HjLVP564MON+OmM5vZQMmdqTo6vci0hERVFy4RPg/+L7e2X53C+B5c1sOTNrByxzEI67zwNOJeqOzNSVqIjONrMViApUNt4BVq+3bkMAM1sF6Et0CaN/A7uaWXsz61yX2cyMqPU3Jcv9Ls0aNKI7V0SaTt2nkgvnALeY2XnA89k80d1rzOxc4GVgBvBuI593VwPrXjez1+LX+JyoeGWT5V0z62ZmXeIuT4AeZvYK0YCe+4GL3H1vM5tI1C37JVFX6GyiDwavedMmGL7RzK6M73/u7psAvyFq7YpIjmhCcJFfYGbHA3Pd/eZ4FOwP7n5pA9t1dvcfzKwj8E9gBNFljqY1VLCbmOWfwND4ahEikgNqKYr8sutp3HXkxsSX2WkPjHf3V4FXWyqEmfUGLldBFMkttRRFRERiGmgjIiISU1EUERGJqSiKiIjEVBRFRERiKooiIiKx/wdVgw+eXeM1UwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 504x504 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "sub_data = np.zeros((0,2))   #due to the large dataset, we subsample the data for the plotting\n",
    "doses = [0, 1, 4, 12, 37, 111, 333, 1000]\n",
    "for j, i in enumerate(np.unique(input_output[:,0])):\n",
    "    temp = input_output[input_output[:,0] == i]\n",
    "    temp = temp[np.random.choice(np.arange(0, len(temp)), 100) ,:]  #subsampling\n",
    "    temp[:,0] = doses[j]\n",
    "    sub_data = np.vstack((sub_data, temp))\n",
    "sub_data = sub_data.astype(int)\n",
    "\n",
    "\n",
    "input_binning = len(pd.unique(input_output[:,0]))  #input binning number is the same as the number of input doses\n",
    "output_binning = 10  #output binning number\n",
    "\n",
    "Pxy, x_margin, y_margin = np.histogram2d(input_output[:,0], input_output[:,1], bins = [input_binning, output_binning])\n",
    "\n",
    "x_margin = np.arange(-0.5, len(x_margin)-1)   # for the alignment of the plot\n",
    "\n",
    "\n",
    "###entropy calculation\n",
    "Pxy = Pxy/Pxy.sum(axis = None)\n",
    "\n",
    "Z = Pxy/Pxy.sum(axis = None) + 1E-10\n",
    "x_marginal = Z.sum(axis = 1) \n",
    "y_marginal = Z.sum(axis = 0) \n",
    "x_entropy = round((-x_marginal*np.log2(x_marginal)).sum(), 2)\n",
    "y_entropy = round((-y_marginal*np.log2(y_marginal)).sum(), 2)\n",
    "joint_entropy = round((-(Z) *np.log2(Z)).sum(axis = None), 2)\n",
    "mutual_information =  round((x_entropy + y_entropy) - joint_entropy, 2)                                                              \n",
    "                                                                       \n",
    "\n",
    "###plotting\n",
    "fig, axs = plt.subplots(2, 2, gridspec_kw={'width_ratios': [2, 0.7], 'height_ratios': [0.7, 2]}, figsize = (7,7))\n",
    "#plt.rcParams['figure.figsize'] = [10, 10]\n",
    "\n",
    "fig.subplots_adjust(hspace=0.07)\n",
    "fig.subplots_adjust(wspace=0.07)\n",
    "\n",
    "axs[0,0].get_xaxis().set_visible(False)\n",
    "axs[0,0].bar(np.arange(Z.shape[0]), x_marginal, color = \"gray\", edgecolor  = \"k\")\n",
    "axs[0,0].set_xlim((-0.5, Z.shape[0] - 0.5))\n",
    "axs[0,0].set_ylim([0,max(x_marginal) + 0.01])\n",
    "axs[0,0].set_xticks([round(0,1), round(x_marginal.max(),1)])\n",
    "axs[0,0].set_ylabel(\"Px\")\n",
    "\n",
    "axs[0,1].get_xaxis().set_visible(False)\n",
    "axs[0,1].get_yaxis().set_visible(False)\n",
    "axs[0,1].axis('off')\n",
    "axs[0,1].annotate(\"Hx: \" + str(x_entropy)  , (0,0.8))\n",
    "axs[0,1].annotate(\"Hy: \" + str(y_entropy) , (0,0.55))\n",
    "axs[0,1].annotate(\"Hxy: \" + str(joint_entropy)  , (0,0.3))\n",
    "axs[0,1].annotate(\"MI: \" + str(mutual_information)  , (0,0.05))\n",
    "\n",
    "c = axs[1,0].pcolormesh(x_margin, y_margin, Pxy.T,cmap ='PuBu', edgecolors='gray', linewidths=0.5)\n",
    "#cbaxes = inset_axes(axs[1,0], width=\"30%\", height=\"3%\", loc= 2 ) \n",
    "cax = axs[1,0].inset_axes([0.1, 0.9, 0.3, 0.03], transform=axs[1,0].transAxes)\n",
    "\n",
    "cbar = fig.colorbar(c, cax = cax, orientation = \"horizontal\", ticks = [0.0, round(Pxy.max()-Pxy.max()/7, 2)])\n",
    "cbar.set_label('Pxy')\n",
    "\n",
    "#cbar.ax.set_xticklabels([0, round(Pxy.max(), 2)])\n",
    "axs[1,0].set_xlabel(\"FurFurMan (\\u03BCg/mL)\")\n",
    "sns.swarmplot(x=sub_data[:,0], y=sub_data[:,1], size = 1, color = \"k\", ax = axs[1,0] )\n",
    "axs[1,0].set_ylabel(\"GFP expression\")\n",
    "axs[1,0].set_xlim(-0.5,Z.shape[0]-0.5)\n",
    "\n",
    "axs[1,0].set_xticklabels(doses, rotation=45, ha=\"right\",  rotation_mode=\"anchor\") \n",
    "#axs[1,0].set_xticks([0, 1, 4, 12, 37, 111, 333, 1000])\n",
    "\n",
    "axs[1,1].get_yaxis().set_visible(False)\n",
    "axs[1,1].barh(np.arange(Z.shape[1]), y_marginal, color = \"gray\", edgecolor  = \"k\")\n",
    "axs[1,1].set_xlabel(\"Py\")\n",
    "axs[1,1].set_xticks([round(0,1), round(y_marginal.max(),1)])\n",
    "axs[1,1].set_ylim((-0.5, Z.shape[1] - 0.5))\n",
    "print(\"Entropy of input = \" + str(x_entropy))\n",
    "print(\"Entropy of output = \" + str(y_entropy))\n",
    "print(\"Joint Entropy = \" + str(joint_entropy))\n",
    "print(\"Mutual Information = \" + str(mutual_information))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3a83567",
   "metadata": {},
   "source": [
    "The plot shown above describes the joint probability distribution of the input and output dataset."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "cfdfd724",
   "metadata": {},
   "source": [
    "# Mutual information optimization function\n",
    "To optimize (maximize) mutual information using weighting values (see SI section 1 of the Supplemental Information of the paper), we use Sequential Least Squares Programming (SLSQP) (https://docs.scipy.org/doc/scipy/reference/optimize.minimize-slsqp.html#optimize-minimize-slsqp)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "765cc68e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def channel_capacity(data, binnum_mem, binnum_GFP):\n",
    "    global history\n",
    "    global entropy_list\n",
    "    global Pxy_origin\n",
    "    np.random.seed()\n",
    "    entropy_list = np.zeros((0,4))    \n",
    "    p = np.repeat(1/binnum_mem, binnum_mem)\n",
    "    history = []\n",
    "\n",
    "    #binning X\n",
    "    starting_data_input = np.min(data[:,0])    #Low cut-off from blue laser\n",
    "    ending_data_input =np.max(data[:,0])   #High cut-off from blue laser\n",
    "    binning_ratio_output =      (ending_data_input - starting_data_input ) /   (binnum_mem)\n",
    "\n",
    "    next_value_output = starting_data_input\n",
    "    bin_x =[starting_data_input]  #initial value for y binning\n",
    "    #generation of y binning interval. it is log scale\n",
    "    for i in range(0,binnum_mem):\n",
    "        next_value_output = next_value_output + binning_ratio_output\n",
    "        bin_x.append(next_value_output) \n",
    "\n",
    "    bin_x = np.array(bin_x)\n",
    "    bin_x[-1] = bin_x[-1] + 0.1\n",
    "\n",
    "    #binning Y\n",
    "    starting_data_input = np.min(data[:,1])    #Low cut-off from blue laser\n",
    "    ending_data_input =np.max(data[:,1])   #High cut-off from blue laser\n",
    "    binning_ratio_output =      (ending_data_input - starting_data_input ) /   (binnum_GFP)\n",
    "\n",
    "    next_value_output = starting_data_input\n",
    "    bin_y =[starting_data_input] \n",
    "\n",
    "    for i in range(0,binnum_GFP):\n",
    "        next_value_output = next_value_output + binning_ratio_output\n",
    "        bin_y.append(next_value_output) \n",
    "\n",
    "    bin_y = np.array(bin_y)\n",
    "\n",
    "    \n",
    "    x = data[:,0]\n",
    "    x = x.astype(int)\n",
    "    probability_x = []\n",
    "    temp_x_bin = np.copy(bin_x)\n",
    "    temp_x_bin = temp_x_bin[1:]\n",
    "    temp_x_bin[-1] = temp_x_bin[-1] + 1\n",
    "\n",
    "    bincount_x = np.bincount(np.digitize(x,temp_x_bin))\n",
    "    \n",
    "    for i in range(0, len(bincount_x)):\n",
    "        initial_repeat = list(itertools.repeat(p[i],bincount_x[i]))\n",
    "        probability_x.extend(initial_repeat)\n",
    "    \n",
    "    probability_x = np.array(probability_x)\n",
    "    data = np.column_stack((data,probability_x.T))\n",
    "    Pxy_origin, dummy_a, dummy_b = np.histogram2d(data[:,0], data[:,1],  bins = [bin_x, bin_y])     \n",
    "    \n",
    "    def MI(p):    \n",
    "        global entropy_list\n",
    "        global Pxy_origin\n",
    "        Pxy=np.copy(Pxy_origin)\n",
    "        Pxy = Pxy * p[:,np.newaxis]\n",
    "        Px = Pxy.sum(axis =1)\n",
    "        Px = Px/sum(Px)\n",
    "        entropy_Px = -Px*np.nan_to_num(np.log2(Px  + 0.0000001))\n",
    "        entropy_Px = entropy_Px.sum()\n",
    "        Py = Pxy.sum(axis =0)\n",
    "        Py = Py/sum(Py)\n",
    "        entropy_Py = -Py*np.nan_to_num(np.log2(Py + 0.0000001))\n",
    "        entropy_Py = entropy_Py.sum()\n",
    "        Pxy = Pxy/Pxy.sum()\n",
    "        Jentropy = 0\n",
    "        Pxy_list = Pxy.flatten()\n",
    "        Jentropy = -Pxy_list*np.nan_to_num(np.log2(Pxy_list  + 0.0000001))\n",
    "        Jentropy = Jentropy.sum()\n",
    "        MI = entropy_Px + entropy_Py - Jentropy\n",
    "        p = p/p.sum()\n",
    "        temp_entropy_list = np.array([entropy_Px, entropy_Py, Jentropy, MI])\n",
    "        entropy_list = np.vstack((entropy_list, temp_entropy_list))             \n",
    "        values = Px\n",
    "        history.append(values)\n",
    "        return -MI\n",
    "\n",
    "    \n",
    "    bin_vs_MI = []\n",
    "    for i in range(2,1):\n",
    "        value = MI(i)\n",
    "        bin_vs_MI.append(value)\n",
    "\n",
    "    history = []\n",
    "    \n",
    "    def constraint1(p):\n",
    "        return sum(p)- 1\n",
    "        \n",
    "    b = np.array([[0.0, 10.0]])\n",
    "    bnds = np.repeat(b, len(p),axis=0)\n",
    "\n",
    "    con1 = {\"type\": \"eq\", \"fun\": constraint1}\n",
    "    \n",
    "    cons = [con1]\n",
    "    minimize(MI, p,  bounds = bnds, constraints = cons,options={ 'disp': False},  method = \"SLSQP\")\n",
    "    history = np.array(history)\n",
    "    pd.history = pd.DataFrame(history)\n",
    "\n",
    "    return entropy_list[-1,3] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32fd3cb2",
   "metadata": {},
   "source": [
    "# Mutual information maximization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3b70c775",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Channel capacity:  0.43404026041191335\n"
     ]
    }
   ],
   "source": [
    "CC = channel_capacity(input_output, input_binning , output_binning) #we already defined input binning and output binning\n",
    "print(\"Channel capacity: \", CC)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18655bd0",
   "metadata": {},
   "source": [
    "# Estimation of channel capacity at infinite sample size using bootstrapping\n",
    "Note that the number of sample size of our input and output dataset is 124488.\n",
    "Using bootstrapping, we can estimate channel capacity at infinite sample size.\n",
    "Please see the detailed information in SI section 1.7.\n",
    "For short, we will do subsampling from the original data (i.e. input_output) by 50, 60, 70, 80, 90 and 100% with replacement. Then plot CC with the inverse of the sample size and extrapolate CC at infinite sample size (i.e.,0)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "955b2d51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Sample size</th>\n",
       "      <th>Trial</th>\n",
       "      <th>Channel capacity</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62242</td>\n",
       "      <td>0</td>\n",
       "      <td>0.43294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62242</td>\n",
       "      <td>1</td>\n",
       "      <td>0.427722</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62242</td>\n",
       "      <td>2</td>\n",
       "      <td>0.435864</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62242</td>\n",
       "      <td>3</td>\n",
       "      <td>0.429958</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>62242</td>\n",
       "      <td>4</td>\n",
       "      <td>0.430207</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Sample size Trial Channel capacity\n",
       "0       62242     0          0.43294\n",
       "0       62242     1         0.427722\n",
       "0       62242     2         0.435864\n",
       "0       62242     3         0.429958\n",
       "0       62242     4         0.430207"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#subsampling procedure\n",
    "data_list = pd.DataFrame(columns = [\"Sample size\", \"Trial\", \"Channel capacity\"]) #making an empty dataset\n",
    "\n",
    "subsample_percentage_list = [50, 60, 70, 80, 90, 100]\n",
    "trial_number = 10   # we will draw the samples 10 times for each subsampling percentage\n",
    "\n",
    "for subsample in subsample_percentage_list:\n",
    "    for trial in np.arange(10):\n",
    "        subsampled_input_output = np.zeros((0,2))\n",
    "        \n",
    "        for dose in np.unique(input_output[:,0]):\n",
    "            temp = input_output[input_output[:,0] == dose]\n",
    "            index_subsample = np.random.choice(len(temp), int(len(temp)*subsample/100), replace=True)\n",
    "            temp = temp[index_subsample,:]\n",
    "            subsampled_input_output = np.vstack((subsampled_input_output, temp))\n",
    "            \n",
    "        CC = channel_capacity(subsampled_input_output, input_binning , output_binning) #we already defined input binning and output binning\n",
    "        data_temp = [len(subsampled_input_output), trial, CC]\n",
    "        temp_CC_list = pd.DataFrame([data_temp], columns = [\"Sample size\", \"Trial\", \"Channel capacity\"])\n",
    "        data_list = pd.concat([data_list, temp_CC_list])\n",
    "\n",
    "data_list.head()      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d4481289",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Channel capacity at infinite sample size:  0.4331101876753283\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm0AAAJeCAYAAADm5TK9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABMYUlEQVR4nO3dd5xddZ3/8ddnMumNkgAhJCSBUAIiZQQUbNhQVpDFAqy6tsXG2gvu7m9V0FXURd3VXRf7uioq6BqkWHbFtogMSAvNSBFQpKe3mfn8/jjnhsswmZmUm3PPzOvJ4z4y59xzz3wOk9x53287kZlIkiSpvXVUXYAkSZKGZmiTJEmqAUObJElSDRjaJEmSasDQJkmSVAOGNkmSpBrorLqA7WHGjBk5b968qsuQJEka0lVXXfVAZs7sv39UhLZ58+bR3d1ddRmSJElDiog7B9pv96gkSVINtDS0RcSxEXFLRCyNiDMGOe6kiMiI6Oq3f25ErIyId5XbcyLipxFxY0QsiYi3trJ+SZKkdtGy0BYRY4DPAs8HFgGnRMSiAY6bCrwVuGKA05wDXNK03QO8MzMXAUcCbx7onJIkSSNNK1vaDgeWZuZtmbkeOA84YYDjzgLOBtY274yIFwG3A0sa+zLzT5l5dfn1CuAmYHZLqpckSWojrQxts4G7mrbvpl/AiohDgTmZeVG//VOA9wIf3NTJI2IecAgDt9BJkiSNKJVNRIiIDoruz3cO8PQHgE9m5spNvHYKcAHwtsxcvoljTouI7ojovv/++7dR1ZIkSdVo5ZIf9wBzmrb3KPc1TAUOBC6LCIDdgMURcTxwBPDiiPgYsAPQFxFrM/MzETGWIrB9PTO/u6lvnpnnAucCdHV15Ta7KkmSpAq0MrRdCSyMiPkUYe1k4NTGk5m5DJjR2I6Iy4B3ZWY38NSm/R8AVpaBLYAvAjdl5jktrF2SJKmttKx7NDN7gNOBH1JMGPh2Zi6JiDPL1rQtcRTwCuCYiLimfLxgG5UsSZLUtiJz5PccdnV1pXdEkCRJdRARV2VmV//93hFBkiSpBgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJkqQaMLRJkiTVgKFNkiSpBgxtkiRJNWBokySpzWT5n9TM0CZJUhv5AT9gIhOZzWzu5M6qy1EbMbRJktRGPsSHWMc6HuABvsN3qi5HbcTQJklSG3kpL2UiExnLWJ7BM6ouR23E0CZJUht5B+/gGq7h9/yeLrqqLmfr/Ply+MnJsPSbVVcyInRWXYAkSXqsfdin6hK2XiZc/FzYsBLu/D7scgRMW1B1VbVmS5skSWqNGNv4AjpsJ9pa/h+UJNVeb/bx0Zt+zV1rVvDBA45m1wmTqy5JEfDCy+Dmc2HOcTBlbtUV1Z6hTZJUe1+7cwn/dPPlrO/r5Z41K7jw6BdXXZIAdj4IjvpM1VWMGHaPSpJqr4Mov4qmr6WRxZY2SVLtvXzPA7hv3SruWr2C/7foKVWXI7WEoU2SVHsdEbxr3yOqLkNqKbtHJUmSasDQJkmSVAOGNkmSpBowtEmSJNWAoU2SJKkGDG2SJEk1YGiTJEmqAUObJElSDRjaJEmSasDQJkmSVAOGNkmSpBowtEmSJNWAoU2SJKkGDG2SJEk1YGiTJEmqAUObJElSDRjaJGmU6utLXvGK+9hhhzv4yEcerrocSUMwtEnSKHXjjRu44IJVLFvWxz/8w8P09WXVJUkaRGfVBUiSqjF3bieTJgVjxsCCBWPp6IiqS5I0CEObJI1S06Z1cOONc7j66nU89akTqi5H0hAMbZI0iu2yyxiOPXZS1WVIGgbHtEmSJNWAoU2SJKkGDG2SJEk1YGiTJEmqAUObJElSDRjaJEmSasDQJkmSVAOGNkmSpBowtEmSJNWAoU2SJKkGDG2SJEk1YGiTJEmqAUObJElSDRjaJEmSasDQJkmSVAMtDW0RcWxE3BIRSyPijEGOOykiMiK6+u2fGxErI+Jdm3tOSZKkkaRloS0ixgCfBZ4PLAJOiYhFAxw3FXgrcMUApzkHuGRzzylJkjTStLKl7XBgaWbelpnrgfOAEwY47izgbGBt886IeBFwO7BkC84pSZI0orQytM0G7mravrvct1FEHArMycyL+u2fArwX+ODmnrPpHKdFRHdEdN9///1bdgWSJEltorKJCBHRQdH9+c4Bnv4A8MnMXLml58/MczOzKzO7Zs6cuaWnkSRJagudLTz3PcCcpu09yn0NU4EDgcsiAmA3YHFEHA8cAbw4Ij4G7AD0RcRa4KohzilJkjQitTK0XQksjIj5FMHqZODUxpOZuQyY0diOiMuAd2VmN/DUpv0fAFZm5mcionOwc0qSJI1ULesezcwe4HTgh8BNwLczc0lEnFm2pm2zc26rmiVJktpVZGbVNbRcV1dXdnd3V12GJEnSkCLiqszs6r/fOyJIkiTVgKFNkiSpBgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJkqQaMLRJkiTVgKFNkiSpBgxtkiRJNWBokyRJqgFDmyRtgb6+vqpLkDTKGNokaTOdffbZjB07lkMOOYRVq1ZVXY6kUcLQJkmb6SMf+Qh9fX0sXbqUX/7yl1WXI2mUMLRJ0mY67rjjmDhxIuPGjePggw+uuhxJo0Rn1QVIUt187WtfY8mSJcydO5fp06dXXY6kUcLQJkmbqaOjgyc84QlVlyFplLF7VJIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1UBn1QVIkiS1s8zktttuY8WKFeyxxx7MmDGjkjpsaZMkSRrEypUrWb58Ob29vdx1112V1WFokyRJGsT48eMB6OjoYNKkSZXVYWiTpM10/vnnM3fuXE499VR6enqqLkcjzFrW8jk+x3/z31WXotK4ceM44IADWLBgAQsXLqysDkObJG2mN7zhDdx1110sXryYyy67rOpyNMK8iTfxdt7OX/FXfI/vVV2OSuPGjWP69Ol0dFQXnQxtkrSZDjjgACZPngzAXnvtVXE1Gmn+wB9Yy1r66OMe7qm6HLURZ49K0ma6+OKLufjiiznooIOYP39+1eVohPk3/o3TOI3ZzOY1vKbqctRGIjOrrqHlurq6sru7u+oyJEmShhQRV2VmV//9do9KkiTVgKFNkiSpBgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJkqQaMLRJkkaEh9ev5a7Vy6suQ2oZQ5skqfaue+Q+5vzg31h4ybn86++uqrocqSUMbZKk2vvxn+9gfV8v6/p6+dqdN1RdjtQShjZJUu2dOHsfpo0dx9jo4N37HlF1OWroWQ13/wTW3F91JSNCZ9UFSJK0tRZM2YH7jn8LPdnHuI4xVZejhsVPg0dugTHj4eSlMH6HqiuqNUObJGlE6IhgXBjY2kYmPPhbyL7isfJOQ9tWsntUkiRtexHQ9SEYOw32fCHs9ISqK6q9loa2iDg2Im6JiKURccYgx50UERkRXeX24RFxTfm4NiJObDr27RGxJCJuiIhvRsSEVl6DJEnaQoe8D169DJ59HoTtRFurZf8HI2IM8Fng+cAi4JSIWDTAcVOBtwJXNO2+AejKzIOBY4H/iIjOiJgNvKV87kBgDHByq65BkiSpXbQy9h4OLM3M2zJzPXAecMIAx50FnA2sbezIzNWZ2VNuTgCy6fhOYGJEdAKTgD+2onhJkqR20srQNhu4q2n77nLfRhFxKDAnMy/q/+KIOCIilgDXA2/IzJ7MvAf4BPAH4E/Assz8UasuQJIkqV1U1sEcER3AOcA7B3o+M6/IzAOAJwHvi4gJEbEjRWvdfGB3YHJEvHwT5z8tIrojovv++10fRpIk1VsrQ9s9wJym7T3KfQ1TgQOByyLiDuBIYHFjMkJDZt4ErCyPfTZwe2ben5kbgO8CTxnom2fmuZnZlZldM2fO3EaXJEmSVI1WhrYrgYURMT8ixlFMGFjceDIzl2XmjMycl5nzgF8Dx2dmd/maToCI2BPYD7iDolv0yIiYFBEBPAu4qYXXIEmS1BZatrhuZvZExOnADylmeX4pM5dExJlAd2YuHuTlRwNnRMQGoA94U2Y+ADwQEecDVwM9wG+Bc1t1DZIkSe0iMnPoo2quq6sru7u7qy5DkiRpSBFxVWZ29d/vSneSNMqNhg/v0khgaJOkUaq3N3ne8/5EZ+ftvOMdD1RdjqQhGNokaZS68cYN/PKXa+nrg099ajl9fba4Se2sZRMRJEntbcGCTnbcsQPo47DDxtHREVWXJGkQhjZJGqUmT+7g5pvncMstGzjooHFVlyNpCIY2SRrFpkzp4LDDxlddhqRhcEybJElSDRjaJEmSasDQJkmSVAOGNkmSpBowtEmSJNWAoU2SJKkGDG2SJEk1YGiTJEmqAUObJElSDRjaJEmSasDQJkmSVAOGNkmSpBrwhvGSJEmDyEzuvPNOVqxYwezZs9lpp50qqcOWNkmSpEGsWrWKhx9+mPXr13PnnXdWVoehTZIkaRDjxo0DoKOjgwkTJlRWh92jkiRJgxg3bhz7778/q1evZvr06ZXVYWiTJEkawoQJEyptZQO7RyVJkmrB0CZJklQDhjZJkqQaMLRJkiTVgKFNkiSpBgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJktQavevhz5fDukeqrmRE8N6jkqTa6+nr48wbf8Vdq5fzT094OrMmTqm6JAH84Jnw4LUwdgq87FYYN63qimrN0CZJqr3/vPMGPnHrb1jf28t961Zz0VNfUnVJyoQ//xroK75efhvMOLjqqmrN7lFJUu2N6xhDAB0RjB9je0RbiIBD3gcd42H2M2GnJ1RdUe35N1uSVHunzl3Ew+vXcPealbx3vyOqLkcNT/pQ8dA2YWiTJNVeRwR/u7Cr6jKklrJ7VJIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSWoz93EfK1hRdRlqM4Y2SZLayOf5PHOYw+7szg3cUHU5aiOGNkmS2sgX+SLrWc861nEJl1RdjtqIoU2SpDbyFt7CWMYyiUmcwAlVl6M2YmiTJKmNnMqpPMzD3Md97MM+VZejNuIdESRJajOTmVx1CWpDtrRJkiTVgKFNkiSpBgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJkqQaaGloi4hjI+KWiFgaEWcMctxJEZER0VVuHx4R15SPayPixKZjd4iI8yPi5oi4KSKe3MprkCTVw4oN67hv7aqqy5BapmV3RIiIMcBngecAdwNXRsTizLyx33FTgbcCVzTtvgHoysyeiJgFXBsRF2ZmD/Bp4NLMfHFEjAMmteoaJEn1sGTZ/Rz5v19jfV8f/3rwszltr4OrLkna5lrZ0nY4sDQzb8vM9cB5MOCdb88CzgbWNnZk5uoyoAFMABIgIqYDTwO+WB63PjMfadkVSJJq4dJ7b2ddby/r+3r5wu3XVl2O1BKtDG2zgbuatu8u920UEYcCczLzov4vjogjImIJcD3whjLEzQfuB74cEb+NiC9EhDdok6RR7vjd92Zy51jGRgdvWdhVdTlSS1R2w/iI6ADOAV410POZeQVwQETsD3w1Ii6hqPdQ4G8z84qI+DRwBvD/Bjj/acBpAHPnzm3JNUiS2sPCqTvx5+P/lvV9vUzpHFd1OVJLtLKl7R5gTtP2HuW+hqnAgcBlEXEHcCSwuDEZoSEzbwJWlsfeDdxdBjqA8ylC3ONk5rmZ2ZWZXTNnztwGlyNJamfjOsYY2DSitTK0XQksjIj55YSBk4HFjSczc1lmzsjMeZk5D/g1cHxmdpev6QSIiD2B/YA7MvNe4K6I2Lc8zbOAx0xskCRJGolaFtrKMWinAz8EbgK+nZlLIuLMiDh+iJcfTTFj9Brge8CbMvOB8rm/Bb4eEdcBBwP/1Ir6JWlTPvWpTzFu3DiOOOIIVq9eXXU5kkaJyMyqa2i5rq6u7O7urroMSSPEjjvuyCOPPMKUKVM4//zzed7znld1SZJGkIi4KjMfN6PGOyJI0mZ67nOfy6RJk+js7OSJT3xi1eVIGiUqmz0qSXX1jW98g2uuuYZ58+ax8847V12OpFHC0CZJm2nMmDEcdthhVZchaTvKTHp7e+nsrC46GdokSZIG0dvby4033sj69euZOXNmZeu/OqZNkiRpEKtXr6anp7i75oMPPlhZHYO2tEXEBOAvgKcCuwNrKG7mflFmLml9eZIkSdVqTDzasGEDO+20U2V1bDK0RcQHKQLbZcAVwH0UN2/fB/hoGejemZnXbYc6JUmSKjFmzBgOPPBAenp6GDt2bGV1DNbS9pvMfP8mnjsnInYBvKmnJEka8SKi0sAGg4xpy8yL+u+LiI6ImFY+f19mumKtJEka2PWfhv+cCZe9BkbBYv6tNuREhIj4RkRMi4jJFOPZboyId7e+NEmSVFuZ8Ot3wNoHYOnX4aHrq66o9oYze3RRZi4HXgRcAswHXtHKoiRJ0gjQaFzrWw/ZW2kpI8FwQtvYiBhLEdoWZ+aG1pYkSZJGhLFTiz87xsP46mZdjhTDWVz3P4A7gGuBn0fEnsCyVhYlSZJqLgKO/xncdC7MPQ6m7ll1RbU3nNB2YWb+S2MjIv4AvKZ1JUmSpBFh5yfC0Z+tuooRYzjdoxc0b2RmAue1phxJkiQNZLDFdfcDDgCmR8RfNj01jWKRXUmSJG0ng3WP7ktxR4QdgBc27V8B/E0La5IkSVI/mwxtmfl94PsR8eTMvHw71iRJkqR+BusefU9mfgw4NSJO6f98Zr6lpZVJkiRpo8G6R28q//RWVZIkSRUbrHv0wvLPrwKU9xzNzFyxnWqTJElSaTj3Hu2KiOuB64AbIuLaiDis9aVJkiSpYTiL634JeFNm/gIgIo4Gvgwc1MrCJEmS9KjhLK7b2whsAJn5S6CndSVJkiS1n+L+AtUZbPbooeWXP4uI/wC+CSTwMuCy1pcmSZJUvb6+Pm655RZWr17N7rvvzqxZsyqpY7Du0X/ut/3+pq+rjZqSJEnbyapVq1i7di0A9957b/uFtsx85vYsRJIkqR1NnDiRjo4OMpNp06ZVVsdg3aMvB76em+jAjYi9gFnlGDdJkqQRqbOzkwMPPJB169YxceLE6uoY5LmdgWsi4irgKuB+ihvF7w08HXgAOKPlFUqSJFVszJgxTJo0qdIaBuse/XREfAY4BjiKYomPNRR3SnhFZv5h+5QoSZKkQddpy8xe4MflQ5IkSRUZzjptkiRJqpihTZIkqQaGc+/RMdujEEmSJG3acFrafhcRH4+IRS2vRpKkUe4qrmJf9uUZPIOHebjqctRGhhPangjcCnwhIn4dEadFRHUry0mSNIK9nbdzK7dyOZfzVb5adTlqI0OGtsxckZmfz8ynAO+luJ3VnyLiqxGxd8srlCRpFDmEQ5jMZMYwhkXYyaVHDbrkB2wc03Yc8GpgHsU9Sb8OPBW4GNinhfVJkjSqnMM5HMMx7MquHMmRVZejNjJkaAN+B/wU+Hhm/l/T/vMj4mmtKUuStD309SUPPdTHzjt3EBFVlyPgCq7gdE5nN3bjYi5mJjOrLkltYjhj2l6Zma9tDmwRcRRAZr6lZZVJklpq/frk0EPvYbfd7uSUU+6ruhyV3sN7uJu7uY7r+Bpfq7octZHhhLZ/GWDfv27rQiRJ29fNN6/n1ls30NsL3/rWKvr6suqSBBzFURvHtB3CIVWXozayye7RiHgy8BRgZkS8o+mpaYBrt0lSze277zjmz+/kd7/bwAtfOImODrtH28FH+SjP5bnMZCYHcVDV5aiNDNbSNg6YQhHspjY9lgMvbn1pktSeLr74Yvbbbz9OO+00ent7qy5ni40fH1x77R7cccdczj9/16rLUelqruZNvIm38BbXadNjRObgzeERsWdm3rmd6mmJrq6u7O7urroMSSPEzJkzeeCBB5g8eTKLFy/mmGOOqbokjSBP42n8gl8wjnF8lI/ydt5edUlbbs398IcfwK5HwQ4uNjFcEXFVZnb13z9Y9+inMvNtwGci4nHJLjOP37YlSlI97L333qxevZrMZO7cuVWXoxHmUA7laq4mSQ7ggKrL2Tr/fSSs/hPEGDjlNpjoTNitMdiSH40pK5/YHoVIUl1ceumlfO973+OQQw5h771dY1zb1j/zzxvXaTuCI6ouZ8tlwsq7IDfAmImw5j5D21YaTvfoZGBNZvaV22OA8Zm5ejvUt03YPSpJUgVu/RpcfRbMPxEO/yi4FuCwbHb3aJP/AZ4NrCy3JwI/ophZKkmSNLB9XlE8tE0MZ522CZnZCGyUX09qXUmSJG2+9X29rOxZX3UZUssMJ7StiohDGxsRcRiwpnUlSZK0eX634iF2Xfyv7PTfn+a/7lxSdTlSSwyne/RtwHci4o9AALsBL2tlUZIkbY7Ff1zKqp4NbMg+/uV33bx8z5rPupQGMGRoy8wrI2I/YN9y1y2ZuaG1ZUmSNHzH7jafD9z4S6IveN38J1ZdjtQSw2lpgyKwLQImAIdGBJn5n60rS5Kk4Ttg+kz++BdvZk1vD7tMmFx1OVJLDBnaIuL9wDMoQtvFwPOBXwKGNklS25g6djxTx46vugypZYYzEeHFwLOAezPz1cATgektrUqSJEmPMZzQ1lhYtycipgH3AXNaW5YkSZKaDWdMW3dE7AB8HriKYpHdy1tZlCRJkh5rOLNH31R++bmIuBSYlpnXtbYsSZIkNRtO9ygR8ZcRcQ7wt8Bewz15RBwbEbdExNKIOGOQ406KiIyIrnL78Ii4pnxcGxEn9jt+TET8NiJ+MNxaJEmS6mw4s0f/Ddgb+Ga56/UR8ezMfPMQrxsDfBZ4DnA3cGVELM7MG/sdNxV4K3BF0+4bgK7M7ImIWcC1EXFhZvaUz78VuAmYNuQVSpIkjQDDaWk7BnheZn45M78MvKDcN5TDgaWZeVtmrgfOA04Y4LizgLOBtY0dmbm6KaBNALLxXETsARwHfGEYNUiSJG21NWvW8OCDD9Lb21tZDcMJbUuBuU3bc8p9Q5kN3NW0fXe5b6PynqZzMvOi/i+OiCMiYglwPfCGphD3KeA9QN8wapAkSdoq69at4+abb+bOO+/k1ltvrayO4cwenQrcFBG/oWjxOpxiRuligMw8fku+cUR0AOcArxro+cy8AjggIvYHvhoRlwDPBu7LzKsi4hlDnP804DSAuXPnDnaoJEnSRpnJqlWrWLZsGcuXL+fee+/l5ptvZocddmD//fevrK7hhLZ/3MJz38Nj13Pbo9zXMBU4ELgsIqC4Ef3iiDg+M7sbB2XmTRGxsjz2KOD4iHgBRbfptIj4r8x8ef9vnpnnAucCdHV1Zf/nJUnSyJKZrFy5kuXLlz/m0QhfjUdfXx+Zj0aDMods3NfR0cHkyZOZNm0a06ZNY+rUqeyyyy5Mnz6dPffcs5Jrg+Et+fGzLTz3lcDCiJhPEdZOBk5tOu8yYEZjOyIuA96Vmd3la+4qJyLsCewH3JGZ7wPeVx7/jPL4xwU2SZJUH/3DVv+Q1Xg0xpP1D1nN+6ZMmbIxbE2fPp1p06ax2267PSaAdXYO99br7WU4s0ePBP4V2B8YB4wBVmXmoDM3y8B1OvDD8jVfyswlEXEm0J2Ziwd5+dHAGRGxgWLs2psy84FhXZEkSdou+vr6HtONOFDr1ooVK+jrK4ahb6p1KyKYMmXKxpDVeMyaNWvjvilTptQ2bG0r0fw/cMADIropWsm+A3QBrwT2KVu9aqGrqyu7u7uHPlCSpFGgr6/vcd2I/YPXihUrHjNTcqDWrY6Ojse0bDW3bjVataZNm8aYMWO2+zXWWURclZld/fcPK7Jm5tKIGJOZvcCXI+K3lN2UkiRp+2i0bG2qC3HZsmWPadnalEbLVnPQmjp1KrNnz964PWXKFMNWmxlOaFsdEeOAayLiY8CfGOadFCRJ0vBbtgYKWxGxsXUrIja2Xk2dOnVjq1YjbDW6EQ1bI9NwQtsrKELa6cDbKWaEntTKoiRJagf9w9amBsj39fU9rvuwsd34uhG2mlu39thjD8NWTSxfvpyVK1cyY8YMxo0bV0kNwwltDwDrM3Mt8MHy9lTjW1uWJElbrq+vjxUrVgy69MOKFSvIzAFDVsNAYWvatGkbw9bUqVOZOnWqYWuEW7t2LUuXLiUzeeihhzjwwAMrqWM4oe1/KBa1XVluTwR+BDylVUVJkkan3t7eIVu2BupGbO5ChGKA/GBhqxG4Ojoc7aOhNU/IqPI2VsMJbRMysxHYyMyVETGphTVJkmqmf8tW/7DVGCA/1IoFjbDVGK/VeMyZM8ewpcpMnjyZ3XffnRUrVjBr1qzK6hhOaFsVEYdm5tUAEXEYsKa1ZUmStoclS9bwk588wFOf2sO4casGHK+1cuXKYc1GbB4Y3xizNWfOnMess2XYUl3ttttu7LbbbpXWMJzQ9jbgOxHxRyAobjf1slYWJUkaXHM34kCtWv3D1kBrbK1Zk3z60yvInMzUqdP493+fx/Tp09lhhx2YO3fuxhBm2JLaw5CL6wJExFhg33Lzlszc0NKqtjEX15XULhpha1NBqxG2htuN2LzGVv9WrsmTJw8atm65ZT2HHHIPa9YknZ2wbt18OjoePxhf0va1tYvrbgBu2OZVSVJN9Pb2bhyztalFTVeuXDnkecaMGfOYMVvTpk1jxx13ZM8999wYuIYKW9vKvvuO413vms63v72Kf/iHHQxsUpsbVktb3dnSJo1ejbC1qfW1NjdsDXSrnu0dtiSNbFvV0iZJ21tPT8+gsxH7dyMOtMYWFN2I/Rc03WmnnZg3b97G1i7DlqQ62GRoi4hDB3thYzapJDVrDltb2o0YEY8LW9OmTWPnnXdm/vz5jxmztamwJkkjzWAtbf88yHMJHLONa5FUoZ6enk12Hza3bA2l0Y3YWF+rOWw1tg1bkrT5NhnaMvOZ27MQSVtmOC1bq1atAh5/T8TmMa2dnZ2Pa9maMWMGCxYseMyYLcOWJFVjyDFt5d0P3gHMzczTImIhsG9m/qDl1Ukj2IYNGwYNW8uXL98YtgYzZsyYx43ZmjFjhi1bkjTCDGciwpeBq3j0XqP3AN8BDG0alfp3Iw4UuIYKW5lJZ2fn42Yg7rLLLuy9994btydNmmTYkiQBwwtte2XmyyLiFIDMXB3+FlENNVq2NtWFuHz5clavXj3keTo7Ox+3kKlhS6reH1Yv5+H1a3niDrtUXYrUEsMJbesjYiLF5AMiYi9gXUurkpps2LBhyAHyww1b/dfYMmxJI8NvHvojz7jsmwD83X5P5h8WPWWIV0j1M5zQ9n7gUmBORHwdOAp4VSuL0sjQHLY21bq1evXqIUPSpsLWXnvttXHbsCWNbr984G56+/pYn3384E9LDW0akYYMbZn544i4GjiS4obxb83MB1pemSqzfv36AVu2+ncjRsSA90ds7B87duzjBsjvuuuuLFy4cGPYmjhxomFL0lZ72Zz9+ezSq7lv3Wo+cMDRVZcjtcRw74gwAXi4PH5R+Uv5560rS1uif9gaqHVrON2IY8eOfdwA+V133ZV99tln47ZhS1I7mT1xKr9/wRuqLkNqqeEs+XE28DJgCdBX7k7A0LaNNLoRhztAfqA1tiJiwJatWbNmse+++xq2JEmqueG0tL2IYl02Jx9swtq1a7n55psHDF1r1qwZ8vWNsNXcutVo2WrsmzBhgmFLkqRRbDih7TZgLM4Y3aQ1a9Zwxx13MG3aNGbNmsV+++23MXwZtiRJ0rYwnNC2GrgmIv6HpuCWmW9pWVU1s+OOO/KiF72o6jIkSdIINpzQtrh8SJJKa9euZfz48bakS9puhrPkx1e3RyGSVBfvf//7+dCHPsTChQu58sormTp1atUlSRoFOoY6ICKOiogfR8StEXFbRNweEbdtj+IkqR194hOfoK+vj9///vf8/OdOpJdGg56eHlauXDng+qTby3C6R78IvJ3ipvG9rS1Hktpfb2/xVtjT08PatWsrrkZSq23YsIElS5aQmUyZMoWFCxdWUseQLW3Assy8JDPvy8wHG4+WVyZJberwww9n4sSJTJw4kUMOOaTqciS12Nq1a8lM+vr6WLlyZWV1DKel7acR8XHguzx29ujVLatKktrYJZdcwsUXX8wTnvAEFixYUHU5klps8uTJTJo0iVWrVjFr1qzK6oih+mYj4qcD7M7MPKY1JW17XV1d2d3dXXUZkiRJQ4qIqzKzq//+4cwefWZrSpIkSdJwDeuG8RFxHHAAxY3jAcjMM1tVlCS1s/Xr1/Pzn/+c/fbbjz322KPqciSNEsO5YfzngEnAM4EvAC8GftPiuiSpbZ144on87Gc/o6OjgyVLljBnzpyqS5I0Cgxn9uhTMvOVwMOZ+UHgycA+rS1LktrX5ZdfzqpVq8hMbrrppqrLkTRKDCe0rSn/XB0RuwMbgOqmTkhSxT7+8Y+zww47cNRRR/H0pz+96nIkjRLDGdP2g4jYAfg4cDWQFN2kkjQqvfa1r+W1r31t1WVIGmWGM3v0rPLLCyLiB8CEzFzW2rIkSZLUbLizR58CzGscHxFk5n+2sC5JkiQ1Gc7s0a8BewHX8Oi9RxMwtEmSJG0nw2lp6wIWZZW3tZckSRrlhjN79AZgt1YXIknS1rjkT7/nS7dfx7renqpLkVpiky1tEXEhRTfoVODGiPgNj71h/PGtL0+SpKFd+MelnPzr75PArx64my8+6QVVlyRtc4N1j35iu1UhSarEPff08JvfrOOYYyYyffpwOl/a0x9WL6cvk7V9vfx+1SNVlyO1xGCh7R5g18z8VfPOiDga+FNLq5Iktdwjj/Sy995/YN06mDWrg3vumVd1SVvsVfMO5LL7/8A9a1bwb4c+t+pypJYYLLR9CnjfAPuXlc+9sAX1SJK2k//7v7WsXVt8/cc/9tHT00dnZz1b2yZ3juM7T35R1WVILTXYv85dM/P6/jvLffNaVpEkabvo6hpPZ/nRffr0DsaMiWoLkjSowVradhjkuYnbuA5J0na2yy6dXHPNHvzsZ2t40YsmE2Fok9rZYKGtOyL+JjM/37wzIl4HXNXasiRJ28MBB4zjgAPGVV2GpGEYLLS9DfheRPwVj4a0LmAccGKL65IkSVKTTYa2zPwz8JSIeCZwYLn7osz83+1SmSRJkjYa8jZWmflT4KfboRZJkiRtQj3ndkuSNEKtYhWf5JN8k2+SeNtvPWo4N4yXJEnbyet5PedzPmPK/17KS6suSW3CljZJktrIfdzHetbTRx8P8mDV5aiNtDS0RcSxEXFLRCyNiDMGOe6kiMiI6Cq3D4+Ia8rHtRFxYrl/TkT8NCJujIglEfHWVtYvSdL29h/8B8/n+bySV/JqXl11OVvnoevh/94Kd/2w6kpGhMhsTX95RIwBbgWeA9wNXAmckpk39jtuKnARxVIip2dmd0RMAtZnZk9EzAKuBXYHZgKzMvPq8nVXAS/qf87+urq6sru7extfoSRJ2qRM+OrOsP5hGDMRXnYzTJlbdVW1EBFXZWZX//2tbGk7HFiambdl5nrgPOCEAY47CzgbWNvYkZmrM7On3JwAxUjMzPxTZl5dfr0CuAmY3bpLkCRJWyw3NL6Avp5BD9XQWhnaZgN3NW3fTb+AFRGHAnMy86L+L46IIyJiCXA98IamENd4fh5wCHDFNq5bkiRtrQh4/g9h/ovh6V+EaQuqrqj2Kps9GhEdwDnAqwZ6PjOvAA6IiP2Br0bEJZm5tnztFOAC4G2ZuXwT5z8NOA1g7lybYyVJ2u52e0rx0DbRypa2e4A5Tdt7lPsaplLcaeGyiLgDOBJY3JiM0JCZNwEry2OJiLEUge3rmfndTX3zzDw3M7sys2vmzJnb4HIkSZKq08rQdiWwMCLmR8Q44GRgcePJzFyWmTMyc15mzgN+DRxfTkSYHxGdABGxJ7AfcEdEBPBF4KbMPKeFtUuSJLWVloW2cgza6cAPKSYMfDszl0TEmRFx/BAvPxq4NiKuAb4HvCkzHwCOAl4BHNO0JMgLWnUNkiRJ7aJlS360E5f8kCRJdVHFkh+SJEnaRgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJkqQaMLRJkiTVgKFNkiSpBgxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJUhtZznLO4iy+wBdIsupy1EY6qy5AkiQ96jRO47t8l7GMZRKTOJVTqy5JbcKWNkmS2sgyltFLL0myghVVl6M2YkubJElt5PN8nnfyTmYzm1fz6qrLURsxtEmS1Eb2YA++xbeqLkNtyO5RSZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZJUe/euXcmiS7/Azt//NP97351VlyO1hKFNklR7X779em5e8SAPrV/L2377k6rLkVrC0CZJqr371q0iy68f2bCu0lqkVjG0SdIo1dubvPCFf2L8+Ns444wHqy5nqxy72wImdHQyoWMMx83aq+pytsrt3M5xHMfreT1rWVt1OWojkZlDH1VzXV1d2d3dXXUZktRWrr9+PUcccQ9r1iQR0NMzn46OqLqsLfaL++/i3rWrOHH2PnR21LdN4nk8jx/zY8YznnM4hzfyxqpL0nYWEVdlZlf//Z1VFCNJqt78+Z1Mn95BR0cfBx00rtaBDeCpM+dUXcI2MZOZjGMcHXSwMztXXY7aiKFNkkapKVM6uPnmOdx003oOPXR81eWo9B/8B4dyKLuxGy/hJVWXozZiaJOkUWz69A6OPHJC1WWoyWQm8w7eUXUZakP17fSXJEkaRQxtkiRJNWBokyRJqgFDmyRJUg0Y2iRJkmrA0CZJklQDhjZJkqQaMLRJkiTVgKFNkiSpBgxtkjRK9fYmL3jBnxg79jbe/e4Hqy5H0hAMbZI0St144wYuu2wtPT3wz/+8jL6+rLokSYPw3qOSNEotWNDJjjt20NHRx8EHj6OjI6ouSdIgDG2SNEpNntzBTTfN4eab13PwweOrLkfSEAxtkjSKTZvWweGHT6i6DEnD4Jg2SZKkIfT29rJ27Voyqxv7aUubJEnSIHp6eliyZAm9vb3ssMMOLFiwoJI6bGmTpC2wbNkyent7qy5D0nawZs0a+vr6yEyWLVtWWR2GNknaTGeccQY777wzCxcurPQNXNL2MXnyZMaPLybr7LLLLpXVYfeoJG2mz33uc/T29nL//fdz+eWXc+yxx1ZdkqQW6ujoYP/99ycz6eiorr3LljZJ2kwvf/nLGTduHNOmTePwww+vuhxJ20FEVBrYwNAmSZvtM5/5DLfffju33XYbO+20U9XlSBol7B6VpC2w++67V12CpFGmpS1tEXFsRNwSEUsj4oxBjjspIjIiusrtwyPimvJxbUScuLnnlCRJGkla1tIWEWOAzwLPAe4GroyIxZl5Y7/jpgJvBa5o2n0D0JWZPRExC7g2Ii4EcjjnlCRJGmla2dJ2OLA0M2/LzPXAecAJAxx3FnA2sLaxIzNXZ2ZPuTmBIqxtzjklSZJGlFaGttnAXU3bd5f7NoqIQ4E5mXlR/xdHxBERsQS4HnhDGeKGPGfT60+LiO6I6L7//vu37kokqcn3v/99FixYwF//9V+7wK6k7aay2aMR0QGcA7xzoOcz84rMPAB4EvC+iNisOxpn5rmZ2ZWZXTNnztz6giWp9LrXvY7bb7+dCy64gMsuu6zqciSNEq0MbfcAc5q29yj3NUwFDgQui4g7gCOBxY3JCA2ZeROwsjx2qHNKUsvtu+++TJo0icxk3rx5VZcjaZRo5ZIfVwILI2I+RbA6GTi18WRmLgNmNLYj4jLgXZnZXb7mrnIiwp7AfsAdwCODnVOStodLL72UCy+8kIMOOoi99tqr6nIkjRItC21l4Dod+CEwBvhSZi6JiDOB7sxcPMjLjwbOiIgNQB/wpsx8AGCgc7bqGiRpIFOmTOGUU06pugxJo0xk5tBH1VxXV1d2d3dXXYYkSdKQIuKqzOzqv9/bWEmSJNWAoU2SJKkGDG2SNIr19SUPPeRac1IdGNokaZRavz457LB72GWXOzn11D9XXY6kIRjaJGmUuvnm9dxyywZ6e+G881bR1zfyJ6ZJdWZok6RRat99xzFvXifjxsEJJ0yioyOqLkkjTV8vPLQEelZXXcmI0MrFdSVJbWz8+OC66/bg3nt7mT17TNXlaCS65Plw7y9hwgx46U0wdnLVFdWaLW2SNIp1dgZ77NFJhK1s2sYy4Z6fQO8aWPsALPtd1RXVnqFNkiRtexEwZW650QfT5ldaztbq6+tj3bp1VHlTAkObJEna9jJh5R+Kr2MMLL+t2nq2Qk9PDzfccANLlizhjjvuqKwOQ5skSdr2ImD/04CAHZ8AOx5QdUVbbM2aNfT29pKZPPLII5XV4UQESZLUGk/9HDz5UzBmfBHiamrSpEmMGzeOtWvXMnPmzMrqMLRJkqTW6ZxQdQVbbcyYMSxatIjMpKOjuk5KQ5skSdIQIqLyWdaOaZMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSW0kSX7Lb7mDO6ouRW3G0CZJUhv5MB/maI7mAA7gci6vuhy1kZaGtog4NiJuiYilEXHGIMedFBEZEV3l9nMi4qqIuL7885imY08p918XEZdGxIxWXoMkSdvTxVzMalazgQ38il9VXc7W+f234DtPgO4PVF3JiNCy0BYRY4DPAs8HFgGnRMSiAY6bCrwVuKJp9wPACzPzCcBfA18rj+0EPg08MzMPAq4DTm/VNUiS6uOiP/2eL9x2LWt7e6ouZat8kA8yjWnMZS6ncErV5Wy5TPjpK+HhG+Daj8HDN1ZdUe11tvDchwNLM/M2gIg4DzgB6P9TOws4G3h3Y0dm/rbp+SXAxIgYD/QBAUyOiAeBacDSll2BJKkWFv/xd5zy68Vkwq8euJsvH35c1SVtsefwHJaxrOoytl4ETJoFa/4M0QET7BjbWq3sHp0N3NW0fXe5b6OIOBSYk5kXDXKek4CrM3NdZm4A3ghcD/yRogXvi9u0aklS7dy1egV9mazp6+H21SMg8IwUL/o1PPmTxZ8Td6m6mtqrbCJCRHQA5wDvHOSYAyha4V5fbo+lCG2HALtTdI++bxOvPS0iuiOi+/7779/G1UuS2smr5z2BE2Yv5KidZ/Pvhz636nLUMGk3WPQG2OnAqisZEVrZPXoPMKdpe49yX8NU4EDgsogA2A1YHBHHZ2Z3ROwBfA94ZWb+vnzNwQCN7Yj4NjDgBIfMPBc4F6Crqyu30TVJktrQpM6xnHfkCVWXIbVUK1vargQWRsT8iBgHnAwsbjyZmcsyc0ZmzsvMecCvgUZg2wG4CDgjM5unztwDLIqImeX2c4CbWngNkiRJbaFloS0zeyhmdv6QIlh9OzOXRMSZEXH8EC8/Hdgb+MeIuKZ87JKZfwQ+CPw8Iq6jaHn7p1ZdgyRJUruIzJHfc9jV1ZXd3d1VlyFJkjSkiLgqM7v67/eOCJIkSTVgaJMkSaoBQ5skbabbbruNd7zjHVxwwQVVlyJpFGnlkh+SNCIde+yx/P73v+dzn/scCxYs4JBDDqm6JEmjgC1tkrSZ1q1bR19fHxHB+vXrqy5H0ihhS5skbaYLL7yQD3/4wzztaU/jiCOOqLocSaOEoU2SNtNBBx3Et771rarLkDTK2D0qSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SRrF/vznHi6+eDUrV/ZVXYqkIXRWXYAkqRrLlvWxaNHdrFuXLFjQyXXXzam6JEmDsKVNkkapu+/uYc2aZNWqZMmSDfT1ZdUlSRqEoU2SRqlFi8bykpdMZscdO/jIR3akoyOqLknSIOwelaRRKiL46ld3qboMScNkS5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA957VJIkaRCZyR/+8AdWrFjB7rvvzk477VRJHba0SZIkDWLVqlU89NBDrFu3jjvvvLOyOgxtkiRJgxg7diwAEcGECRMqq8PuUUmSpEGMHz+e/fbbj9WrV7PDDjtUVoehTZIkaQgTJ05k4sSJldZg96gkSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1YChTZIkqQYMbZIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGjC0SZIk1UBLQ1tEHBsRt0TE0og4Y5DjToqIjIiucvs5EXFVRFxf/nlM07HjIuLciLg1Im6OiJNaeQ2SJEntoLNVJ46IMcBngecAdwNXRsTizLyx33FTgbcCVzTtfgB4YWb+MSIOBH4IzC6f+3vgvszcJyI6gJ1adQ2SJEntopUtbYcDSzPztsxcD5wHnDDAcWcBZwNrGzsy87eZ+cdycwkwMSLGl9uvAT5SHteXmQ+06gIkSZLaRStD22zgrqbtu3m0tQyAiDgUmJOZFw1ynpOAqzNzXUTsUO47KyKujojvRMSuA70oIk6LiO6I6L7//vu3/CokSZLaQGUTEcquzXOAdw5yzAEUrXCvL3d1AnsA/5eZhwKXA58Y6LWZeW5mdmVm18yZM7dp7ZIkSdtbK0PbPcCcpu09yn0NU4EDgcsi4g7gSGBx02SEPYDvAa/MzN+Xr3kQWA18t9z+DnBoqy5AkiSpXbQytF0JLIyI+RExDjgZWNx4MjOXZeaMzJyXmfOAXwPHZ2Z32Q16EXBGZv6q6TUJXAg8o9z1LOAxExskSZJGopaFtszsAU6nmPl5E/DtzFwSEWdGxPFDvPx0YG/gHyPimvKxS/nce4EPRMR1wCsYpHtVkiRppIii8Wpki4j7gTtb/G1mUCxVImn0GCn/7kfKdYDX0q5GyrVsr+vYMzMfNyB/VIS27SEiujOzq+o6JG0/I+Xf/Ui5DvBa2tVIuZaqr8PbWEmSJNWAoU2SJKkGDG3bzrlVFyBpuxsp/+5HynWA19KuRsq1VHodjmmTJEmqAVvaJEmSasDQJkmSVAOGNknSJkVEVF2DHsufyehlaGuRiJgdEXP67fMfmqRaiIhZsPH2gWoDEbET+DMZzQxtLRARL6a4qf0F5W27ToTiH5rBTRqZIuLAiDgkIvaqupatFRHPBv4rIg4st2v/vhURnVXXsDUi4jjgwxGxZ9W16FERsWNEbLcsZWjbxiJiOsX9UN8InAjcCzw3Il4HfkKSRqKIeCHwn8CZwPsi4oSKS9piEfFc4DPAXOAIqO/7VkScGBE/jYjxmdlT1+AWEUcB/w58JzPv7Pdc7X6PR8TzIuJfI+KNEXF0ua92Hwwi4i+ALwJTt9f3rN0PuwbGAD3AI5l5D/At4H+BQ8sfsKQRJCKOBM4GXgG8FLgJeFKlRW2h8j3qk8AJwIuAv4uIJ1Za1BYq6/4QsANwZc2DWxfwhcz834iYFRHPKFveyMy+OgW3iHgq8K/ALcAk4KsRcXLdeqIi4ljgLODTmblse33f2vyg6yIzHwL+B/hQROySmQ8CPwV+BxxWaXGSWmEC8MXMXJKZa4BLgadHxMw6/TItzQBOz8xbKN6zvgs0ukjHVFnYFugDzszMQ4CrgatqHNz+APRExBTgIooPB++LiB9AEdyqLG4z7Q5ckpmfycx/Bl4H/FtEvKwuLboRMZfiA8G3MvNnEbFTRJwcEcc2hhS0St3eUNpSRPxFRPxTRHw2ImYA36D4FPGeiNg1Mx+gePN7TkTsWmmxkraJxi/+zLwM+FK5byzwILAWWFm2gsytrMhhapp08JXM/GlERGauB24D3h4R4zKzt9oqhycinhIRB2bm9cBPADLzVcCVwNVNwW1BlXUOpV+wvB94CfA24N8y802ZeTQwLSJeXkV9W+Ehig8HAGTmT4G/BM4sW63bXmb+AbiYR////xh4HkWY/mhEtKyBxtC2lcofzueAXwMTgU8Di4DLgRXAZyNiX6Dxl3FdFXVK2nYi4vnAJyLiPyNi98x8uAw6G4CHgZ7MXBMRr6D4ZTSt2oo3reyu+klEvKlp9xiAzPx3iha3v69D11VEPAf4JWWIzswHGy2Emflq4DfAL8tr/VC7/lzKrs9PlOO+9svMX1J0Kf4jRUtVwy+owe+UiNijMUEnM38M7BwR/9V4vvzg859AW0/iiYjDI+IYgMz8R2A98A6KrutXA+8Gbgb2bVUNhrattw/wo8xcnJmvoXjDeC5Fl8lXgOuATwFvAN6SmY9UU6akbaER2IAfUoxf/XJEdPbr2lkVEZ+ieEP/eGYu3/6VDtsKitaPhRHxDoCyJaoR0r4JzKcMcu2qDDofAo4Drm9MBsnM3n7BbTzFWKSz2/HnEhGHU4wr/A7FRLb3R8RbMvNLFL9HTo6Il0bEe4AXAtdWV+3Qykk6FwDfiIhPR8TUzDwW2D8ivl62TgNMBg6qrNBBRGFXius4u7wmMvMDwDvLDzeUw6HGArNaVYuhbetdAeweEU+BjZ9MbwBOBR7OzDMpmrWfn5lXV1empK1Vzg5/DfD3mXlJ+UFtFXBS+XwHxViq+cALgJdl5pKq6t2URiArw8waYCXFuK+9IuJVEbE3j7bo/B8wD9i5glKHJSLmA28H3puZl1B0UT+18XyjazcingZ0Ak/PzHYNO4cCl2bmrzLzw8B9wIkR8YYyuJ0BTAHmUPz9urXCWgcVEU+nmKRzGsW/h90pfk5k5mEUk0S+HBFfo1ht4avVVDq4LPyZYujT94DXRrG0V6N7F4CIOBl4OrC4VbV4w/gtEBEHUzRJR2beGBEfBpYD38/Mm8tj/gu4KzPfV12lkraViJiemcsi4lCKD2Z9ZYvUZ4DfZeanm459LXBVZl5TUbmDKsfa/rlp+x+AbwO7AO+lGOLxl41gExFj2nlMW9nNOT0z7yq3F1BMCHt7Zv5303EHAcv6L5vRTqJYAuPNFJNbfhIRZwK7UQy/eXs5RroWIuIvgR0z84vl9pMolsR6VWauLfcdTnF9N2bm0sqKHUREdJTjUz9O8YHgNuBlFMOg1mXmv0bEiyhacE9u5Qc1W9o2U9k1ciHFP6pvR8RJFOu0LABOKD9ZQDF2YkU1VUraliLiecCXoliRfklmrs/MnvLpGym6SYmI48tB8F9s48B2LHBnFAvoNsyn+MXZSbG8xC3AU5qeb8vZieWkgxcDT8rMu8purM7MvI2iheeIiBhbtoCSmde1Y2CLiCdGxN4RsbAcv3YdxTjC71Nc22kULbqvrbTQzfcj4HzYOLGij2L9P8p9kzPzN+XworYMbPCY2bnfAR7KzG9T/Bs5k6K1kPLDwbGtblk3tA1T+WYwBfhb4M2ZeTrwN8DHKNZk+gjFJ6EPRMR5wLsowp2kGis/iP0HcG5mPpSZ68r9jTFfY4HxEfFS4BxgdTWVDq0MbB+i6OI5tGk80ZeAV1MMBn87xZiqvaONb5sUxSLAXwEOAL4bEU8ru7EaYfom4DnAvHZeEqP8mXwTeDlFQ8DJmfkRip/HxyjWzIOidaftW9kiYv+I2B8gM1dmuYZZ+XO5h2LY0Npyks57I2JCheVuUkQ8OyL+LiI+FBGTy90dQFfZgvhSigWPn9Y0xu2elheWmT4240GRrF8OjC23DwfuoOhKANiDYnDo3Kpr9eHDx9Y/KMbjvLH8ejeKFqinNT3/SuBu4GfAoqrrHeQ6ng5cQzGT/TCK8bgzyueeQLEs0bHl9mRgWtU1D3ItMylm7L+g3H4X8AxgYb/jPgd8HQjK4UDt9KBY+uL/gGeV28+kaI36m37HvZ1iwsH+Vdc8xPW8kGKo0D8Dhw7w/ETgaxRd8NcCB1Zd8yau4ziK1s43Uoyz+z+KyYUB/AvwJ+D48tiXbs/f97a0bb57gWdR/OUjM39DsRL630XEXpl5d2ZemMU6LpLqbz1wWDnY/WLgZOBrEfHu8vlbKH5RvTEzb6yoxuGYB5yWmb/OzKsounU/GRFjs1jT7FWZeWk5fmdVtuHMyiaPUKy7tj6KxUz/jqLr8Gfx2KVLPg68K0vbv8whrQCWALfCxkHtF1D02PwVQERMpGgceGVm3lRVoUOJiEnAscC5FOvKnViO/2w830lxB4QXAK8CXpKZN1RQ6qCiWLPwzRSrPfx7Zv41xbI3C8q/Q/8LvCgzG5MNLtiev++diDBMERGNf/QR8S2KMSxvAlZn5oaI+BLF6tt3VFimpG0sIg6gWGrhz8CqzPxkue+/KYZL/AKYkMV0/7ZXjvnqKX+hvhk4KzPvaOyvur7hioiPUEyceDLFyvQfjIguinFUx2cxNqztRcRXgXEULVDHUqzzdxnF363XZrEGYEe2cRdvQ0TMyWJs4UEUKyisBxZnZnfTMR8Ezi8/KLSdclLLUzPzonJ2dVLMBj0/M7/SdNzYLNZl3K5saRtEROwbEU8ux31s/H+VmS8rtz8FvCYi3kzR9VCbNzxJw5PFwOJHgOOBWRExpdx3PjC5bJVq+8DWGIOXjx3zNZsikFKXwNY0qeB9FN1X/04RnCnDwX9RTKhoa03X8dcUt6k6kiK8fYRiVuJyiskH1CGwle6BYsIHcB7F9ZwQEbtGxHERMRP4QLsGNoCyhfl/ys2+8v/9NRRhmihuVTW+isAGNfiLXZVyoOE/UfwlvAfojoivNLoMMvOUiHgNxbozT6T4ZHd3ZQVL2mpR3L1kJ6Cb4g27FyAz/19ErKeY+faWiFhJ0U36hcqKHcJA19LoMShbbtZExFsoFj09vBzq0Zb6X0u5LzJzfRlGXxsRd1OszfZsiokUbWeg6wDIzPf2O+6VFOOjJ1O0VrWtaFoOpjlcZuY15eeE51GMKzwKeGJm3l9JoUPodx1ryz8bXZGN2eEvppgc8izg9krqtHv08cqWtf8C/iUzf1Uu63EkxT+ej2U5G6bp+PFZziiTVE8DfVADNn5QK485huJWO4soZpO25Rijwa4lHl1zqgPYkeJ+lv/Sxr9Mh/Nz+QbF2pkLKcbttd3YwiF+JmPKUD2B4o46nwROLFus2lJE7JPlwr4xyDp+EXEOxQLzLV8OY0sM5zqiWMfwVRRj9V5b5d8vu0c3bRrFGwAU0+N/QDG1/xTYeA+yxiDLtv4kJGlw5Qe1l1G8IT8L+D7FivPvjeIuCABk5v9m5ueBd7dxYBv0WhqtIZnZV3br/lMbB7bh/lxOpZiI8Pw2DWxD/UyaW3hupbhjQzsHtr8ArinDMtl0q7CmYzrKEDobOKFNA9uQ11FaCmygmKxT6d8vQ9sAyr7qc4C/jIinlm9yv6To135aOZvnKOCP5fE2V0r1N9QHtSc1fVBr27sDlDbnWtZu//I2y1DXckREHFyG0HZe0Hw4DQFPyMyb23moTblm2ekULbTro7zxexl4modc7VCG0FOzDW/hONzriGJ91v8FnpuZt1RRazND26b9gmIW0iuiWLSxNzO/QTGGbffM/GRm3lttiZK2hWF+UDuaGnxQG4XX8hSKpZja1mY0BLRli2ezzFxFcf/db1CsjzehKfA0xn49ETg1IiZsqtu0asO8joPLYx7K8hZpVXNM2yAiYkeKact/QfHJaB3wHuCYbLpvn6T6K7tyXgccBPxXZv683H8ZRbfW7yssb7N4Le1npFxHfxGxM8XabGsy8+Xlch8LgV9k5n3VVjd8dbkOZ48OIov1cT5PsQjl6ym6EV5uYJNGnixurfN1inWZ3hcR+1F8UJsJrKy0uM3ktbSfkXId/WXmgxHxeuDjEXELRQ/e09op6AxHXa7DlrZhKgcnZtZnvRxJWyAixlF0VTU+qH06M39bbVVbxmtpPyPlOvqLiLdT3J7qOdnG67ANpd2vw9AmSQMYSR/UvJb2M1KuAzYOJfo28M52nvU6lDpch6FNkiRtlXLSQbvPRB5Su1+HoU2SJKkGXPJDkiSpBgxtkiRJNWBokyRJ2gYi4ksRcV9E3DCMY18VEfdHxDXl43VDvcbQJmm7GOzNLCKOjIjPR8QzIiIj4oVNz/0gIp6xjWq4IyJmbItzDfF9Ph4RSyLi463+XoPUMG84vzgkbVNfAY7djOO/lZkHl48vDHWwoU3S9vIVNv1m9nzg0vLru4G/3x4FbY5+91UcymnAQZn57lbVI6n9lHe6eKh5X0TsFRGXRsRVEfGLcmHlLWJok7RdDPRm1uRZwE/Kr68FlkXEc/of1NxSFhFd5S2AiIgPRMRXyzfEOyPiLyPiYxFxfflmObbpNO8p9/8mIvYuXz8zIi6IiCvLx1FN5/1aRPwK+Fq/WqJsUbuhPN/Lyv2LgSnAVY19Ta95elNXyG8jYmpETImI/4mIq8vznFAeOy8ibo6Ir0TErRHx9Yh4dkT8KiJ+FxGH96vx8nL/3wzw/21MWeuVEXFdufI7ETErIn5e1nNDRDx1Uz8/SVvsXOBvM/Mwivuc/lvTcyeV/ybPj4g5Q53I21hJqlQZwjZk5rKIaOz+MHAW8OPNONVewDOBRcDlwEmZ+Z6I+B5wHPDf5XHLMvMJEfFK4FMU9xb+NPDJzPxlRMwFfgjsXx6/CDg6M9f0+35/CRwMPBGYAVwZET/PzOMjYmVmHjxAje8C3pyZv4qIKRQr4gOcmJnLy/8Xvy6DH8DewEsoblp9JcW9kI8Gjgf+DnhRedxBwJHAZOC3EXFRv+/72vK6nxQR44FfRcSPymv4YWZ+uFzsddKA/2clbZHy3/lTgO80vb+NL/+8EPhmZq4rP0h9FThmsPMZ2iRV7bnAj5p3ZObPI4KIOHozznNJZm6IiOuBMTza3Xo9MK/puG82/fnJ8utnA4ua3lSnlW+2AIsHCGxQhKdvZmYv8OeI+BnwJGDxAMc2/Ao4J4p7UH43M+8uWwH/KSKeBvQBs4Fdy+Nvb9xKJyKWAP+TmVleY/M1fb+scU1E/BQ4HLim6fnnAgdFxIvL7ekUN8O+EvhSWcN/Z2bzayRtvQ7gkYE+xGXmg02bXwA+NpyTSVKVmsezNfsw8A/99vXw6PvWhH7PrQMobwu0IR9dObyPx35AzQG+7gCObBoQPDszGzfxXjXsKxlCZn4UeB0wkaK1az/gryhuGn5Y+cb+Zx69tnVNL+9r2h7smgbaDorumcb1zc/MH5Vd1k8D7gG+UrY+StpGMnM5cHtEvAQ2Dqt4Yvn1rKZDjwduGup8hjZJlYmiaesgHtsqBEBm/gjYsXy+4Q7gsPLrk7bw276s6c/Ly69/BPxtU10HD+M8vwBeVo4Xm0kRfn4z2AsiYq/MvD4zz6Zo5dqPotXrvrKV8JnAnptzMaUTImJCROwMPKM8d7MfAm9sjO2LiH0iYnJE7An8OTM/T/FJ/9At+N6SShHxTYr3lX0j4u6IeC3FB7PXRsS1wBLghPLwt0Qxy/xa4C3Aq4Y6v92jkraL8s3sGcCMiLgbeD/FpIPfNrWK9fdh4PtN2x8EvhgRZwGXbWEpO0bEdRStVqeU+94CfLbc3wn8HHjDEOf5HvBkimtI4D2Zee8Qr3lbGcz6KN68LwGmAheWXZ7dwM2bf0lcB/yUYmzdWZn5x4iY1/T8Fyi6U68ug/L9FOPhngG8OyI2ACsBW9qkrZCZp2ziqcfNnM/M9wHv25zze+9RSZWJiH8AlmbmeVXXUlcR8QFgZWZ+oupaJLWWLW2SKpOZH6q6BkmqC1vaJEmSasCJCJIkSTVgaJMkSaoBQ5skSVINGNokSZJqwNAmSZJUA4Y2SZKkGvj/UeuZtsKZzAoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#estimating CC at infinite sample size using linear regression\n",
    "data_list[\"Sample size\"] = 1/data_list[\"Sample size\"]\n",
    "data_list = data_list.astype(float)\n",
    "fig, ax = plt.subplots(1, 1, figsize = (10,10))\n",
    "\n",
    "max_x = data_list[\"Sample size\"].max()\n",
    "\n",
    "res = stats.linregress(data_list[\"Sample size\"], data_list[\"Channel capacity\"])\n",
    "std_cc_at0 = res.intercept_stderr\n",
    "y0 = res.intercept\n",
    "slope = res.slope\n",
    "\n",
    "plt.plot(np.linspace(0,max_x), np.linspace(0,max_x) * slope + y0, color = \"k\", linewidth = 0.5)\n",
    "\n",
    "\n",
    "scatter = ax.scatter(data_list[\"Sample size\"], data_list[\"Channel capacity\"], marker=\".\", \n",
    "                     c = data_list[\"Sample size\"],  s = 20, cmap=plt.get_cmap(\"nipy_spectral\"))\n",
    "\n",
    "ax.set_xlabel(\"1/Number of samples\")\n",
    "ax.set_ylabel(\"Channel capacity (bits)\")\n",
    "ax.set_xticks(ticks = np.append(0,np.unique(data_list[\"Sample size\"])))\n",
    "ax.ticklabel_format(style='sci', axis='x', scilimits=(0,0))\n",
    "plt.xticks(rotation = 45)\n",
    "\n",
    "print(\"Channel capacity at infinite sample size: \", y0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69658fae",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}