Compatibility with DifferentialEquations.jl

[stevengj]

Right now, this package is not compatible with the RecursiveArrayTools.jl package used by DifferentialEquations.jl.

e.g. here is a simple example from the DifferentialEquations tutorial, modified to use a TupleVec(dx, [dy, dz]) of a scalar and a vector:

using DifferentialEquations, TupleVectorSpaces

function lorenz(u, p, t)
    dx, (dy, dz) = Tuple(u)
    dx = 10.0 * (u[2] - u[1])
    dy = u[1] * (28.0 - u[3]) - u[2]
    dz = u[1] * u[2] - (8 / 3) * u[3]
    TupleVec(dx, [dy, dz])
end

u0 = TupleVec(1.0, [0.0, 0.0])
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
solve(prob, Tsit5());

This gives the error:

MethodError: no method matching ndims(::Type{TupleVec{Tuple{Float64, Vector{Float64}}}})
The function `ndims` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  ndims(::Type{Union{}}, Any...)
   @ Base abstractarray.jl:276
  ndims(::Type{<:AbstractVectorOfArray{T, N}}) where {T, N}
   @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/Y3i0V/src/vector_of_array.jl:595

Although it would technically be possible to support indexing into a TupleVec, i.e. to make it act like an AbstractVector, I'd really rather keep it more general and represent an element of an arbitrary normed vector space. I don't see any fundamental technical reason why this shouldn't be compatible with numerical ODE solvers (just as it's already compatible with numerical integration).

@ChrisRackauckas, is there a way to tell RecursiveArrayTools not to try to "look inside" a TupleVec, and instead treat it as an opaque mathematical object? As long as we're not trying to work in-place, at least for explicit solvers I don't see why you should need to look inside the vector.

[ChrisRackauckas]

I think things are mostly close. I'd add:

Base.ndims(u::Type{<:TupleVec}) = 1
Base.zero(u::TupleVec) = TupleVec(zero.(Tuple(u)))

But the main issue I'm seeing is:

julia> u0 = TupleVec(1.0, [0.0; 0.0])
TupleVec(1.0, [0.0, 0.0])

julia> u0 .+ u0
3-element Vector{Float64}:
 2.0
 0.0
 0.0

While it's generally not a requirement of broadcast, we do generally require that broadcast is type-preserving if it's a linear combination of the same types + scalars, i.e. that the vector space property is held and the linear combination is in the vector space.

[stevengj]

There's already a zero function. ndims seems a bit weird because it's not indexable … why is that required? But I can certainly add it.

Would it suffice to just define

 Broadcast.broadcastable(v::TupleVec) = Ref(v)

so that it is not treated as a container for broadcasting?

[stevengj]

Disabling broadcasting then causes DiffEq to be stuck at this line: https://github.com/SciML/OrdinaryDiffEq.jl/blob/0ea61dc630ea0d6159e70aa134bbb674ae13fa25/lib/OrdinaryDiffEqCore/src/solve.jl#L239

(Aside from the question of why there is a default nonzero abstol at all, is there a way to get DiffEq to not treat the errors elementwise and instead just treat the solution as an opaque normed object?)

If I override the default tolerances with solve(prob, Tsit5(), abstol=0, reltol=1e-6) (since the default reltol is also computed elementwise: https://github.com/SciML/OrdinaryDiffEq.jl/blob/0ea61dc630ea0d6159e70aa134bbb674ae13fa25/lib/OrdinaryDiffEqCore/src/solve.jl#L252), I then get an error at: https://github.com/SciML/OrdinaryDiffEq.jl/blob/0ea61dc630ea0d6159e70aa134bbb674ae13fa25/lib/OrdinaryDiffEqCore/src/solve.jl#L520

because zero(uBottomEltypeNoUnits) ends up calling zero(Any). I'm not sure what the right answer is here because I don't know what last_event_error represents.

[ChrisRackauckas]

using OrdinaryDiffEq, TupleVectorSpaces

using RecursiveArrayTools
function RecursiveArrayTools.recursive_unitless_bottom_eltype(a::Type{
        <:TupleVec,
})
    eltype(a)
end

function RecursiveArrayTools.recursive_unitless_eltype(a::Type{
        <:TupleVec,
})
    eltype(a)
end

Base.eltype(v::Type{<:TupleVec}) = @show promote_type(map(eltype, v.parameters[1].parameters)...)
Base.ndims(u::Type{<:TupleVec}) = 1
Base.zero(u::TupleVec) = TupleVec(zero.(Tuple(u)))
Base.one(u::TupleVec) = TupleVec(one.(Tuple(u)))

function lorenz(u, p, t)
    x, (y, z) = Tuple(u)
    dx = 10.0 * (y - x)
    dy = x * (28.0 - z) - y
    dz = x * y - (8 / 3) * z
    TupleVec(dx, [dy, dz])
end

u0 = TupleVec(1.0, [0.0, 0.0])
tspan = (0.0, 100.0)
prob = ODEProblem(lorenz, u0, tspan)
solve(prob, Tsit5(), abstol=1e-6, reltol=1e-3);

almost does it. But now the issue is that broadcast isn't defined:

Broadcast.broadcastable(v::TupleVec) = Ref(v)
u0 .* u0

that wasn't sufficient.

[stevengj]

I might just need to bite the bullet and make it indexable. I was trying to see how far I could get without that as it's not strictly required for an abstract vector space, but since we support iteration we might as well, and it should make it a lot easier to work with with tools that expect AbstractVector-like containers.

[ChrisRackauckas]

I think it's fine with just a recursive broadcast if the types are all matching.