notebooks/score_vs_R_k.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance comparison for various values of R, k and N\n",
    "\n",
    "This notebook reproduces Figs. 2, S3, S4, S5."
   ]
  },
  {
   "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.pyplot as plt\n",
    "import matplotlib as mpl\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",
    "sns.set_theme(context='paper', style='ticks', font_scale=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the values of sensitivity s<sub>e</sub>, specificity s<sub>p</sub> and number of contacts N."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=100\n",
    "se=0.9\n",
    "sp=0.99"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "name=\"score_vs_R_k\"\n",
    "k_values = np.arange(0.05, 1.01, 0.05)\n",
    "width_pt = 469\n",
    "palette = sns.color_palette('husl', 6)\n",
    "type_dict = {\n",
    "    'false_negatives' : int,\n",
    "    'false_positives' : int,\n",
    "    'k' : float,\n",
    "    'lambda_1' : float,\n",
    "    'lambda_2' : float,\n",
    "    'method' : str,\n",
    "    'avg_group_size' : float,\n",
    "    'num_of_infected' : int,\n",
    "    'num_of_tests' : int,\n",
    "    'r' : float,\n",
    "    'N' : int,\n",
    "    'score' : float,\n",
    "    'se' : float,\n",
    "    'sp' : float,\n",
    "    'seed' : int\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 20/20: Reading file 192/192\n"
     ]
    }
   ],
   "source": [
    "results_list = []\n",
    "\n",
    "for ind, k in enumerate(k_values):\n",
    "    files = glob.glob(\"./outputs/{name}_*_N_{N}_*_k_{k:.2f}_se_{se}_sp_{sp}_*.json\".format(name=name, N=N, k=k, se=se, sp=sp))\n",
    "    assert files!=[],'No files found'\n",
    "    \n",
    "    dicts = []\n",
    "    for fl_id, fl in enumerate(files):\n",
    "        clear_output(wait=True)\n",
    "        print('Iteration ' + str(ind+1) + '/' + str(len(k_values)) + ': 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)\n",
    "\n",
    "    total_df = input_df.astype(type_dict) # Fix types\n",
    "\n",
    "    temp_results = total_df.pivot(index=['lambda_1', 'lambda_2', 'se', 'sp', 'r', 'k', 'seed', 'N'], columns='method', values='score').reset_index()\n",
    "    temp_results['Percentage'] = (temp_results['binomial'] - temp_results['negbin'])/temp_results['binomial']*100\n",
    "    temp_results = temp_results[['r','k','Percentage']]\n",
    "\n",
    "    results_list.append(temp_results)\n",
    "\n",
    "results = pd.concat(results_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# results = total_df.pivot(index=['lambda_1', 'lambda_2', 'se', 'sp', 'r', 'k', 'seed', 'N'], columns='method', values='score').reset_index()\n",
    "# results['Percentage'] = (results['binomial'] - results['negbin'])/results['binomial']*100\n",
    "# results = results[['r','k','Percentage']]\n",
    "results = results.groupby(['r', 'k']).mean().reset_index()\n",
    "results = results.pivot('r', 'k', 'Percentage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mapper(r_quer, k_quer, r, k, Z):\n",
    "    ind_r = 0\n",
    "    while r[ind_r] < r_quer:\n",
    "        ind_r += 1\n",
    "    ind_k = 0\n",
    "    while k[ind_k] < k_quer:\n",
    "        ind_k += 1\n",
    "\n",
    "    return int(np.around(Z[ind_r, ind_k]))\n"
   ]
  },
  {
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV8AAADeCAYAAACaLgNIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABc1klEQVR4nO2dd3iT1fv/31lN0r1YLW0ZhbJBNhQQBVRQ2SooKKKiKBv1pyh+QBHnF4UiAgKyUaasAiKjZW9aRltoS/cu6Uqzc35/1DwmbcaT1abteV1XL5tnnHM/wb5zcp97cAghBBQKhUKpVbh1bQCFQqE0Rqj4UigUSh1AxZdCoVDqAH5dG1DbEEKwe/duxMXFoV27dpgyZQoEAoHF+8rLy3HkyBGkpaXh1VdfRVhYmNVzHz58GFeuXEFwcDDefPNNiEQiWx6BQqE0ADiNbcNt/fr1aN26NYYPH447d+5g06ZNeP755zFixAiT9xBCMHfuXLzzzjsICQnBb7/9BrVajdmzZ8PT05PVvHv27AEhBBMnTkRKSgrWrl2LQYMGYezYseBwOI56PAqFUk9oVOJbUFCAtWvX4osvvmCOEUJw6NAhxMTE4P3330d4eHiN+3bu3ImwsDBERkYyx7KyshAVFYVu3bph8uTJ4HJNe3AqKirw1Vdf4bvvvjM4furUKRw+fBhvvfUWunbt6oAnpFAo9YVGJb6ff/45PvzwQ/j6+tY4p1Ao8Ouvv6KiogILFy6EWCwGAEilUixduhTff/+90TEvXbqEP/74A1OmTEGfPn2MXrN8+XK88cYbCA4OrnFOrVZj48aNyMrKwocffggfHx+DsWNjYzFkyBAMGDDAhieuiVqtxpo1a5CXl4egoCC8/fbbdrk/rl69irNnzyIyMtLgw8kZZGdn4+jRo3Bzc8OkSZOo24ZSvyGNhJs3b5INGzZYvC4rK4vMmzePef3NN9+QjIwMs/doNBqyevVqcuzYsRrnUlJSyI8//mhx3sLCQvLee+8xr7du3Uo2b95MHj9+THbt2kUWLFhAUlJSLI5jiW+++Ybcu3ePEEJIcnIyWbhwIdm/fz/RarVWj7Vr1y6yYcMG8vjxY7J7924yf/588vDhQ7ttNEZubi6ZO3cuefToEblz5w6ZP3++zXZTKK5AoxBfrVZL5s6dS9RqNSGkSugyMzNJXl4eSUtLq3H9999/Tx4/fkwePXpEvvvuO+Z4cXExkUgkpKioiKhUqhpzLFq0qMZYCxcuJJWVlcxrc3Nv2LCBpKSkkKKiIrJ48WKDc5WVleSHH34gy5YtI2VlZda/CYSQpKQk8vPPP9c4fvLkSTJnzhwSFxfHeiyJRFLjeeVyOVmxYgX56quvSGlpqU02muKTTz6p8dw6u+Pj4x06F4VSGzSKaIc9e/Zg/Pjx4PF4AICcnBzk5+fDzc0NCoUCzZo1M/gKW1JSArFYjK+//hpffvklczwpKQlKpRI5OTnw8fHBiBEjmEgJpVJZY+Ps9OnTGDBgAOPCsDR3Xl4efH198fPPP2Pu3LkGY4nFYnz44YfIzs7G119/jc6dO2PKlClWbdatXbsW33zzTY3jw4cPx9ChQ7Fx40bs2bMHH3/8Mby8vMyOtXLlSsyZM8fgmFAoxPz585Gbm4tvvvkGEREReOONN+zeULx69So6dOhQwyad3Rs2bMDu3bsxb948BAQE2DUXhVJbNPg4X5lMhmvXrmHIkCEAgAcPHqB58+bQaDRISUmBTCYzEN6EhAQ0b94cV65cQZ8+feDu7g4ASElJQVBQENRqNSoqKlBZWWkQovbrr7/ijTfeYF6r1WocPHgQ48ePZ46Zmzs7OxuEEOTk5KBFixYmRSQ4OBjffvstWrdujWXLlrF+H44fP46hQ4dCKBQaPc/n8/Huu+9i9uzZ+N///md2rKSkJPj5+aFZs2ZGz7do0YIR3yVLlrC20RiEEOzYsQNTp041afd7772H+fPn44svvsDjx4/tmo9CqTXqeuntbH744Qfy6NEjQgghCQkJZPny5eT+/fvk4cOHRC6XE4VCQQghjO9w7ty5pLKyksyZM4c5lpycTFatWkUSEhJIamoqUSgURC6XM/fl5uaSL7/80mDe9evXk1u3bjGv9edOTk6uMffHH39MSktLydy5c4lSqSSEVLk5srOzSXp6Orl7926NZ/viiy+Ya82hVCrJ3LlzmefJysoiycnJJCEhgTx48KDG9cuWLSMVFRUmx5s3bx5juyUbly5dyrxXtrBjxw5y8eJF5rW5+W7fvk3+/PNPm+eiUGqTBr3yzczMhFqtRqtWrQAAgYGBGD58OEJCQkAIQUZGBtzc3AAAHA4H0dHReOaZZ7Bz506Dr8sikQj9+vVDkyZNUFFRgYKCAmYFyeFwsGrVKgM3gUQiQXp6Onr06MEca9GiBZ5++mk0b94cKpXKYO6rV6+iU6dOOHPmDJ5//nlmRa1SqXD//n1IpVJcuHABqampBs+n0Wig1Wotvg+//fYb3nnnHeZ5AgMDkZ2djYcPH2L37t24e/euwfUSicTAVaJPdHQ0hg8fztjuKBuNUVFRgbi4OINID3PzyeVyxrVEobg6DVp8V69ejdmzZwMAiouLUV5ejmbNmiErKwutWrViRBmo8tmePHkSAwcOREpKCnr27AkAiImJQWpqKrp3746AgACEhoaiadOmzH3Xr19Hu3bt4O3tzRxbuXKlgRhfvnwZWVlZ6NWrF/z8/BAeHs7MrdVqsX37drz88suIiYlhkj2KioogFAqhUCigUqkQFhaG5s2bM2MmJycjICDApBtBR1FREfLz89G5c2cAwN27d/Hw4UNoNBq0a9cOXbp0MXgfzp07h169ehmNW9a9R6NGjWJlY1paGjw9PU0KuSX0//0szUcIwa5duzBmzBib5qJQapsGK74XLlxA9+7d4eHhAQC4d+8e9uzZgxs3buDRo0fg8/kGPtt169bh3XffNRBOlUqFGzduICYmBlFRUbhx4wZ8fHyYVR8hBFu3bsXrr7/OjHP//n0EBgaiSZMmAKoEq6ioCHfv3sWGDRtw9OhRaLVaZu6dO3fi1Vdfxdq1a/H+++8DqFrxnThxAiUlJfDw8IBQKERkZCTjfwaqfMwzZ860+D7oP49CoUBRUREAwNvbG02bNsWYMWOYLD2tVou9e/di0qRJRsdav349ZsyYAQ6Hw9hYXFxs0sZffvkFH3zwgUUbjZGWlgYej4eWLVsavCcSicTofAcPHsQLL7wAPr9R7CFTGgANUny1Wi3+/PNPTJ48mTnWunVrREREoE2bNsjPzzfYgS8oKMDjx4/B4XDg6+vLbCQJBAIMGTIEQ4cOxZtvvgmZTGYwzx9//IGXX36Z+apLCMH69evx7rvvMtdIpVJ4eXmhRYsWGDp0KBQKBSPeFRUViI+PR6tWrVBRUcFk1yUlJUEgEMDT0xNdunSBp6enQRrzP//8g8GDB1tMMoiPj0fLli3h7+/PvBYKhWjVqhXatm2LwsJCg+t1HyTGohMKCwtRVFSEjh07AqjagPTw8ECLFi2M2hgTE2OwYWktv/zyC2bNmsW8Tk9Ph0AggJ+fHzp37gwPDw9mPrlcjnPnzmH48OE2zUWh1AUNUny3bduGqVOngsPhIC0tDf/88w88PDyQk5MDlUqFp556CkCVWAL/rQ7XrVvHrCYLCwuRl5eH3r17Y9CgQXjw4AGCgoKYOaRSKW7duoVBgwYxx44ePYpnn32WWdXm5eVBJpPhySefRKtWrZCbm4t+/foxc0dFRWHOnDkGPuO8vDzI5XK88MIL0Gg0zMpdh1qtxpEjRyx+vSaEYOPGjXjrrbcAALm5uUhISEBZWRkkEgl8fX0NXARlZWV48OABevXqZXQ8/RV0aWkpKisrERcXh9zcXHh7exsItkajwYEDB/DSSy+ZtdEUMTEx6Nu3r0GWYWBgILp06YLCwkIEBAQYuEXYfgugUFyJBie+ZWVlSExMZFJ9r1y5glu3buHmzZsQi8Vo164dU5GMw+Hg1q1baNWqFS5duoSnn36aWZXevHkTW7ZswYkTJyCTydCtWze0adOGmaf6ykypVOLUqVMYOXIkc+zu3bv44YcfsGPHDvj7+yMyMpJJMU5LSwOfz0dubi7at2/P+Ixv3bqF7du3IyYmBpcvX4ZYLDZIS96wYQPeeusti7Gzf/31F0aPHs18DSeEICgoCPfv38fKlSuRmppqkMpsLG5Xx+3btxEaGgo/Pz8AVS6dkJAQjB8/HjExMVAqlQYfTJs3b7Y5vlej0WD//v2YOHEic6ygoAD//PMPysvLsW3bNhQVFTHz5ebmQiqVGq3JQaG4Mg1OfHWrSR0jRozARx99hB49eqB9+/YGG1SEEGzevBlTp07FiRMn8PzzzzPn+vTpgylTpkAikeCrr75iVskAkJGRAUIIQkNDmWNr167Fe++9Z2DL8OHDsWLFCgDA7NmzIZVKmXOrV6/GBx98gG3bthnEB7dq1Qr9+vXD008/DX9/fyiVSuZccXExsrOzLRbhkcvlOH/+PIYNG8YcCwoKAp/Ph0KhwHPPPYfAwEDm3MOHD+Hl5WWwEtZ/j37//XdmBQ0AXl5eOHv2LLy9vdG+fXuD60tKSpCamoonnnjCrI2m+P333zFt2jQD4eZyubhy5QpEIhFGjx5t4MpYtWqVyQ8NCsWlqZsIN+fw8OFDsmLFCqPnEhMTyenTpw2O7d69m5w5c4ZERUUx9Q6qU1xcXCPt9sMPPyRSqZR5nZeXR5YuXVrjXo1Gw/yuizUmhJDTp0+T3bt3kx07dpDz588bnTc5OZkcPXrU4NgXX3xBioqKjF6vz//93/+R5OTkGsczMjJIQkJCjeMLFiwwGYu7Z8+eGu+bXC4nFy9eJPv27ath/9KlS0leXp5FG41hLGVZR0pKCtm8eTM5c+YMc+zq1avk999/t2kuCqWuaVDiu2DBAiKTyZjXukQAHfpiWFlZSRYuXEgKCgrIF198wXqO2NhYsmvXLoNjixYtIhKJxOj11Qu/qNVqMnv2bFJeXk4++ugjs3PpalEQQsjdu3fJmjVrLNqXnZ1Nli1bZjB3SUmJyWSMEydOkIMHDxo9p3uPLNmmIzExkaxcudKijaZYsmQJyc/PNzimUCgM/t10z6TVasns2bMNzlEo9YkGI74nT54kBw4cIIRU/WGqVCoSHx9PtFptDREmpKp4zqNHj8hnn31GHj9+zBzXarXk0qVLRufQCae+oBqrlqZWq8n169eNjrF+/Xpy8+bNGtXSzN2j1WrJvHnzWGWz6RegUSqVRKlUkqSkJKJWq0lJSYnBtdUz36rzww8/kNTUVAMbr127ZnLu+fPnG32v2ZCQkEBWrVrFvLb0b7h9+3Zy4cIFm+aiUFyBBuHzrR4BEBsbi2vXroHD4UAulyM9Pd3g+qysLKjVapSWlhpsJN2+fRscDgetW7dGSUkJEhMTDe7bsmWLwUYS+ddnPG3aNOaaS5cuMQVyysvL8eDBA+acRCLBo0eP4O/vDy6Xi5CQEIv3AFXth0aOHGmx3dHVq1fRsWNHpgDNmTNn8Ndff6GsrAyEEBQUFBhcv2HDBrz99ttGN8ays7OhVCrRunVrAxv9/PyM2njs2DGDDUtrqR6iZ+7fUCqVIi4uDgMHDrRpLgrFFWgQ4qsfAaBUKqFSqUAIQXFxMe7evVtDtKKiojBr1ixs2rQJ06dPB1AljAcPHsTWrVtRUVEBT09Pg0y20tJSJCcnG4RiVa+WJpVKkZCQgNu3b6OsrAwlJSUGsa8rV67EvHnzsHr1aiZSQiqV4sGDB7h9+zbKy8tr3KNQKHD27Fk888wzZt8DXabclClTmNfNmzdH165dkZOTg+PHjxuIbFFREXJzc9GlSxej4+lvZOmeKy4uzuhzqVSqGhuW1hAdHY1hw4Yxwq1UKpmYamP/hlFRUQaZbxRKfaTei2/1CIAjR46gd+/e6NmzJ5o3b47CwkLk5OQw11+8eBHdunXD33//bRCKJRKJ8Pnnn+PZZ5/FhQsXsHv3boN5qu+qV1ZW4tq1a3jyySeZYwkJCZg2bRoiIyPx+PFjXL9+HRUVFQCAxMREBAQEIDExEb1792Z27Pft24fRo0ejf//+KCgowI0bN5h7gKooCjYxrLt27TJoZxQbG4vw8HB06NABw4cPh7u7O+7du2fwPPPmzTM61qVLl9C1a1dGYO/du4c33ngD/fv3h0QiMXgu4L9Vqy2hZdVTlgHg2rVrePrpp9G/f380b94cRUVFzL+hLvNN962BQqmv1Hvx1RcRiUSC/Px8fPnll9iwYQPatWuHpk2bMiFUWq0Wf/zxB8aOHYuLFy8yoVjp6ek4d+4c/vjjD5w7dw6vv/46lEolsrKyAFSVgvT29jYIxaqeOhsfH48DBw7gzJkzuHHjBrOS0xWVWbduHd5++23s27cPL7/8MoCqOgpJSUm4fv06PD098cwzz4DP5zP35OXloby8HO3atTP7Hugy5XQFaOLj43Ht2jUcPXoUBw4cgLu7O1q1aoXevXsDAO7cuYOgoCAm800frVaLXbt24dVXXwVQtUJ+/Pgx7t27hzt37uDpp5+GUChkbCwsLERhYSGT+WYt+inLOtt2796NI0eO4MSJE4iIiIC/vz9atGhh9H2nUOor9ToR/t69e2jevDlT+9bPzw9vvfUWdu/ejZs3b+L48eMGK6pt27ZhypQpWLt2rcEfcEBAALKzszFu3Djs3r0bSUlJBn7c6kXIMzMzodFoDArS6Py4QUFBSEpKgkQiYb6G6yqB6VdL02q1yMzMxOeff84khnTv3h3PP/88I0SrVq3CJ598YvF9qP413MPDg1kFHzlyBJmZmUyCCPk38+3HH380OtaOHTvw2muvgcvlQqvV4tSpU9BqtQgLC8ODBw/g6+tr8J6uWrUKCxYssGijMQoLC1FcXGwg3LqCR506dcKdO3dQUFDAZAXam7JMobgSLiW+0dHR+OOPP1BRUYHWrVvDzc0Nbdq0QVZWFhYuXGhQOYwQgg0bNjCNLRMTEyEWi5GTkwO1Ws0kN+goKyvD/fv3MXz4cCgUCkaMMjMzERISgsLCQuzduxcdO3Y0yJY6ceJEjSLkq1evNuiArBtj5syZCAoKQkZGBmOrUqnE33//jSVLluCHH35gkhUKCwuZ1aVWq0V+fj7UajXjBjFWLc0Yjx49gkAgYArQPHjwAEKhEC1btoRcLkeXLl2YVSNgvgBNeXk57t+/zxQuLywsxCuvvILr16+jpKQEPXr0MEjwiIuLQ0hICLNhaS3VhTsxMRHu7u549dVX0axZM2RnZzOrc41Gg3379mHlypU2zUWhuBx1GmvxLyqVisyePZu0b9+eREVFGcRu7ty5k0RERNToCXbw4EGmYWV6ejp5++23mde6RAL9PmvLli0j2dnZ5OOPP2ZCscrKysj06dPJ7t27ydmzZ0lCQgLJzc1l7jEWinX+/HmyY8cO5rVujO3bt5Po6GhSWFhokOCwatUqcv/+fYPkg9LSUjJ9+nSyefNmEh0dXeP90MWwGoulrc6HH37I9IjTvQ9HjhwhCoWC3L171yCcTS6XkwULFpgc6+uvvybZ2dnMc7355ptk586d5OTJkyQ5Odmg8LpWqyVz5syp0cuOLbdu3SLr169nXuv/GyoUChIfH2/QMHTDhg3kxo0bNs1FobgiLuHz/emnn3DixAn07t0bs2bNMiiaMnnyZEyYMMHgekIIYmJi8NxzzwGo2o2fOXMm+vfvj0OHDjFhULrVXUpKCsRiMSQSCdq1a8eEYnl5eeGzzz7DrVu3sGrVKkgkEgO/bvUi5ACwf/9+g2ppujESExPx22+/ITk5GW3btgVQVZOguLgYPj4+cHd3Z6qleXt747PPPkNKSgp+++03nD171uD5/vzzT4Nqaaa4ceMGevToYVCAZubMmRg8eDCOHz/OrMh1mCtAU1BQwLhNdM/1+eef4969e1izZg2Ki4sNfM/79+/HmDFjbC7huGPHDibSRN/2AQMG4Pjx48jJyWG+nehSlnU1limUhkCdux1KSkqwfft2ADDwJerz3nvvGZRPrKioYPy8eXl58Pb2ZvyGly9fNqivCwBr1qzB119/jT179jCxoXl5eeBwOAgICECfPn0wffp0A19i9SLkQJV7QCQSMWKsP0bPnj3xxhtvGBQOX7VqFRYuXIirV68yhX7073niiSdq3FNZWYmbN2+arKmrz40bN5gQNEvvg6UCNPfu3WPC6PRt7NWrF6ZNm2Zgo1wux6VLl0z6jdkgEAiYDxdLttP6DZSGSJ2vfOPj4yGXywHAZPhQSEiIQfC+l5cX8vPzUV5ejosXLyIjIwOVlZVQKpWYPHkyOnTowFz7zz//YNCgQRCJRHjqqafw559/ghCCixcvIiUlBUqlEp07d0Z4eLhBZa7q3SiAqgIvSqUSxcXFAGB0DF0Fstu3byMsLAx+fn4YOHAgDhw4AK1WW+Oetm3bGsyrHwNsiaeffhq7du1ibDH3PlQvOFQd3bcGYzbqPxdQ9WGmK/xuD3l5eRZtT0pKMqixTKE0FOpcfCsrK5nfiV7lMEt89tlnWLRoEQoKClBaWooLFy4gLy8PXC6XGUetVuPw4cMYO3YsAKBly5YYNmwYFixYgJCQEGaDSSQS4f79+9BoNACqPhCCg4ONhmJ9/vnn+L//+z9s2bIFTz31lNExCCEGCRxeXl54/fXXMX/+fAQGBpqcNyMjAwAMqqWZIzw8HH379sWCBQvQunVrlJWVGX0frl27hoiICLPt4MViMd555x3Mnz8f/v7+Jm3Mzs6GXC43KK9pC5999hl++eUX/Pbbbxg6dKhJ241Vi6NQGgJ1Lr76fsTMzEzW9zVt2hSrVq0CUJXaClTVmc3NzWXcArpi4vo+2yFDhuD777/H1atXcfnyZQQGBuLChQvIyckBj8djQrHefvtto/P6+Phg+fLl6NixI5YsWQIvLy+0adPGYIy9e/di3LhxBj7bnj17YsWKFUhOTkZsbCy8vb0N7gGsW/XqGDZsGL755hvExMTg6tWraNmypcH7QAjB9u3bTbZe16dbt2746aefkJ6ejrNnz8LLy6uGjY7KLvPw8MBXX32FPn364Msvv4RQKESHDh0MbD927JhB5huF0qCos60+PaZNm0bat29Ppk+fbvR8cXGx0eO5ublk9OjRRKVSkfXr15NFixYxu/UZGRnkf//7n9l5CwoKyOeff07WrVvHlIjcvn07OXXqlMl7CgsLmapmGo2G7Nixg3z00UckPT2dKBQKUlRUZLRa2enTp5lW5yUlJeTLL78kP/30EykvLyeEEHLs2DGyZ88es/YaY+/evSQnJ4cQUvV+fPrpp2Tjxo1MBMSGDRsMWq+zpbS0lCxbtoz8+OOPTHTIP//8Q3bu3Gn1WMYoKSkhW7duJYRURU7s2bOHLFy4kCQnJxOFQkFKS0vJvHnzTBb9oVDqOxxCrPiu7yQKCgowZcoUpKenY9myZQbtZ27cuIFjx47h888/r3FfTk4O3n33XRw+fBhAVYbbr7/+ioqKCnh4eGDhwoVG+5xt27YNnp6eGDduHICqeNVdu3aBx+Oha9euZje7pk2bhrFjxzKuDKBqA3Dt2rUoLi6GSCTCvHnzDLpElJeXY/DgwTh16hSzUQhUxeRu3rwZHA4Hbdq0MShYzoaMjAyMHz8eFy5cMIhDvn79Ovbu3Qsej4eePXvWiBYxx+nTp5GYmMj4dFNSUrBx40ZwOByEhYVhxowZVtloilmzZqFfv34GK3KZTIZ169YhPz8fbm5umD17tkHRdwqlIeES4gtU+X43b96Mv//+GwKBAMHBwRAKhejUqRNee+01oyFN1cWXLVFRUfDy8jLIYmNDaWkpRowYgatXr1p13/bt25GWlmb0A8Qevv76a4SGhrJyKbDlwIEDuHfvnsNt1Ucmk2HAgAG4deuWTfUgKJSGgFNDzW7cuIE5c+aAw+Fgy5YtTPyrMdzd3fH+++9btYseEBBg0jfrDB49emQQm8qWgoICvPLKKw63RyqV4oUXXnD4uM4mMzPT5h5vFEpDwanie/XqVZw/f95pf2RSqRTR0dEWO/lWp3379jbVB7h8+bJBggUbCCFQKpUWi+NYS3l5OYKDg21O7TVFaGioVVEntnDhwgUmtZpCaaw4TXyLi4vxzz//YM+ePfjqq68QGRlpcD4qKgqrV6+ucZ+bz2hoPP6r3aoR//e72uB3HtSVZXh8NRbR73xnMAZXbBjEIRAZin9lejJEHgK4nzR0H4hENT8kxHrH7q5fh9/zC+Dm5WVw3BRiEQdquQKXN/2GS03FFq83hntNlzUAoDQrD/cP/YPzfmqj5z1N3GeJ4rRsaEtLcF7IPvLEWg5v+B3ntKnwbW7oz/UR0pWws/ARWr7GlZg32P6kmogOPQAiM3ouODgYp0+ftnsOe3Ca+AYEBGDfvn14+PAhZs+ejd27dxsUiZk9e3aNkKWIiAhnmWOAqjgPHLkI7i1bWXWfVq0G10I3iepoVCrw3Ky7py7HLc8vhiQzD236drZ8sY1oVCrwjdheqiBUgB1EfRNbp0BkcPMZbfRUdvahWjamJk6P823Xrh0mTJhgVQwvW7gCIYSB7BISHAFfLAaHZ91bxuXxIBDbtuo1B89NAL6wfsa/urmLwDNRE6JUQVCqcIk9YArFqTht5UsIYXy9AoHAZE0Bu+BwwOVbL0AcgRBcAbv7ZHLCuBhaDIwEl8evcdwcPKEbwp4cYLWNlnD390Vovx4OH5cvdIOnp3NFvcuI/nAz5U/5F7oKto9ShelzdFXsGjht5Xv8+HFMmTIFmzZtQt++fQ3iUB2FVimHLC/Z6vtEQWEQt2hZ47hcbnzFJfv3eOapf6BRKg2Oy0zcozuvqpQh+bjjfUvleYVIjbUu5I0NfiHNEdKtvcPH1efmwbOQlVVYvpDiFEoV1v1QnIPTVr4jR47EyJEjnTW8VajkxGDTTZb+ABoPEdz8AszcZYg5kW1IFCZnoDInD01aB1m+2MnQ1S+lIVPntR3sgcPjgyuove9QHMDqsDkOhwPYEWpXKTc9rrPCZBUa54z7HxzW7yP1Adc9dAXsHOq8nq89cAVCeLbtZfnCGnBgi26FPf88eFa6TwQe7ogY/ZwNs5nHO7g5wodFmjxfIbc93AwAyhWAl5M+1/q9MgLuvqYrrBlDX4DparhuqC7A1HdsH/V65auRS1GWeNHq+9zbdoJHG+vD2tKjo6FRWLcEUEkrkXTouNVzWaIsOw/Jpy44fNzmndqizcAnAFQJsP6Po7jy50lUlpQ7bkAKxQQaD4HRH7bcuHEDkZGRGDRoEFJSUrBp0yb89ddfTAOIpKQkbN68GT/99BNzT3l5Ofbv329x7HotvqbQyrRmz1em3Ic0NcnqcQkhRrO/zPmDCSGAEzLGqmxx+LDIu5+C1Iu3jJ5znAgbfx/ZQt0QrkFj2JjTZemeP38eEokEJSUlGDt2LEpLSxEXF4f9+/dj6NChyMrKgkQiAQAcOnSIVdZtvRZfDgc2uQ8AAlN/vqYiHij/4chVsK1QX7DrUd8iJSIiIgx+oqKiDM7rsnSHDRuGCxcuIDY2lmkiEB4ejpiYGPj5+UEul0Or1UIsFiMmJgaDBg2y2H8RqO8+X5EXvLs8VWvzNR04FDwrC3sL3MVoM/xJu+atlNdMM/ZsFohWkeb93fb6fZ1FjxcGQ+zl4ZCxaEQExVaSksx/+62epdunTx+mVKxQKERRURE++eQTnD17FhMmTEBpaSnkcjk8PDwQHR2N0NBQdOnSxeT49Xrlq5WXo+zuGavvE4e1g0dIa6vvK759HVqN8VoKptAoVci6csPquSwhLylDzq37Fq+rMBEtYYom4aEI7e281GIASIq9CZXcccsiugJ2XVx99csGXZZubm4uZLKqWhFSqRS+vr5wd3fHqFGjEBkZiZMnT+LZZ5/FypUrERkZifXr15sdt16LLyEw6j6wVFhHU5gBWV52jftEIo7R4jrMfZUVIBrz/uTqVEpVUJY7PqFArVBCWWm8aIg9SDLzkHvX+sQVa5CVSaFWWfchpo+PkFPjh+K61AcXhDH09yUEAgFmzpzJrJaTk5MxePBg5vyxY8eYvAaJRAKxWIzycvObyvVafG1BJOJAq1RCq1LWON5QsWb1q1YoHboqdQY6fy9d8dYv6psAV8/S7dWrF4RCIfbu3Qtvb2/06dMHAJCamgo/Pz+mS82kSZOwb98+9Opl3i1Yr32+PKE7vNr3r7X5/NpHgGuiIIypWg88Nzf4OriWLwCIfLwQ0JZ9USFH+3/tiQNu2TUcQgu1HdhiTIDpSth1KVXUbnywfklaazGWpWus2UObNm0MunkPGjSI1fj1Wny1aiUq0+MhbNvBqvuETZqBy7f+H0WlMr/SMiXAhNSsAcGmKI8+xjbdTIZs2IFXU39wuJa/ENmVgOGc6DuG6oJMxZjiitRrtwPRqKFRVFp/o1YLQqzz3QJAxaMH0Kqt3XBTQvKgpg/V3loR8tJyFKdmWHVPhdyyC4IQAqK1/N7YE26WdTcZSpmVO4EUSgOjXouvrSiKC6GUFDOvLfl7xSKOxZUqm2tcBXMCXFEogSQzr/aMcSJ0M841qS9xwM6mXrsduAIhRM0MQ8YsRTpYi76gunl51yimbklwuTw+3Ly8WF1rDXyhG9w8HF+k3Rps9fuKfTzBFzjvfz0quPUHfQFubLUi6o34Vu/fpoPoPUF14TUFTyQG18pkCQBo0rMnU0ydLR7ebgjt/4Tdwlvd3yvy9UaL7h2tGoPNhptALIKQpajb6vdtH9kDApFj/9Ko4NYfGpvImqLeiK8xtCoFFLmP4NXNugwyUfMgcKwUUQDIjjkLn7ZtwRex36lXVcqQeioWzbqzT1xgEwhQkV+E9Is30axTVSSFoyIZvJsH1ugy4ejqZnHR59GkbTA8/X1Y39NQxZUKUeOlXouvrUjTUsAVCiEMaFLXptTA2ggsR6cPFz/KgiQzD03CQ51WUtIWdBEMDU2ETfk9qSg3fOq1+HJ4fHDFnrU3H4d9EXB77rFmXFet32AeDjicRrnXS6ll1CbifF1B+Gz6C8jPz3e0HTbB5btB3Np04QqTcDg2lVMPGzXK6mLqap4Y7Z1WTH2gw8fV/7BwVvWyqmLqtfehSaG4IjaJ77p16xxth01oFJWQ3r9s9X1e7TrCo81/WWfmykjqx+OmHz0KtVxu8rwxVFIpEg6wL6Zuqm1QdUqzcp1STL1ZxzZoPbCHw8fV58qff9tcTL2xpBTXp9KMFNswK75HjhxB3759MXnyZFRUVEAikWDOnDnIzc2tLfucQvmD+6hIfWD1fbb/2TtHMJxWTP2C8WLqjsSeYuqNESrEDQ+zro8jR45g06ZNSEtLw+LFi3H//n3MmDEDEyZMqC37zFK9OaVWpmUdbkZhhzN7uVFsg/ZSaxiYVapBgwahS5cueOGFFxAYGIjNmzdjwoQJSEtLYz1BSkoKZsyYYa+dRuGKPOHVo/aKqYcOH2FDMXV3tBoxzOG2eDVvgjZP9nX4uLVBr3FPwd3Hdp8vrWhmCF0V10/Mrnx3796NW7eqvoJmZ2dj6dKlAIC0tDQcP27Zj6lUKnHhwgWmALGj0cjKUX7rNITPTrPqPs/wDqyKx1Qn5/w5eLVuBR6XvQCrFQpkxp5HcPf2Vs9njsrHJci4fAuB7Vo7NOKhafsw+LQwDMEzt/Fmy6o4/thFBLYKgtjbvm4WtKIZpT5jVnyffPJJg4LBOq5cucJq8H379uGll17CyZMna5yLiorC6tWraxzX7yxqTzk4U8jlBJVpyeC5CeHm62/VvWq5vKqYejWzTFUzAwCi0UCjUJi9pjpGK5hVQ6NUQa1QsRrPGoofZUOSmQe/0BasrtcXZrZCrJIroLGyQBFbGntFM2Or38bsltDPhtXHFULNzNrwwQcfQGQkm6tbt24WB7548SJ69+4Nsdh4qurs2bMxe/Zsg2MREda3c9eHjb9XJOJAqlZD+2+Du4ZcRN0WNCo1NDZ2mXBF/3BDTc6wBkvuiMYsznWJWfE1Jrzmjuuze/duFBdXVQ5LSEjAr7/+ipkzZ9pgoml4Ig94dq25MtdRvaiOvtAKBOyEV3+1GtC5M7gC61bjPKEQgZ2tq8HAJstN7OuNJhHW96FzBcJ6dnBYMXWK/dR2gXNKFaxW31KpFB4e1vnnfv75Z+b3qVOnOlx4AUCrlEOaeBWioLZW3Sdq1gIid+sL66grZdbHdxECVaXM4eUmtVot1DLH77B4t2gCnpvj3T36KKVyEG3tbZg15lUvW6gAGyclJQXfffcd1q9fj02bNsHf3x8VFRWYMmUKkpKScOnSJUgkEsyfPx8AUF5ejpMnT2L8+PEWx2YlvpMmTcKhQ4eYzCelUgk3G6qCORqi1YCorBcgrVIBLQt9qS6YpY9SEaxWM62E2NT4VahUKE9PsziXNQtBTxEgV1SgIjvH4anFGqUKKr1MD2e4EXKT0tBZ3q9GAR99qGCahopk7aAfMHD9+nWUlJRg+vTp+OWXXxAXF4fo6GhMnjwZUVFRkEgk8PPzw6FDhzBp0iRW47MS34EDB2LLli3w9fUFIQTXr1/H119/zfohtm3bxvpaZ6HvYlBKHoOjEsLLxLXmRFUs4oBv4TwbrBVba7BHMPNLSyAvKHS46OqLKZ8LeLtx4OWCAkuFrfFQfY9p1qxZBvtQ+gEDsbGxTJ+28PBwxMTEwM/PD3K5HFqtFmKxGDExMRg0aBB4POObfNVhvenn6enJZCWx8fk6ElO1fOHtDlHLqhAuRxdRN4bIzx8clm+sDh6fD7G/n13zGgslE4hEEPsa//gwFhrmShth3oG+4Aud69qwFVplrPGgawNvjOoBAxKJBD4+VSVQhUIhioqK8Mknn+Ds2bOYMGECSktLIZfL4eHhgejoaISGhqJLF/N1Z1iJ77vvvovU1FT06tULd+7cwdixY1k+npMhBFqlzOqsNr67B/hi6/+a/Dt1BNeC+FYPKeMKBGjaxboNN2NUF2CRtwead2xj+oZqsBVkoacYYicXvQnv3QECF3BbWYO5iAEqzA2P6gEDV69exYABAwBU7YH5+vrC3d0do0aNAiEEO3bswJQpU7B48WJ8+OGHWLx4MVatWmV2DlaqtXjxYly+fBkcDgfBwcFYvny5nY/mGIhKAWVBptX3uQU0hTAgwOr7ci5cQEWZdT5mlUyG9NiLNY6zLaCjj34DzPKCYjy6FGf9IBbwauqHJq2CHD6uPjePX4asXOrUOWoTmmHmumjFXKM/lvj555+xbds2bNu2DR07dsTOnTuZlXJycrJB/sOxY8eYFvMSiQRisRjl5ZYLR7ES3+HDh6NHjx5VN3C5OHbsGJvbag2tzLpOxJWZjyDNMt35197Ows7GUgdithhbDRel5iDjtvVFhyxB04Ep9ZlevXpBKBRi79698Pb2Rp8+fQAAqamp8PPzQ8C/i7lJkyZh37596NWrl8UxWbkdhEIhzp49i/j4eOzfvx+jR4+24zEcB4fPB8/T1+HFdExtmnG4XNuKqduQymxxXC4XHK7rbViZo1RB4CPkgMvlOuU9qWuo+6FhogsYeP/992uca9OmDbMRB1TVw2ELK/EdNWoU2rZti9TUVKxcuRKdO7PvR+ZMOHw3iILCra5mxuFyoFabvt5UKnDYyJFWN94UeLij/YvPWnUPG/gBzdBhhBOKqfOcK+qlCoInJgyHWuzR4FKB6WYdxRpYKVZcXBy+/fZbrFmzBn/++SdKS0udbRcrtPJKSB9ct/o+z/CO8GgdbvV9aUePQqOw7NjTd1soK6RIOmjcTWOL31dHaVYu4qMdX0w9qEMrtIvs4fBx9bny59+QSspqHNdVK2toLgpaGJ1iDFbi+9FHH2HKlCnYuXMnpk6dih07djjbLrswl1YMAPJH96HINO7XFIs4Ds9Gcyb2+n+r+32z76XiwTnnF1NvzFAhpgAs3Q5PPfUUhg2rqknr5eWFq1evOtWo6vBlKuON8DgcwErfoUjEQfU1lzViy+ba2hbv+tlE0zw6/3BDhxa9abyYFN+5c+dCIpEAAGQyGSZOnMgEHLds2bJ2rLMAV+wJr+5DnTqHvpC2e/E5k8XUTQmum4c72j//jMPt8m7RFG2f6u+QsWq7GlnvCcPg7mMqv5Cij06cqQjbhit3tjEpvuPGjUNkZCQERqp4KZVKpxnEk6oMavqaQ1tpWzF1n/Ydrc5UA4D0k6fhExZqIMCWVrlquRwpf59BQERbq1KKAfOrWUnBY2Rfuo6wzmFWjWlJZFt0bA2/kKZWjWkttw7FIiCkGcQ+ng12dUvFkmIJk+I7dOhQ5vfExETExsYyopuSkoKffvrJ6cY5i7KUh+AJhfD3tS7tV6NUgmj/iym2JLzuIkCh0gIalUXhtdZtoFFbrrtry2o2PzkDjzPyERDSvMY5RwklR6OCB09rd20HKnCU+gwrn+9XX32Fd999l6np4MyKZmxXvfZANGpoNdavfBsDWpXGpKi7mh+Wdm2g1GdYie+QIUMwZMgQ5nVoaKjTDLIGrtgDXt2frLX5mvXsYVUx9Uo5IBSJ0LKP5c4f1uLu74ugbu0cPm5t0GFgVwg9jHc4sRfa2ZdSX2AlvmKxGBMmTICHhwcIISgsLGTVQNPZaBUyVNy7CGFz6zo6uAe1BJdv/QpbVvwY0BLAikUz0WggLZJYPZclNEoVKh/XjJW1F7+Qpmbr7DqiLU9JvgRajXUp4RRKQ4OV+MbExGD79u1MtMPt27edaRN7tFoQdc0mkio5MVlSUi4nUEul4NrQraEsPQNajRpcQdXbxqYppkalguSR9cV/LKGokKIkK98J48qMJkBUxx4XRF5KFlRyBUQezo+Pox0aKK4KqziMl156CVq9jSZuPc/LV5ZKoCz7L0vP1Qvp2Iq5lu+mqJSUoyz/MatrG2I2GoVSW7Ba+a5duxbbtm1jRLegoAAnTpxwqmFs4AhFEId1qrX53Js1tTpEjSsQwKt5E4fb4uYuhmegr8PHtQV9AWazGvZvEQiBqPbq+er7gekquHHhjKYKjoKV+K5ZswZBQf/VeNUVGXYG1sT5QquFuqLE6jkEXt7gCqz/4/dr28bqalxcPg8B7RzfZdjNQ4yA1paTXfRXv2xCz0TeHvAM8LHZLnPuCLVSjajpy+Dh64WfXl2CjoO6YexHr9k8l232Gb6mYkypK1iJb1RUFPO7bsNt48aNTjOKLUSlhKo4x+r73Hz9wOGz7qDEkH3xMpp062Kxm4U+apkcaeeuolWk+fqe1qYIVxQ+Rvq1uwjtzb7CHBsh9vDzZm+ECUwJMN+ND79mAUiLT4bIyx3Dpr9g91z2QsPVKHUFKwUaMGAAUxy4pKQEN2/edKpRzkaamQ6uUAj35s7t2ODKmEopLk7PRXFGHlp2aWvX+KaiIkbNeQnr3vsebXu2h5cdK2xnQjfpKLUBq+/Qo0ePRnBwMIKDg9G5c2ecP3/e2XaxgsMXgO9jfTsga9DfjOPweFYXU5cpuVatlNnC4fLA5dm38WnLhpy1VN+UaxrWHJ2HPoFh0190/uR2UL0MJK1ERnE0rFa+L774Ivz8/EAIgVarZRrJmaOkpATffvst7t69i/fffx+jRo1iZZA1GW4cHh9uTUJYX8/cx+WBw7VeENs+/5z1xdTdxeg0dgSra61xPfgENUHEcMcU1tGHw+OCy3d8NIu+K2L8/5vq8PFrE1o0veFjTL82bdoEf39/VFRUYMqUKUhKSsKlS5cgkUgwf/58AEB5eTlOnjyJ8ePHW5yD1V/ZihUrsHXrVmzbtg07duzAG2+8YfGex48fY/ny5di4cSOio6PZTGOAxqBdvHFB1ipkqEy+bfXY3u0i4GVDMfXkQ+yKqeujrJDi7j72CSm6JpnVf6pTkpWHhBM1G3PaS4sOrdC2f1eHjwvQ0DRK/aG6fl2/fh0lJSUYO3YsSktLERcXh/3792Po0KHIyspiKkAeOnQIY8aMYTUHq5VvTk4Otm3bBpWqKqEhMzMT27dvN3uPrq9RXl4epk6tudKJiorC6tWrWRnpaEqT7oMrFCKwZ986md8WdALs7Lq9OfdSUZyRh36THF8GU4er1YigNFyqN1LQoQIQERFhcGzWrFmYPXs2gJr6FRsbyxwLDw9HTEwM/Pz8IJfLodVqIRaLERMTg0GDBoHH0s3ISnx//fVXzJw5k8lwu3fvHqvBMzMz8eOPP8Lf3x/9+vUzODd79mzmQXVUfzMswuVaXUzd0bDKcqOZtDVoiAJM44nrF7pW8KbQ1y9vb2/4+FRtEAuFQhQVFeGTTz7B2bNnMWHCBJSWlkIul8PDwwPR0dEIDQ1Fly5dzI7PSrlefPFFDBw4EH379kXfvn2ZrhaWCAkJwZYtW5CSkoLHj9llTVkDV+QOr26OKazDJsut7QujwBNa91fl5uAGmroVsG9wM3QYYdn3bona2HRrjJjbsKObdvUDff3i8/mQyWQAAKlUCl9fX7i7u2PUqFGIjIzEyZMn8eyzz2LlypWIjIzE+vXrLY7PauUbHR2N/fv3M4V1MjMzcfbsWVYPwOVy0b17d+ZTwx7U4v+W81oxF1ppGcrvnIW41ZvMcTYZLT4RnWwqpv7o+El0n/GmyW4WxlDL5Ug+dgohXdtYvtgKyvKK8ODMVTQJt7/CnE6AvYRAUOc28A9tZveYbLA2M64hwUaA6eq57tHp1wsvvIDY2FiMGjUKycnJGDx4MHPNsWPHMHLkSACARCKBWCxGeXm5xbFZie+XX36Jtm3/i/tMT0+3eM/mzZuRkpKCnj174pVXXmHtB7EXc0V1dJQ9TAJXKERAD/OJD9XRqlUGxdRZ3aPRQqtWo1IOqztZmKJCDmg1GmjVGscM+C/lCiD9ThpkeXnwn1g7AqyjMQuxKWhURd1RXb+6deuGK1euYO/evfD29kafPn0AAKmpqfDz80NAQFXI66RJk7Bv3z4mL8IcrMRXX3gBICzMcuuaadOmsRna4VjqXAwARKsB0RoKlym/rTl/LttGmTyuY4RXf7PNzqbFJjPcRFwtpGrnOampsNoOFd3aw5h+vf/++zWOtWnThtmIA4BBgwaxnsP6HFsXgiv2hE/vp0yery68utdlANz4plvEmzoeFtkHHl5u4PKMnzcmsEKeCKEDepq0EbAtgsEjwA8hPTsCqN3ml6ZoyKJKRY/iDFhtuGVkZOCff/4BAFy+fBlqtfneYbWFVlGJsusxVt/nERIGTxMdmM2tZkszs0G01sWpajUalGaarj9ha+iYWq5AaXaBw4U3IKw5gm1ILW4I8bs+QuM/FIozYCW+P/74IzIzqwqCd+/eHZ9++qlTjWKNVguisd7vqSyRQFFaavnCapRl5YBY+cGjValQmpVn8ryxBAo2KKSVKMt3fHW5ypIKlOYVOXxcV4eKbMNEJOIY/XEFWLkdhg0bxjiUi4qKcPGi4zOrrEX68CZKLh0Ghy9ARtTHaPneMtZlIlXlZeBprVtyEq0WhBDmv9bWeHAWpgrk2IqstAIVRdZ/MAHGV7+u7I6ggkupS1iJr6+vLw4ePIhjx47hzJkzeO+995xtl1H4Mg0TbiYKbgcOhwNtZRlEoRE21edli1omx9n/9zncvDxxbM4itH56MLq+NoHVvTyBAD4tWzjcJqGHO7ybO7eokCMw5Y5wBVGmSRGUuoSV+Pbs2RPdunVDZmYmFixYwKyC6xKeuxc8OvVHeVwM/AYbr5AllxOjXzHcfHzh5sH+r40vFsG/fTsUJyTCzcsTHSeyr8jF4fHgE2JefK2t5QsAfJEQ3i0c3yHD3c8LKoXS4eO6OjSsi1LbsPL5fvfdd7h+/ToiIiJcQnh1eHUfCnGb7hD4N7XqPr6HJwQeHkbPmcp06zhpIjh8Hlo9ORB8ofFVdqUR/61aLkfGJcfXP5YWS5B1K9Hh44o83eHhb39BdVP4CDkusepli6VMNZrBRrEVVivfxYsXQ6PRICYmBvn5+ejWrRu6d+/ubNsswuW7IWDoy1bfV5mTBa2HyGTEgzHEgQFoN3I42j5rOrQNgM3JFLasfnXoZ6jZy+PMfBRn5CG0e3v7BzNCffML2wpdSVMswUp8U1JSkJ2djX379kGpVMLf398lxLe26TB2JKvr9FfAKhUXhMsunFo/8sGSEHN5PHD5hlmD5uo0uEIssCnYhqk1RJGmNF5YqcL06dPx8ssvY9GiRayy2+oSNunFHB4fHH7tpDvzRSKEjxzGCDLbVbGllbB380C0f4p9SUy2q2OuoKaouwoNoRIabeBZu5iK27dcecH5sBLfHTt2GKQYK5VKuFnZ0cER6BfW0cEV/+e2ZpNaDADe4e2t7kJsK0ppJR4cPoEBC2Y6dNyS7HwknrzkkMI6+jQLD4Vv80CHjmkP9V1sLVEfGnjSDwznYFJ8v/vuO3z88cfgcDg4deoUjh07xpxLTk7Gzz//XBv2Ge1iof1XcK0RXv3XpUn3IfIQwad5zZWjLTUe2GCNH9jRBdOtcTnkJjxyWDH1+iqcVFwMoe+HczApvmPGjGESCTgcDnr16sVUJnN3d3eqURojgmts1WsKc8IrFnEg4Bu+NoY5sbVGSHmqqsI61bFHYL2EQKUAEPDs9+VWF0h3AVDBr7/CaQkqJBRXwaT4dujQgfm9R48eTAk1ACgrK3OuVS6ONRENbp4e6DLhObvmMyawAaHN0W3kQJvHbKjiagndV2gqwpS6xqzPNzc3F1u2bMGdO3cQGlrlWySEICUlBcOHD68VA23BVHKFDv/OnS0WU7fUHoitAKsqZbj310k8+dEMyxdbgSSrAHf/voxm4dZ3bzZHaNe2aNLK8Rl5roalWFwqzhRnY1Z8W7Rogddeew23bt1iigNzOByXSrSwBUlCAnhCIZr27u30uYiNxX/0MVa/QavVQuuE5nCZ9x6hMD0Pgyeza3ffUKGbTBRnY3HLPyQkBP7+/khOTkZwcDCSk5NRUlJSC6ZZRiuzTXyIVst0pDDXu41NX7fawtG91ky1cSdaLbRWdutoDNAMNoqjYRVvdejQIaZa++DBg/G///3PqUbZgspJQimTE4MffSrlxn/04YtFaD3E8S3qvZv6o21/891R2aATYd1PpQpQqE2LM6UKmlJcPxCLOEZ/LFFRUYE5c+Zg2LBhWLJkCQBgz5492LdvHzZs2ACtVov8/HysW7cO3333HbNgUalU2LFjByvbWIlv//79ERJS5VtMSkpCSkoKq8HrErmeUMrlhHktkxMIg1vBK/S/+FhzAlsdNtfoC7BWrUFh0qMa19hSx1d/9auQylCQnGn9IBYIbB2EkO7tmNf6wkwxD6350HCIi4vDt99+iyNHjuDSpUuIj4/HtWvXMGHCBAQGBuL48eM4ceIEunTpAj6fj4SEBADAwYMHMWbMGFZzsEqy6NSpEz7++GNkZmZCoVDg66+/tv2papHqG2+61/KiQigFbuD4/tckUv/TUF9c7Y3v1apUqMgvtGsMfXQCrJDKUG5j3V1zVBSVQJJdiKCOrR0+NsW8AFO/cu0SERFh8HrWrFmYPXs2ACAyMpI53r59e8TGxjLZveHh4di2bRsGDRoEuVwOuVwOb29vxMfHIzw8HJ6enqzmZyW+HTp0wPfff8+8Tk5OZjW4szGXZAGYznBTSyvAFdr3f7qlaIj6iry8EpUlxpMvG0J6r6tChbf2SUpKsnhNRUUFWrRoAT6fz4iqUChEUVERRo0ahRMnTqBfv37w9/fH1atXMWrUKBw9ehRNmzY1CM81BivxHTNmDHx8fAAAMpkMAoEAO3fuZHOrXeiy2/QTLKpnt5mq4+AqrUJ4AgH82zg2BdhLCBBvDwSE1m57d4pzoMLruhw8eBBz5sxBdHQ0Sv9tPSaVSuHr6wsej4dRo0YBAHbt2oWJEydi48aNGDJkCDZu3IgnnngCfL5piWUlvitWrGBqOygUCosO5YqKCixatAj37t3D4MGDGYe1reh3sGCDMeHVP+bmFwCuQGB1dpstWW8cLhe+TX1YZbRZk63GEwjg7uvF/gYW+Ag5aNbCF0KOxukrXCo4FFfnn3/+wfDhw+Hp6YnIyEisWrUKQNU3/8GDBzPXXbp0CX379oVAIIBEIoFYLAYhBGVlZfD39zc5PivxFYvFyMmp6sArk8lw8uRJTJ8+3eT1Omc1h8PB2LFjER8fj27durF6YHuxJLwA4O4lAldQ89GNias9acaeIkCuVCD96h10em6Q+YthOpzMmChXFJcgMz4ZHZ7sxRxzhGAKhAKIPMR2jwNQgaXUX3bs2IGNGzfC19cXKpUKb7zxBrp27Yo9e/agqKgIM2ZUJU09fvwYEokEAwYMAACMGzcOZ8+eRbNmzcwKL8BSfOfPn8+EmgkEArz//vtmr6/urG7SpGa7m6ioKKxevZrN9A6nJCMbIg8RvEItl8d0Bd8u2yaZjvDJPs4uREF6HsL7dLRrnCp7TJ+jwkypDWz9233ttdfw2muvWbzO39+fcT0AVftj+qUZzMFKfNesWWOQ1cZ2w03nrG7Roma66uzZs5mdRR3Vdx9twVJqcW3D5fFMth2yB56AD76R1Xt9gXZ6oDR2TP71Tp06FSqVCgKBACUlJRCLxRCJRNBoNJDJZNi/f7/FwXXOapuNk6mMlpTUh03x9Opw+AJwBebHZQOb+g58kRBdXnjS7rmq49PMHx2fqpkebe/ql+8mAN/N/veGQqGYx6T4fvLJJ+jcuTMAYPPmzZg2bRpzbvPmzRYH1ndWFxUVITDQdQp0e7UJh9iKDTx7UEorcXvfSbToHO7QcR9nFeDO35fRPMKxnUWatw2GTzPzvipnQl0VlMaCyQw3nfACVSlzCkXVX0VKSgqOHDlidtAdO3Zg+fLlmDlzJl588UWcPXvWJuMsrXoB06Fm5ih7kIDSWohVtiWLra7JSkhDwvm4ujbDKDSDjNKQYOU0fPbZZzFnzhwkJibC3d0dixYtMns9W2c1W8yFmbFxO7D1A7vC5pop2G66NXZoqUhKfYGV+IaGhmLZsmVIS0tDu3bt4Ovr6zSDjHWxqGvsEWWhlwd6TWLX9dgS+qFowuYtED7ceDF1U3UYaHYau9UxFWhKbcBKfPfs2YOtW7ciPDwcGo0GzzzzDF544QVn2+Y0fDp0BpfnmEgBS5tuiopK3Nx9As99/q5D5tNRkpWHxJOXENimZjF1UytkNqIc2rUtmjaCYurmoAJNqQ1YKVBxcTEOHz7MvN6zZ4/TDKoNSpPugycUwbtpzdxrU6tcm1e/hDC1gx0JIQREa1xM7XFRZN57hIK0XAx51f4Gmg0ZKtD1A2v6LdY2rEpKikQiEFL1h56eno7r16871SinQwgIoQXDjUG0WubfmmIfdAOQYg5WK98BAwbg9ddfR2pqKlq0aIGlS5c62y5WLeOtaRevQ9e9mGdjB2Nzn6TG6jdo+WJ0fPa/jD9HbZoJQgLAfbIH6/Fc1d9LV4eUxgor8T158iS2bNkCLpfVQtmhsCmow6aeg76Q+oS3BYfHt6qWA5uvL8ZCy1QyFTKu30OrflW1LYzVb7BFkGWlUmTdTUbYE/ZnBerTol0IvJv4mTxPxZJiDT4u/LW/rmElvqWlpQZfRR88eID27ds7zShnI83OAVfoZrQko6PDzbRqNcoKS8xeY4uPVimTo7KkgvX1+ptt5lbBj3OKUJSZj+AI42Uwaet1ijWU1sNY99qClfjeu3cPEydOZIoJP3r0COfPn3eqYc5ELZOB14ibROqE2JgIK6QyyMorWYzx3+9UiCkU62Elvl988YVBpZ60tDRn2VMDNrV82SRR2LuiZVPHwRg8NzcEtm+NCrlxn7AOa1e/7j6eaNom2HqD9GC7GrY8juVrqEBTKIawEt+KigpMnToV+fn56Ny5MxYuXOhsu5yKuEkTcN3caiWjjcPhQCB2vOOLw+OCL3JctTSdEPP8/NDECcEONPOMUhewaWJQV7DaQVu8eDHee+89HDp0CJ988onF2g6OgC9T2XW/3EyHYQ6PB46ZzUNL3YmtQa1QIDcuwWHj6ZA+LkPO/Zpdke2Fy+VCAV6tdyumdRsojQ1W4jt8+HBERkZCJBKhWbNmaNmyJQAgPz/fqcY5i8q8PMgKCmp9XkuFdkx1sqhNSvOKUZiaXfW7i7WLpwJNaUiwcjskJydj7Nix8Pb2hkajQVZWFv78808UFRXh2LFjTjGMTUUzazHmZrDG9WCL35fD40Hg/l9bHlMCrPt6ZE6A9X3CfKEAAge6HUxhbnPOFaH+Z4qjuHbtGlavXo0tW7ZAq9VizZo1CAkJgVarxbhx43DhwgVkZ2ejpKSEaSuUnZ2NhIQEDB8+3OL4rMT3zTffRN++fWscv337tnVPwxJ94dVttmnFjosx5gmF4Dmhu4Qx+EI3tH16gMXrqouyMV+VgTB7+aJl3x522WYMgcgNAnHN96a+ibA5aBcNChv69OnDlNI9fPgwmjRpgjFjxuCLL75A//79sWvXLvz000949913MWPGDGi1Wpw+fRpTp05lNT7rle+pU6fw6aef4syZMwgMDETXrl3Ro0cPmx+MDcaE11h2G9vMNh1eYWEQe7BPsrAnP9xNI8PD6NMI7+2YZAjd6jc/sxCPYq8ivHtro9fZKpKt2rVA02beDUJkLUHFtnFTvW3ZrFmzarQ2E/zb8SY2NhaTJ08GALRq1QoXL16En58fFAoFhMKq/5GOHDmCF198kfX8rMQ3NTUV/fv3BwA89dRTGD9+PKs2QrZgzt1QXXhF1f6rT3UR1X9d+OgBFEIRAlrWbOxpDWx2UuXK/36vjXq81opmdQFKT0lH3qM8RHS0L4yNUv9obNloSUlJrK+VSCTw9vYGAAiFQhQVFWHOnDk4d+4c3nnnHSQlJaFZs2aQSCS4ePEiOnfujLAw811mWIlv9+7dIRKJoNVqsXXrVmYp7gzY9G2rbSqruQR0K2Fj/tsagszhgMNx/CqSw+HA2LDGNsnMCTJNlqDooNlopvHz84NcXvUGSaVS+Pr6okmTJhg5ciSUSiX++usvvPzyy5g7dy5WrFiBhQsX4ueffzY7JitHardu3XDs2DFMmDABV65cwcqVK+1+mLqFA3BsDymrLsbmEHq6o+crzwGo8tnqfuzFP7Q5uo2MtHyhlXCc9GFBodQFXkLjP9by5JNPMivlR48eYeDA/xoZHDx4EKNHjwYAyOVy8Hg8lJaWWhyTlfiGhYXh66+/xoEDB/Drr7+iuLjYeutdiCZdOyOgo2ML0phCUVGJm38er3HcXgF+nJGH+GMXWF1bqiAmfwyvAwI7tEXbyB40jIvS6ElKSkJGRgYePHiA559/HpmZmdi7dy969uyJkJCqJgY3btxAly5dIBJVfeV97rnnsGfPHvujHTQaDWJiYhAWFoa2bdsCqIrt/fbbb3HgwAF7n81hWJteXHjnHvgiIYL693V4lluNNGJCTNbH1QmwLZ/EhBA4ouxudQFOvpmM4ow89Hvlv2LqOrcFzVKjNCYiIiIQExPDvF6wYEGNa3r16mXwety4cazHNyu+ixYtwv3791FWVoZ169YhMTERv/76a71PLwYIYKdw2VrrwRgu1RyToMZ7wybZwkfIoTG2FIoVmBVfX19fHD58GOXl5Xj11VcRFhaGnTt3IiAgoLbsq/e4eYjR5fknLV5nrQD7tghEx6G9LF/oYpQqqABTKIAF8W3WrBlUKhWEQiFeeeUVpmnmgQMHWC2v9TNE2KIWC2rE97IJMWMTs6t7HdQtAlx+VZyvPatXU6Fm+iKqqFQh/eINdBrYifW4bMLFigrKkX3rPjr3bsd6XDaEdwlDcGigzXG+zhRWKtyUhoRZ8V21ahV27NjB+Cx///13EEJQUlLCSnz1M0SciSmfrX40g1jEYV7nJ6WBLxSC6x0AmZx9YgXbCkn6m2lyqRrlkgpW/l1rBE8pV0JeIWN9PVsK0/NQlJGH5m1ti/Ot7npwpFi6qvA2tvhYimMwK74bN26s4VAGgPj4eNYT6DJEqhMVFYXVq1ebvZcr00Ir5kIr04Ir5kJlIjTM2IabsRoOOjQKBTjgGJyrfr0xn66lmrxscCn/rhFUcgUUMsd9YDaGVF4aH+u6H0CunKlpVnyNCS9QFfdrL7Nnz66Rylc93U8fnQADgEpOIBBxmLKRIr3fda+rr3r1UakAfi398fOFbgjqZthySX9lrC/E1Te2TP2PU1ZUgpL8x2jaqgUSzseh46DuDrO3tqDJHQ0L+gFkPawy3Boa7s2aQexZO4V1AABa09EC5uJ9y/8V4+or5WvHr+HyrhPgCnhQSuVYctJxSS++zQPA5VluWupIXE2Ija3WXcEuSsPCJcWXTesgeyBaLUgt9HBTK5QoTstG9p0HaNmzE/zDWtiUPVZdoNuMGISre09BJVdg0HuvWAwFs+arFyEEWo3GahsdhSOSOcwJpa3ju9oHBKX+49Re8PoZIq6ErLAQskLnZ+mlXYlH9JI1KHiQjgMLv4dS6pgNMjexCL1fHQUen4/2T/c1SFs2lr5sKbPN4Nr8xyhIz3OInXWFs4ut02w/iiNw6sq3eoYIG/RXvMZKSQKmy0la2nTTvdbtATo6yqE6rQf0wNWtB1EpKUfX5wYiMMDdtoGM0Gtkf4R1DIOPO7tvCI7eeKCrPwrFPlzO7WCP4JorI6n/u7u3Ozy8hHbF+BqLWKgpcAJM/Ox1/PnFBjw7fRQ8HSmAQjcEdDVfss6mYd1FEHua/5BojMLrSrv5xja3XMk+CjtcTnxrA49mTeDf3Jf19ca6TBjfKKv5db5Vl3C8MP8VePp5WWVjXeEf3ARCD7HZa+rqa3ddij7dzac4GpcTX11sr6PRj+VN+usY+s142eaxrIn35fK46FqP0oBzH2aiID0PwRGhdW1KDcyJfmNcjevD9sOhsa2QXfn/C5cTX330Y3vZYCxZgmi1Bm3iK/ILoVYq4dnEvvoUxgqp647pC3O53uaWqeQKVw4Ery+46iaYq/3x0xW86+DS4quPLrGCLWq5HElHTyDr3EVELvkMQm8vSAslyDx/Bd1em+BES013uDDmqvASWt99wplwuFxaTN2B0JhhiilcTnxNFdMB2G+0ldy5hYQ/9kCjUEDo7YmEHbugLCuHf/s26P/WRJNz25s6zAZXWf2aEoAuPdugsl1zo+fp134KxXG4nPgK/f4LnWIjtsZ+9+zVCYrcPsi4eBOqikoMeWcsfIObMplbtVVboS5Wr/aKYEpcMvIe5eGZN56tcS7Ux76xa5u6/IrdkHyrjdVVsWnTJvj7+6OiogJ9+vTBpUuXIJFIMH/+fABAeXk5Tp48ifHjx9s0vsuJrzWYqmbm5i5Gx5cnoN3oF3B741bA0xeVKh6gqvn135lCrHMn1KYI2+v7LJWqoOYJjGZ0NdY/Qh0NSVCtoTE+9/Xr11FSUoLp06fjl19+wbJly/DVV18hKioKEokEfn5+OHToECZNmmTzHC4nvtb6dvWpXkyHLxKi9wfvQA1A/a9w6De/dBcZDyPTx15xrm+baQWPchHcwTB+2FU3s2zF1m8HbD58GqNQNURiY2PRpk0bAEB4eDjWrFkDuVwOrVYLsViMmJgYDBo0CDw76qC4nPjqY6k3m7n+azohNtefjU3ZSHtLQLKtVOYKaDVatOoWjlY9HFug3dVwZpnLxv7twNUw92FYvYrirFmzmEqLEokEPj5VfjahUIhx48YhNTUVEyZMQGlpKeRyOTw8PBAdHY3Q0FB06dLFatucWtvBFbDUHp5NG3hHtHrXYU2dheqkxSVj26e/OrTerg6tRotDK3ahOLsQfIFLfyY7DdqpuXGRlJRk8KNf4tbf3x8yWVUtFqlUCn9/f4waNQqRkZE4efIknn32WaxcuRKRkZFYv369TfO73F8Zm1ZBpjbczGEuldhYlIMjfMGOWuVmJ6XjxNq/kJeSDR6fh9vHL0MgcoNAJERQ+xCkxydXXcjhILx3ByRfT4SutXFYt3DkPsyE8l/BDghpCq1GC0lOEQBAIBKCL+DhzNZj6PZ0b/QdM9hue2nkA6W+M2TIEMTGxmLUqFFITk7G4MFVfxfHjh3DyJEjAVStjsViMcrLy22aw+XE15jo6tC5GfRdCsbcC/pCa0v4mKM24UoVxCECLCuX4XFOEUAI3MRCEABqlRpcPg9Eq4VapQYAcDgcEC2BRqUGIQRiPuDJ18KNqEFI1TXuXC00RIuKf197CkXoNKAj+jzZBTx+7dbxpVBclV69euHKlSvYu3cvvL290adPH6SmpsLPz49pIDxp0iTs27fPZNMJS3CIrkGbCxAREQHNmHeMCq8p36+lsDPAdcTYHgghyLidhP2L12LmH8shslD8pjrmPgToSpXiaswbPMfuMSIiIvDJ9g+Mnvt2yi9ISkqyew57cLmVrypHBUW1RAsAkBtZEevaBemE1tTv+htrtvRhM5WZVptwOByEPdEB84/8bNP9Zmv41lmhHNfdfKRQnE2D33AzhrH0X2twhZVwQ8CazUYKpaHhcitf+PPBhemNt+qbbqbcC5YwtwKuLq76KzT6FZ1CoTgCl/P5UigUCgC7fbJPP/00srOzjZ4LDg7G6dOn7RrfXlxKfOsDERERde6odwYN9bkA+mz1kYb6XPo0Sp8vhUKh1DVUfCkUCqUOoOJLoVAodQAVXyuZNWtWXZvgFBrqcwH02eojDfW59KEbbhQKhVIH0JUvhUKh1AFUfCkUCqUOcL0MNxdCv4fTlClTmONHjhzBli1bUFFRge+//x5du3atQyutx9Rz6Vi0aBHGjBmDfv361YF19mHu2VJSUnDjxg20a9cOTzzxRB1ZaDumnu3kyZMoKSkBAIhEIrz44ot1ZKFtXLt2DatXr8aWLVsMjt+8eRM3b96EVqvFhAkTmGpiDQW68jWBrofT2LFjUVpairi4OABV1cVEIhH27NmD6dOnIyoqqo4ttQ5Tz6Xj9OnTqKysrCPr7MPcs6WmpmLPnj14+eWX66Xwmnu2rVu34qWXXsJLL72Effv21aGVttGnTx8oFDWrO61YsQJvvfUWXnjhBaxataoOLHMuVHxNUL2HU0xMDICq6mLDhw8HAHTt2hVNmjSpMxttwdRzAUBmZiY0Gg3atm1bV+bZhbln+/rrrxEUFIRly5bh+vXrdWWizZh7ts6dO2PlypW4c+cOXn311boy0S4EAoHB60ePHoHH44HD4SAoKAg3btyoI8ucBxVfE1Tv4VRUVFTjmkuXLuHNN9+sbdPswtRzqdVqxMbGYsSIEXVpnl2YerbKykpkZWVhypQpmD59OubNmwelUlmXplqNuf8f586di4yMDHz//ffo3bt3XZnoUPSfFwBKS0vr0BrnQMXXBNV7OPn6+hqcT09PR1BQEMLDw+vAOtsx9VzXrl3D4cOHMXXqVBw4cADLly9Hfn5+HVpqPaaeTaVSQSQSgcvlIigoCE2bNjX6YerKmPv/8YcffsDSpUvx1ltvYcGCBXVkoWPRf14AcHNzq0NrnAMVXxMMGTKEKeyRnJyMQYMGMb2aioqKkJiYiGeffRZSqbRe+UhNPdeAAQPwxx9/YNu2bRg3bhwWLVqEZs2a1bG11mHq2Xx8fODm5gapVAoA8PPzazDPBgAPHz6Ep6cnhg4dCrVaXZdm2o1Go0FFRQVatWrF+IEzMzPRt2/fOrbM8VDxNUGvXr0gFAqZHk7e3t743//+B4lEgrfeegvr1q3D+PHjMXXqVIjF4ro2lzWmnqshYO7ZFi9ejPXr1+Po0aN45513wOPVr3515p5t+vTp2LZtG/7++2+88sordWyp9SQlJSEjIwMPHjxAbGwsfv31VwDAzJkz8dtvv+HgwYOYP39+HVvpeGiGG4VCodQBdOVLoVAodQAVXwqFQqkDqPhSKBRKHUDFl0KhUOoAKr4UCoVSB1DxbaQsXLjQaYkGiYmJ+Oabb5wytquQlZWFMWPG1LUZlHoMDTVroFy/fh1vvfUW3n77bXh4eOD27dsYMmQIJk6cCKAqmN2Zsa7OHt/R5OTkICgoyKp7Jk+ejF27djnJIkpDh5aUbKD07t0bAQEBmDFjBoRCISorKzF+/Hj4+PhgxIgRThfG+iS8CoUCP/zwA3766Ser7uPz6Z8PxXbo/z2NBHd3d7z++uvYuXMnOnbsiPfffx+HDh3C9evXkZ6ejpiYGMyaNQt//fUXCCHIzs5GcXExNm7ciMrKSpw9exbZ2dlQq9WYOHEiPvjgA4wePRopKSkYNmwYZDIZdu/ejQ0bNmDLli3QarV47733EB0djcrKSiQkJGDMmDHw8PDA4sWLMWrUKBw+fBgzZ87E0KFDAVQV91mxYkWN+a9du4br169DoVDgueeeg0gkwtKlS9G5c2d4e3vDzc0NPj4+iI+Px6effoozZ87g3LlzaN26NU6ePIn58+fj3LlzyMjIwLp16yCRSAye57nnnkNCQgKOHj2KYcOG4ciRI6ioqMCDBw+wYMECfPjhh+jfvz+uXbuGjRs31nhvv//+e4SEhGDs2LH1KtuRUrdQn28jIjg4GPn5+WjZsiU8PDwAAJcvX4ZQKMSSJUvg4+ODtm3bokmTJli1ahWaN2+OI0eOYP369fDx8UH37t2Rn5+P4OBgeHp6YsKECVi2bBliY2PRokULLFmyBFwuF126dIFKpUJ5eTkOHDiAiRMnYsaMGfjss8/Qtm1byOVyvPzyy3j//fdx6tQpxj4+n290/uLiYrz00kvo1q0bLly4gA4dOkAmk+Hjjz/GBx98ALVajcmTJ0MoFOLBgwfo0aMH5HI5PvjgA7z00ktIT0/H559/Dnd3dxQUFNR4noiICDRp0gTPP/889uzZA61WizZt2kAgEIDL5aJJkyYYNGgQ1q9fX+M9PXfuHJ566ilMnjyZCi/FKqj4NiJyc3MREhIC4L+vzK+88gr279+P+fPnQyAQgMPhMBWzevbsifz8fCQnJ6NTp04YMmQIli9fDi6XCx6PBz8/P/B4PEybNg3ffPMNoqKiDHy9qampzO/NmjWDRCIBIQQeHh5wc3ODSCSCSqUysNHY/AMGDMCJEycgkUig1WoBAF5eXnBzcwOPx8OAAQMY36tufk9PTwCAWCxmftfNZ+x5dCQnJ6NNmzYYMmQIli5dCl9fX/B4POa/1fntt9+o6FJsgopvI0GhUGDHjh144403DI7n5eVh06ZN6NmzJw4ePAigSsAAoLy8HF26dEFYWBjzdTsmJgbV92ilUin27t0LkUiE2NhY5nhoaChSUlKg1Wqh1WrRunVrcDgci7ZWn3/58uXo168fU0y8Ov/v//0/TJ48Gb6+vjVsM4ax59H9hIWF4ffff4dGo8G1a9eYSmimWLJkCRYvXmxQ/pBCYQMV3wZKXFwcHj9+jC1btmDr1q344osvMHPmTAwcOBCFhYXIy8tDfHw8jh49irVr14LL5WLIkCEAgPPnz+PYsWPQaDR48sknMWPGDMTHx+O5554Dh8NBfn4+8vLycP78eQDA5s2bsX37dgQGBqJXr164e/cuHj58CLFYjHfffRfffPMNtm7dio8//hhpaWlIT09HWloa7t27h4yMjBqFsqvPHxERgWXLluHOnTuIi4tDXFwcUlJSkJiYCKBK5BcuXIiSkhKcOXMGd+/eRUZGBvLz83H//n3cv38feXl5SE9Px927d2s8D4fDQUBAANasWYNXXnkFcrkcI0aMQHp6OtRqNdLS0nDhwgUDG/Py8pCdnQ2pVIquXbs6NXSP0jChoWYUA/bv3w8AGD9+fKOcn0KpLWi0A4VBqVTi/v374HK5UCqVtd49oK7np1BqE7rypVAolDqA+nwpFAqlDqDiS6FQKHUAFV8KhUKpA6j4UigUSh1AxZdCoVDqACq+FAqFUgf8f8si6I9IS+dQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 373.798x231.02 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "utils.latexify(font_serif='Times New Roman', font_size=10, usetex=False)\n",
    "\n",
    "Z = results.to_numpy()\n",
    "k = results.columns.to_numpy()\n",
    "r = results.index.to_numpy()\n",
    "X, Y = np.meshgrid(k, r)\n",
    "\n",
    "nbins=20\n",
    "levels=mpl.ticker.MaxNLocator(nbins=nbins).tick_values(Z.min(),Z.max())\n",
    "\n",
    "fig_height, fig_aspect = get_fig_dim(width_pt, fraction=0.8)\n",
    "fig, ax = plt.subplots(figsize=(fig_height*fig_aspect, fig_height))\n",
    "\n",
    "cs = ax.contourf(X,Y,Z,levels=levels, cmap=sns.color_palette(\"crest\", as_cmap=True))\n",
    "fig.colorbar(cs,format=mpl.ticker.PercentFormatter(),ticks=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0])\n",
    "sns.despine(ax=ax)\n",
    "\n",
    "cit=23\n",
    "r_quer=0.74\n",
    "k_quer=0.33\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(k_quer + 0.01, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-150,angleB=90,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    ")\n",
    "\n",
    "cit=24\n",
    "r_quer=2.5\n",
    "k_quer=0.1\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(k_quer, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-150,angleB=90,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            )\n",
    "\n",
    "cit=25\n",
    "r_quer=3\n",
    "k_quer=0.5\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(k_quer + 0.01, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-150,angleB=90,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            )\n",
    "\n",
    "cit=45\n",
    "r_quer=2.47\n",
    "k_quer=0.2\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(k_quer + 0.01, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-150,angleB=90,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            )\n",
    "\n",
    "cit=46\n",
    "r_quer=0.61\n",
    "k_quer=0.11\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(0.15, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-90,angleB=0,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            )\n",
    "\n",
    "cit=47\n",
    "r_quer=0.67\n",
    "k_quer=0.25\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(k_quer + 0.01, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-150,angleB=90,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            )\n",
    "\n",
    "cit=48\n",
    "r_quer=0.4\n",
    "k_quer=0.3\n",
    "ax.annotate(r'[{cit}]: {perc}%'.format(cit=cit, perc=mapper(r_quer=r_quer, k_quer=k_quer, r=r, k=k, Z=Z)), xy=(k_quer, r_quer),  xycoords='data',\n",
    "            xytext=(0.41, 5.2), textcoords='data',\n",
    "            size=7,\n",
    "            arrowprops=dict(color='black', arrowstyle=\"-|>\", connectionstyle=\"angle, angleA=-90,angleB=0,rad=5\", linestyle='--', relpos=(0., 0.)),\n",
    "            bbox=dict(boxstyle=\"round\", fc=\"w\", ec=\"black\", linewidth=0.5),\n",
    "            rotation=60,\n",
    "            horizontalalignment='left', verticalalignment='bottom',\n",
    "            )\n",
    "\n",
    "ax.set_yticks([1,2,3,4,5])\n",
    "ax.set_ylabel(r'Reproductive number R')\n",
    "ax.set_xlabel(r'Dispersion parameter k')\n",
    "\n",
    "# ax.text(x=0.01, y=1, s=r'A', size=15, transform = ax.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "# ax.text(x=0.01, y=1, s=r'B', size=15, transform = ax.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "ax.text(x=0.01, y=1, s=r'C', size=15, transform = ax.transAxes, horizontalalignment='center', verticalalignment='bottom', weight='bold')\n",
    "\n",
    "fig.tight_layout()\n",
    "# fig.savefig('figures/fig2.tif', dpi=300)\n",
    "# fig.savefig('figures/figs3a.tif', dpi=300)\n",
    "# fig.savefig('figures/figs3b.tif', dpi=300)\n",
    "# fig.savefig('figures/figs4a.tif', dpi=300)\n",
    "# fig.savefig('figures/figs4b.tif', dpi=300)\n",
    "# fig.savefig('figures/figs4c.tif', dpi=300)\n",
    "# fig.savefig('figures/figs5a.tif', dpi=300)\n",
    "# fig.savefig('figures/figs5b.tif', dpi=300)\n",
    "fig.savefig('figures/figs5c.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
}