Numerical blocks as choice sites (the spindel side of the spindel ↔ raster
block contract, spindel-raster doc/contract.md).

A block is a fixed-shape group of latent variables with the density factors
they touch, given as a description (data) and capabilities (functions over
primitive arrays):

  (block {:block/id :gauss
          :block/latents [{:name :mu :shape [2] :support :real}]
          :block/target :complete-conditional}
         {:log-density (fn [^doubles theta inputs] lp)
          :value+grad  (fn [^doubles theta inputs] [lp ^doubles grad])})

`(block-dist b inputs)` is the block at its inputs as a distribution whose
density is the block's log target, so a block site is an ordinary choice
site: `(sample (block-dist b inputs) :id :mu :init [0.0 0.0])`. Its value is
θ, the latents flattened in declared order, as a vector of doubles.

The target includes every factor the latents touch: their priors and the
observations that depend on them, which must then not be observed again as
sites of their own. A block without `:sample` cannot be drawn from; start
it at an `:init`. `:sample` need not draw from the target (it cannot know
its normalizer); a block drawn from under importance sampling or SMC, or
proposed from by single-site MH, also needs `:sample-log-density`, the log
density of what `:sample` draws, which weighs the draw:
log target(θ) − log sample-density(θ).

Draft 0 supports unconstrained real latents only; transforms of constrained
ones (`:constrain`, `:unconstrain` with their Jacobians) come with raster's
bijectors.

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.

Conjugate parameters, carried as posteriors instead of sampled.

A static parameter of a conjugate pair need never be drawn: the program
holds its posterior as a value, scores each datum under the posterior
predictive and updates the posterior in closed form —

  (loop [c (conjugate/beta-bernoulli 1 1) [y & more] ys]
    (if y
      (do (observe (conjugate/predictive c) y)
          (recur (conjugate/update c y) more))
      (sample (conjugate/posterior c) :id :p)))   ; the parameter, if needed

— which is exact marginalization (Rao-Blackwellization; the conjugate case
of delayed sampling, Murray et al. 2018). Under SMC the parameter cannot
degenerate, since it is never resampled, and a streaming program pays the
same for it at every step. Every operation is a plain function of plain
values, so it composes with any foerster program and inference method.

Families: `normal-mean` (Normal mean with known observation sd),
`beta-bernoulli`, `gamma-poisson` (Gamma shape and SCALE, as
`foerster.dist/gamma`) and `dirichlet-discrete`.

Compositional probabilistic inference algorithms.

Every method runs a probabilistic program as a savepoint handler: SMC
(`foerster.smc`) for the particle methods — smc-infer,
importance-sampling, pimh-infer, pgibbs-infer, pgas-infer, ipmcmc-infer,
bbvi-infer and kernel-infer with a PInferenceKernel — and replay plus
accept over traces (`foerster.trace`) for the Markov-chain kernels.

Pure inference (`:world-policy :fresh`, the default) runs in fresh worlds;
`:world-policy :fork` in canonical forks of the caller's world (see
`in-canonical-worlds`). Particle measures hold `Sample`s (result, trace,
and a canonical particle's world descriptor).

All functions return Spin<EmpiricalMeasure> for composability;
post-processing is measure-centric (query, predict).

Counterfactual queries by abduction, action and prediction (Pearl).

A program's sample and observe sites are structural equations
(`foerster.mechanism`). Given evidence and interventions:

1. abduction — condition the program on the evidence (importance sampling
   by `gfi/generate`) and read each weighted factual trace's exogenous
   noise off its sites;
2. action — install the interventions;
3. prediction — run the program again in a world of its own with every
   other site fed its factual noise, aligned by address.

The result pairs each factual world with its counterfactual twin, under
the factual weight: a joint measure, so quantities that need both worlds
(probability of necessity, effect of treatment on the treated) are
expectations over it. A site that exists only in the counterfactual world
draws fresh noise and is listed as `:unaligned`: its answer is
interventional, not counterfactual.

Checking an inference result.

Markov chains: `kernel-infer` keeps each chain's draws together (the
measure's `:chain-lengths`), so the draws of a quantity split back into
chains (`chains`) and the convergence diagnostics of Vehtari, Gelman,
Simpson, Carpenter & Bürkner (2021) apply — rank-normalized split R-hat,
bulk and tail effective sample sizes, and the Monte Carlo standard error of
the mean, computed as Stan and ArviZ compute them:

  (diagnostics/summary measure :mu)
  ;; => {:mean … :sd … :quantiles … :rhat 1.002 :ess-bulk 1830.0
  ;;     :ess-tail 1590.0 :mcse 0.004 :chains 4 :draws 4000}

An R-hat above 1.01 or an effective sample size below about 100 per chain
says the chains have not mixed: run longer or change the kernel.

Particle methods: weights carry the information, so `summary` reports the
weighted moments with the weight-based ESS; after resampling that counts
particles, not distinct histories — `distinct-count` of an early site
says how many histories survive.

Model comparison: `pointwise-log-likelihood` gives each draw's log density
of each observation, the input of PSIS-LOO and WAIC; particle methods also
estimate the evidence (`measure/log-marginal`).

Probability distributions, portable between the JVM and JavaScript.

A distribution is a value (a record) with

  (draw d)          a sample, from the current generator
                    (`foerster.random`: a world's stream at a site)
  (logpdf d x)      the log density, or log mass for a discrete law;
                    ##-Inf outside the support
  (cdf d x) (quantile d p) (mean d) (variance d)   where defined

Names and parameterizations follow raster's distributions
(`raster.sci.distributions`, and Distributions.jl): Normal[mu sigma],
Uniform[a b], Exponential[lambda] (rate), Gamma[alpha beta] (shape, SCALE:
mean αβ), Beta[alpha beta], Poisson[lambda]. A model's site laws and a
raster block's compiled densities therefore mean the same thing.

Probabilistic programming effects: unified choose primitive.

This provides the fundamental primitive for compositional probabilistic programming:
- choose: Unified effect for both sampling and observation

A choose site is a savepoint when its world handles `:inference/choose`
(inference: `foerster.smc`, `foerster.trace`); otherwise it is forward
simulation.

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.

Gen's generative function interface over savepoint traces.

Each operation is `spindel.trace/run` or `replay` under a policy of
`foerster.trace`; what this namespace adds is the weight each operation
returns and the discard of `update`, with Gen's identities (Cusumano-Towner
et al. 2019, q the model's internal proposal, here the prior):

  generate(c)       w = log p(t) − log q(t; c)     = Σ log p of constrained and observed sites
  assess(choices)   w = log p(choices)             = the log joint
  update(t, c)      w = log p(t') − log p(t) − log q(fresh choices of t')
  regenerate(t, s)  w = log p(t')/p(t) + log q(t | t')/q(t' | t)

so importance sampling is `generate`, and MH with a selection is
`regenerate` accepted on its weight. A model is a spin whose sample sites
are named with `:id` (or addressed structurally, see `addressing/site-address!`);
constraints are keyed by those addresses.

Every operation returns a CPS operation `(fn [resolve reject])`.
`simulate`, `generate` and `assess` run in a root world and savepoint session
of their own (a session runs one computation, and a world hosts one
session); the worlds live until `close!`.
`update` and `regenerate` replay inside the session of the trace they are
given and leave that trace intact: the caller keeps one of the two and gives
the other back with `spindel.trace/release!`.

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)

Hamiltonian Monte Carlo on block sites (`foerster.block`), within Gibbs.

A block supplies the gradient of its log target; the move is spindel's. The
momentum is drawn from the site's stream, leapfrog integrates with the
block's `:value+grad`, and the endpoint is proposed by replaying the
computation from the block site. It is accepted on the change of the FULL
trace log joint plus the kinetic energy, not on the block's own value:

  log α = [log p(trace') − K(p')] − [log p(trace) − K(p)]

Leapfrog is volume preserving and reversible for any position-dependent
force, so the move is exact whatever the block's target covers. A block
whose target misses factors its latents affect (a downstream observe, a
dependent site) mixes worse, and the step says so: `:incomplete-target?`
is true when the replay changed the log probability of anything outside the
block.

Involutive MCMC over savepoint traces (Neklyudov et al. 2020; Gen's
`involutive_mh`).

A move draws auxiliary variables u ~ q(· ; x) from the current choices x,
maps (x, u) through an involution h to (x', u'), replays the program with
x', and accepts with

  log α = log p(x') − log p(x) + log q(u' ; x') − log q(u ; x) + log |det ∂h/∂(x,u)|

Random-walk, scale, swap and data-driven moves are all such pairs (q, h).
The Jacobian term is the caller's (`:log-jacobian` from the involution);
`fd-log-jacobian` computes it numerically for real-valued maps — raster's AD
is the intended replacement.

Scope: moves that keep the set of sites (fixed dimension). A site the
replay reaches afresh is drawn from its prior and one it no longer reaches
is dropped, both scored as in single-site MH; moves that create or remove
sites through the involution itself (reversible jump) are future work.

Inference kernels: what decides a particle's latent sites.

A `PInferenceKernel` decides a particle's latent sites during execution:
`inference/kernel-infer` runs savepoint SMC whose sample sites take the
value the kernel's `step` gives, which makes importance sampling vs SMC a
choice of kernel, not of engine. The Markov-chain kernels (single-site and
random-walk MH, block Gibbs, HMC) are descriptions that `kernel-infer` runs
as replay plus accept over traces.

Training data from inference: the trajectories SMC drew, with their
rewards and weights, for learning the value estimates (twists) and
proposals that make the next search cheaper.

A steered program (`foerster.steer/model`) records its states under
`[:steer/state t]` and its reward under `:steer/reward`. `trajectories`
reads them back from a measure's particles with their normalized weights;
`draws` resamples them by weight into unweighted draws from the target
p · exp(reward), which is what a model trained on plain examples needs.
Training itself (value heads, proposals over a language model's hidden
states) lives with the models: finetune-rstr's typed-decision records take
a state and the probability of success as a soft target.

  (learn/draws (await (infer/smc-infer (steer/model …) 64)) 256)

Measure abstraction for probabilistic programming.

Measures represent probability distributions over execution traces.
This is the foundation for compositional inference algorithms.

Sample sites as structural equations x = f(u), u exogenous noise.

A distribution is read as a mechanism: `noise` draws u, `push` computes the
value from it, `abduct` goes back from a value to its noise. For an
invertible mechanism (location-scale, inverse CDF of a continuous law)
abduction is exact. For a discrete one, many u give the same value, and
`abduct` draws u from its posterior given the value — here under the
inverse-CDF convention x = F⁻¹(u), u ~ U(0,1), under which that posterior is
uniform on the value's CDF interval. Which convention a discrete site uses
changes counterfactual answers and cannot be told from data, so it is fixed
here and documented rather than guessed.

Parameters come from the distribution at the site, so replaying a site with
its old u under new parents is the counterfactual value.

Memoization and random processes whose state lives in the particle's
world, so it forks with the particle and is discarded with it — Anglican's
`mem` and its Chinese restaurant process, on spindel worlds.

`mem` turns a spin-returning function into one that computes each
argument list once per world: the first call samples, later calls (in the
same particle, and in every particle forked after it) return that value.
Nonparametric models are written with it, as a stick per index or a
mean per cluster:

  (let [mean-of (process/mem (fn [k] (spin (sample (dist/normal 0 10) :id [:mean k]))))]
    (spin … (await (mean-of 3)) …))

Name the sample sites inside a memoized function by its arguments, as
above: a value is drawn once, at the first call.

`crp-draw` seats a customer in a Chinese restaurant process: table k with
probability ∝ its customers, a new table ∝ α. The tables live in the world
under the process's name; the draw is an ordinary sample site over the
existing tables and a new one, so every method — exact enumeration
included — handles it.

Where inference draws its randomness from.

Every draw of a seeded inference run must depend only on WHAT is drawn,
not on when: particles, chains and proposals run concurrently, and one
generator shared in arrival order makes a seeded run reproducible only on
a serial executor. So a draw made while deciding something in a world —
a sample site, the site an MH move selects, its acceptance — reads a
STREAM keyed by that world's seed and what is being decided:

  stream(world, key) = generator seeded by hash(seed(world), key)

World seeds are derived per fork from the parent's seed, the site address
and the fork index (`savepoint/fork`), so a replay's fresh proposal draws
from a fresh stream and the same run draws the same numbers. Draws outside
any world — the seeds of an inference's sessions, SMC's resampling at a
barrier — come from the process generator in program order, which
`set-seed!` seeds.

The generator is xoshiro128** on 32-bit words, so a seed gives the same
numbers on the JVM and in JavaScript.

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.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.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.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.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.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.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.smc`, `:anchors` on the static
parameters).

Probabilistic programs as traced computations.

`sample`, `observe` and `factor` publish savepoints (sites
`:inference/choose` and `:inference/factor`), so a probabilistic program is
run and replayed by `spindel.trace` like any other computation. This
namespace adds what is specific to inference and nothing else: a policy that
scores, pure functions over scored traces, and Metropolis-Hastings as replay
plus an accept step.

Every entry's `:note` carries

  :dist          the site's distribution (nil for a factor)
  :log-prob      log density of the entry's value under :dist, for EVERY
                 site, sampled ones included (a factor's weight)
  :log-proposal  log density under whatever drew the value; absent when the
                 value was not drawn (observed, constrained, kept)
  :observed? :constrained? :kept? :symmetric? :factor?

A `deterministic` site is recorded with `{:deterministic? true}` only; it is
not among `entries`, so it neither scores nor moves.

and the world accumulates the importance weight at `[:inference
:log-weight]`: `log p` of what was observed, constrained or factored, and
`log p - log q` of what a proposal drew. A policy writes it to the world it
is about to resume, so the weight of a fork starts from the weight at its
site.