diff --git a/lib/ModelingToolkitTearing/Project.toml b/lib/ModelingToolkitTearing/Project.toml index 0676749..ac0d1ce 100644 --- a/lib/ModelingToolkitTearing/Project.toml +++ b/lib/ModelingToolkitTearing/Project.toml @@ -1,6 +1,6 @@ name = "ModelingToolkitTearing" uuid = "6bb917b9-1269-42b9-9f7c-b0dca72083ab" -version = "1.14.1" +version = "1.14.2" authors = ["Aayush Sabharwal "] [deps] @@ -13,6 +13,7 @@ ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839" Moshi = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -34,6 +35,7 @@ ModelingToolkitBase = "1.37" Moshi = "0.3" OffsetArrays = "1" OrderedCollections = "1.8.1" +Random = "1" SciMLBase = "2.108, 3" Setfield = "0.7, 0.8, 1" SparseArrays = "1" @@ -47,7 +49,8 @@ julia = "1.10" [extras] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "ModelingToolkit", "ForwardDiff"] +test = ["Test", "ModelingToolkit", "ForwardDiff", "SparseArrays"] diff --git a/lib/ModelingToolkitTearing/src/ModelingToolkitTearing.jl b/lib/ModelingToolkitTearing/src/ModelingToolkitTearing.jl index 2efef62..8c2c18d 100644 --- a/lib/ModelingToolkitTearing/src/ModelingToolkitTearing.jl +++ b/lib/ModelingToolkitTearing/src/ModelingToolkitTearing.jl @@ -27,6 +27,7 @@ using SymbolicUtils: BSImpl, unwrap using SciMLBase: LinearProblem using SparseArrays: nonzeros import LinearAlgebra +import Random import UUIDs: UUID, uuid4 const TimeDomain = SciMLBase.AbstractClock diff --git a/lib/ModelingToolkitTearing/src/reassemble.jl b/lib/ModelingToolkitTearing/src/reassemble.jl index ebdb4cd..1dc7ef3 100644 --- a/lib/ModelingToolkitTearing/src/reassemble.jl +++ b/lib/ModelingToolkitTearing/src/reassemble.jl @@ -426,22 +426,15 @@ function generate_system_equations!(state::TearingState, neweqs::Vector{Equation end digraph = DiCMOBiGraph{false}(graph, var_eq_matching) - for (i, scc) in enumerate(var_sccs) - # note that the `vscc <-> escc` relation is a set-to-set mapping, and not - # point-to-point. - vscc, escc = get_sorted_scc(digraph, full_var_eq_matching, var_eq_matching, scc) - var_sccs[i] = vscc - if length(escc) != length(vscc) - isempty(escc) && continue - escc = setdiff(escc, extra_eqs) - isempty(escc) && continue - vscc = setdiff(vscc, extra_vars) - isempty(vscc) && continue - end + # Emit one block (a single SCC, or the union of a jointly-solved family) as either an + # inline linear solve or, failing that, regular per-equation codegen. Returns `true` + # if the block was emitted; with `allow_regular = false` it emits only via the inline + # linear path and returns `false` when that path does not apply (so the caller can + # fall back to emitting the family's members individually). + emit_block! = function (vscc, escc; allow_regular::Bool = true) # Inline linear SCCs pass is only valid on continuous systems. We check if the - # current SCC is algebraic and if the algebraic equations are linear in the - # algebraic variables. + # block is algebraic and if the algebraic equations are linear in the variables. linsol_result = nothing if !is_disc && inline_linear_sccs linsol_result = get_linear_scc_linsol(state, escc, vscc, neweqs, var_eq_matching, total_sub, analytical_linear_scc_limit, simplify) @@ -489,12 +482,71 @@ function generate_system_equations!(state::TearingState, neweqs::Vector{Equation var = eq_var_matching[ieq]::Int codegen_equation!(eq_generator, neweqs[ieq], ieq, var; simplify) end - else + return true + elseif allow_regular for ieq in escc iv = eq_var_matching[ieq] neq = neweqs[ieq] codegen_equation!(eq_generator, neq, ieq, iv; simplify) end + return true + else + # Inline path inapplicable: surface the reason (`NonlinearBlockEqs` carrying + # the offending equations, or `nothing`) so the caller can peel and retry. + return linsol_result + end + end + + # Sort each SCC (and trim the extra equations/variables) into the emittable blocks, + # preserving the block-triangular (topological) order. + prepared = NTuple{2, Vector{Int}}[] + for (i, scc) in enumerate(var_sccs) + # note that the `vscc <-> escc` relation is a set-to-set mapping, and not + # point-to-point. + vscc, escc = get_sorted_scc(digraph, full_var_eq_matching, var_eq_matching, scc) + var_sccs[i] = vscc + if length(escc) != length(vscc) + isempty(escc) && continue + escc = setdiff(escc, extra_eqs) + isempty(escc) && continue + vscc = setdiff(vscc, extra_vars) + isempty(vscc) && continue + end + push!(prepared, (vscc, escc)) + end + + # Group coupled, runtime-rank-deficible inline-linear SCCs into families that are + # solved jointly (issue #98): when a sequentially-solved block becomes rank-deficient + # at a degenerate parameter point, its gauge choice can make a downstream block + # inconsistent even though the union of the family's equations is satisfiable. Only + # active when `maybe_zeros` is set and inline linear SCCs are enabled, so the default + # behaviour (one block per SCC) is unchanged. + groups = _group_inline_linear_families( + state, prepared, MTKBase.maybe_zeros(state.sys), neweqs, graph, is_disc, inline_linear_sccs) + if _famgrp_log_enabled() + fams = [grp for grp in groups if length(grp) > 1] + isempty(fams) || @info "Inline-linear-SCC families formed" nblocks=length(prepared) nfamilies=length(fams) family_equation_counts=[sum(length(prepared[k][2]) for k in grp) for grp in fams] + end + for grp in groups + if length(grp) == 1 + vscc, escc = prepared[grp[1]] + emit_block!(vscc, escc) + continue + end + # Solve the family jointly. The grouping pre-screens members for linearity + # in the union's variables and keeps families path-closed (convex), so a + # failure here is rare (e.g. linearity lost through `total_sub` + # substitution); fall back to emitting each member as its own block. + vscc = reduce(vcat, (prepared[k][1] for k in grp)) + escc = reduce(vcat, (prepared[k][2] for k in grp)) + res = emit_block!(vscc, escc; allow_regular = false) + if res !== true + if _inline_scc_check_enabled() + @warn "Inline-linear-SCC family joint solve fell back to per-member emission" family_blocks=length(grp) reason=res + end + for k in grp + emit_block!(prepared[k]...) + end end end @@ -538,6 +590,280 @@ function generate_system_equations!(state::TearingState, neweqs::Vector{Equation length(solved_vars_set) end +""" + $TYPEDSIGNATURES + +Group the prepared inline-linear blocks (`prepared[k] = (vscc, escc)`, in block-triangular +order) into the units in which they are emitted, returning a `Vector{Vector{Int}}` of index +groups in a valid (topological) emission order. Most groups are singletons `[k]`; a +multi-element group is a *family* that must be solved as one joint linear system. + +The inline-linear-SCC pass solves blocks sequentially, substituting each block's solution +into the downstream blocks. At a degenerate parameter point a block can become +rank-deficient; its (gauge) solution choice is then substituted downstream and can make a +dependent block inconsistent even though the *union* of the family's equations is +satisfiable (issue #98). To avoid this we keep such families unsubstituted and solve each +jointly. + +Families are computed on the block-level dependency DAG (edge `j → k` iff block `k`'s +equations structurally reference block `j`'s variables; members of a family are typically +*not* adjacent in the block-triangular order): + +1. A block is *runtime-rank-deficible* if its equations reference a `maybe_zeros` + parameter, so a coefficient can vanish at runtime and drop the rank. +2. Rank-deficible blocks that reach each other through the DAG (possibly through + intermediate blocks) belong to the same family. +3. Each family is closed over the blocks lying on dependency paths between its members, so + the merged system is self-contained: every variable it references is either solved + upstream or part of the merged block. + +The groups are returned in a topological order of the DAG with each family contracted to a +single node (families are path-closed, so the contraction cannot create cycles). When +`maybe_zeros` is empty no grouping happens at all, leaving the default one-block-per-SCC +behaviour — including the emission order — unchanged. +""" +function _group_inline_linear_families( + state::TearingState, prepared::Vector{NTuple{2, Vector{Int}}}, maybe_zeros, + neweqs::Vector{Equation}, graph, is_disc::Bool, inline_linear_sccs::Bool) + n = length(prepared) + singletons() = [[k] for k in 1:n] + (is_disc || !inline_linear_sccs) && return singletons() + (maybe_zeros === nothing || isempty(maybe_zeros)) && return singletons() + mzs = collect(maybe_zeros) + + references_maybe_zeros(ex) = any(Symbolics.get_variables(unwrap(ex))) do v + base = MTKBase.split_indexed_var(v)[1] + any(z -> isequal(base, unwrap(z)), mzs) + end + # A block's rank can drop iff one of its equations involves a `maybe_zeros` parameter. + droppable = falses(n) + for k in 1:n + _, escc = prepared[k] + droppable[k] = any(escc) do ieq + eq = neweqs[ieq] + res = SU._iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs + references_maybe_zeros(res) + end + end + any(droppable) || return singletons() + + # Block-level dependency DAG: edge `j → k` iff block `k`'s equations structurally + # reference block `j`'s variables. + var_block = Dict{Int, Int}() + for (k, (vscc, _)) in enumerate(prepared), v in vscc + var_block[v] = k + end + + # Whether all equations of block `b` are linear in every incident variable + # belonging to block `j` — i.e. whether `b` tolerates `j` in a jointly-solved + # linear family. Linearity is *set-dependent*: a block may be nonlinear in some + # other block's variables (e.g. cos of another block's angle) and still be + # perfectly linear in its own and its family's variables, so this is checked + # pairwise over a candidate family's closure rather than globally. Memoized. + (; fullvars, sys) = state + pairlinear_cache = Dict{Tuple{Int, Int}, Bool}() + blockpair_linear = function (b, j) + get!(pairlinear_cache, (b, j)) do + _, escc_b = prepared[b] + for ieq in escc_b + eq = neweqs[ieq] + res = SU._iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs + for v in 𝑠neighbors(graph, ieq) + get(var_block, v, 0) == j || continue + lex = MTKBase.get_linear_expander_for!(sys, fullvars[v], true) + _, _, islin = lex(res) + islin || return false + end + end + return true + end + end + closure_linear = function (cl) + for b in cl, j in cl + b == j && continue + blockpair_linear(b, j) || return false + end + return true + end + dag = Graphs.SimpleDiGraph(n) + for (k, (_, escc)) in enumerate(prepared), ieq in escc, v in 𝑠neighbors(graph, ieq) + j = get(var_block, v, 0) + (j == 0 || j == k) && continue + # Forward edges only: the prepared blocks are in BLT order, which is a valid + # linearization of the true (matching-based) dependencies. Raw residual + # incidence also contains backward references (e.g. through torn variables + # of later blocks); treating those as dependencies would make this graph + # cyclic even though sequential emission is perfectly well-defined. + j < k || continue + Graphs.add_edge!(dag, j, k) + end + + bfs = function (sources, neighborfn) + seen = falses(n) + stack = collect(Int, sources) + seen[stack] .= true + while !isempty(stack) + u = pop!(stack) + for w in neighborfn(dag, u) + seen[w] && continue + seen[w] = true + push!(stack, w) + end + end + return seen + end + + # Rank-deficible blocks that reach one another (possibly through intermediate blocks) + # belong to the same family. + dropidx = findall(droppable) + reach = Dict{Int, BitVector}(d => bfs((d,), Graphs.outneighbors) for d in dropidx) + conn = Graphs.SimpleGraph(n) + for d1 in dropidx, d2 in dropidx + d1 < d2 || continue + (reach[d1][d2] || reach[d2][d1]) && Graphs.add_edge!(conn, d1, d2) + end + + # Close each family over the blocks on dependency paths between its members, so the + # merged system is self-contained. Closures of distinct families cannot overlap: a + # rank-deficible block on a path between two members is reachability-connected to + # them and thus in the same component, and a shared non-deficible block would chain + # the two components together. + # Full-graph BFS (ignores mergeability): used to compute the convex (path) + # closure of a candidate family in the *full* DAG. Convexity there is a + # correctness requirement — a family must absorb every block lying on a + # dependency path between its members, since such blocks both consume family + # outputs and feed family inputs, and a joint solve that treats them as + # external inputs would be circular. + bfs_full = function (sources, neighborfn) + seen = falses(n) + stack = collect(Int, sources) + seen[stack] .= true + while !isempty(stack) + u = pop!(stack) + for w in neighborfn(dag, u) + seen[w] && continue + seen[w] = true + push!(stack, w) + end + end + return seen + end + full_closure(members) = bfs_full(members, Graphs.outneighbors) .& + bfs_full(members, Graphs.inneighbors) + + family_of = zeros(Int, n) + nfam = 0 + finalize_family! = function (cur) + length(cur) > 1 || return + closure = full_closure(cur) + # Closures of distinct families must not overlap (every block is emitted + # exactly once); skip if a previous family already claimed part of it. + if any(j -> family_of[j] != 0, findall(closure)) + _famgrp_log_enabled() && + println("FAMGREEDY family skipped (closure overlap): members=", cur) + return + end + nfam += 1 + family_of[closure] .= nfam + end + for comp in Graphs.connected_components(conn) + length(comp) > 1 || continue + # Grow families greedily along the topological order: extend while the + # full-DAG closure of the candidate set remains entirely mergeable, so + # each finalized family is convex AND jointly linearly solvable. Members + # whose inclusion would drag in a non-mergeable block (e.g. the chassis + # block far downstream, with nonlinear kinematics on the intervening + # paths) start a new family instead. + members = sort(comp) + cur = Int[members[1]] + for m in Iterators.drop(members, 1) + cand = vcat(cur, m) + closure = full_closure(cand) + cl = findall(closure) + if all(k -> family_of[k] == 0, cl) && closure_linear(cl) + cur = cand + else + _famgrp_log_enabled() && + println("FAMGREEDY cannot extend ", cur, " by ", m, ": closure=", length(cl)) + finalize_family!(cur) + cur = Int[m] + end + end + finalize_family!(cur) + end + _famgrp_log_enabled() && + (println("FAMGREEDY formed nfam=", nfam, " sizes=", [count(==(f), family_of) for f in 1:nfam]); flush(stdout)) + iszero(nfam) && return singletons() + + # Emit in a topological order of the DAG with each family contracted to one node. + # Even disjoint convex families need not be mutually orderable: contracting two + # families that interleave in the block order can create a cycle through the + # blocks between them. When Kahn's algorithm stalls, dissolve one family that is + # stuck in the residual cycle and retry, keeping a maximal orderable subset. + while true + unit_id(k) = family_of[k] == 0 ? nfam + k : family_of[k] + unit_blocks = Dict{Int, Vector{Int}}() + for k in 1:n + push!(get!(Vector{Int}, unit_blocks, unit_id(k)), k) + end + indeg = Dict{Int, Int}(u => 0 for u in keys(unit_blocks)) + outs = Dict{Int, Set{Int}}(u => Set{Int}() for u in keys(unit_blocks)) + for e in Graphs.edges(dag) + uj = unit_id(Graphs.src(e)) + uk = unit_id(Graphs.dst(e)) + (uj == uk || uk in outs[uj]) && continue + push!(outs[uj], uk) + indeg[uk] += 1 + end + # Kahn's algorithm, preferring the unit containing the smallest original block + # index to keep the emission order stable. + minidx = Dict{Int, Int}(u => first(bs) for (u, bs) in unit_blocks) + ready = [u for (u, d) in indeg if d == 0] + groups = Vector{Int}[] + emitted_units = Set{Int}() + while !isempty(ready) + bestpos = 1 + for p in 2:length(ready) + minidx[ready[p]] < minidx[ready[bestpos]] && (bestpos = p) + end + u = ready[bestpos] + deleteat!(ready, bestpos) + push!(groups, unit_blocks[u]) + push!(emitted_units, u) + for w in outs[u] + (indeg[w] -= 1) == 0 && push!(ready, w) + end + end + sum(length, groups; init = 0) == n && return groups + # Stalled: some families are part of a contracted cycle. Dissolve the stalled + # family with the largest first-block index (preserving the earliest ones). + stalled = [f for f in 1:nfam if !(f in emitted_units) && any(==(f), family_of)] + if isempty(stalled) + # Cycle without any family involved cannot happen on a BLT block list; + # be defensive anyway. + return singletons() + end + drop = argmax(f -> minidx[f], stalled) + _famgrp_log_enabled() && + println("FAMGREEDY dissolving family ", drop, " (unorderable against earlier families)") + family_of[family_of .== drop] .= 0 + any(!=(0), family_of) || return singletons() + end +end + +""" + $TYPEDEF + +Returned by [`get_linear_scc_linsol`](@ref) instead of `nothing` when the inline +linear solve is inapplicable because equations are nonlinear in the block's +variables. Carries the global indices of all offending equations so a +jointly-solved family can peel the member blocks owning them and retry in a +single step (issue #98). +""" +struct NonlinearBlockEqs + ieqs::Vector{Int} +end + const INLINE_LINEAR_SCC_OP = (\) """ @@ -591,12 +917,19 @@ function get_linear_scc_linsol(state::TearingState, alg_eqs::Vector{Int}, b[eqidx] = resid end + bad_eqs = Int[] for (varidx, var) in enumerate(vars) lex = MTKBase.get_linear_expander_for!(sys, var, true) for (eqidx, resid) in enumerate(b) Graphs.has_edge(graph, BipartiteEdge(alg_eqs[eqidx], alg_vars[varidx])) || continue + alg_eqs[eqidx] in bad_eqs && continue p, q, islinear = lex(resid) - islinear || return nothing + if !islinear + # Record and keep scanning so a jointly-solved family can peel all + # offending member blocks in one step rather than one per retry. + push!(bad_eqs, alg_eqs[eqidx]) + continue + end if !SU._iszero(p) # We're iterating in increasing `varidx` (column index) so we can just `push!` push!(A.row_cols[eqidx], varidx) @@ -605,6 +938,51 @@ function get_linear_scc_linsol(state::TearingState, alg_eqs::Vector{Int}, b[eqidx] = q end end + isempty(bad_eqs) || return NonlinearBlockEqs(bad_eqs) + + # The `has_edge` gate above relies on the structural incidence `graph`, which is + # mutated during reassembly and can desync from the `total_sub`-substituted residual. + # A false-negative edge leaves a live SCC variable buried inside `b[eqidx]`; if that + # row is later eliminated by `__reduce_linear_system!` (which only inspects the + # coefficient row, never `b`), the buried variable is smuggled onto the RHS of + # retained equations, producing an inconsistent runtime `A \ b` (issue #98). Restore + # the invariant "`b` is free of all SCC variables" by re-expanding any leftover term. + var_atoms = [Set{Any}(Symbolics.get_variables(unwrap(var))) for var in vars] + scc_atoms = union(Set{Any}(), var_atoms...) + repaired = 0 + if !isempty(scc_atoms) + for eqidx in 1:N + bsyms = Symbolics.get_variables(unwrap(b[eqidx])) + any(in(scc_atoms), bsyms) || continue + present = Set(A.row_cols[eqidx]) + did_repair = false + for (varidx, var) in enumerate(vars) + varidx in present && continue + # Only invoke the (cached) expander if this variable actually occurs. + isempty(intersect(var_atoms[varidx], bsyms)) && continue + lex = MTKBase.get_linear_expander_for!(sys, var, true) + p, q, islinear = lex(b[eqidx]) + islinear || return NonlinearBlockEqs([alg_eqs[eqidx]]) + b[eqidx] = q + if !SU._iszero(p) + push!(A.row_cols[eqidx], varidx) + push!(A.row_vals[eqidx], p) + did_repair = true + end + bsyms = Symbolics.get_variables(unwrap(b[eqidx])) + end + if did_repair + # The re-expanded columns may be out of order; `A.row_cols` must be sorted. + perm = sortperm(A.row_cols[eqidx]) + A.row_cols[eqidx] = A.row_cols[eqidx][perm] + A.row_vals[eqidx] = A.row_vals[eqidx][perm] + repaired += 1 + end + end + end + if repaired > 0 && _inline_scc_check_enabled() + @info "Inline-linear-SCC construction re-expanded buried SCC variables in `b`" block_n=N repaired_equations=repaired alg_vars + end # `-` is important! `b` is on the other side of the equality. for i in eachindex(b) @@ -706,8 +1084,167 @@ function get_linear_scc_linsol(state::TearingState, alg_eqs::Vector{Int}, return (INLINE_LINEAR_SCC_OP(A, b), eqs_mask, vars_mask) end +# --------------------------------------------------------------------------- +# Opt-in self-checks for the inline-linear-SCC pass (issue #98). +# +# These are disabled unless the `MTKTEARING_CHECK_REDUCTION` environment variable +# is set to a non-empty value. When off, the only added cost is a single `get(ENV, +# ...)` per SCC and the snapshots below are skipped, so production pays nothing. +# --------------------------------------------------------------------------- + +""" + $TYPEDSIGNATURES + +Whether the inline-linear-SCC self-checks are enabled. Controlled by the +`MTKTEARING_CHECK_REDUCTION` environment variable (any non-empty value enables it). +""" +_inline_scc_check_enabled() = !isempty(get(ENV, "MTKTEARING_CHECK_REDUCTION", "")) +# Separate, cheap gate for family-grouping progress logging (no snapshots/rank checks). +_famgrp_log_enabled() = _inline_scc_check_enabled() || !isempty(get(ENV, "MTKT_FAMGRP_LOG", "")) + +# Numerically evaluate a symbolic expression under a substitution of *all* its free +# symbols to numbers. Deliberately avoids `iszero`/`simplify`/`expand`, which can OOM +# on large multibody coefficient expressions (see StateSelection.jl#95). +_evalnum(x, subs::AbstractDict) = + Float64(Symbolics.value(Symbolics.substitute(unwrap(x), subs; fold = Val(true)))) + +_free_syms_into!(s::AbstractSet, x) = + (union!(s, Symbolics.get_variables(unwrap(x))); s) + +# Build a deterministic random substitution for all free symbols appearing in the +# given symbolic containers, plus an RNG seeded from the symbol *names* (so a reported +# failure reproduces across runs). Returns `(subs, rng)`. +function _deterministic_subs(containers...) + syms = Set{Any}() + for c in containers, x in c + _free_syms_into!(syms, x) + end + symvec = sort!(collect(syms); by = string) + seed = foldl((h, s) -> hash(string(s), h), symvec; init = UInt(0x5eed)) + rng = Random.MersenneTwister(seed % typemax(UInt) + one(UInt)) + # Bounded draws in (0.15, 0.85): generic enough for rank/identity probing while + # keeping common expression domains valid (sqrt(1 - x^2), log(x), 1/x, ...). + # Array-shaped symbols get an array of draws so indexing still folds. + draw() = 0.15 + 0.7 * rand(rng) + subs = Dict{Any, Any}() + for s in symvec + sh = SU.shape(unwrap(s)) + if sh isa SU.Unknown || isempty(sh) + subs[s] = draw() + else + dims = map(length, Tuple(sh)) + subs[s] = [draw() for _ in CartesianIndices(dims)] + end + end + return subs, rng +end + +""" + $TYPEDSIGNATURES + +Verify that the reduced linear system `A_red x_ret = b_red` produced by +`__reduce_linear_system!` is the *exact* substitution-projection of the full system +`A0 x = b0`, i.e. for every retained equation the reduced residual equals the full +residual after replacing each eliminated variable `v` by `aliases[v]·x_ret + +constants[v]`. Returns `true` iff the identity holds at a deterministic pseudo-random +probe point. Pure and positional: it does not need to know which symbols are the SCC +variables (those are represented by column position, not as free symbols). +""" +function _reduction_identity_ok( + A0::AbstractMatrix, b0::AbstractVector, A_red::AbstractMatrix, b_red::AbstractVector, + aliases::AbstractDict, constants::AbstractDict, eqs_mask::BitVector, + vars_mask::BitVector, old_to_new_eq::Vector{Int}; rtol::Float64 = 1e-7) + N = length(b0) + aliasvals = (SparseArrays.nonzeros(sv) for sv in values(aliases)) + subs, rng = _deterministic_subs(A0, b0, A_red, b_red, values(constants), aliasvals...) + + # New index of each retained variable (mirrors `cumsum(vars_mask)` in the caller). + old_to_new_var = zeros(Int, N) + let c = 0 + for j in 1:N + vars_mask[j] || continue + old_to_new_var[j] = (c += 1) + end + end + nret = count(vars_mask) + xret = randn(rng, nret) + + # Full-length solution: retained vars from `xret`, eliminated vars from their alias. + xfull = Vector{Float64}(undef, N) + for j in 1:N + vars_mask[j] && (xfull[j] = xret[old_to_new_var[j]]) + end + for j in 1:N + vars_mask[j] && continue + acc = _evalnum(constants[j], subs) + I, V = SparseArrays.findnz(aliases[j]) + for (k, coeff) in zip(I, V) + # After Phase-2 flattening, `k` references only retained variables. + acc += _evalnum(coeff, subs) * xfull[k] + end + xfull[j] = acc + end + + nr = length(b_red) + ok = true + for i in 1:N + eqs_mask[i] || continue + ired = old_to_new_eq[i] + rfull = sum(_evalnum(A0[i, j], subs) * xfull[j] for j in 1:N; init = 0.0) - _evalnum(b0[i], subs) + rred = sum(_evalnum(A_red[ired, j], subs) * xret[j] for j in 1:nr; init = 0.0) - _evalnum(b_red[ired], subs) + if abs(rfull - rred) > rtol * (1 + abs(rfull)) + ok = false + @warn "Inline-linear-SCC reduction identity violated" full_row=i reduced_row=ired residual_full=rfull residual_reduced=rred mismatch=abs(rfull - rred) + end + end + return ok +end + +# Report the rank/consistency of the full and reduced blocks at a deterministic probe +# point. Distinguishes "full block already bad (upstream construction)" from "reduction +# broke it". Logs via `@info`; only called when the self-check is enabled. +function _reduction_rank_report(A0::AbstractMatrix, b0::AbstractVector, + A_red::AbstractMatrix, b_red::AbstractVector, alg_vars::Vector{Int}) + subs, _ = _deterministic_subs(A0, b0, A_red, b_red) + N = length(b0) + A0n = [_evalnum(A0[i, j], subs) for i in 1:N, j in 1:N] + b0n = [_evalnum(b0[i], subs) for i in 1:N] + nr = length(b_red) + Arn = [_evalnum(A_red[i, j], subs) for i in 1:nr, j in 1:nr] + brn = [_evalnum(b_red[i], subs) for i in 1:nr] + # Rank-tolerant residual: min-norm least-squares via the pseudoinverse, so a + # (legitimately) rank-deficient block does not throw and we can still tell whether + # `b` lies in the range of `A` (relres ≈ 0 consistent, relres ≈ 1 inconsistent). + relres(A, b) = isempty(b) ? 0.0 : LinearAlgebra.norm(A * (LinearAlgebra.pinv(A) * b) - b) / max(LinearAlgebra.norm(b), eps()) + consistent(A, b) = LinearAlgebra.rank(A) == LinearAlgebra.rank(hcat(A, b)) + @info "Inline-linear-SCC block diagnostics" alg_vars full_n=N full_rank=LinearAlgebra.rank(A0n) full_consistent=consistent(A0n, b0n) full_relres=relres(A0n, b0n) reduced_n=nr reduced_rank=(nr == 0 ? 0 : LinearAlgebra.rank(Arn)) reduced_consistent=(nr == 0 ? true : consistent(Arn, brn)) reduced_relres=relres(Arn, brn) + return nothing +end + + +# Memoized unique-node count of a symbolic expression DAG, with early abort once +# `budget` is exceeded (returns a value > budget in that case). Used to bound the +# size of composed elimination expressions in `__reduce_linear_system!`. +function _node_count_capped(x, budget::Int, seen::IdDict{Any, Nothing} = IdDict{Any, Nothing}()) + x = unwrap(x) + SU.iscall(x) || return 1 + haskey(seen, x) && return 0 + seen[x] = nothing + n = 1 + for a in SU.arguments(x) + n += _node_count_capped(a, budget, seen) + n > budget && return n + end + return n +end + function __reduce_linear_system!(A::StateSelection.CLIL.SparseMatrixCLIL{Num, Int}, b::Vector{SymbolicT}, var_eq_matching::StateSelection.VarEqMatchingT, alg_eqs::Vector{Int}, alg_vars::Vector{Int}) N = length(b) + # Snapshot the pre-reduction system for the opt-in self-check. `A`'s rows are + # mutated in place below and `b` is reassigned, so these must be copies. + _check = _inline_scc_check_enabled() + A0_check = _check ? collect(A)::Matrix{Num} : nothing + b0_check = _check ? copy(b) : nothing # Identify rows (equations) not worth involving in the linear solve. # # The current heuristic is to find all rows with constant coefficients @@ -721,75 +1258,89 @@ function __reduce_linear_system!(A::StateSelection.CLIL.SparseMatrixCLIL{Num, In # `∑_k aliases[i][k] * fullvars[alg_vars[k]] + constants[i]` constants = Dict{Int, SymbolicT}() aliases = Dict{Int, SparseArrays.SparseVector{Num, Int}}() + # Pass 1: collect elimination candidates. A matched row is a candidate when its + # PIVOT coefficient is a known-nonzero constant; off-pivot coefficients may be + # symbolic and fold into the alias expression (StateSelection.jl#95). Requiring + # the whole row to be constant (the previous gate) keeps hundreds of + # sequentially-solvable rows inside the numeric solve. + cand_row = Dict{Int, Int}() # ivar => row index i + cand_piv = Dict{Int, Num}() # ivar => pivot coefficient for (i, coeffs) in enumerate(A.row_vals) - all(SU.isconst ∘ unwrap, coeffs) || continue eq_var_matching[alg_eqs[i]] isa Int || continue - - eqs_mask[i] = false var = eq_var_matching[alg_eqs[i]]::Int - ivar = var_to_idx[var] - vars_mask[ivar] = false - new_N -= 1 - - eqvars = A.row_cols[i] - idx_in_eq = findfirst(isequal(ivar), eqvars)::Int - var_coeff = coeffs[idx_in_eq] - deleteat!(eqvars, idx_in_eq) - deleteat!(coeffs, idx_in_eq) - # Negation to move variables to the other side of the equality - coeffs ./= -var_coeff - aliases[ivar] = SparseArrays.SparseVector(N, eqvars, coeffs) - # `b` is already on the other side of the equality - constants[ivar] = b[i] / var_coeff - end - - # We could have eliminated a variable in terms of other eliminated variables. While - # `get_new_mm` can handle this, it makes updating `b` much more difficult. We can - # topologically sort the dependency graph and use this information to update `aliases` - # and `constants` to fix this issue. + ivar = get(var_to_idx, var, nothing) + ivar === nothing && continue + idx_in_eq = findfirst(isequal(ivar), A.row_cols[i]) + idx_in_eq === nothing && continue + var_coeff = A.row_vals[i][idx_in_eq] + SU.isconst(unwrap(var_coeff)) || continue + SU._iszero(unwrap(var_coeff)) && continue + cand_row[ivar] = i + cand_piv[ivar] = var_coeff + end + + # Pass 2: topologically order the candidates along their dependency DAG (a + # candidate row may reference other candidate variables; tearing guarantees + # acyclicity). dep_graph = Graphs.SimpleDiGraph(length(alg_vars)) - for (var, coeffs) in aliases - I, _ = SparseArrays.findnz(coeffs) - for other_var in I - # Avoid unnecessary edges. - haskey(aliases, other_var) || continue - # Edge from dependency to dependent - Graphs.add_edge!(dep_graph, other_var, var) + for (ivar, i) in cand_row + for other in A.row_cols[i] + other == ivar && continue + haskey(cand_row, other) || continue + Graphs.add_edge!(dep_graph, other, ivar) end end - - # We know there won't be cycles because we're using the results of tearing, which - # partitioned this SCC into a lower var_order = Graphs.topological_sort(dep_graph) - for var in var_order - haskey(aliases, var) || continue - iszero(Graphs.indegree(dep_graph, var)) && continue - coeffs = aliases[var] - cst = constants[var] - I, V = SparseArrays.findnz(coeffs) + # Pass 3: greedy commit in topological order, composing references to already + # eliminated variables inline so committed aliases are always closed over + # retained variables. A variable whose composed expression exceeds the node + # budget stays in the numeric core — this bounds symbolic growth (the failure + # mode reported in StateSelection.jl#95) while still eliminating the cheap + # sequential rows. + elim_budget = something(tryparse(Int, get(ENV, "MTKTEARING_ELIM_SIZE_BUDGET", "")), 1000) + for ivar in var_order + i = get(cand_row, ivar, nothing) + i === nothing && continue + var_coeff = cand_piv[ivar] new_I = Int[] new_V = Num[] - sizehint!(new_I, length(I)) - sizehint!(new_V, length(V)) - for (other_var, coeff) in zip(I, V) - other_coeffs = get(aliases, other_var, nothing) - if other_coeffs === nothing - push!(new_I, other_var) - push!(new_V, coeff) - continue + cst = b[i] / var_coeff + for (other, coeff) in zip(A.row_cols[i], A.row_vals[i]) + other == ivar && continue + c = coeff / (-var_coeff) + other_alias = get(aliases, other, nothing) + if other_alias === nothing + push!(new_I, other) + push!(new_V, c) + else + oI, oV = SparseArrays.findnz(other_alias) + append!(new_I, oI) + append!(new_V, Iterators.map(Base.Fix2(*, c), oV)) + cst += c * constants[other] + end + end + sv = SparseArrays.sparsevec(new_I, new_V, N) + svI, svV = SparseArrays.findnz(sv) + keep = findall(!StateSelection.CLIL.cheap_iszero, svV) + svI = svI[keep]; svV = svV[keep] + # Budget check over the composed coefficients and constant (shared DAG + # nodes counted once via the common `seen` set). + seen = IdDict{Any, Nothing}() + sz = _node_count_capped(cst, elim_budget, seen) + if sz <= elim_budget + for v in svV + sz += _node_count_capped(v, elim_budget, seen) + sz > elim_budget && break end - other_coeffs = other_coeffs::valtype(aliases) - other_cst = constants[other_var] - other_I, other_V = SparseArrays.findnz(other_coeffs) - append!(new_I, other_I) - append!(new_V, Iterators.map(Base.Fix2(*, coeff), other_V)) - cst += coeff * other_cst end + sz > elim_budget && continue # too big: keep this variable in the core - # `sparsevec` sums duplicate indices but keeps explicit zeros; drop them. - aliases[var] = SparseArrays.dropzeros!(SparseArrays.sparsevec(new_I, new_V, length(coeffs))) - constants[var] = cst + eqs_mask[i] = false + vars_mask[ivar] = false + new_N -= 1 + aliases[ivar] = SparseArrays.sparsevec(svI, svV, N) + constants[ivar] = cst end # First we update `b`, since doing so requires the unmodified `A`. @@ -817,6 +1368,23 @@ function __reduce_linear_system!(A::StateSelection.CLIL.SparseMatrixCLIL{Num, In old_to_new_eq[.!eqs_mask] .= 0 A = StateSelection.get_new_mm(aliases, old_to_new_eq, old_to_new_var, A) + if _check + # Best-effort diagnostics: a probe point can still violate an expression's + # domain (e.g. sqrt of a negative subexpression); skip rather than abort. + try + A_red_check = collect(A)::Matrix{Num} + _reduction_identity_ok(A0_check, b0_check, A_red_check, b, aliases, constants, + eqs_mask, vars_mask, old_to_new_eq) + _reduction_rank_report(A0_check, b0_check, A_red_check, b, alg_vars) + catch err + err isa InterruptException && rethrow() + # The self-check is best-effort, opt-in diagnostics; a probe point can + # violate an expression's domain or an exotic term can resist numeric + # evaluation. Skip rather than abort compilation. + @warn "Inline-linear-SCC self-check skipped" block_n=length(b0_check) exception=err + end + end + return A, b, eqs_mask, vars_mask end diff --git a/lib/ModelingToolkitTearing/src/stateselection_interface.jl b/lib/ModelingToolkitTearing/src/stateselection_interface.jl index 7995305..9eb5e39 100644 --- a/lib/ModelingToolkitTearing/src/stateselection_interface.jl +++ b/lib/ModelingToolkitTearing/src/stateselection_interface.jl @@ -105,6 +105,13 @@ function StateSelection.linear_subsys_adjmat!(state::TearingState; kwargs...) return mm end +# Structural zero check for symbolic CLIL values: `Base.iszero(::Num)` performs +# a semantic (expansion-based) zero test that can OOM on large coefficient +# expressions (e.g. multibody models), while explicit stored zeros produced by +# duplicate-index summation are always structural `Const(0)`. +StateSelection.CLIL.cheap_iszero(x::Num) = SU._iszero(Symbolics.unwrap(x)) +StateSelection.CLIL.cheap_iszero(x::SymbolicT) = SU._iszero(x) + function maybe_zeros_descend(ex::SymbolicT) @match ex begin BSImpl.AddMul(; variant) => return variant === SU.AddMulVariant.MUL diff --git a/lib/ModelingToolkitTearing/test/runtests.jl b/lib/ModelingToolkitTearing/test/runtests.jl index 8b983ac..0bae3c3 100644 --- a/lib/ModelingToolkitTearing/test/runtests.jl +++ b/lib/ModelingToolkitTearing/test/runtests.jl @@ -11,6 +11,7 @@ import SymbolicUtils as SU using SymbolicUtils: unwrap using Setfield using ForwardDiff +import SparseArrays @testset "`InferredDiscrete` validation" begin k = ShiftIndex() @@ -171,6 +172,124 @@ end ) reassemble_alg = MTKTearing.DefaultReassembleAlgorithm(; inline_linear_sccs = true) end +@testset "`__reduce_linear_system!` preserves the full-system residual" begin + SymT = Symbolics.SymbolicT + MVT = StateSelection.MatchedVarT + + # Build the 4×4 SCC described in issue #98's plan. Variables x1..x4, equations e1..e4: + # e1: 2*x2 = 0 # eliminates x2 (matched, const coeffs) + # e2: -3*x2 + x3 = 0 # eliminates x3 via x2 -> transitive chain + # e3: x1 + x2 + x3 + x4 = p # RETAINED, references two eliminated vars, symbolic RHS + # e4: x1 + x4 = p # RETAINED, makes the reduced block rank-deficient + # Matching: x2->e1, x3->e2 (eliminated); x1,x4 unassigned => e3,e4 retained. + @variables p + mkA() = StateSelection.CLIL.SparseMatrixCLIL{Num, Int}( + 4, 4, collect(1:4), + [[2], [2, 3], [1, 2, 3, 4], [1, 4]], + [Num[2.0], Num[-3.0, 1.0], Num[1.0, 1.0, 1.0, 1.0], Num[1.0, 1.0]]) + mkb() = SymT[unwrap(Num(0)), unwrap(Num(0)), unwrap(p), unwrap(p)] + vem = BipartiteGraphs.complete( + BipartiteGraphs.Matching{MVT}(Union{MVT, Int}[BipartiteGraphs.unassigned, 1, 2, BipartiteGraphs.unassigned]), + 4) + + Ar, br, em, vm = MTKTearing.__reduce_linear_system!(mkA(), mkb(), vem, collect(1:4), collect(1:4)) + + @test em == Bool[0, 0, 1, 1] + @test vm == Bool[1, 0, 0, 1] + + # The reduction is exact: x2=0, x3=0, so both retained rows become `x1 + x4 = p`. + subs = Dict{Any, Float64}(unwrap(p) => 3.7) + ev(x) = MTKTearing._evalnum(x, subs) + @test ev.(collect(Ar)) == [1.0 1.0; 1.0 1.0] # rank-deficient (rank 1), as expected + @test ev.(br) ≈ [3.7, 3.7] # consistent: b in range(A) + + # Exercise the opt-in self-check code path end-to-end (snapshot + identity + rank report). + local res + withenv("MTKTEARING_CHECK_REDUCTION" => "1") do + res = MTKTearing.__reduce_linear_system!(mkA(), mkb(), vem, collect(1:4), collect(1:4)) + end + @test res[3] == Bool[0, 0, 1, 1] +end + +@testset "`_reduction_identity_ok` detects reduction errors" begin + SymT2 = Symbolics.SymbolicT + @variables p + # Full 2×2 system: e1: 2*x1 = p (eliminate x1), e2: x1 + x2 = 0 (retain x2). + # Correct reduction: x1 = p/2, so e2 becomes x2 = -p/2. + A0 = Num[2.0 0.0; 1.0 1.0] + b0 = SymT2[unwrap(p), unwrap(Num(0))] + aliases = Dict{Int, SparseArrays.SparseVector{Num, Int}}(1 => SparseArrays.spzeros(Num, 2)) + constants = Dict{Int, SymT2}(1 => unwrap(p / 2)) + eqs_mask = BitVector([false, true]) + vars_mask = BitVector([false, true]) + old_to_new_eq = [0, 1] + + A_red = Num[1.0;;] + b_red_good = SymT2[unwrap(-p / 2)] + @test MTKTearing._reduction_identity_ok( + A0, b0, A_red, b_red_good, aliases, constants, eqs_mask, vars_mask, old_to_new_eq) + + # A wrong RHS (off by a constant) must be caught. + b_red_bad = SymT2[unwrap(-p / 2 + 1)] + bad = @test_logs (:warn,) match_mode = :any MTKTearing._reduction_identity_ok( + A0, b0, A_red, b_red_bad, aliases, constants, eqs_mask, vars_mask, old_to_new_eq) + @test bad == false +end + +@testset "`_group_inline_linear_families` merges coupled rank-deficible blocks" begin + @variables a(t) b(t) c(t) d(t) + @parameters p + # Equations indexed 1..4. Blocks 1,2,4 reference the `maybe_zeros` parameter `p` + # (their rank can drop); block 3 does not. + neweqs = [ + p * a ~ 0, # eq1 + 0 ~ b - p * a, # eq2 (couples to block 1's variable below) + 0 ~ c - a, # eq3 (no `p`) + 0 ~ d - p * c, # eq4 + ] + prepared = NTuple{2, Vector{Int}}[([1], [1]), ([2], [2]), ([3], [3]), ([4], [4])] + g = BipartiteGraph(4, 4) + add_edge!(g, BipartiteEdge(2, 1)) # block 2's equation references block 1's variable + # (block 4 is intentionally NOT coupled to block 3) + mz = Symbolics.SymbolicT[unwrap(p)] + + # Coupled + rank-deficible blocks 1,2 merge; block 3 (not deficible) and block 4 + # (not coupled to blocks 1,2) stay separate. + @test MTKTearing._group_inline_linear_families(prepared, mz, neweqs, g, false, true) == + [[1, 2], [3], [4]] + + # No `maybe_zeros` => no grouping at all (default one-block-per-SCC behaviour). + @test MTKTearing._group_inline_linear_families(prepared, Symbolics.SymbolicT[], neweqs, g, false, true) == + [[1], [2], [3], [4]] + # Discrete or inline-linear disabled => singletons. + @test MTKTearing._group_inline_linear_families(prepared, mz, neweqs, g, true, true) == + [[1], [2], [3], [4]] + @test MTKTearing._group_inline_linear_families(prepared, mz, neweqs, g, false, false) == + [[1], [2], [3], [4]] + + # Non-adjacent family members are connected through the dependency DAG. Blocks 1 and 3 + # are rank-deficible; block 2 is not, but lies on the dependency path 1 → 2 → 3, so the + # path closure pulls it into the family. + neweqs2 = [ + p * a ~ 0, # block 1 (deficible) + 0 ~ b - a, # block 2 (not deficible, on the path between 1 and 3) + 0 ~ c - p * b, # block 3 (deficible) + ] + prepared2 = NTuple{2, Vector{Int}}[([1], [1]), ([2], [2]), ([3], [3])] + g2 = BipartiteGraph(3, 3) + add_edge!(g2, BipartiteEdge(2, 1)) # block 2's equation references block 1's variable + add_edge!(g2, BipartiteEdge(3, 2)) # block 3's equation references block 2's variable + @test MTKTearing._group_inline_linear_families(prepared2, mz, neweqs2, g2, false, true) == + [[1, 2, 3]] + + # An unrelated block between two family members is NOT pulled in (no path through it), + # and the family is emitted before it (contracted topological order, stable by index). + g3 = BipartiteGraph(3, 3) + add_edge!(g3, BipartiteEdge(3, 1)) # block 3's equation references block 1's variable + @test MTKTearing._group_inline_linear_families(prepared2, mz, neweqs2, g3, false, true) == + [[1, 3], [2]] +end + @testset "`system_subset(::SystemStructure)` subsets `.state_priorities`" begin @variables x(t) y(t) [state_priority = 2] z(t) [state_priority = 5] @named sys = System([D(x) ~ x, D(y) ~ y, D(z) ~ z], t) diff --git a/src/math/sparsematrixclil.jl b/src/math/sparsematrixclil.jl index ad5a5a1..2d9b15e 100644 --- a/src/math/sparsematrixclil.jl +++ b/src/math/sparsematrixclil.jl @@ -91,6 +91,18 @@ zero!(a::SparseVector) = (empty!(a.nzind); empty!(a.nzval)) zero!(a::CLILVector) = zero!(a.vec) SparseArrays.dropzeros!(a::CLILVector) = SparseArrays.dropzeros!(a.vec) +""" + cheap_iszero(x) + +Structural zero check used by [`SparseArrays.dropzeros!`](@ref) on +`SparseMatrixCLIL`. Defaults to `Base.iszero`. Downstream packages whose CLIL +value type is symbolic should overload this with a cheap *structural* check: +`Base.iszero` on e.g. `Symbolics.Num` performs a semantic (expansion-based) +zero test that can be arbitrarily expensive on large expressions, while +explicit stored zeros are always structural zeros. +""" +cheap_iszero(x) = iszero(x) + # Remove explicitly-stored zeros from each row, in place. function SparseArrays.dropzeros!(S::SparseMatrixCLIL) for r in eachindex(S.row_vals) @@ -98,7 +110,7 @@ function SparseArrays.dropzeros!(S::SparseMatrixCLIL) vals = S.row_vals[r] j = 0 for k in eachindex(vals) - iszero(vals[k]) && continue + cheap_iszero(vals[k]) && continue j += 1 cols[j] = cols[k] vals[j] = vals[k] 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 diff --git a/src/utils.jl b/src/utils.jl index 5960c36..b0dfbee 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -345,7 +345,10 @@ function get_new_mm( # entry: a prior cancellation may have `pop!`ed the matching entry. if !isempty(final_row_cols) && col == final_row_cols[end] final_row_vals[end] += new_row_val_i[indices[i]] - if iszero(final_row_vals[end]) + # Syntactic zero test: semantic `iszero` on symbolic coefficients + # can OOM via polynomial expansion (#95); an uncancelled exact zero + # is pruned later by CLIL `dropzeros!`. + if CLIL.cheap_iszero(final_row_vals[end]) pop!(final_row_cols) pop!(final_row_vals) end 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..ed58396 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,17 @@ using Test include("bareiss.jl") include("carpanzano_tearing.jl") +include("dummy_derivative_blocks.jl") + +# A value type whose `Base.iszero` is "semantic" and must never be consulted by +# `dropzeros!` — stands in for symbolic value types (e.g. `Symbolics.Num`), +# where the semantic zero test is arbitrarily expensive. `dropzeros!` must go +# through the `cheap_iszero` hook instead. +struct SemanticZero + x::Int +end +Base.iszero(::SemanticZero) = error("semantic iszero must not be called by dropzeros!") +SSel.CLIL.cheap_iszero(v::SemanticZero) = v.x == 0 @testset "`get_new_mm`" begin mm = SSel.CLIL.SparseMatrixCLIL( @@ -72,3 +83,11 @@ include("carpanzano_tearing.jl") @test mm2.row_cols == [[2]] @test mm2.row_vals == [[1]] end + +@testset "`dropzeros!` uses the `cheap_iszero` hook" begin + mm = SSel.CLIL.SparseMatrixCLIL{SemanticZero, Int}( + 1, 2, [1], [[1, 2]], [[SemanticZero(0), SemanticZero(3)]]) + SparseArrays.dropzeros!(mm) + @test mm.row_cols == [[2]] + @test map(v -> v.x, only(mm.row_vals)) == [3] +end