Skip to content

Instantly share code, notes, and snippets.

@sellisd
Last active April 6, 2020 09:19
Show Gist options
  • Save sellisd/8e2ae694641d6f56db9bc02273090f3b to your computer and use it in GitHub Desktop.
Save sellisd/8e2ae694641d6f56db9bc02273090f3b to your computer and use it in GitHub Desktop.
a generic SIR model
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# An extensible SIR model\n",
"\n",
"We could add some more population compartments (states) to get a more complex model and modify the corresponding equations"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#!/usr/bin/env python\n",
"import pandas as pd\n",
"\n",
"class Population:\n",
" \"\"\"A class representing a human population randomly mixing\"\"\"\n",
" def __init__(self, states, parameters):\n",
" self.states = states\n",
" self.parameters = parameters\n",
"\n",
" def time_step(self):\n",
" \"\"\"Discrete time step (days?). Update number of people in each state\"\"\"\n",
" S = self.states['Susceptible']\n",
" I = self.states['Infected']\n",
" R = self.states['Recovered']\n",
" N = S+I+R\n",
" a1 = self.parameters['beta'] * I * S / N\n",
" a2 = self.parameters['gamma'] * I\n",
" self.states['Susceptible'] += -a1\n",
" self.states['Infected'] += a1 - a2\n",
" self.states['Recovered'] += a2\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Initialize the model with some arbitrary parameters"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# example of states\n",
"states = {\n",
" \"Susceptible\": 1000,\n",
" \"Infected\": 1,\n",
" \"Recovered\": 0}\n",
"\n",
"# example of parameters\n",
"parameters = {\n",
" \"beta\": 0.2, # average number of contacts per person per time\n",
" \"gamma\": 0.1 # rate of recovery or mortality\n",
"}\n",
"\n",
"# setup model:\n",
"france_simulation = Population(states, parameters)\n",
"\n",
"simulation_time = 180 # in days\n",
"historical_data = pd.DataFrame(columns = states.keys(),index = range(simulation_time))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Run simulation for three months and plot evolution"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.axes._subplots.AxesSubplot at 0x7fb8f881a240>"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3gU1frA8e/Zkmx6pSaB0HsIELoFAgo2/NmwoBdsWK7t2kWRa716saFiv4oNUKyIqCiCIkgLvdcICSVAIJC+5fz+mE1IMLSUnU3yfp5nnpk5M7v75rC8M3vmzBmltUYIIUT9YDE7ACGEEL4jSV8IIeoRSfpCCFGPSNIXQoh6RJK+EELUIzazAziR2NhYnZiYaHYYQghRq6Slpe3XWjeoaJtfJ/3ExESWLl1qdhhCCFGrKKX+Ot42ad4RQoh6RJK+EELUI5L0hRCiHpGkL4QQ9YgkfSGEqEdOmvSVUu8rpbKUUmvKlEUrpX5WSm32zqO85Uop9apSaotSapVSqnuZ14z07r9ZKTWyZv4cIYQQJ3IqZ/qTgKHHlD0MzNZatwFme9cBzgPaeKfRwJtgHCSAcUBvoBcwruRAIYQQwndO2k9fa/27UirxmOKLgQHe5Q+BucBD3vKPtDFe80KlVKRSqol335+11tkASqmfMQ4kU0702XsPF/LG3C0E2qwE2izGZC+zbLMSaC+zbLPgsFsJDbThsFtQSp1iNQghRP1Q2ZuzGmmtd3uX9wCNvMtxwM4y+2V4y45X/jdKqdEYvxIIaNya//64sVIBWi2K0EAboYE2whzGPNRxdD0yOICYkABiQgOIDgkssxxAoM1aqc8UQgh/V+U7crXWWilVbU9i0Vq/A7wDkJKSov94aihFLg9FLjdFTs/RZZfHu+5ddnkocropdLrJLXKTW+Qkt9DFkSIXuYUucotcZOcVs+NAPocLXRzKL8blqTjscIeNppFBNIlw0DQyyDs5aBYdTMvYUKJCAqrrzxVCCJ+qbNLfq5RqorXe7W2+yfKWZwIJZfaL95ZlcrQ5qKR87ql8kMNuxWG3AvZKhloxrTWHC1wcyCsiO6+Y/bnFZOcVk51XRNaRInYdKmR3TgErdh7iYL6z3Gujgu20bBBKi9gQWjUIpVPTcLrERcjBQAjh9yqb9KcDI4HnvPNvy5TfoZSainHRNsd7YPgJeLbMxdtzgUcqH3bVKaWICLYTEWynZYXDEh1VUOwm81ABO7Lz2LYvj6378ti2L5ffN+3ji7SM0v3iIoPoEhdBl/gIkhMi6d4siqAAaSoSQviPkyZ9pdQUjLP0WKVUBkYvnOeAz5VSNwJ/AcO9u88Ezge2APnA9QBa62yl1FPAEu9+T5Zc1K0NggKstG4YSuuGoaS2L78tJ9/J2l05rM40pjWZOfy4dg8Adquia3wkfVrGcGabWHo0j8JmlVsjhBDmUf78YPSUlBRdG0fZzClwsnzHQRZuy2bhtgOszszB7dGEO2wMaNeQQR0aktq+IWGO6m2yEkIIAKVUmtY6paJtfj20cm0VEWRnQLuGDGjXEIAjhU7mb9nP7PVZzNmYxfSVuwi0WRjcsREXd23K2e0aSI8hIYRPSNL3gTCHnaGdmzC0cxM8Hs3ynQf5dsUuZqzazferdhMZbGd4SgLX9m5Os5hgs8MVQtRh0rxjIqfbwx9b9vP5kp3MWrcXj9YMbNeQ6/snckbrWLm5TAhRKSdq3pGk7yd25xQwedEOpizeyf7cIro1i+TuQW04u20DSf5CiNMiSb8WKXK5mbY0gzfnbiXzUAFdEyJ5aEg7+rWONTs0IUQtcaKkL/0H/Uygzcq1fZoz5/4B/OfSLuw/UsQ17y3ilo+XsjM73+zwhBC1nCR9PxVgs3B1r2bMvu9sHhjSjt837WfQS78x/qcNFDrdZocnhKilJOn7OYfdyj8HtmbO/QM4v3NjJs7ZyvmvzmP5joNmhyaEqIUk6dcSjSMcvHJVNz65sTeFxW4ue3MB43/aQLHLY3ZoQohaRJJ+LXNGm1h+/NdZXNY9nolztnLxxPmk788zOywhRC0hSb8WCnfYGX9FV977Rwq7cwq46PU/+HXDXrPDEkLUApL0a7HBHRvx3R1nkBAVzI0fLmXCL5vxHOcZAUIIAZL0a72E6GC+ur0fl3SL4+VfNnHrJ2nSu0cIcVyS9OsAh93Ki1d0ZdxFHfl5/V6u+98ico558IsQQoAk/TpDKcX1/Vvw+tXdWbkzh+Fv/8menEKzwxJC+BlJ+nXMBUlNmHR9TzIPFXDZmwvYti/X7JCEEH5Ekn4d1K91LFNH96HQ6eaadxex44AM3yCEMEjSr6M6x0Xw6c29KXS5uea9hew6VGB2SEIIPyBJvw5r3zicj2/oTU6Bk2veXUjWYWnjF6K+k6Rfx3WJj2DS9b3IOlLEiPcWcTCv2OyQhBAmkqRfD/RoHsX7o3ryV3Y+t3ycRpFL+vELUV9J0q8n+rSM4YUrurI4PZtHvlyNPz88RwhRc+TB6PXIsK5N+Wt/Hi/+vIkWsSHcOaiN2SEJIXxMkn49c0dqa7Z7E3+zmGAuTo4zOyQhhA9J8049o5TiP5d1oVeLaB74YhVrd+WYHZIQwock6ddDgTYrb47oTnRwALd/uozDhTJOjxD1hST9eiomNJDXr+lGxsECHpi2Ui7sClFPSNKvx1ISo3nkvPb8tHYv//tju9nhCCF8QJJ+PXfjGS0Y0qkRz/2wgaXp2WaHI4SoYZL06zmlFP+9vCtxUUHcPXWFtO8LUcdJ0hdEBNl5+cpkducU8MT0dWaHI4SoQVVK+kqpfyml1iql1iilpiilHEqpFkqpRUqpLUqpz5RSAd59A73rW7zbE6vjDxDVo3uzKP45sDVfLsvgxzV7zA5HCFFDKn1zllIqDrgL6Ki1LlBKfQ5cBZwPvKy1nqqUegu4EXjTOz+otW6tlLoKeB64ssp/gag2dw1qw5yNWYz5ejXdm0fSMMxhdkhCmMqjPbg8LlweF06PE6fHaSy7nTi1E6fbWbrNrd14tMeYezzl1yuae8qvl0zH7ufRHrTWePDOtQeNLi1DG3GWbNecuCdeVe/ItQFBSiknEAzsBlKBa7zbPwT+jZH0L/YuA3wBvK6UUlr6CvoNu9XCy8OTufC1P3jky9W8NzIFpZTZYQlRjtaaAlcBec48cp25R+fFxrzQVUihu5BCVyFF7qLyy2Xmhe6jy06PN4Frb0L3JneXdpn955ajUFiUBYVCKXV0vewyJ/4/W+mkr7XOVEq9AOwACoBZQBpwSOvSmsoASu7zjwN2el/rUkrlADHA/nJ/lFKjgdEAzZo1q2x4opLaNArjoaHteXLGOqYtzWB4zwSzQxJ1lNPj5FDhIQ4VHTMdU5ZTlFMusee58vBozyl9hkVZcFgdOGwOHFYHgbbA0vUgWxBRgVEE2gIJsARgt9qxW+zYLDbslhMsW48ul51blRWrsmJRFqwW71ydYG4x5gp1wv2VUqXzU6VGHH/fqjTvRGGcvbcADgHTgKGVfb8SWut3gHcAUlJS5FeACUb1S+THtXt4+vt1DGjfQJp5xGnLLc4lMzeTXbm7yMrPYl/BPvYV7CMrP4v9BfvJys/iYOHB4zZFBNmCiAyMJDIwkvDAcBoENSDEHkJoQKgxtxvzssuhAaGE2EIIsgcRaDWSu81ik1+rx6hK885gYLvWeh+AUuoroD8QqZSyec/244FM7/6ZQAKQoZSyARHAgSp8vqghFoviP5d24bwJ83jyu3W8fk13s0MSfkZrzd78vWw7tI30w+nsyt1FZm5m6XS4+HC5/S3KQowjhtigWBoFN6JTTCcaBjck2hFNpCOSqMCo0iQf6Ygk0Bpo0l9W91Ul6e8A+iilgjGadwYBS4E5wOXAVGAk8K13/+ne9T+923+V9nz/1apBKHcObM2LP2/i0u57SW3fyOyQhAk82kPmkUy2HNrCtpxtxnTImOe78kv3C7QG0jS0KXGhcSQ1SCIuNK50vVFwI6Id0VgtVhP/ElFCVSXvKqWewOiB4wKWAzdhtN1PBaK9ZddqrYuUUg7gY6AbkA1cpbXedqL3T0lJ0UuXLq10fKJqil0eLnxtHrmFLmbdezahgTISd13m0R52HtnJugPrSqf1B9ZzxHmkdJ+GQQ1pGdmSlhHeKbIlieGJxAbFSjOKH1FKpWmtUyrc5s8n25L0zZf210Euf2sBo/olMu6iTmaHI6qR0+Nk3YF1LN+7nLSsNJZnLSenyBhqO8ASQNuotnSM6UiHmA60iWpDy4iWhAWEmRy1OBUnSvpy6iZOqEfzKK7r05xJC9K5pFscSfGRZockKklrzaaDm5iXOY+Fuxayct9KCt2FADQPb05qQipdG3SlU2wnWkW2wm6xmxyxqAmS9MVJPTCkHTNX7+Hxb9fy1W39sFjkZ3xtcbj4MAt3LeSPzD+YnzmfrIIsANpGteXSNpfSo1EPujfqTmxQrMmRCl+RpC9OKsxhZ8z57bn385V8sSyD4SnSd9+fHSk+wtydc/kx/UcW7FqAy+MizB5G36Z9OSPuDPrH9adhcEOzwxQmkaQvTskl3eKYvGgHz/+wgSGdGhMRJD/9/Um+M7800c/PnE+xp5jGIY0Z0X4Eqc1SSWqQhM0i/92FJH1xipRS/HtYJy56/Q9e/nkT/x4mF3XNprVm7YG1fLn5S37Y/gN5zjwaBjVkeLvhDEkcQlKDJCxKBtIV5UnSF6esc1wEI3o346M/07myZwIdmoSbHVK9VOAqYMa2GUzdMJVNBzfhsDo4N/Fc/q/1/9GjUQ9J9OKEJOmL03L/ue34ftVuxk1fy2ej+0jfbB/ak7eHzzZ+xrRN08gpyqF9dHvG9hnLeS3Ok66U4pRJ0henJTI4gPvObcdj36zhp7V7Gdq5sdkh1Xk7D+/kvTXvMX3LdDx4SE1I5dqO19K9YXc56IrTJklfnLareibw4YJ0nvthPantGxJgk+aEmrA9ZzvvrnqXmdtnYlVWLm97OSM7jSQ+LN7s0EQtJklfnDab1cKY8ztw/aQlfLroL67v38LskOqUffn7eGPlG3y9+WsCrAGM6DCCUZ1G0SC4gdmhiTpAkr6olAHtGnBG61gmzN7Mpd3iiQiWLpxVlefM44M1H/DRuo9wup0MbzecW5JuISYoxuzQRB0iv8tFpSilGHN+B3IKnEycu8XscGo1rTU/bv+RYV8P4+1Vb3NW/Fl8+3/fMqb3GEn4otrJmb6otI5Nw7m8ezyT5qdzXZ/mJEQHmx1SrbMtZxv/WfQfFu5eSIfoDrw44EWSGyabHZaow+RMX1TJfee2w2pRPPfjBrNDqVWcHidvrXyLy6Zfxtr9axnTewxTLpgiCV/UOEn6okoaRzi4+ayWfL9qN2l/HTQ7nFph88HNjPh+BBNXTOScZucw/ZLpXN3+annIiPAJSfqiym45qyUNwgJ55vt1+PPzGczm9rh5b/V7DJ8xnL35e3lpwEv89+z/ygiXwqck6YsqCwm0cd85bVm24xAzV+8xOxy/tL9gP7f8fAsTlk0gNSGVry/+mnOan2N2WKIekqQvqsUVKQm0bRTK+J824HR7zA7HryzYtYDLpl/Gyn0rebLfk7xw9gtEO6LNDkvUU5L0RbWwWhQPDW1P+oF8pi7ZaXY4fsGjPUxcMZFbf76VaEc0Uy6YwiVtLpGhE4SpJOmLapPaviG9WkQz4ZfN5BW5zA7HVHnOPO6Zcw9vrXyLYa2GMfmCybSOam12WEJI0hfVRynFw+e1Z39uEe/N2252OKbZeXgn1868lt8zfufhXg/zVP+nCLIFmR2WEIAkfVHNujeLYminxrzz+1b25xaZHY7PLdmzhKtnXk1WfhZvDn6TER1GSHOO8CuS9EW1e2BoOwpdHl7/tX4NzzArfRa3/HxLaft936Z9zQ5JiL+RpC+qXasGoQxPSeDTRX/x14E8s8PxiSkbpnD/b/fTKaYTH5/3Mc3Cm5kdkhAVkqQvasS/BrfBZrHwwqxNZodSo7TWvLrsVZ5d9CxnJ5zNu+e+S0RghNlhCXFckvRFjWgY7uDGM1rw3cpdrM7IMTucGqG15j+L/8O7q9/l8raX8/KAl3HYHGaHJcQJSdIXNWb02S2JCrbzfB0cjM2jPTyz6BmmbJjCyI4jebzP49gsMmit8H+S9EWNCXfYuSO1DX9s2c+8zfvMDqfaeLSHpxY+xWcbP+OGzjdwX8p90kNH1BqS9EWNurZPM+Kjgnjuhw14PLV/MDaP9vDEn0/wxaYvuLnLzdzT/R5J+KJWkd+jokYF2qzcf2477vlsBd+t2sXFyXFmh1RpWmvGLxnPV5u/YnTSaO5IvkMSfhlOp5OMjAwKCwvNDqXecDgcxMfHY7ef+uNKq5T0lVKRwHtAZ0ADNwAbgc+ARCAdGK61PqiM/x0TgPOBfGCU1npZVT5f1A7Dujbl7d+38cKsjQzt3JhAW+0cN/6tlW/xyfpPuLbDtZLwK5CRkUFYWBiJiYlSNz6gtebAgQNkZGTQokWLU35dVZt3JgA/aq3bA12B9cDDwGytdRtgtncd4DygjXcaDbxZxc8WtYTFYgzPsDO7gMmLdpgdTqV8su4T3lj5Bhe3upgHej4gSa0ChYWFxMTESN34iFKKmJiY0/5lVemkr5SKAM4C/gegtS7WWh8CLgY+9O72IfB/3uWLgY+0YSEQqZRqUtnPF7XLWW1i6dcqhtd+3cKRQqfZ4ZyW77Z+x/NLnmdQs0H8u9+/sSi5FHY8kvB9qzL1XZVvbwtgH/CBUmq5Uuo9pVQI0Ehrvdu7zx6gkXc5Dig75m6Gt6wcpdRopdRSpdTSffvqTo+P+q5kMLbsvGLe/X2b2eGcskW7F/H4/Mfp3bg3z5/1vHTLFLVeVZK+DegOvKm17gbkcbQpBwBtPDvvtLpsaK3f0VqnaK1TGjRoUIXwhL9Jio/kwqQmvDtvO1lH/P9i37ZD2/jXnH/RPLw5Lw18iUBroNkhiVPwzDPP0KlTJ5KSkkhOTmbRokWmxLFixQpmzpxZuj59+nSee+45AEaNGsUXX3zxt9fMnTuXCy+8sEbjqkrSzwAytNYlNfoFxkFgb0mzjXee5d2eCSSUeX28t0zUI/ef2w6n28OEXzabHcoJ7S/Yz+2zbyfAGsDEwRMJDwg3OyRxCv78809mzJjBsmXLWLVqFb/88gsJCQknf2ENODbpDxs2jIcffvgEr/CNSid9rfUeYKdSqp23aBCwDpgOjPSWjQS+9S5PB/6hDH2AnDLNQKKeSIwN4ZrezZi6ZCfb9uWaHU6FCl2F3P3r3RwoOMBrqa8RF1p7u5nWN7t37yY2NpbAQONXWWxsLE2bNiUxMZH9+/cDsHTpUgYMGADAb7/9RnJyMsnJyXTr1o0jR44A8Pzzz9OlSxe6du1amqi3bt3K0KFD6dGjB2eeeSYbNhh3mo8aNYpbb72VlJQU2rZty4wZMyguLubxxx/ns88+Izk5mc8++4xJkyZxxx13lMb6yy+/lHvNsfLy8rjhhhvo1asX3bp149tvv/3bPpVR1QbKO4FPlVIBwDbgeowDyedKqRuBv4Dh3n1nYnTX3ILRZfP6Kn62qKXuTG3DF2kZvDBrI2+M6GF2OOVorRm3YByr96/mpQEv0aVBF7NDqpWe+G4t63Ydrtb37Ng0nHEXdTrhPueeey5PPvkkbdu2ZfDgwVx55ZWcffbZx93/hRdeYOLEifTv35/c3FwcDgc//PAD3377LYsWLSI4OJjs7GwARo8ezVtvvUWbNm1YtGgRt99+O7/++isA6enpLF68mK1btzJw4EC2bNnCk08+ydKlS3n99dcBmDRpUrnPrug1ZT3zzDOkpqby/vvvc+jQIXr16sXgwYMJCQk53aorp0pJX2u9AkipYNOgCvbVwD+r8nmibmgQFsjNZ7ZkwuzNLN9xkG7NoswOqdTH6z5m5vaZ3NntTgY3H2x2OOI0hYaGkpaWxrx585gzZw5XXnllaTt6Rfr378+9997LiBEjuPTSS4mPj+eXX37h+uuvJzg4GIDo6Ghyc3NZsGABV1xxRelri4qOPiRo+PDhWCwW2rRpQ8uWLUt/BZzIyV4za9Yspk+fzgsvvAAYXWJ37NhBhw4dTqtOjiVdEYQpbj6rJZ8u+ovnftjA1NF9/KKr36Ldi3gp7SUGNxvMzV1uNjucWu1kZ+Q1yWq1MmDAAAYMGECXLl348MMPsdlseDwegHL92h9++GEuuOACZs6cSf/+/fnpp58qfE+Px0NkZCQrVqyocPux399T+T6f7DVaa7788kvatWtHdZIOx8IUoYE27hrUhkXbs5m7yfyuubtyd3H/b/eTGJ7I02c87RcHIXH6Nm7cyObNRzsJrFixgubNm5OYmEhaWhoAX375Zen2rVu30qVLFx566CF69uzJhg0bOOecc/jggw/Iz88HIDs7m/DwcFq0aMG0adMAIyGvXLmy9H2mTZuGx+Nh69atbNu2jXbt2hEWFlZ6jaAiFb2mrCFDhvDaa69hNJLA8uXLq1g7Bkn6wjRX9WxG85hgnv9hA24TB2MrdBVyz5x7cHvcTEidQIi9am2mwjy5ubmMHDmSjh07kpSUxLp16/j3v//NuHHjuPvuu0lJScFqPToMyCuvvELnzp1JSkrCbrdz3nnnMXToUIYNG0ZKSgrJycmlzSuffvop//vf/+jatSudOnUqd2G1WbNm9OrVi/POO4+33noLh8PBwIEDWbduXemF3GNV9Jqyxo4di9PpJCkpiU6dOjF27NhqqSNVchTxRykpKXrp0qVmhyFq0Hcrd3HnlOW8eEVXLusRb0oMJaNmvp76OmcnHP+inzix9evXV7m9uTYaNWoUF154IZdffrkpn19RvSul0rTWFV1vlTN9Ya4LujShS1wEL/28iUKn2+efP3PbTL7Y9AU3dL5BEr6oFyTpC1OVDMaWeaiATxb+5dPPTs9J54k/nyC5QTJ3dLvj5C8QogKTJk0y7Sy/MiTpC9P1bx3LWW0b8PqcLeQU+GYwtiJ3Eff/dj92q53xZ4/Hbjn18ciFqM0k6Qu/8NDQdhzKd/LWb1t98nnjl4xn48GNPNP/GRqHNPbJZwrhDyTpC7/QqWkE/5fclPf/2M6enJodjG3uzrl8tvEzRnYcKe34ot6RpC/8xn3ntkNr+O+PJ7+bsbL2F+xn3IJxtItqx13d76qxzxHCX0nSF34jITqYm85swVfLM1m242C1v3/JuDq5xbk8d+ZzBFgDqv0zhLlCQ0NPus+8efPo1KkTycnJFBQUnNb7f/PNN6xbt65G4vIVSfrCr/xzYGsahQfyxPS1eKr5hq1pm6bxe8bv3JtyL62jWlfre4va49NPP+WRRx5hxYoVBAUFndZrK5v0/YkkfeFXQgJtPHxee1Zm5PDFsoxqe9/tOdsZv2Q8/Zr24+r2V1fb+wr/NHfuXAYMGMDll19O+/btGTFiBFpr3nvvPT7//HPGjh3LiBEjABg/fjw9e/YkKSmJcePGlb7HRx99RFJSEl27duW6665jwYIFTJ8+nQceeIDk5GS2bt163OGWt2/fTt++fenSpQuPPfaYKXVwPDLgmvA7/5ccx8d//sV/f9zA0M6NCXdUrTuly+Pi0T8eJdAWyFP9n5Jn3PrCDw/DntXV+56Nu8B5xx8x81jLly9n7dq1NG3alP79+zN//nxuuukm/vjjj9I7aGfNmsXmzZtZvHgxWmuGDRvG77//TkxMDE8//TQLFiwgNjaW7OxsoqOjGTZsWLm7bwcNGlThcMt33303t912G//4xz+YOHFi9dZDFcm3X/gdpRT/HtaJA3nFvDa76k/Y+njdx6zev5oxvcbQMLhhNUQoaoNevXoRHx+PxWIhOTmZ9PT0v+0za9YsZs2aRbdu3ejevTsbNmxg8+bN/Prrr1xxxRXExsYCxvDKxyo73HJycjK33HILu3cbz4WaP38+V19t/KK87rrrau6PrAQ50xd+KSk+kuE9EvhgfjrDUxJo0yisUu+zPWc7ry9/nYEJAzmvxXnVHKU4rtM4I68pJU/PAmO4ZZfL9bd9tNY88sgj3HLLLeXKX3vttZO+/+kOt+wv5Exf+K0Hh7YjJNDGo9+soTIDA7o9bh6f/zgOm4Oxfcb67X9CYZ4hQ4bw/vvvk5trPLozMzOTrKwsUlNTmTZtGgcOHAAofXpW2eGSTzTccv/+/Zk6dSpgXDj2J5L0hd+KCQ3kkfPas3h7NtPSTv+i7pQNU1ixbwUP9XqIBsENaiBCUdude+65XHPNNaUXXS+//HKOHDlCp06dePTRRzn77LPp2rUr9957LwBXXXUV48ePp1u3bmzduvW4wy1PmDCBiRMn0qVLFzIzM838E/9GhlYWfs3j0Qx/+0+27stl9n0DiA45tb71Ow/v5NLpl9KzcU8mDpooZ/k+UF+HVjabDK0s6hSLRfHspV04Uuji2ZnrT+k1Hu3h8QWPY7PYeLzv45LwhShDkr7we20bhXHzWS35Ii2DhdsOnHT/zzd+ztK9S3mg5wMymJoQx5CkL2qFu1LbkBAdxJivV5/wYSu7cnfxUtpL9Gvaj0taX+LDCIWoHSTpi1ohKMDKs5d0Ydu+PF7+edNx93tusdFVcFzfcdKsI0QFJOmLWuPMNg24ulcC787bVuGAbHN3zmXOzjnc2vVWmoY2NSFCIfyfJH1Rq4w5vwNNIoK4f9rKcs08Ba4Cnlv8HK0iWnFdB/+6A1IIfyJJX9QqYQ47z11mNPO8VKaZ591V75KZm8mjfR7FbpVHH9ZXVquV5ORkOnfuzEUXXcShQ4fMDqlSBgwYQE11V5ekL2odo5mnGe/O20baXwfZnrOdD9Z+wEUtL6Jn455mhydMFBQUxIoVK1izZg3R0dF+NdhZRcNAmEGSvqiVHr2gA00jgvjX58t58s+nCbIGcW/KvWaHJfxI3759y90Ne6pDKAOkp6eTmppKUlISgwYNYseOHeTk5NC8eXM8Hg8AeXl5JCQk4HQ6jzvE8qhRo7j11lvp3bs3Dz74IHl5edxwww306tWLbt26ld7BW1BQwFVXXUWHDh245JJLTvvhLqdDBlwTtVJooI0JVyVz9eS3OLh3MY/2fpTYoFizwxJezy9+nrzl2LIAACAASURBVA3Z1fvYy/bR7Xmo10OntK/b7Wb27NnceOONAKc1hDLAnXfeyciRIxk5ciTvv/8+d911F9988w3Jycn89ttvDBw4kBkzZjBkyBDsdjujR4+ucIhlgIyMDBYsWIDVamXMmDGkpqby/vvvc+jQIXr16sXgwYN5++23CQ4OZv369axatYru3btXa92VJUlf1FrtmtqJSviB3Lw4bHl9zQ5H+IGCggKSk5PJzMykQ4cOnHPOOUD5IZTBGBZ58+bNrFy5ssIhlP/880+++uorwBga+cEHHwTgyiuv5LPPPmPgwIFMnTqV22+/vdwQyyWKiopKl6+44gqsVmtpHNOnT+eFF14AoLCwkB07dvD7779z113GM5uTkpJISkqqsTqSpC9qrTdWvEGhJ4dWltsZ9+16ejSPpUVsiNlhCTjlM/LqVtKmn5+fz5AhQ5g4cSJ33XVXlYZQLmvYsGGMGTOG7Oxs0tLSSE1NJS8v74RDLIeEHP1Oaq358ssvadeu3en/cdWkym36SimrUmq5UmqGd72FUmqRUmqLUuozpVSAtzzQu77Fuz2xqp8t6q8N2RuYvGEyw9sN5+3h/4fNauGuKcspdnnMDk34geDgYF599VVefPFFXC7XaQ+h3K9fv3JDI5955pmA8YDznj17cvfdd3PhhRditVpPOMTysYYMGcJrr71WOlT48uXLATjrrLOYPHkyAGvWrGHVqlU1US1A9VzIvRsoOxLW88DLWuvWwEHgRm/5jcBBb/nL3v2EOG0e7eGphU8RGRjJnd3upGlkEM9flsTqzBye/r52P7RaVJ9u3bqRlJTElClTTnsI5ddee40PPviApKQkPv74YyZMmFD6vldeeSWffPIJV155ZWnZ8YZYPtbYsWNxOp0kJSXRqVMnxo4dC8Btt91Gbm4uHTp04PHHH6dHjx41VzFa60pPQDwwG0gFZgAK2A/YvNv7Aj95l38C+nqXbd791Inev0ePHlqIY3256UvdeVJn/e2Wb8uVPz1jrW7+0Az9+ZIdJkVWv61bt87sEOqliuodWKqPk1ereqb/CvAgUPKbOgY4pLUu6ZCaAcR5l+OAnd4DjQvI8e5fjlJqtFJqqVJq6b59+6oYnqhrDhYe5KW0l+jesDsXtbyo3LaHhranf+sYHv1mDSt31s6bcoSoaZVO+kqpC4EsrXVaNcaD1vodrXWK1jqlQQN52pEob8KyCeQV5/FYn8f+NqCazWrhtau70yA0kFs/SWN/btFx3kWI+qsqZ/r9gWFKqXRgKkYTzwQgUilV0isoHii5OyITSADwbo8ATj44uhBeK7JW8OXmL7mu43W0iWpT4T7RIQG8fV0PsvOKuf3TZXJh18e0Hz+Jry6qTH1XOulrrR/RWsdrrROBq4BftdYjgDnA5d7dRgIlVzSme9fxbv9VyzdEnCKXx8XTC5+mUXAjbu166wn37RwXwX8vT2Lx9mwe/nKVJCIfcTgcHDhwQOrbR7TWHDhwAIfDcVqvq4l++g8BU5VSTwPLgf95y/8HfKyU2gJkYxwohDglUzdMZePBjbw04CWC7cEn3f/i5Dh2HMjnxZ830STSwQND2vsgyvotPj6ejIwM5Fqc7zgcDuLj40/rNdWS9LXWc4G53uVtQK8K9ikErji2XIiTycrP4vUVr9M/rj+Dmw0+5dfdkdqaXTmFTJyzlSYRQVzbp3kNRinsdjstWrQwOwxxEnJHrvB7Lyx5AafbyZheY07raVhKKZ66uBNZhwt5/Ns1NAp3cE7HRjUYqRD+T0bZFH7tz11/8kP6D9zU5SaahTc77dfbrBZeu6YbXeIi+OfkZczbLE0Pon6TpC/8VrG7mGcXPUtCWAI3dLmh0u8THGBj0vW9aNUglJs+XMqCLfurMUohahdJ+sJvTVo7ifTD6YzpPYZAa2CV3isqJIBPb+pNi9gQbvhwCX9uld7Con6SpC/8UsaRDN5Z9Q7nND+HM+LOqJb3jA4J4JObepMQFcwNkyTxi/pJkr7wO1pr/rP4P1iUhQd7Plit7x0bGsjkm/sQFxXEyA8W89PaPdX6/kL4O0n6wu/M2TmH3zN+55/J/6RxSONqf/8GYYF8fktfOjYJ57ZP0pi8aEe1f4YQ/kqSvvAr+c58nlv8HK0jW3NNh2tq7HOiQwKYfHNvzm7bgDFfr+bV2ZvlTlJRL0jSF37lnVXvsDtvN2P7jMVusdfoZwUH2HjnHylc2j2Ol37exH3TVlLodNfoZwphNrk5S/iNbYe28eG6D7m41cV0b1RzD4Yuy2618OIVXWkeHcLLv2xiS1Yub13bg6aRQT75fCF8Tc70hV/QWvPMomcItgVzb8q9Pv1spRR3D27Du/9IYdu+PIa9/geLt2f7NAYhfEWSvvALM7fPZPGexdzd/W6iHdGmxHBOx0Z8889+hDvsXP3uQl7/dTNuj7Tzi7pFkr4w3eHiw4xfMp4usV24rM1lpsbSumEY39zRn/O7NOGFWZu4+p2FZB4qMDUmIaqTJH1huglpEzhYdJDH+jyG1WI1OxzCHXZevSqZl4Z3Zd3uwwx95Xe+XZEpvXtEnSBJX5hqRdYKPt/0OSM6jKBjTEezwymllOLS7vHMvOtM2jQM5e6pKxj1wRJ2ZuebHZoQVSJJX5jG6XHyxJ9P0DikMXck32F2OBVqFhPMtFv7Me6ijixNz+acl3/j7d+24nLLYxhF7SRJX5jmo7UfseXQFsb0GnNKT8Myi9WiuL5/C36+92zOaN2A//ywgaET5jF7/V5p8hG1jiR9YYqdR3by1sq3GNRsEAObDTQ7nFPSNDKId//Rg3eu64Hbo7nxw6WMeG8RazJzzA5NiFMmSV/4nNaaZxY+g0VZeLjXw2aHc1qUUpzbqTGz/nUWTwzrxPrdh7nwtT+47ZM01u6S5C/8n9yRK3zup/SfmL9rPg/3erhGBlTzBbvVwsh+iVzSPY53f9/GpPnp/LBmD4M7NOSO1DYkJ0SaHaIQFVL+3CaZkpKily5danYYohodLj7MsK+H0SikEZPPn+wXXTSrQ06Bk0nz03l//nZyCpz0aB7FqH6JDO3cGLtVflAL31JKpWmtUyraJmf6wqdeTnuZg0UHeWPwG3Um4QNEBNm5e3Abbjgjkc+XZvDhgnTunLKcxuEOru7VjMt6xBEf5b8Xq0X9IWf6wmcW7l7IzbNuZlSnUdyXcp/Z4dQot0czd2MWkxakM2/zfpSCfq1iuKJHAud2akRwgJxviZpzojN9SfrCJ/Kd+Vw6/VJsFhtfXPQFDpvD7JB8Zmd2Pl8uy+CLtAwyDhYQZLeS2r4h53dpwsD2DeQAIKqdNO8I001YNoFdubuYNHRSvUr4AAnRwdwzuC13pbZh0fZsvl+9ix/X7OH71btx2C30bRnDgHYNObttAxJjQ8wOV9RxkvRFjUvbm8bkDZMZ0WGEz8bJ90cWi6Jvqxj6torhiWGdWbw9m5/W7mHuxizmbFwLQGJMMGe3bcAZbRqQ0jyKqJAAk6MWdY0074gaVegq5LLpl+HWbr4a9pVf33lrpvT9efy2aR+/bdrHgq37KXQawzy0aRhKSmI0Kc2j6JkYTUJ0EEopk6MV/k6ad4RpXln2CjuO7OC9c9+ThH8CibEhJMaGMLJfIoVON6sycliSns3S9GxmrNrFlMXGw9sjg+10ahpO56YRdGwaTqemEbSIDcFqkQOBODWS9EWNWbBrAZ+u/5QRHUbQu0lvs8OpNRx2K71aRNOrhfEwGY9HszkrlyXp2azdlcOazMN8MD+dYu+gbw67hVYNQmnVIJTWDY/Om8cE47DXnW6xonpI0hc1Iqcoh7F/jKVlREvu6X6P2eHUahaLol3jMNo1Distc7o9bMnKZe2uw6zffZgtWbks23GQ6St3lXttw7BA4qOCSIgONuZRwSREBxMXGUSjcAdBAXJQqG8qnfSVUgnAR0AjQAPvaK0nKKWigc+ARCAdGK61PqiMhsgJwPlAPjBKa72sauELf6S15umFT5NdmM1rg16rd711fMFutdChSTgdmoSXKy8odrN1Xy5b9+Wy40A+Ow/mszO7gGU7DjJj1e6/Pf4xLNBGw/BAGoY5vHNjOSY0gKjgACKD7UQFG8thDhsWaUaq9apypu8C7tNaL1NKhQFpSqmfgVHAbK31c0qph4GHgYeA84A23qk38KZ3LuqYmdtn8mP6j9zV7S6/ejBKfRAUYKVzXASd4yL+ts3l9rDncCE7swvIPFRA1pFCsg4Xlc6X7zhE1pHC0ovIx7IoiCxzIIgMshMSaCPUYSM00EZIQMmyldBAOyGBVkK920MCbDjsVhx2Cw67VYamMFGlk77Wejew27t8RCm1HogDLgYGeHf7EJiLkfQvBj7SRnehhUqpSKVUE+/7iDoiMzeTZxY+Q3KDZK7vfL3Z4YgybFYL8VHBJxwOQmvN4UIX2XnFHMwv5lB+MQfznN7l8vM9hwvJK3KR652Od7CoiNWicNgs3gOBlUC7hUCb96BgK39wMCaFzarKr1uMud1qwWa1EGBV2KwWbBZFgM2CzWLBZlVYlcJqUVgsCosCqzKWrd51S8l2pUqXrRbKrVssxvtYFEeXve+hAKVAoSjpWFV23dheZj+Te19VS5u+UioR6AYsAhqVSeR7MJp/wDgg7CzzsgxvWbmkr5QaDYwGaNasWXWEJ3zE6XbywG8PAPDsmc9is8glo9pGKUVEkJ2IIDstOL0bxVxuD3lFbnKLXeQWHj0YlBwYipxuCp0eCp1uCl1llp0eCl3uctuPFDkpdHpwuT043Rqn24PLo3G6PDg9Hlxujcvjv93NT0WFBwSMwmMPJGX3o+x6JQ4qVf5fqZQKBb4E7tFaHy77gVprrZQ6rX8ZrfU7wDtg9NOvanzCd15e9jKr96/m5QEvkxCWYHY4wsdsVgsRwRYigu0++TyPx0j8TrdxEHB6PKXLxSVlbg8erXF7NB7N0WXvulsby8Z27d1+tPzoa43yo/scLdfauKhpzI11MH41VbRNGxsrLC9Zp3T9JO9/nPc40d1NVUr6Sik7RsL/VGv9lbd4b0mzjVKqCZDlLc8EymaCeG+ZqAN+3fErH6/7mGvaX8Pg5oPNDkfUAxaLIsDblCPKe+4E2ypdW97eOP8D1mutXyqzaTow0rs8Evi2TPk/lKEPkCPt+XXDrtxdPDb/MTrGdKzzo2cKUdtV5Uy/P3AdsFoptcJbNgbjIPO5UupG4C9guHfbTIzumlswumzKVb46oMhdxP2/3Y/WmhfOeoEAq4wVI4Q/q0rvnT8wriFUZFAF+2vgn5X9POF/tNY8+eeTrN6/mlcGvEJCuLTjC+HvpDFMVNrH6z5m+tbp3N71dgY1/9txXgjhhyTpi0pZsGsBL6a9yKBmg7il6y1mhyOEOEXSkVqcth2Hd/DAbw/QMqIlz57xLBZVyXMHreFwJmQshV3L4OBfcGQ3HN4NeVmgS272UWALhLDGENYEwptCdCto0hWaJkNow2r724So6yTpi9NyoOAAt/1yG0opXk199fSHSy46Apt/hg3fQ/ofkLvHKLcGQGQzI6k372skcouN0k7JzgLjgHBkN2yfByunHH3PsKbQvB+0HgStUo2DgxCiQpL0xSnLd+Zzx+w72Ju/l/fOfe/Ub8ByFcOG72DlVNg2F9zFEBxrJOj4nhDXAxp3Ns7mT1XhYdizGnavgMxlsP03WPOFsa1RF+g4DDpdCrGtT/vvFKIukydniVPi9Di569e7WLBrAa8MeIWBzQae/EXZ2yFtEiz/BPL3Q0Qz6HARdLgQEnqDpRqH9fV4YO9q2PorbPoJdvxplDdOgi6XQ9drILRB9X2eEH7sRE/OkqQvTkprzWPzH2P61umM6zuOy9tefuIX7F4J816EddNBWaDdeZByPbRMBYuP+g7kZMK6b2DNV5C5FCx2aH8B9BgFLc72XRxCmEAelygqTWvN80ueL+2aecKEv2MR/D4etvwMgeFw5r3Q8ybjwquvRcRB338a075NsOxDWPGpcSCIagE9RkLytXL2L+odOdMXx1WS8D9d/ynXdbyOB1IeqHgEvwNb4efHYcMMo62+7+1Gsnf8fUx3UzkLYf13RpPTX3+ANRCShkOf26GRjPsv6g450xen7ZQSfn42/PZfWPKe0ftm4GNGwg84vSF5fcbugKQrjGnfRlj0FqyYAss/hpYDoM8/ofVgafoRdZqc6Yu/8WgPzy9+nskbJlec8F1FsPhd+P2/RhfMbtfBwEchrNHx39Rf5WdD2gfG33NkN8S2hd63QterIeA0u6MK4SfkQq44ZcXuYsbOH8vM7TP/nvC1hnXfwi/j4GA6tBoE5z4FjTqZGnO1cDth7TewcCLsWg5BUdDjeuh1sznXJISoAkn64pQcKT7Cv+b8i0V7FnF397u5sfONRxN+xlL46VHYuRAadjSSfes6OG6+1rBjoZH8N3xv9D7qdInR7h/X3ezohDgl0qYvTmpv3l5un3072w5t49kznuWiVhcZGw7tgF+eMG58CmkIF00wer1Y6+hXRynjjuDmfY1fM4vehmUfw+pp0KyvkfzbX1C99xgI4UNypi9YkbWC++beR54rj5cGvES/pv2gMAfmvQQL3zTOdvvdAf3vhsAws8P1vcLDxsXeRW8ZB8HIZka7f7frwBFudnRC/I0074gKaa2ZsmEK45eMp0loE14e8DLtIloaXRrn/gfyDxgXNFPHGv3e6zuP22jyWfiGccdvQBh0vw563wJRiWZHJ0Qpad4Rf5PvzOfJhU/y/bbvGRA/gGfOeJrwbb/DpyPgwGZIPBPOfdoYxVIYLFZjTJ+Ow4zxfha+AYvfMX4BtDnXuPDb5hxp+hF+Tc7066EVWSt4bP5j7Di8gzu73cmNYe2wzH4SMpZAbDsY/G9j6ISKbsQS5eVkwtL3YdlHxnDQEQnQfaTxC0BG+xQmkeYdARjdMSeumMiktZNoHNyYpzveSM/l02DzLGN44oGPGAOT1dWLtDXJ7TSaftI+MEYStdig7VBIvsb4FWC1mx2hqEck6QuWZy3nyT+fZMuhLVwWn8oD2YcIWfutcSHyjHuNdml7kNlh1g0HthrJf+VUyNsHwTHQ5Qrj+kiTrvILStQ4Sfr12N68vby87GW+3/Y9DQOjGOcO56zN8yAg1LjxqP/dxo1Iovq5ncZQzysmw8aZxnMEGnY0hnru+H8Q08rsCEUdJUm/Hsp35jN5w2TeWfUOLreTUTqMm9JXERwYAX1ug16jITja7DDrj4KDxjDPK6dCxmKjrFEX6HQxdLxEHvYiqpUk/XqkwFXA5xs/5/0175NdmM0At40Hd+0gITAK+t4BPW+sn33t/UlOhvGsgXXfwM5FRlnDTtBuqNH+H99TegCJKpGkXw/kFOXw9eavmbTmfxwoOkSfIje3H9hHt/BW0Hs0JF0lA4j5o5xMWD8d1s8w+v5rt9Hc1nowtBliPFIyJMbsKEUtI0m/Dtt6aCuT10ziu23fU6Cd9C4o5LZDR+iROMhowmlxllw4rC0KDhnXADbPMh4en7/fKG/U2fh3TDzTeAB8UKS5cQq/J0m/jjlcfJhZW6bz3brJLMvbSYDWXJCbxzWWaNonXQddr5I+4rWdx22M9rltDmyfZzQDuQqNITGadIWEPhCfAgm9jHsD5MAuypCkXwccKT7CH5u+5pdNXzP3yFaK0bQodjLMaeGyxPOI6nKVkQTkP3/d5Cw0nvW7fR6k/wGZaeAqMLaFNjKuA8T1MB4E37izUSbfhXpLhmGohbTWbN+1mD83fslvexaxpDgbl4Jot5vLXTaGNe5Px64jUQm95UlP9YHdAYlnGBMY3UH3rjXuos5YAjsXG4+rLBEcayT/xl2MXkKNO0NMa7AFmhO/8Btypu8nnMV5bN0+m9U7fyMtawWLi7LY583lzV1uUm2xDEw4i6Qu12Ft0M7cYIV/KjhoHAj2rIY9a2DvasjaAO4iY7uyGE1BsW2MA0BMa+NegZjWEB4vJw91iJzp+xOtKTqYTnrmArbsXsaa7LWsyd/DBuWi0GL8HI92e+hli6R3g670bn0R8S0Go2xyG784iaCo8r8GANwuYwC9PWuM+YEtxrRjIRTnHt3PGgDhcRARbwwdHRFvHCBK5mGNITDU93+TqHaS9GtC0RGKs9PZs28Nmdkb2X34L9KP7GRb0QG26SIyrRY83vZWh9Z0sAZxRUgLOjfsSpfmg4iP74eySj9tUQ2sNmjYwZjK0hpy9x49CGRvM+4fyMmArXOM5wVzTCuAPRhCGxoP0wltaFw3CG14tCwk1jjwOCKNHkbSlOSXfJ70lVJDgQmAFXhPa/2cr2M4LR43FB2GwsN4Cg+Rl7ePg4d3cuBIBgfy9pBdsJ8DRQc5UHyEbHcBWZ4idllgv9WKLnMhza6huTWADoHxXBCWQMuYjrSM603Lpr2wWwNM/ANFvaSUcfYe1rj8L4MSbicczjQOAod2GgeIvH3GPHevMb7QXwugIPv4n2EPNg4CZQ8EjkjjF0NA6NF56XKI8YyC0uVQYzwoa6A0PVUjnyZ9pZQVmAicA2QAS5RS07XW6yp8gfYYTy3yuIzJ7Ty67J20qxi324nbXYSrZHIV4/YU43IX4XQVUlicR5EzjyJnPoWufIpchRS6CihyF1LoKqLYXUShx0mhu4hcdxFHPMXkaie52s0Rpcm1WMi1WMhTqlwiLysSC9F2Ow1sEfR3xNA0NI6mES1oGtOeJg060TgsDptFfliJWsJqNx4Mc7KHw7idRw8G+QeMew0KDhrzwjLLBQeNA0VhjtGsVJxr/P8+5XgCweYwfj3YHUeXbUHeuaNMuXey2o07my1277Ld+OVjsRujoJYsl2yzWP++n9W7r7KUnyzWMutW4yBaYXnJ/pbjlJcs+66nla+zUC9gi9Z6G4BSaipwMVBh0t9yYB0XfNILFwqXApdSuCk/d1VjZQVoCLVaCbXZCLUEEWZ10NwaRKg9mLCAUEIDwgkNjCAqpDEx4c2IiUgkOrQRUY4o7BZpcxf1kNUO4U2N6XRoDc4CI/kXHYHiPO/BIM+77l12FoCryLhHoXQq+nt5/n6jW2vJdleBcT3D4wKP92TR3ykLoLwHAHX0QHJsWbk5x3nN8fk66ccBO8usZwC9y+6glBoNjAaIaRZC5+gO2JQVm8WGtWRusWFXNqxWGzaLHavFjs1iLNusJevGss3qwBEQQmBAGI6AMAIDw3HYQwi0OXDYHARaA0snq4x3IoRvKGUMCxLgvU5Q07Qu01rgNJptS5bLtiCUlpU5WLhdxq+S0sl9dNnjNt67wnLP36dy5WVeW1KONsrKzT3eZSrYdpw5rxy3KvyuvUFr/Q7wDhhdNp+/9GuTIxJC1HpKGb9K6s3DbI6f9H19dSQTSCizHu8tE0II4QO+TvpLgDZKqRZKqQDgKmC6j2MQQoh6y6fNO1prl1LqDuAnjC6b72ut1/oyBiGEqM983qavtZ4JzPT15wohhPB9844QQggTSdIXQoh6RJK+EELUI5L0hRCiHvHr8fSVUkeAjWbHUYFYYL/ZQRzDH2MC/4xLYjp1/hiXP8YE/hVXc611g4o2+N0ducfYeLwHAZhJKbXU3+Lyx5jAP+OSmE6dP8bljzGB/8Z1LGneEUKIekSSvhBC1CP+nvTfMTuA4/DHuPwxJvDPuCSmU+ePcfljTOC/cZXj1xdyhRBCVC9/P9MXQghRjSTpCyFEPeK3SV8pNVQptVEptUUp9bBJMSQopeYopdYppdYqpe72lkcrpX5WSm32zqNMiM2qlFqulJrhXW+hlFrkra/PvENX+zqmSKXUF0qpDUqp9UqpvmbXlVLqX95/uzVKqSlKKYcZdaWUel8plaWUWlOmrMK6UYZXvfGtUkp192FM473/fquUUl8rpSLLbHvEG9NGpdSQmojpeHGV2XafUkorpWK966bVlbf8Tm99rVVK/bdMuU/qqlK01n43YQy7vBVoCQQAK4GOJsTRBOjuXQ4DNgEdgf8CD3vLHwaeNyG2e4HJwAzv+ufAVd7lt4DbTIjpQ+Am73IAEGlmXWE8nnM7EFSmjkaZUVfAWUB3YE2ZsgrrBjgf+AHjCah9gEU+jOlcwOZdfr5MTB29/w8DgRbe/59WX8XlLU/AGJb9LyDWD+pqIPALEOhdb+jruqrU32J2AMep4L7AT2XWHwEe8YO4vgXOwbhLuIm3rAnGTWS+jCMemA2kAjO8X/j9Zf6zlqs/H8UU4U2w6phy0+qKo89kjsa4EXEGMMSsugISj0kaFdYN8DZwdUX71XRMx2y7BPjUu1zu/6A3+fb1VV15y74AugLpZZK+aXWFcfIwuIL9fFpXpzv5a/NORQ9QjzMpFgCUUolAN2AR0Ehrvdu7aQ/QyMfhvAI8CHi86zHAIa21y7tuRn21APYBH3ibnd5TSoVgYl1prTOBF4AdwG4gB0jD/Loqcby68Zfv/w0YZ9FgckxKqYuBTK31ymM2mRlXW+BMb1Phb0qpnn4Q00n5a9L3K0qpUOBL4B6t9eGy27RxKPdZv1el1IVAltY6zVefeYpsGD9/39RadwPyMJosSplQV1HAxRgHpKZACDDUV59/OnxdNyejlHoUcAGf+kEswcAY4HGzYzmGDeNXZB/gAeBzpZQyN6ST89ek7zcPUFdK2TES/qda66+8xXuVUk2825sAWT4MqT8wTCmVDkzFaOKZAEQqpUrGUjKjvjKADK31Iu/6FxgHATPrajCwXWu9T2vtBL7CqD+z66rE8erG1O+/UmoUcCEwwnswMjumVhgH7pXe7308sEwp1djkuDKAr7RhMcYv71iTYzopf036fvEAde9R+3/Aeq31S2U2TQdGepdHYrT1+4TW+hGtdbzWOhGjXn7VWo8A5gCXmxGTN649wE6lVDtv0SBgHSbWFUazTh+lVLD337IkJlPrqozj1c104B/enil9gJwyzUA1Sik1O9F9bgAAAQ9JREFUFKPpcJjWOv+YWK9SSgUqpVoAbYDFvohJa71aa91Qa53o/d5nYHSw2IOJdQV8g3ExF6VUW4zOC/sxsa5OidkXFU5w0eR8jN4yW4FHTYrhDIyf3KuAFd7pfIw29NnAZoyr99EmxTeAo713WmJ8sbYA0/D2KPBxPMnAUm99fQNEmV1XwBPABmAN8DFGjwqf1xUwBeO6ghMjad14vLrBuDA/0fvdXw2k+DCmLRjt0SXf97fK7P+oN6aNwHm+rKtjtqdz9EKumXUVAHzi/W4tA1J9XVeVmWQYBiGEqEf8tXlHCCFEDZCkL4QQ9YgkfSGEqEck6QshRD0iSV8IIeoRSfpCCFGPSNIXQoh65P8B1s5L61nOI5oAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"for day in range(simulation_time):\n",
" france_simulation.time_step()\n",
" historical_data.loc[day] = france_simulation.states\n",
"\n",
"historical_data.plot.line()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment