notebooks/kr08/2_1/c_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program 2.1 (SIR model) - original C code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting parameters.txt\n"
     ]
    }
   ],
   "source": [
    "%%writefile parameters.txt\n",
    "Transmission_rate,_beta\n",
    "\t\t\t1.4247\n",
    "Recovery_rate,_gamma\n",
    "\t\t\t0.14286\n",
    "Initial_Susceptibles\n",
    "\t\t\t0.999999\n",
    "Initial_Infecteds\n",
    "\t\t\t0.000001\n",
    "Run_until_time\n",
    "\t\t\t70"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting Program_2_1.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile Program_2_1.c\n",
    "/* ****************************************************************\n",
    "\n",
    "This is the C version of program 2.1 from page 19 of\n",
    "\"Modeling Infectious Disease in humans and animals\"\n",
    "by Keeling & Rohani.\n",
    "\n",
    "It is the simple SIR epidemic without births or deaths.\n",
    "\n",
    "This code is written to be simple, transparent and readily compiled.\n",
    "Far more elegant and efficient code can be written.\n",
    "\n",
    "This code can be compiled using gnu or intel C compilers:\n",
    "icc -o Program_2_1  Program_2_1.c  -lm\n",
    "g++ -o Program_2_1  Program_2_1.c  -lm\n",
    "\n",
    "If you have any difficulties, we suggest you set DEBUG to 1 or -1\n",
    "\n",
    "******************************************************************/\n",
    "\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <math.h>\n",
    "#include <string.h>\n",
    "\n",
    "// Set DEBUG to 1 to output some debugging information\n",
    "// Set DEBUG to -1 to output more debugging information\n",
    "#define DEBUG 0\n",
    "\n",
    "// Set up basic parameters\n",
    "double beta=520.0/365.0;\n",
    "double gamm=1.0/7.0;\n",
    "double S0=1.0-1e-6;\n",
    "double I0=1e-6;\n",
    "\n",
    "double MaxTime=70;\n",
    "\n",
    "// Set up variables and rates of change\n",
    "double t,S,I,R,Pop[3];\n",
    "double dPop[3];\n",
    "\n",
    "// Initialise the equations and Runge-Kutta integration\n",
    "void Diff(double Pop[3])\n",
    "{\n",
    "  // Set up temporary variables to make the equations look neater\n",
    "\n",
    "  double tmpS, tmpI, tmpR;\n",
    "\n",
    "\n",
    "  tmpS=Pop[0]; tmpI=Pop[1]; tmpR=Pop[2];\n",
    "\n",
    "  /* The differential equations */\n",
    "\n",
    "  dPop[0] = - beta*tmpS*tmpI;              // dS/dt\n",
    "  dPop[1] = beta*tmpS*tmpI - gamm*tmpI;   // dI/dt\n",
    "  dPop[2] = gamm*tmpI;                    // dR/dt\n",
    "\n",
    "  return;\n",
    "}\n",
    "\n",
    "\n",
    "void Runge_Kutta(double);\n",
    "void Read_in_parameters(FILE *);\n",
    "void Perform_Checks();\n",
    "void Output_Data(FILE *);\n",
    "\n",
    "\n",
    "// The main program\n",
    "\n",
    "int main(int argc, char** argv)\n",
    "{\n",
    "  char FileName[1000];\n",
    "  double step, Every;\n",
    "  FILE *in, *out;\n",
    "   \n",
    "  /* Find the given parameter  file name and open the file \n",
    "     otherwise use parameters set at top of program */\n",
    "  if(argc>1)\n",
    "    {\n",
    "      strcpy(FileName, *(argv+1));  \n",
    "      \n",
    "      if((in = fopen(FileName,\"r\"))==NULL)\n",
    "\t{\n",
    "\t  printf(\"Cannot read file %s\\n\",FileName);\n",
    "\t  exit(1);\n",
    "\t}\n",
    "\n",
    "      /* Read in the parameters */\n",
    "      \n",
    "      if(DEBUG) printf(\"\\nReading parameters from file %s\\n\",FileName);\n",
    "\n",
    "      Read_in_parameters(in);\n",
    "\t  \n",
    "    }\n",
    "  else\n",
    "    {\n",
    "      if(DEBUG) printf(\"\\nUsing default parameters\\n\");\n",
    "    }\n",
    "\n",
    "  if(DEBUG) printf(\"beta=%g\\ngamma=%g\\nInitial S=%g\\nInitial I=%g\\nRuns until time %g\\n\\n\",beta,gamm,S0,I0,MaxTime);\n",
    "\t\n",
    "\n",
    "  /* Check all parameters are OK & set up intitial conditions */\n",
    "  \n",
    "  Perform_Checks();\n",
    "\n",
    "  S=S0; I=I0; R=1-S-I;\n",
    "\n",
    "  /* Find a suitable time-scale for outputs */\n",
    "\n",
    "  step=0.01/((beta+gamm)*S0);\n",
    "  \n",
    "  Every=pow(10,floor(log10((1.0/((beta+gamm)*S0)))));\n",
    "  while(MaxTime/Every>10000)\n",
    "    {\n",
    "      Every*=10.0;\n",
    "    }\n",
    "  if(DEBUG) printf(\"Using a time step of %g and outputing data every %g\\n\\n\",step,Every);\n",
    "   \n",
    "\n",
    "  if((out = fopen(\"Program_2_1_c.out\",\"w\"))==NULL)\n",
    "    {\n",
    "      printf(\"Cannot open file Output to write data\\n\");\n",
    "      exit(1);\n",
    "    }\n",
    "\n",
    "  /* The main iteration routine */\n",
    "\n",
    "  t=0;\n",
    "\n",
    "  Output_Data(out);\n",
    "\n",
    "  \n",
    "\n",
    "  do\n",
    "    {\n",
    "      Runge_Kutta(step);  t+=step;\n",
    "      \n",
    "      /* If time has moved on sufficiently, output the current data */\n",
    "      if( floor(t/Every) > floor((t-step)/Every) )\n",
    "\t{\n",
    "\t  Output_Data(out);\n",
    "\t}\n",
    "      else\n",
    "\t{\n",
    "\t  if(DEBUG<0)  Output_Data(stdout);\n",
    "\t}\n",
    "      \n",
    "    }\n",
    "  while(t<MaxTime);\n",
    "\n",
    "\n",
    "  Output_Data(out);\n",
    "\n",
    "  fclose(out);\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "void Read_in_parameters(FILE *in)\n",
    "{\n",
    "  char str[200];\n",
    "  fscanf(in,\"%s\",str);\n",
    "  fscanf(in,\"%lf\",&beta);\n",
    "  \n",
    "  fscanf(in,\"%s\",str);\n",
    "  fscanf(in,\"%lf\",&gamm);\n",
    "  \n",
    "  fscanf(in,\"%s\",str);\n",
    "  fscanf(in,\"%lf\",&S0);\n",
    "  \n",
    "  fscanf(in,\"%s\",str);\n",
    "  fscanf(in,\"%lf\",&I0);\n",
    "  \n",
    "  fscanf(in,\"%s\",str);\n",
    "  fscanf(in,\"%lf\",&MaxTime);\n",
    "  \n",
    "  fclose(in);\n",
    "}\n",
    "\n",
    "\n",
    "\n",
    "void Perform_Checks()\n",
    "{\n",
    "  if(S0<=0) {printf(\"ERROR: Initial level of susceptibles (%g) is less than or equal to zero\\n\",S0); exit(1);}\n",
    "  \n",
    "  if(I0<=0) {printf(\"ERROR: Initial level of infecteds (%g) is less than or equal to zero\\n\",I0); exit(1);}\n",
    "  \n",
    "  if(beta<=0) {printf(\"ERROR: Transmission rate beta (%g) is less than or equal to zero\\n\",beta); exit(1);}\n",
    "  \n",
    "  if(gamm<=0) {printf(\"ERROR: Recovery rate gamma (%g) is less than or equal to zero\\n\",gamm); exit(1);}\n",
    "\n",
    "  if(MaxTime<=0) {printf(\"ERROR: Maximum run time (%g) is less than or equal to zero\\n\",MaxTime); exit(1);}\n",
    "  \n",
    "  if(S0+I0>1) {printf(\"WARNING: Initial level of susceptibles+infecteds (%g+%g=%g) is greater than one\\n\",S0,I0,S0+I0);}\n",
    "\n",
    "  if(beta<gamm) {printf(\"WARNING: Basic reproductive ratio (R_0=%g) is less than one\\n\",beta/gamm);}\n",
    "}\n",
    "\n",
    "\n",
    "void Output_Data(FILE *out)\n",
    "{\n",
    "  if(out!=stdout) fprintf(out,\"%g\\t%g\\t%g\\t%g\\n\",t,S,I,R);\n",
    "\n",
    "  if(DEBUG) printf(\"%g\\t%g\\t%g\\t%g\\n\",t,S,I,R);\n",
    "}\n",
    "\n",
    "\n",
    "void Runge_Kutta(double step)\n",
    "{\n",
    "  int i;\n",
    "  double dPop1[3], dPop2[3], dPop3[3], dPop4[3];\n",
    "  double tmpPop[3], initialPop[3];\n",
    "\n",
    "  /* Integrates the equations one step, using Runge-Kutta 4\n",
    "     Note: we work with arrays rather than variables to make the\n",
    "     coding easier */\n",
    "\n",
    "  initialPop[0]=S; initialPop[1]=I; initialPop[2]=R;\n",
    "\n",
    "  Diff(initialPop);\n",
    "  for(i=0;i<3;i++)\n",
    "    {\n",
    "      dPop1[i]=dPop[i];\n",
    "      tmpPop[i]=initialPop[i]+step*dPop1[i]/2;\n",
    "    }\n",
    "\n",
    "  Diff(tmpPop);\n",
    "  for(i=0;i<3;i++)\n",
    "    {\n",
    "      dPop2[i]=dPop[i];\n",
    "      tmpPop[i]=initialPop[i]+step*dPop2[i]/2;  \n",
    "    }\n",
    "\n",
    "  Diff(tmpPop);\n",
    "  for(i=0;i<3;i++)\n",
    "    {\n",
    "      dPop3[i]=dPop[i];\n",
    "      tmpPop[i]=initialPop[i]+step*dPop3[i]; \n",
    "    }\n",
    "\n",
    "  Diff(tmpPop);\n",
    "\n",
    "  for(i=0;i<3;i++)\n",
    "    {\n",
    "      dPop4[i]=dPop[i];\n",
    "\n",
    "      tmpPop[i]=initialPop[i]+(dPop1[i]/6 + dPop2[i]/3 + dPop3[i]/3 + dPop4[i]/6)*step;\n",
    "    }\n",
    "\n",
    "\n",
    "  S=tmpPop[0]; I=tmpPop[1]; R=tmpPop[2];\n",
    "\n",
    "  return;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gcc -O3 -o Program_2_1_c  Program_2_1.c  -lm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "!./Program_2_1_c parameters.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "sir_out = pd.read_csv(\"Program_2_1_c.out\",sep=\"\\t\",header=None,names=[\"t\",\"S\",\"I\",\"R\"],index_col=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdMAAAENCAYAAABUyo/yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8FHX++PHXbEk2nSQbCKGFHroVpNgbAgc2RuEQva93lhP1zjvL/c6v+uX8eqd4p1hOz/O+omf9WLAgduBQFMRG7zWEBEgjvezu/P6Y3ZCEtE02W8j7+XjMY2dnPjPz3rDknU+Zz2iGYSCEEEKI9rOEOgAhhBAi0kkyFUIIITpIkqkQQgjRQZJMhRBCiA6SZCqEEEJ0kCRTIYQQooMkmQohhBAdJMlUCCGE6CBJpkIIIUQH2UIdQAfJ9E1CCNE+WqgDOJFEejLl4MGD7TrO6XSSn58f4Gg6TyTFG0mxQmTFG0mxQmTFG0mxQsfizcjICHA0Qpp5hRBCiA6SZCqEEEJ0kCRTIYQQooMkmQohhBAdJMlUCCGE6KCgjObVdf3/gGnAYaXUyCb2a8BCYApQAVynlPohGLEJIYQQHRWsmukiYHIL+y8BBnuXG4BnghCTEEIIERBBqZkqpVbqup7ZQpEZwEtKKQNYret6N13XeyqlcjsjnsQHH8S2YQOpLhdo2nGLYbE0ub3VfXY7RkICnoQEjIQE3E4n7v79cQ0YgMfpNMsJ0YUYBrhcUFurUV1tvtbUmK++bW63httd/5UG7z0ecLmOrbvdGi7XsfWWyhuGb9HweI6993jM+Mxt5j6Hw0pFRWKT+46d59jxhqHVWz9+n+/zG8axn0X91/o/o+bK1D9P47Lnn2/huusC+I8lOiRcJm3oBWTXe3/Au+24ZKrr+g2YtVeUUjidTr8vZtu9G8tXXxHdvljbxejXD88FF+CZOxdj3Di/E6vNZmvXZw2FSIoVIiveQMXqdkNZGRw9CqWlGkePQkkJlJVpVFZCRYW5mOta3Xrj9xUVUF0NNTW+JAk1Nb5t5mIYkTRBQHyoA2gzp9OImO9tVxAuybTNlFLPAc953xrtmQHEdtddJN95J0eLihr8OalB4z8vG/wZqzW1r/7+6mosZWVopaVYSkuxHDqEbc8ebDt3Ytm3D+u//oX1X/+ieuJEjj70EK5Bg9occyTNzhJJsUJkxds41poaKCy0HLcUFTVcLy21UFqq1b2Wlwdv7KHVamC3G0RFUfcaFWVus9vBZjOw2cBiMctarXgXc91iMcuY+xuuNy5vlj227mtIMl+NBo1JjffFxcVRWVneYJ+5blYFG5ZvfZ+mHauC+v52bum1+X0Nq7K+skOHJsoMSGEkXJJpDtCn3vve3m2dwpWVheF0UhOsX6BuN/aNG3F8+CFxr7xC9KpVOC++mOInnqBq6tTgxCAiimFAQYGFAwes5OVZycuzkJdn5ehRK/v2pZCXZ+XQISvFxe1PivHxHhISDBISjr3GxRnExDRcHI7m3zscBtHR9RPlsXW73SAjw0lRUaT8oRJDfn55qMNoM6fTIEL+BuwSwiWZvg/M03X9dWAccLSz+ktDwmqldswYaseMoeymm0h64AFi336b5Jtvpuif/6Tq4otDHaEIAbcb9u2zsmuXjf37bezbZyU728r+/Tb277dSUdFcorTWrdlsBikpHlJSPCQne+rW6y/dunlISPCQmGjUJdD4eLPm1tms1tbLCHEiCNatMa8B5wBOXdcPAPcDdgCl1LPAUszbYnZi3hrzi2DEFQpGSgrFCxfi7tWLhCeeoNutt3Lkk09w9+8f6tBEJ/F4YPduK9u22dm+3caOHTa2b7eze7eN6urm+867dfPQu7eL9HQP6elu0tPdDBoUS1zcUXr0cNOzp5lAZVybEKEXrNG8s1rZbwC3BCOWsKBplN51F7Y9e4j54AO6/f73FLz1loz2PQH4EueGDVGsW2dnwwY7GzfaKStruhrYs6ebQYNcZGa66NfPRZ8+bvr1c9Onj4tu3Y5/wqDT6SA/v7qzP4YQwk/h0szb9WgaxX/5C1GrVhG9ejWOpUul/zQC1dbChg12vv02itWro1m7NqrJfsz0dDfDh9cyZIiLIUNqGTzYxeDBLhIS5JG8QpwIJJmGkNGtG6W/+x3d/vhH4p98kqopU6R2GuYMA3butLFsWTQrVpjJs7KyYfLs0cPNmDE1jB5dW7ekpXlCFLEQIhgkmYZYxVVXkfC3vxG1YQNRX39NzcSJoQ5JNFJdDStXRvPZZw5WrIgmJ6fhf5sBA1yccUY1Y8fWMG5cDX36uOVvIiG6GEmmoRYTQ8XcuSQ89hixb7whyTRMVFXBf/7jYMkSB5995qC09FjtMzXVzdlnV3PuudWceWa11DqFEJJMw0HFFVeQ8NhjOD76CK2iAiM2NtQhdUmGAWvXRvH667EsWeJoMLHBiBG1XHJJJeefX83IkbVBua1ECBE5JJmGAXf//tSccgpRP/xA9LJlVE2bFuqQupS8PHj22Xhefz2WPXuO/ZcYNaqGqVOrmDq1kgED3CGMUAgR7iSZhomqiy82k+ny5ZJMg2T9ejv//GccH3xgp7Y2CjAHD82cWcHMmRUMGiQJVAjRNpJMw0TVeeeR+Oc/41i+nKOGIaN6O4lhwGefRfPss/GsWWM+6sBiMZg8uZJZsyo455xqbPK/QgjhJ/m1ESZcw4bhTk/HmpeHbds2XFlZoQ7phGIY8OmnDv72t3g2bjRroQkJHmbPruCOO6KIjy8KcYRCiEgmyTRcaBrV48YR+957RK1dK8k0QAwDli2L5pFHEuqSaI8ebm6+uYxZsyqIjze8T2IJcaBCiIgmYxLDSM3ppwMQtXZtiCM5MWzdamP27BTmzk1l48YoevRwM3/+UVatOsSvflVOfLzMPiSECAypmYaRumT63XchjiSyFRZaeOSRBF55JRaPRyMpycPtt5cyd245MTGhjk4IcSKSZBpGXMOG4YmNxbZvH5bCQjwpKaEOKaIYBrzzTgz3359IUZEVq9XgF78o4447ykhJkYkVhBCdR5p5w4nVimvYMABsmzaFOJjIkp1tZc6cFG67LZmiIiuTJlXzxRdHePDBEkmkQohOJ8k0zNSOHAmAXZJpmxgGvPFGDOedl8aKFQ66dfPw2GNFvP56AYMHu0IdnhCii5Bm3jBTO2IEIMm0LYqLNe65pxsffGB2hE6bVsmDDx6VuXKFEEEnyTTMSDJtm7Vr7dxySzI5OTbi4jw89NBRrriiUua6EEKEhCTTMOMaMgQA2+7d4HIh0/Ec7+WXY7n33iRqazVOPrmGp54qIjNTpv4TQoSO9JmGGSM2FlevXmi1tVj37w91OGGlpgbuvjuJu+/uRm2txvXXl7F4cb4kUiFEyEm1Jwy5Bg3ClpODbdcu3AMGhDqcsFBYqHH99Sl8+2000dEGDz9czMyZlaEOSwghAKmZhiXXwIEA2HbtCnEk4SE728qllzr59ttoevZ0s3hxviRSIURYkZppGKpLpjt3hjiS0Nu40cY116Ry+LCVYcNqefnlAtLTZbSuECK8SM00DNUl0927QxxJaK1da+fKK50cPmxlwoRq3nknXxKpECIsSTINQ+5+/QCwdeEBSN9+G8XPf55KaamFqVMrefnlAhITZWJ6IUR4kmQahtwZGRhWK5a8PKiuDnU4Qbd6dRQ//3kK5eUWLrusgr//vYjo6FBHJYQQzZNkGo5sNtwZGWiGgfXAgVBHE1TffWdnzpwUKiosXH55BQsXFsuttkKIsCfJNEy5+/QBwJadHeJIgmfbNhvXXptKZaWFK66o4PHHi7FaQx2VEEK0TpJpmHL17QvQZSZuyMmxMnt2KsXFFi66qJK//U0SqRAickgyDVN1NdMukEyLijRmzUohL8/KuHHV/P3vRdK0K4SIKJJMw5QvmZ7ofaa1tXDjjSns2mVn2LBaXnihkJiYUEclhBD+Cdrf/7quTwYWAlbgeaXUXxrt7wu8CHTzlrlHKbU0WPGFG3fPngDmiN4T2P33J7FqVTTdu7t56aUCkpLk9hchROQJSs1U13Ur8DRwCTAcmKXr+vBGxe4FlFLqZOBq4O/BiC1cudPTAbCewMl00aJYXnwxjuhog3/9q5CMDJmQQQgRmYLVzDsW2KmU2q2UqgFeB2Y0KmMAid71JOBgkGILSx5vzdR66BAYJ15tbc2aKO67LwmABQuKOeWU2hBHJIQQ7ResZNoLqH+PxwHvtvoeAOboun4AWArcGpzQwpMRE4MnKQmtpgZLYWGowwmoggILv/51Mm63xs03l3HFFTJpvRAisoXTmMlZwCKl1F91XR8P/FvX9ZFKqQZtf7qu3wDcAKCUwul0tutiNput3ccGTa9ecPQoqZWVWCMhXq+WfrYeD1x3nY28PAsTJnh49NEobLbQfq6I+C54RVKsEFnxRlKsEHnxnuiClUxzgD713vf2bqvvemAygFLqG13XHYATOFy/kFLqOeA571sjPz+/XQE5nU7ae2ywpHTvjmPzZkq2bCHhpJPCPl6fln62Tz4Zz2efJZKc7GbhwiMUF4e+nzQSvgs+kRQrRFa8kRQrdCzejIyMAEcjgpVM1wKDdV3vj5lErwZmNyqzHzgfWKTr+jDAARwJUnxh6UQbhLR2rZ1HHkkA4IknimXAkRDihBGUPlOllAuYB3wCbDE3qU26rs/XdX26t9jvgF/pur4OeA24Til14o288YPnBEqmFRUav/lNMh6P2U963nldbwJ/IcSJK2h9pt57Rpc22nZfvfXNwMRgxRMJfDXTE+Fe04ceSmDvXhvDhtVy550loQ5HCCECSmZACmN1zby5uSGOpGO++iqKF16Ix2YzePxxeZyaEOLEI8k0jPlmQYrkZt7SUo077ugGwG9+U8rIka4QRySEEIEnyTSMNZi4IUI98kgCOTk2Ro+uYd68slCHI4QQnUKSaRjzJCdjWCxYiouhpibU4fht/Xo7ixbFYbUaPPpoMXZ7qCMSQojOIck0nFkseFJTzfUjkXWXkNsN99yThMejcf315YwYIc27QogTlyTTMOfxznCiHT7cSsnw8u9/x7JuXRQ9e7r53e9KQx2OEEJ0KkmmYc6XTImgZHroEDz8sPnMgvnzjxIf36VvFxZCdAGSTMOcOy0NiKya6f33WykpsXDeeVVccklVqMMRQohOJ8k0zNX1mUZIMt282caiRRZsNoMHHjiKpoU6IiGE6HySTMOcJ8Jqpv/7v4kYhsbcueUMHOgOdThCCBEUkkzDnNvXZxoBo3lXrIhmxQoHSUkGv/2t3FMqhOg6JJmGubrRvGE+cYPbDX/6kzno6O673aSkyBNhhBBdhyTTMOdr5g33mumbb8awdaud3r1d3HKLJFIhRNciyTTMub0DkMK5z7SmBh57zHxO6d13l+JwhDggIYQIMkmmYc5Tv8/UE541PqViOXDAxuDBtcyYURnqcIQQIugkmYa76Gg8iYloLhfa0aOhjuY41dWwcGE8AHfcUYrVGuKAhBAiBCSZRgBf7dSanx/iSI732muxHDxoIyurlmnTZIIGIUTXJMk0Avhuj7GEWTKtroYnnzT7Su+4oxSLfJuEEF2U/PqLAJ7kZAAsRUUhjqSht96KJS/PyrBhtTJtoBCiS5NkGgE8KSkAWAoLWy27pXALn+77lIraik6Nye2GZ54x+0pvuaVMaqVCiC7NFuoAROvaWjN9+qeneWjtQwD0TejLa1NeIzMxs1Ni+vhjB3v22OjTx8XPfiYjeIUQXZvUJyJAW2qmGws28ue1f8aiWegV34v9pfu55uNrqHQFPtEZBjz9tFkrvemmMmzyJ5kQoouTZBoB2lIzfXbdsxgY/GL4L1h2xTKGdBvC7qO7efzHxwMez6pVUaxbF0VKipurrpJaqRBCSJ0iAtTVTJtJpvmV+by/+32smpUbR99IfFQ8C85awIz3Z/D8hue5fsT1dI/tHrB4/vEPs1b6X/9VTkyMPPhbiC5sAvAnICnUgQTBUeC/ga+b2ik10whg+GqmzTTzLs9ejttwMyljEr3iewFwWo/TuLjfxVS5q3huw3MBi2XXLivLljlwOAyuvbY8YOcVQkSkrpJIwfycf2pupyTTCOBupZl3WfYyAM7ve36D7beedCsAr217LWB9p4sWxQFw2WUVpKRIrVSILq6rJFKfZj+vJNMIYLTQzGsYBitzVgJwXp/zGuw7ufvJjHGOobi6mPd3vd/hOEpLNZSKBeAXv5BaqRBC+EgyjQCepCQMTTPn5nW5GuzbU7KH4upiesT2oH9S/+OOvXbEtQC8uPnFDsfx5puxlJVZGD++mhEjXK0fIITocqxW66lZWVnDBw8ePOKSSy4ZUFpaGtA888QTT6TOnTu3b0tllixZkvDZZ5/F+d4/8sgjaU899VRqIONoTJJpJLBaITkZzTCwNJrs/qcjPwEwJm1Mk4dOHzCdxKhE1uWvY0fRjnaH4PHACy+Y302plQohmhMdHe3ZunXr5h07dmyy2+3GX//617Rgx7Bs2bKEL7/8Mt73/q677joyb968gs68piTTSOF9rmnjpt66ZOpsOpnG2GKY2n8qAO/sfKfdl//yy2h277aRkeHi4otl6kAhROsmTZpUtnPnzmiABx54oMfgwYNHDB48eMT8+fO7A2zbti2qf//+I6ZPn95/wIABIyZPnlxXk+3Vq9eo3NxcG8DKlStjx44dO7Tx+V999dWk0aNHZw0bNmz4hAkThmRnZ9u2bdsW9dJLL6U9++yzPbKysoZ//PHH8XfccUfGfffd1wPg66+/jhkzZkzWkCFDhl944YUDjxw5YgUYO3bs0JtvvrnXqFGjhmVmZo78+OOP4xtfryVtSqa6rtt0Xf9C1/VmRzKJzmX4kmmjEb3rj6wH4KS0k5o99vJBlwNmMvUY7Xsm6ssvm32lc+ZUyCQNQohW1dbW8sknnySOGjWq8ssvv4x99dVXU7///vst33333ZaXXnopbdWqVTEAe/fudcybN+/w7t27NyUkJHgWLFjQ5prshRdeWPbTTz9t3bJly+Yrr7yycP78+elDhw6tmTt37pGbbrrp0NatWzdPnjy5rP4x1113Xf+HHnrowPbt2zePGDGi8u67787w7XO5XNqGDRu2PPzww9nz58/POP6KzWtTMlVKuYAsIN2fk9en6/pkXde36bq+U9f1e5opo+u6vlnX9U26rr/a3mudkJqomRqGwbaibQAMTx3e7KFn9DyDjLgMDpQdYG3eWr8vfeSIhU8/dWC1Glx1VefO+SuEiGzV1dWWrKys4aNGjRreu3fvmttvvz1/xYoV8VOmTClOTEz0JCUleaZOnVq0fPnyBID09PSaiy66qBzgmmuuKfj666/bXCPcs2dP1Jlnnjl4yJAhw5944on0rVu3xrRUvqCgwFpaWmqdOnVqGcCvfvWrgtWrV9ddb+bMmUUAEyZMKD9w4ECUP5/bn2be+cCluq6freu63Z+L6LpuBZ4GLgGGA7N0XR/eqMxg4A/ARKXUCOA3/lzjRNdUzfRI5RFKakpIikoiLab5P+YsmqVB7dRfb74Zi8ulcf75VaSnt69mK4ToGnx9plu3bt384osvZjscjhbvodM0rcn3VqvV8HjM3zeVlZVN5qp58+b1/fWvf314+/btm5966ql91dXVHeq69MVqs9lwu91aa+Xr8+fCzwApwDKgStd1t3dpy7DOscBOpdRupVQN8Dowo1GZXwFPK6WKAJRSh/2I7cTXRM10R7E5oGhQt0HHfSEbmzHQ/HF/tPcj3B53my9rGPDqq2YT7+zZUisVQvjv3HPPLVu6dGm30tJSS0lJiWXp0qXJ5557bilAbm5u1Oeffx4H8Morr6RMmDChDKB37941q1atigVQSiU3dd7S0lJr3759awEWLVpUN1o3ISHBXVpaam1cPjU11Z2YmOj29Yf+61//Sh0/fnxZ43Lt4W/vV1O/sduSvXsB2fXeHwDGNSozBEDX9VWAFXhAKfVx4xPpun4DcAOAUgqn98HZ/rLZbO0+NhS0NLPmGVdVhcMbd+6+XABGpY9q9bOcmXomA5MHsqtoF1srt3J2v7PbdN3//Edjzx4bGRkGM2cmYLMltHpMpP1sIyneSIoVIiveSIoVIiveSZMmVcyePbvglFNOGQZwzTXXHJk4cWLltm3bojIzM6uefPLJ7jfccEPs4MGDq37/+98fAbjvvvsO3nTTTZnz5893T5gwobSp8/7xj388OGvWrIFJSUmuSZMmle7fvz8a4Iorrii+8sorB3700UfdHn/88f31j3nhhRf23Hzzzf1uu+02S9++fatfe+21vYH4jP4k0+NvYgwsGzAYOAfoDazUdX2UUqq4fiGl1HOAb348Iz8/v10XczqdtPfYUOjerRsWoConh6PeuH86YI7k7e3o3abPcknfS3iq6Cle/elVRsSNaNN1n322G2Bn5swyioub/D4fJ9J+tpEUbyTFCpEVbyTFCh2LNyPDr7E1fqmoqPixqe0PPPDAoQceeOBQ4+02m4333ntvT+PtkydPLtu7d+/Gxttvu+22AqAAYM6cOcVz5swpblxm9OjR1du3b99c/1y+9QkTJlSuW7dua+Njvv32222+9Z49e7pycnI2NPkBm9HmZl6l1D6l1D4gGhjhe+/d1pocoE+997292+o7ALyvlKpVSu0BtmMmVwEY3r9A6zfz7j66G4CB3Qa26Ry+W2SW7lnaplG9JSUaS5ea/fmzZkkTrxBCNKfNNVNd11MABZwLGN4BQzuAh5RS97Vy+FpgsK7r/TGT6NXA7EZl3gVmAS/ouu7EbPbd3db4TnhNPNN0X6n5d0y/hH5tOsUo5yj6xPchuyyb7w59x9j0sS2WX7rUQVWVxvjx1fTp0/Z+ViGEaIuhQ4fW7NixY1Oo4wgEfwYgLQDOA2oAzVt7XA38rLUDvbfWzAM+AbaYm9QmXdfn67o+3VvsE6BA1/XNwHLgTqVUp85YEUka10zdHjc5ZWblvk9Cn2aPq0/TNKYOMGunS/YsabX8W2+ZA4+uuEKeWSqEEC3xp890MvAxsA24zbttM3BVWw5WSi0Fljbadl+9dQO4w7uIxhrVTPMq8qj11OKMcRJrj23zaab2n8qz659l6Z6lPHDGA1i0pv+eysmx8s030TgcBlOnSjIVQoiW+FMzjQEaP7bEiVlTFZ3Nl0yLi8HjIbvUHBzd1lqpz8lpJ5MRl0FueW7dVIRNWbzY7Cu98MIqEhPlUWtCCNESf5LpemAa3ltadF1/FLOJd10nxCUas9vxJCSgeTxoJSXsLzVHe7e1v9RH0zSm9J8CwJLdTTf1Gga8/baZTK+4QgYeCSFEa/xJpvdijuQ9A/Pe0jsAD/BA4MMSTfHUe65pe2umANP6TwPMUb2GcXytc9MmG9u320lJcXPOOdUdiFgIIULv7rvvTh80aNCIIUOGDM/Kyhq+bNmyuNaP8k+b+0yVUl/pun4acDPQD9gLPKeUWh/ooETTPMnJsG8flqKiusFHveJ7+X2eU3ucSo/YHmSXZbMhfwOj00Y32L94sdkHO2NGJXa/Jo4UQojw8vnnn8d98skn3TZs2LA5JibGyM3NtVVXV/s1VWBb+DWPoVJqI3An8FvgLkmkweVJNmfUshQVkVeeB0DPuJ5+n8eiWbgk8xIAPtz7YYN9hgEffOAAzGQqhBCRLCcnx56SkuKKiYkxwJyQITMzszbQ1/HnPtNkzJmHLq+37R3gRqVUYbMHioCpS6aFheRiTiXYnmQK5qjeRZsXsWT3Eu457Z66uX1//NFOTo6N9HQ3p54a8O+bEKIr07RTO+W8hvF9c7suvfTSkj//+c8ZmZmZIydNmlQya9asQt9TYwLJn5rpv4ArMPtLfcvlwPOBDko0rUHNtKL9NVOAcenjSHWksrdkL1sKt9RtX7LEHHg0dWolFnl0vBAiwiUlJXk2bty4+amnntqXlpbmuvbaawc+8cQTqa0f6R9/7jO9ANgDXAZsBYYBi73bRRD4kmlF0SFKUkqItkaTHN3kwxRaZbVYuSTzEl7e+jIf7vmQ4anDMQxYssRs4v3Zz6SJVwgRYC3UIDuTzWZj2rRppdOmTSsdPXp05b///e9U7xy/AeNP3WMXsFwptV4pVaOUWoc5U9HOQAYkmudLpnllBwFIj01v9dFrLfHN1fvhHrPfVJp4hRAnmnXr1kVv2LAh2vf+xx9/jOndu3fA50dosWaq6/pZ9d6+BNyv6/pPHKuZXoHcGhM0dcm00nzUa3pceofONz5jPN2iu7GjeAfbi7bzwQenAdLEK4Q4cZSUlFhvu+22viUlJVar1WpkZmZWv/jii215QItfWmvmXQE0vhFxYb11DXgUeDyAMYlmGN5kmltrtk50NJnaLXYm95vM69tfZ8nuD1myZBIgTbxCiBPHmWeeWfHjjz8e98i1QGstme7n+GQqQsRXMz3oOQq0f/BRfVMHTOX17a/z9n/2cfCgNPEKIUR7tJhMlVKZQYpDtIEvmeZazFHd6bEdq5kCTMqYRGJUInu/NSduuPjiKmniFUIIP/kzmhcAXddTgQZTMSml9gcsItEsXzLNiaoCOt7MCxBljeLCvhfy9rYZgJlMhRBC+MefSRsuxLzXtPH8dYY/5xHtZ8TGYkRFkRNnDkQLRDMvwNjoq3n78CgsjlLGj5e5eIUQwl/+NOg9C/Sm4aQNmp/nEB2haXhSUshNMN8GopkXoHTDOQB4Bn5IXpU0MgghhL/8qVGmAp8AM5VSAZ+KSbSNK7kbh+LM2Y+cMc6AnPOLz7zZeej7LN07iJtG3xSQ8wohRFfhT63ySaAv0EvX9YDPuC/apiAtHpcVEi2xOGyODp+vqEjj22+jsFg9MHhps884FUKISBUbG3tyZ1/Dn5rp28A8YDOAruu+7YZSSvpMgyTPac6dm6bFB+R8X3zhwO3WmDCxkp8SavnxyI/sK9lHv0T/HjouhBBdmT8101eAJKTPNKTyUqIA6O6JDcj5Pv3UrN1eMrmWyZmTAXh317sBObcQQnQV/tQo+wJrgbuA4s4JR7TmUKIVgO61UR0+V20trFxpTll5wQVVZGqX8c7Od1i8czG3nXRbh+ZrWYClAAAgAElEQVT9FUKIxrT/6ZxHsBn3h2YC/fr8SabPAWOBb5RSAZ8kWLRNnrd1t0dVx1vWf/ghitJSCwMH1tK3r5uenjNJcaSwo3gHmwo3MTJ1ZIevIYQQXYE/v5HPB0YCR3Rd3wu4vdsNpVTnPPBVHOewwwU1kB6A8dTLl5u10nPPNe8ttVvsTB8wnUWbF/HuznclmQohAiocapCdxZ/+ztHe8gnAKOCkeosIksN2s1GgR4m7lZKta5xMAS4ddClg9pt6DE+HryGEEF2BPzXTX3RaFKLNDlvNJ7r0KOrYZPSHD1vYuDEKh8PDGWccS6andT+NPvF9yC7LZk3eGsb3HN+h6wghRFfQ5mSqlHqxMwMRbXPY8E5yf6Rjj0lbscKslU6YUIOj3u2qmqZx6aBLefKnJ1m8c7EkUyFExKuoqPixs6/hz9y8/9fMLkMpdX2A4hGtyHeZA6l7Hu5Yp6kvmdZv4vW5bOBlPPnTkyzZvYT54+cHZHIIIYQ4kfnTzHsd5qT2vvslfOsGIMk0CAzD4Eh1IQDpuaUUu91gtfp9Hrcb/vMfM0Gec87xT4kZmjKU0c7RrM9fzyf7PmHGwBkdC1wIIU5w/gxA+h9gvvf1QeB9zETaXI1VBFhxdTG1nloSqyG2FixHj7brPD/9ZKe42EJmposBA5oeyHTVkKsAeGPbG+2OVwghugp/+kz/p/E2Xdef4/hHsolOkl+ZD/juMXWhFRZCSorf51m+3Fcrbf5xazMGzmD+mvmszFlJTlkOveLln1kIIZrjT59p30abEoGhwJg2Hj8ZWAhYgeeVUn9pptwVwFvA6Uqp79oaX1dwpPIIAN1d0YALa1ER7blBxtdf2lQTr0+yI5nJmZN5b9d7vLn9TX5zym/acSUhhOga/Gnm3dNoWQdMAva2dqCu61bgaeASYDgwS9f14U2USwBuB9b4EVeX4UumaYY5L69WVOT3OUpKNNats2OzGUyY0PJEVr6mXrVdyT2nQgjRAn+SaeMJ7iuBrzEHJrVmLLBTKbXbOxXh60BTo1r+BDwMNF9l6sLqmnm1RAAs7Uimq1dH4fFonHxyDXFxRotlJ2VMIiMug32l+1idu9r/gIUQIgxYrdZTs7Kyhg8ePHjEeeedNyg/P9//kZut8KfPtCNPh+kFZNd7fwAYV7+AruunAH2UUh/qun5nB651wqqrmdq7Ae1Lpl99ZTbxTprU+vTKVouVmUNmsvDHhby+7XUmZEzw+3pCCBFq0dHRnq1bt24GuPzyyzMXLFiQ9vDDD+cF8hqtJlNd11vrluvw80x1XbcAf6MNtVxd128AbgBQSuF0Ott1TZvN1u5jQ8Fms1HmnbAhIzED+J74qipi/PwMq1eb/1RTpjhwOqNbLX/juBtZ+ONClu5dytNxT5Mck9ymWCPtZxsp8UZSrBBZ8UZSrBB58YaLM844o3z9+vUxgT5vW5JgIJ7DlQP0qfe+t3ebTwLmJPorvA8dTwfe13V9euNBSEqp5zCfYANg5Ofntysgp9NJe48NBafTyf7C/QAk2cyaadXBgxz14zMcOWJh06Z0HA4PAwceoS2HJpHEWb3OYmXOSp5Z/Qw3jLqhTbFG2s82UuKNpFghsuKNpFihY/FmZGQEOJq20TQ65xFsBm2aQN/lcrF8+fKE66+/PuD/0G1Jpic3eh8PzAN0zET7UxvOsRYYrOt6f8wkejUw27dTKXUUqPsTS9f1FcDvZTRvQ74+09TEnoD/zbxff20+A3Xs2BqiW6+U1rl2+LWszFnJi5tf5Jcjf4lFk+fBCyEiR3V1tSUrK2v4oUOH7AMHDqy69NJLSwJ9jVaTqVJqHYCu6w7gFuBOoDuwAfgfpdQ7bTiHS9f1ecAnmLfG/J9SapOu6/OB75RS73fgM3QZdX2myWYl399kumpV2/tL67ug7wVkxGWwt2QvX+Z8ydm9z/breCGEgLbXIAPN12daWlpqOeeccwb/5S9/6X7vvfceDuQ12tJnGg3cDNwF9AC2ALcqpd7050JKqaXA0kbb7mum7Dn+nLsrMAzjWM3UmQn4n0x9g48mTmx+soam2Cw25gybwyPfPcKizYskmQohIlJCQoLniSee2D9z5sxBd99992G73R6wc7elvW438FfM2qgC7gWqdV2f7lsCFo1o1tHqo9R4aoi3xxPtTAf8S6bZ2Vb27bORlORh1Cj/H982e+hs7BY7n+//nAOlB/w+XgghwsHEiRMrs7KyKp977jn/p49rQVv6THtybFJ73bvUZ7TxPKIDDpUfAsAZ48STbI6otRQVgWGA1voYsVWrzP7S8eOr2zM3PmmxaUztP5V3d73Lv7f+mz+c/gf/TyKEECHQ+BFsy5Yt2xnoa7SlZrq/lSW7+UNFoBwuN5v302LSICYGj8OBVlODVlHRpuOPNfH6119a33XDrwPg1a2vUunq2PNUhRDiRNKWAUiZQYhDtMJXM02LSQPASE6G3FwshYW44+JaPNYw4Jtv2tdfWt9pPU5jZOpINhZs5O0dbzNn2Jx2n0sIIU4kco9DhKhLprFmMvV4nxbTln7T7GwreXlWkpPdDBniancMmqZx0+ibAPjHhn/IfL1CCOElyTRCNGjmhYb9pq349luzv/T002va0r3aomkDptErvhe7j+7ms32fdexkQohI176HKkeuZj+vJNMI0biZtz3JdOzY9veX+tgtdn418lcAPLP+mQ6fTwgR0f6brpNQj2J+3ibJKNwIkVdmzsncPbY7cCyZtuUxbIFMpgCzs2bz2A+PsfbQWr479B2n9TgtIOcVQkScr4HzQx1EOJCaaYQ4rpk3NRUAaytzcxYWWtixw47D0b77S5sSZ4/jmuHXAPDs+mcDck4hhIhkkkwjRONmXrf3aRGWI0daPG7tWrNWevLJtURFBS6e60dcT5Qlio/3fsz2ou2BO7EQQkQgSaYRwGN4GkzaAOBJM5OqpZWaaaCbeH26x3bn6qFXY2Dw+I+PB/TcQggRaSSZRoDi6mJqPbUkRSXhsDmAY8nU2krNdM2azkmmAPNOmkeUJYr3d73PtsJtAT+/EEJECkmmEeBIhZkwfbVSqNfM20LNtLJSY8MGOxaLwamnBj6Z9orvxaysWVI7FUJ0eZJMI4Dv0Wu+kbzQtmbeH3+043JpDB9eS0KC0SmxzRtj1k4/2P2B1E6FEF2WJNMIUPccU+/gIwAjPh7D4cBSUdHs/Lyd2cTrkxGfweys2RgYPPbjY512HSGECGeSTCPA4QrvbTGxx5IpmtbqiF7fSN7OTKYAt4y5pa52ujF/Y6deSwghwpEk0wjQVM0U6jX1NpFMXS747rvgJNOM+AyuG3EdAPPXzMcwOqdJWQghwpUk0wjgq5l2j+neYLvHWzNtauKGLVvslJdbyMx00aNH509If9tJt5EUlcSqg6v4eNfHnX49IYQIJ5JMI0B+pZksGzTzAu4Waqb1J7cPhmRHMreffDsA/2/5/8Plaf/TaYQQItJIMo0Ahytbrpk2NaI3GIOPGrtuxHX0TejL5vzNqO0qaNcVQohQk2QaAXx9pvXvM4XmJ24wjPqDj9r/MHB/RVujuef0ewBY8N0CymrKgnZtIYQIJUmmYc7lcVFQWYCGRmpMaoN9zY3m3bfPyuHDVlJS3Awc6A5arADTB0xnbMZYDlce5m8//C2o1xZCiFCRZBrmCqoKMDBwxjqxW+wN9nm6m82+jZNp/Sbejj4M3F+aprHwooVoaDy/8Xm2FG4JbgBCCBECkkzDXN3sR3Hdj9vnTk8HwJqX12B7sO4vbc4pPU/h2uHX4jbc/OGrP+AxOn80sRBChJIk0zCXV24myp7xPY/b5+7RA/AmU8+xhLVmTTQQumQKcNdpd+GMcbL20Fre3PFmyOIQQohgkGQa5nLLcwHondD7+J0xMbiTk9Fqa7EUFACQn29h924bMTEeRo4MzMPA2yMpOon7xt0HwINrHqSwqjBksQghRGeTZBrmfDXTjISMJvd7epo1VmuumXR9TbynnFKL3d7kIUFz+aDLmdBzAoVVhfz31/8d2mCEEKITSTINc76aaa+EXk3u9/WbWrz9pqG4v7Q5mqbx6FmPEmuL5d1d7/Lhng9DHZIQQnQKSaZhrrWaqdtXMz14EAjN/aUt6ZfYjz+O+yMAf/jqDxRUFoQ4IiGECDxJpmGuxT5T6iXT3FwqKsyHgVutBqecErr+0sbmDpvLxIyJFFQVcM+qe2QifCHECUeSaZjLq2ilZpphbrfm5fH993bcbo0RI2qJjw+fhGXRLPz1rL8SZ49j6Z6lqB0y1aAQ4sRiC9aFdF2fDCwErMDzSqm/NNp/B/BLwAUcAf5LKbUvWPGFo/LackpqSoi2RpMak0pBxfFNpB7fvaa5uSG/v7QlfRL68Kfxf+KOlXfwx1V/5JS0UxicPDjUYQkhREAEpWaq67oVeBq4BBgOzNJ1fXijYj8CpymlRgNvAY8EI7Zw5mviTY9NR2tmKqP6zbzhcH9pS/QhOpcPupxKVyU3L7uZSldlqEMSQoiACFYz71hgp1Jqt1KqBngdmFG/gFJquVKqwvt2NdB0J2EX4kumPeOOn7DBx5dM3QeP8MMP5r0w4ZpMNU3jzxP/TP/E/mwp3MID3zwQ6pCEECIggtXM2wvIrvf+ADCuhfLXAx81tUPX9RuAGwCUUjidzqaKtcpms7X72GApzy0HoF9Kv+bjdToxEhLYUDqECiwMGmQwbFhKkCNtqKWfrRMnb1z5Bme+eCYvb32Zcwedy5xRc4IcYUOR8F3wiaRYIbLijaRYIfLiPdEFrc+0rXRdnwOcBpzd1H6l1HPAc963Rn4Tz/JsC6fTSXuPDZbtedsBSLGl4HK5mo03rU8fvtx8JgCnnlpJfn5x0GJsSms/217WXswfP5+7v7qbm5feTJoljVN7nBrECBuKhO+CTyTFCpEVbyTFCh2LNyOj6QGNov2C1cybA/Sp9763d1sDuq5fAPwRmK6UCo8bJUPoQNkBADLiWv7iu/r140vMZBou95e2Zs6wOcwdNpcaTw2//OyXdU3aQggRiYKVTNcCg3Vd76/rehRwNfB+/QK6rp8M/AMzkR4OUlxhbX/pfgD6JvZtsZyr77FkOm5cePaXNmX+hPmM7zmew5WHuf7T62VAkhAiYgUlmSqlXMA84BNgi7lJbdJ1fb6u69O9xRYA8cCbuq7/pOv6+82crsvILjW7mfsl9Gux3ObYU8gnjXRHIZmZwX0YeEfYLXaeu+A5+sT3YV3+On697Ne4PK5QhyWEEH4LWp+pUmopsLTRtvvqrV8QrFgigdvj5kCp2czb3OxHPl+VnATAmXHfo2kjOj22QEpxpPDS5Je47IPL+HTfp/zhqz/wyJmPNHsrkBBChCOZASlM5Zbn4jJc9IjtQYwtpsWyX2dnAnBW7bIgRBZ4Q5KHsOjiRTisDl7d9iqPfv9oqEMSQgi/SDINU/tKzcmf+iT0abGcYcDX65MBOLtkCVRVdXpsneH0HqfzzPnPYNWsPP7j4zy7/tlQhySEEG0myTRM+fpL+ya0PPgoO9tKXp6NFEsRI9iELTu7xfLh7KJ+F7HgrAUA/GnNn/jH+n+EOCIhhGgbSaZhal+JWTPtl9jy4KPVq835eCckb8SCgXXv3s4OrVNdNeQqFpxpJtT5a+ZLQhVCRARJpmHKVzNtrZnX9zDwiQPNwUr2HTs6N7AgmJ01u0FCffyHx+WxbUKIsCbJNEz5+kxbuy1m9WpzcvszvPeX2rZs6dzAgsSXUDU0Fny/gHu/vhe3J3Ju+xFCdC2STMOQYRjsProbgMzEzGbL5eVZ2LvXRny8h6HnmvPx2rduDUaIQTE7azbPnv8sUZYoFm1exK+X/Zpqd2TM8CSE6FokmYah/Mp8iquLSbAn0CO2R7Plvv7arJWefnoNjBgKgG3nTqitDUqcwTBtwDReueQVEuwJLNmzhCuXXMmhikOhDksIIRqQZBqGdhSb/Z6Dkge1OHnBl1+ayXTSpGqM+Hhcffui1dRg27MnKHEGy4SMCbz9s7fJiMvgh8M/MOXdKfx05KdQhyWEEHUkmYYhXzId3G1ws2UMA1auNJPpWWeZTZ+1WVnAidNvWt+I1BEsvXQpY3uMJa88j8s/uJw3t78Z6rCEEAKQZBqWdhbvBFpOpjt32sjLs+J0usnKMuezdXmTqf0ETKYAabFpvDH1DX6e9XOq3dX85j+/4fYVt1NWUxbq0IQQXZwk0zC0tdAcRNRSMvXVSs88sxqL91+xdtgwAOwbN3ZugCEUZY3i4UkPs+DMBTisDt7a8RaTF09m/ZH1oQ5NCNGFSTINM4ZhsKlgE2A2bTanfjL1qTntNACivv8ePJ5OjDK0NE1jdtZsPrrsI4alDGNPyR6mvz+dx394nBp35DyCTghx4pBkGmYOlB3gaM1RUhwp9Izr2WSZmhr45htzsgZffymAJyMDV69eWEpKsG3bFpR4Q2lI8hCWzFjCdcOvo9ZTy4LvFzDl3SmsO7Iu1KEJIboYSaZhZmO+2UQ7KnVUsyN5V6+OprzcQlZWLT17NqyB1owdC0DUt992bqBhwmFz8L8T/5c3prxBv4R+bCncwrT3pnH/N/dTUlMS6vCEEF2EJNMws7HATKYjnSObLfPJJw4ALrro+CfE1DX1fvddJ0QXvib1msQXV37BjaNuBOD5jc8z6Y1JvLr1VZk5SQjR6SSZhpnvDplJcEzamCb3GwZ8+qnZX9pkMvXVTNesMQt3ITG2GO474766W2gKqgq488s7mfLuFFZkr5D5fYUQnUaSaRip9dTyw+EfAPP5nk3ZtMnGwYM2und3M2bM8TMduYYOxZ2cjC0nB9uuXZ0ab7ga5RzFOz97h7+f93cy4jLYWLCRn3/8cy7/4HK+Pvh1qMMTQpyAJJmGkc0Fm6lwVZCZmEn32O5Nlnn//RjArJVamvrXs1qpPu88AKI//7yzQg17mqYxY+AMVuoruXfsvSRHJ/PtoW+Z+eFMZi6ZyX8O/EdqqkKIgJFkGka+yf0GgLHpY5vc7/HA4sVmMr3ssspmz1N1/vkAOD79NMARRp4YWww3j7mZb67+ht+f+nsSoxL5OvdrZn80mwvfuRC1XcntNEKIDpNkGkaWZy8H4OxeZze5f+3aKA4etJGR4WLs2OYTQPW552I4HESvWYM1J6dTYo00CVEJ/PaU3/LN1d9wz+n30D2mO1sKt/Db//yWwU8P5s9r/1z3QHYhhPCXJNMwUVZTxpq8NVg0C2f1PqvJMm++eaxW2mQTr5eRmEjVhRcCELN4ccBjjWTdortx60m3snrWav529t/ISs4irzyPp356iglvTOCqD6/ivV3vUelqvuYvhBCNSTINE8sPLKfWU8tJaSeR4kg5bn9h4bEmXl2vaPV8FVdeCUDsq6+CW24NaSzaGs1VQ67i8ys+54s5X3Dl4CtxWB18dfArfr3s14x5eQzzls3j032fyjNUhRCtsoU6AGFavNOsQf5swM+a3P/CCxaqqiycc04Vgwa1nhyrzz0XV79+2Pbtw/HJJ1RNmRLQeE8UmqYxqc8ksmKymD9+Pot3LkZtV6zLX8fiXYtZvGsxiVGJXNTvIi7seyFn9T6LxKjEUIcthAgzUjMNA4VVhSzLXoZFs3DpwEuP219VBc88YwXgF78ob9tJrVbKf/lLABIefRRcroDFe6JKik7iuhHXsfSypay6ahX3nH4Pw1OGU1JTwls73uLGL25k1EujuHLJlTyz7hk2F2zGY5y4cyALIdpOaqZh4KXNL1HrqeX8Puc3eUvMSy/FkZ2tkZVVy3nntb3JsfznPyfun//Evm0bsa++SsXcuYEM+4SWmZjJrSfdyq0n3crO4p18vPdjlmcvZ+2htXyT+w3f5H7Dg98+SHJ0MuPSx3FGzzMY33M8w1KGYbVYQx2+ECLItAi/1844ePBguw50Op3k5+cHOBz/ldeWM+GNCeRX5vPGlDeY1GtSg/1Hj2pMnNidoiIrixYVcOGF/vXfOZYsIeXGG/EkJXHk449x9+0byPCbFC4/27byJ97i6mJWHljJF9lf8NXBr8grz2uwP8GewOi00ZyUdhJj0sZwUtpJZMRlNDvPcmfGGg4iKd5IihU6Fm9GRgZAYL6UApCaacg9ve5p8ivzOTntZCZmTDxu/733JlFUZGXSJA8XXOD/QJiqqVOpvPhiYj75hOQbb6Tg7bcxYmMDEXqX1C26G9MHTmf6wOkYhsG+0n2szl3NN7nfsDp3NQfKDrDq4CpWHVxVd4wzxsnI1JEMSR7C0OShDEkewpBuQ4iPig/hJxFCBJIk0xBaf2Q9f1/3dwDuP+P+42ovb74ZwzvvxBIT4+GZZ1y0q3KjaRQ/9hj2yZOJWr+elDlzKHzxRYyEhAB8gq5N0zQyEzPJTMzk6qFXA5Bbnsv6I+v56chPrM83X/Mr81lxYAUrDqxocHzv+N4MSR5CZmIm/RL70S+hH/0S+9EnoQ8xtpgQfCIhRHtJMg2RnLIcfvn5L6n11HLd8Os4Pb3hXLyffhrN73/fDYD77y9hyJBY2tsCZSQlUfDKKzh1neg1a0ibPJmihQup9T5hRgROz7ie9IzrycWZFwPmw973l+5na+FWthZtZXvRdrYVbWNX8S4OlB3gQNmBJs+THpdOv4R+9IzrSXpcOumx6aTHpTOsahgxrhi6x3bHbrEH86MJIVoQtGSq6/pkYCFgBZ5XSv2l0f5o4CXgVKAAuEoptTdY8QXT2ry13PjFjRyqOMSp3U/l3nH31u2rrYW//z2eRx9NwOPRuOWWUq65pgLoWNOse8AA8t95h5Rf/hL7pk2kzZhB5ZQplP/Xf1Fzxhm0r9orWqNpmlnrTOxXl2ABXB4Xe0v2sr1oO/tK9rGvdB/7Svaxv3Q/2aXZ5JXnHdcf2+C8aKTGpJLqSCXFkUKKI6VuvW5bjLk90Z5IQlQCCVEJ2Czy97MQnSEo/7N0XbcCTwMXAgeAtbquv6+U2lyv2PVAkVJqkK7rVwMPA1cFI75gcHvcfHfoO17c8iLv7XoPgDPSz+D5C58nxhZDQYGFJUscPP98PLt3m/8sv/1tKb/7XWngYujblyPvv0/C448T/89/ErN0KTFLl+Lu2ZPqs86ietw4XMOGUTt4MMRIM2NnsllsDOo2iEHdBh23z+VxkVOWw/7S/WZSrTATa255LvnV+RwoOcDhisPkV+aTX+lfc0WsLZbEqGPJ1beeGJVInD2OGFtM3RJri23yvcPmqNsWbY0myhqFTbMFbJCVEJEoWH+mjgV2KqV2A+i6/jowA6ifTGcAD3jX3wKe0nVdU0oFfLhxTn4Zh8ugoLAQj2Hg8XgXDAyDY+8NAwPwGB48nmP7DLxlDA8eD3XHud0eKl2VlNVUUFZbTkl1KblleRwozWFL4VbKK2uhOglbzWzOdE5jVNEl/PGjKLZvt7Nly7Emu/79XTz44FHOOacTZt5xOCi95x7K584l7pVXiH39day5ucS+8Qaxb7wBgKFpeJxO3D164OnRA09qKp74eIy4OIz4eDxxceYgJrsdw2Y77lVzOrGXlYHNZtZ4vYvh+2VrsTTYflyZlhZ/tLV8RQWWoqLOOXc7ykcB/bVo+lsHQ+JgqDdHREpKCoWFhbg8LvJriiisKaaguojC2mIKaorM9ZpiCmqKKawporC6mBJXGaW1ZZS6yqlwVVDhqiCvovlab3tFW6KwW+xE1S1ROOwObIbFuz2q3j47dosdq2bBqtmwahZsmrXBukWzYrNYsWrmYm4z99k0W926VbNitXj3Y8GiaYCGRbOgoWHRNO9rvX1Y0DQNC5r31UJCQgLlZeVNHgc02G7B0uQ+oO6PCvNKjbc1fG28r63n0NDw1GRiiUoO+L+jaJ9gJdNeQHa99weAcc2VUUq5dF0/CqQCAR+rfsmcIgo2DAEyAn3qNnEBy72Lj8NhMG5cNbpeyZQplURFdW4MnowMSu+8k9Lf/Q7b5s1Ef/klUevWYdu6Fdvu3ViPHMF65Ahs3Niu86cFON7Olh7qAPzgi7W3n8d5NCiLgqPRcNRx/Gu5HSrsUB5lvvqW8vrrjfZVW6HGCi4rVHtqqPbIE3iCRT/ah8d+vzrUYQiviOtA0XX9BuAGAKUUTqfT73PEJuyhMM7M0ZpmgGYAvtfG20DzvoJhVjC8+5pat2gWrBYNq8WC1WIhymonyhZFrD2GhFg7iYmQmAhJSQbp6TBokMGgQQYnn2zgcFiAOO/SkM1ma9dnbZPu3eGcczA/IdS6XHDoEFpuLhw8iFZQAKWlUF6OVlpqrldUQG0tWm2t2dHrW1wuNJcLo6bGnBPYMJpePB4wDLQW9tUtcOy1rfwp35nnDqPyGpAAJBiGmYhrvEsAuDWDWgtUWw1q6l4Nb7L1vloMc7sVqi3mq1szcGvgsoDL4l3XwG0Bl2bgsuDdX6+cZnj3m8f7ytRaDAwNDA08eF81MDC8r+b7+uuGZjQq27iMUa9sw/N6MOqte/8pvBVOo9G6uc84tq8t5Zs8h/l5ANKsCZ33O0H4LVjJNAfoU+99b++2psoc0HXdBiRhDkRqQCn1HPCc963RnpuWV7/ZH6czsQ03PGuNXtvD5V1anpy+rMxcmhP0G8qjoyEz01z81JVufg+2SInVCsQAfSIkXoicn61PACZtEAEUrGS6Fhis63p/zKR5NTC7UZn3gWuBb4ArgWWd0V8qhBBCBFpQJrpXSrmAecAnwBZzk9qk6/p8Xdene4v9C0jVdX0ncAdwTzBiE0IIITpK5uaNEJEUbyTFCpEVbyTFCpEVbyTFCjI3b7iRR7AJIYQQHSTJVAghhOggSaZCCCFEB0kyFUIIITpIkqkQQgjRQRE/mjfUAQghRISS0bwBFOk1U629i67r33fk+GAvkRRvJMUaafFGUqyRFm8kxRqgeEUARXoyFUIIIUJOkqkQQgjRQV05mT7XepGwEknxRlKsEAVgCt4AAAYfSURBVFnxRlKsEFnxRlKsEHnxntAifQCSEEIIEXJduWYqhBBCBETEPRw8EHRdnwwsxHzs4vNKqb+EOKQGdF3/P2AacFgpNdK7LQV4A8gE9gK6UqooVDH66LreB3gJ6IF5q9JzSqmF4RivrusOYCUQjfndf0spdb/30YCvA6nA98A1SqkAPTa7Y3RdtwLfATlKqWlhHuteoBRwAy6l1Gnh+D3w0XW9G/A8MBLzu/tfwDbCLF5d14d6Y/IZANyH+f8urGLtyrpczdT7y+lp4BJgODBL1/XhoY3qOIuAyY223QN8oZQaDHxB+DyizgX8Tik1HDgDuMX78wzHeKuB85RSY4CTgMm6rp8BPAw8ppQaBBQB14cwxsZux3xsoU84xwpwrlLqJKXUad734fg98FkIfKyUygLGYP6cwy5epdQ278/0JOBUoAJYTBjG2pV1uWQKjAV2KqV2e/+ifx2YEeKYGlBKrQQKG22eAbzoXX8RuDSoQTVDKZWrlPrBu16K+QupF2EYr1LKUEqVed/avYsBnAe85d0eFrEC6LreG5iKWXtC13WNMI21BWH3PQDQdT0JOAvzOcoopWqUUsWEabz1nA/sUkrtI/xj7VK6YjNvLyC73vsDwLgQxeKPHkqpXO96HmazaljRdT0TOBlYQ5jG622Z+B4YhNlCsQso9j7AHszvQ68QhdfY48BdQIL3fSrhGyuYf5h8quu6AfxDKfUcYfo9APoDR4AXdF0fg/mduJ3wjdfnauA173q4x9qldMWaacRTShmE2VSKuq7HA28Dv1FKldTfF07xKqXc3uay3pitFFkhDqlJuq77+sy/D3UsfpiklDoFswvlFl3Xz6q/M5y+B5gViVOAZ5RSJwPlNGomDbN40XU9CpgOvNl4X7jF2hV1xWSaA/Sp9763d1u4O6Trek8A7+vhEMdTR9d1O2YifUUp9c7/b+9uQrSsogCO/2sRNCKFKBHKZIS0aCVhICJIZIGBi4hDKSnhIoIiaBO1ScGFuBBqnxiUySGQJMI2o5vAz4zJj1zYh02RgTAR7RRc3Ps6I0gWz+s8d3j/v83zMe/ivHBnzjP33uecervZeAHqlN4RYDXwYEQMZmlaGQ9rgI11U88ByvTuB7QZKwCZ+Vs9/klZ03uKdsfBFDCVmcfr9eeU5NpqvFAeUr7NzCv1uuVYR84oJtOTwIqIeLQ+6b0EHOo5pv/iELC1nm8FvugxlpvqOt5HwIXM3DPrR83FGxFL6g5OIuJ+YD1ljfcI8GL9WBOxZua7mbksM5dTxuhEZm6mwVgBImJBRCwcnAPPAmdpcBwAZOYfwK91pyyUtcjzNBpv9TIzU7zQdqwjZ+TWTDPzWkS8AXxNeTVmb2ae6zmsW0TEZ8A6YHFETAHvA7uAjIhtwC9A9BfhLdYArwDfR8R39d57tBnvw8DHdd30XiAz88uIOA8ciIidwBnqppRGvUObsT4EHIwIKH9X9mfm4Yg4SXvjYOBN4NP6UP0j8Cp1XLQWb31AWQ+8Nut2i79jI8sKSJIkdTSK07ySJA2VyVSSpI5MppIkdWQylSSpI5OpJEkdjdyrMdIwRMQ+yrt9OzJze7/RSOqbyVS6jVp56JF/+chmYBo4NicBSWqayVS6vb3Aonr+OnAfpWTiVL13IjP39xGYpPZYtEG6g4iYBh6g9Oo8Wu/tY9Y0b0Rsp1SqOgxcBV6glNPbROn8sonSwHnLoGVdRIwDu4G1wALgBPB2Zp6do68maUjcgCQN13OUVmmXgVXAKUpbukngCeBDgIgYAyYoJeAmKXVW1wETEbF4zqOW1InJVBquS8AGSt1UgDHgGWZqqq6sx+eBx4DfgYuUZvCXgSXMFLKXNE+YTKXh+qH2lpyu11cy8y/g73o9Vo/L63EppSn1W5TkCqVxuaR5xA1I0nBdv8P1wM/1eBpYVRMwtUXcPXcnNEl3i8lU6sdXwE/Ak8A3ETEJjFPWTTcAR3uLTNL/5jSv1IPM/Ad4mtLseZyyM/hx4BPKGqqkecRXYyRJ6sj/TCVJ6shkKklSRyZTSZI6MplKktSRyVSSpI5MppIkdWQylSSpI5OpJEkdmUwlSeroBsYwXKUC3rlTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sline = plt.plot(\"t\",\"S\",\"\",data=sir_out,color=\"red\",linewidth=2)\n",
    "iline = plt.plot(\"t\",\"I\",\"\",data=sir_out,color=\"green\",linewidth=2)\n",
    "rline = plt.plot(\"t\",\"R\",\"\",data=sir_out,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
}