Feedback before putting in the mass action update

[isaacsas]

I see a few ways I can add in a mass action direct method. My current code

1. Adds a MassActionJump type and updates JumpProblem and all the JumpSet stuff to handle this fourth type. This new jump is passed down as a parameter to aggregate and is stored in the JumpAggregation as a new field.
2. I then have a new DirectMA aggregator algorithm which supports both constant jumps (stored as function wrappers internally) and mass action jumps. It works with a mix of both, or either on their own. (So it would also be useful to people who won't have mass action jumps, but will have > ~10 constant jumps.) It conforms to the new SSA interface @alanderos91 put through yesterday.

I could just put this through as is. Alternatively, I could revise the existing direct method SSA jump routines @alanderos91 updated yesterday to support mass action jumps and function wrappers (the latter through dispatch on the time to next jump function as discussed in the other thread). However, that would mean for now that users don't really have any way to use the function wrappers since the current JumpSet processing code will always result in a tuple of constant jumps being passed into aggregate.

My instinct is to just submit the DirectMA aggregator algorithm as is, and if down the line someone wants to update the JumpSet code to allow other ways to pass in constant jumps (general iterable containers?) then we can look into merging Direct and DirectMA. That would be a good update anyways since the current tuple-based code may have issues with larger numbers of constant jumps... More generally, I will probably just stick with all future SSAs that support both mass action and constant jumps using function wrappers internally to store constant jumps. This will be a bit slower for small numbers of constant jumps, but faster once there are ~10 or more...

Any thoughts?

[alanderos91]

I agree that Direct and DirectMA should be merged. However, I have to own up to not fully understanding the distinction between constant and mass action jumps. The only reason for keeping things separate would be performance. Is this correct? If so, I think this boils down to whether there is a performance penalty for adding a new field to DirectJumpAggregation. The only place I see this cropping up is in building the JumpAggregation but that should be negligible.

It should be possible to get the function wrappers working by adding the aggregate method in your gist to DirectJumpAggregation that specializes on a new aggregator algorithm. The Direct aggregator algorithm ought to be renamed, and the new aggregate method will specialize on the DirectFunWrappers type. All that is left is replacing the helper function for the next reaction with your implementation.

This brings up the issue of how to name these aggregator algorithms, so that a user can appreciate the difference solve(jump_prob, DirectTuples()) and solve(jump_prob, DirectFunWrappers()).

In this way, slight variations on the basic Gillespie algorithm (e.g. FunctionWrappers, sparse/dense stoichiometric matrix representations...) can be implemented by providing (1) a new type that implements AbstractAggregatorAlgorithm in aggregators.jl and (2) a new aggregate method that feeds in the correct types into a type constructor.

[isaacsas]

The only reason for keeping things separate would be performance. Is this correct?

Yes that is correct. But it is a huge difference. For example, consider a system with 16 different A->B reactions (each having a different rate and ungrouped). Running the reaction 8000 times gives:

Direct method with tuples: 10.359 s (45989311 allocations: 1.03 GiB) Direct method with function wrappers 2.928 s (29098591 allocations: 451.39 MiB) Direct method with mass action rate function 146.208 ms (174991 allocations: 10.06 MiB)

As the number of reactions increases the function wrapper version gets much faster than tuples, but is still always much slower than the mass action jump version.

It should be possible to get the function wrappers working by adding the aggregate method in your gist to DirectJumpAggregation that specializes on a new aggregator algorithm.

That's what I did in the PR. Thanks for the suggestion!

This brings up the issue of how to name these aggregator algorithms

I used DirectManyJumps. If you have a better suggestion I'm happy to change it... I admit it's not great.

[ChrisRackauckas]

If so, I think this boils down to whether there is a performance penalty for adding a new field to DirectJumpAggregation. The only place I see this cropping up is in building the JumpAggregation but that should be negligible.

Yeah, it would be negligible.

[isaacsas]

Closed by 22