augustinas1 / MomentClosure.jl

Tools to generate and study moment equations for any chemical reaction network using various moment closure approximations
https://augustinas1.github.io/MomentClosure.jl/dev/
MIT License
44 stars 6 forks source link

Failing example model #2

Closed sdwfrost closed 3 years ago

sdwfrost commented 3 years ago

Hi, I'm trying to write an example for my SIR epidemiological model repository:

https://gist.github.com/sdwfrost/8394921f50fd6b4da0e77a72ac3f96b5

I've written the model with both polynomial and nonpolynomial rates, but I get the following errors:

For the polynomial model:

ocprob1 = ODEProblem(closed_central_eqs1, u0map1, tspan, p1)
ERROR: ArgumentError: Term{Real}[M₂₀₀(t)] are missing from the variable map.
Stacktrace:
 [1] varmap_to_vars(::Dict{Term{Real},Any}, ::Array{Term{Real},1}; defaults::Dict{Any,Any}) at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/variables.jl:18
 [2] #varmap_to_vars#11 at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/variables.jl:28 [inlined]
 [3] process_DEProblem(::Type{T} where T, ::ODESystem, ::Array{Pair{Term{Real},Any},1}, ::Array{Pair{Num,Float64},1}; 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{}}}) at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/systems/diffeqs/abstractodesystem.jl:277
 [4] process_DEProblem at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/systems/diffeqs/abstractodesystem.jl:271 [inlined]
 [5] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Term{Real},Any},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Num,Float64},1}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/systems/diffeqs/abstractodesystem.jl:323
 [6] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Pair{Term{Real},Any},1}, ::Tuple{Float64,Float64}, ::Array{Pair{Num,Float64},1}) at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/systems/diffeqs/abstractodesystem.jl:323
 [7] ODEProblem(::ODESystem, ::Array{Pair{Term{Real},Any},1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/systems/diffeqs/abstractodesystem.jl:303
 [8] ODEProblem(::ODESystem, ::Array{Pair{Term{Real},Any},1}, ::Vararg{Any,N} where N) at /home/simon/.julia/packages/ModelingToolkit/Df3sc/src/systems/diffeqs/abstractodesystem.jl:303
 [9] ODEProblem(::MomentClosure.ClosedMomentEquations, ::Array{Pair{Term{Real},Any},1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/simon/.julia/packages/MomentClosure/BRP8O/src/moment_equations.jl:89
 [10] ODEProblem(::MomentClosure.ClosedMomentEquations, ::Array{Pair{Term{Real},Any},1}, ::Vararg{Any,N} where N) at /home/simon/.julia/packages/MomentClosure/BRP8O/src/moment_equations.jl:89
 [11] top-level scope at none:1

For the non-polynomial model:

ocsol2 = solve(ocprob2,Tsit5())
ERROR: StackOverflowError:
Stacktrace:
 [1] recursive_unitless_bottom_eltype(::Type{Any}) at /home/simon/.julia/packages/RecursiveArrayTools/69f7F/src/utils.jl:86 (repeats 79984 times)
augustinas1 commented 3 years ago

Hi @sdwfrost, thanks for the issue!

The first error you get was due to a bug I overlooked in the code, should be fixed now. The second error you got is due to a variable map being passed to deterministic_IC: only a simple vector of initial values is supported, e.g. u0 = [990.0, 10.0, 0.0]. The function itself should be enforcing this constraint correctly now. Please let me know if you run into any other problems!

sdwfrost commented 3 years ago

Thanks! Runs great now. Updated gist too to do some plotting - not great plotting mean +/- 2se for these parameter sets, but I'll work on some further updates before committing to my repository.