diff --git a/R/E_loo.R b/R/E_loo.R index b449c182..e360c3ef 100644 --- a/R/E_loo.R +++ b/R/E_loo.R @@ -219,7 +219,10 @@ 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)) + # 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)) } .wsd <- function(x, w, ...) { sqrt(.wvar(x, w)) diff --git a/R/helpers.R b/R/helpers.R index 38b401dd..e46831bc 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -26,6 +26,52 @@ 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) { + out <- -exp(x) * expm1(y - x) + # which() drops the NA comparisons that NA or NaN inputs would produce + out[which(x == y)] <- 0 + out +} + +#' 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)`. 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 >= 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 +} + #' Compute point estimates and standard errors from pointwise vectors #' #' @noRd @@ -58,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/importance_sampling.R b/R/importance_sampling.R index f8dfd4d3..46e61662 100644 --- a/R/importance_sampling.R +++ b/R/importance_sampling.R @@ -12,6 +12,16 @@ importance_sampling <- function(log_ratios, method, ...) { UseMethod("importance_sampling") } +validate_log_ratios <- function(x) { + validate_ll(x) + # 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) +} + #' @rdname importance_sampling #' @inheritParams psis @@ -24,8 +34,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 +50,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 10b1bdc7..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), @@ -488,29 +492,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) { - lik <- exp(ll) - w2 <- exp(lw)^2 - E_epd <- exp(E_elpd) - 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(w2)), - 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 - # https://en.wikipedia.org/wiki/Log-normal_distribution#Arithmetic_moments - log(1 + var_epd_i / E_epd[i]^2) - } + # 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 diff --git a/R/loo_model_weights.R b/R/loo_model_weights.R index 946dc7c3..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 @@ -272,11 +275,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) } @@ -317,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/loo_moment_matching.R b/R/loo_moment_matching.R index 110eff93..2d14e20b 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, @@ -531,11 +531,14 @@ 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) + # 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) scaling <- sqrt(mii / matrixStats::colVars(upars)) # transform posterior draws upars_new <- sweep(upars, 2, mean_original, "-") diff --git a/R/loo_subsample.R b/R/loo_subsample.R index bcac4b17..4aad6d46 100644 --- a/R/loo_subsample.R +++ b/R/loo_subsample.R @@ -1185,7 +1185,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/R/psis.R b/R/psis.R index dd8b0176..303da755 100644 --- a/R/psis.R +++ b/R/psis.R @@ -254,7 +254,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/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/_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_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 a3bf5c81..913ca9aa 100644 --- a/tests/testthat/test_E_loo.R +++ b/tests/testthat/test_E_loo.R @@ -224,6 +224,15 @@ 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) @@ -232,4 +241,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_and_waic.R b/tests/testthat/test_loo_and_waic.R index 9bbd16ea..f32595f8 100644 --- a/tests/testthat/test_loo_and_waic.R +++ b/tests/testthat/test_loo_and_waic.R @@ -34,6 +34,82 @@ 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("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) + 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("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 + 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 2c756533..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 @@ -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") @@ -513,6 +518,35 @@ 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", { + 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", { 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 diff --git a/tests/testthat/test_model_weighting.R b/tests/testthat/test_model_weighting.R index 96c075b8..344d7254 100644 --- a/tests/testthat/test_model_weighting.R +++ b/tests/testthat/test_model_weighting.R @@ -20,6 +20,41 @@ 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) + expect_equal( + exp_diff_over_exp(c(-Inf, 0, Inf), c(-Inf, 0, Inf), c(0, 0, 0)), + c(0, 0, 0) + ) +}) + +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), + 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", { expect_error( loo_model_weights(log_lik1), @@ -111,6 +146,26 @@ 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.R b/tests/testthat/test_psis.R index 93d7501a..9e033715 100644 --- a/tests/testthat/test_psis.R +++ b/tests/testthat/test_psis.R @@ -152,6 +152,25 @@ 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) + ) +}) + 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)) 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..581d4898 100644 --- a/tests/testthat/test_tisis.R +++ b/tests/testthat/test_tisis.R @@ -137,6 +137,39 @@ test_that("tis throws correct errors and warnings", { }) +# 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 + + out <- suppressWarnings(method(log_ratios, r_eff = NA)) + 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(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." + + expect_error(method(rep(-Inf, 10), 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'", { lw <- 1:16 expect_silent(tis_true <- tis(log_ratios = lw, r_eff = NA))