Fix GaussAdjoint on immutable (SVector) ODE state

[ChrisRackauckas-Claude]

Important

This PR should be ignored until reviewed by @ChrisRackauckas.

Fixes #1461.

Problem

GaussAdjoint (the default sensealg choice for many out-of-place problems) errored with AssertionError: sensealg isa QuadratureAdjoint in ReverseLossCallback when the ODE state is immutable (e.g. a StaticArrays.SVector). Fixing the assertion exposed three further bugs in the immutable-state path, all fixed here:

1. The asserts themselves (src/adjoint_common.jl): the immutable-state branches of ReverseLossCallback are valid for any adjoint whose integrator state is exactly λ with no parameter-gradient augmentation — exactly what the existing isq flag encodes. isq is widened to AbstractGAdjoint (it previously missed GaussKronrodAdjoint) and the asserts now check isq.
2. Integrand callback arity (src/gauss_adjoint.jl): DiffEqCallbacks' IntegratingSumCallback/IntegratingGKSumCallback call the integrand as integrand_func(u, t, integrator) (allocating) when the adjoint problem is out-of-place, but the 4-arg mutating form was always passed, giving a MethodError. The allocating form is now passed for immutable state.
3. Checkpointing in the out-of-place split_states (src/gauss_adjoint.jl): the OOP split_states ignored checkpointing — it interpolated the (possibly non-dense) forward solution and never refreshed checkpoint_sol/GaussInt.sol, so the quadrature integrand extrapolated a stale checkpoint solution and silently produced garbage dG/dp whenever checkpointing was active (e.g. a saveat forward solution; ischeckpointing(::GaussAdjoint, sol) = checkpointing || !sol.dense). The checkpoint-refresh logic is factored into Gaussupdate_checkpoint_sol! and the OOP split_states is now checkpoint-aware. (Caught by the new tests: du0 looked fine while dp was wrong by orders of magnitude.)
4. Immutable accumulation buffers (src/gauss_adjoint.jl, src/parameters_handling.jl): with static parameters (e.g. SVector p), the quadrature accumulation buffers derived via allocate_zeros/allocate_vjp were themselves static, and DiffEqCallbacks' recursive_zero! then hit setindex!(::SVector, ...). A new mutable_buffer helper (identity for mutable arrays) makes the accumulation buffers mutable.

The mutable-state paths are unchanged (all new branches reduce to the previous behavior for mutable state/parameters).

Verification (all run locally)

• Issue MWE on master reproduces the AssertionError at src/adjoint_common.jl:753; after this PR the assertion no longer fires.
• SVector Lotka–Volterra, Zygote.gradient through solve vs ForwardDiff.gradient:
  • GaussAdjoint(EnzymeVJP()) relerr ≈ 5.5e-12
  • GaussAdjoint(ZygoteVJP()) relerr ≈ 5.4e-12
  • GaussAdjoint() (automatic vjp) relerr ≈ 5.4e-12
  • GaussKronrodAdjoint(EnzymeVJP()) relerr ≈ 5.6e-12
• The issue's SimpleChains + SVector neural-ODE MWE now works with GaussAdjoint(ZygoteVJP()) and with the default automatic sensealg choice, under both outer Zygote.gradient and outer Enzyme.gradient (outer Enzyme ≡ outer Zygote exactly; both match ForwardDiff to Float32 tolerance, maxreldiff ≈ 1e-6).
• test/adjoint_oop.jl passes in full locally (all pre-existing tests plus the new regression tests, which compare discrete/continuous adjoints — including the checkpointing path — against ForwardDiff references).
• test/gauss_zygote_inplace.jl passes (mutable-path smoke test).
• Runic formatting applied to the touched files.

Out of scope (pre-existing, shared with QuadratureAdjoint)

• The MWE's literal GaussAdjoint(autojacvec = EnzymeVJP()) still fails on SimpleChains models, but for an unrelated upstream reason: Enzyme cannot differentiate SimpleChains' TurboDense SIMD kernels (MethodError: from_tape_type(::Type{VecElement{Float32}})). Verified that QuadratureAdjoint(EnzymeVJP()) fails identically on the same MWE, so this is an Enzyme.jl/SimpleChains limitation, not an adjoint-path issue. EnzymeVJP on immutable state works fine for plain Julia rhs functions (tested).
• ReverseDiffVJP has no out-of-place _vecjacobian method, so GaussAdjoint(ReverseDiffVJP()) and QuadratureAdjoint(ReverseDiffVJP()) both fail on immutable state with the same MethodError (verified on master). This PR brings GaussAdjoint to parity with QuadratureAdjoint; the ReverseDiff gap is a separate issue.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Re: mutable_buffer being just similar — agreed; replaced it in e04bb13 with a mutable_zeros helper built on similar (so static tunables get a zeroed MVector instead of a heap Vector from collect). test/adjoint_oop.jl re-verified locally with the change.

On whether the allocating out-of-place integrand path costs anything — timed it (10-state SVector system, 100 parameters):

Integrand level (GaussIntegrand, EnzymeVJP):

3-arg allocating  integrand(t, λ):       76.9 μs (1263 allocations: 37.03 KiB)
4-arg in-place    integrand(out, t, λ):  77.2 μs (1261 allocations: 36.12 KiB)

End-to-end adjoint_sensitivities: 36.0 ms / 549,423 allocs vs 36.4 ms / 549,159 allocs — i.e. for AD vjps the per-node output allocation is noise (<1%); the Enzyme.gradient/Zygote.pullback internals dominate each quadrature node.

Where it does matter is cheap integrands. With a trivial 100-element integrand on an out-of-place SVector ODE, the forced allocation path costs 1.6× in time and 3.7× in allocations for the whole solve:

allocating 3-arg integrand:   196.8 μs (865 allocations: 380.70 KiB)
in-place 4-arg integrand:     123.3 μs (253 allocations: 103.38 KiB)

So I opened SciML/DiffEqCallbacks.jl#313 adding an opt-in integrand_inplace = true to IntegratingSumCallback/IntegratingGKSumCallback that forces the in-place 4-arg integrand even for out-of-place problems (the dG/dp buffer is parameter-shaped and mutable regardless of state immutability). Once that's merged/released, a small follow-up here can drop the 3-arg closure and pass integrand_inplace = true unconditionally — I benchmarked this PR's code against the patched DiffEqCallbacks locally and results are identical.

(Two pre-existing gaps noticed while benchmarking, not touched here: vec_pjac! only calls the in-place vjp_p(out, λ, y, p, t) form, but SciMLBase's conformance check requires out-of-place vjp_p for out-of-place f, so analytical vjp_p is currently unusable with OOP problems; and ReverseDiffVJP lacks an out-of-place _vecjacobian, as already noted.)

[ChrisRackauckas-Claude]

Updated in dbfa6ee now that DiffEqCallbacks 4.18 is out: the Gauss integrand always uses the in-place 4-arg form via integrand_inplace = true (the dG/dp buffer is parameter-shaped and mutable regardless of state immutability); compat bumped to DiffEqCallbacks = "4.18".

Also hunted down the remaining removable allocations in OOP mode (same commit):

1. Adjoint rhs (_vecjacobian, every solver stage): the out-of-place EnzymeVJP/ZygoteVJP methods differentiated w.r.t. both u and p, but the Gauss/Quadrature rhs calls them with dgrad === nothing — the parameter gradient was computed and discarded on every stage. Now p is Enzyme.Const / closed over for Zygote when dgrad === nothing, and the extra primal eval f(y, p, t) is skipped when dy === nothing. This was the biggest win.
2. Quadrature integrand (vec_pjac!, both Gauss and Quadrature): the OOP Enzyme branch used Enzyme.gradient, allocating a gradient vector + result tuple per node; it now writes directly into out via Enzyme.autodiff with Duplicated(tunables, out) when the buffer types match (the common Vector-parameters case; static tunables keep the gradient fallback).

Seeded A/B (10-state SVector ODE, 100 params, discrete-cost adjoint_sensitivities, EnzymeVJP; identical dp checksums before/after):

	before	after
GaussAdjoint	35.86 ms, 549,423 allocs, 15.87 MiB	31.37 ms, 462,095 allocs, 13.18 MiB (−12.5% / −16%)
QuadratureAdjoint	42.90 ms, 651,695 allocs, 19.20 MiB	37.68 ms, 564,475 allocs, 16.56 MiB (−12% / −13%)

Verified locally on the final tree: test/adjoint_oop.jl passes in full, and the SimpleChains SVector MWE gives outer-Enzyme ≡ outer-Zygote gradients matching ForwardDiff (maxreldiff ≈ 3e-6, Float32).

What's left allocates inside the Enzyme/Zygote runtimes per vjp call, not in our plumbing. Two further candidates noted but out of scope here: QuadratureAdjoint integrates with allocating quadgk (the S(t) integrand allocates similar(tunables) per evaluation — could move to quadgk!), and analytical vjp_p is currently unusable with OOP f (SciMLBase requires conforming OOP vjp_p but vec_pjac! only calls the 5-arg in-place form).