From 5bf7a417bfc3957468696492289de937f4425e24 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 27 Apr 2026 13:19:46 +0200 Subject: [PATCH 1/6] add recursivefill! and recursivecopy! for ragged arrays --- .../src/RecursiveArrayToolsRaggedArrays.jl | 38 ++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl b/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl index 009cd044..df70a429 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl @@ -2,7 +2,8 @@ module RecursiveArrayToolsRaggedArrays import RecursiveArrayTools: RecursiveArrayTools, AbstractRaggedVectorOfArray, AbstractRaggedDiffEqArray, VectorOfArray, DiffEqArray, - AbstractVectorOfArray, AbstractDiffEqArray, AllObserved + AbstractVectorOfArray, AbstractDiffEqArray, AllObserved, + recursivefill!, recursivecopy! using SymbolicIndexingInterface using SymbolicIndexingInterface: ParameterTimeseriesCollection, ParameterIndexingProxy, ScalarSymbolic, ArraySymbolic, NotSymbolic, Timeseries, SymbolCache @@ -1725,4 +1726,39 @@ end # Re-export has_discretes and get_discretes for the non-ragged types has_discretes(::TT) where {TT <: AbstractDiffEqArray} = hasfield(TT, :discretes) +function recursivecopy!(b::AbstractRaggedVectorOfArray, a::AbstractRaggedVectorOfArray) + @inbounds for i in eachindex(b.u, a.u) + if ArrayInterface.ismutable(b.u[i]) || b.u[i] isa AbstractRaggedVectorOfArray + recursivecopy!(b.u[i], a.u[i]) + else + b.u[i] = copy(a.u[i]) + end + end + return b +end + +function recursivefill!( + b::AbstractRaggedVectorOfArray{T, N}, + a::T2 + ) where {T <: Union{Number, Bool}, T2 <: Union{Number, Bool}, N} + return fill!(b, a) +end + +function recursivefill!( + b::AbstractRaggedVectorOfArray{T, N}, + a::T2 + ) where {T <: StaticArraysCore.SArray, T2 <: Union{Number, Bool}, N} + @inbounds for arr in b.u, i in eachindex(arr) + arr[i] = map(_ -> a, arr[i]) + end + return b +end + +function recursivefill!(b::AbstractRaggedVectorOfArray{T, N}, a) where {T <: AbstractArray, N} + @inbounds for arr in b.u + recursivefill!(arr, a) + end + return b +end + end # module RecursiveArrayToolsRaggedArrays From 6609672ab7fa97ac814eccca7e4cd5c3c0819d66 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 27 Apr 2026 13:26:02 +0200 Subject: [PATCH 2/6] add tests --- .../test/runtests.jl | 30 +++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl index b3489c41..6519afb0 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl @@ -986,4 +986,34 @@ using Test @test rd isa RecursiveArrayTools.AbstractRaggedVectorOfArray @test !(rd isa AbstractArray) end + + @testset "recursivefill! for RaggedVectorOfArray" begin + # Bool argument — the pattern used by ODE solver cache initialisation + r = RaggedVectorOfArray([ones(2), ones(3)]) + recursivefill!(r, false) + @test r[:, 1] == [0.0, 0.0] + @test r[:, 2] == [0.0, 0.0, 0.0] + + # Numeric argument + r2 = RaggedVectorOfArray([zeros(2), zeros(3)]) + recursivefill!(r2, 1.0) + @test r2[:, 1] == [1.0, 1.0] + @test r2[:, 2] == [1.0, 1.0, 1.0] + + # Ragged sizes are preserved + @test length(r[:, 1]) == 2 + @test length(r[:, 2]) == 3 + end + + @testset "recursivecopy! for RaggedVectorOfArray" begin + src = RaggedVectorOfArray([ones(2), 2 * ones(3)]) + dst = RaggedVectorOfArray([zeros(2), zeros(3)]) + recursivecopy!(dst, src) + @test dst[:, 1] == [1.0, 1.0] + @test dst[:, 2] == [2.0, 2.0, 2.0] + + # Verify deep copy — modifying src must not affect dst + src[:, 1] .= 99.0 + @test dst[:, 1] == [1.0, 1.0] + end end From 08d34500412a8a540525cd628ba122d8badf03f7 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 27 Apr 2026 14:17:08 +0200 Subject: [PATCH 3/6] fix mapreduce --- .../src/RecursiveArrayToolsRaggedArrays.jl | 6 +++++- lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl | 11 +++++++++++ 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl b/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl index df70a429..384dd281 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl @@ -1521,7 +1521,11 @@ end Base.map(f, A::AbstractRaggedVectorOfArray) = map(f, A.u) function Base.mapreduce(f, op, A::AbstractRaggedVectorOfArray; kwargs...) - return mapreduce(f, op, view(A, ntuple(_ -> :, ndims(A))...); kwargs...) + # For full reduction (no kwargs): safely recurse over u to handle ragged inner shapes. + # The view-based approach uses size(A.u[1]) for all columns, which fails when inner + # arrays are themselves ragged with different column counts. + isempty(kwargs) || return mapreduce(f, op, view(A, ntuple(_ -> :, ndims(A))...); kwargs...) + return mapreduce(u -> mapreduce(f, op, u), op, A.u) end function Base.mapreduce( f, op, A::AbstractRaggedVectorOfArray{T, 1, <:AbstractVector{T}}; kwargs... diff --git a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl index 6519afb0..923cc63f 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl @@ -1016,4 +1016,15 @@ using Test src[:, 1] .= 99.0 @test dst[:, 1] == [1.0, 1.0] end + + @testset "mapreduce over nested ragged arrays" begin + # Outer array whose inner RaggedVoA elements have different column counts. + # mapreduce must recurse over A.u rather than building a fixed-shape view. + inner1 = RaggedVectorOfArray([ones(3), ones(3)]) # 2 columns + inner2 = RaggedVectorOfArray([ones(3), ones(3), ones(3)]) # 3 columns — ragged! + u = RaggedVectorOfArray([inner1, inner2]) + + @test mapreduce(identity, +, u) == 15.0 # (2+3)*3 + end + end From e1e0c075585b14aacc7a6357f5e92c09d9c3521b Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 27 Apr 2026 16:53:41 +0200 Subject: [PATCH 4/6] add SciMLBase extension --- .../Project.toml | 9 ++++++- ...siveArrayToolsRaggedArraysDiffEqBaseExt.jl | 27 +++++++++++++++++++ .../test/runtests.jl | 14 ++++++++++ 3 files changed, 49 insertions(+), 1 deletion(-) create mode 100644 lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl diff --git a/lib/RecursiveArrayToolsRaggedArrays/Project.toml b/lib/RecursiveArrayToolsRaggedArrays/Project.toml index 7a900689..f295266f 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/Project.toml +++ b/lib/RecursiveArrayToolsRaggedArrays/Project.toml @@ -10,6 +10,12 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" +[weakdeps] +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" + +[extensions] +RecursiveArrayToolsRaggedArraysDiffEqBaseExt = "DiffEqBase" + [compat] Adapt = "4" ArrayInterface = "7" @@ -20,9 +26,10 @@ SymbolicIndexingInterface = "0.3.35" julia = "1.10" [extras] +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["SparseArrays", "SymbolicIndexingInterface", "Test"] +test = ["DiffEqBase", "SparseArrays", "SymbolicIndexingInterface", "Test"] diff --git a/lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl b/lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl new file mode 100644 index 00000000..61495948 --- /dev/null +++ b/lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl @@ -0,0 +1,27 @@ +module RecursiveArrayToolsRaggedArraysDiffEqBaseExt + +import RecursiveArrayTools: AbstractRaggedVectorOfArray +import DiffEqBase + +# Mirror the AbstractVectorOfArray dispatch in DiffEqBase so that adaptive ODE +# solvers compute the correct RMS-normalized norm instead of the unnormalized +# Euclidean norm. Without these methods, ODE_DEFAULT_NORM falls through to +# `norm(u)` = sqrt(sum_abs2), which is sqrt(n_elements) times larger than the +# intended RMS norm, making the adaptive controller target a stricter tolerance +# than requested (abstol/reltol). + +function DiffEqBase.UNITLESS_ABS2(x::AbstractRaggedVectorOfArray) + return mapreduce(DiffEqBase.UNITLESS_ABS2, +, x.u; + init = zero(real(eltype(x)))) +end + +function DiffEqBase.recursive_length(u::AbstractRaggedVectorOfArray) + return sum(DiffEqBase.recursive_length, u.u; init = 0) +end + +function DiffEqBase.ODE_DEFAULT_NORM(u::AbstractRaggedVectorOfArray, _) + return Base.FastMath.sqrt_fast( + DiffEqBase.UNITLESS_ABS2(u) / max(DiffEqBase.recursive_length(u), 1)) +end + +end # module diff --git a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl index 923cc63f..514e4492 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl @@ -2,6 +2,7 @@ using RecursiveArrayTools, RecursiveArrayToolsRaggedArrays using RecursiveArrayToolsRaggedArrays: RaggedEnd, RaggedRange using SymbolicIndexingInterface using SymbolicIndexingInterface: SymbolCache +import DiffEqBase: ODE_DEFAULT_NORM, UNITLESS_ABS2, recursive_length using Test @testset "RecursiveArrayToolsRaggedArrays" begin @@ -1027,4 +1028,17 @@ using Test @test mapreduce(identity, +, u) == 15.0 # (2+3)*3 end + @testset "ODE_DEFAULT_NORM: RMS-normalised for RaggedVectorOfArray" begin + # Loading OrdinaryDiffEqTsit5 (which depends on DiffEqBase) triggers the weakdep + # extension, giving the correct RMS-normalised norm instead of the unnormalised + # Euclidean norm used by the generic fallback. + r = RaggedVectorOfArray([ones(3), ones(3)]) # 6 ones + @test UNITLESS_ABS2(r) ≈ 6.0 + @test recursive_length(r) == 6 + # RMS norm of 6 ones = sqrt(6/6) = 1 + @test ODE_DEFAULT_NORM(r, 0.0) ≈ 1.0 + # Unnormalised Euclidean norm would be sqrt(6) ≈ 2.449 — make sure we don't get that + @test ODE_DEFAULT_NORM(r, 0.0) < 2.0 + end + end From 48afea0bca34c4618d995a1c6eb97e66a9010516 Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Mon, 27 Apr 2026 17:19:21 +0200 Subject: [PATCH 5/6] fix inference issue on julia v1.10 --- .../src/RecursiveArrayToolsRaggedArrays.jl | 26 +++++++++++++++---- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl b/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl index 384dd281..c5bf97e4 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/src/RecursiveArrayToolsRaggedArrays.jl @@ -1520,12 +1520,28 @@ end Base.map(f, A::AbstractRaggedVectorOfArray) = map(f, A.u) -function Base.mapreduce(f, op, A::AbstractRaggedVectorOfArray; kwargs...) - # For full reduction (no kwargs): safely recurse over u to handle ragged inner shapes. - # The view-based approach uses size(A.u[1]) for all columns, which fails when inner - # arrays are themselves ragged with different column counts. +# Named functor used by the nested-ragged mapreduce to ensure type-stable dispatch. +struct _RaggedMapReduce{F, Op} + f::F + op::Op +end +@inline (w::_RaggedMapReduce)(u) = mapreduce(w.f, w.op, u) + +# When inner elements are themselves ragged, the view-based approach fails: view uses +# size(A.u[1]) for every column, causing BoundsErrors when inner shapes differ. +# We recurse element-by-element instead. Dispatching on the type of A.u (rather than +# using an if-check at runtime) keeps inference type-stable down to Julia 1.10. +function Base.mapreduce( + f, op, + A::AbstractRaggedVectorOfArray{T, N, <:AbstractVector{<:AbstractRaggedVectorOfArray}}; + kwargs... + ) where {T, N} isempty(kwargs) || return mapreduce(f, op, view(A, ntuple(_ -> :, ndims(A))...); kwargs...) - return mapreduce(u -> mapreduce(f, op, u), op, A.u) + return mapreduce(_RaggedMapReduce(f, op), op, A.u) +end + +function Base.mapreduce(f, op, A::AbstractRaggedVectorOfArray; kwargs...) + return mapreduce(f, op, view(A, ntuple(_ -> :, ndims(A))...); kwargs...) end function Base.mapreduce( f, op, A::AbstractRaggedVectorOfArray{T, 1, <:AbstractVector{T}}; kwargs... From 546b1abab99e4b1d1009f19afbdf6b9b0194d7aa Mon Sep 17 00:00:00 2001 From: Joshua Lampert Date: Tue, 28 Apr 2026 14:58:54 +0200 Subject: [PATCH 6/6] remove extension again --- .../Project.toml | 9 +------ ...siveArrayToolsRaggedArraysDiffEqBaseExt.jl | 27 ------------------- .../test/runtests.jl | 14 ---------- 3 files changed, 1 insertion(+), 49 deletions(-) delete mode 100644 lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl diff --git a/lib/RecursiveArrayToolsRaggedArrays/Project.toml b/lib/RecursiveArrayToolsRaggedArrays/Project.toml index f295266f..7a900689 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/Project.toml +++ b/lib/RecursiveArrayToolsRaggedArrays/Project.toml @@ -10,12 +10,6 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" -[weakdeps] -DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" - -[extensions] -RecursiveArrayToolsRaggedArraysDiffEqBaseExt = "DiffEqBase" - [compat] Adapt = "4" ArrayInterface = "7" @@ -26,10 +20,9 @@ SymbolicIndexingInterface = "0.3.35" julia = "1.10" [extras] -DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["DiffEqBase", "SparseArrays", "SymbolicIndexingInterface", "Test"] +test = ["SparseArrays", "SymbolicIndexingInterface", "Test"] diff --git a/lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl b/lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl deleted file mode 100644 index 61495948..00000000 --- a/lib/RecursiveArrayToolsRaggedArrays/ext/RecursiveArrayToolsRaggedArraysDiffEqBaseExt.jl +++ /dev/null @@ -1,27 +0,0 @@ -module RecursiveArrayToolsRaggedArraysDiffEqBaseExt - -import RecursiveArrayTools: AbstractRaggedVectorOfArray -import DiffEqBase - -# Mirror the AbstractVectorOfArray dispatch in DiffEqBase so that adaptive ODE -# solvers compute the correct RMS-normalized norm instead of the unnormalized -# Euclidean norm. Without these methods, ODE_DEFAULT_NORM falls through to -# `norm(u)` = sqrt(sum_abs2), which is sqrt(n_elements) times larger than the -# intended RMS norm, making the adaptive controller target a stricter tolerance -# than requested (abstol/reltol). - -function DiffEqBase.UNITLESS_ABS2(x::AbstractRaggedVectorOfArray) - return mapreduce(DiffEqBase.UNITLESS_ABS2, +, x.u; - init = zero(real(eltype(x)))) -end - -function DiffEqBase.recursive_length(u::AbstractRaggedVectorOfArray) - return sum(DiffEqBase.recursive_length, u.u; init = 0) -end - -function DiffEqBase.ODE_DEFAULT_NORM(u::AbstractRaggedVectorOfArray, _) - return Base.FastMath.sqrt_fast( - DiffEqBase.UNITLESS_ABS2(u) / max(DiffEqBase.recursive_length(u), 1)) -end - -end # module diff --git a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl index 514e4492..923cc63f 100644 --- a/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl +++ b/lib/RecursiveArrayToolsRaggedArrays/test/runtests.jl @@ -2,7 +2,6 @@ using RecursiveArrayTools, RecursiveArrayToolsRaggedArrays using RecursiveArrayToolsRaggedArrays: RaggedEnd, RaggedRange using SymbolicIndexingInterface using SymbolicIndexingInterface: SymbolCache -import DiffEqBase: ODE_DEFAULT_NORM, UNITLESS_ABS2, recursive_length using Test @testset "RecursiveArrayToolsRaggedArrays" begin @@ -1028,17 +1027,4 @@ using Test @test mapreduce(identity, +, u) == 15.0 # (2+3)*3 end - @testset "ODE_DEFAULT_NORM: RMS-normalised for RaggedVectorOfArray" begin - # Loading OrdinaryDiffEqTsit5 (which depends on DiffEqBase) triggers the weakdep - # extension, giving the correct RMS-normalised norm instead of the unnormalised - # Euclidean norm used by the generic fallback. - r = RaggedVectorOfArray([ones(3), ones(3)]) # 6 ones - @test UNITLESS_ABS2(r) ≈ 6.0 - @test recursive_length(r) == 6 - # RMS norm of 6 ones = sqrt(6/6) = 1 - @test ODE_DEFAULT_NORM(r, 0.0) ≈ 1.0 - # Unnormalised Euclidean norm would be sqrt(6) ≈ 2.449 — make sure we don't get that - @test ODE_DEFAULT_NORM(r, 0.0) < 2.0 - end - end