sanity_mcmc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b835d47b-1dd0-4696-8235-42e332bee47a",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "\n",
    "import numpyro\n",
    "import numpyro.distributions as dist \n",
    "from numpyro.infer import SVI, Trace_ELBO\n",
    "from numpyro.infer import MCMC, NUTS\n",
    "from jax import random\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "numpyro.set_host_device_count(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a6797a48-6eaa-4e7b-81cd-ab1cf9f00f8f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cpu\n"
     ]
    }
   ],
   "source": [
    "from jax.lib import xla_bridge\n",
    "print(xla_bridge.get_backend().platform)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "66db2b3c-820d-41b5-a48b-97b443731f12",
   "metadata": {},
   "outputs": [],
   "source": [
    "zeisel = pd.read_csv(\"Zeisel_UMI_counts.txt\", sep = \"\\t\")\n",
    "zeisel = zeisel.set_index(\"GeneID\").T\n",
    "cnt_mat = zeisel.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "65e85258-0a2c-4058-87a6-77a2d4b2ca38",
   "metadata": {},
   "outputs": [],
   "source": [
    "gene = \"Zbed3\"\n",
    "\n",
    "n_g = zeisel[gene].values\n",
    "N_c = cnt_mat.sum(axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c7a3f3d0-f940-447d-b82d-a9335551ca4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGxCAYAAACeKZf2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA0UklEQVR4nO3de3xTZZ7H8W9aaFqUFhB7ASoXQW6WclFqK0JxkA6i0pnZlUEEBlEXB1gqOg6gctGRegNckQFnHWVFEHBWLgsIW6qIhYJyHVBEUATUtniBFAoUaM7+kW0hkJQk5NbD5/165VXynOc853eenCRfTk5ai2EYhgAAAEwiItQFAAAA+BPhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmArhBgAAmEqtUBcQbHa7XT/88IPq1q0ri8US6nIAAIAHDMPQsWPH1KhRI0VEVH9u5ooLNz/88IOSk5NDXQYAAPDBoUOH1KRJk2r7XHHhpm7dupIckxMbGxviagAAgCdKS0uVnJxc9T5enSsu3FR+FBUbG0u4AQCghvHkkhIuKAYAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKYS0nAza9YsdejQoeq3Baenp+uDDz6odp333ntPbdq0UXR0tFJSUrRy5cogVQsAAGqCkIabJk2a6Pnnn9eWLVu0efNm3X777erXr58+//xzl/03bNigAQMGaNiwYdq2bZuys7OVnZ2tXbt2BblyAAAQriyGYRihLuJ8DRo00EsvvaRhw4ZdtKx///4qKyvT8uXLq9puueUWdezYUbNnz/Zo/NLSUsXFxclms/G3pQAAqCG8ef8Omz+cWVFRoffee09lZWVKT0932aewsFBjxoxxasvKytKSJUuCUCEAmNDZs9KUKVJBgdStmzR+vFSrluvl6emSxSJt2OC6LxAmQn5U7ty5U+np6Tp16pSuvvpqLV68WO3atXPZt7i4WAkJCU5tCQkJKi4udjt+eXm5ysvLq+6Xlpb6p3AAMIMpU6RJkyTDkNascbRNmOB6eV7euXZXfYEwEfJvS7Vu3Vrbt2/Xpk2b9Mgjj2jIkCH64osv/DZ+bm6u4uLiqm7Jycl+GxsAaryCAkdwkRw/CwrcLz+fq75AmAh5uImKilLLli3VpUsX5ebmKjU1Vf/xH//hsm9iYqJKSkqc2kpKSpSYmOh2/HHjxslms1XdDh065Nf6AaBG69bN8VGT5PjZrZv75edz1RcIEyH/WOpCdrvd6WOk86Wnpys/P185OTlVbXl5eW6v0ZEkq9Uqq9Xq7zIBwBzGj3f8PP+aG3fLXV1zA4ShkIabcePGqU+fPrruuut07NgxzZ8/X2vXrtXq1aslSYMHD1bjxo2Vm5srSRo9erR69OihqVOnqm/fvlqwYIE2b96sv/3tb6HcDQCouWrVqv66mUstB8JQSMPN4cOHNXjwYBUVFSkuLk4dOnTQ6tWrdccdd0iSDh48qIiIc5+cZWRkaP78+Xrqqac0fvx4tWrVSkuWLNGNN94Yql0AAABhJux+z02g8XtuAACoebx5/w75BcUAAAD+RLgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmQrgBAACmEtJwk5ubq5tvvll169ZVfHy8srOztWfPnmrXmTNnjiwWi9MtOjo6SBUDAIBwF9Jw8/HHH2vEiBHauHGj8vLydObMGfXu3VtlZWXVrhcbG6uioqKq24EDB4JUMQAACHe1QrnxVatWOd2fM2eO4uPjtWXLFnXv3t3tehaLRYmJiYEuDwAA1EBhdc2NzWaTJDVo0KDafsePH1fTpk2VnJysfv366fPPPw9GeQAAoAYIm3Bjt9uVk5OjW2+9VTfeeKPbfq1bt9abb76ppUuX6p133pHdbldGRoa+++47l/3Ly8tVWlrqdAMAAOZlMQzDCHURkvTII4/ogw8+UEFBgZo0aeLxemfOnFHbtm01YMAAPfvssxctnzRpkiZPnnxRu81mU2xs7GXVDAAAgqO0tFRxcXEevX+HxZmbkSNHavny5froo4+8CjaSVLt2bXXq1En79u1zuXzcuHGy2WxVt0OHDvmjZAAAEKZCekGxYRgaNWqUFi9erLVr16p58+Zej1FRUaGdO3fqzjvvdLncarXKarVebqkAAKCGCGm4GTFihObPn6+lS5eqbt26Ki4uliTFxcUpJiZGkjR48GA1btxYubm5kqRnnnlGt9xyi1q2bKmjR4/qpZde0oEDB/Tggw+GbD8AAED4CGm4mTVrliQpMzPTqf2tt97SH/7wB0nSwYMHFRFx7tOzI0eO6KGHHlJxcbHq16+vLl26aMOGDWrXrl2wygYAAGEsbC4oDhZvLkgCAADhocZdUAwAAOAvhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqIQ03ubm5uvnmm1W3bl3Fx8crOztbe/bsueR67733ntq0aaPo6GilpKRo5cqVQagWAADUBCENNx9//LFGjBihjRs3Ki8vT2fOnFHv3r1VVlbmdp0NGzZowIABGjZsmLZt26bs7GxlZ2dr165dQawcAACEK4thGEaoi6j0448/Kj4+Xh9//LG6d+/usk///v1VVlam5cuXV7Xdcsst6tixo2bPnn3JbZSWliouLk42m02xsbF+qx0AAASON+/fYXXNjc1mkyQ1aNDAbZ/CwkL16tXLqS0rK0uFhYUBrQ2XcPas9MwzUu/ejp9nz4a6Iu8Es/6aPldmVvnY9Ool3X67dMcdwTkeJk6Urr/ecZs0yf/bC7djLhzq4TlvbkaYqKioMPr27Wvceuut1farXbu2MX/+fKe2mTNnGvHx8S77nzp1yrDZbFW3Q4cOGZIMm83mt9phGMbkyYZhsRiG5Pg5eXKoK/JOMOuv6XNlZuc/NpW3YBwP529P8v/2wu2YC4d6eM7XODabzeP377A5czNixAjt2rVLCxYs8Ou4ubm5iouLq7olJyf7dXz8v4ICx8uy5PhZUBDaerwVzPpr+lyZ2fmPTaVgHA+etF3uNsLpmAuHenjOm1pYhJuRI0dq+fLl+uijj9SkSZNq+yYmJqqkpMSpraSkRImJiS77jxs3Tjabrep26NAhv9WN83TrJlksjn9bLI77NUkw66/pc2Vm5z82lYJxPHjSdrnbCKdjLhzq4TlvarVCuXHDMDRq1CgtXrxYa9euVfPmzS+5Tnp6uvLz85WTk1PVlpeXp/T0dJf9rVarrFarv0qGO+PHO34WFDieuJX3a4pg1l/T58rMKh+Ldesku12KjJRuuy3wx0NFhfTOO477gwb5f3vhdsyFQz08500tpN+W+uMf/6j58+dr6dKlat26dVV7XFycYmJiJEmDBw9W48aNlZubK8nxVfAePXro+eefV9++fbVgwQJNmTJFW7du1Y033njJbfJtKQAAap4a822pWbNmyWazKTMzU0lJSVW3hQsXVvU5ePCgioqKqu5nZGRo/vz5+tvf/qbU1FT94x//0JIlSzwKNgAAwPzC6vfcBANnbgAAqHlqzJkbAAAAfyPcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAUyHcAAAAU/FLuCktLdWSJUu0e/dufwwHAADgM5/Czb333qvXXntNknTy5EnddNNNuvfee9WhQwf993//t18LBAAA8IZP4WbdunW67bbbJEmLFy+WYRg6evSoXn31Vf3lL3/xa4EAAADe8Cnc2Gw2NWjQQJK0atUq/e53v1OdOnXUt29f7d27168FAgAAeMOncJOcnKzCwkKVlZVp1apV6t27tyTpyJEjio6O9muBAAAA3qjly0o5OTkaOHCgrr76ajVt2lSZmZmSHB9XpaSk+LM+AAAAr/gUbv74xz8qLS1NBw8e1B133KGICMcJoBYtWui5557za4EAAADe8OljqWeeeUZt27bVb37zG1199dVV7bfffrvWrFnjt+IAAAC8ZTEMw/B2pcjISBUVFSk+Pt6p/eeff1Z8fLwqKir8VqC/lZaWKi4uTjabTbGxsaEuBwAAeMCb92+fztwYhiGLxXJR+44dO6q+RQUAABAKXl1zU79+fVksFlksFt1www1OAaeiokLHjx/X8OHD/V4kAACAp7wKN6+88ooMw9ADDzygyZMnKy4urmpZVFSUmjVrpvT0dL8XCQAA4Cmvws2QIUMkSc2bN1dGRoZq164dkKIAAAB85dNXwXv06CG73a6vvvpKhw8flt1ud1revXt3vxQHAADgLZ/CzcaNG3XffffpwIEDuvDLVhaLJay/LQUAAMzNp3AzfPhw3XTTTVqxYoWSkpJcfnMKAAAgFHwKN3v37tU//vEPtWzZ0t/1AAAAXBaffs9NWlqa9u3b5+9aAAAALptPZ25GjRqlxx57TMXFxUpJSbnoW1MdOnTwS3EAAADe8unPL1T+oUyngSyWqt9cHM4XFPPnFwAAqHm8ef/26czN/v37fSoMAAAg0HwKN02bNvV3HQAAAH7hU7h5++23q10+ePBgn4oBAAC4XD5dc1O/fn2n+2fOnNGJEycUFRWlOnXq6JdffvFbgf7GNTcAANQ83rx/+/RV8CNHjjjdjh8/rj179qhbt2569913fSoaAADAH3wKN660atVKzz//vEaPHu2vIQEAALzmt3AjSbVq1dIPP/zgzyEBAAC84tMFxcuWLXO6bxiGioqK9Nprr+nWW2/1S2EAAAC+8CncZGdnO923WCy69tprdfvtt2vq1Kn+qAsAAMAnPoUbu93u7zoAAAD84rKvuTEMQz58mxwAACAgfA43b7/9tlJSUhQTE6OYmBh16NBBc+fO9WdtAAAAXvMp3EybNk2PPPKI7rzzTi1atEiLFi3Sr3/9aw0fPlzTp0/3eJx169bp7rvvVqNGjWSxWLRkyZJq+69du1YWi+WiW3FxsS+7AQAATMina25mzJihWbNmOf2ZhXvuuUft27fXpEmT9Oijj3o0TllZmVJTU/XAAw/ot7/9rcfb37Nnj9NvJ4yPj/e8eAAAYGo+hZuioiJlZGRc1J6RkaGioiKPx+nTp4/69Onj9fbj4+NVr149r9cDAADm59PHUi1bttSiRYsual+4cKFatWp12UVdSseOHZWUlKQ77rhD69evD/j2AABAzeHTmZvJkyerf//+WrduXdUv7Vu/fr3y8/Ndhh5/SUpK0uzZs3XTTTepvLxcb7zxhjIzM7Vp0yZ17tzZ5Trl5eUqLy+vul9aWhqw+gAAQOj5FG5+97vfadOmTZo+fXrVRcBt27bVp59+qk6dOvmzPietW7dW69atq+5nZGTo66+/1vTp091+Uys3N1eTJ08OWE0AACC8+BRuJKlLly565513/FmLT7p27aqCggK3y8eNG6cxY8ZU3S8tLVVycnIwSgMAACHgU7hZuXKlIiMjlZWV5dS+evVq2e12ny4S9tX27duVlJTkdrnVapXVag1aPQAAILR8uqB47NixqqiouKjdMAyNHTvW43GOHz+u7du3a/v27ZKk/fv3a/v27Tp48KAkx1mX879u/sorr2jp0qXat2+fdu3apZycHH344YcaMWKEL7sBAABMyKczN3v37lW7du0uam/Tpo327dvn8TibN29Wz549q+5Xfnw0ZMgQzZkzR0VFRVVBR5JOnz6txx57TN9//73q1KmjDh06aM2aNU5jAACAK5vF8OEPQyUmJmr+/Pm6/fbbndrXrFmj++67T4cPH/Zbgf5WWlqquLg42Ww2p18ECAAAwpc3798+fSzVr18/5eTk6Ouvv65q27dvnx577DHdc889vgwJAADgFz6FmxdffFFXXXWV2rRpo+bNm6t58+Zq27atrrnmGr388sv+rhEAAMBjPl1zExcXpw0bNigvL087duyo+qvg3bt393d9AAAAXvHpmhtPpaSkaOXKlWH1e2W45gYAgJon4NfceOrbb7/VmTNnArkJAAAAJwENNwAAAMFGuAEAAKZCuAEAAKZCuAEAAKZCuAEAAKbiU7g5deqUR/1ef/11JSQk+LIJAAAAn/j0S/zq1aunrl27qkePHsrMzFRGRoZiYmIu6nffffdddoEAAADe8OnMzZo1a/TrX/9amzZtUr9+/VS/fn1169ZNTz75pPLy8vxdIwAAgMcu+zcUnz17Vp999plef/11zZs3T3a7XRUVFf6qz+/4DcUAANQ83rx/+/SxlCR99dVXWrt2bdWtvLxcd911lzIzM30dEgAA4LL5FG4aN26skydPKjMzU5mZmfrzn/+sDh06yGKx+Ls+AAAAr/h0zc21116rEydOqLi4WMXFxSopKdHJkyf9XRsAAIDXfAo327dvV3FxscaOHavy8nKNHz9eDRs2VEZGhp588kl/1wgAAOCxy76g+Oeff9batWu1dOlSvfvuu1xQDAAA/C7gFxS///77VRcSf/HFF2rQoIG6deumqVOnqkePHj4VDQAA4A8+nbmJj49X9+7dlZmZqR49eiglJSUQtQUEZ24AAKh5An7m5vDhwz4VBgAAEGg+/56bSqdOndLp06ed2jgjAgAAQsWnb0uVlZVp5MiRio+P11VXXaX69es73QAAAELFp3DzxBNP6MMPP9SsWbNktVr1xhtvaPLkyWrUqJHefvttf9cIAADgMZ8+lvqf//kfvf3228rMzNTQoUN12223qWXLlmratKnmzZungQMH+rtOAAAAj/h05uaXX35RixYtJDmur/nll18kSd26ddO6dev8Vx0AAICXfAo3LVq00P79+yVJbdq00aJFiyQ5zujUq1fPb8UBAAB4y6dwM3ToUO3YsUOSNHbsWM2cOVPR0dF69NFH9ac//cmvBQIAAHjD62tuzpw5o+XLl2v27NmSpF69eunLL7/Uli1b1LJlS3Xo0MHvRQIAAHjK63BTu3Zt/fOf/3Rqa9q0qZo2beq3ogAAAHzl08dS999/v/7+97/7uxYAAIDL5tNXwc+ePas333xTa9asUZcuXXTVVVc5LZ82bZpfigMAAPCWT+Fm165d6ty5syTpq6++clpmsVguvyoAAAAf+RRuPvroI3/XAQAA4Bc+XXMDAAAQrgg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVAg3AADAVEIabtatW6e7775bjRo1ksVi0ZIlSy65ztq1a9W5c2dZrVa1bNlSc+bMCXidAACg5ghpuCkrK1NqaqpmzpzpUf/9+/erb9++6tmzp7Zv366cnBw9+OCDWr16dYArBQAANYVPfzjTX/r06aM+ffp43H/27Nlq3ry5pk6dKklq27atCgoKNH36dGVlZQWqTAAAUIOENNx4q7CwUL169XJqy8rKUk5OTmgKcuXsWWnKFKmgQOrWTRo/XqpVy3n5s89K77zjuD9okPTUU44+7tb1tt3TWrzZn08+kSoqpIgI6dZbJYtF2rDBMe4TT0gvvuj9di6sz5NxXO2T5Nw2Zox0zz3Sjh1SSorUvbu0cWP12zh/3IwMyTCkwkLv6rjUPrvb33XrJLtdiox0zK3dLs2bJx05ItWr5zhGnn764v10te/uav3LX6S5c8+NOXiwNHas9Pzz59rtdunUKUf/unWlESMcj/e8eY5xrrvOUWP37s7bOXtWeuYZ6bXXpJMnpYQE6frrpR49Lp7vyvsffyx984109KijnhYtpNtuc8x7QYH09ddSSYlktUpxcdLhw45tRUdL9eufqz8311H/0aOO9vvvlx59VEpNlQ4edKxTr540cqQ0caLj/jPPSDNmSMePS1dffW4/586Viosd+xMZ6fhptzv2+/77HcfDN984+khSUpJ0332O58Znn0lRUY72sjLHuB07OubAMBzPFXfPn8rHcfJkaeZMxxxarY7HIibGcQzv2iWdOOHoV17uWD821rGNWrUcj8mYMdJdd0mffupYPzXVsT27XTpwwHnfoqMdY1ksju2UlzvqrFXL8bN2bccY5eXSmTOO42HUKEetU6acm79ataT4eMdYhw87+kqObUqOGuPipGPHHI/Dddc5ajIM6dtvHY+R3e7o8+23jnlz9RypPF5sNsc+L1vmaJ8xwzF2rVqOuWrQwPFYPf206+e13e547a08Xs5//ZUcc/HrXzvm0DAcj2nEeR9wVB57Tz3lGPvOOx2vM6mp0sqVjnHOf31PTj73+FQ+zpXPRUkaOPDcsVB5fJz//KruNdKT1yl/8sf7STAZYUKSsXjx4mr7tGrVypgyZYpT24oVKwxJxokTJ1yuc+rUKcNms1XdDh06ZEgybDabv0p3NnmyYVgshiE5fk6efPFyxyF57lbZx9263rZ7Wou3++PqZrEYRs+evm3nwvo8GcfVPl3Y1qyZ97W6209v6ric/b3UzdV+umtztV1XY/bs6dm2LzUn7sZ3Nd/e7POlbu7qr1ev+jn0x7b9dTv/cbzccdwd9/68+XrMeHKrV8/z58il9vVSz2tXfQ3D8/2bPPnivj17Vv9c8PRxru413t18+Pr67g1/vJ9cJpvNZnj6/h3Gscs/cnNzNXny5OBtsKDAcbhJjp8FBRcvd7VOdet62+5pLd7ujyuG4fifiy/bubA+T8Zxt0/nt333nfe1uttPb+uoTnX768m6letduE1P5syVHTs82/aFLtyOu/Fdzbc3+3wp7uo/etR1uy/Hf6D5+rx0NY67496ffD1mPHH0qOfPkUvt66We1676Sp7vX0HBxX137HB/JsObx7m613h38+Gv46g6/ng/CaIa9VXwxMRElZSUOLWVlJQoNjZWMTExLtcZN26cbDZb1e3QoUOBLbJbN8dpRsnxs1u3i5e7Wqe6db1t97QWb/fHFYvFcUrWl+1cWJ8n47japwvbmjTxvlZ3++lNHZezv56s68m+V3d8XSg11bNtX+jC7bgb39V8e7PPl+Ku/nr1XLdXzlc4Of9xvNxx3B33/uTrMeOJevU8f45cal8v9bx21VfyfP+6dbu4b2pq9c8FTx/n6l7j3c2Hr6/v3vDH+0kQ1agzN+np6Vq5cqVTW15entLT092uY7VaZbVaA13aOZWfq154XcT5yysqnK+5qezjbl1v2z2txZv98faaG2/Gru6am0utc34fX6+5uXBcV59le1OHt/vryTU3rvbzUm3nb9du9/81N+ePf/Zs+F9zU1nz2bPhec3NmTNX9jU3+/c75s3Vc8Tba27cPa9dXXNz/rG8cqVn19yMH++ozdU1N+e/vru65qbyuShVf83NhfV7cs1NIPnj/SSILIbhr3PE3jt+/Lj27dsnSerUqZOmTZumnj17qkGDBrruuus0btw4ff/993r77bclOb4KfuONN2rEiBF64IEH9OGHH+rf//3ftWLFCo+/LVVaWqq4uDjZbDbFxsYGbN8AAID/ePP+HdKPpTZv3qxOnTqpU6dOkqQxY8aoU6dOmjBhgiSpqKhIByv/9yWpefPmWrFihfLy8pSamqqpU6fqjTfe4GvgAACgSkjP3IQCZ24AAKh5asyZGwAAAH8j3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMh3AAAAFMJi3Azc+ZMNWvWTNHR0UpLS9Onn37qtu+cOXNksVicbtHR0UGsFgAAhLOQh5uFCxdqzJgxmjhxorZu3arU1FRlZWXp8OHDbteJjY1VUVFR1e3AgQNBrBgAAISzkIebadOm6aGHHtLQoUPVrl07zZ49W3Xq1NGbb77pdh2LxaLExMSqW0JCQhArBgAA4Syk4eb06dPasmWLevXqVdUWERGhXr16qbCw0O16x48fV9OmTZWcnKx+/frp888/D0a5AACgBghpuPnpp59UUVFx0ZmXhIQEFRcXu1yndevWevPNN7V06VK98847stvtysjI0Hfffeeyf3l5uUpLS51uAADAvEL+sZS30tPTNXjwYHXs2FE9evTQ+++/r2uvvVavv/66y/65ubmKi4uruiUnJwe5YgAAEEwhDTcNGzZUZGSkSkpKnNpLSkqUmJjo0Ri1a9dWp06dtG/fPpfLx40bJ5vNVnU7dOjQZdcNAADCV0jDTVRUlLp06aL8/PyqNrvdrvz8fKWnp3s0RkVFhXbu3KmkpCSXy61Wq2JjY51uAADAvGqFuoAxY8ZoyJAhuummm9S1a1e98sorKisr09ChQyVJgwcPVuPGjZWbmytJeuaZZ3TLLbeoZcuWOnr0qF566SUdOHBADz74YCh3AwAAhImQh5v+/fvrxx9/1IQJE1RcXKyOHTtq1apVVRcZHzx4UBER504wHTlyRA899JCKi4tVv359denSRRs2bFC7du1CtQsAACCMWAzDMEJdRDCVlpYqLi5ONpuNj6gAAKghvHn/rnHflgIAAKgO4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJhKWISbmTNnqlmzZoqOjlZaWpo+/fTTavu/9957atOmjaKjo5WSkqKVK1cGqVIAABDuQh5uFi5cqDFjxmjixInaunWrUlNTlZWVpcOHD7vsv2HDBg0YMEDDhg3Ttm3blJ2drezsbO3atSvIlQMAgHBkMQzDCGUBaWlpuvnmm/Xaa69Jkux2u5KTkzVq1CiNHTv2ov79+/dXWVmZli9fXtV2yy23qGPHjpo9e/Ylt1daWqq4uDjZbDbFxsb6b0cAAEDAePP+XStINbl0+vRpbdmyRePGjatqi4iIUK9evVRYWOhyncLCQo0ZM8apLSsrS0uWLAlkqZd2/Lh0443SgQMXL6tbV0pNlT77TCovd14WGysdOyZdmDEtlovbvFkOAPBeRITjtdXVa3Lt2pLdLp09e669dm3HOhe+tlssUmSkVKuWo8/p0859IiIcy86ccWzLYnG0V243NlYaMULatEn6+muppMTRbrWe61u/vjRokDRunPTii9LHH5/rGxMjjRwpPfGEdM890o4djvehZcukadOkggKpWzdp/HjHWFOmONoyMhw1FBaeW17Li6hw9uy5sXxZ309CGm5++uknVVRUKCEhwak9ISFBX375pct1iouLXfYvLi522b+8vFzl5x1QpaWll1m1GykproON5AgvBQWul7mr51LBhWADAP5nt7tuNwxHQLnQmTPu+58967idOuV6O+ePd+FremmplJt78Xrnj3X0qPTMM9Inn0hr1zqPceqU9Oyz0ty50rffOto++ujce5VhSGvWnOs/aZKjLS/vXFvl8gkTXO+jK1OmnBvLl/X9JOTX3ARabm6u4uLiqm7JycmB2dB33wVmXAAAqrNjh/v/8F743vTdd+f6GobjP94FBa7Xr1zujfPH8mV9PwlpuGnYsKEiIyNVUnm67f+VlJQoMTHR5TqJiYle9R83bpxsNlvV7dChQ/4p/kJNmgRmXAAAqpOaeu6jqgtd+N7UpMm5vhaL46Ojbt1cr1+53Bvnj+XL+n4S0o+loqKi1KVLF+Xn5ys7O1uS44Li/Px8jRw50uU66enpys/PV05OTlVbXl6e0tPTXfa3Wq2yWq3+Lv1iO3dyzQ0A1HRX6jU3kvtrbrxR2d/V+EEU8m9LLVy4UEOGDNHrr7+url276pVXXtGiRYv05ZdfKiEhQYMHD1bjxo2V+/+fPW7YsEE9evTQ888/r759+2rBggWaMmWKtm7dqhtvvPGS2+PbUgAA1Dw15ttSkuOr3T/++KMmTJig4uJidezYUatWraq6aPjgwYOKiDj36VlGRobmz5+vp556SuPHj1erVq20ZMkSj4INAAAwv5CfuQk2ztwAAFDzePP+bfpvSwEAgCsL4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJgK4QYAAJhKyP+2VLBV/rWJ0tLSEFcCAAA8Vfm+7clfjbriws2xY8ckScnJySGuBAAAeOvYsWOKi4urts8V94cz7Xa7fvjhB9WtW1cWiyXU5XiltLRUycnJOnToEH/08zzMi3vMjWvMi2vMi3vMjWvBnBfDMHTs2DE1atRIERHVX1VzxZ25iYiIUJMmTUJdxmWJjY3lyeUC8+Iec+Ma8+Ia8+Iec+NasOblUmdsKnFBMQAAMBXCDQAAMBXCTQ1itVo1ceJEWa3WUJcSVpgX95gb15gX15gX95gb18J1Xq64C4oBAIC5ceYGAACYCuEGAACYCuEGAACYCuEmTFRUVOjpp59W8+bNFRMTo+uvv17PPvvsJX/N9Nq1a9W5c2dZrVa1bNlSc+bMCU7BQeLLvKxdu1YWi+WiW3FxcRArD45jx44pJydHTZs2VUxMjDIyMvTZZ59Vu47ZjxnJ+3kx4zGzbt063X333WrUqJEsFouWLFnitNwwDE2YMEFJSUmKiYlRr169tHfv3kuOO3PmTDVr1kzR0dFKS0vTp59+GqA9CJxAzM2kSZMuOn7atGkTwL3wv0vNy/vvv6/evXvrmmuukcVi0fbt2z0a97333lObNm0UHR2tlJQUrVy50v/FX4BwEyZeeOEFzZo1S6+99pp2796tF154QS+++KJmzJjhdp39+/erb9++6tmzp7Zv366cnBw9+OCDWr16dRArDyxf5qXSnj17VFRUVHWLj48PQsXB9eCDDyovL09z587Vzp071bt3b/Xq1Uvff/+9y/5XwjEjeT8vlcx0zJSVlSk1NVUzZ850ufzFF1/Uq6++qtmzZ2vTpk266qqrlJWVpVOnTrkdc+HChRozZowmTpyorVu3KjU1VVlZWTp8+HCgdiMgAjE3ktS+fXun46egoCAQ5QfMpealrKxM3bp10wsvvODxmBs2bNCAAQM0bNgwbdu2TdnZ2crOztauXbv8VbZrBsJC3759jQceeMCp7be//a0xcOBAt+s88cQTRvv27Z3a+vfvb2RlZQWkxlDwZV4++ugjQ5Jx5MiRAFcXWidOnDAiIyON5cuXO7V37tzZePLJJ12ucyUcM77Mi9mPGUnG4sWLq+7b7XYjMTHReOmll6rajh49alitVuPdd991O07Xrl2NESNGVN2vqKgwGjVqZOTm5gak7mDw19xMnDjRSE1NDWClwXXhvJxv//79hiRj27Ztlxzn3nvvNfr27evUlpaWZvzbv/2bH6p0jzM3YSIjI0P5+fn66quvJEk7duxQQUGB+vTp43adwsJC9erVy6ktKytLhYWFAa01mHyZl0odO3ZUUlKS7rjjDq1fvz7QpQbd2bNnVVFRoejoaKf2mJgYt/9jvBKOGV/mpZLZj5lK+/fvV3FxsdOxEBcXp7S0NLfHwunTp7VlyxandSIiItSrVy9THT++zE2lvXv3qlGjRmrRooUGDhyogwcPBrrcsBeq15wr7m9LhauxY8eqtLRUbdq0UWRkpCoqKvTcc89p4MCBbtcpLi5WQkKCU1tCQoJKS0t18uRJxcTEBLrsgPNlXpKSkjR79mzddNNNKi8v1xtvvKHMzExt2rRJnTt3DmL1gVW3bl2lp6fr2WefVdu2bZWQkKB3331XhYWFatmypct1roRjxpd5uVKOmUqV1xK5OhbcXWf0008/qaKiwuU6X375ZWAKDQFf5kaS0tLSNGfOHLVu3VpFRUWaPHmybrvtNu3atUt169YNaM3hzN1rTqCvZyPchIlFixZp3rx5mj9/vtq3b191PUSjRo00ZMiQUJcXMr7MS+vWrdW6deuq+xkZGfr66681ffp0zZ07N1ilB8XcuXP1wAMPqHHjxoqMjFTnzp01YMAAbdmyJdSlhZS383IlHTMIjPPPJnfo0EFpaWlq2rSpFi1apGHDhoWwsisTH0uFiT/96U8aO3asfv/73yslJUWDBg3So48+qtzcXLfrJCYmqqSkxKmtpKREsbGxpvgfuOTbvLjStWtX7du3L0BVhs7111+vjz/+WMePH9ehQ4f06aef6syZM2rRooXL/lfCMSN5Py+umPWYkRzHgSSXx0Llsgs1bNhQkZGRXq1TE/kyN67Uq1dPN9xwg2mPIU+5e80J9DFDuAkTJ06cUESE88MRGRkpu93udp309HTl5+c7teXl5Sk9PT0gNYaCL/Piyvbt25WUlOTP0sLKVVddpaSkJB05ckSrV69Wv379XPa7Eo6Z83k6L66Y+Zhp3ry5EhMTnY6F0tJSbdq0ye2xEBUVpS5dujitY7fblZ+fb6rjx5e5ceX48eP6+uuvTXsMeSpkrzkBvVwZHhsyZIjRuHFjY/ny5cb+/fuN999/32jYsKHxxBNPVPUZO3asMWjQoKr733zzjVGnTh3jT3/6k7F7925j5syZRmRkpLFq1apQ7EJA+DIv06dPN5YsWWLs3bvX2LlzpzF69GgjIiLCWLNmTSh2IaBWrVplfPDBB8Y333xj/O///q+RmppqpKWlGadPnzYM48o8ZgzD+3kx4zFz7NgxY9u2bca2bdsMSca0adOMbdu2GQcOHDAMwzCef/55o169esbSpUuNf/7zn0a/fv2M5s2bGydPnqwa4/bbbzdmzJhRdX/BggWG1Wo15syZY3zxxRfGww8/bNSrV88oLi4O+v5djkDMzWOPPWasXbvW2L9/v7F+/XqjV69eRsOGDY3Dhw8Hff98dal5+fnnn41t27YZK1asMCQZCxYsMLZt22YUFRVVjTFo0CBj7NixVffXr19v1KpVy3j55ZeN3bt3GxMnTjRq165t7Ny5M6D7QrgJE6Wlpcbo0aON6667zoiOjjZatGhhPPnkk0Z5eXlVnyFDhhg9evRwWu+jjz4yOnbsaERFRRktWrQw3nrrreAWHmC+zMsLL7xgXH/99UZ0dLTRoEEDIzMz0/jwww9DUH3gLVy40GjRooURFRVlJCYmGiNGjDCOHj1atfxKPGYMw/t5MeMxU/n19gtvQ4YMMQzD8ZXnp59+2khISDCsVqvxq1/9ytizZ4/TGE2bNjUmTpzo1DZjxgzjuuuuM6KiooyuXbsaGzduDNIe+U8g5qZ///5GUlKSERUVZTRu3Njo37+/sW/fviDu1eW71Ly89dZbLpefPw89evSo6l9p0aJFxg033GBERUUZ7du3N1asWBHwfeGvggMAAFPhmhsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAGAqhBsAAZeZmamcnBy/j9u9e3fNnz/f7+O6c/r0aTVr1kybN28O2jYBeI9wA6BGWrZsmUpKSvT73/++qq1Zs2ayWCzauHGjU9+cnBxlZmZe9jajoqL0+OOP689//vNljwUgcAg3AGqkV199VUOHDr3or8ZHR0cHNHwMHDhQBQUF+vzzzwO2DQCXh3ADIKiOHDmiwYMHq379+qpTp4769OmjvXv3OvX5z//8TyUnJ6tOnTr6zW9+o2nTpqlevXpVy3/88Ud9+OGHuvvuuy8a/+GHH9bGjRu1cuVKn2t888031b59e1mtViUlJWnkyJFVy+rXr69bb71VCxYs8Hl8AIFFuAEQVH/4wx+0efNmLVu2TIWFhTIMQ3feeafOnDkjSVq/fr2GDx+u0aNHa/v27brjjjv03HPPOY1RUFCgOnXqqG3btheN37x5cw0fPlzjxo2T3W73ur5Zs2ZpxIgRevjhh7Vz504tW7ZMLVu2dOrTtWtXffLJJ16PDSA4aoW6AABXjr1792rZsmVav369MjIyJEnz5s1TcnKylixZon/913/VjBkz1KdPHz3++OOSpBtuuEEbNmzQ8uXLq8Y5cOCAEhISLvpIqtJTTz2lt956S/PmzdOgQYO8qvEvf/mLHnvsMY0ePbqq7eabb3bq06hRIx04cMCrcQEED2duAATN7t27VatWLaWlpVW1XXPNNWrdurV2794tSdqzZ4+6du3qtN6F90+ePKno6Gi327n22mv1+OOPa8KECTp9+rTH9R0+fFg//PCDfvWrX1XbLyYmRidOnPB4XADBRbgBUOM0bNhQR44cqbbPmDFjdPLkSf31r3/1eNyYmBiP+v3yyy+69tprPR4XQHARbgAETdu2bXX27Flt2rSpqu3nn3/Wnj171K5dO0lS69at9dlnnzmtd+H9Tp06qbi4uNqAc/XVV+vpp5/Wc889p2PHjnlUX926ddWsWTPl5+dX22/Xrl3q1KmTR2MCCD7CDYCgadWqlfr166eHHnpIBQUF2rFjh+6//341btxY/fr1kySNGjVKK1eu1LRp07R37169/vrr+uCDD2SxWKrG6dSpkxo2bKj169dXu72HH35YcXFxXv2iv0mTJmnq1Kl69dVXtXfvXm3dulUzZsxw6vPJJ5+od+/eXuw5gGAi3AAIqrfeektdunTRXXfdpfT0dBmGoZUrV6p27dqSpFtvvVWzZ8/WtGnTlJqaqlWrVunRRx91usYmMjJSQ4cO1bx586rdVu3atfXss8/q1KlTHtc3ZMgQvfLKK/rrX/+q9u3b66677nL6qnphYaFsNpv+5V/+xcs9BxAsFsMwjFAXAQDVeeihh/Tll186ff26uLhY7du319atW9W0adOg1dK/f3+lpqZq/PjxQdsmAO9w5gZA2Hn55Ze1Y8cO7du3TzNmzNB//dd/aciQIU59EhMT9fe//10HDx4MWl2nT59WSkqKHn300aBtE4D3OHMDIOzce++9Wrt2rY4dO6YWLVpo1KhRGj58uF/Gvvrqq90u++CDD3Tbbbf5ZTsAQodwA+CKsm/fPrfLGjdu7PHXwQGEL8INAAAwFa65AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApkK4AQAApvJ/7OUZqC3JgyYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log(N_c), n_g, c = \"r\", s = 5)\n",
    "plt.ylabel(\"raw_counts\")\n",
    "plt.xlabel(\"log(N_c)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "885f48b3-8e40-459d-b7e2-f2f3978cca80",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def sanity_model(n_g, N_c):\n",
    "    \"\"\"\n",
    "    n_g - is the vector of counts of gene g across all cells is the observed data\n",
    "    N_c - is the total count of gene transcripts in a given cell.\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    α_g = numpyro.sample(\"α_g\", dist.Uniform(0, 1))\n",
    "    v_gi = numpyro.sample(\"v_gi\", dist.Uniform(jnp.log(0.001), jnp.log(50)))\n",
    "    v_g = jnp.exp(v_gi)\n",
    "    cell_plate = numpyro.plate(\"cell_plate\", size = n_g.shape[0])\n",
    "    \n",
    "    with cell_plate:\n",
    "        δ_gc = numpyro.sample(\"δ_gc\", dist.Normal(loc = 0., scale = v_g))\n",
    "        fold_change = numpyro.deterministic(\"fold_change\", jnp.exp(δ_gc))\n",
    "        \n",
    "        α_gc = numpyro.deterministic(\"α_gc\", α_g*fold_change)\n",
    "        rate = numpyro.deterministic (\"rate\", N_c*α_gc)\n",
    "        numpyro.sample(\"n_cg\", dist.Poisson(rate = rate), obs=n_g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "23f04d84-2f37-4e35-a14f-2937610d706b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"228pt\" height=\"443pt\"\n",
       " viewBox=\"0.00 0.00 228.30 443.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 439)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-439 224.3,-439 224.3,4 -4,4\"/>\n",
       "<g id=\"clust1\" class=\"cluster\">\n",
       "<title>cluster_cell_plate</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"64.3,-8 64.3,-371 212.3,-371 212.3,-8 64.3,-8\"/>\n",
       "<text text-anchor=\"middle\" x=\"170.3\" y=\"-15.8\" font-family=\"Times,serif\" font-size=\"14.00\">cell_plate</text>\n",
       "</g>\n",
       "<!-- α_g -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>α_g</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" cx=\"27.3\" cy=\"-273\" rx=\"27.1\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"27.3\" y=\"-269.3\" font-family=\"Times,serif\" font-size=\"14.00\">α_g</text>\n",
       "</g>\n",
       "<!-- α_gc -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>α_gc</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" stroke-dasharray=\"5,2\" cx=\"121.3\" cy=\"-201\" rx=\"32.49\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"121.3\" y=\"-197.3\" font-family=\"Times,serif\" font-size=\"14.00\">α_gc</text>\n",
       "</g>\n",
       "<!-- α_g&#45;&gt;α_gc -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>α_g&#45;&gt;α_gc</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M44.98,-258.83C58.89,-248.47 78.5,-233.87 94.37,-222.06\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"96.8,-224.61 102.73,-215.83 92.62,-218.99 96.8,-224.61\"/>\n",
       "</g>\n",
       "<!-- v_gi -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>v_gi</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" cx=\"138.3\" cy=\"-417\" rx=\"29.5\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"138.3\" y=\"-413.3\" font-family=\"Times,serif\" font-size=\"14.00\">v_gi</text>\n",
       "</g>\n",
       "<!-- δ_gc -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>δ_gc</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" cx=\"138.3\" cy=\"-345\" rx=\"31.4\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"138.3\" y=\"-341.3\" font-family=\"Times,serif\" font-size=\"14.00\">δ_gc</text>\n",
       "</g>\n",
       "<!-- v_gi&#45;&gt;δ_gc -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>v_gi&#45;&gt;δ_gc</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M138.3,-398.7C138.3,-390.98 138.3,-381.71 138.3,-373.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"141.8,-373.1 138.3,-363.1 134.8,-373.1 141.8,-373.1\"/>\n",
       "</g>\n",
       "<!-- fold_change -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>fold_change</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" stroke-dasharray=\"5,2\" cx=\"138.3\" cy=\"-273\" rx=\"65.79\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"138.3\" y=\"-269.3\" font-family=\"Times,serif\" font-size=\"14.00\">fold_change</text>\n",
       "</g>\n",
       "<!-- δ_gc&#45;&gt;fold_change -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>δ_gc&#45;&gt;fold_change</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M138.3,-326.7C138.3,-318.98 138.3,-309.71 138.3,-301.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"141.8,-301.1 138.3,-291.1 134.8,-301.1 141.8,-301.1\"/>\n",
       "</g>\n",
       "<!-- fold_change&#45;&gt;α_gc -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>fold_change&#45;&gt;α_gc</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M134.1,-254.7C132.2,-246.9 129.92,-237.51 127.81,-228.83\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"131.21,-228 125.45,-219.1 124.41,-229.65 131.21,-228\"/>\n",
       "</g>\n",
       "<!-- rate -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>rate</title>\n",
       "<ellipse fill=\"white\" stroke=\"black\" stroke-dasharray=\"5,2\" cx=\"121.3\" cy=\"-129\" rx=\"29.8\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"121.3\" y=\"-125.3\" font-family=\"Times,serif\" font-size=\"14.00\">rate</text>\n",
       "</g>\n",
       "<!-- α_gc&#45;&gt;rate -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>α_gc&#45;&gt;rate</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M121.3,-182.7C121.3,-174.98 121.3,-165.71 121.3,-157.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"124.8,-157.1 121.3,-147.1 117.8,-157.1 124.8,-157.1\"/>\n",
       "</g>\n",
       "<!-- n_cg -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>n_cg</title>\n",
       "<ellipse fill=\"grey\" stroke=\"black\" cx=\"121.3\" cy=\"-57\" rx=\"31.7\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"121.3\" y=\"-53.3\" font-family=\"Times,serif\" font-size=\"14.00\">n_cg</text>\n",
       "</g>\n",
       "<!-- rate&#45;&gt;n_cg -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>rate&#45;&gt;n_cg</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M121.3,-110.7C121.3,-102.98 121.3,-93.71 121.3,-85.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"124.8,-85.1 121.3,-75.1 117.8,-85.1 124.8,-85.1\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7f233418e760>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpyro.render_model(sanity_model, model_args=(n_g,N_c,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1a109061-7415-4649-aeda-456ab267375f",
   "metadata": {},
   "outputs": [],
   "source": [
    "zsan = pd.read_csv(\"Zeisel_Sanity_normalization.txt\", sep = \"\\t\")\n",
    "zsan = zsan.set_index(\"GeneID\").T\n",
    "cnt_mat_san = zsan.to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1aacbf58-4053-4e5d-8dcb-8bdfb428bfd5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f2240466280>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5H0lEQVR4nO3de3hU9YH/8U8uZJIASRYJNwlgoFyLQBIuSS9ci5c+Vrf+EK2rIDeXR9eibNvwVKGsy4/VXtxq+Wn3qQUqWlHWdtvaYqkgFhKIRIY7gZDQUCASxSSEhFzP74/vziQTZjKZkEkyJ+/X85xncma+5+Q7h4H58L2dMMuyLAEAANhIeGdXAAAAoL0RcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO1EdnYFOkNDQ4MuXLig3r17KywsrLOrAwAAWsGyLF25ckWDBg1SeHjLbTTdMuBcuHBBSUlJnV0NAADQBufOndPgwYNbLNMtA07v3r0lmQsUFxfXybUBAACtUV5erqSkJPf3eEu6ZcBxdUvFxcURcAAACDGtGV7CIGMAAGA7BBwAAGA7BBwAAGA7QQ0469atU0ZGhmJjY5WQkOC1zBNPPKHU1FQ5HA5NnDixVeedMWOGwsLCPLZ//ud/br+KAwCAkBbUgFNTU6N58+Zp+fLlLZZbtGiR5s+fH9C5ly5dqosXL7q3559//kaqCgAAbCSos6jWrl0rSdq0aZPPMi+++KIkqaSkRIcPH271uWNjYzVgwIAbqh8AALCnkB2D8/rrr6tv37764he/qFWrVqmysrKzqwQAALqIkFwH51vf+paGDh2qQYMG6fDhw/re976nvLw8vfPOO17LV1dXq7q62r1fXl7eUVUFAACdIOCAk5mZqeeee67FMidOnNDo0aPbXCl/li1b5v55/PjxGjhwoGbPnq0zZ85o+PDh15Vfv369u7sMAADYX8ABZ+XKlVq4cGGLZZKTk9tanzaZOnWqJCk/P99rwFm1apWeeuop975rqWcAAGBPAQecxMREJSYmBqMubeZ0OiVJAwcO9Pq6w+GQw+HowBoBAIDOFNRBxkVFRXI6nSoqKlJ9fb2cTqecTqcqKircZfLz8+V0OlVcXKyqqip3mZqaGknS+fPnNXr0aOXk5EiSzpw5o2effVa5ubk6e/asfve73+nhhx/WV7/6Vd16663BfDsAYF91ddKpU+YRsIGgDjJevXq1Nm/e7N6fNGmSJGnXrl2aMWOGJGnJkiXavXv3dWUKCws1bNgw1dbWKi8vzz1LKioqSn/5y1/0n//5n7p69aqSkpJ077336umnnw7mWwEA+6qrk9LTpQMHpLQ0KTtbigzJOSiAW5hlWVZnV6KjlZeXKz4+XmVlZdxNHABOnZJGjWrcz8uTRo7svPoAPgTy/R2y6+AAANpJcrJpuZGkyZPNvjd0YyGEEHAAoLuLjDTdUnl5UlaW9+4pVzfWqFHm8do1wg66NAIOAMCEmpEjfY+9KSgwY3Qk8zh5cmPYIeSgCyLgAAD8a9qNNX68dPSo+fnAARN+gC6GgAMA8K9pN1ZOTuvG7ACdiHmAAIDWcXVjSSbsFBSYcMOUcnRBfCoBAIFrGnaALoguKgAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDsEHAAAYDtBCzjr1q1TRkaGYmNjlZCQ4LXME088odTUVDkcDk2cOLHV587OztasWbPUs2dPxcXF6atf/aqqqqrap+IAACDkBS3g1NTUaN68eVq+fHmL5RYtWqT58+e3+rzZ2dm6/fbbNXfuXOXk5Oijjz7S448/rvBwGqMAAIARGawTr127VpK0adMmn2VefPFFSVJJSYkOHz7cqvM++eSTeuKJJ5SZmel+btSoUW2vKAAAsJ2Qava4dOmS9u/fr379+ikjI0P9+/fX9OnTtWfPnhaPq66uVnl5uccGAADsK6QCTkFBgSTpBz/4gZYuXart27crJSVFs2fP1unTp30et379esXHx7u3pKSkjqoyAADoBAEFnMzMTIWFhbW4nTx5Mlh1VUNDgyTp0Ucf1SOPPKJJkybphRde0KhRo/TLX/7S53GrVq1SWVmZezt37lzQ6ggAADpfQGNwVq5cqYULF7ZYJjk5+Ubq06KBAwdKksaOHevx/JgxY1RUVOTzOIfDIYfDEbR6AQCAriWggJOYmKjExMRg1cWvYcOGadCgQcrLy/N4/tSpU7rjjjs6qVYAAKCrCdosqqKiIl2+fFlFRUWqr6+X0+mUJI0YMUK9evWSJOXn56uiokLFxcWqqqpylxk7dqyioqJ0/vx5zZ49W7/61a80ZcoUhYWF6Tvf+Y7WrFmjCRMmaOLEidq8ebNOnjypbdu2BeutAACAEBO0gLN69Wpt3rzZvT9p0iRJ0q5duzRjxgxJ0pIlS7R79+7ryhQWFmrYsGGqra1VXl6eKisr3WVWrFiha9eu6cknn9Tly5c1YcIE7dixQ8OHDw/WWwEAACEmzLIsq7Mr0dHKy8sVHx+vsrIyxcXFdXZ1AABAKwTy/R1S08QBAABag4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4ADAABsh4CD9lVXJ506ZR4BAOgkBBy0n7o6KT1dGjXKPBJyAACdhICD9lNQIB04YH4+cMDshxpaoADAFgg4aD/JyVJamvl58mSzH0pogQIA2yDgoP1ERkrZ2VJenpSVZfZDSWe0QNFiBABBQcBB+4qMlEaODL1wI3V8CxQtRgAQNAQcwKWjW6DsMGYJALooAg7QVEe2QIX6mCUA6MJCsB8BsAlXi1FBgQk3oditBwBdFP+iAp3J1WIEAGhXdFEBAADbIeAAAADbIeAAAADbIeAAAADbIeAAAADbIeAAQHOuW2hcu8atNIAQFdSAs27dOmVkZCg2NlYJCQleyzzxxBNKTU2Vw+HQxIkT/Z7z7NmzCgsL87q9/fbb7fsGAHQ/TW+hcdNN3EoDCFFBDTg1NTWaN2+eli9f3mK5RYsWaf78+a06Z1JSki5evOixrV27Vr169dIdd9zRHtUG0J01vYVGZaV55FYaQMgJ6kJ/a9eulSRt2rTJZ5kXX3xRklRSUqLDhw/7PWdERIQGDBjg8dxvfvMb3XffferVq1fbKwsAUuMtNA4ckGJjTcjhVhpAyAn5lYxzc3PldDq1YcMGn2Wqq6tVXV3t3i8vL++IqgEIRU1voTFkiFRU1DG30nCN+5E67n5ogI2F/CDjV199VWPGjFFGRobPMuvXr1d8fLx7S0pK6sAaAgg5rltoREd3TNioq5OmTZPGjTPbtGmM+QFuUMABJzMz0+cgX9d28uTJYNT1OlVVVXrjjTe0ePHiFsutWrVKZWVl7u3cuXMdUj8AaJWCAik3t3E/N5cxP8ANCvi/JStXrtTChQtbLJPcQX3V27ZtU2VlpR5++OEWyzkcDjkcjg6pEwAELDlZSk1tDDlpaYz5AW5QwAEnMTFRiYmJwahLwF599VV94xvf6DL1AYA2iYyU9u0L/hicujrTMtQRY4qAThbUMThFRUVyOp0qKipSfX29nE6nnE6nKioq3GXy8/PldDpVXFysqqoqd5mamhpJ0vnz5zV69Gjl5OR4nDs/P18ffvihlixZEsy3AAAdIzJSGjvWbMEKN671fVjXB91AUCP86tWrtXnzZvf+pEmTJEm7du3SjBkzJElLlizR7t27rytTWFioYcOGqba2Vnl5eap0rUfxv375y19q8ODBmjt3bjDfAgDYQ9P1fVzr+owc2bl1AoIozLIsq7Mr0dHKy8sVHx+vsrIyxcXFdXZ1ACD4XC04Bw6YdX2ysuimQsgJ5PubTzcAdAdN1/dhDA66AT7hANBduNb3AbqBkF/oDwAAoDkCDgAAsB0CDgCg87juwcW0dbQzAg4AoHOwNg+CiIADAOgc3tbmAdoJAQcA0DmSk819tySzNg/330I7Ypo4AKBzsDYPgohPEwCg87A2D4KELioAAGA7BBwAAGA7BBwAQPfFOjy2RcABAHRPrMNjawQcAED3xDo8tkbAAQB0T6zDY2tMEwcAdE+sw2Nr/GkCALov1uGxLbqoAACA7RBwAADoCExJ71AEHAAAgo0p6R2OgAMAQLAxJb3DEXAAAAg2pqR3OGZRAQAQbExJ73BcYQAAOgJT0jsUXVQAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2CDgAAMB2ghZw1q1bp4yMDMXGxiohIcFrmSeeeEKpqalyOByaOHFiq85bXFyshx56SAMGDFDPnj2VkpKi//7v/26/igMAgJAXtIBTU1OjefPmafny5S2WW7RokebPn9/q8z788MPKy8vT7373Ox05ckTf/OY3dd999+ngwYM3WmUAAGATkcE68dq1ayVJmzZt8lnmxRdflCSVlJTo8OHDrTpvVlaWXn75ZU2ZMkWS9PTTT+uFF15Qbm6uJk2adGOVBgAAthByY3AyMjK0detWXb58WQ0NDXrzzTd17do1zZgxw+cx1dXVKi8v99gAAIB9hVzAeeutt1RbW6ubbrpJDodDjz76qH7zm99oxIgRPo9Zv3694uPj3VtSUlIH1hgAAHS0gAJOZmamwsLCWtxOnjwZrLpKkp555hmVlpbqL3/5iw4cOKCnnnpK9913n44cOeLzmFWrVqmsrMy9nTt3Lqh1BAAAnSugMTgrV67UwoULWyyTnJx8I/Vp0ZkzZ/Szn/1MR48e1bhx4yRJEyZM0F//+ldt2LBBr7zyitfjHA6HHA5H0OoFAAC6loACTmJiohITE4NVF78qKyslSeHhng1PERERamho6IwqAQCALihoY3CKiorkdDpVVFSk+vp6OZ1OOZ1OVVRUuMvk5+fL6XSquLhYVVVV7jI1NTWSpPPnz2v06NHKycmRJI0ePVojRozQo48+qpycHJ05c0Y//vGPtWPHDt1zzz3BeisAACDEBG2a+OrVq7V582b3vmsK965du9wznpYsWaLdu3dfV6awsFDDhg1TbW2t8vLy3C03PXr00B//+EdlZmbqrrvuUkVFhUaMGKHNmzfrzjvvDNZbAQAAISbMsiyrsyvR0crLyxUfH6+ysjLFxcV1dnUAAEArBPL9HXLTxAEAAPwh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4AAAANsh4ARLXZ106pR59PX68eNm81UGAAC0CQEnGOrqpPR0adQo89g8wNTVSdOmSePGmW3aNM8yvsJRoM+39vXWvqdTp6Rr1278XAAABBkBJxgKCqQDB8zPBw6Y/eav5+Y27ufmNpbxFY4Cfd7F3+ut0fQcN93UeK7mYaetQar5ca05T3uENgCAbRFwgiE5WUpLMz9Pnmz2m7+emtq4n5bWWMZXOAr0eRd/r7dG03NUVjaea/Jkz7DTliDVPIC15jwthb2moefaNem998xja0NU0+cJUQAQuqwg+fd//3crPT3diomJseLj472W+Zd/+RcrJSXFioqKsiZMmNCq8+bn51v33HOP1bdvX6t3797WvHnzrOLi4oDqVlZWZkmyysrKAjouILW1lpWXZx59vX7smNmalqmttay0NMuSLGvy5MbXAn3e3/kCfS+uc8TGmsfx482ja9u+3XM/L691587LC/w8zY9xXWdXHdPSLOvKlca6xsRYVmpq42tVVZ5lvV3L1FTPY1r6c2zpz7k1x7TlHADQDQXy/R20FpyamhrNmzdPy5cvb7HcokWLNH/+/Fad8+rVq5o7d67CwsK0c+dO7d27VzU1NbrrrrvU0NDQHtVuP5GR0siR5tHX62PHmq1pmchIKTtbysuTsrIaXwv0eX/nC/S9uM7x2WfmMSfHs5Vq+vSWW618ad7a1ZrzeGsha95StXVrY2tTVVVjl+CBA9Lu3f5bw3JzPY/x1vLVlu6/llqspk1rHHTuGoR++PD1A9G9tSzR2gQAnoKdtjZu3OizBcdlzZo1rWrBee+996zw8HCP5FZaWmqFhYVZO3bsaHWdOqQFpztor5aItpzH2zFNW6p8teBMnuzZguOrNSwtzfMYb3Xx1pLkj78WK1frUUrK9c/V1l5fR1/P+bpmx45Z1qFDni2HTZ8/fNg8un5u3sIIAJ0okO/vNvx3vvNUV1crLCxMDofD/Vx0dLTCw8O1Z88ezZkzx+dx1dXV7v3y8vKg17VbcLVS+doP5nm8HZOdbVpakpPN/mefmdaa6dPNftPXmpf1dg7p+jJNuVqSXOORWtNq1fwYV4uVq+VI8hyA3vQ5VyuSt9an5s81v36umXtNz52aKu3ZI335y95/p0tKirRxo7kGyclSUZHva1JRYVrP7r1XunTJdzkACLKQGmQ8bdo09ezZU9/73vdUWVmpq1ev6l//9V9VX1+vixcv+jxu/fr1io+Pd29JSUkdWGt0mObdgtHR0m23mcfmr/nqQmz6fGu6GQPt/mt+THS02T92rHHgeVqaCRVNuQaie+ue8zeoXbp+5p5k9nfvbjncSNLHH0sTJpglDZrOomveHVZRIcXFSUuWSP/wD57lfA3yrqgwA8ErKhjcDaBdBRRwMjMzFRYW1uJ28uTJYNVViYmJevvtt/X73/9evXr1Unx8vEpLS5WSkqLwcN9vZdWqVSorK3Nv586dC1od0c34C0GtOcY1HmvfPhN8srOl/ftN6Dl0yDxmZzeGruahqjVBq/nMPcmEounTr3++JU1n0TUfl7R1q+lQa+rAARNWfI07iouTbr/dPLrGIU2bZn6eOlX6/e+lI0dMAPrjH6WDB68fl8SimQC8CKjteOXKlVq4cGGLZZJbO7i0jebOnaszZ87o008/VWRkpBISEjRgwIAWf6/D4fDo1gK6pObdbmPHtq6cr+eav75vX2PLSNOQ1fT5sLDGkOL6+ZFHTCuOJMXGmpDjraVo/nxp6VLPkDN5snls2oXWdJC3q6zrsWlr0scfS9/4hu/35K2Lbdw4c+7oaBOKXntN6tdP+sIXPAf0u0LR3/8uzZplygOwlYACTmJiohITE4NVl4D07dtXkrRz505dunRJ32jpH0IAjS1FrX3eZf9+E4Cklsfg9OollZdfPwZH8j3uyBWiXI9paebRX7eZ5L2L7dgx8zv++lfTTdZUSop5L5JpHXKFtuhoqaTEPLpW6z5/Xrr5ZvMe6+ulixcJQkCICdrov6KiIl2+fFlFRUWqr6+X0+mUJI0YMUK9evWSJOXn56uiokLFxcWqqqpylxk7dqyioqJ0/vx5zZ49W7/61a80ZcoUSdLGjRs1ZswYJSYmKjs7W9/+9rf15JNPatSoUcF6K0D31jwAtdRS1KuXtHix+TkhofH55oO6XfuDBkl790pf+pJ04UJjIDp+3HRVVVX5/l1Nu9iahpyjR6UNG64v//HHjd1qrnAjmUCTkiL17u35fHPR0WZ5hC98QTp9Wjp7Vho2TBozprHOtAgBXUewpnItWLDAknTdtmvXLneZ6dOney1TWFhoWZZlFRYWXnfM9773Pat///5Wjx49rC984QvWj3/8Y6uhoSGgujFNHAgBVVWW9e67Zrr6lSvm548/9j7N/dAhyxo3rnFa/+efe59+75pW33wafiCba/kB15aSYlmTJjXux8RY1ssvmzq4ljJwLWtQVcWijsANCOT7O8yymo8KtL/y8nLFx8errKxMcXFxnV0dAO2hrs6zlcjfGJzDh6X77zetMWlpUkNDyy04beUaxJ2b2ziGKS3NdKM17+5rOjboq1+VCgulc+doFQL+VyDf3wQcAg7QfTUNRZLvMTh//7u0dq1nSHFJSTFtNwcPBva7v/hF052Wlma67CTPsUFNB3xHR5tAdOmSNGWKtG2b1LevmW02ZgxrDaHbIOD4QcABEDBXGBoyxPcYnNOnpXnzPGeSNR047QpH48eb6e8ueXnmsS1jCVNSTEuVa1bctWvSzp1SUpJpuWppYUYgxBBw/CDgAAgab6s5S43hqKjIPH7lK40zy7KyTBlfLTitkZYm7dpluuRcg7NdgSo1VfrlLxtbo1ytU2FhUkRE4Gs5AZ2EgOMHAQdAp2s+Zsj1XPMxOGfPSmvWtK4L7Be/MCtJB2rMGOn//l8TfP7+dxN8hgyRevQgAKFLIeD4QcABEFLq6qQTJ8yA46ZjcIYPN4sr5uaalqCdO7234Nyo4cOlt982XWuRkSyUiE5DwPGDgAPANpq3BDUfg1NQID30UPvMEHOtHv2lL3kulLh3r2nl+fBDE7AuXpSGDvWcuQa0AwKOHwQcAN2K6wam1dWeY3Dq66V/+ifTGtNa27eb+4e1xpAhpnXp0qXGW4T4uyM90AICjh8EHAD4X67ur7NnpQEDPMfgNDSY1pqaGlPWWwtOoJoOfN64UbpyxbQ4Pf645+rXgBcEHD8IOADQSk27vFxr7rgWSpw/X8rPb7/f9fnnJuR4G4ANKLDv7/AOqhMAIBRFR0t33tk4wFgyjykppuXn2DGzXbkivfuu9NFH0ptvmlYgb2Jjff+uDRtMuElPN2sCpaebfVcXW0WF9N57JnQBftCCQwsOALS/pusBXbjgOQanoEB6+OHr7xr/+edmvE7TBQ+PHZMWLPC8+3xsrPTZZ+Z8ru61wYPNPtPabS2Q728+AQCA9ufrzvKSmV21b59plSkt9RyD06uXWbTQtQiiZH6WGhc+rKw0x6xefX1Iksz4nn37PNcXosur26EFhxYcAOhamt8jLD39+hac7GxpwgTf58jLMy05ri6vAwcab3JaUNDYokRrT0hhDA4AIHQ1DR6RkSbM5OVJ5eVmmvpnn5lWINed2ptLS/O8RYarBcgVcsaNM+Fo3Dhp2jQTdmA7tODQggMAoanpFHdfY3CatuA0v8mpi6u1x8U1c2zgQHM+idaeLoIxOAAA+4uMNKFl/PiWy2RnN97s9Mtf9hy307S1RzLh5qabvN/iYtw46Y03WKE5RPAnBACwN1fri9Q4uNnXGJzdu33fv+vYMdO11XQQ87Vr5pjp07knVxdDwAEAdB+RkaYFxpfp0/3fpDQ3t7FFyNXa45q67go5ruDzpS+ZafLM4OpwXG0AAFyio01QaToGp/k9u1zdWu+/3xiEKitNoLntNs9uLtfMr7Q001VGyOkwDDJmkDEAwB/XaspSY7dW0yDTtAXnvfe835C06dT1pueSWKenlZgmDgBAe3J1bTUdYOxq7XFNXXd1T7m6uaTGW1ZMnmwCTF2dmZo+bpzZpk41+6NGmcfjx5m23k5owaEFBwDQ3nyNwTl1yvNWFN40785iJWY3WnAAAOhM0dFmPE6vXp4ztZKTPRcoTE29fsHCAwdMoJFavvkoLT0tIuAAANBRIiPNFHPXXdj37WvcdwUdV3eWdP1KzKdOeQaea9cIOz4QcAAA6EjNx/O49vftMwORs7I8W3zS0szPzW8+6rohadPWHbgRcAAA6Aq8LTzY9F5cWVnmdVfgGT9eOnrU/Ny0WwuSCDgAAHRtvm4+mpPj2brT9JYTYKE/AABCStNbT7jus8UMq+twNQAACFVNww480EUFAABsh4ADAABaJ4TW4CHgAAAA/7wtOtiFEXAAAIB/zRcd7OLT0gk4AADAv+aLDnbxaenMogIAAP651uAJkWnpXbt2AACg6wihael0UQEAANsh4AAAANsh4AAAANsJasBZt26dMjIyFBsbq4SEhOteP3TokB544AElJSUpJiZGY8aM0U9/+lO/5718+bIefPBBxcXFKSEhQYsXL1ZFRUUQ3gEAAAhFQR1kXFNTo3nz5ik9PV2vvvrqda/n5uaqX79+2rJli5KSkpSVlaVly5YpIiJCjz/+uM/zPvjgg7p48aJ27Nih2tpaPfLII1q2bJneeOONYL4dAAAQIsIsy7KC/Us2bdqkFStWqLS01G/Zxx57TCdOnNDOnTu9vn7ixAmNHTtWH330kdL+dz7+9u3bdeedd+rvf/+7Bg0a5Pd3lJeXKz4+XmVlZYqLiwvovQAAgM4RyPd3lxuDU1ZWpj59+vh8PTs7WwkJCe5wI0lz5sxReHi49u/f7/WY6upqlZeXe2wAAMC+ulTAycrK0tatW7Vs2TKfZYqLi9WvXz+P5yIjI9WnTx8VFxd7PWb9+vWKj493b0lJSe1abwAA0LUEHHAyMzMVFhbW4nby5MmAK3L06FHdfffdWrNmjebOnRvw8S1ZtWqVysrK3Nu5c+fa9fwAAKBrCXiQ8cqVK7Vw4cIWyyQHeH+K48ePa/bs2Vq2bJmefvrpFssOGDBAly5d8niurq5Oly9f1oABA7we43A45HA4AqoTAAAIXQEHnMTERCUmJrZbBY4dO6ZZs2ZpwYIFWrdund/y6enpKi0tVW5urlJTUyVJO3fuVENDg6ZOndpu9QIAAKErqGNwioqK5HQ6VVRUpPr6ejmdTjmdTveaNUePHtXMmTM1d+5cPfXUUyouLlZxcbFKSkrc58jJydHo0aN1/vx5SdKYMWN0++23a+nSpcrJydHevXv1+OOP6/7772/VDCoAAGB/QV0HZ/Xq1dq8ebN7f9KkSZKkXbt2acaMGdq2bZtKSkq0ZcsWbdmyxV1u6NChOnv2rCSpsrJSeXl5qq2tdb/++uuv6/HHH9fs2bMVHh6ue++9Vy+++GIw3woAAAghHbIOTlfDOjgAAISekF4HBwAA4EYRcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0QcAAAgO0ENeCsW7dOGRkZio2NVUJCwnWvHzp0SA888ICSkpIUExOjMWPG6Kc//ekNnxcAAHRvQQ04NTU1mjdvnpYvX+719dzcXPXr109btmzRsWPH9P3vf1+rVq3Sz372sxs6LwAA6N7CLMuygv1LNm3apBUrVqi0tNRv2ccee0wnTpzQzp072/W8TZWXlys+Pl5lZWWKi4sL6FgAANA5Avn+juygOrVaWVmZ+vTp067nrK6uVnV1tXu/vLy8Xc8PAAC6li41yDgrK0tbt27VsmXL2vW869evV3x8vHtLSkpq1/MDAICuJeCAk5mZqbCwsBa3kydPBlyRo0eP6u6779aaNWs0d+7cgI9vyapVq1RWVubezp07167nBwAAXUvAXVQrV67UwoULWyyTnJwc0DmPHz+u2bNna9myZXr66acDrZJfDodDDoej3c8LAAC6poADTmJiohITE9utAseOHdOsWbO0YMECrVu3rt3OCwAAuq+gjsEpKiqS0+lUUVGR6uvr5XQ65XQ6VVFRIcl0S82cOVNz587VU089peLiYhUXF6ukpMR9jpycHI0ePVrnz59v9XkBAED3FtRZVKtXr9bmzZvd+5MmTZIk7dq1SzNmzNC2bdtUUlKiLVu2aMuWLe5yQ4cO1dmzZyVJlZWVysvLU21tbavPCwAAurcOWQenq2EdHAAAQk8g399dapo4AABAe+hyC/2FvIoKacsWKSZG+vBDybKkOXOkL35RGjlSOnFC+uAD6cIFqUcP6a67pJ49pcpK6b33pFmzpEuXpIsXpWnTpPJyafNm6cknpR/+UBo7VoqONsempUknT0obN0rTp0uffCJ94QvSsGFSRIR05Yr0+99LsbFSv37SzTdL6enmdx8/bn7flCnSP/6jOfbCBal/f/M+CgulvDypvl5yOKTJk6V9+8w5H3nE/DxunJSQIL3wghQVJU2caN7fuHHmOpSXm/f+jW+YOgMA0EHoomrPLqqKCikuzoQab2JipKqq9vt9oeTJJ00YysmRzp2Thg+Xioulr3xF2rPHBKGePU0o+vKXpbfeMuVSU6W9e6WiImnUKBMSr12TRoww5z1xQvr8cyksTBo8WCopkWbPlk6fNqHtyhUT0EpKTCgsKzPB8R//UTp82Py+M2ekmhqptNSc67PPpAceMOHz1lulSP4fAABdQSDf3wSc9gw4r74qLVnSfudD54uLkx580LSU9etnWtoiIkxQGzbMtKgdPGjC7aJFUna2CUhhYVJ4uJSRITmdUkGBOV9YmGlly8szge30aRN6XX8N58+X3nzTtKwdOmR+V+/e5vw/+YkJggQuAN0UAccPWnAQ0oYMkXr1Mq1gR4+a0NXQIA0YYFqpevQwLV7XrpluyYICE8YGDzatVDfdJB05IlVXm/B05Yo5X22tCWAJCeZ8kyebbtG+fTv5DQOAQcDxI6izqBiDY54vK5M2bTJfqAhtCQnmszNtmrR1qwk/lZXmtWHDpE8/lVJSTGhydQkOGSIlJkp9+pjPqOv2KNXVJkwNHmzC1RtvmHMAQCsQcPxgmngHqasz41x27zbjbXr2DI0xOJ98Yp5Dx4iKagw9iYkm3Ofnm89D374mCOXlSf/wD6Z8fb0J+bfcYj5DPXqY/1iUlkobNkj33ks3HmBTBBw/CDjwq7RUev558wUaG2tavDprDE6PHqZ1BK3Xs6cJrT16mD+HigpzLQcOlC5fNq2fp0+bltbKStNq2tBgupgdDun//T/pnnsISkAXQ8Dxg4CDkHP2rHTffdLf/mZaL4I9BodAZfTvb1r1IiJM66GrlW/QIBN+S0tNC1RMjLnukZEmLN15p+k6ZvwS0K4IOH4QcAA/6uqkP/1Jeuwx08LUv3/7jcHp2dN0FXYHiYnmWiYmmmsQFdU4K66uzuyXl5vrFR1tXnvzTcYlAT4QcPwg4ACdrKJCWr1aeuUV04XkcNzYGBzXwHy7GDLEdFuWlZmuzFtvNd2iN91kWuRiY837TUw0XW7vvmsCKGBzBBw/CDiADTmdpmvo889NV1F1deBjcHr2lK5e7ex30ja9e5uWtKoqExAnTTLv/W9/M11pPXtK48ebAfd3321a53r16uxaAwEh4PhBwAHg09mz0v/5P40BKJAxOA0NJkCFit69zaD36GjTWjRwoAmGly+bsVgJCea59esZT4QugYDjBwEHQNCcPSvNm2e61Xr2NIHB3xicm282LS1d2f33m/FW/fubtbxGjzbrZd1xh5SVJY0ZY1qGuPccgoiA4wcBB0CXU1Fhlib4xS9MQKir8z8GJz7ezJjraubPNy0/DoeUnCxt22bqPGSI9M470m9/a7oJgQARcPwg4ACwjaNHzYrolZWmFaWlMTijRknbt3d2jY1vfcusLeWq2yefmEU+337bLOS5bJlZD4q1iNAEAccPAg6Abqu0VHruOTOtf88es56StzE4MTHSgQOdW9eePaVVq8waT9u3m8erVxu7+b74RWnpUgZLdyMEHD8IOADQCsXFJkDk55tVuYuKvI/B+fDDzl0c8r77zLIBX/uamTLvdJoB0itXSlOn0gpkIwQcPwg4ANDOKirMLS7+/GfT0pKYeP0YnL59pVdf7dh6DR4sLV8uffSRuRXHmDGm6+7yZdNitXgxLUAhhIDjBwEHADrJp59KzzxjusAs6/oxOG++Ke3b17F1Sk83rT8JCWYgdFKSqeesWcwI62IIOH4QcACgC6uokH7+c2n/fjNYul8/zzE4UVHm9WBzOKQnnzQLJEZFma6uoUPN8yNH0vXVCQg4fhBwACDEVVRIL79susT69TOtLq4xODk50t69wf39SUnSa6+Z9Y6++U0pO9uEr9JS6Z/+iW6vICHg+EHAAQCb+/RT6d/+zbQAXbpkur1cY3Cqq6Uf/CC4v//RR6XMTOn9980CiLm50vTpdHndIAKOHwQcAOjmKiqkLVvMz2+8Ic2c2TgGZ+BA0wpTWNi+vzM62kzR79/f3P4jIsKsYv3JJ6bra+xYur38IOD4QcABALSork76+GPpD38wU+JdY3ASE6WvfMUMkG5vt9wivf66ua3H8ePSihUmdMGNgOMHAQcA0GYVFdKvf226u44fbxyDc/my9Mgj7Rt+SkrMeQcNMusNuW6B0U0HORNw/CDgAACC4to16Xe/k/76V7PQoGsMzp49JggF+pU7aJB04YIUFuZ57Nixpoutm4UdAo4fBBwAQIerqDAzr/r0aRyD07+/9PDD0pkzbT9vaqoZRB0ZabrWjh83d7UfNsx243oIOH4QcAAAXUZdnXTihAkl/fubx+PHpccfN7fCOHDg+hac5vLyzKrRU6easUMuKSlmPaG6Omn37pCfyUXA8YOAAwAICXV1UkFB4xicxERp4UITgFzS0swYoIICc8f45g4flqZNM3ecj42VPvvMtOoUFJhQFEItPIF8f4d3UJ0AAECgIiPNGJtevaQ775QmT5YOHZKOHTPB5dgxE24iI01YSUnxPD411czKqqw0+5WV0s6d5vYUo0aZx7o681pdnXTqlHls+nOICp3YBgAATJgZO9b78/v3Xz8Gp67OtNy4WnAGDzbdXpJ5dLXkpKeb/dRU81pubmPrUAi18rjQRUUXFQDA7q5daxyDExnZGGYmT5aysnx3b0lmfM/IkR1bXx/oogIAAI2io6XbbjOPkZGmVSYvz4QbV/dWWpopm5bW2IozebJ5LVBdoIsr9NqcAADAjXGN7Wm67xqo7Ao0bR2EXFfX2ELUiV1ctOAAAIDG0BMZ6flzoAoKrh/j0wkIOAAAoP007e5qaxdXO6CLCgAAtJ/m3V2dNAOLgAMAANpX8zE+nSCoXVTr1q1TRkaGYmNjleDllu+HDh3SAw88oKSkJMXExGjMmDH66U9/2uI5z549q8WLF+uWW25RTEyMhg8frjVr1qimpiZI7wIAAISaoLbg1NTUaN68eUpPT9err7563eu5ubnq16+ftmzZoqSkJGVlZWnZsmWKiIjQ448/7vWcJ0+eVENDg37+859rxIgROnr0qJYuXaqrV6/qRz/6UTDfDgAACBEdstDfpk2btGLFCpWWlvot+9hjj+nEiRPauXNnq8//wx/+UC+//LIKWjlSm4X+AAAIPYF8f3e5MThlZWXq06dPux5TXV2t6upq9355eXmb6wcAALq+LjVNPCsrS1u3btWyZctafUx+fr5eeuklPfrooz7LrF+/XvHx8e4tKSmpPaoLAAC6qIADTmZmpsLCwlrcTp48GXBFjh49qrvvvltr1qzR3LlzW3XM+fPndfvtt2vevHlaunSpz3KrVq1SWVmZezt37lzA9QMAAKEj4C6qlStXauHChS2WSQ5wUZ/jx49r9uzZWrZsmZ5++ulWHXPhwgXNnDlTGRkZ+q//+q8WyzocDjkcjoDqBAAAQlfAAScxMVGJiYntVoFjx45p1qxZWrBggdatW9eqY86fP6+ZM2cqNTVVGzduVHh4l+ppAwAAnSyoyaCoqEhOp1NFRUWqr6+X0+mU0+lURUWFJNMtNXPmTM2dO1dPPfWUiouLVVxcrJKSEvc5cnJyNHr0aJ0/f16SCTczZszQkCFD9KMf/UglJSXu4wAAAKQgz6JavXq1Nm/e7N6fNGmSJGnXrl2aMWOGtm3bppKSEm3ZskVbtmxxlxs6dKjOnj0rSaqsrFReXp5qa2slSTt27FB+fr7y8/M1ePBgj9/XATPeAQBACOiQdXC6GtbBAQAg9ATy/c3gFQAAYDtdbqG/juBqtGLBPwAAQofre7s1nU/dMuBcuXJFkljwDwCAEHTlyhXFx8e3WKZbjsFpaGjQhQsX1Lt3b4WFhXV2dQJSXl6upKQknTt3jvFDTXBdfOPaeMd18Y7r4hvXxruOvC6WZenKlSsaNGiQ3yViumULTnh4+HUzsEJNXFwcf8G84Lr4xrXxjuviHdfFN66Ndx11Xfy13LgwyBgAANgOAQcAANgOASfEOBwOrVmzhntrNcN18Y1r4x3XxTuui29cG++66nXploOMAQCAvdGCAwAAbIeAAwAAbIeAAwAAbIeAAwAAbIeA04XU19frmWee0S233KKYmBgNHz5czz77rN97bnzwwQdKSUmRw+HQiBEjtGnTpo6pcAdpy3X54IMPFBYWdt1WXFzcgTXvGFeuXNGKFSs0dOhQxcTEKCMjQx999FGLx9j9MyMFfl3s+pn58MMPddddd2nQoEEKCwvTb3/7W4/XLcvS6tWrNXDgQMXExGjOnDk6ffq03/Nu2LBBw4YNU3R0tKZOnaqcnJwgvYPgCMZ1+cEPfnDd52f06NFBfBftz991eeeddzR37lzddNNNCgsLk9PpbNV53377bY0ePVrR0dEaP368/vjHP7Z/5Zsh4HQhzz33nF5++WX97Gc/04kTJ/Tcc8/p+eef10svveTzmMLCQn3961/XzJkz5XQ6tWLFCi1ZskTvvfdeB9Y8uNpyXVzy8vJ08eJF99avX78OqHHHWrJkiXbs2KHXXntNR44c0dy5czVnzhydP3/ea/nu8JmRAr8uLnb7zFy9elUTJkzQhg0bvL7+/PPP68UXX9Qrr7yi/fv3q2fPnrrtttt07do1n+fcunWrnnrqKa1Zs0Yff/yxJkyYoNtuu02XLl0K1ttod8G4LpI0btw4j8/Pnj17glH9oPF3Xa5evaovf/nLeu6551p9zqysLD3wwANavHixDh48qHvuuUf33HOPjh492l7V9s5Cl/H1r3/dWrRokcdz3/zmN60HH3zQ5zHf/e53rXHjxnk8N3/+fOu2224LSh07Q1uuy65duyxJ1ueffx7k2nWuyspKKyIiwvrDH/7g8XxKSor1/e9/3+sx3eEz05br0h0+M5Ks3/zmN+79hoYGa8CAAdYPf/hD93OlpaWWw+Gwfv3rX/s8z5QpU6zHHnvMvV9fX28NGjTIWr9+fVDqHWztdV3WrFljTZgwIYg17VjNr0tThYWFliTr4MGDfs9z3333WV//+tc9nps6dar16KOPtkMtfaMFpwvJyMjQ+++/r1OnTkmSDh06pD179uiOO+7weUx2drbmzJnj8dxtt92m7OzsoNa1I7XlurhMnDhRAwcO1Ne+9jXt3bs32FXtcHV1daqvr1d0dLTH8zExMT7/59gdPjNtuS4udv/MNFVYWKji4mKPz0N8fLymTp3q8/NQU1Oj3Nxcj2PCw8M1Z84c23yG2nJdXE6fPq1BgwYpOTlZDz74oIqKioJd3S6vs/7N6ZY32+yqMjMzVV5ertGjRysiIkL19fVat26dHnzwQZ/HFBcXq3///h7P9e/fX+Xl5aqqqlJMTEywqx10bbkuAwcO1CuvvKK0tDRVV1frF7/4hWbMmKH9+/crJSWlA2sfXL1791Z6erqeffZZjRkzRv3799evf/1rZWdna8SIEV6P6Q6fmbZcl+7ymWnKNb7I2+fB19ijTz/9VPX19V6POXnyZHAq2sHacl0kaerUqdq0aZNGjRqlixcvau3atfrKV76io0ePqnfv3kGtc1fm69+cYI9vI+B0IW+99ZZef/11vfHGGxo3bpx7fMSgQYO0YMGCzq5ep2nLdRk1apRGjRrl3s/IyNCZM2f0wgsv6LXXXuuoqneI1157TYsWLdLNN9+siIgIpaSk6IEHHlBubm5nV61TBXpdutNnBsHRtFX51ltv1dSpUzV06FC99dZbWrx4cSfWrHuii6oL+c53vqPMzEzdf//9Gj9+vB566CE9+eSTWr9+vc9jBgwYoE8++cTjuU8++URxcXG2+J+41Lbr4s2UKVOUn58fpFp2nuHDh2v37t2qqKjQuXPnlJOTo9raWiUnJ3st3x0+M1Lg18Ubu35mXAYMGCBJXj8Prtea69u3ryIiIgI6JtS05bp4k5CQoJEjR9r6M9Qavv7NCfbnhYDThVRWVio83POPJCIiQg0NDT6PSU9P1/vvv+/x3I4dO5Senh6UOnaGtlwXb5xOpwYOHNieVetSevbsqYEDB+rzzz/Xe++9p7vvvttrue7wmWmqtdfFG7t/Zm655RYNGDDA4/NQXl6u/fv3+/w8REVFKTU11eOYhoYGvf/++7b5DLXlunhTUVGhM2fO2Poz1Bqd9m9OUIcwIyALFiywbr75ZusPf/iDVVhYaL3zzjtW3759re9+97vuMpmZmdZDDz3k3i8oKLBiY2Ot73znO9aJEyesDRs2WBEREdb27ds74y0ERVuuywsvvGD99re/tU6fPm0dOXLE+va3v22Fh4dbf/nLXzrjLQTV9u3brT/96U9WQUGB9ec//9maMGGCNXXqVKumpsayrO75mbGswK+LXT8zV65csQ4ePGgdPHjQkmT95Cc/sQ4ePGj97W9/syzLsv7jP/7DSkhIsP7nf/7HOnz4sHX33Xdbt9xyi1VVVeU+x6xZs6yXXnrJvf/mm29aDofD2rRpk3X8+HFr2bJlVkJCglVcXNzh76+tgnFdVq5caX3wwQdWYWGhtXfvXmvOnDlW3759rUuXLnX4+2srf9fls88+sw4ePGi9++67liTrzTfftA4ePGhdvHjRfY6HHnrIyszMdO/v3bvXioyMtH70ox9ZJ06csNasWWP16NHDOnLkSFDfCwGnCykvL7e+/e1vW0OGDLGio6Ot5ORk6/vf/75VXV3tLrNgwQJr+vTpHsft2rXLmjhxohUVFWUlJydbGzdu7NiKB1lbrstzzz1nDR8+3IqOjrb69OljzZgxw9q5c2cn1D74tm7daiUnJ1tRUVHWgAEDrMcee8wqLS11v94dPzOWFfh1setnxjX9vfm2YMECy7LMlOhnnnnG6t+/v+VwOKzZs2dbeXl5HucYOnSotWbNGo/nXnrpJWvIkCFWVFSUNWXKFGvfvn0d9I7aRzCuy/z5862BAwdaUVFR1s0332zNnz/fys/P78B3deP8XZeNGzd6fb3pdZg+fbq7vMtbb71ljRw50oqKirLGjRtnvfvuu0F/L2GW5WeZXAAAgBDDGBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7BBwAAGA7/x8pXIvla9jAAAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log(N_c), zsan[gene] , c = \"r\", s = 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "505549fb-01c0-42b3-81c3-7482f5e3f523",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((3005, 19972), (3005, 19972))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cnt_mat_san.shape, cnt_mat.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "df0ec063-34b4-4f6e-8539-fd1633c0b91b",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng_key = random.PRNGKey(0)\n",
    "rng_key, rng_key_ = random.split(rng_key)\n",
    "\n",
    "# Run NUTS.\n",
    "kernel = NUTS(sanity_model)\n",
    "num_samples = 4000\n",
    "num_chains = 8\n",
    "mcmc = MCMC(kernel, num_warmup=1000, num_samples=num_samples, num_chains = num_chains,progress_bar = False)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "36a48542-6fd3-44d1-802a-af8e192590e9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min 19s, sys: 5.43 s, total: 2min 25s\n",
      "Wall time: 33.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "mcmc.run(\n",
    "    rng_key_, n_g, N_c\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "61014336-78c5-4947-84ac-2c11321d762c",
   "metadata": {},
   "outputs": [],
   "source": [
    "estimates = mcmc.get_samples()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "96ffe5ed-0ac9-48f5-869a-d1d81d8dbd92",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha_g = estimates[\"α_g\"].mean()\n",
    "delta_gcs = estimates[\"δ_gc\"].mean(axis = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "e44b46a6-e6bd-4dd9-a4a7-b1412f9803fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "vals = np.log(alpha_g*np.exp(delta_gcs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "7567a955-f80f-4a65-a5e7-c61403c33f0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "R_mine = np.corrcoef(np.log(N_c), vals)[0, 1]\n",
    "R_paper = -0.035"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "75f5c8f8-01d8-4a30-8c74-472a33210d19",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(8, -11.3, 'R_sanity : -0.035')"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGdCAYAAAAfTAk2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdSUlEQVR4nO3de1xUdf4/8NfAyE0EUlFEQUW8AF4BL9jN26rlt9WtNSszycytzV9Z7fbVtrS2zNUuZq3f2rZSKyvbLuZma7l5SzENBC+oeAHDG0peQOIiA5/fH+8Oc2GGmQFmBobX8/GYxzBnzjnzOYfB8/JzOzqllAIRERGRF/HxdAGIiIiIGhsDDhEREXkdBhwiIiLyOgw4RERE5HUYcIiIiMjrMOAQERGR12HAISIiIq/DgENEREReR+/pAnhCdXU1zpw5gzZt2kCn03m6OEREROQApRSuXLmCyMhI+PjUXUfTIgPOmTNnEBUV5eliEBERUT2cPHkSXbp0qXOdFhlw2rRpA0BOUEhIiIdLQ0RERI4oLi5GVFRUzXW8Li0y4GjNUiEhIQw4REREzYwj3UvYyZiIiIi8DgMOEREReR0GHCIiIvI6LbIPDhERkTOqqqpQWVnp6WK0CK1atYKvr2+D98OAQ0REVIeSkhKcOnUKSilPF6VF0Ol06NKlC4KDgxu0HwYcIiIiG6qqqnDq1CkEBQUhPDyck8O6mFIKhYWFOHXqFHr27NmgmhwGHCIiIhsqKyuhlEJ4eDgCAwM9XZwWITw8HCdOnEBlZWWDAg47GRMREdnBmhv3aaxzzYBDREREXocBh4iIiLwOAw4REVELNmLECMyZM8fTxWh0DDhERAQYDMCRI/JMzV5qaip0Oh0eeOCBWu899NBD0Ol0SE1NBQB8/vnneO6559xcQtdjwCEiaukMBiAlBejdW54ZcrxCVFQUPv74Y5SVldUsKy8vx4cffojo6OiaZW3btnXo7tzNDQMOEVFLl5sLpKfLz+np8toa1vI0K4mJiYiKisLnn39es+zzzz9HdHQ0Bg0aVLPMsomqW7dueOGFFzBjxgy0adMG0dHReOutt8z2ffLkSdx+++0ICwtD27ZtMXHiRJw4ccLVh+QUBhwiopYuJgZITpafBw+W15ZYy9MszZgxAytWrKh5/e677+Lee++1u93LL7+M5ORkZGZm4o9//CMefPBB5OTkAJC5gcaNG4c2bdrg+++/x44dOxAcHIzx48fj6tWrLjsWZzHgEBG1dHo9sHMnkJMDpKXJa0uWtTxHjrA2x1keqAG7++67sX37dvz000/46aefsGPHDtx99912t7v55pvxxz/+EbGxsfjf//1ftG/fHps3bwYArFmzBtXV1Xj77bfRr18/xMXFYcWKFcjPz8eWLVtcfESO40zGREQkoaZXL9vva7U86enyfM89QEaG/Lxzp/VQREZaDZh2/tx0zsLDwzFhwgSsXLkSSilMmDAB7du3t7td//79a37W6XSIiIjA+fPnAQB79+7FsWPHavXbKS8vx/Hjxxv3ABqA30giIrJPq+XJzZWLdUKCLNf67NQVjsh6Pyc3nbMZM2Zg9uzZAIDly5c7tE2rVq3MXut0OlRXVwOQm48mJSVh9erVtbYLDw9vYGkbDwMOERE5RqvlMRiMtTm2+uyQOdMaMDefM61vjE6nw7hx4xq8v8TERKxZswYdOnRASEhII5TQNdgHh4iInONInx0y58Fz5uvri0OHDuHgwYMNunmlZurUqWjfvj0mTpyI77//Hnl5ediyZQsefvhhnDp1qhFK3DgYcIiIyHlabQ7DjeM8eM5CQkIarbYlKCgI27ZtQ3R0NG699VbExcXhvvvuQ3l5eZOq0dEppZSnC+FuxcXFCA0NRVFRUZP6ZRARUdNSXl6OvLw8dO/eHQEBAZ4uTotQ1zl35vrNGhwiIiLyOgw4RERE5HUYcIiIiMjrMOAQERGR12HAISIiIq/DgENERERehwGHiIiIvA4DDhEREXkdBhwiIiKql9TUVEyaNMnTxbCKAYeck5oKNNEvMxERGRUWFuLBBx9EdHQ0/P39ERERgXHjxmHHjh2N9hnLli3DypUra16PGDECc+bMabT9NwQDTnOVmgrodPJo1Qro3h144gmgvNy1n7tsGWDyZW5y8vOBCROAoCCgQwfgz3+WOx/X5eJFYOpUICQECAsD7rsPKCkxvp+TA4wcCXTsCAQEyF2An3oKqKy0vr+PP5bfi2UQ/PxzYOxYoF07eT8rqwEHSkRUt9tuuw2ZmZlYtWoVjhw5gnXr1mHEiBG4cOFCo31GaGgowsLCGm1/jYkBpzkbPx44exbIzQWWLgX+8Q9gwQLXfmZoqISApqiqSsLN1atyt95VqySMzZ9f93ZTpwLZ2cDGjcBXXwHbtgGzZhnfb9UKuOce4NtvJey8+irwz39aP9cnTgB/+hNw/fW13/vlF+C664DFixtwkERE9l2+fBnff/89Fi9ejJEjR6Jr164YMmQI5s2bh9/+9rcAgFdeeQX9+vVD69atERUVhT/+8Y8oMfnP3cqVKxEWFoZvvvkGcXFxCA4Oxvjx43H27NmadUybqFJTU7F161YsW7YMOp0OOp0OeXl5iI2NxUsvvWRWvqysLOh0Ohw7dsxl54ABpznz9wciIoCoKKktGDNGLtKOOHFCahE++UQuxoGBwODBwJEjwI8/AsnJQHAwcNNNQGGhcTvLJqoRI4CHH5bao7ZtpTzPPGP+WZcvAzNnAuHhUksyahSwd28DDtyGb78FDh4EPvgAGDhQyv7cc8Dy5RJ6rDl0CNiwAXj7bWDoUAkgr78utTBnzsg6MTHAvfcCAwYAXbsCv/2thKLvvzffV1WVLH/2WdnG0rRpErbGjGnUwyYishQcHIzg4GCsXbsWFRUVVtfx8fHBa6+9huzsbKxatQqbNm3CE088YbZOaWkpXnrpJbz//vvYtm0b8vPz8ac//cnq/pYtW4aUlBTcf//9OHv2LM6ePYvo6GjMmDEDK1asMFt3xYoVuOGGGxAbG9s4B2zt+Fy2Z3KvAwek1sLPz7ntFiyQ5pY9ewC9HrjrLgkry5bJBfzYMfs1IKtWAa1bA7t2AUuWAH/9q3nQmjwZOH8e+M9/gIwMIDERGD1amoZsGTFCwpQzdu4E+vWTpiTNuHFAcbHU0NjaJixMAp1mzBjAx0eOx5pjxyQU3Xij+fK//lWaxe67z7lyE1GLYDDI/yHttZo3Br1ej5UrV2LVqlUICwvDtddeiyeffBL79u2rWWfOnDkYOXIkunXrhlGjRuH555/HJ598YrafyspKvPnmm0hOTkZiYiJmz56N7777zupnhoaGws/PD0FBQYiIiEBERAR8fX2RmpqKnJwc7N69u2afH374IWbMmOG6EwAGnObtq6+kliUgQC7s589LnxNn/OlPEgLi4oBHHpEA8vTTwLXXAoMGycV68+a699G/vwSlnj2lKSc5GdD+ALZvB3bvBv71L1nesyfw0ksSKj791PY+o6OBTp2cO5aCAvNwAxhfFxTY3qZDB/Nler3URlluM3y4nOuePaXW669/Nb63fTvwzjvSdEVEZMFgAFJSgN695dkdIee2227DmTNnsG7dOowfPx5btmxBYmJiTafg//73vxg9ejQ6d+6MNm3aYNq0abhw4QJKS0tr9hEUFIQePXrUvO7UqRPOnz/vVDkiIyMxYcIEvPvuuwCAf//736ioqMDkyZMbfpB1YMBpzkaOlI6qu3YB06dLM8pttzm3j/79jT9rYaBfP/Nl9r7MpvsAJJho2+zdKx1227WTMKY98vKA48dt7/O994BFi2y/f9NNxn0lJNRdvsayZo3UdH34IbB+vQQ1ALhyRZqf/vlPoH1795SFiJqV3FwgPV1+Tk+X1+4QEBCA3/zmN3j66aeRlpaG1NRULFiwACdOnMD//M//oH///vjss8+QkZGB5cuXAwCumjTpt2rVymx/Op0OSimnyzFz5kx8/PHHKCsrw4oVKzBlyhQEBQU17ODs0Lty5wsXLsT69euRlZUFPz8/XL58udY6Dz/8MHbs2IEDBw4gLi4OWXZGlly8eBELFizAt99+i/z8fISHh2PSpEl47rnnEBoa6poDaapatwa09st335U+Iu+841wTiemXV6ezvqy62vF9WG5TUiKBZ8uW2ts1pLPy228DZWXmnx8RIbVFps6dM75nTURE7QBnMEjzmeU2UVHyHB8v/W1mzQIef1yC2okTwC23GNfVjl+vl47JJv8DIqKWJyZGKrHT06W7o7Vueu4QHx+PtWvXIiMjA9XV1Xj55Zfh4yN1HZbNU/Xh5+eHqqqqWstvvvlmtG7dGm+88QY2bNiAbdu2Nfiz7HFpwLl69SomT56MlJQUvPPOOzbXmzFjBnbt2mXWNmjLmTNncObMGbz00kuIj4/HTz/9hAceeABnzpzBp3U1eXg7Hx/gySeBxx6TfjSBgZ4ukUhMlKYevR7o1q3x9tu5c+1lKSnAwoUSWLRmp40bpWNzfLz1/aSkSCfojAwgKUmWbdokAWXoUNufX10tw8Srq4E+fYD9+83ff+opqdlZtswYjIioxdLrpctfbq6EG71Lr77AhQsXMHnyZMyYMQP9+/dHmzZtkJ6ejiVLlmDixImIjY1FZWUlXn/9ddxyyy3YsWMH3nzzzQZ/brdu3bBr1y6cOHECwcHBaNu2LXx8fGr64sybNw89e/ZESkpKIxxl3Vx6ip999lkAMJsEyNJrr70GQCYkciTg9O3bF5999lnN6x49emDhwoW4++67YTAYoHf1t6YpmzxZ+uAsXy59a5qCMWMkREyaJB2Qe/WS0Unr1wO/+515515T99wjIaauZipLY8dKkJk2TT6roECCxkMPyYgzQGp47rlH+gh17ix9j8aPB+6/H3jzTQkts2cDd9wBREbKNqtXSy1Rv36yn/R0YN48YMoUWd6qFdC3r3lZtNop0+UXL8o8PdrorJwceY6IsF3DREReQ6+XfwLdITg4GEOHDsXSpUtx/PhxVFZWIioqCvfffz+efPJJBAYG4pVXXsHixYsxb9483HDDDVi0aBHuueeeBn3un/70J0yfPh3x8fEoKytDXl4euv36n9v77rsPL7zwAu69995GOEL7vCINFBUVISQkxGa4qaioMBsmV1xc7K6iuZdeLxfnJUuABx+UJixP0+mAr78G/vIX6SNUWCgX8xtuqN0h2FR+vtRKOcPXVzpeP/ighKrWraVvkmln4NJSCRamk/StXi3nbfRo+czbbgN+Dd4A5LwuXizDH5SSoeKzZwOPPupc+datk3OgueMOeV6woPbQeiKiBvD398eiRYuwqI7/JD766KN41OLfsWnTptX8nJqailSL0ayTJk0y64NjWYHRq1cv7Ny50+rnnT59Gq1atWpwiHKUTtWnt5CTVq5ciTlz5ljtg6N55plnsHbtWrt9cCz9/PPPSEpKwt13342FCxfa3LdWm2RKC0ZERETWlJeXIy8vD927d0dAQICni9MsVVRUoLCwENOnT0dERARWr15d5/p1nfPi4mKEhoY6dP12ehTV3Llza2YotPU4fPiws7utl+LiYkyYMAHx8fF4po7/Ac+bNw9FRUU1j5MnT7qlfERERC3dRx99hK5du+Ly5ctYsmSJ2z7X6Saqxx9/vFaVlaUYN3QPv3LlCsaPH482bdrgiy++qDWUzZS/vz/8tT4YLcULL8jDmuuvl0n3iIiIXMxaU5c7OB1wwsPDER4e7oqyOKy4uBjjxo2Dv78/1q1bx2pDax54ALj9duvvNZURVkRERC7i0k7G+fn5uHjxIvLz81FVVVXTvyY2NhbBwcEAgGPHjqGkpAQFBQUoKyurWSc+Ph5+fn44ffo0Ro8ejffeew9DhgxBcXExxo4di9LSUnzwwQcoLi6u6TQcHh4OX19fVx5S89G2rTyIiIhaIJcGnPnz52PVqlU1rwcNGgQA2Lx5M0aMGAFAZjfcunVrrXW0oWWVlZXIycmpmTp6z5492PXrPYIsb9JlOhyNiIiosbhhPA79qrHOtVtGUTU1zvTCJiKilquyshLHjh1DZGRky5st30OKiopw5swZxMbG1upf68z12yvmwSEiInIFvV6PoKAgFBYWolWrVjW3NSDXqK6uRmFhIYKCgho8cS8DDhERkQ06nQ6dOnVCXl4efvrpJ08Xp0Xw8fFBdHQ0dNr9EeuJAYeIiKgOfn5+6Nmzp9ldtsl1/Pz8GqWmjAGHiIjIDh8fH05J0sywMZGIiIi8DgMOEREReR0GHHKd1FRg0iRPl4KIiFogBpzmKjUV0Onk0aoV0L078MQTQHm5p0tmtGwZsHKl8fWIEcCcOR4qzK+WLwe6dQMCAoChQ4Hdu+1v869/AX36yDb9+gFff23+/jPPyPutWwPXXAOMGQP8OhlljW7djL8v7fG3vzXSQRERkSUGnOZs/Hjg7FkgNxdYuhT4xz+ABQs8XSqj0FAgLMzTpTBaswZ47DE5R3v2AAMGAOPGAefP294mLQ24807gvvuAzEypkZo0CThwwLhOr17A3/8O7N8PbN8uYWbsWKCw0Hxff/2r/L60x//7fy44SCIiAhhwmjd/fyAiAoiKkovumDHAxo2ObXv1KjB7NtCpk9RMdO0KLFpkfP+VV6S2onVr2f8f/wiUlBjfX7lSwss33wBxcUBwsDFwaUybqFJTga1bpVZHq8HIywNiY4GXXjIvW1aWvH/smJMnxI5XXgHuvx+4914gPh54800gKAh4913b2yxbJsf15z/LcT73HJCYKIFGc9ddcu5jYoCEBPmc4mJg3z7zfbVpI78v7dG6deMeHxER1WDA8RYHDkhtg5+fY+u/9hqwbh3wySdATg6werXUPGh8fGSd7Gxg1Spg0yZpAjNVWirh5P33gW3bgPx84E9/sv55y5YBKSkSMLQajOhoYMYMYMUK83VXrABuuEHCjzVbtkgAOnHCsWMFJNBlZEgQMT3GMWOAnTttb7dzp/k2gNT62Nrm6lXgrbek9mrAAPP3/vY3oF07YNAg4MUXAYPB8fITEZFTOA9Oc/bVV1JzYjAAFRVywTatWahLfj7Qsydw3XUSFrp2NX/ftK9Mt27A888DDzwA/N//GZdXVkotSI8e8nr2bGmGsSY0VMJXUJDUXmhSU4H586UvzJAhss8PP6xdq2MqKAjo3Vv6Hjnq55+BqiqgY0fz5R07AocP296uoMD6NgUF5su++gq44w4JfZ06SU1a+/bG9x9+WGp+2raVIDpvnoS8V15x/BiIiMhhrMFpzkaOlOacXbuA6dOl6eW22xzbNjVVtu3dWy6+335r/v5//wuMHg107ixNK9OmARcuyAVcExRkDDeAXNjr6s9iTWQkMGGCsZno3/+WsDZ5su1thgyRUNK5s/X3v/9egp/2WL3auTLVh/a7SEuTJq3bbzc/F489Jp2s+/eXoPjyy8Drr8uxEhFRo2PAac5at5ZmnAEDJCDs2gW8845j2yYmSh+Y554Dysrkgvz738t7J04A//M/cjH+7DNp2lm+XN4znarcsgZFpwPqc3P6mTOBjz+WcqxYAUyZIuGpvpKTJWxoj9/+VmpTfH2Bc+fM1z13zrxGyVJEhGPbaL+LYcPkd6DX1/27GDpUat6caWYjIiKHMeB4Cx8f4MkngaeekqDgiJAQCRP//KeMMPrsM+DiRQk01dVSyzBsmIwSOnOm4WX085NmIks33ywB4Y03gA0bpF9OQwQGStjQHm3ayGcnJQHffWdcr7paXqek2N5XSor5NoA0P9W1jbbvumpnsrLkd9ahg93DISIi5zHgeJPJk6WWQqttqcsrrwAffSRNPUeOyFwvEREyMio2VvrCvP66DEF//33pa9NQ3bpJLdOJE9Inprpalvv6SpPZvHnSL8heeNi9W+adOX3auc9/7DEJc6tWAYcOAQ8+CPzyizTtae65R8qheeQRCV0vvyzn6plngPR06W8EyPZPPgn88APw008SDmfMkLJpzWw7dwKvvgrs3Svnc/Vq4NFHgbvvlnlziIio0THgeBO9Xi68S5bIhbcubdrIesnJwODBEjq+/lpqFQYMkAC0eDHQt69ckE2HkNfXn/4kYSY+HggPl47Omvvuk+Yv07BhS2mpjPyqrHTu86dMkc7L8+cDAwdKLcqGDeadiPPzzYe6Dx8unZ7fekvOy6efAmvXynkB5HgOH5a+T716AbfcIn2Vvv9ehowDMpz/44+BG2+UZQsXSsB56y3nyk9ERA7TKVWfThPNW3FxMUJDQ1FUVISQkBBPF4cACQSjRwMnT9YetURERATnrt8cJk6eVVEhM/4+84w06TDcEBFRI2ATlbd64QXzodKmj5tu8nTpjD76SObguXxZmsyIiIgaAZuovLWJ6uJFeVgTGGh7DhkiIqImik1UJDPmtm3r6VIQERF5BJuoiIiIyOsw4BAREZHXYcAhIiIir8OAQ0RERF6HAYeIiIi8DgMOEREReR0GHCIiIvI6DDjUuAwGuTu5weDpkhARUQvGgEONx2AAUlKA3r3lmSGHiIg8hAGHGk9uLpCeLj+np8trIiIiD2DAocYTEwMkJ8vPgwfL6+aGTWxERF6BAYcaj14P7NwJ5OQAaWnyujlhExsRkddgwKHGpdcDvXo1v3ADsImNiMiLMOAQaTzRxMYmMSIil3BpwFm4cCGGDx+OoKAghIWFWV3n4YcfRlJSEvz9/TFw4ECn9q+Uwk033QSdToe1a9c2uLzUwrm7iY1NYkRELuPSgHP16lVMnjwZDz74YJ3rzZgxA1OmTHF6/6+++ip0Ol19i0dUmzub2NgkRkTkMi79V/zZZ58FAKxcudLmOq+99hoAoLCwEPv27XN431lZWXj55ZeRnp6OTp06NaicRB6hNYmlpzffUWdERE1UM+wJCpSWluKuu+7C8uXLERERYXf9iooKVFRU1LwuLi52ZfGIHKM1ieXmSrhpjh2ziYiaqGbZyfjRRx/F8OHDMXHiRIfWX7RoEUJDQ2seUVFRLi4hkYOa86gzIqImzOmAM3fuXOh0ujofhw8fdkVZAQDr1q3Dpk2b8Oqrrzq8zbx581BUVFTzOHnypMvKR0RERJ7n9H8bH3/8caSmpta5TowL+xJs2rQJx48frzUq67bbbsP111+PLVu21NrG398f/v7+LisTERERNS1OB5zw8HCEh4e7oiwOmTt3LmbOnGm2rF+/fli6dCluueUWD5WKiIiImhKXNvzn5+fj4sWLyM/PR1VVFbKysgAAsbGxCA4OBgAcO3YMJSUlKCgoQFlZWc068fHx8PPzw+nTpzF69Gi89957GDJkCCIiIqx2LI6Ojkb37t1deThERETUTLg04MyfPx+rVq2qeT1o0CAAwObNmzFixAgAwMyZM7F169Za6+Tl5aFbt26orKxETk4OSktLXVlUIiIi8iI6pZTydCHcrbi4GKGhoSgqKkJISIini0NEREQOcOb63SyHiRMRERHVhQGHiMiSdhPU8nLeDJWomWLAISIyZXoT1HbteDNUomaKAYeIyJTpTVC1wQ28GSpRs8OAQ0RkSrsJKgAEBckzb4ZK1OzwBjhERKZMb4IaHQ3k5/NmqETNEP9iiYgsaTdBBYzPrqZ1bNY+k4GKqEHYREVE5GkGAzBsGJCQII9hw9ipmaiBGHCIiDwtNxfIyDC+zshgp2aiBmLAISLytJgYICnJ+Do5mZ2aiRqIjbxERJ6m1wM//OD6PjgGg9QMsdM0tQCswSEiagr0eiA+Xh6uCjfaBIacuJBaAAYcIqKWwHQCQ05cSC0AAw4RUUtgOoEhJy6kFoCNsERELYHpBIbsg0MtAL/hREQthekEhkRejk1URERE5HUYcIiIiMjrMOAQEVHNrbA4epy8BQMOEVEL5+gUOaYhqNECEZMVuQgDDhFRC+fIFDmmIWjYMHk0eM5ATj5ILsSAQ0TUwjkyRY5pCMrIMN4btEFzBnLyQXIhBhwiohZOmyInJwdIS7M+RY5pCEpONt4btEFzBnLyQXIhzoNDRER2p8ixnCcQaIQ5Azn5ILkQv01EROQQyxDUKHMGcvJBchE2UREREZHXYcAhIiIir8OAQ0RERF6HAYeIiFouTjTotRhwiIioZeJEg16NAYeIiFomTjTo1RhwiIioZeJEg16N8+AQEVHLxIkGvRp/m0RE1HJxokGvxSYqIiIi8joMOERERO7AIelu5bKAs3DhQgwfPhxBQUEICwuzus7DDz+MpKQk+Pv7Y+DAgQ7ve+fOnRg1ahRat26NkJAQ3HDDDSgrK2ucghMRETU2Dkl3O5cFnKtXr2Ly5Ml48MEH61xvxowZmDJlisP73blzJ8aPH4+xY8di9+7d+PHHHzF79mz4+LAyioiImigOSXc7l3UyfvbZZwEAK1eutLnOa6+9BgAoLCzEvn37HNrvo48+iocffhhz586tWda7d+/6F5SIiMjVtCHp6ekcku4mzara4/z589i1axc6dOiA4cOHo2PHjrjxxhuxffv2OrerqKhAcXGx2YOIiMhttCHpOTlAWhqHpLtBswo4ub9W6T3zzDO4//77sWHDBiQmJmL06NE4evSoze0WLVqE0NDQmkdUVJS7ikxERCS0IekMN27hVMCZO3cudDpdnY/Dhw+7qqyorq4GAPzhD3/Avffei0GDBmHp0qXo3bs33n33XZvbzZs3D0VFRTWPkydPuqyMRERE5HlOxcjHH38cqampda4T48J2xU6dOgEA4uPjzZbHxcUhPz/f5nb+/v7w9/d3WbmIiIioaXEq4ISHhyM8PNxVZbGrW7duiIyMRE5OjtnyI0eO4KabbvJQqYiIiKipcVlDYH5+Pi5evIj8/HxUVVUhKysLABAbG4vg4GAAwLFjx1BSUoKCggKUlZXVrBMfHw8/Pz+cPn0ao0ePxnvvvYchQ4ZAp9Phz3/+MxYsWIABAwZg4MCBWLVqFQ4fPoxPP/3UVYdCREREzYzLAs78+fOxatWqmteDBg0CAGzevBkjRowAAMycORNbt26ttU5eXh66deuGyspK5OTkoLS0tGadOXPmoLy8HI8++iguXryIAQMGYOPGjejRo4erDoWIiIiaGZ1SSnm6EO5WXFyM0NBQFBUVISQkxNPFISIiIgc4c/1uVsPEiYiIiBzBgENERERehwGHiIiIvA4DDhEREXkdBhwiIiLyOgw4RERE5HUYcKhRGQzAkSPyTERE5CkMONRoDAYgJQXo3Vuem2PIYUAjIvIODDjUaHJzgfR0+Tk9XV43J54IaAxURESuwYBDjSYmBkhOlp8HD5bXzYm7A5o31HgRETVVDDjUaPR6YOdOICcHSEuT182JuwNac6/xIiJqyhhwqFHp9UCvXs0v3ADuD2jNvcaLiKgpa4aXISLX0QKauz5r506puYmJaZ6hkIioqWINDpEHNecaL2+mdf4uL2cncKLmigGHiMiEaefvdu3YCZyouWLAISIyYdr5u7RUnt01qu7gQXkwTBE1HAMOEZEJ087fQUHy7OpO4AYDMGwYkJAgj2HDGHKIGooBh4jIhOlougsX3DOqLjcXyMgwvs7IcE2NESeWpJaEAYeIyILW+TsgwD2dwGNigKQk4+vk5MavMeLEktTScOwGEZGH6fXADz9I7QrgmlBlbWJJd02JQOQJrMEhImoC9HogPl4erqgx4sSS1NIw4BARtQBN9VYq7BdErsKAQ0TUQjS1iSXZL4hciQGHiIg8gjecJVdiwCEiIo9gvyByJQYcIiLyiKbQL4h9gLwXAw4REXmMJ/sFsQ+Qd2PAISKiFol9gLwbAw4REbVI7u4DxOYw92LAISKiFsmdfYDYHOZ+DDhERNRiuasPEJvD3I8Bh4iIyMU4JN79msh8lkRERN5Law7LzZVw01Rmk/ZmPMVERERuoDWHkXuwiYqIiIi8DgMOEREReR2XBZyFCxdi+PDhCAoKQlhYmNV1Hn74YSQlJcHf3x8DBw50aL8FBQWYNm0aIiIi0Lp1ayQmJuKzzz5rvIITERFRs+eygHP16lVMnjwZDz74YJ3rzZgxA1OmTHF4v/fccw9ycnKwbt067N+/H7feeituv/12ZGZmNrTIRERE5CVc1sn42WefBQCsXLnS5jqvvfYaAKCwsBD79u1zaL9paWl44403MGTIEADAU089haVLlyIjIwODBg1qWKGJiIjIKzS7PjjDhw/HmjVrcPHiRVRXV+Pjjz9GeXk5RowY4emiERERURPR7IaJf/LJJ5gyZQratWsHvV6PoKAgfPHFF4iNjbW5TUVFBSoqKmpeFxcXu6OoRERE5CFO1eDMnTsXOp2uzsfhw4ddVVYAwNNPP43Lly/jv//9L9LT0/HYY4/h9ttvx/79+21us2jRIoSGhtY8oqKiXFpGIiIi8iydUko5unJhYSEuXLhQ5zoxMTHw8/Oreb1y5UrMmTMHly9ftrnNM888g7Vr1yIrK6vOfR8/fhyxsbE4cOAAEhISapaPGTMGsbGxePPNN61uZ60GJyoqCkVFRQgJCanzM4mIiKhpKC4uRmhoqEPXb6eaqMLDwxEeHt6gwjVEaWkpAMDHx7ziydfXF9XV1Ta38/f3h7+/v0vLRkRERE2HyzoZ5+fnIysrC/n5+aiqqkJWVhaysrJQUlJSs86xY8eQlZWFgoIClJWV1axz9epVAMDp06fRp08f7N69GwDQp08fxMbG4g9/+AN2796N48eP4+WXX8bGjRsxadIkVx0KERERNTMu62Q8f/58rFq1qua1NoR78+bNNSOeZs6cia1bt9ZaJy8vD926dUNlZSVycnJqam5atWqFr7/+GnPnzsUtt9yCkpISxMbGYtWqVbj55ptddShERETUzDjVB8dbONOGR0RERE2DM9fvZjcPDhEREZE9DDhERETkdRhwPMlgAI4ckWciIiJqNAw4nmIwACkpQO/e8syQQ0RE1GgYcDzlyBEgPV1+Tk8HcnON7zXFmp2mWCYiIiIbGHBcpa5AYDAA99xjfJ2cDMTEGN+zVbNT35DR0HDC2iYiImpmGHBcwV4gyM0FMjKMr1etAvR643vWanbqG3waI5xYlum772x/licCWGPvh4iImj0GHFewFVI0MTFSawMAgwcDvXrZfk+r2alP8HGkLI4wLVNQEDB+vHxWebkxUNQ3SFnbzpGgYrmOtf2UlwPffCPP9ra3t5yIiJoX1QIVFRUpAKqoqMg1H1BZqVRyslKAUoMHy2tr6+TkOP6erX3m5Mgy7ZGT43xZHD2mDRvMP6tvX3lOTlYqO7vucthiWf7sbGN5k5Ntnx/LdSz3s2+fUkFB8nNQkFJXrhjPqbXtre23rMz274iIiNzOmes3a3BcQa8Hdu4EcnKAtDRj85PlOr16Of6erX3aqvFxpiyOHtPo0cbP6tcPOHBAftZqiOoqhy2W5Tfdn60aJ2u1Upb7OXkS+PUWHygtBRITjbU7tjp4W+538GD7NVL1qfGxVvvkzD5Yy0REZJ8bAleT4/IaHHerqzbIVZ9VVla7Zqi+5TDdztHaL2vrmO6nrMxYgxMQYLuWyHJ7bXm/fvZrpGzVBNk7VstaIu11UpKUTTsP2dlK7d1rXFbXZzpz7t35fSEiakTOXL8ZcKj+XHWhdGS/jqxTVibNaleuOB7G6gpwluw1D1pjuY1ls58WdBITay+z1hSnHYMjQauyUpruEhLM17UMmAw/RNREOXP95s02ebPNlsFgMDZlOdpMZ28brWOz1pzlSBOg5TbbtgHXX29sGqtLTo6UxfIzc3OlKc10PdOO69rnDhtmPnoPALKzgenTZX8JCRKbDh6Un99/H/D3l8/Mz3fs3JWXA1u3AjfeCAQE2D8ma+rzuyKiFoE322wC7HWTMBjkOnLwoOMDeeo78Kcxumxo+7AcONWQriS2yueSLiZ19Xmq7zb16d9kuU1AgLzOzgaSkmSd5GTpM2RKmyvJ2mfa64cF1J6aAJB+VIAxXGVnyxdS+zkxUYJOu3YSoAYNAkpKbP9yystl3fHj5XnfPtu/1Lq+zKaj4Uy/cEREznB5fVIT5OomKnstBpWV0uJg2fpQ17bOLne0LM4ej9atJSnJeAyWXUmc+Zy6uqTYa22xNtDMdJnWQlVWVvs9ey1UjrbW1KdFx2ZZyiw+3FofnPoWxPJL17ev8cRoo+EceQQG2v7lrF9fe31rv1TT14mJsl1ZmezDsgmuZ0/r/ZPYjEbUIrEPjh2uDjj2umZYvm+6jq1tnV3uaFnqczzWHpZdSeo7StyR/VgLbZbLrlwxhrHAQMfCmOk+LAOco6PV7WloH2NtH/bCndXPLqtU2V8eVXu/OKqy91Uaj/tKmcrucYvaiwS1Dwlqrz6x5uds9FGV/q2t/9Itfzl79zr25bDW7ygwUPoHmZ4Qf//a68XHy8PaSbc8CVeuKPWPfyi1Z4/tXyCDElGzwoBjhztrcKz1T7X8z7Sti6ytAT6OLHe0LM4ejxYakpONxzB4sGN9ch05V/Xt22u57O23nQ9jdQU5a0HLE32MHQl3tq7lpt85bZ9lZbWXWz4SBxrU3h4TVTb6qLKAMJWDnqoyeZjVL9uVAdeqt5GqLrUKN65XVqYqk4ZafX0FgWoDfqOuIFDeTxqqKq+UqZy3ttS8Vwa/un8xlZUSrkw7T1+6ZL5eXJyxlsjyhFhWoWZnO1ZrRkRux4BjhztGUTnSYmDr31FHmk+c/ayG/kfVdHCRrSac+n6Os/uxFtosl9mqwakrRFkGBdNt7NXgOBrq6gp09h7WgpytZZZshTdrAauuR1BQtZyfpOpax3vlilI6XfWv6xrXKyuTZ2uvdTCYPSdht0pKKDNbFoASVYhQtR7j1B70U3u1mqW+AyUsJQ5R2ehjrHGCr1LPP6+uIFD9H2aqT/FbtRcJqjK+f80vqXJvttqLBLUON6t9SFCVe7Nrp8D4eAlOlhM+OvLHQEQuwYBjB4eJN3/e1gcnO9u8ScyyBsdabY2tcOdoDY7WRGavBsfRWi1btWaOtFDZf1TVWpaUWKXK9uaoJOwyXx6wX13Ku6Twa0CqqYnCLlWZnaMqK5VKjPvFfJu4X1RZRrbKRh+Vgf5qHW5WGeiv9iJB7dEnq/UYp64MuFblfHlIlQ1KUfuQoNb5/FbCUeIQVZk4RO1Fglof/QdVePSSenvOPnXpbJnx98wARNQoOEzcDg4Tp6bIdHQ0YBw8ZDmYy9ooakdGVpsOXDLdp+lynU4u+YDx53vvBfbskWVBQTIxtLVR8SUlQEiIcXug9kh4y9faZ+h0CkrpkJykoKD7dcCXAqCr85xtWG/A+Am1D/j554Gnnqq9fk62jMbqnVB7mwT/Y8iuiK3j06oA+CIIv6AUrWuWJiIdCtXIxJBfl2jllufkpGp8r65H/p5CRCeGI//9rYjpJSc+d+tJxNwYBX2AlEcbXXnqFHDDDUBenkzKbfrzqFH1H4FP1Nw5c/1mwGHAIaqTFoAA+1PilJQAa9YAt90GnD9vXM8ygGmvIyOBHTuAa68Fzpwxv7dshw7AkCHA0aPWy5WcDHz/PXDdtQoZe3RmyzduBK65xnz9pESFH3bpgCNHMDShCHswuOa9BOxHNvo15DTVqS/24QD6IwglKEUwkhKrgcOHkVEaj+SgbOy80BvQ6zF0iMKeTDkWa2ETAAICFJY8XoBO/cLRO16Pnj3Nfyfl5cCmTUBUFGrei46Wc1pRIeGpa1cgPp7TDFHzw4BjBwMOUfNgGoS2bQM6dQJ8fa3XQFVVyXva8pISmauwQwe50Ndc0A0GGIZdh4MZJfhJ3wvdDEfQMzEU1x3+JzJK4+soza81OLoylKrAmqWJA6qgKiqQeTjo1yXmNTj9sA/70b/O48zZkAdERVmtWbInKEihtFSH5GRg82Y53rIy7T2pcQsMNC7TxMUBq1fLOauqkuATESHPOp2Eolat5PyeOgV06ybbMBSRJzHg2MGAQ9TCackpOrqm+sNgAI5sOony9lE4fU6Pzp0Bvc6AquMncPbnVrjh9x1xJuMsoq+NwtHDVTjx7/3odks/xA30BwwGHEy8C6cOXcEQ30x8WXUTbmu1HucrwxCd1AHXYxvSM3xqanCSsRsqIAgZ5X0xOCgbaRd6A7m5GJpwGXt+berSoQoKvrV+rsvbbwMzZ7rutCUlAT/8ID8fPAicOCE1Rf7+ciq3bZOAdfYs0L693I938mSpnbPW3ErkLAYcOxhwiKjRWYYm0/AEPXKPGBA97Ubk7ylETHI7YPNm5O44a+yDYzDAENcPB4/54hSicAO2Ig8xOIkuuAHbkIcYHEEPTMbnZmFHC02DkxU2baxCh2uuogxSm6TV7lirwamv7Gxg2jRjvyxnJSUB770nTWq5ucDly9KkNns2EBYm6/BuHWQLA44dDDhE5BH2rtwlJeZtTFaUIBAfYCra4zx64zh6+p1E/tWOiElqC/27b6F8wGBswkhE4RR67ngP+VeuQfTQTshdvRMVsx/BT+iKJ7EQh+w0m1mjhZOEBKc3rUVrPjN16RIQHGy83VpystyZBDDvs6Xd6sxgAA4dkpqkLl3klJo2U5L3YcCxgwGHiJosrZfw009LNUmvXnJVHzFCqk8CAmQda9avByZMML7W1rWowjHAF0fQE1XwhQJQBR+cQhdE4CxOoQt00CG6uy9affIxDMoXp3afRrfISsTd1E06Qw9V2LOn7hFu9fH889KkZXrvWNP7wWqdrYOCgHPn5JRY3mINMDalWRtlCNgeoUhNHwOOHQw4RNTk1TX0bNgw8xu0ZmQYx+AnJ8t7jSEmBvDzAw4fltdJScD27TBcNwIHM0pwIuoGRH2yFP4h/tIHZ7MBHa4cx9mDl9A+NgQHfgrF5Ic64sx5fc00BPfeK8V1pAZn8GBg5UrrNUb2+hvl5Eh40e7fmp4uxVfKvHnNMgwBxoypdWoHpOysHfI8Bhw7GHCIqFmzrJIwDULl5ZIMDhywPowqPl6u8ocOyeu6aoSs2bBB7hiv6dlT7hwPWA9XiYkynO3XZGDaVcmRPjiAMaA4WoOjNW3p9VJbY1ojZI0WhgA5Fe3a1Q5fmoQEKYs2F1F5ObB1q7HZjFyLAccOBhwi8mrWOjzn5sp72pXcdHKj3FzZxnRWR2v69gV+/NEYoDSxsVLTc/Cg7W0TEoCPPjKONdfG95eXAwUFdc5gaDlvkjN9cExrcJKTgepq80M0DUMA8M035vnN1mnIzJR9a2EoKAi4cKF28DGd48m0TNrpZ42Qcxhw7GDAISKywjR03HOPeW2MadVFSYmMA6+ocP4zkpKAd98FZswwr34JDAQuXjQml0YcRuVMHxx7NTianByZXdo0DG3YAIwbZ74PrdbJtMP0sGHGQ/+11c9y4B1Djw0MOHYw4BAR2WHr3h6A/XafP/4R+L//c/4z168Hxo61PYzKTVd+a31wqqpkeHx2tvFWJbZqcGzVAuXkyLPlqevb17xF0bJWyVZ3LGstlN6OAccOBhwiogYwbfexHKmlTancsaPxyt+1q7HPT1327AFat7Y9jCoxEXjuOWNzlpsnzLH2cdb64FirwdFCEWBeg6OFG0vWOklrtyfR7uVm2sfcMhR5K2eu3z5uKhMREXkLvV6upjk5MvQpJ0eal3JyZHlwsFRnbNggz/v2SVDZt894VU5OBr74wny//v6SHpKT5fXgX+8Xlp4uz3v2yDD4du2kmSwlRcJQfLx0iikvN9Y6ueiwLZu0AgKkWcq0+1BAgPHwi4vltGg3p9XrZdRWdrY8fvzReLhBQcbDNq2d0Q4/PV3ClPY6I8MYlNLTjd2sSLAGhzU4RETuY9m+olVnmFZB2BpGZcraOHGtNqmu6owmOE2ylTuHWO0krc0EoNXgJCdL7ZA2S4AWorwZm6jsYMAhImoiHAkcBoOM0EpJMTZ7nTsnTVa2bjdvOvbbNEFo6cByVFcT1pT64Hg6H7KJioiImgdHphPW64H+/c2bvYKDpcnK37/2+v36Ga/+WhVI795SzaHVBGVnyz6HDavdpFVeLj2FnZkfyIUsT5Hpa3fOxmx6KlNSXNYS2GgYcIiIqHmw7PASHCwzBa5bZ5zuuF8/YPdu4xXftBPLgQMyMaGpjAzzzitaD+Hx4+V5zx6pPWrqV3M3sOwP1NT7/Lg04CxcuBDDhw9HUFAQwrQpKk3s3bsXd955J6KiohAYGIi4uDgsW7bM7n4vXryIqVOnIiQkBGFhYbjvvvtQUlLigiMgIqImLSAAuOUWICtLmqX27DHv8WvZaXnPHvN7PyQnG2t7AOnFq02CU1oqnaITEoCBA1t82LE8laanrSlyaaXW1atXMXnyZKSkpOCdd96p9X5GRgY6dOiADz74AFFRUUhLS8OsWbPg6+uL2bNn29zv1KlTcfbsWWzcuBGVlZW49957MWvWLHz44YeuPBwiImqqtLYaa8t37jTvOJKVZXsq4RtvtH6jLNN7f5newMrTnVLcyNqpbMrc0sl45cqVmDNnDi5fvmx33YceegiHDh3Cpk2brL5/6NAhxMfH48cff0Tyr1Fyw4YNuPnmm3Hq1ClERkba/Qx2MiYiIpv27QMGDKh7nZwcucpro8Cs3bUTaFEByB2adSfjoqIitG3b1ub7O3fuRFhYWE24AYAxY8bAx8cHu3btsrpNRUUFiouLzR5ERERWxccb22ISEuS1Ka1Z68gR40Q0GRnGWiFA+vJ8/TUwdGjz6ZXrZZpUnExLS8OaNWuwfv16m+sUFBSgQ4cOZsv0ej3atm2LgoICq9ssWrQIzz77bKOWlYiIvJRlWwxQ9w2sLFm7oZXWK1ebnpi1Oi7ndA3O3LlzodPp6nwcPnzY6YIcOHAAEydOxIIFCzB27Fint6/LvHnzUFRUVPM4efJko+6fiIi8jOVY7Ph4GVYeH28MJb16mc/MrPUBMu2orNF65ZqOtR40SMKQ6X2/qNE4HR0ff/xxpKam1rlOjJNdqw8ePIjRo0dj1qxZeOqpp+pcNyIiAufPnzdbZjAYcPHiRURERFjdxt/fH/7W5kogIiKqL+2+C5a1MaYdlYOCpDZIC0ZHjpgPW09OllFfWj+e996rXUPEGp96cfpMhYeHIzw8vNEKkJ2djVGjRmH69OlYuHCh3fVTUlJw+fJlZGRkIOnX5Lxp0yZUV1dj6NChjVYuIiIiu6yN3tJuRmV5F05AQorpHTazs43vZWRInx/L21Z4+O7qzZVLOxnn5+cjKysL+fn5qKqqQlZWFrKysmrmrDlw4ABGjhyJsWPH4rHHHkNBQQEKCgpQWFhYs4/du3ejT58+OH36NAAgLi4O48ePx/3334/du3djx44dmD17Nu644w6HRlARERG5nLW7cAISSH78UUIOIKFFa+bSmM6iZzm73pEj5tMJu/gGo82acqHp06crALUemzdvVkoptWDBAqvvd+3atWYfmzdvVgBUXl5ezbILFy6oO++8UwUHB6uQkBB17733qitXrjhcrqKiIgVAFRUVNdKREhEROaGyUqmcHHmurFQqO1uppCSlAKUGD5Zl2nrJycbl2dnys/bo21eek5ON23gxZ67fvNkm58EhIqKmwFZfG1t3V+/XD9i/37ie6Q1GvVSzngeHiIioRbI1BN1yRNfOnRJmdu9uXvdOcDP2TiIiImpOTDs2N6d7J7gZzwYREVFzZeseXMQmKiIiIvI+DDhERETkdRhwiIiIyDHN6LYSDDhERERkn+l9tJrB3dEZcIiIiMg+y1mVtdmWmygGHCIiIrIvJqZZzbvDYeJERERknzbJYDOZd6dpl46IiIiajmY07w6bqIiIiMjrMOAQERGR12HAISIiIq/DgENERERehwGHiIiIvA4DDhEREXkdBhwiIiLyOgw4RERE5HU40V8jKykBPvgACAwEtm0DlALGjAH69pW5kQ4dArZsAc6cAVq1Am65BWjdGigtBb75Bhg1Cjh/Hjh7Fhg2DCguBlatAh59FHjxRSA+HggIkG2Tk4HDh4EVK4AbbwTOnQN69gS6dQN8fYErV4B//xsICgI6dAA6d5b7o505Axw8KJ83ZAjwu9/JtmfOAB07yrb+/sBvfiOfRURE1NzolFLK04Vwt+LiYoSGhqKoqAghISGNtt+SEiAkREKNNYGBQFlZo32cW9x1lwSkwYOBH36Q0HTvvfJzQgIQFgYsXQr4+QEDB0qAS0iQc3Hpkszm3aoVcOGCBLndu4GTJ4EePYCCAuD664Ht2yUEtm4tge6664BPPpH1kpKAHTuA/Hy5gW2vXkB5ORAbK+U7dEg+R6cDunQBCguB0aOBo0cltF25ImGtsFDKUVQkwfF3vwP27ZPPO34cuHoVuHxZ9nXhgjH0KSXlOnhQQmJpqQTA3bslSMbHA5mZcrwzZsgs5hcuSHl8fIDhw4GsLOM96XQ6CaE5OXI8R4/Kd0L7zkyZAqxZAzz2GDBrFgMmEZEpZ67fDDiNGHDeeQeYObPRdkeE3r0lNP/8swSmK1ck9AQHS0g8cEBCV3U1EBEhIa5VKwmE5eVSa5ebK2GsSxcJce3aAfv3AxUVQJs2ss/gYKCyUgJYWJjsr107CWCVlVKzN2yYhK/qavlMQPb7889AYqLsUwuU0dFAeDjQtq3UMp48Ketrnzl+PPDaa0D79p45r0TUPDHg2MEaHKKmoVUrqf3z95cAFh4uzbPHjknNWfv2EohycoBrrpFtqqok5HXvLrWArVrJ39758xK+qqrkdjkVFfKej4+87+sLdOoEXLwozbtHj8rfZGmpcduQEFk/Kgr4/HNjTSERNQ0MOHa4KuAAzb8PTvv2wNy5xv9xE7Vk11wD/PKL/D316mVsxoyMlHB1+bIEtMBA+ZvU62V5YKBsExwsQSs8XP6m/PwksJWUAAaDvC4ulhovvV7+Zg0GqUXr2lX+trt08fBJIGpCGHDscGXA8QYGg/QxWbFC/nH292cfHHf1wbntNuCVV9z+K6cmbOBA+TvYvVuaDXNz5T8tZ89KcDp3Tv7+fH2BuDj5Ph07BgwaJN/Bn36S73br1kC/fsDp07J9YCCwa5f8vfzjH2wupOaBAccOBhxqysrLgfffB156SS5Qbdp4Vx8cg8Ez55Xq1ru3NM1VVEgzXliYfJeOHZPv3A03yPeqY0ephe7TB8jLA266CUhLk/90aDVYkZHyfYyKAu67T75fRI2BAccOBhwizykpAZYskU75hYVNrw9OZSX7yjW2UaPk9xgdLbWzw4ZJ7a5Sxtqlc+ek9vZf/5L1rrlGfi+JifLdyM0FHnlEghe1XAw4djDgEFFdfvgBmDBBAtKlS+7tg3P1qjRzknV33SXB99IlCcjl5fL7iIqSubvWr5dmYX9/CcITJsj5/OEHCUfjxgGTJnEKhuaKAccOBhwiaspKSoAXXpC+MT16AP37N34fHH9/4L//9fSReoZeDzzwgITR4cOBL74ARo6UABQdLWHp3DkgNFTOUZ8+xppGToDqWQw4djDgEBFJ7cennxqbCzt1apw+OB07At9+K/vxNn5+EnoOHJDAePAgcOutMsDgl1+kWS00FEhPBx5/HPjuO2DiRCAjQ0a7Mhw1DAOOHQw4RESuZTDIxT8zE1i9WmqN6tMHp18/YMECTx9N4wgIABYvNt4Sx9dXBgKcOSPnS3t94oScuzlz2OfIEgOOHQw4RETNR3m59K354QfpT9OmjfN9cNLSZEqK5larVFgoZY6MlLnVOnWSIATIqEetf5i+hdxZkgHHDgYcIqKWx2CQ6Qy2bJHmtlOnHO+D4+sLPPWUNNm5U2Sk1PDodLZnyU9KkvCn1xtrzk6cMM7VpddLENy61dhMZjBIf66YmOYVjhhw7GDAISIiZxkMwJEj0gk8I8O5Pjjbt8t6rrri5uRIWBk6FNizx7g8MVGCTceOMiVCUJCEt5EjpYzJyVJ2LRxpoQdomgGIAccOBhwiInK3khKZxLNtW8f64MyeLZ2409PrrsHRQkpurkzYaOntt81vBG35WgtHKSnyWUlJsjwjwzwANQUMOHYw4BARUXOg1ao40gfHYKhdg5OUJE1ytmpwBg+W/km2whEgAahXr/qVu7FrgJy5fvs03sfWtnDhQgwfPhxBQUEIs9IVfO/evbjzzjsRFRWFwMBAxMXFYdmyZXXu88SJE7jvvvvQvXt3BAYGokePHliwYAGuXr3qoqMgIiLyDL1ewkVwMHDzzTK/Uf/+8ujXz9jHRlt31y5g717gyy/l+YcfZNsLF4ANG+Q5OFhqZXJyJNzo9RJEkpNlP8nJxlqcwYONTVaOMhikNqh3b3n21O1ZXFrpdPXqVUyePBkpKSl45513ar2fkZGBDh064IMPPkBUVBTS0tIwa9Ys+Pr6Yvbs2Vb3efjwYVRXV+Mf//gHYmNjceDAAdx///345Zdf8NJLL7nycIiIiJo0vd4YgEwFBMgszqbrmdbK6PXGZq6G9sHJzZXaIUCec3OdrwFqDG5polq5ciXmzJmDy5cv2133oYcewqFDh7Bp0yaH9//iiy/ijTfeQK52y2Y72ERFRETkGloNjmkTWGM1Uzlz/W4i3YaMioqK0LZt20bdpqKiAhUVFTWvi4uL610+IiIiss2yNshTHZRd2gfHWWlpaVizZg1mzZrl8DbHjh3D66+/jj/84Q8211m0aBFCQ0NrHlFRUY1RXCIiIrJCawLz5OgrpwPO3LlzodPp6nwcPnzY6YIcOHAAEydOxIIFCzB27FiHtjl9+jTGjx+PyZMn4/7777e53rx581BUVFTzOHnypNPlIyIioubD6Wz1+OOPIzU1tc51Ypzscn3w4EGMHj0as2bNwlNPPeXQNmfOnMHIkSMxfPhwvPXWW3Wu6+/vD39/f6fKRERERM2X0wEnPDwc4eHhjVaA7OxsjBo1CtOnT8fChQsd2ub06dMYOXIkkpKSsGLFCvj4NKmWNiIiIvIwlyaD/Px8ZGVlIT8/H1VVVcjKykJWVhZKSkoASLPUyJEjMXbsWDz22GMoKChAQUEBCgsLa/axe/du9OnTB6dPnwYg4WbEiBGIjo7GSy+9hMLCwprtiIiIiAAXj6KaP38+Vq1aVfN60KBBAIDNmzdjxIgR+PTTT1FYWIgPPvgAH3zwQc16Xbt2xYkTJwAApaWlyMnJQWVlJQBg48aNOHbsGI4dO4YuXbqYfV4LnJSZiIiIrOCtGjgPDhERUbPQZG7VQEREROQJDDhERETkdRhwiIiIyOsw4BAREZHXYcAhIiIir8OAQ0RERF6HAYeIiIi8DgMOEREReR0GHCIiIvI6DDhERETkdRhwiIiIyOsw4BAREZHXYcAhIiIir8OAQ0RERF6HAYeIiIi8DgMOEREReR0GHCIiIvI6DDhERETkdRhwiIiIyOsw4BAREZHXYcAhIiIir8OAQ0RERF6HAYeIiIi8DgMOEREReR0GHCIiIvI6DDhERETkdRhwiIiIyOsw4BAREZHXYcAhIiIir6P3dAG8Tnk5sGkTEBEB6PXy6NVLng0GIDcXiIwEtm0DunQB4uPlPY22TnQ0kJ9v/Tk3V9bt1UuejxyRzy0oAEaNAgICau8vJsb4OQaDbFNVBfj6GstnjbXtiYiImjhesRpTeTnQrh1QWmq+PD4e2LULGDkSSE83f2/gQGDVKgkbP/0E/OUvwMGDgL8/UFEhYaW8HAgMBMrKjEEJAPr3B0pKjIEHAFq1Ar77DrhwQfb59NPAoUNAXBzwwgsSrlJTZZlm0CBg5UpApwOUMj5XVABTpgB5ecCAAcD778tyQMoREyNl3bkT6NwZuO464MsvgdtuA86frzuk5efL9gYDsHUrcOON5sFMC4rt2wPnzgHduskx1BWyTMNhbq681sqpfZ5p0LQMfXUFOQY9IqJmRaeUdsVqOYqLixEaGoqioiKEhIQ03o6//hqYMMH6ezEx5kHEG2hByJZWrYDKSmNI0541QUHA1asSHvz8gOPHgW++AUaPllqlykrz/XXuDHzyiYS9n36SbQsLgSFDJHTcfbcELi0MarTP7dNHwtfnn0vAS0oC/vpXYPhwCWfZ2UBsLLBmjXz+jh2yzpo1wNKlUr7kZOA//wFeew0IDQV69ABGjJAaucpKOR/+/sD110uovfZa4MwZ81AH2A5LWu0aYB7MAOPyumrciIi8mDPXbwacxgw4+/ZJTYctkZFysaOWR6s5at0aCA6WWqmuXYHf/x4oLpaQ1LMncM89wOHD5ttERUlo0wJOnz4S9Lp3BzZulH3dfbfsF5Awt3WrMVzZqrlqzFop1nARkRsw4NjhsoBjMEhtQmam9fd79wZychrv84hM3XUX0KED8Oqr5stDQ4HHHpNmyLw8CUj33AO8/LI0ZXbqBMyZI0GqoEACyrXXSlj517+AhARZr3Nnea+8HNizR2qz/PykKXTaNKkB69sX+PHHuvuBaQHMslmSiMgOBhw7XBZwAPnHfP9+ICVF+rBY2rVL/td+8qRxWVycXHi0ZpyuXSUIJSZK08iaNcA//2neZGOveahDB+kHY8nHB6iurv/xEdnTujUwaxYweLB83+bNk+93hw7Aiy8CDzxg7E+2aJF8jwMCJFTFx5v3l7LGVm0Ra5GIvB4Djh0uDTia8nJpPrh6FZg/X/qGDB4MpKXJ+5ajmEw72+r1tf+hLi8Hvv1WQo6fn3RYzs8Hpk4FsrIkJK1eLfvT6yVcJSaalykwUJoszpwBLl2STsPr1wOnThnX6dFD+pr4+UnZAbn4dOggnwdIuHr5Zen4GxYmI7eIGlNQkNQ6lZcDp09Ln6sZM4B+/YDp06UZr0sXqWEKCZE+TsOGGftRrV4tTXbR0dKXauhQ+dvTRi4Cxr8x058ZjIiaNKeu38qFnn/+eZWSkqICAwNVaGhorfezsrLUHXfcobp06aICAgJUnz591Kuvvurw/svLy9WAAQMUAJWZmenwdkVFRQqAKioqcnibBqmsVConR57dte/KSqWSkpQClBo0SKn165UqK6u9fVmZUn37ynqDB8vrnBx53rfPuF1lpVJ79ij11ltKXbli/XNiY5XavVupjAylPv9cHrt2KbVunVJpaUrFxMh6UVFK/fCDUtnZsq9165Tq00fe8/eX5x49lPr735W6dEmpLVuUuuYaWd65s1Lz5im1fbtScXGyrD6PwEB57t1bqd/9zv76vr71/yw+mtYjNFS+R4B8ZyMi5Oc+fZT64gv5Dn/xhVJffinfvw0bjH87ln9vlZXyPc7ONl/mqr93ohbOmeu3S2twFixYgLCwMJw6dQrvvPMOLl++bPb+u+++i7179+LWW29FVFQU0tLSMGvWLCxZsgSzZ8+2u/9HHnkER48exX/+8x9kZmZi4MCBDpXLLTU4TYGjVfYNrdpvjM+xnP/HkeYHbcSRwWBssquqAs6eBW64QfZTUQGcOCHvR0cb5yayHKE0dKj0K9FGfsXFAQsWABcvSv8Svd44bP3MGenvUlYmI6bCw6UWTKv16tYNuP12oG1b4K23vG/0XEv1wQfAww/LdyImRqZbWLpUakMBqR1askRG5h0+LNMvvP221EB16yadyPPzpTb0X/+SEXr9+7PWiMgJTa6JauXKlZgzZ06tgGPNQw89hEOHDmHTpk11rvef//wHjz32GD777DMkJCQw4FDD2AtYlkw7ympNita2NQ1hWrDS5ujR3s/Pl7CUnS0XviNHJGClpcnjmmtk/qTSUuDnn4EPP5SQ1bWrXFyLi11/fsg12reXDt+dOsn3JzAQ+M1vpFnY0fma2PeIWpBmHXDuvvtulJeX49NPP7W5zrlz55CUlIS1a9eiffv26N69OwMOeQdHLlam6wAyp8/JkzKfT1qazKJ99apM0pifLyOb3noLOHAAuOUWYO1aCU3PPCM1Dd27y37y8mTiybfeAnbvltFXc+dKDUT37nLxPXjQDSeB8Pvfy+8gP1/OfVqa/N6WLZPfWVwcsHixDGYYNkz6zfXoIf3xgoMZeshrNduAk5aWhhtvvBHr16/H2LFjra6jlMLNN9+Ma6+9Fk899RROnDhhN+BUVFSgwmREU3FxMaKiohhwyLs5MzszYH9kEmDsHH/pktQ2TZggF9/YWODyZRlK3r27NAlqQ8cTE+XnEyckiEVEyH63bzcPTBERMkyd6k+vBx56SGr5CgtlaorRo+X8Jydbv2VMRITUGHECSWoGXNrJ+H//938VgDofhw4dMttmxYoVVjsZm9q/f79q3769eu655+pcb9myZeraa69VBoNBKaVUXl6eAuruZLxgwQKr5XRbJ2OilsheZ1vtfdOO7VqH9fh46XxeViaPdetkGaCUn588JyYq9cYb0okeUEqnM3/mw/zh76/U8uVKdehg/f0uXZT68EPpYL1vn5z37GwZKGBtcAE7UpMHuLSTcWFhIS5cuFDnOjExMfDz86t5ba8G5+DBgxg5ciRmzpyJhQsX1rnvSZMm4d///jd0Ol3NsqqqKvj6+mLq1KlYtWpVrW1Yg0PUTNSnI7q1G9Rq/Zy0zuclJdJJfOZM6UzeoYPUamizRoeHAx07Sg1TUZFMt1BV5f7jb0os59rS6WTeovBwabosLJRaoIwMmXPr+HH5fYwcWfcM2kQN0KyaqLKzszFq1ChMnz4dS5Yssbuv/Px8FJt0qjxz5gzGjRuHTz/9FEOHDkWXLl3s7oN9cIjI7L5fWvOMtRu26nTyrDW5dekiF/7Tp+Xnv/1Nmufi4qTZ7s03PXpYTUZ8PLBwoYwqy8yUPkIrVsg5HDpU5jRi4CEnNZmAk5+fj4sXL2LdunV48cUX8f333wMAYmNjERwcjAMHDmDUqFEYN24cXnzxxZrtfH19ER4eDgDYvXs37rnnHnz33Xfo3Llzrc9wpA+OJQYcInKaI32WABlhp00poN1zLCZG+rs89ZRc7AcOBCZOlAv9Bx9IrUdpqUwt8MQTMtWBJiRE7jW2b5/0W/IWXbvKZI7XXy9D6LWb2375pXSQDwhgzQ/V0mQCTmpqqtUmo82bN2PEiBF45pln8Oyzz9Z6v2vXrjhx4gQAYMuWLRg5ciTy8vLQrVu3Wusy4BBRs+HoKLmDB2X25q5dpSZEq13S7nUXHS0dvEeNkhmer1xx73G4Q7du0jF90CDg6aclOO7bB7RrJ/dN8/eX82CrgzSbxrxSkwk4TRUDDhE1S3XNg6M1q1VUmE92aTAAv/wCrFsnN1wtKJDapcuXjfMxNXedOsk9z+68U+anysqS29BkZkqt0IoVxtviMOw0aww4djDgEFGLZBqQDAZpSuvUSS7+M2ZIh+HERKkxmTtXbvrrTWJipJ+UdlNXRyf2pCaDAccOBhwiIguWtUOWs3CXl0ttiEYLB4GBMhN3z54yOeHRozJRoXZz3qZIu62KdpwdO0qgmznTfKJEBqAmhwHHDgYcIiInGQwyc3J6OjB4sHSatlULYjpCLTJS1j92zBiGEhKA55+XprKvvgLeeUduQ9IUrFgBvPii9IPSyhsXJ6Plxo6VY9SOzfSedgxAbsGAYwcDDhFRPdS3464jN9M9eFBuYpudLYGivFxuHwLITUlvvRX47DNg//7GPSZn+PlJ81ZWlrz295c+T0lJMsJNu7Eu+/q4DAOOHQw4RERNkLVbgwDm8xQdPAicOiUjyj79VIbW+/nJbUJSU+X9Vq2Aykr3ll0Lb4CMAFu6FBgzxnzSQ9Ob9AYEuLd8XoIBxw4GHCIiL2RtokZAgk9mpnSoNp1jyF0GDQLmzwemTpX5jvz9gXffBdq0Md493hSHuNvEgGMHAw4RUQtiGhjKy2VyRa3mp21b6U+Umir9gcLDZZ6de+81vxmsq7RqJeXx8ZHRbJ07Gz87KQn44QeGHBMMOHYw4BARUZ20jtJVVRI8IiOlWezoUQkef/mLPA4dqn3frsa0b5/c1sK0XC14hJcz1++Wc1aIiIgcpddLh2JTBw+aNx3dcou8joyUWqD27eWmo3fd1XjlOH5cmrQiI2XeonnzzEd49e0L/PijeTMXm7gAsAaHNThERNS4SkqANWvkfmMZGXLPsV27gGuvleBx/LiMvnrySeNIMVvi4yXQ1FVL1Lev9DHSOmIPGyafm5Agw/oDArwm9LCJyg4GHCIi8jhtVJhph+j586XZCzCGG0fk5Mhos/37ZVi9Rgs5118vz8nJwM6d8l4zDDwMOHYw4BARUZNkOkliTIwxmJjW4AQESGdp7XnwYCAtTd4bOFDmEjK1YQMwfrzxdXY2MH26eeDRan+aeOhx5vrt46YyERERkT1a35/4eAkwO3dK7UxxsdxAdO9e4NIlWaY9p6XJdrm5tcNNcrLMu6PdZiM5WZ7T043PWg1SSgrQu7c8e8GNWBlwiIiImiq9XpqegoOBm2+W5qeAAFmmPWu1LTExxgDTt6+MwNJqZzRKma83eLC8zs2tHXqauaZZB0VERETO0esl0Fg2Mx05Ip2OAXnOz6+9nhZ6tHuNRUfLdk24ucoe1uAQERF5C63GxzSUWKuxsVxPC0c5OTLk/frrG9ZcZXo3eg9hwCEiIvJmpuFF669ja71evaSGpyHNVU2kPw8DDhERkbezVrNji7UaH2c0kf48DDhERERk5GiNjy0NDUiNpHn2HCIiIiLX0Wp86ruttc7ObsaAQ0RERI2rIQGpkbCJioiIiLwOAw4RERF5HQYcIiIi8joMOEREROR1GHCIiIjI6zDgEBERkddhwCEiIiKvw4BDREREXocBh4iIiLwOAw4RERF5HQYcIiIi8joMOEREROR1WuTNNpVSAIDi4mIPl4SIiIgcpV23tet4XVpkwLly5QoAICoqysMlISIiImdduXIFoaGhda6jU47EIC9TXV2NM2fOoE2bNtDpdJ4ujlOKi4sRFRWFkydPIiQkxNPFaTJ4XmzjubGO58U6nhfbeG6sc+d5UUrhypUriIyMhI9P3b1sWmQNjo+PD7p06eLpYjRISEgI/8Cs4HmxjefGOp4X63hebOO5sc5d58VezY2GnYyJiIjI6zDgEBERkddhwGlm/P39sWDBAvj7+3u6KE0Kz4ttPDfW8bxYx/NiG8+NdU31vLTITsZERETk3ViDQ0RERF6HAYeIiIi8DgMOEREReR0GHCIiIvI6DDhNSFVVFZ5++ml0794dgYGB6NGjB5577jm799zYsmULEhMT4e/vj9jYWKxcudI9BXaT+pyXLVu2QKfT1XoUFBS4seTuceXKFcyZMwddu3ZFYGAghg8fjh9//LHObbz9OwM4f1689Tuzbds23HLLLYiMjIROp8PatWvN3ldKYf78+ejUqRMCAwMxZswYHD161O5+ly9fjm7duiEgIABDhw7F7t27XXQEruGK8/LMM8/U+v706dPHhUfR+Oydl88//xxjx45Fu3btoNPpkJWV5dB+//Wvf6FPnz4ICAhAv3798PXXXzd+4S0w4DQhixcvxhtvvIG///3vOHToEBYvXowlS5bg9ddft7lNXl4eJkyYgJEjRyIrKwtz5szBzJkz8c0337ix5K5Vn/OiycnJwdmzZ2seHTp0cEOJ3WvmzJnYuHEj3n//fezfvx9jx47FmDFjcPr0aavrt4TvDOD8edF423fml19+wYABA7B8+XKr7y9ZsgSvvfYa3nzzTezatQutW7fGuHHjUF5ebnOfa9aswWOPPYYFCxZgz549GDBgAMaNG4fz58+76jAanSvOCwAkJCSYfX+2b9/uiuK7jL3z8ssvv+C6667D4sWLHd5nWloa7rzzTtx3333IzMzEpEmTMGnSJBw4cKCxim2doiZjwoQJasaMGWbLbr31VjV16lSb2zzxxBMqISHBbNmUKVPUuHHjXFJGT6jPedm8ebMCoC5duuTi0nlWaWmp8vX1VV999ZXZ8sTERPWXv/zF6jYt4TtTn/PSEr4zANQXX3xR87q6ulpFRESoF198sWbZ5cuXlb+/v/roo49s7mfIkCHqoYceqnldVVWlIiMj1aJFi1xSbldrrPOyYMECNWDAABeW1L0sz4upvLw8BUBlZmba3c/tt9+uJkyYYLZs6NCh6g9/+EMjlNI21uA0IcOHD8d3332HI0eOAAD27t2L7du346abbrK5zc6dOzFmzBizZePGjcPOnTtdWlZ3qs950QwcOBCdOnXCb37zG+zYscPVRXU7g8GAqqoqBAQEmC0PDAy0+T/HlvCdqc950Xj7d8ZUXl4eCgoKzL4PoaGhGDp0qM3vw9WrV5GRkWG2jY+PD8aMGeM136H6nBfN0aNHERkZiZiYGEydOhX5+fmuLm6T56l/c1rkzTabqrlz56K4uBh9+vSBr68vqqqqsHDhQkydOtXmNgUFBejYsaPZso4dO6K4uBhlZWUIDAx0dbFdrj7npVOnTnjzzTeRnJyMiooKvP322xgxYgR27dqFxMREN5betdq0aYOUlBQ899xziIuLQ8eOHfHRRx9h586diI2NtbpNS/jO1Oe8tJTvjCmtf5G174Otvkc///wzqqqqrG5z+PBh1xTUzepzXgBg6NChWLlyJXr37o2zZ8/i2WefxfXXX48DBw6gTZs2Li1zU2br3xxX929jwGlCPvnkE6xevRoffvghEhISavpHREZGYvr06Z4unsfU57z07t0bvXv3rnk9fPhwHD9+HEuXLsX777/vrqK7xfvvv48ZM2agc+fO8PX1RWJiIu68805kZGR4umge5ex5aUnfGXIN01rl/v37Y+jQoejatSs++eQT3HfffR4sWcvEJqom5M9//jPmzp2LO+64A/369cO0adPw6KOPYtGiRTa3iYiIwLlz58yWnTt3DiEhIV7xP3GgfufFmiFDhuDYsWMuKqXn9OjRA1u3bkVJSQlOnjyJ3bt3o7KyEjExMVbXbwnfGcD582KNt35nNBEREQBg9fugvWepffv28PX1dWqb5qY+58WasLAw9OrVy6u/Q46w9W+Oq78vDDhNSGlpKXx8zH8lvr6+qK6utrlNSkoKvvvuO7NlGzduREpKikvK6An1OS/WZGVloVOnTo1ZtCaldevW6NSpEy5duoRvvvkGEydOtLpeS/jOmHL0vFjj7d+Z7t27IyIiwuz7UFxcjF27dtn8Pvj5+SEpKclsm+rqanz33Xde8x2qz3mxpqSkBMePH/fq75AjPPZvjku7MJNTpk+frjp37qy++uorlZeXpz7//HPVvn179cQTT9SsM3fuXDVt2rSa17m5uSooKEj9+c9/VocOHVLLly9Xvr6+asOGDZ44BJeoz3lZunSpWrt2rTp69Kjav3+/euSRR5SPj4/673//64lDcKkNGzao//znPyo3N1d9++23asCAAWro0KHq6tWrSqmW+Z1Ryvnz4q3fmStXrqjMzEyVmZmpAKhXXnlFZWZmqp9++kkppdTf/vY3FRYWpr788ku1b98+NXHiRNW9e3dVVlZWs49Ro0ap119/veb1xx9/rPz9/dXKlSvVwYMH1axZs1RYWJgqKChw+/HVlyvOy+OPP662bNmi8vLy1I4dO9SYMWNU+/bt1fnz591+fPVl77xcuHBBZWZmqvXr1ysA6uOPP1aZmZnq7NmzNfuYNm2amjt3bs3rHTt2KL1er1566SV16NAhtWDBAtWqVSu1f/9+lx4LA04TUlxcrB555BEVHR2tAgICVExMjPrLX/6iKioqataZPn26uvHGG82227x5sxo4cKDy8/NTMTExasWKFe4tuIvV57wsXrxY9ejRQwUEBKi2bduqESNGqE2bNnmg9K63Zs0aFRMTo/z8/FRERIR66KGH1OXLl2veb4nfGaWcPy/e+p3Rhr9bPqZPn66UkiHRTz/9tOrYsaPy9/dXo0ePVjk5OWb76Nq1q1qwYIHZstdff11FR0crPz8/NWTIEPXDDz+46YgahyvOy5QpU1SnTp2Un5+f6ty5s5oyZYo6duyYG4+q4eydlxUrVlh93/Q83HjjjTXraz755BPVq1cv5efnpxISEtT69etdfiw6pexMk0tERETUzLAPDhEREXkdBhwiIiLyOgw4RERE5HUYcIiIiMjrMOAQERGR12HAISIiIq/DgENERERehwGHiIiIvA4DDhEREXkdBhwiIiLyOgw4RERE5HUYcIiIiMjr/H/Qor8CIgNzBAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.log(N_c), vals, c= \"r\", s = 2, label = \"Mine\")\n",
    "plt.scatter(np.log(N_c), zsan[gene], c= \"b\", s = 2, label = \"Sanity\")\n",
    "plt.legend()\n",
    "plt.text(8, -11.2, f\"R_mine : {round(R_mine, 4)}\", {\"color\":\"r\"})\n",
    "plt.text(8, -11.3, f\"R_sanity : {R_paper}\", {\"color\": \"r\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "4c22ed12-db3b-4a17-99d6-ec213ce049ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAf7klEQVR4nO3de3CU5dnH8V8OJKBmEyMkITXIQRFQRAolRvH1lJoARR1oFWQYYBBaDbaSVg0eiFQrlDLqyEQZrYqd4VDoCFVALA0CogE1wpQCokAsIG5QaA5gCQm53z867HQhQDbZzV6bfD8zzwx59tnNdSeZzZdnd7NRzjknAAAAQ6LDPQAAAMDpCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYExvuAZqivr5eBw8eVEJCgqKiosI9DgAAaATnnKqrq5Wenq7o6HOfI4nIQDl48KAyMjLCPQYAAGiC/fv369JLLz3nMREZKAkJCZL+u0CPxxPmaQAAQGNUVVUpIyPD93v8XCIyUE49rOPxeAgUAAAiTGOensGTZAEAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCJQI1LVgZbhHAAAgpAgUAABgDoECAADMiQ33AGg8HtoBALQVnEEBAADmECgRgDMnAIC2hkABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQIkTXgpXhHgEAgBZDoBhHmAAA2iICBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwJyAAmXmzJn60Y9+pISEBKWkpOiuu+7Srl27/I45fvy48vLydMkll+iiiy7SyJEjVV5e7nfMvn37NGzYMF1wwQVKSUnRww8/rLq6uuavBgAAtAoBBcr69euVl5enTZs2ac2aNaqtrdXtt9+uY8eO+Y6ZOnWq3nnnHS1dulTr16/XwYMHNWLECN/lJ0+e1LBhw3TixAl99NFHevPNNzV//nxNnz49eKtqA7oWrAz3CAAAhEyUc8419crffvutUlJStH79ev3f//2fKisr1alTJy1cuFA//elPJUmff/65evfurZKSEl133XV699139ZOf/EQHDx5UamqqJGnevHl69NFH9e233youLu68n7eqqkqJiYmqrKyUx+Np6vgR4Vwh8tWsYS04CQAAzRPI7+9mPQelsrJSkpScnCxJKi0tVW1trbKzs33H9OrVS126dFFJSYkkqaSkRH379vXFiSTl5OSoqqpK27dvb/Dz1NTUqKqqym8DAACtV5MDpb6+Xg899JBuuOEGXX311ZIkr9eruLg4JSUl+R2bmpoqr9frO+Z/4+TU5acua8jMmTOVmJjo2zIyMpo6NgAAiABNDpS8vDz985//1OLFi4M5T4OmTZumyspK37Z///6Qf04AABA+sU250pQpU7RixQpt2LBBl156qW9/WlqaTpw4oYqKCr+zKOXl5UpLS/Md8/HHH/vd3qlX+Zw65nTx8fGKj49vyqgRiyfBAgDasoDOoDjnNGXKFC1btkxr165Vt27d/C4fMGCA2rVrp+LiYt++Xbt2ad++fcrKypIkZWVladu2bTp06JDvmDVr1sjj8ahPnz7NWQsAAGglAjqDkpeXp4ULF+qvf/2rEhISfM8ZSUxMVIcOHZSYmKiJEycqPz9fycnJ8ng8evDBB5WVlaXrrrtOknT77berT58+Gjt2rGbPni2v16snnnhCeXl5be4sCQAAaFhAgfLyyy9Lkm6++Wa//W+88YbGjx8vSXr++ecVHR2tkSNHqqamRjk5OXrppZd8x8bExGjFihW6//77lZWVpQsvvFDjxo3Tb3/72+atBAAAtBrN+jso4dIW/g5KY56Dwt9BAQBEkhb7OygAAAChQKAAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkEikFdC1aGewQAAMKKQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUY3iCLAAABAoAADCIQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwIOlA0bNmj48OFKT09XVFSUli9f7nf5+PHjFRUV5bfl5ub6HXPkyBGNGTNGHo9HSUlJmjhxoo4ePdqshQAAgNYj4EA5duyY+vXrp6KiorMek5ubq2+++ca3LVq0yO/yMWPGaPv27VqzZo1WrFihDRs2aPLkyYFPDwAAWqXYQK8wZMgQDRky5JzHxMfHKy0trcHLdu7cqdWrV+uTTz7RwIEDJUlz587V0KFDNWfOHKWnpwc6EgAAaGVC8hyUdevWKSUlRVdeeaXuv/9+HT582HdZSUmJkpKSfHEiSdnZ2YqOjtbmzZtDMQ4AAIgwAZ9BOZ/c3FyNGDFC3bp10549e/TYY49pyJAhKikpUUxMjLxer1JSUvyHiI1VcnKyvF5vg7dZU1Ojmpoa38dVVVXBHhsAABgS9EAZNWqU7999+/bVNddcox49emjdunW67bbbmnSbM2fO1IwZM4I1IgAAMC7kLzPu3r27OnbsqN27d0uS0tLSdOjQIb9j6urqdOTIkbM+b2XatGmqrKz0bfv37w/12AAAIIxCHigHDhzQ4cOH1blzZ0lSVlaWKioqVFpa6jtm7dq1qq+vV2ZmZoO3ER8fL4/H47cBAIDWK+CHeI4ePeo7GyJJZWVl2rp1q5KTk5WcnKwZM2Zo5MiRSktL0549e/TII4/o8ssvV05OjiSpd+/eys3N1aRJkzRv3jzV1tZqypQpGjVqFK/gAQAAkppwBuXTTz9V//791b9/f0lSfn6++vfvr+nTpysmJkb/+Mc/dMcdd6hnz56aOHGiBgwYoA8++EDx8fG+21iwYIF69eql2267TUOHDtXgwYP1yiuvBG9VAAAgogV8BuXmm2+Wc+6sl7/33nvnvY3k5GQtXLgw0E8NAADaCN6LBwAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoFiSNeCleEeAQAAEwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BEoE61qwMtwjAAAQEgQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHNiwz0AeDUOAACn4wwKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMCcgANlw4YNGj58uNLT0xUVFaXly5f7Xe6c0/Tp09W5c2d16NBB2dnZ+vLLL/2OOXLkiMaMGSOPx6OkpCRNnDhRR48ebdZCAABA6xFwoBw7dkz9+vVTUVFRg5fPnj1bL774oubNm6fNmzfrwgsvVE5Ojo4fP+47ZsyYMdq+fbvWrFmjFStWaMOGDZo8eXLTVwEAAFqVgN/NeMiQIRoyZEiDlznn9MILL+iJJ57QnXfeKUn605/+pNTUVC1fvlyjRo3Szp07tXr1an3yyScaOHCgJGnu3LkaOnSo5syZo/T09GYsBwAAtAZBfQ5KWVmZvF6vsrOzffsSExOVmZmpkpISSVJJSYmSkpJ8cSJJ2dnZio6O1ubNmxu83ZqaGlVVVfltAACg9QpqoHi9XklSamqq3/7U1FTfZV6vVykpKX6Xx8bGKjk52XfM6WbOnKnExETflpGREcyxAQCAMRHxKp5p06apsrLSt+3fvz/cIwEAgBAKaqCkpaVJksrLy/32l5eX+y5LS0vToUOH/C6vq6vTkSNHfMecLj4+Xh6Px28DAACtV1ADpVu3bkpLS1NxcbFvX1VVlTZv3qysrCxJUlZWlioqKlRaWuo7Zu3ataqvr1dmZmYwxwEAABEq4FfxHD16VLt37/Z9XFZWpq1btyo5OVldunTRQw89pGeeeUZXXHGFunXrpieffFLp6em66667JEm9e/dWbm6uJk2apHnz5qm2tlZTpkzRqFGjeAUPAACQ1IRA+fTTT3XLLbf4Ps7Pz5ckjRs3TvPnz9cjjzyiY8eOafLkyaqoqNDgwYO1evVqtW/f3nedBQsWaMqUKbrtttsUHR2tkSNH6sUXXwzCciJP14KV4R4BAABzopxzLtxDBKqqqkqJiYmqrKyM+OejNDdQvpo1LEiTAAAQWoH8/o6IV/EAAIC2hUABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQRKGHUtWBnuEQAAMIlAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoYdK1YGW4RwAAwCwCBQAAmEOgAAAAcwgUAABgDoECAADMIVAAAIA5BAoAADAn6IHy1FNPKSoqym/r1auX7/Ljx48rLy9Pl1xyiS666CKNHDlS5eXlwR4DAABEsJCcQbnqqqv0zTff+LaNGzf6Lps6dareeecdLV26VOvXr9fBgwc1YsSIUIwBAAAiVGxIbjQ2VmlpaWfsr6ys1GuvvaaFCxfq1ltvlSS98cYb6t27tzZt2qTrrrsuFOMAAIAIE5IzKF9++aXS09PVvXt3jRkzRvv27ZMklZaWqra2VtnZ2b5je/XqpS5duqikpOSst1dTU6Oqqiq/DQAAtF5BD5TMzEzNnz9fq1ev1ssvv6yysjLdeOONqq6ultfrVVxcnJKSkvyuk5qaKq/Xe9bbnDlzphITE31bRkZGsMcGAACGBD1QhgwZop/97Ge65pprlJOTo1WrVqmiokJLlixp8m1OmzZNlZWVvm3//v1BnLjl8P47AAA0TshfZpyUlKSePXtq9+7dSktL04kTJ1RRUeF3THl5eYPPWTklPj5eHo/HbwMAAK1XyAPl6NGj2rNnjzp37qwBAwaoXbt2Ki4u9l2+a9cu7du3T1lZWaEeBQAARIigv4rnN7/5jYYPH67LLrtMBw8eVGFhoWJiYjR69GglJiZq4sSJys/PV3Jysjwejx588EFlZWXxCh4AAOAT9EA5cOCARo8ercOHD6tTp04aPHiwNm3apE6dOkmSnn/+eUVHR2vkyJGqqalRTk6OXnrppWCPAQAAIljQA2Xx4sXnvLx9+/YqKipSUVFRsD81AABoJXgvHgAAYA6BAgAAzAnJn7rH2fG3UAAAOD/OoAAAAHMIFAAAYA6B0kJC9dBO14KVPGwEAGh1CBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6AAAABzCBQAAGAOgQIAAMwhUAAAgDkECgAAMIdAAQAA5hAoAADAHAIFAACYQ6C0gK4FK8M9AgAAEYVAAQAA5hAoAADAHAIFAACYQ6AAAABzCJQQ4YmxAAA0HYECAADMIVAAAIA5BAoAADCHQAEAAOYQKAAAwBwCBQAAmEOgAAAAc2LDPUBrxt9CAQCgaTiDAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAwh0ABAADmEChBxhsEAgDQfAQKAAAwJzbcA7QWnDkBACB4OIMCtFFdC1YS1gDMIlAAAIA5BAoASTxMCcAWAgUAAJhDoAAAAHMIlGZq6ImGnCqHFfwsAohUBAoQ4U6P5HO9OifQYCFwAIQLgQIAAMwhUIBW5FxnTgI9q3Kusyf8DRUAoUagNAN30AiVcP1sBfPhIQBoDgIFAACYQ6AArcTZXk3WmDMfnB0BYA2BEiDuyGFJqH4ez/f8k1B/fgAgUAAAgDmx4fzkRUVF+sMf/iCv16t+/fpp7ty5GjRoUDhHAlpU14KV+mrWsEYdd7rGXK+xt9XY4xu6bmPXAACBCNsZlD//+c/Kz89XYWGhPvvsM/Xr1085OTk6dOhQuEY6r/99TN/aqW1r8yAwp36mIvX7GKlzo2ma8nJ2IFBhC5TnnntOkyZN0oQJE9SnTx/NmzdPF1xwgV5//fVwjQQAAIwIy0M8J06cUGlpqaZNm+bbFx0drezsbJWUlJxxfE1NjWpqanwfV1ZWSpKqqqpCMt/Vhe9Jkv45I+eMfZaF6uuBs7u68D2/n5PzHduQLlOX+n1cVVWl+prvz9h/vuuF26mfv0C+JmhYY7+GTflaB/oz29Cx9TXfN3h/U1/zvaTm3Rc19Dn/d19D989outO/nv97PxWKr/Gpnw3n3PkPdmHw9ddfO0nuo48+8tv/8MMPu0GDBp1xfGFhoZPExsbGxsbG1gq2/fv3n7cVwvok2caaNm2a8vPzfR/X19fryJEjuuSSSxQVFSXpv1WWkZGh/fv3y+PxhGvUsGjLa5fa9vrb8tqltr3+trx2qW2vP5LX7pxTdXW10tPTz3tsWAKlY8eOiomJUXl5ud/+8vJypaWlnXF8fHy84uPj/fYlJSU1eNsejyfivmHB0pbXLrXt9bfltUtte/1tee1S215/pK49MTGxUceF5UmycXFxGjBggIqLi3376uvrVVxcrKysrHCMBAAADAnbQzz5+fkaN26cBg4cqEGDBumFF17QsWPHNGHChHCNBAAAjAhboNxzzz369ttvNX36dHm9Xl177bVavXq1UlNTm3R78fHxKiwsPOOhoLagLa9datvrb8trl9r2+tvy2qW2vf62svYo5xrzWh8AAICWw3vxAAAAcwgUAABgDoECAADMIVAAAIA5ER0oR44c0ZgxY+TxeJSUlKSJEyfq6NGj5zz+wQcf1JVXXqkOHTqoS5cu+uUvf+l7bx/LioqK1LVrV7Vv316ZmZn6+OOPz3n80qVL1atXL7Vv3159+/bVqlWrWmjS0Ahk/a+++qpuvPFGXXzxxbr44ouVnZ193q+XZYF+709ZvHixoqKidNddd4V2wBALdP0VFRXKy8tT586dFR8fr549e0bsz3+ga3/hhRd8928ZGRmaOnWqjh8/3kLTBs+GDRs0fPhwpaenKyoqSsuXLz/vddatW6cf/vCHio+P1+WXX6758+eHfM5QCXT9b731ln784x+rU6dO8ng8ysrK0nvv2X//uPMKzrvrhEdubq7r16+f27Rpk/vggw/c5Zdf7kaPHn3W47dt2+ZGjBjh3n77bbd7925XXFzsrrjiCjdy5MgWnDpwixcvdnFxce71119327dvd5MmTXJJSUmuvLy8weM//PBDFxMT42bPnu127NjhnnjiCdeuXTu3bdu2Fp48OAJd/7333uuKiorcli1b3M6dO9348eNdYmKiO3DgQAtP3nyBrv2UsrIy94Mf/MDdeOON7s4772yZYUMg0PXX1NS4gQMHuqFDh7qNGze6srIyt27dOrd169YWnrz5Al37ggULXHx8vFuwYIErKytz7733nuvcubObOnVqC0/efKtWrXKPP/64e+utt5wkt2zZsnMev3fvXnfBBRe4/Px8t2PHDjd37lwXExPjVq9e3TIDB1mg6//Vr37lfv/737uPP/7YffHFF27atGmuXbt27rPPPmuZgUMkYgNlx44dTpL75JNPfPveffddFxUV5b7++utG386SJUtcXFycq62tDcWYQTFo0CCXl5fn+/jkyZMuPT3dzZw5s8Hj7777bjds2DC/fZmZme7nP/95SOcMlUDXf7q6ujqXkJDg3nzzzVCNGDJNWXtdXZ27/vrr3R//+Ec3bty4iA6UQNf/8ssvu+7du7sTJ0601IghE+ja8/Ly3K233uq3Lz8/391www0hnTPUGvML+pFHHnFXXXWV37577rnH5eTkhHCyltGY9TekT58+bsaMGcEfqAVF7EM8JSUlSkpK0sCBA337srOzFR0drc2bNzf6diorK+XxeBQba/N9E0+cOKHS0lJlZ2f79kVHRys7O1slJSUNXqekpMTveEnKyck56/GWNWX9p/v+++9VW1ur5OTkUI0ZEk1d+29/+1ulpKRo4sSJLTFmyDRl/W+//baysrKUl5en1NRUXX311Xr22Wd18uTJlho7KJqy9uuvv16lpaW+h4H27t2rVatWaejQoS0yczi1pvu8YKivr1d1dXXE3eedzuZv5Ubwer1KSUnx2xcbG6vk5GR5vd5G3cZ3332np59+WpMnTw7FiEHx3Xff6eTJk2f8hd3U1FR9/vnnDV7H6/U2eHxjvy6WNGX9p3v00UeVnp5+xh2YdU1Z+8aNG/Xaa69p69atLTBhaDVl/Xv37tXatWs1ZswYrVq1Srt379YDDzyg2tpaFRYWtsTYQdGUtd9777367rvvNHjwYDnnVFdXp1/84hd67LHHWmLksDrbfV5VVZX+85//qEOHDmGaLDzmzJmjo0eP6u677w73KM1i7gxKQUGBoqKizrk19hfTuVRVVWnYsGHq06ePnnrqqeYPDpNmzZqlxYsXa9myZWrfvn24xwmp6upqjR07Vq+++qo6duwY7nHCor6+XikpKXrllVc0YMAA3XPPPXr88cc1b968cI8WcuvWrdOzzz6rl156SZ999pneeustrVy5Uk8//XS4R0MLWrhwoWbMmKElS5ac8Z/4SGPuDMqvf/1rjR8//pzHdO/eXWlpaTp06JDf/rq6Oh05ckRpaWnnvH51dbVyc3OVkJCgZcuWqV27ds0dO2Q6duyomJgYlZeX++0vLy8/6zrT0tICOt6ypqz/lDlz5mjWrFn6+9//rmuuuSaUY4ZEoGvfs2ePvvrqKw0fPty3r76+XtJ/zy7u2rVLPXr0CO3QQdSU733nzp3Vrl07xcTE+Pb17t1bXq9XJ06cUFxcXEhnDpamrP3JJ5/U2LFjdd9990mS+vbtq2PHjmny5Ml6/PHHFR1t7v+jQXO2+zyPx9Omzp4sXrxY9913n5YuXRpxZ4wbYu4ntlOnTurVq9c5t7i4OGVlZamiokKlpaW+665du1b19fXKzMw86+1XVVXp9ttvV1xcnN5++23z/6uOi4vTgAEDVFxc7NtXX1+v4uJiZWVlNXidrKwsv+Mlac2aNWc93rKmrF+SZs+eraefflqrV6/2e55SJAl07b169dK2bdu0detW33bHHXfolltu0datW5WRkdGS4zdbU773N9xwg3bv3u0LM0n64osv1Llz54iJE6lpa//+++/PiJBToeZa+Vuutab7vKZatGiRJkyYoEWLFmnYsGHhHic4wv0s3ebIzc11/fv3d5s3b3YbN250V1xxhd/LjA8cOOCuvPJKt3nzZuecc5WVlS4zM9P17dvX7d69233zzTe+ra6uLlzLOK/Fixe7+Ph4N3/+fLdjxw43efJkl5SU5Lxer3POubFjx7qCggLf8R9++KGLjY11c+bMcTt37nSFhYUR/zLjQNY/a9YsFxcX5/7yl7/4fY+rq6vDtYQmC3Ttp4v0V/EEuv59+/a5hIQEN2XKFLdr1y63YsUKl5KS4p555plwLaHJAl17YWGhS0hIcIsWLXJ79+51f/vb31yPHj3c3XffHa4lNFl1dbXbsmWL27Jli5PknnvuObdlyxb3r3/9yznnXEFBgRs7dqzv+FMvM3744Yfdzp07XVFRUUS/zDjQ9S9YsMDFxsa6oqIiv/u8ioqKcC0hKCI6UA4fPuxGjx7tLrroIufxeNyECRP8fgmVlZU5Se799993zjn3/vvvO0kNbmVlZeFZRCPNnTvXdenSxcXFxblBgwa5TZs2+S676aab3Lhx4/yOX7JkievZs6eLi4tzV111lVu5cmULTxxcgaz/sssua/B7XFhY2PKDB0Gg3/v/FemB4lzg6//oo49cZmami4+Pd927d3e/+93vTP8H5FwCWXttba176qmnXI8ePVz79u1dRkaGe+CBB9y///3vlh+8mc52X31qvePGjXM33XTTGde59tprXVxcnOvevbt74403WnzuYAl0/TfddNM5j49UUc618nN/AAAg4ph7DgoAAACBAgAAzCFQAACAOQQKAAAwh0ABAADmECgAAMAcAgUAAJhDoAAAAHMIFAAAYA6BAgAAzCFQAACAOQQKAAAw5/8B7OHnpuQXzzMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(delta_gcs, bins=\"auto\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af4ccdfd-8c80-490c-b33d-a18edb44975f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2529e9b-d5e3-41a7-8012-a3d7ba513f78",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b2abf55-20ab-4482-a098-0e24af91d406",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce87a3d7-39a6-4698-bfb8-f7d6553461cd",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "scvi",
   "language": "python",
   "name": "scvi-env"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}