Support tuples? #566
Comments
|
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! |
|
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. |
|
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.
If you broadcast static arrays then you lose a ton of performance right now, so the out-of-place functions do not broadcast. |
|
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) |
|
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) |
That’s roughly 4x slower. Also that was just a model example. A simpler approach might be to just create a |
|
Using an automatic conversion to 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) |
|
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 |
|
Closing as we now have good error messages on this. |
|
It's pretty frustrating not to have some kind of 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 |
|
If the eltypes are the same can’t you use ComponentArrays.jl? |
|
My problem with ComponentArrays.jl is that it converts StaticArrays to 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 (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.) |
|
Yes I agree 100% (hence the creation of this issue!). But in your case since the eltypes are the same it seems like it would be very quick to whip up your own type: struct SVectorAndSMatrix{N,T} <: AbstractVector{T}
vec::SVector{N,T}
mat::SMatrix{(N,N),T}
end
size(A:: SVectorAndSMatrix{N}) where N = (N+N^2,)
getindex(A:: SVectorAndSMatrix{N}, k::Int) where N = k ≤ N ? A.vec[k] : A.mat[k-N] # don't actually use in your time stepper just to implement interface |
|
I created a TupleVecs.jl package to represent arbitrary direct sums of vector spaces (with norms and inner products if available). So far it works with numerical integration, differentiation, interpolation, and extrapolation, but it's not working yet with DifferentialEquations: JuliaMath/TupleVectorSpaces.jl#1 Mathematically, I'd really like to see how far I can get without implementing |
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).
The text was updated successfully, but these errors were encountered: