`sort_eqs = true` leads to surprising results of `structural_simplify`

[devmotion]

IMO the default behaviour of structural_simplify caused by sort_eqs = true is very surprising, and actually leads to problems in our downstream application.

Take the following example:

t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits
ps = @parameters Ka CL Vc
dvs = @variables Depot(t) Central(t)
eqs = [
  D(Depot) ~ -Ka * Depot
  D(Central) ~ Ka * Depot - CL / Vc * Central
]
@mtkbuild sys = ODESystem(eqs, t, dvs, ps)

In contrast to what was desired (and presumably expected) the order of state variables changed from

julia> dvs
2-element Vector{Num}:
   Depot(t)
 Central(t)

to

julia> ModelingToolkit.unknowns(sys)
2-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 Central(t)
 Depot(t)

Arguably changing the order of states/unknowns is much more surprising than changing the order of parameters, in particular in a case where @mtkbuild does not perform any actual simplification.

In our application, users may refer to states/compartments by integers when defining events/doses, so the reordering is not only surprising but actually breaks existing code (arguably, referring to compartments by integers is not ideal but that's how it's done in common data formats in the field, and telling people they can't use their existing datasets anymore doesn't seem feasible).

I think one could avoid these surprises in many (all?) cases by making sort_eqs an internal implementation detail but not user-facing. I'm not sure if it's possible since I'm not familiar with the implementation, but in principle I would assume that the tearing algorithm could operate on the sorted equations but the states of the resulting simplified system could be permuted such that they are in the same order as the original system, before returning the result of structural_simplify to the user.

[ChrisRackauckas]

It has been repeated and documented for years: do not rely on the exact ordering of variables because you cannot know what they will be. Do not use integers to refer to outputs. Etc. This has been mantra of using the library since about day 1. structural_simplify can and will reorder equations, do not rely on the order. So no, it's not breaking if structural_simplify changes the order of equations: it always has, users have always been told that, and it always will. We do not integer index in the documentation, everything points to using the SymbolicIndexingInterface (SII), and nothing relies on equation order. It's a core tenant of the library.

I don't know if we need to put this into bigger and bolder big red letter text in the documentation and docstrings or something but this has been one of the core tenants of the library since the day structural simplification and DAE support was added back in 2021.

And what you get out of this is arguably a more predicable sort order than what existed before: you previously could get different equation orders out of structural_simplify based on the connection ordering for the same set of equations, now with the sorting you always get the same ordering for any equivalent system. So this actually makes the final equation ordering less surprising in the majority of cases.

The odd case out of course is specifically if you have only ever written ODEs before. Specifically for (some subset) of ODEs structural_simplify didn't do anything on them and the equation order was passed through. But that isn't something to be relied on because, of course, it's more of a current limitation of the current compiler. A few things that will change this in the near future are implicit ODE solver optimizations (which should be released in a month or so?) and integrating analytical solvers into structural simplification (which hopefully some of this comes from this summer's GSoC?). So we take the interface very seriously which has always been that structural_simplify can and will change not only the ordering of the equations but it will also change the number of equations. This is demonstrated in the very first tutorial https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/ and the user is greeted with the right ways to interact with the solution via SII from the get go https://docs.sciml.ai/ModelingToolkit/stable/tutorials/ode_modeling/#Named-Indexing-of-Solutions.

I think the question here is then, is it worth the effort for something mostly focused on simple ODEs to update to SII? I think for something like Pumas, in this round, it might as well set sort_eqs=false and continue on is okay-ish. While it's relying on UB (and has for like 4 years now), at least for now there is nothing else that touches the ordering in simple ODE cases so it'll get the next version bump working immediately. That said, https://docs.sciml.ai/SymbolicIndexingInterface/stable/usage/ sol[SII.allvariables] will do this pretty simply, and we should really just update libs at this point that are still relying on UB.

But, Pumas in particular should be anticipating the two parts I mentioned earlier, the implicit solver optimizations and the analytical solution pieces, which will not only change ordering but also the number of equations actually solved for. For example:

t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits
ps = @parameters Ka CL Vc
dvs = @variables Depot(t) Central(t)
eqs = [
  D(Depot) ~ -Ka * Depot
  D(Central) ~ Ka * Depot - CL / Vc * Central
]
@mtkbuild sys = ODESystem(eqs, t, dvs, ps)

this should in the fairly near future give you an ODEProblem where isempty(prob.u0) and sol[Depot] just gives you the analytical solution. This is of course the kind of transformation that is already done by structural_simplify in the vast majority of cases (DAEs do this equation changing for almost every set of equations, there's usually something to eliminate), but again a specific subset of cases has traditionally avoided it because our compiler wasn't smart enough to do anything there. It's now getting smart enough, we plan to use it.

Pumas can and should just turn off some optimizations for now in order to get a bump done, it specifically won't see any improvement from the sorting for now so why not, but we should then discuss separately how to best use the tooling here, in particular the full observed function in order to fully get the benefits in the near future. Note it has already been setup for this path for awhile now. Pumas does not expose the ODESolution object to the user, this is intentional (because of some previous attempts at analytical solving built into Pumas, which then noted that it needed bigger and better guns to work properly, so we're finally getting around to completing that story arch), so it just needs to SII to the matrix and then the user is fine. I'll just make the PR so that it's using the "modern" (2021) interface that is safe for any future change.

[ChrisRackauckas]

Diving into this a bit, one other case that will come up is integrals. If D(u) ~ 1 we currently do not replace it with an observed u ~ t but we should. More generally, D(u) ~ f(t) expressions should use the symbolic integration tooling (being developed this summer in a GSoC) and if a symbolic expression exists, transform u ~ g(t). Otherwise it can at least remove the equation from the set of ODEs into a callback https://docs.sciml.ai/DiffEqCallbacks/stable/integrating/. So 3 cases:

1. Analytical solutions to subsets of equations
2. Equation sorting
3. Removal of integral expressions

can and will effect non-DAE ODEs in the near future. So... it's coming more and more. But we should definitely help ramp the purely ODE users in however we can because while this stuff has always been true for DAEs, purely ODE use cases have traditionally gotten away with some UB just because we didn't have the machinery to do anything fancy.