BREAKING: Make AbstractVectorOfArray <: AbstractArray

Project.toml

@@ -1,7 +1,7 @@
 name = "RecursiveArrayTools"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
+version = "4.0.0"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "3.51.0"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
@@ -15,6 +15,7 @@ StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
 [weakdeps]
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
@@ -29,6 +30,7 @@ Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extensions]
+RecursiveArrayToolsCUDAExt = "CUDA"
 RecursiveArrayToolsFastBroadcastExt = "FastBroadcast"
 RecursiveArrayToolsForwardDiffExt = "ForwardDiff"
 RecursiveArrayToolsKernelAbstractionsExt = "KernelAbstractions"
@@ -46,6 +48,7 @@ RecursiveArrayToolsZygoteExt = "Zygote"
 Adapt = "4"
 Aqua = "0.8"
 ArrayInterface = "7.16"
+CUDA = "5"
 DocStringExtensions = "0.9.3"
 FastBroadcast = "1.1"
 ForwardDiff = "0.10.38, 1"
@@ -61,7 +64,6 @@ Random = "1"
 RecipesBase = "1.3.4"
 ReverseDiff = "1.15"
 SafeTestsets = "0.1"
-SciMLBase = "2.103"
 SparseArrays = "1.10"
 StaticArrays = "1.6"
 StaticArraysCore = "1.4.2"
@@ -86,7 +88,6 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -97,4 +98,4 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "FastBroadcast", "ForwardDiff", "KernelAbstractions", "Measurements", "NLsolve", "Pkg", "Random", "SafeTestsets", "SciMLBase", "SparseArrays", "StaticArrays", "Statistics", "StructArrays", "Tables", "Test", "Unitful", "Zygote"]
+test = ["Aqua", "FastBroadcast", "ForwardDiff", "KernelAbstractions", "Measurements", "NLsolve", "Pkg", "Random", "SafeTestsets", "SparseArrays", "StaticArrays", "Statistics", "StructArrays", "Tables", "Test", "Unitful", "Zygote"]

docs/Project.toml

@@ -4,4 +4,4 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 
 [compat]
 Documenter = "1.3"
-RecursiveArrayTools = "3"
+RecursiveArrayTools = "4"

docs/pages.jl

@@ -2,7 +2,9 @@
 
 pages = [
     "Home" => "index.md",
+    "breaking_changes_v4.md",
     "AbstractVectorOfArrayInterface.md",
     "array_types.md",
+    "plotting.md",
     "recursive_array_functions.md",
 ]

docs/src/assets/Project.toml

@@ -0,0 +1,7 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+
+[compat]
+Documenter = "1.3"
+RecursiveArrayTools = "3"

docs/src/breaking_changes_v4.md

@@ -0,0 +1,123 @@
+# Breaking Changes in v4.0: AbstractArray Interface
+
+## Summary
+
+`AbstractVectorOfArray{T, N, A}` now subtypes `AbstractArray{T, N}`. This means
+all `VectorOfArray` and `DiffEqArray` objects are proper Julia `AbstractArray`s,
+and all standard `AbstractArray` operations work out of the box, including linear
+algebra, broadcasting with plain arrays, and generic algorithms.
+
+## Key Changes
+
+### Linear Indexing
+
+Previously, `A[i]` returned the `i`th inner array (`A.u[i]`). Now, `A[i]` returns
+the `i`th element in column-major linear order, matching standard Julia `AbstractArray`
+behavior.
+
+```julia
+A = VectorOfArray([[1, 2], [3, 4]])
+# Old: A[1] == [1, 2]  (first inner array)
+# New: A[1] == 1        (first element, column-major)
+# To access inner arrays: A.u[1] or A[:, 1]
+```
+
+### Size and Ragged Arrays
+
+For ragged arrays (inner arrays of different sizes), `size(A)` now reports the
+**maximum** size in each dimension. Out-of-bounds elements are treated as zero
+(sparse representation):
+
+```julia
+A = VectorOfArray([[1, 2], [3, 4, 5]])
+size(A)    # (3, 2) — max inner length is 3
+A[3, 1]    # 0      — implicit zero (inner array 1 has only 2 elements)
+A[3, 2]    # 5      — actual stored value
+Array(A)   # [1 3; 2 4; 0 5] — zero-padded dense array
+```
+
+This means ragged `VectorOfArray`s can be used directly with linear algebra
+operations, treating the data as a rectangular matrix with zero padding.
+
+### Iteration
+
+Iteration now goes over scalar elements in column-major order, matching
+`AbstractArray` behavior:
+
+```julia
+A = VectorOfArray([[1, 2], [3, 4]])
+collect(A)  # [1 3; 2 4] — 2x2 matrix
+# To iterate over inner arrays: for u in A.u ... end
+```
+
+### `length`
+
+`length(A)` now returns `prod(size(A))` (total number of elements including
+ragged zeros), not the number of inner arrays. Use `length(A.u)` for the number
+of inner arrays.
+
+### `map`
+
+`map(f, A)` now maps over individual elements, not inner arrays. Use
+`map(f, A.u)` to map over inner arrays.
+
+### `first` / `last`
+
+`first(A)` and `last(A)` return the first/last scalar element, not the first/last
+inner array. Use `first(A.u)` / `last(A.u)` for inner arrays.
+
+### `eachindex`
+
+`eachindex(A)` returns `CartesianIndices(size(A))` for the full rectangular shape,
+not indices into `A.u`.
+
+## Migration Guide
+
+| Old Code | New Code |
+|----------|----------|
+| `A[i]` (get inner array) | `A.u[i]` or `A[:, i]` |
+| `length(A)` (number of arrays) | `length(A.u)` |
+| `for elem in A` (iterate columns) | `for elem in A.u` |
+| `first(A)` (first inner array) | `first(A.u)` |
+| `map(f, A)` (map over columns) | `map(f, A.u)` |
+| `A == vec_of_vecs` | `A.u == vec_of_vecs` |
+
+## Interpolation Interface on DiffEqArray
+
+`DiffEqArray` now has `interp` and `dense` fields for interpolation support:
+
+```julia
+# Create with interpolation
+da = DiffEqArray(u, t, p, sys; interp = my_interp, dense = true)
+
+# Callable syntax
+da(0.5)                    # interpolate at t=0.5
+da(0.5, Val{1})            # first derivative at t=0.5
+da([0.1, 0.5, 0.9])       # interpolate at multiple times
+da(0.5; idxs = 1)          # interpolate single component
+da(0.5; idxs = [1, 2])    # interpolate subset of components
+```
+
+The interpolation object must be callable as `interp(t, idxs, deriv, p, continuity)`,
+matching the protocol used by SciMLBase's `LinearInterpolation`, `HermiteInterpolation`,
+and `ConstantInterpolation`.
+
+When `dense = true` and `interp` is provided, `plot(da)` automatically generates
+dense interpolated output instead of plotting only the saved time points.
+
+## Ragged Arrays Sublibrary
+
+`RaggedVectorOfArray` and `RaggedDiffEqArray` are available via:
+
+```julia
+using RecursiveArrayToolsRaggedArrays
+```
+
+These types preserve the true ragged structure without zero-padding, and do **not**
+subtype `AbstractArray`. See the `RecursiveArrayToolsRaggedArrays` subpackage for details.
+
+## Zygote Compatibility
+
+Some Zygote adjoint rules need updating for the new `AbstractArray` subtyping.
+ForwardDiff continues to work correctly. Zygote support will be updated in a
+follow-up release.

docs/src/plotting.md

@@ -0,0 +1,154 @@
+# Plotting
+
+`AbstractVectorOfArray` and `AbstractDiffEqArray` types include
+[Plots.jl](https://github.com/JuliaPlots/Plots.jl) recipes so they can be
+visualised directly with `plot(A)`.
+
+## VectorOfArray
+
+A `VectorOfArray` plots as a matrix where each inner array is a column:
+
+```julia
+using RecursiveArrayTools, Plots
+A = VectorOfArray([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+plot(A)   # 3 series (components) vs column index
+```
+
+## DiffEqArray
+
+A `DiffEqArray` plots as component time series against `A.t`:
+
+```julia
+u = [[sin(t), cos(t)] for t in 0:0.1:2pi]
+t = collect(0:0.1:2pi)
+A = DiffEqArray(u, t)
+plot(A)   # plots sin and cos vs t
+```
+
+If the `DiffEqArray` carries a symbolic system (via `variables` and
+`independent_variables` keyword arguments), the axis labels and series names
+are set automatically from the symbol names.
+
+## Dense (Interpolated) Plotting
+
+When a `DiffEqArray` has an interpolation object in its `interp` field and
+`dense = true`, calling `plot(A)` generates a smooth curve by evaluating the
+interpolation at many points rather than connecting only the saved time steps.
+
+```julia
+using SciMLBase: LinearInterpolation
+
+u = [[1.0, 0.0], [0.0, 1.0], [-1.0, 0.0], [0.0, -1.0]]
+t = [0.0, 1.0, 2.0, 3.0]
+interp = LinearInterpolation(t, u)
+A = DiffEqArray(u, t; interp = interp, dense = true)
+plot(A)   # smooth interpolated curve with 1000+ points
+```
+
+### Plot Recipe Keyword Arguments
+
+The `AbstractDiffEqArray` recipe accepts the following keyword arguments,
+which can be passed directly to `plot`:
+
+| Keyword | Type | Default | Description |
+|---------|------|---------|-------------|
+| `denseplot` | `Bool` | `A.dense && A.interp !== nothing` | Use dense interpolation for smooth curves. Set `false` to show only saved points. |
+| `plotdensity` | `Int` | `max(1000, 10 * length(A.u))` | Number of evenly-spaced points to evaluate when `denseplot = true`. |
+| `tspan` | `Tuple` or `nothing` | `nothing` | Restrict the time window. E.g. `tspan = (0.0, 5.0)`. |
+| `idxs` | varies | `nothing` | Select which components to plot (see below). |
+
+Example:
+
+```julia
+plot(A; denseplot = true, plotdensity = 5000, tspan = (0.0, 2.0))
+```
+
+## Selecting Variables with `idxs`
+
+The `idxs` keyword controls which variables appear in the plot. It supports
+several formats:
+
+### Single component
+
+```julia
+plot(A; idxs = 1)      # plot component 1 vs time
+plot(A; idxs = 2)      # plot component 2 vs time
+```
+
+### Multiple components
+
+```julia
+plot(A; idxs = [1, 3, 5])   # plot components 1, 3, 5 vs time
+```
+
+### Phase-space plots (component vs component)
+
+Use a tuple where index `0` represents the independent variable (time):
+
+```julia
+plot(A; idxs = (1, 2))      # component 1 vs component 2
+plot(A; idxs = (0, 1))      # time vs component 1 (same as default)
+plot(A; idxs = (1, 2, 3))   # 3D plot of components 1, 2, 3
+```
+
+### Symbolic indexing
+
+When the `DiffEqArray` carries a symbolic system, variables can be referenced
+by symbol:
+
+```julia
+A = DiffEqArray(u, t; variables = [:x, :y], independent_variables = [:t])
+plot(A; idxs = :x)           # plot x vs time
+plot(A; idxs = [:x, :y])     # plot both
+plot(A; idxs = (:x, :y))     # phase plot of x vs y
+```
+
+### Custom transformations
+
+A function can be applied to the plotted values:
+
+```julia
+plot(A; idxs = (norm, 0, 1, 2))   # plot norm(u1, u2) vs time
+```
+
+The tuple format is `(f, xvar, yvar)` or `(f, xvar, yvar, zvar)` where `f`
+is applied element-wise.
+
+## Callable Interface
+
+Any `AbstractDiffEqArray` with an `interp` field supports callable syntax for
+interpolation, independent of plotting:
+
+```julia
+A(0.5)                    # interpolate all components at t=0.5
+A(0.5; idxs = 1)          # interpolate component 1 at t=0.5
+A([0.1, 0.5, 0.9])       # interpolate at multiple times (returns DiffEqArray)
+A(0.5, Val{1})            # first derivative at t=0.5
+A(0.5; continuity = :right)  # right-continuity at discontinuities
+```
+
+The interpolation object must be callable as
+`interp(t, idxs, deriv, p, continuity)`, matching the protocol used by
+SciMLBase's `LinearInterpolation`, `HermiteInterpolation`, and
+`ConstantInterpolation`.
+
+When no interpolation is available (`interp === nothing`), calling `A(t)`
+throws an error.
+
+## ODE Solution Plotting
+
+ODE solutions from DifferentialEquations.jl are subtypes of
+`AbstractDiffEqArray` and inherit all of the above functionality, plus
+additional features:
+
+- **Automatic dense plotting**: `denseplot` defaults to `true` when the solver
+  provides dense output.
+- **Analytic solution overlay**: `plot(sol; plot_analytic = true)` overlays the
+  exact solution if the problem defines one.
+- **Discrete variables**: Time-varying parameters are plotted as step functions
+  with dashed lines and markers.
+- **Symbolic indexing**: Full symbolic indexing through ModelingToolkit is
+  supported, including observed (derived) variables.
+
+These advanced features are defined in SciMLBase and activate automatically
+when plotting solution objects.

ext/RecursiveArrayToolsCUDAExt.jl

@@ -0,0 +1,13 @@
+module RecursiveArrayToolsCUDAExt
+
+using RecursiveArrayTools: AbstractVectorOfArray
+import CUDA: CuArray
+
+# Disambiguate CuArray(::AbstractVectorOfArray) vs CuArray(::AbstractArray{T,N}) from CUDA.jl.
+# This is the exact signature Julia's ambiguity error requests.
+# Uses stack to stay on GPU (avoids GPU→CPU→GPU round-trip).
+function CuArray(VA::AbstractVectorOfArray{T, N}) where {T, N}
+    return CuArray{T, N}(stack(VA.u))
+end
+
+end