How to extract the dependency (also called causal) graph of a ModelingToolkit.jl ODE model?

I am trying to create the dependency (sometimes called causal) graph of a complex ODE System I have created symbolically with MTK. To do this I am using the dependency graph functionality of MTK [found here], however I cannot interpret it correctly to create what I need. I will present a MWE and what I want to create and perhaps someone can help me interpret and process the asgraph output to get what I need?

Here is the simple model:

using ModelingToolkit
using ModelingToolkit: t_nounits as t

@variables x(t) = 1.0 y(t) = 1.0 z(t)

eqs = [
    Differential(t)(x) ~ -2x + z,
    Differential(t)(y) ~ -2y,
    z ~ x*y,
]

sys = ODESystem(eqs, t; name = :test)

ssys = structural_simplify(sys)

go = asgraph(sys; variables = [x, y, z])
gs = asgraph(ssys; variables = [x, y, z])

the output of go (which is what I believe I must use) is:

julia> go
BipartiteGraph with (3, 3) (đť‘ ,đť‘‘)-vertices
  #  src       dst
 1     [1, 3]    [1, 3]
 2     [2]       [2, 3]
 3     [1, 2]    [1]

What I want to create is:

(I know how to plot graphs if I create the graph above).

Ultimately, I want to process the graph above to create the true and full causal graph of the system which is:

However just getting the first picture is more than enough for me!

1 Like

Your first picture is just the src from go with two changes where you change indices for a variable to its derivative if it’s a differential variable. Did you try that?

Thanks that helps a bit. I realized that with this approach I need to keep track of also the order of equations that the sys stores for later processing, which for my actual use case is rather tedious.

Looking later down in the same documentation page, I actually found something that almost automatically creates what I need:

variables = [x, y, z]
varvardep = varvar_dependencies(asgraph(sys; variables), variable_dependencies(sys; variables))

using CairoMakie, GraphMakie
graphplot(varvardep; nlabels = string.([x, y, z]))

I think from the above result, I only need to process all self links to be new links pointing to new variables, the derivatives.

Is it correct to assume that all self links are from variables to their derivatives? Can you imagine any other scenario where a self-link would exist? (I can’t from the top of my head). I understand that given that I know a-priory which variables of my system are differential I can manually make this distinction on the variables vector above, but I am hoping to get a programmatic way to distinguish this. Or is there perhaps a function offered by MTK is_differential(sys, variable) that spits true if that’s the case?

ah actually my approach above is simply incorrect.

What I have to do is first identify which of the variables are time-derivative’d. Then, for each, I have to make all incoming links point to a new variable (the derivative). So a programmatic way to check if a variable is derivative’d (i.e., a state variable) in the sys would be great!

I’m using something like this:

function get_diff_variables(sys)
    only.(get_variables.(getproperty.(diff_equations(sys), :lhs)))
end

where get_variables is from Symbolics and diff_equations from ModelingToolkit.