mod integrator

module integrator

Symplectic integrators for HMC dynamics.

The three-stage update mirrors Stan’s base_leapfrog.hpp:17-22: half-kick momentum update, position drift, then the closing half-kick momentum update.

The Omelyan minimum-norm update keeps the same reversible, volume-preserving HMC map while reducing Hamiltonian error for a comparable trajectory length (doi:10.1016/S0010-4655(02)00754-3).

Variables

const OMELYAN_LAMBDA: f64

Omelyan minimum-norm coefficient for the second-order PQPQP update.

Traits

trait HmcIntegrator

Common interface for fixed-step reversible HMC integrators.

Functions

fn evolve<G, M, Obj>(&self, x0: Array1<f64>, p0: Array1<f64>, u0: f64, temp: f64, gradient: &G, momentum: &M, objective: &Obj) -> LeapfrogResult
where
    G: Gradient<f64>,
    M: Momentum + ?Sized,
    Obj: Fn(&Array1<f64>) -> f64

Evolves (x, p) for the configured number of steps at temperature temp. Returns the final state plus the integrator energy error.

Structs and Unions

struct LeapfrogIntegrator

Explicit-leapfrog integrator. Cooling-aware: epsilon_eff = epsilon * sqrt(temp / temp_ref) keeps the linear-stability margin epoch-invariant under SA cooling.

For q-Gaussian momentum the drift dK/dp is computed by the momentum kernel from the current p. The kinetic contribution to the Hamiltonian likewise comes from momentum.kinetic(&p). The integrator stays explicit because dK/dp only depends on p (constant during the drift step).

epsilon: f64

Base step size; rescaled by sqrt(temp / temp_ref) per call.

l_steps: usize

Number of leapfrog steps per evolve call.

temp_ref: f64

Reference temperature used to normalise the cooling rescaling (typically T_0, the initial temperature).

max_delta_h: f64

Divergence threshold; trajectory marked as diverged when abs(delta_h) exceeds max_delta_h. Mirrors Stan base_nuts.hpp line 113 default of 1000.

Implementations

impl LeapfrogIntegrator

Functions

fn new(epsilon: f64, l_steps: usize, temp_ref: f64) -> Self

Constructs a leapfrog integrator. Requires epsilon > 0, l_steps >= 1, and temp_ref > 0.

Traits implemented

impl HmcIntegrator for LeapfrogIntegrator

struct LeapfrogResult

One leapfrog trajectory (L steps) at fixed temperature temp. Returns the final (x, p) and the integrator energy error delta_H = H_new - H_old (small positive on stable trajectories).

x: Array1<f64>

Final position.

p: Array1<f64>

Final momentum.

delta_h: f64

Hamiltonian change: (U_new/T + K_new) - (U_old/T + K_old).

diverged: bool

true if the trajectory diverged (|delta_h| > max_delta_h).

struct OmelyanIntegrator

Omelyan minimum-norm second-order integrator.

One step applies the sequence P(lambda eps), Q(eps/2), P((1 - 2lambda) eps), Q(eps/2), P(lambda eps), where P is a momentum kick and Q is a position drift. It has the same reversible Metropolis correction as leapfrog and is the default HMC trajectory map for BGSA-style runs.

epsilon: f64

Base step size; rescaled by sqrt(temp / temp_ref) per call.

l_steps: usize

Number of Omelyan steps per evolve call.

temp_ref: f64

Reference temperature used to normalise the cooling rescaling.

max_delta_h: f64

Divergence threshold; trajectory marked as diverged when abs(delta_h) exceeds max_delta_h.

lambda: f64

Minimum-norm kick coefficient.

Implementations

impl OmelyanIntegrator

Functions

fn new(epsilon: f64, l_steps: usize, temp_ref: f64) -> Self

Constructs an Omelyan integrator. Requires epsilon > 0, l_steps >= 1, and temp_ref > 0.

Traits implemented

impl HmcIntegrator for OmelyanIntegrator