PDMPFlux

ForwardECMCAD(dim, U; kwargs...)

Create ForwardECMC sampler with automatic differentiation.

Arguments

• dim::Int: Dimension of the state space
• U::Function: Potential function

Keywords

• grid_size::Int=10: Number of grid points for upper bound
• tmax::Union{Float64, Int}=2.0: Maximum time horizon
• signed_bound::Bool=true: Use signed bound strategy
• adaptive::Bool=true: Use adaptive time horizon
• AD_backend::String="ForwardDiff": Automatic differentiation backend
• ran_p::Bool=true: Use random orthogonal refresh
• mix_p::Float64=0.5: Mixture probability for refreshment

Returns

• ForwardECMC: Configured sampler instance

RHMCAD(dim::Int, U::Function; kwargs...)

Convenience constructor that builds ∇U via the selected AD backend, like ZigZagAD etc.

A buffer function under development
supports faded color animation

Create a contour plot of the potential function U (2D).

Arguments

• U: potential function taking a 2D vector
• x_range: range for the x-axis (default: range(-5, 5, length=100))
• y_range: range for the y-axis (default: range(-5, 5, length=100))
• levels: number of contour levels (default: 20)
• color: color palette (default: :viridis)
• fill: whether to fill contours (default: true)
• show_grid: show grid (default: false)
• show_axis: show axis ticks/labels (default: false)
• show_title: show title (default: false)
• xlabel: x-axis label (default: L"x_1")
• ylabel: y-axis label (default: L"x_2")
• background: background (default: :transparent)
• linewidth: contour line width (default: 1)
• filename: output filename to save (default: nothing)

sample(): draw samples from a `AbstractPDMP` sampler.

This is a thin wrapper that calls `sample_skeleton()` and then `sample_from_skeleton()`.

Args:
    N_sk (Int): Number of skeleton points to generate.
    N_samples (Int): Number of final samples to generate from the skeleton.
    xinit (Array{Float64, 1}): Initial position.
    vinit (Array{Float64, 1}): Initial velocity.
    seed (Int): Seed for random number generation.
    verbose (Bool, optional): Whether to print progress information. Defaults to true.

Returns:
    Array{Float64, 2}: Array of samples generated from the PDMP model.

Sample along a skeleton and return a `Matrix{Float64}` whose rows store per-dimension time series.

Args:
    dt (Float64): The time step.
    output (PdmpOutput): The PDMP output containing the trajectory information.

Returns:
    Array{Float64, 2}: The sampled points from the PDMP trajectory skeleton.

Sample along a skeleton and return a Matrix{Float64} whose rows store per-dimension time series.

Args: N (Int): The number of samples to generate. output (PdmpOutput): The PDMP output containing the trajectory information.

Returns: Array{Float64, 2}: The sampled points from the PDMP trajectory skeleton.

sample_skeleton(sampler::AbstractPDMP, T::Float64, xinit, vinit; seed, verbose, init_capacity)

Time-horizon variant: advance the PDMP up to time T (instead of generating n_sk events) and return the skeleton.

• Since the number of events is not known a priori, PDMPHistory is allocated with init_capacity and grown by doubling as needed (with copying).
• The return value is forced to satisfy t[end] == T by adding one final point at t=T via deterministic flow (because sample_from_skeleton uses t[end] as the time scale).

failsafe dispatch of sample_skeleton(), admitting scalar initial values, used mainly for 1d case.

sample_skeleton(): generate a PDMP skeleton (event times/states) from a PDMP sampler.

Parameters:
- n_sk (Int): The number of skeleton samples to generate.
- xinit (Array{Float64, 1}): The initial position of the particles.
- vinit (Array{Float64, 1}): The initial velocity of the particles.
- seed (Int): The seed value for random number generation.
- verbose (Bool): Whether to display progress bar during sampling. Default is true.

Returns:
- output: The output state of the sampling process.

sample_skeleton_with_diagnostic(): generate a PDMP skeleton (event times/states) from a PDMP sampler with diagnostic information.

Parameters:
- T (Float64): The time horizon to sample up to.
- xinit (Array{Float64, 1}): The initial position of the particles.
- vinit (Array{Float64, 1}): The initial velocity of the particles.
- B (Int64): The number of batches to split the time horizon into.
- U (Function): The potential function to use for the diagnostic.
- seed (Int): The seed value for random number generation.
- verbose (Bool): Whether to display progress bar during sampling. Default is true.

Returns:
- output: The output state of the sampling process.

PDMPState <: Any

A state in the state space of an `AbstractPDMP`, implemented as a struct with the fields below.

Attributes:
    x (AbstractVector{T}): position
    v (AbstractVector{T}): velocity
    t (Float64): time
    is_active (AbstractVector{Bool}): indicator for the freezing state. Used in the sampling loop for Sticky samplers.

    upper_bound_func (Function): upper bound function
    upper_bound (Union{Nothing, NamedTuple}): upper bound box

    lambda_bar (Float64): upper bound for the Poisson process
    exp_rv (Float64): exponential random variable for the Poisson process
    
    lambda_t (Float64): rate at the current time
    horizon (Float64): horizon
    tp (Float64): time to the next event
    ts (Float64): time spent
    tt (Vector{Float64}): remaining frozen time for each coordinate, with dimension `dim`
    ar (Float64): acceptance rate for the thinning

    adaptive (Bool): adaptive indicator
    accept (Bool): accept indicator for the thinning
    stick_or_thaw_event (Bool): indicator for the sticking or thawing event
    
    errored_bound (Int): count of the number of errors in the upper bound
    rejected (Int): count of the number of rejections in the thinning
    hitting_horizon (Int): count of the number of hits of the horizon
    
Note:
    Keep `flow/∇U/rate/velocity_jump/rng` on the sampler side, and keep `PDMPState` focused on
    the evolving state variables `(x, v, t, ...)`.

RHMC(dim::Int, ∇U;
     mean_duration=nothing,
     refresh_rate=1.0,
     phi=pi/2,
     step_size=0.05,
     tmax=10.0,
     adaptive=false,
     AD_backend="FiniteDiff")

Randomized Hamiltonian Monte Carlo (RHMC) sampler (Bou-Rabee & Sanz-Serna, 2017) as a PDMP.

Between events it follows Hamiltonian dynamics

ẋ = v,    v̇ = -∇U(x)

and at Poisson event times (rate refresh_rate = 1/mean_duration) it refreshes momentum via

v ← cos(phi) v + sin(phi) ξ,   ξ ~ N(0, I)

Notes:

• This implementation uses a velocity-Verlet integrator with step size step_size for the flow.
• It is rejection-free in the exact-flow limit; with a finite step_size this is an approximation.

SpeedUpZigZag(dim::Int, ∇U::Function; grid_size::Int=10, tmax::Float64=1.0, 
    vectorized_bound::Bool=true, signed_bound::Bool=true, adaptive::Bool=true, kwargs...)

arguments for constructor

• dim::Int: dimension of the state space
• ∇U::Function: gradient of the potential energy function
• grid_size::Int: number of grid points used for upper-bound discretization (default: 10)
• tmax::Float64: bound horizon (default: 1.0). If set to 0, tmax is chosen adaptively.
• vectorized_bound::Bool: whether to use a vectorized bound strategy (default: true)
• signed_bound::Bool: whether to use signed-rate bound strategies (default: true)
• adaptive::Bool: whether to adapt the horizon during sampling (default: true)
• kwargs...: additional keyword arguments

attributes of a ZigZag construct

• dim::Int: dimension of the state space
• refresh_rate::Float64: refresh rate
• ∇U::Function: gradient of the potential
• grid_size::Int: number of grid points used for upper-bound discretization
• tmax::Float64: bound horizon
• adaptive::Bool: whether to adapt the horizon during sampling
• vectorized_bound::Bool: whether a vectorized bound strategy is used
• signed_bound::Bool: whether signed-rate strategies are used
• flow::Function: deterministic flow / integrator
• rate: scalar (unsigned) event rate
• rate_vect: vectorized (unsigned) event rate
• signed_rate: scalar signed event rate (if used)
• signed_rate_vect: vectorized signed event rate (if used)
• velocity_jump::Function: velocity update at events
• state: sampler state

StickyZigZag(dim::Int, ∇U::Function; grid_size::Int=10, tmax::Float64=1.0, 
    vectorized_bound::Bool=true, signed_bound::Bool=true, adaptive::Bool=true, kwargs...)

arguments for constructor

• dim::Int: dimension of the parameter space

• ∇U::Function: gradient of the potential function (= negative log-likelihood function)

• κ::Vector{Float64}: thawing rate defined from prior inclusion probability. default is fill(0.5, dim).

• grid_size::Int: number of the grid points for discretization of the parameter space. default is 10.

• tmax::Float64: bound horizon (default: 1.0). If set to 0, tmax is chosen adaptively.

• vectorized_bound::Bool: whether to use a vectorized bound strategy (default: true)

• signed_bound::Bool: whether to use signed-rate bound strategies (default: true)

• adaptive::Bool: whether to adapt the horizon during sampling (default: true)

• kwargs...: additional keyword arguments

attributes of a ZigZag construct

• dim::Int: dimension of the state space
• refresh_rate::Float64: refresh rate
• ∇U::Function: gradient of the potential
• grid_size::Int: number of grid points used for upper-bound discretization
• tmax::Float64: bound horizon
• adaptive::Bool: whether to adapt the horizon during sampling
• vectorized_bound::Bool: whether a vectorized bound strategy is used
• signed_bound::Bool: whether signed-rate strategies are used
• flow::Function: deterministic flow / integrator
• rate: scalar (unsigned) event rate
• rate_vect: vectorized (unsigned) event rate
• signed_rate: scalar signed event rate (if used)
• signed_rate_vect: vectorized signed event rate (if used)
• velocity_jump::Function: velocity update at events
• state: sampler state

ZigZag(dim::Int, ∇U::Function; grid_size::Int=10, tmax::Float64=1.0, 
    vectorized_bound::Bool=true, signed_bound::Bool=true, adaptive::Bool=true, kwargs...)

arguments for constructor

• dim::Int: dimension of the state space
• ∇U::Function: gradient of the potential energy function
• grid_size::Int: number of grid points used for upper-bound discretization (default: 10)
• tmax::Float64: bound horizon (default: 1.0). If set to 0, tmax is chosen adaptively.
• vectorized_bound::Bool: whether to use a vectorized bound strategy (default: true)
• signed_bound::Bool: whether to use signed-rate bound strategies (default: true)
• adaptive::Bool: whether to adapt the horizon during sampling (default: true)
• kwargs...: additional keyword arguments

attributes of a ZigZag construct

• dim::Int: dimension of the state space
• refresh_rate::Float64: refresh rate (currently unused for ZigZag)
• ∇U::Function: gradient of the potential
• grid_size::Int: number of grid points used for upper-bound discretization
• tmax::Float64: bound horizon
• adaptive::Bool: whether to adapt the horizon during sampling
• vectorized_bound::Bool: whether a vectorized bound strategy is used
• signed_bound::Bool: whether signed-rate strategies are used
• flow::Function: deterministic flow / integrator
• rate: scalar (unsigned) event rate
• rate_vect: vectorized (unsigned) event rate
• signed_rate: scalar signed event rate (if used)
• signed_rate_vect: vectorized signed event rate (if used)
• velocity_jump::Function: velocity update at events
• state: sampler state