|
| 1 | +# V-nodes `[x_1, x_2, x_3, ..., dx_1, dx_2, ..., y_1, y_2, ...]` where `x`s are |
| 2 | +# differential variables and `y`s are algebraic variables. |
| 3 | +function get_vnodes(sys) |
| 4 | + dxvars = Operation[] |
| 5 | + edges = map(_->Int[], 1:length(sys.eqs)) |
| 6 | + for (i, eq) in enumerate(sys.eqs) |
| 7 | + if !(eq.lhs isa Constant) |
| 8 | + # Make sure that the LHS is a first order derivative of a var. |
| 9 | + @assert eq.lhs.op isa Differential |
| 10 | + @assert !(eq.lhs.args[1] isa Differential) # first order |
| 11 | + |
| 12 | + push!(dxvars, eq.lhs) |
| 13 | + # For efficiency we note down the diff edges here |
| 14 | + push!(edges[i], length(dxvars)) |
| 15 | + end |
| 16 | + end |
| 17 | + |
| 18 | + xvars = (first ∘ var_from_nested_derivative).(dxvars) |
| 19 | + algvars = setdiff(states(sys), xvars) |
| 20 | + return xvars, dxvars, edges, algvars |
| 21 | +end |
| 22 | + |
| 23 | +function sys2bigraph(sys) |
| 24 | + xvars, dxvars, edges, algvars = get_vnodes(sys) |
| 25 | + xvar_offset = length(xvars) |
| 26 | + algvar_offset = 2xvar_offset |
| 27 | + for edge in edges |
| 28 | + isempty(edge) || (edge .+= xvar_offset) |
| 29 | + end |
| 30 | + |
| 31 | + for (i, eq) in enumerate(sys.eqs) |
| 32 | + # T or D(x): |
| 33 | + # We assume no derivatives appear on the RHS at this point |
| 34 | + vs = vars(eq.rhs) |
| 35 | + for v in vs |
| 36 | + for (j, target_v) in enumerate(xvars) |
| 37 | + if v == target_v |
| 38 | + push!(edges[i], j) |
| 39 | + end |
| 40 | + end |
| 41 | + for (j, target_v) in enumerate(algvars) |
| 42 | + if v == target_v |
| 43 | + push!(edges[i], j+algvar_offset) |
| 44 | + end |
| 45 | + end |
| 46 | + end |
| 47 | + end |
| 48 | + |
| 49 | + fullvars = [xvars; dxvars; algvars] # full list of variables |
| 50 | + vars_asso = [(1:xvar_offset) .+ xvar_offset; zeros(Int, length(fullvars) - xvar_offset)] # variable association list |
| 51 | + return edges, fullvars, vars_asso |
| 52 | +end |
| 53 | + |
| 54 | +print_bigraph(sys, vars, edges) = print_bigraph(stdout, sys, vars, edges) |
| 55 | +function print_bigraph(io::IO, sys, vars, edges) |
| 56 | + println(io, "Equations:") |
| 57 | + foreach(x->println(io, x), [i => sys.eqs[i] for i in 1:length(sys.eqs)]) |
| 58 | + for (i, edge) in enumerate(edges) |
| 59 | + println(io, "\nEq $i has:") |
| 60 | + print(io, '[') |
| 61 | + for e in edge |
| 62 | + print(io, "$(vars[e]), ") |
| 63 | + end |
| 64 | + print(io, ']') |
| 65 | + end |
| 66 | + return nothing |
| 67 | +end |
| 68 | + |
| 69 | +function match_equation!(edges, i, assign, active, vcolor=falses(length(active)), ecolor=falses(length(edges))) |
| 70 | + # `edge[active]` are active edges |
| 71 | + # i: equations |
| 72 | + # j: variables |
| 73 | + # assign: assign[j] == i means (i-j) is assigned |
| 74 | + # |
| 75 | + # color the equation |
| 76 | + ecolor[i] = true |
| 77 | + # if a V-node j exists s.t. edge (i-j) exists and assign[j] == 0 |
| 78 | + for j in edges[i] |
| 79 | + if active[j] && assign[j] == 0 |
| 80 | + assign[j] = i |
| 81 | + return true |
| 82 | + end |
| 83 | + end |
| 84 | + # for every j such that edge (i-j) exists and j is uncolored |
| 85 | + for j in edges[i] |
| 86 | + (active[j] && !vcolor[j]) || continue |
| 87 | + # color the variable |
| 88 | + vcolor[j] = true |
| 89 | + if match_equation!(edges, assign[j], assign, active, vcolor, ecolor) |
| 90 | + assign[j] = i |
| 91 | + return true |
| 92 | + end |
| 93 | + end |
| 94 | + return false |
| 95 | +end |
| 96 | + |
| 97 | +function matching(edges, nvars, active=trues(nvars)) |
| 98 | + assign = zeros(Int, nvars) |
| 99 | + for i in 1:length(edges) |
| 100 | + match_equation!(edges, i, assign, active) |
| 101 | + end |
| 102 | + return assign |
| 103 | +end |
| 104 | + |
| 105 | +function pantelides(sys::ODESystem; kwargs...) |
| 106 | + edges, fullvars, vars_asso = sys2bigraph(sys) |
| 107 | + return pantelides!(edges, fullvars, vars_asso; kwargs...) |
| 108 | +end |
| 109 | + |
| 110 | +function pantelides!(edges, vars, vars_asso; maxiter = 8000) |
| 111 | + neqs = length(edges) |
| 112 | + nvars = length(vars) |
| 113 | + assign = zeros(Int, nvars) |
| 114 | + eqs_asso = fill(0, neqs) |
| 115 | + neqs′ = neqs |
| 116 | + for k in 1:neqs′ |
| 117 | + i = k |
| 118 | + pathfound = false |
| 119 | + # In practice, `maxiter=8000` should never be reached, otherwise, the |
| 120 | + # index would be on the order of thousands. |
| 121 | + for _ in 1:maxiter |
| 122 | + # run matching on (dx, y) variables |
| 123 | + active = vars_asso .== 0 |
| 124 | + vcolor = falses(nvars) |
| 125 | + ecolor = falses(neqs) |
| 126 | + pathfound = match_equation!(edges, i, assign, active, vcolor, ecolor) |
| 127 | + pathfound && break # terminating condition |
| 128 | + # for every colored V-node j |
| 129 | + for j in eachindex(vcolor); vcolor[j] || continue |
| 130 | + # introduce a new variable |
| 131 | + nvars += 1 |
| 132 | + push!(vars_asso, 0) |
| 133 | + vars_asso[j] = nvars |
| 134 | + push!(assign, 0) |
| 135 | + end |
| 136 | + |
| 137 | + # for every colored E-node l |
| 138 | + for l in eachindex(ecolor); ecolor[l] || continue |
| 139 | + neqs += 1 |
| 140 | + # create new E-node |
| 141 | + push!(edges, copy(edges[l])) |
| 142 | + # create edges from E-node `neqs` to all V-nodes `j` and |
| 143 | + # `vars_asso[j]` s.t. edge `(l-j)` exists |
| 144 | + for j in edges[l] |
| 145 | + if !(vars_asso[j] in edges[neqs]) |
| 146 | + push!(edges[neqs], vars_asso[j]) |
| 147 | + end |
| 148 | + end |
| 149 | + push!(eqs_asso, 0) |
| 150 | + eqs_asso[l] = neqs |
| 151 | + end |
| 152 | + |
| 153 | + # for every colored V-node j |
| 154 | + for j in eachindex(vcolor); vcolor[j] || continue |
| 155 | + assign[vars_asso[j]] = eqs_asso[assign[j]] |
| 156 | + end |
| 157 | + i = eqs_asso[i] |
| 158 | + end # for _ in 1:maxiter |
| 159 | + pathfound || error("maxiter=$maxiter reached! File a bug report if your system has a reasonable index (<100), and you are using the default `maxiter`. Try to increase the maxiter by `pantelides(sys::ODESystem; maxiter=1_000_000)` if your system has an incredibly high index and it is truly extremely large.") |
| 160 | + end # for k in 1:neqs′ |
| 161 | + return edges, assign, vars_asso, eqs_asso |
| 162 | +end |
0 commit comments