Use Bareiss to select the basis vectors

Author: YingboMa

src/structural_transformation/partial_state_selection.jl

@@ -140,13 +140,14 @@ function partial_state_selection_graph!(structure::SystemStructure, var_eq_match
     var_eq_matching
 end
 
-function dummy_derivative_graph!(state::TransformationState)
+function dummy_derivative_graph!(state::TransformationState, jac=nothing)
     var_eq_matching = complete(pantelides!(state))
     complete!(state.structure)
-    dummy_derivative_graph!(state.structure, var_eq_matching)
+    # TODO: remove state when done
+    dummy_derivative_graph!(state.structure, var_eq_matching, jac, state)
 end
 
-function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching)
+function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, jac, state)
     @unpack eq_to_diff, var_to_diff, graph = structure
     diff_to_eq = invview(eq_to_diff)
     diff_to_var = invview(var_to_diff)
@@ -182,6 +183,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching)
     var_sccs = find_var_sccs(graph, var_eq_matching)
     eqcolor = falses(neqs)
     dummy_derivatives = Int[]
+    col_order = Int[]
     for vars in var_sccs
         eqs = [var_eq_matching[var] for var in vars if var_eq_matching[var] !== unassigned]
         isempty(eqs) && continue
@@ -195,18 +197,35 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching)
             iszero(nrows) && break
             eqs_set = BitSet(eqs)
 
-            structural_rank = 0
-            for var in vars
-                pathfound = construct_augmenting_path!(rank_matching, invgraph, var, eq->eq in eqs_set, eqcolor)
-                pathfound || continue
-                push!(dummy_derivatives, var)
-                structural_rank += 1
-                structural_rank == nrows && break
+            # TODO: making the algorithm more robust
+            # 1. If the Jacobian is a integer matrix, use Bareiss to check
+            # linear independence. (done)
+            #
+            # 2. If the Jacobian is a single row, generate pivots. (Dynamic
+            # state selection.)
+            #
+            # 3. If the Jacobian is a polynomial matrix, use Gröbner basis (?)
+            if jac !== nothing && (_J = jac(eqs, vars); all(x->unwrap(x) isa Integer, _J))
+                J = Int.(unwrap.(_J))
+                N = ModelingToolkit.nullspace(J; col_order) # modifies col_order
+                rank = length(col_order)-size(N, 2)
+                for i in 1:rank
+                    push!(dummy_derivatives, vars[col_order[i]])
+                end
+            else
+                rank = 0
+                for var in vars
+                    pathfound = construct_augmenting_path!(rank_matching, invgraph, var, eq->eq in eqs_set, eqcolor)
+                    pathfound || continue
+                    push!(dummy_derivatives, var)
+                    rank += 1
+                    rank == nrows && break
+                end
+                fill!(rank_matching, unassigned)
             end
-            if structural_rank != nrows
+            if rank != nrows
                 @warn "The DAE system is structurally singular!"
             end
-            fill!(rank_matching, unassigned)
 
             # prepare the next iteration
             eqs = map(eq->diff_to_eq[eq], eqs)

src/structural_transformation/symbolics_tearing.jl

@@ -291,6 +291,10 @@ the system is balanced.
 """
 function dummy_derivative(sys)
     state = TearingState(sys)
-    dds = dummy_derivative_graph!(state)
+    function jac(eqs, vars)
+        symeqs = EquationsView(state)[eqs]
+        Symbolics.jacobian((x->x.rhs).(symeqs), state.fullvars[vars])
+    end
+    dds = dummy_derivative_graph!(state, jac)
     EquationsView(state), state.fullvars[dds]
 end