From b1cf1bdef79d03edfd7aa2248294562a681e9380 Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 22 Oct 2026 16:08:11 +1100 Subject: [PATCH 1/6] Keep mcmc()'s one_by_one advice when a numerical error stops sampling before its first draw, rather than "subscript out of bounds" (#577) - R/sampler_class.R: trace_values() returns early when the trace has no draws, so the clean-up on exit no longer replaces the error that stopped sampling - R/sampler_class.R: abort_numerical_error() raises its error with call = NULL, so the header reads "Error:" rather than naming an internal method - tests/testthat/test_inference.R: a singular solve() stops mcmc() with the one_by_one advice; the one_by_one test also checks it keeps all 40 draws; the numerical-error tests share one model helper - tests/testthat/helpers.R: singular_at_one_model(), whose solve() errors whenever round(r) is 1 or -1 - tests/testthat/_snaps/inference.md: snapshot of the one_by_one advice - tests/testthat/test_inference_numerical_issues.R: removed, as it was all commented out and its test now lives in test_inference.R - NEWS.md: bullet for #577 Co-Authored-By: Claude Opus 5.5 --- NEWS.md | 1 + R/sampler_class.R | 9 ++++- tests/testthat/_snaps/inference.md | 13 ++++++++ tests/testthat/helpers.R | 10 ++++++ tests/testthat/test_inference.R | 33 +++++++++++++------ .../test_inference_numerical_issues.R | 30 ----------------- 6 files changed, 55 insertions(+), 41 deletions(-) delete mode 100644 tests/testthat/test_inference_numerical_issues.R diff --git a/NEWS.md b/NEWS.md index af5085516..5f672ccd6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -65,6 +65,7 @@ * `%*%` now multiplies a data matrix used directly in the product (an R matrix, or a greta array from `as_data()` or `as_data_mutable()`) once, instead of once per chain, so models such as `X %*% beta` run faster in `mcmc()` with several chains (#858). * `[<-` on a greta array fills a matrix or array replacement column by column, as R does (#844). * `mcmc()` and `extra_samples()` keep `n_samples %/% thin` draws, regardless of values of `pb_update` and `one_by_one`. A burst shorter than `thin` used to error, as in a final partial burst, `one_by_one = TRUE` with `thin` above 1, or `extra_samples()` with `pb_update` at or below `thin`. They run all `n_samples` iterations even when `thin` does not divide it, `one_by_one = TRUE` rejects only the proposal that caused a numerical error, and a `thin` larger than `n_samples` is an informative error (#241, #318, #567, #609). +* `mcmc()` gives its error suggesting `one_by_one = TRUE` when a numerical error stops sampling before its first draw, which with `verbose = FALSE` is any numerical error, rather than "subscript out of bounds" (#577). * `mcmc()`'s progress bar reaches its total when `pb_update` does not divide the number of iterations, and fills in step with the iterations run; a bar could read 1010/1010 while still empty (#318). * `opt()`'s `convergence` comes from the optimiser: `bfgs()` and `nelder_mead()` report it, and a Keras optimiser such as `adam()` has converged when its last step changed the objective by no more than `tolerance`, rather than being guessed from whether it stopped before `max_iterations - 1` (#569). diff --git a/R/sampler_class.R b/R/sampler_class.R index 56fae41ea..87a8dd9cc 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -516,13 +516,20 @@ sampler <- R6Class( "This will slow down the sampler slightly.", "The error encountered can be recovered and viewed with:", "{.code greta_notes_tf_num_error()}" - ) + ), + call = NULL ) }, # convert traced free state to the traced values, accounting for # chain dimension trace_values = function(trace_batch_size) { + # on.exit() calls this after an error too; tracing no rows would hide + # that error + no_draws <- nrow(self$traced_free_state[[1]]) == 0 + if (no_draws) { + return() + } self$traced_values <- lapply( self$traced_free_state, self$model$dag$trace_values, diff --git a/tests/testthat/_snaps/inference.md b/tests/testthat/_snaps/inference.md index 55b51d36d..b0da9a279 100644 --- a/tests/testthat/_snaps/inference.md +++ b/tests/testthat/_snaps/inference.md @@ -229,6 +229,19 @@ ! `Lmin` must not be larger than `Lmax` x `Lmin` is 10 and `Lmax` is 5 +# a numerical error stops mcmc() with the one_by_one advice + + Code + mcmc(m, sampler = rwmh(epsilon = 1), warmup = 0, n_samples = 40, chains = 1, + initial_values = initials(r = 0.2), verbose = FALSE) + Condition + Error: + ! TensorFlow hit a numerical problem that caused it to error + greta can handle these as bad proposals if you rerun `mcmc()` with the argument `one_by_one = TRUE`. + This will slow down the sampler slightly. + The error encountered can be recovered and viewed with: + `greta_notes_tf_num_error()` + # thin larger than n_samples is an informative error Code diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index e7ef81434..9cf378176 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -464,6 +464,16 @@ mock_mcmc <- function(n_samples = 1010) { iterate_progress_bar(pb, n_samples, rejects = 10, chains = 1) } +# a model of r whose solve() errors on the singular matrix whenever round(r) is +# 1 or -1, so mcmc() hits numerical errors. r is made by the caller, where +# initials(r = ) looks for it +singular_at_one_model <- function(r) { + singular_at_one <- diag(2) + (1 - diag(2)) * round(r) + y <- as_data(0.5) + distribution(y) <- normal(sum(solve(singular_at_one)), 1) + model(r) +} + # random lkj draws, code from the rethinking package (can't load the package # because of stan*Travis*compiler issues) rlkjcorr <- function(n, eta = 1, dimension = 2) { diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index e128628c8..367ba0422 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -705,15 +705,31 @@ test_that("a numerical error at the end of warmup still tunes", { expect_equal(sampler$parameters$epsilon, exp(log_epsilon_bar)) }) +test_that("a numerical error stops mcmc() with the one_by_one advice", { + skip_if_not(check_tf_version()) + local_greta_seed() + r <- normal(0, 1) + m <- singular_at_one_model(r) + + expect_snapshot( + error = TRUE, + mcmc( + m, + sampler = rwmh(epsilon = 1), + warmup = 0, + n_samples = 40, + chains = 1, + initial_values = initials(r = 0.2), + verbose = FALSE + ) + ) +}) + test_that("a numerical error with one_by_one repeats the draw before it", { skip_if_not(check_tf_version()) local_greta_seed() - # solve() errors on the singular matrix whenever round(r) is 1 or -1 r <- normal(0, 1) - singular_at_one <- diag(2) + (1 - diag(2)) * round(r) - y <- as_data(0.5) - distribution(y) <- normal(sum(solve(singular_at_one)), 1) - m <- model(r) + m <- singular_at_one_model(r) draws <- mcmc( m, @@ -727,6 +743,7 @@ test_that("a numerical error with one_by_one repeats the draw before it", { ) sampler <- get_model_info(draws)$samplers[[1]] expect_gt(sampler$numerical_rejections, 0) + expect_equal(coda::niter(draws), 40) chain <- as.numeric(draws[[1]]) repeated <- which(duplicated(chain)) @@ -736,12 +753,8 @@ test_that("a numerical error with one_by_one repeats the draw before it", { test_that("numerical errors in warmup with one_by_one still tune", { skip_if_not(check_tf_version()) local_greta_seed() - # solve() errors on the singular matrix whenever round(r) is 1 or -1 r <- normal(0, 1) - singular_at_one <- diag(2) + (1 - diag(2)) * round(r) - y <- as_data(0.5) - distribution(y) <- normal(sum(solve(singular_at_one)), 1) - m <- model(r) + m <- singular_at_one_model(r) warmup <- 30 chains <- 2 diff --git a/tests/testthat/test_inference_numerical_issues.R b/tests/testthat/test_inference_numerical_issues.R deleted file mode 100644 index fdcd5e743..000000000 --- a/tests/testthat/test_inference_numerical_issues.R +++ /dev/null @@ -1,30 +0,0 @@ -# set.seed(2020 - 02 - 11) -# -# test_that("numerical issues are handled in mcmc", { -# skip_if_not(check_tf_version()) -# -# # this should have a cholesky decomposition problem at some point -# alpha <- normal(0, 1) -# # x <- matrix(rnorm(6), 3, 2) -# x <- matrix(rnorm(12), 4, 3) -# # y <- t(rnorm(3)) -# y <- t(rnorm(4)) -# z <- alpha * x -# sigma <- z %*% t(z) -# # distribution(y) <- multivariate_normal(zeros(1, 3), sigma) -# distribution(y) <- multivariate_normal(zeros(1, 4), sigma) -# m <- model(alpha) -# -# # running with bursts should error informatively -# expect_snapshot( -# error = TRUE, -# draws <- mcmc(m, verbose = FALSE) -# ) -# -# # setting one_by_one = TRUE should handle those errors as bad samples -# expect_no_error(draws <- mcmc(m, -# warmup = 100, n_samples = 10, -# one_by_one = TRUE, -# verbose = FALSE -# )) -# }) From fa7c6180756d4e213e730cb101f898f7d16ebf2e Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 22 Oct 2026 16:20:04 +1100 Subject: [PATCH 2/6] Say why trace_values() returns early and why the numerical error names no call - R/sampler_class.R: trace_values()'s comment says dag$trace_values() fails on a trace with no rows, and that on.exit() failing would replace the error being unwound - R/sampler_class.R: abort_numerical_error()'s comment says its error names no call because the methods it is raised from mean nothing to the user Co-Authored-By: Claude Opus 5.5 --- R/sampler_class.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/sampler_class.R b/R/sampler_class.R index 87a8dd9cc..6bad1d5e1 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -505,7 +505,8 @@ sampler <- R6Class( # In a call of more than one iteration, the sampler cannot tell which # proposal failed, and the chain would not be valid if it carried on, so - # it stops and says how to run one iteration per call + # it stops and says how to run one iteration per call. The error names no + # call, since the methods it is raised from mean nothing to the user abort_numerical_error = function(error) { greta_stash$tf_num_error <- error cli::cli_abort( @@ -524,8 +525,8 @@ sampler <- R6Class( # convert traced free state to the traced values, accounting for # chain dimension trace_values = function(trace_batch_size) { - # on.exit() calls this after an error too; tracing no rows would hide - # that error + # dag$trace_values() fails on a trace with no rows, and on.exit() calls + # this after an error too, where failing would replace that error no_draws <- nrow(self$traced_free_state[[1]]) == 0 if (no_draws) { return() From 0723026da9d17a4ad2c9b3df8ddf7bb28198bf8d Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 22 Oct 2026 16:20:04 +1100 Subject: [PATCH 3/6] Make stashed_samples() return no draws when several chains stop before their first, rather than a dimnames error (#577) - R/inference_class.R: n_traced counts the columns of the traced values rather than every element, which multiplied by the number of chains and made the empty trace too wide for its names - tests/testthat/test_inference.R: two chains stopped by a numerical error before their first draw give the one_by_one advice, and stashed_samples() then has no draws - NEWS.md: bullet for stashed_samples() Co-Authored-By: Claude Opus 5.5 --- NEWS.md | 1 + R/inference_class.R | 2 +- tests/testthat/test_inference.R | 21 +++++++++++++++++++++ 3 files changed, 23 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 5f672ccd6..895a51ba4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -68,6 +68,7 @@ * `mcmc()` gives its error suggesting `one_by_one = TRUE` when a numerical error stops sampling before its first draw, which with `verbose = FALSE` is any numerical error, rather than "subscript out of bounds" (#577). * `mcmc()`'s progress bar reaches its total when `pb_update` does not divide the number of iterations, and fills in step with the iterations run; a bar could read 1010/1010 while still empty (#318). * `opt()`'s `convergence` comes from the optimiser: `bfgs()` and `nelder_mead()` report it, and a Keras optimiser such as `adam()` has converged when its last step changed the objective by no more than `tolerance`, rather than being guessed from whether it stopped before `max_iterations - 1` (#569). +* `stashed_samples()` returns no draws when `mcmc()` with more than one chain stopped before its first draw, rather than erroring with "length of 'dimnames' [2] not equal to array extent" (#577). ## Deprecated and defunct diff --git a/R/inference_class.R b/R/inference_class.R index ad7d89b74..7253e120c 100644 --- a/R/inference_class.R +++ b/R/inference_class.R @@ -43,7 +43,7 @@ inference <- R6Class( self$set_initial_values(initial_values) # maybe inefficient- potentially speed up? - self$n_traced <- length(model$dag$trace_values(self$free_state)) + self$n_traced <- ncol(model$dag$trace_values(self$free_state)) self$seed <- seed }, diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index 367ba0422..e26630b41 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -725,6 +725,27 @@ test_that("a numerical error stops mcmc() with the one_by_one advice", { ) }) +test_that("stashed_samples() has no draws when several chains stop before one", { + skip_if_not(check_tf_version()) + local_greta_seed() + r <- normal(0, 1) + m <- singular_at_one_model(r) + + expect_error( + mcmc( + m, + sampler = rwmh(epsilon = 1), + warmup = 0, + n_samples = 40, + chains = 2, + initial_values = list(initials(r = 0.2), initials(r = -0.2)), + verbose = FALSE + ), + "one_by_one = TRUE" + ) + expect_identical(stashed_samples(), list(NULL)) +}) + test_that("a numerical error with one_by_one repeats the draw before it", { skip_if_not(check_tf_version()) local_greta_seed() From 48976ff5fb93d93a2a5773daff378b1e822fa4fa Mon Sep 17 00:00:00 2001 From: njtierney Date: Thu, 22 Oct 2026 16:21:44 +1100 Subject: [PATCH 4/6] Put the reason abort_numerical_error() names no call next to call = NULL - R/sampler_class.R: the comment explaining call = NULL moves from above the method to the line above the argument Co-Authored-By: Claude Opus 5.5 --- R/sampler_class.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/sampler_class.R b/R/sampler_class.R index 6bad1d5e1..fa6db36b8 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -505,8 +505,7 @@ sampler <- R6Class( # In a call of more than one iteration, the sampler cannot tell which # proposal failed, and the chain would not be valid if it carried on, so - # it stops and says how to run one iteration per call. The error names no - # call, since the methods it is raised from mean nothing to the user + # it stops and says how to run one iteration per call abort_numerical_error = function(error) { greta_stash$tf_num_error <- error cli::cli_abort( @@ -518,6 +517,7 @@ sampler <- R6Class( "The error encountered can be recovered and viewed with:", "{.code greta_notes_tf_num_error()}" ), + # the methods this is raised from mean nothing to the user call = NULL ) }, From 126ecd8cc4bd1fc7782263544cccdddc07827641 Mon Sep 17 00:00:00 2001 From: njtierney Date: Fri, 23 Oct 2026 11:53:28 +1100 Subject: [PATCH 5/6] Put back the samplers and TensorFlow error the numerical-error tests stash, so no later test runs on them - tests/testthat/helpers.R: local_sampler_stash() restores greta_stash's samplers and tf_num_error when the test exits - tests/testthat/test_inference.R: the two tests whose mcmc() aborts call it Co-Authored-By: Claude Opus 5.5 --- tests/testthat/helpers.R | 14 ++++++++++++++ tests/testthat/test_inference.R | 2 ++ 2 files changed, 16 insertions(+) diff --git a/tests/testthat/helpers.R b/tests/testthat/helpers.R index 9cf378176..28c2c03f0 100644 --- a/tests/testthat/helpers.R +++ b/tests/testthat/helpers.R @@ -474,6 +474,20 @@ singular_at_one_model <- function(r) { model(r) } +# put back the samplers and TensorFlow error that an aborted mcmc() stashes, +# when the caller exits +local_sampler_stash <- function(.local_envir = parent.frame()) { + fields <- c("samplers", "tf_num_error") + old <- mget(intersect(fields, ls(greta_stash)), envir = greta_stash) + withr::defer( + { + rm(list = intersect(fields, ls(greta_stash)), envir = greta_stash) + list2env(old, envir = greta_stash) + }, + envir = .local_envir + ) +} + # random lkj draws, code from the rethinking package (can't load the package # because of stan*Travis*compiler issues) rlkjcorr <- function(n, eta = 1, dimension = 2) { diff --git a/tests/testthat/test_inference.R b/tests/testthat/test_inference.R index e26630b41..cd2475455 100644 --- a/tests/testthat/test_inference.R +++ b/tests/testthat/test_inference.R @@ -708,6 +708,7 @@ test_that("a numerical error at the end of warmup still tunes", { test_that("a numerical error stops mcmc() with the one_by_one advice", { skip_if_not(check_tf_version()) local_greta_seed() + local_sampler_stash() r <- normal(0, 1) m <- singular_at_one_model(r) @@ -728,6 +729,7 @@ test_that("a numerical error stops mcmc() with the one_by_one advice", { test_that("stashed_samples() has no draws when several chains stop before one", { skip_if_not(check_tf_version()) local_greta_seed() + local_sampler_stash() r <- normal(0, 1) m <- singular_at_one_model(r) From db0143c671c531a5c4b67a741b309fef097e148b Mon Sep 17 00:00:00 2001 From: njtierney Date: Fri, 23 Oct 2026 16:19:52 +1100 Subject: [PATCH 6/6] Return NULL explicitly from trace_values() when there are no draws, as greta's other early returns do - R/sampler_class.R: return(NULL) rather than a bare return() Co-Authored-By: Claude Opus 5.5 --- R/sampler_class.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/sampler_class.R b/R/sampler_class.R index fa6db36b8..efa7dedfc 100644 --- a/R/sampler_class.R +++ b/R/sampler_class.R @@ -529,7 +529,7 @@ sampler <- R6Class( # this after an error too, where failing would replace that error no_draws <- nrow(self$traced_free_state[[1]]) == 0 if (no_draws) { - return() + return(NULL) } self$traced_values <- lapply( self$traced_free_state,