feat: jointly solve coupled rank-deficible inline-linear SCC families (#98)

[baggepinnen]

The family-grouping half of the original PR, now stacked on #106 (the get_linear_scc_linsol correctness fix). The base of this PR is fbc/inline-linear-scc-buried-b-repair, so the diff below is only the grouping work; it will auto-retarget to main once #106 merges.

Mechanism (the actual #98 cause)

The diagnostics added in #106 disproved the construction-bug theory: every emitted block is exact at generic parameter values and the repair pass never fires on HalfCar. The real #98 mechanism is cross-block indeterminacy at the degenerate parameter point — blocks are individually exact but solved sequentially, so an upstream rank-deficient block's gauge choice (min-norm, or garbage from a plain LU) is substituted downstream and makes a dependent block inconsistent, even though the union of the family's equations is satisfiable.

Fix

Group coupled, runtime-rank-deficible inline-linear SCCs into families along the block-level dependency DAG and solve each family as one joint linear system, so the gauge is resolved consistently across the whole family. When maybe_zeros is empty the grouping is inert and all existing behaviour is unchanged. Later commits refine this with set-dependent linearity + greedy convex growth, peel-on-failure for nonlinear members, and an equation-count budget on family closures.

Known dependency / why it's separated

The merged family is still emitted as INLINE_LINEAR_SCC_OP(A, b); a legitimately rank-deficient-but-consistent merged block needs a rank-tolerant runtime solve (the sparse direction of #95) to actually be solved. That rank-tolerant solve is closer to the true root cause than the grouping is — see the analysis in #105. Keeping this separate from #106 lets the correctness fix land independently while the grouping is evaluated against that direction.

[baggepinnen]

Ran the pending end-to-end confirmation on HalfCar (this branch + #99's cheap_iszero, MTKTEARING_CHECK_REDUCTION=1, with the rank-revealing runtime-solve override from the #97/#98 measurements). Three findings:

1. Self-check needed a fix (pushed, cec2073)

Symbolics.substitute defaults to fold = Val(false), so fully-numeric expressions like -0.63*sin(-0.54) stayed symbolic and _evalnum's Float64() threw on every block with trig coefficients. With fold = Val(true) the check runs cleanly on the full model.

2. The check is valuable — it disproves the construction-bug theory

Every emitted block, including the big 396 → 295 reduction, passes at generic (random) symbol values:

full_n = 396  full_rank = 396  full_consistent = true  full_relres = 7.9e-14
reduced_n = 295  reduced_rank = 295  reduced_consistent = true
# all 13/15/12 blocks: full_rank == reduced_rank, relres ≤ 2e-15

The repair pass never fires on this model (0 re-expansions) — there is no buried SCC variable in b. Block construction is exact; the hypothesized has_edge desync does not occur here.

3. The real #98 mechanism: cross-block indeterminacy at the degenerate parameter point

At the true parameter/state values (where the joint-axis components are 0 and rank drops), the per-corner block chains behave differently:

corner 1:  n=13 rank=13           → unique solve, |x|=1.25
           n=15 rank=14 relres=2e-16   → rank-deficient but CONSISTENT ✓
corner 2:  n=13 rank=12, b=0      → 1-dim nullspace, min-norm picks x=0 (a gauge choice)
           n=15 rank=14 relres=0.945   → INCONSISTENT ✗

The blocks are individually exact as symbolic functions, but they are solved sequentially, and at the degenerate point the static-reaction indeterminacy spans blocks: the upstream rank-deficient block's solution choice (min-norm here; garbage LU in production) fixes a gauge that no solution of the downstream block is compatible with. The union of the corner's equations is satisfiable (the state satisfies the model to 4e-10; the full 396-block is consistent), so splitting into sequentially-solved sub-blocks is what invalidates the system — not how each block is built.

Implications

• This PR's repair pass is sound hardening but does not fix Inline linear SCC emits an inconsistent runtime system (rank 14/15, ‖Ax−b‖/‖b‖ ≈ 1) at a consistent state — integration goes Unstable #98; the diagnostics (with cec2073) are worth keeping.
• The actual fix needs one of: (a) solve runtime-rank-deficient block families jointly (the sparse joint-solve direction of Inline linear SCCs: 296×296 dense runtime solve for an ~18-unknown core; non-rank-tolerant; symbolic shrinkage intractable — solve sparse+rank-tolerant instead #95 gets this for free), or (b) use maybe_zeros information at emission time — the axis parameters whose zeros cause the rank drop are exactly the ones markable as maybe_zeros, so the pass could refuse to split SCC families whose runtime rank can drop.
• HalfCar integration is unchanged on this branch (Unstable), as expected given the above.

I'll update #98 with the corrected mechanism.

[baggepinnen]

Pushed the joint-solve approach (2e16a5f), per your finding that #98 is cross-block indeterminacy rather than a construction bug.

What it does: in generate_system_equations!, coupled runtime-rank-deficible inline-linear SCCs are grouped into families and each family is solved as one joint linear system, so the upstream gauge choice is no longer frozen and substituted into a dependent block.

A family is a maximal run of consecutive blocks that are:

• rank-deficible — their equations reference a maybe_zeros parameter (so a coefficient can vanish and drop rank), and
• coupled — each block's equations structurally reference the previous block's variables.

This should merge each corner's 13+15(+12) chain while leaving the full-rank 295 chassis block alone (it isn't rank-deficible). When maybe_zeros is empty the grouping is completely inert, so default behaviour and all existing tests are unchanged.

Mechanics: the reassembly loop now emits one block at a time via an emit_block! helper; a merged family passes the union of its equations/variables to a single get_linear_scc_linsol, with a fallback to per-member emission if the joint inline solve doesn't apply. Added a unit test for the grouping decision; full ModelingToolkitTearing + StateSelection suites pass.

Still needed on your side: the joint family is emitted as INLINE_LINEAR_SCC_OP(A, b) = \, and the merged block is itself legitimately rank-deficient-but-consistent at the degenerate point — so it still needs the rank-tolerant runtime solve (your #97/#98 override / the sparse direction of #95) to not throw and to pick a consistent min-norm solution. With that override active, please re-run HalfCar with MTKTEARING_CHECK_REDUCTION=1:

• the per-corner blocks should now appear as a single larger joint block in the diagnostics,
• reduced_consistent = true / relres ≈ 0 for it, and
• integration should no longer go Unstable.

A couple of points I'd value your read on:

1. The grouping heuristic ("consecutive + coupled + each references a maybe_zeros param") is conservative. If the real coupling on HalfCar skips a block or spans non-adjacent SCCs, the run-based grouping may need to follow the dependency DAG instead of adjacency — easy to extend if the diagnostics show a family being split.
2. If you'd rather pursue option (b) from your comment (refuse to split such families during tearing) instead of merging at emission time, say so and I'll move it upstream.

[baggepinnen]

Re-ran HalfCar against 2e16a5f (+ #99, rank-tolerant runtime override, MTKTEARING_CHECK_REDUCTION=1): the grouping is inert on this model — runtime blocks are unchanged (13/15/12 per corner, the 15×15 still relres = 0.945, integration still Unstable).

Instrumented _group_inline_linear_families to see why:

FAMGRP n=681 sizes=[1, 1, ..., 13, ...(26 singletons)..., 13, ...(45)..., 15, ...(21)..., 12, ..., 396, ...]
FAMGRP droppable=[159, 160, 161, 165, ..., 655, 671, 673, 675]   # 60 blocks, incl. all the big ones
FAMGRP coupled=Int64[]                                            # empty!

So, answering your question 1 directly: adjacency-run grouping cannot work here. The prepared block list has 681 entries, almost all singletons, and the members of each corner family are separated by 20–45 interleaved singleton blocks in the topological order. No two droppable blocks are ever adjacent, hence coupled is identically empty and every group is a singleton.

What the grouping needs instead

Follow the block dependency DAG:

1. Build the block-level DAG (edge j → k iff block k's equations reference block j's variables — the same has_edge test you already use, just not restricted to k-1).
2. A family = a connected component of droppable blocks under reachability through that DAG.
3. The merged system must be closed: include the (mostly size-1, cheap) non-droppable blocks lying on dependency paths between family members — otherwise the merged block's equations reference variables that are neither solved upstream nor part of the block. Interval closure in topological order is a simple way to get this, with your existing per-member fallback when the joint solve doesn't apply.

One caveat to decide on: the 396 chassis block is itself droppable (it references the axis parameters too) and is downstream of the corner families, so a pure reachability closure will pull it in and produce one ~500-unknown joint block. That is correct (and the measured 295 reduction was consistent at the true point, so it is not necessary) — if the size is a concern, the grouping could stop at blocks whose own rank was never observed to drop, but there is no static information to distinguish "deficible and actually drops" from "deficible but stays full rank", so I'd accept the big merge and lean on the #95 sparse/rank-tolerant runtime solve for cost.

On question 2: emission-time merging is the right place, IMO. Refusing to split during tearing would forfeit the block-triangular structure for all operating points to guard a degenerate one, and tearing doesn't know the runtime solve strategy; generate_system_equations! has exactly the information needed (maybe_zeros, the block list, the graph).

Side note from this run: both the 5× and the warm-started 30× init LM solves converged this time (resid 1.6e-9 / 4.9e-9) — the init story is in place; the joint-family emission is the last blocker for HalfCar integration.

[baggepinnen]

Pushed the DAG-based grouping (e77da0d), following your instrumented findings exactly:

• Block-level dependency DAG: edge j → k iff block k's equations structurally reference block j's variables — built from the bipartite incidence over all prepared blocks, not restricted to k-1.
• Family = connected component of rank-deficible blocks under reachability through the DAG (intermediate blocks allowed).
• Path closure: each family absorbs the blocks lying on dependency paths between its members, so the merged system is self-contained. (I went with exact path closure rather than interval closure — with 20–45 unrelated singletons interleaved, the interval would drag in blocks that may be nonlinear in the union's variables, and one such block makes get_linear_scc_linsol return nothing for the whole family → fallback → inert again.)
• Emission order: topological order of the DAG with each family contracted to one node (stable Kahn, preferring smallest original block index). Families are path-closed, hence convex, so contraction can't create cycles; a defensive check falls back to per-SCC emission if coverage ever fails.

On the chassis caveat: agreed — a pure reachability component will pull the droppable 396 block into the corner families' component (they all reach it), producing one big joint block. I left that as-is per your "accept the big merge" call; the cost story is #95's sparse/rank-tolerant solve.

One thing to watch in your instrumentation: with the corner families now merged into a single large block, the analytical_linear_scc_limit path is irrelevant (way over the limit) and everything rides the runtime INLINE_LINEAR_SCC_OP — so the rank-tolerant override remains load-bearing for both correctness (gauge) and SingularException avoidance.

Same FAMGRP-style probe should now show a non-trivial family; if the closure still splits something (e.g. a coupling that only appears post-total_sub substitution rather than in the structural graph), paste the trace and I'll chase it.

Tests: 6/6 on the grouping unit tests (incl. non-adjacent members merged through a path, and an unrelated in-between block left out), full ModelingToolkitTearing + StateSelection suites pass; grouping remains inert when maybe_zeros is empty.

[baggepinnen]

Tested e77da0d on HalfCar (same setup). The DAG grouping does form a family now — but it over-collects, the joint solve fails on a nonlinear member, and the fallback silently restores the old behavior (runtime blocks unchanged: 15×15 still relres = 0.945, integration still Unstable).

The trace you asked for:

FAMGRP groups=245 families=1 famsizes=[980]        # one family of 980 equations
LINSOL-FAIL N=980 reason=nonlinear
  var = (excited_suspension_l₊suspension₊r123₊jointUSR₊rod1₊frame_a₊R(t))[3, 3]
  eq#59: 0 ~ ifelse(point_gravity == true, ((-mu) / ((mass₊body₊r_0(t))[1] + mass₊body₊r_cm[1]*R[1,1] + ...)^...), ...)

Two compounding problems:

1. The reachability component over-collects. All droppable blocks reach the (droppable) 396 chassis block, so the family becomes "essentially everything", including kinematic blocks whose variables are rotation-matrix entries. The Inline linear SCC emits an inconsistent runtime system (rank 14/15, ‖Ax−b‖/‖b‖ ≈ 1) at a consistent state — integration goes Unstable #98 gauge coupling lives purely in the reaction-force network, which is linear in the union's force variables — but once R-entries are union variables, any equation nonlinear in them poisons the family.
2. A single nonlinear member kills the whole family. Here it's a gravity equation: the point_gravity branch divides by powers of position, which is nonlinear in the R-entries. (Side note for us, not for this PR: point_gravity being a runtime-ifelse parameter rather than a structural switch is what drags this nonlinearity into an otherwise-linear block — worth a look in the component library.)

Suggested refinement: peel-on-failure instead of all-or-nothing fallback

The failure site already identifies the offending equation and variable. Rather than falling back to per-member emission for the whole family:

1. attempt the joint solve for the family;
2. on islinear == false, map the offending equation back to its member block, remove that block from the family (re-running path closure on the survivors; members that become disconnected revert to singletons);
3. retry. Iterate until the joint solve succeeds or the family collapses to singletons.

This converges (each iteration removes ≥1 block), needs no new analysis machinery, and automatically prunes the kinematic/gravity blocks while keeping the force-network core — which is the part that actually exchanges gauge. The per-attempt cost is the existing symbolic b-build loop, compile-time only.

An equivalent pre-filter (test each candidate block's equations for linearity in the union variables before forming the family) would avoid the retries, but it duplicates what get_linear_scc_linsol already computes, so peeling seems the cheaper patch.

[baggepinnen]

Pushed the grouping rework in 590d510, after a long instrumented session on HalfCar. The #98 inconsistency is now fixed at runtime. Headline measurement (rank-revealing runtime solve, true parameter point):

before:  n=13 rank=12 relres=0.0   |  n=15 rank=14 relres=0.945  ← INCONSISTENT
after:   n=131 rank=129 relres ≤ 1.8e-10                          ← rank-deficient (2 gauge dims) but CONSISTENT
         n=12 / n=12 / n=301: all full-rank, relres ≤ 1e-14

The corners' gauge-coupled reaction chains are merged into one joint block; min-norm resolves both corners' static indeterminacy with a single consistent gauge.

Getting there required four corrections to e77da0d, each driven by an instrumented failure:

1. The dependency "DAG" wasn't acyclic. Raw residual incidence contains backward references (through torn variables of later blocks), so contraction always stalled Kahn and the coverage fallback silently degraded everything to singletons — this is why e77da0d's families never reached emission. Prepared blocks are in BLT order, a valid linearization of the matching-based dependencies, so the graph now keeps forward edges only.
2. Linearity is set-dependent, not global. A first attempt classified blocks as globally "mergeable" (linear in all incident prepared variables) — which disqualified exactly the rank-deficient force blocks the family must absorb (they're nonlinear in other blocks' angles via cos(phi), which only matters if those blocks join). Mergeability is now a memoized pairwise check validated over the candidate family's closure.
3. Convexity is a correctness requirement, peeling is not viable. Emission-time peeling (9f12f16) produced real evaluation cycles: a peeled interior block both consumes family outputs and feeds family inputs, so topsort_equations rightly rejected the result. Families now grow greedily in BLT order while the full-DAG closure stays pairwise-linear, so they are convex by construction; the joint-solve fallback is a plain per-member emission.
4. Even disjoint convex families can be mutually unorderable (interleaved chains). Kahn now dissolves the latest stalled family and retries, keeping a maximal orderable subset.

Two caveats / follow-ups:

• Integration is still Unstable on HalfCar — but this is now demonstrably not a linear-algebra problem (all blocks consistent, clean rank gap 0.11 → 1e-16, FD Jacobian finite and step-size independent). It's the separate integrator-level mystery (dt collapses to 1e-21 from a fully consistent state with a smooth RHS; FBDF reports a singular 22×22 W) that I'll pursue independently — possibly related to the 2 gauge dimensions now living inside the runtime block.
• The rank-deficient-but-consistent 131-block requires the rank-tolerant runtime solve (min-norm pivoted QR; stock \ would throw). That's DyadCompilerPasses#41 / the Inline linear SCCs: 296×296 dense runtime solve for an ~18-unknown core; non-rank-tolerant; symbolic shrinkage intractable — solve sparse+rank-tolerant instead #95 sparse direction — load-bearing for this PR's benefit to materialize in production.

Suites: ModelingToolkitTearing + StateSelection pass locally; the grouping remains inert when maybe_zeros is empty. Note 590d510 also relies on #99's cheap_iszero to avoid the dropzeros OOM on this model — worth merging #99 first.

Side fix landed in MultibodyComponents (04f7359): World.point_gravity is now a structural parameter of the consuming components (Body etc.), so the unused gravity branch — previously a runtime ifelse nonlinear in positions/orientations — never enters the equations. That removed the gravity blocks from the non-mergeable set here, and surfaced two more MTK bugs along the way (SciML/ModelingToolkit.jl#4615, #4616).

[baggepinnen]

End-to-end confirmation on HalfCar (the pending item in the description), run with this branch merged together with #99 into fix-state-priority-merged-blocks, and ModelingToolkitTearing dev'd from lib/:

• MTKTEARING_CHECK_REDUCTION=1 over the HalfCar compile: every per-block report consistent (reduced_relres = 0.0), full = reduced rank throughout. The small blocks report full rank at the compile-time probe points — which is exactly the parameter-degeneracy signature: at the real axis-aligned parameter values two of them are rank-deficient (details and the physical identification posted on Inline linear SCC emits an inconsistent runtime system (rank 14/15, ‖Ax−b‖/‖b‖ ≈ 1) at a consistent state — integration goes Unstable #98).
• Runtime probe at the real init point, pre-model-fix: the genuinely deficient blocks remained (n=13 rank 12 relres 0; n=15 rank 14 relres ≈ 0.93 — upstream-gauge contamination). No false "inconsistent construction" was observed; what was left was model-level, as the 2e16a5f analysis concluded.
• After the model-level state-priority fix (JuliaComputing/MultibodyComponents.jl#67, 65febda): every emitted block full rank, relres ≈ 1e-14 at the real point, no SingularException, no Unstable — HalfCar (0,2) s and FullCar (0,1) s integrate to Success with Rodas5P(autodiff=AutoFiniteDiff()) + BrownFullBasicInit.

[AayushSabharwal]

I'm convinced that get_linear_scc_linsol can be wrong. I'm not yet convinced of this grouping. Is this PR trying to do multiple things at once? Can it be split up?

[baggepinnen]

You're right on both counts, and yes — it's two separable things in two different functions, I'll split them.

PR 1 — the get_linear_scc_linsol fix. This is the part you're convinced of, and it's self-contained: the buried-b repair pass + the safe non-inlined fallback + the synthetic tests + the opt-in self-check. It touches only get_linear_scc_linsol/__reduce_linear_system!. Worth noting this is a latent correctness bug on its own — but it is not what fixes HalfCar. In fact the self-check it adds is what disproved my original construction-bug theory: the repair pass never fires on HalfCar, every block is exact at generic parameters. So I'd like to land this purely as correctness + diagnostics.

PR 2 — the family grouping. This is the actual #98 mechanism (sequential per-block solves freeze an upstream rank-deficient gauge that then makes a downstream block inconsistent, even though the family's union is satisfiable). It reworks the emission loop and is the larger, more speculative change — and it's incomplete: the merged block still needs a rank-tolerant runtime solve (the #95 direction) to actually be solved. Given that dependency, it probably belongs in the same conversation as #105 rather than here. I'll pull it into its own PR so we can evaluate it against that direction without blocking the correctness fix.

[baggepinnen]

Split done:

• Inline-linear-SCC: assert b-free invariant + opt-in reduction self-check (#98) #106 — fix: ensure inline-linear-SCC b is free of SCC variables — the get_linear_scc_linsol correctness fix you're convinced of, off main, self-contained (repair pass + safe fallback + synthetic tests + opt-in self-check).
• This PR (feat: jointly solve coupled rank-deficible inline-linear SCC families (#98) #100) — now scoped to just the family grouping, rebased on top of Inline-linear-SCC: assert b-free invariant + opt-in reduction self-check (#98) #106 (base = fbc/inline-linear-scc-buried-b-repair). The diff here is grouping-only and will retarget to main automatically when Inline-linear-SCC: assert b-free invariant + opt-in reduction self-check (#98) #106 merges.

So the correctness fix can land on its own, and the grouping stays open for evaluation against the rank-tolerant-runtime-solve direction (#95 / #105) without blocking it.