fix: fix remake for parameters dependent on observed variables

[AayushSabharwal]

@isaacsas @TorkelE this should fix the remake issue you were having with conservation laws.

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.

[ChrisRackauckas]

why is this preferred? This would be a slower construction.

[AayushSabharwal]

Is it? I thought it would be faster since it doesn't create a new dictionary

[ChrisRackauckas]

oh I didn't see that it was extending an old one.

It's still adding elements one by one. That might amortize?

[AayushSabharwal]

julia> function oop(dict)
           return Dict{Any, Any}(k => 2v for (k, v) in dict)
       end
julia> function iip(dict)
           for (k, v) in dict
               dict[k] = 2v
           end
       end
julia> dict = Dict{Any, Any}(rand() => rand() for _ in 1:10000);
julia> oop(dict); iip(dict);
julia> @be oop(dict)
Benchmark: 25 samples with 1 evaluation
 min    2.599 ms (50092 allocs: 1.183 MiB)
 median 2.696 ms (50092 allocs: 1.183 MiB)
 mean   3.925 ms (50092 allocs: 1.183 MiB, 3.60% gc time)
 max    29.143 ms (50092 allocs: 1.183 MiB, 90.08% gc time)

julia> @be iip(dict)
Benchmark: 86 samples with 1 evaluation
 min    1.004 ms (20085 allocs: 313.828 KiB)
 median 1.052 ms (20085 allocs: 313.828 KiB)
 mean   1.140 ms (20085 allocs: 313.828 KiB, 1.77% gc time)
 max    5.962 ms (20085 allocs: 313.828 KiB, 81.95% gc time)

In-place seems faster

[ChrisRackauckas]

Okay, that's mainly because the generator isn't optimizing there 😅 . Both are slow, but it's fine.

[isaacsas]

Thanks for the quick fix!

I still think this is likely to lead to silent errors for many users if

remake(prob; u0 = [...])

and

remake(prob; u0 = [...], p = Dict())

are not equivalent. I strongly suspect most people will just use the former and not even realize they might need to use the latter. (Which as we see currently leads to silent correctness issues.)

Is there no way to make the latter the default behavior, or alternatively have a system opt into it? If it is made the default there could be a kwarg to turn it off for more advanced users who know their system doesn't need such behavior (and hence can get better performance avoiding it)? That feels like a safer convention from a correctness perspective.

[ChrisRackauckas]

Well not quite. It's specifically a non-symbolic remake, remake(prob, u0 = [x => 5]), is fine. remake(prob, u0 = [2.0,3.0], p = [3.0,2.0]) reverts back to the normal non-symbolic dispatch though, which then just makes those arrays prob.u0 and prob.p. That of course means it doesn't hit any of the initialization stuff etc... and I hate it. Indeed it's a bit of a footgun. I do think that if you have a symbolic system and you try to do that, you need to do some override boolean and we should error by default. I have seen too many people do non-symbolic things assuming an order etc. and it's just wrong.

There is a use case with optimization where, if an optimizer is choosing new parameters, you want to just remake(prob, p=newp) (and indeed with MTKParameters it's doing something a little special only replacing the tunables, but still), but we should definitely guard that with safety because I see too many people accidentally hit the non-symbolic dispatch thinking it means something that it doesn't. Basically, you should have to say "yes, I'm very sure that I understand the risks here, I know my parameter ordering, I know what SymbolicIndexingInterface, getsym, setsym etc. is, but I want to just put this array here" and we can go "sure, it sounds like you read and know the error with this, good luck".

[isaacsas]

Can one not detect if they are in the symbolic remake context, and in that case enforce that the two forms are equivalent (modulo an explict kwarg to opt out)?

[ChrisRackauckas]

That's what we do right now, without the keyword argument to make it explicit. Dict() just means it's an empty symbol map.

[isaacsas]

Yes, but I’m saying to make passing an empty map equivalent to not passing the kwarg at all. These are different now in the symbolic case.

[ChrisRackauckas]

That's what we do. That's why remake(prob; u0 = [2.0,3.0], p = Dict()) hits the symbolic dispatch while remake(prob; u0 = [2.0,3.0]) does not.

[AayushSabharwal]

I think @isaacsas means remake(prob; u0 = [x => 1.0]) is not equivalent to remake(prob; u0 = [x => 1.0], p = Dict()).

[isaacsas]

Yes, that is what I am saying should be changed.

[ChrisRackauckas]

Yes that should definitely by changed.

[AayushSabharwal]

Isn't it breaking, though?

[AayushSabharwal]

The reason it was implemented this way was (iirc) for the limited cases where symbolic remake worked prior to SII, that's what it did.

[AayushSabharwal]

Yeah, see https://github.com/SciML/SciMLBase.jl/blame/21e74609d18ba86f56fbac7ef4530566524d7504/src/remake.jl#L85-L92

[isaacsas]

This feels like it is a worthwhile reason for a breaking change.

[ChrisRackauckas]

Yeah but what I mean is, is there even any ODE that is solvable where u0 = [x => 2.0]? I would think that would just error in the solver since zero(u0) isn't defined or something. I can't think of a case you'd actual break.

[AayushSabharwal]

u0 = [x => 2.0] hits the symbolic path, which updates the appropriate index of the u0 vector and creates the initialization problem. However, it doesn't update any parameter values which have a default dependent on the new u0.

[ChrisRackauckas]

oh... that's a bug. It's got to be interpreted as a bug. I can't see someone actually relying on that behavior?

[AayushSabharwal]

I mean, it allows you to avoid the overhead of creating the symbolic map, substituting and calling remake_buffer while getting the same solve result because initialization.

[isaacsas]

I mean, it allows you to avoid the overhead of creating the symbolic map, substituting and calling remake_buffer while getting the same solve result because initialization.

But only for ODEs/DAEs. For all other problem types it is a potential bug no?

[isaacsas]

Is there no notion of an underlying dependency map that could be tested to decide if " creating the symbolic map, substituting and calling remake_buffer while getting the same solve result because initialization" could be bypassed?

[isaacsas]

Again though, it feels like the default here should be the safest option, with ways to opt out for users who know what they are doing and need better performance.

[ChrisRackauckas]

Yeah this seems to be trying to hard to save a little bit of performance. The setsym interface is required for full performance anyways, so this isn't necessarily fast, just a little not slower.

[isaacsas]

Why was this merged and released? The Catalyst downstream test failures that appear here are pretty clearly remake related. Several Catalyst test sets are now failing too, see everything I'm having to disable in: SciML/Catalyst.jl#1098

[ChrisRackauckas]

That wasn't clear...