Skip to content

Commit d1189ed

Browse files
YingboMashashi
andcommitted
Add bigraph conversion and pretty printing for debugging
Co-authored-by: "Shashi Gowda" <[email protected]> Co-authored-by: "Yingbo Ma" <[email protected]>
1 parent e5621bc commit d1189ed

File tree

2 files changed

+51
-0
lines changed

2 files changed

+51
-0
lines changed

src/systems/diffeqs/index_reduction.jl

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,3 +12,52 @@ function sys2bigraph(sys)
1212
push!(data, es)
1313
end
1414
end
15+
16+
function get_vnodes(sys)
17+
diffnodes = []
18+
diffedges = Tuple{Int, Int}[]
19+
for (i, eq) in enumerate(sys.eqs)
20+
if !(eq.lhs isa Constant)
21+
# Make sure that the LHS is a first order derivative of a var.
22+
@assert eq.lhs.op isa Differential
23+
@assert !(eq.lhs.args[1] isa Differential) # first order
24+
25+
push!(diffnodes, eq.lhs)
26+
# For efficiency we note down the diff edges here
27+
push!(diffedges, (i, length(diffnodes)))
28+
end
29+
end
30+
31+
diffvars = (first var_from_nested_derivative).(diffnodes)
32+
algvars = setdiff(states(sys), diffvars)
33+
return diffnodes, diffedges, algvars
34+
end
35+
36+
function sys2bigraph2(sys)
37+
diffvars, edges, algvars = get_vnodes(sys)
38+
varnumber_offset = length(diffvars)
39+
40+
for (i, eq) in enumerate(sys.eqs)
41+
# T or D(x):
42+
# We assume no derivatives appear on the RHS at this point
43+
vs = vars(eq.rhs)
44+
for v in vs
45+
for (j, target_v) in enumerate(algvars)
46+
if v == target_v
47+
push!(edges, (i, j+varnumber_offset))
48+
end
49+
end
50+
end
51+
end
52+
vcat(diffvars, algvars), edges
53+
end
54+
55+
print_bigraph(sys, vars, edges) = print_bigraph(stdout, sys, vars, edges)
56+
function print_bigraph(io::IO, sys, vars, edges)
57+
println(io, "Equations:")
58+
foreach(x->println(io, x), [i => sys.eqs[i] for i in 1:length(sys.eqs)])
59+
println(io)
60+
for (i, j) in edges
61+
println(io, "Eq $i has $(vars[j])")
62+
end
63+
end

test/index_reduction.jl

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -30,4 +30,6 @@ eqs = [D(x) ~ w,
3030
0 ~ x^2 + y^2 - L^2]
3131
pendulum = ODESystem(eqs, t, [x, y, w, z, T], [L, g], name=:pendulum)
3232

33+
# V-nodes D(x), D(y), D(w), D(z), T
34+
# E-nodes
3335
sys2bigraph(pendulum)

0 commit comments

Comments
 (0)