{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program 4.4 Multi-host SIR model (R using deSolve)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Author*: Timothy M Pollington @t-pollington\n",
    "\n",
    "*Date*: 2018-10-03"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(deSolve) #lsoda()\n",
    "library(reshape2) #melt()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "checkSize = function(param, L, W){\n",
    "  param_name = deparse(substitute(param))\n",
    "  param_size = dim(param)\n",
    "  if(is.null(param_size)){\n",
    "    print(paste0(\"Warning: \", param_name,\" is a scalar value, expanding to size \", L, \"x\", W, \".\"))\n",
    "    param = param*matrix(1,L,W)\n",
    "    param_size = dim(param)\n",
    "  }\n",
    "  else if(param_size[1]== W & param_size[2]==L & W!=L){\n",
    "    print(paste0(\"Warning: \", param_name,\" was given in reverse dimension order, so transposing it before use...\"))\n",
    "    param = t(param)\n",
    "    param_size = dim(param)\n",
    "  }\n",
    "  else if(param_size[1]!=L | param_size[2]!=W){\n",
    "    print(paste0(\"Error: Parameter \",param_name,\" is of size \",param_size[1],\"x\",param_size[2],\" and not \",L, \"x\",W))\n",
    "    stop(\"See above message\")\n",
    "  }\n",
    "return(param)\n",
    "}\n",
    "\n",
    "checkGreaterOrEqual = function(param, bound.l){\n",
    "  if (sum(param<bound.l)>0) {\n",
    "    param_name = deparse(substitute(param))\n",
    "    print(paste0(\"Error: At least one of the values of \",param_name,\" is less than \",bound.l))\n",
    "    stop(\"See message above\")\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sirODE = function(times,init,params){\n",
    "  with(as.list(c(params,init)), {\n",
    "    dXh = nu1-r*(Tr12*Ym+Tr11*Yh)*Xh-mu1*Xh\n",
    "    dXm = nu2-r*(Tr22*Ym+Tr21*Yh)*Xm-mu2*Xm\n",
    "    dYh = r*(Tr12*Ym+Tr11*Yh)*Xh-mu1*Yh-gamma1*Yh\n",
    "    dYm = r*(Tr22*Ym+Tr21*Yh)*Xm-mu2*Ym-gamma2*Ym\n",
    "    return(list(c(dXh,dXm,dYh,dYm)))\n",
    "  })\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "program4_4 = function (r=0.5/1e3,Tr=matrix(c(0,0.5,0.8,0),2,2,byrow = TRUE),gamma=c(0.033,0),mu=c(5.5e-5,0.143),nu=c(5.5e-2,1.443e3),X0=c(1e3,1e4),Y0=c(1,1),max_time=1000){\n",
    "  # Function contains default arguments if run as program4_4()\n",
    "    \n",
    "  # Check lengths/dimensions of inputs\n",
    "  r = ifelse(length(r)==1,r,stop(\"r is not a scalar.\"))\n",
    "  max_time = ifelse(length(max_time)==1,max_time,stop(\"max_time is not a scalar.\"))\n",
    "  Tr = checkSize(Tr,2,2)\n",
    "  gamma = checkSize(gamma,1,2)\n",
    "  mu = checkSize(mu,1,2)\n",
    "  nu = checkSize(nu,1,2)\n",
    "  X0 = checkSize(X0,1,2)\n",
    "  Y0 = checkSize(Y0,1,2)\n",
    "  \n",
    "  # Check parameter values are valid\n",
    "  checkGreaterOrEqual(r,0)\n",
    "  checkGreaterOrEqual(max_time,0)\n",
    "  checkGreaterOrEqual(Tr,0)\n",
    "  checkGreaterOrEqual(gamma,0)\n",
    "  checkGreaterOrEqual(mu,0)\n",
    "  checkGreaterOrEqual(nu,0)\n",
    "  checkGreaterOrEqual(X0,0)\n",
    "  checkGreaterOrEqual(Y0,0)\n",
    "  \n",
    "  if (Tr[1,1]!=0 | Tr[2,2]!=0) {\n",
    "    print(\"Warning: Transmission probability between human-human or mosquito-mosquito species (i.e. T's diagonal) is non-zero'\") # We print this warning as it is unusual to have non-zero T_HH & T_MM for basic vector models.\n",
    "  }\n",
    "  if (Tr[1,1]<0 | Tr[1,1]>1) {\n",
    "    stop(\"Human-to-human transmission probability is not between zero and 1.\")\n",
    "  }\n",
    "  if (Tr[2,2]<0 | Tr[2,2]>1) {\n",
    "    stop(\"Mosquito-to-mosquito transmission probability is not between zero and 1.\")\n",
    "  }\n",
    "  if (Tr[1,2]<0 | Tr[1,2]>1) {\n",
    "    stop(\"Human-to-mosquito transmission probability is not between zero and 1.\")\n",
    "  }\n",
    "  if (Tr[2,1]<0 | Tr[2,1]>1) {\n",
    "    stop(\"Mosquito-to-human transmission probability is not between zero and 1.\")\n",
    "  }\n",
    "\n",
    "  # ODE solver\n",
    "  params = c(nu[1],nu[2],r,Tr[1,1],Tr[1,2],Tr[2,1],Tr[2,2],mu[1],mu[2],gamma[1],gamma[2])\n",
    "  names(params) = c(\"nu1\",\"nu2\",\"r\",\"Tr11\",\"Tr12\",\"Tr21\",\"Tr22\",\"mu1\",\"mu2\",\"gamma1\",\"gamma2\")\n",
    "  init = c(X0[1],X0[2],Y0[1],Y0[2]) \n",
    "  names(init) = c(\"Xh\",\"Xm\",\"Yh\",\"Ym\")\n",
    "  times = seq(0,max_time,length.out=(max_time+1))\n",
    "  sir_out = lsoda(init,times,sirODE,params,rtol = 1e-5) # rtol matches the original MATLAB code   \n",
    "  return(sir_out)\n",
    "}    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Warning: gamma is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: mu is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: nu is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: X0 is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: Y0 is a scalar value, expanding to size 1x2.\"\n"
     ]
    }
   ],
   "source": [
    "sir_out = program4_4() #Reproduces Keeling & Rohani's results as seen in their original MATLAB code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"Warning: gamma is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: mu is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: nu is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: X0 is a scalar value, expanding to size 1x2.\"\n",
      "[1] \"Warning: Y0 is a scalar value, expanding to size 1x2.\"\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   user  system elapsed \n",
       "  0.278   0.000   0.278 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tic = proc.time()\n",
    "invisible(program4_4()) # invisible() hides function output, apart from warnings\n",
    "proc.time() - tic # ~0.1 secs on Dell Precision M2800 laptop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ggplot2)\n",
    "library(gridExtra)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAFVBMVEUAAAAA/wAzMzNNTU3r\n6+v/AAD////hNSoYAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO2djXbbuA4GkzrN\n+z/yvdv8kKIICgRAiZQG52yzTegRjE8Ty7Ktvn1SFOWut6sboKg7FCJRVEAhEkUFFCJRVEAh\nEkUFFCJRVEAhEkUFFCJRVED5RPrY1u4bUinXaXnBOCVv9L0NSjgiMdWqJ7LSYBHJzEMkWGmw\niGTmIRKsNFhEMvMQCVYaLCKZeYgEKw0Wkcw8RIKVBotIZh4iwUqDRSQzD5FgpcEikpmHSLDS\nYBHJzEMkWGmwiGTmIRKsNFhEMvMQCVYaLCKZeYgEKw02VqT3/1dY04jUUa+vP18v8ctAkd6D\nStdX1NZMdaJIGpMQScHr8uj1bdNL+lJPTBeXdA8u3aWvqTNFUkSDSApel0ffyvz3R/VLPTFl\nXLV7cOHOfGUNEqmon41Rp5dJpO+8+kW6aje+vs55RPreWH8w0jIekSJF+vP/2tzod+/oqot2\n4TlKnAoimXnri2Q5hLhqD56kxLkgkpm3nEj7xLSBpZtdtQPPUhzadeAQSboH8g5W6f+AtauO\n3fmQtdngeq8jIZIDF/A6Uvfp706RGrt0pf8262DRIixECuDNJ1L3C7J9IskOVftvso4WLcJC\npADeVCLZEusSqWlRpf8W63DRIqw0WEQy8x4m0oFG6+z8kaw0WEQy854l0pFG6+z8kaw0WEQy\n8x4l0rFHy+z8kaw0WEQy854pkr5/kaVZtAgrDRaRzLwniaTwaJmdP5KVBotIZt6DRNJ4tMzO\nH8lKg0UkM+85Iqk8Wmbnj2SlwSKSmfc8kfr6l5apFi3CSoNFJDPvMSLpHpCW2fkjWWmwiGTm\nPU6kzv6lZapFi7DSYEd81LznbfldbZ+EQ6TaIkSqsNJgEcnMe5pIvf1Ly1SLFmGlwSKSmfcQ\nkbQeLbPzR7LSYBHJzEOkdv/SMtWiRVhpsIhk5j1DJLVHy+z8kaw0WEQy8xCp3b+0TLVoEVYa\nLCKZeY8SydC/tEy1aBFWGiwimXmPEEn/gLTMzh/JSoNFJDPvPiI1AkOkJisNFpHMvBuIdPiQ\n1OHRMjt/JCsNFpHMPERq9y8tuxMrDRaRzLwHiXT5DjsrKw0Wkcy8B4j0jkhtVhosIpl5iNTu\nX1p2J1YaLCKZec8R6foddlZWGiwimXmI1O5fWnYnVhosIpl59xfpHZEOWGmwiGTmIVK7f2nZ\nnVhpsIhk5j1GpAl22FlZabCIZObdXqR3RDpipcEikpmHSO3+pWV3YqXBIpKZ9xSRtPdglZ0/\nkpUGi0hmHiK1+5eW3YmVBotIZt7dRXpHpENWGiwimXkPEUl9D1bZ+SNZabCIZOYhUrt/admd\nWGmwiGTmIVK7f2nZnVhpsIhk5t1cpHdEOmalwSKSmfcMkfT3YJWdP5KVBotIZh4itfuXlt2J\nlQaLSGYeIrX7l5bdiZUGi0hmHiK1+5eW3YmVBotIZt69RXpHJAUrDRaRzDxEavcvLbsTKw0W\nkcy8R4jUcQ9W2fkjWWmwiGTmIVK7f2nZnVhpsIhk5iFSu39p2Z1YabCIZObdWqR3RNKw0mAR\nycx7gkg992CVnT+SlQaLSGYeIrX7l5bdiZUGi0hmHiK1+5eW3YmVBotIZh4itfuXlt2JlQaL\nSGbenUV6RyQVKw0Wkcy8B4jUdQ9W2fkjWWmwiGTmIVK7f2nZnVhpsIhk5l0l0tvbf/8Zgyu2\niUg+VpaKLY96LIjkwGlF+s+htzerScU2EcnHymIxxSHEgkgOnF6k9J87MUTysbJYLGlIsSCS\nA4dIC7KyWBTDf339+XplXxDpziK9I5KO1SXSt0BfQr1+vEKk6042OE43FNtsi9R3D1bZ+SNZ\nPSJ9PQR9y/TjFCJdJ9Lnm/0BCZFiWT0i/T4SFSL9+X+VK78DMEVMnVHFjoFIPlYaLCI9q4od\nA5F8rDRYu0iVWDi0c+A4tFuQhUgBvKtPNowTqTxpd/0OOysLkQJ4V57+/hx6+nv33ct32FlZ\n/SJx+ns0DpEWZBlE4gXZwThEWpDVJ5I6FkRy4DpFGvmCLCJpWVksljSkWBDJges62TD23d+I\npGVlqZQhdf2uKzeDSHZc1+nvoZ9H2p20u36HnZWVhVJG1PVZl3IziGTH6UXyVLHNlki992CV\nnT+SlamzGXPvO4vLzSCSHYdIC7IydTZjRqQO3nUi2Y/sECmWlWVSJIRIat5tTzYgkpolidT7\nWZdyM4hkx+lPf79lX3qr2CYi+VhZLLWYONlwAa7ndaT8a18V20QkHyuLxZKGFAsiOXDTiLQ/\n+339DjsrK4vFkoYUCyI5cLOJ1H0PVtn5I1lZLGVKXc9ky80gkh03zXMkRNKzsljqHiHS6bhp\nztohkp4li+SKBZEcuGleR0IkPSvLpEjIFQsiOXDTvLMBkfSsNNjyOZIrFkRy4HpPNtiq2CYi\n+VhZLGVMnlgQyYHjrN2CLEkcTjZ08K4Q6S0vRLqchUgBvKtFsnikEOkdkfQsSSRnLIjkwM3y\nHAmROlhZLJ5UECkQN5lI/fdglZ0/kpXFUqbEod3kIvmq2CYi+ViZOZsx8xypg3eBSN+fcRl6\nsgGROliySF9HDnyM4nwcIi3IEkX6+Q6PSKfjOLRbkJUGO+Sj5ocmIZKCF6iPmBgi+ViZOpsx\n/16OyyiS9iEJkRQ8YeaDD+2qvwsv32FnZUki/R6E22JBJAcOkRZkiSL5YkEkB67z0G7UlVar\nCV6+w87KcudRjwWRHLje50iD3iKESD0sMQ7n60iIZMd1izT00M5wD1bZ+SNZWRxFOIik5iFS\nu39p2Z1YWRyucHabQSQzbq6TDYZ7sMrOH8nKYilScsWCSA5ct0h9UdUTQyQfSzbH9wlZRLLj\n5nhBtnr2+/oddlZWGiwnG8w8RGr3Ly27EyszZzNmTjZ08K4+tDMd3xXblESy3INVdv5IliyS\nKxZEcuAQaUFWFkuRkisWRHLg1Id2P5cs7ouqnhgi+ViyOZxsUPMue0Tafu2rYpuI5GNJ4vAc\nqYOHSO3+pWV3YmWxFCEhkpqHSO3+pWV3YmWxWNKQYkEkB67jZEP2xZkYIvlYWSyWNKRYEMmB\n6zjZMPCdDYjUxcpSqWbEod35uI4XZM0aIVIwK8tklxAiTS+So4ptIpKPlalTm7b6191uM4hk\nxk0h0jsidbHSYOvOcM2G83FdL8iaj+2KbSKSj3WkDId25+N6ztoNO9kg5Hf5DjsrK4tFCMsU\nCyI5cD2vI/Vce7CZGCL5WFksRUicbFDzEKndv7TsTqxDkWyxIJIDh0gLsiSRnLEgkgPX9YLs\nqHc2IFIfK0vFkoYUCyI5cPrT32/2ByREimVloewy4jnS7CJ5qtgmIvlYmTmbMXOyoYOHSO3+\npWV3YmXqbMbMR807eDd8r907IvWxskyKhFyxIJID13WyYdALsojUycpSKUJyxYJIDlzPOxs+\nB521k+K7fIedlSWbwzUb1LwLX0fKv/ZVsU1E8rFq4lgu87TbDCKZcYi0IAuRAniI1O5fWnYn\nVk2kgFgQyYHjOdKCrCwWSxpSLIjkwE101s52D1bZ+SNZWSqmOIRYEMmBm+h1JNs9WGXnj2Rl\nmRgDqcaCSA7cDO9sQKROVhqsT6RdfScRC6Xiqtgxtua8I1InKw2WRyQz7+rnSPEnGxCpl5Wl\nYklDigWRHLiOs3aDRTLeg1V2/khWFss+po6XKHabQSQzrvd1JFsV20QkHyuLpUwJkSYXyXcQ\nUWwTkXwsKZV/Fn1yFaErcPqzdh6Tim0iko8lhYJIHbzbnmww3oNVdv5IFiIF8DjZ0O5fWnYn\n1oFIvGn1AhwnGxZkZbGUKXFdu8lF4mTDRCw5lb6Dht1mEMmM42TDgqyYUBApEnf9c6R3ROpl\nZbFY0pBiuVSk9015cA8XyXoPVtn5I1lZLJWAVvuE7LtUNtzkIvmq2CYi+VhpsKuLJEpUkQmR\nECmYlQa79qHdoUabbhBpn1g+IkTqZqXBLiySzqKsIURCpGBWGuxOpK5nsbvNnCdSj0bfPSES\nIgWz0mAXfUG2W6N/hUiIFMvKzNmMufdST7vNnCOS+JjT/rmmM117B8L24r6XIdJyrDTY+kX0\npxZJ64ddJfVm9WW4t+LU+z4z1k4MkVysLJQiou3X3lhOEKm2j8o422694fkE6imlSJ0fvmwn\nlg9EHs7lO+ysrCyVIqPt195Yhou02/mOccJee9xesCbHpROp96Mu7cTyaciTuXyHnZWVxbJP\naeKPmpe7nhIn7rp9q0cXIi3HymLJApr9nQ3FfteDO8cEZyHScqwlRdrudd24k2xwVJdIMW9a\nRSQfqyZSRCwjRdrsczZc4E5fdif02VUdJxvMV/8WRpJ3XpucaryL7PyRrCwVUxxSLONEKvdi\nK86+m++2q39BNlakzg9fthPL7pk03Y8JdthZWVkoWTz9n3XZbWaUSE2POnFegfbdHW23q73s\n7wOq3OZeJPM9WGXnj2Slwa4h0sFe7XyA6zVojzvYbnd7338fUOU2EcnFSoNd4tDuaNcOOQmo\nFkjANRaa2quJZDkh1E4MkVysLJoiKF8sY0Q63MOtL0uJC6/AIdKCrCyaIqjt1ylEOn6keIxI\n/iq3iUguVhpsZsycryNpDrgQSV3lNhHJxcrsyf63/6hht5lwkVRPXJ4mkv3sNyLFsrJMioR8\nsYSLpDsB8CyRRr0gi0j9rCyVIiNfLNEi6Tx6lki9H75sJ4ZILlYWiyUNMZZgkZQePU2k7Vdf\nYojkYmWx7HOa5uInSo0QqaPKbSKSi5XFUqY00ckGrUeIpK9ym2nCrWFfvsPOyspiKUKa6OIn\nao+eJVLvhy/bie1Est+DVXb+SFaWSiWjOV6Q1Xv0HJFGvrMBkQysBUTq8AiREOki1vwi9Xj0\nHJH8VW4TkVysNNhJL37S5REi6avcJiK5WGmwu5MNXYcOu80EidTn0XNEsnxmDJEQSXuxbkRC\npGtYkkjeWIJEekekDw7tVmClwc4oUvIIkfZlfuN3NTFEcrGyWHYx9VxXereZCJHSEyREqtTg\ndzbY78EqO38kK4ulTOlykbITDYhUFq8jzcWSROq9HO5uM4hkxnWLZPEIkWJZE4uUeYRIdZdM\nBgmJ/Q67mdzlO+ysrCyWXUjXipR7hEi1QqSJWFksu5B6LtC+24xXpI1HiBRf5TYRycXK1NnO\n+erPI208QiShRlz8BJEsrCyTSkb6lHabcYr0jkj53+s15mQDIllYWSqmOKRYokQS2u7F7drT\nrZtbpEEXP0EkCyuLxZKGGItTpHdE2vy9WoNekEUkCyuLpZbSRc+RSo8QqVaINBEri6UI6cpr\nNpQeIVKtEGkiVhaLK6TdZjwi7R6QEKlaYy5+gkgWVpZKJaNrRNp7hEi18l23odwmIrlYM4vU\naLsHV1uGSIgUy5JFuuw5UuUBCZHiq9wmIrlYabDTvLOh4hEixVe5TURysTJzykl3HS/sNmMW\nqeYRItWr7zNj7cQQycXKQumO4tWIBZHsuK63CCHSHKwsld4kXkNEqnqESLXq/ajLtsptIpKL\nlcVSSakV0WuISLUzDZW2tThpGSIhUixLEuk/h45ONowU6aBtLU5ahkiiSO3gLt9hZ2XJIqX/\njkT68//a/fA7j850fx+QOm/3yOr98OW2yh0DkVysLJZ9SNeJ1Hmzh9aQzyMh0mUi1WKxHdoJ\nz5A4tGuaZPIIkWJZkki/p1YvEemwbR1OXnYTkTxVbhORXKzMnGLQx69RhIskeoRI8VVuE5Fc\nrEyc7igQaQyu4wXZ9Gd3ldtEJBcrS8UUhxSLRSTZI0Sq1ZhrNiBS9OtIXaeEdptBJDOu53Wk\n/GtfldtEJBcri6UI6XSRGh4hUq0QaSJWFks1K2ssiGTHIdKCrCNnzvs80jsiSTxh5jxHmoh1\npMx5h3YtjxBJCId3NszCylIRsrLF0i1S8wEJkYQacO1vRFr7ZEPTI0SKr3Kb75ty3YNVdv5I\nVqbOZsy9F6fZbaZTpIMEESm8ym0ikouVqRMai1UkZdtHuKNl9xDp8N2QPYkhkouVpWLLQ4il\nU6SjABGpUpoPX+oTQyQXK4tlE9F3UGe9jnSYHyLtS/NRF0Q6iVUV6UugE58jHeeHSIg0NUsQ\nKc/pRJHUbbdxx8sQCZFiWTWR8nfnn3P6G5GaPGnoig9fItJJrClEOvQIkUSToj+PhEjBIn3H\ng0jn43gdaUFWGux1Ih17hEjxVW4TkVysNFhEMvOuex3J8a7VcpuI5GLVRPr890LfT1aIdDYO\nkRZkCSJ1f9Rltxm9SAqPEKlRvEVoBlY9j7fOAztEisT1PkfirN0ELG8cUix6kTQeIVKrBhza\n+e7BKjt/JCuLw5KGGMvnwQfEUtOIdMg7mD0iTcDK4rCkIcaiFknlESLVauDJBt89WGXnj2Rl\nsVjSEGNBJDuuW6SQxBDJxcpiMcUhxaIVSecRIsVXuU1EcrHSYBHJzEOkdv9PYKXBIpKZd4lI\n/R++bCeGSC5WFowtDyEWpUhKjxBpV4YPX7YTQyQXK0vGFIcUCyLZcUqRfv5EpBlYWTKmOKRY\n+kTqbltcp1y2vEiWz4y1E0MkFyuLxpKGGItOJO0DEiKVhUiTsbJoLGmIsSCSHacWqfsdka3E\nEMnFyqKxpCHGohJJ7REilYVIk7GyaCxpiLEgkh2HSAuysmgsaYix9IhkaFtcp1y2vEiWD1+2\nE0MkFytLxpKGGItGJP0DEiLtyvDhy3ZiiORiZcmY4pBiQSQ7TvvOBs+BXVsk5z1YZeePZGXB\n2PIQYlGI1OERIsVXuU1EcrHSYBHJzEOkdv9PYKXBni1Sj0eIFF/lNhHJxUqDRSQzD5Ha/T+B\nlQZ7kUiqu4ZI8VVuE5FcrDTYk0XqekBCpPgqt4lILlYaLCKZeYjU7v8JrDTYc0X6iSx4z0ck\nc2KI5GKlwSKSmYdI7f6fwEqDvUSk6D0fkcyJIZKLlQZ7qkjviNTBcyqjSwyRXKw0WEQy8xCp\n3f8TWGmwZ4qUAkMkBc+pjC4xRHKx0mARycxDpHb/T2ClwSKSmYdI7f6fwEqDPVGkLC9EUvCc\nyugSQyQXKw0Wkcw8RGr3/wRWGux5IuVxIZKC51RGlxgiuVhpsIhk5iFSu/8nsNJgTxNpkxYi\nKXhOZXSJIZKLlQaLSGYeIrX7fwIrDRaRzLy7ieS9B6vs/JGsNNizRNqmhUgKnlMZXWKI5GKl\nwSKSmYdI7f6fwEqDPVmknqYRaUCV20QkFysN9iSR3hGpm+dURpcYIrlYabCIZOYhUrv/J7DS\nYH0i7es7lPp3y29Tp1e5YyCSi5UGe84jUpkVj0gKnlMZXWKI5GKlwZ4i0i4qRFLwnMroEkMk\nFysNFpHMPERq9/8EVhrsGSLtk0IkBc+pjDIxRPKw0mARycxDpHb/T2ClwZ4gUiUoRFLwnMoo\nE0MkDysNFpHMPERq9/8EVhrseJFqOSGSgudURpkYInlYabCIZOYhUrv/J7DSYIeLVI0JkRQ8\npzLKxDQeXb/DzspKg0UkMw+R2v0/gZUGe5ZIhqYRaUDttolIHlYa7GiR6jEhkoLnVEaZGCJ5\nWGmwg0USYkIkBc+pjDIxRPKw0mBPEsnSNCINqN02EcnDSoMdK5KUEiIpeE5llIkhkoeVBnuO\nSKamEWlA7baJSB5WGuxQkcSQEEnBcyqjTAyRPKw02FNEsjWNSANqt01E8rDSYEeKJEeESAqe\nUxllYojkYaXBDhSpEREiKXhOZZSJIZKHlQZ7gkjWphFpQO22iUgeVhrsOJFaCSGSgudURpkY\nInlYabDjRTI3jUgDardNRPKw0mCHidQMCJEUPKcyysQ0Hl2/w87KSoMdJVL7Fx0iKXhOZZSJ\nIZKHlQY7WiRH04g0oHbbRCQPKw12kEgHR96IpOA5lVEmhkgeVhrsYJE8TSPSgNptE5E8rDTY\nsSK5mkakAbXbJiJ5WGmwQ0XyNY1IA2q3TUTysNJgEcnMQ6R2/09gpcGOFMnZNCINqN02EcnD\nSoONFulD4xEiaXhOZZSJIZKHlQY7TiR304g0oHbbRCQPKw12mEj+phFpQO23iUgOVhpsuEgf\nilwQScNzKqNOLCqvVXb+SFYabLxIH8e5IJKG51RGn1hQXqvs/JGsNNgBIkU1jUgDynYfLt9h\nZ2WlwSKSmYdIsNJgEcnMQyRYabCIZOYhEqw0WEQy8xAJVhosIpl5iAQrDRaRzDxEgpUGi0hm\nHiLBSoNFJDMPkWClwSKSmYdIsNJgEcnMQyRYabCIZOYhEqw0WEQy8xAJVhqsT6Si/vyJpIXz\n5saF887e5rNZiDQLDpGWZiHSLDhEWpqFSLPgEGlpVqhIFPXUQiSKCihEoqiAQiSKCihEoqiA\nChXp9XoFcl6v369eWs4Jwb3C2nv9QPdfxlfIhgLTiowqMCdNRpEivX42GcBJLDd1y4lpMqy9\n71Q2vUXN8XjjERuKTCs8qpDGVBkFivT6/SMA9Nuon7rhxDT5imrv+5flprewOR5uPGRDkWlF\nRxWSky6jCUX66fv7kdlN3XACd56Y9lYXKTSt6KiCclpUpH+Y188xQ8gj0it2L/19qA/grS7S\nP0ZUWsFRReW0rki/T+TiHkLCRQri3UCk2LSiRQqALSvSK/s/RBpWQRsKTguR4o4UclzInhq7\nl76yPxApMK3gOYTldLJIn69Ij34Obv3ULSegyd9dJ6K9n0heuy/jK2JDoWnFRhWWkyajCV+Q\nzV9Ji6C+2q+kGYA1rBu2/zK+AjYUmlZsVGE5aTLiLUIUFVCIRFEBhUgUFVCIRFEBhUgUFVCI\nRFEBhUgUFVCIRFEBhUgUFVD3Fentt/77/6u7oRp1h6gWbVtRd0jnIXWHqBZtW1mrpvLAWj2q\nxds/qNXTeVCtHtXi7R/UTzo/xwz//v79za8jCWqSWj2q6Rt0VZHOv0C+U/k5JqfmqNWjmr0/\nX23TSX/8/GT6eJ5Tq0c1eXvOKo8X0h9rpPOgWj2qydtzViudJQ4YnlOrRzV5e846+jVHTVOr\nR7VEk+ZaPZ0H1epRLdGkuZrHC5+rZPSIWj2qydtzlpzO5xIH3g+q1aOavD1nNdJZ4lW+B9Xq\nUU3fIEWtUIhEUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhE\nUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhEUQGFSBQVUIhE\nUQEVItLHtnbfEEq3LJQ2cWs7WkQyZHRWRohkhyHSENrErSHSENhSIlFjC5HssKVEOrV1F23i\n1nhEGgJDpCG0iVtzivT6+vP12n8hpHG0OInIaBQtzVch0rc5X0JtvxDSQBoiBcAmEunrsedb\npu0XQhpJC/WIjIbQ0ny1h3aIdDot0CIyGkRL83WL9Of/FRs5NaJO3b9ctIlbQyTq1P3LRZu4\nNQ7thsA4tBtCm7g1RBoCQ6QhtIlbi3gdqe/099//16DeXbQ1Q4qoysb+KlIioyYtzXfMC7J/\nNRkRUpsWJ1E1oy+PjlIioyYtzXfIW4T+qjIipDYtIhk5o9+Q2imRUZOW5otIdthNRGrGREZN\nWpovItlhdxGplRMZNWlpvohkhyGS1LqLtmZGiGSH3UakRlBk1KSl+Q4ViTNCHlpEMoqM2kGR\nUZOW5otIdtiNRBKTIqMmLc0XkewwRJJad9HWzAiR7LA7iSRFRUZNWprvWJF4sc9Bi0jmMKMP\nRPLQ0nwRyQ67lUhCVGTUpKX5IpIddguRDh6SyKhJS/NFJDvsXiLVsyKjJi3NF5HssHuI1H5I\nIqMmLc0Xkeywm4jUNImMmrQ038Ei8c5iOy0imcOM/vv/Vlhk1KSl+SKSHXYXkVoPSWTUpKX5\njhDpLyJF0CKSETNCpBBami8i2WG3EalhEhk1aWm+iGSHIZLUurjudNhtROJDY2ZaRDJiRluR\nZJPIqElL80UkOwyRpNbFdafDEEns3UVbM6QBGRUiiSaRUZOW5otIdhgiSa2L606HIZLYu4u2\nZkgDMipFkkwioyYtzXe4SFxYA5H6aGtmhEh2GCJJrYvrTochkti7i7ZmSAMy2okkmERGTVqa\nLyLZYYgktS6uOx12I5G4ZhoiddHWzAiR7LB7iVQ3iYyatDTfgSId/gM8hNSkRSQjZoRIIbQ0\nX0Syw24mUtUkMmrS0nwRyQ5bWaRaNojUTUvzPUEkLoe7iEi1b5JRk5bmO1Kko4ckQmrSIpI5\nzqjyzfy7ZNSkpfkikh2GSFLr4rrTYYgk9u6irRnSyIxq382+TUZNWpovItlhiCS1Lq47HYZI\nYu8u2pohjczo4Ntk1KSl+Z4hEv9kiIkWkcxxRgffJqMmLc13qEgHD0mE1KRFJKPIqPbt9H0y\natLSfBHJDkMkqXVx3ekwRBJ7F9edDruhSDuTyKhJS/M9RST+ESsLLSIZRUbV7yOSjpbmi0h2\n2B1FKn9ARk1amu9YkdrHdoTUpEUko8mo+oOfn5BRk5bmi0h2GCJJrYvrTochkti7uO502C1F\n+kCkDlqa7zki8e+TGmgRyWgyqv/kr6d1cd3psLuI1HxIIqQmLSIZVUbVnyCShpbmi0h22D1F\n+kAkPS3NF5HssHuL9NfRurjudNjdRKomRkhNWkQyqozqP0IkBS3NF5HssJuK9IFIalqa72iR\nWokRUpMWkYwuo/rP/tpbF9edDkMksXdx3ekwRJJaF9edDkMksXdx3emwu4r0gUhaWprvWSLV\nIiOkJi0iGV1GZWWZkVGTluY7XKRGZITUpEUko8yo/kNEOqSl+SKSHXZbkT4QSUlL80UkO+z2\nIv0lowNami8i2WGIJLUurjsdhkhi7+K602H3FekDkXS0NN/TRKpkRkhNWkQyyox29RsaGTVp\nab6IZIchktS6uO502I1EkjMjpCatK4LX15+v1/6LJqNdIZKKluaLSHbYRCJ9m/Ml1PYLIg2k\nIVIEbB6Rvh57vmXaftFltKufYzsyatJSBBQpVoEAABWeSURBVOeJtA+NkJq0rggQ6RJaCuAE\nkcTQCKlJ64qgKdKf/1ex/scTsQ4XUJtCJDvsCSJhkrIQyQ5bRiRFRrsSD8g1rYvrTofd8dBu\nlwkiNWldEYSLpFghty6uOx2GSGLv4rrTYbOJFHn6u+8hiYxOEUlKDZGatK4Iol+QRSQVLc0X\nkeywmURyZ7QvRDqmpfkikh32DJE0JpHRqSKVkSBSkyZN++3tv/96gys2hkghtCwVozt9ISGS\ngSYM+z+H3t66TSo2ppEEkQ5pWSxme3pCQiQDTRj21wPS13+RGe1L/5BERojkgCGS1Lq4TrVo\nyYwQyQ67u0j6YzsyQiQH7LKTDZbTDcXGECmElqVitqcnpPpBAiI1aeK43wwPSIg0hJaFYnSn\nMyRE6qdFJNOT0b7UT5LICJEcMESSWhfXqRYtmREi2WG3P7RTH9uR0ckibSNBpCZNmvbb6SId\nmkRGiOSAXXj6+/OU09+IdETLYjHb0xcSInXThGGfKZJ2GRkhkgP2HJGO1pERIjlgl4p0xguy\niHREy2Kx69MVUi0RRGrSpGm/fVd4RrVCpCYtS8UqT29IiNRLE8dt8cgp0sFCMipFGvehMUTq\npXWGEJBRpRCpSUvz3Toz8ENjiNRLC7CnM6NaU4jUomXqbKY98C36lV9tiNSkyfM2HNkh0hBa\nlkmRECLpYU842fBZye24dXGdatGSGZ0mUuW7iNSkyR5lXxDpUpok0sgPjSFSJ00Y9k82J4mk\nWklGu9PfpjdEFvU9etV3qd66SKTmUkQ673WkSiA8IjVpwrARqQd2Q5H230akJk0Y9snPkVRL\nyaj2HKn/lFCxNUQKoUnTPvesneo1WTLavyCLSGrYI15HQqQWLcukSKg3n46QEKmPZorCmVG1\ndUQSaWm+iGSHXXuyYUhG1daPH5LIaPccaWBIuzwQqUkThn32WTtEatBkc0wmFVuT5o5IXbTa\npN/yGpFRtXVEkmiSOGNDQqQu2pFIXRFFiCSvJiNEcsCe8hwJkURaFospFVtIZR6I1KQJw75A\npMPVZIRIDthTTn8jkkhL8z3xnQ27HyBSk1YZ9Pdb888+2YBIEk0SaXBIiNRDm08kcTkZ7V+Q\n/ZfPoDNCiNRD68sgKqNq64gk0NJ8C5F+vjP0EennJ4jUpLndMWVUbx2R6rQ03/M+ar7/CSI1\nacKwrzi0O3pIIiPpclyIdD0NkQJg1521s7xJv9gaIoXQ2jM/6UqrP7dDpCrNnocvpO0vNkRq\n0g6GfuJbhD6OHpLICJEcsItFOvPQDpHqtCyOMp2xx9+IpKcdDP1ckdo3IKNzX5BFpB6aKNAF\nJxsOHpLI6MxPyO5+hEhNmjDsCz5Gkd0AkYq//8ZSpNSXTndIiKSnmaIIyEhoHZEqtDTfMz8h\n+1H8XkOkJs1hjSsjoXVEqtDSfK852YBICpow7M2nzfU5FRszi1S9CRmdfbLhA5HUNEQKgN31\nZAMi6WnStH8uWTwsI6H1hklkdPbJBkTS04Rhn385ru1NEGnzd8mcwScbNmkgUpMmDPsykRq3\nIaPTnyN9IJKWJgz7cpEqNyIjRHLALjvZkH1BpEtpkkjGKraGSCE0adrXvLOheSMyukykv+be\nxXWnwy57QbZfo1CR9rcio/MP7T4QSUnrVSUwI6F1REKkEbDHiSTeiozqh3ZDP8aMSEqaNO3v\nywGMzEhqXXpIIiPBmZFPZFMYiNSkCcP+vUDNwIyk1hGpV6Txh3aIdEQThh1yyTStR2XriFT+\n/TcWMSxzSEcpIZKOJgwbkXpgNz7ZgEhK2sQilTckI0mkrowQaQhNmvZXPr5fdogUQ5NEMtYW\nLj0j3f0ckZo0cdxvhgekIJGEcMlooEhyN4ikokUkY85Iah2Rir//1ICPmiNSDC3AHntGYuuI\ntP37rzmbaYecbECkGJo8b/977dwibW9KRiM+ao5IMTRp2hEnhBAphpalUmTUl44tJERS0YRh\nh3weyS5S9aZkNOKaDVqR/iJSkyYMO+QTsn6RNrcloxHXbDgOCZE0NGHYiNQDu0Ak4/XSEGkU\nbVKRaiaRESI5YM98joRIm7/vRXLUFq4WSZFjpXdx3emwpc7abevnV9n5N75rIZId9sjXkao3\nJqNrRNLnWOldXHc67JHvbPio/R4kI0RywBBJar2Htl+0ZEaXiqQyiZCKOueEUKv1XXxkhEgO\n2GVn7RBJDbu5SD3HdoS0rXPextVuvTSJjK55HQmRNDRp2ojUAUMksXdx3emwq042nPI2rmbr\niPT7dyGUiH8NrkMkjUmEVNT1z5F2+ZGR8DGK8SEh0iFNGPYEJxsQKf39N5YipO3XcSEh0iFN\nGPYMJxtKABkhkgP23JMN5UMSGe2v2ZB9GRmS/kkSIZV1/ckGRPr9uxDKedcDQKQjmjDsKZ4j\nFSaR0f6321nvLEakI5ow7AeI9PfvX0WDc2V00TsbEOmYFpGML6Nm6/m+HjiIv5Uyw2Jbq9DS\nfBHJDkOkaJFqFomdzpXR7uInZ11XGpGOaG533Bm1W8/28ZhByBpVm50ro8rJhlNFUvybcYgU\nUVv2hCIdaLTvd66M9q8jvX2ec4UaRDqgBdjjzUgl0t+IQWg0KlqeKyNEssMeLlL2kOQehNqj\nrOu5MkIkOwyRgkQSZDlQaa6MqiKd8hrFX+U6RMrzOfGjLset/0boG0RhSQETXYrJSLa1eqaw\n52TDOe9sUK9DpJuLVO68e5h+R69ssrbNhj3t6jr93f/WhmqXh3cRkdo0adpnfWZM0bq0f3XQ\n9mrUYVaZNjSzQMmkkv6bSqczgSEhUpsmDPvya3/vMQ6RKlJIMHHnbm/ys3nb7irpv7FUQ+qs\n6h0+nKpqDohU1kwiJZNU93BPq+2ijamqd3DNLWxV3KEUy3UhHYxAMdftOtUiRPpX2oy0Iume\nr+xo1R20PdVwM/pKc2h3+hNZ3UJEKkW6+F+jqIFsImU7aL5Iu80rqlekrowQaQhNmvapGal3\nat1drN92e2vNVEPlqN9Ru0ifZ18PQJcBIu3qrM+MaVq3i1TXSD3VCHuk1uQtItIA2NPf2bBB\n6a5RWLnd7pb6qUYIVGutsc0ZT39/fKpWIlJIbdmXi9R6bOiaqtugXWvNVVqRzn2xD5FaNHHc\np33URdW6dk/d0Fo7uGGqLYFOyuiqqwj96wKRWjRp2ud9ZkzXeodJP7TmA8WaGV34OhIitWnC\nsM98h76u9W6R2sdba2Y0g0ia93goaKpFS4a0y2hOkRRvfNuur99gzYwuFUm1FJG2NZ9IepM+\nN6ulqwMtmdGlz5EQqUUThn3mZ8a0rfeIdOTRohld+HkkRGrTpGlP9c6GDU5zjvnQo0UzuuxK\nqx+ZSM21iFSL6KTPjGlbV5qUXUNVXrlmRpe+IKtai0ghtWVHi6Q0SePRohkhkh127ccoTstI\n17rGJJVHi2Z02ZVWPxDpgCYMe6rPI6XWjx3JNDov8CtPNlwg0mm/oNYMaZvPPBc/2bR+ZInW\no0Uzuu66doh0QDsSqSui8SK1Tco1OvOtLE8RSbEYkUqX+gzyZ6RtXValQ6NVM+oX6fX15+uF\nSIi0bV3QZaNR5wW0olobRDsQqXXc8CXQ60eoIJHOOhu6ZkgRtWWPEan8OMP+O6pNLplR7zsb\nvh6JfmTyh4RIIk2KYKILRO5a33nT79GiGfVfabUQ6c//a7vge2TKgPtWU5NdRajcv9oa3fiX\nXf8BNyJdW5O+jlSQ649GjxBJ+wJS7KHd8bEdh3ZLiSSrpKWtmREi2WGIJLQuaHTrjHKRdK/2\nBYt0uByRSpFmfo60wRca3Tqj/ksWx57+RiSRJgUw6zsbiiolUtPWzKj/ApGxL8gikkiTI5jp\nSqvV1l20NTO68kqrG5HO+cDkmiFF1JaNSDG0NN+LP4/0gUgiLSKZoIx6W3fR1sxIeGfDmU9k\nEalOE8c92QUiK627aGtmtLuKECKpYReebECk2TK6/jnS0bEdIW1rwuvakREieWCIJLXuoq2Z\nUfkcqSucoJAQqUoTho1IPbDLztqZTNrCESmGJgx7xiutktEMJxsOju0IaWfSEu9sMNPWzAiR\n7LALT3/PdqVVMrr+BdmjWxBSSG3ZiBRDS/NFJDvsfJEsLyAFZdTVuou2ZkYziXTCvzu1Zkgp\nLUTqhF0gkvbzSCNCQqQKrSrSvFdaVa07HYZIYu/HtNaiJUMSTEKky2k1kRy1hSNSDE0YNod2\nPTBEEns/prUWLRnSphCpB/ZMkaq3IaSQ2rIRKYaW5juHSK3bEFJRPEfqgCGS2LuC1li0ZEib\n4mRDD+yhItVuREjb4jlSDwyRxN4VtMaiJUPa1BVXeupt3UVbM6NJRGrciJCKOo4s+J/e6W/d\nRVszI0Syw646a3eUWfQ/vdPfuou2ZkaTiVS5FSFtS3s1XEQ6gZbFEuARIg2hCcM2iOT9p3eo\n45pFJPlWiNRfiHR6zSbS/maI1F8c2p1ESyNHJDtsXZEarzaEtO6irZnRNCKJNyOkLC3lxyiO\nTn+rI0KkNg2RImCTi9R4QRaRgmgTi7S7HSGF1AaNSEG0NF9EssMQSWrdRVszo3lEkm5HSCG1\nQSNSEC3NF5HsMESSWnfR1sxoPpHKGxJSSG3QiBRES/NFJDsMkaTWXbQ1M5pIJOGGhBRSGzQi\nBdHSfBHJDkMkqXUXbc2MJhSpuCUhhdQGjUhBtDRfRLLDEElq3UVbM6OZRKrfkpBCaoNGpCBa\nmu+MIm1vSkghtUEjUhAtzReR7DBEklp30dbMaCqRqjclpJDaoBEpiJbmO6VIm9sSUkht0IgU\nREvzRSQ7DJGk1l20NTOaS6TabQkppDZoRAqipfnOKVJ+Y0IKqQ0akYJoab6IZIchktS6i7Zm\nRpOJVLkxIYXUBo1IQbQ0X0SywxBJat1FWzOjSUXKbk1IIbVBI1IQLc0XkewwRJJad9HWzGg2\nkfa3JqSQ2qARKYiW5jurSOnmhBRSGzQiBdHSfONFcoaESOnvkbVBI1IQLc13OpF2tyekkNqg\nESmIluY7rUi/AEIKqQ0akYJoab6IZIchktS6i7ZmRvOJVAIIKaQ2aEQKoqX5zivSD4GQQmqD\nRqQgWprv00TanRR0wBBJat1FQ6Sv8oc0SKS/u3LAYlur0iKSGZZRu3UXDZG+alKR9hq12pwr\npAkzarfuoiHSV4WJ9Lfeey9tg9S6NFdIE2bUbt1FQ6SvmlCklkb1XucKacKM2q27aIj0VQEh\nbXZw/ySOPKq0O1dIM2bUbN1FQ6Svmk4kjUe76yTrNolI9dZdNET6qoiQAkWqGnPo0lwhTZlR\nq3UXDZG+KlCkv7Xe+2iFKjmspdJcIU2ZUat1Fw2RvmoqkUpNtjBZpblCmjKjVusuGiJ9VUhI\n4r7fR9spsoMJKsWFJDzkKWgRyYzMqNG6i4ZIXzWRSJVHmv2q6sOSK6QKUawGLSKZkRkdD8JG\nQ6SvigkpQqTK/lqH7fdv21h7DKoUIsXDEOkXY55E7de+BKvs1ppNJprTIUQaBEMkt0hVJWRY\nbc8+3OSncENEctAQ6auCQvKKVNehBavv3cctIhIizS/SX+MkBBPaMHEXVy90FSLFwxDpdx82\nTUJy4Ag2wo/f1o4WIlI8DJF+dzvLJPK9eLvoGBYkj+aOItJwGCJ5RBI90o01zqB9a607gEjx\nMETKTFLdxapI+0VKWJRCu9ZaqxApHoZIDpEa+3bHWGMcKlprrkKkeBgifWTnG3T3sbxZ7Xad\nY3UrVLTWXoVI8TBE+rCK1NzJfWcuSoN497fQuouGSF8VGFKXSbsHkeqiJUOaOaNq6y7amhnd\nTaSDo641Q5o5o2rrLtqaGU0tUpdJhUjCoiVDmjqjWusu2poZ3UykA48WDWnqjGqtu2hrZrSG\nSFra4fo1Q5o6o1rrLtqaGc0tUs9D0qfGuzVDmjujSusu2poZTS5Sh0mfmsevNUOaPCNE+lxH\npGPe77/YgkgnZ4RIn/OLpH9I0ni0aEizZ4RIC4l0BFR5tGhIs2eESAuIpHxI0nm0aEjTZ4RI\nC4nUJCo9WjSk6TNCpAVE+tQ4oj0AXDOk+TNCpAVEOjZJ+3i0akgLZKRahkjXhnTkid6jRUNa\nICPVMkS6uH5Faf+0/mPqX23yR6QgWppv/CPSiN5bjzkdj0erhhSeESIF0dJ81xBJtuVvl0eL\nhhSeESIF0dJ8FxFJMCnTSHeNlDVDCs8IkYJoab6riLRx5vsnfzffu3FI4RkhUhAtzXcZkQ6u\n6HPrkMIzQqQgWprvOiI1TdLS1gwpPCNECqKl+S4kUsMkNW3NkMIzQqQgWprvSiJJJnXQ1gyJ\nmr+WEqlxrUZEmiUjN23i1u7yiPRRUamPtmZIi2XkpE3c2o1EKk96d9LWDGm5jFy0iVu7lUgu\n2sStIdIIGCKNoU3cGiKNgCHSGNrErSHSCBgijaFN3BoijYAh0hjaxK0h0ggYIo2hTdwaIo2A\nIdIY2sStIdIIGCKNoU3cGiKNgCHSGNrErSHSCBgijaFN3BoijYAh0hjaxK0h0ggYIo2hTdwa\nIo2AIdIY2sStIdII2FoiFfXnz7S0iVsLpp26MTJCpLvSTt0YGSHSXWmnboyMEOmutFM3RkZr\nXESfomYvRKKogEIkigooRKKogEIkigqoASK9Xq84zOv1+9UJyzE+2uubFtLb6we5/zKwyKgL\n9oPcf/mteJFePxv2YxLKC91iIloM6u07lU1fQRNsbZWMejCajMJFev3+4ef8duuGbjARLb5i\nevv+dbnpK2qCja0GbYGMUs0q0k/z3w/QXugGExJSCXWBNn0tIxIZZTWrSP8or5+jh4jfdq/I\nPfX3wd5NW1ikfxQy+qqJRfp9Nhe0b4WHFEJbWyQy+ql5RXpl//fwkKKLjPpJy4r0ymleaPSe\n+sr+eLBIZJRq1tPf6ZDh9zDc21TCuFv83YP8vf1E8tp9GVhk1I86yGjWF2TzF9QCoK/2q2nd\nuBrUiQpsT7FVMrKgGu3xFiGKCihEoqiAQiSKCihEoqiAQiSKCihEoqiAQiSKCihEoqiAQiSK\nCqhbi/T2W//9/9XdULW6S0YLt35cdwnpznWXjBZuXVcrh/OUukNGN7gL7bpDSHevO2R0g7vQ\nrp+Qfg4d/v39+5tfBxTU1XWHjJZo0lNFSP9y+Q7n59CcurjukNEKPbpqG1L64+cnS6R087pD\nRgu06KvysCH9sU5Id687ZLRAi75qhbTMccPN6w4ZLdCir45+21HX1x0yWqZRa90hpLvXHTJa\nplFrNQ8bPleK6r51h4wWaNFXckifyxx/373ukNECLfqqEdIyL/bdve6Q0RJNUtTshUgUFVCI\nRFEBhUgUFVCIRFEBhUgUFVCIRFEBhUgUFVCIRFEBhUgUFVCIRFEBhUgUFVCIRFEBhUgUFVCI\nRFEB9T+H+cmUVO7uJgAAAABJRU5ErkJggg==",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g1 = ggplot(as.data.frame(sir_out[,c(1,2)]),aes(x=time,y=Xh))+\n",
    "  geom_line(lwd=2,colour=\"green\")+\n",
    "  xlab(\"Time\")+ylab(\"Susceptible humans\")+\n",
    "  scale_y_log10()\n",
    "g2 = ggplot(as.data.frame(sir_out[,c(1,3)]),aes(x=time,y=Xm))+\n",
    "  geom_line(lwd=2,colour=\"green\")+\n",
    "  xlab(\"Time\")+ylab(\"Susceptible mosquitoes\")+\n",
    "  scale_y_log10()\n",
    "g3 = ggplot(as.data.frame(sir_out[,c(1,4)]),aes(x=time,y=Yh))+\n",
    "  geom_line(lwd=2,colour=\"red\")+\n",
    "  xlab(\"Time\")+ylab(\"Infected humans\")+\n",
    "  scale_y_log10()\n",
    "g4 = ggplot(as.data.frame(sir_out[,c(1,5)]),aes(x=time,y=Ym))+\n",
    "  geom_line(lwd=2,colour=\"red\")+\n",
    "  xlab(\"Time\")+ylab(\"Infected mosquitoes\")+\n",
    "  scale_y_log10()\n",
    "grid.arrange(g1, g2, g3, g4, ncol=2) #Reproduces Keeling & Rohani's plots as seen in their original MATLAB code."
   ]
  }
 ],
 "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
}