Merge pull request SciML#276 from gzagatti/queue-method-ii

`Coevolve` aggregator for `VariableRateJump`

Author: isaacsas

HISTORY.md

@@ -0,0 +1,9 @@
+# Breaking updates and feature summaries across releases
+
+## JumpProcesses unreleased (master branch)
+- Support for "bounded" `VariableRateJump`s that can be used with the `Coevolve`
+ aggregator for faster simulation of jump processes with time-dependent rates.
+ In particular, if all `VariableRateJump`s in a pure-jump system are bounded one
+ can use `Coevolve` with `SSAStepper` for better performance. See the
+ documentation, particularly the first and second tutorials, for details on
+ defining and using bounded `VariableRateJump`s.

docs/make.jl

@@ -1,7 +1,9 @@
 using Documenter, JumpProcesses
 
-cp("./docs/Manifest.toml", "./docs/src/assets/Manifest.toml", force = true)
-cp("./docs/Project.toml", "./docs/src/assets/Project.toml", force = true)
+docpath = Base.source_dir()
+assetpath = joinpath(docpath, "src", "assets")
+cp(joinpath(docpath, "Manifest.toml"), joinpath(assetpath, "Manifest.toml"), force = true)
+cp(joinpath(docpath, "Project.toml"), joinpath(assetpath, "Project.toml"), force = true)
 
 include("pages.jl")
 

docs/src/api.md

@@ -15,13 +15,16 @@ reset_aggregated_jumps!
 ConstantRateJump
 MassActionJump
 VariableRateJump
+RegularJump
 JumpSet
 ```
 
 ## Aggregators
 Aggregators are the underlying algorithms used for sampling
-[`MassActionJump`](@ref)s and [`ConstantRateJump`](@ref)s.
+[`ConstantRateJump`](@ref)s, [`MassActionJump`](@ref)s, and
+[`VariableRateJump`](@ref)s.
 ```@docs
+Coevolve
 Direct
 DirectCR
 FRM
@@ -36,4 +39,4 @@ SortingDirect
 ```@docs
 ExtendedJumpArray
 SSAIntegrator
-```
+```

docs/src/faq.md

@@ -1,16 +1,20 @@
 # FAQ
 
 ## My simulation is really slow and/or using a lot of memory, what can I do?
-To reduce memory use, use `save_positions=(false,false)` in the `JumpProblem`
-constructor as described [earlier](@ref save_positions_docs) to turn off saving
-the system state before and after every jump. Combined with use of `saveat` in
-the call to `solve` this can dramatically reduce memory usage.
+Exact methods simulate every jump, and by default save the state before and
+after each jump. To reduce memory use, use `save_positions = (false, false)` in
+the `JumpProblem` constructor as described [earlier](@ref save_positions_docs)
+to turn off saving the system state before and after every jump. Combined with
+use of `saveat` in the call to `solve`, to specify the specific times at which
+to save the state, this can dramatically reduce memory usage.
 
 While `Direct` is often fastest for systems with 10 or less `ConstantRateJump`s
-or `MassActionJump`s, if your system has many jumps or one jump occurs most
-frequently, other stochastic simulation algorithms may be faster. See [Constant
-Rate Jump Aggregators](@ref) and the subsequent sections there for guidance on
-choosing different SSAs (called aggregators in JumpProcesses).
+and/or `MassActionJump`s, if your system has many jumps or one jump occurs most
+frequently, other stochastic simulation algorithms may be faster. See [Jump
+Aggregators for Exact Simulation](@ref) and the subsequent sections there for
+guidance on choosing different SSAs (called aggregators in JumpProcesses). For
+systems with bounded `VariableRateJump`s using `Coevolve` with `SSAStepper`
+instead of an ODE/SDE time stepper can give a significant performance boost.
 
 ## When running many consecutive simulations, for example within an `EnsembleProblem` or loop, how can I update `JumpProblem`s?
 
@@ -22,8 +26,9 @@ internal aggregators for each new parameter value or initial condition.
 ## How can I define collections of many different jumps and pass them to `JumpProblem`?
 
 We can use `JumpSet`s to collect jumps together, and then pass them into
-`JumpProblem`s directly. For example, using the `MassActionJump` and
-`ConstantRateJump` defined earlier we can write
+`JumpProblem`s directly. For example, using a `MassActionJump` and
+`ConstantRateJump` defined in the [second tutorial](@ref ssa_tutorial), we can
+write
 
 ```julia
 jset = JumpSet(mass_act_jump, birth_jump)
@@ -42,8 +47,8 @@ vj1 = VariableRateJump(rate3, affect3!)
 vj2 = VariableRateJump(rate4, affect4!)
 vjtuple = (vj1, vj2)
 
-jset = JumpSet(; constant_jumps=cjvec, variable_jumps=vjtuple,
-                 massaction_jumps=mass_act_jump)
+jset = JumpSet(; constant_jumps = cjvec, variable_jumps = vjtuple,
+                 massaction_jumps = mass_act_jump)
 ```
 
 ## How can I set the random number generator used in the jump process sampling algorithms (SSAs)?
@@ -66,16 +71,19 @@ default. On versions below 1.7 it uses `Xoroshiro128Star`.
 ## What are these aggregators and aggregations in JumpProcesses?
 
 JumpProcesses provides a variety of methods for sampling the time the next
-`ConstantRateJump` or `MassActionJump` occurs, and which jump type happens at
-that time. These methods are examples of stochastic simulation algorithms
-(SSAs), also known as Gillespie methods, Doob's method, or Kinetic Monte Carlo
-methods. In the JumpProcesses terminology we call such methods "aggregators", and
-the cache structures that hold their basic data "aggregations". See [Constant
-Rate Jump Aggregators](@ref) for a list of the available SSA aggregators.
+`ConstantRateJump`, `MassActionJump`, or `VariableRateJump` occurs, and which
+jump type happens at that time. These methods are examples of stochastic
+simulation algorithms (SSAs), also known as Gillespie methods, Doob's method, or
+Kinetic Monte Carlo methods. These are all names for jump (or point) processes
+simulation methods used across the biology, chemistry, engineering, mathematics,
+and physics literature. In the JumpProcesses terminology we call such methods
+"aggregators", and the cache structures that hold their basic data
+"aggregations". See [Jump Aggregators for Exact Simulation](@ref) for a list of
+the available SSA aggregators.
 
 ## How should jumps be ordered in dependency graphs?
 Internally, JumpProcesses SSAs (aggregators) order all `MassActionJump`s first,
-then all `ConstantRateJumps`. i.e. in the example
+then all `ConstantRateJumps` and/or `VariableRateJumps`. i.e. in the example
 
 ```julia
 using JumpProcesses
@@ -99,15 +107,15 @@ The four jumps would be ordered by the first jump in `maj`, the second jump in
 `maj`, `cj1`, and finally `cj2`. Any user-generated dependency graphs should
 then follow this ordering when assigning an integer id to each jump.
 
-See also [Constant Rate Jump Aggregators Requiring Dependency Graphs](@ref) for
+See also [Jump Aggregators Requiring Dependency Graphs](@ref) for
 more on dependency graphs needed for the various SSAs.
 
-## How do I use callbacks with `ConstantRateJump` or `MassActionJump` systems?
+## How do I use callbacks with jump simulations?
 
-Callbacks can be used with `ConstantRateJump`s and `MassActionJump`s. When
-solving a pure jump system with `SSAStepper`, only discrete callbacks can be
-used (otherwise a different time stepper is needed). When using an ODE or SDE
-time stepper any callback should work.
+Callbacks can be used with `ConstantRateJump`s, `MassActionJump`s, and
+`VariableRateJump`s. When solving a pure jump system with `SSAStepper`, only
+discrete callbacks can be used (otherwise a different time stepper is needed).
+When using an ODE or SDE time stepper any callback should work.
 
 *Note, when modifying `u` or `p` within a callback, you must call
 [`reset_aggregated_jumps!`](@ref) after making updates.* This ensures that the

docs/src/index.md

@@ -1,9 +1,10 @@
 # JumpProcesses.jl: Stochastic Simulation Algorithms for Jump Processes, Jump-ODEs, and Jump-Diffusions
 JumpProcesses.jl, formerly DiffEqJump.jl, provides methods for simulating jump
-processes, known as stochastic simulation algorithms (SSAs), Doob's method,
-Gillespie methods, or Kinetic Monte Carlo methods across different fields of
-science. It also enables the incorporation of jump processes into hybrid
-jump-ODE and jump-SDE models, including jump diffusions.
+(or point) processes. Across different fields of science such methods are also
+known as stochastic simulation algorithms (SSAs), Doob's method, Gillespie
+methods, or Kinetic Monte Carlo methods . It also enables the incorporation of
+jump processes into hybrid jump-ODE and jump-SDE models, including jump
+diffusions.
 
 JumpProcesses is a component package in the [SciML](https://sciml.ai/) ecosystem,
 and one of the core solver libraries included in
@@ -78,31 +79,33 @@ versioninfo() # hide
 ```
 ```@example
 using Pkg # hide
-Pkg.status(;mode = PKGMODE_MANIFEST) # hide
+Pkg.status(; mode = PKGMODE_MANIFEST) # hide
 ```
 ```@raw html
 </details>
 ```
 ```@raw html
-You can also download the 
+You can also download the
 <a href="
 ```
 ```@eval
 using TOML
-version = TOML.parse(read("../../Project.toml",String))["version"]
-name = TOML.parse(read("../../Project.toml",String))["name"]
-link = "https://github.com/SciML/"*name*".jl/tree/gh-pages/v"*version*"/assets/Manifest.toml"
+projtoml = joinpath("..", "..", "Project.toml")
+version = TOML.parse(read(projtoml, String))["version"]
+name = TOML.parse(read(projtoml, String))["name"]
+link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version * "/assets/Manifest.toml"
 ```
 ```@raw html
 ">manifest</a> file and the
 <a href="
 ```
 ```@eval
 using TOML
-version = TOML.parse(read("../../Project.toml",String))["version"]
-name = TOML.parse(read("../../Project.toml",String))["name"]
-link = "https://github.com/SciML/"*name*".jl/tree/gh-pages/v"*version*"/assets/Project.toml"
+projtoml = joinpath("..", "..", "Project.toml")
+version = TOML.parse(read(projtoml, String))["version"]
+name = TOML.parse(read(projtoml, String))["name"]
+link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version * "/assets/Project.toml"
 ```
 ```@raw html
 ">project</a> file.
-```
+```

docs/src/jump_solve.md

@@ -6,19 +6,37 @@ solve(prob::JumpProblem,alg;kwargs)
 
 ## Recommended Methods
 
-A `JumpProblem(prob,aggregator,jumps...)` comes in two forms. The first major
-form is if it does not have a `RegularJump`. In this case, it can be solved with
-any integrator on  `prob`. However, in the case of a pure `JumpProblem` (a
-`JumpProblem` over a  `DiscreteProblem`), there are special algorithms
-available.  The `SSAStepper()` is an efficient streamlined algorithm for running
-the  `aggregator` version of the SSA for pure `ConstantRateJump` and/or
-`MassActionJump` problems. However, it is not compatible with event handling. If
-events are necessary, then `FunctionMap` does well.
-
-If there is a `RegularJump`, then specific methods must be used. The current
-recommended method is `TauLeaping` if you need adaptivity, events, etc. If you
-just need the most barebones fixed time step leaping method, then `SimpleTauLeaping`
-can have performance benefits.
+`JumpProblem`s can be solved with two classes of methods, exact and inexact.
+Exact algorithms currently sample realizations of the jump processes in
+chronological order, executing individual jumps sequentially at randomly sampled
+times. In contrast, inexact (τ-leaping) methods are time-step based, executing
+multiple occurrences of jumps during each time-step. These methods can be much
+faster as they only simulate the total number of jumps over each leap interval,
+and thus do not need to simulate the realization of every single jump. Jumps for
+use with exact simulation methods can be defined as `ConstantRateJump`s,
+`MassActionJump`s, and/or `VariableRateJump`. Jumps for use with inexact
+τ-leaping methods should be defined as `RegularJump`s.
+
+There are special algorithms available for efficiently simulating an exact, pure
+`JumpProblem` (i.e. a `JumpProblem` over a `DiscreteProblem`).  `SSAStepper()`
+is an efficient streamlined integrator for time stepping such problems from
+individual jump to jump. This integrator is named after Stochastic Simulation
+Algorithms (SSAs), commonly used naming in chemistry and biology applications
+for the class of exact jump process simulation algorithms. In turn, we denote by
+"aggregators" the algorithms that `SSAStepper` calls to calculate the next jump
+time and to execute a jump (i.e. change the system state appropriately). All
+JumpProcesses aggregators can be used with `ConstantRateJump`s and
+`MassActionJump`s, with a subset of aggregators also working with bounded
+ `VariableRateJump`s (see [the first tutorial](@ref poisson_proc_tutorial) for
+the definition of bounded `VariableRateJump`s). Although `SSAStepper()` is
+usually faster, it only supports discrete events (`DiscreteCallback`s), for pure
+jump problems requiring continuous events (`ContinuousCallback`s) the less
+performant `FunctionMap` time-stepper can be used.
+
+If there is a `RegularJump`, then inexact τ-leaping methods must be used. The
+current recommended method is `TauLeaping` if one needs adaptivity, events, etc.
+If ones only needs the most barebones fixed time-step leaping method, then
+`SimpleTauLeaping` can have performance benefits.
 
 ## Special Methods for Pure Jump Problems
 
@@ -28,9 +46,10 @@ algorithms are optimized for pure jump problems.
 
 ### JumpProcesses.jl
 
-- `SSAStepper`: a stepping algorithm for pure `ConstantRateJump` and/or
-  `MassActionJump` `JumpProblem`s. Supports handling of `DiscreteCallback`
-  and saving controls like `saveat`.
+- `SSAStepper`: a stepping integrator for `JumpProblem`s defined over
+  `DiscreteProblem`s involving `ConstantRateJump`s, `MassActionJump`s, and/or
+  bounded `VariableRateJump`s . Supports handling of `DiscreteCallback`s and
+  saving controls like `saveat`.
 
 ## RegularJump Compatible Methods