BREAKING: Make AbstractVectorOfArray <: AbstractArray

[ChrisRackauckas-Claude]

Summary

This PR completes the long-planned migration (documented in the codebase since 2023) to make AbstractVectorOfArray{T, N, A} <: AbstractArray{T, N}. This is a breaking major update that changes fundamental indexing and iteration behavior.

Breaking Changes

Linear Indexing

• A[i] now returns the ith element in column-major order (standard AbstractArray behavior)
• Previously A[i] returned A.u[i] (the ith inner array)
• Migration: Use A.u[i] or A[:, i] to access inner arrays

Ragged Arrays (Sparse Rectangular Interpretation)

• size(A) reports maximum size in each dimension across all inner arrays
• Out-of-bounds ragged elements return zero(T) (sparse representation)
• Array(A) on ragged data produces a zero-padded dense array
• This means all standard linear algebra just works on ragged VectorOfArrays

Iteration & Length

• iterate(A) goes over scalar elements (column-major), not inner arrays
• length(A) returns prod(size(A)), not length(A.u)
• first(A) / last(A) return scalar elements, not inner arrays
• Migration: Use A.u for column-wise iteration, length(A.u) for count

map

• map(f, A) maps over elements, not inner arrays
• Migration: Use map(f, A.u) to map over inner arrays

Other

• eachindex(A) returns CartesianIndices(size(A))
• convert(AbstractArray, A) is identity (A is already an AbstractArray)
• Removed all deprecated linear indexing warnings
• parameter_values(::AbstractDiffEqArray, i) now dispatches correctly (fixes ambiguity with parameter_values(::AbstractArray, i) from SymbolicIndexingInterface)

Implementation Details

Core AbstractArray Interface

• size(VA) uses maximum(size(u, d) for u in VA.u) for each leading dimension
• getindex(A, I::Vararg{Int, N}) checks ragged bounds and returns zero(T) for out-of-storage indices
• setindex! throws ArgumentError when trying to set non-zero values outside ragged storage
• IndexStyle remains IndexCartesian()
• CartesianIndex getindex/setindex! handle ragged zero-padding

Removed/Changed Methods

• Removed: length, eachindex, iterate, first, last, firstindex, lastindex overrides (inherited from AbstractArray)
• Removed: eltype, ndims, axes overrides (inherited from AbstractArray)
• Removed: any, all overrides (inherited from AbstractArray)
• Removed: sum, prod, mapreduce overrides (inherited from AbstractArray)
• Removed: map override (inherited from AbstractArray)
• Removed: All @deprecate linear indexing methods
• Removed: ==(::AbstractVectorOfArray, ::AbstractArray) (AbstractArray default)
• Removed: convert(::Type{AbstractArray}, ...) (identity)
• Fixed: append! uses new_item.u instead of iterating
• Fixed: recursivefill! uses bs.u instead of iterating
• Fixed: zero uses [zero(u) for u in VA.u]

Zygote Extension

• Removed ambiguous broadcast adjoint overrides (Zygote's built-in AbstractArray rules now apply)
• Some Zygote gradient tests marked @test_broken pending Zygote-specific updates
• ForwardDiff continues to work correctly

Test Plan

• All core tests pass (qa, utils, partitions, basic_indexing, interface_tests, table_traits, static_arrays, linalg, measurements, symbolic_indexing)
• Adjoint tests: ForwardDiff passes, Zygote tests marked broken pending update
• Ragged array zero-padding verified
• Linear indexing matches standard AbstractArray behavior
• sum, prod, mapreduce with dims work correctly
• Broadcasting with plain Arrays works correctly

Note on Version Bump

This PR does not bump the version number so that CI tests can run on the current compat bounds. The version bump to v4.0.0 should be done when merging.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Invalidation Analysis & Cleanup (second commit)

Invalidation Results

8 invalidation trees, all minimal:

• 4 from Colon methods for RaggedEnd type (inherent, max 20 mt_backedges)
• 1 from VA[xs...] syntax constructor
• 1 from setindex!(A, v, i::Int, ci::CartesianIndex) (needed for sum(; dims=d))
• 1 from similar(VA, dims::Integer) (needed for VectorOfArray API)
• 1 from any(f, ::ArrayPartition) (unrelated to this PR)

Removed Redundant Overrides

All of these are now inherited from AbstractArray:

• IndexStyle(::AbstractVectorOfArray) (instance method, type method kept)
• size(VA, i), lastindex(VA, d), checkbounds, isassigned
• isapprox, CartesianIndices, adjoint, reshape, vec
• convert(Array, ...), maybeview
• +, -, *, / operator overrides
• 2-arg show for AbstractVectorOfArray

Zygote Fixes

• All 12 adjoint tests now pass (was 4 pass + 8 broken)
• Root cause: VectorOfArray(u) adjoint was returning a VectorOfArray gradient, but downstream vect_pullback indexed it with ȳ[i] which now returns scalars (linear indexing). Fixed by returning .u (plain Vector) as the gradient.
• Removed all getindex/view adjoint overrides — Zygote's built-in AbstractArray rules handle these correctly.

[ChrisRackauckas-Claude]

Subpackage for VA/AP constructors (cc @rokke-git)

The VA[...] and AP[...] shorthand constructors from #546 have been moved to a separate subpackage lib/RecursiveArrayToolsShorthandConstructors. The reason: Base.getindex(::Type{VA}, xs...) supersedes Base.getindex(::Type{T}, vals...) from Base, causing invalidations for every user even if they never use the syntax.

The VA and AP types are still exported from RecursiveArrayTools — only the getindex overrides that enable VA[...]/AP[...] are deferred. Users who want the shorthand syntax opt in with:

using RecursiveArrayToolsShorthandConstructors

This also moves the any/all optimizations for ArrayPartition (partition-level short-circuiting, 779 invalidation children) into the same subpackage.

Result: using RecursiveArrayTools now causes only 1 invalidation tree (from NamedArrayPartition.setindex!, pre-existing), down from 8 at the start of this PR.

[mateuszbaran]

Why is this getting deleted again? I thought that idea was rejected in #540 .

[ChrisRackauckas]

It's a 1.5x optimization that leads to major compile time issues. I think the good solution here is to make it a sublibrary that people can opt into, but not force it to always exist. If people need this optimization, then they can take the compile time hit, but this allows most libraries to avoid the major invalidations.

[ChrisRackauckas]

Are you okay with that solution? It seems best for the mix of use cases.

[mateuszbaran]

OK, if it's still in a subpackage then it's good enough for me.

[JoshuaLampert]

You probably already planned it like this, but I want to mention it anyways: Currently, there are some unreleased commits on master. It would be good to create a (probably final) release of v3, such that the change to v4 only includes the breaking change(s), but not the other non-breaking changes.

[ChrisRackauckas]

Oh good point. I'll make sure to do that.

[JoshuaLampert]

I am testing this branch right now in JoshuaLampert/SimpleDiscontinuousGalerkin.jl#74 to see what breaks because in that package I used some complicated VectorOfArray constructs including ragged arrays. The first thing I noticed is that the following fails, which worked before:

julia> f = VectorOfArray([[1.0], [2.0, 3.0]])
VectorOfArray{Float64,2}:
2-element Vector{Vector{Float64}}:
 [1.0]
 [2.0, 3.0]

julia> f .= zero(f)
ERROR: DimensionMismatch: destination axes (Base.OneTo(2),) are not compatible with source axes (Base.OneTo(1),)
Stacktrace:
 [1] throwdm(axdest::Tuple{Base.OneTo{Int64}}, axsrc::Tuple{Base.OneTo{Int64}})
   @ Base.Broadcast ./broadcast.jl:1086
 [2] copyto!
   @ ./broadcast.jl:983 [inlined]
 [3] copyto!
   @ ./broadcast.jl:947 [inlined]
 [4] copyto!
   @ ~/.julia/dev/RecursiveArrayTools/src/vector_of_array.jl:1174 [inlined]
 [5] materialize!
   @ ./broadcast.jl:905 [inlined]
 [6] materialize!(dest::VectorOfArray{Float64, 2, Vector{Vector{Float64}}}, bc::Base.Broadcast.Broadcasted{RecursiveArrayTools.VectorOfArrayStyle{2}, Nothing, typeof(identity), Tuple{VectorOfArray{Float64, 2, Vector{Vector{…}}}}})
   @ Base.Broadcast ./broadcast.jl:902
 [7] top-level scope
   @ REPL[24]:1
Some type information was truncated. Use `show(err)` to see complete types.

[JoshuaLampert]

I am also not sure about the zero-padding. If I have a ragged array, I usually would not want to have padded zeros. Note that ragged array are the natural data structures for finite element methods with p-adaptivity, where you have different polynomial degrees (and therefore different number of DOFs) in each element. When I want to index into an element, e.g., via f[:, 1] to get all the DOFs in this element, I do not want to have zeros padded.
I see where this change comes from and I think they are necessary for the subtyping of AbstractArray, but I feel like it breaks applications using ragged arrays, which is not so easy to fix. How should I change my code to only get the non-padded arrays? I could of course turn back to using f.u[1], but I want to profit from the multidimensional indexing of RecursiveArrayTools.jl (in my case also because my solution array is allowed to by of type VectorOfArray or could be a regular array (if the user does not use p-adaptivity and we do not have to deal with ragged arrays) and I want to use the same code for both versions rather).
I put some effort into making ragged arrays usable for me in RecursiveArrayTools.jl and I feel like I am loosing the ability to use this feature elegantly with this new breaking release, which would be sad. I would be happy if we could keep an interface, which supports (non-padded) ragged arrays or at least have a solution, which does not require a lot of rewriting and special casing on my end.

[ChrisRackauckas]

I think we need zero-padding as the default way to do things for the ODESolution because (a) 99.99% or something wild of use cases are non-ragged, and (b) the ragged indexing adds major invalidations which slows down everybody by like 10 seconds while making automatic differentiation have a bunch of small bugs that are difficult to work out (since AD is generally written to assume indexing matches the AbstractArray interface and its unwritten assumptions, one of which is the sizing assumption.

That last part then compounds the issue, because what happens is that size either needs to be defined small (i.e. the minimum size, or the size of the first vector), in which case AD overloads build a smaller matrix which succeeds but can give silent incorrect results, or it's defined using the maximum without zero padding, which then makes AD just error. We can then fix this by overloading every operation, which over time we did get a decent amount of coverage for people to not really realize how crazy this is, but it never got full coverage (for example, Cholesky still does weird things) and the more overloads we add, the more invalidations we get from it. We had to make it all not an AbstractArray to ensure we got more correctness but, ultimately it's not the right path for most people.

That said, there are some nice use cases for it. And we did battle test it enough that it ended up being quite a useful interface. So what I can do is keep it alive as a sublibrary as well, again with the other pieces a sublibrary to prevent the invalidations from being in the main path, but make it easy to convert the type into this form and it would just give it different actions. Indeed adaptive meshing is probably the one major use case of that form which always wants to keep the ragged form.

[JoshuaLampert]

I understand and agree that ragged arrays are not that common and am fine with zero-padding by default to avoid the issues you mention as long as there exists a way to opt-in the current interface. Keeping the possibility of non-padded ragged arrays alive by adding a sublibrary and making it easy to convert the type sounds good to me.

[ChrisRackauckas-Claude]

New commit: Ragged arrays, interpolation interface, and plotting delegation

Ragged Arrays (RecursiveArrayToolsRaggedArrays sublibrary)

Per @JoshuaLampert's feedback — RaggedVectorOfArray and RaggedDiffEqArray types that preserve true ragged structure without zero-padding. Isolated in lib/ sublibrary to avoid invalidations.

using RecursiveArrayToolsRaggedArrays
r = RaggedVectorOfArray([[1.0, 2.0], [3.0, 4.0, 5.0]])
r[:, 1]  # [1.0, 2.0] — no zero-padding
r[:, 2]  # [3.0, 4.0, 5.0]
r .= zero(r)  # works correctly (the exact bug reported)

Interpolation Interface on DiffEqArray

New interp and dense fields + callable syntax:

da = DiffEqArray(u, t; interp = LinearInterpolation(t, u), dense = true)
da(0.5)               # interpolate at t=0.5
da([0.1, 0.5, 0.9])  # multiple times → DiffEqArray
plot(da)              # dense interpolated plot

Plotting Delegation

Full-featured AbstractDiffEqArray recipe with idxs, phase plots, denseplot, plotdensity, tspan. Exports interpret_vars, diffeq_to_arrays, solplot_vecs_and_labels etc. so SciMLBase can delegate its recipe (companion PR: SciML/SciMLBase.jl#1297).

Bug Fix

Fixed ragged VectorOfArray .= zero broadcast (DimensionMismatch) by using dest.u[i] instead of dest[:, i] in copyto!.

Invalidations

2 total (pre-existing trivial: show(IOBuffer, Function) + cache). 0 logedges. No regression.

Version bumped to v4.0.0.

[JoshuaLampert]

Suggested change

SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

Should RecursiveArrayTools.jl really depend on SciMLBase.jl? As far as I can see it is not needed and it leads to circular dependencies since SciMLBase.jl depends on RecursiveArrayTools.jl.

[JoshuaLampert]

I just tried some of the tests specifically for ragged arrays we had before using the new RaggedVectorOfArray and some of them failed. For example:

julia> u = RaggedVectorOfArray([[1.0], [2.0, 3.0]])
RaggedVectorOfArray{Float64,2}:
2-element Vector{Vector{Float64}}:
 [1.0]
 [2.0, 3.0]

julia> @test length(view(u, :, 1)) == 1
Error During Test at REPL[19]:1
  Test threw exception
  Expression: length(view(u, :, 1)) == 1
  MethodError: no method matching view(::RaggedVectorOfArray{Float64, 2, Vector{Vector{Float64}}}, ::Colon, ::Int64)
  The function `view` exists, but no method is defined for this combination of argument types.
  
  Closest candidates are:
    view(::AbstractVectorOfArray{T, N, <:AbstractVector{T}}, ::Any...) where {T, N, M}
     @ RecursiveArrayTools ~/.julia/dev/RecursiveArrayTools/src/vector_of_array.jl:986
    view(::AbstractVectorOfArray, ::Any...) where M
     @ RecursiveArrayTools ~/.julia/dev/RecursiveArrayTools/src/vector_of_array.jl:1001
    view(::AbstractArray, ::Any...) where M
     @ Base subarray.jl:211
    ...
  
  Stacktrace:
   [1] top-level scope
     @ REPL[19]:1
   [2] macro expansion
     @ ~/.julia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/Test/src/Test.jl:677 [inlined]
ERROR: There was an error during testing

julia> @test length(view(u, :, 2)) == 2
Error During Test at REPL[20]:1
  Test threw exception
  Expression: length(view(u, :, 2)) == 2
  MethodError: no method matching view(::RaggedVectorOfArray{Float64, 2, Vector{Vector{Float64}}}, ::Colon, ::Int64)
  The function `view` exists, but no method is defined for this combination of argument types.
  
  Closest candidates are:
    view(::AbstractVectorOfArray{T, N, <:AbstractVector{T}}, ::Any...) where {T, N, M}
     @ RecursiveArrayTools ~/.julia/dev/RecursiveArrayTools/src/vector_of_array.jl:986
    view(::AbstractVectorOfArray, ::Any...) where M
     @ RecursiveArrayTools ~/.julia/dev/RecursiveArrayTools/src/vector_of_array.jl:1001
    view(::AbstractArray, ::Any...) where M
     @ Base subarray.jl:211
    ...
  
  Stacktrace:
   [1] top-level scope
     @ REPL[20]:1
   [2] macro expansion
     @ ~/.julia/juliaup/julia-1.12.5+0.x64.linux.gnu/share/julia/stdlib/v1.12/Test/src/Test.jl:677 [inlined]
ERROR: There was an error during testing

julia> u = RaggedVectorOfArray([zeros(1, n) for n in (2, 3)])
RaggedVectorOfArray{Float64,3}:
2-element Vector{Matrix{Float64}}:
 [0.0 0.0]
 [0.0 0.0 0.0]

julia> u[1, :, 2] .= [1.0, 2.0, 3.0]
ERROR: MethodError: no method matching getindex(::RaggedVectorOfArray{Float64, 3, Vector{Matrix{Float64}}}, ::Int64, ::Colon, ::Int64)
The function `getindex` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  getindex(::RaggedVectorOfArray{T, N}, ::Int64...) where {T, N}
   @ RecursiveArrayToolsRaggedArrays ~/.julia/dev/RecursiveArrayTools/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl:360
  getindex(::RaggedVectorOfArray, ::Int64, ::Colon)
   @ RecursiveArrayToolsRaggedArrays ~/.julia/dev/RecursiveArrayTools/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl:418
  getindex(::RaggedVectorOfArray{T, N}, ::Int64) where {T, N}
   @ RecursiveArrayToolsRaggedArrays ~/.julia/dev/RecursiveArrayTools/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl:321
  ...

Stacktrace:
 [1] maybeview(::RaggedVectorOfArray{Float64, 3, Vector{Matrix{Float64}}}, ::Int64, ::Function, ::Vararg{Any})
   @ Base ./views.jl:148
 [2] dotview(::RaggedVectorOfArray{Float64, 3, Vector{Matrix{Float64}}}, ::Int64, ::Function, ::Vararg{Any})
   @ Base.Broadcast ./broadcast.jl:1250
 [3] top-level scope
   @ REPL[22]:1

I stopped trying out more tests. So there are likely more. I would just suggest to keep the tests, which used ragged arrays before and now test them for RaggedVectorOfArray.

[ChrisRackauckas-Claude]

@JoshuaLampert Thanks for testing — fixed in the latest push:

• view(r, :, i) and view(r, I, i) now work
• axes(r, d) and lastindex(r, d) added so end works in indexing expressions
• Range indexing r[1:3, i], r[:, 2:end] etc. now work

One difference from v3: end in the first dimension now uses the max inner array length (not per-column like v3's RaggedEnd), so r[end, i] will throw BoundsError if column i is shorter than the max. Use r[:, i] for safe column access, or index explicitly with length(r.u[i]).

243 tests pass now including ported v3 ragged behavior tests.

[ChrisRackauckas]

Alright, I think this did it all. @JoshuaLampert and @mateuszbaran, do follow up if the ragged versions aren't matching the previous behavior. But I think we hit every requirement here using the sublibraries.

[ChrisRackauckas-Claude]

@JuliaRegistrator register

Release notes:

RecursiveArrayTools v4.0.0

Breaking Changes

• AbstractVectorOfArray{T, N, A} now subtypes AbstractArray{T, N}. Linear indexing A[i] returns the ith element in column-major order (previously returned the ith inner array). Use A.u[i] or A[:, i] for the old behavior.
• length(A) returns prod(size(A)) (total elements), not length(A.u) (number of inner arrays). Use length(A.u) for the old behavior.
• iterate(A) iterates over scalar elements, not inner arrays. Use A.u for column-wise iteration.
• map(f, A) maps over elements, not inner arrays. Use map(f, A.u) for the old behavior.
• Ragged arrays: size(A) reports maximum size; out-of-bounds elements return zero (sparse interpretation). Non-zero-padded ragged arrays are available via using RecursiveArrayToolsRaggedArrays.
• DiffEqArray has new interp and dense fields (new type parameter I).
• Optimized any/all for ArrayPartition moved to RecursiveArrayToolsArrayPartitionAnyAll subpackage to avoid invalidations. Required for GPU arrays.

New Features

• DiffEqArray callable interface: da(t), da(t; idxs=1), da(t, Val{1}) for interpolation.
• Dense plotting: plot(da) with denseplot=true generates smooth interpolated curves.
• Full-featured plot recipe with idxs, phase plots, tspan, plotdensity, plotat.
• RaggedVectorOfArray and RaggedDiffEqArray in RecursiveArrayToolsRaggedArrays sublibrary for non-zero-padded ragged arrays.
• CUDA extension for CuArray(::AbstractVectorOfArray) disambiguation.
• Plotting helpers exported for SciMLBase delegation: interpret_vars, diffeq_to_arrays, solplot_vecs_and_labels.

[JuliaRegistrator]

Error while trying to register: Register Failed
@ChrisRackauckas-Claude, it looks like you are not a publicly listed member/owner in the parent organization (SciML).
If you are a member/owner, you will need to change your membership to public. See GitHub Help

[ChrisRackauckas-Claude]

@JuliaRegistrator register subdir=lib/RecursiveArrayToolsArrayPartitionAnyAll

[JuliaRegistrator]

Error while trying to register: Register Failed
@ChrisRackauckas-Claude, it looks like you are not a publicly listed member/owner in the parent organization (SciML).
If you are a member/owner, you will need to change your membership to public. See GitHub Help

[ChrisRackauckas-Claude]

@JuliaRegistrator register subdir=lib/RecursiveArrayToolsRaggedArrays

[JuliaRegistrator]

Error while trying to register: Register Failed
@ChrisRackauckas-Claude, it looks like you are not a publicly listed member/owner in the parent organization (SciML).
If you are a member/owner, you will need to change your membership to public. See GitHub Help

[ChrisRackauckas-Claude]

@JuliaRegistrator register subdir=lib/RecursiveArrayToolsShorthandConstructors

[JuliaRegistrator]

Error while trying to register: Register Failed
@ChrisRackauckas-Claude, it looks like you are not a publicly listed member/owner in the parent organization (SciML).
If you are a member/owner, you will need to change your membership to public. See GitHub Help

[JoshuaLampert]

Alright, I think this did it all. @JoshuaLampert and @mateuszbaran, do follow up if the ragged versions aren't matching the previous behavior. But I think we hit every requirement here using the sublibraries.

The ragged version does not match the previous behavior regarding the treatment of end citing the bot from above:

One difference from v3: end in the first dimension now uses the max inner array length (not per-column like v3's RaggedEnd), so r[end, i] will throw BoundsError if column i is shorter than the max. Use r[:, i] for safe column access, or index explicitly with length(r.u[i]).

What is the reason to let it throw a BoundsError if the column is shorter than the max? This gave the last element of that column (which is the intuitive behavior I would say) in the previous version.

[ChrisRackauckas]

Okay easiest way to square it was to just make sure all old code and tests was revived #559