C++ API

Sampling functionality

Top-level call

This function will spawn threads to perform end-to-end sampling

template<std::uniform_random_bit_generator RNG, ChainHandler H, GlobalHandler GH, InterruptCallback IC, LogpGrad F>
inline void walnutpie::walnuts(std::size_t seed, std::vector<H> &chain_handlers, GH &global_handler, const IC &interrupt_callback, const F &log_p_grad, const WalnutsConfig &config)

Return the chain records from running Walnuts with the specified seed, sampling event handlers, and configuration.

Template Parameters:

Handler – The type of the event handlers.

Parameters:

• seed – [in] The seed for the pseudo-random number generator.

• chain_handlers – [in] The collection of chain-specific handlers, which are called back.

• global_handler – [in] The handler for global cross-chain events.

• interrupt_callback – [in] The callback for stopping.

• log_p_grad – [in] The log density and gradient function, called back.

• config – [in] The configuration for Walnutpie.

Throws:

std::invalid_argument – If the number of handlers doesn’t match the initialization configuration’s number of chains.

Iterator-style samplers

These classes implement Walnuts in an iteration-per-call style iterator.

template<LogpGrad F, std::uniform_random_bit_generator RNG, SampleHandler H>
class WalnutsSampler

The Walnuts Markov chain Monte Carlo (MCMC) sampler.

The sampler is constructed with a base random number generator, a log density and gradient function, an initialization, and several tuning parameters. It provides a no-argument functor for generating the next element of the Markov chain.

Template Parameters:

• F – The type of the log density and gradient function.

• RNG – The type of the base random number generator.

• Handler – The type of the sampling event handler.

Public Functions

inline WalnutsSampler(RNG &rng, H &sample_handler, const F &logp_grad, const Eigen::VectorXd &theta, const Eigen::VectorXd &inv_mass, double macro_time, std::size_t max_nuts_depth, std::size_t max_step_halvings, std::size_t min_micro_steps, double max_error)

Construct a Walnuts sampler from the specified RNG, target log density/gradient initialization, and tuning parameters.

Parameters:

• rng – [inout] The base random number generator.

• sample_handler – [inout] The sampling and on-stop event handler.

• logp_grad – [in] The target log density and gradient function (see the class documentation.

• theta – [in] The initial position.

• inv_mass – [in] The diagonal of the diagonal inverse mass matrix.

• macro_time – [in] The macro time discretization interval.

• max_nuts_depth – [in] The maximum number of trajectory doublings for Nuts.

• max_step_halvings – [in] The maximum number of times the step size is halved.

• min_micro_steps – [in] The minimum number of micro steps per macro step.

• max_error – [in] The log of the maximum error in joint densities allowed in Hamiltonian trajectories.

Throws:

• std::invalid_argument – If inv_mass_matrix has non-positive or infinite entries.

• std::invalid_argument – If macro_time is not positive or not finite.

• std::invalid_argument – If max_nuts_depth is not positive.

• std::invalid_argument – If max_step_halvings is not positive.

• std::invalid_argument – If min_micro_steps is not positive.

• std::invalid_argument – If max_error is not positive or not finite.

WalnutsSampler(const WalnutsSampler &sampler) = default

Construct a sampler by copying the specified sampler.

Parameters:

sampler – [in] Sampler to copy.

WalnutsSampler(WalnutsSampler &&sampler) = default

Construct a sampler by moving the specified sampler.

Parameters:

sampler – [in] Sampler to move.

inline double operator()()

Generate the next draw and send it to the handler and return its log density.

Returns:

The unnormalized log density of the next draw.

inline const Eigen::VectorXd &inverse_mass_matrix_diagonal() const noexcept

Return a constant reference the diagonal of the diagonal inverse mass matrix.

The value of the inverse mass matrix will change on subsequent calls to operator()(S&).

Returns:

The diagonal of the inverse mass matrix.

inline double macro_time() const noexcept

Return the macro time discretization interval for Nuts.

Returns:

The time discretization interval for Nuts.

inline double max_error() const noexcept

Return the maximum error allowed among Hamiltonians.

Returns:

The maximum error allowed among Hamiltonians.

inline std::size_t dim() const noexcept

Return the number of dimensions.

Returns:

The number of dimensions.

template<LogpGrad F, std::uniform_random_bit_generator RNG, ChainHandler H>
class AdaptiveWalnuts

The adaptive Walnuts sampler.

The adaptive Walnuts sampler is configured in the constructor, then provides a functor method operator()() for returning the next state in warmup. Warmup re-estimates step size and mass matrix each iteration, exponentially discounting the past.

Template Parameters:

• F – Type of log density/gradient function.

• RNG – Type of base random number generator.

• Handler – Type of adaptation and sampling event handler.

Public Functions

inline AdaptiveWalnuts(RNG &rng, H &handler, const F &logp_grad, const InitChainConfig &init_chain_cfg, const WarmupConfig &warmup_cfg, const SamplingConfig &sampling_cfg)

Construct an adaptive Walnuts sampler.

The configuration objects, the base random number generator, and the log density/gradient function are held by reference. The RNG changes state every time a random number is generated. The target depth specifies the expected Nuts tree depth, which is controlled through the minimum number of micro steps per macro step and adjusted with a mean estimator to achieve this average.

Parameters:

• rng – [in] The base random number generator, stored by reference and modifed.

• handler – [inout] Event handler for adaptation and sampling, stored by reference and called back.

• logp_grad – [in] The target log density and gradient function.

• init_chain_cfg – [in] The initialization configuration for a single chain.

• warmup_cfg – [in] The warmup configuration.

• sampling_cfg – [in] The sampling configuration.

inline void operator()()

Generate the next state for adaptation and the handler.

This method should be called a number of time equal to the number of warmup iterations desired. These warmup draws are not drawn from a Markov chain and are not valid for inference. After warmup, call sampler() to return a sampler that fixes the tuning parameters and provides a proper Markov chain.

inline WalnutsSampler<F, RNG, H> sampler()

Return a Walnuts sampler with the current tuning parameter estimates.

The returned sampler forms a proper Markov chain. The method passes along the compound random number generator and log density function and is hence not marked const.

Returns:

The Walnuts sampler with current tuning parameter estimates.

inline Eigen::VectorXd inv_mass() const

Return the diagonal of the diagonal inverse mass matrix.

Returns:

The diagonal of the inverse mass matrix.

inline double step_size() const

Return the step size.

Returns:

The step size.

inline std::size_t min_micro_steps() const

Return the minimum number of micro steps per macro step.

Returns:

The minimum number of micro steps per macro step.

inline std::size_t dim() const noexcept

Return the number of dimensions of the position.

Returns:

The number of dimensions.

inline double log_step_size() const noexcept

Return the natural logarithm of the step size.

Returns:

The log of the step size.

inline Eigen::VectorXd log_mass() const

Return the natural logarithm of the diagonal of the diagonal mass matirx.

Returns:

The log of the diagonal of the mass matrix.

inline std::size_t iter() const noexcept

Return the current iteration.

Returns:

The iteration.

Configuration

The following classes (and their builders) are used to configure Walnuts.

class WalnutsConfig

Encapsulated configuration for Walnutpie.

Walnuts configurations include initialization, warmup, and sampling configurations.

Public Functions

inline WalnutsConfig(InitConfig init, WarmupConfig warmup, SamplingConfig sampling)

Construct a Walnuts configuration given the component configurations.

The arguments will be moved if they are rvalues and copied if lvalues.

Parameters:

• init – [in] The initialization configuration.

• warmup – [in] The warmup configuration.

• sampling – [in] The sampling configuration.

inline const InitConfig &init() const noexcept

Return the initialization configuration.

Returns:

The initialization configuration.

inline const WarmupConfig &warmup() const noexcept

Return the warmup configuration.

Returns:

The warmup configuration.

inline const SamplingConfig &sampling() const noexcept

Return the sampling configuration.

Returns:

The sampling configuration.

class InitConfigBuilder

The builder for initialization configurations.

The usage to return an InitConfig is InitConfigBuilder(4, 20).step_sizes(0.5).build(); with any number of config methods chained between the construction and call to build.

Public Functions

inline InitConfigBuilder(std::size_t num_chains, std::size_t dims)

Construct an initialization builder of the given sizes.

Parameters:

• num_chains – [in] The number of Markov chains.

• dims – [in] The dimensionality of each chain.

inline InitConfigBuilder &step_sizes(double v)

Set the step sizes to all be the specified value.

Parameters:

v – [in] The step size.

Throws:

std::invalid_argument – If the step size is not finite and positive.

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &step_sizes(const std::vector<double> &v)

Set the step sizes to all be the specified values.

Parameters:

v – [in] The step sizes.

Throws:

• std::invalid_argument – If any of the step sizes are not finite positive.

• std::invalid_argument – If the number of chains doesn’t match the number specified in the constructor.

Returns:

A reference to this builder for chaining.

template<std::uniform_random_bit_generator RNG>
inline InitConfigBuilder &positions(RNG &rng, double init_scale)

Randomly initialization the positions.

Initialization is independent in each dimension with values drawn from a zero-centered normal distribution with the specified scale.

Template Parameters:

RNG – The type of the base random number generator.

Parameters:

• rng – [inout] The base random number generator.

• init_scale – [in] The scale of the normal initial values.

Throws:

std::invalid_argument – If the initial scale is not finite and positive.

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &positions(const Eigen::VectorXd &v)

Initialize the positions all to the same value.

Parameters:

v – [in] The initial position.

Throws:

• std::invalid_argument – If the dimensionality doesn’t match that specified during construction.

• std::invalid_argument – If any of the initial positions contains non-finite values.

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &positions(const std::vector<Eigen::VectorXd> &vs)

Initialize the positions to the specified values.

Parameters:

vs – [in] The initial positions.

Throws:

• std::invalid_argument – If the number of initial positions doesn’t match the number of chains specified in the constructor.

• std::invalid_argument – If any of the initial positions contains non-finite values.

• std::invalid_argumet – If any of the initial positions has a dimensionality that does not match the dimensionality specified in the constructor.

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &positions(std::vector<Eigen::VectorXd> &&vs)

Initialize the positions to the specified values via move.

Parameters:

vs – [in] The initial positions.

Throws:

• std::invalid_argument – If the number of initial positions doesn’t match the number of chains specified in the constructor.

• std::invalid_argument – If any of the initial positions contains non-finite values.

• std::invalid_argumet – If any of the initial positions has a dimensionality that does not match the dimensionality specified in the constructor.

Returns:

A reference to this builder for chaining.

template<LogpGrad F>
inline InitConfigBuilder &masses(const F &logp_grad, double mass_smoothing, bool average_masses = false)

Initialize the masses using the Nutpie outer product strategy.

Following Nutpie, the initialization uses a smoothed negative outer product of gradient, which is the absolute value of the outer product of gradients linearly interpolated with a unit matrix with weight mass_smoothing on the unit matrix and 1 - mass_smoothing on the regularized outer product.

If the flag average_masses is true, then each chain’s mass matrix is set to the geometric average of the per-chain mass matrixes.

See: Seyboldt, Adrian and Carlson, Eliot and Carpenter, Bob. 2026. [Preconditioning Hamiltonian Monte Carlo by minimizing Fisher divergence](https://arxiv.org/abs/2603.18845v1). arXiv 2603.18845.

Template Parameters:

LPG – The type of the log density and gradient function.

Parameters:

• logp_grad – [in] The log density and gradient function, called back.

• mass_smoothing – [in] The additive smoothing for mass matrices.

• average_masses – [in] Set to true to geometrically average mass matrices.

Throws:

std::invalid_argumet – If the mass smoothing is not in (0, 1).

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &masses(const Eigen::VectorXd &v)

Initialize the mass matrices all to the same value.

Parameters:

v – [in] The initial diagonal mass matrix.

Throws:

• std::invalid_argument – If the dimensionality doesn’t match that specified during construction.

• std::invalid_argument – If any of the initial mass matrix diagonals contains non-finite or non-positive values.

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &masses(const std::vector<Eigen::VectorXd> &vs)

Initialize the mass matrices to the specified values.

Parameters:

vs – [in] The initial mass matrices.

Throws:

• std::invalid_argument – If the number of initial mass matrices doesn’t match the number of chains specified in the constructor.

• std::invalid_argument – If any of the initial mass matrices contains non-finite values.

• std::invalid_argumet – If any of the initial mass matrices has a dimensionality that does not match the dimensionality specified in the constructor.

Returns:

A reference to this builder for chaining.

inline InitConfigBuilder &masses(std::vector<Eigen::VectorXd> &&vs)

Initialize the mass matrices to the specified values via move.

Parameters:

vs – [in] The initial mass matrices.

Throws:

• std::invalid_argument – If the number of initial mass matrices doesn’t match the number of chains specified in the constructor.

• std::invalid_argument – If any of the initial mass matrices contains non-finite values.

• std::invalid_argumet – If any of the initial mass matrices has a dimensionality that does not match the dimensionality specified in the constructor.

Returns:

A reference to this builder for chaining.

inline InitConfig build()

Return the initialization configuration.

Returns:

The initialization configuration.

template<std::uniform_random_bit_generator RNG, LogpGrad F>
inline InitConfig adapt_step_build(RNG &rng, const F &logp_grad)

Heuristically adapt the initial step sizes, then return the initialization configuration.

Template Parameters:

• RNG – Type of the base random number generator.

• F – Type of the log density and gradient function.

Parameters:

• rng – [in] The base random number generator.

• logp_grad – [in] The log density and gradient function.

class WarmupConfigBuilder

The builder for WarmupConfig objects.

Public Functions

inline WarmupConfigBuilder &min_max_iter(std::size_t min_iter, std::size_t max_iter)

Set the minimum and maximum number of warmup iterations.

Parameters:

• min_iter – [in] The minimum number of warmup iterations.

• max_iter – [in] The maximum number of warmup iterations.

Throws:

std::invalid_argument – If min_iter > max_iter.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_size_converge_tol(double v)

Set the step size convergence tolerance.

Parameters:

v – [in] The step size convergence tolerance.

Throws:

std::invalid_argument – If the tolerance is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &mass_converge_tol(double v)

Set the mass matrix L2-norm convergence tolerance.

Parameters:

v – [in] The mass matrix L2-norm convergence tolerance.

Throws:

std::invalid_argument – If the tolerance is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &mass_init_count(double v)

Set the mass matrix estimator initial count.

Parameters:

v – [in] The mass matrix estimator initial count.

Throws:

std::invalid_argument – If the initial count is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &mass_additive_smoothing(double v)

Set the mass matrix estimator additive smoothing.

Parameters:

v – [in] The mass matrix estimator additive smoothing.

Throws:

std::invalid_argument – If the smoothing is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &max_macro_steps_target(double v)

Set the target number of macro steps.

Parameters:

v – [in] The target number of macro steps.

Throws:

std::invalid_argument – If the target is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_accept_rate_target(double v)

Set the accept-rate target for step-size estimation.

Parameters:

v – [in] The accept-rate target.

Throws:

std::invalid_argument – If the accept rate is not in (0, 1).

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_learning_rate(double v)

Set the step size learning rate.

Parameters:

v – [in] The step size learning rate.

Throws:

std::invalid_argument – If the learning rate is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_gradient_decay(double v)

Set the gradient decay for the step-size estimator.

Parameters:

v – [in] The gradient decay for the step-size estimator.

Throws:

std::invalid_argument – If the decay is not in (0, 1).

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_sq_gradient_decay(double v)

Set the squared gradient decay for the step-size estimator.

Parameters:

v – [in] The squared gradient decay for the step-size estimator.

Throws:

std::invalid_argument – If the decay is not in (0, 1).

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_stabilization(double v)

Set the step-size estimator stabilization term.

Parameters:

v – [in] The step-size estimator stabilization term.

Throws:

std::invalid_argument – If the stabilization is not finite and positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &step_learn_rate_decay(double v)

Set the learning rate decay exponent for step size.

Parameters:

v – [in] The learning rate decay exponent.

Throws:

std::invalid_argument – If the decay exponent is not in (0, 1).

Returns:

This builder for chaining.

inline WarmupConfigBuilder &publish_stride(std::size_t v)

Set the stride between publishing statistics for convergence monitoring.

Parameters:

v – [in] The stride for publishing statistics for convergence monitoring.

Throws:

std::invalid_argument – If the stride is not positive.

Returns:

This builder for chaining.

inline WarmupConfigBuilder &yield_period(std::size_t v)

Set the iteration period between chain threads yielding.

Parameters:

v – [in] The iteration period between chain threads yielding.

Throws:

std::invalid_argument – If the yield period is not positive.

Returns:

This builder for chaining.

inline WarmupConfig build()

Return the warmup configuration.

Returns:

The warmup configuration.

class SamplingConfigBuilder

The builder for sampling configurations.

An example use would be:

SampleConfigBuilder(50u, 100u)
    .max_step_halvings(4u)
    .min_micro_steps(2u)
    .build();

Public Functions

inline SamplingConfigBuilder &min_max_iter(std::size_t min_iter, std::size_t max_iter)

Set the minimum and maximum number of iterations.

Parameters:

• min_iter – [in] The minimum number of iterations.

• max_iter – [in] The maximum number of iterations.

Throws:

std::invalid_argument – If the minimum number of iterations is greater than the maximum number of iterations.

Returns:

A reference to this builder for chaining.

inline SamplingConfigBuilder &max_trajectory_doublings(std::size_t v) noexcept

Set the maximum number of trajectory doublings.

Parameters:

v – [in] The maximum number of trajectory doublings.

Returns:

A reference to this builder for chaining.

inline SamplingConfigBuilder &max_step_halvings(std::size_t v) noexcept

Set the maximum number of step size halvings.

Parameters:

v – [in] The maximum number of step size halvings.

Returns:

A reference to this builder for chaining.

inline SamplingConfigBuilder &max_hamiltonian_error(double v)

Set the maximum error in the Hamiltonian for Walnutpie.

Parameters:

v – [in] The maximum error in the Hamiltonian for Walnutpie.

Throws:

std::invalid_argument – If the error is not finite and positive.

Returns:

A reference to this builder for chaining.

inline SamplingConfigBuilder &min_micro_steps(std::size_t v)

Set the minimum number of micro steps per macro step.

Parameters:

v – [in] The minimum number of micro steps per macro step.

Throws:

std::invalid_argument – If the minimum number of steps is not positive.

Returns:

A reference to this builder for chaining.

inline SamplingConfigBuilder &rhat_converge_tol(double v)

Set the R-hat convergence tolerance.

Parameters:

v – [in] The R-hat convergence tolerance.

Throws:

std::invalid_argument – If the tolerance is not finite and > 1.

Returns:

A reference to this builder for chaining.

inline SamplingConfig build()

Return the sampling configuration.

Returns:

The sampling configuration.

class InitConfig

The initialization configuration for multiple Markov chains.

Rather than a public constructor, it is built using an InitConfigBuilder instance.

The initialization configuration specifies a step size, initial position, and initial mass matrix.

class WarmupConfig

The warmup configuration object. The object supplies methods for all of the tuning parameters for warmup.

class SamplingConfig

A class to hold the configuration for the Walnuts sampler.

Concepts

The following concepts describe the types expected by walnutpie.

template<typename F>
concept LogpGrad

Concept for a log density and gradient function.

A type F satisfies LogpGrad if an object of type const F& can be called with arguments (const Eigen::VectorXd&, double&, Eigen::VectorXd&) and the call returns void. The first argument is the position at which to evaluate, and the second and third are output parameters set to the log density and its gradient, respectively.

tparam F:

The callable type to constrain.

template<typename H>
concept ErrorCallback

Concept for handler for errors The following member is required.

• on_logp_exception(const Eigen::VectorXd&, std::exception&) called when the log density function throws an error

template<typename H>
concept SampleHandler

Concept for a handler of sampling events.

A type H satisfies SampleHandler if it satisfies ErrorCallback and additional provides the following member functions, each callable on a non-const instance and returning void:

• on_sample(const Eigen::VectorXd&, double) called once per draw with the position and log density.

template<typename C>
concept ChainHandler

An extension of the SampleHandler concept for additionally handling warmup events.

A type C satisfies ChainHandler if it provides the following member functions, each callable on a non-const instance and returning void:

• on_warmup(const Eigen::VectorXd&, double, double, const Eigen::VectorXd&) called once per warmup draw with the position, log density, step size, and diagonal inverse mass matrix.

• on_warmup_complete(double, const Eigen::VectorXd&) called once when warmup finishes, with the final step size and diagonal inverse mass matrix.

• on_sample(const Eigen::VectorXd&, double) called once per post-warmup draw with the position and log density.

template<typename H>
concept GlobalHandler

Concept for a handler of cross-chain events.

This only handles R-hat updates now.

A type H satisfies Handler if it provides:

• on_r_hat(double) callable on a non-const instance, returning void,

template<typename H>
concept InterruptCallback

Concept for an interrupt callback.

A type H satisfies Handler if it provides:

• received_interrupt() will return true if the process should be interrupted.