WIP: robust dummy_derivatives and pantelides

Author: YingboMa

src/structural_transformation/pantelides.jl

@@ -83,30 +83,39 @@ function pantelides!(state::TransformationState, ag::Union{AliasGraph, Nothing}
     ecolor = falses(neqs)
     var_eq_matching = Matching(nvars)
     neqs′ = neqs
-    nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)), 1:neqs′)
+    nnonemptyeqs = count(eq -> !isempty(𝑠neighbors(graph, eq)) && eq_to_diff[eq] === nothing,
+                         1:neqs′)
 
     # Allow matching for the highest differentiated variable that
     # currently appears in an equation (or implicit equation in a side ag)
     varwhitelist = falses(nvars)
     for var in 1:nvars
-        if var_to_diff[var] === nothing
+        if var_to_diff[var] === nothing && !varwhitelist[var]
             while isempty(𝑑neighbors(graph, var)) && (ag === nothing || !haskey(ag, var))
                 var′ = invview(var_to_diff)[var]
                 var′ === nothing && break
                 var = var′
             end
             if !isempty(𝑑neighbors(graph, var)) || (ag !== nothing && haskey(ag, var))
-                varwhitelist[var] = true
+                if haskey(ag, var)
+                    # TODO: remove lower diff vars from whitelist
+                    c, a = ag[var]
+                    iszero(c) || (varwhitelist[a] = true)
+                else
+                    varwhitelist[var] = true
+                end
             end
         end
     end
+    @show varwhitelist
 
     if nnonemptyeqs > count(varwhitelist)
         throw(InvalidSystemException("System is structurally singular"))
     end
 
     for k in 1:neqs′
         eq′ = k
+        eq_to_diff[eq′] === nothing || continue
         isempty(𝑠neighbors(graph, eq′)) && continue
         pathfound = false
         # In practice, `maxiters=8000` should never be reached, otherwise, the

src/structural_transformation/partial_state_selection.jl

@@ -157,11 +157,12 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
     var_eq_matching
 end
 
-function dummy_derivative_graph!(state::TransformationState, jac = nothing; kwargs...)
+function dummy_derivative_graph!(state::TransformationState, jac = nothing, ag = nothing;
+                                 kwargs...)
     state.structure.solvable_graph === nothing && find_solvables!(state; kwargs...)
-    var_eq_matching = complete(pantelides!(state))
+    var_eq_matching = complete(pantelides!(state, ag === nothing ? nothing : first(ag)))
     complete!(state.structure)
-    dummy_derivative_graph!(state.structure, var_eq_matching, jac)
+    dummy_derivative_graph!(state.structure, var_eq_matching, jac, ag)
 end
 
 function compute_diff_level(diff_to_x)
@@ -181,8 +182,10 @@ function compute_diff_level(diff_to_x)
     return xlevel, maxlevel
 end
 
-function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac)
+function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac,
+                                 (ag, diff_va))
     @unpack eq_to_diff, var_to_diff, graph = structure
+    display(structure)
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
     invgraph = invview(graph)
@@ -194,6 +197,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
     dummy_derivatives = Int[]
     col_order = Int[]
     nvars = ndsts(graph)
+    @info "" var_eq_matching
     for vars in var_sccs
         eqs = [var_eq_matching[var] for var in vars if var_eq_matching[var] !== unassigned]
         isempty(eqs) && continue
@@ -245,6 +249,16 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
             vars = [diff_to_var[var] for var in vars if diff_to_var[var] !== nothing]
         end
     end
+    n_dummys = length(dummy_derivatives)
+    needed = count(x -> x isa Int, diff_to_eq) - n_dummys
+    n = 0
+    for v in diff_va
+        c, a = ag[v]
+        n += 1
+        push!(dummy_derivatives, iszero(c) ? v : a)
+        needed == n && break
+        continue
+    end
 
     dummy_derivatives_set = BitSet(dummy_derivatives)
     # We can eliminate variables that are not a selected state (differential
@@ -254,6 +268,9 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
         dummy_derivatives_set = dummy_derivatives_set
 
         v -> begin
+            if ag !== nothing
+                haskey(ag, v) && return false
+            end
             dv = var_to_diff[v]
             dv === nothing || dv in dummy_derivatives_set
         end
@@ -269,7 +286,8 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
                                        Union{Unassigned, SelectedState};
                                        varfilter = can_eliminate)
     for v in eachindex(var_eq_matching)
-        can_eliminate(v) && continue
+        dv = var_to_diff[v]
+        (dv === nothing || dv in dummy_derivatives_set) && continue
         var_eq_matching[v] = SelectedState()
     end