notebooks/kr08/3_4/r_desolve.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Program 3.4 from Keeling and Rohani \n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Author: Theresa Stocks @theresasophia\n",
    "\n",
    "Date: 2018-10-02\n",
    "\n",
    "*NOTE*: This does not exactly reproduce the figures from the Python code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: ‘reshape’\n",
      "\n",
      "The following objects are masked from ‘package:tidyr’:\n",
      "\n",
      "    expand, smiths\n",
      "\n",
      "\n",
      "Attaching package: ‘magrittr’\n",
      "\n",
      "The following object is masked from ‘package:tidyr’:\n",
      "\n",
      "    extract\n",
      "\n",
      "\n",
      "Attaching package: ‘plyr’\n",
      "\n",
      "The following objects are masked from ‘package:reshape’:\n",
      "\n",
      "    rename, round_any\n",
      "\n",
      "\n",
      "Attaching package: ‘dplyr’\n",
      "\n",
      "The following objects are masked from ‘package:plyr’:\n",
      "\n",
      "    arrange, count, desc, failwith, id, mutate, rename, summarise,\n",
      "    summarize\n",
      "\n",
      "The following object is masked from ‘package:reshape’:\n",
      "\n",
      "    rename\n",
      "\n",
      "The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Loading all necessary libraries \n",
    "library(deSolve)\n",
    "library(tidyr)\n",
    "library(reshape)\n",
    "library(magrittr)\n",
    "library(plyr)\n",
    "library(ggplot2)\n",
    "library(deSolve)\n",
    "library(tidyr)\n",
    "library(dplyr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m <- 4  #number of age classes\n",
    "mu <- c(0,0,0,1/(55*365)) #death rate in each age group; it is assumed that only adults die\n",
    "nu <- c(1/(55*365),0,0,0) #is the birth rate into the childhood class; it is assumed only adults give birth.\n",
    "n <- c(6,4,10,55)/75 # fraction in each age class (assumption that life expectancy is 75 years)\n",
    "S0 <- c(0.05,0.01,0.01,0.008) # inital value for number of susceptible\n",
    "E0 <- c(0.0001,0.0001,0.0001,0.0001) # inital value for number of exposed\n",
    "I0 <- c(0.0001,0.0001,0.0001,0.0001) # inital value for number of infectious\n",
    "R0 <- c(0.0298, 0.04313333, 0.12313333, 0.72513333) # inital value for number of recovered\n",
    "ND <- 365 # time to simulate\n",
    "beta <- matrix(c(2.089, 2.089, 2.086, 2.037, 2.089, 9.336, 2.086, 2.037, 2.086, 2.086, \n",
    "                 2.086, 2.037, 2.037, 2.037, 2.037, 2.037), nrow=4, ncol=4) # matrix of transmission rates\n",
    "gamma <- 1/5 # recovery rate\n",
    "sigma <- 1/8 # rate at which individuals move from the exposed to the infectious classes\n",
    "TS <- 1 # time step to simualte is days\n",
    "\n",
    "# combining parameter and initial values\n",
    "parms <- list(nu=nu, beta=beta, mu=mu, sigma=sigma, gamma=gamma)\n",
    "INPUT <- c(S0, E0, I0, R0)\n",
    "\n",
    "# constructing time vector\n",
    "t_start <- 0 # starting time\n",
    "t_end <- ND - 1 # ending time\n",
    "t_inc <- TS #time increment\n",
    "t_range <- seq(from= t_start, to=t_end+t_inc, by=t_inc) # vector with time steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Differential Equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# differential equations\n",
    "diff_eqs <- function(times, Y, parms){\n",
    "  dY <- numeric(length(Y))\n",
    "  with(parms,{\n",
    "    # creates an empty matrix               \n",
    "    for(i in 1:m){\n",
    "      dY[i] <- nu[i]*n[4] -  beta[,i]%*%Y[2*m + seq(1:m)] * Y[i] - mu[i] * Y[i] # S_i\n",
    "      dY[m+i] <-  beta[,i] %*% Y[2*m + seq(1:m)] *Y[i] - mu[i] * Y[m+i] - sigma * Y[m+i] #E_i\n",
    "      dY[2*m+i] <- sigma * Y[m+i] - gamma * Y[2*m + i] - mu[i] * Y[2*m+i] #I_i\n",
    "      dY[3*m+i] <- gamma * Y[2*m+i] - mu[i] * Y[3*m + i] #R_i\n",
    "    }\n",
    "    list(dY) \n",
    "  })\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ageing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "RES2=rep(0,17) #initalizing the result vector \n",
    "number_years <- 100 #set the number of years to simulate \n",
    "\n",
    "# initialize the loop\n",
    "k=1\n",
    "# yearly ageing\n",
    "for(k in 1:number_years) {\n",
    "  RES = lsoda(INPUT, t_range, diff_eqs, parms)\n",
    "  #taking the last entry as the the new input that then is propagated accoring to the aging\n",
    "  INPUT=RES[366,-1]\n",
    "  INPUT[16]=INPUT[16]+INPUT[15]/10\n",
    "  INPUT[15]=INPUT[15]+INPUT[14]/4-INPUT[15]/10\n",
    "  INPUT[14]=INPUT[14]+INPUT[13]/6-INPUT[14]/4\n",
    "  INPUT[13]=INPUT[13]-INPUT[13]/6\n",
    "\n",
    "  INPUT[12]=INPUT[12]+INPUT[11]/10\n",
    "  INPUT[11]=INPUT[11]+INPUT[10]/4-INPUT[11]/10\n",
    "  INPUT[10]=INPUT[10]+INPUT[9]/6-INPUT[10]/4\n",
    "  INPUT[9]=INPUT[9]-INPUT[9]/6\n",
    "  \n",
    "  INPUT[8]=INPUT[8]+INPUT[7]/10\n",
    "  INPUT[7]=INPUT[7]+INPUT[6]/4-INPUT[7]/10\n",
    "  INPUT[6]=INPUT[6]+INPUT[5]/6-INPUT[6]/4\n",
    "  INPUT[5]=INPUT[5]-INPUT[5]/6\n",
    "\n",
    "  INPUT[4]=INPUT[4]+INPUT[3]/10\n",
    "  INPUT[3]=INPUT[3]+INPUT[2]/4-INPUT[3]/10\n",
    "  INPUT[2]=INPUT[2]+INPUT[1]/6-INPUT[2]/4\n",
    "  INPUT[1]=INPUT[1]-INPUT[1]/6\n",
    "  RES2 <- rbind(RES2,RES)\n",
    "  k=k+1\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAJFBMVEUAAAAAv8QaGhozMzNN\nTU18rgDHfP/Z2dnr6+vy8vL4dm3///+evM4gAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElE\nQVR4nO2djZqqOrOEh/nEzYH7v9+zhIT8EGJoCux2qp69XUBDTQn9ThAd+Zkoijqtn08HoKhv\nEEGiKIAIEkUBRJAoCiCCRFEAESSKAoggURRABImiAAKC9B9FXS9cw0JFkChbwjUsVASJsiVc\nw0JFkChbwjUsVASJsiVcw0JFkChbwjUsVASJsiVcw0JFkChbwjUsVASJsiVcw0JFkChbwjUs\nVATpWv2+1Ljq/J+bovaEa1ioCNKlmiFqIemXIDUK17BQEaRL9Rs9tq9IkCrCNSxUBOlS/Ub/\nLoPTr3vw//z+5xf9+pmwJrUVrmGhIkiXyuPg8HDkRNMvcn7XBckMB6aicA0LFUG6Vu5iQwRS\nClX0Tz5DkIrCNSxUBOl6Ladt//mTOv9CKEcsnWm/2vfHhGtYqAjSpfIndusAE/iI/tkZkaiS\ncA0LFUG6VBlI68VwntqJhWtYqAjSpfLvIy0XEqILd/+ti3YvNvDUrihcw0JFkK6Vf63jXydt\nL3//JrDx8vc74RoWKoL0YRGXg8I1LFQE6cMiSAeFa1ioCNKHRZAOCtewUBEkypZwDQsVQaJs\nCdewUBEkypZwDQsVQaJsCdewUBEkypZwDQsVQaJsCdewUBEkypZwDQsVEKRhX9ViuyAuirIo\nimInC65hoSJIUhuEh54odrLgGhaqQyA9Ho90KiyYCJLAQ08UO1lgrY/VEZAe83/xVFgyESSB\nh54odrIgux+oAyA9Ng+PpC7dMwd2MMRETxZFUexkgbU+VidB4qndKQ89UexkQXY/UOdA8iT9\n76W3BqMgH0WZ0OlTu3aQRoJEfa0wIM16N1aP47A+Cod8+aZZFoANwkNPFDtZkN0P1N0gjQQp\n8tATxU4WZPcDJbj8/QhTydWGd3uGIG089ESxkwVNAEiCN2Qf69ShN2RnhEaCFHvoiWInC7L7\ngbrvI0IEaeuhJ4qdLLiGhYogSW0QHnqi2MmCa1ioCJLUBuGhJ4qdLLiGhepekEaClHjoiWIn\nC65hobobJL6PFHvoiWInC65hoboNpJEgbT30RLGTBdewUBEkqQ3CQ08UO1lwDQvVrSCNBCn1\n0BPFThZcw0J1M0j8rF3ioSeKnSy4hoXqLpBGglTw0BPFThZcw0J1H0gDQdp46IliJwuuYaEi\nSFIbhIeeKHay4BoWKoIktUF46IliJwuuYaG6F6SBICUeeqLYyYJrWKhuAmksgnSUqW/rGEVR\n7GTBNSxUt4E0bEE6PDh9W8coimInC65hoSJIUhuEh54odrLgGhaqe0AaCVLRQ08UO1lwDQsV\nEKSK5m/iGt1jujRd75Y0FAXXXSNSYTIaprarlX8ZvflN1iRFv3oVRbGTBdewUH0ApDGaGt+A\nlC74to5RFMVOFlzDQqUApHFbKmw378O241CXoo5RFMVOFlzDQnU7SEMGUn6ul25GkFpsICZW\nsuAaFiqCJLVBeOiJYicLrmGh+hxI27/0OwaS9LNGijpGURQ7WXANC9WnQdpcxks3Sz9PNO3V\njkhRxyiKYicLrmGh+gRICz0bkDZf1ZUsGGeQShct1iXln7c9EJWg7SJIZRuICUHaUQ7SuAtS\nyk42Wk3Fq3/DwtBebV2hkOVIbXMp/sTH2L2HleY9YAMxIUg7Shq0BlJgZzwA0vazR28hG/Pa\n+KY2ZLWpUksn92vTeYtXbRJud2WWExYEaV8FkFZo1t0aLxzGEkjpfFrbA6JcG/PaeLQ2Cbfb\nqZ2xmABPB1Wbdlc7/FR3hGtYqD4IUrpzI5DGeKV1owikCLIxqzUcNUxtB6TCpLTWHOtY856N\nVa9ts0if6p5wDQuVHpBWiMKCaDt/MujdxvxwNNbiBfmkstruapvmFW53RW3aXe3oU90TrmGh\n+ghIYxmkaMEYg+TZaABpDLXdA2il1sziTvN+5OnIsuRPla+RrgDJH4hxDDt7jI5MDNmRTiit\ndlNtbKk1209qIheyiJ4qQToI0lprA2mKQBrSwzGOUS0syGtDuVZe8OHa3mrfUHuzGkFqAcmT\nFIE0FkDy+3sM8zlImyNTqwVDI7WTLGp7OgeeKi9/t4HkR4dQG/2eThYE5MogpStFB2LIF6S1\n3SO61goHvlQbT9amSq3UhAdrdz8daW3zdL4epHAT83Wq7a7mTSD51VYQ4rFr1pQvi6HJQYpq\n2RFNPHdqQ1ttOFWbKrXdH31l7dzTkdY2sb4dpMfKzTr1uBykIezvIQIpryXcpH2RHktUbQJ4\nTkkT7fbcp2vI3df0o78cpMf24YEZkdze3C4Yx+EYSENTbe/gHqlNwu2ShqmuBo+so9bwVAnS\nDSCNoZaMVvlRGypHdChOli8g7S2YKrXadknDtDWa1N5KbfNU/xxI6xne/16qbLr54q3kS7mc\n8gXJ/PJlXuvkGL7iK65Fj+VJO7X91a6o6XqqFnUGpMd1IE35fOCmvODioz1dW7vYXnlts5pF\nnQPp8YiuNhw7tUtP88ZowZgtmOfnLSY/OXqLaIFfeq427k3mCybhdvHkVKmN7bVJuJ0k8rta\nnqVlu02sP3dqtz7Oqjz5KkjDurMjcHZBSlcOtSIs5SMKqU0AzwnSvEOleS/fDYgsm9W+HCT/\nkih+cXQRSAVY5n1YAqlwkIbKAUTV9kFq9BxykMSxJuF20NpQyTLUt/tzILm3YR/r1CQGabvC\nBqRCbdqsHB+qSu1AJ7yr+cmpUqttVwKp8Ucjm1ca+V1tEmy3We3rQaqr8uQPgpRuUQCpuF1j\nbQ4aHbXCdi21Upb9BeXJAkjRZK2W2O/9hey77a6oSbIQpFSVJ98A0v7KcpCKP0TRt3woimIn\nC65hoVICUk0bkM5JUccoimInC65hoTIB0roPT7isUtQxiqLYyYJrWKgMgBTtQ4iJno5RFMVO\nFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VO\nFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VO\nFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VO\nFlzDQkWQpDYIDz1R7GTBNSxUQJAqsvodSxTVKI5IUhuEh54odrLgGhYqgiS1QXjoiWInC65h\noSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65hoSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65h\noSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65hoSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65h\noSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65hoSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65h\noSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65hoSJIUhuEh54odrLgGhYqgiS1QXjoiWInC65h\noRLc1TyaCgsmgiTw0BPFThZY62N1BKTH/F809fonkFR58gSp7KEnip0s0PbH6QBIj83DY1rR\nmgiSwENPFDtZkN0P1BmQHiKQ+hM7WL5pZKKnYxRFsZMF2f1AnRqRJn9q97+XKptGX8fV9wcT\nUpQBnQTJvUQSgrSBipRRRnUOpEdcrwzH0ald3/fRZHqet1mQnQd+2zmMoih2ssBaH6vTFxuC\nKk9+H6SElAJZ6T5sOAhvpahjFEWxkwXW+lgJLn8/ov9iVZ58FaQ+rm2GqGQfpsXasdiXoo5R\nFMVOFmT3AyV4Q/bhpx6z1mrlye+BlA1KpSEqmk6P0uY0MCnWDkQlaLsIUtkGYvLtINVVefIA\nkPoZpHTlMO3sIudU0YKpUtsA2O9MvzrmxHV872GleQ/YQEwI0o52QRoykIYaSP3Ul0HqnVNW\n6+OfMgwrb1FtyGvRY702xbUh225oq01TZbs8y35tEm6XZjn/dHazvNsuj0WQ9rSANEPTAFLo\n96gag9SvbtGkW7fvo41Tn15SG3YWTJVabbt4csLEmgQ/+qqaJMvmqRKkKki9Vw5SvMCx5JZn\nIPXZLu+zydba247t39cm4Xbx5FSp9e21SbidJPK7WjlLfbtNLIL0BiRHQwkk1+/JvF87XZDs\n+L5eK0xuyJJt10/C7eLJ6byFttpUqTXbE6R3IK2crCBF+2/dqX2gruvCfl4xdGvGtWz7vut3\nauXDKKhNAM+pUitafGdtsxpBegtSBEp0DhcwGcI/r/+7CJZ/kwtYjo2um5ckNe+ylLreMZXU\n+nKtuOBIbThemyo1PymtnYh1d23zdAhSM0gLCwlIL0VszLUVlpmNRa+dPawgzXt/nvS1zte6\neq1zvHXhcUhrXV7brlasbTz3alOlVrbQVDv0VGu1zWoE6T1IQw2kbtml6YI+DEBd6OMuWhAm\n3dLNgjtqO6tJaxfba6ptVusJ0nuQ/FmbL7nzrRgkN9DHu/gFy7Dd+9R3iiAdB8lfGFhAcrD0\nJViGcHo3RI+hNuS18uQltcpqktqRpwP+0Z+obZ8qQToCkttvZZCGaBcve3iYQqPNnbY+LAui\nybSWTt5f21+tVmu2n77xqRKkdyB1K0idVziTi0edRW6y8yAlv7D6ZTIs2NT6TW0OGjzXyxrr\nY1RLlqa1qVKrbddeG1prk3C7hlqX14Zq5DXLuafKiw17ikCaufEgOWzCmdwWpJiN5Uf08Xlh\nmHxbcwvKQTcrvxM//V22gZgQpB1lIPUrSJ6k5FdVepocsTH5JUMoYo/SARuEh54odrLgGhaq\nT4DURSAN+QC01Qakc1LUMYqi2MmCa1io7gQpfk10AKRoHx44GJUDgXAhSDs2EBOCtKMxRqbr\n1vdXXZkgnbOBmFjJgmtYqIAgVfQPlNfDPOFeE8265YdT1PW6Z0SKTt+iiwvHf1NJNtqY6PnV\nqyiKnSy4hoXqdpD8nGgHJ3PPbPIZL3iuC9zkumCq1NIFz2Qyq02VWm27eHKaTlu8Hifhdklt\nOv90DmSp2xOkPXX+Q92vPXUEpBSWKYbl6Y6Bn3ymR+a5ruA1RI/55GaBtLZRrTadt/iG2mY1\ngrQL0tO3/ry39rl55vPPaHLdOjsyjQeJsiOCtAdIcS8t0+kg81xP1J5h1CmMLE2HIzwOlQXD\nZoG0lg+F4tpw3sJSbdguIEgVkAIrm0EmX/DcqnyWttIXDsdQWBCwXWtDuTa01ZaXN7srA2ql\n1Yq1SbidpFZYuZTl3FPlqd2e4t8zEVARSFtYtiANETdD5DVEv8Gqv8zcgXi3RpN41a5sAzEh\nSDtK2jviZTPybIepmJup4HZYijpGURQ7WXANC9UHQFoHlWRoiqtDtmDdh0ePSPlAIFwI0o4N\nxIQg7agMUrKgaQe3rPTWRE/HKIpiJwuuYaH6BEjyHQwx0dMxiqLYyYJrWKgIktQG4aEnip0s\nuIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0s\nuIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaFiiBJbRAeeqLYyYJrWKgOgfR4PNKpsGAi\nSAIPPVHsZIG1PlZHQHrM/0VTYcFLlSdPkMoeeqLYyYLsfqAOgPTIH8L0rMqTJ0hlDz1R7GQB\nNj9SGJD+91Jl09ef4VPUF+sekCjqy3XPqZ2Z84YDNggPPVHsZAE2P1IESWqD8NATxU4WYPMj\nRZCkNggPPVHsZAE2P1KCy9+P45e/zRylAzYIDz1R7GQBA4CS4A3ZxzrV/IasmaN0wAbhoSeK\nnSzI7gfqno8ImTlKB2wQHnqi2MmCa1io7rmtC0V9uQgSRQFEkCgKIIJEUQARJIoCiCBRFEAE\niaIAIkgUBRBBoiiAgCD9R1HXC9ewUBEkypZwDQsVQaJsCdewUBEkypZwDQsVQaJsCdewUBEk\nypZwDQsVQaJsCdewUBEkypZwDQsVQaJsCdewUBEkypZwDQsVQbpFv7NaV740inXhGhYqgnSL\nDrFBkGrCNSxUBOkWESSYcA0LFUG6Rb/h339neL//zad5/mzv18+5R4JUE65hoSJItygGaWHp\n18+4/3/jKrUvXMNCRZBu0Xqx4ff3N0bHw5Qvo3aFa1ioCNItWtn4TcegZQTy53kEqUG4hoWK\nIN2iGkjxCyOC9Fa4hoWKIN2iCKTiqd0vT+1ahWtYqAjSLdq92BCDFC2jdoVrWKgI0i1yFxt+\nPUW8/C0WrmGhIkgfEEk5IVzDQkWQPiCCdEK4hoWKIH1ABOmEcA0LFUGibAnXsFARJMqWcA0L\nFUGibAnXsFARJMqWcA0LFUGibAnXsFARJMqWcA0LFUGibAnXsFARJMqWcA0LFRCkYV/VYrsg\nLoqyKIpiJwuuYaEiSFIbhIeeKHay4BoWKoIktUF46IliJwuuYaEiSFIbhIeeKHay4BoWKoIk\ntUF46IliJwuuYaEiSFIbhIeeKHay4BoWqk+ANIp3sHTDxERPxyiKYicLrmGh+gBII0HyHnqi\n2MmCa1ioCJLUBuGhJ4qdLLiGhep+kEaCtHroiWInC65hoSJIUhuEh54odrLgGhYqgiS1QXjo\niWInC65hoQKC1Kh/IN3+MynqYnFEktogPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPRE\nsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPRE\nsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPRE\nsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPREsZMF17BQESSpDcJDTxQ7WXANC1UTSI/HY29J\nVGjcMwQpeOiJYicLqPHRagMpASZe8CBIZzz0RLGTBdf7ULWA5GFJh6HH/D9BOuOhJ4qdLND2\nx6kNpMdCzToOTRNBQnjoiWInCxwBjJpAWkjagrSe4f3vpcafyD81p75Qrad2j/nM7uH+WZY+\nCBJFLToCUrJgBukRX89rHKt5ahc89ESxkwVNAEgHQcpeI/Hy9zkPPVHsZAEDgNKx10jRZTu+\njwTw0BPFThY4AhgdeUM2Oo/jG7IEqWwDMflWkNrUuGcIUvDQE8VOFlzDQkWQpDYIDz1R7GTB\nNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUHwFJStK3dYyiKHay4BoWqg+A1BMk76En\nip0suIaF6n6Q+r4nSM5DTxQ7WXANCxVBktogPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktog\nPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktog\nPPREsZMF17BQAUFq1Auk238oRV0rjkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SH\nnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SH\nnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SH\nnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhagTpkc/L\n70ZBkIKHnih2skDaHq82kB4ZSOv9kR4E6YyHnih2sqA6H6y2+yNN6f2RTt3VnCAFDz1R7GQB\nA4DSAZAKt770IB25GXP+F7K8NzP1BWq8h2wZJNFdzTOQeJNz6hvUANLDg7Tcw9zfyXy9saxX\n41idndoduTfFt53DKIpiJ8tFIJxVE0iOoLBkefAFp7Y9MxKk4KEnip0seAYgOnD5u/AaSTAi\nFUBKSapw9W0doyiKnSzI7gfqGEjRZTvp+0jvQKqNUN/WMYqi2MkCbX+cjrwhG53Hyd+QJUiR\nh54odrLgeh+q2z8itAWp74dkniAdtYGYWMmCa1ioPg9ST5DO2kBMrGTBNSxUt4PUE6TgoSeK\nnSy4hoWKIEltEB56otjJgmtYqAiS1AbhoSeKnSy4hoWKIEltEB56otjJgmtYqLSBtHl/NrF5\ndwRapKhjFEWxkwXXsFCpACma34IUzX5bxyiKYicLrmGhUg/SSJDe2kBMrGTBNSxU2kEaCdJ7\nG4iJlSy4hoWKIEltEB56otjJgmtYqD4N0kiQzttATKxkwTUsVB8BKbpMR5AANhATK1lwDQvV\nB0AaCJL30BPFThZcw0IFBKlNff/6b9UCUjw/Zl+Nwq90oAyII5LUBuGhJ4qdLLiGhYogSW0Q\nHnqi2MmCa1ioDIAU5r+tYxRFsZMF17BQ6QOpJ0hHbSAmVrLgGhYqdSD9m9v9Sodv6xhFUexk\nwTUsVNpA6nOQeoL01gZiYiULrmGhUgFSNJ+DlGD2bR2jKIqdLLiGherjIPUE6bQNxMRKFlzD\nQqUcpJEgNdhATKxkwTUsVHeDNL4HqW8HSXKjJUUdoyiKnSy4hoXKNEgHvn+/lOWUCFLZBmJC\nkN7umRykfjufXcSrgiQgSVHHKIpiJwuuYaG6CaS14ReQQv+/IJKCVP2ilMqBOL5NyQbhoSeK\nnSy4hoXqHpBGglT00BPFThZcw0LVdg/Z/K7mR+9GQZDKHnqi2MkCany0Wu8h+9gsWSZuByli\nhyA5G4iJlSzvW/r2v7Gbf2rDOuvt+pJhKL6r+azKk8eBNBKkrQ3ExEqW9y2tFaSHv13ftLn1\npQfp3V3Nw5+5jvNfyIa/en39tWz8F7N9+gez2eyY/AFt/te0FPXi6CMktd6xrwjSeoYHA2l5\nWykux3PjO5DGnWnq70g1SA93Zjfz5M7wHlN6q/Mzp3ZxuU/L8dyLnL5yapd8uGjvtE/ROYyi\nKHayvGvon3ByNzO1zFyOV9tVu+hxnXosZIWllSePAymZTb9qMgFr//WToo5RFMVOlncNHUD6\ncVonj4BxVK0XG6bia6T2Ecl39RuQ+j79yFA6t/zxUrzy7h8BEqSjJlayvOnnn9JjPH+RWi42\nhH+jy3YH30dqBmnYgjQGl2SEOgaSn1fUMYqi2Mnypp8DMuu4pAYkfwYXnccdfkM2AqnrjoA0\n7II0HgJpjEASXTbPRJDKNhCTMyBFp3NuSbL0Kt31EaHQ/V0M0ngXSGME0v5odUAEqWwDMRGD\n9PMTmNmMSJfKMEj99vsd+gSkMmQFkNJ3eTfZi0+JIJVtICYnQAr/EqR8vgJSckmvfwPS6A5E\nDaRmyKZ8wfY579e8h5XmPWADMYGA5B9+3MPHXyO1qvLkUSAN+yCNySW9sfS1XnsgxdC5yeKb\nVb42z09DCbKEyGItnpym/Vq23d5qsFqWZW+1W2pikAIq2SujP/EaqRODVBitxtU3AunfVADL\nfUDCw/KqbWAZ1+O6V/OTvgl2avnkkE5Oldpbi/f2h2oT6EcjahiQ9L0h26bKk98FqZ8v4iXg\nxCNUn8HSi0Aak9HKgxQO4/JfBMu4ra1HexyTyaQ1NrWxWAsLpkqttt0VtWNZro2F/DOKmz4w\ndDNI4wakoQ7SGMPSvwXJzzuoxrEA0r+Zpc8dLOtpoO//MTm4nqEtIf6/47VC825qte3EtXF/\nwVSp1bYD1DaRMSBd/8Io/mE4q8tBch97CLDkIA0JSP1yYFaQeg+SL40rR72HJWLJTTggE8jW\n6ZS9Wm04XxuO18ZtbbPaFbWzT3Xth5MgXf7BoOiH4azaQepykPpofjmhC7MrHKMQpHGMRqvX\nb94VJHck+7BgASc6onEpHOocskItmxwLfQivbXLuWlxWm4TbZTUISP5q3S26G6S+AFJ3BKQh\nBykadPLackhex8fNv9abHDf+3C8Faa82JLWlY1ZWlx8easVJ5xMtyGv5Akht3K1d/qNP1EAg\n3agPgPR8vgFpOALSOPqXN9vaXI5B8qzMtX6tBViWyT4c0d1av9bCpJ9PFyBrhdWkNdWRCdJ7\nkJ7/QIrByUHqhoSH8Ipq4WgtLlCtv+Y9ZL7oXlEtpGxB8oDEtSGqueMaTY59ccGm1qcLlNZG\nfbU4MkHaAenZDFLXRVcf+vjSRN9Pyz+eldG/FIpHqzG8XnLFPtQcLP6yhQdnSIDooyPqR7Zy\nLTv8tmuV1aQ1aaxv/Rah83o+3R+Fe5B84d9Om8lZ5zsHki/Pc+u6yz9jmJ33/bj8O/8Rez+/\n3HXTc3GMamM/RbW1E+L5yS3Ia6vS2nhNbWrYbnxbK1jUatc9nZZa/HTu6EmwMpDCB5OOq/Jb\n5AWSG5G6wogUzXfxiDQ6kHZGpGG9OudO5fyp3ricAvrjM2tw04OvDX7ACSNUn9b6ttqQTWqq\nDUprb56O+RFp+UiFkCQZSAs5YzqbgNSt5HiQxqUWQHJELMUxhizAEo7gFqS3tXxy2yW11a6o\nRZPTzT+6VsNk+QKQTnzsvB2kZ/Qea9cNFZC6Eki9u7owrFfuloPgKAvTQwWWQXaI08khnfxQ\nbSrVPhRrA5LIkyDhQVo/+dC/9m/3+ohr79+5mU/7+gWWFaR+gSAGacgP1VAAae/gbmr7K3+o\ntgOS1LP23N9aTJinKgXp/3YlaOhDuhmkPgVpHnJKIM17clxB6pfTtWnFbH5ZtLy68iA5gz4D\nqY8gW49atGDIj2i5tpmcKrWd1jhb24s1VWq17VprRyLvgCR4qqZBmn5+5JcbGkAaM5DmAen1\nmikFyY1BMUhjAtIYgTQEkOZV/SDTFdkowZI2U2NtEm5X41RYmy7wlD6dCfaj91TtQTUgTT/i\nAakOUu9BeuYgPTOQhgBSv4K0vPSZHGbukvZretn5855fB7N5fqmF325+nCt3wm5r7NYm4Xbv\na4e7eQvSFbHaatOh7aqHZEfVHtQD0glVnnwC0rABKcx3r/kEpBdmHqTOg7Ry1HULSIPjqPOH\no8tAimvDUooh66LDmNbSQxxqOyANlQ4C1ob4cbr1R9drsixD+ih/Q2IbfEoAACAASURBVPYv\ngOTQmEHqAkjzCOWvh7sBannR9MJlmX0dhtHDsZ6/LYPXepASkDxkvv9DzReH2XbwxTWP+xnr\nb8W1Nk8Pa20q17rC5G7ttc8qq91aw2YZWmq7FrZB+kkksKo8+XaQllkP0jJezSDNo8yyn1c4\nZuYcLA6OzrOxgrUOZK7rp3W2S2vL/9G60UyXFONSvuBQLcqSl+IljbXBVq0rLVmmCNJRkF6g\nRFcbAkjL6yA3u5y+vXbzsFQXWJZVAyzDPHp5kFaS3BEa3GGaSkcyHMbdGnWjLIN0WlWQunF5\nL6ibwXlGI9AQQPLkPN0FhWFZezm3+wfE9IzGr24FycOyzDsSfM0rnl9GsuWI5bV1aFsPaVeq\nDcUFn6ltxrXNdvfVTmUJkwRpD6Sn+4MGD1IYgYYwP1/TGxZW5k8ruLUXkjpHmQdpmfFHx9fi\n2cHNO+ieUa1zteG51oa49uwcb0stmXTT+WSt1rjaFbUP/mh57WKQfuKvG8IgUHgf6ZJTOw/S\ncyEnGoGW07XlxG+enefdJWwPUueuQ7i1Xc8PAR4PzjPBaj0+fvr18IyKxcnrah+xv6J29d66\nFqSfqO9/QIPJ5kOrG5C2tzQPCx+NN2N+zud2bkAa1hdJjqvl3O5VX8qOpGUcGVY6nqG4zrpD\n0Pm9Hx+K9IgltfxwSmvl1WqSWtRq0+d+9FW1S0H6iRr/J18YvlBynm/+eMLmI0LOICCT3HEi\nXdh8x77X/plf9TxXcPoUHHddLmLFcxTBMq1myx6PZ4dNLcAzpJC52rBMJ5Ob2jOv+QVTpba/\nXeNqx2pTzR5cG95sN51/OufeRxrL2gfpJyD1M4MzBQh+3HjVQNLPZu4nWRzuu/zY3rFvnW4B\n6ZmCMoPzTOf9oLO8TFr5WAkogjTE0/EhDvPZgYhiVTK/Eb9Ev2wDMbkRpJWkn/U07ycMT20Q\nFUDafGg1uu9yfDNzCUhR92ezu/OpgT9KWziOSFHHKIpiJ8sbkA6c2v2sI1AE0k98brd+C/9x\nkH6y2222gfTurubhPN7PP8vz0RbRDEW1SfQayYH0Ey48rKd2bu7wqd36Jf7rghWkcDNzAUgU\ndYtOgBSuDmxP7Y6DtFE8Ivm7YD4On9qZOW84YIPw0BPFTpZqux6+/P0TXW0IVxh+pvRCAxgk\n+WskM0fpgA3CQ08UO1kAIJXfkF3v/zL51zjTT7ZyRW/fR4qudAdwjl7+NnOUDtggPPREsZMF\nAVIDCkfV/Iasv/z9KC2cJd0zB3YwxERPFkVR7GSptrMWkE5ZSvfMgR0MMdGTRVEUO1mqPXgC\npJOfuStvjf5TczNH6YANwkNPFDtZqj14ZkQ6px2QgB8Kp6jbRJAoCiAtIJ36C1mK+rTUgXT1\nj6WoK6QFJIqiRCJIFAVQAOn0twhR1Kel4dSOIFHmpQGkZfYnPFKUNWkBKflE+VH9R1HXq9qD\nBImi2lTtQYJEUW2q9qAekA78KVOuT+9h6k+o2oNaQJrOXLT79B6m/oSqPagGpOyrTw7p03uY\n+hOq9qAikOT69B6m/oSqPdgIUjxYHPjKrZoIEmVL1R5sBSn5GqF0gVDxJxtO/hnFp/cw9SdU\n7cE2kKLu9l+Wnz8c/i59gkTZUrUH/+//+rIykLJTu+grI6Mvtjv0Xfo8tbtBv2Hqt1KkWlTt\nwUaQYpIKIO194WpNBOkG/Ram1nmCdEzVHmw/tQunb1MJpMPfpV98Q1amT+9hvaqCRB1TtQcP\ngZTN+Mt3su/SL/2peRM2W316D+vVPO68zup+14cwufzzu6z231qidlTtwfMgSb9Lv/iGrAyl\nT+9hvfpdXhyt//33W/ovqnCk2le1BxtByl8jba6HR5ccpCAtXrxqh1SOz3+/yfzvZh1qX9Ue\nbAOp+oas8Lv0CdINSkH6dSd0ayUHiWd2NVV7sBGkmoSvbHhqd4M2I9J/NZCIUlXVHtQCEi82\nXKLjp3YEaVfVHjwPkrT9efn7Bm0uMfymFx+iixEE6Z2qPQgYkYTiG7I3KBttfuPH6PK3m+OZ\nXVXVHiRIFNWmag9qAIkfWqUsqNqDBImi2lTtQQ0gndan9zD1J4RrWKgIEmVLuIaFKj61O/nd\n35/ew9SfULUHNZzaESTKgqo9qAGkZfYnPB7Vp/cw9SdU7UEtIPEriyntqvYgQaKoNlV78CtA\nGvZVLbYL4qIoi6IodrJUe1APSC1/w7Qj6Z45sIMhJnqyKIpiJ0u1B7WAdOpL9KV75sAOhpjo\nyaIoip0s1R5UA5L/AiKJpHvmwA6GmOjJoiiKnSzVHlQEklzSPXNgB0NM9GRRFMVOlmoPtoK0\n+YKtUy9qgiNE0j1zYAdDTPRkURTFTpZqDzaCtPnm/J/164JAIJ35ZEPznhnFO1i6YWKip2MU\nRbGTpdqDbSDl35yff49d9B1dzd+lX/rOhotBGgmS99ATxU6Wag/+3/91ZRVO7dKvV01BOv5d\n+uX3kURq3TMEafXQE8VOlmoPCkFyvZ+BtPeFq0UBQWrVON7wQ6g/qSMXG/LTt6kAUvt36eev\nkU48idZfMRyRVg89UexkqfagBKTQ+OnJ3k80PcVzDad2p0hq3TMEafXQE8VOlmoPSkH6CY/x\np+TKp3ZNr5Guv9gwEqTVQ08UO1mqPXgEpOhid3wVfEVhSi45EKRaljM2CA89UexkqfbgIZDC\nWdpP9GJoRSGs0fJd+ve/IUuQgoeeKHayVHuwFaT3OgoGQZLaIDz0RLGTpdqDBKltB0s33Mty\nxgbhoSeKnSzVHoSBdPjFTfyG1MkvP2ncMwQpeOiJYidLtQdxI9JRESSpDcJDTxQ7Wao9qAGk\n02rcMwQpeOiJYidLtQfVg/TI5x+PTaFxzxCk4KEnip0s1T7VDtIjA+nhCXoQpDMeeqLYyYLq\nfLCaQHos2KTD0CMUFjXuGYIUPPREsZMFDABKB0Dy/xdA+t9LjT9x5Ke/qauk+9TuMe2AtJ7h\nESRKhVSD9PAgPeZzu4c7w3usBafGsZqndsFDTxQ7Waq9qhwkR1BYsjz4glPjniFIwUNPFDtZ\nqr2qGqSXdl4jcUQ656Enip0s1T61AlJ02Y7vIwE89ESxk6XapyZAmqLzOL4hS5DKNhCT7wWp\nRY17hiAFDz1R7GSp9iBBatvB0g33spyxQXjoiWInS7UHCVLbDpZuuJfljA3CQ08UO1mqPdgI\n0ulv+t6KIEltEB56otjJUu3BNpDCt6YSJLkUdYyiKHayVHuwGaT0C4ojFqJvjHz/hd/5xhg1\n7hmCFDz0RLGTpdqDTSBtvuk7/jauQ1/4HZtKkCmrcc8QpOChJ4qdLNUe/L//e5aVX2z42QNp\n71sh3+kjIElJ+raOURTFTpZqD7aClH7TdwZS+xd+x46HedlV454hSMFDTxQ7Wao92HrVLnqI\nvqTk8Bd+55YYNe4ZghQ89ESxk6Xag+0XG0oj0uEv/N543qmRf5BEXaW2iw3h359kweEv/C6Y\nAtT4K4YjUvDQE8VOlmoPNl61y77pOwXpwBd+x6aHedlV454Ze4LkPfREsZOl2oONr5H2dIKG\n+0Hq+54gOQ89UexkqfYgQWrbwcLtdrOcsUF46IliJ0u1B8+BdOajdwRJaoPw0BPFTpZqD54c\nkU6IIEltEB56otjJUu1BgtS2g4Xb7WY5Y4Pw0BPFTpZqDxKkth0s3G43yxkbhIeeKHayVHuQ\nILXtYOF2u1nO2CA89ESxk6XagwSpbQcLt9vNcsYG4aEnip0suIaFqu27v/O7UZz5FiGCFDz0\nRLGTBdT4aLV+9/djs2SZIEhnPPREsZMF2PxINYLkv/w7LBDf1oUgBQ89UexkQRMAUutrpEfx\nK4slt3V5gXQsI0Wp14E79oFu60KQqC9U61cWP2C3deGpXfDQE8VOlisoAOjgqd06I7+tC0EK\nHnqi2MkCbX+cGm80tvMaiSPSOQ89UexkQRMA0oHL36DbuhCk4KEnip0scAQwOvKGLOa2LgQp\neOiJYicLsvuB4keEpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElq\ng/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElq\ng/DQE8VOFlzDQvV5kA58gfG3dYyiKHay4BoWqo+DdOQr9b+tYxRFsZMF17BQ3X9bl+zvkXiT\nF+obxBFJaoPw0BPFThZcw0KlAaSUpApX39YxiqLYyYJrWKjUgVQbob6tYxRFsZMF17BQESSp\nDcJDTxQ7WXANC5UCkPp+SOYJ0lEbiImVLLiGherzIPUE6awNxMRKFlzDQnU7SCNBCh56otjJ\ngmtYqAiS1AbhoSeKnSy4hoWKIEltEB56otjJgmtYqAiS1AbhoSeKnSy4hoXq9tu6lEBq/aTD\nt3WMoih2soAaH63bb+vSE6TgoSeKnSzA5kfq9tu6vECKz+UykLafGIpmv61jFEWxkwVNAEjn\nvrJYcFuXvn/9FzSmnwYfs0+D88PhlAkduRvFhLitC0GivlFH7o8Eua3Lv/O6I6d2I0/t3ttA\nTKxkuYICgNqu2kWP65Twti4EKfLQE8VOFmj749R6sWFC3dYlA2kkSOdtICZWsiC7H6iWiw3h\nX8BtXQhS5KEnip0s0PbHqemqnTuDg9zWhSBFHnqi2MkCBgCl2z8idBykMP9tHaMoip0suIaF\nSiNIu39A+20doyiKnSy4hoVKH0g9QTpqAzGxkgXXsFCpACme7wnSYRuIiZUsuIaFiiBJbRAe\neqLYyYJrWKi0gdTnIMUnft/WMYqi2MmCa1ioPg5SXwepJ0jvbSAmVrLgGhaqz4AU2CiA1BOk\ngzYQEytZcA0LFUGS2iA89ESxkwXXsFApB2msgyS5Y5mijlEUxU4WXMNCdfdtXcb575HGbD7U\n++zPk7KbwGRm/GMlSonuHpHG7Yg0iEekI7eEKWQ5JY5IZRuIyd8ekWpPfm34HKS+DlJfA2nz\n/Q5NUtQxiqLYyYJrWKjuAWmsgDQQpLM2EBMrWXANC5U5kCJ2CJKzgZhYyYJrWKisgTQSpK0N\nxMRKFlzDQmULpPSPLAiSs4GYWMmCa1iodICUXHyIPyOUgjQSpKINxMRKFlzDQqURpJyryOYN\nSOnnxvcORCVouwhS2QZiQpB2FDr+DUhjAaQ+thlzrso/pvIek6KOURTFThZcw0KlEqQ+Ksdc\nvTjqk9Ft92+X9k/7FHWMoih2suAaFqrGb1rdLDj2LUInQEqqY8LV5m+XqiD5eUUdoyiKnSyQ\ntser8f5I+YKDt3VJQer2QeprII1nQBpjkCQfds1EkMo2EJMvBckPP/LbuiQgdVWQhmaQRjFI\nu6PVARGksg3E5EtBSsjJZtvuRjE+/ce0x38cddEHul+f/I4+4D1/FjyfX7d1IK1WL5DinxJt\nOGa3tYhn8ztepJ8i5wfKKYFOgdR8W5fnHkjjCZD69E8uxuRPLkogRT77IG3+NKMGGf+Mg3I6\nApL8ti7/QHJTDqRhnS2e2qXz4dQuKeYXIpJTvX8Tm9O+0Z0alE77xu3kuqA4PW1q2XOOF4zr\nQ1KL/rqkvF2uco2ndgp0cERKZttv61IFqUvB2YIlAil5/TRGYE3ZB43S96PcZPHFVlIb39cG\nv2CeHvPJf/9MlVpwK1lsapvJo7Upq5VWu6JWivX9IAlv61IDqYOA5E/7xtU3fsMpfgPKgxQ6\nNh69xqybPWRjVvOT43ayVhvTx6lS29tufzVpbZksZ2mOJa0VV/t6kKS3dZGC1Kew9P07kPy8\nA8nzEI9WL5AWPjwcOUhRzUOWHvhxM3moFhZMlVptu2O1tpzbLIefHyzyd4N04rYuGUhdBM4M\nUh/mM5DGFKsYli1IQ1SbyRlXWMJ20zI6RTU/Wq2AhSM6epb80S7XNt1Vq8WrXVErrDZUarXt\nrqzVns7XgtSkFpDGDUhDCaR+nb0CJF/yx87Vlvl+SGqjG9kiRV05fLQ27NVK6FxfQz5VgiQB\nqctAGsoguW/uagRpWEEaR3/aN89MfQLS4F9cLVom+6FUW0YyR9Y0riPbsP4cfz6Z18b0ca+W\nrXb17/0rapNwu+zpEKR3IL1O5Z5dDE4dpOEtSNGgswGpH11ju/kIqj40fQTLMlmoLf3va0vH\n9IG6pOYmPXVjNDlmq42V1d7UhkqtZFGrHfzR19cIUgNIzzpI3ZDw8AakfnQg+doYgRSACCOU\nA6mPawGWCKQNEGvNr1GalNcutpf96LtrfpIgvQfpWQep66LLeO4llJ9xw9W4znoeIsjGcJrn\nilHNwbLWZkJ8zR9GX/SHPD7ErtaHWtYfewvU12rP4e5Y3/uGbJPqII0VkLp9kLoYpCkFaQx/\nmT4fhGWECqd5rhjVUpBWOAZ3BJNaH2r5gs1kodY318ZCbajUCp3XWivEktYOPZ3G2vpIkKQg\nrRfmugMgDevVuRykfgcWB8qQdcIKy4la3l1DvgamNua1/ljtoliwp0OQjoEUrmh3FZBGN0B5\ncjxI41JLQFrO/Prl/C3UwlELR3AYssMYrqzv1YZyrd/WemAtb8WdyemCH31t5Hc1gvQepO4f\nSM89kPomkPoNSH6kWUCK3o6Sw1JY0FgbrquVf/R0xOLiGiYLQZKANBRA6v3aC0ijB6nzIM3s\nzKd9/TLSeHjiaTc/DPmhGlawBLAMu5P316LmLdfujlXMctzzT4NU0QukeWIFyRVe5EzhD5QC\nSH7tudq/tn4t9OPV8m7oPPNaPq43hum9Jv83g737Z5me9mtTXJv68BjV0qWbBcprteeArAEi\nX9CDF+vmEanPR6SuNCK5c7txdOOVfxm0gjTM8+56nxt1XupWkJJRx49Q6y++aEH8S/BIbTqy\nXX9pbZI+nQtiTRjPPz0iVUFaUBkzkOYzu+jigzvTewvSuII0xCC5s8IZlrmUnCbsdtdQWPC2\nNpVXLjfvtbUtSHdn2YX6hP2ecA0L1e0gPXOQXvNdPD+sIC3leZiZXxNNnX/NNI803UJS70YZ\nR9lyOBxkAaSoNgwbyOLDGNXqIO3WpHAKa1P7dqXVoLXp/XbvnyrfkN3TCx2HxvM5ZOBEbywt\ns+5C3ejGLzck9fNbTjMQEUcOpMFxlIDUrUcmhWyoQeZGsiGqpQvc5LQUi7VWmfnz7gM2EBOC\ntKMKSM8UpGV2Aal35e7V4aMfZWauFnJcJy+DjNPi5CDzFwA7D918IHytWwanqObWdPR0oebm\nk9qwTnZLbdiv7SyYKrV8u5pno0W1Ngm3u6JGkCQgRS+SAkj96M/sHEjzKDNb+c86+NHLDzIO\nDseDh8fXVlqmZDaa6cJG6eTOAkBtqtQMSrobNiJIuyB1M0genGc8Am1Benbuoz2eq8Gdyr2c\nnissS20edZZjtazrIXMHcEj/LR5Uaa28GnVaBGkXJPfZUg9SAGcI8/0K0vI6aHDjlQdp8pit\nA9KCizu7cdDNR2IZycIp3vD0ZHVuJHOD2ODn15pfuqzhRr1ibYhq4eTO17r92gCtTRd4grMU\nLSo1grQH0tP9QcNzIScHqXMgPf388jrIYbeQ5MhZRp3Oz/mR4LkW47lnNjPkxdejm/YPz+59\nbeE0X/AMC2q1jdGZWndd7arI72s8tauCNC4vkQJIywC1DFiD+/zQ4E/fFo6efihxI9AQhhY3\n8/RNv647r+pK/viUVatRRQ3vV0H8lC8E6bG5F0VY+Gi8G8Vr3/TucrYDZwwj0OBGIF92Q9LC\ngz9jm1nxxS6iai77vR8fimg6nx+eviOG0rrPik+83WbBsFlwcW2enCq12nZva4JYE+apfiFI\n6dcSpwsPfGXxfLkhgPNc3hpKwOl9OUPnuWpazZadP0TzQ7E2ZNPDsFsb8ppfWq5Nm1qyAFAb\nWmoviwn+o+U1WZbsqX7j+0jhO1Yf268sXqdnVZ780twpKO6i3Drvz9bC+muj++kEpNS98sM3\nB6J93ZoNwkNPFDtZ0ASA1ArSY5qyLy6OQXp/N4p1UGman9K5ZcGBZ0VRN+sekF5KyMi42cxT\nlC01gxTu5iI4tTNz3nDABuGhJ4qdLJdgcF5HRiSHU7g/0rrCLOmeObCDISZ6siiKYifLJRic\nF+bUbpZ0zxzYwRATPVkURbGT5QoKAGq9/L3e0yUjq+3yt5mjdMAG4aEnip0s13BwWs1vyK63\nNi8tnCXdMwd2MMRETxZFUexkgSOA0T0fETJzlA7YIDz0RLGTBdewUN3zLUIU9eUiSBQFEEGi\nKIAIEkUBRJAoCiCCRFEAESSKAoggURRABImiAAKC9B9FXS9cw0JFkChbwjUsVASJsiVcw0JF\nkChbwjUsVASJsiVcw0JFkChbwjUsVASJsiVcw0JFkChbwjUsVASJsiVcw0JFkChbwjUsVASJ\nsiVcw0JFkC7W7yzJhvAo3yFcw0JFkC7WzIOEJIJUFq5hoSJIF+s3ehRsSOXCNSxUBOliBZCW\nM7z18dePVuuyaPL3lyDtCNewUBGki7We2i2sJP8t5ERT8UoEqSxcw0JFkC7WerFhHXTc4xap\nZCWOSDvCNSxUBOlihRHpN5ywZSCFE71oJYJUFq5hoSJIF2sl5zedjc7rfsOaaZkqCNewUBGk\ni5WAlJ/a/Red9PHUrk24hoWKIF2sGKTXABSd0P0XQIovNPBiQ1W4hoWKIF2s2uVvj0u4GsHL\n32+Fa1ioCBJlS7iGhYogUbaEa1ioCBJlS7iGhYogUbaEa1ioCBJlS7iGhYogUbaEa1ioCBJl\nS7iGhYogUbaEa1ioCBJlS7iGhYogUbaEa1iogCAN+6oW2wVxUZRFURQ7WXANCxVBktogPPRE\nsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPREsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPRE\nsZMF17BQfRFIz+f8OKyPz/AYTU5uwTPUdhZUa5Nwu3hy2s2STr6pTXtPdbMgqaWrTXu1obLd\ndrdPlVrx6JRWI0j7T75SHMq7vHAYy837XBY8n9ljPvnnVHvutV1UrtV2bVTbWS2vVVcjSLtP\nvvBr+hlPRnuxvJ+pPyWCtDfeXLfL88md34A5pfFgVllQr+2v3Ox5erXN75/Nr6JqrbhadQA5\nNLiUa+9XI0hSkHZOMLadsNPNy5neUF6QHohK0HZ988WG7f6MJuu1aXe10gu8ndV4arend7+m\n092aLYj2Ya0TWqWoexVFsZMF17BQ3QcSYgdDTPR0jKIodrLgGhYqgiS1QXjoiWInC65hoQKC\nVNHzecuPoahPiSOS1AbhoSeKnSy4hoWKIEltEB56otjJgmtYqAiS1AbhoSeKnSy4hoWKIElt\nEB56otjJgmtYqAiS1AbhoSeKnSy4hoWKIEltEB56otjJgmtYqAiS1AbhoSeKnSy4hoWKIElt\nEB56otjJgmtYqAiS1AbhoSeKnSy4hoWKIEltEB56otjJgmtYqAiS1AbhoSeKnSy4hoWKIElt\nEB56otjJgmtYqAiS1AbhoSeKnSy4hoVKAtLj8dhMTQRJ4KEnip0siK6/QAKQHvN/+SRBEnjo\niWInC6bv4ToO0mN9iCB6qfLkCVLZQ08UO1kATX+FToLEU7tTHnqi2MmC6Xu4zoHkSfrfS5Vt\n+Bey1Jfr9KkdQaIoFEizKsMxT+3KHnqi2MmC6Xu4CJLUBuGhJ4qdLJi+h0t++fsRvUaaVXny\nBKnsoSeKnSy43odK/obsY+Ibsmc99ESxkwXU+GjxI0JSG4SHnih2suAaFiqCJLVBeOiJYicL\nrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicL\nrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicL\nrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicL\nrmGh4l3NKQogjkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVB\neOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVB\neOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVB\neOiJYicLrmGhIkhSG4SHnih2suAaFiqCJLVBeOiJYicLrmGhkt+Nwk2HxZUnT5DKHnqi2Mly\nvucvkfz+SMs0QTrjoSeKnSyIrr9Ap+7Y97qteShUnjxBKnvoiWInC6Tt8cKAxJsxU39cZ+8h\nS5AoajoH0uPBiw3nPPREsZMF1PhonQPp8eDNmM946IliJwuq88E69xqJI9I5Dz1R7GQBNP0V\nkl/+Xl8meVWePEEqe+iJYicLpu/hkr8hS5DOe+iJYicLpO3x4keEpDYIDz1R7GTBNSxUBElq\ng/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElq\ng/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElq\ng/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUvKs5\nRQHEEUlqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTB\nNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTB\nNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTB\nNSxUBElqg/DQE8VOFlzDQkWQpDYIDz1R7GTBNSxU5+5q/uDNYgbybgAABfJJREFUmE956Ili\nJwui6y/QqbuaP5LbmleePEEqe+iJYicLqPHROnsP2YkgnfDQE8VOFkzfw3X2HrIE6YyHnih2\nsoAaH63Tt76cp3hXc+qP6yRI7iUSQaL+uM6B9IgLleGYp3ZlDz1R7GQBNP0VOn2xIajy5AlS\n2UNPFDtZAE1/hc7c1TzliCAd99ATxU4WSNvjdeau5o9Z6+LKkydIZQ89UexkgbU+VvyIkNQG\n4aEnip0suIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaFiiBJbRAeeqLYyYJrWKgIktQG\n4aEnip0suIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaFiiBJbRAeeqLYyYJrWKjuAqlT\npGF96IZhfsgeh/A4hMfseIpaJO8JgAlBUqG/CBKlXQRpF6R1sutKv9/bJP511yVDjR+Thk/2\nClURQdoRXyOVPfREsZMF17BQESSpDcJDTxQ7WXANCxVBktogPPREsZMF17BQAUGqiH9qTn25\nOCJJbRAeeqLYyYJrWKgIktQG4aEnip0suIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaF\niiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaFiiBJbRAeeqLYyYJrWKgIktQG4aEnip0suIaF\n6i6Q+pfcXDztZ+P5ZVFUX/fhseOxdyAQLgRpxwZiQpB25EB6r+FIbcgmhz579Mj6h3otLI4W\nXFz74I++vHbCgiDtgjT4keelcewT5fP/luQLqL8kgrQL0ujVv7gZwuSsYa31Y6xofq7F9WiB\nn+zH+DFMLnpb68ds8u7aB3/01U/n0FMdCdJ7kGZuFmWzyYJxzLaIpkfqy/VHQBLe1Xz0hFRW\nrOvbXlUrimInC6LrL9Dpu5qH5ZUnz8vfZQ89UexkwfQ9XKduxrw+zKo8eYJU9tATxU4WTN/D\nRZCkNggPPVHsZMH0PVwYkHgzZuqP6x6QKOrLdc+pnZnzhgM2CA89UexkwfQ9XARJaoPw0BPF\nThZM38N15q7m7Ze/zRylAzYIDz1R7GSBtT5WZ+5q3v6GrJmjdMAG4aEnip0soMZH656PCJk5\nSgdsEB56otjJgmtYqAiS1AbhoSeKnSy4hoXqnm9a1XRtXFEWRVGY5aQI0gelKAqznBRB+qAU\nRWGWkyJIH5SiKMxyUveARFFfLoJEUQARJIoCiCBRFEAEiaIAugWk5CN5n5T7lODj84F8hs8n\nWaI8lOwWnyj6x4buACn9kPgHNXeLjjDFr4/5oHQEefi/zAl/XmBFN4CU/dnSB6Xo8DzCo4o4\nKnZL+LOCfw9q9k2b/hJI/jApOIVxGfTsmknJbpkI0q5U7RHXvh9vmSWDll2znlF9OgtBqkjT\nHgnfEfvRGIv0NMtjM/EhEaSKFO2RR2Hqc9LTLATptP4USHpe4etqltJX2XxGBKkmFVeEphBC\nyYsB/3Ltw0mm6BeMht3Cy9/7UnA56KX4nUcNWdy/n44SvQ368Sx8Q5ai/rQIEkUBRJAoCiCC\nRFEAESSKAoggURRABImiACJIFAUQQaIogAgSRQFEkID6WfWablj/oL0sFXWHeHCAOghSC2tn\n1qduFI8NWO3dfpwLkqRXPDRgEaS/KR4asHyz+9O7ed4tXE76khXTx2g9t5n7dy3c8hyo4+KR\nASsDaWZg5SImaQPST+Bn3mQFKWxHkNSKRwasFKTwkJIzbUHKptZN4zGNIKkVjwxY+aldjEVc\nD6d7RVR+ClsQJMXikQGrBlJybrcFKT63S0DiayT94pEB692IlK+YX5CYNieD4eIDQdIrHhmw\njoL0OgIZJKUtCJJy8ciAVT21i+rJcOP+zTdLLunF21DqxCMD1j5IU/E1UnYC54njayRj4pEB\nqwLSFGFUAGmKePmJRqRoM4KkVjwyH9NxOsiRXvHQfEwE6ZvEQ/M5+RdOh9andIrH5oM6CtKF\nUaiT4sGhKIAIEkUBRJAoCiCCRFEAESSKAoggURRABImiACJIFAUQQaIogAgSRQFEkCgKoP8H\nRZj+IWtFl2wAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#rescaling time to years\n",
    "time <- seq(from=0, to=100*(ND+1))/(ND+1)\n",
    "#changing time to the rescaled time\n",
    "RES2[ ,\"time\"] <- time\n",
    "\n",
    "#labeling of the output from ODE solver\n",
    "label <- c(\"S1\", \"S2\", \"S3\", \"S4\",\"E1\", \"E2\", \"E3\", \"E4\", \"I1\", \"I2\", \"I3\", \"I4\", \"R1\", \"R2\", \"R3\", \"R4\")\n",
    "label1 <- substr(label, 1, 1)\n",
    "Age <- substr(label, 2, 2)\n",
    "\n",
    "df <- data.frame(time = RES2[, 1],\n",
    "                 label1 = rep(label1, each =  nrow(RES2)),\n",
    "                 Age = rep(Age, each =  nrow(RES2)),\n",
    "                 value = c(RES2[, -1]))\n",
    "\n",
    "#plotting  the data\n",
    "df$label1 <- factor(df$label1, levels = c(\"S\",\"E\",\"I\",\"R\"))\n",
    "df$Age <- factor(df$Age)\n",
    "df %>% mutate(label1 = recode(label1, S = \"Susceptible\")) %>% \n",
    "       mutate(label1 = recode(label1, E = \"Exposed\")) %>% \n",
    "       mutate(label1 = recode(label1, I = \"Infectious\")) %>%\n",
    "       mutate(label1 = recode(label1, R = \"Recovered\"))  %>%\n",
    "       mutate(Age = recode(Age, \"1\" = \"0-6 years \")) %>%\n",
    "       mutate(Age = recode(Age, \"2\" = \"6-10 years \")) %>%\n",
    "       mutate(Age = recode(Age, \"3\" = \"10-20 years \")) %>%       \n",
    "       mutate(Age = recode(Age, \"4\" = \"20+ years \")) %>%  \n",
    "ggplot() +\n",
    "   geom_line(aes(x = time, y = value, color = Age)) +\n",
    "   facet_wrap( ~label1, ncol=1, scales =  \"free_y\")+ \n",
    "   xlab(\"Time (years)\") + ylab(\" Individuals\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}