API

initialize(rng, logdensity, initial; move=StretchMove(), executor=SerialExecutor(),
           walker_ids=1:nwalkers, logdensities=nothing)

Create one ensemble from a vector of finite coordinate vectors or a matrix with axes (coordinate, walker). The coordinates must have full affine rank and initial log densities other than NaN or +Inf. The state owns copies of coordinates and RNG. Supported RNGs are Random123 Philox4x/Threefry4x with UInt64 counters. Use independent RNG seeds or streams for independent ensembles.

Supply logdensities to reuse cached values without evaluating the target. The vector must match walker order and contain real values other than NaN or +Inf. Values are copied and promoted independently of coordinate precision. The caller must ensure they equal the target at the initial coordinates.

Walker IDs default to 1:nwalkers and remain fixed throughout the run. The density must not mutate its input. A scalar density must support concurrent calls with ThreadedExecutor; a batched callback controls its own parallelism. Do not mutate state fields. Use current_state for borrowed inspection and snapshot for copies.

step!(state)
step!(state, nsweeps)

Advance one full ensemble sweep, or nsweeps complete sweeps without storing a history. Mutate and return the state. Zero sweeps leave a valid state unchanged. Each sweep selects one move and updates every walker once, in ordered groups against frozen complements. Use current_state to consume each sweep without a copy, or snapshot to retain it. An exception during a sweep invalidates the state. Initialize a fresh state after correcting the target instead of resuming a partially completed sweep.

step!(state, rng; proposal_index=1)

Advance one sweep using an externally addressed RNG without changing rng. The RNG must match the initialization RNG type and leave two partition levels for purposes and walkers. proposal_index selects disjoint proposal streams. Inner mixtures use purpose 2 for selection, leaving purpose 1 for an outer mixture. Cyclic mixtures count calls to this state. The supplied address replaces the standalone cycle/step address for this sweep.

sample!(state, nsweeps)

Advance and collect complete sweeps. Positions have axes (coordinate, walker, sweep). Rejections repeat the current state. Returned arrays own their storage. Each state represents one coupled ensemble, not independent walker chains.

Return a named tuple with positions, logdensities and accepted (the latter two indexed by walker and sweep), move_indices (one move index per sweep), and walker_ids (one ID per walker). The initial positions are not included.

current_state(state)

Inspect a valid ensemble without copying arrays or collecting a history. Return a named tuple with borrowed positions (a vector of coordinate vectors), logdensities, accepted and walker_ids (one entry per walker), and attempts and acceptances (cumulative counts per move). Scalar move_index names the move selected for the latest sweep; sweep_count counts completed sweeps. Before the first sweep, both scalars and all counts are zero, and accepted is false for every walker. candidates, candidate_logdensities, and acceptance_probabilities describe the last attempted transition, including rejections. Before stepping or after synchronization they describe the current points with zero acceptance probabilities and move_index == 0.

Treat all borrowed arrays as read-only. Consume them before the next mutation of the state. Scalar metadata is captured at query time, not updated live. Do not inspect a state concurrently with stepping. Use snapshot to retain an owned copy. A state invalidated by a failed sweep cannot be queried.

snapshot(state)

Copy the fields returned by current_state, including every nested array. Later steps and edits to the snapshot cannot affect each other. A snapshot contains sample data and counts, not a resumable sampler checkpoint.

synchronize!(state, positions, logdensities)

Copy externally owned positions and their matching cached log densities into a valid state. Inputs must not alias arrays borrowed from this state. Accept a coordinate-by-walker matrix or a vector of coordinate vectors. Walker identities, RNG, sweep count, mixture phase, and cumulative counters stay unchanged. Clear the last-transition metadata. No target calls or affine-rank check run here. The caller must preserve density consistency and affine rank. Use validate_positions after replacing positions through initialization retries.

validate_positions(positions)

Check that a vector of coordinate vectors or a coordinate-by-walker matrix has positive, matching dimensions, finite real coordinates, and full affine rank. Use the floating-point coordinate type for the rank tolerance. Return nothing without changing the input. Rank validation uses an owned host matrix. Call this after retry initialization, before synchronizing replacement positions. Move-specific walker requirements are checked by initialize.

StretchMove(; scale=2)

Stretch a walker away from a companion in the frozen complement. The stretch factor lies in [1/scale, scale]; scale must be finite and greater than one. Uses two groups and requires at least 2d walkers for dimension d.

DEMove(; gamma0=nothing, sigma=1e-5)

Add a scaled difference of two distinct companions in the frozen complement. The scale is gamma0 * (1 + sigma * randn()), with default gamma0 = 2.38 / sqrt(2d) for dimension d. gamma0 must be finite and positive when supplied. sigma must be finite and nonnegative. Uses two groups and requires at least max(2d, 4) walkers.

DESnookerMove(; scale=1.7)

Project a companion difference onto the line from a reference walker to the active walker, then multiply that displacement by scale. The scale must be finite and positive. Degenerate reference directions produce a rejection. Uses four groups and requires at least max(2d, 4) walkers for dimension d. This move is not generally affine-equivariant.

GaussianReplacementMove(; shrinkage=0.5)

Replace a walker with an independent Gaussian draw fitted to the frozen complement, with the independence-proposal Hastings correction. The unbiased covariance is shrunk toward tr(C)/d * I by finite shrinkage in [0, 1]. Uses two groups and requires at least max(2d, 4) walkers, or 2(d+1) when shrinkage is zero. A non-positive or non-finite fit rejects the whole group.

MoveMixture(moves, weights; schedule=:random)

Select one move per complete sweep. schedule=:random gives fixed random selection probabilities after normalization, regardless of weight types. schedule=:cycle gives a repeating weighted cycle in component order and requires integer-valued weights. Weights must be finite and nonnegative, with at least one positive weight. The schedule and weights do not adapt.

BatchedLogDensity(scalar, batch!)

Wrap scalar and batched log-density callbacks. Initialization uses scalar(x) unless cached logdensities are supplied. During sampling, batch!(values, positions) receives borrowed SubArray views of output storage and a coordinate-by-candidate matrix containing only valid proposals. The callback must treat positions as read-only, fill every entry of values with the same result as scalar, and must not retain either array. It is called once per nonempty frozen proposal group and controls any parallelism used to evaluate that group.

SerialExecutor()

Evaluate walker proposals serially within each group of a sweep.

ThreadedExecutor()

Use Julia threads when proposal groups contain enough work to benefit. The first few groups calibrate task size and compare serial and threaded costs for each move. Calibration uses normal proposals, without extra density calls. Groups still advance in order. The log-density function must support concurrent calls. A deterministic target gives the same trajectory as SerialExecutor for the same initial state, move, RNG, and walker IDs.

KernelExecutor()

Execute proposals and acceptance with KernelAbstractions on the initial matrix's backend. Requires BatchedLogDensity. Initialization evaluates the scalar callback on host vectors; sampling calls the batch callback with backend arrays. The caller places target data explicitly. CPU and CUDA arrays are supported.