breaking change in behavior of complete in non-breaking release

[isaacsas]

Calling complete(sys) now changes the underlying parameters of a system, interjecting lots of Initial(unknown)s into the parameter list. This seems breaking, as anything that previously worked with completed systems, like various Catalyst functionality, suddenly has lots of non-parameter parameters to deal with when calling get_ps(sys).

Older MTK v9s (like 9.59):

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

@parameters k1 
@variables X1(t) 
eqs = [
    D(X1) ~ k1 * X1
]
@mtkbuild osys = ODESystem(eqs, t, [X1], [k1])

such that

julia> ModelingToolkit.get_ps(osys)
1-element Vector{Any}:
 k1

Newer MTK v9s:

julia> ModelingToolkit.get_ps(osys)
2-element Vector{Any}:
 k1
 Initial(X1(t))

[ChrisRackauckas]

Those are all parameters? There has never been anything guaranteeing there won't be other parameters, and in fact there are many cases where this is already happened?

[ChrisRackauckas]

Also get_ps is not for the parameters, it's for the internal ps array which is the internal representation of the parameters. But everything there is still a parameter so it's not like it broke any contract.

[isaacsas]

Anywhere that previously assumed that in a non-hierarchical or flattened system issetequal(get_ps(sys), parameters(sys)) == true, and used them interchangeably, is now broken. This behavior had been true since the first versions of MTK. Likewise, the behavior of complete has now also changed, for I believe the second(?) time since v9 was released.

It is impossible to build stuff on MTK without using get_ps, parameters, or complete. If the way that such core API functions work keeps changing in non-breaking releases, peoples' codes built on MTK will be continually breaking...

[ChrisRackauckas]

Where did we say those were equal? They were made as different calls because they don't do the same thing... why would two different functions that do two different things always be the same? There are many systems for which that already wasn't true anyways.

[isaacsas]

• parameters(sys): All parameters of the system and its subsystems.
• get_ps(sys): Parameters that are in the current-level system.

from

https://docs.sciml.ai/ModelingToolkit/stable/basics/AbstractSystem/#Composition-and-Accessor-Functions

That is the long-standing documentation for them in the AbstractSystem interface. That has no description of them being different in any way beyond that the former also includes parameters from subsystems. Nowhere does it state that get_ps(sys) could be different from parameters(sys) for a system that has no subsystems.

[TorkelE]

Whatever the new intended behaviour is, in that case it seems that we need a new, public version, of get_ps? I.e. in Catalyst we have a nice section discussing how to access stuff in models: https://docs.sciml.ai/Catalyst/stable/model_creation/reactionsystem_content_accessing/#model_accessing_hierarchical_symbolic_variables-2. Here we describe how while parameters can be used generally, for hierarchical models there is get_ps and parameters, both which have different uses. If get_ps is now some internal thing I need to replace get_ps in Catalyst documents with something else, and then create a new description of get_ps (and similar ones) in the Catalyst devdocs.

[ChrisRackauckas]

Let's be clear here. get_xxxxx is just getfield(sys, :xxxxx) where sys.x is overloaded for hierarchical modeling. It's a simplified getfield call. It's intended to give whatever the internal representation of the type is doing, and it's been that way since hierarchical modeling was added 4 years ago. "is now some internal thing", all of the get_xxxxx have always been just a getfield call and are all just generated functions created by looping over the fields to create one accessor function per field. This hasn't ever changed, and it still hasn't changed. We should keep it like that, since if we have one get_xxxx function that is different from the others that would be fairly odd and can lead to some unintended issues with the function generation. If we need another function we should make another functions, but keep the get_xxxx as just getfield aliases.

[ChrisRackauckas]

For reference, he's the code https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/abstractsystem.jl#L799-L874

and it has been just a generated alias for getfield since af37b7f#diff-287e017b6c265948a34862d11b63bade810cb38f11173644350028d4ba836661R132

[isaacsas]

Yes, but the meaning of the field has now changed.

Moreover, previously in a flattened/non-hierarchical system get_ps(sys) was an allocation-free way to get the list of parameters. parameters(sys) calls unique so will always allocate a new array. Because of that we preferred using get_ps in many places where we knew that a system was flattened.

[ChrisRackauckas]

Well it still gives you that, since the initial values are parameters?

[isaacsas]

I'll close this as it doesn't feel like there is anything productive to continue discussing about this and it sounds like this was an intentional change that won't be modified. At the end of the day though, this change has broken stuff in Catalyst, as have many of the other changes in API function behavior over the last 5+ months. This has resulted in our being unable to have a fully working version of Catalyst master for most of that time, and prevented us from making a release going back September.

[ChrisRackauckas]

If the way that such core API functions work keeps changing in non-breaking releases, peoples' codes built on MTK will be continually breaking...

Just to be clear here, two things. One, if you rely on the exact values in a field because you're using a getfield call directly, then you're prone to breakage because then you're relying on implementation not interface. This is why we tend to only document subsets of fields across the board for all of SciML, any change is breaking if someone is relying not on interface but implementation details. I can see why you'd skip the main interfaces to eek out a bit of performance, but that's always going to have a trade-off unless we make a function for the thing you're looking for.

But also, the FUD is very misplaced here. This is the change to uniformity, and it's the change you requested 😅. The steps to get here were:

1. Initialization systems added for non-trivial DAE systems. The issue though is that if it's only for DAEs, we cannot handle things like y => 2 for y ~ x an observed variable in an ODE system
2. Initialization systems added for non-trivial DAEs and a subset of ODEs for which initial conditions are placed on non-ODE values. But the issue here is that we also wanted to support parameter equations p2 ~ 2*p1. If we want to support those, they need to resolve at the same time as the initialization equations.
3. Then there's lots of edge cases because it's hard for a user to understand when there will be an initialization problem or not. So the latest is a major simplification: everything gets an initialization problem. Many times that initialization problem is just trivial, so when it's trivial we do some things eagerly.

The hope though is that now that everything has an initialization problem, which is well-motivated since it ends up being required for correctness across the board, that means there's no more problems that should ever switch from "not having an initialization problem" to "having an initialization problem": there is no choice anymore, they are always on the same side.

One thing required to do that correctly though is that we have to really understand the initial values as parameters. In that sense, (1) and (2) become the same problem, simplifying the mechanism a lot.

I don't think we're at the final state yet, mostly because of some interactions like https://discourse.julialang.org/t/modelingtoolkit-changes-with-remake/126004 where if you directly modify values non-symbolically it can be a bit unintuitive, but we're hopefully pretty close now.

The reason we didn't roll out initialization systems to everything at the start was that we couldn't run an initialization system on non-DAEs in the first place, so we had to update the whole ecosystem to allow for intiailization problems everywhere (not just ODEs, but optimization, nonlinear solves, etc.), and now with that it can be standard.