fix: fix remake with u0 dependent on Symbol parameter

[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]

I have verified that this passes the tests in https://github.com/SciML/Catalyst.jl/blob/0b4092e7b6992440b07fc0d58516ad2340b0f88e/test/upstream/mtk_problem_inputs.jl#L1-L176

[AayushSabharwal]

@isaacsas @TorkelE what are your opinions on the new API?

[isaacsas]

So the default here (missing_is_symbolic=true) is the behavior we discussed with Chris (i.e. not having to pass p = Dict()), but setting it to false goes back to the old behavior? If so that sounds fine to me.

I am running the full set of Catalyst tests against this PR locally and will report back. (The stability computation tests were also failing after the remake changes yesterday.)

[isaacsas]

I'm now getting a failure: https://github.com/SciML/Catalyst.jl/blob/0b4092e7b6992440b07fc0d58516ad2340b0f88e/test/network_analysis/conservation_laws.jl#L265 locally.

I'll upcap SciMLBase on master so this PR can test against it.

[isaacsas]

OK, it is uncapped so hopefully downstream here will just run against it. It was passing tests fine last night when restricted to 2.71.1 so any failures are likely due to this and/or the previous PR.

If we need to update a test due to incorrectly assumed behavior please let me know. I probably won't have any time free today, but can take a look tonight or tomorrow morning.

[AayushSabharwal]

So the main failure here is MTK and it's a very annoying one

    @variables x(t) y(t)
    @parameters p q
    @mtkbuild sys = ODESystem(
        [D(x) ~ x, p ~ x + y], t; defaults = [p => missing], guesses = [x => 0.0, p => 0.0])
    prob = ODEProblem(sys, [x => 1.0, y => 1.0], (0.0, 1.0))
    prob2 = remake(prob; u0 = [x => 1.0, y => 2x + exp(t)])

this is the code that gets run. The u0 provided to remake promotes to a Vector{Pair{Num, Num}}, so it tries to create an ODEProblem where u0 == Num[1.0]. The problem is, SciMLBase doesn't know about unwrapping Num to get the number inside.

[AayushSabharwal]

The best solution I have right now is to make SII.symbolic_evaluate (which is used to allow SciMLBase to perform symbolic substitution) unwrap everything it returns. Not the prettiest way to do things, though.

[isaacsas]

Yeah, that sounds a bit messy. But if it fixes this issue it is a starting point at least that can hopefully be improved in the future.

[AayushSabharwal]

JuliaSymbolics/Symbolics.jl#1329

[AayushSabharwal]

https://github.com/SciML/Catalyst.jl/blob/0b4092e7b6992440b07fc0d58516ad2340b0f88e/test/network_analysis/conservation_laws.jl#L265

This failure is also annoying. How remake works is it builds dictionaries mapping unknowns/parameters to values. It ends up with the dictionaries u0 = [X1 => 1.0] and p = [k1 => 0.3, k2 => 0.4, Gamma[1] => X1 + X2]. Since X2 is an observed, it doesn't show up in the u0 dictionary and hence Gamma[1] can't get a value. If I make it so u0 is also populated with observed variables, we get X2 => Gamma[1] - X1 which leads to a circular substitution.

My solution would be to make it populate observed variables in u0, then if we end up in a case with circular rules (basically if one of the values in the map doesn't resolve to a non-symbolic quantity) use SII.getsym to get the existing value and substitute that in. The downside is that this requires us hitting maxiters on all circular substitutions. It can also lead to ordering issues, where we might get different results if it finds X2 to be circular first and substitutes its existing value versus finding Gamma first.

[AayushSabharwal]

I could also not put observed variables in the u0 map, and if we don't have enough information to symbolically evaluate a value then use getu. This has the advantage of not hitting maxiters but could potentially run into similar ordering issues and have incorrect values, since getu won't take into consideration the newly provided values.

On a sidenote, this is why we have initialization systems 😅 regardless of what remake does here, initialization will do the right thing.

[ChrisRackauckas]

My solution would be to make it populate observed variables in u0, then if we end up in a case with circular rules (basically if one of the values in the map doesn't resolve to a non-symbolic quantity) use SII.getsym to get the existing value and substitute that in. The downside is that this requires us hitting maxiters on all circular substitutions. It can also lead to ordering issues, where we might get different results if it finds X2 to be circular first and substitutes its existing value versus finding Gamma first.

This is what I keep saying. It's inherently an implicit system the moment any simplification is made. Cycling substitutions to a fixed point is just fixed point iteration of the nonlinear system in a slow way. We should not do this since it's not robust and it can be really really slow (like hours 😅). I know someone will look at this and say "well I as a human would first substitute this into this and then get this, then calculate..." I know you would! How do we come up with that algorithm? That algorithm for the proper order to do the substations in order to best resolve the system is structural_simplify!!!

So this is why I'm saying it's just not a fruitful path. At its core, we're saying, can we avoid structural simplify and the initialization system by instead substituting the right values? Okay, how do we substitute that right? Well, we just need to detect where the cycles are in the graph, perform tearing, ... congrats, that's just structural simplify again with extra steps 😅.

If you do want to enforce some kind of partial solve correctness, then what you could do is just have some way of doing a no-solve solve on the initialization problem (i.e. an algorithm that takes u0 and calls that the solution and gives you back a NonlinearSolution to play with, which should have the initialization analytical solutions in the observed due to structural simplification) so that you get the observed setup on the u0 and that would tell you how to resolve the linear things that can be analytically solved for. In Catalyst's case this seems like a good idea because the conservation laws, in theory and if supplied correctly, should always show up in the analytical solution. In practice though, this would not check for inconsistency, so you'd need to assert that any relationship pulled is not an unknown of the initialization system.

[AayushSabharwal]

The problematic part of this PR is the direct result of the decision at the end of #832 to re-evaluate symbolic defaults... I don't see how we can have it both ways.

[AayushSabharwal]

It's also why I mentioned that this is why initialization exists 😅 we're trying too hard to make remake calculate the final u0/p when this is exactly what OverrideInit does. remake(; u0 = [x => 1, y => 2], p = Dict()) should not be equivalent to remake(; u0 = [x => 1, y => 2]). remake can do exactly what the user asks, no more no less, and let MTK figure out the implications on the rest of the values in the system.

[isaacsas]

Sorry this is such a mess. From the perspective of someone that doesn't know the difference between what remake is doing vs. *Problem, it is hard to see why this has always worked for the latter but can't be made to work easily for the former.

If you all want to go back to the old remake behavior we can just proceed with wrapping it in Catalyst. This seems like it would be a lot less complicated for everyone since we are perfectly happy to special case the treatment of these constants in Catalyst, exploiting that they

1. only appear via a linear relationship with the initial conditions of unknowns and observables
2. are assumed to either be defined via this relationship or to be set explicitly

We can use setp to just update gamma, and setu for any unknown that might need updating, after remake finishes directly (since that doesn't impact other parameters/unknowns right?).

[AayushSabharwal]

The core of the problem is that the *Problem constructors live in MTK, and remake is a SciMLBase-level abstraction. MTK has a lot more information to work with and can make more assumptions to allow handling these cases. Also, remake is fundamentally different in that you can specify a new value for a subset of the variables and only those will be updated. You can't do that with problem constructors, since that is an underspecification of the problem. MTK handles this by letting remake do its thing and solving for the actual values during initialization since in general there is no way to explicitly evaluate everything.

[isaacsas]

But every problem now wraps a SciMLFunction that wraps a system when coming from MTK right? Can't this then be dispatched to have a custom MTK remake instead of trying to make one remake that works for everything? That would also have the advantage of MTK-related updates to remake no longer needing SciMLBase modifications.

I would think there may be other cases in the future where we'd want to exploit information that only an MTK system might have to improve remake anyways.

In any case, it sounds like this has become too complicated to handle via remake currently, so we should just go back and have our own Catalyst wrapper than can exploit the extra information we have about this (linearity of the relationship, each gamma corresponding to a specifical observed, etc).

[ChrisRackauckas]

We already have been speicalizing remake for awhile for other reasons https://github.com/SciML/SciMLBase.jl/blob/master/src/remake.jl#L109-L204

[AayushSabharwal]

But every problem now wraps a SciMLFunction that wraps a system when coming from MTK right?

Not every problem. People still use the non-MTK interface. We also can't assume it comes from MTK, since people can both build their own systems or have an alternative DSL that hooks into this interface. I know of at least one such non-MTK project, though it isn't public. MTK also doesn't codegen to every SciMLFunction, and some SciMLFunctions don't wrap a system.

MTK also can't solve the circular substitute problem above, its solution is to create the initialization system. The problem is that the intent is not immediately clear from the values passed to remake because remake allows you to pass an incomplete problem specification. You both want to retain old values that you haven't touched, and update ones that depend on the variables you provided new values for. The algorithm that figures out exactly how to update everything is structural_simplify + a nonlinear solve, as Chris mentioned, by treating the new values as constraints during initialization that override any old constraints for those values and merge with the existing ones.

[isaacsas]

@AayushSabharwal, a problem passed to remake, and which stores a SciMLFunction that stores an MTK.AbstractSystem, can not be assumed to be from an MTK-based library? Isn't that what a dispatch would handle?

People who store non-MTK.AbstractSystems in the SciMLFunction, or no SciMLFunction at all, would not get routed to the dispatch, so still use what is in SciMLBase.

In any case, that is tangential to the bigger issue. The conclusions I've come to from all our discussion is this can't currently be done at the level of remake, requires initialization with how the parameter/u0 modification process is being envisioned here, but initialization is only supported for ODEs currently (and even for ODEs, we might have to tell users to do non-standard stuff with remake or problem construction like passing p = Dict() or guesses to get initialization to work?).

So it seems like the most immediate solution is adding a wrapper in Catalyst where we can exploit the additional knowledge we have there. I'm happy to consider alternative solutions, but I'd like a solution we can implement immediately and include in Catalyst 15 (which I want to release soon as it has at least one dramatic performance fix). Something that requires waiting an indeterminate amount of time for functionality to be added across all possible solvers Catalyst targets isn't a useful solution right now (but of course would be great long-term).

[AayushSabharwal]

a problem passed to remake, and which stores a SciMLFunction that stores an MTK.AbstractSystem, can not be assumed to be from an MTK-based library? Isn't that what a dispatch would handle?

Sure, but that dispatch won't be able to do much more than what we do here in SciMLBase anyway for the reasons mentioned previously. In fact, I'm not sure what new thing it'll do anyway.

(and even for ODEs, we might have to tell users to do non-standard stuff with remake or problem construction like passing p = Dict() or guesses to get initialization to work?).

Passing p = Dict() has the effect of updating the ODEProblem. Regardless of whether you do this, initialization will update the dependent values appropriately.

So it seems like the most immediate solution is adding a wrapper in Catalyst where we can exploit the additional knowledge we have there.

I think so, yeah. In Catalyst you assume that the problem always contains the final value of u0 and p, which is not a valid assumption anywhere else. I'm not sure how else to resolve this. It's worth noting, though, that if you add a rebuild which passes p = Dict() to remake, you can still run into the circular dependency issue.

[isaacsas]

Passing p = Dict() has the effect of updating the ODEProblem. Regardless of whether you do this, initialization will update the dependent values appropriately.

Except the many examples I put in the rebuild issue show that is not currently happening? The only time solutions were correct was when building a new problem directly (i.e. calling ODEProblem), or when explicitly passing a value for gamma to remake. The ODE solutions were wrong with every variant of remake I tested that didn't give an explicit value to gamma.

It's worth noting, though, that if you add a rebuild which passes p = Dict() to remake, you can still run into the circular dependency issue.

Sure, but that is a potential issue even when not using the conservation law functionality (i.e. a user could define initial conditions and parameters to introduce a circular dependency). Right now, that will impact all non-ODE problems anyways since they also don't include initialization.

[AayushSabharwal]

Except the many examples I put in the rebuild issue show that is not currently happening? The only time solutions were correct was when building a new problem directly (i.e. calling ODEProblem). The ODE solutions were wrong with every variant of remake I tested except where I explicitly passed a value for gamma.

I don't think those cases had a guess for Gamma? Parameters have to opt in to initialization (https://docs.sciml.ai/ModelingToolkit/stable/tutorials/initialization/#Initialization-of-parameters) unlike unknowns since usually a missing parameter value means a fault in the model.

[AayushSabharwal]

Sure, but that is a potential issue even when not using the conservation law functionality (i.e. a user could define initial conditions and parameters to introduce a circular dependency). Right now, that will impact all non-ODE problems anyways since they also don't include initialization.

That is correct.