Low-level Sampler Interface

SSA{KeyType,TimeType}

This abstract type represents a stochastic simulation algorithm (SSA). It is parametrized by the clock ID, or key, and the type used for the time, which is typically a Float64. The type of the key can be anything you would use as a dictionary key. This excludes mutable values but includes a wide range of identifiers useful for simulation. For instance, it could be a String, but it could be a Tuple{Int64,Int64,Int64}, so that it indexes into a complicated simulation state.

enable!(sampler, clock, distribution, enablingtime, currenttime, RNG)

Tell the sampler to start a clock.

• sampler::SSA{KeyType,TimeType} - The sampler to tell.
• clock::KeyType - The ID of the clock. Can be a string, integer, tuple, etc.
• distribution::Distributions.UnivariateDistribution
• enablingtime::TimeType - The zero time for the clock's distribution, in absolute time. Usually equal to when.
• when::TimeType - The current time of the simulation.
• rng::AbstractRNG - A random number generator.

These times are absolute since the start of the simulation. The current time should be when. If you want to shift the distribution so that this event cannot happen for a little while then choose enablingtime > when. If you want to modify the distribution by shifting it left, then choose enablingtime < when. Usually, enablingtime == when. The truncated() function from Distributions.jl can make the left-shift conditioning explicit, but pair it with the shifted enabling time — enable!(sampler, clock, truncated(dist; lower=age), when - age, when, rng) — because a truncated distribution with enablingtime == when measures the truncated sample from now, pushing every firing at least age into the future.

reset!(sampler)

After a sampler is used for a simulation run, it has internal state. This function resets that internal state to the initial value in preparation for another sample run.

reset!(sampling)

Clear all enabled clocks (and delayed state if present) and reset time to the fixed start time.

disable!(sampler, clock, when)

Tell the sampler to forget a clock. We include the current simulation time because some Next Reaction methods use this to optimize sampling.

next(sampler, when, rng)

Ask the sampler for what happens next, in the form of (when, which)::Tuple{TimeType,KeyType}. rng is a random number generator.

The returned (when, which) is a reservation, not a commitment. It is valid only until the next enable!, disable!, or fire! call changes the sampler's state. Calling next twice without an intervening state change returns valid reservations, but they are not guaranteed to be identical across samplers: FirstReaction redraws on each call, while CombinedNextReaction returns the same cached reservation. Act on the most recent call's result. The two supported responses are to fire it—call fire! with the returned clock and time—or to decline it and stop the simulation (the fixed-horizon pattern). There is deliberately no peek.

The when argument must not decrease from one call to the next, and must never advance past a pending firing time without that event being fired. The time of the most recently fired event is always a safe choice. Within this invariant, re-querying next at a later when is legal — this is how MultiSampler composes sub-samplers: the global minimum event time cannot exceed any sub-sampler's pending minimum, so losing sub-samplers are re-queried in-contract at the advanced global time. Outside the invariant, samplers legitimately disagree: FirstReaction re-conditions every clock on survival to when, while CombinedNextReaction returns putative times cached at enabling. Within it they agree, because every surviving clock's putative time is at least when.

timetype(sampler)

Return the type of clock times.

timetype(sampling)

Time type T, usually Float64.

enabled(sampler)

Returns a read-only set of currently-enabled clocks.

length(sampler)::Int64

Return the number of stored clocks.

getindex(sampler, clock::KeyType)

Return stored state for a particular clock. If the clock does not exist, a KeyError will be thrown. Different samplers have different stored state.

keys(sampler)

Return all stored clocks as a vector.

clone(sampler)

Given an existing sampler, make a copy that has the same type and same constructor options but has no data in it. Use this to initialize an array of samplers.

clone(ds::DelayedState)

Deep copy of a DelayedState.

copy_clocks!(dst, src)

Copy enabled clocks and delayed state from src into dst, jittering the destination sampler so it will generate different next() samples.

jitter!(sampler, when, rng)

Takes a sampler that is at a statistically-valid stopping time (this is a technical term) and resamples all clocks currently-enabled so that this copy of the sampler will yield different results from another copy of the sampler despite its internal cache of clock data.

CombinedNextReaction{KeyType,TimeType}()

This combines Next Reaction Method and Modified Next Reaction Method. The Next Reaction Method is from Gibson and Bruck in their 2000 paper called "Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels". The Modified Next Reaction Method is from David F. Anderson's 2007 paper, "A modified Next Reaction Method for simulating chemical systems with time dependent propensities and delays". Both methods reuse draws of random numbers. The former works by accumulating survival of a distribution in a linear space and the latter works by accumulating survival of a distribution in a log space.

Each enabled clock specifies a univariate distribution from the Distributions package. Every distribution is more precise being sampled in the manner of the Next Reaction method (linear space) or the manner of the Modified Next Reaction method (log space). This sampler chooses which space to use depending on the type of the UnivariateDistribution and based on performance timings that are done during package testing. Defaults are set for those distributions included in the Distributions.jl package. If you want to add a distribution, then define:

sampling_space(::MyDistribution) = LogSampling

If you want to override a choice in the library, then create a sub-type of the given distribution, and specify its sampling space.

struct LinearGamma <: Distributions.Gamma end
sampling_space(::LinearGamma) = LinearSampling

If you want to test a distribution, look at tests/nrmetric.jl to see how distributions are timed.

DirectCall{KeyType,TimeType,TreeType}()

DirectCall is responsible for sampling among Exponential distributions. It samples using the Direct method. In this case, there is no optimization to that Direct method, so we call it DirectCall because it recalculates everything every time you call it.

The algorithm for the Direct Method relies heavily on what data structure it uses to maintain a list of hazard rates, such that it can know the sum of those hazards and index into them using a random value. This struct has a default constructor that chooses a data structure for you, but there are several options.

Example

If we know that our simulation will only use a small number of different clock keys, then it would make sense to use a data structure that disables clocks by zeroing them out, instead of removing them from the list. This will greatly reduce memory churn. We can do that by changing the underlying data structure.

prefix_tree = BinaryTreePrefixSearch{T}()
keyed_prefix_tree = KeyedKeepPrefixSearch{K,typeof(prefix_tree)}(prefix_tree)
sampler_noremove = DirectCall{K,T,typeof(keyed_prefix_tree)}(keyed_prefix_tree)

FirstReaction{KeyType,TimeType}()

This is the classic first reaction method for general distributions. Every time you sample, this goes to each distribution and asks when it would fire. Then it takes the soonest and throws out the rest of the sampled times until the next sample. It can also be very fast when there are only a few clocks to sample.

One interesting property of this sampler is that you can call next() multiple times in order to get a distribution of next firing clocks and their times to fire.

FirstToFire{KeyType,TimeType}()

This sampler is often the fastest for non-exponential distributions. When a clock is first enabled, this sampler asks the clock when it would fire and saves that time in a sorted heap of future times. Then it works through the heap, one by one. When a clock is disabled, its future firing time is removed from the list. There is no memory of previous firing times.

This uses a DataStructures.MutableBinaryMinHeap which is a Fibonacci heap. It has been tested against many other heaps and rarely loses by more than a few percent, so we are sticking with it.

PSSACR{KeyType,TimeType}(; ngroups::Int=64)

Exact, continuous-time sampler using composition–rejection over groups.

Assumptions:

• Each enabled clock is Exponential with rate λ (i.e., distribution Exponential(θ) where λ=1/θ).
• The sampler maintains group-wise sums and maxima of rates.
• Time to next event is drawn from Exp(∑λ), and the firing key is drawn by a two-stage selection: group by probability ∝ group-sum, then within-group by rejection from a uniform proposal with acceptance λ/λ_max_group.

Performance notes:

• Choose ngroups so that groups remain small and rate magnitudes similar. If you possess "owner" metadata (partial-propensity style), call assign_group! before enable! to place clocks with common owners into the same group.

References:

• R. Ramaswamy & I. F. Sbalzarini, "A partial-propensity variant of the composition-rejection SSA", J. Chem. Phys. 132, 044102 (2010).
• A. Slepoy, A. P. Thompson, S. J. Plimpton, "A constant-time kinetic Monte Carlo algorithm...", J. Chem. Phys. 128, 205101 (2008).

RSSA{KeyType,TimeType}(; bound_factor=1.05)

Rejection-based SSA with global Fenwick tree for candidate selection. This is for exponential distributions only not time-dependent rates.

• bound_factor ≥ 1.0 controls default upper bounds: \bar ai ← max(\bar ai, boundfactor * ai). Set to 1.0 for no rejections (reduces to direct-method timing with tree selection).

Algorithmic core:

• Maintain per-clock true rate ai and a certified upper bound ar ai >= a_i.
• Maintain Ā = sumi ar ai and a Fenwick tree over { ar a_i } for O(log N) sampling.
• Draw candidate times from Exp(Ā). Select candidate clock i by Categorical(ar a_i/Ā).
• Accept with probability ai / ar ai; otherwise reject and continue thinning.

Exactness: standard thinning of a Poisson process with rate Ā, with acceptance ai/ar ai, yields the target Markov jump process (homogeneous propensities). See Thanh et al. (2014, 2015).

Notes:

• This implementation targets time-homogeneous propensities (Exponential only). For time-dependent rates (tRSSA), you need piecewise-time envelopes and an integral sampler; these can be added as a thin extension without changing the public interface.

Petri{KeyType,TimeType}(dt=1.0)

A chaos-monkey sampler for smoke-testing. It ignores the clocks' distributions entirely and picks the next clock to fire uniformly at random among all currently enabled clocks, advancing time by a fixed step dt (default 1.0). It's called "Petri" because a Petri net model always chooses the next event at random.

Because it drives a simulation through improbable event orderings that a distribution-weighted sampler would rarely produce, it is useful for checking that a simulation's enable/disable bookkeeping doesn't crash under unlikely sequences of events. See PetriMethod to select it through the builder interface.

MultipleDirect{SamplerKey,K,Time,Chooser} <: SSA{K,Time}

MultipleDirect is a sampler that contains multiple prefix-search data structures. Here SamplerKey is an identifier for the sampler. K is the type of the Clock key. Time is a Float64 or other clock time. Chooser is a function that selects a prefix-search given a clock key. This differs from a MultiSampler because this works only for exponential distributions, and it combines them by summing propensities.

MultiSampler{SamplerKey,Key,Time}(which_sampler::Chooser) <: SSA{Key,Time}

A sampler returns the soonest event, so we can make a hierarchical sampler that returns the soonest event of the samplers it contains. This is useful because the performance of a sampler depends on the type of the event. For instance, some simulations have a few fast events and a lot of slow ones, so it helps to split them into separate data structures.

The SamplerKey is the type of an identifier for the samplers that this MultiSampler contains. The which_sampler argument is a strategy object that decides which event is sampled by which contained sampler. There is an example of this below.

Once a clock is first enabled, it will always go to the same sampler. This sampler remembers the associations, which could increase memory for simulations with semi-infinite clocks.

Examples

Let's make one sampler for exponential distributions, one for a few clocks we know will be fast and one for slower clocks. We can name them with symbols. The trick is that we need to direct each kind of distribution to the correct sampler. Use a Float64 for time and each clock can be identified with an Int64.

using CompetingClocks
using Distributions: Exponential, UnivariateDistribution

struct ByDistribution <: SamplerChoice{Int64,Symbol} end

function CompetingClocks.choose_sampler(
    chooser::ByDistribution, clock::Int64, distribution::Exponential
    )::Symbol
    return :direct
end
function CompetingClocks.choose_sampler(
    chooser::ByDistribution, clock::Int64, distribution::UnivariateDistribution
    )::Symbol
    if clock < 100
        return :fast
    else
        return :slow
    end
end
sampler = MultiSampler{Symbol,Int64,Float64}(ByDistribution())
sampler[:direct] = OptimizedDirect{Int64,Float64}()
sampler[:fast] = FirstToFire{Int64,Float64}()
sampler[:slow] = FirstToFire{Int64,Float64}()

DebugWatcher()

For debugging, it helps to have visibility into the simulation. This Watcher records everything that is enabled or disabled as a list of all enables and all disabled. It's the complete event history, and you can think of it as the filtration for the process going forward.

watcher = DebugWatcher{String}()
# enable and disable some things.
(watcher.enabled[1].clock,
watcher.enabled[1].distribution,
watcher.enabled[1].te,
watcher.enabled[1].when)

TrackWatcher{K,T}()

This Watcher doesn't sample. It records everything enabled. You can iterate over enabled clocks with a for-loop. If we think of the model as providing changes in which transitions are enabled or disabled, this Watcher accumulates those changes to provide a consistent list of all enabled transitions. Together, a model and this Watcher provide the Semi-Markov core matrix, or the row of it that is currently known.

for entry in tracker
    entry.clock
    entry.distribution
    entry.te
    entry.when
end

TrajectoryWatcher{K,T}

This doesn't sample but calculates the likelihood of the path of samples from start to finish. It has many of the same interface functions as a sampler, but the core value is in the pathloglikelihood function.

PathLikelihoods{K,T}

Calculates the likelihood of one path of events and times according to a vector of different distributions. This accepts a vector of distributions in each enable! call. This can be useful for tuning distribution parameters, but its main use is for importance sampling with mixtures of distribution parameters in order to stabilize importance sampling.