# Causal modeling in MTK?

**URL:** <https://discourse.julialang.org/t/causal-modeling-in-mtk/117161>\
**Category:** Modelling & Simulations\
**Created:** [July 17, 2024, 8:07pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161 "2024-07-17T20:07:02Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [July 17, 2024, 8:07pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/1 "2024-07-17T20:07:02Z")

</div>

After seeing @ChrisRackauckas slides from JuliaCon 2024, I found out that (it seems) ModelingToolkit supports causal models, through passing input and outputs around. Is this really possible? I ask because Causal.jl struggles with a model I’m trying to implement - see [this issue](https://github.com/zekeriyasari/Causal.jl/issues/71). Any tips on porting this Causal.jl model to MTK would be appreciated! I’m not sure how to define the component models in MTK.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 17, 2024, 10:03pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/2 "2024-07-17T22:03:19Z")

</div>

I whipped this together really quick and haven’t run it yet, but yeah you can define two systems, one with an input and and one with an output, and connect them. It would look something like this:

```julia
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

@mtkmodel InfectionSystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        S
        u [output = true]
    end
    @equations begin
        D(S) ~ rate_to_proportion(β*I,δt)*S
        u ~ rate_to_proportion(β*I,δt)*S
    end
end

@mtkmodel RecoverySystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        R
        u [input = true]
    end
    @equations begin
        D(R) ~ u - rate_to_proportion(γ,δt)*I
    end
end

@mtkmodel SIR begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @components begin
        infection = InfectionSystem()
        recovery = RecoverySystem()
    end
    @equations begin
        connect(infection.u, recovery.u)
    end
end

```

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [July 17, 2024, 11:16pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/3 "2024-07-17T23:16:14Z")

</div>

Thanks @ChrisRackauckas! The line `u [output=true]` throws a parse error (as it expects no gap between `u` and `[`.

What’s going on under the hood? Is it a true causal system where `InfectionSystem` would be run, then `RecoverySystem`?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 18, 2024, 2:28am UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/4 "2024-07-18T02:28:15Z")

</div>

> [@sdwfrost](#):
>
> Thanks @ChrisRackauckas! The line `u [output=true]` throws a parse error (as it expects no gap between `u` and `[`.

```julia
    @variables begin
        R
        u, [input = true]
    end

```

the comma was missing.

> [@sdwfrost](#):
>
> What’s going on under the hood? Is it a true causal system where `InfectionSystem` would be run, then `RecoverySystem`?

That’s orthogonal. That’s separate compilation.

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [July 19, 2024, 5:31pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/5 "2024-07-19T17:31:09Z")

</div>

The following code throws an error when I try to instantiate the model that I can’t easily debug:

```Julia
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D

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

@mtkmodel InfectionSystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        S(t)
        I(t)
        u(t),[output=true]
    end
    @equations begin
        D(S) ~ rate_to_proportion(β*I,δt)*S
        u ~ rate_to_proportion(β*I,δt)*S
    end
end

@mtkmodel RecoverySystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        I(t)
        R(t)
        u(t),[input=true]
    end
    @equations begin
        D(R) ~ u - rate_to_proportion(γ,δt)*I
    end
end

@mtkmodel SIR begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @components begin
        infection = InfectionSystem()
        recovery = RecoverySystem()
    end
    @equations begin
        connect(infection.u, recovery.u)
    end
end

@mtkbuild sir=SIR()

```

Output:

```Julia
ERROR: None Sym BasicSymbolic doesn't have a name
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] nameof(s::SymbolicUtils.BasicSymbolic{Real})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/qyMYa/src/types.jl:249
  [3] nameof(n::Num)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/num.jl:88
  [4] unique(f::typeof(nameof), C::Tuple{Num, Num}; seen::Nothing)
    @ Base ./set.jl:261
  [5] unique
    @ ./set.jl:247 [inlined]
  [6] connect(::Num, ::Num)
    @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/equations.jl:18
  [7] __SIR__ (; name::Symbol, β::Nothing, γ::Nothing, δt::Nothing)
    @ Main ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/model_parsing.jl:620
  [8] __SIR__
    @ ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/model_parsing.jl:132 [inlined]
  [9] #_#321
    @ ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/model_parsing.jl:25 [inlined]
 [10] top-level scope
    @ ~/.julia/packages/ModelingToolkit/cAhZr/src/systems/abstractsystem.jl:1445

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 20, 2024, 1:32am UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/6 "2024-07-20T01:32:46Z")

</div>

Using connectors is probably nicer. Here’s a version that builds:

```julia
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary.Blocks: RealInput, RealOutput

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

@mtkmodel InfectionSystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        S(t)
    end
    @components begin
        output = RealOutput()
    end
    @equations begin
        D(S) ~ rate_to_proportion(β*I,δt)*S
        output.u ~ rate_to_proportion(β*I,δt)*S
    end
end

@mtkmodel RecoverySystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        I(t)
        R(t)
    end
    @components begin
        input = RealInput()
    end
    @equations begin
        D(I) ~ rate_to_proportion(γ,δt)*I
        D(R) ~ input.u - rate_to_proportion(γ,δt)*I
    end
end

@mtkmodel SIR begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @components begin
        infection = InfectionSystem()
        recovery = RecoverySystem()
    end
    @equations begin
        connect(infection.output, recovery.input)
    end
end

@mtkbuild sir=SIR()

```

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [July 23, 2024, 2:03pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/7 "2024-07-23T14:03:58Z")

</div>

That doesn’t build for me; I had to add in `I` as a variable in `InfectionSystem`, but the following gives me an error.

```julia
ERROR: ExtraVariablesSystemException: The system is unbalanced. There are 6 highest order derivative variables and 5 equations.
More variables than equations, here are the potential extra variable(s):
 infection₊S(t)
 recovery₊I(t)
 recovery₊R(t)
 recovery₊input₊u(t)
Stacktrace:
  [1] error_reporting(state::TearingState{…}, bad_idxs::Vector{…}, n_highest_vars::Int64, iseqs::Bool, orig_inputs::Set{…})
    @ ModelingToolkit.StructuralTransformations C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\structural_transformation\utils.jl:50
  [2] check_consistency(state::TearingState{ODESystem}, orig_inputs::Set{Any})
    @ ModelingToolkit.StructuralTransformations C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\structural_transformation\utils.jl:85
  [3] _structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systemstructure.jl:691
  [4] structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systemstructure.jl:634
  [5] __structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkit C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systems.jl:85
  [6] __structural_simplify
    @ C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systems.jl:66 [inlined]
  [7] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, split::Bool, allow_symbolic::Bool, allow_parameter::Bool, conservative::Bool, fully_determined::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systems.jl:24
  [8] structural_simplify
    @ C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systems.jl:20 [inlined]
  [9] structural_simplify(sys::ODESystem)
    @ ModelingToolkit C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\systems.jl:20
 [10] top-level scope
    @ C:\Users\sdwfr\.julia\packages\ModelingToolkit\3AmNp\src\systems\abstractsystem.jl:1881
Some type information was truncated. Use `show(err)` to see complete types.

```

```julia
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkitStandardLibrary.Blocks: RealInput, RealOutput

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

@mtkmodel InfectionSystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        S(t)
        I(t)
    end
    @components begin
        output = RealOutput()
    end
    @equations begin
        D(S) ~ rate_to_proportion(β*I,δt)*S
        output.u ~ rate_to_proportion(β*I,δt)*S
    end
end

@mtkmodel RecoverySystem begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @variables begin
        I(t)
        R(t)
    end
    @components begin
        input = RealInput()
    end
    @equations begin
        D(I) ~ rate_to_proportion(γ,δt)*I
        D(R) ~ input.u - rate_to_proportion(γ,δt)*I
    end
end

@mtkmodel SIR begin
    @parameters begin
        β = 0.5/1000
        γ = 0.25
        δt = 0.1
    end
    @components begin
        infection = InfectionSystem()
        recovery = RecoverySystem()
    end
    @equations begin
        connect(infection.output, recovery.input)
    end
end

@mtkbuild sir=SIR()

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 23, 2024, 2:06pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/8 "2024-07-23T14:06:14Z")

</div>

You cannot have I in both or you would be missing a connection. This error is saying to use what I shared

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [July 23, 2024, 3:04pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/9 "2024-07-23T15:04:16Z")

</div>

…but then the InfectionSystem doesn’t compile (as it doesn’t know about `I`)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 23, 2024, 9:09pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/10 "2024-07-23T21:09:36Z")

</div>

The models are incomplete without connections (which is known because of the inputs/outputs), but the connected model is complete and it compiles fine (the code I shared)

---

<div class="post-metadata">

**Author:** ![sdwfrost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdwfrost/32/2831_2.png) [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)\
**Post date:** [July 23, 2024, 9:33pm UTC](https://discourse.julialang.org/t/causal-modeling-in-mtk/117161/11 "2024-07-23T21:33:09Z")

</div>

Sorry for not being clear - the code you shared throws this error when building the final system:

```julia
ERROR: UndefVarError: `I` not defined
Stacktrace:
 [1] __InfectionSystem__ (; name::Symbol, β::Nothing, γ::Nothing, δt::Nothing, S::Nothing)

```

My version of MTK:

```julia
(models) pkg> status ModelingToolkit
  [961ee093] ModelingToolkit v9.26.0

```