Nested LabelledArray for DiffEq

[JinraeKim]

I'm trying to find a way to use a "nested LabelledArray" for *DiffEq.

For example, the following code is what I wanted.

    xa = @LArray [1, 2] (:x, :y)
    xb = @LArray [3, 4] (:z, :w)
    x = @LArray [xa, xb] (:a, :b)
    dynamics = function (dx, x, p, t)
        dx.a = -x.a
        dx.b = -x.b
    end
    tspan = (0.0, 1.0)
    prob = ODEProblem(dynamics, x, tspan)
    sol = @time solve(prob, Tsit5())

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?

[JinraeKim]

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

[JinraeKim]

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

[JinraeKim]

ComponentArrays.jl seems an appropriate package for this. So I closed this issue.