Select dummy derivatives on merged Mattsson-Söderlind blocks so state_priority can act across matching-SCCs

src/partial_state_selection.jl

@@ -187,6 +187,86 @@ struct DummyDerivativeSummary
     state_priority::Vector{Vector{Float64}}
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Merge the SCCs of the Pantelides matching into the blocks on which dummy-derivative
+selection must operate.
+
+The Mattsson–Söderlind dummy-derivative selection problem is posed on the subproblem of
+differentiated equations and their highest-derivative candidate variables. The SCCs of
+the full matching can be strictly finer than the blocks of that subproblem: a
+differentiated equation matched to a candidate in one SCC may be incident to candidate
+variables in other SCCs. A common example is a twice-differentiated connection alias
+`0 ~ D(D(x)) - D(D(y))` matched to `D(D(x))`, with `D(D(y))` belonging to a
+kinematic-loop SCC: `D(D(x))` then sits in a singleton SCC where it is demoted
+unconditionally, and a high `state_priority` on `x` cannot prevent it even though
+demoting `D(D(y))` instead would be structurally valid (see issue #101). Selecting per
+merged block restores the full selection freedom of the subproblem, and the
+priority-sorted greedy selection inside `dummy_derivative_graph!` then maximizes the
+total priority of the kept states.
+
+Returns the merged list of variable blocks; SCCs without coupling are returned
+unchanged (in particular the result is `===` the input when nothing merges).
+"""
+function merge_dummy_derivative_blocks(
+        structure::SystemStructure, var_eq_matching, var_sccs::Vector{Vector{Int}})
+    (; eq_to_diff, var_to_diff, graph) = structure
+    diff_to_eq = invview(eq_to_diff)
+    diff_to_var = invview(var_to_diff)
+
+    # SCC index of every variable that is a dummy-derivative candidate of its SCC,
+    # mirroring the candidate filter in `dummy_derivative_graph!`.
+    scc_of_candidate = zeros(Int, ndsts(graph))
+    for (i, vars) in enumerate(var_sccs), var in vars
+        var_eq_matching[var] isa Int || continue
+        (diff_to_var[var] !== nothing && is_present(structure, var)) || continue
+        scc_of_candidate[var] = i
+    end
+
+    # Union-find over SCC indices, merging along differentiated equations that are
+    # incident to candidate variables outside the SCC they are matched in.
+    parent = collect(1:length(var_sccs))
+    function root(i::Int)
+        while parent[i] != i
+            parent[i] = parent[parent[i]]
+            i = parent[i]
+        end
+        i
+    end
+    merged_any = false
+    for (i, vars) in enumerate(var_sccs), var in vars
+        eq = var_eq_matching[var]
+        eq isa Int || continue
+        diff_to_eq[eq] === nothing && continue
+        for var2 in 𝑠neighbors(graph, eq)
+            j = scc_of_candidate[var2]
+            (j == 0 || j == i) && continue
+            ri = root(i)
+            rj = root(j)
+            ri == rj && continue
+            # union by min keeps roots at the first SCC of each block, which
+            # preserves the original SCC order in the output
+            parent[max(ri, rj)] = min(ri, rj)
+            merged_any = true
+        end
+    end
+    merged_any || return var_sccs
+
+    buckets = Dict{Int, Vector{Int}}()
+    order = Int[]
+    for (i, vars) in enumerate(var_sccs)
+        r = root(i)
+        b = get!(buckets, r) do
+            push!(order, r)
+            Int[]
+        end
+        append!(b, vars)
+    end
+    sort!(order)
+    return [buckets[r] for r in order]
+end
+
 """
     $TYPEDSIGNATURES
 
@@ -227,6 +307,7 @@ function dummy_derivative_graph!(
     end
 
     var_sccs = find_var_sccs(graph, var_eq_matching)
+    var_sccs = merge_dummy_derivative_blocks(structure, var_eq_matching, var_sccs)
     var_perm = Int[]
     var_dummy_scc = Vector{Int}[]
     var_state_priority = Vector{Float64}[]
@@ -281,6 +362,9 @@ function dummy_derivative_graph!(
                 sortperm!(var_perm, sp)
                 permute!(vars, var_perm)
                 permute!(sp, var_perm)
+                # keep the Jacobian columns aligned with the permuted variable
+                # order; `col_order` below indexes into `vars` (#102)
+                J === nothing || (J = J[:, var_perm])
                 push!(var_dummy_scc, copy(vars))
                 push!(var_state_priority, sp)
             end

test/dummy_derivative_blocks.jl

@@ -0,0 +1,132 @@
+# Tests for merge_dummy_derivative_blocks (#101) and the Jacobian column
+# permutation in dummy_derivative_graph! (#102).
+
+using BipartiteGraphs
+import Graphs: add_edge!
+
+struct DDBlockTestStructure <: StateSelection.SystemStructure
+    graph::BipartiteGraph{Int, Nothing}
+    solvable_graph::BipartiteGraph{Int, Nothing}
+    var_to_diff::StateSelection.DiffGraph
+    eq_to_diff::StateSelection.DiffGraph
+end
+
+# A tearing algorithm probe that simply returns the selected dummy-derivative
+# set, so the tests can assert on the selection without running tearing.
+struct DummySetProbe <: StateSelection.TearingAlgorithm end
+(::DummySetProbe)(structure, dummy_derivatives) = (dummy_derivatives, (;))
+
+function ddblock_structure(neqs, nvars, edges, var_diffs, eq_diffs)
+    graph = BipartiteGraph(neqs, nvars)
+    solvable_graph = BipartiteGraph(neqs, nvars)
+    for (eq, var) in edges
+        add_edge!(graph, eq, var)
+        add_edge!(solvable_graph, eq, var)
+    end
+    var_to_diff = StateSelection.DiffGraph(nvars, true)
+    for (v, dv) in var_diffs
+        var_to_diff[v] = dv
+    end
+    eq_to_diff = StateSelection.DiffGraph(neqs, true)
+    for (e, de) in eq_diffs
+        eq_to_diff[e] = de
+    end
+    DDBlockTestStructure(graph, solvable_graph, var_to_diff, eq_to_diff)
+end
+
+@testset "dummy derivative selection on merged blocks (#101)" begin
+    # Distilled from a multibody arm: a high-priority coordinate chain
+    # x -> D(x) -> D²(x) rigidly aliased to a low-priority chain
+    # y -> D(y) -> D²(y) (a frame angle), with the dynamics formulated in y.
+    #
+    # Variables: 1: x, 2: D(x), 3: D²(x), 4: y, 5: D(y), 6: D²(y)
+    # Equations: 1: 0 ~ x - y         (alias)
+    #            2: 0 ~ D(x) - D(y)   (alias′)
+    #            3: 0 ~ D²(x) - D²(y) (alias″)
+    #            4: 0 ~ D²(y) - f(y)  (dynamics)
+    #
+    # Pantelides-style matching: D²(x) ↔ eq 3, D²(y) ↔ eq 4. The matching-induced
+    # SCCs are then singletons; per-SCC selection demotes D²(x) through eq 3
+    # unconditionally and the priority of x can never act, even though demoting
+    # D²(y) through eq 3 instead is structurally valid. The blocks must be merged
+    # so that the selection sees both candidates.
+    structure = ddblock_structure(
+        4, 6,
+        [(1, 1), (1, 4), (2, 2), (2, 5), (3, 3), (3, 6), (4, 6), (4, 4)],
+        [1 => 2, 2 => 3, 4 => 5, 5 => 6],
+        [1 => 2, 2 => 3])
+
+    make_matching = function ()
+        m = Matching(6)
+        m[3] = 3
+        m[6] = 4
+        complete(m, 4)
+    end
+
+    # x (and via the derivative chain, D²(x)) has high priority
+    state_priority_x = v -> v == 1 ? 100.0 : 0.0
+
+    # the singleton SCCs of the matching must merge into one block
+    var_eq_matching = make_matching()
+    sccs = StateSelection.find_var_sccs(structure.graph, var_eq_matching)
+    merged = StateSelection.merge_dummy_derivative_blocks(structure, var_eq_matching, sccs)
+    blocks = [b for b in merged if 3 in b || 6 in b]
+    @test length(blocks) == 1
+    @test 3 in blocks[1] && 6 in blocks[1]
+
+    # end-to-end: selection must demote the low-priority chain (D(y), D²(y)),
+    # keeping the priority-100 chain of x as states
+    dummys, _ = StateSelection.dummy_derivative_graph!(
+        structure, make_matching(), nothing, state_priority_x, Val(false);
+        tearing_alg = DummySetProbe())
+    @test dummys == BitSet([5, 6])
+
+    # with the priorities swapped, the selection must flip
+    state_priority_y = v -> v == 4 ? 100.0 : 0.0
+    dummys2, _ = StateSelection.dummy_derivative_graph!(
+        structure, make_matching(), nothing, state_priority_y, Val(false);
+        tearing_alg = DummySetProbe())
+    @test dummys2 == BitSet([2, 3])
+end
+
+@testset "Jacobian columns follow the priority permutation (#102)" begin
+    # Three integrator chains x, y, z whose derivatives Dx, Dy, Dz form one SCC
+    # through the matching, with two differentiated equations carrying an
+    # all-integer Jacobian:
+    #
+    # Variables: 1: x, 2: D(x), 3: y, 4: D(y), 5: z, 6: D(z)
+    # Equations: 1: a0(x, y)            (integral of eq 2)
+    #            2: a(D(x), D(y))       differentiated, ∂/∂D(x) = 1, ∂/∂D(y) = 0
+    #            3: b0(z, x)            (integral of eq 4)
+    #            4: b(D(z), D(x))       differentiated, ∂/∂D(z) = 1, ∂/∂D(x) = 0
+    #            5: c(D(y), D(z))       algebraic
+    #
+    # Matching: D(x) ↔ eq 2, D(y) ↔ eq 5, D(z) ↔ eq 4 puts {2, 4, 6} in one SCC.
+    # y has high state priority, so the two demotions (for eqs 2 and 4) must fall
+    # on D(x) and D(z). Sorted candidate order: [D(x), D(z), D(y)]; the permuted
+    # Jacobian has pivots in columns 1 and 2. Without permuting the Jacobian
+    # alongside the candidates, bareiss pivots on the unsorted columns
+    # ([D(x), D(y), D(z)] with a zero column for D(y)) and reports column order
+    # (1, 3, 2), so the selection demotes sorted[3] = D(y) — the high-priority
+    # variable — and keeps the structurally unsolvable D(z) as a state.
+    structure = ddblock_structure(
+        5, 6,
+        [(1, 1), (1, 3), (2, 2), (2, 4), (3, 5), (3, 1), (4, 6), (4, 2), (5, 4), (5, 6)],
+        [1 => 2, 3 => 4, 5 => 6],
+        [1 => 2, 3 => 4])
+
+    m = Matching(6)
+    m[2] = 2
+    m[4] = 5
+    m[6] = 4
+    var_eq_matching = complete(m, 5)
+
+    state_priority = v -> v == 3 ? 100.0 : 0.0
+    jac = (eqs, vars) -> [((eq == 2 && var == 2) || (eq == 4 && var == 6)) ? 1 : 0
+                          for eq in eqs, var in vars]
+
+    dummys, _ = StateSelection.dummy_derivative_graph!(
+        structure, var_eq_matching, jac, state_priority, Val(false);
+        tearing_alg = DummySetProbe())
+    @test dummys == BitSet([2, 6])
+end

test/runtests.jl

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