Sampler failure when num of samples = 1

[yizhang-yiz]

Not sure if this has been reported. Setting n_samples = 1 in the minimal example stops the sampler from working properly.

using AdvancedHMC, Distributions, ForwardDiff
using LinearAlgebra

function debug()
    # Choose parameter dimensionality and initial parameter value
    D = 10; initial_θ = rand(D)

    # Define the target distribution
    ℓπ(θ) = logpdf(MvNormal(zeros(D), I), θ)

    # Set the number of samples to draw and warmup iterations
    n_samples, n_adapts = 1, 1000

    # Define a Hamiltonian system
    metric = DiagEuclideanMetric(D)
    hamiltonian = Hamiltonian(metric, ℓπ, ForwardDiff)

    # Define a leapfrog solver, with initial step size chosen heuristically
    initial_ϵ = find_good_stepsize(hamiltonian, initial_θ)
    integrator = Leapfrog(initial_ϵ)

    # Define an HMC sampler, with the following components
    #   - multinomial sampling scheme,
    #   - generalised No-U-Turn criteria, and
    #   - windowed adaption for step-size and diagonal mass matrix
    proposal = NUTS{MultinomialTS, GeneralisedNoUTurn}(integrator)
    adaptor = StanHMCAdaptor(MassMatrixAdaptor(metric), StepSizeAdaptor(0.8, integrator))

    # Run the sampler to draw samples from the specified Gaussian, where
    #   - `samples` will store the samples
    #   - `stats` will store diagnostic statistics for each sample
    samples, stats = sample(hamiltonian, proposal, initial_θ, 1, adaptor, n_adapts; progress=true)
end

Checked on v1.6.2 & v1.6.3, I'm getting

julia> debug()
ERROR: ArgumentError: reducing over an empty collection is not allowed
Stacktrace:
  [1] _empty_reduce_error()
    @ Base ./reduce.jl:299
  [2] mapreduce_empty(f::Function, op::Function, T::Type)
    @ Base ./reduce.jl:342
  [3] reduce_empty(op::Base.MappingRF{AdvancedHMC.var"#22#23"{Vector{Float64}}, typeof(Base.add_sum)}, #unused#::Type{Int64})
    @ Base ./reduce.jl:329
  [4] reduce_empty_iter
    @ ./reduce.jl:355 [inlined]
  [5] mapreduce_empty_iter(f::Function, op::Function, itr::UnitRange{Int64}, ItrEltype::Base.HasEltype)
    @ Base ./reduce.jl:351
  [6] _mapreduce(f::AdvancedHMC.var"#22#23"{Vector{Float64}}, op::typeof(Base.add_sum), #unused#::IndexLinear, A::UnitRange{Int64})
    @ Base ./reduce.jl:400
  [7] _mapreduce_dim
    @ ./reducedim.jl:318 [inlined]
  [8] #mapreduce#672
    @ ./reducedim.jl:310 [inlined]
  [9] mapreduce
    @ ./reducedim.jl:310 [inlined]
 [10] #_sum#682
    @ ./reducedim.jl:878 [inlined]
 [11] _sum
    @ ./reducedim.jl:878 [inlined]
 [12] #sum#680
    @ ./reducedim.jl:874 [inlined]
 [13] _mean(f::AdvancedHMC.var"#22#23"{Vector{Float64}}, A::UnitRange{Int64}, dims::Colon)
    @ Statistics /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Statistics/src/Statistics.jl:170
 [14] #mean#1
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Statistics/src/Statistics.jl:104 [inlined]
 [15] mean
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/Statistics/src/Statistics.jl:104 [inlined]
 [16] EBFMI(Es::Vector{Float64})
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/yd6UP/src/diagnosis.jl:2
 [17] sample(rng::Random._GLOBAL_RNG, h::Hamiltonian{DiagEuclideanMetric{Float64, Vector{Float64}}, var"#ℓπ#1"{Int64}, AdvancedHMC.var"#∂ℓπ∂θ#66"{var"#ℓπ#1"{Int64}}}, κ::HMCKernel{AdvancedHMC.FullMomentumRefreshment, Trajectory{MultinomialTS, Leapfrog{Float64}, GeneralisedNoUTurn{Float64}}}, θ::Vector{Float64}, n_samples::Int64, adaptor::StanHMCAdaptor{WelfordVar{Float64, Vector{Float64}}, NesterovDualAveraging{Float64}}, n_adapts::Int64; drop_warmup::Bool, verbose::Bool, progress::Bool, pm_next!::typeof(AdvancedHMC.pm_next!))
    @ AdvancedHMC ~/.julia/packages/AdvancedHMC/yd6UP/src/sampler.jl:199
 [18] #sample#26
    @ ~/.julia/packages/AdvancedHMC/yd6UP/src/sampler.jl:118 [inlined]
 [19] debug()
    @ Main ./REPL[4]:29
 [20] top-level scope
    @ REPL[5]:1

Am I missing something here?