Nested vector variables in ODE problem definitions

I am trying to define an ODE in which the state variables are nested in arrays. Here is a toy example:

using DifferentialEquations

function f!(du, u, p, t)
    du[1] = u[1] - u[2][1]
    du[2] .= -5 * u[2]
    return nothing

u₀ = [1.0, [2.0, 3.0]]
tspan = (0.0, 10.0)

prob = ODEProblem(f!, u₀, tspan)

sol = solve(prob, Tsit5())

I assume this type of thing is supported since the solve function makes a call to the RecursiveArrayTools function recursive_unitless_bottom_eltype to get the base element type of the nested array. But here it throws a stack overflow error because the methods of recursive_unitless_bottom_eltype look like:

recursive_unitless_bottom_eltype(a) = recursive_unitless_bottom_eltype(typeof(a))
recursive_unitless_bottom_eltype(a::Type{T}) where T = recursive_unitless_bottom_eltype(eltype(a))
recursive_unitless_bottom_eltype(a::Type{T}) where {T<:AbstractArray} = recursive_unitless_bottom_eltype(eltype(a))
recursive_unitless_bottom_eltype(a::Type{T}) where {T<:Number} = eltype(a) == Number ? Float64 : typeof(one(eltype(a)))

Since the typeof(u₀) is Array{Any,1}, the second pass goes into the second method. Since eltype(Array{Any,1}) is Any and eltype(Any) is Any, we recurse indefinitely.

Is my use case supported or am I doing something that just shouldn’t be done?

This is one thing I haven’t been able to solve without a type because you cannot declare broadcast as recursive. u₀ = VectorOfArray([[1.0], [2.0, 3.0]]) works though, and so does ArrayPartition.


Oh, cool, I’ll have to give that a try. Thanks!

This is the same problem I am having on the Slack discussion.

julia> recursive_unitless_bottom_eltype([xyzVec(),xyzVec(),xVec(1.0)])
ERROR: StackOverflowError:
 [1] recursive_unitless_bottom_eltype(::Type{Any}) at C:\Users\bakerar\.julia\packages\RecursiveArrayTools\DBKJE\src\utils.jl:86 (repeats 65195 times)
 [2] recursive_unitless_bottom_eltype(::Type{Array{xyzVec,1}}) at C:\Users\bakerar\.julia\packages\RecursiveArrayTools\DBKJE\src\utils.jl:87
 [3] recursive_unitless_bottom_eltype(::Array{xyzVec,1}) at C:\Users\bakerar\.julia\packages\RecursiveArrayTools\DBKJE\src\utils.jl:85
 [4] top-level scope at none:0```

My composite type is simple:
struct xyzVec <:AbstractVec

I have lots of math for vectors defined…
and the thing that is causing the problem is an ill defined eltype function which I can’t find a good example of.

Base.iterate(v::AbstractVec, state=1) = state > length(v) ?
                                        nothing :
                                        (v[state], state+1) #Return the positional number

Base.getindex(v::AbstractVec, i1::Int64) = toTuple(v)[i1]

Base.eltype(::xyzVec) = Float64
Base.eltype(::xyVec) = Float64
Base.eltype(::xyVec) = Float64
Base.eltype(::AbstractVec) = Float64

I think I need my syntax corrected on the eltype and it will all work out, but so far no luck.

ComponentArrays.jl is basically what came out of this topic. Check it out and see if it fits your needs.

1 Like

ComponentArrays does work for my case, but VectorOfArray integrators are not resize!()able or deleteat!()able within callbacks (no methods). Any way to have a nested structure that allows for easily changing dimensionality?

@levasco Yep, MultiScaleArrays.jl was created specifically for this purpose (assuming the resizing and deleting you’re doing is of whole branches of the tree)

1 Like