{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Coronavirus Scratch Pad\n", "\n", "## The Basic Model\n", "\n", "### The Framework\n", "\n", "Assume a continuum of people, indexed by j ∈ [0, 1]. \n", "\n", "* For each j there is a measure of how gregarious/infective the people at that j are: how many others they instantaneously contact: $ g_j $. \n", "\n", "* For each j, there are measures of how many of the people at that j are currently, at each time t, susceptible, infected, and recovered: $ s_{jt} + i_{jt} + r_{jt} = 1 $.\n", "\n", "* People recover from being infected at rate $ \\rho $ (including dying): $ \\frac{dr_{jt}}{dt} = \\rho i_{jt} $\n", "\n", "* A fraction $ \\delta $ of the recovered are dead: $ d_{jt} = \\delta r_{jt} $\n", "\n", "* People get infected if they are susceptible and if one of the $ g_j $ they contact is currently infected: $\\frac{di_{jt}}{dt} = g_j s_{jt} I_t - \\rho i_{jt} $\n", "\n", "* The total number of currently infected: $ I_t = \\int_0^1 {i_{jt} dj} $ \n", "\n", "* The total number of currently recovered: $ R_t = \\int_0^1 {r_{jt} dj} $\n", "\n", "* The total number of currently susceptiable: $ S_t = \\int_0^1 {s_{jt} dj} $\n", "\n", "* The total number of currently dead is: $ D_t = \\delta R_t $\n", "\n", " \n", "\n", "### Initial Conditions\n", "\n", "* Initially, before the epidemic: $ s_{j0} = S_0 = 1 - I_0 $\n", "\n", "* Initially, before the epidemic: $ r_{j0} = R_0 = 0 $\n", "\n", "* Initially, before the epidemic: $ i_{j0} = I_0 $\n", "\n", " \n", "\n", "### Parameters\n", "\n", "There is one initial infection rate $ I_0 $, one recovery rate $ \\rho $, one share of the recovered who are dead $ \\delta $, and a continuum of gregarious rates $ g_j $.\n", "\n", "For simplicity let us cut our number of parameters to five and assign the $ g_j $ linearly based on two end parameters $ g_{min} $ and $ g_{max} $: $ g_j = j \\left( g^{max} \\right) + (1-j) \\left( g^{min} \\right) $\n", "\n", " \n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwc9X3/8ddnd3Vbt+RTkuVDBhuwLVs2JpwNBAwhkOYgQA74lTY/2tK0TdOENm1JSdom6Y+0OWgLTUgTAgFCEuIkpIQbkoCxfIJtbMu3LFmSdd/Saj+/P2bWrGXJWlmrHWn383x4H9qdaz87u37v7He+MyOqijHGmMTl87oAY4wxk8uC3hhjEpwFvTHGJDgLemOMSXAW9MYYk+As6I0xJsFZ0Ju4EJGPisivva5jqhCRL4jID+LwPBeLyD4R6RKR90/2852hjri8XjMyC/okIyKHRKTX/Y8fvs2N8XOUi4iKSCA8TFUfUdWrY/k8Ec+XIyL/LiJH3NdT4z4umoznGy8RuUJEamO0rNPW7RjuBb6lqjNU9akJPO8hEbnqbOc33rKgT07vc//jh291wycYR5B4SkRSgeeB84D1QA7wLqAZWHsWyzvtdU+XdTGK+cBOr4sw3rKgN8ApW4p3iMgR4AV3+I9E5LiItIvIKyJyXsQ8GSJyn4gcdsf/RkQygFfcSdrcLeyLROR2EflNxLzvEpFN7nybRORdEeNeEpEvishvRaRTRH59hq3zTwBlwO+r6i5VDalqo6p+UVWfdpenIrI4Yvn/IyJfcu9fISK1IvI5ETkOfHekYe6014vINhFpE5HficjyiGUeEpHPiMgO9zU9LiLpIpIF/AqYO9ovKBH5pYj82bBhO6JpanFfy/3uMjpFZKOILHLH7QcWAj93nzdNRHJF5DsiUi8ix0TkSyLij1jeH4nIbndZu0RklYg87K7j8HI+6067zl0PbSKyXUSuiFjOAhF52V3Os8CU+HWVtFTVbkl0Aw4BV40wvBxQ4PtAFpDhDv8DIBtIA/4d2BYxz/3AS8A8wI+zJZ0WsaxAxLS3A79x7xcArcDHgQBwi/u40B3/ErAfWAJkuI+/PMrreQz43hivWYHFEY//B/iSe/8KIAh8xa09Y5Rhq4BG4EL3td7mrsu0iPX6BjDXfX27gTsjnqN2WE1fAH7g3r8J2BgxbgXOL5LUM7xPgYjX0oLz6yUAPAI8Ntr7DTwFPOC+xzPdmv+vO+7DwDFgDSDAYmD+KMuZ59Z4Hc4G43vcx8Xu+NeAr7nr7zKgM/x67Rb/m23RJ6en3K2wNhEZ3m77BVXtVtVeAFV9SFU7VbUfJ5xWuFuFPpwvgT9X1WOqOqSqv3OnG8t7gX2q+rCqBlX1h8DbwPsipvmuqu5163gCWDnKsgqB+qhf+chCwD2q2h9+3SMM+yPgAVXd6L7W7wH9wLqI5XxDVetUtQX4+RlqHu5nQIWIVLiPPw48rqoDUc7/E1V9Q1WDOEE/4vOKyCzgWuAv3Pe4Efg34GZ3kj8Evqqqm9RRo6qHR3nOjwFPq+rT6vyKehaoBq4TkTKcL4u/d9ffKzjrw3jEgj45vV9V89zb8OaBo+E7IuIXkS+LyH4R6cDZqgPnZ3gRkI6z5T1ec4HhAXIYZysx7HjE/R5gxijLagbmnEUNkZpUtW+MYfOBv4r4gmwDSnFey3hrPoX75fgE8DH3C/QW4OFx1B/t884HUoD6iNfwAM6WPTivJ9r3cz7w4WHr4xKc92Iu0Kqq3RHTj/aFYeLAgt4MF3k601uBG4GrgFycZgNwftafAPqARWMsYyR1OEERqQyn2WC8ngOucdvCR9MDZEY8nj1s/Ej1Dh92FPiniC/IPFXNdH+NjCWaU8R+D/gocCXQo6qvRTHPeB3F+RVSFPEaclT1vIjxI72fMPL6eHjY+shS1S/j/MLKH/aelMXyhZjxsaA3Z5KNEwzNOEH5z+ERqhoCHgK+JiJz3a3/i0QkDWjCafpYOMpynwaWiMitIhIQkY8Ay4BfnEWND+OEzo9F5FwR8YlIoYj8rYhc506zDbjVrXE9cPlZPM9/A3eKyIXiyBKR94pIdhTzNgCFIpI72gRusIeA+xjf1nzUVLUe+DVwnzhdUn0iskhEwuvj28BnRGS1+xoXi0j4C7mBU9/PHwDvE5Fr3PWa7u7ELnGbe6qBfxSRVBG5hFOb5UycWdCbM/k+zk/uY8Au4PVh4z8DvAlswtkh+BXAp6o9wD8Bv3V/1ke2Y6OqzcD1wF/hfIl8FrheVU+Mt0C32eMqnDb+Z4EOnB2MRcBGd7I/xwmaNpyt5nH3J1fVapx2+m/h7DiuwdnBHM28bwM/BA6462O04xa+D1yAE6KT5RNAKs772Qo8idv0pao/wnnfHsXZefoUzo5lgH8B/s6t/zOqehTn197f4nyxHwX+mncy5VacHdctwD3uazMeEVW78IgxU4GIfAL4pKpe4nUtJrHYFr0xU4CIZAJ/AjzodS0m8VjQG+MxEbkGp/mjAafZxJiYsqYbY4xJcLZFb4wxCS6qkzW5XdK+jnPo97fdvrKR4z+Nc1RdEOcn6B+Ej6gTkSGcnhkAR1T1hjM9V1FRkZaXl4/nNRhjTNLbvHnzCVUtHmncmEHvnvDofpxzWdQCm0Rkg6ruiphsK1Clqj0i8sfAV4GPuON6VTXaQ8EpLy+nuro62smNMcYAIjLq0cfRNN2sBWpU9YB77o3HcPrPnqSqL7p9p8Hpa11ytsUaY4yJrWiCfh4R5z/B2aqfN8q0AHfgnJY1LF1EqkXk9dFOuyoin3SnqW5qaoqiJGOMMdGKpo1eRhg2YlcdEfkYUMWph5iXqWqdiCwEXhCRN1X1lBMnqeqDuP2Hq6qqrBuQMcbEUDRb9LU4Z7ULK8E5KdUpxLnM2OeBGyJPVavu1YtU9QDOecUrJ1CvMcaYcYom6DfhnCt7gTiXbbsZ2BA5gYhU4pzu9Ab3HNfh4fnuSa4Q5wpBF+OcY8MYY0ycjNl0o6pBEbkLeAane+VDqrpTRO4FqlV1A/CvOOfA/pGIwDvdKJcCD4hICOdL5cvDeusYY4yZZFPuyNiqqiq17pXGGDM+IrJZVatGGjedr25vpom3j3ew5XAb6Sk+blgxl4DfDsg2Jp4s6M2kemTjYT7/07dOPv72qwe5/6OrWFB0pgtCGWNiyTatzKT56dZaPv/Tt3j3uTN59bO/x/23rqKuvZc7H95M3+CQ1+UZkzQs6M2k6Ogb5N6f72JNeT7/+bFVlBZk8t7lc/j6zZXsaejk3l/YPnlj4sWC3kyKB18+QGvPIPe87zzSAv6Twy9fUswfXbqARzceYXd9h4cVGpM8LOhNzJ3o6ufbvznADSvmcv6806+HfdfvVZCdFuCbL+zzoDpjko8FvYm5p7Yeo28wxJ+9e/GI43MzU7j94nKefvM4e453xrk6Y5KPBb2JuSc317KiNI+KWdmjTnPHJQvISPHzP787GMfKjElOFvQmpnbWtfP28U4+tOpMJziFvMxUrrtgDj/fXk/vgPXAMWYyWdCbmPrJlmOk+n28b8XcMae9qaqErv4gv3qrPg6VGZO8LOhNzKgqz+w8zqUVReRlpo45/doFBZQXZvKj6to4VGdM8rKgNzGzv6mb2tZefu/cmVFNLyJ8cFUJrx1o5nh73yRXZ0zysqA3MfPSHucM1VecM+L1iUd07QWzAXh21/FJqckYY0FvYuilPU1UzJxBSX5m1PMsnpnNouIs/nenBb0xk8WC3sREd3+QNw62jGtrPmz9+bN5/UALrd0Dk1CZMcaC3sTEG4daGBgKcdmS8Qf9NefNZiikPLe7YRIqM8ZY0JuYqD7Ugt8nrJ6fP+55L5iXy+ycdF7c0zj2xMaYcbOgNzGx6WAr58/NITN1/Jc4EBEuX1LMq/tOEBwKTUJ1xiQ3C3ozYf3BIbbVtrGmvOCsl3HZkmI6+4Jsr22LYWXGGLCgNzHwZm07A8EQaxacfdBfsrgIn8DLe5piWJkxBizoTQy8cagFgKqzaJ8Py81MobIsn5f3WtAbE2sW9GbCqg+1sqg4i8IZaRNazmUVxew41k5bj3WzNCaWLOjNhKgq24+2UVl29lvzYRctKkQV3jjYEoPKjDFhFvRmQura+2juHmBFyelXkhqvFaW5pAV8bLSgNyamLOjNhOw46vSSuaAkb8LLSgv4qSzLY+PB5gkvyxjzDgt6MyE7jrWT4heWzhn9alLjceGCQnbVddDRNxiT5RljLOjNBO2obeOc2dmkBfwxWd66hYWE1DnS1hgTGxb05qyFQsqO2naWx6DZJqyyLI9Uv4+NByzojYkVC3pz1g639NDZF4zJjtiw9BQ/K0vzeP2AtdMbEysW9OasvXWsHYDz5sYu6AEuXFjAW3UddPUHY7pcY5KVBb05a7vrOwj4hIpZM2K63AsXFDIUUmunNyZGLOjNWdtd38HimTNitiM2bNX8PAI+sf70xsRIVEEvIutFZI+I1IjI3SOM/7SI7BKRHSLyvIjMjxh3m4jsc2+3xbJ4461d9R0snZMT8+VmpgZYXpLLRmunNyYmxgx6EfED9wPXAsuAW0Rk2bDJtgJVqroceBL4qjtvAXAPcCGwFrhHRCZ+rLzxXEv3AA0d/SybhKAHuHBhITtq2+kZsHZ6YyYqmi36tUCNqh5Q1QHgMeDGyAlU9UVV7XEfvg6UuPevAZ5V1RZVbQWeBdbHpnTjpd31HQCTskUPsLa8gKDbfdMYMzHRBP084GjE41p32GjuAH41nnlF5JMiUi0i1U1Ndpra6eCdoI/NEbHDVZY5ffM3H26dlOUbk0yiCXoZYZiOOKHIx4Aq4F/HM6+qPqiqVapaVVw8/otLm/jbVd/BzOy0CZ+aeDR5maksLM5i6xELemMmKpqgrwVKIx6XAHXDJxKRq4DPAzeoav945jXTz57jnZw7Sc02YavL8tlypA3VEbcrjDFRiiboNwEVIrJARFKBm4ENkROISCXwAE7IN0aMega4WkTy3Z2wV7vDzDQWCin7m7qomBnb/vPDrZqfT0v3AIeae8ae2BgzqjGDXlWDwF04Ab0beEJVd4rIvSJygzvZvwIzgB+JyDYR2eDO2wJ8EefLYhNwrzvMTGPH2nrpGwyxeJKDfrV7acIt1k5vzIQEoplIVZ8Gnh427B8i7l91hnkfAh462wLN1FPT1AUw6UG/uHgG2ekBNh9p5YOrS8aewRgzIjsy1ozb/kY36IsnN+h9PmFlaZ5t0RszQRb0ZtxqGrsozEolPyt10p9r9fx89jR00mkXIjHmrFnQm3Graexi0SRvzYetKstHFbYftQOnjDlbFvRm3PY3dbFoktvnw1aW5SECW6w/vTFnzYLejEtzVz+tPYOTviM2LCc9hSUzs+0IWWMmwILejEtNY3x63ERaNT+PrUdaCYXswCljzoYFvRmXeHWtjFRZlk9HX5D97nMbY8bHgt6MS01jFxkpfubkpMftOVeVOQdObT3aFrfnNCaRWNCbcalp7GLRzCx8vpHOVzc5FhZlkZMeYOsRC3pjzoYFvRmXA03dk36g1HA+n7CyLN/OZGnMWbKgN1Hr7g9yrK03ru3zYZWleext6KSr3644Zcx4WdCbqB1o6gbiuyM2rLIsj5DCjlprvjFmvCzoTdRqmjoBb4J+ZalzxSlrpzdm/CzoTdRqGrvw+4Sygqy4P7ddccqYs2dBb6K2v7Gb+YWZpAa8+dhUluaz1a44Zcy4WdCbqNU0dcW9x02kyrI8mrsHONrS61kNxkxHFvQmKoNDIQ6d6PakfT6sssxtpz9qzTfGjIcFvYnK4eYegiH1NOjPmZVNZqrfdsgaM04W9CYq4ZOZxes89CMJ+H0sL8m1HbLGjJMFvYlK+IRi8ToP/Wgqy/LZWddB3+CQp3UYM51Y0Juo7G/sYk5uOjPSorqe/KSpLM0jGFJ21tkVp4yJlgW9iUpNU5en7fNhK8vswCljxsuC3oxJVdkfx+vEnsnM7HRK8jMs6I0ZBwt6M6b69j66B4Y8b58Pq7QzWRozLhb0ZkwnLx84BbbowWmnr2vv43h7n9elGDMtWNCbMXlxndgzCR84tc0OnDImKhb0Zkz7m7rIzUihaEaq16UAsGxuDql+n7XTGxMlC3ozpppGp8eNSPwuH3gmaQE/583LYYu10xsTFQt6M6b9TV0sKo7/qYnPpLI0nx217QwOhbwuxZgpz4LenFFbzwAnugamTPt8WGVZHv3BEG/Xd3pdijFTngW9OaOptiM2zM5kaUz0LOjNGYXPcbO4ONvjSk41Ly+D4uw02yFrTBSiCnoRWS8ie0SkRkTuHmH8ZSKyRUSCIvKhYeOGRGSbe9sQq8JNfNQ0dpEW8DEvP8PrUk4hIlSW5tmBU8ZEYcygFxE/cD9wLbAMuEVElg2b7AhwO/DoCIvoVdWV7u2GCdZr4qymsYuFxTPw+6ZGj5tIlWX5HGruoaV7wOtSjJnSotmiXwvUqOoBVR0AHgNujJxAVQ+p6g7AukAkmJop2OMmzA6cMiY60QT9POBoxONad1i00kWkWkReF5H3jzSBiHzSnaa6qalpHIs2k6lvcIja1t4ptyM2bHlJLn6fWDu9MWOIJuhH+s2u43iOMlWtAm4F/l1EFp22MNUHVbVKVauKi4vHsWgzmfY3daE69XrchGWmBjh3drYFvTFjiCboa4HSiMclQF20T6Cqde7fA8BLQOU46jMe2t/UDUzdoAen+Wbb0TaGQuPZ9jAmuUQT9JuAChFZICKpwM1AVL1nRCRfRNLc+0XAxcCusy3WxFdNYxc+gfLCqdlGD84Rsl39wZPdQI0xpxsz6FU1CNwFPAPsBp5Q1Z0icq+I3AAgImtEpBb4MPCAiOx0Z18KVIvIduBF4MuqakE/Texv7KK0IJP0FL/XpYzq5IFT1s3SmFFFdQFQVX0aeHrYsH+IuL8Jp0ln+Hy/Ay6YYI3GIzWNXVPmHPSjWVCURW5GClsOt/GRNWVel2PMlGRHxpoRBYdCHDzRPaXb58E5cGpVWR6bbYvemFFZ0JsR1bb2MjAUmjKXDzyTqvICahq7aLUDp4wZkQW9GdFUPZnZSKrm5wOw+bBt1RszEgt6M6IatxfLoineRg+wojSPFL+w6XCL16UYMyVZ0JsR1TR2UZydRm5GiteljCk9xc8F83KpPmRb9MaMxILejGg69LiJtKa8gB21bfQNDnldijFTjgW9OY2qst+9Tux0UVVewOCQsqO23etSjJlyLOjNaZo6++nsD06roF/t7pCttnZ6Y05jQW9OE+5xMx12xIYVZKWyeOYMa6c3ZgQW9OY04R4302mLHmBNeT7Vh1oI2QnOjDmFBb05TU1jFzPSAszKSfO6lHGpml9AR1+QfY12gjNjIlnQm9PUNHaxaOYMRKbe5QPPpKrcaaffdMja6Y2JZEFvTjPdulaGlRVkUpydRrUFvTGnsKA3p2jvHaSxs5+KWdMv6EWEteUFvHGwBVVrpzcmzILenGJfQycAS6Zh0AOsW1hAXXsfR1t6vS7FmCnDgt6cYm+DsyOzYma2x5WcnYsWFQLw2oETHldizNRhQW9Osbehk8xUP/PyMrwu5awsKp5B0YxUXj9g7fTGhFnQm1Psa+xk8cwZ+HzTq8dNmIhw4cJCXj/QbO30xrgs6M0p9jZ0Tdtmm7B1Cwupb+/jcHOP16UYMyVY0JuT2noGaOrsn7Y7YsMuWui0079+oNnjSoyZGizozUnhHbFLZk3vLfpFxVkUzUizoDfGZUFvTtrX6HStnI596COJCOsWFvCatdMbA1jQmwj7GrrImsY9biJdtKiQho5+Dlk7vTEW9OYdexucHjfT7Rw3I1nnttO/tt+ab4yxoDcn7W3oomKat8+HLSzKojjb2umNAQt642rtHuBE1/TvcRMmIly0sJDf7T9h56c3Sc+C3gBOsw2QMFv0AJdWFHGia4C3j3d6XYoxnrKgNwDsbUyMrpWRLq0oBuDVfU0eV2KMtyzoDeCctTIr1c/c3HSvS4mZ2bnpnDMrm1f32QnOTHKzoDeA2+NmVnZC9LiJdGlFEW8caqF3YMjrUozxjAW9AZw+9Eum2cXAo3HpkmIGgiE2HrTeNyZ5WdAbGjv7aO4e4Nw5OV6XEnNrywtIDfis+cYktaiCXkTWi8geEakRkbtHGH+ZiGwRkaCIfGjYuNtEZJ97uy1WhZvY2VXXAcCyBAz6jFQ/Fy4osB2yJqmNGfQi4gfuB64FlgG3iMiyYZMdAW4HHh02bwFwD3AhsBa4R0TyJ162iaVd9Ykb9OC00+9t6OJ4e5/XpRjjiWi26NcCNap6QFUHgMeAGyMnUNVDqroDCA2b9xrgWVVtUdVW4FlgfQzqNjG0q66DeXkZ5GameF3KpAh3s3zFtupNkoom6OcBRyMe17rDohHVvCLySRGpFpHqpib7zxhvu+o7WDY3MbfmAc6dnU3RjDRrpzdJK5qgH6m/XbTHlEc1r6o+qKpVqlpVXFwc5aJNLPQMBDl4ojthm23AOR3CZUuKeHVfE8Gh4T86jUl80QR9LVAa8bgEqIty+ROZ18TBnuOdqJLQW/QAV547i7aeQbYcafO6FGPiLpqg3wRUiMgCEUkFbgY2RLn8Z4CrRSTf3Ql7tTvMTBGJviM27LIlRaT4hed2N3hdijFxN2bQq2oQuAsnoHcDT6jqThG5V0RuABCRNSJSC3wYeEBEdrrztgBfxPmy2ATc6w4zU8Suug6y0wOU5E//i42cSXZ6CusWFlrQm6QUiGYiVX0aeHrYsH+IuL8Jp1lmpHkfAh6aQI1mEu2q72DZnJyEO/XBSK5aOot7NuzkQFMXC4sT7yhgY0ZjR8YmsaGQ8nZ9Z8K3z4dduXQmAM/vbvS4EmPiy4I+iR1q7qZ3cCjh2+fDSvIzOXd2tjXfmKRjQZ/Ewqc+WJokQQ9O80314Vbaega8LsWYuLGgT2K76zsI+ISKBLl8YDSuXDqToZDy0h47MM8kDwv6JLarvoPFM2eQFvB7XUrcrCjJo2hGGs9a841JIhb0SWxXXWKf+mAkPp9w1dKZvLynib5BuxiJSQ4W9EmqoaOPxs5+zpub63UpcXfdBXPo6g/y8l5rvjHJwYI+SW0/6pwKYGVp8gX9RYsKyc9M4Zc76r0uxZi4sKBPUjtq2/H7hGVzki/oU/w+1p8/h+d2N1jzjUkKFvRJanttG+fMyiYjNXl2xEa6fvkcegaGePFtO3jKJD4L+iSkquyobWdFEjbbhF24oICiGan84k1rvjGJz4I+CR1u7qG9d5AVJXlel+KZgN/H+vNn88LuRnoGgl6XY8yksqBPQttrnR2xy5M46AGuXz6X3sEhXrDmG5PgLOiT0LajbaSn+FiSREfEjmRNeQHF2Wn8Yrs135jEZkGfhLYcbmVlaR4Bf3K//X6f8N4L5vDCnkbaewe9LseYSZPc/9OTUO/AEDvrOlhVlu91KVPCB1eVMBAM8fPtdoVLk7gs6JPMjto2giFl9XwLeoDz5+Vw7uxsfrS51utSjJk0FvRJZvORVgAqbYseABHhw1WlbD/axt6GTq/LMWZSWNAnmS2HW1lYnEVBVqrXpUwZ7185l4BP+FH1Ua9LMWZSWNAnEVVl8+FWVtvW/CkKZ6Rx5dKZ/HTrMQaHQl6XY0zMWdAnkf1NXbT2DFJVbkE/3IdXl3Kia8BOiWASkgV9Enn9QAsAFy4o9LiSqeeKc4opzk6znbImIVnQJ5HXDzQzKyeN+YWZXpcy5QT8Pj6wah4vvN3I8fY+r8sxJqYs6JOEqrLxYAvrFhYiIl6XMyV9dO18Qqo8svGw16UYE1MW9Eni4Ilumjr7rdnmDMoKM7ny3Fk8uvGInafeJBQL+iSx8aDTPr9uYYHHlUxt/+ficpq7B/iFXX3KJBAL+iTxu/3NzMxOY0FRltelTGnvWlRIxcwZfPe3B1FVr8sxJiYs6JNAKKT8Zl8Tl1QUWfv8GESE2y8uZ2ddB5sPt3pdjjExYUGfBHbWddDaM8ilFUVelzIt/H7lPHLSA3z3d4e8LsWYmLCgTwKv1jQBcPFiC/poZKYGuHltGf/71nGOtvR4XY4xE2ZBnwRe3XuCc2dnMzM73etSpo0/uHgBfhH+46X9XpdizIRZ0Ce4noEg1YdbuGxJsdelTCuzc9O5aU0JT24+Sl1br9flGDMhUQW9iKwXkT0iUiMid48wPk1EHnfHbxSRcnd4uYj0isg29/ZfsS3fjOXVfScYHFIut6AftzsvX4QqPPCybdWb6W3MoBcRP3A/cC2wDLhFRJYNm+wOoFVVFwP/BnwlYtx+VV3p3u6MUd0mSs/vbiA7PcDaBdZ/frxK8jP54KoSfrjpKI0ddloEM31Fs0W/FqhR1QOqOgA8Btw4bJobge+5958ErhTrx+e5UEh54e0mLl9STEqSXx/2bP3J7y1iKKQ8+MoBr0sx5qxF879/HhB5RYZad9iI06hqEGgHwsfaLxCRrSLysohcOtITiMgnRaRaRKqbmprG9QLM6LbXtnGiq5+rls7yupRpa35hFjeumMsjG4/YVr2ZtqIJ+pG2zIcfMjjaNPVAmapWAp8GHhWRnNMmVH1QVatUtaq42NqSY+X53Y34fcIV59g6nYhPXVlBMBTia8/u9boUY85KNEFfC5RGPC4B6kabRkQCQC7Qoqr9qtoMoKqbgf3AkokWbcamqjz9Vj1rywvIy7TLBk5EeVEWH19XzhPVR3n7eIfX5RgzbtEE/SagQkQWiEgqcDOwYdg0G4Db3PsfAl5QVRWRYndnLiKyEKgArLEzDvY0dHKgqZv3Lp/jdSkJ4VNXLiY7PYV/+uVur0sxZtzGDHq3zf0u4BlgN/CEqu4UkXtF5AZ3su8AhSJSg9NEE+6CeRmwQ0S24+ykvVNVW2L9IszpfrmjHp/A+vNne11KQsjLTOXP3r2YV/ed4KU9drlBM73IVDtDX1VVlVZXV3tdxrSmqlz5tZeZnZPOo3+0zutyEsZAMMR7/u1l0gI+nv7UpQSsJ5OZQkRks6pWjTTOPqkJaGddBweaurnuAmu2iaXUgI+/uXYpexu6eOi3B70ux5ioWdAnoB9vqSXV7+N6a5+PuWvOm8VVS2fxtWf3cmZhsH0AAA74SURBVLi52+tyjImKBX2CGQiG+Nm2Oq5aNtN620wCEeGL7z+PgM/H3/70Tbs4iZkWLOgTzEt7GmnpHuBDq0u8LiVhzcnN4HPXnstva5r58ZZjXpdjzJgs6BPME9W1FM1I47IKO0hqMn10bRlV8/P54i920WBHzJopzoI+gdS29vDC2w18ZE2J9QiZZD6f8JUPLWcgGOIvH9/GUMiacMzUZWmQQB7ZeASAWy+c73ElyWFR8Qz+8cbz+N3+Zv7zpRqvyzFmVBb0CaJvcIjHNx3lqqWzmJeX4XU5SePDq0u4ceVcvvbsXjYdsmMBzdRkQZ8gfrr1GC3dA9z+rnKvS0kqIsKX3n8+pQWZfOqHW2nu6ve6JGNOY0GfAIZCygMv72d5SS4XLSocewYTU9npKdx/6ypauge48web6Q8OeV2SMaewoE8Av3qrnkPNPfzx5Yuw67144/x5udx30wo2HWrl7h9b/3oztQS8LsBMzFBI+ebzNSwsyuLq8+wEZl66fvlcDjZ1c9+ze1lYlMWfXVnhdUnGABb0096G7cfY09DJN2+pxO+zrXmv3fXuxRw84YT9rJx0blpTOvZMxkwyC/ppbCDoXPVo2Zwc3msnMJsSRIR/+eAFNHcP8Lmf7MDvEz5oRykbj1kb/TT23d8e5GhLL59dfw4+25qfMtICfh74+GouXlTEXz+5nZ9ts9MkGG9Z0E9Tx9v7+Prz+7jy3Jlccc5Mr8sxw6Sn+PnvT1SxdkEBf/n4Np6oPup1SSaJWdBPU1/8xS6CIeWe953ndSlmFBmpfr5z2xretaiIzz65g68/t8964xhPWNBPQ7/cUc8v36znU+9eTFlhptflmDPISgvw0O1r+MCqefzbc3u5+8dvMjgU8rosk2RsZ+w009jRx9899SYrSnK58/JFXpdjopAa8HHfh1dQkpfBN16o4cCJLr55yypm56Z7XZpJErZFP40Eh0Lc9cOt9A2GuO+mFXaGymlERPj01efw9ZtXsrOug/d+41V+s++E12WZJGFJMY189Zk9vHGwhX/6/fNZPDPb63LMWbhx5Tw23HUxBVmpfPyhjfzz07vpG7RTJpjJZUE/TTy68QgPvnKAj60r4wOrrF/2dLZ4ZjY/u+tibl5TyoOvHODar79qZ740k8qCfhp4Zudx/v5nb3HFOcV8wXrZJITM1AD/8oHl/OCOCxkcCnHTA6/xuSd30GhXqzKTwIJ+intm53H+9JEtLC/J5Vu3rrJ2+QRzSUURz/zFZfzhJQv4ydZarvh/L/GtF/bRO2DNOSZ2ZKr1662qqtLq6mqvy5gSfvjGEf7uqbe4YF4u379jLTnpKV6XZCbRwRPdfPlXu3lmZwOFWan84aUL+fhF85mRZp3jzNhEZLOqVo04zoJ+6hkcCvEvT7/NQ789yOVLivmPj64iy/6zJ43qQy1844UaXtnbRG5GCh9bV8atF863K4eZM7Kgn0YOnujm009sY+uRNm5/Vzmff+9SUqy5JiltO9rGf7xYw3O7GwC4auksbllbxiUVRfaZMKexoJ8GegeG+ParB7j/pRpS/T7++QMXcP3yuV6XZaaA2tYeHtl4hMc3HaWle4DCrFTet2Iu110wh1VlebbfxgAW9FPaUEj5yZZa7vv1Xo539HHNebP4xxvOt6MmzWkGgiFe3tvEU1uP8ezuBgaCIfIyU7hiSTHvXjqLyyuKyc20/TjJ6kxBbw2/HmntHuDx6qM8/NphjrX1sqIkl2/cUsnaBQVel2amqNSAj/csm8V7ls2is2+QV/ed4Pndjby4p5GnttXhE1g6J4c15QXObUE+M7Ntg8HYFn1cdfYN8uKeJp556zjP7W6gPxhi3cICbn9XOVcvm23nlDdnZSikbDvaxst7m6g+1MLWI230ukfbzsvLYNncHJbNyTn5d15ehn3WEpBt0Xukuz/Ijtp2Nh5sZuOBFjYfaWUgGKJoRho3VZXy0XVlnDs7x+syzTTn9wmr5+ezen4+4PTa2lnXwaaDLew41s6uunae291AeJsuNeBjfkEm5UVZLCjKYn5hJgsKs5iTl8HsnHQyUv0evhozGaIKehFZD3wd8APfVtUvDxufBnwfWA00Ax9R1UPuuL8B7gCGgE+p6jMxq34KCA6FaOjsp66tl2Otvexr7GTP8U72NHRytKUXABFYOjuHT6ybzzXnz2ZVWb5d39VMmhS/j5WleawszTs5rGcgyJ7jneyu7+RQczcHT3RzuLmbl/c2MRA89bTJOekBZuemMysnndk56RTOSCM/M4X8zFTyMlPIz0olPzOFvMxU8jJSbGfwNDBm0IuIH7gfeA9QC2wSkQ2quitisjuAVlVdLCI3A18BPiIiy4CbgfOAucBzIrJEVeN+2J+qEgwpQyH375AypEowFGJwSOkbHKJ3YIj+4BC9AyF6B4foHRyib3CIzr4g7T0DtPUO0tYz6P4doKmzn4aOPkIRrV8Bn7CwOIsVJXl8pKqUZXNzWD2/gNwM20lmvJOZGqCyLJ/KsvxThodCSn1HH4dPdFPf3sfxjj4aOvo47t7fc7yT1p4BBodGb+JNC/iYkRYgM81PVmqArDT3luo/+TctxU+q30dqwLmluX9T/b6T49Iihvt9QsDn/I28BXyCz/3r9wl+Efx+92/EcBHbkIoUzRb9WqBGVQ8AiMhjwI1AZNDfCHzBvf8k8C1x1vSNwGOq2g8cFJEad3mvxab8d7R0D7D+3185GeShUGSwh04J47PhE8jNcLZicjNSKMhKpWJmNvPy0pmTl8HcvAzm5qYzvzCL1IBt4ZjpwecT5uVlnPFgLFWle2CI1u4B2noGae0ZoLXHud/WM0jPQJCu/iA9A0Pu3yDtvYPUt/XS3R+k292AGghO/P9htERAAJ+Ie98ZcOow5/TREp4+Ynh4Gjh9WHg+AJ/PWXbk94qcUoecNowzTLt0Tg7fvKUyRmvhHdEE/Twg8oKXtcCFo02jqkERaQcK3eGvD5t33vAnEJFPAp8EKCsri7b2U6QFfFy5dOY73/I+HwH/qd/yzl8ffh/O+IjhGal+0lP8ZKT4yUh1/qan+ElPcbZWctJTbAeWSUoiwoy0ADPSApROsFNYcChEfzDEQDDEwJDztz84dHJY+O9Q5K/vkPPreygUYigEQ6HQaRtzJ6cZcjf0VFEFRQkpJ++7/wiFFOWd4arOF1rksPB8nByPs1xOXV4ookNL5PdYePCpw0aeNvygNH9yjn6OJuhHSrfh38ujTRPNvKjqg8CD4PS6iaKm02SlOWcDNMZMXQG/j4DfR1aa15Ukl2jaGGqB0ojHJUDdaNOISADIBVqinNcYY8wkiiboNwEVIrJARFJxdq5uGDbNBuA29/6HgBfU+Y2yAbhZRNJEZAFQAbwRm9KNMcZEY8ymG7fN/S7gGZzulQ+p6k4RuReoVtUNwHeAh92drS04Xwa40z2Bs+M2CPypFz1ujDEmmdmRscYYkwDOdGSs9QM0xpgEZ0FvjDEJzoLeGGMSnAW9McYkuCm3M1ZEmoDDE1hEEXAiRuXEktU1PlbX+Fhd45OIdc1X1eKRRky5oJ8oEakebc+zl6yu8bG6xsfqGp9kq8uabowxJsFZ0BtjTIJLxKB/0OsCRmF1jY/VNT5W1/gkVV0J10ZvjDHmVIm4RW+MMSaCBb0xxiS4hAl6EVkvIntEpEZE7vawjlIReVFEdovIThH5c3f4F0TkmIhsc2/XeVDbIRF5033+andYgYg8KyL73L/5Yy0nxjWdE7FOtolIh4j8hVfrS0QeEpFGEXkrYtiI60gc33A/cztEZFUca/pXEXnbfd6fikieO7xcRHoj1tt/TUZNY9Q26nsnIn/jrq89InJNnOt6PKKmQyKyzR0el3V2hmyY/M+Xqk77G87pk/cDC4FUYDuwzKNa5gCr3PvZwF5gGc41dT/j8Xo6BBQNG/ZV4G73/t3AVzx+H48D871aX8BlwCrgrbHWEXAd8CucK6mtAzbGsaargYB7/ysRNZVHTufR+hrxvXP/H2wH0oAF7v9Zf7zqGjb+PuAf4rnOzpANk/75SpQt+pMXMFfVASB8AfO4U9V6Vd3i3u8EdjPCdXKnkBuB77n3vwe838NargT2q+pEjoyeEFV9BeeaCpFGW0c3At9Xx+tAnojMiUdNqvprVQ26D1/HuXpb3I2yvkZzI/CYqvar6kGgBuf/blzrEhEBbgJ+OBnPfYaaRsuGSf98JUrQj3QBc8/DVUTKgUpgozvoLvcn2EPxbiJxKfBrEdkszgXZAWapaj04H0Rgpgd1hd3Mqf/5vF5fYaOto6nyufsDnC2/sAUislVEXhaRSz2oB0Z+76bK+roUaFDVfRHD4rrOhmXDpH++EiXoo7oIeTyJyAzgx8BfqGoH8J/AImAlUI/z0zHeLlbVVcC1wJ+KyGUe1DAicS5TeQPwI3fQVFhfY/H8cycin8e5etsj7qB6oExVK4FPA4+KSE48a2L0987z9eW6hVM3KOK6zkbIhlEnHWHYWa2vRAn6KXURchFJwXkjH1HVnwCoaoOqDqlqCPhvJukn65moap37txH4qVtDQ/jnoPu3Md51ua4Ftqhqg1uj5+srwmjryNPPnYjcBlwPfFTdRl23WaTZvb8Zpx18Sbxqcp93tPfO8/+nIhIAPgA8Hh4Wz3U2UjYQh89XogR9NBcwjwu3/e87wG5V/VrE8Mi2td8H3ho+7yTXlSUi2eH7ODvz3uLUC7vfBvwsnnVFOGUry+v1Ncxo62gD8Am3d8Q6oD38E3yyich64HPADaraEzG8WET87v2FQAVwIB41RdQw2nu3AbhZRNJEZIFb2xvxrA24CnhbVWvDA+K1zkbLBuLx+ZrsPc3xuuHsod6L8238eQ/ruATn59UOYJt7uw54GHjTHb4BmBPnuhbi9HjYDuwMryOgEHge2Of+LfBgnWUCzUBuxDBP1hfOl009MIizRXXHaOsI56f1/e5n7k2gKo411eC034Y/Y//lTvtB9/3dDmwB3ufB+hr1vQM+766vPcC18azLHf4/wJ3Dpo3LOjtDNkz658tOgWCMMQkuUZpujDHGjMKC3hhjEpwFvTHGJDgLemOMSXAW9MYYk+As6I0xJsFZ0BtjTIL7/0f8v70mZjuEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import scipy as sp\n", "import pandas as pd\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "# paramaters\n", "\n", "time = 200 # time periods we run the simulation for\n", "bins = 101 # number of bins into which the population is divided\n", "\n", "R_zero_min = 2.5 # R_0 for least gregarious/infective\n", "R_zero_max = 2.5 # R_0 for most gregarious/infective\n", "\n", "rho = 0.1 # recovery rate\n", "g_min = R_zero_min * rho # minimum gregariousness/infectiveness\n", "g_max = R_zero_max * rho # maximum gregariousness/infectiveness\n", "delta = 0.01 # death rate\n", "\n", "I_zero = 0.0001 # initial infection rate\n", "\n", "# averages over the population at the current moment of S, I, R—\n", "# susceptible, infected, recovered (or dead), initialized to their\n", "# time-zero values\n", "I_tavg = I_zero \n", "R_tavg = 0\n", "S_tavg = 1 - I_zero\n", "\n", "# time series for the entire population of the fractions S, I, R;\n", "# initialized with their time-zero values\n", "T_pop = [0]\n", "R_pop = [R_tavg]\n", "I_pop = [I_tavg]\n", "S_pop = [S_tavg]\n", "\n", "# heterogeneity across the population in gregariousness/infectiveness\n", "G = np.linspace(g_min, g_max, 101)\n", "\n", "# initial conditions across the population, for all j:\n", "R_0 = 0*G\n", "I_0 = 0*G + I_zero\n", "S_0 = 1 - I_zero\n", "\n", "# arrays to hold current-time population heterogeneity—all j at the current\n", "# moment of time—initialized to their values at time zero:\n", "R = R_0\n", "I = I_0\n", "S = S_0\n", "\n", "# arrays to hold all the numbers—for all j, and for all t—for when we want\n", "# to look at them later:\n", "R_array = [R]\n", "I_array = [I]\n", "S_array = [S]\n", "\n", "# loop to calculate all the numbers:\n", "for t in range(0, time): \n", " R = R + rho * I\n", " I = (1-rho) * I + I_tavg * G * S\n", " S = - (R + I - 1) \n", " \n", " # subloop to calculate average S, I, R across all j at the current\n", " # moment of time:\n", " I_tavg = 0\n", " R_tavg = 0\n", " S_tavg = 0\n", " for j in range(0, bins):\n", " I_tavg = I_tavg + I[j]/bins \n", " R_tavg = R_tavg + R[j]/bins\n", " S_tavg = S_tavg + S[j]/bins\n", "\n", " # update the full arrays with the current numbers for all j:\n", " R_array = R_array + [R]\n", " I_array = I_array + [I]\n", " S_array = S_array + [S]\n", " \n", " # update the time series of population averages:\n", " T_pop = T_pop + [t+1]\n", " R_pop = R_pop + [R_tavg]\n", " S_pop = S_pop + [S_tavg]\n", " I_pop = I_pop + [I_tavg]\n", "\n", "pd.DataFrame(I_pop).plot(title = \"Fraction Currently Infected\", legend = False)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1wUZ/7A8c+zS5MiimDFggULCBgVaxTsLaZZUkwuZ6KJiTG/u0vPmUQvPaaZmDPmzpjiaXo0RqMx9t57L4gIVpRelt3n98cgAiIsCCzl+/a1r5nZeWbmu4Bfhmdmvo/SWiOEEKLyMzk6ACGEEKVDEroQQlQRktCFEKKKkIQuhBBVhCR0IYSoIiShCyFEFSEJXYgyoJTSSqmWdrSLUErFlEdMouqThC4qPKVUT6XUBqVUglIqXim1XinVOXvdQ0qpdbnaRiml0pRSyUqps0qpOUopz0L2vSo7+Ybme/+X7PcjyuyDCVHKJKGLCk0pVRNYBHwM+ACNgClARiGb3aa19gTCgA7AC0Uc5gjwYK5j1gG6AhdKHrkQ5U8SuqjoAgG01vO01latdZrWepnWek9RG2qtzwJLMRJ7YeYCo5VS5uzle4GfgcyrDZRSrkqpD5VSsdmvD5VSrrnWP6OUisteNzb3zrO3naaUilZKnVNKzVRK1bDr0wtRDJLQRUV3BLAqpb5USg1WStW2d0OllD8wGDhWRNNY4AAwIHv5QeCrfG1ewjhrDwNCgXDgn9nHGQQ8DfQHWgH98m37NsYvpjCgJcZfGS/b+zmEsJckdFGhaa0TgZ6ABj4HLiilFiql6hWy2S9KqSTgNHAeeMWOQ30FPKiUag3U0lpvzLf+fmCq1vq81voCRrfPA9nrRgFfaK33aa1TgFevbqSUUsA44G9a63itdRLwBnCPHTEJUSyS0EWFp7U+qLV+SGvtDwQDDYEPC9nkDq21FxABtAF87TjMT0Af4Eng6wLWNwRO5Vo+lf3e1XWn8627yg9wB7Yrpa4opa4Av2e/L0SpkoQuKhWt9SFgDkZiL6rt6uy20+xomwosASZQcEKPBZrmWm6S/R5AHNA437qrLgJpQJDWulb2yzv7oq0QpUoSuqjQlFJtlFL/yO4PRynVGOOi5SY7d/Eh0F8pVdSFUYAXgd5a66gC1s0D/qmU8lNK+WL0gX+Tve474CGlVDullDu5uni01jaMrqIPlFJ1sz9DI6XUQDvjF8JuktBFRZcEdAE2K6VSMBL5PuAf9myc3d/9FTDZjraxWut1N1j9GrAN2APsBXZkv4fWegnGL44VGBdgV+Tb9rns9zcppRKB5UBre+IXojiUDHAhhBBVg5yhCyFEFSEJXQghqghJ6EIIUUVIQhdCiCrCyVEH9vX11c2aNXPU4YUQolLavn37Ra11gQ+mOSyhN2vWjG3btjnq8EIIUSkppU7daJ10uQghRBUhCV0IIaoISehCCFFFSEIXQogqQhK6EEJUEUUmdKXUbKXUeaXUvhusV0qp6UqpY0qpPUqpW0o/TCGEEEWx5wx9DjCokPWDMYbdagWMB/5982EJIYQoriLvQ9dar1FKNSukye3AV9oo27hJKVVLKdVAax1XSjHmdWojHF8BJjMoEyiVPTWBMueaz36ZTNe/l6edApMTOLmBk4sxNbvmms+eOruBi5exPyGEqIBK48GiRuQdfism+73rErpSajzGWTxNmjTJv9o+MVtgzTsl27Y0uHiBW01wrWlM3bzBoy545np5Nwaf5uBex/iFIYQQ5aA0EnpBGavAIuta61nALIBOnTqVrBB7j6eMl9agbddeNmve5Ru98rTL3ofNAlkZxsuacW0+ZzkdLGmQkQTpiZCR/UpPhORzcHYfpJwHW1beWF1rgk8A+LSA+u2hQSg07ADuPiX66EIIUZjSSOgx5B1P0Z9rYy2WHaWMrhPMZX4ou9hskHbZSPBXouHySYg/YbzObIP9P11rW7cdBPQyXs0jwcXdcXELIaqM0kjoC4GJSqn5GEOFJZRZ/3lFZjKBRx3jVa/d9evTLkPcHiO5R62D7V/C5png7A6BgyDoTmPq5FL+sQshqoQih6BTSs0DIgBf4BzGALjOAFrrmUopBXyCcSdMKvBXrXWRVbc6deqkq3VxrqwMiN4E+3+Ggwsh9RJ4NYDOj0CnsdItI4QokFJqu9a6U4HrHDWmaLVP6LlZs+DYcuOM/cRKcKoBXSfArX8HVy9HRyeEqEAKS+hyD15FYHaC1oPgwV9gwkZoexusex8+7gg7vjb654UQogiS0Cuaeu3g7s/hkRVQqyksnAhf3mZcaBVCiEJIQq+o/DvCw8vg9hkQtxv+3RMOL3F0VEKICkwSekWmFHQYAxPWGfezz7sH1kwz7p8XQoh8JKFXBrWbwdil0H4krPgX/PYP6VcXQlzHYWOKimJydoM7Z0HNhrD+I7Bmwm3TpbaMECKHJPTKxGSCflOM4mFr3jGKhg15V+rFCCEASeiVj1IQ+SJYUmHjJ+DdCHr+zdFRCSEqAEnolZFS0P9fkBQHy18FvzbQerCjoxJCOJh0wFZWJpNxS2ODUPjpUbh03NERCSEcTBJ6ZeZcA0Z/YyT3H8ZCVqajIxJCOJAk9MquVhPjbpe4XbDqTUdHI4RwIEnoVUG74cYDSOs+gDPbHR2NEMJB5KJoVTHwDTi6HBY+BeNXgtnZ0RGJSkRrjcWqyciykpllIyP7lZllI8tmw2aDLJsNq01fe2lNlk1jtRrzVpuxbMueWm02rDawao3W2hggTGs0xsPOtuwnnrUGjcamr83ntNXkaa+zN9BkL+daf3Xf+T9X3uUCPvt12+RfX/Q+Cvp6FnaMYSENCQ8o/RLZktCrCjdv45707x4wyvB2f9LREYkyZLNpEtMtxKdkcjk1k8spFhLSLKRkZpGckUVKRhYpGdac+avT1EwrmVYbGZarCduak7wrE5MCpRSK7MHLcs/nGxUz/2MaBT21ofI1uq5NEfvIv31Rx23fyFsSuihC29ug1QBY/Q6E3AOefo6OSBSTzaY5l5ROXEI6ZxOuTtOIS0jnXGJ6dgK3cCU1E1shZ4pmk8LDxYynqxMe2S9PVyfqeLri6mTC1cmMq7MJF7MJV2cTrmYTrs7ma8tOJlycTDibTTiZTDiZFObsl5NJYcqemvO8b8JsAnN2+6ttlAJTTsI1piYj82Yn4Oz12clYKfLO514vD9EVShJ6VaIUDHgd/t0NVr0Bwz5wdETiBhLTLRw5m8SJiymcvJjCyQspRF0y5vOfLbs6mWjg7Ua9mm60ru9FbXcXfDxcqO3uQm0P55xl7xrOOYnb1ckkya8akoRe1fgFQqeHYevn0GWCsSwcKjkji53Rl9kTk8D+2AT2xyZy6lJqznpns6KJjzsBvp7c2sqXZr4eNPSuQX1vNxp4u+Fdw1mSs7CLJPSqqPezsPMb4zbGkV84OppqJzHdwqbjl9gaFc+Wk/Hsi03Emt0/0sTHnaCGNRnZ0Z92DWvSws+TRrVq4GSWG87EzZOEXhV5+ELXx2Dte9DraagX5OiIqrzT8an8efAcyw+eZ/PJS1isGhcnE2GNa/F4RAs6N/MhtHEtvGvI3Uei7EhCr6q6TYQtn8Pqt2HUV46Opko6l5jOgl1n+HlnLAfjEgFo4efB2B4BRLapS1jjWrg5mx0cpahOJKFXVe4+0PkR42GjS8ehTgtHR1QlWG2aPw+e45vN0aw7egGbhrDGtXhpSFv6tq1Lcz9PR4coqjFJ6FVZ1wmwcYYxIMbw6Y6OplJLSrcwb0s0X208RczlNBp4u/FEZEvu7NBIkrioMCShV2WedaHD/cYF0siXwKueoyOqdBLSLHy5IYr/rjtJQpqFLgE+vDSkLf3b1ZMLmaLCkYRe1XV9ArbNhu1zIOI5R0dTaWRm2fhyQxQfrzhKYnoW/drWY1LfloT413J0aELckCT0qs63JbTsZyT1nn8DJxdHR1Shaa3548A53lh8kKhLqUS09uPpAa0JbuTt6NCEKJL8zVgdhI+H5LNwcKGjI6nQziemM/7r7Yz/ejtOZhNz/tqZOX8Nl2QuKg05Q68OWvaH2gHGWXr7EY6OpsLRWvPzzjNM+fUA6RYrLw5pw9geAdJHLiodSejVgckEtzwIf06Bi8eMbhgBGI/lv/DTXn7dHUunprV5Z0SI3LUiKi05Bakuwu4DZYadXzs6kgrj8Nkkhn+yjt/2xPLMwNZ8+2g3SeaiUpOEXl141YfAgbDrf2C1ODoah1uyN447ZqwnMS2LuY905YnIlphNUgBLVG52JXSl1CCl1GGl1DGl1PMFrG+ilFqplNqplNqjlBpS+qGKm9bhAUg5D8f+dHQkDvWftSd4/H87aNvAi8WTetKtRR1HhyREqSgyoSulzMAMYDDQDrhXKdUuX7N/At9prTsA9wCflnagohS07AdutWDfD46OxCGsNs2UX/fz2m8HGRRUn/+N60rdmm6ODkuIUmPPGXo4cExrfUJrnQnMB27P10YDNbPnvYHY0gtRlBonF2h3OxxaDJmpRbevQqw2zTPf7+aL9VE83DOAGffdIoWzRJVjT0JvBJzOtRyT/V5urwJjlFIxwGKgwAEtlVLjlVLblFLbLly4UIJwxU1rPwIsKXBkiaMjKTdWm+bZH/bw084z/KN/IJOHtcMk/eWiCrInoRf0k59/NMN7gTlaa39gCPC1Uuq6fWutZ2mtO2mtO/n5yXiXDtG0B3jWh70/OjqScmGzaZ7/cQ8/7ojhb/0CebJvK0eHJESZsSehxwCNcy37c32XysPAdwBa642AG+BbGgGKUmYyQ/BdcOwPSLvi6GjK3NRFB/h+ewxP9W3FU/0kmYuqzZ6EvhVopZQKUEq5YFz0zP8MeTTQF0Ap1RYjoUufSkUVPAKsmXDwV0dHUqZmrzvJnA1RjO0RwP9JMhfVQJEJXWudBUwElgIHMe5m2a+UmqqUGp7d7B/AOKXUbmAe8JDWOn+3jKgoGt1ilAKowne7LNt/ln/9doCBQfV4aWhbGWRZVAt2PfqvtV6McbEz93sv55o/APQo3dBEmVEKgu+Gde9D0rkqVyd9b0wCk+bvJMS/Fh+O7iAPDIlqQ54Ura6C7wJtq3J3u1xJzeSxb7ZTx8OV/zzYiRoucmuiqD4koVdXddtBrabGPelVhM2m+cd3uzmflM6M+2/Bz8vV0SEJUa4koVdXSkHrIXBiFWSmODqaUvH52hP8eeg8Lw1pS1hjGVlIVD+S0Kuz1oPBmgHHVzo6kpu2LSqed5YeZkj7+vylezNHhyOEQ0hCr86adgc3bzhcubtdUjOz+Pt3u2lUqwZv3R0id7SIaksSenVmdoZWA+DI72CzOjqaEnvn98NEx6fy7ogQaro5OzocIRxGEnp113oIpF6C01scHUmJbDpxiTkbonioezO6NJcyuKJ6k4Re3bXsBybnStntkpqZxbM/7KFpHXeeHdTa0eEI4XCS0Ks7t5oQcCscrnz3o7+37AjR8am8c3cI7i4yPK4QktCF0e1y6ShcPOroSOx25FwSczZEcW94E+lqESKbJHRhjDUKcPQPx8ZhJ601ry7cj6erE88OlK4WIa6ShC6gVhPwDYRjyx0diV0W7z3LhuOXeHpAILU9XBwdjhAVhiR0YWjRF06tB0uaoyMpVGpmFq//doC2DWpyX5emjg5HiApFErowtOwHWelwaoOjIynUzNUniE1IZ8rwIKmiKEQ+ktCFoWl3MLvCsT8dHckNXUjK4D9rTzC0fQPCA3wcHY4QFY4kdGFwcYdmPeB4xU3on646RkaWjX8MCHR0KEJUSJLQxTUt+sKFQ5AQ4+hIrhNzOZW5m6IZcYs/zf08HR2OEBWSJHRxTcu+xrQCdrtM/9O4R14GehbixiShi2v82kDNRhWu2+XY+WR+2B7DmK5NaVirhqPDEaLCkoQurlEKWvSB46vAmuXoaHJ89OdR3JzNPB7ZwtGhCFGhSUIXebXsBxkJcGa7oyMBIOpiCr/tieWBbk3x9ZQh5YQojCR0kVfz3oAyhqarAD5bcwIns4mHewQ4OhQhKjxJ6CKvGrWhQSicXO3oSDiXmM6P22MY0dGfujXdHB2OEBWeJHRxvea9jQEvHDx49Ox1J8my2Xi0V3OHxiFEZSEJXVwvoBfYLBC9yWEhJKRa+GbTKYaFNKRpHQ+HxSFEZSIJXVyvSTdjFCMHdrt8vSmKlEwrj/WWO1uEsJckdHE9Fw9oHA4nHJPQM7NsfLnxFL0D/WjXsKZDYhCiMpKELgoW0BvidkNqfLkfesm+OC4kZfDXHs3K/dhCVGaS0EXBmvcGNEStK/dDz9kQRYCvB71a+ZX7sYWozGRkXVGwRh3B2cPoR283vNwOuyfmCjujr/DKbe0wlUK9c4vFQkxMDOnp6aUQnRDlx83NDX9/f5ydne3exq6ErpQaBHwEmIH/aK3fKqDNKOBVQAO7tdb32R2FqHjMzkaN9HLuR5+zIQoPFzMjOvqXyv5iYmLw8vKiWbNmKCUDYojKQWvNpUuXiImJISDA/ofqiuxyUUqZgRnAYKAdcK9Sql2+Nq2AF4AeWusg4P+KE7yooJr3hktHITG2XA53MTmDRbvjuLujP15u9p+VFCY9PZ06depIMheVilKKOnXqFPsvS3v60MOBY1rrE1rrTGA+cHu+NuOAGVrrywBa6/PFikJUTAG9jenJNeVyuHmbo8m02niwW7NS3a8kc1EZleTn1p6E3gg4nWs5Jvu93AKBQKXUeqXUpuwumoICHK+U2qaU2nbhwoViByvKWb1gcK9TLt0uVptm3pZoerb0pWVdGcBCiJKwJ6EX9GtC51t2AloBEcC9wH+UUrWu20jrWVrrTlrrTn5+cgdDhWcyQbNbjQujOv+3vHStPXqB2IR07g1vUqbHcQSz2UxYWBjBwcHcdtttXLlyxdEhlUhERATbtm0rcN2IESM4ceIEAM2aNePixYvF3n9UVBTBwcHF2mbOnDlMnDgRgJkzZ/LVV18V+7gFiYiIoHXr1oSEhNCmTRsmTpxYat+3VatWMWzYMAAWLVrEK6+8Uir7BfsSegzQONeyP5C/UzUGWKC1tmitTwKHMRK8qOwCekHiGbh0vEwP8+3W0/h4uNCvXd0yPY4j1KhRg127drFv3z58fHyYMWOGo0PKkZV183Xv9+/fj9VqpXlzx9bceeyxx3jwwQdLbX9z585lz5497NmzB1dXV26/PX9P880bOnQoCxcuJDU1tVT2Z09C3wq0UkoFKKVcgHuAhfna/AJEAiilfDG6YE6USoTCsZpHGNOTq8rsEBeTM/jjwDnu6tAIVydzmR2nIujWrRtnzpzJWX733Xfp3LkzISEhec7UvvrqK0JCQggNDeWBBx4A4NSpU/Tt25eQkBD69u1LdHQ0CQkJNGvWDJvNBkBqaiqNGzfGYrFw/PhxBg0aRMeOHbn11ls5dOgQAA899BB///vfiYyM5LnnniMlJYWxY8fSuXNnOnTowIIFCwBIS0vjnnvuISQkhNGjR5OWllbgZ5o7d26ByS4qKoq2bdsybtw4goKCGDBgQM4+jh07Rr9+/QgNDeWWW27h+PG8Jwy5z7wBhg0bxqpVqwD44osvCAwMpHfv3qxfvz6nzauvvsq0adMA4wz7ueeeIzw8nMDAQNauXZvz9Rk1alTOZ+rSpcsN/+q4ysXFhXfeeYfo6Gh2794NwDfffEN4eDhhYWE8+uijWK1WACZMmECnTp0ICgrK8/38/fffadOmDT179uSnn37KeV8pRUREBIsWLSo0BnsVedui1jpLKTURWIpx2+JsrfV+pdRUYJvWemH2ugFKqQOAFXhGa32pVCIUjuXTHGr6GxdGOz9SJof4eccZsmya0Z0bF934Jkz5dT8HYhNLdZ/tGtbklduC7GprtVr5888/efjhhwFYtmwZR48eZcuWLWitGT58OGvWrKFOnTq8/vrrrF+/Hl9fX+Ljjad1J06cyIMPPshf/vIXZs+ezaRJk/jll18IDQ1l9erVREZG8uuvvzJw4ECcnZ0ZP348M2fOpFWrVmzevJnHH3+cFStWAHDkyBGWL1+O2WzmxRdfpE+fPsyePZsrV64QHh5Ov379+Oyzz3B3d885S73lllsK/Fzr16/n3nvvLXDd0aNHmTdvHp9//jmjRo3ixx9/ZMyYMdx///08//zz3HnnnaSnp2Oz2Th/vuh7KeLi4njllVfYvn073t7eREZG0qFDhwLbZmVlsWXLFhYvXsyUKVNYvnw5n376KbVr12bPnj3s27ePsLCwIo8JRrdZaGgohw4dwsXFhW+//Zb169fj7OzM448/zty5c3nwwQd5/fXX8fHxwWq10rdvX/bs2UNgYCDjxo1jxYoVtGzZktGjR+fZd6dOnVi7di2jRo2yK5bC2HUfutZ6MbA433sv55rXwN+zX6IqUcq4ffHwErDZjH71UqS1Zv7WaDo2rU2rel6luu+KIi0tjbCwMKKioujYsSP9+/cHjIS+bNmynISUnJzM0aNH2b17NyNGjMDX1xcAHx8fADZu3JhzdvfAAw/w7LPPAjB69Gi+/fZbIiMjmT9/Po8//jjJycls2LCBkSNH5sSRkZGRMz9y5EjMZnNOHAsXLsw5u01PTyc6Opo1a9YwadIkAEJCQggJCSnw88XFxXGja2IBAQE5SbNjx45ERUWRlJTEmTNnuPPOOwHjARp7bd68mYiIiJzjjR49miNHjhTY9q677spzXIB169bx1FNPARAcHHzDz1QQnX0d6c8//2T79u107twZML6/desaXYXfffcds2bNIisri7i4OA4cOIDNZiMgIIBWrYxe6DFjxjBr1qyc/datW5fY2NK5NVieFBVFC+gNu+bCub3G4BelaPupyxy/kMI7I8q+qqK9Z9Kl7WofekJCAsOGDWPGjBlMmjQJrTUvvPACjz76aJ7206dPt+uWtatthg8fzgsvvEB8fDzbt2+nT58+pKSkUKtWLXbt2lXgth4e10oSa6358ccfad269Q2PUdTnu9H90q6u14YNNJvNpKWl5STGwjg5OeV0IwF59m/v7XxXj202m3OuFdhz7IJYrVb27t1L27ZtOX/+PH/5y194880387Q5efIk06ZNY+vWrdSuXZuHHnooJ+7CYk5PT6dGjdIZ/FxquYiiBfQypmVw++L8rafxdHViaPsGpb7visbb25vp06czbdo0LBYLAwcOZPbs2SQnJwNw5swZzp8/T9++ffnuu++4dMnotbza5dK9e3fmz58PGP3WPXv2BMDT05Pw8HCeeuophg0bhtlspmbNmgQEBPD9998DRiK72v+b38CBA/n4449zkt3OnTsB6NWrF3PnzgVg37597Nmzp8Dt27Zty7Fjx+z+OtSsWRN/f39++eUXwPjLIf9FwWbNmrFr1y5sNhunT59my5YtAHTp0oVVq1Zx6dIlLBZLzuezV8+ePfnuu+8AOHDgAHv37i1yG4vFwgsvvEDjxo1zrl/88MMPOV1E8fHxnDp1isTERDw8PPD29ubcuXMsWbIEgDZt2nDy5Mmc6wTz5s3Ls/8jR44U++6eG5GELopWswH4BpZ6ffTUzCwW741jWEgDPFyrxx+LHTp0IDQ0lPnz5zNgwADuu+8+unXrRvv27RkxYgRJSUkEBQXx0ksv0bt3b0JDQ/n7342ezOnTp/PFF18QEhLC119/zUcffZSz39GjR/PNN9/k6Z+dO3cu//3vfwkNDSUoKCjnYmd+kydPxmKxEBISQnBwMJMnTwaMC3zJycmEhITwzjvvEB4eXuD2Q4cOzblgaa+vv/6a6dOnExISQvfu3Tl79mye9T169CAgIID27dvz9NNP5/TfN2jQgFdffZVu3brRr1+/G/br38jjjz/OhQsXCAkJ4e233yYkJARvb+8C295///05X5OUlJScr1+7du147bXXGDBgACEhIfTv35+4uDhCQ0Pp0KEDQUFBjB07lh49egBGl9KsWbMYOnQoPXv2pGnTpnmOs3LlSoYOHVqsz3FDWmuHvDp27KhFJbLoH1q/Vl9rS0ap7fLnHTG66XOL9KbjF0ttn/kdOHCgzPYtDKmpqbpLly46KyvL0aEUKSsrS6elpWmttT527Jhu2rSpzsgovZ/p4jp79qzu06fPDdcX9POLcTNKgXm1epwWiZvXvDds/RzObIem3Upllz/tPEOjWjXo3MynVPYnHKNGjRpMmTKFM2fO0KRJxX4wLDU1lcjISCwWC1pr/v3vf+Pi4uKweKKjo3nvvfdKbX+S0IV9mvUEZYITq0oloZ9PSmfd0QtMiGhRKmVyhWMNHDjQ0SHYxcvLq8j7zsvT1TtlSov0oQv71KgNDTsYCb0ULNwVi03DnR3ylwUSQpSUJHRhv+YRELMV0m/+4Zyfd56hfSNvWtatmveeC+EIktCF/ZpHgrbCqfVFty3EkXNJ7I9NlLNzIUqZJHRhv8bh4FQDjq+8qd38vPMMZpPittCGpRSYEAIkoYvicHLNHpZuVYl3YbNpFu6KpWdLX/y8XIveoAp4/fXXCQoKIiQkhLCwMDZv3gzAhx9+mOeBmiFDhuSUaPX0NGrC36ikrM1mY9KkSQQHB9O+fXs6d+7MyZMny+HTFGzVqlVs2LAhZzl3Kdsbld3NX4BL3Dy5y0UUT4tIWPZPSDgD3sXvMtl5+jJnrqTxjwGBZRBcxbNx40YWLVrEjh07cHV15eLFi2RmZgJGQh8zZgzu7u4ALF68uLBd5fHtt98SGxvLnj17MJlMxMTE5Hmcv7ytWrUKT09PunfvDhilbEX5kzN0UTzNI4xpCZ8a/XV3HC5OJvq3q1dqIVVkcXFx+Pr65tQV8fX1pWHDhkyfPp3Y2FgiIyOJjIwEijcwRFxcHA0aNMCUXSzN39+f2rVrA9fO7gF++OEHHnroIQC+//57goODCQ0NpVcvo5yD1Wrl6aefpn379oSEhPDxxx8DsH37dnr37k3Hjh0ZOHAgcXFxgHG2/X//9390796d4OBgtmzZQlRUFDNnzuSDDz4gLCyMtWvX5illC0a52dzb5HfhwgXuvvtuOnfuTOfOnfOUxRX2kzN0UTx1g8DDz+hHD7uvWJtabZrFe+OICPQrtUGgi2XJ840HpFUAACAASURBVHC26NodxVK/PQx+64arBwwYwNSpUwkMDKRfv36MHj2a3r17M2nSJN5//31WrlyZU1WxOEaNGkXPnj1Zu3Ytffv2ZcyYMTcsI3vV1KlTWbp0KY0aNcrp2pk1axYnT55k586dODk5ER8fj8Vi4cknn2TBggX4+fnx7bff8tJLLzF79mwAUlJS2LBhA2vWrGHs2LHs27ePxx57DE9PT55++mnAqEiYW0Hb5PbUU0/xt7/9jZ49exIdHc3AgQM5ePBgsb8u1Z0kdFE8JpNRffHEKmNYumIMZLs1Kp7zSRnV6mKop6cn27dvZ+3ataxcuZLRo0fz1ltv5Zw1l5S/vz+HDx9mxYoVrFixgr59+/L999/Tt2/fG27To0cPHnroIUaNGpVTWnb58uU89thjODkZqcDHx4d9+/axb9++nDK/VquVBg2uFU+7Wvu8V69eJCYm2jU0W1HbLF++nAMHDuQsJyYmkpSUhJeX3NZaHJLQRfG1iIR9P8D5A1DP/pK0i/bEUsPZTN+2DhpmrpAz6bJkNpuJiIggIiKC9u3b8+WXX950QgejPOzgwYMZPHgw9erV45dffqFv3755SrXmLjs7c+ZMNm/ezG+//UZYWBi7du1Ca31daVetNUFBQWzcuLHA4+ZvX5xSvzdattlsbNy4sdTKyFZX0ocuiq95hDEtxt0uWVYbS/aepU/buri7VJ/ziMOHD3P06NGc5V27duVU2/Py8iIpKalE+92xY0fOoAg2m409e/bk7LdevXocPHgQm83Gzz//nLPN8ePH6dKlC1OnTsXX15fTp08zYMAAZs6cmVMvPD4+ntatW3PhwoWchG6xWNi/f3/Ofr799lvAGCzC29sbb2/vIj9LQdvkNmDAAD755JM8XydRfNXnf5YoPd7+UKeV0Y/e7Qm7Ntl44hKXUjK5LaTq1z3PLTk5mSeffJIrV67g5OREy5Ytc0arGT9+PIMHD6ZBgwasXFm8e/vPnz/PuHHjckYhCg8Pz7kF8K233mLYsGE0btyY4ODgnHrrzzzzDEePHkVrTd++fQkNDSU4OJgjR44QEhKCs7Mz48aNY+LEifzwww9MmjSJhIQEsrKy+L//+z+Cgoy/xmrXrk337t1JTEzM6Ve/7bbbGDFiBAsWLMi5sJpbQdvkNn36dJ544glCQkLIysqiV69ezJw5s1hfEwFKl3AEj5vVqVMnXZGK5Ihi+u1pYxSj506BU9HV6p77YQ+/7Y1j2z/74eZcfgNBHzx4kLZt25bb8aq6iIgIpk2bRqdOnRwdSrVQ0M+vUmq71rrAb4B0uYiSaREJllSIuf4WtPwsVhtLD5ylf7t65ZrMhahupMtFlEzucrrNehbadPOJeK6kWhgcXL98YhNlprgjE4nyJWfoomTcvKFRR7vquizZF4e7i5legQWPDC+EKB2S0EXJtegDsTsgNf6GTaw2zdL954hsXVe6W4QoY5LQRcm1GgDaBsf+vGGTHdGXuZicwUDpbhGizElCFyXXsAO414Gjy27YZMnes7g4mejTxkEPEwlRjUhCFyVnMkPL/nBsOdis163WWrN0/1l6tfLF07X6Xn8/e/Ys99xzDy1atKBdu3YMGTKEI0eOlGsMV65c4dNPP81ZvlFZ3htZtWoVw4YNK7Ld9OnTadu2Lffff3+xY8xfTrg046ouJKGLmxM4ANLiIeb6Zwr2nkngzJU0BgZV3+4WrTV33nknERERHD9+nAMHDvDGG29w7tw5u/dhteb9ZXn1qc7iyJ/Qy8qnn37K4sWLmTt3brG3LUlCF3lJQhc3p0UfUGY4uvS6Vb/vO4vZpKpNqdyCrFy5Emdn5zz1wcPCwrj11luvO7ucOHEic+bMAYxSulOnTqVnz558//33RERE8OKLL9K7d28++uijG5abffXVVxk7diwRERE0b96c6dOnA/D8889z/PhxwsLCeOaZZ/LEeOutt+Z51L5Hjx7s2bPnhp/pRsd47LHHOHHiBMOHD+eDDz4gJSWFsWPH0rlzZzp06MCCBQuAgkv2FlROeNmyZXTr1o1bbrmFkSNH5jzx+vvvv9OmTRt69uzJTz/9VKLvS1VVff8OFqWjRm1o3MXoR+/7cs7bWmt+33eWbs3rUMu96CdJy8PbW97mUPyhUt1nG582PBf+3A3X79u3j44dO5Zo325ubqxbtw4wCmtduXKF1auNOvT33XffDcvNHjp0iJUrV5KUlETr1q2ZMGECb731Fvv27ctJ3FFRUTnHeeSRR5gzZw4ffvghR44cISMjg5CQkEJjK+gYM2fO5Pfff88pCfziiy/Sp08fZs+ezZUrVwgPD6dfv3589dVX15Xs9fHxyVNO+OLFi7z22mssX74cDw8P3n77bd5//32effZZxo0bx4oVK2jZsiWjR48u0de2qpKELm5e4ABY/iokxkJNozTukXPJnLiYwtieAY6NrRLLn6xyL9+o3CzA0KFDcXV1xdXVlbp16xbZvTNy5Ej+9a9/8e677zJ79my7KkEWdAx/f/88bZYtW8bChQtzBrpIT08nOjq6wJK9+W3atIkDBw7Qo0cPADIzM+nWrRuHDh0iICCAVq1aATBmzJic2jhCErooDa0GGgn96DLo+BBgdLcoBQOCKk53S2Fn0mUlKCiIH374ocB1Tk5O2Gy2nOXcpW6B64aUy71cWLnZq6MjgVG6t6g+d3d3d/r378+CBQv47rvvChz/syTH0Frz448/0rp16+veL6rkrtaa/v37M2/evDzv79q1y65yvdWVXX3oSqlBSqnDSqljSqnnC2k3QimllVJSuac6qdsWvBvDkWu3Ly7ZF0enprWp6+XmwMAcr0+fPmRkZPD555/nvLd161ZWr15N06ZNOXDgABkZGSQkJFw3yk9hiltutqjyto888giTJk2ic+fOBZ4xl8TAgQP5+OOPuVoAcOfOnQAFluzNH2PXrl1Zv349x44dAyA1NZUjR47Qpk0bTp48yfHjxwGuS/jVXZEJXSllBmYAg4F2wL1KqXYFtPMCJgGbSztIUcEpZTxkdGIVWNKJupjCobNJDAquXqVyC6KU4ueff+aPP/6gRYsWBAUF8eqrr9KwYUMaN27MqFGjCAkJ4f777y9yCLncpk+fzrZt2wgJCaFdu3ZFlpqtU6cOPXr0IDg4+LqLogAdO3akZs2a/PWvfy32Z7yRyZMnY7FYCAkJITg4mMmTJwPGL48mTZoQEhJCaGgo//vf/4Br5YQjIyPx8/Njzpw53HvvvYSEhNC1a1cOHTqEm5sbs2bNYujQofTs2TOnBrwwFFk+VynVDXhVaz0we/kFAK31m/nafQgsB54GntZaF/p3m5TPrWKOLYdv7oZ75jHzXGveWnKIdc9F4l/b3aFhSflc+8TGxhIREcGhQ4dyBp4WjlcW5XMbAadzLcdkv5f7AB2AxlrrRYXtSCk1Xim1TSm17cKFC3YcWlQazXoZBbsOLmTJvrOE+Hs7PJkL+3z11Vd06dKF119/XZJ5JWfPd6+gKxA5p/VKKRPwAfCPonaktZ6lte6kte7k5yeV96oUJxdoPQTbocUcOH2xWj9MVNk8+OCDnD59mpEjRzo6FHGT7EnoMUDjXMv+QGyuZS8gGFillIoCugIL5cJoNdR2OKaMBLqZ9leo2ueOGpVLiJtRkp9bexL6VqCVUipAKeUC3AMszHXQBK21r9a6mda6GbAJGF5UH7qoglr0IU3V4B6PnTT383R0NIDxcM6lS5ckqYtKRWvNpUuXcHMr3l1iRd6HrrXOUkpNBJYCZmC21nq/UmoqsE1rvbDwPYjq4mKGYmNWGH2dthrFukyOr3/u7+9PTEwMcs1GVDZubm7XPaxVFLseLNJaLwYW53vv5Ru0jShWBKLKWLb/HGut4dxm2QinNkDArY4OCWdnZwIC5GlVUT3IJW1Ran7ff5YTtbqhnWrAgQWODkeIakcSuigVCakWNhy7SET7ZqiWfeHgr5DrsXYhRNmThC5KxZ+HzpFl0wwKqg9Bd0LyWYha6+iwhKhWJKGLUrF471kaeLsR6l8L2gwF15qwW+psCFGeJKGLm5ackcWaoxcYFFwfk0mBcw0IugMOLISMZEeHJ0S1IQld3LQVh86TmWVjSPtcxbhC7wNLChyUu1qFKC+S0MVNW7I3jrpernRsUvvam026Qu0A2PU/xwUmRDUjCV3clNTMLFYePn+tu+UqpSD0XuPC6JVoxwUoRDUiCV3clFWHL5BusTG4oNrnofcY093flm9QQlRTktDFTVm8N446Hi6EBxQwyk3tptC0J+z+H0gtFSHKnCR0UWLpFisrDp1nYHB9zKYbjPPY4X6IPwEn15RvcEJUQ5LQRYmtPnKB1EwrQwobai7oLnCvA5s/K7/AhKimJKGLEluyN47a7s50aV7IoMLObtDxr3B4McSfLL/ghKiGJKGLEsnIsvLnwfMMaFcfZ3MRP0adHzZK6W75vHyCE6KakoQuSmT9sYskZWQxuL0dIxPVbAjtboedX0NGUtkHJ0Q1JQldlMjivWep6eZE9xa+9m3QZQJkJMLu+WUbmBDVmCR0UWyZWTb+OHCOfu3q4eJk549Q487QqCNsnilldYUoI5LQRbGtPXqBhDQLQ9sXcndLQbpMgEvHjAukQohSJwldFNvC3bHUcnfm1lZ+xdsw6E7waQ6r35YHjYQoA5LQRbGkZmaxbP85hrRvYH93y1VmJ7j1aTi7B44uK5sAhajGJKGLYll+8DxpFivDQxuWbAcho6BWU1j5uvSlC1HKJKGLYlm46wz1a7oR3qyQh4kKY3aGyBchbjcc+Ll0gxOimpOELux2JTWT1UcucFtog7ylcour/UioGwR//gusltILUIhqThK6sNuSfWexWDW3hzW6uR2ZzNDvVbh8ErbMKo3QhBCAk6MDEJXHgl1naO7rQVDDmnneT7WksvvCbo5dOUZaVhq13WrTqlYr2vu2x2wyF7yzVv2hZX9Y9ZZxxu5Ztxw+gRBVmyR0YZe4hDQ2n4xnUp9WKGV0t5xNOct/9/6XBccXkJaVdt02Pm4+3NHyDh5s9yB1atTJu1IpGPQmfNoN/ngF7vx3eXwMIao0SejCLj9uj0FruPsWf7TWfH/ke97b9h6ZtkyGBAxhcMBgguoE4eHswaW0S+y+uJulJ5cyZ/8c5h2axxNhT3B/2/txMuX6kfNtBd2fhHXvG3e/tIh03AcUogpQ2kEPeHTq1Elv27bNIccWxaO1JmLaKhp4u/H1w52YsnEKC44voGuDrrzS7RX8vfxvuO3JhJNM2zaNNTFrCPML493e71LfI1dBL0sa/LsHaCtM2AAuHuXwiYSovJRS27XWnQpaJxdFRZG2nIzn1KVU7uhQjydXPMmC4wuYEDqBWf1nFZrMAQK8A/ikzye8eeubHLl8hFG/jmLX+V3XGjjXgOEfw+VTsGxyGX8SIao2SeiiSN9vj8HTVbHyyrusj13P1O5TeTzs8Zy+9KIopRjWfBjzh83Hy8WLccvGsSYm15B0zXpA94mw7b9w+Pcy+hRCVH12JXSl1CCl1GGl1DGl1PMFrP+7UuqAUmqPUupPpVTT0g9VOEJyRha/7YmlSeslbIhdx8vdXubOVneWaF8B3gF8NfgrmtdqzqQVk/j5aK4Hi/pMhvrt4ZfHZGQjIUqoyISulDIDM4DBQDvgXqVUu3zNdgKdtNYhwA/AO6UdqHCMxXviyPJazWnLKsaHjGdk4Mib2l+dGnWYPXA24fXDeXnDy3x76FtjhZMrjPramJ9/H2Qk32TkQlQ/9pyhhwPHtNYntNaZwHzg9twNtNYrtdap2YubgMI7VkWl8eWO1bjVXULfJn15IuyJUtmnh7MHM/rOIMI/gtc2v3btTN0nAEZ8ARcOGWfqUutFiGKxJ6E3Ak7nWo7Jfu9GHgaWFLRCKTVeKbVNKbXtwoUL9kcpHGL3mTiiTLPwcvZlao+pmFTpXXJxNjszLWIa3Rt255UNr/Dbid+MFS0iYcDrcPBXWP1WqR1PiOrAnv+hBV35KvBeR6XUGKAT8G5B67XWs7TWnbTWnfz8illLW5QrrTXPrZ6Mck7grZ5vU9OlZtEbFZOr2ZUPIz+kU/1OvLTuJVZErzBWdJ0AYWOMuukbPi714wpRVdmT0GOAxrmW/YHY/I2UUv2Al4DhWuuM0glPOMrcA99xxrKZQOeR9Gpa4C2vpaKGUw0+6fMJQXWCeGb1M2w9u9V4ivS2jyDoLlj2T9gkT5EKYQ97EvpWoJVSKkAp5QLcAyzM3UAp1QH4DCOZny/9MEV5ik2O5f3t75GV3JKXez1e5sdzd3ZnRt8ZNPZqzJMrnmT/pf3GYBh3zYK2t8Hvz8OmmWUehxCVXZEJXWudBUwElgIHge+01vuVUlOVUsOzm70LeALfK6V2KaUW3mB3ooLTWvPqhilYrDZamv5KWOMS1j0vplputfis/2d4u3gz4Y8JnEw4adROv3s2tB4Kvz8HS18Cm7Vc4hGiMrLrKpfWerHWOlBr3UJr/Xr2ey9rrRdmz/fTWtfTWodlv4YXvkdRUf1y7Bc2xm0g/fwgHunWsVyPXc+jHrMGzEIpxfg/xnM25Sw4ucCoryD8Udj4SfYtjUnlGpcQlYU8KSpynE89z7vb3sVTB+Kd1YvBwQ3KPYamNZsys99MkjOTGf/HeC6nXza6X4a8A0OmwdE/4D/94ezeco9NiIpOEroAjK6W1za9RkZWBmePD2dMl2bFHwS6lLSt05aP+3xMbHIsjy9/nBRLirEifByM+RHS4uHzPrD+I+mCESIXSegCgKVRS1l5eiVNTXfhSl3GdHVs9YZO9Tsxrfc0DsYf5KmVT5FpzTRWtIiECRuh1QD442WYM1TO1oXIJgldEJ8ezxub36CVdzt27w9hTJem+Hq6OjosIhpHMLXHVDbHbea5Nc9hvXo27lEHRn8Dd/wbLhyGz3rBor9ByiXHBiyEg0lCF7y5+U2SLEnUy3wAs8mJ8b2aOzqkHMNbDOfZzs+yPHo5UzdNJad+v1IQdh9M2mFcMN3+JXwUCsunQMpFxwYthINIQq/m/jz1J79H/c69rR7mj10m7gtvQt2abo4OK48H2j3A+JDx/HT0Jz7Y8UHelTVqw+C34PGNxjil6z6AD4Jh8bPG2bsQ1Ygk9GrsSvoV/rXpX7T1aUt8bA9MSvFo74pzdp7bxLCJjAocxRf7vmD2vtnXN/BrDSO/gCe2QNCdsG02zAiHL4bA7vlyq6OoFmRM0Wrs7a1vk5CRwL+6fsTYWWcY3bkxDbxrODqsAimleLHLiyRmJvLB9g+o5VqLu1rddX1Dv0BjwOn+U2DXXNg+B35+FJzcoGU/aHcHtOhj9MMLUcVIQq+mVp9ezaITi3gs9DF+3qxRSjEhoqWjwyqU2WTmjZ5vkJSZxJSNU6jhVIPBAYMLbuxZF3r+Dbo/Bac3wf5f4MACOLQIUNAgBJpHGnfNNO5iDIUnRCUng0RXQ4mZidz5y514u3nzcodZ3PnpZh7t1YLnB7dxdGh2SbWk8vifj7Pz/E6mdJ/CHS3vsG9DmxXO7IATK+H4SojZArYsUGao1w4a3gKNbjGmvoHgXLGuJQgBhQ8SLQm9Gvrnun+y6MQi5g6Zy79+SuTIuWRWPRNBTTdnR4dmt7SsNJ5a8RQb4zbyUpeXuKfNPcXfSUYSRK03EvuZHRC7E9KvGOuUCWo1NfrmfQOhTgvwbgy1moC3v5zRC4cpLKFLl0s1s+TkEhYcX8C49uM4FlOLTSeieP3O4EqVzMEou/tx3495etXTvL75dZIyk3ik/SN2D1wNgKsXtB5kvAC0hvgTELfbuEPm4mG4cASOr4CrDzZd5eFnJPaajcDDF9x9jfc8fMG9zrV5N2+j/744cQlRQpLQq5HTSaeZsnEKoX6hjGnzCAM/2ECovzf3dG7i6NBKxNXsyvuR7zN5/WSm75xOXEocL3Z5ESdTCX+slTLOxOu0yPu+NQuSYuHKaUjIfl2dv3QcTm+G1EugbzBknsnJ+OXh4mVMXb3A1fPavLO7kfSd3IyxVXNPnXO9b3Y1KlCanMFkNvZrdjamV5dN+Zavrlcm+aVSDUhCryYsNgvPrXkOEybe7vU27/5+nEvJGfznwU6YTZX3P7qzyZk3er5BA48G/Gfvfzifep53er2Du7N76R3E7GR0tdQq5BefzQZplyH1ovFgU8oFYz490ejayUw2phlJkJEIqfFw+ZQxb0mHrHSwlvG4MMpsJHVlArKnVxN9znuqgPfsbafIM8DZdb9A8i3nWV/Yunzrr/txLY9jlvL/kR5PGbX+S5kk9Gri4x0fs/fiXt7r/R4n4lyYtyWaR3s1J7RxLUeHdtNMysRTtzxFA48GvL75dR5Y8gAfRn5IY6/GRW9cakGYjFshPeoY/e4lYbMZXTtZaZCVYST53FNLmnER12bNnlryLltzL1+dz162Woy/ILQN0NnzuabXvVfCdjnyXZu77lqdLmSdndtdt/5mjlmM/ZYGU9l0cUpCrwYWn1jMF/u/YGTgSMLrRjD4o7U09/Pgb/0DHR1aqRrVehSNPBvx7JpnuWfRPbzT6x16NOrh6LDsZzKByU3urhElJk+KVnF7L+xl8vrJ3FL3Fp7v/DzP/LCHi8kZfDS6A27OZkeHV+p6NOrB/GHzqe9RnwnLJ/DvXf8my5bl6LCEKBeS0KuwsylnmbRyEn7ufnwQ+QFfbojhjwPneG5QG9r7ezs6vDLT2KsxXw/+mqHNh/Lp7k8Zu3QsZ5LPODosIcqcJPQqKiEjgSdXPElaVhqf9PmEvdFW3lxykMHB9Xm4Z4Cjwytz7s7uvHnrm7x565scvXyUEQtH8MuxX3DUcxdClAdJ6FVQUmYSj/3xGMevHOe93u+RmVaXiXN3EFjPi2kjQ4t3r3YlN6z5ML6/7XsCawcyef1kxv0xjujEaEeHJUSZkIRexaRYUpiwfAKH4g/xfsT7NHbrwF++2IKnmxOzH+qMh2v1uw7u7+XPF4O+YHLXyey/uJ+7Ft7FzN0zSctKc3RoQpQqSehVyOX0yzz2x2Psu7iPd3u/S4B7Z+6ZtRGL1caXY8NpWKv6Pq5uUiZGtR7FgjsW0Mu/FzN2zWDYz8P49fiv2G70QJAQlYwk9CriRMIJ7l98PwcuHeDd3u/SyCWcUZ9tJM1i5X+PdCWwnpejQ6wQ6rrX5f2I95kzaA5+Nfx4cd2LjPx1JMuilkliF5WeJPQqYGPsRsb8NoYUSwr/HfhfzGkhjJy5EYVi3viutGtY09EhVjgd63Xkf0P/x1u3vkWmNZN/rP4Hdy+8m1+P/4rFanF0eEKUiFRbrMQyrZnM2DWDOfvn0Ny7OR/2/phv1ify+dqTtGtQk//8pVO17maxl9VmZWnUUj7b8xknEk7gW8OXUa1HMTJwJL41fB0dnhB5SPncKuhQ/CFeXPciRy8f5e5WdzMiYAL//OkYu05f4YGuTXlpaNsq+eBQWbJpGxtiNzD34FzWnVmHWZnp1rAbtzW/jcgmkdRwkl+OwvEkoVchF1Iv8OnuT/n56M/4uPnwXKfJ7DhUn9nrT+LmbObtu0MY0r6Bo8Os9KISolhwfAGLTizibMpZPJw96NekH/2a9iO8fnjpFv8SohgkoVcBF9MuMu/QPL4+8DUWm4W7Wo7E1zKU/665wMXkDEZ29OeZQa2p6yV1QEqTTdvYfm47C48v5I9Tf5BiScHF5EJ4g3B6+feiZ8Oe+Hv5V6t7+4VjSUKvpLTW7L+0n3mH5rHk5BIsNgu9GvbDJ+N2FmzLIDE9i/BmPrw0tG2VqJpY0VmsFraf387q06tZE7OG6CTjAaW6NepyS71b6FC3Ax3rdaRFrRYlr8kuRBEkoVciNm3jcPxh/jj1B0ujlhKdFE0NJ3faevbhclxn9kS5ohQMCqrPI7c2p2PT2o4OudqKSohic9xmtp/fzo5zOziXeg4wBt4IrB1Ia5/WtPVpS2DtQAK8A/B2rbr1c0T5kYRegVmsFo5eOcq+i/vYenYrm89u5nL6ZUyY8HUKIu1KMLFnWoPNjVZ1PbmjQyOGhzaksY/04VYkWmtiU2LZcW4HB+MPcij+EIfiD5GUmZTTxtvVm6Y1m9KsZjOa1mxKE68m1PeoTz33evi6++JcRjWyRdVy0wldKTUI+AgwA//RWr+Vb70r8BXQEbgEjNZaRxW2z+qU0G3axuX0y5xJPkN0UjTRidEcvxLF8csniEo6gVUb9z076VrYUluSkhCANbkNLsqL8AAferb0pXdrP1rX85K+2krkapI/evkopxJP5byiEqM4n3o+T1uFwreGL3Xd61LXvS4+bj7Ucq1FbbfaeaeutfFy8cLD2QNns/wCqI5uapBopZQZmAH0B2KArUqphVrrA7maPQxc1lq3VErdA7wNjL750MufTduw2CxkWjNzpumWTNKyMknLyiAtM4N0q4VUSwbJmSkkW5JJykgmMTOJhIxkkjOTSbYkk2JJIskST4r1Mum2BDTXnkLUWqEt3tgy/bBldMOa1ogauimNazellZ8XoR28CfGvRdsGXrg6ya2HlZVSikaejWjk2ei6damWVGKSYzifep5zKec4l5r9SjlHTHIM+y7u43LG5UJrubuYXPBw9sDd2R0PZ49r807G1MXkgovZBVezK65mV5zNzjnzLmYXXEwu1+bNLpiVGbPJnHeaPe+knDApE04mJ8zKXOC8SZlQKJRSOVNRvuy5chMOHNNanwBQSs0HbgdyJ/TbgVez538APlFKKV0G/TkvLPucJTHz0GQPh5Xrlf89jQaVvx3ZU1v2kIG51qmbC1drM9rqBjZXtM0NneWFsrbEGW+8nGpT26Uu9d398fdqTD0vT/y8XGlWx50AXw98PFzkP0A14u7sTmDtQAJr33jUKK01KZYULmdc5kr6FS5nXOZy+mWSMpNIsaSQkpVCqiXVmLcY8wnpCcRmxZJqtRSr8wAABuFJREFUScVis5Bhzcg5OXGEq4ndhMkYdjT7n0mZcn7ec34RGA3yLOf+5VDQNP+xrjt+Ef+n7NnmuuMUsM+C9lPYfieETmBwwOAitykuexJ6I+B0ruUYoMuN2mits5RSCUAd4GLuRkqp8cB4gCZNSjbSfD0PP3ycWuT6Acj9zTdhDFNruvaNz/nmGwPZGj9cKnvehCl7YNvsHzNMSmFWzjkvF7MTztlnOi5m5+wzHmPe3ckdTxdPPJw9+P/2zi5EqjKM47+/rrt+rB+ZFkKla1jgVS5dCKY3RamU2icbQULdBAktEWQIEt1Z1EUQiaFkYiVR0l4YGRV1paW26wdqahmY2xoGGfSx7c7TxXlnmR3nrDO7zntmpucHh/POM++Z8+d53/PMe97z8cycNI1pLZOYOnECU1rG09rSxOTmJpqb/O0KzuiQRGtzK63NrWPOj5qzHP2D/fTn+pP1YP9QsM+Xc5ZjwAYYzA0OKw9aWEJ5IDdAznJD5UEL9XMDySDKkuFUznLJsCqM64Y+hzqFn82Gr/Pv1TEzcuSGfmPod4vGilZG3s9ytim2lbWf4nSkZfxutS6QlxPQS/31FCsupw5mtgXYAskcehn7vozOJWvoXLJmNJs6zv+WcRrHxKaJTMSfU2hkyhk+ngMKhwc3AOfT6khqAqYDv10NgY7jOE55lBPQvwUWSGqT1Ax0AF1FdbqAtaH8EPBFNebPHcdxnHSuOOUS5sTXAZ+S3La4zcyOSXoJOGBmXcBWYIek0yQj845qinYcx3Eup6znk81sD7CnyLaxoPw38PDVleY4juNUgt+C4TiO0yB4QHccx2kQPKA7juM0CB7QHcdxGoTM3rYo6Vfgp1FuPouip1BrBNdVGa6rcmpVm+uqjLHommtms0t9kVlAHwuSDqS9bSxLXFdluK7KqVVtrqsyqqXLp1wcx3EaBA/ojuM4DUK9BvQtWQtIwXVVhuuqnFrV5roqoyq66nIO3XEcx7mceh2hO47jOEV4QHccx2kQ6i6gS1ou6aSk05LWZ6jjRklfSjou6ZikZ4L9RUk/S+oOy8oMtJ2VdCTs/0CwzZT0maRTYX1NZE23FvikW9IlSZ1Z+EvSNkkXJB0tsJX0jxJeD/3tsKT2yLpekXQi7Hu3pBnBPk/SXwV+2xxZV2q7SXoh+OukpHsi69pVoOmspO5gj+mvtNhQ/T5mZnWzkLy+9wwwH2gGeoCFGWmZA7SH8lTge2AhSW7V5zL201lgVpHtZWB9KK8HNmXcjr8Ac7PwF7AMaAeOXsk/wErgE5KsXIuB/ZF13Q00hfKmAl3zCutl4K+S7RaOgR6gBWgLx+v4WLqKvn8V2JiBv9JiQ9X7WL2N0IcSVptZP5BPWB0dM+s1s0Oh/AdwnCS3aq2yGtgeytuBLPP43QmcMbPRPik8Jszsay7PqJXmn9XAO5awD5ghaU4sXWa218wGwsd9JBnDopLirzRWA++b2T9m9iNwmuS4japLkoBHgPeqse+RGCE2VL2P1VtAL5WwOvMgKmkesAjYH0zrwqnTtthTGwED9ko6qCQxN8D1ZtYLSYcDrstAV54Ohh9oWfsL0v1TS33uCZKRXJ42Sd9J+krS0gz0lGq3WvHXUqDPzE4V2KL7qyg2VL2P1VtALysZdUwktQIfAp1mdgl4E7gZuA3oJTnti80SM2sHVgBPS1qWgYaSKEljuAr4IJhqwV8jURN9TtIGYADYGUy9wE1mtgh4FnhX0rSIktLarSb8BTzK8EFDdH+ViA2pVUvYRuWzegvo5SSsjoakCSQNttPMPgIwsz4zGzSzHPAWVTrdHAkzOx/WF4DdQUNf/jQurC/E1hVYARwys76gMXN/BdL8k3mfk7QWuBd4zMKka5jSuBjKB0nmqm+JpWmEdqsFfzUBDwC78rbY/ioVG4jQx+otoJeTsDoKYY5uK3DczF4rsBfOfd0PHC3etsq6pkiami+TXFQ7yvBE3muBj2PqKmDYyClrfxWQ5p8u4PFwJ8Ji4Pf8aXMMJC0HngdWmdmfBfbZksaH8nxgAfBDRF1p7dYFdEhqkdQWdH0TS1fgLuCEmZ3LG2L6Ky02EKOPxbjqe5WvIK8kuWp8BtiQoY47SE6LDgPdYVkJ7ACOBHsXMCeyrvkkdxn0AMfyPgKuBT4HToX1zAx8Nhm4CEwvsEX3F8kfSi/wL8no6Mk0/5CcDr8R+tsR4PbIuk6TzK/m+9jmUPfB0L49wCHgvsi6UtsN2BD8dRJYEVNXsL8NPFVUN6a/0mJD1fuYP/rvOI7TINTblIvjOI6Tggd0x3GcBsEDuuM4ToPgAd1xHKdB8IDuOI7TIHhAdxzHaRA8oDuO4zQI/wE7UZzEw1HlywAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "Results = [R_pop, S_pop,I_pop]\n", "Results_df = pd.DataFrame(Results).transpose()\n", "Results_df.columns = [\"Recovered (Including Dead)\", \"Still Susceptible\", \"Currently Infected\"]\n", "\n", "Results_df.plot(title=\"SIR Model\")" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.10087844165170001" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "S_tavg" ] } ], "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.7" } }, "nbformat": 4, "nbformat_minor": 4 }