{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic SIR model (discrete state, continuous time) in R/Rcpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(reshape2)\n",
    "library(Rcpp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cppFunction('\n",
    "  List sirc(double beta, double gamma, double N, double S0, double I0, double R0, double tf){\n",
    "    double t = 0;\n",
    "    double S = S0;\n",
    "    double I = I0;\n",
    "    double R = R0;\n",
    "    std::vector<double> ta;\n",
    "    std::vector<double> Sa;\n",
    "    std::vector<double> Ia;\n",
    "    std::vector<double> Ra;\n",
    "    do{\n",
    "      ta.push_back(t);\n",
    "      Sa.push_back(S);\n",
    "      Ia.push_back(I);\n",
    "      Ra.push_back(R);\n",
    "      double pf1 = beta*S*I;\n",
    "      double pf2 = gamma*I;\n",
    "      double pf = pf1+pf2;\n",
    "      double dt = rexp(1,pf)[0];\n",
    "      t += dt;\n",
    "      double r = runif(1)[0];\n",
    "      if(r<pf1/pf){\n",
    "        S--;\n",
    "        I++;\n",
    "      }else{\n",
    "        I--;\n",
    "        R++;\n",
    "      }\n",
    "      if(I==0){break;}\n",
    "    } while (t<=tf && (I>0));\n",
    "  return List::create(_[\"time\"] = ta, _[\"S\"] = Sa, _[\"I\"] = Ia, _[\"R\"]=Ra);\n",
    "  }'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sir_out_list <- sirc(0.1/1000,0.05,1000,999,1,0,200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sir_out <- as.data.frame(sir_out_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "if(dim(sir_out)[1]==1){\n",
    "    sir_out_list <- sirc(0.1/1000,0.05,1000,999,1,0,200)\n",
    "    sir_out <- as.data.frame(sir_out_list)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th scope=col>time</th><th scope=col>S</th><th scope=col>I</th><th scope=col>R</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><td> 0.000000</td><td>999      </td><td>1        </td><td>0        </td></tr>\n",
       "\t<tr><td> 1.891201</td><td>998      </td><td>2        </td><td>0        </td></tr>\n",
       "\t<tr><td> 3.470562</td><td>997      </td><td>3        </td><td>0        </td></tr>\n",
       "\t<tr><td> 4.169704</td><td>997      </td><td>2        </td><td>1        </td></tr>\n",
       "\t<tr><td> 8.149657</td><td>996      </td><td>3        </td><td>1        </td></tr>\n",
       "\t<tr><td>11.145904</td><td>995      </td><td>4        </td><td>1        </td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|llll}\n",
       " time & S & I & R\\\\\n",
       "\\hline\n",
       "\t  0.000000 & 999       & 1         & 0        \\\\\n",
       "\t  1.891201 & 998       & 2         & 0        \\\\\n",
       "\t  3.470562 & 997       & 3         & 0        \\\\\n",
       "\t  4.169704 & 997       & 2         & 1        \\\\\n",
       "\t  8.149657 & 996       & 3         & 1        \\\\\n",
       "\t 11.145904 & 995       & 4         & 1        \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "time | S | I | R | \n",
       "|---|---|---|---|---|---|\n",
       "|  0.000000 | 999       | 1         | 0         | \n",
       "|  1.891201 | 998       | 2         | 0         | \n",
       "|  3.470562 | 997       | 3         | 0         | \n",
       "|  4.169704 | 997       | 2         | 1         | \n",
       "|  8.149657 | 996       | 3         | 1         | \n",
       "| 11.145904 | 995       | 4         | 1         | \n",
       "\n",
       "\n"
      ],
      "text/plain": [
       "  time      S   I R\n",
       "1  0.000000 999 1 0\n",
       "2  1.891201 998 2 0\n",
       "3  3.470562 997 3 0\n",
       "4  4.169704 997 2 1\n",
       "5  8.149657 996 3 1\n",
       "6 11.145904 995 4 1"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "head(sir_out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "sir_out_long <- melt(sir_out,\"time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAG1BMVEUAAAAAujgzMzNNTU1h\nnP/r6+vy8vL4dm3////njUASAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2diXbb\nSrIE7aH4rP//4nclbli6garqBLqajDxzLEsW4yYKFQNusv98E0Ka86d3AULeIYhEiCCIRIgg\niESIIIhEiCCIRIggiESIIIhEiCASka6b2ftzV3SspLXo5UVtsxT7bQkiaWBCFL18KEQKjU1H\nylmLXl4UIoXGpiPlrEUvLwqRQmPTkXLWopcXhUihselIOWvRy4tCpNDYdKSctejlRSFSaGw6\nUs5a9PKiECk0Nh0pZy16eVGIFBqbjpSzFr28KEQKjU1HylmLXl4UIoXGpiPlrEUvLwqRQmPT\nkXLWopcXhUihselIOWvRy4tCpNDYdKSctejlRSFSaGw6Us5a9PKiECk0Nh0pZy16eVGIFBqb\njpSzFr28KEQKjU1HylmLXl4UIoXGpiPlrEUvLwqRQmPTkXLWopcXhUihselIOWvRy4tCpNDY\ndKSctejlRSFSaGw6Us5a9PKiECk0Nh0pZy16eVGIFBqbjpSzFr28KEQKjU1HylmLXl4UIoXG\npiPlrEUvLwqRQmPTkXLWopcXhUihselIOWvRy4saR6TL7dfLZf3hlqYjdY5NR8pZi15e1DAi\n3c25CTX/cE/TkTrHpiPlrEUvL2oUkW7XnrtM8w+P7B3pv3+6selIYy6GEyZEjdnrEGsKsd61\nq4r0v/+yB/i3SrgvISnTSSScIu+VZpF+snPt3fRoGsuF3H7N3yONeVfFCROixux1gDPFpBLJ\noBIi+WBC1Ji9DnCmmGQi7bqESD6YEDVmrwOcKcb8OlLD098+kbZ9QiQfTIgas9ch1hRy0guy\nMpMQyQcTosbspVemnD5vEYrLhEg+mBA1Zi/FflvS/b12TpUQyQcTosbspdhvS7qL9IxNJ0Ty\nwYSoMXsp9tuSPCJdN116wqys3Qy6GE6YEDVmL8V+W5JKpJ+EHjlFMuhiOGFC1Ji9FPttSTqR\nbjnBpkEXwwkTosbspdhvS5KKZFApAm2vVYMJUfTyoRBpL8eaNOhiOGFC1Ji9FPttSWaR9lQK\nY1trFWBCFL18KEQy5iCZBl0MJ0yIGrOXYr8tGUCkq+l5cXcGXQwnTIgas5divy0ZQ6QJbGVV\nmDTmYjhhQtSYvRT7bcmAIolUGnQxnDAhasxeiv22ZEiRJCYNuhhOmBA1Zi/FflsypkjX4sMm\nJ2nMxXDChKgxeyn225JhRWo2adDFcMKEqDF7KfbbknFFKqnkIo25GE6YEDVmL8V+WzKySG0m\nDboYTpgQNWYvxX5bMrRITSYNuhhOmBA1Zi/Fflsytkjll2qNpDEXwwkTosbspdhvS0YXKW7S\noIvhhAlRY/ZS7Lclw4v0iNulQRfDCROixuyl2G9L3kYkt0mDLoYTJkSN2Uux35a8j0hekwZd\nDCdMiBqzl2K/LXkjka6ItIYJUWP2Uuy3Je8k0hWRVjAhasxeiv225K1Eus5c2iONuRhOmBA1\nZi/FflvybiJdzRelQRfDCROixuyl2G9L3k6kq/WqNOhiOGFC1Ji9FPttyTuKZDNp0MVwwoSo\nMXsp9tuStxTJZNKgi+GECVFj9lLstyVvL1LVpEEXwwkTosbspdhvS95TJMsladDFcMKEqDF7\nKfbbkjcVyXBNGnQxnDAhasxeiv225F1Fuk5dKpPGXAwnTIgas5divy15Y5F2TBp0MZwwIWrM\nXor9tuQzRCqZNOhiOGFC1Ji9FPttyTuLtG3SoIvhhAlRY/ZS7Lclby3S5p27QRfDCROixuyl\n2G9L3lukiUpr0piL4YQJUWP2Uuy3Je8uUv158EEXwwkTosbspdhvS95epKpJgy6GEyZEjdlL\nsd+WvL9INZN616qi6OVDIVJobIHbIJIENWYvxX5b8gkilU3qX6uCopcPhUihsYVuVTIpQa0y\nil4+FCKFxha7WeFxUoZaRRS9fChECo0tesOVSjlqFVD08qEQKTS28C2XJiWptUbRy4dCpNDY\n4jddmJSl1gpFLx8KkUJja7gtIjWhxuyl2G9LPkmk+ZtYE9Wao+jlQyFSaGxtN5+YlKnWDEUv\nHwqRQmNrvD0ihVFj9lLstyWfJtIVkaKoMXsp9tuSjxPpaVKyWi8UvXwoRAqNrZnwvHOXq9YL\nRS8fCpFCY2smFN8K3hoW1gkTohApNLZ2BCLFUGP2Uuy3JR8o0hEmsbBOmBCFSKGxKSB6k1hY\nJ0yIQqTQ2BQQ/cMkFtYJE6IQKTQ2CWX900mNYWGdMCEKkUJj02DUJrGwTpgQhUihsYk4YpNY\nWCdMiEKk0NhUIETyosbspdhvSz5WJK1JLKwTJkQhUmhsOpLSJBbWCROiECk0Nh3pG5FcqDF7\nKfbbkk8WSfjCLAvrhAlRiBQam470LXyLAwvrhAlRiBQam470g1KZxMI6YUIUIoXGpiPNRGo0\niYV1woQoRAqNTUf6RYlMYmGdMCEKkUJj05EWIjWZxMI6YUIUIoXGpiMtRWoxiYV1woQoRAqN\nTUe6oRDJihqzl2K/Lfl4kTQmsbBOmBCFSKGx6UgvVLtJLKwTJkQhUmhsOtIE1WwSC+uECVGI\nFBqbjjRFtZrEwjphQhQihcamI81QiLSLGrOXYr8tkYj0BnleknoXIWOGK9I9bXfu+H9+J0yI\neqMrUtOROsemIy1QiLSDGrOXYr8tQaRnWkxiYZ0wIQqRQmPTkRDJhxqzl2K/LUGkVxpMYmGd\nMCEKkUJj05EQyYcas5divy1BpFcQaRM1Zi/FfluCSJPE377KwjphQhQihcamI22J5DWJhXXC\nhChECo1NRyqgwiaxsE6YEIVIobHpSCXUv38xl1hYJ0yIQqTQ2HSkXZEcJrGwTpgQhUihselI\nRVTQJBbWCROiECk0Nh2pjEKkGmrMXor9tgSR1gmYxMI6YUIUIoXGpiPVUX6TWFgnTIhCpNDY\ndKQNFCIVUGP2Uuy3JYhUCSItUWP2Uuy3JYhUie+SxMI6YUIUIoXGpiMhkg81Zi/FfluCSLUg\n0gI1Zi/FfluCSLUg0gI1Zi/FfluCSLW47tuxsE6YEIVIobHpSIjkQ43ZS7HfliBSNYg0R43Z\nS7HfliBSNYg0R43ZS7HfliBSNZ77diysEyZEIVJobDqSVSSLSSysEyZEIVJobDrSLgqRZqgx\neyn22xJEqgeRZqgxeyn22xJE2ggiTVFj9lLstyWItBH7JYmFdcKEKEQKjU1HQiQfasxeiv22\nBJE2gkhT1Ji9FPttCSJtxWwSC+uECVGIFBqbjoRIPtSYvRT7bQkibcZqEgvrhAlRiBQam46E\nSD7UmL0U+20JIm0HkZ6oMXsp9tsSRNoOIj1RY/ZS7LcliLQdRHqixuyl2G9LEGk7xgdJLKwT\nJkQhUmhsOhIi+VBj9lLstyWItBNEeqDG7KXYb0sQaSeI9ECN2Uux35Yg0k5s9+1YWCdMiEKk\n0Nh0JETyocbspdhvSxBpJ7a/uYGFdcKEKEQKjU1HMqIQ6Y4as5divy1BpL2YLkksrBMmRCFS\naGw6khVlMYmFdcKEKEQKjU1Hcou0YRIL64QJUYgUGpuOZEYZTGJhnTAhCpFCY9OR7Kh9k1hY\nJ0yIQqTQ2HSkiEg1k1hYJ0yIQqTQ2HQkB2rXJBbWCROiECk0Nh3Jg0KkUXsp9tsSRDJm2yQW\n1gkTohApNDYdyYlCJB0LkcppOlLn2HQkRPKhxuyl2G9LEMmcLZNYWCdMiEKk0Nh0JETyocbs\npdhvSxDJng2TWFgnTIhCpNDYdCRE8qHG7KXYb0sQyZ6NZ8BZWCdMiEKk0Nh0JD+qbhIL64QJ\nUYgUGpuOFEAhkgomRCFSaGw6EiL5UGP2Uuy3JYjkCiKJYEIUIoXGpiM1iLQyiYV1wkScry9E\nigWRfDAhKlGvr1cQKZbOItVMetOFnaGy9PqaB5FCQSQfTIjK0OurlO3/0klBJF8QSQOL3azo\nESJF0lukiklvtrBFVM9eZYMQKZ40Is1NepeF3UJ16rXp0O+DpO3/0klBJGfKf3vDGyzsLqpL\nr12NeLIhlu4ilS9Jwy+sAdWjl+H+HCKF0l+kokmjL6wFdXov20MiRAolgUil5xuGXlgj6uRe\nNo0QKRhE8sGEqHN7WT1CpFgyiFQwadyFtaNO6bX7kMjbS7HfliBSIIjUDCt/2SuRoZdivy2x\ninT5zf3j7fPXHzYdqS+pRHqZNNjChlAH99p/fi7U6whpSnFdkS6//5v/7idNR+pLCpHWl6SB\nFjaMOrZXzKL9XkpZtuIR6fLU5/L85TdNR+oLIvlgQtShvaIajSnS9+2u3Uyk//2XA3olz12k\n3jWGTfV+3Cy9W3riEOlhz38mIRIiNeT9NPKLdPvNp9+1W923G+YuVANKxtr2R9xL6MpmECkU\nRIpHadF+L6Erm7GLdJn88vEiLU3KuLC/MCFKwpJejCy9tLrU4xTp/hjp05/+Xr2UlG5hHzAh\nSsDS3qkz9VILU4tbpPsrsR/9guxPEMmfAzQaUKSNNB2pc2w6UiMKkTyRebMKIsXGpiNpRLqZ\nlGZhlzAhqoGlu/64eyn22xJEigaRbCnejTuvl2K/LUGkcBBpJxuPhhCpmKYjdY5NRxKJ9GsS\nIq2z9aQCIhXTdKTOselIKpF+TEKkefaem0OkYpqO1Dk2HakZhUi17HmESOU0HalzbDoSIvlQ\nVtaOQyf3Uuy3JYjUkJdJiPTMvkWn9lLstyWI1JDXoyREemRfonN7KfbbEkRqCSIts+/QDSar\nhUjBselIQpH+IdJvDNeis3sp9tsSRGoKIj1iuUPXo5divy1BpKYg0m9Mj4s69LoiUhWmI0lQ\nyx+VbU/ScW30sj2/cH6v+x+fE0RqCyKF3tqNSMU0HalzbDoSIvlQFVZAI0SqpOlInWPTkTSo\nTxcp5BEildN0pM6x6UiI5EMVWTGPEKmcpiN1jk1HQiQfqsQKOHSDqVohUnRsOpJUJJ1JScdV\n7BWz6IpIlTQdqXNsOhIi+VBrVtgjRCqn6UidY9ORtCLJTEo6rkKvuEeIVE7TkTrHpiOJUB8r\nUoNHiFRO05E6x6YjqVBik5KOa9WrxSNEKqfpSJ1j05EQyYeasYLP1j1holJXRIqOTUeSoT5Q\npDaNEKmSpiN1jk1HQiQfasJq9QiRymk6UufYdCS1SCKTko5r2qvZI0Qqp+lInWPTkRDJh3qw\nGh8e3WCiUldEio5NR9KhpCYlHdezl8IjRCqn6UidY9OREMmHurEEFl0RqZKmI3WOTUdCJB/q\nl6W4HF0RqZKmI3WOTUcS1voYkUQeIVI5TUfqHJuOpBdJYlLScf32UnmESOU0HalzbDoSIvlQ\n35KnGe4wSaUbCpFCY9ORlLU+QKQvoUeIVE7TkTrHpiMdIJLCpKTjknqESOU0HalzbDrSESIJ\nTMo5LqlGiFRJ05E6x6YjIZI9Yo8QqZymI3WOTUeS1npvkcQaIVIlTUfqHJuOpK0lMynfuNSX\noysiVdJ0pM6x6UjHiNRsUrpxHeARIpXTdKTOselI4lrvKtIBGiFSJU1H6hybjnSQSK0mJRvX\nIR4hUjlNR+ocm46kriUyKdm4Xh7l6vVCIVJobDoSIhky8ShVrwkKkUJj05HktTQm5RrX625d\nrl4TFCKFxqYj6Wu9n0iTh0epek1RiBQam46ESLv5QiRrEEkDu31QmJRoXLPn6xL1mqMQKTQ2\nHQmR9oJI9iCSBnb78F4ifSGSPYikgd0/vqdIN1SaXgsUIoXGpiMdKFKLSWnGtXhHQ5peSxQi\nhcamIyHSdhbvDErTa4lCpNDYdCRE2szyHXZZeq1QiBQam450SK12k5KMa/VW1SS91ihECo1N\nR0Kkraze8p2k1xqFSKGx6UjH1HoTkdY/O5GjVwGFSKGx6UiItJH1zyDl6FVAIVJobDrSoSLF\nTUoxrsIP86XoVUIhUmhsOhIi1VP4odgUvUooRAqNTUc6VqSwSRnGVbggpehVRCFSaGw60kG1\n3kmkGSpBryIKkUJj05GOqjW+SCWPMvQqoxApNDYd6WCRoib1HtfrPd+I5AgiaWCT3w8t0lfN\no6ynEZGCY9OREGmdL0QKBpE0sOknbSZ1HNfXhkdZTyMiBcemIx0uUsykfuPa9CjraUSk4Nh0\nJESaZ9ujrKcRkYJj05GOqzWiSF87HmU9jYgUHJuOdGCtFpP6jGvHomva04hIwbHpSIj0yr5H\nWU8jIgXHpiOdIFLEpB7jMniU9TQiUnBsOtKRtRpMOn9cBouuaU8jIgXHpiMh0i02j7KeRkQK\njk1HOkUkv0lnj8voUdbTiEjBselIh9aKm3TyuIwapT2NiBQcm450bK3xRNpD5TyNiBQcm450\ncK0RRLLerftF5TyNbyUSKeUuUu8aG5l61LvL6OGKpIEVvpb+iuS5HqU9jW91RWo6UufYdCRE\n8niU9TQiUnBsOtJJInlNOm1cvgtS1tOISMGx6UiI5PEo62lEpODYdKQPF8npUdbTiEjBselI\nh9eKmXTSuHwWXdOeRkQKjk1HQiSPR1lPIyIFx6YjnSaSz6RzxuW+IGU9jYgUHJuOdHytvCL5\nPcp6GhEpODYd6TyRXCadMS7n8ww3VM7TiEjBselIJ9SKmHRCr4hHWU8jIgXHpiOdKZLDpDNF\ncqFynkZECo5NRzqjVsCk43uFPMp6GhEpODYd6VNFinmU9TQiUnBsOtIptfKJFHqAdE17GhEp\nODYd6ZxabpMO7hX1KOtpRKTg2HSkk2p5TTpLJDcq52lEpODYdKSzaqUSKexR1tOISMGx6Uif\nKFLco6ynEZGCY9ORTqvlM+nIXuEHSNe0pxGRgmPTkT5PpBaPsp5GRAqOTUf6YJFCqJynEZGC\nY9ORzquVTaQYKudpRKTg2HSk00WymXRcryaPsp5GRAqOTUc6XySTSYf1arpjl/Y0IlJwbDrS\nibU8Jh3Vq9GjrKcRkYJj05HOrOUw6aBerR5lPY2IFBybjoRIPlTO04hIwbHpSKfWspt0TK9m\nj7KeRkQKjk1HOrdWEpEaUDlPIyIFx6YjnVzLatKhIrWgcp5GRAqOTUdCJB8q52lEpODYdKRP\nEqn9nl3W04hIwbHpSGfXMpqESE7UASL9+VP+5E9dF0TSwAzfk0CkJlTO04hIwbHpSL1E2jEJ\nkZyoo+/aIdIO6fRaNpOO6CW4Z5f1NCJScGw60ueI1P5q7DXtaWwR6ebF769/fnL75Pa/xZcm\n3/r8DJFupPNrmUzS92r6wdgXKudp1Ij052bLVKTJl/5MnHp8jkhPUodanUVqQ+U8jSqRJkZV\nv/T6+vpOHiJpYLZvs1ySjhOpEZXzNDY9RnqJ8f39vCJNNJl86fm9iLQk9ahlMOkwkVpROU+j\nSKTH/bXlw6aSSMX7doikgRm/79++SYjkRClE+rO+3iy+tLgirVGIJIFZv7GDSJp7dllPY9vT\n33/uBhQeAS2+hEhVUp9auyape2meakh7GhtF+vPHINLiyYbnp4j0S+pUa88kcS+VR1lPY+ML\nss/Lzevp74kwG09/I9KT1KvWjkmI5EQpRLprs3yy4fElXpDdIn2ESDKPsp5GfowiODYdqVut\nHiIJUDlPIyIFx6Yj9au1aZK0l86jrKcRkYJj05H6i1Q0CZGcKEQKjU1H+gCRdI+Q0p5GRAqO\nTUfqWOtskRSspKcRkYJj05F61towSdhL6VHW04hIwbHpSIjkStLTiEjBselIXWvVTdL1knqU\n9TTmFGnjZ9K30nSkzrHpSO8ukvKZhmvveW2gECk0Nh0phUhrk1S9xB5lPY05RQre1Ws6UufY\ndKS+tU4TSUPrPq86KqNIf/78Kb+3dTtNR+ocm46UQ6SVSaJe6gtS73nVUXGR/q8eRLKTOtc6\nViS5R93nVUVlFCmYpiN1jk1HQiRXes+rikKk0Nh0pN61KiZJeuk96j+vGiqnSNO/N8WcpiN1\njk1H6l2r8ihJ0evp0XkL64QJUSlF+oNIMVjgNseJ9LoeIdLW3h75ZMPtKzzZ4IVFbnSUSF+I\nNP1jRCrCdKT+tYomNff6QqTZHyNSEaYj9a9VvHPX2mvqESIpRDK9HlQUideR3LDQrf6VTGrs\nNfMIkQQilf8eu22Ran+xMSLtwWI3K5nU1mvuESIpRCqJsidS5R9/QaQ9WPB2B4rU1KuQFPMq\noVpE+ldOQaS98IKsBha94dqkpl4LjxBJIZLpwoJIGlj4lseI1N5rlSTzWqOOFinyZMN35e9j\nRaQ9WPiWiNSKOlykmxY7DvBkgwYWv+kRIil6LZNlXivUOY+RfCI9/65wnv52wuI3VYq0eqsq\nIvUSaf4Rkayw+E2X9+0aeq3f841IApECryMhUhQWv+m/hUnxXoWfnUCkbZGOfGcDIkVg8Zuq\nRFq8FNvca5ks81qhMr7XrvZvn++k6UidY9ORstRamBTsVfQIkTqI9GcWRPLBWm48NynWq+wR\nIiFSDaYjpaklFUnYa4FKM68FKptIDWk6UufYdKQ8tWYmhXqVr0eIlFCky09uH79fHxDpDmu7\n+dSkSK+aR4jU8+9sKN+1uzw/XF4fEOkBa7x9m0hVjxCp20/I7oh0uxo9PiDSE9Z4+8klCZGc\nqJwiVb/xds9uLdL//ov7v0qWeYgUue3LI3UrYs7yrl31G38suiDSYVGIpO5E7HH9axQlkX7S\ndO11Xsh1pFy1nvft/L1q9+skvV6oXPN6oTLetdt4sgGRNmHNhLhIGx4hUronG+YGIdIS1kwI\ni1R9okHT64XKNa8XKqdI1W+8PUbi6e8KrB0RFGnTI0TqdteuLtL8lVhekJ3D2hGv5+1cN0Ok\nlCLxT18GYe2ImEjbHiFSusdIiLQNa0c8X5MNiXRcrycq2byeKEQKjU1HylYLkWKojCIF03Sk\nzrHpSNlqRUTa8QiRBCLxN63ukNLVWvzMuSE7j5AQqZdI3LWLwhSQ5d/esJtdjxBpR6S/5SCS\nh5SuFiKFUBlFct30laYjdY5NR8pXy2nSvkeI1Fckr0lNR+ocm46UsJbPpH2PEKm3SNy1c8I0\nGJdIhgsSIiFSDaYjZazlMcngESLxZEMNpiOlrGUXyXJBQqTeIrk8QiRhLatIJo8QiRdkazAd\nKWct6yXJ5BEiIVINpiPlrGUUyeYRInV4rx1/ZXELTEZCJC8KkUJj05Fy1jKaZPMIkbq++5sn\nG/wwHQqRnKi0Irk9QiTpYlhMMt6zQ6R+Ivk1QqRjRNowyfbct7hX1nmlFMn/GtJPmo7UOTYd\nKWetiUh1k6weIVK/HzV3E74RSbwY+yYh0gSVT6SgRoikXgyrSGf3EsKEqLhIykzUCXuESOrF\nQCQHKptIvI7UAhOifljbIpnv2SESL8jWYDpSzlozkSomIdIUlU2khjQdqXNsOlLOWg6ROvSS\nwYQoRAqNTUfKWevea9MkRJqiECk0Nh0pZ62FSCWT7PfsEAmRajAdKWetRy9EsqIQKTQ2HSln\nraVIBZPsHiESItVgOlLOWs9eVZMcFyREQqQaTEfKWQuRvKijRbK9HoRIGpgQ9WBVRPJ4hEiq\nv7Nh1yRE0sCEqCeraJL5JyiO66WACVEtIn2VU/jLTxCpSspZa0ckn0eIhEg1mI6Us9a0V8Ek\nn0eIpBKJu3Z1Us5a2yI5L0iIpBCJJxu2STlrlUR6meT0CJG4ItVgOlLOWrNeS5O8HiESj5Fq\nMB0pZ615r7lJ3jt2iIRIVZiOlLNWRaRfk9weIZLqdSREqpNy1qqJ9GMSIpVQiBQam46Us9ai\nV0mkDL1aYUIU77ULjU1Hyllr2QuR9lCIFBqbjpSz1qrXyyT/PTtEQqQaTEfKWWvd62kSIhVR\niBQam46Us1ahV9wjREKkGkxHylmrKpL/RaSje7XAhChECo1NR8pZq9QrfEFCJESqwXSknLWK\nvaIXJERCpBpMR8pZqyJSzCNE2hRJGUTSwISokkhfiFRFIVJobDpSzlrlXkGPEAmRajAdKWet\nbZHcMEGhB2qkeU3/+JwgkgYmRBVYD4/2/q3zNUxS6YYaaF6zPz4niKSBCVFr1hcibaAQKTQ2\nHSlnrUKvL0TaQiFSaGw6Us5a615Tj7wmIZJivy1BJA1MiFqwvhBpG4VIobHpSDlrbYhU/qvA\ne/VqgwlRiBQam46Us9ay19Ijn0mIpNhvSxBJAxOiZqy1Ry6TEEmx35YgkgYmRBVFui7/RqHO\nvVphQhQihcamI+WsZRHJbhIiKfbbEkTSwISoKWv21iBEKqIQKTQ2HSlnrQ2RNv+p83N7NcOE\nKEQKjU1Hyllr1mv1XlWvSoik2G9LEEkDE6JMIhlNQiTFfluCSBqYEPVirX8KyfuEAyIp9tsS\nRNLAhKgnq/TTfE6TEEmx35YgkgYmRG2K5DQJkRT7bQkiaWBC1INV9MhpEiIp9tsSRNLAhKil\nSMtvQKQ5CpFCY9ORctZ69qp55HoOHJEU+20JImlgQtSuSI5rEiIp9tsSRNLAhKgba8Mjh0mI\npNhvSxBJAxOiflmVZxoesZqESIr9tgSRNDAhaiZS5ZusT90hkmK/LUEkDUyIsohkNQmRFPtt\nCSJpYELUVKT6tyHSA4VIobHpSDlr3XrtXZCuRpMQSbHflkhEIuo8RNr6nqdIZ5UiG+GKpIEJ\nUd/Wf3yCK9Iv6o2uSE1H6hybjpSz1kyknW/dNwmRFPttCSJpYELUt/lfQ9p/lIRIiv22BJE0\nMCHq2+qRwSREUuy3JYikgQlREZFqJiGSYr8tQSQNTIj6Nnu0bxIiKfbbEkTSwISob7tHu29w\nQCTFfluCSBqYEOW4IO2ahEiK/bYEkTQwIcol0s6dO0RS7LcliKSBCVEuj3ZMQiTFfluCSBqY\nDuW7IP1kwyREUuy3JYikgelQfpE2TEIkxX5bgkgamIxkfxHplfq9O0RS7LcliKSByUgBjzYu\nSYik2G9LEEkDU4FCHtXfvopIiv22BJE0MBUoKFLtmoRIiv22BJE0MBEn8gjpN5WHSYik2G9L\nEEkDE3GiHtXe4YBIiv22BJE0MBEnLFLFJERS7LcliMw84k4AABYzSURBVKSBaTBxjyomIZJi\nvy1BJA1Mg2kRqWgSIin22xJE0sA0mCaRSk84IJJivy1BJA1MQmnzqGQSIin22xJE0sAklFaR\n1iYhkmK/LUEkDUwBCb+I9MrSJERS7LcliKSBKSDtHq3e4YBIiv22BJE0MAFDcEG6Lk1CJMV+\nW4JIGpiA8fCokTUzCZEU+20JImlgAoZIpJlJiKTYb0sQSQNrRzw9kon0XxBJsd+WIJIG1o6Q\niTQzafffPTcn27yeKEQKjU1HylXrSyfSMSYlm9cLhUihselIqWq9nrJT9DrCpFzzmqAQKTQ2\nHSlVrddT35JeiPT643OCSBpY4+2/xCIdYFKqeU1RiBQam46UqNb0tVhRL7lJmeY1QyFSaGw6\nUqJaB4h0Nf675+ZkmtcMhUihselIiWpN3xwk7CU1KdO8ZihECo1NR8pTa/YmO2UvREKkGkxH\nylPrMJGUJiWa1xyFSKGx6Uhpas3f9a3tpTMpz7wWKEQKjU1HylLr6xSRmlVKM68lCpFCY9OR\nktRa/hiSuJfMpCzzWqEQKTQ2HSlJrYNFkr2glGVeKxQihcamI+Wotfq5WHkvkUlJ5rVGIVJo\nbDpSilrrny/X99KYlGNeBRQihcamI2Wo9XWGSJrHSSnmVUIhUmhsOlKGWueIJDEpxbxKKEQK\njU1HSlCr4NExvRDp6CCSBha7WcGjg3q1m5RhXkUUIoXGpiMlqFXw6LBerSZlmFcRhUihselI\n/WuVPDpepKBJCeZVRiFSaGw6Uv9ap4rU+ix4gnmVUYgUGpuO1L1W0aNzRAr41H9eFRQihcam\nI3WvdbJIjSb1n1cFhUihselIvWuVnrK7HtqryaTu86qhECk0Nh2pd62yR8f2QqSjgkgamP8m\nFY8O7hU3qfe8qihECo1NR/pIkeJ373rPq4pCpNDYdKS+tSqPkE7pFTIp6WlEpODYdKQcIq1R\nJ/XyupT0NCJScGw60oeL5DUp6WlEpODYdKSutaoendjLdw8v6WlEpODYdKSPF6n0xEPdp6Sn\nEZGCY9ORMohUQp3Zy2FS0tOISMGx6Ug9a9UvSGf3MpuU9DQiUnBsOhIi/cZqUtLTiEjBselI\niPSIyaSkpxGRgmPTkTrW2vCoSy+LSUlPIyIFx6Yj9atVfVfDL6pXrx2Tkp5GRAqOTUfqL1IZ\n1XVcdZeSnkZECo5NR+pWa/OC1Htc1ctS0tOISMGx6Ui9am171H1cNZN696qiECk0Nh2pU62v\n5CLVTOreq4ZCpNDYdKQ+tb7GEullUvdeNRQihcamI/UWqYbqPq6ySf17VVCIFBqbjtRZpCqq\n/7iKJiXoVUYhUmhsOlJXkTZQSca1VClLrxUKkUJj05G61Nq7HiUa18KkNL2WKEQKjU1HQqSd\nzE3K02uBQqTQ2HSkniJtofKMC5HsQSQNzPh9+xekVOOampSp1wyFSKGx6UiItJ/ps3eZek1R\no4l0uVxuvz4+vv6o6UidY9OROtTafe77mm1ctXfeNeaDRbr8yPP74fHZ68+ajtQ5Nh2po0ib\nqFzjOsakzxbppc/l+ctvmo7UOTYdqZ9I26hk46q89a4tnyvS5SHS5Xlduon0v/9yULe3y+OC\n1LuHL/+W6V0oZTxPNtwV+u8DIoUypkiYZIlDpOfTCxfu2q1glm+yPNWQclwrk5qJn3vXbuoN\nIq1hlm8yeZRzXGqTPlik2bMMiLSEWb7J5FHScYlN+lyRXnfrng+VXn/YdKTOselIZ9eyXZCy\njqvyF4UHjfpgke6vxN5fj+UF2QXM8D2Di/RfL51LnyvSZpqO1Dk2HenkWkaPso7r3kukEiIV\n03SkzrHpSIjkQ91YGpMQqZimI3WOTUc6t5btue9r2nE9e1Xu3/lcQqRimo7UOTYdqY9I+6ic\n45r0qqnk8AmRimk6UufYdKRTa5kvSFnHNeu1qZLFJUQqpulInWPTkc6sZfco67iKvcImIVIx\nTUfqHJuO1EMkCyrnuGq9YtclRCqm6UidY9OREMmHqrEid/MQqZimI3WOTUc6sZbDo6zj2ui1\n/YipKBUiFdN0pM6x6Ujn1XI8Qko7rp1eVpn+ndwLkWowHem0Wl8fINI9Tp9O6KXYb0sQSQPb\n+sMPEsmhkkonRIqNTUc6q5bPo6zjsvdyqdTuFyKFMp5ITo+yjsvb6zSTECmU4UTyepR1XLFe\nJ6iESKGMK5IZlXNczb2aLlEbniFSKKOJ5PYo67gkvY6QCZFCGVUkByrnuKS9lD4hUiiDieT3\nKOu4DumFSPPsDC3nZiCSE3VWL0SqHop/0HWYjnRCrYBHWcd1Zi9EqhyKd8xbMB0JkXyoM3sh\nUvFQfEPehulIx9eKeJR1XKl6IVJzEMkHE6KS9UKkpiCSDyZEjdlLsd+WIJIGVvpiyKOs4xq1\nl2K/LUEkDaz0RUSqw4QoRAqNTUc6ulbMo6zjGrWXYr8tQSQNrPA1RNqACVGIFBqbjoRIPtSY\nvRT7bQkiaWCFr8U8yjquUXsp9tsSRNLA1l8KXpCyjmvUXor9tgSRNLD1lxBpCyZEIVJobDrS\nsbX8P9H3QOUc16i9FPttCSJpYKuvRD3KOq5Reyn22xJE0sCWXwh7lHVco/ZS7LcliKSBLT4P\n37FLO65Reyn22xJE0sDmnzZ4lHVco/ZS7LcliKSBzT9t8CjruEbtpdhvSxBJA5t/uhLp70/6\n92pCjdlLsd+WIJIGNvvs4dHfhz9/n+naqw01Zi/FfluCSBrY7LO5SH8nHhlcSjquUXsp9tsS\nRNLAZp/NPVqljvn506TjynoaESk4Nh3psFpfOyLVTLr/6ffV/HjK16sRlfM0IlJwbDrSUbW+\ndkUqa7L/LX7BEEmx35Ygkgb2+u2XQaSFLL+/3Zdt537hTq/WJD2NiBQcm450TC2HR7enIdYK\nVUTauprt92pO0tOISMGx6UiH1HJ5tC3ZIlXBTL3ak/Q0IlJwbDrSEbW2PLrWrj0bLj3E2TLM\n0ktxiDlPIyIFx6YjHSpS5QrjNOlx0zXI2UtxiDlPIyIFx6YjHVBrfUFaPbSJqbQ0y9lLcog5\nTyMiBcemI+lrlTy6Lp61lohkVgmRFPttCSJpYL+2zD0qf6dIJOs7YIWHmPM0IlJwbDqStNbP\nbk89qn/rth5ikxBJsd+WIJIkS4+21vyuQUmgq+OKZeqFSIr9tgSRFJl5tH1Buj4fM1Wk+P42\nymQphkiK/bYEkQT56xNpcrPC235+elXdQaQ1CpFCY9ORtCL5Park2at0DfJekhBJsd+WIFJz\n/s4vSObHL4Zec3le/zFEmqAQKTQ2HUmE+jsXyXyvy95rRvSZhEiK/bYEkRqj9sjSy24SIin2\n2xJEakxBpMN7IdIUhUihselICtRzpV8indDLfucOkRT7bQkiteTvUqR2pKnX4nm8DZigkKOX\nHSZEIVJobDpSO2rl0fki7VyYEEmx35YgUjx/+4u0dxcPkRT7bQkihfMXkeIwIQqRQmPTkYQi\nCT2y9VqZ9Hiz0UIqRFLstyWIFE5PkTbf1zqFKRp5ellhQhQihcamIzWippt7vkhGkxBJsd+W\nIFIwR3lk7rVl0tMnRFLstyWIFMt0W7+6iGQyCZEU+20JIoUyXVatR45eBpMQSbHfliBSKIX7\ndeeLdN13CZEU+20JIkVy4AXJ32vn3p0qSU8jIgXHpiPFUbM9FXsU6HWOSUlPIyIFx6YjhVGH\nehTqdW9zqEpJTyMiBcemIwVRixVNIdIkh6mU9DQiUnBsOlIMtVhQuUfKF4qlLiU9jYgUHJuO\nFEItlzOfSI+oTUp6GhEpODYd6c1FUpuU9DQiUnBsOlIEtVzNh0cpX69BpOtgIn1MZpv53+fP\n15B6Fyvn77IvOSxckRypXY+SXpG016Skp/GtrkhNR+ocm47UItLt89ebGpKKdP3WmZT0NCJS\ncGw6UsM7CG6fT94clFak+D8/u0LlPI2IFBybjmRHLZ7/WnuUWKTqPx/jRuU8jYgUHJuOFP5x\nhdtXv0YRyfH3SW6icp5GRAqOTUeyogwe5RbJ8Vccb6FynkZECo5NRwqKdPvi4qeQcoskuSYl\nPY2IFBybjmREFS9Iy5/mSy6S4nFS0tOISMGx6Ug21NKj2RuDnm8Nyi6SwKSkpxGRgmPTkWIi\n/Xxt5dFQIgVNSnoaESk4Nh0pJNLPl9YejSVSTKWkpxGRgmPTkSIi/Xzla3yRIiYlPY2IFByb\njmRCrTev5NFwIgVMSnoaESk4Nh0pJlLRo/witV+Skp5GRAqOTUeyoNaLV9JoBJGaTUp6GhEp\nODYdySPS8wvF69EQIrWalPQ0IlJwbDqSAWW9IA0h0i1hk5KeRkQKjk1Hcoh0/7T8+Ehb6/Bx\nRa9KSU8jIgXHpiPZRbp/VvdoJJGiJiU9jYgUHJuOtI+ar9vXm4gUNCnpaUSk4Nh0pF3UfNu2\nPBpLpNiLSklPIyIFx6YjmUX6/WTTo8FECpmU9DQiUnBsOtIearZm2x6NJlLpX0VP0SuCQqTQ\n2HSkqEjH1jppXH6Tkp5GRAqOTUfyiLSt0YAi+U1KehoRKTg2HWkH5fFoQJG2/v2XolZJTyMi\nBcemIxlFuho8GlGkuknlC1TS04hIwbHpSNWnhh8f7R4NKVJNpcpdvaSnEZGCY9OR6i9WPn/z\nd/avlm/82y1jivT7/xolmxDJHUSaZb1OJo9GFekWm0lJTyMiBcemI+28odPj0dgi2UxKehoR\nKTg2HWmNKuySzaPBRbpa/uXZpKcRkYJj05FWqLhH44u0b1LS04hIwbHpSAaRrB69gUg/2TQp\n6WlEpODYdKQlakOj3X9r+T1Eej1VWTAp6WlEpODYdKQFasujTxHp+Vps4aqU9DQiUnBsOtIc\n1eTR24i0ej36ZVLS04hIwbHpSHfU+lHBXCKLR+8j0jOrqSTptUYhUmhsOtINVXp47fbo/UWa\nP17q2GuNQqTQ2HSkX9S+RyfXSjOukkkylxCpmKYjdY5NRyqKdA159I4i1UyS2IRIxTQdqXNs\nOtIPSuTRW4pUfztr+6UJkYppOlLn2HSkb51H7ynSgSYhUjFNR+ocm470XfDI/zSDulaqcR1m\nEiIV03SkzrHpSN+bl6NetdKNC5FM+ViRZivR6tE7i3SISohUTNOROscmoRT2IXifTlrrhso3\nrhtqeQVvUQmRimk6UufYFBC5R58h0rX8sluGXlt/fE4+UaTCLjR69DEiqd7wgEjFNB2pc2zN\nhNIatHr0OSLdXo9tdgmRimk6UufYWgGHePRBIt3SahIiFdN0pM6xNdy2crck+tqRqtYSlWVc\nS9SCVZqm2SdEKqbpSJ1jC9+ydtoFGn2kSLVnxbv3Wv7xOfkMkYpn/Ht5LYpaFK5VRvUfVxm1\nYlVMssiESMU0HalzbJEbFc/197fOo88UqcEkRCqm6UidYwvcpnyqhRp9qkjX7Xfj9ew1+eNz\n8gEiFU/yUqM2jz5XpJ33tdaUQqRimo7UOTb3Lc7w6INF2japdnFCpGKajtQ5Nu8NCqd0ZVGr\nRh8t0k9MMk2FQqRimo7UOTbn9xfO4wEefbpIj5hVQqRimo7UOTbXd5fuXKgdCtTaRg0s0tV6\nTw+Rimk6UufYPN98mkeINM+uSYhUTNOROsfm+N49jzrV2kO9gUiG69Ls2w7spdhvS95VpMJZ\nuy4uR++wsLswISr0PI/jiYigTogUG5vt2ywevc3CbsGEqAaW2ySHVogUSuQR7+1ry8dG77iw\na5gQ1c6yqeS8RCFSKCbWrkeIFEAJWMYrk8skRAplk/V3/a/dr180ejzN8M4L+4IJUSqWWCVE\nCqXOKpyFskX3p+vefmF/YUKUjKW9LCFSKFWW2aLH097vv7DXrL2WjjSZhEihzFn14e9YdGyt\nNtRn9tqTqWoUIoUyY7ktmr4G+6ELG0cd3GvHpOrVCZFCebG8Di08+tiFDaNO6DW1xaLSz+eI\nFMqTFdBo/p6gT17YEOqsh7qT39ZM+jt9ehaRInmwHBadWUuBotczNZfq16j1f+mkjCmSWaNz\na2lQ9HoFkWaHYpyaYaY/Z8B2Hdq2SFfrDhOi6DUJIk0PxTq1wBDDHrGwXlTPXlsrwGOk8Oh2\nLbL8B1hYJypBr6JHiGSc1O6TcW6H3LX2YUIUvepZ7wavI+3PKKCR46deUyxGCUWv7TwVuqMQ\nqTieBol8PzueZjGWKHr5UJ8n0uv/RqpXn4g+IY2GXQwnTIgas5divy2JinS5XF6f7B/pjjkN\n+vgNsp0AJ0yIopcPNbRIl9//PbJzpGFDjrNo/wQ4YUIUvXyokUW6PH/5zc6RJrRo/wQ4YUIU\nvXyoNxHpf/9l57vbhQmVJOS8DCFSqCMhJyb9XbslrO1+wOyaL0NxF8oNE6Le5K7dT3aOtFWe\nGcw/6eqEZSgW1g0TohDJI9ATZp/wHmnMxXDChKgxe2k02c8ZT3/nPANJa9HLixpaJO8Lsrqx\n6Ug5a9HLixpbpFmajtQ5Nh0pZy16eVGIFBqbjpSzFr28KEQKjU1HylmLXl4UIoXGpiPlrEUv\nLwqRQmPTkXLWopcXhUihselIOWvRy4tCpNDYdKSctejlRSFSaGw6Us5a9PKiECk0Nh0pZy16\neVGIFBqbjpSzFr28KEQKjU1HylmLXl4UIoXGpiPlrEUvLwqRQmPTkXLWopcXhUihselIOWvR\ny4tCpNDYdKSctejlRSFSaGw6Us5a9PKiECk0Nh0pZy16eVGIFBqbjpSzFr28KEQKjU1HylmL\nXl4UIoXGpiPlrEUvLwqRQmPTkXLWopcXhUihselIOWvRy4tCpNDYdKSctejlRSFSaGw6Us5a\n9PKiECk0Nh0pZy16eVGIFBqbjpSzFr28KEQKjU1HylmLXl4UIoXGpiPlrEUvLwqRQmPTkXLW\nopcXhUihselIOWvRy4t6I5G2s/tvzPZJ0lr0ciZJL0TKFnr5kqQXImULvXxJ0guRsoVeviTp\ndYJIhLx/EIkQQRCJEEEQiRBBEIkQQY4X6XK5HP7fcOZyuZVKVu1W5l4qT7dHnWQjm8+pe6/D\nRbo8FiRRLs8PmardN+JWKk+3Z637hyy1bs3SjOtokaZnIU0ur1/zVLv/P/7vbxN1u0xFylPr\nrk+ecX2mSJfLJWG1lCJN7tplGtkFkRLkdrcgXbXUIiUc2QWREiTF9OfJLNLtN5lqfU8ukQl6\nIVKeIJInl8kvCXp9okh5pj9PYpHy1Zo9+5Gg16c+/T156jRNcj79PX2MlKfW5fUxx7g+9gXZ\n+8feVaa5/59stm6TZ+3y1LosXiHu3ou3CBEiCCIRIggiESIIIhEiCCIRIggiESIIIhEiCCIR\nIggiESIIIvXNn2d+ft+7DQmHc9c3iPQm4dz1DwK9QTiH/YNIbxDOYf88RHrcvfv9/P7F250+\nkj6cpf5ZiPTrzl2gx8Mnkj2cpP6Zi/T65fEnmDRAOEf9s7xr9/oFkYYJ56h/tkTivt0g4Rz1\nz94ViQwQzlT/INIbhDPVP5t37b7RaYhwjvqnLtI3j5FGCeeofzZE4gXZUcJZIkQQRCJEEEQi\nRBBEIkQQRCJEEEQiRBBEIkQQRCJEEEQiRBBEIkQQRCJEEEQiRBBEIkQQRCJEEEQiRJD/B3I6\n5yJ9V17GAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ggplot(sir_out_long,aes(x=time,y=value,colour=variable,group=variable))+\n",
    "  # Add line\n",
    "  geom_line(lwd=2)+\n",
    "  #Add labels\n",
    "  xlab(\"Time\")+ylab(\"Number\")"
   ]
  }
 ],
 "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
}