diff --git a/NEWS.md b/NEWS.md index 608dc615..bcc1e16c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,6 @@ # loo (development version) +* Handle errors of the user functions in `loo_moment_match()` by @avehtari in #410 * Allow different numbers of draws in `loo_model_weights()` by @avehtari in #409 * Improve numerical stability in `loo()`, `psis()`, model weighting, subsampling, and moment matching in #395 diff --git a/R/loo_moment_matching.R b/R/loo_moment_matching.R index 91d17ee2..08d6ced6 100644 --- a/R/loo_moment_matching.R +++ b/R/loo_moment_matching.R @@ -139,8 +139,38 @@ loo_moment_match.default <- function(x, loo, post_draws, log_lik_i, } # update results + failed <- integer(0) + # observations where log_lik_i() failed, a subset of failed + log_lik_failed <- integer(0) + stopped <- integer(0) + split_failed <- integer(0) + error_msg <- NULL for (ii in seq_along(I)) { - i <- mm_list[[ii]]$i + i <- I[ii] + # parallel::mclapply() returns NULL if a worker dies, for example when it + # runs out of memory, and a "try-error" string if a worker stops with an + # error. Errors of the user functions are caught in loo_moment_match_i(), + # so a worker error is an error in loo itself, as it is with cores = 1. + if (is.null(mm_list[[ii]])) { + mm_list[[ii]] <- list(i = i, failed = TRUE, + error = "a parallel worker did not return a result") + } else if (inherits(mm_list[[ii]], "try-error")) { + stop(try_error_message(mm_list[[ii]]), call. = FALSE) + } + if (!is.null(mm_list[[ii]]$error)) { + error_msg <- c(error_msg, mm_list[[ii]]$error) + if (isTRUE(mm_list[[ii]]$failed)) { + failed <- c(failed, i) + if (isTRUE(mm_list[[ii]]$log_lik_failed)) { + log_lik_failed <- c(log_lik_failed, i) + } + next + } + stopped <- c(stopped, i) + } + if (isTRUE(mm_list[[ii]]$split_failed)) { + split_failed <- c(split_failed, i) + } loo$pointwise[i, "elpd_loo"] <- mm_list[[ii]]$elpd_loo_i loo$pointwise[i, "p_loo"] <- mm_list[[ii]]$p_loo loo$pointwise[i, "mcse_elpd_loo"] <- mm_list[[ii]]$mcse_elpd_loo @@ -154,6 +184,14 @@ loo_moment_match.default <- function(x, loo, post_draws, log_lik_i, loo$psis_object$log_weights[, i] <- mm_list[[ii]]$lwi } } + # log_lik_i() is evaluated at the original draws, so an error for every + # observation is probably an error in the function, for example wrong + # arguments. An error at the first transformation can come from extreme + # draws, and keeps the PSIS results with a warning. + if (length(I) > 0 && length(log_lik_failed) == length(I)) { + stop("'log_lik_i' failed for all observations. The first error ", + "message:\n", error_msg[1], call. = FALSE) + } if (!is.null(loo$psis_object)) { attr(loo$psis_object, "norm_const_log") <- matrixStats::colLogSumExps(loo$psis_object$log_weights) loo$psis_object$diagnostics <- loo$diagnostics @@ -173,11 +211,14 @@ loo_moment_match.default <- function(x, loo, post_draws, log_lik_i, loo$se_p_loo <- loo$estimates["p_loo","SE"] loo$se_looic <- loo$estimates["looic","SE"] + throw_moment_match_error_warning(failed, stopped, error_msg) # Warn if some Pareto ks are still high throw_pareto_warnings(loo$diagnostics$pareto_k, k_threshold) # if we don't split, accuracy may be compromised if (!split) { throw_large_kf_warning(kfs, k_threshold) + } else if (length(split_failed) > 0) { + throw_large_kf_warning(kfs[split_failed], k_threshold, split_failed = TRUE) } loo @@ -252,7 +293,13 @@ loo_moment_match_i <- function(i, uparsi <- upars ki <- k kfi <- 0 - log_liki <- log_lik_i(x, i, ...) + log_liki <- tryCatch(log_lik_i(x, i, ...), error = function(e) e) + if (inherits(log_liki, "error")) { + # the user function failed for this observation, which keeps its PSIS + # results + return(list(i = i, failed = TRUE, log_lik_failed = TRUE, + error = conditionMessage(log_liki))) + } S_per_chain <- NROW(log_liki) N_chains <- NCOL(log_liki) dim(log_liki) <- c(S_per_chain, N_chains, 1) @@ -274,6 +321,9 @@ loo_moment_match_i <- function(i, mapping = diag(npars) ) + # message of an error that stopped moment matching for this observation + error_msg <- NULL + # try several transformations one by one # if one does not work, do not apply it and try another one # to accept the transformation, Pareto k needs to improve @@ -295,11 +345,13 @@ loo_moment_match_i <- function(i, r_eff_i = r_eff_i, cores = 1, is_method = is_method, - ...) + ...), + silent = TRUE ) if (inherits(quantities_i, "try-error")) { - # Stan log prob caused an exception probably due to under- or - # overflow of parameters to invalid values + # log_prob_upars() or log_lik_i_upars() failed, for example because the + # transformed draws are extreme + error_msg <- try_error_message(quantities_i) break } if (quantities_i$ki < ki) { @@ -325,11 +377,13 @@ loo_moment_match_i <- function(i, r_eff_i = r_eff_i, cores = 1, is_method = is_method, - ...) + ...), + silent = TRUE ) if (inherits(quantities_i, "try-error")) { - # Stan log prob caused an exception probably due to under- or - # overflow of parameters to invalid values + # log_prob_upars() or log_lik_i_upars() failed, for example because the + # transformed draws are extreme + error_msg <- try_error_message(quantities_i) break } if (quantities_i$ki < ki) { @@ -356,12 +410,14 @@ loo_moment_match_i <- function(i, r_eff_i = r_eff_i, cores = 1, is_method = is_method, - ...) + ...), + silent = TRUE ) if (inherits(quantities_i, "try-error")) { - # Stan log prob caused an exception probably due to under- or - # overflow of parameters to invalid values + # log_prob_upars() or log_lik_i_upars() failed, for example because the + # transformed draws are extreme + error_msg <- try_error_message(quantities_i) break } if (quantities_i$ki < ki) { @@ -386,19 +442,28 @@ loo_moment_match_i <- function(i, # if we don't do split transform, or # if no transformations were successful # stop and collect values + split_done <- FALSE + split_failed <- FALSE if (split && (iterind > 1)) { # compute split transformation - split_obj <- loo_moment_match_split( + split_obj <- try(loo_moment_match_split( x, upars, cov, total$shift, total$scaling, total$mapping, i, log_prob_upars = log_prob_upars, log_lik_i_upars = log_lik_i_upars, cores = 1, r_eff_i = r_eff_i, is_method = is_method, ... - ) - log_liki <- split_obj$log_liki - lwi <- split_obj$lwi - lwfi <- split_obj$lwfi - r_eff_i <- split_obj$r_eff_i + ), silent = TRUE) + if (inherits(split_obj, "try-error")) { + # use the result without the split transformation + error_msg <- try_error_message(split_obj) + split_failed <- TRUE + } else { + log_liki <- split_obj$log_liki + lwi <- split_obj$lwi + lwfi <- split_obj$lwfi + r_eff_i <- split_obj$r_eff_i + split_done <- TRUE + } } - else { + if (!split_done) { dim(log_liki) <- c(S_per_chain, N_chains, 1) r_eff_i <- loo::relative_eff(exp(log_liki), cores = 1) dim(log_liki) <- NULL @@ -420,7 +485,12 @@ loo_moment_match_i <- function(i, n_eff = min(1.0 / sum(exp(2 * lwi)), 1.0 / sum(exp(2 * lwfi))) * r_eff_i, lwi = lwi, - i = i) + i = i, + error = error_msg, + # an error before any transformation was accepted leaves the PSIS + # results + failed = !is.null(error_msg) && iterind == 1, + split_failed = split_failed) } @@ -635,6 +705,46 @@ shift_and_cov <- function(x, upars, lwi, ...) { } +#' Message of a "try-error" object without the call and the newline +#' @noRd +try_error_message <- function(x) { + cond <- attr(x, "condition") + if (!is.null(cond)) conditionMessage(cond) else trimws(as.character(x)) +} + +#' Warning message if a function call failed during moment matching +#' @noRd +#' @param failed Observations for which moment matching failed completely. +#' @param stopped Observations for which moment matching stopped early. +#' @param error_msg Error messages, the first one is shown. +throw_moment_match_error_warning <- function(failed, stopped, error_msg) { + if (length(failed) == 0 && length(stopped) == 0) { + return(invisible(NULL)) + } + msg <- NULL + if (length(stopped) > 0) { + msg <- c(msg, paste0( + "Moment matching stopped early for observation(s) ", + paste(stopped, collapse = ", "), ". The results for these ", + "observations use the last transformation that succeeded.\n" + )) + } + if (length(failed) > 0) { + msg <- c(msg, paste0( + "Moment matching failed for observation(s) ", + paste(failed, collapse = ", "), ". The results for these ", + "observations are the PSIS results without moment matching.\n" + )) + } + warning( + msg, + "The cause was an error in 'log_prob_upars', 'log_lik_i_upars' or ", + "'log_lik_i', or a parallel worker that did not return a result. The ", + "first error message:\n", error_msg[1], + call. = FALSE + ) +} + #' Warning message if max_iters is reached #' @noRd throw_moment_match_max_iters_warning <- function() { @@ -649,11 +759,19 @@ throw_moment_match_max_iters_warning <- function() { #' Warning message if not using split transformation and accuracy is #' compromised #' @noRd -throw_large_kf_warning <- function(kf, k_threshold) { +throw_large_kf_warning <- function(kf, k_threshold, split_failed = FALSE) { if (any(kf > k_threshold)) { + if (split_failed) { + advice <- paste0( + "The split transformation failed for these observations, ", + "see the warning about the error." + ) + } else { + advice <- "Setting the argument 'split' to 'TRUE' will likely improve accuracy." + } warning( "The accuracy of self-normalized importance sampling may be bad.\n", - "Setting the argument 'split' to 'TRUE' will likely improve accuracy.", + advice, call. = FALSE ) } diff --git a/tests/testthat/test_loo_moment_matching.R b/tests/testthat/test_loo_moment_matching.R index c69f3401..1ea822ff 100644 --- a/tests/testthat/test_loo_moment_matching.R +++ b/tests/testthat/test_loo_moment_matching.R @@ -734,3 +734,172 @@ test_that("passing arguments works", { passed_arg = TRUE )) }) + +# log_lik_i_upars() that fails in the calls listed in `fail_calls`, as a model +# can fail for the extreme draws of a transformation +failing_log_lik_i_upars <- function(fail_calls) { + calls <- 0 + function(x, upars, i, ...) { + calls <<- calls + 1 + if (calls %in% fail_calls) { + stop("overflow in log_lik_i_upars") + } + log_lik_i_upars_test(x, upars, i) + } +} + +mm_with <- function(loo_x, log_lik_i = log_lik_i_test, + log_lik_i_upars = log_lik_i_upars_test, split = TRUE, + k_threshold = 0.7, cores = 1) { + loo_moment_match( + x, loo = loo_x, post_draws = post_draws_test, log_lik_i = log_lik_i, + unconstrain_pars = unconstrain_pars_test, + log_prob_upars = log_prob_upars_test, log_lik_i_upars = log_lik_i_upars, + k_threshold = k_threshold, split = split, cores = cores + ) +} + +test_that("an error in the split step keeps the result without the split", { + loo_manual <- suppressWarnings(loo(loglik)) + # observation 1 needs two transformations; the third call is the split step + warnings <- character(0) + mm <- withCallingHandlers( + mm_with(loo_manual, log_lik_i_upars = failing_log_lik_i_upars(3)), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + expect_match(warnings, "stopped early for observation(s) 1", fixed = TRUE, + all = FALSE) + expect_match(warnings, "The split transformation failed", all = FALSE) + mm_no_split <- suppressWarnings(mm_with(loo_manual, split = FALSE)) + expect_equal(mm$pointwise, mm_no_split$pointwise) + expect_equal(mm$diagnostics, mm_no_split$diagnostics) +}) + +test_that("the accuracy warning names a failed split transformation", { + expect_warning(throw_large_kf_warning(0.8, 0.7), "Setting the argument 'split'") + expect_warning( + throw_large_kf_warning(0.8, 0.7, split_failed = TRUE), + "The split transformation failed" + ) + expect_no_warning(throw_large_kf_warning(0.5, 0.7, split_failed = TRUE)) +}) + +test_that("an error in the first transformation keeps the PSIS results", { + loo_manual <- suppressWarnings(loo(loglik)) + # k_threshold = 0 moment-matches all observations; the first call of + # log_lik_i_upars is the first transformation of observation 1 + warnings <- character(0) + mm <- withCallingHandlers( + mm_with(loo_manual, log_lik_i_upars = failing_log_lik_i_upars(1), + k_threshold = 0), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + # no transformation was accepted, so the observation counts as failed + expect_match(warnings, "failed for observation(s) 1.", fixed = TRUE, + all = FALSE) + expect_match(warnings, "overflow in log_lik_i_upars", all = FALSE) + expect_equal(mm$pointwise[1, ], loo_manual$pointwise[1, ]) + expect_equal(mm$diagnostics$pareto_k[1], loo_manual$diagnostics$pareto_k[1]) +}) + +test_that("first-transformation errors everywhere keep the PSIS results", { + loo_manual <- suppressWarnings(loo(loglik)) + # only observation 1 is above k_threshold = 0.7, and its first + # transformation fails, as with an overflow in a *_rng function (#274) + warnings <- character(0) + mm <- withCallingHandlers( + mm_with(loo_manual, log_lik_i_upars = failing_log_lik_i_upars(1)), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + expect_match(warnings, "failed for observation(s) 1.", fixed = TRUE, + all = FALSE) + expect_match(warnings, "overflow in log_lik_i_upars", all = FALSE) + expect_equal(mm$pointwise, loo_manual$pointwise) + expect_equal(mm$diagnostics, loo_manual$diagnostics) +}) + +test_that("moment matching stops if log_lik_i fails for all observations", { + loo_manual <- suppressWarnings(loo(loglik)) + # a log_lik_i with the wrong arguments fails for every observation + wrong_log_lik_i <- function(x) x + expect_error( + suppressWarnings(mm_with(loo_manual, log_lik_i = wrong_log_lik_i, + k_threshold = 0)), + "'log_lik_i' failed for all observations" + ) +}) + +test_that("errors in loo itself are not reported as failed observations", { + loo_manual <- suppressWarnings(loo(loglik)) + # log_lik_i returns without an error, and loo's own code then fails + non_numeric_log_lik_i <- function(x, i, ...) { + if (i == 1) "not numeric" else log_lik_i_test(x, i) + } + cores_list <- if (os_is_windows()) 1 else c(1, 2) + for (cores in cores_list) { + expect_error( + # with cores = 2, mclapply() also warns that a worker had an error + suppressWarnings(mm_with(loo_manual, log_lik_i = non_numeric_log_lik_i, + k_threshold = 0, cores = cores)), + "non-numeric argument to mathematical function" + ) + } +}) + +test_that("a parallel worker that dies gives failed observations", { + skip_on_cran() + skip_on_os("windows") # the worker is a forked process + loo_manual <- suppressWarnings(loo(loglik)) + parent <- Sys.getpid() + dying_log_lik_i <- function(x, i, ...) { + if (i == 1 && Sys.getpid() != parent) { + tools::pskill(Sys.getpid(), tools::SIGKILL) + } + log_lik_i_test(x, i) + } + warnings <- character(0) + mm <- withCallingHandlers( + mm_with(loo_manual, log_lik_i = dying_log_lik_i, k_threshold = 0, + cores = 2), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + expect_match(warnings, "a parallel worker did not return a result", + all = FALSE) + expect_equal(mm$pointwise[1, ], loo_manual$pointwise[1, ]) +}) + +test_that("an error for one observation keeps its PSIS results", { + skip_on_cran() + loo_manual <- suppressWarnings(loo(loglik)) + failing_log_lik_i <- function(x, i, ...) { + if (i == 1) stop("log_lik_i failed") + log_lik_i_test(x, i) + } + mm_ok <- suppressWarnings(mm_with(loo_manual, k_threshold = 0)) + # On Windows the workers are PSOCK processes that load the installed loo, so + # with devtools::test() they would not run this version of the code + cores_list <- if (os_is_windows()) 1 else c(1, 2) + for (cores in cores_list) { + # observation 1 keeps its high Pareto k, which gives a second warning + suppressWarnings(expect_warning( + mm <- mm_with(loo_manual, log_lik_i = failing_log_lik_i, + k_threshold = 0, cores = cores), + "failed for observation(s) 1.", + fixed = TRUE + )) + expect_equal(mm$pointwise[1, ], loo_manual$pointwise[1, ]) + expect_equal(mm$pointwise[-1, ], mm_ok$pointwise[-1, ]) + } +})