Add Pantelides algorithm

Author: YingboMa

src/systems/diffeqs/index_reduction.jl

@@ -66,8 +66,7 @@ function print_bigraph(io::IO, sys, vars, edges)
     return nothing
 end
 
-
-function matching_equation!(edges, i, assign, active, vcolor=falses(length(active)), ecolor=falses(length(edges)))
+function match_equation!(edges, i, assign, active, vcolor=falses(length(active)), ecolor=falses(length(edges)))
     # `edge[active]` are active edges
     # i: equations
     # j: variables
@@ -88,7 +87,7 @@ function matching_equation!(edges, i, assign, active, vcolor=falses(length(activ
         # color the variable
         vcolor[j] = true
         if match_equation!(edges, assign[j], assign, active, vcolor, ecolor)
-            assign[v] = i
+            assign[j] = i
             return true
         end
     end
@@ -98,7 +97,66 @@ end
 function matching(edges, nvars, active=trues(nvars))
     assign = zeros(Int, nvars)
     for i in 1:length(edges)
-        matching_equation!(edges, i, assign, active)
+        match_equation!(edges, i, assign, active)
     end
     return assign
 end
+
+function pantelides(sys::ODESystem)
+    edges, fullvars, vars_asso = sys2bigraph(sys)
+    return pantelides!(edges, fullvars, vars_asso)
+end
+
+function pantelides!(edges, vars, vars_asso)
+    neqs = length(edges)
+    nvars = length(vars)
+    assign = zeros(Int, nvars)
+    eqs_asso = fill(0, neqs)
+    neqs′ = neqs
+    for k in 1:neqs′
+        i = k
+        pathfound = false
+        iter = 0
+        while !pathfound && iter < 3
+            iter += 1
+            # run matching on (dx, y) variables
+            active = vars_asso .== 0
+            vcolor = falses(nvars)
+            ecolor = falses(neqs)
+            pathfound = match_equation!(edges, i, assign, active, vcolor, ecolor)
+            if !pathfound
+                # for every colored V-node j
+                for j in eachindex(vcolor); vcolor[j] || continue
+                    # introduce a new variable
+                    nvars += 1
+                    push!(vars_asso, 0)
+                    vars_asso[j] = nvars
+                    push!(assign, 0)
+                end
+
+                # for every colored E-node l
+                for l in eachindex(ecolor); ecolor[l] || continue
+                    neqs += 1
+                    # create new E-node
+                    push!(edges, copy(edges[l]))
+                    # create edges from E-node `neqs` to all V-nodes `j` and
+                    # `vars_asso[j]` s.t. edge `(l-j)` exists
+                    for j in edges[l]
+                        if !(vars_asso[j] in edges[neqs])
+                            push!(edges[neqs], vars_asso[j])
+                        end
+                    end
+                    push!(eqs_asso, 0)
+                    eqs_asso[l] = neqs
+                end
+
+                # for every colored V-node j
+                for j in eachindex(vcolor); vcolor[j] || continue
+                    assign[vars_asso[j]] = eqs_asso[assign[j]]
+                end
+                i = eqs_asso[i]
+            end # if !pathfound
+        end # while !pathfound
+    end # for k in 1:neqs′
+    return assign, vars_asso, eqs_asso
+end