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

[baggepinnen]

Problem

Fixes #101. Fixes #102.

dummy_derivative_graph! consults state_priority only within each SCC of
find_var_sccs(graph, var_eq_matching), but the SCC partition of the (priority-blind)
Pantelides matching can be strictly finer than the blocks of the Mattsson–Söderlind
selection subproblem (differentiated equations × highest-derivative candidates). When a
differentiated equation matched in one SCC is incident to candidate variables in another
— e.g. a twice-differentiated connection alias 0 ~ D²(x) - D²(y) matched to D²(x)
while D²(y) belongs to a kinematic-loop SCC — the prioritized variable sits in a
singleton SCC where demotion is unconditional. state_priority is then silently
ineffective, and the resulting state realization can be singular.

Concretely (#101): in a planar 2-link arm with joint coordinates at state_priority = 10
(or 100 — the value cannot matter) and body Cartesian variables at priority 1–2, the
Cartesian variables were selected as states together with an algebraic equation for a
joint angle whose Jacobian cos(phi) is singular exactly at the example's target
configuration; solvers abort with dt → ε as tracking improves.

Fix

1. merge_dummy_derivative_blocks: after find_var_sccs, merge (union-find) SCCs
   that are coupled through differentiated equations incident to dummy-derivative
   candidates in other SCCs, using the same candidate filter as the selection loop.
   The existing priority-sorted greedy demotion with augmenting-path feasibility then
   runs per merged block — on a merged block that greedy selects a minimum-priority
   basis of the transversal matroid, i.e. it provably keeps the highest-priority states
   subject to structural validity. When nothing merges, the input is returned unchanged
   (===), so the common path is unaffected.
2. Jacobian column permutation (dummy_derivative_graph!: Jacobian columns not permuted after priority sort — vars[col_order[i]] indexes mismatched orders #102): when the candidates are priority-sorted, the
   previously-computed integer Jacobian is now column-permuted alongside (J = J[:, var_perm]); bareiss' col_order indexes into the sorted candidate list, so the two
   orders must agree. This latent mismatch becomes reachable in practice once blocks can
   merge (a merged block may have an all-integer Jacobian and non-uniform priorities).

Validation

On the motivating model (TwoJointPlanarRobotPID, MultibodyComponents.jl twodof
branch; joint priorities 10, body priorities 1–2):

	before	after
forward-arm states	body.phi/w, body.r[2]/v[2]	revolute1.phi/w, revolute2.phi/w
algebraic equations	1 (singular at target)	0
Rodas5P	Unstable at t ≈ 2.2 (always)	solves, 72 steps
FBDF	solves (fragile, config-dependent)	solves, 257 steps

Sweep of 12 PlanarMechanics models (pendula, wheels, slip models, sensors, balancing
robot, parallel kinematic robot): 11 pass with identical or better results.
ParallelKinematicRobot previously failed with a stable_eachindex MethodError and now
simulates successfully (the merged-block selection avoids the offending downstream code
path). TwoTrackModelTest fails identically with and without this patch (pre-existing,
unrelated: Infs/NaNs in runtime LU during solve).

New regression tests distill the pattern into 6-variable structural fixtures: the merge
test, end-to-end selections for both priority assignments (the swapped-priority case
passes on the unpatched code precisely because the old selection is priority-blind), and
a Jacobian-path case where the unpermuted col_order demotes the priority-100 variable.

Behavioral note

In models with non-uniform priorities (MTK always passes priorities), state selection may
change where SCCs now merge — by design, in favor of higher-priority states. For uniform
priorities the candidate ordering within a merged block is the concatenation of the
former per-SCC orders, so different-but-equally-valid dummy choices are possible in
models that previously had coupled blocks.

🤖 Generated with Claude Code

[baggepinnen]

NOTE: see correction to this comment in the comment after

---

Tested this on the MultibodyComponents suspension models (stacked with #99 + #100's family grouping):

• HalfCar (22 states): behavior-neutral. Same state selection, same runtime blocks, init converges, and the (separate, pre-existing) integration instability is unchanged. No regressions.

• FullCar (36 states, ~1400 prepared blocks): compile is OOM-killed (kernel kill at 13.8 GB RSS, 31 GB machine; --heap-size-hint=10G doesn't help). Bisect on otherwise-identical stacks:

  • with 90fc5f8/f19b11a: OOM during mtkcompile, reproducibly;
  • without (same branch minus these two commits): compiles in ~20 min, 36 states / 8125 observed, and the downstream pipeline runs.

  FullCar is presumably the worst case for merge_dummy_derivative_blocks: a free-floating chassis whose orientation couples to four corner assemblies through twice-differentiated connection aliases — exactly the pattern that triggers merging — so the merged Mattsson–Söderlind blocks (and their integer Jacobians / candidate lists) may grow very large. I haven't profiled where the memory goes; happy to re-run with instrumentation if useful. A size guard with fallback to the unmerged per-SCC behavior would make this safe to land even before the scalability question is settled.

For calibration: the no-#103 stack takes FullCar further than it has ever been — compiles, initialization residuals all finite, and every runtime linear block consistent (220×220 with rank 216 = the four corners' gauge dimensions, relres 1e-13; 568×568 full-rank, relres 1e-15). The remaining wall is the same integrator-level failure as HalfCar (FBDF: SingularException(36) on W from a near-consistent state), which is independent of this PR.

[baggepinnen]

Correction to my previous comment: the FullCar OOM was not this PR.

After adding an equation budget to #100's family grouping (22b1ab8 there — an unbounded greedy could chain the four corner families through the free chassis into one enormous joint block) and re-running the bisect on an otherwise idle machine, all three configurations compile FullCar cleanly (~20 min, 36 states / 8125 observed, identical family formation):

stack	FullCar compile
main + #99 + #103 only	✓
main + #99 + #100 (no #103)	✓
main + #99 + #100 + #103	✓

The two failures I attributed to #103 were (1) the unbudgeted #100 grouping (13.8 GB kernel kill — genuinely reproducible, but #100's bug, now fixed there) and (2) a re-run on a machine that at the time had only ~13 GB available. Classic bisect-under-confounders mistake — apologies for the noise.

So: #103 is compile-neutral on both car models, behavior-neutral on HalfCar, and the merged dummy-derivative selection on FullCar produces the same state count (36). No scalability concern from these models.

[AayushSabharwal]

While this seems reasonable, I'm concerned it's masking a problem in how Pantelides ignore state priorities. I also think this might cause us to blow up the size of some SCCs and make tearing's job that much harder? Is this really the best solution, or can we:

• do something better in alias elimination so that the alias equations don't exist
• make Pantelides better?

Alias elimination would remove this equation, but both sides have a positive state priority so it doesn't. In such cases, is it valid to unconditionally replace (substitute) the lower priority variable for the higher priority one? I seem to recall doing so ran into issues, but can't recall what those issues were.

[baggepinnen]

but both sides have a positive state priority so it doesn't

I think this is problematic, I often use varying levels of positive priority to indicate preference and always want the highest one to win and any lower one to be optimally eliminated, even though it has positive priority

[AayushSabharwal]

Okay, I can make that change. What if we get two with equally high state priority? Arbitrary choice/tiebreaker with the canonical_rank PR?

[baggepinnen]

Yeah if equal then canonical would be good. if we have the concept of StateSelect.Always, then it should perhaps error if there are two of those to choose from?