Help with running simple example

[sdwfrost]

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?

using Causal

@inline function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end

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

function recovery(t,u,p)
    (S,I) = u
    (β,γ,δt) = p
    rfrac = rate_to_proportion(γ,δt)
    rec = rfrac*I
    return [S,I-rec]
end

@def_static_system struct InfectionSystem{IP, OP, RO} <: AbstractStaticSystem
      β::Float64 = 0.5/1000
      γ::Float64 = 0.25
      δt::Float64 = 0.1
      p = (β,γ,δt)
      input::IP = Inport(2)
      output::OP = Outport(2)
      readout::RO = (t,u) -> infection(t,u,p)
end

@def_static_system struct RecoverySystem{IP, OP, RO} <: AbstractStaticSystem
      β::Float64 = 0.5/1000
      γ::Float64 = 0.25
      δt::Float64 = 0.1
      p = (β,γ,δt)
      input::IP = Inport(2)
      output::OP = Outport(2)
      readout::RO = (t,u) -> recovery(t,u,p)
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) # run for single step

The output is:

[ Info: 2021-04-16T14:42:40.211 Started simulation...
[ Info: 2021-04-16T14:42:40.211 Inspecting model...
┌ Info:     The model has algrebraic loops:[[1, 2]]
└       Trying to break these loops...
[ Info:     Loop [1, 2] is broken
[ Info: 2021-04-16T14:42:41.039 Done.
[ Info: 2021-04-16T14:42:41.039 Initializing the model...
[ Info: 2021-04-16T14:42:41.081 Done...
[ Info: 2021-04-16T14:42:44.434 Running the simulation...
┌ Error: Failed for Causal.LoopBreaker{Outport{Outpin{Float64}}, Causal.var"#926#929"{Causal.var"#951#952"{Vector{Function}, Int64}}, Inpin{Float64}, Outpin{Bool}, Nothing, Base.UUID}(nothing, Outport(numpins:2, eltype:Outpin{Float64}), Causal.var"#926#929"{Causal.var"#951#952"{Vector{Function}, Int64}}(Core.Box(2), Causal.var"#951#952"{Vector{Function}, Int64}(Function[Causal.var"#func#938"{BitVector, InfectionSystem{Inport{Inpin{Float64}}, Outport{Outpin{Float64}}, var"#8#13"{Tuple{Float64, Float64, Float64}}, Inpin{Float64}, Outpin{Bool}, Nothing, Base.UUID}}(Bool[1, 1], InfectionSystem{Inport{Inpin{Float64}}, Outport{Outpin{Float64}}, var"#8#13"{Tuple{Float64, Float64, Float64}}, Inpin{Float64}, Outpin{Bool}, Nothing, Base.UUID}(0.0005, 0.25, 0.1, (0.0005, 0.25, 0.1), Inport(numpins:2, eltype:Inpin{Float64}), Outport(numpins:2, eltype:Outpin{Float64}), var"#8#13"{Tuple{Float64, Float64, Float64}}((0.0005, 0.25, 0.1)), Inpin(eltype:Float64, isbound:true), Outpin(eltype:Bool, isbound:true), nothing, Symbol(""), UUID("fa11882b-390d-4141-a53d-3059090e747b"))), Causal.var"#func#937"{RecoverySystem{Inport{Inpin{Float64}}, Outport{Outpin{Float64}}, var"#18#23"{Tuple{Float64, Float64, Float64}}, Inpin{Float64}, Outpin{Bool}, Nothing, Base.UUID}}(RecoverySystem{Inport{Inpin{Float64}}, Outport{Outpin{Float64}}, var"#18#23"{Tuple{Float64, Float64, Float64}}, Inpin{Float64}, Outpin{Bool}, Nothing, Base.UUID}(0.0005, 0.25, 0.1, (0.0005, 0.25, 0.1), Inport(numpins:2, eltype:Inpin{Float64}), Outport(numpins:2, eltype:Outpin{Float64}), var"#18#23"{Tuple{Float64, Float64, Float64}}((0.0005, 0.25, 0.1)), Inpin(eltype:Float64, isbound:true), Outpin(eltype:Bool, isbound:true), nothing, Symbol(""), UUID("d78e5f28-2174-45da-af02-c054edb68f44")))], 2)), Inpin(eltype:Float64, isbound:true), Outpin(eltype:Bool, isbound:true), nothing, Symbol(""), UUID("2a9ba428-100f-4702-96f1-687044c07531"))
└ @ Causal ~/.julia/packages/Causal/vsleZ/src/models/taskmanager.jl:39

[zekeriyasari]

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;

[sdwfrost]

Hi @zekeriyasari

Thanks for your detailed response. Perhaps I can clear up some vagueness in my question!

1. I understand that the static system doesn't have any internal dynamics - what I am trying to do is run a single 'step' of the difference equation, such that there is not an internal state except that which is in input and output.
2. My full example has three states, but my test case in Causal only considers two: S and I. For some parameter values, this is numerically equivalent to the example I give. I used this smaller system in demonstrating AlgebraicDynamics.jl, an acausal framework: https://github.com/epirecipes/sir-julia/blob/master/markdown/ode_algebraicdynamics/ode_algebraicdynamics.md.
3. You're right in that I want to model the two processes separately. For the 'full' system, I would introduce some coupling system like you have in some of your examples to couple infection (S,I) with recovery (I,R).

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

[zekeriyasari]

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.Interpolantand 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])

[sdwfrost]

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

1. The number of S can go up in a timestep (it should only ever go down)
2. The rows are grouped in threes - one unique row followed by duplicates.

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.