SIR on networks with LightGraphs

New to Julia
1

I am creating the simplest and lowest level SIR simulation with LightGraphs and ModelingToolkit (I am a beginner). For now, I have a simple network of 10 nodes and 7 edges. And it is a directed graph.

Here is the top part of my code up until the line that crashes:

using ModelingToolkit
using DiffEqBase
using DiffEqJump
using LightGraphs

nedges = 7
nverts = 10
dg = DiGraph(nverts, nedges)

# Allow for different beta on each edge
@parameters t β[1:nedges]  γ[1:nedges] ;
@variables S[1:nverts](t);
@variables I[1:nverts](t);
@variables R[1:nverts](t);

rxsS = [Reaction(β[i],[I[src(e)], S[dst(e)]], [S[dst(e)]])
    for (i,e) ∈ enumerate(edges(dg))]

rxsI = [Reaction(γ[i],[I[src(e)]], [R[src(e)]])  # note: src to src, yet there is no edge
    for (i,e) ∈ enumerate(edges(dg))]

rxs = vcat(rxsS, rxsI)
rs = ReactionSystem(rxs, t, u, k)
js = convert(JumpSystem, rs)

S0 = Array{Float64,1}(undef, nverts)
I0 = Array{Float64,1}(undef, nverts)
R0 = Array{Float64,1}(undef, nverts)
fill!(S0, 1.);
fill!(I0, 0.);
fill!(R0, 0.);
S0[1] = 0.; # One person is infected
I0[1] = 1.;
R0[1] = 0.;

# Two column vectors
γ = Array{Float64,1}(undef, nedges)'
β = Array{Float64,1}(undef, nedges)'
fill!(γ, 0.25)
fill!(β, 0.50)
# Row vector
pmap = hcat(β,γ)
u0 = vcat(S0,I0,R0)

tspan = (0.0,10.0)
dprob = DiscreteProblem(js, u0, tspan, pmap)

DiscreteProblem throws an error:

ERROR: MethodError: Cannot `convert` an object of type 
  Float64 to an object of type 
  Variable
Closest candidates are:
  convert(::Type{Variable}, ::Operation) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/ModelingToolkit.jl:72
  convert(::Type{T}, ::T) where T at essentials.jl:171
  Variable(::Any) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/variables.jl:50
  ...
Stacktrace:
 [1] varmap_to_vars(::Array{Float64,1}, ::Array{Variable,1}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/variables.jl:251
 [2] DiscreteProblem(::JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}, ::Array{Float64,1}, ::Tuple{Float64,Float64}, ::LinearAlgebra.Adjoint{Float64,Array{Float64,1}}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/jumps/jumpsystem.jl:143
 [3] DiscreteProblem(::JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}, ::Array{Float64,1}, ::Tuple{Float64,Float64}, ::LinearAlgebra.Adjoint{Float64,Array{Float64,1}}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/jumps/jumpsystem.jl:143
 [4] top-level scope at none:0

which I assume is due to the way I combine my reactions:

rxs = vcat(rxsS, rxsI)

Alternatively, it may have to do with my unconventional reactions equations, which I am not sure are allowed:

rxsS = [Reaction(β[i],[I[src(e)], S[dst(e)]], [S[dst(e)]])
    for (i,e) ∈ enumerate(edges(dg))]

rxsI = [Reaction(γ[i],[I[src(e)]], [R[src(e)]])  # note: src to src, yet there is no edge
    for (i,e) ∈ enumerate(edges(dg))]

Notice that rxsS and rxsI have reactions defined as

I(src) + S(dst) --> I(dst)
I(src) --> I(src)

which is not the usual way of writing reaction equations on a network, but the ones above reflect what is really happening. The infected person at a node eventually recovers at that same node. There is interaction with another node.

Any help would be greatly appreciated.

2

I don’t have anything topical to add, but I wanted to let you know you can simplify your code.

You don’t need to preallocate the vectors if you’re filling them:

S0 = Array{Float64,1}(undef, nverts)
I0 = Array{Float64,1}(undef, nverts)
R0 = Array{Float64,1}(undef, nverts)
fill!(S0, 1.);
fill!(I0, 0.);
fill!(R0, 0.);

becomes

S0 = ones(nverts)  # Float64 is implied
I0 = zeros(nverts)
R0 = zeros(nverts)

and then

# Two column vectors
γ = Array{Float64,1}(undef, nedges)'
β = Array{Float64,1}(undef, nedges)'
fill!(γ, 0.25)
fill!(β, 0.50)

becomes, simply,

γ = fill(0.25, nedges)'
β = fill(0.50, nedges)'
1 Like
3

Thank you very much, Seth! Regarding the use of fill! versus zeros, I was aware, but I did some timings and concluded that my approach, while heavier, was more efficient for large arrays. Doing my timings this morning, that does not appear to be the case. I will change my code to take the more elegant approach.

4

The way we decide which approach to take in LightGraphs is as follows: if we’re going to reuse the vector many times (that is, if we need the same type of vector for many iterations), we’ll reuse it (whether we preallocate or not is dependent on what the first iteration needs) by calling fill! on it at the end of each iteration. If we need the vector once inside a function, there’s no benefit to preallocating.

1 Like
5 
using ModelingToolkit
using DiffEqBase
using DiffEqJump
using LightGraphs
using SpecialGraphs

nedges = 7
nverts = 10

#=
# Simpler than using callbacks. Different results though
sir_model = @reaction_network sir_rn begin
  β/tot_pop, s + i --> 2i
  γ, i --> r
end tot_pop β γ
=#

dg = DiGraph(nverts, nedges)

# Allow for different beta on each edge
@parameters t β[1:nedges]  γ[1:nedges] ;
@variables S[1:nverts](t);
@variables I[1:nverts](t);
@variables R[1:nverts](t);
rxsS = [Reaction(β[i],[I[src(e)], S[dst(e)]], [S[dst(e)]])
    for (i,e) ∈ enumerate(edges(dg))]

rxsI = [Reaction(γ[i],[I[src(e)]], [R[src(e)]])  # note: src to src, yet there is no edge
    for (i,e) ∈ enumerate(edges(dg))]

rxs = vcat(rxsS, rxsI)
vars = vcat(S,I,R)

rs = ReactionSystem(rxs, t, vars, [β,γ])
js = convert(JumpSystem, rs)

I get the following error when calling ReactionSystem and have no idea what to do about it:

julia> rs = ReactionSystem(rxs, t, vars, [β,γ])
ERROR: MethodError: Cannot `convert` an object of type 
  Array{Operation,1} to an object of type 
  Variable
Closest candidates are:
  convert(::Type{Variable}, ::Operation) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/ModelingToolkit.jl:72
  convert(::Type{T}, ::T) where T at essentials.jl:171
  Variable(::Any) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/variables.jl:50
  ...
Stacktrace:
 [1] _broadcast_getindex_evalf at ./broadcast.jl:631 [inlined]
 [2] _broadcast_getindex at ./broadcast.jl:614 [inlined]
 [3] getindex at ./broadcast.jl:564 [inlined]
 [4] copy at ./broadcast.jl:854 [inlined]
 [5] materialize at ./broadcast.jl:820 [inlined]
 [6] ReactionSystem(::Array{Reaction{Variable,Int64},1}, ::Operation, ::Array{Operation,1}, ::Array{Array{Operation,1},1}; systems::Array{ReactionSystem,1}, name::Symbol) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/reaction/reactionsystem.jl:141
 [7] ReactionSystem(::Array{Reaction{Variable,Int64},1}, ::Operation, ::Array{Operation,1}, ::Array{Array{Operation,1},1}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/reaction/reactionsystem.jl:141
 [8] top-level scope at none:0

Any help would be much appreciated. The documentation could provide an example of a system with multiple variables, where each variable is an array, as it is in this case. Thanks.

6 
rs = ReactionSystem(rxs, t, vars, [β;γ])

β is an array of variables, so you need to concatente it rather than make an array of arrays. This is something we’re improving though.

7

Very good observation, Chris. So now, I have an error in DiscreteEquations. Another conversion problem. Here is some code:

using ModelingToolkit
using DiffEqBase
using DiffEqJump
using LightGraphs
using SpecialGraphs

nedges = 7
nverts = 10
#dg = SpecialGraphs.CompleteGraph(nverts);
#nedges = ne(dg)
dg = DiGraph(nverts, nedges)

# Allow for different beta on each edge
@parameters t β[1:nedges]  γ[1:nedges] ;
@variables S[1:nverts](t);
@variables I[1:nverts](t);
@variables R[1:nverts](t);

#@parameters t k[1:nedges];
#@variables u[1:nverts](t);

# Specify one or more Reaction equations per edge
# u's are numbered 1 through nedges
# I need to create one list for S, one list for R, one list for I,
#  and concatenate the lists

# Two people interact: Node A (Infeced) and node B (Susceptible)
# Node B might get infected
# I, S, R will be either 1 or 0 at each node.

rxsS = [Reaction(β[i],[I[src(e)], S[dst(e)]], [S[dst(e)]])
    for (i,e) ∈ enumerate(edges(dg))]

rxsI = [Reaction(γ[i],[I[src(e)]], [R[src(e)]])  # note: src to src, yet there is no edge
    for (i,e) ∈ enumerate(edges(dg))]

rxs = vcat(rxsS, rxsI);
vars = vcat(S,I,R);
params = vcat(β,γ);
rs = ReactionSystem(rxs, t, vars, params);
js = convert(JumpSystem, rs);

S0 = ones(nverts)
I0 = zeros(nverts)
R0 = zeros(nverts)
S0[1] = 0.; # One person is infected
I0[1] = 1.;
R0[1] = 0.;
vars = vcat(S0, I0, R0);

# Two column vectors
γ = fill(0.25, nedges);
β = fill(0.50, nedges);
params = vcat(β,γ);
pmap = hcat(β,γ)
u0 = vcat(S0,I0,R0)

# I create two dictionaries; one for initial conditions, and 
# another for the problem parameters. Quite inconvenient 
# and wasteful of resources if the problem were very large. 

u0 = Dict(convert(Variable,state) => vars[i] for (i,state) in enumerate(states(js)));
params0 = Dict(convert(Variable,par) => params[i] for (i,par) in enumerate(parameters(js)));

tspan = (0.0,10.0)
dprob = DiscreteProblem(js, u0, tspan, params0)

I print out params0 and u0 (I guess the order is immaterial since the data is encoded in a dictionary):

julia> params0
Dict{Variable{ModelingToolkit.Parameter{Number}},Float64} with 14 entries:
  β₆ => 0.5
  β₇ => 0.5
  β₂ => 0.5
  γ₂ => 0.25
  γ₃ => 0.25
  γ₇ => 0.25
  γ₆ => 0.25
  β₅ => 0.5
  ⋮  => ⋮

and

julia> u0
Dict{Variable{Number},Float64} with 30 entries:
  R₄  => 0.0
  I₄  => 0.0
  I₁  => 1.0
  I₁₀ => 0.0
  I₇  => 0.0
  R₁  => 0.0
  R₆  => 0.0
  S₄  => 1.0
  S₆  => 1.0
  I₉  => 0.0
  I₅  => 0.0
  R₃  => 0.0
  S₉  => 1.0
  S₈  => 1.0
  I₃  => 0.0
  S₇  => 1.0
  S₃  => 1.0
  R₇  => 0.0
  S₅  => 1.0

Here is the error (still related to convert):

dprob = DiscreteProblem(js, u0, tspan, params0)
ERROR: MethodError: no method matching similar(::Dict{Variable{Number},Float64}, ::Type{Float64})
Closest candidates are:
  similar(::JuliaInterpreter.Compiled, ::Any) at /Users/erlebach/.julia/packages/JuliaInterpreter/owKmt/src/types.jl:7
  similar(::Array{T,1}, ::Type) where T at array.jl:358
  similar(::Array{T,2}, ::Type) where T at array.jl:359
  ...
Stacktrace:
 [1] varmap_to_vars(::Dict{Variable{Number},Float64}, ::Array{Variable,1}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/variables.jl:249
 [2] DiscreteProblem(::JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}, ::Dict{Variable{Number},Float64}, ::Tuple{Float64,Float64}, ::Dict{Variable{ModelingToolkit.Parameter{Number}},Float64}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/jumps/jumpsystem.jl:143
 [3] DiscreteProblem(::JumpSystem{RecursiveArrayTools.ArrayPartition{DiffEqJump.AbstractJump,Tuple{Array{MassActionJump,1},Array{ConstantRateJump,1},Array{VariableRateJump,1}}}}, ::Dict{Variable{Number},Float64}, ::Tuple{Float64,Float64}, ::Dict{Variable{ModelingToolkit.Parameter{Number}},Float64}) at /Users/erlebach/.julia/packages/ModelingToolkit/MohAp/src/systems/jumps/jumpsystem.jl:143
 [4] top-level scope at none:0

I solved the problem: I replaced the dictionary by an Array of Pairs.

Question: What is the functional difference between a dictionary and an array of pairs? Can one access elements in an array of pairs via their name with the efficiency of a dictionary ?

When does one use an array of pairs and when does one use a dictionary. I come from Python, and before that, C++ and Fortran.

Thanks.

8

Looks like this was solved in Lightgraph + ModelToolkit optimization .