notebooks/score_vs_N.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance comparison for various values of N\n",
    "\n",
    "This notebook reproduces Figs. 1, S2 and Tables 1, S1, S2, S3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.chdir(\"../\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import glob\n",
    "import json\n",
    "from src import utils\n",
    "from src.utils import get_fig_dim\n",
    "from IPython.display import clear_output\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sns.set_theme(context='paper', style='ticks', font_scale=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the values of sensitivity s<sub>e</sub> and specificity s<sub>p</sub>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "se = 0.99\n",
    "sp = 0.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "name=\"score_vs_N\"\n",
    "width_pt = 469\n",
    "palette = sns.color_palette('husl', 6)\n",
    "type_dict = {\n",
    "    'N' : int, \n",
    "    'false_negatives' : int,\n",
    "    'false_positives' : int,\n",
    "    'k' : float,\n",
    "    'lambda_1' : float,\n",
    "    'lambda_2' : float,\n",
    "    'method' : str,\n",
    "    'num_of_infected' : int,\n",
    "    'num_of_tests' : int,\n",
    "    'avg_group_size' : float,\n",
    "    'r' : float,\n",
    "    'score' : float,\n",
    "    'se' : float,\n",
    "    'sp' : float,\n",
    "    'd' : float,\n",
    "    'seed' : int\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading file 8/8\n"
     ]
    }
   ],
   "source": [
    "files = glob.glob(\"./outputs/{name}_*_se_{se}_sp_{sp}_*.json\".format(name=name, se=se, sp=sp))\n",
    "\n",
    "dicts = []\n",
    "for fl_id, fl in enumerate(files):\n",
    "    clear_output(wait=True)\n",
    "    print('Reading file ' + str(fl_id+1)+'/'+str(len(files)))\n",
    "    with open(fl,\"r\") as f:\n",
    "        js = json.load(f)\n",
    "        sizes = []\n",
    "        for group_id in js['groups']:\n",
    "            sizes.append(int(js['groups'][group_id]))\n",
    "\n",
    "        avg_group_size = np.mean(sizes)\n",
    "        \n",
    "        for seed in js['seeds']:\n",
    "            seed_js = js.copy()\n",
    "            seed_js.pop('seeds', None)\n",
    "            seed_js.pop('groups', None)    \n",
    "            seed_js = seed_js | js['seeds'][seed]\n",
    "            seed_js['seed']=seed\n",
    "            seed_js['avg_group_size'] = avg_group_size\n",
    "\n",
    "            dicts.append(seed_js)\n",
    "\n",
    "input_df = pd.DataFrame(dicts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_df = input_df.astype(type_dict) # Fix types\n",
    "total_df['tests_per_contact'] = total_df['num_of_tests']/total_df['N']\n",
    "total_df['fns_per_contact'] = total_df['false_negatives']/total_df['N']\n",
    "total_df['fps_per_contact'] = total_df['false_positives']/total_df['N']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {}\n",
    "temp = total_df[['method', 'N', 'tests_per_contact', 'fns_per_contact', 'fps_per_contact']]\n",
    "binomial = temp[temp['method']=='binomial'].groupby(['N'])\n",
    "negbin = temp[temp['method']=='negbin'].groupby(['N'])\n",
    "results['binomial'], results['negbin'] = {}, {}\n",
    "results['binomial']['avg'], results['binomial']['std'] = binomial.mean().reset_index(), binomial.std().reset_index()\n",
    "results['negbin']['avg'], results['negbin']['std'] = negbin.mean().reset_index(), negbin.std().reset_index()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the values of the table containing the average number of tests, false negatives and false positives per contact, set the first field to 'binomial' for classic Dorfman's method (Dorf-Cl) or 'negbin' for our method (Dorf-OD). Set the second field to 'avg' for the average values or 'std' for the standard deviation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>tests_per_contact</th>\n",
       "      <th>fns_per_contact</th>\n",
       "      <th>fps_per_contact</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20</td>\n",
       "      <td>0.281900</td>\n",
       "      <td>0.003560</td>\n",
       "      <td>0.001760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>50</td>\n",
       "      <td>0.241802</td>\n",
       "      <td>0.002800</td>\n",
       "      <td>0.001448</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>100</td>\n",
       "      <td>0.199127</td>\n",
       "      <td>0.001900</td>\n",
       "      <td>0.001176</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>200</td>\n",
       "      <td>0.159320</td>\n",
       "      <td>0.001236</td>\n",
       "      <td>0.000990</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     N  tests_per_contact  fns_per_contact  fps_per_contact\n",
       "0   20           0.281900         0.003560         0.001760\n",
       "1   50           0.241802         0.002800         0.001448\n",
       "2  100           0.199127         0.001900         0.001176\n",
       "3  200           0.159320         0.001236         0.000990"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results['negbin']['avg']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "size_df = total_df[['method', 'N', 'avg_group_size']]\n",
    "size_df = size_df.rename(columns={'avg_group_size' : 'Average pool size'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc_df = total_df.pivot(index=['lambda_1', 'lambda_2', 'se', 'sp', 'r', 'k', 'seed', 'N'], columns='method', values='num_of_tests').reset_index()\n",
    "perc_df['Perc'] = 100*(perc_df['binomial'] - perc_df['negbin'])/perc_df['binomial']\n",
    "\n",
    "df_list = []\n",
    "means = []\n",
    "for N in [20, 50, 100, 200]:\n",
    "    temp_df = perc_df[(perc_df['N']==N)].sort_values(by=['Perc'])\n",
    "    means.append(temp_df['Perc'].mean())\n",
    "    num_of_rows_to_skip = int(0.05 * temp_df.shape[0])\n",
    "    temp_df = temp_df.iloc[num_of_rows_to_skip:-num_of_rows_to_skip]\n",
    "    \n",
    "    df_list.append(temp_df)\n",
    "\n",
    "perc_df = pd.concat(df_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on the figure, comment/uncomment the lines that determine the subplot letter identifier and the file name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAACmCAYAAAAf6h86AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxSElEQVR4nO3deVhUZf/H8feAwCCyC7ghuIK7qZSGS5ptaO5aPpmZ9rPUsKwsNxTJXFofRXkyl1xCyy3NtUtyIRF3JZcHBfcFF3BAQBCU+f3Bz/lFIA7IzJkD39d1eV1w5sw5n8GZmy/nvs99a/R6vR4hhBBCCAtjpXQAIYQQQoiiSJEihBBCCIskRYoQQgghLJIUKUIIIYSwSFKkCCGEEMIiSZEihBBCCItUSekARYmPj2fBggVs2rQJgD59+uDg4MDt27fZtWsXbdu2JSIiQuGUQojyIj4+nkWLFvHbb78B+W2Ovb09N27cIDo6mpYtW7J8+XKFUwpR8VhkkeLv70///v0NRcqoUaOoVasWAJcuXSI8PFzJeEKIcsbf358+ffoYipS/tzmnTp3iq6++UjKeEBWWxXb3aDSaQtvWrVtH7dq1ee211xRIJIQoz6ysim4OGzduTLt27cycRggBFlyk/FNMTAzXrl0DoE2bNgqnEUJUBBs3bgSge/fuCicRomJSRZHy1VdfMWLECGQGfyGEuSQnJxuKlBo1aiicRoiKySLHpPzT2LFjefXVV/nvf/+rdBQhRAUQHh7On3/+iYeHh9JRhKjQVHElBaBr1654enoCoNfriY+PVziREKK8Cg4OZv369VSrVk3pKEJUaKq4kvJQt27dANi3bx95eXkKpxFClGeenp706NFD6RhCVGgWeyXl70XI/fv3AahSpQq3bt0iNDTUcHugEEKUhb+3OQ/Hvz38w0gIoQyLvJISHx/P2rVrDd+HhITQuHFjrl+/TnR0NNnZ2VSvXl3BhEKI8uSfbc4vv/xCv3798PX1VS6UEAKNXm6ZEUIIIYQFstjuHiGEEEJUbFKkCCGEEMIiWXyR4ufnp3QEIUQFIe2NEJbF4osUIYQQQlRMUqQIIYQQwiJJkSKEEEIIiyRFihDCJHQ6ndIRhBAqZ5GTuT1OXl4e165dIzc3V+koFsfGxoYaNWpgZSX1pzC9adOmkZmZWWh7Xl4eycnJLFq0SIFU5jdjcijpZizKHF1dGR8WarbzlReTwsYDMG3yDIWTCGOpski5du0aTk5OODk5KR3F4ty5c4dr167JsgHCLJo0aUKLFi2ws7NjyZIlvP7662i1WgA2b96scDrzSdfpmGxT02znC9NdNdu5ypOUtBSlI4gSUuWf27m5uVKgPIKTk5NcYRJm07t3b+rWrUvNmjXx8/OjXr161KxZEy8vL/bv3690PCGEyqnySooQwvK4u7szfPhwbG1tOXXqFAMGDFA6ktk4urqW6urGvXv3sLOzK9X5ysLnU8eTlmq+qwvOLu6ETFGuq6WSxgZnFf6BOyEsDIDpkycrnMT8pEhRQExMDHFxcVy9epUvvvjisfvHxsayYsUKwsPDzZBOiNLp3Lkzbdu25ezZs1StWpVq1aopHclsSjs+ZOLEiUa1AaaSlprCy/6HzXa+bfGtzXauojjZuTBtsnI/79JKSUtVOoJiykWRkpeVjf7+/TI7nqZSJazstY98/NChQwwbNox33nkHBwcHjh07RseOHenXr99jj52bm8vvv/9OWFgYMTExBR7LyMjghx9+wNPTEwcHB65cuYKPjw/dunUjIiLiiV+XEKb03XffcfPmTWbMmMGuXbs4ffo0nTp1UjqWEELFykWRor9/n5vfLSmz43mOGVLs423atDFc2razs+Pu3bv06dMHZ2dnXnjhhWKfm5KSQk5ODgCBgYEFHnv//ff54IMPeOqppwC4ceMGhw8fxtrauvQvRggzsbGxISgoCIDnnnuOoKCgMi1SDh8+zOjRo9FoNCxdupTdu3fj5uZGRkYGgwYN4vTp08TGxqLT6RgzZgwA6enpbN++nT59+pRZjrJUo0YNRc/v7OJeqqsbpe2mcnZxL/FzypLSP+/Scnd24c6dNKVjKMJkRUpqaiozZ87kxIkTjBw5kqCgIBYvXlygUSkvKleuzODBg1mxYgU6nQ4bGxuOHj3KsGHDyMrKYurUqTRp0gQnJye0Wi0XLlzgjz/+4Pnnnzcc4+TJk1y5csVQoAB4eXnRpUsXJV6SECVWs2ZN7t+/z40bN1iwYAGOjo5levwDBw6wZ88eNBoNhw4dIjU1laFDhzJv3jzi4uLYsmULAwcOJDw8HJ1Oh6urK7/99huvv/56meYoS9euXVP0/KUdH6J0N1VpKf3zLq3pkyczceJEpWMowmR399y+fZvp06ezaNEitmzZYmhUevXqRVpaGnFxcaY6tSJq1qzJ4cOHOXr0KL1796Z3795Mnz4df39/srKy+PTTTxk1ahRBQUH4+voWKFAALly4gIeHR6HjPrydUwhL16lTJw4cOMCkSZMAynQMVUpKClFRUTz//PPExMQQHR1N3bp1Aahfvz67d+/G1dWV7Oxs8vLysLe3Z/fu3bRv316uRAqhYia7kvKwAbl+/TpvvvlmkY1KixYtTHV6s0tKSgLAwcEBAH9/fxISEgBwdHTE1ta20HOmTZvGkSNHqFOnDgMHDiQlRe7hF+q1c+dOPvvsM5Mc293dnbVr15KQkEBwcDABAQE4OzsDYGdnR3JyMuPGjWPXrl307duXtLQ0srOzcXBwYMuWLdSuXZumTZsWOGZ4eDhz5841SV5hWpPCxpdqzpOcu7mM+Hh4iZ/n7uwuE8ApxKRjUi5fvszXX3+Nm5sbTk5OhRqVf1Jro3Hv3j0iIyMZO3asYQKrrKwsmjRpUuzzHv7FCfkzdGq1WmJjY2nXrh0Ad+/e5fz58489jhCWoEqVKqxfvx4rKyueeeYZvLy8yvwcDRo0oG/fvuzfv5+srCwAMjMzcXFxoXLlygQFBaHX64mMjGTQoEGEhITwySefEBISwpw5cwocKzg4mODg4ALb/Pz8SpUrL1uPPlv//xs0YO1shf6Bnrx0fYF9rRw0aGw05KXrqetRjwepeflPsdFg5aBBf09PXtY/nuOsAT3k3cnfrtFqsNJqSpW1LCg9tiMlLYXrnQ6Z74S725TZoSaEhZXqbp2cu3d59+OPSvw8d2cXVd+6bNIixdvbm6VLl9KjRw8CAgIKNSr/VJaNhinFxcVx+/Ztli5dilar5eTJk4wYMYKgoCDu3LnDnDlz0Gq1fPrpp5w/f56zZ88SHx+Pv78/e/fu5cyZM1y6dInatWsbjmllZcX333/PrFmz2LlzJz4+Pri7u/PSSy9x/PhxkpKSuHHjhkkafiHKwosvvohGoyEuLo733nuPNm3alFk/ul6vR6PJ/6VsY2PDiBEjiI6OJigoiMTERDp06GDYd+vWrbzyyitA/vpB9vb2pKenl0mOR7m7PZfMTf8/iaJGC56zHXiQoiclJKvAvi4j7bBrUYnUedl0P/8ayZ/lP27X2hqX4Vqy9t4n/eecAs/xjKiMPgfDvg7dbajyauGrs+ai1rEdliAlLZVjgea7FbtljPluMTcFk9/dY2VlRYsWLejevfsjG5UnpalU6bF35JT0eMVp0aIFx44dK/KxUaNGFdq2d+9ew9cDBgx45CRXNWrUYPbs2YW2N2vWjKioqGIzCaG0999/n1u3buHr68uECRMICAgos2Nv27aNyMhIunTpQtu2bWncuDH79+9nzZo1ODk5Gc517tw5XF1dcXfPv4vk9ddfZ+3atbRubdpfCpVfsMG+w9/ajf+7yGHtrqHqLPsC+1o55D/oMkrLzC9m8tm4/C4yjU3+dvtnK2H31D/G0ViBxg7DsTQKXkURwpyMLlIejpbPyMigSpUqj91/yZIlnD17llatWvHaa6/RvHnzIhuVslDcnCZCCPOoVq0aM2bMMMmSFa+88orh6shDI0eOLLRf3bp1DWPfANq3b1/mWYpipdVAEYWDxlqDtUvRBYWVowZnbyesXQrev6Cx02BtV8RzNDzyWOamdHePu7N7qbpgSnvrtLtz2d067e7sUqqrG6XP7lLi51gSo4qUiRMnYmtry5QpU4iPj+fSpUuPnXdgyJAhhbYV1agIIcqH8ePHk5KSQkZGBnq9nj/++IPBgwcrHcuiKd1tUtqxNGmX7yg6lqa0g1gt4dbp0o4PsYTsSjCqSGnQoAH169cHoFWrVowbN85iJ0cSQigjJCQEnU6Hra0tWq2W+/fvS5Fi4Uo7lmaYbrQqx9II9TGqSHF0dCQpKYnDhw/z448/4uvra+JYQgi1eeaZZ+jZsyfr16+nd+/eTJs2TelIFk/pbpPSjqX5cdGPhqvlSoylqWh3U4Hy7xWlGFWkdO/endWrV7Nx40Zat25t0TM4CiGUcerUKa5du0bjxo3p3bs3eXl5SkeyeEp395R2LM25W2cVHUuj5rup1NrFphSjipTp06fz6aef4uDgwO3bt5k/fz4ffvihiaMJIdRkwoQJhq8bNmyIp6engmlEeabmu6mki61kjCpS0tLSePPNNxk7dizt2rWT7p4nFBMTQ1xcHFevXq2QA6FE+fThhx8ycuRIdu3axcaNG+nSpYthoT9RNLVewlc6t5rvplJrF5tSjCpSnnvuOTp06MCUKVPYunUrjRo1MnWuEsm4d4ecBzmP39FItta2VLF79G2Uhw4dYtiwYbzzzjs4ODhw7NgxOnbsSL9+/R577NzcXH7//XfCwsKIiYkp8FhqaipLly6lZs2aZGVlce/ePYYOHUp6ejrTp0/nwoULBAUFcf36dW7evMmnn34qk7sJi9G5c2ecnJxYtWoVGzZskLl9jKB0d09pSe7SU2sXm1KMKlIOHjxIp06dmDt3LuvXr+err75i4MCBps5mtJwHOYxaU/QEaaUxr9+qYh9v06YN7u7uDB8+HDs7O+7evUufPn1wdnbmhRdeKPa5KSkp5OTkF1SBgYGG7Xl5ebz77rtMnjzZMA3+woULiYiI4P333+eZZ57Bw8ODt956C4DVq1fz7rvv8uuvvxpm4hRCSXfv3mXSpElMmjSJ8+fPs3z5cnr27Kl0LCGEihlVpPy9S6JXr15lOltseVC5cmUGDx7MihUr0Ol02NjYcPToUYYNG0ZWVhZTp06lSZMmODk5odVquXDhAn/88UeBlZCPHDlCSkpKgXV6Xn75ZXr06FHkLLb9+vXj3//+N0eOHDH5bJpCGGPgwIEF/nhZs2aNgmmEKSnd3VNa6VlpTAobr7rFAieEhZGelfX4HcuhRxYpH330Ed988w0ajYawsDDDujsAly5dIjIy0iwB1aJmzZocPnzYMOumr68v06dPZ/78+WRlZfHpp59ibW1NUlIS586dK1CgANy8edMwlfdDnp6eZGZmotPpCp1Po9FQo0YNbty4YdLXJYQwHbX+sreEbpPSyNXklGr1ZKWlpKUaxq5UNI8sUkJCQgzdCC1btqRly5ZYW+cP0Dlx4oR50qlIUlISAA4ODgD4+/uTkJAA5M8zY2tbeHT1tGnTOHLkCHXq1GHIkCHcunWrwOM6nQ5HR0fc3NweeU5vb++yfBlCCDNS6y97IczlkUWKq6ur4etmzZphb2+PjY0N69evp1OnTmYJpxb37t0jMjKSsWPHsnnzZgCysrIKdN0UZdKkSYav9Xo9Xl5exMXF0aJFCwC2b9/OsGHDinzuhg0b8PHxoVmzZmX0KoR4Mnv27KF69erY2tqybNkyXn75ZemKLKfUegUoK1OdXSbuzi7cuZOmdAxFWD1+F5g9ezY2NjaMGTOGO3fuGH4RV1RxcXHcvn2bpUuXsmzZMiZPnsyIESN488036dChA3PmzGHNmjV8+umnnD9/nrNnzxIfHw/kr4h85swZLl26VOCYGo2GefPmsW3bNlatWsXy5cvJycnh3XffJS0tjSNHjnD06FFWrlzJ7NmzOXDgAOHh4Uq8fCGKtGPHDnx9ffnggw/w9/cnLi5O6UgWbcbkUO7cvMWMyaFKRykxtV4Bsnewx97B/vE7WpjpkyfjbFcxF9I1auBsx44d2bt3L1ZWVowePZr58+ebOleJ2FrbPvaOnJIerzgtWrTg2LFjRT5W1CDXvXv3Gr4eMGAAAwYUfSeSm5sbn332WaHtzs7OMsW4sHj16tXjxx9/pE6dOvTt25eZM2cqHcmipet0fOVYnzDdVaWjCGGxjCpSateuzf79+5k1axYnTpwgNTXVxLFKprg5TYQQ5hEYGMixY8f4/PPPOXHihGq7BET5pdbuHlBvF9uTMqpIadOmDW3atAGgatWqNG3a1KShhBDqc+rUKXr16gVA06ZN+fPPP5UNZOEcXV2ZeP0cVatVUzpKiXw+dTy3k6/z+dTxhExR1628auzqeUitXWxPyqgiRQghHiUpKYmvvvqK//73v/z888/o9fkLm1WtWtXk5168eDFubm5kZGQQEBBAbGwsOp3OMB1/eno627dvp0+fPibPUlLjw0KZOHEi48NCFU5SMmmpKfRt9Rfb4m2UjiIqgFIVKdnZ2Wi1FXMQjxCioOrVqzNz5kzOnTuHv7+/2c576NAhUlNTGTp0KPPmzWPatGl8/vnnhIeHo9PpcHV15bfffpNV24WBmrt7KqpHFik//PCDYfr2fzpz5gxz5swxWajHsbGx4c6dOzg5yViUf7pz5w42NvIXjjAvW1tbVq5cyRtvvMHJkyf5z3/+Q79+/Rg+fLjJzhkdHU3dunUBqF+/PhEREWRnZ5OXl4e9vT27d++mffv2hvmdLJEaxxk4u7izIa4pnl7uj9/Zwqitu2dCWFj+RG5Azt27vPvxR0D+LcnTJ09WMJn5PLJIsbKyonXr1kV+wCtXrmzSUI9To0YNrl27RkqK+mYONDUbGxtVNnxC/Zo3b46Xlxfvvfceq1atIjY21qTn0+l0ODs7A2BnZ0fv3r05d+4cffv2JS0tjezsbBwcHNiyZQu1a9cuNJYuPDycuXPnmjTj46hxnEHIlBlMnDiRkCmygruppaSlciyw8FxDLWMOK5BGGY8sUt55550C3588eZKEhAQaN27M0KFDTR6sOFZWVtSqVUvRDEKIgtLS0ggODmbUqFGkpaWxbNkyXn31VZOdz83NzbBcR2ZmJm5ubgQFBaHX64mMjGTQoEGEhITwySefEBISUujqb3BwMMHBwQW2+fn5mSyvEKLkjBqTEhERweHDh/H39+f48ePUr1/folZBFkIob+jQobz55pvY2NiQnZ3N4sWLTXq+jh07Eh0dTVBQEImJiYaFT7du3corr7wC5F9tsbe3Jz093aRZSkutVz3VmlttY1LcnV1oGXOYO//3/nVydDRsryiMKlIcHR1ZtGiR4fvVq1ebLJAQQp2+++47bt68yYwZM9i3bx8ajcakS2i0bt2a/fv3s2bNGpycnAgICODcuXO4uroaFut8/fXXWbt2rcVOz6/G7h5Qb261jUl5OO7k4ViU+d98q2QcRRhVpGRlZZGUlISDgwMnTpxg37599O/f39TZhBAqYmNjQ1BQEADPPfccQUFBJl/na+TIkQW+r1u3rmEwLUD79u1Nen4hhGkZVaT06NGDr776ivj4eOrXr88nn3xS7P4ZGRlMmDCBkydP0qFDB0JDQ1m9ejVWVlbodDqGDh2KlZVRywYJIVSiZs2a3L9/nxs3brBgwQIc/+/StCh/1NjdMylsPOl30g1fT5usnonosjIzlY6gGKMqhWrVqvH222/zP//zP4waNQpvb+9i94+Li2PmzJls2rSJ2NhY/vrrLw4ePEjfvn2pWrUq27ZtK5PwQgjL0alTJw4cOGBY3VvpO2csnSwwaF4paSlkvnqazFdPk5KmrjtD7R0csHdwUDqGIoweOHvkyBH8/PyMGjgbGBho+Lphw4ZER0fj4+MD5M9nsHz5csNlYSFE+XD8+HHDApl6vZ41a9ZIt3AxZIFBIR7P6IGzCxcuNHxv7MDZjIwMqlevTqVKlahSpQqQP59BcnJykfsrPW9Bxr075DwoegK7J2VrbSsLIYpSM9V7syzel0lJScydO5cjR44YBtg/nBpfipTySY3dPe7O7qRvapj/dS11TURXkbt7TDpwdsOGDYwePZotW7aQlpYG5M9n4OLiUuT+Ss9bkPMgh1FrBpjk2PP6rTLJcUXFYKr3Zlm8L6tXr87YsWM5e/asxd5FY4nUusAgqLO7Z9rkGYz4eLjhazWpqF09YOSYlB49evD1118zcOBAfvnlFz788MPHPicqKoquXbtSpUoVAgMDSUxMBCgwn4EQonxwcXGRAqWExoeFovVwV90Cg0KYk1FXUqpVq0b//v25ceMGTZo0eezA2cjISBYtWoSLiwu5ubm89dZbNGvWjNWrV5OcnGzS9TyEEOYzb948WrduTdu2bZWOIsxIjd09oL7J3B7Kyszk/v0HTAgLqzBr9jxkVJEyevRokpOTeeqppzh8+DD169dn8ODBj9z/jTfe4I033iizkEIIy3Tnzh1DgXLgwAGefvppw2O3bt3Cw8NDqWiqoNZf9mrs7gH1Teb2kL2DA8cCW+NYgdbsecioIiUzM5MVK1YYvt+yZYvJAgkh1CMgIIA5c+ZgZWXFqVOnOHDggOGx48ePM3/+fAXTWT61/rIXwlyMKlJatWrF9u3bcXFx4cGDBxw6dAgPDw+OHz+u+GKDQgjldO3alRYtWnDp0iXu3r1b4ErKwzt8hLAUau3ucXd2odGOvbh7eikdxeyMKlIuX75MXl4eGo0GyF99dP/+/Zw9e9ak4YQQls/DwwMPDw8aNWqEra0tt27dwsvLSwbSGkGN3T2fTx3P7eTrfD51PCFT1HaXjDq7e6ZPnszEiRP5ooKNRwEji5TJkydTuXLlQtszMjLKPJAQQp0OHjzI559/jpubG5UrV2bYsGFyJ99jqLG7Jy01hb6t/mJbvI3SUUQFYFSRUlSBAhgmaBNCiIMHD7J582bs7OwA+PnnnxVOJERBau3umRAWRnqWOrM/KaOKFFE+WPKspWWhvL8+S1erVi1DgZKZmcnJkycVTiRMITMzi1X762NdSX2/NNXa3ZOSlgoapVMow6giZcWKFVy8eJHx48ezc+dOqlatSrNmzUydTZQxS561tCyU99dn6Xx8fHjjjTdIT08nIyOD0NBQpSNZPDWOSXFwsOdl/1Nsi5cxR8L0jCpSzp07Z5gLoXPnzvTp04d169aZNJgQQl3atWtH27ZtuX37Nm5uboaB9uLR1DgmRQhzMqpIadGiBVqtlry8PJYtW8a9e/dMnUsIoUIajQZ397JfvG3atGls3boVf39/Fi1aREpKCj/99BNVq1alUaNGtGrVip9++gmNRkOtWrXo1KkTAH/++Sc1atSgXr16ZZ6ponJ2cWdDXFM8vdS1SB9AJY0Nzk7q67p1d3bhzp00pWMowqi1e5o3b87WrVvp27cv+/fvZ/bs2abOJYQQANy4cYNGjRoRExNjWGX522+/pUePHrzxxhv88MMP6PV69uzZQ//+/Vm5ciWQP+Ntenq6xRYoMyaHcufmLWZMDlU6SomETJmBk2st1d1+DOBk56K6xQUh/xZkZzut0jEUYVSRkpyczBdffMGvv/7Kf/7zHzQaDZs3b+bBgwemzieEUIljx46xalX++J3o6GjDyudPKjY2loiICN59911u374NQExMDL6+voZ9rly5glarJSsry7DK+u+//05QUFCZZDCFdJ2Orxzrk67TKR1FCItlVJGybNkygoKCOHXqFADTp08nMDDQ0CAJIcTy5cvJy8sDIDAwkM8++6zEx/j+++8ZN25cgX9VqlQhKiqKdu3aMXPmTADu379vGPNia2tLcnIyQ4cOJSYmhjFjxrBt2zZefPFFDhw4wObNm8nOzi50rvDwcPz8/Ar8E8ZR44BfgC5duigdodTUnP1JGDUmJSAggFmzZrFjxw4aN27M7du3cXFxITEx0dT5hBAq0blzZ8NVjJMnT3LixIkSH+O999575GNDhgxhzJgxQMG5mzIzM3FxcaFOnTo0b96cCxcu4ODggKurK7NmzWLo0KEsW7as0OrrwcHBBAcHF9hmzkLF0dWVidfPUbVaNbOds6yodcDv888/r3SEUlNz9idh1JWUzMxMduzYwd69e7lw4YLhMu7Vq1dNGk4IoR5eXl6sWrWKIUOG8O677zJx4sQyOe7DNYBycnIMUx88/fTTXLx40bC9Tp06QP4VltjYWDp06GC4emJvb8/169fLJEtZGh8WitbDnfFhoQonEcJyGXUl5aWXXmLFihU8++yz5OXlMXr0aBYtWmRoGIQQolatWsyePZvbt2/j6upKTk4OFy9exMfH54mO+8EHH+Ds7EyjRo14/fXXgfyrIA/v7hk9erRh302bNvHqq68C4OjoSIMGDfj999/p2bPnE2UQBam1u0eoj1FFiq+vLxMmTDB8f+7cOXr16mWqTEIIFZo/fz6QX0BYWVkREhJC7969SUxMfKJL1XPmzCm0zcvLi48//rjQ9n+2S8V1H1kCtf6yV2t3j1Afo4qUqVOnEhcXh62tLdbW1mi1Wrp27WrqbEIIFdFqtbRr146dO3fSr18/EhMTefbZZ/nyyy8rbH/648gveyGKZ9SYlKeeeopffvmFt956i8jISJo2bWrqXEIIlWnYsCFNmjQxjAXJzc0FZOxaeVRR7zQR5mdUkXLgwAEiIyPJzc3lvffeY+vWrabOJYRQGScnJ9555x3OnTvHrl27aNmyJT/99JPhtmRRmFq7e+TKmDAXo7p7unXrRsuWLbG3t6dx48a4ubmZOpcQQmW6du1aoBvY29ubWrVqyfi1Ykh3jxDFM6pIWbhwoWE66vr16xtmfRRCiIciIyOJioriwYMH6PV67t69y9q1a7Gzs1M6mhBCpYwqUpycnBgxYoRhoqbTp0/LKshCiAL0ej1ff/01CQkJPP3006xYsULpSEIIlTOqSOnUqROtW7fGyip/CEtcXJxRBz948CBz585l6dKl5OXlERERgbe3N3l5efTu3bv0qYUQFuf06dN4e3tz4cIF4uPjWblyJYMGDVI6lkVT65gUIczFqIGzPj4+xMbGUrNmTc6ePUtgYKBRBw8ICODevXsAbNy4EQ8PD3r27MnRo0dJSkoqfWohhMUZNWoUVatW5V//+he2trZMnjxZ6UgWT8akCFE8o4qUv4/QL+nCYTY2NkD+qqgPl0z39fVl7969Jc0qhLBg+/bto0mTJtjY2PCvf/2LXbt2KR1JCKFyRnX3lMXCYTqdDicnJwDs7OxITk4utE94eDhz584t8bGFsBR5Wdno798v+wNb6cv+mGUkKSmJGTNmkJCQwLp169Dr9VhZWdG4cWOlo1k86e4RonhGFSleXl4sX76chQsXcvr0aaZMmVLiE7m6uhomeXq4auk/Kb0qqRBPSn//Pje/W1Lmx7X7+LUyP2ZZqV69OrNnz+by5cvUrl3bsL2oP0REQdLdI0TxjCpSateuXWDhsIcDaEuiU6dOnD59mubNm3P+/HlGjhxZ4mMIISyTRqNh//79TJ061TDT7I0bN/j9998VTiaEUDOjipR58+bRqlUrnnnmmRIVKKdPn+bSpUucOXOGbt26MXv2bNasWUOrVq3w9vYudWghhOXZtGkTkydPRqvVAvldw6J40t0jRPGMXmBQo9Fw8OBBli9fTtOmTQkKCnrs8/z8/Ni9e7fh+48++qj0SYUQFq1///54eXlRpUoVAFJSUhROZPmku0eI4hlVpERFRXH16lU2bdpEgwYNpPoXQhSyePFifv75Z6ysrNDr9Vy/fp3t27crHUsIoWJGFSnffPMNb7/9NkuWLKFKlSrcN8XdC0IIVYuIiMDW1hZXV1cyMzPJyspSOpLFkz/4hCieUUXKkiVLqFatGpA/9XV4eDhjxowxaTBRPpnsFl2w6Nt0K4Lw8HBsbW2ZMmUK8fHxXLp0iT59+pToGNnZ2SxZsgQrKyuGDx8OwPnz59myZQv29vZ07tyZOnXqFLnthx9+AKB9+/aG2583btxIYGCgxS6KKt09QhTPqCKlWrVqJCcns2bNGtasWUNWVpYUKaJUTHWLLlj2bboVQYMGDahfvz4ArVq1Yty4cSUuUrRaLU2bNuXo0aOGbV988QWzZ8/GxsaGjz76iLlz5xbaNmPGDM6cOcPIkSNZsWIFjRs3JiEhAXd3d4stUIQQj/fYImXv3r2sXLmS+Ph4GjVqRGRkJHq9/MUqhCjI0dGRpKQkDh8+zI8//oivr2+pjvNwlmrIv7Jy+fJlHBwcALhy5QoZGRmFtmm1WjQaDdnZ2Tg7O5OTk8OhQ4cYOHDgE78uIYRyii1Spk2bxrZt25gyZQqzZ89mw4YNeHl5mSubEEJFunfvzurVq9m4cSNt2rRhwIABxe7//fffc+HChQLbunbtiqOjo+H7tLQ0QzECYG1tTWZmZqFtOp2Ovn37cv78eYYMGcKGDRvo2bMn0dHRZGRk8PLLLxeaPsESZriWMSlCFK/YImXSpEkMGTKEHTt2sHLlSsMkTTk5Odja2poloBBCHc6dO0fTpk0ZNGgQsbGxj93/vffeK3L7/v37DV+7urqSk5Nj+D47O5vKlSsX2ubk5ETbtm0BOHr0KP7+/ty5c4c9e/bQpk0btm7dSrdu3QqcxxJmuJYxKUIU77Ezs9WqVYvBgwfTv39/PD09Wbx4MdOnTzdHtgorLyubB+kZZf4P6aYTJhQREWHoqmndujWTJk164mPa2tpSo0YNsrKyuHfvHtWqVcPR0bHQtocTyGVkZHDx4kWaNWtGamoqWq0We3t7rl+//sRZhBDmZ9TAWchvLB5O4LZq1SqTBRIVc/0XoX4vvPACTZo0AfIXFD127FiJj3H//n2OHj1KQkICaWlpODs788knn7Bw4UJsbW0ZP348QJHbIP9unv79+wP5A3ltbGw4cODAY7uelCLdPUIUz+gi5e8s9QMvhFCOp6cn48aNIzMzk1OnTvHJJ5+U+BiVKlUq1A3UsGFDGjZs+NhtQIGBstbW1oW6cyyNdPcIUbxSFSlKM9lcGzLPhhCllpSUxPTp09HpdKVeiFQIIf5OlUWKdIcIYXmio6Pp3bs37u7uAFy9epWaNWsqnMqySXePEMVTZZEihLA8d+/epWfPnjg7O6PX67l06VKBBUZFYdLdI0TxpEgRQpSJt956i7Zt2xq6eeLi4hROJIRQO+k0FkKUibt377JgwQIAdu3aZej2EY8m3T1CFE+KFCFEmfjjjz+oXr06AM899xzjxo1TOJHlk+4eIYon3T1CiDLRoUMHnJycANi+fbv8AhZCPDEpUoQQZcLHx4fFixcza9YsNBoNs2bNUjqSRZsxOZScjEylYwhh0aRIEUI8scTERDw9Pfnmm2+UjqIa6TodNo/fTYgKTcakCCGeSHh4OP3796dbt25GLSwohBDGkiJFCPFEzp49y/79+4mKiuLo0aNKxxFClCNSpAghnkjDhg3RaDRotVrq1atHTk4OOTk5bNq0SeloFs3R1RVHV1elYwhh0WRMihDiiSxYsIC1a9ei1+evffXll1+i1+tJTU2le/fuCqezXOPDQpWOIITFM3uRsnjxYtzc3MjIyGDQoEHmPr0QoowtXLiQ1q1bF9r+119/KZBGCFGemLW759ChQ6SmptKrVy/S0tJk2mwhyoGiChSA5s2bmzmJEKK8MWuREh0dTd26dQGoX7++LD4mhBBCiEfS6B92JJtBSEgIXbp0oXPnzuzatYsdO3YQFhZmeDw8PJy5c+eaK44QohinT59WOoLZ+fn5KR1BiArpUe2NWcekuLm5kZWVBUBmZiYuLi4FHg8ODiY4ONickQrw8/Mr1w2zvD71Ks+vzZKY+2es1v9XyW1+as7+JMza3dOxY0fDDzkxMZEOHTqY8/RCCCGEUBGzFimtW7fGzs6ONWvW4OTkREBAgDlPL4QQQggVMfstyCNHjjT3KYUQQgihQjLj7N+8//77SkcwKXl96lWeX1tFptb/V8ltfmrO/iTMenePEEIIIYSx5EqKEEIIISySFClCCCGEsEgVeoHBjIwMJkyYwMmTJ+nQoQOhoaGsXr0aKysrdDodQ4cOxcpK3XXc4cOHGT16NBqNhqVLl7J7925Vr5108OBB5s6dy9KlS8nLyyMiIgJvb2/y8vLo3bt3kdvUwtj3Y3l7j1YUam9v1NSWqLGdkM9/0cr3q3uMuLg4Zs6cyaZNm4iNjeWvv/7i4MGD9O3bl6pVq7Jt2zalIz6xAwcOsGfPHvbs2YNOp1P92kkBAQHcu3cPgI0bN+Lh4UHPnj05evQoSUlJRW5TC2Pej1euXCl379GKQu3tjZraEjW2E/L5L1qFLlICAwOpXLky9vb2NGzYkOjoaHx8fIDysbZQSkoKUVFRPP/888TExJSbtZNsbGyA/LWg6tWrB4Cvry979+4tcptaGPN+jImJKVfv0YpEze2NGtsStbUT8vkvWoXu7nkoIyOD6tWrU6lSJapUqQKAnZ0dycnJCid7Mu7u7qxdu5aEhASCg4MJCAjA2dkZKB+vT6fT4eTkBPz/6ylqm9oU937U6XTl6v+wIlJje6PmtkRt7YR8/guq0FdSHtqwYQOjR49+7NpCatWgQQP69u1LUlJSuXp9rq6uZGdnA///eorapjbFvR/L63u0IlFze6PGtkRt7YR8/guq8EVKVFQUXbt2pUqVKgQGBpKYmAiUj7WF/j4Fjo2NDSNGjChXayd16tTJ8HrOnz/Ps88+W+Q2NXnc+7F9+/bl6j1a0ai1vVFzW6KmdkI+/4VV6MncIiMjWbRoES4uLuTm5vLWW2+RnZ1tuIw2fPhwrK2tlY5Zalu3biUyMpIuXbrQtm1bGjduTEREBJ6enqSnp/P2228rHbHETp8+zfDhw1mwYAH16tVj9uzZ1K5dG71eT//+/Xnw4EGhbWph7Pvxp59+Kjfv0YpEze2N2toSNbYT8vkvWoUuUoQQQghhuSp8d48QQgghLJMUKUIIIYSwSFKkCCGEEMIiSZEihBBCCIskRYoQQgghLJLMOCsMVq9eTXx8PCEhIWV63KioKA4cOICTkxPvv/9+mR5bCKFO0t4IY8iVFBXLysril19+4YUXXuD27dsA3L17l2+++Yaff/6ZBw8elOh4bdq0ITMzs0wz3rx5kxMnTjBu3Dhat279xMfLy8vj+vXrRu9/7dq1Jz6nEELaG2NIe1P2pEhRMXt7e1577TWqVq3Khx9+yIMHD6hcuTLt27enffv2JZ7k5+GCXGXp0qVLWFtbY2VlRbt27Z74eMuWLePixYtG7avT6Zg7d+4Tn1MIIe3N40h7YxrS3VMO9OvXj127dvHll18yfvx4IH9Z9fHjx7Ns2TLmz59Pbm4uL774IqGhofTq1YsNGzYwYMAALl68SGxsLP/+97+B/NVOJ0yYwLFjx4iIiMDHx4dff/2VnJwc9u7dyxdffEFoaCje3t7s3LmTDRs2GHLExsZy4cIFLl++zNNPP01AQAB//PEH58+fZ9++fbRt29aw7/r169FoNGzatInvvvuO9PR0fvvtN+zt7UlPT2fIkCFMmzYNX19fTp8+jZ2dHZMmTeLgwYMkJyfj7e3NvHnzCAgI4NixY4SGhnL16lWio6M5c+YMTZs2xcfHhyNHjrBz504cHR25ePEiu3fv5oMPPjCsgCqEKBlpb6S9MSe5klIOaDQaZs2axd69e9m4cSMATz/9tOGxp556CgA/Pz/u3btH7969mTp1KlFRUXzwwQf06NGDgwcPAqDVapk+fTqvvfYaCxYsYOfOnVy/fp0aNWpQs2ZNbt68ibe3N40aNWLdunWGDHl5eYSHhzNw4EBGjx7N559/jo2NDc899xxNmjQp0GAkJCRw8+ZNevbsSe/evcnKyuLrr7+me/fuDB48mCNHjnDu3Dm8vLyoV68e3377LXv27MHBwQF/f386dOiAi4sL9erVIygoiLi4OADmzp1L//79CQ4OxtbWljZt2uDh4UHnzp3Zt28fdnZ2hIaGGlYWFUKUnLQ30t6YkxQp5UTlypWJiIjgm2++IT4+/pH7OTg4YGtri729veHDo9Vqyc3NNTwO+Y1OSkoKiYmJ1KhRg44dO/LZZ59Rp04drK2tcXFxKXB59/bt26SnpxuO5+rqyq1bt4rMkJCQgFarBSAoKAgPDw9OnTplOLefnx+JiYlYW1vj6OgIFL40XLlyZdzd3YmKisLW1hbIXyisUqVKuLm58eqrrxbY/7XXXmPdunWMGTPGJJeZhahIpL2R9sZcpEgpR7y9vZkxYwbfffcdAPfv3wcgNTWVvLy8Eh3r+vXrtGzZEh8fH1asWEF2djZnz57l8uXLRe7v5uZGRkYGGRkZQP6H2svLq8h9fXx82Lx5M7m5uVy7do3Tp0/ToEEDTpw4AUB2djaNGjUq8rkajQa9Xk9CQgJbt24lKCgIa2tr9Ho9Li4u7NixA4Bdu3YB+au3Phz8tnjxYlq1alXgkrEQonSkvZH2xhykSFG5v/76i127dnHjxg0A2rVrx4cffghAvXr1GD16NAkJCVy5coX4+HguXbrEhQsXOHbsGJcuXeL69eucPHmSU6dO4eLiwv379/nll184deoUQ4cOpWvXrtSuXZuXXnqJ7du34+npyZkzZ9i7d2+B5dutrKwIDQ1lxowZrFixgmHDhgGwf/9+/vrrL0M+gCZNmtCyZUu6d+/Ozz//jJ+fH2PHjmXdunWsWrWKevXqUadOHeLj44mLi+PixYukp6dz/Phx/Pz8WLJkCXfv3iUxMZFvv/0WrVbL9u3b+fjjj/n6668ZPHgw1apVM5xv9erVbN68me+//x4rKys6duxohv8ZIcofaW+kvTE3WQVZCCGEEBZJrqQIIYQQwiJJkSKEEEIIiyRFihBCCCEskhQpQgghhLBIUqQIIYQQwiJJkSKEEEIIiyRFihBCCCEskhQpQgghhLBI/wtN+y4DIzSYJQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 560.697x173.265 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "utils.latexify(font_serif='Times New Roman', font_size=10, usetex=False)\n",
    "\n",
    "fig_height, fig_aspect = get_fig_dim(width_pt, fraction=0.6)\n",
    "fig, [ax1, ax2] = plt.subplots(nrows=1, ncols=2, figsize=(2*fig_height*fig_aspect,fig_height))\n",
    "\n",
    "# Subplot 1\n",
    "sns.barplot(data=size_df, x=\"N\", y=\"Average pool size\", hue=\"method\", ci=None, \n",
    "            palette=[palette[0]]+[palette[2]], hue_order=['binomial', 'negbin'], ax=ax1)\n",
    "sns.despine(ax=ax1)\n",
    "\n",
    "# ax1.text(x=0.01, y=1, s=r'A', size=15, transform = ax1.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "# ax1.text(x=0.01, y=1, s=r'C', size=15, transform = ax1.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "ax1.text(x=0.01, y=1, s=r'E', size=15, transform = ax1.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "\n",
    "ax1.legend(labels=['Dorf-Cl', 'Dorf-OD'], loc='upper left')\n",
    "ax1.set_xlabel(\"Number of contacts\")\n",
    "ax1.set_ylim([0,35])\n",
    "\n",
    "\n",
    "# Subplot 2\n",
    "labels = perc_df['N'].unique()\n",
    "data_sets = [perc_df[perc_df['N']==N]['Perc'].to_numpy() for N in labels]\n",
    "\n",
    "hist_range = (np.min(data_sets), np.max(data_sets))\n",
    "number_of_bins = 30\n",
    "binned_data_sets=[]\n",
    "for d in data_sets:\n",
    "    hist = np.histogram(d, range=hist_range, bins=number_of_bins)[0]\n",
    "    hist_max = np.max(hist)\n",
    "    binned_data_sets.append(hist/hist_max)\n",
    "\n",
    "binned_data_sets = np.array(binned_data_sets)\n",
    "\n",
    "x_locations = np.linspace(0, len(binned_data_sets)+1, len(binned_data_sets))\n",
    "\n",
    "i=0\n",
    "for x_loc, binned_data in zip(x_locations, binned_data_sets):\n",
    "    non_zero = np.nonzero(binned_data)[0]\n",
    "    effective_range_ids = non_zero[0], non_zero[-1]\n",
    "    effective_bins = non_zero[-1] - non_zero[0] + 1\n",
    "    effective_range = hist_range[0]+non_zero[0]*(hist_range[1]-hist_range[0])/number_of_bins, \\\n",
    "                        hist_range[0]+non_zero[-1]*(hist_range[1]-hist_range[0])/number_of_bins\n",
    "\n",
    "    bin_edges = np.linspace(effective_range[0], effective_range[1], effective_bins+1)\n",
    "    centers = 0.5 * (bin_edges + np.roll(bin_edges, 1))[1:]\n",
    "    \n",
    "    heights = np.diff(bin_edges)\n",
    "    lefts = x_loc - 0.5 * binned_data[effective_range_ids[0]:effective_range_ids[1]+1]\n",
    "    rights = x_loc + 0.5 * binned_data[effective_range_ids[0]:effective_range_ids[1]+1]\n",
    "    \n",
    "    ax2.barh(centers, binned_data[effective_range_ids[0]:effective_range_ids[1]+1], height=heights, left=lefts, color=palette[i], edgecolor='#595959')\n",
    "    \n",
    "    ax2.hlines(y=means[i], xmin=min(lefts), xmax=max(rights), linestyles='dashed', colors=palette[5])\n",
    "\n",
    "    i += 1\n",
    "\n",
    "sns.despine(ax=ax2)\n",
    "\n",
    "# ax2.text(x=0.01, y=1, s=r'B', size=15, transform = ax2.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "# ax2.text(x=0.01, y=1, s=r'D', size=15, transform = ax2.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "ax2.text(x=0.01, y=1, s=r'F', size=15, transform = ax2.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "\n",
    "ax2.set_xticks(x_locations)\n",
    "ax2.set_xticklabels(labels)\n",
    "\n",
    "ax2.set_ylabel(r'Percentage of tests saved', loc='center')\n",
    "ax2.set_xlabel(\"Number of contacts\")\n",
    "ax2.get_yaxis().set_major_formatter(matplotlib.ticker.PercentFormatter())\n",
    "ax2.set_xlim([-1, 7])\n",
    "\n",
    "fig.tight_layout()\n",
    "# fig.savefig('figures/fig1.tif', dpi=300)\n",
    "# fig.savefig('figures/figs2a.tif', dpi=300)\n",
    "# fig.savefig('figures/figs2b.tif', dpi=300)\n",
    "fig.savefig('figures/figs2c.tif', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "e312f55c3c97e43ef328a86651ae8cf61ffd8a19c7888f1a2843b6f17cf5ecdd"
  },
  "kernelspec": {
   "display_name": "Python 3.9.1 64-bit ('env')",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  },
  "metadata": {
   "interpreter": {
    "hash": "c425b3d18afb7ef01beae38734bbc4c3f66a0926b3e7f0a803887fb6dc8ff96a"
   }
  },
  "orig_nbformat": 2
 },
 "nbformat": 4,
 "nbformat_minor": 2
}