notebooks/khan2015/python_scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EVD transmission model from Khan et al 2015 using Python/SciPy\n",
    "\n",
    "Author: Verity Hill @ViralVerity\n",
    "\n",
    "Date: 2018-10-03"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.integrate import ode, solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def LIB(times, init, paramsLIB):\n",
    "    beta, eta, psi, pi, p, alpha, thetaI, thetaH, tau, mu, deltaI, deltaH = paramsLIB \n",
    "    SL,SH,E,I,H,R = init\n",
    "    \n",
    "    #N = ST + E + I + R\n",
    "    N = SL + SH + E + I + H + R\n",
    "    \n",
    "    lamb = beta * ((I + (eta*H))/N)\n",
    "    \n",
    "    dSL = pi*(1-p) - lamb*SL - mu*SL \n",
    "    dSH = pi*p - lamb*psi*SH - mu*SH\n",
    "    \n",
    "    dE = lamb*(SL + psi*SH) - (alpha + mu)*E \n",
    "    \n",
    "    dI = alpha*E - (tau + thetaI + deltaI + mu)*I\n",
    "    dH = tau*I - (thetaH + deltaH + mu)*H\n",
    "    \n",
    "    dR = thetaI*I + thetaH*H - mu*R\n",
    "        \n",
    "    return(dSL, dSH, dE, dI, dH, dR)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SL(times, init, paramsSL):\n",
    "    beta, eta, psi, pi, p, alpha, thetaI, thetaH, tau, mu, deltaI, deltaH = paramsSL \n",
    "    SL,SH,E,I,H,R = init\n",
    "    \n",
    "    #N = ST + E + I + R\n",
    "    N = SL + SH + E + I + H + R\n",
    "    \n",
    "    lamb = beta * (I + (eta*H))/N\n",
    "    \n",
    "    dSL = pi*(1-p) - lamb*SL - mu*SL \n",
    "    dSH = pi*p - lamb*psi*SH - mu*SH\n",
    "    \n",
    "    dE = lamb*(SL + psi*SH) - (alpha + mu)*E \n",
    "    \n",
    "    dI = alpha*E - (tau + thetaI + deltaI + mu)*I\n",
    "    dH = tau*I - (thetaH + deltaH + mu)*H\n",
    "    \n",
    "    dR = thetaI*I + thetaH*H - mu*R\n",
    "    \n",
    "    return(dSL, dSH, dE, dI, dH, dR) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For Liberia\n",
    "paramsLIB = [0.371, 0.7, 1.6, 1.7, 0.2, 0.1, 0.1, 0.2, 0.16, (1/63)/365, 0.1, 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For Sierra Leone\n",
    "paramsSL = [0.361, 0.7, 1.6, 1.7, 0.2, 0.1, 0.1, 0.2, 0.16, (1/63)/365, 0.1, 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#For both\n",
    "init = [1000000, 20000, 15, 10, 0, 0]\n",
    "times = np.linspace(0, 200, 200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_LIB = solve_ivp(fun=lambda t, y: LIB(t, y, paramsLIB), t_span=[min(times),max(times)], y0=init, t_eval=times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_LIB = pd.DataFrame({\"t\":sol_LIB[\"t\"],\"SL\":sol_LIB[\"y\"][0],\"SH\": sol_LIB[\"y\"][1],\"E\":sol_LIB[\"y\"][2],\"I\":sol_LIB[\"y\"][3], \"H\":sol_LIB[\"y\"][4], \"R\":sol_LIB[\"y\"][5]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol_SL = solve_ivp(fun=lambda t, y: SL(t, y, paramsSL), t_span=[min(times),max(times)], y0=init, t_eval=times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_SL = pd.DataFrame({\"t\":sol_SL[\"t\"],\"SL\":sol_SL[\"y\"][0],\"SH\": sol_SL[\"y\"][1],\"E\":sol_SL[\"y\"][2],\"I\":sol_SL[\"y\"][3], \"H\":sol_SL[\"y\"][4], \"R\":sol_SL[\"y\"][5]})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdAAAAENCAYAAAC//TTxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VFX+//FXGj3UgUDoSCRUASmKoIiAqCyu7ShYwPUnogvWVVGXou5i72VdFBBWLMdVdlmlKUWkK0VaQFpoCZBQQ0lCkvn9MQPfiAlkhszcSeb9fDzymLl37p375maYT+69554T4Xa7EREREd9EOh1ARESkJFIBFRER8YMKqIiIiB9UQEVERPygAioiIuIHFVARERE/qICKiIj4QQVURETEDyqgIiIifoh2OoAf1HWSiIh/IpwOUJqUxAJKSkqKX+u5XC7S09OLOc35Uy7fKJdvlMs3pTVXfHx8MaYR0ClcERERv6iAioiI+EEFVERExA8qoCIiIn5QARUREfGDCqiIiIgfVEBFREKc2+0mJy/H6RhyBhVQEZEQ9+22b+kwrgPzd813OorkowIqIhLCsnOzeeGnF0hKT2LbkW1Ox5F8VEBFRELYv5L+RfKRZC6sfiEDEgc4HUfyUQEVEQlRh7MO88aKNwD4+5V/JyYyxuFEkp8KqIhIiHrvl/c4mHWQzrU784eEPzgdR86gAioiEoJ2H93NR2s/AmBE5xFERGgglVCjAioiEoJe+uklsnKz6NekH+1qtXM6jhRABVREJMSsTV/L15u/JiYyhuEdhzsdRwqhAioiEkLcbjfPL30eN24GtRhEw8oNnY4khQjKgNrGmGbAF/lmNQFGApO88xsByYCx1h4MRiYRkVD0/Y7vWZCygCplqvBQu4ecjiNnEZQjUGvtRmttW2ttW+Bi4DgwBRgOzLbWJgCzvdMiImEpOzebZ5c8C8BD7R6iWrlqDieSs3HiFO5VwBZr7XbgemCid/5E4I8O5BERCQnj141n25FtNKnShLtb3u10HDkHJwrobcBn3udx1tpU7/M9QJwDeUREHJd+Ip03V7wJwKhLRlEmqozDieRcgnIN9BRjTBmgH/DUma9Za93GGHch6w0GBnuXw+Vy+bX96Ohov9cNJOXyjXL5Rrl841SukdNHknEyg95NenNru1t/d99nqO6vcBbUAgpcA6yw1u71Tu81xtSx1qYaY+oA+wpayVo7FhjrnXSnp6f7tXGXy4W/6waScvlGuXyjXL5xItfa/WsZv2o8URFRPN3+afbv31/sueLj488nohQg2Kdw+/N/p28BpgIDvc8HAv8Nch4REUe53W5GLx7tuW2l5SASqiU4HUmKKGgF1BhTEegFfJ1v9otAL2PMJqCnd1pEJGxMS57G4tTFVCtbjUfbP+p0HPFB0E7hWmuPATXOmLcfT6tcEZGwk5mTyfNLngfg8Q6PU7VsVYcTiS/UE5GIiEPGrhnLzqM7SayWyO2JtzsdR3ykAioi4oA9x/bwzqp3ABh96WiiI4PdplPOlwqoiIgDnl3yLMdzjtOnYR+61e3mdBzxgwqoiEiQzd81n6lbp1IuqhzPXvqs03HETyqgIiJBlJWbxdMLnwbgkfaPUC+2nsOJxF8qoCIiQfSPX/7BtiPbaFq1KYNbD3Y6jpwHFVARkSDZfmT76YZDYy4bo/5uSzgVUBGRIHC73fx10V/JzM3kxqY3cln8ZU5HkvOkAioiEgQzt89kzs45xMbEMqLzCKfjSDFQARURCbDjJ48zYpGnaD7R4QlqVajlcCIpDiqgIiIB9ubKN0k5lkKrGq24q8VdTseRYqICKiISQL8e/JV/rv4nEUTwQtcX1ONQKaICKiISIHnuPJ748Qly3DkMSBxA+1rtnY4kxUgFVEQkQD5J+oSf9v5EzfI1ebrT007HkWKmAioiEgCpx1IZs2wMAM93eV5DlZVCKqAiIgEwctFIMk5m0LNBT/o27ut0HAkAFVARkWI2I3kG05KnUTGmImMuG0NERITTkSQAVEBFRIrRkewjPLPwGQCGdxhO3Up1HU4kgaICKiJSjF5Y9gJ7ju+hXa12DGwx0Ok4EkAqoCIixWRx6mImJU0iOiKaV7q9QlRklNORJICCdkevMaYq8BHQCnADfwI2Al8AjYBkwFhrDwYrk4hIcTl+8jiP/fAYAMPaDaN59eYOJ5JAC+YR6FvADGttInARkAQMB2ZbaxOA2d5pEZES58WfXmR7xnaaV2/Og20fdDqOBEFQCqgxpgpwOTAOwFqbba09BFwPTPQuNhH4YzDyiIgUpyWpSxi3bhxREVG8ecWbGuczTATrFG5jIA2YYIy5CFgOPATEWWtTvcvsAeKClEdEpFgcP3mcx+Z7T922HUYrVyuHE0mwBKuARgPtgWHW2qXGmLc443SttdZtjHEXtLIxZjAw2LscLpfLvxDR0X6vG0jK5Rvl8o1y+cbXXI999xjJR5JpVbMVz/d6PmBHn6G6v8JZsAroLmCXtXapd/rfeAroXmNMHWttqjGmDrCvoJWttWOBsd5Jd3p6ul8hXC4X/q4bSMrlG+XyjXL5xpdcS1OX8t7P7xEVEcWrXV/lyMEjIZGrIPHx8cWYRiBI10CttXuAncaYZt5ZVwHrganAqRulBgL/DUYeEZHzdSLnBI/OfxQ3boa2HUprV2unI0mQBXNgumHAZGNMGWArcDeeAm6NMfcA2wETxDwiIn578acXST6STPPqzXm43cNOxxEHBK2AWmtXAR0KeOmqYGUQESkOi1IWMW6tp9Xt65e/rla3YUo9EYmI+OBw1mEemvcQbtw82O5B2tRs43QkcYgKqIiID55Z+Awpx1JoW7MtD7V7yOk44iAVUBGRIvrP5v8wZcsUykeX5+3ubxMTGeN0JHGQCqiISBHsPrqbpxY+BcDoS0ZzQdULHE4kTlMBFRE5hzx3Hg/Pe5gj2Ufo2aAntyfe7nQkCQEqoCIi5zB2zVgWpS7CVd7Fa5e/RkREhNORJASogIqInMW6/et46aeXAHi126u4yqs7PfFQARURKURmTibD5g4jOy+bOxLvoFfDXk5HkhCiAioiUojnlj7HxoMbaVy5MaMuGeV0HAkxKqAiIgX4dtu3TFw/kZjIGN7v8T4VYio4HUlCjAqoiMgZdhzZwV/m/wWAEZ1HqLchKZAKqIhIPidzT/LA3Ac4kn2EqxtezZ9a/snpSBKiVEBFRPIZ+cNIVu5bSXzFeN2yImelAioi4jVn5xxeX/o6URFRvN/jfaqVq+Z0JAlhKqAiIkDqsVQemufpHP7xDo/TsXZHhxNJqFMBFZGwl5uXy7C5wziQeYCrGl3Fny/6s9ORpARQARWRsPfq8ldZnLqYmuVrMqHfBCIj9NUo56ZPiYiEtVnbZ/H2qreJjIjk3SvfJa5inNORpIRQARWRsLX18FYenPsgAE91fIqudbs6nEhKEhVQEQlLx08eZ/D3g8k4mcE1ja7h/jb3Ox1JSpjoYG3IGJMMZAC5QI61toMxpjrwBdAISAaMtfZgsDKJSHhyu9088eMTJB1IokmVJrxxxRu631N8Fuwj0CuttW2ttR2808OB2dbaBGC2d1pEJKA+Xv8xU7ZMoUJ0BT7q+RGxZWKdjiQlkNOncK8HJnqfTwT+6GAWEQkDP+39idGLRwPw6uWv0qx6M2cDSYkVzALqBmYZY5YbYwZ758VZa1O9z/cAav4mIgGTdjyNId8PIcedw72t7uX6C653OpKUYEG7Bgp0tdbuNsbUAr4zxmzI/6K11m2McRe0orfgDvYuh8vl34jw0dHRfq8bSMrlG+XyjXJ5ZOdmc8v0W9hzfA9d63fljWvfICYqxvFcRRWqucJZhNtdYM0KKGPMaOAocC/Q3VqbaoypA8yz1p7rfIo7JSXFr+26XC7S09P9WjeQlMs3yuUb5fI0GnpywZNM3jCZOhXrMO2P06hVoZbjuXxxvrni4+MB1FKqGAXlFK4xpqIxJvbUc6A3sBaYCgz0LjYQ+G8w8ohIeJmYNJHJGyZTLqoc43qNK7R4ivgiWNdA44AFxphfgGXAt9baGcCLQC9jzCagp3daRKTYLEpZxKhFowBPo6GLal7kcCIpLYJyDdRauxX43afWWrsfuCoYGUQk/OzM2Ml9s+8jx53DA20e4IamNzgdSUqRIhVQY0w0MBNYZK0dEdhIIiLn72j2Ue6edTcHMg9wZb0rGd5Rt5lL8SrSKVxrbQ6QCNQObBwRkfOXm5fLA3MeON3T0Hs93iMqMsrpWFLK+HIN9Dngj8aYK4wxv2/7LSISIp5d8iyzd86matmqTLx6IlXKVnE6kpRCvlwD/QeezhDmABhjTs13W2uDeT+piEihJqybwLh14ygTWYbxvcbTpEoTpyNJKeVr4SvoHiLdVyQiIWH2jtmMXDwS8LS47Vyns8OJpDTzpYA2DlgKEZHztG7/Ou6fcz957jweaf8INyXc5HQkKeWKXECttdsBjDEXAk2ttdMClkpExAd7ju1h4MyBHDt5jBsuuIHH2j/mdCQJA0UuoN6xOy1wJeA2xiQAm4Ax1tqRAconInJWx04e4+5Zd5N6LJWOcR159fJXNbanBIUvp3BfAXoAWUAZa+02Y8wS4A+ACqiIBN3JvJMM/n4wq9NX0zC2IeN6jaNcdDmnY5V2XYDngdLetPkwMAJYVNgCvtzG0geYAXyQb956QE3cRCTo3G43j/3wGPN2zaNGuRp8cs0n1Chfw+lY4SAciid4/o3Pn20BXwpoeeDgGfNcQLaPoUREztuYZWP4avNXVIiuwKQ+k3S7SvCEQ/E85az/Vl8K6GqgL9AZwBjzKp7Tt7/4HU1ExA8frvmQ91e/T3RENB/2/JC2Nds6HUnCkC8F9K9AWeASPPd+PgrkAaOLP5aISMH+s/k/jF4yGoDXr3id7vW7O5pHICoq6uLExMQWCQkJLa+55pomGRkZxTrS19tvv13jrrvuanC2Zb755pvY7777ruKp6Zdffrnmu+++G9Bz+kX+R1prFwAd8PRINA14H+jonS8iEnDzd8/n4R8eBmBE5xG61zNElC1bNm/Dhg3rN23atC4mJsb92muv1Qx2hjlz5sT++OOPlU5NP/HEE2lDhw7dH8ht+vRXgrV2LfA48AjwhLV2dUBSiYic4Ze0X/h/3/0/T8vb1oMZ0maI05GkAF27dj26efPmsgCjR4+OS0hIaJmQkNDyueeeqwWwcePGMo0bN27Zr1+/xk2aNGnZp0+f00esdevWbZ2amhoNMH/+/AqdOnVqdub7f/rpp1XatGmT2Lx58xZdunS5cOfOndEbN24sM2nSpJoffPBBXGJiYosZM2ZUevTRR+NHjhwZB7Bo0aLyF110UeKFF17YolevXhekpaVFAXTq1KnZ/fffX7d169bNGzVq1GrGjBmVztze2RS5gBpjqhljvgQygA1AhjHmS+/9oSIiAbPhwAYGTB9wuqOEEZ01qmIoOnnyJDNnzqzcunXrEz/++GOFTz/9tMby5cuTfv7556RJkybVXLhwYXmA5OTkckOHDt23devWdbGxsXmvvPJKkY9Ye/XqdXTVqlUbkpKS1t98880HnnvuudrNmjXLvuuuu9KGDBmyd8OGDev79OlzNP86gwYNajxmzJhdv/766/qWLVueePLJJ+NPvZaTkxOxZs2apJdeemnnc889F//7LRbOlyPQccBNeK5/nvq5EfjIlw2KiPhi6+Gt3DbtNg5lHaJXg1680f0NIiOK9RKbnKesrKzIxMTEFq1bt25Rr1697Iceeih93rx5la699tpDlStXzqtSpUreddddd3Du3LmxALVr187u3bv3MYA777xz/6JFi4p85Ldt27Yy3bp1S7jwwgtbvP3227U3bNhQ/mzL79+/PyojIyPquuuuOwpw77337l+yZMnp7d1yyy0HAbp06XJs165dZXz5d/vyKewJbAPaAuWAdsB273wRkWK3++hubv32VtJOpNE1visfXPUBMZEaTTHUnLoGumHDhvUTJ07cWa5cOffZlj+zp6hT01FRUe68vDwATpw4UWB9Gjp0aIMHHnhg36+//rr+3Xff3Z6VlXVef02dyhodHU1ubq5PXVj5suEtwFxr7Wprbba19hdgLrDZlw2KiBTF3uN7Md8aUo6l0CGuA+N7j1cvQyXIlVdeeXTatGlVMzIyIo8cORI5bdq0aldeeWUGQGpqapnvv/++IsDkyZOrd+nS5ShAvXr1shcuXFgBwFpbraD3zcjIiGrQoMFJgI8//vh0K9vY2NjcjIyM342aXqNGjdzKlSvnnrq+OW7cuBqXXnrp0TOX88dZu/Izxlyeb3ISMMoYswrPNdDmeE7pji6OICIipxzIPED/af1JPpJMqxqtmHT1JCrGVDz3ihIyunbtenzAgAH727dv3xzgzjvvTLvssstObNy4sUyjRo0y33nnnVqDBw+ukJCQkPmXv/wlDWDkyJEpQ4YMafTcc8/ldunSJaOg933mmWdS+vfvf0GVKlVyunbtmrFjx46yADfddNOhm2+++YLp06dXffPNN3fkX2fChAnb7r///oYPPvhgZIMGDbI+++yz5OL4N0a43YUfaRtj8vAMol3o+kBeUQfUNsZEAT8Du621fY0xjYHPgRrAcuBOa+25ejZyp6SkFGVzv+NyuUhPT/dr3UBSLt8ol29KWq6M7Axu/fZWfkn/hYSqCXzV96ugdtFX0vZXUcXHx0PxjN/88/msvHHjxjJ9+/ZN2LRp07piyBIMHQp74VyncHec42c7sNOHIA8BSfmmXwLesNY2xdNN4D0+vJeIlDIZ2RkMmD6AX9J/oWFsQz6/9nP1bysh66xHjtbaRsW1IWNMPeA64O/Ao8aYCDyjuwzwLjIRz+ngfxTXNkWk5DhVPFfsW0G9SvX44rovqF2xttOxpJg1a9YsuwQdfZ6VL8OZAWCMqQH85mKEtXZHIYvn9ybwBBDrna4BHLLW5nindwF1fc0jIiXfmcXz333/Tf3Y+k7HEjkrXwbU7oXnXtAzi5z7XO9jjOkL7LPWLjfGdPc1pDFmMDAYwFqLy+Xy9S0ATzNlf9cNJOXyjXL5JtRzHck6wsDPB7Ji3woaVmnIrAGzaFS1keO5Qk2o5gpnvhyBfgDUK2B+US5KXwb0M8Zci+ce0srAW0BVY0y09yi0HrC7oJWttWOBsd5Jt78X0ktr44BAUS7fKJdvXC4X21K2/fa07TVfUCmnkqN5Q3l/FUMjIilGvhTQGsBM4BZrrU/30FhrnwKeAvAegf7FWnu7t2vAm/G0xB0I/NeX9xWRkutI1hGdtpUSzZeOFN4BGgB1vQ2AisOTeBoUbcZToMcV0/uKSAg7kHmAPp/2UfGUgKtQoUK7QL23L0egXwFDgfUAxphT891FvQ8UwFo7D5jnfb4V6ORDBhEp4fYd30f/af3ZcHADDWMb8sV1X6h4SonkSwGdDFQpYH5xHY2KSCm3K2MXt067leQjySTWSGTy1ZN1q4qUWL4U0AbAT3huRTkUmDgiUlptPbyVW7+9lZRjKbSq0YoZd8wg4rj+/g4XEc9GXByI93WPci8PxPsWhS8FdCye062Li9DdnojIaUkHkug/rT9pJ9LoENeBSVdPomaFmqQfD73WriJF5UsBvQpoBaQZY5KBXO98t7U2IH9ZiEjJtyptFbdPv51DWYfoGt+VCb0nUCGmgtOxJMicPFIMFF8KaBvvYyzQOt/8s477JiLha+7OuQz+fjDHc47Tu2Fv/tHjHxqSTEoNXwro3QFLISKlzr83/ZvHfniMHHcONzW9ideueE2DYUup4svtJxMDGURESge3280Hqz/gb8v+BsADbR7g6U5PExGhBkMSfMePH18ZqPf2pS/c8YW85LbWahgyESHPncezS57lo7UfATD6ktHc2/peh1OJBIYvp3AH4bneeerPyFPP3WgcT5Gwl5WbxcPzHmbq1qnERMbwVve3uP6C652OJRIwvhTQZ/M9j8LTkKgfUNiRqYiEiUNZh7j3u3tZlLqISjGV+KjXR3Sr283pWCIB5cs10GfPnGeMGYvG8BQJa9sOb2PgzIFsObyFmuVr8kmfT2jlauV0LAmcwxTcK11pdPhsL/pyDbTBGbMqA82Ai/wIJSKlwLI9y/jTrD9xMOsgzas3Z+LVE6lbSX9Tl3IjgOcp/UX0MJ5/a6F8OYW7rZD5a3x4DxEpJb7e/DWP/fAY2XnZ9Kjfg/d7vE9smVinY0ngLcLTsU7Y82U4s4gzfk7g2ZGDij+WiIQqt9vNa8tfY9jcYWTnZTOoxSAm9J6g4ilhx5droL4UWxEphU7knODx+Y8zZcsUIiMiGX3JaO5ppUb4Ep7OWUCNMbnnWMSn8UBFpGTalbGLe767h7X711IhugLv93ifXg17OR1LxDFFKXzqPkQkzC1KWcR9s+/jQOYBGsY2ZHzv8SRWT3Q6loijilJA250xXQkYChg8xXVVcYcSkdDgdrv5eP3HjFo8ilx3LlfUvYL3erxHtXLVnI4m4rhzFlBr7S8AxphywJ+Bx4FaeFrfPmut/TqgCUXEEZk5mTy98Gm++PULwNOn7fCOw4mKjHI4mUhoKMo10LLA/cATQByQBAyz1n4Z4Gwi4pDUY6nc+929rExbSbmocrx+xevqlk/kDEU5hbsVqI2nz1sLfAHkGWP6nVrAWjv1bG/gPXqdD5T1bvPf1tpRxpjGwOdADWA5cKe1Ntuff4iIFI+lqUu5b/Z9pJ1Io16leozrPY5WNdSzkMiZinJrSh3vYwSe655fAVPy/RTlFG4W0MNaexHQFuhjjLkEeAl4w1rbFDiIOqUXcUyeO493Vr3DLd/eQtqJNLrU6cL0G6areIoUoihHoDvwHH36zVrrBo56J2O8P26gBzDAO38iMBr4x/lsS0R8dyDzAA/Ne4g5O+cA8OeL/szjHR7XANgiZ1GURkSNimNDxpgoPKdpmwLvAVuAQ9baHO8iu1DH9CJB99Oen7h/zv2kHkulatmqvNX9LXo26Ol0LJGQF7QOEKy1uUBbY0xVPKd+i3wTmTFmMDDY+z64XC6/MkRHR/u9biApl2+UyzeF5cpz5/HG0jcYMW8Eue5cLql7Cf+6/l80qHLmuBHBzeU05ZKiinC7z+vsrF+MMSPx9KX7JFDbWptjjLkUGG2tvfocq7tTUlL82q7L5SI9Pd2vdQNJuXyjXL4pKNfBzIM8/MPDfL/jewCGtBnC8I7Dg3rKtiTtr1Bwvrni4+NBHeMUq6D0b2uMqek98sQYUx7ohed2mLnAzd7FBgL/DUYekXC2JHUJV0+5mu93fE/VslWZ0HsCIzqP0PVOER8Fq4P4OsBcY8xq4CfgO2vtN3iOQB81xmzGcyvLuCDlEQk72bnZvLDsBW7+5mZ2H91Nu5rtmHnDTHo37O10NJESKSjXQK21q/l9l4BYa7cCnYKRQSScbT60maFzh7ImfQ2REZE82PZBHmn/iI46Rc6DRlERKcXcbjcfLP+AJ2c/SWZuJvUr1eedK9+hY+2OTkcTKfFUQEVKqbTjaTw6/9HT93benHAzf+vyNw18LVJMVEBFSqGZyTN5/MfH2Z+5n2rlqjGmyxj6XdDv3CuKSJGpgIqUIgcyDzBq8Si+3uzpYbNrfFcm3jCRctnlHE4mUvqogIqUEjOSZzB8wXDSTqRRLqocwzsO555W91Crcq2QvK9RpKRTARUp4Q5kHmDEohH8Z8t/AOhcuzOvXf4ajas0djiZSOmmAipSQrndbqYlT+PphU+TfiKd8tHlebrj0wxqOYjIiGDd4i0SvlRARUqg3Ud389dFf2XW9lkAXFrnUl69/FUaVW7kbDCRMKICKlKC5OTlMH7deF75+RWO5xwnNiaW4Z2Gc1fzu3TUKRJkKqAiJcTqtNU8seAJ1qSvAaBv4748e+mz1K5Y2+FkIuFJBVQkxB3NPsory19h/Lrx5LnzqFupLn/v8nd6NezldDSRsKYCKhKi3G4305OnM3LxSFKPpRIVEcV9re/jsYsfo2JMRafjiYQ9FVCRELTp4CZGLB7Bj7t/BOAi10W83O1lWrlaOZxMRE5RARUJIRnZGby+4nXGrx1PjjuHqmWr8niHx7kz8U6iIqOcjici+aiAioSAPHce/970b8YsG0PaiTQiiOCOxDt4suOTVC9X3el4IlIAFVARhy1NXcrzS59nZdpKADrEdeBvXf5Ga1drh5OJyNmogIo4ZMuhLbzw0wtMT54OQK3ytXim8zPc1PQmIiIiHE4nIueiAioSZAcyD/DGijeYtH4SOe4cykeX54E2D3Bfm/vUulakBFEBFQmSzJxMJqybwNur3uZI9hEiiKB/s/785eK/qDMEkRJIBVQkwHLzcvlq81e8vvx1dh7dCcAVda/gr53/SosaLRxOJyL+UgEVCZA8dx7Ttk3j1eWvsunQJgASqyUyovMIutfv7mw4ETlvQSmgxpj6wCQgDnADY621bxljqgNfAI2AZMBYaw8GI5NIoLjdbubumsvLP798ut/aBrENeLT9o9zY9EbdzylSSgRr+IYc4DFrbQvgEuDPxpgWwHBgtrU2AZjtnRYpsZakLuGG/93AnTPuZE36GuIqxDHmsjH8cMsP3HLhLSqeIqVIUI5ArbWpQKr3eYYxJgmoC1wPdPcuNhGYBzwZjEwixWnZnmW8+927zE6eDUC1stUY2nYoA1sMpHx0eYfTiUggBP0aqDGmEdAOWArEeYsrwB48p3gLWmcwMBjAWovL5fJr29HR0X6vG0jK5ZtQyeV2u5m7fS4vLHyB+TvmAxBbJpZHOj/CsI7DqFy2ssMJPUJlf51JuXwTqrnCWYTb7Q7axowxlYAfgL9ba782xhyy1lbN9/pBa221c7yNOyUlxa/tu1wu0tPT/Vo3kJTLN07ncrvdzNk5h7dWvsXyfcsBqFymMsM6DqN/k/5UK3euj3BwOb2/CqNcvjnfXPHx8QDqoaMYBe0I1BgTA3wFTLbWfu2dvdcYU8dam2qMqQPsC1YeEV/lufOYtX0Wb618i9XpqwHPqdrBrQczqOUgmsQ3CckvXhEJjGC1wo0AxgFJ1trX8700FRgIvOh9/G8w8oj4Iis3iymbp/DB6g9O345Ss3wTl7dEAAAQE0lEQVRNhrQZwp3N71TvQSJhKlhHoJcBdwJrjDGrvPOexlM4rTHmHmA7YIKUR+ScDmUd4pOkTxi/bjx7j+8FoE7FOtzf5n4GJA5Q4yCRMBesVrgLKPzc+1XByCBSVLsydvHh2g/5bONnHDt5DIDm1ZszpM0Q+jXpR5moMg4nFJFQoJ6IRLxW7FvB+LXjmbp1KrnuXAC61e3GkNZDuKLeFRohRUR+QwVUwlpWbhZTt0zl4/UfsyrNc3UhKiKKGy64gSFthtDK1crhhCISqlRAJSylHE3hX0n/YvKGyezP3A9A1bJV6d+sP4NaDKJebD2HE4pIqFMBlbDhdrtZsmcJE9ZNYEbyjNOnaVtUb8GfWv6JPzb9oxoGiUiRqYBKqXc46zBTtkzhk6RPSDqQBEB0RDT9mvTj7pZ30zGuo65viojPVEClVHK73Szbs4zJGybz7bZvyczNBDz3b97R/A7uSLxDg1iLyHlRAZVSJf1EOl/++iWfbvyUrYe3np7frW43+jfrT59GfSgbVdbBhCJSWqiASomXm5fLj7t/ZPLGycxKnkWOOweAuApx3HrhrdzW7DYaVm7ocEoRKW1UQKXESjqQxFebvmLKlinsObYHgMiISHo16MWAxAH0qN+D6Eh9xEUkMPTtIiVK6rFUJm2ZxKRVk043CAJoGNuQW5vdirnQUKdiHQcTiki4UAGVkHc0+yjTk6fz1eavWLB7AW48Q/BVLVuVPzT5Azc1vYkOcR3UklZEgkoFVELSiZwTzN05l2+2fcOs7bM4kXMCgDKRZbg24Vr61u9LjwY91CBIRByjAioh40TOCebtnMc3277hux3fne7IHaBTXCduTLiRvo37klAvQeNuiojjVEDFUZk5mczbNY//bf3f74rmRa6L6NukL30b96VB5QYOphQR+T0VUAm6jOwM5uycw6zts1Q0RaTEUgGVoEg5msKsHbOYlTyLRamLOJl38vRrp4rmdY2v0/2aIlJiqIBKQLjdbpIOJDFz+0xmbZ/F6vTVp1+LjIjkktqX0Lthb/o06qOiKSIlkgqoFJtjJ4+xMGUhc3fOZe7Ouew8uvP0a+Wjy9O9Xnd6N+xNzwY9qV6uuoNJRUTOnwqo+M3tdvPrwV+Zu8tTMJftWUZ2Xvbp113lXfRu0JveDXvTtW5XDRUmIqWKCqj4JCM7gwW7F5wuminHUk6/FkEE7Wu1p0f9HnSv1502rjZERUY5mFZEJHCCUkCNMeOBvsA+a20r77zqwBdAIyAZMNbag8HII0WXnZvNyn0rWZiykAUpC1i+d/npztrBc5TZvV53etTvQbe63XRqVkTCRrCOQD8G3gUm5Zs3HJhtrX3RGDPcO/1kkPJIIXLzclm3fx0LUxayNG0pC3YuON0LEEBURBSd4jrRvb6naLas0ZLIiEgHE4uIOCMoBdRaO98Y0+iM2dcD3b3PJwLzCFABPZR1iDHLxrD92HY+v/pz9ZmaT25eLhsObmBp6lIWpS5icepiDmUd+s0yzao147L4y+ga35VL6lxClbJVHEorIhI6nLwGGmetTfU+3wPEFbagMWYwMBjAWovL5fJpQ1XzqjI9eToHMg+QHpFOc1dzfzMHRHR0tM//Jn9l5WSxfM9yFu5cyIKdC1i8azGHsw7/ZpmGVRrSo1EPrmpyFd3qdaN2pdpByVZUwdxfvlAu3yiXb0I1VzgLiUZE1lq3McZ9ltfHAmO9k25/+kG9qv5VfLnpSz5f9TnD2g7zM2lguFyugPXtejT7KD/v/Zmle5aybM8yVqWtIjM38zfL1KtUj061O3FpnUvpGt/1dA9Ap3KlZ4ZWv7OB3F/nQ7l8o1y+Od9c8fHxxZhGwNkCutcYU8dam2qMqQPsC+TG+jTqw5ebvmRm8syQK6DFJc+dx9bDW1mxb8Xpn6QDSeS5836zXLNqzehUuxOda3emU+1O1K1U16HEIiIll5MFdCowEHjR+/jfQG7sinpXUD66PCvTVpJ6LLVUDLp8MPMgK9NWsnLfSlbsW8HKfSs5nP3b07HREdFcVOsiOtfuTOfanekQ10EtZUVEikGwbmP5DE+DIZcxZhcwCk/htMaYe4DtgAlkhvLR5enZuCf/2/Q/ZibPZFDLQYHcXLE7fvI46/avY3X6an5J+4WVaSvZenjr75arXaE27Wu1P/3T2tWaCjEVHEgsIlK6BasVbv9CXroqGNs/5abEm/jfpv8xbt047mh+B9GRIXEJ+HeOZh89XSxXp69mTfoaNh/ajJvfXiYuF1WONq42tKvV7nTBjK+k6xwiIsEQmhUkQG5ufjOj5o1i6+GtTNk8hVsuvMXpSKQdT2PVtlUs2baEtelrWZ2+mq2Ht/6uWEZHRNOsejNa12hN65qtaV+zPc1rNCcmMsah5CIi4S2sCmhMVAwPt3+YR354hDdWvEHfJn2D1j9rZk4mmw5tYv2B9STtTyLpQBIbDm4g/cTvW9XFRMaQWD2RNq42tHa1prWrNYnVEikXXS4oWUVE5NzCqoAC3Nj0Rt7/5X02HdrEg/Me5J9X/bNYe9LJzctlR8YONh3aRNIBb6E8sIEth7f8rjUsQGxMLG3i2nBB7AW0crWidY3WNKvejLJRZYstk4iIFL+wK6DRkdF82PND+k3tx7Rt03hw7oP8/bK/+9y7TmZOJlsPb2XToU1sPrT59OPWw1vJys363fKREZEkVE2gefXmNK/enMTqibSo3oK6lepSs2bNkLzvTEREChd2BRQgoVoC/+z5T+6eeTdTtkxh/u753JRwE5fFX0ajyo2oFFOJ7NxsMk5mkJGdwZ5je9h5dCc7M3ayK2MX2zO2syNjR4FHlAC1K9YmoWoCidUSaV6jOS2qt6Bp1aYazktEpBQJywIKcHndy5l540we+eERVuxbwdg1Yxm7Zuy5V/SKjIikceXGJFRLIKFqAk2rNj39GFsmNoDJRUQkFIRtAQVoWrUpU/tNZcW+FXyz7RvW7V/HroxdnMg5QZmoMsSWiaViTEXiKsRRP7Y+9SvVp15sPepXqk+jKo10nVJEJIyFdQEFiIiI4OK4i7k47mKno4iISAmigRxFRET8oAIqIiLiBxVQERERP6iAioiI+EEFVERExA8qoCIiIn5QARUREfGDCqiIiIgfItxu97mXCi0lLrCISIiIcDpAaVISj0Aj/P0xxiw/n/UD9aNcyqVcofNTynNJMSqJBVRERMRxKqAiIiJ+CLcCWvTxyoJLuXyjXL5RLt8olxRJSWxEJCIi4rhwOwIVEREpFmEzHqgxpg/wFhAFfGStfdGBDPWBSUAcnttxxlpr3zLGjAbuBdK8iz5trZ3mQL5kIAPIBXKstR2MMdWBL4BGQDJgrLUHg5ipmXf7pzQBRgJVCfI+M8aMB/oC+6y1rbzzCtw/xpgIPJ+3a4HjwCBr7Yog5noF+AOQDWwB7rbWHjLGNAKSgI3e1ZdYa4cEMddoCvm9GWOeAu7B8/l70Fo7M4i5vgCaeRepChyy1rYN8v4q7PvB8c+YFCwsCqgxJgp4D+gF7AJ+MsZMtdauD3KUHOAxa+0KY0wssNwY8533tTesta8GOU9BrrTWpuebHg7Mtta+aIwZ7p1+MlhhrLUbgbZw+ve4G5gC3E3w99nHwLt4vuROKWz/XAMkeH86A//wPgYr13fAU9baHGPMS8BT/N/vbYu1tm2AspwrFxTwezPGtABuA1oC8cD3xpgLrbW5wchlrb01X5bXgMP5lg/W/irs+2EQzn/GpADhcgq3E7DZWrvVWpsNfA5cH+wQ1trUU38hWmsz8PxlWzfYOXx0PTDR+3wi8EcHs1yF58tsuxMbt9bOBw6cMbuw/XM9MMla67bWLgGqGmPqBCuXtXaWtTbHO7kEqBeIbfua6yyuBz631mZZa7cBm/H8vw1qLu9RnQE+C8S2z+Ys3w+Of8akYGFxBIrnQ7gz3/QuHP5LzXtqqB2wFLgMGGqMuQv4Gc9foUE7TZqPG5hljHED/7TWjgXirLWp3tf34Dm95JTb+O0XWyjss8L2T0GfubpAKsH3J357GryxMWYlcAT4q7X2xyDnKej3VhdPoT/l1P4Ktm7AXmvtpnzzgr6/zvh+KAmfsbAULkegIcUYUwn4CnjYWnsEz6mXC/CcqkwFXnMoWldrbXs8p4b+bIy5PP+L1lo3DnWlaIwpA/QDvvTOCpV9dpqT+6cwxphn8JwanOydlQo0sNa2Ax4FPjXGVA5ipJD7vZ2hP7/9Iy3o+6uA74fTQvEzFs7CpYDuBurnm67nnRd0xpgYPP85Jltrvwaw1u611uZaa/OADwnQqatzsdbu9j7uw3OdsROw99RpIe/jPiey4SnqK6y1e70ZQ2KfUfj+cfwzZ4wZhKexzO3eL168p0j3e58vx9PA6MJgZTrL7y0U9lc0cCP5jtaDvb8K+n4ghD9j4S5cCuhPQIIxprH3SOY2YGqwQ3ivr4wDkqy1r+ebn/+6xQ3AWgeyVfQ2XMAYUxHo7c0xFRjoXWwg8N9gZ/P6zZFBKOwzr8L2z1TgLmNMhDHmEuBwvtNwAedtdf4E0M9aezzf/JrexlgYY5rgaYCyNYi5Cvu9TQVuM8aUNcY09uZaFqxcXj2BDdbaXadmBHN/Ffb9QIh+xiRMroF6WyIOBWbiuY1lvLV2nQNRLgPuBNYYY1Z55z0N9DfGtMVzaiYZuM+BbHHAFGMMeD4Xn1prZxhjfgKsMeYeYDueBhZB5S3ovfjtfnk52PvMGPMZ0B1wGWN2AaOAFyl4/0zDc3vBZjy3GNwd5FxPAWWB77y/01O3X1wOPGeMOQnkAUOstUVt6FMcuboX9Huz1q4zxlhgPZ5Tzn8OUAvcAnNZa8fx+2vsEMT9ReHfD45/xqRg6olIRETED+FyCldERKRYqYCKiIj4QQVURETEDyqgIiIiflABFRER8UNY3MYiUhyMMR/juQ/vWWvtaGfTiIjTVEBFvLzDuTU8yyK3A4f4bZ+tIhKmVEBF/s94oLr3+f1AGTzdqp3qmWaZtfZTJ4KJSOhRRwoiBTDGHAKq4BkfdZ533sfkO4XrHRx6FDAD2I+nH9W1wAA83egNwNPbzl2nhqkyxjQAXsYz6kdFPN3VPWqtdaorQhHxkxoRiZy/q4EawA6gI55hutoBq/EMEP02gDGmAjAHT1dsq/H0ZdodmGOMcQU9tYicFxVQkfO3BU+fpC96pyvg6Zj8VP+87byP1+EZyisF2IhnUOcdQE3g5mCFFZHioQIqcv42eIcLO+Sd3mutPQxkeKcreB8beR/rAg95fy7wzmsahJwiUozUiEjk/J05akhho4gkex+XAx1PjdFpjKkKRAQmmogEigqoSPBMA7YBFwMLjTGrgQZ4roNeC8xzLJmI+EyncEWCxFp7DOiBZ8zJBnha9DYDPsFzTVREShDdxiIiIuIHHYGKiIj4QQVURETEDyqgIiIiflABFRER8YMKqIiIiB9UQEVERPygAioiIuIHFVARERE/qICKiIj44f8DjDDk1b6l/hAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#sLline = plt.plot(\"t\",\"SL\",\"\",data=output_SL,color=\"red\",linewidth=2)\n",
    "#sHline = plt.plot(\"t\", \"SH\", \"\", data = output_SL, color = \"red\", linewidth = 2)\n",
    "iline = plt.plot(\"t\",\"I\",\"\",data=output_SL,color=\"green\",linewidth=2)\n",
    "#rline = plt.plot(\"t\",\"R\",\"\",data=output_SL,color=\"blue\",linewidth=2)\n",
    "plt.xlabel(\"Time\",fontweight=\"bold\")\n",
    "plt.ylabel(\"Number\",fontweight=\"bold\")\n",
    "legend = plt.legend(title=\"Population\",loc=5,bbox_to_anchor=(1.25,0.5))\n",
    "frame = legend.get_frame()\n",
    "frame.set_facecolor(\"white\")\n",
    "frame.set_linewidth(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAENCAYAAACy40S2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOWhxvFfFghbQoAgkEAIhk0WBVREQIsoFsGqV/EVtKIWNyp1Qeq+LxTFBZRai4JCxeWl2ityAUHUglgUAQFZwhoCIRACCSSQfeb+MQMNGCATZuZMkuf7+eSTzJlzZh5OhnlyzpzznjC3242IiIj4R7jTAURERKoTFauIiIgfqVhFRET8SMUqIiLiRypWERERP1KxioiI+JGKVURExI9UrCIiIn6kYhUREfGjSKcD+JGGkBIRqZwwpwNUJ9WpWNm1a1ellouLiyMrK8vPaU5fqOaC0M2mXL5RLt9Ux1zx8fF+TiPaFSwiIuJHQdliNcZMBa4EMq21XbzTxgO/A4qALcBt1toc732PAiOAUuBea+2XwcgpIiJyuoK1xfo+MPC4aQuALtbas4GNwKMAxphOwFCgs3eZt4wxEUHKKSIiclqCUqzW2kXA/uOmzbfWlnhvLgVaen++GvjYWltord0GbAZ6BiOniIjI6QqVz1j/AMz1/pwA7Chz307vNBERkZDn+FHBxpjHgRJgRiWWvRO4E8BaS1xcXKUyREZGVnrZQArVXBC62ZTLN8rlG+WSinC0WI0xt+I5qOlSa+2R81DTgVZlZmvpnfYr1trJwGTvTXdlDzevjofQB1qoZlMu3yiXb0IxV6mrlGZnNNPpNiHEsV3BxpiBwEPAVdbaw2XumgUMNcZEGWPaAO2AH53IKCISyvKK8uhr+/Lkt09S7Cp2Oo54BaVYjTEfAf8BOhhjdhpjRgCTgGhggTHmZ2PM2wDW2rWABdYB84B7rLWlwcgpIlKVfLDhA9Jy0/hux3dEhjn+yZ54BeU3Ya0dVs7kKSeZ/0XgxcAlEhGp2gpKCvj76r8D8NCFDxEWplEJQ0WoHBUsIiI+sBstmfmZdGrciYHJxw8TIE5SsYqIVDHFrmLeWvUWAKO6jdLWaohRsYqIVDGfbvqUHXk7aBvblivbXOl0HDmOilVEpAopdhXzxso3ALiv+31EhGvE11CjYhURqUL+tflfbM/dTpuYNlx15lVOx5FyqFhFRKqIElfJ0a3Ve7vfS2S4TrEJRSpWEZEqYtbWWWw7uI3W0a25tu21TseRE1CxiohUAaWuUiasmABoazXUqVhFRKqA2dtms+XAFlo2aMl17a5zOo6chIpVRCTEudyuo1urf+r2J2qF13I4kZyMilVEJMTN3jqbjTkbia8fj2lvnI4jp6BiFREJYaWuUl5d8Srg2VqtHVHb4URyKipWEZEQ9q8t/2JzzmZaNWjF0A5DnY4jFaBiFREJUcWuYl5b/hoAD5z7gLZWqwgVq4hIiJq5cSbbc7dzZsMzua6tjgSuKlSsIiIhqLC0kNdXvA7Agz0e1HmrVYiKVUQkBH244UN2HdpFx0YduSpZYwJXJSpWEZEQk1+Sf3RM4DHnjiE8TG/VVYl+WyIiIWbaumlk5mfSNa4rA5MGOh1HfKRiFREJIQeLDjLp50kA/PncPxMWFuZwIvGVilVEJIS8teotsguz6dmsJ/1b9Xc6jlSCilVEJERkHMrgnTXvAPD4BY9ra7WKUrGKiISIV5e/SkFpAYPaDOK8Zuc5HUcqScUqIhICUvan8MnGT4gIi+CR8x5xOo6cBhWriEgIGLtsLC63i9+f9XuSY5OdjiOnQcUqIuKwpRlL+SrtK+rXqs8D3R9wOo6cJhWriIiD3G43L/zwAgB3d72bpvWaOpxITldQBp80xkwFrgQyrbVdvNMaA58ASUAqYKy12caYMGAiMAg4DNxqrV0RjJwiIsE2e9tsVu5dSdO6Tbnr7LucjiN+EKwt1veB44cPeQRYaK1tByz03ga4Amjn/boT+FuQMoqIBFWxq5hxy8YBMLrHaOrXqu9wIvGHoBSrtXYRsP+4yVcD07w/TwOuKTN9urXWba1dCsQaY1oEI6eISDDNWD+D1IOpJDdMZljHYU7HET9x8jpEzay1Gd6fdwPNvD8nADvKzLfTOy2D4xhj7sSzVYu1lri4uEoFiYyMrPSygRSquSB0symXb5TLN/7MdaDgABN+ngDAXy79Cy3OqPz2Q6iur5oqJC7wZ611G2PclVhuMjDZe9OdlZVVqeePi4ujsssGUqjmgtDNply+US7f+DPX8z88z97Dezm/2fn0btz7tB73dHLFx8dX+nmlfE4eFbznyC5e7/dM7/R0oFWZ+Vp6p4mIVAtbD2xlyi9TAHj2wmc1dGE14+QW6yzgFmCc9/vnZaaPMsZ8DFwAHCizy1hEpMp7/ofnKXYVc0P7Gzin6TlOxxE/C9bpNh8B/YA4Y8xO4Gk8hWqNMSOA7YDxzj4Hz6k2m/GcbnNbMDKKiATDop2LmL99PvVr1eeR8zV0YXUUlGK11p7ocLdLy5nXDdwT2EQiIsFX4irhmaXPAHBvt3s5o94ZzgaSgNDISyIiQfLB+g9IyU4hMTqR27vc7nQcCRAVq4hIEGQXZDN++XgAnrrgKepE1nE4kQSKilVEJAheW/EaOYU59G7Rm4FJxw9EJ9WJilVEJMB+2fcL7697n4iwCJ1eUwOoWEVEAsjldvHEkidwuV3c2vlWOjXp5HQkCTAVq4hIAP1z0z9ZtmcZTes2Zcy5Y5yOI0GgYhURCZADhQeOXmv1iQueIKZ2jMOJJBhUrCIiAfLK8lfYV7CPns16cl3b65yOI0GiYhURCYCyByy92OdFHbBUg6hYRUT8zOV28fiSx3XAUg2lYhUR8bOPUz7mpz0/6YClGkrFKiLiR5mHM48esPRMr2d0wFINpGIVEfGjp/7zFAeKDtC/VX+uTr7a6TjiABWriIifLNi+gC+2fkHdyLqM7TNWByzVUCpWERE/yCvK47EljwHw0HkP0Sq6lcOJxCkqVhERP3j5p5fZdWgX58Sdw4jOI5yOIw5SsYqInKYVmSuYunYqEWERvHzxy0SERzgdSRykYhUROQ3FrmIeWvwQbtzc1fUuujTp4nQkcZiKVUTkNLy9+m3W719P6+jWjD53tNNxJASoWEVEKmnrga28vuJ1AMb1HUfdyLoOJ5JQoGIVEamEUlcpo/89msLSQoa0G8LFLS92OpKECBWriEglvPvLuyzbs4xm9ZrxTK9nnI4jIUTFKiLio805m3n5p5cBeKnvSzSq08jhRBJKVKwiIj4odZXywL8foKC0gOvbXc+A1gOcjiQhRsUqIuKDyWsmsyJzBc3rNefZC591Oo6EIBWriEgFrc9az/jl4wEYf/F4GkY1dDiRhKJIpwMYYx4AbgfcwBrgNqAF8DHQBFgO3GytLXIspIjUeMWuYu6YcweFpYUMbT+U/q36Ox1JQpSjW6zGmATgXuA8a20XIAIYCrwEvG6tbQtkAxp4U0Qc9cbKN1iWsYz4+vE8feHTTseREBYKu4IjgbrGmEigHpAB9Af+6b1/GnCNQ9lERFiRuYKJKycSRhgT+k3QxcvlpBzdFWytTTfGvAKkAfnAfDy7fnOstSXe2XYCCeUtb4y5E7jT+1jExcVVKkdkZGSllw2kUM0FoZtNuXyjXKeWV5TH/f+8n1J3KQ9e+CBXnx16Fy8PpfUlDherMaYRcDXQBsgBZgIDK7q8tXYyMNl7052VlVWpHHFxcVR22UAK1VwQutmUyzfKdWoPLX6ILdlbOKvxWTzd9+mQyVXW6ayv+Ph4P6cRp3cFXwZss9butdYWA58BfYBY765hgJZAulMBRaTmmr99PjM2zKB2eG3evORNoiKjnI4kVYDTRwWnAb2MMfXw7Aq+FPgJ+AYYgufI4FuAzx1LKCI1UlZ+FmMWjQHg0Z6PclbjsxxOJFWFo1us1tof8ByktALPqTbheHbtPgyMNsZsxnPKzRTHQopIjeNyu7j/2/vZV7CPvvF9ub3L7U5HkirE6S1WrLVPA8cfu74V6OlAHBER3l79Nt/s/IZGUY14/TevEx7m9KdmUpXo1SIiUsayPcsYt2wcABP6TSC+gQ7uEd+oWEVEvLILsrnn63sodZdy99l3c1niZU5HkiqoQruCvUfofgl8b619MrCRRESCz+12M2bRGNLz0unetDsPn/ew05GkiqrQFqt3sIaOQPPAxhERccZ7a99j3vZ5xNSO4a3+b1E7orbTkaSK8mVX8HPANcaY3xhjagUqkIhIsK3eu5rnf3gegFcufoXEmESHE0lV5stRwX/DcwWarwGMMUemu621jh9dLCJSGblFuYz8eiRFriJu6XQLg9sMdjqSVHG+FmJYBaeJiIS8I5+rph5MpVPjTjx1wVNOR5JqwJdibROwFCIiDvj7mr8ze9tsGtRqwNuXvk2dyDpOR5JqoMLFaq3dDmCMaQ+0tdbOCVgqEZEA+y79O1788UUAJvabSHJsssOJpLqocLEaYxoDFrgEcBtj2gGbgLHWWu0/EZEqIz0vnZFfj8TldnFvt3sZmFThi2qJnJIvu4LH47kAeSFQ21q7zRizFPgdoGIVkSqhoKSAOxbcwf6C/fRr2Y8x545xOlJ10Rt4HmjodJAgOAA8CXxf3p2+nG4zEJgHvF1m2jrgzEpHExEJsie+f4JVWato1aAVb17yJhHhEU5Hqi5qSqmC59/5/Inu9KVY6wLZx02LA4oqEUpEJOhmbJjBRykfUSeiDu8OeJfGdRo7Hak6qSmlesQJ/72+FOtq4ErgAgBjzCt4dgOvOq1oIiJBsCJzBU8seQKAly56iS5xXRxOJNWVL8X6BBAF9MJz7upowAU84/9YIiL+k3Eog9sX3E6Rq4jbOt3GkHZDnI5UI0RERJzbsWPHTu3atet8xRVXnJmbm+vXC7+88cYbTYYPH37SYbJmz54dvWDBgvpHbr/88stNJ02a1MSfOY5X4X+ktfY74Dw8IzDNAd4CzvdOFxEJSfkl+YyYP4I9h/fQq3kvnuqlYy2DJSoqyrVhw4Z1mzZtWlurVi33q6++2jTYGb7++uvoxYsXNzhy+6GHHto7atSofYF8Tp/+erDW/gL8GXgAeMhauzogqURE/MDtdjP636NZlbWKxOhE3hnwjgbXd0jfvn3zNm/eHAXwzDPPNGvXrl3ndu3adX7uuefOAEhJSandpk2bzldddVWbM888s/PAgQOPbuEmJCR0zcjIiARYtGhRvZ49e3Y4/vE//PDDhmeffXbHs846q1Pv3r3b79ixIzIlJaX29OnTm7799tvNOnbs2GnevHkNRo8eHf/UU081A/j+++/rnnPOOR3bt2/facCAAcl79+6NAOjZs2eHkSNHJnTt2vWspKSkLvPmzWtw/POdTIWL1RjTyBgzE8gFNgC5xpiZ3vNbRURCzoSVE5i1dRYNajXg/cvf18FKDikuLubLL7+M6dq1a/7ixYvrffjhh02WL1++/qefflo/ffr0pkuWLKkLkJqaWmfUqFGZW7duXRsdHe0aP358hbdwBwwYkPfzzz9vWL9+/bohQ4bsf+6555p36NChaPjw4XvvvvvuPRs2bFg3cODAvLLL3HrrrW3Gjh27c+PGjes6d+6c//DDDx+9qn1JSUnYmjVr1r/00ks7nnvuOZ+udu/LFusU4Do8n68e+boWeNeXJxQRCYbZW2fzyvJXCCOMv/b/Kx0a/2ojRwKssLAwvGPHjp26du3aqWXLlkX33Xdf1rfffttg0KBBOTExMa6GDRu6Bg8enP3NN99EAzRv3rzo8ssvPwRw88037/v+++8rvKW4bdu22hdddFG79u3bd3rjjTeab9iwoe7J5t+3b19Ebm5uxODBg/MA7rjjjn1Lly49+nzXX399NkDv3r0P7dy506fdHL4U62XANqAbUAfoDmz3ThcRCRmr9q7ivm/vA+CJC57gskS9TTnhyGesGzZsWDdt2rQdderUcZ9s/rCwsHJvR0REuF0uFwD5+fnl9taoUaMS//jHP2Zu3Lhx3aRJk7YXFhae1oFSR7JGRkZSWlrq08VmfHniLcA31trV1toia+0q4Btgsy9PKCISSGkH0xj+5XAKSgu4of0N3NX1LqcjSRmXXHJJ3pw5c2Jzc3PDDx48GD5nzpxGl1xySS5ARkZG7a+++qo+wIwZMxr37t07D6Bly5ZFS5YsqQdgrW1U3uPm5uZGJCYmFgO8//77R4/6jY6OLs3Nzf3VKCBNmjQpjYmJKT3y+emUKVOaXHjhhXnHz1cZJx3S0BhzcZmb04GnjTE/4/mM9Sw8u4af8UcQEZHTlV2Qze/n/Z6s/CwuSriIcX3H/WorSJzVt2/fwzfeeOO+Hj16nAVw88037+3Tp09+SkpK7aSkpII333zzjDvvvLNeu3btCsaMGbMX4Kmnntp19913Jz333HOlvXv3zi3vcR9//PFdw4YNS27YsGFJ3759c9PS0qIArrvuupwhQ4Ykz507N3bChAlpZZd57733to0cObL1vffeG56YmFj40Ucfpfrj3xjmdp94y9wY48JzcfMTLg+4QuRC5+5du3ZVasG4uDiysrL8HOf0hWouCN1syuWb6pSroKSAYXOG8eOeHzmr8Vl89rvPiKkd43iuYDidXPHx8eCf62r/dDoLp6Sk1L7yyivbbdq0aa0fsgTLeeVNPFUhpnHyYhURcZzL7eL+f9/Pj3t+pHn95kz/7XS/l6pIRZ20WK21SUHKISJSaS/++CJfbP2C6FrR/OO3/yC+gU9nR0gI6NChQ1EV21o9IZ934RpjmgD1y06z1qadYHYRkYB6b+17vL36bSLDInlnwDt0atLJ6UhSw/lyofMBeM5lTTjuLrcvj1PO48biORe2i/ex/gCkAJ8ASUAqYKy1x19ZR0RquDnb5vDk908C8MrFr3BRwkUOJxLx7XSbt4GWHDtARJiPj1GeicA8a21H4BxgPfAIsNBa2w5Y6L0tInLUovRF3PP1PbhxM+bcMVzf/nqnI4kAvm1pNgG+BK631vrlXB9jTEPgYuBWAGttEVBkjLka6OedbRrwLfCwP55TRKq+5XuWM2L+CIpcRYzoPIL7u9/vdCSRo3wp1jfxDGGYYIzZaK31x9HCbYC9wHvGmHOA5cB9QDNrbYZ3nt1AMz88l4hUA+v3r2f4l8M5XHKYIe2G8MyFz+hcVfFZvXr1uh8+fHhlIB7bl2L9FBgFrAMwxhyZ7j6N81gjgR7An6y1PxhjJnLcbl9rrdsYU26JG2PuBO70zkdcXFzlQkRGVnrZQArVXBC62ZTLN1Ut15bsLfx+3u/JKczhd+1+x7RrpxEZHrzT6Kva+hJn+PKKnAE0LGf66fypuBPYaa39wXv7n3iKdY8xpoW1NsMY0wLILG9ha+1kYLL3pruyJ0hXx5O+Ay1UsymXb6pSrt2HdvM/X/wPuw/tpk98Hyb0nUDO/hzHc4UCPwwQIX7kS7EmAsuAhwC/vJqttbuNMTuMMR2stSnApXi2iNcBtwDjvN8/98fziUjVtL9gPzfOvZG03DS6Ne3G1AFTqRNZx+lY4gdhz4adG4jHdT/tXh6Ix60IX4p1MtAT+I/3ICN/+RMwwxhTG9gK3IbnSGNrjBmB5wo65iTLi0g1ll2QzY1zbyQlO4X2se35x8B/0KC2T9edFgkqX4r1Ujznmu41xqQCpd7pbmttpf/isNb+TPnjLV5a2ccUkeohpzCHYXOHsSZrDUkxSXw46ENdrLyacXLLMlB8Kdazvd+jga5lpmssYRHxu5zCHIbN+W+pzhw8kxb1WzgdS+SUfCnW2wKWQkSkjOz8bG6ccyOrs1YfLVWN/ytVRYWL1Vo7LZBBREQADhQe4OYvbmZV1ipaR7fGDrYqVfG7QJ3DCr6NFTz1BHe5rbUj/JRHRGqwA4UHGDZn2NFSnXnlTBIaHD88uUho82VX8K14Pk89ct7qkZ/dgIpVRE7L/oL93DT3Js/u39gk7BVWpSpVki/F+myZnyPwHMB0FXCiLVkRkQrJPJzJ0DlDSclOISkmiQU3LqBecT2nY4lUii+fsT57/DRjzGR+fRk5EZEKS89Lx/yfIfVgKu1j2/PRoI9IbJgYkiMcyUkdoPzR+aqrAye6w5fPWBOPmxQDdMBzqTcREZ9tO7CNG+bcQHpeOl2adOHDKz6kSd0mTseSynkSeJ6aUa4H8Px7y+XLruBtJ5i+xqc4IiJ4rlJz45wbyczP5NwzzuUfA/9Bw6ia8J5cbX2PBvYBfLtI+fEXOM/HsyJv9X8sEanO/pPxH6794loy8zPp3aI3Hw36SKUq1YYvn7H6UsIiIuWas20Oo74ZRWFpIYPaDOLNfm9qQH2pVk5ZrMaY0lPMcjrXYxWRGmT6uuk8/v3juNwuhp81nBd6v0BEeITTsUT8qiKFeDrXWxURwe1289qK13htxWsAjDl3DPd3v5+wML29SPVTkWLtftztBsAoPJdyCwN+9ncoEak+Sl2lPLbkMT7Y8AHhYeGM6zuOmzre5HQskYA5ZbFaa1cBGGPqAPcAfwbOwHM08LPW2s8CmlBEqqy8ojz++PUfWbhjIXUi6vDX/n9lYNJAp2OJBFRFPmONAkYCDwHNgPXAn6y1MwOcTUSqsPS8dG758hbW719PbFQs713+Hj2b93Q6lkjAVWRX8FagOZ4xgS3wCeAyxlx1ZAZr7azAxBORqmhl5kpum38be/P3cmbDM5n+2+m0adjG6VgiQVGRYm3BfwfcN96vstwVfBwRqQFmb53Nfd/eR0FpAb1b9OadAe8QGxXrdCyRoKlIIabhKU8RkRNyu928+fObvPTTSwAM6zCMsX3GUjuitsPJRIKrIgcvJQUhh4hUYfkl+fx50Z/515Z/EUYYj/d8nLvPvlun00iNpF24InJa0g6mcftXt7N231rqRdbjjX5vcEWbK5yOJeIYFauIVNqi9EWMXDiSnMIckmKSmDJgCh0bd3Q6loijNP6viPjM7Xbz9uq3uWnuTeQU5tC/VX/+75r/U6mKoC1WEfHR4eLDPLjoQWZt9Zxld3/3+3nw3AcJD9Pf6SKgYhURH2zM3shdX93FxpyN1K9Vn4m/majPU0WOo2IVkQr5ZOMnPL7kcfJL8mkb25Z3L3uXdo3aOR1LJOSERLEaYyKAn4B0a+2Vxpg2wMdAE2A5cLO1tsjJjCI11eHiwzy25DFmbvKMYjqk3RDG9hlL/Vr1HU4mEppC5UOR+/CMQXzES8Dr1tq2QDYwwpFUIjVcyv4UBv/vYGZumkmdiDq89pvXmNhvokpV5CQcL1ZjTEtgMPCu93YY0B/4p3eWacA1zqQTqZncbjcfbfiIQf87iI05G2kX244518zhhvY3OB1NJOSFwq7gCXiunBPtvd0EyLHWlnhv7wQSnAgmUhPtL9jPw4sfZk7qHACub3c9Y/uMpV6teg4nE6kaHC1WY8yVQKa1drkxpl8llr8TuBPAWktcXFylckRGRlZ62UAK1VwQutmUyzfH51qwdQF3/N8dZORlEBMVw4QBE7ipa/AvSl5V1leoCNVcNZXTW6x9gKuMMYOAOkAMMBGINcZEerdaWwLp5S1srZ0MTPbedGdlZVUqRFxcHJVdNpBCNReEbjbl8s2RXPkl+fzlx78wZe0UAC5ofgET+02kVXQrR3KH+voKNaeTKz4+3s9pxNHPWK21j1prW3oH+h8KfG2tvQn4Bhjine0W4HOHIopUe2v3rWXw/w5mytopRIZF8sj5jzBz8ExaRbdyOppIleT0FuuJPAx8bIx5AVgJTHE4j0i1U1RaxAuLX+AvS/5CibuEMxueyaRLJnFO03OcjiZSpYVMsVprvwW+9f68FejpZB6R6uyXrF944N8PsG7/OgBu6XQLT/R8QgcoifhByBSriAReUWkRE1dOZNLPkyhxl5AUm8T4PuPpHd/b6Wgi1YaKVaSGWL13NaMXjWb9fs9YLH/o/AdeGfgK+QfzHU4mUr2oWEWqubyiPMYvH8/UtVNxuV0kxSTx6sWv0qtFL+rXrk8+KlYRf1KxilRTbrebeanzePI/T5JxKIPwsHDu6HIHD5//MHUj6zodT6TaUrGKVEM7c3fyxPdPsCBtAQDnxJ3DSxe9RNe4rg4nE6n+VKwi1UhRaRFT107lleWvkF+ST4NaDXjk/EcYftZwIsIjnI4nUiOoWEWqiW92fMPT/3maLQe2AHBlmyt59sJnaV6/ucPJRGoWFatIFbftwDaeXfrs0d2+bWLa8Fzv5+jfqr/DyURqJhWrSBV1qPgQb6x8g8lrJlPkKqJ+rfo80P0BRnQZQe2I2k7HE6mxVKwiVUyJq4SZG2cyfvl49hzeA3gu7fZoz0dpVq+Zw+lERMUqUkW43W6+SvuKsT+OZWPORgC6Ne3G872fp8cZPRxOJyJHqFhFqoCVmSt54YcXWLp7KQCtGrTikfMf4arkqwgPc/QiVSJyHBWrSAhLPZjKuGXj+GLrFwDERsVyf/f7Gd5pOFERUQ6nE5HyqFhFQtDuQ7t58+c3+WD9B5S4S6gTUYcRXUZwzzn30DCqodPxROQkVKwiISTzcCaTVk3ig/UfUFhaSBhhmPaGMeeOIaFBgtPxRKQCVKwiISArP4u3Vr3FtHXTKCgtAGBQm0E82ONBOjbu6HA6EfGFilXEQZmHM3lnzTu8v+59DpccBmBg64GMPnc0nZt0djidiFSGilXEAWkH0/jb6r/xycZPKCwtBOCyxMt4sMeDnN30bIfTicjpULGKBFHK/hQmrZrE51s+p9RdCni2UEd1G0X3M7o7nE5E/EHFKhIEP6b/yAv/foEvt38JQERYBNe1vY57zrmHDo07OJxORPxJxSoSIKWuUuZvn8+7v7x7dGCHqIgohnYYysizR9IqupXDCUUkEFSsIn6WW5TLxykfM3XtVNJy0wCIiYpheMfh3N7ldprWa+pwQhEJJBWriJ9sP7idKWun8EnKJ+QV5wHQOro1I7qMYOSFIynKLXI4oYgEg4pV5DS43C4Wpy9m2rppzN8+HzduAC5scSF3dLmDyxIvIyI8gpioGLJysxxOKyLBoGIVqYR9+fv4ZOMnfLD+A7bnbgev06UCAAAPYUlEQVSgdnhtrml7DSO6jKBLky4OJxQRp6hYRSrI7XazdPdS/rH+H8zZNodiVzEALRu05KaONzGswzB9fioiKlaRU8nKz+KzzZ8xY8MMNudsBiA8LJwBiQO4+ayb6deyHxHhEQ6nFJFQ4WixGmNaAdOBZoAbmGytnWiMaQx8AiQBqYCx1mY7lVNqnqLSIhamLcRusnyd9jUl7hIAmtVrxrAOw7ix440aFF9EyuX0FmsJ8KC1doUxJhpYboxZANwKLLTWjjPGPAI8AjzsYE6pAdxuN2v3rcVutHy2+TOyCz1/y0WERTAgcQCmvWFA6wHUCq/lcFIRCWWOFqu1NgPI8P6ca4xZDyQAVwP9vLNNA75FxSoBkp6Xzqwts/h086es37/+6PSOjTpi2huubXutPjsVkQpzeov1KGNMEtAd+AFo5i1dgN14dhWXt8ydwJ0A1lri4uIq9dyRkZGVXjaQQjUXhG62iubKPJTJpxs+Zea6mSzZueTo9MZ1GjO081CGnz2cbs26ERYWFtRcwaZcvlEuqYgwt9vtdAaMMQ2AfwMvWms/M8bkWGtjy9yfba1tdIqHce/atatSzx8XF0dWVuidYxiquSB0s50s18Gig8xNncvnmz/nu13fHR0Ev05EHS5LvIxrkq+hf2J/oiKigprLScrlm+qYKz4+HsA/f0EKEAJbrMaYWsCnwAxr7WfeyXuMMS2stRnGmBZApnMJpSrbl7+P+dvnMzd1LovTF1Pk8ox+FBkWyaWtLuWattdweeLlNKjdwOGkIlJdOH1UcBgwBVhvrX2tzF2zgFuAcd7vnzsQT6qo9Lx05qXOY27qXH7Y/QMutwuAMMK4sMWFXJN8DYPaDKJxncYOJxWR6sjpLdY+wM3AGmPMz95pj+EpVGuMGQFsB4xD+aQKcLvdbDmwhbmpc1mwYwHLdy8/el+t8Fr0a9mPK5Ku4PLWlxNXV59DiUhgOX1U8HeceN/+pcHMIlVLQUkBSzOW8vWOr1m4YyGpB1OP3lc3si79W/VnUNIg+if2J6Z2jHNBRaTGcXqLVaTC0vPSWZi2kIU7FrJk1xLyS/KP3hcbFcuAxAHccM4NdIvuRt3Iug4mFZGaTMUqISu/JJ+f9vzEop2LWLhjISnZKcfc3zWuK/1b9ad/q/50b9qdiPCIkD1qU0RqDhWrhIxSVylr9q1hcfpivkv/jmV7llFYWnj0/ga1GnBxwsVcmngp/Vr2o3n95g6mFREpn4pVHON2u9l6YCuLd3mK9Ptd33Og6MDR+8MIo0uTLvRN6Ev/Vv05v9n51I6o7WBiEZFTU7FK0LjcLlKyU/gh4wd+2O352nN4zzHzJMUk0Se+DxclXESf+D46JUZEqhwVqwRMsauYNVlr+HH3jyzNWMqyPcvIKcw5Zp4mdZrQN6EvfeM9X4kxiQ6lFRHxDxWr+M3ew3tZuXclyzOXs2LPClbuXXnMkbsA8fXj6dWiFxc0v4BeLXqR3DDZb+PxioiEAhWrVEphSSHL9yxn5d6VrMhcwYo9K9iRt+NX8yU3TOaC5hdwQYsL6NW8Fy2jWzqQVkQkeFSsckolrhK25Gxhzb41rN67mhV7V7B231qKSouOma9eZD26Ne1Gj2Y96NG0Bz3O6KHLrYlIjaNilWMUlRaxMXsja7LWsGbfGtZkrWHdvnUUlBYcM18YYbSPbU+PM3rQo1kPujftTodGHYgIj3AouYhIaFCxApv3b6ahu2GN+6wvpzCHlP0prN+/nl/2/cKarDWkZKdQ7Cr+1bytGrSiS1wXujTpQo9mPbi0w6UU5/16PhGRmq7GF+uT3z/Je2vf44OBH9CvVT+n4wREfkk+m7I3sSF7Axv2byAlO4UN+zew+/Ducuc/s+GZdI3rStcmXY+WaaM6x14Ot2GdhmTlaYQjEZHj1fhija8fjxs3L//0Mr9p+ZsqvdVaUFJA6sFUNuVsIiU75ejWaOrBVNz8+oL2dSLq0KFRBzo07kCXJp4C7dSkE9G1ox1ILyJSPdT4Yr218628s/YdVmWtYv72+fw26bdORzopt9vNnsN72JyzmS0HtrDlwBa25mxly4Et7MjdUW6BRoRFkNwwmQ6NO9CxUUc6Nu5Ih0YdSIxO1GeiIiJ+VuOLtW5kXR7u/TCjF4xm7LKx9E3oS/1a9R3N5Ha7yS7MZmv6VlbvWM22g9vYkrPlaJEeKj5U7nIRYREkRieSHJt8TIEmxyYTFREV5H+FiEjNVOOLFeD2brfz1x//yuaczYz6ZhTvXvZuwLfkikqLSM9LJy03jdSDqaTlppF28L8/5xXnnXDZ2KhY2sa2Jblhsucr1vO9dUxrjaUrIuIwFSsQFRnFe5e/x1WzrmL+9vnc8uUtvPqbV2lWr1mlH7OwtJCMQxmk56WzK28X6XnppOelsz13O2kH00g/lI7L7Trh8tG1oklunEx8vXiSopOOlmdybLLGzxURCWEqVq/k2GSmDJjCiAUj+GbnN/T+uDeD2wzm4pYX0z62PU3qNgE8lzYrdhWTU5jDvoJ97C/Yz778fWTmZ7Irbxe7DnlKdG/+3pM+XxhhJDRIIDE6kaSYJBKjE2kd05rWMa1JjE6kUVQjmjZtqmuLiohUMSrWMnq16MXC6xby2JLH+HL7l3y6+VM+3fxppR4rIiyC5vWbk1A/gYQGCcQ3iCe+QTytoz3lmdAgQZ97iohUQyrW4zSv35ypl08l9WAqs7fO5ue9P7P94HayC7Nx4yYyLJLI8Ega1m5Ik7pNaFynMU3qNKFp3aZHyzOhfgLN6jXTEbciIjWQivUEkmKSGNVtlNMxRESkigl3OoCIiEh1omIVERHxIxWriIiIH6lYRURE/EjFKiIi4kcqVhERET9SsYqIiPiRilVERMSPwtzuX1+/s4qqNv8QEZEgC3M6QHVSnbZYwyr7ZYxZfjrLB+orVHOFcjblUi7lqtSX+FF1KlYRERHHqVhFRET8SMXqMdnpACcQqrkgdLMpl2+UyzfKJadUnQ5eEhERcZy2WEVERPyoxl+P1RgzEJgIRADvWmvHOZSjFTAdaIbn1KHJ1tqJxphngDuAvd5ZH7PWzglytlQgFygFSqy15xljGgOfAElAKmCstdlBzNTB+/xHnAk8BcQS5PVljJkKXAlkWmu7eKeVu36MMWF4Xm+DgMPArdbaFUHMNR74HVAEbAFus9bmGGOSgPVAinfxpdbau4OY6xlO8HszxjwKjMDz+rvXWvtlEHN9AnTwzhIL5FhruwV5fZ3ovcHx15iUr0YXqzEmAvgrMADYCSwzxsyy1q5zIE4J8KC1doUxJhpYboxZ4L3vdWvtKw5kKusSa21WmduPAAutteOMMY94bz8crDDW2hSgGxz9PaYD/wJuI/jr631gEp43vyNOtH6uANp5vy4A/ub9HqxcC4BHrbUlxpiXgEf57+9ti7W2W4CynCoXlPN7M8Z0AoYCnYF44CtjTHtrbWkwcllrbyiT5VXgQJn5g7W+TvTecCvOv8akHDV9V3BPYLO1dqu1tgj4GLjaiSDW2owjf1Vaa3Px/DWc4ESWCroamOb9eRpwjYNZLsXzJrfdiSe31i4C9h83+UTr52pgurXWba1dCsQaY1oEK5e1dr61tsR7cynQMhDP7Wuuk7ga+NhaW2it3QZsxvP/Nqi5vFuBBvgoEM99Mid5b3D8NSblq9FbrHhenDvK3N5JCPxl593N1B34AegDjDLGDAd+wvOXa9B2uXq5gfnGGDfwd2vtZKCZtTbDe/9uPLupnDKUY9/wnF5fcOL1U95rLgHIIPj+wLG709sYY1YCB4EnrLWLg5ynvN9bAp4/AI44sr6C7SJgj7V2U5lpQV9fx703VIXXWI1U07dYQ44xpgHwKXC/tfYgnt04yXh2e2YArzoQq6+1tgeeXUz3GGMuLnuntdaNQ0NKGmNqA1cBM72TQmF9HcPJ9XMixpjH8exinOGdlAEkWmu7A6OBD40xMUGMFHK/t+MM49g/3oK+vsp5bzgqFF9jNVlNL9Z0oFWZ2y290xxhjKmF5z/ODGvtZwDW2j3W2lJrrQt4hwDtBjsZa22693smns8xewJ7juxe8n7PDHYuryuAFdbaPd6Mjq8vrxOtH8dfc8aYW/EcpHOT9w0Z767Wfd6fl+M5sKl9sDKd5PcWCusrEriWMlv3wV5f5b03EMKvsZquphfrMqCdMaaNd8tnKDDLiSDez3CmAOutta+VmV72s5H/AX4Jcq763gMmMMbUBy73ZpgF3OKd7Rbg82DmKuOYLQmn11cZJ1o/s4DhxpgwY0wv4ECZ3XkB5z0K/iHgKmvt4TLTm3oPAsMYcyaeA1+2BjHXiX5vs4ChxpgoY0wbb64fg5XL6zJgg7V255EJwVxfJ3pvIERfY1LDP2P1Hhk5CvgSz+k2U621ax2K0we4GVhjjPnZO+0xYJgxphue3TypwF1BztUM+JcxBjyvlw+ttfOMMcsAa4wZAWzHc2BHUHmLfgDHrpOXg72+jDEfAf2AOGPMTuBpYBzlr585eE6D2IznVIjbgpzrUSAKWOD9nR45TeRi4DljTDHgAu621lb0ACN/5OpX3u/NWrvWGGOBdXh2Xd8ToCOCy81lrZ3Crz/DhyCuL0783uD4a0zKp5GXRERE/Kim7woWERHxKxWriIiIH6lYRURE/EjFKiIi4kcqVhERET+q0afbiPiDMeZ9POcRPmutfcbZNCLiNBWryCl4L5vX+iSz3ATkcOyYtiJSQ6lYRU5tKtDY+/NIoDae4eWOjMTzo7X2QyeCiUjo0QARIj4wxuQADfFcn/Zb77T3KbMr2HvR7qeBecA+POPM/gLciGc4wRvxjC40/MjlwIwxicDLeK6iUh/PsH2jrbVODckoIpWkg5dEAue3QBMgDTgfz+XQugOr8Vy4+w0AY0w94Gs8Q9KtxjPWaz/ga2NMXNBTi8hpUbGKBM4WPGO2jvPerodnQPcj4xd3934fjOeSabuAFDwX204DmgJDghVWRPxDxSoSOBu8l2XL8d7eY609AOR6b9fzfk/yfk8A7vN+JXuntQ1CThHxIx28JBI4x1+F5URXZUn1fl8OnH/kGqnGmFggLDDRRCRQVKwizpsDbAPOBZYYY1YDiXg+Zx0EfOtYMhHxmXYFizjMWnsI6I/nmp+JeI4w7gB8gOczVxGpQnS6jYiIiB9pi1VERMSPVKwiIiJ+pGIVERHxIxWriIiIH6lYRURE/EjFKiIi4kcqVhERET9SsYqIiPiRilVERMSP/h/3hX+wFRh+bAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#sLline = plt.plot(\"t\",\"SL\",\"\",data=output_SL,color=\"red\",linewidth=2)\n",
    "#sHline = plt.plot(\"t\", \"SH\", \"\", data = output_SL, color = \"red\", linewidth = 2)\n",
    "iline = plt.plot(\"t\",\"I\",\"\",data=output_LIB,color=\"green\",linewidth=2)\n",
    "#rline = plt.plot(\"t\",\"R\",\"\",data=output_SL,color=\"blue\",linewidth=2)\n",
    "plt.xlabel(\"Time\",fontweight=\"bold\")\n",
    "plt.ylabel(\"Number\",fontweight=\"bold\")\n",
    "legend = plt.legend(title=\"Population\",loc=5,bbox_to_anchor=(1.25,0.5))\n",
    "frame = legend.get_frame()\n",
    "frame.set_facecolor(\"white\")\n",
    "frame.set_linewidth(0)"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}