SMC throws BoundsError when discard_initial > 1 or thinning > 1

[vandenman]

This following works on Turing v0.21.1:

using Turing
# Define a simple Normal model with unknown mean and variance.
@model function gdemo(x, y)
	s² ~ InverseGamma(2, 3)
	m ~ Normal(0, sqrt(s²))
	x ~ Normal(m, sqrt(s²))
	y ~ Normal(m, sqrt(s²))
end

# works
chn = sample(gdemo(1.5, 2), SMC(), 10; discard_initial = 1)

changing discard_initial (or thinning) to a value larger than 1 and an error is thrown:

# fails
chn = sample(gdemo(1.5, 2), SMC(), 10; discard_initial = 2)

ERROR: BoundsError: attempt to access 10-element Vector{AdvancedPS.Trace{Turing.Essential.TracedModel{AbstractMCMC.AbstractSampler, DynamicPPL.AbstractVarInfo, DynamicPPL.Model, Tuple}}} at index [11]
Stacktrace:
  [1] getindex
    @ ./array.jl:861 [inlined]
  [2] step(::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(gdemo), (:x, :y), (), (), Tuple{Float64, Int64}, Tuple{}, DynamicPPL.DefaultContext}, spl::DynamicPPL.Sampler{SMC{(), AdvancedPS.ResampleWithESSThreshold{typeof(AdvancedPS.resample_systematic), Float64}}}, state::Turing.Inference.SMCState{AdvancedPS.ParticleContainer{AdvancedPS.Trace{Turing.Essential.TracedModel{AbstractMCMC.AbstractSampler, DynamicPPL.AbstractVarInfo, DynamicPPL.Model, Tuple}}}, Float64}; kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nparticles,), Tuple{Int64}}})
    @ Turing.Inference ~/.julia/packages/Turing/uMoX1/src/inference/AdvancedSMC.jl:147
  [3] macro expansion
    @ ~/.julia/packages/AbstractMCMC/0eT8o/src/sample.jl:163 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/ProgressLogging/6KXlp/src/ProgressLogging.jl:328 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/AbstractMCMC/0eT8o/src/logging.jl:9 [inlined]
  [6] mcmcsample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(gdemo), (:x, :y), (), (), Tuple{Float64, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{SMC{(), AdvancedPS.ResampleWithESSThreshold{typeof(AdvancedPS.resample_systematic), Float64}}}, N::Int64; progress::Bool, progressname::String, callback::Nothing, discard_initial::Int64, thinning::Int64, chain_type::Type, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:nparticles,), Tuple{Int64}}})
    @ AbstractMCMC ~/.julia/packages/AbstractMCMC/0eT8o/src/sample.jl:111
  [7] sample(rng::Random._GLOBAL_RNG, model::DynamicPPL.Model{typeof(gdemo), (:x, :y), (), (), Tuple{Float64, Int64}, Tuple{}, DynamicPPL.DefaultContext}, sampler::DynamicPPL.Sampler{SMC{(), AdvancedPS.ResampleWithESSThreshold{typeof(AdvancedPS.resample_systematic), Float64}}}, N::Int64; chain_type::Type, resume_from::Nothing, progress::Bool, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol}, NamedTuple{(:discard_initial,), Tuple{Int64}}})
    @ Turing.Inference ~/.julia/packages/Turing/uMoX1/src/inference/AdvancedSMC.jl:91
  [8] #sample#2
    @ ~/.julia/packages/Turing/uMoX1/src/inference/Inference.jl:145 [inlined]
  [9] #sample#1
    @ ~/.julia/packages/Turing/uMoX1/src/inference/Inference.jl:135 [inlined]
 [10] top-level scope
    @ ~/github/juliaPkgs/issues/Turing/smc_issue.jl:14

perhaps that in this function nparticles

https://github.com/TuringLang/Turing.jl/blob/715dd985ac331f5df26338e6b30ac686015610c4/src/inference/AdvancedSMC.jl#L91-L98

needs to be adjusted for discard_initial and thinning?

The issue can also be replicated by calling mcmcsample directly:

import Random
n_iter = 10
AbstractMCMC.mcmcsample(
	Random.GLOBAL_RNG,
	gdemo(1.5, 2),
	DynamicPPL.Sampler(SMC()),
	n_iter;
	chain_type=MCMCChains.Chains, 
	progress=Turing.PROGRESS[],
	nparticles=n_iter,
	discard_initial = 1
)

AbstractMCMC.mcmcsample(
	Random.GLOBAL_RNG,
	gdemo(1.5, 2),
	DynamicPPL.Sampler(SMC()),
	n_iter;
	chain_type=MCMCChains.Chains, 
	progress=Turing.PROGRESS[],
	nparticles=n_iter,
	discard_initial = 2
)

No error is thrown if nparticles is adjusted for discard_initial:

AbstractMCMC.mcmcsample(
	Random.GLOBAL_RNG,
	gdemo(1.5, 2),
	DynamicPPL.Sampler(SMC()),
	n_iter;
	chain_type=MCMCChains.Chains, 
	progress=Turing.PROGRESS[],
	nparticles= n_iter + 1,
	discard_initial = 2
)

If this is as simple as modifying AbstractMCMC.sample(..., spl::Sampler{<:SMC}; ...) so that nparticles = N * thinning + discard_initial - 1, I'd be happy to open a PR!

[devmotion]

SMC is a bit special, these options don't make sense. The returned samples are just a set of particles that were propagated together through the model. So there's not any burn-in phase that you might want to discard.

[vandenman]

SMC is a bit special, these options don't make sense.

In that case, shouldn't a DomainError be thrown for the sake of user-friendliness whenever these options are set?

[devmotion]

IMO an ArgumentError would be more appropriate, usually DomainError is used for actual domain errors such as sqrt(-1). However, I think it would be more consistent with other samplers and AbstractMCMC to just ignore unsupported keyword arguments (possibly showing a warning in this case but I don't think the keyword arguments are commonly used with SMC).