Agent Based Model Augmentation

SemanticModels supports model augmentation, which is the derivation of new models from old models with different (usually more advanced) capabilities. The approach is to define a type to represent a class of models and then a set of transformations that can act on that type to change the capabilities of the model.

We can apply our model augmentation framework to models that are not defined as an analytical mathematical expression. A widely used class of models for complex systems are agent based in that they have an explicit representation of the agents with states and functions to represent their behavior and interactions. This notebook examines how to apply model transformations to augment agent based simulations.

We are going to use the simulation in examples/agentbased.jl as a baseline simulation and add capabilities to the simulation with SemanticModels transformations. The simulation in question is an implementation of a basic SIRS model on a static population. We will make two augmentations.

  1. Add un estado de los muertos or a state for the dead, transforming the model from SIRS to SIRD
  2. Add vital dynamics a represented by net population growth

These changes to the model could easily be made by changing the source code to add the features. However, this notebook shows how those changes could be scripted by a scientist. As we all know, once you can automate a scientific task by introducing a new technology, you free the mind of the scientist for more productive thoughts.

In this case we are automating the implementation of model changes to free the scientist to think about what augmentations to the model should I make? instead of how do I implement these augmentations?

In [1]:
using SemanticModels
In [2]:
using SemanticModels.Parsers
using SemanticModels.ModelTools
using SemanticModels.ModelTools.ExpStateModels
import Base: push!
In [3]:
samples = 7
nsteps = 10
finalcounts = Any[]
Out[3]:
0-element Array{Any,1}
In [4]:
println("Running Agent Based Simulation Augmentation Demo")
println("================================================")
println("demo parameters:\n\tsamples=$samples\n\tnsteps=$nsteps")
Running Agent Based Simulation Augmentation Demo
================================================
demo parameters:
	samples=7
	nsteps=10

Baseline SIRS model

Here is the baseline model, which is read in from a text file. You could instead of using parsefile use a quote/end block to code up the baseline model in this script.

Agents progress from the susceptible state to infected and then recovered, and get become susceptible again after recovery. See the file ../examples/agentbased.jl for a full description of this model

In [5]:
expr = parsefile("agentbased.jl")
m = model(ExpStateModel, expr)
#ModelTools.funclines(m.expr, :main)
Out[5]:
ExpStateModel(
  states=:([:S, :I, :R]),
  agents=Expr[:(a = sm.agents), :(a = fill(:S, n))],
  transitions=Expr[:(T = Dict(:S => (x...->begin
                      #= none:120 =#
                      if rand(Float64) < stateload(x[1], :I)
                          :I
                      else
                          :S
                      end
                  end), :I => (x...->begin
                      #= none:121 =#
                      if rand(Float64) < ρ
                          :I
                      else
                          :R
                      end
                  end), :R => (x...->begin
                      #= none:122 =#
                      if rand(Float64) < μ
                          :R
                      else
                          :S
                      end
                  end)))]
)
In [6]:
println("\nRunning basic model")
AgentModels = eval(m.expr)
for i in 1:samples
    newsam, counts = AgentModels.main(nsteps)
    push!(finalcounts, (model=:basic, counts=counts))
end
Running basic model
newsam.agents = Symbol[:R, :S, :S, :S, :R, :S, :R, :S, :I, :S, :S, :S, :I, :S, :S, :I, :S, :S, :R, :R]
newsam.agents = Symbol[:S, :R, :R, :S, :I, :R, :S, :S, :I, :R, :R, :S, :S, :R, :S, :S, :S, :S, :S, :S]
newsam.agents = Symbol[:I, :I, :R, :S, :R, :I, :R, :R, :I, :R, :R, :I, :R, :R, :I, :S, :I, :I, :S, :S]
newsam.agents = Symbol[:S, :I, :R, :I, :S, :S, :R, :I, :R, :S, :I, :R, :S, :I, :S, :S, :S, :I, :R, :R]
newsam.agents = Symbol[:R, :R, :S, :S, :S, :S, :S, :S, :R, :S, :I, :S, :S, :S, :I, :S, :I, :S, :I, :S]
newsam.agents = Symbol[:S, :S, :I, :R, :S, :S, :S, :I, :S, :I, :S, :R, :R, :S, :S, :S, :S, :R, :R, :R]
newsam.agents = Symbol[:I, :R, :S, :I, :I, :S, :S, :I, :S, :R, :I, :S, :R, :S, :S, :I, :R, :S, :I, :I]
In [7]:
m
Out[7]:
ExpStateModel(
  states=:([:S, :I, :R]),
  agents=Expr[:(a = sm.agents), :(a = fill(:S, n))],
  transitions=Expr[:(T = Dict(:S => (x...->begin
                      #= none:120 =#
                      if rand(Float64) < stateload(x[1], :I)
                          :I
                      else
                          :S
                      end
                  end), :I => (x...->begin
                      #= none:121 =#
                      if rand(Float64) < ρ
                          :I
                      else
                          :R
                      end
                  end), :R => (x...->begin
                      #= none:122 =#
                      if rand(Float64) < μ
                          :R
                      else
                          :S
                      end
                  end)))]
)

Adding the Dead State

We are going to add an additional state to the model to represent the infectious disease fatalities. The user must specify what that concept means in terms of the name for the new state and the behavior of that state. D is a terminal state for a finite automata.

In [8]:
println("\nThe system states are $(m.states.args)")
println("\nAdding un estado de los muertos")

put!(m, ExpStateTransition(:D, :((x...)->:D)))

println("\nThe system states are $(m.states.args)")
# once you are dead, you are dead forever
println("\nThere is no resurrection in this model")
println("\nInfected individuals recover or die in one step")

# add a transition rule for infected -> recovered, dead, or infected
m[:I] = :((x...)->begin
        roll = mod(rand(Int),3)
        if roll == 1
            return :R
        elseif roll == 2
            return :D
        else
            return :I
        end
    end
)
@show m[:I]
The system states are Any[:(:S), :(:I), :(:R)]

Adding un estado de los muertos

The system states are Any[:(:S), :(:I), :(:R), :(:D)]

There is no resurrection in this model

Infected individuals recover or die in one step
m[:I] = :(x...->begin
          #= In[8]:12 =#
          #= In[8]:13 =#
          roll = mod(rand(Int), 3)
          #= In[8]:14 =#
          if roll == 1
              #= In[8]:15 =#
              return :R
          elseif #= In[8]:16 =# roll == 2
              #= In[8]:17 =#
              return :D
          else
              #= In[8]:19 =#
              return :I
          end
      end)
Out[8]:
:(x...->begin
          #= In[8]:12 =#
          #= In[8]:13 =#
          roll = mod(rand(Int), 3)
          #= In[8]:14 =#
          if roll == 1
              #= In[8]:15 =#
              return :R
          elseif #= In[8]:16 =# roll == 2
              #= In[8]:17 =#
              return :D
          else
              #= In[8]:19 =#
              return :I
          end
      end)
In [9]:
println("\nRunning SIRD model")
AgentModels = eval(m.expr)
for i in 1:samples
    newsam, counts = AgentModels.main(nsteps)
    push!(finalcounts, (model=:sird, counts=counts))
end
Running SIRD model
newsam.agents = Symbol[:D, :I, :I, :D, :I, :D, :D, :S, :I, :D, :D, :D, :I, :R, :D, :I, :I, :R, :S, :I]
newsam.agents = Symbol[:S, :S, :S, :D, :D, :D, :D, :S, :S, :S, :D, :D, :D, :D, :S, :S, :D, :R, :S, :D]
newsam.agents = Symbol[:D, :S, :R, :S, :D, :D, :S, :I, :D, :D, :R, :D, :D, :I, :D, :S, :D, :D, :S, :D]
newsam.agents = Symbol[:R, :I, :S, :D, :S, :I, :D, :S, :D, :D, :S, :S, :R, :I, :D, :D, :I, :S, :D, :R]
newsam.agents = Symbol[:D, :D, :S, :D, :R, :D, :D, :D, :D, :D, :D, :D, :D, :S, :D, :S, :S, :R, :D, :S]
newsam.agents = Symbol[:D, :S, :I, :R, :S, :D, :S, :S, :I, :S, :D, :I, :D, :S, :D, :D, :I, :D, :S, :S]
newsam.agents = Symbol[:D, :D, :D, :S, :S, :S, :S, :D, :S, :S, :D, :S, :D, :D, :D, :D, :D, :D, :S, :D]
WARNING: replacing module AgentModels.

Some utilities for manipulating functions at a higher level than expressions.

In [10]:
struct Func end

function push!(::Func, func::Expr, ex::Expr)
    push!(bodyblock(func), ex)
end
Out[10]:
push! (generic function with 26 methods)

Population Growth

Another change we can make to our model is the introduction of population growth. Our model for population is that on each timestep, one new suceptible person will be added to the list of agents. We use the tick! function as an anchor point for this transformation.

In [11]:
println("\nAdding population growth to this model")
stepr = filter(x->isa(x,Expr), findfunc(m.expr, :tick!))[1]
@show stepr
push!(Func(), stepr, :(push!(sm.agents, :S)))
println("------------------------")
@show stepr;
Adding population growth to this model
stepr = :(function tick!(sm::StateModel)
      #= none:59 =#
      sm.loads = map((s->begin
                      #= none:59 =#
                      stateload(sm, s)
                  end), sm.states)
      #= none:60 =#
      return sm.loads
  end)
------------------------
stepr = :(function tick!(sm::StateModel)
      #= none:59 =#
      sm.loads = map((s->begin
                      #= none:59 =#
                      stateload(sm, s)
                  end), sm.states)
      #= none:60 =#
      return sm.loads
      push!(sm.agents, :S)
  end)
In [12]:
println("\nRunning growth model")
AgentModels = eval(m.expr)
for i in 1:samples
    newsam, counts = AgentModels.main(nsteps)
    push!(finalcounts, (model=:growth, counts=counts))
end
Running growth model
newsam.agents = Symbol[:I, :D, :S, :D, :I, :D, :I, :D, :D, :S, :I, :D, :S, :D, :D, :D, :D, :R, :R, :D]
newsam.agents = Symbol[:D, :S, :D, :R, :D, :D, :R, :S, :S, :D, :D, :R, :D, :D, :D, :S, :S, :S, :D, :D]
newsam.agents = Symbol[:S, :D, :I, :S, :R, :S, :I, :R, :R, :I, :D, :S, :S, :D, :D, :S, :I, :D, :I, :D]
newsam.agents = Symbol[:S, :S, :D, :D, :R, :S, :D, :D, :S, :D, :S, :D, :D, :D, :S, :D, :S, :D, :S, :D]
newsam.agents = Symbol[:I, :S, :D, :S, :D, :S, :S, :D, :S, :I, :D, :D, :D, :D, :S, :D, :D, :D, :S, :D]
newsam.agents = Symbol[:D, :S, :D, :S, :D, :D, :S, :S, :D, :D, :I, :S, :S, :D, :S, :D, :S, :S, :D, :R]
newsam.agents = Symbol[:D, :S, :D, :S, :S, :D, :I, :D, :D, :D, :D, :S, :S, :I, :D, :I, :S, :S, :S, :D]
WARNING: replacing module AgentModels.

Presentation of results

We have accumulated all of our simulation runs into the list finalcounts we process those simulation runs into summary tables describing the results of those simulations. This table can be used to make decisions and drive further inquiry.

In [13]:
println("\nModel\t Counts")
println("-----\t ------")
for result in finalcounts
    println("$(result.model)\t$(result.counts)")
end

function groupagg(x::Vector{Tuple{S,T}}) where {S,T}
    c = Dict{S, Tuple{Int, T}}()
    # c2 = Dict{S, T}()
    for r in x
        g = first(r)
        c[g] = get(c, g,(0, 0.0)) .+ (1, last(r))
    end
    return c
end

mean_healthy_frac = [(r.model,
                  map(last, filter(x->(x.first == :R || x.first == :S), r.counts))[1] / sum(map(last, r.counts))[1])
                 for r in finalcounts] |> groupagg

num_unhealthy = [(r.model,
                  map(last,
                      sum(map(last, filter(x->(x.first != :R && x.first != :S),
                             r.counts)))))
                 for r in finalcounts] |> groupagg

println("\nModel\t Count \t Num Unhealthy \t Mean Healthy %")
println("-----\t ------\t --------------\t  --------------")
for (g, v) in mean_healthy_frac
    μ = last(v)/first(v)
    μ′ = round(μ*100, sigdigits=5)
    x = round(last(num_unhealthy[g]) / first(num_unhealthy[g]), sigdigits=5)
    println("$g\t   $(first(v))\t  $(rpad(x, 6))\t   $(μ′)")
end
Model	 Counts
-----	 ------
basic	Pair{Symbol,Int64}[:S=>12, :I=>3, :R=>5]
basic	Pair{Symbol,Int64}[:S=>12, :I=>2, :R=>6]
basic	Pair{Symbol,Int64}[:S=>4, :I=>8, :R=>8]
basic	Pair{Symbol,Int64}[:S=>8, :I=>6, :R=>6]
basic	Pair{Symbol,Int64}[:S=>13, :I=>4, :R=>3]
basic	Pair{Symbol,Int64}[:S=>11, :I=>3, :R=>6]
basic	Pair{Symbol,Int64}[:S=>8, :I=>8, :R=>4]
sird	Pair{Symbol,Int64}[:S=>2, :I=>8, :R=>2, :D=>8]
sird	Pair{Symbol,Int64}[:S=>9, :I=>0, :R=>1, :D=>10]
sird	Pair{Symbol,Int64}[:S=>5, :I=>2, :R=>2, :D=>11]
sird	Pair{Symbol,Int64}[:S=>6, :I=>4, :R=>3, :D=>7]
sird	Pair{Symbol,Int64}[:S=>5, :I=>0, :R=>2, :D=>13]
sird	Pair{Symbol,Int64}[:S=>8, :I=>4, :R=>1, :D=>7]
sird	Pair{Symbol,Int64}[:S=>8, :I=>0, :R=>0, :D=>12]
growth	Pair{Symbol,Int64}[:S=>3, :I=>4, :R=>2, :D=>11]
growth	Pair{Symbol,Int64}[:S=>6, :I=>0, :R=>3, :D=>11]
growth	Pair{Symbol,Int64}[:S=>6, :I=>5, :R=>3, :D=>6]
growth	Pair{Symbol,Int64}[:S=>8, :I=>0, :R=>1, :D=>11]
growth	Pair{Symbol,Int64}[:S=>7, :I=>2, :R=>0, :D=>11]
growth	Pair{Symbol,Int64}[:S=>9, :I=>1, :R=>1, :D=>9]
growth	Pair{Symbol,Int64}[:S=>8, :I=>3, :R=>0, :D=>9]

Model	 Count 	 Num Unhealthy 	 Mean Healthy %
-----	 ------	 --------------	  --------------
growth	   7	  11.857	   33.571
sird	   7	  12.286	   30.714
basic	   7	  4.8571	   48.571