Support tuples?

[dlfivefifty]

Fairly regularly (in teaching) I need to turn a second order ODE to a first order one. I'm surprised tuples can't be used to do this, see below. I also don't understand the conceptional reason why: tuples do not support addition or multiplication but they do support broadcasting addition and multiplication, which is what the "time-stepping" algorithms should be using.

I realise StaticArrays gives a work around but this is less than ideal for pedagogical usages, and also practical reasons when one has mixed types (say first argument is complex and the second real).

julia> using DifferentialEquations, Plots

julia> V = x -> x^2/2  # Potential
#5 (generic function with 1 method)

julia> Vp = x -> x     # Force
#7 (generic function with 1 method)

julia> λ_0 = 2.3  # initial location
2.3

julia> v_0 = 1.2  # initial position
1.2

julia> λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)); reltol=1E-6);
ERROR: MethodError: no method matching vec(::Tuple{Float64,Float64})
Closest candidates are:
  vec(::Transpose{T,#s627} where #s627<:(AbstractArray{T,1} where T) where T) at /Users/sheehanolver/Projects/julia-1.3/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/adjtrans.jl:201
  vec(::AbstractSparseArray{Tv,Ti,1} where Ti where Tv) at /Users/sheehanolver/Projects/julia-1.3/usr/share/julia/stdlib/v1.3/SparseArrays/src/sparsevector.jl:913
  vec(::FillArrays.Ones{T,N,Axes} where Axes where N) where T at /Users/sheehanolver/.julia/packages/FillArrays/Aj0C4/src/fillalgebra.jl:3
  ...
Stacktrace:
 [1] _vec(::Tuple{Float64,Float64}) at /Users/sheehanolver/.julia/packages/DiffEqBase/avuk1/src/utils.jl:224
 [2] build_J_W(::Rodas4{0,false,DefaultLinSolve,DataType}, ::Tuple{Float64,Float64}, ::Tuple{Float64,Float64}, ::DiffEqBase.NullParameters, ::Float64, ::Float64, ::Function, ::Type, ::Val{false}) at /Users/sheehanolver/.julia/packages/OrdinaryDiffEq/nV9bA/src/derivative_utils.jl:534
 [3] alg_cache(::Rodas4{0,false,DefaultLinSolve,DataType}, ::Tuple{Float64,Float64}, ::Tuple{Float64,Float64}, ::Type, ::Type, ::Type, ::Tuple{Float64,Float64}, ::Tuple{Float64,Float64}, ::ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Float64, ::Float64, ::Float64, ::DiffEqBase.NullParameters, ::Bool, ::Val{false}) at /Users/sheehanolver/.julia/packages/OrdinaryDiffEq/nV9bA/src/caches/rosenbrock_caches.jl:389
 [4] (::OrdinaryDiffEq.var"#194#195"{Tuple{Float64,Float64},Tuple{Float64,Float64},DataType,DataType,DataType,Tuple{Float64,Float64},Tuple{Float64,Float64},ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64,Float64,DiffEqBase.NullParameters,Bool})(::Rodas4{0,false,DefaultLinSolve,DataType}) at /Users/sheehanolver/.julia/packages/OrdinaryDiffEq/nV9bA/src/caches/basic_caches.jl:15
 [5] map(::OrdinaryDiffEq.var"#194#195"{Tuple{Float64,Float64},Tuple{Float64,Float64},DataType,DataType,DataType,Tuple{Float64,Float64},Tuple{Float64,Float64},ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,Float64,Float64,DiffEqBase.NullParameters,Bool}, ::Tuple{Vern7,Rodas4{0,false,DefaultLinSolve,DataType}}) at ./tuple.jl:140
 [6] alg_cache(::CompositeAlgorithm{Tuple{Vern7,Rodas4{0,false,DefaultLinSolve,DataType}},AutoSwitch{Vern7,Rodas4{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}, ::Tuple{Float64,Float64}, ::Tuple{Float64,Float64}, ::Type, ::Type, ::Type, ::Tuple{Float64,Float64}, ::Tuple{Float64,Float64}, ::ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Float64, ::Float64, ::Float64, ::DiffEqBase.NullParameters, ::Bool, ::Val{false}) at /Users/sheehanolver/.julia/packages/OrdinaryDiffEq/nV9bA/src/caches/basic_caches.jl:14
 [7] #__init#328(::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Nothing, ::Bool, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Float64, ::Float64, ::Float64, ::Bool, ::Bool, ::Rational{Int64}, ::Nothing, ::Float64, ::Rational{Int64}, ::Int64, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Nothing, ::Nothing, ::Int64, ::typeof(DiffEqBase.ODE_DEFAULT_NORM), ::typeof(opnorm), ::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), ::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), ::Nothing, ::Bool, ::Bool, ::Bool, ::Bool, ::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol,Symbol},NamedTuple{(:default_set, :second_time),Tuple{Bool,Bool}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern7,Rodas4{0,false,DefaultLinSolve,DataType}},AutoSwitch{Vern7,Rodas4{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}, ::Array{Tuple{Float64,Float64},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at /Users/sheehanolver/.julia/packages/OrdinaryDiffEq/nV9bA/src/solve.jl:278
 [8] (::DiffEqBase.var"#kw##__init")(::NamedTuple{(:default_set, :second_time, :reltol),Tuple{Bool,Bool,Float64}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern7,Rodas4{0,false,DefaultLinSolve,DataType}},AutoSwitch{Vern7,Rodas4{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}, ::Array{Tuple{Float64,Float64},1}, ::Array{Float64,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./none:0 (repeats 5 times)
 [9] #__solve#327(::Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol,Symbol},NamedTuple{(:default_set, :second_time, :reltol),Tuple{Bool,Bool,Float64}}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern7,Rodas4{0,false,DefaultLinSolve,DataType}},AutoSwitch{Vern7,Rodas4{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}) at /Users/sheehanolver/.julia/packages/OrdinaryDiffEq/nV9bA/src/solve.jl:4
 [10] (::DiffEqBase.var"#kw##__solve")(::NamedTuple{(:default_set, :second_time, :reltol),Tuple{Bool,Bool,Float64}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::CompositeAlgorithm{Tuple{Vern7,Rodas4{0,false,DefaultLinSolve,DataType}},AutoSwitch{Vern7,Rodas4{0,false,DefaultLinSolve,DataType},Rational{Int64},Int64}}) at ./none:0
 [11] #__solve#1(::Bool, ::Base.Iterators.Pairs{Symbol,Real,Tuple{Symbol,Symbol},NamedTuple{(:second_time, :reltol),Tuple{Bool,Float64}}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Nothing) at /Users/sheehanolver/.julia/packages/DifferentialEquations/9ez1L/src/default_solve.jl:7
 [12] #__solve at ./none:0 [inlined]
 [13] #__solve#447 at /Users/sheehanolver/.julia/packages/DiffEqBase/avuk1/src/solve.jl:179 [inlined]
 [14] #__solve at ./none:0 [inlined]
 [15] #solve_call#442(::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:reltol,),Tuple{Float64}}}, ::typeof(DiffEqBase.solve_call), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at /Users/sheehanolver/.julia/packages/DiffEqBase/avuk1/src/solve.jl:38
 [16] #solve_call at ./none:0 [inlined]
 [17] #solve#443 at /Users/sheehanolver/.julia/packages/DiffEqBase/avuk1/src/solve.jl:61 [inlined]
 [18] (::DiffEqBase.var"#kw##solve")(::NamedTuple{(:reltol,),Tuple{Float64}}, ::typeof(solve), ::ODEProblem{Tuple{Float64,Float64},Tuple{Float64,Float64},false,DiffEqBase.NullParameters,ODEFunction{false,var"#9#10",UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}) at ./none:0
 [19] top-level scope at REPL[7]:1

[rvignolo]

Hi,

If the first one is real and the second one is complex, you can just use complex types.

Your example works if you just change it to the following:

julia> solve(ODEProblem(((λ,v),_,t) -> [v,-Vp(λ)], [λ_0,v_0], (0.0, 10.0)); reltol=1E-6)

which is funny, because the part (λ,v) in the function definition works. I suspect it does (λ,v) = u internally (in Julia).

I am trying to think about this but I do not see anything wrong. However, I see your point. I didn't know about the support of broadcasted functions in tuples.

Finally, try:

julia> function f(du, (λ,v),_,t) 
         du[1] = v
         du[2] = -Vp(λ)
       end
julia> solve(ODEProblem(f, [λ_0,v_0], (0.0, 10.0)); reltol=1E-6)

Thanks!

[ChrisRackauckas]

Tuples don't have vector semantics which is what causes the issue. Here's the amount of type piracy that you'd need to do in order to do this:

using DifferentialEquations, Plots

V = x -> x^2/2  # Potential
Vp = x -> x     # Force
λ_0 = 2.3  # initial location
v_0 = 1.2  # initial position

Base.zero(x::NTuple{N}) where {N} = ntuple(i->zero(x[i]),N)
Base.size(x::Tuple) = (length(x),)
Base.:*(x::Float64,y::Tuple) = x .* y
Base.:+(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] + y[i],N)

λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)), Tsit5(); reltol=1E-6)

Honestly, I like those definitions so I'd add them to DiffEqBase, but I am sure that would cause a frenzy haha.

[ChrisRackauckas]

However, to support stiff ODE solvers... you're asking for linear algebra on tuples. Not impossible, but 🤷‍♂ . Do you have a proposed fix for:

u = (λ_0,v_0)
J = false .* u .* u'

? That would be required to support stiff ODEs. What is the tuple equivalent to a matrix? What factorization code should it call? That all seems to be begging for static arrays.

tuples do not support addition or multiplication but they do support broadcasting addition and multiplication, which is what the "time-stepping" algorithms should be using.

If you broadcast static arrays then you lose a ton of performance right now, so the out-of-place functions do not broadcast.

[ChrisRackauckas]

We can also get auto-switch algorithms to work, if they never switch to the implicit method:

using StaticArrays
OrdinaryDiffEq._vec(x::Tuple) = SVector(x)
Base.adjoint(x::Tuple) = adjoint(SVector(x))
Base.:-(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] - y[i],N)

λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)); reltol=1E-6)

[ChrisRackauckas]

For future reference to myself, playing around:

using DifferentialEquations, Plots

V = x -> x^2/2  # Potential
Vp = x -> x     # Force
λ_0 = 2.3  # initial location
v_0 = 1.2  # initial position

Base.zero(x::NTuple{N}) where {N} = ntuple(i->zero(x[i]),N)
Base.size(x::Tuple) = (length(x),)
Base.:*(x::Float64,y::Tuple) = x .* y
Base.:+(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] + y[i],N)

λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)), Tsit5(); reltol=1E-6)

u = (λ_0,v_0)
J = false .* _vec(u) .* _vec(u)'

using StaticArrays
OrdinaryDiffEq._vec(x::Tuple) = SVector(x)
Base.adjoint(x::Tuple) = adjoint(SVector(x))
Base.:-(x::NTuple{N},y::NTuple{N}) where {N} = ntuple(i->x[i] - y[i],N)

using ForwardDiff
ForwardDiff.derivative(f,x::Tuple) = ForwardDiff.derivative(f,SVector(x))

λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)); reltol=1E-6)

λv = solve(ODEProblem(((λ,v),_,t) -> (v,-Vp(λ)), (λ_0,v_0), (0.0, 10.0)), TRBDF2(autodiff=false); reltol=1E-6)

[dlfivefifty]

If the first one is real and the second one is complex, you can just use complex types

That’s roughly 4x slower. Also that was just a model example.

A simpler approach might be to just create a VectorTuple type that has vector semantics but allows different types for the entries. Then a Tuple would auto convert to this type at the beginning. Linear algebra would require some work though... as my actual usage is just in lecture notes this is not urgent at all.

[ChrisRackauckas]

Using an automatic conversion to TupleVec, I can get it to be good for non-stiff ODEs. It's not robust though, and fails at the stiff ODE solver in some deep ways that really need StaticArray types of properties.

using OrdinaryDiffEq, StaticArrays, BenchmarkTools

struct TupleVec{T}
    tup::T
end
Base.:+(x::TupleVec,y::TupleVec) = TupleVec(map(+,x.tup,y.tup))
Base.:-(x::TupleVec,y::TupleVec) = TupleVec(map(-,x.tup,y.tup))
Base.:-(x::TupleVec) = TupleVec(map(-,x.tup))
Base.zero(x::TupleVec) = TupleVec(map(zero,x.tup))
Base.size(x::TupleVec) = (length(x.tup),)
Base.length(x::TupleVec) = length(x.tup)
Base.iterate(x::TupleVec,i) = iterate(x.tup,i)
Base.:*(x::Number,y::TupleVec) = TupleVec(x .* y.tup)
Base.:*(x::TupleVec,y::Number) = TupleVec(x.tup .* y)
Base.:/(x::TupleVec,y::Number) = TupleVec(x.tup ./ y)
Base.eltype(x::TupleVec) = promote_type(map(typeof,x.tup)...)
RecursiveArrayTools.recursive_unitless_bottom_eltype(x::TupleVec) = promote_type(map(x->typeof(one(x)),x.tup)...)
ArrayInterfaceCore.zeromatrix(x::TupleVec) = false .* SVector(x.tup) .* SVector(x.tup)'
Base.vec(x::TupleVec) = x

(f::SciMLBase.ODEFunction)(u::TupleVec,p,t) = TupleVec(f.f(u.tup,p,t))

DiffEqBase.function get_concrete_u0(prob, isadapt, t0, kwargs)
  if eval_u0(prob.u0)
    u0 = prob.u0(prob.p, t0)
  elseif haskey(kwargs, :u0)
    u0 = kwargs[:u0]
  else
    u0 = prob.u0
  end

  if u0 isa Tuple
    u0 = SciMLBase.TupleVec(prob.u0)
  end

  isadapt && eltype(u0) <: Integer && (u0 = float.(u0))

  _u0 = handle_distribution_u0(u0)

  if isinplace(prob) && (_u0 isa Number || _u0 isa SArray || _u0 isa SciMLBase.TupleVec)
    throw(IncompatibleInitialConditionError())
  end

  _u0
end


V(x) = x^2/2  # Potential
Vp(x) = x     # Force
λ_0 = 2.3  # initial location
v_0 = 1.2  # initial position
function ftup((λ,v),_,t)
    (v,-Vp(λ))
end

function fsa(u,_,t)
    SA[u[2],-Vp(u[1])]
end

@benchmark solve(ODEProblem(ftup, (λ_0,v_0), (0.0, 10.0)), Tsit5(); reltol=1E-6)
@benchmark solve(ODEProblem(fsa, SA[λ_0,v_0], (0.0, 10.0)), Tsit5(); reltol=1E-6)

using ForwardDiff
ForwardDiff.derivative(f,x::SciMLBase.TupleVec) = ForwardDiff.jacobian((args...)->SVector(f(args...)),SVector(x.tup))
solve(ODEProblem(ftup, (λ_0,v_0), (0.0, 10.0)), TRBDF2(); reltol=1E-6)

[ChrisRackauckas]

For now, I added an informative error:

julia> solve(prob,Tsit5())
ERROR: Tuple type used as a state. Since a tuple does not have vector
properties, it will not work as a state type in equation solvers.
Instead, change your equation from using tuple constructors `()`
to static array constructors `SA[]`. For example, change:

function ftup((a,b),p,t)
  return b,-a
end
u0 = (1.0,2.0)
tspan = (0.0,1.0)
ODEProblem(ftup,u0,tspan)

to:

using StaticArrays
function fsa(u,p,t)
    SA[u[2],u[1]]
end
u0 = SA[1.0,2.0]
tspan = (0.0,1.0)
ODEProblem(ftup,u0,tspan)

This will be safer and fast for small ODEs. For more information, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Further-Optimizations-of-Small-Non-Stiff-ODEs-with-StaticArrays   


Stacktrace:
 [1] get_concrete_u0(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool, t0::Float64, kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:u0, :p), Tuple{Tuple{Float64, Float32}, SciMLBase.NullParameters}}})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:907
 [2] get_concrete_problem(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool; kwargs::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:u0, :p), Tuple{Tuple{Float64, Float32}, SciMLBase.NullParameters}}})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:801
 [3] solve_up(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Tuple{Float64, Float32}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:722
 [4] solve_up(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Tuple{Float64, Float32}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:713
 [5] solve(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:708
 [6] solve(prob::ODEProblem{Tuple{Float64, Float32}, Tuple{Float64, Float64}, false, SciMLBase.NullParameters, ODEFunction{false, typeof(f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False})
   @ DiffEqBase C:\Users\accou\.julia\dev\DiffEqBase\src\solve.jl:703
 [7] top-level scope
   @ REPL[3]:1

SciML/DiffEqBase.jl#767

[ChrisRackauckas]

Closing as we now have good error messages on this.

[stevengj]

It's pretty frustrating not to have some kind of TupleVec type that supports this, because tuples (unlike svectors) support heterogeneous elements.

I just ran into a problem where I want to co-evolve ODEs for a 4-component vector and a 4x4 matrix; not clear how to do this without reshaping, but a TupleVec would make it a snap.

[dlfivefifty]

If the eltypes are the same can’t you use ComponentArrays.jl?

[stevengj]

My problem with ComponentArrays.jl is that it converts StaticArrays to Array:

julia> dump(ComponentArray(v=SA[1,2,3]))
ComponentVector{Int64, Vector{Int64}, Tuple{Axis{(v = 1:3,)}}}
  data: Array{Int64}((3,)) [1, 2, 3]
  axes: Tuple{Axis{(v = 1:3,)}}
    1: Axis{(v = 1:3,)} Axis(v = 1:3,)

which really hurts performance.

I think there is a place in the ecosystem for a type which represents a direct sum of arbitrary vector spaces, i.e. a tuple of vectors, that leaves them as-is rather than trying to squash them into a single flat Array.

(Obviously, such a type would not work with code that really requires an array, but for quadrature and ODEs often all you really need is a normed vector space.)