code/chap14.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 14\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from previous chapters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(beta, gamma):\n",
    "    \"\"\"Make a system object for the SIR model.\n",
    "    \n",
    "    beta: contact rate in days\n",
    "    gamma: recovery rate in days\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(S=89, I=1, R=0)\n",
    "    init /= np.sum(init)\n",
    "\n",
    "    t0 = 0\n",
    "    t_end = 7 * 14\n",
    "\n",
    "    return System(init=init, t0=t0, t_end=t_end,\n",
    "                  beta=beta, gamma=gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Update the SIR model.\n",
    "    \n",
    "    state: State (s, i, r)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (sir)\n",
    "    \"\"\"\n",
    "    s, i, r = state\n",
    "\n",
    "    infected = system.beta * i * s    \n",
    "    recovered = system.gamma * i\n",
    "    \n",
    "    s -= infected\n",
    "    i += infected - recovered\n",
    "    r += recovered\n",
    "    \n",
    "    return State(S=s, I=i, R=r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "        \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    unpack(system)\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t0] = init\n",
    "    \n",
    "    for t in linrange(t0, t_end):\n",
    "        frame.row[t+1] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_total_infected(results):\n",
    "    \"\"\"Fraction of population infected during the simulation.\n",
    "    \n",
    "    results: DataFrame with columns S, I, R\n",
    "    \n",
    "    returns: fraction of population\n",
    "    \"\"\"\n",
    "    return get_first_value(results.S) - get_last_value(results.S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_beta(beta_array, gamma):\n",
    "    \"\"\"Sweep a range of values for beta.\n",
    "    \n",
    "    beta_array: array of beta values\n",
    "    gamma: recovery rate\n",
    "    \n",
    "    returns: SweepSeries that maps from beta to total infected\n",
    "    \"\"\"\n",
    "    sweep = SweepSeries()\n",
    "    for beta in beta_array:\n",
    "        system = make_system(beta, gamma)\n",
    "        results = run_simulation(system, update_func)\n",
    "        sweep[system.beta] = calc_total_infected(results)\n",
    "    return sweep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SweepFrame\n",
    "\n",
    "The following sweeps two parameters and stores the results in a `SweepFrame`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep_parameters(beta_array, gamma_array):\n",
    "    \"\"\"Sweep a range of values for beta and gamma.\n",
    "    \n",
    "    beta_array: array of infection rates\n",
    "    gamma_array: array of recovery rates\n",
    "    \n",
    "    returns: SweepFrame with one row for each beta\n",
    "             and one column for each gamma\n",
    "    \"\"\"\n",
    "    frame = SweepFrame(columns=gamma_array)\n",
    "    for gamma in gamma_array:\n",
    "        frame[gamma] = sweep_beta(beta_array, gamma)\n",
    "    return frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0.1</th>\n",
       "      <th>0.3</th>\n",
       "      <th>0.5</th>\n",
       "      <th>0.7</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.10</th>\n",
       "      <td>0.084693</td>\n",
       "      <td>0.005444</td>\n",
       "      <td>0.002736</td>\n",
       "      <td>0.001827</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.18</th>\n",
       "      <td>0.708623</td>\n",
       "      <td>0.015914</td>\n",
       "      <td>0.006118</td>\n",
       "      <td>0.003783</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.26</th>\n",
       "      <td>0.900780</td>\n",
       "      <td>0.055380</td>\n",
       "      <td>0.011639</td>\n",
       "      <td>0.006427</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.34</th>\n",
       "      <td>0.956888</td>\n",
       "      <td>0.267864</td>\n",
       "      <td>0.022115</td>\n",
       "      <td>0.010191</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.42</th>\n",
       "      <td>0.977045</td>\n",
       "      <td>0.524563</td>\n",
       "      <td>0.047816</td>\n",
       "      <td>0.015946</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           0.1       0.3       0.5       0.7\n",
       "0.10  0.084693  0.005444  0.002736  0.001827\n",
       "0.18  0.708623  0.015914  0.006118  0.003783\n",
       "0.26  0.900780  0.055380  0.011639  0.006427\n",
       "0.34  0.956888  0.267864  0.022115  0.010191\n",
       "0.42  0.977045  0.524563  0.047816  0.015946"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta_array = linspace(0.1, 0.9, 11)\n",
    "gamma_array = linspace(0.1, 0.7, 4)\n",
    "frame = sweep_parameters(beta_array, gamma_array)\n",
    "frame.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how we can plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXecnNV5v309U3e2S9qmLerSUa9IQkIFAzIQDCZgMATDS4wNb+yExPmR/BLbARtj494owcYtsXBwbOyXJkAUIQkQkjAqqOxRA0mr1fY2vT7vH8/s7Mw2zUq7U3bP9fGaOWeembklzc53zjn3fX81XddRKBQKhSLTMKU7AIVCoVAo+kMJlEKhUCgyEiVQCoVCochIlEApFAqFIiNRAqVQKBSKjEQJlEKhUCgyEiVQCoVCochIlEApFAqFIiNRAqVQKBSKjEQJlEKhUCgyEku6AxhuhBB2YDlwFginORyFQqFQGJiBicBuKaU/mQeMOoHCEKft6Q5CoVAoFP2yFngrmQtHo0CdBXjqqaeoqKhIdywKhUKhABoaGrjtttsg+hmdDKNRoMIAFRUVVFdXpzsWhUKhUCSS9NGLSpJQKBQKRUaSlhWUEGIF8IKUsmyA+ycBvwQuBpqAf5BSbkphiAqFQqFIMyldQQkhNCHE54DNgG2QS58G9gMTgM8DTwshpqUgRIVCoVBkCKne4vs68HfAQwNdIISYBVwE3C+lDEgp3wCeA+5KTYgKhUKhyARSvcX3hJTyfiHEpYNcMxc4JaV0x83VAiuGI4BIJEJdXR1ut/vcFyuyGqvVSllZGYWFhekORaFQnAcpFSgpZX0Sl+UDnl5zHiB3OGJoaWlB0zSEEJhMKkdktKLrOl6vlzNnzgAokVKg63qvcdzt/ibj5/t9vn5nhxBP0pcO4VlHDotJw2xO7WdmJqaZuwFHr7lcwDUcT97R0cGUKVOUOI1yNE0jNzeXqqoq6uvrs0agIhGdYDhCMBQhGAoTCulEdJ1IpJ//RnR0HWMcNxfRdfQIhPX4cfQancTroteEI8TGekSPPpboaxgfjzo9H/K6Djo60f8ZH7a6HrudeJ+e8OeL3Y7e1HU99mHd53F6z4dz/PPofW70O+wjSorzx241c+myambWjEvZa2aiQB0CJgkhHFJKb3RudnT+ggmHw1it1uF4KkUW4HA4CAaDI/Lc4W4hiQlKhEAwTChuHAyFCYR6xqFQ2Liuey4YTnh8KBwZkVgVigvFHwxz5GT72BYoKaUUQuwDvimE+HdgNfBJYNVwvYamacP1VIoMJ9l/63BEx+UJ0OUO4Iz+1+UJ4A/2CE1MhKK341cDiuyi9/tC62fQ+52T7HtpwKuG8LGTiZ9ReTlWFs4sTelrZoRACSFuA34mpcyPTt0I/ByjBqoFuEtKeSBd8Smyn0hEx+UN0uX243QHoyLkp8ttzLl9oYzYDtI0DavFhNVswmoxYbGYMGkaJpMW918wmTS0uHmziYRxz397zUfnuseapmE2aWga/TzW+NHo+cDUtLjb0f/T0KLzPX+G3vfF5rXuxxmTprjPYZPWowzdz6FFn6j353WioGh953qmY6+tyD7SIlBSyjeB4rjxU8BTcePTwNWpj0yRLmpra3nggQeQUlJTU8M3v/lNFi5cOOhj9u/fzz333MOOHTvizlKMH+N2xFgZeYM88af9RIZZgEzdYmIxYbWYY7dtFhOWuLEx1zO29DPX/RwWs6Y+TBWKKBmxglKMbQKBAF/4whe444472LhxI5s3b+auu+5iy5Yt5Ocbi+p4AQqFI/zpmT/yox/+AICWTi+R8MDi050gMBiappGXY6Eg10Zhno2CPBsFuTYcdksfEem+bTYpMVEoRhIlUBnM5s2b+cEPfkBHRweXX345J06c4Oabb+aGG25ASsl3vvMdjhw5gtPpZPHixTz88MNUVlbyyCOPUF9fT2dnJzt27KC6uppvfOMb/OpXv2L79u1UV1fz/e9/HyHEkK4NBAJ897vfZdu2bTQ1NVFaWso//dM/cc011/SJ/f777+f555/vM3/ttdfy4IMPJszt2rWLYDDInXfeCcBVV1/Nf//3b3nmz89y3fU3xoSpO0XrFz9/nB3vvMUdd36O//7Nk4OKUzx5OdaY8BTm9fwU5NooyLWmPIVWoVAMzpgXqD2yiV2HGgiGRj57ymoxsWJuBUtEvy0IE/jwww+57777eOSRR1i9ejW//OUveeaZZ7j55psBuPfee7npppv4xS9+gcvl4otf/CJPPvkkDzzwAADPPfccjz32GI888ghf/OIX+cxnPsNjjz3G97//ff7t3/6NRx99lEceeWRI1/7qV7/iwIED/OEPf6CgoICnnnqK+++/nyuvvBKLJfGt9OCDD/YRooE4duwY06ZNx+sP4Q+ECQTDVNVMRsojBIN9/12u/+tP8bm7v8Ce999LmO8+MzHH/ZjMJvIcVv7fG2ZjUQKkUJwXejgEgGZOrWSMeYHae6Q5JeIEEAxF2HukOSmBevHFF1m9ejXr168H4J577uF3v/td7P4nn3ySyspKgsEgDQ0NjBs3jqamptj9Cxcu5NJLLwVgxYoVnD17NvZcq1ev5qmnnhrytbfccgs333wzhYWFNDY24nA4cLlceL1eCgoKhvz3oes6gWCYtg4nZosVpzsQuy8nJwefz5dwvRYVnaqqiZhNGg67GU3TGF+UEzvY7w+TSVPipFAMgh4JE/G6CHudhD1dRLxOwh6n8V9vFxG/F81kJn/+OnKqRcriGvMCtXhWaUpXUItnJZem2dTURHl5eWysaVqCAeOBAwe45557cDqdzJw5E6/Xy/jx42P3FxfHclAwmUwJhaomk4lIJDLka10uFw8++CD79u2jqqqKqVOnAv0XQ37ta1/jhRde6DP/iU98gi9/5T/wBYzVkq6D1WrH7090gA4E/BQV5lOUb8ds1voVILvNePsq8VEoBkfXI0S8bsJeJxFvV5z4RAXJ5+Fc/Sr0SBh/wwklUKlkiShLakWTaioqKti3b19srOs6jY2NADQ2NnLfffexceNGli5dCsBDDz1EfX1PJ6mhHN4ne+0DDzzA5MmTefzxx7FYLBw6dKhfEQJDoL72ta/FYg+FI/gCYfyBMB3ORDGaMnUa//v7p7BYTNitZnJsZs6cPsn111+P3WZO+s+hUIxVdF0n4vdEVz49K6Cwt4uIx0nY5xpab6XeaBqWggnkTl8yfEEnwZgXqEzl2muv5Wc/+xnbt29n1apVbNy4kYaGBsBYyei6Tk5ODgA7duzg2WefZfHixSMak9PpxG63YzKZaGpq4gc/MLLoBurUYIiSsVIKD5DIYDJrrFmziu9r8Nyfnua2227j5ZdfQkrJhg0bRuzPolBkE7quowd8PYITXf1EPNHVkM+FHknaqLYfNEw5uZhzCzE7CjA5CjDnFmBydI/z0LTU71QogcpQampqePjhh3nggQdwuVxceeWVVFZWYrVamT59Ovfeey933XUXoVCIqVOncuutt/LCCy+MaLHpV77yFb761a/y9NNPM27cOG6++WYOHjzIkSNHWLXKaPQR7l4pBcOEBtg2NZk07DZjpWQxm9A0LZbg8ZOf/ITq6moee+yx2JblE088wfPPP8+LL744Yn82hSJT0MNBQs52Ql0thLtaCXW1EHK2oYcvrGWXye7A7CiME5+COEHKRzNl3m6FlgnV88OJEGIK8OHrr79OdXV1n/sPHz7MnDlzUh7XUKmvr8fj8TBjxozY3OrVq/nud7/LmjVr0hhZXyIRPbZSGugsT9OIipJRV5TK+qFs+TdXjD0ifo8hQF2txo+zlbCrg/PpX26y2jF1C05uAWZH/O2ClGfg9aauro7LL78cYKqU8qNkHqNWUBlKU1MTX/jCF/j9739PdXU1v//97wkEAiO+jZcskYiRgeeLpoX3i2Z0QLbbzNitZlXUqhiz6HqEsLuLcLcYOVsIdbUQ8XvP/eAomsVqiE7c6sfkKIiJkMkymEl5dqIEKkNZvHgxd999N7fffjudnZ1MmzaNJ554ItZZIR10p4V3b+EN9CXPZjVht1mwW82YTEqUFGMLY4uuLbYyCkdXRt21ROdGw5xXhKWwBEvhBOOnYAKaPXfMfclTApXB3HnnnbHuCulC13WCoWiyQzCMPkA2vtViMlZKNjNm5bWlGCMkbtEZq6Kwu4tkt+g0swVLwQTM3UJUWIIlfzyaRVkCgRIoRT/0TgsfyFbCbDaRE012UG2CFKMZY4uuM7oiaiHkbB3yFp3J7jAEqMBYGZkLSzDnFaYlOy5bUAKlSEDXdZyeAD5//+dKJrNGjtVIdrBY1C+WYnSi6zqhrmaCzacJtNQR6mgaQhq3hjm/GEtB3KqocAIme+6IxjwaUQKlSMDtC/YRJ82kxVZK3WnhCsVoI+xzE2w5TaC5jmDLaSJB/zkf071FZyksiW7TlWApGIdmVlt0w4ESKEUMnz+Ex9tzkGuzmnHkWLClOC1coUgFejhEsL2BQMtpgs2nCTnbBr3elJPXsyoqUFt0qUAJlAKAYChMl6enWavNaqIo36aESTFq0HWdsLsjtm0XbD0z6LadyebAWlKNrbQG24QqTDl5KYxWAUqgFBjdHzpdgVjikdmsUZhnV+KkyHoiQT/BljPGKqnlNGGva8BrNc2EZVyFIUilNZgLJqjfgTSjBGqME9F1Ot2BWKaeZtIoyrenvH5pKJbvW7Zs4Yc//CF1dXVMmDCBz33uc9xyyy0pjVeRmeh6hFBnM4Hm0wRb6gh2NA7aJNWcW4SttMZYKU2oUundGYYSqDGMrus43YGEnnlFebaU21ckY/neTVNTE/feey+PPvoo69ev5+DBg9x6660sWLCAefPmpTRuRWYQ9roIttRFV0l1gyY3aBYrtglVWEtrsJXUYM4tHPBaRfpRApXBjLTle9Wkqfzn44/R0HAWp7OLv+zeSU1NTdot36+55ho2btzIpk2bYg7C3ZSVlbFjxw7y8/OJRCJ0dHRgNpvJy1PnA2MFPRwi2FZvnCM11xFyDZbcoGEpKsFWUoOttBpLcXlGNkVV9M+YFyjPiX14jr53wZ2Ck0EzW8mdeRG50xad89qRtnz/yU9+ytce+h4Am1/exA9+9BN+/sTjabN8nz59esLctGnTOHLkSL/X5+fn4/V6ueiiiwiFQnz+859nypQpSb2WIvvQdZ2wqz2WbRdsOzt4coPdga2kBmvpJGwlVZhsjhRGqxhOxrxAeT/clxJxAqNHl/fDfUkJ1Ehavq+8+GJ++9sey/d58+dz1ccvR9O0lFq+d+PxeGLeVt04HA683oGr9O12O3v27EFKyd13383kyZO56aabzjsGReYR9nThO3UIX/1RIj73gNdpJnNPckNJDeaC8Sq5YZQw5gXKMXVRSldQjqnnFicYOct3w68pgh5tqqdpUDKh5xc6FZbv3U673eTm5vaxfPd6veTmDlx5bzKZsNlsLFiwgJtvvpnXX39dCdQoQNd1gs2n8Z48QKD5NAP1tDPnFWErnWSslMZPVMkNo5QxL1C50xYltaJJNSNh+d6dsdetJ0aHCMuIW76fi+nTp/Ob3/wmYe7EiRNcf/31fa7dtWsX3/72t/nTn/4UmwsEAgmiqsg+IgEfvrpafCcPEvY6+9yvWWy9khvOf8WuyB5UCXSGcu211/Luu++yfft2QqEQ//Vf/3VOy/eBrNdh4Iy9oaSTD9XyPVlWrlyJruv85je/IRgM8uKLLw5o+T5nzhwaGxv59a9/TTgc5v333+eZZ57hU5/61AXFoEgPwY4mnPu30PbGb3HXvttHnGylNRQuu4oJG+6kcNmVOCbNVeI0hhjzK6hMZbgt3yMRHX+g52DZZNKwWYeWzZSM5fv5YLPZkrZ8Lygo4Oc//zkPPfQQjz76KBMnTuShhx5ixYoV5/36itSih0P4zx7De/IQoc6mPvdrVjs51bMNMcorSkOEikxBWb5nKMNp+e7zh+hy97QxcuRYKMgdfe6bA5Et/+ajnbCnC++pg/hP1/Zbq2QpKsUxeR72idNVs9VRiLJ8H0UMl+V7fz328h3ql1+RGnQ9Ek16ONhv0oNmMmOfOJ2cyfOxFJWq7DtFAkqgMpThsHxXPfYU6SIS8OI7XYvv1KF+kx7MjgJyJs0lp2a2qlNSDIgSqAzmQizfM6XHnmJsEexownfyAP6zx/spptWwldbgmDwPa2mNsqlQnJOUCpQQYhHwBLAQOAF8Vkq5u5/rLgZ+CgigGfi2lPIXqYw1m8mUHnuKsUFP0sNBQp3Nfe43We3Ya2bjqFFJD4qhkTKBEkLYgGeBHwPrgBuBzUKIyVLKrrjrTNHr/o+UcqMQYjmwXQixW0q5r7/nViTi9gUTMvYK8mxDzthTKM5F2N2J99QhfHW16IMmPcxAM6vNGsXQSeW75lLAKqX8cXT8tBDi74FPA0/GXTcOKAM0IYSGcYISAgIozklvV1xHjgWHXX04KIYHXY8QaDqF71R30kMi8UkP1uKyNESoGE2k8pNrLnC411wtsCB+QkrZKoR4FPgv4NeAGfhHKWXvxyp60Tdjz6wy9hTDQlJJD5PnkVM9G5Mtp59nUCiGTioFKh/w9JrzAAkN16JbfD7gb4BngNXAn4QQtVLKzakINBvpP2NPWbYrLozkkh7mR5Me1HtNMbykUqDcQO980lygtwfzDcAlUsp/iY63CiF+CdwDKIHqB5WxpxhuIgEvrkPv4K8/2ue+WNLDpHnK8E8xoqQyresQRlZePLOj8/HUAPZecyEgNZ4YWcZoydirra3l05/+NIsXL+baa69l//79A177xz/+kXnz5rFkyZLYz5///OcURjt60XUd35kjtG99uo84WYrKKFj4McZfdjv5s1cpcVKMOKlcQW3BSHz4EvAoRhbfQqD3J8tm4GEhxN0YyRNLgc8Dn0thrFmD25v9GXtDsXwHOHToEH/7t3/Lfffdl4ZoRy9hjxPXwW19kh/sFdNwTF+MtUglPShSS8oESkoZEEJcjVEH9SDwEXC9lLJZCPFl4DYp5Twp5UEhxA3AN4DvAQ3Av0kpn01VrJnCuSzfv/Xwtzl69Agup4t58xfy4Dceomz85KQt34UQQ7o2EyzfAQ4ePMgdd9wxPH/JCnQ9gvejD/Ac2Y0e7skANTvyyZ+3DlvZpDRGpxjLpDT/WEp5AOjT6VRK+S3gW3HjTcCmVMS0v+Ew79XvJxT3izlSWMwWLqpcyMKKczcuTcby/epPXM93vv8Ibrebr/77/+F3G3+TtOV7t437UK7NBMv3cDiMlJJnn32Whx9+GIfDwU033cTnP/95dUh/HoS6WnF+8GavAlsNx5T55M1aoYwAFWkluw4qRoD9jYdTIk4AoXCI/Y3JZcvHW75brVbuueeemMNuOBzhOz94hJs/fRvBYJDW1iZKS8b3a/luNptZsWIF06dPZ/369dhsNlavXk1dXd2Qr73lllt4/PHHKSwspKmpKcHy/UIYiuV7W1sb8+fP5/rrr+eNN97gpz/9Kf/zP//D7373uwuKYayhh0O4a3fS/vYfE8TJUjCe4tXXkz/3EiVOirQz5is4F5bPSekKamF5crYPA1m+69GMvcOHDvKv/3wvLreLWTNn4PP5krJ87x5327gP5dpMsHwvLS1l48aNsfGcOXP4zGc+w+bNm7ntttv6XK/oS6D1DK4PthH2dMbmNJOZ3BnLcExbhGbKrjNMxcjT7u3EFXBTVViBKYU9FJVAVcxJasst1Qxk+e7zhzlb38A3vvYVHvnPX7B29UpsVnNSlu8DkU2W70ePHuWll17i3nvvjc0Fg0Hs9t6Jn4reRIJ+3Id34KurTZi3jp9I/vz1WPKLB3ikYiyi6zqnO+v5oLGWM12Gm/fCijlcXLM0ZTGM+S2+TGUgy/dgOILH40bXdcYXFWCzmpOyfB8OMsHyvbCwkF//+tf87//+L5FIhAMHDvDb3/6WG2644YJiGM3ouo7/7HHatz2dIE6axUb+/HUUrbxOiZMiRigSprb5OH84+CIvH30zJk4AnqAvpbGM+RVUptLb8v2KKzZQXlGB1WJl8pSp3PN3X+QLf3d30pbvw0EmWL6Xl5fz+OOP873vfY+HH36YcePG8YUvfIGrrrpquP6Yo4qw14Xr4HYCTScT5u0V08ifewmmnLw0RabINHxBH4eaj3Kg6Qi+3kKkaUwbV8PqSctSGtOAlu9CiHXJPomUctuwRXSBjEbL92AoTLvTz3VXX85XH3iINWvWUJSv2hglS7b8mw8nuq7jO3UQt9yJHupZ4Zpy8sifuwZ7xdQ0RqfIJDp8XXzQUMuR1hOEe7WzspgtzC6ZwfxyQaE9ebPU/hhuy/c3e411QAMiQBiwRm8H6NVPT3HhdFu+/8//PI2jsIRn//RHgsEACxcuVD32FIMScrbjOrCVYHtDwnzOpLnkiZWYrOq8bqyj6zpnXU3sbzjMqY4zfe7Ps+Uyv1wwu2QGdostDREaDCZQBXG3PwV8CaMf3ntSyrAQYiHwM+A3Ixfe2GXx4sV87vOf5/bbb6fL2cXkyVP4zvd/SmXFBNVjT9EveiSM59j7eI/vQdd7sjTNecUULFiHdXxlGqNTZAKRSIQT7afY33iYFndbn/tL8sazsHwO08ZNwmRKf4rCgAIlpXR33xZCfAO4Md79Vkq5XwjxReAlDKFSDCO6rnPjTX/DJz75aWNCg+J8e9b12FOkhmB7A84PthJ2tcfmNM2EY/oScmcsVanjY5xAKMDhlmMcaJS4A71NJWBScRULK+YwMb8so3Znkk2SKMTwZepNwRCeQzEE+vTYy82+HnuKkScSCuCRO/GePJgwby0uJ3/BeiwF4wd4pGIs4PS7+KBRUttyrE+tp9lkZtaEaSwoFxQ7itIU4eAkKy5/BH4dbfS6B+MsaiXwQ+C/Ryi2MYvPH8LjU664isHxN36E6+B2Ir7YZgea2UqeWEHO5HloKSyoVGQWTe5W9jcc5kT7KeiVCJdjzWFe2Szmls7EYc1sc8lkP/X+AfhP4Nm4xwSBXwD/MtCDFEMnGIooV1zFoET8HlwH38LfcCJh3lY2mfx5azE7LizbSpGdRPQIpzrOsL/xMA3O5j73FzuKWFg+hxkTpmDJki3fpARKSukF7hRC/AM9nk61UsreZoOKC8TrDylXXEW/6LqOv64WV+276MGe1lAmm4P8uZdgmzhdvVfGIMFwkCOtH/JBYy1dPmef+ysLK1hYPpuaosqse38kvW8khCgAbscQqG8A64QQh6WUH45UcGORYCjx3Ell7CkAwu5OnB9sJdhWnzCfUz2bvNkXY7Jl9laNYvjxBLwcaJIcbj6GP5TYy9KkmZg+fjILKmZTkpu955BJCZQQYhZGXVQXMA34CXAr8EkhxFVSyndGLMIxRDgcIRyOLp80sFrUGcJYR4+E8Z7Yh+fYX9DjiijNuYXkz1+HraRvMbpidNPqaeeDxlqOtX5EJK6cAMBmsTG3dCbzymaRZ8v+8tRkV1A/Af5XSvlPQggngJTydiHEYximgpeMVIBjiUCcbbvNYsq65fiFUFtbywMPPICUkpqaGr75zW+ycOHCPtc999xzMc+rbnw+H6tWreJXv/pVqsJNCcHOJlz7txJytvZMahq5UxeRO3MZmlmdTY4VdF2nruss+xsOJ/TG66bAns+C8tmIkmlYR9H7IlmBWgX8Yz/zPwT29TOvOA/it/esluw4xBwOhmL5ft1113HdddfFxocOHeKzn/0s//qv/5rqsEeUQEsdXe+9lLBqshSVUrBgPZbCkjRGpkgluq5zsuMMu8/spd3b2ef+8vxSFlbMYXJxVUptMFJFsgLlBcqA3hanMzG2/bKWjr37aNv1HpHQyHYCBzBZrIxfcRHFixf1uU/XdQLBuBWU1XROy/fvfOc7HDlyBKfTyeLFi3n44YeprKwc9Zbv3QSDQe677z7+4R/+gdmzZyf7z5Dx9BYnzWQmd9YKHFMXqNTxMUSHt5N3Tv+Fus6ziXdoGlOLa1hYMZvy/NL0BJcikn23/xfwmBCiu2V1qRDiOozU86dGJLIU0bF3X0rECSASCtKxt/8FZziiE4kY50+aBqdPneS+++7jy1/+Mm+99RaTJk1iz549sevvvfdeVq9ezbZt29i+fTuRSIQnn3wydv9zzz3HzTffzHvvvUdVVRWf+cxnuPHGG9m5cyczZ87k0UcfHfK18Zbv77//PnfccQf3338/oVBfs8cHH3yQPXv29PnpzwZ+KJbv8Tz11FPk5OTwN3/zN4Nel00EWs8kiJMpJ4/iNTeRO22REqcxQiAc5N3T7/OHgy8miJPFbGF++WxuWXAdG2asHfXiBMmvoL6Kkfz8BmAHdgAh4HHgKyMTWmooXrwopSuo/lZPYNQ/dWO1mNm0aVPM8h3gnnvuSbA1f/LJJ6msrCQYDNLQ0MC4ceP6tXwHWLFiBWfPno091+rVq3nqqaeGfO0tt9zCzTffTGFhIY2NjQmW7wUF8a0bh8ZQLN+7CQQC/PKXv+TrX//6qDmrC7TW07V7U6I4rbwOc15mVvkrhhdd1zna+iE76/bgjbe70DTmls7goqpF5FjGVqPfZOugQsC/CyG+DkyPPu6YlNIthCgDmgZ9ggymePGiAUUjlQSCPWcNNqtpQMv3bg4cOMA999yD0+lk5syZeL3eMWH53s327dsxmUwxYc12gm31dL2nxGms0uxu5e1T79HkakmYrygoZfWki7I6VfxCSDbNPAxUSCmbgYNx85OAQ4AqXb8AdF3vs4IayPIdoLGxkfvuu4+NGzeydKlhvzxWLN+7ef3117n66qszouPyhRJsq6dz90vo0V5pSpzGDt6gj91n9lHbcjyhJVGuLZeLq5cwffzkUbNDcD4MKFBCiFuBv44ONeAXQgh/r8smA317tiuGRML5k0nDYta49tpr+dnPfsb27dtZtWoVGzdupKHBSC91uVzouh7bFuu2fF+8ePGIxpkKy/fbbruNzZs3D2j53s2+ffv4x3/sL7E0uwi2nY2Kk/F3aLLnUrTiWiVOo5yIHuFQ01Heq99PINTT2sykmVhYMYclE+eNqnTx82Wwr5+vAi6guxOlN3q7+8cF7AQG/praqVGwAAAgAElEQVSrSIqE7b1o/VO85fvq1as5fvw4lZWVWK1Wpk+fzr333stdd93F8uXL+dGPfsStt97K8ePHR9zy/a233mLZsmXccsstLF++nHHjxp0zmeFcdFu+v/LKK6xYsYInnniij+V770zBM2fOUFZWdkGvm26C7Q107t4UJ04OilZehyW/+ByPVGQz9c5Gnjm4iXdOvZcgTjVFldw0/xpWVC9W4hRlQMv3eIQQDwDfk1J64uYs0bOpjCIbLd87Xf6YtUZ+rpXcHGuC5Xs3q1ev5rvf/S5r1qxJV6hZSSb+mxvi9GLMjr1HnMalOTLFSOEKuHn39B5OtJ1MmC/MKWB1zTImFVelKbLUcD6W78lu4P8IeEII8eW4uRNCiF8KIbK/n0Ya0XWdQCg+QcIo0G1qauKOO+7g9OnT6LrO008/TSAQGPFtPMXI00ecbA6KVihxGq2EImHerz/A7z94PkGcLGYLK6oXc9O8a0a9OJ0vyaaZPwLMA34cN3c78H3gB8DfDXNcY4ZQWKe7nZbJpGGONoddvHgxd999N7fffjudnZ1MmzaNJ554ok9nBUV2Eexo6itOK6/DUqDEabTR3QVix+m/4PQnGj9MHz+ZlTVLyLflpSm67CBZgfoE8DEp5f7uCSnlViHEPcDLKIE6bxLbGyX237vzzjtj3RUU2U+wo4nOXS/EiVMORSuvVeI0CunwdfHOqff6dIEYnzuOSyZdxMSC7D4/TRVDsWkdqJ+/bTgCGasktjcaO/33xhrBzm5xMg7FTVY7RSuuVZbso4xAOMj79R/wQWNtQsKSzWJjedUi5pTOGJU980aKZAXqeYxWR38rpTwAIISYg7H1t2mkghspdF3PiNoCo/4pcQWlGF7iC5LTRaizmc5dLyaK08rrsBROSHNkiuFisC4Qc0pmsLx67HWBGA6SFagvAf8fsD9aC6VjtDzajGEHnzXk5OTQ2trKhAkT0i5SoXAkVpsXf/6kuHB0XScYDNLY2EheXvr2+UNdLXTseiHmgGuI07VKnEYRLe423jq1u08XiPL8Ui6ZPHa7QAwHybY6agfWCyHmAnOBAHBESlk7lBcTQiwCngAWAieAz0opd/dz3USMRrQfA3zAz6WU/zGU1xqI6upq6urqaG5uHo6nuyD8wTCBaHq51WKi2T6UHVfFubBYLBQVFVFSkh57ilBXC507n4+JkxYTJ2WXMRrwBX3s6rcLhIOV1UuYMX5K2r8EZztDsXw3AVOAGuDXwCwhRKGUMim7DSGEDXgWIxNwHXAjsFkIMbmf53gW+AtQDkwEtkbt5X/HBWK1WmM95NLNs9uOc7rRCcAVK6qYPVl90xothLpa6dz5PJE4cSpe8QklTqOAwbpALKiYzZKJ87GpQtthIdlefBOBVzAEKhdDQP4NWCmEuEJKeTiJp7kUsEopu1PVnxZC/D3waSDmEyGEWIlhK3+JlDIIfCiEuBSjk8WoIRyOcLbFHRtXl6r08dGCIU7PJYhT0YprsBSNfnuE0U69s5F3Tr1Hm6cjYb6mqJJVk5ZRnFM4wCMV58NQLN8PAcuB7o3WzwC/jd738SSeYy7QW8hqgQW95pYBHwBfE0LcibHF97iU8gdJxpoVNLR5CIWNA/zifDv5uSoZcjQQcrbRuStu5WSxUbT8GqxFKq04mxm4C0Q+q2qWMamoSm3njQDJCtRlwDoppV8IAYCU0iOE+ApGP75kyAc8veY8GCuyeMYDa4GtGCup2cDLQoizw7HFlymcaeop3KsuU6un0UDI2WZs6wWMLC7NYqVoxSewFitxylYieoR9DYfYU3+QUKSns5vFZGFJ5TwWlM/BYlLlISNFsgKl0X8dVAlGwkQyuAFHr7lcjKaz8fiBLinl16LjfUKIXwA3AKNGoOriBKpKCVTWE3K2R8XJ2IlW4pT9eIM+Xj/xNvVdDQnzqgtE6ki28OYZ4AfRsygdQAixEHgM4zwqGQ4Botfc7Oh8PLVAbjSpoptRld4WDIVpaOs5f6pS509ZTcjZTueu5xLFafk1WIvLz/FIRabS6GrmT4deShCn8bnFfGL2FVw+fY0SpxSR7Af/P2MkMpyJjg8BVowC3n9O8jm2AJoQ4kvAoxhZfAuBP/e67lWgGUMQ/w+GqN3FKGqndLbFHfN/mlCYQ26OyvjJVkKuqDj5o+JkjorTuIpzPFKRiei6zqHmo+w49Rciek+R99LK+SytXKC6QKSYAf+2hRB3CyHyAaSULinlrcAM4FqMzLs5Usrrk00zl1IGgKsxhKkN+ApwvZSyWQjxZSHEweh1PmA9xvnTWYxef9+VUj5zvn/ITKMu4fypII2RKC6EkKvDyNaLF6cVSpyylWA4yJYPd/D2yd0xcbJZbFw181IuqlqkxCkNDLaC+jHwEuCKs3w/gVFge15E2yT1MTOSUn4L+Fbc+ARwTe/rRgtnmtX5U7bTrzgtv1qJU5bS4evi1WPbaPd2xuYm5I5jw4x1FNrV72i6GEyg6jE8oHZhJEn8ixCid0IDAFLKB0ciuNGIPximqT36oaZpVJaqvexsI+zuNFLJ/UZSqma2GOI0vjLNkSnOhw/bT/PmhzsIRp2NAUTJdC6ZvFxl6KWZwQTqTuB+jJWMDmwA+nPQ1QElUElS3+yKdTkuLXaQYxtV+R+jnrC7k46dzxHxGUkumtlC4UV/pcQpC4noEXaf2ce+sz15WiaTiTWTVjC7dHoaI1N0M+Cno5TyLaIFuEKID4ENUsqWga5XJIdKL89ewp6uRHEymSm86GpsE5Q4ZRueoJc3TrxNfVdjbC7fnseG6WspzVONfDOFZJvFTgUQQlijj9F63d+7AFcxAPHnT6q9UfYQ9nQZZ059xElZdWcbDa5mXjv+Fp5Az8dWTVElH5u2WlliZBjJ9uJbCfwcmN/rLg1ji09t1CaB1x+ipcM4fzJpGhNL1PlTNhD2OOnc+Txhr/HlIiZOJdVpjkwxFHRd52CTZMfp93vMBDWNiyoXsGTifNWqKAMZSi++TuB6IKm0ckVf4ldP5eNzlYNuFmCI03OEvUbXec1kpnDZVUqcsoxgOMi2j3ZyPK6Xnt1i57Jpq6kpUlu0mUqyArUAuFhK+cFIBjPaOaPOn7KKsLc/cboSW2lNmiNTDIUObyebj2+nIy6FvCRvPBumr6VApZBnNMkK1GGgEqPLuOI8SUiQUOdPGU3Y66Lz3V7itPRKbKWT0hyZYiicaDvFmx/tIBTuSUCeXTqD1ZMuUinkWUCyAvUI8KQQ4hHgKL0axEopNw13YKMNlzdIu9Pocm02qfOnTEbXdZx7X+srTmVKnLKFiB5hV91e9jf0OPyYTWbWTF6OKFEp5NlCsgL16+h/v9PPfSpJIgnONDljtyeW5GExq7YpmYrv9GGC7dEmoZpGwZINSpyyCE/Qy2vH36LB2RSbK7Dns2HGWkpylWt1NpFsmrn6NL1AEtobqe29jCXi9+CWPRZnudMWYy+fkr6AFEOiwdnEq8e34w36YnOTiqv42NTV2C3KFDTbGFCghBC53fVNQojepoIJqDqoc6MaxGYHrsPvoEfdcM2OAnJnLEtzRIpk0HWdDxpr2Vm3JyGFfHnVQhZXzFMp5FnKYCsopxBiopSyCcNUUO/nGlUHlQSdLj9dbuPYzmoxUTaut2+jIhMItNThrz8WG+fPX4dmVq2oMp1ANIU83o49x2LnsumXUF04MY2RKS6UwX77LsOwxQD4WApiGbXEb+9NLMnDrM6fMg49HMJ1YFtsbK+codLJs4B2byevHt9Gh7enPLM0bwIbZqxVpoKjgMF68W3t77Zi6MTXP1WXqu29TMRz7H3CHuNDTrPayZ+zOs0RKc7F8baTbP3o3YQU8rlls1hVsxSzSiEfFaj9ixFG1/XE/nuqQDfjCDnb8J7YGxvniZWY7IMeuyrSSCQS4d26PRxorI3NmU1m1k1ZycwJU9MYmWK4UQI1wnS4/Li8hs+M3WqmpFidP2USuq7jOrANPeqgah1XQU7NnDRHpRgId8DDa8ffotHVHJsrzClgw/S1TMgdl8bIFCOBEqgRJn57r7I0H5NJZRNlEv46mVDzlD9/ncr4ylDqnY28fvythBTyKeOquXTKKmwqhXxUogRqhElIL1f1TxlFxO/FVbsjNs6duhhLgSrkzDR0XWd/42F21e1NSCFfUbWYRRVz1BeKUUyydhulwFeBZYCVvn5QK4Y/tOynz/lTuRKoTKJPzdPMpWmOSNGbQDjI1g/f5cP2U7G5HGsOl0+7hKrCijRGpkgFya6gfgmsBDai7DaSprXTh9dvZBg57BbGF+akOSJFN0bN09HYOH/+WjSzNY0RKXrjC/p44cgbtHnaY3Nl+SVsmL6WPJtKYhkLJCtQ64BPqnTzoXGmV/dytRWRGRg1T9tjY/vEGapLeYbRnzjNK5/FxdUqhXwskaxAtWMYFiqGQJ1KL89IPMf3EPYYb2fNYiN/rqp5yiR8IX+iOGka66esVF3IxyDJCtT9wKNCiC/Rv92G6sXXi0gk8fxJGRRmBiFnO97je2LjvNmq5imT8IX8vCBfTxCnS6dczKySaekNTJEWkhWoHwLFwLsD3K/W3L1o7vASCIYByHdYKc63pzkiha7ruA7G1TwVl5NTMzfNUSm68YX8vKjESRFHsgL1qRGNYhSizp8yD3+dJNh21hiomqeMolucWpU4ZRx6OEzHvv34W1oYt2wZ9gmpK8VI1g9qK4AQwgHMBEzAcSmlc9AHjmHqmnv+apS9RvqJ+L24E2qeFmEpnJDGiBTd9CdO66esVOKUAQQ7O2nY/Br+ZqNzhx6OMPHqK1P2+snWQZmBbwH/SE8dVEAI8Rvg76WUoUEePuYIhyOcbXbHxur8Kf24a3cQSah5Uj5PmYAv5GfTkTf6iJNKiEg/zqPHaH5zK5FgMDbnmJja2rNkfR++CXwGuAOoif7cAfwVRgKFIo7Gdg/BsHHOUZhnozBPtWFJJ4GWOnxnjsTG+fNUzVMm4A8F2HTkDVrcUVcfTWPdZCVO6SYSCtH05jYaX30tJk6ayUTJJaspWrQwpbEkewZ1B/A5KeWmuLn/FUI4gZ+jRCqBBHsNtXpKK31rnqZjK1M1T+nGHwrw4pHXe8QJWDd5JbNLlTilk0B7Ow2vvEqgreffxVpQQPmVG8gpK0t5PMkKVD5wrJ/5E0DJ8IUzOlD27plDYs2TVfk8ZQD9itOUi5U4pZmuWknLtu1EQj0nNvkzplO6fh1me3qykJMVqN3AFzHOoOL5e+Avyb6YEGIR8ASwEEPcPiul3D3I9cXAfuB+KeVvkn2ddBIKR2hojTt/Ug1i00bI1d7L5+liTDnKZTWd9C9OauWUTiLBIM1bt+M80rMNrpnNlKxZTeHcuWnNdE1WoP4v8KYQ4lJ6aqEuBqYAVyXzBEIIG/As8GOM1kk3ApuFEJOllAP193sCqEoyxozgbIubcMTouDyuIIc8hzrrSAeGz9N29IhRi2YtLiNnkqp5Sid9zpyAtVNWMLt0RhqjGtv4W1ppeGUzwc6eRkG24mLKP74Be0n6s1yTSpKQUr4HLAVexUiQKAGeB2ZLKXcM9tg4LgWsUsofSymDUsqngYPAp/u7WAjx/wCFwAdJPn9GoNxzMwP/GUmwrd4YaBr589ermqc0EoiKU7O7NTa3dsoK5pTOTGNUYxdd1+k8eJC6P/4pQZwKxCyqb7oxI8QJhuAHJaU8Atx3Aa81Fzjca64WWND7QiHEVOABYDXw8gW8ZsqJP39S6eXpIRLw4j7c873JMXWhqnlKI4FQgBd7idOayUqc0kXY76f5za24jp+IzZksVkrXr6VAzEpjZH0ZUKCEELuAK6WU7UKI3YA+0LVJ+kHlA7179nmAhEZo0ZqrjcB9UsoGIUQST50ZBIJhmtp6/ojKoDA9uGvfTah5yptxUZojGrsEQgE2Hd3SR5zmlilxSge+xiYaN79K0NnTSMA+YQLlH78C27hxaYysfwZbQb0I+KO3XxiG13IDjl5zuYCr19x/AFJK+adheM2UcrbFTSTq+FlS7CDHrgyLU02gtR5fnYyN8+etRbOoc8B0EAgH2XR0C02ultjcmsnLlTilAV3X6dy3n9Z3d6JHIrH5onnzmHDJKkyWzPysGjAqKeXX44ZbgB1SymD8NUIIO0axbjIcAr7Ua2428N+95m4BKoUQN0THBcDjQogVUsovJPlaaUHZa6QXPRLGdaDHssxeMU3VPKWJQDjIpiNvJIjTJZOXM7css7aQxgJhn4+m17fgPnkyNmeyWim77GPkT8/sdlLJyuYWoAJo7jU/DfgdfVdGAz2HFrXseBQji28h8Of4i6SUs+PHQoi9wI+zIc28rqln2azSy1OP5/gewu6emqe8uZekOaKxyUDiNE+JU8rx1tfT+OrrhNw9pS85ZWWUf/wKrIWFaYwsOQY7g/o7oHsVpQGHhBC9z6HygT0kgZQyIIS4GiN1/EHgI+B6KWWzEOLLwG1SynlDjD9j8PlDtHT4ANA0jUolUCkl5OpI9HkSKzGrmqeUEwgHeenIFiVOaUaPRGh/fw9tu94jPn2geNEiJly8As2cHQ5Jg62gnsQ4NzIBvwK+QaKrro5xfvR6si8mpTwArOln/lsYzWj7e8ziZJ8/ndS3uNGj509l4xzYrdnxBhgNxHyeojVPliJV85QOusWp0dWz0bJ60kVKnFJMyOOh6bU38NTVxebMdjtll3+MvClT0hfYeTDYGVSI6PmQEOJD4B2gQErZFp1bDuxRncwN4rf31PlTavHXHyXY2lPzVLBgHZqWbB9kxXAQCAd5+WhfcZpfnj1ZuKMBz+k6Gl97nbDXG5tzTJxI+YbLseRn3+dSsr/FzYAE/j1u7nngAyGE6lFC7waxqv9eqogEfLgPvxMbO6YsxFKo2kOmkmBUnBqc8eK0TIlTCtEjEVrf3Un98y/GiZPGuGVLqfzktVkpTpB8ksQjwHZ6zqQApgL/iZHwcPUwx5VVeHxBWruM8yeTSaNigjr7SBXu2neJBIy/e7Mjn7yZquYplQTDQV46+maCOK2atIz55bMHeZRiOAm5XDRsfg1fQ0NszuxwUH7F5eTWVKcxsgsnWYFaAdwtpYwtE6SUXiHEQ8D7IxJZFhHf3qhifB5Wi9peSgVGzVNtbJw3d42qeUohPeLUFJtbNWkZC5Q4pQz3Rx/R9PoWwn5/bC63upqyKy7Dkps7yCOzg2QFqg2YDxzvNT8LGPO273XK/ynlGDVP22Jje8U07OVT0hfQGKM/cbq4ZqkSpxShh8O0vruTjn37Y3OapjF+xXKKly4ZNX0nkxWoXwA/F0LUAN15i0sxuj78coRiyxqUQWHq8R7fS9jdAaiap1RjnDlt7SNOCyvmpDGqsUOws5OGza/hb+7ZVrXk5VH+8StwTJyYxsiGn2QF6pvRa/8DKI3ONQE/Ar43AnFlDS5PgA6Xsby2mE2Uj8/+ZXWmE3Z34jnes7OcN2uFqnlKEd3idNbZGJtT4pQ6nEeP0fzm1pgVO0De5MmUXf4xzDk5aYxsZEhKoKSUEYzu4g8IIUqAwCAeTmOK+PZGE0vyMJvV+dNIous6zgPxNU+l5EzO2vrurCIUDvHKsURxWlmzRIlTCoiEQrS89Q5dhw7F5jSTiQmrLqZo4YJRs6XXm6Q7BAohlmBYZpijYw2wA8uklPeMTHiZT11jnL2G6h4x4hg1T2eiI42CBetVzVMKCIVDvHzsTeq7EsVpUYUqiB5pAh0dNLy8mUBbj9GjtaCA8is3kFNWlsbIRp6kBEoI8RWMThIuIA+jo0RR9O5NIxNa5qPruirQTSF9ap6mLlA1TynAEKetCeK0olqJUypwf/QRja++nrCllz9jOqXr12G229MYWWpI9qvnPcC/SCkLgbMYTV6rMOzfd49QbBlPlzuAy2u8cWxWM2Xj1PnTSBJf82TKySNv5vI0RzT66RGnnhqbFdVLWDxRidNIous6be/9hbObXo6Jk2Y2U7p+LeUbrhgT4gTJC1QF8Ez09l5glZSyAfhX4PaRCCwbiE8vryzJw2QanfvAmUCwLbHmSfk8jTyhSJhX+ojTYiVOI0wkEKDh5Vdo29Xz3d9aUED1DX9N0bx5o/a8qT+G0uqo2zP7CLAoevsMUDncQWULqv4pNeiRMM4D22Nje/lUVfM0woQjYV49to0zceK0vHoRiyeqhJSRJNDRSd0zf8b94UexOUdVJdWfugF76djbzk42SeJZjDqouzB8nR4TQryK4el0ctBHjlJ0XU/oIFFVqvrvjRTeE3sJu9oB0Myq5mmkiegRtny4g9Od9bG5i6oWsWTi/DRGNfpxnzxF46uvEQkEYnPFCxcwYfUqNNPYTARKVqDuA36I0U3iKQxheh0jaeJvRia0zKbd6cfjM/aG7TYzJcWjrwYhEwi7O/Eci695Wo7ZoVarI4Wu67x1cjcn2nq+dy6pnM/SSiVOI4Wu63S8v4fWnbvp9m7SzGbKLl1PgRjbViXJCtSNwFeklK3R8Z1CiL8HfGPVbiMhe680f0ztC6cKw+dpe2LN0xT1QTlS6LrOzro91DYfi83NK5/FRZUL0xjV6CYSDNL0+hZcJ07E5ix5eVRcfeWoTyFPhmQF6qfALqBboIhvHDsWUfYaI4+//hiBlm7TNY38+WtVzdMIsrfhIPsbDsfGMydMZXXNRerL1wgR7Ozk7EuvJNQ3OSorKf/4Biy5jjRGljkk+9u+E/jrkQwkm9B1PaGDRJVKkBh2+vo8zcdapL5RjhQHm46wu25fbDxlXDXrp16sxGmE8Jw6xek//ClBnIrmz6fy2muUOMWR7AoqAnxLCPFV4EPAG3+nlHLFcAeWybR0+PAHjG2n3Bwr4wrGRk1CKnHLnUQCxtvMlJNH7ixV8zRSHG39kLdP9qQ0VxZWcNm0NZjUanXY0XWdjj17aX13F7HzJpNR31Q4R3WC702yArUz+qMg8fypSp0/DTuBljp8p3u2mvLnrsFksaUxotHLR+11vPnhjti4LL+EK2esw2IypzGq0UkkGKRpy5u4jvW4Flny8qi46kpyytXuQH8MKFBCCFO0SSxSyq8PdN1YRNU/jRwRnxvn3tdjY1vZFOwVU9MY0ejlTFcDr53Yjq4b3+TH5xZz1cxLsZpVAfRwE+zqouGlV/C3xo7xyamooOKqj48KY8GRYrA1fFAIkSDrQoh1QogxvZ8VjujUtyiBGgl0PULX3td7tvbsDgoWrEtzVKOTJlcLrxzbSiQSAaAwJ5+/mnUZOZYx/es9InhO11H3h2cSxKlo/jyqPnmtEqdzMNgWX3/7Vi8Ai4ET/dw3Jmhu9xAMGb/UBbk2CvPU1tNw4TnyHsG27uJQjYLFV2Cyq1/g4abN08Gmo1sIhY0KkTxbLtfMupxcqzqcH050Xadj7z5ad+xEnTedH0nbbUQZ84ctvbf31PnT8BBoPoXn+J7YOG/mRdgmVKUxotFJp6+LF4+8TiBkdCvIsdi5ZtZlFNjVTsBwYpw3bcV1rKemzJKXR8WVG8ipqEhjZNnFUAVqzBMvUCq9fHgIe104971B97dMW0k1jhlL0hvUKMQVcPPikTfwBo2O8Fazlb+adRnFjqJzPFIxFIJdThpefgV/S0tsLqeigoorN2DJU87PQ0EJ1BAIhyM0tLpj42plUHjB6HoE597XEmw0ChZdrgpyhxlv0MemI2/g8hvvX7PJzFUzL6Ukb3yaIxtdeOrqaHzlVcJ+f2yucO5cStdegmZWmZFD5VwCdacQIr5jhAX4jBCiJf4iKeXjwx5ZBtLQ5iEUNs6fivPt5Oeq86cLxSN3EWzv7pitUbj4ckx2dRYynARCAV46uoUObxcAJs3Ex2esY2KBSm0eLnRdp3Pfflp3vBvLitRMJkrWrqFonrInOV8GE6hTwN/1mmsA/rbXnA6MCYE6o9LLh5VA00k8J/bGxnmzlmMdP2bdW0aEbsPBFne0Y4Gmcdm01dQUqb/n4SISCtH85lacR47G5swOBxVXXYljojpvuhAGFCgp5ZQUxpEVJBToKoG6IMJeZ/TcycBWOgnHdHXuNJyEI2FePb6dBmdTbG7d5JVMGz85jVGNLoJOp1HfFH/eVF5OxZUfx5KvzpsuFHUGlSTBUJiGNk9sXKXOn84bPRLGuec1IkFjn944d/qYyogcRvrzdLq4ZimzS6enMarRhfdMPQ2vbCbs88XmCufMpnTdWnXeNEykVKCEEIuAJ4CFGLVUn5VS7u7nug3At4GZQBPwPSnlz1IZa2/OtriJRIy95QmFOeTmqGr788UtdxHsaDQGmkbhkg2YbOrcabjoz9NpaeV8FlbMSWNUowdd1+n84ACtb7/Tc96kaZSsW0Ph3Lnqi9YwkrJUKSGEDcOZ9/dAMfBNYLMQorDXdTXAM8BD0etuBR4WQlyZqlj7Q6WXDw/+xo/wftjTNTtPrMQ6Tu3TDxcDeTotU55Ow0IkFKLpjS20vPV2TJzMDgeV119H0bx5SpyGmVSuoC4FrFLKH0fHT0dNDz8NPBl33RTgd1LKP0fHu4UQbwKXAK+kJtS+xNu7K/+n8yPsceLcvyU2tpVNxjF1URojGn3sOas8nUaKkMvF2Zdewd/cHJvLKSsz+unlqy+tI0EqBWoucLjXXC2wIH5CSrkd2N49FkKMB9YCvx3pAAfCHwzT1G70h9M0jcpSdfg5VPRImK49r6JHz53MjnwKFqpzp+HkQKPkvTPK02kk8NbX0/DKq4S9PU5DhbMFJevWYrKoo/yRIpV/s/mAp9ecBxiw2ZoQogh4DsPq49mRC21w6ptdseV8abGDHJt6Qw4Vd+27hDqj2WRatM+eLSe9QY0ijrSc4J1T78XGVYUVXK48nS6YAc+b1lxC4Xy1pTfSpPKT1g30PgnPBfq1jhdCzMIQpUPAbd3WH+lAnT9dGP6GE3g/+iA2zpt9sTp3GkY+aj/N1o/ejY3L8kv4+GTNA9gAAB93SURBVIx1mJWn0wURaGunedt2vPU9mZBmh4OKKzfgqFR1ZKkglQJ1CPhSr7nZwH/3vlAIsQ5DnJ4Aviyl1Ec+vIFJOH9S6eVDIuzpwrn/zdjYXj4FxxR1YD9cGJ5ObyV4Ol0982PK0+kCiIRCtP/lfTr27EWP9HwvtpeWMvHqK9V5UwpJpUBtATQhxJeAR4EbMdLN/xx/kRBiOoatx1eklI+kML5+8fpDtHRE/Yk0jYkl6vwpWWLnTtHO2WZHAfnq3GnYaHQ19/J0KuCvZl2GXbkPnzeeU6do3rqdoLOnKF/TNIoWLWT8iuXqvCnFpOxvW0oZEEJcjbEqehD4CLheStkshPgyxjbePOCLQAFGavnDcU/xmJTy/6Yq3m7iV0/l43OxWdW2SbK4D79DqNPIeNI0EwVLNmCyKkO84aDV085LR9/s5el0mfJ0Ok9Cbjctb7+TYMcORleI0vXrsJdMSFNkY5uUfh2QUh4A1vQz/y3gW9Hb/wz8cyrjGgx1/nR++OuP4T15MDbOm7MKa7FqTjocdPq62HTkDeXpNAzokQidBw7StnMXkWAwNm+y2ZmwaiWFc+eoFX8aUevVcxDfIFa1N0qOsLsT54GtsbG9Yho5k+enMaLRg/J0Gj78zc00vbktoa4JoGDWLCasXoUlV61G040SqEFweYO0O40PArNJnT8lgx4ORc+djG+j5txC8hesV99ChwFv0MeLMtHT6Wrl6TRkIoEArTt30/nBAbpNMgGsRUWUrl9LbnV1+oJTJKAEahDOxHUvr5iQh8WsakrOhevwO4S6jM7Omsls9NlT504XTCAUYNORN+j0JXo6VShPp6TRdR33iQ9peettQu4e41HNZGbcRUsZt2SxavKaYSiBGoTE9kZqe+9c+M4cxXfqUGycN2c1lqLSNEY0OjA8nd6k1dNuTChPpyET7OqiedtbeE6dSpjPra6mZN1abMVqizQTUQI1CHVNqv9esoRc7bgObIuN7ROnkzNJOYleKOFImM3Ht9Hg7DknUZ5OyaOHw3Ts3Ufbe39BD4dj82aHg5I1l5A/Y7rafs5glEANQKfLT5fbyJKyWkyUjVMHpgOhh4M4338VPdx97lSkzp2GAcPT6R3qOs/G5pSnU/J46+tp3rqdQHt73KxG0fy5jF+5ArNdbT1nOkqgBiB+e29iSR5mdf40IK6DbxNyGZbimslM4dKPY1LFoheEJ+hlx6m/cKKtZ0tqaeUC5emUBGGvl9Yd79JVKxPm7SUllK5fS055eZoiUwwVJVADEJ9eXl2qtvcGwnfmCL662tg4b+4lWApVUeP54g8F2N94mA8aa2NFuADzywXLKhcM8kiFrus4ayWt7+wg7PfH5k0WK+NXLqdowXw0k/qimU0ogeoHXddVgkQShJy9zp0qZ5JTo77hnw+hcIgDTUfY23AwVoDbjSiZzqqaZWrLdBD8rW20bNuO9+zZhPm8qVMpXXuJ6p+XpSiB6ocOlx+X1zhPsVvNlBSr86fe6KEgXXs2o0e/5ZvziimYv059iA6RSCRCbctx3j97AE8g0Y1mfO44llctZFJRlfp7HYBIMNjT2FWPq2kqKKBk7SXkTZmSvuAUF4wSqH6Iz96rLM3HZFIfDvHouo7r4HbCLuPwufvcSbOoDtrJous6x9tO8l79frp8zoT7Cuz5LK9axPTxk5UwDYL75ClatvVt7Fq8eBHjLlqGyarej9mOEqh+SDx/UlsDvfHXSXxnjsTG+fPXYilQ3QySQdd1TnfWs+vMPto87Qn35docLJ24gNkl0zGps5IBCbnctLz1Nq4TJxLmcyoqjMauE9R7cbSgBKoXfc6fypVAxRPqasV1cHtsnFM1i5zq2WmMKHtocDax68zehJomAJvFxuKKecwvm4XFrH4lB0KPROj84ABtO3cTCfU0djXb7UxYvYqC2UKtOEcZ6rehF62dPrx+41zFYbcwvlDZkncTCQWMc6eIUfBozh9H/ry1aY4q82nxtLG7bh+nO+sT5i0mC/PLBYsq5ioPp3Pga2qi+c1t+FtaEuYLxCwmrFKNXUcrSqB60bt7ufpGZqDrOq4D2wi7OwHQzBYKl6hzp8Ho9HXx3pn9HG87mTBv0kzMLp3B0sr5yr/pHITcbtr/sofOAweJb+xqKy6mdP06HFWq3dNoRglUL+pUenm/+E4fxl9/LDbOn78OS8G4NEaUubgDHt6vP0Bty7GEzDI0jRnjp3BR1UIKlXdTv+iRCL6GBjynTuM5dbrPikkzmxm3TDV2HSsogYojEkk8f1IGhQahrhbch96OjXOqZ5NTNSuNEWUmvpCffQ2HONAoCUfCCfdNLq5medUixucWpym6zCXkcuM5fRrPyVN46uqIBAL9Xpf7/7d37nF2VFW+/1adc7rTr7wIQickPBxZmJAwMYA4Ojz8XEVHrzCiFzUfr4iMDuJwZRDHgYuCMMqAMAwyIyJccVAvjuKDy4xzYdTLOOqMARQyBNcYMCSQBzHP7iTd51X3j7Wru3LSnZwOfU6fdK/v53M+p2rXrqpV1X32r9beu9aaP5/DT38dhRke2HWq4AKVIUkSSuUqAN0dBWZ2e6yuarnIzscfHhp3yvfMpnvRPkmRpzSlSomVm5QnN66iWCntta235whOPeokjuj2qO4pSaXCwKZN7H5uHbvXrmVwy5ZR60ZRxLTeXmYsPpGu4471LvcphgtUhlwu5rQTj+TZF3Zw6qIjp/yPIUkS+lc+QmV3Ou5UoGfpG4l8phlgkcaf3ryaxzf8BwMhw23KnK7ZnDLvJI6a3jvl/4+gxktat26v9Oq15Lu66Dx6AZ0LFtB51DziNp9AMlXxlqaGZSccwbITPJgkwMDapxjc8MzQevfiM8h3exdVNamyessaHl3/5FB225QZ06ZzyrwlHDtrwZQWpqRSYWDjJnavXWtjSfvzkuKYab1HmiAtWEDb7FlT+t45w7hAOSNS2vEiu57+2dD6tAULmTb3dybQooknSRKe2/48K154gm17duy1rautk2VzF3P8nOOIo6n5km25fxe7165l13Nr2fP88+4lOS8ZFyhnH6qlQcvvlI47TZ9D9yt/b4KtmljW79zIz194ghf7955V1p5vZ2nvIha+7Hjy8dSaVZb1knY9t5bi1q2j1jUvqZfOBfPdS3LqxgXKIUkSKn1bKG1ZT2nrBkpb11MtWbqCKF9g+tI3TNlxp827tvDz53/JCzs37lWez+VZcsQrWXLkK2nLTZ13wcr9/ex6zrrtDugldXfTuWABXUcvoGPeXPeSnDEzNVudKU6SVCnv/K2J0ZYNlLZtICkNjli3Z/GZ5LqmxrTeJEnYXdrD9oGdbNuzg/V9m1izbd1edeI4ZuHhx7O0dxEdhckfZeRgvKSu0HVXmDXTvSTnJeECNQVIqhUTpNRD2raBpDz6ky9A3DaNzlecTHvv5EsvXk2q9A/uYtvAjiEx2j6wk+17duwzTXyIKEIOO45lcxfT3d7VXIMbSLVcpty/i3JfH+X+fko7+yj391Hu67eyXbtIqtVR9y/09NAxf755SUfN8wjizrjiAjUJSaoVyjs2B0FaT2nbJpLRGt5A3N5BYfZcCrN7KcyeS6770B8jqFQr7Bjo20eIdgzs3OdF2v1x7KwFnDxvCbM6Dj1PsjI4aGLT3xfEJwhPXz+lvj4qe/aM6XhRnKNj7vBYkntJTiNxgZoEJJUy5R0vUtqygeLW9ZS3bRya4DAa8bQuCrN6KRw21wSpa8Yh29AUK6UhD2jbwA6279nJtoEd7Bzsh2yooToo5ArM6pjBzGnTmTltBkfNOJI5na2ZviFJEiq7d1Pu66PU17+P+JT7+vY7RlQvhZ4em2139Hw65rmX5DQPF6hDkKRSorRtU5jQsIHy9k0HFKRcR3fwkMxLijunH1KClCQJe8oDIwjRzn0y0dZDR2FaECITo1SUOgsdLXFfkiQhKZWoDAyY59OXej9BfMLygf7uByYi391FvrubwvQe8t3d5Ht6KPTYd7672wXJmTBcoA4BknKJ0raN1l23dQPl7S+SJKOPCwDkOnqGvKPC7LnkOnuaZO3YqSZVipUSxXKRwYp9iuUi/cXdQ11z2wZ2UCyPHKNtVKKInraujADNYGbHdGZNm9Hw9BZJpUK1WKQyWKRaHKRaLFIdLIbvwVBuy9VicfiT2ZaN3n2wRLncsPgEwSn09JDv6Sbf3UO+u4vIkyM6LYoLVIuRVCskpaKNIQVBKu148YBdVbmuGXuPIXU0N9BtpVoZEpbBSpHBcpFipZRZtu/BSs1yqPdSiKOYGdN6mDltBrM6pg8J0cz26WNOAJgkCUm5TFKpkJTLw0KSEZBUVEYUmXS5XH5J11Qvufb2IeGp9XzyPT3kOqa1hEfoOAeDC9Q4kyQJVCtUS4Mk4VMtDZKUi8PLQ+VFknKmTql4wMkMKbnuWbQNCVIv8bSxzyxLkoRKUqVSrVCulu07qYT1yrCw7FdkzPMpVxvUICcJVBOiapUCMTPyXczMdzG90Mn0fCc9+U664jYoV0iKFaq7yyTlbSTlzeyoVKiWyiSVMkm5TLVc2Ut8quXy0PpQvcpL7TIbP+J8gbi9LeP5BPHJrPu7Rc5kpqkCJSInAXcAS4BngQtVdcXB1msUSZKQVEokxUGSclZUisPL5YyoDIlQWB/DuECSJFSBKglJ+K6SUEmgQjL0oXM60fQ5RNNnQ/csqvl8EJIBKltWDwtMdVhgKkl2vUylWqZcqVKulKhWylQqZaKEYRFIEkgSomr6TfiuDpVH1VQwQp1KlfYkoT3sH1Wrw9uHyhKoVm3/JBlhe5UcMXkickn4JiKfxLRFOdpzbbTl2ijEeYj2jn3XHz6tR0Tc1kauvY24vZ24rY24vS2UhfW2sK22vL2duFDwfEfOlKdpAiUibcD3gFuB04HzgIdE5GhV3TnWeo1i88pH+PXKn1Aul6km1dB2J0FI7DtbliSJLTO8TLodhuokJCRV2zetm26Pwj4RhO+IKMoR5wtE+TbiuEDEdkjWDjfwqZiEdRL2LksYEpFckpAPy40iiiLiKCYX5cjF8V7LuShHHOfIRTG5KM4s54hzVreVuqGiOEeUzxPnc0SFgglHVmhGFZl2E6S2NqJCoaWuyXEORZrpQZ0JFFT11rB+n4h8BDgf+NJB1Bt3ysUBfvydbxD3jRxVoR7S4eY4ioAIohhbjMN6ZA1X+iEeLiMtS6kCA+HTGKIoIiYiCiIRR2GZaEhchgUnJhfnaoTIxCaO4qCwDbAxjonyeaJcjjift+V0vZAnyuWHBWVoW3Y9FXvbx9bzRPnC8HIus68Li+O0BM0UqIXA0zVlvwIWH2S9caeaQKGYY6iDLoqDYNSIShQRkRWZVIhqBaY+UnGKbYUoHD8rFnEoT8VkeN3OO9I22y8jPETEsXkwUZwjzuWI4xji2EQgXY6i4bJcbuibOCbKZcritCxn5aPVj4fL9i5P62eOkz1H+vFZZo4zJWmmQHUDtS+s7AY6D7LeuNPWPo3fPf8C1v1ihY2ZxDmiODTyYTmK0gY+zjT2w+XDZfGQANRuG1rO2bIJgh2DKKoRirAeZUUkyohIzvbJpQKZEZqhfWuO6TiOcwjQTIHaBXTUlHWy7xh3vfUawvzFy5i/eFkzTuU4juPsh2b2nawCpKbshFB+MPUcx3GcSUwzPagfAZGIXAbcjs3OWwJ85yDrOY7jOJOYpnlQqloE3owJzlbgKuBcVd0sIleKyFMHqtcsWx3HcZyJp6kv6qrqfwCvG6H8M8BnDlTPcRzHmTr4/F3HcRynJXGBchzHcVqSyRgsNgewcePGibbDcRzHCWTa5LqDTE5GgeoFWL58+UTb4TiO4+xLL/BMPRUno0CtAH4f2AC0Tu4Ex3GcqU0OE6e6M1NESQMjXDuO4zjOweKTJBzHcZyWxAXKcRzHaUlcoBzHcZyWxAXKcRzHaUlcoBzHcZyWxAXKcRzHaUlcoBzHcZyWxAXKcRzHaUkmYySJuhCRk4A7sGSIzwIXquqobziLyKnAg6r6slayT0TeANwAvAJ4EbhJVb/YYja+FUuncmyw8cZWszFTfybwJPBJVb2nlWwUkQuBLwKDmeJLVPUrLWJfL/AF4CxgALhTVa9upG1jsVFElmP3L0sH8ANVfWMr2BjqnQbchmUW3wzcoKp3Ndq+Mdr4euAm4Hjg18AnVPWh8bZnSnpQItIGfA/4BjAT+AvgIRGZPkLdSEQuAh4C2lrJPhGZD9wPXB/qvRv4rIic3UI29gLfAv5MVXuAdwK3isirWsXGGu4A5jXatpQx2vgq4GZV7c58Gi1OY7Hve1iIsSOA04D3ich7GmnfWGxU1a9l7x0WEm0rcEWr2Cgicah3m6rOwH7TtwfhaBUbjwH+D3APMBv4CPB1EVk03jZNSYECzgQKqnqrqpZU9T7gKeD8EepeC1yMiUCr2XcM8HVV/Y6qVsOTzv8DXtsqNqrqBuBwVf1++PEdBpSBvlaxMUVE3gdMB1Y2wbaUM6nfxmXAL5toG9Rpn4i8GjgOuFRVB1T1N2HfH7WKjTX2FoCvAdeo6hMtZOMs4GVAJCIRkGC/l2IL2fhm4Feq+vlQ76fAN4ELxtugqSpQC4Gna8p+BSweoe4dqroMeLThVg1Tl32q+mNV/eN0XURmY0+Fv2i4hWO4h6raJyKdWNfUQ8DfqOqvG29i/TaKyLHAp4ALm2BXlrpsFJEc1u3yXhFZLyKrReQToRGbcPsw8VwJXCMiL4jIM8AfhgeURjOW33PKJcAe4G8bZVQN9f6mtwC3A18BSlhg1StVtXbfCbMR041dNWUVrLtvXJmqAtUN7K4p2w101lZU1fVNsWhv6rYvRURmAA8A/4656Y1mrDYOAF3AKcCFIvKBBtqWUpeNofH/KvAxVW12IrF67+Ph2EPSV7CxvHdgnv3FLWJf+nBUwjyptwMfa0YXH2P8XwxdWVdg3lOzomXX+78YY7+V92DjY2cCnxKRho+R1Wsj8H3gVSKyXEQKYczsXZi948pUnSSxi31vZifQPwG2jMSY7BOR4zFRWgUsV9VqY80DxmhjsKkIPCoidwLnAHc31ML6bbzaTNRvN9iekajLxiCcZ2SKfikinwfOo7FeQL33cBDYqarXhPUnROQuTKi+3kD7YOy/5zcBVeAfGmlUDfXa+Hbgtaqajos9IiJ3Ax/Ceh8aSb3/i8+KyDnAX2KTOX6CPTgdNd4GTVUPahU2QybLCaG8FajbPhE5HfOavgu8Q1UHGm8eUKeNInKGiDxWU68d2N5A21LqvY/vAt4hIttFZDvWpfG3ItKM7p967+MiEbm2pl4b9rTdSOq9h78COoN3ktKsB+Cx/p7PAf6+SQ9yKfXaOB/7fWQpY55po6n3f7EH2KyqJ6vqYar6NmAuUPs7f8lMVQ/qR9gg5GVYf+95WP/+dybUqmHqsk9EXg48CFylqp9vRRuxQf15IvKnwF8DrwY+APxhq9ioqidk10Xkl8CtTZpmXu993A5cLiLPY57nUuBSbAZVK9j3MDYl+mYRuRxr6D5A47sgx2JjymmY19xM6rXxIWwm7geBL2EzN/8IuKiFbDwM+Pcw1XwFJvhvBC4bb4OmpAelqkVsJsp52DTTq4BzVXWziFwpIk8dIvZdAvRg/9D9mc9ftoqNqroD+AOs62IrcCdwkao+0io2TiRjuI8vAG/Dunp2Yq8XXKeq32oR+wawLsjjsKnm/4S973Z/I+0bi40ZjgGaOrY8hvv4FPZb+RD2UPJ17B2jho8rj8HGNdjDx73ANuBy4A8aMX7rGXUdx3GclmRKelCO4zhO6+MC5TiO47QkLlCO4zhOS+IC5TiO47QkLlCO4zhOS+IC5TiO47QkLlDOhCEiXSJynYj8p4jsEZHfiMjNISfTeJ2je7zi/oW4Yx8ej2M1CxF5VESuaeDxTxaRH4bla0TkoIMqh9Q2fyQi08Zw7h8c7Pmc1scFypkQQo6ZfwP+CxYRYSHwx8DZWA6a8Qo8eTnjF83gPcCnx+lYhzwhyO6d2Aud48Hp4Xh1RbhR1UeBbSFNijMJmaqhjpyJ5wbsAen1qronlP1GRFYBq4H3Mz5BUMczHUWjU1scapwLVFX1Z+N0vIO5v7cBXxaRe5scW89pAi5QTtMRkXZgOZZld092m6quE5GzAA11I+BPsJhzC7CgpFeq6j+G7fdgKQE6sBQUW4C7VPV6EbkAy/GEiCRYmorNWKrqc7EUFhuBL6rq9aFeDrgSi382G4s19pFQ98uZY50F/CeWQvx0LLHcw1gK9hdHuOYLgI9i2UqvwPLnfAn4n2nDGpL+3YLlVloXtn9OVath/0uxwJ1vw8II7ZNEM8TBuwxLvHgbmUZfRPJYAs7lWNbgrcB9wJ9iSfLWAWer6g8y+6wG/kJVv1x7LuzvUhsBPhaRv8byavUDn1XV2zLHezPwWSxW37Ph+r4csrSmyQ37ROT9qnpPiAt3MRaeqB+LQH6xqqYRtv8VSxPxX2lOmhmniXgXnzMRHIc1oCtG2qiqPw2J28DE4lrgk1jgyu8CD8jeKbAvAp4HTgbuAq4TkWWYGNwMPAH0Yg3wLcBrMIESrBFP6xPOcykmJkvDPv8A/CyUbQ3H+inwBUxoTsHi0B0TzjcarwTeCrwe8xA/SAhaKiIvA/4vFsNucbDhEuDjmf2XYpmIXwX8Xe3Bg4hdC3wMC4j68lA35QrgvcD7gFeE9UuweGsbgB9iKcbT452GCdk+aUhCF+3rgr1ZlmL359XBjhtEZHnYZxEWQ/AO4ESsu/RmEXkXdp/PC8d4OfANEXk3cA0moK/AMraeg8WpA4bSuDyMxZBzJhnuQTkTwazwvWN/lYL39FHsCf6+UHxN8DQ+jnkCAM+oahqd+vrw1L1MVR8TkX6gnAayFJGfYB7T46H+TSLySWChiDyOPa1fn+aGEpFLsEayJ9ibZI51DJa9eI2qFkOD2rOfS8pj+brWYPmcrgc+ISKfxoRihapeF+r+WkT+HBPQGzLHuHY/QTkvBr6Q3qswOSSb6O4p4IJMoN41InIFNv53P5a08VYRuVhVS9iY24Mh4G8tS7EH3Np0EduA96vqLmCViJwMfBhLr/5x4Guqekeo+0yIyP8xVb1PRLaG8hdVdY+IrA/2PhjKnxORR4K9WVZhAVadSYYLlDMR/DZ8z9pvLet2moN5L1n+FXhnZr02fXwfUBjlmF8F3ioi78VSVP8u1kWUC+c6nIxnp6p92EQLRGpT5XA15qW9Pcwm+zbwv/dzPeuCOKX8HDginHcRcGYQ1JQY6BCRw8L67gNEjD4RuDVj++5sJG9VfUBEzhKRG8O1L8G8vlyo8m3MKzxbRL4P/Dds4spIHAHsGiH/2MogTikrMG+RcI2Lg5Cn5Bkl15GqPhJm6l2P5SVahHm9td7jFux/xZlkeBefMxGsxhqVU0baKCK3hCf7PSNtx8ZVsv+7xVHqjMT/wiZf7MEautcwnDwxPU5dIf5V9QEswdxlWFffHcA/7meXcs16eg0VrKG+HxPM9LME69pKPZjR7kdKwr7XPXRvwnTzb2OCdD/WLfZE5nr6sS7U87FuyPb9XE91hHOl15IlztiQBz7P3td4IuaN7UPosvwxNhb4fSyx5AMjVM2NcF5nEuAC5TSdMG5wL/Ante+8hG6zDwGDqroTy9vzmppD/B42WaIehsRGLBPoe4H/rqpXquo3sMZzBhCFrqwXyYzbiEiHiGwK3YrZY0UichMwT1XvVtV3YEkY3xDGk0ZigYjMyayfinlVW4GngRNUdXX6wcasPomJQT2sxMZ+UhvbMAFIuQS4XFUvV9V7gTXA0ewtNF/F8nedC3wz5AgaiY1A1wivAywME01STmO4G/Bp4HdqrvFMrAsQ9n0wuAS4SVU/rKp3A09igl0rjHOCPc4kw7v4nIniOuAtwI9E5FNYN91JwI1YQ3RnqHcD8GmxTLKPYU/3Z2MNWz30A0eKyHHYQPwurEvuN1ia6puwBi9Ns/1XwNUisgabpXcV5sH8AjgK6BaRhdgMtIXA7SJyKdat+G6s0U+7MGspAPeIyJ9hEwGuAj4Ttv0NcKmI3BaW52MzBL8XZvHVc623AH8XxtL+DZsEcXhm+xbgLSLyL9gklWuwbtZsivGHsC63C4E37edcT4Z6J4VzpRwO3B26EU/FJrCcG7Z9DsvEeiXw92HfvwLS2Yhp9+YyEXks2HtWuN8R5qkuxIQ4y0mMMuHGObRxD8qZEILX8DrgUawhXoXNgPsu8KbM2MbtmGjdiDVM5wBvVdUf13mqb2HisQrrUnoP9nLwKqyL7yGsGyudxXcTllL9LkyUeoG3BE/iB6HsF5i4fgB7cv9nrMGeH+qO5vFswcad0hmAN2Gigqo+jwnvyVi3273Y+NZH67xOVPWbwP/AZvI9jnleP8xUuQCbar8Su89rwrUuyxyjgonHb4F/2c+5dmJjgWfUbHoYE64V2EPIpar6T2Gfx7BXAc7HJmzcgj2A3Bj2XQk8iP1NPhiuJcH+R/4ZE9LPsreHGwOvDfs5kwzPqOs4TSCMp3xOVeccqO5EIyLfBFar6p8foN47gatVdUlzLBvRhrOxhxjxF3UnH+5BOY4DgIicISIfxbzDL9Wxy/1AJCK/31jL9suHgc+4OE1OXKAcx0l5F/by7MdV9dkDVQ6icBHWldd0ROQUbILLPRNxfqfxeBef4ziO05K4B+U4juO0JC5QjuM4TkviAuU4juO0JC5QjuM4TkviAuU4juO0JP8fbDzqYlujjbMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for gamma in gamma_array:\n",
    "    label = 'gamma = ' + str(gamma)\n",
    "    plot(frame[gamma], label=label)\n",
    "    \n",
    "decorate(xlabel='Contacts per day (beta)',\n",
    "         ylabel='Fraction infected',\n",
    "         loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's often useful to separate the code that generates results from the code that plots the results, so we can run the simulations once, save the results, and then use them for different analysis, visualization, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Contact number"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running `sweep_parameters`, we have a `SweepFrame` with one row for each value of `beta` and one column for each value of `gamma`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(11, 4)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "frame.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following loop shows how we can loop through the columns and rows of the `SweepFrame`.  With 11 rows and 4 columns, there are 44 elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1 0.1 0.0846929424381071\n",
      "0.18 0.1 0.7086227853695759\n",
      "0.26 0.1 0.9007802517781114\n",
      "0.33999999999999997 0.1 0.9568878995442757\n",
      "0.42000000000000004 0.1 0.9770452570735504\n",
      "0.5 0.1 0.9845958628261559\n",
      "0.58 0.1 0.9874003453175401\n",
      "0.66 0.1 0.9884042490643622\n",
      "0.74 0.1 0.9887434214062726\n",
      "0.82 0.1 0.9888495150524135\n",
      "0.9 0.1 0.9888795705171926\n",
      "0.1 0.3 0.0054435591223862545\n",
      "0.18 0.3 0.015914069144794984\n",
      "0.26 0.3 0.055379762106819386\n",
      "0.33999999999999997 0.3 0.2678641677332422\n",
      "0.42000000000000004 0.3 0.5245629358439001\n",
      "0.5 0.3 0.6860504839161878\n",
      "0.58 0.3 0.7883785563390235\n",
      "0.66 0.3 0.8550657464101674\n",
      "0.74 0.3 0.8994791356903035\n",
      "0.82 0.3 0.9294693026191699\n",
      "0.9 0.3 0.9498533103273188\n",
      "0.1 0.5 0.0027357655411521797\n",
      "0.18 0.5 0.006118341358324897\n",
      "0.26 0.5 0.011639469321666152\n",
      "0.33999999999999997 0.5 0.022114766524234164\n",
      "0.42000000000000004 0.5 0.04781622666891572\n",
      "0.5 0.5 0.13243803845818214\n",
      "0.58 0.5 0.30326419264834004\n",
      "0.66 0.5 0.4641102273186152\n",
      "0.74 0.5 0.5884769725281787\n",
      "0.82 0.5 0.6827496109784223\n",
      "0.9 0.5 0.7545952983288148\n",
      "0.1 0.7 0.001826769346999102\n",
      "0.18 0.7 0.003782561608418833\n",
      "0.26 0.7 0.0064266722107564345\n",
      "0.33999999999999997 0.7 0.010190551933453973\n",
      "0.42000000000000004 0.7 0.015945826561526877\n",
      "0.5 0.7 0.025707925046422053\n",
      "0.58 0.7 0.04500775311679983\n",
      "0.66 0.7 0.09069406882939202\n",
      "0.74 0.7 0.18979521165624091\n",
      "0.82 0.7 0.3183431867354656\n",
      "0.9 0.7 0.4369993744563133\n"
     ]
    }
   ],
   "source": [
    "for gamma in frame.columns:\n",
    "    series = frame[gamma]\n",
    "    for beta in series.index:\n",
    "        frac_infected = series[beta]\n",
    "        print(beta, gamma, frac_infected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can wrap that loop in a function and plot the results.  For each element of the `SweepFrame`, we have `beta`, `gamma`, and `frac_infected`, and we plot `beta/gamma` on the x-axis and `frac_infected` on the y-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sweep_frame(frame):\n",
    "    \"\"\"Plot the values from a SweepFrame.\n",
    "    \n",
    "    For each (beta, gamma), compute the contact number,\n",
    "    beta/gamma\n",
    "    \n",
    "    frame: SweepFrame with one row per beta, one column per gamma\n",
    "    \"\"\"\n",
    "    for gamma in frame.columns:\n",
    "        series = frame[gamma]\n",
    "        for beta in series.index:\n",
    "            frac_infected = series[beta]\n",
    "            plot(beta/gamma, frac_infected, 'ro')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap06-fig03.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3XmcXFWZ//FPdXUnnV7Ckg1CgiwmD0aBccedhFGHGQcQ+A0QFBhcEMEFBx0HBEQWFxRR0EGWEbdMmN84iCKjYBIQfrjEFYbAk7A5CQlJkwDp7kovVV2/P86tUF10Nzedqro3Xd/369Uvuk7duveprnCfOuc+95xMsVhEREQkbZqSDkBERGQkSlAiIpJKSlAiIpJKSlAiIpJKSlAiIpJKSlAiIpJKSlAiIpJKSlAiIpJKSlAiIpJKSlAiIpJKzUkHUG1mNhl4LbABKCQcjoiIBFlgb2Clu/fHecGES1CE5HRP0kGIiMiI3gLcG2fDiZigNgD84Ac/YK+99ko6FhERAZ566ilOPvlkiM7RcUzEBFUA2GuvvZgzZ07SsYiIyHCxL72oSEJERFIpkR6Umb0OuM3dZ47y/L7AjcBhwCbgI+5+ex1DFBGRhNW1B2VmGTN7P3AHMGmMTZcC9wPTgA8AS83sgDqEKCIiKVHvIb6LgTOBS0fbwMzmA68BLnT3AXdfDvwYeF99QhQRkTSo9xDfte5+oZkdPsY2C4D/dffesraHgdfVNDIRkZi2rl5D1/IV9G3cROusmcxYtJCp8+dNyOMm9V6hzgnK3dfH2KwDyFW05YC26kckMjHpBFq7Y25dvYa1S5aSbW9j0vRpDHZ3s3bJUuYuPnHCHTep91qSxiq+XmBKRVsb0JNALCK7nNJJZbC7e9hJZevqNRPuuEkcs2v5CrLtbTS3t5PJZGhubyfb3kbX8hU1O2ZSx03qvZakMUGtAvY1s/IkdVDULrLL2bp6DY9eex0PXnwpj157Xc0ThU6gtT1m38ZNZNuGD+hk29ro27ipZsdM6rhJvdeS1CUod3fgz8BlZjbZzBYCRwNLko1MZMcl8Q1fJ9DaHrN11kwKueFXIQq5HK2zRrxrZpc+blLvtSQVCcrMTjaz8iG844CXEe6BugF4n7v/TyLBieyEJL7h6wRa22POWLSQQm+OfG8vxWKRfG8vhd4cMxYtrNkxkzpuUu+1JJEE5e53ufvuZY9/4O4dZY/XuvuR7r6bux/o7v+RRJwiOyuJb/g6gdb2mFPnz2Pu4hNp6exk4OnNtHR21qVoIInjJvVeSybiXHwiI0qiwqx11kwGu7tpbm/f3lbrb/ilk0r5e5199FF1O4HW87hJvtd6naSTPm5S7xWUoGQXNJ5Ek1S57IxFC1m7ZCkQek6FXI5Cb47ZRx9Vs2OCTqAyMaTiGpRIXOMtOkiqsi3pIRKRXZl6ULJLKU80wPb/di1fMeZJv2/jJiZNnzasrV7lsvqGLzI+6kHJLmW8RQdJl8uKyI5TgpJdyngTTdLlsiKy4zTEJ6kRp/hhvEUHSVV7icj4KUFJKsStstuZRKNrQSK7FiUoSYUdKX5QohFpDLoGJamQ9KSUIpI+SlCSCqqyE5FKGuKTRFQWRLS99EC23HsfUN8ZF0QkvZSgpOZGS0blBRFb7r2PPd/8RnKPPKoqOxEBlKCkxkaqzlv37zfTOns2kysKInKPPMqBH/pgkuGKSIroGpTU1Ehz4BXzBQa3PDNsOxVEiEglJSipqRGr8zrayff2DGtTQYSIVFKCkpoaqTpv8p57kslmNe2QiIxJCUpqaqQ58DLZLHNOOkFLUIjImFQkITU15tREf31E0uGJSIopQUnNaWoiERkPDfGJiEgqqQcl4xZneQwRkfFSD0rGpXQD7mB397DlMbauXpN0aCIyQShBybiMdANutr2NruUrkg5NRCYIJSgZFy2PISK1pgQl46LlMUSk1pSgZFxGugFXs0GISDUpQcm4lG7A1WwQIlIrKjOXcdMNuCJSS+pBiYhIKilBiYhIKtV1iM/MDgWuBQ4BHgNOd/eVI2x3GPB1wIAu4AvufkM9YxURkWTVrQdlZpOAW4Gbgd2By4A7zGxqxXZN0XZfd/fdgJOAa6LkJiIiDaKeQ3yHAy3ufpW7D7r7UuBB4ISK7fYAZgIZM8sARSAPDNQxVhERSVg9h/gWAA9VtD0MHFze4O6bzewa4DvAt4Es8DF3r3yt1IgmgRWRNKhnD6oDyFW05YBh8+VEQ3x9wGJgCqHndZGZvaMOMTY8TQIrImlRzwTVS0g45dqAnoq2Y4E3ufvSaCjwbuBG4Iw6xNjwNAmsiKRFPRPUKkJVXrmDovZyc4HJFW15YLBGcUkZTQIrImlRz2tQKwiFD+cA1wDHEcrNb6nY7g7g82b2QeB64FXAB4D31zHWhtU6ayaD3d00t7dvb9MksCKShLr1oNx9ADiSkJi2AOcDx7h7l5mdZ2YPRts9SBjmOwN4FlgCfNrdb61XrI1Mk8CKSFpkisVi0jFUlZntBzy+bNky5syZk3Q4uyRV8YlIta1bt44jjjgCYH93fyLOazRZbIMbLRkpIYlI0jQXXwNTSbmIpJkSVANTSbmIpJkSVANTSbmIpJkSVANrnTWTQm745B4qKReRtBi1SMLM3hp3J+7+y+qEI/U0Y9FC1i5ZCoSeUyGXo9CbY/bRRyUcmYjI2FV8d1U8LgIZYAgoAC3R7wNUzKcnu4ap8+cxd/GJw6r4Zh99lCr4RCQVxkpQnWW/Hw+cQ7h59nfuXjCzQ4BvATfVLjypNZWUi0hajZqg3L239LuZXQIcV776rbvfb2ZnAf9NSFQiIiJVE/dG3amEdZkqde7APiSFNGuEiKRV3OTyn8C3o4le/0i4FvV64ErguzWKTWqsdKNutr1t2I26cxefqCQlIomLW2b+EeA3wK3AeuBJ4GbC8N4naxOa1Jpu1BWRNIvVg3L3bcBpZvYRnl/T6WF3r1xsUHYhfRs3MWn6tGFtulFXRNIi9o26ZtYJvDf6eQJ4q5ntX6O4pA50o66IpFmsHpSZzSfcF7UVOAD4GnAScLSZ/Y2731ezCKVqKgsi2l56IFvuDR+dbtQVkbSJ24P6GvAf7n4Q0A/g7u8FvgdcUaPYpIo2/GIZ/sUr2Pzr39C/qYvck+vZcu997PnmN9LS2cnA05tp6exUgYSIpEbcKr43AB8bof1K4M/VC0dqYevqNaz795spEnpKQ4MD9D25ntZ9ZpN75FEO/NAHkw5RROQF4iaobcBMYHVF+zzCsJ+kWNfyFRTzBbId7WSATEsLQ8Dglmfoa9ZtbCKSTnGH+L4DfMPM3hA9nmFmRwH/CvygJpFJ1fRt3ES2o51iPr+9LdPcTL63RwURIpJacb8+f4YwWexyYDLwKyAPfBM4vzahSbW0zppJMZ9n25PrGSIkp0J/P5lslhmLFiYdnojIiGL1oNw97+7/AuwBHAy8EtjD3T8O7F7D+KQKZixaSCabZco+s0Ny6uklA8w56QQVRIhIasUtMy8Ae7l7F/BgWfu+wCqgozbhSTWUL6uRaW6m9WDNuSci6TfWgoUnAe+OHmaAG8ysv2KzlwBbahSbVJGW1RCRXc1YQ3x3Aj1AadmNbdHvpZ8ewvx8x9QyQBERaUxjrQf1NHA6gJk9AVzh7tvnxTGzZnfPj/xqERGRnRO3zPyrwLVmdl5Z22NmdqOZabl3ERGpurgJ6mrg5cDPytreCxwCfKXaQYmIiMRNUO8C/tHd/1BqcPe7gTOA42oRmIiINLbYy20AraO0T6pGICIiIuXiJqifEKY6ekWpwcxeRhj6u70WgYmISGOLO9XROcCPgPuje6GKhCmP7iAsBy8iIlJVcZd8fwZ4m5ktABYAA8Bqd394Rw5mZocC1xKKKx4DTnf3lSNstzdhItqFQB9wnbtfsCPHEhGRXduOLPneBOwHzAV+CUw1s6k78PpJwK3AzYT5+y4D7hhlH7cCG4BZwGHAqWa2OO6xRERk1xcrQUU9mj8BSwkr6O4JfBp4KLoWFcfhQIu7X+Xug+6+lDCv3wkVx3o9YVn5j7p7n7s/Hr12RczjiIjIBBD3GtTXCJPCvhZ4Omp7D2HJ968B74ixjwXAQxVtDxNmRy/3auAB4LNmdhphiO+b7q77rWLaunoNXctX0LdxE62zNDGsiOya4g7xLQI+5+7bJ4uNpj06H3h9zH10ALmKthxQORPFnsBbgEFCT+pY4FwN8cWzdfUa1i5ZymB3N5OmT2Owu5u1S5aydfWapEMTEdkhcRNUhpHvg5pOKJiIoxeYUtHWRph0tlw/sNXdP+vu/e7+Z+AGQqKSF9G1fAXZ9jaa29vJZDI0t7eTbW+ja7lGSEVk1xI3Qf0Q+Ep0LaoIYGaHAN8gFDTEsQqwiraDovZyDwNtUVFFSdyhyIbXt3ET2bbhndJsWxt9GzclFJGIyPjEPfF/ArgeeDJ6vApoIdzA+4mY+1gBZMzsHOAawhRJhwC3VGx3J9BFSIj/REhq7wPOjHmchtY6ayaD3d00t7dvbyvkcrTOmplgVCIiO27UHpSZfdDMOgDcvcfdTwJeCvw9ofLuZe5+jLtvjXMgdx8AjiQkpi2E61fHuHuXmZ1nZg9G2/UBbyNcf9pAmKD2S+7+w/G+yUYyY9FCCr058r29FItF8r29FHpzzFi0MOnQRER2SKZYLI74hJnlAHP3tRVLvqeame0HPL5s2TLmzJmTdDiJUBWfiKTNunXrOOKIIwD2d/cn4rxmrCG+9YQ1oH5LKJL4pJlVFjQA4O6f28FYpYa0vLuITARjJajTgAuBvyMURrwdGGkF3SKgBCUiIlU11pLv9xLdgGtmjwNvj5aBFxERqbm4k8XuD2BmLdFrMhXPV96AKyIislNiJahofrzrgFdUPJUhDPFlqxyXiIg0uB2Zi+854BggVlm5iIjIzoiboA4GDnP3B2oZjIyfSstFZKKJO9XRQ8DsWgYi46cJYkVkIorbg7oauN7MrgbWUDFBrLvfXu3AJL7yCWKB7f/tWr5CvSgR2WXFTVDfjv77xRGeU5FEwvo2bmLS9GnD2jRBrIjs6uKWmcdeGl7qTxPEishENGqCMrO20v1NZla5qOAwug8qWTMWLWTtkqVA6DkVcjkKvTlmH31UwpGJiIzfWD2jbjMrfQXvAbpH+Cm1S4Kmzp/H3MUn0tLZycDTm2np7GTu4hN1/UlEdmljDfEtIiyLAaC1GlJOE8SKyEQz1lx8d4/0u4iISD2o+EFERFJJCUpERFJJCUpERFIp7o26kkKaf09EJrK4y23MAD4DvBpo4YXrQb2u+qHJWErz72Xb24bNv6fychGZKOL2oG4EXg98Hy23kQqaf09EJrq4CeqtwNEqN08Pzb8nIhNd3CKJZwgLFkpKtM6aSSE3fIYpzb8nIhNJ3B7UhcA1ZnYOIy+3obn46kzz74nIRBc3QV0J7A78epTntdxGHZWq9/K5HP2bN9PUOoXOA/dn9tFH6fqTiEwYcRPU8TWNQmIrr95re8m+23tOKjEXkYkm7npQdwOY2RRgHuHa1aPurpnM60zVeyLSKOLeB5UFLgc+xvP3QQ2Y2U3A2e6er1mEMoyq90SkUcSt4rsMeA9wCjA3+jkF+FtCAYXUiar3RKRRxL0GdQrwfne/vaztP8ysG7gOJam6UfWeiDSKuD2oDuCREdofA6ZXLxx5MVo9V0QaRdwe1ErgLMI1qHJnA7+PezAzOxS4FjiEkNxOd/eVY2y/O3A/cKG73xT3OBOdVs8VkUYQN0H9M3CXmR3O8/dCHQbsB/xNnB2Y2STgVuAqwtRJxwF3mNlL3H20+f2uBfaJGaOIiEwgsYb43P13wKuAOwkFEtOBnwAHufuvYh7rcKDF3a9y90F3Xwo8CJww0sZmdiowFXgg5v5FRGQCib0elLuvBs7diWMtAB6qaHsYOLhyQzPbH7gIeCPws504poiI7KJGTVBm9lvgne7+jJmtBIqjbRtzPagOoHLOvhzQVnHcLGFZj3Pd/Skzi7FrERGZaMbqQf0U6I9+v60Kx+oFplS0tQE9FW0XAO7u/1WFY4qIyC5q1ATl7heXPVwB/MrdB8u3MbPJhJt141gFnFPRdhDw3Yq2E4HZZnZs9LgT+KaZvc7dPxzzWCIisouLew1qBbAX0FXRfgCwhBf2jEbbRyZasuMaQhXfIcAt5Ru5+0Hlj83sT8BVKjMXEWksY12DOhMo9aIywCozq7wO1QH8Mc6B3H3AzI4klI5/DngCOMbdu8zsPOBkd3/5DsYvIiIT1Fg9qOsJ142agH8DLmH4qrpFwvWjZXEP5u7/A7x5hPbLCZPRjvSav4q7fxERmTjGugaVJ7o+ZGaPA/cBne6+JWp7LfBHzWQuIiK1EHcuvi7AgX8pa/sJ8ICZHVj1qEREpOHFTVBXA/fw/DUpgP2B3xAKHkRERKoqboJ6HfA5d99+z5K7bwMuBd5Ui8BERKSxxU1QW4BXjNA+H9Cy7yIiUnVx74O6AbjOzOYCvyNU8L2KMOvDjTWKTUREGljcBHVZtO0FwIyobRPwVeCKGsQlIiINLlaCcvchwuziF5nZdGBgjDWcpMq2rl5D1/IV9G3cROusmcxYtFALForIhBd7uQ0zeyVhyYxs9DgDTAZe7e5n1CY82bp6DWuXLCXb3sak6dMY7O5m7ZKlWuZdRCa8WAnKzM4nzCTRA7QTZpTYLXr69tqEJgBdy1eQbW+jub0dYPt/u5avUIISkQktbhXfGcAn3X0qsIEwyes+hOXfV9YoNgH6Nm5iaGCQrase4pnf/4Gtqx5iaGCQvo2bkg5NRKSm4iaovYAfRr//CXiDuz8FfAp4by0Ck6Bp8mR6Vq9haHCQptZWhgYH6Vm9hqbJk5MOTUSkpnZkqqNp0e+rgUOj358EZlc7KClXDD/Fsp9Sm4jIBBY3Qd1KuA/qrwjrOp1iZm8DPgH8pVbBCQz1D9A+fz5NkyYx1N9P06RJtM+fz1D/QNKhiYjUVNwqvnOBKwmzSfyAsNjgMkLRxOLahCYArbNmMtjdzdQFL9velu/tpWXanglGJSJSe3F7UMcB57v799296O6nAbsD091dVXw1NGPRQgq9OfK9vRSLRfK9vRR6c8xYtDDp0EREaipugvo6ML28wd17tBZU7U2dP4+5i0+kpbOTgac309LZqXugRKQhxE1QvwHeXctAZGSaRUJEGlXca1BDwOVm9hngcWBb+ZPu/rpqByaaRUJEGlvcBPWb6EfqSLNIiEgjGzVBmVlTNEks7n7xaNtJ7fRt3MSk6dOGtWXb2jSLhIg0hLGuQQ2a2czyBjN7q5lpCoM6aZ01k0IuN6ytkMvROmvmKK8QEZk4xkpQmRHabiPMwSd1oBJzEWlkcav4SkZKWlIjKjEXkUYWez0oScbU+fOUkESkIe1oD0pERKQuXqwHdZqZ9VRs/x4ze7p8I3f/ZtUjExGRhjZWgvpf4MyKtqeAf6xoKwJKUCIiUlWjJih336+OcYiIiAyja1AiIpJKda3iM7NDgWuBQ4DHgNPdfeUI270d+AIwD9gEXOHu36pnrCIikqy69aDMbBJhZd6bCWtJXQbcYWZTK7abC/wQuDTa7iTg82b2znrFKiIiyatnD+pwoMXdr4oeLzWzs4ETgOvLttsPWOLut0SPV5rZXcCbgJ/XJ1QREUlaPRPUAuChiraHgYPLG9z9HuCe0mMz2xN4C/C9WgeYFloDSkSkvkUSHUCuoi0HtI32AjPbDfgxYamPW2sXWnqU1oAa7O4etgbU1tVrkg5NRKSu6pmgeoEpFW1tQM8I22Jm84FfAxuB40tLf0x05WtAZTIZmtvbyba30bV8RdKhiYjUVT0T1CrAKtoOitqHMbO3EnpNPyIkp77ah5cOfRs3kW0b3qnUGlAi0ojqeQ1qBZAxs3OAa4DjCOXmt5RvZGYHEpb1ON/dr65jfKnQOmsmg93d21fPBa0BJSKNqW49KHcfAI4kJKYtwPnAMe7eZWbnmdmD0aZnAZ2E0vKesp8v1ivWJGkNKBGRIFMsFpOOoarMbD/g8WXLljFnzpykwxkXVfGJyESzbt06jjjiCID93f2JOK/RelAppDWgREQ0F5+IiKSUEpSIiKSSEpSIiKSSEpSIiKSSiiRSRNV7IiLPUw8qJTQHn4jIcEpQKaE5+EREhlOCSgnNwSciMpwSVEq0zppJITd8NRLNwScijUwJKiU0B5+IyHBKUCkxdf485i4+kZbOTgae3kxLZydzF5+oKj4RaVgqM08RzcEnIvI8JagU0P1PIiIvpASVsNL9T9n2NjLNzTzzxz/x9L33MfUVC9jn2HcrUYlIw9I1qISV7n8qDubpffQxAJraptD7+BO6UVdEGpoSVMJK9z9tW7+eTEszTc3hp5jP60ZdEWloSlAJK93/VNi2jUxzGHEt5vNkp0zRjboi0tCUoBJWuv8p09zM0OAgQ4ODFPMFWmfvrRt1RaShKUElrHT/U/v++zOU2waZDG0HHkBTS4tu1BWRhqYqvoRUlpbvc+wx7HPsMdvbWqbtyeyjj1IVn4g0LCWoBJSXlpcvrTF38Ykc+KEPJh2eiEgqKEEl4InvfJfeRx+nmM/TNGkSU16yL5P33IOu5SvUYxIRiegaVJ1t+MUyuh9yikOFUBiRH6R3zSP0bX5GFXsiImWUoOpsw60/IdPcTKYpS6YpQ1NzMzRl2PbEE6rYExEpowRVZ4PPPktzZycMDVEcKobGTIahgQFV7ImIlNE1qBqrrNbLTpnCUKFAy25TyffmKBbyUITJ06fr+pOISBklqBra8ItlrPv3mynmC2Q72inm82QmTaKweQu0TaF5990Y6u+n2D/APv9wXNLhioikihJUlW34xTI23PoT+p/ezFB/H9nWKbTssTvFfJ5tT65nyj6zadltN/JbtzL47LO07L47e590Anv/9RFJhy4ikipKUFVQGsZ79sFV9G94imxbW3hiqEghl6OppYXmjnaGgP4tW2idOZNDrr4q0ZhFRNJOCWoHbF29hif/6xZ61jwKQMe8l9I6Zx+eXnEXxUKBfHcPxQwM9fVRLA6RyWZDe08PzR3tZJqbKfT00nqwqvVERF5MXROUmR0KXAscAjwGnO7uK8e7XbWVekI9jz1OYds2mlqnMHnaHkCGbevXs239eigM0TR5EtnWKTx3/wNs+e1Kmjs6aO5oZ/DZ5wAotjTBUJFiUwaGMhQLBYrFIoX+fjLNWVXriYjEULcyczObBNwK3AzsDlwG3GFmU8ezXbWVph/a9uR6+rueJp/L0ffUBp574EGee+ABtq3fAIUhKBYpDg1RyOUo5vOhXDyfJ5PJkMlmw86GhqApQ1OmCaL2Qi5HBphz0gmq1hMRiaGePajDgRZ3L118WWpmZwMnANePY7uqKq1sm/vL/9I0qYWm5mb6e3PPb1AohP9mMiFRtbQw1N8PQHFwEIBsRzv557ZSLBRomjSFppYWirkck/fei91fvoAZixYqOYmIxFTPBLUAeKii7WHg4HFuV1V9Gzcxafq0aGivNTQODVGkSIZMeJzJAKEH1dQU2jLRTBBD+TzNHR0UB/MU+rYB0LLbVPY95WRV6ImIjEM9E1QHkKtoywFt49yuqlpnzWSwuzvcSDs4GCWeplJqItPSTLEwFHpPTU0M5fOQydDU0hIWF+zro9DTS3ZKKy85/VQlJRGRnVTPBNULTKloawN6xrldVc1YtJC1S5YyaY89yK17kqF8nkxLMwwVgSKZlmYymSbyvb0haQ0VmbrgIKYf/jZyjzwaZoo4eKaG8UREqqSeCWoVcE5F20HAd8e5XVWVVrbtWr6CQn9/GOrbc9r2Kr6BLVsobNvG5Fmz6Dxw/+GJSL0lEZGqq2eCWgFkzOwc4BrgOEIZ+S3j3K7qps6fp96PiEhK1K3M3N0HgCMJCWcLcD5wjLt3mdl5Zvbgi21Xr1hFRCR5db1R193/B3jzCO2XA5e/2HYiItI4tB6UiIikkhKUiIik0kScLDYL8NRTTyUdh4iIRMrOydm4r5mICWpvgJNPPjnpOERE5IX2Bh6Ns+FETFArgbcAG4BCwrGIiEiQJSSn2CtTZIrFYu3CERERGScVSYiISCopQYmISCopQYmISCopQYmISCopQYmISCopQYmISCopQYmISCopQYmISCpNxJkkYjGzQ4FrCYshPgac7u6x73BuBGb2duALwDxgE3CFu38r2ajSy8x2B+4HLnT3mxIOJ3XMbG/gX4GFQB9wnbtfkGxU6WJmhwFfBwzoAr7g7jckG1VyGrIHZWaTgFuBm4HdgcuAO8xsaqKBpYiZzQV+CFxK+BudBHzezN6ZaGDpdi2wT9JBpNithCnIZgGHAaea2eJkQ0oPM2si/I2+7u67Ef6fuyb6Mt2QGrUHdTjQ4u5XRY+XmtnZwAnA9YlFlS77AUvc/Zbo8Uozuwt4E/DzpIJKKzM7FZgKPJB0LGlkZq8HDgDe5O6DwONmdjiwLdHA0mUPYCaQMbMMUATywECiUSWoIXtQwALgoYq2h4GDE4glldz9Hnf/UOmxme1JmIT3j8lFlU5mtj9wEXB60rGk2KsJyfuzZvakmT0KvNvdNyQcV2q4+2bgGuA7wCBhUtXz3L3yXNUwGjVBdQC5irYc0JZALKlnZrsBPwZ+QxiCkIiZZYHvA+e6uxYhG13pC84goSd1LHCuhvieFw3x9QGLgSmEkZ6LzOwdScaVpEYd4usl/AMo1wb0JBBLqpnZfEJSWgWc7O5DCYeUNhcA7u7/lXQgKdcPbHX3z0aP/2xmNxAS1ZLEokqXYwlDoJ+MHt9tZjcCZwB3JBdWchq1B7WKUCVT7qCoXSJm9lZCr+lHwPHu3pdwSGl0InC8mT1rZs8Shom/aWbfTDiutHkYaIsKlEoa9QvyaOYCkyva8oReZ0Nq1H8gKwgXIs8hjPkeRyg3v2XMVzUQMzsQuA04392vTjqetHL3g8ofm9mfgKtUZv4CdxLKpr9iZv9E+IL4PuDMRKNKlzsIlbIfJBRrvQr4APD+RKNKUEP2oNx9ADiSkJi2AOcDx7h7V6KBpctZQCfhf5iesp8vJh2Y7Hqi3vfbCNefNgA/A77k7j9MNLAUcfeVAQvQAAAJH0lEQVQHCcN8ZwDPEoY+P+3uDXvdVyvqiohIKjVkD0pERNJPCUpERFJJCUpERFJJCUpERFJJCUpERFJJCUpERFJJCUrGZGbtZnaJma02s21m9riZfSVa+6hax+gws/dVaV8tZvbhauyrGszsNDN7uk7HOsfMLo5+v8vMvrwT+xrXZ2JmvzKzN4/3uGlhZq8xs2VJx9HolKBkVNH6WL8G/hr4KGEW+A8B7ySsn1U5n+F4/RPVm1FgMfC5Ku1rlxGt3/Ux4Ioq7XKHP5NoUuGDCP9mdmnu/jvgmWgZFUlIo051JPF8gfAlZpG7l9btedzMVgGPAP8IVGPOuUwV9lGLfe1KPgX8p7tXa8Lj8fwdjwDucfd8lWJI2teBb5vZ9zRJcjKUoGREZjYZOBn457LkBIC7rzWzhYBH22aAjwBnA/sSJgY9z91vj56/ibCcyRTgeGAzcIO7X2pmpxHWUsLMisD+hDnbrgCOAWYATwHfcvdLo+2ywHmEecr2JKybc3a07bfL9rXQ3e+qeF+fJcy7+Bhh/aZ+4P8C57h7IYq1w92PL3vNE8CX3f2a6PUHA2sIPYwB4BLgD4TlzA8AfkmY+X1L2T7OI/RKiLa7yN0L0XOvB64krJm0ljAP25fdfSj6+3yUMJHxUYTpgS6teE8dwGmEnm25GWb234Ql1h8DznL3FWWvO5fwuU0D/kRYMuTX4/lMIm8nmnU7+jdxMWHanlbg34C/Ar7j7jeZWftY+4s+h62EpXFOAJ4BPkGYPPUKwsJ+twGnuvvAeD6XaC7OMwmLc/YAPwXOLEvy90bH/3u0zEwiNMQnozmAsELsypGedPf7ogXWICSLi4ELCSf/HwE/rliq+v3AOuA1wA3AJWb2auBm4CvAn4G9CSfoK4E3EE5eRvgmW9qe6DgfBT4OvDJ6zU+BX0VtW6J93TfKe3sXYRn7N0T7Ois6VlzvIpy4XkVY5v1KwqTDHwbeQUg0Hy/bfhoheSwCTiGctD8NYGYzCSsU/4xwgv1oFM+nyl7/SqA7Ot53R4jnbUCB8P7LnQLcAxxKmAj559HiipjZGYQhwQ9H+78dWB49P57PhOi93xn9/s+Ek//7CaswWxRnSZz9nQE8Gv1dfg7cCHwS+D+EWeSPBt5Ttn3sz8XMTgI+S0h68wgJ/ujomABEvaY7CfN2SgLUg5LR7BH997mxNoq+KX8cuMzdl0bNn416BZ8i9MIAHnX3C6LfL42+vb7a3X9vZj1AvrTgn5n9P8K36T9E219hZhcCC8zsD4QT36WlNZjM7CzCyaYzirf4IosHbgM+HE0a7GZ2JuHkFXfi0j7g4+6eN7Orgc8A17j7PVE8twMvL9u+ALzH3dcS1kG6NPrbXEZIRivd/ZJo2zVm9i+EE/YXyvZx8Rjv6TXAQ+5eObHmne5+efT7+WZ2JCFhnE/4UvFpd/9p9Pzl0RLsZ7n7uTvymQC/N7MDgGZ39+j5jwCXuPtPote/h5DoSsbcX9T2iLtfFr3+W4Qe7wWl15jZbxn+d96Rz2U9cJq73xY9/ouZ3R0dv9wqwgSukgAlKBlNqfJsjzG3CkMt03nht/d7Cd90S9ZUPN8NtIyyz+8D7zKz9wLzCUNDHUA2OtYMynp27t5NNHxmVrnM14j+EiWnkq1jxDKSJ8qus5RWZn6s7Pm+KM6SJ6PkVLISmB1VQr4cODxKCCVNwBQzm1Y6xosk3Fk8/3mVq/xMfge8PBoS3Be4Pjrxl0wmDHmOZKzPBMp6T2Y2HZjN8M9oi5mt3oH9QbjOWTLa37l8/aTYn4u73x1V6l1KKOx4OaEnV9lD3Uz4Ny4J0BCfjOYRwv+crx3pSTO70sw+SeiNjCTD8H9fA6NsM5J/IxRfbCOcMN5AWH6gfD87Mw3/WLGMtN/KL3IjLSA31kX0QsXj0t9lMNr3Dwkn6NLPIYRhp1LvdbS/cfmxR/pbjnTcAZ5PAqdWHPdlhCrNkYz1mUDZ9See//uMdX55sf2V76fcWH/n2NtH19nuIVzD/G/CkOGPR9g0ywv/jlInSlAyomj8/XvAR8ystfw5M9uPMFbf7+5bCcMlb6jYxRsJxRJxbE8KZtYJvBc4xd3Pc/ebCSfV3YCMuz8HbCJcZyi9ZoqZbYyGFXd2/ZjSsUr77mDnv0HPMbPynuhhhF5cL/AQcJC7P1L6ISSKCxn7ZFzuKUKvstIhpV+iodjXAauiv+FTwD4Vxz2b5wstYn8mUdHK4cAygGj/6wjDpqV97EZIui+6v5jveWedBVzh7h929xuB+6P4Ko8/nfC3kgRoiE/Gcgnwd8AKM7uIMEx3KPAlwv/Q10XbfQH4nJmtI1w/OIFwojs85nF6gL2i6xhrgV7gWDN7nDBUdAXhxFEazvkqcEFUXbeacE3lOeCPwBygw8wWAI+NY5n6lcD7zOwYQvK4mJ3/Bp0FlpjZpwkVcedFPwDfAD5qZl+Pfp8LfAu4Narii7P/3wOfNrNsqTIwcrSZfYJQQPIBQrVa6baALwEXmtkGwnteTEhQpUKGHflMXkv4W5eKZiB8Rp8xs8cIw2yXAu2ExNf3Ivurh83AwujfSQY4h3D96YGK7Q5llEIhqT31oGRUUTnumwnXLr5FuGD8FUKV3t+UnfyvIZzwvkT4H/xo4F2li9Mx/CfhmtQqwlDTYsLNwasIwz93EKrMSt/IryBUdN1ASEp7A38XXVdaFrX9kZBcd9T3gJui495DKL/+f+PYT7mHCEnkl4Sk/lXC3xN3X0dI5q8hVM19j1BF9/ER9zSyFYST7Gsq2q8mlEjfT7hH6W/dfVP03NeALxM+s1XAPwDHu3vpve7IZ/J2nq/eK7mK8Hf8DuFa2CPAX4ABdx98kf3Vw8cIyfJ3wC8IifHzDO+ZNxEqEG8baQdSe1pRV2QCMLNvAAV3/2jSsQCY2d8CfyirAmwmFHIc5e6/TDS4mMzsnYQvX6YbdZOhHpTIxPAl4LjoWk8avA/4vpkdbGbzCD22Z4HfJBvWDvkwcLmSU3KUoEQmAHf/C2Ho8JNJxxI5m3Bd8JeE2RwOBN7h7qOVsaeKmb2WULRxU8KhNDQN8YmISCqpByUiIqmkBCUiIqmkBCUiIqmkBCUiIqmkBCUiIqn0/wF9gM2AvTO1YgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_sweep_frame(frame)\n",
    "\n",
    "decorate(xlabel='Contact number (beta/gamma)',\n",
    "         ylabel='Fraction infected',\n",
    "         legend=False)\n",
    "\n",
    "savefig('figs/chap06-fig03.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns out that the ratio `beta/gamma`, called the \"contact number\" is sufficient to predict the total number of infections; we don't have to know `beta` and `gamma` separately.\n",
    "\n",
    "We can see that in the previous plot: when we plot the fraction infected versus the contact number, the results fall close to a curve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the book we figured out the relationship between $c$ and $s_{\\infty}$ analytically.  Now we can compute it for a range of values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_inf_array = linspace(0.0001, 0.9999, 101);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_array = log(s_inf_array) / (s_inf_array - 1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`total_infected` is the change in $s$ from the beginning to the end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "frac_infected = 1 - s_inf_array\n",
    "frac_infected_series = Series(frac_infected, index=c_array);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the analytic results and compare them to the simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap06-fig04.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8XGW9+PHPmSWZmezN0tANuj5tWQVvhQpKq4grlEXZFEVQ5OLGFa8IClLW64L8pCiyXEC0t1UECrKI0qogIlUKAi3frli6JWnTJpNZklnO748zSabTNJ22mSWZ7/v1yqszz5w55zsJnO88z3nO97Fs20YppZQqNq5CB6CUUkoNRBOUUkqpoqQJSimlVFHSBKWUUqooaYJSSilVlDRBKaWUKkqaoJRSShUlTVBKKaWKkiYopZRSRUkTlFJKqaLkKXQAQ80YUw78B7AVSBQ4HKWUUg43cAiwXES6s3nDiEtQOMnp+UIHoZRSakAnAS9ks+FITFBbAX71q1/R3Nxc6FiUUkoB27Zt44ILLoDUOTobIzFBJQCam5sZN25coWNRSim1u6wvvegkCaWUUkVJE5RSSqmipAlKKaVUUdIEpZRSqigVZJKEMWYW8DsRadrL6xOA+4DjgVbgKyLy1FAcO5lMsmnTJkKh0FDsTu0Hr9dLU1MT1dXVhQ5FKTUM5DVBGWMs4GLgh/vYdBHwN+BjwInAY8aYY0Rk/cHGsH37dizLwhiDy6UdyHyxbZtIJMLmzZsBNEkpVcRs28a2nX+TqX/dLgu3O7/nzHz3oK7HSTo3At8ZaANjzDTg3cApItIDLDXGPI6T2K452AB27drFYYcdpskpzyzLIhAIMHbsWLZs2aIJSuXEQCfWpG2DDcm9vNbbZu+xTaotaQMZryXTt7exATuZ9jzttaRN3/F3j6U/jt32S+/j/ljB+TfpvCEt/szYB//8/Z9lz8/f+5l6j5+p3Ovm5OPGMXV8XZ7+mvlPUHeJyLXGmJMH2WYmsFFE0sfg3gJmDUUAiUQCr9c7FLtSB8Dv9xOLxQodRtFLJm1iiSR20h7wxLrHSTTjRGvbNom9nVjJeJxxYrWxsZOQyDiZZh4ntHkrHW+8QU9HJ57qGipnzsDXPLrvpJ5+YrVtO/VZnP33nnwzP0fvyT6Z3DOxRHe0E/r3RuLhMK5AAN+48Xhra/aIayjFOjqIbN5CIhLB7ffjHzsGb03NkB5juByzO5Zg9b93jtwEJSJbstisEghntIWBwFDFYVnWUO1K7adS/d3btk13T4JId5xId5xwd5xINE6kx/m373nq9WhPfL/2n++TWqyjg66167A8HixPPXZnHPuF16mcEs7JcXc7XnkAuztORNZQOWVyzj5n+jFdPh/JWIyutetK6JhTKK+rwWVZWBZU+L0cNbUxJzHsTTFWkggB/oy2ANBVgFiU2qtYPEE4Lan0Pm5/+x3a/rWSro4QiYpqXOMOJVlR5fQIchFHAU5qkc1bnON5nVOI5fWQTLXn4pj7czzLsrAAl8s5sVqpE6xzorVwZbSxl9d2rnkDv8/GXWZhkcAqA7sniWf7Rppmnti3PfS/z+Xqff/u+0p/7tqtrff4zvOtT6wgWRXB4/NhWTEsIBmJ4O1ax/hTz3Ridlm7x5/2WXs/R//jfX/+Dff9L/FRXZRVOH0AC0iEQpTF32LyWV8c8r/l/ijGBLUSmGCM8YtIJNU2PdVesq699lp+/etf88QTTzB16tSD3t9VV11FVVUV11xzYJf1Hn/8cRYuXMiiRYsOOpZilkgk6Qj1sCvY7fx0RdkV7GFXVzfh6J5Dlbt90y+rwY7EsVet3u9kYVkWXo/LOcFg4XL1ngz3PEHtXP0mdeX0n0i9FnbMxtP2Dg1Tx/Vt57KA3U5QA5xYSTuJuvY8sfaezDf+60+U1VRjWXFczm7Btol3bGfSuz+0x0m498QKu5849zixknmid16TN56mvKmu/704x4vt2MHhZ5yzRzIYCm++8GvKGuuxrP6/s23b9GzfxuGzPzMkx8jkCW6mrCHjmOVuerZvYcr42pwcM9nWiq+hnvRfm1URINrSmpPj7Y+iS1AiIsaY14CbjDHfBmYDpwMnFDaywgmFQjz99NOcccYZPPTQQ8yfP7/QIXHaaadx2mmnFTqMIWHbNsFwjI6u7r5EtLMryq5gN8FwbL+uawz2Tb+iYRSBcg++cg/+cg8Bn/Nv+k9vm6/M03dC35c3X3qYsqZ6LKtnt8/Us30rh7//0/v1u8hW3dgAseBOPIGKvrZ4KIR3bBWTJ9YP+fHqD6knFgziqdj9eL7mRrye3Ex48o1u2uOYiXAY3+gB747RY+ZAUSQoY8wFwM9FpDLVdBZwN849UNuBi0XkjULFV2i/+93vmDFjBpdccglnn3023/jGN6ipqeGRRx7hiSeeoLm5mWeffZa6ujouueQSzj33XABefvllbr/9djZs2EA0GmX27NnceuutVFVV9e1727ZtzJkzh9///vdMmDABgMcee4xf/OIXPPLII9x999388pe/pLu7mylTpvCtb32Lo446ikceeYQHH3yQJUuW0N7ezre//W1eeeUVKioqOOGEE7juuuvw+XwF+X3tSygSo6U9TEt7iJb2CK07w/TE9n/pMJfLIlDuIeDz9iWY7a+1UVlTic8do9yyKXfZlFsJrPbtHDXvvBx8msKcYBrnzuGdhU7v2R0IkAiHSYTCjDk9N19a8n08PWZuj5mtgiQoEfkTUJv2/FfAr9KevwN8JB+xrJBWXl65jVg8mfNjeT0uZs1s5l1m/04cixcv5qKLLmLy5MkcccQRPPzww1x88cUAvPjii1x//fXccMMNPProo8yfP5+PfvSjeDweLr/8cq6//no++tGP0traymc/+1kWL17MJZdc0rfv5uZmZs2axZNPPslll10GOAnxtNNO44033uC+++5jyZIlNDU1cccdd3DbbbfxwAMP7BbfnXfeSVVVFS+++CKhUIgLL7yQZ555hnnz5h3cL2wI7FglrPnDC2xtCxKsrCc6ZiI9vsp9vzHFsiyqAl5qK8uprXJ+airLqa0spypQtkcvZ90rPqdn4cvoWTSPnGQBUD1tKuPPP5e2pcuItrTiG93EmNNPo3rawQ8/F8Px9Ji5PWa2iqIHVUivrm7LS3ICiMWTvLq6bb8S1Ouvv87WrVs59dRTATjvvPP44Q9/yEUXXQRAY2NjX49p3rx5fOc732Hr1q1MmTKFhx9+mEMPPZRwOExbWxujRo2itXXPceXTTjuNBx54gMsuu4z29nZefvllbr31VoLBIKFQiEceeYRTTjmFr3zlKwPeP1ZZWclf//pXnn32WWbPns1jjz025PeZda5es9v/QI1z5wz4P5Bt27TtjLCxJciaN9az7p+rwOPF8jRhB+PYbw58PchX5qGuqj8J9SakmspyPPtxc2IpJIv04+bzJJbv4+kxC6/kE9Qx0xrz2oM6Ztr+TdNcvHgxwWCQOXPmAE6ppvb2dpYuXQpAfX3/eH/v/V3JZBK3281f/vIX7r//fpLJJNOnT6ezs3PA6ymnnnoq8+fPZ82aNbz88svMmjWLhoYGGhoauPPOO7n//vv52c9+Rn19PZdffjmf/OQnd3v/5ZdfjsvlYsGCBVx55ZUcd9xxzJ8/n0mTJu3XZ92bztVreGfhItwVAcoanGsR7yxcxPjzz6V62lTC0RjvtATZuC3IxpYgkW5ninbnyrdhgOtBPVu2cOjksYweFXB+6gNU+r1DcnG9VJKFUvlQ8gnqXaZpv4fc8qWrq4snn3ySe+65hylTpvS133XXXTz00EOcfvrpe33vihUruP322/nNb37Tlyh6h/AyVVZWMnfuXH7/+9+zfPlyzjrrLABaW1sZNWoUDzzwAJFIhGeeeYarrrqKE088cbf3r169mnPPPZevfe1rbN26lVtuuYX58+fvMRR4oNqWLsNdEei7xuKpqKA7afHSky/Qtdli647QgIk3EYng8vmodicZ5UlQ701Q54nh37mRI+fk5noQaLJQaqhovZ8i9vjjj9Pc3MwJJ5xAY2Nj388555zDSy+9RCQS2et7g8EgLpeL8vJykskkTz/9NM8///xeqzicfvrpPP3006xatYoPfvCDAKxbt45LLrmE1atX4/f7qa+vp6ysDL9/99vUHnzwQW666SZCoRD19fX4fD5qhvBemGhLK+5AgFgS/h318peOAE9Fmvl7i82W7V17JCd/uQczoY7ZzfAx/zY+PKqLWdURJvt7qO7pIpDD60FKqaFT8j2oYrZ48WI+/vGP79E+bdo0Dj/8cG655RYmT5484HtPOukkPvGJTzBv3jxcLhczZszgU5/6FCIy4PYnnngiV199NSeffDKBgHPD3gknnMCll17KpZdeys6dOxkzZgw//vGPqa3d/X6Mb3/721x77bXMmTOHeDzOrFmzuP766w/y0/eLjGpmxY4km+3KvrWik/E47lSitCyLQ+oDTGiuZsLoKhrr/FiWRWfdSbyzcBFxl110s5OUUvtmDXXtqkIzxhwGbHjuuecYN27cHq+vWrWKGTNm5D2u4eC0007jqquuYvbs2Tk9TvrfYG+THxKJJOs2d/D62u1sXL8lrayOBzsex47HmXTsDI48dgpTxtUS8A1cXzHbyRVKqdzatGkTH/jABwAmisjb2bxHe1CKjRs38tJLLxEKhTj++OPzdtyBJj9s+NVigu/7GKvDZX2VGrw1NVROmezc7NodZFJ9Gcd+8ATGHr3vLxp6PUip4UsTlOL73/8+K1as4Ic//GFelyFJn/wQS8JaRvFW0iL2wltUz+xPPi6XxeFHHsaR895Nc32gZAvOKlVqNEEpFixYUJDjRltacdfXI+EyVobLidkWttspjglQ6fdyxOQGZk4ctdchPKXUyKUJShWEbdtsrxnDilaLsLu8vz0ep7rSx5zjxjP90Lq8r+CplCoeJZmgbNvWYaICiUci9HQG+fk1d7HdCpCIduP2u7E8HgLJbozVzuyzP0HtpKEvOKqUGl5KLkG53W5isRhlZWWFDmXES0SjxDqD2PEYlseLy1dOe3sHu8IJdpTVQNyp+OAliUm0MHNMBc0fKI4aYEqpwiu5BFVbW0tLSwtjx47N64SAUpOIRulpbweXC8vjIRZPsKutgy2trby1ssVZ58frYUqZzVENPcz8zysKHbJSqsiUXIJqaGhg06ZNe71hVQ2NWDAISRtcFnHboicJka4oK1/ZSEfUos6T5N1VYWrdSXradhQ6XKVUESq5BOVyufrWPVK58+b1N8KoBpZ3BdjS4/xn1tO+Cyvew3ETq5nm78FlQTxUHAujKaWKj45xqZyI1jfzh7byvuQEUOe3OKl7HVPsnVjYxEMhEqEwjXPnFDBSpVSx0gSlhty6Tbt4MTCNYI9NMhbHtmGyFWRO+XZmnDsPb1UVPdt34K2q6lsyQymlMpXcEJ/KHdu2eW1NG3/911asqmoqp0wmtmULx8Q3M7m5msa5qWT0wQ8UOlSl1DCgCUoNCdu2eeG1Lby2pq2vrWFsEx/95Czqa/yDvFMppQamCUodtGTSZtk/32HV2+19bYfUV/DR907EX67/iSmlDoyePdR+S1/CoqypiTfHHMOmnv4bnyePq+WUWRPwaJkipdRB0DOI2i+9S2TEgkE89fX8rRVe/9sbxDo6AJg5cRSnvudQTU5KqYOmZxG1X3qXyHAHKnglFGATlVgeD5HNWzhqSgNzjhuPy6V1DpVSB08TlNov0ZZW3IEAb4TL2RB1lsCwPB4OjW3npGPGahFepdSQ0QSl9otvdBOrd9qsCvcvkXGoq4v3jC3X5KSUGlKaoNR+6T7mBFZ0+fpuwB1thzk60ULTB7QahFJqaGmCUlnb2RnlhW1QMWUyLq+Xmp4gJzYlOfQCrQahlBp6Os1cZaUnluCpF9+mJ5bAW1NDXXMDn/rgNF2KXSmVM9qDUvtk2zZL//EOO4NRADxuFx977yRNTkqpnNIEpfbpjfU7WLtpV9/zOceNo7FOyxcppXIrr0N8xpijgbuAo4D1wOdFZPkA2x0P/AQwQBtwq4jcm89YlWNHR4QXXt3c9/yIyQ2YQ0cVMCKlVKnIWw/KGFMGLAEWA7XATcCzxpjqjO1cqe1+IiI1wHnAglRyU3mUSCT5w8sbSSRtABpq/Zx49JgCR6WUKhX5HOI7GfCKyO0iEhORRcCbwDkZ29UBTYBljLEAG4gDPXmMVQEvr9zG9l0RwLnu9CEtYaSUyqN8DvHNBFZltL0FHJneICI7jDELgAeB+wE38DURyXyvyoHeQrBbt+7kedcEfGPH4K2pYfZRhzCq2lfo8JRSJSSfX4crgXBGWxgIpDekhviiwPmAH6fndZ0x5kN5iLGk9RaC7e4M8pp3DIlYjK6162iwohw5uaHQ4SmlSkw+E1QIJ+GkCwBdGW1nAu8VkUWpocA/A/cBl+YhxpLWWwh2vauOjoQbl9eDx+tmRsvrWsZIKZV3+UxQK3Fm5aWbnmpPNx4oz2iLA7EcxaVSoi2tRMsreDPUP5R3ZHUCz/aWAkallCpV+bwGtQxn4sMVwALgLJzp5o9mbPcscIsx5ovAPcCxwBeAS/IYa0nyjW7i5RaLROp5jTvJJHsn5aObChqXUqo05a0HJSI9wEdwElM7cA0wT0TajDFXG2PeTG33Js4w36XALmAhcJWILMlXrKUqfuxsNkXcfYVgj3a3Y4fDNM7VQrBKqfzL6426IvIGcOIA7TcDN6c9fwp4Ko+hlTzbtnm1w0vllMlENm9hTE87hzT7aJz7ES0Eq5QqCC0WW8J6p5RHW1rZVj2WLVWT8NbU4B9Vx9kfnk5loKzQISqlSpjedVmieqeUx4JB3PX1vLbTomvtOmIdHRwzrVGTk1Kq4DRBlajeKeWeigo2RMuJusuwPB7srZs51uikCKVU4e11iM8Y875sdyIifxmacFS+RFtaKWuoJ27Tt3y75fEwuaeFMq+7wNEppdTg16D+lPHcBiwgCSQAb+pxDxnVIFTx841uIhYMss6qI2o7N+H6kj3MGFNR4MiUUsox2BBfVdrPRcDrwAlAuYiUA8cAy4Erch2kGnqNc+fQ0xVmVYcL24ZkLM7UZDvNH9Ap5Uqp4rDXHpSIhHofG2NuAM5KX7tJRP5ljLkceBr4eU6jVEOuetpUInM+Qc+fV5GMRKis9HHS2R/XKeVKqaKR7TTzapyq4pmq9mMfqogkkzYSLqN65gwATjp6LHXTGgsclVJK9cs2uTwM3J8qU7QC51rUe4DbgF/kKDaVQ+s276Iz5Cyx5SvzMHOSrpKrlCou2SaorwA/w1nptvc9MeBe4Js5iEvlkG3bvLq6jVhHB5HNWxgda2Xj1hdpnDtHh/iUUkUjq/ugRCQiIp8DGnB6Tu8BGkTkyyISz2F8Kge27QizacNWutaug3gP00eXEwsGeWfhIjpXryl0eEopBezHjbrGmCrgM6mft4H3GWMm5igulUP/WttGZPMWLI+HiZXgd4OnogJ3RYC2pcsKHZ5SSgFZJihjzDRAgK8Cl+FMmjgPeM0YMzt34amh1hWJsW5TB4lIBMvjYaq/u+81dyBAtKW1gNEppVS/bHtQ/w/4tYhMB7oBROQzwEPAD3IUmxpCnavXsO6uu3nmpp+y682VWG439USo9ST7tkmEw/h07SelVJHIdpLECcDXBmi/DXht6MJRubD1j8+x6f8Wk0wkWFdzNHFfGDueYHx3C/GKAO5AgEQ4TCIUZszppxU6XKWUArLvQUWAgb5aTwU6hy4cNdQ6V69h0/8txga2VzQRTrpIhEL4PHBovQ9vVRU923fgrapi/Pnn6iw+pVTRyLYH9SBwpzHmS6nnjcaYI3CG/n6Vk8jUkGhbugw7nsBdWcE77losy4WdhLHhVqwATP7SFwsdolJKDSjbBPUdnGKxS4Fy4G9AHPgpztLtqkhFW1pxV1YQjdu0eCoBsFwWYzq34Dv8mAJHp5RSe5ftfVBxEfk2UAccCbwLqBORrwO1OYxPHSTf6CbKR41ik11J0rYBqI0FqXLFaZyrhWGVUsUrqx6UMSYBNItIG/BmWvsEYCVQmZvw1MFqnDuHdxYuYtuoCVgRGzsWY3x8J+POO0evNymlitpgCxaeB5yRemoB9xpjujM2OxRoz1FsaghUT5tKxelnEXr8NdxEKKuv4/1nf5iGmabQoSml1KAG60H9ATgFJzmBM5Mvkva6DfwdeCAnkakhs8Vd21e1fOr4WhpmHlbYgJRSKguDrQe1Hfg8gDHmbeAHIhLufd0Y49E6fMXPtm3WbtrV93zahLoCRqOUUtnL9j6oHwN3GWOuTmtbb4y5zxijy70XsW07wn3LapSXuZkwuqrAESmlVHayTVB3AIcDz6S1fQY4CvjRUAelhk5672ny2Brc7qzrAyulVEFle7b6OHCRiLzS2yAifwYuBc7KRWDq4Nm2zbq0BDV1vA7vKaWGj/35Ou3bS3vZUASihl5Le5iuSAxwVs0d06h3Ayilho9sE9QTOKWOjuhtMMbMwBn6eyoXgamDt25zR9/jSWOrcbusQbZWSqnikm2poyuAx4B/pe6FsnFKHj2Lsxy8KjK2bbN+twSlBT+UUsNLVglKRHYC7zfGzARmAj3AahF5K5fBqQPX3hmlo8u5r9rrcTGuSYf3lFLDS7Y9KIwxLuAwYDxwPzDNGFMtIlkvt2GMORq4C2f233rg8yKyfIDtDgF+BswBosDdIvLdbI+jYMOW/j/LhOZqPDp7Tyk1zGS75PshwKvAIpwVdEcBVwGrUteistlHGbAEWIxTYPYm4FljTPUAmy8BtgKjgeOBzxpjzs/mOMrx9tb+BDVpzEC/YqWUKm77s+T7SqCR/nJHnwZeSr2WjZMBr4jcLiIxEVmEU3j2nPSNjDHvASYBXxWRqIhsSL13WZbHKXnhaIyWdqfoh2VZHNqsCUopNfxkO8Q3F3ifiHQb4xQZFZGwMeYanHp82ZgJrMpoewtn+Y50xwGvA98zxnwOZ4jvpyKiNwTvQ+fqNbQtXcbqzUE6PGPwjx3DhElj8JVnPZKrlFJFI9selMXA90E14EyYyEYlEM5oCwOZpZJGAScBMZye1JnAlTrEN7jO1Wt4Z+EiYsEg2wMNJGMxutauozEZLHRoSil1QLJNUL8FfpS6FmUDGGOOAu7EuV6UjRDgz2gLAF0Zbd1Ap4h8T0S6ReQ14F6cRKX2om3pMtwVAVyBClpiHlxeD5bHQ2DVK/t+s1JKFaFsE9R/AduAzTg9oZXACmBD6rVsrAQyFyGanmpP9xYQSE2q6KVjVPsQbWnFHQiwK+6mx3ZuyPWXufHt2FbgyJRS6sAMtmDhF4GFItIlIl3AealrTjNS71slIqv341jLAMsYcwWwAKeG31HAoxnb/QFow+mxfQMnqV0MXLYfxyo5vtFNxIJBtjGqr63JDuFvbipgVEopdeAG60HdDtSBs+S7MaZRRNaLyJMismQ/kxMi0gN8BCcxtQPXAPNEpM0Yc7Ux5s3UdlHg/TjXn7biVFD/voj8dn8/XClpnDuHRCjMllAS24ZkLE5DvJPGuXMKHZpSSh2QwYbOtuCsAfUyziSJbxpjMq8XASAi87M5mIi8AZw4QPvNwM1pz9cDH8tmn8pRPW0qzed8il2/eYVkOILb7+ddZ32c6mlTCx2aUkodkMES1OeAa3EShY2z/PtAK+jaQFYJSuVWV+0hVM5w7ptuqPXTfETmJT+llBo+Blvy/QXgQwDGmA3AKall4FWR2tTW38HV2ntKqeEu22KxEwGMMd7Ue6yM1zPvb1IFsKm1/56nsbr2k1JqmMsqQaXKD90NHJHxkoUzxOce4rjUfor2xNm+Kwo45Y10cUKl1HCX7f1F/w/oAOYBWVcvV/mzdXsI27YBaKz1U+7V7wxKqeEt2wR1JHC8iLyey2DUgdvSFup7PFavPymlRoBsE9QqYAxOEVdVZDpXr+H1J//Ozs4e3H4/NePcOH8upZQavrJNUHcA9xhj7gDWkFEgVkSeGurAVHY6V69hw68WsyMxAZfPhx2LEXtmCZ31Pr0HSik1rGWboO5P/fs/A7ymkyQKqG3pMjp81Vg9zp+y2ufCX+6nbekyTVBKqWEt22nmul54kYq2tNLu7x/Oa/TGcQcCRFtaCxiVUkodvMGKxQZ6728yxmSu2bQbvQ+qcHyjm2jbluy7M63BmyARDuMbrUVilVLD22A9o6Axpvcs1wUEB/jpbVcFUj/nZLZHLZKxOLYNtbFOEqGwFolVSg17gw3xzcWpOg6gZ7silWgej2/KFCKbt+CNhqgZXUHTGVokVik1/A1Wi+/PAz1WxWXbjjDemhq8NTVMHlvDlNkTCx2SUkoNCZ38MMxt29F/g+7o+ooCRqKUUkNLE9Qw19LePz+ledSgc1mUUmpY0QQ1jPXEEuwMdgPgsiwa6/wFjkgppYaOJqhhrHVnuK9A7KgaH16P3i+tlBo5sl1uoxH4DnAc4GXP9aBmDX1oajCdq9fw6hN/ZWdbErffz8TjphQ6JKWUGlLZljq6D3gP8Et0uY2C61y9hncWLqLNasblqyAZi2G/9Bc6p9Xq9HKl1IiRbYJ6H3C6TjcvDm1Ll+GuCNAZDWAlwfJ6qC/3aP09pdSIku01qJ04CxaqIhBtaSXuqyCUdEZaXUBdZZnW31NKjSjZ9qCuBRYYY65g4OU2tBZfHvlGN7FpZ7Tvea0nARGtv6eUGlmy7UHdhnMN6iVgB3vW5FN51Dh3Du2hRF/9vepEWOvvKaVGnGx7UGfnNAqVtc7Va2hbuowdUYjHu7DcHhrGlTH+Y+fq9Sel1IiS7XpQfwYwxviBqTg9r3Uior2nPOqdveeuCBCuGY2nO4kdjzP5fbM0OSmlRpyshviMMW5jzP/gTJZYAbwCtBlj7jLGZNsLUwepd/ae7a+gK+nC5fXg8nhIvPzXQoemlFJDLttrUDcBnwYuBManfi4EPoozgULlQbSlFXcgQEeiv2JEjc8i1qqz95RSI0+2vZ8LgUtE5Km0tl8bY4LA3WiSygvf6CZiwSC7XHV9bdWJiM7eU0qNSNn2oCqBtQO0rwcahi4cNZjGuXNIhMLsCMWwbUgFuBjeAAAaEklEQVTG4lTFQjp7Tyk1ImWboJYDlw/Q/mXgn9kezBhztDHmb8aYkDHmdWPMf+xj+1pjzEZjzOeyPcZIVj1tKuPPP5cuTwXJaBSX14v5xCk6QUIpNSJlO8T3LeBPxpiTce6FAjgeOAz4cDY7MMaUAUuA23FKJ50FPGuMOVRE9lbf7y5gbJYxloSqqVNgRoS6eBKAw442BY5IKaVyI6selIj8AzgW+APOBIkG4Alguoj8LctjnQx4ReR2EYmJyCLgTeCcgTY2xnwWqAZez3L/JaEz1EMslZz85R4CPm+BI1JKqdzIeoq4iKwGrjyIY80EVmW0vQUcmbmhMWYicB0wG3jmII454uzo6C9xVF+jCxQqpUauvSYoY8zLwKkistMYsxyw97ZtlutBVQKZNfvCwG7rlBtj3DjLelwpItuM0SGsdO2d6QnKV8BIlFIqtwbrQT0JdKce/24IjhUCMr/yB4CujLbvAiIijwzBMUecHR2RvseaoJRSI9leE5SIXJ/2dBnwNxGJpW9jjCnHuVk3GyuBKzLapgO/yGg7FxhjjDkz9bwK+KkxZpaI/GeWxxqx2tOG+EZVa4JSSo1c2V6DWgY0A20Z7ZOAhezZM9rbPqzUkh0LcGbxHQU8mr6RiExPf26MeRW4XUQeyDLWESuRtNnZ1d33XBOUUmokG+wa1GVAby/KAlYaYzKvQ1Xi1ObbJxHpMcZ8BGfq+HzgbWCeiLQZY64GLhCRw/cz/pLS0dVNMun8CSr9Xsq87n28Qymlhq/BelD34Fw3cgH/C9zA7qvq2jjXj57L9mAi8gZw4gDtNwM37+U9x2S7/5EufYLEKL3+pJQa4Qa7BhUndX3IGLMBeBGoEpH2VNt/ACtS26k82Nmp15+UUqUj21JHbYAA305rewJ43RgzecijUgNq79TrT0qp0pFtgroDeJ7+a1IAE4G/40x4UHmwM9jfg6qr0gSllBrZsk1Qs4D5ItJ3z5KIRIAbgffmIjC1u2TSZlewvwdVV1VewGiUUir3sk1Q7cARA7RPA3TZ9zwIhnuIJ/pr8PnKdSFjpdTIlu1Z7l7gbmPMeOAfODP4jsWp+nBfjmJTaXbvPenwnlJq5Ms2Qd2U2va7QGOqrRX4MfCDHMSlMuyWoKp1eE8pNfJllaBEJIlTXfw6Y0wD0DPIGk5qCHWuXkPb0mW8vqmHTm89/rFjqK0cU+iwlFIq57K+kGGMeRfOkhnu1HMLKAeOE5FLcxNeaetcvYZ3Fi7CXREg4m8kGYnRtXYd3sNHAU2FDk8ppXIqqwRljLkGp5JEF1CBU1GiJvXyU7kJTbUtXYa7IoCnooKuqBuX1yIJxP/5EpywxzJaSik1omQ7i+9S4JsiUg1sxSnyOhZn+fflOYqt5EVbWkn2xGhf+Ra7WnbQs3MXVjKJe3tLoUNTSqmcyzZBNQO/TT1+FThBRLYB/w18JheBKXCVl9O1eg3BuAvL7cFOJvHubMXj00kSSqmRb39KHdWnHq8Gjk493gzoFfucsQGbEF7nsQ0VyW4GWdxYKaVGjGwT1BKc+6COwVnX6UJjzPuB/wL+navgSl2yu4eKadOIlAWwEwkst4tRo0eR7O4pdGhKKZVz2c7iuxK4DaeaxK9wFht8DmfSxPm5CU35RjcRCwZJjDmE8m4vANXuXfgadAafUmrky7YHdRZwjYj8UkRsEfkcUAs0iIjO4suRxrlzSITCdEbj2DYkY3H83V00zp1T6NCUUirnsk1QPwEa0htEpEvXgsqt6mlTGX/+uUQ8fpLRKC6vl2lnf4LqaVMLHZpSSuVctgnq78AZuQxE7alz9Rq2PreMrq4obr+fwLixjDlyeqHDUkqpvMj2GlQSuNkY8x1gAxBJf1FEZg11YKWut4pEl68Kl6+JZCyGvfYtutYeoj0opVRJyDZB/T31o/Kkt4pE1FOJ1Q2W10OVx0Xb0mWaoJRSJWGvCcoY40oViUVErt/bdio3oi2tlDXU0xXpH4Wt8nuItmwqYFRKKZU/g12DihljdpvPbIx5nzFGyxjkgW90E4lwmFCi/0/kj4XxjdYp5kqp0jBYgrIGaPsdTg0+lWN9U8wj/VPMfd0hnWKulCoZ2c7i6zVQ0lI50DfF3Ns/xXzqGR/T609KqZKR9XpQKv+qpk7BMz1MXdKpvTf2SFPgiJRSKn/2twel8igUjZNIJSdfmYcyr7vAESmlVP7sqwf1OWNMV8b2nzbGbE/fSER+OuSRKTpD3X2PqyvKChiJUkrl32AJaiNwWUbbNuCijDYb0ASVA8FQf9VyTVBKqVKz1wQlIoflMQ41gM60BFWlCUopVWL0GlQRC4a1B6WUKl15ncVnjDkauAs4ClgPfF5Elg+w3SnArcBUoBX4gYj8PJ+xFoPOUKzvcXVAE5RSqrTkrQdljCnDWZl3Mc5aUjcBzxpjqjO2Gw/8Frgxtd15wC3GmFPzFWuxSO9B6RCfUqrU5LMHdTLgFZHbU88XGWO+DJwD3JO23WHAQhF5NPV8uTHmT8B7gd/nJ9TCs22brvQEFfAWMBqllMq/fCaomcCqjLa3gCPTG0TkeeD53ufGmFHAScBDuQ6wGHSuXkPb0mXs3LadndZE/GPHUNVYj9ej90AppUpLPidJVALhjLYwENjbG4wxNcDjOEt9LMldaMWhdw2oWDBIvLaBZCxG19p1lEU6Cx2aUkrlXT4TVAjwZ7QFgK4BtsUYMw14CWgBzu5d+mMk610DylNRQTjpxuX1YHk88Pa6QoemlFJ5l88EtRLILCY3PdW+G2PM+3B6TY/hJKdo7sMrvGhLK+6A06EMJ526vJbHg6drVyHDUkqpgsjnNahlgGWMuQJYAJyFM9380fSNjDGTcZb1uEZE7shjfAXnG91ELBh0elCpdaDseJzaUVUFjkwppfIvbz0oEekBPoKTmNqBa4B5ItJmjLnaGPNmatPLgSqcqeVdaT//k69YC6V3Dah4KEQ4YZGMxbHjcca+59hCh6aUUnmX1xt1ReQN4MQB2m8Gbk49/i/gv/IZV7HoXQOqbekyOjdGcZVX4z/sUJpnTCl0aEoplXe6HlSRqZ42leppUwkseQNXTxyASr/eA6WUKj1ai68IxeJJoqnk5HJZBHz6PUIpVXo0QRWhrkh/BYlKvxfLsgoYjVJKFYYmqCLUFe4vEqvDe0qpUqUJqgiFIv0JqsKvRWKVUqVJL24Uid4afNGWVtZXjiNWMxFvTQ2VWiRWKVWitAdVBNJr8JU11BMMddO1dh2xjg4d4lNKlSxNUEUgvQafZVn0eP1YHg+RzVs0QSmlSpYmqCKQXoMPIJx0YXk8JCIRKnUlXaVUidIEVQR8o5tIhPtXIokkLOx4HLffT4X2oJRSJUoTVBFIr8EXT9pEehLY8Tj+sWMJlOs8FqVUadIEVQR6a/B5q6robNuFy+ulcspk6prrcbn0Jl2lVGnSr+dFoq8G3/YQ1cvWAOjwnlKqpGmCKrD0+598o5sIHnl832uaoJRSpUwTVAH13v/krghgeTzsXPEqa19ew65DjqBi4mFUTG4odIhKKVUweg2qgHrvf7JjcULr1gPQ46sgFgzStXYdVntrgSNUSqnC0QRVQL33P0W2bMHyenB5PHS7yyGZxPJ46Hn9tUKHqJRSBaMJqoB6739KRCJYHme0NZq0sDweLI8H187tBY5QKaUKRxNUAfXe/2R5PCRjMZKxGFHbjTsQwI7HqWmqK3SISilVMJqgCqj3/qeKiRNJhiNgWcRrG7BcLux4nPFzTip0iEopVTA6i68AMqeWjz1zHmPPnMfW55YR/7eNu9xLxcTDaJw5rdChKqVUwWgPKs8yl9aIBYO8s3ARAM2f+Sx1xx1L9cwZ1DY36FLvSqmSpj2oPHv7wV8QWrcBOx7HVVaG/9AJlI+qo23pMvxnn9+3XcCnfxqlVGnTs2Aebf3jcwRXCZbH7UyMiMcIrVmLPWUKiXAEOxLv21arSCilSp0O8eXR1iVPOFPIXW4sl4XL4wGXReTtt/GNbiLcHevbVquYK6VKnSaoPIrt2oWnqgqSSeyk7TRaFsmeHhrnziGU1oMKaA9KKVXi9Gt6jqXP2LOTSYjF8NZUEw+FsRNxsKG8oYHqaVOJ/GNj3/u0B6WUKnV6FsyhrX98jk3/txg7nsBdWYGnpobubS1AJZ7aGpLd3djdPYz91FkAu/Wg9BqUUqrUaYIaYlv/+BxblzxB9/YdJLujuH1+vHW12PE4dixOefNo4h0dJIJBvLW1HHLeORzywQ8A7H4NyqcJSilV2jRBHaTO1WvY/MhjdK1ZSzwSIRmN4qmsdF5M2iTCYVxeL57KCpKAZVlUTZvG4dd9Z499RaJp16B0mrlSqsTpWTBLTiJ6lI6Vb2F3d+PyleMfN454ZyfxziBWeRnJaBSSSRLhMFhgud3YiQTxri48lRVYHg+JrhC+I5v22L9t24TTEpRfr0EppUpcXs+CxpijgbuAo4D1wOdFZPmBbjfUeic0dK3fQCISIWmDywIsi9iuDhLRKNhJsCwSoTBdbwm2ZeGtqsLt9fa9ZieTWJaF7bIgaWEnEti2TaK7G8vjpnHunD2OHe1JkLSdmX3lXjcet06wVEqVtrydBY0xZcASYDFQC9wEPGuMqT6Q7YZabwmiyOYtdLdtp6ejg+5tW4l1dBLdto1EOIwdiwHO/UuWx+PMyovHSXZ3A2C53M7ObCezuSwXWBaW200iHMYCxp13DtXTpu5x/HC0//qTX4f3lFIqrz2okwGviNyeer7IGPNl4BzgngPYbkj1rm4b/vdGXGVeksHuvqrilg3JRAIAO5EArwfL1V8nz0lc4K6sIL6rA1wWlseLy+vFDocpP6SZ2sNn0jh3zoDJCdhteC9QrhMklFIqnwlqJrAqo+0t4MgD3G5IRVtaKWuoJxGJ4PL5nMTkcmMn4lheD8TjYAHJJIBzo63H49x0i00iFsNVXo7l80HSSWbemmomXHhB3yy9wZSXufsej6rx5eQzKqXUcJLPBFUJhDPawkDgALcbUr7RTcSCQdx+P8lYrK9WnsvjxV1eTjIWw04kwYJkPA5JG3d5Gd66Ospqa4hubQGg7l1HM/bMM/baU9qbxlo/s48cw66uKO+eMToXH1EppYaVfCaoEODPaAsAXQe43ZBqnDuHdxYuoqyujvCmzeBxY/d0O9eagMCECXS3bScRj2MlErgqyqmeMYOxZ87b72Q0EMuyOHb6nrP7lFKqVOUzQa0Ershomw784gC3G1K9q9u2LV1GorvbmcVXXYPLArffT+WkiUy69AtDkoyUUkrtWz4T1DLAMsZcASwAzsKZRv7oAW435KqnTdUEpJRSRSJv08xFpAf4CE7CaQeuAeaJSJsx5mpjzJv72i5fsSqllCq8vN5wIyJvACcO0H4zcPO+tlNKKVU6tFyBUkqpoqQJSimlVFEaiTV13ADbtm0rdBxKKaVS0s7J7sG2SzcSE9QhABdccEGh41BKKbWnQ4B12Ww4EhPUcuAkYCuQKHAsSimlHG6c5JT1yhSWnVriQSmllComOklCKaVUUdIEpZRSqihpglJKKVWUNEEppZQqSpqglFJKFSVNUEoppYqSJiillFJFSROUUkqpojQSK0lkxRhzNHAXzmKI64HPi0jWdziPdMaYU4BbgalAK/ADEfl5YaMqPsaYWuBfwLUi8kCBwykaxphDgJ8Bc4AocLeIfLewURUPY8zxwE8AA7QBt4rIvYWNqviUZA/KGFMGLAEWA7XATcCzxpjqggZWJIwx44HfAjfi/H7OA24xxpxa0MCK013A2EIHUYSW4JQbGw0cD3zWGHN+YUMqDsYYF87v5yciUoPz/9eC1JdmlaZUe1AnA14RuT31fJEx5svAOcA9BYuqeBwGLBSRR1PPlxtj/gS8F/h9oYIqNsaYzwLVwOuFjqWYGGPeA0wC3isiMWCDMeZkIFLQwIpHHdAEWMYYC7CBONBT0KiKUEn2oICZwKqMtreAIwsQS9ERkedF5Eu9z40xo3AK8K4oXFTFxRgzEbgO+HyhYylCx+Ek7e8ZYzYbY9YBZ4jI1gLHVRREZAewAHgQiOEUT71aRDLPSSWvVBNUJRDOaAsDgQLEUtSMMTXA48DfcYYlSp4xxg38ErhSRHThsT31fqGJ4fSkzgSu1CE+R2qILwqcD/hxRnSuM8Z8qJBxFaNSHeIL4fyHkS4AdBUglqJljJmGk5RWAheISLLAIRWL7wIiIo8UOpAi1Q10isj3Us9fM8bci5OoFhYsquJxJs7w5zdTz/9sjLkPuBR4tnBhFZ9S7UGtxJk9k256ql0Bxpj34fSaHgPOFpFogUMqJucCZxtjdhljduEMDf/UGPPTAsdVLN4CAqnJSL1K9cvwQMYD5RltcZwep0pTqv/RLMO5QHkFzljwWTjTzR8d9F0lwhgzGfgdcI2I3FHoeIqNiExPf26MeRW4XaeZ9/kDztTpHxljvoHzZfBi4LKCRlU8nsWZFftFnElZxwJfAC4paFRFqCR7UCLSA3wEJzG1A9cA80SkraCBFY/LgSqc/4m60n7+p9CBqeKX6m2/H+f601bgGeD7IvLbggZWJETkTZxhvkuBXTjDnleJiF7jzaAr6iqllCpKJdmDUkopVfw0QSmllCpKmqCUUkoVJU1QSimlipImKKWUUkVJE5RSSqmiVKo36qphzBhTAVyFU31+PLANeAS4QUR2DdExKoFzROS+IdiXF/iCiBRFpQljzOeAH4pIQx6OdQVQKyLX7WO7K4EKEbk+1zGp4UN7UGpYSa3Z9RLwQeCrOJXpvwScirOmV2aNxQP1DYau8sH5wPwh2tewkVpX7GvAD7LY/A7gwlT9R6UATVBq+LkV57/buSLyjIhsEJHf41QGORq4aIiOYw3RfoZ6X8PJfwMPi8g+izCLSDfO8hNX5zwqNWxoJQk1bBhjynGWn/+WiNw1wOuzcaqM70gtBPcV4MvABJwCpleLyFOpbR/AWWLFD5wN7ADuFZEbU0Ng96fteiJObbkfAPOARpxhxZ+LyI2p/blxTq5fwFluYnnq2I04tR97zRGRP2XE/T2cWpDrcdaX6gZ+A1whIolUrJUicnbae97GGaZbkHr/kcAanF5fD3AD8ArOsuuTgL/gVKRv7x3iA27D6SmS2u46EUmk9v+e1OvHAe/g1Iz7oYgkU+//Kk5x5dNwyhjdmPGZKnHKHJ0qIi+m2vzAzTgryPqBPwH/KSKbU68fAfwTGKdlxxRoD0oNL5NwVrBdPtCLIvJiajE4cJLF9cC1OCf/x4DHM5bVvgTYBLwbuBe4wRhzHLAY+BHwGnAIzgn6NuAEnARlgJ+kbU/qOF8Fvg68K/WeJ4G/pdraU/t6cS+f7eNAbeoY1+LUQ5yXxe8k/f2VOIVH70rFuwD4T+BDOInm62nb1+MMi84FLsSpC3cVgDGmCWfl5GdwEt9XU/H8d9r73wUEU8f7xQDxvB9I4Hz+XncBpwOfxlkG3ofzuwZARN7A+SLwgf343GoE0wSlhpO61L8dg22U6j19HbhJRBaJyOrU2kR/YPeT7DoR+a44bsRJIseJSARnbbC4iGxL9Sr+ClwsIi+LyHoR+UFqm5mp410G3Cgij4jIGpwT+m9xiu52AHZqX3tb1juC05sQEbkHJzket5dtBxIFvi4i63Cu57iBBanVkf8KPAUcnrZ9Avi0iLwmIk8CN+IkM1KxLxeRG0RkjYg8DXwbuDLjmNeLyFoR2ThAPO8GVomIDX0LX56P0yv8Y2r12C8BL2Qsy7Ey9V6ldBafGla2p/6tG3QraAIa2P3bO8ALwCfTnq/JeD0IePeyz18CHzfGfAaYBhyD02Nxp47VSFrPTkSCpIbPjMlcemxA/85IXp2DxDKQt0Uknnrcu1r0+rTXo6k4e20WkXfSni8HxhhjanES2cnGmPRrRy7Ab4yp7z3GPlYTHk3/3wuc35mH3X9HG0j12tLswPn7KaUJSg0ra3FOYP/BAMN8xpjbcK57/Hwv77fYfdRgoN7M3iY0/C/OzMFfpH4uw+nlpO/nYC7oDhbLQPvN/H93oMXuBlsBOZHxvPf3Ekvt+7fAdwZ4X2/vNTLIvnuPnf67zPZ35MZJpkrpEJ8aPlJLzj8EfMUY40t/zRhzGM51lG4R6QS24FzPSTcbZ7JENvpOpMaYKuAzwIUicrWILMY54dYAloh04EzeODbtPX5jTEtqssHBzkTqPVbvvis5+F7GOGNMek/0eJxeXAhYBUxPDd+tFZG1wAyca2ODJb1023B6lb3W4yTF9N/RocaYdmPM2LTtGlLvVUp7UGrYuQH4GLDMGHMdzjDd0cD3gX8Bd6e2uxWYb4zZhDMz7BycSQEnZ3mcLqDZGDMJZ8JDCDjTGLMBGIMzo8+if+nuHwPfTc2uW42zCGYHsAIYB1QaY2YC61ML+u2P5cDFxph5OMnjevbsAe0vN7DQGHMVzizFq+mf4n0n8FVjzE9Sj8fj9EqXpGbxZbP/fwJXGWPcIpIQkaAx5l7gNmNMEGf470fAv3pn8aUchTO5QyntQanhRUTagROBf+CcNFfinOgeAz6cdvJfgJO0vg+8jjN77OMi8nyWh3oY55rUSpzrTefjDPGtxBniexZn4kHvRIYfAPfhzAZcgTNj72Op60rPpdpW4CTX/fUQ8EDquM8Dr+JM2jgYq3CSyF9wkvqPSQ2NisgmnGT+bpxhzIdwZtt9fcA9DWwZTgJPn/DwjVT7ozizGYM4XxwAMMYcjjOp5A8H8oHUyKP3QSmlcsIYcyeQEJGvZrn9LUCTiFyc28jUcKE9KKVUrnwfOCs1xXxQqZt4z0+9RylAE5RSKkdE5N84Q4ffzGLzLwMPiojkNio1nOgQn1JKqaKkPSillFJFSROUUkqpoqQJSimlVFHSBKWUUqooaYJSSilVlP4/h2zGxT3KpycAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_sweep_frame(frame)\n",
    "plot(frac_infected_series, label='Analysis')\n",
    "\n",
    "decorate(xlabel='Contact number (c)',\n",
    "         ylabel='Fraction infected')\n",
    "\n",
    "savefig('figs/chap06-fig04.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The agreement is generally good, except for values of `c` less than 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  If we didn't know about contact numbers, we might have explored other possibilities, like the difference between `beta` and `gamma`, rather than their ratio.\n",
    "\n",
    "Write a version of `plot_sweep_frame`, called `plot_sweep_frame_difference`, that plots the fraction infected versus the difference `beta-gamma`.\n",
    "\n",
    "What do the results look like, and what does that imply? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotSweepFrameDiff(frame):\n",
    "    for gamma in frame.columns:\n",
    "        series = frame[gamma]\n",
    "        for beta in series.index:\n",
    "            frac_infected = series[beta]\n",
    "            plot(beta-gamma, frac_infected, 'ro')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap06-fig04.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VPW5+PHPbFlmkhAISYAAshS+gLtUXGot0MVbrxep9l4RW1urrfXazXvtvS6trbh1ubW2Ymtdftp69Qe/27q1LrW30NbWWmnrhuATZJGdhARIMksyy/n9cSbDZAhhgMyck8nzfr3yYubMmTNPApwn3+95zvP1WJaFUkop5TZepwNQSiml+qMJSimllCtpglJKKeVKmqCUUkq5kiYopZRSrqQJSimllCtpglJKKeVKmqCUUkq5kiYopZRSrqQJSimllCv5nQ5gsBljyoFTgR1A0uFwlFJK2XzAWGCViHTn84aSS1DYyelFp4NQSinVr/cDf8xnx1JMUDsAHn30UcaMGeN0LEoppYCdO3dyySWXQPocnY9STFBJgDFjxjB+/HinY1FKKdVX3pdetEhCKaWUK2mCUkop5UqaoJRSSrmSJiillFKu5EiRhDFmDvArEWk4yOsTgQeB04EW4Isi8mwRQ1RKKeWwoo6gjDEeY8wVwAtA2QC7LgPeAOqAzwLLjDFTihCiUkoNe5ZlkUpZJJMp4okUPfEkyWSq6HEUewR1M/CPwK3A1/rbwRgzHXgv8GER6QFWGGOeBi4HbixWoEqp0mJZFinL/jPzOGWRsiwsa//rqZSFhb0tlUr/mX5P7+NUygLs11NW//tYqX629Ymhd7/92/rEkrKwyHrc+9n2xvRnH+T7yj7WwY7f7377PydXecDH3NnjmTZhZNH+zoqdoO4VkZuMMXMH2GcWsFlEwlnb3gbmFDQyF+poXkfripXEdrXgLS8DPKS6uzOPu9vaScWi+CorKRs1MvN6RWMD9fPnUTN9msPfgSoWy7JIJFMkU7knz/5OfP2drLJPogeeFK2c1zPbUpC0LMicPLMep0+w9nFyTtKpPE6sB421/5P5ASf+rO+h90Svjlx3PEnzu3tKN0GJyPY8dqsCIjnbIkBw8CNyp47mdWx7/Ak6Vq/BFwoSGDmK7g0bAA/lY8fSvWEDqXgCj8+HpyxAz569RLdtx+v3UzV9GvHOTrY8towJixdpkhri4okkkVgi/RUnEksQ7U4QTj+OxOJEu+3XEw5MwajS5fF48HjAm/4zVBnghGn1RY3BjZ0kwkBlzrYg0OVALEXTO1rqXL+ReHsbWBbeoP1jiGzYgC8UxFtWRnTzZvxVIZLdPVipFOU11XSH7XzurQoR27GDmlkzAWhdsVITlAslkyki3X2TTubP7gSRaDz9epx4QpPOYPJmn3S92SdgD15P7kk5vc3rybzP3pbeJ/1+Dx68Xntf0q/3Hsvr7T1m1rY+x+rv+P18Dh483r7xe70e+3tKvx/A5/Uc+D14D/J95fws9m/b/9hpbkxQa4CJxphKEYmmt81Iby9JHc3r2PLYMnyhIKlIBAtIdHQSqB2Bt7wcy7JIxrrxBYOkenrw+EdAKoU9gULmscfvJxm1f2S+YJDYrhbnvilFdzzJno4Yezq6ae+MsbcjRltHjM5IvCDTTX6ft58TVNYJLOvE5qH3BJvzuHef9Imz7wnzICdbb/8n1t7HXi85J9is9/UTR5/35ySRPt9X9jZv/yfe3u/Vk3VMNXS4LkGJiBhjXgduM8ZcD5wJnA+c4WxkhdO6YiW+UBB/KEQyFsVXXk7S7yfRFcZXXo43ECAV78FKJPCWlWElEuD1kvmvln5sJRL4Ku1RVzISoaKx3yp+NYgsyyLanWBPZzftHTH2dMRo7+hmb2eMrmj8qI/v83oIVgQIVvgJlvsJVgaoLPcTqghQWeFPb7dfD/i9egJWJcUVCcoYcwnwExGpSm+6ELgP+x6o3cDlIrLaqfgKLbarhbLRdQD4KitJ9fTgqwqR2LuPVDwOZQGIx0nFuqmcOJHuHdvxeL14fD4SsRgev9+u6ol1E5w4kUQ4TDIcYdz5Cxz+zkpPKmWxe2+UHW1hdraF2bE7fNiJyOPxEKrw2wkmnVyCFYF0Atr/uLLCT3nAp0lHDVuOJCgR+R1Qm/X8UeDRrOdbgI8WP7Li6r3uFNm8hej27YSOOYaKcWMJv7MBPOCvHQEeD55Ekqrp0whUV5Pq7qZs5gz6VPHV1fWp4gvUjWLc+Qv0+tMg6Ikn2dkWZmdbhO27w+xqD+d9Xcjn9VBbVc7ImgpG1VQwsqackdUV1FaX4/dpExelDsUVI6jhKPu6U3DKZMLNzXRKM1XTp1HRNI7Y9u2UjRpF1ZTJmZLx7LLzisYGmi5YqEmoAPZ0xti4vYNN2zvY0RY+5PWigN9rJ6DqrGRUXUFNqCxzIVspdfg0QTkk+7qTPwQeM53Iu5sJb9jIqFNnc8wnL+mTfLITWtnoOi0lH0SplMXOtjAbd3Swcfs+9nYOvBp1VWWAsaOrGDc6xJi6EHUjKjQRKVUAmqAckn3dCaCstpbAiBH07G5j6uc/d8D+2QkNyPyppeRHrm1flLWb2pF39xDtTvS7j8fjYfSICsbUhRg7OsS40SGqggN16VJKDRZNUA6paGwg3tmZSTQwcOVdbkIDLSU/ErHuBM1b9rB2Uzute6L97hPweZkwpprJY0dwzNhqghWBIkeplAJNUI6pnz+PLY8tA+xEk4xEBqy8O9yEpvpq2xflteZWmjfvIZk68JpSqCLA5HE1TB43gqaGKi1iUMoFNEE5pGb6NCYsXtSn6OFglXcdzevo2bs30/ooOHEi3rKAlpIfgmVZbG3p4tXmFjbv7DzgdZ/Xw5SmWmZOGsn4hmq9jqSUy2iCclDN9GmHvH6UXRxRNXMG0c2b6Vz7NjXHHasFEgPYsquTl1fvYFd7bltHaBwVZOakUbxnQi0VZfpfQCm30v+dLpdb7Vc+spZEOEygulqTUz92tUf485s72NrSd8Tk8XiY0jSCk6fXM6YudJB3K6XcRBOUy2lxRH66onH+9Po21m3Z22e73+dl5qRRnDitntrqcoeiU0odCU1QLqfFEQNLpSzefGc3f1mzk554MrPd6/Ewa/Io3jtrDFWVWoWn1FCkCcoBuR0hBlpc8HCr/YaT3Xuj/HbVZlr39i0XnzahljnHjmFkdYVDkSmlBoMmqCLLtyNEn9V0K8qxkil6drcNWO03XKRSFq81t/KXt3b0KRkfVVPBB04ZT1N91QDvVkoNFZqgiiyfjhC5Sax31KRVe/a1phdefpftu/evX+n3eTl1ViMnTavHp/cvKVUyNEEVWT5FD9rWqH/bWrv49cvvEontX96iYWSQD8+ZyMganc5TqtRogiqyfIoetHKvL8uyeH1dK396Y0ems7jH4+HUmY3MntmIT2+wVaok6XxIkdXPn0cyHCERtpdx6F1csH7+vMw+FY0NJCN9bzAdrpV7qZTF7/++lT++vj2TnCrL/Zx/9hTmHDtGk5NSJUwTVJH1tjgKVFfTs7uNQHX1AdeW8kliw0E8keSZP21k9Ya2zLYxdSEu+rBhfEO1g5EppYpBp/gccKgWR4fTp69UxXoSPPPHjexoC2e2TZswkg+eOkEbuSo1TGiCcql8+vSVqlh3gqf+sL7P/U3vndnIaceOwePRKT2lhgtNUMpVYj0Jnnqxb3J6/4lNnDi93sGolFJO0ATlIofTYaIUxRNJfvXHjZmFBD0eD/Nmj2fW5LpDvFMpVYp0Mt8lem/OjXd29ukw0dG8zunQiiKZsnjuz5vYmXXNae4pmpyUGs40QblE9s25Ho8HfyiELxSkdcVKp0MrOMuy+N3ftvRZVPADJ4/n2CmanJQazjRBuURsVwu+YLDPtuFyc+6r0sraTe2Z53NmjeH494x2MCKllBtognKJ4Xpz7sbt+/jz6h2Z57Mmj+LUWY0ORqSUcgtNUC4xHG/O3dvZzf++sjnTIWLc6Co+cPJ4LSVXSgGaoFwjnw4TpSSRTPH8y5voTi8yWBMq46NnTtJu5EqpDC0zd5HhdHPun17fzu70vU4+r4d/OH0SleX6z1EptZ/+uqqKbuP2fby5fnfm+VknNdEwKjjAO5RSw5EmKFVU0e4EK/66JfN8atMIjtNycqVUP4o6p2KMORG4FzgB2AB8RkRW9bPf6cAPAQO0At8SkQeKGasqjN//fSvR7gQAoYoA82ZP0KIIpVS/ijaCMsaUAU8By4Fa4DbgBWNMTc5+3vR+PxSREcDFwNJ0clND2IZt+3hn697M8/nvnUCFXndSSh1EMaf45gIBEblLROIisgx4C7goZ7+RQAPgMcZ4AAtIAD1FjFUNsu54kj+8ujXzfNbkURwztmaAdyilhrti/vo6C1ibs+1t4PjsDSLSZoxZCvwUeAjwAV8Wkdz3lozh0CT2lbd20hWNAxCsCHDmCeMcjkgp5XbFHEFVAZGcbRGgT/lWeoovBiwGKrFHXt8wxnykCDEW3XBoEtu2L8qb72RV7Z04jooyndpTSg2smAkqjJ1wsgWBrpxtFwDvE5Fl6anA3wMPAlcWIcaiK/UmsZZl8eJr20ilu0WMb6hm2oRah6NSSg0FxUxQa7Cr8rLNSG/PNgEoz9mWAOIFistRpd4kdtOODra22L+DeD0ezj65Sav2lFJ5KeY8y0rswodrgKXAhdjl5k/k7PcCcIcx5nPA/cApwGeBK4oYa9FUNDYQ7+zEHwpltpVKk9hUyuJPb2zPPD92Sh2jaiocjEgpNZQUbQQlIj3AR7ETUztwI7BQRFqNMTcYY95K7/cW9jTflcBe4DHgOhF5qlixFlMpN4ldu6mdvZ3dAJQFfMw5dozDESmlhpKiXqkWkdXAWf1svx24Pev5s8CzRQzNMb1NYrOr+Madv2DIV/ElkilWrdmZeX6KadBee0qpw6JnDBfIbRLb0byO9ffeN6TLzt/a0NanrPzEafUOR6SUGmq0F5/LlELZeSKZ4u9v7y/ymD2jgYBf/6kppQ6PnjVcphTKztduaiccs0dPoYoAx2ozWKXUETjoFJ8x5ux8DyIifxiccFRsVwtlo/ue0IdS2XkqZfGq7I/1FNOAXxchVEodgYGuQf0u57kFeIAUkAQC6cc95HSDUEduqJedv7N1Lx1hu21iRZmfWVNGORyRUmqoGuhX2+qsr8uAN4EzgHIRKQdOAlYB1xQ6yOFkKJedW5bFq837R08nTBtNwO9zMCKl1FB20AQlIuHeL+AW4LMi8hcRSaZffwO4GlhSnFCHh96y80B1NT272whUVzNh8aIhUcW3fXeY1j32Mu5+n5fjp452OCKl1FCWb5l5DXZX8VzVh3EMlafcsvOh4o11rZnH5piRet+TUuqo5HsG+TnwULpN0avY16JOA+4Eflag2NQQ0hnpYcP2jsxzve9JKXW08k1QXwR+jL3Sbe974sADwFcLEFfJKfU1n97a0IaV1bFce+4ppY5WXvW/IhIVkU8Do7FHTqcBo0XkCyKSKGB8JaEUbr4dSDJlsWZje+b58VP1viel1NHL+wYVY0w18Mn01ybgbGPM5ALFVVJK4ebbgWzcvo9I+sbcqsoAk8eNcDgipVQpyCtBGWOmAwJ8CbgKu2jiYuB1Y8yZhQuvNJT6mk9rNrZlHs+YNAqvV9d7UkodvXxHUD8A/p+IzAC6AUTkk8AjwHcLFFvJqGhsIBnpu9p9bOcuunfv5q2bb2X9vfcN2em+rkgPW3btXxR55iS9MVcpNTjyLZI4A/hyP9vvBF4fvHBKU/38eWx5bBmQHjnt3EVk4yZCkydlrkltuP9ByuvqSHV3D6kiirff3dOnOGJEVe5iyEopdWTyHUFFgf567UwDOvrZrrLk3nwb37uX0ORJVIwdg8fjIRWP09PSQnjjxiFVRGFZFvLunszzGZNGOhiNUqrU5DuC+ilwjzHm8+nn9caY47Cn/h4tSGQlJvvm27duvrVPQ9jY9h14ysuxEolMEQXYxRVuHkW17omypzMGQMDvZWqTFkcopQZPvgnqa9jNYlcA5cCfgQTwI+yl29VhyG0Im4xGwefDV1mZ2WcoFFE0b9k/epraNEL77imlBlW+90ElROR6YCRwPHAyMFJEvgLUFjC+kpTbENbj92N191A5blxmH7d3MLcsi3e27M08nzZRp/eUUoMrrxGUMSYJjBGRVuCtrO0TgTVAVWHCK02916R6O0uEJk+ip60dT8CPZVkkIxGS4Qjjzl/gdKgHtWN3OLOke2W5nwkN1Q5HpJQqNQMtWHgx8LH0Uw/wgDGmO2e3Y4B21GHLbQib2wpp3PkLXH396Z2t+0dPU5tG6L1PSqlBN9AI6jfAh7GTE9iVfNGs1y3gL8DDBYlsmBlKHcwty2L91n2Z5++ZoLO8SqnBd9AEJSK7gc8AGGM2Ad8VkczdpsYYv/bhG552tkUIx/ZP740brTO8SqnBl+99UN8H7jXG3JC1bYMx5kFjjC73Psxs2LZ/9DR5nE7vKaUKI98EdTdwLPB81rZPAicA3xvsoJR7WZbFhu37E5Te+6SUKpR8E9R5wGUi8vfeDSLye+BK4MJCBKbcaU9nN/u67FqZgN/L+Aad3lNKFUbey20AB1uBrmwwAlFDw8as0dPEMTX4fIfzT0gppfKX79nll9itjo7r3WCMmYk99fdsIQJT7rQpa1n3yeNqHIxEKVXq8m11dA3wJPBG+l4oC7vl0QvYy8GrYSDWnWBnu13I6fF4mNioN+cqpQonrwQlInuADxhjZgGzgB6gWUTeLmRwyl027+rMLK3RMLKSYEXA4YiUUqUs3xEUxhgvMAmYADwETDfG1IhI3sttGGNOBO7Frv7bAHxGRFb1s99Y4MfAPCAG3CciX8/3c1RhbN7ZmXl8zFid3lNKFVa+S76PBV4DlmGvoDsKuA5Ym74Wlc8xyoCngOXYDWZvA14wxvR3pnsK2AE0AqcDnzLGLM7nc1RhWJbF5l1ZCWqMJiilVGEdzpLva4B69rc7+gTwcvq1fMwFAiJyl4jERWQZduPZi7J3MsacBkwBviQiMRHZmH7vyjw/RxVAe0eMSLp7REWZn/raykO8Qymljk6+U3zzgbNFpNsYA4CIRIwxN2L348vHLGBtzra3sZfvyDYbeBP4pjHm09hTfD8SkZK7ITi3Qaybl3nPnt6b0Fil3SOUUgWX7wjKQ//3QY3GLpjIRxUQydkWAXJbJY0C3g/EsUdSFwDXltoUX0fzOrY8tox4Z+eQWOZ9a0tX5vF4XVpDKVUE+SaoXwDfS1+LsgCMMScA92BfL8pHGMidFwoCXTnbuoEOEfmmiHSLyOvAA9iJqmS0rliJLxTEHwpllnn3hYK0rnDfTGYymWL77uwEpd0jlFKFl2+C+jdgJ7ANeyS0BngV2Jh+LR9rAJOzbUZ6e7a3gWC6qKJX3tWGQ0VsVwu+YN/Bo1uXeW/ZEyWeSAFQEypjRFW5wxEppYaDgRYs/BzwmIh0iUgXcHH6mtPM9PvWikjzYXzWSsBjjLkGWIrdw+8E4Imc/X4DtGKP2P4dO6ldDlx1GJ/lehWNDcQ7O/GHQpltbl3mfWvL/utPOnpSShXLQCOou4CRYC/5boypF5ENIvKMiDx1mMkJEekBPoqdmNqBG4GFItJqjLnBGPNWer8Y8AHs6087sDuof0dEfnG435yb1c+fRzIcIREOY1kWiXCYZDhC/fx5Tod2gG2t4czjpnpNUEqp4hho6mw79hpQr2AXSXzVGJN7vQgAEVmSz4eJyGrgrH623w7cnvV8A/CP+RxzqKqZPo0Jixe5fpn3ZDLFzjZNUEqp4hsoQX0auAk7UVjYy7/3t4KuBeSVoFRfQ2GZ95Y9URLJ/defqoLavF4pVRwDLfn+R+AjAMaYjcCH08vAq2Fkx24dPSmlnJFvs9jJAMaYQPo9npzXc+9vUiUiu7x87OjQAHsqpdTgyitBpdsP3Qccl/OSB3uKzzfIcSkXsCyLHVnXn8aN1hGUUqp48r2/6AfAPmAhkHf3cjW0tXfE6O5JAlBZ7mdElV5/UkoVT74J6njgdBF5s5DBKHfZ2bZ/5nZMnd3xQimliiXfBLUWGIfdxFUdoaHUHBboU14+tk6vPymliivfBHU3cL8x5m5gHTkNYkXk2cEOrNT0Nof1hYJ9msNOWLzItUlqV3v2CCq3p69SShVWvgnqofSf3+7nNS2SyEN2c1gg82fripWuTFCxngTtHTEAvB4P9SM1QSmliivfMvN8m8qqg4jtaqFsdF2fbW5tDgvQuieaeVxXW0HAr/8ElFLFNVCz2GDv/U3GmAF/fdb7oA5tKDWHhb7Te406elJKOWCgX4s7jTG9Z88uoLOfr97t6hCGUnNYgF1ZBRKNo7RAQilVfANN8c3H7joO4M6z6BAyVJrDgn2D7q6sKb6GUbnrTCqlVOEN1Ivv9/09VkduKDSHBQjHEkRicQACfi8jqyscjkgpNRzplW91gNY9+68/1dcG8Xr1Bl2lVPFpglIHaNXpPaWUC2iCUgfoO4LSBKWUcoYmKHWA1r37R1B6g65Syin5LrdRD3wNmA0EOHA9qDmDH1rpcnNPvkgsTlc0XSDh81JbVe5wREqp4SrfVkcPAqcB/40ut3FU3N6TL3v0VFdbqQUSSinH5JugzgbO13Lzo+f2nny7sxLU6BFaXq6Uck6+16D2YC9YqI5SbFcLvmDf6zpu6sm3e28s81ivPymlnJTvCOomYKkx5hr6X25De/Hlye09+dr2ZU3x6QhKKeWgfEdQd2Jfg3oZaOPAnnwqT27uyZdIptjT2Q2Ax+PRBKWUclS+I6iPFzSKYaS3J9+2x59g799fA6Bq2nscjsrWvi+GZVkAjAiVEfDrMl9KKefkux7U7wGMMZXANOyR13oR0dHTEUrFuqmeMR1fMEgyEnFFJV/bvv3Xn3T0pJRyWl5TfMYYnzHm29jFEq8CfwdajTH3GmPyHYWptOxKPo/Hgz8UwhcK0rpipaNxtXVkX3/SDhJKKWflew3qNuATwKXAhPTXpcC52AUU6jC4tZKvPWsENUpHUEoph+U7+rkUuEJEns3a9v+MMZ3AfWiSOixureRr78ia4qvRBKWUcla+I6gq4J1+tm8ARg9eOMODGyv5Yj2JTIsjn9fDCG1xpJRyWL4JahVwdT/bvwD8Ld8PM8acaIz5szEmbIx50xhz6iH2rzXGbDbGfDrfzxgKeiv5AtXV9OxuI1Bd7XiBxJ6O7szjkTUV2uJIKeW4fKf4/hP4nTFmLva9UACnA5OAf8jnAMaYMuAp4C7s1kkXAi8YY44RkYP197sXaMozxiHFbavrZk/vjdLpPaWUC+Q1ghKRvwKnAL/BLpAYDfwSmCEif87zs+YCARG5S0TiIrIMeAu4qL+djTGfAmqAN/M8vjoKmqCUUm6Td4m4iDQD1x7FZ80C1uZsexs4PndHY8xk4BvAmcDzR/GZKk/ZCWpktV5/Uko576AJyhjzCnCOiOwxxqwCrIPtm+d6UFVAbs++CNCn3toY48Ne1uNaEdlpjMnj0Opo7clOUDqCUkq5wEAjqGeA3ivnvxqEzwoDuXd/BoGunG1fB0REHh+Ez1R5iCeSmQo+r0cr+JRS7nDQBCUiN2c9XQn8WUTi2fsYY8qxb9bNxxrgmpxtM4Cf5WxbBIwzxlyQfl4N/MgYM0dE/jXPz1KHIbuCb0RVOT6t4FNKuUC+16BWAmOA1pztU4DHOHBkdLBjeNJLdizFruI7AXgieycRmZH93BjzGnCXiDycZ6zqMO3pzJ7e09GTUsodBroGdRXQO4ryAGuMMbnXoaqwe/Mdkoj0GGM+il06vgTYBCwUkVZjzA3AJSJy7GHGrwZB7xIbALU6vaeUcomBRlD3Y1838gL/B7iFvqvqWtjXj36b74eJyGrgrH623w7cfpD3nJTv8dWR2ZuVoEZWa4GEUsodBroGlSB9fcgYsxF4CagWkfb0tlOBV9P7qSFsb1d2FwkdQSml3CHfVketgADXZ237JfCmMWbqoEelisayLPbpFJ9SyoXyTVB3Ay+y/5oUwGTgL9gFD2qICkfjxJMpACrK/FSU6/JeSil3yDdBzQGWiEjmniURiQK3Au8rRGCqOPoUSGgHCaWUi+SboNqB4/rZPh3QZd+HsH1d2dN7ZQ5GopRSfeU7n/MAcJ8xZgLwV+wKvlOwuz48WKDYVBHs6+rJPK7VCj6llIvkm6BuS+/7daA+va0F+D7w3QLEpYoku4KvJqQjKKWUe+SVoEQkhd1d/BvGmNFAzwBrOKmD6GheR+uKlcR2tVDR2ED9/HmOrwnV0aUVfEopd8q7ZMsYczL2khm+9HMPUA7MFpErCxNe6ehoXseWx5bhCwUpG11HvLOTLY8tc3QlXcuy2BfeP8VXo9eglFIukleCMsbciN1JogsIYXeUGJF++dnChFZaWlesxBcK4g+FADJ/tq5Y6ViCCkfjJLJLzMu0xFwp5R75VvFdCXxVRGqAHdhNXpuwl39fVaDYSkrXho2EN73Lnr/9nY41a4nv3YcvGCS2q8WxmDqyRk8jdPSklHKZfBPUGOAX6cevAWeIyE7gP4BPFiKwUtLRvI6etjaSsRjeigpS8Thd76yne+cuKhobHIsru4KvJqTXn5RS7nI4rY7q0o+bgRPTj7cB4wY7qFLTumIlFePG4bHAisfx+HxYHohu3079/HmOxbUvrBV8Sin3yjdBPYV9H9RJ2Os6XWqM+QDwb8C7hQquVMR2tVAxppHQe6bgLSsj1d2Nv6KSwKg6R6v4skdQOsWnlHKbfK+KXwvcid1N4lHsxQZ/i100sbgwoZWOisYG4p2dlNXWUlZbC0AiHCZQXe1oXB06glJKuVi+I6gLgRtF5L9FxBKRTwO1wGgR0Sq+Q6ifP49kOEIiHMayLBLhMMlwxNHpPcgtktBrUEopd8k3Qf0QGJ29QUS6dC2o/NRMn8aExYsIVFfTs7uNQHW1o/c/AcQTSaLd9l+f1+shVBFwLBallOpPvgnqL8DHChlIKXNlB4ms0VN1sAyv1+NgNEopdaB8r0GlgNuNMV8DNgLR7BdFZM5gB1Yq3NhBAqAzEs88rg7q9SccLOvsAAAZBUlEQVSllPvkm6D+kv5Sh8mNHSQAOrNbHGmBhFLKhQ6aoIwx3nSTWETk5oPtpwYW29VC2ei6Ptuc7iAB0BHRBKWUcreBrkHFjTF92hwYY842xmi512GoaGwgGYn02ZaMRBztIAG516C0QEIp5T4DJaj+rpr/CrsHn8qTW0vMs6f4qnUEpZRyoXyr+HppqddhcmOJOUBn9hSfFkkopVxI11cogprp0xxPSNniidT+e6A8HoJ6D5RSyoUOdwSlSkBX1uipKhjQe6CUUq50qBHUp40xXTn7f8IYszt7JxH50aBHpgome3qvqlKn95RS7jRQgtoMXJWzbSdwWc42C9AENYRk36RbE9LpPaWUOx00QYnIpCLGoYooe4ovpCMopZRL6TWoYagrmj2C0gSllHKnolbxGWNOBO4FTgA2AJ8RkVX97Pdh4FvANKAF+K6I/KSYsZay7Cm+qkqd4lNKuVPRRlDGmDLslXmXY68ldRvwgjGmJme/CcAvgFvT+10M3GGMOadYsZa63Co+pZRyo2KOoOYCARG5K/18mTHmC8BFwP1Z+00CHhORJ9LPVxljfge8D/h1cUItXZZl9Zniq9KbdJVSLlXMBDULWJuz7W3g+OwNIvIi8GLvc2PMKOD9wCOFDrBQ3LQeVHdPkkQyBUBZwEd5wOdIHEopdSjFLJKoAiI52yJA8GBvMMaMAJ7GXurjqcKFVji960HFOzv7rAfV0bzOkXj0+pNSaqgoZoIKA5U524JAVz/7YoyZDrwM7AI+3rv0x1CTvR6Ux+PBHwrhCwVpXbHSkXi6otk36WqCUkq5VzET1BrA5Gybkd7ehzHmbOxR05PYySlW+PAKI7arBV+w7yDRyfWgwlnXn0KaoJRSLlbMa1ArAY8x5hpgKXAhdrn5E9k7GWOmYi/rcaOI3F3E+AqiorGBeGdnZiVdcHY9qOwCCV3qXSnlZkUbQYlID/BR7MTUDtwILBSRVmPMDcaYt9K7Xg1UY5eWd2V9fbtYsQ4mt60H1RXREZRSamgo6o26IrIaOKuf7bcDt6cf/xvwb8WMq5B614PKruIbd/4Cx6r4wjEtklBKDQ26HlQRuGk9KL0GpZQaKrQX3zDTpQlKKTVEaIIaRuKJJD3xJAA+r4eKMr1JVynlXpqghpHc0ZPHoyvpKqXcSxPUMBKJJTKPQxU6vaeUcjctkiggN/XgAy2QUEoNLTqCKhC39eCDnCk+HUEppVxOE1SBuK0HH0AkpiMopdTQoQmqQNzWgw8gHN1/DSpYqbO7Sil30wRVIBWNDSQjfVcXcbIHH+Rcg9IpPqWUy2mCKhC39eCDvlN8wQodQSml3E0TVIH09uALVFfTs7uNQHU1ExYvcraKL6YjKKXU0KG/RheQm3rwxRNJ4gl7zUef10O5dpFQSrmcJqgCcd89UFkFEhXaRUIp5X6aoAqg9x4oK5mku72d8IaNtK/6K+MvvoixH/qgIzFFuvX6k1JqaNEzVQG0rliJlUwS3bYdT8CPrypEMhZj6/9dTmjiREdGUtltjoJ6/UmVsI6ODlpaWojH44feWQ2aQCBAQ0MDNTU1g3ZMTVAFENvVQnd7O56AH6/f/hH7ystJRiK0rljpSIKK9unDp3/tqjR1dHSwa9cumpqaqKys1KnsIrEsi2g0yrZt2wAGLUlpFV8BVDQ2kOwK4/HvTwRWIoE/VOXYjbrZFXxB7SKhSlRLSwtNTU0Eg0FNTkXk8XgIBoM0NTXR0jJ45zhNUAVQP38eHr+PZCyGZVmk4nGsRJLAqJGO3ajbZ4qvXEdQqjTF43EqKyudDmPYqqysHNSpVU1QBVAzfRrjL74ID3b3CG+gjIqmcXh9Psdu1I32uUlXR1CqdOnIyTmD/bPXBFUgoYkTqZo2DY/PTzIWwx8KOnqjbrhPkYSOoJRS7qcJqgB6y8w9Pi+1p5xE9YzppGLdjsaU3eaoUqf4lHKVm266iRkzZrBu3eAsx3Pddddx2223HfH7n376aRYtWjQosRwNTVCDqKN5HevvvY+1t91Bx9vC3jdX07n2bVLxuKNLbViWRbRbR1BKuVE4HOa5557jYx/7GI888ojT4QCwYMECli1b5nQYmqAGS++oKbptO4muMFYqSSrWTbyri/A7G0j2xB2r4OuOJ0mmLAACfi8Bv7Y5UsotfvWrXzFz5kyuuOIKfvnLX7Jv3z4AHn/8cS677DKuv/56Zs+ezYc+9KE+SeOVV15h8eLFnHHGGZx88slcffXVdHZ29jn2zp07mTlzJps3b85se/LJJ7ngggsAuO+++zj77LM57bTTuOSSS3jjjTcyn33++ecD0N7ezpVXXsmpp57K3Llzuf7664nFYgX9mfTSX6UHSe8ChZF3N+MN9BYhWFg9cTxVZUQ3b2bkySc5Elv26Emn99Rw8qq08MqanZk+lIUW8HuZM2sMJ5v8q3WXL1/OZZddxtSpUznuuOP4+c9/zuWXXw7ASy+9xM0338wtt9zCE088wZIlSzj33HPx+/1cffXV3HzzzZx77rm0tLTwqU99iuXLl3PFFVdkjj1mzBjmzJnDM888w1VXXQXYCXHBggWsXr2aBx98kKeeeoqGhgbuvvtu7rzzTh5++OE+8d1zzz1UV1fz0ksvEQ6HufTSS3n++edZuHDh0f/ADkFHUIOkd4HCZDSKv6oKUiksy77/yQJHl9roO72nFXxq+HitubVoyQkgnkjxWnNr3vu/+eab7Nixg3POOQeAiy++mEcffZRUyo65vr6eRYsW4ff7WbhwIT09PezYsYPy8nJ+/vOfc+655xKJRGhtbWXUqFH93oO0YMECnn32WcAeDb3yyiucd955hEIhwuEwjz/+OOvXr+eLX/ziAckJoKqqitWrV/PCCy9gWRZPPvlkUZIT6AhqUHQ0r6N7927CGzaSSsTxlpUTGFFDvLMLsPB4PNQcd6xjFXzZ90DpCEoNJydNry/6COqk6fV57798+XI6OzuZN8/+5TWVStHe3s6KFSsAqKur23/s9MxMKpXC5/Pxhz/8gYceeohUKsWMGTPo6OjAsqwDPuOcc85hyZIlrFu3jldeeYU5c+YwevRoRo8ezT333MNDDz3Ej3/8Y+rq6rj66qv553/+5z7vv/rqq/F6vSxdupRrr72W2bNns2TJEqZMmXLYP5/DpWero9R77ckbCJCIRrESCZLhCL5gEF+wksqmJrw+H00XFOc3jv5ogYQark42DYc13VZMXV1dPPPMM9x///285z3vyWy/9957eeSRRzLXgPrz6quvctddd/E///M/mUTRO4WXq6qqivnz5/PrX/+aVatWceGFFwJ2141Ro0bx8MMPE41Gef7557nuuus466yz+ry/ubmZRYsW8eUvf5kdO3Zwxx13sGTJkn5HW4NNp/iOUG/Fnnz7vwi/u5nu3bvxlpXZN6pZFsmI3eoo2DTO8YUK9RqUUu7z9NNPM2bMGM444wzq6+szXxdddBEvv/wy0Wj0oO/t7OzE6/VSXl5OKpXiueee48UXXzxoF4fzzz+f5557jrVr1/KhD30IgPXr13PFFVfQ3NxMZWUldXV1lJWVHdCJ46c//Sm33XYb4XCYuro6KioqGDFixOD9IAagZ6vD0LvGU+f6jcTb2wjU1pLo7CTV02Pv4PPhDQTwlAWw4glIJh1fBwr6NorVBKWUOyxfvpzzzjvvgO3Tp0/n2GOP5Y477mDq1Kn9vvf9738///RP/8TChQvxer3MnDmTf/mXf0FE+t3/rLPO4oYbbmDu3LkEg0EAzjjjDK688kquvPJK9uzZw7hx4/j+979PbW1tn/def/313HTTTcybN49EIsGcOXO4+eabj/K7z4+nvznLocwYMwnY+Nvf/pbx48cPyjE7mtex7fEn6Fi9Bl8oSLK7m2RXGCuRAI8HUlnz2z4fHp8Pb8CPt6yMutNPY+rnPzcocRyp5/68ifVb9wLwkdOOYfrEkY7Go1ShrF27lpkzZzodhistWLCA6667jjPPPLOgn3Owv4OtW7fywQ9+EGCyiGzK51hF/XXaGHMicC9wArAB+IyIrDrS/QZb9iq43vIywEN3Wzvx9jZSPXFSiTjJ3W37E1NucgJIJrFSKSgvc7R7ebaYTvEpNWxt3ryZl19+mXA4zOmnn+50OIelaNegjDFlwFPAcqAWuA14wRhTcyT7DbbeYod4Zyf4/XSufZvOtW+T2LsHK54g0dmJlcxKRgONPH0+rJ64o93Ls+k1KKWGr+985zv84Ac/4NZbb8XrHVplB8U8W80FAiJyV/r5MmPMF4CLgPuPYL9B1XujrT8UIvLuWrzl5eDxEN+7F7w+8HqwkknAskdOlrV/9OTzQTJpP/Z68ZaX4cHjaPfybFrFp9TwtXTpUqdDOGLFPFvNAtbmbHsbOP4I9xtUsV0tlI227zlIRqOZBAVgJeLgD0BPD3i9QNboqTdZlZXhr6zAW16OFY/jr652vHqvV3nAR7Q7QWW5n4oyTVBKqaGhmGerKiCSsy0CBI9wv0FV0dhAvLMTfyiEr7LSrszzePDX1JDY1wFWCqusDC+Q6unBUxbAWxnE5/fRs3cfgeoqQlOn4isLkAxHXJOcAOa9dwKr1+9mxqRReL26Vo4qbZZl6ZpQDhnsortiJqgwkLvUZRDoOsL9BlX9/HlsecxuxFg+dizh5mbAQ9X0afS0tRHZvAVfWYCy2lo85WUk9u2jrK6OqimTCb5nKpF31hPb1UKgbhTjzl/gmuQE0FRfRVN9ldNhKFVwgUCAaDSaKaVWxRWNRjMdLwZDMRPUGuCanG0zgJ8d4X6Dqmb6NCYsXkTripUkd7VQPXMG4CHV3U319Gk0nvsPmSRU0dhw4P1NH/pgIcNTSuWhoaGBbdu20dTURGVlpY6kisSyLKLRKNu2baOxsXHQjlvMBLUS8BhjrgGWAhdil5E/cYT7Dbqa6dMGHvloElLK1Wpq7GLf7du3H7SrgiqMQCBAY2Nj5u9gMBQtQYlIjzHmo9j3Ny0BNgELRaTVGHMDcImIHDvQfsWKVSk1dNXU1AzqSVI5p6glXSKyGjirn+23A7cfaj+llFLDx9C6a0sppdSwoQlKKaWUK5XiXZs+gJ07dzodh1JKqbSsc7Iv3/eUYoIaC3DJJZc4HYdSSqkDjQXW57NjKSaoVcD7gR1A0uFYlFJK2XzYySnvlSlKbj0opZRSpUGLJJRSSrmSJiillFKupAlKKaWUK2mCUkop5UqaoJRSSrmSJiillFKupAlKKaWUK2mCUkop5Uql2EnCNYwxJ2Kva3UCsAH4jIjkfRd1qTDGfBj4FjANaAG+KyI/cTYqZxljaoE3gJtE5GGHwyk6Y8xY4MfAPCAG3CciX3c2KmcYY04HfggYoBX4log84GxU7qAjqAIxxpQBTwHLgVrgNuAFY8ywWknNGDMB+AVwK/bP4WLgDmPMOY4G5rx7gSang3DQU9jtyBqB04FPGWMWOxtS8RljvNg/ix+KyAjs/x9L07/cDns6giqcuUBARO5KP19mjPkCcBFwv2NRFd8k4DEReSL9fJUx5nfA+4BfOxWUk4wxnwJqgDedjsUJxpjTgCnA+0QkDmw0xswFoo4G5oyRQAPgMcZ4AAtIAD2ORuUSOoIqnFnA2pxtbwPHOxCLY0TkRRH5fO9zY8wo7Ga+rzoXlXOMMZOBbwCfcToWB83GTs7fNMZsM8asBz4mIjscjqvoRKQNWAr8FIhjN1K9QURyzx3DkiaowqkCIjnbIkDQgVhcwRgzAnga+Av2tMawYozxAf8NXCsiw3nBst5fUuLYI6kLgGuH8RRfDFgMVGLPvHzDGPMRJ+NyC53iK5ww9j+4bEGgy4FYHGeMmY6dlNYAl4hIyuGQnPB1QETkcacDcVg30CEi30w/f90Y8wB2onrMsaiccQH2VOdX089/b4x5ELgSeMG5sNxBR1CFswa7KifbjPT2YcUYczb2qOlJ4OMiEnM4JKcsAj5ujNlrjNmLPd37I2PMjxyOq9jeBoLpQqJew/WX5QlAec62BPboctgbrv8oimEl9oXPa7DnmC/ELjd/YsB3lRhjzFTgV8CNInK30/E4SURmZD83xrwG3DUMy8x/g11O/T1jzL9j/yJ3OXCVo1E54wXsqtbPYRdPnQJ8FrjC0ahcQkdQBSIiPcBHsRNTO3AjsFBEWh0NrPiuBqqx/xN2ZX192+nAlDPSI+gPYF9/2gE8D3xHRH7haGAOEJG3sKf5rgT2Yk9xXiciw+4abX90RV2llFKupCMopZRSrqQJSimllCtpglJKKeVKmqCUUkq5kiYopZRSrqQJSimllCvpjbpqyDHGhIDrsDvDTwB2Ao8Dt4jI3kH6jCrgIhF5cBCOFQA+KyKu6BhhjPk08F8iMroIn3UNUCsi3zjEftcCIRG5udAxqaFDR1BqSEmvp/Uy8CHgS9hd4z8PnIO93lZu/8Mj9e8MXmeDxcCSQTrWkJFeC+zLwHfz2P1u4NJ0z0alAE1Qauj5Fva/2/ki8ryIbBSRX2N37TgRuGyQPsczSMcZ7GMNJf8B/FxEDtkgWUS6sZecuKHgUakhQztJqCHDGFOOvWT8f4rIvf28fiZ2t/C29OJvXwS+AEzEblB6g4g8m973YezlTyqBjwNtwAMicmt6CuyhrENPxu4d911gIVCPPa34ExG5NX08H/bJ9bPYy0msSn92PXZfxl7zROR3OXF/E7tP4wbsdaK6gf8BrhGRZDrWKhH5eNZ7NmFP0y1Nv/94YB32qK8HuAX4O/ay6lOAP2B3kW/vneID7sQeKZLe7xsikkwf/7T067OBLdh94v5LRFLp938Ju/HxAuw2RbfmfE9V2G2MzhGRl9LbKoHbsVeNrQR+B/yriGxLv34c8Ddg/DBsCab6oSMoNZRMwV6JdlV/L4rIS+kF4MBOFjcDN2Gf/J8Ens5ZSvsKYCvwXuAB4BZjzGxgOfA94HVgLPYJ+k7gDOwEZYAfZu1P+nO+BHwFODn9nmeAP6e3taeP9dJBvrfzgNr0Z9yE3cNwYR4/k+z3V2E3G703He9S4F+Bj2Anmq9k7V+HPS06H7gUuxfcdQDGmAbs1Y6fx058X0rH8x9Z7z8Z6Ex/3s/6iecDQBL7++91L3A+8AnsZd4rsH/WAIjIauxfBD54GN+3KmGaoNRQMjL9576BdkqPnr4C3CYiy0SkOb320G/oe5JdLyJfF9ut2ElktohEsdftSojIzvSo4k/A5SLyiohsEJHvpveZlf68q4BbReRxEVmHfUL/BXaj3H2AlT7WwZbyjmKPJkRE7sdOjrMPsm9/YsBXRGQ99vUcH7A0vaLxn4BngWOz9k8CnxCR10XkGeBW7GRGOvZVInKLiKwTkeeA64Frcz7zZhF5R0Q29xPPe4G1ImJBZrHKxdijwv9Nrxj7eeCPOcturEm/Vymt4lNDyu70nyMH3AsagNH0/e0d4I/AP2c9X5fzeicQOMgx/xs4zxjzSWA6cBL2iMWX/qx6skZ2ItJJevrMmNxlwfr1bk7y6hgglv5sEpFE+nHvSs4bsl6PpePstU1EtmQ9XwWMM8bUYieyucaY7GtHXqDSGFPX+xmHWBW4kf1/X2D/zPz0/RltJD1qy9KG/fenlCYoNaS8g30CO5V+pvmMMXdiX/f4yUHe76HvrEF/o5mDFTT8H+zKwZ+lv67CHuVkH+doLugOFEt/x839v9vfAncDrVqczHne+3OJp4/9C+Br/byvd/QaHeDYvZ+d/bPM92fkw06mSukUnxo60svEPwJ80RhTkf2aMWYS9nWUbhHpALZjX8/JdiZ2sUQ+MidSY0w18EngUhG5QUSWY59wRwAeEdmHXbxxStZ7Ko0xu9LFBkdbidT7Wb3HruLoRxnjjTHZI9HTsUdxYWAtMCM9ffeOiLwDzMS+NjZQ0su2E3tU2WsDdlLM/hkdY4xpN8Y0Ze03Ov1epXQEpYacW4B/BFYaY76BPU13IvAd4A3gvvR+3wKWGGO2YleGXYRdFDA3z8/pAsYYY6ZgFzyEgQuMMRuBcdgVfR72L9f9feDr6eq6ZuwFKvcBrwLjgSpjzCxgwxEseb8KuNwYsxA7edzMgSOgw+UDHjPGXIddpXgD+0u87wG+ZIz5YfrxBOxR6VPpKr58jv834DpjjE9EkiLSaYx5ALjTGNOJPf33PeCN3iq+tBOwizuU0hGUGlpEpB04C/gr9klzDfaJ7kngH7JO/kuxk9Z3gDexq8fOE5EX8/yon2Nfk1qDfb1pMfYU3xrsKb4XsAsPegsZvgs8iF0N+Cp2xd4/pq8r/Ta97VXs5Hq4HgEeTn/ui8Br2EUbR2MtdhL5A3ZS/z7pqVER2YqdzN+LPY35CHa13Vf6PVL/VmIn8OyCh39Pb38Cu5qxE/sXBwCMMcdiF5X85ki+IVV69D4opVRBGGPuAZIi8qU8978DaBCRywsbmRoqdASllCqU7wAXpkvMB5S+iXdx+j1KAZqglFIFIiLvYk8dfjWP3b8A/FREpLBRqaFEp/iUUkq5ko6glFJKuZImKKWUUq6kCUoppZQraYJSSinlSpqglFJKudL/B70jyVo3IrHqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotSweepFrameDiff(frame)\n",
    "plot(frac_infected_series, label='Analysis')\n",
    "\n",
    "decorate(xlabel='Contact number (c)',\n",
    "         ylabel='Fraction infected')\n",
    "\n",
    "savefig('figs/chap06-fig04.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose you run a survey at the end of the semester and find that 26% of students had the Freshman Plague at some point.\n",
    "\n",
    "What is your best estimate of `c`?\n",
    "\n",
    "Hint: if you print `frac_infected_series`, you can read off the answer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000050    0.000100\n",
      "1.005083    0.010098\n",
      "1.010185    0.020096\n",
      "1.015356    0.030094\n",
      "1.020598    0.040092\n",
      "1.025914    0.050090\n",
      "1.031305    0.060088\n",
      "1.036772    0.070086\n",
      "1.042317    0.080084\n",
      "1.047943    0.090082\n",
      "1.053651    0.100080\n",
      "1.059444    0.110078\n",
      "1.065323    0.120076\n",
      "1.071291    0.130074\n",
      "1.077350    0.140072\n",
      "1.083503    0.150070\n",
      "1.089751    0.160068\n",
      "1.096099    0.170066\n",
      "1.102547    0.180064\n",
      "1.109099    0.190062\n",
      "1.115758    0.200060\n",
      "1.122527    0.210058\n",
      "1.129409    0.220056\n",
      "1.136407    0.230054\n",
      "1.143524    0.240052\n",
      "1.150765    0.250050\n",
      "1.158132    0.260048\n",
      "1.165630    0.270046\n",
      "1.173263    0.280044\n",
      "1.181034    0.290042\n",
      "              ...   \n",
      "1.743384    0.709958\n",
      "1.767898    0.719956\n",
      "1.793487    0.729954\n",
      "1.820238    0.739952\n",
      "1.848249    0.749950\n",
      "1.877628    0.759948\n",
      "1.908499    0.769946\n",
      "1.941002    0.779944\n",
      "1.975299    0.789942\n",
      "2.011573    0.799940\n",
      "2.050040    0.809938\n",
      "2.090947    0.819936\n",
      "2.134590    0.829934\n",
      "2.181315    0.839932\n",
      "2.231541    0.849930\n",
      "2.285771    0.859928\n",
      "2.344627    0.869926\n",
      "2.408879    0.879924\n",
      "2.479505    0.889922\n",
      "2.557767    0.899920\n",
      "2.645332    0.909918\n",
      "2.744467    0.919916\n",
      "2.858363    0.929914\n",
      "2.991711    0.939912\n",
      "3.151808    0.949910\n",
      "3.350924    0.959908\n",
      "3.612133    0.969906\n",
      "3.987365    0.979904\n",
      "4.642296    0.989902\n",
      "9.211261    0.999900\n",
      "Length: 101, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(frac_infected_series)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.158096819542062"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Alternative solution\n",
    "\n",
    "\"\"\"We can use `np.interp` to look up `s_inf` and\n",
    "estimate the corresponding value of `c`, but it only\n",
    "works if the index of the series is sorted in ascending\n",
    "order.  So we have to use `sort_index` first.\n",
    "\"\"\"\n",
    "\n",
    "frac_infected_series.sort_index(inplace=True)\n",
    "np.interp(0.26, frac_infected_series, frac_infected_series.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}