GitHub

Continuous-Time Sampling Algorithms

Continuous-time Monte Carlo in MonteCarloX is event-driven. You sample when the next event happens and which event happens.

Why this differs from Metropolis-style sampling

  • Importance sampling focuses on stationary distributions.
  • Continuous-time sampling focuses on trajectories in real/simulation time.

Conceptually, trajectories can be viewed as states in path space. The continuous-time interface is a computational specialization for explicit event-time dynamics.

Path-space notation:

  • trajectory $\omega = (x_t)_{t \in [0,T]}$
  • target path distribution $\pi(\omega) = p(\omega \mid \lambda)$

Gillespie workflow

At each step:

  1. compute total rate
  2. sample waiting time $\mathrm{d}t \sim \mathrm{Exp}(\lambda_\mathrm{tot})$
  3. sample event index proportional to rates
  4. advance time and apply event update

Gillespie stores steps and current time.

Minimal loop

using Random
using MonteCarloX

alg = Gillespie(MersenneTwister(10))
rates = [0.5, 1.0]  # e.g. birth, death
for _ in 1:10000
    t, event = step!(alg, rates)
    # apply event to your state
    # update rates from new state
end

step! vs advance!

  • step!: one event at a time (full manual control)
  • advance!: run until total_time, with optional callbacks (measure!, update!)

Example with explicit state + event source:

using Random
using MonteCarloX

sys = Dict(:N => 30)
event_source(sys::Dict{Symbol,Int}) = [0.2 * sys[:N], 0.1 * sys[:N]]
alg = Gillespie(MersenneTwister(11))

modify_cb = (state, event, t) -> begin
    if event == 1
        state[:N] += 1
    elseif event == 2
        state[:N] = max(0, state[:N] - 1)
    end
end

advance!(alg, sys, 20.0; modify!=modify_cb)

Event sources supported

  • raw vectors of rates
  • StatsBase weights
  • rate-based event handlers (AbstractEventHandlerRate)
  • time-ordered event queues (AbstractEventHandlerTime)
  • callback rates_at_time(t)

Low-level building blocks

  • next_time: waiting-time sampling (homogeneous or thinning-based)
  • next_event: event-index sampling
  • next: one combined draw $(\mathrm{d}t,\, \text{event})$

Example map (algorithm ↔ application)

  • Birth-death dynamics (Gillespie): examples/stochastic_processes/gillespie_birth_death.ipynb
  • Reaction-network dynamics (Gillespie): examples/stochastic_processes/gillespie_dimerization.ipynb
  • Poisson primitives / low-level event draws: examples/stochastic_processes/kmc_poisson.ipynb

API reference

MonteCarloX.Gillespie — Type
Gillespie <: AbstractKineticMonteCarlo

Gillespie algorithm for continuous-time event-driven sampling.

Fields

  • rng::AbstractRNG: Random number generator
  • steps::Int: Number of events sampled
  • time::Float64: Current simulation time
source
MonteCarloX.next — Function
next(alg::AbstractKineticMonteCarlo, rates::AbstractVector)

Draw next event waiting time and event id from raw rates.

For zero total rate, returns (Inf, 0).

source
next(alg::AbstractKineticMonteCarlo, rates::AbstractVector)

Draw next event waiting time and event id from raw rates (StatsBase.AbstractWeights is also an AbstractVector).

For zero total rate, returns (Inf, 0).

source
next(alg::AbstractKineticMonteCarlo, event_handler::AbstractEventHandlerRate)

Draw next event waiting time and event from an event handler.

source
next(alg::AbstractKineticMonteCarlo, event_handler::AbstractEventHandlerTime)

Draw next event from a time-ordered event queue.

source

next(alg::AbstractKineticMonteCarlo, ratesattime::Function)

Draw next event waiting time and event from a time-dependent rates function.

rates_at_time is called as rates_at_time(alg.time) and should return either an AbstractVector, AbstractWeights, or AbstractEventHandlerRate.

source
MonteCarloX.step! — Function
step!(alg::AbstractKineticMonteCarlo, event_source)

Perform one kinetic Monte Carlo step from an event_source (rates or an event handler).

Updates alg.time and alg.steps internally and returns (t_new, event).

If dt=Inf then this is reflected in the algorithm.time and returned correspondingly If event is nothing, then this is returned as well.

source
MonteCarloX.next_time — Function
next_time(rng::AbstractRNG, rate_generation::Number)::Float64

Draw next waiting time for a homogeneous Poisson event stream.

source
next_time(rng::AbstractRNG, rate::Function, rate_generation::Real)::Real

Draw next waiting time for an inhomogeneous Poisson stream using thinning.

source
MonteCarloX.next_event — Function
next_event(rng::AbstractRNG, rates::Union{AbstractWeights, AbstractVector})::Int

Draw one event index proportional to rates.

source
MonteCarloX.advance! — Function
advance!(alg::AbstractKineticMonteCarlo, sys, total_time; t0=0, measure!=measure!, modify!=modify!)

Advance system sys using alg until total_time.

Loop order per step:

  1. step! — draw next event from event_source(sys)
  2. measure! — observe before modification: measure!(sys, event, t)
  3. modify! — apply event to state: modify!(sys, event, t)

Stops early if the event source is exhausted or total time is reached.

source
MonteCarloX.reset! — Method
reset!(alg::AbstractKineticMonteCarlo)

Reset kinetic Monte Carlo statistics (steps, time) to zero.

source