From 90fc5f8966cbf31ece4f75b374cf2961628d0f6c Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 11 Jun 2026 13:42:18 +0000 Subject: [PATCH 1/2] Select dummy derivatives on merged Mattsson-Soderlind blocks, not matching-SCCs The SCCs of the full Pantelides matching can be strictly finer than the blocks of the (differentiated equations) x (highest-derivative candidates) subproblem on which Mattsson-Soderlind selection is posed. A differentiated equation matched in one SCC but incident to candidate variables in another (e.g. a twice-differentiated connection alias) creates a singleton SCC whose candidate is demoted unconditionally, making state_priority silently ineffective and potentially forcing a state realization with singularities. Merge such SCCs with union-find before selection so the priority-sorted greedy demotion sees the full block. Also permute the integer Jacobian's columns when the candidates are priority-sorted, so bareiss col_order indexes the sorted variable list consistently. Fixes #101. Fixes #102. Co-Authored-By: Claude Fable 5 --- src/partial_state_selection.jl | 84 ++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/src/partial_state_selection.jl b/src/partial_state_selection.jl index f4eb3ea..d927328 100644 --- a/src/partial_state_selection.jl +++ b/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 From f19b11a80590fce7cdb4dd9be6194aa0d7ce81f7 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Thu, 11 Jun 2026 14:12:18 +0000 Subject: [PATCH 2/2] test: dummy-derivative block merging and Jacobian column permutation Co-Authored-By: Claude Fable 5 --- test/dummy_derivative_blocks.jl | 132 ++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 133 insertions(+) create mode 100644 test/dummy_derivative_blocks.jl diff --git a/test/dummy_derivative_blocks.jl b/test/dummy_derivative_blocks.jl new file mode 100644 index 0000000..0ed117f --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 1a03238..3a0dd88 100644 --- a/test/runtests.jl +++ b/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(