notebooks/lloydjansen2004/r_odin.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Metapopulation SEIR model with cross-coupling and/or migration\n",
    "    \n",
    "*Author*: Constanze Ciavarella @ConniCia\n",
    "\n",
    "*Date*: 2018-10-02"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Requirements\n",
    "To use odin, we need to install the package along with its dependencies. This is done using the following commands:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: drat\n",
      "Installing package into ‘/home/simon/R/x86_64-pc-linux-gnu-library/3.4’\n",
      "(as ‘lib’ is unspecified)\n",
      "Installing package into ‘/home/simon/R/x86_64-pc-linux-gnu-library/3.4’\n",
      "(as ‘lib’ is unspecified)\n"
     ]
    }
   ],
   "source": [
    "if (!require(\"drat\")) install.packages(\"drat\")\n",
    "drat:::add(\"mrc-ide\")\n",
    "install.packages(\"dde\")\n",
    "install.packages(\"odin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also define the following helper functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# beta matrix: random initialisation -----------------------------------------------------\n",
    "beta.mat <- function (nr_patches) {\n",
    "  ## beta: positive, non-symmetric matrix of effective contact rates\n",
    "  ## values are higher on the diagonal (transmission is higher *within*-patch),\n",
    "  ## values decrease gradually further away from the diagonal (*between*-patch transmission)\n",
    "  beta <- matrix(0, nrow=nr_patches, ncol=nr_patches)\n",
    "  for (i in 0:(nr_patches-1)) {\n",
    "    ## superdiagonal: scale vector s.t. numbers decrease further away from diagonal\n",
    "    rand.vec <- rnorm(nr_patches-i, 10 * (nr_patches-i)/nr_patches, 2)\n",
    "    beta[row(beta) == col(beta) - i] <- rand.vec\n",
    "    ## subdiagonal\n",
    "    rand.vec <- rnorm(nr_patches-i, 10 * (nr_patches-i)/nr_patches, 2)\n",
    "    beta[row(beta) == col(beta) - i] <- rand.vec\n",
    "  }\n",
    "  return(beta)\n",
    "}\n",
    "\n",
    "# mobility matrix: random initialisation -----------------------------------------------------\n",
    "mob.mat <- function (nr_patches) {\n",
    "  ## mobility matrix (origin-destination matrix of proportion of population that travels) #TODO trip counts or relative?\n",
    "  C <- matrix(0, nrow=nr_patches, ncol=nr_patches)\n",
    "  if (nr_patches > 1) {\n",
    "    diag(C) <- -sample(1:25, nr_patches, replace=TRUE)\n",
    "    for (i in 1:nr_patches) {\n",
    "      C[i, which(C[i, ] == 0)] <- rmultinom(n = 1, size = -C[i,i], prob = rep(1, nr_patches-1))\n",
    "    }\n",
    "    C <- C/100\n",
    "  }\n",
    "  return(C)\n",
    "}\n",
    "\n",
    "# plotting --------------------------------------------------------------------#\n",
    "plot.pretty <- function(out, nr_patches, what) {\n",
    "  \n",
    "  # plot total densities by disease status --------------------------------------#\n",
    "  if (what == \"total\") {\n",
    "    ## compute total densities by disease status\n",
    "    S_tot <- rowSums(out$S) / nr_patches\n",
    "    E_tot <- rowSums(out$E) / nr_patches\n",
    "    I_tot <- rowSums(out$I) / nr_patches\n",
    "    R_tot <- rowSums(out$R) / nr_patches\n",
    "    \n",
    "    ## plot total densities by disease status\n",
    "    par(mfrow=c(1, 1), las=1, omi=c(1,0,0,0), xpd=NA)\n",
    "    plot( t, S_tot, col=\"green\",\n",
    "          type=\"l\", xlab=\"Days\", ylab=\"Densities\",\n",
    "          main=\"Total densities by disease status\")\n",
    "    lines(t, E_tot, col=\"orange\")\n",
    "    lines(t, I_tot, col=\"red\")\n",
    "    lines(t, R_tot, col=\"blue\")\n",
    "    legend(-8.5, -0.3, title=\"Disease statuses\", horiz=TRUE,\n",
    "           legend=c(\"Susceptible\", \"Exposed\", \"Infectious\", \"Recovered\"),\n",
    "           col=c(\"green\", \"orange\", \"red\", \"blue\"), lty=1)\n",
    "  }\n",
    "  \n",
    "  # plot densities of some patches by disease status ----------------------------#\n",
    "  if (what == \"panels\") {\n",
    "    ## define the plot panels\n",
    "    if (nr_patches >= 6) {\n",
    "      mfrow=c(2, 3)\n",
    "      panels <- as.integer( seq(1, nr_patches, length.out=6))\n",
    "    } else if (nr_patches >= 4) {\n",
    "      mfrow=c(2, 2)\n",
    "      panels <- as.integer( seq(1, nr_patches, length.out=4))\n",
    "    } else {\n",
    "      mfrow=c(1, nr_patches)\n",
    "      panels <- 1:nr_patches\n",
    "    }\n",
    "    par(mfrow=mfrow, las=1, omi=c(1,0,0.3,0), xpd=NA)\n",
    "    \n",
    "    ## plot the disease statuses of some patches\n",
    "    ymax <- max(out$N[, panels])\n",
    "    for (i in panels) {\n",
    "      plot (t, out$S[, i], col=\"green\", ylim=c(0, ymax),\n",
    "            type=\"l\", xlab=\"Days\", ylab=\"Densities\",\n",
    "            main=paste(\"Patch \", i))\n",
    "      lines(t, out$E[, i], col=\"orange\")\n",
    "      lines(t, out$I[, i], col=\"red\")\n",
    "      lines(t, out$R[, i], col=\"blue\")\n",
    "    }\n",
    "    title(\"Densities by disease status\", outer=TRUE)\n",
    "    legend(-130, -1.2, title=\"Disease statuses\", horiz=TRUE,\n",
    "           legend=c(\"Susceptible\", \"Exposed\", \"Infectious\", \"Recovered\"),\n",
    "           col=c(\"green\", \"orange\", \"red\", \"blue\"), lty=1)\n",
    "  }\n",
    "  \n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model is specified in odin as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Compiling shared library\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "gcc -I/usr/local/lib/R/include -DNDEBUG   -I/usr/local/include   -fpic  -g -O2  -c odin.c -o odin.o\n",
      "g++ -shared -L/usr/local/lib/R/lib -L/usr/local/lib -o odin_b0e443ee.so odin.o -L/usr/local/lib/R/lib -lR\n"
     ]
    }
   ],
   "source": [
    "SEIR_cont <- odin::odin({\n",
    "  nr_patches <- user()\n",
    "  n <- nr_patches\n",
    "  \n",
    "  ## Params\n",
    "  lambda_prod[ , ] <- beta[i, j] * I[j]\n",
    "  lambda[] <- sum(lambda_prod[i, ]) # rowSums\n",
    "  \n",
    "  mob_prod[ , ] <- S[i] * C[i, j]\n",
    "  mob_S[] <- sum(mob_prod[, i])     # colSums\n",
    "  mob_prod[ , ] <- E[i] * C[i, j]\n",
    "  mob_E[] <- sum(mob_prod[, i])\n",
    "  mob_prod[ , ] <- I[i] * C[i, j]\n",
    "  mob_I[] <- sum(mob_prod[, i])\n",
    "  mob_prod[ , ] <- R[i] * C[i, j]\n",
    "  mob_R[] <- sum(mob_prod[, i])\n",
    "  \n",
    "  N[] <- S[i] + E[i] + I[i] + R[i]\n",
    "  output(N[]) <- TRUE\n",
    "  \n",
    "  ## Derivatives\n",
    "  deriv(S[]) <- mu - mu*S[i] - S[i]*lambda[i]         + M[1] * mob_S[i]\n",
    "  deriv(E[]) <- S[i]*lambda[i] - (mu + sigma) * E[i]  + M[2] * mob_E[i]\n",
    "  deriv(I[]) <- sigma*E[i] - (mu + gamma)*I[i]        + M[3] * mob_I[i]\n",
    "  deriv(R[]) <- gamma*I[i] - mu*R[i]                  + M[4] * mob_R[i]\n",
    "  \n",
    "  ## Initial conditions\n",
    "  initial(S[]) <- 1.0 - 1E-6\n",
    "  initial(E[]) <- 0.0\n",
    "  initial(I[]) <- 1E-6\n",
    "  initial(R[]) <- 0.0\n",
    "  \n",
    "  ## parameters\n",
    "  beta[,] <- user()   # effective contact rate\n",
    "  sigma   <- 1/3      # rate of breakdown to active disease\n",
    "  gamma   <- 1/3      # rate of recovery from active disease\n",
    "  mu      <- 1/10     # background mortality\n",
    "  C[,]    <- user()   # origin-destination matrix of proportion of population that travels\n",
    "  M[]     <- user()   # relative migration propensity by disease status\n",
    "  \n",
    "  ## dimensions\n",
    "  dim(beta)        <- c(n, n)\n",
    "  dim(C)           <- c(n, n)\n",
    "  dim(M)           <- 4\n",
    "  dim(lambda_prod) <- c(n, n)\n",
    "  dim(lambda)      <- n\n",
    "  dim(mob_prod)    <- c(n, n)\n",
    "  dim(mob_S)       <- n\n",
    "  dim(mob_E)       <- n\n",
    "  dim(mob_I)       <- n\n",
    "  dim(mob_R)       <- n\n",
    "  dim(S)           <- n\n",
    "  dim(E)           <- n\n",
    "  dim(I)           <- n\n",
    "  dim(R)           <- n\n",
    "  dim(N)           <- n\n",
    "  \n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the model\n",
    "\n",
    "We first define the user-specified parameters. The default values provided below all satisfy model requirements, but can be changed if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set seed\n",
    "set.seed(1)\n",
    "\n",
    "# SEIR free parameters --------------------------------------------------------#\n",
    "## total number of patches in the model\n",
    "nr_patches = 20\n",
    "## relative migration propensity by disease status (S, E, I, R)\n",
    "M <- c(1, 0.5, 1, 1)\n",
    "## matrix of effective contact rates\n",
    "beta <- beta.mat(nr_patches)\n",
    "## mobility matrix\n",
    "C <- mob.mat(nr_patches)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we run the model. Note that each patch is initialised with the same population size. We also run a test to make sure that the total population size has remained constant throughout the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run SEIR model --------------------------------------------------------------#\n",
    "mod <- SEIR_cont(nr_patches=nr_patches, beta=beta, C=C, M=M)\n",
    "t <- seq(0, 50, length.out=50000)\n",
    "out <- mod$run(t)\n",
    "out <- mod$transform_variables(out)\n",
    "# error check -----------------------------------------------------------------#\n",
    "if ( ! all( abs(rowSums(out$N) - nr_patches) < 1E-10 ) )\n",
    "  warning(\"Something went wrong, density is increasing/decreasing!\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the disease dynamics in the total population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAElBMVEUAAAAAAP8A/wD/AAD/\npQD////SDJfdAAAACXBIWXMAABJ0AAASdAHeZh94AAAdVklEQVR4nO3di3aqOBiAUXpa3/+V\nZ1pvIQRE+SWY7L3WdLwgQeU7odba4QRsNtTeAGiBkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCA\nkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCA\nkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkNYY7h4v+eCCh7d7dIsVa/xb\nZOXIhPBYryEkHvBYrzEb0gsXzI5wXfaIIb16t/rh4ViruOuEhfTEsmtDiiSkRzwcawnprSN8\nOg/HWumucz3Gux3tDckl46WSJW+ruR8ijm53WVt+i2SpfA2Lg6bbON7s6aXTBUZXTe7naFOT\nu9ZvXt3e8adlO2G2PyWXnMaLjffEdNlsTZO9c7pUYb2n2UHT9YyzH21EcYHsjk7uZymkIb95\nV3q938+bhHTfhyZ772ixIb/F6HS6M45PFpbK13u7xak8aD5EYSPKC9zWNWplcj/zKzs+4uv1\nfj/vto+U9qHb5ePckq/XxWZuXQgpv8VpbsX5vj+9bjraad1ApavmQxrdujO93u/nZSFN9uFC\nSMml2QLJP+KTGxduMXqW0hWfpiGN1zYaYnQ3Hi+Qb97Mpk7vTI96vu/PWQopP+TJFrsdEmV7\n4Ol+KJVdl94iWWoy/uke3cy2ZWvJzswucF3baLDJ/Zy5Mz3q9o4/bSGk8x70QkijvXg2pGwX\nzUIalkM6JQMPo9NLC9xXlxUzvp8zd6ZHvd7v582HVNh7J4tlq8n+2Z87tBuNXlrxipCSIUpL\nFReYDDxzP2fuTI86vdsvyEKaTCTlkNJddGYPXAhpcotTsspkPsk2cqnN8VQ0N2Z+R1eENImq\nL53e7ReM9+RhvJcN1x06+9f5funoFukOP7fXp7dIf4abrvd67Wgrk0HzIZLNOeWnR2PeV11a\nV7qps5vZmU7v9gsme/Lt5GTvHS2W7nGnaUjpxbe15beY7N9PhTTa2DyA8gKjdc3cz0Lv2cb0\npNf7/bx0H7nvMPf9KT+0Sy9NbpGGlByY3XfGyS1Kh2/jK7PNXPEWoXygyQL5wPn9vC6arCfb\nzM50e8eb0e++eyiehU8npEPwLHy2no+mDsWz8Nl0dBCehs8mpIPwNEAAIUEAIUEAIUEAIUEA\nIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEA\nIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUEAIUGAHUIa4MO8sJfHh1NhCIgk\nJAggJAggJAggJAggJAggJAggJAggJAggJAggJAiwX0ijdyM9eHOSkPgwO4WUva3v0dv8hMSH\nqRVSNkMFDAEVVTm0O59cKElIfBghQQAhQYDjhLTq1w3/JV7bDHiH44T07BCC4kA+N6QrLXEA\nnx/S6bel1zYJotR5Z0P0z5GURGV7h3T9GvzOBiVRV52Q4t9rpySqaubd30qipmZCOgmJihoK\nSUnU005IpiQqEhIEaCkkJVFNQyGZkqhHSBCgqZCURC0thWRKohohQQAhQYC2QlISlTQVkimJ\nWoQEAYQEAYQEARoLSUnU0VZIpiQqERIEEBIEEBIEEBIEaC0kJVFFYyGZkqhDSBBASBBASBBA\nSBBASBCguZCURA2thWRKogohQQAhQQAhQQAhQQAhQQAhQQAhQYD2QlISFTQXkimJGoQEAYQE\nAYQEAYQEAYQEAYQEAYQEAYQEARoMSUnsr72QTElUICQIICQIICQIICQIICQIICQIICQIICQI\nICQI0GJISmJ3DYZkSmJ/QoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIATYakJPbW\nYkimJHYnJAggJAggJAggJAggJAiwV0jD/+bPhgxxJyT2tlNIwzBKJzsbMkRCSOxtt5DO/13P\nDZcvgUMkhMTe9gnp3MytnOxsyBApIbE3IUGAKiFlR3ohQ6SExN5qhTR9sWFIPT9ESkjsrdKh\nnRcbaEuT3yMpib0JCQI0+XMkIbG3vd/ZcJ2M3vrOBiGxt73fa3fJ573vtRMSe2vy3d9CYm9C\nggBCggBCggBCggBCggBCggBCggBCggBCggCNhqQk9tVmSKYkdiYkCCAkCCAkCCAkCCAkCCAk\nCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCNBqSEpiV42GZEpiX0KCAEKCAEKCAEKCAEKCAEKC\nAEKCAEKCAEKCAEKCAEKCAEKCAEKCAEKCAEKCAM2GpCT21GpIpiR2JSQIICQIICQIICQIICQI\nICQIICQIICQIICQIICQIICQIICQIICQIICQI0G5ISmJHzYZkSmJPQoIAQoIAQoIAQoIAQoIA\nQoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIAQoIADYekJPbTbkimJHYkJAggJAiwV0jD/+bP\nhgyRExL72SmkYRilM2TnI4aYEBL72S2k83/J2ZOQaMc+IZ2buZWzeFj34hBTQmI/QoIAtULy\nYgNNqRRS4bWGIfX8EFNCYj8VD+282EA7hAQBhAQBqvwcafBzJBqz9zsbrpPR8isKMSEpid3s\n/V67Sz57vNfOlMR+Gn73t5DKvthg7lEV0md6x77ABhtDOr+IEJ6WkDJqOLptIQ3nN/uElySk\nUxZP7Y3hga0hTX5DIkTXIYnnEwnpOAT0wYR0CAr6dJtfbHjLyw09haShJmx9+Xt4x4TUTUga\naoafI9UioqYIqQoRtcah3f5U1KCgFxvCtmc6xOsOGZKK2hTw8vfJy99rqahZLYd0tJJU1DAh\n7cRk1LaQkA76A9kDhaSi1gW82HDYd38XQ/r5X8za15NR+wJe/j7s7yOVQvqtaOeUZNSDln8g\nWwjpktCOIcmoD32FdJ2K9puSZNSJ10O6vO37Ld8kvSmkez87lWQ66kZfIRVPvo2MOtLToV06\nDe0wJcmoJ9t/jpT+P8qbQpo98wamo750FNJ4EnrzlCSjzmwIKfwvg02H2CgL6bR0NpTpqDtB\nIYVu03tCyqegN05JMupPzKFdtKiQ0pIm3bwtJB11qOmQ0ilpOgG9aUpyWNelpl/+HoU0vfYt\nIcmoT03/QPZRSG8oSUed6iWkYjTxIemoV70c2hWbCQ9JR93qO6TgknTUr82/IfuGXzR/Q0gz\nycSGpKOObX35+8h/aCwNqXx9aEg66lnAh5+84RMidwspsCQdda2PkGaDiQtJR33rJKS5BcKm\nJB11Luizvz81pKgpSUe9C/g4rje85S4spHNJC/NOTEg66l7bP0e6TEkLtYQc2+mI3kOKmJJ0\nRNRHFh/0e6RdQtIRAT+Q/YCQFg/fNh/b6YhT67/YdwlpcZGNIemIX0LaGtK2m9MIIW07tjMh\n8Sfg3d9vEBrSo1K2hKQjzjp4seFRKBtC0hEXQtpQko646uAHsm8M6dUb0hwhvf5ygwmJm7bf\ntPr7rtUVmbwWko64C/o1irDtmQ6xzb81lbwUko5IBPxi3+m4v9i3NqRXStIRCSH9eiEkExIp\nIf16PiQdMRIS0nF/jnT6t+qw7fljOx0xEvT7SGHbMx1im3/rEnk2JBMSYwEvfx/4k1bfFZKO\nyDT+A9nVIT1Xko7ICOnsqZBMSOQi3rQatzWFITZ6R0g6YmJLSG/7o+Y1QnqiJB0xsSGkW0BH\nftVudSDrQzIhMfV6SEk+B/450s+/x8ucF1xdko6Y2hJSctHnh7R6SjIhUSCk25LrFtMRJUK6\nLbmuJB1RIqT7omsWMiFRJKT7omsW0hFFjYf080xIK0oyIVG26eXv4V0/ko0L6bQ+pDVTko4o\nE1K68CMmJGY0/qbV50J6VJKOmNN+SIFTko6Y03xIgVOSCYlZQsoWX6AjZu0V0vQFiaUXKA4Z\nkgmJeTuFNH1pb/GVvqiQfp4NaakkHTFvt5CyKWj5JfOwkE7PhbQ0JZmQWLBPSOdoknT+P3nM\nkOZL0hEL6oT06JP3q4U0PyWZkFgipOktinTEkioh/Z2YvooX/46jF0KaKcmExKIaIU2+Y4oY\nouj5kOamJB2xqE5Ij+ad0JAiju1MSCxrP6SQKUlHLKv1c6RdDu3OTQSEZELigb3f2TB+5S5y\niILXQiqUpCMe2Pu9dscPqTAlmZB4pOl3f78Y0qQkHfGIkOZudmdC4iEhzd3sTkc8JKTCzbKS\nhMRDQpq93ZUjOx7rIaSNU5KOeKyDkDZOSSYkVmg5pGsOL4SUlKQjVhDS4k1NSKwjpPJNbyXp\niDWEtHxbExKrCGnmtpcb64hVhLR8YyGxipDmbvx3a0d2rNNFSK9PSTpinR5Cen1KMiGxkpCW\nbq4jVhLS/M1/TEis1XBI93cnvBbS6UdHrCWkhRUIibWENM+RHasJad7Xwz9zDhdCmvX/hKQk\nVuojpJdK+jqdJp+DAmVdhPTalPQXkpJYRUhzzi81CIlVhDTn/JKdKYlVhDTj+tq3klhDSDNu\nP0MSEiu0G9JPVEhK4rGGQ0pOPx9S8q4GJfGYkMrSdwcJiYc6CWlVSd9/zqdHISmJR/oIacWU\ndG3oHNP4/apK4hEh/bnNRJdz2fu+lcQDQjrlGZ1+X2rILhESy4T021F+yVcelymJZUIqdHR+\nqWGUkpJYJKRSR1/Tq5TEku5DKnR0f+07nZSUxAIhFS5LXrNTEqv0EtJcScUJ6au8gJKY1WxI\n2V4/E1Kpo/zzVe/LCIk57YY0PlsOqdxR/ilct6VMSczpPKTShdNPs1MSj3QdUnFCKn1y/reU\nWNZ3SKULy5+vqiQW9RzS6glptLCSKOg6pOIN5z7wW0ks6CakaUkzE9LsJ+ffvlFSEhO9hFSY\nkp6bkP5uoSRm9BtSeUJa/muX95KkxEjHIRVv9uBvIvlGibJuQ3plQjqNvlGSEnf9hlS+2eM/\n0mdSoqDXkOYmpBV/7TIpSUpctBrSZBfPQyrfbNVfjU3eDq4kzvoJaVTShgnpNH7rnZT41U1I\n4ylpy4T0d/tkUpISQhpbHdL4F2elRJ8hbTuyu6xCStx1GlL5Rk90dMo/rEtLfesypJkJ6cmQ\nstVIqWt9hlS+zTNHduf1TFLSUq+EdPdsR6fp3KalXnUU0q2koCO7y7rylf2IqUf9hHSfkqKO\n7C5rK6xOTN0R0s1rHZUmpfMWiKknHYYUeWS3vMZzTHLqQaMhlXbeW0jlm7x4ZHde5WxK542R\nU/NaDalw2aOQNg24WNKfnx9BNay/kOKP7C6rfZjS2U9i25AcSIchlW+x5cjusuK1Kd38FGzd\nCuroKaRzSe+ZkM6rfrakqVJbHMjcE9dRSOcp6U1HdmfPT0o0or+QyjfYfmR3JqVOCeksqKOT\nlDrVW0hvPbK7kFKHugupvHzUkd2FlLojpD+xHZ2k1J29Qhr+N382ZIjUXEh7HNldfGupJzuF\nNAxpOsP4bMwQI+WX+//tdGR3JaV+7BbS+b/7ufR8yBAjz4a0dbw5pqVe7BNSsZw3hjTzA+h/\nOx7Z3UipC/VCWjq22xxS+eJ/+x7ZXZmWOlAtpLd+j/RsSBuHe0xLrasV0rSjIfX8ECOHC+mk\npcZVCulBK28KKeCTijf5FlOz6oT0aM55V0jTP23+a6+OztsgpiZVCenhsVvDIf0SU3sq/hwp\neIhUOaTvcki7HdmNfKupKXu/s+Hy5cFLCm8KafIXmf9U6ehMTc3Y+712Qpr4llMD2nz39zMh\n1Tmym/jW00frKKTvmZA2jhXr+1tRH6mnkPI/bX52rJBuvlO1N4aHmgyp/J7V391xGtJBjuwW\nfU/V3iQybYZUuvB7JqRtQ1VSSKuo9nb2o6OQfr+0EtJKa3tjtbmHuveQPuHIjg/QfUjbRoKz\nbkL6FhJv1E9If1/z178d2RGjs5DyKUlHxBASBOglpO9iSI7sCNJNSJf/ZyFtGweuhAQBWgyp\n8Fa720+kRy/bObIjSpMhTS+6v7NjFNKmYeBOSBCg55Ac2RGmj5C+yyFtGgUSnYR0Pykk3qG/\nkG4lObIjTnch3ackHRGni5C+hcSb9RFSeuYakiM7AnUc0qYxYKTbkExIRGowpMlb7cYf/XJ5\n2U5HRGoxpPyC7COUhES8XkNyZEeobkPaMgLkOggp/3RMIRGvh5Cy87+vNjiyI1aHIf1OSToi\nlpAgQPshTf+AwD9HdkTrIKTJAv9MSERrL6T8jQ2FkP4JiWANhpSdL/xpKB0RrfmQSn9jTUhE\naz+k6RJeaiBclyHlfyUJtuoypMLfZIZNWg+p8C3Sl5AI13xI0wW+Cn9KFrbpL6QvIRGvuZDG\nP48tHdmdpn+TGTZqL6TRubkfIgmJWN2F9CUk3qC/kP6+OrYjVtshTb9Fur6rQUiEajykydXX\ndwcJiVBCggCthfSzHNLt/aq+SSJUcyGlZwrfIt1OCYlIbYeUX/slJN6js5DuJx3bEamrkEa/\n0SckArUc0uRbpNFvxgqJQE2HlF03/hVzx3YEaiykn8WQxmeFRJzWQkpO50d2+WeemJKI03JI\n2XWTzw4SEmH6CWn6IVxCIky7IU2O7CYLO7YjTFshpa81PJyQTEnEaSyk5PTDCcmURJxmQ8qO\n7MofUywkgrQb0via8sd9m5II0klIc5+bryRiNBXS7JHd/N+fEBIhmg1pdMX833ExJRGii5CW\n/iCSkIjQUkjJT5GyI7ulWymJAE2FdD+5ekJycEeINkMaTUiP/tKlktiuoZDSI7vk4sd/MVZJ\nbNZSSLdT6YS05i8vK4mt2glpbkJac1slsVFDId1OPTshnX5LkhJb7BXS8L/5swFD3CekFzo6\nmZTYZqeQhmGUTnY2YIhiR1/rOzrypPSP45h7knYL6fxf+WzAEOWOnlvH4gO1XfyTx3HsE9I5\nmls62dntQ/zcJqR7R09NRzdrdl89MNFESElG145eyyilCJ7w+SH9XKej71tGX9szgqccJ6Qh\nNbeir3We30LY5DghbRwCahISBBASBGjl50hQ1d7vbLhORtHvbICq9n6v3SWf+PfaQU3tvPsb\nKhISBBASBBASBBASBBASBBASBBASBBASBDhoSPBhXtjL48P5iLEPsQHGb2d8IRnf+Adb1yeN\nfYgNMH474wvJ+MY/2Lo+aexDbIDx2xlfSMY3/sHW9UljH2IDjN/O+EIyvvEPtq5PGvsQG2D8\ndsYXkvGNf7B1fdLYh9gA47czvpCMb/yDrQu6JSQIICQIICQIICQIICQIICQIICQIICQIICQI\nICQIICQIICQIICQIICQIICQIICQIUDGk1z71P2rsmttwH7XSY3D/iwvVnoSK46d/cSJs/Hoh\nvfr3M4JGrrcNyfNY6TG4b0G1J6HmA/CWJ6BmSPeJYe+R7yFV2IbrU1hr/PobcJ8IKj4Bt5Mf\nHlLybFYevVbN1cbP9+EqKd863n/8+3iB4wupzjYMNffj26gVS64aUjohCilq9EqHNgcYv+uQ\ngscXUo1tqLsfn27fo1Q8tKo/fuwhgZCqfItwf62j4isulTYge62jkW9ShVTjyGaYbsXuKoZ0\nO7YSUsDA3YZ0/8FFn6/aVQ/pLfe/WkjJHak0eK1tSAas8xjcv0eo+CTke/PeQ0ff/5ohVXyP\n0C2kej9YH2qNP379t97bS6qN39g7G061nsLz2PW2YcieyL3HHw1b7UlI/imrMHb8/a8YErRD\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBA\nSBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASBBASJ9hSP9gJsfj\nifkMg5KOzfPyGZI/Hlx3QyjztHyGWz9COiZPy2e49/N36n6UN/4qs1o88J9hHNJw+5KG5Huo\nijzunyEP6XbZ/WRylt153D9Ddmh3PqGh4/Dgf4bi90h5SI7s6vHAf4bC90j5Ud0lpTrb1z2P\n+2cYv/ydvlT3+xQO0+XYl8f9M4x/IDsOKX0hXEiVeNw/w/gtQteXv68lJYvU28S+eeA/Q/ZW\nu2EYfWuUXkgVHvnPJ58D8Bx8PiEdgOfg0zmeOwRPwqcT0iF4EiCAkCCAkCCAkCCAkCCAkCCA\nkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCAkCCA\nkCCAkCCAkCCAkCCAkCBAHtIAPD+/TEJ6OUlohpAggJAggJAggJAggJAggJAggJAggJAggJAg\ngJAggJCqGr9T64U3bHEQQqpq/J7HiJAm61DnLoRU1XUvf+Xdw8trnL+AdxBSVbe9XEgfTkhV\n3ffy31PXA7zbt0zJiWTuyq5NL78uebvF5UR66FhYw0u/T0NKSFUVQsobuBz3DfcKxtdeb34J\n5vq/+82ykAprePE300gIqappSJd9/JTs+6MKppdNr8wmoNNkZYX1C2kbIVU1E9L9quTscFq6\ntrTG2ZBO5Yt4nZCqKoU0OvJKT57K157Sy5Pvg+ZCmqzBkV0AIVU192LDcJpMLekLD8UX55LF\nhqyeyYsN2RqUtJmQqip8k3M9lx+tlZorrOzR90hzaxDSNkKqKnn9+TR5MaDwwkLh2tP0ypmQ\nhqX1C2kbIVU13JzPjb/buX8Pcz5XvDZZz/VUeovbBZNviHyPFElIVY06Ot13/NuVyYm5a7Nz\no+zGFyyvny2EBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGEBAGE\nBAGEBAGEBAGEBAGOF9L2DxDIPmxq188jGH+ayetridqerQMv35GDfdRDyGNfXO+KZZ5f6+Y1\nLK5920ORfYTV/bK9NBbS/B054md4xTz4pfWuWOb5tW5ew9LKkw++ee32+YndQzrOWgIGnt+Q\ngzX0Z/wpgfHrXVzm+bVuXsPSyofx/1+9vZBiBv7MkKI37iNDGtIz6dfkM9jGn8Z2+fzD63XD\n9dAuuT5dwXuVtz/ZmvQelD5Z7v4BdDVcD9iyT46cbHzySCd3qMLjPd365ETpMV/c2JnP/vvA\nkNK9LL2ro6d0KF40fnqHdMF06TfLPrTx/iV5ArOTo61Ll6xg/IBmj2uy8WlI2SO97+M93fpT\n+g3C5DFf3tjRzpPc8ANDGv/jcf2aBpFedL3bp+xxGV04XvrNhvShTza9sInFrUsWryHf2uXt\nLn7JnrB9tz5/9Ndv7PyutO7peEtI/9aZHyJ/JPI5u3Rlad+8nXwqpJ915jZ9SP8pGA1c2LCF\n+/Ky73XKmz/anPIuN15ycoeyBZ/2tU5x68cP/nRr8n1lNqmnn47jzUh/65zcwdF8fFnkduBx\nui63ENJOxxrZEONnYrJh4+ONqJA2mPsXebx5o6tnQqpwZHfbvOt2TDd7JqTxgi89HccMafKE\nnuYektP0wZmfkXbwTEjz9+poIY2uzpYshnS7a7vKQ8ovX5yRFm74gSGNHorsbn1gSMuHdvll\np6eeubeIC+lU4V6MR360r/zuucVw2ghpuP/vOu9OntTLwcPkVPHqyTJvNRoiPXiYbuLsxp72\n/7c83eLCdk43eenL9F/6Pbc+2fbF7Uz/mRv/o/DS03GwkE7jl13ux6rZZemp/Can7MLx0m82\nDEO6Yfcvs/egeK8OFVLpQb8/0sWNr3MnsrlwbrNvW3b//22zT8UbfmBIp/EdHYo/Icv/QRlf\nff23Jbl+tIa3God0HT7fxOweTE4dLKS5jU932Pty2RO2oyyk2cd8mC6f7iovPR2HC+l51fa6\n9T5gE9lGSHv4gE1kGyHt4QM2kW2EtIcP2ES2aSAkqE9IEEBIEEBIEEBIEEBIEEBIEEBIEEBI\nEEBIEEBIEEBIEEBIEEBIEEBIECAgJGB7SMALhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQB\nhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQB\nhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQBhAQB/gNoEeM5B0AwlAAAAABJRU5E\nrkJggg==",
      "text/plain": [
       "Plot with title “Total densities by disease status”"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot total densities by disease status --------------------------------------#\n",
    "plot.pretty(out, nr_patches, \"total\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the disease dynamics of some of the patches. Dynmics between patches may vary due to cross-coupling and migration rates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAElBMVEUAAAAAAP8A/wD/AAD/\npQD////SDJfdAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2di3ajMAwF3Tb8/y/v\nEt6JBDYIJMPM2W3TxrHhRgOGJDQ1AHCY5L0AAHcAkQAMQCQAAxAJwABEAjAAkQAMQCQAAxAJ\nwABEAjAAkQAMQCQAAxDJmPRm7f7ld/le4ddrnYI7PDvGpA2TBicuEGn2CCw8GwI25l2yG3Wr\n370i0s5FOdIB5EPAxnQl2+vUfx1+1+9Y3jdmLZppF9bfsbBx+TvhcctfjncM481HnLcHS0jU\nmEmkNFXy4sd5Wc/v6R43thpEWv4qffes9DIpNLZbdAKmkKgxM5Gaxd5ktGj247JFMxb6eG+/\nn/nso5nthpa/nPUy3CcsBFM9cwjUmFEkYUcw3D/bM/V7jcWjZw7MZnlTH7PzGdPD5yc5ZkON\nbT47AVvI1BhVpOkY6VOk5THSmkiL6dpi0rgQaTHUh83NfDgwg0SNmYk0/ULeN8zKOX1rI4o0\nP8IZ75EOe5I+tZsNB2YQqDFTeUsHReLRyscx0nL+9TG1mz3g8xjp80hLE2nWEMwgUGOmKdbi\nbNvHJOv7rN1sf9M/eNHhRx/CDO9jZjc1GI6j5jNBRLKGQI1ZWtF8HZrMRRp0+VCgv7XocXl4\nM31d7JwWD5nORPTfxs7x6ARINCg8MXXB8xUS9hm1wfMVEkSqDZ4vAAMQCcAARAIwAJEADEAk\nAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIw\nAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQ\nCcAARAIwIIJIKal/sP7Rf0qVXGRC5hLhCUkfycg31d/clsJcHhNMSS6qcuYLdcko67xXNT+Y\naxYqACW5PCmYglxWdl7WC3XFIBtMwbzXulv1NP4iLZo+p15KcnlUMEX1ctUcJkL6YzBpZEjk\ne4vynHopyiU9KJiCXJ4p0uz7EMB3EM+pl5Jc0hNFyqmXhx0jTaubpmC6Hx8tUnYul214I1BS\nL88UadxNI1JLfi7XzWAiUFAvD5zaDd8RaSQ/lzTfSN+egnpBJI6RSnJBJP3Y8ZEijWdgmvEs\nDCLl5fKgYApyuW77EiH9cVXHCW93IrNBpOHGdi4PCqYkl8v2049JH+BMEAnAAEQCMACRAAxA\nJAADEAnAAEQCMACRAAxAJAADEAnAAEQCMACRAAxAJAADEAnAAEQCMACRAAxAJAADEAnAAEQC\nMACRAAxAJAADEAnAAEQCMACRAAxAJAADEAnAgMMipTqxyI5cyGVar8PBHO3AhQsK5vQRzoBc\nZBBJgYKRIRcZRFKgYGTIRQaRFCgYGXKRQSQFCkaGXGQQSYGCkSEXGURSoGBkyEUGkRQoGBly\nkUEkBQpGhlxkEEmBgpEhFxlEUqBgZMhFBpEUKBgZcpFBJAUKRoZcZBBJgYKRIRcZG5HGd5FL\n7ye/azA5nZCL3MkTc8lYrzQ2S0L7uwaT1we5yH08L5cSkaYvZR0ExLRgyOWrj+flkrlezwum\noB9ykft5VC5565Wmr08JpqAbcpG7eVQuWeuVZt/GB2R/mj0iNktNLqu9PCqXrLN28+9P2cIU\ndEIuciePyiXzZMN46zFnYUr6IBe5jyflkiNSt0NOz3pdIKMPcpH7eGQuh9frrsHEH+EMyEUG\nkRQoGBlykUEkBQpGhlxkEEmBgpEhFxlEUqBgZMhFBpEUKBgZcpFBJAUKRoZcZBBJgYKRIRcZ\nRFKgYGTIRQaRFCgYGXKRQSQFCkaGXGQQSYGCkSEXGURSoGBkyEUGkRQoGBlykUEkBQpGhlxk\nEEmBgpEhFxlEUqBgZMhFBpEUKBgZcpFBJAUKRoZcZBBJgYKRIRcZRFKgYGTIRQaRFCgYGXKR\nQSQFCkaGXGQQSYGCkSEXGURSoGBkyEUGkRQoGBlykUEkBQpGhlxkEEmBgpEhFxlEUqBgZMhF\nBpEUKBgZcpFBJAUKRoZcZBBJgYKRIRcZRFKgYGTIRQaRFCgYGXKRsRJpaCX8LcO7BlPUDbnI\n3TwolyKRpMZ3DaakG3KRu3lSLjnrNW5WnhRMTifkInfyxFzK9kgP+ivVJd2Qi9zNk3IpndpN\nyaSefUvmjPkUhlyEbp6US9nJBuEBdw2mtBtykbt5Si6IZNQNucjdPCWX/VO7og6iceIUxnaE\niyEXGUuR3v+fc/CY3Q25yN08LJfD63XXYOKPcAbkIoNIChSMDLnIIJICBSNDLjKIpEDByJCL\nDCIpUDAy5CKDSAoUjAy5yCCSAgUjQy4yiKRAwciQiwwiKVAwMuQig0gKFIwMucgUiCS9L8pm\niIhkL/XeWMjl6AihyBcpDf/sh4hI7lLvjoVcDo4QC0RSoGBkyEUGkRQoGBlykSk5RiIYseHO\nWMjl6AihKDprt28V7xrMnqYGD3OGXGSK9khp13mYuwYzNNwZC7kcHSEUHCMpcCwgQy4yiKRA\nwciQi0zRC7KJYKSG+2Ihl6MjhKLkZMPOa/fdNZix5c5LGpLL0REikStSanZfBvOuwXSt9l8d\nlFyOjBCNotPfZw0RkfOXmly8RjiDopMNWe13DBGRgoPqgtZ7RogFuchkT+32X+H8rsG8Wx24\n8Du5HBkhGqV7pFOGiEjZlvfMEWJBLjKcbFDgoFqGXGQ42aDAQbUMucgUvSDLe6fkhry+Jjck\nF6kFbxHS2vFWGLkduYgtEElrR8HI7chFbIFIWjsKRm5HLnILPiErN9wZC7kcHSEUVmft+lZP\n+nPvJd2Qi9zNk3IZp3arLwyk8dtjgunbrb9eQi6r3Twpl/TxXWyT0tTms+Fdg8loRy7KnU/M\nJUek4d5HBZPVjlzW7n5ULiVTu0cF07fLmcKQi9zNo3IpOdmwDGb5DuDfjuIldML0oJpc5G4y\ncqklmCKRVhrnbmFqyaWkYFaPHqev67nUEszFuVQTTIFIqxcqyz0LU0suJa+XrF2/LfvsFLko\nI1QSTNHJBmm1Z61Sk/G6QCW5lBxUr8RCLivdZOVSSzBWImUPUUkuVgWTPUItu+qrc6klmJKp\nXWNxXbs6YimZwuyM5bZb3rGh0XXt6gim5GSDyXXtKtnAFNSA0fXbyEUZoY5grE5/53dQRy4X\nvJpBLnkj1BFMwTHSzo/g37xgdsdCLpkj1BFMrkjzN0cdG6KOXHJfh94fC7lkjlDHwUC2SNnt\nt4aoIpbsgilqvTpCHfVyfS6VVAwiKVxfMHUEg0gyiKSASDKIJJN/jGR1gcgqYsk/FrC7EGId\nwVyfSx3BcPpbgdO8MtfnUkcwiKSASDKIJOMgUhW5IJICIslcL1IduSCSAiLJIJICIskgkoyL\nSK/X0T7Px0WkCoLxEKkGkxBJAZFkHESqYpfkJBIFg0j5I7TJnD7sQXxEip+Ll0jhg/ESKXow\niKTgIFIVwfiIFH8Lg0gKLmenKgjGSaTwwSCSAiLJIJKMl0jRc3ER6RW/XvxECp6Mk0gUjCjS\nLQ6q7UeooWIQScFJpPDBIJIMIikgkgwiyXiI9IofCyJpIJKMi0gNxwLSCIikjNCJFDsZL5Eo\nGE2k4MG4iRS8YhBJwUuk6MEgkgwiKSCSjItIL0QSOujrJXYuiKTgI1L7NXbFuIlEwWgixQ7G\nT6TYFeMgUlNBLI4iBQ8GkWRsRJouBShcFPDBIpXlUkO9eORSRcUYiTQ0kxprIoXOxahginL5\nrSEYh1yqqBgTkdLyy1YHr8W3oFgUTGEuvzUE45BLFRVjLdL2X6muIhbzgsnJ5Yki5eRSRcXY\n75GmZLQLqVcQywlb3u1cnihSTi5VVIyxSNIDVJEi52JcMFKPukiRg3HIpYqKcRQp9AbGU6TQ\nwXiKFDkY+7N2uae/Q8dywtmpjFxqCMYhlyoqxkik98w2dbc2OxhjCZyL5esl2bkMVxSVgvlZ\nYrF0u3DIpYqKsRGprIPX142AmBRM2QijSN/BfKrjp5NDLlVUDCIp+Iq0TEZTxmP3hEgyviLF\nzcVXpEXBbJpy6YzPWaSwFeMqUuANjLdIUzIldvycirzU5qyJFLdiEEnBV6QpGcfzCiKIJOMg\n0mxrG3dP7S5Sn0w0j1xEqqFiPESSNrzhcBapTyacRz4iVVAx3iJ95LI2O78Wd5HaZOJ55C9S\nUJOcRfo8PbV+mvcibFZ7k3WR2mQCeuQuUtRdkrtIs6PqUGXjL9IrViA9iCTjLVIzHlQHqxpH\nkf5a/v8YLJEOf5FimuQpUlcu71yiaeQpUptKm83PT8SScRXpHU3EVDxFajVq//83KdzuqHEU\n6V0sTbtt+Yu48fURqU/mry+YgLiJ1NdLG0xAjTxF6n/6+Z9NwJpxEan9u+Z9xUQ1yU+k/qd2\nCnN0hDPwEmnaITUhTXIUqfuhNen0ZSjHRaTfqV7aign5t9/dROp/eG9d/uKp5CbS37TpbU9o\nur3AqOAh0juW8cefgNXSuIk03yG1hMvGT6Txx7/RoEAuuYvEFGbg9fr92CG1vDczgdIJIFJ7\nGmb2QwyVnEQa99M//RQmVLG0OIg0K5d5dbwPr8PE4yXS/FigrZfp3hgqeYnU355mMJGKpcVZ\npPkd/ettMQJyE2n4odvyvrLeW3Yh3iINv5wVS4BqabxE6re7n3Xx1//+NWP+yGsQl9qe72PH\n703v18GAu0o+In0eUrcMxXJVWYQrmP/j/ko7pDd/09Smb/21xNfgIVLzXTHfh9XOu6XrRWq3\nL4JHjVAsnniKJFbEX4x4vETqbv7MpjDfGxDNpSs2OA4irWx4Y9TKG1eRlEdEcMlZpOm33WH1\nZ9vhBSZhjnHax3Hkpd5er1JUkaTNR4RaeeMiUrfua1OUP++AXEVaJDMcVU/tJm/WSv6ipd7R\noqyDzQ3vn3uxtPiI9L6x9Xz/ZSb0Z4y41PasifSxfu97XQ4Vv3ARaXvD23w8feI9Z2Gz2pto\nImVuOHNXxBgfkeSK8d/cTviI9L6xt16OLk8ejiKdPvIRnER636j7mPoskbxP/G9wecG8+mDI\n5fMXbTDt+3nlY+rTlycPR5GO9nsu14vUVwW5fPz86kWStzAbOyXjWb8+g/YQ6a+CDa+PSGq5\nxMHpbObv+osCf9+/ufRAwEGk/9uXGurFT6TTxz2G07Hjb+65qUvtmXG9SG29/IavFy+R2MCo\nIp0+8CFsRBr/FvX2n3vvYglfLzYFU5DLe8r7+0sugkh/FVSMkUhDsyS0/xTpHUve0jliUzAF\nubzaczDRy+X6XGoJxkSkNHwZb6gdVBKLTcEU5NJtYMjl++fuLO/PbxOby0Xq6iV6LJcXzHsD\n8z+X8MFcvoF5vS80Fb5izhQp9XQ/Ld48+BuZ3NU+IZfQwfjlEjuYzFxs90jVcPkUphLIRQaR\nFCgYGXKRufysXS1cfnaqEshFxkik98w2NVmvC9SBTcGQi9zJE3M5vF53DSb+CGdALjKIpEDB\nyJCLDCIpUDAy5CJzhUh1cnS1yYVclutllc8VLS4ZxBZy2TtcZbkg0rmQy97hKssFkc6FXPYO\nV1kuiHQu5LJ3uMpyQaRzIZe9w1WWCyKdC7nsHa6yXBDpXMhl73CV5YJI50Iue4erLJeL8wO4\nJ4gEYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYICNSJsfxk0bjbq7NltsjLPV4n1n\n1geHjSCXlQFXG2w0ipeLkUhbHaX1RtPdaov3fevjpPUu+js2l9UQclFHvFsuJtmljZ5SSuuN\numA2W2yMk9a7mN15UcWQy9aAyv0V5nKJSJPWa+ll5LvaIG110e2pN5fVDnLRlmhrrPpyCSPS\nZouM7UdGi62tlCnkog54u1yiiJQyullfqe1gtvswhlyUwTbHqi+XICKl5Re5wUYw44X8blQw\n5LLVQ5RcbKJLWx2l9UZps5uug41x0vqSpGF/f1G9kMv6Aq03qCwXI5EOvi7QbxwueF0gY1kN\nIRd1wLvlcll2AHcGkQAMQCQAAxAJwABEAjAAkQAMQCQAAxAJwABEAjAAkQAMQCQAAxAJwABE\nAjAAkQAMQCQAAxAJwABEAjAAkQAMQCQAAxAJwABEAjAAkQAMQCQAA8KJlK68vlpFkItMlFwi\nLMOC/hqb8AG5yETJJcIyLBiutTz+jZzUxNjkOEMuMlFyCfdU9MFMF39OUbY5vpCLTJRcwj0T\nafwyDwbIRSZKLuGejdmuOvWbGqYwDbloRMkl3FPxuatO028fDbnIRMkl3DMxTXb7v+0R5ryM\nL+QiEyWXcM/E+GfU0nQehikMuWhEySX6UxF9+bwgFxm3XKI/IdGXzwtykUEkmeCL5wa5yPjl\nwjMCYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYAAiARiA\nSAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgE\nYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYAAiARiASAAGIBKAAYgEYAAiARiASAAG\nIBKAAe4ipTfKXRcvSyTIRSZqLu7PSfpIRr65+I0W5K3Yl8v9gynJJV1YL+7Jv1eyKBh9k3Qn\nyEWmIJc+kGtycU9+Cua9vuO6N98JpDQ2vn/F7MjFZTmvpiiX9/9rsnFPf+bGwGJTkr4bN48S\nKTeX9/1OC3shBbmMtl1RL+7RD+s5+z7scr53PeNP9y+Z4lxWjsLvRFEuzxJpev7TFEz3oybS\nIwqmMJeLCsabolzmxp29XGcPsLkAadxopFyR7l8uO3J5mkjBcnEPfrGTyRPp/tXS7MjlMSLN\nvm/kMkXyXJFWjpHuXywtO3JpnpBNQS6zze75sbgHP5vvDmdgmvEsjCLSE46qd+bisKTXUpDL\n8pze2ct1+ghbC5DGG/3mpTsD1ayfnXJf7rMpzmW49+bk5zIWyiW5PCB6gPNBJAADEAnAAEQC\nMACRAAxAJAADEAnAAEQCMACRAAxAJAADEAnAAEQCMACRAAxAJAADEAnAAEQCMACRAAxAJAAD\nEAnAAEQCMACRAAxAJAADEAnAAEQCMACRAAxAJAADDouU6sQiO3Ihl2m9DgdztAMXHnBR9V2Q\niwwiKVAwMuQig0gKFIwMucggkgIFI0MuMoikQMHIkIsMIilQMDLkIoNIChSMDLnIIJICBSND\nLjKIpEDByJCLDCIpUDAy5CKDSAoUjAy5yCCSAgUjQy4yiKRAwciQiwwiKVAwMuQiYyPS+C5y\n6f3kdw0mpxNykTt5Yi4Z65XGZklof9dg8vogF7mP5+VSItL0payDgJgWDLl89fG8XDLX63nB\nFPRDLnI/j8olb73S9PUpwRR0Qy5yN4/KJWu90uzb+IDsT7NHxGapyWW1l0flknXWbv79KVuY\ngk7IRe7kUblknmwYbz3mLExJH+Qi9/GkXHJE6nbI6VmvC2T0QS5yH4/M5fB63TWY+COcAbnI\nIJICBSNDLjKIpEDByJCLDCIpUDAy5CKDSAoUjAy5yCCSAgUjQy4yiKRAwciQiwwiKVAwMuQi\ng0gKFIwMucggkgIFI0MuMoikQMHIkIsMIilQMDLkIoNIChSMDLnIIJICBSNDLjKIpEDByJCL\nDCIpUDAy5CKDSAoUjAy5yCCSAgUjQy4yiKRAwciQiwwiKVAwMuQig0gKFIwMucggkgIFI0Mu\nMoikQMHIkIsMIilQMDLkIoNIChSMDLnIIJICBSNDLjKIpEDByJCLDCIpUDAy5CKDSAoUjAy5\nyCCSAgUjQy4yiKRAwciQiwwiKVAwMuQig0gKFIwMucggkgIFI0MuMlYiDa2Ev2V412CKuiEX\nuZsH5VIkktT4rsGUdEMucjdPyiVnvcbNypOCyemEXOROnphL2R7pQX+luqQbcpG7eVIupVO7\nKZnUs2/JnDGfwpCL0M2Tcik72SA84K7BlHZDLnI3T8kFkYy6IRe5m6fksn9qV9RBNE6cwtiO\ncDHkImMp0vv/cw4es7upKpcfnff9T81lCyuRzuzAhfOX+vpcViT50EXnjrlYgEgK9ymYfEly\nuE8utiCSwg0KxlCfiRvkcgqIpFB5wdgb1FN5LqeBSAoVF8xpErVUnMupIJJCrQVzpkQtteZy\nNoikUGfBnGxRU2su54NIChUWzNk7ozcV5nIJiKRQX8FcYFFTYy7XgEgKtRXMFXujltpyuYoC\nkaT3RdkMEZHspd4bi2kuV2lUWy7XkS9SGv7ZDxGR3KXeHYtlLpdpVFkuF4JIChUVzHW7o6aq\nXC4FkRTqKZgrNaopl2spOUYiGLHhzliscrnWo3pyuZiis3b7VvGuwexpavCwDy72qJpcrqZo\nj5R2nYe5azBDw52x2ORytUe15HI5HCMp1HEscLlHleRyPYikUEXBXO9RHbk4UPSCbCIYqeG+\nWAxycfCoilw8KDnZsPPafXcNZmy585KGh3Px8KiGXFzIFSk1uy+Deddgulb7rw56NBcXjyrI\nxYei099nDRGR85f64Ag+HsXPxYmikw1Z7XcMEZGCg+qC1ntGUPDxKH4uTmRP7fZf4fyuwbxb\nHbjw+7FcnHZI4XPxonSPdMoQESnb8p45goiXR9FzcYOTDQrBD6q9PIqeixucbFCIfVDt5lHw\nXPwoekGW907JDS9/fc1tYhc8F0d4i5BC6LfC+HkUOxdHEEkhcsE47pBC5+IJIikELhhPjyLn\n4krJMRLBiA13xnJApL0PtCBwLq5YnbXrWz3pz72XdGOZi+sOKXAuvhRM7VZfGEjjt8cE07db\nf73khFxcPQqciy9Fx0h6m5SmNp8N7xpMRrszcvHdIcXNxRkbkYZ7HxVMVrsTRNr3MCvC5uKM\n7dTuUcH07XKmMIa5OO+Qwubije3JhmUwy3cA/w4UL6MLpgfVWbnkJePskUsuNVRNkUgrjfO3\nMPEzeVNQMKtHj9PXzS1vTjLeOySPXFrCV02BSKsXKis4CxM9ko7sglm/flvR2amMavH2yCWX\nlugmFZ1skFZ71io1Wa8LRI+ko+CgeiWWklyajG2M+w7JJ5eW4FVjJVLJEMEj6TAqmLIRtkXa\nMY4tLrm0BK+akqldY3Rdu+CRdORPYXbGIoq0EY3/Dsknl5bgE5mSkw1W17WLnUhP/qpaXr9t\nS6RdA5nik0tL7LKxOv1d0kHwbUvH+a9mlIsUYIfklEtL7KopOEba+RF86ezUnm4uJvtYYG8s\ne0TaNY4tPrm0xK6aXJHmb446PETsSDry1vRALOUiRdghOeXSEnsiky1SdvuMIUIn0pNZMEWt\nt0dYLZYIHjnl8iZ02SCSglPBIBIi5Q8ROpGeeCKFmNkhkkL+MZLhBSJDJ9KTeyxgfCHENZF2\njGGPUy4tocvG4/R38MPGDqfTvHoyMXZIjqe/Een7V6Ej6Ygn0mkLUoSjSKHLBpEUEEkGkWQQ\nSSGaSEFmdoikgEgK4UQ6b0GKQCQZRFLwEkmLBpFilw0iKXgVjBJNlJkdIikgkkI0kU5ckCIQ\nScZNpNfraL/n4iiSEE2YHRIiKTiKFNukYCKdvTDZuIoU2CQ/kRpEkpC3MYj0BpE+QCR1BDGa\nODM7RFJwEukXkYpEOn1hsgkt0k+H6eJk4iNS8xYptEm+In1Eg0gdGx8g7hVykclRpNi7JD+R\nXl/RBJrZ+YuknKX6cOdylxBJwVGk5kuk05clH3+RxLoRtLnWJS+Rvje7wUAkGW+RXuIxgaLM\nhSq5idQgkvjbXqR5NpFmdv4iNULhrCR0lUuuIkU2yVekRa1E8iikSOuuXDPF8xQp9C7JTaTv\naBBpYBCp9FzMBefxEEnBrWD6aKZsQs3snEUS9te5+ZwsEyIpOIs0yyaUR74iSRPfknxOfMHW\nV6TAJiGSjLNIr6/CKffiZ8l2iw02l7rJbLGjg69qCYe7SEM2sWZ27iJ13w0nvrobRdiINF0K\nULgo4INFKs6l+a6VWB755fJGEMlieQwwEmloJjV+skiluTSzWum/R6mUHrdc3nyGE2d/bSJS\nWn7J6OBz/hIPi4Ipz6WZZXJfkXbl8ub3I5w46ViLlPlXqr930dEwLpj8v949RtIZFWaT2+OW\ny5tJpO5GnHTs90hTMisXUn+cSJm5NPNIOpEMlsMSt1ze/E6bmffXOOkYiyQ9YFWksCbZFozU\n47ZI7c04pdLhlsub3+X+Os4OyVukuLukACK1xRKoVDqiiNQE28rYn7UrOf19c5HKc2mWibxe\ngUqlwy2XN6/l/jpQOkYivWe2qbuV08Fr+NBw2Lmd4eslBbk0H5uWu4q0I5c3r9mHzWPtrm1E\nKu5gCuTWIu0Z4eMNmdHyccvlzUKkUFsZRFIIIlK4C2nGESnUDimASLHqZCSISH9LTl+oTSKJ\nFKl03EWKukuKINLfz89bnder3y8F8CmUSIFqJ4BIgdKY4SfSGM3f32zu8urp7vDTKY5I/2d2\ngUrHX6RIaczwF+m/Kd8HATOffHZPQUT6v97ttDdO7SCSgrtIrSHq0fRr5GqfYojUbmR+3jqd\nvjiZRBDp44Sv3/Wb53iL9PZoO4VRqKvOSLiK9JqHM94KQQCRpl3SzKBjH1c0wKtgxlppil+4\n/1bqg9cnZkttyZpIzVc4UUwKIdLwlnj1cpkOWrmJ1NXK3w6RttANy2NlqS3ZFmmxtw4yvYsg\nUve2qZKrKq1ehsIGb5GaYC84DoQQqZlvZF4hTPIVadjSGb1tquziL+s6+op0xg7JBm+Rfr/D\nCWGSq0hdAO3kPd7G11mk9ma4TFqcRWrT+dpbRzDJV6T+p794781EJIUoIi1/73+g5CNSd2pq\nXPmfn0CvrPW4ivQX9hApgEjitNd9p+Qk0rTZbd4VE+k16g5fkdpbIT2KIFIjbWS8TXIUadoh\nNU08kxBJxlOk1ySScG/0N/OeJNLvfIfUxDPJU6TAMztfkf7//105odm+5HzOMm3jJVIzivQz\nvkQd6n3griK1N2J6FECkRt/IdO/eOM1liycAAAnISURBVGWptggg0vBLvxAkEEkmhkhqk/GN\nUGof9mws9fZ6ZaKI1M9oZxuX/7+Jo5K3SEFndtFFGg6VNsreHkeRupvzTP46lWbre8L2Yxub\n1d5EHOHVHj6GfVtD4y9SxrviXc46eIv0kckYwcnbj228RGoivz+oCSDSb0Y2Dp8fdhfp4x7/\nl6h7PEVq4s7sIoiUlU3WJ7OOvh++6F3xJ4nUjS9kEkQlb5FOH34njiK9BpFye8p0wAQXkdpE\n5B3SG4cd8zeIJOMpUtPWRoFIV+IjUjO8eK9lYry52IGfSJkfMnfCW6TmN2Y2ziKtPc5mArsD\nm9XeRBUp53DajQAinb4Ee/ASKfK7YN4gkoyvSP3VgwLiJlL79u+YkXT4FMzrlX9eygc/kfoj\n659f5X5fPEUKXC5uIgU+nO5wFKlpRfpBpJFusxu6XDxFino43eEq0vvis4g00CUSulxcRQq9\nhXETqZvZ/S+a35Am2Yg0/i3qrD/33ldL5HIxKpiyXN7FEvdwusMll5Zhh9Tk75J+r2FtqbfX\n67NNWt5Y6aCrltCb3caqYEpyacZNTMhNbo9LLi3DDqlp1vZJQolfg4lIafgy3ljroIodkknB\nlOUy7ZBiTl46PHJpabMZLwG6tXtw4HqRaqiWxqVgpondWkGcOUlZxyuXLpu/v9Db3jNFSj3d\nT1+XM/WriA1yV/uUXPpq2Vo8N/xyCe2Rx9SuDrymMNEhFxlEUqBgZMhFxuGsXR24nZ0KDrnI\nGIn0ntmmpuR1geDYFAy5yJ08MZfD63XXYOKPcAbkIoNIChSMDLnIIJICBSNDLjJXiFQnR1eb\nXMhluV5W+VzR4pJBbCGXvcNVlgsinQu57B2uslwQ6VzIZe9wleWCSOdCLnuHqywXRDoXctk7\nXGW5INK5kMve4SrLBZHOhVz2DldZLoh0LuSyd7jKcrk4P4B7gkgABiASgAGIBGAAIgEYgEgA\nBiASgAGIBGAAIgEYgEgABtiItPlh3LTRqLtrs8XGOFst3ndmfXDYCHJZGXC1wUajeLkYibTV\nUVpvNN2ttnjftz5OWu+iv2NzWQ0hF3XEu+Vikl3a6CmltN6oC2azxcY4ab2L2Z0XVQy5bA2o\n3F9hLpeINGm9ll5GvqsN0lYX3Z56c1ntIBdtibbGqi+XMCJttsjYfmS02NpKmUIu6oC3yyWK\nSCmjm/WV2g5muw9jyEUZbHOs+nIJIlJafpEbbAQzXsjvRgVDLls9RMnFJrq01VFab5Q2u+k6\n2BgnrS9JGvb3F9ULuawv0HqDynIxEung6wL9xuGC1wUyltUQclEHvFsul2UHcGcQCcAARAIw\nAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQ\nCcAARAIwAJEADEAkAAMQCcAARAIwIJxI6crrq1UEuchEySXCMizor7EJH5CLTJRcIizDguFa\ny+PfyElNjE2OM+QiEyWXcE9FH8x08ecUZZvjC7nIRMkl3DORxi/zYIBcZKLkEu7ZmO2qU7+p\nYQrTkItGlFzCPRWfu+o0/fbRkItMlFzCPRPTZLf/2x5hzsv4Qi4yUXIJ90yMf0YtTedhmMKQ\ni0aUXKI/FdGXzwtykXHLJfoTEn35vCAXGUSSCb54bpCLjF8uPCMABiASgAGIBGAAIgEYgEgA\nBiASgAGIBGAAIgEYgEgABiASgAGIBGAAIgEYgEgABiASgAGIBGAAIgEYgEgABiASgAGIBGAA\nIgEYgEgABiASgAGIBGAAIgEYgEgABiASgAGIBGAAIgEYgEgABiASgAGIBGDADURKAC2+Veg6\nugk3WAUwAJEOcoNVAAMQ6SA3WAUwAJEOcoNVAAMQ6SA3WAUwAJEOcoNVAAMQ6SA3WAUwAJEO\ncoNVAAMQ6SA3WAWd4ZXGvJVMXzceBCId5AaroJMv0dB6eeNBINJBbrAKOqn/kpp+z9R9mX/t\nG6Y0/u9udI/7fNR9QaSD3GAVdEaRul1T6r/Mvw5N3j/1e7Dh39ejbgwiHeQGq6CzEGkmQ2qW\nv0jN4utcpJl4twaRDnKDVdBZTu36dzmnaf42vu25+1Uz3Df+WzzKd11OBpEOcoNV0JlEasbj\nodl+KS3bfk3tPhrdP6qHjm7CDVZBZ5CGY6RNEOkgN1gFnXHmNj9VN5+qDRO2/o7xvtn8b3qU\n33pcACId5AarAAYg0kFusApgACId5AarAAYg0kFusApgACId5AarAAYg0kFusApgACId5Aar\nAAYg0kFusApgACId5AarAAYg0kFusApgACId5AarAAYg0kFusApgACId5AarAAYg0kF2rcL2\nnwGZPmF6Qkh5f4fkxCdn7PrzI01XVUTxH2LxTCuDh4o0fsRno9vTRDJsdWwBkvzr05nFW/CA\nQw1O5cEijR+LW3xwp+n3Vt2FeKYPbtsGNX12dfbBIWEpTiPNVnO6yNDnSp/4Iabp44iLj1mN\nn0ycPpKYGwYiHWSXSOOupn/q0nAzjfXdFXkzfPDUNKk0fBuGm1k0LtXZIi3+TSs9fEJwzOas\nBfgY/GPc8YnIDQORDnLgGGk6Dpo2kMO3WSEVXOu0YPhhczxNM+eXCjqvhJt+kPn6z44Gx0/k\nTr89ZQHGZ+DjKZh2R4vngqnd6cxX4VdHeuA0fekLeTa7WIiUfWj80pGWeaiQ+efJp9H2PDl/\nKh/rPtsdTGv3KVL51O5H53MBPgYfx02zJ6IgDEQ6yK6pXfdlnLEtL1o63TfbI5kyjPSxRxoX\nZVLsJD72SPOV/srmrAX4GHwu0nwhm8wwEOkgu0T6lGU2P29mt4ebyTipsVLG7f649V8cnZ1G\n+vw3rbTHMdIUQTMdGQ2jZ4aBSAfZtQpfp8o+LtIzCDSdtbMVaZyyyGftxl+cxbJYlwf73wt0\nzgJ8Dz7N4cZ9dX4YiHSQs+ZdV3GDpyAEiHQQRIIWRDoIIkELIh3kBqsABiDSQW6wCmAAIh3k\nBqsABiDSQW6wCmAAIh3kBqsABiDSQW6wCmAAIh3kBqsABiDSQW6wCmAAIh3kBqsABiDSQW6w\nCmAAIh0kAbT4VqHr6AA3AZEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQ\nCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEA\nDEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAA\nRAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIwAJEADEAkAAMQCcAARAIw4B8GLqGm\nevtwmAAAAABJRU5ErkJggg==",
      "text/plain": [
       "Plot with title “Densities by disease status”"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot densities of some patches by disease status ----------------------------#\n",
    "plot.pretty(out, nr_patches, \"panels\")"
   ]
  }
 ],
 "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
}