StackOverFlow Error in SIR model using ModelingTookit

substitute.(eqs, Dict([S => s*N, I => i*N, R => r*N]))

I tried that too but got this error:

ArgumentError: broadcasting over dictionaries and `NamedTuple`s is reserved

Stacktrace:
 [1] broadcastable(#unused#::Dict{Num, Num})
   @ Base.Broadcast .\broadcast.jl:683
 [2] broadcasted(::Function, ::Vector{Equation}, ::Dict{Num, Num})
   @ Base.Broadcast .\broadcast.jl:1313
 [3] top-level scope
   @ In[24]:3
 [4] eval
   @ .\boot.jl:360 [inlined]
 [5] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base .\loading.jl:1094

substitute.(eqs, (Dict([S => s*N, I => i*N, R => r*N]),))

2 Likes

Thanks! So I now have the following:

eqs = substitute.(eqs, (Dict([S => s*N, I => i*N, R => r*N]),))
eqs = expand_derivatives.(eqs)

s_eq = Symbolics.solve_for(eqs[1],D(s))
i_eq = Symbolics.solve_for(eqs[2],D(i))
r_eq = Symbolics.solve_for(eqs[3],D(r))

eqs_scaled = [D(s) ~ s_eq, D(i) ~ i_eq, D(r) ~ r_eq]

Which gives me:

\begin{align} \frac{ds(t)}{dt} =& - \beta i\left( t \right) s\left( t \right) \\ \frac{di(t)}{dt} =& - \frac{N \gamma i\left( t \right) - N \beta i\left( t \right) s\left( t \right)}{N} \\ \frac{dr(t)}{dt} =& \gamma i\left( t \right) \end{align}

Now for some reason it doesn’t simplify i(t) by canceling out the N's. I tried simplify(i_eq) but it just gives back the same thing…is there any way to fully simplify it?

Can you make an MWE out of that? I would assume that would further simplify.

I managed to figure it out! Here’s a MWE:

using ModelingToolkit, Symbolics

@parameters t β γ N
@variables s(t) i(t) r(t)

D = Differential(t)
eq = D(i) ~ (N*β*s*i - N*γ*i)/(N)     
eq = Symbolics.solve_for(eq,D(i))
Symbolics.simplify(eq)

produces

\begin{equation} \frac{N \beta i\left( t \right) s\left( t \right) - N \gamma i\left( t \right)}{N} \end{equation}

Changing the last line to Symbolics.simplify(eq;expand=true) makes it cancel out the N:
\begin{equation} \beta i\left( t \right) s\left( t \right) - \gamma i\left( t \right) \end{equation}

1 Like

@shashi is there a reason why it doesn’t expand that by default?

Almost there… here’s what I have so far:

@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) 

#1) Define the (unscaled equations)
D = Differential(t)
SIR_eqs = [D(S) ~ -β/N*S*I,
           D(I) ~ β/N*S*I - γ*I,
           D(R) ~ γ*I]

SIR = ODESystem(SIR_eqs, name = :SIR)

#2) Rewrite equations in terms of s = S/N, i = I/N, r = R/N
sir_eqs = substitute.(SIR_eqs, (Dict([S => s*N, I => i*N, R => r*N]),))
sir_eqs = expand_derivatives.(sir_eqs)

s_eq = Symbolics.solve_for(eqs_new[1],D(s))
i_eq = Symbolics.solve_for(eqs_new[2],D(i))
r_eq = Symbolics.solve_for(eqs_new[3],D(r))

#cancel out N 
i_eq = Symbolics.simplify(i_eq;expand=true) 

#Define the new scaled system
sir_eqs = [D(s) ~ s_eq, D(i) ~ i_eq, D(r) ~ r_eq]
SIR_scaled = ODESystem(eqs_new, name = :SIR_scaled) 

#3) Specify the initial condition and parameter values
u0 = [S => 0.99,
      I => 0.01,
      R => 0.00
     ]

#N omitted from `pars` b.c. it doesn't appear in the scaled system.
pars = [β => 0.20, γ => 0.15]   

#4) Scale the initial condition 
function scale_u0(u0::Vector{Pair{Num, Float64}}, x::Vector{Float64})
    #x = the scaling factors (constants)
    vars = first.(u0)
    vals = last.(u0)
    
    u0_scaled = [Pair(vars[i], x[i]*vals[i]) for i=1:length(vars)]
    return u0_scaled
end    

N = 1e3
u0_scaled = scale_u0(u0,repeat([N],3))

#4) Now create the ODEProblem and solve the scaled system
prob = ODEProblem(SIR_scaled,u0_scaled,(0.0,200.0),pars)
#sol = solve(prob,Tsit5())

But when I try to create the ODEProblem, it gives the following error:

ArgumentError: Term{Real, Nothing}[s(t), i(t), r(t)] are missing from the variable map.

Stacktrace:
  [1] throw_missingvars(vars::Vector{Term{Real, Nothing}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\variables.jl:71
  [2] _varmap_to_vars(varmap::Dict{Num, Float64}, varlist::Vector{Term{Real, Nothing}}; defaults::Dict{Any, Any}, check::Bool, toterm::typeof(Symbolics.diff2term))
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\variables.jl:63
  [3] varmap_to_vars(varmap::Vector{Pair{Num, Float64}}, varlist::Vector{Term{Real, Nothing}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function)
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\variables.jl:38
  [4] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, parammap::Vector{Pair{Num, Float64}}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:461
  [5] process_DEProblem
    @ C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:437 [inlined]
  [6] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Num, Float64}}; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:550
  [7] (ODEProblem{true, tType, isinplace, P, F, K, PT} where {tType, isinplace, P, F, K, PT})(sys::ODESystem, u0map::Vector{Pair{Num, Float64}}, tspan::Tuple{Float64, Float64}, parammap::Vector{Pair{Num, Float64}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:550
  [8] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any, N} where N; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:529
  [9] ODEProblem(::ODESystem, ::Vector{Pair{Num, Float64}}, ::Vararg{Any, N} where N)
    @ ModelingToolkit C:\Users\Michael\.julia\packages\ModelingToolkit\o7Ahw\src\systems\diffeqs\abstractodesystem.jl:529
 [10] top-level scope
    @ In[156]:46
 [11] eval
    @ .\boot.jl:360 [inlined]
 [12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base .\loading.jl:1094

A quick search brought this up: Missing from the variable map error - Usage - JuliaLang. Not sure if mine is exactly the same issue. I wasn’t able to figure it out…

I think you need to specify initial values u0 for lower case [s, i, r] not upper case since you did a substitution. Those are the values mentioned in the error message, the variable map is the list of initial values you provided.

2 Likes

Thanks! That was indeed the issue!

Alright, here is my completed little toy problem!

using ModelingToolkit, DifferentialEquations, Symbolics

@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) x y 

#1) Define the (unscaled equations)
D = Differential(t)
SIR_eqs = [D(S) ~ -β/N*S*I,
           D(I) ~ β/N*S*I - γ*I,
           D(R) ~ γ*I]

SIR = ODESystem(SIR_eqs, name = :SIR)

#2) Rewrite equations in terms of s = S/N, i = I/N, r = R/N
sir_eqs = substitute.(SIR_eqs, (Dict([S => s*N, I => i*N, R => r*N]),))
sir_eqs = expand_derivatives.(sir_eqs)

s_eq = Symbolics.solve_for(sir_eqs[1],D(s))
i_eq = Symbolics.solve_for(sir_eqs[2],D(i))
r_eq = Symbolics.solve_for(sir_eqs[3],D(r))

i_eq = Symbolics.simplify(i_eq;expand=true)   #cancel out the N 


sir_eqs = [D(s) ~ s_eq, D(i) ~ i_eq, D(r) ~ r_eq]
SIR_scaled = ODESystem(sir_eqs, name = :SIR_scaled)  #Define the new scaled system

#3) Specify the initial condition and parameter values
u0 = [S => 990.0,
      I => 10.0,
      R => 0.00
     ]

pars = [β => 0.20, γ => 0.15]   #N is omitted because it doesn't appear in the scaled system.


function scale_u0(u0::Vector{Pair{Num, Float64}},vars_new::Vector{Num},x::Vector{Float64})
    
    #x = the scaling factors for the variables
    vars = first.(u0)
    vals = last.(u0)
    
    u0_scaled = [Pair(vars_new[i], x[i]*vals[i]) for i=1:length(vars_new)]
    return u0_scaled
end    

N = 1e-3
u0_scaled = scale_u0(u0,[s,i,r],repeat([N],3))

#4) Now solve the ODE and plot it
prob = ODEProblem(SIR_scaled,u0_scaled,(0.0,200.0),pars)
sol = solve(prob,Tsit5(), saveat = 1:1:200)  #works!

Thanks for the help, @ChrisRackauckas! (And also @contradict!).

One follow-up question: Suppose I want to write a function or two that generalizes the above process of scaling the variables to any ODE system whose equations are polynomials in the dependent variables. Ideally I’d like some function that takes in a list of ODE equations, and a vector of variable-scale factor pairs (e.g. [S => 1e3, I => 1e3, R => 1e3], though the scaling factors don’t need to be the same for each). Then I’d want the function to return a new ODE system written in terms of the scaled variables. Is this possible the way that ModelingToolkit is written? I something like this could ever be added to the MKT package, that would really be awesome.

You can add observed variables, which are just scalings of other state variables, and the system will automatically eliminate them from the solve. You then just do sol[r_scaled] or whatnot and it’ll give you those values in the solution.

Hmm, I’m not sure how to implement that…Do you think you could illustrate this with a simple example? I’ve checked this page for more info but this particular topic doesn’t have any code…

I’m not sure how to add them back after the fact, but you can do this:

using ModelingToolkit, DifferentialEquations
using Plots

@parameters t β γ N
@variables S(t) I(t) R(t) s(t) i(t) r(t) x y
D = Differential(t)

SIR_eqs = [D(s) ~ -β*s*i,
           D(i) ~ β*s*i - γ*i,
           D(r) ~ γ*i,
            0 ~ s*N - S,
            0 ~ i*N - I,
            0 ~ r*N - R,
]

SIR = ODESystem(SIR_eqs, name = :SIR)

# After this, the system has only 3 states [s,i,r]
SIR = structural_simplify(SIR)

u0 = [s => 0.99, i => 0.01, r => 0.00]
pars = [β => 0.20, γ => 0.15, N=>1e3]

prob = ODEProblem(SIR, u0, (0.0, 200.0), pars)

sol = solve(prob; saveat=1:1:200)

To plot scaled

plot(sol; vars=[s,i,r])

Or to plot unscaled

plot(sol; vars=[S,I,R])
2 Likes

@contradict Thanks, this helps a lot! Let me see if I understand correctly: by saying

SIR = structural_simplify(SIR)

the variables s,i,r are the variables which are actually being integrated, whereas as the variables S,I,R are the “observables” which are calculated from the solution for s(t),i(t),r(t)?

Yes, in this case that is what happens. structural_simplify is kind of a magic box that applies several tools to try to make a simpler system to solve. As such, it isn’t very controllable so you sometimes have to play with the input equations if you want a particular shape. The fun starts here if you want to see what it is doing.

2 Likes

Awesome, this looks to be exactly what I originally wanted!

Another question: How does structural_simplify “know” that s,i,r are the observables? Is it the 0 on the LHS of the equations? Or the fact that the s,i,r equations don’t involve derivatives? Or did we just happen to get lucky in this case and structural_simplify somehow inferred what the observables would be?

…I’ll check out the source code in more detail (though based on a quick glance, it seems a bit over my head, currently :slight_smile: ).

It uses a graph algorithm to identify the minimal equation to solve and solves that, generating explicit representations of the other equations from the elimination equations.

2 Likes

Hmm, very interesting!

I have some follow-up questions on this. I think I’ll probably start a new thread though since it doesn’t really have to do with my original title.

Also, would it be ok if I made a PR for the documentation to include the change of variables example with the SIR equations?

1 Like

That would be nice.