Closed JinraeKim closed 3 years ago
Sorry, it actually works! (but why...?)
My test code is:
function test3()
_myvec = @SLVector (:a, :b)
myvec = @SLVector (:x, :y)
dynamics = function (X, p, t)
dx = -X.x
dy = -X.y
myvec(dx, dy)
end
tspan = (0.0, 1.0)
x = _myvec(1.0, 2.0)
y = _myvec(1.0, 2.0)
X = myvec(x, y)
prob = ODEProblem(dynamics, X, tspan)
sol = @time solve(prob, Tsit5())
end
The result is:
julia> include("test/test_APIs.jl"); test3()
1.151163 seconds (3.47 M allocations: 226.099 MiB, 3.33% gc time)
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 5-element Vector{Float64}:
0.0
0.10001264591185442
0.34206322751658697
0.655359775820945
1.0
u: 5-element Vector{SLArray{Tuple{2}, SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)}, 1, 2, (:x, :y)}}:
2-element SLArray{Tuple{2}, SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)}, 1, 2, (:x, :y)} with indices SOneTo(2):
:x => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 1.0
:b => 2.0
:y => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 1.0
:b => 2.0
2-element SLArray{Tuple{2}, SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)}, 1, 2, (:x, :y)} with indices SOneTo(2):
:x => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.9048259756770471
:b => 1.8096519513540943
:y => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.9048259756770471
:b => 1.8096519513540943
2-element SLArray{Tuple{2}, SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)}, 1, 2, (:x, :y)} with indices SOneTo(2):
:x => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.7103033089908293
:b => 1.4206066179816585
:y => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.7103033089908293
:b => 1.4206066179816585
2-element SLArray{Tuple{2}, SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)}, 1, 2, (:x, :y)} with indices SOneTo(2):
:x => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.5192552947803908
:b => 1.0385105895607816
:y => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.5192552947803908
:b => 1.0385105895607816
2-element SLArray{Tuple{2}, SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)}, 1, 2, (:x, :y)} with indices SOneTo(2):
:x => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.3678795934784995
:b => 0.735759186956999
:y => 2-element SLArray{Tuple{2}, Float64, 1, 2, (:a, :b)} with indices SOneTo(2):
:a => 0.3678795934784995
:b => 0.735759186956999
But I still cannot find a way for "in-place" methods.
An example is:
function test4()
dynamics = function (dX, X, p, t)
dX.x = -X.x
dX.y = -X.y
end
tspan = (0.0, 1.0)
x = @LArray [1.0, 2.0] (:a, :b)
y = @LArray [1.0, 2.0] (:a, :b)
X = @LArray [x, y] (:x, :y)
prob = ODEProblem(dynamics, X, tspan)
sol = @time solve(prob, Tsit5())
end
which yields
julia> include("test/test_APIs.jl"); test4()
ERROR: MethodError: no method matching zero(::Type{LArray{Float64, 1, Vector{Float64}, (:a, :b)}})
Closest candidates are:
zero(::Union{Type{P}, P}) where P<:Dates.Period at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Dates/src/periods.jl:53
zero(::ForwardDiff.Partials) at /Users/jinrae/.julia/packages/ForwardDiff/m7cm5/src/partials.jl:39
zero(::SparseArrays.AbstractSparseArray) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/SparseArrays.jl:55
...
Stacktrace:
[1] zero(x::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)})
@ Base ./abstractarray.jl:1085
[2] alg_cache(alg::Tsit5, u::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, rate_prototype::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, uEltypeNoUnits::Type, uBottomEltypeNoUnits::Type, tTypeNoUnits::Type, uprev::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, uprev2::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, f::ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, t::Float64, dt::Float64, reltol::Float64, p::SciMLBase.NullParameters, calck::Bool, #unused#::Val{true})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yCczp/src/caches/low_order_rk_caches.jl:349
[3] __init(prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta1::Nothing, beta2::Nothing, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yCczp/src/solve.jl:293
[4] __init(prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, alg::Tsit5, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yCczp/src/solve.jl:66
[5] __solve(::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yCczp/src/solve.jl:4
[6] __solve(::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Tsit5)
@ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/yCczp/src/solve.jl:4
[7] solve_call(_prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5; merge_callbacks::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/krS36/src/solve.jl:61
[8] solve_call(_prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5)
@ DiffEqBase ~/.julia/packages/DiffEqBase/krS36/src/solve.jl:48
[9] solve_up(prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, p::SciMLBase.NullParameters, args::Tsit5; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/krS36/src/solve.jl:82
[10] solve_up(prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, p::SciMLBase.NullParameters, args::Tsit5)
@ DiffEqBase ~/.julia/packages/DiffEqBase/krS36/src/solve.jl:75
[11] solve(prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ DiffEqBase ~/.julia/packages/DiffEqBase/krS36/src/solve.jl:70
[12] solve(prob::ODEProblem{LArray{LArray{Float64, 1, Vector{Float64}, (:a, :b)}, 1, Vector{LArray{Float64, 1, Vector{Float64}, (:a, :b)}}, (:x, :y)}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, var"#51#52", LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5)
@ DiffEqBase ~/.julia/packages/DiffEqBase/krS36/src/solve.jl:68
[13] macro expansion
@ ./timing.jl:210 [inlined]
[14] test4()
@ Main ~/.julia/dev/NestedEnvironments/test/test_APIs.jl:76
[15] top-level scope
@ REPL[3]:1
ComponentArrays.jl seems an appropriate package for this. So I closed this issue.
I'm trying to find a way to use a "nested LabelledArray" for *DiffEq.
For example, the following code is what I wanted.
However, it does not work since
x
is not an Array of Numbers. The nested LabelledArray would be very helpful when dealing with dynamical equations of nested and complicated systems, e.g., multi-agent and interactive simulation.Any thoughts?