Simulations of covid-19 pool tests

pooltests.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "mighty-seating",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "incredible-handling",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext blackcellmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "german-china",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_persons = 9438\n",
    "n_pools = 701\n",
    "n_positive_pools = 101\n",
    "\n",
    "n_samples = 50\n",
    "positivities = [0.001, 0.005, 0.0075, 0.01, 0.02, 0.03, 0.05, 0.075, 0.1, 0.135, 0.2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "furnished-prospect",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD6CAYAAACh4jDWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXV0lEQVR4nO3dcayd9X3f8fendkIcUoihBnm2I9NhZQO0kHLlectUZXNabpuqZhNMN1KLN3lyhWiXTJM203/aTrIEU1dapIHkhQzD0oDnJMVqShfLNOoqMZObhNYxxOMuMLi1Z18CIWQVtCbf/XF+dz2+Offcc6997z3XvF/So/Oc73l+j3/PI9uf8/s9zzknVYUkST+y3B2QJA0HA0GSBBgIkqTGQJAkAQaCJKkxECRJwICBkORfJTme5JtJPpfkPUmuSHI4yfPtcW3X9nclmUhyIsnNXfWbkhxrr92XJK1+SZLHWv1oks0X/EglSX3NGQhJNgD/EhipqhuAVcAYsAc4UlVbgCPtOUmua69fD4wC9ydZ1Xb3ALAb2NKW0VbfBbxWVdcC9wL3zNWv0dHRAlxcXFxc5rfMatApo9XAmiSrgfcCJ4EdwP72+n7glra+A3i0qt6qqheACWBrkvXAZVX1VHU+DffwjDbT+zoIbJ8ePczmlVdeGbDrkqRBzBkIVfXnwG8CLwGngNer6svA1VV1qm1zCriqNdkAvNy1i8lW29DWZ9bPaVNVZ4HXgSsXdkiSpIUYZMpoLZ138NcAfwO4NMkv9GvSo1Z96v3azOzL7iTjScanpqb6d1ySNC+DTBl9DHihqqaq6q+ALwB/HzjdpoFoj2fa9pPApq72G+lMMU229Zn1c9q0aanLgVdndqSq9lXVSFWNrFu3brAjlCQNZJBAeAnYluS9bV5/O/AccAjY2bbZCTze1g8BY+3OoWvoXDx+uk0rvZFkW9vP7TPaTO/rVuDJ8lv3JGlJrZ5rg6o6muQg8HXgLPANYB/wPuBAkl10QuO2tv3xJAeAZ9v2d1bV2213dwAPAWuAJ9oC8CDwSJIJOiODsQtydJKkgWWlvhEfGRmp8fHx5e6GJK00s97B6SeVJUmAgSBJagwESRIwwEVlaTls3vOlvq+/ePfHl6gn0juHIwRJEmAgSJIaA0GSBBgIkqTGQJAkAQaCJKkxECRJgJ9D0ArmZxWkC8tA0JLzP3JpODllJEkCDARJUmMgSJIAA0GS1BgIkiRggEBI8sEkz3Qt30vyqSRXJDmc5Pn2uLarzV1JJpKcSHJzV/2mJMfaa/clSatfkuSxVj+aZPOiHK0kaVZzBkJVnaiqG6vqRuAm4C+ALwJ7gCNVtQU40p6T5DpgDLgeGAXuT7Kq7e4BYDewpS2jrb4LeK2qrgXuBe65IEcnSRrYfKeMtgP/q6r+N7AD2N/q+4Fb2voO4NGqequqXgAmgK1J1gOXVdVTVVXAwzPaTO/rILB9evQgSVoa8w2EMeBzbf3qqjoF0B6vavUNwMtdbSZbbUNbn1k/p01VnQVeB66c+Ycn2Z1kPMn41NTUPLsuSepn4EBI8m7g54H/OtemPWrVp96vzbmFqn1VNVJVI+vWrZujG5Kk+ZjPCOFngK9X1en2/HSbBqI9nmn1SWBTV7uNwMlW39ijfk6bJKuBy4FX59E3SdJ5mk8gfIK/ni4COATsbOs7gce76mPtzqFr6Fw8frpNK72RZFu7PnD7jDbT+7oVeLJdZ5AkLZGBvtwuyXuBnwJ+qat8N3AgyS7gJeA2gKo6nuQA8CxwFrizqt5ube4AHgLWAE+0BeBB4JEkE3RGBmPncUySpAUYKBCq6i+YcZG3qr5D566jXtvvBfb2qI8DN/Sov0kLFEnS8vCTypIkwECQJDUGgiQJMBAkSY0/oamLlj/VKc2PIwRJEuAIQYug3ztz35VLw8sRgiQJMBAkSY2BIEkCDARJUmMgSJIAA0GS1BgIkiTAQJAkNQaCJAkwECRJjYEgSQIGDIQk709yMMm3kjyX5O8luSLJ4STPt8e1XdvflWQiyYkkN3fVb0pyrL12X5K0+iVJHmv1o0k2X/AjlST1NegI4XeAP6yqvwV8CHgO2AMcqaotwJH2nCTXAWPA9cAocH+SVW0/DwC7gS1tGW31XcBrVXUtcC9wz3kelyRpnuYMhCSXAT8JPAhQVX9ZVd8FdgD722b7gVva+g7g0ap6q6peACaArUnWA5dV1VNVVcDDM9pM7+sgsH169CBJWhqDjBB+HJgC/nOSbyT5dJJLgaur6hRAe7yqbb8BeLmr/WSrbWjrM+vntKmqs8DrwJUzO5Jkd5LxJONTU1MDHqIkaRCDBMJq4CeAB6rqw8D/pU0PzaLXO/vqU+/X5txC1b6qGqmqkXXr1vXvtSRpXgYJhElgsqqOtucH6QTE6TYNRHs807X9pq72G4GTrb6xR/2cNklWA5cDr873YCRJCzdnIFTV/wFeTvLBVtoOPAscAna22k7g8bZ+CBhrdw5dQ+fi8dNtWumNJNva9YHbZ7SZ3tetwJPtOoMkaYkM+hOavwJ8Nsm7gW8D/5xOmBxIsgt4CbgNoKqOJzlAJzTOAndW1dttP3cADwFrgCfaAp0L1o8kmaAzMhg7z+OSJM3TQIFQVc8AIz1e2j7L9nuBvT3q48ANPepv0gJFkrQ8/KSyJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJGDwD6ZJF6XNe77U9/UX7/74EvVEWn6OECRJgIEgSWoMBEkSYCBIkhovKmtevAgrXbwcIUiSAANBktQYCJIkwECQJDUGgiQJGDAQkryY5FiSZ5KMt9oVSQ4neb49ru3a/q4kE0lOJLm5q35T289EkvuSpNUvSfJYqx9NsvkCH6ckaQ7zGSH8w6q6saqmf1t5D3CkqrYAR9pzklwHjAHXA6PA/UlWtTYPALuBLW0ZbfVdwGtVdS1wL3DPwg9JkrQQ5zNltAPY39b3A7d01R+tqreq6gVgAtiaZD1wWVU9VVUFPDyjzfS+DgLbp0cPkqSlMWggFPDlJF9LsrvVrq6qUwDt8apW3wC83NV2stU2tPWZ9XPaVNVZ4HXgypmdSLI7yXiS8ampqQG7LkkaxKCfVP5IVZ1MchVwOMm3+mzb65199an3a3NuoWofsA9gZGTkh16XJC3cQCOEqjrZHs8AXwS2AqfbNBDt8UzbfBLY1NV8I3Cy1Tf2qJ/TJslq4HLg1fkfjiRpoeYMhCSXJvnR6XXgp4FvAoeAnW2zncDjbf0QMNbuHLqGzsXjp9u00htJtrXrA7fPaDO9r1uBJ9t1BknSEhlkyuhq4IvtGu9q4Her6g+TfBU4kGQX8BJwG0BVHU9yAHgWOAvcWVVvt33dATwErAGeaAvAg8AjSSbojAzGLsCxSZLmYc5AqKpvAx/qUf8OsH2WNnuBvT3q48ANPepv0gJFkrQ8/KSyJAkwECRJjYEgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCTAQJEnNoF9/Lb1jbd7zpb6vv3j3x5eoJ9LicoQgSQIMBElSYyBIkgADQZLUGAiSJMBAkCQ1BoIkCZhHICRZleQbSX6/Pb8iyeEkz7fHtV3b3pVkIsmJJDd31W9Kcqy9dl/aDzUnuSTJY61+NMnmC3iMkqQBzGeE8Engua7ne4AjVbUFONKek+Q6YAy4HhgF7k+yqrV5ANgNbGnLaKvvAl6rqmuBe4F7FnQ0kqQFGygQkmwEPg58uqu8A9jf1vcDt3TVH62qt6rqBWAC2JpkPXBZVT1VVQU8PKPN9L4OAtunRw+SpKUx6Ajht4F/A/ygq3Z1VZ0CaI9XtfoG4OWu7SZbbUNbn1k/p01VnQVeB66c2Ykku5OMJxmfmpoasOuSpEHM+V1GSX4OOFNVX0vy0QH22eudffWp92tzbqFqH7APYGRk5Ide1/nxO3ukd7ZBvtzuI8DPJ/lZ4D3AZUn+C3A6yfqqOtWmg8607SeBTV3tNwInW31jj3p3m8kkq4HLgVcXeEySpAWYc8qoqu6qqo1VtZnOxeInq+oXgEPAzrbZTuDxtn4IGGt3Dl1D5+Lx021a6Y0k29r1gdtntJne163tz3AEIElL6Hy+/vpu4ECSXcBLwG0AVXU8yQHgWeAscGdVvd3a3AE8BKwBnmgLwIPAI0km6IwMxs6jX5KkBZhXIFTVV4CvtPXvANtn2W4vsLdHfRy4oUf9TVqgSJKWh59UliQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkSYCBIElqDARJEnB+X38tqen3a3P+0pxWCkcIkiTAQJAkNQaCJAkwECRJzZyBkOQ9SZ5O8qdJjif5jVa/IsnhJM+3x7Vdbe5KMpHkRJKbu+o3JTnWXrsvSVr9kiSPtfrRJJsX4VglSX0MMkJ4C/hHVfUh4EZgNMk2YA9wpKq2AEfac5JcB4wB1wOjwP1JVrV9PQDsBra0ZbTVdwGvVdW1wL3APed/aJKk+ZgzEKrj++3pu9pSwA5gf6vvB25p6zuAR6vqrap6AZgAtiZZD1xWVU9VVQEPz2gzva+DwPbp0YMkaWkMdA0hyaokzwBngMNVdRS4uqpOAbTHq9rmG4CXu5pPttqGtj6zfk6bqjoLvA5c2aMfu5OMJxmfmpoa6AAlSYMZKBCq6u2quhHYSOfd/g19Nu/1zr761Pu1mdmPfVU1UlUj69atm6PXkqT5mNddRlX1XeArdOb+T7dpINrjmbbZJLCpq9lG4GSrb+xRP6dNktXA5cCr8+mbJOn8DHKX0bok72/ra4CPAd8CDgE722Y7gcfb+iFgrN05dA2di8dPt2mlN5Jsa9cHbp/RZnpftwJPtusMkqQlMsh3Ga0H9rc7hX4EOFBVv5/kKeBAkl3AS8BtAFV1PMkB4FngLHBnVb3d9nUH8BCwBniiLQAPAo8kmaAzMhi7EAcnSRrcnIFQVX8GfLhH/TvA9lna7AX29qiPAz90/aGq3qQFiiRpefhJZUkSYCBIkhoDQZIE+AM57xj9fsAF/BEXSY4QJEmNgSBJAgwESVJjIEiSAANBktQYCJIkwECQJDUGgiQJMBAkSY2fVJaWgJ8U10rgCEGSBBgIkqTGQJAkAQaCJKmZMxCSbEryR0meS3I8ySdb/Yokh5M83x7XdrW5K8lEkhNJbu6q35TkWHvtviRp9UuSPNbqR5NsXoRjlST1McgI4Szwr6vqbwPbgDuTXAfsAY5U1RbgSHtOe20MuB4YBe5Psqrt6wFgN7ClLaOtvgt4raquBe4F7rkAxyZJmoc5A6GqTlXV19v6G8BzwAZgB7C/bbYfuKWt7wAeraq3quoFYALYmmQ9cFlVPVVVBTw8o830vg4C26dHD5KkpTGvawhtKufDwFHg6qo6BZ3QAK5qm20AXu5qNtlqG9r6zPo5barqLPA6cGWPP393kvEk41NTU/PpuiRpDgMHQpL3AZ8HPlVV3+u3aY9a9an3a3NuoWpfVY1U1ci6devm6rIkaR4GCoQk76ITBp+tqi+08uk2DUR7PNPqk8CmruYbgZOtvrFH/Zw2SVYDlwOvzvdgJEkLN8hdRgEeBJ6rqt/qeukQsLOt7wQe76qPtTuHrqFz8fjpNq30RpJtbZ+3z2gzva9bgSfbdQZJ0hIZ5LuMPgL8InAsyTOt9qvA3cCBJLuAl4DbAKrqeJIDwLN07lC6s6rebu3uAB4C1gBPtAU6gfNIkgk6I4Ox8zssSdJ8zRkIVfUn9J7jB9g+S5u9wN4e9XHghh71N2mBIklaHn5SWZIEGAiSpMZAkCQBBoIkqTEQJEmAP6F50ej3E43+PKOkQRgI0pDwd5e13JwykiQBBoIkqTEQJEmAgSBJagwESRJgIEiSGgNBkgQYCJKkxkCQJAEGgiSpMRAkScAAgZDkM0nOJPlmV+2KJIeTPN8e13a9dleSiSQnktzcVb8pybH22n1J0uqXJHms1Y8m2XyBj1GSNIBBRggPAaMzanuAI1W1BTjSnpPkOmAMuL61uT/JqtbmAWA3sKUt0/vcBbxWVdcC9wL3LPRgJEkLN2cgVNUfA6/OKO8A9rf1/cAtXfVHq+qtqnoBmAC2JlkPXFZVT1VVAQ/PaDO9r4PA9unRgyRp6Sz066+vrqpTAFV1KslVrb4B+B9d20222l+19Zn16TYvt32dTfI6cCXwysw/NMluOqMMPvCBDyyw69LK5VdkazFd6IvKvd7ZV596vzY/XKzaV1UjVTWybt26BXZRktTLQgPhdJsGoj2eafVJYFPXdhuBk62+sUf9nDZJVgOX88NTVJKkRbbQQDgE7GzrO4HHu+pj7c6ha+hcPH66TS+9kWRbuz5w+4w20/u6FXiyXWeQJC2hOa8hJPkc8FHgx5JMAr8G3A0cSLILeAm4DaCqjic5ADwLnAXurKq3267uoHPH0hrgibYAPAg8kmSCzshg7IIcmSRpXuYMhKr6xCwvbZ9l+73A3h71ceCGHvU3aYGi3ryQKGkp+EllSRJgIEiSGgNBkgQYCJKkZqGfVJY0hLwBQefDEYIkCTAQJEmNgSBJAgwESVJjIEiSAANBktR42+ky8zZBLTX/zmk2jhAkSYCBIElqDARJEuA1BEk9eJ3hnckRgiQJGKIRQpJR4HeAVcCnq+ruZe7SefNdlqSVZCgCIckq4D8CPwVMAl9Ncqiqnl3enkmaTb83PL7ZWZmGIhCArcBEVX0bIMmjwA7AQJBWKEfIK8+wBMIG4OWu55PA312mvgxsrr/wkvobJDQMlqWTqlruPpDkNuDmqvoX7fkvAlur6ldmbLcb2N2efhA40WN3Pwa8sojdXQz2eWmstD6vtP6CfV4q59PnV6pqtNcLwzJCmAQ2dT3fCJycuVFV7QP29dtRkvGqGrmw3Vtc9nlprLQ+r7T+gn1eKovV52G57fSrwJYk1yR5NzAGHFrmPknSO8pQjBCq6mySXwb+G53bTj9TVceXuVuS9I4yFIEAUFV/APzBBdhV3ymlIWWfl8ZK6/NK6y/Y56WyKH0eiovKkqTlNyzXECRJy+yiCYQko0lOJJlIsme5+zOIJC8mOZbkmSTjy92fXpJ8JsmZJN/sql2R5HCS59vj2uXs40yz9PnXk/x5O9fPJPnZ5ezjTEk2JfmjJM8lOZ7kk60+tOe6T5+H9lwneU+Sp5P8aevzb7T6UJ7nPv1dlHN8UUwZta+++J90ffUF8Ilh/+qLJC8CI1U1tPdAJ/lJ4PvAw1V1Q6v9e+DVqrq7he/aqvq3y9nPbrP0+deB71fVby5n32aTZD2wvqq+nuRHga8BtwD/jCE91336/E8Z0nOdJMClVfX9JO8C/gT4JPBPGMLz3Ke/oyzCOb5YRgj//6svquovgemvvtB5qqo/Bl6dUd4B7G/r++n8JzA0ZunzUKuqU1X19bb+BvAcnU/wD+257tPnoVUd329P39WWYkjPc5/+LoqLJRB6ffXFUP/FbAr4cpKvtU9hrxRXV9Up6PynAFy1zP0Z1C8n+bM2pTQUUwK9JNkMfBg4ygo51zP6DEN8rpOsSvIMcAY4XFVDfZ5n6S8swjm+WAIhPWorYS7sI1X1E8DPAHe2qQ4tjgeAvwncCJwC/sOy9mYWSd4HfB74VFV9b7n7M4gefR7qc11Vb1fVjXS+EWFrkhuWuUt9zdLfRTnHF0sgDPTVF8Omqk62xzPAF+lMfa0Ep9v88fQ88pll7s+cqup0+4f1A+A/MYTnus0Rfx74bFV9oZWH+lz36vNKONcAVfVd4Ct05uOH+jzDuf1drHN8sQTCivvqiySXtgtxJLkU+Gngm/1bDY1DwM62vhN4fBn7MpDpf+zNP2bIznW7ePgg8FxV/VbXS0N7rmfr8zCf6yTrkry/ra8BPgZ8iyE9z7P1d7HO8UVxlxFAu+3qt/nrr77Yu7w96i/Jj9MZFUDnE+O/O4x9TvI54KN0vl3xNPBrwO8BB4APAC8Bt1XV0FzEnaXPH6UzvC7gReCXpueMh0GSfwD8d+AY8INW/lU6c/JDea779PkTDOm5TvJ36Fw0XkXnDfGBqvp3Sa5kCM9zn/4+wiKc44smECRJ5+dimTKSJJ0nA0GSBBgIkqTGQJAkAQaCJKkxECRJgIEgSWoMBEkSAP8P1SRc4pUqqR0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "persons_per_pool = []\n",
    "for _ in range(1000):\n",
    "    _, sample = np.unique(\n",
    "        np.random.choice(range(n_pools), size=n_persons), return_counts=1\n",
    "    )\n",
    "    persons_per_pool += list(sample)\n",
    "\n",
    "plt.hist(persons_per_pool, bins=max(persons_per_pool), rwidth=0.9)\n",
    "sns.despine()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "norman-reporter",
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = pd.DataFrame()\n",
    "for positivity, _ in product(positivities, range(n_samples)):\n",
    "    sample = pd.DataFrame(\n",
    "        {\n",
    "            \"person\": np.random.uniform(size=n_persons) <= positivity,\n",
    "            \"pool\": np.random.choice(range(n_pools), size=n_persons),\n",
    "        }\n",
    "    )\n",
    "    n_positive_pools_sample = (\n",
    "        sample.groupby(\"pool\").apply(lambda pool: pool[\"person\"].any()).sum()\n",
    "    )\n",
    "    samples = samples.append(\n",
    "        pd.Series(\n",
    "            {\"positivity\": positivity, \"n_positive_pools\": n_positive_pools_sample}\n",
    "        ),\n",
    "        ignore_index=True,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "previous-installation",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWMAAADGCAYAAAAHZOEIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5cElEQVR4nO2deXhV1dX/PyvzHJIQAgQQEJkJo4jVCg51QOtYW6x1bn191apttc6tHXyr1V+1Dq21jq2IY1XqVAXFeSggMkNAwkwGIHNucof1++OcJDfJTXIDuUkg6/M898nd55y9zzon537vumvvvbaoKoZhGEb3EtXdBhiGYRgmxoZhGD0CE2PDMIwegImxYRhGD8DE2DAMowdgYmwYhtEDMDE2jA4iIpUiMry77egORGSWiGzrbjsORkyMD2JEZIyIvCciZSKyQUTOauW4X4uIisgJQduuE5FvRKRcRHaIyH0iEhO0f5KIfOS2vU1EftUV19QTUNUUVf2mo/VE5NuukAe/VETOCTrmZyKyy72vT4hIfNC+q0VksYjUishTnXQ5Rg/BxPggxRXO14DXgUzgcuAZERnZ7LhDge8BO5s18W9giqqmAeOBicA1QfufBT50254J/K+InB6BS2lC8BfCgYaqfuQKeYqqpgCnAZXA2wAichJwE3A8MBQYDvwmqIkdwO+BJ7rSbqNrMDE+eBkNDATuU1W/qr4HfAJc0Oy4h4Abgbrgjaq6UVVL3aIAAWBE0CFDgblu2xuBj4FxoQwRkTtE5CUReV5EKkRkqYhMDNo/UEReFpFiEdkkIteEqPuMiJQDF4vIdNdDLBeRQhH5U9Dxp4vIKhEpFZFFIjImaF+BiFwvIstdz/N5EUlw9/UVkdfdentcrz/k58P1Zke4758SkYdF5A332r5wv+DC4SLgJVWtCio/rqqrVHUv8Dvg4vqDVfVfqvoqsLu9hkXkYhH5REQedK91rYgcH7R/oIjMd691g4j8JGhfvIjc7/4i2uG+j2/lPDeKyHb32tcFn8PoGCbGBy/SyrbxDQWRc4E6VX0zZAMiP3QFsATHM/5b0O77gQtFJFZERgFHAgvasOcM4EUcT/pZ4FW3bhSOF/41kIvjFV7neonBdV8C+gBzgT8Df3a99kOBF1x7RwLzgOuAbOBN4N8iEhfU1veBk4FhQB6NYvcLYJtbLwe4BQg3V8B5OB5sBrABuLO9CiKShPOL5OmgzeNw7kM9XwM5IpIVph3NOQL4BugL/Br4l4hkuvvm4VzvQNeO/wsS0luBGcAknP/7dOC2ENcwCrgaOFxVU4GTgIJ9tLXXY2J88LIWKAJucEXvRJxwQhKAiKQA/4cjXCFR1WddwRsJPAIUBu1+HedDXOOe63FV/W8b9ixR1ZdU1Qv8CUjA+cAfDmSr6m9Vtc6Nxf4dmBNU9zNVfVVVA6paA3iBESLSV1UrVfVz97gfAG+o6rvuee4FEoFvBbX1gKruUNU9OF8Ck9ztXmAAcIiqet2QQrhi/C9V/VJVfThfFpPaOR7gHJwvuQ+CtqUAZUHl+vepYdrRnCLgfvd6ngfWAaeKyGDgaOBGVfWo6jLgMRp/NZ0P/FZVi1S1GOeLpvkvKgA/EA+MFZFYVS1wfyUZ+4CJ8UGKK0ZnAqcCu3A8vxdwvCFwPmD/VNVNYbSVD6wC/gLgeldvA7/FEdXBwEkicmUbzWwNai9Ao1d2CDDQDQ+UikgpjleaE6quy2U4XxBrReS/InKau30gsLnZebbieNz17Ap6X40jgAD34Hi174jTcXlTG9fSnNbabIuLgH80E/xKIC2oXP++ogO2BLO9Wfubce7RQGCPqlY021d/n5rcx6B6TVDVDThf5ncARSLynIi0OM4IDxPjgxhVXa6qM1U1S1VPwukQ+tLdfTxwjdtzvwtHUF8QkRtbaS4GJySA245fVf+hqj5V3QY8B8xuw5zB9W/c0MQgnA6prcAmVe0T9EpV1eC2mnioqpqvqucB/YC7gZdEJNlt75Cg84h73u1t2FXfZoWq/kJVhwPfBX4eqfin65nOAv7RbNcqnLBAPROBQlVtN0bcCrnuPahnCM492gFkikhqs33196nJfQyq1wL319PR7vGK8/8w9gET44MYEckTkQQRSRKR63F+hj/l7j4eJ348yX3tAP4HeNit+2MR6ee+HwvcDCx06653NssPRSRKRPrjhAiC453NmSoiZ7ujIa4DaoHPcb4cyt2OoEQRiRaR8SJyeBvX9SMRyXY931J3sx/H8z9VRI4XkVicXwO1wKdh3KvTRGSEK17lbnv+9urtIxcAn4b4Sf8P4DIRGSsiGThx2qeCbIxxOxyjgWj3f9vW6JJ+OF+4sW7/wBjgTVXdinNP/uC2kYfza2OuW28ecJuIZItIX+BXwDPNGxeRUSJynNu558EJWUXqnh30mBgf3FyAM2StCEd8v6OqtQCqultVd9W/cD5Ee1W10q17FLBCRKpwOsLexAkfoKrlwNnAz4C9wDJgJW13XL2GI9h7XbvOdmOZfhxPdBKwCSeO+hiQ3kZbJwOrRKQSpzNvjhv7XAf8CHjQbee7wHdVta71pho4DKcDshL4DPiLqi4Ko96+cCFNO+4AUNW3gT8C7+OEBjbjdLzVcxuO4N2Ec501hOhYC+ILnOsqwfnffC/Iyz4PZ0TMDuAV4Neq+q677/fAYmA5sAJY6m5rTjxwl9v+Lhzxv6UNe4w2EEsub0QaEbkDGKGqP+puW3oLInIx8GM3hGAcAJhnbBiG0QMwMTYMw+gBWJjCMAyjB2CesWEYRg/AxNgwDKMHcMBmwAI4+eST9e233+5uMwzDMMIlVM4Y4AD3jEtKSrrbBMMwjE7hgBZjwzCMgwUTY8MwjB6AibFhGEYPwMTYMAyjB3BAj6YwDMOINIGAUucPUOsN4PH5qfMFGJyZ1OnnMTE2DKNXo+qIrccboNbnp9YboNbnvPd4A3j9Afx+ZXdVHYXlHgorPEwc1IfvTuzcPPomxoZhHPR4/QE8Xr8rsgFqvf4G8a3zBQgo+ANKSWUtu8o8FJZ72FXe+LeovBZfoDF1xPRhmSbGhmEYzfEHNEhs/Q0hhXov1+8KqdcfoKiilsKyRrF1XrUUV9TiD8rVEx8TRU5aAoP6JDF1SAY56Qn0T0ugf3oCs8cP6PRrMDE2DKPHEwhoo9D6grxc17v1+htFtNbnp6i8ttGzDfJ0d1fWNVnDKzE2mv7pCQzrm8yM4ZnkpDmCm5OeQJ/EWJquWuUgAlFRrU6k22dMjA3D6HZUtTGE0ODR1ocSAtT5Ak2Or67zUVhe2xhOcD3doopa9lQ1XdglJT6G/ukJjOqfRv+0+CaCmxofE1JwmxMlEB8bTUJsFImx0Z167fWYGBuG0SXU+YIFttGzrR+h0Dybb6XH1ySU0BjDraW8xtvk2D6JseSkJTAhN90VW0d0c9ISSI4PT+ZEIMEV3ISY6Mb3sdHEx0SFJdr7g4mxYRidgs8faBpCCBZfr59AM7FVVco9vhahhPq/VbVN1zbNTI6jf1oCU4dkOGLrxnBz0hJICNNbFXFiwYlx0d0iuG0RUTEWkT44i0uOx1nG+1JgHfA8zmKIBcD3VXWve/zNOKvU+oFrVPU/kbTPMIzwCQS0SadY8/itz99yoYqAKqXV3iahhOBOsxpvo+CKQHaK49EeObyvG0qIp39aAv1SE4iLCW+OWr3gJsRGkxjbswS3LSLtGf8ZeFtVvycicUASzuqxC1X1LhG5CWel2xvd5eDnAOOAgcACERnprh5sGEaEaYjbekN0lPn81PlCrwoUCCi7q2opDNFpVlheS52/Md4bLUI/N4TQPIabnRpPTHTHBTchKJbb0wW3LSImxiKSBhwDXAzgLpdeJyJnALPcw54GFgE3AmcAz7lLyW8SkQ3AdJxl0w3D6ASCY7XN47eh4rb1+AIBSirqmoQRCl3BLapoOgY3Nlrol+oMAZswqE8Twc1KiSc6zJEIoQQ3ISaaxLgDV3DbIpKe8XCgGHhSRCYCS4BrgRxV3QmgqjtFpJ97fC7weVD9be42wzDCxB9Qarx+Z1KDO7kheKJD87htMHW+AMUVtU0F130VV9Q2qRsfE0X/tAQGZSYxbWim21nmhBQykuOIClMoRSAuplFkgzvP4mOiIjKErKcSSTGOAaYAP1XVL0TkzzghidYIdddbPDoicjlwOcCQIUM6w07DOGCp8wWo8Hip8Pio8PioqvO16t0CeLz+FpMe6v82H4ObFBdNTloCw/umODHc9EYPN72VMbihCBbchNimnWe9TXDbIpJivA3YpqpfuOWXcMS4UEQGuF7xAKAo6PjBQfUHATuaN6qqjwKPAkybNs2WtjZ6FR6vn/Ig8a2pa9mlUj8GN9QIhdLqpkPCUhNi6J+WwOj+aY7QpjcOC0sJcwwutBTc+tBCoglu2ERMjFV1l4hsFZFRqroOOB5Y7b4uAu5y/77mVpkPPCsif8LpwDsM+DJS9hnGgUB1nSO65TVeyj2+FpMfAEqr61i5o5wV20pZuaO8xaSHPkmx5KQmMHFQn30eg1tPXEx9R1nLWK4J7v4R6dEUPwXmuiMpvgEuwcmh/IKIXAZsAc4FUNVVIvICjlj7gKtsJIXRm1BVKmtd8fV4qfT4mkzzrcfj9bN2VzkrtpezYnsZW/dUA85Ms3ED0xieneKOv43v0BjceuJiopqMTjDB7RpE2wow9XCmTZumixcv7m4zDGOf8AeUSld4Kzw+Kmt9DQltggkElG9Kqli5vYwV28tYX1iBL6DERgujclIZn5vOhNx0hvZNDrvjrF5wg8MJJrhdQqs312bgGUYX4fUH3Fhvo/i25gsVlntYvq2MlTvKWLWjrGE22iFZSZw8vj8TctMZ1T+V+JjWvd7mgpvgzjyLj4kOe3iZ0XWYGBtGhKj1+RvivRUeH9UhOtvqqfT4WLmjrMH7LaqoBSArOY7DD8lkwqB0xg1MJz0xNmT9KIH0pFgyk+JIjo8hIdYE90DDxNgwOomaOj8VHqejrdzjpdbbsrOtHq8/wLpdFazcUcaKbWVsKqlCcVI6jhuYxqkTBjAhN53+6QmtjmiIEuiTFEdmchwZSbFhz14zeiYmxoaxD6gqVa741oceWpsuDE6Ohq17qlnher5rd1ZQ5w8QLcKIfimcM3UQE3LTOTQ7pU2PNjpKyEiKJTM5jj5Jceb9HkSYGBtGGAQCSkWtr8kEi1CdbcHsrqxlxXY39LCjvCHtY26fRI4b3Y8JuemMGZBGYlzbox1iousFOJ4+ibHWwXaQYmJsGCHw+QNU1voor3FCDlW1vjanEoMzJnj1znJWuB1vO0o9AKQnxjLBHfEwITedzOS4ds8fGy1kJMeRlRxHWoIJcG/AxNgwXDxeP8UVteytrqO6zt/mtGJwkudsLKpixfZSVm4vJ7+ogoA6eRtG90/luFE5TBiUzuCMxLBmssXFCBlJcWQlx5OWGP7sN+PgwMTY6NWom2+3sMKZKtyWAKsqO8o8DZ7v6h3l1Hj9iMDwvsmcPnEgE3LTOSwnldgwO9PiYqLISo4jMyUu7CWAjIMTE2OjV1Lrc7zgooraNkc9lNV4G4abrdhe1jDVOCctnqNGZDE+1xlyltKBacUJsVFkJceTkRxLakLooWpG78PE2OhVlFbXUVjuhCJCecG1Pj9rd1Y0iO+WZlON6+O+/dISOnTexLhoxwNOjutwPgijd2BPhXHQ4/UHKKqopajcg6cVL/ib4koWri3ikw0l1PoCxEQJo/qnMufwwYzPTWdYVnKHO9GS46PJdAU4Kc4+akbb2BNiHLSU1XgpKvewp6ou5EgIj9fPpxt3s3BNId+UVBEXHcWRh2Zx5PAsRg9oe6pxa6TEx5CZ4oyC6GiCHqN3064Yi8ihOHmJa0VkFpAH/ENVSyNrmmF0HJ8/QHFlLUXlta1OP968u4oFaxwvuMbrZ3BGIhd/ayhHj+i7TyGE1IQYslLiyEgyATb2nXCevJeBaSIyAngcN+8wMDuShhlGR6jweCksr2VPVV3IyRi1Pj+ff7ObhWuKyC+qJDZamDEsixPG5nBYv5QOjWIQcQU4OZ7M5LiwVy02jLYIR4wDquoTkbOA+1X1QRH5KtKGGUZ7+ANKiesFV9b6Qh6zdU81C9cW8XF+MVV1fgamJ3DBjEM45rBsUhLC94JFnMkbWclxZCTHhT10zTDCJZyn0Ssi5+GsyvFdd5uNxzG6japan7NmW1UdvhDJ1+t8Ab7Y5HjB6woriIkSpg/L5PgxOYzpnxq2F1yfiCcj2cmGZol4jEgSjhhfAlwB3Kmqm0RkGPBMZM0yjKYEAkpJleMFV3hCe8E7SmtYuKaQD/NLqKz10T8tgfOPGMIxh2WT1krqyeZERwl93EQ8GZaIx+hC2hVjVV0NXBNU3oSzfp1hRJyaOj+F5R5KKmtDLkHk9Qf4b8EeFq4pYvXOcqJFmDY0gxPG5DB2YFpYK1+IOOvE1ceATYCN7qBVMRaRFUCoyaECqKrmRcwqo1ejquyuqqOw3EN5TWgvuLDcw8I1hXywvphyj4/slHh+cPhgZo3Mpk9S+4l4wOmEy051BNhiwEZ305ZnfFqXWWEYOON+i8prKa70hMwN7AsEWFKwlwVri1i5vYwogamHZHD8aCchTzhecHJ8NFkp8fRNidunccSGESlaFWNV3Vz/XkRygMPd4peqWhRpw4zeQ1Wtj617q1tN1FNc4eG9tUUsWldMaY2XrOQ4zp06iFmj+oWVjjIhNoq+KfH0TYlvN3ewYXQX4Uz6+D5wD7AIJ0TxoIjcoKovRdg24yDH4/WzbW8NJZW1LUTYH1CWbtnLwjWFLN9WBgKTB2dw/Jh+TBrUp92pyXExQlZyPFkpcZaMxzggCGc0xa3A4fXesIhkAwsAE2Njn/D5A2wvrWFXmafFNOXdlbW8t87xgvdU1ZGRFMtZU3I5blQ/slLi22zXWREjjuwUywdsHHiEI8ZRzcISuwHr7TA6TCCg7Cr3sL20psn4YFXlq62lLFxTyFdbS0Ehb1A6F39rKFOGZLQ5uiFKaFgRIyMpzlbEMA5YwhHjt0XkP8A8t/wD4M3ImWQcbKgqxZW1bNtb0yJ38Nqd5TzzxWY2FleRnhjLGRMHcuyofm2mqBSBtIRY+qbG2WQM46Ch3adYVW8A/oaTIGgi8Kiq3hhpw4yDg9LqOlZsL2NjUVUTId5ZWsP/e2cdv3l9NXuq6vifY4bz0A8n84PDh7QqxKkJMQzrm8yUIRmMHZhGv9QEE+JuRES44IILGso+n4/s7GxOO61jA7FmzZrF4sWLAZg9ezalpaWdaeYBQ7iT8z8BvDjjjr+MnDnGwUJlrY/Nu6tajBMur/Hy8tJtLFxTRGyMcO7UQZyaN6DVYWZJcdFkpcTRNyXeMqL1MJKTk1m5ciU1NTUkJiby7rvvkpubu19tvvlm7/3R3a5b4Y6m+BL4HvB94AsR+V6kDTMOTDxeP/mFFazYVtZEiOt8AeYv2851zy9jwZpCjh2dzX3fn8TZUwa1EOL42Chy+ySSNyidiYP7MCgjyYS4h3LKKafwxhtvADBv3jzOO++8hn1VVVVceumlHH744UyePJnXXnsNgJqaGubMmUNeXh4/+MEPqKmpaagzdOhQSkpKADjzzDOZOnUq48aN49FHH204JiUlhVtvvZWJEycyY8YMCgsLu+JSI044v/HqR1NcpKoXAtOB2yNrlnGg4fUH2FRSxbKtpZRU1jVsD6jy8YYSfvHiMub9dytjBqRy9zl5XHb08CYz5WKjhZy0eMblpjFlSAZDspJseaIDgDlz5vDcc8/h8XhYvnw5RxxxRMO+O++8k+OOO47//ve/vP/++9xwww1UVVXx17/+laSkJJYvX86tt97KkiVLQrb9xBNPsGTJEhYvXswDDzzA7t27AUfkZ8yYwddff80xxxzD3//+9y651kgT8dEUIhINLAa2q+ppIpIJPA8MBQqA76vqXvfYm4HLAD9wjar+J9zzGN2DP6DsLKthZ5mnRQa11TvKeOaLLWwqqWJY32SumHko4wamNzlGBAakJzAoI8lyQuwjv/n3KlbvKO/UNscOTOPX3x3X7nF5eXkUFBQwb948Zs9umuL8nXfeYf78+dx7770AeDwetmzZwocffsg111zTUD8vL3RmhQceeIBXXnkFgK1bt5Kfn09WVhZxcXENcempU6fy7rvv7vN19iS6YjTFtcAaIM0t3wQsVNW7ROQmt3yjiIwF5gDjgIHAAhEZqaqhl2swuhVVpbiilq17a6jzNR0hsX1vDc9+uYWlW/aSlRzHlbMO5agRfVtMV05NiGF4drKtD3eAc/rpp3P99dezaNGiBu8VnGfk5ZdfZtSoUS3qtDcGfNGiRSxYsIDPPvuMpKQkZs2ahcfjASA2NrahfnR0ND5f6PwlBxrhZG27QUTOBo7GmYH3qKq+Ek7jIjIIOBW4E/i5u/kMYJb7/mmcmX03utufU9VaYJOIbMAJiXwW7sUYXcOeqjq27KmmptmyRmU1Xl5aso331hYSHxPNnMMHc8r4AS1WwoiNFoZkJdEvtWMrLBuhCceDjSSXXnop6enpTJgwgUWLFjVsP+mkk3jwwQd58MEHERG++uorJk+ezDHHHMPcuXM59thjWblyJcuXL2/RZllZGRkZGSQlJbF27Vo+//zzLryi7iFcl+RTnNBBAPhvB9q/H/glkBq0LUdVdwKo6k4R6eduzwWC7/g2d5vRQyj3eNmyu7pFPuFan5+3Vuxi/tc7qPX5OWFMDudMGdQih7AI9EuNZ0hmkg1JO4gYNGgQ1157bYvtt99+O9dddx15eXmoKkOHDuX111/nf//3f7nkkkvIy8tj0qRJTJ8+vUXdk08+mUceeYS8vDxGjRrFjBkzuuJSuhXRUJlZgg8Q+THwK+A9HM94JvBbVX2inXqnAbNV9Up3IdPr3Zhxqar2CTpur6pmiMjDwGeq+oy7/XHgTVV9uVm7lwOXAwwZMmTq5s2bMSJLTZ2fLXuq2VNV12R7QJWP8kt4YfFW9lTVMe2QDM6bPoSBfRJbtJESH8Ow7GRSrFPO6N20Gp8J55NxAzBZVXcDiEgWjqfcphgDRwGni8hsIAFIE5FngEIRGeB6xQOA+s7BbcDgoPqDgB3NG1XVR4FHAaZNm9b2N4mxX9T5AmzbW01RRctEPiu3l/HMF5vZvLuaQ7OTufrYEYwZkNaijZhoYXBGEjlp8ZYrwjDaIBwx3gZUBJUrgK3tVVLVm4GbAYI84x+JyD046+nd5f59za0yH3hWRP6E04F3GDbBpFvwB5Qdpc4IieYrLW/dU82zX25h2dZS+qbEcfWxIzjy0KyQuYSzU+MYkplsqycbRhiEI8bbcSZ6vIYzA+8M4EsR+TmAqv6pg+e8C3hBRC4DtgDnuu2sEpEXgNWAD7jKRlJ0LapKYXkt2/ZWt1jiqLS6jheXbOP9dUUkxkZz/hFDOHFs/5BCmxQXzdC+yaSHue6cYRjhxYx/3dZ+Vf1Np1rUAaZNm6b1c9qN/aOkspate6rxNEvk4/H6eWPFTv799Q58fuU743I4e3JuyBzB0VFCbkYiA9MTLCRhGKHZ95hxd4qtEXnKapwREpW1TUdIBALKB/nFvLB4K6XVXqYPy2TO4YMZkN6ycw4gMzmOoX2TbCkjw9hHrGu7lxIIKJv3VLOrzNNi3/JtpTzzxRa27qlmRL8Urjt+JKP6p4ZoxVnSaGhWMhlhLH9kGEbrmBj3QpxkPpUtvOHNu6t49ostLN9eRr/UeK49/jCOGJYZMuQQJTCwTyK5fRItoXsvJTo6mgkTJuDz+Rg2bBj//Oc/6dOnT3ebdcDSbsy4J2Mx445TUlnLN8VVTUZJ7Kmq48XFW/lgfTFJ8dGcPXkQ3xmb0+ry9X2SYhnWN9kyqfUgPtu4u/2DOsCRh2a1e0xKSgqVlZUAXHTRRYwcOZJbb721U+04CNn3mLGIjAT+ijNzbryI5AGnq+rvO9FAI8IEAsqm3VUUldc2bPN4/fx7+Q7eWL4TX0A5ZcIAzpqUS0pC6MciLiaKoVlJ7a5FZ/Q+jjzyyIZpzV9++SXXXXddQ57jJ598klGjRvHUU08xf/58qqur2bhxI2eddRZ//OMfAXj88ce5++67GThwIIcddhjx8fE89NBDFBcXc8UVV7BlyxYA7r//fo466qhuu85IEk6Y4u84Ez/+BqCqy0XkWcDE+AChus5HfmEl1UG5JDbvruKBhfnsKPMwY3gmcw4fQk4rK2xYZjWjLfx+PwsXLuSyyy4DYPTo0Xz44YfExMSwYMECbrnlFl5+2ZlIu2zZMr766ivi4+MZNWoUP/3pT4mOjuZ3v/sdS5cuJTU1leOOO46JEycCcO211/Kzn/2Mo48+mi1btnDSSSexZs2abrvWSBKOGCep6pfN4oYHR5qkXkBRuYeC3dUNYQlV5f11xTz16SaS42K47dQxLdJaBmOZ1YzWqKmpYdKkSRQUFDB16lS+853vAE6Sn4suuoj8/HxEBK/X21Dn+OOPJz3ded7Gjh3L5s2bKSkpYebMmWRmZgJw7rnnsn79egAWLFjA6tWrG+qXl5dTUVFBamroDuUDmXA+YSUicijOhA/cVT52RtQqY7/xB5RNJZUUVzTmk/B4/Tz+8SY+3lDC+Nx0rpp1aJME78FYZjWjPRITE1m2bBllZWWcdtppPPzww1xzzTXcfvvtHHvssbzyyisUFBQwa9ashjrx8Y0hrvr0l231WwUCAT777DMSE0MPqTyYCGee6lU4IYrRIrIduA64IpJGGftHVa2P5dtKmwjxlj3V3PrKCj7ZWMK5Uwdx88mjQwqxCOSkxTNxcB8TYiMs0tPTeeCBB7j33nvxer2UlZU1rIX31FNPtVt/+vTpfPDBB+zduxefz9cQ0gA48cQTeeihhxrKy5Yt62zzewzhiPFmVT0ByAZGq+rRqmqp0noou8o8rNxe1jCTTlV5f20Rt726guo6P7fNHsPZUwaFHI6WEh/DuIFpDM9OaXUkhWGEYvLkyUycOJHnnnuOX/7yl9x8880cddRR+P3tZzTIzc3llltu4YgjjuCEE05g7NixDaGMBx54gMWLF5OXl8fYsWN55JFHIn0p3UY406G3AG/jLJX0nvagsXA2tK0Rnz/AxuKqJmkuww1LWGY1o7uprKwkJSUFn8/HWWedxaWXXspZZ53V3WZFgv1KoTkK+C5OuOJxEXkdZ0WOjzvJOGM/qfB4yS+qpDYor8TWPdXcv3A9O8s8nDt1EGdOyg3pDVtmNaMncMcdd7BgwQI8Hg8nnngiZ555Zneb1OV0aNKHiGQAfwbOV9VuH/FvnjFsL61h657qhnzDqsqi9cU89UkBSXHRXH3ciJCjJaKjhOHZyfS1McOG0ZXsl2eMiMzEWYj0FJxll77fOXYZ+4rXH2BDUSWl1Y3DhjxeP098vImPNpQwfmAaVx07ImRYIjk+mpE5qTaDzjB6EOHMwNsELANeAG5Q1apIG2W0TVmNlw1FlU1WZW4IS5R6+N7UQZzVSlgiJy2eoVnJlk/CMHoY4XjGE1W1POKWGO2iqmzbW8P20pomYYkP1hfz5CcFJMZFc8vsMYzPbRmWiIkWhvW1sIRh9FRaFWMR+aWq/hG4U0RaBJZV9ZqIWmY0odbnZ0NRJeU1jZMfPV4/T3yyiY/y2w5LpMTHcFhOioUlDKMH01YXev0E8MXAkhAvo4sora5jxbayJkK8dU81t726ko/zS/je1EHcfMqYkELcPz2BcQPTTIiNiPDKK68gIqxdu7ZhW3FxMUcccQSTJ0/mo48+4i9/+ct+n2fo0KGUlJSEffxTTz3F1VdfDcAjjzzCP/7xjxbHFBQUMH78+P22rbNoVYxV9d/u22pVfTr4BVR3jXm9G1Vly+5q1uysaLIm3aJ1Rdz26koqa33cMnsM54SYxBETLYzMSWFYX4sPG5Fj3rx5HH300Tz33HMN2xYuXMjo0aP56quvGDx4cIfFWFUJBALtHxgmV1xxBRdeeGGntRcpwhlcenOY24xOxOP1s2pHOdtLa5ps++uiDfztw28Y0S+Fu86eEDI+nBIfw4TcdEt1aUSUyspKPvnkEx5//PEGMV62bBm//OUvefPNN5k0aRI33ngjGzduZNKkSdxwww0A3HPPPRx++OHk5eXx6187S2wWFBQwZswYrrzySqZMmcLWrS0XoH/wwQeZMmUKEyZMaPDE9+zZw5lnnkleXh4zZsxoSOMZzB133MG9994LwJIlS5g4cSJHHnkkDz/8cMMxBQUFfPvb32bKlClMmTKFTz/9FIALLriA1157reG4888/n/nz57Nq1SqmT5/OpEmTyMvLIz8/f7/vZ1sx41OA2UCuiDwQtCsNy9oWUfZU1bGxuBJfkDe8dU81f16Yz47SGs6ZksvZk0NPaR6QnsCQzCTzho2I8+qrr3LyySczcuRIMjMzWbp0KVOmTOG3v/0tixcv5qGHHqKgoIBVq1Y15JR45513yM/P58svv0RVOf300/nwww8ZMmQI69at48knn2zVk+7bty9Lly7lL3/5C/feey+PPfYYv/71r5k8eTKvvvoq7733HhdeeGGb+SsuueQSHnzwQWbOnNnw5QDQr18/3n33XRISEsjPz+e8885j8eLF/PjHP+a+++7jjDPOoKysjE8//ZSnn36an/3sZ1x77bWcf/751NXVhTXtuz3a8ox34MSLPTSNFc8HTtrvMxsh2bK7mnW7KpoI8QfrnbBEhRuW+N7UwSHDEqP6pzLUwhK9kydPha/mOu/9Xqf89fNOua7aKa90E/B4ypzy6vlOuWq3U173llOuKAzrlPPmzWPOnDkAzJkzh3nz5rVb55133uGdd95h8uTJTJkyhbVr1zZ4lYcccggzZsxote7ZZ58NwNSpUykoKADg448/5oILLgDguOOOY/fu3ZSVlYWsX1ZWRmlpKTNnzgRoqAfg9Xr5yU9+woQJEzj33HMb0nbOnDmTDRs2UFRUxLx58zjnnHOIiYnhyCOP5P/+7/+4++672bx5c6dklWvVM1bVr4GvRWSuqponHGECAWVDcSW7K5vmlnjyk018mF/C2AFpXH3cCDJstITRA9i9ezfvvfceK1euRETw+/2ISMPKHa2hqtx88838z//8T5PtBQUFJCcnt1m3Pv1mferN+vaa01p+FVVtdd99991HTk4OX3/9NYFAgISExoyFF1xwAXPnzuW5557jiSeeAOCHP/whRxxxBG+88QYnnXQSjz32GMcdd1yb9rdHW2GKF1T1+8BXzYa2iXNdmrdfZzYa8PoDrNtVQYWn8Ttv295q7l/ghCXOnpLLOW2EJQ7JSrIEP72dS95ofB8d27Qcl9S0nJDetJyc1bScmtPu6V566SUuvPBC/va3vzVsmzlzJh9/3DRlTWpqKhUVFQ3lk046idtvv53zzz+flJQUtm/fTmxsbPvX1wrHHHMMc+fO5fbbb2fRokX07duXtLS0kMf26dOH9PR0Pv74Y44++mjmzp3bsK+srIxBgwYRFRXF008/3STscPHFFzN9+nT69+/PuHHjAPjmm28YPnw411xzDd988w3Lly+PnBgD17p/T9uvMxhtUlPnZ+2u8oaUl4A7iWMT8bHR3Dx7DBNamcRxaHYKmcmhk8MbRiSZN28eN910U5Nt55xzDs8++yxHHHFEw7asrCyOOuooxo8fzymnnMI999zDmjVrOPLIIwFnUdNnnnmG6Oh9+1V3xx13cMkll5CXl0dSUhJPP/10m8c/+eSTXHrppSQlJXHSSY3R1iuvvJJzzjmHF198kWOPPbaJl56Tk8OYMWOaJC96/vnneeaZZ4iNjaV///786le/2if7gwknhWYyUKOqAXdx0tHAW6rqbbNiF3CgJwoqq/GyvrAxPuwLBHj8o00sWl/cZlgiNSGGEf0sLGEYXUF1dTUTJkxg6dKlDXmW94NWf8KGM7TtQyBBRHKBhcAlwFP7a1Fvp7iilrU7yxuE2OsP8OcF+SxaX8yZk3K5dfaYkEI8sI9N4jCMrmLBggWMHj2an/70p50hxG0STm4KUdVqEbkMeFBV/ygiX0XUqoOcrXuq2ba3cfxwrc/Pn95Zz/LtZVzyraGcOK5/izqxblgiw8IShtFlnHDCCWzZsqVLzhWWGIvIkcD5wGUdqGc0Q1XZWFxFcUVtw7aaOj/3vLOWtTsruPyY4Rw7ql+LehaWMIyDn3DCFNfhzLh7RVVXichw4P32KonIYBF5X0TWiMgqEbnW3Z4pIu+KSL77NyOozs0iskFE1onIQTWW2ecPsGZnRRMhrqr18Ye31rBuVwVXHTsipBBbWMIwegdhr/QhIqk4Q9oqwzx+ADBAVZe6dZcAZwIXA3tU9S4RuQnIUNUbRWQsMA+YDgwEFgAjVbXVqS0HSgeex+tn7a4KauoaL6Xc4+UPb65h694arj3uMA4fltmkjoUlDOOgZN878ERkghsjXgmsFpElIjKuvXqqulNVl7rvK3CywOUCZwD140+exhFo3O3PqWqtqm4CNuAI8wFNhcfLqh1lTYR4b3Udv3t9NdtLa7j+xFEthDg1IYYJg9JNiA2jFxFO7PdvwM9V9X0AEZkF/B34VrgnEZGhwGTgCyBHVXeCI9giUv/bPBf4PKjaNnfbAcueqjryCysIBP342F1Zy+/fWMPe6jpuPHl0i/XpBvZxckvYJA7D6F2EI8bJ9UIMoKqL3LHHYSEiKcDLwHWqWt6GyITa0SKGIiKXA5cDDBkyJFwzupySylo2FFUSHAUqLPfw+zdWU1Xr55bZYxiZk9qwLzpKODQ72TKtGUYvJZwOvG9E5HYRGeq+bgM2hdO4iMTiCPFcVf2Xu7nQjSfXx5WL3O3bgMFB1QfhJCtqgqo+qqrTVHVadnZ2OGZ0OWXVXjY2E+LtpTX85t+r8HgD3HZqUyFOiou2lJeG0csJR4wvBbKBf7mvvjgTP9pEHBf4cWCNqv4paNd84CL3/UXAa0Hb54hIvIgMAw4DvgznInoSlbU+1jULTWzeXcVvX19NQOH208YyPDulYV/flDjG56aTGGejJQyjN9NWoqAE4ApgBLAC+EUHp0AfBVwArBCRZe62W4C7gBfcSSRbgHMB3GFzLwCrcfIlX9XWSIqeSE2dn7U7y/EHKfHG4kr+8NYa4mOiuXX2GAb2cVLticCQzKSGsmEYvZtWh7aJyPOAF/gIOAUoUNXrus609ulJQ9tqfc7KHLVBCX/W7arg7rfXkhIfw22njqFfmpOWLy5GGNEvlfTEfc9WZRjGAUmrnWZtdeCNVdUJACLyOAdgyKCr8PkDrN1Z0USIV+0o457/rCMzOY5bZ49piAenJji5h+NjLCxhGEYjbYlxQ0hCVX021Co0gYCydlcF1UHjiNfuKuee/6wjOzWeW2c3rtrcPz2BoZZ72DCMELQlxhNFpNx9L0CiW65PLh86g3MvQlVZX9Q0KfzG4kr++HajR9wnKY7oKGFY32SyU220hGEYoWlr2SX7Hd0OG4ur2FvV2Ke5dU81d73lxIjrhTghNopR/VNJirPcSoZhtI4pxD6yZXd1k6Q/O8tquPPNNcRGC7ee6sSIM5JjGZGdQkx0OCMIDcPozZgY7wM7SmvYXtqYj7i4opY731iDqnLraePISUsgOzWeEf1S2mjFMAyjEXPZOkhxRS2bd1c3lPdW13Hnm6vxeP3cPHsMuX0SyUqJ49DssGeMG4ZhmBh3hL1VdWwsbswgWu7xcucbayir8XLjyaMZmpVMRnIsh/VLsREThmF0CBPjMKnweMkPyjdRVevjD2+uoajCww0njuKwnFTSEmMY2S/VhNgwjA5jYhwG1XU+1u2qaJjm7PH6+eN/1rJ1bw0//85Ixg5MJzUhhtH904iKMiE2DKPjmBi3g6qyvrASr7uKc50vwL3vrGNDUSXXHHcYkwZnkBQXzaj+qUSbEBuGsY+YGLdDUUVtwyodPn+A+xesZ/WOcq6YeSjTh2WSEBvFmAFpxNrwNcMw9oPeqSBPngpfzXXe+71O+evnnXJdtVNe+TL+gLKjcBdj3zmP9M1v89D7G/hqaym3Zi7k9ITlxMdGMSa1hrh/fhfyFzj1y7Y59Te6+fj3bHLKBR875ZJ8p7zlC6dcuNopb1/ilHcud8o7lzvl7UuccuFqp7zlC6dcku+UCz52ynvcFNMb33fKZduccv4Cp1xR6JTXveWUq3Y75dXznbKnzCmvfNkp17kjRr5+3in73cktX811yvUseQqePr2x/OXf4ZlzGsuf/xWendNY/uQBeP5HjeWP/gQvBmVk/eCP8PJPGsvv3QmvXtlYXnAHzL+msfyfW+GNXzSW37rJedXzxi+cY+qZf43TRj2vXumco56Xf+LYUM+Llzg21vP8j5xrqOfZOc411vPMOc49qOfp0517VE+Yzx7g/E+ePNX5H4HzP3vyVOd/CM7/9MlT7dmrp6ufvU6md4pxmOwsq8HrUwIK961M5ItNe7hwajbnpq4kJhrGDkgjwRL+GIbRCYS9OnRPJJIpNL3+AMu2luL1BXjy0wLeXV3I96YO4pwpg4iNFsYOTLMpzoZhdJR9Xx26t7KjtAavL8CzX27h3dWFnJY3gLMn5xIdJYweYEJsGEbnYmIcAo/Xz64yDy8v3c7ry3dywpgcfjh9CNFRwqj+qaTEmxAbhtG5mBiHYGNxJa8t28HLS7cxc2Q2lxw1tEGIbXUOwzAigYlxM3aW1fDSkm08++UWjhyexeXfHk50lDCiX0pDknjDMIzOxn5vB+Hx+nlnVSFPfVrAlCF9uPLYQ4mNiWJY3+SGZZMMwzAigYmxi6qyZPNe/rwwn5zUBK46dgSZyXEcmp1CQqwNXzMMI7KYGLtsL63h3v+so7zGy+/OHM/43HRy3NWcDcMwIo3FjHESAT3ywUa+2lrKj2YcwvGj+5kQG4bRpfR6MVZVXv96J/O+2Mr0oZn8cPpg+pkQG4bRxfR6MV66ZS93vrmG7NR4fnr8CIZn21JJhmF0Pb1ajPdU1vLLl5ZT6/Nz/UkjmTS4jy0eahhGt9BrlcfvD/CLF79mY3EVV80awVEj+pKaYBM6DMPoHnrtaIr7FuTz/rpizp06iO8fPtg67AzD6FZ6nGcsIieLyDoR2SAiN7Vfo+O8s2oXDy/awPShmdxw8igTYsMwup0eJcYiEg08DJwCjAXOE5GxnXmOzbur+PkLXzMoI5H750ykX6oJsWEY3U+PEmNgOrBBVb9R1TrgOeCMzmo8EFCunLsUgL/9aCoD+yR1VtOGYRj7RU+LGecCW4PK24AjOqvxqCjhtlPHUOP1M3Zgemc1axiGsd/0NDEOlQW/yVIkInI5cDnAkCFDOnyCIw/tu0+GGYZhRJKeFqbYBgwOKg8CdgQfoKqPquo0VZ2WnZ3dpcYZhmFEip4mxv8FDhORYSISB8wB5nezTYZhGBGnR4UpVNUnIlcD/wGigSdUdVU3m2UYhhFxepQYA6jqm8Cb3W2HYRhGV9LTwhSGYRi9ElHV9o/qoYhIMbC5g9X6AiURMKejmB1N6Sl2QM+xxexoSU+xZV/tKFHVk0PtOKDFeF8QkcWqOs3sMDtao6fYYna0pKfYEgk7LExhGIbRAzAxNgzD6AH0RjF+tLsNcDE7mtJT7ICeY4vZ0ZKeYkun29HrYsaGYRg9kd7oGRuGYfQ4Dmgxbi8RvTg84O5fLiJT2qsrIpki8q6I5Lt/MyJlh4gMFpH3RWSNiKwSkWuD6twhIttFZJn7mh3h+1EgIivccy3en/uxn/dkVNA1LxORchG5LoL3ZLSIfCYitSJyfTh1I/SMhLSjs5+RTrgnnfac7Mc96epn5Hz3GV0uIp+KyMT26u7T50ZVD8gXznTpjcBwIA74Ghjb7JjZwFs42eBmAF+0Vxf4I3CT+/4m4O4I2jEAmOK+TwXWB9lxB3B9V9wPd18B0DdEux26H51hS7N2dgGHRPCe9AMOB+4MbrsbnpHW7Oi0Z2R/benM52R/7ejiZ+RbQIb7/hQioCOqekB7xuEkoj8D+Ic6fA70EZEB7dQ9A3jaff80cGak7FDVnaq6FEBVK4A1ODmd94X9uR9t0dH70Zm2HA9sVNWOTuwJ2w5VLVLV/wLeDtTt9GekNTs6+RnZL1vaocvuSTO64hn5VFX3usXPcbJJtle3w5+bA1mMQyWib/6QtnZMW3VzVHUnOB8EnG/nSNnRgIgMBSYDXwRtvtr9afREGD9z9tcOBd4RkSXi5Iyup6P3ozNsqWcOMK/Zts6+J/tSNxLPSLt0wjPSGbZ01nPSKfeErn9GLsP5Rdde3Q5/bg5kMW43EX0bx4RTtyvscHaKpAAvA9eparm7+a/AocAkYCfw/yJsx1GqOgXnZ9hVInJMO+eLpC2Ik0L1dODFoP2RuCeRqNvpbXXSM9IZtnTWc9IZ96RLnxERORZHjG/saN1wOJDFuN1E9G0c01bdwvqfy+7fogjagYjE4nzI5qrqv+oPUNVCVfWragD4O85PoojZoar1f4uAV4LO19H7sd+2uJwCLFXVwvoNEbon+1I3Es9Iq3TiM7LftnTic7Jfdrh02TMiInnAY8AZqro7jLod/twcyGIcTiL6+cCF4jADKHN/MrRVdz5wkfv+IuC1SNkhIgI8DqxR1T8FV2gWPz0LWBlBO5JFJNU9bzJwYtD5Ono/9suWoP3n0eznZ4Tuyb7UjcQzEpJOfkb215bOfE46YyGJLnlGRGQI8C/gAlVdH2bdjn9u2uvh68kvnB759Tg9mre6264ArnDfC/Cwu38FMK2tuu72LGAhkO/+zYyUHcDROD9rlgPL3Ndsd98/3WOXu//YARG0YzhOT/DXwKr9vR+d8L9JAnYD6c3ajMQ96Y/j4ZQDpe77tG54RkLa0dnPyH7a0qnPyX7+b7ryGXkM2Bt0/xdHQkdsBp5hGEYP4EAOUxiGYRw0mBgbhmH0AEyMDcMwegAmxoZhGD0AE2PDMIwegImxcUAjIn5xsnOtFJEXRSSpg/UHishL7vtJEpTlS0ROlxBZvJrV/62InOC+v66j5zeMemxom3FAIyKVqprivp8LLNFmkyM60NbFOOOdr97H+gVu/Z6werFxgGGesXEw8REwQpxcsq+6yWI+d6eyIiIzpTHP7VcikioiQ12vOg74LfADd/8PRORiEXlIRNLFyeMb5baTJCJbRSRWRJ4Ske+JyDXAQOB9cfIPXyYi99UbJiI/EZF9+pIwegcmxsZBgYjE4OQqWAH8BvhKVfOAW4B/uIddD1ylqpOAbwM19fXVSYH4K+B5VZ2kqs8H7SvDmXU20930XeA/quoNOuYBnLwEx6rqsTjpFE9380oAXAI82akXbRxUmBgbBzqJIrIMWAxswcnjcDTOtFhU9T0gS0TSgU+AP7lebB9V9XXgPM8DP3Dfz3HLraKqVcB7wGkiMhqIVdUVHTif0cuI6W4DDGM/qXE93Qbc5DrNUVW9S0TewMkn8Lnb8eYJ8zzzgT+ISCYwFUdo2+MxHM98LeYVG+1gnrFxMPIhcD6AiMwCSlS1XEQOVdUVqno3jic9ulm9CpyljVqgqpXAl8CfgddV1R/isCb1VfULnBSLP6RlAnTDaIKJsXEwcgcwTUSWA3fRmMrwOrez7mucePFbzeq9D4yt78AL0e7zwI9oPUTxKPCWiLwftO0F4BNtXLbHMEJiQ9sMI4KIyOvAfaq6sLttMXo25hkbRgQQkT4ish4npm1CbLSLecaGYRg9APOMDcMwegAmxoZhGD0AE2PDMIwegImxYRhGD8DE2DAMowdgYmwYhtED+P9wktgeRyYJuwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 363.6x205.2 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[5.05, 2.85])\n",
    "plt.plot(\n",
    "    samples.groupby(\"positivity\")[\"n_positive_pools\"].median(),\n",
    "    label=\"Median\",\n",
    ")\n",
    "plt.fill_between(\n",
    "    samples[\"positivity\"].unique(),\n",
    "    samples.groupby(\"positivity\")[\"n_positive_pools\"].min(),\n",
    "    samples.groupby(\"positivity\")[\"n_positive_pools\"].max(),\n",
    "    alpha=0.25,\n",
    "    label=\"Range\",\n",
    ")\n",
    "plt.hlines(\n",
    "    n_positive_pools,\n",
    "    min(positivities),\n",
    "    max(positivities),\n",
    "    ls=\":\",\n",
    "    color=\"C1\",\n",
    "    label=\"After holidays\",\n",
    ")\n",
    "plt.title(f\"{n_persons} persons in {n_pools} pools\")\n",
    "plt.xlabel(\"Positivity\")\n",
    "plt.ylabel(\"Positive pools\")\n",
    "sns.despine()\n",
    "plt.legend(frameon=False, loc=\"center right\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"pools.png\", transparent=False, dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "extraordinary-folder",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}