DiffEqJump.jl API
Core Types
DiffEqJump.JumpProblem
— Typemutable struct JumpProblem{iip, P, A, C, J<:Union{Nothing, DiffEqJump.AbstractJumpAggregator}, J2, J3, J4, R} <: SciMLBase.AbstractJumpProblem{P, J<:Union{Nothing, DiffEqJump.AbstractJumpAggregator}}
Defines a collection of jump processes to associate with another problem type.
Constructors
JumpProblem
s can be constructed by first building another problem type to which the jumps will be associated. For example, to simulate a collection of jump processes for which the transition rates are constant between jumps (called ConstantRateJump
s or MassActionJump
s), we must first construct a DiscreteProblem
prob = DiscreteProblem(u0, p, tspan)
where u0
is the initial condition, p
the parameters and tspan
the time span. If we wanted to have the jumps coupled with a system of ODEs, or have transition rates with explicit time dependence, we would use an ODEProblem
instead that defines the ODE portion of the dynamics.
Given prob
we define the jumps via
JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps::JumpSet ; kwargs...)
JumpProblem(prob, aggregator::AbstractAggregatorAlgorithm, jumps...; kwargs...)
Here aggregator
specifies the underlying algorithm for calculating next jump times and types, for example Direct
. The collection of different AbstractJump
types can then be passed within a single JumpSet
or as subsequent sequential arguments.
Fields
prob
The type of problem to couple the jumps to. For a pure jump process use
DiscreteProblem
, to couple to ODEs,ODEProblem
, etc.aggregator
The aggregator algorithm that determines the next jump times and types for
ConstantRateJump
s andMassActionJump
s. Examples includeDirect
.discrete_jump_aggregation
The underlying state data associated with the chosen aggregator.
jump_callback
CallBackSet
with the underlyingConstantRate
andVariableRate
jumps.variable_jumps
The
VariableRateJump
s.regular_jump
The
RegularJump
s.massaction_jump
The
MassActionJump
s.rng
The random number generator to use.
Keyword Arguments
rng
, the random number generator to use. On 1.7 and up defaults to Julia's builtin generator, below 1.7 uses RandomNumbers.jl'sXorshifts.Xoroshiro128Star(rand(UInt64))
.save_positions=(true,true)
, specifies whether to save the system's state (before,after) the jump occurs.spatial_system
, for spatial problems the underlying spatial structure.hopping_constants
, for spatial problems the spatial transition rate coefficients.
Please see the tutorial page in the DifferentialEquations.jl docs for usage examples and commonly asked questions.
DiffEqJump.SSAStepper
— Typestruct SSAStepper <: SciMLBase.DEAlgorithm
Highly efficient integrator for pure jump problems that involve only ConstantRateJump
s and/or MassActionJump
s.
Notes
- Only works with
JumProblem
s defined fromDiscreteProblem
s. - Only works with collections of
ConstantRateJump
s andMassActionJump
s. - Only supports
DiscreteCallback
s for events.
Examples
SIR model:
using DiffEqJump
β = 0.1 / 1000.0; ν = .01;
p = (β,ν)
rate1(u,p,t) = p[1]*u[1]*u[2] # β*S*I
function affect1!(integrator)
integrator.u[1] -= 1 # S -> S - 1
integrator.u[2] += 1 # I -> I + 1
end
jump = ConstantRateJump(rate1,affect1!)
rate2(u,p,t) = p[2]*u[2] # ν*I
function affect2!(integrator)
integrator.u[2] -= 1 # I -> I - 1
integrator.u[3] += 1 # R -> R + 1
end
jump2 = ConstantRateJump(rate2,affect2!)
u₀ = [999,1,0]
tspan = (0.0,250.0)
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, Direct(), jump, jump2)
sol = solve(jump_prob, SSAStepper())
see the tutorial for details.
DiffEqJump.reset_aggregated_jumps!
— Functionreset_aggregated_jumps!(integrator, uprev = nothing; update_jump_params=true)
Reset the state of jump processes and associated solvers following a change in parameters or such.
Notes - update_jump_params=true
will recalculate the rates stored within any MassActionJump that was built from the parameter vector. If the parameter vector is unchanged this can safely be set to false to improve performance.
Types of Jumps
DiffEqJump.ConstantRateJump
— Typestruct ConstantRateJump{F1, F2} <: DiffEqJump.AbstractJump
Defines a jump process with a rate (i.e. hazard, intensity, or propensity) that does not explicitly depend on time. More precisely, one where the rate function is constant between the occurrence of jumps. For detailed examples and usage information see the
Fields
rate
Function
rate(u,p,t)
that returns the jump's current rate.affect!
Function
affect(integrator)
that updates the state for one occurrence of the jump.
Examples
Suppose u[1]
gives the amount of particles and p[1]
the probability per time each particle can decay away. A corresponding ConstantRateJump
for this jump process is
rate(u,p,t) = p[1]*u[1]
affect!(integrator) = integrator.u[1] -= 1
crj = ConstantRateJump(rate, affect!)
Notice, here that rate
changes in time, but is constant between the occurrence of jumps (when u[1]
will decrease).
DiffEqJump.MassActionJump
— Typestruct MassActionJump{T, S, U, V} <: DiffEqJump.AbstractMassActionJump
Optimized representation for ConstantRateJump
s that can be represented in mass action form, offering improved performance within jump algorithms compared to ConstantRateJump
. For detailed examples and usage information see the
Constructors
MassActionJump(reactant_stoich, net_stoich; scale_rates=true, param_idxs=nothing)
Here reactant_stoich
denotes the reactant stoichiometry for each reaction and net_stoich
the net stoichiometry for each reaction.
Fields
scaled_rates
The (scaled) reaction rate constants.
reactant_stoch
The reactant stoichiometry vectors.
net_stoch
The net stoichiometry vectors.
param_mapper
Parameter mapping functor to identify reaction rate constants with parameters in
p
vectors.
Keyword Arguments
scale_rates=true
, whether to rescale the reaction rate constants according to the stoichiometry.nocopy=false
, whether theMassActionJump
can alias thescaled_rates
andreactant_stoch
from the input. Note, ifscale_rates=true
this will potentially modify both of these.param_idxs=nothing
, indexes in the parameter vector,JumpProblem.prob.p
, that correspond to each reaction's rate.
See the tutorial and main docs for details.
Examples
An SIR model with S + I --> 2I
at rate β as the first reaction and I --> R
at rate ν as the second reaction can be encoded by
p = (β=1e-4, ν=.01)
u0 = [999, 1, 0] # (S,I,R)
tspan = (0.0, 250.0)
rateidxs = [1, 2] # i.e. [β,ν]
reactant_stoich =
[
[1 => 1, 2 => 1], # 1*S and 1*I
[2 => 1] # 1*I
]
net_stoich =
[
[1 => -1, 2 => 1], # -1*S and 1*I
[2 => -1, 3 => 1] # -1*I and 1*R
]
maj = MassActionJump(reactant_stoich, net_stoich; param_idxs=rateidxs)
prob = DiscreteProblem(u0, tspan, p)
jprob = JumpProblem(prob, Direct(), maj)
Notes
DiffEqJump.VariableRateJump
— Typestruct VariableRateJump{R, F, I, T, T2} <: DiffEqJump.AbstractJump
Defines a jump process with a rate (i.e. hazard, intensity, or propensity) that may explicitly depend on time. More precisely, one where the rate function is allowed to change between the occurrence of jumps. For detailed examples and usage information see the
Fields
rate
Function
rate(u,p,t)
that returns the jump's current rate.affect!
Function
affect(integrator)
that updates the state for one occurrence of the jump.idxs
rootfind
interp_points
save_positions
abstol
reltol
Examples
Suppose u[1]
gives the amount of particles and t*p[1]
the probability per time each particle can decay away. A corresponding VariableRateJump
for this jump process is
rate(u,p,t) = t*p[1]*u[1]
affect!(integrator) = integrator.u[1] -= 1
crj = VariableRateJump(rate, affect!)
Notes
VariableRateJump
s result inintegrator
s storing an effective state type that wraps the main state vector. SeeExtendedJumpArray
for details on using this object. Note that the presence of anyVariableRateJump
s will result in allConstantRateJump
,VariableRateJump
and callbackaffect!
functions receiving an integrator withintegrator.u
anExtendedJumpArray
.- Must be used with
ODEProblem
s orSDEProblem
s to be correctly simulated (i.e. can not currently be used withDiscreteProblem
s). - Salis H., Kaznessis Y., Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions, Journal of Chemical Physics, 122 (5), DOI:10.1063/1.1835951 is used for calculating jump times with
VariableRateJump
s within ODE/SDE integrators.
DiffEqJump.JumpSet
— Typestruct JumpSet{T1, T2, T3, T4} <: DiffEqJump.AbstractJump
Defines a collection of jumps that should collectively be included in a simulation.
Fields
variable_jumps
Collection of
VariableRateJump
sconstant_jumps
Collection of
ConstantRateJump
sregular_jump
Collection of
RegularJump
smassaction_jump
Collection of
MassActionJump
s
Examples
Here we construct two jumps, store them in a JumpSet
, and then simulate the resulting process.
using DiffEqJump, OrdinaryDiffEq
rate1(u,p,t) = p[1]
affect1!(integrator) = (integrator.u[1] += 1)
crj = ConstantRateJump(rate1, affect1!)
rate2(u,p,t) = (t/(1+t))*p[2]*u[1]
affect2!(integrator) = (integrator.u[1] -= 1)
vrj = VariableRateJump(rate2, affect2!)
jset = JumpSet(crj, vrj)
f!(du,u,p,t) = (du .= 0)
u0 = [0.0]
p = (20.0, 2.0)
tspan = (0.0, 200.0)
oprob = ODEProblem(f!, u0, tspan, p)
jprob = JumpProblem(oprob, Direct(), jset)
sol = solve(jprob, Tsit5())
Aggregators
Aggregators are the underlying algorithms used for sampling MassActionJump
s and ConstantRateJump
s.
DiffEqJump.Direct
— TypeGillespie's Direct method. ConstantRateJump
rates and affects are stored in tuples. Fastest for a small (total) number of ConstantRateJump
s or MassActionJump
s (~10). For larger numbers of possible jumps use other methods.
Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. Journal of Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
DiffEqJump.DirectCR
— TypeThe Composition-Rejection Direct method. Performs best relative to other methods for systems with large numbers of jumps with special structure (for example a linear chain of reactions, or jumps corresponding to particles hopping on a grid or graph).
A. Slepoy, A.P. Thompson and S.J. Plimpton, A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks, Journal of Chemical Physics, 128 (20), 205101 (2008). doi:10.1063/1.2919546
S. Mauch and M. Stalzer, Efficient formulations for exact stochastic simulation of chemical systems, ACM Transactions on Computational Biology and Bioinformatics, 8 (1), 27-35 (2010). doi:10.1109/TCBB.2009.47
DiffEqJump.FRM
— TypeGillespie's First Reaction Method. Should not be used for practical applications due to slow performance relative to all other methods.
Gillespie, Daniel T. (1976). A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. Journal of Computational Physics. 22 (4): 403–434. doi:10.1016/0021-9991(76)90041-3.
DiffEqJump.NRM
— TypeThe Next Reaction Method. Can significantly outperform Direct for systems with large numbers of jumps and sparse dependency graphs, but is usually slower than one of DirectCR
, RSSA
, or RSSACR
for such systems.
M. A. Gibson and J. Bruck, Efficient exact stochastic simulation of chemical systems with many species and many channels, Journal of Physical Chemistry A, 104 (9), 1876-1889 (2000). doi:10.1021/jp993732q
DiffEqJump.RDirect
— TypeA rejection-based direct method.
DiffEqJump.RSSA
— TypeThe Rejection SSA method. One of the best methods for systems with hundreds to many thousands of jumps (along with RSSACR
) and sparse dependency graphs.
V. H. Thanh, C. Priami and R. Zunino, Efficient rejection-based simulation of biochemical reactions with stochastic noise and delays, Journal of Chemical Physics, 141 (13), 134116 (2014). doi:10.1063/1.4896985
V. H. Thanh, R. Zunino and C. Priami, On the rejection-based algorithm for simulation and analysis of large-scale reaction networks, Journal of Chemical Physics, 142 (24), 244106 (2015). doi:10.1063/1.4922923
DiffEqJump.RSSACR
— TypeThe Rejection SSA Composition-Rejection method. Often the best performer for systems with tens of thousands of jumps and sparse depedency graphs.
V. H. Thanh, R. Zunino, and C. Priami, Efficient Constant-Time Complexity Algorithm for Stochastic Simulation of Large Reaction Networks, IEEE/ACM Transactions on Computational Biology and Bioinformatics, Vol. 14, No. 3, 657-667 (2017).
DiffEqJump.SortingDirect
— TypeThe Sorting Direct method. Often the fastest algorithm for smaller to moderate sized systems (tens of jumps), or systems where a few jumps occur much more frequently than others.
J. M. McCollum, G. D. Peterson, C. D. Cox, M. L. Simpson and N. F. Samatova, The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior, Computational Biology and Chemistry, 30 (1), 39049 (2006). doi:10.1016/j.compbiolchem.2005.10.007
Private API Functions
DiffEqJump.ExtendedJumpArray
— Typestruct ExtendedJumpArray{T3<:Number, T1, T<:AbstractArray{T3<:Number, T1}, T2} <: AbstractArray{T3<:Number, 1}
Extended state definition used within integrators when there are VariableRateJump
s in a system. For detailed examples and usage information see the
Fields
u
The current state.
jump_u
The current rate (i.e. hazard, intensity, or propensity) values for the
VariableRateJump
s.
Examples
using DiffEqJump, OrdinaryDiffEq
f(du,u,p,t) = du .= 0
rate(u,p,t) = (1+t)*u[1]*u[2]
# suppose we wish to decrease each of the two variables by one
# when a jump occurs
function affect!(integrator)
# Method 1, direct indexing works like normal
integrator.u[1] -= 1
integrator.u[2] -= 1
# Method 2, if we want to broadcast or use array operations we need
# to access integrator.u.u which is the actual state object.
# So equivalently to above we could have said:
# integrator.u.u .-= 1
end
u0 = [10.0, 10.0]
vrj = VariableRateJump(rate, affect!)
oprob = ODEProblem(f, u0, (0.0,2.0))
jprob = JumpProblem(oprob, Direct(), vrj)
sol = solve(jprob,Tsit5())
Notes
If
ueja isa ExtendedJumpArray
withueja.u
of sizeN
andueja.jump_u
of sizenum_variableratejumps
then# for 1 <= i <= N ueja[i] == ueja.u[i] # for N < i <= (N+num_variableratejumps) ueja[i] == ueja.jump_u[i]
In a system with
VariableRateJump
s all callback,ConstantRateJump
, andVariableRateJump
affect!
functions will receive integrators withintegrator.u
anExtendedJumpArray
.As such,
affect!
functions that wish to modify the state via vector operations should useueja.u.u
to obtain the aliased state object.
DiffEqJump.SSAIntegrator
— Typemutable struct SSAIntegrator{F, uType, tType, tdirType, P, S, CB, SA, OPT, TS} <: SciMLBase.DEIntegrator{SSAStepper, Nothing, uType, tType}
Solution objects for pure jump problems solved via SSAStepper
.
Fields
f
The underlying
prob.f
function. Not currently used.u
The current solution values.
t
The current solution time.
tprev
The previous time a jump occured.
tdir
The direction time is changing in (must be positive indicating time is increasing)
p
The current parameters.
sol
The current solution object.
i
tstop
The next jump time.
cb
The jump aggregator callback.
saveat
Times to save the solution at.
save_everystep
Whether to save everytime a jump occurs.
save_end
Whether to save at the final step.
cur_saveat
Index of the next
saveat
time.opts
Tuple storing callbacks.
tstops
User supplied times to step to, useful with callbacks.
tstops_idx
u_modified
keep_stepping