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
136 stars 35 forks source link

function wrapper type instability fix #309

Closed isaacsas closed 1 year ago

isaacsas commented 1 year ago

@ChrisRackauckas what do you think about this approach? Basically, when not using tuples for affects we store them as affects!::Any (but they are concretely a Vector{Any}). When init is called the Vector{Any} that is stored is concretized into a vector of type-stable function wrappers using the integrator type as input. Then I add a manual dispatch when affects! isn't a tuple to handle pulling out this vector from the affects!::Any field.

Preliminary testing is that overall I'm seeing a half-percent or less speed decrease when using Direct with tuples (so still much faster than before the recent generated function approach I put in), but a ~17% or more speed up when using function wrappers. Moreover, the memory allocation issues go away, so I think we now have type-stability outside the highest level call into the jump aggregator callback.

If this approach looks reasonable I'll cleanup and refactor the code from Direct to work with all the SSAs.

One nice thing about the setup here is that it wouldn't be too bad to support affect! tuples for any method now, which can give a nice speed up (Direct vs. DirectFW on a small system is still like a 40% speed difference). We might want to add a similar manual dispatch for rate functions so that we can also support them as tuples with other SSAs. I’m thinking we could then switch based on the number of jumps to automatically use tuples vs. function wrappers for any method.

coveralls commented 1 year ago

Pull Request Test Coverage Report for Build 4745185989

Warning: This coverage report may be inaccurate.

This pull request's base commit is no longer the HEAD commit of its target branch. This means it includes changes from outside the original pull request, including, potentially, unrelated coverage changes.

Details


Changes Missing Coverage Covered Lines Changed/Added Lines %
src/aggregators/ssajump.jl 33 34 97.06%
<!-- Total: 82 83 98.8% -->
Files with Coverage Reduction New Missed Lines %
src/aggregators/prioritytable.jl 1 86.86%
src/extended_jump_array.jl 1 27.03%
src/jumps.jl 1 84.46%
src/coupling.jl 2 78.18%
src/spatial/spatial_massaction_jump.jl 3 85.45%
src/problem.jl 51 64.4%
<!-- Total: 59 -->
Totals Coverage Status
Change from base Build 4549721708: 0.3%
Covered Lines: 2136
Relevant Lines: 2420

💛 - Coveralls
ChrisRackauckas commented 1 year ago

This seems reasonable. Though I do think that we should be moving to a scheme where anything cached is in the init phase and when restructuring I think the Any can be avoided by directly declaring the type of the integrator. But that may never happen, so this is very good in the meantime.

isaacsas commented 1 year ago

I'd like to switch to such a setup too as that would also get rid of @gaurav-arya's issue, however, I'm not in favor of doing so until we have a nice way to rerun the same simulation many times for sampling without a user having to drop into the integrator interface. (i.e. I think EnsembleProblems would need updating to also then use the integrator interface and reinit between runs on a thread.)

But in any case I think we'd still have this issue. Fundamentally, the problem is a circular reference in that the integrator stores the affects, but the wrapped affects need to know the type of the integrator. We'd need some other object to cache the JumpAggregation, which would have to be created after the integrator to avoid this.