intro_abm_jl_0.master.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "In this worksheet we are going to compare an analytical version of an SI model (Susceptible-Infected, a very simple epidemiological model) with an agent-based solution. This will show why it is in some situations necessary to use agent-based models while at the same time demonstrating the implementation details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An analytical SI model\n",
    "\n",
    "The SI (susceptible-infected) model is a simple epidemiological model with a [long tradition](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). It is based on the assumption that the modelled population is large and homogeneous enough for stochasticity, individual differences and waiting times between events to become irrelevant. This makes it possible to treat the healthy and infected portions of the population as \"compartments\" with infections and recoveries as \"flows\" between them which can be modelled using some [very well understood mathematics](https://en.wikipedia.org/wiki/Multi-compartment_model).\n",
    "\n",
    "We implement the model using the [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) package. We start by defining the update function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we need this for the solution\n",
    "using DifferentialEquations\n",
    "\n",
    "# f updates the derivative du dependent on the \n",
    "# current state u using parameters p \n",
    "# TODO what about t?\n",
    "function f(du, u, p, t) \n",
    "    # parameters\n",
    "    β, γ = p\n",
    "    # susceptible\n",
    "    du[1] = -β * u[2]*u[1] +   # infection\n",
    "             γ * u[2]          # recovery\n",
    "    # infected\n",
    "    du[2] = β * u[2]*u[1] -    # infection\n",
    "            γ * u[2]           # recovery\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To numerically solve the system we have to define boundary conditions and parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = (0.1, 0.03)       # parameters beta, gamma\n",
    "u0 = [0.99, 0.001]    # initial conditions\n",
    "tspan = (1.0, 500.0)  # simulation time\n",
    "# prepare the problem\n",
    "problem = ODEProblem(f, u0, tspan, p)\n",
    "# and solve it\n",
    "deqsolution = solve(problem);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DifferentialEquations provides some infrastructure to make plotting a solution (using the Plots package) very straightforward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots; gr() # Using the Plotly Backend\n",
    "Plots.plot(deqsolution, xlabel=\"time\", ylabel=\"N\", lab=[\"susceptible\" \"infected\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to try different parameter values we have to solve again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = (0.1, 0.06)       # parameters beta, gamma\n",
    "u0 = [0.99, 0.001]    # initial conditions\n",
    "tspan = (1.0, 500.0)  # simulation time\n",
    "# prepare the problem\n",
    "problem = ODEProblem(f, u0, tspan, p);\n",
    "# and solve it\n",
    "deqsolution = solve(problem);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the two solutions by plotting them in the same graph. We do this by simply calling `plot!` instead of `plot` (by convention the '!' indicates that a function modifies its arguments, in this case the previous graph)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Plots.plot!(deqsolution, xlabel=\"time\", ylabel=\"N\", lab=[\"susceptible2\" \"infected2\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An agent-based version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now going to implement an agent-based version of the same model. We will start off with a very simplistic implementation that is as close to the analytical model as possible and then add complications that will clarify the advantages of using an ABM."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prelude\n",
    "\n",
    "Some modules we are going to need down the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Type declarations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to tell Julia about how the objects look like that we are going to use. \n",
    "\n",
    "An `@enum` is a type that can only have a handful of named values, e.g.:\n",
    "```Julia\n",
    "@enum WeekDay Monday Tuesday Wednesday Thursday Friday\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all possible states a person can be in\n",
    "@enum Status susceptible infected\n",
    "\n",
    "# this is our agent type\n",
    "mutable struct SimplePerson\n",
    "    # state\n",
    "    status :: Status\n",
    "    # other agents this one can infect or be infected by\n",
    "    contacts :: Vector{SimplePerson}\n",
    "end\n",
    "\n",
    "# how we construct a person object\n",
    "SimplePerson() = SimplePerson(susceptible, [])   # default Person is susceptible and has no contacts\n",
    "SimplePerson(state) = SimplePerson(state, [])  # default Person has no contacts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation itself is mainly used for bookkeeping. We might want to switch out the Agent type later, therefore we keep it as a parameter to the Simulation type (using a so-called [parametric type](https://docs.julialang.org/en/v1/manual/types/#Parametric-Types))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is a parametric type\n",
    "# we can specify which type AGENT is going to be replaced with\n",
    "# when constructing our Simulation\n",
    "mutable struct Simulation{AGENT}\n",
    "    # model parameters:\n",
    "    # infection rate\n",
    "    inf :: Float64\n",
    "    # recovery rate\n",
    "    rec :: Float64\n",
    "\n",
    "    # and this is our population of agents\n",
    "    pop :: Vector{AGENT}\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we implement the update functions that change the agents' state. In this version of the model we use simple updating in time steps.\n",
    "\n",
    "Note that none of these include type information and can therefore be used for any agent type with the right struct members."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# update function for susceptible individuals\n",
    "function update_susceptible!(person, sim)\n",
    "    if length(person.contacts) == 0\n",
    "        return\n",
    "    end\n",
    "    # if random contact is infected check for transmission\n",
    "    if rand(person.contacts).status == infected && rand() < sim.inf\n",
    "        person.status = infected\n",
    "    end\n",
    "end\n",
    "\n",
    "# update function for infected individuals\n",
    "function update_infected!(person, sim)\n",
    "    if rand() < sim.rec\n",
    "        person.status = susceptible\n",
    "    end\n",
    "end\n",
    "\n",
    "# update an agent\n",
    "function update!(person, sim)\n",
    "    if person.status == susceptible\n",
    "        update_susceptible!(person, sim)\n",
    "    else\n",
    "        update_infected!(person, sim)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have to call these functions for every agent in our simulation at every time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function update_agents!(sim)\n",
    "    # we need to change the order, otherwise agents at the beginning of the \n",
    "    # pop array will behave differently from those further in the back\n",
    "    order = shuffle(sim.pop)\n",
    "    \n",
    "    for p in order\n",
    "        update!(p, sim)\n",
    "    end\n",
    "end   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With that out of the way we can create the simulation and run it. The only thing left to do is to setup our population. \n",
    "\n",
    "Note that if we *only* wanted to do a fully mixed population there would be no need to explicitly track contacts per agent and to set up the population like this. However, we might also look at situations where agents know only part of the population therefore we do it this way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up a mixed population\n",
    "# p_contact is the probability that two agents are connected\n",
    "function setup_mixed(n, p_contact)\n",
    "    pop = [ SimplePerson() for i=1:n ]\n",
    "\n",
    "    # go through all combinations of agents and \n",
    "    # check if they are connected\n",
    "    for i in eachindex(pop)\n",
    "        for j in i+1:length(pop)\n",
    "            if rand() < p_contact\n",
    "                push!(pop[i].contacts, pop[j])\n",
    "                push!(pop[j].contacts, pop[i])\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "\n",
    "    pop\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point we can already run the simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "source": [
    "<span style='color:Red'>See if you can complete the code below on your own. You should be able to produce a spreading infection.</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO add arguments\n",
    "#pop = setup_mixed(...) \n",
    "\n",
    "# is there anything else we need to do here?\n",
    "\n",
    "# TODO add arguments\n",
    "#sim = Simulation(...)\n",
    "\n",
    "for i in 1:20\n",
    "    # TODO update sim here\n",
    "    \n",
    "#    println(count(p -> p.status == infected, sim.pop)) # prints current number of infected\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "outputs": [],
   "source": [
    "pop = setup_mixed(1000, 1.0) \n",
    "\n",
    "pop[1].status = infected\n",
    "pop[2].status = infected\n",
    "pop[3].status = infected\n",
    "\n",
    "sim = Simulation(0.2, 0.1, pop) \n",
    "\n",
    "for i in 1:20\n",
    "    update_agents!(sim)\n",
    "    println(count(p -> p.status == infected, sim.pop)) # prints current number of infected\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make things a bit easier down the line I have wrapped both, the setup and the model run into their own separate functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#                     function parameters after ; have to be \n",
    "#                     named on call \n",
    "function  setup_sim(;p_i, p_r, N, p_c, n_inf, seed)\n",
    "    # for reproducibility\n",
    "    Random.seed!(seed)\n",
    "\n",
    "    # create a population of agents, fully mixed\n",
    "    pop = setup_mixed(N, p_c)\n",
    "\n",
    "    # create a simulation object with parameter values\n",
    "    sim = Simulation(p_i, p_r, pop)\n",
    "    for i in 1:n_inf\n",
    "        # one percent of agents are infected\n",
    "        sim.pop[i].status = infected\n",
    "    end\n",
    "    \n",
    "    sim\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function run_sim(sim, n_steps, verbose = false)\n",
    "    # we keep track of the numbers\n",
    "    n_inf = Int[]\n",
    "    n_susc = Int[]\n",
    "\n",
    "    # simulation steps\n",
    "    for t in  1:n_steps\n",
    "        update_agents!(sim)\n",
    "        push!(n_inf, count(p -> p.status == infected, sim.pop))\n",
    "        push!(n_susc, count(p -> p.status == susceptible, sim.pop))\n",
    "        # a bit of output\n",
    "        if verbose\n",
    "            println(t, \", \", n_inf[end], \", \", n_susc[end])\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    # return the results (normalized by pop size)\n",
    "    n = length(sim.pop)\n",
    "    n_inf./n, n_susc./n\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With that in place running the model and displaying the result is quite easy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = setup_sim(p_i=0.5, p_r=0.3, N=1000, p_c=1.0, n_inf=10, seed=42)\n",
    "# inf and susc are arrays containing the proportion of infected/susceptible over time\n",
    "inf, susc = run_sim(sim, 500)\n",
    "\n",
    "Plots.plot([susc, inf], labels = [\"susc (ABM)\" \"inf (ABM)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison with analytical model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to compare the ABM and the analytical model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "source": [
    "<span style='color:Red'> Let's run the model using the same parameter values as we used for the analytical version and plot both in the same graph.</span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO add code to run the ABM with the same parameters as the diff eq and plot it\n",
    "# code goes here\n",
    "\n",
    "\n",
    "p = (0.1, 0.03)       # parameters beta, gamma\n",
    "u0 = [0.99, 0.001]    # initial conditions\n",
    "tspan = (1.0, 500.0)  # simulation time\n",
    "problem = ODEProblem(f, u0, tspan, p)\n",
    "deqsolution = solve(problem);\n",
    "\n",
    "#Plots.plot!(deqsolution, xlabel=\"time\", ylabel=\"N\", lab=[\"susceptible\" \"infected\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "outputs": [],
   "source": [
    "sim = setup_sim(p_i=0.1, p_r=0.03, N=1000, p_c=1.0, n_inf=10, seed=42)\n",
    "inf, susc = run_sim(sim, 500)\n",
    "\n",
    "Plots.plot([susc, inf], labels = [\"susc (ABM)\" \"inf (ABM)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "source": [
    "<span style='color:Blue'> That does look pretty similar to the previous results, but to make sure we can plot both sets of curves in the same graphs. </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "outputs": [],
   "source": [
    "p = (0.1, 0.03)       # parameters beta, gamma\n",
    "u0 = [0.99, 0.001]    # initial conditions\n",
    "tspan = (1.0, 500.0)  # simulation time\n",
    "# prepare the problem\n",
    "problem = ODEProblem(f, u0, tspan, p);\n",
    "# and solve it\n",
    "deqsolution = solve(problem);\n",
    "\n",
    "Plots.plot!(deqsolution, xlabel=\"time\", ylabel=\"N\", lab=[\"susceptible\" \"infected\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "source": [
    "<span style='color:Blue'> As we can see the ABM produces very similar results to the analytical model. There is still a difference in the speed with which both models converge, however.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "source": [
    " <span style='color:Red'> Can you think of a reason for the difference in results?\n",
    " </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "source": [
    "<span style='color:Blue'> The reason for the difference between the analytical and the first agent-based version of the model lies in the update rules we used: For the sake of simplicity we updated the agent population in place. If we instead updated the population based on a complete copy of the previous state the curves of the two models would overlap fully (save for stochastic fluctuations in the ABM).</span>\n",
    "\n",
    "<span style='color:Blue'> This is a somewhat important point. In models with stepwise updating update rules can sometimes make a huge difference. When building a new (time stepped) model it is therefore always a good idea to think about what is going to be updated in which order and based on which information and how that relates to reality. </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytical as well as the first agent-based version of the model assumed that our population is well-mixed, i.e. that any two individuals/agents meet with the same probability. This is obviously unrealistic as soon as we go beyond the level of a single social group. Does it matter, though?\n",
    "\n",
    "Let's implement a very simple spatial version of the model and see if we get different results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is our agent type\n",
    "# in principle we could have re-used the Person type from before, however,\n",
    "# it is much easier to produce a nice visualization of the model if we\n",
    "# keep the coordinates in the agents (as we will se below)\n",
    "mutable struct SpacePerson\n",
    "    # state\n",
    "    status :: Status\n",
    "    # other agents this one can infect or be infected by\n",
    "    contacts :: Vector{SpacePerson}\n",
    "    # location in the world\n",
    "    x :: Float64\n",
    "    y :: Float64\n",
    "end\n",
    "\n",
    "# how we construct a person object\n",
    "SpacePerson(x, y) = SpacePerson(susceptible, [], x, y)   # default Person is susceptible and has no contacts\n",
    "SpacePerson(state, x, y) = SpacePerson(state, [], x, y)  # default Person has no contacts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We kept the agent type as a parameter to the Simulation type therefore we don't need to change it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since our update functions (see above) didn't contain any type information we actually don't have to do anything here (yay generic programming!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spatial setup is a tad fiddly, therefore I put it in a separate file (take a look if you are interested) and include it here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"setup_world.jl\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sets up the simulation and the agents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function setup_space_sim(; p_i, p_r, xs, ys, n_inf, seed)\n",
    "    # create a population of agents living on a \n",
    "    # 50x50 grid (von Neumann)\n",
    "    pop = setup_grid(SpacePerson, xs, ys)\n",
    "    \n",
    "    # create a simulation object with parameter values\n",
    "    sim = Simulation(p_i, p_r, pop)\n",
    "\n",
    "    # for reproducibility\n",
    "    Random.seed!(seed)\n",
    "\n",
    "    # pick n_inf random agents and set them to infected\n",
    "    for p in rand(sim.pop, n_inf)\n",
    "        p.status = infected\n",
    "    end\n",
    "    \n",
    "    sim\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The run function hasn't changed either, so we are done with the preliminaries and can proceed with running the simulation and plotting it against the simple, non-spatial version. \n",
    "\n",
    "This time we have to specify the size of the world (as `xs` and `ys`) instead of the number of agents, otherwise the function calls are identical. As before we try to keep the parameters as close to the analytical version as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "source": [
    "<span style='color:Red'> Try to run the spatial model using the same parameters as the simple version. Plot the results of both in the same graph to make them easier to compare. </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO run spatial model here\n",
    "\n",
    "# Plots.plot([susc, inf], labels = [\"susc\" \"inf\"])\n",
    "\n",
    "# TODO run simple model here\n",
    "\n",
    "# Plots.plot!([susc, inf], labels = [\"susc (mixed)\" \"inf (mixed)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "outputs": [],
   "source": [
    "sim = setup_space_sim(p_i=0.1, p_r=0.03, xs=50, ys=50, n_inf=10, seed=42)\n",
    "inf, susc = run_sim(sim, 500)\n",
    "\n",
    "Plots.plot([susc, inf], labels = [\"susc\" \"inf\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "source": [
    "<span style='color:Blue'> Clearly this is not in any way similar to what happened before. Just for confirmation let's give it a try. </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "outputs": [],
   "source": [
    "sim = setup_sim(p_i=0.1, p_r=0.03, N=2500, p_c=1.0, n_inf=10, seed=42)\n",
    "inf, susc = run_sim(sim, 500)\n",
    "\n",
    "Plots.plot!([susc, inf], labels = [\"susc (mixed)\" \"inf (mixed)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "source": [
    "<span style='color:Blue'> Yes, that's pretty much as expected. </span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But wait a moment. In the spatial model we drastically reduced the number of contacts per agent, didn't we? What happens if we do the same in the simple model?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "source": [
    "<span style='color:Red'> Re-run the simple model with reduced numbers of contacts and plot the results. </span>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "questions"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO run low-contact version of the simple model here\n",
    "\n",
    "# Plots.plot!([susc, inf], labels = [\"susc (low contacts)\" \"inf (low contacts)\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "outputs": [],
   "source": [
    "sim = setup_sim(p_i=0.1, p_r=0.03, N=2500, p_c=4/2500, n_inf=10, seed=42)\n",
    "\n",
    "inf, susc = run_sim(sim, 500)\n",
    "\n",
    "Plots.plot!([susc, inf], labels = [\"susc (low contacts)\" \"inf (low contacts)\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "solutions"
    ]
   },
   "source": [
    "<span style='color:Blue'> As we can see the changed behaviour in the spatial model is not caused by the lower number of contacts alone, but by the spatial structure itself. Why this is the case will become clear once we see the output of the visualization below.</span>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is generally a good idea to have a proper visualisation for a simulation. It makes it a lot easier to catch bugs, over time it provides an intuition of the system's behaviour and it is fun.\n",
    "\n",
    "Usually it would be good to have a dynamic display of the system state over time, but for this simple example we will make do with static output. If you want to track the progress of the simulation you can just execute `runsim(sim, 500)` again and then rerun the visualisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this uses the Luxor package, which allows for very straighforward \n",
    "# generation of SVG (PDF, etc.) drawings\n",
    "using Luxor\n",
    "\n",
    "sim = setup_space_sim(p_i=0.1, p_r=0.03, xs=50, ys=50, n_inf=10, seed=42)\n",
    "inf, susc = run_sim(sim, 50)\n",
    "\n",
    "@svg begin\n",
    "    # scale from model coordinates to screen coordinates\n",
    "    f_scale = 11\n",
    "    # size of agents\n",
    "    p_size = 5\n",
    "    \n",
    "    # otherwise we'll only have to lower right quarter\n",
    "    # (don't ask...)\n",
    "    origin(0, 0)\n",
    "\n",
    "    # paint it black\n",
    "    sethue(\"black\")\n",
    "    \n",
    "    # ** first draw all connections\n",
    "    # (we don't want to worry about agent size here,\n",
    "    #  so we just draw the lines from x1,y1 to x2,y2\n",
    "    #  and paint the agents on top of that)\n",
    "    \n",
    "    # for each agent\n",
    "    for p in sim.pop\n",
    "        x1 = p.x * f_scale\n",
    "        y1 = p.y * f_scale\n",
    "        \n",
    "        # draw all lines to other agents\n",
    "        for p2 in p.contacts\n",
    "            x2 = p2.x * f_scale\n",
    "            y2 = p2.y * f_scale\n",
    "            \n",
    "            setline(1)\n",
    "            setdash(\"solid\")\n",
    "            line(Point(x1, y1), Point(x2, y2), :stroke)\n",
    "        end\n",
    "    end\n",
    "    # yes, this draws all lines twice\n",
    "    # who cares\n",
    "\n",
    "    # ** now draw all agents\n",
    "    \n",
    "    for p in sim.pop\n",
    "        x = p.x * f_scale\n",
    "        y = p.y * f_scale\n",
    "        \n",
    "        # show the status as colours\n",
    "        if p.status == infected\n",
    "            sethue(\"red\")\n",
    "        elseif p.status == susceptible\n",
    "            sethue(\"blue\")\n",
    "        end\n",
    "        \n",
    "        circle(Point(x, y), p_size, :fill)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Things to try\n",
    "\n",
    "* Run for longer, is there an equilibrium or an end state?\n",
    "* Try different values for `inf` and `rec`.\n",
    "* Does changing the shape of the world (`x` vs. `y` size) have an effect?"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  },
  "toc-autonumbering": true,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 4
}