JumpProcesses.jl

When you call using JumpProcesses, the library automatically re-exports several essential types and functions from SciMLBase and DiffEq ecosystems. This allows you to build, solve, and inspect JumpProblems without importing additional packages.

Key Re-exported Components:

Problem Definitions

• DiscreteProblem, EnsembleProblem (Primary interfaces for jump processes)
• ODEProblem, SDEProblem (For hybrid or paired problems)
• remake, NullParameters

Solving and Execution

• solve, solve!, init, step!
• EnsembleSerial, EnsembleThreads, EnsembleDistributed, EnsembleSplitThreads (For parallel/distributed ensemble solving)

Integrator Interface

• add_tstop!, add_saveat!, savevalues!, set_proposed_dt!, set_t!, set_u!, reinit!, terminate!, u_modified!, derivative_discontinuity!

Callbacks

• DiscreteCallback (Supported by SSAStepper)
• ContinuousCallback, VectorContinuousCallback (Used with ODE/SDE integrators)
• CallbackSet

Solutions and Status

• ODESolution, EnsembleSolution, EnsembleSummary, EnsembleAnalysis
• ReturnCode, successful_retcode

Note: If you need BVP, DAE, DDE, or optimization-specific classes, you must import them directly from SciMLBase as they are not re-exported.

JumpProcesses.jl models stochastic processes where jumps occur at random times. The core components are:

• Rate Function (rate(u, p, t)): A user-provided function that returns the intensity $\lambda_i$ (the probability per unit time) of a jump occurring, given the current state $u$, parameters $p$, and time $t$.
• Affect Function (affect!(u, p, t)): A user-provided function that defines how the state $u$ changes when a jump occurs. This allows for general state updates beyond simple increments, such as random changes or marking.
• Aggregators: Algorithms used to determine the time of the next jump and which specific jump $N_i$ occurs.
• Time-stepping methods: Methods used to advance the state from one jump to the next.

Mathematically, the library simulates processes of the form:

$$du = \sum_{i=1}^I h_i(u,p,t) , dN_i(t)$$

where $h_i$ is the change in state. It also supports coupling jumps with continuous dynamics (ODEs or SDEs) to create jump-diffusion or piecewise deterministic Markov processes.

To bridge the gap between jump process simulation and temporal point process (TPP) theory, you can implement the SciMLPointProcess structure. This structure acts as a concrete implementation of PointProcesses.AbstractPointProcess that leverages solvers from JumpProcesses and SciML.

A SciMLPointProcess is modeled as a multivariate TPP where each subcomponent is a marked TPP on a continuous space. It requires:

• jumps: A vector of jump processes (Vector{J <: JumpProcesses.AbstractJump}).
• mark_dist: A vector of mark distributions (Vector{D}).
• g: A sparsely connected graph representing the relationships between processes.
• p: Parameters.
• tmin, tmax: The temporal boundaries of the process.

using JumpProcesses
using PointProcesses

struct SciMLPointProcess{M, J <: JumpProcesses.AbstractJump, G, D, T <: Real} <:
       AbstractPointProcess{M}
    jumps::Vector{J}
    mark_dist::Vector{D}
    g::G
    p::Any
    tmin::T
    tmax::T
end

When exactness of stepping to every jump time is not required, you can use RegularJump to implement τ-leaping. This method pools jumps together to perform larger, more efficient updates, introducing a controlled time-discretization error.

To define a RegularJump, you must provide:

1. A rate function: Computes the rates of each jump process as a vector.
2. A change function: Calculates the change in state (du) given the number of times each jump occurred during a time-step (counts).
3. The number of jumps being encoded.

Note: When a JumpProblem contains a RegularJump, you must use a τ-leaping algorithm (such as SimpleTauLeaping) for simulation.

# 1. Define the rate function
function rate(out, u, p, t)
    out[1] = p[1] * u[1] * u[2]   # β * S * I
    out[2] = p[2] * u[2]          # ν * I
    nothing
end

# 2. Define the change function
# 'counts' is a vector of how many times each jump occurred
function change(du, u, p, t, counts, mark)
    # Example using a stoichiometry matrix 'c'
    mul!(du, c, counts)
    nothing
end

# 3. Create the RegularJump
rj = RegularJump(rate, change, 2)

# 4. Build the JumpProblem with a leaping aggregator
u₀ = [1000.0, 50.0, 0.0]
prob = DiscreteProblem(u₀, tspan, p)
jump_prob = JumpProblem(prob, PureLeaping(), rj)

# 5. Solve using a τ-leaping algorithm like SimpleTauLeaping
sol = solve(jump_prob, SimpleTauLeaping(); dt = 0.001)

When transition rates depend on time or state variables in a way that is not constant between jumps, use the VariableRateJump type. There are two ways to implement this:

1. General VariableRateJump: Used when you do not have known bounds on the rate or the time interval. This requires coupling the jump process to a continuous problem (like an ODEProblem) and using an ODE solver (e.g., Tsit5()) as the stepper. All jumps in the system must be redefined as general VariableRateJumps if they depend on variables modified by a general jump.

2. Bounded VariableRateJump: More performant but more restrictive. It requires providing functions that define an upper bound (urate), an optional lower bound (lrate), and a time interval (rateinterval) over which the rate is guaranteed to stay below the upper bound. This can be used with the Coevolve() aggregator and SSAStepper().

Requirements for Bounded Jumps:

• rate1(u, p, t) <= urate(u, p, t) must hold for all t <= s <= t + rateinterval(u, p, t).
• You must specify a dependency graph (dep_graph) to indicate which jumps need their rates or bounds recalculated when a specific jump occurs.

# Bounded VariableRateJump example
rate1(u, p, t) = p.λ * (sin(pi * t / 2) + 1)
affect1!(integrator) = (integrator.u[1] += 1)

rateinterval(u, p, t) = typemax(t)
urate(u, p, t) = 2 * p.λ
lrate(u, p, t) = p.λ

vrj1 = VariableRateJump(rate1, affect1!; lrate, urate, rateinterval)

# Dependency graph: dep_graph[i] contains indices of jumps that depend on jump i
# Example: jump 1 (death) affects jump 2 (birth), and jump 2 affects jump 1
dep_graph = [[1], [1, 2]]

jprob = JumpProblem(dprob, Coevolve(), vrj1, deathcrj; dep_graph)
sol = solve(jprob, SSAStepper())

JumpProcesses simulates TPPs by treating the count of events $N(t)$ as the state of a dynamical system.

• ConstantRateJump: Used for processes with a constant rate between points (e.g., Homogeneous Poisson Process).
• VariableRateJump: Used for processes where the rate $\lambda(t)$ can change arbitrarily based on time or the current state $u$.
• DiscreteProblem: The base problem type used when the system evolves in discrete steps (typical for TPPs).
• JumpProblem: Wraps a DiscreteProblem with a simulation algorithm (aggregator) and one or more jump processes.
• Aggregators:
  • Direct(): A thinning algorithm for constant rates.
  • Coevolve(): A thinning algorithm for multivariate processes that uses a dependency graph to handle interacting rates.
• SSAStepper: A discrete time-stepper that advances the system specifically to the times proposed by the simulation algorithm.