de_sim/examples/jupyter_examples_for_talk/3. An epidemic model using DE-Sim.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Requirement already satisfied: de_sim in /usr/local/lib/python3.7/site-packages (0.1.9)\n",
"Requirement already satisfied: matplotlib in /usr/local/lib/python3.7/site-packages (from de_sim) (3.1.2)\n",
"Requirement already satisfied: progressbar2>=3.39 in /usr/local/lib/python3.7/site-packages (from de_sim) (3.47.0)\n",
"Requirement already satisfied: wc-utils[git]>=0.0.16 in /root/.wc/wc_sandbox/packages/wc_utils (from de_sim) (0.0.24)\n",
"Requirement already satisfied: configobj in /usr/local/lib/python3.7/site-packages (from de_sim) (5.0.6)\n",
"Requirement already satisfied: pympler in /usr/local/lib/python3.7/site-packages (from de_sim) (0.9)\n",
"Requirement already satisfied: logging2 in /usr/local/lib/python3.7/site-packages (from de_sim) (0.1.2)\n",
"Requirement already satisfied: numpy in /usr/local/lib/python3.7/site-packages (from de_sim) (1.18.1)\n",
"Requirement already satisfied: setuptools in /usr/local/lib/python3.7/site-packages (from de_sim) (44.0.0)\n",
"Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/site-packages (from matplotlib->de_sim) (0.10.0)\n",
"Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/site-packages (from matplotlib->de_sim) (1.1.0)\n",
"Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.7/site-packages (from matplotlib->de_sim) (2.4.6)\n",
"Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/site-packages (from matplotlib->de_sim) (2.8.1)\n",
"Requirement already satisfied: six in /usr/local/lib/python3.7/site-packages/six-1.12.0-py3.7.egg (from progressbar2>=3.39->de_sim) (1.12.0)\n",
"Requirement already satisfied: python-utils>=2.3.0 in /usr/local/lib/python3.7/site-packages (from progressbar2>=3.39->de_sim) (2.3.0)\n",
"Requirement already satisfied: abduct in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (2.0.1)\n",
"Requirement already satisfied: attrdict in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (2.0.1)\n",
"Requirement already satisfied: dataclasses in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.6)\n",
"Requirement already satisfied: diskcache in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (4.1.0)\n",
"Requirement already satisfied: humanfriendly in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (4.18)\n",
"Requirement already satisfied: mendeleev in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.5.1)\n",
"Requirement already satisfied: natsort in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (6.2.0)\n",
"Requirement already satisfied: objsize in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.3.2)\n",
"Requirement already satisfied: openpyxl<=3.0.1,>=2.6.1 in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (3.0.1)\n",
"Requirement already satisfied: pint>=0.10 in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.10)\n",
"Requirement already satisfied: pronto>=1 in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (1.1.3)\n",
"Requirement already satisfied: pyexcel>=0.5.9.1 in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.5.15)\n",
"Requirement already satisfied: pyexcel_io>=0.5.9.1 in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.5.20)\n",
"Requirement already satisfied: pyyaml>=5.1 in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (5.3)\n",
"Requirement already satisfied: qualname in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (0.1.0)\n",
"Requirement already satisfied: requests in /usr/local/lib/python3.7/site-packages/requests-2.22.0-py3.7.egg (from wc-utils[git]>=0.0.16->de_sim) (2.22.0)\n",
"Requirement already satisfied: xlsxwriter in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (1.2.7)\n",
"Requirement already satisfied: gitpython in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (3.0.5)\n",
"Requirement already satisfied: pygithub in /usr/local/lib/python3.7/site-packages (from wc-utils[git]>=0.0.16->de_sim) (1.45)\n",
"Requirement already satisfied: contextlib2 in /usr/local/lib/python3.7/site-packages (from abduct->wc-utils[git]>=0.0.16->de_sim) (0.6.0.post1)\n",
"Requirement already satisfied: pandas in /usr/local/lib/python3.7/site-packages (from mendeleev->wc-utils[git]>=0.0.16->de_sim) (0.25.3)\n",
"Requirement already satisfied: pyfiglet in /usr/local/lib/python3.7/site-packages (from mendeleev->wc-utils[git]>=0.0.16->de_sim) (0.8.post1)\n",
"Requirement already satisfied: colorama in /usr/local/lib/python3.7/site-packages (from mendeleev->wc-utils[git]>=0.0.16->de_sim) (0.4.3)\n",
"Requirement already satisfied: sqlalchemy in /usr/local/lib/python3.7/site-packages (from mendeleev->wc-utils[git]>=0.0.16->de_sim) (1.3.12)\n",
"Requirement already satisfied: et-xmlfile in /usr/local/lib/python3.7/site-packages (from openpyxl<=3.0.1,>=2.6.1->wc-utils[git]>=0.0.16->de_sim) (1.0.1)\n",
"Requirement already satisfied: jdcal in /usr/local/lib/python3.7/site-packages (from openpyxl<=3.0.1,>=2.6.1->wc-utils[git]>=0.0.16->de_sim) (1.4.1)\n",
"Requirement already satisfied: fastobo~=0.6.0 in /usr/local/lib/python3.7/site-packages (from pronto>=1->wc-utils[git]>=0.0.16->de_sim) (0.6.1)\n",
"Requirement already satisfied: networkx~=2.3 in /usr/local/lib/python3.7/site-packages (from pronto>=1->wc-utils[git]>=0.0.16->de_sim) (2.4)\n",
"Requirement already satisfied: frozendict~=1.2 in /usr/local/lib/python3.7/site-packages (from pronto>=1->wc-utils[git]>=0.0.16->de_sim) (1.2)\n",
"Requirement already satisfied: contexter~=0.1.4 in /usr/local/lib/python3.7/site-packages (from pronto>=1->wc-utils[git]>=0.0.16->de_sim) (0.1.4)\n",
"Requirement already satisfied: nanoset~=0.1.2; platform_python_implementation == \"CPython\" in /usr/local/lib/python3.7/site-packages (from pronto>=1->wc-utils[git]>=0.0.16->de_sim) (0.1.3)\n",
"Requirement already satisfied: chardet~=3.0 in /usr/local/lib/python3.7/site-packages/chardet-3.0.4-py3.7.egg (from pronto>=1->wc-utils[git]>=0.0.16->de_sim) (3.0.4)\n",
"Requirement already satisfied: texttable>=0.8.2 in /usr/local/lib/python3.7/site-packages (from pyexcel>=0.5.9.1->wc-utils[git]>=0.0.16->de_sim) (1.6.2)\n",
"Requirement already satisfied: lml>=0.0.4 in /usr/local/lib/python3.7/site-packages (from pyexcel>=0.5.9.1->wc-utils[git]>=0.0.16->de_sim) (0.0.9)\n",
"Requirement already satisfied: idna<2.9,>=2.5 in /usr/local/lib/python3.7/site-packages/idna-2.8-py3.7.egg (from requests->wc-utils[git]>=0.0.16->de_sim) (2.8)\n",
"Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /usr/local/lib/python3.7/site-packages (from requests->wc-utils[git]>=0.0.16->de_sim) (1.24.2)\n",
"Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.7/site-packages/certifi-2019.9.11-py3.7.egg (from requests->wc-utils[git]>=0.0.16->de_sim) (2019.9.11)\n",
"Requirement already satisfied: gitdb2>=2.0.0 in /usr/local/lib/python3.7/site-packages (from gitpython->wc-utils[git]>=0.0.16->de_sim) (2.0.6)\n",
"Requirement already satisfied: pyjwt in /usr/local/lib/python3.7/site-packages (from pygithub->wc-utils[git]>=0.0.16->de_sim) (1.7.1)\n",
"Requirement already satisfied: deprecated in /usr/local/lib/python3.7/site-packages (from pygithub->wc-utils[git]>=0.0.16->de_sim) (1.2.7)\n",
"Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.7/site-packages (from pandas->mendeleev->wc-utils[git]>=0.0.16->de_sim) (2019.3)\n",
"Requirement already satisfied: decorator>=4.3.0 in /usr/local/lib/python3.7/site-packages (from networkx~=2.3->pronto>=1->wc-utils[git]>=0.0.16->de_sim) (4.4.1)\n",
"Requirement already satisfied: smmap2>=2.0.0 in /usr/local/lib/python3.7/site-packages (from gitdb2>=2.0.0->gitpython->wc-utils[git]>=0.0.16->de_sim) (2.0.5)\n",
"Requirement already satisfied: wrapt<2,>=1.10 in /usr/local/lib/python3.7/site-packages (from deprecated->pygithub->wc-utils[git]>=0.0.16->de_sim) (1.11.2)\n",
"\u001b[31mERROR: Error while checking for conflicts. Please file an issue on pip's issue tracker: https://github.com/pypa/pip/issues/new\n",
"Traceback (most recent call last):\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 3021, in _dep_map\n",
" return self.__dep_map\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 2815, in __getattr__\n",
" raise AttributeError(attr)\n",
"AttributeError: _DistInfoDistribution__dep_map\n",
"\n",
"During handling of the above exception, another exception occurred:\n",
"\n",
"Traceback (most recent call last):\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 3012, in _parsed_pkg_info\n",
" return self._pkg_info\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 2815, in __getattr__\n",
" raise AttributeError(attr)\n",
"AttributeError: _pkg_info\n",
"\n",
"During handling of the above exception, another exception occurred:\n",
"\n",
"Traceback (most recent call last):\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_internal/commands/install.py\", line 520, in _determine_conflicts\n",
" return check_install_conflicts(to_install)\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_internal/operations/check.py\", line 108, in check_install_conflicts\n",
" package_set, _ = create_package_set_from_installed()\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_internal/operations/check.py\", line 50, in create_package_set_from_installed\n",
" package_set[name] = PackageDetails(dist.version, dist.requires())\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 2736, in requires\n",
" dm = self._dep_map\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 3023, in _dep_map\n",
" self.__dep_map = self._compute_dependencies()\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 3032, in _compute_dependencies\n",
" for req in self._parsed_pkg_info.get_all('Requires-Dist') or []:\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 3014, in _parsed_pkg_info\n",
" metadata = self.get_metadata(self.PKG_INFO)\n",
" File \"/usr/local/lib/python3.7/site-packages/pip/_vendor/pkg_resources/__init__.py\", line 1895, in get_metadata\n",
" raise KeyError(\"No metadata except PKG-INFO is available\")\n",
"KeyError: 'No metadata except PKG-INFO is available'\u001b[0m\n"
]
}
],
"source": [
"!pip install de_sim"
]
},
{
"cell_type": "markdown",
"metadata": {
"codehighlighter": [
[
0,
0
],
[
0,
0
]
]
},
"source": [
"<!-- :Author: Arthur Goldberg <Arthur.Goldberg@mssm.edu> -->\n",
"<!-- :Date: 2020-07-13 -->\n",
"<!-- :Copyright: 2020, Karr Lab -->\n",
"<!-- :License: MIT -->\n",
"\n",
"# A stochastic epidemic model\n",
"\n",
"<font size=\"4\">\n",
"Epidemics occur when infectious diseases spread through a susceptible population.\n",
"Models that classify individuals by their infectious state are used to study the dynamics of epidemics.\n",
"The simplest approach considers three infectious states:\n",
"<br>\n",
"\n",
"* *Susceptible*: a person who can become infected if exposed\n",
"* *Infectious*: a person who is infected, and can transmit the infection to a susceptible person\n",
"* *Recovered*: a person who has recovered from an infection, and cannot be reinfected\n",
"\n",
"Dynamic analyses of epidemics are called Susceptible, Infectious, or Recovered (SIR) models.\n",
"SIR models are described by the initial population of people in each state and the rates at which they transition between states.\n",
"\n",
"\n",
"*SIR model states and transitions*\n",
"\n",
"S and I represent the number of individuals in states Susceptible and Infectious, respectively. β and γ are model parameters.\n",
"\n",
"We present a stochastic SIR model that demonstrates the core features of DE-Sim.\n",
"The SIR model uses DE-Sim to implement a continuous-time Markov chain model, as described in section 3 of Allen, 2017. Infectious Disease Modelling.\n",
"\n",
"Let's implement and use the SIR model.\n",
"\n",
"\n",
"\n",
"First, define the event messages\n",
"</font>"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"\"DE-Sim implementation of an SIR epidemic model\"\n",
"\n",
"import enum\n",
"import numpy\n",
"\n",
"import de_sim\n",
"\n",
"\n",
"class StateTransitionType(enum.Enum):\n",
" \" State transition types \"\n",
" s_to_i = 'Transition from Susceptible to Infectious'\n",
" i_to_r = 'Transition from Infectious to Recovered'\n",
"\n",
"\n",
"class TransitionMessage(de_sim.EventMessage):\n",
" \"Message for all model transitions\"\n",
" transition_type: StateTransitionType\n",
"\n",
"\n",
"MESSAGE_TYPES = [TransitionMessage]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"<font size=\"4\">\n",
"\n",
"Next, define a simulation object. It has these attributes:\n",
"<br>\n",
" \n",
"* s (int): number of susceptible subjects\n",
"* i (int): number of infectious subjects\n",
"* N (int): total number of susceptible subjects, a constant\n",
"* beta (float): SIR beta parameter\n",
"* gamma (float): SIR gamma parameter\n",
"* random_state (numpy.random.RandomState): a random state\n",
"\n",
"</font>"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"codehighlighter": [
[
14,
15
],
[
34,
35
],
[
49,
51
],
[
52,
53
],
[
44,
45
]
]
},
"outputs": [],
"source": [
"class SIR(de_sim.SimulationObject):\n",
" \"\"\"Implement a Susceptible, Infectious, or Recovered (SIR)\n",
" epidemic model\"\"\"\n",
" def __init__(self, name, s, i, N, beta, gamma):\n",
" \" Initialize an SIR instance \"\n",
" self.s = s\n",
" self.i = i\n",
" self.N = N\n",
" self.beta = beta\n",
" self.gamma = gamma\n",
" self.random_state = numpy.random.RandomState()\n",
" super().__init__(name)\n",
"\n",
" def init_before_run(self):\n",
" \" Send the initial events \"\n",
" self.schedule_next_event()\n",
"\n",
" def schedule_next_event(self):\n",
" \" Schedule the next SIR event \"\n",
" rates = {'s_to_i': self.beta * self.s * self.i / self.N,\n",
" 'i_to_r': self.gamma * self.i}\n",
" lambda_val = rates['s_to_i'] + rates['i_to_r']\n",
" if lambda_val == 0:\n",
" # no transitions remain\n",
" return\n",
" tau = self.random_state.exponential(1.0/lambda_val)\n",
" prob_s_to_i = rates['s_to_i'] / lambda_val\n",
" if self.random_state.random_sample() < prob_s_to_i:\n",
" self.send_event(tau, self, TransitionMessage(StateTransitionType.s_to_i))\n",
" else:\n",
" self.send_event(tau, self, TransitionMessage(StateTransitionType.i_to_r))\n",
"\n",
" def handle_state_transition(self, event):\n",
" \" Handle an infectious state transition event \"\n",
" transition_type = event.message.transition_type\n",
" if transition_type is StateTransitionType.s_to_i:\n",
" self.s -= 1\n",
" self.i += 1\n",
" elif transition_type is StateTransitionType.i_to_r:\n",
" self.i -= 1\n",
" self.schedule_next_event()\n",
"\n",
" event_handlers = [(TransitionMessage, 'handle_state_transition')]\n",
"\n",
" # register the message types sent\n",
" messages_sent = MESSAGE_TYPES"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from de_sim.checkpoint import AccessCheckpoints, Checkpoint\n",
"from de_sim.simulation_checkpoint_object import (AccessStateObjectInterface,\n",
" CheckpointSimulationObject)\n",
"\n",
"class AccessSIRObjectState(AccessStateObjectInterface):\n",
" \"\"\" Get the state of the simulation\n",
"\n",
" Attributes:\n",
" sir (`obj`): an SIR object\n",
" random_state (`numpy.random.RandomState`): a random state\n",
" \"\"\"\n",
"\n",
" def __init__(self, sir):\n",
" self.sir = sir\n",
" self.random_state = sir.random_state\n",
"\n",
" def get_checkpoint_state(self, time):\n",
" \"\"\" Get the simulation's state\n",
"\n",
" Args:\n",
" time (`float`): current time; ignored\n",
" \"\"\"\n",
" return dict(s=self.sir.s,\n",
" i=self.sir.i)\n",
"\n",
" def get_random_state(self):\n",
" \" Get the simulation's random state \"\n",
" return self.random_state.get_state()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"<font size=\"4\">\n",
"The next cell defines code to run the SIR model and visualize its predictions.\n",
"</font>"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"code_folding": [
22
]
},
"outputs": [],
"source": [
"import pandas\n",
"\n",
"class RunSIR(object):\n",
"\n",
" def __init__(self, checkpoint_dir):\n",
" self.checkpoint_dir = checkpoint_dir\n",
"\n",
" def simulate(self, recording_period, max_time, **sir_args):\n",
" \"\"\" Create and run an SIR simulation\n",
"\n",
" Args:\n",
" recording_period (`float`): interval between state checkpoints\n",
" max_time (`float`): simulation end time\n",
" sir_args (`dict`): arguments for an SIR object\n",
" \"\"\"\n",
" # create a simulator\n",
" simulator = de_sim.Simulator()\n",
"\n",
" # create an SIR instance\n",
" self.sir = sir = SIR(**sir_args)\n",
" simulator.add_object(sir)\n",
"\n",
" # create a checkpoint simulation object\n",
" access_state_object = AccessSIRObjectState(sir)\n",
" checkpointing_obj = CheckpointSimulationObject('checkpointing_obj',\n",
" recording_period,\n",
" self.checkpoint_dir,\n",
" access_state_object)\n",
" simulator.add_object(checkpointing_obj)\n",
"\n",
" # initialize simulation, which sends the SIR instance an initial event message\n",
" simulator.initialize()\n",
"\n",
" # run the simulation\n",
" event_num = simulator.simulate(max_time).num_events\n",
"\n",
" def last_checkpoint(self):\n",
" \"\"\" Get the last checkpoint of the last simulation run\n",
"\n",
" Returns:\n",
" `Checkpoint`: the last checkpoint of the last simulation run\n",
" \"\"\"\n",
" access_checkpoints = AccessCheckpoints(self.checkpoint_dir)\n",
" last_checkpoint_time = access_checkpoints.list_checkpoints()[-1]\n",
" return access_checkpoints.get_checkpoint(time=last_checkpoint_time)\n",
"\n",
" def history_to_dataframe(self):\n",
" fields = ('s', 'i', 'r')\n",
" hist = []\n",
" index = []\n",
" access_checkpoints = AccessCheckpoints(self.checkpoint_dir)\n",
" for checkpoint_time in access_checkpoints.list_checkpoints():\n",
" state = access_checkpoints.get_checkpoint(time=checkpoint_time).state\n",
" state_as_list = [state['s'], state['i'], self.sir.N - state['s'] - state['i']]\n",
" hist.append(dict(zip(fields, state_as_list)))\n",
" index.append(checkpoint_time)\n",
" return pandas.DataFrame(hist)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"<font size=\"4\">\n",
"Use the model to view an epidemic's predictions.\n",
"We use parameters from Allen (2017), and print and plot the trajectory of a single simulation.\n",
"Since the model is stochastic, each run produces a different trajectory.\n",
" </font>"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUZdrH8e+T3isJ6QUICb2FUEURpdiwUiwosrC76toV7GXftXdUlI6KYlnXithFBASSgNTQCYSEAOkJCWnP+8cZQlAgCWTmzCT357rmmpmTmTk3AfLLU4/SWiOEEEKcjpPZBQghhLB/EhZCCCEaJGEhhBCiQRIWQgghGiRhIYQQokEuZhdwNtq0aaPj4uLMLkMIIRxKWlraYa11SFPe49BhERcXR2pqqtllCCGEQ1FKZTb1PdINJYQQokFWCwul1Fyl1EGl1MZ6x4KUUt8rpbZb7gMtx5VS6jWl1A6l1HqlVG9r1SWEEKLprNmymA+M/NOxacCPWusE4EfLc4BRQILlNgWYYcW6hBBCNJHVwkJr/SuQ/6fDo4EFlscLgMvrHX9HG34HApRS4daqTQghRNPYesyirdY6x/L4ANDW8jgS2FfvdVmWY3+hlJqilEpVSqUeOnTIepUKIYSoY9oAtzZ2MGzyLoZa65la62StdXJISJNmfgkhhDhDtg6L3GPdS5b7g5bj+4Hoeq+LshwTQghhB2wdFl8AN1oe3wh8Xu/4BMusqP5AUb3uqlPKL6u0TpVCCCFOYM2psx8AK4FEpVSWUmoS8AxwoVJqO3CB5TnAYmAXsAOYBdzSmHPsLyznu00Hmr12IYQQJ1KOfPGjgJgkHXbjyyyaMoCe0QFmlyOEEA5BKZWmtU5uynscegV3XBtvQnzd+duCNezNO2J2OUII0WI5dFi4OCnmT0yhulZz0/zVFMgYhhBCWIVDhwVA+xAfZk1IJqugnCnvplJRVWN2SUII0eI4fFgA9I0L4sVrerBmTwH3fvwHtbWOOw4jhBD2yKG3KK/v0h4RZBeW8/Q3GUQGevLAqE5mlySEEC1GiwkLgClD2pFVUM7bS3cRFejFDf1jzS5JCCFahBYVFkopHru0M9mF5Tz2+UYi/D0Y1qltw28UQghxWi1izKI+F2cnpl/biy4R/tz2/lrWZxWaXZIQQji8FhcWAF5uLsy5KZkgbzdunp/KvnxZgyGEEGejRYYFQKivBwtu7ktldQ03zVtN0ZEqs0sSQgiH1WLDAqBDqC+zJiSzL7+cye+mcrRa1mAIIcSZaNFhAdCvXTDPX9Od1bvzue/j9bIGQwghzkCLmg11KqN7RrK/sJznlmwlKtCT+0cmmV2SEEI4lFYRFgD/PLc9WQXlvPnLTiIDPbmun6zBEEKIxmo1YaGU4snLupBTWM4jn20kwt+ToUmhZpclhBAOocWPWdTn4uzE69f2pnOEH7e+n87G/UVmlySEEA6hVYUFgLe7C3Nv7EuglxsT568hq0DWYAghRENaXVgAhPp5MG9iXyqqapg4bw1F5bIGQwghTqdVhgVAx7a+vH1DH/bklfF3WYMhhBCn1WrDAmBg+zY8f3UPft+Vz7T/bsCRr0cuhBDW1GpmQ53K5b0iySo4wgvfbSMq0JN7hieaXZIQQtidVh8WALcO7UBWQTnTf9pBZIAn41JizC5JCCHsioQFxhqMf1/eleyiCh76bCNh/h6clyhrMIQQ4phWPWZRn6uzE29e15vEtr7cujCdTdmyBkMIIY6RsKjHx92FeRP74u/pysR5a9hfWG52SUIIYRckLP6krZ8H8yamUF5Zw83z1lBcIWswhBBCwuIkEsN8eeuGPuw8VMo/30ujsrrW7JKEEMJUEhanMKhDG569qjvLd+Qx7dP1sgZDCNGqyWyo07iqTxRZBeW8/MM2ogK9uPvCjmaXJIQQppCwaMDtwzqQVXCE137cTlSgJ2OSo80uSQghbE7CogFKKZ66shsHiit48NMNhPl5MKRjiNllCSGETZkyZqGUuksptUkptVEp9YFSykMpFa+UWqWU2qGU+lAp5WZGbSdzbA1Gh1AfblmYzubsYrNLEkIIm7J5WCilIoHbgWStdVfAGRgHPAu8rLXuABQAk2xd2+n4ergyb2JffNxduHn+GnKKZA2GEKL1MKsbygXwVEpVAV5ADnA+cK3l6wuAx4EZplR3CuH+nsyb2Jdr3lrJxHlr+OgfA/DzcDW7LCHESdTqWsqqyiitLKW4sphqXW12SQ7N5mGhtd6vlHoB2AuUA98BaUCh1nV/m1lApK1ra4xO4X68eV1vbp6/hlsXpjP3pr64OssMZCGak9aaozVHKa0yftCXVpYaP/SrjMcllSV1t9Kq0rpAKK0yvlZaWUppVSkamfLeXGweFkqpQGA0EA8UAh8DI5vw/inAFICYGHN2hx3SMYSnruzG/Z+s58FPN/Dc1d1RSplSixD2SmtN4dFCco/kUlJZUvdDv6SyhJKqkrrH9QPh2POSyhKqak+/e4KTcsLH1QdfN1983XzxcfUh0ifyhOf1H7s5280wqOmGMrTJ7zGjG+oCYLfW+hCAUupTYBAQoJRysbQuooD9J3uz1nomMBMgOTnZtF8bxiRHk1VQzms/bmdA+2Cu7B1lVilCmKKmtoZD5YfILs0muyybnNKcuvucMuNWXn3qsT1PF098XS0/zN18CPQIJMY3Bh+34z/kfV19T3hePwC8XLzklzQbMiMs9gL9lVJeGN1Qw4BU4GfgamARcCPwuQm1NcmdwxL4OeMgL363jYu7h+Pu4mx2SUI0m4rqCg6UHfhLEBy7zz2SS40+8XLEge6BhPuE086/HQMjBhLhE0GYdxj+bv7GD31LOHi7eePqJON9jsSMMYtVSqlPgHSgGliL0VL4GliklPo/y7E5tq6tqZycFFNHJnH9nFW89/teJg2ON7skIRpFa01xZbERBidpGWSXZZNfkX/Ce5yUE6FeoUR4R9CrbS8ivCMI9wkn3DucCG8jFLxcvUz6EwlrM2U2lNb6MeCxPx3eBaSYUM5ZGZzQhsEd2vDGzzsYkxyFr8yOEnai6GgRe4r31P3wzy7NJqcsp+6+rKrshNe7O7sbP/h9IkgMSqx7fOw+xCtEWgOtmKzgbgZTRyZx6eu/MevXXdwt1/AWJjl45CDpuemk5qaSlpvGjsIdJ3zdz82PCJ8Ion2j6Rfej3Dv8BMCIcgjSMYAxClJWDSDblH+XNw9nNm/7eaGAXGE+LqbXZJo4bTWZJVkkZqbSvrBdNJy09hXsg8ALxcveob2ZFT8KJKCkuq6i7xdvU2uWjgyCYtmcu/wRL7deIDpP23nydFdzS5HtDC1upadhTtJy00jPdcIh4PlBwEIcA+gd2hvxiaOJbltMolBibg4yX9t0bzkX1QziW/jzbiUaN5ftZebB8UT10Z+ixNnrrq2moz8DNJy00jNTWXtwbUUHTWuCx/qFUqfsD70Ce1Dn7Z9aBfQDiclC0OFdUlYNKPbhyXw37T9vPj9NqaP72V2OcKBHK05yoZDG0jLTSMtN411h9bVrVGI8Y3h/Ojz6dO2D73b9ibKJ0rGFoTNSVg0o1BfDyYNjuf1n3fw9yHt6Brpb3ZJwk6VVpay7tC6ui6lDYc3UFVbhUKREJjA6Paj61oPIV6yJb4wn4RFM5tybjsWrsrk2SUZvDupn9nlCDtRUFFgBMNBo+WQkZ9Bra7FWTnTJbgL13e6nt5te9MrtBf+7vJLhrA/EhbNzM/DlVuHduD/vt7Cih2HGdihjdklCRMcKDtwwmD0zqKdgLGWoXtId6Z0n0Lv0N70COkhC9mEQ5CwsILr+8cy97fdPLskg89uHST9y63I+kPrmbV+Fr9k/QKAj6sPPUN7ckn7S0hum0zn4M6yoZ1wSBIWVuDh6sxdF3bkvk/W883GA1zULdzskoQVaa1JzU1l5vqZ/J7zO/7u/vyzxz8ZGj2UjoEdcXaSPcOE45OwsJIre0cxa9kuXvh2K8M7t8VFrnnR4mit+W3/b8xcP5N1h9YR7BHMPX3uYUziGOlaEi2OhIWVODsp7huRxOR3UvkoNYtr+5lz7Q3R/Gp1LT/t/YmZ62eyJX8L4d7hPNjvQa7ocAUeLh5mlyeEVUhYWNEFnUJJjg3klR+2cUWvSDzdpDvCkVXXVrNkzxJmr5/NzqKdxPrF8uTAJ7mk3SW4OssGe6Jlk7CwIqUUU0clcc1bK5m3Yje3nNfB7JLEGaisqeSLnV8wZ8Mcskqz6BDQgeeGPMfw2OEyHiFaDQkLK+sbF8SwpFBm/LKTa1NiCPCSmTCOory6nE+3f8q8jfPIPZJL1+Cu3Nf3Ps6LPk+21xCtjoSFDdw/MomRr/7KjF928sBFncwuRzSgtLKUD7d+yDub3yG/Ip8+bfvw5MAnGRAxQKZBi1ZLwsIGEsN8ubJXFPNW7OHGgXFEBHiaXZI4iaKjRby35T0WbllISWUJgyIGMbn7ZPq07WN2aUKYTsLCRu66MIEv/8jmlR+28dzVPcwuR9RzuPww72x6hw+3fsiR6iOcH30+U7pPoUubLmaXJoTdkLCwkahAL24YEMu85buZfE47Etr6ml1Sq5dTmsO8TfP4dPunVNVWMSJuBJO7TSYhMMHs0oSwOxIWNnTr0A58tGYfz3+7lZkTks0up9XaW7yXORvn8MXOL0DDZR0u4+auNxPrF2t2aULYLQkLGwrydmPKkHa8+P020jIL6BMbaHZJrcr2gu3M3jCbJXuW4KJcuKbjNUzsMpFwH9mORYiGSFjY2KRz4lmw0tjC/MMp/WV2jQ1sOryJmetn8tO+n/B08eTGzjcyocsE2njKjsBCNJaEhY15ublwx7AOPPL5Jn7ZeoihSaFml9RipeWmMWv9LJZnL8fXzZd/9PgH1yVdR4BHgNmlCeFwJCxMMC4lhtmWLczP7RiCk5O0LpqL1pqV2SuZuWEmablpBHkEcWfvOxmbOBYfNx+zyxPCYUlYmMDV2Yl7hidy+wdr+fyP/VzRK8rsklqEzOJMHl/xOKm5qYR6hTItZRpXJlyJp4usaxHibElYmOSSbuG8vXQnL363jYu6hePuInsMnama2hoWblnI9LXTcXV25eF+D3NFwhVykSEhmpFscGMSJyfF1JFJZBWU8/6qvWaX47B2F+3mpiU38Xzq8/QL78dnoz9jbNJYCQohmpm0LEx0TkIbBrYP5vWfdnBNcjQ+7vLX0Vg1tTW8u/ldXl/3Ou7O7jw1+CkuaXeJzC4TwkqkZWEipYzWRV5ZJbN+3WV2OQ5jV+EuJiyZwItpLzIoYhCfjf6MS9tfKkEhhBXJr7Im6xEdwMXdwpm9bBfX948lxNfd7JLsVnVtNQs2LeDNdW/i5erFc0OeY2TcSAkJIWxAWhZ24J7hHamoruX1n7abXYrd2lGwgxsW38Ar6a8wJGoI/xv9P0bFj5KgEMJGpGVhB9qF+DC2bzTvr97LpMHtiAn2Mrsku1FdW828jfOY8ccMfFx9eP7c5xkRO0JCQggbM6VloZQKUEp9opTKUEptUUoNUEoFKaW+V0ptt9y3qo2T7hiWgLOT4sXvt5pdit3YVrCN6xZfx2trX+P8mPP57PLPpNtJCJOY1bJ4FViitb5aKeUGeAEPAj9qrZ9RSk0DpgFTTarP5tr6eXDzoHje/GUnU4a0o0uEv9klmaaqtoo5G+bw9vq38XPz46XzXuLC2AvNLkuIv6iqqiIrK4uKigqzSzkpDw8PoqKicHV1PevPUlrrhl+k1CDgcSAWI2AUoLXW7Zp8QqX8gXVAO13v5EqprcB5WuscpVQ48IvWOvF0n5WcnKxTU1ObWoLdKiqvYshzP9MzOoAFN6eYXY4pMvIzeGT5I2TkZzAqfhQPpDxAoEeramQKB7J79258fX0JDg62uxav1pq8vDxKSkqIj48/4WtKqTStdZOuk9DYlsUc4C4gDahpyglOIh44BMxTSvWwfOYdQFutdY7lNQeAtid7s1JqCjAFICYm5ixLsS/+nq7cOrQ9Ty3OYOXOPAa0Dza7JJupqqli1oZZzFo/C393f14Z+grDYoaZXZYQp1VRUUFcXJzdBQUYU/ODg4M5dOhQs3xeY8csirTW32itD2qt847dzvCcLkBvYIbWuhdQhtHlVMfS4jhpk0drPVNrnay1Tg4JCTnDEuzXhAFxhPt78MySDBrT6msJtuRtYdzX45jxxwxGxo/k88s/l6AQDsMeg+KY5qytsWHxs1LqectAdO9jtzM8ZxaQpbVeZXn+CUZ45Fq6n7DcHzzDz3doHq7O3HVBR/7YV8i3mw6YXY5VVdZUMn3tdMZ/PZ6CigKmnz+dp895Gn/31jteI4S9amw3VD/Lff0+Lg2c39QTaq0PKKX2KaUStdZbgWHAZsvtRuAZy/3nTf3sluLK3pHMXLaL577dygWd2uLi3PKWw2w6vImHlz/MjsIdjG4/mvv63ichIYQda1RYaK2HNvN5/wUstMyE2gVMxGjlfKSUmgRkAmOa+ZwOw8XZiftGJPL3d9P4JC2LcSktZ2ymsqaSGX/MYN7GeQR7BvPGsDcYEjXE7LKEEA1oVFhYZjA9Bhz7X70UeFJrXXQmJ9Var+PEVsox0lFtMbxzW3rHBPDKD9u5vFckHq6Ov4X5hkMbeGT5I+ws2skVHa7g3r734ufmZ3ZZQjissrIyxowZQ1ZWFjU1NTzyyCOMHTvWKudqbDfUXGAjx3/bvwGYB1xpjaKEMTA1bVQnxry9kvkr9vCPc9ubXdIZO1pzlDfWvcGCTQsI8QzhrQveYlDkILPLEqJZPfHlJjZnFzfrZ3aO8OOxS7uc8utLliwhIiKCr7/+GoCiojP6/b1RGtsZ3l5r/ZjWepfl9gTQ5DUWomlS4oM4PymUN3/eQdGRKrPLOSPrDq7jmi+vYd7GeVyZcCWfjf5MgkKIZtKtWze+//57pk6dyrJly/D3t964X2NbFuVKqcFa69+gbpFeudWqEnXuH5nIqFeX8ebSHTwwqpPZ5TRaRXUFr699nXc2v0OYdxhvX/g2AyMGml2WEFZzuhaAtXTs2JH09HQWL17Mww8/zLBhw3j00Uetcq7GhsU/gQWWsQsF5AM3WaUicYKkMD+u6BnJ/OV7uGlgHOH+9n896fTcdB5d8SiZxZmMTRzLXX3uwtvV2+yyhGhxsrOzCQoK4vrrrycgIIDZs2db7VyNnQ21DuihlPKzPG/ejjlxWndd2JGv1ufw6g/beeaq7maXc0rl1eW8lv4aC7csJMIngtnDZ9MvvF/DbxRCnJENGzZw33334eTkhKurKzNmzLDauU4bFkqp67XW7yml7v7TcQC01i9ZrTJRJzrIi+v6x7BgxR7+dk47OoT6mF3SX2SXZjP5u8nsLdnL+KTx3Nn7TrxcZat1IaxpxIgRjBgxwibnamiA+1jfge9Jbvb3E6sFu21oB7zcXHjhW/vbwry6tpqpv04lryKPuSPm8mC/ByUohGhhTtuy0Fq/bXn4g9Z6ef2vWQa5hY0E+7gz+Zx2vPzDNtbuLaBXjP3sxDpz/UzWHVrHM+c8Q9+wvmaXI4SwgsZOnZ3eyGPCiv52TjxtfNx41o42GUzLTePt9W9zWfvLuLjdxWaXI4SwkobGLAYAA4GQP41b+AGOv6TYwXi7u/Cv8xN47ItNLN12iPMSQ02tp+hoEdOWTSPSJ5IH+z1oai1CCOtqqGXhhjE24cKJ4xXFwNXWLU2czPiUGGKCvHh2yVZqa81rXWiteWLlExw+cpjnhjwnU2OFaOEaGrNYCixVSs3XWmfaqCZxGm4uTtwzvCN3LFrHl+uzGd0z0pQ6/rfjf3yf+T139r6Trm26mlKDEMJ2GjtmccRyPYvFSqmfjt2sWpk4pUu7R9A53I8Xv9tGZXWtzc+/q2gXz6x+hn5h/ZjYdaLNzy+EOG7gQNvsjNDYsFgIZGBcEvUJYA+wxko1iQY4OSmmjkpib/4RPli916bnrqypZOqvU3F3duepc57CSbW8a20I4UhWrFhhk/M09n96sNZ6DlCltV6qtb6ZM7jwkWg+QxLaMKBdMNN/2k7Z0WqbnfeV9FfIyM/g34P+TaiXuQPsQgjw8bHNkrfG7g11bMvTHKXUxUA2EGSdkkRjKGW0Li5/Yzmzl+3mjgsSrH7O3/b/xrub32Vc4jjOiz7P6ucTwqF8Mw0ObGjezwzrBqOead7PPEONbVn8n2UTwXuAe4HZwF1Wq0o0Ss/oAEZ1DWPmrzs5XHrUquc6XH6Yh357iA4BHbgn+R6rnksIYX8au5HgV5aHRUBzX2JVnIV7RyTy3eZcXvp+G09d0c0q56jVtTy8/GHKqsqYPXw2Hi4eVjmPEA7NTloA1tLQorzpwCkn82utb2/2ikSTtA/x4aaBccz5bTeRAZ7cOrRDs5/jvc3vsXz/ch7q9xAJgdbv7hJC2J+GWhapNqlCnJUHL+pEflklz3+7FVdnxZQhzXcJ1i15W3g5/WWGRg9lbKJ1ru0rhLB/DS3KW2CrQsSZc3ZSPH91d6pqanlqcQauzk5MHBR/1p97pOoI9/96P0HuQTwx8Im6remFEPajtLTUJudp1JiFUupnTtIdpbWW6bN2wsXZiZfH9qS6RvPEl5txcXbihv6xZ/WZz615jsziTGYNn0Wgh/3sciuEsL3GTp29t95jD+AqwHaT+0WjuDo78dr4XtyyMI1HPtuIm7NibN+YM/qs7/Z8x3+3/5dJXSfJ1e6EEI2eDZX2p0PLlVKrrVCPOEtuLk68cV1vpryTxrRPN+Di5MRVfaKa9Bk5pTk8vvJxurXpxq29brVSpUIIR9KodRZKqaB6tzZKqRGAv5VrE2fI3cWZt2/ow6D2bbjvkz/4fN3+Rr+3praGacumUVNbw7PnPIurk6sVKxVCOIrGdkOlYYxZKIzup93AJGsVJc6eh6szsyYkM3H+au7+6A9cnZ24qFt4g++buWEm6QfTeWrwU0T7RdugUiGEI2hsN9TZT60RNufp5sycG/ty49zV3P7BWlycFMO7hJ3y9WsPruWtP97i4nYXc2n7S21YqRDC3jW2G8pDKXW3UupTpdR/lVJ3KqVkGa8D8HZ3Yd7EvnSN9OfW99P5OePgSV9XXFnMtF+nEeEdwcP9HrZxlUIIe9fYvaHeAbpgXHf7dcvjd61VlGhevh6uLLg5haQwP/7+Xhq/bjt0wte11vx75b/JPZLLs0OexcfNNrtYCiGaj9aa2lrrXd+msWHRVWs9SWv9s+U2GSMwhIPw93Tl3UkptA/xYfI7qazYcbjua5/t+Iwle5Zwa89b6R7S3cQqhRBNsWfPHhITE5kwYQJdu3Zl3759VjtXYwe405VS/bXWvwMopfohW4E4nAAvN96blML4Wb8zaUEqC25OITSomKdXP03fsL7c3PVms0sUwmE9u/pZMvIzmvUzk4KSmJoy9bSv2b59OwsWLKB///7Neu4/a2zLog+wQim1Rym1B1gJ9FVKbVBKrT+TEyulnJVSa5VSX1mexyulVimldiilPlRKuZ3J54rTC/ZxZ+Hf+hMR4MHEeSu47Yd7cHN24+nBT+Ps5Gx2eUKIJoqNjbV6UEDjWxYjrXDuO4AtgJ/l+bPAy1rrRUqptzCm5s6wwnlbvRBfd96f3J9L3ptKZuk27uz2H9p6tzW7LCEcWkMtAGvx9va2yXka1bLQWmcCAcCllluA1jrz2K2pJ1VKRQEXY1xECWXsUHc+8InlJQuAy5v6uaLxdpakU+H9E25HBvHqF+5s3F9kdklCCDvW2KmzdwALgVDL7T2l1L/O4ryvAPcDx4bug4FCrfWx/aaygMhT1DJFKZWqlEo9dOjQyV4iGpBXnseDvz1Ie//2LLrqKXw9XLl+ziq25BSbXZoQwk41dsxiEtBPa/2o1vpRoD8w+UxOqJS6BDh4kv2mGkVrPVNrnay1Tg4JCTmTj2jVtNY8svwRSipLeO7c50gIDeL9yf3wcHHm+tmr2J5bYnaJQohGiouLY+PGjTY5V2PDQgE19Z7XWI6diUHAZZaB8kUY3U+vAgFKqWNjKFFA4zc0Eo32fsb7LNu/jHuS76FjYEcAYoO9eX9yP5ycFONnrWLnIdvsjy+EcByNDYt5wCql1ONKqceB34E5Z3JCrfUDWusorXUcMA74SWt9HfAzcLXlZTcCn5/J54tT25q/lRdTX+TcqHMZnzT+hK+1C/Hhg8n90Fpz7azf2XO4zKQqhRD2qLED3C8BE4F8y22i1vqVZq5lKnC3UmoHxhjGGYWROLny6nLu+/U+AtwDeHLQkye96l2HUF8WTu5HZXUt1876nX35R0yoVAjHovVfrgtnN5qzttOGhWVPqDuVUq8DfYE3tdavaa3XNsfJtda/aK0vsTzepbVO0Vp30Fpfo7U+2hznEIbn1zzPnqI9/GfwfwjyCDrl65LC/Hjvb/0oq6xh/KzfyS4st2GVQjgWDw8P8vLy7DIwtNbk5eXh4dE82/g1tM5iAVAFLANGAZ2AO5vlzMJmfsj8gY+3fczErhMZEDGgwdd3ifDn3UkpXDdrFeNn/c6HUwYQ5i/7RgrxZ1FRUWRlZWGvMzM9PDyIimraxc9ORZ0uEZVSG7TW3SyPXYDVWuvezXLmZpCcnKxTU2XXkdM5UHaAq764imjfaN4d9S6uzo2/mFH63gJumL2Ktv4eLJrSn1BfCQwhWgKlVJrWOrkp72lozKLq2IN6ayCEg6ipreGBZQ9QVVvFs0OebVJQAPSOCWT+zSnkFFZw3axV5JVKz6AQrVVDYdFDKVVsuZUA3Y89VkrJCi47N2fjHFJzU3mo30PE+sWe0Wf0jQti7k192VdwhOtmr6KgrLKZqxRCOILThoXW2llr7We5+WqtXeo99jvde4W51h1cx5vr3mRU/Cgua3/ZWX3WgPbBzJqQzK7DZVw/ZxVFR6oafpMQokVp7DoL4UBKKkuYtmwaYd5hPNL/kZNOk22qcxJCePv6PmzLLWHC3FUUV0hgCNGaSFi0MFpr/v37vzlQdoBnznkGXzffZvvsoUmhvHFtbzZlFzNx3hpKj8owlhCthYRFC/Plri/5Zvc3/LPHP+kZ2rPZP394lzCmj+/Fun2F3Dx/DUcqJTCEaA0kLFqQvcrDC6UAABkJSURBVMV7+c/v/yG5bTJ/6/Y3q51nVLdwXhrTg9Q9+fxtQSoVVTUNv0kI4dAkLFqIqpoq7v/1flycXHj6HOtf9W50z0iev7oHK3flMfkdCQwhWjoJixZi+rrpbMrbxBMDnyDMO8wm57yqTxTPXNmNZdsPc8vCdCqraxt+kxDCIUlYtAArs1cyb+M8ru54NRfEXmDTc4/tG8P/Xd6VnzIOcuv76ew8VGqX++QIIc5OY6/BLezU4fLDPLDsAdr5t+P+vvebUsP1/WOpqqnliS838/3mXNr4uJESH0RKXBB944NICvPD2ensp+8KIcwjYeHAamprmPrrVMqqypg1fBaeLp6m1TJxUDznJYby+648Vu/OZ/XufBZvOACAr4cLfeOC6BsXREp8EN0i/XFzkUatEI5EwsKBvbX+LVYfWM2TA58kITDB7HKIb+NNfBtvxqfEAJBVcIQ1e4zgWLU7n58yDgLg4epE75hA+sYF0S8+iF4xgXi6WXdAXghxdiQsHNSK7BW8/cfbXNb+Mq5IuMLsck4qKtCLqEAvruhlbJF8qOQoqXuM4FizJ5/XftqO1uDipOgW5U9KvBEefWKD8Pds2qaHQgjrOu0W5fautW5RfvDIQa758hoC3QN5/+L38XL1MrukM1JUXkV6ZkFdeKzPKqSqRqOUcRGmfvFGt1XfuCBCfN3NLleIFuNMtiiXsHAw1bXV/O27v7E5bzOLLl5Eu4B2ZpfUbMora1i3r9AY89iTR3pmIeWW9Rvt2ngbg+aWW1SgYwakEPbgTMJCuqEczBvr3iAtN42nBj/VooICwNPNmQHtgxnQPhhIoKqmlo37i+oNmOewaM0+ACL8PSzBEUxKfCDtQ3yaZcNEIcTJSVg4kGVZy5i9YTZXJVzFpe0vNbscq3N1dqJXTCC9YgL5+7ntqa3VbM0tqQuP5Tvz+GxdNgDB3m7GjKv4IHrFBJAU5ouXm/zzFqK5SDeUgzhQdoBrvryGUK9QFl60EA8XucSp1po9eUdYvTuP1bsLWL0nj3355QAoBbFBXnQK9yMpzI+kcF86hfkRFeiJk6z5EK2cdEO1UFW1Vdy39D4qayp58dwXJSgslFJ103XH9jWm6+YUlbMhq4iMAyVsySkm40AJSzYd4NjvRD7uLiSG+ZIU5ktSuB+dw33p2NYXXw+ZfSXE6UhYOIDp6dNZd2gdzw15jjj/OLPLsWvh/p6E+3syvMvx/bGOVFaz9UAJGQdKyMgpZsuBEr74I5uFq/bWvSY6yJOkMD86hfvRyRIksUFe0goRwkLCws4t3beUeZvmMabjGEbFjzK7HIfk5eZSN/ZxjNaa7KIKMiytj805xWTkFPPjllxqLa0QT1dnEsN86RTuWxckiWG+sgZEtEoyZmHHskuzuebLa4j0ieTdi97F3VnWGlhbRVUN23JLyMgpYcuB4rr7wnrXHY8M8CQpzNcYD7EESXwbb9n/SjgMGbNoQapqqrh36b3U6lpePPdFCQob8XB1pntUAN2jAuqOaa3JLT56PDxyisk4UMwv2w5RY2mGuLs4HR8LsbRCekT7y4ws0WLIv2Q79VLaS2w4vIGXznuJaL9os8tp1ZRShPl7EObvwdDE0LrjR6tr2HGwlC05JXXdWT9uOchHqVnAiduYpMQFkRwbhL+XdGEJxyRhYYd+zPyR97a8x7VJ13Jh7IVmlyNOwd3FmS4R/nSJ8K87prXmUOlRNu0vrttEce5vu3l76a66bUxS4gJJiQ+mb3wgob4ys004BhmzsDP7SvYx9suxxPrFsmDUAtyc3cwuSZyliqoa1u4trAuPtMyCum1M4tt4kxJXfxsTT1mJLqxOxiwcXGVNJfcuvRcUPH/u8xIULYSHa/1tTKCqppZN2cWWxYT5LNl0gA9TjW1Mwi3bmBzbvr1DqGxjIuyDhIUdeSH1BTbnbeaVoa8Q5RtldjnCSlydnegZHUDP6ACmDDG2Mdl28Pg2Jit35vG5ZRuTIG83kmMDLdu3B9Mp3BcXZ7lwlLA9m4eFUioaeAdoC2hgptb6VaVUEPAhEAfsAcZorQtsXZ9Zvt3zLR9kfMCEzhMYFjOsaW8+WgI1VeAVZJ3ihFU5OSljS5IwPyYMiENrTWbeEcvuu0aAfLc5FwBvN2f6WFodKfFBdI/yx91FLhwlrM/mYxZKqXAgXGudrpTyBdKAy4GbgHyt9TNKqWlAoNZ66uk+q6WMWWQWZzL2q7G0D2jP/JHzcXVq5IyZqnL4fQb89gpUlkL7odB9HCRdBG7e1i1a2NSBogpLcBhdV9tySwFwczFaKf0sXVe9YwPxcZcOA3F6Dnk9C6XU58Drltt5WuscS6D8orVOPN17W0JYHK05yvWLryenLIePL/mYcJ/wht9UWwPrFsLPT0NJNnQcCSFJsOETKM4CNx/odCl0HwvxQ8BJfvNsaQrKKusGzFfvyWdTdjE1tRpnJ0XXCL+6651HB3khQx7izzqF+ztWWCil4oBfga7AXq11gOW4AgqOPT+VlhAWT658ko+3fcwbw95gSNSQ079Ya9i2BH54HA5lQGQfuPBJiBtsfL22FvaugD8WwebP4Wgx+IZDt6uN4AjrZvU/jzBH6dFq0jML6sJj3b5CKqtrzS5L2KnMZy9xnLBQSvkAS4H/aK0/VUoV1g8HpVSB1jrwJO+bAkwBiImJ6ZOZmWmzmpvb4l2LmbpsKhO7TuTuPnef/sX71sD3jxphENQehj0KnUdzyl8bqypg2zew/iPY/h3UVkNoF+g+BrpdA/6Rzf8HEnajoqqG9VlF5JUeNbsUYYcu6h7hGGGhlHIFvgK+1Vq/ZDm2lVbUDbW7aDfjvhpHYlAic0bMOfU4xeEd8OMTsOUL8A6Bc6dCn5vAuQkrgcvyYNOnsP5DyFoDKKN7qvtYo7vKw685/khCCAfhEGMWli6mBRiD2XfWO/48kFdvgDtIa33/6T7LUcOivLqc6xZfx+Ejh/no0o8I8w7764tKcmHpM5C2AFw8YNDtMOA2cPc5u5Pn7TRaG+s/hILd4OJpDIh3Hwvtz29aCAkhHJKjhMVgYBmwATjWqfogsAr4CIgBMjGmzuaf7rMcNSweW/EY/9v+P9684E0GRw4+8YtHS2DFdFjxOtQcNVoR504Fn9CTftYZ09poZfyxyGh1lBeAVxvL+MYYiOh96i4uIYRDc4iwaE6OGBZf7PyCh357iMndJnN779uPf6G6EtIXwNJnoewQdL7cGJcIbm/9oqorYcf3Rmtj6xIjpIITjNZG9zEQGGv9GoQQNiNhYed2Fu5k/Nfj6RLchVnDZ+Hi5GL8hr/pf/Djk0a3UOxgY4ZTVB9ziiwvNGZSrf8QMpcbx2IGGMHR5XLw/MucAyGEg5GwsGNHqo5w7dfXUnC0gE8u/YQQrxDYvcyY4ZSdDqGd4YLHIWG4/XT/FO49Pr5xeBs4u0HHEUZwJAwHF7nGhhCOSDYStFNaa/6z6j/sKtrFzOEzCSk5CJ/+0+j68YuE0W9Cj3H2t3guIAaG3Avn3AM564zg2PAxbPkSPAKgyxVG3dH97CfghBBWIS0LG/jf9v/x6IpH+WfitdySsxf++MCYrjr4buj3d3D1NLvExquphl2/wPpFsOUrqC6HgFjL+MZYaNPB7AqFEA2Qbig7tK1gG9d+PZ6ezv68vWMjzgD9phhB4egb/x0tMQJj/SLYtRTQxvqNlCnQcRQ4S8NVCHskYWFnyo7kMe6z0ZRWFPBxVg5tuo2BoQ8a3TstTXG20WJKnQdF+8AvCpInQu8bwSfE7OqEEPVIWNiL2hr0H4uYtuZplrjBbOdo+l74XOvYm6mm2ti/avVM2L3UGBTvcqXR2jBrhpcQ4gQywG02rWHHD/D9Y3xcnsniNkH8K/YS+p73tNmV2Y6zC3S6xLgd2gqrZxktjvWLjIV+KVOMgXFXufa0EI5EWhbNZX8afP8Y7FnGljaxXO/nRN/wfrx5wQycVCu/sllFsTH9dvVMYwquV7DRPZV8MwREm12dEK2OdEOZIW8n/PRvY2GdVxtKB9/F2AOLqag5yseXfkyQh4MPYjcnrY2uqdWzYOti41jiRUZrI36ITL8VwkakG8qWjpbCr8/DyjeMzfeG3I8ecBuPrXqS/aXZzB0xV4Liz5SCducZt8K9kDrX2Cgx4ytokwgpk411G+6+5tYphPiLVt4/cga0hk2fwRspsPwVY++k29fC+Q+xKHMx32V+x796/YvebXubXal9C4gxVqzfvQUun2GsNVl8L7zYCRbfD4e3m12hEKIe6YZqisM74Jv7YOdP0LYbXPwixPQDYFPeJm5YfAP9w/vz+rDXZZyiqbQ2xn1Wz4SNn0JtFbQbalmzMcL+VrcL4cBkzMJaKstg2Yuw/DXjN+DzH4bkSXWLzoorixnz5RhqdA0fX/IxAR6nvRqsaEjpQWMH3jVzjWuM+8dA30nQe4LjL2QUwg5IWDQ3rY3+9CUPGAvNeow3doStd20JrTV3/XIXS/ctZf6o+fQI6WG9elqbmmrY+rUxIL5nmXERqK5XG2MbET3Nrk4IhyUD3M0pbyd8M9XY7C+0C0z8BmIH/uVl7215jx/3/si9yfdKUDQ3ZxfjOuOdR0PuZlgzy7hY07r3ICrF6KLqPBpc3MyuVIgWT1oWf1Z5BH572Ri8dnY3tudImfKXfY5KKkuYu3Eu8zfO55yoc3h16KsomfppfeWFxiK/1bMgfyd4hxpXE0yeCH4RZlcnhEOQbqizlbEYlkw1pnV2uwaG/x/4nnh97MqaShZlLGLmhpkUHS3ioviLeKj/Q/i5+TVfHaJhtbWw6ycjNLZ9C8oJOl1qBHvsQFmzIcRpSDfUmcrfDUumGXsahSTBjV9B/DknvKRW17J492JeX/s6+0v30z+8P3f1uYvOwZ1NKrqVc3KCDhcYt/zdkDoH0t+FzZ9BUDtjsV/HkRDT31gHI4Q4K627ZVFVDstfhWUvGT9QzpsG/f7xlx8uK7JX8EraK2zJ30JSUBJ39bmLgRF/Hb8QJqs8Ahv/awTG7l+hphI8/KHDhZA4CjoMk8vCCoF0QzXNtm/hm/uhYI+xK+qI//ylz3tL3hZeTnuZlTkrifSJ5LZet3FR/EWyhsIRHC2FXT/D1iWw/VsoOwTK2eii6jjSCI/g9mZXKYQpJCwaoyDTmAq79Wto0xEuet7YfqKerJIspq+dzuLdi/F392dKtymMSxqHm7PMunFItbXGgr9t38DWb+DgZuN4cAIkjjS6rKJS5GJNotWQsDidqgpYMR2WvWAMhp47FfrfcsK0y4KKAmaun8mHWz/ESTlxQ+cbmNh1ogxetzQFmcb41NZvYM9vxmpxz0BIGG60OjoMM7qvhGihJCxOZccPsPg+yN9lzMsf8RT4R9V9uby6nIVbFjJnwxyOVB/h8g6Xc0uPW2jr3daK1Qu7UFFsbN+ybYnRNVmeD04uEDvI6KrqOBKC4s2uUohmJWHxZ4X74NsHYMuXENTe6HLqMKzuy9W11Xy+43PeXPcmB8sPcl70edzZ+07aB0hfdqtUWwP7Vlu6q5bA4a3G8ZAkyzjHRRCVLPtUCYcnYXFMdSWsfN3YQlxrGHIvDPwXuLgDxhYdv+z7hVfTX2Vn0U66h3Tn7j5306etXPZT1JO/ywiNbd9A5gqorTYu3JQwwhjraH++bKcuHJKEBcDOn40up7ztkHQJjHza2A7bYt3Bdbyc9jLpB9OJ84vjjt53MCxmmKy+FqdXXgg7fzTGObZ/DxWFxvXF4wZDx1FGeNT7dyaEPWvdYVG0H7590JhjHxhvdDklXFj32t1Fu3kt/TV+2PsDwR7B3NLzFq5IuAJXJ1mwJZqophr2/W4Ex7YlkLfDOB7a5fjsqojexsJBIexQ6wyL31fAqhnwy7Oga+Cce2Dg7eDqAcChI4eY8ccMPt3+Ke7O7kzsOpEJnSfg5eplcvWixTi84/g4x96Vxr9DD39w8zG7MnD3M65z7h99/P7YY58wCbRWqvWFRbdEnfp3f2MgMvEio8spMA6Asqoy5m2cxzub36GqpooxiWOY0n0KwZ7B5hYtWrYj+cbsu8wVxpRcM2mM7rLCvcYW++UFJ37dyRX8Iy3hEWMJkqh6oRJVN84nWpbWFxYRzjr1/i4w6jmj+Q9U1VTx8baPeXv92+RX5DMybiS397qdaL9ok6sVwmRHS43QKNxn3Nd/XLgPSnIwEqYenzBLeESdGCrHAsVD1iA5otYXFomROnXjDnD1pFbX8t2e73ht7WvsK9lHSlgKd/W5i65tuppdphCOoboSivdbgiTLEiR764VLlrHfVn3u/n/t5qrf3eUTKjsA2yGH33VWKTUSeBVwBmZrrZ857Rt8w8HVk9U5q3kp7SU25W0iITCBN4e9yeDIwTLDSYimcHEzFiCeahFibS2UHTxFy2QvZC6Ho8UnvsfZ3dhzzTPQaIW4+xnjOcdudc/96j2v91jWtNgNu2lZKKWcgW3AhUAWsAYYr7XefKr3dO3ZVQ9+YTDL9y8n3Duc23rdxsXxF+Ms/8CEMEd54V9bJsXZUFFkuRUb90eLoepIw5/n5nvqIPlLyJwkgFy9pGVzEo7eskgBdmitdwEopRYBo4FThsWuol34H/Ln3uR7GZc0DndnGYwTwlSeAcYtrFvDr62pMsLj6EmC5ITHRccfl+bC4e3Hn9dWn/4cTi7HQ8bFo3n+jK2UPYVFJLCv3vMsoN+fX6SUmgJMAQiODWbxlcbOsEIIB+PsCt7Bxu1MaG20TiqK64VKsTEDrP7zY4//PN7Sqq1u8jvsKSwaRWs9E5gJxjoLCQohWimlwM3buBFudjWOZey7TX6LPa3I2Q/Un98aZTkmhBDCZPYUFmuABKVUvFLKDRgHfGFyTUIIIbCjbiitdbVS6jbgW4yps3O11ptMLksIIQR2FBYAWuvFwGKz6xBCCHEie+qGEkIIYackLIQQQjRIwkIIIUSDJCyEEEI0yG72hjoTSqkSYKvZddiJNsBhs4uwE/K9OE6+F8fJ9+K4RK11ky4gb1ezoc7A1qZuhtVSKaVS5XthkO/FcfK9OE6+F8cppVKb+h7phhJCCNEgCQshhBANcvSwmGl2AXZEvhfHyffiOPleHCffi+Oa/L1w6AFuIYQQtuHoLQshhBA2IGEhhBCiQQ4bFkqpkUqprUqpHUqpaWbXYxalVLRS6mel1Gal1Cal1B1m12QmpZSzUmqtUuors2sxm1IqQCn1iVIqQym1RSk1wOyazKCUusvyf2OjUuoDpVSrur6qUmquUuqgUmpjvWNBSqnvlVLbLfeBDX2OQ4aFUsoZeAMYBXQGxiulOptblWmqgXu01p2B/sCtrfh7AXAHsMXsIuzEq8ASrXUS0INW+H1RSkUCtwPJWuuuGJc/GGduVTY3Hxj5p2PTgB+11gnAj5bnp+WQYQGkADu01ru01pXAImC0yTWZQmudo7VOtzwuwfiBEGluVeZQSkUBFwOzza7FbEopf2AIMAdAa12ptS40tyrTuACeSikXwAvINrkem9Ja/wrk/+nwaGCB5fEC4PKGPsdRwyIS2FfveRat9AdkfUqpOKAXsMrcSkzzCnA/UGt2IXYgHjgEzLN0y81WSnmbXZStaa33Ay8Ae4EcoEhr/Z25VdmFtlrrHMvjA0Dbht7gqGEh/kQp5QP8F7hTa11sdj22ppS6BDiotU4zuxY74QL0BmZorXsBZTSiq6GlsfTFj8YIzwjAWyl1vblV2RdtrJ9ocA2Fo4bFfiC63vMoy7FWSSnlihEUC7XWn5pdj0kGAZcppfZgdEuer5R6z9ySTJUFZGmtj7UyP8EIj9bmAmC31vqQ1roK+BQYaHJN9iBXKRUOYLk/2NAbHDUs1gAJSql4pZQbxoDVFybXZAqllMLol96itX7J7HrMorV+QGsdpbWOw/j38JPWutX+Bqm1PgDsU0olWg4NAzabWJJZ9gL9lVJelv8rw2iFA/0n8QVwo+XxjcDnDb3BIXed1VpXK6VuA77FmN0wV2u9yeSyzDIIuAHYoJRaZzn2oOV65qJ1+xew0PIL1S5gosn12JzWepVS6hMgHWPm4Fpa2bYfSqkPgPOANkqpLOAx4BngI6XUJCATGNPg58h2H0IIIRriqN1QQgghbEjCQgghRIMkLIQQQjRIwkIIIUSDJCyEEEI0yCGnzgpha0qpYIwN1wDCgBqM7TQAjmitZaGXaNFk6qwQTaSUehwo1Vq/YHYtQtiKdEMJcZaUUqWW+/OUUkuVUp8rpXYppZ5RSl2nlFqtlNqglGpveV2IUuq/Sqk1ltsgc/8EQjRMwkKI5tUD+AfQCWNlfUetdQrGtun/srzmVeBlrXVf4CpkS3XhAGTMQojmtebY1s9KqZ3Ase2wNwBDLY8vADobWxUB4KeU8tFal9q0UiGaQMJCiOZ1tN7j2nrPazn+/80J6K+1rrBlYUKcDemGEsL2vuN4lxRKqZ4m1iJEo0hYCGF7twPJSqn1SqnNGGMcQtg1mTorhBCiQdKyEEII0SAJCyGEEA2SsBBCCNEgCQshhBANkrAQQgjRIAkLIYQQDZKwEEII0aD/B2JTnOufdJZBAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import tempfile\n",
"\n",
"sir_args = dict(name='sir',\n",
" s=98,\n",
" i=2,\n",
" N=100,\n",
" beta=0.3,\n",
" gamma=0.15)\n",
"with tempfile.TemporaryDirectory() as tmpdirname:\n",
" run_sir = RunSIR(tmpdirname)\n",
" run_sir.simulate(10, 100, **sir_args)\n",
"\n",
" # print and plot an epidemic's predicted trajectory\n",
" sir_data_frame = run_sir.history_to_dataframe()\n",
" axes = sir_data_frame.plot()\n",
" axes.set_xlabel(\"Time\")\n",
" rv = axes.set_ylabel(\"Population\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"<font size=\"4\">\n",
"An important prediction generated by the SIR model is the severity of the epidemic, which can be summarized by the fraction of people who became infected.\n",
"We run an ensemble of simulations and examine the predicted distribution of severity.\n",
" </font>"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAXTklEQVR4nO3de5gldX3n8fdHUBBBBKedICLjhaDERCQT48YbBnVBI+KuYXXFgFExUWNc0c1EiaKJEe9usiaKQkBEVsQbCkQBUaIr4oCD3LzioFyEAVQuojDwzR9VrU3blxqm65zprvfrec4zdTv1+57qns+p/lWd30lVIUkajruNuwBJ0mgZ/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvzZKkuVJzkpyY5J3LvC+L0qy50Luc1ORpJI8tJ1+X5K/W6D9PjDJTUk2a+e/mORFC7Hvdn+nJjlwofan8dh83AVo05NkLfCiqjq9w+YHA9cC966N+FBIkqOBy6vq0MllVfU7d3V/C6l98/lwVT2gj/1X1V90rGMt8/xcquqHwNYLUVeSw4CHVtUBU/a/z0LsW+PlGb821s7AxRsT+uOUxpL4f5DEEzl1U1U+fNzpAawFntxOHwR8GXgH8BPgB8A+7bqjgduAW4GbgCfTnEysAr4PXAecAGw/Zd+PA/4/8FPgR+3+D562n8/MUMcWwHuAK9vHe4At2nV7ApcDhwDXAFcBL5jj9X0ReDPwFeAW4KHAC4BLgBuBS4GXtNveq93mjra2m4D7z/c6Z2jzNW1dVwJ/DhTN2fTkcfyHdnoZ8Nn2+FwP/Efb1rFtDbe0NfxvYEW7nxcCPwTOmrJs8ymv9S3AOcANwKcn65w8bjP97IG925/HbW1750/Z34va6bsBhwKXtcf9Q8C27brJOg5sa7sWeN24f7d9NI8lcaaj3v0h8G2aUHobcGSSVNVBwHHA26pq62q6IP4K2A94Ik1A/gR4L0CSnYFTgX8GJoDdgTVVdcS0/TxjhhpeBzymfc4jgUfThM6k3wK2BXakCcL3Jtlujtf0fJo3nG34dXD9CXBvmjeBdyfZo6puBvYBrmxr27qqrpzrdU6XZG/g1cBTgF1ognU2h9C8iU0Ay4HXAlVVz6cJ0Ge0NbxtynOeCDwc+K+z7PPPaN5sdgDWA/80R/vQNPjvwD8CH23be+QMmx3UPp4EPJimi+n/TtvmccCuwF7A65M8fL621T+DX11cVlUfqKrbgWNoAmT5LNv+Bc2Z3eVV9UvgMODZbTfE/wROr6rjq+q2qrquqtZ0rOF5wJuq6pqqWge8kSa8J93Wrr+tqk6hOUvddY79HV1VF1XV+vY5J1fV96vxJeDzwOPneP5cr3O6/YF/q6oL2zeSw+bY7200x3fntq7/qKr5utEOq6qbq+qWWdYfO6XtvwP2n7z4u5GeB7yrqi6tqpuAvwWeM+0YvLGqbqmq84Hzad60NWYGv7r48eREVf28nZztAuLOwCeT/DTJT2m6T26neaPYiaZr5K64P82Z+aTL2mWTrquq9VPmfz5HjdB0M/1Kkn2SnJ3k+rbup9H8hTObuV7nTLVPbe+yGbaZ9Hbge8Dnk1yaZNUc2874WuZZfxlwd+Z+bV3N9DPZnDsfgx9PmZ7vZ6IRMfi10H5Ecw3gPlMeW1bVFe26h8zyvPnOaq+kCdtJD2yX3VW/ai/JFsDHaa5jLK+q+wCnAJmjtrle53RX0bzpTa195qKqbqyqQ6rqwcC+wKuS7DVHHXMtnzS97dto+txvBraaXNH+FTCxAfud6WeyHrh6nudpzAx+LbT3AW9u+/NJMpHkme2644AnJ9k/yeZJ7ptk93bd1TT9xLM5Hji03d8y4PXAhxeo5nvQXDxeB6xPsg/w1Cnrrwbum2TbKcvmep3TnQAclGS3JFsBb5itkCR/kuShSQL8jOaviDum1DHXMZrNAVPafhNwYttt9x1gyyRPT3J3mmsmW0x53tXAijnuejoe+F9JHpRka359TWD9LNtrE2Hwa6H9H+Akmq6KG4GzaS4OU8095k+juYB5PbCGX/f5Hgns1nadfGqG/f4DsBr4JnABcF67bKNV1Y3AK2gC+ic01yJOmrL+WzQhd2lb3/3nep0z7P9UmruQvkDTjfOFOcrZBTid5hrFV4F/qaoz23VvoXnz+2mSV2/ASzyW5s6hHwNbtq+VqvoZ8FLgg8AVNH8BXD7leR9r/70uyXkz7Peodt9n0dzt9Quai97axGX+60aSpKXEM35JGhiDX5IGxuCXpIEx+CVpYBbFoE7Lli2rFStWjLsMSVpUzj333GuramL68kUR/CtWrGD16tXjLkOSFpUkM35K3K4eSRoYg1+SBsbgl6SBMfglaWAMfkkaGINfkgbG4JekgTH4JWlgDH5JGphF8cndjbFi1clja3vt4U8fW9uSNBvP+CVpYAx+SRoYg1+SBsbgl6SBMfglaWAMfkkaGINfkgbG4JekgTH4JWlgDH5JGhiDX5IGxuCXpIEx+CVpYAx+SRqYJT8ss6SlY1zDrC+1IdY945ekgTH4JWlgDH5JGpjegj/JTknOTHJxkouS/HW7/LAkVyRZ0z6e1lcNkqTf1OfF3fXAIVV1XpJtgHOTnNaue3dVvaPHtiVJs+gt+KvqKuCqdvrGJJcAO/bVniSpm5H08SdZATwK+Fq76OVJvpnkqCTbzfKcg5OsTrJ63bp1oyhTkgah9+BPsjXwceCVVXUD8K/AQ4Ddaf4ieOdMz6uqI6pqZVWtnJiY6LtMSRqMXoM/yd1pQv+4qvoEQFVdXVW3V9UdwAeAR/dZgyTpzvq8qyfAkcAlVfWuKct3mLLZs4AL+6pBkvSb+ryr57HA84ELkqxpl70WeG6S3YEC1gIv6bEGSdI0fd7V82UgM6w6pa82JUnz85O7kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvyQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA2PwS9LAGPySNDAGvyQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sD0FvxJdkpyZpKLk1yU5K/b5dsnOS3Jd9t/t+urBknSb+rzjH89cEhV7QY8BnhZkt2AVcAZVbULcEY7L0kakd6Cv6quqqrz2ukbgUuAHYFnAse0mx0D7NdXDZKk3zSSPv4kK4BHAV8DllfVVe2qHwPLZ3nOwUlWJ1m9bt26UZQpSYPQe/An2Rr4OPDKqrph6rqqKqBmel5VHVFVK6tq5cTERN9lStJg9Br8Se5OE/rHVdUn2sVXJ9mhXb8DcE2fNUiS7qzPu3oCHAlcUlXvmrLqJODAdvpA4NN91SBJ+k2b97jvxwLPBy5IsqZd9lrgcOCEJC8ELgP277EGSdI0vQV/VX0ZyCyr9+qrXUnS3PzkriQNjMEvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA1Mp+BP8rt9FyJJGo2uZ/z/kuScJC9Nsm2vFUmSetUp+Kvq8cDzgJ2Ac5N8JMlTeq1MktSLzn38VfVd4FDgb4AnAv+U5FtJ/ltfxUmSFl7XPv7fS/Jumu/N/WPgGVX18Hb63T3WJ0laYF2HZf5n4IPAa6vqlsmFVXVlkkN7qUyS1Iuuwf904Jaquh0gyd2ALavq51V1bG/VSZIWXNc+/tOBe06Z36pdJklaZLoG/5ZVddPkTDu9VT8lSZL61DX4b06yx+RMkt8Hbplje0nSJqprH/8rgY8luZLme3R/C/gfvVUlSepNp+Cvqq8neRiwa7vo21V1W39lSZL60vWMH+APgBXtc/ZIQlV9qJeqJEm96RT8SY4FHgKsAW5vFxdg8EvSItP1jH8lsFtVVZ/FSNKmaMWqk8fW9trDn77g++x6V8+FNBd0JUmLXNcz/mXAxUnOAX45ubCq9u2lKklSb7oG/2F9FiFJGp2ut3N+KcnOwC5VdXqSrYDN+i1NktSHrsMyvxg4EXh/u2hH4FPzPOeoJNckuXDKssOSXJFkTft42l0tXJJ013S9uPsy4LHADfCrL2W53zzPORrYe4bl766q3dvHKV0LlSQtjK7B/8uqunVyJsnmNPfxz6qqzgKu34jaJEk96Br8X0ryWuCe7Xftfgz4zF1s8+VJvtl2BW0320ZJDk6yOsnqdevW3cWmJEnTdQ3+VcA64ALgJcApNN+/u6H+leYTwLsDVwHvnG3DqjqiqlZW1cqJiYm70JQkaSZd7+q5A/hA+7jLqurqyekkHwA+uzH7kyRtuK5j9fyAGfr0q+rBG9JYkh2q6qp29lk0nwiWJI3QhozVM2lL4E+B7ed6QpLjgT2BZUkuB94A7Jlkd5o3kbU03UaSpBHq2tVz3bRF70lyLvD6OZ7z3BkWH7kBtUmSetC1q2ePKbN3o/kLYEPG8pckbSK6hvfUu2/W03TT7L/g1UiSete1q+dJfRciSRqNrl09r5prfVW9a2HKkST1bUPu6vkD4KR2/hnAOcB3+yhKktSfrsH/AGCPqroRmlE2gZOr6oC+CpMk9aPrkA3LgVunzN/aLpMkLTJdz/g/BJyT5JPt/H7AMf2UJEnqU9e7et6c5FTg8e2iF1TVN/orS5LUlw35ENZWwA1V9W9JJpI8qKp+0FdhkjZNK1adPO4StJG6fvXiG4C/Af62XXR34MN9FSVJ6k/Xi7vPAvYFbgaoqiuBbfoqSpLUn67Bf2tVFe3QzEnu1V9JkqQ+dQ3+E5K8H7hPkhcDp7ORX8oiSRqPrnf1vKP9rt0bgF2B11fVab1WJknqxbzBn2Qz4PR2oDbDXpIWuXm7eqrqduCOJNuOoB5JUs+63sd/E3BBktNo7+wBqKpX9FKVJKk3XYP/E+1DkrTIzRn8SR5YVT+sKsflkaQlYr4+/k9NTiT5eM+1SJJGYL7gz5TpB/dZiCRpNOYL/pplWpK0SM13cfeRSW6gOfO/ZztNO19Vde9eq5MkLbg5g7+qNhtVIZKk0eg6Vo8kaYkw+CVpYAx+SRqY3oI/yVFJrkly4ZRl2yc5Lcl323+366t9SdLM+jzjPxrYe9qyVcAZVbULcEY7L0kaod6Cv6rOAq6ftviZwOTwD8cA+/XVviRpZqPu419eVVe10z8Gls+2YZKDk6xOsnrdunWjqU6SBmBsF3enfofvLOuPqKqVVbVyYmJihJVJ0tI26uC/OskOAO2/14y4fUkavFEH/0nAge30gcCnR9y+JA1en7dzHg98Fdg1yeVJXggcDjwlyXeBJ7fzkqQR6voNXBusqp47y6q9+mpTkjQ/P7krSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA9PbJ3cl9WvFqpPHXYIWKc/4JWlgDH5JGhiDX5IGxuCXpIEx+CVpYAx+SRoYg1+SBsbgl6SBMfglaWAMfkkaGINfkgbG4JekgTH4JWlgDH5JGhiDX5IGxuCXpIEx+CVpYAx+SRqYsXz1YpK1wI3A7cD6qlo5jjokaYjG+Z27T6qqa8fYviQNkl09kjQw4wr+Aj6f5NwkB8+0QZKDk6xOsnrdunUjLk+Slq5xBf/jqmoPYB/gZUmeMH2DqjqiqlZW1cqJiYnRVyhJS9RYgr+qrmj/vQb4JPDocdQhSUM08uBPcq8k20xOA08FLhx1HZI0VOO4q2c58Mkkk+1/pKr+fQx1SNIgjTz4q+pS4JGjbleS1PB2TkkaGINfkgbG4JekgTH4JWlgDH5JGhiDX5IGxuCXpIEx+CVpYAx+SRqYcX4Ry5K3YtXJY2l37eFPH0u7QzSun7G0MTzjl6SBMfglaWAMfkkaGINfkgbG4JekgTH4JWlgvJ1TS4K3VUrdecYvSQNj8EvSwBj8kjQwBr8kDYzBL0kDY/BL0sB4O6cWjLdUSouDZ/ySNDAGvyQNjMEvSQNj8EvSwIwl+JPsneTbSb6XZNU4apCkoRp58CfZDHgvsA+wG/DcJLuNug5JGqpxnPE/GvheVV1aVbcC/w945hjqkKRBGsd9/DsCP5oyfznwh9M3SnIwcHA7e1OSb8+z32XAtQtS4eK2LG/1OODvwySPQ2PRHoe8daOevvNMCzfZD3BV1RHAEV23T7K6qlb2WNKi4HFoeBwaHoeGx+HOxtHVcwWw05T5B7TLJEkjMI7g/zqwS5IHJbkH8BzgpDHUIUmDNPKunqpan+TlwOeAzYCjquqiBdh1526hJc7j0PA4NDwODY/DFKmqcdcgSRohP7krSQNj8EvSwCyq4J9vqIckWyT5aLv+a0lWjL7K/nU4Dq9KcnGSbyY5I8mM9/IuBV2H/0jy35NUkiV5S1+X45Bk//b34qIkHxl1jaPQ4f/GA5OcmeQb7f+Pp42jzrGrqkXxoLkQ/H3gwcA9gPOB3aZt81Lgfe30c4CPjrvuMR2HJwFbtdN/uRSPQ9dj0W63DXAWcDawctx1j+l3YhfgG8B27fz9xl33mI7DEcBfttO7AWvHXfc4HovpjL/LUA/PBI5pp08E9kqSEdY4CvMeh6o6s6p+3s6eTfNZiaWo6/Affw+8FfjFKIsboS7H4cXAe6vqJwBVdc2IaxyFLsehgHu309sCV46wvk3GYgr+mYZ62HG2bapqPfAz4L4jqW50uhyHqV4InNprReMz77FIsgewU1Ut5e+F7PI78dvAbyf5SpKzk+w9supGp8txOAw4IMnlwCnAX42mtE3LJjtkgzZekgOAlcATx13LOCS5G/Au4KAxl7Ip2Jymu2dPmr8Az0ryu1X107FWNXrPBY6uqncm+S/AsUkeUVV3jLuwUVpMZ/xdhnr41TZJNqf5U+66kVQ3Op2GvEjyZOB1wL5V9csR1TZq8x2LbYBHAF9MshZ4DHDSErzA2+V34nLgpKq6rap+AHyH5o1gKelyHF4InABQVV8FtqQZwG1QFlPwdxnq4STgwHb62cAXqr2Ks4TMexySPAp4P03oL8W+3ElzHouq+llVLauqFVW1guZ6x75VtXo85famy/+NT9Gc7ZNkGU3Xz6WjLHIEuhyHHwJ7ASR5OE3wrxtplZuARRP8bZ/95FAPlwAnVNVFSd6UZN92syOB+yb5HvAqYMl9u1fH4/B2YGvgY0nWJFmSYyF1PBZLXsfj8DnguiQXA2cCr6mqJfXXcMfjcAjw4iTnA8cDBy3Bk8N5OWSDJA3MojnjlyQtDINfkgbG4JekgTH4JWlgDH5JGhiDX4takps6bPP4dkTKNUnuuYH73y/JblPm39R+OG4kktwnyUtH1Z6GweDXEDwPeEtV7V5Vt2zgc/ejGcURgKp6fVWdvpDFtZ8yn819aEadlRaMwa8lIcmeSb6Y5MQk30pyXBovAvYH/j7Jce22r0ny9XY89jdO2ceftcvOT3Jskj8C9gXe3v618JAkRyd5drv9Xu247hckOSrJFu3ytUnemOS8dt3DZqj3oCQnJfkCcEaSrdvvTph8zuSokocDD2nbf/tc9UtdOUiblpJHAb9DM9TuV4DHVtUHkzwO+GxVnZjkqTRj1DwaCM3YPU+gGdPpUOCPquraJNtX1fXtp54/W1UnAkyO8p1kS+BoYK+q+k6SD9F898F72lqurao92m6aVwMvmqHePYDfa9vZHHhWVd3QDqlwdtv2KuARVbV72+6M9VfVWQt1ELX0ecavpeScqrq8HWlxDbBihm2e2j6+AZwHPIwmSP8Y+FhVXQtQVdfP09auwA+q6jvt/DHAE6as/0T777mz1AFw2pR2Avxjkm8Cp9MMJ7x8A+qXOvOMX0vJ1FFIb2fm3+/Q9Pe//04Lk4Uel32yltnqALh5yvTzgAng96vqtnY00S1neM6M9UsbwjN+Dc3ngD9PsjVAkh2T3A/4AvCnSe7bLt++3f5GmuGdp/s2sCLJQ9v55wNf2oi6tgWuaUP/ScDk9yRPb3+2+qXOPOPXoFTV59vheL/a9tffBBzQjuL4ZuBLSW6n6Uo5iObr+z6Q5BU0Q31P7ucXSV5AMwLq5jRDAr9vI0o7DvhMkguA1cC32nauS/OtWRcCp1bVa2aqH1jKw29rgTk6pyQNjF09kjQwBr8kDYzBL0kDY/BL0sAY/JI0MAa/JA2MwS9JA/OfAzz56mkSRr0AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import math\n",
"\n",
"num_sims = 100\n",
"infection_rates= []\n",
"\n",
"for _ in range(num_sims):\n",
" with tempfile.TemporaryDirectory() as tmpdirname:\n",
" run_sirs = RunSIR(tmpdirname)\n",
" run_sirs.simulate(recording_period=10, max_time=60, **sir_args)\n",
" # infection rate = infectious + recovered\n",
" # N = s + i + r => i + r = N - s\n",
" final_state = run_sirs.last_checkpoint().state\n",
" N = sir_args['N']\n",
" infection_rate = (N - final_state['s'])/N\n",
" infection_rates.append(infection_rate)\n",
"\n",
"import matplotlib.pyplot as plt\n",
"fig = plt.figure()\n",
"ax = fig.add_subplot(111)\n",
"rv = plt.hist(infection_rates)\n",
"ax.set_title('Infection rate distribution')\n",
"ax.set_xlabel('Infection rate')\n",
"rv = ax.set_ylabel('Frequency')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<font size=\"4\">\n",
"As predicted by Allen's (2017) analysis, for the parameters in `sir_args` the infection rate distribution is bimodal. Most of the epidemics infect a majority of the population and a small fraction of them (Allen predicts 25%) burn out and infect only a minority of the population.\n",
"\n",
"</font>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"<font size=\"4\">\n",
"The simple model above only touches the surface of epidemic modeling. Many extensions are possible:\n",
"\n",
"* A spatial model with multiple geographic areas: each area would be represented by an instance of SIR.\n",
"* An extension of the spatial model that also represents travel between geographic areas\n",
"* A model that represents individuals in more states, such as multiple infectious states which distinguish between asymptomatic and symptomatic individuals, with a lower transmission parameter β for symptomatic individuals who would likely isolate while recovering\n",
"* A model that can model both small and large populations: it would use the stochastic approach above to integrate small populations and ODEs to integrate large populations. Models and simulators that use multiple integration methods are called *multi-algorithmic*.\n",
"\n",
"We encourage you to experiment with different parameters for this model and build your own models!\n",
" </font>"
]
},
{
"cell_type": "markdown",
"metadata": {
"codehighlighter": [
[
2,
4
]
]
},
"source": [
"**References**\n",
"\n",
"Allen, L.J., 2017. A primer on stochastic epidemic models: Formulation, numerical\n",
"simulation, and analysis. Infectious Disease Modelling, 2(2), pp.128-142."
]
}
],
"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.7.6"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": false,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {
"height": "calc(100% - 180px)",
"left": "10px",
"top": "150px",
"width": "165px"
},
"toc_section_display": true,
"toc_window_display": false
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}