{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Modeling and Simulation in Python\n",
"\n",
"Chapter 13\n",
"\n",
"Copyright 2017 Allen Downey\n",
"\n",
"License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
]
},
{
"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": "markdown",
"metadata": {},
"source": [
"`make_system`, `plot_results`, and `calc_total_infected` are unchanged."
]
},
{
"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 plot_results(S, I, R):\n",
" \"\"\"Plot the results of a SIR model.\n",
" \n",
" S: TimeSeries\n",
" I: TimeSeries\n",
" R: TimeSeries\n",
" \"\"\"\n",
" plot(S, '--', label='Susceptible')\n",
" plot(I, '-', label='Infected')\n",
" plot(R, ':', label='Recovered')\n",
" decorate(xlabel='Time (days)',\n",
" ylabel='Fraction of population')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"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": "markdown",
"metadata": {},
"source": [
"Here's an updated version of `run_simulation` that uses `unpack`."
]
},
{
"cell_type": "code",
"execution_count": 5,
"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": "markdown",
"metadata": {},
"source": [
"**Exercise:** Write a version of `update_func` that uses `unpack`."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# Original\n",
"\n",
"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": 7,
"metadata": {},
"outputs": [],
"source": [
"# Solution\n",
"\n",
"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",
" unpack(system)\n",
" s, i, r = state\n",
"\n",
" infected = beta * i * s \n",
" recovered = 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": "markdown",
"metadata": {},
"source": [
"Test the updated code with this example."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" S \n",
" I \n",
" R \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0.988889 \n",
" 0.011111 \n",
" 0.000000 \n",
" \n",
" \n",
" 1 \n",
" 0.985230 \n",
" 0.011992 \n",
" 0.002778 \n",
" \n",
" \n",
" 2 \n",
" 0.981296 \n",
" 0.012929 \n",
" 0.005776 \n",
" \n",
" \n",
" 3 \n",
" 0.977071 \n",
" 0.013921 \n",
" 0.009008 \n",
" \n",
" \n",
" 4 \n",
" 0.972541 \n",
" 0.014970 \n",
" 0.012488 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plot_results(results.S, results.I, results.R)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sweeping beta"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Make a range of values for `beta`, with constant `gamma`."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.25"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"beta_array = linspace(0.1, 1.1, 11)\n",
"gamma = 0.25"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run the simulation once for each value of `beta` and print total infections."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.1 0.0072309016649785285\n",
"0.2 0.038410532615067994\n",
"0.30000000000000004 0.33703425948982\n",
"0.4 0.6502429153895082\n",
"0.5 0.8045061124629623\n",
"0.6 0.8862866308018508\n",
"0.7000000000000001 0.9316695082755875\n",
"0.8 0.9574278300784942\n",
"0.9 0.9720993156325133\n",
"1.0 0.9803437149675784\n",
"1.1 0.9848347293510136\n"
]
}
],
"source": [
"for beta in beta_array:\n",
" system = make_system(beta, gamma)\n",
" results = run_simulation(system, update_func)\n",
" print(system.beta, calc_total_infected(results))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wrap that loop in a function and return a `SweepSeries` object."
]
},
{
"cell_type": "code",
"execution_count": 12,
"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": [
"Sweep `beta` and plot the results."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" values \n",
" \n",
" \n",
" \n",
" \n",
" 0.1 \n",
" 0.007231 \n",
" \n",
" \n",
" 0.2 \n",
" 0.038411 \n",
" \n",
" \n",
" 0.3 \n",
" 0.337034 \n",
" \n",
" \n",
" 0.4 \n",
" 0.650243 \n",
" \n",
" \n",
" 0.5 \n",
" 0.804506 \n",
" \n",
" \n",
" 0.6 \n",
" 0.886287 \n",
" \n",
" \n",
" 0.7 \n",
" 0.931670 \n",
" \n",
" \n",
" 0.8 \n",
" 0.957428 \n",
" \n",
" \n",
" 0.9 \n",
" 0.972099 \n",
" \n",
" \n",
" 1.0 \n",
" 0.980344 \n",
" \n",
" \n",
" 1.1 \n",
" 0.984835 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"label = 'gamma = ' + str(gamma)\n",
"plot(infected_sweep, label=label)\n",
"\n",
"decorate(xlabel='Contacts per day (beta)',\n",
" ylabel='Fraction infected')\n",
"\n",
"savefig('figs/chap06-fig01.pdf')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Sweeping gamma"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the same array of values for `beta`"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1])"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"beta_array"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And now an array of values for `gamma`"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[0.2, 0.4, 0.6, 0.8]"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gamma_array = [0.2, 0.4, 0.6, 0.8]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For each value of `gamma`, sweep `beta` and plot the results."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Saving figure to file figs/chap06-fig02.pdf\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXd4XFmZ5/+pqFLOOVhyOs6Ws7tttzu57aYDDZMWGHZ3fsMzMPwGJrDzW3qBIQ08MDA9s7DAsEwzMzDApIWls7vbNO3U7SxnH0u2ck4lVUmqfH9/XKmCLEtlWxWkOp/n0aM6556q+6pUdb/3vOc972vQNA2FQqFQKJINY6INUCgUCoViJpRAKRQKhSIpUQKlUCgUiqRECZRCoVAokhJzog2YL4QQacA2oBvwJ9gchUKhUESHCSgHTkkp3eEHFo1AoYvTkUQboVAoFIq7Yg9wNLxjMQlUN8BPfvITysrKEm2LQqFQKKKgp6eHD33oQzB5DQ9nMQmUH6CsrIyqqqpE26JQKBSKO+OWpRkVJKFQKBSKpEQJlEKhUCiSkri6+IQQfwT8V2A98DMp5X+dZeyfAv8dSAf+D/CH0yM8FAqFQrF4ifcMqgv4S+CHsw0SQuwHPg08AtQCS4Evxto4hUKhUCQPcRUoKeXPpZT/FxicY+h/AZ6XUl6WUg4DX0afeSkUCoUiRUjWNai1wPmw9nmgVAhRmCB7FAqFQhFnkjXMPAsYCWtPPc5m7tmXQqFQJBRN0wgENAKaRkCDQEAL9vmn+gMa2uSxQNgxLew5gfDXCR8f0PBrGmigofcDBCb7AtpUn/5bm7RJ04h4jqZpk8fCxmrhz9H7A5MNDdACoedoaJQVZrJlVSkmo2He38dkFSgnkBPWnnrsSIAtCoViAaJf8AN4fQF8fg2fPxD88U+29WN62zt13Dc1Tps8FtCP+SLbEeIyJSQB/WKeSnX22nocFOelU1eRO++vnawCdRnYCPzbZHsj0CulVLMnhWIRo2kaHl8Al9uHy+Of/O3D5fbj8viY8Pjxev0RAhIUlYA2+Vtv+wOpIxKJxGQ0kJeVFpPXjneYuXnynCbAJISwAT4ppW/a0B8B/yiE+Al6+ovPAv8YT1sVCsW9EQhourgEhcYfITYuj48Jtx/35JgJtw+3x6+7qRYBRqMBo8Ew7fe0/snHBoN+oTdMHz/1nLC2wWDAZDBgMILBYMDA5G8DwcdMPgcmfwePgYHJ5xLebwgdMzBj/1Sf0Rj2ukBxfjoZNktM3sN4z6A+C3w+rP27wBeFED8ErgBrpJRtUsrXhBB/BbxFaB/U5295NUXUdHR08Oyzz3LhwgXKy8v5i7/4C+6///4Zx37961/n0KFD9Pf3U1paysc+9jGeeeaZOFusSCY0TWPM5WPC5QuKS0hoIgVoSmjc3sQWFTAYDJiNBsxmI6bJ3xaTEbPJiNlsDB4zm8J/JvuMk2NMhtAxc2iMyWiMFBCjISQwUxd1w/yvyaQacRUoKeUXgC/c5nDWtLHPAc/F2KSU4VOf+hT19fX84Ac/4O233+aTn/wkr7/+OgUFBbeMTU9P53vf+x51dXVcvHiRj3zkI9TU1LB58+YEWK6IJy63D7vTjd3hZtjhxu50M+J0M+Jw4/UH4mKDxWwkPc1MmtWEzWqe/DGRnmbGlmbCajZFiEWEwISJypRgKBYuyboGlTK8//3vp7Kykm9961vBL9M3vvEN7HY7X/nKV+blHM3NzVy+fJnnn38em83G/v37+ad/+icOHjzIBz7wgVvGf/KTnww+3rhxI1u2bKGhoUEJ1CLB6/Njd3iwO12MOD3YHS6GHW5GnB5cnune9rvHYDCQZjFhs5qwpZlJt5pIs+oiMyU6tjSzLjyTx9KtJkymZN39oog3KSlQ52QfJ6/04PXF7o7QYjayfU0Zm0TJrON+8pOf8Mwzz3D69Gm2bduG3+/nhRde4LnnZp48fvSjH+XMmTMzHtuyZQvf//73b+lvamqiurqarKzQJHXVqlU0NTXN+Xe4XC4uXbrEBz/4wTnHKpIHvz/A6LgHu0OfDU3NikacbpwT3rt6zTSriSybBVtaSFyCQhMuPGkm0q1mrBZTcL1CobgbUlKgGq73x1ScALy+AA3X++cUqPT0dJ5++mkOHjzItm3bOHr0KBaLha1bt844fiYBmouxsTGys7Mj+rKzs+nt7Z3zuZ///OcRQrBnz547Pq8itmiahmPcy4gzUoTsTjeOMc9dBRtYTEbystPIzUojL3vyJ0v/saWl5OVCkUBS8hNXv7I4LjOo+pXFUY3dt28ff/AHf8BnPvMZfvGLX/De9753Xn3nmZmZOJ3OiD6n00lmZuasz/v6179OY2MjP/rRj5QvP4F4fQH67ePBGdDUrGhkzIPvLtaFjAYDOVlW8rPSyMu2kZtlnRQjG5k2s/pfK5KGlBSoTaJkzplNPFm5ciVWq5WjR4/yq1/9ihdffPG2Yz/ykY/M6uL7+7//+1v6ly9fTnt7O06nM+jmu3btGk8++eRtz/Otb32LI0eO8OMf/zjCNaiIDyNON209Dlq6R+nsd96xEBkMBrLSLREzoKmZUU6mVbneFPfMhNeF3TVKTloWmdaMmJwjJQUqGXnkkUf43Oc+x5o1a1iyZMltx80kQHNRV1fH6tWr+c53vsOf/MmfcPjwYaSUfPvb355x/Pe//31eeukl/vmf/5n8/Pw7Pp/izvEHNHoGx2jpHqW1e5ShUVdUz0tPM5Mf7pILEyKzCjZQzANun4dh1whD43aGXXaGJ0YYmhjB5dU/owaDgd9Y+x4K0vPm/dxKoJKEffv28cMf/pCPfexjMXn95557jmeffZZt27ZRXl7Ot771rWCI+QsvvMD3v/99Xn755eBYi8XC/v37g8//6Ec/GjPbUpVxlzc4S2rvdcy6b6ggx0ZhbvqkGFmDrjmbVX2FFfOD1+9l2DXK8MSUCNkZmhhh3DM+6/M0TcPhdiqBWszU19cDsGPHjpi8flVVFT/+8Y9nPPb000/z9NNPB9tSypjYkOpomkb/8AQtPfosqW944rY528wmI5XFWdSW57CkPIecTGucrVUsVvwBP3bXaJgI6YLkcDvnfnIYZqOZ/PRcavOrqM6piImtSqCSBKPRGPFbsTjweP209zpo7RmlpdvBuOv2Id7ZGVaWlOdQW55DZXEWFrP6LCjunoAWYMTlCArR8MQIwy47Iy7HHSWzNRqN5NlyKUjPJT89l4L0PPLT88i2ZsY8oEYJlEIxzww7XLR264LUNeAkcJukpUaDgbLCzMlZUjYFOTYVQae4YzRNw+EZu2WNyO4aIRCIPrjGYDCQa8sh35ZLQUae/js9lxxbNkZDYm6WlEAlEcq1tjDx+wN0DYQCHOxO923H2qxmlpRls6Q8h5qybLWGpLgjptZ7uhx99Dj7GJqwY58YxRe4gwwgBgPZ1szJmVBu8HeeLQeT0RQ74+8C9e1QKO4C54SX1u5RWnv0AIfZ9tQV56UHXXcl+RkqxFsRNZqmMep20O3oo8vRR5ejd86ghXAyrRkRIlSQnkeeLQeLKTbZx+cbJVAKRRQEAhp9w+PBWVK/feK2Yy1mI1Ul2cEAh6z0hXExUCSeKUHqcvQGRSkaQbJZbJMilBNcI8q35ZJmXtjBNUqgFIrboGkaLd2j3Oiw09rjYMJ9ezdKblYatWX6WlJlcZZKeKqICk3TGHE76Hb00jXaS7ezj3HP7W9+ACwmC2VZxZRnl1KSVUi+LZd0iy1OFscXJVAKxQwMjbo4fK6Djr6ZQ2+NBgMVxVnUluvrSXlZaSrAQTEnmqYx4hqly9Gni5Kjlwnv7JuyLSYL5dkllGeXUJFdSmFGfsKCFuKNEiiFIgyvz8/JK72cv95/S7LVDJslFOBQmo3VklwLyorkQ9M07K7RSXddL91RCJLVZKFsUozKs0tSSpCmowRKoUC/kDR12Dl2viuiHIXRYGDdskJWLSmgOD9dzZIUszIlSKE1pN5gSqDbYTVbKc8KzZAKMvJSVpCmowQqRbiTku9T2O12Dhw4QF1dHT/72c/iZGn8GR518fa5Tjr6HBH9FUVZ7N1cSWFueoIsUyQ7mqYx7BqhO+iy64takCpy9BlSQboSpNuhBCpFuJOS71N885vfZNmyZXe02W8h4fX5OXWll4bG/ojNtBk2C7s2lLOyJl/NmBS3MOp20mbv1EXJObcgpZnTKM/WgxoqskspSM9Tn6soUQKVYJKx5DvAuXPnaGxs5Ld/+7f5j//4j3mxI1nQNI0bnSMcbeiMcOcZDAY2LC9i+9oy0tT6kiIMt8/DzeFWGgeb6XH0zzpWF6QSKrJLKFeCdE+kpEBd6LnK6a4L+Px3sPv6DjGbzGyt2MCGstWzjkvGku9+v58vfelLfPnLX+b69euz2r/QGHa4OHyuk/beSHdeeWEmezdXUZSn3HkKnUAgQPtoF9cHmmkd6bitJyEkSPoMKT89VwnSPJGaAtV7NabiBODz+7jQe3VOgUrGku8//vGP2bBhA+vWrVs0AuX1+Tl9tY9z1/si3HnpaWZ2bahALFHuPMVkxvnxIRoHb3JjsBWX79a0VQaDgaqccqpzK6jIKSXfpgQpVqSkQG0oXR2XGdSG0tnFaYpkKvne29vLj370I37+85/P2/kTiaZp3Owc4ej5LhzjnmC/wWBg/bJCtq8tU/nwFDjdYzQONdM42Ix9YnTGMUWZBawsrGNZQe2i3RibbKTkN3ND2eo5ZzbxJJlKvl+8eJH+/n6eeOIJAFwuF263m127dnH48GFMpoWzNmN3uDnc0EFbT6Q7r6wwk72bqijOV+68VMbj93JzqI3GwWa6HTN7EzKtGaworGNFYR356blxtlCRkgKVjCRLyfcHHniAX/3qV8H2K6+8wksvvcR3v/vdBSNOXl+AM9d6OSf78E9z592/voJVtcqdl6oEtAAdI900DjbTYu/AH7i1irHZZGZpfg0rCuuoyC5Vn5UEogQqSUiWku9Wq5Xi4uLg87KzszGbzRF9ycpU7rwjDZ2MjkW689YtLWTHOuXOS0U0TWNwYpjGwWaaBltmzuRgMFCVU8aKwjpq86oWTLbvxY7hTiorJjNCiFqg+dChQ1RVVSXanDsmEAiwevVqXnvtNerq6hJtzoJjxOnm8LlOWnsi1w9KCzLYu6mKkoKMBFmmSBRjnnEaB1toHLzJ8MTIjGMKMvJZWVjH8oJaMqzK5ZsIOjo6eOSRRwDqpJQt4cfU7WSSoEq+3x0+f4Cz1/o4c603wp1ns5q5f0M5q2sLlIsmhfD6vTQPt9M42EynoxdmuAFPt9iC60qFGfkJsFIRLUqgFAuW5q6RGd15a+sK2LmuHFua+ninAgEtQJejl8aBZprt7TNG55qMJurya1hRWEtlTplKLbRAUN/gJEKVfI+OEaebow2dNHcrd14qMzRh5/rATZqGWmcu6mcwUJFdysrCOmrzq7GqdaUFhxIoxYLB5w9wVvZx5uqt7rz71pezpk658xY7494JmgZbaBxsZnB8eMYxeem5+rpSYS1Z1lv3+ikWDkqgFAuClu5RDp/ruMWdt6augPuUO2/R4/K5Odd9ict912dMOWSz2FheUMuKwlqKMtSNymIhrt9qIUQB8DzwGDAAPCul/OkM49KA/wm8D7AAx4CPSSk742iuIgkYcbo5er6L5q7IKKyS/Az2bq6iVLnzFjW+gJ9LvZKG7kt4/N6IY0ajkdq8alYW1lGVU64CjBYh8b7t/A7gAUqBeuBlIcR5KeXlaeP+GLgP2ACMAD8Avg28P462KhKIP6Bx9lovZ6714fOH7pjTrCbuW1fOmrpCjEZ1l7xYCWgBGgebOd15gbFp60tFmQWsKV5BXX4NaWZrgixUxIO4CZQQIhP4DWCdlNIJHBVCvAB8GPj0tOF1wEEpZe/kc/8FmDm9t2LRoWkavzrVhmyLXGNYU1fAfesrSFfuvEWLpmm0j3RxouPcLXuXcmzZbK+spy6/WrnwUoR4ftNXAn4pZXh67PPA3hnGPg/8TyFEBWAHPgS8GnsTFcnAleahCHEqzktn7+YqygrVgvdipm9skBPt527Ji5dusbGlYj2ripYrN16KEc//dha6uy6cESB7hrHXgTagExgFVgNfiql1i5yOjg4+/OEPs3HjRg4cOMDx48dnHX/8+HHe9773UV9fz969e3nllVfiYueAfYIjDaGlxtW1BfzWIyuVOC1iRl0O3rxxlP975bUIcTKbzGyp3MB/Wv80a0pWKnFKQeI5g3ICOdP6cgDHDGO/B9iAQmAM+P/QZ1A7YmngYuZOSr43NTXxqU99iq997Wvs2rULh8OBwzHTv2l+8fr8vPZuS3DNqTDHxgObqtRa0yJlwuvibPdFrvQ1Ep5yzWAwsLp4BZsr1pFhUemHUpl43pJcB8xCiBVhfRuB6QESU/3/KKUcklK60QMktgshiuJgZ1x5//vfzyc+8YmIL+g3vvENPvOZz8zbOaZKvn/iE58IlnxfuXIlBw8enHH89773PX7nd36HvXv3Yjabyc/Pp6amZt7smQlN03j7bAd2h14gzmIysv++Wixmdde82PD6vZztusi/XHyBy73XIz77dfk1/Na6J9m9ZJsSJ0X8ZlBSyjEhxM+BLwkhPoIexfde4P4Zhp8C/rMQ4tfAOPBxoEtKOTAfttgbzjN08jQBn3fuwXeJ0WyhYPtW8uo3zjouGUu+NzQ0UF1dzVNPPcXw8DA7d+7ks5/9LHl5ebP+LffCtZZhrrWG1p32bq6iIEcVhVtMBLQA1/pvcKbrwi0ZxcuyS9hRVU9pVvJnzVfoeEcdOJuaMKWnk71KxCRwJd7hUB8Hfgj0AYPAH0opLwsh9gCvSimnrqD/DfgW0AhYgUvoe6LmBXvD+ZiKE0DA58XecH5OgUrGku+9vb288MILPP/885SUlPDpT3+aL3/5y/z1X//1HZ87GgZHJnj7XEewvWpJAatqb3U9KhYmmqbRau/gZGfDLdVq89Jz2VFVT01upYrMWyC4+wewNzTgbLoRnP2aMzPIiIGXJa4CJaUcAp6Zof8IehDFVHsQPXIvJuTVb4zLDGoucZoimUq+A6SlpfH+978/WPbjox/9KL/3e783b/aE4/X5Ofhua3DdqSDHxt7NlTE5lyL+9Dj7OdF+jl5nf0R/hjWDrRUbWFlUpxK3LgA0TWOivQN7w3nGOzoijhkMBoxpsfF2pOSGkrz6jVGLRzxIppLvAEKIu/gr7o4jDZ0MjeruHrPJyP6dS7CYF0blXsXtsU+McLKzgZbhyIuZxWShvnwN60tWYTal5OVnQaH5/TibbmBvOI97cPCW4+mVFRRs3YqttCQm51efkCQhWUq+gx648d3vfpf3vve9FBUV8YMf/IAHH3zwjs87F7J1iCvNQ8H2A5sqKcxVC+MLmXHPBKe7LiAHbkQEPxgNRtaUrGRz+VpsFrW2mOwEPB5Gr1zFfv4CvrGxaUcNZC1fSl79RmwlsRGmKZRAJQnJUvId4Dd/8zfp6urit37rtwDYs2cPn/3sZ+fVnuFRF78+G7q7FjX5rFbrTgsWj9/L+Z4rXOy5hi8QWY9pWcEStlXVk5OWdZtnK5IF39gYIxcuMXL5MgGPJ+KY0Wwme9Uq8uo3YMmZvmMoNqiS70lCKpV89/kD/PuhRgZHJgDIy07jtx9ZidWiXHsLDX/Az9X+Js52X8I1LTKvIqeMHVX1FGcWJsg6RbR4hoaxN5zHcb0RLeCPOGZKTyd3/Tpy163FZJv/2e9dlXwXQnww2hPMlJFccWekUsn3Iw2dQXEyGQ0c2FmrxGmBoWkaN4fbONXZwKgrMgCnICOfHVX1VOWUq8i8JEbTNFzd3djPnWestfWW45bcXPLqN5AtBEZzYpxts53169PaBUA6ocwP2cAEeri4EihFVFxvG+byzdBi6576Sory1LrTQqLL0cuJ9nP0j0UummelZbK1YgMrCuuUMCUxWiDAWHML9nMNuPr6bjluKy0lb9NGMmtrMST4hvm2AiWlrJ56PDmb+ih6Taark32r0VMS3fmqvWJGFnvJd7vDzVtn2oPtFdV5rF2q3D8LhaFxOyc6ztE+0hXRbzVb2Vy+jjUlKzEb1Uw4WQl4vTiuSeznL+AdHb3leGZtLXmb6kkvL0uAdTMT7bztq8DTU+IEIKW8KoT4E+AF4J9jYZxi8eDzBzj4bgten77fKTcrjYe2qLIJCwFN0zjddYFz3ZchPDLPaGR96Srqy9aqukxJjH9igpFLlxm5eAm/K3Kd0GA0kS1WkFe/EWt+foIsvD3RClQRuotvOnlA8v1ViqTj2Pku+u2hdaf9O5eodacFgMfv5a2bx2m1h+1nMhhYWVjH1ooNZKWpLPPJindkBPv5C4xevYbmjwx8MFrTyF23htwN6zFnJG9V6mgF6l+BHwkhPoOeJ08DtgN/CfxbjGxTLBKaOuxcvBFKo7hrYwUl+cn7pVDoONxODja9zdC4PdhXkVPKfdVbKMxQ96XJiqu3D/u5Bpw3m9Ev1SEs2dnkblhPzprVGC2WxBh4B0QrUH8IfA74G0IzqSHg74Avx8AuxSJhxOnmV6dD607LKnNZv2zRJaVfdHQ7+nij6TAunzvYt6FsNdur6lVqoiRE0zTGW1uxnzvPRHf3LcfTiorI21RP1rKlCQ98uBOiEigppQddoD4nhCgEDPOVWVyxePH7Axx8txWPV3cv5GRaeWirWndKdq72N3K09VQwE4TRaGTPku2IomUJtkwxHc3vx3G9EXvDeTzDw7ccz6iuJm/TRtIrF2Yy3jsKbhdCbAaWAq9MtjMBt5TSN+sTFSnJ8Yvd9A2PA2A0Gti/sxabVSUvSVYCWoB3289yqTcUTWqz2Hhs+QOUqTIYSUXA49EDHy5cxDc+HnHMYDCQtUIPfEgrWthRslFdLYQQlcAvgbWABVgJ3AS+AQSAP4qVgYr5oaOjg2effZYLFy5QXl7OX/zFX3D//TOV4gK73c4XvvAF3nnnHQB2797NF7/4xYh6UnNxs3OE842hDNa71ldQWqDWnZIVt8/DmzeO0DnaE+wrzMhn//K9KhAiidA0DWfTDQaPv3NLjjyj2ULO2tXkbdyA+Q6+q8lMtM7I7wAXgVz0zblT/Buwf76NUsw/n/rUp1izZg0nTpzgT//0T/nkJz/J0NDQjGP/9m//ltHRUQ4dOsSbb77J4ODgbRPLzsTomIdDp9uC7bqKXDasUOtOyYp9YoRfXH0tQpzq8mt4etU+JU5JhMdup/vFl+l9480IcTJnZFC4cwdL/svvUrTr/kUjThC9QO0Fvjq5FhVOC7DwEt8lEclY8n0qN1ZWVhbZ2dns27fvttV3p+Of3O/k9ujrTtkZVh5R605JS/tIF//36kFGXY5g3+aK9Ty6bDcWU/JHeaUCAa+XwXdP0P4v/xZRi8mUnk7Jg3tZ8uEPkb95E6a0tARaGRuiXRDwAjPdSi1HT3W0oLgh+2m80ovP55978F1iNptYsaaUZWJ2330ylnz/0Ic+xE9/+tNgvaiDBw/y8MMPz/p3TPHu5R56hybXnQz6fidbmlp3SjY0TeNi7zXe7TgX3HxrMpp4qO4+lhbcvtyLIn5omsZYcwsDR4/hiyg4aiB3/ToKtm9dlKIUTrRXjv8DfFkI8TuTbU0IsQz4JvoeqQXFzev9MRUnAJ/Pz83r/XMKVDKWfF+zZg1er5cdO3YAcN999/HBD86dO7ile5RzMpTba+f6csoKlYso2fAH/BxpPcn1gZvBvkxrBvtX7KUoQ5U8SQa8IyP0HznGeFtbRL+ttJTiB/aQVpwaLvNoXXyfQk8SOwBkACeB60AzMH++qDixdGUx5hhXbTWbTSxdGV3k0759+3jzzTfRNC0pSr7/8R//MbW1tZw9e5YzZ85QXV3Nn//5n896Due4hzdPhr5MS8py2BTl36+IH+PeCV6ShyLEqSSriPetOaDEKQkI+HwMnTxF28/+LUKcTDYbJQ89SOX7n0kZcYLo90GNAx8UQiwF1gFZwHkp5eVYGhcrloniOWc28STZSr5LKfnCF75AxmQKlA984AOzzqD8AY2D77bi8ui7DbLSLTy6vUatOyUZA+NDHGx8mzFPKCx5ZdFS9izZjkkleU04Y61tDBw+gtfhCOs1kLNmNYU7t8ekFlOyE22Y+f8A/kZKeRM9vHyq3wb8mZTyqzGyL2VIppLv69at49///d+Ds6Z//dd/RQhx29c/ebmb7kE9qshoMPDYziWkq3WnpOLmUBu/bn4nVO3WYGBn1SbWl65SNxIJxjvqYODYMcaaWyL604qLKX5gD7bS2JZVT2aidfF9Gb3+03SygC/Nnzmpy759++ju7uaZZ56Jyes/99xzXLp0iW3btvHNb37zlpLvTzzxRHDsV7/6VTo7O9m7dy8PPPAA7e3tfO1rX5vxdVt7RjlzLbTutH1tGRVFiyfMdaGjaRpnui7y5o0jQXGymCwcWL6XDWWrlTglEM3vZ/jMWdp+9i8R4mS0plG8dw9Vv/G+lBYnmGMGJYSomHxoAMqEEOHjTcA+4NaKV4o7pr6+HiAYmDDfVFVV8eMf/3jGY08//TRPP/10sF1dXc3f/d3fzfmazglvxLpTTVk2W1al9hcqmfD5ffy65R1uDoX+Rzm2LPYvf5D89NwEWqYYb++g//ARvCMjEf05q1ZRsHMH5gxVxBPmdvF1oKfD1YBz044ZADcw++q5IioWWsn3QEDj9XdbmXDrd+WZNguPblPrTsmC0zPGwca3GRwP5WeryCnj0WW7sZkXd2hyMuNzjjFw/DjOphsR/WmFhRQ9sCepigUmA3MJ1Ap0IboO3IcexTeFF+iZYfOuIgU4eaWHrgE9MtAwue6UYVMbO5OBXmc/rzcdZsIbKk63tnQl91VtWTA3QIsNze9n5OIlhk6eJuDzBvuNFgsFO7aTu27tgsoyHi9mFSgp5Q0AIYRFShnbjUOKBVPyvb3XEbnutKaUymK17pQMXB+4yeHWEwQCeuVig8HArpptrClZkWDLUpeJri76Dx/FMy21WPbKFRTetxPzbbZ7KKLfqPsFIUSrlDIihEwI8RGgWkr5+fk3TZGMjLu8vHGyLZiaqaokmy31HqINAAAgAElEQVSrShNslSKgBTjZ0cCFnqvBvjRzGvuW76EiW/1/EoFvfJzB4+/guN4Y0W/Nz6f4gT2kV1bc5pmKKaIVqN8D3jtD/xng85M/ikVOIKDx+ok2xl26iyLDZuGxHTUYjWrdKZF4fB4O3TxG+0hXsC8/PZf9Kx4kJ03NbOONFggwcukyQydOEvCGufPMFvK3bSFvw3oMJrXvLBqiFahCYHSG/jEgdbY1pzinr/XS0advIjQYDOzbXqPWnRLMiGuUg02HsU+EosFq8ip5eOkurCrZa9xx9fTQf/go7oHIeq5Zy5Yuukzj8SBagboM/DbwlWn9vw1cm1eLFElJZ7+TU1dCufu2riqhunSmrXGKeNE52sMbN47g8YXilOrL17KtcqOKpowz/okJBt85wei1yMuhJTeX4j27yaipTpBlC5toBeqzwC+FEOuBw+hh53uB9wHvj5FtiiRh3OXl9Xdbg+tOlcVZbFujwmETyeW+6xxvOx1Rln1v7U5WFNYl2LLUQtM0Rq9cYfCdkwQ87mC/wWQif8tm8jfVK3fePRBtLr7XhBBbgf8O/AF66PklYJuU8kIM7VMkGE3TePNkG2OT607paWb27Vii1p0SRCAQ4Fjbaa72hxbeM6zpPLZ8LyWZC7u890LD1dfHwOGjuPoicxVk1tZStHsXlhzlYbhXok6YJqW8CPzuvZxMCFEAPA88hr6n6lkp5U9vM3Yz8LfAZvS1rq9KKf/nvZxfceecudZHW28oeeWj22vISldrG4nA5XXxxo2jdDtCrtaizAL2L99LpjUjgZalFn63m6F3TzJy+Qq6M0nHkp1N0Z5dZNbWJsy2xUbUAiWEKAc+ANQBX5BSDgoh7ge6pZTNUb7MdwAPUArUAy8LIW7Jii6EKAJeA/4U+A/AiqrcG3e6+p2cuBwqA75lVSlLynISaFHqMjRh52Dj2zjcobIpywqWsLfuPswqE3nccPX20vPa6xEl1w1GE/mb68nbvAmjWSVJnk+izWa+C3gFOA3sBv4GvZLuQ8BG9GCJuV4jE/gNYJ2U0gkcFUK8AHwY+PS04X8GHJRS/mSy7QauoogbE24fr58IrTtVFGWyY61ad0oEbfZO3rx5FJ/fF+zbVrWR+rK1KhgijoxevUb/20fQAqGcBRnV1RTt2Y01T+U2jAXRyv1fA5+TUn5LCBFerOR14ONRvsZKwC+lvB7Wdx492GI6O4GLQojj6GXlTwD/r5SybYaxinlmat3JOaGvO9msZh5T604Joc85wOtNhwloemYIs8nMw3W7qM1XDoV4ofn9DBx7h5FLl4J9RmsaJQ/tJXNpnbpJiCHRCtQ64IUZ+gfQ90hFQxYwMq1vhJnLeFShrz3tAy4CfwX8DNgV5bkU98C56/209oS2vT26vYasDGsCLUpNJrwuXr9xJChOWWmZ7F++l8KM/ARbljr4xifoPfg6E93dwT5rQQHlj+/HkqtmTbEmWoHqBZYBLdP6d6OXfY8GJzB9ASMHvZT8dCaAX0gpTwEIIb4IDAghcqWU00VOMY94vH5Ohq07bRIl1Jardad4E9ACvHnjKOOT1W+tZitPikdVZog44urro+fVgxHrTVlLl1LyyEMYLSpQKB5Emz73fwHfEUI8PNleJoT4feA59Ei7aLgOmIUQ4VkrN6JvAp7OBcLDY0KP1Vw6xjR3jeDz63fsBTk2dq4rT7BFqcnJjoZQtJ7BwMN19ytxiiMOeZ3OX/wyTJwMFO7YTun+fUqc4ki0+6D+RggxBvwDkAkcBPqBr0gpvx/la4wJIX4OfGkyyWw9en6/+2cY/g/A/xFCfAtdwD4HHJVS2qM5l+LuaWoPvcViST4mte4Ud24OtUYkfd1SsZ6avMoEWpQ6aIEAg++8i/18aHun0WqldN+jZC6pSaBlqcltBUoI8R7gDSmlF0BK+b+B/y2EyAGypJRdt3vuLHwc+CF6Fd5B4A+llJeFEHuAV6WUWZPn+pUQ4n8ALwMZwFHgg3dxPsUd4PL4IvY8La/KS6A1qcnwxAi/bnk32K7Jq2Rz+boEWpQ6+F0ueg6+wURnZ7DPmp9P2eP7seap70IimG0G9SJQBvQLITxApZSyX0o5ysyJY+dESjkEPDND/xH0IIrwvu8B37ub8yjujpauUfwB3Ztakp9BbpaqvBpPPH4vrzcdDoaT59iyeKjufhUlFgfcA4P0vPoaXkfoBi2ztpbSRx/GaFUBQoliNoHqAZ4WQrw2Oa5cCDGj8/UuZ1OKJKMxzL2nZk/xRdM0ft38DiMu/d7PZDSxb9kDpJnVxTHWOBqb6H/r1wR8oX1mBVu3kL9tq7o5SDCzCdTngW8DaehBCudmGGOYPKa2si9wXB4f7eHuvWolUPHkfM8VWobbg+29tTtVOHmM0QIBhk6cYvhc6NJmNFsoefRhspaqpLvJwG0FSkr590KIHwEVwE3gPvR9T4pFSHPnKIHJrBGlBRnkZKo793jRMdrNyc7zwfa60lUsL6xNnEEpgN/tpvf1NxlvD90UWHJzKX/8ANYCdWOQLMwaxSel9AAtQoh9wBkppW+28YqFS2PHcPCxcu/FD4fbyaEbx2Dy5qAsu4SdVZsSbNXixj04RM9rB/GOhLZUZtTUULrvEUxpat01mYg2zPyQECJXCLEFKGHa/qnbZSRXLAxcbh8dvWFJSJVAxQVfwM8bN47g9ul1hDKs6Ty6bDdGY7TbExV3ivPmTfrefIuAL1SKPX/zJgq2b8Og3vekI9pkse8D/gmwoIeHT99EqwRqAXOjcyTo3isrzFTuvTigaRrHWk8xMDYEgMFg4NFle8iwpCfYssWJpmkMnzrN0OkzwT6j2UzJww+RtXxZAi1TzEa0qY6+iR4w8Xnl5lt8NHWEovdWqNlTXLg20IQcuBFs31e9hbKs4gRatHgJeDz0vnGIsdbWYJ8lO5uyxw+QVqSKPCYz0QpUCfBDJU6Lj3GXl86+cPeeSoAZa/rGBjnWdjrYXl5Yx9qSlQm0aPHisdvpeeU1PPbQTVhGVRWljz2KyWZLoGWKaIhWoP4VvQqu2ji7yGjuCkXvlRdmqqzlMWbC6+KNpsMEApP5DjPyeWDJdrXfJgaMtbTS+8abBLyh9aa8+o0U7tyh1pvmCU3T0DRiVoonWoHqBv5SCPEIevkLb/hBKeVX59swRXyI2Jyr9j7FlIAW4NDNo4yFZSh/bNkezCZVhXU+0TSN4TNnGTp5mqnlcoPJRMlDD5K9csXsT1ZEhaZptN4c5PrlXswWE/ftXUp6DG5uo/1mPAhcQS/VXjrtmAYogVqAjLu8dPbr7j2DwaCi92LMyY4GukanZSi3zVQOTXG3BLxe+g69hfPmzWCfOSuL8sf3k1as1vjmA6/Hx/nTHfR06mH6HreP4cHxxAmUlHLPvJ9ZkXBudI4ES7qXF2aSla7KCMSKm0NtKkN5jPGOjND96kE8Q0PBvvSKCkof24c5Q0VHzgdDA2OcO9HGxLgn2JeTl05pRWxqxinfQgoTXlpjhXLvxQw9Q/k7wbbKUD7/jLe10fP6IQIed7Avb8N6Cu/bicGkMrHdK1pAo0n2cf1yb/CmFqB2WRGrN5ZjMsVmTW+2chuvA78lpRwRQrxB5N6nCKSUj8XCOEXsGHd56RrQi7Hp7j0VvRcLVIby2KJpGvaG8wy+c4LgepPRRPHePeSsXpVY4xYJrgkv5062MRgW7WuxmtiwpZryGF83ZptBnSAUDPHuLOMUC5AbHSH3XkVRJhk25d6bb1SG8tgS8Hrpe+ttnE1NwT5zZiZlB/ZjKy1JoGWLh77uURpOteNxh2V6L8pk046amKw5TWe2ZLGfm+mxYnGgovdij8pQHju8ow56XjuIeyCUv9pWVkbZgccwZ2Qk0LLFQcAf4NqlHm5e7w/2GQwGlq8uYeXqUgxxqrSt1qBSEOeEl+7BMPdepXLvzTe3ZigXKkP5PDHe0UHvwTfwu0PrTTlr1lC8Z5dab5oHxhxuzp5oY2R4PNhnS7ewaUcNhcVZszxz/lEClYLc6LAH3XuVxVnKvTfP3JqhvJidVZsTbNXCR9M0Ri5cZPD4O8HPr8FopGjPbnLXrkmwdYuDjtZhLp3txOfzB/tKy3PYuK0aa1r85UIJVAqiovdix/QM5ekWG48u3aMylM8D05O9mtLTKTuwn/TysgRatTjw+fxcOttFR2soRN9oNLJqfRl1K4oSFtSjBCrFcI57gu49o8FAXYz2L6Qq0zOU71v+ABlWtQfnXhk+1xAhTraSEsoO7MeclZlAqxYHI8MTnD3Rypgj5DLNzEpj884acvMTu56nBCrFuNERKtJWWaLce/PJ1f5GlaE8BoxcuszgO6FA4oyaGsoOPIbRrC5f94KmabQ0DXL1QncwNyRA1ZJ81m2qxGxJ/HpetPWgsoD/Duxl5oKFKhXzAqExrLSGqpw7f6gM5bHBIa/Tf/hIsJ1eUaHEaR7wuH2cP9VOb/dosM9sNrFucyVVS5In0jTa//I/ANuA59ETx952064ieRkd89AT5t5T0Xvzg8pQHhucN27Se+itYNtWUkL5ew4ocbpHBvudnDvRhmsilPM7Nz+dzTuWkJmdXCXvo/1PPwY8IqU8PedIRdJyI2z2VFWahS0BUTmLDT1D+TGVoXyeGWtto/eNN5m6F04rLKT8qScwWtUm57tFC2hcv9pL09W+iHRFdSuKWb2+DGOM0hXdC9F+i7oA/5yjFElNZOXc5JnGL2ROdZ6na7RHb6gM5fPCRGcXPa8dRJuckVpycyl/6glMacl1d7+QmBj3cO5EO0MDoXRF1jQzG7dVU1qevIFS0QrUHwFfF0L8KXBFSqlcfAuM0TEPvUP6Xb7RaKCuMnk/lAuFm0NtnO++EmyrDOX3jqu3l+6XX0Xz6/fDluxsKt/7lMoOcQ/0dI1w/lQHXk8oXVFhcRabdtRgS/IKBtEK1GvogREXAE0IEQg/KKVU8+4kJ3zvU3VJNjarckHdCypD+fzjHhik68VXCPj0tRFzRgYVTz+JOSu+2QsWC35/gKsXumlpCqWDMhgMrFxbynJRErd0RfdCtFepAzG1QhFzItx7anPuPTE9Q3l2mspQfq947Ha6XnwpWC7DZLNR8fRTWHJVIM/d4Bx1cfZEG6P2iWCfLd3K5p01FBQtnL1j0RYsPBRrQxSxY8Tppm8yr5bJaKBWbc69a2bKUP7YcpWh/F7wjjro+uWL+Cf0i6nRaqXiqSewFqh10jtF0zQ6Woa5dK4Tvz/k6CqrzGXDlqqEpCu6F6K2VghRBnwMmCqycg34vpSyOxaGKeaP8NlTTaly790L53uuRmQof6B2h8pQfg/4nGN0/fIFfGOT2x/MFsqfeI8qz34X+Lx+Lp7tpLNtONhnNBpZU1/OkqWFC3KGH1VcoRDiIeAG8ATQN/nzBNA0eUyRxDSp0hrzQudoDyc7G4LttaUrWVFYl0CLFja+8Qm6XnwJr8MB6IUGy96jcuvdDfahcQ6/2RghTlk5NnY/spzaZYnLpXevRHsr/dfA30gpPxveKYT4CvAcsGm+DVPMD3aHm/5JP7TJaKCuQvn07wane4xDN44GM5SXZqkM5feC3+2m+6WX8QzrF1SDwUDZgX1kVFUl2LKFhaZp3Lw+wLWL3RF7m6rrClhbX4HZnPh0RfdCtAK1BvidGfr/EfhUtCcTQhSgZ6N4DBgAnpVS/nSW8Vb0yMEsKaX65N4FEe69shysSZBfa6HhC/h5/cZhXGEZyvct24PJqN7LuyHg9dL90ithxQYNlDz6CJm1tYk0a8HhdvloONVGf48j2Ge2mNiwpYqKReIpiVag2oGHgMZp/Q8DHXdwvu8AHqAUqAdeFkKcl1Jevs34P0d3J6o407tERe/dO8fbTqsM5fNEwOej+5XXcPX2BvtKHtpL9orlCbRq4THY7+Tsu224XaF0RXkFGWzeUUNG1uLZ0BytQH0Z+IEQ4n5gKq3wTuA/AR+N5gWEEJnAbwDrpJRO4KgQ4gXgw8CnZxhfB/wu8GfAD6K0UxHG8KiLgUn3ntlkpDaJd4wnK9f6m7jW3xRsqwzld4/m99N78A0mOjuDfUW7d5GzetUsz1JMp6N1mAun2wkEQi69ZaIEsa4M4wLY23QnRBtm/iMhxE3gk8DHAQN6FN+jUsqjUZ5rJeCXUl4P6zuPniF9Jr4N/A9g4jbHFXMQPntaUpat3Ht3yOD4MEfbTgXbKkP53aMFAvQe+hVjra3BvsId28nbsD6BVi0sNE3j+pVeGq+EZp/WNDP122soKVuc6bWijjeeFKJoxWgmsoCRaX0jwC3vrBDifYBZSvkLIcSD93DOlEZF7909mqZxrO2UylA+D2iaRv+v38bZFKqVlb9pE/lbVJBJtAT8AS6c6aCjNRSll52bzvbdtaRnLN49eLcVKCFEhZSya+rxbC8yNW4OnMB0H1MO4AjvmHQF/hXwniheU3EbhkZdDI66AOXeuxsaB5vpcfQDYDQYeWTpLpWh/C7QNI2BI8cYvSaDfbnr1lGwc3sCrVpYeD0+Th9vZbA/lOi1uCybLTuXJEVRwVgy2zeuQwhRJqXsQw+EmClBrGGyP5p36TpgFkKskFJOBVtsBKYHSKwAaoEjQggAK5ArhOgBdkopW6I4V8oTPntaUp6DZYGHm8YTj8/DiY5zwfb6slXkp6vw/Lth6MRJRi5dCrZzVgmK9uxSM9EoGXe6OXm0BafDFeyrWVrIuk2Vi269aSZmE6gVQH/Y43tCSjkmhPg58CUhxEfQo/jeC9w/begloDqsfT/wv4DNYfYoZkHTNBW9dw+c7rrIhFe/IGRYM1QS2Ltk+MxZhs+GhD5r2VKKH9yrxClKhgfHOXWsGY87lIV81fpylonilHkPbytQUsobYc1S4ISUMqImlBDCBOxAzzIRDR8HfogeOj4I/KGU8rIQYg/wqpQyS0rpA3rCzjEEBKSUPTO+ouIWhkZdDE269ywmI0vKlHsvWobG7VzuC7mj7qvejMWU3CUJkhH7hYsMnjgZbGcuqaH00UcwGJOvKF4y0t0xQsPJtmA+PaPRSP326kWzvylaonWqHwHK0YUlnLzJY1H5j6SUQ8AzM/Qf4TZ7naSUvwbUJt07INy9V1uRg8WsLgrRoGkaR9tOBXfkV+SUsTS/JsFWLTxGr15j4OixYDu9soLS/Y9hMCk381xomkZz4wBXL4QyQ1isZrbtql1QWcjni2gFamqtaTrVwOj8maO4VzRNozHMvbe8KrXuuO6FG0Ot9Dj0ezCDwcCumq0p40qZLxyNTfS99XawbSstpfw9j2M0qwCTudACGpcbumi5EarflJmVxvbddWRmL57Nt3fCrJ8aIUQjujBpwLtCiHAXnwmoAP41duYp7pTBERd2h56Sx2I2skRF70WFx+/l3Y6zwfb6UhUYcaeMtbTQ9+Yhpu5l04qKKH/yPRgtykU6Fz6fn7PvttHXHbrfLyjKZOv9tQuuRMZ8Mtdf/jX02dP/Rg9UCJ8teYFW4HBsTFPcDY3h7r3yXMwm5d6LhrNdFxn36HvCM6zpbK5QG0jvhPGODnpeeyPolrLm5VH+5BOY0lLzzv9OcE14OXm0OaK4YEV1Hhu3VWNK8e/vrAIlpXweQAjRDLw9GcCgSFI0TeNGmHtvZY1y70XD8MQIF3uvBds7qzZjVYERUePq6aHnlYNoAd3BYsnOpuLppzBnqHyFczFqn+Dk0RZcE55g3/JVetoi5V6Osh4UejDE/umdQognJ7M+KJKAAbsLu1N371ktJqpLF2f6k/lkKmPE1J1/eXYpywqWJNiqhYO7f4Cul14h4NOTlpozM6l471OYs1JvQf9O6etxcPzXN4LiZDAY2LC1mlXry5U4TRKtQP0l4Jqh3wl8Zf7MUdwLTR2hNCh15TnKvRcFN4fb6BrVc5sZDAZ2LVGBEdHiGRqm68WXCHj0C6wpPZ2Kp5/EkqPWPeei9eYgp4424/Pqs06zxcT2PXXU1BUk2LLkItrVt1rg5gz9rYAqKZoEaJoWsf6kcu/Njdfv5Z32UGDE2hJBQbp636LBOzpK14sv4Xfp961GaxoVTz6BNT8/wZYlN5qmce1iDzdkaMeOLd3K9j215OQql+h0ohWoLvQNuc3T+u8Dem8drog3/cMTjI7pd7JpFhM1yr03J2e7LzHuGQf0IoRbVWBEVPicTrp++SK+sTEAjGYLFU8+TlpxUYItS278/gDnT7XTFXYjmZufzrZdddjS1ZrnTEQrUN8B/pcQIhd9Y66GXibjL9ETuyoSTPjep7qK3JSP/pkL+8QIF3quBts7qjZhNS/erNDzhW98gs5fvojXoed4NphMlL3nALaysgRbltx43D5OHWtheHAs2FdansOmnTULvix7LIm2HtRzQogJ4HPoe58AuoHPSym/EyvjFNExPXpP5d6bHT0w4nQwMKIsu5gVhcpTPRd+t5vuF1/CO6JXzTEYDJTtf4yMqsoEW5bcjDncnDzazNhkABNA7fIi1m6swJACCV/vhTupB/U94HtCiDzAIKUcnus5ivjQF+7es5qoKpkxa5RikubhdjpHJ1M7GgzsqtmmAiPmIODx0P3iy7gHByd7DJTue5TMWhXxOBtDA2OcOtaC16Pv0DEYDKzeUE7diiL1mYuCO96iLKW0zz1KEU8a20P3CkuVe29W9MCIM8H22pIVFGaohf3Z0Px+ul89iKsvtLBf8vCDZC1fljijFgCdbXbOn2oPFr00mYzUb6+hvEplKImWqARKCGEAfh/4TfT8exErelJKVQc7QWiapirn3gHnui8zNhkYYbPY2FqxMcEWJTeaptH31q+Z6OwM9hXv2U3OKpFAq5IbTdNoutaHvBQqwGBNM7N9dx15BRkJtGzhEe2t9ueBLwJvAUuBn6AHSxQA34+NaYpo6B0axzmhb5K0Wc1Ulajovdthd41OC4yoJ00FRszK0ImTOK43BtsF27aSu17Vx7odgYDGhTMdEeKUlWNj98PLlTjdBdEK1H8Bfl9K+XXAB/xESvn76EETO2JlnGJuwvc+La3MxaQWXWdE0zSOt50moOnulpKsIlYWLk2wVcnNyKXLEQUHc1avIn/rlgRalNx4vX5OHm2mvXko2FdYksWuh5aRkaVyEt4N0a5BlQBTt56jwJTT/iAqzDxhqOi96Gmxd9Ax0q03DAZ2L1GBEbMx1txC/+GjwXZGTQ3Fex9Q79ltGB/zcOpoM47RUMKdqiX5bNhShVGtCd810b5zTeiuPYDLwIeEEFbgNwAVzZcgugfHItx7lcUqem8mfH5fRGDEmuIVFGWolDK3w9XTQ8/rbxAsm1FcTNn+faoa7m2wD41z7FdNEeK0cm0ZG7dVK3G6R6KdQX0bWIa+BvUF4CXgE+juvo/FxDLFnIQHRyyrysWo3Hszcq7nMk63vkHSZk5jW6UKjLgdHvsI3a+8huYPZSYvf+JxVdPpNvR2jXL23daw0ux6wteqJSoydD6IVqB+KKUMAEgpjwshlgBrgFYpZc/sT1XEgkBA40bHSLCtKufOzIhrlPM9V4Lt7VWbVGDEbfCNT9D90svB/HqmtDTKn3oCc4Za3J+J5qYBrjR0hZVmN7H1/loKlSdj3phToIQQFmBcCLFRSnkFQErpAE7E2jjF7ekZHGPMpbv30tOUe28m9MCIM8F9KCVZRYgiFRgxEwGvl55XXsU7qtckNZhMlD/xONY8deMzHS2gceVCN82N/cG+jEwr23fXkZVjS6Bli485HaRSSi/6GpQq8JJENEa49/KUe28GWu2dtI906Q2DgV01qpTGTGiBAL2vvxm2EVfPEqHy692Kz+fn9DutEeKUV5DBrodXKHGKAdGu4P0Z8E0hxNbJGZUigQQCGjc6Q+49Fb13K76An+Ptp4Pt1UXLKc4sTKBFyYmmafQfPsJYa2uwr3jPLrKWqtyE05kY93D8rRv0doW+e+VVudz34DLSbHeclEcRBdG+qy+gi9kJACGEP/yglFI59eNI14CT8Un3XobNQnmhmtxOp6E7FBiRZk5jW5UKjJiJ4TNnGb0S2rycv6lebcSdAfvQOKeOteCe/N4BLF1ZzOoNqvptLIlWoA7E1ArFHRERvVepovemM+p20tBzOdjeVrkRm1ltlJzO6DXJ0MlTwXb2yhUU7FT77qfT3WGn4WR7MFLPYDCwbnMlS5aqGXmsibbcxqFYG6KIDuXem5vjbaeDgRFFmQWsKlZJTacz3tZO/1u/DrbTKyspeehBNRsIY6acehariS331VKkKgbEhVnXoIQQr0wWKZxqv0cIoeoSJ5DOficTbj11f6bNQnmRcu+F02bvpM0eSmy6u2YbRoPaLBmOu3+AntdeD4ZHWwsKKDvwGAaTKpw3RcAfoOFke4Q4ZWalsevhFUqc4shcM6j9QLhv5F+AeuBmzCxSzEpTWGqj5VV56o43DF/Az7G2UGCEKFpGSZYqQx6Od9RB10svE/DpaynmzEwqnnwPpjTlAp3C7fJx+nhk9dvCkiy27FyCNU0FQ8STud7t6Vc/dTVMIP7pm3OVey+CCz1XcLidAFjNVrZX1SfYouTC73LR/fIr+CcmADBarZQ/+QTmLDUjmMIx4uLksWYmJguAAtQsLWTdpkq11psA1O3AAqKzz4FrsjJnVrqFskK1w38Kh9vJue7IwIh0i9qXMkXA56Pn1YN4hvXUmQajifLHD5BWqHISTtHX4+Dsu634vHqQsqp+m3iiEaj/JoRwTj62Ap8UQgyFD5BSfmneLVPcQrh7b5ly70XwTvsZ/AH9wlKYkc/q4uUJtih50DSNvjd/xUR3d7Cv5JGHSK+sSKBVyYOmabTcGIxIW2Q2m9i0o4bSipwEW5fazCVQh4FtYe3jwPQNJRqgBCrG+FX03m1pH+miZbgj2N69RAVGhDN4/B2cN0PLxoX37SR7hRJw0KNiLzd00npjMNhnS7eyfXctOXkqHizRzCpQUsoH4z76VacAACAASURBVGSHYg46+hy4PfoMITvDSqmqzgmAf1pgxMqipZRmFSfQouTC3nAe+/kLwXbuunXk1atNywBej48z77Yx0OsI9uUVZLBtVy1pNpUwJxmI6xqUEKIAeB54DBgAnpVS/nSGcX+OXsV3yeS470opvxFPW5ON8M25KnovxIXeq4y69AuM1WRhR9WmBFuUPDgamxg4/k6wnVlXR9Hu+9VnBxhzujl1tAWnI1TDqaI6j43bqjGpGk5JQ7z/E98BPEAp8CHge0KItTOMMwD/Gb1y7wHgj4QQ/yluViYZfn+Am50qem86TvcYZ7suBdtbVWBEkImuLvoOvRVs28rKKN33iCo6CAz2Ozl6qClCnFasKWXTjholTklG3GZQQohM9Aq866SUTuCoEOIF4MPAp8PHSinDy8hLIcQvgV3o+7BSjvY+J+7JyKKcTCsl+co3DvBO+9lgYERBRj5rSlYk2KLkwDM0TPcrB9Em3xtLbi7l7zmA0ayCdtubh7h4toNAQA+GMBqNbNxWTWWNuulLRuJ5u7AS8Espr4f1nQdmmkEFEUIYgD3opeZTkqb24eBjFb2n0zHaTfNwW7C9u2arCowAfM4xfSOuxw2AKT2diqeewGRL7ZmlpmlcvdDN+dPtQXGyppm578GlSpySmHjeUmUBI9P6RoDsOZ73BXQh/YcY2JT0+P0BbnaNBtsrVOVcPTCiNRQYsaKwjrLskgRalBwEPB66X34Fn1PfFWI0Wyh/4nEsOakdKu3z+Wk42U5PmJs8OzedbbtqychUhRiSmXgKlBOY/k3JARwzjAVACPFH6GtRe6SU7hjalrS09TrwhLn3ipV7j4u9khGXLtoWk4Ud1SowQvP76XntddyDeri0wWCg7MA+bCWpLdwT4x5OHWth1D4R7Cspz2HzjhrMFpV7MNmJp0/kOmAWQoQvFGzkNq47IcT/g7429YiUsmOmMalAePTeimrl3nN6xjjbfTHY3lq5gQxLaou2pmn0vfU24x2hr0nxg3vJqKlJoFWJxz40ztFDTRHitHRlMdvur1XitECI2wxKSjkmhPg58CUhxEfQk86+F7h/+lghxIeArwIPSSlTNjGtzx/gZlj1zuVV+Qm0Jjl4t/0cPr+e7ik/PZe1JSsTbFHiGTp5Csf10NJuwbat5KxelUCLEo+q4bQ4iHdYz8eBHwJ9wCDwh1LKy0KIPcCrUsqprJV/CRQCp4QQU8/9Zynlx+Jsb0Jp63Hg9elfsLysNIryUnuhu2u0h5tDodLku5dsT/nAiJHLlxk+czbYzlm1ivytWxJoUWK5bQ2nnUsoKp1ruVuRbMRVoKSUQ8AzM/QfQQ+imGrXxdOuZKUxfHNuirv3AoEAR8MyRiwvrKU8xQMjxlpa6H/7aLCdUVND8d49Kfs5CfgDnD/dQWdbKOo1MyuNbbvryMpW5UQWImpjRJLi8wdo6Q5376V29N6lPol9Qn8/zCZzymeMcPX20XPwDfRUmJBWXEzZ/n0pW3RQ1XBanKj/XJLS0j0acu9lp1GYm7ruvTHPOKe7QvnktlZsINOaurkIPfYRul95Fc0/uRE3O5vyJx7HaEnN/HGOERenjjUzHlbDqbqugPWbKjGqzBALGiVQScqNsNIaK1J8c+6JjlBgRF56LutKxBzPWLz4xiciig6a0tIof+oJzBmpKdgz1XBatb6cpStVDafFgBKoJMTr89MStjk3lXPvdTl6aRpsCbZ312zFmKL55AJeLz2vvIp3RHd1Gkwmyp94HGtean4+mpsGbqnhVL+jmrKK3ARbppgvlEAlIa3dDryT4bEFOTYKclLTvRcIBDjWeirYXlqwhIqcsgRalDi0QIDe19/E1dc32WOgdN+j2MpS7/3QAhqXG7pouTEQ7FM1nBYnSqCSkMYOFb0HcLn/OsNhgRE7UzRjhKZp9B8+ylhrKMS+eM//396Zh8lxVAn+l5nV1dX3pcMtS61bT7Iky9Zt+cLmGBZ2B7CHHc8Yxiwws8CaHWOWawYGGBgPxwDDssbArgcP+xlYdowxs3g5bAw+JbUkW7JlK2Tral1tdavP6qOOzJw/Iqu7qtTqQ1YdXR2/r+urzIjIzFdR1fHyRbx872qql8w8Z9dEPMme7W10mBxOMwKjoIqMRNLl2Om06b0Z6r03mBhi18lRx4j1zWupDlcVUKLC4LsuXTtb6XvxxZGyhiuvoG7tmgJKlX983+fVU3289PxpBvpHo56ZHE6ljVFQRcaRU30kg+m9phk6vXei7zRPHN1Jwk0AUF9Ry9q5M8sxwvd9oq8comv7DhL9o9ZCzfJlNG7dUkDJ8ovn+Zxs6+aQ6iDaN5xRt/yyuay4bO6MnWGYCRgFVWQcyprem0kMJ2NsP76Hg52Z0a22tWzEsWfO8z1DJ0/R+fQzxDo6Msor589nzo03zIgB2U16tB3t4rDqYGgwnlHnODaXb5zPpS0m9FepYxRUETE4nODoDJze832fw91tPNW2i+HE6F1yOBTmmpZNzK9tLqB0+SN2touu7Tsy1ppAu5I3bNxA3ZrVJf8gbiKe5Oihsxx5uZN4LJlRFypzWLikicXLZxGpMOtNMwGjoIqEI6d6eWz3Cdwgmdqs+goaZsD0XjQ+wJPHWmnrOZlRvqRxIdtaNsyISOXJ6ABdra30vaRIRYYA7UZef/nl1K+/Aqe8tEP1xIYTHD7YybHDZ0eeaUoRLg+xePksFi1toixshqyZhPm2C0w84fLk3pO8eKQro3y9lHacOd/3ebHj5YyHcAGqwpVcs3ATC+vnF1C6/ODF43Q/u5fevXvxkunWgkWNLKdpy2ZC1dXnPb4UGIzGOKQ6OH60G8/zMuoqKsMskdm0LGrECRkniJmIUVAF5GRHlEdb2+hLC9FSGSnjhg3zWVzCDxt2D/Xy+6PbORMdfY4Fy+Ky2cvZPP8Kwk5pT9/4rkvfSwfoat01EhEiReWCBTRdtZXyWaWdFqKvd4hDBzo4dbxn5EHbFNW1EZatnMO8BfXYdumvtxnOj1FQBSDpejzz/Gn2vpy5CL58QT3XXzmfSIkGt3Q9l+faX+TZUy/g+aN3y/UVdVy3aAuXVM8uoHS5x/d9Bo4c4ewzO0aiQaQob2qiadtVVC4obcuxq3OAQwfO8GraWmuK+sZKlq2cw9x5tTPCEcQwMaU5EhYxr3YN8sjONrr7R50BysMO1185nxUl7JXUHu3g8aM7RiKSA9iWzRXNq7myeXXJe+kNt7fT+fR2htvbM8pDVVU0bd1C9YrlJTso+75PR3s/rxzooKszek79rLk1LJPZNM2pLtk+MFwYRkHlCdf12PXSq+w+cAYvbUqj5ZIabtzYQnWJeiXF3QQ7TzzHix0vQ9rnnlM9i+sWbaGxorQ9FeM9vXRt30H0cKbrvB0O07D+SuouX4sdKs1/Q9/zOX2yl1cOnMlIuw46qOvcebUsWzmH+saZGejWMDGl+Z9RZJztHeKRnW10pP2TloVsrll3KZctbizZu8ZjPSd48lgrA/HBkbKQE2LzpVdw2ZzlJZ0NNzk4RPeu3fTt35+xxmLZNnVrVtOwcQNOpDS9ND3X4/ixbg6rDgaisYw6y7KYv7CBJTKbmhngpWp4bRgFlUM8z+e5lzvY8cLpEfdxgHmzqnn9pgXUVZem6/BgYoin23ZnpGcHWFA3j2sXbqa6vHRDFnmJBL37nqd7z7N4iURGXfWyZTRt3UxZbW2BpMstyYTLscNdHHm5g+GhzM/uODYtixtZIrOpqAwXSELDdMMoqBzRG43xaGsbpzpHM3w6tsXWNc2sWz67JL2TfN/n4NnDPHN8D/HkqGdipCzCtgUbWNq4sGStRd/z6FcH6drZSnJgIKOuormZpm1XEZlbmo8OxGNJjrzSydFXOknEM59hKgs7LFo6i0XLZlEeMcONYWqYX8xFxvd99h8+y1P7To1kxAWY3VDBGzcvLNnYen2xKI8f3cGpvkwngBWzlrB1wXoiodK0Fn3fZ7DtOGef2U68K/NZtnBDA01XbaFyYWkq5qHBOIcPdtJ2+Cyum/kMU3mkjCUrZrFwSROhstJ2gDHkDqOgLiLRoQS/3dVGW/tocE/bsti4ai4bVs3FKUGryfM9nn/1ALtO7sP1Ru+eq8uruG7RlpIOUxTr6KDz6e0MncyMguFUVNC4eRO1q1ZilWByxWh/jEMHznCyrRvPy3yGqbK6nKUym/kLG0yEccNrxiioi4Dv+xxs6+bx504SS5viaKyN8IZNLcwpUS+lzsEuHj+6g86BNMvBsrh87ko2zFtLWYk+cJvo66dr5076D76cUW6Hyqi/ch31V6zDLiutz+77Pr3dQxxSZ2g/2XfOw7W19RUslTnMm1+HVYI3YobCYBTUa2QoluR3e05kRCG3LIt1y2exdU0zoRK8i0x6LrtP7WNf+0sZA1VjZQPXL9rC7KrSjILgxmJ0795D774X8NOsRcuyqL3sMho2bSBUOf1vRpJJl/7eYfp6hvV77xB9vcPnxMgDaJxVxbKVc5h9SU1JTmMaCotRUK+BVIDXweFRj6XaqjCv39TCpbNLM4baqb52Hj+2k77htGlM22bjvMu5fO4q7BKc0vJdl97nX6Br1x68eKbbdNXiRTRt3UK4Yfo9ZO37PoPROH2BEtLKaJjBLNfwsZhzSS3LVs2hcVbpemQaCo9RUBdAPOHyxHMneelo5qL46iVNXH35PMIluCgcS8bZfnwPqvNQRnlzzVyuXbSZ+khpuU67w8PEu3uIdXTQu3dfRtJAgMjcuTRt20pF8/RYY0vEk/T1Do8oof7glUyeaxWdj7Kww5xLalkqs6mtL/0o84bCYxTUFDlxpp9HW4/Tn5ZErSpSxo0bF7CwubQGadB32Ue6j/NUWytD6bmanDK2LliPzFo6bad2fN8nGY2S6O4m3t1DvLtHb/f0nBPENUVZbS1NV22lasniovzcvuczEI2NKqOeIfr6hhkaiE98cIBlWVTXlFNTF6G2voLaugg1dREiFWVF+ZkNpYtRUJMkkfTY/vxp9r4ycwK8DsQHefJYK8d6TmSUL25o4eqWjVSGp8ddtO+6JHr7iPd0E+/qJtGTUkY9eMnExCcAnEhEJw1cfVnRJA2Mx7RV1NczRH+fXjOK9g2f4/I9HuHyUIYSqq2roKa2HLsE104N04/SG1VzQPvZAR5pbaOnf3RuPhIOcf36S1m+YPqtPUyE53sc6DjEjhPPknBHB/DKcAVXt2xiccOCAkp3frx4nHiPVjzx7m7i3d0kuntI9Pae43U2EZbjEG5ooKy+nsic2dSsWlmQpIG+75NMuAwNJUYcF1LrRdnRGsbDti2qayOBIqqgtl5vl0dKy9vQUFoYBTUOruvRGgR4TR/gFl5Sy40bF1A1TQO8up5LND5AND5If2yAaDxKf2yA/vgA0fgAA/HBcwb0VUGupvJQYcPU+L6POzSklVBPN/GuYFquu/ucCA6TwSkvp6yhgXBDvVZIwXaoJjdeaZ7rEY+7xGNJYrEk8ZhLPJ4kHkuSiLm6LJ4kPpwcaTdV5RqpKBuxhrQiqqCqprwko5cYShujoM7DdA7wmnSTRONa4WgFNEB/LBoooEEGE0MZkcXHoy5Sy3WLttBck78wPb7r4sXjuLEYiZ7eEUso3tNDvKv7HE+6yRCqqiLc2BAooXrC9VoZORWRC/4uU9ZNSsmMKJxYkkQ8SWx4VPnEY1rhjOWqfaE4jq2tovoINbWj60XhEpxuNsxMzC85C8/zee5gBzv2ZwZ4vXR2NTduLI4Ar3E3QTQ2QH9g+UTTFVF8gOE0Z4YLpTJcycpZS7iieQ2hKeZq8n0fP5HAjcXwYnG8WEwrnOFhvLjed1PlqbrUdiw+6XWhbCzLoqyuLssi0u/ZD876no/resSGk7hJD9f1SAbvbjJzO5FwL5p1cyGEQg7hSIiamnJq6ysC6yhCVXW5eSjWUNIYBZVGT78O8Hr6bGaA16vW6gCv+bCafN8n5sYDBaQtn2h8MGMaLj0Q6wVhWVSVVVBTXk11uIqa8ipqwlVUB++VZZWEsLTS6OtnOFAs6UrFDRSLN5ytYLTygdc+cPs+eD64voXnW7hYeD74dhl2TS12dTV2VQ1WZRVWRRWEIwx5jCqbMx7u6X6Syd4MxeO53pQcCS4mlmURLg8RDjv6vTx4D4cIR4L3codwOER5JERZ2DEhgwwzFqOgAjp7hnjgty+TSBu45jRU8obNLWMGePV9n6SXJOElSbr6PeEldZmb9e4lSbiJkbJz2o3UuyS8ZOb0m+9juZ5+eT64HiHPw3J9LG+03HI9CN5tz6fCCVNphYnYZVRYZZRbISJWGWEcwjjgunjJftxkN17SxU8mGUwmiSY9vKSLpy+ND3h+oBiwtNLIej+3TXlaG/DT6jLPE9Sd08YC28ZyHKyyMuxwOU55GLu8XL/Kgp+tD0SB6MhGbn4c4xAqc0aVTTiUqXTK05RNsB0qc4p6ethgKCbyqqBEpBG4F3gT0Al8Sin1wzHaWcCXgPcHRfcCn1BK5Ww+5aVjZ+kaPAqD3Vi+yyVNEar8CE89/SzJRBLXTZJMurhuEjfp4rlJfA+s1Ejr6XfL90fKM+r9oI0Plufj+2D5ZLR1ACdVF5xT6yorsEcsfN8Cy8axHGzLwbYdHCuEldq3bCzLxsciCvT7+ljfd/Hxgm1GzsfIe1nwukhYllYwjoOVUjbBNo6D7Tjg2Fi2g+M4WMG2FWznGsuycBwbJ2TjOBZOyCGUtu04lt4PytOtm/KU8gk7xh3bYMgh+bag7gbiwFzgCuAXIrJXKbU/q91fAG8H1qHH0t8Ah4Hv5EqwwdMvUPncMXxChByLgSMJBug/T2sbyJ03m4WFbVnYto2NHbwH+5YuI+0mfGTZ3edizKyNypGlYMhWNiNKZbQMOxi0L4KVkKFEQnagQFJKJdhP27Ydi1CZM1rnnP8427GMJWMwFDl5U1AiUgXcDKxRSkWBJ0Xk58C7gU9mNb8N+JpS6kRw7NeAPyeHCir56iChKbpQW1joPwsLPaDqsmA/vT6os1JHBm3S26YGzHMGTiuos21tgaTvW9qiwrawbAvLCpRD2rZl21n7+hjL1nW2o1+WrS0GK9i37UAp2ha2beE4+ljb1pZGqty29YBv21bacWmvVJmTeb7sstFz6n3jAGAwzGzyaUGtAFyl1MG0sr3A9WO0XR3UpbdbnUPZ2PrG64j1P85A/zB2oAhs28YJtp3UtmMTchwc2xkd4G2tlPQgr5VHasC3gzaWk7YdDMJWqo0zer0R5RBycEKh4LjUMcG57UAZ2trS0joqvd4alWOsOitVP4YyNBgMhiIhnwqqGujNKusFaibRtheoFhErV+tQTfNmcdPtN+Xi1AaDwWC4APK5whsFsqOp1sKYCz3ZbWuBaC6dJAwGg8FQXORTQR0EQiKyPK1sHZDtIEFQtm4S7QwGg8FQouRtik8pNSAiPwX+VkTej/biexuwbYzmPwDuFJGH0X5pHwW+lS9ZDQaDwVB48v0Qx4eACuAM8CPgg0qp/SJyrYikP2X5XeBfgeeBF4BfBGUGg8FgmCHk9TkopVQX+vmm7PIn0I4RqX0f+HjwMhgMBsMMxDwGbzAYDIaixCgog8FgMBQlpRQs1gFob28vtBwGg8FgmCRpY/Y5eX1KSUE1A9x6662FlsNgMBgMU6cZOJReUEoKqhW4FjhNWvxUg8FgMBQ1Dlo5tWZXWPnICGowGAwGw1QxThIGg8FgKEqMgjIYDAZDUWIUlMFgMBiKEqOgDAaDwVCUGAVlMBgMhqLEKCiDwWAwFCWl9BzUtEFEGoF7gTcBncCnlFI/HKPdx4DbgIVBu28rpb6aT1lzzWT7Iq19GNgHVCul5udHyvwwlb4QkfXAPwLrgQHgLqXUN/Mlay6Zwv9HOfBN4B1AGfAU8AGl1Mk8iptTROR24D3AWuBHSqn3jNP2I8An0BkjHkBni4jlQcycYSyownA3EAfmArcC94jI6jHaWcCfAQ3Am4HbReSWvEmZHybbFyk+hk7XUopMqi9EZBbwS3QKmiZgGfDrPMqZayb7m/hL4CrgcmAe0EPp5Y07BXwR+KfxGonIHwCfBF4PLAKWAJ/PtXC5xlhQeUZEqoCbgTVKqSjwpIj8HHg3+gc2glLqK+m7IvIQcDXw43zJm0um0hdB+8XAu4A7gf+ZT1lzzRT74k7gV0qp+4P9GPBS3oTNIVPsh8Xofng1OPbHwNfzKW+uUUr9FEBENgLjzRjcBtyrlNoftP8CcD9j/B9NJ4wFlX9WAK5S6mBa2V5gPKsBEbHQoZz251C2fDPVvvgW8FfAUK4FKwBT6YutQJeIPC0iZ0TkX0WkJS9S5p6p9MO9wNUiMk9EKtHW1v/Pg4zFyGp0P6XYC8wVkaYCyXNRMAoq/1QDvVllvUDNBMd9Dv19fT8HMhWKSfeFiLwDCCmlHsyHYAVgKr+L+eg75r8EWoAj6AzVpcBU+uEg0AacBPqAVcDf5lS64iW731LbE40rRY1RUPknCtRmldUC/ec7IFgo/TPgrdN90TOLSfVFMO3zFeDDeZKrEEzldzEEPKiUalVKDaPXGraJSF2OZcwHU+mHe4AIeh2uCvgpM9eCyu631PZ5x5XpgFFQ+ecgEBKR5Wll6zjP1J2IvJdg8VMpdSIP8uWTyfbFcvTC7xMi0o4eiJpFpF1EFuVD0Dwwld/FPiA9ynNq28qRbPlkKv2wDrhPKdUV3Lh9C9gcOJHMNPaj+yPFOuBVpdTZAslzUTDRzAtAsJjrA+8HrgAeBralFjjT2t0KfA24QSlVEovg2UymL0QkBKQPOtuA/4F2se5QSpVEepUp/C5uRLsR34AemL4CbFRKXZtfiXPDFPrh+2hL4b3AINrD878opS7Nr8S5I/jth4DPoqd2/xxIKqWSWe3eDNwH3IhOOfQAsFMpZZwkDFPmQ+hnFc6g1w4+qJTaLyLXikg0rd0X0dMXrSISDV7fKYC8uWTCvlBKJZVS7akX0AV4wX5JKKeASf0ulFK/RTuL/CJouwz40wLImysm+//x34Bh4GWgA3gL+pmoUuLT6CndT6I9WIeAT4tISzAetAAopX6JvlF5DDgWvD5bGJEvHsaCMhgMBkNRYiwog8FgMBQlRkEZDAaDoSgxCspgMBgMRYlRUAaDwWAoSoyCMhgMBkNRYhSUwWAwGIoSo6AMhhmAiNwnIvfl+BpbReSFILBxXq55Hjk+ISL35Pu6houPUVCGvCMiS0TkByJyWkSGROSAiPyDiMy9iNc4KiLvuYjne52ImIcGx+eLwFeUUhetny7we/w28M4SCoM1YzEKypBXRGQVsAuoB94JCDqkTRPwgQKKNu0QEUdEnELLASAiK9DJAx8otCxKqX50eKT3FloWw2vDJCw05Ju70aFp3pZ2p92GTkxXn2okIncCdwBzgD3A7UqpPUHde9DpRz4BfBmoQw+MH1JKxUXkd8BC4PtBvLbfK6VeJyJ/CPw1cBk6yvPPgI8ppQbSrvunQZtlQDvwjaDdY0F9Sub/BPwzcBc60nwTOvvpl5VS3x3rg4vIUXQMwRvRcfSOAn+hlHoirc0bgC+h8/scB76hlLonqFuETq3xTuBT6Eyya4EDY1zrw+gwORF03iQb8NLqP4VO2dGCTlfxLaXUf0+re7tSakta+zlBu01KqefG+Hh/BDye3pcBtoh8A60shoAvKKXuTjvvEnTa9hvQKTMeAD6ulBp6Ld8jWkH9TfAyTFOMBWXIG0GU6dcBXx9rGkgp1RO0+2N0ColPAFeiA6I+HKTdSDEHHX/uP6AzsN4MvC+ouwk4gVZwzcE+6MH6i+hIz38MXE9avDIReRM6tfZ3gTXoBHidaEVxc9CsOXj9H7SiuBW4BW0Jvg94dYJu+Cv0gHol8AjwkIjUBNcX9AB9N1pB3Ql8XkRuzjpHqm9WB58zAxG5Hh1k+DPAluBzZ8eoi6Et19WBTH8nIm8J6n4AbAysohR/Ahw4j3ICHcB3zxjl70DfCG9BK4tviMg1gZxh4FdoBbsBeBuwGR1TDi7wewzYDaya7gn7ZjrGgjLkk6XolBBqgnZ3oO/ofwQgIh8E3oxWBt8L2pQD71NKdQZt/i96oLpHKdUlIi7QGwSXBUAp9ZO0axwWkc+iB8OPB2V/DXw7ZUmgLb0ng/N3BecYOZ+ILNBFIxbQsUn0wcNKqe8Fx98B/GHwub6DVjr3KKVSSSkPi8g/oiNYp0+d3aWUemSca3wI+EnadT4M/Pv0Bkqp9NToRwKl9keBfCdF5FF0cNKUBfIu4H+Pc80W4JdjlHcBdwRBfQ+IyHXA7eh+vQU4q5T6WKqxiHwEeERE/utr+B5BW78AC4BpnXJiJmMUlKEYWQn8fWpHKZUUkV1BeYr2lHIKOI2eljsvIrISPSW3EWhA//7T/wfWoKf0JssDwEdF5CV0oryfKaUen+CYnakNpZQrIrvR1hfo6bq1QYLKFCH0FGg6z05wDUFbgunX2ZXRQOTfoS2nFehkf2Hg92lN/hmdnfZvAstuPfD2ca4ZQVtl2ezOiji/Ez09CvrzbsiKUG6hI5k3o6dMz/1wE3+PoKOcE5zLME0xU3yGfHIoeJdxW02ORNa+z8S/558DLtpi2Yi2NNIHNovMRIDjopQ6ik6m+Bl0au3/F6y3jMd4568GvorOgZR6rQHemNVucIJrjPU5RpIZishi9DTjo8Bb0dON/wSUpbV/EJgtItvQ1tNvlVInx7nmWbTjSzYTfd7HyPy869B9emac4yb6HkErLtBTtIZpirGgDHlDKdUZLHzfISI/yV6HEpE6pVQvegpwK3ogSiVt28jU0nkngBEPt2D9azlwk1LqhaDspqxjXkBPEz50nvMhIk66RRAszP8L8C8i8hvgfwEfGUeuzWky2WjL5NdB0V5ghVLqlcl8wHFQWddxgus8FhRtAPqVUp9La7M44wRKDQbTSBdPjQAAAoxJREFUpu8C/gC97jUe+8i0cFNsEBFbKZVy0NjE6BTvXrSCbFNKxc9z3gv5Hglk6QcOTyC3oYgxCsqQb1LrD78WkS8DrwDz0F5ebehppW8C3xORfcBzwEfRa04/nMJ1jgHXiMgv0N5j3cHrAyLydbTC+2DWMXcBPxORI2hlOAtYqpS6n9FptjeLyA704HcL2jLZgb6jfzs6Zfl4vFVE3g88gb7zr0v7XF8FnhKRzwM/RluEW4CwUmoqiSrvAX4lIr8PrvNhoDGt/hDQKCK3AU8B/xHtIp7t5HAf2sqKAz+d4JqPor+7bJrQjhHfBq5FOzW8Pqi7H73u9mMRuQvoQXvmXaOUSq0nXcj3CNpp43clltByxmGm+Ax5JUjbvQnt7XY/2oPr++jF9O8EbX6EHuz+AX2XvQZ4i1IqOtY5z8Pn0VbYceChYKC6FW0N7Ec/c/WZLNl+CfxntBJ9Ea0kmoK642gFdh86e+ufAL1oJbMzeDUycWbbv0d7/+0NZHmHUqovuMZu4E3oAXwPWrnchnZHnzRKqcfQ6c/vAlrRyvPBtPpn0Q4hXw2usxT9cGs2T6Jdyx+cRN8/hJ4SXJtV/iB6nGkF/g74aMqpJHhe6XWBfI+g++QL6PXEFFP+HgPeiV5HM0xjTEZdgyFPBM9BfU4pdV9hJZkcIhJBOyrcopT69STafwa4VClV0AeuReQq9E3PamNBTW+MBWUwGM4hCDv1abRl+5tJHvZ14GgqFl8BaUA/gmCU0zTHrEEZDIaxaEdbT++ebGy9wGHkSzmVanJyPFxoGQwXBzPFZzAYDIaixEzxGQwGg6EoMQrKYDAYDEWJUVAGg8FgKEqMgjIYDAZDUWIUlMFgMBiKkn8Dwx4kahA6py0AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for gamma in gamma_array:\n",
" infected_sweep = sweep_beta(beta_array, gamma)\n",
" label = 'γ = ' + str(gamma)\n",
" plot(infected_sweep, label=label)\n",
" \n",
"decorate(xlabel='Contacts per day (beta)',\n",
" ylabel='Fraction infected',\n",
" loc='upper left')\n",
"\n",
"savefig('figs/chap06-fig02.pdf')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"** Exercise:** Suppose the infectious period for the Freshman Plague is known to be 2 days on average, and suppose during one particularly bad year, 40% of the class is infected at some point. Estimate the time between contacts."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" values \n",
" \n",
" \n",
" \n",
" \n",
" 0.1 \n",
" 0.002736 \n",
" \n",
" \n",
" 0.2 \n",
" 0.007235 \n",
" \n",
" \n",
" 0.3 \n",
" 0.015929 \n",
" \n",
" \n",
" 0.4 \n",
" 0.038603 \n",
" \n",
" \n",
" 0.5 \n",
" 0.132438 \n",
" \n",
" \n",
" 0.6 \n",
" 0.346765 \n",
" \n",
" \n",
" 0.7 \n",
" 0.530585 \n",
" \n",
" \n",
" 0.8 \n",
" 0.661553 \n",
" \n",
" \n",
" 0.9 \n",
" 0.754595 \n",
" \n",
" \n",
" 1.0 \n",
" 0.821534 \n",
" \n",
" \n",
" 1.1 \n",
" 0.870219 \n",
" \n",
" \n",
"
\n",
"