Downstream Updates from MTK v9

ChrisRackauckas commented 8 months ago

Upgrade bug issues:

(None so far!)

Upgrade explanation issues:

ChrisRackauckas commented 8 months ago

Also any downstream community packages I may have missed, please post here. I tried to weed out unmaintained projects and things that just look like archives of project code.

@ctessum @cadojo @jake484 @hexaeder @jedforrest @TorkelE @isaacsas @TS-CUBED @wsphillips @baggepinnen @Vaibhavdixit02 @AlCap23 @sumiya11 @michakraus @RafaelArutjunjan @xiaomingfu2013 @bgctw @augustinas1 @voduchuy @JordiBolibar @facusapienza21 @sebapersson @iliailmer @orebas @mattar13 @killah-t-cell @nathanaelbosch @david-pl @paulflang @pogudingleb @YichengDWu @RXGottlieb @meyde @paulobruno @piluc @xtalax please see https://github.com/SciML/ModelingToolkit.jl/blob/master/NEWS.md and let me know if you have any questions. Some general updating notes:

Any PDE package should be able to update without a hiccup. PDESystem was unchanged.
Many packages which just use MTK for ODEs will update straightforwardly. We would recommend using the canonical t and D, i.e. using ModelingToolkit: t_nounits as t, D_nounits as D instead of defining it in file, but your tests shouldn't be broken from this. Most test breakage is likely just because ODESystem(eqs) was deprecated for ODESystem(eqs,t) and that the system needs to be completed before generating the problem (i.e. use @mtkbuild or change structural_simplify(sys) to complete(structural_simplify(sys)), though @mtkbuild is now the documented way). It's not possible to automatically upgrade these so it had to be made breaking (t detection did not always work and there were subtle edge cases) and so this will need a manual intervention. I presume this will be 99% of the upgrade work for most packages but should be straightforward and mechanical.
states(sys) was changed to unknowns(sys), that should be mechanical.
Symbolic arrays have had many changes so that they can now keep structure. If you were previously using u[1:5] and you meant an array of symbolic variables, use u = Symbolics.variables(:u, 1:5) instead. To clarify the difference between them, see https://symbolics.juliasymbolics.org/stable/manual/arrays/. I can see that the former was used in a lot of places that meant the latter. If you're not being setup for JuliaSimCompiler there won't be a performance difference (and in fact there will likely be a performance improvement due to improved tracing performance)
If your package assumes parameters from the generated ODEProblem are a vector (which IIUC is a very small subset, but includes PETab.jl), use the SciMLStructures interface https://github.com/SciML/SciMLStructures.jl to interact with the parameters structure. For these cases you may need to come talk to me since there are not that many examples of handling this right now. Also, some of these may need https://github.com/SciML/SciMLSensitivity.jl/pull/1010

Again if any of this is confusing and you need help upgrading, please feel free to ask questions, ask for help. We want help to get the whole ecosystem updated.

(@AlCap23 you were going to move DynamicOED.jl into SciML? This would be a good time to do it)

AayushSabharwal commented 8 months ago

change structural_simplify(sys) to complete(structural_simplify(sys)),

structural_simplify calls complete, so this is not necessary

JumpProcesses? Catalyst?

yes and yes

SturcturalIdentifiability?

yes

ChrisRackauckas commented 8 months ago

I was just making a note that I haven't gotten those CompatHelper PRs yet. We're definitely covering any non-archived SciML package!

ChrisRackauckas commented 8 months ago

A more complete description of the parameter splitting if anyone needs help figuring out what it's all about: https://github.com/SciML/ModelingToolkit.jl/issues/2482#issuecomment-1959330130

sumiya11 commented 8 months ago

Thanks for this summary @ChrisRackauckas ! We will try to update MTK in StructuralIdentifiability.jl

cadojo commented 8 months ago

Thank you for the detailed information! I am trying to update AstrodynamicalModels.jl, and I'm running into an issue with ODEFunction. I think I have traced the issue to different parameter ordering inside functions written by Symbolics.generate_function for ODESystem instances created by MTK Version 9. **Is the usage in PlanarEntrySystem() still supported? I had thought that ODESystem(eqs, t, [γ, v, r, θ], [R, P, H, m, A, C, μ]; name=name) was supported as part of the non-DSL syntax.

I can see the word-level differences in the generated function string in git diff --no-index, but those differences are shown by colors and I don't know how to show that in GitHub markdown or HTML. I've included an image of the terminal output below. Alternatively, you can run the following command.

git clone git@github.com:cadojo/AstrodynamicalModels.jl.git
cd AstrodynamicalModels.jl
git diff --no-index --word-diff=color v8/entry.jl v9/entry.jl

pogudingleb commented 8 months ago

@ChrisRackauckas What happened to ModelingToolkit.DiscreteSystem ? Now the subtypes of ModelingToolkit.AbstractTimeDependentSystem seem to be only JumpSystem and ModelingToolkit.AbstractODESystem.

pogudingleb commented 8 months ago

@ChrisRackauckas What happened to ModelingToolkit.DiscreteSystem ? Now the subtypes of ModelingToolkit.AbstractTimeDependentSystem seem to be only JumpSystem and ModelingToolkit.AbstractODESystem.

My bad, found in the release notes

ChrisRackauckas commented 8 months ago

It will be returning, we just need to redo it a bit on the new discrete-continuous infrastructure

ChrisRackauckas commented 8 months ago

Thank you for the detailed information! I am trying to update AstrodynamicalModels.jl, and I'm running into an issue with ODEFunction. I think I have traced the issue to different parameter ordering inside functions written by Symbolics.generate_function for ODESystem instances created by MTK Version 9. **Is the usage in PlanarEntrySystem() still supported? I had thought that ODESystem(eqs, t, [γ, v, r, θ], [R, P, H, m, A, C, μ]; name=name) was supported as part of the non-DSL syntax.

Ordering isn't guaranteed. In fact, parameters aren't even a vector so that doesn't make a lot of sense. See the discussion on parameter splitting for more details: https://github.com/SciML/ModelingToolkit.jl/issues/2482#issuecomment-1959330130.

But anything you were looking to do that needed a parameter ordering can and should be done by using the symbolic indexing interface: https://docs.sciml.ai/SymbolicIndexingInterface/stable/. For example, f = getp(prob, [R, P, H, m, A, C, μ]) gives you a function f(prob) that returns the parameters in that chosen ordering. This decouples the ordering and representation in the solver from how you use and interact with it. There's also setp as well.

pogudingleb commented 8 months ago

It will be returning, we just need to redo it a bit on the new discrete-continuous infrastructure

Does it mean that right now MTK does not have an interface for discrete-time systems?

cadojo commented 8 months ago

But anything you were looking to do that needed a parameter ordering can and should be done by using the symbolic indexing interface: https://docs.sciml.ai/SymbolicIndexingInterface/stable/. For example, f = getp(prob, [R, P, H, m, A, C, μ]) gives you a function f(prob) that returns the parameters in that chosen ordering. This decouples the ordering and representation in the solver from how you use and interact with it. There's also setp as well.

Are the ODEFunction instances intended for users to call, or only the solvers? To clarify, my motivation is this: I have been calling ODEFunction types inside other functions to get the Jacobian and vector field. Is there any way for users to call ODEFunction or the outputs of Symbolics.generate_function themselves? I can't seem to use Symbolics.generate_function as a way to build functions that users can call.

Edit

My goal is to use Symbolics.build_function or Symbolics.generate_function on an ODESystem to print the dynamics to a fast, non-allocating function for users (and solvers) to execute. This notebook is a full example for what I'm trying to do. I use my own subtypes of AbstractTimeDependentSystem (within GalacticPotentials.jl) to represent scalar potentials used in galactic dynamics, and then I use Symbolics code generation to print the dynamics to code.

AayushSabharwal commented 8 months ago

The best thing to support this for a general workflow would be for me to create a build_custom_function method in MTK, which is a wrapper over build_function that does the parameter ordering magic. It can be done manually, but it would rely on undocumented features and could break anytime.

To be able to use this build_custom_function (which I'll get to soon :D) you will need to ensure your systems can be completed, and require that this be done before any codegen. In your systems, you'll need:

complete::Bool
index_cache::Union{Nothing, MTK.IndexCache}

If you need to create a parameter object, use MTK.MTKParameters(sys, parammap). Don't rely on parameter order, or the internals of IndexCache or MTKParameters. If there's any other guidelines I come across while adding the function, I'll mention them here.

AayushSabharwal commented 8 months ago

Any functions inside ODEFunction should be callable as f(u, p, t) or f(du, u, p, t) as before. If they're not, please open an issue.

cadojo commented 8 months ago

Thank you @AayushSabharwal! Sounds good. To make sure I understand: the method signature for all ODEFunction types should still be (u, p, t) and/or (du, u, p, t), but the order of the elements in du, u, or p is unknown to the user at runtime. Is that right?

I don't think I understand how / why generate_function is a public method then. What is the use case to external packages if the generated function has unknown state (unknown) and parameter ordering? In my understanding, users can generate functions, but cannot call those generated functions. Is there some supported use case where users call ODEFunction with a variable-map-like structure? Or should ODEFunction and generate_function only for use in internal SciML codes? (I hope I'm not coming across as argumentative — I'm trying to understand how best to use SciML tools.)

ChrisRackauckas commented 8 months ago

Thank you @AayushSabharwal! Sounds good. To make sure I understand: the method signature for all ODEFunction types should still be (u, p, t) and/or (du, u, p, t), but the order of the elements in du, u, or p is unknown to the user at runtime. Is that right?

Again, the first thing to do is to stop assuming that there is such a thing as an "ordering" since p is not even a vector, so there is no such thing as p[1]. There are Cartesian indices and such, see the parameter splitting.

But also, see what I wrote above. If you use the standard setters and accessors https://docs.sciml.ai/SymbolicIndexingInterface/stable/ then you have everything that you need in order to use the p object.

I don't think I understand how / why generate_function is a public method then. What is the use case to external packages if the generated function has unknown state (unknown) and parameter ordering? In my understanding, users can generate functions, but cannot call those generated functions. Is there some supported use case where users call ODEFunction with a variable-map-like structure? Or should ODEFunction and generate_function only for use in internal SciML codes?

It probably should be considered internals now, that is a good point. We'd either need to make sure the MTKParameters construction is also documented API, or we need to assume that you built the problem so that prob.p is how you'd get it. I think I'd prefer just doing the latter as it decreases the API surface. I don't think you're losing anything either since building the problem to get prob.f is not much slower and it gives you all of the required pieces for the indexing. So @AayushSabharwal in the documentation update we may want to simply remove those and stop exporting the direct function generators. Since otherwise, if someone cannot build p, they cannot use those functions very well, and building p is pretty much the other thing that the XProblem constructors do, so we might as well say you have to build the problem.

This also fixes a few other buggy things because things which went around the problem construction have always made a few assumptions which don't quite hold true.

But again, I think you get all of the functionality you need back via getu, getp, setu, and setp, so I don't think anything is lost now that the SymbolicIndexingInterface is completed. It just means the better symbolic indexing tooling is more central to the story and we should highlight it better in the MTK docs.

The best thing to support this for a general workflow would be for me to create a build_custom_function method in MTK, which is a wrapper over build_function that does the parameter ordering magic. It can be done manually, but it would rely on undocumented features and could break anytime.

What would be the difference from just using build_function manually? It already lets you set the order through its definition.

ChrisRackauckas commented 8 months ago

My goal is to use Symbolics.build_function or Symbolics.generate_function on an ODESystem to print the dynamics to a fast, non-allocating function for users (and solvers) to execute. This notebook is a full example for what I'm trying to do. I use my own subtypes of AbstractTimeDependentSystem (within GalacticPotentials.jl) to represent scalar potentials used in galactic dynamics, and then I use Symbolics code generation to print the dynamics to code.

I think showing the generated code is a fine use case, especially to other targets. That one reason to keep it around as public API. We would need to expose the two parameter construction for the separate targets better in order to make this good though. I'm not sure the C code would work well for other parameter types at the moment, so we would need to make that error whenever we see any declared symtype.

ChrisRackauckas commented 8 months ago

Does it mean that right now MTK does not have an interface for discrete-time systems?

Yes, though I hope we can get around to getting something back next week. I think from other various discussions this has come up as a high priority. @AayushSabharwal let's discuss this and get something together.

cadojo commented 8 months ago

I think showing the generated code is a fine use case, especially to other targets. That one reason to keep it around as public API. We would need to expose the two parameter construction for the separate targets better in order to make this good though. I'm not sure the C code would work well for other parameter types at the moment, so we would need to make that error whenever we see any declared symtype.

If you do decide to keep the function generation methods in the public API, I'd be interested in contributing to the _build_function implementations for each target. I worked on those implementations a bit in Symbolics.jl right after it was released.

One way I could see keeping the function generation capabilities is to unpack all of the parameters as individual arguments. For CTarget, we could have a type signatures map to assign a C type to every parameter, based on their Julia types: Integer => "int", Bool => "uint8_t", Float32 => "float", Float64 => "double", AbstractArray => "double*". And we could use the AbstractSystem name as the function name and the parameter symbols as the argument names by default.

AayushSabharwal commented 8 months ago

What would be the difference from just using build_function manually? It already lets you set the order through its definition.

Yeah but we don't want to let people set the order 😅. The custom function would do something along the lines of build_function(expr, unknowns(sys), reorder_parameters(sys, parameters(sys)), get_iv(sys); kwargs...). Then probably it would also change the function definition to accept an MTKParameters and destructure it into subarrays internally, so it has a (u, p, t) signature.

AayushSabharwal commented 8 months ago

the order of the elements in du, u, or p is unknown to the user at runtime. Is that right?

As @ChrisRackauckas mentioned, there is no concept of an order for p since it isn't a vector anymore. u and du should still follow the order of unknowns(sys).

pogudingleb commented 7 months ago

Does it mean that right now MTK does not have an interface for discrete-time systems?

Yes, though I hope we can get around to getting something back next week. I think from other various discussions this has come up as a high priority. @AayushSabharwal let's discuss this and get something together.

@ChrisRackauckas Is there an issue/PR I could subscribe to in order to know the status of this?

ChrisRackauckas commented 7 months ago

https://github.com/SciML/ModelingToolkit.jl/pull/2507

SciML / ModelingToolkit.jl

Downstream Updates from MTK v9 #2483

Edit