Samplers

Continuation/Functional Pattern

Key Insight: CommonRandom provides RNG transformation as a higher-order function.

struct CommonRandom{K,R} recorder::Dict{K,Vector{R}} mode::Symbol end

function withcommonrng( crc::CommonRandom{K,R}, clock::K, baserng::R, f::Function ) where {K,R} # Determine which RNG to use if crc.mode == :replaying && haskey(crc.recorder, clock) # ... retrieve saved RNG ... activerng = savedrng else activerng = copy(base_rng) if crc.mode == :recording # ... save RNG snapshot ... end end

  # Call continuation with prepared RNG
  return f(active_rng)

end

function enable!(ctx::SamplingContext, clock, dist, te, when) if ctx.commonrandom !== nothing # Wrap sampler call in RNG transformation withcommonrng(ctx.commonrandom, clock, ctx.rng) do rng enable!(ctx.sampler, clock, dist, te, when, rng) end else enable!(ctx.sampler, clock, dist, te, when, ctx.rng) end

  # Observer watchers use the base RNG (fanned out over the watcher tuple).
  watch_enable!(ctx.watchers, clock, dist, dist, te, when, ctx.rng)

end

Call this before replaying common random numbers. Call it before each replay.

freeze_crn!(ctx::SamplingContext)

Switch to using only previously recorded common random numbers.

reset_crn!(ctx::SamplingContext)

Reset stored CRN draws and time.

split!(dst, src)

Split a src context into multiple copies in dst, adjusting split_weight.

How many times the sampler looked for a random number and found no previous value.

The Pairs of clock keys and values that weren't found.

This updates the survival for a transition in the linear space, according to Gibson and Bruck. Transition was enabled between time record t0 and tn. Divide the survival by the conditional survival between t0 and tn. te can be before t0, at t0, between t0 and tn, or at tn, or after t_n.

$u=\exp\left(-\int_{t_e}^{t_n}\lambda_0(s-t_e)ds\right)\exp\left(-\int_{t_n}^{\tau'}\lambda_{n}(s-t_e)ds\right)$

This updates the survival for a transition in log space, according to Anderson's method.

$\ln u=-\int_{t_e}^{t_n}\lambda_0(s-t_e)ds - \int_{t_n}^{\tau'}\lambda_{n}(s-t_e)ds$

This function decides whether a particular distribution can be sampled faster and more accurately using its cumulative distribution function or using the log of its cumulative distribution function, also called the integrated hazard. The former is used for the Next Reaction method by Gibson and Bruck. The latter is used by the Modified Next Reaction method of Anderson. We are calling the first a linear space and the second a logarithmic space.

steploglikelihood(tw::TrackWatcher, now, when_fires, which_fires)

Calculate the log probability density of a single step in which the which_fires transition fires next. now is the current time. when_fires is the time when which_fires happens so when > now. You have to call this before the transition fires so that it is before transitions are enabled and disabled from the previous step.

One way to compute a marginal likelihood of a particular clock firing $P[K]$ is to integrate:

Using QuadGK
quadgk(t -> exp(steploglikelihood(tw, t0, t, clock)), t0, Inf)[1]

It would be slow but could be done.

steploglikelihood(ctx, when, which)

Step log-likelihood of an event which at when. In a delayed context, which is the public event identity (clock, phase)—the same identity that fire! takes and next_delayed returns.

See sampler.jl for the MultiSampler to understand how we're making a chooser. We would like to implement this: "Compiling Pattern Matching to Good Decision Trees" - Luc Maranget (2008). It looks like MLStyle.jl has a good example of how to do it.