Skip to content

Fix "cannot evaluate norm recursively" in LowStorageRK / SSPRK under RAT v4#3513

Closed
ChrisRackauckas-Claude wants to merge 1 commit intoSciML:masterfrom
ChrisRackauckas-Claude:fix-norm-recursive-lowstorage-ssprk
Closed

Fix "cannot evaluate norm recursively" in LowStorageRK / SSPRK under RAT v4#3513
ChrisRackauckas-Claude wants to merge 1 commit intoSciML:masterfrom
ChrisRackauckas-Claude:fix-norm-recursive-lowstorage-ssprk

Conversation

@ChrisRackauckas-Claude
Copy link
Copy Markdown
Contributor

Summary

Root cause

The failing assertion in both suites is @test sol_SA ≈ sol_SV, where sol_SA and sol_SV are ODESolutions whose internal u storage is a Vector{VectorOfArray{Float64, 2, StructVector{SVector{1,Float64},…}}} vs Vector{VectorOfArray{Float64, 2, Vector{SVector{1,Float64}}}}. ODESolution <: AbstractVectorOfArray <: AbstractArray, so the generic LinearAlgebra.isapprox(::AbstractArray, ::AbstractArray) fallback calls norm(sol). On Julia 1.12 + RecursiveArrayTools v4, norm now routes through a stricter norm_recursive_check, which throws for this container shape because AbstractArray iteration yields the solution itself (iterate(sol)[1] === sol) when the leaf storage is a StructArray{<:SVector}-backed VectorOfArray. No OrdinaryDiffEq solver code appears on the failing stack — solve completes fine for both storages; the failure is purely in at the test site.

Fix

  • Add Base.isapprox(x::AbstractTimeseriesSolution, y::AbstractTimeseriesSolution) in lib/OrdinaryDiffEqCore/src/misc_utils.jl. It:
    • checks trajectory lengths and compares x.t ≈ y.t,
    • iterates x.u / y.u snapshot-by-snapshot and delegates to isapprox on leaves,
    • materialises VectorOfArray / StructArray snapshots (via parent(...) + collect) into plain Arrays before the inner isapprox, so LinearAlgebra.norm never sees a container whose element type equals itself.
  • Import LinearAlgebra.norm (and the LinearAlgebra module) in OrdinaryDiffEqCore.
  • Bump OrdinaryDiffEqCore to 4.0.1.

The scalar and Vector{Float64} comparison paths still dispatch to isapprox on the raw snapshot, so the existing @test sol_old.u[end] ≈ sol_new.u[end] assertions elsewhere in the LowStorageRK / SSPRK tests are untouched.

Tests were not removed, disabled, marked @test_broken, or commented out.

Test plan

Locally on Julia 1.12.6 (matching CI), monorepo Pkg.developed against this branch:

  • OrdinaryDiffEqLowStorageRK VectorOfArray/StructArray compatibility testset: 16/16 pass (was 15 pass / 1 error on master).
  • OrdinaryDiffEqSSPRK VectorOfArray/StructArray compatibility testset: 2/2 pass (was 1 pass / 1 error on master).
  • Pre-fix reproducer (same test code, on master tip 72052c030) confirms the ArgumentError: cannot evaluate norm recursively … with stack trace identical to the CI failure.
  • Full CI on this PR (LowStorageRK + SSPRK + rest).

🤖 Generated with Claude Code

…RAT v4

The VectorOfArray/StructArray compatibility testsets in
`lib/OrdinaryDiffEqLowStorageRK/test/ode_low_storage_rk_tests.jl` and
`lib/OrdinaryDiffEqSSPRK/test/ode_ssprk_tests.jl` assert
`sol_SA ≈ sol_SV` where both sides are `ODESolution`s whose `u` is a
`Vector{VectorOfArray{Float64, 2, StructVector{SVector{1,Float64},…}}}`
vs `Vector{VectorOfArray{Float64, 2, Vector{SVector{1,Float64}}}}`.
`ODESolution <: AbstractVectorOfArray <: AbstractArray`, so with no
`isapprox` override the generic `LinearAlgebra.isapprox(::AbstractArray,
::AbstractArray)` fallback calls `norm(sol)`, which — under
RecursiveArrayTools v4 / Julia 1.12 — hits a stricter
`norm_recursive_check` and throws

  ArgumentError: cannot evaluate norm recursively if the type of the
  initial element is identical to that of the container

because `AbstractArray` iteration of the solution yields the solution
itself (`iterate(sol)[1] === sol`) when the leaf storage is a
`StructArray{<:SVector}`-backed `VectorOfArray`. No OrdinaryDiffEq
solver code is on the failing stack; `solve(...)` completes normally for
both storages. The failure is purely in the `≈` path at the test site.

Fix it by adding an `isapprox(::AbstractTimeseriesSolution,
::AbstractTimeseriesSolution)` method in `OrdinaryDiffEqCore` that
compares `t` and iterates `u` snapshot-by-snapshot, materialising
`VectorOfArray` / `StructArray` snapshots to plain `Array`s before the
inner `isapprox`. This avoids the pathological recursive `norm` path,
sidesteps the cross-container-type issue when the two solutions use
different backing arrays, and preserves the semantics callers expect
(element-wise approximate equality of the trajectory). The scalar /
`Vector{Float64}` paths still dispatch to the existing `isapprox` on the
leaf snapshot, so existing `@test sol_old.u[end] ≈ sol_new.u[end]`
assertions are untouched.

- Adds `Base.isapprox(x, y)` for `SciMLBase.AbstractTimeseriesSolution`
  in `lib/OrdinaryDiffEqCore/src/misc_utils.jl`.
- Imports `LinearAlgebra.norm` (and the `LinearAlgebra` module) in
  `lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl`.
- Bumps `OrdinaryDiffEqCore` to 4.0.1.

Locally verified both offending testsets now pass on Julia 1.12.6
(same tool version as CI) with the monorepo `Pkg.develop`ed against
this branch:

  LowStorageRK VectorOfArray/StructArray compatibility: 16/16 passes
  SSPRK         VectorOfArray/StructArray compatibility:  2/2  passes

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
@ChrisRackauckas-Claude
Copy link
Copy Markdown
Contributor Author

Diagnosis looks right, but the fix lands Base.isapprox(::SciMLBase.AbstractTimeseriesSolution, ::SciMLBase.AbstractTimeseriesSolution) in OrdinaryDiffEqCore — that's type piracy (neither function nor type is owned here), same class of issue as the get_fsalfirstlast piracy fixed in #3505. The natural home for this dispatch is SciMLBase, right next to AbstractTimeseriesSolution where the existing array-interface methods live.

Suggested path:

  1. Move the Base.isapprox(::AbstractTimeseriesSolution, ...) method (plus _isapprox_snapshot, _snapshot_array, _isapprox_materialised, _cmp_eltype) into SciMLBase — probably under src/solutions/ode_solutions.jl or a new src/solutions/comparisons.jl. Use RAT's exported AbstractVectorOfArray / recursive_bottom_eltype (RAT is already a direct SciMLBase dep). Bump SciMLBase patch.
  2. Once that's registered, this PR can shrink to just the OrdinaryDiffEqCore version bump (or be closed entirely — no source change needed here).

Alternative: if the real bug is in RAT's norm_recursive_check being over-eager for StructArray{<:SVector}-backed VectorOfArrays, the cleanest fix is upstream in RecursiveArrayTools (teach norm_recursive_check to see through a StructArray parent) — that would unblock every downstream that trips over the same recursion check, not just OrdinaryDiffEq's tests. Worth confirming which layer owns the bug before picking the landing point.

@ChrisRackauckas-Claude
Copy link
Copy Markdown
Contributor Author

Reproduced and tracked down: the root cause is in SciMLBase, not RecursiveArrayTools. VectorOfArray's own norm works fine.

The broken iterator is SciMLBase/src/solutions/solution_interface.jl:188-192:

function Base.iterate(sol::AbstractTimeseriesSolution, state = 0)
    state >= length(sol) && return nothing
    state += 1
    return (solution_new_tslocation(sol, state), state)
end

For a solution typed <: AbstractArray{Float64, 3} with eltype = Float64, this yields another ODESolution of the same concrete type (same fields, different tslocation). Under Julia 1.12, LinearAlgebra.norm_recursive_check iterates the container, checks v == itr, and throws because the first element has the container's own type — and AbstractVectorOfArray's == compares entries, so a "copy with tslocation rewound to 1" compares equal to the full container. That's the exact preconditions the check is guarding against.

This iterate has been like this for years (for plotting/tslocation-stepping semantics), so I don't want to change its behavior. The upstream fix is to add LinearAlgebra.norm(::AbstractTimeseriesSolution) in SciMLBase that operates on sol.u directly — bypassing both the Base iterate-based fallback and the recursive-check trip. One method, one file, no test piracy, no downstream churn.

I'll close this PR once the SciMLBase patch is in — or I can repurpose it to just be the version bump of OrdinaryDiffEqCore once SciMLBase ships the fix. WDYT?

@ChrisRackauckas-Claude
Copy link
Copy Markdown
Contributor Author

SciMLBase upstream fix: SciML/SciMLBase.jl#1322 — replaces the custom iterate(::AbstractTimeseriesSolution) (which returned solution_new_tslocation(sol, state)) with a linear-indexing form that yields scalars per the AbstractArray contract. That lets LinearAlgebra.norm(sol) + isapprox(sol_SA, sol_SV) run on Julia 1.12 without hitting norm_recursive_check.

Once #1322 merges and registers, this PR should close: the Base.isapprox piracy in OrdinaryDiffEqCore becomes unnecessary.

@ChrisRackauckas-Claude
Copy link
Copy Markdown
Contributor Author

Superseded by SciML/SciMLBase.jl#1322 (registering as SciMLBase v3.4.0). The Base.isapprox piracy in OrdinaryDiffEqCore is no longer needed — with 3.4.0, iterate(::AbstractTimeseriesSolution) defers to the AbstractArray fallback, so LinearAlgebra.norm(sol) / sol ≈ sol stop tripping norm_recursive_check on Julia 1.12. Closing.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants