diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 00000000..1e72b507 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1 @@ +style="blue" diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml new file mode 100644 index 00000000..40d66912 --- /dev/null +++ b/.github/workflows/Format.yml @@ -0,0 +1,28 @@ +name: Format + +on: + push: + branches: + - master + pull_request: + +jobs: + format: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@latest + with: + version: 1 + - name: Format code + run: | + using Pkg + Pkg.add(; name="JuliaFormatter", uuid="98e50ef6-434e-11e9-1051-2b60c6c9e899") + using JuliaFormatter + format("."; verbose=true) + shell: julia --color=yes {0} + - uses: reviewdog/action-suggester@v1 + if: github.event_name == 'pull_request' + with: + tool_name: JuliaFormatter + fail_on_error: true diff --git a/README.md b/README.md index 6fe8d052..e0b5da75 100644 --- a/README.md +++ b/README.md @@ -3,6 +3,7 @@ [![Build Status](https://github.com/TuringLang/AdvancedPS.jl/workflows/CI/badge.svg?branch=master)](https://github.com/TuringLang/AdvancedPS.jl/actions?query=workflow%3ACI%20branch%3Amaster) [![Coverage](https://codecov.io/gh/TuringLang/AdvancedPS.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/TuringLang/AdvancedPS.jl) [![Coverage](https://coveralls.io/repos/github/TuringLang/AdvancedPS.jl/badge.svg?branch=master)](https://coveralls.io/github/TuringLang/AdvancedPS.jl?branch=master) +[![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle) ### Reference diff --git a/src/AdvancedPS.jl b/src/AdvancedPS.jl index 95d0974b..ce2ecf18 100644 --- a/src/AdvancedPS.jl +++ b/src/AdvancedPS.jl @@ -1,10 +1,10 @@ module AdvancedPS -import AbstractMCMC -import Distributions -import Libtask -import Random -import StatsFuns +using AbstractMCMC: AbstractMCMC +using Distributions: Distributions +using Libtask: Libtask +using Random: Random +using StatsFuns: StatsFuns include("resampling.jl") include("container.jl") diff --git a/src/container.jl b/src/container.jl index 7fdaf230..bb0fe76f 100644 --- a/src/container.jl +++ b/src/container.jl @@ -6,7 +6,7 @@ end const Particle = Trace function Trace(f) - ctask = let f=f + ctask = let f = f Libtask.CTask() do res = f() Libtask.produce(nothing) @@ -34,7 +34,7 @@ reset_model(f) = deepcopy(f) delete_retained!(f) = nothing # Task copying version of fork for Trace. -function fork(trace::Trace, isref::Bool = false) +function fork(trace::Trace, isref::Bool=false) newtrace = copy(trace) isref && delete_retained!(newtrace.f) @@ -48,7 +48,7 @@ end # Create new task and copy randomness function forkr(trace::Trace) newf = reset_model(trace.f) - ctask = let f=trace.ctask.task.code + ctask = let f = trace.ctask.task.code Libtask.CTask() do res = f() Libtask.produce(nothing) @@ -105,7 +105,7 @@ end function Base.push!(pc::ParticleContainer, p::Particle) push!(pc.vals, p) push!(pc.logWs, 0.0) - pc + return pc end # clones a theta-particle @@ -116,7 +116,7 @@ function Base.copy(pc::ParticleContainer) # copy weights logWs = copy(pc.logWs) - ParticleContainer(vals, logWs) + return ParticleContainer(vals, logWs) end """ @@ -183,9 +183,9 @@ of the particle `weights`. For Particle Gibbs sampling, one can provide a refere function resample_propagate!( rng::Random.AbstractRNG, pc::ParticleContainer, - randcat = resample_systematic, - ref::Union{Particle, Nothing} = nothing; - weights = getweights(pc) + randcat=resample_systematic, + ref::Union{Particle,Nothing}=nothing; + weights=getweights(pc), ) # check that weights are not NaN @assert !any(isnan, weights) @@ -232,24 +232,24 @@ function resample_propagate!( pc.vals = children reset_logweights!(pc) - pc + return pc end function resample_propagate!( rng::Random.AbstractRNG, pc::ParticleContainer, resampler::ResampleWithESSThreshold, - ref::Union{Particle,Nothing} = nothing; - weights = getweights(pc) + ref::Union{Particle,Nothing}=nothing; + weights=getweights(pc), ) # Compute the effective sample size ``1 / ∑ wᵢ²`` with normalized weights ``wᵢ`` ess = inv(sum(abs2, weights)) if ess ≤ resampler.threshold * length(pc) - resample_propagate!(rng, pc, resampler.resampler, ref; weights = weights) + resample_propagate!(rng, pc, resampler.resampler, ref; weights=weights) end - pc + return pc end """ @@ -290,9 +290,13 @@ function reweight!(pc::ParticleContainer) # The posterior for models with random number of observations is not well-defined. if numdone != 0 - error("mis-aligned execution traces: # particles = ", n, - " # completed trajectories = ", numdone, - ". Please make sure the number of observations is NOT random.") + error( + "mis-aligned execution traces: # particles = ", + n, + " # completed trajectories = ", + numdone, + ". Please make sure the number of observations is NOT random.", + ) end return false diff --git a/src/contrib/AdvancedSMCExtensions.jl b/src/contrib/AdvancedSMCExtensions.jl index 921662b6..243c217e 100644 --- a/src/contrib/AdvancedSMCExtensions.jl +++ b/src/contrib/AdvancedSMCExtensions.jl @@ -26,12 +26,12 @@ Arguments: - `parameters_algs::Tuple{MH}` : An [`MH`](@ref) algorithm, which includes a sample space specification. """ -mutable struct PMMH{space, A<:Tuple} <: InferenceAlgorithm +mutable struct PMMH{space,A<:Tuple} <: InferenceAlgorithm n_iters::Int # number of iterations algs::A # Proposals for state & parameters end function PMMH(n_iters::Int, algs::Tuple, space::Tuple) - return PMMH{space, typeof(algs)}(n_iters, algs) + return PMMH{space,typeof(algs)}(n_iters, algs) end function PMMH(n_iters::Int, smc_alg::SMC, parameter_algs...) return PMMH(n_iters, tuple(parameter_algs..., smc_alg), ()) @@ -40,7 +40,7 @@ end PIMH(n_iters::Int, smc_alg::SMC) = PMMH(n_iters, tuple(smc_alg), ()) function Sampler(alg::PMMH, model::Model, s::Selector) - info = Dict{Symbol, Any}() + info = Dict{Symbol,Any}() spl = Sampler(alg, info, s) alg_str = "PMMH" @@ -51,15 +51,15 @@ function Sampler(alg::PMMH, model::Model, s::Selector) for i in 1:n_samplers sub_alg = alg.algs[i] - if isa(sub_alg, Union{SMC, MH}) + if isa(sub_alg, Union{SMC,MH}) samplers[i] = Sampler(sub_alg, model, Selector(Symbol(typeof(sub_alg)))) else error("[$alg_str] unsupport base sampling algorithm $alg") end if typeof(sub_alg) == MH && sub_alg.n_iters != 1 warn( - "[$alg_str] number of iterations greater than 1" * - "is useless for MH since it is only used for its proposal" + "[$alg_str] number of iterations greater than 1" * + "is useless for MH since it is only used for its proposal", ) end space = union(space, sub_alg.space) @@ -80,10 +80,13 @@ function step(model, spl::Sampler{<:PMMH}, vi::VarInfo, is_first::Bool) old_θ = copy(vi[spl]) @debug "Propose new parameters from proposals..." - for local_spl in spl.info[:samplers][1:end-1] + for local_spl in spl.info[:samplers][1:(end - 1)] @debug "$(typeof(local_spl)) proposing $(local_spl.alg.space)..." propose(model, local_spl, vi) - if local_spl.info[:violating_support] violating_support=true; break end + if local_spl.info[:violating_support] + violating_support = true + break + end new_prior_prob += local_spl.info[:prior_prob] proposal_ratio += local_spl.info[:proposal_ratio] end @@ -98,9 +101,9 @@ function step(model, spl::Sampler{<:PMMH}, vi::VarInfo, is_first::Bool) @debug "Decide whether to accept..." accepted = mh_accept( - spl.info[:old_likelihood_estimate] + spl.info[:old_prior_prob], - new_likelihood_estimate + new_prior_prob, - proposal_ratio, + spl.info[:old_likelihood_estimate] + spl.info[:old_prior_prob], + new_likelihood_estimate + new_prior_prob, + proposal_ratio, ) end @@ -114,13 +117,13 @@ function step(model, spl::Sampler{<:PMMH}, vi::VarInfo, is_first::Bool) return vi, accepted end -function sample( model::Model, - alg::PMMH; - save_state=false, # flag for state saving - resume_from=nothing, # chain to continue - reuse_spl_n=0 # flag for spl re-using - ) - +function sample( + model::Model, + alg::PMMH; + save_state=false, # flag for state saving + resume_from=nothing, # chain to continue + reuse_spl_n=0, # flag for spl re-using +) spl = Sampler(alg, model) if resume_from !== nothing spl.selector = resume_from.info[:spl].selector @@ -134,8 +137,8 @@ function sample( model::Model, time_total = zero(Float64) samples = Array{Sample}(undef, sample_n) weight = 1 / sample_n - for i = 1:sample_n - samples[i] = Sample(weight, Dict{Symbol, Any}()) + for i in 1:sample_n + samples[i] = Sample(weight, Dict{Symbol,Any}()) end # Init parameters @@ -148,22 +151,24 @@ function sample( model::Model, # PMMH steps accept_his = Bool[] - PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0)) - for i = 1:n - @debug "$alg_str stepping..." - time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, i==1) - - if is_accept # accepted => store the new predcits - samples[i].value = Sample(vi, spl).value - else # rejected => store the previous predcits - samples[i] = samples[i - 1] - end - - time_total += time_elapsed - push!(accept_his, is_accept) - if PROGRESS[] - haskey(spl.info, :progress) && ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1) - end + PROGRESS[] && + (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0)) + for i in 1:n + @debug "$alg_str stepping..." + time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, i == 1) + + if is_accept # accepted => store the new predcits + samples[i].value = Sample(vi, spl).value + else # rejected => store the previous predcits + samples[i] = samples[i - 1] + end + + time_total += time_elapsed + push!(accept_his, is_accept) + if PROGRESS[] + haskey(spl.info, :progress) && + ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1) + end end println("[$alg_str] Finished with") @@ -172,18 +177,17 @@ function sample( model::Model, println(" Accept rate = $accept_rate;") if resume_from !== nothing # concat samples - pushfirst!(samples, resume_from.info[:samples]...) + pushfirst!(samples, resume_from.info[:samples]...) end c = Chain(-Inf, samples) # wrap the result by Chain if save_state # save state - c = save(c, spl, model, vi, samples) + c = save(c, spl, model, vi, samples) end - c + return c end - #### #### IMCMC Sampler. #### @@ -212,41 +216,45 @@ Arguments: A paper on this can be found [here](https://arxiv.org/abs/1602.05128). """ -mutable struct IPMCMC{T, F} <: InferenceAlgorithm - n_particles::Int # number of particles used - n_iters::Int # number of iterations - n_nodes::Int # number of nodes running SMC and CSMC - n_csmc_nodes::Int # number of nodes CSMC - resampler::F # function to resample - space::Set{T} # sampling space, emtpy means all +mutable struct IPMCMC{T,F} <: InferenceAlgorithm + n_particles::Int # number of particles used + n_iters::Int # number of iterations + n_nodes::Int # number of nodes running SMC and CSMC + n_csmc_nodes::Int # number of nodes CSMC + resampler::F # function to resample + space::Set{T} # sampling space, emtpy means all end IPMCMC(n1::Int, n2::Int) = IPMCMC(n1, n2, 32, 16, resample_systematic, Set()) -IPMCMC(n1::Int, n2::Int, n3::Int) = IPMCMC(n1, n2, n3, Int(ceil(n3/2)), resample_systematic, Set()) -IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int) = IPMCMC(n1, n2, n3, n4, resample_systematic, Set()) +function IPMCMC(n1::Int, n2::Int, n3::Int) + return IPMCMC(n1, n2, n3, Int(ceil(n3 / 2)), resample_systematic, Set()) +end +function IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int) + return IPMCMC(n1, n2, n3, n4, resample_systematic, Set()) +end function IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int, space...) - _space = isa(space, Symbol) ? Set([space]) : Set(space) - IPMCMC(n1, n2, n3, n4, resample_systematic, _space) + _space = isa(space, Symbol) ? Set([space]) : Set(space) + return IPMCMC(n1, n2, n3, n4, resample_systematic, _space) end function Sampler(alg::IPMCMC, s::Selector) - info = Dict{Symbol, Any}() - spl = Sampler(alg, info, s) - # Create SMC and CSMC nodes - samplers = Array{Sampler}(undef, alg.n_nodes) - # Use resampler_threshold=1.0 for SMC since adaptive resampling is invalid in this setting - default_CSMC = CSMC(alg.n_particles, 1, alg.resampler, alg.space) - default_SMC = SMC(alg.n_particles, alg.resampler, 1.0, false, alg.space) - - for i in 1:alg.n_csmc_nodes - samplers[i] = Sampler(default_CSMC, Selector(Symbol(typeof(default_CSMC)))) - end - for i in (alg.n_csmc_nodes+1):alg.n_nodes - samplers[i] = Sampler(default_SMC, Selector(Symbol(typeof(default_CSMC)))) - end - - info[:samplers] = samplers - - return spl + info = Dict{Symbol,Any}() + spl = Sampler(alg, info, s) + # Create SMC and CSMC nodes + samplers = Array{Sampler}(undef, alg.n_nodes) + # Use resampler_threshold=1.0 for SMC since adaptive resampling is invalid in this setting + default_CSMC = CSMC(alg.n_particles, 1, alg.resampler, alg.space) + default_SMC = SMC(alg.n_particles, alg.resampler, 1.0, false, alg.space) + + for i in 1:(alg.n_csmc_nodes) + samplers[i] = Sampler(default_CSMC, Selector(Symbol(typeof(default_CSMC)))) + end + for i in (alg.n_csmc_nodes + 1):(alg.n_nodes) + samplers[i] = Sampler(default_SMC, Selector(Symbol(typeof(default_CSMC)))) + end + + info[:samplers] = samplers + + return spl end function step(model, spl::Sampler{<:IPMCMC}, VarInfos::Array{VarInfo}, is_first::Bool) @@ -254,16 +262,16 @@ function step(model, spl::Sampler{<:IPMCMC}, VarInfos::Array{VarInfo}, is_first: log_zs = zeros(spl.alg.n_nodes) # Run SMC & CSMC nodes - for j in 1:spl.alg.n_nodes + for j in 1:(spl.alg.n_nodes) reset_num_produce!(VarInfos[j]) VarInfos[j] = step(model, spl.info[:samplers][j], VarInfos[j])[1] log_zs[j] = spl.info[:samplers][j].info[:logevidence][end] end # Resampling of CSMC nodes indices - conditonal_nodes_indices = collect(1:spl.alg.n_csmc_nodes) - unconditonal_nodes_indices = collect(spl.alg.n_csmc_nodes+1:spl.alg.n_nodes) - for j in 1:spl.alg.n_csmc_nodes + conditonal_nodes_indices = collect(1:(spl.alg.n_csmc_nodes)) + unconditonal_nodes_indices = collect((spl.alg.n_csmc_nodes + 1):(spl.alg.n_nodes)) + for j in 1:(spl.alg.n_csmc_nodes) # Select a new conditional node by simulating cj log_ksi = vcat(log_zs[unconditonal_nodes_indices], log_zs[j]) ksi = exp.(log_ksi .- maximum(log_ksi)) @@ -276,51 +284,53 @@ function step(model, spl::Sampler{<:IPMCMC}, VarInfos::Array{VarInfo}, is_first: end nodes_permutation = vcat(conditonal_nodes_indices, unconditonal_nodes_indices) - VarInfos[nodes_permutation] + return VarInfos[nodes_permutation] end function sample(model::Model, alg::IPMCMC) + spl = Sampler(alg) - spl = Sampler(alg) - - # Number of samples to store - sample_n = alg.n_iters * alg.n_csmc_nodes - - # Init samples - time_total = zero(Float64) - samples = Array{Sample}(undef, sample_n) - weight = 1 / sample_n - for i = 1:sample_n - samples[i] = Sample(weight, Dict{Symbol, Any}()) - end - - # Init parameters - vi = empty!(VarInfo(model)) - VarInfos = Array{VarInfo}(undef, spl.alg.n_nodes) - for j in 1:spl.alg.n_nodes - VarInfos[j] = deepcopy(vi) - end - n = spl.alg.n_iters - - # IPMCMC steps - if PROGRESS[] spl.info[:progress] = ProgressMeter.Progress(n, 1, "[IPMCMC] Sampling...", 0) end - for i = 1:n - @debug "IPMCMC stepping..." - time_elapsed = @elapsed VarInfos = step(model, spl, VarInfos, i==1) - - # Save each CSMS retained path as a sample - for j in 1:spl.alg.n_csmc_nodes - samples[(i-1)*alg.n_csmc_nodes+j].value = Sample(VarInfos[j], spl).value + # Number of samples to store + sample_n = alg.n_iters * alg.n_csmc_nodes + + # Init samples + time_total = zero(Float64) + samples = Array{Sample}(undef, sample_n) + weight = 1 / sample_n + for i in 1:sample_n + samples[i] = Sample(weight, Dict{Symbol,Any}()) end - time_total += time_elapsed + # Init parameters + vi = empty!(VarInfo(model)) + VarInfos = Array{VarInfo}(undef, spl.alg.n_nodes) + for j in 1:(spl.alg.n_nodes) + VarInfos[j] = deepcopy(vi) + end + n = spl.alg.n_iters + + # IPMCMC steps if PROGRESS[] - haskey(spl.info, :progress) && ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1) + spl.info[:progress] = ProgressMeter.Progress(n, 1, "[IPMCMC] Sampling...", 0) end - end + for i in 1:n + @debug "IPMCMC stepping..." + time_elapsed = @elapsed VarInfos = step(model, spl, VarInfos, i == 1) + + # Save each CSMS retained path as a sample + for j in 1:(spl.alg.n_csmc_nodes) + samples[(i - 1) * alg.n_csmc_nodes + j].value = Sample(VarInfos[j], spl).value + end - println("[IPMCMC] Finished with") - println(" Running time = $time_total;") + time_total += time_elapsed + if PROGRESS[] + haskey(spl.info, :progress) && + ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1) + end + end + + println("[IPMCMC] Finished with") + println(" Running time = $time_total;") - Chain(0.0, samples) # wrap the result by Chain + return Chain(0.0, samples) # wrap the result by Chain end diff --git a/src/resampling.jl b/src/resampling.jl index 4bb21046..b2797461 100644 --- a/src/resampling.jl +++ b/src/resampling.jl @@ -7,13 +7,13 @@ # - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf # Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering -struct ResampleWithESSThreshold{R, T<:Real} +struct ResampleWithESSThreshold{R,T<:Real} resampler::R threshold::T end -function ResampleWithESSThreshold(resampler = resample_systematic) - ResampleWithESSThreshold(resampler, 0.5) +function ResampleWithESSThreshold(resampler=resample_systematic) + return ResampleWithESSThreshold(resampler, 0.5) end # More stable, faster version of rand(Categorical) @@ -30,9 +30,7 @@ function randcat(rng::Random.AbstractRNG, p::AbstractVector{<:Real}) end function resample_multinomial( - rng::Random.AbstractRNG, - w::AbstractVector{<:Real}, - num_particles::Integer = length(w), + rng::Random.AbstractRNG, w::AbstractVector{<:Real}, num_particles::Integer=length(w) ) return rand(rng, Distributions.sampler(Distributions.Categorical(w)), num_particles) end @@ -40,7 +38,7 @@ end function resample_residual( rng::Random.AbstractRNG, w::AbstractVector{<:Real}, - num_particles::Integer = length(weights), + num_particles::Integer=length(weights), ) # Pre-allocate array for resampled particles indices = Vector{Int}(undef, num_particles) @@ -67,7 +65,6 @@ function resample_residual( return indices end - """ resample_stratified(rng, weights, n) @@ -81,9 +78,7 @@ i.e., `xᵢ = j` if and only if ``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``. """ function resample_stratified( - rng::Random.AbstractRNG, - weights::AbstractVector{<:Real}, - n::Integer = length(weights), + rng::Random.AbstractRNG, weights::AbstractVector{<:Real}, n::Integer=length(weights) ) # check input m = length(weights) @@ -130,9 +125,7 @@ normalized `weights`, i.e., `xᵢ = j` if and only if ``uᵢ \\in [\\sum_{s=1}^{j-1} weights_{s}, \\sum_{s=1}^{j} weights_{s})``. """ function resample_systematic( - rng::Random.AbstractRNG, - weights::AbstractVector{<:Real}, - n::Integer = length(weights), + rng::Random.AbstractRNG, weights::AbstractVector{<:Real}, n::Integer=length(weights) ) # check input m = length(weights) diff --git a/src/smc.jl b/src/smc.jl index da5523ae..1beab676 100644 --- a/src/smc.jl +++ b/src/smc.jl @@ -31,17 +31,14 @@ function AbstractMCMC.sample(model::AbstractMCMC.AbstractModel, sampler::SMC; kw end function AbstractMCMC.sample( - rng::Random.AbstractRNG, - model::AbstractMCMC.AbstractModel, - sampler::SMC; - kwargs... + rng::Random.AbstractRNG, model::AbstractMCMC.AbstractModel, sampler::SMC; kwargs... ) if !isempty(kwargs) @warn "keyword arguments $(keys(kwargs)) are not supported by `SMC`" end # Create a set of particles. - particles = ParticleContainer([Trace(model) for _ in 1:sampler.nparticles]) + particles = ParticleContainer([Trace(model) for _ in 1:(sampler.nparticles)]) # Perform particle sweep. logevidence = sweep!(rng, particles, sampler.resampler) @@ -83,13 +80,10 @@ struct PGSample{T,L} end function AbstractMCMC.step( - rng::Random.AbstractRNG, - model::AbstractMCMC.AbstractModel, - sampler::PG; - kwargs..., + rng::Random.AbstractRNG, model::AbstractMCMC.AbstractModel, sampler::PG; kwargs... ) # Create a new set of particles. - particles = ParticleContainer([Trace(model) for _ in 1:sampler.nparticles]) + particles = ParticleContainer([Trace(model) for _ in 1:(sampler.nparticles)]) # Perform a particle sweep. logevidence = sweep!(rng, particles, sampler.resampler) @@ -105,7 +99,7 @@ function AbstractMCMC.step( model::AbstractMCMC.AbstractModel, sampler::PG, state::PGState; - kwargs... + kwargs..., ) # Create a new set of particles. nparticles = sampler.nparticles diff --git a/test/container.jl b/test/container.jl index 2b2f83c8..7313ac79 100644 --- a/test/container.jl +++ b/test/container.jl @@ -27,8 +27,8 @@ # Initial state. @test pc.logWs == zeros(3) - @test AdvancedPS.getweights(pc) == fill(1/3, 3) - @test all(AdvancedPS.getweight(pc, i) == 1/3 for i in 1:3) + @test AdvancedPS.getweights(pc) == fill(1 / 3, 3) + @test all(AdvancedPS.getweight(pc, i) == 1 / 3 for i in 1:3) @test AdvancedPS.logZ(pc) ≈ log(3) @test AdvancedPS.effectiveSampleSize(pc) == 3 @@ -36,21 +36,26 @@ AdvancedPS.reweight!(pc) @test pc.logWs == logps @test AdvancedPS.getweights(pc) ≈ exp.(logps) ./ sum(exp, logps) - @test all(AdvancedPS.getweight(pc, i) ≈ exp(logps[i]) / sum(exp, logps) for i in 1:3) + @test all( + AdvancedPS.getweight(pc, i) ≈ exp(logps[i]) / sum(exp, logps) for i in 1:3 + ) @test AdvancedPS.logZ(pc) ≈ log(sum(exp, logps)) # Reweight particles. AdvancedPS.reweight!(pc) @test pc.logWs == 2 .* logps @test AdvancedPS.getweights(pc) == exp.(2 .* logps) ./ sum(exp, 2 .* logps) - @test all(AdvancedPS.getweight(pc, i) ≈ exp(2 * logps[i]) / sum(exp, 2 .* logps) for i in 1:3) + @test all( + AdvancedPS.getweight(pc, i) ≈ exp(2 * logps[i]) / sum(exp, 2 .* logps) for + i in 1:3 + ) @test AdvancedPS.logZ(pc) ≈ log(sum(exp, 2 .* logps)) # Resample and propagate particles. AdvancedPS.resample_propagate!(Random.GLOBAL_RNG, pc) @test pc.logWs == zeros(3) - @test AdvancedPS.getweights(pc) == fill(1/3, 3) - @test all(AdvancedPS.getweight(pc, i) == 1/3 for i in 1:3) + @test AdvancedPS.getweights(pc) == fill(1 / 3, 3) + @test all(AdvancedPS.getweight(pc, i) == 1 / 3 for i in 1:3) @test AdvancedPS.logZ(pc) ≈ log(3) @test AdvancedPS.effectiveSampleSize(pc) == 3 @@ -58,7 +63,9 @@ AdvancedPS.reweight!(pc) @test pc.logWs ⊆ logps @test AdvancedPS.getweights(pc) == exp.(pc.logWs) ./ sum(exp, pc.logWs) - @test all(AdvancedPS.getweight(pc, i) ≈ exp(pc.logWs[i]) / sum(exp, pc.logWs) for i in 1:3) + @test all( + AdvancedPS.getweight(pc, i) ≈ exp(pc.logWs[i]) / sum(exp, pc.logWs) for i in 1:3 + ) @test AdvancedPS.logZ(pc) ≈ log(sum(exp, pc.logWs)) # Increase unnormalized logarithmic weights. diff --git a/test/resampling.jl b/test/resampling.jl index f72f07bf..8f7eb8cb 100644 --- a/test/resampling.jl +++ b/test/resampling.jl @@ -5,12 +5,12 @@ resSystematic = AdvancedPS.resample_systematic(rng, D, num_samples) resStratified = AdvancedPS.resample_stratified(rng, D, num_samples) - resMultinomial= AdvancedPS.resample_multinomial(rng, D, num_samples) - resResidual = AdvancedPS.resample_residual(rng, D, num_samples) + resMultinomial = AdvancedPS.resample_multinomial(rng, D, num_samples) + resResidual = AdvancedPS.resample_residual(rng, D, num_samples) AdvancedPS.resample_systematic(rng, D) - @test sum(resSystematic .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples - @test sum(resStratified .== 2) ≈ (num_samples * 0.4) atol=1e-3*num_samples - @test sum(resMultinomial .== 2) ≈ (num_samples * 0.4) atol=1e-2*num_samples - @test sum(resResidual .== 2) ≈ (num_samples * 0.4) atol=1e-2*num_samples + @test sum(resSystematic .== 2) ≈ (num_samples * 0.4) atol = 1e-3 * num_samples + @test sum(resStratified .== 2) ≈ (num_samples * 0.4) atol = 1e-3 * num_samples + @test sum(resMultinomial .== 2) ≈ (num_samples * 0.4) atol = 1e-2 * num_samples + @test sum(resResidual .== 2) ≈ (num_samples * 0.4) atol = 1e-2 * num_samples end diff --git a/test/runtests.jl b/test/runtests.jl index 7a3cd3c7..5497a17e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,13 @@ using Random using Test @testset "AdvancedPS.jl" begin - @testset "Resampling tests" begin include("resampling.jl") end - @testset "Container tests" begin include("container.jl") end - @testset "SMC and PG tests" begin include("smc.jl") end + @testset "Resampling tests" begin + include("resampling.jl") + end + @testset "Container tests" begin + include("container.jl") + end + @testset "SMC and PG tests" begin + include("smc.jl") + end end diff --git a/test/smc.jl b/test/smc.jl index 2726b718..1f98eb37 100644 --- a/test/smc.jl +++ b/test/smc.jl @@ -7,12 +7,12 @@ sampler = AdvancedPS.SMC(15, 0.6) @test sampler.nparticles == 15 @test sampler.resampler === - AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_systematic, 0.6) + AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_systematic, 0.6) sampler = AdvancedPS.SMC(20, AdvancedPS.resample_multinomial, 0.6) @test sampler.nparticles == 20 @test sampler.resampler === - AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_multinomial, 0.6) + AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_multinomial, 0.6) sampler = AdvancedPS.SMC(25, AdvancedPS.resample_systematic) @test sampler.nparticles == 25 @@ -39,7 +39,7 @@ m.b = b = rand(Normal(a, 1)) # Second observation. - AdvancedPS.observe(Normal(b, 2), 1.5) + return AdvancedPS.observe(Normal(b, 2), 1.5) end sample(NormalModel(), AdvancedPS.SMC(100)) @@ -87,7 +87,7 @@ m.c = rand(Beta()) # Second observation. - AdvancedPS.observe(Bernoulli(x / 2), 0) + return AdvancedPS.observe(Bernoulli(x / 2), 0) end chains_smc = sample(TestModel(), AdvancedPS.SMC(100)) @@ -104,12 +104,12 @@ sampler = AdvancedPS.PG(60, 0.6) @test sampler.nparticles == 60 @test sampler.resampler === - AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_systematic, 0.6) + AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_systematic, 0.6) sampler = AdvancedPS.PG(80, AdvancedPS.resample_multinomial, 0.6) @test sampler.nparticles == 80 @test sampler.resampler === - AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_multinomial, 0.6) + AdvancedPS.ResampleWithESSThreshold(AdvancedPS.resample_multinomial, 0.6) sampler = AdvancedPS.PG(100, AdvancedPS.resample_systematic) @test sampler.nparticles == 100 @@ -141,7 +141,7 @@ m.c = rand(Beta()) # Second observation. - AdvancedPS.observe(Bernoulli(x / 2), 0) + return AdvancedPS.observe(Bernoulli(x / 2), 0) end chains_pg = sample(TestModel(), AdvancedPS.PG(10), 100) @@ -208,4 +208,3 @@ end # check_MoGtest_default(chain2, atol = 0.2) # end # end -