Allow getu to be used directly on the state vector?

[Datseris]

Is your feature request related to a problem? Please describe.

I am creating a dynamical system that iterates a vector of vectors, that is, many different states in parallel. I do this because it is important that the states are iterated at exactly the same time. The system is in essence using an MTK model and gets the function f and uses it in a for loop to create a new newf that is the ODE rule for the system with the parallel states.

I would like to access "observed" variables within the MTK definition and symbolic indexing.

Describe the solution you’d like

getu would allow this, if it could be applied directly on the state instead of an ODEproblem/sol/integ. Because the ODEInteg of the parallel system is not created from an MTK system, I can't pass it to getu.

Describe alternatives you’ve considered

I am welcoming alternative ways to achieve that.

[Datseris]

Sorry, to be clear: above I mean to allow the function generated by getu to be applied directly on the state. I am keeping track of the original MTKSystem and hence I can create the function by getu, I just can't apply it.

[Datseris]

wait a moment, maybe https://docs.sciml.ai/SymbolicIndexingInterface/stable/api/#SymbolicIndexingInterface.observed is what I need. I will test this now!

[Datseris]

Hm, observed doesn't seem to work as documetned:

MethodError: no method matching symbolic_container(::ODESystem)

  Closest candidates are:
    symbolic_container(::SciMLBase.AbstractSciMLFunction)       
     @ SciMLBase C:\Users\gd419\.julia\packages\SciMLBase\dpafx\src\scimlfunctions.jl:4220
    symbolic_container(::SciMLBase.AbstractSciMLProblem)        
     @ SciMLBase C:\Users\gd419\.julia\packages\SciMLBase\dpafx\src\problems\problem_interface.jl:1
    symbolic_container(::SciMLBase.AbstractOptimizationSolution)     @ SciMLBase C:\Users\gd419\.julia\packages\SciMLBase\dpafx\src\solutions\optimization_solutions.jl:102
    ...

  Stacktrace:
   [1] observed(sys::ODESystem, sym::Int64)
     @ SymbolicIndexingInterface C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\65Q0o\src\interface.jl:100
   [2] observe_state(u::Vector{Float64}, p::Vector{Float64}, prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#k#503"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x855e1841, 0xb46b4383, 0x84212e78, 0x4808e69a, 0xd156bdd9), Nothing}, RuntimeGen

where sys is the system generated by the first tutorial of MTK:

sys = structural_simplify(connected)
SymbolicIndexingInterface.observed(connected_simp, fol_1.τ)

[AayushSabharwal]

observed isn't valid for systems, since the observed function isn't actually generated yet. It's valid for SciMLFunctions, integrators, problems and solutions.

[AayushSabharwal]

That said, I think it's a nice idea to allow the function returned by getu to be used on state vectors directly. Thoughts @ChrisRackauckas ?

It would look something like

getter = getu(sys, sym)
# assume sym is the second variable
getter([1.5, 2.5, 3.5]) # returns 2.5

And the best part is I think all we will need is SII.state_values(arr::AbstractArray) = arr for it to work:

julia> using SymbolicIndexingInterface

julia> SymbolicIndexingInterface.state_values(arr::AbstractArray) = arr

julia> sc = SymbolCache([:x, :y])
SymbolCache{Vector{Symbol}, Nothing, Nothing}([:x, :y], nothing, nothing)

julia> getter = getu(sc, :x)
(::SymbolicIndexingInterface.var"#getter#8") (generic function with 1 method)

julia> getter([1.5, 2.5])
1.5

julia> getter = getu(sc, :y)
(::SymbolicIndexingInterface.var"#getter#8") (generic function with 1 method)

julia> getter([1.5, 2.5])
2.5

[Datseris]

thanks, yeah that would be very nice to have!

[AayushSabharwal]

This will be supported once #33 is merged and tagged

[Datseris]

@AayushSabharwal I have checked out your PR in my local code, where I am working in integrating DynamicalSystems.jl with ModelingToolkit.jl. It works wonderfully! It does exactly what I wanted, and it significantly simplifies my code as well! Coincidentally, it turns out I wanted the same thing for parameters as well, for getp to be applicable on the resulting parameter vector, so thanks for including that!

[Datseris]

@AayushSabharwal I am having some difficulty using the PR that I cannot create a MWE for (because I have not understood why the problem exists at all). It appears that for some ODESystem this PR works, while for some others it does not. I am always passing an ODESystem in the getu(sys, index) function. Sometimes it works as expected, sometimes I get the error

ERROR: MethodError: no method matching symbolic_container(::ODESystem)


Stacktrace:
 [1] observed(sys::ODESystem, sym::Num)
   @ SymbolicIndexingInterface C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\rvbno\src\interface.jl:100
 [2] _getu(sys::ODESystem, ::SymbolicIndexingInterface.ScalarSymbolic, ::SymbolicIndexingInterface.ScalarSymbolic, sym::Num)
   @ SymbolicIndexingInterface C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\rvbno\src\state_indexing.jl:149
 [3] getu
   @ C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\rvbno\src\state_indexing.jl:113 [inlined]

Above you said that observed does not allow for ODESystems. So in the end, what is the correct way to use getu? Do you pass in an ODESystem or an ODEProblem generated from the system?

[Datseris]

A!!!! I Was able to reproduce the problem!!! The syntax doesn't work for observed variables!

@parameters σ=28.0 ρ=10.0 β=8/3
@variables t x(t)=5.0 y(t)=0.0 z(t)=0.0 w(t)
D = Differential(t)

eqs = [D(x) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ w - β * z,
    w ~ x*y]

@named lorenz = ODESystem(eqs)
sys = structural_simplify(lorenz)

julia> SymbolicIndexingInterface.getu(sys, sys.x)
(::SymbolicIndexingInterface.var"#getter#24"{SymbolicIndexingInterface.var"#_getter#23"{Int64}}) (generic function with 2 methods)

julia> SymbolicIndexingInterface.getu(sys, sys.w)
ERROR: MethodError: no method matching symbolic_container(::ODESystem)

Stacktrace:
 [1] observed(sys::ODESystem, sym::Num)
   @ SymbolicIndexingInterface C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\rvbno\src\interface.jl:100
 [2] _getu(sys::ODESystem, ::SymbolicIndexingInterface.ScalarSymbolic, ::SymbolicIndexingInterface.ScalarSymbolic, sym::Num)
   @ SymbolicIndexingInterface C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\rvbno\src\state_indexing.jl:149
 [3] getu(sys::ODESystem, sym::Num)
   @ SymbolicIndexingInterface C:\Users\gd419\.julia\packages\SymbolicIndexingInterface\rvbno\src\state_indexing.jl:113
 [4] top-level scope

The above works fine if I give the ODE problem as the first input, and not the system.

so in the end, getu cannot be used with the system, only the problem instance. So perhaps the docstring should rename sys to prob/sol/integ instead.

[AayushSabharwal]

It works for sys, just not for observed variables. It's a bit of a weird thing to document, since getu says it is valid if the first argument supports observed, which is an implementation detail. Since this is an interface designed to be independent of implementation, it can't assume that all systems do/do not implement observed. In general for MTK, if support for observed quantities is required then problems must be used.