Fix non-deterministic behavior and incorrect variable selection in Carpanzano tearing algorithm

src/carpanzano_tearing.jl

@@ -52,8 +52,8 @@ function (alg::CarpanzanoTearing)(structure::SystemStructure)
     full_var_eq_matching = copy(var_eq_matching)
     var_sccs = find_var_sccs(graph, var_eq_matching)
 
-    active_vars = Set{Int}()
-    active_eqs = Set{Int}()
+    active_vars = OrderedSet{Int}()
+    active_eqs = OrderedSet{Int}()
     for vars in var_sccs
         for var in vars
             # Identify variables and equations in this SCC
@@ -87,7 +87,7 @@ In the context of the paper, `structure.graph` is the associated bipartite graph
 """
 function find_single_solvable_eq!(
         structure::SystemStructure, var_eq_matching::MatchingT,
-        active_vars::Set{Int}, active_eqs::Set{Int}, condition::F = _ -> true;
+        active_vars::AbstractSet{Int}, active_eqs::AbstractSet{Int}, condition::F = _ -> true;
         nbors_buffer::Vector{Int} = Int[]
     ) where {F}
     (; graph, solvable_graph) = structure
@@ -116,7 +116,7 @@ all of `active_vars` to `unassigned`, and will be modified to match solvable var
 """
 function carpanzano_tear_scc!(
         alg::CarpanzanoTearing, structure::SystemStructure, var_eq_matching::MatchingT,
-        active_vars::Set{Int}, active_eqs::Set{Int}
+        active_vars::AbstractSet{Int}, active_eqs::AbstractSet{Int}
     )
     # TODO: This is an implementation of algorithm A1 in the paper. Find an efficient
     # way to implement algorithm A2 and analyze the benefits.
@@ -211,6 +211,8 @@ function carpanzano_tear_scc!(
             solvable_cnt = count(in(active_eqs), 𝑑neighbors(solvable_graph, ivar))
             if iszero(alg_var) || cnt > max_incidence_cnt || cnt == max_incidence_cnt && solvable_cnt < min_solvable_cnt
                 alg_var = ivar
+                max_incidence_cnt = cnt
+                min_solvable_cnt = solvable_cnt
             end
         end
 

test/carpanzano_tearing.jl

@@ -0,0 +1,61 @@
+# Tests for carpanzano_tearing.jl - determinism and Heuristic 2 correctness
+
+using BipartiteGraphs
+import Graphs: add_edge!
+
+struct TestSystemStructure <: StateSelection.SystemStructure
+    graph::BipartiteGraph{Int,Nothing}
+    solvable_graph::BipartiteGraph{Int,Nothing}
+    var_to_diff::DiffGraph
+    eq_to_diff::DiffGraph
+end
+
+@testset "carpanzano_tearing" begin
+    @testset "Heuristic 2 picks variable with maximum incidence" begin
+        # Graph: 3 equations (source vertices), 3 variables (destination vertices).
+        # Every edge is also present in the solvable_graph (all variables solvable in
+        # all equations they appear in), so Heuristic 1 never finds a non-solvable
+        # variable and falls through to Heuristic 2.
+        #
+        # Incidence of each variable (number of equations it appears in):
+        #   v1 (idx 1): eq1, eq2        => 2
+        #   v2 (idx 2): eq1, eq2, eq3   => 3   # maximum -> must be torn (algebraic)
+        #   v3 (idx 3): eq1, eq3        => 2
+        #
+        # Before the fix, max_incidence_cnt was never updated inside the Heuristic 2
+        # loop, so every variable satisfied cnt > typemin(Int) and the last variable
+        # in hash-iteration order was chosen -- non-deterministic and wrong.
+        graph = BipartiteGraph(3, 3)
+        for (eq, var) in [(1,1),(1,2),(1,3),(2,1),(2,2),(3,2),(3,3)]
+            add_edge!(graph, eq, var)
+        end
+        solvable_graph = BipartiteGraph(3, 3)
+        for (eq, var) in [(1,1),(1,2),(1,3),(2,1),(2,2),(3,2),(3,3)]
+            add_edge!(solvable_graph, eq, var)
+        end
+
+        structure = TestSystemStructure(
+            graph, solvable_graph,
+            DiffGraph(3), DiffGraph(3)
+        )
+        alg = StateSelection.CarpanzanoTearing()
+        # complete() makes invview(var_eq_matching) available, which the algorithm uses
+        # internally to find the variable matched to a solved equation.
+        var_eq_matching = complete(
+            Matching{Union{Unassigned, SelectedState}}(3), 3
+        )
+        active_vars = Set{Int}([1, 2, 3])
+        active_eqs  = Set{Int}([1, 2, 3])
+
+        StateSelection.carpanzano_tear_scc!(
+            alg, structure, var_eq_matching, active_vars, active_eqs
+        )
+
+        # v2 has the highest incidence (3 equations) -- Heuristic 2 must choose it as
+        # the algebraic (torn) variable, leaving it unmatched.
+        @test var_eq_matching[2] === unassigned
+        # v1 and v3 must each be matched to some equation.
+        @test var_eq_matching[1] isa Int
+        @test var_eq_matching[3] isa Int
+    end
+end

test/runtests.jl

@@ -4,6 +4,7 @@ using SparseArrays
 using Test
 
 include("bareiss.jl")
+include("carpanzano_tearing.jl")
 
 @testset "`get_new_mm`" begin
     mm = SSel.CLIL.SparseMatrixCLIL(