{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Forced SIR model using R and simecol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(deSolve)\n",
    "library(simecol)\n",
    "library(reshape2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sirforcedode <- new(\"odeModel\",\n",
    "  main = function(time, init, parms, ...){\n",
    "    with(as.list(c(init,parms)),{\n",
    "      # ODEs\n",
    "      N <- S+I+R\n",
    "      dS <- mu*N-beta(beta0,beta1,omega,time)*S*I/N-mu*S\n",
    "      dI <- beta(beta0,beta1,omega,time)*S*I/N-gamma*I-mu*I\n",
    "      dR <- gamma*I-mu*R\n",
    "      list(c(dS,dI,dR))\n",
    "    })},\n",
    "  equations = list(\n",
    "    beta = function(beta0,beta1,omega,time){beta0*(1+beta1*sin(omega*time))}\n",
    "  ),\n",
    "  parms = c(beta0=10./7,beta1=0.05,omega=2*pi/365,gamma=1./7,mu=1./(70*365)),\n",
    "  times = c(from=0,to=100*365,by=1),\n",
    "  init = c(S=99999,I=1,R=0),\n",
    "  solver = \"lsoda\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate until equilibrium\n",
    "sirforcedode <- sim(sirforcedode)\n",
    "# Reset initial values\n",
    "init(sirforcedode) <- unlist(out(sirforcedode)[100*365,2:4])\n",
    "# Look at 10 years\n",
    "times(sirforcedode) <- c(from=0,to=10*365,by=1)\n",
    "# Simulate\n",
    "sirforcedode <- sim(sirforcedode)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sirforced_out <- out(sirforcedode)\n",
    "sirforced_out_long <- melt(sirforced_out,\"time\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {},
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAElBMVEUAAAAzMzNNTU3r6+v4\ndm3///+pfg9ZAAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3di3brNrJF0bRJ//8v\nd2L5IZIA9gYKlEBw1Rg3JzGBehCcliwnt//5JAgiHP+8uwGCmCGARBAdAkgE0SGARBAdAkgE\n0SGARBAdAkgE0SGARBAdohHSosNaFInT859d4OwJGEBXiO4HkpF/+GNU6RlAVojuB5KRf/hj\nVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4Y1Tp\nGUBWiO4HkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlA\nVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI\n7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4H\nkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlAVojuB5KR\nf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4\nY1TpGUBWiO4HkpF/+GNU6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU\n6RlAVojunwnS+m+ckn/4Y1TpGUBWiO6fCNK6niRp/GNU6RlAVojunwfSup4lafxjVOkZQFaI\n7p8G0gqkfHoGkBWi+4Fk5B/+GFV6BpAVovuBZOQf/hhVegaQFaL754PUX9L4x6jSM4CsEN0P\nJCP/8Meo0jOArBDdPwukFUiF9AwgK0T3A8nIP/wxqvQMICtE9wPJyD/8Mar0DCArRPcDycg/\n/DGq9AwgK0T3A8nIP/wxqvQMICtE908Iqbuk8Y9RpWcAWSG6H0hG/uGPUaVnAFkhuh9IRv7h\nj1GlZwBZIbofSEb+4Y9RpWcAWSG6H0hG/uGPUaVnAFkhuh9IRv7hj1GlZwBZIbofSEb+4Y9R\npWcAWSG6fxJIK5BK6RlAVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlA\nVojuB5KRf/hjVOkZQFaI7geSkX/4Y1TpGUBWiO4HkpF/+GNU6RlAVojunxFSb0njH6NKzwCy\nQnQ/kIz8wx+jSs8AskJ0P5CM/MMfo0rPALJCdD+QjPzDH6NKzwCyQnQ/kIz8wx+jSs8AskJ0\nP5CM/MMfo0rPALJCdP8ckFYgFdMzgKwQ3Q8kI//wx6jSM4CsEN0/FaQFSOn0DCArRPcDycg/\n/DGq9AwgK0T3A8nIP/wxqvQMICtE9wPJyD/8Mar0DCArRPcDycg//DGq9AwgK0T3A8nIP/wx\nqvQMICtE9wPJyD/8Mar0DCArRPcDycg//DGq9H3yF27t6wfofNBAWhKQOksC0leUbu3LB+h9\n0EBaniGd85IEpP+i+E3q1QN0/5dYgLQASacHkqwQ3Q8kIz+QFiCJ/UAy8gPp+T9USVaIFygF\nkB41gk3KAkAS6YEkK0T3A8nID6TBIK3lbloqRPdPAWkFUjk9kGSF6H4gGfmBBCS1fyZIC5Ay\n6YEkK0T3A8nIDyT17AIJSEZ+IKn/aVEg1UH6+PrLxweQOlc4OT2QZIXo/ipIX4I+vj0BqWOF\nk9NPDKnTWb8U0scvoj9JVo1gk7IAkET6ySCtk0H637+hd70gvm/p5u+IrvH86L67l00z7+9m\nFxrSxyeQbhtDPbrXhvRxgPRfWK96wZdNWYC3diI9b+1khej+Gkj/BZDOCCCVYypI34AuAamv\npKYBqrq4FqTEXECqhjT2x9+nvCS1DFDXBpDKMSWkoX8hOwqkylMeHtIqnl0g1UE6hlUj2KQs\nACSRHkiyQnQ/kIz87ZDMRoBUDiA9agSblAWAJNIDSVaI7geSkb96gNpjBlI5gPSoEWxSFgCS\nSA8kWSG6H0hGfiANDKnPWQNp2X5EBqRU+rkgraqblgrR/RNBOvx9r4hA8joBUjmSkLqeNZAW\nIOn0QJIVovuBpKN6gOo3HkAqB5AeNYJNqgCSSt8Z0nEoIAFJB5CAJPcDSQeQVvXsAglIOoAE\nJLkfSDqA9DsLkHL7gaQDSECS+4GkIwTJagVI5QDSo0awSRVAUumBJCtE9wNJB5CAJPfPBukM\nSUACktwPJB21A2wcAalDFCF1OWsgLUDS6aeC9HtjgdQ5gKTSzwmp51kDaQGSTj8jpAVIvQNI\nKj2QZIXofiDpABKQ5H4g6QDSHtJhKCABSQeQgCT3A0lHDJLTymUg5W4vkICkoxFSRStAKgeQ\nHjWCTYrY3k4gJdIDSVaI7geSDiABSe4Hkg4gAUnuB5KOygHWjpD6jBI9gfVvEiBl9gNJx/sg\ndZoFSLpCdD+QdLRBqmklU2DtNEwnSLu/3VSIFVABpEeNYJMigCTTA0lWiO4Hko53QfpxBCQg\nPWoEmxQxP6ToNJND6nHYQAKSkX5SSB0PG0hAMtLPBGnV3bRUiO4Hkg4gDQjp8PfBCtH9QNIB\nJCDJ/dNBOkFSEJLRSrLAnyMgAWkBEpDCAaRHjWCTIoaF5LeiIAWnAZKuEN0PJB1AOkDaNwQk\nIOkAknx2gQQkHUACktwPJB1AApLcDyQddQM8PftA6hRAetQINiliTEhVrQCpHEB61Ag2KWJS\nSM+OgFTTTFOF6H4g6QASkOR+IOkYAFJsnOAJbHoAUnr/9SHtbyaQjum7QEr8w1+FUAEZQHrU\nCDZZDiDp9ECSFaL7gaQDSECS+4GkA0hAkvuBpANIg0PqcNhAApKRflZI/Q4bSEAy0gNJVoju\nB5IOIA0EaVsfSP1iUkhbR0BKNgOkjgEknR5IskJ0P5B0vBNSl3GApCtE9wNJRxSS7kVDCs0D\nJF0huh9IOpoh2b0AqRxAetQINlmOySH1mAdIukJ0/7SQOkqaAFIoAZCM/fNB6v+SVDXAxvEo\nkGIZgGTsB5KOFkh1rZwMKZgCSMZ+IOm4OqQVSLpCdD+QdMwCqTWH8ewCCUg6Lg7p1xGQjvGz\nG0hAMjsCUiJ+twMJSGZHQDrG340BEpDMjppz7DZPBOnpxgBpUkh/R9wRUluSG0BagTQ5pIoc\nqiMgHQJIzwEkqyEgHQNIzwEkqyEgHQNIzwEkqyEgHeL5xgAJSF5Dc0IK3R0gbQJIVkNAOgSQ\nNnETSOFfOgJpH0DaBJCshiaBtK/eCdIKpCtAkr2UIPX612CakwDJ6vDqkI6PyECQ3F6AVA4g\nPWoEmyzG8U4CKZEISMcA0nMAyeoHSIfY3BggAcnqB0iHANImgGT1A6RDAGkTQLL6AdIhgLQJ\nIFn9AOkQQNrEnJCeH3wgvQBS9HEBkhHvg1STo9wPkA4BpE0MBsl56hIVSmkGg5TsB0hA0jER\npJYsvSE1tAGkR41gk8UAktfQMJBa+gDSo0awyWIAyWtoFEhND+4ZkLaOgAQkr6FBILU9uCdC\n6vS4AMmIe0M6bOwDqaoTID1qBJssRhZSP0lAApLucD5I3V+SxoDUNs/vZiClOwHSd8wPKTTP\n714gpTsB0ncAyesHSOlOgPQd3SEdtwKpG6T1KSp6eAGkoCQgpfNtvgQkIOkOgZRKt9kMpJkh\n9XlvB6RkNiBlWsm1AyQgpbJtdwcgmb28DlJ9lqEgWc1UtQKk7wCS1w+Q0q0A6TuAZPYzJ6QO\nb1mB9BVAMvsBUrIVIH1HV0jpY349pG0Tk0KqaOUGkN4e37dBfKkuWePuVO22bNtdkY7+9rZm\nOewLtLNxFDig4lfqWond3kTwirTk3nm87RWpKofqpzXLYV8qkXmHtpD8Hk54RfprInJ7nzq8\nOKTUoQAplQlIqVYCGXYdzgFJfq0q2zyQnoYBUqqVQIZdh0A67gRStpdMIiABKbFzOkifgSxA\n8joE0nFnOyTrqTsGkMrxCkgxSUBKbARSvpcl/eABCUjHjUDK95Lp58qQury3A1Ji46CQAu9V\ngZRrJZJi2yGQDhv3+98NKf5D35iQ/F6A9KgRarEcQHL6GQTSzhGQgGS3Mxyk2iyJbVFI9QMB\n6VEj1GI5gOT089meJbELSKkOgXTYCKR8K5kvAQlIx41AyreS+RKQgHTcCKR8K5kvAWlGSK13\nZusISLnaM0B6PuR5IAXGOAHSYTuQpoUUSbHrcABIkTmA5PQDpGwr7Sl2Hb4f0hoZZHRIXi9A\nKgeQHjWsmdomORHS7/7rQvqZBEjZVtpT7DoE0nHbbJA+gZRtpT3FrkMgHbcBKdtK5kszQApJ\nGgDSryMg5TZNAunvnPtDiow1HaTAT9PGFytyASnbSuZL74PUYywgTQppv2cuSMszKS+A9Kjh\njASk/B4gAWkBEpCO+xohpR4UIG1HAlJ+D5C2kFSDdb1MAunJEZAye4AEpAVIg0GqTAMks8OR\nILU/KsYXK3I1Q/KeumO8BFJbGiCZHV4bUnpj6N4OCKn5pIFk9AKkJXcLOkH6SQAkIOkOgXTc\nBaRy6ePXgASk1C4gFUtfH9L2gQMSkKx+gFRopTXFvkMgJXYBqVgaSIkO3wxp4whImT1Aqmmm\nrpXWFPsOgZTYBaRdhh6QnlMBCUiV7bwT0s9tjUNSXwTSWJBaHxXrqxW55oK0DAWppRUgPWoY\nEwGptGdWSG4vQHrUMCYCUmnPdJAqewHSo4YxEZBKe4CU3wWkzQzTQrKaAVI5gPSoYUx0BqT2\nmwukYuVJIUUkzQkpfHOBVKw8G6QOL0lA2u8BUrmV9BeB9GZIv0/vMJAS+18N6bgDSPldQNpM\nMCSk7wRAApLuEEipPUAqVR4RUmis+SAFPrI2vuynAlKp8lCQeowFJCA57QCp2MtUkBpnAZLT\nDpCKvQAJSF47xefQSyG+CCQgJfcMBil4d2aBtPaCtN8CJCA57UwGyWmxopW2FIcOgZTaA6RS\n4SZIK5CAFOkGSKlNQAJSZTdASm0CEpAquwFSatNa1QuQHjXUQEBSO5ruzu/dPQnS01erIdX1\nAqRHDTUQkNSOCCTxaxgzR/mrQBoCUvMsuU1A2mwCktgCJCA57QBJbAFSP0grkIo5yl8FEpCS\nW4BULAykY4dASm1Zn/8RSDeBFJA0HKTAzzXG1/1UQCptmQ1S/CUJSJktQCptAdKxw1EgxT9p\nM77upwJSaQuQjh0CKbmlHySnGSCVA0iPGtZAQCpsAFJFM9W9AGlKSKn174K0AsntEEjJLUAq\n7QDSsUMgJbcAqbQDSMcOZ4YUeGRmhlSRBkhmvBfS5mB7Qop/750SUn2aOSGthy1AApLRDpDU\njjkgtc8CJKudW0AKTQUkIBntDAzJ66UIqcdUQAKS0c6IkKp6AdKjhjUQkArrrw/psAVIQKpr\nB0jJLUACUl07QEpuqYe0AskZKHebStEP0r76BJD+9swFyW2yohcg9YaU/GcXUmaCHpBaThpI\ndi9TQYo9K+4FN1MAUkMvQCoHkB41vIGAVFgPJCABCUhAAhKQsiP9fX0SSO2SgJTZAKTihtkg\nhV+SgJTb8PQFIMmvAwlImQ1AKm0A0qHDK0NagWS1AyS5AUj/RiI/kIBU0QuQgOS0AyS5AUhA\n0u0ASW4AEpB0O0CSGy4NaQXSdzwXWJPrgQQkBal9lpkh1SfJ1/0boDYNkMwAUm4DkErrLw0p\ntR5IQNLtNENaT4RU0wuQHjW8gYBUWP5WSPoCkMaDVP/0XwNSqZn5ISUyAekkSM1P/+CQjGaA\nVA4gPWqYA70R0gqkcgpxAUhzQwo8MUAqLQfSocM5IcWfGCCVlo8IKTQVkICk25kcUo+pgAQk\n3Q6Q5HIgAUm3AyS5HEhA0u0ASS4H0k0gNZw0kPxegHQapORzWJ3E7AVI5fgaYAWSORCQCqtb\nb85X+pkg1bRp97J8BiVVQfr4+Hj6A0huTSABaePoP0OPP4BkJjF7AVI5poP07/99/PwtkNya\nQALSEdMTpP/9G9audKw/DnZfqM5QeclcXt1NdkddpszqPjcnflv6dVObI7++x1T1x50JB9LH\nz9u6npAKX6jO0JwLSDW9AKkQ1ivSR8+3dt+tl19yi/EzvJPcypT+kvnWbs3V1L3w1q4ck7+1\nA5JXEkhA2r4aAakyh9sLkMoxFaTuH3/nIPnDAMlrB0hy+ZV/IXtsHUj51X0gVaY5FVJFL5NB\nSoVVw52n7YbMBSm3eEhIv1eANAGk5CkC6VRI+yvqDiUTAQlIZg63FyCVA0iPGu48QAJSMoD0\nqOHOAyQgJQNIjxruPEACUjK6QkquBhKQdDutkFYguQGk59UnQdLNjA3JuAQkIOVWAwlIfgAp\nuxpIQPIDSNnVQAKSH0DKr/79IpDkJSBNDinyxADpKpA6TAWkDKQO9xZIQLIjB+kf15VVozQP\nkJYaSP5JA6mml8/671S7Dt8NSY6Yj1tBar05X+mBpBafBcl+q2fVcOcBEpCScV1I//wGkICk\nLwEJSNnFQBoK0npJSHZYNdx5SveqkAFIouZMkGobNRfPBalyGCCZ3QBJLj7xw4b/3tXpD++s\nGvY8QAJSKi4M6R8gAQlIduQgfSn6NCRZNex5gASkVADpUcOeB0hnQ6pKAyQ7gJRfXAlpBVJl\nohtBevHvkWaA1NQLkMpxXUif5u9jgeQUrICUXzotJKOXC0P69BwBySkIpHSiW0Ayw6phzwMk\nIKUCSI8a9jxAugqk3bULQUqvPRGS9c5uDkjptUACkh1ZSLf6sAFIOkf5GpAykL4J3eTjbyDp\nHOVrQMpC2v55Y0grkPQ1IAEpv7YbJNnMmJBWINkBpMJaIBWWAmnXYRrS58/PSMJRO6Rk30AC\nUiquCemfTYwLaQWS2Q6Q5FogrZlTBBKQKtaKAkaHR0hVYdWw5wESkFIBpEcNex4g9YK0WQ0k\nufY0SLy1Gw1S4835Sg8ktfYsSPyMtAAJSH7kIQlAQKopCKS7QrJ/ZLJq2PMACUipuDCk0//X\nKNJtN90RIInFd4EUmGqyDxuABKRUKEjxqSb7sKE/pMgTAyQguZGHJAAND6nDEzMlpIo0QPIj\nB+lNHzYAaUpIuUR3gHT5DxuABKSKpfyMBCRdcVRIdi9AetTw5wESkBJxXUh2WDX8ea4LaQVS\nfaIKSMX7C6QhIOVWtkBquLgAScWn+u/uqofirV1iLZDEYiDJlUACkq4IJLny5Ld2J/6vUQDp\nKYBUjstD0pKsGv5AQ0H6ufBCSOvVIP1UaEt0I0hnvbXLdA2k/MKxIO0qtCUCEpCqygHprpBO\n/7DhlpCyzQCpHBNAEo6A5JUD0k0h2WHV8OcB0jiQViD5AaTSSiAVFgJp2+ER0kv+XxYD6TmA\nVA4gPWr48wCpD6TtWiDJlee+tXvbhw0VdwRIei2Q5MpTIRmOekNquiP3gmSdNJAqezkRksMI\nSGa5HpBqThpIlb18ViVJdZiB5PwOCUh2OSDdFJLJCEhmOSDdEpLNCEhmOSDdEVKFIyB55YB0\nQ0jv/D0SkICUCBdS+1RAWoAklw4LyW1GQgpPddqndm5YNfx5gASkRADpUcOfpz+kyBMDJCCZ\nMTEkO5U85rkgdbgvx8tAAhKQ2nr2IRUSAQlIVdWANAik3DogLUCSS4Ek110VUrbnoSD9SrK7\nAVJdIhuSuL9Asi8UVg4EqfEykMrxC6m1WbkOSAuQ5FIgyXVAWmaCVF4GJCC58wDpapAe14EE\nJCA1lgPSpkMgAamtHJA2HQLpFZDEdSCVA0iPGv4823fexQCS3Q6Q5LrpIPnTAMluB0hyHZAW\nIMmlTZDkMiBtOgTShJB2K4Ek1wFpAZJcCSS5DkjL7SBV3ZtHeiCpdZ81SZIdAukykFruzSM9\nkNS6+0JagWR3AyS57vaQvgrIFW6mzCUgqQVAAhKQGosBadMhkIDUVqwHJLMZID1q+PMACUiJ\nANKjhj/PGyAZxwwkteDtkKJTAWkBkly4OwEvDZBqAkhAai3mQiomAtJlIEUeGCBNBim7CkhL\n/hSB9JMeSGoVkBYgyYVAkquAtABJLgSSXAWkBUhyIZDkqotCWoG0CSCVQ/6//r45pNTqESFV\nddOy4NPLcntIkXbLq4C0nAmppZuWBUAqB5AeNfx5/BcBIPntAEmumg+SOw2Q/HZOhVQ8ApkI\nSECqWjchJOcIZCIgheK74+prdcu8TMVVZjN6nZVILXK7EQs73Jeadnrc4F7N9LgxInhFCjTD\nK1JrIq+XG7wiWTX8efxpgCSLAclfBaQFSHIdkOQqIC3XgFReAaRyvAhSQBKQgNRcDEhP+4EE\npNZiQHraDyQgtRYD0tN+IAGplONsSF4zQHrU8Ofxh+kFqfzA3B6SuQRIV4X0vKoHpEgKax2Q\nmusA6aeGP489zC0hNbxaP9IDSa0C0gIkuWx/Ak6abpBEos6QWk8bSAuQ5LI7QAqeNpAWIMll\nQJJrgLQASS4DklwDpAVIchmQ5BogLYVTjD8wQPKWAOktkIr9AilSJrkMSHINkJZXQKp+dCtX\nAKnYC5C+a/jzLO6z+yJIrY9u5Qr3GIEUavhekOo/IwCSWAYkuQZIywiQoucMpGIvvRoGUnkR\nkMQyIMk1QFqAJJcBSa4B0gIkuewkSE9rgAQkIDWXAtLffiABqbkUkP72AwlIzaVmhdQgCUiv\nglRcA6RinZdBan9JAhKQ2ksB6W8/kIDUXApIf/uBBKTmUkD62w+k+SAdVgFJrgHSAiS5Ckhy\nDZAWIMlVGUjhSkD62w+kG0LqVcmBpBIZqIH0U8Ofxx1msyYCSZ3jWJAav3+PDOkFSJw1QFq6\nQGq+7q8CUnultrFrlgBpAZJcBCS5BEgLkOQiIMkll4QkugVSrBCQGpYAaQGSXASkzYq7QKp+\noAKn2OecgXQpSKlLQFqmgXTeEwUkuQJIC5DkIiDJFUBagCQXAUmuANICJLkISHIFkBYgyUVX\nh+Q93kDaB5A6dwMkuQJIC5DkIiDJFUBagCQXAUmuANICJLkISHLFp15SDCC9GlL+GIEUagZI\nh7gcpPBBAwlITg1/HncYINUsApJcAaQFSHLRWZD+VgEJSNNBOq4BklwBpAVIck3uBKKfsQDp\nb//NIfV5oi4HqccRbJad+d9WAum3hj+PO0xnSG6lwBogtS8B0ncNfx53mLtCanswgSRXAGkB\nklxyK0htp31FSHLg2jsLJLHkHpBCp31hSPkNQAqVAlLDAiAtQJJLgPR8HUi5FW//JQaQrgQp\neQ1ISwlS9881mpeUVwEJSE4Nf57nBTY1IIklQJLXp4RU+5oFJLEESPI6kBYgySVAkteBtMwC\naQVSr24argNpmQxSe5L8kotDWoH0XcOfx1wBpKolc0BSzQDpEEDq2w2Q5HUgLUCSS4Akr+8g\nVUsC0usgFZYBaRRIrS9JQAJSY5nNMiABKa4kftBAApJTw5/HXAGkqiVnQ1qBBKTpIK3HJYlj\nEnlqRwLSvSElHrr6ZsaEtEn/NkjODe718gikQ7waUqGVe0Kynv/nNAKSl6T1eo88QFqAJFe0\nQhK9LED62w+kS0Kq+4kCSPL6BSEZnQIpUAxILdeBtABJLgCSvA6kBUhyAZCeLgMpu+BcSE2P\nbuU6IL0OUvri1JDcWwsksQBI8vKckCpvLZDEgptAigwFpAVIcgGQ5NV7QtpdngLSCqTW6+ZC\nIKnLM0EKJMkvAJK8CqQFSHIBkORVIC1FSD0eGCDpYkACEpAaiuwXAglIL4eUWAgkIDk1/Hm8\nNReFlF8IJCA5Nfx5vDVAAlL9QiCpy0ASC64NaQXSTw1/Hm9NL0jeEQFJF3sBJN1MT0i1koCk\n75la1eGgO0JKPQinQ1pDkHo10wVS40sSkCaFtE1/GqSflUAC0q0hpfMULyaXAglI94RUzFPx\nLAHpez+QgFRd47g0cwQrkIBkr3otpJbPuN4NyUzSdLVTHiAtQJKXgSTTXA/SWgHJvLVAEpeB\nJNNcFpLYA6TWckACkr0KSEAqrGwZCkjLHJCKz0BNuTtDCgwFpOVVkKoe3cqVQAKSU8Ofx1q1\nf+hOhtTy6Fau/OyRJH8ZSDLNnSH9FeiVqCUHkIC03AdS6IEBkuwGSDWQPj4+nv64BqQeD0xv\nSMelQLoTpI//DD3+AFJ9ltJSIN0M0r//9/Hzt0Cqy1JaCqQbQfqYDtIKJCB9x4s/bPjYQPrf\nv+Hs2sTPnQ+s8lLoxauZSCyraCe/1E/SqZvSui4TVSQqHoR7Sqpe5Y0xH69jOJA+fn4+mgpS\nJEVlO0A6tZk+ec6H9LH9y1dYr3r+K6y1an8x+NZOtGK+8xjkrV2ym8Jbu9CNEWn63OGKZrq8\ntWt7b1cD6WOnCUhVWUpru0Papa+DVMZR0w6Qso5m+vj71pAKeaqeJCA99ldA+orPiX4hC6ST\nIa1AssOq4c9jrQKSf28e6d8L6UVJymsF6skhebf2XpAa3kwBSaa5HCTxjcGaBkhAqs4zKSS5\nCUht7QAJSJu4GaQVSB0hNQwFpOVlkGoe3cpmPnskyV+8DaT2oYC0AEleBBKQnGtnQyo/5kAC\n0jUhrS2QIg8MkBq7AdIVID0VaMwEJCD9BJAWIMmLQAKSc6kR0gokIP0EkJYLQjosBhKQnBr+\nPM66zpBkKx0h5RYDqeK427oBkr4EJHGxBClyY9RiP5EE8HpIdZKABKS6AjWLO0LymukCqekl\nCUi3hdThxjwtDr4rA5JVo6HJy0AqPUh2zQtDqjqphiRAeqrR0OTVIOkspdVzQlqBNAwk65YA\nSVysg1Tx/JfbqcgDJKvGoUf/E1HnCpD214BUmwdIywyQqp7d+m6aIFm9FJePBEnd4YtCMnYB\nqTJH4drtINW/zAJpmQhSKEfh2n0gNQ8FpOV1kLzPEFuaaYFU0Q2QgORcOR1S6SXp9ZBKq4EE\npF30hhR4YIDU1g2QhoZ0hAEkcQlIQHIuiCYGhLRfDiQgOTVaegSSSuJfAhKQnAsXgpRZDiQg\nOTVaejwd0gqk4o0xmwHSYz+QdCtAaoXn6/YAAA1+SURBVOumByT/lMqJGiBVSQKS0cqkkPp8\nhNgVUhbASyG1vCQByWgFSG3ddITkPGfFmreH5NwSIIlL74GUPcGaLEB6rtHS4kUgVT0t2aoj\nQeozUW0iIDk1mlq8FCQjS2F95bM7JKT6Z9daDKTnGvsOgZT8YhxSJlELJLOXQSFVf3e4JiSr\ntftBCuUoXTkVUvOLQE2STpDy+4C0AEleuRGk1qGAtLwQUv4naiNLYX0bJL8bIAEJSNkcQPIT\nAakfpPYH5t8BhoJU3Q2QgJR7IamAFH5gToC0pr5Ym8S+ACQgzQYpuQFIQLJqNHVYhrQp0JQJ\nSGtyOZCAVJMp/RQlA0hN3QAJSNuYFVKfTz7mhVQj6faQnE6A1NRNC6QMACA9avjzHKcxvgwk\ncQFIdZAa3tsByQggtXSToVGX5ccRkBYg6Q0jQap8/vNFgQSkXQDJ6ab62XWWd4Kkh7oWpKpD\n8s/n6pCqn90XQLJ7GRRS7XcHIC1AkonOhdT6bqoiSSdIhW2XhGT2dgVI1QYKkEI5il+/E6TG\noYC0vBJS9CfqdF0gqSRAetRo6+9CkJws2x3r8UuhHMVEQALSIJAqPkNUAaT6JEB61GjrD0gV\nOYqJgASkjpCaPykDUks3QJoOUuEHfCAB6StuBykFIADJamR2SLGPEIH0tR9IOqaFFL0xxS1t\nkJIA3gWp4i5cFdK2QEMmIGU2AMnbsN8PJB13gpR8nlu6qU5UAgCkR4229u4CqeHZPR1SRS8X\ngGT0cgNI+nsLkEQ3r4EU/PTvfEilXZeCVPnE3AxSKEc5E5CmhGQ3NzWkdf+Fhhwq80+cDKnx\nx5KKJA2QKocC0nI1SIc9QJJJgPSo0djdgJCiP1GnSveGdFx+K0htQwFpuTEk4xPNrwASkHpC\nav2k7GuAukfX6CaSBEhViYDUCVLwk7LxIFV2oyDFPkIE0gIkK4DU0A2QhoZ0PPkFSDJR/pjG\ng5T8SRZIjxqN3QHJzlFO5ENK3t+mboB0BUi7AvWZBoLU9OxOCKl03DFIzlBAWiaBFMohMlVD\nqmtmXkh++1eAJF+krwlp3f5jUw6zm6tC+r1TXSAVN00NyX23ezFIu13zQEoaqM6SSvJySNU3\nAkhGAMnr5ixIpQGcboAEJC/F5mtvgNT2QVlFkipITUNdCVL+oHNxGqS6VoBU3w6QTodk99Yf\n0v6TsiikujT7bcl/akqhEklIsY8QgXQbSPt9da2MC2lNfC2xGkhAugIkM01m21shxW5MoTSQ\ngHSIZ0jmo+s3k8pbm0N1AyQgbXekn7qXQKp7dP1muj26s0A6HncQUvqh2cU9Ie0LVGcCUnqT\n9cxZ3bQkOhdSec8tIKkX6YtCWjvk8LqphVTbDJCA5ASQ6rvpAenvPr0Nkj3A2JDMHxuvBmmz\nbyZIKQP1WRJJmiBVDXUoACQnEZDSm3pDasqS+Oc6SC1DAWmpgdT2wNwAUuyTD/egmpMA6VHj\nubP3QEp+UjYCpJa7ki5dyFQ6JiB9xfUguZ0l9/SF5HbxAkiRFLoZIAFpZEjr+yGth6+kFgMJ\nSBeA5KY5bFwDOU6BdExa2Y7TjEyS+EgISI8aLZ2l9mQO+qWQvEfXaiYC6bjx+pAyx70A6adG\nS2epPZkUr4FU9eja3YwIqb6ZkSGtQDrcBCAVigNpn6dmqMkheS/SF4a0NufoDCn0ycc5kJ46\nq4TU8N0hB8kd4W6Qmr7zAsntJmGgIcvxn94AqfJm3AbS897aYz4fUuVzUttN8Zi6QPq0DkoF\nkESNhsaS0/SFZDdxGqQVSLlWFiAlajz1BaSnnYNBajuiRwUgvRiS2Vd6mmkg/W4FUrqV3T8B\n6VGjvq/0NGNBCjx3Q0Oq72VkSOYp3QpS7p5cF9JPdIBU6saBtI4DaXvcX+nbIFV8d7gjpGOB\nikwt33kz/x1I4Lk7D1Jy7e0g1b/M3gRS+d1uJaQ9CBmZ/zKxHyRvgnQGp5vyMY0DKfWzDZCe\natT3tZumA6TA83ICpK2kFkg13XiQ9jRrYgupOVPutN8IyRxidEjW5y9AEt28CVJ9mpEg1U0B\nJCPuBamhlT9INZ84F1rZ/307pNXs5e2Q7PhuK7KtMcU+U3Oa3cZQO8+OQhni3fzs7HR7I710\nO+39YUfuS2Wc/IrkfmcofVvIpbjiK9LzS5I7QTKB0431itR6Qo8K1kHJyCVpfUXyX2bf/opk\n1ajtKj1N9qSrIDU8MECS0RfS83E/0gNp6Q4pUaAmUzdIoQfvuZ2lB6RiN+KYJoRU/YPfXSAV\nP38x73Lz87KDFP2Md9/O0hFSeullIB0/bQPSc43artLjvB3StoOgo02CACSnmxpIDZ0sO0iB\nVJnTboZktwKk5bKQlr/9LZAqunkhpIpf3RRa2Sd5JyRvDiAZcRakpwotm7pBWhrvy3OFsSHp\n9akCQCpmqm/kTpBaGlm+BhgLUu1pXwVS+zn1g9T8jRdIMjpDWoGUrNHhc5zSvwoGpCCkpe2+\nPFfoDWnbD5AeNXpCShWoTAWkTKr2aZ6POJTrFEjG8vkhGb8RsO9ym6MDpOgvSxIVWjYlmmmF\nFB7mLEg/6QP/wRaQduMMBin0uBwqtGxKdJNLr/MHZ3n+MTh0Y4BUqjEKpMan/2/KeSEFozek\n3d2NQXJWFyB5+y8DqTBVxV1uaqMEqTJVrkLTrrEgGd/xjOgHqe4Hv4tAqmkpN04XSEtLF0CS\n0Q1S+jepbQPUNJIsMCwka6TcOPkUZz8mQJLxBCnyTfOQ5Sd94AZ5S4G0AEl2c1lIv+lfMkCu\nF2s/kHQcIAV/65io0LTr0MxIkFozAalQY0JI0edlV6Fp16GbbPrXPIddICV/k/peSNbH5xeA\nJH8j8DpInX7reKzQtAtI/eIakII3d2BIvSq0bZsd0l/690CqOGUgGQEkGQdIgVxAytbo9k0q\nU6A5sRfDQ1JP76uew86QntK/F5IxzTUgqU932xNbMTok+fReC1Li3gLpq0av71LZAoHEThwh\ndXleniu0bXObedlz2Om+HFK8C5L//fIqkMoFzkr8k39USO6vbq4GadlneDckPQ+QnPxPBU5x\nNA2kzu94/9IDaVn2/yNU/QNIw0Cq+dfbatK/G5Ic6RWQzn5BAtJAkE5K/7YBBoLU/ak7DnFa\n5u/814BUSA8kWSHzdfeYXwipZqy6eCWkXr913FVo3AekbhVyF8aB9PMzkj9UbdwXktkNkHSF\n7JXRINkj1QeQgBSukL80FiR3oJZ4KaRev77fVmjd6DUDJF0huv9FkIJtigKnZs9D6liheSeQ\nOlWI7n8FpNMf9LPzDwzJ+ixn/OdQpR9+ACA5+bcFTvgYMjKB0cz4z6FKP/wAQHLy7wr0/xQy\nNIFuZvznUKUffgAgOfmHP0aVngFkheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmH\nP0aVngFkheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFkheh+IBn5hz9G\nlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFkheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWe\nAWSF6H4gGfmHP0aVngFkheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFk\nheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFkheh+IBn5hz9GlZ4BZIXo\nfiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFkheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4g\nGfmHP0aVngFkheh+IBn5hz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFkheh+IBn5\nhz9GlZ4BZIXofiAZ+Yc/RpWeAWSF6H4gGfmHP0aVngFkheh+IBn5hz9GlZ4BZIXo/igkI/73\nv9NSvyYY4N1xpQGAlA0GeHdcaQAgZYMB3h1XGgBI2WCAd8eVBjgPEkHcKIBEEB0CSATRIYBE\nEB0CSATRIU6D9PHxcVbqU+PR9Xf32z8uEcnOGeD0OAvSx88jebH4PrdH99s/LhEfjxEY4NVx\nEqSP379cKx7f+b673/5xjXi0zwAvDyDt4tqQPq76HG7iigMAaRfXhvRffFzxOXyKj48rngCQ\ndnF5SNd8Dp/jkt8JgLSLq0P6ePrLJQf45K3dU1xh9HRcHNLH318vOcBlf8jj4+9dXPvj74+/\nP685wCcff+9i/F+hpeP7W/rFfh34Ex9f8XndATKdjz8A/4oQQXQIIBFEhwASQXQIIBFEhwAS\nQXQIIBFEhwASQXQIIBFEhwASQXQIIF0j/vmN//7+3d0Qh+BMrhFAGjw4k+sEgAYOzuY6AaSB\ng7O5TvxA+nl79/XP3198vOkj3hbc/evEDtKXnW9APz8+Ee8Kbv51Ygvp7y8/V5D0xuDeXyf2\nb+3+/gKktwf3/jpRgsR7uzcH9/46oV6RiDcGJ3CdANLAwQlcJ4pv7T7h9Nbg3l8n8pA++Rnp\n3cG9v04UIPEL2XcHd58gOgSQCKJDAIkgOgSQCKJDAIkgOgSQCKJDAIkgOgSQCKJDAIkgOgSQ\nCKJDAIkgOgSQCKJDAIkgOgSQCKJDAIkgOsT/AZwVIGMpxhnzAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ggplot(sirforced_out_long[sirforced_out_long$variable==\"I\",],aes(x=time,y=value,colour=variable,group=variable))+\n",
    "  # Add line\n",
    "  geom_line(lwd=2)+\n",
    "  #Add labels\n",
    "  xlab(\"Time\")+ylab(\"Number\")+\n",
    "  theme(legend.position=\"none\")"
   ]
  }
 ],
 "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
}