# SIR on networks with LightGraphs

**URL:** <https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497>\
**Category:** New to Julia\
**Created:** [May 31, 2020, 2:17am UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497 "2020-05-31T02:17:33Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [May 31, 2020, 2:17am UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/1 "2020-05-31T02:17:33Z")

</div>

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:

```julia
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:

```julia
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:

```julia
rxs = vcat(rxsS, rxsI)

```

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

```julia
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

```julia
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.

---

<div class="post-metadata">

**Author:** ![anon94023334](https://avatars.discourse-cdn.com/v4/letter/a/e274bd/32.png) [@anon94023334](https://discourse.julialang.org/u/anon94023334)\
**Post date:** [May 31, 2020, 4:14am UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/2 "2020-05-31T04:14:47Z")

</div>

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 `fill`ing them:

```julia
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

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

```

and then

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

```

becomes, simply,

```julia
γ = fill(0.25, nedges)'
β = fill(0.50, nedges)'

```

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [May 31, 2020, 2:45pm UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/3 "2020-05-31T14:45:05Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![anon94023334](https://avatars.discourse-cdn.com/v4/letter/a/e274bd/32.png) [@anon94023334](https://discourse.julialang.org/u/anon94023334)\
**Post date:** [May 31, 2020, 3:34pm UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/4 "2020-05-31T15:34:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [June 3, 2020, 2:29pm UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/5 "2020-06-03T14:29:36Z")

</div>

```julia
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
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 3, 2020, 2:37pm UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/6 "2020-06-03T14:37:18Z")

</div>

```julia
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.

---

<div class="post-metadata">

**Author:** ![erlebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/erlebach/32/12973_2.png) [@erlebach](https://discourse.julialang.org/u/erlebach)\
**Post date:** [June 3, 2020, 3:17pm UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/7 "2020-06-03T15:17:24Z")

</div>

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

```julia
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
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
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):

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [June 4, 2020, 12:51am UTC](https://discourse.julialang.org/t/sir-on-networks-with-lightgraphs/40497/8 "2020-06-04T00:51:15Z")

</div>

Looks like this was solved in [Lightgraph + ModelToolkit optimization](https://discourse.julialang.org/t/lightgraph-modeltoolkit-optimization/40696) .