Black-box variational inference (Ranganath, Gerrish & Blei 2014): a
mean-field q learned by score-function ascent on the ELBO, each iteration
importance sampling from q on savepoint SMC.

The particle cascade (Paige, Wood, Doucet & Teh 2014): asynchronous SMC
without barriers. Each particle runs on its own; at an observation (or a
barrier factor) it compares its weight W with the running mean W̄ of the
weights that have arrived there so far, itself included, and branches on
R = W/W̄ — so no particle ever waits for another (Eq. 14 of the paper):
below the mean it survives with probability R, as one copy weighted W̄;
at or above it takes ⌈R⌉ copies while the children produced at that stage
so far are at most min(K₀, the arrivals before it), ⌊R⌋ otherwise, each
weighted W/M. That feedback keeps the population near K₀; children drawn
independently of each other would let its variance grow without bound.
Copies are forks of the particle's world.

The evidence estimate (1/K₀) Σ W over the particles that reach the end (K₀
launched) is unbiased whatever order the particles arrive in, so slow and
fast particles (an uneven model call, a long simulation) cost no waiting.
`:cap` bounds the particles alive at once: copies beyond it collapse into a
multiplicity carried by the particle (Anglican's pcascade), which counts in
the running means and multiplies its final weight.

Inference ends when every particle has finished — by count, not by wall
clock, which would favour fast particles (Murray, Singh & Lee 2021). The
branching draws come from each particle's world stream, but the running
means depend on arrival order, so runs on a multi-threaded executor are not
reproducible from their seed.

A model whose structure is static given its data, compiled to plain
functions of its latent values, and inference on them without worlds,
savepoints or continuations.

`compile-model` reads the model's instance graph (`analysis/instance`):
every site with a concrete address, a distribution written over constants
and the latent sites before it, and no site whose existence depends on a
value. It returns its latent addresses, a log density of a double array
of their values (the sum of every site's log density, as a run of the
model would score it) and a draw from the prior. `rmh` runs `:rmh`'s
moves (one latent drawn uniformly, a Gaussian step of `:step-size`) on
that density: each move evaluates the density once, as Gen's dynamic MH
re-runs its program once, but as compiled code.

A model that does not qualify is refused with the reasons, so a caller
can fall back to the general path.

Scope: an opt-in (`infer … {:method :rmh :compiled? true}`), on the JVM
only, for models whose latents are all continuous, and single-site
random-walk MH only, with a fixed `:step-size` and the whole density
evaluated per move (no incremental density). Positive and interval
latents step in their own coordinates, a proposal outside the support
being rejected, where the general path's block moves adapt scales in
unconstrained coordinates. The chains are tasks of the run's executor.

Exact inference by enumeration, for programs whose latent choices all have
finite support (`dist/support`): at every latent sample site the particle's
world forks once per value, each branch weighted by that value's
probability, and observations and factors score as usual. Every branch that
ends is one outcome with its exact joint probability, so the measure is the
exact posterior and its `m/log-marginal` the exact evidence.

Branches share the program up to the site where they part — the world fork
is a copy-on-write of everything computed so far — so the cost is the
number of branches, not that number times the program's length. That
number grows exponentially with the latent sites: enumeration is an oracle
for small discrete models, and the reference other methods are checked
against.

  (infer/infer (model) {:method :enumerate})

A latent site with infinite or continuous support is refused
(`::infinite-support`); `:max-branches` (default 100000) bounds the work.

Gradient protocol for variational inference.

Provides grad-log and grad-step for distributions, enabling:
- Black Box Variational Inference (BBVI)
- Amortized inference with learned proposals

Gradients are w.r.t. unconstrained parameters for numerical stability:
- Normal: (mean, log-std)
- Gamma: (log-shape, log-scale); Beta: (log-alpha, log-beta)
- Flip: logit(p)

Point estimates and the Laplace approximation of a block's target.

  (optimize/map-estimate (block/block-dist b inputs) {:init θ0})
  ;; => {:theta θ* :latents {:mu … :sigma …} :log-density lp
  ;;     :iterations 23 :grad-norm 3e-9 :converged? true}

  (optimize/laplace (block/block-dist b inputs) {:init θ0})
  ;; => {:theta θ* :covariance Σ :log-evidence … :draw (fn [] {:mu … …})}

`map-estimate` maximizes the log target by L-BFGS (Nocedal & Wright 2006,
Alg. 7.4, with a backtracking Armijo line search). For a block with
constrained latents it returns the mode in the latents' own coordinates
(σ, p): the transforms' log-Jacobian is taken out of the target, as Stan's
`optimize` does by default; `:jacobian? true` keeps it, giving the mode of
the unconstrained θ.

`laplace` approximates the posterior by a Gaussian in θ at the mode of
the unconstrained target (Jacobian included, as Stan's `laplace`): mean
θ*, covariance the inverse of the negative Hessian, which it takes by
central differences of the block's gradient. Draws are mapped back
through the transforms, so a positive latent stays positive. Its
`:log-evidence` is log p(θ*) + d/2 log 2π − ½ log det(−H), exact for a
Gaussian target. It is cheap and good when the posterior is close to
Gaussian in θ; compare it with NUTS when that is in doubt.

Non-reversible parallel tempering (Syed, Bouchard-Côté, Deligiannidis &
Doucet 2022) on a model that is one block (`foerster.population.theta`).

N + 1 chains sit on the path π_β ∝ q^(1−β)·t^β from the block's draw
density q (the reference) to its target t, at 0 = β_0 < … < β_N = 1. A
scan explores every chain at its own β (the reference chain by an exact
draw, the others by random-walk Metropolis), then offers swaps to
neighbours: even pairs on even scans, odd pairs on odd scans. A swap of
chains i and i+1 is accepted with min(1, exp((β_{i+1} − β_i)·(V(x_i) −
V(x_{i+1})))), V = log t − log q.

The schedule is tuned in rounds of doubling length: the swap rejection
rates r_i give the cumulative barrier Λ(β_i) = Σ_{j<i} r_j, and the next
round spaces its β evenly in Λ. Λ(1) is the global communication barrier.

The last round's states of the target chain are the posterior sample; the
log evidence (q normalized) is the stepping-stone estimate
Σ_i log mean exp((β_{i+1} − β_i)·V) over chain i's states.

Particle Markov chains over an injected sweep operation. World execution is
the sweep's responsibility; normalization, retention and accept/reject live
here, shared by savepoint and canonical-world drivers.

Particle MCMC drivers: PIMH and particle Gibbs (Andrieu, Doucet &
Holenstein 2010), PGAS (Lindsten, Jordan & Schön 2014) and iPMCMC
(Rainforth et al. 2016). In fresh worlds PIMH, particle Gibbs and PGAS
run on savepoint SMC (`foerster.inference.smc`); in canonical worlds, and
iPMCMC always, they run over SMC sweeps (`foerster.inference.particle-mcmc`).

Reading a posterior measure: weighted statistics of the program's value
(`query`), draws passed to a function (`predict`), and posterior
predictive replays of a model (`predictive`).

Sequential core-set Monte Carlo (Beronov, Weilbach, Wood & Campbell, UAI
2021): SMC over data arriving in batches whose rejuvenation moves cost a
bounded amount however much data has arrived.

The target after t batches is p_t(θ) ∝ p₀(θ) · Π_{s≤t} Π_b p(x_b⁽ˢ⁾ | θ).
Plain resample-move SMC (IBIS; Chopin 2002) must score every past datum in
each move, so step t costs O(t). SCMC keeps at most `:memory` M past data
points with weights — a core-set — whose weighted log-likelihood
approximates that of all data so far, and moves the particles against

  log p₀(θ) + Σ_j w_j log p(u_j | θ).

After each step's moves, when the core-set and the new batch exceed M, it
is recompressed: a sparse nonnegative least-squares fit (GIGA; Campbell &
Broderick 2018) of the core-set's log-likelihood functions, centred and
evaluated at the current particles, to their weighted sum. GIGA fixes the
weights' proportions greedily and fits only their common scale; its
weights are then refitted by nonnegative least squares on the points it
selected, the fully corrective update of Campbell & Beronov (NeurIPS
2019) (`:refit :nnls`, the default; `:refit :none` keeps GIGA's weights,
as in the UAI 2021 algorithm). The fit is best where the population has
its mass.

`:projection-inflation k` is a research option and a deviation from the
UAI 2021 algorithm, off by default: the projection is fitted at the
particles spread by the factor k about their mean, μ + k·(θ − μ), rather
than at the particles themselves, so the core-set's log-likelihood is
matched over a wider region than the current population (the posterior
can move outside it as data arrive). It costs N·|E| likelihood
evaluations per compression, and the likelihood must be finite at the
spread points. k = 1, or no factor, is the UAI 2021 projection.

Each new batch reweights the particles by its own exact likelihood, but
the particles it reweights are not distributed as the previous target:
they were moved against the previous step's core-set before that step
compressed it (the core-set plus that step's batch), while the posterior
they stand for from then on uses the compressed core-set. No importance
correction exp(Σ_C' w'·ℓ − Σ_C w·ℓ) is applied for that switch, as in
the paper's Algorithm 1; the next moves target the compressed core-set.
So the particles carry the core-set's approximation error: small when the
log-likelihoods compress well (exchangeable data, few parameters), larger
for a small `:memory`. With M at least the amount of data, SCMC is
resample-move SMC and nothing is approximated.

The engine is `foerster.population.theta/theta-smc` on θ vectors: a step
resamples, then moves every particle by population-scaled random-walk
Metropolis.

Sequential Monte Carlo as a savepoint handler.

The program runs once in a session. The first savepoint it reaches is
forked into N worlds, one per particle; every sample site is decided by
the scoring policy of `foerster.runtime.trace` and recorded in its world's trace;
every observe is scored and then PARKED — its savepoint left pending. When
every particle is parked or has returned, the population is resampled
(when its ESS is below the threshold, folding the log mean weight into the
evidence): the parked savepoints of the chosen ancestors are forked into
the next generation's worlds, the old ones are abandoned, and the new ones
resume. A particle that returned early is carried along with its final
weight. No coordinator, no barrier thread: the barrier is the arrival of
the last particle.

Worlds are forks of savepoints, so particles are frozen, coherent worlds
(see docs/forking.md) and inherit whatever world state the program holds.

Static parameters of state-space programs: particle marginal
Metropolis-Hastings (PMMH; Andrieu, Doucet & Holenstein 2010) and SMC²
(Chopin, Jacob & Papaspiliopoulos 2013).

A program's parameters are named sample sites (`:params`, a set of
addresses). Run with those sites constrained to θ (`foerster.runtime.trace/policy
:constraints`), SMC's evidence estimate is p(θ)·p̂(y | θ) — a constrained
site adds its prior density to the weight — so the posterior ratio of two
parameter values is the difference of two SMC log-evidences.

- `pmmh` is a Markov chain over θ whose every proposal runs an SMC over
  the program's other sites: exact for any number of particles.
- `smc2` is online: θ-particles each carry an inner streaming SMC
  (`foerster.inference.smc/stream`); a pushed observation reweights each by its
  inner evidence increment; when the θ-particles' ESS falls they are
  resampled — a duplicate forks its inner population's worlds, so copies
  evolve independently — and moved by a PMMH step that runs a fresh inner
  filter on the data so far.

Priors for θ come from the program itself: a θ-particle starts from a
simulation of the program (`foerster.runtime.gfi/simulate`). Correlated PMMH, which
needs SMC driven by explicit noise, is not provided.

Steering a sequential process by SMC: a program that takes steps — an
agent's turns, a model's sentences — whose proposals come from the process
itself (a language model, a simulator: randomness without a sample site,
so its density cancels as a prior draw's does), scored by a value estimate
at each step and by a reward at the end.

The target is p(trajectory) · exp(reward): the process's own law tilted by
the reward. A value estimate ψ (a verifier, a process reward model, a
judge) twists the intermediate targets without changing the final one
(twisted SMC; Whiteley & Lee 2014, Zhao et al. 2024): step t adds
log ψ_t − log ψ_{t−1} as a barrier factor, so SMC resamples on it, and the
end adds reward − log ψ_T. With no value estimate every intermediate
factor is 0 and SMC is best-of-N weighted by the reward.

  (infer/smc-infer (steer/model {:init s0
                                 :step (fn [state] (spin …))   ; the next state
                                 :value (fn [state] (spin …))  ; log ψ
                                 :reward (fn [state] (spin …)) ; log potential
                                 :done? (fn [state] …)})
                   16 {:resampling :stratified})

`:step`, `:value` and `:reward` may return a spin or a plain value.

A step drawn from a proposal q other than the process p — a learned
proposal, a smaller model — returns `(weighted state log-w)` with
log-w = log p(state | previous) − log q(state | previous); the weight joins
that step's factor and the target stays p · exp(reward).

Tempered SMC: an SMC sampler (Del Moral, Doucet & Jasra 2006) over whole
program traces, from the prior to the posterior through the targets

  π_β(x) ∝ p(x) · L(x)^β,   0 = β_0 < β_1 < … < β_T = 1,

L the product of the program's observations and factors. Every particle is
a complete run of the program (`foerster.runtime.gfi`), recorded at temperature 0
with each observation's untempered log-likelihood in its note. A step picks
the next β adaptively, so that the conditional ESS of the incremental
weights L^(β'−β) is `:ess-target`·N (Zhou, Johansen & Aston 2016), reweights,
resamples and moves every particle by Metropolis-Hastings targeting π_β'
(`foerster.runtime.trace/mh-step` at that temperature; random-walk proposals scaled
by the population's spread on continuous sites, prior proposals on the
others). The evidence estimate is the product of the steps' mean
incremental weights.

With `:waste-free P` (Dau & Chopin 2022) a step resamples N/P particles
and keeps every state of their P-step chains as particles.

For data that arrive one observation at a time (IBIS), use SMC with
resample-move instead (`foerster.inference.smc`, `:anchors` on the static
parameters).

The worlds a particle method runs in, and the options every driver
checks. `:world-policy :fresh` (the default) runs pure inference on
savepoints in fresh worlds; `:fork` runs it in canonical forks of the
caller's world (`in-canonical-worlds`), from one frozen fork of the caller
when a method runs several times (`foerster.runtime.canonical`).

Errors and option keys keep `foerster.core`'s keywords: they are what the
public entry points raise and take.