ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

List of missing interface function (where numeric indexes must be used instead of `something[x]`)

Open TorkelE opened this issue 1 year ago • 15 comments

Compiled a list of various of these that I have encountered that should probably fixed (so that one is not required to fetch the correct index of the system and insert that instead).

This following list should be fairly self-explanatory, I'm not sure whether we wish to support all of these after the latest MTK update, but there are definitely missing cases.

using ModelingToolkit, NonlinearSolve

# Cases for `NonlinearSystem`s
begin 
    # Generates system with states: x, parameters: p, observables: y.
    @parameters p
    @variables t x(t) y(t)
    eqs = [0 ~ p - x^3 + 2x^2,
           0 ~ p - y]    
    @named nsys = NonlinearSystem(eqs, [x, y], [p])
    nsys = structural_simplify(nsys)

    # Creates NonlinearProblem, Integrator, and Solution

    nl_prob = NonlinearProblem(nsys, [x => 1.0, y => 3.0], [p => 2.0])
    nl_integrator = init(nl_prob)
    nl_sol = solve(nl_prob)
end

nl_prob[x] # ERROR: type NonlinearSystem has no field iv
nl_prob[y] # ERROR: type NonlinearSystem has no field iv
nl_prob[p] # ERROR: type NonlinearSystem has no field iv
nl_prob[nsys.x] # ERROR: type NonlinearSystem has no field iv
nl_prob[nsys.y] # ERROR: type NonlinearSystem has no field iv
nl_prob[nsys.p] # ERROR: type NonlinearSystem has no field iv
nl_prob[:x] # ERROR: type NonlinearSystem has no field iv
nl_prob[:y] # ERROR: type NonlinearSystem has no field iv
nl_prob[:p] # ERROR: type NonlinearSystem has no field iv

nl_prob[x] = 5.0 # Works.
nl_prob[p] = 5.0 # Works.
nl_prob[nsys.x] = 5.0 # Works.
nl_prob[nsys.p] = 5.0 # Works.
nl_prob[:x] = 5.0 # Works.
nl_prob[:p] = 5.0 # Works.

nl_integrator[x] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[y] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[p] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[nsys.x] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[nsys.y] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[nsys.p] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[:x] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[:y] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...
nl_integrator[:p] # ERROR: MethodError: no method matching getindex(::NonlinearSolve.FastShortcutNonlinearPolyalgCache{...

nl_sol[x] # Works.
nl_sol[y] # Works.
nl_sol[p] # Works.
nl_sol[nsys.x] # Works.
nl_sol[nsys.y] # Works.
nl_sol[nsys.p] # Works.
nl_sol[:x] # ERROR: UndefVarError: `x` not defined
nl_sol[:y] # ERROR: UndefVarError: `y` not defined
nl_sol[:p] # Works.

# Cases for `ODESystem`s
using OrdinaryDiffEq
begin 
    # Generates system with states: x, parameters: p, observables: y.
    @parameters p
    @variables t x(t) y(t)
    D = Differential(t)
    eqs = [D(x) ~ p - x^3 + 2x^2,
           D(y) ~ p - y]    
    @named osys = ODESystem(eqs)
    osys = structural_simplify(osys)

    # Creates NonlinearProblem, Integrator, and Solution

    ode_prob = ODEProblem(osys, [x => 1.0, y => 3.0], (0.0, 10.0), [p => 2.0])
    ode_integrator = init(ode_prob, Tsit5())
    ode_sol = solve(ode_prob, Tsit5())
end

solve(ode_prob, Tsit5(); save_idxs=1) # Works
solve(ode_prob, Tsit5(); save_idxs=x) # ERROR: ArgumentError: invalid index: x(t) of type SymbolicUtils.BasicSymbolic{Real}
solve(ode_prob, Tsit5(); save_idxs=y) # ERROR: ArgumentError: invalid index: y(t) of type SymbolicUtils.BasicSymbolic{Real}
solve(ode_prob, Tsit5(); save_idxs=osys.x) # ERROR: ArgumentError: invalid index: x(t) of type SymbolicUtils.BasicSymbolic{Real}
solve(ode_prob, Tsit5(); save_idxs=osys.y) # ERROR: ArgumentError: invalid index: y(t) of type SymbolicUtils.BasicSymbolic{Real}
solve(ode_prob, Tsit5(); save_idxs=:x) # ERROR: ArgumentError: invalid index: :x of type Symbol
solve(ode_prob, Tsit5(); save_idxs=:y) # ERROR: ArgumentError: invalid index: :y of type Symbol

TorkelE avatar Nov 01 '23 02:11 TorkelE

It would be good to add an array variable to this list, since they are often broken in my experience, i.e.,

@variables z(t)[1:2]

and then trying to access the solution or problem at z.

baggepinnen avatar Nov 01 '23 05:11 baggepinnen

Previously, basically all these cases has to be hard-coded, almost case by case. Is any united internal interface for setfield! and getfield planned/in the works? That would make it more avoidable to miss edge cases and stuff.

TorkelE avatar Nov 01 '23 13:11 TorkelE

The unified way is the SymbolicIndexingInterface.jl overhaul + removing the Symbol dispatches and requiring systems are complete (MTK v9)

ChrisRackauckas avatar Nov 01 '23 13:11 ChrisRackauckas

Could we keep teh symbol interfaces, e.g. sol[:x]. Removing that is going to permanently break Catalyst.

TorkelE avatar Nov 01 '23 13:11 TorkelE

sol[sys.x] is better for Catalyst since it makes less assumptions which break the general interfaces. It only requires that Catalyst models are made complete.

ChrisRackauckas avatar Nov 01 '23 16:11 ChrisRackauckas

We shouldn't be using Symbols in Catalyst except in the u0 and parameter mappings (which was done as syntactic sugar to avoid having to unpack every variable/parameter from a DSL generated system). There shouldn't be anywhere else we are using symbols or telling users to use them (and if there is the associated tutorial should get updated).

isaacsas avatar Nov 01 '23 17:11 isaacsas

Has there been any consideration to making a macro to construct parameter/u0 mappings given a complete system? i.e. in something like

u0 = [
    rc_model.capacitor.v => 0.0,
]

the macro could be

u0 = @symbolic_mapping rc_model begin
    capacitor.v => 0.0
end

the idea being that if a user has a lot of initial conditions or parameters to specify this would allow them to avoid having to type rc_model.variablename before each variable/parameter?

isaacsas avatar Nov 01 '23 17:11 isaacsas

Maybe we can discuss in person next time we meet to see if we can expand to allow :x usage? I can see some obscure cases where it might not be preferred, but for the vast majority of users and cases it is.

TorkelE avatar Nov 01 '23 17:11 TorkelE

I think a macro that makes the syntax a little nicer is fine. That's going to be more correct than symbols

ChrisRackauckas avatar Nov 01 '23 18:11 ChrisRackauckas

Ideally I would like Symbols to work for non-hierarchial systems (and throw an error if they are not, I think that is how it works currently)

TorkelE avatar Nov 01 '23 18:11 TorkelE

Most of this should work, bar the remake stuff (which is on the list)

AayushSabharwal avatar Feb 23 '24 06:02 AayushSabharwal

Can you show it? Is this all captured in a test?

ChrisRackauckas avatar Feb 23 '24 06:02 ChrisRackauckas

I can add it to a test. The cases that don't work:

  • Indexing nl_integrator. This isn't an integrator at all, and is only defined in NonlinearSolve. Is this something that is commonly indexed? Can it subtype something appropriate in SciMLBase, or would that break things and require SII to be implemented in NonlinearSolve?
  • Indexing using Symbols as parameters (nl_sol[:p]). As of MTKv9, indices are stored based on the hash of the symbolic variable, so this just doesn't work. I could add hashes for Symbol(var) but that would immediately double the size of the Dict.
  • remake. As mentioned above, this will eventually work
  • save_idxs. This is a bit of a tough one to handle, since it'll require changes in multiple libraries (ODE and SciMLBase) and require the solution to manually implement variable_index if save_idxs is used, since the system's index is now invalid.

AayushSabharwal avatar Feb 23 '24 06:02 AayushSabharwal

This isn't an integrator at all, and is only defined in NonlinearSolve. Is this something that is commonly indexed? Can it subtype something appropriate in SciMLBase, or would that break things and require SII to be implemented in NonlinearSolve?

It's an analogous iterator type to the integrator and should get the same interface.

Indexing using Symbols as parameters (nl_sol[:p]). As of MTKv9, indices are stored based on the hash of the symbolic variable, so this just doesn't work. I could add hashes for Symbol(var) but that would immediately double the size of the Dict.

Yup that's understood.

remake. As mentioned above, this will eventually work

Yes separate issue.

save_idxs. This is a bit of a tough one to handle, since it'll require changes in multiple libraries (ODE and SciMLBase) and require the solution to manually implement variable_index if save_idxs is used, since the system's index is now invalid.

Separate issue.

So let's split this into the 3 concrete issues with MWEs?

ChrisRackauckas avatar Feb 23 '24 06:02 ChrisRackauckas

So let's split this into the 3 concrete issues with MWEs?

Sounds good

AayushSabharwal avatar Feb 23 '24 06:02 AayushSabharwal