Open sdwfrost opened 3 years ago
Hi @sdwfrost
Let me answer this first.
I'm also not sure how to feed the initial state
You cannot because the StaticSystem
is used to model systems that are not dynamical. StaticSystem
models systems whose input output relation is of the form
y = g(u, t)
where t
, u
, y
is the time, input and output of the system. So there is not dynamics here. However, the SIR model you simulated here using the FunctionMap
of the DifferentialEquations.jl
is dynamical system represented by difference equations. So, instead of @def_static_system
, you should use @def_discrete_system
to define a DiscreteSystem
whose input-state-output relation is given by
x[k + 1] = f(x[k], u[k], k)
y[k] = g(x[k], u[k], k)
where k
, u
, x
, y
is the time, input, state and output of the system. Below is the case when def_discrete_system
is used.
using Causal
using Plots
@inline function rate_to_proportion(r::Float64,t::Float64)
1-exp(-r*t)
end
function sir_map!(dx, x, u, t, p=nothing)
S, I, R = x
β, c, γ, δt = p
N = S + I + R
infection = rate_to_proportion(β * c * I / N, δt) * S
recovery = rate_to_proportion(γ, δt) * I
@inbounds begin
dx[1] = S - infection
dx[2] = I + infection - recovery
dx[3] = R + recovery
end
nothing
end
@def_discrete_system mutable struct SIRSystem{RH, RO, ST, IP, OP} <: AbstractDiscreteSystem
β::Float64 = 0.05
c::Float64 = 10.
γ::Float64 = 0.25
δ::Float64 = 0.1
righthandside::RH = (dx, x, u, t, p=(β, c, γ, δt)) -> sir_map!(dx, x, u, t, p)
readout::RO = (x, u, t) -> x
state::ST = [990., 10., 0.]
input::IP = nothing
output::OP = Outport(3)
end
@defmodel model begin
@nodes begin
sirsys = SIRSystem()
writer = Writer(input=Inport(3))
end
@branches begin
sirsys => writer
end
end
δt = 0.1
nsteps = 400
tmax = nsteps * δt
sim = simulate!(model, 0., δt, tmax)
t, x = read(getnode(model, :writer).component)
plot(t, x[:, 1], label="S")
plot!(t, x[:, 2], label="I")
plot!(t, x[:, 3], label="R")
Above, the whole SIR model (with all there state variables S, I, R), is used. And I think, this is not what you want. You want to divide the system into sub systems as infection system (with state variables S, I) and recovery system (with state variables I, R) and connect them together. From your equations I could not get how the whole SIR dynamics is divided into two. I could not decouple the variables (S, I, R) into the state variables (S, I) and (I, R) so that two sub-dynamical systems for infection and recovery dynamics can be constructed. For example, from
function infection(t,u,p)
(S,I) = u
(β,γ,δt) = p
ifrac = rate_to_proportion(β*I,δt)
inf = ifrac*S
return [S-inf,I+inf]
end
I understand the state variables and output of the infection system is (S, I) and (S - inf, I + inf), respectively. From
function recovery(t,u,p)
(S,I) = u
(β,γ,δt) = p
rfrac = rate_to_proportion(γ,δt)
rec = rfrac*I
return [S,I-rec]
end
it seems that the state variables and output of the recovery system is (S, I) and (S, I - rec), respectively. (By the way, I think, here the state variable S should be renamed to R to reflect recovery). And from
infsys[1:2] => recsys[1:2]
recsys[1:2] => infsys[1:2]
it seems that the systems infsys
and recsys
drives each other. But, following the signal flow between infsys
and recsys
, I could not reconstruct the whole SIR dynamics which should be like
function sir_map!(du,u,p,t)
(S,I,R) = u
(β,c,γ,δt) = p
N = S+I+R
infection = rate_to_proportion(β*c*I/N,δt)*S
recovery = rate_to_proportion(γ,δt)*I
@inbounds begin
du[1] = S-infection
du[2] = I+infection-recovery
du[3] = R+recovery
end
nothing
end;
Hi @zekeriyasari
Thanks for your detailed response. Perhaps I can clear up some vagueness in my question!
So, back to my original code; I don't want to specify an initial state (as the nodes take input and produce output in a single time step), I want to specify an initial input to get things started.
I first tried a DiscreteSystem, but as the inputs determine the (initial) state of the system, I was using functions like this, where my state is overwritten by the inputs.
function infection(dx, x, u, t, p)
(S,I) = u
x = u
(β,γ,δt) = p
ifrac = rate_to_proportion(β*I,δt)
inf = ifrac*S
dx[1] = S-inf
dx[2] = I+inf
end
function recovery(dx, x, u, t, p)
(S,I) = u
x = u
(β,γ,δt) = p
rfrac = rate_to_proportion(γ,δt)
rec = rfrac*I
dx[1] = S
dx[2] = I-rec
end
Here's the full example using DiscreteSystem:
using Causal
@inline function rate_to_proportion(r::Float64,t::Float64)
1-exp(-r*t)
end
function infection(dx, x, u, t, p)
(S,I) = u
x = u
(β,γ,δt) = p
ifrac = rate_to_proportion(β*I,δt)
inf = ifrac*S
dx[1] = S-inf
dx[2] = I+inf
end
ofunc(x, u, t) = x
@def_discrete_system mutable struct InfectionSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
β::Float64 = 0.5/1000.0
γ::Float64 = 0.25
δt::Float64 = 0.1
p = (β,γ,δt)
righthandside::RH = (dx, x, u, t, p) -> infection(dx, x, u, t, p)
state::Vector{Float64} = [990.0,10.0]
readout::RO = ofunc
input::IP = Inport(2)
output::OP = Outport(2)
end
@def_discrete_system mutable struct RecoverySystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
β::Float64 = 0.5/1000.0
γ::Float64 = 0.25
δt::Float64 = 0.1
p = (β,γ,δt)
righthandside::RH = (dx, x, u, t, p) -> recovery(dx, x, u, t, p)
state::Vector{Float64} = [990.0,10.0]
readout::RO = ofunc
input::IP = Inport(2)
output::OP = Outport(2)
end
@defmodel model begin
@nodes begin
infsys = InfectionSystem()
recsys = RecoverySystem()
writer = Writer(input=Inport(2))
end
@branches begin
infsys[1:2] => recsys[1:2]
recsys[1:2] => infsys[1:2]
infsys[1:2] => writer[1:2]
end
end
simulate!(model, 0.0, 0.1, 0.1)
When running the above model, I get the following output (and Julia hangs there):
[ Info: 2021-04-17T07:39:20.862 Started simulation...
[ Info: 2021-04-17T07:39:21.212 Inspecting model...
┌ Info: The model has algrebraic loops:[[1, 2]]
└ Trying to break these loops...
[ Info: Loop [1, 2] is broken
[ Info: 2021-04-17T07:39:23.078 Done.
[ Info: 2021-04-17T07:39:23.078 Initializing the model...
[ Info: 2021-04-17T07:39:23.591 Done...
[ Info: 2021-04-17T07:39:25.724 Running the simulation...
┌ Warning: `top` is deprecated, use `first` instead.
│ caller = modify_dt_for_tstops!(integrator::OrdinaryDiffEq.ODEIntegrator{OrdinaryDiffEq.FunctionMap{false}, true, Vector{Float64}, Nothing, Float64, Interpolant{Buffer{Cyclic, Float64, 1}, Buffer{Cyclic, Float64, 2}, Vector{Interpolations.Extrapolation{Float64, 1, Interpolations.ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Interpolations.Line{Nothing}}}}, Float64, Float64, Float64, Vector{Vector{Float64}}, DiffEqBase.ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, DiffEqBase.DiscreteProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Interpolant{Buffer{Cyclic, Float64, 1}, Buffer{Cyclic, Float64, 2}, Vector{Interpolations.Extrapolation{Float64, 1, Interpolations.ScaledInterpolation{Float64, 1, Interpolations.BSplineInterpolation{Float64, 1, OffsetArrays.OffsetVector{Float64, Vector{Float64}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line{Interpolations.OnGrid}}}, Interpolations.Line{Nothing}}}}, DiffEqBase.DiscreteFunction{true, var"#18#23", Nothing, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, OrdinaryDiffEq.FunctionMap{false}, OrdinaryDiffEq.InterpolationData{DiffEqBase.DiscreteFunction{true, var"#18#23", Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.FunctionMapCache{Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats}, DiffEqBase.DiscreteFunction{true, var"#18#23", Nothing, Nothing}, OrdinaryDiffEq.FunctionMapCache{Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DEOptions{Bool, Bool, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), DiffEqBase.CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Bool, Nothing, OrdinaryDiffEq.DefaultInit}) at integrator_utils.jl:46
└ @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/LQQYm/src/integrators/integrator_utils.jl:46
Hi @sdwfrost
The function signatures for the right-hand-side of a DiscreteSystem
must be of the form
function righthandside(dx, x, u, t)
# update dx
end
So, to be used in a DiscreteSystem
, the signatures of infection
and recovery
functions should be
function infection(dx, x, u, t, p=nothing)
(S,I) = u
x = u
(β,γ,δt) = p
ifrac = rate_to_proportion(β*I,δt)
inf = ifrac*S
dx[1] = S-inf
dx[2] = I+inf
end
function recovery(dx, x, u, t, p=nothing)
(S,I) = u
x = u
(β,γ,δt) = p
rfrac = rate_to_proportion(γ,δt)
rec = rfrac*I
dx[1] = S
dx[2] = I-rec
end
The parameter container p
must be provided as default argument. Moreover, the input u
is considered as a vector of functions, i.e., u(t)= [u_1(t), u_2(t)]
. The type of u
is Causal.Interpolant
and there is no iterate
method to be applicable for the type Causal.Interpolant.
, So you must explicitly overwrite the state variables (S, I).
With these modifications, I ran the code for a time interval of [0., 40.] seconds. However, the time waveforms of S
and I
differ from the ones you obtained here In order for the time waveforms to settle to their final values, we need a simulation of approximately 100 seconds. I could not understand why, unfortunately, but I wanted to inform you of the current status. I suspect that the inconsistency may come from the algebraic loop-breaking algorithm but I am not a hundred percent sure.
using Causal
using Plots
@inline function rate_to_proportion(r::Float64,t::Float64)
1-exp(-r*t)
end
function infection(dx, x, u, t, p=nothing)
S, I = u[1](t), u[2](t)
(β,γ,δt) = p
ifrac = rate_to_proportion(β*I,δt)
inf = ifrac*S
dx[1] = S-inf
dx[2] = I+inf
end
function recovery(dx, x, u, t, p=nothing)
S, I = u[1](t), u[2](t)
(β,γ,δt) = p
rfrac = rate_to_proportion(γ,δt)
rec = rfrac*I
dx[1] = S
dx[2] = I-rec
end
ofunc(x, u, t) = x
@def_discrete_system mutable struct InfectionSystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
β::Float64 = 0.5/1000.0
γ::Float64 = 0.25
δt::Float64 = 0.1
righthandside::RH = (dx, x, u, t, p=(β,γ,δt)) -> infection(dx, x, u, t, p)
state::Vector{Float64} = [990.0,10.0]
readout::RO = ofunc
input::IP = Inport(2)
output::OP = Outport(2)
end
@def_discrete_system mutable struct RecoverySystem{RH, RO, IP, OP} <: AbstractDiscreteSystem
β::Float64 = 0.5/1000.0
γ::Float64 = 0.25
δt::Float64 = 0.1
righthandside::RH = (dx, x, u, t, p=(β,γ,δt)) -> recovery(dx, x, u, t, p)
state::Vector{Float64} = [990.0,10.0]
readout::RO = ofunc
input::IP = Inport(2)
output::OP = Outport(2)
end
@defmodel model begin
@nodes begin
infsys = InfectionSystem()
recsys = RecoverySystem()
writer = Writer(input=Inport(2))
end
@branches begin
infsys[1:2] => recsys[1:2]
recsys[1:2] => infsys[1:2]
infsys[1:2] => writer[1:2]
end
end
simulate!(model, 0.0, 0.1, 40)
t, x = read(getnode(model, :writer).component)
plot(t, x[: ,1])
plot!(t, x[: ,2])
If you look at the output, there are a couple of differences between the output:
julia> x[1:7,1:2]
7×2 Matrix{Float64}:
990.0 10.0
989.505 10.4949
989.517 10.2358
989.517 10.2358
988.999 10.742
989.024 10.4768
989.024 10.4768
A similar thing happens with plugging recsys
into the writer:
julia> x[1:7,1:2]
7×2 Matrix{Float64}:
990.0 10.0
990.0 9.7531
990.0 9.7531
989.505 10.2358
989.517 9.98304
989.517 9.98304
988.999 10.4768
Here's the expected output (assuming that the models are composed into a single step):
990.0 10.0
989.505 10.248
988.998 10.5018
988.479 10.7617
987.947 11.0278
987.403 11.3001
986.845 11.5788
However, what I actually want is to have a system where for each time step, infsys is run, followed by recsys, followed by the writer taking the final output of recsys - only then should t be incremented.
I'm trying to use Causal.jl to re-implement one of my examples using DifferentialEquations here, modularizing it into two systems. It's a simple epidemiological model, with infection followed by recovery. According to the docs, this should be a static system, not a dynamic one, but although Causal.jl breaks the algebraic loop, it results in an error. Code below, noting that I have deliberately kept the code unnecessarily verbose. I'm also not sure how to feed the initial state (which should be [990.0,10.0]) as an input in a static system - do I need to use a memory component?
The output is: