SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
135 stars 35 forks source link

Ambiguity fix #385

Closed isaacsas closed 6 months ago

isaacsas commented 6 months ago

Let's see if tests pass with this.

isaacsas commented 6 months ago

@ChrisRackauckas were there more plotting changes recently? The doc build crash here is plotting a solution object that previously worked.

ChrisRackauckas commented 6 months ago

There was just the major change https://github.com/SciML/SciMLBase.jl/pull/572. The syntax of the plots isn't changed at all, just the internals. The internals got rid of everything special about plotting and now just uses the solution interface. This revealed that we had an issue with idxs with FunctionMap https://github.com/SciML/OrdinaryDiffEq.jl/pull/2099.

ChrisRackauckas commented 6 months ago

https://github.com/SciML/JumpProcesses.jl/actions/runs/7360900012/job/20037606665?pr=385#step:6:671 parameter indexing test needs to be updated but @AayushSabharwal is syms still supported? I thought that was effectively removed.

isaacsas commented 6 months ago

I think @TorkelE had added parameter indexing stuff here a while back. I guess that needs the updates SciMLbase got too?

TorkelE commented 6 months ago

I'm confused. Has the new MTK version with the new indexing been pushed and is that what we are using? Or has it not, by why is then problem indexing affected?

isaacsas commented 6 months ago

This is about SciMLBase not MTK.

TorkelE commented 6 months ago

SUre, but it is all tied together, right?

My understanding was that a big MTK update was coming that would make indexing easier (MTK v9). Obviously this would change internal stuff as well. But has this been rolled out already? Or it hasn't? Or just in some places? And if big changes have been rolled out in SciML base, why wasn't it made a breaking release? Or this have nothing to do with the new big changes, but it seems like Chris is talking about them.

Basically, it would be useful to learn what is actually planned to happen where and when, and have some kind of announcement before major things happens. But maybe this wasn't a major change, but then exactly what it is that have changed?

ChrisRackauckas commented 6 months ago

The RecursiveArrayTools v3, and SymbolicIndexingInterface v0.3 roll outs have all already happened. See https://docs.sciml.ai/SymbolicIndexingInterface/dev/usage/ for examples on how it's used.

And if big changes have been rolled out in SciML base, why wasn't it made a breaking release?

I mentioned many times before that indexing from Symbol isn't something it truly supported. There were some deprecated functions that were still around and could be hacked to be made used, but they (a) weren't there in all interfaces and (b) aren't always correct. And there's not really a way to make them correct. That's why I said Symbol indexing isn't supported even when this PR was going in. It's fine to track it but that's why I've been heavily against putting any of that into documentation, since we knew it was going away (as if it was ever there, again, it was not very stable and usable beyond simple cases to begin with). We kept it out of the documentation so we could easily remove it.

Basically, it would be useful to learn what is actually planned to happen where and when, and have some kind of announcement before major things happens. But maybe this wasn't a major change, but then exactly what it is that have changed?

We had 2 of the core packages do a major update, and hundreds of PRs to go along with it.

isaacsas commented 6 months ago

@ChrisRackauckas should I just remove the getindex/setindex stuff on JumpProblems (and associated tests)?

TorkelE commented 6 months ago

Last time we chatted you said that it was being discussed, and you generally was no in favour of it, and we both stated why we think differently. Both since then and before I have been in contact with the people handling the implementation, and I have been reassured that Symbol indexing is being accounted for and will be supported.

Just to demonstrate why this is important to end-users, here is an example where you create a bifurcation diagram, both using Symbol indexing, and system. indexing. This is not a contrived example, it is a real one from one of my actual research papers. There is no way anyone can tell me that the new notation is an acceptable alternative. I know alternatives have been discussed, but I don't think anyone have actually looked at even a proposed notation.

Using `Symbol notation that thought we would support:

narula_system = @reaction_network begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    v0*((1+F*σB)/(K+σB)),     ∅ ⟶ σB
    λW*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ w
    λV*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ v
end

p_vals = [:kBw => 3600, :kDw => 18, :kD => 18, :kK1 => 36, :kK2 => 36, :kP => 180,:kDeg => 0.7, 
          :kB1 => 3600, :kB2 => 3600, :kB3 => 3600, :kB4 => 1800, :kB5 => 3600, 
          :kD1 => 18, :kD2 => 18, :kD3 => 18, :kD4 => 1800, :kD5 => 18, 
          :v0 => 0.4, :F => 30, :K => 0.2, :λW => 4, :λV => 4.5, :pInit => 0.001, :pStress => 0.4]

u_guess = [:w => 1.0, :w2 => 1.0, :w2v => 1.0, :v => 1.0,:w2v2 => 1.0,  :w2σB => 1.0, :σB => 1.0,
          :vP => 1.0, :phos => 1.0, :vPp => 1.0]
u0 = [:phos => 0.4, :vPp => 0.0]
bprob = BifurcationProblem(narula_system, u_guess, p_vals, :kK2; plot_var = :σB, u0)

Using the system. notation:

narula_system = @reaction_network begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    v0*((1+F*σB)/(K+σB)),     ∅ ⟶ σB
    λW*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ w
    λV*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ v
end

p_vals = [narula_system.kBw => 3600, narula_system.kDw => 18, narula_system.kD => 18, narula_system.kK1 => 36, 
          narula_system.kK2 => 36, narula_system.kP => 180,narula_system.kDeg => 0.7, 
          narula_system.kB1 => 3600, narula_system.kB2 => 3600, narula_system.kB3 => 3600, narula_system.kB4 => 1800, 
          narula_system.kB5 => 3600,  narula_system.kD1 => 18, narula_system.kD2 => 18, narula_system.kD3 => 18, narula_system.kD4 => 1800, 
          narula_system.kD5 => 18,  narula_system.v0 => 0.4, narula_system.F => 30, narula_system.K => 0.2, narula_system.λW => 4, 
          narula_system.λV => 4.5, narula_system.pInit => 0.001, narula_system.pStress => 0.4]

u_guess = [narula_system.w => 1.0, narula_system.w2 => 1.0, narula_system.w2v => 1.0, narula_system.v => 1.0,
           narula_system.w2v2 => 1.0,  narula_system.w2σB => 1.0, narula_system.σB => 1.0,
          narula_system.vP => 1.0, narula_system.phos => 1.0, narula_system.vPp => 1.0]
u0 = [narula_system.phos => 0.4, narula_system.vPp => 0.0]
bprob = BifurcationProblem(narula_system, u_guess, p_vals, narula_system.kK2; plot_var = narula_system.σB, u0)
isaacsas commented 6 months ago

@TorkelE the latter issue could be fixed by a @valuemap macro as I had suggested:

p_vals = @valuemap narula_system [kBw => 3600, kDw => 18, ...]

(I like explicitly using the desired container type here since then one could potentially give tuples or other input type arrays to store the mapping.)

But I understand that such functionality doesn't exist currently.

ChrisRackauckas commented 6 months ago

You're talking about something different. Did you read the link that I just sent? https://docs.sciml.ai/SymbolicIndexingInterface/dev/usage/.

If necessary, Symbols can be used to refer to variables. This is only valid for symbolic variables for which hasname returns true. The Symbol used must match the one returned by getname for the variable.

What I am referring to is the fact that DiscreteFunction has an undocumented and not public f.syms that is relied on in this test here. There has been no guarantee that is staying. Meanwhile, SII v0,3 makes it very explicit the cases in which indexing using a symbol on a well-defined symbolic structure will work, i.e. when it matches hasname / getname. This is still not a great idea for hierarchical systems (because . is not a supported symbol in Julia), but there is a guarantee that will work.

Right now, f.syms works by building a sys::SymbolCache which implements the SII. However, to simplify the interfaces what we want to do is just tell DSLs to use f.sys, document that f.sys exists, and then there's only one way by which the symbolic indexing is defined and everyone is on the same interface.

From this, other DSLs can implement how to symbolically interface with their systems by either building a SymbolCache and passing that as f.sys or directly defining the interface functions on their system type (like ModelingToolkit does). With that, it's now completed documented what's going on, there's a single interface, and MTK isn't special in any sense, and also f.syms isn't special.

TorkelE commented 6 months ago

You are right Chris, I had another look and get it better, my mistake.

Seems the issue is specific to DiscreteFunction then. I guess long term it would be good to bring it up to date with the rest of MTK, but I will leave to Sam to decide what to do for now.

Also, do SII allow Symbol interfacing for hierarchical systems then, and does not by default throw an error?

Finally, does this mean that the big indexing change has happened? Whatever the outcome I have been looking forward to that one. having this in place has been the last piece for needed for the comprehensive upgrade of the Catalyst docs, looking forward to be able to work on all parts of that one!

TorkelE commented 6 months ago

@isaacsas Yes, that was the one I alluded to. I think there are some stuff that needs to be ironed out.

E.g. in some cases you want to provide not a map, but a list of e.g. parameters or expressions (e.g. funcs_to_check in SI). Should you have a separate macro, or the same one, which also works one these, like:

known_ps = @make_symbolic system [p, d, k]
funcs_to_check = @make_symbolic system [X1+x2+X3, Y + 2*Y2Z]
isaacsas commented 6 months ago

That seems sensible to me.

ChrisRackauckas commented 6 months ago

Finally, does this mean that the big indexing change has happened? Whatever the outcome I have been looking forward to that one. having this in place has been the last piece for needed for the comprehensive upgrade of the Catalyst docs, looking forward to be able to work on all parts of that one!

It has already happened. We're cleaning up a few test errors that we're finding here and there, but I think it's all done now? This might be the last one.

Also, do SII allow Symbol interfacing for hierarchical systems then, and does not by default throw an error?

That's not really up to SII. SII states that in order for fast matching to occur, any DSL needs to define unique names for all of its symbols. getname and hasname are defined by SII https://docs.sciml.ai/SymbolicIndexingInterface/dev/api/#SymbolicIndexingInterface.hasname . DSLs with hierarchical models then just need to make sure their hierarchical naming scheme gives unique names for each symbol. However they do that, it will work and those symbols will then be able to be used as indexes.

The reason why that's not great still though is that :x.y is simply not a valid Julia Symbol, and so any hierarchical scheme can come up with a naming but arguably the most sensible naming is not one that can be used. This is why ModelingToolkit uses the unicode act₊vol₁₊x, which is pretty bad. If I think of something better for the MTK v9 then that's something that will change, I just have never found a better solution.

isaacsas commented 6 months ago

I just posted an issue on SymbolicIndexingInterface. It seems like there would be a benefit to having a formalized interface for getting a display representation of a name (as a string), which can be distinct from the Symbol-based indexing interface. This could be used then by show and plot recipes.

ChrisRackauckas commented 6 months ago

Is there more to that then string(getname(sol,sym))?

isaacsas commented 6 months ago

Yes, it could handle converting unicode plus to a normal dot for MTK systems.

ChrisRackauckas commented 6 months ago

interesting idea. I'm not opposed.

isaacsas commented 6 months ago

If we do that correctly, users would only encounter unicode plus when accessing variables/parameters via flattened systems directly, but all their display output would look nice. (And if people are trained to access variables/parameters via the pre-flattened base system then they might never have to see the unicode plus and it would just be internal.)

ChrisRackauckas commented 6 months ago

Yeah, I could see that. @YingboMa would you be opposed to having the printing show the .?

isaacsas commented 6 months ago

It could always just be used for plots and Latexify, but those are both cases where the unicode plus is definitely not desirable.

ChrisRackauckas commented 6 months ago

I heavily agree there. At that point though, we might as well also remove it from display/show.

TorkelE commented 6 months ago

Is MTK v9 something imminent (as in there is a plan to release when some stuff gets finished and we more or less know what changes it will bring), or still off in the hard to forsee future?

ChrisRackauckas commented 6 months ago

It's very imminent. https://github.com/SciML/ModelingToolkit.jl/issues/2262

ChrisRackauckas commented 6 months ago

I want it done by next week before the class starts.

TorkelE commented 6 months ago

Makes sense to have it ready for the class. Looking forward to the release!

isaacsas commented 6 months ago

@ChrisRackauckas will having heterogenous tuple inputs (i.e. floats and ints) be supported again? There needs to be a way to perfomantly have mixed parameter types.

ChrisRackauckas commented 6 months ago

As parameters? Yes. It's still implemented it's just off by default. That's getting turned back on.

isaacsas commented 6 months ago

Ok, I’d like to get the stuff that broke in catalyst working again. We had to basically add a comment that such support was broken due to a bug in our docs, and have had people asking why their tuples are getting converted to union type vectors that give performance warnings from ScoMLBase.

ChrisRackauckas commented 6 months ago

Yeah sorry, it's a bit of a bumpy month, but I'm trying to get everything onto standard interfaces so it's easier to maintain it all going forward. No secret side channels, closing all of those 😅

isaacsas commented 6 months ago

I'm all in favor of standard interfaces! I had wanted to actually chat with you about dev practices at some point for when one is making big interface changes (whether non-breaking or user facing and breaking).

The current approach is to merge the incremental dev work for such efforts into master if it doesn't appear breaking, which means changes get released piecemeal before the overall effort is completely finished. I was wondering if more ambitious rewrites within core packages like SciMLBase or MTK should perhaps go via PRs to feature branches. When they are complete an announcement could be made to test against them (basically an alpha/beta type release), and only after a chance for people to test would they then get merged into master and released. I feel like this would have (perhaps) avoided some of the bugs that slipped through in the past semester. (Note, for more minor changes, adding features that don't change existing (internal) API functions, or bug fixes I wouldn't think anything needs changing. This is more for things that make big internal changes that have a higher likelihood of inadvertently introducing changes/bugs that impact downstream packages.)

I admitedly have little experience on dev practices for larger software packages, but it seemed like maybe SciML is starting to reach the size where the most core libraries that many things depend might benefit from a modified dev/release workflow.

ChrisRackauckas commented 6 months ago

We did hold this on PR branches for a few months while doing a lot of downstream testing. There were about 12 PRs or something like that prepared and a thread with 500 replies of multiple people testing out a bunch of things. Most tests passed, but there was still quite a bit of clean up to do. Part of that is because v1.10 happened to land right in the middle of the rollout. Also the RecursiveArrayTools v3 happening at the same time didn't make it easier (though it was necessary because that's where the indexing break actually took place, and thus changing the AbstrsctVoA to no longer be a AbstractArray because of the linear indexing had been waiting for a long time and we finally had a good chance to do it). Changing very low level interfaces is just something that takes a significant amount of work.

prbzrg commented 6 months ago

Does this fix https://github.com/SciML/JumpProcesses.jl/issues/384 ?

isaacsas commented 6 months ago

Your comment said this was the missing dispatch.

isaacsas commented 6 months ago

Is there still an issue?

prbzrg commented 6 months ago

Yes, by [ccbc3e58] JumpProcesses v9.10.1, the report didn't change. Is it the latest version?

isaacsas commented 6 months ago

Can you try with master? We haven’t made a release with this fix yet.

prbzrg commented 6 months ago

I think the best way to be sure is to set up Aqua test in CI.

isaacsas commented 6 months ago

@ChrisRackauckas @TorkelE I don't know about the internal field f.syms, but it is the case that both syms and paramsyms are documented as keywords for DiscreteFunction that support symbol inputs, see

https://docs.sciml.ai/DiffEqDocs/stable/types/discrete_types/#SciMLBase.DiscreteFunction

So it would be nice to ensure they keep working within SciMLBase 2.X (perhaps by auto-constructing the SymbolCache now if they are used?).

ChrisRackauckas commented 6 months ago

That should've never gotten documented. If anything, we'll just need to do a fast SciMLBase v3 and merge without looking at CI just to get it out. There's no way we have the manpower to keep that thing around.

TorkelE commented 6 months ago

Basically DiscreteFunctions are getting discontinued? Is the idea that they are pretty much covered by DynamicalSystems.jl anyway?

ChrisRackauckas commented 6 months ago

No, DiscreteFunction isn't going anywhere. If anyone is saying that it's not true. It's f.syms that's not staying. It wasn't intended to be documented or used, and it only barely works in some scenarios.

TorkelE commented 6 months ago

Got it, good to know. There was something in https://github.com/SciML/ModelingToolkit.jl/issues/2262 about removing discrete time stuff, so got a bit uncertain exactly what was happening.

ChrisRackauckas commented 6 months ago

That's a different kind of discrete stuff. It's a form of discrete+continuous model that's used in control systems, and it's removing an older implementation that was deprecated for a newer one years ago (and it was undocumented). The newer tested and documented one would serve as a nice building block for making hybrid solvers (see https://github.com/SciML/ModelingToolkit.jl/pull/2410 for some details on how to use it), but the old one was more of a test run.