{ "cells": [ { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import datetime as dt\n", "import requests as rq\n", "\n", "from scipy import integrate\n", "from scipy import optimize" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Define SEIR" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [], "source": [ "def diff_eqns(y, t, beta, sigma, gamma, N):\n", " St, Et, It, Rt = y\n", " dSdt = -beta * St * It / N\n", " dEdt = beta * St * It / N - sigma * Et\n", " dIdt = sigma * Et - gamma * It\n", " dRdt = gamma * It\n", " return ([dSdt, dEdt, dIdt, dRdt])\n", "\n", "def seir_model(t, beta, sigma, gamma, E0, N):\n", " S0 = N - E0\n", " I0 = R0 = 0\n", " out = integrate.odeint(diff_eqns, (S0, E0, I0, R0), t, args=(beta, sigma, gamma, N))\n", " return out.T # S, E, I, R" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sample SEIR plot" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEWCAYAAACKSkfIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3xUVdrA8d+Zkkx6L4RQEnpv0sSCCgpIc0XFfW1r233dteur67qru66u7rq4oq6KK4IV1gYIgoIQEQSkikjvJEAq6cnU8/4xE0ggCekzyTxfP+PM3Hvuuc9ckmduzj33HKW1RgghhH8weDsAIYQQLUeSvhBC+BFJ+kII4Uck6QshhB+RpC+EEH5Ekr4QQvgRSfpCtAJKKa2U6urtOETrJ0lftApKqYuUUt8rpQqUUnlKqbVKqaGedbcppZxKqeKzHkme9YeVUmOqKVuolPpRKTXRm59NiJYkSV/4PKVUOLAYeAWIBtoDfwaslYqt01qHnvU4XkOV67TWoUAk8G9gnlIqshk/ghA+Q5K+aA26A2itP9JaO7XWZVrrr7XW2xtTqdbaBbwHhADdqiujlIpVSi1WSuV7/sL4Till8Kx7XCl1QClVpJTaqZS6ptJ2t3n+GnnJs+1BpdSFnuXHlFJZSqlbK5Wfo5R6Qym13FPft0qpTjXEFKiUelEpdVQplenZLqgxx0L4D0n6ojXYCziVUnOVUuOVUlFNUalSygj8CrADR2oo9jCQDsQBCcATQMXYJQeAi4EI3H95vK+Ualdp2+HAdiAG+BCYBwwFugI3Aa8qpUIrlf8f4BkgFtgGfFBDTC/g/iIc6KmrPfCnunxmISTpC5+ntS4ELsKdbN8CspVSi5RSCZWKjfCcUVc8DtRS5QilVD5QDrwI3KS1zqqhrB1oB3TSWtu11t9pz4BVWuuPtdbHtdYurfV8YB8wrNK2h7TW72itncB8oAPwF621VWv9NWDDnbQrLNFar9ZaW4E/ACOVUh0qB6OUUsBdwINa6zytdRHwHDC9ls8rxGmS9EWroLXepbW+TWudDPQFkoB/VSqyXmsdWenRpZbq1mutI4EoYBHus/Wa/APYD3ztaaJ5vGKFUuoWpdS2ii8aT1yxlbbNrPS6zPM5zl5W+Uz/WKXPWwzkeT5nZXFAMLC50n6XeZYLcV6S9EWro7XeDczBnWQbU08xcA9ws1JqUA1lirTWD2utU4FJwENKqSs87e1vAb8DYjxfIjsA1YiQTp/Ve5p9ooGzL0bn4P6y6FPpCy7Cc2FaiPOSpC98nlKqp1LqYaVUsud9B+BGYH1j69Za5wL/oYY2caXURKVUV0+zSiHg9DxCcDc3ZXvK/YpGfgkBEzxdUwNwt+1v0Fofq1zAc/H5LeAlpVS8Z9/tlVJXNXLfwk9I0hetQRHui6IblFIluJP9DtwXWSuMrKaf/tA61v8v3Am3fzXrugErgGJgHfBvrXWa1non8E/PskygH7C2IR+ukg+Bp3A36wzBfWG3Oo/hbnJar5Qq9MTXo5H7Fn5CySQqQnifUmoOkK61ftLbsYi2Tc70hRDCj0jSF0IIPyLNO0II4UfkTF8IIfyIydsB1CYyMlJ37SqjyVYoKSkhJCTE22H4DDkeZ8ixqMrfj8fmzZtztNbV3rDn00k/ISGBTZs2eTsMn5GWlsbo0aO9HYbPkONxhhyLqvz9eCilahpLSpp3hBDCn0jSF0IIPyJJXwgh/IgkfSGE8COS9IUQwo+cN+krpWZ7pnbbUWlZtGdat32e5yjPcqWUmqmU2q+U2q6UGlxpm1s95fdVniZOCCFEy6nLmf4cYNxZyx4HvtFadwO+8bwHGI97VMJuwN3A6+D+ksA9euBw3DMLPdVUU94JIYSou/P209dar1ZKdT5r8RRgtOf1XCAN93CvU4B3PdPJrVdKRXrmDB0NLNda5wEopZbj/iL5qLZ951s1Ly3fi0EpDAqMRoXZYMBoUJhNBgKNBgLNBgJNBixmI0FmI0EBRoIDTIQGmgi1mAg2GzEYGjOvhRBCtB0NvTkrQWt9AkBrfaJiMgfcEzRXnvQh3bOspuXnUErdjfuvBAISu/LyN/saGKKnPiDYDKFmRYhZER6gCA9URAQqogIVMUGKWIuBmCCFxeTbXw7FxcWkpaV5OwyfIcfjDDkWVcnxqFlT35FbXdbUtSw/d6HWs4BZAD169NC7/zYBrcGpNU6XxuHSOJwubE4XNocLq8OF1e6i3OGk3Oak1OakxOag2OqgxOqgsMxBQZmd/DI7+aU2copt7CqwkltsxXVWBInhFrrEh9A1LpRe7cLpnxxJ94RQTEbfuN7t73cZnk2OxxlyLKqS41Gzhib9TKVUO89Zfjsgy7M8nUrzfALJuOf4TOdMc1DF8rS67EgphVJgQGE2NjDaajhdmuwiKxn5ZWTkl3Esr5QD2cUcyC7h0y0ZFFvddzFbzAb6J0dyUddYRnWNZUByhM98CQhRVy6bDVdJCbq0FFdZGa6ycnR5Ga7yclzl5WirDW0tR9tsaJsNl9WGttvQdvvpBw6H57UD7XSiHXZwON2vnQ5wutzPDifa5QKn59nlQruc4NLgcoF2oSteu1xoXek12l3Os0yj3aeHFWW0Zx2es8aK9xWjBXtexzmd7DYaqyw7/bqys+urvLwmDV3XkHJ1VY/6Gpr0FwG3As97nhdWWv47pdQ83BdtCzxfDF8Bz1W6eHsl8PsG7rtJGA2KxAgLiREWhnSqek3Z5dIcyStle3o+Px4rYOPhPF5asZcZy/cSZjFxVZ9EpgxMYmRqjHwBiBblKinBkZODIzsbR14ezlP5OE/lEfrTDo4v+wpnYSGuwkKcxcW4iopwFRfjLC0Fu71hOzQaUSYTymx2P0wmMJtQJjPKaESZjGAyowwGMBlRBqNnuQkVGAAGIxgNKGVw12VQoAygFMpoqPoahfsMzwAGhXta4opl7ucqy8D9rJRnUUWDguJYRjodkjucVa7i5VkNDxV1nn7Nua/PVltLcG3bVSnW1M3Jleq7/74aS5036SulPsJ9lh6rlErH3QvneeC/Sqk7gKPAdZ7iXwITcM/fWQr8CkBrnaeUegbY6Cn3l4qLur7IYFCkxIaQEhvClIHuSw95JTa+P5DDqt3ZfLXjJJ9sTic2NJBrh7TnVxemkBhh8XLUorXTWuPIysJ+9Ci29Azsx45hP34ce+ZJHCczsWdmoktLq902KDCQkuhojOHhGMPDMSclYQwLxRASiiEkxP0IDsYQEowhKAhlCcIQZEEFWjBYAlEWCyogEENgACqg0sPYhH9et6BdaWkk+HPzTmOSvtb6xhpWXVFNWQ38toZ6ZgOzz7c/XxUdEsDE/klM7J9Eub0vaXuy+GxLBm+tPsjsNYeYPKA9v740le4JYd4OVfg4rTWOzEzKd+/Gumcv1n37sB08iO3wYVyVk7pSmBISMCcmEtijB6GXXIwpLg5jbCymmFhMsTEYo6MxRkayet06acMWdeLTQyv7KovZyLi+7RjXtx1Hc0uZvfYQ8zce47Ot6Uwf2pFHr+pBdEiAt8MUPsJVUkLp1m2U/biNsu3bKd/+E85Tp06vNyW1IzC1CxFDhhCQ0pmAjp0I6JCMuV07VID8HImmJUm/kTrGBPP05D7cf0U3Xl21nznfH2bJ9uM8fGUP/md4R2nz90Pabqd08xZK1q6l9IcfKNuxA5xOUIqALqmEXnYZlj69sfTsSWD37hjD5K9D0XIk6TeRqJAA/jixN9OHduDpL37mqUU/8+VPJ3jlxkHEh0t7f1vnLC6heNVKilaupOS7NbiKi8FkIqhfP2LuvJPgoUMJGjgAY2iot0MVfk6SfhPrlhDG+3cM59MtGfxxwQ4mzFzDzBsHcmGXWG+HJpqYy2ajeFUahUuWUPztt2irFVNcHGHjriJs9GiCR4zEGOq/U/YJ3yRJvxkopZg2JJn+yRHc88EWbvrPBh69qie/uTS1GbppiZZmO3aM/P/+l/xPP8OZl4cxJobIadMIv/pqggYOcHdfFMJHSdJvRt0Twlj421E89ul2Xli2m6JyO49e1UMSfytVumUrubNmUZyWBkYjoZeNJur66wm58EJ3/3UhWgH5SW1mIYEmZk4fRJjFzL/TDmB1uHjy6l6S+FsJrTUla78n9403KN20CWNUFLG//S2R103DnJjo7fCEqDdJ+i3AYFA8d01fAk0G3l5zCKvDyTNT+kri93HlO3eS+fd/ULp+PaaEBBKe+D2R06ZhCA72dmhCNJgk/RailOKpSb0JNBl4c/VBEsIs3HtFN2+HJaphz8wie8YMChYtwhgRQcITTxA1/QbpMy/aBEn6LUgpxePje5JVZOWfy/fSIzGMK/tIE4Gv0FqT//HHZP3jRbTVSsyddxBz110Yw8O9HZoQTUa6GbQwpRR/+0U/+idH8OD8bezNLPJ2SAKwHTnC0Vtv4+SfnsLSqxepixYS//DDkvBFmyNJ3wssZiNv3jyEoAATd727ifxSm7dD8msFixZx8JpfUL5rF4l/+TMd584hoHNnb4clRLOQpO8l7SKCePPmwRzPL+P3n/3k7XD8kqusjONPPsnx/3uMoN69SV38BVHXXy8X2EWbJknfi4Z0iubBsd1ZuuMkS3864e1w/Irt6FEOX38DBZ9+Rsxvfk3HOe9gTkjwdlhCNDtJ+l5218Wp9EkK548Lf6agtIETXYh6Kdu2jcM3TMeRlUWHWbOIf+ABublK+A1J+l5mNhp44dr+nCq18dclO70dTptX+NXXHLn1NgxhYXSa9xGhF1/k7ZCEaFGS9H1A3/YR/PqSVD7enM53+7K9HU6blffe+2Q88ACWXr3oPO8jAlNSvB2SEC1Okr6PuO+KbqTGhvDE5z9hdTi9HU6bk/fuu2Q++yyhV1xOxznvYIqO9nZIQniFJH0fYTEbeXpyH47llTHvh2PeDqdNyXv3XTKf+xthY8eS/NJLGCwyv4HwX5L0fcjF3WIZkRrNKyv3U2pzeDucNqFywm8/458os9nbIQnhVZL0fYhSikev6kFOsZU53x/2djitXsEXX0jCF+IskvR9zJBO0VzRM5430g5QUCZdOBuqZMMPHH/iDwQPG0b7f74oCV8ID0n6PujhK3tQWO7grdUHvR1Kq2Q9cID0e+8loFNHkl+ZKaNjClGJJH0f1DspnEkDkpi99hDZRVZvh9OqOHJyOHbX3ajAADq++SbGiAhvhySET5Gk76MeHNONMruT99Yd9nYorYZ2OMh46GEceXl0eP0NzO3bezskIXyOJH0flRoXyhU94/lgw1Hpt19H2a+8SukPP9Duz08T1LePt8MRwidJ0vdht12YQm6JjcU/ymBs51OUlkbum28Sed11REyZ4u1whPBZkvR92KiuMXSND2XO94fRWns7HJ9lz8jg+GOPE9irFwlP/sHb4Qjh0yTp+zClFLdd2JmfMgrYcvSUt8PxSdrpJOPhR8DpJPnlf2EIDPR2SEL4NEn6Pu4Xg9sTZjHxztrD3g7FJ+XNmUPZtm0kPvUUAR07ejscIXxeo5K+UupBpdTPSqkdSqmPlFIWpVSKUmqDUmqfUmq+UirAUzbQ836/Z33npvgAbV1wgInpQzuwdMdJ8spd3g7HpxiPnyD75ZmEjR1L+MSrvR2OEK1Cg5O+Uqo9cB9wgda6L2AEpgMvAC9prbsBp4A7PJvcAZzSWncFXvKUE3Vw84jOuLRm1VEZj6eCdjiImDsXQ0gIiU8/JVMcClFHjW3eMQFBSikTEAycAC4HPvGsnwtM9bye4nmPZ/0VSn5T66RjTDCX9YhnTYYDp0su6ALkvj0b85EjJD71J0wxMd4OR4hWo8FJX2udAbwIHMWd7AuAzUC+1rrilDQdqLhDpj1wzLOtw1Neflvr6BeD23PKqll/MNfboXid9eAhcl59lfIhQwgfN87b4QjRqjR4YlClVBTus/cUIB/4GBhfTdGKU9PqzurPOW1VSt0N3A0QFxdHWlpaQ0NsUwKcGotR8/rSzdj7+XEPFa2JnPkKZrOJkxMnys+HR3FxsRyLSuR41Kwxs0GPAQ5prbMBlFKfARcCkUopk+dsPhk47imfDnQA0j3NQRFA3tmVaq1nAbMAevTooUePHt2IENuWoTu/YmsODL/wYoICjN4OxysKv/6ajF27SHjySbLbJSI/H25paWlyLCqR41GzxrTpHwVGKKWCPW3zVwA7gVXANE+ZW4GFnteLPO/xrF+p5Y6jerkwyUSx1cGKXZneDsUrXGVlZD7/PIE9ehA1/QZvhyNEq9SYNv0NuC/IbgF+8tQ1C3gMeEgptR93m/3bnk3eBmI8yx8CHm9E3H6pR7SBdhEWPt+a4e1QvCJn1iwcx0+Q+McnUabG/JEqhP9q1G+O1vop4KmzFh8EhlVTthy4rjH783cGpZg8MIn/fHeI3GIrMaH+07ZvO3KEvP+8TfikSQRfcIG3wxGi1ZI7cluZXwxKxunSLN7uX4OwZf1zBspsJv6RR7wdihCtmiT9VqZHYhi92oX7VRNP2U8/UfT110TffjvmhHhvhyNEqyZJvxW6ZlAS247lczS31NuhtIisGTMwRkURfdtt3g5FiFZPkn4rNK5POwC+3nnSy5E0v5Lvv6d03Xpi//c3GENDvB2OEK2eJP1WqGNMMD0Tw1i+s2133dRakzXjJUxJ7YicPt3b4QjRJkjSb6XG9k5g4+E88kps3g6l2RR9vZzyHTuIu/c+DAEB3g5HiDZBkn4rNbZ3Ai4NK3dneTuUZqGdTrJffpmArl2ImDzJ2+EI0WZI0m+l+rWPIDHcwvI22q5ftHw5toMHifvd71BG/xxyQojmIEm/lVJKMbZ3Aqv35lBud3o7nCaltSbnjTcJSEkhbOxYb4cjRJsiSb8VG9s7gTK7kzX7crwdSpMqWb0a6+7dxNx1l5zlC9HEJOm3YiNSYwgLNLWpXjwVZ/mmpHZETJro7XCEaHMk6bdiASYDo3vGs2JXZpuZUat040bKtm4l5o47UGazt8MRos2RpN/Kje2dQG6Jja1HT3k7lCaR++YsjLGxRF57rbdDEaJNkqTfyo3uEYfJoFixq/V33Szb8TMla9cSc9utGCwWb4cjRJskSb+VC7eYGdwpiu/2ZXs7lEbLmzsXQ0iI3H0rRDOSpN8GXNItlp+PF5JTbPV2KA1mz8qicNkyIq79BcbQUG+HI0SbJUm/DbikexwAa/e33q6b+fPmgcNB9E03eTsUIdo0SfptQJ+kCKKCzXy7t3U28bisVk7Nm0/oZZcR0LGjt8MRok2TpN8GGA2KUV1j+W5fDq1xrvnCxUtw5uURfcvN3g5FiDZPkn4bcUn3OLKLrOzJLPJ2KPWitSbvvfcI7NaN4OHDvR2OEG2eJP024uJusQB8t7d1teuXbtyIdfduom65GaWUt8MRos2TpN9GtIsIolt8KKtbWdfNU++9jzEykohJMnyyEC1Bkn4bckn3ODYcyms1o27as7IoWrmSiGt/ITdjCdFCTN4OQDSdi7vF8vaaQ/xwKO90N05fVvDZ5+B0EnXddd4ORYgq7HY76enplJeXezuUWlksFpKTkzHXY5wqSfptyPCUGAJMBlbvzfb5pK9dLvI//pjgESMI6NzZ2+EIUUV6ejphYWF07tzZZ681aa3Jzc0lPT2dlJSUOm8nzTttSFCAkWGdo/muFYyvX7L2e+wZGURdL2f5wveUl5cTExPjswkf3BMpxcTE1PuvEUn6bcyFXWPYk1lEro8PyZD/3/kYo6MJGzPG26EIUS1fTvgVGhKjJP02ZkRqDAAbDuV5OZKauS/griLimqmogABvhyOEz3r22Wfp06cP/fv3Z+DAgWzYsKHRdUqbfhvTr30EIQFG1h3IZUK/dt4Op1pyAVeI81u3bh2LFy9my5YtBAYGkpOTg81ma3S9kvTbGLPRwAWdo1l/MNfboVRLLuAKUTcnTpwgNjaWwMBAAGJjY5ukXkn6bdDILjE8v3Q32UVW4sICvR1OFaUbNmDPyCDuwQe9HYoQdfLnL35m5/HCJq2zd1I4T03qU2uZK6+8kr/85S90796dMWPGcMMNN3DppZc2et+NatNXSkUqpT5RSu1WSu1SSo1USkUrpZYrpfZ5nqM8ZZVSaqZSar9SartSanCjoxfVOtOu73tn+wULFmAICyNszBXeDkUInxYaGsrmzZuZNWsWcXFx3HDDDcyZM6fR9Tb2TP9lYJnWeppSKgAIBp4AvtFaP6+Uehx4HHgMGA908zyGA697nkUT65sUTmigiXUHcpnYP8nb4ZzmLC6h8OvlREyaJHfgilbjfGfkzcloNDJ69GhGjx5Nv379mDt3Lrfddluj6mzwmb5SKhy4BHgbQGtt01rnA1OAuZ5ic4GpntdTgHe123ogUinlm1caWzmT0cDQzlE+165f9NVX6LIyIq6Zev7CQvi5PXv2sG/fvtPvt23bRqdOnRpdb2PO9FOBbOAdpdQAYDNwP5CgtT4BoLU+oZSK95RvDxyrtH26Z9mJypUqpe4G7gaIi4sjLS2tESG2LcXFxXU+HvHYWZVtY8GylURafKNnbtScORji4/khPx+a4N+1PsejrZNjUVVjj0dERARFRd4dpjwzM5NHH32UgoICTCYTqampzJw585y4ysvL6/dZtdYNegAXAA5guOf9y8AzQP5Z5U55npcAF1Va/g0wpLZ9dO/eXYszVq1aVeeyPx47pTs9tlgv3JbRfAHVg/XoUb2zR0+d/fobTVZnfY5HWyfHoqrGHo+dO3c2TSAtoLpYgU26hrzamFPAdCBda11xt8AnwGAgs6LZxvOcVal8h0rbJwPHG7F/UYs+SRGEedr1fUHBgoWgFBFTJns7FCH8WoOTvtb6JHBMKdXDs+gKYCewCLjVs+xWYKHn9SLgFk8vnhFAgfY0A4mmZzQohqVEs8EH2vW1y0XBwoWEjByBuZ1cxhHCmxrbe+de4ANPz52DwK9wf5H8Vyl1B3AUqLjt8ktgArAfKPWUFc1oZJcYvtmdxcmCchIjvNdbpmzzZuzp6cTdd6/XYhBCuDUq6Wutt+Fu2z/bOZ2wPe1Mv23M/kT9VO6vP2Vge6/FUbBoEYbgYBlcTQgf4BvdOkSz6NXO3V9/42HvDb7mstko/OprwsaOwRAc7LU4hBBukvTbMKNBMbhTFD94ccTNktWrcRUWEj5xotdiEEKcIUm/jRueEs3ezGJOlTR+dL6GKFi8BGN0NCEjR3pl/0K0VkajkYEDB55+PP/8801Srwy41sYN7RwNwMbDeVzZJ7FF9+0sLqZ41Soip01DmeRHTYj6CAoKYtu2bU1er5zpt3H9kyMIMBm80q5ftHwF2molfOLVLb5vIUT15PSrjbOYjQxMjvRKu37h4sWYk5MJGjiwxfctRJNZ+jic/Klp60zsB+Nrb64pKytjYKXfnd///vfccMMNjd61JH0/MCwlmte/PUCJ1UFIYMv8kzuysylZt46Yu+9qFXONCuFrmqt5R5K+HxiaEs2rq/az5egpLu4W1yL7LFy6DFwuIqTXjmjtznNG3tpIm74fGNIpCoOCjS3YxFOwZDGBvXoR2LVri+1TCHF+kvT9QGigiT5JEWxooaRvS0+n/MfthE8Y3yL7E6ItqmjTr3g8/vjjTVKvNO/4iWEp0by//ghWh5NAk7FZ91W4dCkA4eMl6QvRUE6ns1nqlTN9PzG0czRWh4uf0guafV+FS5di6d+fgOTkZt+XEKJ+JOn7iaGdowD4oZn769sOH8a6c5ec5QvhoyTp+4mY0EC6xoc2e3/90007465q1v0IIRpGkr4fGdo5ms2HT+F06WbbR+GXSwkaPFgmSxHCR0nS9yPDUqIosjrYfbKwWeq37t+Pdd8+adoRwodJ0vcjpwdfa6YmnsIvl4JShF11ZbPUL4RoPEn6fiQ5KpikCAsbD59q8rq11hQuXUrw0KGY4+ObvH4h/E1oaGiz1CtJ388MS4lmw6E83LNXNh3r3r3YDh2SG7KE8HGS9P3M0JRocoqtHM4tbdJ6C5cuBYOBsLFjm7ReIUTTkjty/cywSu36KbEhTVKn1pqiZV8RPHwYppiYJqlTCF/xwg8vsDtvd5PW2TO6J48Ne6xJ66wrOdP3M13jQ4kKNjfpTVrWPXuwHT5M+FXjmqxOIUTzkDN9P6OU4oLO0U06k1bhsmXupp0rpWlHtD3eOiNvLnKm74eGp0RzJLeUrMLyRtdVpWknOroJohNCNCdJ+n6oor9+UzTxSNOOEK2LJH0/1CcpnOAAY5OMwyNNO0I0j+Li4mapV5K+HzIZDQzuGNXopK+1pmjpMmnaEaIVkaTvp4alRLMns4iCUnuD67Du2YPtyBHCx8kNWUK0FpL0/dTwlGi0bly7fuGyZWA0EjZ2TBNGJoRoTpL0/dSADpEEmAxsOJjboO1PN+0MGypNO0K0IpL0/ZTFbGRwx0jWH2pY0pemHSFap0YnfaWUUSm1VSm12PM+RSm1QSm1Tyk1XykV4Fke6Hm/37O+c2P3LRpneEoMPx8vpKCs/u36hUulaUeI1qgpzvTvB3ZVev8C8JLWuhtwCrjDs/wO4JTWuivwkqec8KIRqTFoDZvq2a6vtaZw2VJChg+Xph0hmonRaGTgwIH07duXSZMmkZ+f3yT1NirpK6WSgauB/3jeK+By4BNPkbnAVM/rKZ73eNZf4SkvvGRQR3e7/vp6tutbd+/GfuQoYTIPrhDNJigoiG3btrFjxw6io6N57bXXmqTexo698y/g/4Awz/sYIF9r7fC8Twfae163B44BaK0dSqkCT/mcyhUqpe4G7gaIi4sjLS2tkSG2HcXFxU1+PFLCYMWPRxgVklXnbUIXLCDYYGB7cDDai/8+zXE8Wis5FlU19nhERERQVFTUdAE1UEUMgwYNYseOHdXGVF5eXq/P2uCkr5SaCGRprTcrpUZXLK6mqK7DujMLtJ4FzALo0aOHHj169NlF/FZaWhpNfTy22Pfy6sp9DB4xinCL+bzltdYc+NvfCBgxgt4TJzZpLPXVHMejtZJjUVVjj8euXbsIC3Ofy5587jmsu5p2aOXAXj1JfOKJ85YLCwvD6XSydu1a7rjjjtMxVWaxWBg0aFCd992Y5p1RwHcUehoAACAASURBVGSl1GFgHu5mnX8BkUqpii+TZOC453U60AHAsz4CaJ7JWkWdjUiJxqVhcx2nULTu2uVu2hkvY+0I0ZzKysoYOHAgMTEx5OXlMbaJJihq8Jm+1vr3wO8BPGf6j2it/0cp9TEwDfcXwa3AQs8mizzv13nWr9RNPWefqLdBHaMIMLrb9S/ref65bQuXfeXutTNGeu0I/1CXM/LmUNGmX1BQwMSJE3nttde47777Gl1vc/TTfwx4SCm1H3eb/due5W8DMZ7lDwGPN8O+RT0FBRgZ0CGC9XUYh8fda2cZISNGYIqKaoHohBARERHMnDmTF198Ebu94cOmVGiSpK+1TtNaT/S8Pqi1Hqa17qq1vk5rbfUsL/e87+pZf7Ap9i0ab0RqDDsyCigqr/0HqvznndiPHiVcmnaEaFGDBg1iwIABzJs3r9F1yR25guEpMThdmk1Ham/XL/zySzCbpWlHiBZw9tDKX3zxBTfffHOj65WkLxjSyd2u//3+nBrLaJeLwmVLCR01CmNkZAtGJ4RoSpL0BUEBRi7oHMV3+2pO+mXbfsRx/AThE2SsHSFaM0n6AoCLusWy+2QRWUXVz5tb+OWXqMBAQi+/vIUjE0I0JUn6AoCLu8YBsLaaJh7tdFL41TJCL70UY2hoS4cmhFe0hh7lDYlRkr4A3PPmRgWbq23iKd24CWd2DuETJnghMiFansViITc316cTv9aa3NxcLBZLvbZr7Ng7oo0wGBQXdo1lzb4ctNZUHguv8MsvUcHBhF56iRcjFKLlJCcnk56eTnZ2trdDqZXFYiE5Oble20jSF6dd3DWWJdtPsC+rmO4J7jE+tN1O0VdfEXb55RiCgrwcoRAtw2w2k5KS4u0wmoU074jTLuoWC1Cliadk3TqcBQXStCNEGyFJX5yWHBVMSmwIa/ad+ZO2YNEXGCIiCLlolBcjE0I0FUn6ooqLu8Wy4VAeNocLZ3EJRStWED5+HIaAAG+HJoRoApL0RRUXdY2l1OZky9FTFH+zAl1eTsTkyd4OSwjRRCTpiypGdInBaFCs2ZdDwaIvMLdvT1A9JmgQQvg2SfqiinCLmYEdItmydR8l69YRPnkSMpWxEG2HJH1xjimdnVy4eR64XET0j4Uj68Dl9HZYQogmIP30xRlZu+C7Gdy84xMOpUejoiHw23vgWyA6FUbcAwN/CQEh3o5UCNFAcqYvoOwU/PcW+PcI2L0Ea6ebsJ4KYH2/8XDHCrj2bQiKhi8fgRm9YdNsb0cshGggSfr+rigT5kyE3V/CJY/CgzsozE3BZTAwI+wSSuIHQb9pcNc3cMdyaDcAFj8IXzwADpu3oxdC1JMkfX+WfxTeGQd5B+GX8+HyJ9GBERQsWoRz8DCyTSFVB2DrMAxu/hwuehA2vwPvToFi3x6bRAhRlSR9f5V7AGaPg9JcuHkBdL0CgJI1a3BkZtLhpulEBJlZsSuz6nYGI4x52t3kc3wLvD0GCk+0ePhCiIaRpO+PrMXw0Y3gKIfblkDH4adX5X/yCcboaCIvv4zLesSxcncWTlc1w8v2mwa3LoaSHHj/F1Ca14IfQAjRUJL0/Y3WsPgByN0H02ZDYr/Tqxw5ORStSiNiyhRUQABjeieQV2Jjy9EaJkzvMBSmfwC5++HDG8BW0kIfQgjRUJL0/c2m2fDTxzD6CUgdXWVVwcJF4HAQOe1aAC7tHofZqFixM/PceiqkjnY39WRsgvk3y8VdIXycJH1/cnwrLHscuo6Bix+uskprTf6nnxI0aBCBXboAEGYxMyI1huVnt+ufrfdkmPQyHPjGXb8QwmdJ0vcXtlL4+FcQEg/XzAJD1X/6sq1bsR08ePosv8LY3gkczC5hX2ZR7fUPvgVG3Q+b3oYt7zZ19EKIJiJJ31+k/Q1OHYJr3oCQmHNW53/yKYbgYMLHjauyfFzfRAwKFm47fv59XPEUpF4GSx6G9E1NFbkQoglJ0vcHx7fButfcZ+MpF5+z2llcTOGyZYRfPQFDSNUhFuLDLFzULY4F2zLOP0m0wei+OBzWDubf5L7xSwjhUyTpt3VOByy6F4JjYOxfqi1S8PkCdGkpkddfX+36awYlkX6qjM1HaujFU1lwNEz/EMoL4OPb3PsXQvgMSfpt3fp/w8ntMOHvEBR1zmrtcnHq/fcJGjiQoH79qqkAruydSJDZyOdbM+q2z8S+7gu7R7+HVc82JnohRBOTpN+WnToMq56D7uOh99Rqi5SsWYPtyBGibrqpxmpCAk1c2SeBxdtPYHO46rbv/tfD4FthzQzY+3UDghdCNIcGJ32lVAel1Cql1C6l1M9Kqfs9y6OVUsuVUvs8z1Ge5UopNVMptV8ptV0pNbipPoSowbInQCm4+kX3czXy3nsfU1wc4VeOrbWqqYPaU1BmJ21PVt33P/4FSOgHn98NBen1iVwI0Uwac6bvAB7WWvcCRgC/VUr1Bh4HvtFadwO+8bwHGA908zzuBl5vxL7F+exbAXuWuEfOjEiutoj10CFKvvuOyBuno84z8fnFXWOJCQlgwbY6NvEAmIPgujngtHva9+1131YI0SwanPS11ie01ls8r4uAXUB7YAow11NsLlDRrjAFeFe7rQcilVLtGhy5qJnDBsseg+guMPK3NRY79cGHKLOZqBou4FZmMhqYNCCJFbuyKCirR/KO7QqTZ0L6RljxdN23E0I0iyaZOUsp1RkYBGwAErTWJ8D9xaCUivcUaw8cq7RZumdZlSEalVJ34/5LgLi4ONLS0poixDahuLi4Tsejw9HP6JK7n+39/kTemnXVllFlZcR+/DHWIYNZs2NHnfbfUTuxOVz865M0Lu1grkfkMXRLGk/7da/yU1E4ubHDz79JHdT1ePgDORZVyfGohda6UQ8gFNgM/MLzPv+s9ac8z0uAiyot/wYYUlvd3bt31+KMVatWnb9QQYbWzyZp/cENtRbLnTNH7+zRU5du/6nO+3e5XHrsjDQ94eXV2uVy1Xk7rbXW9nKt37hY67910DrvcP22rUGdjoefkGNRlb8fD2CTriGvNqr3jlLKDHwKfKC1/syzOLOi2cbzXHHlLx3oUGnzZKAOt3mKeln+J3fb+bjnaizistnInf0OwRdcQFC/vnWuWinFbRem8PPxQjYerkOf/cpMge72fa3d7fsyMJsQXtGY3jsKeBvYpbWeUWnVIuBWz+tbgYWVlt/i6cUzAijQnmYg0UQOr3WPoHnhve6JzGtQ8PkCHJmZxN7zv/XexTWD2hMZbGb2mkP1jy86Faa86p585asn6r+9EKLRGnOmPwq4GbhcKbXN85gAPA+MVUrtA8Z63gN8CRwE9gNvAfc0Yt/ibE4HLP0/CE8+ZwTNyrTdTu6sWVgG9Cd45Mh67yYowMgvh3Xk650nOZZXWv84e0+Bkb+DjW/Bj/Prv70QolEafCFXa70GqL7zN1xRTXkN1NyVRDTO5ncgcwdcNxcCgmssVrB4CfaMDBKe/AOqhr7753PzyE7MWn2Qd9cd5g9X965/BWOedg/z/MX9kNC7ykQuQojmJXfktgUlObDyGUi51H0mXQPtdJL75psE9upF6OjRDd5du4ggxvdrx7yNxyi2NmBsHaPZ3b4fFOkemK2sntcHhBANJkm/LfjmL+6pCsf/vcY7bwEKly3Ddvgwsb/5TYPP8ivcPqozReUOPt3cwDttQ+Ph+nehIAM+uxtczkbFI4SoG0n6rd2xH9yTlgz7NcT3rLGYdjrJfeMNArp0IWzsmEbvdlDHKAZ1jGT22kPYnXUcj+dsHYa5h2rY97W715EQotlJ0m/NnHZ3u3h4Elz2+1qL5n/6KdZ9+4m7916UoXH/7HaXnYziDMYNKSO9fAszVi/nSOERimznmV2rOkPvcH9hrXsVNs89f3khRKM0yR25wku+fwWydsL0jyAwrMZizuJisl+eSdCQIYRddWW9d+N0OdmatZWVx1ayOn01RwuPonFPqBLcEd4/6n4AdArvxAUJFzAkYQgXJl1ITNC5s3Sd46rnIHc/LHkIolMg5ZJ6xyiEqBtJ+q1V3kH49gXoNQl6Tqi1aO6st3Dm5pLwxuv1assvtBXy/s73mb9nPnnleZgNZoa3G874lPG0C2lHYkgi6bkufr9gA+MGhDEoxcCPWT/y9ZGv+XTfp5iUiUuSL+GabtdwUfuLMBlq+HEzmuC6d+A/Y2H+zXDHcojrXp+jIYSoI0n6rZHWsPghMJjdF29rYc/IIG/OHMInT6pxkpSzVST793e+T5G9iNHJo7m6y9Vc3P5iQsxVp1MkCVb9GMSqTdk8fdll3NnvTpwuJ3tP7WXpoaUsPLCQlcdWEh8cz219bmNa92kEmYLO3aklAn45H96+Et6bCnd8XePooEKIhpM2/dbox3lwcBVc8Sd3e34tsma8BEoR/+CD561Wa81Xh79i0ueTeP3H1xnWbhgfT/qYV654hXGdx52b8D0evaoHVoeLmd/sA8BoMNIrphcPXfAQK65bwcuXvUyn8E78fePfGffpOGbvmE2pvZobu6JT4KZPwVoE710DpXnnPxZCiHqRpN/a5B9133nbcaT7ImgtSjdvpnDJEqJv/xXmdrWPYp1dms2DaQ/yyLePkBiSyPyJ8/nXZf+iZ3TNPYIqpMaFcuOwDnz0w1H2Z1W9mGs2mLm84+XMvmo2c8fNpWd0T17a/BKTPp/E5/s+x3l2V812/eHGee7P+cE0sBafd/9CiLqTpN+aaCd8/ht38841b4DBWGNRV1kZx594AnP79sTeeWet1a46uoqpC6fyXfp3PDjkQT6Y8AG9Y+p3p+39V3Qn1GLi4Y+346ihC+fghMG8OfZN3h3/Lomhifzp+z9xw+IbWH9ifdWCnUfBtHfg+Db4aLokfiGakCT9VqTDsYVwZK27b3tU51rLZv/rX9iPHKXds3/FEFJ9s4zdZeefm/7Jfavuo31oez6Z/Am397295guutYgLC+SvU/vy47F8Xlt1oNayg+IH8f749/nHpf+g2F7MXV/fxUNpD3GiuNL4ez0nwC9mwZHvPWf8DegOKoQ4hyT91uLEdlIOfeDurTPwl7UWLd28mbx33yPyxumEjBhRbZnMkkzu+OoO5vw8hxt63MB7E94jJSKlUSFO7J/ElIFJzFy5j+3p+bWWVUoxrvM4Fk5dyL2D7uW79O+YsnAKb21/C5vTM+xyv2kw7W33rFvvXQPlBY2KTwghSb91KC+ET+/Abg6DiS/XOtTC6WadpCQSHnmk2jLbsrYxfcl0duft5oWLX+DJEU8SaAxsklD/MrkvcaGBPDh/G+X28w+tEGgM5O7+d7Nw6kJGJY1i5taZXLvoWr4//r27QJ9r3IPIHd8GcydDcT0mZhdCnEOSvq9zudxj0+QeYFevhyGk9pudsv7xj1qbdT7f9zm3f3U7QaYgPpzwIRNSa+/jX18RwWZevG4AB7JL+PMXOytmSTuvpNAkXrrsJd4Y8wYu7eLXy3/NI98+QmZJJvSaCNM/gOw98J8xkL23SWMWwp9I0vd1aX+DvUth3PPkR9Xezz7/s8859eFHRN966znNOg6Xg+d/eJ4/ff8nhiQM4aOrP6JrVNdmCfmibrH85tIufPTDUd749mC9th3VfhSfTfmMewbew6qjq5i8YDLv7HgHe5fL4VdLwF4Gb4+Bw2uaJXYh2jpJ+r5s5yJY/XcYeBMMu6vWomXbt3Py6acJHjGC+EerNusUWAv4zYrf8MGuD7ip1028PuZ1IgIjmjNy/u+qHkwakMQLy3bz2Zb6jcQZaAzkfwf8LwumLmBo4lBmbJ7BtV9cyzplgztXQGgivDuVdse/cvdkEkLUmSR9X5Wx2d09s/0FMHFGre34juxs0u+9D1NsLO1fmoEynel9s//Ufm5cciNbMrfwzKhneGzYYw3qnVNfBoPixev6MzI1hv/7ZDvf7cuudx0dwjrw6hWv8urlr2J32rl7+d3cv3UGR69/G1Iupsfef7uPka2kGT6BEG2TJH1fdHIHvPcLCImFG953TypeA1dZGen3P4CzoIDk117FFBV1et03R77hpqU3UeYoY/ZVs5nadWpLRH9aoMnIm7cMoWt8KL9+bzPf7q1/4ge4tMOlLJi6gPsH38+6E+uY8tUtzOg5ip86Xwfb58Nbl0PW7iaOXoi2SZK+r8neC+9OAXMw3LoIwmu+k9ZVXk76b39L2datJD3/Nyy9egHu9vsZm2fwQNoDpISn8NHVHzEwfmBLfYIqwi1m3r19GJ1jQrh9zkY+3HC0QfUEGgO5s9+dLLlmCRNTJzJn57vcafyJuZfdi7UkB968BNbOlMlYhDgPSfq+JO8gvDsZlAFu/aLWG7BcVivpv7uXknXraffcc4SPGwdAblkuv17+a97Z8Q7Xd7+euePnkhiS2EIfoHrx4Rb++5uRXNwtlic+/4kXlu3G5WpYW3xccBzPjHqG+RPn0ymgEy8eXsDEzp34PGUw9uV/hNlXSe8eIWohSd9XHNvoHlrYUQ63LITYmnvWuGw20u+7j5I1a2j312eIvMbdbLMmYw3XLrqWH7N/5NmLnuWPI/9IgDGgpT5BrUIDTfznlgv4n+EdeT3tALfM/oH0U9UMulZHvWJ6cU/CPbx15VvEBMfzJ2cGk3r057PSI9jfGAUrnpa7eIWohiR9X/DzApg7EQJD4favIaHmcW8cOTkcve1XlHy7msRn/kLktddS7ijn+R+e539X/C9Rlig+vPpDJneZ3IIfoG5MRgN/ndqXZ6/py9ajp7jqpdV8sOFInfvyV2dEuxF8dPVHvHL5K0SEJfFUhIWJnTrxwfa3KH1lCGx9332vgxACkKTvXS4XrH4RPr4V2g2AO7+pdfIQ05EjHJp2HeW7dtH+pRlEXXed++7axdNPd8ecN3Ee3aN8dwISpRT/M7wTyx64hAEdIvnD5zu4YdZ6Nh1u+DDKSilGdxjNvKvn8doVrxEf24vnY6IZE2vh5dVPkPnGCNjxmSR/IZCk7z35x9zt9yufgb7T4JZF7t461dBak//Z50S/+E8wKDp/+AH68gv587o/c/PSmylxlPDGmDd4bNhjTTacQnPrEB3MB3cO57lr+nEwu4Rpb6zjltk/8OOx2sfsqY1SikuSL+G9Ce/x3vj3GNHpCt6OjOCq4DLuX/0oa2cNw/XjfHDYmvCTCNG6yMxZLU1rdzfDLx8F7YLJr8Kgm2rsh28/cYITTz9NybersXfrRuo7/2Fx3mpeWfAbCqwF3Nr7Vu4ZeA/B5uAW/iCNp5Til8M7MnVQEu+uO8Kb3x5gymtrGdQxkhuHdWRi/3YEBzTsR3Rg/EAGxg/kWNExPtnzMQv2zGelo5SkjU9z9do/M6nbVFJG3A9hCU38qYTwbaox7anNrUePHnrPnj3eDqPpHN8KX/8RDn/nngTlmjdq7KGjHQ5OzZ9P9j9noLUm7oH7eTc6h5WONA4XHmZA3AD+MPwP9Irp1bKfoRkVlduZv/EY8zYeY39WMWGBJq7sk8jY3glc0j32nC+AtLQ0Ro8eXae6bU4b3xxezoId77A+fw8uoK/VxpXBHRnT+5d06P/LWu+H8HX1ORb+wN+Ph1Jqs9b6gmrXSdJvAXmHYNWz8NPHEBwLox+HC26vdhIU7XRS+OWX5Lz6GrYjRwgaOYKf77iY2XlL2HtqL10ju3LfoPsY3WF0vSY5b0201mw6cop5Pxxjxa5MCsrsBJgMjEyNYXhqNMNTounXPpLv16xu0C92dmk2X/78Hl/u/YydDvdwzT3sTkaHpnBR14n07XcTpsCwJv5Uzcvfk9zZ/P14SNL3Bq3dE56sfx32fAnGABj5Wxj1AFjCzynuKi2lcOlScme/g+3AAQzdUtk+pRevhW4k15pH18iujDCO4JEJj2CsZcastsbudLHp8CmW78xk9b5s9me5Z9EKNBloHwIjeibTu104PRLDSI0NITokoF5fhscLjvLNj2+x4sg3bHMW4lKKMJeLEcZIhsYP4oIe19Cl02gMPn7M/T3Jnc3fj4ck/ZZ06gj8/Bls/xiyfoagaBhym3vAtLMmMddaU75jB/mffkrhF4txlZRQ3imeZZeF81HiIbRSXNT+Im7ufTMj243k22+/9esfZIDcYisbD59i0+E81u48QkapgcJyx+n14RYTKXGhJEcG0T4qiKQIC4kRFuLCAokPcz9bzNUn8IKSHNbvmMvaQ8tZV5rBSU+xSJemvzmKflE96d/xUvp0GU9EcO1DXLc0f09yZ/P341Fb0pcLuY3ldMDxLXDwW9i7DDI2uZe3HwKTZkL/68EcdLq4q6yM0s1bKF75DUUrV+I4mYnTbOTnARF80rOc3cm5dI4I53dd7mNCygSSw5K99MF8U0xoIOP6JjKubyJpoVlceumlpJ8qY39WMQdzSjiUU8zhnFJ2nihk+a5MbI5zu2kGmY1EhwQQFWImIshMuMVMmMVEmMVMSOAkkpOmcqvZgLH8J/JylnG4+Cf2WPNYnbMectbDlhdIdBnoFhBF9/BUusT3pUv7C0lJHECQKaiaqIXwHZL060NrKDoJJ7a5L8pmbIFjG8BaCCh3X/sxT7tne4rqjHa5sKenU757N6XbtpH/wzpcu/ehHE7sZgPbUxQbhhjY1E2R2qEzlydfxJ+TL6VXdK82217f1JRSdIgOpkN0MJedtU5rTU6xjczCcrKLrWQXWskutpJfaiOvxM6pUhsFZXayCospKLNTbHVQaqs8do8RuNrz0CQZMkgN3kxo0EHsgbkcc2TyvT0HZ95G2P0OANEORYQ9AIsjAuWMxelKoFx3oJROmFQsRoMZk0FhMCiMCowGhUF5HgbOvFbu16rKa/czFe8Bg3Ifg8xMK4uytqFwl1PgeVanO4a5nyu9h0qvVaUyZ9ZXHONzj3st/yaoOparYXkT/OgfO2bj+9Jdja+oDWrxpK+UGge8jPs36j9a6+dbOoZa2UqhOBMKM9x96QvS4dQh96xNOfvA6pmnVRkgrie651TskYOwmztTmJ1DwfrdlM37M65jGQQcOYmp3A6AwwD728HuCxR7O5pxDOpF3/aDmRg/kD8nDiXaEu3FD902KaWICwskLqzuvXKcLk2pzUGJ1Umpzf0lUGZ3Um53Um53YXVMotzuwuZwYbU7UMWHsRduobh0FwW2o+TpPHIMpRwPKiXXlAXsrFJ/uBMiXSYsrkAsOgizDkG5wkGH4yQKmyuKchWFVUfh0OG4XBbQ7vMNl9ZoPM/a/aVW8b6szMmRUvcNbhXLtQb3q4rXlacf0KdfVyyq3NR7Zhnn3DFda4OwrvblucVqaFZuisZmrcHpcmLMONIEtbU9LZr0lVJG4DVgLJAObFRKLdJa76y2PBrs5Z6vfoX7p9/lfric4HKceXbazjw7rJ5HufthL0XbSrGXF+IqLcBRVoi9tAB7ST6OknwcpYXYSwuxFxfisFmxOxUOp8LuMOB0KJwE49JB4EhC2ZIxlLowlzgILComuGQFsKJK3MXhcDJKkdHXQFHHOJxdOhLaszfdEvowNao7qZGpreYmKn9jNCjCLGbCLOY6bpEKXH7uYnsZZbn7OJG1g5On9nGi8ChZpVnkWPPJchSTRzHZqoBcg4FyQ833SCqtCUERohUhykiwMhKkTAQbAggymrEYArAYA7AWlxMfHU+gMZAAo4UAUyABRvfDbLJgNgZgNgViNlowGQMwGwMwmdyvjYYAjEYzJkMARqPnYTBjMJo8z2aMBjPKYMJoMGEwmjEoIwZl8Nm/SP2+Tf+vNa9r6TP9YcB+rfVBAKXUPGAKZ58OeeiMo6RdemZIYKU9P2Da86eprljufhhcnuUuMGgwutzLjS73w3De0wgLRixUd5nPaiqjxFJGmUVRFGKmPDYAW2okzqhwdFw0xoR4Atu1JzKlG/FRHeganEB8cLxf9bQRlZiDCErsT2pif1JrKqM1WAspK8ygoDCd/OIT5JecpKgsjyJrPkW2IgrtxZQ6yih2llPitFHqslOmy8jSJZThohwoV+6HNS8D7YUkbNAaA+7fPQOe30fOPCq+0lQ1jyrLddUmH3XOszpnGdWUB0Br/n5QVbvu7G1VLe9q0tCj7AtfkS2d9NsDxyq9TweGVy6glLobuBsgNdxCYXxQpSNVqTHS3bDp/iFXoJUBlPu9NhjAYEAbDGijAZfRhDYY0UYj2hSANpnQRhPaaASz2f0wmSEgEGUOgIBADJZgDJZgjEGhmAJDCQgMIVAFYlJGwoFzO116nID8E/nkk88emrbnUXFxMWlpaU1aZ2vWto5HANAJ6ITZBNEmiD53XvsalRQVEBoUgHaV4XSW4XSV43JacehyXC4bTpcNl7bjcNlwaQdO7cCl7bhw4dIOXNqFUzvQuHBpJy7touI/rV3u5WhcuuK9u+nIvfTMfy5dsQbPEs+zdmf002srmqCASo1QZ5qVTn8yffr/1S3jrGUVr1wujcGgzimnqeXMr9amqcY3PPlKP8mWTvrVfdFVORZa61nALHB32Zy8YGtLxNUq+PufrGeT43FGWloal8qxOM3ffzZeu6vmvylaesC1dKBDpffJwPEWjkEIIfxWSyf9jUA3pVSKUioAmA4sauEYhBDCb7Vo847W2qGU+h3wFe4um7O11j+3ZAxCCOHPWryfvtb6S+DLlt6vEEIImURFCCH8iiR9IYTwI5L0hRDCj0jSF0IIP+LT4+krpYqgiW9rbd1igRxvB+FD5HicIceiKn8/Hp201nHVrfD1oZX31DQRgD9SSm2S43GGHI8z5FhUJcejZtK8I4QQfkSSvhBC+BFfT/qzvB2Aj5HjUZUcjzPkWFQlx6MGPn0hVwghRNPy9TN9IYQQTUiSvhBC+BGfTfpKqXFKqT1Kqf1Kqce9HU9LUkp1UEqtUkrtUkr9rJS637M8Wim1XCm1z/Mc5e1YW5JSyqiU2qqUWux5n6KU2uA5HvM9w3X7BaVUpFLqE6XUbs/PaYyuRwAABI5JREFUyUh//flQSj3o+T35//buLkSqOozj+PdHGvlCG0WGrcS2sGUWpbIXYlFh3lSySSQWGhJFN0UpRfRyEQVdBJFW9nKhmIgKoUt1VRcW4UUZ+EJvGwQWZqyukJqJpev+uvj/px1s3RG0OSP/5wMLc86cszxzeOaZM8858/9/J2mjpItKzo1GWrLo102gficwDXhA0rRqo2qqQeAp29cBs4DH8ut/FthiuwvYkpdL8iTQV7f8KrA8H4+DwMOVRFWNN4BPbE8FbiIdl+LyQ1I78ATQbfsG0pDt91N2boyqJYs+dROo2z4O1CZQL4Ltfts78uMjpDd0O+kYrM2brQXmVxNh80maAtwNrMrLAuYAm/ImxRwPSRcDtwKrAWwft32IcvNjDDBO0hhgPNBPoblxJlq16I80gXp7RbFUSlIHMAPYBlxhux/SBwMwqbrImm4F8AwwlJcvAw7ZHszLJeVIJ3AAWJPbXaskTaDA/LD9G/AasIdU7A8D2yk3Nxpq1aLfcAL1EkiaCGwGltr+o+p4qiJpHjBge3v96hE2LSVHxgAzgXdtzwCOUkArZyT5usU9wNXAlcAEUlv4VKXkRkOtWvSLn0Bd0lhSwV9vuzev3i9pcn5+MjBQVXxNdjPQI+kXUqtvDunM/5L8lR7KypG9wF7b2/LyJtKHQIn5MRf42fYB2yeAXmA25eZGQ61a9IueQD33q1cDfbZfr3vqY2BJfrwE+KjZsVXB9nO2p9juIOXCZ7YXAZ8D9+XNSjoe+4BfJV2bV90B/ECZ+bEHmCVpfH7f1I5FkblxJlr2F7mS7iKdzdUmUH+l4pCaRtItwFbgW4Z72M+T+vofAFeRkn2B7d8rCbIikm4HnrY9T1In6cz/UmAnsNj231XG1yySppMual8I7AYeIp3EFZcfkl4CFpLuetsJPELq4ReZG420bNEPIYRw7rVqeyeEEML/IIp+CCEUJIp+CCEUJIp+CCEUJIp+CCEUpNUnRg+haSSdJN0mO5Z0+99aYIXtoVF3DOE8EkU/hGHHbE8HkDQJ2AC0AS9WGlUI51C0d0IYge0B4FHgcSUdkrZK2pH/ZgNIWifp3xFgJa2X1CPpeklfS9ol6RtJXVW9lhDqxY+zQsgk/Wl74inrDgJTgSPAkO2/cgHfaLtb0m3AMtvzJbUBu4AuYDnwle31eSiRC2wfa+4rCuG/or0Twuhqo3mOBVbm4Q9OAtcA2P5C0tu5HXQvsNn2oKQvgRfyPAC9tn+qIvgQThXtnRBOI4/tc5I0WuUyYD9plqpu0pg3NeuARaTxb9YA2N4A9ADHgE8lzWle5CGcXhT9EEYg6XLgPWClUw+0DejPd/I8SBoIsOZ9YCmA7e/z/p3Abttvkka/vLF50YdwetHeCWHYOEm7GL5lcx1QG9r6HWCzpAWkYXuP1nayvV9SH/Bh3f9aCCyWdALYB7zchPhDaCgu5IZwliSNJ93fP9P24arjCWE00d4J4SxImgv8CLwVBT+cD+JMP4QQChJn+iGEUJAo+iGEUJAo+iGEUJAo+iGEUJAo+iGEUJB/AAhEXzad37y7AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "N = 1000\n", "days = np.arange(100)\n", "beta = 1.5\n", "sigma = 1 / 5.2 # 5.2 days incubation, i.e. from exposed to infectious\n", "gamma = 1 / 2.9 # 2.9 days from infectious to removal\n", "S, E, I, R = seir_model(days, beta, sigma, gamma, 1, N)\n", "\n", "df = pd.DataFrame({\n", " \"Days\": days,\n", " \"S\": S,\n", " \"E\": E,\n", " \"I\": I,\n", " \"R\": R\n", "})\n", "\n", "df.plot(\n", " x='Days',\n", " y=['S', 'E', 'I', 'R'],\n", " grid=True,\n", " title=\"SEIR sample\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fit Italy's data" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
datedeceasedrecoveredconfirmedactivedaily
592020-04-202411448877.01812281082372256
\n", "
" ], "text/plain": [ " date deceased recovered confirmed active daily\n", "59 2020-04-20 24114 48877.0 181228 108237 2256" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "it_df = pd.read_csv(\"csv/italy.csv\", parse_dates=['date'])\n", "it_df.sample()" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Params: [0.73755627]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAE5CAYAAABlF9zIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXwV1fn48c+TnQDZSMKSsIQtsq8CoiguILihVEVcQFwo1lZbtZWufq3a2uqvKm6UCgKtC6hVqIIISEA22dcEJEAggZBANhJC1nt+f8wEQ7hJINvNvXner97emTNnZp6bkfvcmXPmjBhjUEoppZzxcnUASimlGi9NEkoppSqlSUIppVSlNEkopZSqlCYJpZRSldIkoZRSqlKaJFSDEpEOIpInIt6XsM4IEdlfn3HVBftzdXZ1HErVJU0Sql6ISJKInLW/OMte7YwxR40xLYwxpXa9OBF5pMK6RkS6ls0bY74zxsQ29GeoirO47c91yFUxNSQRGSkiKbVY/0ERWVuXMan64ePqAJRHu9UYs8LVQSilak7PJFSDEpFO9pmCj4i8BIwA3rLPNN4SkTV21Z122YSKv1rts5RnRGSXiOSIyAIRCSi3/Dcikioix0XkkYpnJhXimSIiCSKSKyKHROSnFZaPE5EdInJaRA6KyBhncdt1jYh0FZFhInKi/CU1EblDRHbZ014iMt3eXoaILBSRsCr+ZhfEUF3sIhIuIl+KSLaIZIrIdyLiZS9rJyKfichJETksIk+UW2+IiGyx95UmIv9wEk9zYCnQrvxZooj4i8jr9t/9uD3tX9nnUm7CGKMvfdX5C0gCbnBS3gkwgI89Hwc8UqGOAbqWmx8JpFTY9iagHRAGJADT7GVjgBNALyAQ+HfF7VXY181AF0CAa4B8YKC9bAiQA4zC+kEVBVx2MXEDB4FR5ZZ9Aky3p38JbASiAX/gn8BHlcRXVQxVxf5XYCbga79G2PW8gK3AnwA/oDNwCLjRXm8D8IA93QIYVklc5x0Tu+zP9ueKBCKA9cALlaz/ILDW1f+d6qv6l55JqPr0hf1LNltEvqjjbc8wxhw3xmQC/wP62+V3A+8bY/YaY/KB56vaiDHmK2PMQWNZDXyD9YUK8DAwxxiz3BjjMMYcM8bsu8j4PgImAohIS+Amuwzgp8DvjTEpxphC4P+AO0XE2eXfSmOoJvZioC3Q0RhTbKx2HQNcDkQYY/5sjCkyVhvKv4B7yq3XVUTCjTF5xpiNF/l5Ae4D/myMSTfGnMT62z9wCeurRkiThKpPtxtjQuzX7XW87RPlpvOxfvWCdXaRXG5Z+ekLiMhYEdloX5LJxvoyD7cXt8c6I6iJD4Hx9uWW8cA2Y8wRe1lH4POyBIp1JlQKtHaynUpjqCb2V4BE4Bv7UtT0cvtuVy55ZwO/K7fvh4HuwD4R2Swit1zCZ24HHCk3f8QuU25MG66Vq9X1MMSpWJdxyrSvrKL9Bf4ZMAlYZIwpts94xK6SjHU5x5kq4zbGxIvIEWAscC9W0iiTDDxkjFlX1TaqiqG62I0xucDTwNMi0gtYJSKb7e0dNsZ0qyTuA8BEu/1iPPCpiLQyxpypWNXJ6sexktBee76DXabcmJ5JKFdLw7ouXl3ZxVoITBGRHiISiHXtvTJ+WG0CJ4ESERkLjC63fLa9revtxuYoEbnsEmL8EHgCuBqrTaLMTOAlEekIICIRIjKukm1UFkOVsYvILXYjugCnsc5USrHack6LyLMi0kxEvEWkt4hcbq93v4hEGGMcQLa9uVIncaUBrUQkuFzZR8Af7M8TjvW3/081fyPVyGmSUK72Btb1+CwRmWGX/R8wz74ccvelbMwYsxSYAazCutyywV5U6KRuLtaX+EIgC+sX/+JyyzcBU4DXsBqPV2P9Uq4s7oo+wmrg/dYYc6rCZ16MdSkoF6uxd2gln8dpDNXFDnQDVgB59t/gHWNMnLHuT7kVqw3nMHAKeA8o+7IfA+wVkTw7znuMMQVO4tpnf75D9nFqB7wIbAF2AbuBbXaZcmNitWUp5ZlEpAewB/A3xpS4Oh6l3I2eSSiPY9+T4CciocDfgP9pglCqZjRJKE/0U6xr9Qexrqc/5tpwlHJferlJKaVUpfRMQimlVKU0SSillKqUx91MFxISYrp2dTqWW4M6c+YMzZs3d3UYqo7o8XR/egyrtnXr1lPGmIiK5R6XJFq3bs2WLVtcHQZxcXGMHDnS1WGoOqLH0/3pMayaPULABfRyk1JKqUpVmyREZI6IpIvIHifLnhFrDP1we15EZIaIJIo11v/AcnUni8gB+zW5XPkgEdltrzPDHkYAEQkTkeV2/eV2n3ellFIN6GLOJOZi3ap/HhFpjzXG/dFyxWOxhgPoBkwF3rXrhgHPYQ09MAR4rtyX/rt23bL1yvY1HVhpD0S20p5XSinVgKptkzDGrBGRTk4WvQb8BlhUrmwcMN8et36jiISISFus8WuW22P/IyLLgTEiEgcEGWM22OXzgduxnno1zl4PYB7WQ16evaRPZysuLiYlJYWCgguGoKk3wcHBJCQkNNj+qhMQEEB0dDS+vr6uDkUp5UZq1HAtIrcBx4wxO+2rQ2WiOH/8/hS7rKryFCflAK2NMakAxphUEYmsSawAKSkptGzZkk6dOlEh3nqTm5tLy5YtG2Rf1THGkJGRQUpKCjExMa4ORynlRi45SdjDL/+e84dUPrfYSZmpQfmlxjQV65IVERERxMXFnbc8ODiYVq1akZeXd6mbrrHS0lJyc3MbbH/V8fPzIzs7+4K/jbo4eXl5+rdzc3oMa6YmZxJdgBisB9WD9YCXbSIyBOtMoPxDXqKxHjqSwo+XjsrK4+zyaCf1AdJEpK19FtEWSK8sIGPMLGAWQGxsrKnYzS0hIYGgoKBL+Yy11pjOJMoEBAQwYMAAV4fhlrT7pPvTY1iF+MWVLrrkLrDGmN3GmEhjTCdjTCesL/qBxpgTWOPZT7J7OQ0DcuxLRsuA0SISajdYjwaW2ctyRWSY3atpEj+2cSwGynpBTeb8tg+PFhcXx/r168/Nz5w5k/nz57swIqWUxyrKh4WTKl1c7ZmEiJQ9OCVcRFKA54wxsyupvgTrObuJWM8dngJgjMkUkReAzXa9P5c1YmON0DkXaIbVYL3ULn8ZWCgiD2P1oLqrulg9RVxcHC1atGD48OEATJs2zcURKaU81sl9VHWV/2J6N02sZnmnctMGeLySenOAOU7KtwC9nZRnANdXF587uf3220lOTqagoIAnn3ySqVOn8vXXX/O73/2O0tJSwsPDmT17NjNnzsTb25v//Oc/vPnmm6xcuZIWLVpw8803M3nyZDZt2gRAUlISt912G7t27WLr1q089dRT5OXlER4ezty5c2nbtq2LP7FSqtFL21vlYo8blqM6z/9vL/HHT9fpNnu2C+K5W3tVW2/OnDmEhYVx9uxZLr/8csaNG8ejjz7KmjVriImJITMzk7CwMKZNm0aLFi145plnAFi5ciUAPXr0oKioiEOHDtG5c2cWLFjA3XffTXFxMb/4xS9YtGgRERERLFiwgN///vfMmXNBTlZKqfOlx4NPM6xHoV+oySUJV5oxYwaff/45AMnJycyaNYurr776XLfUsLCwardx9913s3DhQqZPn86CBQtYsGAB+/fvZ8+ePYwaNQqwelbpWYRS6qKk7YXIy4A0p4ubXJK4mF/89SEuLo4VK1awYcMGAgMDGTlyJP369WP//v2XtJ0JEyZw1113MX78eESEbt26sXv3bnr16sWGDRvqKXqllMdK2wuxY4DVThfrAH8NJCcnh9DQUAIDA9m3bx8bN26ksLCQ1atXc/jwYQAyM622/JYtW1Z6j0WXLl3w9vbmhRdeYMKECQDExsZy8uTJc0miuLiYvXurvs6olFLkpUP+KYis/MezJokGMmbMGEpKSujbty9//OMfGTZsGBEREcyaNYvx48fTr1+/c1/6t956K59//jn9+/fnu+++u2BbEyZM4D//+Q933303YN0o9+mnn/Lss8/Sr18/+vfvf14XWqWUcirNHre1deVJosldbnIVf39/li5d6nTZ2LFjz5vv3r07u3btOjc/YsSI85Y/88wz5xq1y/Tv3581a9bUUbRKqSYhLd56ryJJ6JmEUko1VWl7oXkkNA+vtIomCaWUaqrS91Z5FgGaJJRSqmkqLYGT+zVJKKWUciLzEJQUaJJQSinlRLrdTT6yZ5XVNEkopVRTlLYXxAsiYquspkmigcyYMYMePXoQGhrKyy+/DMAXX3xBfHz8uTpz587l+PHjlW3CqaSkJHr3vmB8RKWUqlpaPLTqCr7NqqymSaKBvPPOOyxZsoSsrCymT58O1E2SUEqpGknfW+2lJtCb6RrEtGnTOHToELfddhsPPfQQBw8e5N5772Xx4sWsXr2aF198kYkTJ7Jlyxbuu+8+mjVrxoYNG4iPj3c6/PfWrVt56KGHCAwM5KqrrnL1x1NKuZvCXMhKgv73V1u16SWJpdPhxO663WabPjD25UoXz5w5k6+//ppVq1bx5ZdfAjB8+HBuu+02brnlFu68804rtKVLefXVVxk8eHCVw39PmTKFN998k2uuuYZf//rXdftZlFKeL32f9d5azyTcVmXDf+fk5JCdnc0111wDwAMPPFDpcB9KKeVU2ZhNernJiSp+8Tcmxhinw39nZ2djPQ5cKaVqKD0e/FpASMdqq2rDtQtVHBK8/Hxlw3+HhIQQHBzM2rVrAfjggw8aPnCllHtL2wuRPcCr+hSgScKF7rnnHl555RUGDBjAwYMHefDBB5k2bRr9+/entLS00uG/33//fR5//HGuuOIKmjWruvuaUkqdx+GA1F3Qpu/F1TfGVPkC5gDpwJ5yZa8A+4BdwOdASLllvwUSgf3AjeXKx9hlicD0cuUxwPfAAWAB4GeX+9vzifbyTtXFaoyhe/fupqL4+PgLyurb6dOnG3yf1XHF38FTrFq1ytUhqFrSY2g7ecCY54KM2TrvvGJgi3HynXoxZxJz7S/48pYDvY0xfYEf7MSAiPQE7gF62eu8IyLeIuINvA2MBXoCE+26AH8DXjPGdAOygIft8oeBLGNMV+A1u55SSqnaSN1hvbftf1HVq00Sxpg1QGaFsm+MMSX27EYg2p4eB3xsjCk0xhzGOgsYYr8SjTGHjDFFwMfAOLFaYK8DPrXXnwfcXm5b8+zpT4HrRVtslVKqdo5vB29/q03iItRF76aHsC4LAURhJY0yKXYZQHKF8qFAKyC7XMIpXz+qbB1jTImI5Nj1T1UMQESmAlMBIiIiiIuLO295cHBwpc+Mri+lpaUNvs/qFBQUXPC3URcnLy9P/3ZuTo+hpV/CarwDO7Dtu3UXVb9WSUJEfg+UAGVdbJz90jc4P2MxVdSvalsXFhozC5gFEBsba0aOHHne8oSEBFq0aNGgXUdzc3Np2bJlg+2vOsYYAgICGDBggKtDcUtxcXFU/O9KuRc9hoAxsOEo9B5/0X+LGvduEpHJwC3AfXajB1hnAu3LVYsGjldRfgoIERGfCuXnbcteHkyFy14XKyAggIyMDH4Ms2kxxpCRkUFAQICrQ1FKuVLmISjMgXYX1x4BNTyTEJExwLPANcaY/HKLFgMfisg/gHZAN2AT1llBNxGJAY5hNW7fa4wxIrIKuBOrnWIysKjctiYDG+zl35oafstHR0eTkpLCyZMna7J6jRQUFDSqL+WAgACio6Orr6iU8lyX2GgNF5EkROQjYCQQLiIpwHNYvZn8geX2JZyNxphpxpi9IrIQiMe6DPW4MabU3s7PgWWANzDHGGM/8YJngY9F5EVgOzDbLp8N/FtEErHOIO656E9Vga+vLzExMTVdvUbi4uL00o5SqnFJ3Qlevhc1HEeZapOEMWaik+LZTsrK6r8EvOSkfAmwxEn5IazeTxXLC4C7qotPKaXURTq+wxrUz8fvolfRO66VUqopMMY6k7iES02gSUIppZqGrCQoyL6kRmvQJKGUUk1D6k7rXc8klFJKXSB1B3j5QOtel7SaJgmllGoKju+whuLw8b+k1TRJKKWUpzPGOpO4xEtNoElCKaU8X04ynM265EZr0CShlFKe73jZndaXfoOvJgmllPJ0x7eDeFs30l0iTRJKKeXpUjZDmz7ge+mPO66L50kopTzc6YJi9qXmkpB6moTU06SdLsBLBBHBS6BtcAC9o4LpGx1Cl4jm+Hjr789Go7QEjm2D/vfWaHVNEkqpCxhjOJCex/L4NJbHp7EjOfvcstBAX6JCrV+kDgc4jGFd4inmbTgCgL+PF4F+3uee39LM15uru4czpndbhndpha8mkIZ1MgGKz0D7C4bIuyiaJJRq4owxbDqcye5jORzNzOdIRj6J6Xkcyz4LQL/oYH51Q3f6RgfTo20QrYP8L3iAV6nDcPjUGXYfyyYhNZeC4lKMAYMhI6+IxTuO89GmZIICfLiyazgdwgKJCm1GVEgzerQNol3IpV8GURcpeZP1Hn15jVbXJKFUE1Vc6uDLXcf515rDxKeeBqBlgA8dWwXSv30Ij43swg09WtMmuPrnonh7CV0jW9A1sgV3OOlAU1BcyncHTrF0TyrbjmSxMiGdolLHueWdWgVyRZdwhnUOw0uE9NxC0k8XkHO2mPAW/rQJDqBdSAAdwgLpHN4CLy993P1FS9kMzSMgtFONVtckoVQTk51fxILNycxdn0RqTgFdI1vw95/0ZVTP1oQE+tbLY34DfL0Z1bM1o3q2BsDhMJzKKyQ56yw7krPZcPAUX+48zkebjp5bx8/bi6BmvmTlF1Hq+PF5Y8HNfBnYIYTBncIY2CGUPtHBtPDXr7JKJW+C6CFQw+Oqf1mlmoh9J04zb30Sn28/RkGxgys6t+Ivd/Thmu4RDf7L3MtLiAwKIDIogEEdQ3n4qhhKSh3sO5GLn48XkS39CW5mJaxSO6Eczz5LYnoeW49kseVIFqv27wes776uES3o1z6Eadd0pmtk43m2vMudyYDMgzDg/hpvQpOEUh4q80wRGw9lsP7gKdYfzODQyTME+Hpxx4AoJg/vxGVtglwd4nl8vL3oHRV8Qbm3l9A6KIDWQQEM6BDKXYPbA5B1pogdKdnsSs5hZ0o2y/acYO2BUyz6+ZW0Dmo8jw52qWNbrPcaNlqDJgmlPM6x7LO8ufIAn2xNodRhaO7nzZCYMO4b2pHxA6IIbX7xTyVrzEKb+3FtbCTXxkYCkJB6mjvfXc8j87aw8KdX0MzP28URNgLJm6yb6NrV/FHKmiSU8hDpuQW8s+ogH35vXdd/YFhHbu3Xjr7RwU2i22mPtkHMmDiAR+Zv4VcLdvDOfQO1gTtlE7TpDX7Na7wJTRJKeYClu1P5zae7yC8u5e7B0fz8um5ENcFupdf3aM3vb+rBi18l8Mo3+/nNjbH10hDvFhyl1k10/e6p1WaqTRIiMge4BUg3xvS2y8KABUAnIAm42xiTJdbReAO4CcgHHjTGbLPXmQz8wd7si8aYeXb5IGAu0AxYAjxpjDGV7aNWn1YpD1NU4uCvSxN4f10S/dqH8Nrd/egc0cLVYbnUw1fFcOjUGd6NO8hnW1O4vFMYgzqGUppVSpfMfCJa+hPg2wQuRaXHQ1Ge1bOpFi7mTGIu8BYwv1zZdGClMeZlEZluzz8LjAW62a+hwLvAUPsL/zlgMGCArSKy2P7SfxeYCmzEShJjgKVV7EMpBaRk5fP4h9vZmZzNlCs78duxPfDz8fzLStUREZ6/rRd9ooLZeCiDLUlZfLU7FYCXvl8FQAt/H9qHBdK7XRB9ooPpHRVM+9BAQgJ9PefSXMpm6719zW6iK1NtkjDGrBGRThWKxwEj7el5QBzWF/g4YL4xxgAbRSRERNradZcbYzIBRGQ5MEZE4oAgY8wGu3w+cDtWkqhsH0o1eSsT0nhq4U4cDsO79w1kbJ+2rg6pUfH19mLikA5MHNIBgNScsyxYto52MbGczCvkZG4hh06dYeW+dD7ZmnLeui0DfIho6c8Dwzpy/7CO7ps0kjdDYDiExtRqMzVtk2htjEkFMMakikikXR4FJJerl2KXVVWe4qS8qn1cQESmYp2NEBERQVxcXA0/Vt3Jy8trFHGoutFYjmepw/DZgWKWHC6mY5AXj/f3p1nGfuLi9rs6tEava2ABLc4cJFKAIBgZBCbGh8wCb5JOO8gqMJwpNuQVG46cPsvz/4vnX98mcM9lfvSNcL/m2yE/xJEf2Jk9q1fXajt1/cmdtRCZGpRfEmPMLGAWQGxsrBk5cuSlbqLOxcXF0RjiUHWjMRzP5Mx8nl64k01J+dw7tAN/uqVn07i2Xkcu5RgaY1iRkM6LX8Xzj635jOkVylv3DnCf0W3zMyHuOIHDH2HkiJG12lRNk0SaiLS1f+G3BdLt8hSgfbl60cBxu3xkhfI4uzzaSf2q9qFUk3Iqr5C3VyXywcajeHsJr03oxx0DoqtfUdWYiDCqZ2uu7h7OGysO8E7cQb47cIprL6v0gkbjUtYeUctGa6j5Q4cWA5Pt6cnAonLlk8QyDMixLxktA0aLSKiIhAKjgWX2slwRGWb3jJpUYVvO9qFUk5BfVMJry3/gmr+vYv6GI/xkUBTfPnONJogG5O/jzS9v6E5Ycz8+2Zpc/QqNRdJa8PaD6MG13tTFdIH9COssIFxEUrB6Kb0MLBSRh4GjwF129SVY3V8TsbrATgEwxmSKyAuAnd74c1kjNvAYP3aBXWq/qGIfSnm8FfFpPLd4L8eyz3JTnzY8PTqWLk28a6ur+Pl4Ma5/Oz7YeJSsM0Xuccd60lqIGlyjJ9FVdDG9myZWsuh6J3UN8Hgl25kDzHFSvgXo7aQ8w9k+lPJkx7LP8n+L97I8Po3urVvwybQruLxTmKvDavLuGtSe99clsXjncSYP7+TqcKpWcBpSd8CIp+tkc+7XZK+UB0pMz2P22kN8tu0Y3iJMH3sZD18V477dLz1Mz3ZB9GoXxCdbkxt/kkj+HowDOl1VJ5vTJKGUi5SUOth4KJO565NYkZCGn48XPxkYzePXdiE6NNDV4akK7hoUzf/9L56E1NP0aNu4RtA9T9J34OVbJ43WoElCqQblcBg2Hs7gq12pfL3nBBlniggN9OWJ67sx6YqOhLfwd3WIqhLj+kfxlyX7+GRLCn+6taerw6lc0jqIGgR+dfNDQ5OEUg0kNecsT360g01JmTTz9eb6HpHc0rctI2Mj9X4HNxDa3I8bekbyxY5jTB97WeMcAqUwF45vh6t+WWeb1CShVANYmZDGM5/spLDEwV/u6MMdA6L0eQdu6K5B7Vmy+wTf7ktnTO82rg7nQsnfgymts/YI0CShVL0qKXXw8tJ9vLf2MD3aBvHWvQO0K6sbG9EtnMiW/sxee4gR3cJp3tierZ20Frx8oP3QOttkIzxfUsozGGP43ee7eW/tYSZd0ZHPfzZcE4Sb8/H24lejurPlSBa3vLmWPcdyXB3S+ZLWQbuBtXrIUEWaJJSqJ6+vOMDCLSk8cV1X/jyut7Y7eIiJQzrw0aPDKCgu5Y531vGvNYdwOC55yLm6V5gHx7dBpyvrdLOaJJSqBws2H+WNlQe4c1A0vxrV3dXhqDo2rHMrlj45gusui+SlJQnc+95GDp8649qgkr8HR0mdtkeAJgml6tyq/en87vM9jOgWzl/H92m6j8/0cCGBfsy8fxAvj+/D3uOnufH1Nby9KpGiEodrAjqyDsS7TtsjQJOEUnVqR3I2j3+wjdjWLXn3/kF6x7SHExHuGdKBlU9dww09Inll2X5ufXMtSa44q0haC+0GgH/LOt2s/hesVB05eDKPKe9volULP+ZOuZwWja3ni6o3kUEBvHPfIP41aTDpuQU8NHczOfnFDRdA0Rk4VvftEaBJQqk6cSKngEmzN+ElwvyHhhIZFODqkJQLjOrZmn8+MJjkrHwe+2ArxaUNdOkpaS04iqHzyDrftCYJpWopJ7+YyXM2kZ1fxNwpQ4gJr7vuh8r9DIkJ46/j+7L+YAZ/WrQXa3Dsepa4EnyaQYfhdb5pPR9WqhYyzxTx0NzNHD51hvenXE6f6GBXh6QagTsHRXPoZB7vxB2ka2QLHr4qpn53eHCl1avJt+7PYPVMQqkaSs7M58531xOfepq37h3AlV3DXR2SakSeGR3LmF5teOHLeJ5asIO00wX1s6OsJMhIhK718/gdTRJK1cCeYznc8c56Ms4U8cEjQxndqxGO46NcystLeP2e/jx+bRe+3JXKta/G8faqRAqKS+t2R4krrfcu9ZMk9HKTUpdo7YFT/PTfWwhu5stHj15Bt9Z12+VQeY4AX29+feNlTBjcgb8sSeCVZft5f91hRvVsw429WjO8S3jtR5M9+C0Et4fwbnUTdAWaJJS6BF/vOcETH22nc0Rz5k4ZQptg7cWkqtehVSAzHxjE+oOn+OD7oyzacYyPNh2lZYAPt/Zrx9QRnelUkw4PpcVwaDX0Hg/1dNNmrZKEiPwKeAQwwG5gCtAW+BgIA7YBDxhjikTEH5gPDAIygAnGmCR7O78FHgZKgSeMMcvs8jHAG4A38J4x5uXaxKtUbXy+PYVnPtlFn6hg5k0ZQnCgr6tDUm5meJdwhncJp6C4lHWJp/hqdyqfbknh401HGdu7LdOu6XJpnR9SNkNRbr21R0At2iREJAp4AhhsjOmN9UV+D/A34DVjTDcgC+vLH/s9yxjTFXjNroeI9LTX6wWMAd4REW8R8QbeBsYCPYGJdl2lGtzKo8X8asFOhsaE8Z9HhmqCULUS4OvN9T1a84+7+7P22WuZenUX1vxwklvfWsubKw9c/IYSV1pDccRcU2+x1rbh2gdoJiI+QCCQClwHfGovnwfcbk+Ps+exl18v1qA244CPjTGFxpjDQCIwxH4lGmMOGWOKsM5OxtUyXqUu2b83HuHf8UXc0COSOQ/qndSqbkUGBTB97GWs++113N6/Hf9v+Q+8v+7wxa18cCVEXw7NQuotvhonCWPMMeBV4ChWcsgBtgLZxpgSu1oKEGVPRwHJ9roldv1W5csrrFNZuVIN5us9J/jToj30i/Dm3fsH6XDfqt4EBfjy6l39GN2zNc//L57PtqZUvcKZDDi+o14vNT5qbrkAACAASURBVEEt2iREJBTrl30MkA18gnVpqKKy2w2dtaqYKsqdJTCnty6KyFRgKkBERARxcXFVhd4g8vLyGkUcquZ+yCrl75sLiAnyYnLXEtZ9t8bVIalacJd/k3dGGVLSvPj1pzs5kriPga2df01Hpq2hJ4atOSHk1uPnqs158w3AYWPMSQAR+S8wHAgRER/7bCEaOG7XTwHaAyn25algILNceZny61RWfh5jzCxgFkBsbKwZOXJkLT5W3YiLi6MxxKFq5kBaLk+8u54OYc359LHh7Nq8Xo+nm3Onf5NXXlXCfe99z8xdp3nh9u5MuLzDhZU+XwDNQhl0y8PgVX9nuLVpkzgKDBORQLtt4XogHlgF3GnXmQwssqcX2/PYy7811qAmi4F7RMRfRGKAbsAmYDPQTURiRMQPq3F7cS3iVeqi7DtxmklzNuHv6828h4YQ1tzP1SGpJqa5vw/zpgxhSEwYz362m19/spOzReVuwnOUQuJy6HJdvSYIqF2bxPdYDdDbsLq/emH9mn8WeEpEErHaHGbbq8wGWtnlTwHT7e3sBRZiJZivgceNMaX2mcjPgWVAArDQrqtUvflqVyp3vL2eUodh7pTLaR8W6OqQVBMVHOjLvIeG8MR1Xflkawp3vLPux6ffpWyBMych9qZ6j6NW3TSMMc8Bz1UoPoTVM6li3QLgrkq28xLwkpPyJcCS2sSo1MUodRhe/WY/78YdZGCHEGbeP0iH+1Yu5+0lPDU6lgEdQ/nVgh3cNXM9K566hpB9X4KXL3QbVe8x6NhNqslLzsxn8pxNvBt30HrI/dRhmiBUo3JtbCQfPDKUzDNF/L9l+2HfVxAzAgLqf9Rh7fCtmqziUgez1x7m9RU/4C3CX8f3YeIQJw2ESjUCvdoFM+mKTqzduA78DsIVP2uQ/WqSUE3SvhOn+eXHO9h3IpdRPVvz/G29aBfSzNVhKVWlX43qTviOt8ABjm5jG+RSkCYJ1eTkFhTzyLwtFBQ7mHn/IMb01mG+lXsIbubLvUG72ZHZmR8OOLj78vrfp7ZJqCbnhS/jOZ59ln8+MFAThHIvp1MJy95NfNAIXv56H9n5RfW+S00SqklZHp/Gwi0pTLumC4M6hrk6HKUuzX6rs+fQmyaRc7aYZz7ZyemC4nrdpSYJ1WRk5BXy2//uokfbIH55Q3dXh6PUpdu/BMI606XHIH5/Uw9W7T/JTW98x9YjmfW2S00SqkkwxvC7z3dz+mwJr03oV/ungSnV0ApOWw8YuuxmEOGhq2JY+NMrEIG7/7mR11f8QEmpo853q/9SVJPw2bZjLNubxtOju3NZmyBXh6PUpUtcDo5iuOyWc0WDOoay5IkR3NavHa+vOMCEWRs5mpFfp7vVJKE8XtKpM/xp0R6GxoTxyIjOrg5HqZqJXwzNI6znR5TTMsCX1yb05417+vNDWi5j31jDws3JWEPj1Z4mCeXRikocPPHxdny9vXhtQn+8vernOcBK1avCXPjha+h5e6UD+o3rH8WyX15N3+gQfvPZLn76761k5BXWeteaJJRHe23FD+xKyeHl8X30ZjnlvvYtgZIC6HNnldXahTTjg0eG8oebexC3/ySjX1vDsr0narVrTRLKY61PPMXM1QeZOKQ9Y/u0dXU4StXc7k8guANEXzB26gW8vIRHRnTmyyeuok1wAD/991aeXljzrrKaJJRHyjlbzK8W7qBzeHP+eEtPV4ejVM2dyYBDq6D3ePC6+K/s7q1b8vnPruSJ67ryxY5jjPrHahZuTqbUcWltFZoklEd6bfkPnMwt5PUJAwj009FnlBuL/wIcJdVeanLGz8eLp0bH8tljw2kT3IzffLaLMa+vYUV82kU3bGuSUB5n/4lc/r3xCPcO7UCf6PofSlmperXnMwiPhda9a7yJ/u1D+OJnw3nnvoGUOAyPzN/Co/O3kF9UUu26miSURzHG8NziPbQM8OHpUbGuDkep2slJgSPrrbMIqV3PPBHhpj5t+eZXV/P7m3rw7b507v3X92SeqXr8J00SyqN8tTuVjYcyeWZ0LKH6bGrl7vb8FzDQ+yd1tklfby8evboz794/iITU09w5cz0pWZXfgKdJQnmM/KISXvoqgV7tgvThQcoz7PkU2g2EVl3qfNM39mrDvx8eyqncQn7y7vpK62mSUB7j7VWJpOYU8PxtvfSmOeX+TiVC6s4aNVhfrCExYXwybXiVdWqVJEQkREQ+FZF9IpIgIleISJiILBeRA/Z7qF1XRGSGiCSKyC4RGVhuO5Pt+gdEZHK58kEistteZ4ZILS/KKY+VdOoM/1pzmDsGRDG4kw4BrjzAzg9BvKDX+HrdTWyblix5YkSly2t7JvEG8LUx5jKgH5AATAdWGmO6ASvteYCxQDf7NRV4F0BEwoDngKHAEOC5ssRi15labr0xtYxXeagXvozH11v47djLXB2KUrVXWgLbP4BuoyGo/m8EbdXCv9JlNU4SIhIEXA3MBjDGFBljsoFxwDy72jzgdnt6HDDfWDYCISLSFrgRWG6MyTTGZAHLgTH2siBjzAZjdeidX25bSp3z7b40Vu5L58kbuhEZFODqcJSqvcTlkHcCBk5ydSS1esZ1Z+Ak8L6I9AO2Ak8CrY0xqQDGmFQRibTrRwHJ5dZPscuqKk9xUn4BEZmKdcZBREQEcXFxtfhYdSMvL69RxOHpih2GP6w9S5vmQueSo8TFJVe/Ug3o8XR/7nQMe+/+By39Qtl43B9zIs6lsdQmSfgAA4FfGGO+F5E3+PHSkjPO2hNMDcovLDRmFjALIDY21owcObKKMBpGXFwcjSEOT/f2qkTS8vcz/6EhXN09ot72o8fT/bnNMTydCqu3wpVPcM11N7g6mlq1SaQAKcaY7+35T7GSRpp9qQj7Pb1c/fbl1o8GjldTHu2kXCkAjmef5a1vE7mxV+t6TRBKNaidH4IphQEPuDoSoBZJwhhzAkgWkbLbWq8H4oHFQFkPpcnAInt6MTDJ7uU0DMixL0stA0aLSKjdYD0aWGYvyxWRYXavpknltqWaOGMMf1q0F4cx/OFmHcBPeQiHA7b9GzqNqJd7I2qitiOf/QL4QET8gEPAFKzEs1BEHgaOAnfZdZcANwGJQL5dF2NMpoi8AGy26/3ZGFP2VO/HgLlAM2Cp/VKK2WsPsyIhjT/e0pP2YYGuDkepunFkLWQdhmt/5+pIzqlVkjDG7AAGO1l0vZO6Bni8ku3MAeY4Kd8C1HxUK+WRth3N4uWl+xjdszUPXdnJ1eEoVXe2zYeAYOhxq6sjOUfvuFZuJTu/iF98uJ02wQG8cmc/9P5K5THyM63nWPedAL6N5ymKOtC+chvGGJ75ZCfpuQV8Om04wYG+rg5Jqbqz9X0oLYRBU1wdyXn0TEK5DasdIp3f3dSDfu1DXB2OUnWntBg2vQedR0LrxtURQ5OEcgu7U3L429dWO8SDwzu5Ohyl6lb8Isg9DsN+5upILqBJQjV6eYUl/OKjbYS38Ofvd/bVdgjleTa+C2FdoOsoV0dyAU0SqtF7btFejmbm89qE/oQE6oOElIdJ3gzHtsCwx8Cr8X0lN76IlCrni+3H+GxbCj+/rhvDOrdydThK1b2N74B/MPSb6OpInNIkoRqt5Mx8/vDFHgZ3DOWJ67q6Ohyl6l5OitUeMWgS+LdwdTROaZJQjVKpw/DUwh0I8Po9/fHx1v9UlQfa9C/AwJCpro6kUnqfhGqUZq05xOakLF6b0I/oUB12Q3mgghzr3oget0JI430mu/48U43O3uM5/GP5fm7q04bb+zt9hIhS7m/TLCtRXPWUqyOpkiYJ1agUFJfyqwU7CA3046Xb+2h3V+WZCnNhw9vQfQy06+/qaKqkl5tUo/Lqsv38kJbH3CmXE9pcu7sqD7V5NpzNgqt/4+pIqqVnEqrR2HQ4k9nrDnP/sA6MjI2sfgWl3FHRGVj/JnS5HqIHuTqaammSUI3C2aJSfvPpTqJCmvHbsT1cHY5S9WfrXMg/Bdc0/rMI0MtNqpF49Zv9JGXk8+EjQ2nur/9ZKg9VfBbWvQExV0OHYa6O5qLomYRyua1HMpljX2Ya3jXc1eEoVX+2zYe8NLdoiyijSUK5VEFxKb/+ZJdeZlKerzAX1rwCHa+ETle5OpqLpuf1yqVeW/EDh06d0ctMyvOtmwFnTsLEBeBGXbtrfSYhIt4isl1EvrTnY0TkexE5ICILRMTPLve35xPt5Z3KbeO3dvl+EbmxXPkYuyxRRKbXNlbVuCRn5jNn7WHuGhStl5mUZzudChvegl7j3aJHU3l1cbnpSSCh3PzfgNeMMd2ALOBhu/xhIMsY0xV4za6HiPQE7gF6AWOAd+zE4w28DYwFegIT7brKQ7z6zX68vYSnR8e6OhSl6lfcX6ynz13/J1dHcslqlSREJBq4GXjPnhfgOuBTu8o84HZ7epw9j738erv+OOBjY0yhMeYwkAgMsV+JxphDxpgi4GO7rvIAe47lsGjHcR6+KoY2wQGuDkep+pOeANv/A0MehbAYV0dzyWp7JvE68BvAYc+3ArKNMSX2fApQNvhOFJAMYC/PseufK6+wTmXlys0ZY/jr0gRCA3356TVdXB2OUvVr+XPg1xKu/rWrI6mRGrcUisgtQLoxZquIjCwrdlLVVLOssnJnCcw4KUNEpgJTASIiIoiLi6s88AaSl5fXKOJojHafLGFdYiH3XebHto3rXB3ORdHj6f5ccQxDsnbS/8AyDnaeTPKmXQ2677pSm+4kVwK3ichNQAAQhHVmESIiPvbZQjRw3K6fArQHUkTEBwgGMsuVlym/TmXl5zHGzAJmAcTGxpqRI0fW4mPVjbi4OBpDHI1NqcPw8ozv6BDmzXP3X4Ofj3v0wtbj6f4a/BiWFMK7z0BIR7pM/DtdfN3zsmqN/4UaY35rjIk2xnTCanj+1hhzH7AKuNOuNhlYZE8vtuexl39rjDF2+T1276cYoBuwCdgMdLN7S/nZ+1hc03hV4/D59mPsO5HLMzfGuk2CUKpG1s+AjANw8z/ATRME1M99Es8CH4vIi8B2YLZdPhv4t4gkYp1B3ANgjNkrIguBeKAEeNwYUwogIj8HlgHewBxjzN56iFc1kPyiEl5dtp9+0cHc0qetq8NRqv5kHoY1r0LPcdDtBldHUyt1kiSMMXFAnD19CKtnUsU6BcBdlaz/EvCSk/IlwJK6iFG53r/WHObE6QLevHcAXl7uczORUpfEGFjyDHj5wJiXXR1Nren5vmoQaacLmLn6IDf1acPlncJcHY5S9Sd+ESSugOv+AEHtXB1NrWmSUA3i1WX7KXUYnh1zmatDUar+FOTA19OhTR+4/FFXR1MnNEmoerfnWA6fbkvhwSs70bFVc1eHo1T9Wfos5KXDLW+At2eMRaZJQtUrYwwvfZVASDNfHr+2q6vDUar+xC+GnR/BiKfdbnymqmiSUPXqm/g0NhzK4Jc3dCe4ma+rw1GqfuSmwZe/hLb93eaJcxdLk4SqNwXFpbzwZTzdW7fgvqEdXB2OUvXDGPjfE1CYB+Nngbdn/RjyjItmqlH615pDpGSd5cNHhuLjrb9HlIfaNh9++Bpu/CtEeN6IxvovV9WL49lneTsukZv6tNFnRSjPlZ4AX//Wemb10GmujqZeaJJQ9eIvSxIwBn53kz6SVHmogtOw4AHwaw53zAIvz/w69cxPpVxq46EMvtyVymMjuxAdGujqcJSqe8bA4p9D5iG4630I8txhZrRNQtWpvMISnlu0l6iQZkzTZ0UoT7XhbevO6lEvQKerXB1NvdIkoepMXmEJD87ZROLJPN6bPJgAX29Xh6RU3UtaB8v/BD1uheG/cHU09U4vN6k6caawhIfe38z25Gxm3DOAa2MjXR2SUnUv4yAsfABCO8G4t0E8f6BKTRKq1vKLSpgydzNbjmTy+oT+3NzXc6/PqiYsPxM+uMtqj7jvEwgIdnVEDUIvN6laKSpx8Oj8LWxJyuS1Cf25tZ/7j3qp1AWKC+CjiZCTApMXQ6um096mSULVmDGGZz/bxbrEDF65sy/j+ke5OiSl6p7DAV88Bskb4c73ocMwV0fUoPRyk6qxV7/Zz+fbj/H0qO7cNbh99Sso5W6MgWW/hb3/hRueh97jXR1Rg9MkoWrkw++P8vaqg9xzeXt+fp2O7qo8kDFWL6bvZ8Kwn8GVT7o6IpfQJKEu2Td7T/DHRXsYGRvBi7f3RppADw/VBK16CdbPgMsfgRv/0iR6MjmjSUJdkkU7jvHYB9voHRXM2/cO1IH7lGda/XdY8woMnARjX2myCQJqkSREpL2IrBKRBBHZKyJP2uVhIrJcRA7Y76F2uYjIDBFJFJFdIjKw3LYm2/UPiMjkcuWDRGS3vc4M0Z+sLvXB90f45YIdDO4YygePDKW5v/Z7UB7GGFjxvHUW0e9e6wlzHjom08WqzacvAZ42xvQAhgGPi0hPYDqw0hjTDVhpzwOMBbrZr6nAu2AlFeA5YCgwBHiuLLHYdaaWW29MLeJVtTBz9UF+//kero2NZN5DQ2ihCUJ5mtIS67kQa/8Bgx6EcW81+QQBtUgSxphUY8w2ezoXSACigHHAPLvaPOB2e3ocMN9YNgIhItIWuBFYbozJNMZkAcuBMfayIGPMBmOMAeaX25ZqQP9cfZCXl+7jlr5tmXn/IB1uQ3me4gL4ZLL1bIirfwO3vA5e+t851NF9EiLSCRgAfA+0NsakgpVIRKRsfIYoILncail2WVXlKU7Kne1/KtYZBxEREcTFxdXq89SFvLy8RhFHba1JKWbOniKGtPFmfNsc1q9d4+qQXMJTjmdTVtkx9Ck+Te89fyMkZw8Huj7KMa8rYfXqhg+wkap1khCRFsBnwC+NMaeraDZwtsDUoPzCQmNmAbMAYmNjzciRI6uJuv7FxcXRGOKojWV7TzB32VZGdAtn9uTL8fNpuqfennA8mzqnxzBtL3z0BOSdgJ/MplufO+nmkugar1r9qxcRX6wE8YEx5r92cZp9qQj7Pd0uTwHK33EVDRyvpjzaSblqABsOZvCLj7bTr30I/3xgUJNOEMpDxS+G90ZBSSFMWQp97nR1RI1SbXo3CTAbSDDG/KPcosVAWQ+lycCicuWT7F5Ow4Ac+7LUMmC0iITaDdajgWX2slwRGWbva1K5bal6tCslm0fnb6FTq0Def/ByAv20kVp5kNIS+PZFazTXyB4wNQ6iB7k6qkarNv/6rwQeAHaLyA677HfAy8BCEXkYOArcZS9bAtwEJAL5wBQAY0ymiLwAbLbr/dkYk2lPPwbMBZoBS+2XqkcH0nKZPGcTIYG+zH9oKCGBfq4OSam6k30UPnvUGoep//1w8/8D3wBXR9Wo1ThJGGPW4rzdAOB6J/UN8Hgl25oDzHFSvgXoXdMY1aVJzszn/tnf4+PtxQePDKVNsP7jUZ4j/OR6mDnJGrBv/HvQ967qV1I6CqyypJ0u4L73vqeg2MHCn15Bx1bNXR2SUnUjPxOW/Z7eez+EqEHwk9kQFuPqqNyGJglFYnouj8zbQkZeIR88OozYNi1dHZJStWeM9RzqJc/A2SySOt5Np0nvgLevqyNzK5okmrhv9p7gqYU7CfD1Yv7DQ+jfPsTVISlVe9nJ8PV02PcltO0H9/+XpP2ZdNIEcck0STRRDodhxrcHeH3FAfpGB/PPBwbRNriZq8NSqnaK8mHdG7DudWv+hufhip+Dtw/sj3NpaO5Kk0QTdLaolKcW7mDpnhP8ZGA0L93RW4faUO7N4YA9n8GK5+D0Meg1HkY9DyEdXB2Z29Mk0cSk5xbw6Pyt7ErJ5g839+Dhq2L0eRDKfRkD+5fAty9B+l5o0xd+8h50HO7qyDyGJokm5Ie0XKa8v5nMM0X88/5BjO7VxtUhKVUzDgckLoe4v8Lx7RDWxeq11OsOHZivjmmSaCJW7UvniY+2E+DnzcKfXkGf6GBXh6TUpSspgt2fwPo34WQCBHeAcW9D33usdgdV5/Sv6uEcDsOb3yby+sof6NEmiH9NHkxUiDZQKzeTmwbb58Pm2ZCbCpG94I5/Qu+faJfWeqZJwoPlnC3mqQU7WLkvnfEDonjpjj4089NTceUmHA5I+g62zLG6sjpKoPNI62FAXa5v0o8UbUiaJDxIYUkpCam57E7JZvexHL47cIqTuYU8f1svJl3RURuolXs4dQB2fgy7FkLOUQgIgaHTYNAUCO/q6uiaHE0SHsAYw+Kdx3nxqwRO5hYCENbcjz5RwcyYOIDLO4W5OEKlqpF5yBq6O34RHN8G4mWdNVz3B+h5G/jqJVJX0STh5hLTc/njF3vZcCiDvtHBPH9bL/pGBxMV0kzPHFTj5XBA6g448A0kfAlpu63ytv1h9IvQ+04IauvaGBWgScJtncwtZObqg8zfkEQzX29evL03E4d0wNtLE4NqpE6nWm0MB7+FxBVw5iQg0H4IjH4JetwKoR1dHaWqQJOEmzmVV8isNYeYvyGJohIHdw6K5jdjLiO8hb+rQ1PqfNnJkPw9HN0Ah9fAqR+s8mahVsNz9xuhy3XQPNy1caoqaZJwE4dPnWHe+iQWbE6msKSU2/tH8YvruxETrkN6q0ag4DSk7rRubDu+DZI3WcNjAPg2t+6AHvAAxFwNbfroDW9uRJNEI1ZQXMr3hzOZvz6Jb/en4+Ml3NqvHY9f25UuES1cHZ5qihwOq8dRegKc2ANp9ivjIGCsOsHtof1Q6DDMem/dW290c2N65BqRzDNFbD+axeakLLYkZbIrJYeiUgfhLfx44rpu3DesA5Et9WlxqgGczYbMg5BxyH5PhJP7re6pJWd/rBcaA617Qd8J0G4gtOuvl488jCYJFzmVV8gPabn8cCKX3cdOs/1oFodOnQHA11voExXMlCs7MbhTGCO6hesoraruOErhzCnrclBOyo/v2Ucg64j1XpBTbgWB4GgI7w6dRkBELERcBq17gr8+oMrTaZKoJ8YYcvKLSc8tIDWngEMn8ziQnkdiuvWeeaboXN1Wzf0Y0CGUuwa3Z2CHEPpGh+id0erSFBfA2SzIPwX5GVYSyM+AvHQ4kw55JyEvDXJPWO+m9Pz1fQKsYbVDOlq9jUI6QlhnaNUFQjvpfQpNWKNPEiIyBngD8AbeM8a87OKQzlNS6iAzv4hDJ8+w51gOu1Jy2HM8h6MZ+ZQs++a8ukEBPnRv3ZLRPVvTrXVLurduQffWLYls6a/3NDRVjlIoOmO9ivOhKA8K8+z3XOu94LQ1XXjami7Itn7pn822ps9mWes6I17QPAKaR0KLSIjsad1/0LINtGwHwVEQFA2BYTrMhXKqUScJEfEG3gZGASnAZhFZbIyJr6t9GGMoLHFQWOzgbHEpZ4pKOFtUypnCEnILSsgtLCa3oITTZ4s5lVfEybxCMvIKOZVXREZeIVn5xedtr21wAL2jgunevIjBvboSGRRA65b+xEQ0J6KFJoMGYYz15WtKK7w7rPF/HCXW/Ln3Ymu6tNjJvPUekb4TdqVDaZH9KoaSwh/ny6ZLCu3pQig+a88X2NP2e/FZ60u9ON9a56II+AdBQJA1TEVAMITFWO/NQn98Bbay2gQCw63pwDDtSaRqpVEnCWAIkGiMOQQgIh8D44DKk0ReGhte+QkOYzDG4Ch7OaDUGBwOg8NAqXFgHIZSY5xupuyr3AsItl+dvQR/X2/8fbzw9xH8g7zwC7PmA/28CQn0w9/HC4whPTedyBORcMLZ1ivs87wYqlpWzfLK6p4rr2reybILpo2T/dllxuGknrN3h72rsuny5Y4ft1V+WdkX/LlX6Y/Lz1tmJ4OKf6M60Auq+q8OvP2sSzbefuDjb037BPw4HRAEPq3BNwB8moFfIPiWvZqBfwurq6hfIPi1sK71+7Wwyv2DrGkvrzr/XEpVp7EniSggudx8CjC0YiURmQpMBejT1o8OZ3ZZX/L2/8m5ej++ixfgbS+1/me9yqZFzk17lZs+p9R+FZ5fVHbSH2gM+XmHK/1gpsoziorLzp+vfN2yT1G+bsVtVLat8uXiZD/Wtp2XgZVOy9b5cZvWfNkyrwrlP9Y14gXyY/3yy4142euev9wqK7/M67wy5/M+9rt3uXfnL4eXj13fm7yCIpo1D7LLvDHii8PLB4eXL0Z8aneppth+nSfPfqm6kpeXR1xcnKvDcDuNPUk4+5d3wc9EY8wsYBZAbGysifq//fUdV7Xi4uIYOXKkq8NQdSQuLo6hejzdmv6brJnGfv6aArQvNx8NHHdRLEop1eQ09iSxGegmIjEi4gfcAyx2cUxKKdVkNOrLTcaYEhH5ObAMqwvsHGPMXheHpZRSTUajThIAxpglwBJXx6GUUk1RY7/cpJRSyoU0SSillKqUJgmllFKV0iShlFKqUmIqGZbCXYlILlD+brpgIMfJdG3nq6sbDpyqhzjccd3GEkdt1q3qeNbnfutyW01t3YrzrjqG7rJurDHmwrHfjT3Gkae8gC0V5mc5m67t/EXUrZc43HHdxhJHLdet9Hjq365xrttYjqEbrXve36fs1RQuN/2vkunazldXt77icMd1G0sctf0MrtpvY/j87rius3lXxOEu6zrliZebthhjBmscqi7p8XR/egyrVtnfxxPPJGa5OgBbY4lD1Q09nu5Pj2HVnP59PO5MQimlVN3xxDMJp0SkvYisEpEEEdkrIk9WWP6MiBgRCa+DfY0Rkf0ikigi0+2yD+yyPSIyR0R862k/14nINns/80Sk1kOv2PGmi8ieCuW/sPe/V0T+f3v3HmzXeIdx/PuICHULKkQxUdeqRlwad41MOy1VibvUkLR6YaajtIo2M6pUKUo7pq2pW1B1ibjXLTWIa6JCJBWSKCqkoY27QfHrH++72U72Omefs9be50Sez8yZs/a71l7vu953rf1bl73f9/QK8mnYRpJOlvSYpEcl3S5pnZL5LC9pmqQZOZ9f5PQNJE2VNFfSlblTydIKjcfHawAACtlJREFU2kmSTpE0J2/vkSXzWKyNqq63vM6iNlpd0uRcd5MlrdbCvK7M2/SopGckPVpBXou1Ud28cyRVMrhHQTvtn7fvA0mV3A4ryGeYpAdzvf1d0vCmVtboafYn8Q8YDGydp1cG5gCb59frkToRfBb4dMl8+gFPAZ8FlgNmAJsDe/DRKDuXA0e0KJ/ngE3yMicBh1VQd7sCWwOz6tJ2A/4GDMivB7WqjYBV6pY5Eji3ZD4CVsrT/YGpwPbAVcBBOf3csm3URTt9C7gEWKaK+itoo0rrrYs2Oh04PqcfD/y6VXl1WOY3wAmtaKM8b1vgUuCNstvTSTt9DtgUuAvYtoX53A7snqf3AO5qZl1LzZVERCyIiOl5+nVgNmnkO4CzgWOpZtzLD4dcjYh3gSuAURFxc2TANNLYGFXnsy/wTkTMyctMzmmlRMQUYFGH5COA0yLinbzMixXk07CNIuK1usVWpGQ75WaonRn2z38BjASuzukXA6PL5JM13B9I9XdSRHyQy1Sq/hq1UdX1ltdZdByNItUZVFR3XRyzSBJwAOmkq4yGbSSpH3AG6bOhEgXtNDsiKh0preCYDWCVPL0qTY7Ns9QEiXqShgBbAVMl7QU8HxEzKlp9oyFX63fs/sAhwK0tyGdtoH/dJet+fHzQpiptAuySb8/cLemLVa68vo3y61MkPQccDJxQwfr75dsUL5KC6VPAKxHxXl7kY+1WQtH+sCFwYL7sv0XSxhXktZiq663DuofwURutFRELIH24A4NamFfNLsDCiJhbcvVFbfQD4Ibadn0CHAWckfeHM4GfNvOmpS5ISFoJmESqsPeA8VR78HQ15OofgCkRcU8L8vmANDDT2ZKmAa+TtrEVlgVWI92m+QlwVT6zK62+jWpnwxExPiLWAy4jHbylRMT7ETGMdEU3nHTJv9hiZfOheH8YALwd6SuH5wEXVpDX4hlVXG81jdqoVTrJawzlryKgcRsNAPYHzqlg/X3FEcDReX84GrigmTctVUEin8VPAi6LiGtIZ3MbADMkPUP6wJguae0S2RQOuSrp58CawI9KrL/TfCLigYjYJSKGA1OAsmdZneV/Tb51M40UoKp46N+xjTr6CxXcQquJiFdI94K3BwbWPeivaqjcov1hPmk7Aa4FhlaQV2cqq7eCNlooaXCeP5h0hdaqvMjttA9wZQXZNGqjZ4CNgHn5s+FTkuZVkFdvGgvU6nAi6eSoS0tNkMhnuRcAsyPiLICImBkRgyJiSEQMIe0sW0fEv0tk1XDIVUnfAb4KjKndhy6pKJ9BAJIGAMeRHsC2wnWke/hI2oT0wO8/nb6jC43aKKfX34rZC3iiZD5rShqYp1cAvky6330n6RYdpAPq+jL5ZEVD8H5Yf8CXSA9lK1V1veV1Nmwj0jaNzdOV1F0neUFqsyciYn7ZfGjcRtdFxNp1nw1vRcRGFeTVm14g7WuQ9r3mTiCreJK+JPwBO5Mu8x8DHs1/e3RY5hlKfrspPvrmwBzSfe7xOe29/LqWd6lvZHSSzxmkD7wnSZfnVdTd5cAC4H+kQHoYKSj8GZgFTAdGtqqNSGeSs3L6jaSH2WXyGQo8ktc3q9YWpG+3TAPmkc60BlRUf43aaSDwV2Am8ACwZQvaqNJ666KN1gDuIH3w3AGs3qq88rwJwOFVtE9RG3WYX9W3mxq10955+h1gIXBbi/LZGXiY9O2tqcA2zazLP6YzM7NCS83tJjMz6z4HCTMzK+QgYWZmhRwkzMyskIOEmZkVcpAwM7NCDhJmZlbIQcLMzAo5SJiZWSEHCTMzK+QgYWZmhRwkzMyskIOEmZkVcpAwM7NCDhJmZlbIQcLMzAo5SPSApJB0ad3rZSW9JOmm3iyXlSdp79y+m/V2Wax5PiZbx0GiZ94EtsjjIwN8BXi+OyvIA7lb3zMGuJc0znHTJPVrTXGsSaWPSWvMQaLnbgG+nqfHkMaUBUDScEn3S3ok/980p4+TNFHSjcDt7S+ydUbSSsBOpPGAD8ppIyRNkXStpMclnStpmTzvDUknSZoK7NB7JbesJ8fkPZKG1S13n6ShbS11H+cg0XNXAAdJWh4YShpYvOYJYNeI2Ao4AfhV3bwdgLERMbJtJbVmjQZujYg5wCJJW+f04cCPgS8AGwL75PQVgVkRsV1E3Nv20lpHPTkmzwfGAUjaBBgQEY+1rcRLAAeJHso70hDSGcvNHWavCkyUNAs4G/h83bzJEbGoLYW07hpD+qAh/x+Tp6dFxD8j4n3S2enOOf19YFJ7i2hFenhMTgT2lNQf+DYwoS2FXYL4vng5NwBnAiOANerSTwbujIi9JQ0B7qqb92abymbdIGkNYCTpvnYA/YAgfdhEh8Vrr9/OgcP6jm4dkxHxlqTJwCjgAGDbNpZ1ieAgUc6FwKsRMVPSiLr0Vfnoodm4dhfKemQ/4JKI+H4tQdLdpKuG4ZI2AJ4FDgT+1DtFtCb05Jg8H7gRuMdX+Yvz7aYSImJ+RPyuwazTgVMl3Uc6I7W+bwxwbYe0ScA3gQeA04BZwNMNlrM+oifHZEQ8DLwGXNSGIi5xFNHxStrMavLZ6DERsWdvl8VaQ9I6pNtPm0XEB71cnD7HVxJmttSSdCjpW1DjHSAa85WEmZkV8pVEkyStJ+lOSbMl/UPSD3P66pImS5qb/6+W0w+W9Fj+u1/SlnXr+pqkJyXNk3R8b22TmVlXfCXRJEmDgcERMV3SysDDpB9fjQMWRcRp+QN/tYg4TtKOwOyIeFnS7sCJEbFd7r5hDqnbgPnAQ8CYiHi8N7bLzKwzvpJoUkQsiIjpefp1YDbwGdL3qy/Oi11MChxExP0R8XJOfxBYN08PB+blH2e9S/rR1qj2bIWZWfc4SPRA/jHOVqQHXmtFxAJIgQQY1OAth5H6lYEUWJ6rmzc/p5mZ9Tn+MV035U7gJgFHRcRrkrpafjdSkKh15dDoDb7nZ2Z9kq8kuiH37zIJuCwirsnJC/Pzitpzixfrlh9K+jXnqIj4b06eD6xXt9p1gRdaXXYzs55wkGiS0iXDBaSH0WfVzboBGJunxwLX5+XXB64BDsm9itY8BGwsaQNJy5G6pL6h1eU3M+sJf7upSZJ2Bu4BZgK1H938jPRc4ipgfeBfwP4RsUjS+cC+pP5+AN6LiG3zuvYAfkvqHuDCiDilbRtiZtYNDhJmZlbIt5vMzKyQg4SZmRVykDAzs0IOEmZmVshBwszMCjlImFVM0omSjulk/mhJm7ezTGY95SBh1n6jAQcJWyL4dxJmFZA0HjiU1HnjS6Su5F8FvgcsB8wDDgGGATflea+SfnAJ8HtgTeAt4LsR8UQ7y29WxEHCrCRJ2wATgO1InWZOB84FLqr12SXpl8DCiDhH0gTgpoi4Os+7Azg8IuZK2g44NSJGtn9LzBbnXmDNytsFuDYi3gKQVOuLa4scHAYCKwG3dXxj7lV4R2BiXY/CA1peYrMmOUiYVaPRJfkEYHREzJA0DhjRYJllgFciYljrimbWc35wbVbeFGBvSSvkoW2/kdNXBhbkLuYPrlv+9TyPiHgNeFrS/pB6G64fD92st/mZhFkF6h5cP0saM+Rx4E3g2Jw2E1g5IsZJ2gk4D3gH2I/Uq/AfgcFAf+CKiDip7Rth1oCDhJmZFfLtJjMzK+QgYWZmhRwkzMyskIOEmZkVcpAwM7NCDhJmZlbIQcLMzAo5SJiZWaH/Azrrq4iJB9a7AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def seir_model_italy_i(t, beta):\n", " sigma = 1 / 5.2\n", " gamma = 1 / 2.9\n", " E0 = 19\n", " N = 60461826\n", " return seir_model(t, beta, sigma, gamma, E0, N)[2]\n", "\n", "params, covar = optimize.curve_fit(seir_model_italy_i, it_df.index, it_df.active)\n", "it_df['fitted'] = seir_model_italy_i(it_df.index, *params)\n", "\n", "print(\"Params:\", params)\n", "ax = it_df.plot(\n", " x='date',\n", " y=['active', 'fitted'],\n", " grid=True,\n", " title=\"Fitting active cases to I\"\n", ")" ] } ], "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 }