feat!: use SciMLStructures and add new MTKParameters

[AayushSabharwal]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[AayushSabharwal]

Also worth noting that this still scalarizes vector parameters. Not doing that won't be breaking, so it's not blocking the release.

[TorkelE]

You mean that parameters that are length-1 vectors become scalars? If so, would it be possible to not do that and let the user decide? I have a case which heavily depends on parameters that are length-1 vectors.

[YingboMa]

We probably should remove the # !s.

[AayushSabharwal]

Fixed in b20236e.

[YingboMa]

Should we keep the original p to test the fallback case?

[AayushSabharwal]

If the system has an IndexCache it will always use that for indexing, and expect that the parameter buffer is an MTKParameters. As such, collect(4:6) will be turned into MTKParameters(de, parameters(de) .=> collect(4:6))

[isaacsas]

Is this really the best way to do this? Doesn’t this mean every call is reallocating a vector of pairs? BifurcationKit will call this function many many times.

[ChrisRackauckas]

No, this is not a good way to do this. We need to do the interface a bit differently here, but I don't want that to block the v9. It should build the problem and use the generated functions from that.

[isaacsas]

I haven't reviewed these changes, but I hope that condition/affect functions that are generated will have similar performance to the current MTK and don't have such allocations (which would be really bad for jump models).

[isaacsas]

(given how often those functions are called)

[ChrisRackauckas]

There are no allocations, and the split form is type stable so it should be faster for cases that have heterogeneous parameters.

[isaacsas]

Awesome!

[AayushSabharwal]

I haven't reviewed these changes, but I hope that condition/affect functions that are generated will have similar performance to the current MTK and don't have such allocations (which would be really bad for jump models).

This is a BifurcationKit-only hack, which I'm removing in favor of just not using MTKParameters here. Everywhere else in the codebase, this sort of thing does not happen. MTKParameters is only created once for a problem, and all of the codegen leverages it.

[AayushSabharwal]

85e06f6

[AayushSabharwal]

You mean that parameters that are length-1 vectors become scalars? If so, would it be possible to not do that and let the user decide? I have a case which heavily depends on parameters that are length-1 vectors.

No, this only means that parameter values need to be specified as [p[1] => 1.0, ...], similar to how get_u0_p expects them to be. Indexing array parameters still returns an array:

julia> @parameters p[1:1]
julia> @variables x(t)
julia> @mtkbuild sys = ODESystem([D(x) ~ sum(p) * t], t)
julia> prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p[1] => 1.0])
julia> prob.ps[p]
1-element Vector{Float64}:
 1.0

This is something that will be changed, to make specifying [p => [1.0]] possible, but that probably needs some changes outside of just IndexCache/MTKParameters and isn't blocking for MTKv9, so I've pushed it down the road a bit in the interest of releasing sooner.

[ChrisRackauckas]

Open an issue on that.

[TorkelE]

You mean that parameters that are length-1 vectors become scalars? If so, would it be possible to not do that and let the user decide? I have a case which heavily depends on parameters that are length-1 vectors.

No, this only means that parameter values need to be specified as [p[1] => 1.0, ...], similar to how get_u0_p expects them to be. Indexing array parameters still returns an array:

julia> @parameters p[1:1]
julia> @variables x(t)
julia> @mtkbuild sys = ODESystem([D(x) ~ sum(p) * t], t)
julia> prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0), [p[1] => 1.0])
julia> prob.ps[p]
1-element Vector{Float64}:
 1.0

This is something that will be changed, to make specifying [p => [1.0]] possible, but that probably needs some changes outside of just IndexCache/MTKParameters and isn't blocking for MTKv9, so I've pushed it down the road a bit in the interest of releasing sooner.

That makes sense, thanks