From e7d45e7d4d05023d795c65a5a546fb606deb5214 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 12 Sep 2026 22:39:47 +0300 Subject: [PATCH 01/17] fix: stabilize MCSE calculation --- R/loo.R | 27 +++++++++++++++++---------- tests/testthat/test_loo_and_waic.R | 12 ++++++++++++ 2 files changed, 29 insertions(+), 10 deletions(-) diff --git a/R/loo.R b/R/loo.R index 10b1bdc7..3833e04d 100644 --- a/R/loo.R +++ b/R/loo.R @@ -488,24 +488,31 @@ importance_sampling_loo_object <- function(pointwise, diagnostics, dims, #' @return Vector of standard error estimates. #' mcse_elpd <- function(ll, lw, E_elpd, r_eff, n_samples = NULL) { - lik <- exp(ll) - w2 <- exp(lw)^2 - E_epd <- exp(E_elpd) if (length(r_eff) == 1 && !is.null(ncol(ll))) { r_eff <- rep(r_eff, ncol(ll)) } var_elpd <- vapply( - seq_len(ncol(w2)), + seq_len(ncol(lw)), FUN.VALUE = numeric(1), FUN = function(i) { - # Variance in linear scale - # Equation (6) in Vehtari et al. (2024) - var_epd_i <- sum(w2[, i] * (lik[, i] - E_epd[i]) ^ 2) / r_eff[i] - # Compute variance in log scale by match the variance of a - # log-normal approximation + # Numerically stable way to compute + # 1) variance in linear scale. Equation (6) in Vehtari et al. (2024) + # 2) variance in log scale by matching the variance of a log-normal # https://en.wikipedia.org/wiki/Log-normal_distribution#Arithmetic_moments - log(1 + var_epd_i / E_epd[i]^2) + log_lik_ratio <- ll[, i] - E_elpd[i] + positive <- log_lik_ratio > 0 + log_abs_diff <- numeric(length(log_lik_ratio)) + log_abs_diff[positive] <- + log_lik_ratio[positive] + log1p(-exp(-log_lik_ratio[positive])) + log_abs_diff[!positive] <- log(-expm1(log_lik_ratio[!positive])) + log_var_epd_ratio <- + matrixStats::logSumExp(2 * lw[, i] + 2 * log_abs_diff) - log(r_eff[i]) + if (log_var_epd_ratio > 0) { + log_var_epd_ratio + log1p(exp(-log_var_epd_ratio)) + } else { + log1p(exp(log_var_epd_ratio)) + } } ) sqrt(var_elpd) diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R index 9bbd16ea..d7e50e6c 100644 --- a/tests/testthat/test_loo_and_waic.R +++ b/tests/testthat/test_loo_and_waic.R @@ -34,6 +34,18 @@ test_that("loo with cores=1 and cores=2 gives same results", { expect_equal(loo1$estimates, loo2$estimates) }) +test_that("mcse_elpd is stable for extreme log likelihoods", { + ll <- cbind(c(-1000, -1001), c(1000, 999)) + lw <- matrix(log(0.5), nrow = 2, ncol = 2) + E_elpd <- matrixStats::colLogSumExps(ll + lw) + shift <- apply(ll, 2, max) + lik <- exp(sweep(ll, 2, shift)) + E_epd <- exp(E_elpd - shift) + expected <- sqrt(log1p(colSums(exp(lw)^2 * (lik - E_epd)^2) / E_epd^2)) + + expect_equal(mcse_elpd(ll, lw, E_elpd, r_eff = 1), expected) +}) + test_that("waic returns object with correct structure", { expect_true(is.waic(waic1)) expect_true(is.loo(waic1)) From 4aa7c4d452a023f9ec6eb2a7d3fddf37b0a7be4f Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 12 Sep 2026 22:41:16 +0300 Subject: [PATCH 02/17] fix: stabilize weighted variance calculations --- R/E_loo.R | 3 ++- R/loo_moment_matching.R | 9 +++++---- tests/testthat/test_E_loo.R | 3 +++ tests/testthat/test_loo_moment_matching.R | 8 ++++++++ 4 files changed, 18 insertions(+), 5 deletions(-) diff --git a/R/E_loo.R b/R/E_loo.R index b449c182..b0a43f46 100644 --- a/R/E_loo.R +++ b/R/E_loo.R @@ -219,7 +219,8 @@ E_loo.matrix <- # sample size ESS is estimated with the generic target quantity invariant # estimate 1/sum(w^2), see e.g. "Monte Carlo theory, methods and examples" # by Owen (2013). - (sum(.wmean(x^2, w)) - sum(.wmean(x, w)^2)) / (1 - sum(w^2)) + weighted_mean <- .wmean(x, w) + sum(w * (x - weighted_mean)^2) / (1 - sum(w^2)) } .wsd <- function(x, w, ...) { sqrt(.wvar(x, w)) diff --git a/R/loo_moment_matching.R b/R/loo_moment_matching.R index 110eff93..7087ef16 100644 --- a/R/loo_moment_matching.R +++ b/R/loo_moment_matching.R @@ -531,11 +531,12 @@ shift_and_scale <- function(x, upars, lwi) { # compute moments using log weights S <- dim(upars)[1] mean_original <- colMeans(upars) - mean_weighted <- colSums(exp(lwi) * upars) + weights <- exp(lwi) + mean_weighted <- colSums(weights * upars) shift <- mean_weighted - mean_original - mii <- exp(lwi)* upars^2 - mii <- colSums(mii) - mean_weighted^2 - mii <- mii*S/(S-1) + centered <- sweep(upars, 2, mean_weighted) + mii <- colSums(weights * centered^2) + mii <- mii * S / (S - 1) scaling <- sqrt(mii / matrixStats::colVars(upars)) # transform posterior draws upars_new <- sweep(upars, 2, mean_original, "-") diff --git a/tests/testthat/test_E_loo.R b/tests/testthat/test_E_loo.R index a3bf5c81..4837c860 100644 --- a/tests/testthat/test_E_loo.R +++ b/tests/testthat/test_E_loo.R @@ -232,4 +232,7 @@ test_that("weighted variance works", { w <- c(rep(0.1, 10), rep(0, 90)) expect_equal(.wvar(x, w), var(x[w > 0])) + + x <- 1e8 + c(-1, 0, 1) + expect_equal(.wvar(x, rep(1 / 3, 3)), 1) }) diff --git a/tests/testthat/test_loo_moment_matching.R b/tests/testthat/test_loo_moment_matching.R index 2c756533..b444a019 100644 --- a/tests/testthat/test_loo_moment_matching.R +++ b/tests/testthat/test_loo_moment_matching.R @@ -515,6 +515,14 @@ test_that("loo_moment_match_split works", { expect_snapshot_value(split2, style = "serialize") }) +test_that("shift_and_scale is stable for large parameter values", { + upars <- matrix(1e8 + c(-1, 0, 1), ncol = 1) + result <- shift_and_scale(NULL, upars, rep(-log(3), 3)) + + expect_equal(result$scaling, 1) + expect_equal(result$upars, upars) +}) + test_that("passing arguments works", { log_lik_i_upars_test_additional_argument <- function( x, From 9103840b0998661e4d25c3e27c640297630d45ec Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 12 Sep 2026 22:42:26 +0300 Subject: [PATCH 03/17] fix: stabilize subsampling square differences --- R/loo_subsample.R | 7 ++++++- tests/testthat/test_loo_subsampling_approximations.R | 7 +++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/R/loo_subsample.R b/R/loo_subsample.R index bcac4b17..4c1adbdb 100644 --- a/R/loo_subsample.R +++ b/R/loo_subsample.R @@ -1166,6 +1166,11 @@ loo_subsample_estimation_diff_srs <- function(x) { update_psis_loo_ss_estimates(x) } +difference_of_squares <- function(x, y) { + (x - y) * (x + y) +} + + #' Difference estimation using SRS-WOR sampling (Magnusson et al., 2020) #' @noRd #' @param y_approx Approximated values of all observations. @@ -1185,7 +1190,7 @@ srs_diff_est <- function(y_approx, y, y_idx) { t_pi_tilde <- sum(y_approx) t_pi2_tilde <- sum(y_approx^2) t_e <- N * mean(e_i) - t_hat_epsilon <- N * mean(y^2 - y_approx_m^2) + t_hat_epsilon <- N * mean(difference_of_squares(y, y_approx_m)) est_list <- list(m = length(y), N = N) # eq (7) diff --git a/tests/testthat/test_loo_subsampling_approximations.R b/tests/testthat/test_loo_subsampling_approximations.R index e5a02d58..9e20d412 100644 --- a/tests/testthat/test_loo_subsampling_approximations.R +++ b/tests/testthat/test_loo_subsampling_approximations.R @@ -628,6 +628,13 @@ test_that("whhest works as expected", { }) +test_that("difference_of_squares is stable for nearby large values", { + x <- 1e8 + 1.4901161193847656e-8 + + expect_equal(difference_of_squares(x, 1e8), 2.9802322387695312) +}) + + test_that("srs_diff_est works as expected", { set.seed(1234) N <- 1000 From aa65ac6415a1a074feab78b071db3fb411abbf49 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 12 Sep 2026 22:43:23 +0300 Subject: [PATCH 04/17] fix: stabilize PSIS tail differences --- R/psis.R | 7 ++++++- tests/testthat/test_psis.R | 9 +++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/R/psis.R b/R/psis.R index dd8b0176..30a8ffa2 100644 --- a/R/psis.R +++ b/R/psis.R @@ -240,6 +240,11 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { list(log_weights = lw_i, pareto_k = khat) } +exp_x_minus_exp_y <- function(x, y) { + -exp(x) * expm1(y - x) +} + + #' PSIS tail smoothing for a single vector #' #' @noRd @@ -254,7 +259,7 @@ psis_smooth_tail <- function(x, cutoff) { exp_cutoff <- exp(cutoff) # save time not sorting since x already sorted - fit <- posterior::gpdfit(exp(x) - exp_cutoff, sort_x = FALSE) + fit <- posterior::gpdfit(exp_x_minus_exp_y(x, cutoff), sort_x = FALSE) k <- fit$k sigma <- fit$sigma if (is.na(k)) { diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R index 93d7501a..0d814aeb 100644 --- a/tests/testthat/test_psis.R +++ b/tests/testthat/test_psis.R @@ -152,6 +152,15 @@ test_that("do_psis_i throws warning if all tail values the same", { expect_equal(val$pareto_k, Inf) }) +test_that("exp_x_minus_exp_y is stable for nearby values", { + cutoff <- -30 + x <- cutoff + 1e-14 + + expect_equal(exp_x_minus_exp_y(x, cutoff), 9.973486536736925e-28) + expect_equal(exp_x_minus_exp_y(0, -1000), 1) +}) + + test_that("psis_smooth_tail returns original tail values if k is infinite", { xx <- log(c(2, 2, 2, 2, 3, 4, 5, 6)) val <- suppressWarnings(psis_smooth_tail(xx, 0)) From 5b2bee17f205d76928e420143c71ff427a969d30 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sat, 12 Sep 2026 22:44:49 +0300 Subject: [PATCH 05/17] fix: stabilize stacking gradient calculation --- R/loo_model_weights.R | 25 +++++++++++++++++++++---- tests/testthat/test_model_weighting.R | 8 ++++++++ 2 files changed, 29 insertions(+), 4 deletions(-) diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index 946dc7c3..c6ef2460 100644 --- a/R/loo_model_weights.R +++ b/R/loo_model_weights.R @@ -234,6 +234,19 @@ loo_model_weights.default <- } +exp_diff_over_exp <- function(a, b, denominator_log) { + a_is_larger <- a >= b + out <- numeric(length(a)) + out[a_is_larger] <- + exp(a[a_is_larger] - denominator_log[a_is_larger]) * + -expm1(b[a_is_larger] - a[a_is_larger]) + out[!a_is_larger] <- + exp(b[!a_is_larger] - denominator_log[!a_is_larger]) * + expm1(a[!a_is_larger] - b[!a_is_larger]) + out +} + + #' @rdname loo_model_weights #' @export #' @param lpd_point If calling `stacking_weights()` or `pseudobma_weights()` @@ -272,11 +285,15 @@ stacking_weights <- stopifnot(length(w) == K - 1) w_full <- c(w, 1 - sum(w)) grad <- rep(0, K - 1) - # avoid over- and underflows using log weights, rowLogSumExps, - # and by subtracting the row maximum of lpd_point - mlpd <- matrixStats::rowMaxs(lpd_point) + mixture_lpd <- matrixStats::rowLogSumExps( + sweep(lpd_point, 2, log(w_full), "+") + ) for (k in 1:(K - 1)) { - grad[k] <- sum((exp(lpd_point[, k] - mlpd) - exp(lpd_point[, K] - mlpd)) / exp(matrixStats::rowLogSumExps(sweep(lpd_point, 2, log(w_full), '+')) - mlpd)) + grad[k] <- sum(exp_diff_over_exp( + lpd_point[, k], + lpd_point[, K], + mixture_lpd + )) } return(-grad) } diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index 96c075b8..1c525c81 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -20,6 +20,14 @@ loo_list <- lapply(1:length(ll_list), function(j) { tol <- 0.01 # absolute tolerance of weights +test_that("stacking gradient is stable for similar model predictions", { + a <- -30 + 1e-14 + b <- -30 + + expect_equal(exp_diff_over_exp(a, b, b), 1.065814103640156e-14) + expect_equal(exp_diff_over_exp(b, a, b), -1.065814103640156e-14) +}) + test_that("loo_model_weights throws correct errors and warnings", { expect_error( loo_model_weights(log_lik1), From 572ae71172908d2af09fc2403cdad54cad17f220 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Sun, 13 Sep 2026 10:41:54 +0300 Subject: [PATCH 06/17] refactor: move new helper functions to helpers.R --- R/helpers.R | 40 ++++++++++++++++++++++++++++++++++++++++ R/loo_model_weights.R | 13 ------------- R/loo_subsample.R | 5 ----- R/psis.R | 5 ----- 4 files changed, 40 insertions(+), 23 deletions(-) diff --git a/R/helpers.R b/R/helpers.R index 38b401dd..c3086905 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -26,6 +26,46 @@ colLogMeanExps <- function(x) { matrixStats::colLogSumExps(x) - logS } +#' More stable version of `exp(x) - exp(y)` +#' +#' @noRd +#' @param x A numeric vector. +#' @param y A numeric scalar or vector recycled to the length of `x`. +#' Must satisfy `x >= y` elementwise. +#' @return A numeric vector equal to `exp(x) - exp(y)`. +#' +exp_x_minus_exp_y <- function(x, y) { + -exp(x) * expm1(y - x) +} + +#' More stable version of `x^2 - y^2` +#' +#' @noRd +#' @param x,y Numeric vectors of the same length. +#' @return A numeric vector equal to `x^2 - y^2`. +#' +difference_of_squares <- function(x, y) { + (x - y) * (x + y) +} + +#' More stable version of `(exp(a) - exp(b)) / exp(c)` +#' +#' @noRd +#' @param a,b,c Numeric vectors of the same length. +#' @return A numeric vector equal to `(exp(a) - exp(b)) / exp(c)`. +#' +exp_diff_over_exp <- function(a, b, c) { + a_is_larger <- a >= b + out <- numeric(length(a)) + out[a_is_larger] <- + exp(a[a_is_larger] - c[a_is_larger]) * + -expm1(b[a_is_larger] - a[a_is_larger]) + out[!a_is_larger] <- + exp(b[!a_is_larger] - c[!a_is_larger]) * + expm1(a[!a_is_larger] - b[!a_is_larger]) + out +} + #' Compute point estimates and standard errors from pointwise vectors #' #' @noRd diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index c6ef2460..d22ecb6c 100644 --- a/R/loo_model_weights.R +++ b/R/loo_model_weights.R @@ -234,19 +234,6 @@ loo_model_weights.default <- } -exp_diff_over_exp <- function(a, b, denominator_log) { - a_is_larger <- a >= b - out <- numeric(length(a)) - out[a_is_larger] <- - exp(a[a_is_larger] - denominator_log[a_is_larger]) * - -expm1(b[a_is_larger] - a[a_is_larger]) - out[!a_is_larger] <- - exp(b[!a_is_larger] - denominator_log[!a_is_larger]) * - expm1(a[!a_is_larger] - b[!a_is_larger]) - out -} - - #' @rdname loo_model_weights #' @export #' @param lpd_point If calling `stacking_weights()` or `pseudobma_weights()` diff --git a/R/loo_subsample.R b/R/loo_subsample.R index 4c1adbdb..4aad6d46 100644 --- a/R/loo_subsample.R +++ b/R/loo_subsample.R @@ -1166,11 +1166,6 @@ loo_subsample_estimation_diff_srs <- function(x) { update_psis_loo_ss_estimates(x) } -difference_of_squares <- function(x, y) { - (x - y) * (x + y) -} - - #' Difference estimation using SRS-WOR sampling (Magnusson et al., 2020) #' @noRd #' @param y_approx Approximated values of all observations. diff --git a/R/psis.R b/R/psis.R index 30a8ffa2..303da755 100644 --- a/R/psis.R +++ b/R/psis.R @@ -240,11 +240,6 @@ do_psis_i <- function(log_ratios_i, tail_len_i, ...) { list(log_weights = lw_i, pareto_k = khat) } -exp_x_minus_exp_y <- function(x, y) { - -exp(x) * expm1(y - x) -} - - #' PSIS tail smoothing for a single vector #' #' @noRd From a336b76f6681ac30d8186162a15726ab49109faf Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 20:38:39 +0300 Subject: [PATCH 07/17] fix: handle equal infinite PSIS values --- R/helpers.R | 5 ++++- tests/testthat/test_psis.R | 11 +++++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/R/helpers.R b/R/helpers.R index c3086905..fa7d2845 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -35,7 +35,10 @@ colLogMeanExps <- function(x) { #' @return A numeric vector equal to `exp(x) - exp(y)`. #' exp_x_minus_exp_y <- function(x, y) { - -exp(x) * expm1(y - x) + out <- -exp(x) * expm1(y - x) + equal <- x == y + out[!is.na(equal) & equal] <- 0 + out } #' More stable version of `x^2 - y^2` diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R index 0d814aeb..183f868a 100644 --- a/tests/testthat/test_psis.R +++ b/tests/testthat/test_psis.R @@ -152,12 +152,23 @@ test_that("do_psis_i throws warning if all tail values the same", { expect_equal(val$pareto_k, Inf) }) +test_that("psis handles negative infinite log ratios", { + log_ratios <- c(rep(-Inf, 90), seq(-9, 0, length.out = 10)) + + expect_no_error(out <- suppressWarnings(psis(log_ratios))) + expect_s3_class(out, "psis") +}) + test_that("exp_x_minus_exp_y is stable for nearby values", { cutoff <- -30 x <- cutoff + 1e-14 expect_equal(exp_x_minus_exp_y(x, cutoff), 9.973486536736925e-28) expect_equal(exp_x_minus_exp_y(0, -1000), 1) + expect_equal( + exp_x_minus_exp_y(c(-Inf, 0, Inf), c(-Inf, 0, Inf)), + c(0, 0, 0) + ) }) From 9f4341fec13e80b99a084c204aa0b7f45c4fca3e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 20:39:10 +0300 Subject: [PATCH 08/17] fix: handle equal infinite stacking terms --- R/helpers.R | 2 ++ tests/testthat/test_model_weighting.R | 14 ++++++++++++++ 2 files changed, 16 insertions(+) diff --git a/R/helpers.R b/R/helpers.R index fa7d2845..a11553b6 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -66,6 +66,8 @@ exp_diff_over_exp <- function(a, b, c) { out[!a_is_larger] <- exp(b[!a_is_larger] - c[!a_is_larger]) * expm1(a[!a_is_larger] - b[!a_is_larger]) + equal <- a == b + out[!is.na(equal) & equal] <- 0 out } diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index 1c525c81..28b5cf0f 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -26,6 +26,20 @@ test_that("stacking gradient is stable for similar model predictions", { expect_equal(exp_diff_over_exp(a, b, b), 1.065814103640156e-14) expect_equal(exp_diff_over_exp(b, a, b), -1.065814103640156e-14) + expect_equal( + exp_diff_over_exp(c(-Inf, 0, Inf), c(-Inf, 0, Inf), c(0, 0, 0)), + c(0, 0, 0) + ) +}) + +test_that("stacking handles equal infinite log predictive densities", { + lpd <- matrix( + c(-Inf, 0, -Inf, -1, -1, -1), + nrow = 2, + byrow = TRUE + ) + + expect_equal(as.numeric(stacking_weights(lpd)), c(0, 1, 0), tolerance = 1e-6) }) test_that("loo_model_weights throws correct errors and warnings", { From 23005f4b41ec4aa6016b197183ea8dc9aa77b551 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 22:16:00 +0300 Subject: [PATCH 09/17] fix: pass log-scale ELPD to moment matching MCSE --- R/loo_moment_matching.R | 2 +- tests/testthat/_snaps/loo_moment_matching.md | 266 +++++++++---------- tests/testthat/test_loo_moment_matching.R | 5 + 3 files changed, 139 insertions(+), 134 deletions(-) diff --git a/R/loo_moment_matching.R b/R/loo_moment_matching.R index 7087ef16..c31a12bb 100644 --- a/R/loo_moment_matching.R +++ b/R/loo_moment_matching.R @@ -408,7 +408,7 @@ loo_moment_match_i <- function(i, elpd_loo_i <- matrixStats::logSumExp(log_liki + lwi) mcse_elpd_loo <- mcse_elpd( ll = as.matrix(log_liki), lw = as.matrix(lwi), - E_elpd = exp(elpd_loo_i), r_eff = r_eff_i + E_elpd = elpd_loo_i, r_eff = r_eff_i ) list(elpd_loo_i = elpd_loo_i, diff --git a/tests/testthat/_snaps/loo_moment_matching.md b/tests/testthat/_snaps/loo_moment_matching.md index a227a3aa..1a1d41dc 100644 --- a/tests/testthat/_snaps/loo_moment_matching.md +++ b/tests/testthat/_snaps/loo_moment_matching.md @@ -32,44 +32,44 @@ # loo_moment_match.default works - WAoAAAACAAQFAAACAwAAAAMTAAAACgAAAg4AAAAGwFJQKLlLGUdAJeeelRE5/EBiUCi5SxlH - QDIH7c01qv5AJJ/C+xNRK0BCB+3NNar+AAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA + WAoAAAACAAQGAQACAwAAAAMTAAAACgAAAg4AAAAGwFJQKLlLGUlAJeeelRE6DEBiUCi5SxlJ + QDIH7c01qwlAJJ/C+xNRP0BCB+3NNasJAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA AwAAAAIAAAQCAAAAAQAEAAkAAAAIZGltbmFtZXMAAAATAAAAAgAAABAAAAADAAQACQAAAAhl bHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAAAEAAAAAIABAAJAAAACEVzdGlt - YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsAz2SCKL/OPv/uoQkpWQUi//t6gM8RRaL/7ME6e + YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsAz2SCKL/OZv/uoQkpWQUi//t6gM8RRaL/7ME6e 7mn4v/slBOxuIMC//9jJ8p3XcL/7L6Fev0U4wAA2/FFJKpi//ScFgWbouL/8Pocg7R74v/0p - qq1wTGC/+x8Kh2qjyL/7JFoKTeIwv/soK3EUBUC//KBEQBPCGMAAKtiBefwUv/s4omTXOVjA - A0qKEEuTHL/7hlp0sBHYv/xVS8d60Ki//yJy272SKL/7oVoD6+JIv/7iRaQt0cC//VaEo30D - QL/85OxxwadowAHqd1LUWNy/+9Uw/+PdqL/7IYmDCMJQv/+TLUcS2KC//UsYUnCS0D/ADGFe - xkQGP2Gfcpi2J8Y/Y+tB1l78/T9hmK3v3QgHP2GaKR7p1QA/ZQkD2NQUNz9hp3Fk+oG7P2Wd - 0KgLpYI/Yjujzb+UVT9hx0KNqoPZP2J1vTdzw+s/YaMJLQv8lD9hpM5D/++/P2GZy1TTq7Y/ - YfB+8HlL/j9lqO+CzesuP2GpxXMbu9A/cGE1wkkFgT9hvBsU728qP2HQSdm5/Yo/Y/zERK8U - MT9hni0KCSNZP2O2kFJN6R4/Ylo+R+/P/j9iFB/WgCCUP2s9ryP5flA/YdKJGMivnD9hmtGw - ba2EP2SA235CBJ8/YotQDpXP+EAkqqbkVm+LP5IXf6UkPQA/lrxpW1EOAD+SDmssqA4AP5IS + qq1wTGC/+x8Kh2qj0L/7JFoKTeIwv/soK3EUBUC//KBEQBPCGMAAKtiBefwUv/s4omTXOVjA + A0qKEEuTGL/7hlp0sBHYv/xVS8d60LC//yJy272SKL/7oVoD6+JIv/7iRaQt0cC//VaEo30D + QL/85OxxwadowAHqd1LUWNi/+9Uw/+PdoL/7IYmDCMJQv/+TLUcS2KC//UsYUnCS0D/IX47y + Iyd4P2Gfcpi2o14/Y+tB1l8doj9hmK3v3TkPP2GaKR7pqDc/ZQkD2NQ5Jj9hp3Fk+g3LP2Wd + 0KgLB8I/Yjujzb/pcj9hx0KNqzrIP2J1vTd0O6A/YaMJLQunTj9hpM5D/xQmP2GZy1TTNTQ/ + YfB+8HkBPD9lqO+CzYUVP2GpxXMcLEE/cGE1wkkqlD9hvBsU75ARP2HQSdm5R3s/Y/zERK8H + xD9hni0KCRVhP2O2kFJNSG8/Ylo+R/BQoj9iFB/Wf65IP2s9ryP57Cc/YdKJGMftEz9hmtGw + bW4sP2SA235BWtA/YotQDpT7XkAkqqbkVm+fP5IXf6UkPQA/lrxpW1EOAD+SDmssqA4AP5IS 8YLzxAA/mQEBdV99AD+SOdWCpaAAP5otz1Wr6wA/k0U02BITAD+SYx8kcE4AP5PhNpO18QA/ - ki0P3eyAAD+SMjlDn5oAP5IR1Ro8TgA/krLjpKmXAD+aTlBMoagAP5JAN/f7OgA/qwiqhKr0 - AD+SbOl4OJEAP5J03K3aTgA/lsJL+aFkAD+SFPC7p4gAP5Y04S7O5wA/k4CJvuH6AD+S+DEh - 3hkAP6Omxr93f4A/kp2i9JwgAD+SFOjEYZoAP5fSMVoWZAA/lAsiPEe5AEBD2SCKL/OPQAuo + ki0P3eyCAD+SMjlDn5oAP5IR1Ro8TgA/krLjpKmXAD+aTlBMoagAP5JAN/f7OgA/qwiqhKrz + AD+SbOl4OJEAP5J03K3aUAA/lsJL+aFkAD+SFPC7p4gAP5Y04S7O5wA/k4CJvuH6AD+S+DEh + 3hkAP6Omxr93foA/kp2i9JweAD+SFOjEYZoAP5fSMVoWZAA/lAsiPEe5AEBD2SCKL/OZQAuo QkpWQUhADt6gM8RRaEALME6e7mn4QAslBOxuIMBAD9jJ8p3XcEALL6Fev0U4QBA2/FFJKphA - DScFgWbouEAMPocg7R74QA0pqq1wTGBACx8Kh2qjyEALJFoKTeIwQAsoK3EUBUBADKBEQBPC - GEAQKtiBefwUQAs4omTXOVhAE0qKEEuTHEALhlp0sBHYQAxVS8d60KhADyJy272SKEALoVoD - 6+JIQA7iRaQt0cBADVaEo30DQEAM5OxxwadoQBHqd1LUWNxAC9Uw/+PdqEALIYmDCMJQQA+T - LUcS2KBADUsYUnCS0D/2XcffnmqVv5oq/ZSIVXI/rc6o9F9nBr9RT6Q+S52Qv3AJF5hupNw/ - teKv62UxfL+O7vQo7zgIP8HB+bh4VvU/p+nRG8MPFz+j+bKjZAlIv6dNik1Iv/W/pGhy2+ci - f7+iYyxrlH4Ov4shzHMAyAw/r8r08fJQBz+9B+WoAPW6v6EcmkpM9v8/sa5n3E19AL+gl0c6 - bbNbP5Z/I4ovlu4/xB68LmfXyr+KgBHuAEaAP7sysVWUlHQ/tAdcziRrAj+tf7e/7a72P6bM - edRbec6/pOE4IxDMN79Twd0WCzSwP8WPgi1BERa/kp3azvQypgAABAIAAAH/AAAADQAAAAIA + DScFgWbouEAMPocg7R74QA0pqq1wTGBACx8Kh2qj0EALJFoKTeIwQAsoK3EUBUBADKBEQBPC + GEAQKtiBefwUQAs4omTXOVhAE0qKEEuTGEALhlp0sBHYQAxVS8d60LBADyJy272SKEALoVoD + 6+JIQA7iRaQt0cBADVaEo30DQEAM5OxxwadoQBHqd1LUWNhAC9Uw/+PdoEALIYmDCMJQQA+T + LUcS2KBADUsYUnCS0D/2XcffnmqMv5oq/ZSIVTg/rc6o9F9myr9RT6Q+S9cKv3AJF5hurk0/ + teKv62UxTb+O7vQo7zfSP8HB+bh4V7g/p+nRG8MPsz+j+bKjZAlUv6dNik1Iv/q/pGhy2+ck + Sr+iYyxrlH76v4shzHMAxXs/r8r08fJP/D+9B+WoAPYGv6EcmkpM9eg/sa5n3E19L7+gl0c6 + bbJXP5Z/I4ovk/g/xB68LmfW37+KgBHuAEWmP7sysVWUlGw/tAdcziRrIj+tf7e/7a48P6bM + edRbekW/pOE4IxDLd79Twd0WCzszP8WPgi1BEUO/kp3azvQ1TQAABAIAAAH/AAAADQAAAAIA AAAeAAAABQAABAIAAAL/AAAAEwAAAAIAAAD+AAAAEAAAAAUABAAJAAAACGVscGRfbG9vAAQA CQAAAA1tY3NlX2VscGRfbG9vAAQACQAAAAVwX2xvbwAEAAkAAAAFbG9vaWMABAAJAAAAEmlu - Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAeP+fdHDQ/QUa/mir9lIhVcj+t - zqj0X2cGv1FPpD5LnZC/cAkXmG6k3D+14q/rZTF8v47u9CjvOAg/wcH5uHhW9T+n6dEbww8X - P6P5sqNkCUi/p02KTUi/9b+kaHLb5yJ/v6JjLGuUfg6/iyHMcwDIDD+vyvTx8lAHP70H5agA - 9bq/oRyaSkz2/z+xrmfcTX0Av6CXRzpts1s/ln8jii+W7j/EHrwuZ9fKv4qAEe4ARoA/uzKx - VZSUdD+0B1zOJGsCP61/t7/trvY/psx51Ft5zr+k4TgjEMw3v1PB3RYLNLA/xY+CLUERFr+S - ndrO9DKmAAAADgAAAB5ARQcvr5nxRUCurmCYTUMgQK6GZIiNpldArq7KVLrefECurq6LsT47 + Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAeP+fdHDQ/Q3e/mir9lIhVOD+t + zqj0X2bKv1FPpD5L1wq/cAkXmG6uTT+14q/rZTFNv47u9CjvN9I/wcH5uHhXuD+n6dEbww+z + P6P5sqNkCVS/p02KTUi/+r+kaHLb5yRKv6JjLGuUfvq/iyHMcwDFez+vyvTx8k/8P70H5agA + 9ga/oRyaSkz16D+xrmfcTX0vv6CXRzptslc/ln8jii+T+D/EHrwuZ9bfv4qAEe4ARaY/uzKx + VZSUbD+0B1zOJGsiP61/t7/trjw/psx51Ft6Rb+k4TgjEMt3v1PB3RYLOzM/xY+CLUERQ7+S + ndrO9DVNAAAADgAAAB5ARQcvr5nxk0CurmCYTUMgQK6GZIiNpldArq7KVLrefECurq6LsT47 QK5xjnyab6BArq3SHuGPKUCuZrTzmGSJQK6kDS9Qm5FArqvKHptSZ0CuoDIHBkg3QK6uB1vl - FFtArq3uNy14MkCurq95Y/PBQK6pD98CEa1ArmWF3fh0P0CuraLKskjiQK1cwbtJXeNArqxg - HKDKmkCuqx85Cm9rQK6FgdGWcINArq598JLLRECuil5ipQqpQK6iHFktv3tArqazEBI7tkCt + FFtArq3uNy14MkCurq95Y/PDQK6pD98CEa1ArmWF3fh0P0CuraLKskjiQK1cwbtJXeNArqxg + HKDKnECuqx85Cm9rQK6FgdGWcINArq598JLLRECuil5ipQqpQK6iHFktv3lArqazEBI7tkCt 65ndUr6PQK6q4cDxAOlArq6k0f++pUCufAQQaRkWQK6e1Q/LzlsAAAAOAAAAHj/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAA AAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAA @@ -77,8 +77,8 @@ AAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAAAABAIAAAABAAQACQAAAAVuYW1lcwAAABAAAAADAAQACQAAAAhw YXJldG9fawAEAAkAAAAFbl9lZmYABAAJAAAABXJfZWZmAAAA/gAAAP4AAAAOAAAAAcBSUCi5 - SxlHAAAADgAAAAFAJeeelRE5/AAAAA4AAAABQGJQKLlLGUcAAAAOAAAAAUAyB+3NNar+AAAA - DgAAAAFAJJ/C+xNRKwAAAA4AAAABQEIH7c01qv4AAAQCAAAD/wAAABAAAAAKAAQACQAAAAll + SxlJAAAADgAAAAFAJeeelRE6DAAAAA4AAAABQGJQKLlLGUkAAAAOAAAAAUAyB+3NNasJAAAA + DgAAAAFAJJ/C+xNRPwAAAA4AAAABQEIH7c01qwkAAAQCAAAD/wAAABAAAAAKAAQACQAAAAll c3RpbWF0ZXMABAAJAAAACXBvaW50d2lzZQAEAAkAAAALZGlhZ25vc3RpY3MABAAJAAAAC3Bz aXNfb2JqZWN0AAQACQAAAAhlbHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAQA CQAAAAtzZV9lbHBkX2xvbwAEAAkAAAAIc2VfcF9sb28ABAAJAAAACHNlX2xvb2ljAAAEAgAA @@ -88,44 +88,44 @@ --- - WAoAAAACAAQFAAACAwAAAAMTAAAACgAAAg4AAAAGwFJ5ESZa2ndAJy7h/Y9DekBieREmWtp3 - QDKrWN/9QglAJecFy/N4C0BCq1jf/UIJAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA + WAoAAAACAAQGAQACAwAAAAMTAAAACgAAAg4AAAAGwFJ5ESZa2nlAJy7h/Y9DjEBieREmWtp5 + QDKrWN/9QhVAJecFy/N4I0BCq1jf/UIVAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA AwAAAAIAAAQCAAAAAQAEAAkAAAAIZGltbmFtZXMAAAATAAAAAgAAABAAAAADAAQACQAAAAhl bHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAAAEAAAAAIABAAJAAAACEVzdGlt - YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsA0fMI+bvhOv/uoQkpWQUi//t6gM8RRaL/7ME6e + YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsA0fMI+bvhZv/uoQkpWQUi//t6gM8RRaL/7ME6e 7mn4v/slBOxuIMC//9jJ8p3XcL/7L6Fev0U4wAA2/FFJKpi//ScFgWbouL/8Pocg7R74v/0p - qq1wTGC/+x8Kh2qjyL/7JFoKTeIwv/soK3EUBUC//KBEQBPCGMAAKtiBefwUv/s4omTXOVjA - A0qKEEuTHL/7hlp0sBHYv/xVS8d60Ki//yJy272SKL/7oVoD6+JIv/7iRaQt0cC//VaEo30D - QL/85OxxwadowAHqd1LUWNy/+9Uw/+PdqL/7IYmDCMJQv/+TLUcS2KC//UsYUnCS0D+TxTjn - PtpbP2Gfcpi2J8Y/Y+tB1l78/T9hmK3v3QgHP2GaKR7p1QA/ZQkD2NQUNz9hp3Fk+oG7P2Wd - 0KgLpYI/Yjujzb+UVT9hx0KNqoPZP2J1vTdzw+s/YaMJLQv8lD9hpM5D/++/P2GZy1TTq7Y/ - YfB+8HlL/j9lqO+CzesuP2GpxXMbu9A/cGE1wkkFgT9hvBsU728qP2HQSdm5/Yo/Y/zERK8U - MT9hni0KCSNZP2O2kFJN6R4/Ylo+R+/P/j9iFB/WgCCUP2s9ryP5flA/YdKJGMivnD9hmtGw - ba2EP2SA235CBJ8/YotQDpXP+EAl8epM1HkJP5IXf6UkPQA/lrxpW1EOAD+SDmssqA4AP5IS + qq1wTGC/+x8Kh2qj0L/7JFoKTeIwv/soK3EUBUC//KBEQBPCGMAAKtiBefwUv/s4omTXOVjA + A0qKEEuTGL/7hlp0sBHYv/xVS8d60LC//yJy272SKL/7oVoD6+JIv/7iRaQt0cC//VaEo30D + QL/85OxxwadowAHqd1LUWNi/+9Uw/+PdoL/7IYmDCMJQv/+TLUcS2KC//UsYUnCS0D/Zw5bt + Aw0/P2Gfcpi2o14/Y+tB1l8doj9hmK3v3TkPP2GaKR7pqDc/ZQkD2NQ5Jj9hp3Fk+g3LP2Wd + 0KgLB8I/Yjujzb/pcj9hx0KNqzrIP2J1vTd0O6A/YaMJLQunTj9hpM5D/xQmP2GZy1TTNTQ/ + YfB+8HkBPD9lqO+CzYUVP2GpxXMcLEE/cGE1wkkqlD9hvBsU75ARP2HQSdm5R3s/Y/zERK8H + xD9hni0KCRVhP2O2kFJNSG8/Ylo+R/BQoj9iFB/Wf65IP2s9ryP57Cc/YdKJGMftEz9hmtGw + bW4sP2SA235BWtA/YotQDpT7XkAl8epM1HkfP5IXf6UkPQA/lrxpW1EOAD+SDmssqA4AP5IS 8YLzxAA/mQEBdV99AD+SOdWCpaAAP5otz1Wr6wA/k0U02BITAD+SYx8kcE4AP5PhNpO18QA/ - ki0P3eyAAD+SMjlDn5oAP5IR1Ro8TgA/krLjpKmXAD+aTlBMoagAP5JAN/f7OgA/qwiqhKr0 - AD+SbOl4OJEAP5J03K3aTgA/lsJL+aFkAD+SFPC7p4gAP5Y04S7O5wA/k4CJvuH6AD+S+DEh - 3hkAP6Omxr93f4A/kp2i9JwgAD+SFOjEYZoAP5fSMVoWZAA/lAsiPEe5AEBEfMI+bvhOQAuo + ki0P3eyCAD+SMjlDn5oAP5IR1Ro8TgA/krLjpKmXAD+aTlBMoagAP5JAN/f7OgA/qwiqhKrz + AD+SbOl4OJEAP5J03K3aUAA/lsJL+aFkAD+SFPC7p4gAP5Y04S7O5wA/k4CJvuH6AD+S+DEh + 3hkAP6Omxr93foA/kp2i9JweAD+SFOjEYZoAP5fSMVoWZAA/lAsiPEe5AEBEfMI+bvhZQAuo QkpWQUhADt6gM8RRaEALME6e7mn4QAslBOxuIMBAD9jJ8p3XcEALL6Fev0U4QBA2/FFJKphA - DScFgWbouEAMPocg7R74QA0pqq1wTGBACx8Kh2qjyEALJFoKTeIwQAsoK3EUBUBADKBEQBPC - GEAQKtiBefwUQAs4omTXOVhAE0qKEEuTHEALhlp0sBHYQAxVS8d60KhADyJy272SKEALoVoD - 6+JIQA7iRaQt0cBADVaEo30DQEAM5OxxwadoQBHqd1LUWNxAC9Uw/+PdqEALIYmDCMJQQA+T - LUcS2KBADUsYUnCS0D/2XcffnmqVv5oq/ZSIVXI/rc6o9F9nBr9RT6Q+S52Qv3AJF5hupNw/ - teKv62UxfL+O7vQo7zgIP8HB+bh4VvU/p+nRG8MPFz+j+bKjZAlIv6dNik1Iv/W/pGhy2+ci - f7+iYyxrlH4Ov4shzHMAyAw/r8r08fJQBz+9B+WoAPW6v6EcmkpM9v8/sa5n3E19AL+gl0c6 - bbNbP5Z/I4ovlu4/xB68LmfXyr+KgBHuAEaAP7sysVWUlHQ/tAdcziRrAj+tf7e/7a72P6bM - edRbec6/pOE4IxDMN79Twd0WCzSwP8WPgi1BERa/kp3azvQypgAABAIAAAH/AAAADQAAAAIA + DScFgWbouEAMPocg7R74QA0pqq1wTGBACx8Kh2qj0EALJFoKTeIwQAsoK3EUBUBADKBEQBPC + GEAQKtiBefwUQAs4omTXOVhAE0qKEEuTGEALhlp0sBHYQAxVS8d60LBADyJy272SKEALoVoD + 6+JIQA7iRaQt0cBADVaEo30DQEAM5OxxwadoQBHqd1LUWNhAC9Uw/+PdoEALIYmDCMJQQA+T + LUcS2KBADUsYUnCS0D/2XcffnmqMv5oq/ZSIVTg/rc6o9F9myr9RT6Q+S9cKv3AJF5hurk0/ + teKv62UxTb+O7vQo7zfSP8HB+bh4V7g/p+nRG8MPsz+j+bKjZAlUv6dNik1Iv/q/pGhy2+ck + Sr+iYyxrlH76v4shzHMAxXs/r8r08fJP/D+9B+WoAPYGv6EcmkpM9eg/sa5n3E19L7+gl0c6 + bbJXP5Z/I4ovk/g/xB68LmfW37+KgBHuAEWmP7sysVWUlGw/tAdcziRrIj+tf7e/7a48P6bM + edRbekW/pOE4IxDLd79Twd0WCzszP8WPgi1BEUO/kp3azvQ1TQAABAIAAAH/AAAADQAAAAIA AAAeAAAABQAABAIAAAL/AAAAEwAAAAIAAAD+AAAAEAAAAAUABAAJAAAACGVscGRfbG9vAAQA CQAAAA1tY3NlX2VscGRfbG9vAAQACQAAAAVwX2xvbwAEAAkAAAAFbG9vaWMABAAJAAAAEmlu - Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAeP7Z3M/Z+3oW/mir9lIhVcj+t - zqj0X2cGv1FPpD5LnZC/cAkXmG6k3D+14q/rZTF8v47u9CjvOAg/wcH5uHhW9T+n6dEbww8X - P6P5sqNkCUi/p02KTUi/9b+kaHLb5yJ/v6JjLGuUfg6/iyHMcwDIDD+vyvTx8lAHP70H5agA - 9bq/oRyaSkz2/z+xrmfcTX0Av6CXRzpts1s/ln8jii+W7j/EHrwuZ9fKv4qAEe4ARoA/uzKx - VZSUdD+0B1zOJGsCP61/t7/trvY/psx51Ft5zr+k4TgjEMw3v1PB3RYLNLA/xY+CLUERFr+S - ndrO9DKmAAAADgAAAB5AFp345qujjUCurmCYTUMgQK6GZIiNpldArq7KVLrefECurq6LsT47 + Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAeP7Z3M/Z+5by/mir9lIhVOD+t + zqj0X2bKv1FPpD5L1wq/cAkXmG6uTT+14q/rZTFNv47u9CjvN9I/wcH5uHhXuD+n6dEbww+z + P6P5sqNkCVS/p02KTUi/+r+kaHLb5yRKv6JjLGuUfvq/iyHMcwDFez+vyvTx8k/8P70H5agA + 9ga/oRyaSkz16D+xrmfcTX0vv6CXRzptslc/ln8jii+T+D/EHrwuZ9bfv4qAEe4ARaY/uzKx + VZSUbD+0B1zOJGsiP61/t7/trjw/psx51Ft6Rb+k4TgjEMt3v1PB3RYLOzM/xY+CLUERQ7+S + ndrO9DVNAAAADgAAAB5AFp345qujwkCurmCYTUMgQK6GZIiNpldArq7KVLrefECurq6LsT47 QK5xjnyab6BArq3SHuGPKUCuZrTzmGSJQK6kDS9Qm5FArqvKHptSZ0CuoDIHBkg3QK6uB1vl - FFtArq3uNy14MkCurq95Y/PBQK6pD98CEa1ArmWF3fh0P0CuraLKskjiQK1cwbtJXeNArqxg - HKDKmkCuqx85Cm9rQK6FgdGWcINArq598JLLRECuil5ipQqpQK6iHFktv3tArqazEBI7tkCt + FFtArq3uNy14MkCurq95Y/PDQK6pD98CEa1ArmWF3fh0P0CuraLKskjiQK1cwbtJXeNArqxg + HKDKnECuqx85Cm9rQK6FgdGWcINArq598JLLRECuil5ipQqpQK6iHFktv3lArqazEBI7tkCt 65ndUr6PQK6q4cDxAOlArq6k0f++pUCufAQQaRkWQK6e1Q/LzlsAAAAOAAAAHj/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAA AAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAA @@ -133,8 +133,8 @@ AAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAAAABAIAAAABAAQACQAAAAVuYW1lcwAAABAAAAADAAQACQAAAAhw YXJldG9fawAEAAkAAAAFbl9lZmYABAAJAAAABXJfZWZmAAAA/gAAAP4AAAAOAAAAAcBSeREm - Wtp3AAAADgAAAAFAJy7h/Y9DegAAAA4AAAABQGJ5ESZa2ncAAAAOAAAAAUAyq1jf/UIJAAAA - DgAAAAFAJecFy/N4CwAAAA4AAAABQEKrWN/9QgkAAAQCAAAD/wAAABAAAAAKAAQACQAAAAll + Wtp5AAAADgAAAAFAJy7h/Y9DjAAAAA4AAAABQGJ5ESZa2nkAAAAOAAAAAUAyq1jf/UIVAAAA + DgAAAAFAJecFy/N4IwAAAA4AAAABQEKrWN/9QhUAAAQCAAAD/wAAABAAAAAKAAQACQAAAAll c3RpbWF0ZXMABAAJAAAACXBvaW50d2lzZQAEAAkAAAALZGlhZ25vc3RpY3MABAAJAAAAC3Bz aXNfb2JqZWN0AAQACQAAAAhlbHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAQA CQAAAAtzZV9lbHBkX2xvbwAEAAkAAAAIc2VfcF9sb28ABAAJAAAACHNlX2xvb2ljAAAEAgAA @@ -144,44 +144,44 @@ --- - WAoAAAACAAQFAAACAwAAAAMTAAAACgAAAg4AAAAGwFKt/dMmDJVAKNZHY+jUZEBirf3TJgyV - QDN+yjVXeO1AJ45qhpBmBkBDfso1V3jtAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA + WAoAAAACAAQGAQACAwAAAAMTAAAACgAAAg4AAAAGwFKt/dMmDJVAKNZHY+jUaEBirf3TJgyV + QDN+yjVXePBAJ45qhpBmDEBDfso1V3jwAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA AwAAAAIAAAQCAAAAAQAEAAkAAAAIZGltbmFtZXMAAAATAAAAAgAAABAAAAADAAQACQAAAAhl bHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAAAEAAAAAIABAAJAAAACEVzdGlt - YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsA1UHTxm8DDv/uoQkpWQUi//t6gM8RRaL/7ME6e + YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsA1UHTxm8DHv/uoQkpWQUi//t6gM8RRaL/7ME6e 7mn4v/slBOxuIMC//9jJ8p3XcL/7L6Fev0U4wAA2/FFJKpi//ScFgWbouL/8Pocg7R74v/0p - qq1wTGC/+x8Kh2qjyL/7JFoKTeIwv/soK3EUBUC//KBEQBPCGMAAKtiBefwUv/s4omTXOVjA - A0qKEEuTHL/7hlp0sBHYv/xVS8d60Ki//yJy272SKL/7oVoD6+JIv/7iRaQt0cC//VaEo30D - QL/85OxxwadowAHqd1LUWNy/+9Uw/+PdqL/7IYmDCMJQv/+TLUcS2KC//UsYUnCS0D+cqILP - qzBKP2Gfcpi2J8Y/Y+tB1l78/T9hmK3v3QgHP2GaKR7p1QA/ZQkD2NQUNz9hp3Fk+oG7P2Wd - 0KgLpYI/Yjujzb+UVT9hx0KNqoPZP2J1vTdzw+s/YaMJLQv8lD9hpM5D/++/P2GZy1TTq7Y/ - YfB+8HlL/j9lqO+CzesuP2GpxXMbu9A/cGE1wkkFgT9hvBsU728qP2HQSdm5/Yo/Y/zERK8U - MT9hni0KCSNZP2O2kFJN6R4/Ylo+R+/P/j9iFB/WgCCUP2s9ryP5flA/YdKJGMivnD9hmtGw - ba2EP2SA235CBJ8/YotQDpXP+EAnmU+zLgnzP5IXf6UkPQA/lrxpW1EOAD+SDmssqA4AP5IS + qq1wTGC/+x8Kh2qj0L/7JFoKTeIwv/soK3EUBUC//KBEQBPCGMAAKtiBefwUv/s4omTXOVjA + A0qKEEuTGL/7hlp0sBHYv/xVS8d60LC//yJy272SKL/7oVoD6+JIv/7iRaQt0cC//VaEo30D + QL/85OxxwadowAHqd1LUWNi/+9Uw/+PdoL/7IYmDCMJQv/+TLUcS2KC//UsYUnCS0D+iQ4TC + EOOrP2Gfcpi2o14/Y+tB1l8doj9hmK3v3TkPP2GaKR7pqDc/ZQkD2NQ5Jj9hp3Fk+g3LP2Wd + 0KgLB8I/Yjujzb/pcj9hx0KNqzrIP2J1vTd0O6A/YaMJLQunTj9hpM5D/xQmP2GZy1TTNTQ/ + YfB+8HkBPD9lqO+CzYUVP2GpxXMcLEE/cGE1wkkqlD9hvBsU75ARP2HQSdm5R3s/Y/zERK8H + xD9hni0KCRVhP2O2kFJNSG8/Ylo+R/BQoj9iFB/Wf65IP2s9ryP57Cc/YdKJGMftEz9hmtGw + bW4sP2SA235BWtA/YotQDpT7XkAnmU+zLgn7P5IXf6UkPQA/lrxpW1EOAD+SDmssqA4AP5IS 8YLzxAA/mQEBdV99AD+SOdWCpaAAP5otz1Wr6wA/k0U02BITAD+SYx8kcE4AP5PhNpO18QA/ - ki0P3eyAAD+SMjlDn5oAP5IR1Ro8TgA/krLjpKmXAD+aTlBMoagAP5JAN/f7OgA/qwiqhKr0 - AD+SbOl4OJEAP5J03K3aTgA/lsJL+aFkAD+SFPC7p4gAP5Y04S7O5wA/k4CJvuH6AD+S+DEh - 3hkAP6Omxr93f4A/kp2i9JwgAD+SFOjEYZoAP5fSMVoWZAA/lAsiPEe5AEBFUHTxm8DDQAuo + ki0P3eyCAD+SMjlDn5oAP5IR1Ro8TgA/krLjpKmXAD+aTlBMoagAP5JAN/f7OgA/qwiqhKrz + AD+SbOl4OJEAP5J03K3aUAA/lsJL+aFkAD+SFPC7p4gAP5Y04S7O5wA/k4CJvuH6AD+S+DEh + 3hkAP6Omxr93foA/kp2i9JweAD+SFOjEYZoAP5fSMVoWZAA/lAsiPEe5AEBFUHTxm8DHQAuo QkpWQUhADt6gM8RRaEALME6e7mn4QAslBOxuIMBAD9jJ8p3XcEALL6Fev0U4QBA2/FFJKphA - DScFgWbouEAMPocg7R74QA0pqq1wTGBACx8Kh2qjyEALJFoKTeIwQAsoK3EUBUBADKBEQBPC - GEAQKtiBefwUQAs4omTXOVhAE0qKEEuTHEALhlp0sBHYQAxVS8d60KhADyJy272SKEALoVoD - 6+JIQA7iRaQt0cBADVaEo30DQEAM5OxxwadoQBHqd1LUWNxAC9Uw/+PdqEALIYmDCMJQQA+T - LUcS2KBADUsYUnCS0D/2XcffnmqVv5oq/ZSIVXI/rc6o9F9nBr9RT6Q+S52Qv3AJF5hupNw/ - teKv62UxfL+O7vQo7zgIP8HB+bh4VvU/p+nRG8MPFz+j+bKjZAlIv6dNik1Iv/W/pGhy2+ci - f7+iYyxrlH4Ov4shzHMAyAw/r8r08fJQBz+9B+WoAPW6v6EcmkpM9v8/sa5n3E19AL+gl0c6 - bbNbP5Z/I4ovlu4/xB68LmfXyr+KgBHuAEaAP7sysVWUlHQ/tAdcziRrAj+tf7e/7a72P6bM - edRbec6/pOE4IxDMN79Twd0WCzSwP8WPgi1BERa/kp3azvQypgAABAIAAAH/AAAADQAAAAIA + DScFgWbouEAMPocg7R74QA0pqq1wTGBACx8Kh2qj0EALJFoKTeIwQAsoK3EUBUBADKBEQBPC + GEAQKtiBefwUQAs4omTXOVhAE0qKEEuTGEALhlp0sBHYQAxVS8d60LBADyJy272SKEALoVoD + 6+JIQA7iRaQt0cBADVaEo30DQEAM5OxxwadoQBHqd1LUWNhAC9Uw/+PdoEALIYmDCMJQQA+T + LUcS2KBADUsYUnCS0D/2XcffnmqMv5oq/ZSIVTg/rc6o9F9myr9RT6Q+S9cKv3AJF5hurk0/ + teKv62UxTb+O7vQo7zfSP8HB+bh4V7g/p+nRG8MPsz+j+bKjZAlUv6dNik1Iv/q/pGhy2+ck + Sr+iYyxrlH76v4shzHMAxXs/r8r08fJP/D+9B+WoAPYGv6EcmkpM9eg/sa5n3E19L7+gl0c6 + bbJXP5Z/I4ovk/g/xB68LmfW37+KgBHuAEWmP7sysVWUlGw/tAdcziRrIj+tf7e/7a48P6bM + edRbekW/pOE4IxDLd79Twd0WCzszP8WPgi1BEUO/kp3azvQ1TQAABAIAAAH/AAAADQAAAAIA AAAeAAAABQAABAIAAAL/AAAAEwAAAAIAAAD+AAAAEAAAAAUABAAJAAAACGVscGRfbG9vAAQA CQAAAA1tY3NlX2VscGRfbG9vAAQACQAAAAVwX2xvbwAEAAkAAAAFbG9vaWMABAAJAAAAEmlu - Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAeP7Z3M/Z+3oW/mir9lIhVcj+t - zqj0X2cGv1FPpD5LnZC/cAkXmG6k3D+14q/rZTF8v47u9CjvOAg/wcH5uHhW9T+n6dEbww8X - P6P5sqNkCUi/p02KTUi/9b+kaHLb5yJ/v6JjLGuUfg6/iyHMcwDIDD+vyvTx8lAHP70H5agA - 9bq/oRyaSkz2/z+xrmfcTX0Av6CXRzpts1s/ln8jii+W7j/EHrwuZ9fKv4qAEe4ARoA/uzKx - VZSUdD+0B1zOJGsCP61/t7/trvY/psx51Ft5zr+k4TgjEMw3v1PB3RYLNLA/xY+CLUERFr+S - ndrO9DKmAAAADgAAAB5Ak/Dwg4seSECurmCYTUMgQK6GZIiNpldArq7KVLrefECurq6LsT47 + Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAeP7Z3M/Z+5by/mir9lIhVOD+t + zqj0X2bKv1FPpD5L1wq/cAkXmG6uTT+14q/rZTFNv47u9CjvN9I/wcH5uHhXuD+n6dEbww+z + P6P5sqNkCVS/p02KTUi/+r+kaHLb5yRKv6JjLGuUfvq/iyHMcwDFez+vyvTx8k/8P70H5agA + 9ga/oRyaSkz16D+xrmfcTX0vv6CXRzptslc/ln8jii+T+D/EHrwuZ9bfv4qAEe4ARaY/uzKx + VZSUbD+0B1zOJGsiP61/t7/trjw/psx51Ft6Rb+k4TgjEMt3v1PB3RYLOzM/xY+CLUERQ7+S + ndrO9DVNAAAADgAAAB5Ak/Dwg4sez0CurmCYTUMgQK6GZIiNpldArq7KVLrefECurq6LsT47 QK5xjnyab6BArq3SHuGPKUCuZrTzmGSJQK6kDS9Qm5FArqvKHptSZ0CuoDIHBkg3QK6uB1vl - FFtArq3uNy14MkCurq95Y/PBQK6pD98CEa1ArmWF3fh0P0CuraLKskjiQK1cwbtJXeNArqxg - HKDKmkCuqx85Cm9rQK6FgdGWcINArq598JLLRECuil5ipQqpQK6iHFktv3tArqazEBI7tkCt + FFtArq3uNy14MkCurq95Y/PDQK6pD98CEa1ArmWF3fh0P0CuraLKskjiQK1cwbtJXeNArqxg + HKDKnECuqx85Cm9rQK6FgdGWcINArq598JLLRECuil5ipQqpQK6iHFktv3lArqazEBI7tkCt 65ndUr6PQK6q4cDxAOlArq6k0f++pUCufAQQaRkWQK6e1Q/LzlsAAAAOAAAAHj/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAA AAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAA @@ -189,8 +189,8 @@ AAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAAAABAIAAAABAAQACQAAAAVuYW1lcwAAABAAAAADAAQACQAAAAhw YXJldG9fawAEAAkAAAAFbl9lZmYABAAJAAAABXJfZWZmAAAA/gAAAP4AAAAOAAAAAcBSrf3T - JgyVAAAADgAAAAFAKNZHY+jUZAAAAA4AAAABQGKt/dMmDJUAAAAOAAAAAUAzfso1V3jtAAAA - DgAAAAFAJ45qhpBmBgAAAA4AAAABQEN+yjVXeO0AAAQCAAAD/wAAABAAAAAKAAQACQAAAAll + JgyVAAAADgAAAAFAKNZHY+jUaAAAAA4AAAABQGKt/dMmDJUAAAAOAAAAAUAzfso1V3jwAAAA + DgAAAAFAJ45qhpBmDAAAAA4AAAABQEN+yjVXePAAAAQCAAAD/wAAABAAAAAKAAQACQAAAAll c3RpbWF0ZXMABAAJAAAACXBvaW50d2lzZQAEAAkAAAALZGlhZ25vc3RpY3MABAAJAAAAC3Bz aXNfb2JqZWN0AAQACQAAAAhlbHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAQA CQAAAAtzZV9lbHBkX2xvbwAEAAkAAAAIc2VfcF9sb28ABAAJAAAACHNlX2xvb2ljAAAEAgAA @@ -200,53 +200,53 @@ # variance and covariance transformations work - WAoAAAACAAQFAAACAwAAAAMTAAAACgAAAg4AAAAGwFKLcIr2U6lAJ7mJXVQ1mEBii3CK9lOp - QDMBFEV6xQ1AJnj9Bj3oV0BDARRFesUNAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA + WAoAAAACAAQGAQACAwAAAAMTAAAACgAAAg4AAAAGwFKLcIr2U6dAJ7mJXVQ1jkBii3CK9lOn + QDMBFEV6xQhAJnj9Bj3oUEBDARRFesUIAAAEAgAAAAEABAAJAAAAA2RpbQAAAA0AAAACAAAA AwAAAAIAAAQCAAAAAQAEAAkAAAAIZGltbmFtZXMAAAATAAAAAgAAABAAAAADAAQACQAAAAhl bHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAAAEAAAAAIABAAJAAAACEVzdGlt - YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsA00hUKopC/v/uiIMLHLuy//ssyhve2vL/7Jr8K - SbUAv/sa1cItPhC//8POqw85FL/7JoUvY15YwAA7D8sD/c6//SUayysMjL/8O2Quh638v/0W - mUgAlIC/+xKZXXIQsL/7F40nuPLAv/seLsVhXcC//J23rmyzmMAAIARDWDzAv/srA2aD7MjA - A1R95vJwGL/7djZ68aYYv/xSR/sGCki//yevfroaBL/7mxSAeStQv/7m3aUwM2i//VTP19jb - WL/84rONukOUwAHx9fU/2Fa/+8QvunAEmL/7Fx4TRHXUv/+ZdSn1rdS//Tf3NN2ngD+W4GjR - aryTP5P9k5w+sG4/ayNFdQtzdT+UNR65X0TEP5QzA295mvs/bFTXaJrIOz+U3CtVdiGyP5RR - Ot5s8NI/aT33KWiSfT9pDE2N/1sTP5S2SpsPGnU/lHsSPZvdgj+UhiM2HUpuP5QsiYrA7e4/ - aRVFdDks8D9tDHOUss4tP5Sd+c05JPQ/lYiPg3JzFD+UxOlt8wDqP2kMYnIOue0/lBcHGRWx - ez+T/tJ6LYPNP5QNq0mErcY/aVODDQzbTz9pJSlA+yZNP5UHY/IZmJA/lNKlqN7G3T+UO0Qc - JW5GP5QpKO2aIqQ/lK/I/7EmaEAmg6Qr+gHKP5I5hg6AAwA/lTqBKHTFAD+SIcNl86MAP5If - omEDgQA/lvi9WmUhAD+THZI/2f4AP5mPg19ObgA/ktItKsZkAD+SjsiR3esAP5KkmoaHBgA/ - khZKaRj3AD+SFGHT2/EAP5Ig3ZVNGAA/kpwrkaAHAD+YCFgsLCIAP5IPS9yzsgA/qlFcWiEQ - gD+RznOE6nMAP5KPYOKjMwA/lmRIfQIhAD+SN2h/dggAP5Xg9zNtrgA/ku5c/6I+AD+SsXhR - Iu0AP6MUf55iiAA/kd0sR7uFAD+SH4L0gRUAP5dfNWv4+gA/kshKi4uJAEBE0hUKopC/QAui - IMLHLuxADssyhve2vEALJr8KSbUAQAsa1cItPhBAD8POqw85FEALJoUvY15YQBA7D8sD/c5A - DSUayysMjEAMO2Quh638QA0WmUgAlIBACxKZXXIQsEALF40nuPLAQAseLsVhXcBADJ23rmyz - mEAQIARDWDzAQAsrA2aD7MhAE1R95vJwGEALdjZ68aYYQAxSR/sGCkhADyevfroaBEALmxSA - eStQQA7m3aUwM2hADVTP19jbWEAM4rONukOUQBHx9fU/2FZAC8QvunAEmEALFx4TRHXUQA+Z - dSn1rdRADTf3NN2ngD/5Fjgy7IjZP7bHI1nATQa/p51R7PmwGD+sHH7uZvEqP6/qoVFau3K/ - qA5Hs2wKJj+VYzdWUCsiP7mzNB+LFFy/vdDQDBpFx7+VFeEkxRxeP7F7m5KVODU/pVACJbjF - Jj+jdMHsaW1mP65HVhP0fJm/s2b06ISqEb+kQJCWKzfqP6COIa9BZSQ/zvjmbLoMiD9TIc+d - D8Ugv5o1IjVvHMo/s1oN2T4FLD+3OGyWxNwrP7EIIW+h0fm/uyMfiLLHcb+0y/cCrI9DP8sJ - wPl6FQA/o6KWFvbG6j+0vX2EvGQYP7O2T6sntQo/qCNltsfvAgAABAIAAAH/AAAADQAAAAIA + YXRlAAQACQAAAAJTRQAAAP4AAAIOAAAAlsA00hUKopC7v/uiIMLHLvi//ssyhve2vL/7Jr8K + SbT4v/sa1cItPgi//8POqw85FL/7JoUvY15YwAA7D8sD/c6//SUayysMjL/8O2Quh638v/0W + mUgAlIC/+xKZXXIQkL/7F40nuPLQv/seLsVhXbC//J23rmyzkMAAIARDWDzAv/srA2aD7NTA + A1R95vJwGL/7djZ68aYcv/xSR/sGCki//yevfroaIL/7mxSAeStQv/7m3aUwM2i//VTP19jb + WL/84rONukOMwAHx9fU/2Fa/+8QvunAEmL/7Fx4TRHXQv/+ZdSn1rci//Tf3NN2neD/nAgcM + Ar3hP2mZjSnOPhE/ayNFdQtNbz9qIjO8mlf2P2onEdowf9Y/bFTXaJqRrD9rK1dC1mI/P2w9 + mmBlEWQ/aT33KWgVMz9pDE2N/9UDP2pdJT4k8hE/aouFsw7Gnz9qlz3OrlTAP2oc4HHNYSc/ + aRVFdDl95D9tDHOUskvNP2qrErFFyQA/dMRV+0sb3z9qnM+HJMU3P2kMYnIPPmQ/apILe568 + Iz9pniYJ+vRZP2pRL+B2abs/aVODDQyKDj9pJSlA+y+/P3Gae7w4vXo/aoHzLWAumD9qNLl+ + 2CDgP2sSjLsR9GA/al5k0IGIuEAmg6Qr+gHCP5I5hg6ABgA/lTqBKHTFAD+SIcNl86EAP5If + omEDfwA/lvi9WmUhAD+THZI/2f4AP5mPg19ObgA/ktItKsZkAD+SjsiR3esAP5KkmoaHBgA/ + khZKaRjvAD+SFGHT2/UAP5Ig3ZVNFAA/kpwrkaAFAD+YCFgsLCIAP5IPS9yztQA/qlFcWiEQ + gD+RznOE6nQAP5KPYOKjMwA/lmRIfQIoAD+SN2h/dggAP5Xg9zNtrgA/ku5c/6I+AD+SsXhR + IusAP6MUf55iiAA/kd0sR7uFAD+SH4L0gRQAP5dfNWv49wA/kshKi4uHAEBE0hUKopC7QAui + IMLHLvhADssyhve2vEALJr8KSbT4QAsa1cItPghAD8POqw85FEALJoUvY15YQBA7D8sD/c5A + DSUayysMjEAMO2Quh638QA0WmUgAlIBACxKZXXIQkEALF40nuPLQQAseLsVhXbBADJ23rmyz + kEAQIARDWDzAQAsrA2aD7NRAE1R95vJwGEALdjZ68aYcQAxSR/sGCkhADyevfroaIEALmxSA + eStQQA7m3aUwM2hADVTP19jbWEAM4rONukOMQBHx9fU/2FZAC8QvunAEmEALFx4TRHXQQA+Z + dSn1rchADTf3NN2neD/5Fjgy7IhNP7bHI1nATNu/p51R7PmvmD+sHH7uZvExP6/qoVFauie/ + qA5Hs2wLOD+VYzdWUCoGP7mzNB+LFPm/vdDQDBpFTr+VFeEkxR0uP7F7m5KVOH4/pVACJbjE + sz+jdMHsaW26P65HVhP0fb2/s2b06ISpO7+kQJCWKzeLP6COIa9BZf0/zvjmbLoMHD9TIc+d + D90iv5o1IjVvGpE/s1oN2T4Ezj+3OGyWxNxQP7EIIW+h0Zu/uyMfiLLHPL+0y/cCrI2fP8sJ + wPl6FPE/o6KWFvbGlj+0vX2EvGOhP7O2T6sntYs/qCNltsfv5AAABAIAAAH/AAAADQAAAAIA AAAeAAAABQAABAIAAAL/AAAAEwAAAAIAAAD+AAAAEAAAAAUABAAJAAAACGVscGRfbG9vAAQA CQAAAA1tY3NlX2VscGRfbG9vAAQACQAAAAVwX2xvbwAEAAkAAAAFbG9vaWMABAAJAAAAEmlu - Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAev+M6GYBVASq/vGuV8SI2X7+n - nVHs+bAYv6TW5Y0JFQu/kWww4vFX5L+oDkezbAomP5VjN1ZQKyI/lNht5+qt4r+90NAMGkXH - v5UV4STFHF6/oHuKtB6KiL9jFZ1qwuogP3lA+cNvwvi/nu9VS24Wgr+zZvTohKoRv6RAkJYr - N+o/myfZhL4yPT+5EKUE8xZgv4E3FDJaGYy/mjUiNW8cyr8zDtdI8zuAv7tWWF8xxVc/hs8u - D/wJuL+7Ix+Issdxv7TL9wKsj0M/sGq450lfB7+1y0lRJZXnv6J/sQ9joN0/n+DuBpUSPj9T - FE1IlvogAAAADgAAAB5AA4tpeQZcIkCdJ778gk1wQJ6TSnvdXnBAnD6Pn+iHVUCcPOYri38u - QJ6D/vF10JpAmzztybeSzECefRV6RttJQJ6p/nbR4JNAnqxBiqhjm0Cb0Jv6ZpKMQJtvTMQm - Lw9Am1O/s5Mn+0CcUUc/BUohQJ6ryQO6wvNAnnp5yf77kUCbHdP+dPZgQJxwFCHevvhAmt1+ - OmGOAECerENzxOK5QJ6kbKmimTZAnSAFYnwMUUCepBjV4PmdQJ6pF6rZX3pAnqslCMgZy0Cd - UOL9AXOaQJr484ME1q1AnCOQo84P50CendeyiBCCQJvx7is3eLQAAAAOAAAAHj/wAAAAAAAA + Zmx1ZW5jZV9wYXJldG9fawAAAP4AAAITAAAAAwAAAA4AAAAev+M6GYBUtuK/vGuV8SEFY7+n + nVHs+a+Yv6TW5Y0ItKu/kWww4vIr3r+oDkezbAs4P5VjN1ZQKgY/lNht5+Gb3L+90NAMGkVO + v5UV4STFHS6/oHuKtB6H+79jFZ1pm6peP3lA+cM1Zw2/nu9VSy318r+zZvTohKk7v6RAkJYr + N4s/myfZhKlCvz+5EKUE8jDyv4E3FDJY90a/mjUiNW8akb8zDtdIUiGLv7tWWF8w/Zo/hs8u + D/kEfb+7Ix+Issc8v7TL9wKsjZ8/sGq450lfLL+1y0lRJZaqv6J/sQ9lzDo/n+DuBoWlbD9T + FE1Ilwn2AAAADgAAAB5AA4tpeQZbA0CdJ778gk1OQJ6TSnvdXnBAnD6Pn+iHVkCcPOYri39h + QJ6D/vF10JpAmzztybeSzECefRV6RtsvQJ6p/nbR4JNAnqxBiqhjm0Cb0Jv6ZpKMQJtvTMQm + L0pAm1O/s5Mn4kCcUUc/BUoAQJ6ryQO6wvNAnnp5yf77kUCbHdP+dPZRQJxwFCHevvJAmt1+ + OmGOCECerENzxOK5QJ6kbKmimNNAnSAFYnwMPECepBjV4PmnQJ6pF6rZX3pAnqslCMgZy0Cd + UOL9AXOaQJr484ME1q1AnCOQo84P8UCendeyiBCHQJvx7is3eLQAAAAOAAAAHj/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAA AAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAA AAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/w AAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAAP/AAAAAAAAA/8AAAAAAAAD/wAAAAAAAA P/AAAAAAAAA/8AAAAAAAAAAABAIAAAABAAQACQAAAAVuYW1lcwAAABAAAAADAAQACQAAAAhw YXJldG9fawAEAAkAAAAFbl9lZmYABAAJAAAABXJfZWZmAAAA/gAAAP4AAAAOAAAAAcBSi3CK - 9lOpAAAADgAAAAFAJ7mJXVQ1mAAAAA4AAAABQGKLcIr2U6kAAAAOAAAAAUAzARRFesUNAAAA - DgAAAAFAJnj9Bj3oVwAAAA4AAAABQEMBFEV6xQ0AAAQCAAAD/wAAABAAAAAKAAQACQAAAAll + 9lOnAAAADgAAAAFAJ7mJXVQ1jgAAAA4AAAABQGKLcIr2U6cAAAAOAAAAAUAzARRFesUIAAAA + DgAAAAFAJnj9Bj3oUAAAAA4AAAABQEMBFEV6xQgAAAQCAAAD/wAAABAAAAAKAAQACQAAAAll c3RpbWF0ZXMABAAJAAAACXBvaW50d2lzZQAEAAkAAAALZGlhZ25vc3RpY3MABAAJAAAAC3Bz aXNfb2JqZWN0AAQACQAAAAhlbHBkX2xvbwAEAAkAAAAFcF9sb28ABAAJAAAABWxvb2ljAAQA CQAAAAtzZV9lbHBkX2xvbwAEAAkAAAAIc2VfcF9sb28ABAAJAAAACHNlX2xvb2ljAAAEAgAA diff --git a/tests/testthat/test_loo_moment_matching.R b/tests/testthat/test_loo_moment_matching.R index b444a019..f7b18f07 100644 --- a/tests/testthat/test_loo_moment_matching.R +++ b/tests/testthat/test_loo_moment_matching.R @@ -286,6 +286,11 @@ test_that("loo_moment_match.default works", { loo_moment_match_object$pointwise[, "influence_pareto_k"], loo_manual$diagnostics$pareto_k ) + expect_equal( + unname(loo_moment_match_object$pointwise[1, "mcse_elpd_loo"]), + 0.1904162, + tolerance = 1e-6 + ) expect_snapshot_value(loo_moment_match_object, style = "serialize") From ea1cac446015e099dd9a2f6d7c00fc160fbe8321 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 22:18:17 +0300 Subject: [PATCH 10/17] fix: handle infinite log inputs --- R/importance_sampling.R | 12 +++++-- R/loo.R | 3 ++ R/loo_model_weights.R | 8 ++++- R/psis_approximate_posterior.R | 18 +++++++--- R/waic.R | 12 +++++-- tests/testthat/test_0_helpers.R | 9 +++++ tests/testthat/test_E_loo.R | 10 ++++++ tests/testthat/test_loo_and_waic.R | 35 +++++++++++++++++++ tests/testthat/test_loo_moment_matching.R | 23 +++++++++++- tests/testthat/test_model_weighting.R | 22 ++++++++++++ .../test_psis_approximate_posterior.R | 23 ++++++++++++ tests/testthat/test_relative_eff.R | 28 +++++++++++++++ tests/testthat/test_tisis.R | 25 +++++++++++++ 13 files changed, 218 insertions(+), 10 deletions(-) diff --git a/R/importance_sampling.R b/R/importance_sampling.R index f8dfd4d3..b85ad0dc 100644 --- a/R/importance_sampling.R +++ b/R/importance_sampling.R @@ -12,6 +12,14 @@ importance_sampling <- function(log_ratios, method, ...) { UseMethod("importance_sampling") } +validate_log_ratios <- function(x) { + validate_ll(x) + if (any(colSums(is.finite(x)) == 0)) { + stop("Each column of log ratios must contain at least one finite value.") + } + invisible(x) +} + #' @rdname importance_sampling #' @inheritParams psis @@ -24,8 +32,8 @@ importance_sampling.array <- cores <- loo_cores(cores) stopifnot(length(dim(log_ratios)) == 3) assert_importance_sampling_method_is_implemented(method) - log_ratios <- validate_ll(log_ratios) log_ratios <- llarray_to_matrix(log_ratios) + log_ratios <- validate_log_ratios(log_ratios) r_eff <- prepare_psis_r_eff(r_eff, len = ncol(log_ratios)) do_importance_sampling(log_ratios, r_eff = r_eff, cores = cores, method = method) } @@ -40,7 +48,7 @@ importance_sampling.matrix <- cores = getOption("mc.cores", 1)) { cores <- loo_cores(cores) assert_importance_sampling_method_is_implemented(method) - log_ratios <- validate_ll(log_ratios) + log_ratios <- validate_log_ratios(log_ratios) r_eff <- prepare_psis_r_eff(r_eff, len = ncol(log_ratios)) do_importance_sampling(log_ratios, r_eff = r_eff, cores = cores, method = method) } diff --git a/R/loo.R b/R/loo.R index 3833e04d..30773532 100644 --- a/R/loo.R +++ b/R/loo.R @@ -496,6 +496,9 @@ mcse_elpd <- function(ll, lw, E_elpd, r_eff, n_samples = NULL) { seq_len(ncol(lw)), FUN.VALUE = numeric(1), FUN = function(i) { + if (is.infinite(E_elpd[i]) && E_elpd[i] < 0) { + return(NA_real_) + } # Numerically stable way to compute # 1) variance in linear scale. Equation (6) in Vehtari et al. (2024) # 2) variance in log scale by matching the variance of a log-normal diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index d22ecb6c..4c005cd6 100644 --- a/R/loo_model_weights.R +++ b/R/loo_model_weights.R @@ -257,6 +257,9 @@ stacking_weights <- if (K < 2) { stop("At least two models are required for stacking weights.") } + if (any(rowSums(is.finite(lpd_point)) == 0)) { + stop("Each observation must have a finite predictive density for at least one model.") + } negative_log_score_loo <- function(w) { # objective function: log score @@ -321,9 +324,12 @@ pseudobma_weights <- if (K < 2) { stop("At least two models are required for pseudo-BMA weights.") } + elpd <- colSums2(lpd_point) + if (!any(is.finite(elpd))) { + stop("At least one model must have a finite total predictive density.") + } if (!BB) { - elpd <- colSums2(lpd_point) uwts <- exp(elpd - max(elpd)) wts <- structure( uwts / sum(uwts), diff --git a/R/psis_approximate_posterior.R b/R/psis_approximate_posterior.R index ebd9b4a4..37f61dfa 100644 --- a/R/psis_approximate_posterior.R +++ b/R/psis_approximate_posterior.R @@ -34,10 +34,11 @@ psis_approximate_posterior <- function(log_p = NULL, log_g = NULL, log_liks = NU checkmate::assert_flag(save_psis) if (is.null(log_liks)) { - approx_correction <- log_p - log_g - # Handle underflow/overflow - approx_correction <- approx_correction - max(approx_correction) + approx_correction <- validate_approx_correction(log_p, log_g) log_ratios <- matrix(approx_correction, ncol = 1) + log_ratios <- validate_log_ratios(log_ratios) + # Handle underflow/overflow + log_ratios <- log_ratios - max(log_ratios) } else { log_ratios <- correct_log_ratios(log_ratios = -log_liks, log_p = log_p, log_g = log_g) } @@ -65,9 +66,18 @@ psis_approximate_posterior <- function(log_p = NULL, log_g = NULL, log_liks = NU #' @inheritParams ap_psis #' @noRd #' @keywords internal -correct_log_ratios <- function(log_ratios, log_p, log_g) { +validate_approx_correction <- function(log_p, log_g) { approx_correction <- log_p - log_g + if (any(is.nan(approx_correction))) { + stop("The log density ratio is undefined for one or more draws.") + } + approx_correction +} + +correct_log_ratios <- function(log_ratios, log_p, log_g) { + approx_correction <- validate_approx_correction(log_p, log_g) log_ratios <- log_ratios + approx_correction + log_ratios <- validate_log_ratios(log_ratios) # Handle underflow/overflow log_ratio_max <- apply(log_ratios, 2, max) log_ratios <- sweep(log_ratios, MARGIN = 2, STATS = log_ratio_max) diff --git a/R/waic.R b/R/waic.R index 2dac6456..ab15b940 100644 --- a/R/waic.R +++ b/R/waic.R @@ -64,6 +64,14 @@ waic <- function(x, ...) { UseMethod("waic") } +validate_waic_log_lik <- function(x) { + validate_ll(x) + if (any(x == -Inf)) { + stop("All log-likelihood values must be finite for WAIC.") + } + invisible(x) +} + #' @export #' @templateVar fn waic #' @template array @@ -77,7 +85,7 @@ waic.array <- function(x, ...) { #' @template matrix #' waic.matrix <- function(x, ...) { - ll <- validate_ll(x) + ll <- validate_waic_log_lik(x) lldim <- dim(ll) lpd <- matrixStats::colLogSumExps(ll) - log(nrow(ll)) # colLogMeanExps p_waic <- matrixStats::colVars(ll) @@ -108,7 +116,7 @@ waic.function <- S <- length(as.vector(.llfun(data_i = data[1,, drop=FALSE], draws = draws, ...))) waic_list <- lapply(seq_len(N), FUN = function(i) { ll_i <- .llfun(data_i = data[i,, drop=FALSE], draws = draws, ...) - ll_i <- as.vector(ll_i) + ll_i <- validate_waic_log_lik(as.vector(ll_i)) lpd_i <- logMeanExp(ll_i) p_waic_i <- var(ll_i) elpd_waic_i <- lpd_i - p_waic_i diff --git a/tests/testthat/test_0_helpers.R b/tests/testthat/test_0_helpers.R index e17f3ec1..a0add421 100644 --- a/tests/testthat/test_0_helpers.R +++ b/tests/testthat/test_0_helpers.R @@ -48,6 +48,15 @@ test_that("colLogMeanExps(x) = log(colMeans(exp(x))) ", { expect_equal(colLogMeanExps(LLmat), log(colMeans(exp(LLmat)))) }) +test_that("log-mean-exp helpers handle negative infinity", { + x <- c(-Inf, log(2), log(3)) + expect_equal(logMeanExp(x), log(5 / 3)) + expect_equal(logMeanExp(rep(-Inf, 3)), -Inf) + + x <- cbind(x, rep(-Inf, 3)) + expect_equal(unname(colLogMeanExps(x)), c(log(5 / 3), -Inf)) +}) + test_that("validating log-lik objects and functions works", { f_ok <- function(data_i, draws) return(NULL) f_bad1 <- function(data_i) return(NULL) diff --git a/tests/testthat/test_E_loo.R b/tests/testthat/test_E_loo.R index 4837c860..e3e54e36 100644 --- a/tests/testthat/test_E_loo.R +++ b/tests/testthat/test_E_loo.R @@ -224,6 +224,16 @@ test_that("weighted quantiles work", { ) }) +test_that("E_loo handles negative infinite log ratios", { + log_ratios <- c(-Inf, seq(-9, 0, length.out = 99)) + psis_object <- suppressWarnings(psis(log_ratios)) + x <- seq_along(log_ratios) + + expect_no_error(out <- E_loo(x, psis_object, log_ratios = log_ratios)) + expect_equal(out$value, sum(weights(psis_object, log = FALSE) * x)) +}) + + test_that("weighted variance works", { x <- rnorm(100) w <- rep(0.01, 100) diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R index d7e50e6c..c73290a2 100644 --- a/tests/testthat/test_loo_and_waic.R +++ b/tests/testthat/test_loo_and_waic.R @@ -46,6 +46,41 @@ test_that("mcse_elpd is stable for extreme log likelihoods", { expect_equal(mcse_elpd(ll, lw, E_elpd, r_eff = 1), expected) }) +test_that("mcse_elpd returns NA for an all-zero likelihood column", { + ll <- cbind(c(-Inf, -Inf), c(-1, -2)) + lw <- matrix(log(0.5), nrow = 2, ncol = 2) + E_elpd <- matrixStats::colLogSumExps(ll + lw) + + out <- mcse_elpd(ll, lw, E_elpd, r_eff = 1) + expect_true(is.na(out[1])) + expect_true(is.finite(out[2])) +}) + +test_that("elpd handles negative infinite log likelihoods", { + log_lik <- cbind(c(-Inf, 0), c(-Inf, -Inf)) + out <- elpd(log_lik) + + expect_equal(out$pointwise[, "elpd"], c(-log(2), -Inf)) +}) + + +test_that("waic rejects negative infinite log likelihoods", { + log_lik <- matrix(-1, nrow = 10, ncol = 2) + log_lik[1, 1] <- -Inf + error <- "All log-likelihood values must be finite for WAIC." + + expect_error(waic(log_lik), error, fixed = TRUE) + expect_error(waic(array(log_lik, dim = c(5, 2, 2))), error, fixed = TRUE) + + llfun <- function(data_i, draws) log_lik[, data_i$i] + expect_error( + waic(llfun, data = data.frame(i = 1:2), draws = matrix(0, 10, 1)), + error, + fixed = TRUE + ) +}) + + test_that("waic returns object with correct structure", { expect_true(is.waic(waic1)) expect_true(is.loo(waic1)) diff --git a/tests/testthat/test_loo_moment_matching.R b/tests/testthat/test_loo_moment_matching.R index f7b18f07..58325a29 100644 --- a/tests/testthat/test_loo_moment_matching.R +++ b/tests/testthat/test_loo_moment_matching.R @@ -249,7 +249,7 @@ test_that("loo_moment_match.default works", { lwi_x <- lwi_1 lwi_x[which.min(lwi_1)] <- -Inf expect_no_error(suppressWarnings(importance_sampling.default( - lwi_1, + lwi_x, method = "psis", r_eff = 1, cores = 1 @@ -518,6 +518,27 @@ test_that("loo_moment_match_split works", { ) expect_snapshot_value(split2, style = "serialize") + + log_prob_with_inf <- function(x, upars, ...) { + out <- log_prob_upars_test(x, upars, ...) + out[1] <- -Inf + out + } + expect_no_error(split3 <- loo_moment_match_split( + x, + upars, + cov = FALSE, + total_shift = c(0, 0), + total_scaling = c(1, 1), + total_mapping = diag(c(1, 1)), + i = 1, + log_prob_upars = log_prob_with_inf, + log_lik_i_upars = log_lik_i_upars_test, + cores = 1, + r_eff_i = 1, + is_method = "psis" + )) + expect_false(anyNA(split3$lwi)) }) test_that("shift_and_scale is stable for large parameter values", { diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index 28b5cf0f..7e12b4a6 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -133,6 +133,28 @@ test_that("loo_model_weights (stacking and pseudo-BMA) gives expected result", { expect_identical(w3, w3_b) }) +test_that("pseudo-BMA gives zero weight to an impossible model", { + lpd <- cbind(rep(-Inf, 3), c(-2, -1, 0)) + + expect_equal(as.numeric(pseudobma_weights(lpd, BB = FALSE)), c(0, 1)) +}) + + +test_that("model weighting rejects inputs with no finite predictive density", { + stacking_lpd <- rbind(c(-Inf, -Inf), c(-1, -1)) + expect_error( + stacking_weights(stacking_lpd), + "Each observation must have a finite predictive density for at least one model.", + fixed = TRUE + ) + + pseudobma_lpd <- matrix(-Inf, nrow = 2, ncol = 2) + error <- "At least one model must have a finite total predictive density." + expect_error(pseudobma_weights(pseudobma_lpd, BB = FALSE), error, fixed = TRUE) + expect_error(pseudobma_weights(pseudobma_lpd, BB = TRUE), error, fixed = TRUE) +}) + + test_that("stacking_weights and pseudobma_weights throw correct errors", { xx <- cbind(rnorm(10)) expect_error(stacking_weights(xx), "two models are required") diff --git a/tests/testthat/test_psis_approximate_posterior.R b/tests/testthat/test_psis_approximate_posterior.R index 8d0ffe70..a46d1965 100644 --- a/tests/testthat/test_psis_approximate_posterior.R +++ b/tests/testthat/test_psis_approximate_posterior.R @@ -308,6 +308,29 @@ test_that("ADVI meanfield approximation, normal model", { }) +test_that("approximate posterior methods reject undefined log density ratios", { + log_p <- c(-Inf, rep(0, 9)) + log_g <- c(-Inf, rep(0, 9)) + error <- "The log density ratio is undefined for one or more draws." + + expect_error( + psis_approximate_posterior( + log_p = log_p, + log_g = log_g, + cores = 1, + save_psis = FALSE + ), + error, + fixed = TRUE + ) + expect_error( + ap_psis(matrix(seq(-9, 0, length.out = 10), ncol = 1), log_p, log_g), + error, + fixed = TRUE + ) +}) + + test_that("Deprecation of log_q argument", { log_p <- test_data_psis_approximate_posterior$laplace_independent$log_p log_g <- test_data_psis_approximate_posterior$laplace_independent$log_q diff --git a/tests/testthat/test_relative_eff.R b/tests/testthat/test_relative_eff.R index 195a2c47..b728d21a 100644 --- a/tests/testthat/test_relative_eff.R +++ b/tests/testthat/test_relative_eff.R @@ -8,6 +8,34 @@ test_that("relative_eff results haven't changed", { expect_snapshot_value(relative_eff(exp(LLarr)), style = "serialize") }) +test_that("relative_eff handles zero likelihood values", { + likelihood <- c(0, exp(seq(-9, 0, length.out = 99))) + + expect_equal( + relative_eff(likelihood, chain_id = rep(1, length(likelihood))), + posterior::ess_mean(likelihood) / length(likelihood) + ) +}) + +test_that("relative_eff returns NA for all-zero likelihoods", { + likelihood <- rep(0, 100) + chain_id <- rep(1, length(likelihood)) + expect_true(is.na(relative_eff(likelihood, chain_id = chain_id))) + + likelihood_matrix <- cbind(likelihood, likelihood) + expect_true(all(is.na(relative_eff(likelihood_matrix, chain_id = chain_id)))) + + zero_likelihood <- function(data_i, draws) rep(0, nrow(draws)) + out <- relative_eff( + zero_likelihood, + chain_id = chain_id, + data = data.frame(i = 1:2), + draws = matrix(0, nrow = length(likelihood), ncol = 1), + cores = 1 + ) + expect_true(all(is.na(out))) +}) + test_that("relative_eff is equal to ESS / S", { dims <- dim(LLarr) ess <- r_eff <- rep(NA, dims[3]) diff --git a/tests/testthat/test_tisis.R b/tests/testthat/test_tisis.R index edfcbd3f..eecd7008 100644 --- a/tests/testthat/test_tisis.R +++ b/tests/testthat/test_tisis.R @@ -137,6 +137,31 @@ test_that("tis throws correct errors and warnings", { }) +test_that("importance sampling methods handle negative infinite log ratios", { + log_ratios <- c(-Inf, seq(-9, 0, length.out = 99)) + methods <- list(psis = psis, tis = tis, sis = sis) + + lapply(methods, function(method) { + out <- suppressWarnings(method(log_ratios, r_eff = NA)) + expect_identical(weights(out, normalize = FALSE, log = TRUE)[1], -Inf) + expect_identical(weights(out, normalize = TRUE, log = TRUE)[1], -Inf) + expect_identical(weights(out, normalize = TRUE, log = FALSE)[1], 0) + expect_equal(sum(weights(out, normalize = TRUE, log = FALSE)), 1) + }) +}) + +test_that("importance sampling methods reject all negative infinite columns", { + methods <- list(psis = psis, tis = tis, sis = sis) + log_ratios <- cbind(seq(-9, 0, length.out = 10), rep(-Inf, 10)) + error <- "Each column of log ratios must contain at least one finite value." + + lapply(methods, function(method) { + expect_error(method(rep(-Inf, 10), r_eff = NA), error, fixed = TRUE) + expect_error(method(log_ratios, r_eff = NA), error, fixed = TRUE) + }) +}) + + test_that("explict test of values for 'sis' and 'tis'", { lw <- 1:16 expect_silent(tis_true <- tis(log_ratios = lw, r_eff = NA)) From 912929c494650762dc3c58463b6de436ad67038f Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 22:23:07 +0300 Subject: [PATCH 11/17] fix: retain tiny positive MCSE deviations --- R/loo.R | 2 +- tests/testthat/test_loo_and_waic.R | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/R/loo.R b/R/loo.R index 30773532..17837fd2 100644 --- a/R/loo.R +++ b/R/loo.R @@ -507,7 +507,7 @@ mcse_elpd <- function(ll, lw, E_elpd, r_eff, n_samples = NULL) { positive <- log_lik_ratio > 0 log_abs_diff <- numeric(length(log_lik_ratio)) log_abs_diff[positive] <- - log_lik_ratio[positive] + log1p(-exp(-log_lik_ratio[positive])) + log_lik_ratio[positive] + log(-expm1(-log_lik_ratio[positive])) log_abs_diff[!positive] <- log(-expm1(log_lik_ratio[!positive])) log_var_epd_ratio <- matrixStats::logSumExp(2 * lw[, i] + 2 * log_abs_diff) - log(r_eff[i]) diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R index c73290a2..dd74ebce 100644 --- a/tests/testthat/test_loo_and_waic.R +++ b/tests/testthat/test_loo_and_waic.R @@ -46,6 +46,16 @@ test_that("mcse_elpd is stable for extreme log likelihoods", { expect_equal(mcse_elpd(ll, lw, E_elpd, r_eff = 1), expected) }) +test_that("mcse_elpd retains tiny positive deviations", { + z <- 5e-17 + ll <- matrix(c(-z, z), ncol = 1) + lw <- matrix(log(0.5), nrow = 2, ncol = 1) + expected <- sqrt(log1p(sum(exp(lw)^2 * expm1(ll)^2))) + out <- mcse_elpd(ll, lw, E_elpd = 0, r_eff = 1) + + expect_equal(out / expected, 1, tolerance = 1e-12) +}) + test_that("mcse_elpd returns NA for an all-zero likelihood column", { ll <- cbind(c(-Inf, -Inf), c(-1, -2)) lw <- matrix(log(0.5), nrow = 2, ncol = 2) From 9ea4fb59c6e97597aa347dd0599150af102259f7 Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:15:44 +0300 Subject: [PATCH 12/17] perf: vectorize mcse_elpd and cheapen log-ratio validation --- R/helpers.R | 4 +-- R/importance_sampling.R | 4 ++- R/loo.R | 79 +++++++++++++++++++++++++++-------------- 3 files changed, 57 insertions(+), 30 deletions(-) diff --git a/R/helpers.R b/R/helpers.R index a11553b6..0deb140a 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -36,8 +36,8 @@ colLogMeanExps <- function(x) { #' exp_x_minus_exp_y <- function(x, y) { out <- -exp(x) * expm1(y - x) - equal <- x == y - out[!is.na(equal) & equal] <- 0 + # which() drops the NA comparisons that NA or NaN inputs would produce + out[which(x == y)] <- 0 out } diff --git a/R/importance_sampling.R b/R/importance_sampling.R index b85ad0dc..46e61662 100644 --- a/R/importance_sampling.R +++ b/R/importance_sampling.R @@ -14,7 +14,9 @@ importance_sampling <- function(log_ratios, method, ...) { validate_log_ratios <- function(x) { validate_ll(x) - if (any(colSums(is.finite(x)) == 0)) { + # validate_ll() has already ruled out NA and +Inf, so a column without a + # finite value is exactly a column whose maximum is -Inf + if (any(matrixStats::colMaxs(x) == -Inf)) { stop("Each column of log ratios must contain at least one finite value.") } invisible(x) diff --git a/R/loo.R b/R/loo.R index 17837fd2..e11f3466 100644 --- a/R/loo.R +++ b/R/loo.R @@ -488,39 +488,64 @@ importance_sampling_loo_object <- function(pointwise, diagnostics, dims, #' @return Vector of standard error estimates. #' mcse_elpd <- function(ll, lw, E_elpd, r_eff, n_samples = NULL) { - if (length(r_eff) == 1 && !is.null(ncol(ll))) { + if (!is.matrix(ll)) { + ll <- as.matrix(ll) + } + if (!is.matrix(lw)) { + lw <- as.matrix(lw) + } + S <- nrow(ll) + if (length(r_eff) == 1) { r_eff <- rep(r_eff, ncol(ll)) } - var_elpd <- - vapply( - seq_len(ncol(lw)), - FUN.VALUE = numeric(1), - FUN = function(i) { - if (is.infinite(E_elpd[i]) && E_elpd[i] < 0) { - return(NA_real_) - } - # Numerically stable way to compute - # 1) variance in linear scale. Equation (6) in Vehtari et al. (2024) - # 2) variance in log scale by matching the variance of a log-normal - # https://en.wikipedia.org/wiki/Log-normal_distribution#Arithmetic_moments - log_lik_ratio <- ll[, i] - E_elpd[i] - positive <- log_lik_ratio > 0 - log_abs_diff <- numeric(length(log_lik_ratio)) - log_abs_diff[positive] <- - log_lik_ratio[positive] + log(-expm1(-log_lik_ratio[positive])) - log_abs_diff[!positive] <- log(-expm1(log_lik_ratio[!positive])) - log_var_epd_ratio <- - matrixStats::logSumExp(2 * lw[, i] + 2 * log_abs_diff) - log(r_eff[i]) - if (log_var_epd_ratio > 0) { - log_var_epd_ratio + log1p(exp(-log_var_epd_ratio)) - } else { - log1p(exp(log_var_epd_ratio)) - } - } + # Everything is computed relative to the loo predictive density, so that + # 1) exp() of the log likelihood never over- or underflows, and + # 2) expm1() avoids the cancellation in `exp(ll) - exp(E_elpd)`. + # `ll - E_elpd` is bounded above by `-lw`, so the product below cannot + # overflow for consistent (ll, lw, E_elpd); the fallback covers the rest. + # + # Variance in linear scale, relative to E_epd^2. + # Equation (6) in Vehtari et al. (2024) + var_epd_ratio <- + matrixStats::colSums2((exp(lw) * expm1(ll - rep(E_elpd, each = S)))^2) / + r_eff + # Variance in log scale by matching the variance of a log-normal + # https://en.wikipedia.org/wiki/Log-normal_distribution#Arithmetic_moments + var_elpd <- log1p(var_epd_ratio) + undefined <- is.infinite(E_elpd) & E_elpd < 0 + overflow <- !is.finite(var_epd_ratio) & !undefined + if (any(overflow)) { + lvr <- log_var_epd_ratio( + ll[, overflow, drop = FALSE] - rep(E_elpd[overflow], each = S), + lw[, overflow, drop = FALSE], + r_eff[overflow] ) + var_elpd[overflow] <- + ifelse(lvr > 0, lvr + log1p(exp(-lvr)), log1p(exp(lvr))) + } + var_elpd[undefined] <- NA_real_ sqrt(var_elpd) } +#' Log of the relative linear-scale ELPD variance, for the rare case where +#' `exp(lw) * expm1(log_lik_ratio)` over- or underflows +#' +#' @noRd +#' @param log_lik_ratio Matrix of `ll - E_elpd` values. +#' @param lw Matrix of normalized log weights. +#' @param r_eff Vector of relative effective sample sizes. +#' @return Vector of `log(var_epd / E_epd^2)` values. +#' +log_var_epd_ratio <- function(log_lik_ratio, lw, r_eff) { + log_abs_diff <- log(abs(expm1(log_lik_ratio))) + big <- which(log_lik_ratio > 700) + if (length(big)) { + log_abs_diff[big] <- + log_lik_ratio[big] + log(-expm1(-log_lik_ratio[big])) + } + matrixStats::colLogSumExps(2 * (lw + log_abs_diff)) - log(r_eff) +} + #' Warning message if r_eff not specified #' @noRd From 374b8fae43940b9c873e4b4469e3bbd094eab90a Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:16:56 +0300 Subject: [PATCH 13/17] fix: propagate missing values through exp_diff_over_exp --- R/helpers.R | 23 ++++++++++++----------- tests/testthat/test_model_weighting.R | 13 +++++++++++++ 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/R/helpers.R b/R/helpers.R index 0deb140a..5b0ec517 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -55,19 +55,20 @@ difference_of_squares <- function(x, y) { #' #' @noRd #' @param a,b,c Numeric vectors of the same length. -#' @return A numeric vector equal to `(exp(a) - exp(b)) / exp(c)`. +#' @return A numeric vector equal to `(exp(a) - exp(b)) / exp(c)`. Elements +#' with `a == b` are returned as an exact zero regardless of `c`; elsewhere +#' `NA` and `NaN` inputs propagate. #' exp_diff_over_exp <- function(a, b, c) { - a_is_larger <- a >= b - out <- numeric(length(a)) - out[a_is_larger] <- - exp(a[a_is_larger] - c[a_is_larger]) * - -expm1(b[a_is_larger] - a[a_is_larger]) - out[!a_is_larger] <- - exp(b[!a_is_larger] - c[!a_is_larger]) * - expm1(a[!a_is_larger] - b[!a_is_larger]) - equal <- a == b - out[!is.na(equal) & equal] <- 0 + # `a >= b` is NA if `a` or `b` is NA or NaN, and R silently ignores NA + # indices in `[<-`. Seed the result from the inputs and index with which() + # so that missing values propagate instead of leaving a zero behind. + out <- a + b + c + larger <- which(a >= b) + smaller <- which(a < b) + out[larger] <- exp(a[larger] - c[larger]) * -expm1(b[larger] - a[larger]) + out[smaller] <- exp(b[smaller] - c[smaller]) * expm1(a[smaller] - b[smaller]) + out[which(a == b)] <- 0 out } diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index 7e12b4a6..4a03a930 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -32,6 +32,19 @@ test_that("stacking gradient is stable for similar model predictions", { ) }) +test_that("exp_diff_over_exp propagates missing values", { + expect_equal(exp_diff_over_exp(NaN, 0, 0), NaN) + expect_equal(exp_diff_over_exp(0, NaN, 0), NaN) + expect_equal(exp_diff_over_exp(0, -1, NaN), NaN) + expect_equal(exp_diff_over_exp(NA_real_, 0, 0), NA_real_) + # equal numerators short-circuit to an exact zero, whatever the denominator + expect_equal(exp_diff_over_exp(0, 0, NaN), 0) + expect_equal( + exp_diff_over_exp(c(NaN, 0), c(0, -1), c(0, 0)), + c(NaN, exp_diff_over_exp(0, -1, 0)) + ) +}) + test_that("stacking handles equal infinite log predictive densities", { lpd <- matrix( c(-Inf, 0, -Inf, -1, -1, -1), From bf22710bb2d10f8b148a68d2768739374641f83b Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:18:56 +0300 Subject: [PATCH 14/17] fix: report -Inf log-likelihood errors in terms of the loo() input --- R/helpers.R | 28 +++++++++++++++++++++++++++- R/loo.R | 10 +++++++--- tests/testthat/test_loo_and_waic.R | 21 +++++++++++++++++++++ 3 files changed, 55 insertions(+), 4 deletions(-) diff --git a/R/helpers.R b/R/helpers.R index 5b0ec517..e46831bc 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -104,11 +104,37 @@ validate_ll <- function(x) { } else if (anyNA(x)) { stop("NAs not allowed in input.") } else if (any(x == Inf)) { - stop("All input values must be finite or -Inf.") + # classed so that callers which negate a log-likelihood matrix can report + # the error in terms of the input the user actually supplied + stop(errorCondition( + "All input values must be finite or -Inf.", + class = "loo_positive_infinity_error" + )) } invisible(x) } +#' Report `+Inf` log ratios in terms of the log likelihood that produced them +#' +#' `loo()` negates the log-likelihood before importance sampling, so a `-Inf` +#' log-likelihood value reaches [validate_ll()] as `+Inf`. Wrap the importance +#' sampling call so the user sees a message about their own input. The wrapped +#' expression is only forced inside the handler, so there is no cost unless an +#' error is raised. +#' +#' @noRd +#' @param expr Expression that negates a log-likelihood and importance samples it. +#' @return The value of `expr`. +#' +with_log_lik_error_message <- function(expr) { + withCallingHandlers( + expr, + loo_positive_infinity_error = function(cnd) { + stop("-Inf log-likelihood values are not allowed.", call. = FALSE) + } + ) +} + #' Convert iter by chain by obs array to (iter * chain) by obs matrix #' #' @noRd diff --git a/R/loo.R b/R/loo.R index e11f3466..af9fbb7d 100644 --- a/R/loo.R +++ b/R/loo.R @@ -198,7 +198,9 @@ loo.array <- cores = getOption("mc.cores", 1), is_method = c("psis", "tis", "sis")) { is_method <- match.arg(is_method) - psis_out <- importance_sampling.array(log_ratios = -x, r_eff = r_eff, cores = cores, method = is_method) + psis_out <- with_log_lik_error_message( + importance_sampling.array(log_ratios = -x, r_eff = r_eff, cores = cores, method = is_method) + ) ll <- llarray_to_matrix(x) pointwise <- pointwise_loo_calcs(ll, psis_out) importance_sampling_loo_object( @@ -222,13 +224,14 @@ loo.matrix <- cores = getOption("mc.cores", 1), is_method = c("psis", "tis", "sis")) { is_method <- match.arg(is_method) - psis_out <- + psis_out <- with_log_lik_error_message( importance_sampling.matrix( log_ratios = -x, r_eff = r_eff, cores = cores, method = is_method ) + ) pointwise <- pointwise_loo_calcs(x, psis_out) importance_sampling_loo_object( pointwise = pointwise, @@ -371,13 +374,14 @@ loo_i <- if (!is.matrix(ll_i)) { ll_i <- as.matrix(ll_i) } - psis_out <- + psis_out <- with_log_lik_error_message( importance_sampling.matrix( log_ratios = -ll_i, r_eff = r_eff, cores = 1, method = is_method ) + ) structure( list( pointwise = pointwise_loo_calcs(ll_i, psis_out), diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R index dd74ebce..fdd11ce9 100644 --- a/tests/testthat/test_loo_and_waic.R +++ b/tests/testthat/test_loo_and_waic.R @@ -74,6 +74,27 @@ test_that("elpd handles negative infinite log likelihoods", { }) +test_that("loo rejects negative infinite log likelihoods", { + log_lik <- matrix(-1, nrow = 10, ncol = 2) + log_lik[1, 1] <- -Inf + error <- "-Inf log-likelihood values are not allowed." + + expect_error(loo(log_lik, r_eff = NA), error, fixed = TRUE) + expect_error( + loo(array(log_lik, dim = c(5, 2, 2)), r_eff = NA), + error, + fixed = TRUE + ) + + llfun <- function(data_i, draws) log_lik[, data_i$i] + expect_error( + loo(llfun, data = data.frame(i = 1:2), draws = matrix(0, 10, 1), + r_eff = NA, cores = 1), + error, + fixed = TRUE + ) +}) + test_that("waic rejects negative infinite log likelihoods", { log_lik <- matrix(-1, nrow = 10, ncol = 2) log_lik[1, 1] <- -Inf From 4d884ceb150420162082ddae08846554098742df Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:19:11 +0300 Subject: [PATCH 15/17] docs: note the normalized-weight assumption in the two-pass variances --- R/E_loo.R | 2 ++ R/loo_moment_matching.R | 2 ++ 2 files changed, 4 insertions(+) diff --git a/R/E_loo.R b/R/E_loo.R index b0a43f46..e360c3ef 100644 --- a/R/E_loo.R +++ b/R/E_loo.R @@ -219,6 +219,8 @@ E_loo.matrix <- # sample size ESS is estimated with the generic target quantity invariant # estimate 1/sum(w^2), see e.g. "Monte Carlo theory, methods and examples" # by Owen (2013). + # The two-pass form avoids the cancellation in E[x^2] - E[x]^2 and is + # equivalent to it only because `w` sums to one. weighted_mean <- .wmean(x, w) sum(w * (x - weighted_mean)^2) / (1 - sum(w^2)) } diff --git a/R/loo_moment_matching.R b/R/loo_moment_matching.R index c31a12bb..2d14e20b 100644 --- a/R/loo_moment_matching.R +++ b/R/loo_moment_matching.R @@ -534,6 +534,8 @@ shift_and_scale <- function(x, upars, lwi) { weights <- exp(lwi) mean_weighted <- colSums(weights * upars) shift <- mean_weighted - mean_original + # The two-pass form avoids the cancellation in E[x^2] - E[x]^2 and is + # equivalent to it only because `weights` sums to one. centered <- sweep(upars, 2, mean_weighted) mii <- colSums(weights * centered^2) mii <- mii * S / (S - 1) From 2d12c3ff78a6a772acd1d5c392f144f681024e0e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:20:59 +0300 Subject: [PATCH 16/17] style: match the surrounding blank-line spacing in the new tests --- tests/testthat/test_E_loo.R | 1 - tests/testthat/test_loo_and_waic.R | 2 -- tests/testthat/test_model_weighting.R | 2 -- tests/testthat/test_psis.R | 1 - 4 files changed, 6 deletions(-) diff --git a/tests/testthat/test_E_loo.R b/tests/testthat/test_E_loo.R index e3e54e36..913ca9aa 100644 --- a/tests/testthat/test_E_loo.R +++ b/tests/testthat/test_E_loo.R @@ -233,7 +233,6 @@ test_that("E_loo handles negative infinite log ratios", { expect_equal(out$value, sum(weights(psis_object, log = FALSE) * x)) }) - test_that("weighted variance works", { x <- rnorm(100) w <- rep(0.01, 100) diff --git a/tests/testthat/test_loo_and_waic.R b/tests/testthat/test_loo_and_waic.R index fdd11ce9..f32595f8 100644 --- a/tests/testthat/test_loo_and_waic.R +++ b/tests/testthat/test_loo_and_waic.R @@ -73,7 +73,6 @@ test_that("elpd handles negative infinite log likelihoods", { expect_equal(out$pointwise[, "elpd"], c(-log(2), -Inf)) }) - test_that("loo rejects negative infinite log likelihoods", { log_lik <- matrix(-1, nrow = 10, ncol = 2) log_lik[1, 1] <- -Inf @@ -111,7 +110,6 @@ test_that("waic rejects negative infinite log likelihoods", { ) }) - test_that("waic returns object with correct structure", { expect_true(is.waic(waic1)) expect_true(is.loo(waic1)) diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index 4a03a930..344d7254 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -152,7 +152,6 @@ test_that("pseudo-BMA gives zero weight to an impossible model", { expect_equal(as.numeric(pseudobma_weights(lpd, BB = FALSE)), c(0, 1)) }) - test_that("model weighting rejects inputs with no finite predictive density", { stacking_lpd <- rbind(c(-Inf, -Inf), c(-1, -1)) expect_error( @@ -167,7 +166,6 @@ test_that("model weighting rejects inputs with no finite predictive density", { expect_error(pseudobma_weights(pseudobma_lpd, BB = TRUE), error, fixed = TRUE) }) - test_that("stacking_weights and pseudobma_weights throw correct errors", { xx <- cbind(rnorm(10)) expect_error(stacking_weights(xx), "two models are required") diff --git a/tests/testthat/test_psis.R b/tests/testthat/test_psis.R index 183f868a..9e033715 100644 --- a/tests/testthat/test_psis.R +++ b/tests/testthat/test_psis.R @@ -171,7 +171,6 @@ test_that("exp_x_minus_exp_y is stable for nearby values", { ) }) - test_that("psis_smooth_tail returns original tail values if k is infinite", { xx <- log(c(2, 2, 2, 2, 3, 4, 5, 6)) val <- suppressWarnings(psis_smooth_tail(xx, 0)) From 32fad02c41df263f97d89af5cc3690577074015e Mon Sep 17 00:00:00 2001 From: Aki Vehtari Date: Mon, 14 Sep 2026 23:21:48 +0300 Subject: [PATCH 17/17] test: name the importance sampling method in the -Inf log-ratio tests --- tests/testthat/test_tisis.R | 40 ++++++++++++++++++++++--------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/tests/testthat/test_tisis.R b/tests/testthat/test_tisis.R index eecd7008..581d4898 100644 --- a/tests/testthat/test_tisis.R +++ b/tests/testthat/test_tisis.R @@ -137,29 +137,37 @@ test_that("tis throws correct errors and warnings", { }) -test_that("importance sampling methods handle negative infinite log ratios", { - log_ratios <- c(-Inf, seq(-9, 0, length.out = 99)) - methods <- list(psis = psis, tis = tis, sis = sis) +# One test_that() per method, so that a failure names the method that failed +for (method_name in c("psis", "tis", "sis")) { + test_that(paste0(method_name, "() handles negative infinite log ratios"), { + method <- match.fun(method_name) + zero <- c(TRUE, FALSE, TRUE, rep(FALSE, 97)) + log_ratios <- seq(-9, 0, length.out = 100) + log_ratios[zero] <- -Inf - lapply(methods, function(method) { out <- suppressWarnings(method(log_ratios, r_eff = NA)) - expect_identical(weights(out, normalize = FALSE, log = TRUE)[1], -Inf) - expect_identical(weights(out, normalize = TRUE, log = TRUE)[1], -Inf) - expect_identical(weights(out, normalize = TRUE, log = FALSE)[1], 0) - expect_equal(sum(weights(out, normalize = TRUE, log = FALSE)), 1) + w_log <- weights(out, normalize = TRUE, log = TRUE) + w <- weights(out, normalize = TRUE, log = FALSE) + + expect_identical(weights(out, normalize = FALSE, log = TRUE)[zero], c(-Inf, -Inf)) + expect_identical(w_log[zero], c(-Inf, -Inf)) + expect_identical(w[zero], c(0, 0)) + expect_true(all(w[!zero] > 0)) + expect_equal(sum(w), 1) }) -}) -test_that("importance sampling methods reject all negative infinite columns", { - methods <- list(psis = psis, tis = tis, sis = sis) - log_ratios <- cbind(seq(-9, 0, length.out = 10), rep(-Inf, 10)) - error <- "Each column of log ratios must contain at least one finite value." + test_that(paste0(method_name, "() rejects all negative infinite columns"), { + method <- match.fun(method_name) + error <- "Each column of log ratios must contain at least one finite value." - lapply(methods, function(method) { expect_error(method(rep(-Inf, 10), r_eff = NA), error, fixed = TRUE) - expect_error(method(log_ratios, r_eff = NA), error, fixed = TRUE) + expect_error( + method(cbind(seq(-9, 0, length.out = 10), rep(-Inf, 10)), r_eff = NA), + error, + fixed = TRUE + ) }) -}) +} test_that("explict test of values for 'sis' and 'tis'", {