search

epirecipes / sir-julia

Various implementations of the classical SIR model in Julia
MIT License
195 stars 39 forks source link

Add matrix exponential formulation #86

Closed sdwfrost closed 1 year ago

sdwfrost commented 1 year ago

See Keeling and Ross and Sherlock

ChrisRackauckas commented 1 year ago

https://github.com/kaandocal/FiniteStateProjection.jl

sdwfrost commented 1 year ago

Thanks @ChrisRackauckas! I hadn't seen this. Any idea why the following fails at the last step? Is it because I'm using a ReactionSystem rather than a reaction network?

using Catalyst
using FiniteStateProjection
using OrdinaryDiffEq

@parameters t β c γ
@variables S(t) I(t) R(t)

N=S+I+R
rxs = [Reaction((β*c)/N, [S,I], [I], [1,1], [2])
       Reaction(γ, [I], [R])]
@named rs  = ReactionSystem(rxs, t, [S,I,R], [β,c,γ])

fspsys = FSPSystem(rs)
u0 = zeros(3, 1001)
u0[1,991] = 1.0
u0[2,11] = 1.0
p = [0.05, 10.0, 0.25]
tspan = (0.0, 40.0)
fspprob = convert(ODEProblem, fspsys, u0, tspan, p)
fspsol = solve(fspprob, Tsit5())

Output (on solve):

ERROR: BoundsError: attempt to access Tuple{Int64, Int64} at index [3]
Stacktrace:
  [1] getindex
    @ ./tuple.jl:29 [inlined]
  [2] getindex(index::CartesianIndex{2}, i::Int64)
    @ Base.IteratorsMD ./multidimensional.jl:93
  [3] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/6v5Gn/src/RuntimeGeneratedFunctions.jl:137 [inlined]
  [4] macro expansion
    @ ./none:0 [inlined]
  [5] generated_callfunc
    @ ./none:0 [inlined]
  [6] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)})(::Matrix{Float64}, ::Matrix{Float64}, ::Vector{Float64}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/6v5Gn/src/RuntimeGeneratedFunctions.jl:124
  [7] (::ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing})(::Matrix{Float64}, ::Vararg{Any})
    @ SciMLBase ~/.julia/packages/SciMLBase/ys6dl/src/scimlfunctions.jl:2096
  [8] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Matrix{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Matrix{Float64}}, ODESolution{Float64, 3, Vector{Matrix{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Matrix{Float64}}, Vector{Float64}, Vector{Vector{Matrix{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Matrix{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/4OfcV/src/perform_step/low_order_rk_perform_step.jl:766
  [9] __init(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/4OfcV/src/solve.jl:499
 [10] __init(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/4OfcV/src/solve.jl:73
 [11] #__solve#582
    @ ~/.julia/packages/OrdinaryDiffEq/4OfcV/src/solve.jl:5 [inlined]
 [12] __solve
    @ ~/.julia/packages/OrdinaryDiffEq/4OfcV/src/solve.jl:5 [inlined]
 [13] #solve_call#22
    @ ~/.julia/packages/DiffEqBase/egmnd/src/solve.jl:494 [inlined]
 [14] solve_call
    @ ~/.julia/packages/DiffEqBase/egmnd/src/solve.jl:466 [inlined]
 [15] #solve_up#29
    @ ~/.julia/packages/DiffEqBase/egmnd/src/solve.jl:915 [inlined]
 [16] solve_up
    @ ~/.julia/packages/DiffEqBase/egmnd/src/solve.jl:890 [inlined]
 [17] #solve#27
    @ ~/.julia/packages/DiffEqBase/egmnd/src/solve.jl:825 [inlined]
 [18] solve(prob::ODEProblem{Matrix{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.FullSpecialize, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:du, :u, :p, :t), FiniteStateProjection.var"#_RGF_ModTag", FiniteStateProjection.var"#_RGF_ModTag", (0x8ae605f3, 0xb893c5c5, 0xdfddb792, 0xe99de15c, 0x0cc083e6)}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/egmnd/src/solve.jl:817
 [19] top-level scope
    @ Untitled-1:21
ChrisRackauckas commented 1 year ago

I am not sure.

sdwfrost commented 1 year ago

I've mostly worked it out. I should try to replicate Keeling and Ross.

sdwfrost commented 1 year ago

Added FiniteStateProjection.jl example with 2f1b191f88facb892accc88d3ec012544f5b9513. Closing for now.