-
-
Notifications
You must be signed in to change notification settings - Fork 220
sort_eqs = true
leads to surprising results of structural_simplify
#3633
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Comments
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. 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 The odd case out of course is specifically if you have only ever written ODEs before. Specifically for (some subset) of ODEs 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 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 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. |
Diving into this a bit, one other case that will come up is integrals. If
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. |
IMO the default behaviour of
structural_simplify
caused bysort_eqs = true
is very surprising, and actually leads to problems in our downstream application.Take the following example:
In contrast to what was desired (and presumably expected) the order of state variables changed from
to
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 ofstructural_simplify
to the user.The text was updated successfully, but these errors were encountered: