struct ExtendedJumpArray{T3<:Number, T1, T<:AbstractArray{T3<:Number, T1}, T2} <: AbstractArray{T3<:Number, 1}

Extended state definition used within integrators when there are VariableRateJumps in a system. For detailed examples and usage information, see the

• Tutorial

Fields

• u: The current state.

• jump_u: The current rate (i.e. hazard, intensity, or propensity) values for the VariableRateJumps.

Examples

using JumpProcesses, 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 with ueja.u of size N and ueja.jump_u of size num_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 VariableRateJumps all callback, ConstantRateJump, and VariableRateJumpaffect! functions will receive integrators with integrator.u an ExtendedJumpArray.

• As such, affect! functions that wish to modify the state via vector operations should use ueja.u.u to obtain the aliased state object.

• :u and :jump_u always refer to the ExtendedJumpArray fields. Any other property access (and propertynames) is forwarded to the wrapped state ueja.u, so e.g. ueja.x works when ueja.u is an ArrayPartition.

mutable struct JumpProblem{iip, P, A, C, J<:Union{Nothing, JumpProcesses.AbstractJumpAggregator}, J1, J2, J3, J4, K} <: SciMLBase.AbstractJumpProblem{P, J<:Union{Nothing, JumpProcesses.AbstractJumpAggregator}}

Defines a collection of jump processes to associate with another problem type.

• Documentation Page
• Tutorial Page
• FAQ Page

Constructors

JumpProblems 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 ConstantRateJumps or MassActionJumps), we must first construct a DiscreteProblem

prob = DiscreteProblem(u0, tspan, p)

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 ConstantRateJumps and MassActionJumps. Examples include Direct.

• discrete_jump_aggregation: The underlying state data associated with the chosen aggregator.

• jump_callback: CallBackSet with the underlying ConstantRate and VariableRate jumps.

• constant_jumps: The ConstantRateJumps.

• variable_jumps: The VariableRateJumps.

• regular_jump: The RegularJumps.

• massaction_jump: The MassActionJumps.

• kwargs: kwargs to pass on to solve call.

Keyword Arguments

• save_positions=(true,true) when including variable rates and (false,true) for constant rates, specifies whether to save the system's state (before, after) the jump occurs. NOTE, this only controls saving for non-VariableRateJumps. VariableRateJump saving is controlled at the individual jump level via the value of this kwarg passed to the individual VariableRateJumps constructor when it is created.
• spatial_system, for spatial problems the underlying spatial structure.
• hopping_constants, for spatial problems the spatial transition rate coefficients.
• use_vrj_bounds = true, set to false to disable handling bounded VariableRateJumps with a supporting aggregator (such as Coevolve). They will then be handled via the continuous integration interface, and treated like general VariableRateJumps.
• vr_aggregator, indicates the aggregator to use for sampling variable rate jumps. Current default is VR_FRM.
• tstops, time stops to pass through to the solver. Can be an AbstractVector of times or a callable (p, tspan) -> times.
• callback, callbacks to run in every solve of the problem. They are combined with any callback passed to solve or init, which replaces them instead when merge_callbacks = false is also passed. Do not store callbacks on the wrapped problem prob; see the warning below.

Please see the tutorial page in the DifferentialEquations.jl docs for usage examples and commonly asked questions.

PureLeaping()

Request that all jumps in a JumpProblem are handled by the leaping algorithm passed to solve, instead of being converted into callback-based SSA aggregators during problem construction.

Returns

• A stateless aggregator marker for JumpProblem(prob, PureLeaping(), jumps; kwargs...).

Notes

• PureLeaping is currently intended for tau-leaping algorithms such as SimpleTauLeaping and SimpleExplicitTauLeaping.
• Spatial jump problems are not supported by the PureLeaping construction path.

Examples

using JumpProcesses, DiffEqBase

rate!(out, u, p, t) = (out[1] = 0.2 * u[1])
affect!(du, u, p, t, counts, mark) = (du[1] = -counts[1])
rj = RegularJump(rate!, affect!, 1)

prob = DiscreteProblem([10], (0.0, 1.0))
jprob = JumpProblem(prob, PureLeaping(), rj)
sol = solve(jprob, SimpleTauLeaping(); dt = 0.1)

struct SSAStepper{SaveUprev} <: SciMLBase.AbstractDEAlgorithm

Highly efficient integrator for pure jump problems that involve only ConstantRateJumps, MassActionJumps, and/or VariableRateJumps with rate bounds.

Constructor

SSAStepper(; save_uprev = false)

• save_uprev: whether the integrator keeps the state at the start of each step in integrator.uprev, so that it can be evaluated anywhere in the current step (see "Evaluating the integrator" below). This copies the state on every step, at a cost proportional to the size of the state, which can dominate the per-step cost of efficient aggregators on large or spatial systems. SSAStepper() keeps no copy.

Supported state types

The state u can be a Vector or an SVector of integers or floats; a species × sites Matrix{Int} for spatial problems, which NSM and DirectCRDirect keep and other aggregators flatten into a Vector; or, for models without MassActionJumps, an integer or float scalar (not with RSSA or RSSACR). Bounded VariableRateJumps require Coevolve, and affect! functions for SVector and scalar states must assign a new value to integrator.u. Other state types are not supported. See the supported state types for details.

Evaluating the integrator

Callbacks, such as DiffEqCallbacks' SavingCallback and FunctionCallingCallback with save times, may evaluate the integrator at times before the current one with integrator(t) or integrator(out, t). The sampled path is piecewise constant, and jumps and callbacks only change the state at the end of a step, so:

• Evaluating at the current time, integrator.t, returns the current state.
• With SSAStepper(; save_uprev = true), evaluating at any t in [integrator.tprev, integrator.t) returns the state at the start of the step, integrator.uprev. get_tmp_cache(integrator) then also provides scratch space for callbacks that evaluate the integrator in place.
• Any other time raises an ArgumentError: a path cannot be extrapolated, and earlier states are only available from the saved solution. With the default, SSAStepper(), every time other than integrator.t raises this error, as does get_tmp_cache.

integrator(t) returns a copy for array states, and integrator(out, t) writes into out and returns it. Instead of save_uprev = true, the times at which callbacks evaluate the integrator can be passed as tstops, so that it stops exactly at them. This avoids the per-step copy, but adds a stop and a pass through the callbacks at each time, and, when every jump is saved, a saved point at each stop, so pair it with save_positions = (false, false). To save the full state, saveat is exact and needs neither. DiffEqCallbacks' integrating callbacks (IntegratingCallback, IntegratingSumCallback) are not yet supported with SSAStepper.

Notes

• Only works with JumpProblems defined from DiscreteProblems.
• Only works with collections of ConstantRateJumps, MassActionJumps, and VariableRateJumps with rate bounds.
• Only supports DiscreteCallbacks for events, which are checked after every step taken by SSAStepper.
• Only supports a limited subset of the output controls from the common solver interface, specifically save_start, save_end, and saveat.
• Supports rng and seed keyword arguments in solve/init to control the random number generator used for jump sampling. rng accepts any AbstractRNG, while seed creates a Xoshiro generator. rng takes priority over seed.
• Reuses the JumpProblem's jump state (aggregator and jump callback), re-initializing it at every init; it is never copied. One JumpProblem therefore supports one active solve or integrator at a time. For concurrent solves use independent copies, for example deepcopy(jprob), or an EnsembleProblem. See the JumpProblem docstring.
• Supports the common alias keyword argument for its u0, p, and tstops inputs, following SciMLBase's alias-specifier convention: pass nothing (the default), a Bool that applies to all three, a SciMLBase.DiscreteAliasSpecifier (alias_u0, alias_p), or an ODEAliasSpecifier (which also has alias_tstops). true permits aliasing and false requests a copy. The defaults copy u0 (with recursivecopy), reuse p, and never modify a caller's tstops array. With alias_u0 = true the integrator mutates the problem's u0, so a later solve starts from its final state. alias_f and alias_du0 have no effect, and alias never controls jump state.
• tstops may be a number, a tuple, or an AbstractVector of times (such as a range or view), or a callable (p, tspan) -> times. Only a Vector of the time type can be aliased; other supported inputs are copied into a new Vector at init.
• Saved states are independent snapshots of the integrator state made with recursivecopy, as in OrdinaryDiffEq.
• As when using jumps with ODEs and SDEs, saving controls for whether to save each time a jump occurs are via the save_positions keyword argument to JumpProblem. Note that when choosing SSAStepper as the timestepper, save_positions = (true,true), (true,false), or (false,true) are all equivalent. SSAStepper will save only the post-jump state in the solution object in each of these cases. This is because solution objects generated via SSAStepper use piecewise-constant interpolation, and can therefore exactly reconstruct the sampled jump process path with knowing just the post-jump state. That is, sol(t) for any 0 <= t <= tstop will give the exact value of the sampled solution path at t provided at least one component of save_positions is true.

Examples

SIR model:

using JumpProcesses
β = 0.1 / 1000.0;
ν = 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.

David F. Anderson, Masanori Koyama; An asymptotic relationship between coupling methods for stochastically modeled population processes. IMA J Numer Anal 2015; 35 (4): 1757-1778. doi: 10.1093/imanum/dru044

reset_aggregated_jumps!(integrator, uprev = nothing)

Reset the state of jump processes and associated solvers following a change in parameters or such. Rate updates are handled automatically by initialize! via fill_scaled_rates!.

remake(jprob::JumpProblem; u0 = missing, p = missing, tspan, prob)

Construct a JumpProblem with a new initial condition u0, parameters p, time span tspan, or wrapped problem prob.

A new wrapped problem prob must not store callbacks, as for the JumpProblem constructor; an ArgumentError is raised if it does.

The new problem shares the original's jump aggregator and jump callbacks rather than copying them. Solvers owned by JumpProcesses re-initialize this state from the problem being solved at every init, so the original and remade problems can be solved one after another. They must not be solved concurrently, or used for integrators that are alive at the same time; use deepcopy to obtain an independent problem. SortingDirect's learned search order is also shared.

struct ConstantRateJump{F1, F2, B<:Union{Nothing, JumpProcesses.RateBoundFunctions}} <: JumpProcesses.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

• Tutorial

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.

• bounds: Optional RateBoundFunctions holding user supplied bracketing functions.

Examples

Suppose u[1] gives the number 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).

Rate bounds may be supplied separately, for use with bracketing aggregators such as RSSA:

rate(u, p, t) = p[1] * u[1] / (1 + u[2])
affect!(integrator) = integrator.u[1] -= 1
function bounds(ulow, uhigh, u, p, t)
    RateBounds(lrate = p[1]*ulow[1] / (1+uhigh[2]),
        urate = p[1]*uhigh[1] / (1+ulow[2]))
end
crj = ConstantRateJump(rate, affect!; bounds)

Notes

• When rate bounds are not supplied, they are computed by evaluating the rate at ulow and uhigh. These bounds are only correct if the rate is monotonic with respect to state, i.e. increasing in all species or decreasing in all species.

struct MassActionJump{T, S, U, V} <: JumpProcesses.AbstractMassActionJump

Optimized representation for ConstantRateJumps 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

• Main Docs
• Tutorial

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 stored (scaled) fixed reaction rate constants, or nothing for a parameter-mapped definition. Numeric working rates live in solver/aggregator buffers and are refreshed at initialization/reset. Stored definition arrays can be shared across remade problems; do not mutate them during concurrent solves. Prefer parameter indices or a custom mapper for rates that change.

• 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.

• rescale_rates_on_update: Whether the built-in mapper callable should apply stoichiometric scaling to rates.

Keyword Arguments

• scale_rates = true, whether to rescale the reaction rate constants according to the stoichiometry.
• nocopy = false, whether the MassActionJump can alias the scaled_rates and reactant_stoch from the input. Note, if scale_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; a scalar index is accepted for one reaction.
• param_mapper = nothing, an alternative to param_idxs: a callable (dest, maj, params) that fills the final scaled working rates in dest and returns nothing. It must not mutate the definition or input parameters.
• rescale_rates_on_update = scale_rates, the scaling policy used by the built-in parameter mapper. Custom mappers implement their own scaling and should honor this policy when mapping unscaled coefficients; pre-scaled output must not be scaled again.

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, ν = 0.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

• By default, fixed reaction rates are rescaled when constructing the MassActionJump; the built-in parameter mapper rescales rates when filling each solver's working buffer at initialization/reset, as explained in the main docs. Disable this with the kwarg scale_rates=false.
• Also see the main docs for how to specify reactions with no products or no reactants.
• For non-spatial problems the state u must be a Vector or an SVector of integers or floats; scalar and other states are not supported. Spatial problems, built by passing hopping_constants and a spatial_system to JumpProblem, use a species × sites Matrix{Int}. NSM and DirectCRDirect keep the matrix, while other aggregators flatten it into a Vector. See the supported state types.

struct VariableRateJump{R, F, R2, R3, R4, I, T, T2} <: JumpProcesses.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

• Tutorial

Note that two types of VariableRateJumps are currently supported, with different performance characteritistics.

• A general VariableRateJump or VariableRateJump will refer to one in which only rate and affect functions are specified.

  • These are the most general in what they can represent, but require the use of an ODEProblem or SDEProblem whose underlying timestepper handles their evolution in time (via the callback interface).
  • This is the least performant jump type in simulations.

• Bounded VariableRateJumps require passing the keyword arguments urate and rateinterval, corresponding to functions urate(u, p, t) and rateinterval(u, p, t), see below. These must calculate a time window over which the rate function is bounded by a constant. Note that it is ok if the rate bound would be violated within the time interval due to a change in u arising from another ConstantRateJump, MassActionJump or boundedVariableRateJump being executed, as the chosen aggregator will then handle recalculating the rate bound and interval. However, if the bound could be violated within the time interval due to a change in u arising from continuous dynamics such as a coupled ODE, SDE, or a general VariableRateJump, bounds should not be given. This ensures the jump is classified as a general VariableRateJump and properly handled. One can also optionally provide a lower bound function, lrate(u, p, t), via the lrate keyword argument. This can lead to increased performance. The validity of the lower bound should hold under the same conditions and rate interval as urate.

  • Bounded VariableRateJumps can currently be used in the Coevolve aggregator, and can therefore be efficiently simulated in pure-jump DiscreteProblems using the SSAStepper time-stepper.
  • These can be substantially more performant than general VariableRateJumps without the rate bound functions.

Reemphasizing, the additional user provided functions leveraged by bounded VariableRateJumps, urate(u, p, t), rateinterval(u, p, t), and the optional lrate(u, p, t) require that

• For s in [t, t + rateinterval(u, p, t)], we have that lrate(u, p, t) <= rate(u, p, s) <= urate(u, p, t).
• It is ok if these bounds would be violated during the time window due to another ConstantRateJump, MassActionJump or bounded VariableRateJump occurring. However, they must remain valid if u changes for any other reason (for example, due to continuous dynamics like ODEs, SDEs, or general VariableRateJumps).

Fields

• rate: Function rate(u,p,t) that returns the jump's current rate given state u, parameters p and time t.

• affect!: Function affect!(integrator) that updates the state for one occurrence of the jump given integrator.

• lrate: Optional function lrate(u, p, t) that computes a lower bound on the rate in the interval t to t + rateinterval(u, p, t) at time t given state u and parameters p. This bound must rigorously hold during the time interval as long as another ConstantRateJump, MassActionJump, or boundedVariableRateJump has not been sampled. When using aggregators that support bounded VariableRateJumps, currently only Coevolve, providing a lower-bound can lead to improved performance.

• urate: Optional function urate(u, p, t) for general VariableRateJumps, but is required to define a bounded VariableRateJump, which can be used with supporting aggregators, currently only Coevolve, and offers improved computational performance. Computes an upper bound for the rate in the interval t to t + rateinterval(u, p, t) at time t given state u and parameters p. This bound must rigorously hold during the time interval as long as another ConstantRateJump, MassActionJump, or boundedVariableRateJump has not been sampled.

• rateinterval: Optional function rateinterval(u, p, t) for general VariableRateJumps, but is required to define a bounded VariableRateJump, which can be used with supporting aggregators, currently only Coevolve, and offers improved computational performance. Computes the time interval from time t over which the urate and lrate bounds will hold, t to t + rateinterval(u, p, t), given state u and parameters p. This bound must rigorously hold during the time interval as long as another ConstantRateJump, MassActionJump, or boundedVariableRateJump has not been sampled.

• idxs

• rootfind

• interp_points

• save_positions

• abstol

• reltol

Examples

Suppose u[1] gives the number 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
vrj = VariableRateJump(rate, affect!)

To define a bounded VariableRateJump that can be used with supporting aggregators such as Coevolve, we must define bounds and a rate interval:

rateinterval(u, p, t) = (1 / p[1]) * 2
rate(u, p, t) = t * p[1] * u[1]
lrate(u, p, t) = rate(u, p, t)
urate(u, p, t) = rate(u, p, t + rateinterval(u, p, t))
affect!(integrator) = integrator.u[1] -= 1
vrj = VariableRateJump(rate, affect!; lrate = lrate, urate = urate,
    rateinterval = rateinterval)

Notes

• When using an aggregator that supports bounded VariableRateJumps, DiscreteProblem can be used. Otherwise, ODEProblem or SDEProblem must be used.
• When not using aggregators that support bounded VariableRateJumps, or when there are general VariableRateJumps, integrators store an effective state type that wraps the main state vector. See ExtendedJumpArray for details on using this object. In this case all ConstantRateJump, VariableRateJump and callback affect! functions receive an integrator with integrator.u an ExtendedJumpArray.
• 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 VariableRateJumps within ODE/SDE integrators.

struct RegularJump{iip, R, C, MD}

Representation for encoding rates and multiple simultaneous jumps for use in τ-leaping type methods.

Constructors

• RegularJump(rate, c, numjumps; mark_dist = nothing)

Fields

• rate: Function rate!(rate_vals, u, p, t) that returns the current rates, i.e. intensities or propensities, for all possible jumps in rate_vals.

• c: Function c(du, u, p, t, counts, mark) that executes the ith jump counts[i] times, saving the output in du[i].

• numjumps: Number of jumps in the system.

• mark_dist: A distribution for marks. Not currently used or supported.

Examples

function rate!(out, u, p, t)
    out[1] = (0.1 / 1000.0) * u[1] * u[2]
    out[2] = 0.01u[2]
    nothing
end

const dc = zeros(3,2)
function c(du, u, p, t, counts, mark)
    mul!(du, dc, counts)
    nothing
end

rj = RegularJump(rate!, c, 2)

## Notes
- `mark_dist` is not currently used or supported in τ-leaping methods.

struct JumpSet{T1, T2, T3, T4} <: JumpProcesses.AbstractJump

Defines a collection of jumps that should collectively be included in a simulation.

Fields

• variable_jumps: Collection of VariableRateJumps

• constant_jumps: Collection of ConstantRateJumps

• regular_jump: Collection of RegularJumps

• massaction_jump: Collection of MassActionJumps

Examples

Here we construct two jumps, store them in a JumpSet, and then simulate the resulting process.

using JumpProcesses, 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())

RateBounds(; lrate, urate, rateinterval = Inf)

Computed bounds on a jump's rate for a fixed state bracket [ulow, uhigh] at time t.

Fields

• lrate: Lower bound of the rate.

• urate: Upper bound of the rate.

• rateinterval: Time window over which the bounds hold.

Notes

• The rates must satisfy lrate <= rate(v, p, s) <= urate for every v with ulow .<= v .<= uhigh and s in [t, t + rateinterval].
• At least one of lrate or urate must be given. An omitted bound defaults to the trivially valid one: zero below and typemax above.
• For ConstantRateJumps, rateinterval is always Inf, as rates do not explicitly depend on t.

Examples

Increasing rate in one species but decreasing in another:

rate(u, p, t) = p[1] * u[1] / (1 + u[2])
function bounds(ulow, uhigh, u, p, t)
    RateBounds(lrate = p[1] * ulow[1] / (1 + uhigh[2]),
        urate = p[1] * uhigh[1] / (1 + ulow[2]))
end

Extremas lying inside the bracket:

rate(u, p, t) = p[1] * u[1] * (p[2] - u[1])
function bounds(ulow, uhigh, u, p, t)
    L = p[1] * max(abs(p[2] - 2*ulow[1]), abs(p[2] - 2*uhigh[1]))
    δr = L * max(u[1] - ulow[1], uhigh[1] - u[1])
    r = rate(u, p, t)
    RateBounds(lrate = max(r - δr, zero(r)), urate = r + δr)
end

Rate that depends on time:

rate(u, p, t) = p[1] * u[1] * exp(-p[2] * t)
function bounds(ulow, uhigh, u, p, t)
    Δ = p[3]
    RateBounds(lrate = p[1] * ulow[1] * exp(-p[2] * (t + Δ)),
        urate = p[1] * uhigh[1] * exp(-p[2] * t),
        rateinterval = Δ)
end

BracketData(fluctrate, threshold, Δu)
BracketData{T1, T2}()

Configure species-population brackets used by RSSA-based aggregators.

For species population u[i], the bracket is [(1 - fluctrate) * u[i], (1 + fluctrate) * u[i]] when u[i] >= threshold. For smaller populations, the bracket is [max(u[i] - Δu, 0), u[i] + Δu]. Each field may be either a scalar shared by all species or a vector indexed by species.

Fields

• fluctrate: Relative fluctuation width used for populations at or above threshold.
• threshold: Population threshold below which the absolute Δu bracket is used.
• Δu: Absolute bracket half-width used for populations below threshold.

Notes

• BracketData{T1, T2}() constructs BracketData(T1(0.1), T2(25), T2(4)).
• The bracketing rules follow the RSSA construction in Thanh et al., J. Chem. Phys. 142, 244106 (2015).

Examples

using JumpProcesses

bd = BracketData(0.1, 25, 4)
bd.fluctrate == 0.1

A constant-complexity NRM method. Stores next reaction times in a table with a specified bin width.

Kevin R. Sanft and Hans G. Othmer, Constant-complexity stochastic simulation algorithm with optimal binning, Journal of Chemical Physics 143, 074108 (2015). doi: 10.1063/1.4928635.

An improvement of the COEVOLVE algorithm for simulating any compound jump process that evolves through time. This method handles variable intensity rates with user-defined bounds and inter-dependent processes. It reduces to NRM when rates are constant. As opposed to COEVOLVE, this method syncs the thinning procedure with the stepper which allows it to handle dependencies on continuous dynamics.

G. A. Zagatti, S. A. Isaacson, C. Rackauckas, V. Ilin, S.-K. Ng and S. Bressan, Extending JumpProcess.jl for fast point process simulation with time-varying intensities, arXiv. doi:10.48550/arXiv.2306.06992.

M. Farajtabar, Y. Wang, M. Gomez-Rodriguez, S. Li, H. Zha, and L. Song, COEVOLVE: a joint point process model for information diffusion and network evolution, Journal of Machine Learning Research 18(1), 1305–1353 (2017). doi: 10.5555/3122009.3122050.

Gillespie's Direct method. ConstantRateJump rates and affects are stored in tuples. Fastest for a small (total) number of ConstantRateJumps or MassActionJumps (~10). For larger numbers of possible jumps, use other methods.

Daniel T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (4), 403–434 (1976). doi:10.1016/0021-9991(76)90041-3.

The 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.

The Direct Composition-Rejection Direct method. Uses the DirectCR method to determine where on the grid/graph a jump occurs, and the Direct method to determine which jump occurs at the sampled location. The state is a species × sites matrix, with u[i, j] the number of species i at site j.

Kevin R. Sanft and Hans G. Othmer, Constant-complexity stochastic simulation algorithm with optimal binning, Journal of Chemical Physics 143, 074108 (2015). doi: 10.1063/1.4928635.

Gillespie's Direct method. ConstantRateJump rates and affects are stored via FunctionWrappers, which is more performant than Direct for very large numbers of ConstantRateJumps. However, for such large numbers of jump different classes of aggregators are usually much more performant (i.e. SortingDirect, DirectCR, RSSA or RSSACR).

Daniel T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (4), 403–434 (1976). doi:10.1016/0021-9991(76)90041-3.

Gillespie's First Reaction Method. Should not be used for practical applications due to slow performance relative to all other methods.

Daniel T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (4), 403–434 (1976). doi:10.1016/0021-9991(76)90041-3.

Gillespie's First Reaction Method with FunctionWrappers for handling ConstantRateJumps. Should not be used for practical applications due to slow performance relative to all other methods.

Daniel T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (4), 403–434 (1976). doi:10.1016/0021-9991(76)90041-3.

The 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.

The Next Subvolume Method for spatial jump process simulations. Usually slower than DirectCRDirect. Uses an indexed priority queue tree structure to determine where on the grid/graph the next jump occurs, and then the Direct method to determine which jump at the given location occurs. The state is a species × sites matrix, with u[i, j] the number of species i at site j.

Elf, Johan and Ehrenberg, M, Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases,Systems Biology, 1(2), 230-236 (2004). doi:10.1049/sb:20045021.

A rejection-based direct method.

The Rejection SSA method. One of the best methods for systems with hundreds to many thousands of jumps (along with RSSACR) and sparse dependency graphs. Requires an array state, such as a Vector or SVector; scalar states are not supported.

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.

The Rejection SSA Composition-Rejection method. Often the best performer for systems with tens of thousands of jumps and sparse dependency graphs. Requires an array state, such as a Vector or SVector; scalar states are not supported.

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, 14 (3), 657-667 (2017). doi:10.1109/TCBB.2016.2530066.

The 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.

The search order learned during a solve is kept in the problem's jump state and is not reset by init, so repeated solves of the same problem (including problems created with remake) start from the order learned so far. This speeds up many consecutive simulations. Sampling remains exact, but a given seed reproduces a trajectory only from the same learned order: for exact replay, solve a deepcopy of a problem that has never been solved. Ensemble results additionally depend on the seeds, batching, thread count, and ensemble algorithm.

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.

get_num_majumps(jumps) -> Int

Return the number of mass-action jumps represented by a jump container.

Arguments

• jumps: A MassActionJump, SpatialMassActionJump, JumpSet, or nothing.

Returns

• The number of mass-action reaction channels. nothing returns 0.

Examples

using JumpProcesses

maj = MassActionJump([1.0], [[1 => 1]], [[1 => -1]])
get_num_majumps(maj) == 1
get_num_majumps(nothing) == 0

needs_depgraph(aggregator) -> Bool

Return whether an aggregator requires a reaction dependency graph when building a JumpProblem.

Arguments

• aggregator: An AbstractAggregatorAlgorithm value such as DirectCR, SortingDirect, NRM, CCNRM, RDirect, or Coevolve.

Returns

• true when the aggregator needs a dependency graph from each reaction to the reactions whose propensities must be updated after it fires.
• false for aggregators that update all rates or otherwise do not use this graph.

Examples

using JumpProcesses

needs_depgraph(Direct()) == false
needs_depgraph(NRM()) == true

needs_vartojumps_map(aggregator) -> Bool

Return whether an aggregator requires the species-to-reaction dependency map used by RSSA-style bracketing.

Arguments

• aggregator: An AbstractAggregatorAlgorithm value.

Returns

• true for RSSA and RSSACR.
• false for other built-in aggregators.

Examples

using JumpProcesses

needs_vartojumps_map(RSSA()) == true
needs_vartojumps_map(Direct()) == false

abstract type VariableRateAggregator

An abstract type for aggregators that manage the simulation of VariableRateJumps in jump processes.

Notes

• In hybrid ODE/SDE systems with general VariableRateJumps, integrator.u may be an ExtendedJumpArray for some aggregators.

struct VR_Direct <: VariableRateAggregator

A concrete VariableRateAggregator implementing a direct method-based approach for simulating VariableRateJumps. VR_Direct (Variable Rate Direct Callback) efficiently samples jump times using one continuous callback to integrate the total intensity / propensity for all VariableRateJumps, sample when the next jump occurs, and then sample which jump occurs at this time. VR_DirectFW a separate FunctionWrapper mode, which wraps things in FunctionWrappers in cases with large numbers of jumps

Examples

Simulating a birth-death process with VR_Direct (default) and VR_DirectFW:

using JumpProcesses, OrdinaryDiffEq
u0 = [1.0]           # Initial population
p = [10.0, 0.5]      # [birth rate, death rate coefficient]
tspan = (0.0, 10.0)

# Birth jump: ∅ → X
birth_rate(u, p, t) = p[1]
birth_affect!(integrator) = (integrator.u[1] += 1; nothing)
birth_jump = VariableRateJump(birth_rate, birth_affect!)

# Death jump: X → ∅
death_rate(u, p, t) = p[2] * u[1]
death_affect!(integrator) = (integrator.u[1] -= 1; nothing)
death_jump = VariableRateJump(death_rate, death_affect!)

# Problem setup
oprob = ODEProblem((du, u, p, t) -> du .= 0, u0, tspan, p)
jprob = JumpProblem(oprob, birth_jump, death_jump; vr_aggregator = VR_Direct())
sol = solve(jprob, Tsit5())

jprob = JumpProblem(oprob, birth_jump, death_jump; vr_aggregator = VR_DirectFW())
sol = solve(jprob, Tsit5())

Notes

• VR_Direct and VR_DirectFW are expected to generally be more performant than VR_FRM.

struct VR_DirectFW <: VariableRateAggregator

Function-wrapper variant of VR_Direct for simulations with many VariableRateJumps.

VR_DirectFW uses the same direct-method callback strategy as VR_Direct, but stores the rate and affect functions in FunctionWrappers-based containers.

Returns

• A stateless VariableRateAggregator value for the vr_aggregator keyword of JumpProblem.

Examples

using JumpProcesses, OrdinaryDiffEq

u0 = [1.0]
p = [10.0, 0.5]
tspan = (0.0, 10.0)

birth_rate(u, p, t) = p[1]
birth_affect!(integrator) = (integrator.u[1] += 1; nothing)
birth_jump = VariableRateJump(birth_rate, birth_affect!)

death_rate(u, p, t) = p[2] * u[1]
death_affect!(integrator) = (integrator.u[1] -= 1; nothing)
death_jump = VariableRateJump(death_rate, death_affect!)

oprob = ODEProblem((du, u, p, t) -> du .= 0, u0, tspan, p)
jprob = JumpProblem(oprob, birth_jump, death_jump; vr_aggregator = VR_DirectFW())
sol = solve(jprob, Tsit5())

struct VR_FRM <: VariableRateAggregator

A concrete VariableRateAggregator implementing a first-reaction method variant for simulating VariableRateJumps. VR_FRM (Variable Rate First Reaction Method with Ordinary Differential Equation) uses a user-selected ODE solver to handle integrating each jump's intensity / propensity. A callback is also used for each jump to determine when its integrated intensity reaches a level corresponding to a firing time, and to then execute the affect associated with the jump at that time.

Examples

Simulating a birth-death process with VR_FRM:

using JumpProcesses, OrdinaryDiffEq
u0 = [1.0]           # Initial population
p = [10.0, 0.5]      # [birth rate, death rate]
tspan = (0.0, 10.0)

# Birth jump: ∅ → X
birth_rate(u, p, t) = p[1]
birth_affect!(integrator) = (integrator.u[1] += 1; nothing)
birth_jump = VariableRateJump(birth_rate, birth_affect!)

# Death jump: X → ∅
death_rate(u, p, t) = p[2] * u[1]
death_affect!(integrator) = (integrator.u[1] -= 1; nothing)
death_jump = VariableRateJump(death_rate, death_affect!)

# Problem setup
oprob = ODEProblem((du, u, p, t) -> du .= 0, u0, tspan, p)
jprob = JumpProblem(oprob, birth_jump, death_jump; vr_aggregator = VR_FRM())
sol = solve(jprob, Tsit5())

Notes

• Specify VR_FRM in a JumpProblem via the vr_aggregator keyword argument to select its use for handling VariableRateJumps.
• While robust, it may be less performant than VR_Direct due to its integration of each individual jump's intensity, and use of one continuous callback per jump to handle detection of jump times and implementation of state changes from that jump.

EnsembleGPUKernel()
EnsembleGPUKernel(backend)

Ensemble algorithm marker for GPU execution of tau-leaping ensemble simulations.

Arguments

• backend: Optional KernelAbstractions-compatible backend. nothing requests the default backend selected by the extension.

Fields

• backend: Backend object used by the GPU extension.
• cpu_offload: Fraction of trajectories to offload to CPU execution.

Random number generation

For multiple trajectories, kernel ensembles use backend-local RNGs and reject explicit host rng and rng_func keywords. The seed solve keyword is supported with the KernelAbstractions CPU backend and seeds Julia's task-local generator. CPU-kernel replay with the same seed requires the same thread count, workload (including trajectory count), CPU backend settings, and Julia and package versions. Identical trajectories across thread counts are not guaranteed. These conditions apply to CPU kernel execution; CPU ensemble algorithms that explicitly assign RNGs to trajectories have their own RNG behavior.

Device backends reject seed with ArgumentError. Seed the device RNG through the backend's API, such as CUDA.seed!, before calling solve.

With trajectories = 1, execution delegates to EnsembleSerial() on the CPU and accepts its usual rng, seed, and ensemble rng_func keywords.

Aliasing

Kernels always build device-owned state. For multiple trajectories, the removed alias_jump keyword raises its migration error, and the alias_jumps keyword and any alias value other than nothing raise an ArgumentError, whether passed to solve or stored on the problem wrapped by the JumpProblem. Kernels never expose an integrator, so SSAStepper(; save_uprev = true) also raises an ArgumentError. With trajectories = 1, the serial solver's rules apply: SSAStepper supports alias and save_uprev, while the tau-leaping solvers reject both alias keywords.

Returns

• A SciMLBase.EnsembleAlgorithm value for use as the ensemble algorithm argument to solve.

Examples

using JumpProcesses

ensemble_alg = EnsembleGPUKernel()

struct SimpleAdaptiveTauLeaping{T<:AbstractFloat, A<:SciMLBase.AbstractDEAlgorithm} <: SciMLBase.AbstractDEAlgorithm

A tau-leaping method that switches between an explicit and an implicit step according to whether the system currently looks stiff.

Explicit tau-leaping is efficient on non-stiff systems but is limited by the fastest reaction, while an implicit step lifts that limit at the cost of a nonlinear solve per step. This solver measures stiffness at each step and pays for the implicit step only when it is needed, following Cao et al. (2007).

Stiffness is judged from the spread of the propensities, or, when eigenvalue_check is set, from the eigenvalue ratio of the Jacobian of the drift. When a step is taken implicitly the tau-selection tolerance is relaxed by implicit_epsilon_factor, since the implicit step is not restricted by the fast timescale it damps.

The implicit step itself is whichever algorithm is passed as implicit_alg, either SimpleImplicitTauLeaping or SimpleTrapezoidalLeaping.

Fields

• epsilon: Error control parameter used when selecting tau.

• implicit_alg: The algorithm used for steps that are taken implicitly.

• eigenvalue_check: Whether to judge stiffness from the eigenvalues of the drift Jacobian.

• stiffness_ratio_threshold: Eigenvalue ratio above which the system counts as stiff.

• implicit_epsilon_factor: Factor relaxing epsilon when a step is taken implicitly.

Notes

• Only works with JumpProblems defined from DiscreteProblems that contain only a MassActionJump, built with the PureLeaping() aggregator.
• Supports saveat, save_start and save_end.

Examples

using JumpProcesses

maj = MassActionJump([1.0, 1.0], [[1 => 1], [2 => 1]], [
    [1 => -1, 2 => 1], [1 => 1, 2 => -1]])
prob = DiscreteProblem([100, 100], (0.0, 10.0))
jprob = JumpProblem(prob, PureLeaping(), maj)

sol = solve(jprob, SimpleAdaptiveTauLeaping())
sol = solve(jprob, SimpleAdaptiveTauLeaping(implicit_alg = SimpleTrapezoidalLeaping()))

SimpleExplicitTauLeaping(; epsilon = 0.05)
SimpleExplicitTauLeaping(epsilon)

Adaptive explicit tau-leaping algorithm for pure MassActionJump problems.

Use SimpleExplicitTauLeaping with JumpProblem(prob, PureLeaping(), mass_action_jump). The algorithm computes step sizes from the mass-action propensities and the error-control parameter epsilon.

Arguments

• epsilon: Positive floating-point error-control parameter. Smaller values generally produce smaller steps.

Fields

• epsilon: Stored error-control parameter used by the adaptive leaping step selector.

Returns

• A SciMLBase.AbstractDEAlgorithm value.

Examples

using JumpProcesses, DiffEqBase

maj = MassActionJump([0.1], [[1 => 1]], [[1 => -1]])
prob = DiscreteProblem([20], (0.0, 2.0))
jprob = JumpProblem(prob, PureLeaping(), maj)
sol = solve(jprob, SimpleExplicitTauLeaping())

SimpleImplicitTauLeaping(; epsilon = 0.05)

An implicit tau-leaping method for stiff pure-jump problems.

Explicit tau-leaping is limited by the fastest reaction in the system, so a stiff model forces a step size far smaller than the timescale of interest. Each step here instead solves a nonlinear equation for the new state, which lifts that restriction; see Rathinam et al. (2003) and Cao et al. (2004).

The deterministic part of the step is taken implicitly and the fluctuations are then sampled with Poisson random variables, after which the step is rejected and tau halved if it would drive a population negative.

\[X(t + \tau) = X(t) + \sum_j \nu_j a_j(X(t + \tau)) \tau\]

as in Rathinam et al. (2003) and Cao et al. (2004).

Fields

• epsilon: Error control parameter used when selecting tau.

Notes

• Only works with JumpProblems defined from DiscreteProblems that contain only a MassActionJump, built with the PureLeaping() aggregator.
• Supports saveat, save_start and save_end.

Examples

using JumpProcesses

maj = MassActionJump([1.0, 1.0], [[1 => 1], [2 => 1]], [
    [1 => -1, 2 => 1], [1 => 1, 2 => -1]])
prob = DiscreteProblem([100, 100], (0.0, 10.0))
jprob = JumpProblem(prob, PureLeaping(), maj)
sol = solve(jprob, SimpleImplicitTauLeaping())

SimpleTauLeaping()

Fixed-step tau-leaping algorithm for pure RegularJump problems.

Use SimpleTauLeaping with JumpProblem(prob, PureLeaping(), regular_jump) and pass the timestep through the dt keyword to solve.

Keyword Arguments

The algorithm constructor has no fields or keyword arguments. The solve method accepts:

• dt: Required fixed timestep.
• rng: Optional RNG used by this solve; takes precedence over seed.
• seed: Optional seed used to construct a solve-local RNG.
• saveat: Optional scalar interval or collection of save times.
• save_start: Whether to save the initial time. Defaults follow SciML save conventions.
• save_end: Whether to save the final time. Defaults follow SciML save conventions.

Returns

• A stateless SciMLBase.AbstractDEAlgorithm value.

Examples

using JumpProcesses, DiffEqBase

rate!(out, u, p, t) = (out[1] = 0.1 * u[1])
affect!(du, u, p, t, counts, mark) = (du[1] = -counts[1])
rj = RegularJump(rate!, affect!, 1)

prob = DiscreteProblem([20], (0.0, 2.0))
jprob = JumpProblem(prob, PureLeaping(), rj)
sol = solve(jprob, SimpleTauLeaping(); dt = 0.1)

SimpleTrapezoidalLeaping(; epsilon = 0.05)

An implicit trapezoidal tau-leaping method for stiff pure-jump problems.

The method averages the propensities at the current and new states,

\[X(t + \tau) = X(t) + \sum_j \nu_j \frac{a_j(X(t)) + a_j(X(t + \tau))}{2} \tau.\]

This formulation damps the excessive stiffness of the fully implicit step and keeps the equilibrium distribution closer to the exact one.

Fields

• epsilon: Error control parameter used when selecting tau.

Notes

• Only works with JumpProblems defined from DiscreteProblems that contain only a MassActionJump, built with the PureLeaping() aggregator.
• Supports saveat, save_start and save_end.

Examples

using JumpProcesses

maj = MassActionJump([1.0, 1.0], [[1 => 1], [2 => 1]], [
    [1 => -1, 2 => 1], [1 => 1, 2 => -1]])
prob = DiscreteProblem([100, 100], (0.0, 10.0))
jprob = JumpProblem(prob, PureLeaping(), maj)
sol = solve(jprob, SimpleTrapezoidalLeaping())

CartesianGrid(dims)

Construct the default Cartesian grid topology for spatial jump simulations.

CartesianGrid currently returns a CartesianGridRej, which samples random neighbors by rejection from the potential coordinate offsets.

Arguments

• dims: Tuple or vector of side lengths. One-, two-, and three-dimensional grids are supported by the built-in offset tables.

Returns

• A CartesianGridRej with sites indexed by Julia linear indices.

Examples

using JumpProcesses

grid = CartesianGrid((3, 4))
num_sites(grid) == 12
collect(neighbors(grid, 1)) == [2, 4]

CartesianGridRej(dims)
CartesianGridRej(dimension, linear_size::Int)

Cartesian grid topology with rejection-based random neighbor sampling.

Sites are represented by linear indices over dims. Neighbor relations use the nearest Cartesian offsets in one, two, or three dimensions and reject offsets that would leave the domain.

Arguments

• dims: Tuple or vector of side lengths.
• dimension: Number of Cartesian dimensions for a hypercube grid.
• linear_size: Side length used in every dimension when constructing from (dimension, linear_size).

Fields

• dims: Side lengths of the grid.
• nums_neighbors: Number of valid neighbors for each site.
• CI: Cartesian indices for the grid domain.
• LI: Linear indices for the grid domain.
• offsets: Candidate Cartesian offsets used to enumerate or sample neighbors.

Examples

using JumpProcesses

grid = CartesianGridRej((2, 2))
outdegree(grid, 1) == 2
collect(neighbors(grid, 1)) == [2, 3]

SpatialMassActionJump(uniform_rates, spatial_rates, reactant_stoch, net_stoch,
    param_mapper = nothing; scale_rates = true, useiszero = true, nocopy = false)
SpatialMassActionJump(spatial_rates, reactant_stoch, net_stoch, param_mapper = nothing;
    kwargs...)
SpatialMassActionJump(uniform_rates, reactant_stoch, net_stoch, param_mapper = nothing;
    kwargs...)
SpatialMassActionJump(ma_jumps::MassActionJump; scale_rates = false, kwargs...)

Represent mass-action reactions whose rate constants may be uniform across sites, vary by site, or include both uniform and spatially varying reaction channels.

Uniform reactions are ordered before spatially varying reactions. For a spatial rate matrix, rows index reactions and columns index sites.

Arguments

• uniform_rates: Vector of rate constants for reactions that use the same rate at every site, or nothing.
• spatial_rates: Matrix of rate constants for site-dependent reactions, or nothing.
• reactant_stoch: Reactant stoichiometry for each reaction, using the same pair-vector representation as MassActionJump.
• net_stoch: Net stoichiometry for each reaction.
• param_mapper: Optional function mapping problem parameters to rate constants.
• ma_jumps: Existing MassActionJump to reinterpret as spatial mass-action reactions.

Keyword Arguments

• scale_rates: Whether to divide rates by factorial stoichiometry factors. Defaults to true for raw rates and false when constructing from an existing MassActionJump.
• useiszero: Whether a single zero reactant entry is treated as an empty reactant list.
• nocopy: Whether to store input arrays directly instead of copying them.

Fields

• uniform_rates: Uniform reaction rates, or nothing.
• spatial_rates: Site-dependent reaction-rate matrix, or nothing.
• reactant_stoch: Reactant stoichiometry by reaction.
• net_stoch: Net stoichiometry by reaction.
• param_mapper: Optional parameter-to-rate mapping.

Examples

using JumpProcesses

rates = [0.1 0.2 0.3]
reactant_stoch = [[1 => 1]]
net_stoch = [[1 => -1]]
smaj = SpatialMassActionJump(rates, reactant_stoch, net_stoch)
get_num_majumps(smaj) == 1

neighbors(grid, site)

return an iterator over neighbors of site in ascending order. Do not use in hot loops

num_sites(spatial_system) -> Int

Return the number of sites in a spatial system.

Arguments

• spatial_system: A Graphs.AbstractGraph or CartesianGridRej-compatible grid.

Returns

• The number of graph vertices or Cartesian grid cells.

Examples

using JumpProcesses

grid = CartesianGrid((2, 3))
num_sites(grid) == 6

outdegree(grid, site) -> Int

Return the number of valid nearest-neighbor sites adjacent to site.

Arguments

• grid: A CartesianGridRej-compatible grid.
• site: Linear site index.

Returns

• The number of in-domain neighbors for site.

Examples

using JumpProcesses

grid = CartesianGrid((2, 2))
outdegree(grid, 1) == 2

mutable struct SSAIntegrator{F, uType, tType, tdirType, P, S, CB, SA, OPT, TS, R, UP, TC} <: JumpProcesses.AbstractSSAIntegrator{SSAStepper, Nothing, uType, tType}

Integrator 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 time at the start of the current step. A step ends at the next jump or tstop, and solve!'s final advance to the end time also counts as a step.

• uprev: The state at the start of the current step, at time tprev, when solving with SSAStepper(; save_uprev = true); otherwise nothing.

• 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 every time 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

• alias_tstops: If true, will write tstops into the user-passed array

• copied_tstops: If true indicates we have already allocated the tstops array

• rng: The random number generator.

• tmp_cache: Scratch space returned by get_tmp_cache for save_uprev = true and array states other than SVectors; otherwise nothing.