{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "_oaQGMHSBltD"
   },
   "source": [
    "# Network Epidemiology\n",
    "\n",
    "<div style=\"display: flex; align-items: center;\">\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/yy/netsci-course/blob/master/docs/m09-netepi/lab09.ipynb\">\n",
    "        <img src=\"https://www.tensorflow.org/images/colab_logo_32px.png\" style=\"vertical-align: -8px; margin-right: 10px;\" />\n",
    "        Open this notebook in Google Colab\n",
    "    </a>\n",
    "</div>\n",
    "\n",
    "\n",
    "<div style=\"display: flex; align-items: center;\">\n",
    "    <a href=\"https://raw.githubusercontent.com/yy/netsci-course/master/docs/m09-netepi/lab09.ipynb\" download>\n",
    "        <img src=\"http://yyahn.com/netsci-course/images/download_icon.png\" style=\"vertical-align: -8px; margin-right: 10px;\" />\n",
    "        Download this notebook (File -> Save As)\n",
    "    </a>\n",
    "</div>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "et3I3l-iBltF"
   },
   "source": [
    "You may remember the March of 2020. 👀\n",
    "\n",
    "<img src=\"https://yyiki.s3.us-east-2.amazonaws.com/public/imgs/month_of_2020-time_spent_looking_at_exponential_graphs.jpg\" width=\"400\">\n",
    "\n",
    "It was crazy how everything in our lives was suddenly governed by this new pathogen. It was a time when we all became armchair epidemiologists. We were all trying to understand the exponential growth of the virus and how it was spreading across the world. We were all trying to understand the basic reproduction number $R_0$ and how it was affecting the spread of the virus.\n",
    "\n",
    "Let's bring back some of that excitement(?) and learn how network perspectives can help us understand the spread of infectious diseases!\n",
    "\n",
    "## Friendship paradox, again. \n",
    "\n",
    "The first and probably the most important insight that network science can provide to epidemiology is the biased sampling happening in an epidemic. This is the exactly same bias that we saw in the friendship paradox.\n",
    "\n",
    "Imagine a contact network through which a disease is spreading. Imagine that someone in the network just caught the virus and will spread to their neighbors. In this case, who will be the most likely to be the next victim of the virus? If every contact is equally likely to be infected, then the probability is proportional to the likelihood of them _appearing_ in the contact list of the infected. And as we have learned multiple times, this is roughly proportional to the degree of the node due to the friendship paradox! \n",
    "\n",
    "Let's quickly simulate this. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "executionInfo": {
     "elapsed": 408,
     "status": "ok",
     "timestamp": 1648464750012,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "W2gd8rYsBltF"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "index:  254\n",
      "secondary:  [66, 38, 996]\n",
      "degree of index:  6\n",
      "degree of secondary:  [14, 42, 4]\n"
     ]
    }
   ],
   "source": [
    "import networkx as nx\n",
    "import random\n",
    "\n",
    "G = nx.barabasi_albert_graph(1000, 4)\n",
    "\n",
    "# sample a random node as the initial infected node (index case)\n",
    "index_node = random.choice(list(G.nodes()))\n",
    "secondary_infections = [\n",
    "    node for node in G.neighbors(index_node) if random.random() < 0.5\n",
    "]\n",
    "print(\"index: \", index_node)\n",
    "print(\"secondary: \", secondary_infections)\n",
    "print(\"degree of index: \", G.degree(index_node))\n",
    "print(\"degree of secondary: \", [G.degree(node) for node in secondary_infections])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's repeat many times to get the degree distribution (use CCDF). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG1CAYAAAAV2Js8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOOElEQVR4nO3deXxU1f3/8dfMJJM9gRDIRlgVNCwhBIhURVAUsKhAXeqKqFgtWi2VKq3VWq1WEcQlLYoF/Lb2V+peQQRFEYvIKigCIhAgLNnIvk4yM78/bgiEJJBAJneSeT8fj/uYufeeufMZwJm39557jsXtdrsRERER8UFWswsQERERMYuCkIiIiPgsBSERERHxWQpCIiIi4rMUhERERMRnKQiJiIiIz1IQEhEREZ+lICQiIiI+y8/sArydy+Xi8OHDhIWFYbFYzC5HREREmsDtdlNcXExcXBxWa+PnfRSETuPw4cMkJCSYXYaIiIicgYyMDLp27drofgWh0wgLCwOMP8jw8HCTqxEREZGmKCoqIiEhofZ3vDEKQqdx7HJYeHi4gpCIiEgbc7puLeosLSIiIj5LQUhERER8li6NiYhIu+VyuXA4HGaXIR7g7++PzWY76+MoCImISLvkcDhIT0/H5XKZXYp4SIcOHYiJiTmr4W0UhEREpN1xu90cOXIEm81GQkLCKceRkbbH7XZTVlZGdnY2ALGxsWd8LAUhERFpd6qrqykrKyMuLo7g4GCzyxEPCAoKAiA7O5suXbqc8WUyRWQREWl3nE4nAHa73eRKxJOOhdyqqqozPoaCkIiItFuaGql9a4m/X58IQkuWLKFv376ce+65vP7662aXIyIiIl6i3fcRqq6uZvr06Xz++edERESQkpLCxIkT6dSpk9mliYiIiMna/Rmh9evX069fP+Lj4wkNDWXcuHGsWLHC7LJERETqGTlyJA8++OAZv37fvn1YLBa2bNnSYjW1d14fhFavXs1VV11FXFwcFouF999/v16btLQ0evToQWBgIKmpqaxfv7523+HDh4mPj69dj4+P59ChQ61RuoiISLO8++67PPnkk2aX4VO8PgiVlpaSlJREWlpag/sXL17M9OnTefzxx9m8eTNJSUmMGTOmdmyB5qqsrKSoqKjO4gmbZk/g279cxuF9P3jk+CIi0vZERkaedrZ0aVleH4TGjRvHU089xcSJExvcP2fOHKZOncqUKVNITExk3rx5BAcHs2DBAgDi4uLqnAE6dOgQcXFxjb7fM888Q0RERO2SkJDQsh+oRo/izQys2EhlmWeCloiIHOd2uylzVJuyuN3uJtd54qWxHj168PTTT3PHHXcQFhZGt27deO211+q0X79+PcnJyQQGBjJkyBC++eabesfctm0b48aNIzQ0lOjoaG699VZyc3MBWLVqFXa7nS+//LK2/XPPPUeXLl3Iyso6gz/ptqdNd5Z2OBxs2rSJmTNn1m6zWq2MHj2atWvXAjBs2DC2bdvGoUOHiIiIYNmyZfzhD39o9JgzZ85k+vTptetFRUUeC0MiItI6yqucJD623JT33v6nMQTbz+zndvbs2Tz55JP87ne/4+233+bee+/lkksuoW/fvpSUlDB+/Hguv/xy/vnPf5Kens4DDzxQ5/UFBQVceuml3HXXXbzwwguUl5fz8MMPc/311/PZZ5/VBq9bb72VrVu3snfvXv7whz/w1ltvER0d3RIf3+u16SCUm5uL0+ms95cVHR3Nzp07AfDz82P27NmMGjUKl8vFb3/721PeMRYQEEBAQIBH6xYREWmKK6+8kl/+8pcAPPzww7zwwgt8/vnn9O3bl3/961+4XC7+/ve/ExgYSL9+/Th48CD33ntv7etfeeUVkpOTefrpp2u3LViwgISEBHbt2kWfPn146qmn+OSTT7j77rvZtm0bkydP5uqrr271z2qWNh2Emurqq6/2qb9UERGpK8jfxvY/jTHtvc/UwIEDa59bLBZiYmJq+8Du2LGDgQMHEhgYWNtm+PDhdV6/detWPv/8c0JDQ+sde8+ePfTp0we73c6bb77JwIED6d69Oy+88MIZ19sWtekgFBUVhc1mq3cdMysri5iYGJOqEhERb2OxWM748pSZ/P3966xbLBZcLleTX19SUsJVV13Fs88+W2/fiROVfvXVVwDk5eWRl5dHSEjIGVbc9nh9Z+lTsdvtpKSksHLlytptLpeLlStX1kvFzZWWlkZiYiJDhw492zJFRERa3Pnnn8+3335LRUVF7bavv/66TpvBgwfz/fff06NHD84555w6y7Gws2fPHn79618zf/58UlNTmTx5crPCVlvn9fG4pKSE3bt3166np6ezZcsWIiMj6datG9OnT2fy5MkMGTKEYcOGMXfuXEpLS5kyZcpZve+0adOYNm0aRUVFREREnO3HaFRQ1mYILAUsYLE085FT74fG94VEQXCkxz6XiIh41k033cTvf/97pk6dysyZM9m3bx/PP/98nTbTpk1j/vz53Hjjjfz2t78lMjKS3bt38+9//7t2yqlbbrmFMWPGMGXKFMaOHcuAAQOYPXs2M2bMMONjtTqvD0IbN25k1KhRtevH7uiaPHkyixYt4oYbbiAnJ4fHHnuMzMxMBg0axMcff+z1vd3dNUEl5ovfmlOA1R9+sRqiE815fxEROSuhoaF8+OGH3HPPPSQnJ5OYmMizzz7Lz372s9o2cXFxrFmzhocffpgrrriCyspKunfvztixY7FarTz55JPs37+fJUuWAMblstdee40bb7yRK664gqSkJLM+XquxuJszwIEPOnZGqLCwkPDw8BY77gt/+hVjqz+jV1QwATYr4Aa3+9SPcMI2mvaa2r/eE7ZVFIHbCdcuhP6TWuwziYh4i4qKCtLT0+nZs2edzsTSvpzq77mpv99ef0bILGlpaaSlpeF0Oj1y/DetV/GiYwzLrx9B35hWHkV00XjY9+Xp24mIiLRzbbqztCdNmzaN7du3s2HDBrNLEREREQ9REBIRERGfpSAkIiIiPkt9hHxZdSU4Grh132JtYJvF3FpFREQ8QEHIl71/j7E02amC0snbrGC1wYgZMPyXHilfRETkbOnSWCPa9cjSvUdRO+Bis7jB7QJXNbiqwOkAZyVUV0BVGVSVgqMEHMVQWQjlefDt4pauXkREpMXojFAjWmtkaVNc/BsYfr8xlpDbVXfMIbeLOuMP1dvmrv+ahrbt+x98+CuTPqCIiEjTKAj5Kj+7Z4+ft9ezxxcREa+watUqRo0aRX5+Ph06dDC7nGZTEDLZb9/5lhC7DavFUtsf2WKxYKGm681J62A5YTtYatatNRtPbn/iOie0j40IZNqocwj0t3n2A1ZXQNZ2o8/QyX2J6q0fe35yW+vx9ljAP9jzQU5ERHyCgpBJOocFklviYGtGgWk1JHXtwOhED8/JlrMT/ja8ZY8ZGAF3fwGRPVv2uCIiYgqHw4Hdbs7/4CoImeSNKUNZl56HG3C73UZXHGoea6YSc7vdx6cUO7YP6rbF2OByH29/4utpoP2b6/aTkVdORbVnpg8BoOsQiBsMhRkn9TNynTBPmqt+H6Q67U6YY+1EFYWQ+a2CkIi0S2+//TZPPPEEu3fvJjg4mOTkZD744ANCQkJ4/fXXmT17Nunp6fTo0YNf/epX/PKXx+/MPXjwIDNmzGD58uVUVlZy/vnnk5aWRmpqKgB/+9vfeP7558nIyKBnz548+uij3HrrrbWvt1gszJ8/n6VLl7J8+XLi4+OZPXs2V199dW2bjz76iAcffJCMjAwuuOACJk+eXKf+o0ePct9997F69Wry8/Pp3bs3v/vd77jxxhtr24wcOZL+/fvj5+fHP//5TwYMGEDPnj3Jzs6unQAWoKqqivj4eJ555hnuvPPOFv+zBgWhRnl6rrEu4YFclRTnkWOfzqofssnIK/fsmwR1hLs/b5ljuU8ISYt+ChlfQ1W5MXnsiZfRGrq0pvGPRASM75CqMnPe2z+4yd9FR44c4cYbb+S5555j4sSJFBcX8+WXX+J2u3nzzTd57LHHeOWVV0hOTuabb75h6tSphISEMHnyZEpKSrjkkkuIj4/nv//9LzExMWzevBmXywXAe++9xwMPPMDcuXMZPXo0S5YsYcqUKXTt2pVRo0bV1vDEE0/w3HPPMWvWLF5++WVuvvlm9u/fT2RkJBkZGUyaNIlp06Zx9913s3HjRn7zm9/U+QwVFRWkpKTw8MMPEx4eztKlS7n11lvp3bs3w4YNq233xhtvcO+997JmzRrACFAjRozgyJEjxMbGArBkyRLKysq44YYbzuqv4FQ0+/xpeGr2eTP9/LW1fL03j1duSmb8QHPC2BlbMA4OfNXMFzXW78gK/SfCNWmeqFRETFRvVnJHKTxt0vfd7w6DPaRJTTdv3kxKSgr79u2je/fudfadc845PPnkk3XOrDz11FN89NFHfPXVV7z22ms89NBD7Nu3j8jIyHrHvvDCC+nXrx+vvfZa7bbrr7+e0tJSli5dChhnhB599FGefPJJAEpLSwkNDWXZsmWMHTuW3/3ud3zwwQd8//33tcd45JFHePbZZ0/ZWXr8+PGcd955PP/884BxRqioqIjNmzfXadevXz8mT57Mb3/7WwCuvvpqOnXqxMKFCxs8bkvMPq9xhKRtOecymj8GktsYKsBVfcLYR+XGuEff/scTVYqInJGkpCQuu+wyBgwYwHXXXcf8+fPJz8+ntLSUPXv2cOeddxIaGlq7PPXUU+zZsweALVu2kJyc3GAIAtixYwcXXnhhnW0XXnghO3bsqLNt4MCBtc9DQkIIDw8nOzu79hjHLrMdM3x43X6gTqeTJ598kgEDBhAZGUloaCjLly/nwIEDddqlpKTUq/Guu+6qDT1ZWVksW7aMO+64o9E/r5agS2PStox4CC58oIH+RQ31Nzp5/wltijPh76NN/jAi0mr8g40zM2a9dxPZbDY++eQTvvrqK1asWMHLL7/M73//ez788EMA5s+fXy+I2GzG3b9BQUEtU66/f511i8VSe3mtKWbNmsWLL77I3LlzGTBgACEhITz44IM4HI467UJC6p8lu+2223jkkUdYu3YtX331FT179uTiiy8+sw/SRApC0vbY/E/f5nQsNSdD3S7I3V23n1GjSwNtbHaw6T8jEa9nsTT58pTZLBYLF154IRdeeCGPPfYY3bt3Z82aNcTFxbF3715uvvnmBl83cOBAXn/9dfLy8ho8K3T++eezZs2aOp2b16xZQ2JiYpNrO//88/nvf/9bZ9vXX39dZ33NmjVcc8013HLLLQC4XC527drVpPfp1KkTEyZMYOHChaxdu5YpU6Y0ubYzpW9w8W2uanil/unZJgsIh8n/hbjklqtJRHzWunXrWLlyJVdccQVdunRh3bp15OTkcP755/PEE0/wq1/9ioiICMaOHUtlZSUbN24kPz+f6dOnc+ONN/L0008zYcIEnnnmGWJjY/nmm2+Ii4tj+PDhzJgxg+uvv57k5GRGjx7Nhx9+yLvvvsunn37a5PruueceZs+ezYwZM7jrrrvYtGkTixYtqtPm3HPP5e233+arr76iY8eOzJkzh6ysrCYHrrvuuovx48fjdDrr3ZHmCeojJL4pLBbOuRwCO0BABNjDwB5aM1hjoHGmx+p3/MxRYyqLIGNDq5QsIu1feHg4q1ev5sorr6RPnz48+uijzJ49m3HjxnHXXXfx+uuvs3DhQgYMGMAll1zCokWL6NnTGErEbrezYsUKunTpwpVXXsmAAQP4y1/+UnvpbMKECbz44os8//zz9OvXj1dffZWFCxcycuTIJtfXrVs33nnnHd5//32SkpKYN28eTz/9dJ02jz76KIMHD2bMmDGMHDmSmJgYJkyY0OT3GD16NLGxsYwZM4a4OM93cNddY6ehu8YEON636MTlnbtg5xIYNwtS7za7QhE5wanuJhLvVlJSQnx8PAsXLmTSpEmnbNsSd43p0lgjPD2OkLQxFgtYbMAJU5JY9Z+PiEhLcblc5ObmMnv2bDp06FBnEEdP0jd5I9r17PPSsqpPHNyxsU7WGthRRORUDhw4QM+ePenatSuLFi3Cz691IoqCkMjZ+uQxYzmlBu44i0uG25eA1cMT34qItAE9evTAjN466iwtcqbOuQysTb2V/9igjlXHB3Q88BUUHDj9S0VExGN0RkjkTA2+DQbdXL8TdZ2lgU7Wbhe8MtS8eY9ERKSWgpDI2bCe1IG6qY7dlp/+BeTuqnvJzGqreW47Yb2ms3ZwJHTo1qIfQaQ9043R7VtL/P0qCPmwl1b+yL/XZ2CxgNViwVrzaDnhudVKzfqJ+0/R3nK8vb/NwjWD4kmMax/DDrQoS014+vCB5r/25rfh3Mtbth6RdubY2DkOh6PFpp4Q71NWZpxZP3lakOZQEPJBXcKMsRZ2ZZWwK6vEo++1JaOAxb8YfvqGvmbUTNj2bs2lMqfx6HKdtO6su6/sqNG3KOcHBSGR0/Dz8yM4OJicnBz8/f2xWtUltj1xu92UlZWRnZ1Nhw4daoPvmdCAiqfRHgdULKqo4qvduTicbtxuNy63G5cLXG43brfx6HIfWz/+3OXmeHv3Ce1d9dunHy1l6bdHGBAfwYf3X2T2R24f3pkK3/2nZsVSc8nMdvzRYsU4hXfCtoBQ+Ols6DnC1NJFzOBwOEhPT2/WhKHStnTo0IGYmBgsDQxRogEVz1J7HlAxPNCfsf1jPfoen/+QzdJvj3j0PXxO9+Hw3VuA21hc1UA1nOqfaDHGmScFIfFBdrudc889t96s59I++Pv7n9WZoGMUhBqhARXF6wy5AwZcB9WVxy+b1Xl01922cYGxoJO+4rusVqum2JBTUhASaUsCwoylKcJqzvoVZMCu5SfdlWar+9xqNSaa7ZKoAR5FxKcoCIm0V8du0d+z0liaIvkWuCbNczWJiHgZBSGR9qrfRDiwFsoLGrl85jr+3FEKZbmQs8vsqkVEWpWCkEh71ak33PJO09ruXAr/vglKc2DrYuPy2Ml3pVlPuJzW+TwIi/Fs/SIirUBBSESOz5mWnw7v3X369vZQ+M0Pxu35IiJtmIKQiECPiyDldig8aFwqc1XXHdSx9tEFWd+BowTK8xSERKTNUxASEbAHw1UvNq3tUzHGCNe7lkNoF7D61SzHLqHVrPvZISYJbPqaERHvpW8oEWkem78RhD566PRtB91s3IXWwKivIiLeQEFIPOpgfhkz3/22dmJWC8cnZj02easFsFqN9eP7j0/oauHEyV4hJiKIicnx2Kz6cTXF6D/C9veNy2Su6ppLZtU1l9Rqnuf+YLTd8qax1J41qjlzdOy5LQAumQGDbzPzE4mID1MQakR7nmKjNYQHGp1v88uq+H/rM1r8+HERgfzknKgWP640wdA7jeVU8vbCqyOhstBYd1XXTAnSgG/eVBASEdNo0tXTaI+TrrYGt9vN+1sOcSi/vGay1ppJWTk+cauxDdwcn7zVzfHJXI9N4OrGXXuMT7ZnkVtSyd9uHsy4AZ6dL03OUnWlMT7RsRBUu9ScNdrzOSyfCQmpcOcKs6sVkXZGk66KqSwWCxOTu7b4cfdkl5BbUtnixxUP8AswlsYc3WM8HtwAf+l+/JJZbYdrq/EY2Ruu/z/w13xRItLyFIRExBydzwO/QKiugIqCxtsd3Q1/jjbaN3R3WtS5MP4FzZEmImdEQUhEzBF1Djy0C0pza8YsOuGy2bGxi968Fipq+hnl7Gz4OAe+Mh67Dq3pgO1//OySLQB6XAj2kNb5TCLS5igIiYh5AiOMpTEP/QhHvgVXVf0701zVsPhmo93mN4ylIf0mwXULW752EWkXFISkTTpSWMGenBKsNbflWy0WrNbjzy0WsFksNfstWKzH14/dmm+1gJ/NavZHkVPxC4CEoY3vv3ExbP1/RihyVtXtlF2SDUd/hKJDrVeviLQ5CkLSJv1pyXZYcvbHuXJADH+9OeXsDyTm6DvWWBqyY4lxxqgk27hF3+pnjHJt9a95boeuKRDUsXVrFhGvoiAkbcqE5HjSj5ZS7XThdLlrb8t31T4ef96UgSE+3Z7t+aLFHDa78ZifDh/8suE20QPg3v+1Xk0i4nUUhKRNuSm1GzeldmtSW3dNGHLWGbfIjdPl5khhBVe8sNrD1Yqpel4Mw34BhRn1L51VlhiTx2Zvh39MqulgXXOWyOZv3M025A6IG2T2pxARD1MQknardhoP6k/FUVLZyCjH0n74B8GVzzW8rywPnj+3ZmDHlQ232fwG3PDm8XB07LFDdwiL9lzdItKqFIRExPcER8IvVkP2DuNMkdNh3JnmrIbv34WMdUa7Y3elncjqDw9sgYiWHzBURFqfgpCI+KbofsZyskE3wUcPQeEhIyA5HTWX1hyQl24EppV/go49j58p8guAkCg47yrws7f+ZxGRM6YgJCJyosBwmPRaw/teGwWHN8O3ixt/vX9ITT+jACMkhcXCtQugQ4Jn6hWRs6IgJD6tyuXi/v/3DTYLWK0WbBYLNqul7nOLBZuVBrZZiAoNYNLgeAL9Nb2DT7hqLmx7B6odx88WOavg238fb1NVClUnvKYwA+ZdCEPvMjphHzuDFJ8CCcNa+xOIyEkUhMQnhQT44W+zUOV08+HWw2d1rGC7jQnJ8S1UmXi12CRjOdmkV6G8AKrKoLqypt9RJbx9B+TuMqYJ+XJ2/df5h0BwJ+Nymi3AGNNozJ91t5pIK1IQEp8UHujPv+++gG2HinC63LW31TvdblwuN04XOF0unG7jee3+E9qu/jGHjLxyCsurTv+G0v4FdTCWE93yDmz5FzhKas4iVRqPW/9l7K8qhcLSuq/Z+m8FIZFWpCDUiLS0NNLS0nA6nWaXIh6S0j2SlO6RZ/z6aW9uJiOvvAUrknanQzcY+Uj97Ve/ZHS8PhaMqitg0yL47j/GkvG1cQnN6g9WqzGhbP+f1WyrGe8oPA4s9YeGEJHmURBqxLRp05g2bRpFRUVERJxiUkgRkeay+UPnPnW35e01QlDZUWM5Ufrq+pfW+oyFm07RaVtEmkRBSETEGyTfYtzOX55fMwp2ldG36KPfGpfcjo13VFFgtN/1MTzX+/hAj/5BcMEvIWWymZ9CpM1REBIR8QYWC8QPrr89+Za665Ul8OLAmjNHuXX3rX9NQUikmRSERM5SXqmDg/ll+Fmt2KwW/KwW/GyWOutWq/pySAsJCIUHv4OiI8dv4c9YD8tmgMsJLpfRr0hEmkRBSOQsvbjyR15c+eMp21gsGAHJasXPasFmszCwawcW3T5UIUmazx4CUeccXz92uSxnB/ypI1hsxlhFzir4yf1wzujjAzx27GEMGikigIKQyBkbNyCGdel5VFQ5qXa5qHa6qXa5G2zrdkOV003VCXchrt6Vw8H8crp1Cm6tkqW9iu4PoTFQkmmsu53GmEYA/5tjLCfqfZkxZpFfzdhFw6cZd6GJ+CAFIZEzNH5gHOMH1v3xcLvduNxQ7XLhdBnByAhINetOYwyicS9+SXmVhmaQFhISBdO31wzoWHO57MhWWPW0cbmsutJYCg8Y7fesrPt6qx9c/kTr1y3iBRSERFqQxWLBZgGb9dRTbuhqmLQ4qw0CwiCgZj08FvqOrdsmfx/s+9/xYLRrmXFr/sENsHGhMQWIfyD0GAEhnVr7E4iYQkFIRMRXdOxhLMc4So0gtH+NsRzTbTjc8XFrVydiCgUhERFfNfhWKM2G0lzjDFFJFhzaCNk7YG2acYbIL9CYHDbqXLOrFfEIBSEREV8VFgNXzjq+nvU9/O0nxl1oy393fHtgB5ixB2z6yZD2R/+qRUTE0CURrvgz5Ow05j+rLDH6EVUUGCNdKwhJO6R/1SImWvTVPqLC7NhtxvhC/n5W/K1W/P0s+Nus+Fmt2P2M8Yf8bVY6h9k5p0uY2WVLe2WxwE/uO75eWQLPxBvPP33C6IztH2RcKutxkTk1irQwBSEREwTZ/Sh1OFmwJr3Zr51/2xAuT4z2QFUiJ7HZjT5C1RWw7m919/WbCCGdwT8Yeo2E3qNMKVHkbCkIiZjghRuS+GR7FlVON9VOF1VOV82Ai8bzapcbR3Xd54fyyymurGb/0VKzyxdf4WeHm/4D+78yxiiqKocN84193793vN36+TAzw7iFX6SNURASMcHF53bm4nM7N+s1D/77G97fcthDFYk0otclxnJMymRjLKKqMijLg7WvQFUpuF2AgpC0PQpCIiLSdDEDjAWgPN8IQiJtmIKQiIicvT/HGJPB2sOMucuG/9LsikSaREFIRETOTEAExA2Gw5vBVQ0VhcayfKZx6cweasyDdt5PjbvNRLyQgpCIiJwZqxWmfgaVRcZ0HQe+hrenGPs+e/J4u1GPwiUzzKlR5DQUhETamPlf7uXDrYcJ8LNh97Ni97MSUO/R2BcR5M91KV3pFBpw+gOLnAmLBQIjjKXfROOMUM4P4CiBgxshZweU5phdpUijFIRE2ohunUIAyCqqJKuossmvyy91MPPK8z1VlshxFgsMmXJ8feWTRhDa8V/j0RYAfgE14xMFGGMQDZlyvPO1iAl8IghNnDiRVatWcdlll/H222+bXY7IGXnwsnO59LwulFRUU1ntxFHtorLaZTw6XVRWOXE4XVRWuXA4Xazbe5TNBwooqqg2u3TxVRE1o1IXHzGWhhQfgRv/X+vVJHISnwhCDzzwAHfccQdvvPGG2aWInDGr1cKghA5Nbv/Syh/ZfKDAY/WInNbg2yGqD5QdhWoHOB3grDSeH9oI370Fh7fA0oeMM0R+AcZI1jY7BEdCv0kQEGr2p5B2zieC0MiRI1m1apXZZYiI+BartfE5yX78xAhCxYePj1Z9stIcuPg3nqtPBC8IQqtXr2bWrFls2rSJI0eO8N577zFhwoQ6bdLS0pg1axaZmZkkJSXx8ssvM2zYMHMKFhGRs9drFEx8DQozoLqy5kxRzZKxDrK3w86lxrpfYM0SYNyG33UoRJ1r9ieQdsL0IFRaWkpSUhJ33HEHkyZNqrd/8eLFTJ8+nXnz5pGamsrcuXMZM2YMP/zwA126dAFg0KBBVFfX7wexYsUK4uLimlVPZWUllZXHO6IWFRU18xOJeJcDeaUs/z6TQH8bQf42Av2ttc8D/K0122z426xmlyq+xOYHSTc0vO+L54wgdGiTsZwsMAJm7DWOIXKWTP9XNG7cOMaNG9fo/jlz5jB16lSmTDHuRJg3bx5Lly5lwYIFPPLIIwBs2bKlxep55plneOKJJ1rseCJm8bNZAFiz+yhrdh89bXu7zcqMMX2ZOqKXp0sTObWhd4HFakzhUVVec6aowrglf9fHxi36zkoFIWkRXv2vyOFwsGnTJmbOnFm7zWq1Mnr0aNauXeuR95w5cybTp0+vXS8qKiIhIcEj7yXiSRMGxfNDZjG5JZWUO5xUVLmoqHJSUeWkvMpYL69y1rZ3OF0s23ZEQUjMFxwJIx6qv91RCk837yy/yOl4dRDKzc3F6XQSHR1dZ3t0dDQ7d+5s8nFGjx7N1q1bKS0tpWvXrrz11lsMHz68wbYBAQEEBGjwOWn74joE8eLPk0/Zxu12U1nt4uNtmTy4eEvrFCbSErb8CwI7QHQiRPczuxppw7w6CLWUTz/91OwSRLySxWIh0N9GsN1mdikip2exGYvbCR/VnDGyWGH6TgiLPvVrRRrh1UEoKioKm81GVlZWne1ZWVnExMSYVJWIiJjCPxDGvwB7Pzf6Du1eCa4q2Ph3CI83RrbG0sCj9aRtGI8WayPtT7GvUy+I1OXj9sSrg5DdbiclJYWVK1fW3lLvcrlYuXIl9913n0ffOy0tjbS0NJxO5+kbi7QTB/LKeOLD7wmx+xFktxFitxEc4EeI3Y9gu3HmKCTAj+6dggkL9De7XPFFKZONBeCF/sbt918823rvb7HBr7+H8NjWe0/xKNODUElJCbt3765dT09PZ8uWLURGRtKtWzemT5/O5MmTGTJkCMOGDWPu3LmUlpbW3kXmKdOmTWPatGkUFRURERHh0fcSMVuHYDsAuSUOFq7Zd9r2nULsrHnkUgL9dUlNTHTFk/Dd2+B2g9sFuI3ndR5dDWxzG69vdF8jr8vabpyBKj6iINSOmB6ENm7cyKhRo2rXj92xNXnyZBYtWsQNN9xATk4Ojz32GJmZmQwaNIiPP/64XgdqETlzQ7p35OUbkzmQV0ZpZTVlDidljmpKHU7KKmseHcb2vTmlHC11kFtSSdeOwWaXLr6s30RjaS3HzkBJu2Jxu49FY2nIsTNChYWFhIeHm12OiOnO+8MyKqpc/O/hUQpC4luOBaGYARAUaYx0bbMbd62NnHm8/5F4hab+fpt+RshbqY+QiIjUER5vBKHM7+pu37nEmCC2y3nm1CVnRUGoEeojJCIiddy0GA58bYxy7XQYI14v/z1U1ox0LW2SgpCInJGRs1YRGuhHWKAfoQH+hAX6ER7oR2iAH2GB/kQE+TNpcDy9OoeaXapIywjqAH3H1t32+dNGEJI2S0FIRJrlkj6dWf59FtUuNwVlVRSUVQHlDbbdfqSIBbcPbd0CRczww8eQ8wNYbWD1q7/YQyAu2dgvXkVBSESaZd4tKZQ6nJRUVFNcUUVxZTXFNc+NbdVsOVjA0m+PUFJRbXa5Ip5lqxlPa9XTp2970a9h9B89Wo40n4KQiDSLxWIhNMC4BBYTEdhgm2XfHWHpt0dIP1rK0x/tICLIn/BAP8KDjEtmx5aOwXY6hthb+ROItKDLn4Cti8FV3fhSkgMlmZCXbna10gAFoUborjGRM3dsgMac4kpeW733lG2njerNjDG620baqKaMZbR+/vG50cTraByh09A4QiLN53K5WbYtk/TcEgrLq05aqikqryKv1EF5lZOU7h15596fmF2yiOccC0JBkcaYQza7sfjVPIbHwSUPG/2IpMVoHCERMY3VauGnA089BcHy7zP5xT82cbignHlf7CEy2E6HYH8iQ+x0CLYTGWInIsgfm1WD1EkbFx5nPJbnwb4vG24Tl9y6o2RLLQUhETFFWIDx9XOksIK/LNvZYBuLBfpGh/HWPcM1yau0XX2vhDs/gZIsY+whZ5UxDpHTAetfg9xdUFVhdpU+S0FIREyR2qsTz0wawO7sEvJLHeSXOcgrq6KgzEFeqYPiimrcbtiZWcyOI8UM6xlpdskiZ8ZigYRhDe/b9bERhMQ0CkIiYgqb1cKNw7o1ur/K6WL0nC/Yf7SMwwXl5Jc6iAjyx6pLZSLSghSEGqG7xkTM5W+z4lcTeh5cvAUwwlNkiJ2o0ACiQu10CrHTPz6COy/qiUUTXorIGdBdY6ehu8ZEzPOvdQd4/X97OVrioLC8qtF2S391Ef3iNCegtEH//Bns/hT8gyEwwngM6gCXPQ69LjG7ujatqb/fCkKnoSAk4h0c1S7ySh3kllSSW1LJ0RIHTy3dTn6ZEZA6BvvTJSyQLuEBdA4NoHN4ADHhgfx0QCxdwhse+FHEdKuebXhU6v4/g2sXtH497YhunxeRdsXuZyUmIrDOaNZ7c0t4bfVeqpxu8suqyC+r4oes4jqvW7P7KK9PHtLa5Yo0zciHYcgUqCiCqlJjlOqv08ClbhmtRUFIRNqsGWPO46Er+lJQVkV2cSXZxRVkF1WSXVzJpv35fLoji/wyh9llipxaaBdjAchYb24tPkhBSETaNIvFQscQY86yvjFhtduXf5/JpzuyTKxMRNoCq9kFiIh4Uk5xJcu+O8KWjAKyiytwudQtUkSO0xkhEWmXgvxtABzIK+PeNzfXbrfbjL5GsRGBXDckgWtTuppVooh4AQWhRmgcIZG27YJenZgxpi/bjxRxpKCcwwUVZBdX4HC6OJBXxoG8Mtal5/H5D9l07RhEQsdgunYMomvNY2BNkBKR9k23z5+Gbp8XaT+qnC6yiirYlVXMHYs2NtrOYoGZ487j7hG9W7E6EY7PVG8PhfB4sPqB1QY2/5rnNetWP7D61133C4Rhd0F8itmfwivo9nkRkZP426w1Z3yC+XrmZWw9WEBGXhkH88trljIy8soodTj5bGe2gpC0vsiexqOjBHJ/aP7rS3PglrdbtqZ2TkFIRHySMSZRTL3tS749zH3/+obyKhdFFVWEa9Z7aU3njIb7NxuBxlkFrmpjTCFXdc1SVXf9WJtDm2Hrv6Bas9g3l4KQiMgJbDVzlm3NKGDgH1fQKcROj6gQenQKoWdUMCndIxneu5PJVUq71qm3sTTHtneMICTNpiAkInKC4b07cdl5Xdh6sNCYyqPUwdFSB5v259e2eXjseQzr2ZHIkAAiQ+yEB/pp0leRNkpBSETkBB2C7fz99qEAFFdUsf9oGem5pezLLWX2J7sAePbjnXVe42c1BnXsFGInsmZwx2PPO4cF0Cc6jPNiwgjTZTYRr6MgJCLSiLBAf/rHR9A/3pjZPjEunH9+vZ+8sirySivJK3FQ6nBS7XKTU1xJTnHlKY+XEBnEeTHhnB8bTmJsGOfFhNMtMhirVWeTRMyiICQi0kSXnR/NZedH19lWUeUkv8zB0RIH+WUO8kqPPz9a6uBIQTk7M4s5UlhBRl45GXnlfLL9+NQfIXYbfWPCOD/WCEgD4iM4LzaMAD+NYyTSGhSEGqEBFUWkKQL9bcRGBBEbEXTKdvmlDnZmFrPjSBE7jhSxM7OYH7KKKXU42XyggM0HCmrb2m1WzosNY2DXCAZ27cCghA707hyKTWeORFqcBlQ8DQ2oKCKeUu10kZ5byo6agPT94SK+O1hAfllVvbbBdhv94yNIqglHl/TtrFv75bht78Dbd0CPi+H2JWZX4xU0oKKIiJfzs1k5NzqMc6PDuDopDgC3283B/HK2Hizg24OFbM0o4LtDhZQ5nKxPz2N9eh4AA7tG8N/7LjKzfJF2QUFIRMSLWCwWEiKDSYgMZvxAIxw5XW725JSwNaOAr/Yc5b1vDpGRV2ZypSLtg4KQiIiXs1kt9IkOo090GMndOvDeN4fMLkmk3bCaXYCIiDRfqcPJ/63dx66sYtTVU+TM6YyQiEgb0jHYjs1qwVHt4rEPvgcgMsROas9I7rq4JyndI02uUKRt0RkhEZE2pFNoAB/96mIeuqIPF50TRaC/lbxSB8u2ZfKnJTvMLk+kzdEZIRGRNqZvTBh9Y8K471JwVLtYvDGDP7y/jQqHxj0TaS6dERIRacPsflZ6RYWYXYZIm6UgJCIiIj5LQagRaWlpJCYmMnToULNLERFpEje6e0ykuTTFxmloig0R8XZf7cnlpvnrAOgTHcromslhkxM6aGZ7X3Fsio2gSOj+E7DawOpXs/jXXbf5N3+/PQR6Xwr+gWZ/0ibTFBsiIj5iUEIHLj2vC1/symFXVgm7skr466o9dAqx0y8+gnO7hBpLdBjndAklIkhzlLU7IV2Mx/I82OmhucYu+jWM/qNnjm0inRE6DZ0REpG2oqDMwRe7cvh0RzarfsimuKK6wXbR4QGc28UIRYmx4Vw9KI5Af1srVystyu2GPZ9B0WFwVddfnA1sa+r+vHTI2wNJN8LEeWZ/0ibTGSERER/TIdjONYPiuWZQPFVOF98eLOTHrGJ+zC5hV1Yxu7NLOFJYQVZRJVlFlfxvdy4ARworeGD0uSZXL2fFYoFzLvPMsde8BJ/8wTPH9gIKQiIi7ZC/zUpK946kdO9YZ3txRRW7s0v4MbuEtzZmsGFfPvllDpOqFDGfgpCIiA8JC/QnuVtHkrt15MDRMjbsyze7JBFTKQiJiPgoS80NZe99cwirxcINQxPoGxNmblEirUzjCImI+Kgx/WKICQ+ksLyKBWvSGTN3NdekrWHxhgOUa7oO8RHNCkILFiygsrLSU7WIiEgr6h8fwf8eHsWC24cwpl80flYLWzMKePid70h9+lOeXLKd9NxSs8sU8ahmBaGpU6dSWFhYux4XF8e+fftauiYREWklfjYrl54Xzau3DmHtzMt4ZNx5JEQGUVRRzd//l86o51fx0FtbzS5TxGOaFYROHnKouLgYl8vVogWJiIg5OocFcM8lvfnioVEsvH0oI/t2BuDtTQepqNKlMmmf1EdIRETqsFotjDqvC2k3DTa7FBGPa1YQslgsWCyWRtdFRERE2pJm3T7vdrvp06dPbfgpKSkhOTkZq7VunsrLy2u5CkVEREQ8pFlBaOHChZ6qQ0RERKTVNSsITZ482VN1iIiIF/vr57uZOqIXYYGauV7alzMaWdrtdrNp0yb27duHxWKhZ8+eJCcnt6v+QmlpaaSlpeF06k4JEfFNwXYbw3t1Yu3eo7z02W7+7+v93HNJb24b3p1guyYmkPbB4j75nvjT+Pzzz7nzzjvZv39/7e30x8LQggULGDFihEcKNUtRUREREREUFhYSHh5udjkiIq3K5XLz0bYjzPlkF3tzjMEVe0aF8MmvR+Bn043HPuHY7PNJN8LEeWZX02RN/f1u1r/i3bt3M378eHr06MG7777Ljh072L59O2+99RZdu3blyiuvZO/evWddvIiIeAer1cL4gXGseHAEz1+XBEB6bikF5VUmVybSMpoVhObOncsFF1zAZ599xjXXXEPfvn0577zzmDRpEp9//jmpqam88MILnqpVRERM4mezcm1K19r1A3llJlYj0nKaFYRWrVrFgw8+2OA+i8XCgw8+yOeff94SdYmIiBeKjQgE4Lp5a/nTh9spqtCZIWnbmhWEDhw4wIABAxrd379/f/bv33/WRYmIiHd6657hjOkXjdPlZsGadC59fhVvfLWP3BJNyC1tU7OCUElJCcHBwY3uDw4OpqxMp0tFRNqrrh2DefXWIfzfHcPo1TmE3BIHj//3e4b9+VNufO1r/rF2H9lFFWaXKdJkzb7/cfv27WRmZja4Lzc396wLEhER7zeiT2c+fmAEb67bz7ubD/HdoULW7j3K2r1Heey/3zO0eyTjBsQwtn8MsRFBZpcr0qhmB6HLLrus3iz0YPQRcrvd7WosIRERaZzdz8qUC3sy5cKeZOSV8fG2TD7adoRvDhSwfl8e6/fl8cSH24nvEERiXDiJseH0iwsnMS6c+A5B+r0Qr9CsIJSenu6pOkREpA1LiAxm6oheTB3Ri8MF5Xy8LZNl246wcX8+hwrKOVRQzifbs2rbhwf6kRgXTr+4CBJjjXB0TpdQ/DU2kfcqPAi7loNfAPgFNvB4wnOrzexqm6zZAyr6Gg2oKCJy5grLq9hxpIjth4vYfqSI7w8X8WNWMdWu+j89dpuVPjGhJMaGM7BrB342uCtB9rbzg9puff03+PiR5r3G6t9IWKp59D8hOHW/EIZNbfGym/r73awg9OOPP/LYY4/x6quv1jtoYWEh9957L0899RS9evU688q9jIKQiEjLqqx2sju7pE442nG4iOLK6jrtZozpy7RR55hUpdQqPATLZ0JJDlRXQHXlCY/lx9dd1ac/VoMs8Mh+CIxo0bKb+vvdrEtjs2bNIiEhocEDRkREkJCQwKxZs/jb3/7W/IpFRMQnBPjZ6BcXQb+44z98brebg/nlfH+4iIVr0lmXnkdBmcPEKqVWRDxc/3+nb+esBmflCUHp5NB00mNVOXz4K8ANTvPGo2pWEPriiy/45z//2ej+66+/nptuuumsixIREd9isVhIiAwmITKYbw7ksy49z+ySpLlsfsZiD2n6az78lefqaaJmD6jYpUuXRvdHRUWRkZFx1kWJiIiItIZmBaGIiAj27NnT6P7du3erH42IiIi0Gc0KQiNGjODll19udP9LL73ExRdffNZFiYiIVFa7zC5BfECzgtDMmTNZtmwZ1157LevXr6ewsJDCwkLWrVvHz372M5YvX87MmTM9VauIiPiAQH/jlvn/W7ufu97YyLZDhSZXJO1ZszpLJycn8/bbb3PHHXfw3nvv1dnXqVMn/vOf/zB48OAWLVBERHzL7T/pwYG8Mj7YcohPd2Tx6Y4srkiM5oHR59a500ykJZzRgIrl5eV8/PHH7N69G7fbTZ8+fbjiiitOOSFrW6VxhEREzLEnp4SXV/7IB1sPc+yXaky/aH5zRV/6RIeZW5y0jD/WBNsZeyAkqkUP7ZEBFT/77DPuu+8+vv766wYHVPzJT37CvHnz2lU/IQUhERFz7c4u4eXPfuS/NYEoLNCPb/5wOX6ajqPt84Ig1Kx/RXPnzmXq1KmNDqj4i1/8gjlz5jS/WhERkUac0yWUF3+ezJL7LwKguKIah1MdqaVlNCsIbd26lbFjxza6/4orrmDTpk1nXVRLysjIYOTIkSQmJjJw4EDeeusts0sSEZEz0DOqGQP1iTRRszpLZ2Vl4e/v3/jB/PzIyck566Jakp+fH3PnzmXQoEFkZmaSkpLClVdeSUiI/oMSEWmrHv/ge+I7BhEbEUhMRBAx4YHERAQSHuiHxWIxuzxpQ5oVhOLj49m2bRvnnNPwJHjffvstsbGxLVJYS4mNja2tKSYmhqioKPLy8hSERETaGLvNSnigH0UV1by16WCDbYL8bcRGBBIdHmg8RgTWWY8JDyQqNACrVWFJDM0KQldeeSV/+MMfGDt2LIGBgXX2lZeX8/jjjzN+/PhmFbB69WpmzZrFpk2bOHLkCO+99x4TJkyo0yYtLY1Zs2aRmZlJUlISL7/8MsOGDWvW+wBs2rQJp9NJQkJCs18rIiLm8rNZ+e99F/H13qNkFlWQWVhR57GgrIryKid7c0vZm1va+HGsFqLDA4kODyA2IqheaIoJD6RLeAABfrZW/HRilmYFoUcffZR3332XPn36cN9999G3b18Adu7cSVpaGk6nk9///vfNKqC0tJSkpCTuuOMOJk2aVG//4sWLmT59OvPmzSM1NZW5c+cyZswYfvjhh9p5zwYNGkR1dXW9165YsYK4uDgA8vLyuO2225g/f36z6hMREe/RIyqEHo30FSp3OMkqquBIYcVJj+VkFlWSWVhOTnEl1S43hwrKOVRQDhQ0+l5RofbjIanOYxAxEQHERAQRGtCsn1HxQs0eR2j//v3ce++9LF++nGMvtVgsjBkzhrS0NHr27HnmxVgs9c4IpaamMnToUF555RUAXC4XCQkJ3H///TzyyCNNOm5lZSWXX345U6dO5dZbbz1t28rKytr1oqIiEhISdPu8iEg7UO10kVNSaZxFOumM0okBytHE6T1CA/wY3L0jr982BLufbudvNi+4fb7ZUbZ79+589NFH5Ofn1w6oeO6559KxY8ezKrghDoeDTZs21Zm2w2q1Mnr0aNauXdukY7jdbm6//XYuvfTS04YggGeeeYYnnnjijGsWERHv5WezEhsRRGxEUKNt3G43BWVVdYKREZiOn1k6UlhBcUU1JZXVrN6Vw4/ZxRr1uo0643N6HTt2ZOjQoS1ZSz25ubk4nU6io6PrbI+Ojmbnzp1NOsaaNWtYvHgxAwcO5P333wfgH//4BwMGDGiw/cyZM5k+fXrt+rEzQiIi4hssFgsdQ+x0DLGTGNf4mYQyRzWXzFpFTnElzZ+jQbxFu7+4edFFF+FyNX3grYCAAAICAjxYkYiItAfBdj9sulW/zfPqC5pRUVHYbDaysrLqbM/KyiImJsakqkREROqqbGKfIvE+Xh2E7HY7KSkprFy5snaby+Vi5cqVDB8+3KPvnZaWRmJioscv/4mISNsV6G/8jE5esJ45n+yisLzK5IqkuUwPQiUlJWzZsoUtW7YAkJ6ezpYtWzhw4AAA06dPZ/78+bzxxhvs2LGDe++9l9LSUqZMmeLRuqZNm8b27dvZsGGDR99HRETarr/8bCDnxYRRUlnNSyt/5KJnP+PFT3+kqEKBqK1o9u3zLW3VqlWMGjWq3vbJkyezaNEiAF555ZXaARUHDRrESy+9RGpqaqvUp9nnRUTkVFwuN8u/z+SFT3exK6sEgIggf+4e0YvJP+mhsYZOxQtunzc9CHk7BSEREWkKl8vN0u+OMPfTXezJMUa27hjszy8u6c1tw7sTbFcgqscLgpDpl8ZERETaA6vVwlVJcaz49SW8+PNB9IwKIb+sir8s28nFz37O/37MNbtEaYCCUCPUWVpERM6EzWrhmkHxfPLrEcy+LomEyCCOljr494YDZpcmDVAQaoQ6S4uIyNnws1n5WUpXpl7cC0CDLnopBSERERHxWQpCIiIi4rMUhERERMRnKQiJiIh4kM1qzEe27XAh2cUVJlcjJ1MQaoTuGhMRkZYw+vxoOocFsP9oGdfPW8vB/DKzS5ITKAg1QneNiYhIS4gOD+StXwwnvkMQ+46Wcd28tezOLjG7LKmhICQiIuJhPaJCePve4fTuHMKRwgpueHUt2w4Vml2WoCAkIiLSKmIjgvjPL4bTPz6co6UObpz/NRv35Zldls9TEBIREWklnUID+NfUCxjWI5Liimpu+fs6vtiVY3ZZPk1BSEREpBWFB/rzxh3DGNm3MxVVLu56YwOb9uvMkFkUhBqhu8ZERMRTguw2Xrt1CJcnRlPldPPYB9/jdGkODjMoCDVCd42JiIgn2f2s/GXSAMIC/fj+cBFvbcwwuySfpCAkIiJikk6hATxw2bkAPL/iB4oqqkyuyPcoCImIiJjotuE96NU5hNwSB698ttvscnyOgpCIiIiJ7H5W/jA+EYCFa9LZm6PBFluTgpCIiIjJRvXtwqi+nalyuvnz0h1ml+NTFIRERES8wKPjE/GzWli5M5tXPvuRH7OKcbt1J5mn+ZldgLdKS0sjLS0Np9NpdikiIuIDencO5faf9OD1/6Xz/IpdPL9iF51C7AzrGUlqz0hSe3Wib3QY1prZ7KVlWNyKm6dUVFREREQEhYWFhIeHm12OiIi0Y45qFwvXpPPFrhw2H8inospVZ3+HYH+G9jCC0QW9OnF+bDi2thyM/hhhPM7YAyFRLXropv5+KwidhoKQiIiYwVHt4tuDBaxLz+PrvUfZtD+fMkfdqxRhgX61wSi1Vyf6x4XjZ2tDvV4UhLyfgpCIiHiDKqeLbYcKWZeex7q9R9m4L5/iyuo6bULsNlJ6RHJBr0hSe3ZiYNcI/L05GHlBEFIfIRERkTbA32YluVtHkrt15J5LelPtdLHjSDHr0o/y9d481qcfpaiimtW7clhdM5FrkL+NlO4da88YJSVEEOBnM/mTeBcFIRERkTbIz2ZlQNcIBnSN4K6Le+F0udmZWcS6vXmsSz/K+vQ88suq+N/uXP63OxeAAD8ryd06kNqzE6m9IhncrSOB/r4djBSERERE2gGb1UK/uAj6xUVwx0U9cbnc/Jhdwrr0o7XhKLfEwdd78/h6bx6sBLvNSlJCBKk9O3FBr04M7t6BYLtvRQP1EToN9RESEZH2wO12syentE4wyiqqrNPGz2phYNcIUnt1IrVnJEN6RBIa4MFg5AV9hBSETkNBSERE2iO3283+o2UnBKM8DhWU12ljs1roHxfOBb06cefFPekSFtiyRXhBEPKt818iIiICgMVioUdUCD2iQrhhaDcAMvLKau9KW5eex4G8MrYeLGTrwULKq5z86Zr+Jlfd8hSEGqGRpUVExNckRAaTEBnMtSldAThSWM6cFbt4a9NBSk66Vb+98OLBBcw1bdo0tm/fzoYNG8wuRURExBSxEUGcGx1qdhkepSAkIiIiPktBSERERHyWgpCIiIj4LAUhERER8VkKQiIiIuKzFIRERETEZykIiYiIiM9SEBIRERGfpSAkIiIiPktBqBFpaWkkJiYydOhQs0sRERERD1EQaoSm2BAREWn/FIRERETEZykIiYiIiM9SEBIRERGfpSAkIiIiPktBSERERHyWgpCIiIj4LAUhERER8VkKQiIiIuKzFIRERETEZykIiYiIiM9SEBIRERGfpSAkIiIiPktBSERERHyWgpCIiIj4LAUhERER8VkKQo1IS0sjMTGRoUOHml2KiIiIeIiCUCOmTZvG9u3b2bBhg9mliIiIiIcoCImIiIjPUhASERERn6UgJCIiIj5LQUhERER8loKQiIiI+CwFIREREfFZCkIiIiLisxSERERExGcpCImIiIjPUhASERERn6UgJCIiIj5LQUhERER8loKQiIiI+CwFIREREfFZCkIiIiLisxSERERExGcpCImIiIjPUhASERERn6UgJCIiIj5LQUhERER8VrsPQgUFBQwZMoRBgwbRv39/5s+fb3ZJIiIi4iX8zC7A08LCwli9ejXBwcGUlpbSv39/Jk2aRKdOncwuTUREREzW7s8I2Ww2goODAaisrMTtduN2u02uSkRERLyB6UFo9erVXHXVVcTFxWGxWHj//ffrtUlLS6NHjx4EBgaSmprK+vXrm/UeBQUFJCUl0bVrV2bMmEFUVFQLVS8iIiJtmelBqLS0lKSkJNLS0hrcv3jxYqZPn87jjz/O5s2bSUpKYsyYMWRnZ9e2Odb/5+Tl8OHDAHTo0IGtW7eSnp7Ov/71L7Kyslrls4mIiIh3M72P0Lhx4xg3blyj++fMmcPUqVOZMmUKAPPmzWPp0qUsWLCARx55BIAtW7Y06b2io6NJSkriyy+/5Nprr22wTWVlJZWVlbXrRUVFTfwkIiIi0taYfkboVBwOB5s2bWL06NG126xWK6NHj2bt2rVNOkZWVhbFxcUAFBYWsnr1avr27dto+2eeeYaIiIjaJSEh4ew+hIiIiHgt088InUpubi5Op5Po6Og626Ojo9m5c2eTjrF//37uvvvu2k7S999/PwMGDGi0/cyZM5k+fXrtelFRkcKQiIj4rBuGdOPyxBhCAmxml+IRXh2EWsKwYcOafOkMICAggICAAM8VJCIi0oZEBPsTEexvdhke49WXxqKiorDZbPU6N2dlZRETE2NSVSIiItJeeHUQstvtpKSksHLlytptLpeLlStXMnz4cI++d1paGomJiQwdOtSj7yMiIiLmMf3SWElJCbt3765dT09PZ8uWLURGRtKtWzemT5/O5MmTGTJkCMOGDWPu3LmUlpbW3kXmKdOmTWPatGkUFRURERHh0fcSERERc5gehDZu3MioUaNq1491VJ48eTKLFi3ihhtuICcnh8cee4zMzEwGDRrExx9/XK8DtYiIiEhzWdyab+KUjp0RKiwsJDw83OxyRERE2o8/1lxxmbEHQlp21oem/n57dR8hEREREU9SEGqEOkuLiIi0fwpCjZg2bRrbt29nw4YNZpciIiIiHqIgJCIiIj5LQUhERER8loKQiIiI+CwFIREREfFZCkKN0F1jIiIi7Z+CUCN015iIiEj7pyAkIiIiPktBSERERHyWgpCIiIj4LAUhERER8VkKQo3QXWMiIiLtn4JQI3TXmIiISPunICQiIiI+S0FIREREfJaCkIiIiPgsBSERERHxWQpCIiIi4rMUhBqh2+dFRETaPwWhRuj2eRERkfZPQUhERER8loKQiIiI+CwFIREREfFZCkIiIiLisxSERERExGcpCImIiIjPUhASERERn6Ug1AgNqCgiItL+KQg1QgMqioiItH8KQiIiIuKzFIRERETEZykIiYiIiM9SEBIRERGfpSAkIiIiPktBSERERHyWgpCIiIj4LAUhERER8VkKQiIiIuKzFIRERETEZykINUJzjYmIiLR/CkKN0FxjIiIi7Z+CkIiIiPgsBSERERHxWQpCIiIi4rMUhERERMRn+ZldgIiIiPioia8aj/ZQ00pQEBIRERFzJP3c7Ap0aUxERER8l4KQiIiI+CwFIREREfFZCkIiIiLisxSERERExGcpCImIiIjPUhASERERn6UgJCIiIj5LQUhERER8loJQI9LS0khMTGTo0KFmlyIiIiIeYnG73W6zi/BmRUVFREREUFhYSHh4uNnliIiISBM09fdbZ4RERETEZykIiYiIiM/S7POncezKYVFRkcmViIiISFMd+90+XQ8gBaHTKC4uBiAhIcHkSkRERKS5iouLiYiIaHS/Okufhsvl4vDhw4SFhWGxWFr02EOHDmXDhg0tekxve29PvE9LHvNsj3Wmry8qKiIhIYGMjAx1wjeRmf8NthZv/4zt/XvQU+/RUsdtieN46/eg2+2muLiYuLg4rNbGewLpjNBpWK1Wunbt6pFj22w2034EW+u9PfE+LXnMsz3W2b4+PDxcQchEZv432Fq8/TO29+9BT71HSx23JY7jzd+DpzoTdIw6S5to2rRp7f69PfE+LXnMsz2WmX+HcvZ84e/P2z9je/8e9NR7tNRxW+I43v5v7HR0aUzEBBqfSkR8nbd8D+qMkIgJAgICePzxxwkICDC7FBERU3jL96DOCImIiIjP0hkhERER8VkKQiIiIuKzFIRERETEZykIiYiIiM9SEBIRERGfpSAk4oUmTpxIx44dufbaa80uRUSkVWVkZDBy5EgSExMZOHAgb731lkffT7fPi3ihVatWUVxczBtvvMHbb79tdjkiIq3myJEjZGVlMWjQIDIzM0lJSWHXrl2EhIR45P10RkjEC40cOZKwsDCzyxARaXWxsbEMGjQIgJiYGKKiosjLy/PY+ykIibSw1atXc9VVVxEXF4fFYuH999+v1yYtLY0ePXoQGBhIamoq69evb/1CRUQ8oCW/Azdt2oTT6SQhIcFj9SoIibSw0tJSkpKSSEtLa3D/4sWLmT59Oo8//jibN28mKSmJMWPGkJ2d3cqVioi0vJb6DszLy+O2227jtdde82i96iMk4kEWi4X33nuPCRMm1G5LTU1l6NChvPLKKwC4XC4SEhK4//77eeSRR2rbrVq1ildeeUV9hESkzTrT78DKykouv/xypk6dyq233urRGnVGSKQVORwONm3axOjRo2u3Wa1WRo8ezdq1a02sTETE85ryHeh2u7n99tu59NJLPR6CQEFIpFXl5ubidDqJjo6usz06OprMzMza9dGjR3Pdddfx0Ucf0bVrV4UkEWkXmvIduGbNGhYvXsz777/PoEGDGDRoEN99953HavLz2JFF5Ix9+umnZpcgImKKiy66CJfL1WrvpzNCIq0oKioKm81GVlZWne1ZWVnExMSYVJWISOvwxu9ABSGRVmS320lJSWHlypW121wuFytXrmT48OEmViYi4nne+B2oS2MiLaykpITdu3fXrqenp7NlyxYiIyPp1q0b06dPZ/LkyQwZMoRhw4Yxd+5cSktLmTJliolVi4i0jLb2Hajb50Va2KpVqxg1alS97ZMnT2bRokUAvPLKK8yaNYvMzEwGDRrESy+9RGpqaitXKiLS8trad6CCkIiIiPgs9RESERERn6UgJCIiIj5LQUhERER8loKQiIiI+CwFIREREfFZCkIiIiLisxSERERExGcpCImIiIjPUhASkTZp5MiRPPjgg2aXISJtnIKQiIiI+CwFIRGRBjgcDrNLEJFWoCAkIl6vtLSU2267jdDQUGJjY5k9e3ad/ZWVlTz00EPEx8cTEhJCamoqq1atqtNm/vz5JCQkEBwczMSJE5kzZw4dOnSo3f/HP/6RQYMG8frrr9OzZ08CAwMBKCgo4K677qJz586Eh4dz6aWXsnXr1jrH/uCDDxg8eDCBgYH06tWLJ554gurqao/8WYhIy1IQEhGvN2PGDL744gs++OADVqxYwapVq9i8eXPt/vvuu4+1a9fy73//m2+//ZbrrruOsWPH8uOPPwKwZs0a7rnnHh544AG2bNnC5Zdfzp///Od677N7927eeecd3n33XbZs2QLAddddR3Z2NsuWLWPTpk0MHjyYyy67jLy8PAC+/PJLbrvtNh544AG2b9/Oq6++yqJFixo8voh4IbeIiBcrLi522+1293/+85/abUePHnUHBQW5H3jgAff+/fvdNpvNfejQoTqvu+yyy9wzZ850u91u9w033OD+6U9/Wmf/zTff7I6IiKhdf/zxx93+/v7u7Ozs2m1ffvmlOzw83F1RUVHntb1793a/+uqrte/z9NNP19n/j3/8wx0bG3vmH1pEWo2f2UFMRORU9uzZg8PhIDU1tXZbZGQkffv2BeC7777D6XTSp0+fOq+rrKykU6dOAPzwww9MnDixzv5hw4axZMmSOtu6d+9O586da9e3bt1KSUlJ7XGOKS8vZ8+ePbVt1qxZU+cMkNPppKKigrKyMoKDg8/0o4tIK1AQEpE2raSkBJvNxqZNm7DZbHX2hYaGNutYISEh9Y4dGxtbr78RUNu/qKSkhCeeeIJJkybVa3Osn5GIeC8FIRHxar1798bf359169bRrVs3APLz89m1axeXXHIJycnJOJ1OsrOzufjiixs8Rt++fdmwYUOdbSevN2Tw4MFkZmbi5+dHjx49Gm3zww8/cM455zTvg4mIV1AQEhGvFhoayp133smMGTPo1KkTXbp04fe//z1Wq3GvR58+fbj55pu57bbbmD17NsnJyeTk5LBy5UoGDhzIT3/6U+6//35GjBjBnDlzuOqqq/jss89YtmwZFovllO89evRohg8fzoQJE3juuefo06cPhw8fZunSpUycOJEhQ4bw2GOPMX78eLp168a1116L1Wpl69atbNu2jaeeeqo1/ohE5CzorjER8XqzZs3i4osv5qqrrmL06NFcdNFFpKSk1O5fuHAht912G7/5zW/o27cvEyZMYMOGDbVnkC688ELmzZvHnDlzSEpK4uOPP+bXv/71aS9dWSwWPvroI0aMGMGUKVPo06cPP//5z9m/fz/R0dEAjBkzhiVLlrBixQqGDh3KBRdcwAsvvED37t099wciIi3G4na73WYXISLS2qZOncrOnTv58ssvzS5FREykS2Mi4hOef/55Lr/8ckJCQli2bBlvvPEGf/3rX80uS0RMpjNCIuITrr/+elatWkVxcTG9evXi/vvv55577jG7LBExmYKQiIiI+Cx1lhYRERGfpSAkIiIiPktBSERERHyWgpCIiIj4LAUhERER8VkKQiIiIuKzFIRERETEZykIiYiIiM9SEBIRERGf9f8BrDNISWTMHZAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# repeat the simple simulation to obtain the degree distribution of the index and secondary nodes.\n",
    "# The degree distribution of the secondary nodes is expected to be different from the index node.\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "N = 1000\n",
    "infection_probability = 0.5\n",
    "index_degrees = []\n",
    "secondary_degrees = []\n",
    "\n",
    "# YOUR SOLUTION HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we expected, the secondary infection is more likely to happen to the nodes with higher degree. If you are well-connected, then you're more likely to get infected! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating functions, contact tracing, and friendship paradox\n",
    "\n",
    "One fascinating implication of this happens in the context of contact tracing, which I published during the pandemic: [Kojaku et al., The effectiveness of backward contact tracing in networks, Nature Physics, 2021](https://www.nature.com/articles/s41567-021-01187-2). \n",
    "\n",
    "Let's first write down the generating function for the degree distribution of the underlying network. With the degree distribution represented by $\\{ p_k \\}$, the probability generating function is written as:\n",
    "\n",
    "$$ G_0(x) = \\sum_k p_k x^k. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, imagine a node with degree $k$ that has just infected by another node. If we think about the remaining degree that the infected node can further spread the disease, that's $k-1$ (because one of the edges is already used to infect the current node). Then the generating function for the remaining degree (that the node can propagate the disease) distribution is:\n",
    "\n",
    "$$ G_1(x) = \\frac{1}{\\langle k \\rangle} \\sum_k k p_k x^{k-1} = \\frac{G_0'(x)}{G_0'(1)}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that, if we randomly sample those who are infected, then we expect to see this distribution. \n",
    "\n",
    "Interestingly, when we do contact tracing, another factor comes into play. First, think about two directions of contact tracing. When we found an infected person, we can try to identify who infected this person (backward contact tracing) or whom this person infected (forward contact tracing).\n",
    "\n",
    "Forward tracing will produce the same distribution because we are simply following the edge to sample a new node, which is, according to the friendship paradox principle, the same as $G_1(x)$.\n",
    "\n",
    "By contrast, backward tracing will produce a different distribution. This is because the more \"offspring\" a node has created, the more likely it is to be sampled! This is another bias that comes into play. \n",
    "\n",
    "<img src=\"https://yyiki.s3.us-east-2.amazonaws.com/public/imgs/backward_contact_tracing_illustration-fig1.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from this figure (d), the \"parent\" with more \"offsprings\" is more likely to be sampled and the factor is their degree in the transmission tree. If you work out the generating function for this distribution, it simply turns out to be:\n",
    "\n",
    "$$ G_2(x) = \\frac{G'_1(x)}{G'_1(1)} = \\frac{1}{\\sum_k k(k-1)p_k} \\sum_k k(k-1)p_k x^{k-2}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this means that the degree distribution sampled from backward contact tracing process is _even more_ biased towards the nodes with higher degree ($\\sim k^2 p_k$)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an **extra credit** qusetion (10pts): **Can you show these three degree distributions by running simple simulations?** It is essentially Fig. 2a in the paper: https://www.nature.com/articles/s41567-021-01187-2/figures/2 \n",
    "\n",
    "The first degree distribution is easy. It is just the degree distribution of the network. The second one is what we did above (randomly sampling infected nodes). You just need to plot the correct, expected degree distribution. The last one is a bit tricky. You'll need to simulate the transmission at least for a few steps and keep track of \"parents\" of the infected nodes so that we can trace back to them. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "7U8uddN0LPf3"
   },
   "source": [
    "## SI model\n",
    "\n",
    "Let's do a quick and dirty implementation of the simplest model, the SI model. In this model, every node is in one of the two states: susceptible (S) or infected (I). When a susceptible node is infected, it becomes infected and stays infected forever. The infected nodes can infect their neighbors with a certain probability.\n",
    "\n",
    "Given a network $G$, we can start with initializing the state of the nodes. Note that we are using a list of states to keep track of how the state changes over time for every node. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 263,
     "status": "ok",
     "timestamp": 1648465170308,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "xzPd867cLPRQ",
    "outputId": "9adbd7a0-343e-4475-c631-87a52d095c55"
   },
   "outputs": [],
   "source": [
    "def initalize_graph(G):\n",
    "    for node in G.nodes():\n",
    "        G.nodes[node][\"state\"] = [\"S\"]\n",
    "    return G"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "YPA1eSbCLmGg"
   },
   "source": [
    "Now we can randomly choose a certain fraction of nodes and mark them as initially infected using `choice()` function in numpy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 8,
     "status": "ok",
     "timestamp": 1648465374102,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "aIMIziOZL-8w",
    "outputId": "dc2b6bf8-30b8-48dc-8c1f-e7b8dbc32fd3"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{2632, 3187, 5233, 6914, 7019, 7681, 8675, 8924, 9159, 9363}"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init_infected_fraction = 0.001\n",
    "init_infected = set(\n",
    "    np.random.choice(\n",
    "        G.nodes(), size=int(len(G) * init_infected_fraction), replace=False\n",
    "    )\n",
    ")\n",
    "init_infected"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "Da9t1uoCMgEq"
   },
   "source": [
    "We should set these nodes' states to be `I` and then we can begin the simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 510,
     "status": "ok",
     "timestamp": 1648465499752,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "RgO-hI6AMpUq",
    "outputId": "afc43957-b2c4-4fc8-d03b-7444e085c9ed"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['S']\n",
      "['I']\n"
     ]
    }
   ],
   "source": [
    "G = initalize_graph(G)\n",
    "for node in init_infected:\n",
    "    G.nodes[node][\"state\"][0] = \"I\"\n",
    "print(G.nodes[1455][\"state\"])\n",
    "print(G.nodes[2632][\"state\"])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "Cuh5w7hPNULr"
   },
   "source": [
    "But how should we update the states and run the simulation? \n",
    "\n",
    "There is an important issue here. If we immediately change the state of a particular node, that may affect the updating of other nodes! Imagine an extreme case where we are simulating a highly infectious disease on a chain-like network where node 0 is connected to node 1, node 1 is connected to node 0 and 2, and so on. In the simulation we are going through each node one by one from node 0 to node N and update them immediately. But it happens to be that node 0 was infected. Then, it can propagate to many nodes within a single time step because it is possible that node 1 gets updated to be infected by node 0, and then node 2 gets infected by node 1, and so on. This is clearly unrealistic. \n",
    "\n",
    "If we change a node's state, then the next node will be updating with respect to a network that is now in a different state!\n",
    "\n",
    "When modelling discrete-time dynamical systems there are generally two different update strategies: synchronous and asynchronous updating. In the asynchronous updating, a random node is picked and its state is updated according to the current network state. In the synchronous updating, there is a global time clock that all nodes are synced to, so nodes only update according to the state of the network at the _current time-step_ and all nodes move to the next time step simultaneously.    \n",
    "\n",
    "Choosing the updating scheme can have a huge impact on dynamics. We will be using the synchronous updating scheme, which means we need to store the _next state_ of the system before updating everyone all at once. There are many ways to accomplish this, but for the sake of simplicity, we will just keep the whole history. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "executionInfo": {
     "elapsed": 208,
     "status": "ok",
     "timestamp": 1648467910491,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "BunVNgrZBltG"
   },
   "outputs": [],
   "source": [
    "def run_SI(graph, tmax: int, beta: float, initial_inf: float):\n",
    "    \"\"\"Runs the SI model on the given graph with given parameters.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    graph : networkx graph object\n",
    "        The network (graph) on which the simulation will be run\n",
    "    tmax : int\n",
    "        The maximum time step that we will run the model\n",
    "    beta : float\n",
    "        The transmission probability\n",
    "    initial_inf : float\n",
    "        The initial fraction of infected nodes\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    list[float]\n",
    "        the time-series of the fraction of infected nodes\n",
    "    \"\"\"\n",
    "    # First lets generate a set of initially infected nodes.\n",
    "    init_infected = set(\n",
    "        np.random.choice(\n",
    "            graph.nodes(), size=int(len(graph) * initial_inf), replace=False\n",
    "        )\n",
    "    )\n",
    "\n",
    "    # The code below uses a dictionary comprehension to generate a dictionary\n",
    "    # with keys=nodes and values=infection states. The \"I\" is for infected\n",
    "    # and \"S\" is for susceptible. We then give that dictionary to networkx's\n",
    "    # attribute function which then gives all the nodes the 'state' attribute.\n",
    "\n",
    "    # YOUR SOLUTION HERE\n",
    "\n",
    "    # Now we need to loop through for `tmax` time step. One time step equals to\n",
    "    # updating the whole network once.\n",
    "    for t in range(tmax):\n",
    "        for node in graph.nodes():\n",
    "            # Now we check if the node if susceptible to infection\n",
    "            # If it is, we need to determine the probability of it switching\n",
    "            # and then switch it for the next time-step\n",
    "            if graph.nodes[node][\"state\"][t] == \"S\":\n",
    "                # First determine how many infected neighbors the node has at time t:\n",
    "\n",
    "                # YOUR SOLUTION HERE\n",
    "\n",
    "                # Instead of drawing a bunch of random numbers for each neighbor\n",
    "                # we can just calculate the cumulative probability of getting\n",
    "                # infected since these events are independent and then just\n",
    "                # draw 1 random number to check against:\n",
    "                if np.random.random() < (1 - (1 - beta) ** num_inf_neighbors):\n",
    "                    # If infection occurs we add a 1 to the state list of the node.\n",
    "                    # Note that by doing this we don't change how the other\n",
    "                    # nodes update, because they will be using time index t not t+1\n",
    "                    graph.nodes[node][\"state\"].append(\"I\")\n",
    "\n",
    "                else:\n",
    "                    # If no infection occurs, then just append the current state (0)\n",
    "                    graph.nodes[node][\"state\"].append(\"S\")\n",
    "\n",
    "            # Similarly, if the node is already infected it can't change back\n",
    "            # So we append the current state if it wasn't susceptible\n",
    "            else:\n",
    "                graph.nodes[node][\"state\"].append(\"I\")\n",
    "                \n",
    "    # the function returns a time series of the fraction of infected in each time step as a list.\n",
    "    # YOUR SOLUTION HERE"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "Z-DB-9N4BltH"
   },
   "source": [
    "And there we have our SI model. The function is mostly comments, there are only a dozen lines of code involved in the whole process. Lets give it a run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 303
    },
    "executionInfo": {
     "elapsed": 1081,
     "status": "ok",
     "timestamp": 1648468028065,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "3qjDbTvVBltH",
    "outputId": "3bb26a59-fd89-4339-a190-3b62923222ce"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGzCAYAAADT4Tb9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABKsElEQVR4nO3deVwUdeMH8M/uAsshLCg3gqB4g4AoPF6ZhWn18yyPvM20fLRUtMfsSc0uu/Qx07TMIzNTszRN05Q8shSUSzxBRUFuRFjOZdmd3x8kRSKyusvsLp/367WvR2Znhs8847YfZ+Y7IxEEQQARERGRmZCKHYCIiIhIn1huiIiIyKyw3BAREZFZYbkhIiIis8JyQ0RERGaF5YaIiIjMCssNERERmRWWGyIiIjIrLDdERERkVlhuiIiIyKxYiPnLjx8/jo8++gixsbHIysrCrl27MHTo0HqXOXr0KCIjI3H+/Hl4e3vjjTfewKRJkxr8O7VaLTIzM2Fvbw+JRPJwG0BERESNQhAEFBcXw9PTE1Jp/cdmRC03paWlCAoKwvPPP4/hw4ffd/7U1FQ8/fTTeOmll/DNN98gKioKL7zwAjw8PDBgwIAG/c7MzEx4e3s/bHQiIiISQXp6Olq2bFnvPBJjeXCmRCK575Gb+fPnY9++fTh37lzNtNGjR6OwsBAHDhxo0O8pKiqCo6Mj0tPT4eDg8LCxiYiIqBEolUp4e3ujsLAQCoWi3nlFPXKjq5MnTyIiIqLWtAEDBmD27Nn3XEalUkGlUtX8XFxcDABwcHBguSEiIjIxDbmkxKQuKM7Ozoabm1utaW5ublAqlSgvL69zmaVLl0KhUNS8eEqKiIjIvJlUuXkQCxYsQFFRUc0rPT1d7EhERERkQCZ1Wsrd3R05OTm1puXk5MDBwQE2NjZ1LiOXyyGXyxsjHhERERkBkzpy06NHD0RFRdWadujQIfTo0UOkRERERGRsRC03JSUlSEhIQEJCAoDqod4JCQlIS0sDUH1KacKECTXzv/TSS7h27Rr+85//4NKlS/jss8+wY8cOzJkzR4z4REREZIRELTdnzpxBSEgIQkJCAACRkZEICQnBokWLAABZWVk1RQcA/Pz8sG/fPhw6dAhBQUFYtmwZvvzyywbf44aIiIjMn9Hc56axKJVKKBQKFBUVcSg4ERGRidDl+9ukrrkhIiIiuh+WGyIiIjIrLDdERERkVlhuiIiIyKyw3BAREZFZMak7FBMREZFxEgQB5WoNisrV0AqAl2PdTw5oDCw3REREBADQagWUVFahqEyNonI1lOVqKCuq/3znpSyv+sfPf82j1lTfXeZfrZtj2zTxnh7AckNERGSmKtQa5BWrkFtcgVylCrl//rmgtLLOolJcUX3U5WFYSCUQ+w56LDdEREQmplRVVV1UlBXILVYhR1nxZ4mpXWSKytUPtH4rCykUNpZ3vRysLar/9+/T/jGPrZUMEolEz1usG5YbIiIiI1JUpsb5zKKa0pJ7p7T8WWBylBUordQ0eH1WFlK42sv/fFnD1UGOFnZyKGwsoLC1hIP1PwqMjSWsLWUG3ELDY7khIiISkUYrIPFmIY4n5+FYch4S0wsbdGrI1kpWU1hcHORw+7O4/L3EuNrLobCxFP1ISmNjuSEiImpk2UUV1WUmJQ8nUvLvOn3UqoUtPBU2cHWQw83BGq72crj8rbS4OVijmZxf4ffC/2eIiIgMTFWlwenU2ziekodjl/NwOae41vsO1hbo3dYZj7R1wSPtXOAp4jBqc8ByQ0REpGeCICA1vxTHkvNwPDkPp64VoFz913UyEgnQpaUj+rZzQd92zghq6QgLGe+rqy8sN0RERHpQXKHGH1dv1Vw7c/N2ea33Xe3leKSdC/q2c0Fvf2c42VmJlNT8sdwQERE9AK1WwIUsJY79WWbibtxG1d+uBLaSSdHdz6nmVFMHd/smd2GvWFhuiIiIGqCssgoXs5RIulmEhPRCnLiSj/ySylrz+Dnb4ZG2zujb3gX/at0Ctlb8mhUD/18nIiL6h+IKNS5kKpGUUYTzmUqcyyjC1bySu4Zo21nJ0NPfufp0U1sX+LSwFScw1cJyQ0RETdqdm+YlZRTh3J9FJjW/tM553RzkCPBUoLOXAj3btEBXHydYWfBCYGPDckNERE1GQWklzmUU/XlEpvp/0wvK65zXy9EGnT0dEOClQKCXAp29HOBqb93IielBsNwQEZFZyi2uwPmM6lNL5/58ZRZV1DmvT3NbBHg5oLPnn0XG0wEtmskbOTHpC8sNERGZjUvZSuyKz8BPiVnIKKz7iExrZzt09lIg0Muh+hSTpwIKW8tGTkqGxHJDREQmLUdZgR8TMrArPhMXs5Q106USoI1LMwR4Kapfng7o5OkAe2sWGXPHckNERCanRFWFA+eysTs+A79fzYfw5ygmS5kE/dq7YnhXLzzSzoVDsZso7nUiIjIJVRotfruSj11xGfjlQjYq1Nqa97q1csLQEC/8XxcPONryzr9NHcsNEREZLUEQkJRRhF3xGdibmFnrpnl+znYYFuKFocFevL8M1cJyQ0RERie9oOzP62gycDXvr3vONLezwuAgTwwN8UJQSwUfZ0B1YrkhIiKjUFSmxr6kLOyOz0DM9YKa6XILKfp3csPwrl7o09YFlnx6Nt0Hyw0REYlGVaXBkUt52B2fgV8v5aJSU30djUQC9GjdAsNCvDAwwJ0jnEgnLDdERNSo1Botzly/jb1nM7HvbBaKytU173Vwt8ewEC8MDvaEh8JGxJRkylhuiIjI4ApKK3H0ci5+vZSLY8l5KK6oqnnPzUGOIcFeGBbihY4eDiKmJHPBckNERHonCAIuZRfj10vVhSYu7XbNvWgAwMnWEo91cMOwEC/0aNMCMikvDCb9YbkhIiK9KK/U4OS1fERdzMWRS7l3Pcepo4cDHuvggsc6uCHY25GFhgyG5YaIiB5YRmE5fr1UXWZ+v5IPVdVfN9aztpSiVxtnPNbRFf3au8LTkdfQUONguSEiogbTaAUkpN9G1MXq002Xsotrve/laIN+HVzweAc39GjTAtaWMpGSUlPGckNERPUqKlPjWEoejlzKxdHLubhd9tfoJqkE6OrjhMc6uuKxDq5o72bPG+uR6FhuiIjoLjduleLAuWz8eikXZ27chkb719XADtYW6NveFY93cEXfdi5wsuOznMi4sNwQEREAIKuoHPvOZmFPYibO3iyq9V5b12bVR2fauyK0lRMseJdgMmIsN0RETVh+iQo/J2Vhb2JWrUceyKQS9GjdAv07ueGxDq7wbs4HU5LpYLkhImpiisrVOHg+G3sTM/HH1Vu1TjmF+TbHoCAPPBnoAedmchFTEj04lhsioiagrLIKhy/mYk9CJo4n59U8wwkAurRUYHCQJ54K9OBwbTILLDdERGaqQq3BseQ87E3MRNTFXJSrNTXvtXezx6AgD/xfF0/4OtuJmJJI/1huiIjMiFqjxR9Xb2FvYiYOnstGseqvZzi1amGLQV08MSjIE+3d7UVMSWRYLDdERCZOqxUQc70AexMz8fO5bBSUVta856Gwxv918cCgIE8Eeil4DxpqElhuiIhMkCAIOHuzCHsSM/HT2UzkKFU177Wws8JTgdWFplsrJ0j5DCdqYlhuiIhMiFqjxf6kLHz5WyqSMv66F429tQUGdnbHoCBP9GzTgvehoSaN5YaIyAQUllVia0waNv9xA9nK6qdtyy2kGPBnoXmknTPkFnyOExHAckNEZNSu5ZVg4+/XsTP2Zs1oJ+dmckzs0Qpjwn3QgveiIboLyw0RkZERBAEnr93ChhOpiLqUC+HPe+x19HDAlN5+GBTkwaM0RPVguSEiMhKVVVrsTczE+hOpuJClrJn+eAdXTOnthx5tWnC0E1EDsNwQEYmsoLQSW6Nv4KuTN5BXXD3qydpSimdDW2JyLz+0cWkmckIi08JyQ0Qkkiu5xVh/4jp+iLsJVVX14xDcHOSY0MMXY8J84GRnJXJCItPEckNE1IgEQcCJK/lYfyIVRy/n1UwP8HLAC71b46lAD1hZcBg30cNguSEiagQVag32JGRiw++puJRdDACQSICIjm54obcfwvya83oaIj1huSEiMqD8EhW2nLqBLaduIL+k+rEItlYyjOzmjUk9ffnQSiIDYLkhIjKAElUVlu6/iO9ib6Lyz+tpPBTWmNTTF6O7+0BhaylyQiLzxXJDRKRnyTnFeGlLLK7llQIAgloqMKVPazwZ4A5LPhaByOBYboiI9OjHhAy89n0SytUauDtYY9nIIPTk/WmIGhXLDRGRHqiqNHh330VsPnkDANDLvwU+GR0CZz4egajRsdwQET2kjMJy/PubOCSmFwIAZvbzx5z+7SCT8mgNkRhYboiIHsLx5DzM2haP22VqKGws8b9RQXisg5vYsYiaNJYbIqIHoNUK+PTXK1gRlQxBqL4J35qxofBubit2NKImj+WGiEhHt0srMXt7Ao4lV99h+LkwHywe1AnWlnxSN5ExYLkhItJBYnoh/v1NHDIKyyG3kOLdYYF4NrSl2LGI6G9YboiIGkAQBGyJTsPbey+gUqOFbwtbrBkXio4eDmJHI6J/YLkhIrqPssoq/HfXOeyKzwAAPNHJDR+PDIKDNe8yTGSMWG6IiOpxNa8E07fEIjmnBDKpBPMHtsfUPq15Uz4iIyb6fcBXr14NX19fWFtbIzw8HDExMfXOv2LFCrRv3x42Njbw9vbGnDlzUFFR0Uhpiagp2Z+UhSGrfkdyTglc7OXY+kI4pj3ShsWGyMiJeuRm+/btiIyMxNq1axEeHo4VK1ZgwIABuHz5MlxdXe+af+vWrXjttdewYcMG9OzZE8nJyZg0aRIkEgmWL18uwhYQkTlSa7R4/+dLWH8iFQAQ5tccq54LgauDtcjJiKghJIIgCGL98vDwcHTv3h2rVq0CAGi1Wnh7e+Pll1/Ga6+9dtf8M2fOxMWLFxEVFVUzbe7cuYiOjsaJEyfq/B0qlQoqlarmZ6VSCW9vbxQVFcHBgRcCElFt2UUVmLk1Dmdu3AYAvPhIa7w6oD0s+MBLIlEplUooFIoGfX+L9mmtrKxEbGwsIiIi/gojlSIiIgInT56sc5mePXsiNja25tTVtWvXsH//fjz11FP3/D1Lly6FQqGoeXl7e+t3Q4jIbPxxNR//9+lvOHPjNuzlFlg7LhQLnurIYkNkYkQ7LZWfnw+NRgM3t9q3KXdzc8OlS5fqXGbMmDHIz89H7969IQgCqqqq8NJLL+H111+/5+9ZsGABIiMja36+c+SGiOgOrVbA2uNX8fHBy9AKQAd3e6wZFwo/ZzuxoxHRAzCpf44cPXoU7733Hj777DPExcXhhx9+wL59+/D222/fcxm5XA4HB4daLyKiO4rK1Zj2dSw+PFBdbJ7p2hK7/t2LxYbIhIl25MbZ2RkymQw5OTm1pufk5MDd3b3OZRYuXIjx48fjhRdeAAAEBgaitLQU06ZNw3//+19IpSbV1YhIZOczizB9SxzSCspgZSHFksGdMbq7N0dDEZk40dqAlZUVQkNDa10crNVqERUVhR49etS5TFlZ2V0FRiarfpaLiNdFE5EJir52CyPXnkRaQRlaOtng+5d64rkwHxYbIjMg6lDwyMhITJw4Ed26dUNYWBhWrFiB0tJSTJ48GQAwYcIEeHl5YenSpQCAQYMGYfny5QgJCUF4eDiuXLmChQsXYtCgQTUlh4jofk6k5OOFzadRodbiX62bY+24UDjaWokdi4j0RNRyM2rUKOTl5WHRokXIzs5GcHAwDhw4UHORcVpaWq0jNW+88QYkEgneeOMNZGRkwMXFBYMGDcK7774r1iYQkYk5cikXL26JRWWVFn3bueDz8aF8mjeRmRH1Pjdi0GWcPBGZl4PnszFzaxzUGgERHd2wemwI5BYsNkSmQJfvbz5bioiahL2JmZi9PQEarYCnAz2wYnQwLHn/GiKzxHJDRGbv+9ibeHVnIrQCMCzECx8924U35iMyYyw3RGTWvo1Jw+u7kiAIwOju3nh3WCBkUo6IIjJnLDdEZLa++uM6Fu85DwCY0KMV3hzUGVIWGyKzx3JDRGbpi+NX8d7+6ke5TO3jh9ef6sh72BA1ESw3RGR2Po1KwbJDyQCAmf38MfeJdiw2RE0Iyw0RmQ1BELDsl2SsOnIFADC3fzu8/HhbkVMRUWNjuSEisyAIAt7ddxFfnkgFALz+VAdMe6SNyKmISAwsN0Rk8rRaAYv3nMfXp24AAJYM7oyJPX3FDUVEomG5ISKTptEK+O+uJGw7nQ6JBHhvWCCeC/MROxYRiYjlhohMVpVGi1d3nsWu+AxIJcDHI4IwvGtLsWMRkchYbojIJKk1WszeloB9SVmQSSX4ZHQw/q+Lp9ixiMgIsNwQkclRVWkwc2s8Dl3IgaVMglVjumJAZ3exYxGRkWC5ISKTUqHW4MWvY3EsOQ9WFlJ8Pj4U/dq7ih2LiIwIyw0RmYyyyiq88NUZ/HH1FmwsZfhyYjf08ncWOxYRGRmWGyIyCcUVajy/6TROX78NOysZNk4OQ5hfc7FjEZERYrkhIqNXVKbGhI0xSEwvhL21Bb56PgxdfZzEjkVERorlhoiMWkFpJcavj8b5TCUcbS2xZUo4ArwUYsciIiPGckNERiuvWIVxX0bjck4xnJtZYcsL4ejg7iB2LCIyciw3RGSUbt4uw/j1MUjNL4WrvRxbp/4L/q7NxI5FRCaA5YaIjM6V3BKMXx+NrKIKeDna4JsXwuHrbCd2LCIyESw3RGRUzmUUYcKGGBSUVqKNix22vBAOD4WN2LGIyISw3BCR0YhJLcCUTadRrKpCoJcCmyZ3R4tmcrFjEZGJYbkhIqNw5HIuXvo6FqoqLcL8mmP9xG6wt7YUOxYRmSCWGyIS3d7ETMzZnoAqrYDHOrjis7FdYW0pEzsWEZkolhsiEtW3MWl4fVcSBAEYHOSJZSODYCmTih2LiEwYyw0RiWbtsat4/+dLAICx4T54a0gAZFKJyKmIyNSx3BBRoxMEAR8dvIzPjl4FAEx/tA3+M6A9JBIWGyJ6eCw3RNSotFoBi/acw5ZTaQCA+QM7YPqjbURORUTmhOWGiBqNWqPFvO8S8WNCJiQS4J2hARgb3krsWERkZlhuiKhRVKg1mPFNHKIu5cJCKsHyUcEYHOQpdiwiMkMsN0RkcMUVarzw1RlEpxZAbiHF2nGh6NfBVexYRGSmWG6IyKAKSisxaWMMzt4sQjO5BdZP7Ibw1i3EjkVEZozlhogMJruoAuPWR+NKbgma21nhq8lhCGypEDsWEZk5lhsiMojr+aUYtz4aN2+Xw0Nhja+nhMPftZnYsYioCWC5ISK9u5StxLgvY5BfooJvC1tseSEcLZ1sxY5FRE0Eyw0R6VVc2m1M3ngaReVqdHC3x+YpYXC1txY7FhE1ISw3RKQ3J1LyMe3rMyir1KCrjyM2TgqDwpZP9iaixsVyQ0R6ceBcNl75Nh6VGi36tHXG5+NDYWvF/8QQUePjf3mI6KHtjL2J/+xMhFYAngxwx4rRwZBbyMSORURNFMsNET2Ujb+nYsneCwCAEaEtsXR4ICxkUpFTEVFTxnJDRA9s++m0mmIzpbcf/vtUR0ilfLI3EYmrQeVmz549DV7h4MGDHzgMEZmO+LTbWLj7PABgRr82mPdEe0gkLDZEJL4GlZuhQ4fW+lkikUAQhFo/36HRaPSTjIiMVl6xCtO3xKFSo8UTndwwtz+LDREZjwadGNdqtTWvX375BcHBwfj5559RWFiIwsJC7N+/H127dsWBAwcMnZeIRKbWaDFzaxyylRVo7WKHZSODeCqKiIyKztfczJ49G2vXrkXv3r1rpg0YMAC2traYNm0aLl68qNeARGRc3v/5EqJTC2BnJcMX40Nhb8372BCRcdF5SMPVq1fh6Oh413SFQoHr16/rIRIRGasfEzKw/kQqAGDZyGD4u9qLnIiI6G46l5vu3bsjMjISOTk5NdNycnLw6quvIiwsTK/hiMh4XMhUYv73ZwFUX0A8MMBd5ERERHXTudxs2LABWVlZ8PHxgb+/P/z9/eHj44OMjAysX7/eEBmJSGSFZZV4ccsZVKi1eKSdCyL7txc7EhHRPel8zY2/vz/Onj2LQ4cO4dKlSwCAjh07IiIigqMliMyQRitg1rYEpBeUw7u5DVaODoaMFxATkRF7oJv4SSQSPPHEE3jkkUcgl8tZaojM2P8OJeNYch6sLaVYOy4UjrZWYkciIqqXzqeltFot3n77bXh5eaFZs2ZITa2+uHDhwoU8LUVkZg6ez8aqI1cAAO8P74LOngqRExER3Z/O5eadd97Bpk2b8OGHH8LK6q9/wQUEBODLL7/UazgiEs+V3BLM3ZEIAJjcyxdDQ7xETkRE1DA6l5vNmzfjiy++wNixYyGT/fXU36CgoJprcIjItBVXqPHi12dQoqpCmF9zvP5UR7EjERE1mM7lJiMjA/7+/ndN12q1UKvVeglFROIRBAHzvkvE1bxSuDtYY/WYrrDkU76JyITo/F+sTp064bfffrtr+s6dOxESEqKXUEQkns+OXsXB8zmwkkmxZlxXuNjLxY5ERKQTnUdLLVq0CBMnTkRGRga0Wi1++OEHXL58GZs3b8ZPP/1kiIxE1EiOJ+fh418uAwDeHNwZIT5OIiciItKdzkduhgwZgr179+Lw4cOws7PDokWLcPHiRezduxf9+/c3REYiagTpBWV4+dt4CAIwurs3xoT7iB2JiOiBPNB9bvr06YNDhw7pOwsRiaS8UoNpX8eiqFyNIG9HLBnSWexIREQPTOcjN61bt8atW7fuml5YWIjWrVvrJRQRNR5BELDgh7O4mKVECzsrrBnbFXIL2f0XJCIyUjqXm+vXr0Oj0dw1XaVSISMjQy+hiKjxbPrjOnYnZEImlWD12K7wdLQROxIR0UNp8GmpPXv21Pz54MGDUCj+ulOpRqNBVFQUfH199RqOiAwr+totvLPvIgDg9ac64l+tW4iciIjo4TW43AwdOhRA9XOlJk6cWOs9S0tL+Pr6YtmyZXoNR0SGk11UgRlb46DRChgS7Inne/mKHYmISC8aXG60Wi0AwM/PD6dPn4azs7PBQhGRYamqNHhpSyzySyrRwd0e7w/vwgfgEpHZ0Hm01J0HZRKR6Vqy9wIS0gvhYG2BL8Z3g40VLyAmIvOh8wXFr7zyClauXHnX9FWrVmH27Nn6yEREBrT9dBq2RqdBIgFWPhcCnxa2YkciItIrncvN999/j169et01vWfPnti5c6deQhGRYSSkF2Lh7vMAgLn92+HR9q4iJyIi0j+dy82tW7dqjZS6w8HBAfn5+ToHWL16NXx9fWFtbY3w8HDExMTUO39hYSFmzJgBDw8PyOVytGvXDvv379f59xI1NfklKkzfEotKjRZPdHLDvx+9+wG4RETmQOdy4+/vjwMHDtw1/eeff9b5Jn7bt29HZGQkFi9ejLi4OAQFBWHAgAHIzc2tc/7Kykr0798f169fx86dO3H58mWsW7cOXl5eum4GUZNSpdFi5tY4ZBVVoLWLHZaNDIJUyguIicg86XxBcWRkJGbOnIm8vDw89thjAICoqCgsW7YMK1as0Gldy5cvx9SpUzF58mQAwNq1a7Fv3z5s2LABr7322l3zb9iwAQUFBfjjjz9gaWkJAPe9t45KpYJKpar5WalU6pSRyBws/fkSTl0rgJ2VDF+MD4W9taXYkYiIDEbnIzfPP/88li1bhvXr16Nfv37o168ftmzZgjVr1mDq1KkNXk9lZSViY2MRERHxVxipFBERETh58mSdy+zZswc9evTAjBkz4ObmhoCAALz33nt13jH5jqVLl0KhUNS8vL29G76xRGbgx4QMrD9RPcpx2chg+Lvai5yIiMiwdC43ADB9+nTcvHkTOTk5UCqVuHbtGiZMmKDTOvLz86HRaODm5lZrupubG7Kzs+tc5tq1a9i5cyc0Gg3279+PhQsXYtmyZXjnnXfu+XsWLFiAoqKimld6erpOOYlM2cUsJV77PgkA8O9H22BggLvIiYiIDO+BngpeVVWFo0eP4urVqxgzZgwAIDMzEw4ODmjWrJleA/6dVquFq6srvvjiC8hkMoSGhiIjIwMfffQRFi9eXOcycrkccrncYJmIjFVhWSVe/DoW5WoN+rR1xtwn2osdiYioUehcbm7cuIGBAwciLS0NKpUK/fv3h729PT744AOoVCqsXbu2QetxdnaGTCZDTk5Orek5OTlwd6/7X5ceHh6wtLSETPbXDcc6duyI7OxsVFZWwsrKStfNITJLGq2AWdsSkFZQBu/mNvj0uRDIeAExETUROp+WmjVrFrp164bbt2/DxuavpwcPGzYMUVFRDV6PlZUVQkNDay2j1WoRFRWFHj161LlMr169cOXKlZpHQQBAcnIyPDw8WGyI/mbF4WQcS86DtaUUa8eFwtGWnw8iajp0Lje//fYb3njjjbvKhK+vLzIyMnRaV2RkJNatW4evvvoKFy9exPTp01FaWlozemrChAlYsGBBzfzTp09HQUEBZs2aheTkZOzbtw/vvfceZsyYoetmEJmtg+ez8emvVwAA7w/vgs6ed9+XiojInOl8Wkqr1dY5OunmzZuwt9dtFMaoUaOQl5eHRYsWITs7G8HBwThw4EDNRcZpaWmQSv/qX97e3jh48CDmzJmDLl26wMvLC7NmzcL8+fN13Qwis3Q1rwRzdyQCACb38sXQEN4DioiaHokgCIIuC4waNQoKhQJffPEF7O3tcfbsWbi4uGDIkCHw8fHBxo0bDZVVL5RKJRQKBYqKiuDg4CB2HCK9KVFVYejq33EltwRhfs3xzQvhsJQ90IBIIiKjo8v3t85HbpYtW4YBAwagU6dOqKiowJgxY5CSkgJnZ2d8++23DxyaiB6cIAiYtyMRV3JL4OYgx+oxXVlsiKjJ0rnctGzZEomJidi+fTsSExNRUlKCKVOmYOzYsbUuMCaixvPZ0as4cD4bljIJ1owLhYs9b39ARE1Xg05Lde3aFVFRUXBycsJbb72FefPmwdbWtjHy6R1PS5G5OZ6ch4kbYyAIwHvDAjEm3EfsSEREeqfL93eDjltfvHgRpaWlAIAlS5agpKTk4VMS0UNLLyjDy9/GQxCA0d29WWyIiNDA01LBwcGYPHkyevfuDUEQ8PHHH9/zTsSLFi3Sa0Aiqlt5pQbTvo5FUbkaQS0VeHNwZ7EjEREZhQadlrp8+TIWL16Mq1evIi4uDp06dYKFxd29SCKRIC4uziBB9YWnpcgcCIKAOdsTsDshEy3srLD35d7wdOQ1b0RkvnT5/tZ5KLhUKkV2djZcXV0fKqRYWG7IHGz8PRVL9l6ATCrBlinh6NGmhdiRiIgMyqBDwf/+6AMianzR127hnX0XAQALnuzAYkNE9A8P9FTwlJQUHDlyBLm5uXeVHV5zQ2Q42UUVmLE1DhqtgCHBnpjS20/sSERERkfncrNu3TpMnz4dzs7OcHd3h0Ty15OGJRIJyw2RgaiqNHhpSyzySyrRwd0eS4cH1vr8ERFRNZ3LzTvvvIN3332Xz3MiamRv7rmAhPRCOFhb4PPxobC1eqADr0REZk/n+7Pfvn0bI0aMMEQWIrqHbTFp+DYmDRIJsPK5ELRqYSd2JCIio6VzuRkxYgR++eUXQ2QhojokpBdi0Y/nAQBz+7fDo+1Nc6QiEVFj0fm4tr+/PxYuXIhTp04hMDAQlpaWtd5/5ZVX9BaOqKnLL1Fh+pZYVGq06N/JDf9+1F/sSERERk/n+9z4+d17dIZEIsG1a9ceOpQh8T43ZCrUGi3GfRmN6NQCtHaxw48zesHe2vL+CxIRmSGD3ucmNTX1gYMRUcO9//MlRKcWwM5Khi/Gh7LYEBE1kM7X3BCR4f2YkIH1J6r/IbFsZBD8Xe1FTkREZDoadOQmMjISb7/9Nuzs7BAZGVnvvMuXL9dLMKKm6kKmEvO/PwsA+PejbTAwwEPkREREpqVB5SY+Ph5qtbrmz/fCG4oRPZzCskq8uOUMKtRa9GnrjLlPtBc7EhGRyWlQuTly5EidfyYi/dFoBczaloD0gnK0dLLBytEhkEn5DwYiIl3xmhsiI/G/Q8k4lpwHa0spPh8fCic7K7EjERGZJJYbIiNw8Hw2Vh25AgBYOjwQnT0VIiciIjJdLDdEIruSW4K5OxIBAJN6+mJYSEuRExERmTaWGyIRFZZVYtrmMyhRVSHMrzn++3RHsSMREZk8lhsikaiqNJj2dSyu5ZfCU2GN1WO6wlLGjyQR0cPS+Q7FAJCSkoIjR44gNzcXWq221nuLFi3SSzAicyYIAhZ8n4SY1AI0k1tgw+TucLGXix2LiMgs6Fxu1q1bh+nTp8PZ2Rnu7u617m0jkUhYboga4NNfr+CH+AzIpBKsHtsVHdz5nDMiIn3Rudy88847ePfddzF//nxD5CEye7vjM7D8UDIA4K0hndG3nYvIiYiIzIvOJ/hv376NESNGGCILkdmLSS3Af3ZWP1ph2iOtMTa8lciJiIjMj87lZsSIEfjll18MkYXIrKXml2La12dQqdFiYGd3vDawg9iRiIjMks6npfz9/bFw4UKcOnUKgYGBsLS0rPX+K6+8ordwRObidmklnt90GoVlagS1VOB/o4Ih5aMViIgMQiIIgqDLAn5+fvdemUSCa9euPXQoQ1IqlVAoFCgqKoKDAy/iJMNTVWkw/ssYxFwvgJejDXbN6AlXe2uxYxERmRRdvr91PnKTmpr6wMGImhpBEDB/51nEXC+AvdwCGyZ1Z7EhIjKwh7pjmCAI0PHAD1GTsuJwCnYnZEImleCzcV3R3t1e7EhERGbvgcrN5s2bERgYCBsbG9jY2KBLly74+uuv9Z2NyKTtir+JT6JSAADvDA1An7Yc8k1E1Bh0Pi21fPlyLFy4EDNnzkSvXr0AACdOnMBLL72E/Px8zJkzR+8hiUxN9LVbNUO+X+zbGs+F+YiciIio6XigC4qXLFmCCRMm1Jr+1Vdf4c033zT6a3J4QTEZ2tW8Egz/7A8UlavxZIA7Vo/pypFRREQPSZfvb51PS2VlZaFnz553Te/ZsyeysrJ0XR2RWSn4c8h3Ubkawd6OHPJNRCQCncuNv78/duzYcdf07du3o23btnoJRWSKKtQaTNt8BjdulaGlkw3WTegGa0uZ2LGIiJocna+5WbJkCUaNGoXjx4/XXHPz+++/Iyoqqs7SQ9QUCIKA/+w8izM3bsPe2gIbJ/Ep30REYtH5yM0zzzyD6OhoODs7Y/fu3di9ezecnZ0RExODYcOGGSIjkdH736Fk7EnMhIVUgjVjQ9HWjUO+iYjEovORGwAIDQ3Fli1b9J2FyCTtjL2Jlb9eAQC8OywAvds6i5yIiKhpa1C5USqVNVcmK5XKeuflCCRqSk5evYUFP1QP+f73o20wqjuHfBMRia1B5cbJyQlZWVlwdXWFo6MjJJK7R38IggCJRAKNRqP3kETG6EpuCV78+gzUGgFPd/HAvCfaix2JiIjQwHLz66+/onnz5gCAI0eOGDQQkSm4VaLC85tOQ1lRhRAfRywbEcQh30RERqJB5aZv3741f/bz84O3t/ddR28EQUB6erp+0xEZoQq1BlM3n0FaQRm8m3PINxGRsdF5tJSfnx/y8vLuml5QUAA/Pz+9hCIyVlqtgHnfJSIurRAOfw75dm7GId9ERMZE53Jz59qafyopKYG1tbVeQhEZq2WHLuOns1mwkEqwdnwo/F055JuIyNg0eCh4ZGQkAEAikWDhwoWwtbWteU+j0SA6OhrBwcF6D0hkLHacScfqI1cBAEuHB6JnGw75JiIyRg0uN/Hx8QCqj9wkJSXBysqq5j0rKysEBQVh3rx5+k9IZAT+uJKP139IAgDM7OePEd28RU5ERET30uByc2eU1OTJk/HJJ5/wfjbUZFzJLcGLW2JRpRUwKMgTkf3biR2JiIjqofM1NytWrEBVVdVd0wsKCu57gz8iU1NZpcUr38ajuKIKoa2c8NGzXTjkm4jIyOlcbkaPHo1t27bdNX3Hjh0YPXq0XkIRGYtPf03BhSwlnGwtsWZcVw75JiIyATqXm+joaPTr1++u6Y8++iiio6P1EorIGMSn3cbqI9XPjHpnaCBc7TkakIjIFOhcblQqVZ2npdRqNcrLy/USikhs5ZUazN2RCK0ADAn2xNNdPMSOREREDaRzuQkLC8MXX3xx1/S1a9ciNDRUL6GIxPbBgUu4ll8KNwc53hocIHYcIiLSQYNHS93xzjvvICIiAomJiXj88ccBAFFRUTh9+jR++eUXvQckamy/X8nHpj+uAwA+eKYLFLaW4gYiIiKd6HzkplevXjh58iS8vb2xY8cO7N27F/7+/jh79iz69OljiIxEjUZZocar3yUCAMaG++DR9q4iJyIiIl3pfOQGAIKDg/HNN9/oOwuR6N7aewGZRRXwaW6L15/qKHYcIiJ6AA9Ubu6oqKhAZWVlrWm8uR+Zql/OZ2Nn7E1IJMCykUGwkz/Ux4OIiESi82mpsrIyzJw5E66urrCzs4OTk1OtF5EpulWiwuu7qh+vMO2R1uju21zkRERE9KB0Ljevvvoqfv31V6xZswZyuRxffvkllixZAk9PT2zevNkQGYkMShAEvL4rCfkllWjvZs/HKxARmTidj7vv3bsXmzdvxqOPPorJkyejT58+8Pf3R6tWrfDNN99g7NixhshJZDC74jNw8HwOLKQSLBsZBLkF70JMRGTKdD5yU1BQgNatWwOovr6moKAAANC7d28cP35cv+mIDCyzsByL95wHAMyOaIsAL4XIiYiI6GHpXG5at26N1NRUAECHDh2wY8cOANVHdBwdHfUajsiQtFoB/9l5FsUVVQj2dsRLfduIHYmIiPRA53IzefJkJCZW3wfktddew+rVq2FtbY05c+bg1Vdf1XtAIkPZEn0DJ67kw9pSimUjg2Ah0/njQERERkgiCILwMCu4ceMGYmNj4e/vjy5duugrl8EolUooFAoUFRVx2HoTlppfiic/OY4KtRZvDuqESb38xI5ERET10OX7W6d/qqrVajz++ONISUmpmdaqVSsMHz78oYrN6tWr4evrC2tra4SHhyMmJqZBy23btg0SiQRDhw594N9NTU+VRou5OxJQodaiZ5sWmNDDV+xIRESkRzqVG0tLS5w9e1avAbZv347IyEgsXrwYcXFxCAoKwoABA5Cbm1vvctevX8e8efP4yAfS2efHryEurRD2cgt8NCIIUqlE7EhERKRHOl9kMG7cOKxfv15vAZYvX46pU6di8uTJ6NSpE9auXQtbW1ts2LDhnstoNBqMHTsWS5YsqRm5dS8qlQpKpbLWi5quC5lKrDicDABYPLgzvBxtRE5ERET6pvN9bqqqqrBhwwYcPnwYoaGhsLOzq/X+8uXLG7yuyspKxMbGYsGCBTXTpFIpIiIicPLkyXsu99Zbb8HV1RVTpkzBb7/9Vu/vWLp0KZYsWdLgTGS+VFUaRO5IgFojoH8nNzzT1UvsSEREZAA6l5tz586ha9euAIDk5ORa70kkuh3ez8/Ph0ajgZubW63pbm5uuHTpUp3LnDhxAuvXr0dCQkKDfseCBQsQGRlZ87NSqYS3t7dOOck8rDicgkvZxWhuZ4WlwwN1/vtKRESmoUHl5uzZswgICIBUKsWRI0cMnemeiouLMX78eKxbtw7Ozs4NWkYul0Mulxs4GRm72BsF+PzYVQDAe8MC4NyMfyeIiMxVg8pNSEgIsrKy4OrqitatW+P06dNo0aLFQ/9yZ2dnyGQy5OTk1Jqek5MDd3f3u+a/evUqrl+/jkGDBtVM02q1AAALCwtcvnwZbdrwRmxUW1llFebuSIRWAIaHeGFggIfYkYiIyIAadEGxo6NjzV2Jr1+/XlMoHpaVlRVCQ0MRFRVVM02r1SIqKgo9evS4a/4OHTogKSkJCQkJNa/BgwejX79+SEhI4OkmqtPS/Zdw/VYZPBTWWDy4s9hxiIjIwBp05OaZZ55B37594eHhAYlEgm7dukEmq/vhgteuXdMpQGRkJCZOnIhu3bohLCwMK1asQGlpKSZPngwAmDBhAry8vLB06VJYW1sjICCg1vJ3Hvnwz+lEAHA8OQ9fn7oBAPjw2S5Q2FiKnIiIiAytQeXmiy++wPDhw3HlyhW88sormDp1Kuzt7fUSYNSoUcjLy8OiRYuQnZ2N4OBgHDhwoOYi47S0NEilvC0+6a6oTI3/7Ky+L9OEHq3Qp62LyImIiKgx6Pz4hcmTJ2PlypV6KzeNjY9faDrmbE/ArvgM+Lawxf5ZfWBrpfPgQCIiMhK6fH/r/F/7jRs3PnAwosbyc1IWdsVnQCoBlo0MZrEhImpCeL6HzE5esQqv70oCALzUtw1CWzmJnIiIiBoTyw2ZFUEQsOCHJNwuU6ODuz1mRbQVOxIRETUylhsyKztjb+LwxRxYyiT436hgyC3qHtVHRETmi+WGzMbN22VYsvcCAGBO/3bo6MELxomImiKWGzILWq2AV787ixJVFbr6OOLFR3inaiKiporlhszCVyev4+S1W7CxlGHZyGDIpHwoJhFRU8VyQyYvJacY7/9c/RT515/qAD9nO5ETERGRmFhuyKSVqqow/Zs4qKq06NPWGeP+1UrsSEREJDKWGzJZgiDgv7uScCW3BK72ciwfGQyJhKejiIiaOpYbMllbY9KwOyETMqkEq8Z0hYu9XOxIRERkBFhuyCQl3SzCkj3Vw75fHdAeYX7NRU5ERETGguWGTE5RmRr/3hqLSo0WER1dMa1Pa7EjERGREWG5IZMiCALm7UxEekE5WjrZYNmIYEg57JuIiP6G5YZMyrrfruHQhRxYyaT4bGxXKGwtxY5ERERGhuWGTMbp6wX44MBlAMDCQZ3QpaWjuIGIiMgosdyQScgvUWHm1jhotAIGB3liXLiP2JGIiMhIsdyQ0dNoBczeloAcpQptXOywdHgg72dDRET3xHJDRm9lVApOXMmHjaUMa8aFwk5uIXYkIiIyYiw3ZNSOJ+dh5a8pAIB3hwWgnZu9yImIiMjYsdyQ0coqKsfs7QkQBOC5MG8M79pS7EhERGQCWG7IKKk1WszcGo+C0kp08nDA4kGdxY5EREQmguWGjNKHBy4h9sZt2MstsGZcV1hbysSOREREJoLlhozOwfPZWPdbKgDgoxFd0KqFnciJiIjIlLDckFFJu1WGed8lAgCm9PbDwAAPkRMREZGpYbkho1Gh1mD6N7EorqhCVx9HvPZkB7EjERGRCWK5IaPx1k8XcD5TCSdbS6wa0xWWMv71JCIi3fHbg4zC7vgMbI1Og0QCrBgdAk9HG7EjERGRiWK5IdGl5BRjwQ9JAICXH2uLvu1cRE5ERESmjOWGRFWqqsL0b+JQrtagt78zZj3eVuxIRERk4lhuSDSCIOC/u5JwJbcEbg5yrBgdDJmUD8QkIqKHw3JDotkak4bdCZmQSSVYNaYrnJvJxY5ERERmgOWGRHEuowhL9lwAAPxnQHt0920uciIiIjIXLDfU6IrK1Zj+TSwqNVpEdHTDtEdaix2JiIjMCMsNNSpBEDDvu0SkF5SjpZMNlo0IgkTC62yIiEh/WG6oUa377RoOXciBlUyKNWNDobC1FDsSERGZGZYbajSnrxfggwOXAQALB3VCYEuFyImIiMgcsdxQo8gvUWHm1jhotAIGB3liXLiP2JGIiMhMsdyQwWm1AuZsT0COUoU2LnZYOjyQ19kQEZHBsNyQwa377Rp+S8mHjaUMa8aFwk5uIXYkIiIyYyw3ZFBJN4vw0cHq62wWD+qEdm72IiciIiJzx3JDBlOqqsIr2+JRpRXwZIA7RnX3FjsSERE1ASw3ZDBv7b2A1PxSeCiseZ0NERE1GpYbMoj9SVnYfiYdEgnwv1HBcLS1EjsSERE1ESw3pHeZheV47fuzAIB/P9oG/2rdQuRERETUlLDckF5ptAJmb0+AsqIKQd6OmB3RTuxIRETUxLDckF6tPXYVMakFsLOS4ZNRwbCU8a8YERE1Ln7zkN7Ep93G8kPJAIC3hgTA19lO5ERERNQUsdyQXhRXqDFrWwI0WgGDgjwxvKuX2JGIiKiJYrkhvVi85zzSCsrg5WiDd4YGcNg3ERGJhuWGHtqPCRn4IS4DUgnwyehgKGwsxY5ERERNGMsNPZT0gjK8sescAODlx9qim29zkRMREVFTx3JDD6xKo8Xs7QkoVlUhtJUTXn7MX+xIRERELDf04FYduYLYG7dhL7fAilHBsOCwbyIiMgL8NqIHcuZ6AVZGpQAA3hkWAO/mtiInIiIiqsZyQzorKq8e9q0VgOEhXhgSzGHfRERkPFhuSCeCIOCN3eeQUVgOn+a2WDKks9iRiIiIamG5IZ38EJeBvYmZkEkl+GR0MOytOeybiIiMC8sNNdiNW6VY9GP1sO85EW0R4uMkciIiIqK7sdxQg6g1WryyLQGllRqE+TXH9Ec57JuIiIwTyw01yIrDyUhML4SDdfWwb5mUj1cgIiLjxHJD93Xq2i18dvQqAOD9Z7rA09FG5ERERET3xnJD9Sosq8Sc7QkQBGBkt5Z4KtBD7EhERET1YrmhexIEAa/vSkJWUQX8nO2weBCHfRMRkfFjuaF72nEmHfuTsmEpk2Dl6BDYyS3EjkRERHRfLDdUp6t5JXhzzwUAwNwn2iOwpULkRERERA1jFOVm9erV8PX1hbW1NcLDwxETE3PPedetW4c+ffrAyckJTk5OiIiIqHd+0l1llRaztsWjXK1BzzYtMK1Pa7EjERERNZjo5Wb79u2IjIzE4sWLERcXh6CgIAwYMAC5ubl1zn/06FE899xzOHLkCE6ePAlvb2888cQTyMjIaOTk5mvZL5dxLkMJR1tLLB8ZDCmHfRMRkQmRCIIgiBkgPDwc3bt3x6pVqwAAWq0W3t7eePnll/Haa6/dd3mNRgMnJyesWrUKEyZMuO/8SqUSCoUCRUVFcHBweOj85uZESj7GrY8GAHw+PhQDOruLnIiIiEi3729Rj9xUVlYiNjYWERERNdOkUikiIiJw8uTJBq2jrKwMarUazZs3r/N9lUoFpVJZ60V1KyitROSOBADA2HAfFhsiIjJJopab/Px8aDQauLm51Zru5uaG7OzsBq1j/vz58PT0rFWQ/m7p0qVQKBQ1L29v74fObY4EQcD8788it1gFf9dmeOPpTmJHIiIieiCiX3PzMN5//31s27YNu3btgrW1dZ3zLFiwAEVFRTWv9PT0Rk5pGjafvIFDF3JgJZPik9HBsLGSiR2JiIjogYh64xJnZ2fIZDLk5OTUmp6TkwN39/pPiXz88cd4//33cfjwYXTp0uWe88nlcsjlcr3kNVc/JmTgzb3nAQD/GdgenT057JuIiEyXqEdurKysEBoaiqioqJppWq0WUVFR6NGjxz2X+/DDD/H222/jwIED6NatW2NENVsHzmUhckciBAEYE+6DKb39xI5ERET0UES/5WxkZCQmTpyIbt26ISwsDCtWrEBpaSkmT54MAJgwYQK8vLywdOlSAMAHH3yARYsWYevWrfD19a25NqdZs2Zo1qyZaNthiqIu5uDlb+Oh0Qp4NrQl3hkSAImEw76JiMi0iV5uRo0ahby8PCxatAjZ2dkIDg7GgQMHai4yTktLg1T61wGmNWvWoLKyEs8++2yt9SxevBhvvvlmY0Y3ab+l5GH6ljioNQIGBXnig2e68H42RERkFkS/z01j431ugFPXbmHSxhhUqLUY0NkNq8Z0haXMpK8tJyIiM2cy97mhxhd7owDPbzqNCrUWj3VwxafPsdgQEZF54bdaE3L2ZiEmbTiNskoNevs747OxXWFlwb8CRERkXvjN1kRcyFRi/PoYFKuqEObXHOsmdIO1Je9lQ0RE5oflpglIySnG+PXRKCpXI8THERsmdedN+oiIyGyx3Ji51PxSjPkyGrdKKxHopcCmyWFoJhd9kBwREZHBsNyYsfSCMoxZdwp5xSp0cLfH5ufDoLCxFDsWERGRQbHcmKnMwnKM+fIUsooq0MbFDlteCIeTnZXYsYiIiAyO5cYM5SorMPbLaKQXlMO3hS22Tv0XnJvx+VpERNQ0sNyYmVslKoz9Mhqp+aXwcrTBN1P/BTeHup+YTkREZI5YbsxIYVklxq2PQUpuCdwdrPHt1H/By9FG7FhERESNiuXGTCgr1JiwIQYXs5RwbibH1qnh8GlhK3YsIiKiRsdyYwZKVVWYvPE0zt4sQnM7K2ydGo7WLnxCOhERNU0sNyauvFKDKV+dRuyN23CwtsDXU8LQzs1e7FhERESiYbkxYRVqDaZ9fQanrhWgmdwCX08JR2dPhdixiIiIRMVyY6Iqq7SYuTUOv6Xkw9ZKhk2TuyPI21HsWERERKJjuTFBVRotZm2Lx+GLuZBbSPHlxG7o5ttc7FhERERGgeXGxGi0AuZ+l4ifz2XDSibFFxO6oWcbZ7FjERERGQ2WGxOi1QpY8MNZ/JiQCQupBKvHdkXfdi5ixyIiIjIqLDcmQhAELNpzDjvO3IRUAqx8LgT9O7mJHYuIiMjoWIgdgOonCAJ+S8nH8kPJSEgvhEQCLB8ZjKcCPcSORkREZJRYbozYH1eqS82ZG7cBANaWUrw9JABDQ7xETkZERGS8WG6MUExqAZYfuoxT1woAAHILKcb9qxVe6tsGLvZ8ujcREVF9WG6MSOyN2/jfoWScuJIPALCSSfFcmDf+3c+fT/YmIiJqIJYbI5CQXoj/HUrGseQ8AIClTIKR3bwxo58/PPlUbyIiIp2w3IjoXEYRVhxOxuGLuQAAmVSCEaEtMaOfP7yb84neRERED4LlRgQXs5RYcTgZB8/nAACkEmBYSEu88rg/WrWwEzkdERGRaWO5aUQpOcVYcTgF+5KyAAASCTA4yBOzHm+L1i7NRE5HRERkHlhuGsHVvBKsjErBnsRMCEL1tKe7eGD2423R1s1e3HBERERmhuXGgK7nl2LlrynYHZ8B7Z+lZmBnd8zu3xYd3B3EDUdERGSmWG4MIL2gDKt+vYKdcTeh+bPVRHR0xeyIdgjwUoicjoiIyLyx3OhRZmE5Vh25gh2n01H1Z6l5tL0L5kS0Q5C3o7jhiIiImgiWGz35OSkLs7YloFKjBQD0aeuM2RHtENrKSeRkRERETQvLjZ6E+jpBKgXCfZojsn87hLduIXYkIiKiJonlRk9c7a3xy+y+8GnBm+8RERGJSSp2AHPCYkNERCQ+lhsiIiIyKyw3REREZFZYboiIiMissNwQERGRWWG5ISIiIrPCckNERERmheWGiIiIzArLDREREZkVlhsiIiIyKyw3REREZFZYboiIiMissNwQERGRWWG5ISIiIrNiIXaAxiYIAgBAqVSKnISIiIga6s739p3v8fo0uXJTXFwMAPD29hY5CREREemquLgYCoWi3nkkQkMqkBnRarXIzMyEvb09JBKJXtetVCrh7e2N9PR0ODg46HXdxobbar6a0vZyW81XU9reprKtgiCguLgYnp6ekErrv6qmyR25kUqlaNmypUF/h4ODg1n/Bfs7bqv5akrby201X01pe5vCtt7viM0dvKCYiIiIzArLDREREZkVlhs9ksvlWLx4MeRyudhRDI7bar6a0vZyW81XU9reprStDdXkLigmIiIi88YjN0RERGRWWG6IiIjIrLDcEBERkVlhuSEiIiKzwnKjo9WrV8PX1xfW1tYIDw9HTExMvfN/99136NChA6ytrREYGIj9+/c3UtIHt3TpUnTv3h329vZwdXXF0KFDcfny5XqX2bRpEyQSSa2XtbV1IyV+OG+++eZd2Tt06FDvMqa4XwHA19f3rm2VSCSYMWNGnfOb0n49fvw4Bg0aBE9PT0gkEuzevbvW+4IgYNGiRfDw8ICNjQ0iIiKQkpJy3/Xq+plvLPVtr1qtxvz58xEYGAg7Ozt4enpiwoQJyMzMrHedD/JZaAz327eTJk26K/fAgQPvu15j3Lf329a6Pr8SiQQfffTRPddprPvVkFhudLB9+3ZERkZi8eLFiIuLQ1BQEAYMGIDc3Nw65//jjz/w3HPPYcqUKYiPj8fQoUMxdOhQnDt3rpGT6+bYsWOYMWMGTp06hUOHDkGtVuOJJ55AaWlpvcs5ODggKyur5nXjxo1GSvzwOnfuXCv7iRMn7jmvqe5XADh9+nSt7Tx06BAAYMSIEfdcxlT2a2lpKYKCgrB69eo63//www+xcuVKrF27FtHR0bCzs8OAAQNQUVFxz3Xq+plvTPVtb1lZGeLi4rBw4ULExcXhhx9+wOXLlzF48OD7rleXz0Jjud++BYCBAwfWyv3tt9/Wu05j3bf329a/b2NWVhY2bNgAiUSCZ555pt71GuN+NSiBGiwsLEyYMWNGzc8ajUbw9PQUli5dWuf8I0eOFJ5++ula08LDw4UXX3zRoDn1LTc3VwAgHDt27J7zbNy4UVAoFI0XSo8WL14sBAUFNXh+c9mvgiAIs2bNEtq0aSNotdo63zfV/QpA2LVrV83PWq1WcHd3Fz766KOaaYWFhYJcLhe+/fbbe65H18+8WP65vXWJiYkRAAg3bty45zy6fhbEUNe2Tpw4URgyZIhO6zGFfduQ/TpkyBDhscceq3ceU9iv+sYjNw1UWVmJ2NhYRERE1EyTSqWIiIjAyZMn61zm5MmTteYHgAEDBtxzfmNVVFQEAGjevHm985WUlKBVq1bw9vbGkCFDcP78+caIpxcpKSnw9PRE69atMXbsWKSlpd1zXnPZr5WVldiyZQuef/75eh8ia8r79Y7U1FRkZ2fX2m8KhQLh4eH33G8P8pk3ZkVFRZBIJHB0dKx3Pl0+C8bk6NGjcHV1Rfv27TF9+nTcunXrnvOay77NycnBvn37MGXKlPvOa6r79UGx3DRQfn4+NBoN3Nzcak13c3NDdnZ2nctkZ2frNL8x0mq1mD17Nnr16oWAgIB7zte+fXts2LABP/74I7Zs2QKtVouePXvi5s2bjZj2wYSHh2PTpk04cOAA1qxZg9TUVPTp0wfFxcV1zm8O+xUAdu/ejcLCQkyaNOme85jyfv27O/tGl/32IJ95Y1VRUYH58+fjueeeq/fBirp+FozFwIEDsXnzZkRFReGDDz7AsWPH8OSTT0Kj0dQ5v7ns26+++gr29vYYPnx4vfOZ6n59GE3uqeCkmxkzZuDcuXP3PT/bo0cP9OjRo+bnnj17omPHjvj888/x9ttvGzrmQ3nyySdr/tylSxeEh4ejVatW2LFjR4P+RWSq1q9fjyeffBKenp73nMeU9ytVU6vVGDlyJARBwJo1a+qd11Q/C6NHj675c2BgILp06YI2bdrg6NGjePzxx0VMZlgbNmzA2LFj73uRv6nu14fBIzcN5OzsDJlMhpycnFrTc3Jy4O7uXucy7u7uOs1vbGbOnImffvoJR44cQcuWLXVa1tLSEiEhIbhy5YqB0hmOo6Mj2rVrd8/spr5fAeDGjRs4fPgwXnjhBZ2WM9X9emff6LLfHuQzb2zuFJsbN27g0KFD9R61qcv9PgvGqnXr1nB2dr5nbnPYt7/99hsuX76s82cYMN39qguWmwaysrJCaGgooqKiaqZptVpERUXV+pft3/Xo0aPW/ABw6NChe85vLARBwMyZM7Fr1y78+uuv8PPz03kdGo0GSUlJ8PDwMEBCwyopKcHVq1fvmd1U9+vfbdy4Ea6urnj66ad1Ws5U96ufnx/c3d1r7TelUono6Oh77rcH+cwbkzvFJiUlBYcPH0aLFi10Xsf9PgvG6ubNm7h169Y9c5v6vgWqj7yGhoYiKChI52VNdb/qROwrmk3Jtm3bBLlcLmzatEm4cOGCMG3aNMHR0VHIzs4WBEEQxo8fL7z22ms18//++++ChYWF8PHHHwsXL14UFi9eLFhaWgpJSUlibUKDTJ8+XVAoFMLRo0eFrKysmldZWVnNPP/c1iVLlggHDx4Url69KsTGxgqjR48WrK2thfPnz4uxCTqZO3eucPToUSE1NVX4/fffhYiICMHZ2VnIzc0VBMF89usdGo1G8PHxEebPn3/Xe6a8X4uLi4X4+HghPj5eACAsX75ciI+Prxkd9P777wuOjo7Cjz/+KJw9e1YYMmSI4OfnJ5SXl9es47HHHhM+/fTTmp/v95kXU33bW1lZKQwePFho2bKlkJCQUOtzrFKpatbxz+2932dBLPVta3FxsTBv3jzh5MmTQmpqqnD48GGha9euQtu2bYWKioqadZjKvr3f32NBEISioiLB1tZWWLNmTZ3rMJX9akgsNzr69NNPBR8fH8HKykoICwsTTp06VfNe3759hYkTJ9aaf8eOHUK7du0EKysroXPnzsK+ffsaObHuANT52rhxY808/9zW2bNn1/z/4ubmJjz11FNCXFxc44d/AKNGjRI8PDwEKysrwcvLSxg1apRw5cqVmvfNZb/ecfDgQQGAcPny5bveM+X9euTIkTr/3t7ZHq1WKyxcuFBwc3MT5HK58Pjjj9/1/0GrVq2ExYsX15pW32deTPVtb2pq6j0/x0eOHKlZxz+3936fBbHUt61lZWXCE088Ibi4uAiWlpZCq1athKlTp95VUkxl397v77EgCMLnn38u2NjYCIWFhXWuw1T2qyFJBEEQDHpoiIiIiKgR8ZobIiIiMissN0RERGRWWG6IiIjIrLDcEBERkVlhuSEiIiKzwnJDREREZoXlhoiIiMwKyw0RERGZFZYbIjIJR48ehUQiQWFhodhRiMjI8Q7FRGSUHn30UQQHB2PFihUAgMrKShQUFMDNzQ0SiUTccERk1CzEDkBE1BBWVlZwd3cXOwYRmQCeliIiozNp0iQcO3YMn3zyCSQSCSQSCTZt2lTrtNSmTZvg6OiIn376Ce3bt4etrS2effZZlJWV4auvvoKvry+cnJzwyiuvQKPR1KxbpVJh3rx58PLygp2dHcLDw3H06FFxNpSIDIJHbojI6HzyySdITk5GQEAA3nrrLQDA+fPn75qvrKwMK1euxLZt21BcXIzhw4dj2LBhcHR0xP79+3Ht2jU888wz6NWrF0aNGgUAmDlzJi5cuIBt27bB09MTu3btwsCBA5GUlIS2bds26nYSkWGw3BCR0VEoFLCysoKtrW3NqahLly7dNZ9arcaaNWvQpk0bAMCzzz6Lr7/+Gjk5OWjWrBk6deqEfv364ciRIxg1ahTS0tKwceNGpKWlwdPTEwAwb948HDhwABs3bsR7773XeBtJRAbDckNEJsvW1ram2ACAm5sbfH190axZs1rTcnNzAQBJSUnQaDRo165drfWoVCq0aNGicUITkcGx3BCRybK0tKz1s0QiqXOaVqsFAJSUlEAmkyE2NhYymazWfH8vRERk2lhuiMgoWVlZ1boQWB9CQkKg0WiQm5uLPn366HXdRGQ8OFqKiIySr68voqOjcf36deTn59ccfXkY7dq1w9ixYzFhwgT88MMPSE1NRUxMDJYuXYp9+/bpITURGQOWGyIySvPmzYNMJkOnTp3g4uKCtLQ0vax348aNmDBhAubOnYv27dtj6NChOH36NHx8fPSyfiISH+9QTERERGaFR26IiIjIrLDcEBERkVlhuSEiIiKzwnJDREREZoXlhoiIiMwKyw0RERGZFZYbIiIiMissN0RERGRWWG6IiIjIrLDcEBERkVlhuSEiIiKz8v/GEVQfygpk1wAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "G = nx.barabasi_albert_graph(1000, 4)\n",
    "plt.plot(run_SI(G, tmax=20, beta=0.05, initial_inf=0.1))\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"fraction infected\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try aggregating results from multiple runs. Let's do 50 runs and plot them all to see the variation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 303
    },
    "executionInfo": {
     "elapsed": 1697,
     "status": "ok",
     "timestamp": 1648468087684,
     "user": {
      "displayName": "Yong Yeol Ahn",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "00405984523192108505"
     },
     "user_tz": 240
    },
    "id": "mf5smaShBltI",
    "outputId": "a9f6dc8a-f848-40d9-9db3-42042255d08a"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x340b39cd0>]"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACRC0lEQVR4nOz9ebTk913f+T+/a33rW3vV3W/f3rTa2MjGi7ATAyZyjOMYPBMmTsJgxyeT/JIDHAbNnMHOgHMyyeBMGIhnxiZOHMCAIZg4YRkWO0aDIQ4CgY0wtrVY6r3vfmuvb3337++P0uejut0tWy1L6u39OOceWber+tZtS7de+nzei1EURYEQQgghxDViXusXIIQQQohbm4QRIYQQQlxTEkaEEEIIcU1JGBFCCCHENSVhRAghhBDXlIQRIYQQQlxTEkaEEEIIcU1JGBFCCCHENWVf6xfwbOR5zubmJrVaDcMwrvXLEUIIIcSzUBQFo9GItbU1TPOZzz9uiDCyubnJxsbGtX4ZQgghhHgOzp8/z5EjR57x12+IMFKr1YDZN1Ov16/xqxFCCCHEszEcDtnY2NDv48/khggj6mqmXq9LGBFCCCFuMF+rxEIKWIUQQghxTUkYEUIIIcQ1ddVh5A/+4A9429vextraGoZh8Gu/9mtf8zmf+cxn+KZv+iZKpRK33347H/3oR5/DSxVCCCHEzeiqw8hkMuGee+7hQx/60LN6/OnTp3nrW9/KG9/4Rh5++GH+x//xf+R/+B/+Bz71qU9d9YsVQgghxM3nqgtY3/KWt/CWt7zlWT/+wx/+MCdOnOAnfuInAHjJS17CZz/7Wf7Vv/pXvPnNb77aLy+EEEKIm8wLXjPy4IMPct999x363Jvf/GYefPDBZ3xOFEUMh8NDH0IIIYS4Ob3gYWR7e5vl5eVDn1teXmY4HDKdTq/4nPe///00Gg39IQPPhBBCiJvXddlN8973vpfBYKA/zp8/f61fkhBCCCFeIC/40LOVlRV2dnYOfW5nZ4d6vU65XL7ic0qlEqVS6YV+aUIIIYS4DrzgJyOve93reOCBBw597tOf/jSve93rXugvLYQQQogbwFWHkfF4zMMPP8zDDz8MzFp3H374Yc6dOwfMrlje+c536sf/w3/4Dzl16hT/y//yv/Doo4/yUz/1U/zKr/wKP/RDP/T8fAdCCCGEuKFddRj50z/9U175ylfyyle+EoD777+fV77ylbzvfe8DYGtrSwcTgBMnTvBbv/VbfPrTn+aee+7hJ37iJ/h3/+7fSVuvEEIIIQAwiqIorvWL+FqGwyGNRoPBYCCL8oQQQojnqCgK8jynKIrLPkql0tdcaHe1nu379w2xtVcIIYQQh10pUHytj0vleU6e56Rpim3b2Pa1iQUSRoQQQojrjDrBmP/4WsHia/1eaZoShiFxHJMkif5rlmVkWcbLX/5yCSNCCCHEreZKoUN9PBvz1ypFUeiQEccxURQxmUxIkoQoikiShDzP9a+HYUgYhkRRRBiGHD16lEql8kJ9q1+VhBEhhBDiBVYUBVmWXfG045kYhnFZDYcKHNPplOl0ShAETKdToigiiiKCINBBJMsyptMpWZYRRRFpmhLHMXme69OQJElI05Q0TQmC4IX+Y3hGEkaEEEKI58kznXJ8tdBhmiaGYejHZFlGHMdMJhMmkwnj8ZjRaEQYhozHY+I4Jk1THSSyLNOBIk1THWDU5w3D0K8hiGIGic0gsxnlPoHhERhlQqvC2/aG3H77i/LHdBkJI0IIIcRVUtcr6rRD/fWrPV4FDlW/EccxQRAwmUwIgkD/bxU81OkGPB1QDMMgyzIdONTXVicvUZIxLlwOQhgXLoN0FjamZpnIrpLYPjjMPi5xrisnI0IIIcR1af6K5dK/XvoxHwxUfYa6Ppm/UgnDUJ9cRFF0qKg0z3MdXNQJCMyCRxgnTIoSw8xmmLuMc5fArDM1fSKrQuI8VfPhPvXir7BZxcwTvHSCl08oZwE+IX4RctvC216cP9ArkDAihBBCPGU+VKiaivngMX8doj6najhUqFAdK+rv1VWKChvzNRwAlmXpk5MwDEmSlEFm00tdDtIyQzqMjQqRXSG2fDCMK55sKGaeUErH+EVIhZCqEdNyc1puzkLZpOHZ+H4Zz/OwLB/Hmf1mdx0/8mL8EV+RhBEhhBC3pPlwoeovVMeJChCXFnmq6xXVnTLfKquCizoZUc+L4xh4+mpHf+0ceqnDXmRykLj0C4+xWWFq1chN+xmvU8wi1ScadTOh4WQ0rIRO2WC9WWah6lGrVbHtJpZlYRgGpmkeClnqmseyLExzNozddd3Lv9iLRMKIEEKIm54KBfPXIZeGj/lrEhUmLg0iKqTkea6vWlQ4Udcq6npGF6SaDt3U5SBx6aYu/bzEyKgQmD4YJhg8fa3yFKPIqOQTGkxpWRELpYzlikW7VNAqO9TrNVy3geu6mKaJaZqkaXoo8MzPELFtW5/AOI6DZVmzr3OFUHItSBgRQghxU5kvEJ1OpwyHQ309cul1iQoPapCYesx8F4waFjadTvWVjPr1KIqe+qoG08LiIHEYUKOflejnHkPKTA3v6RdncmgrnJXH+OmIWjGhZUUsl2GplLFScyh7JSqVCqVSTYcGwzCwLEsHjiiK9OkHoKeo5nmOZVnYtq2vgEzTxHVdbNumVCrpx9q2jeM4+rrmWpAwIoQQ4oY1HyBUvcZkMtFFolEU6XoOVRiqThJM0zwUUC4NMaPRSJ+gqK8zCy3QTWwO8jI92uznPr2iQnTpncrciBA3m+KnI6r5mDoBC27KasVgqeZRKrn4vo/vt7AsS59SzLfoqs85jnPZoDNAhwl1AqL+qk5BXNc99BgVQNRjr9XkVUXCiBBCiBvCfDutChAqdKhulVkB6NP1H4Zh4LounufpwKJORVTx6Xg8PvTc+TbdNMvpZSW6RYWDokI3r9AtfBKsK71A/CKgko7wsxGVbETDmLJYymmUHarVKtVqFd9folyeFZCq0KAChgoXtm3roWfzgUMFB8uyKJVKOqSYpnkoWKjwMf/4a30V89VIGBFCCHFdmr8ymW+DVa2yKkTM136ox89vp51MJrqmQp2cTKdTJpPJoa6YrDDoFWW6RZ1uUaFHlW5eJuPyN3CzyGgUI6pxj1rWp5GPqBcTquUSZb9MrVajUmnRaBzH8zx839cnE/OtwOqkBmZBRIUQ13X1iYYKFoAOHSp4qF9Tj7/0VOVGIWFECCHENTd/4jFfKKpGng+HQ4bDoe5iUQGkKAp93QKzEwXTNPUJhwof3W5X/94AuTG7ZukbLfpGlR5VenmZnMvfxG1SmsWYRjagmvSopX2qRYDvlbAci9ZSi0pljUajged5VCqVQ6ce6pRFBaT5QKFOPxzH0c9R1yfzpx2lUkmHmfnTkBstdDwTCSNCCCFeNPOL4VQ4mN+noq5c1FTS+c2y6iRBjTdXnSTqBCSKIl2sqpbAJUnCJEroG3W6Ro2BUaNbVBjkHgUGFMw+nlIipWNOaBsBrWJEPR/gRAPyLJudbpQdarVFKpUK1WqVWq1GtVrFtm3K5bKuLYGnx7yr1w1PF6CWSiVKpZI++Zg/6ZgPJPNXODczCSNCCCFeEFcamT4/p0MViqorl/nBYfPhQ11XqI2yqtBUdcrMb6bN85xBlLOdVtgrVtijTtfwyTEvCx5VM2XRDlm0QzrWlFYxwkkmpOms1sR2ZqHAWzxCuVymXC5TqVRoNpuHQoQKSGp+hwpIpmli2zae5+F53mXhYz5w3EynHM+FhBEhhBBftyuNTFf1GCqAqJZadQKilsCp2R4wO02o1Wq4rkupVKIoCv2cXq+nR6oHQfBUQWpCv/DYSsrs5m32qDPML5+BXrdz1ryYJTtk0Y5YtEPsNNDBp8iKWTAoezhOTddhlEolms0mtVqNer2uu1lUsJoPEb7vUyqV8DyPcrl8WR3H/F/FYRJGhBBCXJX5a5Zn+ut8UakaHqaChGq9VW/m6g3cMAxs29ZXNfv7+/rKRp2apIXBQVFlO2uyk1XYzStExeG3MoOCpVLGES/mSClixZ5QM1PSdBaG4jAmznNS08TzPF1AWiqVKJfLVKtVHUB83z+0lE69xnK5rGs4rtS5Mh8+xNcmYUQIIcTXpALFpUvh5geEAYdGoUdRxGQy0R8qpKgrilJpdoJhGAaj0YjpdEoYhozHY12oOklNdvIKu/nS7Oolu7zI1DYK1koR627ImhuyaIxxeGqi6jgmynPCotC1GpVK5dBrqFartFqtp9pufX2iUTz1nEtPNeaHhF1aUCqeGwkjQgghrqgoCn21ojpYVPiYr3NQRZrq+mU8HjOZTAjD8NBuFtU1ArPrGFVoOhgMGI/HhGHEQWKzk1fZzWcnH4MrXLlUzJQVO2DFnrBiBbSMCUWWQgZRPyK0baZPjT2fP8VQtRv1ep16vU6z2aRUKunaDnV9okKH+t4ubbGVE4/nn4QRIYQQh6jdKvO7XAD9ZqxCiro6UcFjNBrpLhbVQmsYBuVyGcMwKJVKJEnCaDRiOBzOnhOmPDn1OJMsspnXL7tygYK2FbFsjlm2Ao6UQrwsoCieGtUepcTMBoKZpkmr1dI1G+VymXq9roeNqSJS9T2qiaaXjkFXJybz8zvEC0vCiBBCCAA9oVTtYMmyTL9hA4eKT9Xpx2Aw0J+D2amCKt6c7zKJooiDgwPiOKYXFnx5aPNkuMRWVp11ujzFJmfRClh1AlbtKRvlBLt4ehsugOk6OtyoIlN1xaL+WqlUKJfLWJalr5PUjhZV26Fenzr1UAHkVu5quVYkjAghxC1ufhHcpYvjVDBRRaX9fl8XoKqJobZt60JPda1hGIY+PZlOQ7YmOV8eOjwxbbGXVw59/ZYZcns54A4/ZKUUQ57pBXS2aWMY9qHuFNd1dehQwWO2UO7pMKFqVpIkuWKdhwogKpSIa0vCiBBC3IJUy6wqGp1fLKeKUwG9MG48HpMkiS7UVFcv83Ujao9KEARMgilnxgaPDB0eD1qXtNsWLFsBt3sT7q7FtJ2U0WhEySyRp7OTCjU+HdAzRtSgsfkWWjVxVRXNBkGga1rUdYt67PxcEHF9kTAihBC3kCzLmEwmDIfDQx0s8yvnAT2ELI5j3VUyX+A5v0dFbbYNooTHBgZf6lt8ZdJmOlf/YZJzxB5zZyXkzmqEV8yKV63MIrdclpeXD12ZlMtlms2mDhHz7bOqgFTVtqhrGDVg7NIAIq5/EkaEEOImdOkskDiOOTg4YDAYEIbhoZAx32Y7v/V2fpdLmqaHVs2rQDKYJnxhP+cvuhZPTMokxdP1H66RcdQecltpzO2VGNecjWy3YgvDdel0OvqUpVqt0ul0KJfLOuSor6e206qrI7W7pigKHVzUZt75ehBx45AwIoQQN4H5IWNq5keSJAwGA13noQKG+dSwr0ajgeM4erR6HMeMx2NGoxFRFOlrGPV4VXexN0n53HbCnx8UnJ64sx0vT6kYCcfsPsedISerGY41a+F1zdkViWqndV2XWq3GwsLCoSsfFYxUASw8PbtEXR2VSiVqtZoOHup1SQC5cUkYEUKIG9T8dFN1YqBqPwaDgS5IVW/mvu/TbrfxPA+AyWRCr9fTrbaj0Yg8z/UsDdu2qdVmo9H3pgW/ey7kz3YjLk4t4Ol22JY55bg94LgzZMVNsO3ZplrLmAWFer2ur1DUjA8VONTXUmFHdbyoOR5q30upVNLBQ4UW9b/FjU/CiBBC3EBU6JjfZKs6XdRMELW3pSgKPWujXq9jmibj8Zi9vT36/b4OIXEc642xlmXh+z7lcpncdPjs2Ql/cGHCqZF607eAghVrwobZ47bShIYV63oSVXyqBoyp65NKpYLnebr4VHW0qOue+Q216npo/jWpgWnzw8nEzUPCiBBC3ABUt4g66YjjWAcQdSIymUx0DYjaraJ2vezu7tLr9RgMBnrXS5qmlMtlfV2jdrM8sp/wmcenfG4vIM4NwAQK1q0Rx80DbvMCSiQ6QLjuLLzUarVDM0YqlQqu61Kv13VhqTpxUfUf88WoaleNqhNRdSPq9ETcvCSMCCHEdUpNOp2vBQmCgDAMD+1/UfteOp2Ovv7I85yDgwOGwyHdbvfQeHbbtmk0Gnorbrlcph8bPHA24g83h+yFqvbCoGaE3GHtc5fbp2Im+qrE8xqUy2V9ilKpVCiKQneyqCVzapaHKnydHz4G6FMZdVIC6L9XzxE3P/l/WgghrjPqGkZ1sqhrGDWQTF3HqKsLtdq+KAp99TIYDPTI9SzLsCyLSqWix6YXRUFhWPx5Fz77pZBHe8VThagGDhnHrS63W/us2lMsa7arpVyu6ZChgoh6DaVSiUajoaeeqhBiWZb+fgB9DTNfrDq/gE46YW5NEkaEEOI6cKVi1CAI9CRSNSVVXWOok4RSqcRkMmFra4vRaKQLV5MkObRb5emC0JzN0OGPduCh3Yxpql6BwZo94aSxywb7VD3nqXqP2awPdQVTLpd1sallWdTrdcrlMqZp6oJU1aGjFuGp8KOuXTzPOxQ+5AREyD8BQghxjTxTMep0OtXtuepqRe1fAXTh52Aw4OzZs3rvy2Qy0dcbanqpKgIdxvBnXZsHd2BzUujXULdS7rAPOFbs0LDipwJHB8dxqNVq+tpFTUBVpxnqfwN6I+58MS2gC09VHcv8lY3UgIh5EkaEEOJF9kzFqGma6muM6XSqB4Kpgk/LsnQtyMHBAePxWLfkqiml6vrD930KDP7iAP7knM0X9lKypxbN2UbBHeWAY/km7WQf17Cp1qpUqy3dDaNqT9RAMrV4Tn0NdRJiWRZhGNLr9fRUVnUSMz8zRNpwxVcjYUQIIV4E6hQkjuPLilHVXhfVGeO6Ls1mU19vAERRxN7eni5KHQ6HZFmG67o0Gg3dvVKtVrk4gd/bzPnDzYRBmAOzu5gj5ZQ7rH0WgrPY0wTf92l2lnUNiGq5tW1btwSr4tT54lLbtkmShNFopMewz++RUZ0zUv8hni0JI0II8QJSHTFq/Lq6xlBv4gBBEOgpp7VaTYcQVTcyGAx0TUgcx7poVD22Wq1iuh5/vFPw+48mfKUb669fcwq+wZ+wnpzHnuzpOo9ms0m1WtWnG6rWQwWaRqOB7/uH9tHkea5bgtXJhwovaraI1H+I50L+qRFCiBfApQWp0+mUyWSip4tmWUYcx7obZn7eRhiGTCYT9vf32dnZYTQa6cJV3/d1y6zneaRY/N75hN96MqAfzupMTANeWk85aexSGZ4l7yd4nkdzZUWHDDULxLIsXY/SaDR0R4z6HtT8jzAMddGpCh2qoFUCiPh6yT9BQgjxPJpvvQX0dYa6iomiCNM09TUGoJfABUFAv99nZ2eH/f19kiTRE1FVYFCDw8Zhyu+czfmdU1OG0SyEtEsFr20GrITnmPb3AKjW69Try7TbbXzfB9AhxPM8qtUqCwsL+loojmO9CVed5qjgoYKIqh0R4vkiYUQIIZ4H6s1bzQbJskxfxxiGQRAEetYHcGj/ynA4ZG9vj93dXYbDod4PU6lUqNfrVKtVPM8jjmN6k4jfu5DyydMR43hWT7LowWtrPdqDJwg3J2Sex/Lysp4JUq1WL2utbbVaLC0t6d9XdcCoa6X5Dhg1xl2W0YkXioQRIYT4OqiiVLVRVp2MqELVMAwpioJKpaLbdV3XJU1T9vb22NnZodfrEUUReZ5j2zaVSoVGo0G1WsVxHOI4Zrs74tNnEz59LiVIZiFkuVzw2mqP1vBJptsTTN9nbW1Nj4Ivl8u6xkRtyW00GrTbbQAdnizL0ltx1fh21VUzvz1XiBeKhBEhhHgOVFGqChiqvkIFERVQVCtsnueUSiXiOOb8+fNsb28znU4PPa7RaNBqtXTh6HQ65dxOj985HfH/nc8I01kIWSkX3Fvr0R6fIdgdk1oWKysrtFotFhcX9XVLlmVUKhV836fT6eiR7XEc670yapKrWmSnilJVd4wQLwYJI0II8SxduitGmS/yVCchakZIURT6dOPUqVN0u109IVWdmHQ6nUObbqfTKU9e3OM3nwz5zPmUaHbowppf8Lpaj870ApODMYlhsLy8rEOM67r6iqVarVKtVnWtiBqqpuaAFEWhT0FUYarqrBHixSZhRAghvgZ1mqACxLwsy3S7rqrLUKchMKsN2d7e5sKFC0wmEz28bHFxUV+HqCuROI559Ow2v/romN+/kJI8lXc2KgV/qTmgPd1k3BuRWBbtdptGo0GlUqFUKununXq9rq9i5gtkPc/Tm3LVXBDVnitXMeJakzAihBDP4NL2XEB3xah6kPF4TJ7nevS5ut4wDIPBYMC5c+fo9XokSYLruiwuLtJsNnUxaa1WI89zHj2/y698ocd/uZjx1G0Mxyo539KZ0JxuEvQnZLat54OUy2Ucx9FL6SqVCs1mk1arpQeTqZChXr+aCaKmp8pQMnG9kDAihBCXuLQ9F57es5JlGVEUMRqNdPGnOmkwDIOiKIiiiLNnz7K3t0cYhti2zeLiIp1OR9djqLqQL57Z5pf+bJ8/3MrIilkwOFnLeePilEa4w2QwhqdOTlRXi5qUOt9ts7CwgG3buvVX7bxRc0TUnBK5ihHXIwkjQgjxlEvbcwE9jEztkxmPxzpgqADiOA55nhNFEbu7u/pKxjAMms0my8vLerOt7/uYpsmXzu3xC5/f44+3MnIMwOD2WsZfWYnpZN3Zrpqnil5VDYgKIbVaTV/TNJuzrbppmuptuUmS6AFp6vWpaxkhrkcSRoQQt7wrhRB1naHe3MfjsZ4Voq4/XNfV1zUHBwdsb2/T7/d1F8vKygrNZlOfqpimyaNbfX7uoW3+eDujeCqE3FFL+Y4jOUvmmOl0ShiG+upHLahTY9wXFhZot9vU63U8z9ML92C24dfzvEOL6mQ2iLgRSBgRQtyy1GnGfAhRdRiqdTcMQz2OvVQqURQFpVIJQG/NVVt0oyjCcRxWV1dZXFykXC7roWd744T/+w/O8dBcCLmzGvPW4yYrzmz8+2AU6LHrakeMbdtUq1VWVlbodDo0Gg19zTKdTsnzXNeD+L6vr2JkRLu4kcg/rUKIW86VakJUcWeSJIRhSJqmDAYDAB1C1FXHcDhkNBoxGo3odruMx2MAOp0Oq6ur1Ot1YDaLJMfkF/7oAv/xsYA4n4WQuyohbztps+rlhGFAtzvUE1Cr1Sr1ep1SqUSlUmFhYYGFhQVarZZuyY2iSJ/SqAV7akS71IOIG5GEESHELeNKIUTVfqRpqk8aVHGq4zi6FVaFkH6/r4PKwcEBeZ5TrVZZWlpicXER27aJ49nW3D863edfP7TPznQWQtbckL91p8Wxmk0YhvT7YwaDgZ430mg0qNVqeJ6nr2TUEDR1HaPGxav2YTVpVepBxI1MwogQ4qanThMuDSFqLPt0OtWL6lRxqnqzdxyHfr9Pr9fTM0XUlYyamrq+vq53zoRhyIXuhA/+1y0+v1cABmUz5S2rEa9fs3Fsg/F4zN7eHmmaUi6XaTabejiZ7/vU63WazSaNRgPTNJlMJoRhyHQ61Vcw6nESQsTNQMKIEOKmpYaVqZMKePokJMsyHULml9qpEGIYBpPJhH6/r0erq6sZwzCo1+usrKzo05DpdEp/OOaXPr/DbzwZk+QGBgWvaU55+2021dJssurW5kXiOKZcLrO4uMjCwoJehOf7PpVKRQcTVTSrvn6pVKJUKulhaULcLCSMCCFuOlcKIWrehjrdUEPMkiQhiiIAHUKm0ymDwUA/Tv29uh5pt9scOXJEX5/0ej1+/7FdfvbPx+xFJmBwpBTx3Sdy7lyqkCQJ29vbjEYjSqUSS0tLeqEdQLlcplQq0Ww2aTabJEnC3t6e7vKZv46RKxlxM5IwIoS4aVxpbLsKIXme65MQRQURdTqSpqk+jVBdKupkBKBer3PkyBGWlpZIkoTBYMCXTm/x0S+M+PMDAzCpmBlvWQv5tuMVLMtkf3+fvb09iqKg0WiwsbHBysoK8HTRrO/7LC0t6amtqrU3jmN9aqJORYS4GUkYEULc8NQCuziOddhQrbiAfnOHp8e5q425auR7FEX69MMwDIIgYDqdkmUZtm2zurrK+vo6juMwHA65uL3LLz+8z6fOFSTF7Erm3taUt9/u0qo2GA6H7O7uMhqN8H2fjY0Njh07Rrlc1jNHLMvSS/JUh4wKRbZts7CwoLtlpFVX3Mzkn24hxA3rq4UQwzCIoogsm628VbtgsixjPB7rYWZZljEcDplMJrq1V/16URQsLCywsbGB7/uEYcju7i6f/uJFfunRlIN4diVz3E/4b4+n3LFYJUkSLly4QK/XoygK1tbWOHbsGCsrK4eGj1WrVZrNJnme6wJV9T1Uq1U970SFFyFuZhJGhBA3JFUTMh9CVPHppSFETSFVpx0qwERRxHA4BGbXOYPBYDYbJM/xfZ8jR47QarUoioJut8sXT2/x05/v86WBDZjU7Jy/fiTmdWsutu2yv79Pv99nOp1Sq9XY2NhgY2ODSqWiT2Ycx6HZbOoiWjV0TbURq+/BdV25lhG3DAkjQogbynydBzwdQizLumyaquu6OI5DFEVMJhN9NaOuROI4plQqEccxo9FI73RZXl7Wi+eCIGBzZ49f/NwO//mCQVrYmBR8y0rGdxwpqJRcgiCg253tk3Ech2PHjrG2tsby8jKmaerR7uVymWq1qhfuxXGsx8vneY7neRiGoce/C3GrkDAihLghqB0w8yceaomcGtuuOI6j3+BVQaqaIRKGIUmSYBgGvu8TBAGj0YgoiqhUKnpmSFEU7O3t8am/uMi/fyyhG8/CwW3VjHfcYbBUml0R7ezsM5lMKIqCVqvF2toaq6urVKtVsizTYalSqeh2XFUcWyqViKJIb9u1bVsHEiFuJRJGhBDXvfnTELWgznGcy0KIelNX9SLT6ZTRaEQQBEwmE7IsoygKvUxOTVTNsoxOp8PCwgKe5xEEAV88s82/+1yfR4Y2YFF3cv7GSYNvWjAJggnD4dNXPL7v02g0OHLkCJ1OB9u29etU80NU4SygtwBHUaRfb6lUktkh4pYlYUQIcd0qikLviYHDs0LUacT85y3L0sPM1ICy8XhMFEVYloVlWVSrVeI4Zmdnh8lkguM4HDlyhEajAcD5rR1+8XO7fPoC+krm21Zz3nabTRKMGQ4TPQ3V8zw8z2NpaYnl5WV839dfx7IsfN/X9SqqeFa17KpQpLpl5FpG3MokjAghrkvzE1JhtqzONM1Dn1P1FqrtVZ2G9Ho9xuMxk8lEnzo4joPneXS7Xfr9PmEY0mg0WFxcxPd9hsMh//kvLvKLj6X04tlQsTvrOe+4w6Tj5Ay7+8CsTdg0TRqNBr7v62266vWpExFV96Gui9TsE7XZ17IsuZYR4ikSRoQQ150oivSgMXVyoHbIqM+pGgtADzQbDof0ej0mkwlxHOuTC/X8Cxcu6DCztramT0MeeeI0H314yB/vz7pkGk7O3zgJ97QLgmBIfxLpLp1KpYJlWTQaDVZWVqjX67oTRhWpuq57qKg2z3M9P0Rdy3ieh+M41+TPV4jrzXOaKfyhD32I48eP43ke9957Lw899NBXffwHPvAB7rrrLsrlMhsbG/zQD/3QoXteIYSAWahQXS4wK0Qtl8uHwokqBlVBJI5jBoMBm5ubbG1t6cFlrVaLZrNJpVKh3+9z7tw5wjCkVCpx9OhR6vU6g8GA3/zjR/jHvz96KojAt65kvOcbE+7yAwaDPnme6/HrtVqNUqnE+vo6J06coNPpYFmWvnKp1Wp62R2gtwCraxl1lVSpVCSICDHnqk9GPv7xj3P//ffz4Q9/mHvvvZcPfOADvPnNb+axxx5jaWnpssf/0i/9Eu95z3v4mZ/5GV7/+tfz+OOP83f/7t/FMAx+8id/8nn5JoQQN74kSfR/pKiTAzUbRBWuzk8iVach3W6XXq+nu2Sq1Sq1Wk0vr9vc3GQ6nWKaJu12m2azSRAEXNjc5mNfGPDZ/RJg0HJz/vbJlBOVdHaKkWWUy2W9Kdc0TSqVCmtra7RaLf0a1PWP2i8DT098VSciqm7EcRx9MiKEeJpRzC9qeBbuvfdeXvOa1/DBD34QmP3LuLGxwQ/8wA/wnve857LHf//3fz+PPPIIDzzwgP7c//Q//U/88R//MZ/97Gef1dccDoc0Gg0GgwH1ev1qXq4Q4jp3aZGqqqO4dIHd/II4tZyu2+0yHo+J41iPT3ddlziO6ff7DAYDsizD8zy9lG5vb48/O3vAL52y2Y9npxOvXUj5rqMpDrMC0/kaFNd1MQyDxcVF1tfXsW2boihwHAfHcfB9X9eHRFFEkiQAugXZsqzLgpQQt4pn+/59Vf9mxHHM5z73Od773vfqz5mmyX333ceDDz54xee8/vWv52Mf+xgPPfQQr33tazl16hS//du/zfd+7/dezZcWQtyE0jQ9NAZdFZrOh5P504SiKBgMBjqEqIV2zWaTer1OURSMx2N6vR5BEGDbNo1GA8/zmEwmnLtwkd/4Ssh/6VbJMajaOX/nDrirmurwoGZ/FEVBuVzGtm1OnDhBtVrVVzYqMNVqNRzHIU3TQ8v11OMMw8CyLN01I4S4sqsKI/v7+2RZxvLy8qHPLy8v8+ijj17xOX/n7/wd9vf3+ct/+S9TFAVpmvIP/+E/5B//43/8jF8niiL9X0SA7uUXQtwcVFeJOkVQJx+AflMHDhV5FkWha0LUWPdKpUKr1dIBZjgc6lbeer1OuVwmTVNOnz7Nl87v86ubNbbiGgCv6BS84zZw8ogkSfVskCAIdIdOu91mY2ODLMv0ZFTbtqnX64euZNT3Md/lA8hIdyGepRf8zPAzn/kMP/ZjP8ZP/dRPce+99/LEE0/wgz/4g/yzf/bP+NEf/dErPuf9738///Sf/tMX+qUJIa4BVeuhAod6w56vGbl09kaapmxubjIajRgOhziOw9LSkh4m1u/39ZRVy7JYWFjAMAz29/c5d/4Cv78JfzDokBYGZavge+62eWllSpomFE+FjDiOSdMU3/dxXZdjx47hui5Zlumvo+pRDMO47FRHUcPNZKS7EM/eVdWMxHGM7/t84hOf4O1vf7v+/Lve9S76/T6//uu/ftlz3vCGN/DN3/zN/PiP/7j+3Mc+9jH+wT/4B4zH4yseXV7pZGRjY0NqRoS4wam9MMChOor504VLZ29EUcTW1paelqq23XqepxfdqQ4cte12Mplw4cIFTu30+c29Fuej2anLN7QNvueOgjIxeZ7rotLJZKJrQ9rtNsvLyxRFobt2XNel2WzqepH5Ux31OtWPUpkdIsTTXpCaEdd1edWrXsUDDzygw0ie5zzwwAN8//d//xWfEwTBZYFD/dfCM+WgUqkkR5tC3ESKomA6neq6DPWGXRQFQRAcqteYH4k+mUzY29uj2+0yGo1ot9s6cKgQoqao1mo1kiRhd3eXzc0t/njP4DODFeLCxLXgb95h85pWhGFAGCaUSiV9SqOuiNbW1nQnjrqSUe26cHmNi2maFEWh/15mhwjx3Fz1Nc3999/Pu971Ll796lfz2te+lg984ANMJhPe/e53A/DOd76T9fV13v/+9wPwtre9jZ/8yZ/kla98pb6m+dEf/VHe9ra3yRGmELeA+Tfw+Wmol37+0rbd8XjM/v4+e3t7RFFEp9OhUqkA0O129emE7/skScJwOGRra4sLB0N+Z6/FqWj22DtaJu+6y6RpJ0RRoq9b1IZdNaBMdcqoFl21YVfNEZkvqlXFqSpEXdrtI4S4OlcdRt7xjnewt7fH+973Pra3t3nFK17BJz/5SV3Ueu7cuUP/Qv7Ij/wIhmHwIz/yI1y8eJHFxUXe9ra38b//7//78/ddCCGuO5deZ8x3lcwPMbu02yRNU0ajEfv7++zs7FAUhR7ZPp1OmU6nh6azqtHum5ubfKFn8bv9Vaa5hW3Cd9/p8ca1jHA6JQgivY03DENdJ9JqtWi325RKJarVKp7nUa1W9enspachtm3r+SEgRapCPB+ues7ItSBzRoS4sWRZpgd/wdNv2Jde16iBYUoURYxGI/b29tjb28NxHFqtFtVqlcFgoBfUwaz+RI2A3zoY8P/1O3xxPLtOOdaw+XsvtejYEWEYHlqwB7OAoQpdK5XKocmpatndlU5DLMsiTdMrnuYIIS73gtSMCCHE13KlvTLz23Qvva6B2bVMGIaMRiN2dnbo9Xp4nke73cb3fXq9HlEUMRgM9NeYTCYMh0OenLj8xs4ao9TCAP6bu33etJoQhxOm00R3taj23CRJqFQqNBoN3R3j+z61Wk3Xq6iBa+q/1eaX3QF6B40UqQrx/JAwIoR4Xlx6kjA/rGy+i+bS+gp1DTIcDtne3mY0GuH7Pp1OB8dxODg40CcmSZIwGo1mpyJJxqf2GjzUm52srFZM/sE9HktWQDSdhYZqtaqHkKVpimEY1Ot1arWaPnGpVqv4vq+LUee/B8uycByHOI71Kc+lRbZCiK+fhBEhxNft0usX1VWiPn/pqPf5tt0oiuj3+2xvbxOGIdVqlVarhWVZdLtdXSeSJAl7e3vkec526vPLZ1y68awI/s0nSrz9pEEWTfQAM3X9o05q1P6YZrNJs9mkVqvpjbvAoZMb4NBQM7h89okQ4vkjYUQI8XWZH2I2P+zr0uFm8ycKKqQkSUK322Vvb0/PCWk2m5imeahjZjweMxgMSHP4/YMan9lxKDBoeQb/v1f43FZJSJLZ0DLVCZMkiR4toEa3LywsUK/Xqdfr+jQEuGzgWqlUIo7jQ7UtsuBOiBeOhBEhxHM2f5owf/1yadvu/DRS9WvqRGRvb48sy6hUKjSbTfI8Z39/nzRNSZKEg4MDJpMJYW7xi2d9Tk9mP7Zev2ryvS/zsbLZOPeiKKjVZqPewzDUhau2bdPpdOh0OjQaDRqNxqFZIPM1LrZt6yFsV6ptEUK8MCSMCCGek/nThPmCzkvbducLPdWvBUHAYDDg4OBATzptNpukaUqv1wNm3TL7+/vEcUw3sfjoKZ+D2MKz4N0vtXndRnnW5pskmKaph56FYaiLVT3PY2lpiVarRafT0WPd4co1LsCh70kW3Anx4pAwIoS4avMFqWpQGBw+ZZifvzF/LTMejxmNRvT7fQzD0J0tYRgyGAywLIvBYEC/3ydJEk6NTH7hTJlpZtJyC77/Gy1uWyjpK5xyuawX5SVJQpZlFEVBo9FgZWWFTqdDs9k8dLpxpRbjLMsua0UWQrw4JIwIIa7K/B6Z+TftOI51EJkfi66ucuI4ZjKZEAQB3W5X12bU63XG4zHT6RTLsjg4OGA8HpPnOX+6Z/CfLvpkhcFGJef7v9Fmoerqtttms0kcx/o1xXGsl+gtLCywvLx86DRk/vWoaxjTNA9tD5YiVSFefBJGhBDPyjN1zMDTczmAQzUWKiioGpHxeEyv19NXKM1mk8FgoH9PVR+SJCmfumjx/+3ORrq/rJny91/mUPNnJyJqWNpkMsE0TYIgIEkSGo0GS0tL+uPSgWRpmuqhaep7kiJVIa49CSNCiK/p0o6Z+cmjatoqoPe8qOCipqRmWXYoiKghY8PhUP8eBwcHTKdTJmHMr5x1+UJ/1nnzrUsR73hJmZLrMp1OqVarGIbBZDLBMAz6/T62bbOyssLi4iLr6+u6kHXe/NWS+j7Uh0xSFeLakn/7hBBf1aXXGvOdMSqkwKx2RI1cn06n+uokz3MGgwGDwQDDMKjVapRKJYIgAGYhQV3NHIwj/v2FOmfGJiYFbz8S8qbbKnoMe6PRYDqd6iV1/X6fRqPBwsICKysrrK2tXbHzRc06yfOcPM918Lh07okQ4tqQMCKEeEbzLbqXTk7N85wgCCiKQnfNqNke6jl5nnNwcMBoNMKyLHzf1xNNsywjiiIODg4YDAZsTwp+4VyDg8igZOZ8z7GAV2/UdFAol8u6Xbfb7ZKmKcvLy7RaLY4fP06r1bosVMxfLSXJbGOvbdvSsivEdUbCiBDiip6pdReefpOfDylJkjCZTHQhahzHdLtdRqMRjuNQLpfxPE+fUERRxP7+PgcHB5wNPX7xTJlpZtC0U773+Jg7l+uHrlFUp8v29jblcpmlpSWWl5c5duzYoWV7ijqhUV/LcRwsy5KWXSGuQxJGhBCXmW/RvXSzLnCofsT3fX0iooLIdDrVxaiO4+D7PuVyWc/0GI/H7O3tcXBwwJfDBv/xnEteGKyXIv7W+oBjSx0dHiqVyqHn1Go1lpeXOXr0KEtLS1fsfFGFqmma6lHwqntH9soIcf2RMCKEOGR+l8yV3rzVtcczBRFVH6KuVFSNiLK/v8/u7i69Xp8/mi7xwOYsTLzEn/CdywOOrq/p4FCpVDAMg16vx+7uLo1Gg9XVVe666y5dyHop1cGj2n9Vjcv8FZMQ4voiYUQIAXz11l0liiIdVNTVzHQ6JQgCXNdld3dXzxKxbZt2u61Hsk8mE3Z3d9nZ2aE/mvDJ/goPH8zCxDfXenxLa8SxjWO6yLVcLmPbNhcvXqTb7dJqtTh69Ci33367HrJ2KdVGHIahLk6VAWZCXP8kjAghnnHZ3bz5oWbzQWQ6nWLbNjs7OwRBwGg0olKp0Gq19AlHt9tlZ2eH7e1tRgn8h90VzozApODba1u8qp1y7NgxLMui2Wzq+pQzZ84wGo3odDqcPHmSEydOXPGaRQWpIAj0tUypVJIBZkLcICSMCHGLe6Zld/MuHWqmajJUm+329jZRFDEYDKjX6zQaDVzXxXEc9vb2OHv2LP1+n/3E4ZcuNDgIwTMz/lrtAne2TI4c2aBUKtFsNvF9nzzPeeyxxwjDkFarxV133cXRo0evOAskz3PG47Hu7PF9X4cRadkV4sYgYUSIW9j8RNJLO2aUS4eaqb+fTCbAbGqq2jnTbDapVCr4vg/AxYsXOXfuHEEQcDos8/HzFaYpNO2E7yif4kS7wtGjR/V+mmq1yng85itf+QppmrK4uMhLXvISVldXn7FQdTgcEoYhpmnqry0DzIS4sci/sULcoi5ddnel4V/zQ80syyLPc73QDmA0GulZJPV6Hc/zqFQqFEXB6dOn2draIkkSHh7X+LXzLnkB686U+7xTrLRrnDx5Et/36XQ6+jrn9OnTpGnK0tISL3vZy1hYWLhi4al6HVmW6ULZK4UpIcT1T8KIELegZ1p2N29+qJkKA5PJRG/WHY/HZFlGmqY4joPjOFQqFbIs47HHHqPb7ZKkGb/Xa/KZrdnz7yoNeEPpHMuLHW677TZqtRrtdhvTNNnb2+PUqVNYlsXa2hove9nLrjjIDGYhaDweUxQFruvSaDRkgJkQNzAJI0LcYr5W6y4cHmpmGAZFUTAcDhmPxziOQxAEmKZ56Pep1+vEccwjjzwy20uTw3/aafOFg9nveW95l3ucbVZXVzh58iStVotGo0GWZWxubnLu3Dlc12VlZYWXvexl1Ov1K57U9Ho9faJTqVSo1WrSsivEDU7CiBC3kDAMD7XmPlNtheqsgafrMqIownVdPT9kPB6TpqnevpskCY888giTyYRJ7vCxs3XOjgoso+CN/gVud/qsrx/h5MmTLC0t4XkeYRiytbXF5uYm5XKZ1dVVXvrSl15x0Z1appemKYZh0Gg0dG2KEOLGJmFEiFtEHMf6auZrBZEsy8jznDiOdbus6qIxTZPBYKBbaBuNBlEU8dhjjzEej9mJHH7hrE8vKvDNjDeVT3PUT1lfP87x48dZWVnBsiyCIGBra4vt7W0qlQrr6+vcfffdVwwYaZpycHCg60M6nY4UqQpxE5F/m4W4BaRpeqg195neyNVQsyiK9P9OkoRKpUKSJBRFwWg0IooiPM+jXq8TRRFf+cpXGI1GnJu6fPR0mSiDjpPwHf5plisWKysb3H777SwtLZFlGePxmM3NTQ4ODqjVajqIXGnHTJIkHBwc6G27nU5HZocIcZORMCLETU7NEYFZseoz7WZRY9Sn0ylJkujTkUajQRAEZFlGv9/XQUQVq546dYrBYMD5yONnT3nEucFRL+Sv+mdpVUqsrKxw991302639WCyra0tBoMBjUaDjY1ZULnS64qiiG63S1EUOI5Dp9OR+hAhbkISRoS4ic235tq2/Yxj0ZMkYTQa6Suaoih0u6zqWun3+0ynU1zX1ZNNn3jiCYbDIRfiMj/7ZIk4NzhWmvLXqufoNOt0Oh2+4Ru+gWq1qsfEb29vEwSBDiInT568YhAJgoB+vw/MTnPa7ba07Qpxk5IwIsRNar4jxrKsK16BwCyIdLtd4jjWNRm2beM4jg4i3W6X6XSqB4uVSiWefPJJBoMBpyY2P/tkiSQ3OO5NeZN3iqXOog4ipVKJyWTCeDxmf3+fKIp0EDl27NgVW3JVCzGA7/s0Gg0JIkLcxCSMCHGTunTXzJXezKfTKf1+X1/LqMeZpqknrO7v7+ug0mq18DyPJ598ktFoxJNjm4+e8mZBpBTwJu80K4sdOp0OL3/5y7Ftm9FoRBAEdLtdoiiiVqtx/Phxjhw5csXaldFoxGg0AqBarVKv11/YPyghxDUnYUSIm1AYhmRZhmEY+L5/xRHvQRAwHo9JkoQ8z6lWq3rA2aVBJI5jOp0Onudx6tQpJpMJj/bgo6c80sLgeGnCXymdYqnTYXFxkXvuuUfPJgnDkF6vR5ZlNBoNjh07xtra2mVBRD1efe1arXbFFl8hxM1HwogQN5koinQLr+d5lxV8qk6ZIAhIkkQvxwMwDIPJZEJRFBwcHBCGIVEUsbi4iOd5bG5uMplM+OJ+xs89FUROliZ8e+kUi50Wq6urvPzlL9eL9YIgYDgcUhQF9Xqdo0ePsrKycsUgompSDMOgXq9TqVRenD8wIcQ1J2FEiJtIkiTEcQzMgsilb/rqlEMVqqraEGUymegpp5PJhCzLWFxcxPd9Ll68SL/f5wu7CT9/ukxaGJwoTfgr3mk6rSZHjx7lJS95CVEU6fkkKtjUajWOHTvG4uLiZa9pfqqqYRg0m00djoQQtwYJI0LcJC7drntpYag6rVC1JKZpYlkWhmHo9t80TRkMBjqItNttfN9nZ2eHXq/Hn+/EOoicLI359tJpGrUqR48e5a677tJBJwxDvdem0Whw9OjRKw4qy/Ncb/01TZNWq/WMHT9CiJuXhBEhbgJfq4VXBRV1haOKVC3LIkkSfVKigkhRFLTbbWq1Gjs7O+zv7/PwTsIvXBJEapUyt912G3fddRdFURBFkQ4iamT7kSNHaLfblwWRNE3p9Xo6iHQ6HVl2J8QtSsKIEDe4S1t4L73iUEFFnYyoYlbHcUiShCAI9P4ZFWiazSb1ep3d3V12d3d5eDfl5097ZIXBbaUx3+aeol6tcPvttx8KIupqxrZt6vU6Gxsb1Ov1y4KIaidWV0WtVkuCiBC3MAkjQtzAVBBR1y6XBpGiKHShqioOLYoCz/OI45jJZEKapkwmEz0uvtlsUqlU2N3dZWtriz/fy/j5UyWywuB273AQueOOOyiKgjzP9VZfy7JotVqsra1Rr9cvCxlRFNHr9cjzHMdxaLfbMt5diFuchBEhbmDzLbyXzhJRQURd0ahfc12XKIp0W28YhsRxTFEUunh0f3+fra0t/uIAHUTuKE/4VvtJapUKt912G7fffjumaWKaJnt7e0wmE0zTZHFxkbW1NWq12mVBZDqdMhgMyPOcUqlEq9WS8e5CCAkjQtyo1CI7mG3hvfRNXQWVIAj0yYlpmkRRxGg0OrR/RrXTuq5Lt9vl4sWLfLlv8dEnLLLC4M7ymG+xn6Req3Ls2DHuuOMObNvGsiy2t7f16cz6+jrLy8tUq9XLgkgQBAwGA4qioFwu02g0JIgIIQAJI0LckFSLLsyCyKXXHGEYkiSJbtVVQUTNDymKQl/ZqHBQKpXo9/tsbW3x5aHDz37FmJ2IeCO+1T1DvVrjyJEj3H333TqI7Ozs6C6c9fV1FhcXqdVqh3bNFEXBZDJhNBpRFIWMdxdCXEb+s0SIG0yaprq+o1QqXXGWiLp+URNV1b6Z3d1dsiwD0OFEtQH3ej02Nzf5Yt/iZx6bBZE7vRHf7p2lWa+xvr7OS1/6UlzXxbIs9vb2dBBZW1uj0+no0xWlKArG4/Gh8e4SRIQQl5KTESFuIPOzRBzHuWzbreqYCcNQX+E4jqOvU9RkVtu2SdMUwzCwbZter8fBwQFfGrr8zGMFeWFwd3nMt7hnqNdqLC8v85KXvATbtjFNk263q+tRlpeXabfbNJvNQy3FeZ4zGo2YTCYYhkG1WqVarUoQEUJcRsKIEDeI+RZe27Yv28Krgoo6GYHZyYlpmnqialEUuqU3yzIcx6Hb7TIYDPiLvs3PPJKTY3C3N+KN5fNUK08HEcdxsG2bwWCgg0inM1uKd6UgMhgM9KbfWq0m492FEM9IwogQNwDVGaOuXS4NImqWyPwVjuu6mKbJdDplf39fzyHJ85w0TXVImU6n/HnX1EHkJeURb/QvUPF9Op2OLlYtlUqHprPW63U6nQ6tVutQS3GWZXpmiWmaNBoNGe8uhPiqJIwIcQMIw1B3vTxTC2+apnpomaoDSdOUCxcuUBQFWZbheZ5+zHQ6ZTqd8ic7OR99tCDH4KXlEd/mX8Ave7RaLe688048z6NSqTAajXSHjud5LC0tXRZE1Dj5KIqwLItGo3FZcBJCiEtJGBHiOqfqP1QQmW+HnZ8lMj8O3nVd8jzn3Llz5HlOkiS02219qqG6bf5kt+BnH4UCg2/wR3xr+QIVv0y73ebEiRNUKhWq1aoeihZFEbZt6yAyf/UyP97dtm0ajYbsmRFCPCsSRoS4js3Xf3ied8UWXhVE1BVOuVymKAq2trYIgoA4jlleXtYj4VXtyJ/sFvzsIwUF8LLKiG/xLlCt+HqfTKvVolqtHiqKLYqCpaUl2u021WpVv4758e6O49BoNC4rrhVCiGciYUSI69TXauFVJybzVziVSgXDMNjf32dvb48kSVhYWCBNU+I4ZjgcPhVE4KOPzoLIy/1ZEKlUfKrVKmtraywtLVGr1XQtihopv7q6qhfoKfPj3V3XpdlsXvZahRDiq5GfGEJch+avXVzXveyUYX6WiBoHX6lUsCyLwWDA2bNnyfOcZrOJ67qMx2P6/T55nvPHOwW/8LhBAXxjdcy3eBfx/VkQWV9f1ztlDMPQM0KiKNItvPNzQqbTKf1+n6IoZLy7EOI5k58aQlxn1GkEoLtY5qlrEzUOXi2+c12XMAx57LHHyPMcz/NoNptEUcT+/j5JkvDlbsHHVBCpjPjW8kU8r0S5XGZtbY0jR47QbDaxLIsgCBgOh0RRxOLiIu12m1arpYPIZDKh1+vpr99utyWICCGeE/nJIcR1RtV/WJZ1WSeKupZJkoQkSXQQ8H2f6XTKo48+SpZlWJbF+vo60+mUra0toijiwijnpx81yIGX+mO+1d+kVHIpl8tsbGxw5MgRFhYWdDtwr9cjCAIajQadToeFhQUdREajEYPBAADf92m32zLMTAjxnMk1jRDXkSiKnrGFN89zfS2jtuy6rovv+8RxzJNPPslkMsGyLE6cOMFkMtFL7IYxfORxlzCDI+6Uby1fwHEcPM/TQWR5eVmHHTUkrVqtsrS0RKfT0bts1NAzgFqtdqh+RAghngsJI0JcJ1SRKcw6Z640S0S15aoprNVqlTzPOX36NL1eD8MwOH78OEmSsLu7y3A4JEwLfvorJXqRQdtOeHPlLLY5O9E4duwYR48eZWVlRdeh9Pt9xuMx5XKZlZUVFhcX9bA0NVVVbfmVqapCiOeDhBEhrgNFUeidM67rHupGUUFEnYyoFl7VWnvhwgV2dnYwDIP19XVM02R7e5vd3V3yAn7xVInzEwPfzPiu1iZumlKrtTh69CjHjx9naWlJj4fv9/sMBgMcx2F1dZWlpSUsyyJN01mwCUOZqiqEeN5JGBHiOjA/J+TSgtXpdKq7a+ZbeC3L4uLFi1y8eBGAxcVFKpUK+/v7nD9/HsMw+M3zNl/qW1hGwXd1dnDjAfVmk42NDW6//Xba7TZJkpCmKf1+n16vh2VZrKyssLy8jG3bxHGsO2osy7psD40QQny9JIwIcY3Fcazbcy89bVA1IqqWBGbXK7Zts7e3x+bmJmma0mw2abfb9Pt9Tp8+TZ7n/Jdtk/+y6wDw1xd7NOIDavU6Kysr3HnnnTSbTdI0JU1T3foLsLS0xMrKCo7jEIahnr4qw8yEEC8UCSNCXEMqaMDTG3aVKIp0G2+WZRRFoYPIcDhkc3OTMAyp1WabdSeTCadPnyZJEr6wn/NbW7N6jjd2RhzJtnHLZb34rl6vk+e5rkE5ODggTVMWFhZYW1vDdV2CINATXEulErVaTYKIEOIFIWFEiGtkvk7Etm0cx9G/liSJLihN01TPDbFtm+l0yubmpi4yXVxcJAxDzpw5QxAEnOrGfGKzTYHBKxtTXm5tURQmrVaLkydP6s6YNE3JskzPIGm1WmxsbFAqlQiCQI+PL5fLVCoVCSJCiBeMzBkR4hpRVy+maR6aJ6Laa1V3TZ7nlEolLMsiyzK2t7cZj8fYts3CwgIA58+fp9frsT0I+eWtNnFucLIS88baDkUxm8R69OhRVldX9TZfgJ2dHX26cvToUTzP00EkTVN835cgIoR4wUkYEeIaUEPL4HAb75Vmidi2jWEYeufMeDwmTVNarRau67K5ucnu7i79ScjHt9oME4OlUsZ/u9QljSPqT9WJrK+v4/u+Hoq2vb1NEAS6xVcFEfX1fd/H930JIkKIF5yEESFeZHmeH6oTUZt4i6LQHTPzy+9M08RxHEajEf1+n+l0Sq1Wo1wus7+/z9bWFoPhiF/dbrEVWlTtnL9zZEA6HeH7Pq3WrI233W6TZRm2bbO1tcVkMsF1XR1EwjDUrcO+71Mulw9dHQkhxAtFwogQLzL1hm9Z1qFTB3Vto9p8i6IAZnNHptMp3W6XyWRCuVymXq8TBAFnzpyh3x/w6W6Lx8cOjlHw3x+bYE57uK5Lo9FgY2ODlZUViqLAcRx9zWNZFseOHaNcLuule+rKyPM8CSJCiBeNFLAK8SJSnTGGYRyqE1HFqiqQZFmGaZr4vk9RFBwcHDAejzHNWSFqkiQ8+eSTjEYjHur7fG7gY1DwN48G1OIDsCx832dlZYVjx45hWRaGYXDhwgVGoxGGYbCxsaGvbaIowrZtbNuWICKEeNFJGBHiRaLqQGBWJ6LaeFUYUAWraZpiGAa+72NZFru7u0wmE+I4ZnFxUQeRwWDAX/QsHug2AHjLWsR6sY9pWViWxdLSErfffjvlcpksyzh79ixhGGIYBqurq3qCq5ohYpom5XL50PRXIYR4Mcg1jRAvAlUPAuA4jn7DV5+f/2ue51SrVUqlEsPhkH6/z3A4pN1uA7POmW63y5O9lN886FBg8M0LCa+qDHQh7NLSEnfccQetVos0TTl37pz++qurqzQaDUzTJI5jXNeVICKEuKbkJ48QL4L5nTLzo9RVAAmCgKIoSJJEDxcLw5C9vT0GgwHNZhPDMNje3mZnZ4etYcSvHayT5AZ31lLe1OkTTWenHq1WixMnTrC+vs54PObixYu6HmR1dZVKpYJt2yRJomtWJIgIIa4l+ekjxAtMXb3A7E1fnV6o+hHVSptlGZ7n4bouRVFw8eJFRqORDgpq50w/iPi17iqj1GC1nPE3j4zIohjTNKlUKhw/fpzbbruNwWDA9vY2YRhiWRarq6v690/TVNeFSBARQlxrz+ma5kMf+hDHjx/H8zzuvfdeHnrooa/6+H6/z/d93/exurpKqVTizjvv5Ld/+7ef0wsW4kZyaRuvqhNRE1bTNNW/bts2lUoF0zR16y3MwkK/3+fs2bNMo5jf2F9iJ7SoOznfe2xMHs02+jqOw9GjR3npS1/KYDBgZ2eH6XSK67qsr69TKpXwPI8kSXT4kCAihLgeXPVPoY9//OPcf//9fPjDH+bee+/lAx/4AG9+85t57LHHWFpauuzxcRzzpje9iaWlJT7xiU+wvr7O2bNnaTabz8frF+K6NV8nYtu2vhJRc0SKomAymWCaJlmW0WjMClF3d3cZDoeEYUiz2WQwGHD27FmKAv5zt82TEwfXLHjX8YBSNiXJc2zb5siRI3zDN3wDg8GAg4MDPdBsYWEBx3HwPI8oivTrUAWyQghxrRmFGmbwLN1777285jWv4YMf/CAw+8G6sbHBD/zAD/Ce97znssd/+MMf5sd//Md59NFHn3O74HA4pNFoMBgMqNfrz+n3EOLFFoYhSZJgGAaVSgXDMHR9SJ7njEYjYNZNU61WsSyLyWTCzs4O/X6fWq3GdDrl3LlzJEnCH+x5PLBfxaDge09MOeGOiePZ9czy8jLf9E3fRFEU9Ho9JpMJlUqFTqdzKIio05lyuSxBRAjxgnu2799XdU0TxzGf+9znuO+++57+DUyT++67jwcffPCKz/mN3/gNXve61/F93/d9LC8v87KXvYwf+7EfI8uyZ/w6URQxHA4PfQhxI0nTVI97n68TUZNV5ztnKpUKlmXp7bmqTmQ8HnPu3DmyLOOLQ5cH9metuG9dnXKbN9uma9s2zWaTu+++mzRN9YTWer1Oq9XCcRxKpZI+EZEgIoS4Hl1VGNnf3yfLMpaXlw99fnl5me3t7Ss+59SpU3ziE58gyzJ++7d/mx/90R/lJ37iJ/jn//yfP+PXef/730+j0dAfGxsbV/Myhbim5rfxuq6r3/jnZ4nEcYxhGLiui+M4RFFEr9djMBjoltsLFy6QZRnnpza/ujn7L4rXL0S8phno2SC+73PHHXdgWRbj8ZggCKjVatRqNV2sqgaaWZYlQUQIcV16weeM5HnO0tIS//bf/lte9apX8Y53vIP/9X/9X/nwhz/8jM9573vfy2Aw0B/nz59/oV+mEM8bdephWZZu41UhJMsygiDANE1s26ZcLhOGIePxmG63q09ULl68CMAwd/ml83XSAu6uxtzXGepw4bouR48exfd9vVemXq9TqVR0Yao6PXEcR4KIEOK6dVUFrAsLC1iWxc7OzqHP7+zssLKycsXnqJXl8z8EX/KSl7C9va0HLl2qVCodmsUgxI1CBY75ce+XFqxalqV3wERRpK9ngiDQ23SzLCPG5ufPVBknsOal/DerffKnFt25rsvi4iLNZlNPcK1UKniepyerJkmCaZq4ritBRAhxXbuqkxHXdXnVq17FAw88oD+X5zkPPPAAr3vd6674nL/0l/4STzzxBHme6889/vjjrK6uympycVNRoQCebuOdn6w6mUwwDEMXtKpdNAcHBwyHQ4qiYG9vjzzPyTH5pXNVdqbQcHL+5mqXkmXo5XqVSoXFxUUMwyCOYyqVCqVSiXq9TlEUeq6J53kSRIQQ172rvqa5//77+chHPsLP/dzP8cgjj/CP/tE/YjKZ8O53vxuAd77znbz3ve/Vj/9H/+gf0e12+cEf/EEef/xxfuu3fosf+7Ef4/u+7/uev+9CiGtsvk7EcRzdOTZfsKoCeaVS0fNFRqMRvV6PLMsYDodPLciz+K3dBl8ZQMks+JvL+zRLxqEaEzXArCgKPM/Dtm0ajYZesgez1l0JIkKIG8FVzxl5xzvewd7eHu973/vY3t7mFa94BZ/85Cd1Ueu5c+f0YCeAjY0NPvWpT/FDP/RDfOM3fiPr6+v84A/+ID/8wz/8/H0XQlxjKnTMj3tXQ83URl61/E6FkyAI2NraIo5jgiDQbbr/tV/jj3ZyTArevrjPWgW9yC7Pc5aXl6nVapimiWEYuqNmPoiouhEJIkKIG8FVzxm5FmTOiLieJUmiT0XUIDFVqJplmb6e8TwPx3EYj2fzQc6dO0cQBIxGI9I0xXVdToc+H/xCRgG8ZbHPa1qR3vA7nU5ZW1tjfX0d13XJnxp2poalZVlGURRUKhXdLiyEENfSCzJnRAhx2KXj3i3LumwTrzq9KJVKjMdjsixje3ub6XTKYDAgTVM8z2NqVfnpL8+CyCtrY17VmOrfM4oiGo0Ga2tr+grINE3q9TqGYZCmKXme4/u+BBEhxA1HwogQX4f5Nl5VkK3qQ9SvmaZJtVplNBqR5zm7u7v0+30dREqlErZX4cNfyglSWHND3tjY11cx0+mUcrnMbbfdhuu6epJrs9nEsiziOKYoCsrlsp7kKoQQNxIJI0I8R1EUkec5hmFQLpcB9AZe9VfDMKhWq0wmE/I8p9frsb+/z2g0Yjqd4vs+tVqNX/yKwflhRsVMeVt7h2a9RpZlJEmC4zicPHlSX9cURUGj0cC2bX095Hke9XpdgogQ4oYkYUSI50AVpsIsCBiGQZIkul13vqNFfX4ymbC9vc1kMmE0GulJqX+w6/LgxQiTgu+oX2SpVsJ1XR1mbr/9dqrVKoZhkGUZtVqNUqlEEAS6w6bRaEgQEULcsCSMCHGVLh33btu2Pg1J05Q0TXXLrWEYTKdT4jjm/PnzBEHAwcEBzWaTarXKZlbjF74wW5j3huoOJ2uFnkECcOzYMV30VRQF1WoV3/cZj8d6imur1ZIgIoS4oUkYEeIqqWmqqo1XFaqqoWdFUeh5IEEQkCQJZ8+eJQgCdnZ2aDab+L6P3Vjixz+7TwHc6fR4ZXXM0tISqsFtbW2NhYUFPZ/E932q1SrD4RDLsrAsi3a7LUFECHHDkzAixFVIkkRPN1V1IqpgVQURy7KwbVsvxtvc3NRXNLVajUqlwvLaEf757+0wjAsWzID7mnssLy8BsyughYUFlpeXdftuuVymVqsxGo2wbRvTNGm329j2VY8KEkKI646EESGepaIoLhv3rupD1LWNauNVVzm7u7v0ej22t7cplUpUKhU2Njb4fx7c48legkfCW5ubLHVa+L7PZDKhVquxsrJClmW6S6fRaDAajfSgs3a7rVt8hRDiRidhRIhnSV3PqICgiljDMMQwDPI81wFBzRDZ3d3l4OBAzwQ5cuQIn3wy4IFTEwwK3ly/yErNZXl5mX6/T61WY3V1VYca13VpNpuMx2MMw9AnIrLXSQhxM5EwIsSzMH8943me3sSrZnyoCaowu2YZj8dsbW3R7XZJkoRGo8Hq6ioXwhL/+o/2AHhtaZMTfsKxY8c4ODigUqmwvr5OuVzW23bb7TZBEOgQ1Gw2ZaO1EOKmI2FEiK/h0usZ1SGjAkpRFDiOo1tx1YnI/v4+QRDoaxer2uZ9nzpLVsAJq8urq0PuvPNOBoMBtm2zuLhIrTabL+I4Du12WxfG2rZNvV7XdSpCCHEzkTAixNdw6fWMauGNogjTNDFNUweGXq/HwcEBFy9eZDKZzIpVl5dpthf40d85Qy/MaRoBb2rscOLEcT2ldWlpieXlZeI4xrIsarWaXrDnOI4ufBVCiJuRhBEhvopLr2dUQJhOp3q+iAoio9GI4XDI+fPnmU6nOI7D4uIiy8vLfPihfb68F+GS8lf9Mxw/MtsxE4YhzWaTtbU1fbJSKpWwbZs4jnFdVy++E0KIm5WEESGewaXXM6p9VwUNNWVV7Y8ZjUacOXOG6XQKwPLyMqurq/z+uYhf/1IXKPi20hlefnyFSqWiO2c2NjZ0K7BhGNRqNcIwxHVdPS7eMIxr9ccghBAvOAkjQjyD+esZtQdGXc1kWUZRFBRFQZZljEYjzp8/z2QyIcsyOp0OS0tLbIUO//L3zgPwSnuTe49WqVQqTKezjbxHjhyh1WoxGAz0zpkkSSSICCFuKRJGhLiC+euZUqmkC1bVaYi6nsnznCAIuHDhAv1+nyiKaLVaLC4uYnh1fvj/fYI0h6NWn+/YgEqloieqHjlyRLf0ZllGpVLRI97VTBLTlH9FhRA3P/lJJ8QlLr2eUSEkiiIsy9J/XxQFeZ5z/vx5BoMBQRDQaDRot9t0Fhb5kd85xUFY0DAjvnOpR71ew7Is0jRlbW2NlZUV4jhmNBrheR6+7+tTmHq9LtNVhRC3DAkjQlxi/npGbeNVBatqU68a065ORAaDAZVKhWazyerqKj/1h5t8cTfGMXL+enOT1YXZMrssy2i326yvr1Mqldje3sZxHH0KYpomjUZDpqsKIW4pEkaEmDN/PaO6XVRhaZIkGIah54Bsbm7qIFIul2k2mywvL/O7j/f41S8PAbivtsVL1pr6FKVSqXDs2DGazSZbW1uYponneXq8fLVaxfO8a/lHIIQQLzoJI0I8Zf56xnVdveguTVOyLMM0TYqiwLZtPeZ9f38fy7KoVqssLCxwqhvxr/7rbMLqqytd7l2bjY0vlUo4jsOxY8fodDrs7+/rItZqtapDSa1Wu5Z/BEIIcU1IGBHiKfPXM6pAVW3iVUHEMAz6/T67u7t0u13yPKdardLpdEgNh3/+mV3iHI57U75jbXaSUi6XybKMo0ePsrKyQhiG7O/v644Z13X1qHchhLgVSRgRgsPXM+oqRoUTwzD0krrpdMrW1haDwYAoiqhWq7Mx7X6F939mm71pQcNO+W9XBzi2RalUIk1Tjhw5wsrKCqZpsrW1hW3beJ5HpVLBMAw6nY50zgghblny00/c8tQwM5gNMFOj3ueHmlmWRRRFbG1tMRqNdAdMpVKhVqvx8392wJ/vpdhGwXcvd6m5pi54XVlZYWVlhVqtxvb2NkmS6CsZ0zRpNpvSOSOEuKVJGBG3PHUVYxgGeZ6TpilhGAJg2zaO4xBFEfv7+wwGA/r9PuVyGd/3qVQq/OG5Cb/2+Ozxb1sZcqSKHhW/sLDA6uoqi4uL7O/vMx6PcV2XarWK4zj4vi/L74QQtzwJI+KWNr95V3W8BEGAaZo4joPjOKRpymAwYG9vj36/j2EYOlBsTwr+9ecnAHxza8o9jZiiKADwfZ8jR46wurpKv9/Xz/V9n1KpJAWrQgjxFAkj4pY1fz2jAsR4PAbQG3oBRqMROzs7DIdD0jTFdd1Z228KP/knY8IMTlRS/uryhKIocF0X0zQ5fvw4y8vLRFHEaDQiiiIqlYo+DVFdNEIIcauTn4TilqUKVNM0xTRNJpOJ7pzxfR/TNBmNRuzt7dHr9QjDEMdxsG0bt1Ti334hZGtSUHdy/tbGhCSa/XqSJBw/fpyFhQVc12UwGDAej6nVanrAWblc1mFHCCFudRJGxC1JjXRXQSSKIuI4xjRN3eEShiEHBwfs7e0RhqGekOq6Lv/5PHx+N8MyCv774wFEs4LWOI45evQoi4uLtNttut0uo9FIX/vU63XK5bIMNhNCiDkSRsQtR13P5Hmua0WCIMC2bXzf18WnBwcH7OzsEASBHnpmmiZnph7/6SsJAG/fiGkXI33Ksbi4yPLysl6AN51OSdOUSqWil9+paatCCCFm5CeiuOWEYUie58RxjGVZ9Pv92dWL6+oBZd1ul52dHQaDAVmWYVnWrK7Eq/HTX04pgHsXM17uD3Unjud5HD16lKWlJcIwJAgCRqMR1WoVy7LodDq6KFYIIcTTJIyIW4q6ngnDENd1GQ6HOI6DaZrU63XiOCYIAnZ3d9nf39cdNnmeU/Yr/MLjBoOoYKVc8NaVCVEU4TgOhmFw8uRJWq0WjuMwHA6ZTCZ4nodpmiwsLGBZllzPCCHEFUgYEbcMdT2jTkQmk1lLrmEYtFotHVS2t7fZ3d3VNSWGYWBZFn/cK/OFvRTHhL91NCAYDXBdlyzLOH78OK1Wi0ajwWAwYDqdArMdN/V6nVKpRLlcxjCMa/lHIIQQ1yUJI+KWEYbhocV3KpQ0m09v1d3d3WVnZ4fpdEqe51iWRRzHTEod/sNjMQBvXQsphQc4jkOe56ytrbGwsEC73WYymTCdTvWU1VKpRK1W0/tnhBBCXE7CiLglqOFmqismCAI9uMx1XZIkYTAYsLW1Ra/XA2YnKWEYUmst8pEvZaQFvKyR8hLnQNeJtFotVldXabfbFEXBcDjUV0CO49BqtbBtm1KpdI3/BIQQ4volYUTc9NT1zHQ6xTRNxuMxjuPgeR7VapUwDInjmO3tbfb29jBNkziOyfOcWq3Gb1wssT3OaLoFf22xRxhO8TyPcrnM+vo6rVYL3/fpdrv6Wmc+6Mi4dyGE+OokjIibXhiG+tpFDS5zHIdGo6EX4m1vb7O1tUVRFERRpGtBTmUdfv9chAG8fblPMhlQr9cBWFtbo9Fo0Gq16Pf7ZFlGkiRUKhU8z9Nj36WNVwghvjr5KSluanEcE4YhSZKQJLPZII7j0Gw2yfOcJEnodrtcuHCBIAjI8xzbtmdbeReO8DNfCAB449KUhaJPqVQiyzJ9IrK8vEwYhoRhSBRF+L6PZVn4vo/rujJlVQghngXZWy5uWnmeM51OmU6nZFlGURS6oNS2bYIgYDqdcu7cOV0nkqYpeZ7TXljiQ1/KmKYFJ6oZr3R3Zu29T23rVUHEMAwGgwFpmuI4jg4gnudJG68QQjxLcjIiblrT6ZQgCEiSRC+wU3thVJ3IhQsX2NvbA9B7aVzX5TPdGk90E3y74K3tPZI4otlsYlmWLlhVdSJqDkm1WsUwDD3uXa5nhBDi2ZGfluKmFMcx4/GYMAyB2dVMuVymUqnoK5u9vT0uXLigr2/yPCfLMgblNX79sdkMkrcu9jHDAY1GA4DV1VVqtRqLi4uMRiPSNCWKIhqNBmmaUi6XKZfLMmVVCCGugoQRcdPJ85zhcKgHj9m2rYNIURSEYchwOOTMmTNMp1MMw9B7avz2Cv/6z2ZB5LWtkPViX3fOdDod6vU6a2trhGHIZDIhTVPq9TpZlunOGbmeEUKIqyNhRNx0giDQQUEtvyuXy9i2rTtrzp49y3A41APQAKq1Or/4pEk/zFnxcl7nb2NZFq1Wi1KpRLvdpt1uY1kWo9GIPM9xHEcPM1MdNDJlVQghro6EEXFTieOY4XBIFEVYlkW5XKZUKlEqlXRnjRr3rvbKFEWB4zg8HLT4/GaIYxa8bWGPIonodDp6t0ylUqHT6TAYDMjzHEDvsymXy7owVgghxNWRMCJuGkVRMBqN9NWL7/v6iiXLMqbTKf1+X1/POI5DkiS4rkvoL/GxvxgC8OalMU6wz8LCArZts7CwgOd5rK2tMRqNdEFss9lkMpnguq6eKSKEEOLqSRgRN40gCBiPx6RpSqlU0rthDMNgOp0yHo85deoUQRBQFIUOIqVKnX/95xFpDt9QjzmeXqTT6eB5HpVKhVKpxOLiImma6jqUarVKkiQYhqEfJ4QQ4rmRMCJuClmWMRwOieMY27Z1iHAchzAMCYKAM2fO0O12CcMQ27Z1YetvXCyzOUxoOjnfXtvG80p6xHutVqNer1Ov1xmPx8CsINZxHKIowvM86vW6LMETQoivg4QRcVMYDAZEUUSe5/rKpFQqkSQJk8mE3d1dtre3CYKAer2uB5g9kbZ54NQYg4K3LR5gJFOWlpawLEvvlllZWdEFq0VR0Gg0mEwmlEoluZ4RQojngYQRccNTw81UIam6nimKgslkwmAw4OzZswRBoAeeVatV0nKLf/u5AQDf2pnQSrusrq7q7hvTNFlZWSGKIl0nooKIZVm4rkutVrvG370QQtz4JIyIG1qe5/pUxLZtHUQcx9E1JKdPn6bf7+urlEqlglPy+DdfiJmmBcf8hFe4OzSbTWq1GqZpUiqVdEuvWrJXKpX0tFW130amrAohxNdPfpKKG5pq483zXO+DUW28o9GIc+fOcXBwoItV1cbe393xeXw/omzl/LXmDn7ZY2lpCdM08X0fx3Fot9u62NWyLDzPIwxDXNel0WjIEjwhhHieyFAEccOKoojJZEIcx7pGxHVdvbxud3eXra0tgiCgVquR5zmu67JttPmPX+4C8JbOAS3PYGFhQdd+FEXB8vKyHhs/fz2jJrn6vn8tv3UhhLipyMmIuCEVRcFgMCCOY12/oQLJeDzWdSLj8RjXdXVhq+k3+X8e6gPwTbUR39DIqFarNJtN8jzHtm0ajYbulimKgmq1ShAE+tRF7akRQgjx/JAwIm5Io9GIKIoumymSJAmDwYDTp08zGMyKU9XgM7dU4qe/mNAPc5bchL/SGepumTAM8TwP27bpdDqEYah/7zRN9RVPq9W6xt+5EELcfCSMiBtOkiSMx2N9PeO6Lq7rYts2w+GQixcvsrOzQ5ZlekS767r8Uc/n81tTHKPgOxf2adWrdDodAEqlElmW6eFmSZLgOA6maWLbtt5RIwWrQgjx/JOfrOKG0+/3SZIE0zRxHEefjIzHY/b29jh37hxJkujleK7r0qXGz/1ZH4A3Ng+4Y7mm98mkaYphGFSrVV2kalkWtm1jGIYOIo7jXNtvXAghblISRsQNZTQaEcex7o5RtSJZltHr9XjiiSeYTCYYhqEno9rlGv/XnwzJCrjLn/CtR5xD3TKu62KaJouLi/q5tm3rsFOv1/E871p/60IIcdOSMCJuGGmaHrqeUe22juPoBXjdbpc0TWm1WuR5TqVS4WOPxGyPMxp2xnetTqjXa9RqNYqiwLZt0jRlaWmJKIoAsCwL0zRxXZdyuUy1Wr3G37kQQtzcJIyIG0a/3ydN00NzP0qlEkEQsLOzw/nz54njmFqtpq9v/nTf5PfPTjEp+K6lLieOrOihZnEcYxgGlUoF13WJogjTNDEMQ3fmSOeMEEK88CSMiBuCGveeJIk+DVFzQfb393niiSeI4xjHcajVahiGwbDw+MifjQB4Q3vEa07MilVrtRphGOI4DkVRsLCwoFuAi6LQQ88ajYYUrAohxItAftKK616e5wyHQ5Ik0UWl6uSi1+vx5JNP0u12ieOYpaUlsizD8yv8X38yJMrgmBfx1hOODhl5nmNZFnEc0+l0dN1IkiRUq1VdJ2LbMhNQCCFeDBJGxHWv3++TZZneCaO6Z8Iw5MKFC1y8eJEkSVhYWNCB4uOPhJwdFpTNnL91ImJpcYEsy/TckDzPdfCwLIskSajX6/pkRQpWhRDixSNhRFzXptMpYRjqjbzqVMSyLDY3Nzl9+jRhGOL7vu6seaRv8MkzCQDfvTHljvVFsizTG3vzPAeg0WiQJAl5nlMulymVSlSrVcrl8rX8loUQ4pYjYURct9T1TJZl2LZ9aBne/v4+p0+fptfrYZom9XqdoiiILY9//WcTAO5tTXnDbU08z8M0TbIs039tt9t6lLxlWXrGiO/7GIZxjb9zIYS4tUgYEdctFUTSNNXbdj3PI45jzpw5o69n1F6ZSrXKh/8sYJwYLJcS/tZLZuFChQ41WbVSqei2XrUEz/M8qtUqlmVd629bCCFuORJGxHUpjmPdQVMulymKQnfRnDlzhnPnzhEEAdVqFcMwKJfLfPJ0wpe7BY5R8K47ctqNmm4DTtOUOI6xbZtSqYTjOHpcvDoRkQmrQghxbTynMPKhD32I48eP43ke9957Lw899NCzet4v//IvYxgGb3/725/LlxW3iKIo6Pf7OkgURaG7Z7a3tzlz5gy9Xo9KpaJnjVyY2vzHJ1IA3ro25a7Vhp7MCuiR7yrQFEVBpVKhUqlQLpelYFUIIa6hqw4jH//4x7n//vv5J//kn/D5z3+ee+65hze/+c3s7u5+1eedOXOG//l//p95wxve8JxfrLg1jEYj0jTVW3PVkLI0TTl16hRbW1vAbBuvZVlkpsO/+UJMXsA31CLedNKnUqmQpimmaerrGbVMzzRNLMuiXq/j+74UrAohxDV21WHkJ3/yJ/n7f//v8+53v5uXvvSlfPjDH8b3fX7mZ37mGZ+TZRnf8z3fwz/9p/+UkydPfl0vWNzckiRhMpkcasNV1yqPP/44m5ubTKdTKpUKlmXhOC6//BXYmxY0nYx3nMxotVr6aibPcz1HxLZtfRXTarXwfV8XtwohhLh2ruqncBzHfO5zn+O+++57+jcwTe677z4efPDBZ3ze//a//W8sLS3x9/7e33tWXyeKIobD4aEPcfObv55RHS2qaHVzc5Pt7W329/ep1WpYloXruvz5wOXBrRSTgu9eH3N0dVEXrQJMJhM9KM33fWzbpl6vU6vVKJVKMthMCCGuA1cVRvb398myjOXl5UOfX15eZnt7+4rP+exnP8tP//RP85GPfORZf533v//9NBoN/bGxsXE1L1PcoCaTCUmSkKapHs3ueR5RFHH27FnOnz+P7/vA7Iqmlzr8/JdnoePbFgJedaxFo9EgiiKyLCNJEt0WPD9Cvt1u65kkQgghrr0X9Hx6NBrxvd/7vXzkIx9hYWHhWT/vve99L4PBQH+cP3/+BXyV4nqQpimj0UgXq6prGtu2eeyxx9ja2tJzQsrlMjkmH33cmI17L8e86UjB0tISSZLowWbT6RTXdQH0qcji4qI+bRFCCHF9uKoz6oWFBSzLYmdn59Dnd3Z2WFlZuezxTz75JGfOnOFtb3ub/pyafqneZG677bbLnqc6J8StoSgKBoMBRVHoKxpVtHru3DkODg44ODigWq3quo9P7fqc7keUrZzvXh9z4vgdGIZBGIZ6NollWeR5TqvVwnEcfSLieZ4MNhNCiOvIVZ2MuK7Lq171Kh544AH9uTzPeeCBB3jd61532ePvvvtu/uIv/oKHH35Yf3znd34nb3zjG3n44Yfl+kUAs5HvURSR57kuWvU8j/F4zObmJufOnaNarc4W4HkeZyKf3zkVAfDXF3vcfWwF3/d14av6UK28lUoF3/f1TBEZbCaEENeXq67eu//++3nXu97Fq1/9al772tfygQ98gMlkwrvf/W4A3vnOd7K+vs773/9+PM/jZS972aHnN5tNgMs+L25NWZbpAmU1iExdrZw+fZrNzU1s2yaKImq1GpHp8TNfns0OeU1jwmvWPJaXlxkMBnqwGcxCsmEYLC4uHrqekcFmQghx/bnqMPKOd7yDvb093ve+97G9vc0rXvEKPvnJT+qi1nPnzkmrpHjWhsOhvrpTM0Fs2+bcuXPs7+8zHA51a2/JK/PR0y7DKGbJTbhvYcQdd3wTQRDoGhHHcYjjGNM0WVpaAmBpaQnLsqRORAghrlNGURTFtX4RX8twOKTRaDAYDKjX69f65YjnyXQ6pdfrURQFruuSJAme5zEcDjl16hRPPPEErusynU5ZWlriT8dN/v0jIbZR8HdXt3nDPXdQqVQYj8d6lojaQ1Ov16lUKiwuLlKr1fB9X0KyEEK8yJ7t+7f8dBbXRJZljEYjAL2wzrIsoihia2uLCxcuYFkW4/GYVqtFlxq/8mgIwF9pdXnZ0QUajQbj8ZjpdIplWSRJok9AqtUq9Xpdb+OVICKEENcv+Qktrgk18h1mYUT9793dXXZ3dwnDkCiKKJfLOH6Nn36kICvgLn/C65fh+PHjjEYjkiTR9SEAlmXRbDbxPE93z8hgMyGEuL5JGBEvujAMmU6nGIZBqVTSVyvdbpednR329/cpioI8z+l0Fvj1C2W2xykNO+WvLfS5/fbbCIKAMAxJ01QXvhqGQbvdxvM8XbAqLeJCCHH9kzAiXlRpmjIcDimKQl/PAIzHY7rdLnt7e6RpShAELCws8GhU5w/OTTEoeGt7j9uPrWPbtm7jNU1TzydR+2ZksJkQQtxYJIyIF01RFPp6xrZtbNsmjmOiKKLX67G3t8dkMiEIglmhU22Jj/7FFIDX13u8bMWn3W4zHo8xDEN3zaRpqutD2u02juNQLpdlsJkQQtwgJIyIF426WlHTVeM4JkkShsMh+/v77O/vE8cxjuOwun6Ej3wxJUwLNkpTvmUhZG1tjTiOKYqCKIqoVqsEQUC1WqXRaNBqtSiXy1KwKoQQNxj5iS1eFHEcMx6P9fK7LMuI45jJZEK/36fb7eqrl5WVFT65WeLJXkLZzHjbQpeV5dmskDRNybIMx3EYDoe6a6bRaFCtVmWwmRBC3ICkzUC84PI8ZzQaHZquGoYhk8mE4XBIt9ul2+1SFAWNRoNuaYXfeHwAwHe09lhvV2g2m6Rpimma5HmuZ5LU63X9obbzCiGEuLHIyYh4wU0mE6IowrZtXNfVQSQIAnq9Hjs7O/q0Y+no7XzwT2fj4V/h97m7nrK4uAjMxsVPp1PyPNctvM1mk1arhW3blMvla/ltCiGEeI4kjIgXVBRFBEGg23jV30+nU0ajEfv7+0ynU4qi4MTJk/ybP58yjAoW7Yhv7wx0HYjruozHY7IsI0kSWq0WzWaTRqMhBatCCHGDkzAiXjBpmuoA4Xme3h+jxrcfHBzQ7XYxDIPV1VX+677HF3YTbCPnLc1tGtUKKysruK5LFEXEcUyWZTqE1Go1yuUypVJJNvEKIcQNTMKIeEEURUEQBERRhOu6WJbFZDJhNBrpVt6trS2KoqBSqZA21vnlLwcAfEt1l+VywfLyMr7vk+c5/X6fPM+pVCq0220qlQq1Wg3HcXQdihBCiBuThBHxggjDkCAIdBvvZDJhMpnorpqtrS2yLMOyLDZO3MH/9dCQrIA7SiNeWuqzvLzM4uIihmGwu7tLnue4rsvCwgLlcplms4llWTJhVQghbgISRsTzTrXs5nmO7/u6YFWNgd/c3GQ0GlEUBcePn+DnvhSyGxTUrZTXu2dYWlpkbW0Nz/PY2toiSRJs22ZpaUnvnFEFq1InIoQQNz4JI+J5lec5QRAQx7GuE1GdM1EUsbe3x97eHqZpsri4yBdGHn+4mWJS8G2lUyy36qyvr7OyssL58+f1Rl51IqKW38lgMyGEuHnIT3PxvFJTVm3b1oPJgiAgSRJGoxHnzp3Dtm1838dqrvLzX4oA+KbSNidqsL6+zu23386pU6f0NU+73aZWq1Gr1fB9n1KpJJt4hRDiJiJhRDxv1DWMmrI6Ho8Jw5AkSYjjmK985SvkeQ7AkaPH+Kk/mxJlsO5MeFX5gKWlJe655x62trZ0fUm73abZbFIul2k0GnpWiRBCiJuHhBHxvEjTVAcPz/N03YgqZH3yySeZTmdL744fP84nvpJydlRQNjO+xTnN4kKHl7/85YzHY/b39wmCQJ+IOI7DwsKCTFgVQoiblIQR8XUrioLpdMp0OsV1XQzDYDgcEscxYRiyt7fHzs4OhmGwvLzMl4cunz6XAfCG0lnWWj533nkntVqNs2fPMp1OaTab1Ot1PM9jcXFRClaFEOImJmFEfN3CMCQMQwzDwHEcRqMRcRzrKatf+cpXsG2barVKXlngZ76UAHBPaY87KxEnTpzg6NGjPPHEE6RpimVZ1Ot1fN+n2WzieZ4UrAohxE1MfrqLr4s6/UiSRI97n06n+nrm0UcfBcC2bVaOHOWn/jwmzGDdnXKvt82RI0e48847uXDhAnEcE0URCwsLeJ5HpVKhXq/juq4UrAohxE1Mwoh4zvI816cinueRpinD4ZAsy5hOp5w5c4YgCHAch42No3zssYwL44KKmfLG0mlWlpe444476PV6jEYjJpOJbt2tVCp0Oh1s25bBZkIIcZOTMCKeE1UnEoYhlmXpeSJpmjKZTNja2mJnZwfTNFleXuaPuy5/tJVjUHCff5b1dpWTJ09imiaDwYA4jvXVTKVSOVQnIoQQ4uYmYUQ8J1EUEUWRHumurljCMGRnZ4fz58+T5zkLCwv0rQa//FgKwOv9HU5Uc06cOEGj0aDX6+kJq57n4fu+Ph2RglUhhLg1SBgRVy1NUx0+XNcljmOCICDLMra2tjh37hyTyYROp4Nba/NvvpiRFXBbacQr/T7Hjh1jcXFRn6TkeU6pVNI1IqqLRjbxCiHErUHCiLgqRVHoOhHbtkmShCiKyPNcj3rv9/u0Wi2qtTofe9KmG0LLirmvusX6+horKys60GRZpk9Bms0m7XYbx3FwHOdaf6tCCCFeJNKiIK5KGIbEcawnqcZxTBzH9Ho9tra22Nra0icc/6Vb5UsHObaR89eam6wsNDl69Oih38c0TTzPo9FoHNo7I4QQ4tYhJyPiWVPBI4oiLMvSY96jKOLcuXO6YLXRaLBJm988Mwss9zV2OVq3OXnyJIZhEEURaZpiGAamaVKtVul0OpRKJQkiQghxC5IwIp6VPM91gapt27pGJE1Tzpw5w8HBAdPplMXFRdJSg599ZBZE7vH7vLwW6iCSpilZlumTlWq1ysLCApVKhXK5LIPNhBDiFiQ/+cXXpNp44zgGYDqd6hOOra0ttre3dZ2I4/l89Cs2QQorTsgbWz2OHTuml9slSUKSJOR5ju/7dDodPWVVClaFEOLWJDUj4mtS1yqqVqQoCrIso9vtcvbsWQaDAZVKhVqtxv+7Web8xMAzMv7GSo+1xSWq1SqAvtYxDINSqcTCwoKuE5FNvEIIceuSkxHxVaVpSpIk+jREnY6MRiNOnz5Nr9fDtm3a7TYPDzwe3DWBgu9aOmClXmJpaUlfy4RhSFEU2LbN0tISnU6HcrksdSJCCHGLkzAinpFq442iCIDJZILjOEwmE06dOkW326UoCprNJjuRw386NzvdeENzyJ31nNXVVeI4xrZtJpMJhmFgWRbLy8u0Wi1835fBZkIIISSMiGc2nU719cx0OsVxHMbjMRcuXGBzc5M0TanVamSWy8fOlEkLg5NewLevJBw5coQkSfA8j+FwCMzCzfLyMp1Oh0qlIpt4hRBCABJGxDOI45g0TQnDkDSdjXJXdSKnT58mz/PZCPeyz38459ONTRp2yndvBLRbTfI8x/M8xuMxRVGQ57m+mlFBRDbxCiGEAAkj4grSNCWKIuI41hNWTdOk1+tx6tQpXTdSq9X4g70Sj44cLKPgb64NWG5WMQwD27b171EUBZ1Oh5WVFWq1mhSsCiGEOETCiDhEFZpmWXZorshkMuH06dP0+32KosD3fU5NHH53Z1Z8+h0LA+5a8nEcR594hGEIQL1e58iRI1SrVRzHkYJVIYQQh0gYEZqaJ6ICSRiGmKZJlmV6wipAuVxmUrj8h4s1CgxeUQ9447GSDiGlUonJZKJniZw4cYJarYbjOFKwKoQQ4jISRoQWhuGhSat5nmMYBjs7O1y4cIE8zzFNE8st8SsXagSZybIb890ncizLwjAMyuUyo9GIJEkolUqcPHlSn4jIhFUhhBBXIu8MAnh6sFmSJIfaefv9PmfPntVXLpVKhU9ueVwIHTwz53tPhFQ8l6IoKJfLTKdTkiTBcRxuu+02Go0GpVKJUqkkE1aFEEJckYQRodt35+eKFEVBkiR6sFmapvi+z+cPTB7qVwD47o2AjbZPHMeUy2W9NM8wDI4fP0673cbzPBzHkYJVIYQQz0jCyC1O1YfAbK5IGIYkSYJhGJw9e5Z+v0+apniex/bU4Lf3mgB8y0LAa9Y8oijC932yLGM6nVIUBevr66ysrOh9M1KwKoQQ4quRMHILUwWrRVHoNtwoirBtm52dHfb29giCANu2yUyH/7jdJilMbqumfOcJiyAIcF0X0zT1PJKlpSWOHTum60PK5fK1/jaFEEJc5ySM3MJUkapq41X7Z3Z2dtjc3NRtvI7j8pu7DQ4Sm7qd8b23p+RZim3bejx8GIYsLCxw++23U6lUdBCRzhkhhBBfi4SRW5QqWIWnr2fSNKXX63Hx4kUODg4oigLP8/iTYYVHxmVMCv7uHRllI9UTVieTCUmS0Gg0OHnypA4i6opGCCGE+FpkHvctKEkSXbCqQkkYhkynU86fP0+32wVm80Iuxh4P7NcB+M6jKatuSFEYlEolPaG1XC5z55130mq19GmJ4zjX8lsUQghxA5GTkVtMlmW6YFWFkvF4TBzHesJqlmWUSiUSy+M/bjXJMbinmXBvK9Sj3g3D0NNZT548SbvdxnEcKVgVQghx1SSM3EJUwSo83c6r6j1UEAmCgFKphGHZfGKzzjizWCpl/HfHE6AAwLZtxuMxACdPntSdM2romRBCCHE1JIzcQlTnTJZlegFeEAScP3+eg4MDJpMJvu9j2w6f3qtxdlrCNQv+7h0pRTI7FXFdl8lkQpZlHD16lJWVFSqV2dwRKVgVQgjxXEgYuUWo5XeqeyaOY/r9PhcuXGBnZ4cgCCiXy9i2zX/t+fxRd3bV8o7jMZVspOtA4jgmTVNWVlbY2NigXp/Vk0jBqhBCiOdKClhvAUmSkCQJeZ6T5zlxHDMcDtnc3GRra4soinS9xxfGFX53d3bS8fajKXeWJzhOCUC3ALfbbU6ePEmj0cA0TVzXlYJVIYQQz5mcjNzkVMFqURQ6iEwmE86cOcP29jZxHGNZFrZt82To8xubsyDy7asp39ye6tMO0zSZTCbU63VOnjypO2csy6JUKl3Lb1EIIcQNTsLITezSglU1+v2JJ57Q01UNw8AwDLbTMp+4WCfH4JtaCW9ZS/TWXsdxmE6nlMtlNjY26HQ6lEolmbAqhBDieSFh5CamClaTJME0TabTKU888QT7+/sMh0Ncd7Ztt5d7/PsLTZLC4M5awt84GlIUub66UWPfjxw5wtrami5UlYJVIYQQzwcJIzcpVbCapimWZTGdTjlz5gxbW1sMBgN8f7ZtN6DEv99sM81N1r2Ev310gm3OAobjOHoL78LCAkePHsXzPD1h1TTlHx8hhBBfP3k3uQmpyahZlmEYBkmScO7cOc6cOcNoNOuMCYKAqLD5ld1FhqnFgpvyPRsDys6sRqRSqehdNc1mkxMnTlAul3EcB9d1sW2pfRZCCPH8kDByk1EdL6pzBuD8+fOcPn2ayWTy9IbdwuTXuyvsRTY1O+MdK/u0/NkG3lqtRpIkFEVBtVrl+PHj1Go1SqUStm1LwaoQQojnlYSRm0ie57pOJE1nW3UvXrzIk08+yWg0wrKs2S6aLOd3hqucmzp4Zs5/t7jLerOMaZpUq1XCMCQMQ3zfZ3V1lXa7reeIyKh3IYQQzzc5a79JqM6ZoiiI4xjXddne3uaJJ55gMBjoXTJxnPDZ+CiPj0vYRsF/s7DD8ZarT0GCICCOY6rVKktLS2xsbFAqlbAsSwpWhRBCvCAkjNwkwjAkz3M9wKzb7fL444/T7/cxDIPpdEoQBDycHeHhoY9BwVsaW7xkoYTruvi+z3Q6JY5jKpUKy8vLHD9+XI+Al4JVIYQQLxQJIzcBNaJdDTCbTCY8+uijDIdDXcgaBAGP5sv84aABwLfXdnnFoolt23ieRxzHRFFErVZjZWWFEydOkOc5nufpWhEhhBDihSD/qXuDS9N0VgeSpsDshOSRRx5hOBwync4mqPZ6Pc7kHX6v3wHgdZV97l1M8X2fUqmkl+ZVq1VWVlY4duwYAL7v6+4ZIYQQ4oUi/7l7A1MTVfM817tnnnzySXq9HqPRCNu2uXDhAntmm08Nlykw+MZyn7/cnlCvtzEMg6IoCMOQRqPB+vo66+vresS7bdsyYVUIIcQLTk5GblCqYFUFkizLuHDhAru7u4zHY1zXZWdnh4FZ57eHR8gKg9tLI97Y2KfZbJBlmZ6u2mw2WVtbY3V1Fdd1KZVmdSQSRIQQQrwY5GTkBjQfRIIgIM9zDg4O2NzcJAgCiqJgb2+PbmzyW8Ex4sLkiBNwX/Ui7dYCaTq7oplMJrTbbVZWVlhbW8P3fd2+K1t4hRBCvFie08nIhz70IY4fP47nedx777089NBDz/jYj3zkI7zhDW+g1WrRarW47777vurjxdemTkKCICDLMvr9PhcuXCAIApIkYTKZsD+O+OT0dia5xYId8ebKWZY6LQDK5TJBENDpdFhaWmJ1dZVKpaKvZSSICCGEeDFddRj5+Mc/zv33388/+Sf/hM9//vPcc889vPnNb2Z3d/eKj//MZz7D3/7bf5vf+73f48EHH2RjY4O/+lf/KhcvXvy6X/ytaDqdkqYpk8mEJEkYj8dsbm4yGo1mk1XTlJ2DPr8b30k/c6hbCW/xT7HSmXXRlMtlwjCk0+mwsLDAkSNHqNfreo6IdM0IIYR4sRlFURRX84R7772X17zmNXzwgx8EZkWUGxsb/MAP/ADvec97vubzsyyj1WrxwQ9+kHe+853P6msOh0MajQaDwYB6vX41L/emEoYhSZIwGo10K++FCxfY399nPB4DcObcBX43uYMLSZWymfFdlSc40pjNESmXy0RRxOLiIu12m2PHjlGtVnUQkTkiQgghnk/P9v37qt594jjmc5/7HPfdd9/Tv4Fpct999/Hggw8+q99DXSW02+1nfEwURQyHw0Mft7ooikiShMFgQJqmZFnG3t4eg8GA8XiMYRhcuHCRzyYnuJBUcYyct5RPsewbeJ6nT0RUEFH7Zmzbxvd9CSJCCCGumat6B9rf3yfLMpaXlw99fnl5me3t7Wf1e/zwD/8wa2trhwLNpd7//vfTaDT0x8bGxtW8zJuOGkjW6/XIsowsy9jd3WUwGNDtdnEch62tbf4oWueJpIlJwZvKp1n3M6rVqg4iy8vLNJtNTp48SbVaxXEcfN+XEe9CCCGuqRf1P4f/xb/4F/zyL/8yv/qrv/pVF669973vZTAY6I/z58+/iK/y+pIkCUEQ0O129Sbe3d1der0eu7u7+L7P9vY2n5+2+WKyCMC3emc5Yo+p1+t6uuri4iLVapU77rhDDzuTpXdCCCGuB1dVrbiwsIBlWezs7Bz6/M7ODisrK1/1uf/n//l/8i/+xb/gd3/3d/nGb/zGr/rYUqkka+qZTVcdjUZ60Z0KIvv7+3S7XXzfZ2triy+MfB6K1gD4ZvcCd7gDGo0G5XKZNE1ZWFig0Whw1113US6XpXVXCCHEdeWqTkZc1+VVr3oVDzzwgP5cnuc88MADvO51r3vG5/3Lf/kv+Wf/7J/xyU9+kle/+tXP/dXeQlTL7ny9zPb2Njs7O+zv71OtVhmPx3ypZ/AH4ewa6+XONi9z96jVatTrdaIoot1u02w2ufvuu/F9X494F0IIIa4XV93Hef/99/Oud72LV7/61bz2ta/lAx/4AJPJhHe/+90AvPOd72R9fZ33v//9APwf/8f/wfve9z5+6Zf+/+3de2xU150H8O8d2zO2Z2yP7WE8Hr8NfoGxkxBwIEvTJoRH2AJK00CEWlDSpoqINmxbif5DTFStkjZVpU02oqkUoFG0aak2JFVJQYZiQgiPFMwzrO0h47fHYxt7PJ7XncfZP6hnMfjBEOzrGX8/kiV77rlH5/Dz0f1yfR//jcLCwvC1JTqdDjqd7j5OJXYEAgEMDAzA5XIhFApBpVKhvb0ddrsdTqcTOp0OLpcLV3vcqHPnQ0DCXFUfHo7vhFabAr1eD1mWYTAYkJ6eHg4ivGOGiIhmoojDyMaNG9Hb24tXX30VNpsNDzzwAA4dOhS+qLWtrW3UAW/37t2QZRnPPPPMqH5qa2uxa9eubzb6GBQIBNDf3x9+wqokSbBarejr6wsHEY/Hg0tdTvxlMAdBxMEsDWBp3NfQalOQnp4efo5IZmYmysrKoNVqkZSUxAtViYhoRor4OSNKmC3PGZFlGTdu3IDX68VIWSwWC/r7++FyuaDT6eBwOHDaLqF+2IQQVDBITjwZ/78w6G8GkVAohJSUFOTm5qK0tBQ6nQ6JiYkMIkRENO3u9vjNx23OECO37o4EkWAwiObmZgwODobPiNjsvajrS8UVXyYAoDBuAI9IzchI1UGr1SIUCiE5ORm5ubkoKSlBamoqLwQmIqIZj2FkBvB4PHA4HPB6vQBuPmnVYrFgaGgITqcTWq0W1u4+HBwwoSugBQAs1tiwQOpAcrIWWq0WCQkJSExMRHFxMYqKiqDX66FWq5WcFhER0V1hGFGQEAJutxvDw8PhMyJDQ0NoaWnB4OAgZFlGfHw8Lrb245CrAMNCjQQE8Z2kNuSpBpCRYUAwGIRGo0FCQgLmzZuHwsJCZGZm8h0zREQUNXjEUkgwGITH44Hb7YbP54MkSbDb7eju7sbAwAA8Hg+EEPhHL3DcMxdBqJAqefHt+Ebka9XIzMyFLMvhO2QqKiqQn5/PIEJERFGHRy0F+P1+eL1eeDwe+Hw+BINBdHZ2ore3Fzdu3IAsy/B4fai/kYJL/psPk8tROfCteAvMc9KhVqsRCASgVqshhAgHkZGH0hEREUUThpFp5vP5IMty+A28siyjra0NN27cQH9//80/2/iC+JvDhM6QHgBQGdeFh+I7kZmRDiEEJElCXFwc1Go1SkpKUFhYCIPBwGeIEBFRVGIYmSZCCHi9XgQCAfh8PgQCAbhcLlitVgwPD6Ovrw/Dw8PodgkcdhfCiSTEIYgayYIFyR7odOkIBALQarVQqVQwGo0oLi6G0WhEZmYmb90lIqKoxTAyDUKhUPghZrIsIxgMor+/H21tbfB4POjo6IDf78fVwTic8BcjIMUjGV4sxzUU6xPg99+81VervXnnTElJCXJzc5GWlgadTscgQkREUY1hZIoFAoHwnTKBQACBQABdXV1ob2+Hx+NBZ2cnfD4ZXzhScRkFgAQY4cBj8c3I0KrhdruRkZEBtVqN7OxslJaWwmg0IjU1le+YISKimMAwMoVkWYbP5wNw8880siyjtbX1nwHEh5aWFjg9Mo6589AVZwQAzA11YFliNyCC8PkE0tPTkZqainnz5qG4uBh6vR7Jyck8G0JERDGDYWSKjFyoCgCSJMHtdsNiscBms2F4eBgdHR3ocQVxLFAKZ1wKVAjhgUAjFiQ54HF7kJKSAoPBAJPJhIqKCpjNZuh0Ot62S0REMYdHting8XgQCAQAAPHx8RgYGEBjYyPsdjsGBgbQ3d2NFm8ivhDz4VepkQgZS+SLyIpzwen0w2g0wmAwoLS0FGVlZUhLS+OL7oiIKGYxjNxHQgh4PB4Eg0EAgFqths1mQ1NTE3p7e9HX14eeHjsu+TJxOa4EUElIDznwSOAy4mQn/PHxyM/PR3Z2NiorK5Gfnw+tVsuzIUREFNN4lLtPbr1jRpIkqFQqdHZ2oqmpCTabDT09PegdcOC0vxAdCTkAgLxAJ6rkr+D3uqFLS0NRURHmzZuHyspKZGRkQKPR8GwIERHFPIaR+2Dk0e4jDyQLBAJoaWlBZ2cnOjo60N3djV5XACfEAgwl6CEJgXLvVeS4LfAEgygsLERRURGqqqpQVFSE5ORkng0hIqJZg0e8b+jWW3cBwOFwwGq1ore3Fx0dHbDZbOjwaXBK9SDkuESohYwFjjPQuTqhSkpC+fz5KC8vR3V1NebMmcOzIURENOswjHwDI++YGfm+p6cHra2t6O7uRnt7O2RZxhVPGi4lVEBIKqQEh1DeexyJQRf0ej2qq6tRWVmJ8vJy3ilDRESzFo9+9+jWW3eHh4fDzw/5+uuv0dfXB7+kxj8CeWhT5wIA5njaUWj/HJo4ICc3F4sXL0ZVVRXy8vJ4NoSIiGY1hpF7MPKSu0AgALvdjra2NlitVly/fh0ur4yv4/PRFF+EQHwCIATyBi/ANHAJ2uRklJeXY9myZSgrK+NbdomIiMAwEpFbX3bn9XrR0tICi8WCa9euocduh02di0ZNGbxxyQAAnX8QOfbT0Mt2pOn1qKmpwcMPP4zi4mK+U4aIiOifGEbu0sgzRPx+P3p7e2GxWNDU1HQziIhUNCU/imF1BgBAHXSjYPAijJ4WQAjkFBRg+fLlmD9/PgoKCnhtCBER0S14VLwLI88QkWUZTU1NaG5uxuXLl3G9dxgWbTUGtXkAAFVIRlb/ReR5LIhHCPr0dOTn5+ORRx5BRUUFjEajwjMhIiKaeRhGJjHyDJGhoSFcuHABFy9exOXmFjSpS9A3518ASQWIENL7r8A8cBGpagkZGRmYP38+CgsLYTKZUF5eDq1Wq/RUiIiIZiSGkQkEAgG4XC50dHTg5MmTOHfhEhrcGeidsxoiTg0ASB6wwGw/i/R4GTl5OeH3yWRnZyM3NxdarZYXqRIREU2AYWQcsixjcHAQDQ0NqD/+GU60+9BteBTB1FQAQMJQJwydJ2COd6NobhEWLlyI0tJS5OXlYc6cOVCr1bxAlYiI6C4wjIzB5/Ohra0Nn376Kf7W0IJmbSWCOSYAgOQZQGb75yiM60fpghIsWrQIZWVlMJvN0Ov1UKlUCo+eiIgoujCM3MblcuHkyZPY8z+HcC6YD79pxc0Nshu6tpOYJzrxYHUlFi9eh8rKSphMJqjVamUHTUREFMUYRv5JCIGenh785+/34UCzH57s70BSqSCCASR1fon5og1Laxbi29/ejHnz5iEtLY1nQYiIiO4DhhHcDCIf/eUgXtv/BRymhyHlaCABiOu6hPn+Zjz1rcV48smtKCgoQEJCgtLDJSIiiimzPoz03RjA+u1voF1fBVXeo5AAoO9rlHuvYcva5Vi58kVkZGTwYlQiIqIpMmvDiBACz/zbLpyVcxBn/hZUAEJDPSgauoSdW/8Vy5e/CI1Go/QwiYiIYt6sDSMujw9fimLEpRsQ8jqRZT+H//r3TahZ/LzSQyMiIppVZm0Y0SUnokJuQov1Gn7/ytN4bNkmpYdEREQ0K83aMAIAh3//H0oPgYiIaNbjvalERESkKIYRIiIiUhTDCBERESmKYYSIiIgUxTBCREREimIYISIiIkUxjBAREZGiGEaIiIhIUQwjREREpCiGESIiIlIUwwgREREpimGEiIiIFMUwQkRERIqKirf2CiEAAENDQwqPhIiIiO7WyHF75Dg+nqgII06nEwCQl5en8EiIiIgoUk6nE2lpaeNul8RkcWUGCIVC6OrqQkpKCiRJum/9Dg0NIS8vD+3t7UhNTb1v/c5Us2m+nGvsmk3z5Vxj12yZrxACTqcTZrMZKtX4V4ZExZkRlUqF3NzcKes/NTU1pn8Zbjeb5su5xq7ZNF/ONXbNhvlOdEZkBC9gJSIiIkUxjBAREZGiZnUY0Wg0qK2thUajUXoo02I2zZdzjV2zab6ca+yabfOdTFRcwEpERESxa1afGSEiIiLlMYwQERGRohhGiIiISFEMI0RERKSomA8j77zzDgoLC5GYmIiamhqcPXt2wvZ//vOfUV5ejsTERCxcuBCffvrpNI30m3n99dexePFipKSkwGg0YsOGDWhsbJxwn3379kGSpFFfiYmJ0zTie7dr1647xl1eXj7hPtFa18LCwjvmKkkStm3bNmb7aKvpZ599hu9+97swm82QJAkff/zxqO1CCLz66qvIzs5GUlISVqxYgebm5kn7jXTdT4eJ5ur3+7Fjxw4sXLgQWq0WZrMZP/zhD9HV1TVhn/eyFqbDZHXdunXrHeNevXr1pP3OxLoCk893rDUsSRLefPPNcfucqbWdKjEdRv70pz/hpz/9KWpra3H+/HlUV1dj1apVsNvtY7b/4osv8Nxzz+GFF15AQ0MDNmzYgA0bNuDKlSvTPPLIHT9+HNu2bcPp06dRV1cHv9+PlStXwuVyTbhfamoquru7w1+tra3TNOJvZsGCBaPG/fnnn4/bNprr+uWXX46aZ11dHQDg+9///rj7RFNNXS4Xqqur8c4774y5/de//jXeeust/O53v8OZM2eg1WqxatUqeL3ecfuMdN1Pl4nm6na7cf78eezcuRPnz5/HRx99hMbGRqxbt27SfiNZC9NlsroCwOrVq0eN+8MPP5ywz5laV2Dy+d46z+7ubuzZsweSJOF73/vehP3OxNpOGRHDlixZIrZt2xb+ORgMCrPZLF5//fUx2z/77LNi7dq1oz6rqakRP/nJT6Z0nFPBbrcLAOL48ePjttm7d69IS0ubvkHdJ7W1taK6uvqu28dSXV955RUxd+5cEQqFxtwerTUVQggA4sCBA+GfQ6GQMJlM4s033wx/Njg4KDQajfjwww/H7SfSda+E2+c6lrNnzwoAorW1ddw2ka4FJYw11y1btoj169dH1E801FWIu6vt+vXrxeOPPz5hm2io7f0Us2dGZFnGuXPnsGLFivBnKpUKK1aswKlTp8bc59SpU6PaA8CqVavGbT+TORwOAEBGRsaE7YaHh1FQUIC8vDysX78eV69enY7hfWPNzc0wm80oLi7G5s2b0dbWNm7bWKmrLMv44IMP8Pzzz0/4wshorentrFYrbDbbqNqlpaWhpqZm3Nrdy7qfqRwOByRJgl6vn7BdJGthJqmvr4fRaERZWRleeukl9Pf3j9s2lura09ODgwcP4oUXXpi0bbTW9l7EbBjp6+tDMBhEVlbWqM+zsrJgs9nG3Mdms0XUfqYKhULYvn07Hn30UVRWVo7brqysDHv27MEnn3yCDz74AKFQCMuWLUNHR8c0jjZyNTU12LdvHw4dOoTdu3fDarVi+fLlcDqdY7aPlbp+/PHHGBwcxNatW8dtE601HctIfSKp3b2s+5nI6/Vix44deO655yZ8iVqka2GmWL16Nd5//30cPXoUv/rVr3D8+HGsWbMGwWBwzPaxUlcA+MMf/oCUlBQ8/fTTE7aL1treq6h4ay9FZtu2bbhy5cqkf19cunQpli5dGv552bJlqKiowLvvvotf/vKXUz3Me7ZmzZrw91VVVaipqUFBQQH2799/V//biFbvvfce1qxZA7PZPG6baK0p/T+/349nn30WQgjs3r17wrbRuhY2bdoU/n7hwoWoqqrC3LlzUV9fjyeeeELBkU29PXv2YPPmzZNeWB6ttb1XMXtmxGAwIC4uDj09PaM+7+npgclkGnMfk8kUUfuZ6OWXX8Zf//pXHDt2DLm5uRHtm5CQgAcffBAWi2WKRjc19Ho9SktLxx13LNS1tbUVR44cwY9+9KOI9ovWmgII1yeS2t3Lup9JRoJIa2sr6urqIn61/GRrYaYqLi6GwWAYd9zRXtcRJ06cQGNjY8TrGIje2t6tmA0jarUaixYtwtGjR8OfhUIhHD16dNT/HG+1dOnSUe0BoK6ubtz2M4kQAi+//DIOHDiAv//97ygqKoq4j2AwiMuXLyM7O3sKRjh1hoeHcf369XHHHc11HbF3714YjUasXbs2ov2itaYAUFRUBJPJNKp2Q0NDOHPmzLi1u5d1P1OMBJHm5mYcOXIEmZmZEfcx2VqYqTo6OtDf3z/uuKO5rrd67733sGjRIlRXV0e8b7TW9q4pfQXtVPrjH/8oNBqN2Ldvn/jqq6/Eiy++KPR6vbDZbEIIIX7wgx+IX/ziF+H2J0+eFPHx8eI3v/mNuHbtmqitrRUJCQni8uXLSk3hrr300ksiLS1N1NfXi+7u7vCX2+0Ot7l9vq+99po4fPiwuH79ujh37pzYtGmTSExMFFevXlViCnftZz/7maivrxdWq1WcPHlSrFixQhgMBmG324UQsVVXIW7eNZCfny927Nhxx7Zor6nT6RQNDQ2ioaFBABC//e1vRUNDQ/gOkjfeeEPo9XrxySefiEuXLon169eLoqIi4fF4wn08/vjj4u233w7/PNm6V8pEc5VlWaxbt07k5uaKCxcujFrDPp8v3Mftc51sLShlork6nU7x85//XJw6dUpYrVZx5MgR8dBDD4mSkhLh9XrDfURLXYWY/PdYCCEcDodITk4Wu3fvHrOPaKntVInpMCKEEG+//bbIz88XarVaLFmyRJw+fTq87bHHHhNbtmwZ1X7//v2itLRUqNVqsWDBAnHw4MFpHvG9ATDm1969e8Ntbp/v9u3bw/82WVlZ4qmnnhLnz5+f/sFHaOPGjSI7O1uo1WqRk5MjNm7cKCwWS3h7LNVVCCEOHz4sAIjGxsY7tkV7TY8dOzbm7+3InEKhkNi5c6fIysoSGo1GPPHEE3f8OxQUFIja2tpRn0207pUy0VytVuu4a/jYsWPhPm6f62RrQSkTzdXtdouVK1eKOXPmiISEBFFQUCB+/OMf3xEqoqWuQkz+eyyEEO+++65ISkoSg4ODY/YRLbWdKpIQQkzpqRciIiKiCcTsNSNEREQUHRhGiIiISFEMI0RERKQohhEiIiJSFMMIERERKYphhIiIiBTFMEJERESKYhghIiIiRTGMEBERkaIYRoiIiEhRDCNERESkKIYRIiIiUtT/ASZj3/aoJCP3AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot individual runs in light grey with small alpha (more transparent) and the average in blue\n",
    "\n",
    "G = nx.barabasi_albert_graph(1000, 4)\n",
    "\n",
    "# YOUR SOLUTION HERE"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "71yaTNYFBltI"
   },
   "source": [
    "As you can see, it shows an S-shaped curve that is typical of the spread of infectious diseases. The number of infected nodes grows exponentially at first, then slows down to saturate as the number of susceptible nodes decreases. \n",
    "\n",
    "This type of average is commonly done, but a care must be taken when aggregating epidemic (forecast) curves. [A very nice paper by Jonas Juul et al.](https://www.nature.com/articles/s41567-020-01121-y) showed that this type of aggregation (fixed-time descriptive statistics) can underestimate extreme cases. \n",
    "\n",
    "<img src=\"https://yyiki.s3.us-east-2.amazonaws.com/public/imgs/juul_lehmann-pitfalls_in_using_fixed_time_desc_stat_in_epidemic_curve_ensembles.png\" width=600>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the median curve (the black solid line) heavily underestimate how high the peak is expected to be. Our simulation does not show this because everything is aligned to the same time step, but this was a very important insight in the early days of the COVID-19 pandemic, when it is difficult to forecast when the epidemic curve would take off. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "id": "K0X_c5sIBltI"
   },
   "source": [
    "## SIR model \n",
    "\n",
    "SI model is useful to understand the early stage behavior of the epidemic. But the most commonly used model is the SIR model. In this model, the infected nodes can recover and become immune. For the sake of simplicity, you can assume that there is a recovery probability parameter that is the same for all nodes. In each time step, each infected node can recover and becomes immune with this probability.\n",
    "\n",
    "**Q: implement the SIR model.**\n",
    "\n",
    "1. Implement `run_SIR()`. \n",
    "2. Make a function that runs the SIR model and then produces a plot of \"S\", \"I\", \"R\" states (number of nodes) over time from the run. Show them in the same plot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_SIR(graph: nx.Graph, tmax: int, beta: float, mu: float, initial_inf: float):\n",
    "    # YOUR SOLUTION HERE"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "collapsed_sections": [],
   "name": "m10-networkepidemiology.ipynb",
   "provenance": [
    {
     "file_id": "https://github.com/yy/netsci-course/blob/master/m10-epidemics/epidemics_assignment.ipynb",
     "timestamp": 1648429000871
    }
   ]
  },
  "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}