Tearing differentiated variables on priority

Author: YingboMa

src/structural_transformation/bipartite_tearing/modia_tearing.jl

@@ -9,34 +9,62 @@ function try_assign_eq!(ict::IncrementalCycleTracker, vj::Integer, eq::Integer)
     end
 end
 
-function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int},
-                        vs::Vector{Int})
+function try_assign_eq!(ict::IncrementalCycleTracker, vars, v_active, eq::Integer,
+                        condition::F = _ -> true) where {F}
     G = ict.graph
-    vActive = BitSet(vs)
+    for vj in vars
+        (vj in v_active && G.matching[vj] === unassigned && condition(vj)) || continue
+        try_assign_eq!(ict, vj, eq) && return true
+    end
+    return false
+end
 
+function tearEquations!(ict::IncrementalCycleTracker, Gsolvable, es::Vector{Int},
+                        v_active::BitSet, isder′::F) where {F}
+    check_der = isder′ !== nothing
+    if check_der
+        has_der = Ref(false)
+        isder = let has_der = has_der, isder′ = isder′
+            v -> begin
+                r = isder′(v)
+                has_der[] |= r
+                r
+            end
+        end
+    end
     for eq in es  # iterate only over equations that are not in eSolvedFixed
-        for vj in Gsolvable[eq]
-            if G.matching[vj] === unassigned && (vj in vActive)
-                r = try_assign_eq!(ict, vj, eq)
-                r && break
+        vs = Gsolvable[eq]
+        #=
+        if check_der
+            # if there're differentiated variables, then only consider them
+            try_assign_eq!(ict, vs, v_active, eq, isder)
+            @show has_der[]
+            if has_der[]
+                has_der[] = false
+                continue
             end
         end
+        =#
+        try_assign_eq!(ict, vs, v_active, eq)
     end
 
     return ict
 end
 
-function tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, eqs, vars)
+function tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, eqs, vars,
+                                 isder::F) where {F}
     ict = IncrementalCycleTracker(DiCMOBiGraph{true}(graph); dir = :in)
-    tearEquations!(ict, solvable_graph.fadjlist, eqs, vars)
+    tearEquations!(ict, solvable_graph.fadjlist, eqs, vars, isder)
     for var in vars
         var_eq_matching[var] = ict.graph.matching[var]
     end
     return nothing
 end
 
-function tear_graph_modia(structure::SystemStructure, ::Type{U} = Unassigned;
-                          varfilter = v -> true, eqfilter = eq -> true) where {U}
+function tear_graph_modia(structure::SystemStructure, isder::F = nothing,
+                          ::Type{U} = Unassigned;
+                          varfilter::F2 = v -> true,
+                          eqfilter::F3 = eq -> true) where {F, U, F2, F3}
     # It would be possible here to simply iterate over all variables and attempt to
     # use tearEquations! to produce a matching that greedily selects the minimal
     # number of torn variables. However, we can do this process faster if we first
@@ -52,14 +80,22 @@ function tear_graph_modia(structure::SystemStructure, ::Type{U} = Unassigned;
     var_eq_matching = complete(maximal_matching(graph, eqfilter, varfilter, U))
     var_sccs::Vector{Union{Vector{Int}, Int}} = find_var_sccs(graph, var_eq_matching)
 
+    ieqs = Int[]
+    filtered_vars = BitSet()
     for vars in var_sccs
-        filtered_vars = filter(varfilter, vars)
-        ieqs = Int[var_eq_matching[v]
-                   for v in filtered_vars if var_eq_matching[v] !== unassigned]
         for var in vars
+            if varfilter(var)
+                push!(filtered_vars, var)
+                if var_eq_matching[var] !== unassigned
+                    push!(ieqs, var_eq_matching[var])
+                end
+            end
             var_eq_matching[var] = unassigned
         end
-        tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars)
+        tear_graph_block_modia!(var_eq_matching, graph, solvable_graph, ieqs, filtered_vars,
+                                isder)
+        empty!(ieqs)
+        empty!(filtered_vars)
     end
 
     return var_eq_matching

src/structural_transformation/partial_state_selection.jl

@@ -84,7 +84,8 @@ function pss_graph_modia!(structure::SystemStructure, var_eq_matching, varlevel,
             end
             filter!(var -> ict.graph.matching[var] === unassigned, to_tear_vars)
             filter!(eq -> invview(ict.graph.matching)[eq] === unassigned, to_tear_eqs)
-            tearEquations!(ict, solvable_graph.fadjlist, to_tear_eqs, to_tear_vars)
+            tearEquations!(ict, solvable_graph.fadjlist, to_tear_eqs, BitSet(to_tear_vars),
+                           nothing)
             for var in to_tear_vars
                 var_eq_matching[var] = ict.graph.matching[var]
             end
@@ -177,7 +178,6 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
     invgraph = invview(graph)
 
     eqlevel, _ = compute_diff_level(diff_to_eq)
-    varlevel, _ = compute_diff_level(diff_to_var)
 
     var_sccs = find_var_sccs(graph, var_eq_matching)
     eqcolor = falses(nsrcs(graph))
@@ -191,7 +191,7 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
         iszero(maxlevel) && continue
 
         rank_matching = Matching(nvars)
-        for level in maxlevel:-1:1
+        for _ in maxlevel:-1:1
             eqs = filter(eq -> diff_to_eq[eq] !== nothing, eqs)
             nrows = length(eqs)
             iszero(nrows) && break
@@ -254,13 +254,10 @@ function dummy_derivative_graph!(structure::SystemStructure, var_eq_matching, ja
     isdiffed = let diff_to_var = diff_to_var, dummy_derivatives_set = dummy_derivatives_set
         v -> diff_to_var[v] !== nothing && !(v in dummy_derivatives_set)
     end
-    should_consider = let graph = graph, isdiffed = isdiffed
-        eq -> !any(isdiffed, 𝑠neighbors(graph, eq))
-    end
 
-    var_eq_matching = tear_graph_modia(structure, Union{Unassigned, SelectedState};
-                                       varfilter = can_eliminate,
-                                       eqfilter = should_consider)
+    var_eq_matching = tear_graph_modia(structure, isdiffed,
+                                       Union{Unassigned, SelectedState};
+                                       varfilter = can_eliminate)
     for v in eachindex(var_eq_matching)
         can_eliminate(v) && continue
         var_eq_matching[v] = SelectedState()