### JumpProcesses: Basic JumpProblem setup and solve

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/HISTORY.md

Sets up and solves a basic JumpProblem using ConstantRateJump and SSAStepper. This example illustrates the fundamental components required for defining and simulating jump processes.

```julia
using JumpProcesses
rate(u, p, t) = u[1]
affect(integrator) = (integrator.u[1] -= 1; nothing)
crj = ConstantRateJump(rate, affect)
dprob = DiscreteProblem([10], (0.0, 10.0))
jprob = JumpProblem(dprob, crj)
sol = solve(jprob, SSAStepper())
```

--------------------------------

### Install DifferentialEquations and Plots Packages

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/jump_diffusion.md

Installs the necessary Julia packages for solving differential equations and plotting. These packages are required for running the examples in this tutorial.

```julia
using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")
```

--------------------------------

### Install JumpProcesses and DifferentialEquations Packages

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Installs the necessary Julia packages for simulating jump processes and differential equations. This includes DifferentialEquations.jl (which bundles JumpProcesses.jl), Plots.jl for visualization, and optionally Catalyst.jl for chemical kinetics modeling.

```julia
using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")
Pkg.add("Catalyst") # optional
```

--------------------------------

### Setup for Tutorial (DifferentialEquations)

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Sets up the environment for the tutorial by loading necessary packages and configuring default plot settings. This code block is intended for use with a documentation generation system like Documenter.jl.

```julia
@setup tut2
using DifferentialEquations, Plots, LinearAlgebra
default(; lw = 2)
```

--------------------------------

### ConstantRateJump Example: Poisson Birth-Death Process

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Demonstrates the usage of ConstantRateJump for a simple Poisson birth-death process. It defines birth and death rates and their corresponding state modification functions. The example sets up a DiscreteProblem, builds a JumpProblem with the Direct algorithm, and solves it using SSAStepper.

```julia
using JumpProcesses

# Define a simple Poisson birth-death process
# Birth rate is constant λ, death rate is μ * population
birth_rate(u, p, t) = p.λ
birth_affect!(integrator) = (integrator.u[1] += 1)
birth_jump = ConstantRateJump(birth_rate, birth_affect!)

death_rate(u, p, t) = p.μ * u[1]
death_affect!(integrator) = (integrator.u[1] -= 1)
death_jump = ConstantRateJump(death_rate, death_affect!)

# Set up the problem
u₀ = [50]                          # Initial population
p = (λ = 2.0, μ = 0.02)            # Birth rate λ, death rate per individual μ
tspan = (0.0, 100.0)

# Create discrete problem (pure jump, no ODE/SDE dynamics)
dprob = DiscreteProblem(u₀, tspan, p)

# Build jump problem with Direct (Gillespie) algorithm
jprob = JumpProblem(dprob, Direct(), birth_jump, death_jump)

# Solve using SSAStepper for pure jump problems
sol = solve(jprob, SSAStepper())

# Access solution at specific times
println("Population at t=50: ", sol(50.0))
println("Final population: ", sol[end])
```

--------------------------------

### JumpProblem Setup and Aggregators in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

This section illustrates setting up a JumpProblem with different aggregator algorithms (Direct, SortingDirect, RSSA) for performance tuning. It covers pure jump problems and hybrid jump-ODE systems.

```julia
using JumpProcesses, OrdinaryDiffEq

# Pure jump problem (DiscreteProblem + SSAStepper)
rate1(u, p, t) = p[1] * u[1] * u[2]
affect1!(integrator) = (integrator.u[1] -= 1; integrator.u[2] += 1)
jump1 = ConstantRateJump(rate1, affect1!)

rate2(u, p, t) = p[2] * u[2]
affect2!(integrator) = (integrator.u[2] -= 1; integrator.u[3] += 1)
jump2 = ConstantRateJump(rate2, affect2!)

u₀ = [999, 1, 0]
p = [0.001, 0.1]
tspan = (0.0, 100.0)

# Different aggregator options for performance tuning
dprob = DiscreteProblem(u₀, tspan, p)

# Direct (Gillespie): Good for small systems (<10 jumps)
jprob_direct = JumpProblem(dprob, Direct(), jump1, jump2)

# SortingDirect: Better for medium systems (10-100 jumps)
jprob_sorting = JumpProblem(dprob, SortingDirect(), jump1, jump2;
                            dep_graph = [[1,2], [1,2]])

# RSSA: Best for large systems with many species updates
jprob_rssa = JumpProblem(dprob, RSSA(), jump1, jump2;
                         vartojumps_map = [[1,2], [1,2], Int[]],
                         jumptovars_map = [[1,2], [2,3]])

# Solve pure jump problems with SSAStepper
sol = solve(jprob_direct, SSAStepper())

# Control saving behavior (reduce memory for large simulations)
jprob_sparse = JumpProblem(dprob, Direct(), jump1, jump2;
                           save_positions = (false, false))
sol_sparse = solve(jprob_sparse, SSAStepper(); saveat = 1.0)

# Hybrid jump-ODE (piecewise deterministic Markov process)
function f!(du, u, p, t)
    du[4] = -0.1 * u[4]  # Exponential decay of 4th component
end
oprob = ODEProblem(f!, [999.0, 1.0, 0.0, 100.0], tspan, p)
jprob_hybrid = JumpProblem(oprob, Direct(), jump1, jump2)
sol_hybrid = solve(jprob_hybrid, Tsit5())  # Use ODE solver for hybrid problems

```

--------------------------------

### Add JumpProcesses and Plots Packages

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/simple_poisson_process.md

Installs the necessary JumpProcesses and Plots packages for simulating and visualizing jump processes. This is a prerequisite for running the examples.

```julia
using Pkg
Pkg.add("JumpProcesses")
Pkg.add("Plots")
```

--------------------------------

### Install JumpProcesses.jl Package

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/README.md

Installs the JumpProcesses.jl library using the Julia package manager. This is a lighter dependency than the meta DifferentialEquations.jl package.

```julia
using Pkg
Pkg.add("JumpProcesses")
```

--------------------------------

### Setup Jump Problem with Different Aggregators in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Demonstrates setting up a `JumpProblem` using various aggregator algorithms provided by JumpProcesses.jl, including Direct, NRM, RSSA, and their variants. It shows how to define reaction systems, initial conditions, and parameters, and how to optionally provide dependency graphs or variable-to-jump mappings.

```julia
using JumpProcesses

# Setup a test system
rs = [[1 => 1, 2 => 1], [2 => 1]]
ns = [[1 => -1, 2 => 1], [2 => -1, 3 => 1]]
maj = MassActionJump(rs, ns; param_idxs = [1, 2])

u₀ = [1000, 10, 0]
p = [0.0001, 0.05]
tspan = (0.0, 200.0)
dprob = DiscreteProblem(u₀, tspan, p)

# Direct: O(n) per step, simple implementation
jprob_direct = JumpProblem(dprob, Direct(), maj)

# DirectFW: Direct with FunctionWrappers (can reduce compilation)
jprob_dfw = JumpProblem(dprob, DirectFW(), maj)

# FRM: First Reaction Method - O(n) per step
jprob_frm = JumpProblem(dprob, FRM(), maj)

# NRM: Next Reaction Method - O(log n) per step with dependency graph
dep_graph = [[1, 2], [1, 2]]
jprob_nrm = JumpProblem(dprob, NRM(), maj; dep_graph = dep_graph)

# SortingDirect: Maintains sorted rates - good for heterogeneous rates
jprob_sd = JumpProblem(dprob, SortingDirect(), maj; dep_graph = dep_graph)

# DirectCR: Composition-Rejection - O(1) average for sparse networks
jprob_dcr = JumpProblem(dprob, DirectCR(), maj; dep_graph = dep_graph)

# RSSA: Rejection SSA - excellent for large networks
vartojumps = [[1, 2], [1, 2], Int[]]  # Which jumps depend on each variable
jumptovars = [[1, 2], [2, 3]]         # Which variables each jump modifies
jprob_rssa = JumpProblem(dprob, RSSA(), maj;
                         vartojumps_map = vartojumps,
                         jumptovars_map = jumptovars)

# RSSACR: RSSA with Composition-Rejection
jprob_rssacr = JumpProblem(dprob, RSSACR(), maj;
                           vartojumps_map = vartojumps,
                           jumptovars_map = jumptovars)
```

--------------------------------

### Set Initial Conditions and Time Span

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Initializes the starting count `u0` to `[0]` and defines the simulation time span `tspan` from 0.0 to 10.0. These parameters are essential for setting up the simulation problem.

```julia
u0 = [0]
tspan = (0.0, 10.0)
```

--------------------------------

### Optimized Species Index Ordering for MassActionJump in Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/jump_types.md

Provides an example of how to order species indices from smallest to largest within stoichiometry vectors for MassActionJump in Julia. This ordering is recommended for performance reasons.

```julia
reactant_stoich = [[1 => 2, 3 => 1, 4 => 2], [2 => 2, 3 => 2]]
```

--------------------------------

### Add OrdinaryDiffEq Package

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/simple_poisson_process.md

Installs the OrdinaryDiffEq package, which is required for simulating jump processes coupled with ODEs. This is a prerequisite for using ODEProblem and its solvers.

```julia
using Pkg
Pkg.add("OrdinaryDiffEq")
```

--------------------------------

### Instantiate JumpProblem with Direct Aggregator

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/jump_types.md

Example of creating a JumpProblem using Gillespie's Direct SSA method as the aggregator. This demonstrates how to pass an instantiated aggregator type to the JumpProblem constructor.

```julia
JumpProblem(prob, Direct(), jump1, jump2)
```

--------------------------------

### VariableRateJump Example: Time-Varying Birth Rate

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Introduces VariableRateJump for handling jumps with rates that change continuously over time. This example outlines a birth-death process with a seasonally varying birth rate, defined by a sinusoidal function. It requires integration with OrdinaryDiffEq for time-stepping.

```julia
using JumpProcesses, OrdinaryDiffEq

# Example: Birth-death with time-varying birth rate (seasonal)
# Birth rate oscillates: λ(t) = λ₀ * (sin(πt/12) + 1.5)


```

--------------------------------

### Solve Jump-ODE Example with ConstantRateJump

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/README.md

Solves an ODE for exponential growth coupled with a constant rate jump process that halves the solution. It defines the ODE function, initial conditions, time span, and the jump's rate and affect functions. The JumpProblem is then solved using the Tsit5 method.

```julia
using DifferentialEquations, Plots

# du/dt = u is the ODE part
function f(du, u, p, t)
    du[1] = u[1]
end
u₀ = [0.2]
tspan = (0.0, 10.0)
prob = ODEProblem(f, u₀, tspan)

# jump part

# fires with a constant intensity of 2
rate(u, p, t) = 2

# halve the solution when firing
affect!(integrator) = (integrator.u[1] = integrator.u[1] / 2)
jump = ConstantRateJump(rate, affect!)

# use the Direct method to handle simulating the jumps
jump_prob = JumpProblem(prob, Direct(), jump)

# now couple to the ODE, solving the ODE with the Tsit5 method
sol = solve(jump_prob, Tsit5())
plot(sol)
```

--------------------------------

### Display Package Status - Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/index.md

Displays the current status of packages in the Julia environment. This is useful for understanding installed packages and their versions for reproducibility.

```julia
using Pkg # hide
Pkg.status() # hide
```

--------------------------------

### Construct ODEProblem and JumpProblem for Simulation

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/simple_poisson_process.md

Constructs an ODEProblem with a trivial ODE and then a JumpProblem using the Direct aggregator. This setup is used to simulate jump processes coupled with a constant state, leveraging ODE solver time steppers.

```julia
u₀ = [0.0, 0.0]
oprob = ODEProblem(f!, u₀, tspan, p)
jprob = JumpProblem(oprob, Direct(), vrj2, deathvrj)
```

--------------------------------

### Setup RegularJump for SIR Model Tau-Leaping in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Sets up a `RegularJump` for simulating an SIR model using tau-leaping. It defines the rate function, stoichiometry matrix, and state change function. The `JumpProblem` is then created with `PureLeaping` and solved using `SimpleTauLeaping` for fixed-step approximation.

```julia
using JumpProcesses

# Define RegularJump for SIR model
# rate: computes all jump rates into output vector
function rate!(out, u, p, t)
    out[1] = p[1] * u[1] * u[2]  # Infection rate
    out[2] = p[2] * u[2]          # Recovery rate
    nothing
end

# Net stoichiometry matrix: change[species, reaction]
# Columns are reactions, rows are species changes
stoich_matrix = [-1  0;   # S changes: -1 for infection, 0 for recovery
                  1 -1;   # I changes: +1 for infection, -1 for recovery
                  0  1]   # R changes: 0 for infection, +1 for recovery

# change!: given jump counts, compute state change
function change!(du, u, p, t, counts, mark)
    # counts[i] = number of times reaction i occurred this step
    mul!(du, stoich_matrix, counts)
    nothing
end

# Create RegularJump (2 reactions)
num_reactions = 2
regular_jump = RegularJump(rate!, change!, num_reactions)

# Setup problem
u₀ = [10000.0, 100.0, 0.0]  # Must be Float64 for tau-leaping
p = [0.00005, 0.1]
tspan = (0.0, 100.0)

dprob = DiscreteProblem(u₀, tspan, p)
jprob = JumpProblem(dprob, PureLeaping(), regular_jump)

# Solve with SimpleTauLeaping (fixed step size)
sol_simple = solve(jprob, SimpleTauLeaping(); dt = 0.01)
```

--------------------------------

### Initialize SciMLPointProcess as a Hawkes Process (Julia)

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/applications/advanced_point_process.md

Demonstrates the initialization of a `SciMLPointProcess` object to represent a Hawkes process. This involves defining the graph structure, mark distributions, jump definitions, and time span. The `SciMLPointProcess` is then constructed with these components and the appropriate parameter initialization function.

```julia
using Graphs
using Distributions
V = 10
G = erdos_renyi(V, 0.2)
g = [[[i]; neighbors(G, i)] for i in 1:nv(G)]
mark_dist = [MvNormal(rand(2), [0.2, 0.2]) for i in 1:nv(G)]
jumps = [hawkes_jump(i, g, mark_dist) for i in 1:nv(G)]
tspan = (0.0, 50.0)
hawkes = SciMLPointProcess{
    Vector{Real},
    eltype(jumps),
    typeof(g),
    eltype(mark_dist),
    eltype(tspan)
}(jumps,
    mark_dist,
    g,
    hawkes_p,
    tspan[1],
    tspan[2])
```

--------------------------------

### Initialize DiscreteProblem (Julia)

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Sets up the initial state and time span for a discrete problem, which serves as the base for a JumpProblem. The state `u0` is a vector representing the counts of the multivariate TPPs.

```julia
u0 = [0, 0]
tspan = (0.0, 10.0)
dprob = DiscreteProblem(u0, tspan)
```

--------------------------------

### MassActionJump Example: SIR Epidemic Model

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Illustrates the MassActionJump type for an SIR epidemic model. This specialized type optimizes performance by defining reactions through stoichiometry rather than individual rate and affect functions. The example defines reactant and net stoichiometries, parameter indices, and then sets up and solves the jump problem.

```julia
using JumpProcesses

# SIR epidemic model: S + I -> 2I (infection), I -> R (recovery)
# Species: S=1, I=2, R=3
# Parameters: β (infection rate), ν (recovery rate)

# Define reaction stoichiometries
# Reactant stoichiometry: which species are consumed (species => count)
reactant_stoich = [
    [1 => 1, 2 => 1],    # S + I -> ... (consumes 1 S and 1 I)
    [2 => 1]             # I -> ...     (consumes 1 I)
]

# Net stoichiometry: net change in each species (species => change)
net_stoich = [
    [1 => -1, 2 => 1],   # S decreases by 1, I increases by 1
    [2 => -1, 3 => 1]    # I decreases by 1, R increases by 1
]

# Parameter indices: which parameters are the rate constants
param_idxs = [1, 2]  # p[1] = β, p[2] = ν

# Create mass action jump (handles all reactions)
mass_jump = MassActionJump(reactant_stoich, net_stoich; param_idxs = param_idxs)

# Set up problem
u₀ = [990, 10, 0]                     # [S, I, R] initial populations
p = [0.1/1000, 0.05]                  # [β, ν] rate constants
tspan = (0.0, 250.0)

dprob = DiscreteProblem(u₀, tspan, p)
jprob = JumpProblem(dprob, Direct(), mass_jump)
sol = solve(jprob, SSAStepper())

# Plot results
using Plots
plot(sol, label=["S" "I" "R"], xlabel="Time", ylabel="Population",
     title="SIR Model via MassActionJump")
```

--------------------------------

### Creating and Using JumpSet in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Illustrates how to create a `JumpSet` in Julia to organize different types of jumps (MassActionJump, ConstantRateJump) for use in a `JumpProblem`. It shows both positional and keyword argument methods for constructing `JumpSet` and how to use it with `JumpProblem` and `solve`.

```julia
using JumpProcesses

# Create individual jumps
rate1(u, p, t) = p[1] * u[1]
affect1!(integrator) = (integrator.u[1] -= 1; integrator.u[2] += 1)
cj1 = ConstantRateJump(rate1, affect1!)

rate2(u, p, t) = p[2] * u[2]
affect2!(integrator) = (integrator.u[2] -= 1)
cj2 = ConstantRateJump(rate2, affect2!)

# Create a MassActionJump
rs = [[1 => 2]]  # Dimerization: 2A -> ...
ns = [[1 => -2, 3 => 1]]  # ... -> A₂
maj = MassActionJump(rs, ns; param_idxs = [3])

# Combine into JumpSet (positional arguments)
jset = JumpSet(maj, cj1, cj2)

# Or use keyword arguments for collections
cj_vec = [cj1, cj2]
jset2 = JumpSet(; massaction_jump = maj, constant_jumps = cj_vec)

# Use JumpSet in JumpProblem
u₀ = [100, 0, 0]
p = [0.1, 0.05, 0.001]
tspan = (0.0, 50.0)

dprob = DiscreteProblem(u₀, tspan, p)
jprob = JumpProblem(dprob, Direct(), jset)
sol = solve(jprob, SSAStepper())

# Dynamically build jump collections
constant_jumps = ConstantRateJump[]
for i in 1:5
    rate_i = (u, p, t) -> p[i] * u[i]
    affect_i! = integrator -> (integrator.u[i] -= 1)
    push!(constant_jumps, ConstantRateJump(rate_i, affect_i!))
end

jset_dynamic = JumpSet(; constant_jumps = constant_jumps)
```

--------------------------------

### Create Graph Topology in Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/spatial.md

Creates an unstructured mesh topology using an AbstractGraph from the Graphs.jl package. An example shows creating a directed cyclic graph.

```julia
using Graphs
graph = cycle_digraph(5) # directed cyclic graph on 5 nodes

```

--------------------------------

### Construct and Solve JumpProblem (Julia)

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Creates a JumpProblem using the DiscreteProblem, the Coevolve algorithm, the defined TPPs, and the dependency graph. It then solves the problem using the SSAStepper and plots the results.

```julia
jprob = JumpProblem(dprob, Coevolve(), poisson_process, seasonal_process; dep_graph)
sol = solve(jprob, SSAStepper())
plot(sol, labels = ["N_1(t)" "N_2(t)"], xlabel = "t", legend = :topleft)
```

--------------------------------

### Load JumpProcesses and Plots Packages

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Loads the necessary JumpProcesses and Plots packages for simulating and visualizing the Poisson process. These are fundamental for the subsequent steps.

```julia
using JumpProcesses, Plots
```

--------------------------------

### Create RegularJump Object

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Instantiates a RegularJump object using the previously defined 'rate' and 'change' functions, specifying the number of jumps being encoded (2 in this SIR model example).

```julia
rj = RegularJump(rate, change, 2)
```

--------------------------------

### Simulate Hybrid System with Multiple Jumps

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Demonstrates simulating a hybrid system by combining a `MassActionJump` (e.g., `mass_act_jump`) with a newly defined `ConstantRateJump` (`birth_jump`) using `JumpProblem`. The solution is then solved and plotted.

```julia
jump_prob = JumpProblem(prob, Direct(), mass_act_jump, birth_jump)
sol = solve(jump_prob, SSAStepper())
plot(sol; label = ["S(t)" "I(t)" "R(t)"])
```

--------------------------------

### Define Jumps and JumpSet with Dependency Ordering (Julia)

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/faq.md

Illustrates the definition of `MassActionJump` and `ConstantRateJump` and their aggregation into a `JumpSet`, highlighting the internal ordering of jumps for SSAs.

```julia
using JumpProcesses
rs = [[1 => 1], [2 => 1]]
ns = [[1 => -1, 2 => 1], [1 => 1, 2 => -1]]
p = [1.0, 0.0]
maj = MassActionJump(rs, ns; param_idxs = [1, 2])
rate1(u, p, t) = u[1]
function affect1!(integrator)
    u[1] -= 1
end
cj1 = ConstantRateJump(rate1, affect1)
rate2(u, p, t) = u[2]
function affect2!(integrator)
    u[2] -= 1
end
cj2 = ConstantRateJump(rate2, affect2)
jset = JumpSet(; constant_jumps = [cj1, cj2], massaction_jump = maj)
```

--------------------------------

### Display Version Information - Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/index.md

Shows detailed information about the Julia version, operating system, and hardware. This is crucial for reproducing computational environments.

```julia
using InteractiveUtils # hide
versioninfo() # hide
```

--------------------------------

### SIR Model Simulation with JumpProcesses.jl

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/README.md

Demonstrates simulating a stochastic chemical kinetics SIR model using JumpProcesses.jl. This example involves three species (S, I, R) and two reaction types represented as jump processes.

```julia
using JumpProcesses, Plots

# here we order S = 1, I = 2, and R = 3

```

--------------------------------

### SSAStepper for Pure Jump Problems in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

This code demonstrates the basic usage of SSAStepper for solving pure jump problems defined with DiscreteProblem and ConstantRateJump. It shows how to use saveat for specific time points and EnsembleProblem for multiple simulations.

```julia
using JumpProcesses

# Basic usage with SSAStepper
rate(u, p, t) = p[1] * u[1]
affect!(integrator) = (integrator.u[1] -= 1)
jump = ConstantRateJump(rate, affect!)

u₀ = [100]
p = [0.1]
tspan = (0.0, 50.0)

dprob = DiscreteProblem(u₀, tspan, p)
jprob = JumpProblem(dprob, Direct(), jump)

# Solve with SSAStepper
sol = solve(jprob, SSAStepper())

# Use saveat to save at specific times
sol_saveat = solve(jprob, SSAStepper(); saveat = 0.5)

# Multiple simulations with EnsembleProblem
using DiffEqBase
ensemble_prob = EnsembleProblem(jprob)
ensemble_sol = solve(ensemble_prob, SSAStepper(), EnsembleThreads();
                     trajectories = 1000)

```

--------------------------------

### Animate Spatial Jump Solutions in Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/spatial.md

Generates an animation of the solution to a spatial jump process. It defines functions to get individual frames and create the animation object. Requires the Plots.jl package for plotting and animation.

```julia
using Plots
is_static(spec) = (spec == 3) # true if spec does not hop
"""
get frame k
"""
function get_frame(k, sol, linear_size, labels, title)
    num_species = length(labels)
    h = 1 / linear_size
    t = sol.t[k]
    state = sol.u[k]
    xlim = (0, 1 + 3h / 2)
    ylim = (0, 1 + 3h / 2)
    plt = plot(xlim = xlim, ylim = ylim, title = "$title, $(round(t, sigdigits=3)) seconds")

    species_seriess_x = [[] for i in 1:num_species]
    species_seriess_y = [[] for i in 1:num_species]
    CI = CartesianIndices((linear_size, linear_size))
    for ci in CartesianIndices(state)
        species, site = Tuple(ci)
        x, y = Tuple(CI[site])
        num_molecules = state[ci]
        sizehint!(species_seriess_x[species], num_molecules)
        sizehint!(species_seriess_y[species], num_molecules)
        if !is_static(species)
            randsx = rand(num_molecules)
            randsy = rand(num_molecules)
        else
            randsx = zeros(num_molecules)
            randsy = zeros(num_molecules)
        end
        for k in 1:num_molecules
            push!(species_seriess_x[species], x * h - h / 4 + 0.5h * randsx[k])
            push!(species_seriess_y[species], y * h - h / 4 + 0.5h * randsy[k])
        end
    end
    for species in 1:num_species
        scatter!(plt, species_seriess_x[species], species_seriess_y[species],
            label = labels[species], marker = 6)
    end
    xticks!(plt, range(xlim..., length = linear_size + 1))
    yticks!(plt, range(ylim..., length = linear_size + 1))
    xgrid!(plt, 1, 0.7)
    ygrid!(plt, 1, 0.7)
    return plt
end

"""
make an animation of solution sol in 2 dimensions
"""
function animate_2d(sol, linear_size; species_labels, title, verbose = true)
    num_frames = length(sol.t)
    anim = @animate for k in 1:num_frames
        verbose && println("Making frame $k")
        get_frame(k, sol, linear_size, species_labels, title)
    end
    anim
end
# animate
anim = animate_2d(solution, 5, species_labels = ["A", "B", "C"], title = "A + B <--> C",
    verbose = false)
fps = 5
gif(anim, fps = fps)

```

--------------------------------

### Define Dependency Graph and JumpProblem

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Specifies the dependency graph for jump recalculations and constructs a `JumpProblem` using the `Coevolve` aggregator. This setup is necessary for correctly handling interactions between jumps, especially with bounded `VariableRateJump`s.

```julia
dep_graph = [[1, 2], [1, 2]]
prob = DiscreteProblem(u₀, tspan, p1)
jump_prob = JumpProblem(prob, Coevolve(), jump3, jump4; dep_graph)
```

--------------------------------

### Simulating Jump Processes with Callbacks in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Demonstrates how to use callbacks to modify jump process parameters during simulation. It shows how to define a condition and an affect function for the callback, and how to integrate it into the `solve` function. This is useful for simulating systems where parameters change at specific time points.

```julia
using JumpProcesses

# Assuming jprob is a pre-defined JumpProblem

# Remake with new timespan
jprob4 = remake(jprob; tspan = (0.0, 200.0))
sol4 = solve(jprob4, SSAStepper())

# Using callbacks with jumps
# Example: Change parameters at t=50
condition(u, t, integrator) = t == 50.0

function affect_callback!(integrator)
    integrator.p[1] = 0.0        # Turn off infections
    integrator.p[2] = 0.2        # Speed up recovery
    reset_aggregated_jumps!(integrator)  # REQUIRED after modifying p!
    nothing
end

callback = DiscreteCallback(condition, affect_callback!)

sol_cb = solve(jprob, SSAStepper();
               callback = callback,
               tstops = [50.0])  # Must specify tstop for exact callback timing

# Parameter sweep using remake
results = []
for β in [0.00005, 0.0001, 0.0002, 0.0005]
    jprob_sweep = remake(jprob; p = [β, 0.05])
    sol_sweep = solve(jprob_sweep, SSAStepper())
    push!(results, (β = β, peak_infected = maximum(sol_sweep[2, :])))
end
```

--------------------------------

### Julia ConstantRateJump and MassActionJump Rate Functions

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Example rate functions for ConstantRateJump and MassActionJump in Julia. These functions must remain constant between consecutive jumps. They depend on species/states (u) and parameters (p) at a given time (t).

```julia
rate1(u, p, t) = p[1] * u[1] * u[2]
rate2(u, p, t) = p[2] * u[2]
```

--------------------------------

### Define Bounded Variable Rate Jump Parameters (Julia)

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Sets up the upper bound, the time interval for the upper bound's validity, and the lower bound for a VariableRateJump. This allows for efficient simulation of TPPs with predictable rate fluctuations.

```julia
urate(u, p, t) = 2           # upper bound
rateinterval(u, p, t) = Inf  # time window bound is valid over
lrate(u, p, t) = 1           # lower bound
```

--------------------------------

### Define Variable Rate Jump Function

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/jump_diffusion.md

Defines a function for the rate of a variable rate jump. The rate is dependent on the solution vector `u`, parameters `p`, and time `t`. In this example, the rate is set to the first element of the solution vector.

```julia
rate(u, p, t) = u[1]
```

--------------------------------

### Simulate Jump Process with SSAStepper

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/simple_poisson_process.md

Solves a JumpProblem using the SSAStepper, which is suitable for discrete jump processes. The solution can then be plotted to visualize the system's state over time.

```julia
sol = solve(jprob, SSAStepper())
plot(sol, labels = ["N(t)" "D(t)"], xlabel = "t", legend = :topleft)
```

--------------------------------

### Load Packages and Set Plot Defaults

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Loads the required Julia packages for the tutorial: DifferentialEquations, Plots, and LinearAlgebra. It also sets a default line width for plots to 2.

```julia
using DifferentialEquations, Plots, LinearAlgebra
default(; lw = 2)
```

--------------------------------

### Compare Tau-Leaping with Exact Solution in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Compares the results of a tau-leaping simulation (`sol_simple`) with an exact simulation (`sol_exact`) of an SIR model. It sets up an exact `JumpProblem` using `MassActionJump` and `Direct` aggregator, then plots both solutions for visual comparison.

```julia
using StochasticDiffEq
# For adaptive tau-leaping, use StochasticDiffEq
# sol_adaptive = solve(jprob, TauLeaping())

# Compare with exact solution
dprob_exact = DiscreteProblem([10000, 100, 0], tspan, p)
rs = [[1 => 1, 2 => 1], [2 => 1]]
ns = [[1 => -1, 2 => 1], [2 => -1, 3 => 1]]
maj = MassActionJump(rs, ns; param_idxs = [1, 2])
jprob_exact = JumpProblem(dprob_exact, Direct(), maj)
sol_exact = solve(jprob_exact, SSAStepper())

using Plots
plot(sol_simple, label=["S (τ-leap)" "I (τ-leap)" "R (τ-leap)"], linestyle=:dash)
plot!(sol_exact, label=["S (exact)" "I (exact)" "R (exact)"])
```

--------------------------------

### Solve and Plot Poisson Process Simulation

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Solves the `JumpProblem` `jprob` using the `SSAStepper` time-stepper and plots the resulting solution `sol`. This generates and visualizes one realization of the homogeneous Poisson process.

```julia
sol = solve(jprob, SSAStepper())
plot(sol)
```

--------------------------------

### Initialize and Reset Parameters for SciMLPointProcess in Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/applications/advanced_point_process.md

Defines methods for initializing and resetting parameters of a SciMLPointProcess. The `params` function retrieves the current parameters, while `params!` updates them. This flexibility is crucial for general TPP modeling and parameter estimation.

```julia
params(pp::SciMLPointProcess) = pp.p(pp)
params!(pp::SciMLPointProcess, p) = pp.p(pp, p)
nothing # hide
```

--------------------------------

### Simulate SIR Model with MassActionJump

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/README.md

Simulates an SIR model using MassActionJump, defining substrate stoichiometry, net change, and rate constants. It then sets up a DiscreteProblem and solves it as a pure jump process using the Direct method.

```julia
# substrate stoichiometry:
substoich = [[1 => 1, 2 => 1],    # 1*S + 1*I
    [2 => 1]]                     # 1*I
# net change by each jump type
netstoich = [[1 => -1, 2 => 1],   # S -> S-1, I -> I+1
    [2 => -1, 3 => 1]]            # I -> I-1, R -> R+1
# rate constants for each jump
p = (0.1 / 1000, 0.01)

# p[1] is rate for S+I --> 2I, p[2] for I --> R
pidxs = [1, 2]

maj = MassActionJump(substoich, netstoich; param_idxs = pidxs)

u₀ = [999, 1, 0]       #[S(0),I(0),R(0)]
tspan = (0.0, 250.0)
dprob = DiscreteProblem(u₀, tspan, p)

# use the Direct method to simulate
jprob = JumpProblem(dprob, maj)

# solve as a pure jump process, i.e. using SSAStepper
sol = solve(jprob)
plot(sol)
```

--------------------------------

### Benchmark Jump Aggregators with BenchmarkTools in Julia

Source: https://context7.com/sciml/jumpprocesses.jl/llms.txt

Benchmarks the performance of different jump aggregator algorithms using the BenchmarkTools.jl package. It solves the jump problem using `SSAStepper` for each specified aggregator and reports the execution time.

```julia
using BenchmarkTools
@btime solve($jprob_direct, SSAStepper())
@btime solve($jprob_nrm, SSAStepper())
@btime solve($jprob_rssa, SSAStepper())
```

--------------------------------

### Create Discrete Problem for TPP

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/point_process_simulation.md

Creates a `DiscreteProblem` using the initial conditions `u0` and time span `tspan`. This serves as the base problem for simulating processes that evolve in discrete time steps, like TPPs.

```julia
dprob = DiscreteProblem(u0, tspan)
```

--------------------------------

### Define Zero Order Reaction with MassActionJump in Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/jump_types.md

Illustrates two methods for defining a zero-order reaction (e.g., \O \overset{k}{\rightarrow} A) using MassActionJump in Julia. It shows how to represent the reactant stoichiometry using either an explicit zero mapping or an empty vector.

```julia
p = [1.0]
reactant_stoich = [[0 => 1]]
net_stoich = [[1 => 1]]
jump = MassActionJump(reactant_stoich, net_stoich; param_idxs = [1])
```

```julia
p = [1.0]
reactant_stoich = [Vector{Pair{Int, Int}}()]
net_stoich = [[1 => 1]]
jump = MassActionJump(reactant_stoich, net_stoich; param_idxs = [1])
```

--------------------------------

### Load OrdinaryDiffEq Package

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/simple_poisson_process.md

Loads the OrdinaryDiffEq package into the current Julia session. This makes its functionalities available for use in defining and solving problems.

```julia
using OrdinaryDiffEq
```

--------------------------------

### Solve Jump Problem with SSAStepper in Julia

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/discrete_stochastic_example.md

Solves the `JumpProblem` using the `SSAStepper` for efficient simulation of pure jump processes. The resulting solution object `sol` can be plotted using standard plotting recipes.

```julia
sol = solve(jump_prob, SSAStepper())
plot(sol, label = ["S(t)" "I(t)" "R(t)"])
```

--------------------------------

### Create SDEProblem for Jump Diffusion

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/jump_diffusion.md

Initializes an `SDEProblem` with a drift function `f`, a noise function `g`, an initial state `[0.2]`, and a time span `(0.0, 10.0)`. This sets up the continuous part of the jump diffusion process.

```julia
prob = SDEProblem(f, g, [0.2], (0.0, 10.0))
```

--------------------------------

### Load Packages and Set Plotting Defaults

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/jump_diffusion.md

Loads the DifferentialEquations and Plots packages into the current Julia session and sets default line width for plots. This prepares the environment for creating visualizations.

```julia
using DifferentialEquations, Plots
default(; lw = 2)
```

--------------------------------

### Set up JumpProblem with Next Subvolume Method

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/tutorials/spatial.md

Configures the JumpProblem using the Next Subvolume Method (NSM) for spatial simulations. It integrates the discrete problem, jump processes, hopping constants, and spatial grid.

```julia
alg = NSM()
jump_prob = JumpProblem(prob, alg, majumps, hopping_constants = hopping_constants,
    spatial_system = grid, save_positions = (true, false))
```

--------------------------------

### JumpProblem Configuration with Dependency Graph

Source: https://github.com/sciml/jumpprocesses.jl/blob/master/docs/src/jump_types.md

Demonstrates how to pass a dependency graph to JumpProblem for aggregators like DirectCR. This is necessary when using ConstantRateJump or bounded VariableRateJump without auto-generation from a Catalyst reaction network.

```julia
JumpProblem(prob, DirectCR(), jump1, jump2; dep_graph = your_dependency_graph)
```