{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stochastic SEIRD model using odin\n",
    "\n",
    "Author: Thibaut Jombart\n",
    "\n",
    "Date: 2018-10-03"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Requirements\n",
    "This code uses [*odin*](https://github.com/mrc-ide/odin), an **R** package for describing and solving differential equations. The commented commands will install the latest version of the package and its dependencies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#if (!require(\"drat\")) install.packages(\"drat\")\n",
    "#drat:::add(\"mrc-ide\")\n",
    "#install.packages(\"dde\")\n",
    "#install.packages(\"odin\")\n",
    "library(odin)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also use the following helper functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## function to make colors transparent, stolen from adegenet::transp\n",
    "transp <- function (col, alpha = 0.5) {\n",
    "    res <- apply(col2rgb(col), 2,\n",
    "                 function(c) rgb(c[1]/255, c[2]/255, \n",
    "                                 c[3]/255, alpha))\n",
    "    return(res)\n",
    "}\n",
    "\n",
    "## x: instance of odin model\n",
    "## t: time steps\n",
    "## n: number of replicates\n",
    "run_model <- function(x, t = 0:100, n = 1, ...) {\n",
    "    res <- x$run(t, replicate = n, ...)\n",
    "    res <- x$transform_variables(res)\n",
    "    res <- cbind.data.frame(t = res[[1]], res[-1])\n",
    "    attr(res, \"n_compartments\") <- length(x$names) - 1\n",
    "    attr(res, \"n_replicates\") <- n\n",
    "    attr(res, \"compartments\") <- x$names[-1]\n",
    "    class(res) <- c(\"pretty_odin\", class(res))\n",
    "    res\n",
    "}\n",
    "\n",
    "\n",
    "## plot function\n",
    "plot.pretty_odin <- function(x, pal = seird_pal, transparency=FALSE,...) {\n",
    "    ## handle colors\n",
    "    n_compartments <- attr(x, \"n_compartments\")\n",
    "    n_replicates <- attr(x, \"n_replicates\")\n",
    "    col_leg <- pal(n_compartments)\n",
    "    alpha <- max(10 / n_replicates, 0.05)\n",
    "    if(transparency){\n",
    "        col <- rep(transp(col_leg, alpha), each = n_replicates)\n",
    "    }else{\n",
    "        col <- rep(col_leg, each = n_replicates)\n",
    "    }\n",
    "    ## make plot\n",
    "    par(mar = c(4.1, 5.1, 0.5, 0.5), las = 3)\n",
    "    matplot(x[, 1], x[, -1], xlab = \"Time\", ylab = \"Number of individuals\",\n",
    "            type = \"l\", col = col, lty = 1, ...)\n",
    "    legend(\"topright\", lwd = 1, col = col_leg, bty = \"n\",\n",
    "           legend = attr(x, \"compartments\"))\n",
    "}\n",
    "\n",
    "                 \n",
    "## colors\n",
    "seird_col <- c(\"#8c8cd9\", \"#e67300\", \"#d279a6\", \"#ff4d4d\", \"#999966\", \"#660000\")\n",
    "seird_pal <- colorRampPalette(seird_col)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model description\n",
    "\n",
    "The model is specified in *odin* using:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "seird_generator <- odin::odin({\n",
    "  ## Core equations for transitions between compartments:\n",
    "  update(S) <- S - n_SE + n_RS\n",
    "  update(E) <- E + n_SE - n_EI + n_import_E\n",
    "  update(Ir) <- Ir + n_EIr - n_IrR\n",
    "  update(Id) <- Id + n_EId - n_IdD\n",
    "  update(R) <- R + n_IrR - n_RS\n",
    "  update(D) <- D + n_IdD\n",
    "\n",
    "  ## Individual probabilities of transition:\n",
    "  p_SE <- 1 - exp(-beta * I / N)\n",
    "  p_EI <-  1 - exp(-delta)\n",
    "  p_IrR <- 1 - exp(-gamma_R) # Ir to R\n",
    "  p_IdD <- 1 - exp(-gamma_D) # Id to d\n",
    "  p_RS <- 1 - exp(-omega) # R to S\n",
    "\n",
    "\n",
    "  ## Draws from binomial distributions for numbers changing between\n",
    "  ## compartments:\n",
    "  n_SE <- rbinom(S, p_SE)\n",
    "  n_EI <- rbinom(E, p_EI)\n",
    "\n",
    "  n_EIrId[] <- rmultinom(n_EI, p)\n",
    "  p[1] <- 1 - mu\n",
    "  p[2] <- mu\n",
    "  dim(p) <- 2\n",
    "  dim(n_EIrId) <- 2\n",
    "  n_EIr <- n_EIrId[1]\n",
    "  n_EId <- n_EIrId[2]\n",
    "  n_IrR <- rbinom(Ir, p_IrR)\n",
    "  n_IdD <- rbinom(Id, p_IdD)\n",
    "\n",
    "  n_RS <- rbinom(R, p_RS)\n",
    "\n",
    "  n_import_E <- rpois(epsilon)\n",
    "\n",
    "  ## Total population size, and number of infecteds\n",
    "  I <- Ir + Id\n",
    "  N <- S + E + I + R + D\n",
    "\n",
    "  ## Initial states\n",
    "  initial(S) <- S_ini\n",
    "  initial(E) <- E_ini\n",
    "  initial(Id) <- 0\n",
    "  initial(Ir) <- 0\n",
    "  initial(R) <- 0\n",
    "  initial(D) <- 0\n",
    "\n",
    "  ## User defined parameters - default in parentheses:\n",
    "  S_ini <- user(1000) # susceptibles\n",
    "  E_ini <- user(1) # infected\n",
    "  beta <- user(0.3) # infection rate\n",
    "  delta <- user(0.2) # inverse incubation period\n",
    "  gamma_R <- user(0.2) # recovery rate\n",
    "  gamma_D <- user(0.4) # death rate\n",
    "  mu <- user(0.05) # CFR\n",
    "  omega <- user(0.1) # rate of waning immunity\n",
    "  epsilon <- user(0.05) # import case rate\n",
    "    \n",
    "  }, verbose = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the model\n",
    "The model is first parsed and compiled using `odin::odin`, and user-provided parameters are passed using the resulting model generator (the object `seird_generator`). All parameters specified in the model description above as `user()` can be set through the model generator. Here, we use 1 initial exposed individual, 1 million susceptibles, and run 100 independent simulations: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAGFBMVEUAAABmAACMjNmZmWbS\neabmcwD/TU3///9Z1nUoAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2di6KjKLNG\nY7YO7//G05FbgVwKARX51jl/t0kIG3e7poqCmI8AAFTzuXsAALwBiARAAyASAA2ASAA0ACIB\n0ACIBEADIBIADYBIADQAIgHQAIgEQAMgEgANuFGkDwBjQ6/m2zxCMASDA5EAaEBXkWS88+Je\ndhgAjEd/kQ4ZZG4YAIxHT5GUR+aQOQwAxgMiAdAAiARAA/qK9FG9ZidJEAmMTediQ2i5KjcM\nAMaj8zoS0yOIBAYHC7IANAAiAdCAa1K7omEAMB4oNgDQgP7lb/eIMQwAxqP7guzhMDsMAMYD\nIgHQAKR2ABQQm/Cj2AAAm/jljPI3AFw+5qOqx5cixxcDkcAAaIEgEgAVJJKrC0RiZHYQCQxB\nfKICkQAoIKYSRAKz8eURff+NVTusI4EXcUPVLiVS7IavLsz/fExA+38dcJJ7qna9U7vGVyv7\nks5c3QVNC86yrhNQx73l70vmSKd1Ob7T684+8p7VPzc3LGZT7lnW9gHOc++C7FXFBoY/vizC\nsUroP7zuhA0vAZVE5touaMo9y+pOwElu2yLEpc2PrlJJkAded/aRIx35uZlh0RG2OMvqPsBJ\nojP6V4mUM0nQoKLaO0FJuCpZvUzn9llRYlLToASVHkh3kXZ/r9v9zTZJaFf0dRkySTgm6bSq\n3KSSpsyzrO8EtKS3SJ+LRaozSZw2ScCkuem+jqTvWcwfRh1pkc7HJNu5eEZMqu8DNOR1IgXK\n2TGTtB1CmWRelU+Z7jyTyLM0DmVNCh5WnGaDTkAz3idSPibZNqq9iULGIccko416ZN9UokdB\nHsg8zQadgFZcI9K1e+0Ci6x+SAqZpHQiAUrol21bQVpgmgQMVxQbrr5nwyGlE1QiG5KsKkK7\nIxjJ3dmQhOTuzVxR/r782ygyIckmd8d6Q2FyJ0r0aJzcISQ9iXctyCpC6VwoJH3jIenLC0mi\nKM4gJL2XV4rUKiTpC1WHJPvIPCsQksCPl4rEDUkiGpKYsySBkDQT8Q/QvVMkYXdyh4zSApSF\nJKeMZ5+lTmRFQuFubO4RifcR2C4/OhmSvuGQFBCJZHMISSC1HvpSkZIhSaTWkqwe1BEvJCmR\nTBpIfmxuWOymrLOESNdy033tOB/p6/SjWeUGRm5HK+BOPZw8W1JuiByfPk2YdCkTinSu3CAy\nuZ15hNxuSuYTSYekoFEktxOBkMTJ7UyL0nJD8PD8WUKkS4nPVF5ateuT233pI/OsQG43FhuP\n0FtnFMk3iZHbWYGQ24EQt95En0OHH83I7UhwqsvtINIkTCxSOrf72jATyO0OIvm53ffEmmzk\n+OxZQqQrmVKkQCRKTZJYa7INJkkISQMzo0gHk8omSUQr3R02N0zPdHvtfkRzO2+SFCqAx3I7\n4YUkFMDnAiJFJ0lBkZxJEi2AH0W6O7eDSM/gzSLFJTL6FE2SzEKuqMntmosEk57Ai0UqmCSp\nkCTss86USXeHSRKI8WqREpOkL2eSlBDJnySRyxkizciUIhVOknR3zElS7sr+Ro7PniVEegLT\ni5SeJFFHgiLdv5KESdIjeLVIcYm4K0letYG2xSQJEKYSya826EmSvqDVY2H0yFQbSEuINDlv\nFumQ2wnnwMyJ7CRJeLmdk9rxJ0nZYbGbss4SIj2AqUQKTpJIvnb9kixEegsTi3ScJIVEQrUB\ncHi3SAmHQtUGWl9AtQGUMJdIsb0N+j/qJAiVVRu+EGkKpvwYhbrWYyLR1C5Qbfii2gAOTCpS\nrtrwNSbFJkknRUK14aVApIRIIiFS92pDCwm+MOkqphUpJRGvbCd8kZ5XbYBIl7GL9O+PqT7Y\nJ4rKdmaSROZIx7KdoG2fIhJyu8uYV6REtcEp20WrDamy3deK9C2a+UCkQVEiBS7XyUSKpXYn\nNwk5wauk2hA8rDlNwGflEXqrjkiBVyLHF9PrR7cq26WqDQGRsEnopUCksyKhbAcI84qUBmU7\nUARECkajZNkuJRLKdpMysUgVZbuQSCjbTc3sIgUiE6dsRxM90h8t26knD2W73LiChzWnCa5g\nWpGale10d63q3yjbvQyIpGNPsmyXrH+jbAdmF+lc2a5epNZluyZxDVTwepFS0eisSI+rfyO3\nux2IpEVyqg1CGGFQ/wYM3i9Sov6dKdvZyBStfxORUP+emklESta/RUQkRtnOFQn173l5u0jl\n9W9RJlKk/l1StmtxlhDpZiCSbjN4/bu+D1DD7CKdKts1EAn175fxfpFS0SgmkkiL9Lz69xch\n6WYgkm7Drn/HynbiRpGQ210D3WbnbrmbRqTIQUIkERYpWf8mPzc3Ln5T3mk26ATkmFokr/5N\nRbLlBH3olO2+TqJH+oNIkzKxSCK1kGTr18YJV6RM/ZssREGkGYBIsWlSYCHpyxfJZH9YSJoD\nJc/nc7i33dwiFSwkta5/f9uGpC9MugIpz0fhvBI5vphbRErUvzMLSQ+sf2Mh6Qq0SOJwm8gJ\nREpFo9eIhNyOz388Qm+lHzWHSDGRsJAEMswtkjDaHA5I/dsr2+m3ioNIqH/PC0Q6u5C0l+2w\nkAQkM4uUXUja21QsJJnpFER6PTMXG8oXkgRfJJr9YSHp/cxc/m6ykCQyImEhaQ5mXpDNiSQY\nIh3q384j8yZRIhIWkl7FbCIdpKL1byKSSIuEhSTgMoNISY4i0WwusZBkHkEkMJVI4YPTIgXq\n3wIiTcsUInEXkoIifQ8iYSEJHJhApMNCUkAkNe3x6t/7E75IWJEFAeYRKTQ/arIi6ywkQaRJ\nmVqkL1mRDde/syK1WUiCSMMzu0gFK7LuQlJApMJqQ+sVWSwk3QlEMqkZVmTBeaYQyQtAQZG+\nrkjp+rdt+xyRkNvdynQixSLSl8ohMiJhRXZOPqHtquqlyPHFXCuSqBQpWv8WEOnVfOImzSGS\niIjEWZE9LCQlREL9+91ogQImTSRSIMXDiiwowfgzp0jxsh1WZEEJECktUqL+nRUJK7ITMXlq\nlxTpRSuyEKk3kxcbciIJhkiHhSTnkXmTKBGpcW6HFVkmfzwC75y8/F28kGTqDAmRsJA0ITa1\nO74UOb6YZ4mUXUiy+kCkmYBIoZSuWqRA/VtApDczfdWOuSKLrQ0gCUQ6JZKeKX09kRIrsgIi\nvZnZRcouJO1ttDQiJFJyRRYiTQJE0gEnIpJK5Y4isbY2+AU+8mNzw2I3ZZ0lRLqPqUSKRqSv\nbkPr33SSlBAJWxuAgEjapoxIzK0N50VqYRJWZG8EIkVXZE9sbSgUqXFuB5FuZBKRkpwTKbQi\ne69IyO1upL9Ikb1J8WF0oFSk7EKS1QdlO/Cjp0j6i80YJl0tErY2gLb0FukT+prA9DA6ULe1\nQf6/IxIWkoBLf5HMIXcYHeCIROoOJkRFRMKKLPCZQ6TsQtLeJijSFyIBBnOJlNnaEBDpzNYG\nvh5YkX0NfUVSdYZ8teEakQJhyd3aYJZiSP37y93a8MXWhom5SqSCYfQgJVLB1oakSNja8Hpm\n/6h5eo5UsCLrihRYSLp5awNCUl9mv/mJaLW1wREJWxum477bcfGWY58nUmRFliRz2NowIXd9\nHumTCIXxYfQgnNIViyQ8kVD/ngt7HR+u6M7FhsNRrGnjH+0TX5H9YmsD4HKbSKHDcNPGP9rn\n5NYGEREJK7JzApGSZTuINBcLj8A7kdppkaJbG/Y2IZF4WxuOS07CdJMeFrsp6ywhUlduEulB\nxYZsRFKTJLOmeah/J0RK1L8h0qu4S6TnlL8Fa2sDqTY02tpQIFILk7C1oS/T344rHZFs/Vuk\nRBJfX6TAQhJEejOzfz+SEE/d2vDlN2WeZoNOQIx7twgxMrvnieSW7axIVp0mWxsg0lDcu2n1\nmSK12dpgHqH+PTnziFSztUEaApFAlGuqdg9YR2JsbSDByRNJ+CI5hQmIBO4U6UNp/6M9sLUB\n9GSW1E6L1GlrgyfSqfo3RBqZ2USKRiSV25mlGP1YP8vb2mCimP2xuWGxm7LOEiLdxGQi9d7a\nUCVSAwm+MOkmZlmQTUYk5tYGspehh0gtQhK2NtzFPCIl0UWCZEQSvkjfg0jHPUKZCxsrsi+h\nu0gfcptI3jD6UCqSaCNSzg6I9BJ6i/R5mUjuimxAJNS/56T7OpLs9vYF2VNbG4hI36BIWJEF\nCoh02NoQEumwImveZzpXf0OkOZlXJBKMKvYImcoERJqba0S6f6+dqX+ntjaERTqztYH+2Nyw\n2E05QKSbuKLY8IR7NuRF0tWD0NaGvEikpVf/zg6L3ZR1lhDpHq4of9//Zcw/GFsbXJGOZbuo\nSDokQaRpmWZBVoQcIiJpEaq2NtSJ1MIkiHQTE4kUwxHp23xrQ3ZY5LjFaTboA5TTVyTm9/Xd\nK9KXiFS0tSG8Inv31oYGfYByZrlBZDOR3BVZbG0Aku63LOZ99+UdIgnvkW7DXZG1+kAk0Fck\n0+vTRRInRMLWBmCBSPkV2eTWBvtA/Q2RpmSe1E74/ny5ImFrA8gyT7FBLyTlRBLfUyJ9sbVh\nauYpfwtsbQD9mGdBNieSFiG5tUF0FamFSRDpHmYSKQZ/a0NWpNDWhsyV3XhFFiLdA0TSNnFW\nZD2RBEeknB0Q6RVApFKRnBXZBiI1niRBpHuYWSQREOkbEymwImv1Qf17eiCSPg6LJBIiYUUW\nGCCSEcKKZBeSEhEJWxsAYSKRsLUB9GMqkbC1AfRiPpGSK7KuSNn6N7Y2AMV8IgXA1gZQC0Ty\nRCILSa22NuSG1TYkQaRbgEjGJu6KLFUHWxuABCJ5In2zIonGIjWeJH1bdAJKmUqkJDmRAiuy\nTvRC/XtqIJKaJDkLSUQkERUJK7LAMqNI5SuySZFMZQIiTcxcIvFWZI8imRXZApFO1b8h0qhM\nKFJ4RVZ8dbnb39qg/8yJhBXZmZlQpPAcKbgiK7yyXUIkrMjOzUwixevfjkjfhEjY2gDCQCRf\nJFEnElZk5wQiWZuOIn09kUIrsoH6980rshDpBuYSKUkgtTuKRCdJWJEFBohk07uDSMIRKVBt\nsPpApMmZUqSShSRsbQAcJhPJW5EVYZG8hST11uxCkisSFpKmYkaRxLdMJFq2S4mEFdmJgUj2\ngFbdykVqs5AEkQZlKpFyC0lf06ZiIQkiTcnEIh2k0td/nUh1K7IQaVAmFukYljgLSbSO0GFF\nFiINymQiJcmJdJgkYUUWaCCSyeuchSQtR1Ik+TJEArOKlFtIMvmaVSe8kGQTOYg0N7OJ1HhF\nNioSynaTMaVIuYUkMkkSJrf7MlZkIdK0TClScIpEREqtyCbr3xBpWuYSKbsi+9VtYvXvtEix\n+nfuysZC0vBAJGqTyIqkFKH1b6caTp7FiuxMzCZSEpPKxUTyJ0lRkVC2mw2I5CZ4wohk69/J\nFVksJIEfs4vEXkgK1r+1SIH6t4BIUzGdSImFJBKRvIUkPS8qEAllu7mYUyQsJIHGzClScHpk\nZkZYSALFTCZSrv6NhSRwjulESpKtfwvPpFYLSQXO8c6yvg9QBEQqFumwkIT6N4BIcZG8+rcW\nBgtJIMD0IjllOysSrX9/TQ0htLUhUf8umPlApMGZT6R8/Vs/lm8QJfVvLCTNyqQihdI6a8/Z\n+jdZiIJIkwGRfJtEptqQEIkuyaL+PReziVRftlPNcmW7wvr3l+8c8yzrOwEFQKS4SGS+k6p/\nP7FsB5GuBiLFRDqU7fRT2P8NjkwokudSbf07VbZD/XsaZhWppv7NKduh/j0Zs4oUyezSH6Rg\n1b+NcoVlu+BhzWmCK5lQpPQkiVX//kZFek7Zrr4PUMB0ItXXv7VOurtI/buw2vBtG5K+MOla\nIFJUpO+h/p0q25lHKNtNCURy0rozZbto/fsLkeZhRpEqynZEK179m/zc3Lj4TZmnCS5kWpFi\nQSlQbRAl9W+32lBQQihYc2KeJriQ+UTKTZLYZTun2vDAsh1MuhKI5IvELdsJR6THVRu+TeIa\n4AKR4iKZmOOkdrFqg3n0EJEQkS4FInlzpqhIqDaABFOKFAtKtmynSwjyDaKsbOcGL/pz0+MK\nHtacJriOCUUSqbJdptrwpdUG3V2HTUItzhIiXcmUIiU5U21oMkn68puyzhLVhiuBSHGRbIXb\nEQnVBnBkYpGsUJFqg7kWzSRJJnaoNoADc4pkPPGmS7lqg/aMIZJ+CiLNwYwixasNJKNLTpJU\n5CLZXONJEkQajTlFYk6SrByXVBuwt2FgZhYpIpQjksnt6Gf6DiJFJ0niRpEQkq5kUpESkyRh\nRaiaJGkPIdIUTClS+SRJGHfMkqx+XhjnyITKPotJ0hxApKNUmCSBYiYV6RCFGCKlJknCDUnP\nEAkh6UJmFSkwP6K5nWxk5ktObsedJOkDiDQDc4qUmyQdQlJgJUmQQoJ20j4yz4pbJ0kw6Spm\nFelMbmdEuiK3a+IARLoMiJTaJURCkhDCFQm5HSD0F+nzo2QYl5CZJGlxQgXwaG5nrtsvMepb\nlLChAD4qpSL9nGCIYdpKj7JvuF6kNrkdEemJBXCYdBmFIv2UYIUY2VhLl33DDSKx6nZ+SFLP\nZnM7vV6L3G4WikWy/2M3Fvn210/PeHU7Z75Dczv9ihuSArmdZ9LlIsGka5hXJJrVpXI7I8ch\nt3NDkivS2dzu29ykBp2APNOKlCs36NVXLyTJI73f7ktEoktJ9oC+Tz+TGRa7KessIdI1lBcb\n+OUGU2fIT6ruEOkgTyAkaWn2NwRDEglCrJCUu7IRksakuPz9YQckR6SCYVxETCH5lCAhicjB\nLzcQkb6iJM5ApCGZdEH2x9EgzyZhShLqDflZkrOEZJ89WW5oYRK2gF8DRApLZEOSk69pd762\ngJ6pgNeFpBYmQaQrKBDp48Dsntn4FoczIclMexIhiUYb4qBASJqNriLxmz9MJGGdoCK5JtmQ\nZIIQCUlkVdYU72wv6XFFjs+fZn0nIEPP1M7688SqneDVG9Qj/QZdq5MR56u10t159QaSBhaY\n1Lhw94VJF9BXpNBhdhgXEjcpltyZuZM26ZDcmcbxyl22BB45PnuWSO76c1Ik9oIss/29IkVN\n+vom2ZAUSu7CJn2LF5OclxuEE5jUnzMLstw50uNTO5GNSQGTvtQkbUvMpK81yY1J2VHRR/Un\nidpdb8p3fxcU7R5ebNjhm2Tau08mYtKXxKSvezVns7umRXCY1J3yvXbycxHc7p9c/paEDXJs\nOqikHDMSucU64b5PiKOOpSY1UKlBJyBOqUi6WeMr/0aRGCaR9I4EGNcnc6F+7cvCmCRCJuUm\nSq5KlRp83Z8OGgORBEulmDrei7HubAbIV8m77qs1cMYIGtPzYxTkPSXDuAWOS0zPwr0ZiXRr\n2y43JvpEi7MUTYrqwOHsR83fJlK9StGCReZNyemS21A0qeB5PYImnLz5ScmFP4hIP+pEcgU5\n+bbjiILjsuM9Eae8wdX/3kDv3d/c+vdDRGpp0nliYwsKVdRF+mzlU41+k7Nxo0gnNpNfRnM5\nquAOTr8ugk/nzle47+C9ESjOLsi+b44UoYcYLThkZ7x3xddmWe9GmSIGROJz4nJ/PO7JRTU5\n1eVUy8DnUru3FhtKoNeIugrBAZGSkPsb837riVb3Wdtz9/fJYbyH2JoSGJX4vzVEuhNSH5AP\n3X+yuy4XECP+T9ldJN7tuyYV6Szhf9Xv9yu83Up+A/tm/wq54jIcn/g/Se9iwwcigRno+nkk\ntRFCPPYTsgA0ovuCLEQCMwCRAGhAgUjl97UzIg2z1w6Ac3QWSTCbQyQwNmc+RlG0jMTTDiKB\nsTnxCVn6d49hADAeEAmABnQVSeeBSO3A2zlxzwb7V761YE6SIBIYm/JiA79o97ErSCh/g3dT\nviBbtogEkcAU9NzZAJHANPQV6SP0DiGIBF5N0c6G0o9R8DdCQCQwNl1F4tcmIBIYm867v7lA\nJDA2EAmABpxbkO06DADG49RHzfsOA4DxOLcg2/y6h0hgbM7NkZqrBJHA2EAkABqA1A6ABqDY\nAEADUP4GoAFYkAWgARAJgAZ03rR6ZhgAjAdEAqABSO0AaABEAqABRandiXt/nxgGAOMBkQBo\nQPl97eyfnYYBwHh0vWXxuWEAMB4QCYAGQCQAGnBu0yqKDQA4lBcbunySAiKBsel5E/2TwwBg\nPLCzAYAGQCQAGnB2joRiAwCEU/dsgEgAuJxcR+o5DADGAyIB0ACIBEADzu3+7joMAMYDxQYA\nGgCRAGgAFmQBaABEAqABuK8dAA2ASAA0AKkdAA2ASAA0ACIB0ACIBEADIBIADYBIADSguPzd\nfxgAjAdEAqAB+DYKABpQMkeCSABEwCdkAWgARAKgASh/A9AA3PsbgAbg2ygAaAC+HwmABpws\nNkAkACgQCYAGQCQAGoA5EgANQNUOgAZgHQmABmBnAwANgEgANAAiAdAAiARAAyASAA2ASAA0\nAB/sA6ABEAmABhSndn0ueYgExqZ8rx3uIgTAAYgEQANQtQOgARAJgAaUFxu63AIcIoGxOfN5\nJIgEgMeJdaSPwEfNAXCBSAA0ACIB0IBzImEdCQAH3PwEgAacKX93uPsJRAJjgwVZABoAkQBo\nAO5rB0ADUGwAoAG49zcADej/bRSsAAaRwNj0FOkj90BwTIJIYGx6i/TRNvGHAcB4lBcb7J+c\ntqolRALvpkCkjwOnb4gEZgEiAdCAnjsbjHJ58SASGJurRCoYBgDjcXZnAxZkASDgBpEANODk\nOhK7e6Z1EAmMTVeR+EU+iATGpqdI1h9U7cDLObezgdn3J3SYHQYA49GzageRwDT0rNohtQPT\ngGIDAA1A+RuABnQW6cwwAOjJUkG8155VuwIgEghTc9m3Jz7OC7YIcdpCpJm4W4fzxM8JIoEe\n3H3FdyJ+wj0/RqG7hUiv5u6LO0fZKM/+FjqLxI1gEGkYel7zJzgxxIpTjbftnNql2pd+ch3c\nRBcBTnBieM1ONftOzJFAiOYWnKF8cK3OtLybc6ld4SbwsmGA2+gkRBGlQ2t0oud/aZKTcyRE\npDfRUQsuZSNrc55nf11BLhCpdBjgKtrrUNq+ZGBtzvN8L2kg0oSc9SQhRGn73MjsYZvzrOiF\nR/diw4fcJpI3DNCNWmECQhS/ITOwFgKJ/A9szlmR2P1DpEdQq8xRCPcx4w3pgdnDJqdZ08kp\nui/Iym6xIHsbjdSxV6j3ONYue0UfemxzmlW9nAcivZcTlqRe84NQvKTGVyjbtOQ863qpBCK9\nknKHUohD4SvWTOSuaK9l7cX/BIUkBSI5W3rYW4T2brHX7jrKNUk7lHpIHcpe0l0cquykHZ1F\n4n4LOkRqQrFEfIeivYsih7Ihi3+alZ005lRqV7LVjqcdRKql3JOkQ/4aTKxZXgzPoTdK9OOM\nSB02a0OkGlo6dJQo2qxYoibnWdtJJ8pF6vKZB4h0lpYSuf1Fe7c/NzswXlP2edZ20pHym5+U\ne4TbcfXieonMK9lxkRE2Oc3aTjpTvrOhoG9+cQIiFdNNonQ6Jy6WKG/tMygSqTAaQaRe1EoU\n1ybUu72Wsxe1a1GT86zu5BpKRCrO6phfac740cBSKFGDUGR/bG5c/lHtaVb3chl915FkM4jU\njkqL4sEn1DlCEZ/uOxs4e78hEo9CiepXiuyPzY6L2ZR7mtWdXE3v+9r9TIJILbjAIn87kH4l\nNy7BbMo9zepObqD/DSJZ4QsipamzSNC/3e7so7MW2aZtTrO+l1voLxKrRgGREtRZZF+ARR25\nQCQOEClGG4uWUosKl4uanGZ1J3cCkR5NoUWJlM7tjjwSi9+qxKImwWh4iwREejLnLRKeRflg\ndNqiKfb/MIBID6XEorMpnXENwagaiPRICi2KpXQJjUQkGDE3LzSz6B0aQaQnUmJRap+2vkqd\ntmmNsuNiNuWeZnUnjwEiPY1SjTxzzAsHycwjYlGZRsyW7NOs7+U5QKRnwZHHM4PUDWpyuuyw\nmE3ZZ1ndy6OASE+CoZFvRuA9R41s50/RqLqTpwGRngPHnqBGzjt/z7vd2Ueo03UDIj2FRhot\nUY3EWY0Ch1VnWd3JE4FIz+CERsQc+wpTI/Jjs8NiNuWeZXUnzwQiPYECi+IaCZ5GCzTqAUS6\nH4ZG+aROUI1IWczTqCwacZtyz7K6k+cCke6mViOrinBf1p0LRyPMjfoAke4lrxGjxCASGlUX\nvAU04gCRbuW0Ru5Sq33ZthXkRU+jnBvQqByIdCN5hxyNhKPRMV9rpFFJU+5Z1vfydCDSbRRo\npOINmRXpR45GqDHcB0S6iYxG2hGiUW5yZNt64cgJLPlwxG7KPMvqToYAIt0Dw6OgRvQl0Xxy\n1CGrq+9kDCDSHeQdCk+OvKwuGI6SGiGr6wVEuh6GRnUrR65H5OfmxsVuyjzL6k7GASJdTsYj\nK0ykVucWGci79EPxkKxuJo8g0tWUhyMahxjhyP6Qkhiz8Jtyz7K6k6GASJeSj0aMpSNxnBz5\nsyMdvOzPzY6L25R5mvWdjAVEuhJGODqqkQxHjnNaOH1Afm5uWNym3LOs7mQ0INKF1IYj0Ssc\nhY/Pn2V1H+MBkS6DEY4Ki3XHcHT/7GhSjyDSZTDCkU3LBAlH2jLHEDccHYt15MfmhhU+rjjL\n6k5GBCJdQ204iqV19pF50407gmYNRwIiXUSBR6q9nhI54SjukX22JFcraMo9y+pOBgUiXUBG\no8UPR6gyjAdE6k9BOApXGZhp3YK07j4gUndK45FTpmOmdaQj8nNzw+I2ZZ5ldR8jA5F6UxiO\nitI64hHSunuBSH3JxaJ0Whf3yD4yzyKtuxOI1JVyj/79z9o0yCLsskzvEUTqCtsj4U6PaFpH\nDDl6JIJpXd4jdlPmWVb3MTwQqSMZj2zgYU2P3LaYHj0LiNQPbjg6eKRfdjw6hqP6zQwC06NW\nQKRekHBzyiOsHg0FROpEOhxFygyLfdZ4ZLoz9jyozACPNBCpDxmPrDHChCPrkS0zEI+8eGSf\nLSsztPaovo93AJG6cMojUeaRfqpgzoMyQzcgUrAZqSIAACAASURBVA8YHnlqpMp11D6RLnvn\nhsVuyjzL+k7eAkTqQNYhm5Wp9iyP7CMBj54GRGpOPhY18wjluscAkVrTvuxNynVuXgePngNE\naowRgu2RuMKjBR71BSK1JRWMOB4d8zrikT6CR88DIjUlZE6zePQkj+r7eBkQqSVn45HwPTLd\nNfIofHz+NOv7eBsQqSHlHi2nPFrg0eOASO241iPyc3PDCh+fP836Pt4HRGpGed37XF4Hj54I\nRGpFuUeDxiOIFAIiNeIqjw51hty4nCE2OE14FAQiteE2j/gBCR71BCI14axHndeP4NFlQKQW\n1HrUaT8DPLoOiNSAxh6JhEfnbxfU4jThURSI1ICURY5H6c+Vw6OBgUj1MDzSAUi153sk9BE8\nejYQqRq2R0vSI+F4RKZLi/0h1Ah49CggUi0Zjw4TpKxHVj9pFXn2Po+aFCveDESqJBSC/IBU\n5pEIxiPhenTx/RngUQ6IVEepR84GO3FYiD0WGtTfNR41Sezq+3g3EKmKVFbH8aig8O2uCuXG\nRYfY4jTr+3g5EKkGPwAVeZRbQDJJX8AjLMQ+DYhUA9ujc4VvUe8RCnYXAZEqSKZ1tGBnSgms\nQoN+BI8GAiKdJxSCIomdeoPxSAQKDUYu07n825QlzI/NDSt8fP406/t4PxDpNOUTpIXcKL/0\nkxPk52aHFTyuOE2QByKdRWdmRR4tnkexCdJ5j/BJvnuASGdJRSOdoalD1d7UvOOFBrq/QZyb\nIEWOz59lfSczAJFOcsojY1PUI5IE6mdRaBgBiHSOkDrkwaFgRwrfdKIUKzTEdthhZ9BTgUin\nKC800MI3o2D3lAlSfR+TAJHOcKpgF/KotNDAD0goNFwLRDpDyiJngmQ9IiIZq/yPINlH9FlM\nkEYAIp0gpI77IOCR/jPu0aHQIFyP8NWWDwYildNsgmRSOW+CRJ49O0Fqc5qADUQqhuORX7HT\nISg6QeIUGvgrSJggXQ5EKiZlUcQjU6NbbHSKTZB0QBK0kUCh4eFApFJC6gQ9Ok6QUoWG4ATp\nvhUkeFRIf5E+P0qG8WwOiZ2gB7ZAZ0UiEyT9t0isIJFnMUEah54i7f58PhyThhHpzATJzImW\nyARpsY/cCVJBYucMscFpQqQyeov00Tbxh/FoUhYxJkhFK0iYIA1Ef5HMIXcYTyakji8STez0\n40Plm3Z3rNjdPEFCYlcORCpBR5RgLHI80gHJTdfKKt8kyuTGFTysOU1QCEQqITNB0jdZyCV2\nsQkSbtIwLn1FUnWGfLVhDJFSFjkBKXazk1xiJ/9GYjcgV4lUMIzHEq/Y+YkdWRXSHi2HCVIm\nsSM/NzcuZ4gtThMUgwVZPox4lErskhOkmsQucnz+LOs7mQ+IxOYQkJyHtkIXSOw4H4o1id2C\nxG5AOovEW44dQqQziZ0oT+xKJ0io2D2DriJ9LAXDeCgtErtYxc68qyqxazNBgkin6FxsOBwx\nhvFMShO7Q+Ubid2r6b7X7nCYHcYjKU/shE3sBDuxw17VUYFIPCoSuyWa2FmRYoldQUBCYncr\nSO1YXFWxQ0AaFRQbOFQldsV77Jyfmx5X8LDmNMFJUP7mUJHYCbvQSgJSJLEjesmXcsMKHdac\nZX0ns4IFWQY28IQSOycgqTcgsZsNiJQnntjFApI4VhrEAIkdRDrPBSIxMruni5T0yH54wuZq\nJBDpBI+X2J1bQmqiADyqASJliSd2TqXBKOEEpOMSUiyxwxLS0ECkLKzEzgSfQKVBxS0vIKkH\n7hKS4EcZVBoexTVVu5HXkULq2AOyhGRii1VHRxk3sXMCUoPEDgHpAdwo0ufDX2e6D06lQR/K\nN4hjQBKRSoPO7GxH9OemxxU8rDlNUANSuwwpi4IVu4olJFQaxgUipWFWGpyARMNM2RISsSM3\nruBhzWmCKiBSmuIlpEClQTATu5OVhhZnCZFqwYJskpA6TkAibWR7WmkQfqVBeJUGLyCRn5sb\nFrsp8zQbdDI3ECnFIbELBCQ/sQsGJNqfE54EKg3voLtIH3KbSN4wHkQqrXMCUjixO5a+bT6I\n0ve76C3SZ2SR4gEpXmkwk6JcQNLpYLlHCEgPpPs6kux2zAXZTDzS2+hypW+n0rB4iZ1R67bS\n9wKRGgCR4qRyOjcgLSYgKcN0pQEBaRYgUpxMQMqXvoVRRb0crzSUlL7pCJucZX0n4BqRhtxr\nF45DwYCk3mACkmAFJPUkKg1v4Ipiw6D3bPDsiSd2x4C00BmS6S4wQ7Id2R+bG1bwsOIsG3QC\nLil/j/llzAeHnIfCimSK24fSNw1I6i3HSoOzYCuydiAgPRMsyMZglL5jM6RoQCKJnVOyuzMg\nQaQmQKQI8RnSYgSJeBTZ9Y2A9Gb6imSSuuGKDQs7INnSt141Kih9L8Yl3S43Ln5T5mmCJvS/\nQaQ64g/jERwFKg5IguRsS+xjSMvpgNTEAYjUit73/lYmjSZSMBZ5AclaJRCQpqf7TfQ/nxEX\nZE8FJHEISMILSPaReVYgIL2C7iLxtn8/TKR46fsQkNQbhIlFrICk/kZAeg39RdpNGk6kZEBK\nlb6LAtL5kl0LBxaI1I7OxQb112BbhA4lO3EUqSAgLZGAJBCQ3sMVIonxRKoOSN4akvWoUUBq\nc5YNegESLMgeqQtISy4gCQSkFwKRjrACEpn1lASkxQtICwLSO4BIB+IBSa8UCWEfCRqQFgSk\nWYFIByoC0lISkCrWkNqcZYNegAYi+aTXkEz8cQLSkgpIAgFpAiCST3RzUMEMSfajXz6W7JYH\nBKQGnQALRPLgBqSFO0NCQJoCiOSBgATOAJFc0gFpCQakghnSgoD0UiCSCysgGaWaBaTchY2A\n9HQgkkNpQDJVbzOBcgKSQECaBYjk0CIgLcIJQZghTQFEoiQDkjUlsqkhEJD8bd8690NAehsQ\niZIJSDazq9j2bTuiPzc9LASkxwORCNmApEt0CEjAAyIRWAEpvYZEY00wIImqgNRCpAUidQAi\nEaIi6TikA5IWSRiRygMSP11DQBoAiGRJhSPmDMkLSGKhjxCQXgxEsvifQ8oGJCuS0UzEA5Jp\ngYD0QiCSoTQg0VwusTvIPHpGQIJIfYBIBmZACmV2JNyEA5LwAhLdHZQdFrsp6yzhURfmFGnb\n/v0/5d/DZdkC/uiDf69v22Y9+j0wHm2/t+6Ptv1q33tXIm17Y25A2jZ7qF6uEGk/Mw+I1Icp\nRdpCSI88mcxTP8+USLtSu3jbop+Vmm3KN/Wyfpd6JLStWyQgWa2FPTrvEemO9AGP+jChSEGN\njgFp80TyPFJPEH+UVfLK1SLpzsWi/v6nwx7MNAsdlx2feeq0SJuwnZA+IFIf5hMp7NEWikVO\nQBK/q3ALBaTf/+12LEqtTTf994owAenf02rStC3m0iYX9kbHaF63z5Y5QPXZznYC2ECkXEAS\nJrNTz8uphwpIQgekfb5EMzuh0zhBAtKiZ0n60g57ZC/9xT5/2iPyCB71YjqR5HW++kTCkeLf\n67rNuuzv3p/758m/P38q/Z5Z96fWvfdfu581uvP197xU6/dY6ChhL2x1qe/t90nS70DGL/lK\naUAy3dHHEKkXs4kU9chTaSUH6/6HkU2+QT7Wf6i/14X2t+q26l2yN7HrpS5tzyN52e9X/z+d\n/v2xrOa1Ex7Z89M/Dh71YkqRDh6tYYmCAUmLYtSRz+k/1cubcnBV4el3sD8jC3G/a9upyRGP\ntv1IXvoh53inabrTB3Lm1ur3CDxmFOno0XqUx0q1mYC0WpH+iaHDjAlIYksHpP2nCFm3kKqY\nC3vT9lCl9pKdlqDCI6HFhUf9mEykmEcFMyQtkn7FBqRfqUEHJBmUhA5Pyi7xS9UWPX0hNTnr\nkT78vUNV93YLikWiHskuIVJHIFIuIC10hvSTRjilBiWXCjqbfln+/yZUeNItZXVC6fM71OP6\nPTRH8vBXlDA1g8WpwmXP8tfF6j65G13QCShiLpF+e4GCASlWatgfbE6lYS/K/eoAe0VcFSL2\nZpuZIakudX1vN0o++0v/THHbHPoBSXivC+sU8zTFuh6e3uvzUKkPs4l0JiDp3E0VDTZaathU\nuJGzpl869ps82RnStql3CdXRr5EpRutDcuGTw70HVbLbZ0r8s/QTO9ndLiZM6sJkIv0CkqfS\npgOSCTpO7VvVF0yKt8mpy6H2vemS3eaKpFeT5LxJpmurqcNJk2QqpkZIAtI+0VFNTfWOcZbR\ngLTCpF5MJVIss0sWG0K1732Twk/BQ+17nwzp6OOWGja1brtpk35Tf3NphwKSfp0e8k4zEJD2\nTnYxYVIHJhNJByTr0yEgJRZjbe1byOlSKCDJ/gStfasAtgclGbT2wLWPSUUafd1vTkAS6tI3\nTVkS7CIdnl10JzCpB3OJFCk1FAekff+P0LMmWfteF1WyVruDzMrtIpdC98rbppLDX3VPl+xW\nGaP0ED2R9uBimnIys3hmp/orqVsAJjOJtP+XOrA9KKGSH5BkyW5fjP3tV3V2B5GNdXK76qpq\n3ntI0dU9GQ6FWwYXq/7Iw8Ej4ZTsfu9knebhWbuGpOddoCVzicSqfTspHglIi97xI58Quqqg\nM7tNzZC2VW3yXvXuoL3ksKqgpBJLUgbf7KUdFImUwRnxZMuJBJN6MJVI50sNOrSsZJvdvutb\nWbVteopkAtK2bzpV+aDZ9y03h69aPPlBdJneySE6pQZ1SPXJmpT1CCWHHswlUrz2be954tS+\nbUBSS0Wr2aWw6SmStkrHN53Zbb+ShNxTtK/q/nxYVm3Ssqd3+8fK920I+6UdCUhuEMnEE4ZI\nZs8saMZEIpHMjuhUXvs221V//ydM7XtTAUnOnPQjsaj933vBYO9CftZo0cU7vYikCgnHgCTU\nJtfNBqWUBByRUHJozmQiHXK7Q+3bfDLW2Wa3WpHsxrrVBCSd2NFtdnKX6GJKDsu+hXuVKpla\n+Wov/N2PcECSL7PSO55HSO9aM5NIOrNzF5GiwSgVkDZ337fcHWQSxU2XJfYPpqt8T23mVvao\n6t1Ga997jFKDDd3zhOgTD0pckWBSW+YRSa73szM7JyARkVYbg3RA+hWoN53L6d1BYt3UnVKE\n7Uh9vEFPqfSKkr7w9516+1Hw3kGsiVJQpPDnJ6BSQ2YSKZjZedDMToeckEh2V8NmJj+69r0t\ncor0e5vc4KOD0qZNWjZbB7d+bMaPyE24HJPCK0q7mT6RDyLBpHbMJFJgNbak9q0yOx2QaO3b\nfoBCry3pO6DITaxiMXvxpEqLnBDp9I4Uv1fvI+GuAvTKDwalgoAkYFJDJhIpkNnFSw3CRCF1\ncAxITu3bz+x0C7mJ1YokL1x9ZIOSKX7Lz4TbYYeqbYnqXZlIWJxtxjQikSkSJyAJG5BoYmc3\nBK0qs1O17816ST7ftydz+77vfaOP8vj38GCSKX7vDtpxB8ptImHSvgzsvyX1EXMEpTZMJFJo\niuQFJOe++aHa90pFIlb9ECSz03neqmvfYg9DcgCb2j9E0zu1OCvr4OQj4QEFSFA61sFDU6T0\nrRpgUhPmEUlokWhmlwpICw1IVCRtkyw1GJE2K5KcIS1qV8O2LuaOQLIbfZcfG5TkNge9GmtD\nTaTcFjNpKxcJ6V0TJhLpkNht6VKD0B8iNwFpXQIBaXNC0GI2eOvd4vIpeRsgEwhN+W412y22\nVV/Q+yqu/jR68GRoeueYtPfmNc7fOwhBqZ5ZRCqbIu03Fxbqg0eBUgPdZreRu9mRzM5GNFMV\nl2FI7moQxiRSbyB3jfRvIOlDgxKR4ExAkp0gKNUxj0ibv60hGY5kQNqcUoMXkGypQXenS3ab\nMDMku6thNbdtVJngMShJnchHWRO3hbRBiUpwViR8craWeUQiF6sTkAI+rToquSItfma3mo9G\nuCLp7UF6d5De1CpNWvQQzC0kzHZa75Oz+XKbaakWogIise8KiaBUwzQiHTO7gz/uIpJIByS/\nZqeV+c2HVmeGtKpNrXrhdSHRZ9MxycyUtnWxSVv6jo40vTMfbTovEoJSDZOIFJoiRZO6NZzZ\nRQOSNIBkdpvO7Oyt9JfV1LgXYePivuzj5HebCV3CKd+FOAQlsn1cU3SbYgSl00wj0nEVKbVf\nVcjMTj82U5tFLy0dApKb2ckIZzM7tftuL4fr0reKS5tn0rIZkxa/KHfAUWmryuz8/kARs4h0\n/HBsNCBJCeRX8qnimxbJC0heqWH3Tj8SpPa96o/x7Y6ZooPJ5nSI0p9klze+22j5LoZj0hoS\nqfA3hfzuHBOJFAxIyVKDnDj9RJL3Tz1kdkswIKml2VVmdPaed7IkoS99mt+ZQrjziT/zXRTM\nmZL6rJND+RdQrFDpDHOI9LtUDwHJc8j5aOy2BxT1vNhvl0UmRYlSg+18kzfO3/QMSQWMRX8s\n1g1Kxina1Ik10ZOzV76Mec6LZ77JJacuCDCLSIf9QanF2F0kYUT6ZXahgLRGApJWytyvS29q\nVRFpC6okH9iWm94AIUSuCqBVkvsanJanvhJphUrFTCJSYIqULDVsat+2/CjRLpK8uRbZZhcp\nNaiZ0brq3Q3y9qpKEL0YK0sDZEyqfmcmU07VITt1USpt8hvKyJ7Xk18ttiK/K2QekTyPkmyb\nn9n9RFp2p9yARHc16MxusyLJuzfIz5WbKZKJjN4NKzehp1CqKTUpU3UwUc1tefo7+hCUCplC\npNAUKZ3ZLUqkTYukY5H8mqFU7Ztus1tWuT1i097YooO0RH+JuY5Ki9k+pD8G6FXl0lFpVfvH\n9efQa77rEioVMYlI/irSodQgnAf7zjj1/KoyOyXSnvfFdzXozs2z+o5B6iO0JB5ZS8i47KoS\nacorOtBbHyvp6r40NmcuIMwh0jGz25YwgmR26/54/zbWRWZ2Qom0MkoNVi0dkNZVf5fzSkof\nVKU9EAo9RNKUo5LauUdUWmq/fRkqsZlTpENAWl2R1O4e9c3kYp/lCB2QxJbL7MwMSa/Ybvqj\nsUYhjb5dg+pnNfU7UqBYj6EmeJb7FIk2XZZqCYiaIMUMIh2nSJFwRDK71WR2YtVfpGwC0nbc\n1aDnRauf2a12smMO9S4G84KSZ7MzJN1UZ3oMldS2BnLpL0uDgAKVWMwhkh8H0qWGPZHb9J67\nVd63RGd2P83kV4yJTKlhlaUGnfDt32xOq3XOBzu2/ceYpps9clodP4HknqZ+WW4zX5rkZmuL\nTl7PFCIdMzsP55ndFLsYu0+RNiXSvpv1Fy8EESlRaiDzqFVPe2wt3JfEaWrLDmSy5H9uwjlL\nuz9IybiINtMcqJRnEpGYAUlvs1vV/Rw2ndmZL7vUIskvhpWVskypge6mCH7VmVZpMTU7tbJL\ntuDlVXJvxCX/c2FbtpgsQaUEE4hkp0ibuUojAWkPO0IlcntIWlVmt5qAtEmr9n3c++XLD0hW\nCu8P9bq+kYPZuup/xGKlscG7tP2d31pG2XRrMllCWIoyhUjef/0j4YhmdoJkdkJ+b6UMQSog\nyR1DQotEjAmVGkhNzuphR7XpgKQuVduU5nUBl6hKx88iLbpeYZo2yPAQliLMINJhirQe5kUm\nMilT9DKTzOw2+zXmQqiAtMpP+hiRSECyuxqOAcl4sHoRabOekVKDq5K7yfWHPTp8ylzPkNg7\nIxjIH13VxVt5v0iH4ncmIAmVyO17V+Vi7P4dLiYgyRvo68xOyIRwSwSkg0gm1Ljj0u6s9GaS\nygTyJtclfRQISMK+vgnatgbaH7DMIJL33/5sqUHuZNCLSDSz2wOSvC34/ukKE5D2zXnCr33b\nTa0HkVa13duasXcqX1i88COjgLbHV0lf24GARF82D5Dh9WACkbzM7hiDqEi7KaueLKnMTom0\nLDYgrUKJJDVTIolcqeHgkv+6XLfdVpP/UWlCKtkULyaSd+mvDSLK2qKTdzGFSF4GlczslkBm\nt+8DF8Jkdov8kvJ9hq9uy7rsC676AxSLmSolRaIu2di1LXpTkZ/KhYp4eorkXduLt83Od6lR\nilfXx5t4vUj+FCm+79sLSCaz22xmtyirbGZnApLQ9/uWm9xWP7MLe2Rs0V9CK5uShdvt2Pbo\n0iZoSeJH+EtjvRSvvvKAsKSZQCT3us0EJJLZ6X3fMqfbbEBa3cxOL+GqzO73Q9kByUgRnkxJ\nofzKgy+VfiBs4PEDkoKY1sQD2t/kvF+kYECyNjn3PLG3V12NSMJkdoeAtH9myRFJKI/od/Tl\nPJKvk6qDaxmVx1hz+FzIum5kf9AW9ki92jrFQ1gS7xfJZHZmVwMzIMnP2B0DkvzqWCuSqn2T\nzE4IVagQjIBkS4kmW/Obmk8ZuS55HThS7LtlY1e317R604N2e27eL5K+1mQ9LHc3O7EJPUPa\nN4GqJ6xIi75PnQhldvKr+fQtUxZz7WcDkjYhOZkS3uPAxj15Ret937HfCmnaRAS49H6R/Mwu\nhpvZyYC0rxztudweiGQBzw1I66JFWuW+CKEyO7Hqm5wWiJRrekjqQntgN7vvO3V1U5n0G2t+\n0/U9DM3LRTrU7GIq/ak//mxmtx0zuz8nIAkbkP79lL+9C/GnAtIfvVvwkhHJq+4d9bBB1Xcp\nFJZMYpfxw/WnXibq5mx0FukjKRpGS7z/gMcDkhLpT+qwK7Uf749/R3/yOdXm70+9bNv+qfer\nZqbRn3uocccl/+I0Vddq8rHa/ur2EvkNufrApbN0FeljKRhGS4oDkhFJu0FFsg+ISNqj/VD8\nWbWsEotrh5Hm3+GfufDl44W8qNqGbfJXmHyP5Jt0ZPpbczJt7Up5a30X49FTJOtP1qROInn7\n7AoCUsijY0AiwYsEpH+p3B8zIMlDdeHr2OS+LB/EQpP70Ipkm5OSQiI0OZGk3oT5AlNfkUKH\n2WE0xM3sgtU64pEKOcyA9Jsh+QHpLxKQxNEkc6HrgKXGuK6eL07TUGjyH/07Tb+dvKStS0GZ\naNOGgWkSmV4vEv0EamFA+jsXkMgMSdCAJII67fI4h74pq5vdBWzyZ0ih4OUU6DKRaW0r0xRp\n3qtTO38OvhppYjOkWEAympFQwwtIMs/7M4+OKunrXvUf8Chj0+b8FfGI2rQf5UOTf3wW6vFr\neXWxgRuQTGK3KY9MJHESu1BA8jxilezccpqpRPyGSGWJ2xSQyRTIl6SM5pJWh1GZjstMbULT\na216c/n7UBROJ3aL8uinhNxHqkTSQWqJByQasvT/Yh6pUBB4/d+7rCJRG/7Iy8eGvkdb4Ija\nZKrkx9+fbbq1CSu6v5o+nsqbF2RLA9KmA9K/q2rd11QDAYmGoENAOiR2h9K38B8todgVdiSm\nk+OdltEhfAuVdbXFmJxNdLdffT3vhTa9WCT3BqsJh9QfQnukAtLfnibtO+ic+p3yww1Ieg3K\nJnfHgGP0IetE8abk+s7YRE0IeqQvXu9vc00bn2I60abtpk2v0unFqZ0XkAJFBs8jKdK/lpv9\n/olVfqeDjE5LNrFbDoldJLPzKnrxpquby2V0WmgpPYO/O0LQkl5aJ//BORp08RzeW2xws5hM\nYrffNNV4tOnPt6rb5v882v45su9isx7ZMrmtNKQTuwN0iSnVlF7gKZv2gJRumHTLRhyGTtan\nBjaN7tN7y99OQMokdtSjX2ha//4p9CfvXic3pu7RafuTt17V2d4hsVsYiZ0rEr+pUiNpCfE2\nH8hSPgnfJ18op20DE4b36bULsu4denIe6QnSsu/03h36W1RUWpRH6qZA//qW6ujwJPspTOwO\nr+eahn1yrn85Yv1ANaBvOu2T/dHeb5m2tm87q4L304cy6q0iOR7FpkeK/X51us26xyPq0Z/2\naFP3AjcVO3tP1ks9IkY5bi3HT2uE31UslKlHcIwi7zr7rzegUW9N7Whix/for9ijVb5Lp4Db\npR753UULDY5G9LgY9R8o24XIK1WRrTXq5gJeWmwo8OjPiUf/QtleATh6tE+ZUh79yVujqBH8\nZWc9fTyi5JzKteVAuvD+DYLNX+vUO8vfdNc3xyN9+EsjrEcyW5OShT2SR7rlHo/0EOhh5MIP\nHZ4m4JEriq9Kqm2lUX+eVJF3nLWhlZpNeeOCLJ0fZTRatDbSiF2jxZjzw/WIdMf2iIzLufBD\nh6fJeeSR0SDZtFgqR6voe05Nhdr00oIXilTi0ZLz6M/86Xq0WY/0k8Sjv4BH9Db3v7cIFaoK\nFQhS2UfwI1DRxtVa/dh/DYn3DKfVBSIxMruGP/q8RqvVSGdvasOqDEdrMBz9qeYZjcRRI+9Q\nxHLALEsLFwNwtHKbMsWKdJkS65wWoV4Ku+DyMpEKNbJN1v+OGtk2q/wgq9WIfKyPSLnjukGG\npQ6jGkXIXcO9NGoF167KbniOMQbC6ufIm0QiFjEkovyXbP+f8/J/Rqk/t9E/5Ev/RfnjN2Xw\n93yLyjgrmddLQrr07WdrLLumatd5Hem8QY0gqqUMki9nmrIUWgr6+Puvz2U/Fau8BwDvCr5Y\npM+Hs850mQyDcd7CBvz9N6mciSs9ctyQa4sNANwBRAKgARAJgAa8cEEWgOuBSAA0oLtIe0Xu\ntpvoA3ANvUX6QCQwA93XkWS3N937G4CLgEgANAAiAdCAa0S664vGALiIK4oNN371JQDXcEX5\nm3HTBogExgYLsgA0ACIB0IB33o4LgIt56Q0iAbiWt96yGIBLeetN9AG4FIgEQAOQ2gHQABQb\nAGgAyt8ANAALsgA0ACIB0ACIBEADniISAGNDr+bbPErwyEG5DDDEEcb4niE+8kweOSiXAYY4\nwhjfM8RHnskjB+UywBBHGON7hvjIM3nkoFwGGOIIY3zPEB95Jo8clMsAQxxhjO8Z4iPP5JGD\nchlgiCOM8T1DfOSZPHJQLgMMcYQxvmeIjzyTRw7KZYAhjjDG9wzxkWfyyEG5DDDEEcb4niE+\n8kweOSiXAYY4whjfM8RHnskjB+UywBBHGON7hvjIM3nkoFwGGOIIY3zPEAc4EwCeD0QCoAEQ\nCYAGQCQAGgCRAGgARAKgARAJgAZAJAAaUe+Z7QAAAs5JREFUAJEAaABEAqABEAmABkAkABrw\nOJGOt957HAMMcYQxvmuIDzuP8F0sH8UAQxxhjG8b4rPO4sP/mrK7GGCII4zxdUN81kkUfHHm\nXQwwxBHG+LohPuskXvfrvYcBxvi6IT7rJF4X8O9hgDG+bogPO4m3TUFvYoAxvm2IjzuLh/9y\nfwwwxBHG+K4hPvo8ABgFiARAAyBSOQPkJCOM8V1DfPR5PJIBZskjjPFtQ3zsWTyVsqLoPQww\nxtcN8aEn8VwGWEkcYYyvG+JDT+K5DHAFjDDG1w3xoSfxXAbISUYY4+uG+NCTeDADzJJHGOPb\nhvjYs3gwD//33xlgjO8a4qPPA4BRgEgANAAineJNScmNPH2MSO36sf9an36ZPn8mP8KvEcWG\nnvx+qR99Gdw9mAhD1JYf/2tE+bsr8gowh49kkNXOh/8asSDblQGuAIjUBIjUlQGugHFSO3P4\nRJDadcXMPl8zTb6FEX6NKDb0hF4Bd48lwbM1GuXXiPI3ANcCkQBoAEQCIILJ6lBsAOA0dn4E\nkebkQ7l7MGGGGKKJSRBpToa4SgcYorDbmPKtu48HXM9Tr03C84eog9EHIk3L86/SAYZo5kcf\niDQrz79KBxoiL/t8/ukAcA+2+A2RALgGiARAAyASAA2ASAA0ACIB0ACIBEADIBIADYBIADQA\nIgHQAIgEQAMgEgANgEgANAAiAdAAiARAAyDSO6Cf2x7goz7vA7/zdwCRbga/8/cAgW4Ev/v3\nAJFuBL/790DvC7rfscP5tofn3vfqFeB3+x5ckbwbxz35DnJvAL/a9+CJJOjtDdl3DAXnwG/2\nPXipnaB/QqTO4Df7HtIiIbfrCn6z7yEfkUA38Pt9DxDpRvD7fQ+Z1E5Ap47gN/seUiI9/tuZ\nRwe/2feQFAkLsn3B7xaABkAkABoAkQBoAEQCoAEQCYAGQCQAGgCRAGgARAKgAf8Da6siV1to\nm7wAAAAASUVORK5CYII=",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "set.seed(1)\n",
    "seird <- seird_generator(S_ini = 1e6, E_ini = 1, epsilon = 0.01, omega = 1)\n",
    "x <- run_model(seird, t = 0:2000, n = 100)\n",
    "plot(x)"
   ]
  }
 ],
 "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
}