{
 "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": 1,
   "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:  768\n",
      "secondary:  [8, 71]\n",
      "degree of index:  4\n",
      "degree of secondary:  [56, 33]\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": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAG1CAYAAAAV2Js8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAATMxJREFUeJzt3Xl0VPX9//HnzCSTPSEhkIUEguxrwi6uLFFEwYJFKVqlWLH6i1bEomJbKdWvVhHESioVK1i7oSguKIqCCAKyCohhJ+yQhUD2dWZ+f0wSEpJAAknuJPN6nHPPzNx75877Ep28cu9nMTkcDgciIiIibshsdAEiIiIiRlEQEhEREbelICQiIiJuS0FIRERE3JaCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRteRhdgKuz2+2cPHmSgIAATCaT0eWIiIhILTgcDrKzs4mMjMRsrvm6j4LQJZw8eZLo6GijyxAREZHLcOzYMaKiomrcriB0CQEBAYDzHzIwMNDgakRERKQ2srKyiI6OLv89XhMFoUsoux0WGBioICQiItLEXKpZixpLi4iIiNtSEBIRERG3pVtjIiLSbNntdoqKiowuQxqAp6cnFovlio+jICQiIs1SUVERycnJ2O12o0uRBtKiRQvCw8OvaHgbBSEREWl2HA4Hp06dwmKxEB0dfdFxZKTpcTgc5OXlkZqaCkBERMRlH0tBSEREmp2SkhLy8vKIjIzE19fX6HKkAfj4+ACQmppK69atL/s2mSKyiIg0OzabDQCr1WpwJdKQykJucXHxZR9DQUhERJotTY3UvNXHz9ctgtCyZcvo0qULnTp14q233jK6HBEREXERzb6NUElJCVOnTuWbb74hKCiIfv36MXbsWFq2bGl0aSIiImKwZn9FaNOmTfTo0YM2bdrg7+/PyJEjWbFihdFliYiIVDFkyBCmTJly2e8/fPgwJpOJ7du311tNzZ3LB6E1a9YwevRoIiMjMZlMfPTRR1X2SUxMJCYmBm9vbwYNGsSmTZvKt508eZI2bdqUv27Tpg0nTpxojNJFRETq5MMPP+S5554zugy34vJBKDc3l9jYWBITE6vdvnjxYqZOncqMGTPYtm0bsbGxjBgxonxsgboqLCwkKyur0tIQts4ew86/DOfk4b0NcnwREWl6QkJCLjlbutQvlw9CI0eO5Pnnn2fs2LHVbp8zZw6TJ09m0qRJdO/enfnz5+Pr68vbb78NQGRkZKUrQCdOnCAyMrLGz3vxxRcJCgoqX6Kjo+v3hErFZG+jd8EWCvMaJmiJiMh5DoeDvKISQxaHw1HrOiveGouJieGFF17g/vvvJyAggLZt2/Lmm29W2n/Tpk306dMHb29v+vfvzw8//FDlmLt27WLkyJH4+/sTFhbGvffeS3p6OgCrV6/GarWydu3a8v1ffvllWrduTUpKymX8Szc9TbqxdFFREVu3bmX69Onl68xmM/Hx8WzYsAGAgQMHsmvXLk6cOEFQUBDLly/nj3/8Y43HnD59OlOnTi1/nZWV1WBhSEREGkd+sY3uz35pyGcn/XkEvtbL+3U7e/ZsnnvuOZ555hmWLFnCww8/zI033kiXLl3Iyclh1KhR3HTTTfzrX/8iOTmZxx57rNL7z507x7Bhw3jggQd49dVXyc/P56mnnuKuu+5i1apV5cHr3nvvZceOHRw6dIg//vGPvP/++4SFhdXH6bu8Jh2E0tPTsdlsVX5YYWFh7NmzBwAPDw9mz57N0KFDsdvtPPnkkxftMebl5YWXl1eD1i0iIlIbt956K//v//0/AJ566ileffVVvvnmG7p06cJ//vMf7HY7//jHP/D29qZHjx4cP36chx9+uPz98+bNo0+fPrzwwgvl695++22io6PZt28fnTt35vnnn+err77iwQcfZNeuXUycOJHbb7+90c/VKE06CNXW7bff7lY/VBERqczH00LSn0cY9tmXq3fv3uXPTSYT4eHh5W1gd+/eTe/evfH29i7fZ/DgwZXev2PHDr755hv8/f2rHPvgwYN07twZq9XKv//9b3r37k27du149dVXL7vepqhJB6HQ0FAsFkuV+5gpKSmEh4cbVJWIiLgak8l02benjOTp6Vnptclkwm631/r9OTk5jB49mpdeeqnKtooTla5fvx6AjIwMMjIy8PPzu8yKmx6Xbyx9MVarlX79+rFy5crydXa7nZUrV1ZJxXWVmJhI9+7dGTBgwJWWKSIiUu+6devGzp07KSgoKF/3/fffV9qnb9++/PTTT8TExNCxY8dKS1nYOXjwII8//jgLFixg0KBBTJw4sU5hq6lz+SCUk5PD9u3byweHSk5OZvv27Rw9ehSAqVOnsmDBAt555x12797Nww8/TG5uLpMmTbqiz01ISCApKYnNmzdf6SmIiIjUu7vvvhuTycTkyZNJSkri888/55VXXqm0T0JCAhkZGUyYMIHNmzdz8OBBvvzySyZNmoTNZsNms/HLX/6SESNGMGnSJBYuXMjOnTuZPXu2QWfV+Fz+OuGWLVsYOnRo+euyHl0TJ05k0aJFjB8/nrS0NJ599llOnz5NXFwcX3zxRZNp7R7z3k1gMgGmy3jkMt9XuvQeD8P+YMBZi4jIlfL39+fTTz/loYceok+fPnTv3p2XXnqJn//85+X7REZGsm7dOp566iluvvlmCgsLadeuHbfccgtms5nnnnuOI0eOsGzZMsB5u+zNN99kwoQJ3HzzzcTGxhp1eo3G5KjLAAduKCsri6CgIDIzMwkMDKy3466aeRPDHJsuvWND8g6Cp48aW4OISAMoKCggOTmZ9u3bV2pMLM3LxX7Otf397fJXhIySmJhIYmIiNputQY7/pHka5KbzvwcH0bGVHzgcgKOWj9Rx/wseM4/De/eCIrCIiLg5BaEaJCQkkJCQUJ4o653JRDpB2HxbQ2MPp+7dAOcjIiLSBLl8Y2kRERGRhqIgJCIiIm5Lt8bcWVEOLBrlfF5tzzJz1XXV9T6rcZvZufT8OXQ2ZkRXERGRi1EQckc+wWD2BHsxHF576f2v1LFNCkIiIuKSFIRq0NC9xgzlGwIPrYXU3Vy8J5r94r3PHPaLvz/zOGyYByWFjX6KIiIitaEgVIMG7zVWatwb67FYTJhwziHjfAQwOe8ulb42VXptAs7fmSrbZi59PxX2u/D9lK43m+BncbE8eEOHBjs3Tu1wBiEREREXpSBkkG4Rgazdn052YYlhNRw/m9+wQUhERJq91atXM3ToUM6ePUuLFi2MLqfOFIQMsmjSQA6fyS29q+Rw3mmi9I5T2esKz7lgm93hKN+/yvsrbHPg3FDx9anMAp5cshO7BhUXERE3pyBkEIvZRIdW/oZ8dnJ6buN+YF46/HNM1Z5oVZ5Ti31Kn5stEDsBYq5r3HMREZF6V1RUhNVqNeSzNY6QNBy/VoAJbEVw6Bs4uAoOfA0HvoL9K2DfF7BvOez9HPYscy67P4WkjyHpI/hpKez6AHYtgR/fh52LYef/YMd/4Yd/wZfPGH2GIiL1bsmSJfTq1QsfHx9atmxJfHw8ubnOP2DfeustunXrhre3N127duVvf/tbpfceP36cCRMmEBISgp+fH/3792fjxo3l29944w06dOiA1WqlS5cuvPvuu5XebzKZeOuttxg7diy+vr506tSJTz75pNI+n3/+OZ07d8bHx4ehQ4dy+PDhStvPnDnDhAkTaNOmDb6+vvTq1Yv//ve/lfYZMmQIjzzyCFOmTCE0NJQRI0Zw//33M2rUqEr7FRcX07p1a/7xj39c1r9lbeiKUA2ada+xxhIYCZNXwZkDNfREq+l5TT3XSp9nHIKtC6G4wMizE5GmxOGA4jxjPtvTl/LeKpdw6tQpJkyYwMsvv8zYsWPJzs5m7dq1OBwO/v3vf/Pss88yb948+vTpww8//MDkyZPx8/Nj4sSJ5OTkcOONN9KmTRs++eQTwsPD2bZtG3a7HYClS5fy2GOPMXfuXOLj41m2bBmTJk0iKiqKoUOHltcwc+ZMXn75ZWbNmsXrr7/OPffcw5EjRwgJCeHYsWPccccdJCQk8OCDD7JlyxaeeOKJSudQUFBAv379eOqppwgMDOSzzz7j3nvvpUOHDgwcOLB8v3feeYeHH36YdevWAc4AdcMNN3Dq1CkiIiIAWLZsGXl5eYwfP/6KfgQXo9nnL6GhZp83UnJ6LkNfWU2Atwc//qkJju+TvBbeGQWhXeCRTUZXIyIuqMqs5EW58EKkMcU8cxKsfrXaddu2bfTr14/Dhw/Trl27Sts6duzIc889x4QJE8rXPf/883z++eesX7+eN998k9/97nccPnyYkJCQKse+9tpr6dGjB2+++Wb5urvuuovc3Fw+++wzwHlF6A9/+APPPfccALm5ufj7+7N8+XJuueUWnnnmGT7++GN++umn8mM8/fTTvPTSSxdtLD1q1Ci6du3KK6+8AjivCGVlZbFt27ZK+/Xo0YOJEyfy5JNPAnD77bfTsmVLFi5cWO1xNfu8uLeSAkj5CWe7odJRrMtHtK7YpujCbeYL2ht5gHfzCLki0rTFxsYyfPhwevXqxYgRI7j55psZN24cVquVgwcP8utf/5rJkyeX719SUlI+xMv27dvp06dPtSEIYPfu3Tz44IOV1l177bW89tprldb17t27/Lmfnx+BgYGkpqaWH2PQoEGV9h88eHCl1zabjRdeeIH33nuPEydOUFRURGFhIb6+vpX269evX5UaH3jgAd58802efPJJUlJSWL58OatWrar2fOqLgpA0PWWXmM8dgTeuqZ9j9v81jJpTP8cSEdfj6eu8MmPUZ9eSxWLhq6++Yv369axYsYLXX3+d3//+93z66acALFiwoEoQsVgsAPj4+NRPuZ6elV6bTKby22u1MWvWLF577TXmzp1Lr1698PPzY8qUKRQVFVXaz8+v6lWy++67j6effpoNGzawfv162rdvz/XXX395J1JLCkLS9ETEQfTVcDb5fLshh73ySNfVrrdXbnNUUfKaxj8PEWk8JlOtb08ZzWQyce2113Lttdfy7LPP0q5dO9atW0dkZCSHDh3innvuqfZ9vXv35q233iIjI6Paq0LdunVj3bp1TJw4sXzdunXr6N69e61r69atW5XG099//32l1+vWreNnP/sZv/zlLwGw2+3s27evVp/TsmVLxowZw8KFC9mwYQOTJk2qdW2XS0FImh4vf/j1l1d2jLIG2UfWOdsbiYi4gI0bN7Jy5UpuvvlmWrduzcaNG0lLS6Nbt27MnDmT3/72twQFBXHLLbdQWFjIli1bOHv2LFOnTmXChAm88MILjBkzhhdffJGIiAh++OEHIiMjGTx4MNOmTeOuu+6iT58+xMfH8+mnn/Lhhx/y9ddf17q+hx56iNmzZzNt2jQeeOABtm7dyqJFiyrt06lTJ5YsWcL69esJDg5mzpw5pKSk1DpwPfDAA4waNQqbzVYptDUUBSFxT6YK4xGBcz60Uzuo3KbIfEH7ohq2YQJPb/BuuKlYRMQ9BAYGsmbNGubOnUtWVhbt2rVj9uzZjBw5EgBfX19mzZrFtGnT8PPzo1evXkyZMgUAq9XKihUreOKJJ7j11lspKSmhe/fuJCYmAjBmzBhee+01XnnlFR577DHat2/PwoULGTJkSK3ra9u2LR988AGPP/44r7/+OgMHDuSFF17g/vvvL9/nD3/4A4cOHWLEiBH4+vry4IMPMmbMGDIzM2v1GfHx8URERNCjRw8iIxu+gbt6jdWgYvf5ffv2NcteY14eZmbe3sM5R5np/Bxk5fOWlc5VZq6w3rlP6XozF+xjOp8XcL7HYjbRs00Q3p4Wg8+6Bkc2wMJb6uFAJrj9deh7bz0cS0Su1MV6E4lry8nJoU2bNixcuJA77rjjovuq11gDaqxJV43gaXE2Ni4ssfP0hz82+OcNah/C4t8MvvSORoiIhZjrISO5hvZEFdoVOapbbweHzfmeYxsVhERELpPdbic9PZ3Zs2fTokULbr/99kb5XAUhNxQV7Mtjwzvx08kswIG9bO6yCo8OHNjtpY9l85eVbreXz13mqPC+C/dxUFBs58S5fI6cMWgQs9qw+sKvll3ZMdbOhpV/rp96RETc1NGjR2nfvj1RUVEsWrQID4/GiSgKQm7q8Zs6N/hn7DqRyajXv2vwz3EZBZmQfqD0nmHZYqncnshsqdDWqHSbxdO5iIi4sZiYGIxoraMgJFJfdn/iXOrKYoUxb0CvcfVfk4iIXJQmXRW5Uh3jIbg9eLcAryCwBjgHUPPwdoYckwW4yDxDtiI4tLqRihURkYp0RUjkSkXEwmPbL71feYNr2/mG1t/NhW//Ajmpzh5sZkuF7vqWC26pXdBt3z/M2cZJRGqkjtHNW338fBWERBpL2dgCFS/Eeng5H/d/6VzqwicEHtuhedJEqlE27URRUVG9TT0hricvz9kZ58JpQepCQUjESF1Hwd7lkH+2cnf8sq76dtsF60u78hdmQX4GZB4D7x5Gn4WIy/Hw8MDX15e0tDQ8PT0xm9USpDlxOBzk5eWRmppKixYtyoPv5VAQqkHFARVFGkyrzvDAV3V/36xOkJsKm98C//DS22cVeqKZLRWel95K8wqELrfqdpq4BZPJREREBMnJyRw5csTocqSBtGjRgvDw8Cs6hoJQDZrzgIrSDFh9IRfY8nbd3jfsj3DD7xqkJBFXY7Va6dSpU5VZz6V58PT0vKIrQWUUhKTBFZTYWPHTacwmExazcxoOi9lUPi2H8zmYS9dZTJX3sZRO5WEp3d9sBg+zmbBAL0ymi/TGas5unwe7P3XeLrPbzt82s1e4hVZxfUoSZByE3HSjKxdpVGazWVNsyEUpCEmD8SidyuNcXjEPvru13o9/Z78oZt0ZW+/HbRLaX+9camvln50jYB9cBR8lnL+NVnYLraxXmtnj/LqO8dB2UMOdg4iIC1AQkgbTuXUAEwa25UBqNja7A1vpNBw2e+m0HnYHttLpOOyl62x2h3MfhwObnfLnFbeX2O0U2xz8cOyc0afYdPgEOx/T9zqX2tj2T/hdLfcVEWmiFISkwZjNJl68o1e9H3fDwTNMWPB9vR+3Wes3CbwCnNOAlN0yK7uNZreBveT8uoJzsP3fzsczB8+PaVTp6lHpVCFmi3PgyLJhAEREmhgFIRF34OUP/X5Vu33PHnYGoZICeL3vpff38IYJ/4UOw66kQhERQ2hgBRGpLCjaGWq8g5xd7j39nGHH7OlsR3ShkgI4qit0ItI06YqQiFRmtsC9S2veXnGqkOVP1r0Lv4iIC1EQEpG6qThViKl0DI9D3zrbGZk9StsRmSs8tzifB7ZxDuioEX5FxIUoCInI5fPydz4e+965XMrEZXXr9i8i0sAUhETk8g16GCxWKMyp0POsrBdahUEdD66CvDPORUTEhSgI1UBzjYnUQkAYDH3m0vu9PRKOrm/4ekRE6kg362uQkJBAUlISmzdvNroUkeYjNw0yj0P2aed0H/lnoTAbivPBVuxshC0i0oh0RUiarCNncrlpzreYS+cmK5uHzGwyYcI5P5m5bH11+5Ru9zCbuKt/NDf3uLIZjKUWPv+dc6mJ1R9+8R+46sbGq0lE3JqCkDQ5UcE+mE1QbHOwPzWnXo55/Gy+glBD6nkHpO2GkiJn+6Gy9kQXKsqB1X+BY5vO9zazeJ6fA83sCWE9oE0tBnoUEakFk8Oha9EXk5WVRVBQEJmZmQQGBhpdjpQ6fjaPE2fzsZfOX2Z34JyzzOHAUfrcUb6u5n0OpObwt9UHuaqVH6ueGGL0abkXh+N8w2p7Cax+ETbMu/T7zB7wu/3gG9LwNYpIk1Xb39+6IiRNUlSwL1HBvld8nE3JGfxt9cF6qEjqzGQCi4dzAbjm0dIeaNlgL3aGI1tpSLIXO0PT3uXO11sXQUBEhStFpVePzB7Qoh207mroqYlI06EgJCKuISAc4mdcfJ+XO0BeOqycefH9xr7pvIVWFo4qBiX/MA3qKCLlFIREpOkY8QIkfVR6lajieEWlV46Ol/byXPpgzceIuR5+taxRyhUR16cgJAKcySnipS/2YCnraWZ29jSzlD6aTWAxO3uaWSo+Nzu3+Vo9GN6tNb5W/S/VoGLHO5eabHkbvptb2hW/Qvsju825zlaoCWJFpBJ9a4tb8/NyzpWVmV/MG1fYVihhaAemjVDbFEP1v9+5VCfrJMzp1rj1iIjLUxASt9Y9IpDnx/TkyJlc7A6w2R04HA5sDgc2u7O3mc1+vseZ87mjwnM4lJbDwbRc0rOLjD4dqQ17CSybWqFrfmm3fIsnBMdA7/Glk8qKiDtQEBK3ZjKZ+OXV7a7oGInfHGDWl3vrqSJpMJ6+zvBjL4Et/6h5v/R90HZwhaDk6ezZ1qId+IU2Xr0i0igUhETEPfi0gLvfg5PbzrcZKuuWbyuGTX937rd2dvXvt3jBlB+d86uJSLOhICQi7qPjcOdSnZhrYePfoaQQbKUjYJeFpXPHnA2tv3neeWXIYnVeLbJ4OgNSx+EQGNm45yIi9UJBSEQEoPvPnEt13rgWUnbBtn9Wv73tNXD/8oarTUQajIKQSD3ZfTqL+d8exKO0y72HpfTRbMJsdj5ayhbT+ed+Xh70bxeMh0WD/LmsUa/Crg+cV4psRc4Rr21FkJMCR9ZBahKs/DOYSke5NlmcgzaaLOcHdKy4rmVHaH+90WclIigI1SgxMZHExERstmomhhSpwMfT2QV/5/FMdh7PvKxjPB7fmcfiO9VnWVKfogc6lwud2gF/vwEKztXctqgmj2yF0I71Up6IXD5NunoJmnRVLuVsbhHzvz3I2bwibHaw2e3YHKWPdkf5UmI/3+2+bDl5roDTWQWM7x/NS+N6G30qUlcOB2x+CzKSSwdwtJ1/rPi84uOBVVCcC51udrYrsniBh9XZ7ijmOugwzOizEmkWNOmqSCMJ9rMy/dbLG6hPXe+bOJMJBk6u23veindOBbJ/RdVta2fD/V+Cp0+F22rVTCxr9gCfYI13JFIPFIRERBrTHQtg35dQUuBsZ1RSCAWZsHmBc/vbI2p/rF/8Bzy8wMMbPHwgpD34hjRM3SLNlIKQiEhjCmkPVz9UzQYHHPga7BUmka00uWyFudPK/O/uqoeJ+yX4tzrfxT+sJ3QZ2WCnI9LUKQiJiLiC22rZ2NrhgG9egKMbnFeVSgqcV5XS9zm3b/9X1fc8ngRBbeqvVpFmREFIxAUkp+fywdbjeFhMeJjNpY8mPCxmPEu72XtYzKXrnPtEtvAmwNvT6NKlsZlMMOz3Vdef2Ap7l1fu3r/tn86BIAsyFYREaqAgJGIgT4uzseumwxlsOpxRp/cGeHuw7ulhBCoMCUCbfs6lop+WQl4hbP83BEWXjoRtdc6lpq77IoCCkIihbo9tw97TOZzLK6LY7sBmt1Nsc1Bic3a9L7Y5KLHbKbE7KLE5StfZScspJLughJPn8gkMVxCSGnj6Oh83zKu83q81TNvf+PWIuCAFIREDhQd5M/uu2Dq/r//zX5OeU9gAFUmzMnouJH18vndaYTYc+ApyU+Hrmc6rQx5W8AqEHmPBL9ToikUanYKQiEhzdeEks4XZ8FKMs+fZd3Mq75ua5JxKRMTNKAiJiLgLrwC4659wbCOUFDkbUp/eBcc3QeZx51JxfrRKAzp6OudKE2lmFIREmrAVP6Ww+1QWnhYznhYzVouzx1nF154ezl5mgT4etA7wNrpkMVrX25xLmU0LnEFo/wp4tUfN7/PwgTsXakwiaXYUhESaIGtpb7M5X+2r0/tm3xnLz/tFNURJ0lS1vwGCYyAn7fyAjY5qJpsuyYddH0LLTuDp7WyIrVGspRnQpKuXoElXxRV9/uMpPt5+gmKbsxeZc3E+Lypx9jIrttkpLrFTZHOQXVBMYYmdSdfGMGP0Rf7qFwHnoI2OCiNcr34R1r9edb9BD8HIlxq/PpFa0KSrIs3Yrb0iuLVXRK33n/XlHhK/OdiAFUmzYjKVthWyAF7Q4w7Y/zXknXGOZF2Y5dxv43zIP+uc68zTF7re6rzCJNKEKAiJuJElW47z7b40rBYzVg/z+ccLnnt5mLmpexjDuoYZXbK4gjZ9IeH7869P7YS/X+98vnPx+fVJH8ETexq1NJErpSAk4gY6tvYHILuwhOy0kkvs7fTlTyls++NNDVmWNFXhveCXH0LGISjOh6yTsPENyMuA9fPA08d5hSjmWmjR1uhqRS5KQUjEDYztE0Wf6GDO5hVRVGKnqLQtUdnzwtLnxTY7qdmFvLH6IPlF1TSYFQHnrbOOw4HSMYrKgpCtEFZUmAct5Cr47Q+GlChSWwpCIm4iJtSPGPwuud+xjDzeWK32RFIHgZEw+jU4sc15hSgvHQ6uguzTkHumdI4zz9KxiCzOICXiIhSERETkyvX7lXMByEiGv8ZBcR7MuuqCHU3OUGT1h1FznFN7iBhIQUhEqpVfbKPzH5bj5WHG29OCV2kj6rLnZY8hfl48cXNnIlv4GF2yuIqgaIga4LxCVGVMIodz7rP8DEj6BLr9TCNWi6EUhESkkrBAb2Ja+nL4TF55O6Lsgos3sL6qlR8JQzs2UoXi8iwe8MDXzud2O9iLneHHVuxctrwN3/4FfvoQflrqbFgdHAO/XOK8zSbSiNxiQMWxY8eyevVqhg8fzpIlS+r0Xg2oKO6oxGYnM985CGNBsY3CEnul52WPizcfZd2BM1x9VQjDurbGx9OCl6cFH08L3qWPPlYzXh4WIoK8aenvZfSpiSs4+QP8cwwUnKu8vk1/6HsfePlDu2shINyI6qSZqO3vb7cIQqtXryY7O5t33nlHQUikHv1l+R7mf1u7htWeFhPLHr2eLuEBDVyVNAl2u3PajqJcWDQK0vdW3t6yEzy6xZjapFnQyNIVDBkyhNWrVxtdhkizc/+1MZhNcC6/mIIiGwUlNvKLbOQX2ygodl45yi+2cepcAUU2OwfTchSExMlsBqufcxn3D9iyEAoyITcVktdA9imjKxQ3YXgQWrNmDbNmzWLr1q2cOnWKpUuXMmbMmEr7JCYmMmvWLE6fPk1sbCyvv/46AwcONKZgESnXOtCbJ2/pesn97pq/gU2HMxqhImmSwns5e5DB+R5nIo3E8CCUm5tLbGws999/P3fccUeV7YsXL2bq1KnMnz+fQYMGMXfuXEaMGMHevXtp3bo1AHFxcZSUVG3MuWLFCiIj69bwrrCwkMLCwvLXWVlZdTwjERG5YiUF8PEjpeMOWS54NFez3lz5tac3dLsd/FsbfSbi4gwPQiNHjmTkyJE1bp8zZw6TJ09m0qRJAMyfP5/PPvuMt99+m6effhqA7du311s9L774IjNnzqy344nIeesOpFNss+Nn9cDXy4Kv1QM/qwVfLw98PS34eXlg9VBXarfmFegMOvYS+OHdKzvW0Y3w8wX1U5c0W4YHoYspKipi69atTJ8+vXyd2WwmPj6eDRs2NMhnTp8+nalTp5a/zsrKIjo6ukE+S8RdeHk6w82/Nx7l3xuP1rift6eZN+/tzw2dWzVWaeJq/FrCPUvg9E6w28BhL320XfBY3Xq7M0BlHILjmyDlJ9jxPzB7lF5F8nBeMSp7NFkqryu7ouTpA6GdNQK2m3DpIJSeno7NZiMsrPIM2GFhYezZU/sZjuPj49mxYwe5ublERUXx/vvvM3jw4Gr39fLywstLXXxF6tPjN3Um1N+L7IIS8opKyC2ykVdYQl6Rrfx1UYmdgmI7m5IzFITcXcfhpXOZXaYflziDUOpPsPQ3l3eMwY/AiP+7/BqkyXDpIFRfvv76a6NLEHFrfdsG07dt8EX3mfHxLt7ZcKSRKpJmrdNNzvGIsk87rxbZSypcNSo5/7psm6PCPoVZzt5rabX/Y1uaNpcOQqGhoVgsFlJSUiqtT0lJITxcA22JNCcm3YaQ+uIdBLe/fnnv3f4f+OhhOLUTlvwaLFbnSNlmTwhuB1cnOF9Ls+HSP02r1Uq/fv1YuXJleZd6u93OypUreeSRRxr0sxMTE0lMTMRmu3CeHBERabZ8Q52Puamwq5oBeCP7QPsbGrcmaVCGB6GcnBwOHDhQ/jo5OZnt27cTEhJC27ZtmTp1KhMnTqR///4MHDiQuXPnkpubW96LrKEkJCSQkJBQPjKliDSOtfvTsDscBHh74u/tQaC3B/5eziXA25N2LX3x8zL8q0uaq47DYfy/nQM6ls2PZi+Gzf9wrivKM7pCqWeGf5ts2bKFoUOHlr8u67E1ceJEFi1axPjx40lLS+PZZ5/l9OnTxMXF8cUXX1RpQC0iTVugjycAO45nsuN4Zo37hfhZ+e6pofhaDf/6kubIbIFuo6qu37vcGYQOrwVbIUQNhMCIxq9P6p1bzDV2JTTXmEjjyMwv5qMfTpCeU0h2QQnZBSXkFBaTU1hCTunrQ+m5AHzzuyG0D/UzuGJxK2+PhKPrz78OiICpu9XF3oVprrErpDZCIo0ryMeTidfEXHSfXn/6kuyCqqPIizS4oc/A5gWQl+G8KpR9ChwOBaFmQFeELkFXhERcR1kQ+tU1MbQN8SXIx9O5+DofA709CfGzanRqaTh5GfBye+fzZ886p/YQl6QrQiLS7AR4eZBdUMKi9Ydr3CfEz8rnv72e8CDvxitM3FPqT87u9WaPyouHFXwuPm6WuA4FIRFpMubd05evklLIyi8ms3TJyi8mq6CEzPxizuYVkZFbxJ7TWQpC0vDmX1fztn6/gtGvNVopcvkUhESkybjUCNW3/XUtP53M4uS5AlKzCmjhq9tkUs98gqHXXXD4O2e3+rIRqW1lz4ud+x3+ztg6pdYUhESk2TCXNlx9ZumPPLPUuc7PaiHYz0qwr5UWvp6M6h3B+AFtDaxSmjST6eIz2h9ZDwtHNl49csX0p1INEhMT6d69OwMGDDC6FBGppYnXxNChlR8hflbMpZ15cotsHD+bz48nMlm7P52ZnyYZW6SIuBRdEaqBRpYWaXrG9YtiXL8oAOx2B1kFxZzNKyYjt4jk9Fx+9/4Oim12g6sUEVeiICQizZLZbKKFr5UWvlbah/rRpoWP0SWJiAvSrTERERFxWwpCIiIi9S0nFZY/Bd/Ph5JCo6uRi9CtMRFxK8U2B4Ne+JpWAV6E+nvRyt+LVgHOpVtEIIPah2DStAlyuXxCnI+FWbBxvvN5UFT1E7mKS1AQqoHmGhNpXkL9rXRs7c+B1BxSsgpJyar+r/TPfnsdPSLVQUIuU+uucM8SSE2CLQvhbDIUZhtdlVyE5hq7BM01JtJ82OwO0rILnUtOAWnZhaTnFJGWXcjSH06QmV/MO/cP5MbOrYwuVZqDd++AgythzHyIm2B0NW5Hc42JiFzAYjYRHuRdOv1G5as+mw9nkJlfbExhImIYBSERkQreXHOQ9QfSCQ/yJiLIm/AgH8IDvWkV4IXFrLZDchmOfAdWX7hqCHjrtqurURASEQHCAr356WQW6w6cYd2BM1W2RwZ5s3zKDQT5eBpQnTRJnqVjV/3wL+fS/kaY+ImxNUkVCkIiIsBrv4hj9d40TmcWcCqzgNNZ+c7H0tcnMws4kJpDv3Y1T/oqUsmQ6eAbAun74egGyD5ldEVSDQUhEREgwNuT0bGR1W67cdY3HDmTx+LNR9mfkk2bYB+ign2JbOGNl4elkSuVJiO8J9z+OhxeB4tuNboaqYGCUA3UfV5EyrTw8eQI8N6W47y35Xilba0DvHjoxg7cf117Y4oTkSui7vOXoO7zInLkTC6f/3iaE+fyOH42nxNn8zl+Np/8YucfSle18mPVE0OMLVJcV9kVodDO8Mhmo6txG+o+LyJST9q19OPhIR0qrXM4HKzem8akRZtBf06KNFkKQiIil8FkMuHn5fwKzSksYdWeFFr4Wgn2tRLs60mgtydmdbcXcXkKQiIil8nT4gw6qdmF3L9oS6VtZhME+XjSwtdKC19Pgis8dosI5Od922hOMxEXoCAkInKZeke1YPL17dlzOpuzeUWczS3mXF4RuUU27A44m1fM2bzqR6vecyqLUbGRdGztj7+XvopFjKLG0pegxtIiUldFJXbO5RdxLq+Yc3nFnM0r4lxeEWfzivnL8j1V9o8M8qZjWAC39Qpn/IC2BlQsDaqssbRvS7jmt2DxBLMHmC2ljxVel28rXSye0Ka/c2RqqRM1lhYRMYjVw0zrAG9aB3hX2XZdx1De33KMA2k57E/JITW7kJOlAzauO5DOz+La4O2psYmalbIRpvPOwNcz6v7+DsPg3qX1W5OUUxASEWlEPdsE0bPN+fmmMvOK+elkJne/tRGb3cGrX+0jIsibVgHO+c3KFj+rRW2KmqqIOIj/E5w5CHYb2EvAXlz6aANb2fMLloJMyDgEZw8bfALNm4JQDTSgoog0hiBfTwa0D8HXaiGvyMbf1xyqdj8fT8v5YOTvRbuWvsRGt6B3VBBtWvgoJLkysxmue7zu7zv6Pbw9ov7rkUrURugS1EZIRBrD1iNnWX8gnbScQtJzCknLPr/kFl38D7JQfyu9o1oQG9WC3tFBxEa1IMTP2kiVS4MpC0IhV8FvfzC6miZHbYRERJqQfu2Ca5zQNbewpHI4yilk7+lsdhw/x55T2aTnFLFqTyqr9qSWv6dtiC+9o4KI7xbGmD5tGus0RJocBSERERfn5+WBn5cH7Vr6VdlWUGwj6VQWO46dY+fxTHYcO8eh9FyOZuRxNCOPZTtP0balL33bVh+yRNydgpCISBPm7Wmhb9vgSkEnM7+YH49n8st/bATgjr+tp22IL/1jghkQE0J8tzBaBXgZVbKISzEbXYCIiNSvIB9PrusUyit3xtI9IhCTCY5m5PHhthNM//BHxiSuM7pEEZehK0IiIs3UuH5RjOsXRVZBMT8cPcd3+9NYsDaZE+fyjS5NxGUoCImINHOB3p7c2LkVPSIDWbA22ehypK7OHobZ3cDDChYv56OH9/nnFi/w8AKfYGc3/eB2RlfcpCgIiYiIuKKQq8DTD4pzIftk7d7jHQQ3zWzYupoZBSERERFX5N8antgNWSehpBBsRVBSACVFYCt0rispdD5P+gQOrnS+ljpREBIREXFV3kHO5VLOHXUGIakz9RoTERERt6UgVIPExES6d+/OgAEDjC5FRKTebT6cgc2uGZZEdGusBgkJCSQkJJTPVSIi0tT5eFrwMJsosTu4c/4Ggn09GdY1jPhurYlr24LwQG9N3ipuR0FIRMRN+Hl58M/7B/LelmN8szeNs3nFfLDtOB9sOw5AgJcHHcP86dTan06tA+gU5k+nsAAigxSQpPlSEBIRcSPXdAzlmo6hlNjsbDlylq+TUvh2XxqH0nPJLizhh6Pn+OHouUrv8bNa6BoRyJ9/1oMekbpCLs2LgpCIiBvysJi5+qqWXH1VS/4AFJXYOXwml30p2exPyWF/qvMxOT2X3CIbW4+c5cNtJxSEpNlREBIREaweZjqHBdA5LKDS+mKbneeWJfHPDUfUuLop+OlDOPkDmC1gMoPJBKay56VLpW3mGrabKqyr63Zz5dcVt1d8rzUAut4KVj9D/8kUhEREpEaeFjMB3s5fFUmnsli1J4WebYJoHeBtcGVSSWCk8zEnxbk0FUN/Dzc+aWgJCkIiInJRQT6eAGxKzmBTcgYArQO86NUmiB5tgujVJoiebQLV68xI/SZBaBcoyASHvXSxgcNx/rXdVmFbHbdfyXur2566G87sh9w0o//lFIREROTi7r06Bi8PC9uPnePHE5kcTMshNbuQlXtSWbkntXy/23pHkHh3XwMrdWNmC7S/3ugqam/V/8Gal42uAlAQEhGRS/CxWph4TQwTS1/nFpaw+1QWu05ksutkFtuPneNAag5r9hr/171IXdVpZOm3336bwkJN6CYi4s78vDzoHxPCr65tzyt3xvKPif0ByC4s4aF3t7L1SIbBFYrUXp2C0OTJk8nMzCx/HRkZyeHDh+u7JhERaUKign0Z2TMcgC9+Os3P39jA2L+tY/Hmoxw5k4vDod5m4rrqdGvswv+Ys7Ozsdvt9VqQiIg0LRaziTd+2Y99Kdn8Y20yS384UWlgxvBAbwa2D2HQVSEMah9Ch1b+alQtLkNthEREpF50DgvgpXG9+d2ILvx301HW7Etjx/FznM4q4JMdJ/lkx0kAWvpZGXRVCFNv6kzH1gGXOKpIw6pTEDKZTJVS/IWvRUREWgV48dvhnfjt8E7kF9n44dhZNh5ydr3fdvQsZ3KL+PzH0/hZPZh1Z6zR5Yqbq/Otsc6dO5eHn5ycHPr06YPZXLmpUUaGGsqJiIizx9k1HUK5pkMoAIUlNuas2Mff1xziXH4xDodDf1CLoeoUhBYuXNhQdbicxMREEhMTsdlsRpciItJseHlYaB/qnFLhq6QU7v3HJp65tRvdIwMNrkzclcmh5vwXlZWVRVBQEJmZmQQG6n9UEZErVVhi49Wv9vP2d8kU2eyYTHBnvyievKUrof5eRpcnjaFsQMWBD8KtsxrkI2r7+/uyGks7HA62bt3K4cOHMZlMtG/fnj59+ujypoiIXJKXh4WnR3blnkFteemLPSzbeYr3thzni12nmTaiC3cPaofFrN8n0jjqNI4QwDfffEOHDh0YNGgQd911F3feeScDBgygU6dOrFmzpiFqFBGRZig6xJd5d/flg4cH0yMykKyCEv748U+MSVzHrhOZlz6ASD2oUxA6cOAAo0aNIiYmhg8//JDdu3eTlJTE+++/T1RUFLfeeiuHDh1qqFpFRKQZ6tcuhE8euY4//6wHAd4e/Hgik8n/3GJ0WeIm6nRrbO7cuVx99dWsXLmy0vquXbsyduxY4uPjefXVV3n99dfrtUgREWneLGYT9w2OoX+7EG7961qy8ouNLkncRJ2C0OrVq3nxxRer3WYymZgyZQrTp0+vl8JERMT9+Hs5fy3lFtkYMusbYkL9aF9hiWnpR2QLH7UhknpTpyB09OhRevXqVeP2nj17cuTIkSsuSkRE3FN4kDedw/zZl5LD4TN5HD6Tx+oLZrW3ephpF+JLTKgfV4X6VQpLrQO81HFH6qROQSgnJwdfX98at/v6+pKXl3fFRYmIiHuyepj54rEbSMkuIDk9l8PpeSSn55CcnsfhM7kcOZNLUYmd/ak57E/NqfJ+X6uFmJbOUBQd4kt0iA9Rwb5EB/sQ2cIHb0+LAWclrqzO3eeTkpI4ffp0tdvS09OvuCAREXFvZrOJiCAfIoJ8uKZD5W02u4OT5/I5lJ7L4fRckkuXw2dyOZaRR16RjaRTWSSdyqr22GGBXuXBKCrYl6hgH6JDnI+RLXzwtNS5M7U0cXUOQsOHD68yCz042whpqHQREWlIFrOp9EqPLzd2blVpW1GJnWNn88oD0vGz+RzLyHM+nnWGpJSsQlKyCtl65GyVY5tNEB7o7QxIFa4kRQU7ryyFB3rjoaDU7NQpCCUnJzdUHSIiIlfE6mGmQyt/OrTyr7LN4XBwNq+4UjA6fjavUlgqLLFzMrOAk5kFbDpc9fih/l58+ui1RAT5NPzJSKOpUxBq165dQ9UhIiLSYEwmEyF+VkL8rMRGt6iy3eFwkJZTWCkYVQxKRzPySM8pJOlkloJQM1Ona3z79+9nwoQJZGVVvfeamZnJ3XffrQEVRUSkyTGZTLQO8KZv22B+FteGhKEdefGO3rz760GsnjaUXlEtjC5RGkidgtCsWbOIjo6udvKyoKAgoqOjmTWrYSZPExERMUpZ69eDaVV7qknTVqcg9O2333LnnXfWuP2uu+5i1apVV1yUiIiIK7mlZzgAL32xl5W7UwyuRupTnYLQ0aNHad26dY3bQ0NDOXbs2BUXJSIi4kp+c8NVjOsXhc3uIOE/29h6JMPokqSe1CkIBQUFcfDgwRq3HzhwoNrbZiIiIk2ZyWTixTt6MaxrawqK7fz2v9urHUpGmp46BaEbbrjhohOq/vWvf+X666+/4qJERERcjafFzF/ucE4zdeJcPspBzUOdgtD06dNZvnw548aNY9OmTWRmZpKZmcnGjRv5+c9/zpdffqlJV0VEpNnSyNPNT53GEerTpw9Llizh/vvvZ+nSpZW2tWzZkvfee4++ffvWa4EiIiIiDaXOU2yMGjWKI0eO8MUXX3DgwAEcDgedO3fm5ptvvuiErCIiIiKupk5BaNWqVTzyyCN8//33jB07ttK2zMxMevTowfz589VOSERERGpmKr3FuHsZdBgOXW4xrJQ63eycO3cukydPrnFAxd/85jfMmTOn3oqrD8eOHWPIkCF0796d3r178/777xtdkoiINAM2tZa+fN1vh8A2kH0S/jcBCo0bqLJOQWjHjh3cckvNqe3mm29m69atV1xUffLw8GDu3LkkJSWxYsUKpkyZQm5urtFliYhIE+RhMZU/v2XuGpb+cJwSm93AipqosB7wyGbnc4cdSgoMK6VOQSglJQVPT88at3t4eJCWlnbFRdWniIgI4uLiAAgPDyc0NJSMDA2EJSIidRfg7cmzo7oT5OPJwbRcHl+8g/g53/L+lmMUKxDVjdXP6AqAOgahNm3asGvXrhq379y5k4iIiDoVsGbNGkaPHk1kZCQmk4mPPvqoyj6JiYnExMTg7e3NoEGD2LRpU50+o8zWrVux2WxER0df1vtFRETuv6493z01lCdv6UKwryeHz+QxbclOhs1ezf82HaWoRIGoKalTELr11lv54x//SEFB1UtY+fn5zJgxg1GjRtWpgNzcXGJjY0lMTKx2++LFi5k6dSozZsxg27ZtxMbGMmLECFJTU8v3iYuLo2fPnlWWkydPlu+TkZHBfffdx5tvvnnRegoLC8nKyqq0iIiIVBTg7cn/G9KR754axjO3diXU38qxjHye/vBHhr6ymq+TNB9ZU2Fy1GGM8JSUFPr27YvFYuGRRx6hS5cuAOzZs4fExERsNhvbtm0jLCzs8ooxmVi6dCljxowpXzdo0CAGDBjAvHnzALDb7URHR/Poo4/y9NNP1+q4hYWF3HTTTUyePJl77733ovv+6U9/YubMmVXWZ2ZmavoQERGpVn6Rjf9sOsr8bw+Sll1Il7AAvnz8BqPLcn1/CnI+TjsIfqH1euisrCyCgoIu+fu7TleEwsLCWL9+PT179mT69OmMHTuWsWPH8swzz9CzZ0++++67yw5B1SkqKmLr1q3Ex8efL9hsJj4+ng0bNtTqGA6Hg1/96lcMGzbskiEInKNnl42YnZmZqUlkRUTkknysFn59XXvmTegDoPZCTUidB1Rs164dn3/+OWfPni0fULFTp04EBwfXe3Hp6enYbLYq4SosLIw9e/bU6hjr1q1j8eLF9O7du7z90bvvvkuvXr2q3d/LywsvL68rqltERNyTyWS69E7iUuochMoEBwczYMCA+qylQVx33XXY7UrmIiIiUpVLzx4XGhqKxWIhJaVyo7OUlBTCw8MNqkpERKR6ZXOynssvprDEZmwxUisuHYSsViv9+vVj5cqV5evsdjsrV65k8ODBDfrZiYmJdO/evUlc9RIREdfQIzKI1gFeZOQW8b9NamPaFBgehHJycti+fTvbt28HIDk5me3bt3P06FEApk6dyoIFC3jnnXfYvXs3Dz/8MLm5uUyaNKlB60pISCApKYnNmzc36OeIiEjz4e1p4bfDOwHw+qr95BaWGFyRXMpltxGqL1u2bGHo0KHlr6dOnQrAxIkTWbRoEePHjyctLY1nn32W06dPExcXxxdffFGvvdNERETqy/gB0SxYe4gjZ/JYuC6ZR4Z1MrokuYg6jSPkjmo7DoGIiEiZj7ef4LH/bSfAy4M1Tw4l2M9qdEmuqamNI+RO1EZIREQu1+jekXSLCCS7sISF6w8bXY5chIJQDdRGSERELpfZbOKRoR0BeHfDYfKK1FbIVSkIiYiINIBbeobTNsSXs3nFvL/luNHlSA0UhERERBqAxWxi8vXtAXjru0OUaNoNl6QgJCIi0kDG9YsmxM85M/3yXaeNLkeqoSAkIiLSQHysFu4b3A6Af31/xOBqpDoKQjVQrzEREakPV1/VEoD0nEKDK5HqKAjVQL3GRESkPmg+etemICQiIiJuS0FIRERE3JaCkIiIiLgtBSERERFxWwpCNVCvMRERkeZPQagG6jUmIiLS/CkIiYiIiNtSEBIRERG3pSAkIiIibktBSERERNyWgpCIiIi4LQWhGqj7vIiISPOnIFQDdZ8XERFp/hSERERExG0pCImIiIjbUhASERERt6UgJCIiIm5LQUhERETcloKQiIiIuC0FIRERkUbgMLoAqZaCUA00oKKIiNSHFr5WAA6l5TJnxV7sdkUiV6IgVAMNqCgiIvWhS3gAv7nxKgD+uuoA/+/f28grKjG4KimjICQiItLApo/sxqxxvfG0mPjip9OMe2MDJ8/lG12WoCAkIiLSKO7sH81/J19NSz8rSaeyuH3eOrYeOWt0WW5PQUhERKSR9I8J4eNHrqVreADpOYVMePN7Ptx23Oiy3JqCkIiISCOKCvblg4ev4ebuYRTZ7Ex9bwd/Wb5HjagNoiAkIiLSyPy8PJj/y34kDO0AwPxvD/Lgu1vIKVQj6samICQiImIAs9nEtBFdee0XcVg9zHy9O5Wf/209xzLyjC7NrSgIiYiIGOhncW1Y/ODVtArwYm9KNj9LXMem5Ayjy3IbCkIiIiIG69M2mE8euZaebQLJyC3inre+573Nx4wuyy0oCImIiLiAiCAf3v/NNdzWK4Jim4MnP9jJ88uSsKkRdYNSEBIREXERPlYL8+7uw5T4TgC89V0yUxZvx+FQGGooCkI10FxjIiJiBJPJxJT4ziTe3RdPi4lPd5xk2c5TRpfVbCkI1UBzjYmIiJFu6x3BI0OdV4ZmfPITZ3IKDa6oeVIQEhERcVEPD+lA1/AAMnKLmPlpktHlNEsKQiIiIi7K6mHm5XG9MZvgkx0n+SopxeiSmh0FIRERERfWO6oFD97gHIH690t/JDO/2OCKmhcFIRERERc3Jb4TV4X6kZpdyAuf7Ta6nGZFQUhERMTFeXtaeGlcb0wmWLzlGGv3pxldUrOhICQiItIEDIgJYeLgGACe/uBHcjVBa71QEBIREWkipo3oQlSwDyfO5TPry71Gl9MsKAiJiIg0EX5eHrx4Ry8AFq0/zObDmpz1SikIiYiINCHXd2rF+P7RADy1ZCcFxTaDK2raFIRERESamGdu60ZYoBeH0nN59et9RpfTpCkIiYiINDFBPp48P8Z5i2zBmkPsPH7O2IKaMAUhERGRJuim7mHcHhuJ3QFPLtlJUYnd6JKaJAUhERGRJmrG6O6E+FnZczqbv60+YHQ5TZKCkIiISBPV0t+LP93eA4DEbw6w93S2wRU1PQpCNUhMTKR79+4MGDDA6FJERERqNLp3BEO7tKLY5uDDH44bXU6ToyBUg4SEBJKSkti8ebPRpYiIiNTIZDLRPTIQQO2ELoOCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQiJiIiI21IQEhEREbelICQiIiJuS0FIRERE3JaCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQiJiIiI21IQEhEREbelICQiIiJuS0FIRERE3JaCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQiJiIiI22r2QejcuXP079+fuLg4evbsyYIFC4wuSURERFyEh9EFNLSAgADWrFmDr68vubm59OzZkzvuuIOWLVsaXZqIiIgYrNlfEbJYLPj6+gJQWFiIw+HA4XAYXJWIiIi4AsOD0Jo1axg9ejSRkZGYTCY++uijKvskJiYSExODt7c3gwYNYtOmTXX6jHPnzhEbG0tUVBTTpk0jNDS0nqoXERGRpszwIJSbm0tsbCyJiYnVbl+8eDFTp05lxowZbNu2jdjYWEaMGEFqamr5PmXtfy5cTp48CUCLFi3YsWMHycnJ/Oc//yElJaXGegoLC8nKyqq0iIiISPNkeBuhkSNHMnLkyBq3z5kzh8mTJzNp0iQA5s+fz2effcbbb7/N008/DcD27dtr9VlhYWHExsaydu1axo0bV+0+L774IjNnzqzbSYiIiEiTZPgVoYspKipi69atxMfHl68zm83Ex8ezYcOGWh0jJSWF7OxsADIzM1mzZg1dunSpcf/p06eTmZlZvhw7duzKTkJERKSBXdshlCnxnbixcyujS2lyDL8idDHp6enYbDbCwsIqrQ8LC2PPnj21OsaRI0d48MEHyxtJP/roo/Tq1avG/b28vPDy8rqiukVERBrTNR1Duaaj2r9eDpcOQvVh4MCBtb51JiIiIu7FpW+NhYaGYrFYqjRuTklJITw83KCqREREpLlw6SBktVrp168fK1euLF9nt9tZuXIlgwcPbtDPTkxMpHv37gwYMKBBP0dERESMY/itsZycHA4cOFD+Ojk5me3btxMSEkLbtm2ZOnUqEydOpH///gwcOJC5c+eSm5tb3ousoSQkJJCQkEBWVhZBQUEN+lkiIiJiDMOD0JYtWxg6dGj566lTpwIwceJEFi1axPjx40lLS+PZZ5/l9OnTxMXF8cUXX1RpQC0iIiJSVyaH5pu4qLIrQpmZmQQGBhpdjoiISPPxp9I7LtMOgl/99nqr7e9vl24jZCS1ERIREWn+FIRqkJCQQFJSEps3bza6FBEREWkgCkIiIiLithSERERExG0pCImIiIjbUhASERERt6UgVAP1GhMREWn+FIRqoF5jIiIizZ+CkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQjVQL3GREREmj8FoRqo15iIiEjzpyAkIiIibktBSERERNyWgpCIiIi4LQUhERERcVsKQiIiIuK2FIRqoO7zIiIizZ+CUA3UfV5ERKT5UxASERERt6UgJCIiIm5LQUhERETcloKQiIiIuC0FIREREXFbCkIiIiLithSERERExG0pCNVAAyqKiIg0fwpCNdCAiiIiIs2fgpCIiIi4LQUhERERcVsKQiIiIuK2FIRERETEbSkIiYiIiNtSEBIRERG3pSAkIiIibktBSERERNyWgpCIiIi4LQUhERERcVsKQjXQXGMiIiLNn4JQDTTXmIiISPOnICQiIiJuS0FIRERE3JaCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQiJiIiI2/IwugARERFxU2P/7ny0+htWgoKQiIiIGCP2F0ZXoFtjIiIi4r4UhERERMRtKQiJiIiI21IQEhEREbelICQiIiJuS0GoBomJiXTv3p0BAwYYXYqIiIg0EJPD4XAYXYQry8rKIigoiMzMTAIDA40uR0RERGqhtr+/dUVIRERE3JaCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQiJiIiI29Ls85dQNrpAVlaWwZWIiIhIbZX93r7UKEEKQpeQnZ0NQHR0tMGViIiISF1lZ2cTFBRU43YNqHgJdrudkydPEhAQgMlkqtdjDxgwgM2bN9frMV3x8xvqc+rzuPVxrMs9RlZWFtHR0Rw7dkyDdhrM6P8nG0tTOE+ja9T3Y/0dy6jvRofDQXZ2NpGRkZjNNbcE0hWhSzCbzURFRTXIsS0Wi6G/+Brr8xvqc+rzuPVxrCs9RmBgoIKQwYz+f7KxNIXzNLpGfT/W37GM/G682JWgMmosbaCEhAS3+PyG+pz6PG59HMvon6dcOXf5GTaF8zS6Rn0/1t+xjP5ZXopujYkYTPPZiYhU1VjfjboiJGIwLy8vZsyYgZeXl9GliIi4jMb6btQVIREREXFbuiIkIiIibktBSERERNyWgpCIiIi4LQUhERERcVsKQiIiIuK2FIREXNjYsWMJDg5m3LhxRpciIuIyjh07xpAhQ+jevTu9e/fm/fffv+xjqfu8iAtbvXo12dnZvPPOOyxZssTockREXMKpU6dISUkhLi6O06dP069fP/bt24efn1+dj6UrQiIubMiQIQQEBBhdhoiIS4mIiCAuLg6A8PBwQkNDycjIuKxjKQiJNJA1a9YwevRoIiMjMZlMfPTRR1X2SUxMJCYmBm9vbwYNGsSmTZsav1ARkUZWn9+PW7duxWazER0dfVm1KAiJNJDc3FxiY2NJTEysdvvixYuZOnUqM2bMYNu2bcTGxjJixAhSU1MbuVIRkcZVX9+PGRkZ3Hfffbz55puXXYvaCIk0ApPJxNKlSxkzZkz5ukGDBjFgwADmzZsHgN1uJzo6mkcffZSnn366fL/Vq1czb948tRESkWbpcr8fCwsLuemmm5g8eTL33nvvZX++rgiJGKCoqIitW7cSHx9fvs5sNhMfH8+GDRsMrExExFi1+X50OBz86le/YtiwYVcUgkBBSMQQ6enp2Gw2wsLCKq0PCwvj9OnT5a/j4+O58847+fzzz4mKilJIEpFmrzbfj+vWrWPx4sV89NFHxMXFERcXx48//nhZn+dxxRWLSIP5+uuvjS5BRMTlXHfdddjt9no5lq4IiRggNDQUi8VCSkpKpfUpKSmEh4cbVJWIiPEa+/tRQUjEAFarlX79+rFy5crydXa7nZUrVzJ48GADKxMRMVZjfz/q1phIA8nJyeHAgQPlr5OTk9m+fTshISG0bduWqVOnMnHiRPr378/AgQOZO3cuubm5TJo0ycCqRUQanit9P6r7vEgDWb16NUOHDq2yfuLEiSxatAiAefPmMWvWLE6fPk1cXBx//etfGTRoUCNXKiLSuFzp+1FBSERERNyW2giJiIiI21IQEhEREbelICQiIiJuS0FIRERE3JaCkIiIiLgtBSERERFxWwpCIiIi4rYUhERERMRtKQiJSJM0ZMgQpkyZYnQZItLEKQiJiIiI21IQEhGpRlFRkdEliEgjUBASEZeXm5vLfffdh7+/PxEREcyePbvS9sLCQn73u9/Rpk0b/Pz8GDRoEKtXr660z4IFC4iOjsbX15exY8cyZ84cWrRoUb79T3/6E3Fxcbz11lu0b98eb29vAM6dO8cDDzxAq1atCAwMZNiwYezYsaPSsT/++GP69u2Lt7c3V111FTNnzqSkpKRB/i1EpH4pCImIy5s2bRrffvstH3/8MStWrGD16tVs27atfPsjjzzChg0b+N///sfOnTu58847ueWWW9i/fz8A69at46GHHuKxxx5j+/bt3HTTTfzf//1flc85cOAAH3zwAR9++CHbt28H4M477yQ1NZXly5ezdetW+vbty/Dhw8nIyABg7dq13HfffTz22GMkJSXx97//nUWLFlV7fBFxQQ4REReWnZ3tsFqtjvfee6983ZkzZxw+Pj6Oxx57zHHkyBGHxWJxnDhxotL7hg8f7pg+fbrD4XA4xo8f77jtttsqbb/nnnscQUFB5a9nzJjh8PT0dKSmppavW7t2rSMwMNBRUFBQ6b0dOnRw/P3vfy//nBdeeKHS9nfffdcRERFx+SctIo3Gw+ggJiJyMQcPHqSoqIhBgwaVrwsJCaFLly4A/Pjjj9hsNjp37lzpfYWFhbRs2RKAvXv3Mnbs2ErbBw4cyLJlyyqta9euHa1atSp/vWPHDnJycsqPUyY/P5+DBw+W77Nu3bpKV4BsNhsFBQXk5eXh6+t7uacuIo1AQUhEmrScnBwsFgtbt27FYrFU2ubv71+nY/n5+VU5dkRERJX2RkB5+6KcnBxmzpzJHXfcUWWfsnZGIuK6FIRExKV16NABT09PNm7cSNu2bQE4e/Ys+/bt48Ybb6RPnz7YbDZSU1O5/vrrqz1Gly5d2Lx5c6V1F76uTt++fTl9+jQeHh7ExMTUuM/evXvp2LFj3U5MRFyCgpCIuDR/f39+/etfM23aNFq2bEnr1q35/e9/j9ns7OvRuXNn7rnnHu677z5mz55Nnz59SEtLY+XKlfTu3ZvbbruNRx99lBtuuIE5c+YwevRoVq1axfLlyzGZTBf97Pj4eAYPHsyYMWN4+eWX6dy5MydPnuSzzz5j7Nix9O/fn2effZZRo0bRtm1bxo0bh9lsZseOHezatYvnn3++Mf6JROQKqNeYiLi8WbNmcf311zN69Gji4+O57rrr6NevX/n2hQsXct999/HEE0/QpUsXxowZw+bNm8uvIF177bXMnz+fOXPmEBsbyxdffMHjjz9+yVtXJpOJzz//nBtuuIFJkybRuXNnfvGLX3DkyBHCwsIAGDFiBMuWLWPFihUMGDCAq6++mldffZV27do13D+IiNQbk8PhcBhdhIhIY5s8eTJ79uxh7dq1RpciIgbSrTERcQuvvPIKN910E35+fixfvpx33nmHv/3tb0aXJSIG0xUhEXELd911F6tXryY7O5urrrqKRx99lIceesjoskTEYApCIiIi4rbUWFpERETcloKQiIiIuC0FIREREXFbCkIiIiLithSERERExG0pCImIiIjbUhASERERt6UgJCIiIm5LQUhERETc1v8HWatRtGaq/wgAAAAASUVORK5CYII=",
      "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** question (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 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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": 3,
   "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": 10,
   "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": [
       "{112, 453, 489, 704, 924}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "init_infected_fraction = 0.005\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": 12,
   "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[0][\"state\"])\n",
    "print(G.nodes[112][\"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": 13,
   "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": 14,
   "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": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGxCAYAAACeKZf2AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAASN9JREFUeJzt3XlcFHXjB/DPLrDLvaDI6Qoo3gcgCh752EEe9ZhXimahZoc+mhrZo/0ejyzLyjQrTcv7qNTMo9I0JS0tFUXAGy8E5EaO5Vxgd35/kJvEIasLs7t83q/XvtJhZvhM48rH2fl+RyIIggAiIiIiMyEVOwARERGRIbHcEBERkVlhuSEiIiKzwnJDREREZoXlhoiIiMwKyw0RERGZFZYbIiIiMissN0RERGRWLMUO0Ni0Wi1SU1Ph4OAAiUQidhwiIiKqB0EQUFBQAE9PT0ildV+baXLlJjU1FUqlUuwYRERE9ACSk5PRsmXLOtcRtdz8/vvvWLJkCaKjo5GWlobdu3dj2LBhdW5z9OhRRERE4OLFi1AqlZg7dy4mTJhQ7+/p4OAAoPJ/jqOj40OkJyIiosaiUqmgVCp1P8frImq5KSoqgr+/P1588UWMGDHivusnJCTg6aefxuTJk/H1118jMjISL730Ejw8PDBw4MB6fc+7H0U5Ojqy3BAREZmY+txSImq5GTx4MAYPHlzv9VevXg1fX18sXboUANCxY0ccP34cn3zySb3LDREREZk3kxotdeLECYSGhlZZNnDgQJw4caLWbdRqNVQqVZUXERERmS+TKjfp6elwc3OrsszNzQ0qlQolJSU1brN48WIoFArdizcTExERmTeTKjcP4q233kJ+fr7ulZycLHYkIiIiakAmNRTc3d0dGRkZVZZlZGTA0dERNjY2NW4jl8shl8sbIx4REREZAZO6ctO7d29ERkZWWXbo0CH07t1bpERERERkbEQtN4WFhYiNjUVsbCyAyqHesbGxSEpKAlD5kVJ4eLhu/cmTJ+PmzZv473//iytXruCLL77Ajh078Prrr4sRn4iIiIyQqOXmzJkzCAwMRGBgIAAgIiICgYGBmD9/PgAgLS1NV3QAwNfXF/v27cOhQ4fg7++PpUuXYu3atRwGTkRERDoSQRAEsUM0JpVKBYVCgfz8fE7iR0REZCL0+fltUvfcEBEREd0Pyw0RERGZFZYbIiIiMissN0RERGQQWq2ADFUpEu8UiZrDpCbxIyIiInGpSsuRdKcYt3OLkZxTgqScYiTnFiMppxi3c0tQVqFFv7Yu2DIpRLSMLDdERESko67QICW3BMm5lcXl9j3lJTmnBPkl5XVubyGVoFyjbaS0NWO5ISIiakK0WgGZBerKwnLnnqsuOSVIzi1GuqoU95skprmdDMpmtpUvZxu00v3aFh5O1rCyEPeuF5YbIiIiM1Wh0eJGVhEupOTjQmo+LqaocDE1H0Vlmjq3s5VZQOn8V2FpZqP7datmtmjpbAM7uXHXB+NOR0RERPWirtDgWkahrshcSFHhcpoK6orqHxFZSCXwcrKpUlzuvQrTzE4GiUQiwlEYBssNERGRiSkp0+ByugoXUypLzIXUfFzNKEC5pvrnSfZyS3TydEQXTwW6eDmii5cCrV3sYCnyR0cNieWGiIjIiBWUluNSqgoXUv8qM6n5uJ5ZCG0N98U42Vqhi6cCnb3ulhkFvJvZQio13aswD4LlhoiIyEgUlJbj/O18nEvJx4WUfFxMVSEhu+Y5Y1zs5ej615WYzn9dlfFysjHpj5MMheWGiIhIBGUVWlxJVyEuOQ+xyfmIu52HG1mFNY5U8nKyQWfPyiLT5a+rMq6O1o0f2kSw3BARETUwQRBw607xX0UmD3G383AxVYWyGm72belsA/+WTroi09lTgWZ2MhFSmy6WGyIiIgPLLlQjLjmvsszczkdccl6Nk98pbKzgr3RCQEsF/JVO6NbSCS0c5CIkNi8sN0RERA+huKwCF1JUiE3ORVxyPmKT85CSV1JtPZmlFJ09HeHf0gkBysqXd3Nb3iPTAFhuiIiI6qGkTIPU/BKk5lU+luD87coiczWjoNrIJYkE8GthD3+l019XZpzQ3t0BMkvzHX5tTFhuiIioySur0CJDVYrUvBKk5ZfqSkxaXilS80uRll+CvOLan6nk7mgNf6VCV2S6tFTA0dqqEY+A7sVyQ0REZk2jFZBdqEZqXglS8yqLiu6/+aVIyytBVqH6vs9TAgA7mQU8nGzg+dfopbsfMbkrOHLJmLDcEBGRWYlNzsP208m4nlmA1LxSZKhKUVHTjHf/ILOUwkNhDQ+FNTydbOCpsIGHk7Xuvx4KGzhaW/IeGRPAckNERCavtFyDfefSsPnELcTdzq/2dQupBG4Ocng42cBDYQ2vv/7rcU+JaW7iz1Oiv7HcEBGRyUrNK8HXpxKxLSoZd4rKAFRegRnSzROPtm9ReQXGyRot7OVm/SwlqorlhoiITIogCDhx8w42/5mIXy6l60YqeSqsMa6XN8b0VKK5PeeKacpYboiIyCQUqSuwOyYFm0/cwtWMQt3y3q2bY3wfH4R2dOXVGQLAckNEREYuIbsIW04k4rvoZBSUVgAAbGUWGNHdC+G9fdDOzUHkhGRsWG6IiMjoaLUCfruahY1/3sJvV7N0y31d7PBCL2+MDGoJhQ3nkaGasdwQEZHRyC8ux3fRydhyMhGJd4oBVM72+1h7V4T39sa/2raAVMoRTVQ3lhsiIhLdlXQVNv2ZiD0xKSgp1wAAHK0tMbqHEi/09oZ3czuRE5IpYbkhIiJRlGu0OHQpA5v+vIVTCTm65R3cHTC+jw+GBnjCVsYfU6Q//qkhIqJGlZZfgu/O3MY3p5KQrioFUDnJ3qDO7gjv7Y1g32acTI8eCssNERE1uHKNFpGXM7H9dBJ+u5qlm5vGxV6GscGt8FxIK3gobMQNSWaD5YaIiBpMQnYRtp1OwvfRKcguVOuWh/g2w9jgVhjc1R1ySwsRE5I5YrkhIiKDKi3X4OcLadgWlVzlXhoXezmeDWqJ0T1aonULexETkrljuSEiIoO4kJKPHWeSsTsmRTfZnlQCPNreFWE9lXi8gyusOIMwNQKWGyIiemCq0nL8EJuKbaeTcCFFpVve0tkGYT2UeLZHS95LQ42O5YaIiPQiCALOJOZiW1Qy9p1PRWm5FgAgs5BiQGc3jOnZCn3aNOdkeyQalhsiIqqX7EI1dp29jW2nk3Ezq0i3vJ2bPcJ6tsLwQC80s5OJmJCoEssNERHVSqMVcOxaFrafTsahSxmo+GsMt63MAkO6eSIsWIlApRPnpSGjwnJDRETVJOcUY2f0bXx3Jhmp+aW65QFKJ4zpqcS//T1hL+ePEDJO/JNJREQAgNyiMvx0Pg17YlIQnZirW+5ka4XhgV4I66lEB3dHERMS1Q/LDRFRE1ZarsHhyxnYE5OCo/FZuo+dpBKgTxsXjO6pxIBObrC24kR7ZDpYboiImhiNVsDJm3ewOyYFBy6ko1BdoftaZ09HDA/0whB/T7g5WouYkujBsdwQETUBgiDgUpoKe2JS8ENcKjJUfz8KwcvJBsMCPTEswAtt3RxETElkGCw3RERm7HZuMfbGpmJPTAquZRbqlitsrPB0Nw8MD/RCUCtnzklDZoXlhojIzOQVl2H/+XTsiUlB1K2/n+0ks5QitKMrhgV4oX/7FnxgJZktlhsiIjNQWq7BkSuZ2B2TgiPxmSjXVN4YLJEAvXybY3igFwZ1dYejtZXISYkaHssNEZGJ0moFnEy4g70xqdh/IU33sEoA6ODugOGBXngmwJPPdqImh+WGiMgEHb+Wjdnfn0NKXolumYfCGkMDvDAs0JPz0VCTxnJDRGRitp9Owv92X0CFVoCDtSWe7uqBYYFeCPZpxhuDicByQ0RkMrRaAR//Eo8vjt4AAAwL8MQHI7txgj2if2C5ISIyAaXlGsz6Lg4/nUsDAEx/3A+vP9mOD6wkqgHLDRGRkcspKsPLm88gOjEXllIJFo/oilE9lGLHIjJaLDdEREbsZlYhJm48jcQ7xXCwtsSXzwehj5+L2LGIjBrLDRGRkYpKyMErW84gr7gcLZ1tsHFiT/i58vEIRPfDckNEZIT2xKTgvzvPoUyjhb/SCWvDe6CFg1zsWEQmgeWGiMiICIKAFb9ex9JDVwEAg7u4Y9noANjIOCKKqL5YboiIjERZhRb/t/s8dkbfBgC88q/WmDOoA+euIdITyw0RkRHILynHlK3R+PPGHUglwDtDu+D5Xt5ixyIySSw3REQiS84pxsSNp3E9sxB2MgusGNcdj7V3FTsWkcliuSEiElFsch5e2nQa2YVlcHe0xvoJPdHJk8+FInoYLDdERCI5cCEdM7fHoLRci04ejlg/oSfcFdZixyIyeSw3RESNTBAErD2WgPd/vgxBAB5r3wKfP9cd9nL+lUxkCHwnERE1ogqNFm//eBFbTyYBAMJ7e2P+vzvB0kIqcjIi88FyQ0TUSArVFZj2zVkcjc+CRAL876mOmPSILx9+SWRgLDdERI0gPb8UEzeexuU0FaytpFgeFohBXdzFjkVkllhuiIga2KVUFV7ceBrpqlK42MuwdnxPBCidxI5FZLZE/5B35cqV8PHxgbW1NUJCQhAVFVXn+suXL0f79u1hY2MDpVKJ119/HaWlpY2UlohIP0euZGLU6j+RripFW1d77P5PXxYbogYm6pWb7du3IyIiAqtXr0ZISAiWL1+OgQMHIj4+Hq6u1Sew+uabbzBnzhysX78effr0wdWrVzFhwgRIJBIsW7ZMhCMgIqqZIAjYejIRC364CK0A9PVrji/GBUFhYyV2NCKzJxEEQRDrm4eEhKBnz55YsWIFAECr1UKpVOK1117DnDlzqq0/bdo0XL58GZGRkbplb7zxBk6dOoXjx4/X63uqVCooFArk5+fD0ZETZRGR4eWXlGPungv4MS4VADAqqCXeG94VMkvRL5YTmSx9fn6L9k4rKytDdHQ0QkND/w4jlSI0NBQnTpyocZs+ffogOjpa99HVzZs3sX//fjz11FONkpmI6H5O38rBU58ew49xqbCQSjBncAd89Gw3FhuiRiTax1LZ2dnQaDRwc3OrstzNzQ1XrlypcZvnnnsO2dnZeOSRRyAIAioqKjB58mT83//9X63fR61WQ61W636vUqkMcwBERPeo0Gjx2a/XseLXa9AKQKtmtvh0TAACWzmLHY2oyTGpf0ocPXoU77//Pr744gucPXsWu3btwr59+/Duu+/Wus3ixYuhUCh0L6VS2YiJiagpSM4pxugvT+CzyMpiM7J7S+yf0Y/Fhkgkot1zU1ZWBltbW+zcuRPDhg3TLR8/fjzy8vKwd+/eatv069cPvXr1wpIlS3TLtm7dildeeQWFhYWQSqt3tZqu3CiVSt5zQ0QGsScmBXP3XEChugIO1pZ4b3hXPOPvKXYsIrNjEvfcyGQyBAUFVbk5WKvVIjIyEr17965xm+Li4moFxsLCAkDlyISayOVyODo6VnkRET0sVWk5Zm6LwcztsShUV6CnjzN+ntGPxYbICIg6FDwiIgLjx49Hjx49EBwcjOXLl6OoqAgTJ04EAISHh8PLywuLFy8GAAwZMgTLli1DYGAgQkJCcP36dcybNw9DhgzRlRwiooYWnZiLGdticDu3BBZSCWY80Rb/ebQNnw9FZCRELTdhYWHIysrC/PnzkZ6ejoCAABw4cEB3k3FSUlKVKzVz586FRCLB3LlzkZKSghYtWmDIkCF47733xDoEImpCKjRarDxyA5/9eg0arQBlMxssDwtEkDfvrSEyJqLOcyMGznNDRA/idm4xZm6LxZnEXADA8EAvvDO0MxysOSkfUWPQ5+c3ny1FRHQfP8Sl4n+7z6OgtAL2ckssGtYFwwK9xI5FRLVguSEiqkWhugLz917ArrMpAIDurZzw6ZhAKJvZipyMiOrCckNEVIOYpFzM3B6LxDvFkEqAaY+3xfTH/XjTMJEJYLkhIrqHRitg1dHr+ORw5U3DXk42WD4mAD19mokdjYjqieWGiOgvqXklmLk9FlEJOQCAIf6eWDSsC5/kTWRiWG6IiADsO5eGt3adg6q0AnYyC7wztAtGdPeCRCIROxoR6YnlhoiatCJ1Bd7+4SK+i74NAAhQOuHTMQHwbm4ncjIielAsN0TUZJ2/nY/p22KQkF0EiQSY+qgfZoS2hRVvGiYyaSw3RNQknbmVg/D1USgu08BTYY1PwgIQ0rq52LGIyABYboioyYlNzsOEDadRXKZBX7/m+OK5IChsedMwkblguSGiJuVCSj7C151CoboCvVo3w9rwnrCR8cG7ROaEHywTUZNxOU2F59edgqq0Aj28nbFuPIsNkTliuSGiJuFqRgHGrT2FvOJyBCidsGFiT9jJefGayByx3BCR2buRVYjn1pxCTlEZunopsOnFYD7Nm8iMsdwQkVm7lV2E59acRHahGh09HLFlUjBnHCYycyw3RGS2knOK8dyak8hQqdHOzR5bJwXDyVYmdiwiamAsN0RkllLzSjB2zUmk5peiTQs7fP1SLzS3l4sdi4gaAcsNEZmdDFUpxq45idu5JfBpbotvXu6FFg4sNkRNBcsNEZmVrAI1xq45icQ7xVA2s8E3L/eCm6O12LGIqBGx3BCR2bhTqMa4tSdxM6sIngprfPNSL3g62Ygdi4gaGcsNEZmF3KIyjFt7ClczCuHmKMe3r/SCspmt2LGISAQsN0Rk8vJLyvHC+lO4kl4AF3s5vnm5F7yb24kdi4hEwnJDRCatoLQc49dH4UKKCs3tZPj25RC0aWEvdiwiEhHLDRGZrCJ1BSZuOI3Y5Dw42Vph60shaOvmIHYsIhIZyw0RmaSSMg1e3HgaZxJz4Whtia2TQtDRw1HsWERkBFhuiMjklJZr8PLmMziVkAN7uSU2TwpBFy+F2LGIyEiw3BCRSVFXaDB5azSOX8+GrcwCm17siQClk9ixiMiIsNwQkckoq9Bi6tdncTQ+C9ZWUmyY0BNB3s3EjkVERoblhohMQrlGi+nfxuDw5UzILaVYN74nQlo3FzsWERkhlhsiMnoarYCIHXE4cDEdMgspvnwhCH39XMSORURGiuWGiIyaRivgze/i8GNcKqwsJFj1fHc82t5V7FhEZMRYbojIaGm1Av5v13nsikmBhVSCz8d2xxMd3cSORURGjuWGiIySRivgrV3nsf1MMqQSYHlYAAZ1cRc7FhGZAEuxAxAR/ZO6QoOZ22Lx84V0SCXA0tH+GOLvKXYsIjIRLDdEZFSK1BV4dUvlPDYyCyk+GxuAQV08xI5FRCaE5YaIjEZuURkmbqx8VpStzAJfvdADj7TlqCgi0g/LDREZhfT8Uryw7hSuZRbCydYKGyb0RGArZ7FjEZEJYrkhItHdyi7C8+tO4XZuCdwc5dgyKQTt+HRvInpALDdEJKpLqSqEr49CdqEaPs1tsWVSCJTNbMWORUQmjOWGiERz5lYOJm48jYLSCnT0cMTmF4PRwkEudiwiMnH1Kjc//PBDvXf4zDPPPHAYImo6jsRnYsrWaJSWa9HD2xnrJvSEwsZK7FhEZAbqVW6GDRtW5fcSiQSCIFT5/V0ajcYwyYjIbO2NTcEbO+JQoRXwaPsWWDUuCDYyC7FjEZGZqNcMxVqtVvf65ZdfEBAQgJ9//hl5eXnIy8vD/v370b17dxw4cKCh8xKRidtyMhEzt8eiQivgGX9PfPVCDxYbIjIove+5mTlzJlavXo1HHnlEt2zgwIGwtbXFK6+8gsuXLxs0IBGZB0EQsPLIdXz8y1UAwAu9vLHwmc6QSiX32ZKISD96l5sbN27Aycmp2nKFQoFbt24ZIBIRmRutVsB7+y9j3fEEAMBrj/sh4sl2VT7SJiIyFL0fnNmzZ09EREQgIyNDtywjIwNvvvkmgoODDRqOiExfhUaL/35/Tlds5v27E94Y0J7FhogajN5XbtavX4/hw4ejVatWUCqVAIDk5GS0bdsWe/bsMXQ+IjJhpeUaTP82Br9cyoCFVIIPR3bDs0EtxY5FRGZO73Lj5+eHc+fO4dChQ7hy5QoAoGPHjggNDeW/xIhIp1BdgVc2n8GfN+5AZinFirGBGNDZXexYRNQESIR7x3TrqbS0FHK53KRKjUqlgkKhQH5+PhwdHcWOQ2SWcorKMHFDFOJu58NOZoE143ugTxs+AJOIHpw+P7/1vudGq9Xi3XffhZeXF+zt7ZGQ8Nfn6PPmYd26dQ+WmIjMRlp+CUZ/eQJxt/PhbGuFb1/pxWJDRI1K73KzaNEibNy4ER999BFkMplueZcuXbB27VqDhiMi03IzqxDPrjqB65mFcHe0xneTe6NbSyexYxFRE6N3udm8eTO++uorjBs3DhYWf0+85e/vr7sHh4iangsp+Ri1+gRS8krg62KHnVN6w8+VT/Ymosan9w3FKSkp8PPzq7Zcq9WivLzcIKGIyLREJeRg0sbTKFBXoJOHIzZPCoaLPR+ASUTi0PvKTadOnXDs2LFqy3fu3InAwECDhCIi0/HrlQy8sO4UCtQVCPZphm2v9mKxISJR6X3lZv78+Rg/fjxSUlKg1Wqxa9cuxMfHY/Pmzfjpp58aIiMRGanDlzIweWs0KrQCnujgipXjusPais+JIiJx6X3lZujQofjxxx9x+PBh2NnZYf78+bh8+TJ+/PFHPPnkkw2RkYiM0NmkXEz79qzuAZirXwhisSEio/BQ89yYIs5zQ/TwbmYVYuSqP5FbXI5H27fAmvAesLLQ+99KRET11qDz3LRu3Rp37typtjwvLw+tW7fWd3dEZGIyC0oxfkMUcovL0a2lAiuf685iQ0RGRe+/kW7dugWNRlNtuVqtRkpKikFCEZFxKlRX4MWNp5GcUwLv5rZYP6En7OR637pHRNSg6v230g8//KD79cGDB6FQKHS/12g0iIyMhI+Pj0HDEZHxKNdoMWVrNC6kqNDcToZNEzncm4iMU73LzbBhwwAAEokE48ePr/I1Kysr+Pj4YOnSpQYNR0TGQRAEzP7+HI5dy4aNlQXWTegJHxc7sWMREdWo3uVGq9UCAHx9fXH69Gm4uPBZMURNxce/xGPX2RRYSCVYOS4QAUonsSMREdVK7w/L7z4ok4iahi0nbmHlkRsAgPeHd8HjHdxETkREVDe9byiePn06Pvvss2rLV6xYgZkzZxoiExEZiQMX0jH/h4sAgNdD2yGsZyuRExER3Z/e5eb7779H3759qy3v06cPdu7caZBQRCS+M7dyMGNbDAQBGBusxPQnqj9TjojIGOldbu7cuVNlpNRdjo6OyM7ONkgoIhLX9cwCTNp0BuoKLZ7o4Ip3h3aBRCIROxYRUb3oXW78/Pxw4MCBast//vnnB5rEb+XKlfDx8YG1tTVCQkIQFRVV5/p5eXmYOnUqPDw8IJfL0a5dO+zfv1/v70tENctQlWL8+tPILylHgNIJnz8XCEtO0kdEJkTvG4ojIiIwbdo0ZGVl4fHHHwcAREZGYunSpVi+fLle+9q+fTsiIiKwevVqhISEYPny5Rg4cCDi4+Ph6upabf2ysjI8+eSTcHV1xc6dO+Hl5YXExEQ4OTnpexhEVIOC0nJM2HAaKXkl8HWxw7rxPWAr4yR9RGRaHujZUqtWrcJ7772H1NRUAICPjw/efvtthIeH67WfkJAQ9OzZEytWrABQOdxcqVTitddew5w5c6qtv3r1aixZsgRXrlyBlZWVvrEB8NlSRLUpq9Bi4sYo/HH9DlzsZdg1pS9aNbcVOxYREQD9fn4/1IMzs7KyYGNjA3t7e723LSsrg62tLXbu3KmbIBAAxo8fj7y8POzdu7faNk899RSaNWsGW1tb7N27Fy1atMBzzz2H2bNnw8Kifk8jZrkhqk6rFRCxIxZ7YlNhK7PA9ld6o2vL6vfWERGJpUEfnAkAFRUVOHz4MHbt2oW73Sg1NRWFhYX13kd2djY0Gg3c3KrOmeHm5ob09PQat7l58yZ27twJjUaD/fv3Y968eVi6dCkWLVpU6/dRq9VQqVRVXkRU1YcHr2BPbCospRJ8Ma47iw0RmTS9P0xPTEzEoEGDkJSUBLVajSeffBIODg748MMPoVarsXr16obICaDyYytXV1d89dVXsLCwQFBQEFJSUrBkyRIsWLCgxm0WL16MhQsXNlgmIlO34Y8EfPnbTQDAByO74dH21e93IyIyJXpfuZkxYwZ69OiB3Nxc2NjY6JYPHz4ckZGR9d6Pi4sLLCwskJGRUWV5RkYG3N3da9zGw8MD7dq1q/IRVMeOHZGeno6ysrIat3nrrbeQn5+veyUnJ9c7I5G5238+De/8dAkA8ObA9ng2qKXIiYiIHp7e5ebYsWOYO3cuZDJZleU+Pj5ISUmp935kMhmCgoKqFCKtVovIyEj07t27xm369u2L69ev655zBQBXr16Fh4dHtTx3yeVyODo6VnkREXDq5h3M3B4LQQCe79UK/3m0jdiRiIgMQu9yo9VqodFoqi2/ffs2HBwc9NpXREQE1qxZg02bNuHy5cuYMmUKioqKMHHiRABAeHg43nrrLd36U6ZMQU5ODmbMmIGrV69i3759eP/99zF16lR9D4OoSbuaUYCXN59BWYUWAzq5YeEznKSPiMyH3vfcDBgwAMuXL8dXX30FAJBIJCgsLMSCBQvw1FNP6bWvsLAwZGVlYf78+UhPT0dAQAAOHDigu8k4KSkJUunf/UupVOLgwYN4/fXX0a1bN3h5eWHGjBmYPXu2vodB1GSl5Zdg/PooqEorEOTtjM/GBsJCymJDROZD76Hgt2/fxsCBAyEIAq5du4YePXrg2rVrcHFxwe+//17j5HvGhEPBqSnLLynH6NUnEJ9RgNYt7PD95D5wtqv5I10iImOiz89vva/ctGzZEnFxcdi+fTvi4uJQWFiISZMmYdy4cVVuMCYi46Ku0ODVLWcQn1GAFg5ybJoYzGJDRGapXlduunfvjsjISDg7O+Odd97BrFmzYGtrmjOX8soNNUVarYDp22Lw07k02Mstsf3VXujsyblsiMh0GHwSv8uXL6OoqAgAsHDhQr0m6yMi8b2//zJ+OpcGS6kEq58PYrEhIrNWr4+lAgICMHHiRDzyyCMQBAEff/xxrY9cmD9/vkEDEtHDWXvsJtYeTwAALBnVDY+0dRE5ERFRw6rXx1Lx8fFYsGABbty4gbNnz6JTp06wtKzeiyQSCc6ePdsgQQ2FH0tRUxJ5OQMvbT4DQQDmDO6Ayf05lw0RmaYGfXCmVCpFenq60Y+Kqg3LDTUV1zMLMWzlHyhUV2BcSCssGsa5bIjIdDXoaKl7ZwcmIuOUX1KOVzafQaG6AsE+zbBgSGcWGyJqMvQuNwBw7do1HDlyBJmZmdXKDu+5IRKXRitg5rYY3MwugqfCGl883x0yS70nIyciMll6l5s1a9ZgypQpcHFxgbu7e5V/DUokEpYbIpEt/SUeR+KzILeU4ssXesDFXi52JCKiRqV3uVm0aBHee+89PvKAyAj9dC4VXxy9AQD46Nlu6NqSQ76JqOnR+1p1bm4uRo0a1RBZiOghXEpV4c3vzgEAXvlXawwN8BI5ERGROPQuN6NGjcIvv/zSEFmI6AHlFJXh5c1nUFKuQb+2Lpg9qIPYkYiIRKP3x1J+fn6YN28eTp48ia5du8LKyqrK16dPn26wcER0f+UaLaZ+fRYpeSXwbm6LFWO78ynfRNSk6T3Pja+vb+07k0hw8+bNhw7VkDjPDZmbt3+4iI1/3oKdzAK7p/ZFOzcHsSMRERlcg85zk5CQ8MDBiMiwdpxJxsY/bwEAloUFsNgQEeEB7rkhIuMQk5SLubsvAABmPNEWAzu7i5yIiMg41OvKTUREBN59913Y2dkhIiKiznWXLVtmkGBEVLtMVSkmb41GmUaLAZ3cMOOJtmJHIiIyGvUqNzExMSgvL9f9ujac3p2o4akrNHh1azQyVGq0dbXHsrAASHkDMRGRTr3KzZEjR2r8NRE1LkEQMH/PRcQk5cHR2hJrwnvAXv5AT1EhIjJbvOeGyIRsOZmI7WeSIZUAnz/XHT4udmJHIiIyOiw3RCbi5M07eOfHSwCA2YM6oH+7FiInIiIyTiw3RCbgdm4x/vP1WVRoBTzj74lX/tVa7EhEREaL5YbIyJWUafDqlmjkFJWhs6cjPhzZjTfvExHVgeWGyIgJgoD/fn8OF1NVaG4nw1fhPWAjsxA7FhGRUXugYRbXrl3DkSNHkJmZCa1WW+Vr8+fPN0gwIgK+/P0mfoxLhaVUgi/GdYeXk43YkYiIjJ7e5WbNmjWYMmUKXFxc4O7uXuXyuEQiYbkhMpCj8Zn48MAVAMCCIZ0Q0rq5yImIiEyD3uVm0aJFeO+99zB79uyGyENEABKyi/DatzEQBGBMTyWe7+UtdiQiIpOh9z03ubm5GDVqVENkISIABaXleHnzGRSUVqB7KycsHNqZNxATEelB73IzatQo/PLLLw2RhajJ02oFvL49DtczC+HmKMfq54Mgt+QNxERE+tD7Yyk/Pz/MmzcPJ0+eRNeuXWFlZVXl69OnTzdYOKKmZnnkNRy+nAGZpRRfvtADro7WYkciIjI5EkEQBH028PX1rX1nEglu3rz50KEakkqlgkKhQH5+PhwdHcWOQ6Rz4EIaJm89CwD4eJQ/ng1qKXIiIiLjoc/Pb72v3CQkJDxwMCKqWXx6ASJ2xAEAJvb1YbEhInoIDzWJnyAI0PPCDxH9Q15xGV7efAbFZRr0adMc/3uqo9iRiIhM2gOVm82bN6Nr166wsbGBjY0NunXrhi1bthg6G5HZ02gFvPZtDJJyitHS2QYrnusOSwtOHE5E9DD0/lhq2bJlmDdvHqZNm4a+ffsCAI4fP47JkycjOzsbr7/+usFDEpmr5Yev4ti1bFhbSfHVCz3QzE4mdiQiIpP3QDcUL1y4EOHh4VWWb9q0CW+//bbR35PDG4rJWERezsCkTWcAAMvDAjAs0EvkRERExkufn996X/9OS0tDnz59qi3v06cP0tLS9N0dUZOUdKcYr2+PBQCE9/ZmsSEiMiC9y42fnx927NhRbfn27dvRtm1bg4QiMmel5RpM3hoNVWkFAls5Ye7TncSORERkVvS+52bhwoUICwvD77//rrvn5o8//kBkZGSNpYeI/iYIAubuuYBLaSo0t5Phi3HdIbPkDcRERIak99+qI0eOxKlTp+Di4oI9e/Zgz549cHFxQVRUFIYPH94QGYnMxrdRydgZfRtSCfD52EB4KGzEjkREZHb0vnIDAEFBQdi6dauhsxCZtbjkPLz9w0UAwKyB7dHHz0XkRERE5qle5UalUunuTFapVHWuyxFIRNXlFJVhytZolGm0GNDJDVP6txE7EhGR2apXuXF2dkZaWhpcXV3h5OQEiURSbR1BECCRSKDRaAweksiUabQCZmyLQWp+KXxd7PDxaP8a30NERGQY9So3v/76K5o1awYAOHLkSIMGIjI3907Ut+r57nC0thI7EhGRWatXuenfv7/u176+vlAqldX+5SkIApKTkw2bjsjERV7OwOe/XgcAfDCiGzq482NbIqKGpvdoKV9fX2RlZVVbnpOTA19fX4OEIjIHnKiPiEgcepebu/fW/FNhYSGsra0NEorI1HGiPiIi8dR7KHhERAQAQCKRYN68ebC1tdV9TaPR4NSpUwgICDB4QCJTw4n6iIjEVe9yExMTA6DyL+7z589DJvv76cUymQz+/v6YNWuW4RMSmRhO1EdEJK56l5u7o6QmTpyITz/9lPPZENXg3on63hzYgRP1ERGJQO9r5cuXL0dFRUW15Tk5Ofed4I/InP1zor7J/VuLHYmIqEnSu9yMGTMG27Ztq7Z8x44dGDNmjEFCEZkaTtRHRGQ89C43p06dwmOPPVZt+aOPPopTp04ZJBSRqeFEfURExkPvcqNWq2v8WKq8vBwlJSUGCUVkSjhRHxGRcdG73AQHB+Orr76qtnz16tUICgoySCgiU3HvRH3jOVEfEZFRqPdoqbsWLVqE0NBQxMXF4YknngAAREZG4vTp0/jll18MHpDIWN07UV/3Vk74HyfqIyIyCnpfuenbty9OnDgBpVKJHTt24Mcff4Sfnx/OnTuHfv36NURGIqMjCAL+t/vvifpWcqI+IiKjofeVGwAICAjA119/begsRCbjm6gkfH+WE/URERmjByo3d5WWlqKsrKzKMk7uR+YuLjkPC3+4BIAT9RERGSO9r6MXFxdj2rRpcHV1hZ2dHZydnau8iMwZJ+ojIjJ+epebN998E7/++itWrVoFuVyOtWvXYuHChfD09MTmzZsbIiORUeBEfUREpkHvj6V+/PFHbN68GY8++igmTpyIfv36wc/PD97e3vj6668xbty4hshJJLq7E/XZWFlwoj4iIiOm95WbnJwctG5deSne0dEROTk5AIBHHnkEv//+u2HTERmJKhP1jezKifqIiIyY3uWmdevWSEhIAAB06NABO3bsAFB5RcfJycmg4YiMwT8n6hsawIn6iIiMmd7lZuLEiYiLiwMAzJkzBytXroS1tTVef/11vPnmmwYPSCSmCo0W07fFcKI+IiITIhEEQXiYHSQmJiI6Ohp+fn7o1q2boXI1GJVKBYVCgfz8fA5bp/v69PA1fHL4KhysLXFg5r/g5cT5bIiIxKDPz2+9rtyUl5fjiSeewLVr13TLvL29MWLECJMoNkT6iE3Ow2e/Vv5ZXzSsC4sNEZGJ0KvcWFlZ4dy5cwYPsXLlSvj4+MDa2hohISGIioqq13bbtm2DRCLBsGHDDJ6Jmrbisgq8vj0WGq2AIf6evM+GiMiE6H3PzfPPP49169YZLMD27dsRERGBBQsW4OzZs/D398fAgQORmZlZ53a3bt3CrFmz+DwrahDv7buMhOwiuDtaY9HQLmLHISIiPeg9z01FRQXWr1+Pw4cPIygoCHZ2dlW+vmzZMr32t2zZMrz88suYOHEiAGD16tXYt28f1q9fjzlz5tS4jUajwbhx47Bw4UIcO3YMeXl5+h4GUa2OXMnE16eSAABLR/tDYcv5bIiITIne5ebChQvo3r07AODq1atVvqbvbK1lZWWIjo7GW2+9pVsmlUoRGhqKEydO1LrdO++8A1dXV0yaNAnHjh2r83uo1Wqo1Wrd71UqlV4ZqWm5U6jGmzsrP3p9sa8v+vK5UUREJqde5ebcuXPo0qULpFIpjhw5YrBvnp2dDY1GAzc3tyrL3dzccOXKlRq3OX78ONatW4fY2Nh6fY/Fixdj4cKFDxuVmgBBEPDWrvPILlSjnZs9/juovdiRiIjoAdTrnpvAwEBkZ2cDqJzE786dOw0aqjYFBQV44YUXsGbNGri41O9f1G+99Rby8/N1r+Tk5AZOSabquzO38culDFhZSPBJWACsrSzEjkRERA+gXldunJyckJCQAFdXV9y6dQtardYg39zFxQUWFhbIyMiosjwjIwPu7u7V1r9x4wZu3bqFIUOG6JbdzWJpaYn4+Hi0adOmyjZyuRxyudwgecl8Jd0pxsIfLwIAIp5sj86eCpETERHRg6pXuRk5ciT69+8PDw8PSCQS9OjRAxYWNf+r9ubNm/X+5jKZDEFBQYiMjNQN59ZqtYiMjMS0adOqrd+hQwecP3++yrK5c+eioKAAn376KZRKZb2/N9FdGq2AiB2xKCrTINinGV75V2uxIxER0UOoV7n56quvMGLECFy/fh3Tp0/Hyy+/DAcHB4MEiIiIwPjx49GjRw8EBwdj+fLlKCoq0o2eCg8Ph5eXFxYvXgxra2t06VJ1WO7d51n9czlRfa3+7QbOJObCXm6JpaP9YSHV78Z4IiIyLvUeLTVo0CAAQHR0NGbMmGGwchMWFoasrCzMnz8f6enpCAgIwIEDB3Q3GSclJUEq1Xs6HqJ6OX87H58cqhz1t/CZzlA2sxU5ERERPayHfraUqeGzpeiu0nINnv7sGG5kFWFwF3d8Ma673tMZEBFR42iwZ0sRmZMPfr6CG1lFcHWQ4/3hXVlsiIjMBMsNNUm/X83Cxj9vAQCWjPKHs51M3EBERGQwLDfU5OQWlWHWd3EAgPG9vdG/XQuRExERkSGx3FCTIggC5u65gMwCNdq0sMOcwR3FjkRERAbGckNNyu6YFOw7nwZLqQTLwwJhI+MsxERE5oblhpqM27nFWLC3chbimaFt0bUlZyEmIjJHLDfUJFTOQhyHAnUFgrydMbl/m/tvREREJonlhpqEtcduIiohB3YyCywb7Q9LC/7RJyIyV/wbnszepVQVPv4lHgAwf0gneDe3EzkRERE1JJYbMmul5RrM3B6Dco2AAZ3cMLoHH65KRGTuWG7IrH18MB5XMwrhYi/D4hGchZiIqClguSGz9ef1bKw9ngAA+OjZbmhuLxc5ERERNQaWGzJL+cXleOOvWYifC2mFxzu4iZyIiIgaC8sNmaV5ey8gLb8Uvi52mPs0ZyEmImpKWG7I7OyNTcEPcamwkEqwbLQ/bGWWYkciIqJGxHJDZiU1rwTz9lwAAEx7zA+BrZxFTkRERI2N5YbMhlYrYNZ3cVCVVsBf6YRpj/uJHYmIiETAckNmY/0fCfjzxh3YWFlgeVgArDgLMRFRk8S//cksxKcX4KODlbMQz/13R/i6cBZiIqKmiuWGTJ66QoOZ22NRVqHF4x1c8VxwK7EjERGRiFhuyOR9fDAel9NUaGYnwwcjOQsxEVFTx3JDJu3YtSysOVY5C/GHI7vB1cFa5ERERCQ2lhsyWTlFZXhjR+UsxONCWuHJTpyFmIiIWG7IRAmCgNnfn0NmgRp+rvaY+3QnsSMREZGRYLkhk/RNVBIOXcqAzEKKT8cEwEZmIXYkIiIyEiw3ZHKuZxbg3Z8uAQD+O6g9OnsqRE5ERETGhOWGTIq6QoPp38aitFyLfm1d8GJfX7EjERGRkWG5IZPy8cF4XPpr2PfSUf6QSjnsm4iIqmK5IZNRbdi3I4d9ExFRdSw3ZBLuHfb9fC8O+yYiotqx3JDREwQB/93597Dv/z3FYd9ERFQ7lhsyel+fSsLhyxz2TURE9cNyQ0btemYBFu3jsG8iIqo/lhsyWuoKDV7jsG8iItITyw0ZrSUH/n7aN4d9ExFRfbHckFH6/WoW1h7nsG8iItIfyw0ZnTuFarzxHYd9ExHRg2G5IaNy92nfWRz2TURED4jlhoxK5bDvTMgspPhsTCCHfRMRkd5Ybsho/HPYdydPR5ETERGRKWK5IaPAYd9ERGQoLDdkFDjsm4iIDIXlhkR377Dvjzjsm4iIHhLLDYnqn8O+Qznsm4iIHhLLDYmGw76JiKghsNyQaLZy2DcRETUAlhsSxbWMAiz6icO+iYjI8FhuqNGpKzSYvi0W6goO+yYiIsNjuaFG9xGHfRMRUQNiuaFG9fvVLKz7a9j3kmc57JuIiAyP5YYazb3Dvl/o5Y0nOnLYNxERGR7LDTUKQRDw352Vw77butrjf093FDsSERGZKZYbahRbTyYi8spfw77HBsLaisO+iYioYbDcUIOLSsjBu/suAwBmD+6Ajh4c9k1ERA2H5YYa1NWMAry06TTKKrQY0MkNE/v4iB2JiIjMHMsNNZi0/BKMXx8FVWkFgryd8dnYQA77JiKiBsdyQw0iv6QcE9afRlp+KVq3sMPa8B68z4aIiBoFyw0ZnLpCg1e3nEF8RgFaOMixaWIwnO1kYsciIqImguWGDEqrFfDGjjicvJkDe7klNk7sCWUzW7FjERFRE8JyQwb1/v7L+OlcGiylEqx+PgidPRViRyIioiaG5YYMZu2xm1h799EKo7rhkbYuIiciIqKmiOWGDOKHuFQs+msumzmDO2B4YEuRExERUVPFckMP7c8b2Zi1o/KZURP6+ODVf7UWORERETVlLDf0UK6kq/Dq5miUabQY3MUd8/7dCRIJ57IhIiLxsNzQA0vJq5ykr0BdgWCfZvgkLAAWnKSPiIhExnJDDyS/uBwT1kchQ1X5lO81nKSPiIiMhFGUm5UrV8LHxwfW1tYICQlBVFRUreuuWbMG/fr1g7OzM5ydnREaGlrn+mR4peUavLz5DK5lFsLNUY6NLwZDYWsldiwiIiIARlButm/fjoiICCxYsABnz56Fv78/Bg4ciMzMzBrXP3r0KMaOHYsjR47gxIkTUCqVGDBgAFJSUho5edOk0QqI2BGLqFs5cJBbYuPEYHg52Ygdi4iISEciCIIgZoCQkBD07NkTK1asAABotVoolUq89tprmDNnzn2312g0cHZ2xooVKxAeHn7f9VUqFRQKBfLz8+Ho6PjQ+ZsSQRCw8MdL2PjnLVhZSLDpxWD0acO5bIiIqOHp8/Nb1Cs3ZWVliI6ORmhoqG6ZVCpFaGgoTpw4Ua99FBcXo7y8HM2aNavx62q1GiqVqsqLHsyaYzex8c9bAIClowNYbIiIyCiJWm6ys7Oh0Wjg5uZWZbmbmxvS09PrtY/Zs2fD09OzSkG61+LFi6FQKHQvpVL50Lmbor2xKXh//xUAwNynO+IZf0+RExEREdVM9HtuHsYHH3yAbdu2Yffu3bC2tq5xnbfeegv5+fm6V3JyciOnNH1/XM/GrO8qJ+mb9IgvXurHSfqIiMh4WYr5zV1cXGBhYYGMjIwqyzMyMuDu7l7nth9//DE++OADHD58GN26dat1PblcDrlcbpC8TdGlVBVe3RKNco2Ap7t54H9PdRQ7EhERUZ1EvXIjk8kQFBSEyMhI3TKtVovIyEj07t271u0++ugjvPvuuzhw4AB69OjRGFGbpNu5xZiwIQqF6gr0at0My0b7Q8pJ+oiIyMiJeuUGACIiIjB+/Hj06NEDwcHBWL58OYqKijBx4kQAQHh4OLy8vLB48WIAwIcffoj58+fjm2++gY+Pj+7eHHt7e9jb24t2HOYmr7gM49dHIbNAjfZuDvjyhR6QW3KSPiIiMn6il5uwsDBkZWVh/vz5SE9PR0BAAA4cOKC7yTgpKQlS6d8XmFatWoWysjI8++yzVfazYMECvP32240Z3WyVlmvw0qYzuJFVBA+FNTa+2BMKG07SR0REpkH0eW4aG+e5qZtGK+A/X0fj4MUMOFpbYueUPmjn5iB2LCIiauJMZp4bMi6CIODtHy7i4MUMyCykWBPeg8WGiIhMDssN6az67Qa2nEyERAJ8EhaAkNbNxY5ERESkN9HvuSHxCYKAraeS8NGBeADA/H93wtPdPERORURE9GBYbpq4gtJyzN1zAXtjUwEAr/yrNSb29RU5FRER0YNjuWnCzt3Ow2vfxiDxTjEspBJEPNkOU/q3ETsWERHRQ2G5aYK0WgHrjifgo4NXUK4R4OVkg8/GBiDIu+aHjxIREZkSlpsmJrtQjVnfxeFofBYAYHAXd3wwohsUtpzHhoiIzAPLTRPyx/VszNwei6wCNeSWUsz7dyeMC2kFiYSPVCAiIvPBctMEVGi0+OTwVXxx9AYEAfBztceK5wLRwZ2TGBIRkflhuTFzt3OLMWNbLKITcwEAY4OVmP/vzrCR8TlRRERknlhuzNjP59Mw+/tzUJVWwEFuicUju+Lf3TzFjkVERNSgWG7MUGm5Bu/+dAlfn0oCAAQonfD52EAom9mKnIyIiKjhsdyYmWsZBZj2TQziMwoAAJP7t8EbA9rByoJP2iAioqaB5cZMCIKA7aeT8faPF1FaroWLvQzLRgfgX+1aiB2NiIioUbHcmAFVaTne2nUe+86lAQD6tXXBstEBaOEgFzkZERFR42O5MXExSbmYvi0GyTklsJRKMGtge7zSrzWkUs5dQ0RETRPLjYnSagV8dewmPj4YjwqtgJbONvh8bCACWzmLHY2IiEhULDcmKKtAjYgdsTh2LRsA8HQ3Dywe0RWO1nyEAhEREcuNiTl2LQuvb49DdqEa1lZSvD2kM8J6KvkIBSIior+w3JiI0nINlh++htW/3QAAtHdzwIrnAtHWzUHkZERERMaF5cbIFZSWY+vJJKw7noDsQjUAYFxIK8z7dydYW/ERCkRERP/EcmOksgvV2PBHAjafSERBaQUAwFNhjXn/7oTBXT1ETkdERGS8WG6MzO3cYqz5/Sa2nU6GukILoPIp3pP7t8HQAE/ONExERHQfLDdG4mpGAVYfvYG9canQaAUAgL/SCf95tA2e7OjGeWuIiIjqieVGZDFJufji6A0cupShW/aInwv+82gb9G7TnKOgiIiI9MRyIwJBEHD8eja+OHIDJ27eAQBIJMCgzu6Y3L8N/JVO4gYkIiIyYSw3jUijFXDwYjpWHb2B8yn5AABLqQTDA73wav828HO1FzkhERGR6WO5aQRlFVrsiUnB6t9u4GZ2EQDAxsoCY4KVeLlfa3g62YickIiIyHyw3DSgInUFvo1KwtpjCUhXlQIAFDZWGN/HBxP6+KCZnUzkhEREROaH5aYB5BWXYeOft7Dxz1vIKy4HALg5yvHSI60xNqQV7OX8305ERNRQ+FPWgNLzS7H22E18E5WE4jINAMCnuS1e7d8GI7p7QW7JGYWJiIgaGsuNgfx8Pg3Tt8WgXFM5R00nD0f857E2GNzFAxaco4aIiKjRsNwYSJCPMyQSCYJ9nfGfR9ugf7sWnKOGiIhIBCw3BuLqYI3IiP5QNrMVOwoREVGTxgcVGRCLDRERkfhYboiIiMissNwQERGRWWG5ISIiIrPCckNERERmheWGiIiIzArLDREREZkVlhsiIiIyKyw3REREZFZYboiIiMissNwQERGRWWG5ISIiIrPCckNERERmheWGiIiIzIql2AEamyAIAACVSiVyEiIiIqqvuz+37/4cr0uTKzcFBQUAAKVSKXISIiIi0ldBQQEUCkWd60iE+lQgM6LVapGamgoHBwdIJBKD7lulUkGpVCI5ORmOjo4G3bex4bGar6Z0vDxW89WUjrepHKsgCCgoKICnpyek0rrvqmlyV26kUilatmzZoN/D0dHRrP+A3YvHar6a0vHyWM1XUzrepnCs97ticxdvKCYiIiKzwnJDREREZoXlxoDkcjkWLFgAuVwudpQGx2M1X03peHms5qspHW9TOtb6anI3FBMREZF545UbIiIiMissN0RERGRWWG6IiIjIrLDcEBERkVlhudHTypUr4ePjA2tra4SEhCAqKqrO9b/77jt06NAB1tbW6Nq1K/bv399ISR/c4sWL0bNnTzg4OMDV1RXDhg1DfHx8ndts3LgREomkysva2rqREj+ct99+u1r2Dh061LmNKZ5XAPDx8al2rBKJBFOnTq1xfVM6r7///juGDBkCT09PSCQS7Nmzp8rXBUHA/Pnz4eHhARsbG4SGhuLatWv33a++7/nGUtfxlpeXY/bs2ejatSvs7Ozg6emJ8PBwpKam1rnPB3kvNIb7ndsJEyZUyz1o0KD77tcYz+39jrWm969EIsGSJUtq3aexnteGxHKjh+3btyMiIgILFizA2bNn4e/vj4EDByIzM7PG9f/880+MHTsWkyZNQkxMDIYNG4Zhw4bhwoULjZxcP7/99humTp2KkydP4tChQygvL8eAAQNQVFRU53aOjo5IS0vTvRITExsp8cPr3LlzlezHjx+vdV1TPa8AcPr06SrHeejQIQDAqFGjat3GVM5rUVER/P39sXLlyhq//tFHH+Gzzz7D6tWrcerUKdjZ2WHgwIEoLS2tdZ/6vucbU13HW1xcjLNnz2LevHk4e/Ysdu3ahfj4eDzzzDP33a8+74XGcr9zCwCDBg2qkvvbb7+tc5/Gem7vd6z3HmNaWhrWr18PiUSCkSNH1rlfYzyvDUqgegsODhamTp2q+71GoxE8PT2FxYsX17j+6NGjhaeffrrKspCQEOHVV19t0JyGlpmZKQAQfvvtt1rX2bBhg6BQKBovlAEtWLBA8Pf3r/f65nJeBUEQZsyYIbRp00bQarU1ft1UzysAYffu3brfa7Vawd3dXViyZIluWV5eniCXy4Vvv/221v3o+54Xyz+PtyZRUVECACExMbHWdfR9L4ihpmMdP368MHToUL32Ywrntj7ndejQocLjjz9e5zqmcF4NjVdu6qmsrAzR0dEIDQ3VLZNKpQgNDcWJEydq3ObEiRNV1geAgQMH1rq+scrPzwcANGvWrM71CgsL4e3tDaVSiaFDh+LixYuNEc8grl27Bk9PT7Ru3Rrjxo1DUlJSreuay3ktKyvD1q1b8eKLL9b5EFlTPq93JSQkID09vcp5UygUCAkJqfW8Pch73pjl5+dDIpHAycmpzvX0eS8Yk6NHj8LV1RXt27fHlClTcOfOnVrXNZdzm5GRgX379mHSpEn3XddUz+uDYrmpp+zsbGg0Gri5uVVZ7ubmhvT09Bq3SU9P12t9Y6TVajFz5kz07dsXXbp0qXW99u3bY/369di7dy+2bt0KrVaLPn364Pbt242Y9sGEhIRg48aNOHDgAFatWoWEhAT069cPBQUFNa5vDucVAPbs2YO8vDxMmDCh1nVM+bze6+650ee8Pch73liVlpZi9uzZGDt2bJ0PVtT3vWAsBg0ahM2bNyMyMhIffvghfvvtNwwePBgajabG9c3l3G7atAkODg4YMWJEneuZ6nl9GE3uqeCkn6lTp+LChQv3/Xy2d+/e6N27t+73ffr0QceOHfHll1/i3XffbeiYD2Xw4MG6X3fr1g0hISHw9vbGjh076vUvIlO1bt06DB48GJ6enrWuY8rnlSqVl5dj9OjREAQBq1atqnNdU30vjBkzRvfrrl27olu3bmjTpg2OHj2KJ554QsRkDWv9+vUYN27cfW/yN9Xz+jB45aaeXFxcYGFhgYyMjCrLMzIy4O7uXuM27u7ueq1vbKZNm4affvoJR44cQcuWLfXa1srKCoGBgbh+/XoDpWs4Tk5OaNeuXa3ZTf28AkBiYiIOHz6Ml156Sa/tTPW83j03+py3B3nPG5u7xSYxMRGHDh2q86pNTe73XjBWrVu3houLS625zeHcHjt2DPHx8Xq/hwHTPa/6YLmpJ5lMhqCgIERGRuqWabVaREZGVvmX7b169+5dZX0AOHToUK3rGwtBEDBt2jTs3r0bv/76K3x9ffXeh0ajwfnz5+Hh4dEACRtWYWEhbty4UWt2Uz2v99qwYQNcXV3x9NNP67WdqZ5XX19fuLu7VzlvKpUKp06dqvW8Pch73pjcLTbXrl3D4cOH0bx5c733cb/3grG6ffs27ty5U2tuUz+3QOWV16CgIPj7++u9rameV72IfUezKdm2bZsgl8uFjRs3CpcuXRJeeeUVwcnJSUhPTxcEQRBeeOEFYc6cObr1//jjD8HS0lL4+OOPhcuXLwsLFiwQrKyshPPnz4t1CPUyZcoUQaFQCEePHhXS0tJ0r+LiYt06/zzWhQsXCgcPHhRu3LghREdHC2PGjBGsra2FixcvinEIennjjTeEo0ePCgkJCcIff/whhIaGCi4uLkJmZqYgCOZzXu/SaDRCq1athNmzZ1f7mimf14KCAiEmJkaIiYkRAAjLli0TYmJidKODPvjgA8HJyUnYu3evcO7cOWHo0KGCr6+vUFJSotvH448/Lnz++ee639/vPS+muo63rKxMeOaZZ4SWLVsKsbGxVd7HarVat49/Hu/93gtiqetYCwoKhFmzZgknTpwQEhIShMOHDwvdu3cX2rZtK5SWlur2YSrn9n5/jgVBEPLz8wVbW1th1apVNe7DVM5rQ2K50dPnn38utGrVSpDJZEJwcLBw8uRJ3df69+8vjB8/vsr6O3bsENq1ayfIZDKhc+fOwr59+xo5sf4A1PjasGGDbp1/HuvMmTN1/1/c3NyEp556Sjh79mzjh38AYWFhgoeHhyCTyQQvLy8hLCxMuH79uu7r5nJe7zp48KAAQIiPj6/2NVM+r0eOHKnxz+3d49FqtcK8efMENzc3QS6XC0888US1/wfe3t7CggULqiyr6z0vprqONyEhodb38ZEjR3T7+Ofx3u+9IJa6jrW4uFgYMGCA0KJFC8HKykrw9vYWXn755WolxVTO7f3+HAuCIHz55ZeCjY2NkJeXV+M+TOW8NiSJIAhCg14aIiIiImpEvOeGiIiIzArLDREREZkVlhsiIiIyKyw3REREZFZYboiIiMissNwQERGRWWG5ISIiIrPCckNEJuHo0aOQSCTIy8sTOwoRGTlO4kdERunRRx9FQEAAli9fDgAoKytDTk4O3NzcIJFIxA1HREbNUuwARET1IZPJTOaJzUQkLn4sRURGZ8KECfjtt9/w6aefQiKRQCKRYOPGjVU+ltq4cSOcnJzw008/oX379rC1tcWzzz6L4uJibNq0CT4+PnB2dsb06dOh0Wh0+1ar1Zg1axa8vLxgZ2eHkJAQHD16VJwDJaIGwSs3RGR0Pv30U1y9ehVdunTBO++8AwC4ePFitfWKi4vx2WefYdu2bSgoKMCIESMwfPhwODk5Yf/+/bh58yZGjhyJvn37IiwsDAAwbdo0XLp0Cdu2bYOnpyd2796NQYMG4fz582jbtm2jHicRNQyWGyIyOgqFAjKZDLa2trqPoq5cuVJtvfLycqxatQpt2rQBADz77LPYsmULMjIyYG9vj06dOuGxxx7DkSNHEBYWhqSkJGzYsAFJSUnw9PQEAMyaNQsHDhzAhg0b8P777zfeQRJRg2G5ISKTZWtrqys2AODm5gYfHx/Y29tXWZaZmQkAOH/+PDQaDdq1a1dlP2q1Gs2bN2+c0ETU4FhuiMhkWVlZVfm9RCKpcZlWqwUAFBYWwsLCAtHR0bCwsKiy3r2FiIhMG8sNERklmUxW5UZgQwgMDIRGo0FmZib69etn0H0TkfHgaCkiMko+Pj44deoUbt26hezsbN3Vl4fRrl07jBs3DuHh4di1axcSEhIQFRWFxYsXY9++fQZITUTGgOWGiIzSrFmzYGFhgU6dOqFFixZISkoyyH43bNiA8PBwvPHGG2jfvj2GDRuG06dPo1WrVgbZPxGJjzMUExERkVnhlRsiIiIyKyw3REREZFZYboiIiMissNwQERGRWWG5ISIiIrPCckNERERmheWGiIiIzArLDREREZkVlhsiIiIyKyw3REREZFZYboiIiMissNwQERGRWfl/c6XS4ZXxGTYAAAAASUVORK5CYII=",
      "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": 15,
   "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 0x12c8eea50>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAjrpJREFUeJzs/Xmcpedd33l/7v0+99lP7VVd6m5JlmxjbNmyLWQDwSDiMGCGJMzjSQImnpAMPMAQK5lgJ2CGIYMzk4Hx87zsxIMJA1kAZzxAFowxKDhgMBjb8YK1S92q7q7tVJ39nHu/7/nj6LpU1ZJstaxWb7/361UvSdV1qk61uvt8+7p+i1GWZYkQQgghxBViXuknIIQQQogbm4QRIYQQQlxREkaEEEIIcUVJGBFCCCHEFSVhRAghhBBXlIQRIYQQQlxREkaEEEIIcUVJGBFCCCHEFWVf6SfwXBRFwfb2NvV6HcMwrvTTEUIIIcRzUJYl4/GY9fV1TPPZzz+uiTCyvb3N5ubmlX4aQgghhHgezp07x4kTJ571x6+JMFKv14H5N9NoNK7wsxFCCCHEczEajdjc3NSv48/mmggj6mqm0WhIGBFCCCGuMV+pxEIKWIUQQghxRUkYEUIIIcQVdclh5A/+4A94y1vewvr6OoZh8Ju/+Ztf8TEf//jHec1rXoPnedx666380i/90vN4qkIIIYS4Hl1yGJlOp7zqVa/i/e9//3P6+DNnzvDt3/7tvOlNb+Jzn/scf/fv/l2+//u/n9/5nd+55CcrhBBCiOvPJRewftu3fRvf9m3f9pw//gMf+ACnT5/mZ3/2ZwF42ctexic+8Qn+j//j/+DNb37zpX55IYQQQlxnLnvNyCc/+UnuueeeY+9785vfzCc/+cnL/aWFEEIIcQ247K29u7u7rKysHHvfysoKo9GIMAypVCpPe0wcx8RxrP97NBpd7qcphBBCiCvkquymec973kOz2dRvMn1VCCGEuH5d9jCyurrK3t7esfft7e3RaDSe8VQE4F3vehfD4VC/nTt37nI/TSGEEEJcIZf9mubuu+/mIx/5yLH3/e7v/i533333sz7G8zw8z7vcT00IIYQQV4FLPhmZTCZ87nOf43Of+xwwb9393Oc+x9bWFjA/1Xjb296mP/4HfuAHePzxx/kH/+Af8OCDD/LP/tk/49/+23/LO97xjhfmOxBCCCHENe2Sw8inP/1pXv3qV/PqV78agHvvvZdXv/rVvPvd7wZgZ2dHBxOA06dP81u/9Vv87u/+Lq961av42Z/9WX7hF35B2nqFEEIIAYBRlmV5pZ/EVzIajWg2mwyHQ1mUJ4QQQnwVyrKkLEuKojj2757nfcWFdpfqub5+XxNbe4UQQgjxlalwcTRsXPzPZzuDcBwHy7Je5Gc8J2FECCGEuEY8W8A4GjRU2Lj434++zzAM/ZgkSYiiiJtuuknCiBBCCHEjuzhYZFlGnucURUGWZc85bByl/jtNU8IwZDqdkiQJk8mE0WjEbDZjPB4ThiFve9vbWF5efhG/46dIGBFCCCEus4vrM1TYUEEjz3P9Y+rHjz7u6D+PUp8zSRLyPCeOY6IoYjKZEEURg8GA2WzGbDYjDENms5kOJkmSkKap/vpvetObJIwIIYQQ16JnuzZRQUOdaqhTDvUGT4WNo2FEfVye5zpgqLBRlqUOE2EYEkUR4/GYJEn0mpXpdKpDRxzHpGn61NcxbRK7RmJXyb02qVsn8Wtkbp3dUfwVvtPLR8KIEEII8RypgKBOFNTpQpqmOkSoUw/g2EnHxdcqeZ4f+7xlWRKGIXEc66ChAoZ6O3rCEUURaZoSRdH8eRUFqelTVFqkTp3Mq5O1G2ROfR46nBq55T/r9zbMrlwkkDAihBBCPEkFhosDh1rgevHVysXXKoAuDjUM41j4OHraocJLFEU6YIxGI33KMZvNmEwm+jqlKAriLCexqswMn9TpkFUapK06qTMPGoldpTS+cgGqQ0aViCoJNTOhbmU07ZyXn3jpZft5/UokjAghhLihHL0GUWHjaOhQQeFo2FAMw8A05/NCVeBwHEf/+NFTkaOhQ121TKdTJpMJ4/GYyWSig8h0OiWKIpI0JSwsplaDpLJIFLSI2m0ip0liPfM+t6MMSgIjpWok1IzkyX/GNOycppPTdudhBMA0TQzDwDAMLMtisX7l1rBIGBFCCHHdSpLk2Js67VBvR4OJOsU4+gLteR6maR4LIBefdqjPozpVVMBQYWM6nTIej5nNZjqURHFMWLqMzSqh0yRyN4kWOsycBpn57KHAIqdmJNTMlLo5Dx11cx48AmLqVoZtmdi2rYPTU282lUpF73/zfR/P87BtG8uyqFarl/9/yLOQMCKEEOK6oIo71Qu+ulZRdRXqiqUsSyzL0i/Snufhui6m+dSLuDrhOFr7oT7X0boNVSyqrlRms5me2xGGIUmaMi0cJladqdVhZjeZ1ZvM2g0y03m274SGmdKxYjpWTNOY0TRCakaMU6SY5jwoOY6jn7PjuHheHdu29ZvrugRBgG3b+L6vP159z0fffN8nCIIX73/WRSSMCCGEuCap8KGCh+o4UW9hGFIUBaZp4routVpNvyBblnUsdKjrlNlsduz65uK2WBU0VBeL+rG8KBimFmOzxtRqMbHqTKo1Jmad3Hjml1qDkpaV0rFjOnZCy4hoMKWaT7CNUocG27bxPJ9KpX0sWFiWpUOJZVn6pMNxHFzX1f9UYUudgKhTnjzP9c/daDSSMCKEEEJ8JWqextE3VWyqAoS6QjFNk2q1qus51ImAqt9QIUYFEXXioU80nvyxLMv0+yaTyXx+Rxwzzm0GZpOhvcbQaTM0G2SVZy4eNSlpmvNTjo6T0jJCWsaMWhlCkenTDdM0cRyHIFjC9319peK6rg4gjuPg+z6WZenAoU431OdQIQWeqltRwero9dTFA9KSJMG2r0wskDAihBDiqqSuRtSph5qXARyr91DXEqZpHutiUQHlaABRJx4qYBRFof+p3n90jsdsNmOawWFZY2i1GFin6Pl1YsN92vM1KWiZMW19vRLSNGbUyxDKQj/H+UmHh2HUn3ZyUa1WqdfrBEFApVLRAcR1XVzX1UFDhSvLsvT3qd6ODjK7eEiaqiNRn+Piz3+lSBgRQghxVSiK4ti1i6rxgKe6VFS9h6qJUC++SZIA6BMNdVqiAkYcx8f+++LuGRVCstJgJ7I5KKsclsv0jDoT9+ldLAYlC+aMRXPGAmM6jGkaEbZlUpYlju08GSY8TLOiQ4RhGHieR71ep1arEQQBtVoN3/f1CYgKLerNtu1j4+HVz83RgWrq50kV36rHAfrnSoUOFd6u1B6aZyJhRAghxBWjAoKaGHr0b/JHh4Wpug/DMEjTlOl0euyKRYWNozNCjp6AHC1iVUGlxKAbW2zHLntZk25RZVhWKQ1j/gSMp55nw4hYNCYsmlMWmLBghVjMn6s6VfD9uq7pUC2/tm1Tq9Wo1WpUKhX9pk5DVL2HarM9OjRNXUWpUKZ+TIUU1fGjTjrUKcfRsKGCzdVOwogQQogXlQoLR08sYP63ekAHkqMvolmW6Y4VFTTCMDw2yVQFE9Viqz53nueYpsW4cDg389iJPXZTj8OySs5FpwMG+CQsGhOWzClL5oxlJyKwn5qSOm/5rVKpVKjXn7pq8X1fn3DUajVdTDrvdnkqJKhrpKNXSc90tXLx6YgKGCpkqACiAof6+bsWSRgRQghx2akXXRVAVMsscOwq5mihqarZUEvfVO1IlmVPqws5WgOi5obMjAqPTzzORh7beZ2wfHorrU1GpxyzZM1YsSNW3ZimUwJPjW637fkJRr1ep9Fo6BOIer2uw4c65ThaAHr0uuXiVmEVbNTb0auVowHk6GnH0ROR642EESGEEJfFxTUgaZoeO8lQ4UNdT8RxrFfaq3Chgot68XYcR1/tHD1JSJKESWrwROSyFQecT6uMy+PDw4yyoMWERWPCshWy7qcs+yWmOX9xnweCp9pkVU3HU9cwvm4PfqaCWRUSjl6NPFt4mBewGscGk118xXIjkTAihBDiBaMCyNGNsepEQL342rb9VIB4sl322HTSJxe/HT0xyPMc27Z1oWkcxwxnCVuhw7kk4EK2SL88XmhqUNIpR6waIzbdGZtBRq0yDwHz+hNP13Ycnc8RBMGxmo5KpYLrujpIHb1OUlclRwtD1fcJT01pPTp6HdD1HhfP/rhRSRgRQgjxVVEB5GjRZRRFuuDyaOeLunIZjUZP7WO5qG3X87xj7avqZOWgP+Shw5TtrMZuucBBWaXk+KlDmymrxpBVc8TJIKVdqzx50tHA9/1jxbDqBEQVnapAocKJeg5Hu1XU9YlqvVUfd3QR3jPtsjk6/0PCx9NJGBFCCHHJjo5HV3M51JXK0b/1J0lCv9/XW2mPTjJVJx/qpELVW6jumsFwxBOjgnNpwG7ZpFuuknP8RbxuRGzYE1aNERvOlIZ7dFS6S6VSIQgCPZ1Ufa2j3ScXTyw1TVNfAamTE3ViUqlUsG372AyTNE2f9vOjwoc6+bge6zxeSBJGhBBCPCdHA8jR+R1Jkhw7QVCzMI4uihuNRnrCpzoZqFQqGIZBlmVMp1OGwxHnRilPRD67ZZN9biG96GUqMFJOOFM2vZBNZ4pfhPoKx3U9HT4ajQbValVfj6guGFX/cfG49KNL8lQ77tEuGFXXor7vix099ZDwcekkjAghhHhGR7tV1CTUOI71vI6jHSKqi2U2mzEajRiPx3oWiNqx0ul0jk0NPTw8ZPdwwKMTmyfSBrucJLposqlLxqYXcdKPOOFMaRCS5/OOFNdx8byGnt3RarWoVqu6vuPoLA61k+ZoG64aFQ8cCxDqakWFmDiOn/Zzo65rJHy8MCSMCCGE0I7uebl4fLpaIqd+XA0gU0vjhsOhvrYxDINKpUKz2dSDvaIo4vDwkCe6Qx6eeGzlLfZ5KYVh6gFjNgUbbsTN1YRTQcqSHZGlib7OcZyaDhS1Wo1qtUqr1aJSqRx7XkdrNNRpieqMUac3qqsHODZQ7GjHjyLXLpeXhBEhhLiBqRdmVf9wdMx4HMeEYfjUqPQng8nR6wpVK5KmqT4BaTabegrpYDDg4OCQRw8jHpq4nC/a9I31+Rd/8vW8aSbcXov5mk7J6VoJxfxrz0NMFctq6NoNdf0SBIE+7VDBQgUgFVbUqHXLsnSRrbpiORpGjo5OV46GGXVNIy4fCSNCCHEDuXjL7dF/T5JE13uoEALof1dXLqoLRp1CeJ5Huz1fb59lGcPhkO29fR48zHlkVuFcsUho+PMnYACUbHgxL62nvKJdslZVba/mk903Ho1GA5jP7KhUKjQaDV3DARyb7aG21qrTErXVVgUQNYlVfb+Avro5ejVz9NpFul1eXBJGhBDiOndx+Di6+0TVQ6h6D3UKoiagqgV06lQhiiI97tyyLN2Wu7u7yyDMeHBk8eiswnaxTmY8+RJjzK9fTlciXtZIeWkzp115aombOnk4uk1W7XJR9R/qY4+eWKgOF3X9olpsn2nGiRqY5nnesVoSmfNxdZAwIoQQ15mLr15Ul8vRZXFqDHmWZXrg2MV1IoZh6FoRNVvD8zzG4zHdbpckSehGBg9PfR6Z1eiWDTiyZK5qpNxcmXFbNeKlLQPfVcHD0SccqqNFBQ8VMlQNhzr5UNc0vu9TrVaPbbdV38fRUfNHJ7Z6nndNLo+7kUgYEUKI60BZlsc2vAL6hCCO42P7X9TodTVwLM9zfSWjPk8YhmRZpq9GhsMhOzs7FCXsZhUenvqciduMODL11IAFM+S0N+H2WswtbRff9zCMCkVR6BoP3/f1zA41Bv1oa+3FVy/q49UsEnW1kuf5005zgCcX2XnH9rpI0enVTcKIEEJc49TUU1XvcXQPDHBs/ocqPFXBRYUXdbqgHquCzGQyIcsLziUB908X2MoaxDzVEmtSsG5NOe2Oub2estGuPHnqYevZHmqnS7Va1dcjgA4J6mvbtq2LUyuVij79OFpcWhSFrmFRxbSALp5VM0RuxP0u1zIJI0IIcY1SJxp5njOZTBiPx/rF/eibOglRQUNd1agiziiKGA6HuiZEXeMcRAZfnDR5KGkxLZ+a/+EZOTfZI067Y05XYhab1SeLRudhw/d9KpWKXiqngoY6sVATTPM8x7IsXfehrmouDhLq+avOnaOnKOrqRn1euX65NkkYEUKIa4waxKXGsPd6PV1foa4q1L+rawxVaAroXSy9Xo/BYKD3wqRpSpgWPBLVuD9qspNV9df0jIyXuCNucYfcVEmp16pUq41jVypq4Jj6+upr1et13WmjTktUi646+VCzSI5+j6p1+OKJpyqAqJMTuX659kkYEUKIa8jRRXSHh4ckSaInms430RrMZjPCMCQMQx0yVG1GFEUcHBwwmUx0sWocJ+zmAQ+nyzwc1UhLdbpQcsIac7vb4xZ/RrtRo1Kp6aVyR0emq3kf6vRDnVQcXUzXbrfx/XmLrzqVORomjtayRFGkT0DgqRZftRtGXF/k/6gQQlwDVD1HFEX0ej3CMMQ0TV0LYpomg8FAh5U0TfWsDcMw9H6Y8Xisiz1npctjxRJ/Pq3RS596OWgYMbc5h9zmDlipz0OH5zV13YfqhFHXImorr9r3osLH0foNeOpaRX0MPHsNiPr8EkBuDPJ/VwghrmLqxVptvZ1MJnqWRpIkuK6ra0aybL6z5Wir69FtuWVZkpUG5/IO90dNHps6lE+OQXWMgtNWn1utA054Ma1WkyBYIggCWq3WsUJS1Y6rTkXUSYgKH2qOh3J04ZwKTyqAXLz3RV0hVSoVHXrE9U/CiBBCXKVmsxmTyYQwDJlMJiRJwmw2I0kSPbxrOBzqUehhGOqV9qpgdT4vxKRPjfvjJn8+rhDmTwWFNWvKrVaXm4xDOo0qtVqNen2Vdrutt95mWaZHrKvuGNXtotpzj16pwFMj1tWPqfoWFaIuPgFRNSBHp6KKG4eEESGEuIqUZXksfEwmE/r9vu4kqdfrVCoVptOpfsFXgUXtjFE1Ionhcn/e4bN9l93oqeLQmplyi3nAqXKXRb+k2WzSbJ6k1WpRr9epVqskSUJZlnrHi+u6NBoN6vW6rk0B9HAx9dzVJFVAT29VAUmd3KhrnKNFqOLGJr8ChBDiKqBmhKjC0slkwvb2tj5F8H2fRqOh54n4vk+/39fXNtPpdF48Wgk4m9T4dNfhgaFFXs5Dg0XBKWfELcY+KwwIKj6t1gILCwu0Wi2q1eqx56NqQ1zXpdlsHqsVURNO1ZULPLVYDjg2t0Ttr1GTVdWVjZyAiKMkjAghxBWkri3UycdoNGJ7e1uPYFcj0lW3iuu6jMdjJpOJbt3N85zSb/CJA5fPPGYzzp6atbFiR9xqdTlt9XDJqNVqtFonWVpaol6v4/s+RVHomhB1ZaL2vhztipnNZjqAqImqR4tRFRVGji62Ux8rg8jEM5EwIoQQLzLVbquCyHg8ZjgccnBwwHQ6PbZPRY0xNwyDJEno9/t6AFiapiR2lT849Ph0zyV78hQkMHNusXvcZh/SNkMMwyAIAtbX11laWsL3ff05VS3IxZNLVSHq0cFqgA5IqiVXBSW1DbgoCh1OVDuxnIKIr0TCiBBCvEjUOHa1kG4wGDAajRgMBkwmE7235ehmWvU4dYWjZnAMMotPDjr8l75L/mRHzJoT8kpnn5ucMUY5b/dttRbZ2NhgYWFBz/xQQ8NUEap6/1HPNGpd7XtRz+1oqDpawCqnIOJSSRgRQojLTE1MVXNCRqORbtNVNR+u6+qrjaNzObIs06Pep9Mp3cjgT4YNvjiuUDwZQk44M17j73PCmT25NNdkZWWNzc1NPM8D0KcfqmjU8zx9mqHmjqhBZGo2yTONWle1IKpzR5FTEPHVkDAihBCXkZqpobpi1HTU8XhMURTHrmHU9YmqsxiNRvR6PUajERfGOX82W+SBaUXPBrnJnXGnv8+yMQKgWm2wvr7O4uKiPr1QhaeO41CtVvE8jzRN9Zj4LMt04FAnMkfnfBydjnp0b40ipyDihSBhRAghLhO1gE5dxcxmMz37Q7XllmWp52uoK5rBYMDu7i79fp9z44JPh4s8HFbhyRBy2pvyGm+fZWtKnufUanU2NzfpdDp4nqfDBMzDQqPRwPd9JpMJ+/v7xHFMnudUq1WazSaWZelJp0fbdgF9DSOnIOJykjAihBAvsCzLGAwGOoiMRiM9bdSyLH0N43mertsAGAwG7Ozs0Ov12BqX/Fm4xGPRUy23t3hjXu3usWjNZ4oEwVMhRHXGmKapa08ajQa2bevPqxblVSoVOp0OlUpFd8wcDRVlWep5JRfvh1HdPUK8kCSMCCHEC0TVdxwtTJ1Op8dmcaigoIaHFUVBv99nb2+PXq/H46OSP5stcjZWIaTkJe6IO4NDFu356PRarcX6+jqdTkcvn1N1HipgGIZBv99nMBjoQtRKpcLCwgK1Wo1KpfK0U42jBbaK6riRUxBxOUkYEUKIr4IadR7HsT4JOTw8ZDweH5s4ahgG1ep83LoKIYeHh+zv79Pv93l4UPKp6SLnkgAAg5KXuENeGxywUpl/jmq1w8rKCktLS/rqJcsyiqLQc0EMw9CnMiqgVKtVlpeXaTQaTzvVUM8/TdNjI9qPDikT4nKTX2VCCPE8HG1rDcOQ4XBIt9ul1+vpIk/1wl+pVGg0GjqU9Pt9ut0u/f6AB3o5fzLpsJ3OQ4hJye3egNdV+6xUTRxnXsuxvLysTzXUSUiSJFiWRRAEWJbFeDxmPB7rrpwgCFheXqZarT7tVEM9/mhBqmEY+irm4lZfIS4nCSNCCHEJjl5lqCuW/f19Dg4OiONYv4irYWK1Wk0PBpvNZuzu7jIajflSr+QT/SZ7uQohBS/zBry+PmC96RMEHX2do5bWBcH8Y6Mo0lcyrusynU6ZzWYAepfMwsKCPik5KssyPVJeUbti1CAzIV5sEkaEEOI5uHjb7GQyYW9vj729PcIwpCxLyrI8NlBMzfgoioJz585x2OvzxcOSP+zXOHgyhFgUvMIfcHd7yno7oNm8SZ9ONJtNfN/XgWY2m1GWJUEQYNs2YRjS7/f1x6vC1ItPQqQgVVztJIwIIcSzeKYX8SiK2N/fZ2dnh+l0eqxDpdFoUJYllUoFy7Ioy1LXhfz5fszHunV65bww1SbnlcGQr1+MObXaoV5fBeanFOoqJggCHMdhOp1SFIWegKrmhKjTEd/3n/Ek5NkKUtUCPDkFEVcLCSNCCHGRZ3oRV3thtre3GQwGJEmig0ez2QTmQUJ1qUwmE3Z3d3l0p8/v7AWczedhwyHnjuqAN63D6fX5srokSQB0m616UyPg1ZWPZVnEcaw7XFzXZWFh4djGXTUh9egpDkhBqri6ya9KIYR4ktoBo+opiqIgyzLG4zH7+/t0u109Bt11XTqdDpZlkSQJrVYLwzBI05Rut8sTF/b42AWDL8RrFJgYlLzCH/BdtzqcXj9BEATEcUyaprr2Q72pRXhqGqraH1MUhR6Q1mg0qNVq+nRDjZyXglRxLZIwIoQQzGtCoigC0MWp0+lUF6iq8e22bdNut3FdlzRNdZFpnucMh0POnb/A7z0+5U9mS8S4AGw6U95yIuaO0yu02209hVXtfFF7X9TOF9UJo/bJZFmmP04NM1N1HuoqSZ2ugExIFdceCSNCiBueCgdq7HmSJAwGAw4PD+n3+2RZhmVZulNFFaq2Wi2yLCOOY86dO8cfPrTLx0eLDNgAoGXFfPtaxOtOVFldvUkHB3XSAU+12ML8KqVWqxEEgf4YtUTP8zw9KA2eOYSo5ygFqeJaI2FECHHDKsuSKIr0rJA4jplOp4zHY7rdrm7VVbUgjuPodl3DMMjznL29Pf74C4/w+70G5zgFgGfkfPNKxJtO2Cy0V/RwMt/3yfNcn8CoxXQqhDQaDd2Bo2pS1II7VYvyTEW1EkLEtU7CiBDihlSWJWEYEkWR3mA7HA71KPeyLMnzHMuyaDQaVKtVvd22LEt2dnb47Bfv5z/tujxYnKIw5nUhr+/EfMcpk3YQ0Gg0APT8jvF4DKALSW3bplar6WsfFTTyPNdXNUc7c9RJyNEQIkWp4nogv4KFEDccVQ8ShiGz2YzZbKanoqqx6KZp0mq1aDQaVCoVqtUqRVFw4cIFHnr4Ef7zVsxn0g1iwwUDXlLP+Ks3w2rF0jM/1ImFmkNi2zaO41CpVKjVarRaLX2aoaa5WpalB5w5jvOMIcQ0TV1vIsT1QH4lCyFuKGqZ3XQ6ZTKZ6BbcyWRClmUYhkGr1dIts7VaDYCtrS0eeeQRPnNuzB+HawyMGhiw5BX8ldMlt9VzbNsiCKr6SqcsS2azmS5Q9TxPb9iFebdLURS6XbdSqehhaYZhSAgRNwz5FS2EuGGoZXaz2UzXhRwcHJDnOXmeEwQB6+vrNJtNarUaaZqytbXF448/zqO7A/7zcIFz5UvAAN8s+I7TJm9cKinyFNd19VWO4zjkea6npZqmSaPRoF6v6yBhmiZxHFOWJb7v47ouvu9jWZZevHc0hKjNuUJcjySMCCFuCJPJhNFoxHA41O26URTpyaadToeFhQVarRZ5nnP27FnOnz/P2fO7/NGwwf3ZLRSYmJS8cbXkO2+2sPOYsij1iYaqEYnjGN/3MQxDt/6q0xHLsnQHjmrdVW+qvVgNKzMMQ09dFeJ6JmFECHFdK4qC4XCot+oeHBwwGAwAdHHq8vIy9Xod0zQ5e/Ys29vb7Ozu8Zm+y6fj00TMw8DtzYLvvsVkvWowm01xnrxSUe24YRjqgWQAtVpNF6+qmR9RFOE4jv4xz/PI81yPfAeOTVgV4kYgYUQIcd1KkoTDw0MODw/Z2dlhNBqRJAmWZVGpVHRdSKPRIAxDvvSlL7Gzs8MDhxl/HG3QKyoALPklf+VUwR0rzpNTU9H7Yzqdjm7XVScganKq2ltj2zZJkuA4DkEQ6MJWgDAMnxZCZFiZuNFIGBFCXHfKsmQ6nbK3t8fu7i4HBwdEUYRpmrootdFo6OFivV6PL37xizyy3eOPZiuczVoAVKySb9ssedMJE9ty9MmH2kfjOPP3BUGgT1Z83ydJEmazGZ7n6Tklqh5ELaiL41hCiBBPkjAihLiuZFlGr9fjwoULerOu2ulSq9X0FNNqtUqapjz++OPcf/8DfKoLf5reSoo1rwtZKfj2UyatikOSJIRJrIePNZtNPXSs3W7rDhc1uVVt3lUL744GjaO7b2SDrhBzEkaEENeFoigYj8f0+322trbodrt6p4vqjlGnGsB8j8y5c3zu/of5/ckqW+UiAKdrOX/9JSYbdVu316p2X3UNk2UZtVpN15yEYajHxtdqtWMj3Y+GEHVtIyFEiOOe1xrH97///Zw6dQrf97nrrrv41Kc+9WU//r3vfS+33347lUqFzc1N3vGOd+hxyEII8dVQszy63S5nzpzh85//PPv7+8C8rmN9fV2fXqhrlV6vx5e+9CV++zOP8X+PX8JWuYhJyVtOlrzjVSbrNUPPAAmCgE6no0fAqx0x9XqdSqXC7u4uvV6PoihYWFhgZWWFpaUlms3mseJUFUTUeHc1S0QI8TxORj70oQ9x77338oEPfIC77rqL9773vbz5zW/moYceYnl5+Wkf/yu/8iu8853v5Bd/8Rd5wxvewMMPP8zf/Jt/E8Mw+Lmf+7kX5JsQQtx41GTSyWTCeDzmzJkz9Pt9vbG21WpRrVb1lUiWZUynU6Io4v6HHuZjOx73Fy8FDJb8kr/1couNSobneRRFQVmW1Go1fZKihpHBPOQMh0PCMATQ80lqtRqmOf87XlEURFGkv/7RuhEhxHFGqabqPEd33XUXr3vd63jf+94HzH/DbW5u8iM/8iO8853vfNrH//AP/zAPPPAA9913n37f3/t7f48//dM/5ROf+MRz+pqj0Yhms8lwONQtc0KIG1NZlnqx3XQ6pdvtcv78eT091XVdVlZWSNOULMv0x6dpynQ65TOPXOA/9pbol1UA3rBS8ldvhprv6H00ruvSaDQoyxLP8/SCO9UJMxgM9P6YtbU1FhcXj51yqC3AIG264sb2XF+/L+lkJEkSPvOZz/Cud71Lv880Te655x4++clPPuNj3vCGN/Cv//W/5lOf+hSvf/3refzxx/nIRz7C937v9z7r14njmDiOj30zQgiRZRlRFBGGIaPRiJ2dHQ4PD7EsC9M0abfbBEHAeDymLEuKotC7ZkajMf/xoSF/NL2JHJPAKvie201etTAPC1mW6X00juNQFAXVapWyLHWh6nQ6ZTQa6TqUtbU1fVoC6BZf1SWj5oio0xIhxDO7pDCixiavrKwce//KygoPPvjgMz7mr//1v87BwQFf//VfT1mWZFnGD/zAD/AP/+E/fNav8573vIef+qmfupSnJoS4zsVxTBRFDAYDBoMB+/v7enaH67q0223SNKXX6wHoq5YoitjujfnVxyy2sjUAbqtnvP3lNgtVB9M0yfP82DWMmspaFAWNRgPbtvUm33q9Trvdpt1u69OQsiyfVqDq+77skBHiObrscf3jH/84P/MzP8M/+2f/jM9+9rP8+q//Or/1W7/FT//0Tz/rY971rnfpiYmq4l0IcWMqioLZbKbnhqjZIVmW6WLQdrvNYDCg3+/rlts0TRmPx/zZhZD//8N1trIGllHynSdi/u5rXDrBUztiVJdMHMd6D4zv+ywsLOj3q8mq6+vrdDodHUTU9Y8KIqptWIKIEM/dJf1uWVxcxLIs9vb2jr1/b2+P1dXVZ3zMT/zET/C93/u9fP/3fz8AX/u1X8t0OuXv/J2/wz/6R//oGY8vVW++EOLGphbGzWYz+v2+HuXuuq6eYlqWJTs7O1iWhW3bpGnKZDIhzko+/Bh8ejDfkLvsJnzPLRm3rzb0yYXarKvad9WI9k6nQxAERFFEmqZ6fHuj0dB7YqRAVYgXziWdjLiuy5133nmsGLUoCu677z7uvvvuZ3zMbDZ7WuBQv1kvsXZWCHGDUNcrarvuwcEB29vbDIdDLMvS1yhq8V2lUsEwDKbTKXmeszUu+adfsPj0YN4J89r6iL/3ypLbVupYlqWnpwJ6YZ3qnLnpppsIgoDZbEaWZQRBoK9ljj5GfS3V7lutViWICPE8XfI54r333sv3fd/38drXvpbXv/71vPe972U6nfL2t78dgLe97W1sbGzwnve8B4C3vOUt/NzP/RyvfvWrueuuu3j00Uf5iZ/4Cd7ylrfIb1whxNMURUEYhiRJwnQ65eDggIODA9I0Jc9zLMvCMAz6/T6+7+s5I2VZYlo2HzmT8LEdl7w0qJoZ377Y4w03t6lUKnqfjCoyTZKEen0eUFqtFgsLC0RRRJZlen+NmrpqGIYUqApxmVxyGHnrW99Kt9vl3e9+N7u7u9xxxx189KMf1UWtW1tbx35j/viP/ziGYfDjP/7jXLhwgaWlJd7ylrfwv/wv/8sL910IIa4LaZoSRRFRFDGZTOj3+3S7XfI8pygKLMvS9RtqB0xRFBiGwcEs5988ZvDYeH7F+xJ/ylvWxpxaW6JardJsNnWhaZZlALo4dWNjA8dxiKKIsiz1ortqtYrruvqkRgpUhbg8LnnOyJUgc0aEuL4dfbEfjUaEYchwOOTg4ECPZLcsizRNsSwLx3H0HA/TNPlMF371EQhzcIyCb24dcmc75cSJDT0tVdV3qGsZdTWzurpKlmUURaFPXiqVCpVKRdegxHGsr5Udx5HpqUI8R5dlzogQQrzQ1NWHCiBlWdLr9RiNRhiGQRRF+kTEtm0sy2I2m81PRrD41w9m/OnePCisORHfvnDAiZbH+vpJms0mlUqF2WxGnuckSYLv+3iex+LiIq1WS5+SZFmG67o6iAD6cSAFqkJcThJGhBBXTJIkeoCZaqvd29sjDEM9/NBxHNI0pVqtEkURSZLQaDR4dFjy81+IOYzAoOTrGkPe0Byy2GmztLREu93GNE2m0ylZlulZIkEQsLKyQqVSoSgKPRjN930qlQq+7+vTEEVt3hVCXB4SRoQQL7qyLAnDUO+VUVce+/v7hGHIYDDQHTOGYdBsNpnNZvNJq7bDvztT8B8eTSiBlpPzX3W63FwvaTQWWVpaYmlpiSRJ9NVPWZZ64Z2aHVIUBXEc62V4atT7bDaTAlUhXmQSRoQQL6osy3RdSBiGOI5DlmX0ej36/T6j0UjPEHFdlyAI9Aj27iznFz+fcWY0v5a5oxnzLe1Dqq5FEFTZ3NykVqvpExQVRFSnTLPZ1LUnYRhimqZu6VVdOSAFqkK82OR3mhDiRaN2uyRJQpZlVKtVwjCk2+2ys7Oja0HUBFPP85hOp7iuywMHKR+8v2SSQcUqecvqhJfXI/LcpFqtcuutt+qOmDRNddGrurKp1+sYhkEYhkRRhOM4BEFApVIhyzJ9LWPbNr7vS4GqEC8iCSNCiMsuyzIGg4G+OlEj2MfjMefOnWN3d1dvxa3VavpkZDqdYlkW/+lsxL993CAvDU4EJd9zaopfzE8/lpaWuOWWW3Tbrnqr1+ssLy/rbposy5hMJqRpiud5BEGA7/vHdsq4rivTn4W4AiSMCCEum6IomE6nTCYT3c0SBAGu67Kzs8P58+d1+26r1SIIAj3RdDKZUAL/6otT/nB/Pvn0VZ2C/89NEbPxiLgoOH36NKurq/pzq4Fky8vLtFot/Tlns5muBVGDzBzHIQxD3S2jxsMLIV58EkaEEC84te9lPB7rK5myLOl0OkRRxBNPPMG5c+cYDAb4vk+r1cJxHBxnvkV3NpsxTXLe/19CHpnMu1jevJHzpuWIw4Muvu9z22230W639YnIbDbDdV3W19ep1+t0Oh1M02Q8HhOGIYZhUK1WCYJAfw01MK1SqUjLrhBXkIQRIcQLShWHqlOHPM/xPA/f9xkMBuzt7XHmzBmiKKLVatFoNPSSO5hf6Wz1Qt7/xZyDxMUxSr7ntpLbKzMG/QGtVouTJ09Sq9Uoy5LxeEyaptTrdZaWlqjX67TbbT1ALYoibNvWHTPHxsebJpVKRbplhLjCJIwIIV4QaoqqegP0Jtw0Tdnf32d7e5tz585hWZae9WFZlt43k+c5n3hkn3/1mEdY2DTsgh96lU0QHTCZxCwuLrK6ukqz2SRNU4bDIY7jsLCwQLvdPrZtV9WHuK5LtVrF932yLNPPTU1alUJVIa48CSNCiK/a0eV2YRjqUw7HcfQskSeeeIJer4fv+zSbTX1KYVkWrusyGo34jS92+Y/bAQUGJ4KcH3qlTTHtgWnqGSGNRoPZbEYYhlSrVer1OgsLC7RaLWzbZjqd6o26vu9Tq9VwXZc4jvUIebXbRghxdZAwIoT4quR5rq9k1Mh2ddrQ7/cZj8c89thjZFlGrVaj2WzieR5FUeC6LrZts7u3zy99Ycyn+lUAXrOQ8z23mWTxBMtxsG2bVqtFvV5nOBzqbpxWq8Xi4iK1Wg2A8XisZ4VUq1Wq1SqWZRGGoR77LtNUhbj6SBgRQjxvasuuKiI1TVMvnQvDkIODA3Z3d4F5t8rKygp5nlOWpT6ZeOSJC/zzLyScDQMAvuNkyV/cMCiKDMdxMAyDRqNBvV5nNBrpTpjFxUXa7TZBEJBlGePxWO+xUYWqaq6IdMwIcXWTMCKEeF6OXntkWXbshX86ndLtdhmPx5imieM4LC8v62Fjvu8ThiGffmiLf/GIQz/zcYyS/+7lJl/bLgEDMMnznEajodtzG40GjuOwurpKo9HA8zxdH5IkCY7jUKvV8H2fsiyZTqeUZSkdM0Jc5SSMCCEuidoro04bDMOgKApmsxlxHDMejxkMBkRRpNtpW60W0+mUarWKaZr0+33+05e2+dD5GnFp0XQKfviVNidqAPNx7UVR0Gw2cRyHsix1+FhbW6NWq+n6EDXDxPM8arUanufpQlXpmBHi2iBhRAjxnKnrFzWfw3EcxuMx0+mUMAwZDoeMRiPSNMU0TV0fEkURzWaTMAw5PDzkN7/U5/d6LUoMTtUKfvgOl6pVYJomURSR5zntdvvYiYZq3a3VahiGoetD1BI8FVCSJNGj3aVjRohrg4QRIcRzcvFpg+d5DAYDJpMJ0+mU8XjMeDwmjmNc16Ver2NZli5cHY/H7Ozu86/uj/jCrA3AXSvwvS91sJgHEbVFt9ls6loRwzBot9ssLi4SBAFFUehFe2rRnTpxUXtpQDpmhLiWSBgRQnxFz3Ta0O/3ddvucDjUcz1836der2PbNmVZ4nke3W6Xx8/v8S8fdzmfNDEo+e6XOHzTaobJvL13NpvpcOG6Lo1GA4CFhQU6nQ6VSuXYVFfbtvXGXUA6ZoS4hkkYEUI8KzXITL3Iq9OGwWDAaDRiOBwyGAwIw1AHETVRFeanKfv7+3zxiS6/eqHJMHfxzJIfuMPn1kqIadqYpqlnk7iuS61W0yPbFxcX6XQ6uK6r60NUS3C9Xtctwkevjnzf13NOhBDXBvkdK4R4Rkdf5OGptlhVoNrr9RgMBqRpSpIk1Go1vXFXjYQ/PDzkT7bG/PvuIklpslgx+B/u8GgZMyxr/sdPmqbYto1t27pI9WjrrmEY9Pt9PTm1Uqnokxc140Q6ZoS4tkkYEUI8zbO9yE8mEw4ODnQYUYPOms2mHjw2mUyIooj9/S6/vw2/31+ixOD2jsXfebmBmY6xnrxCUZ/fdV1arZaeEaImqqZpymAwoCgKPehMzQ/JsowwDAEwTVO/Xwhx7ZEwIoQ4Rg0yA461xc5mMw4ODuj1ehwcHFAUBXme02w2qdfrlGXJcDhkPB7TPezxH3aqfH48n6j6zSd9/usTEWWWYpgmlmUdO3FRxa7tdlsvzzs6TdVxHL3ZF47XsNi2je/7EkSEuIZJGBFCaEcHmR19kQ/DkG63S7fb1UEEoNls0mg0yPOc4XDI4eEh57sDfqO7yFboYgDf98oar64OyZKUsiyp1+vAPOi4rksQBLiuS6fTodVq4fs+h4eHuk6lVqtRr9d12DjaMeO6Lp7nvcg/S0KIF5qEESHE0waZHe1GUUFke3ubXq8HoNtu2+02SZLoIPLY7oBf219mkNoEtsGPvLbKCWvEdDrT+2VgXtjq+z6e5xEEAe12m3Z73u57eHioF+i1Wi0dNi5+jjLaXYjrh4QRIW5wz9aNol78e70e58+f5/DwEEDP/eh0Onriarfb5cz+kF/bW2GQWaxULe59bQUv7jOdRniep4NImqYEQYBlWTrQ1Ot1kiTR10O+79NqtfTUVOmYEeL6Jr+bhbiBPdvYdBVEhsMhZ86cYTgc6tOMVqvFwsICaZoyHA7Z3d1l62DEr+0uM8gs1moW73y9TzmbzyFR4+AtyyKOY331o6apBkGgZ4SoSatHp6bKaHchrn8SRoS4QT3b2HQVRAaDAQ8//LDujqnX69TrdRYXF8nznF6vx+7uLhd6E/7N9iL9zGalavI/vtYlHR0QhiH1ep1Op6Nniaggsr6+rus9VNDwPI96vX5sWJmMdhfixiBhRIgbzLMNMlM/puaDPP7443rnTLPZpFqtsrKyogeZ7ezssN2b8K8uLNDPHBYrJj/2+grZcJ8oimi1WiwtLek2Ydd1sSyL1dVVbNvWE1Uty8L3farV6rEZITLaXYgbh4QRIW4gX64IVP1Yt9vliSee0K217XabIAhYW1ujLEv29vbY2dlhbzDlX1/o0MscFnyDH39DjaS/w2QyYWlpSQcXNbrddV0WFxePBRHP8/B9/9iJx5crphVCXJ8kjAhxA1EbcS+eVlqWJbPZjG63y7lz55hMJs8YRHZ2dtje3manN+bfbC9wmLm0PYOf/MY28eF5ptMpi4uLrK2t6TbhsiypVCp6v4wKFtVqFc/zjrXmSqGqEDcm+V0uxA1CFYk+WxDZ3t5mZ2eH6XTKdDql1WpRrVbZ2NigLEvOnTs3PxHpj/nVnQUOMo+mZ/DT37JE3N1iMpnQbDZZW1vTJxtqY2+r1dIj3B3H0TNMjrbmSqGqEDcuCSNC3ACO1oj4vv+0IHL+/Hn29/eZzWZPCyIAFy5cYG9vj+1un1/dWWQ/82i48I+/eYn08DyTyYR6vc7a2hrT6RSYD1Cr1+u0220ajQau6+q6kYuDhhSqCnFjkzAixHUuSRJdCFqpVPS1R1mWjMdjtra2GAwGTKdTZrOZ3v9y4sQJAM6fPz9v393Z59/ur7CX+dQd+JlvXSXvnWc8HlOv11leXiaOY72vpt1us7CwoFt1Lct6xtHtUqgqhJAwIsR1LE1TfeJwtP5CBZEzZ84wHo8Jw5AwDAmCgHq9zokTJyjLku3t7XkQ2d7jw90VdlOfqgM//S3LGKNdxuMxtVqNhYUF4jimKAriOGZpaYlWq0Wn09E1IRePbpdCVSGEImFEiOuUqsGAeRBQ9RlFUTAajTh79izT6ZQ4jpnNZlQqFYIg4KabbiLLMnZ3d9nb2+Ps+R1+/WCF7cQnsOGnvmmRIOlzMBxSrVZpNBqk6XzvTJIkLC8v0+l06HQ6OI7zjIWoUqgqhDhKfvcLcR3KsowwDIH51cfR/S6DwYAzZ87o65EoivB9nyAIOH36NFmW0e12uXDhAue29/jNwxXOJxUqFrz7Gzu0mXDQ7+vwkmUZtm2TpinLy8tUq1UWFhawbfsZC1GlUFUIcTEJI0JcZ/I81yciqkZDUVczqsV3Npvpjzl9+jRpmupdNOd39viNgyW24gqeBf/wjU3WvIRut3/sJMOyLLIsY2lpCd/39VCzZ6oPkUJVIcQzkTAixHVEXX+UZamHiymz2YzHHnuM2WxGlmV6F4zv+5w6dYo4jhkOh5w9e5YLu3v8RneRJ+IAz4J3vaHBTdWcfr+vr15M09RdOQsLC3iex8rKij6Jubj+QwpVhRDPRsKIENcJ1aZblqUOGUqapjzyyCNMJhPCMMQwDIqioFKpsLm5SZqmjMdjzp07x+5+l9/cX+BsXMUx4Z1vaHJrEwaDMaZpkmUZnufpeSG1Wk1PV/V9X7fwHn1eUqgqhPhyJIwIcR04GkRUHYa6/sjznEcffZTRaKTrSIqioFqtsrq6SpqmzGYztra22Nnb59d32zwe13DMkn/wdQ1ubZSMxxNgHmpUMazneQRBgG3bdDodHUpkoqoQ4lLJnwpCXOPUyYN6wT8aRIqi4MyZMxweHupaDZiPYldL7NT01e7BIf/Pdp3H4hqWUXLv62q8tG2Qpqn+OBU2qtWqDhbNZpNms3msUBakUFUI8dxJGBHiGnd030wQBPoFX41w39vb06FA1Xp0Oh2yLCOOY/b399nb7/J/n6vwSFTXQeRrFkwMwyCKImazmS44bTabeqz70c6Zo9dCUqgqhLgUEkaEuIZdvG9GBZGiKNjf32dnZ0d31qhi03a7TZqmpGlKv99nZ3ePD2/5PBTVMY2S/+HOKl+7aOJ5Hv1+n/F4XitSrVZpt9tYlqULUJeWlp5WKCuFqkKISyVhRIhrVBzHz7hvJs9zBoMB586dYzqdkmUZvu9TFAW1Wo2iKHSdyPkL2/z6OY/7wzomJT/86oDXrNhUq1W63S7D4ZCyLPW1jmVZeqz7ysrK005EjgYRKVQVQjxXEkaEuAYlSUKSJMDxMe95njOZTNja2mI8HpNlGdVqlaIoCIKAoih0W+/29g6/ed7ji9M6BiU/cEeF16051Go1Dg4OGAwGlGVJvV5nY2MDwzD029raGo7jHLt+ORpEju7AEUKIr0T+tBDiGnN034xqsYV5wajawNvr9UiShGazSZ7neJ6nr3PSNGVnZ5cPnzH4/GQeRP77V1V4wwmPRqPB/v4+BwcHANRqNdbX149t+V1dXZUgIoR4QUlpuxDXkIv3zahrkDRNCcOQ3d1d9vf3SZKEVqsFzAPLeDwmz3OSJGF/v8uHHy/57HgeRL7/lT7feHJemNrr9eh2u5RlqTf3Ht1ps7KyosfASxARQrxQ5E8NIa4ReZ4/474Z1blyeHjI9vY2URRRr9exbRvDMOh2u/i+T5ZljEYjfv3RjE+P6gD8za/x+JabawRBMC9m3dkB5qHi9OnT2LZNURQURcHi4iJBEBw7EVEFtOoxEkSEEM+HnIwIcQ04GkSOFo3GcUwcx3p6ahRFuv3Wsiy63S6O4+gTld9+POKTw3kQedvLHf7SbQ2q1Sqj0YgLFy4A86Bz+vRpHMehLEuKoqDdbtNoNI61DksQEUK8UCSMCHGVO7pvRs3sgPn1SJIkhGHIuXPnmM1mOI7D2toaWZaxv7+vH18UBf/5sSG/s18D4K/eavOdL+9QqVTo9/ucP39e77O5+eab9elHnufU63Xa7fax1mEJIkKIF5KEESGuYmq66tEgot6XpilJkrC1tcVoNMIwDE6ePEkcx3S7XV1bYhgGnzs/5P/ZrgEGf2Hd4K/dsYDruoxGI/b29ijLEsdxuOWWW6jX5ycnSZIQBAGLi4tUKhVdxCpBRAjxQpMwIsRVSu2bKYpCj1OHp8LAvD13m36/T57nnD59mizL2N3dZTQa4TgOruvy2P6EX3rMJytNXtEp+TuvbeM4DsPhkL29PbIsw3EcTp06Rb1epygKkiTB8zyWlpYkiAghLjsJI0JchZ5p3wzAbDYjz3N9DbO/v0+e52xubmJZFhcuXKDX62Ga8wmqe4MZ/+eDFmFhsVkteMddLWrVgOFwyMHBAWma4jgO6+vrtNttYN6ZY5omS0tLVKtVHTgkiAghLhcJI0JchS7eNwMwnU719NTBYMDOzg5ZlrG6ukqz2WRra0u35TYaDUazmPf/eUE/tei4Of/j62u0ahVGoxG9Xo84jnWNyfLyMmVZ6vklKysruiMHJIgIIS4vCSNCXGWiKDq2b0Zd15Rl+dQY9/PnSdOUdrvN0tISZ86c4eDggCRJWFhYIEoSPvilnAuhTWAVvONOn+VmhSRJ6PV6pGmK67osLS2xvLyMaZokSUJZliwtLektvCBBRAhx+UkYEeIqkmWZHiDm+/6xAlZVsHr+/HniOKZarbK+vs758+fpdruMRiMWFxfJ85xffSjngaGFbZT8wNcYnOrMQ8RgMNA1KJ1Oh9XVVSzLIo5j8jyn0+nQbrcliAghXlQSRoS4SpRleWy6KqDrRtI0Jc9ztre3CcMQ3/dZXV1lf3+fbrdLv99nYWEB27b56BMFf7RnAiV/4+aEV6wG1Go1RqORDjStVouVlRVc19VD09rtNgsLC8e+tgQRIcSLQcKIEFeJKIooyxLTNDEM49iJiGEYbG9vM5lMsCyLpaUlJpMJ+/v77Ozs0G63cV2XP9nJ+Hdn55/vLRsxbzxZm1/bRBFxHDOdTqnX6ywvLxMEAXEcE4ahDiJqmJoEESHEi0nCiBBXgSRJ9Iu/ujZRJyKWZbG3t8d0OsUwDFqtFlmW0e122d7eptFo4Hkej40tfvnBAoA3LIR860mH1dVV8jxnNBrR7/fxfZ8TJ07QbDZJkoTZbEar1dKj3kGCiBDixSdhRIgrrCgK3cXiOA5pmlIUBVmWYdu2rgfJskzP/Oj1emxtbeG6LtVqlV7u877PJ+SlwcvrEd99q83p06cxTZPBYMD+/j6WZXHy5Ek6nQ5pmjKZTGg0GiwsLFCtVgEJIkKIK0PCiBBXmNo5Y1kWWZaR57me/9Hv95lOp0RRpHfSjMdjLly4gG3bVKtVUqfKz306JMoNNisJ3/uSkpfceguO4zAajdje3sY0TU6cOMH6+rpemFerza9wGo2Gfh4qiARBIEFECPGikTAixBWkrmMMw6AsS/I8J45jPM9jOBwyHo8Zj8cAVKtVoihie3ubJEnwfR+/3uKf/umUfmLQcVL+1ktSXnbbrXiex3Q65ezZs5RlydraGrfeeitxHDMYDKhUKnQ6HZrNJvD0IKImrgohxItBwogQV0iWZSRJov9bXdf4vs9kMmE2mzEcDvX1TJ7nunbEcRzqzRb/v09P2Z4ZVM2c77815mtvv4VarUYYhjz66KOkacri4iKveMUr9LA0y7J0C68qlJUgIoS4kiSMCHEFHG3jPfrfjuMwnU6ZzWb0ej2SJMF1XVzX5eDggNFohGmatFotfvGLEQ8ODByj4G/eEnLnS0/RarV0EJnNZjSbTe644w7yPKfX6wGwtLTEwsICpmlKEBFCXBUkjAhxBag2XpgHkSRJMAyDNE2J45her8dkMsG2bTzPo9/vMx6PyfOcRqPBv3ss41NdE4OSv34y5I0v22RhYYE4jjl79iyDwUAHEcdxODw8pCgKlpaWWFxcxLIsCSJCiKvG8woj73//+zl16hS+73PXXXfxqU996st+/GAw4Id+6IdYW1vD8zxuu+02PvKRjzyvJyzEtS5NU7Is02GkKAq9h2Y8HtPv9+n3+7pgNUkSXTsSBAF/uJ3z0fPz37rfdSLiW1+xrjtkzp07R7fbpdFocPvtt9Nqtdjb2yNJEtrtNsvLy9i2rYOI2n0jQUQIcSVdcrn8hz70Ie69914+8IEPcNddd/He976XN7/5zTz00EMsLy8/7eOTJOFbv/VbWV5e5sMf/jAbGxs88cQTtFqtF+L5C3FNOdrGW5alXk7nOA6DwYAwDOl2u9i2jeM4mKapA0oQBDwwMPjwmflv229eSfiur12i1WphGAZbW1vs7u4SBAG33HILN910ExcuXCCOYzqdDmtra08LIqpVWAghriSjVH89e47uuusuXve61/G+970PmP/hurm5yY/8yI/wzne+82kf/4EPfIB/+k//KQ8++KDed3GpRqMRzWaT4XCo2xCFuBbNZjPyPCfPcyzLYjab4TiOrg/Z2dnRU1gXFxfpdrvs7OxgGAb7mcc/f8gnKQxe00754dc19NTUra0t9vb2sCyL06dP89KXvpT9/X2m0ym1Wo3NzU19yqLCkJyICCEut+f6+n1J1zRJkvCZz3yGe+6556lPYJrcc889fPKTn3zGx/z7f//vufvuu/mhH/ohVlZWeMUrXsHP/MzPkOf5s36dOI4ZjUbH3oS41qlldGpRnZodMhwOdadLURSUZUmr1WI4HNLtdknTlFHu8IuPeCSFwa21jP/v65p0Oh1c1+XcuXMcHh7qWSI333wzvV6P2WxGpVJhfX0d3/d12zDMl/BJEBFCXC0uKYwcHByQ5zkrKyvH3r+yssLu7u4zPubxxx/nwx/+MHme85GPfISf+Imf4Gd/9mf5x//4Hz/r13nPe95Ds9nUb5ubm5fyNIW46uR5TpIkOmyozbyTyYQ4jplMJnoXTa1WI01T9vb25icplsf/9XiFcWay6hf86GurLHba2LbN7u4u/X6fLMs4ceIEJ0+eJI5jxuMxnuextrZGEAR6+y/Mp7w+31NKIYS4HC57N01RFCwvL/PzP//z3Hnnnbz1rW/lH/2jf8QHPvCBZ33Mu971LobDoX47d+7c5X6aQlw2R9t41YlglmXHQsh0OgXQbbznzp2bnwiaNr+yVWM/Mmk4Bfe+1mdzdRHTNOl2uxweHpJlGRsbG6yururx747jsLy8TL1eP7Z0zzRNvQxPCCGuFpdUwKpaAvf29o69f29vj9XV1Wd8zNraGo7jHDsSftnLXsbu7q6eoXAxz/PwPO9SnpoQVy01ZTVJEl1AWpalvpZRdSQAjUaDxx57jMlkgmFa/Lu9Jo9PLDyz5N47PW47sURZlhwcHHB4eAjMTyZXVlbwfZ/hcKiHmjWbTUzT1NdDqmBVCCGuNpd0MuK6LnfeeSf33Xeffl9RFNx3333cfffdz/iYN77xjTz66KMURaHf9/DDD7O2tvaMQUSI60mWZaRpqrfvqhOSg4MDfWKS5zlZltFsNtne3iYMQ4qi4A+HbT7Xt7GMkh++w+WOU8sURcFgMKDX62FZFq1Wi6WlJXzfJwxDDMOg2WzSbs+vcY5OefV9H9OU0UJCiKvPJf/JdO+99/LBD36QX/7lX+aBBx7gB3/wB5lOp7z97W8H4G1vexvvete79Mf/4A/+IL1ejx/90R/l4Ycf5rd+67f4mZ/5GX7oh37ohfsuhLgKFUVBFEUURUFRFHob78HBgS5gVScjzWaTwWDAYDAgiiL+PF7gP+3MTxPf/jUuX/+SJWAeYnq93nwcfL3OwsICnudRliVFUVCr1Wi323iep78+zP8iIYvvhBBXq0v+0+mtb30r3W6Xd7/73ezu7nLHHXfw0Y9+VBe1bm1tHfvb1+bmJr/zO7/DO97xDl75yleysbHBj/7oj/JjP/ZjL9x3IcRVSIWNJEn0qchwOMQwDD36/fDwkCAIiKKIfr9PHMc8kTX5za35b83vutniO16xjGVZbG9vMxqNcF2XZrNJvV7H8zy97bdSqdBqtQiCQH/9siyxLEuuPYUQV7VLnjNyJcicEXGtUfM8oijSo9cnkwlJkpDnOY7jsL29jW3buK7L3t4ew+GQvcjgF88tkhQGX79m8D9+0waWZXHu3DmGwyG+7+vTkCAIqFQqFEWBbdv6esYwDOI41iPmgyCQ6xkhxBVxWeaMCCG+MjXPI0kSvYxOBZEsy3Ach263i2maWJbFaDRiPB4zTQo+vNOZzxJplPy9bzqBbducP3+efr+P53m6PqRSqWDbNoZhYJom1WpVT2KVOhEhxLVG/pQS4gWkilKzLCPPc9I0ZTqdEsexvkoZjUbEcayDw3z6asq/77Y5SCxabsGPf/MaFc/lwoULHBwcUKlUdBBxXVd3xhRFQaVSodPpYJqm1IkIIa5JEkaEeAGp0BHHMWVZMhwOdWtvpVIhiiK9jddxHA4ODojjmD/oVXl44mEbJe/8+kXWO3UuXLhAt9vVQUSNb7dtm3q9TpIk+L5Pp9PRoUO1DUudiBDiWiJhRIgXiGrjVXUi/X5fb+j1fZ+yLOn1epRlied59Pt9oijiC4fwh/35Xer331HjVZstzp07p2tKlpeXaTQaGIaBZVlUq1XiOMb3fZrNph5ipgpmZZ6IEOJaI2FEiBeAup5RVyRqMqoKDbZt6z0zlUqF6XTKeDzmzMGM3+7Pt13fc9LhO75mkb29PXZ2djBNk/X1dRqNhi56VaPifd+nWq1Sr9cB9CwTgEqlgmEYV+YnQgghngcJI0K8AKIo0gWqw+FQB5FqtYrruhwcHBCGIY7j6OubC/uH/IfBGnFh8JKWwQ/ctcxkMmFnZ4csy1hbW6NWq2FZFq7r4vs+aZriOA6+79NqtQCO1YmoVl8hhLiWSBgR4quUJAlJkhBFEdPpVP97pVLBcRwGgwHD4VDvhTk4OKDbPeBjo1UOEpumC+/6xmUso+TcuXPMZjPW1tZotVp4nodhGDiOQ1EUOpSogtWjC/BUm7AQQlxrJIwI8VVQpxJq2V0YhvpExPM8kiTRdSKtVoter0e32+UTvYBHwyqWUfLOb1igE9hsbW0xHo9ZWFig1WpRrVYxDEN3z6gTEjXqHY7XicgCPCHEtUrCiBBfhTAM9YlIGIYkSaJPRAzDYHt7m6Io9NCf/f19vtDN+ZPpIgB/+zUNXrbks7OzQ6/X050znU4HwzCwbRs1l9D3fRqNhu6SUcWxIHUiQohrm4QRIZ4nNWFVDS1TRaaqYPXcuXPkeY7v+zqYPLY/4vdmm5QYfMspj2+7rcnBwYHehL2xsaG3Y1uWpXfOBEFAEATUajVgPlhN6kSEENcLCSNCPA95njObzXQQUacX9Xody7LY2dnRm3rr9Tpnz55lp9vjo5ObiAuTW1smP/R1y/R6PXZ2dkiShLW1NZaXl/XiuyzLyLKMWq2G67q6YFV17oDUiQghrg8SRoS4RGVZMp1OGY1GDAYDfXqh9sIcHBwwnU4py5Ll5WXOnDnD/n6X3xkscZh5ND2DH3/TKrPJiN3dXSaTCWtra5w4cYIgCHQQSdOUarWK4zj62gaeqhNRBbFCCHGtk1nRQlyiJEkYjUb0ej2KogCg0+kA6IBimiadTocLFy6wt7fHJw58Hk9bWAb8w29cxs1Ddvb2GAwGdDod1tbWaDabwLwWRNWeqBMRdQ2j2ocBff0jhBDXOjkZEeISFEVBr9fj8PCQNE0py5Jms4lpmkRRRK/XwzRN6vU64/GYs2fP8oX9lD9L1gD423e2uKUJOzs7HB4eUqvVOHnyJCsrK9i2TRzHxHGM67pUq9VjBatqAR/Mg4jUiQghrhcSRoS4BOPxmF6vx3g81qHD8zziOKbb7VKWpd6o++CDD7J1OOX349OUGNxzOuDNtwRcuHCBfr+P67psbGywvr6O67qEYUgYhpimSa1W00WrwNPmiTiOcyV/GoQQ4gUlYUSI5yhJErrdLoPBANu2CYKAarVKGIb0ej2yLNMdL5///Oc5HI756PQkcWnxko7DD961wIULFzg8PMQwDJaXlzl9+jS+75MkCbPZjLIsaTQaVCoV6vX6sTqRsiylTkQIcV2SMCLEc1CWJd1uVy+/8zyPZrNJHMeMRiPCMCQIAur1Ol/60pfo9fr8Tm+RXlGh6Zn8+DetcLC3w8HBAXmes7CwwC233EIQBBiGwWg0Is9z6vW6DjSmOf/tebROROaJCCGuRxJGhHgORqMRh4eHDIdDgiCg3W4TxzGDwYDpdEqj0aBarc5beHd2+ONDnzPFApYBP/HNa6Sjrm73bTQa3HzzzTSbTVzX5fDwkCRJdJjxfV+3615cJ6ICihBCXE/kTzYhvgJVD3J4eKhPRAAGg4E+EfE8j16vx6OPPsoDvYLP5JsA/ODXLbNQ9PUskUqlwqlTp1hYWCAIAg4PD5nNZvi+T7PZxPM8fQ1ztE7EcRypExFCXLckjAjxZaRpSr/fZ39/n6Io9DXK4eEhWZbhui5BEDCdTrn//vu5MAj5z+mtlBi8+bYWdy8mbG9v69CysbHB6uoq9XpdD0xT7buqlVd1yYRhqOtEVEeNEEJcjySMCPEsiqJgOp2yvb2tTy+WlpY4ODigLEvyPKfZbBKGIQ888AD7h33ui24hKi1uW/T53pfZnD9/nul0SqVSodPpsLm5SaPRII5jDg4OMAyDdrtNtVo9FjriOCbPcwzDkDoRIcR1T8KIEM9AXZEcHBwwGAywLIv19XUmkwlFURCGoQ4i8zqRXf4gPMFBUaHpW/zYG9t0d7cZDAZUKhVarRabm5u0Wi0cx2FnZ0dPbVVdM+p6JssykiQBpE5ECHFjkD/lhHgGURQxmUy4cOECZVmysLCAbdtMJhOiKKJer5NlGefOneOJJ57gc7Mmj2YdTAP+4TetMDu4wMHBAZVKhWq1ytraGktLS9RqNc6dO0eWZTQaDT3m3XEcLMuiKAq9d8ZxHGxbhiQLIa5/EkaEuEiSJMRxzNbWFlmW4fs+y8vLuk5EbdTd39/n7NmzPDoy+FSqClZXaEZ7dLtdXU+yvLzM6uoq7XabCxcuEIYhruuysrICgGEYejmeqhOxLEvmiQghbhgSRoQ4Issy4jjWC+wMw2Bzc5PhcEie5yRJQq1W4+DggLNnz3K+N+Xj6S2UGHzrS1rcEcyLXWF+xdJut1lfX2dhYYFer8doNMKyLDY3N/XXVDtm1AI8VScihBA3CgkjQjxJXZH0+316vR55nrO2tqZrRKIoolarMRgMOH/+PNt7XX4/vZWwsLl1weNv3GbQ7Xb1JNZms8ny8jLLy8vEccz+/j5lWbKxsaGvX9RVjAw2E0LcyCSMCMFTBathGNLtdpnNZjSbTarVKuPxmDiOMQyDNE3Z3t5me3uHT0Qn2M8q1D2Te++qs78z75ypVqtUq1UWFhY4ceIEANvb2+R5zuLiItVqVZ+AeJ6nT2NAFuAJIW5MEkaEYF6wmmUZe3t7hGGI53ksLi4ym810DYnneU8GkW0+P2vwcDovWP2xr18i6e0wmUyo1Wp4nqfbeH3fZ29vTxe9Li8v6xMQ3/cpy/JYwaoMNhNC3IgkjIgbXpIkpGnK7u4ucRwTxzGLi4s6KMxmMzzPY29vj52dHR4flfxJPD/xePtrOixkXXq9HpVKBcdxaLVarK+v02q1ODw8ZDQa4fu+vvKB+eZdy7KkYFUIIZAwIm5wavfLcDjUu2YWFxexbZswDInjmLIsmc1m7OzssNufcF90mgKDr7+pwtcvRnS7XYIgwLIsms0mCwsLrKysMBqN6Pf7OI7D8vIyjuPoQWa+70vBqhBCPEnCiLhhqTqR2WzGaDRiMplQr9fxPE8XrU6nUwAuXLjAaDTmvukJpoXDRt3m7a+Yn5bYto1t2zQaDVqtFjfffDNxHHN4eEhZlrTbbVqtlh5k5nkeaZpKwaoQQjxJwoi4YYVhSJIk9Ho9fQJSr9f1jw2HQ2zbZnd3l8FgwB/2As7nDVwL7v26BqPDffI8x/M8HMehXq9z8803k+c5g8GAJEloNBosLCzoIGJZFoZhSMGqEEIcIWFE3JCiKNJBJM9z3T2jds6MRiNgvpl3MBhw/2HOf0nXAfjv72xRTYd6+Z0a6765uYnjOPqUpVqt0ul0sCyLPM8BcF1XClaFEOIiEkbEDSdNU+I4ZjQakSQJYRjqExHTNBkMBsRxrGeOXOhN+Hh0ihKDe26p8upWrAtW4zhmaWmJpaUl2u02URQxHo/xPI9Go0GtVtOnIK7r6hMYKVgVQoinSBgRN5Q8z3WdiBpkpsa7W5bFcDhkMpmQpinD4ZBef8DvTk4wK21Othy+56Uuh4eH2LbNbDZjYWGBWq3G+vq6foxpmgRBQKfTORY+iqKQglUhhHgGEkbEDUMVrEZRxHQ6JU1TAL0XZjqdMh6PybKMyWRCv9/nD/p1trMqvm3w97+uRa+7S5Ikur6kVqtx8uRJyrJkPB5TFAWe57GwsACgi1RN05SCVSGEeBYSRsQNQwWR2WxGmqakaYrv+2RZRpZlDIdD0jQlDMN5nUjf4L/E82V273jjMtbsgPF4rEe4LywssLq6iud5uhjWcRyazaZu3YV50aoKPlKwKoQQTydhRNwQVA1IFEWkaUqSJLrFtigKBoMBeZ4zmUwYj8dc6M+4bzofbPaWlza51R1ycHCgg8Ty8jKtVovFxUU9j8QwDIIgoNVq6esZQA86k4JVIYR4ZhJGxHUvyzK9d0YFEdOc/9JXbbhlWTIajcjznO5hj49NNolKi1s7Lv/NrSY7Ozu63mNpaYlqtcrGxgZxHJNlmT5l6XQ65HlOmqaUZanfpGBVCCGenYQRcV0rikIXq6ogouo6ptMps9lM13tkWUa32+VPozV2Up/AMfh7X9eku7dDHMdYlsXy8jK1Wo2bbrpJBw01Lr7ZbOI4jr6eUeFFClaFEOLLkzAirlsqKMxmM10Xok4wVLGpumKZzWb0ej0ejap8ZtIE4O/evUA23GMwGGAYBp1Oh3q9zsbGBo7j6BDjOA5BEFCv1/X1TJZl+vRFClaFEOLLkzAirluqWDXLMn11opbTqYJVy7IYDAaMRiMOIvhofwmA//r2Krf6Uw4PD0mShGazSavVYmVlhVqtRpZlRFFEWZb4vk+73dZfQ+2fUTtopGBVCCG+PAkj4rqUJIlu3y2K4lgNx2QyIYoiXNdlf3+fMAyZhDG/NVgjLkxuX3D4K7fa7OzsEIYhzWaTpaUlFhYWWFhY0J0x0+kU3/dpNptYlqUX32VZhmVZUrAqhBDPkYQRcd1RXTFqH4y6ojEMg/F4rN83HA71bJE/DtfYS1zqrsGP3Fllf3eHwWBAEASsrq6yuLjIysqKDhr9fh/f96lWq3oSa57nOuRIwaoQQjx39pV+AkK8kNTwMnWFkue5LlpV71NzRg4ODpjNZjyWtfnMMMAAfvi1dcppj8PDQ1zXZXl5maWlJZaXl4H58LLhcIjjOHrku/oacRzjui6maUrBqhBCXAI5GRHXlTiOCcNQz/hQrbdHW3CLotBXMP3c46OHHQC+67YKp7wZ+/v7ZFnGwsICa2trdDodHMfBMAzdfeM4Dp1OB9M09dK9siyxbVsKVoUQ4hJJGBHXjTRNmU6neux6HMd60qrqcInjmL29vfk01rTgP/ZXSQqDly/YvHkjZX9/X++cOXHiBK1Wi1qtpmtBwjDEdV3a7Ta2bZMkiW4Z9n0fz/OkYFUIIS6RhBFxXTjapluWJWma6omrakFdmqYcHh4yGAxI04xPRCfYjy2ansHffoXNaDik3+/rxXedTod2u63rTUaj0bE6EfU1wzDE8zxc18V13Sv9UyGEENccCSPiuhBFEWEYkue5HnSmgohhGCRJwnA4ZG9vjzzPebRY5HNDDwP4gVd6FLMBh4eHeJ7H6uoqKysrul3Xtm36/b4OG83mfA6JCiKqWFUKVoUQ4vmRMCKueUmS6AmrhmHopXVFUeirlPF4zLlz5wDoFRV+52AeKP7qbS4rxlCPhF9eXubEiRO0221M08S2bUajkW7V7XQ6+muq+hHf96VgVQghvgoSRsQ1Lc9zvY1XXc8kSUKWZdj2vFlsPB5z/vx58jxnlpb89miDtDB45aLF3a0Jk8mEMAzpdDpsbGzQ6XRwXVcXrKrTkXa7jWVZFEXBZDIhTVM8z6NarUrBqhBCfBUkjIhrVlmWui4kyzLKstQnJGVZ6jbcbrfLcDgkzwv+OD1JNzLo+Ab/7amYKJwxmUxoNBpsbm6yvLxMEATAU0HHdV0ajQau6x4bmmZZFrVaTQpWhRDiqyRhRFyz4jjWYUSNeVcnI5ZlMR6P6fV6bG9vA/BwscIXBg6WAf/dS0tIpkwmEzzPY2VlheXlZRqNBmVZ6oJVz/OoVCrUajVgXicynU4xDEMHFCGEEF8dCSPimqTadMMwxDRN0jQlyzKSJNF1IqPRiK2tLcqypGfU+d3DJ+tEbrVo5wP92KWlJdbX12m1WgDYts1wOMR1XWzb1u9XtSdlWeqOGiGEEF89mcAqrjnqekZ1zwB6nohaUheGIU888cR8Dojh8NvDdfISXrNs8iq/pwtcm80mKysrLC4uYlkWtm0zHo8xTRPLsuh0OhiGocfH53l+7KRECCHEV09ORsQ1R008VScb6kQkiiL939vb24xGI9I04w+TU/Rig+XA4DtXRqTpPIjU63XdxquGlamNvpZl0Wq1sG2boigYj8ekaYrjONTrdUxTfusIIcQLRf5EFdeUJEn0dYlt2/p6ZjabAfPrm263y97eHkmS8Lh7igdGDrYJf+NUBGmoTzcajQYrKyvU63Usy9KL7hzHoVar4fu+LlhV80Sazabu0hFCCPHCkDAirhlFURDHMdPpFEBf0UynU9I0BZ5q451Op0z8Zf7TQR2Av3q6pFXOTzeCIKBarep5IqobZjKZ4DgOvu9Tr88fd7RgtV6vS8GqEEJcBhJGxDVDDTNTI9/VBl51PRNFETs7O/MBZm6V3x6uUWDw+hWDr60MiaKIarWK4zgsLy/rBXimaTIej7EsS1/PGIZBHMdSsCqEEC8CCSPimqC27o5GI2zb1ovvptMpeZ5TliXdbpeDgwPKEv5zcppharJeM/mOlTGz2ZQgCHAch4WFBVZWVgiCANu29RWPZVl68mqWZYzH42MFqzLYTAghLg8JI+KqpwpU1emFWlynho+Zpkmv12N/f58wDHnEPsljUxfXgu85HZOGEyzLIggCarUaq6ur1Go1bNsmiiLSNMU0TX0NoyasJkmC67pSsCqEEJeZ/AkrrmpHp6yq1l3XdRmPx0wmE90Bs7+/z2g0Yuwu8IlhG4C33mpQzYZkWaYHlK2trdFoNPA8TxesqqBSrVYpy5LpdEoYhti2Tb1el4JVIYS4zORPWXFVi6KIPM919wzAbDYf4a429B4cHDAYDEhKk9+dbFBg8LpVi5e6PabTWAeRxcVF6vW6nhEym82wLAvP82g0Gvp9qmC1Wq3ied4V+96FEOJGIScj4qql2nYHgwGO45DnOaZpMhqNmM1m+t97vR5RFPGn2Sn6qc1ixeA7lkdEUYjnedRqNWq1Gu12WxenzmYzDMPQ7bqmaepOHVWwqnbUCCGEuLwkjIir0tE2XrWNNwgChsMh4/EYmLfi9vt9ZrMZZ1nm/lkNA/ieWzLS2QjDMGi329i2zfLyMs1mE9d1CcNQ759Rc0PSNGUymZBlGb7vEwSBFKwKIcSL5HmFkfe///2cOnUK3/e56667+NSnPvWcHvdrv/ZrGIbBd33Xdz2fLytuIKpGRHW6VCoVHT6SJNHj2UejEYexyccnKwB8xymDRnKgg4hlWSwvL+vrGbXhF9CDzfI8ZzKZEMcxruvKJl4hhHiRXXIY+dCHPsS9997LT/7kT/LZz36WV73qVbz5zW9mf3//yz7u7Nmz/P2///f5hm/4huf9ZMWNQbXxDodDHMcB5oWsg8FAh5MwDBkMBowmMz4enyIpTF7ShDvcvWP1Hu12m2q1SqvVIssy0jTFMAzdrqsKVqMowrZtarWa/ppCCCFeHJccRn7u536Ov/23/zZvf/vbefnLX84HPvABgiDgF3/xF5/1MXme8zf+xt/gp37qp7j55pu/qicsrm95nus2XsMwSNOUWq1Gr9ej3+9TFAV5ntPv95lMJny+2GQn8Qjskr+8NoKywPd9Go2G/me73dZDzABc16XRaOjaEVU/Uq1W8X3/Cv8MCCHEjeeSwkiSJHzmM5/hnnvueeoTmCb33HMPn/zkJ5/1cf/z//w/s7y8zN/6W3/rOX2dOI4ZjUbH3sT1T7XxxnFMmqYURUG1WtVFqqpodTgcMp1OeSKu8OnpvI33r5wIqZQRnufpdtx2u02z2aRSqRBFEUVRYFkWjUYDy7KIoojZbEZZlgRBIAWrQghxhVxSGDk4OCDPc1ZWVo69f2Vlhd3d3Wd8zCc+8Qn+xb/4F3zwgx98zl/nPe95D81mU79tbm5eytMU1ygVQiaTCYZh4DiObt1Vk1ejKGI8HrM/nPH74U0A3LWQcNoZ4TgOjUZDB5F6vU69XieOY4qiOLZfJk1TvdPG932q1aoUrAohxBVyWbtpxuMx3/u938sHP/hBFhcXn/Pj3vWudzEcDvXbuXPnLuOzFFeDNE319YxpmuR5ThAEHBwccHBwQFEUmKbJYDBgMBjyh8kpprnFspfzLZ0B1WpVb99VrbytVktf6xzdL5PnOdPplDiO8TyParUqBatCCHEFXdLQs8XFRSzLYm9v79j79/b2WF1dfdrHP/bYY5w9e5a3vOUt+n1FUcy/sG3z0EMPccsttzztcZ7nybCpG4hq4w3DUA8yazQaDAYDut2uXnA3HA4ZDAZ8PmxzJq5iGSXfvT6iUa3gui62bVOpVGi1WjSbTSzLIkkSimJeR6IKVmezGWEY4jgO1WpVNvEKIcQVdkknI67rcuedd3Lffffp9xVFwX333cfdd9/9tI9/6Utfyhe/+EU+97nP6bfv/M7v5E1vehOf+9zn5PpFAOg6kTAMgXkYLYqC3d1der0eruuSJAmj0YitUcEnZ/Nrwr+0MmWzYVGpVLBt+9j1jDoBUePj6/U6hmEQhqGuPQmCQEKvEEJcBS55HPy9997L933f9/Ha176W17/+9bz3ve9lOp3y9re/HYC3ve1tbGxs8J73vAff93nFK15x7PGtVgvgae8XN6YkSUiSRHe0AARBwNbWFt1uF8uyME2Tfr9Ptz/kvvA0OQa3VSNe14poNBb1dUuz2aRWq1GtVimKgjRN9bWN4zg6iKjC2EqlInUiQghxFbjkMPLWt76VbrfLu9/9bnZ3d7njjjv46Ec/qotat7a2ZMOpeE7UorowDCmKgqIoqNfr9Ho9dnZ2CMOQTqfDZDKh2+3yh5MV+rlHzcr5rvUpq6srjMdjqtWqHveuOmlUG69q102ShDAM9STXIAjk16kQQlwljLIsyyv9JL6S0WhEs9lkOBzqhWbi2qZqN9SG3LIs8TwP0zR58MEH2dnZ0bVDFy5c4M92Mz42vQko+b7NIW+4dZHxeKw/ZmNjg3q9TqPRIM9z4jjWRa1FUTCZTJjNZrp2ROpEhBDi8nuur9/yV0NxRcRxrGtFAF18ev78efb397EsiyAIODw8ZOtgwsdnGwC8sT3lzs06SZLgOI4e965ChvrcnufpuSFhGBJFEY7jUKlUJIgIIcRVRsKIeNGlaUoURURRBIBhGPi+z8HBAefPnyfLMur1OtPplP3uAR+bbJCUFutuzLefnNeUJEmCZVl0Oh19AuK6LnEcY9s2QRBg27auEzEMgyAIZMKqEEJchSSMiBdVnuf6pKIsS93tUhQFjz32GJPJhFqtRpZlbG9v88fDJnt5Fdco+G9OTNhYW+Xg4EAPKmu321QqFarVqh5uVqlU8H3/WD1KEARSsCqEEFcpCSPiRVMUxbEgUhQFjuPgOA5nzpyh3+9jWRae57G3t8djI4PPxssA/FdLQ15xeo3Dw0NqtRqmabK6uorjOHoJXhzHOogkSUIURSRJQqVSoVKpSMGqEEJcpeRPZ/GiKMuSMAxJkkTP/zAMA8/zjl3PtNtter0e+4MJHxuvU2LwiuqEb3lJkziOMU2TsixZX1/HcRwWFhb05z46LE9dA/m+r+eQCCGEuDpJGBEviiiKyLKMJEmA+SmJqvF4/PHHCcOQVqtFFEXs73f52GCJaenSthK++/R8LkgURdi2zeLiop60qqasmqaJ53l6noiasOr7vgw2E0KIq5yEEXHZxXFMlmV6wqrqhDFNk8cee0xPWS2KgsPDQz47rHAmbWJS8lfXR9y0vsJgMMDzPHzfZ2Fh4diemSRJ8H0fx3H0iYhpmvpURAghxNVNwoi4rNRpSBzHGIZBmqbYtq3rQnZ2dijLEtM0mU6nnDmM+MPpfIDeN7b6vPaWFSaTiR73vra2huu6NJtNXYOirmFU3UhZljqISMGqEEJc/SSMiMtGhYUsy8jznCzLgPnumdFoxJkzZ3SLrmEYdHtDPjJYIcfkJnfKd7ykqjtubNtmeXlZX8+UZUkcxziOg23blGWpR8urglXZxCuEENcGCSPislATVsuyJE1TPe7d8zziOObs2bMMBgPKssR1XabTKb93UOMw96kYGf/tqZhGvaYLU1utFq1Wi1qthmVZ+vOpAWbqVCQIAlzXxXGcK/wzIIQQ4rmSMCIuCzXiPUkSiqIgyzIcx6EoCra3t7lw4QJZluG6LlmW8V/2c74QtgH4zpUBt2ws64LVSqWi60TU0DLVKVMUxbFuGlW0KoQQ4tohYUS84KIoIs9zfSKSZRmWZeE4Dt1ulzNnzhwLD+d7E377sAPAa6p97j7V1CcfQRDoIOJ5HoZh6HHvRVFgGIYe9S5BRAghrk0SRsQLKkkS0jTVNSKqTkSNez9z5gzj8RjXdedj3dOM39iuE5U2i1bIX77ZwnVdXfvRbDZptVr4vo9pmno+CcyvgrIswzAMXNeVwWZCCHGNkj+5xQtGbctVRadZllEUBb7vMxwOeeKJJzg8PMRxHJrNJpZl8VuPJ5xPq9jk/DcbYzqthp7M2mg0WFxcxHVdbNvGNE2yLMO2bdI0xTAMsizTw86kYFUIIa5NEkbEC0J1zgDHrmgcxyGOY86dO8fu7i6maeqldp95oscfjeZ1It/a6XPbWguYL86r1+ssLCzo+SFquJmaJaIGpvm+j23bsolXCCGuYRJGxFdNFZCqfTPqmsayLEzTZGtri+3tbfI81223Owd9fn27QYnBS9wh33Ry3opbliVBENButwmCQNeaqM8Xx7EOIqp9V+pEhBDi2iZhRHzVoijSXS1pmuo6Edd1OXfuHDs7O0ynU6rVKrVajSiK+dBjJqPCpWbEfPepDMex9ZK8TqdDrVbDtm0cx8EwDP25Ab3p1zRNGWwmhBDXAQkj4quiRr0DejR7WZZ4nsfu7i57e3sMh0NqtZouMP0PD414KKpjUPCXV4e0a/P9MbZt02639TWOmrqa5/mxrhz1popahRBCXNvkT3LxvKVpqhfflWVJFEXA/ESk3++zvb3N3t6ervvwfZ8/O9vj9/vzOpFvqB3w0kVXj3NvNps0Gg1c19UdMurqR13PALpGRDbxCiHE9UHCiHhe8jzX4cM0TT1t1bIsoihia2uLg4MDHMfBdV2q1Spb+wN+batKgcEtzoBvPmHqTpharcbCwoIOIur6RXXoqOChQohs4hVCiOuHhBFxyVTBKoBlWcxmMz2ADNCj3pMkIQgCKpUK0zDi/3rIYFrYtIyQ774pxvNcfSXT6XQwTVOHF1XMqgKPZVm6TkQKVoUQ4voiYURcMtU5Y5qmHnIG4DgOTzzxBP1+n9FoRKPR0PUdv3b/jPNJBYeM717rU/Pns0Ycx2FpaUnvkzEMQ4eNo6ctateMFKwKIcT1R8KIuCRq1LsKBLPZDJhfn1y4cIFer8doNNKbdYMg4PceOuRTowYA39rc50TDodPp4DgO7Xb72OTUarWKaZq6Q8c0TT0GXgpWhRDi+iR/sovn7OJTkNFoBMyDyMHBAfv7+/T7fV3P0Wg0eOBCj9/cqQPwSnef1ywZtNttHMdhYWGBRmM+cdWyLIIgOHbaUpalDiCyiVcIIa5fEkbEc5JlGXEcA/PwMRqN9BXKaDTi/PnzDAYDTNPUdR398YxfetQlLU1WGPJt6wm1Wk0PNWs0GpRlqVt4Pc87NkDN9319zSMFq0IIcf2SMCK+oqIodCGp4ziEYah3w6Rpyvnz5xmPxxRFcaz24xf/POEwdaiUMd+10qdWDVhcXKRardJqtbAsC8MwME2TarVKWZb62ke176rOGiGEENcvCSPiyzp6UqE6XFRgMAyDCxcuMB6PSZIE3/cxDIMgCPh3D4x4YFrBKAv+Yu0ca+0qq6uruK7LwsICjuPojbu1Wk1fzxRFAaADiBSsCiHE9U/CiPiyVCGpYRjYts1gMADmrbY7Ozv0+32m0ym+75PnOdVqlc9sDfjoXgDAa+0tXrEasLy8jO/7tFotHUQsy9IDz7Is0wPUVABR1zRCCCGubxJGxLM6Ourd930Gg4Fu6T08POTg4IDpdIrrumRZRrVaZWcw418+5lJicFO5zzeuliwtLdFoNAiCQNeFqL0zQRDowWYwvwZSPyYFq0IIcWOQMCKe0cUnFZPJRAeTyWTC7u4u0+kU27Z1sWmWF/zC/SXT3KJRTPivlvqsrCyzvLysx8GrTb6WZVGv18nznDzPdcjxPE8GmwkhxA1Gwoh4mov3zCRJogeQZVnG9vY2k8lEf7zqhvnVByK2Qge7zPjW4CybaytsbGzoThk1QVWdiKjrnyRJMAxDBxEpWBVCiBuLhBHxNFEU6YJV0zR1G29Zluzu7jKbzXQ3TVEUeJ7HH22FfKI7X2T3ButRbl9vs7m5eawjRp14qPcZhkEcx7p917ZtGWwmhBA3IPlTXxyTpqm+jnFdl8FgoK9SDg8PGQwGRFGE4zjkeU4QBJztx/zq4/NC09vyJ3jdmsvp06epVqtkWUatVjsWSBzH0acsRwON53myiVcIIW5AEkaEVhTFscFmk8mENE1J05Qoiuh2u4RhiGmapGlKrVZjEmf8/JcK0tJkMT/kW5amnDx5kk6noxflHW3f9TxPd9IURUFRFPi+rxfkCSGEuPFIGBHa0euZLMuYzWbEcUxRFOzs7JAkCZZl6RORLMv4pfszDhIbv4z4lmCL06dOsrGxQZZlulakWq3qrbtZluE4DnEc664a13WlYFUIIW5gEkYEMN87oxbgmabJeDzWpyS7u7vEcUyapuR5riesfmwr54sDG6Ms+HrjQV5+8yYnTpygLEvSNKVer+O6rt45k+c5lmURRZEeZOb7vgw2E0KIG5yEEXHsekYtwFOnJL1ejyiKjrXxuq7L/Qcp/+7sPEC8qniE155a4JZbbsHzPGazGZ1OB9M0abfbpGmqi1LjONZ1IkEQHNvYK4QQ4sYkrwJCt/GapslsNiOKIn1NMxqNGI/H+L5PkswX3XUnKf/iASgxOJFv88blgpe97GVUq1Wm0yn1eh3LslhcXCRJElzX1YWxpmlSliVBEOD7vhSsCiGEkDByozt6PVMUBZPJ5NipyGw2w3VdwjCkXq8zDWN+/ks509ykUYz5lvoeL3vZS6nX68xmM3zfx/d9FhYWjtWHqJ0zAEEQ4LquFKwKIYQAJIzc0I5ezxiGwXg81vUc3W5XDzpL05RKpUKapnz40YwnphZ2mfIXrId5+e0vYWVlRX+eWq1Gs9nENE1M0ySOY/I8B9ADzdQ0ViGEEAIkjNzQjl7PjMdjZrMZeZ5zcHBAGIZ6sJlhGFiWxR9fSPjDvfm1yuuLB3jVLeucOnVKzwxpNpvU63WCINDFrkmS6GCiluLJhFUhhBBHSRi5QanrGZiHktlsRlEUOpSMx2M8z9OnIo92Z/zaY/NfLrdlZ3jdus/p06f1FNV6vU6tVqNerxPHsa4/sW0bwzAIgkBv6RVCCCGOkjByAzp6PZNlGePxmDRNdbHqeDym2WwSRRGNRoO93oh/8ZBJWpos5Ye8sdHn9ttvp9VqEUURvu/TaDSo1+v6JGQ0Gumx79VqVZ+MSAuvEEKIi0kYuQGFYQhAnud6nshkMmEymehumMlkQq1WYzye8KuP2xwmFpUy4k3+WV56+20sLS0xnU7xfZ9ms0m1WtVBYzgc4rouZVnqyatq1ogQQghxMXl1uMGozpaiKJhOp3qGyHg8ZjQa4XkecRxTqVRIkoSPnSv584GFWRa8kQd46ekT3HTTTSRJgu/71Go1Go2G3u47mUz01Uy9Xj+2k0YIIYR4JjLk4QaiCkoBptMps9mM6XTKaDRiOp3ieZ4e0V6WJV/Yi/nI+fkCvDuKR/ia9SonT54EwLIsgiCg3W7jOA5RFOmCWNu29YmI7/tYlnVlvmEhhBDXBPnr6g1EhYUwDHUQGQ6H+tpGdc7Yts32IOSXHzYpMThZ7PCaZsgtt9xCtVrV01NbrRaWZZEkCVEUkee5Hv9uWRae5+E4zpX8loUQQlwDJIzcINT1TJIkjMdjwjDk8PBQb+X1fV8PKdsfTPjA/QbT3KRVjnmDu8Wtt97C4uIieZ5Tq9V0EIF5yEnTFM/zjm3glaFmQgghngsJIzcAdT2T5zmDwYAoiuj1ehRFwWw2IwgCkiShUqnQH035P++HndDEJ+GNPMDNJzdZW1vTY9yr1Sq2bWOaJtPplCzL8DxPv6ltvUIIIcRzITUj17myLPX1TK/XI0kS9vf39e6ZRqOhC1aH4wkfvD/nzMTGJeMvlF/k9o0FNjY29ElHo9HQLbpqSJoKH6o+RGaJCCGEuBRyMnKdU9cz/X7/WBCZTqcEQaCDyGQy5V8+kPHA0Mai4OvLP+elq3VWV1dpt9vYtk2r1dJFrmp/jWEY+L5PpVKRICKEEOJ5kZOR61ie56Rpymw2Yzab0e12CcOQKIoIgoA8z6lUKkynU/7tQzGf6XkYlLyBB3jFSoXV1VVWV1cxTZNOp4PjOBRFQZ7nlGVJWZZ6MZ4KIjLUTAghxKWSMHKdKsuSMAzJsozDw0PdvquKVcuyxPM8oijiPzwy4w+68xON1/MIdyxZnDhxgrW1NVzX1UGkLEsdRABc19UFqzLUTAghxPMlYeQ6pa5ndnd3mUwm9Pt98jzXHTBqSNnHHhny2ztVAF7FGe7spKytnWBjYwPf96lWq7rAtSgKffLhOA6e5+G6rgw1E0II8VWRV5DrUJZlpGnK3t4e0+mUw8NDsizDNE09/yPLMj7x+IBfPx8A8FLO8/r6iLW1NTY3N/F9nyAIaDabxHFMlmX689u2jeM4x2pFhBBCiOdLwsh1RhWXDodDhsMhh4eHJElCWZaYpolt2xRFwZ8+fsCvPBFQYnCKPd5Y7XLixAabm5tUq1Xq9TqdTocwDEnTFADTNHWgCYIA3/exbTlcE0II8dWRMHKdieOYMAzZ3t6m3+8zm82wLAvHcXAcB8Mw+MITB/yrs1VyTNY55E3BBTY3T7C5uamX3rXbbdI0JYoiDMPAsix9AlKtVvF9X4aaCSGEeEFIGLmOZFlGGIZsbW0xGo2YTCY6RLiui2ma3L+1xy886hOXFgvlgHsqT3BiY53NzU1arRbVapVOp0OWZYxGI2zb1tcyavqq7/sy1EwIIcQLRsLIdUJ1z1y4cIF+v89oNMI0zWMnIo+d3+cXHvGZFjaNcsI97qOcPLHOiRMnaLVa1Ot1Wq0WWZYxGAywbVufqiRJQhAEBEEgs0SEEEK8oOTC/zoRhiG7u7tsb2/rxXfqVMMwDM7tdvn5h2wGuUtQRvwl/1Fu2VhjY2ODpaUlqtUqjUaDPM/p9/vYto3ruti2TRRFVCoVqtWqBBEhhBAvODkZuQ6kacrBwQFbW1tEUUSWZXr+h2VZ7HYP+YUHLfYzH69M+FbnQU6vdjhx4gQLCwvUajWazaaeSeI4DpVKhUqlQhzHOI5DvV6nWq3KUDMhhBAvOAkj17iyLOn3+5w9e5bJZEIYhnoiKkD34JBferDkXFLBLjO+2X6Q29fbnDhxgsXFRRqNBvV6XV/NqGV3vu8zmUwwDINms0mtVpMgIoQQ4rKQa5pr3HA45LHHHmMwGBCGoQ4NhmEwHI74Nw9mPBI1MMuCv2A9yNes1XUQUZ0zaZoyHo9xXVcHkdFohGEY1Ot16vW6DDUTQghx2UgYuYZFUcSjjz5Kt9tlNpvpkey2bTMajfh/Ho74wqwNZcndxkPcsRZw8uRJFhYW6HQ6VCoViqLQQUQVu6oumlqtRqvVkqFmQgghLisJI9eoLMt0EJlOp3ieh2maGIbBaDTio49H/Mm4DcBrjUd53arN5uYm7XabhYUFvX13MpnoQlXbtplMJniep+tIJIgIIYS43J7X2fv73/9+Tp06he/73HXXXXzqU5961o/94Ac/yDd8wzfQbrdpt9vcc889X/bjxVeW5zmPP/64HmymTjVs22Y2m/EHT8z4vd48iLyiPMsblgtuuukmOp0Oi4uLuK5LWZZMp1P9OEDXm7RaLdrttgQRIYQQL4pLDiMf+tCHuPfee/nJn/xJPvvZz/KqV72KN7/5zezv7z/jx3/84x/nr/21v8bv//7v88lPfpLNzU3+4l/8i1y4cOGrfvI3ojzPOXPmDFtbW3oWiHobj8d8ZjvkP3TnQeSW/Bzfsppy00030Wq1WFpa0tt3wzDUM0SKoqAoCoIgoNPp0Gg0pFhVCCHEi8Yo1T745+iuu+7ida97He973/sAKIqCzc1NfuRHfoR3vvOdX/HxeZ7Tbrd53/vex9ve9rbn9DVHoxHNZpPhcEij0biUp3tdyfOcra0tnnjiCR3+PM/Dsiym0ykPdiP+1fk2GRYn8l2+a2XA5ua8WHVlZQXXdbEsiziO9UC0PM/10rt2u43v+1f4uxRCCHG9eK6v35d0MpIkCZ/5zGe45557nvoEpsk999zDJz/5yef0OWazGWma0ul0nvVj4jhmNBode7vR5XnO9va2DiJlWeL7vr5uOduL+JXzTTIslvID3rLcZ2NjncXFRZaXl/U01SiKALAsizRNcV2XarXK4uKiBBEhhBBXxCWFkYODA/I8Z2Vl5dj7V1ZW2N3dfU6f48d+7MdYX18/Fmgu9p73vIdms6nfNjc3L+VpXneyLGN3d5fHHntMB5FKpUJZlsRxzPYg5F9u1YlxaOUD/srKIWsryywuLtLpdPR1jJrMapomeZ5TrVap1WosLi7iOM4V/i6FEELcqF7U4RH/5J/8E37t136N3/iN3/iyfwt/17vexXA41G/nzp17EZ/l1SXLMrrdru6cKYpCL6mL45juaMYvn60xLT1qxYTvXu6yvrTA2toanU5Ht+tOp1PKstS1IGp+SKfTkRkiQgghrqhLau1dXFzEsiz29vaOvX9vb4/V1dUv+9j//X//3/kn/+Sf8Hu/93u88pWv/LIfq6aA3uiyLKPf7+sgkue5DnHT6ZTBJOSXz9QYlBUqZcRfbm+z2q6zurpKo9HAcRxc12U2m2EYBqZpYlkWjUaDWq0mU1WFEEJcFS7pr8Su63LnnXdy33336fcVRcF9993H3Xff/ayP+9/+t/+Nn/7pn+ajH/0or33ta5//s72BpGnKYDDgoYceotvt6n0zWZYxHo8ZTqb8yhMB3aKKW6Z8R+0sJzoBGxsbVKtVfSIymUwoy1J33LTbbR1GJIgIIYS4Glzy0LN7772X7/u+7+O1r30tr3/963nve9/LdDrl7W9/OwBve9vb2NjY4D3veQ8A/+v/+r/y7ne/m1/5lV/h1KlTurZE/c1cPF2apoxGIx555BH29/dJkgTXdYnjmDAM6Q5n/EZ3ie0swCpz/lL1DDcvzoNIrVbTI93DMNRdM7Zts7CwoBfgSRARQghxtbjkMPLWt76VbrfLu9/9bnZ3d7njjjv46Ec/qotat7a2jtUg/PN//s9JkoTv/u7vPvZ5fvInf5L/6X/6n766Z38dSpKE2WzGo48+eiyIzGYzJpMJ2+OM3+yvMyo8HDK+NTjLSxc9NjY29NWM53nMZjNs29b/vbi4iOd5VCqVK/0tCiGEEMdc8pyRK+FGmTOSJAlhGPLII49w4cIFfbIxnU4ZDoecmVh8dHITMTZVYr69/gQnWy6rq6ssLCxgWRb1ep0kSfA8D9d18X2fhYUF/e9CCCHEi+W5vn7LbpqrRBzHRFHEY489xvb2NmEYYhgG/X6fwWDAn0+r/FFykgKTRWPCdy3scmKxTRAENJtNDMOgUqkQRRGVSgXf9wmCgFarpbfxCiGEEFcjCSNXgTiOieP/t707j426/PMA/p5p5746vWY6nZ7IJUIFFhswrr8fEigeFHU9WNcjnnEhkVU3mGy0GpOfB8Y/NAT9JQK67Kr4i0CirgQqoCCXtIRCL1qHWmhLmYG2w0yncz37B+lXC52WqbTfzvT9Spq0M8/3yfP00yfz7vfsQ0tLC86ePQu/349oNAq3242LF7twNGjH8UguAKBA4UF5Tg9sGZkwmUzS/UY0Gg1CoRCMRqMUREwmk7SHhIiIaLxiGJFZIBBAMBhEa2srWltb4fP5EA6H0dnZifOei/g5UowWZAMAZqa2oywnCKs1HRaLZUAAAQCz2SydoGowGKDVahlEiIho3GMYkVF/EOm/zfulS5fQ09MDt9uNc91+/CSm44IyDQpEcav6N/y1QAut1iLdVbX/BmhqtRoWiwVarVbaK6LT6XhXVSIiSggMIzLp7e1FKBRCe3s7XC4XfD4fPB4P3G43WruC2J86C36lASqEUWZsxZxcA9RqNTIyMqBQKKBQKKRLeNPS0qTwodfrodfrkZrK0hIRUWLgJ5YMent7EQ6H4Xa7cfr0afh8Ppw9e/bybd97lDiqm4uwUg0jArg3owNFGXqo1Wrp4YIpKSnQ6XTSCar9V8no9XoYDAakpKTIOT0iIqK4MIyMsf4g4vF40NTUhJ6eHjQ3N+P8+fOo7TWj3jQHQqFElsKLf7F5YNWroNVqYTKZpAfk9V9BY7VaEYlEAAAGg4FBhIiIEhLDyBgRQqC3txeRSAQXLlxAc3MzPB4PGhsb4blwAVWhXJwxzwAAFKVcQLndC506FWazGXq9HiqVSrprrdVqhdFoRCgUAnB5j4jRaOQD74iIKCExjIyBPwaRnp4eNDU1oa2tDfX19eg470G1+iZcMBcBAGarz2Fhlh8pKUqYzWZotVrodDpYLBbo9XrpYYV9fX1ISUlhECEiooTHMDLKrgwijY2NcLlcqK+vxxl3N2osC+DTZUOJKP5iaMdNhktQKjWwWq0wGAywWn+/sZnFYkE4HEYkEpEu4dVoNHzODBERJTSGkVEkhJBuYObz+dDY2IiGhgbU19fjV08v6rMWIqg2Q4Uw7ra0IVflg1Z7+WZl/eeE6HQ6ZGZmQqVSIRwOS7d112q1PD+EiIiSAsPIKIlGo+jt7UU0GoXf78fJkydRW1uLuro6nOpJwa85SxBJ0cCEAO62tMIk/DAYLp8L4nQ6YTAYpL0h0WhUOiTDO6oSEVGyYRgZBeFwGIFAAEII+Hw+HD9+HDU1Nairq0NdMB1nc28DFEpkK7xYoj8NPQQsaWmwWq1wOBwwGo1IT0+HRqORQkj/82V4SIaIiJINw8h1FgwG0dfXJwWRqqoqHD9+HLV1dahLnQy3Yw4AoEjpxj9rW6FWKmCxWGG32+F0OmE2m2E0GqXwkZqaykMyRESU1BhGrqNAIIBQKIRoNAqPx4Njx46hvr4ex07Uos48D960GwAAMxWtmJ3aDq1Ki4yMDDidTjgcDphMJumZMikpKTwkQ0REEwLDyHXwxytmAoEAfv31V9TW1qK1tRVVtU1oyF6IXr0NChFBqaIJk1MvwmAwwmazIT8/H3a7HSaTCUajEampqUhNTYVGo+HlukRENCEwjPxJkUhEOlHV7XajpqYGLpcLp0+fRtN5PxocSxFSm5Ea7cP8cA2KjBGYzRYUFBTA6XQiMzMTFosFOp0OSqVSOjRDREQ0UfBT708Ih8PS7d0bGxtRW1uL5uZmnHK1oEU/FedyZ0IoU6EN9WB+6BjyrVpkZmZh8uTJyMrKQlZWFiwWC5RKJdRqNdRqNU9QJSKiCYdhZIT6+voQDAbR3d2NqqoqNDc3Xz5J1W9Eh/1uhFVGAICltx3zwjUosGciNzcXkydPhtV6+YTV/qtltFotD8kQEdGExTASJyEEAoEAgsEgXC4Xjhw5gvr6elT/dhFttgUIpOcAADQhL6YG6lCk6kZ+cT6cTieKi4thtVqRmZkpnReiUqlknhEREZG8GEbi0H8jM6/Xi19++QXV1dWobnChUTMNlyYtBAAoIkHkeWsxS+OGzZ6O3NxZyMnJQW5urnRXVd4zhIiI6HcMI9coEonA7/fjzJkz2LdvH37cfwA1IRu6HPcAKZcvv03vbsSMSDNuyM1CQUEJsrKyYDKZkJWVBbPZDKvVKh2aISIiossYRq5BMBiEz+fDoUOHsKuyEvtO+9DpWIiozgoA0F5qw/Tek7jRbkBR0c2w2+3IyMiAxWKBWq2GyWRCeno6D8kQERENgmFkGIFAAO3t7di6dSt2/VKPRlMJwpMKAQDKQA8Ku49hdkYEhVMnIT8/HzabDWlpaVCpVFCpVLBYLDCbzTwkQ0REFAPDSAz9t3M/ePAgPtuyDYd6s9GXdw+gUAKRILLOH8NcwwVMnVWM4uJi5Obmwmw2Q61Ww2AwSHdS5VUyREREQ2MYGUQkEkFbWxs+2/y/+MdxDy46/gKk6wAA+vO1uFlxGjffWIDp0/8JeXl5MJlM0Ol0MJvN0Gq1vGkZERFRHPipeYVQKITKykq8vfn/cNoyGyi4CQCg7GrFVF8Nbpueizlz7pEeatf/PBk+Q4aIiGhkGEb+wO1247/WrsOOThPgXHz5xd4uONy/oGxqOm699T44nU5YrVZpLwgRERH9OQwjuHx+yP/8Yzve3FqFgHMeFDYlRDgI49nDKJ+sRdn996KwsBBWqxUGg4HngRAREV1HEz6MuC9cRNmqv6HTNg/K/FIoACjOHMOyvDD+7T/uxQ033ACz2czzQIiIiEbJhP6EvX/VazgSzoMy/y9QAoh6fsO81Ba89p8rMGXKFOh0OrmHSERElPQmbBjxB4I4ghugTMtAtLcbuZ4q/Pcb/47iokIehiEiIhpDEzaM6LVq3Bish8sVxd9X34fb5/+r3EMiIiKakCZsGAGA7//+N7mHQERENOHxeAQRERHJimGEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrBhGiIiISFYMI0RERCQrhhEiIiKSFcMIERERyYphhIiIiGTFMEJERESyYhghIiIiWTGMEBERkawS4qm9QggAQE9Pj8wjISIiomvV/7nd/zkeS0KEEa/XCwDIy8uTeSREREQUL6/XC4vFEvN9hRgurowD0WgUbW1tMJlMUCgU163fnp4e5OXlobW1FWaz+br1O15NpPlyrslrIs2Xc01eE2W+Qgh4vV44HA4olbHPDEmIPSNKpRJOp3PU+jebzUn9x3CliTRfzjV5TaT5cq7JayLMd6g9Iv14AisRERHJimGEiIiIZDWhw4hGo0FFRQU0Go3cQxkTE2m+nGvymkjz5VyT10Sb73AS4gRWIiIiSl4Tes8IERERyY9hhIiIiGTFMEJERESyYhghIiIiWSV9GFm3bh0KCwuh1WpRWlqKw4cPD9n+q6++wrRp06DVajFz5kx89913YzTSP+ett97CvHnzYDKZkJ2djeXLl6OhoWHIbTZt2gSFQjHgS6vVjtGIR+7111+/atzTpk0bcptErWthYeFVc1UoFFi5cuWg7ROtpj/++CPuueceOBwOKBQKbNu2bcD7Qgi89tpryMnJgU6nw6JFi3Dq1Klh+4133Y+FoeYaCoWwZs0azJw5EwaDAQ6HA4899hja2tqG7HMka2EsDFfXJ5544qpxl5WVDdvveKwrMPx8B1vDCoUCa9eujdnneK3taEnqMPLll1/ixRdfREVFBaqqqlBSUoIlS5ags7Nz0PY///wzVqxYgaeeegrV1dVYvnw5li9fjhMnTozxyOO3d+9erFy5EgcPHsTOnTsRCoWwePFi+Hy+Ibczm81ob2+XvlpaWsZoxH/OjBkzBox73759Mdsmcl2PHDkyYJ47d+4EADzwwAMxt0mkmvp8PpSUlGDdunWDvv/uu+/igw8+wEcffYRDhw7BYDBgyZIlCAQCMfuMd92PlaHm6vf7UVVVhVdffRVVVVX4+uuv0dDQgGXLlg3bbzxrYawMV1cAKCsrGzDuzz//fMg+x2tdgeHn+8d5tre3Y8OGDVAoFLj//vuH7Hc81nbUiCR2yy23iJUrV0o/RyIR4XA4xFtvvTVo+wcffFDcddddA14rLS0Vzz333KiOczR0dnYKAGLv3r0x22zcuFFYLJaxG9R1UlFRIUpKSq65fTLV9YUXXhCTJk0S0Wh00PcTtaZCCAFAbN26Vfo5Go0Ku90u1q5dK73W1dUlNBqN+Pzzz2P2E++6l8OVcx3M4cOHBQDR0tISs028a0EOg8318ccfF+Xl5XH1kwh1FeLaalteXi4WLlw4ZJtEqO31lLR7RoLBII4ePYpFixZJrymVSixatAgHDhwYdJsDBw4MaA8AS5Ysidl+POvu7gYApKenD9nu0qVLKCgoQF5eHsrLy3Hy5MmxGN6fdurUKTgcDhQXF+ORRx7Bb7/9FrNtstQ1GAxi8+bNePLJJ4d8YGSi1vRKLpcLHR0dA2pnsVhQWloas3YjWffjVXd3NxQKBdLS0oZsF89aGE/27NmD7OxsTJ06Fc8//zw8Hk/MtslU13PnzuHbb7/FU089NWzbRK3tSCRtGHG73YhEIrDZbANet9ls6OjoGHSbjo6OuNqPV9FoFKtXr8att96Km266KWa7qVOnYsOGDdi+fTs2b96MaDSKBQsW4MyZM2M42viVlpZi06ZN+P7777F+/Xq4XC7cdttt8Hq9g7ZPlrpu27YNXV1deOKJJ2K2SdSaDqa/PvHUbiTrfjwKBAJYs2YNVqxYMeRD1OJdC+NFWVkZPvvsM1RWVuKdd97B3r17sXTpUkQikUHbJ0tdAeDTTz+FyWTCfffdN2S7RK3tSCXEU3spPitXrsSJEyeGPb44f/58zJ8/X/p5wYIFmD59Oj7++GO8+eaboz3MEVu6dKn0/axZs1BaWoqCggJs2bLlmv7bSFSffPIJli5dCofDEbNNotaUfhcKhfDggw9CCIH169cP2TZR18LDDz8sfT9z5kzMmjULkyZNwp49e3DHHXfIOLLRt2HDBjzyyCPDnlieqLUdqaTdM5KZmYmUlBScO3duwOvnzp2D3W4fdBu73R5X+/Fo1apV+Oabb7B79244nc64tlWpVJg9ezaamppGaXSjIy0tDVOmTIk57mSoa0tLC3bt2oWnn346ru0StaYApPrEU7uRrPvxpD+ItLS0YOfOnXE/Wn64tTBeFRcXIzMzM+a4E72u/X766Sc0NDTEvY6BxK3ttUraMKJWqzF37lxUVlZKr0WjUVRWVg74z/GP5s+fP6A9AOzcuTNm+/FECIFVq1Zh69at+OGHH1BUVBR3H5FIBDU1NcjJyRmFEY6eS5cuobm5Oea4E7mu/TZu3Ijs7GzcddddcW2XqDUFgKKiItjt9gG16+npwaFDh2LWbiTrfrzoDyKnTp3Crl27kJGREXcfw62F8erMmTPweDwxx53Idf2jTz75BHPnzkVJSUnc2yZqba+Z3GfQjqYvvvhCaDQasWnTJlFbWyueffZZkZaWJjo6OoQQQjz66KPilVdekdrv379fpKamivfee0/U1dWJiooKoVKpRE1NjVxTuGbPP/+8sFgsYs+ePaK9vV368vv9Upsr5/vGG2+IHTt2iObmZnH06FHx8MMPC61WK06ePCnHFK7ZSy+9JPbs2SNcLpfYv3+/WLRokcjMzBSdnZ1CiOSqqxCXrxrIz88Xa9asueq9RK+p1+sV1dXVorq6WgAQ77//vqiurpauIHn77bdFWlqa2L59uzh+/LgoLy8XRUVFore3V+pj4cKF4sMPP5R+Hm7dy2WouQaDQbFs2TLhdDrFsWPHBqzhvr4+qY8r5zrcWpDLUHP1er3i5ZdfFgcOHBAul0vs2rVLzJkzR0yePFkEAgGpj0SpqxDD/x0LIUR3d7fQ6/Vi/fr1g/aRKLUdLUkdRoQQ4sMPPxT5+flCrVaLW265RRw8eFB67/bbbxePP/74gPZbtmwRU6ZMEWq1WsyYMUN8++23YzzikQEw6NfGjRulNlfOd/Xq1dLvxmaziTvvvFNUVVWN/eDj9NBDD4mcnByhVqtFbm6ueOihh0RTU5P0fjLVVQghduzYIQCIhoaGq95L9Jru3r170L/b/jlFo1Hx6quvCpvNJjQajbjjjjuu+j0UFBSIioqKAa8Nte7lMtRcXS5XzDW8e/duqY8r5zrcWpDLUHP1+/1i8eLFIisrS6hUKlFQUCCeeeaZq0JFotRViOH/joUQ4uOPPxY6nU50dXUN2kei1Ha0KIQQYlR3vRARERENIWnPGSEiIqLEwDBCREREsmIYISIiIlkxjBAREZGsGEaIiIhIVgwjREREJCuGESIiIpIVwwgRERHJimGEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrP4fvVoudEHL1C8AAAAASUVORK5CYII=",
      "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": ".venv",
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}