[ChrisRackauckas-Claude]

SCC-lowered SteadyStateProblem solves rebuild sol.u in the original unknown order via variable_symbols, recompute sol.resid from the original residual at that u (the lowering's residual can be empty or shorter when states are observed/eliminated — remapping alone cannot recover eliminated residuals), and apply save_idxs after that restore. ModelingToolkit test compat is raised to 11.45 (true minimum 11.45.0: MTK #5163's empty-SCC schedule fix; Downgrade CI failed on 11.44.0 with SCCDecomposition(...; combine_sccs) when eq_sccs was Vector{Any}).

Failing before (src unfixed, resid/save_idxs assertions):

SteadyStateProblem SCC solutions use the original state order | 10 passed, 5 failed, 5 errored, 20 total
Some tests did not pass: 10 passed, 5 failed, 5 errored, 0 broken.

Passing after:

SteadyStateProblem SCC solutions use the original state order | 20 passed, 20 total
DynamicSS on a SteadyStateProblem with an SCC lowering | 4 passed
SICNM on a SteadyStateProblem with an SCC lowering | 4 passed

Validation: focused test/scc.jl run (before/after with src stashed); Runic on changed files; typos; git diff --check.

Not verified: full GROUP=Core Pkg.test after this follow-up; docs; downstream; Downgrade CI re-run on the new floor.

Reviewer pushback: confirm recomputing resid from prob.f (vs attempting to permute the lowering resid) is the right contract; non-lowered solve(...; save_idxs=...) still ignores save_idxs (only the SCC restore path applies it).

Fixes #171

Please ignore this draft until reviewed by @ChrisRackauckas.

Risk assessment

• Risk: low
• Blast radius: solutions of SCC-lowered SteadyStateProblems. sol.u is reconstructed in the original state order (via SymbolicIndexingInterface variable_symbols), sol.resid is recomputed from the original prob.f at the restored u, and save_idxs is applied after the restore. Plain (non-MTK) SteadyStateProblems give identical results to base. The test-only MTK compat floor is raised to 11.45; 11.45.0 is the first release with Fix solve(::SteadyStateProblem) regressions from the deferred SCC lowering ModelingToolkit.jl#5163 combine_sccs support. No public API change.
• Evidence: on base, the LV sol.u comes back reversed and the chain and mix cases throw DimensionMismatch. On the PR, sol.u, sol[x] and resid match the non-lowered solve for all three, and save_idxs=[2] returns only element 2. All 13 CI checks are green, including Downgrade, which resolves MTK 11.45.0.
• Independent review: Claude Code (Opus 5.5), two rounds. Round 1: medium, CHANGES (resid layout, save_idxs, compat floor). Round 2 on 680dcc2: low, medium-high confidence, MERGE. SSRootfind(NewtonRaphson()) on the chain system throws a ForwardDiff chunk-size error on base too; that's not caused by this PR.
• Merge: auto-merge candidate. The authors were Codex gpt-6-luna, then Cursor Auto.

Risk assessment and review coordination by 🤖 Claude Code (model: claude-opus-5-5[1m]), https://claude.ai/code/session_01NdFVhANTF4nepTbqYxjfRP

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Risk assessment (independent review, Codex gpt-6.1-sol, head 680dcc2)

• Verdict: CHANGES
• Rating: Codex gpt-6.1-sol rated it high (H1) at head 680dcc2
• Risk: high (H1: changes the explicitly documented returned problem/state ordering for SCC-lowered SteadyStateProblems; downstream positional access and sol.prob consumers observe a different contract). Also matches L2 for added tests and ML4 for the raised compat floor; H1 is the highest matched criterion.
• Blast radius: MTK/symbolic SteadyStateProblems with a stored lowering, through bare nonlinear algorithms, SSRootfind, DynamicSS and SICNM. The new save_idxs handling can silently assign one variable's value to another variable. No new public names or dependencies are introduced; variable_symbols is exported by its owning dependency, and no new copyleft dependency or skipped/loosened assertion appears in the diff.
• Evidence: Read issue 171, the full git diff origin/HEAD...HEAD, algorithm docstrings, relevant solution-indexing implementation, and current CI metadata. All 13 reported checks are successful, so there is no failing CI check requiring a default-branch comparison. Ran the PR's new 20-assertion testset unchanged in scratch copies: head 20 passed; default branch 471a4a039f637459ff8625cc832f20a7b85f4f49 7 passed, 7 failed, 6 errored. Ran an additional symbolic-indexing probe: default branch 4 passed; head 2 passed, 2 failed, with both algorithm variants swapping symbolic values. The two passing head assertions independently check the analytic LV equilibrium [u,v] = [3,1.5]. Julia 1.12.4; ModelingToolkit 11.45.4; ModelingToolkitBase 1.77.4; SciMLBase 3.58.0; NonlinearSolve 4.32.0. git diff --check origin/HEAD...HEAD passed, and final git status --porcelain was empty. G1 pass: analytic equilibrium checks and comparison with the monolithic, non-lowered solve provide independent numerical evidence. G2 not applicable: no GPU kernel, AD rule, Reactant tracing or MPI implementation changes. G3 not applicable: no shared-state/cache/threading/default change is introduced by this diff.
• Merge: needs changes (repair symbolic indexing after save_idxs; update the affected documentation and handle the documented behavior change under the H1 major-release policy).

Targeted commands run (S denotes <local path>):

timeout 3600 julia --project="$S/env" "$S/regression.jl"
timeout 3600 julia --compiled-modules=existing --project="$S/env" "$S/adversarial.jl"
timeout 3600 julia --compiled-modules=existing --project="$S/base-env" "$S/base-probes.jl"

All used JULIA_NUM_THREADS=4, JULIA_PKG_PRECOMPILE_AUTO=0, TMPDIR="$S/tmp", and JULIA_DEPOT_PATH=<local path>. The third script runs the same adversarial probe followed by the same extracted PR testset; its nonzero exit is the expected unfixed regression failure. Scratch copies came from git archive HEAD and git archive origin/HEAD, with identical dependency versions and only the developed SteadyStateDiffEq source path differing.

---

🤖 Posted by an AI agent — harness: Codex CLI 0.162.1 · model: gpt-6.1-sol
Conversation: local headless review run (no public session URL); transcript kept by the fleet master.

[ChrisRackauckas]

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Independent review (Devin CLI 3000.11.3, model fusion-claude-opus-5-5-high-sidekick-swe-2-medium): CHANGES, risk medium.

Full review

VERDICT: CHANGES
RISK: medium

The core fix works. I ran it, and on MTK-built problems sol.u now comes back in prob.u0 order and sol.resid is the original residual at that u. This holds for NewtonRaphson, SSRootfind, DynamicSS, SICNM and a bare solve. It also covers a completed system whose lowering has fewer unknowns and a 4-block SCC system. On base, all of these return the lowering's order. A ForwardDiff derivative through solve(...).u[i] is wrong on base (0.0588 vs the analytic 0.2100) and right on the PR (0.20999). The restore path also has two regressions against base, both confirmed by running them. Both need small changes.

Blocking findings

1. save_idxs on the NewtonRaphson/SSRootfind SCC path silently returns the wrong value for symbolic indexing. Base gave the right value.
   src/solve.jl:32-38, reached from src/solve.jl:55-58 (SSRootfind __solve) and src/solve.jl:195-197 (top-level solve).

   • Confirmed by running (scratch/probe3.jl): Lotka–Volterra with unknowns(lv) == [v, u], true v = 1.5, u = 3.0. I ran s = solve(prob, NewtonRaphson(); save_idxs = [2]).

     • PR: s.u == [3.0], s[v] == 3.0 (wrong, and no error), s[u] throws BoundsError. Same for SSRootfind(NewtonRaphson()).
     • Base: save_idxs was ignored on this path, so s[v] == 1.5 and s[u] == 3.0 (right).

   • Cause: the restored solution is built on the full SteadyStateProblem, so its symbolic container still maps v → 1 and u → 2, but the u vector holds only the saved entry.

   • It also makes the same call give different shapes depending on whether a lowering exists. remake(prob; lowered_problem = nothing) with NewtonRaphson and save_idxs = [2] returns a 2-vector, and the lowered problem returns a 1-vector.

   • The new test (test/scc.jl:441-446) checks only saved.u/saved.resid, so it cannot catch this.

   • DynamicSS has the same symbolic-indexing defect with or without a lowering, on base too. That part is pre-existing. The PR's new part is extending it to the nonlinear-solver paths, which used to be correct.

   • Required change, pick one:

     • (a) Don't apply save_idxs in the restore when the call came through the NonlinearSolve/SSRootfind paths. That matches the non-lowered behaviour and base.
     • (b) If save_idxs must be honoured, build the subset solution so symbolic indexing can't misread it, e.g. without the full-system symbolic container, or with a saved-subsystem mapping like ODE solutions use.

     In both cases, add a test asserting that saved[states[k]] matches reference[states[k]], or that it throws. It must not silently return another state's value.

2. A non-MTK symbolic sys plus a hand-built SCC lowering now throws on every algorithm. Base solved it.
   src/solve.jl:26-30.

   • The guard hasproperty(prob.f, :sys) is always true for SciMLFunctions, because the field exists even when it is nothing. When sys is nothing, the length(states) == length(prob.u0) check is what saves the hand-built-lowering tests.
   • When prob.f.sys is a SymbolCache, variable_symbols(prob) has the right length. The restore then calls sol[states] on an SCC solution that has no symbolic information.
   • Confirmed by running (scratch/probe4.jl): SteadyStateProblem(ODEFunction(f; sys = SymbolCache([:a,:b,:x,:y])), zeros(4); lowered_problem = sccprob).
     • PR: ERROR SciMLFunction does not have observed for NewtonRaphson, SSRootfind, DynamicSS and SICNM.
     • Base: all four return Success, u ≈ [1,2,1,2].
   • Required change: before indexing, check that the lowered solution resolves every state, e.g. all(s -> is_variable(sol, s) || is_observed(sol, s), states) (SymbolicIndexingInterface public API). Return nothing (the fallback) otherwise. Drop the vacuous hasproperty guard. Add this case to the "SteadyStateProblem with an SCC lowering" testset.

Non-blocking findings

• Overlaps with sibling PR Map SCC-lowered steady-state solutions back onto the problem's state (stacked on #173) #177.
  • Map SCC-lowered steady-state solutions back onto the problem's state (stacked on #173) #177 is by the same agent account, still a draft, and stacked on Add __init for DynamicSS/SICNM so init() doesn't fall through to NonlinearSolve's default #173. Its __scc_solution_on_problem maps SCC solutions back onto prob in the same functions (__solve_scc_lowering and the SSRootfind __solve). It also recomputes resid from prob.f(u, p, Inf) and raises the same MTK floor to 11.45.
  • The two will conflict textually. Map SCC-lowered steady-state solutions back onto the problem's state (stacked on #173) #177 explicitly leaves the top-level solve(prob, NewtonRaphson()) path, i.e. solve(::SteadyStateProblem) returns the lowered problem's solution: sol.u reordered/empty, sol.prob.f.sys time-independent #171, unfixed. This PR covers that path.
  • Map SCC-lowered steady-state solutions back onto the problem's state (stacked on #173) #177 already does the resolvability check that finding 2 asks for.
  • One of the two should absorb the other before either merges. Merging both means rewriting the same code twice.
• Suspected from reading: the silent fallback at src/solve.jl:29 (length(states) == length(prob.u0) || return nothing) brings back the reordered lowering solution with no warning. I could not construct an MTK problem that reaches it.
• init(prob, NewtonRaphson()); solve! on the LV SCC problem errors with type ODEFunction has no field resid_prototype, on base and on the PR alike. This is not caused by this PR. Add __init for DynamicSS/SICNM so init() doesn't fall through to NonlinearSolve's default #173/Map SCC-lowered steady-state solutions back onto the problem's state (stacked on #173) #177 handle init for DynamicSS/SICNM only.
• solve(prob, NewtonRaphson(); p = remake(prob; p = [k => 8.0]).p) errors with MethodError: current_time(::GeneralizedFirstOrderAlgorithmCache…) on base and on the PR alike. This is pre-existing. Passing the new value through remake instead gives the right answer on the PR (a = 2.0 = cbrt(8)).
• The PR body's "failing before" counts (10 pass / 5 fail / 5 error on the new testset) don't match exactly what I see: 8 pass / 6 fail / 6 error, with MTK 11.45.3 and SciMLBase 3.57.0. The likely cause is the dependency versions. The claim that the testset fails on base holds.
• The MTK compat floor 11.44 → 11.45 is test-only ([extras]) and justified. The Downgrade log for the first commit (https://github.com/SciML/SteadyStateDiffEq.jl/actions/runs/36347793019) shows MTK 11.44.0 failing at test/scc.jl:404 with MethodError: no method matching ModelingToolkit.SCCDecomposition(...; combine_sccs::Bool). The rerun (https://github.com/SciML/SteadyStateDiffEq.jl/actions/runs/36371936324) resolves MTK 11.45.0 and passes.
• The CI job Core (julia pre) reports "pass" in 28 s, so it likely ran no tests. I did not investigate further.
• Added comments are about 7 of ~108 added lines and are present-tense. No new dependency (SymbolicIndexingInterface was already a dep, and variable_symbols is public SII API). No new exports.
• sol.prob changes from SCCNonlinearProblem to SteadyStateProblem for MTK problems, and the PR has no version bump or docs note. It is a bug-fix restoration of pre-2.18 semantics, so a patch release is defensible, but it is user-visible.

What I ran

All runs used Julia 1.11 in scratch/env: PR head 680dcc2 dev'd, MTK 11.45.3, NonlinearSolve 4.32.0, SciMLBase 3.57.0, SCCNonlinearSolve 1.15.4, OrdinaryDiffEq 7.8.1. Base runs used scratch/envbase, the same environment with origin/master 735cb01 dev'd (extracted via git archive).

• test/scc.jl on the PR: 431/431 pass.
• The PR's test/scc.jl against base src: 416 pass, 9 fail, 6 error. The failures are the new testset and the updated sol.prob isa SteadyStateProblem assertions.
• scratch/probe.jl, on PR and base:
  • 4-block SCC system (a^3, b^3, …), with NR, SSRootfind, DynamicSS and SICNM.
  • Full solves and save_idxs ∈ {[2], [4], [1,3], 3}. PR values match the non-lowered reference to ≤2.5e-13 (NR) / ≤1.5e-11 (DynamicSS). Base is reversed, with error 0.58.
  • maxiters = 1: on the PR, resid equals prob.f(u) exactly. On base, resid was in lowering order.
  • Out-of-place problem, u0 kwarg, remake p, the issue solve(::SteadyStateProblem) returns the lowered problem's solution: sol.u reordered/empty, sol.prob.f.sys time-independent #171 reconstruction (SteadyStateProblem(sol.prob.f.sys, op) now works), plain non-MTK save_idxs, and ForwardDiff.
• scratch/probe2.jl: a completed (uncompiled) 2-state system whose lowering has fewer unknowns, from SciMLBase's comprehensive_indexing.jl. Under all 5 algorithm choices, sol.u == [5.0, 1.5], resid == [0, 0], sol[XY] == 6.5.
• scratch/probe3.jl: finding 1, on PR and base.
• scratch/probe4.jl: finding 2, on PR and base.
• julia +1.12 --project=@runic -m Runic --check on the 3 changed Julia files: clean. git diff origin/master...HEAD | typos -: clean.
• gh pr checks 175: all green. I also read the Downgrade logs for both commits.

What I did not verify

• The full GROUP=Core suite and GROUP=QA locally. CI reports both green; I read only the job status and the Downgrade logs.
• The docs build, and downstream packages: SciMLBase GROUP=SymbolicIndexingInterface, and MTK's test/scc_nonlinear_problem.jl.
• Julia LTS (1.10) and the Downgrade minimum environment locally.
• Whether ODEProblem(sys, ...; u0 = sssol.u) round-trips. My probe gave [5.0, 0.5] for a steady state of [5.0, 1.5], most likely because initialization overrides u0. I did not investigate this, and it is probably unrelated to the PR.
• Performance cost of the extra residual evaluation and of sol[states] (observed-function codegen) on large systems.

---

🤖 Posted by an AI agent — harness: Devin CLI 3000.11.3 (review), Claude Code 2.1.285 (posting) · model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Conversation: local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master)

[ChrisRackauckas]

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Hold: Chris's direction (2026-10-14) is to replace the post-solve reversal with an SCC form of SteadyStateProblem whose solution is stored in the compressed SCC layout, as SCCNonlinearProblem solutions are, with symbolic indexing mapping names onto it. A worker is designing that across SciMLBase, ModelingToolkit and SteadyStateDiffEq. This PR and #177 are expected to be superseded; their test cases (LV order, chain, mix, resid, save_idxs) are being carried over.

---

🤖 Posted by an AI agent — harness: Claude Code 2.1.285 · model: claude-opus-5-5
Conversation: local session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master)