Nested vector variables in ODE problem definitions

Specific Domains Modelling & Simulations
1

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
end

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?

2

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.

2 Likes
3

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

4

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

julia> recursive_unitless_bottom_eltype([xyzVec(),xyzVec(),xVec(1.0)])
ERROR: StackOverflowError:
Stacktrace:
 [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
    x::Float64
    y::Float64
    z::Float64
end

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.

5

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

1 Like
6

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?

7

@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