Skip to content

Commit b56663d

Browse files
committed
improve coverage and error messages
1 parent 07525bc commit b56663d

2 files changed

Lines changed: 184 additions & 37 deletions

File tree

src/MarkovChainMonteCarlo.jl

Lines changed: 85 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -210,14 +210,7 @@ function log_transition_density(
210210
θ_to;
211211
stepsize::FT = 1.0,
212212
) where {MHS <: AdvancedMH.MHSampler, FT <: AbstractFloat}
213-
throw(
214-
ArgumentError(
215-
"log_transition_density not implemented for $(MHS). Every MHSampler in this module " *
216-
"must implement its own log_transition_density(sampler, model, θ_from, θ_to; stepsize) " *
217-
"= log q(θ_to | θ_from); there is no generic/symmetric fallback, since silently assuming " *
218-
"one caused the prior-double-counting bug this interface replaces.",
219-
),
220-
)
213+
_throw_log_transition_density_not_implemented(MHS)
221214
end
222215

223216
function AdvancedMH.logratio_proposal_density(
@@ -819,11 +812,8 @@ $(TYPEDSIGNATURES)
819812
Fraction of MC proposals in `chain` which were accepted (according to Metropolis-Hastings.)
820813
"""
821814
function accept_ratio(chain::MCMCChains.Chains)
822-
if :accepted in names(chain, :internals)
823-
return mean(chain, :accepted)
824-
else
825-
error("MH `:accepted` not recorded in MCMC chain — available internals: $(names(chain, :internals)).")
826-
end
815+
:accepted in names(chain, :internals) || _throw_accepted_not_recorded(chain)
816+
return mean(chain, :accepted)
827817
end
828818

829819
function _find_mcmc_step_log(mcmc::MCMCWrapper)
@@ -886,12 +876,7 @@ function optimize_stepsize(
886876
n_evals = 0
887877
function acc_at(stepsize)
888878
n_evals += 1
889-
if n_evals > max_iter
890-
error(
891-
"optimize_stepsize: acceptance rate did not reach target $(target_acc) ± $(tol) " *
892-
"within $(max_iter) iterations — last stepsize tried: $(round(stepsize; sigdigits = 3)).",
893-
)
894-
end
879+
n_evals <= max_iter || _throw_max_iter_exceeded(n_evals, max_iter, target_acc, tol, stepsize)
895880
trial_chain = sample(rng, mcmc, N; stepsize = stepsize, sample_kwargs...)
896881
acc_ratio = accept_ratio(trial_chain)
897882
_find_mcmc_step_log(n_evals, stepsize, acc_ratio, trial_chain)
@@ -919,14 +904,8 @@ function optimize_stepsize(
919904
n_expand = 0
920905
while acc_hi > target_acc
921906
n_expand += 1
922-
if n_expand > max_expansions
923-
error(
924-
"optimize_stepsize: acceptance rate stayed above target $(target_acc) ± $(tol) " *
925-
"after $(max_expansions) doublings (up to stepsize $(round(hi; sigdigits = 3))) " *
926-
"without crossing it — the acceptance-vs-stepsize relationship may not be monotonic " *
927-
"in this range rather than simply needing a larger stepsize.",
928-
)
929-
end
907+
n_expand <= max_expansions ||
908+
_throw_max_expansions_exceeded(:above, n_expand, max_expansions, target_acc, tol, hi)
930909
lo = hi
931910
hi *= 2
932911
acc_hi = acc_at(hi)
@@ -938,15 +917,8 @@ function optimize_stepsize(
938917
n_expand = 0
939918
while acc_lo < target_acc
940919
n_expand += 1
941-
if n_expand > max_expansions
942-
error(
943-
"optimize_stepsize: acceptance rate stayed below target $(target_acc) ± $(tol) " *
944-
"after $(max_expansions) halvings (down to stepsize $(round(lo; sigdigits = 3))) " *
945-
"without crossing it — the acceptance-vs-stepsize relationship may not be monotonic " *
946-
"in this range (e.g. a numerically noisy log-density at very small stepsizes) rather " *
947-
"than simply needing a smaller stepsize.",
948-
)
949-
end
920+
n_expand <= max_expansions ||
921+
_throw_max_expansions_exceeded(:below, n_expand, max_expansions, target_acc, tol, lo)
950922
hi = lo
951923
lo /= 2
952924
acc_lo = acc_at(lo)
@@ -1058,6 +1030,83 @@ Suggestion:
10581030
"""))
10591031
end
10601032

1033+
@noinline function _throw_log_transition_density_not_implemented(sampler_type)
1034+
throw(ArgumentError("""
1035+
`log_transition_density` is not implemented for sampler type $(sampler_type).
1036+
1037+
Expected:
1038+
Every `AdvancedMH.MHSampler` used with this module must implement its own
1039+
`log_transition_density(sampler, model, θ_from, θ_to; stepsize)`, returning
1040+
`log q(θ_to | θ_from)`. There is no generic/symmetric fallback, since silently
1041+
assuming one previously caused a prior-double-counting bug.
1042+
1043+
Got:
1044+
typeof(sampler) = $(sampler_type)
1045+
1046+
Suggestion:
1047+
Define `log_transition_density(sampler::$(sampler_type), model, θ_from, θ_to; stepsize = 1.0)`
1048+
for your sampler type.
1049+
"""))
1050+
end
1051+
1052+
@noinline function _throw_accepted_not_recorded(chain::MCMCChains.Chains)
1053+
throw(ArgumentError("""
1054+
Chain is missing the `:accepted` internal required by `accept_ratio`.
1055+
1056+
Expected:
1057+
A `Chains` object produced by this module's `sample`/`optimize_stepsize`, which always
1058+
records `:accepted` as an internal.
1059+
1060+
Got:
1061+
available internals = $(names(chain, :internals))
1062+
1063+
Suggestion:
1064+
Only pass `Chains` objects returned by this module's `sample`/`optimize_stepsize`; do
1065+
not construct or filter the internals section by hand.
1066+
"""))
1067+
end
1068+
1069+
@noinline function _throw_max_iter_exceeded(n_evals, max_iter, target_acc, tol, stepsize)
1070+
throw(ArgumentError("""
1071+
optimize_stepsize did not reach the target acceptance rate within the iteration budget.
1072+
1073+
Expected:
1074+
Acceptance rate within target_acc ± tol = $(target_acc) ± $(tol), reached within
1075+
max_iter = $(max_iter) `sample()` calls.
1076+
1077+
Loop context:
1078+
n_evals = $(n_evals) (max_iter = $(max_iter))
1079+
last stepsize tried = $(round(stepsize; sigdigits = 3))
1080+
1081+
Suggestion:
1082+
Increase `max_iter`, or widen `tol` if a looser acceptance-rate window is acceptable.
1083+
"""))
1084+
end
1085+
1086+
@noinline function _throw_max_expansions_exceeded(direction::Symbol, n_expand, max_expansions, target_acc, tol, bound_stepsize)
1087+
above = direction == :above
1088+
action = above ? "doublings" : "halvings"
1089+
stayed = above ? "stayed above" : "stayed below"
1090+
monotonic_note =
1091+
above ? "rather than simply needing a larger stepsize" :
1092+
"(e.g. a numerically noisy log-density at very small stepsizes) rather than simply needing a smaller stepsize"
1093+
throw(ArgumentError("""
1094+
optimize_stepsize: acceptance rate $(stayed) target $(target_acc) ± $(tol) after $(max_expansions) $(action) without crossing it.
1095+
1096+
Expected:
1097+
Acceptance rate to cross target_acc ± tol = $(target_acc) ± $(tol) within
1098+
max_expansions = $(max_expansions) $(action).
1099+
1100+
Loop context:
1101+
n_expand = $(n_expand) (max_expansions = $(max_expansions))
1102+
stepsize reached = $(round(bound_stepsize; sigdigits = 3))
1103+
1104+
Suggestion:
1105+
Increase `max_expansions` to search farther. Otherwise, the acceptance-vs-stepsize
1106+
relationship may not be monotonic in this range, $(monotonic_note).
1107+
"""))
1108+
end
1109+
10611110

10621111

10631112

test/MarkovChainMonteCarlo/runtests.jl

Lines changed: 99 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@ using GaussianProcesses
55
using Test
66
using AdvancedMH
77
using AbstractMCMC
8+
using MCMCChains
89

910
using CalibrateEmulateSample.EnsembleKalmanProcesses
1011
using CalibrateEmulateSample.MarkovChainMonteCarlo
@@ -693,6 +694,90 @@ end
693694
end
694695
end
695696

697+
@testset "optimize_stepsize: branch coverage" begin
698+
# optimize_stepsize's bracket-and-bisect search has several distinct code paths
699+
mcmc = MCMCWrapper(RWMHSampling(), mcmc_params[:obs_sample], prior, em_1; init_params = vec(collect(mcmc_params[:init_params])))
700+
701+
@testset "returns immediately if the initial stepsize is already within tolerance" begin
702+
step = optimize_stepsize(Random.MersenneTwister(1), mcmc; init_stepsize = 0.125, N = 800, target_acc = 0.25)
703+
@test isapprox(step, 0.125; atol = 1e-8) # no expansion/bisection needed at all
704+
end
705+
706+
@testset "phase 1, 'acceptance too high' branch (expands stepsize upward)" begin
707+
step = optimize_stepsize(Random.MersenneTwister(2), mcmc; init_stepsize = 0.001, N = 800, target_acc = 0.25)
708+
@test 0.15 <= accept_ratio(MCMC.sample(Random.MersenneTwister(2), mcmc, 800; stepsize = step)) <= 0.35
709+
end
710+
711+
@testset "phase 1, 'acceptance too low' branch (expands stepsize downward)" begin
712+
step = optimize_stepsize(Random.MersenneTwister(3), mcmc; init_stepsize = 4.0, N = 800, target_acc = 0.25)
713+
@test 0.15 <= accept_ratio(MCMC.sample(Random.MersenneTwister(3), mcmc, 800; stepsize = step)) <= 0.35
714+
end
715+
716+
@testset "phase 2: bisection is actually reached" begin
717+
718+
step = optimize_stepsize(
719+
Random.MersenneTwister(4),
720+
mcmc;
721+
init_stepsize = 0.0625,
722+
N = 1500,
723+
target_acc = 0.33,
724+
tol = 0.05,
725+
)
726+
@test 0.0625 < step < 0.125 # strictly inside the bracket -> bisection ran, not just phase 1
727+
end
728+
729+
@testset "max_expansions exceeded: 'stayed above target' (expand-up) branch" begin
730+
let thrown = @test_throws ArgumentError optimize_stepsize(
731+
Random.MersenneTwister(5),
732+
mcmc;
733+
init_stepsize = 0.25,
734+
N = 300,
735+
target_acc = -0.5,
736+
tol = 0.05,
737+
max_expansions = 3,
738+
)
739+
@test contains(thrown.value.msg, "stayed above")
740+
@test contains(thrown.value.msg, "doublings")
741+
@test contains(thrown.value.msg, "max_expansions")
742+
@test contains(thrown.value.msg, "-0.5")
743+
end
744+
end
745+
746+
@testset "max_expansions exceeded: 'stayed below target' (expand-down) branch" begin
747+
# target_acc set above the achievable range (>1), symmetric to the case above.
748+
let thrown = @test_throws ArgumentError optimize_stepsize(
749+
Random.MersenneTwister(6),
750+
mcmc;
751+
init_stepsize = 0.25,
752+
N = 300,
753+
target_acc = 1.5,
754+
tol = 0.05,
755+
max_expansions = 3,
756+
)
757+
@test contains(thrown.value.msg, "stayed below")
758+
@test contains(thrown.value.msg, "halvings")
759+
@test contains(thrown.value.msg, "max_expansions")
760+
@test contains(thrown.value.msg, "1.5")
761+
end
762+
end
763+
764+
@testset "overall max_iter budget exhausted" begin
765+
# A reachable target, but with only 1 sample() call allowed in total: the very first
766+
# (non-converging) phase-1 expansion step must hit the n_evals > max_iter guard.
767+
let thrown = @test_throws ArgumentError optimize_stepsize(
768+
Random.MersenneTwister(7),
769+
mcmc;
770+
init_stepsize = 4.0,
771+
N = 300,
772+
target_acc = 0.25,
773+
max_iter = 1,
774+
)
775+
@test contains(thrown.value.msg, "iteration budget")
776+
@test contains(thrown.value.msg, "max_iter")
777+
end
778+
end
779+
end
780+
696781
@testset "Sampler transition-density interface (pCN Hastings-term fix)" begin
697782
# Regression tests for: pCN's (and Barker's) log_transition_density must reflect the
698783
# *actual* asymmetric proposal kernel, not the symmetric additive random-walk kernel.
@@ -819,7 +904,20 @@ end
819904
# A hypothetical new sampler that forgets to implement log_transition_density must
820905
# fail loudly (ArgumentError) rather than silently falling back to a symmetric,
821906
# possibly-wrong correction — this is what the shared interface guards against.
822-
@test_throws ArgumentError MCMC.log_transition_density(_DummyMHSampler(), dummy_model, a, b)
907+
let thrown = @test_throws ArgumentError MCMC.log_transition_density(_DummyMHSampler(), dummy_model, a, b)
908+
@test contains(thrown.value.msg, "_DummyMHSampler")
909+
@test contains(thrown.value.msg, "log_transition_density")
910+
end
911+
end
912+
end
913+
914+
@testset "accept_ratio errors on a chain missing the :accepted internal" begin
915+
# A Chains object not produced by this module's sample()/optimize_stepsize() (e.g. hand-built,
916+
# or with internals filtered out) must fail loudly rather than silently mis-reporting.
917+
bad_chain = MCMCChains.Chains(rand(5, 3, 1), [:a, :b, :other_internal], (parameters = [:a, :b], internals = [:other_internal]))
918+
let thrown = @test_throws ArgumentError accept_ratio(bad_chain)
919+
@test contains(thrown.value.msg, "accept_ratio")
920+
@test contains(thrown.value.msg, "other_internal")
823921
end
824922
end
825923
end

0 commit comments

Comments
 (0)