SciML / RecursiveArrayTools.jl

Tools for easily handling objects like arrays of arrays and deeper nestings in scientific machine learning (SciML) and other applications
https://docs.sciml.ai/RecursiveArrayTools/stable/
Other
212 stars 57 forks source link

Getting the `sol` of an observed variable (using SymbolicIndexingInterface) of a trivial ODE (with no reduced states) fails #365

Closed avinashresearch1 closed 6 months ago

avinashresearch1 commented 6 months ago

If I have a trivial ODE after structural simplification (no states), only observed vars, assuming I want to solve such a problem, it appears getting sol[] breaks e.g.,

MWE:

using ModelingToolkit, OrdinaryDiffEq, Test, SciMLBase
using ModelingToolkit: t_nounits as t, D_nounits as D

vars = @variables begin
    x(t)
    y(t)
end

eqs = [x ~ 1,
    y ~ x]

@named sys = ODESystem(eqs, t, vars, [])
sysRed = structural_simplify(sys)

unknowns(sysRed)

observed(sysRed)

prob = ODEProblem(sysRed, [], (0.0, 1.0))
sol = solve(prob, Rodas5())

sol[x]

gives the following stacktrace

ERROR: BoundsError: attempt to access Tuple{Vector{Float64}, Float64} at index [3]
Stacktrace:
  [1] getindex(t::Tuple, i::Int64)
    @ Base ./tuple.jl:31
  [2] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:162 [inlined]
  [3] macro expansion
    @ ./none:0 [inlined]
  [4] generated_callfunc
    @ ./none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
  [6] (::ModelingToolkit.var"#699#712"{…})(u::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/gIWqH/src/systems/diffeqs/abstractodesystem.jl:463
  [7] _broadcast_getindex_evalf
    @ ./broadcast.jl:709 [inlined]
  [8] _broadcast_getindex
    @ ./broadcast.jl:682 [inlined]
  [9] getindex
    @ ./broadcast.jl:636 [inlined]
 [10] copy
    @ ./broadcast.jl:942 [inlined]
 [11] materialize(bc::Base.Broadcast.Broadcasted{…})
    @ Base.Broadcast ./broadcast.jl:903
 [12] _getindex
    @ ~/.julia/packages/RecursiveArrayTools/Vclmo/src/vector_of_array.jl:367 [inlined]
 [13] getindex(::ODESolution{…}, ::Num)
    @ RecursiveArrayTools ~/.julia/packages/RecursiveArrayTools/Vclmo/src/vector_of_array.jl:433
 [14] top-level scope
    @ ~/Desktop/JULIAHUB/JULIASIM/HVAC.jl/tmp/mtkv9_rat_bug.jl:22
Some type information was truncated. Use `show(err)` to see complete types.

It appears to be related to the broadcast of observed here: https://github.com/SciML/RecursiveArrayTools.jl/blob/e78cb7443ebb68d910606ab03b2b915316eb6e98/src/vector_of_array.jl#L367

so:

julia> sol
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 1.0
u: 2-element Vector{Vector{Float64}}:
 []
 []

julia> RecursiveArrayTools.observed(sol, x)(sol.u, (RecursiveArrayTools.parameter_values(sol),), sol.t)
1.0

julia> RecursiveArrayTools.observed(sol, x).(sol.u, (RecursiveArrayTools.parameter_values(sol),), sol.t)
ERROR: BoundsError: attempt to access Tuple{Vector{Float64}, Float64} at index [3]
Stacktrace:
  [1] getindex(t::Tuple, i::Int64)
    @ Base ./tuple.jl:31
  [2] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:162 [inlined]
  [3] macro expansion
    @ ./none:0 [inlined]
  [4] generated_callfunc
    @ ./none:0 [inlined]
  [5] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{…})(::Vector{…}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/Yo8zx/src/RuntimeGeneratedFunctions.jl:150
  [6] (::ModelingToolkit.var"#699#712"{…})(u::Vector{…}, p::ModelingToolkit.MTKParameters{…}, t::Float64)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/gIWqH/src/systems/diffeqs/abstractodesystem.jl:463
  [7] _broadcast_getindex_evalf
    @ ./broadcast.jl:709 [inlined]
  [8] _broadcast_getindex
    @ ./broadcast.jl:682 [inlined]
  [9] getindex
    @ ./broadcast.jl:636 [inlined]
 [10] copy
    @ ./broadcast.jl:942 [inlined]
 [11] materialize(bc::Base.Broadcast.Broadcasted{…})
    @ Base.Broadcast ./broadcast.jl:903
 [12] top-level scope
    @ REPL[25]:1
Some type information was truncated. Use `show(err)` to see complete types.
> ] st
  [fbb218c0] BSON v0.3.9
  [49dc2e85] Calculus v0.5.1
  [e084ae63] CoolProp v0.2.0
  [ffbed154] DocStringExtensions v0.9.3
  [06fc5a27] DynamicQuantities v0.13.1
  [615f187c] IfElse v0.1.1
  [033835bb] JLD2 v0.4.46
  [961ee093] ModelingToolkit v9.7.0
  [16a59e39] ModelingToolkitStandardLibrary v2.6.0
  [2774e3e8] NLsolve v4.5.1
  [77ba4419] NaNMath v1.0.2
  [1dea7af3] OrdinaryDiffEq v6.74.1
  [731186ca] RecursiveArrayTools v3.12.0
  [189a3867] Reexport v1.2.2
  [3059ac8e] RefpropSplines v0.5.4
⌃ [0bca4576] SciMLBase v2.30.1
  [efcf1570] Setfield v1.1.1
  [90137ffa] StaticArrays v1.9.3
  [2efcf032] SymbolicIndexingInterface v0.3.11
  [d1185830] SymbolicUtils v1.5.1
  [0c5d862f] Symbolics v5.25.2
  [3a884ed6] UnPack v1.0.2
  [37e2e46d] LinearAlgebra
  [44cfe95a] Pkg v1.10.0
  [2f01184e] SparseArrays v1.10.0
  [4607b0f0] SuiteSparse
  [bea87d4a] SuiteSparse_jll v7.2.1+1
ChrisRackauckas commented 6 months ago

works with latest versions

ChrisRackauckas commented 6 months ago
julia> prob = ODEProblem(sysRed, [], (0.0, 1.0))
ODEProblem with uType Nothing and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: nothing

julia> sol = solve(prob, Rodas5())
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 1.0
u: 2-element Vector{Vector{Float64}}:
 []
 []

julia> 

julia> sol[x]
2-element Vector{Float64}:
 1.0
 1.0