From ba04ad413b6d1a24004931a20d65befa434549ff Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 09:46:55 +0000 Subject: [PATCH 01/13] Add function for constructing target distribution from formula --- R/bridges.R | 46 +++++++++++++++++++ ...t_distribution_from_log_density_formula.Rd | 31 +++++++++++++ 2 files changed, 77 insertions(+) create mode 100644 man/target_distribution_from_log_density_formula.Rd diff --git a/R/bridges.R b/R/bridges.R index d705ec1..d1d5022 100644 --- a/R/bridges.R +++ b/R/bridges.R @@ -111,3 +111,49 @@ example_gaussian_stan_model <- function(n_data = 50, seed = 1234L) { fileext = ".stan" ) } + + +#' Construct target distribution from a formula specifying log density. +#' +#' @param log_density_formula Formula for which right-hand side specifies +#' expression for logarithm of (unnormalized) density of target distribution. +#' +#' @return A list with entries +#' * `log_density`: A function to evaluate log density function for target +#' distribution given current position vector. +#' * `value_and_gradient_log_density`: A function to evaluate value and gradient +#' of log density function for target distribution given current position +#' vector, returning as a list with entries `value` and `gradient`. +#' +#' @export +#' +#' @examples +#' target_distribution <- target_distribution_from_log_density_formula( +#' ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 10), +#' ) +#' target_distribution$value_and_gradient_log_density(c(0.1, -0.3)) +target_distribution_from_log_density_formula <- function(log_density_formula) { + variables <- all.vars(log_density_formula) + deriv_log_density <- deriv(log_density_formula, variables, func = TRUE) + value_and_gradient_log_density <- function(position) { + names(position) <- variables + value <- rlang::inject(deriv_log_density(!!!position)) + gradient <- attr(value, "gradient") + attr(value, "gradient") <- NULL + list(value = value, gradient = gradient) + } + log_density <- function(position) { + value_and_gradient_log_density(position)$value + } + trace_function <- function(state) { + trace_values <- state$position() + names(trace_values) <- variables + trace_values["log_density"] <- log_density(state$position()) + trace_values + } + list( + log_density = log_density, + value_and_gradient_log_density = value_and_gradient_log_density, + trace_function = trace_function + ) +} diff --git a/man/target_distribution_from_log_density_formula.Rd b/man/target_distribution_from_log_density_formula.Rd new file mode 100644 index 0000000..4fc6f5d --- /dev/null +++ b/man/target_distribution_from_log_density_formula.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bridges.R +\name{target_distribution_from_log_density_formula} +\alias{target_distribution_from_log_density_formula} +\title{Construct target distribution from a formula specifying log density.} +\usage{ +target_distribution_from_log_density_formula(log_density_formula) +} +\arguments{ +\item{log_density_formula}{Formula for which right-hand side specifies +expression for logarithm of (unnormalized) density of target distribution.} +} +\value{ +A list with entries +\itemize{ +\item \code{log_density}: A function to evaluate log density function for target +distribution given current position vector. +\item \code{value_and_gradient_log_density}: A function to evaluate value and gradient +of log density function for target distribution given current position +vector, returning as a list with entries \code{value} and \code{gradient}. +} +} +\description{ +Construct target distribution from a formula specifying log density. +} +\examples{ +target_distribution <- target_distribution_from_log_density_formula( + ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 10), +) +target_distribution$value_and_gradient_log_density(c(0.1, -0.3)) +} From 76cdafb2349b3a4fd5550eda6036a438e7a97c1d Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 18:08:04 +0000 Subject: [PATCH 02/13] Wrap trace function into target distribution arg --- R/chains.R | 21 ++++++++++++++++++--- man/sample_chain.Rd | 29 +++++++++++++++++------------ 2 files changed, 35 insertions(+), 15 deletions(-) diff --git a/R/chains.R b/R/chains.R index bd2c127..51ba236 100644 --- a/R/chains.R +++ b/R/chains.R @@ -4,7 +4,21 @@ #' target distribution and proposal (defaulting to Barker proposal), optionally #' adapting proposal parameters in a warm-up stage. #' -#' @inheritParams sample_metropolis_hastings +#' @param target_distribution Target stationary distribution for chain. A list +#' with named entries `log_density` and `gradient_log_density` corresponding +#' to respectively functions for evaluating the logarithm of the (potentially +#' unnormalized) density of the target distribution and its gradient (only +#' required for gradient-based proposals). As an alternative to +#' `gradient_log_density` an entry `value_and_gradient_log_density` may +#' instead be provided which is a function returning both the value and +#' gradient of the logarithm of the (unnormalized) density of the target +#' distribution as a list under the names `value` and `gradient` respectively. +#' The list may also contain a named entry `trace_function`, correspond to a +#' function which given current chain state outputs list of variables to trace +#' on each main (non-adaptive) chain iteration. If a `trace_function` entry is +#' not specified, then the default behaviour is to trace the position +#' component of the chain state along with the log density of the target +#' distribution. #' @param initial_state Initial chain state. Either a vector specifying just #' the position component of the chain state or a list output by `chain_state` #' specifying the full chain state. @@ -78,7 +92,6 @@ sample_chain <- function( n_main_iteration, proposal = barker_proposal(), adapters = list(scale_adapter(), shape_adapter()), - trace_function = NULL, show_progress_bar = TRUE, trace_warm_up = FALSE) { progress_available <- requireNamespace("progress", quietly = TRUE) @@ -90,8 +103,10 @@ sample_chain <- function( } else { stop("initial_state must be a vector or list with an entry named position.") } - if (is.null(trace_function)) { + if (is.null(target_distribuion$trace_function)) { trace_function <- default_trace_function(target_distribution) + } else { + trace_function <- target_distribution$trace_function } statistic_names <- list("accept_prob") warm_up_results <- chain_loop( diff --git a/man/sample_chain.Rd b/man/sample_chain.Rd index 0793c2a..c29faa2 100644 --- a/man/sample_chain.Rd +++ b/man/sample_chain.Rd @@ -11,21 +11,26 @@ sample_chain( n_main_iteration, proposal = barker_proposal(), adapters = list(scale_adapter(), shape_adapter()), - trace_function = NULL, show_progress_bar = TRUE, trace_warm_up = FALSE ) } \arguments{ \item{target_distribution}{Target stationary distribution for chain. A list -with named entries \code{log_density} and \code{gradient_log_density} corresponding to -respectively functions for evaluating the logarithm of the (potentially -unnormalized) density of the target distribution and its gradient. -As an alternative to \code{gradient_log_density} an entry -\code{value_and_gradient_log_density} may instead be provided which is a function -returning both the value and gradient of the logarithm of the (unnormalized) -density of the target distribution as a list under the names \code{value} and -\code{gradient} respectively.} +with named entries \code{log_density} and \code{gradient_log_density} corresponding +to respectively functions for evaluating the logarithm of the (potentially +unnormalized) density of the target distribution and its gradient (only +required for gradient-based proposals). As an alternative to +\code{gradient_log_density} an entry \code{value_and_gradient_log_density} may +instead be provided which is a function returning both the value and +gradient of the logarithm of the (unnormalized) density of the target +distribution as a list under the names \code{value} and \code{gradient} respectively. +The list may also contain a named entry \code{trace_function}, correspond to a +function which given current chain state outputs list of variables to trace +on each main (non-adaptive) chain iteration. If a \code{trace_function} entry is +not specified, then the default behaviour is to trace the position +component of the chain state along with the log density of the target +distribution.} \item{initial_state}{Initial chain state. Either a vector specifying just the position component of the chain state or a list output by \code{chain_state} @@ -59,15 +64,15 @@ coerce the average acceptance rate to a target value using a dual-averaging algorithm, and adapting the shape to an estimate of the covariance of the target distribution.} -\item{trace_function}{Function which given current chain state outputs list -of variables to trace on each main (non-adaptive) chain iteration.} - \item{show_progress_bar}{Whether to show progress bars during sampling. Requires \code{progress} package to be installed to have an effect.} \item{trace_warm_up}{Whether to record chain traces and adaptation / transition statistics during (adaptive) warm-up iterations in addition to (non-adaptive) main chain iterations.} + +\item{trace_function}{Function which given current chain state outputs list +of variables to trace on each main (non-adaptive) chain iteration.} } \value{ A list with entries From 2bc057ecb241e4a270f58e5afed968f12da2c74e Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 18:12:38 +0000 Subject: [PATCH 03/13] Ensure gradient a vector and add to NAMESPACE --- NAMESPACE | 1 + R/bridges.R | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 7d62b18..f7f8639 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -13,6 +13,7 @@ export(sample_chain) export(scale_adapter) export(shape_adapter) export(stochastic_approximation_scale_adapter) +export(target_distribution_from_log_density_formula) export(target_distribution_from_stan_model) export(trace_function_from_stan_model) export(variance_shape_adapter) diff --git a/R/bridges.R b/R/bridges.R index d1d5022..62bcba7 100644 --- a/R/bridges.R +++ b/R/bridges.R @@ -112,7 +112,6 @@ example_gaussian_stan_model <- function(n_data = 50, seed = 1234L) { ) } - #' Construct target distribution from a formula specifying log density. #' #' @param log_density_formula Formula for which right-hand side specifies @@ -138,7 +137,7 @@ target_distribution_from_log_density_formula <- function(log_density_formula) { value_and_gradient_log_density <- function(position) { names(position) <- variables value <- rlang::inject(deriv_log_density(!!!position)) - gradient <- attr(value, "gradient") + gradient <- drop(attr(value, "gradient")) attr(value, "gradient") <- NULL list(value = value, gradient = gradient) } From 1ec12aacf6ca735260fc47c771b4b55df5fc6a29 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 18:18:02 +0000 Subject: [PATCH 04/13] Combine BridgeStan interface functions for constructing target distribution and trace function --- NAMESPACE | 1 - R/bridges.R | 73 +++++++++------------- man/example_gaussian_stan_model.Rd | 5 +- man/target_distribution_from_stan_model.Rd | 21 ++++++- man/trace_function_from_stan_model.Rd | 43 ------------- tests/testthat/test-bridges.R | 73 ++++++++++++---------- 6 files changed, 94 insertions(+), 122 deletions(-) delete mode 100644 man/trace_function_from_stan_model.Rd diff --git a/NAMESPACE b/NAMESPACE index f7f8639..2bb2300 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,5 +15,4 @@ export(shape_adapter) export(stochastic_approximation_scale_adapter) export(target_distribution_from_log_density_formula) export(target_distribution_from_stan_model) -export(trace_function_from_stan_model) export(variance_shape_adapter) diff --git a/R/bridges.R b/R/bridges.R index 62bcba7..32514bc 100644 --- a/R/bridges.R +++ b/R/bridges.R @@ -1,6 +1,13 @@ #' Construct target distribution from a BridgeStan `StanModel` object. #' #' @param model Stan model object to use for target (posterior) distribution. +#' @param include_log_density Whether to include an entry `log_density` +#' corresponding to current log density for target distribution in values +#' returned by trace function. +#' @param include_generated_quantities Whether to included generated quantities +#' in Stan model definition in values returned by trace function. +#' @param include_transformed_parameters Whether to include transformed +#' parameters in Stan model definition in values returned by trace function. #' #' @return A list with entries #' * `log_density`: A function to evaluate log density function for target @@ -8,6 +15,9 @@ #' * `value_and_gradient_log_density`: A function to evaluate value and gradient #' of log density function for target distribution given current position #' vector, returning as a list with entries `value` and `gradient`. +#' * `trace_function`: A function which given a `chain_state()` object returns a +#' named vector of values to trace during sampling. The constrained parameter +#' values of model will always be included. #' #' @export #' @@ -18,46 +28,13 @@ #' 876287L, state <- chain_state(stats::rnorm(model$param_unc_num())) #' ) #' state$log_density(target_distribution) -target_distribution_from_stan_model <- function(model) { - list( - log_density = model$log_density, - value_and_gradient_log_density = function(position) { - value_and_gradient <- model$log_density_gradient(position) - names(value_and_gradient) <- c("value", "gradient") - value_and_gradient - } - ) -} - -#' Construct trace function from a BridgeStan `StanModel` object. -#' -#' @param model Stan model object to use to generate (constrained) parameters to -#' trace. -#' @param include_log_density Whether to include an entry `log_density` -#' corresponding to current log density for target distribution in values -#' returned by trace function. -#' @param include_generated_quantities Whether to included generated quantities -#' in Stan model definition in values returned by trace function. -#' @param include_transformed_parameters Whether to include transformed -#' parameters in Stan model definition in values returned by trace function. -#' -#' @return A function which given `chain_state` object returns a named vector of -#' values to trace during sampling. The constrained parameter values of model -#' will always be included. -#' -#' @export -#' -#' @examplesIf requireNamespace("bridgestan", quietly = TRUE) -#' model <- example_gaussian_stan_model() -#' trace_function <- trace_function_from_stan_model(model) -#' withr::with_seed(876287L, state <- chain_state(rnorm(model$param_unc_num()))) -#' trace_function(state) -trace_function_from_stan_model <- function( +#' target_distribution$trace_function(state) +target_distribution_from_stan_model <- function( model, include_log_density = TRUE, include_generated_quantities = FALSE, include_transformed_parameters = FALSE) { - function(state) { + trace_function <- function(state) { position <- state$position() trace_values <- model$param_constrain( position, include_transformed_parameters, include_generated_quantities @@ -68,15 +45,25 @@ trace_function_from_stan_model <- function( } trace_values } + list( + log_density = model$log_density, + value_and_gradient_log_density = function(position) { + value_and_gradient <- model$log_density_gradient(position) + names(value_and_gradient) <- c("value", "gradient") + value_and_gradient + }, + trace_function = trace_function + ) } #' Construct an example BridgeStan `StanModel` object for a Gaussian model. #' #' Requires BridgeStan package to be installed. Generative model is assumed to -#' be of the form `y ~ normal(mu, sigma)` for unknown `mu` and `sigma`. +#' be of the form `y ~ normal(mu, sigma)` for unknown `mu ~ normal(0, 3)` and +#' `sigma ~ half_normal(0, 3)`. #' #' @param n_data Number of independent data points `y` to generate and condition -#' model against. +#' model against from `normal(0, 1)`. #' @param seed Integer seed for Stan model. #' #' @return BridgeStan StanModel object. @@ -88,20 +75,22 @@ trace_function_from_stan_model <- function( #' model$param_names() example_gaussian_stan_model <- function(n_data = 50, seed = 1234L) { rlang::check_installed("bridgestan", reason = "to use this function") - model_string <- "data { - int N; - vector[N] y; + model_string <- " + data { + int N; + vector[N] y; } parameters { real mu; real sigma; } model { + mu ~ normal(0, 3); + sigma ~ normal(0, 3); y ~ normal(mu, sigma); }" withr::with_seed(seed, y <- stats::rnorm(n_data)) data_string <- sprintf('{"N": %i, "y": [%s]}', n_data, toString(y)) - model_file <- tempfile("gaussian", fileext = ".stan") withr::with_tempfile("model_file", { writeLines(model_string, model_file) diff --git a/man/example_gaussian_stan_model.Rd b/man/example_gaussian_stan_model.Rd index cb38e01..8202fad 100644 --- a/man/example_gaussian_stan_model.Rd +++ b/man/example_gaussian_stan_model.Rd @@ -8,7 +8,7 @@ example_gaussian_stan_model(n_data = 50, seed = 1234L) } \arguments{ \item{n_data}{Number of independent data points \code{y} to generate and condition -model against.} +model against from \code{normal(0, 1)}.} \item{seed}{Integer seed for Stan model.} } @@ -17,7 +17,8 @@ BridgeStan StanModel object. } \description{ Requires BridgeStan package to be installed. Generative model is assumed to -be of the form \code{y ~ normal(mu, sigma)} for unknown \code{mu} and \code{sigma}. +be of the form \code{y ~ normal(mu, sigma)} for unknown \code{mu ~ normal(0, 3)} and +\code{sigma ~ half_normal(0, 3)}. } \examples{ \dontshow{if (requireNamespace("bridgestan", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} diff --git a/man/target_distribution_from_stan_model.Rd b/man/target_distribution_from_stan_model.Rd index aa222ef..22ae35c 100644 --- a/man/target_distribution_from_stan_model.Rd +++ b/man/target_distribution_from_stan_model.Rd @@ -4,10 +4,25 @@ \alias{target_distribution_from_stan_model} \title{Construct target distribution from a BridgeStan \code{StanModel} object.} \usage{ -target_distribution_from_stan_model(model) +target_distribution_from_stan_model( + model, + include_log_density = TRUE, + include_generated_quantities = FALSE, + include_transformed_parameters = FALSE +) } \arguments{ \item{model}{Stan model object to use for target (posterior) distribution.} + +\item{include_log_density}{Whether to include an entry \code{log_density} +corresponding to current log density for target distribution in values +returned by trace function.} + +\item{include_generated_quantities}{Whether to included generated quantities +in Stan model definition in values returned by trace function.} + +\item{include_transformed_parameters}{Whether to include transformed +parameters in Stan model definition in values returned by trace function.} } \value{ A list with entries @@ -17,6 +32,9 @@ distribution given current position vector. \item \code{value_and_gradient_log_density}: A function to evaluate value and gradient of log density function for target distribution given current position vector, returning as a list with entries \code{value} and \code{gradient}. +\item \code{trace_function}: A function which given a \code{chain_state()} object returns a +named vector of values to trace during sampling. The constrained parameter +values of model will always be included. } } \description{ @@ -30,5 +48,6 @@ withr::with_seed( 876287L, state <- chain_state(stats::rnorm(model$param_unc_num())) ) state$log_density(target_distribution) +target_distribution$trace_function(state) \dontshow{\}) # examplesIf} } diff --git a/man/trace_function_from_stan_model.Rd b/man/trace_function_from_stan_model.Rd deleted file mode 100644 index c148281..0000000 --- a/man/trace_function_from_stan_model.Rd +++ /dev/null @@ -1,43 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/bridges.R -\name{trace_function_from_stan_model} -\alias{trace_function_from_stan_model} -\title{Construct trace function from a BridgeStan \code{StanModel} object.} -\usage{ -trace_function_from_stan_model( - model, - include_log_density = TRUE, - include_generated_quantities = FALSE, - include_transformed_parameters = FALSE -) -} -\arguments{ -\item{model}{Stan model object to use to generate (constrained) parameters to -trace.} - -\item{include_log_density}{Whether to include an entry \code{log_density} -corresponding to current log density for target distribution in values -returned by trace function.} - -\item{include_generated_quantities}{Whether to included generated quantities -in Stan model definition in values returned by trace function.} - -\item{include_transformed_parameters}{Whether to include transformed -parameters in Stan model definition in values returned by trace function.} -} -\value{ -A function which given \code{chain_state} object returns a named vector of -values to trace during sampling. The constrained parameter values of model -will always be included. -} -\description{ -Construct trace function from a BridgeStan \code{StanModel} object. -} -\examples{ -\dontshow{if (requireNamespace("bridgestan", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} -model <- example_gaussian_stan_model() -trace_function <- trace_function_from_stan_model(model) -withr::with_seed(876287L, state <- chain_state(rnorm(model$param_unc_num()))) -trace_function(state) -\dontshow{\}) # examplesIf} -} diff --git a/tests/testthat/test-bridges.R b/tests/testthat/test-bridges.R index e7fc7f0..f8ac4ea 100644 --- a/tests/testthat/test-bridges.R +++ b/tests/testthat/test-bridges.R @@ -25,10 +25,11 @@ test_that("Creating target distribution from Stan model works", { expect_type(target_distribution, "list") expect_named( target_distribution, - c("log_density", "value_and_gradient_log_density") + c("log_density", "value_and_gradient_log_density", "trace_function") ) expect_type(target_distribution$log_density, "closure") expect_type(target_distribution$value_and_gradient_log_density, "closure") + expect_type(target_distribution$trace_function, "closure") position <- rep(0, model$param_unc_num()) log_density <- target_distribution$log_density(position) value_and_gradient_log_density <- ( @@ -46,37 +47,43 @@ test_that("Creating target distribution from Stan model works", { }) for (include_log_density in c(TRUE, FALSE)) { - test_that("Creating trace function from Stan model works", { - model <- cached_example_gaussian_stan_model() - trace_function <- trace_function_from_stan_model( - model, - include_log_density = include_log_density - ) - expect_type(trace_function, "closure") - position <- rep(0, model$param_unc_num()) - state <- chain_state(position) - trace_values <- trace_function(state) - if (include_log_density) { - expected_length <- model$param_num() + 1L - } else { - expected_length <- model$param_num() - } - expect_identical(length(trace_values), expected_length) - expected_parameter_values <- model$param_constrain(position) - expected_parameter_names <- model$param_names() - if (include_log_density) { - expected_names <- c(expected_parameter_names, "log_density") - } else { - expected_names <- expected_parameter_names - } - expect_named(trace_values, expected_names) - for (i in seq_along(expected_parameter_values)) { - name <- expected_parameter_names[[i]] - expect_identical(trace_values[[name]], expected_parameter_values[[i]]) - } - if (include_log_density) { - expected_log_density <- model$log_density(position) - expect_identical(trace_values[["log_density"]], expected_log_density) + test_that( + sprintf( + "Creating trace function from Stan model works (include_log_density = %i)", + include_log_density + ), + { + model <- cached_example_gaussian_stan_model() + trace_function <- target_distribution_from_stan_model( + model, + include_log_density = include_log_density + )$trace_function + expect_type(trace_function, "closure") + position <- rep(0, model$param_unc_num()) + state <- chain_state(position) + trace_values <- trace_function(state) + if (include_log_density) { + expected_length <- model$param_num() + 1L + } else { + expected_length <- model$param_num() + } + expect_identical(length(trace_values), expected_length) + expected_parameter_values <- model$param_constrain(position) + expected_parameter_names <- model$param_names() + if (include_log_density) { + expected_names <- c(expected_parameter_names, "log_density") + } else { + expected_names <- expected_parameter_names + } + expect_named(trace_values, expected_names) + for (i in seq_along(expected_parameter_values)) { + name <- expected_parameter_names[[i]] + expect_identical(trace_values[[name]], expected_parameter_values[[i]]) + } + if (include_log_density) { + expected_log_density <- model$log_density(position) + expect_identical(trace_values[["log_density"]], expected_log_density) + } } - }) + ) } From 3fe2a84d0c285710cfcdc8773e953dc42d0e4762 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 18:37:41 +0000 Subject: [PATCH 05/13] Allow passing formula or Stan model directly to sample_chain --- R/chains.R | 56 +++++++++++++++++------- README.Rmd | 26 +++-------- README.md | 35 +++++---------- man/figures/README-banana-samples-1.png | Bin 32852 -> 37401 bytes man/sample_chain.Rd | 38 ++++++++++------ 5 files changed, 83 insertions(+), 72 deletions(-) diff --git a/R/chains.R b/R/chains.R index 51ba236..8131791 100644 --- a/R/chains.R +++ b/R/chains.R @@ -4,21 +4,31 @@ #' target distribution and proposal (defaulting to Barker proposal), optionally #' adapting proposal parameters in a warm-up stage. #' -#' @param target_distribution Target stationary distribution for chain. A list -#' with named entries `log_density` and `gradient_log_density` corresponding -#' to respectively functions for evaluating the logarithm of the (potentially -#' unnormalized) density of the target distribution and its gradient (only -#' required for gradient-based proposals). As an alternative to -#' `gradient_log_density` an entry `value_and_gradient_log_density` may -#' instead be provided which is a function returning both the value and -#' gradient of the logarithm of the (unnormalized) density of the target -#' distribution as a list under the names `value` and `gradient` respectively. -#' The list may also contain a named entry `trace_function`, correspond to a -#' function which given current chain state outputs list of variables to trace -#' on each main (non-adaptive) chain iteration. If a `trace_function` entry is -#' not specified, then the default behaviour is to trace the position -#' component of the chain state along with the log density of the target -#' distribution. +#' @param target_distribution Target stationary distribution for chain. One of: +#' * A one-sided formula specifying expression for log density of target +#' distribution which will be passed to +#' [target_distribution_from_log_density_formula()] to construct functions +#' to evaluate log density and its gradient using [deriv()]. +#' * A `bridgestan::StanModel` instance (requires `bridgestan` to be +#' installed) specifying target model and data. Will be passed to +#' [target_distribution_from_stan_model()] using default values for optional +#' arguments - to override call [target_distribution_from_stan_model()] +#' directly and pass the returned list as the `target_distribution` argument +#' here. +#' * A list with named entries `log_density` and `gradient_log_density` +#' corresponding to respectively functions for evaluating the logarithm of +#' the (potentially unnormalized) density of the target distribution and its +#' gradient (only required for gradient-based proposals). As an alternative +#' to `gradient_log_density` an entry `value_and_gradient_log_density` may +#' instead be provided which is a function returning both the value and +#' gradient of the logarithm of the (unnormalized) density of the target +#' distribution as a list under the names `value` and `gradient` +#' respectively. The list may also contain a named entry `trace_function`, +#' correspond to a function which given current chain state outputs list of +#' variables to trace on each main (non-adaptive) chain iteration. If a +#' `trace_function` entry is not specified, then the default behaviour is to +#' trace the position component of the chain state along with the log +#' density of the target distribution. #' @param initial_state Initial chain state. Either a vector specifying just #' the position component of the chain state or a list output by `chain_state` #' specifying the full chain state. @@ -103,7 +113,21 @@ sample_chain <- function( } else { stop("initial_state must be a vector or list with an entry named position.") } - if (is.null(target_distribuion$trace_function)) { + if (is(target_distribution, "formula")) { + target_distribution <- target_distribution_from_log_density_formula( + target_distribution + ) + } else if (is(target_distribution, "StanModel")) { + target_distribution <- target_distribution_from_stan_model( + target_distribution + ) + } else if ( + !is.list(target_distribution) || + !("log_density" %in% names(target_distribution)) + ) { + stop("target_distribution invalid - see documentation for allowable types.") + } + if (is.null(target_distribution$trace_function)) { trace_function <- default_trace_function(target_distribution) } else { trace_function <- target_distribution$trace_function diff --git a/README.Rmd b/README.Rmd index 29e851b..e7763cc 100644 --- a/README.Rmd +++ b/README.Rmd @@ -72,36 +72,22 @@ cat( ``` As a second example, the snippet below demonstrates sampling from a two-dimensional banana shaped distribution based on the [Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function) and plotting the generated chain samples. -Here we use the default values of the `proposal` and `adapters` arguments to `sample_chain`, +Here we use the default values of the `proposal` and `adapters` arguments to `sample_chain()`, corresponding respectively to the Barker proposal, and adapters for tuning the proposal scale to coerce the average acceptance rate using a dual-averaging algorithm, and for tuning the proposal shape based on an estimate of the target distribution covariance matrix. +The `target_distribution` argument to `sample_chain()` is passed a formula specifying the log density of the target distribution, which is passed to `target_distribution_from_log_density_formula()` to construct necessary functions, +using `deriv()` to symbolically compute derivatives. ```{r banana-samples, fig.width=6, fig.height=4} library(rmcmc) set.seed(651239L) -target_distribution <- list( - log_density = function(x) -sum(x^2) / 8 - (x[1]^2 - x[2])^2 - (x[1] - 1)^2 / 10, - gradient_log_density = function(x) { - c( - -x[1] / 4 + 4 * x[1] * (x[2] - x[1]^2) - 0.2 * x[1] + 0.2, - -x[2] / 4 + 2 * x[1]^2 - 2 * x[2] - ) - } -) results <- sample_chain( - target_distribution = target_distribution, + target_distribution = ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 100), initial_state = rnorm(2), n_warm_up_iteration = 10000, - n_main_iteration = 10000, -) -plot( - results$traces[, "position1"], - results$traces[, "position2"], - xlab = expression(x[1]), - ylab = expression(x[2]), - col = "#1f77b4", - pch = 20 + n_main_iteration = 10000 ) +plot(results$traces[, "x"], results$traces[, "y"], col = "#1f77b4", pch = 20) ``` diff --git a/README.md b/README.md index b7aa09b..fe9cb7c 100644 --- a/README.md +++ b/README.md @@ -76,39 +76,28 @@ As a second example, the snippet below demonstrates sampling from a two-dimensional banana shaped distribution based on the [Rosenbrock function](https://en.wikipedia.org/wiki/Rosenbrock_function) and plotting the generated chain samples. Here we use the default values of -the `proposal` and `adapters` arguments to `sample_chain`, corresponding -respectively to the Barker proposal, and adapters for tuning the -proposal scale to coerce the average acceptance rate using a +the `proposal` and `adapters` arguments to `sample_chain()`, +corresponding respectively to the Barker proposal, and adapters for +tuning the proposal scale to coerce the average acceptance rate using a dual-averaging algorithm, and for tuning the proposal shape based on an -estimate of the target distribution covariance matrix. +estimate of the target distribution covariance matrix. The +`target_distribution` argument to `sample_chain()` is passed a formula +specifying the log density of the target distribution, which is passed +to `target_distribution_from_log_density_formula()` to construct +necessary functions, using `deriv()` to symbolically compute +derivatives. ``` r library(rmcmc) set.seed(651239L) -target_distribution <- list( - log_density = function(x) -sum(x^2) / 8 - (x[1]^2 - x[2])^2 - (x[1] - 1)^2 / 10, - gradient_log_density = function(x) { - c( - -x[1] / 4 + 4 * x[1] * (x[2] - x[1]^2) - 0.2 * x[1] + 0.2, - -x[2] / 4 + 2 * x[1]^2 - 2 * x[2] - ) - } -) results <- sample_chain( - target_distribution = target_distribution, + target_distribution = ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 100), initial_state = rnorm(2), n_warm_up_iteration = 10000, - n_main_iteration = 10000, -) -plot( - results$traces[, "position1"], - results$traces[, "position2"], - xlab = expression(x[1]), - ylab = expression(x[2]), - col = "#1f77b4", - pch = 20 + n_main_iteration = 10000 ) +plot(results$traces[, "x"], results$traces[, "y"], col = "#1f77b4", pch = 20) ``` diff --git a/man/figures/README-banana-samples-1.png b/man/figures/README-banana-samples-1.png index c46b203a09ee0c6d9f8791f88671e2046ecf7fd0..c6304115cec450cf562befd618a3e96ba28c2014 100644 GIT binary patch literal 37401 zcmeEu1ydYd)Mf)gf&_;EL4v!xyGw9~;10oE65NA31b1g}OR&J;?j*Pq++q9Ot*>f- z!`4zoQAJO~boagIJa#7Xld=pNG9fYm0BCZul4<||D+&NGK}fLRZ$|FqWx*eBzsTyk z0>Hbz|Gr^{oXX4rfDDk66w~m^K3=mDRky+i{6nVd0=sKJQ6xkb2jbb!QZtyDm4 ztC#hvD{bd&EA#YDzBHZc-mb&pN)nJqG~A+cylFUn4tOI$4O@UAG*P?l;(A#?xz~O4 z;adT%DALD}E41o1q96=$BrXI5)({MFpSN)D#gRUCg}|VG4Eb*K3kbsymt`>lB#=H9 zrLhKm48fRIc;5@Y{(uVg^MB`K{r}VcpSvY&nH?}r6HoV-P8;28tJR`N`JTvWP@#aQ z*4WswmB6E0+Mu3My2nkta8o3DZidU|?d(5_i+b1k)+$OTxfC&~D8x@~4kKE2+RhpUoN(XE46Dq+}*oec1CJ@UfMm+ z)^{DZvanJ#vAecky0fCdB6U?#d|5dgha{<1ZPjefL>>gb&Cswg$K@uQwRZP`fq{!) z^-^W{pqrDGvyGnD!gyiZd4z5utVG*WYO;`Z|Gv(XJuKDa>070sM%_!U(gg4h!vUV2 zw!$==zDfivdSLZAy zp4|o#WDcmX2!@MVB8m^LXVqZSKO_yixx4%C1*bL`=D`QF3m8bN0lj8y$`gm{t4!iZ zUHGG@iz6zR(8?r~H%1;zj`^a*dRl-l`t_GaStAZg{>uK*?+2zOqzr-Jf(IBVQkO9m zD4ORv^cnpRnYtBB8>7ByZR}q{(eg2X-iayHLVWRk2g0eJKL>t^p1T?RNTO1ifPH_{ zoB1={0iB$~ZoZ;P-#frD;&JtbWNd{pWb$z^qK}uOa_DC`L0hpV4D)WEz#^AtG$AUS z4#y`+mbS44EH<$4uM1-{UskR=WOm@jhD$_@#s5mV1h4(e``%0#nEpx?uxn`0pne|f z4XVcKT|K!XJK_*Wf^dR|#a2~<1$c?zFhns(?*xpVOJrPI{%alR%$wFErB=0fdudZt zzXmB1#GYW_SIDXh-F%}1SvH(6Q{=rxaPw-~uq^OxjItmVo`^N$I*>11OhZod#k{-( zI#s{rI5M%92NRU@_O>OS-_6c!$~;XPi3>9!IlKPsd~6RMm>ik@vbU}5IIojl&|Qak zqle(sMw65P0HTqU=6iEIgjW|x22mF?-@GTTTWC-jfA z$1m+UjPXXmcP)*!J9GQXW|^MPz)U0yMmRNYF`nkUbDUmUkND~zvV>+2YY`FsZwN2B zolWi%M62?|-ZWGxqI5QxlzU2(@82a_FX);)^6UArt)>1kF%eqGKxo3NK5db zLr*jEm0fAZ=Xd9ud#}QdB_#(iAH-)*H8gHkDW)!lwt}-4;?Ac6a+klm_gQ8AQo9Qc zol~|%STDh&eA`6bURbwXzMuzDqt{&U?WWpXPs7-{;}`4nVawK_SCW3 z6TrU#FosGH-!wOzKOeYwr{Dnn><1%_=W|#tmU-O3r+l_o4ISN+ezk+a*vCgPgq{)# zO(v?LGo--o3Zvk+g9f{(#d?(;66)?be7{$kh4Ceky1Jw5y+>5AFGiTWk6C=0wS;9c z-W2nbp05G7EmOSyg%-OM%*@Ja=uSP=PO%vODo?>~_WR9i)9;{QY$Xt-fu z)*Pk<0OX*8!I%3oJ^B=|c~mc8@1Xm5>(G<#u^%W~9tMAB)1`OqG<{2FeW!*AU?|8{ z4v@gCe!KXk14>`3>_XPP1)&=<;~hElZw>*4ZcI{5`&EQqd=w}!4@K1K2PQC zgZ{xl3Lpx{%4LqL9fd2xzRNmNNg2-}kD7iJq_vAm2m^2uewXLTR)RDitre2=f~`w? zLVNyIQE2DP{_xCR@0QQE)Z^t(yss228^D`Ebo|0t))&QK+oItX_AqaWf4K}ZY{t>9 z89n@K#CEJZZMDCKIdD}`OAx~Lmyt~MnfImc?=N@9ANI>=&T05kkmBT#tWZ!FBHFXo&*v-`; zHMPa>5at*}J0M=z{lXg%5RgQpLL%5|b2Han2`)JjrnhvAwIv+Zy9kT`nRSQ@KP_y; zxaiQi0c>as?li$)5FgT(@2?n?*`0b>9&swu!0cd#3 z*;cs)JOS(&U}pooD zoxg--Au`jSK|o)6=3P7*`e~%JU=%W`PfkJiRo1^+OJ9H_=bZt6PQf1=kBoVIY0Fe* zqV@HmpBk0MpVbtZkhAC$x25!*XBJL{PhRTC9l z^R7#Sposr{ZvBOU0VKqH%4Q?1HE?c7Dd6d$@rizp(Q(iqEZRb;S>hXjPE?35uj_&n zCo;nSLW=@~JuW*v_xE5b(PKqCzWAtj&)%0dI@^i-du~lqr>L1pELV9K78Vvq$YneE zL(I|95yZp}soWr~BN=AIWEOgE(i4Xb(vfy-eIlo>EHN&AkEO9FoL_h`Q{1;a=YQU) z;!@MUP}PYaVFAaz?$~(uZ6zMxiX@`yG1QpUt7@@KLalDUhla0M^ZGILD$EBU;)BpPd zdLOQjhT=)jwy$d3`Lc38F|K)CO9f@t2 z_|;hW9~X~CUUUQ3Zp_K{8z8Bp8N7dv)Hng2g3|YI_R>(oE}>E*L$!nDPE)Fu4*hgE zW!ZT;KP;-$Y7a(3TP?>jNdxZk!OO#@R=oSyhdwqowqe=0t|i z2I-O~x9VR8d|0lDFV%EydG~SlJ*kJAf$K|#9}*_;-fcNKw;D&pkh^NlYE91q2c^qC zmNCV}*3M$sEO2tl`MLOw)G?d+OwZv?$@O---CWJ%q}AQL?+C%rMwh%9KhN)-y}dj^ z@9FjRfb~xA=H_M!=@?#@?TAkia+gNX&7^wrVuCfS{5Jp=%A36B%TKOVA3hf0qw2&& zyZU;3;bh}5hqQ3;aW!V)s#w(5H@1nP^FXzwS@;OiPG2US41T!!3M63?ACod7#|<;L z;bb8K(jQWo!}HV)ZB+d~DBk#!^qWirL9lfkf0@1ZTs)4Bzt%4uCs4@z*EDHyas@n3 z|NZUD3wWSQ8vY>Ub48QqhXwAw&!4Fu5ft{0s+c`k%7638S}}hP`J}LG`_ryM4GZWs znVk3k0|)R~Xnw7%jVhu9CXq}}5ZCR-MOZ4auz?#(11pr^vv?7Lle0-pnXc>H(^Qk_ zZxTsS`m&TQE3<>iI+<{L->(7-M}JZv+OE98{#f0o{S z!=L0t(CHeekCeZ_>GSA61!XLulwuMfCT$d?bXc3L&bzZp(;)Uu_P$cxy(2bB+q!Tw z{5c!MH@e_h*K*Q<-}U1ZxY?g7q*6$2a@ts*nQ@&dQCjQp?D4s_b#``^mmdMIzf?IV zGtYPdhkUxpF~;enS+e*l>Ye3kH!ot?cM>d(VUG=!uImB}fcg~qR?lO~cVOq)gy7oS zpuX{#%W|=jwZiD7DG^JHlj8;~A zY-WQ`*9*oV>nZFSrCqAftc)e%IbE!$`1tYX;Nb2U*Xh;440yD$!N;?Ei4S$Tkm_S; zVZ%%m1K9SAeKPU040mm9cpFuS%xZi_`UkxejL&a#lF7f8HC7HJpH}mZN?RyPBF*)s z5=9c#c77M3wp~#8uAn@SswHoE?NoS~0zk(F*E~fGUlIdJ%t>;ln7}4`p^1<~O{IAw zw)Za-_(Z3*G^I(8)8{K^1!QsBOoOzV-~F)Yy3~M~hmLM|U|=Vf3l;c0Tc)~DV|WDe zH6s7C7N??~;HR(P!W;f30=H}yfDuAyGXvd+|6zn@aa7&k4!5}ecpyJf<{Y@DqN|_P zzb7j>n0h2Ac@q|!rN8nG)QaE$qusnMmI73FGk}K%E88S~=ZP;EB|fS`ncn@a>QYmz|HrU?5hd+x{n zNNn+ZRx+}nVU7FdvEkwJs@XP|vz%mY#?$tfOg*L4)dpL`z1!Cl50`vUd!*gjr=r18 zNP$b2gx!(1_(2zhTbgGINR)|j8+r#EJEuNTR8;y1W&M2L_1yh*UYnK|{ydr2r}ow~ zT>5qb&K{r`M!~ru1+Q4VZ?RC#oRIokXTKerv z7;{@{Vl23sFMkFjDx@)^At5b-jXQ>bgU}fh5iz{pZ0PXeV5!ln&7h_xula+f8k592 z_u%z1ES;6ruVJRj(^Jd)OKT?xDLjN-o2ZS|lxJ)cGeob7FT7w}sQh5&Ltr-K_L7Y7 zT#rrhoK-zya?IpNcNOuEPk)x( zh~DI5^}R=h1;Tvzye1&;8=g~TrEe}9UVIa8;E~U8uK_$lLVm|pUVxR8)1brSWMy;@ zYRJQWabw*vx|7}`k#mK! ztp8NtHS&2f#5AMrV@G_rCIgDZ*d=0fOVe|3^-$DH zS~sOAtwU3`gQs;8C!dFjxzqp?q_v#o`vL4d;aNLJ2A9FW5OGV=PQC)4d^P_gqmMTd zfHa^6K|081ik%Trhl3csK8hLj7*rY3vLvWAqDzDLW?~=}yTxVlsOgEg@FkRaLVCio z5d<1wJFcKMmd%7dZo`}z*(7oEAKl*g(t*yEiRh%qCJ-6k=>E621P;Ov|9`zU@ z`4UyFzUDTF7&D$m7*sCh=BA#f>79#ZZGT(7uX|65z)KRCO(o*kZ({u}eD*JhHyKg+ zHb}Km->X^WnFoH@p6x>u@uzpQEOe4CDJcm@Clv-2OfRoH-;{T+&U%r{AmyjV?>T@0 z=G54ujdv{7X9^$w_Ry3Tx^Bm_4DA@E&ot(AIZM~qEj(n6u!vyf<`eV+z2p6VTH>F9 zpn^a0e(P1i30MU10M4iQ83oH@e~q#xf0w;A90DmAQJ4HbV(hOMEnR{RF`0+j()&s? zjMkfbRsN6iDUzwGj7=Q7vsUW6dDw9Bp;$dB+R5n4?_E)sxQr`$;!arG3;WvkKUksk zvG}&J6E2dr{^OK6`JMY{^S*bsY4sz-d24VOkz<3Dq*>Y%D>^yeo_G%(03UQt^jI-% zlo&d(<`X`?=ZS_nt4nrVj~9VSC+k?XI9T~kB}xiEOjR5-U~sY*=*zH9vG z{xCJ44A)MDkIY$!52O>+y&*5(b>_?0xw$>rp6o9`am{3OrxSOY9r~&);diVS^@e#j z%yxx*XQYC1=2l2c@Oei~9Q-F2f^ekI*MCfG*#J?$mzRGr5mG-#D*wtn@`W@(?srUo zg_e5zecf+im|#NHbvSV=5u6KU!EjCIlsZEL@X#*!2@|BctCHksB0}L-xUU+{I&b}P zEu(*#PLij-3B*1Yo?Qr-843vRK86KA)SWVskdQ1})m;`KDnJDM`_xriG!>z4d160B9fM-0itRF7D^z&42q&3eqflG9gmq{+Q3g8$BzKW=nLI8Yr;Vdf$NMq z>4KqhKu=oIMZQwn5d}@hAM$(w-FkVKBx@(UMr|;^aL^cOATmQ%-6pzd>!m871tmK-ji zbFqlneQZBq=FpbgbZ$E~xtZ6JRKv1=Vy>mW2Fb}`D`SRrbQpfLF@FS+W|tr%KI&Cf zS?c0@*E`tr2k2Yy;a8j}m8Z({F$HMS_=r$#YmM@KNKx0b7&ogX+S0!rDU1{#)zYX$ z=UdGD0+6w*`)`<=1~`~wMOnDMlQQ%zNu0@M!KfC*M+aUX7oC!O;DGFyk(U>*7{AQ6 zmhidUa>Mj;};bI8{9m7VgHT7O?f#xl)h(Ugp4nC(Qy+fH-ZEiOjHmji`?mj>tI0U z-VBlBbk4hrl=8Bey&)@$Dt-(Is^+nEg+y9xmpR3n%_l!qjy5%J8IPV+8$Y5b?-&UV zAcTQRH$7HWRxS^R?#yLQK^LYdi>lFptD8>~Y>RE~N{oVcWLrZV^_0aH0*y~4U#{=1 zf;%0Yx@fZlk0VVf|J_9Y(h|wen<~9V4_<~9wNtT*z{+pWRC|oHSpPL zr0vO}mQ-bMH(Ar_B9y+(3ZO>S)jDe1L`ucAOIFTE71iV{uCrczLq_ka#>V#B`|D=%80Dii^9)u@ zTfNc(mx|)UV0TJ*hVg#%3%5EhD>_o1&AQofunUu9O(}Kj^bMw8%6*6a4Cm7I6ta)8 zymp5RN5eoo6d~Z+RcV&(wVxQ9&y5eML(7tb=bCa>od=e#MSpREdc;>g@DH){an4B3u zku3Dr0nhM@vCCian2`chIJ{uFgF>Q_GGMILHtr@pLPomoZ+T93Ns|C zCQHhX3F35VKzNPvn)xrkmb;JRoOXWF9jVQ>JI{I{u^aMQ=xW8xZ#fXRKRuRE-5?Dl zH}-%53i%fgMsPH19K25sGK9kJ7In`T$MFz);?H>N3{PZseZ<$Fi^!#>*M|>5jN=*@{LG5SY@k!Cz4;*k`@TW@2&KMFHKW*e2=eZw7vN(hgR2n zJ}FG#Ps^dLjfJ6<9vy+ju)6NXV&GPF@;7~!YQui~;hh??uw{KS*sXQIbV1PdUMYN;Oh`@7TM~Z=l+^08E zqIIR3w_3|6BA~|+9vC8LP6^0XP(mmA+m`;dx73d^VOQknsxq;gvoS~>=AXOQ5i}{0 zKEk&^RM)9u4<&?<!U@QT`e6|-kv7vdA;40#fJX% ziIh;_iy#*l7cK4Y)zzn1?)|O-i>Y|!GoL;0pBEhGjMOI+kvSkmYRW0Ys-q<{sBSBV zhulW43_gZZ;dk49BdChy1wK(F#6DjJ6bzOqynQTcU=gq zcq1xije48&|4A;kHgl3{OFW2$SX$@^LEn(4%FB-FBdhY#Ao0@Dk`^sk$aCmZ0I-Yp z<@Jsw#p2XPYWe7h-G-A+%68mWwSoZ`T3O66F4V}XmM6uGX*!Y5=60h4n)HaM@BgiN zZou?LV$&?QIEi@u`>die0cw3Mi=DU1EGkD?8`o`uEopK4@PerqLWpct!F^-i{w@By zRwln)Alwl~Q-QLG;mi(Iv#GbX%4XNsP0Jn}2R%c&nNui0Q%xN`ZWLD`KKzfiScERP z>@m$3Sbbevy0MvKtr+WC-P2M%d1`!Q(1Mp}Sv7S1yHrxi66S=a8SQGzNingvW&ra8 zltQ%q*A6*QR)=QzeRO3>(HZz--@OqNRmwN7X8BKm6nM-Gq7U5aGOZbMltBk{Ahw!czI>6FkLu}pO{7v-H@*DcUdn(7>Y+v0y_6B|W?g}#8!-P- zJV7Hu7gr-=b>*>~obtodRH@*V3jgo#ghv;@+}|gvt3bL@v{zz;!i$T0l!3(z&d;i* zko3;yPTEFEe_I#`XX9g=E`!C_`@h%R%7PQMcaO>=pvSwTsWdb9UA8xHZYQp?#wd1T zp{nlqv+2Lv9Gsm3{37?Ry)w1m4z2=X&N4A8KX?MAh`GyooL$OrYX^nRkQMV;@IqUI0A>9q`Yc zRfwa*SI9VGUD(Psm|^3l`$HN{YNEejb@nI({q5|~QhN4Mc)@mji;;@pz5DILM%mHp zrNvKZO^;q~CVSg6~BSuGP@gDdy*^jF73fhuY%LIbCr=_d_3OTi0y2I&W zg%9s!)^R!5*()l1*7UxVO;r?cH>lS;FX1{=p=rn8uVcz&o zNSOum#--)SCgkRqzgnV!$6>i?)vlb!0ZSSd-n(CDQ%q^)Fqd1s)ah#ml^O=TEmaR^ z0(8(Y8#?s`vP^hveh51QuKa`d5+ljxHDk-FynF8O#8Qkj1`a&Kgy16&BQ9;q@1~Q( z=zpXS&jB1u>Cnu%7PQ+P<*LK>(<=b5M%B|7z`l>NuTOs3pFm{O<$^|zp=n5qfjxRQ z0Yh0?cpC|TDi5srB`horY%at<8=c>r3FT}AiHV8nH`~`YHLW(;%=CVR1J z_;Roly1_dZOqoB#var*UyX69{q1hIBaSn2^K(we1Z2h6d;Bg+g*qy6oRWNG>d%;j) zK`8~ZovdY%rlF}Y9II_qgr{|N-Ed#j*tk55z&KbloMU+$IEo(u(>op_AkyDeCBN!C z$K~Ba52?Elw>RpK|Cy;2rO2R~WrL)&&C-9df5u-PqgEN2t-JaL?)RJG=DSp0;X)!H zqIkRC5M2Tr1tv%uF7XXew^^TOnbvINmMc~o$>2*DjN(gxqG`e}5q~ekkNMbOg z@O>O;4tRd}3Wtb*h)6_41S)BW0Kdm^dTQ!_dSutYZoZzCSwo7ehenin;D%Dh%!bnC zXptaH`3c1>XNSQh-X{{TZ8!#$SKo$1`%Y0t_@K%OWMwHU(wue8d-KEs(BK;%KlK(V zj~M%Lec8Zm&LO<={Herk?H97DHI{Lso4uXgwV@<8^`XuFcalh1hHi{3Xwl4|o7+{x zT$@yk7)>X!Vy@)^v+41lvMw^0;JS3qd;G2s;fAg>w+m}Hpu1QWX#{+tEz3ZiA?yNz zYRsBga#-vgEYb~V$&6}fstVDIKXh#mdbl-lT58LcBeGneFlbW0>Y%pk=y60@oBY-j znJ?{0J2||ZNDidShte^ zHcT%ukOtVjZ8LDXA0F1OOqm#bc#H$<@AV3HosoYSFmju-(GRFi>a10geX1b1VTJIv zeVQDMq6zdl(S6y(u7}4ZClq6~327;mOXJ=RvzKm{?mMN%vqoT5Ac1ZSYt$QOiIp-D z>&3@9mPx7Q&OdL&~ zGcq!glD1AxgglNHTjm+x0HCRnng`|}YCY;%bvTx-zW542esk)W?QPXAU#dB4unj+|)H+|J$Qb4$6 z_50XA`vNzTRO-?GIr=`kZiX7`-M9+$Z)Agr9m{BQwkG6n zz4sowFNSu*6s6Xu(yONF>#JI_AbM1z%{S?~70m9{No{`j4LLbEjg5`*@#w%CL`0{x z_G<^OWB)W=#_aAhXCrSepy)ExHC}Qhy{`1;i>AFC!Ok8=l?I=bW7^?DRHiN6uvxv@ z!aE_Jj+JFLVT9$>fy#putt6j4IW65OsPp@QQl6ep{tn637YSomi4@qo<1YMkF=;~&;-ICHT|2C2e z(K#{;Z1{uM=!)OkApZelIMbL88S;f0oUrf8w0b0M-2y_vcw^v5G974Xoh&u(70D$5 zKfn$Pc2l&hv0RYy>D&MHJv7Obr&|Q~tefRQgUQJ*ysu;t{-Z!fX1BqXq$;ba{(j@C z3%^*>??ZWxiuAGz57tW$R@Y|M^yWS1Vx7IsHI)uX1N$x@YWnr_82q*Tp^K&|G7zVH zhSz!Cd~xn?a**i>flat_G{pq2zBF=HxW1*b{oCHu%!{r5_`Fg)o;c(rMsIu+@TBMV z4S+SL1o3(y(vIBbc}|XcH2tgRiIj_pg}%U+_s3H(AHP(oCpG;gF1b`NKUs!=VT1KnXNo-IMiIZAk(9Z) zIq*IV05^E=UDIu46k8R6KhE#{g*5q%UMg#nk*NT((Ds6~r{Jc8&4GmtMYb#rTBqocmMAO9=D(lM;HJ>de)2z)U@CsX{5|O z4GUZqsyH!MjmCmS4iYoVShYHzVd8cXxRM|^v;}DpOuLt$W>D=waO!3-roZRxI5wd- zuD5R(Z)5%@7SuGZG=4+9IzJzl9|^gR1T~qy=GCe8?0|hJ!B+=t*uz$8wGxh?;^N}a zkTR>jqm`C3FgUXDB?$u%MfX(jw96%*P!EdEsBMPm`+aR~3ecV1&g%oieYfi>9rZur z--bhWO`n2`OarZEeC7P=HM)2gPbZs#|LT%7ccxG*S4y#bBcuIs9NQMC@{Ww|FwfWK zk)){7-(Oaip1P5dy7BlaX*InV&fDOjUE%EoC=uN6%=H+^FS4~^14tFcdB>>HP^MBC zHyelZiDXgf905^Ph0j0Dh;!$QTh&H?U=A8wXH|ziJX=nWchFr7y{@Zlu^*ruCERW4 zI^Kh_916XpSoNsPqreAW1FJka{UW+DjQFd@j)YG?IwwU8ZN#;;ugC!epA{|SBF`Dl zMA~l7&fw?y^B17*(s6`EJ2w~=L{U+{X|?>2X#qlv>S|Ofc9*i+oXv5pxx8^;7tAe_ zkO+I#jc^_m%ET`>{nWc~f=iPV;V(Qoot*85{ThyTwmBm>L52S{hzvhkV#b!;IT^W# zsHV$$KSd!Vu`RstVnhlwqyREz{B+LD+D2VP>BtKk*aQu-L_8|Z@Ul2GSU^g3d>4Vn z%A@5^y+?RS4L`o|(8{o^?cj^dW6yE6#t#D|1e|*BU*mN&x)uN#Tx*sN=aJ}DUeaOK zG8`Anj`(`*R{NoRc;&=}F+uP2{V&kz97&QAlHCQPFz2q^vEyF4NwN?T%V9GGcYE@2 z=IZHdifYvKraMipji@Xq>RHTQt@sLE^GB(PdE+a4rodJA6SRT1Zvgn&gW8TuUev^b z7f+`1#X?3pxp%h^W8OstS2UdIB95T3Tp?dDR*Q;I82I`EeA3bJ0@KEGpl_$&a>~Sg z)@5ngn<>(Z4fxdd6zj0Pk`-HT#{GLX_AmCHiG*795#Zlb*(y z{*CLP6G6Sk_W!DaLdzG-7_HpQqFXcs=eDp_W%lOpJXE!FfNcdM0l{* zYqs#j7G1pE*y)0ovuk4)GRJ!js|`V@pV!ECflc5PE?xaL)c(Y&!Z-JK?(S(wL42RH zHE(pH5~D^i>BkfOzyk#rD&g4DD#&6uIX=(RkovuPG-|Tv69JW$|i&$T$9poAdyvRwz33+ z^VwVAeqx0#M9qQH5F0^DkQ;bnT|5$VQT2x407yaoPZuOTtE}c6v(q0Y-D61cw4>7U zB0DY%CPSdv^Oe{A%*+C3B~`-4jZNFND|Dq+?mQfy!Jw$}J&(B2<*$ihw~1M@qy>^c z;z|slTc_a`ry(8M1i5M8KnRM5%SmSV%a=jq!`wxg@Ec5sC~fhFbTXd0^VEM?|C|0x z9usJ_Q~dGnX(LD1e+>%-eu&)2U!97EYB*ZOg|z7WNOufOuttS$>LQXH@eh5M=(+UY zXcGThFO!kmJlGuvu`qZSQcAO(%5V|PdeEPmHgHqZUw$P9npGx1WqiBVxBpfa=tY95 z2{5*`_$CTZ#p^Vn6hX%ElyB z%H%|xSh6<-nHJPIqP7V)P5cs_khE+4FO0j=;xIt6rd(>akF8DrOWt%`3%KS$3#j(7 zy_y%hv}_S4N4Nc6Z?oZVZI%Z*SN`5TE*#HVg2)SqSbAvSQLUXG#bji-g9FIfql;*r z5Vtm*W{(jtBZ_BkgabM$`?H4%Xw$vPDmHiYyvNu@rC_I`8lBM`CbB7}uK8kt{pEH` zh2h!2Ux&QvI!q{f{)P!RU-W9{qr;m(;C+5Cl3T^}@@LQ7$IC~19O^6%fsyziQdQdi zlEKcw;d!#Wb9D4LB}2-`%X_@m!6P8>SPem?)49+XB=5A3U9Fm{ytE-6E2EL&I4#IG z%583MR^lomt=kBb;PyPTM>eiI)I(J0q zqhr@vhIb_k-(T10s&FuY8oV&8?o+wa;8;wa9&tLEiqpHfkyTbD{b{<;s>BKc-yI2Z zJp`<@cXk~jj6$)-aui77g>IoAf~7xX&JPUy`f@(JXm@^ZD?$Z)P;6_aHZVl{!@mD_ z`Qnut4{7&9luxRD%@QXhGyzA{QTH0(mvo8HHM!EUTNx9hi#j0z+fxWtrF;pFqE_l3 zH_<1i1meiwip)4(%Y^>3!*cZ^Jlc>EIb|D7!k!B-rg2d|8PCr*OndDXO3wBcP-%qb ztNIZu@+oFsU|3KTrOpP-Q?j4D#dTLDB1z^j=h&CWZlZWTKH%iVN}!?`HX%M9LhPP< zWqdAfMC5i_7X-shO--B>6coU}%S-6U$RP;*Zb8O}mlX~GRbEvQR6ju}rlf}0$4rcr zpBNx}8VZo`KSo_bRUmoXp!;ecHqzNb-V)Ri`!_1#T7|ESiq|j-s=!o1N%uV=x#n8$ z#&We+Yjc@=C;s1GwuVddp0V@y{XE^& zme8>NnBF{}(?o>|g8o_iN|~Y~rl{iC4l>|4AItdG0Tfd*=C60Z3D?eG3*h9I|7QBU zgk#5|)H?bz>ugknurOl3&+p{cpTzASlsWGW-_~-W&-4;pbVS5rvqM&n;YrAw@V*Q# zJ3LECd%k8m%|rvi{-e^y^NZTI65P-G+J-~Vm|7nbsdu%N*Rank1ee`G9vh~O)(5)L zEJ;5}Td<*$4J_?#e#_M<8q$U;zb9Dm;oCBy>M2AG1l-iwPwY&%sWE&p+u~yU#*5XE z;sVjppLH?I`^8Z}9k*pA?!k9Evw@Cl!#r@g*(3Nz#OJ`D2`PEiww()BN|ybGuJVFo z2PZs9(|6*eh)Eg#E9nILPc8o6N8>e?Q=D1KY55zQ3L6Q4(Ze&mEJM6ZjW>Wj!rI)o zmZRJemm;FK@NZNb?V`m>tBJSF%znWKy^FJaSF(ESu-;t^WH9q^qP6ktgs zyF2Pgf#Y->URz&lVk?*ZCR^ezPJ9#z4hc&&qtK^Xo_bM?G<6S04;2Jr2`#gI)X!l3 zEiOL1L-azcot@VA+ z)@Q~tC*SX4L(Ru65kz9GCrl>6L`j^`xL=}CO(#?T#M120RQmKCkxve;&KX~hP~8qR zZEt(q@9BJqlW=T$x}mjIhY>HHSis}@aL%dGzYKHwUOfL4r|bHhyeI!WUd{L7`s;56 zS1>Y}TEmOYf(0G*KBli%lopp)JxNJJb>Ly1w3%&w4E?s8KAxOIzJ{iZ#n$@h^K5ZS=*O&c!Jf^E z!$d?xj=0VGtU-)2-1zFH6*z!7$%v|)fvfD>@{ezYbZN$PE|5RqSfuOU@_p=hyLm%_ z41%i1Iym4T{3bYoqs`DKoSVvnJje(n==%TKfoe``8Me5$;lwX^6hPd ze=duT+h+NN!fKWG$#I>fJjDU(XE%QTz406rd}dS>6zK2Ybr1*|D(cnswSQKuL~P?9 z9Yx5t(%2$4=Q0B#72SWn&o68h9AjzSV!e`Xk|G@csTr=yaoj`~aBQg5%`<0u&3!z- zj3RzE+FaM^GvE~%lIt%?~yH+-uNDo9nnQjd{{ti@p zBpP~HQ=xY>8#DA5&VO!z;=qU}-~fzXRhkrC(2AmjgF`y|%?9YEE3{u8?FT@KGq(Z7 z1Hzup;zrg;#GWp7$<(#Qtgq;6X>If#qJSN+vV~WJrvgA_77}= z+hg^%jLksvVD_eb3bxvP%15HJDw4Tw*P{nK81HPOD!)==osT~0T7nvdKddkLZz-+j ze`eATBldB;WMGwPwTjlrk}U@?+}-s;>c4e<_Uz_Z<#c={&-gXK7BE-MQuEAM#9_s+q}0D1-HMg=fc{p`dk72 ze=@nTryBF%Zjz{Fm%;D#sbS7omhFap4{8Hf3B2&JGOaaC*~D2pG;An~xawCsM@gom zg|RWe`JcaCRWy)=piaT&%V96sOsJgEf4hEHw*{)5;i}F0cFYvYfsG4C3XFNsQd5U% z5LU;5Y99(pA24sKijxhR-qSYM8x_CEBH9(itt$k4xR{~}jJ8p>B{C6zLGa%nPwOc# z{i#KlTn|v?_+{3s^;)6CSeEWidpOe4b*Lj>8dIve7B=MQd_6i#{mwJI#1?@S9QdoF zc=8T_2UZp6+n)>mJ!R)2u-poMi2OYpW0Wd*ll^iEOs@1*w1lVZKFXVKk5A4LO(3^2 z`gkE5vlUW&9OerzWbvO-f8nGryb{AQzTB~=He{rVKT@s$14T^qjW~dONZszQ#?Fh1 zw_nFQ#ECL+bY3#G|0fH8b7(N8Ltgc;0Vy)fAI6>2lctD<;C?*(W6fVuNI?}G!I_e) zTf2NOn|)?0&o2^E2_oxmzqZm+E6{vQNl9sJboKG^@%Fwa4SYHe%Op=4W@l$VK0fCC z5O;A7=5QMt*ola8IwHRRGWJYdFjVxq_GXDa-T%U_*3(!yigvSq;7El71*_Qpo4+TY zvkleYMQOo@j3mX5BEfA8A=T3P7&StBQmQ(qR|&dmEc%5ZotYsIr}LfTHM;;trbg45 zSobjpLXg_cqj`Id zOQ-$D^w;U9lW51}Gxt-fa?iT$s5hIMd$}3&XUgJZr<0r#sSj&%&)SQ@(;05xCXZDe z+HX2v(g;bnn{12FgO&ZJ=e#dtDW!5@D3~wzRE0qS*#OU&?jpmqV0@v__4pbD);R*6 zj^g6KKnV@>4@d;PzGr8zcKiFfx-xlo)XLGPFzYvufxvxA3_eRiZ9h-x zLQ_6wIjXD@2nWcdzm8M}b8iqg-N!8t=ec?07OT3cw$ZzpgqFp@N7*BcNzlz{8LQa6 zO9~Ll5OKD;+iC6Wh?IgVi>!G;eil1qYVq)-U`~+`LKVS?+C0%=tcKB9q2Ma{0<=#J zvx|9P*PFR*@i6?dMkM?7cq3D48wWm1A+-?OUUwuEO51WP^R?g<*< zZCW!Mf0E!$cP>mrxBe|a=Fb0G3IsuSn2kEU<|;Im)+3a%xWPQV z2WXOm(>!_L7{g#JQP0x_lsYi&R;>7Nvrn0M*XAf}{kgl=(jq*txg;F~130M1Yajg}-d&#E?}_HrZ4PL;uRpyQ6N$5dm+1Vx^~?Q6Pu~ zBk;*n=g-6fe@#38^+tre0}5}S`khsQT}K7zn1P=MK;mDdCU*WMIVol3Xgt`Wf(`8YaxJMXzAx)(dUj%k{j4st~V~*btDHr6(fJ!cw}s+T$?41@5EUB!CTxg zxUj5uYnpb#a5Pu2hA=b9?YFnLcel}l0`uL$VQ650mY%)_Y+o*}h)~P#ZM&+AjSjmm zYo7HcBLpya>%JJ^TcDOn{cJH^@RmZ9{Yn9V#>q?ctEA9YtJrqs$JXUW`yXy|zX(rh z1TyFmO0Tm}VAfBqn#_Xc7E#dfol3qoUT8J(UH1__wqr^azGPJ~^=_|ejBDy(D85UP zdDohX0$MBQWwP*qhw^Znwp9>T#oo3K`Xw*n^}TnzPmdp3w~W^gMInFlii4%t>C`FB zDeJ`PFqH~dA44fACFGqv;Mb%W7riJNCFvfUA~av>zGYx`1a;lHEButBr$7YkMM=Mn z#{oJ;xt)NL(oJY-+^m8Z(U#+tUEao;gNh)R$iUzslB>m1+;346Pu`9gQB1O9E15l# zpt9&lu-^m=cY9baa3RcSg!{X@617TlZf@7Z+0x&?e}l&!gjy1QH=1;yiO>f^REOQs(z4mA zn?zq~fbgPfQ#8oyPqk_HY zo)t?70h-lLBM9I$QsyZ4v) zwz%56KdXt-WNNC+jF5%7snde`X#~n^SVbt~I38>jNm@133HFfru!4_z*%(4?hl+DR zLTE8$+@-L3UY=e5r%Z%MU0oflf0}^jL!I>!kP6m&J~~s!xUP5nqP7(V;)pWcb^)V) z%$}$7~TC(DpZqD;>VIxM5h08m0;m1g4PSrSz*IL12?N}LOCgV(AQ`NP~ z-%|C)gZ0~Pk`oU0ta?qxW7*)wkP4F&w<_N|p}IRda!w7u%Fv>jSBjwhniM9vtvZQl zCeB%-VMPsnsP}X$X`hPrrxiTB1ZvxE(3%7>ogz*$5RGuQSe1!0#oF2$ z3aH9q_$*mjSrG0_=1F5!vx6ds^Fgt6lZ0=}@7qA5{O-dpC`j#y*^Uj~-_NCC^PDq= zH3OB)eML-{k0A>rD!zbv?FV7TWk?axlzT1g&mmCUE-b$^qyf#$I#Qap0`lusbZQ%> zAz1H~x#PpO1Iw>9SRuI<#hbBQDIf8a1!O)UnJKS%_F=zf)X?Kp}q@|Wbkjn4?? zA|p^0#|sxkH^TW0!c9QXTNOUh)w;+v4wL2oNJO32A-r)o$+Vb4-%7W zYu=zbGwcnEB4uG?LqS4v{TiOg@9wOwz7T*$SZB2m9UD8H_l*z+pVg?pP%%GGHX}VI z+xV1DVI8HLo#4qq_K{~t@H7XIcme;3-J9|S@&N`tz|;sA1?SB*7lX5#gt9jiVd=zB z_s>R)SS}cGr_;!4PJ`1Vg6Y{&d>!YxSYuoic&*{sPu&Ts62M<#KIuEL6Q!OKXf|yp z%qi3JGM1O5!ydKwmt>3yek#MRuVU7E$t<^j{odH3mLGOq#H%Ob>r#3c`}0kKq@;0A z`+QtOi@7`L3*X}Iy9|&tGJB6H>k} zsS!})&5n}oAFxdzc255^@5n+SNIWQ+a3~DE#~D9sKTu7#?KdfDj>U{or;D2@duqt& z7JQe+Hn{dLgZ9@UBQ>^)3cLisZJ~UN&9zh*FhI&$!^VECZqt~8b2_9`T6^+J3@wZQ zCy=v(RXhQb&a1tNWbRLEen9H!@dCXCkp&&S7Nr?H*YWA;Uf@gV?xq%}hMv%tVS9Tr zUkI_1!Z;p(p`$6*@e+M(BJ4pKjOEPRv%P}*x%FfZn8_RwDTK`BOna6l9j{<)!)w;s ze&=C}sGq43rwwB(r^i+M4z-UimnvhOYyTe7^k9u=!i4Y0DPt@aoyCnFWh}&!Gy+uo z!GXmWXvt+osf#f-Rl5eErB0PXr?JmZOeRej8mu50v0bHLV`xjs(t@Q#6aX(+Ud+Ot z>_HsdrgB{}%fK-;psy@@=tZf7O>!l8eEgRdT36nii4Hor3i>wg{B&t2e4V_A6mb zSaly(WtD|$|0-IBm1f7&l1Hc=W>w1XlAFtCUse0YW+$1U}YkCaAd>wkevnlg!G z{WS{iUKS02nYzB(;h-Y>Pqu@OCRD$r5V=@Zu=RXnBnUq4vEgz-cG?u}qY_u$*>1(} zO)bq)CZ$@d>}s6QARWX^5)zdd_a>+_ztqElZ;K=f?y5!`)XLiA~UP`b00(@;xiv0a|>{!wT3qnOTi=%lwbNtc|xPH89H7k1YD zR+xt0+T&yV-HV>2c-TdPn>Anunhi%hn_Rzmr818(XCJ-bkoh=aBz}x zISMU)JFN|FWM$AcL6t03py9a?8pF6rHTJ*i28d z-pbA5x#C?o8NRYYk1_UpQm*Kz>Y?wo^4k*?4>khG-c&ZmO9Rw|HHk38Q{Kd=FR4!? zusE^M>Zt&W@XUkMw=25Pb#d;Sf1wZ3s=iw z~BH#Qy%RE7n?YYX%AjGL(^t1rc$#9Q?kZW`~_j#qkPnm zy>t((TVC~T10=i;IXd6?)ZA@!~#s%3NSA?k|4jR8+zC=^8WBmcJWbjwyL+ z^K30`p8tNbO}8VQKun~oFy*Ho1nK-vM3onxgFE+Q>HE#t-ub{pVB_|+h%y=xLT#DE z{E@gIZupasV7b_9WgSkFsjQXYz3Z%+n}nhFJu}~;NaEjRFHegw5hPKL^9Y(gN-~t~ zw;h!1+#3-*6VT&N-)3tOzS>w_^~lm6Pa8Y=r%y$^*=0JFFWT(7nAr2MP^9_q-Cocl zjETO!ypmE*sV~>-$7|ZEq#-7s?#u0<7HKHb`h{$TJVugQU(HS8 zxKE6P(p2?Kt>)Q_^Hik4Z=l5f!qk4%!eE7-b0wAb))UpfF7o-OM+y56*8TyH9~GN- zBxU#kPe+Z;TB7o`i>@)<7-{!!PSSqxQ%G`VO|v7bqe`~+#v5_^R#_VS@RNRi7L2*% z{Ft{JJM->aVQ?W;#6@wZmg|7zr?3xxpQVu}bmuK5S;Q8P3M*2+gcUEBCE=uy?Al(Q z>*?-W*}w&462I&6Zo49<%I{o==4dYeSVKe*VgMs0DXG!w#&P}wII5YQgb4u4Adcv*jlP`a`LS8~?X=#BV?+fBF=-`u){ThfR zB$dup_$ZS)U#=_YiCfv~*doA8OC+XWJUY1g1rOtlI=)0*IzSKF1(k72_pt47rA z>-ER`MJ*4NteB1c4*~l4kRpY-{f$-ws1@Cz-vBoq9WJo+qtKo%z4XDXGMZA_K|I+p z6X5DCE)t^QRhGgHzzom_TkJ7UEH7dD*`A8+T>X4#PXAA6xH^~E<$q|SK zIJOAmu3RR*#4JXo?@8DO(ak(s2s80gvP`8@)bBN`om>cBzY7X;u_t>^5o~Ae-62mq ze=_ivKinP72A(xbpkjS}o!|ZDgWu%TR6TITs8m#S#GeU=_%|ClFQlk2Md80dcT(EQ zLU@Gyi!?2AQjBv_wqKx20aHF zbYRV1!eKRzWR`YjZ`|4Q$1JkJR^dYCIHkyCbsKngV;OtkQ4M%f$E)plhju11;SQUdT zlKq^%jbOyq?~Kz@hDmp_f5bb5VZYifF77COWK#B_ppK&_?7&yL82_d^G|kjGSz>j+ zKS{R+ZSU);QqC#?bQPu|^rUYWiioin8w3>m7`6Fw;&96E7?PEzXrwwp0+Ai#`qzuxc5C(CzmxrdzBC!M5go8NkwNJ4ik z=Z33jNoK7uEN?n}X{#neC8;M09|+%-)Gt{die{Gz#%RzBN|NWd=xARZ92$QeCYX(G z<(?Y*o0@0ws?OdweMedFN`l+LLu$WKbgF4j^TEsIJa!HO~y} z5kx|9-)aJAY47k*6ZDZaHG%$P8CX9|PfydNy}?&wVtT%hcz71z*ywKAT;aoRGX?L+AcblST(}6r17= z`NT^H(v?{J(2&cu$8{!b)U4*;b_k-;o9MTDzM0;^>MAz!L$hdlbm#YKUw?EMbAQn_ z$Q$nb^gh=K#KgtOb&{H30~lWK=N4dTs_PA%$(JSd8wa#M5&@3}`)x@F2M2=iI`GS< zPsN9G<%C%Uf8Os%Taelw+DH~gV?tKzQeSYd8!y$DuVo+H2i8WB=?Jr0OFpB93}(C# z*dU5N#?5%BTHmBM8ycv^=Bun+4}L3UXAKqs^(U8g^CY?}=^{j>h8ty!I zw+B*LCjq_Ap7GWuc=pyVsX2JAU^Df=+#JZ0S4yuCqDgN)WrP&yd?0DxRCz(e3w-{; zandYhP8IU07nVQr5T4ib;RcJNvl|lCu87IxjJ`|0ujRBBpq{USV8)slC1*dToHjvZ-zi(f5po zuF;m7=gmmrakO(z_UI-sBt$9B4#x`k$(etfsuF4oee;7W6y3s%;NJ(K7GO4@^Wmyp49}J)Y{rw23p#}g8zQ@fTc#e7qqCT zh)$=vzFyGl>4%8OYhk8$RJ@__UyRgTDvh9=R0iLFmRc^wV@{_}%>D~gCnOYKI;V+o z5L158_QTH3Rr6sIaP&U^`-l~cRBS`;#$^HhFpy8y_J^D zTNF{-)@KlrC{Z8jS7mQzfV~VRMCKxtbt#4Is>u6a$sWFnXXxYZ&JHd%_SwA$!Lc${MMp=+yN@!Cr>kjWf%|*J zu;CDU<*EE*g8>cvT`071N7dsrl1HZly6nk=w!69#m_M znpjgWn$vxK@t1S0J~6jT4hvd{+KX(R@ild3 zH~S*8v$NkqfX{DOXeelxBO)QGR~q<}guYn}|3X<9HTF4OMd9%#-ZJ-MfvSxk3Y<3? zhvo=y-VvZ0@J$qJ3$PnD5IIdqt-Qv1CzF%Ppc$$3dM{>?U>-0b+Mzc~0Gr3g#m&pf zN#%E+U0Lx2q<<7#`f=0yFH6BxcX)jy?<`es&MF?fv?_=beqr)l|M}v&JH0{RCx@cC z$Z8>s)`4%6-TO!6I=nMpd$RhE^{;Hvg=NvYOm&WDs9($6VZiyG18ONWn$WGFDgu_ONNlk^cqDQ z0o$Jxsx%<&-{)BQ(!Bf}c7|ZlvDB2lruyJKm{9WX>-hEV*G?rP`O58aLBfO6QPX#g zR*a0tE;kY^X%MgClr+b5BxX*F6if**Ol&Gw~@T|3@;JAz0Mg+^3VYY@PEyC0}A1~rD z>2&w>7@L}U-)|=%W{!`K1DOF<#4n_Om8!$J!$^tY! zV8{ph0sue-q8u_ZGQ>>fVin-3ogsdKB8?dJ`R7j%pd$jW3cpY$RLMFQg5HxFWt63% z(cTYp9BUi3krTrX8iSgs8N`1C8CZ)Ku5}OdAo*WBtNXGxee+9QO${3d2beWW)C2D# zYT&K4u+UgkB*|n9f1I6?eeWImK-7??$b^Y;Rv&&@9oL3^WECZ^@$pgVAyI~Gf_95K z2@&tmREk>v<0;;r&EFvtUcuy679ngFq{QBQYq9)jLjTf%e0pp`Z5_iuzy9wU0Ttov z%?zL`0`KHk(%}&dyF>cgarNFcMQa4VIyD#~ ziwWDU`bpjcBbq5UJgk%3X;{4=owxid2>bvM@j;>4ms_jPS234a;K{&*wA~O*9y;cr zpyy!{BB9dlPd}}u&zk6Dq~i`uw~OP?PC6YbqiY55`u(oCs^$Nq+ad-*GspjRzR`CMaV+U#jh9#D#2WV(@NN`2TiHI1rzkBu9hT`O1Z1 z8cYq~qH46w78|N^AkMxdxV<v1XJ?i9#@ljWp|}IdLk}c{U0uXvcWi= zMf0}pas%u{@ze(Mme+zFewMs=?Y9v4`m-DHxOsEiDvim@>DN=aRZmXeHlB(QG#S6n zgwkEz)X5CQND+xhCzKT6(Q|tUUK*z*Z~oM$ne70?)8!);(3qO}rvI~24S$!-`$cc)o=gu|DW)jr%}pClPV?j0C`;wfwAe_$;-OnpozMuHqQ|K0a!>o0pCM z1y9RE=V7l{2xP%jpJlQmF(CjRO7&DIu{RZXL(M^m*8`)XpV9&)U8%{`+D}t@e4r-)irw$L0MSYC+H&_>xZY@fW8l?r&ZIMrYNpxDHntoxGlQW)q@Ub*#D9i3BxY~UUlBT|D`5@4n^Jqi#1 z#o|0OHrqwFw5;kw-H>XRa#%mccXnHaE=-v~7@STkmh`#1cFidux5)U*1Mp%b^okTH zF2}()#?thhB~8$Z2@y<{C46^s)UBwc=aa}q;Tb~FofTzzbFMc4$X+1vScn>{Ild<6cKaOm3 z(gLIb8;*=G37DVT0xd<&wSQ8~O=(m<_+v(oAfjIV*>tk0-NYH7-ogo@L>^LI{Eq*d zPt-b_rgv^!C?+aeUGX2@VtA58J6`vhJzw>FAn`I%hlc z2JgZ!YlanF{Wl|~FxL88wOtFI{Ak**#%-n%DgqF!r*!L9THSB6va?G|N&?V{4-OA4 z&CP>1e*l5^+C>JEl{FRn_CEXvjpgl!>-4oCB)H}T2I$w^7?-@eTZ zMQsOj98Q2L4p@|8tJn%5S}!-)fxf#`JUzGjjaBt7Ee*|ly={hYR@0aY4K^o7VImw{ zjqd=OVZ6C0y4?=#Q9c@V+%iqHnamrjVZkl=W>J-mT;~Wb(L2}!_B=+yBf>+|OPY2B zvRc>@`49S2jw*0}bS*wgZ|&@$6Mg=^Kbeb|2?|V*V{ixvtoHTof`vZT2?A3OWKuZ4 zal4WE(O?h0A@zuijm;TxwXn1VG=IA34m&MM*0TK?>YU~jt=N$3>L`I5mvt|e>Gs@@ zR+u46R7sv<27c$FR17amN+)?!mh?q7L{HFdey-~2W6TW$ePe`w#hY@P@G1W*Wesj= z{a*>|yzny}|8?4PAOX9@ih}thBhYrJyPD zYnboG1|>D0E{_&1q`eod?W$-!ld(nndNyjU+Mb>qwX=q2UlNCyj?a_>0Uu! zc{6~mP@U;;58FOQJPfmn%XIi?^L!@wW9bhHJSi0$5l^W@ljBv{JT8K)*mG=!U$7_v z`r<2Vy=NcZ#0h2wAS`xug#VVbhKUiU)M5m7@u%%CtoUKTy2lSQuxp}c{kNZ=A1H#X zfE^?d!Bs7>poRtqbK7kQdw6&>x6_(Y5Q)`(BS(m1cilan-cU9FF`P6FaCXVc{f!7BrQW#uIH3$E;@a8zfjeO|bC;B5Mj;{}YU zfyhgX(FM3r{?FvPq@)Bm2Je-z^N~dKZ8{wn8PuZd@ulOacst$+oPSc|bbBj|)iYkb zzY;0?;`Re7@;*Z3k-4|Hyq*JMCLGk&mcotfUbVM%2h%LkZ1X9#v*6%W(^Mr>wM0Q8 zcD6b`L_BX~d?X|;ce=-)g1xSpEwn6AtZVHz{6uk=B4Gz;$N;4kuuF7ww7Int=%1l_^w$>XsQi*Bq`>4P&;i~KzE z6R;0uwHN`qQ(%nx^dZFWvWN5*n&T!E=N*>PnVHbwzfy-b*u(5EDjj5W`geT2M zdG)Z_p)63K7P`5)VPs?kG5#we7>oi1qir7O@R@b&p6hH>Y^f=)E-P2L7%q15E8m#j zAUaqJ_-Qh0uDRU{PCFh#B4m0{y+m%D(`_YKT(@1=_w6)_aQU+gexNZh4sBMKC7A_S)CMr88^4H3|!^78Dva|*Uj?wxfh$8JH0@xJyX6y)7bi)po^PJTuMeE}fP3s=&i7wx_H5b9sS@EAHCd6&evw}&)gN^a&SO1Tk)#kEG_ZHkh_wT7 zu(4}^4S7e$7ifF!B7vHU%Ir^Acx=wz?ry7vDv>@BIO5cH{#z5#q$_jPH#XIQ9xl{Q zrz2Jym48o_5tzR81^ub^l69^!N@Ag@4W~0vod5p5T^eRUD>6#)!N0(NCQw90%^pF93VA zxL%r!rb=YW85kPAg#b<>$YX(BQe1$M?;l0lG8>}Cdjb}B$Eo|87qnBJ#{3k{HxKjn zqu$fEIx?M8%lhv@)hV64lZRr`^L%EVP zB^FN)-J475z!|>RXP)sKmJZ?Wr;OCGr#|R#0WU3G&+D1p4Da8VRh5;QQ~Mxs&}sOd z34oojxf8;Lsz!`9!CMc7*p&Ee?nEVfS=r*9TQDH)>O=Wp*u=2?DLX4QCSl%hBg%=P zyUjiKC`Q#Oot@ZEqr&kVEQ8AWLr1Lxw+^4)LX2yDxZPtSMvif|tTMY4;9|0Pej@pI zf}r%@*_Xz1eZJWTxF%ri5ZJSum-lXrq|6#U7fJ4PhyG4*QSY^vaCT}OsY)F?%j;8C zM&u@j#PPpgO=<-QX&>YEx=r=xtd*G-F8}R|uLb}~CtrKkfhf6~M4GDtDlDpaVcUv0 zWprd@B#_vI5pd9f1BD#I5)KEGS%xcivLlfJ7I%>!q$PE_qh_}bbxO^CZ_;}4iPIWI zj8z&<3s6qePDZtf8!#@5$!Aj)?XUydNC{cV=Y_Vp<944W7MU>!9Qle~2=*KwrVACl zU!X2(_Lv21Wz9oNWaaiFjcQe#E$?zdU1*Enu{TGQy@)PR@!_`wOe0&y| z2=#B}-ki!V+Fkv8o_LUj2kDRgc6e$@(JZ1bOuo0pZk*ljxUV2C-oO6W7jVlYa=0>O z@T7cTW4==B6ZRo^6P;?Ve7IUqk?-ee6uCXCP>>0y*&Z}CZ7;FxAtd~oq8X^F&ZPEz zodja`v%R#*({H2}YU9SWm5;Jt&pK4X9ftAX>ZXxW>p=U97MvD2@*(0PiM-8aF zjM+GERJagT)7{*j%~=A4k|u|=R}gz65Zrc^*;8G$axP|ES#)^Ij4&6@j5r80a>@zl zqMZ^EL(Jvj1s-$MR-|HgdjBL(&?0YQ!A61uJxHVxW;_aHi?x6NhUp<`5?#g zUR}sWgrEB5K9Q|36yXQcZ%Z!%UYLWLmkJh`WRXs0|2hi)31cU?iOnYqHl7|`>Iv*L z$RaoG#>QndKyXMkM-LQj>5M-AMiFk{2=sFEI!&h5mz0T1!LeLPng(q>ITR-+c3U$n zPfyZ(UJ3tG09cPb{;kIOlJ_ZVA zNEzBVHDBN|;^Oe*PZue7MRW5QrH`j9)nG;5BVdO5ku!wJA z5oi#zn^hzKsrXrTXl&y1T&oAUdd0i)p3p!tznz{RXy_Hi*4Sz8Do){V= zdI|>LhJUXQJFmR=v=GLK;0opzc)o;v3##?GN!^o{_`(ZR;0-cJ=%PCoo!2y0EQiXs z_`TME%x@S!HtO$hQnV}QN~2GwQ%!OJ3~1hwrs_rxdwUX}4M}4^HvK$VnMY61SWqO- z?Iz!^23oucwcjQgz{7~DqOtzXnHvs91IvAq@U2>MThvT2XaQ9qLOTNCV8q#Bq(OhPHeOA*%JC#~&Z!QU= zfedb<0d!3L?d};>==V2yeu-;vUsi5YG6pXQu|CKZ%3KyC4ky+_`f^|@n*9^4~?@t1JW}XbR^OX7*lJj{$f8bw) zC-d9f^ZrP3+LrE*c~QIhbfr^}5dJ<5La=AcD`%2^1S7aqExT<(d_Qlu+xvt2T&+7t zqp%EpxvwZsJU?!_`(Z(e3lWD}?)JKl`&Q}LgW69HbU*`GgtF^(=BoA;?4;iWm|TH$ zY-OZS%@tzT2@{=R=TFxCCptwKi#!s*gbt7H=xB~VOqbBmF55H_x-fq$N>FAB+jv3n zG|?gb3(;|kd9)1WK7drC`Oy!`%!sc@&<_8lSO^UjKTTfXelS3(!Xb+nUihJTG8GG&T?Uq8l`z_SoQvq39$!TO&C$VqTr3m&D}o4ZD**- z=JyWi&23gB%SV#xDj~+XjfxL%vfOvd-S4YDsj?JqVLZ6y;mlA1elpDsSAGVG%&6l0 zXU=8UFRDv0Vm1-BjO%wPba(;d7m`g(IC z+G}zpq==9!O-hQP?7VJ!IO3YF(0tLILi1mCVl~&oP_NZX7)hzshu2?_UNP$t+?lX1 zAGqRw?j zP*z$GkW9IA$8i!pRA&)v=aypIEI1;&N;gcs1UOM4occGNowWV-#)kw71|hpt#-2Pq zmvhB%h>D>RyxN#t6t!gQ1u!ICb^>7r)VY8H5+vPO=Cy1)Qkef3V~<_zT11zK51E1c z$E3>o_P@9`Sx9Ms?hD7uXgE?}Gns$0GyI1I_f)jWN=CXpPG0BHt-=7%q(C9RM*WHo z4B1op34%_wZu;5c9_Q-jcSY5lTD>lYK3yw6z#Sg7re$yUeRNC9E-q0y>#6Kh zjW0yhqXP@%l{U?9aNVepG;t#v+e-sCZUK%Q^Y3*=X^C@(5t`f}{%fc?sr+bVhQW83)XI+ymk1I6eEmn9 zvJ(&EB0iw2kbforMyZ6L$d3$r%PNYzJvF;2@++00sS?$&@nkU9`_qe_rq&lhvhLmjXBnY{8`M#a@aeuU{=S1sQFfVWLEu&WEUTqv{w;(eMz za=)^$Yu>~0jJ_;lw>o)_2P#A8mChb%ebwRXjXRqgWVGn9P|K5TrNtz+=ft;-`yOuk z%DxN4X-gZDP#ah4Oq7rpCQ@B@ zf)Avl8RJV@>#KTLHy8b}2L9!w{)N>jEi&%OOPCqBjt42F1~xcBGrAu~@AwN#RDNQT z8|Hw%A=VO8UOcWc7b6KwBc{-y5Da|Y?lgQOdQr7(2!KK2t&*Fdxpg{j=f{F$L4DrK zqexDbLFcvFygZSG2YGz#<2uGZPug=;Z9M0RKZy4AYi%_s{lk!iqwDLn6uBAOt8^;v z4IZw{BCbspm;{+>_2dD%9)=+&EmLV~5nzC(*Xi3Lo1VZ(WxA$axK7ot+2@pcb8+Hp_PFU7dYCgT49L z!a?r&-uOj^Qdxpnsb9wD|D40f zX@W>5Ygt{J?b{$thtx%%I8C~SX17&YwLSGNwEGUmyZLk#2@vN_`VD`h+8u58>{goH z?AT{mR8;#97F(t#BRN7nt^dp7{)TS<1a zR54CLn@UwGz7w~7uJPSl$lp=A3zU>Y`I#a4*n|bT8VoK%vC-eixxW&$t3*95SBWYg@}r}8DBjgXFK&yGVnyiAa*iy$otZKXfFqO7z^P`b4d(O1wxus-D73igs&*ow8tJc>0WQahQ*uCkvKYz(| zCx2bz=g5+vw0%RC#BVkx+dR*XWF~79L zhHu|EO+ZLkTvTLbZ9PMb_$73lHfJJJCVQ!3C!RP#KYz!WXx+@lT9nC?|6j77KD>vu z%;0pE);+!db!E_1S`=TcZO$ zUE7?*0*=&+#}~e~H1+2NT@3FGf(5{zog^;dzLL9VkA}hL&fbH67c8vnasE26>-_Aj zqNe5vAb~v_T)u9izm63#89v9+o%onk$Eb49J7r+%IcE%9MjMrz6P_#ZT9oqe2t^#7WgT5)zJiQMdhcXA{J! z_!fqlL{uJr3CuJFj?vY6e|-R33{V~aL*D>I1o(yQj*T#8Xl?J0JBDc%7yzAcqCzlE zn&br?K^Cmx~L-u=N5DnTW?zkM^7oh7(E1^aPV z>GMT$5}BMllia+5&PTuvZo8)*!TCMDU`kc>Nu_;eWc9#_RWwMaL-@T2%+T0a#P8qn z0-iv51pcl0qu_yll3b=BD$FKj{Q5hW)%Z~iN1`+sn71hs!a_>(X}HJF_j2BT<562J zy5+{RoFNa7i}b`>+1#etaSABFBh-y-ly9Nm4>QEE9)d`V%4&rrDKirKylpS4uWXnT zxD7B>4Q+fF(P0`d*ZarE3_))l=OLgm1Fit^^y&r$;M$s%Mb8=iF0fmbHAJzp;L?q$X3Z`oBo>zDsoNR%CL!JF zHYCs0deou@N)jEmqNOt071P^mD|5-OFkzAM(;d{(iOk;59#P`#fEWj4oa6hBgoG)e zi-1Nd(2|$eyQ)*%osF0Y_Wbytb6o*JFTG zM)tzv?{}O;a^F{)n(sd&;G7Eh(!u`S$unbp54{)9ao3QA;IRz!>&J&XMRGKR0R48E zNlxu}o?MP^FZ_()EB!3EWJ1UZI}VT9$59?r_VGVIg?bxpU}x^;)|9_nY&-#`!*PS@ z?>;_oevzEp*vS47B>d3@8ovS8(tI)pJT_N*nV%6*f7Ug}$#XhCep_HzNgZyr=B}bx z=SCfbGQFwue9v~gt-@|4A#FH88;Wq3Sg8J@C-qQvu?TF*EqXJM+Cu6{-3Cqci z1J^t_c(e+cLaao^p&VrqiS-(BG#IfLFW$9YDZtTv%-47~+EY;|(}Kq&vz~ z#;PCw90jN1?Av9N?% z2+L|uuXL9lnqrM(VG}mYPt|@(X|hKFm;mDzSb9#**>MptB~P#A53tmN z`6?>Z*^?#4EwzjqL zOqKL5b!kIGqGYayhI3X9hKPuWBTri2p|GVCE-o%>k~B)F3_E)}h(K9cogWrHEows` zJ}D_FLaU2v3KFo^|CuQP8um~;1Dm7JmFYy5NLg7Kfr6GZANa7qLH7`}oAvg5_Qn85 zjGTf3=$>u<`wNT@1d|w;^*djE#p|vSBE<(EJLVNE>VF^J2L=HKrh^is&|pQPKm3G+ zfnk=w>;SJwqy)|`h6e`R*Th)II-VK5`>R_BUazP%;0A7 z>+0eqb0z2YPf&UOQO!+HgT*2UHvyNIb#mf}8SulO75(3DYX*7s|F{;oBsCJ;|Nkvg z?*J~y!9?~OaQa>8KvFL$NtJ+3Ciq^T7t zvEElGHN~f-G`n8zfav$Z@4wKyx3MQDEj+8|cQ7~b(P(L)sGHq8dVVt|Y zljk)Db4*l1f~vMQzvG_VempoY!5HAwR3i9z$+`dN?cOh5!NQZ3szSvD1>(}uL~JIg z=dlAu=KJodTiXtW zeGSn1_VD0?6o5gA0gR#Opd#5-64@{kC@=evJ_q|VX0<;0ru@^2wQeVY>XY8o0FpLM_nXn_23M$rA|{04UT zDmFeII9Av_0fKv9gl@G-cxEX4BJA%?fK7K{*F_E(`Nkpbp?df)p}b<=bb5YrJ*ovX z+INJq83M$ZnBRei7~r`x zce){Dq@9O^gn*Xh`}6t9Nh~tQhpYX6U4di3p0X@WsjBVaS| zgI(qqDNaeA=l{M*>kq$EKzT*F5DeS<1SVZZ8k&mw`n6W~`lA~#c(PNr5nXQ}nhM;) zmy~Nj~9pI_+ok1JQEs&A`2-QDev%(4Cb3%HGEdDEl< z&o)~;P)u_0^Y=G(OtE!w`PlWN1ID`cakfl2S$UnO6)hpQWMwp{t&&SAX{o8NeyK10FxpTTF3H4;*ovduC=@U-6YzF)7(^l*nW@-3 zTMxnyKF#ZRK2NbDR0B9o&(c2A*UiY>DtX=^zZ04Av)}O3$S=uX69kCH3 ziRb3#HU!AH+v%V}10J1l!@E!NK1{ z@bK)frW-bNvBb+XG;sEfPu0UraVxxdlalIu99j^t<6Gf{Q6E3zEQK{MkQm-)G8kgcXbWa6 z+X}rfrU(4xQTsi&C_Vzk+qXArYfUf?9I;yF7A=dSN+pt04_Mlonwpr-r>$*~j%D9v zo(B^k_m`aK_Y91)=bznKHsM6<(Kb1e)7I5x$nq`3HiDR<|Y7eX_@cb}#nzE@XD)B;qnH{rzFKMi7`3T;$+=?hYPl;|7It z2|&60?i=J>)a+P~b%us;jIVI?XE4@e?+m^kp_wxd#IMOXbE)z$F~1|uoyMro-!S?~LP3YZZS6T>mF(%M_(keoB10)K-NSd*SoxtkzKNCggpH`AAT1t4dUu0=Ar z7jYeKGM2%6WTU!K?Fo)#ayY`h*9o`j*{#8u928d*DhRG&yxh!vLC|WQd`g1E&|(Tv9tEoj9pPQU{|G5sSiv`OQX?fVPVhU zr{GVm@sR$Y|6bAUVnBxqNG6kYbaa3XzncFEnX7oQdo_O!Idx4VjE`@nD)T+doKs%Y zR-N;BrcjQc@x$6!uJ(wCFn_qHdEzfLGo#59Xsh|!cvr3?LQyOay&7gQOAH-!IIDhJ;)bb?S)FVd7pz zN(PGxon7nDP~pHufk5(`;`1|K-FEoegmmgoPj{3jqzYY~ot;wYRm6C#qah)dNLY(7jw3F7YGCv zKRUaa=fM>ja-?DuOAwwN+?g&_WA9)wqlb_1l)`ddfKu#x_Ilmkq zX9NV?jr_BE>=yRoblX^Gr&qc@-__IG+iJeLx;o&~_{a)}$B1G~bsy~OvwSLOi^y}g zB=Jf|V*qn>5eMg@Qt`5;y}kYA%b2mU*RMH$qj9GCJ|8=Pt7I~N!F=A}WMP${u<$*m zlf&mQ3MmW)Lz+?^b=6ywwjTB6|$=NLH?7G(M?2sa&mG&z?1l9<<6&0g`Q=X%gQ!dStU>R@1Cfz zT)q5sThNW@$M*}q;P2RXjd_5@>Tf!lJSZ`2XuZ|cRN*%zpOF66U5t%yW<4mxsv{em zDJa11XM=?S%?&!bgyiIpvhreA^M?c1zr)qk)WwfnJb(TjZl#)Uto4ema>lWrS=rg9 z+oJtaQYB7D+3Z@vFj5^kX-S~1 z%^AmI7gvXnV84KXk%pv$80~djrQFZu*1nkEk#@G<&Lh^+HxHCB~WIvAH2SWSesi5Um<#uQf9pnHS6MVXuW=ti_!a0#{Effl#9EbTVyD>cT_9 zVpB1USU|a|8yE6p=FIbcGVNY|LIX`fi6EZ@d+_d`AM37IfpI?ubP|cw{wP{^yV1Wp o5uEaWwI_dgb@E>&lZpGcT#Y%B^cX#(ubavs}z?X}hpRge=$MZ!k{0031|LPQAwpalQ`${zt5`~-3@B>}!5+Dm9U0RUR> ze^01EyJ8anAO$2v1XbMAk5|mSl}&MgbEi1ggH^8yQEj9Gf#FqUcBEejKY#uHX8(g6 z>Gw|=s2x=S^cXVX7+$SgTbhEn(~di-spIWV5?g}Uw26AG#qI~c^`zzH{HD)CEq5Mc zPH2L$6r_RZLczfk6Y&1%LI@48aLnL)I7Ap+Aq293KqwTlz-YlgKrp({Us6Lr7=g?r zn%SQ$5S>K=w;Q}3j{?O2fh-W^!~g%<|8rAP!}SPcW)w&#t+spAmovu9X=SBSt2;hDcKo{D7xOhCON=iyH29y=m)zvjM4~~w84eFP^ znHv}weE znTULx!Og{$A&1f3Xt%~^x61P|Y-fKT0l3;5j|vN0YH`}@55vJ>(rdR~s11OB?|S`( zJP`d*hpQibbDaOJ-KhV!%h4>9fQQEeI1JmBRyOB-1pt@dmnUv$xzT2+-Q$Yukr_j1 z--F$Is^KR-Q~)CPUPuJ;zCUNh>>(CZ`0!3}P>^C19UUE1wtN=fZ}AXRg4-p#u7fKs z8X6j|0Rm+uB_t#yI27FMap{>k;VM*|9~6?LP`~O(P|$AcC-*|Z1M|zfvy~M717O0U zQO~ryg!ltFv#FypKs*ksp4H`gqGe}yCY{9iay=yvk4E@!&m}!Q{Ud``{>cP5 zKlpff2VOE1P^E7f2xQmoiyoiS8v98*U%5m9WE7Nkr@b+7%Gj;v6_k~QZFHLL$pGiG z)faFD>&!+y2SP+jU+?hx-admtWOVU>_cd}%{UP@g*P(`1j`BNZV1KdUD3XZV_hBa~ zH#Zkt?9^0z0s;bab93;XA~{JHVR3Wj{B~!hQ$q5t@*Q!vt?sU@I_Aq8sXWv!1JnLY zm3eF?d9}%-I#^sCkp@Qd(T8j5exaIBhTHYon&2Y+3c2Q zwLjn?0&|nJl|?N!UKQ#Dewcyilfs*aNpb^N&+cD(X9l)M#2NgyT;MY3{9TCE)S{bw zjmk0|Sz8qzY@KZx0FfP)@g)L=m?5a{KlA&OMZG%Z5iae28>J5>JY4#3plLl7<@#~# zJQ~!Wuqgx)$WZkVue?8xfJZDBo7umCphTuuW+B9HP;mkozz>WSg(u0>&Api}+`CI^ zh7NsnjSc*x<%`2Tb<3p^i0+BDuIbZ;hi_2jJvY3H1z0{9?L@Sw^4Q3)5ys?v@q_}x z1P!jiGl!32aAI~eq9^=(oPzR`w0nHwdA@Xt#fG_)yYi{j^MlG<{m4A5ZHtTMF+ zH-JzNS)XgKL3fVZ;}m-*{~x=f9p^)K!R|u~-g;i2=%|%|cC8%-WO%QqZ>dpdr`C33 zws(1=2xP(`P{#+4?)CcBZ9DpRm*41ji;SNGbm;3(9u}JYUlbo)o(mO4aICiMiq?a0 zgb*x(ppHuj@YcVvbSq}@Zgmkg9f^&t9G6b45Pn2Cyl;4iKsKFX2ps8rj5^ok`BU^? zBr=geD)rE=GD&erE~9k&YStmJcF6Rv>NmR3c$O@N1$QNfD$6r>>3UZ`jBtlF_ic10 z-iVM0LJnG@5`?Jb!Ly2lv7_^;Q-jiK#>q`jy`=YPAnX0&%|9--|I*oX49l=bT;tkU z#qQk2=NL}a4!`Ef(xNFOY6OakFh4ld|B|CiKGUnT!XY0*1i6JtVrK?2)V}Ff#B~|Z zt9LhGX0Yw}m)b5GvP0xUijZ@wIEp#ZNP#tT^mc9NOUvf3CWB-hV$6b`+-E4kE}nnZ z9>YK-2`AXj)@LdjSCu=>8+P z+E_u~d3CB#o-F2;~RF-RRc>v8ae@bLCu<*i@GOwG35+Tx8VSGBcCa8pkiBzFutsjs^wFmg+LKW>O>`7}*}4jA4zwV%%ijl= z`LY{F(EaDR%BX;ml_!DKXI#Ydm4)pVnAB7cKWB;!Z=o-U6al2J3(ke%g~$B6 zcQ$AGIN+&LOPApJu!7a<`&*)`&KpbHd{(3ESApY>G$TM=j!kX5lv*<<|6b35Jhb&` zOd;!)?|SNP3sPf!EC>du3WrHp%oIdcG1?U7guN8Uc(zg@Eb20?vcl3YqI7F%bcYVPaKsMo4Q70+$!}_Ta9h(k+-&1@Et(FSom^bHlJ)Ohkv(? zL~(U1mNbnBkmj|M?|qQw=OUC|;Y-t~by52p&UdP(pP{~{ ztbLt}HGN){pV*cAd0a1Tgk%s#Ku1sZY1@Ymj_1Pl*>K`D=q1FpVeZk}KY!Cx z^cvcKTIqD9yj1Z=WCa&W;J>M1_;USHUl4)c^-R2p zB@CM}kQsJ77BF?pu$IbrI>0Icxa)ZvxgoBaS0~$~=F$mB@&Vp1y0izOOfj@;s%KpS zn47jI+4xdVoy@ue7%W(2l<@uKgiaT6IUHTMIbJfNnoMp9(oK4$78ek>`}xVvNIu?= zexLGcQ9~9v!AM)4Ew)2g6wCzepBMJ{_N6PUs3N7QV&Si?+$U(l+pZbEvou5Q-P;u0 zW0sbpqaQ&pg~83QTuo6?0eHZCYE9GFA$D=&YzBEE1oJltvPuTT6A=3y7eiChT#%g% z8qJPJzY6P=heLI*re7VUuQJG`f%Ouvmy;X^j>_S2HFTw?$vJ zug{8j45281wakJa<34xypxOnj(s&a6JWfS>;zITYVX_(dOye@@pF-xZlz$Cn{vj-m zVFu_@!nIMK2D{X>G&+mqNmv&0UfAtLlxH!T{@tbkWDr~zIGds(10HsrV=rC!p{ieQ z60bgpzxyUqiK{z`35HnIQekKz37oHJNAtNqG!uA!Xf$4l-Ih+t)@>OFq`QtI=lRSB z=rHiwjK-J6!60fg7Q8%bX*;HI6SPB5h2K(s7xsc8nkb8q>HyY_ln(rA`nmLI z8&pgUKOrhG{c{r<9wW@d5gI4>>0cSSOChP}ecg;buoHio{EUB-;X(e${9b+jvO7w! z8z@Zaijxu@{_-O=&U=xPyy2&74N0a-uCbrKSAw9#Pg2hoV*=(tJ!esu>zHnGS>u4P zV6Ul0Z|3r&DXv0uWH344B%%CkZ-L$mb5^8V z*t9Mwo85?-3K;zs74p|Sxb+;Bu8z=FV#^>pOz~?Yl3f!X8JV|wTDpl^%&Nu()T{p1dR1D#OoGnZ z$l?8j+r8Kmzm|*LkBsQwekqm#0`uj2$7x8Maq>7hAT%+qGrI7k+(d>QXx_Y3ZY&EJEpDN~IVDq|tv92}9{Bu>nAx^x@Uo z(c4g$n-~TJ$?X4)%zz}ph|uP5%N$urbEjglcaS*a3zeWxQ}|`@Z@7d@!kE-eLLmfR z4lWmf+vm~Qv;azlSCh9Qtw3^0Ki0TUM)zwDk^Zy04u%bKMH6!WYPep&#H%QgssKU< zDfk^Pu%yJq#N>XY){530N^5^XQ~#$Dr$-TU(vmQ>h|`jydC+zFbNu&J4!_QBsOm== zFO7*$aRYB(ISAs0cy+PeI*vyusj=Y&5pwqAizGI`d0cLXVp6g4)%(6YcD+6#18BrN z_cBqEc3)k^kP2vBij0sgv74$2RJ||kj4}kQ2h4oF()_agC*#CacRAWIF@0j$`R%@K zxaEk7hS&xrBhk zdG%LrVs%HJN?ZVe)R99*<#@mdz3$0#+)znUneAJr#N*>>@pZ!q z=kPiC5f+s7W?yJcO%2Egn6jii-JRXt*=e1Np!dDK-gQmni61?1O-xLzd0#~|io-$+ z`uKD-iZ2aFR1AIzlcdMr6<*`fxuTqQr;WU%J(#L- zrrzwxQR<B?>L(J<* zuivS&l;?J~+E8CVefxcBKDAc2LPqM4WRB#?)iFBPA@tGG2wH$V&d>WB`%cKk?}%44 zN*Yt1ot4hCN#P~=d^VZrxMl;qoJcaNA*LuQVn2;{K;*~tEb}|!ZL6YAf-ph;QBx}< zm##I4rXN{CYsZt?SiyMk1GmL#UmzBjGyns9`jm%Cz%G-_YzW?&`B?hH)i}TQ6H;Me zq47ZY#4UB)5Jjet;q8<(90S)mNbdB^|9aR!=z4(ZzIzT-0eU0L(}cQx~6AKxptAZo+(UGfypC?JY9p zYq#;F3*|peb6#9% z1f`a9_z7o6H>xtQ%nITRVWS)AT@h2aFbGs=lb}~B!B!ZfdSCLb@QZbvVDpJ*vFGIMf=6ja&BkD-ev(LCVGqSPq0g1Cb@u<}< zzqQ`pYUK|?%Qm*YUs!01JjXJhe&e!p&a3*2_voPobP#HJ_wcVG`0Od50025*a{7n} zRKMx1x{bo}HtCPHw>mdBre=bb`=(jNhq4h)dU)i=TpUPhl2^Hh=21HURO9xM>bH3< zTyMf42DyN>xi?qg4p%W#;TUc=b)ES*nEGACNkpG;*NxqF^LoHytQmTUu4DNlN>0<> z#M(ZwTnoXzCPA0(xKd}^neD1pM%Kr`e{o(FReKxzKhq(FtOuX-+G=NNTd$ie%>KrH zs&Rn5MzSv{K(*oT+tb}cH5gAJwpD{)JAGr-)1;-c5Y0yu#~}=d6v2)8d2^ba0;pl* z=tA;;|0KKE4*mOz_!(PggOvZNno+*Kb$Df98(cUm{QfE54sDI6=8Ehnn8Ht%GbMBm zq(0TSfn9YF{~R3}j*GM@d?YIJz{Mb%G|slF!NI(?Ms`m&go1||@09e+KvS!MS6j>I z&B~^Z6eCSEr4V(4U|HIr)3uN=3aJ7Cgqc%DifuK=4kE24Uja%gW$6~j4?lJ2Z!8#G zSvv2F&Ygb>G7{k2J2H)~=%KUxhh<+JsMKjc4Z?9&WKL%E_`7VBc-ihkOG2Zs$sl=~ z7wh$ytprW4G&us@D+4A^M;#VHoVf2MMi+96<^Ih~crRIXSj*HDm#E=(qb$*J+m(qh zDV)i|?Glb+4lho#{@y$3$kO=(3xGu%4DwS6l^TtY~=c^_fL`*j~*I=iFr_db2IPJf8PvUBEq$uGI*v4c zV(suDf#(cQ--rxdB)p`NwSW+OCzt7#?Iwk&GCZn3Fn{R%b&9d#+P#mDUQw-UWg zBlDL_H`b_UskOrxp!&sfmr4Dr%kjG5(x2Mzl!Dao4A1RkFywJ@Z+!HvX{x5imh}eF zO_A|kuX^klG&{fBY}qM|TbIt-bh9rlN68qDi>}t>QGA_}iyn`7{#GciQBz5-RWZSF zrM>ZRn{4&7_TTDd&Tb5u3Nj|bc(Y=a2)VIi_imFG1acy%y>GuI&RB4BZy`!{t_qa%ttb!LVkPEVSMWg?cE-;UBD|T3pexzQ z4pm{4XWZy?W4^r(2>^F8Y?~%*i^c%X^ou>a>FsH;!)VfX2|T?YcuqR<_l(MGT!aIo>YYrULg&jlX)0%F|xI$r_uHjDLBpJo*!d z2k6V5+E%cahAQk)#0=f7j4CfzY=*vu(&+Hwd==92A;OE68CO7*2>|QynHHu6W0!YnL(NS_`5i46+J61}zx1NA zro6Yz)IOi&cxlcg#SUX5srdp#XS(cBmt%%N{<0ETeag^-xT+V-CN6P0s9Qo^(&$vrNTEhqr-S5FKO4v`}nB{mIkWsXiTda|BcJ%tTewUZdw!%?K|= zKQ9)3hBTHqlF|QrbEaC2mmco!Jvof2sS(zuVTQtev**;%ypTCv9+N&2*04#_SF1m* z8k#BJb%Srj-WUgsf7jBbd68iRdb%Oh$naSXSM;3?Ou(%+Y%t@ss79>bnsxB*)8+Lm zwUB#4b(c;J;dX;V7VU3<{;&EUL0aw4uRZY5(%VbZAN9Mr`-@Z=4Intz@H}i77=O&8 z@nj>LOCR_4`aW0B!ie30MebBO1EcXqQWekf3s@IYI%uTBYwE9f376ODH#j=8}2PEuB<8jvFeiQ$XUSC7Q46iD<~i)=h4))YC_wVmewd*;Iut+ zg*2KN%u9{$zG;hjbE#}OF@$UG@e@0em6R7Mvi_ap`#V!Rc!%NrBRnx_7+V79h6Jki z+HG8-0U|69f3s}qiEye-;`Ni0WM%Xi8dlPK9h=x9lz%tu=#Z(C*05R(_w_iz#C^*o z*9sK^ioc^x?0-^O zt+0%sK2%QM)PDV6FTjSWIXmsY3#}_cRZs#-nij*hHJ|iiE$7T~WTv$nqVB2@yv`V> zh~IY##DYhbb5k!%yb~_QY);j>xP2`o$7y1r1=BuBI zah9|Vhb!F4NA<;NQ4UNF?jjSDeKM39QEK8pA!e0PNzcF1cZ zi-U5>qJHEvH$3&|TPcOIg+vi5dUHnXem~pIxF)MC`IAapE10J!Yw#hz*7HYSKg4}2 zJo96+(=U|18r%fzPvRt4N>ImTp-}wkGpe@1L;eTcr;1HY`a{1*N?+YBNP%y7%J&Nw z6#X6r;S@MT?!|m+j`tCT!A&o-dOJu30BI~FF-zw7!1ihlKH$_;75@v~)9G$P$rx;L z)CkD=Q}N+qB`E6-CtY7Px=c@c9S8kfB-~Z44fp-L`4r4 z?^M`|M|63=XC^&9yOZ1K#~kKjeefs!b-KC>FyG8^qY5SPME2z_g0m6q$HLkYC>^F( zb$4~$I9HvYIz$Thv|ziH?bMQP2vp&7+O%qp)qAT$vhHO5GXq``i?uJMTf{HNKs2)wv=`~`}ON>XEJe6`*^567BTRNu6V8E{(<--cMDymII4Ju~*|>Fm3ZGFF4xu*IDnfD+^NZOg;oAAhyMpk@=LA z#aWoOK+Ihye&EV$Pi=DjXv~vh*pgG`rE{reJ>8vbeHCW{iG8Aya5`^uKN-^i1mLoi zpqC^8x_djbvTNkD-vs@r2!8yPR|yUyw}c$L57QwfXRj zm}wClsA^kJ`F}<;+t~{*JarI{{+pi-rLYGxgQ$L82^HQo8iP1g*mLl0hQ5PNbT!4V z9Rcq_r3zqPUNgLVW{$i~Z!OxE$%QaK^^2SG#Hje7Y5WjG&x*Ktvo7x%5Eovquwx=g z(Soq;RX=CY7f0a4ruqhM^2(!E8wrsD*2%L+COS2jfk|~O7Z0C%%}!gz+vY*lE#S82 z_KB12tJ`Yb7h@;b4=L%_wO^N>E#+U-!#dZQBRW5=$W)c9NH_2&X$U0@?4Q}V-Dg07 z^qbm%Sno3tv3JErYOwE@eD-2Pqx?(#<-sY+gQp}u!hiTy;6^^&)mJCr2yY5 z%7?4)sL4D&OfRk%I(@7kM^bXL?2GSN_;ub%*YTQ&C9yVjj7>-ShK1Fm<5?nI0R8gY zx~Ju5M2a}lKB*)%e*K|jtmDahIy(;F`6d#02ysbH?;Gn9Go0F7>@_t)UFL_b9 zt*!7daZeGR+#KIbY;t|fN2%)Oij6C{ScUZwhx@ijlfmh^kT!Zyl@k#yR&aAYNjzwb z6sXjqLa9vWNFvKnF!?n#^7}zZ4%}#H&lYCNj=c=HKqpDLD%b}=32LL}r z8n4W5sl>N@J7YChs6EdNP}pWc6Tzp6ApO#Uewc2$j+MwR{{7`xGB)pv#9=&;>%4(R zy!n#tgS$=f5C}^4g*T2i9#qr+4N9XrZKkv75n4dVaAl>OmfLA89$o!KJ)CA-F2?cz zyUsUiOp8O%A*r8&UN$^lp4er?osQ72xXchYsjBr{^yVvrvNWwli~{*qb#iled7M{E zlXiC4PQ*2mG-e(RuP*nJ=WRMf46Hr%Kd0b-F?GBg<)Z;U-Xab%l-Ur_e*aN}j4hwjqfPJQRzmz|-~+lA}l% zmi*wv%Qd%jD*AELkZ99zWzk#~(c$G#lwvr5X?VVk>*6 zB$2&D*bvok;i3;D>URle+j_nE@K#jc_Hxac{ww>!Mlo^ZPmO3vZ`}bakK;?pz3^9gXa*W&qxAkXGr3Uq_6K`2#hF)!6JX9;#3Z zUfLHYH;M#^>6oGJvvYAsL~u&D{J_1Recwj1Zwy8q^Uoldidl=PAY`Zyrrqp4GC zPW^oXaDqBgjN0D&WDUY9ayMs9wMX|K^x^9g;(o$CE-EhFEkhK2i6SwaX7aQ7&vyT* zsajHeN|WnS%kXxt<{Umyv94*LSsGpZmY1H1Xk_*uA$cf#gin1kHMQdC_-pC@bS#r= z4-Jf);E7jgDs#_${}zd=$#;P_SKlvwQ4P)!r@)<$UABUuV?%uo4%F-*^2Zl zPTYInoymabvPJKPay{z!zp>KIBqR~a6bb)ob45NZK$ImPMVJ5s;Am#$4&rXKcx^+Z zU;axE(j!Hk@K%#eaw6*4tjVc)SSk{y(WmRi=Q3Y--rUzQ$*z?e9IAsw&MOF<5yM+462> zSUKHlV0-4(gG0SVWvd{*wjbWOM0NbF|KrQGWfh*!GT?9A$GTM=pZ3a2Tj{5w^LCXm z;$wG1gtuONQmq*BNm9f4eB&kE68J?|+?J*7a~RNLYq=8x2vEn2&U?sDBwV~~bPgxP z3W9v07iU3ws}E&obE{%O)Of`Rg!6@0)Y7VpIsr9zcxVTdi=g~{R~sS1K^5em3f?7m!ws}vx1Tx zIg`k+VBoEaNMiYpqsFf*In}(^OiHzbRsN5W-Az5bcH92ZW2s?dj}~MlRj&HddZ*<2 zvFTFz{^o1t;Q$DquUXgA4roIo%j7=OS-M73l=IGExoM&p$?=cIK`m98O*TGbDt|j^ zYQt<=TlHmAo_yrNnC?WfdQiE}Rbl@{9y`o>>cLo}(o%FjmbTC^YS;J^zY|0`0*V#> zi9^PwM(=>FhC9+#N6x%_IRGGuxc2&^*8e!}+M3UjLYz|Y)w=c2{#MccQSHlgWC03-C<`vzSI|}7CY2^S!l{Wd#g0XoHLgpq z;B@hpWIJ;}XQ)utq_wB5)QCAvSNenf%M@aMvA>{P)k_}z-L|S>gZu@hyje{^S9rR_ zg%0h+$Vu5EWsm@-73cb%>9jxcc+&rvu^irSl8hCOkWBXvkCiLZH#9NRR(d?2 z43HwcCEz6JGqedUjxRN)-nk=Ck4P7U;#chKmfU|pF<=FWDbTGuIZWV9G?EmF?7rCY z6C|o$56S?6ewKTJ9U=$%NIV`C)3Ik!xH+KG2` ziku*T;VF4+ya;OWTbpZbt`25;JrPshxi4Cw)A5->0)Ewt*!~WG*QsF%F{9yX0SoVp z8?Ss7Il1GF+T#M?#xwR43?Zl{f;aAbrzD_S-7bDS5&AR{TQ;S0!cW2Sdlo8jiBvh(xV^0hUL=EnN~giZQKOQEUy=l!)DZ{du^8!!TLrp8K& zFa7c%d^0H@C|awSV`NZja))AGe|`$;!LmJr*F6|56blTpm99LY!&Q&bmHaB>=2e&bGj_DlwtJ(v z+WcZU+qg>}=w-K;E=~QqS>S&%S1a-%1tZSj3m!e(1ZjRq1c?4a`&|~l@=SNSr@FY+ z9TqftpIkZUN*^x6*)&(qEA?Z7R%h0X*2zYNA@xErnNATkB`Wu$s{N!^?i0^X;i;G@ zlT$arlkJNjIII-79L!j*uMGyMu9`ybHPR*zg&VuX`){Uv+oVt0yrJ1;9;7bB7XeKdj7*QIbJ2(dk^-U`m1-$ zF7l2Mv2h289UkNJETZ+dE-44Jvv7a_<=MhEJM@RMO)5dbiQ1yX4a~^z7|f@a-sc)uG!rc zsiZ+2Zp^|AB-F;j6!p%jw4~%`QIV31uc%yX0#8nWN6~komBFcj-r4%<2lA1HN%v?w z?N~(`(Afr|Ie33%MY27MGJXAXMTb7*iFdIpAF$o_xvJLZ|SqH@WK4qL~Mkf}4{Sk3T<+ zD2$aw9wU!ri*i7s>glq;Ka#t63=Kx3$kqCy8`e*Z<5xb)PX)PFa$uf*ms_}l>RLe( z#|%BNEQFP&P;>>ZUj@dV~JJQ4iD=J>5Uz> zW-2+RJ)QFD@g?g|-4AVOa3hZyo;=_~d^<8tr)C|DMKxE8FotjgE7_=6 z69>eAXtmk@EY;Opp0o9Ol~8i>9VSUX{s6NjHy8DiIxinwxEr%?uZ%BHg%$+#m-P3p zCvM}XNYIF@UrUYrhvz2ftF)}RomMnZf>qq;1U`UT%v~A|WfKOXnf2^WeC;resSF(1 z)Mcn&vOYX$)O5j#l{!juy4XTGJiJdq`$)LcnQf@`2nCU3=}H|ZMF2@JD*BJ@Ck=WR zunj5@!my!v)()qc%Zo;irJI4Rcoj}b5K9v$XrCBb)p|<3_qtN+t6R;ywr)~5d9V|c zflD|>1E!vuK|Od6i6Rf0NS-SGT!eORw(bwl7-jLba$!e+(dAL;jM7J!V%6b4b_qq! zL}rG6D))}KUDQ`e-vK)&a#C6rFmTK^&alis1gb5j3d_sOQ@)zZrLjjWHXHQ(1kFZA znuPdDGlNsR12~^tj=?&{QQ6&7vBKb2mJbeh-Otnk$)3up;FkBR_7xh+D@rRP;hGT; z-2}u^>B8XsT)*o{_sn4f^n=kOmwHhq)0n1Yxs$lD5HwStmY4jN z$CkgoeZ?(3Gy0B`#=nY^_jxyA78z_tpu(PUtNHrmUG5iM_POCj(cz`*eOP(;LYTOF zT5TF1lNgE-*vBDFvok;s!t6FqEsKoT1C#;VmK5C?1OPNcBq>AAyg7f*dWgW|&~2^l z-ag!3=UHaXkEOr7OrKWH;!9ZR21$OzSW?$k2(mXSEMJT!&F=~88L@<0q33f&!z?W=%|?=Iz_uphC=Rz% z9pL_~>n&3Wh#NA2g@t8gWE>cv)fW)TsqXAt9U78c!o&dnbLB1($gXszC?bk{j*iw) zIWxA~J^QY!M3qer6W6m|jrY9!SHwq-6%MuYkmf25jU5{3mhc-S(LZq5_uU|Asw!>os|o zb_!!7?LPIYH?u3@-1+zw7e%P@K%}Cnhz)qg(a!fa&^;aj3$v25aO&3B21eahXV7XB zgp3m$5>lnt(ONF!GVeeE7Uiv$_!y}KTA z`T%4XD)Jc5bsXUpOb-SD z5cYu$c(YadZ$ngBD1_h6!OlVLCfkbg@@LR8mNhW$saJC&K}hujMc1pLt0^Mg``K_R z`Cc6TT>z;`sW8bRpX8F>QNx9hRv4X3)F4hdPhLy5Zy2cYB+Vj%QKGF)%a?D=?{#$_ zE(c;LzN=#8#n`_9(eOZ}QPq$af4-Qd@TeEDS4~D=!CYB?urt>PQXgYg2 zje+*_R*opiMnTulOl`5mA^821;0}Mupl75=P*Qy-s&oNTHMSu96((USC*GB3V!Ixp z{j@ebd!98pIJG!Fp+iZvcM=h4`KnKs^ga3X5XIQP;?Hii$mrI`MKnE1sxGb+=j)Xp z1fP^Pky0(CMp0bWc)kWAS?b*uLhMwCf(p_1(G0-$Mk3bl|KLfytHGvFj#c+h~$^8G0W|-YR zQEr{TKqqP?NDzfZQ-gj=F&fyGKU$p{kIfLp=&*vQfRx;=-~d?j=B|J*X@-3!)N@D| z)oif-m2w{=!JL&z1T#Ge@oSB*oCh5x;0|;(AnP~?Zaq4{2HBbN_1|D^?#)}1misv* zGHiG$MfI*$WwUQ)8i`Y~Fo~@iS41>Gg8b9qojsbLikMFU+I55Ww|ti!eNgjVsx;}t zo^H{=W5SzhC~E8cq}7#TrSLC)sMYW7 zRk=tG7QV^vkE(;iMKqZ(mfPM~1_3eg>p@XgBmqYxJ4xZvXc`9&qfYUOB43ui^3(t& zat4$>eMSd2Ehs|9{)c4F+~9HPbX_>L?9ssd1(NKDpKl(xv`G#$tY8I5s3vqjDbgoI zyzRIuqIRL`6I5_@$$Q!fIY_Rz&$+ zQJIWB-9JV#zXFl%!a{+K4HShpGP<;aS_iXg@AG0$r)%SGx~Oo6c6s}&r{B2o``?G} zqTWaL_(7h}RlWn%bWFTl@#I#{x4rEnMk0!8gF zaZ@r|H4t+u^I&aTdK}al&P~z!`-9CNOzutbA17cmjEsyxTOzg_XeI_9Q0oH`1fRsh zllgoeuMeh{UPyob{8_2j!R4?i(x1qn-S`e3o=Lm$cuwDM=FKG-jmUEcf^E}(;dy&( zVO7T2vq(C!M3A8|?Tm57+~qp)3JHpV7Q{60%64} z9a)<{N64ag|6fKqtqqMmND&NzZ9M7__DURy_(8-OE2n>r&NW|G9^g9gFTHww#NoL} zLWi0-jJw$W4mP($>;WWgPSU~0%`u|7W{96zj6a#>NwjW>-GbyAm9#OARh)ByBwHzZ zg>6!NYA7I}Bd9Bx;s5C4)X>L|$?jz9Cyy3(kv0h-i!g2(B>8BazVNj&Ta7%rs$*+j$=e>KAs;UiO$a;CYN}piY znq2wqeZugMs+SR6d8L;*!^U6?^w?Q5XDI)%$&M|*&JJU_tc#v8^AA#LS(23&SZzzm z*2xo2zp{9Na1>|fMpiZNmJoN1r=zCEf$<^hnS9bku}+QK!%7KDy+axAp@?fpfZg5G zwMM(6IGwuHRxB48(7cfB4Q3E05l<7>h7MTs3;QfS%dj%Qza)`7Xa0k7w#T+VQ{s*q zn=jw&+UUTxh_#}SVn?In2*@@3Uq6Pr{7I*P%}o>g;z#yq3e6rLo?oeI5TC86;yMyk zh6)_8i>v>Mq)?p|5CF9|kY`Z$wAWLUpTr}Nmzt{7tMowo_h(7RwN9TOKYldXu5kIi zK2LG~`YfGTu35KR)%A+cVXLF4I90OH7m6tlCgL4PLIeQ^Uk8beKyvAW#NnFTos-pf zWrai4y{sRJTOXie8#AS4{!epn85CC=g$bg8;1UQf2^O5-?(XjHPLSXdf;+(-g1bA7 zy9RBX;7$i=T;_heTk~sccK=M(RCUosb@!F`p7;3k9JFtyY@DT)qq0{61DMJ~DSm!b zadOWL6?lGr9C3QwR?gE8PaCfQ=<{^f_Ik;r)*tn;;tdw+`57d}Te;(;wfK^qIZm{9TX%2%OYhS8hn}7Cu(O!Bm%xQyBBxVVr@)fgs)1NkAD! z;&4ggK*L)59~qu~tPxE9`luH)!(vsB&k_imA>=+h&dxW}#JuIj2HeEEFfSeTMJZSB zs(XYicAFEkoK=7t-ETL0tsFbM0t(tO;~I+-LQrzE?@jP`IZSN1{e6L{#_5waF}E3nKr*Vhx%(?%^0E8{5)iVAw% z?wfiH+BMGr!9#RZ)cQJ^vXw%v0I#s{+fucDMh4MI0mw$`9gmTWk3cJ10g}I|0RW-e zAeZ;lpXVx4Y2_kR08b;NH)fa-mT$tTl4n`H>jAMyjcA0&Akfbn7KG6%(N&dxy#&&a zhd0;H4|2n?H*z&1kNfv9wDP-r-S<7=LqLV#M#a-y@GEa4y;tQ7)kvAu_D?+z37BO! z{mqy(M-$|&hIa_Nvp&?ap zMG=ONnXwXubfCD6%QSR}!q_>d+n3|P)IpS{Gu>>jxD3Ca{;AK1Vp=Fy)tth8BglOO z2oVm@&vs1uc=I&IP1%f2PiM8|w+V_c(Pc?UFh-V^Qu|CaG&BGXjgYV~q4Oh2)5aAc zA)%XH^LnXRB2K%7GHVx@ob}q;yN6#nF4WC%GOlV7r|~TF{O~JTjs7H)9($;CX^?s2 zwQdTVEyjfnpm#}DcI96z_u9W;`;vFkR^S%X-$v_QC9JI#-} zYmeBJx<2E}SOAhOL6w$4z;<|hx4xBKV8IY~T*H zn?I@QIpGs8f9Ef5U(XNCXn#d4x|Q90WfAR=>s6^fiww+}?Y~(<7*Bc7+RW*~2+u4H9SUO)mkq-y7CIHCdu4>Qw<)sd1^G~QN(zcn1xS{7)$ zLL*E1YCg{*0`2_VP4?=tPjf!olZJ__sdQKfs`GRgt*>&<%31ZIO5E)9`d!@`$%{xq zK+Xax0U20N@xw^p5rdXYJomA~Bo*qf3rDKU>Eveb#(S4C$C_uRh`qApl>NG2Hc zpX3106d=@SZft}P0yLt}a10@&`^gDle<3F3mt5m668GI-k+Hb)>BF7MvwfBUUk)V& zb}Pb&wz$$AH+MS$1PWEpImL(C?~+8j+99IeIe&)NN5?~Q6DRR(UT~iOxm1PLV=L1^ zo7mC6H_qXMfN=2aye)VYQ~Gfwz2PnLOVK%?G2cU*yc3Mgou_rbjw|z?=V}d^zYJ*@ zOE7A#qC9Af-*x{cSw1#`Zi&!DXE5o!4h@3dKwcn(<$ip7pFr^OX zr`w><;|XK36^oR$QhG3OFW@m(L^b*@YP`kyE8&#cgOE6hLomnO^=BXH>b0K&KmIF4^?NfB> z8831y&ZCu1KueRaM%8xOZx48{tFnQ#GK~$tzVbW|R&)W*j)7Tt0wy0-@en(u$l1w$?8=8GeHY#0)gw zPpf@^W?(ZM9Kwf?W>$yT*b0PmWh4FQ*Ce2tdE-O{k zqE?o^ZO7;_npkyXu=c||Y_7dwHE*1zO<@N5w7>#lz;?hzD9FCVUK%PXQO(vH*Q}E5 z$c=etb%_ut?fw|r$7dpPii|8Dx_FK zO$|w%%DMQt%&=%2QafWKn==xBC!}B}C5+0%82J#*Pkq~+oL-cg#xjLrrwVPpv6X?I zIa=n$oIg$*W7<4vUrCG&8A!4le8w+ylszZi9~)_&jLC=<14(J&36afv=xB#K&j_di z3?j?=E&yl=BP+vbe33w8pXy7(C_o8Pn(C^paCW=BRo_IKT$uI~u{4GO#n3mSC4Lx) z`%F}{pKN$$9c?je-Bu|NK4NlEP^{ zkhXQ=QiSrB_lRizl%Y>5O|ywJ8dsYj+I(i$yA#qx-}xjcu)(Gd=n|@G1RQc`!%X3C>QnY z<{|o67{$)g%$=|M$lWh;d|uYe(bM8j1v32uD@Z|nEPm+|EYOA(dtTCBSQY0mrc^of z7D4Ohe0OY?)M~opsoqik=DeCg2~7b$wm-;pJ!s_M~Mj@Foneprd1s|{k=DUHBIL_+C~f`=l0gQeG# z2vZ9~+eSo57Sq?0yW8xPW$zYNo7Mg^ohLweyt)5mKmaZY7G^v_g=8&*t=V0(5fU2a^FOjbGSp|A_D6FrwoFZ_ zNb*?@zB=J?TOogpqh#j?$z0N4!vir`45h0r-?|nXQk<{F%QZFyZ(hK|- z7Qn|3mD8pn7jzCFMJc+bM)S{u-P~2Dq|&YkP-Zw+KcRs8nZk%#I&V{^SRx@KDL&fN z`YG@2i(;oYl?fZpXG499@yUiUvNAPDbB>xL3gThw;?cN zE};7SeC2%9Z)F)Ds^%99r|&9+Zzv`e)1gaEo=Ir&KY#mZoGQI_r*%YCIj0n#6pN6U zz!temW+@bGI$ZtW2$=@a*pAqbwfe>x*;%|o74=(;;)mHdCJmnq{kWV4tZbxfsm;S) za;o`ll9ZB=b!LU&fi}Go-#mhdzRW9Ig@ny zk5hQI;byGNJxs2W&R2T(r>J3Ij^n<$`hhlutJi33B~?%Pz>HKm$PXcATw~U`4YbhX zQ0&7(QRmd)Qo(N^6=NhavC)a#(wE|4uQf&FqwD2QA~1KoVZm;^nkW0`ix0Zd9f3F;~Rm8G+T6>?yU2lukk>a3(XmqTZpCeOM_2_lQr$y z8;)V?edYrN+$~T?nU8nPtk#;|{ z-t@0(ebk@MiljV&9#>{r?rsUs%QgE#@xkD6sR9+IHQ(aES^K{!SC0rit#Ov0EV_w9F@rL}*e_jCqVzhJxf ztkWs&Q~W&_mz(r z6#*QdtE%_erbjOPoMOx4tX9;eugEdPBc=w)<6MMYG_7%lelSQP)qI$>ZQ8!*uf>;+ zhT?m`AHc#+>aEjN8AH#TPprO(jSi){9c2(i1{rFs`FsLZf$zA|za^FiEZC%0BTSJC z%^WVJx3GX$YhYqZ#UeyUemP^|u9SUA9$yIVkmY~v{%024K@6WxnC)(Vc{>|(f)pK| zveXLt<;P#yRo+&k*mHpMEhB`nK2w61ru;Ze$HO)eVz%17UU1`fqggxuUF5FmI{cpAbGe3$7bUr}s;+C`~x7YssyKs%myj5$u-GuD9c1)hb;Ez>C zf>8CKnqw=)J%y~-`WB`;A!oTH5E9=pTd5t|Tn+5goVhSnqn((Bl7m}}kUzi=pY8fB zx9Y&wNbE;BNr~UHCtq(A4!t%Ck1dh%h}enrtbU@wE&urdVflttM^Dbt|4edE$;_pY z%|*S8O?h)!o~>4Gm6|@%wEuM>9`~~Wg2V6?NV|joYL0-#T7Ri?geki`VlM4%DJ)8O z302tBQt`>A?_XlBS!B@u+xWM0f4j`$%%f1O)xctvr^TP(D2gC-uJlm>4q1dM=hXO2 zlD(tT-aO*t4(A~{noX$w7W3lRY?|gK6{h7@xBz!-q zOSML<$twA&KHWVh>T*c@*mImHorkGcm6qK(y;9{`r(6yDf}vIBsf7sU2nkqFt-9zb zuzaC9*C51z)gUL_Xu!~uT96wXsq|Fi>w8!X*tW;Ho zETzU6QbBB#*8YAyBU7pDo=#qATYWTRyVWe*3XK@pQAfwI+qp;%_dh2=Hv|r z1Z8P-HtYr5GeQp1T>6wW@+9wvR-8djPO<}^Yyr}3(Zq_*_AL)fXwzlO_%pbtVzuezuTnfDN2Xs=YKjC>o3-^ zQzd9W&>Mcdhjxc8CdSAg#r5Y$?CR3^X*31mdIjh_&vM<8@)fQt!SI=Z-q!8Qe*fkj zycGg0n!A0`Mals*n6Dp$#ESei;F$yg1OX^j%`LgTz*j|f)M&*veE)FmEAmPjRQWtZ zPJwL!H8x{EI;DBS=GtU-WSymjpAUk0*Fn7~-*&fJg$0G@7rMTdRAo`-{+G8!`a+oB z>s1t0@*>sG+VwYjg4QW5p1?|ene3SP`Hm`(Exz3Fz|^dUcjBE7;$5#uYY77X#e4@k z{-J%(r|IiRHnL)0M?v`O19PrfjPkH8Gsd;qxLDndb zGzQCTn}fwS_f{wk0$#5?8~;!xrgwZ=p5G~le|lNo3>>`lh&TJgUVfaWtg4f#Q-AcG zorJ+~lL{>td}=M#M`V+H+Uvp%Ujp>5(;E=xYlr>xRlaCjS*cugh14t86c*N@y|QmY z%+!|ZUx>}wxxR^&3nto4l3ZvvU*Ohrp8U~LBiBNdk>HyMeZy<;JFo&Ql|%;&P10CJ zW##Ef1pFq;p_irwOY*681H>(yGle;-MhxXz6QoIa`Fm5x-$4^ZSD6?QP}*IzBDg z%O#jw=@y2k6C+QtoKQ%@tks!{`}sWsRJuGQ7~91;{GQCZ4Hmn@v9!v$90w{;C{)}8 zhMQW0OIDaaX>GG0;ySS^@2oNdI%$uhhfAwd$NvcQ@cACjlPxW2{~gqYCIR$rBVk9!M)!?Z;t}UG4>}Fxf{rDOY^CeOfA)GEdLS+bA9LT<%IEjXlO^9* zDDak>MoHdrJkI)2x;>#mTOHu=Nht@dw~d+O0OyD;gW$Ki6w^NjID?$EyG7 zy2`tp58Ai++C|7Fva+($aU>7_1~9d249Y>I19iYs0F?xVJfRIWOZayWa`!)!f~&4* zrc(22(~=iY*9mIve&K|Dj)$gXw;673uub8DWQ?S_tZ;-Y`iE#ycFO$ zH}wN8$^dqHfa`_1}n8=XraA^77`N z&A!AeIj{7A6T zb*kzglR4cjpN3gEXy-ioMmM$fOD{y;U6ck8!-W{UC$)w_u2jAMs#5->BMA>%~k3@=MN_qaM3pu%q=%Y zYWI_@hy0eY-8kJ56X|zD)Gy{03}4yubk`JO6G?v-Y&G*{WYzst<>EtMam}b!YrmBe zLVep4rTrnLA3n|c`+J62pNuCG@lyTG6<5Y{ygr}*dkEPK$)Po-?2#-g906;r4CY#1 zUYx41kXN(I>1@^Te}30ed%@~Ntxe2b8t|e3p@f;4ihD zKxRz;SWCANEKd(^DWK{8w}N(f#d0ZHYI&+(2erK+&Fn6;tq1b}Kw{^=838gab68@q zN^r2mrSmcLR`Q(&fxFh$5Nbb{+p`V&^H_V@$Lg4Coa^*Y%x;p|;78ejRiKUV{d*+a zKNBmkqwsnry@TINEVinTNQ{g(|J+XAfg!@9u`^6cMxXK$7->*A9yFIKXN zOHASAlyR4^zp}>Hku6{y@+ZK1Ew#`RXou+aV5)5}`fsN- zp|a2zOw^vLl>HWz^uvEWsrn0uEb+Q<|CH7z6HQK8PJJogHR#iPpUF>{cM1xff$u1j zjjNN??I1&4j&1|U&f8y)SrpK{OznCL7IpxTwRJKVXVcg(=VYoaI@4Q*Ml(-54r+QQ z{}40GXMGb59y}nlD#Bl|O6G)iGXNCW(!^Slg&Kjrnn^bWq=%?>`dDIdTxXyOsdt_3 z{B&h@92aVn8U~cYpr<{ai%sQMwI{+Wrg(KkTFAf(iWRb(!Pgra4Bsy zyLd%ZR!j8uldsc0-F*lQQHfT_HQY>o$tPcv9i6+#oTj&zkP14Ie2>)|6zLt&PWvYK zYGDdaf20RqXK{}F!B!vi8@Q`^o)q2o zLo*|D{ShE#)`R1hB-LhX^l0vAV40LWQFN^fdIL?xx$ycm-b`*xr?4Cmlz6eCB)Ztn zSp8XjTI@r9w`elR#V^Fp)*mWbDoAwxYSZ+T(zV(eZ}lmJEY5|EXEZo;{y?Imf?Sw; z1>y&+AEG?qgQjGJg*(I02LrP8_v(dcihlB#8XYfo?(Q(u9+IC|TF8IG7!5|Tn;>5O zEGY2?atVYC`9bB=mtuXiv88)0EHnugNf#Kd;%z9Z`}jR}U#Jfo_@#RpH0h_j9+L?1 zD!b1~((ZTPPfZNKS7w3PZYKokit2%=NtQh^Vxu$aC>#L0PTA{Aqc9`gR1Z>p-ws~ zLH9ugsXW_^Ox>4G#ZZ3aN2TP6WY?=NrLynLo6+bgFEF0Ni_uT3?b;+6R|g+~)dcR5 z0k%SerYvM_Kgw!jPo$aYYlc|!z5FLG3)*1JEH;Le#c5z)>3&zoDp zEu*)FQAWE43mL>8B#@b;h5)VptXm)6!z@~Qyk7oPUcg6^Ui_Hww`;7NOrZr+> zSqi;YR0w>2?R50elrL|FIF7$C6{l+bUk^D#5)!>*s*=#NRu(pWGpWhnI|s&J*PwOi zYLK6q74p8wkNDVK#R>{9<`wb}gTgh{`u@B#jOgV!knWs`zwl2??!Ic;rG8^*sjvPy zG5C$IM7LE|yhbH=tfI?EZx^3et#SR3UJD*uZ{AHmCniwEjx>a|NP8Mb-FY6( zcP_UnU8;X&uy4{@%=oxp@^bQQlT+I3BZ7*Spdn{9blS7^rQ9CST#h-o746rTBOTxD zyOlxetD{zb-ba=*Z&o^eO-@;N@Ugi`r^;TwBy64{-Lu6rBRH`t4`P;Vb_8Cl;wbNR zcyfmFIa@rSQo|3c{L6JsjT99c6^7QXMtK|UT2J$rLIU^!j3Bg9Q&cG*Dk8lp7#(==xksZo1C8hP??$h7mzAHG5z?;cB~In{-q)dS9^X ztU^5vefUvVHmXTVOeZWENlzbM##Ote_HcLVU9_>m_eYkFn=)votiBZ@tQ+;`+fMyo z`r>e4?L)$WGDAQjS<5EQW)x5Dv8mFc9@8;K=v9En1Ttn~yaa7|I_4e$yiY#jyHS}y zzZUe7-mR04$Zl?wQ@NA1ba%#YZ4A=RgSR{N1Km%hjfw}hGyGQ$*+H?NLvjfO&-Ss_HpDzqj9>&5n2LI&-}W* zuPy=t+M0B~$rb7*`e8D2-B_0_$wO{TgU|ZgLl`GSu>(UnQ@t{5j8#>t~dWZ{yP1 zp|t`&ymqC}5;t(r8uTgGS6a2pv{ZG@Ti9%%vRb~WH>;roxF0){L#!Op zV?XwMyGG2DtyP!)?Lw{ZWjqy^aDQXU8c1lI4{y-3=O(hnPv!Sy5a689WIyCxmoNn- z>hLHEA}mla;~v&W{Io$nkZkv774A=;0aE<`rY|RbQ>K8}WQsi)E~nGb1t&p8LShhT z;{HfRw)s?3@{Wns#Oc4}+T*i>Q(2{pRr(JBk5A0kq2;iGB;mCm=OheD(gbIQQa65c zu6R!*@;CP9;K?n{Zc8o9~PRr{W;hbJqfqQX>1!c zfDNzE*>egxtB6-?7y$*+RCYy5i8JhWi5Z0OrYI);B(Vy?`}%)( zDT|&OTF=b1I)2DL*~_ASeDaED_rtHz%}e`R>2|l3OWa{{N3E!X?ukoN8@XtqmHly_ zi)EbYc4K}B>iE*Ru{hsZXSAcQE5258ED9l+;kc(KE^WSBR+=-rH%B0m>Qsosk{V$Z zrhPS6C2N%Ro<#5Sd(gegG;9+j3t_N}|0c5C|2pk|XI>7ID0fx3s;3vwm`BemKj zA-Ty@gX^)}88up~8xJf3FJUpf*z+`3`zGSh)+pkdk2uyF)6uTP$$dni=4l%a#?|@W zBfOA5;*?5nYXaUt$SY>nJI!sDUKufedg#0I&8=m6bcIjwy4v5H`jPdDI$@`n(Y4?6%3klmFA@vFk1CXh& zX$@DaJaKclBXz7xr(ciN#*Zw{sX4MR{*I`&UB7c3p7Shu#7vPIvC;XsPWx?SU!l{U z!|U#D7vm?UZj5J^kJ1D5v?eN|Fz^RqrZ!q_jkX^NVlG?wYf^ssRmJ-kg;(EV@Sv}? z`Yr4~aRNzJ3?qgP?Ufj;Oro8eNrb#{njC~aF^ukG@bnk6-0xk+;26B?ba$fL#b0Gj zzaKr2oIWgI(T2y~d`RPN0_qeP6VkQJ78{9Ibrtz^?I#}n}-TN~} zg|$8K&Zq0FeCSz3Cw_FzWx5sV1XP6kwuh7l6}Y#x5!X}*ti<>jv#4lX$R;LqUf0`y zR+r-5HflwY1;6$~kG-5<{8`b4?w5g=An$Z%Sxw47zB~eBu^bfYor4$9OiTylD~xk` zY6*`Z5UO?Hj zKXY#N8}(K>o!B7>3w68+4!{t+ha(_B0qP(>m1E1D&rSXk#IH;kbB0l3vVEjyzon2LAW)Drq zE1z(vu!a~o@D2z{kWbHNd>QJub*2gH${+~F%7l&~_UoTL0E<`KT%=ZeD4TRjGi$xu zWTQO05{ru^_O*TPEoGO_Se;%Pnay1rc2QvER%ZFo+O64GZmcYweE)UZ*ppL#+EV)^ zMU_o(LH6Y)ZLyxA2b8J5GFddQHSCQ(yvhf#tHQvGhd+K#OH8I;>H(o;4Nh#4tr z9qsVh3no*j(4^Cp_DCWddR%qR19N)M8&2YV%YarJRwuA~GQ4rldQKrgw38xjD{VN= zcI0t=nG0&+-@eO&lVP}iPMhx(F!vKn4mKqqQ?`qJG1*qs!~F|{mevq0C1e*7WVUaa zOyR1v9<3~P>By8EU;g|U$Zr))uTmVSdlMkzWmB1>)CRdJ94--s%YrY~G&+qJIq!$& zbA9c^{q)0uD9PrJFY+Xu6eS#krpmnD?~foXF?llq>fmY7n9l+dmBSvlQYUw_*r}a^ zw#%$+X>~L5;ZEfmN5}661p=I;(69~i*CKRz#^DG4Ff<19Zg8+=6@_FBRcI$U+vh%e z9?K}MEo`a&%R?N%76UZtTBeRMl2Ej$XPs%_v(Da`i!~YU*0VmJscGhjl(*L!?#Lk0 zaq;9|)}xnbBPO+nky6*(M`mb_IRII^hY|TNMdqululI3)L8iS(R>t%I3RmgBR(O*#MN9v4@4-03 z_sjmBW(%9q4@_Sl21UYCz8N7G015lGC4CgMSgmUJIMv+W{Q4yaTAXW_=`HIu?VdlT zgOkxOS0QS2n${lS_EMO*<&DJ-S;Qr>!2*EUq3q=-+ym_Om2-{J*r`GA-@%-#!qZ_} z9BsHir1!Y6T16$gw-eBvgWI@!nD8Uh-^J#FcP~k`0b`nj59W^95fc)~?|Eb1!@k3z zaLgg5*<)3*)8*%(n)ay72K(NmQ;ISqF7>&)GJNA|n;`sYRo#C}Y3WmzNM?W(h)0U>?MlL5FaNxGSj zJtvgJ+7$AV=l@$6FP@bg1Lmj9H%tZQw4U2ld$=FHYis(9_&ork3q=xacms?8*YKKC%ds8~ zAf9CxvGs?g_vBvmRv{Nft#4UG33nHhyQd?BZj)1Ep8Z?*>^reNEhNWs`p^AMyiB z!MawVk@G?U0u6`o)Auvpica^X-HICm>=uhWGmW}@KRM&JzKGl0i$rZWe-<19a*(Z_ zbSB$!CIV{il&Q(dK#|EdSLkvOHyWSG=c4^^b#2Yp>m!VR15b(DFx!8g#$BHPdv%ZY z3f0TA!v&C4Ey2ffL-0J@0oi!~^Y9iqCX8O0Kto4~;Tx1y4p*;v8OW5y7dhIESGtW>vE;(JLuR}@JZBr7+-z+1<>g~5E6Y{2Uo#j5bJDLw#LGEEyIn*) zl(-x$h7nKi3`Xa3Q=e4jrgsZ{lyQ#CRfeqcR%Of4V3@v@w& zd>8s#eMP29oLGd?PI+eX*9v?=_&x182@X^HL_ z?{}-LE6y)wGeJ&h%osJdUJT#sz?w4N585R9no-r(w!FBQGIa82oHDTwAPFbN#&khG zP?&gfT--)rgRz+pqy6K7A%Ij~-?fsR%hB`ni{G4P)(>~7(rCZtYUvq(`TE6+1k=1J zPS*Cx-%`#QcWQ_Zd3(h{uTP{?BrYKdrGm;ioCG=r5fKqUXFIdNx&n7fbK-Brl3|Ol zAW)}P-*S`fd>A@070ws{1EN(rqdxxvv)ANTrLCv>Lf$%#B_@RP`Xm)kV!G1=G(nXx z#V`>EGSv&KXF>6fOa9H&I;mD@Bn*5^Cyaz6=<9@ZcRH+mjU40v=fdb9D2;q)BDp-N zy2UT=ABc;d50k+oKcN@M(Jc2S+IB*nMNKFF0vKu5Cf?OL5Oo(NE?}nt*UM(kt!>%0 zWeUB~!8^#~lQSwBjEJ5A&C6ePwhYTyN75*)kHqOch>XEqu#D8yLmQ;{9*99jtZ0ax z>EEQ_V0=c2dubFiEsj3{(}yWo+K}GCzy$poU4tV=3&$ky_Bc%~mY?m6iHUiB_WL1X zavZFqCMz2m9gS)=WWHwcRaad6XA0(EbwvdZN~l5!^q@_xTmu~)9TxNeD&-%K7wzl=Y zXD|Q@lK6z3xAi18C!S)!%tA-!@{JPmxB%1zb8~a?A8c?ecVLK|yDr~_0X^2Q9mS|- zvH;w_yu5tax}U@+j)eDWQQP?C^$C*796?P@LsO&I95s0!->0Reg^!14=_e&v33t~E z+_nl1Fq#wwCK+gb0N$v95CjHBL>BA+|zmpCOQv< zBFNi3({XB936*PV+OW$fLqj1S2qWyoF1#`{1Pj$v9oJvZf;@|A6@ac|HG&OZ}f8$=HEw3;`?_HZlT$ zZ6o*l_rrj>)jg7thMIaGcnq!XUnb9i`K^UgvA{*;@O;|=;4VytB+>QL#>Vg`KXOH3 zJ{Pub1^Aw?HUk4+fQhl`{GJyjidg`Z%j{T8H%zU39T=heaDB*Rh}FpFb940ka9vnf z2+&yA<45Ct{HU#=VF7$E$m8{qNCopo2~kFhK4o=zxk|rPsd4||Fz?WX$F+SuxcF-i zFcVf;_{IJC{t~$WK-cMY@K~iJB(OTJazJ;Nmb5%Q?*NeYN`vHAcKcW8Vm=w`yQW^5 z(zdp?`ucB6Rl1o>fOU0rbOb>2_>3tJ#igZ?tsdiYjfNEm`Lv;`larI_Y+m-jS3jsa zMmW3M!6bmP)Nj$SvTE+>ezyQW926=fJF>V4-30kMIXPuA0fu_dg=czt`oE(o=KWt7 zS!7CKVd4LdCMPeCEhZ`|n#lzG>GXu8j#((s%)9P?WvG0Ca{fQp7CvP5aDQK?!w8)x z54hh>QWo;QU9qmdRa;n1PD}e7BJnBFDq4q8t6U>CHa3;v11onKQVc9?OjkukMJl6i zgUu{pV(}f==qVc0Eyj@_eE-s@m2zzYV!Sb3K!<@A#Qswv55x88gN zfU~n0v?XL`cX>jLNXU0`K<|`;J%H-91o9#nY~Qi~1~=+3o^uh_;Vfxk-r!|pOhiXU zj*XA23FAi(As3jB#QnEfj6jbAn8HoJ$lza+OqvZ04cX8D03R?er^!cNNh$6Ygk?fq zRY&$;zj5lwzAHQ&9EbA0>qCWRB|fvkJpl36tJAQO;(5QMz=%&vMMXtR3ruJxprlHh zSe)+%ydD#+C0*d`2&{9w0BeV>ZFytkHDGp{T3K}JU~2lS1_5L;FhK=mIgu)xMB@)E zn}>%-rFLzLdqI${sGtcNvPn@=QqoT>9{|Pjeyvn7>s?@*mKgwJ21N8f0|qco{4Oj6 z1~#AX?^_;=!8Qi_^b(}|1mOZ(ON5CP1siz(U+wk(xq*Q1{r}ex|2uxw+Yh)~HXN%f zjn<3Gn52Az7GIq=IxJ_!U`1fl<%%2J31RBK1C$3K*L^uXy-W>ybQBc5Zodb8XV8B) zF7h2vWp!PE_~+;618&CR(o!z4%@3xs81X0OLg9kWT;x9jTL0|)+!NqPuv?@rR`WiX zBow=`v$4TK2~DLZhA9FPVK6cUg@pq^S5IHx_H+(EL6~g|Hy0OxZIfToE&@+f4k!WsXX)}$|8F?}KK?&j9B$s!%>2&04l`7M3x$!DRFbF>GYS45_XSzA diff --git a/man/sample_chain.Rd b/man/sample_chain.Rd index c29faa2..8ca9498 100644 --- a/man/sample_chain.Rd +++ b/man/sample_chain.Rd @@ -16,21 +16,33 @@ sample_chain( ) } \arguments{ -\item{target_distribution}{Target stationary distribution for chain. A list -with named entries \code{log_density} and \code{gradient_log_density} corresponding -to respectively functions for evaluating the logarithm of the (potentially -unnormalized) density of the target distribution and its gradient (only -required for gradient-based proposals). As an alternative to -\code{gradient_log_density} an entry \code{value_and_gradient_log_density} may +\item{target_distribution}{Target stationary distribution for chain. One of: +\itemize{ +\item A one-sided formula specifying expression for log density of target +distribution which will be passed to +\code{\link[=target_distribution_from_log_density_formula]{target_distribution_from_log_density_formula()}} to construct functions +to evaluate log density and its gradient using \code{\link[=deriv]{deriv()}}. +\item A \code{bridgestan::StanModel} instance (requires \code{bridgestan} to be +installed) specifying target model and data. Will be passed to +\code{\link[=target_distribution_from_stan_model]{target_distribution_from_stan_model()}} using default values for optional +arguments - to override call \code{\link[=target_distribution_from_stan_model]{target_distribution_from_stan_model()}} +directly and pass the returned list as the \code{target_distribution} argument +here. +\item A list with named entries \code{log_density} and \code{gradient_log_density} +corresponding to respectively functions for evaluating the logarithm of +the (potentially unnormalized) density of the target distribution and its +gradient (only required for gradient-based proposals). As an alternative +to \code{gradient_log_density} an entry \code{value_and_gradient_log_density} may instead be provided which is a function returning both the value and gradient of the logarithm of the (unnormalized) density of the target -distribution as a list under the names \code{value} and \code{gradient} respectively. -The list may also contain a named entry \code{trace_function}, correspond to a -function which given current chain state outputs list of variables to trace -on each main (non-adaptive) chain iteration. If a \code{trace_function} entry is -not specified, then the default behaviour is to trace the position -component of the chain state along with the log density of the target -distribution.} +distribution as a list under the names \code{value} and \code{gradient} +respectively. The list may also contain a named entry \code{trace_function}, +correspond to a function which given current chain state outputs list of +variables to trace on each main (non-adaptive) chain iteration. If a +\code{trace_function} entry is not specified, then the default behaviour is to +trace the position component of the chain state along with the log +density of the target distribution. +}} \item{initial_state}{Initial chain state. Either a vector specifying just the position component of the chain state or a list output by \code{chain_state} From 75c3c18825dba7774880f8777757651c2b250e64 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 20:51:14 +0000 Subject: [PATCH 06/13] Use dummy variable declaration to avoid check note --- R/bridges.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/bridges.R b/R/bridges.R index 32514bc..fd8fb30 100644 --- a/R/bridges.R +++ b/R/bridges.R @@ -91,6 +91,7 @@ example_gaussian_stan_model <- function(n_data = 50, seed = 1234L) { }" withr::with_seed(seed, y <- stats::rnorm(n_data)) data_string <- sprintf('{"N": %i, "y": [%s]}', n_data, toString(y)) + model_file <- NULL # to avoid 'no visible binding for global variable' note withr::with_tempfile("model_file", { writeLines(model_string, model_file) From d1de18519a032ffa3c99d374505d1749ed3ad942 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 20:54:13 +0000 Subject: [PATCH 07/13] Qualify deriv call with stats package name --- R/bridges.R | 4 ++-- README.Rmd | 2 +- README.md | 2 +- man/target_distribution_from_log_density_formula.Rd | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/R/bridges.R b/R/bridges.R index fd8fb30..cb484e6 100644 --- a/R/bridges.R +++ b/R/bridges.R @@ -118,12 +118,12 @@ example_gaussian_stan_model <- function(n_data = 50, seed = 1234L) { #' #' @examples #' target_distribution <- target_distribution_from_log_density_formula( -#' ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 10), +#' ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 10) #' ) #' target_distribution$value_and_gradient_log_density(c(0.1, -0.3)) target_distribution_from_log_density_formula <- function(log_density_formula) { variables <- all.vars(log_density_formula) - deriv_log_density <- deriv(log_density_formula, variables, func = TRUE) + deriv_log_density <- stats::deriv(log_density_formula, variables, func = TRUE) value_and_gradient_log_density <- function(position) { names(position) <- variables value <- rlang::inject(deriv_log_density(!!!position)) diff --git a/README.Rmd b/README.Rmd index e7763cc..e306d7e 100644 --- a/README.Rmd +++ b/README.Rmd @@ -76,7 +76,7 @@ Here we use the default values of the `proposal` and `adapters` arguments to `sa corresponding respectively to the Barker proposal, and adapters for tuning the proposal scale to coerce the average acceptance rate using a dual-averaging algorithm, and for tuning the proposal shape based on an estimate of the target distribution covariance matrix. The `target_distribution` argument to `sample_chain()` is passed a formula specifying the log density of the target distribution, which is passed to `target_distribution_from_log_density_formula()` to construct necessary functions, -using `deriv()` to symbolically compute derivatives. +using `stats::deriv()` to symbolically compute derivatives. ```{r banana-samples, fig.width=6, fig.height=4} diff --git a/README.md b/README.md index fe9cb7c..7ce016d 100644 --- a/README.md +++ b/README.md @@ -84,7 +84,7 @@ estimate of the target distribution covariance matrix. The `target_distribution` argument to `sample_chain()` is passed a formula specifying the log density of the target distribution, which is passed to `target_distribution_from_log_density_formula()` to construct -necessary functions, using `deriv()` to symbolically compute +necessary functions, using `stats::deriv()` to symbolically compute derivatives. ``` r diff --git a/man/target_distribution_from_log_density_formula.Rd b/man/target_distribution_from_log_density_formula.Rd index 4fc6f5d..ac0f09b 100644 --- a/man/target_distribution_from_log_density_formula.Rd +++ b/man/target_distribution_from_log_density_formula.Rd @@ -25,7 +25,7 @@ Construct target distribution from a formula specifying log density. } \examples{ target_distribution <- target_distribution_from_log_density_formula( - ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 10), + ~ (-(x^2 + y^2) / 8 - (x^2 - y)^2 - (x - 1)^2 / 10) ) target_distribution$value_and_gradient_log_density(c(0.1, -0.3)) } From b677e13372bb4a9538b6482af8958289d8d05cb4 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 20:58:33 +0000 Subject: [PATCH 08/13] Test target distribution from formula function --- tests/testthat/test-bridges.R | 44 +++++++++++++++++++++++++++++------ 1 file changed, 37 insertions(+), 7 deletions(-) diff --git a/tests/testthat/test-bridges.R b/tests/testthat/test-bridges.R index f8ac4ea..425b18f 100644 --- a/tests/testthat/test-bridges.R +++ b/tests/testthat/test-bridges.R @@ -19,9 +19,7 @@ test_that("Creating example Gaussian Stan model works", { expect_identical(model$param_unc_num(), 2L) }) -test_that("Creating target distribution from Stan model works", { - model <- cached_example_gaussian_stan_model() - target_distribution <- target_distribution_from_stan_model(model) +check_target_distribution <- function(target_distribution) { expect_type(target_distribution, "list") expect_named( target_distribution, @@ -30,19 +28,33 @@ test_that("Creating target distribution from Stan model works", { expect_type(target_distribution$log_density, "closure") expect_type(target_distribution$value_and_gradient_log_density, "closure") expect_type(target_distribution$trace_function, "closure") - position <- rep(0, model$param_unc_num()) +} + +check_log_density_and_gradient <- function( + position, target_distribution, true_log_density) { log_density <- target_distribution$log_density(position) value_and_gradient_log_density <- ( target_distribution$value_and_gradient_log_density(position) ) expect_type(log_density, "double") - expect_identical(log_density, model$log_density(position)) + expect_equal(log_density, true_log_density(position)) expect_type(value_and_gradient_log_density, "list") expect_identical(value_and_gradient_log_density$value, log_density) expect_equal( value_and_gradient_log_density$gradient, - numerical_gradient(target_distribution$log_density)(position), - tolerance = 1e-6 + numerical_gradient(true_log_density)(position), + tolerance = 1e-6, + ignore_attr = TRUE + ) +} + +test_that("Creating target distribution from Stan model works", { + model <- cached_example_gaussian_stan_model() + target_distribution <- target_distribution_from_stan_model(model) + check_target_distribution(target_distribution) + position <- rep(0, model$param_unc_num()) + check_log_density_and_gradient( + position, target_distribution, model$log_density ) }) @@ -87,3 +99,21 @@ for (include_log_density in c(TRUE, FALSE)) { } ) } + +test_that("Constructing target distribution from log density formula works", { + log_density_formula <- ~ -(x^2 + y^2 + z^2) / 2 + target_distribution <- target_distribution_from_log_density_formula( + log_density_formula + ) + check_target_distribution(target_distribution) + withr::with_seed(seed = default_seed(), position <- rnorm(3)) + check_log_density_and_gradient( + position, target_distribution, function(x) -sum(x^2) / 2 + ) + trace_values <- target_distribution$trace_function(chain_state(position)) + expect_named(trace_values, c("x", "y", "z", "log_density")) + expect_equal( + trace_values, c(position, -sum(position^2) / 2), + ignore_attr = TRUE + ) +}) From 4a8e4f416d89c5c470ed4d5f0edac42377ccbd1e Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 21:00:32 +0000 Subject: [PATCH 09/13] Use base::inherits in place of methods::is --- R/chains.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/chains.R b/R/chains.R index 8131791..4d380c4 100644 --- a/R/chains.R +++ b/R/chains.R @@ -113,11 +113,11 @@ sample_chain <- function( } else { stop("initial_state must be a vector or list with an entry named position.") } - if (is(target_distribution, "formula")) { + if (inherits(target_distribution, "formula")) { target_distribution <- target_distribution_from_log_density_formula( target_distribution ) - } else if (is(target_distribution, "StanModel")) { + } else if (inherits(target_distribution, "StanModel")) { target_distribution <- target_distribution_from_stan_model( target_distribution ) From db07fd221d3d430414cab88e2ffd24008308ca3b Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Wed, 11 Dec 2024 21:02:04 +0000 Subject: [PATCH 10/13] Remove removed trace_function argument from sample_chain docs Clarify trace_function allowable output types --- R/chains.R | 12 +++++------- man/sample_chain.Rd | 13 +++++-------- 2 files changed, 10 insertions(+), 15 deletions(-) diff --git a/R/chains.R b/R/chains.R index 4d380c4..59c5388 100644 --- a/R/chains.R +++ b/R/chains.R @@ -24,11 +24,11 @@ #' gradient of the logarithm of the (unnormalized) density of the target #' distribution as a list under the names `value` and `gradient` #' respectively. The list may also contain a named entry `trace_function`, -#' correspond to a function which given current chain state outputs list of -#' variables to trace on each main (non-adaptive) chain iteration. If a -#' `trace_function` entry is not specified, then the default behaviour is to -#' trace the position component of the chain state along with the log -#' density of the target distribution. +#' correspond to a function which given current chain state outputs a named +#' vector or list of variables to trace on each main (non-adaptive) chain +#' iteration. If a `trace_function` entry is not specified, then the default +#' behaviour is to trace the position component of the chain state along +#' with the log density of the target distribution. #' @param initial_state Initial chain state. Either a vector specifying just #' the position component of the chain state or a list output by `chain_state` #' specifying the full chain state. @@ -56,8 +56,6 @@ #' coerce the average acceptance rate to a target value using a dual-averaging #' algorithm, and adapting the shape to an estimate of the covariance of the #' target distribution. -#' @param trace_function Function which given current chain state outputs list -#' of variables to trace on each main (non-adaptive) chain iteration. #' @param show_progress_bar Whether to show progress bars during sampling. #' Requires `progress` package to be installed to have an effect. #' @param trace_warm_up Whether to record chain traces and adaptation / diff --git a/man/sample_chain.Rd b/man/sample_chain.Rd index 8ca9498..4e7c774 100644 --- a/man/sample_chain.Rd +++ b/man/sample_chain.Rd @@ -37,11 +37,11 @@ instead be provided which is a function returning both the value and gradient of the logarithm of the (unnormalized) density of the target distribution as a list under the names \code{value} and \code{gradient} respectively. The list may also contain a named entry \code{trace_function}, -correspond to a function which given current chain state outputs list of -variables to trace on each main (non-adaptive) chain iteration. If a -\code{trace_function} entry is not specified, then the default behaviour is to -trace the position component of the chain state along with the log -density of the target distribution. +correspond to a function which given current chain state outputs a named +vector or list of variables to trace on each main (non-adaptive) chain +iteration. If a \code{trace_function} entry is not specified, then the default +behaviour is to trace the position component of the chain state along +with the log density of the target distribution. }} \item{initial_state}{Initial chain state. Either a vector specifying just @@ -82,9 +82,6 @@ Requires \code{progress} package to be installed to have an effect.} \item{trace_warm_up}{Whether to record chain traces and adaptation / transition statistics during (adaptive) warm-up iterations in addition to (non-adaptive) main chain iterations.} - -\item{trace_function}{Function which given current chain state outputs list -of variables to trace on each main (non-adaptive) chain iteration.} } \value{ A list with entries From c667b55693595374a7c1a196af81c7a0d1c57692 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 12:20:17 +0000 Subject: [PATCH 11/13] Test using invalid target distribution with sample_chain raises error --- tests/testthat/test-chains.R | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/tests/testthat/test-chains.R b/tests/testthat/test-chains.R index 944c2c9..4cc9be0 100644 --- a/tests/testthat/test-chains.R +++ b/tests/testthat/test-chains.R @@ -89,3 +89,15 @@ test_that("Sample chains with invalid initial_state raises error", { "initial_state" ) }) + +test_that("Sample chains with invalid target_distribution raises error", { + expect_error( + sample_chain( + target_distribution = list(), + initial_state = c(0., 0.), + n_warm_up_iteration = 1, + n_main_iteration = 1, + ), + "target_distribution" + ) +}) From 2c23ae13165609454c68cae03dedb021246fcad8 Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 12:21:00 +0000 Subject: [PATCH 12/13] Test passing explicit trace function to sample_chain --- tests/testthat/test-chains.R | 130 +++++++++++++++++++---------------- 1 file changed, 71 insertions(+), 59 deletions(-) diff --git a/tests/testthat/test-chains.R b/tests/testthat/test-chains.R index 4cc9be0..38c2809 100644 --- a/tests/testthat/test-chains.R +++ b/tests/testthat/test-chains.R @@ -4,70 +4,82 @@ for (n_warm_up_iteration in c(0, 1, 10)) { for (trace_warm_up in c(TRUE, FALSE)) { for (show_progress_bar in c(TRUE, FALSE)) { for (wrapped_initial_state in c(TRUE, FALSE)) { - test_that( - sprintf( - paste0( - "Sampling chain with %i warm-up iterations, ", - "%i main iterations, dimension %i, ", - "trace_warm_up = %i, show_progress_bar = %i ", - "and wrapped_initial_state = %i works" + for (explicit_trace_function in c(TRUE, FALSE)) { + test_that( + sprintf( + paste0( + "Sampling chain with %i warm-up iterations, ", + "%i main iterations, dimension %i, ", + "trace_warm_up = %i, show_progress_bar = %i ", + "wrapped_initial_state = %i and ", + "explicit_trace_function = %i works" + ), + n_warm_up_iteration, + n_main_iteration, + dimension, + trace_warm_up, + show_progress_bar, + wrapped_initial_state, + explicit_trace_function ), - n_warm_up_iteration, - n_main_iteration, - dimension, - trace_warm_up, - show_progress_bar, - wrapped_initial_state - ), - { - target_distribution <- standard_normal_target_distribution() - adapters <- list( - scale_adapter("stochastic_approximation", initial_scale = 1.) - ) - withr::with_seed(default_seed(), { - position <- rnorm(dimension) - }) - if (wrapped_initial_state) { - initial_state <- chain_state(position) - } else { - initial_state <- position - } - results <- sample_chain( - target_distribution = target_distribution, - initial_state = initial_state, - n_warm_up_iteration = n_warm_up_iteration, - n_main_iteration = n_main_iteration, - adapters = adapters, - trace_warm_up = trace_warm_up, - show_progress_bar = show_progress_bar - ) - expected_results_names <- c( - "final_state", "traces", "statistics" - ) - if (trace_warm_up) { + { + target_distribution <- standard_normal_target_distribution() + if (explicit_trace_function) { + target_distribution[["trace_function"]] <- ( + default_trace_function(target_distribution) + ) + } + adapters <- list( + scale_adapter( + "stochastic_approximation", + initial_scale = 1. + ) + ) + withr::with_seed(default_seed(), { + position <- rnorm(dimension) + }) + if (wrapped_initial_state) { + initial_state <- chain_state(position) + } else { + initial_state <- position + } + results <- sample_chain( + target_distribution = target_distribution, + initial_state = initial_state, + n_warm_up_iteration = n_warm_up_iteration, + n_main_iteration = n_main_iteration, + adapters = adapters, + trace_warm_up = trace_warm_up, + show_progress_bar = show_progress_bar + ) expected_results_names <- c( + "final_state", "traces", "statistics" + ) + if (trace_warm_up) { + expected_results_names <- c( + expected_results_names, + "warm_up_traces", + "warm_up_statistics" + ) + } + expect_named( + results, expected_results_names, - "warm_up_traces", - "warm_up_statistics" + ignore.order = TRUE, ) + expect_nrow(results$traces, n_main_iteration) + expect_ncol(results$traces, dimension + 1) + expect_nrow(results$statistics, n_main_iteration) + expect_ncol(results$statistics, 1) + if (trace_warm_up) { + expect_nrow(results$warm_up_traces, n_warm_up_iteration) + expect_ncol(results$warm_up_traces, dimension + 1) + expect_nrow(results$warm_up_statistics, n_warm_up_iteration) + expect_ncol(results$warm_up_statistics, 2) + } } - expect_named( - results, - expected_results_names, - ignore.order = TRUE, - ) - expect_nrow(results$traces, n_main_iteration) - expect_ncol(results$traces, dimension + 1) - expect_nrow(results$statistics, n_main_iteration) - expect_ncol(results$statistics, 1) - if (trace_warm_up) { - expect_nrow(results$warm_up_traces, n_warm_up_iteration) - expect_ncol(results$warm_up_traces, dimension + 1) - expect_nrow(results$warm_up_statistics, n_warm_up_iteration) - expect_ncol(results$warm_up_statistics, 2) - } - } - ) + ) + } } } } From 6efa4f1b98f1b17594743e6adee3997bf97e84bd Mon Sep 17 00:00:00 2001 From: Matt Graham Date: Mon, 6 Jan 2025 12:21:29 +0000 Subject: [PATCH 13/13] Test using sample_chain with Stan model and log density formula works --- tests/testthat/test-bridges.R | 44 +++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/tests/testthat/test-bridges.R b/tests/testthat/test-bridges.R index 425b18f..403b1c5 100644 --- a/tests/testthat/test-bridges.R +++ b/tests/testthat/test-bridges.R @@ -100,6 +100,28 @@ for (include_log_density in c(TRUE, FALSE)) { ) } +test_that("Using sample_chains with Stan model as target_distribution works", { + model <- cached_example_gaussian_stan_model() + dimension <- model$param_unc_num() + n_warm_up_iteration <- 50 + n_main_iteration <- 50 + results <- sample_chain( + target_distribution = model, + initial_state = rep(0, dimension), + n_warm_up_iteration = n_warm_up_iteration, + n_main_iteration = n_main_iteration + ) + expect_named( + results, + c("final_state", "traces", "statistics"), + ignore.order = TRUE, + ) + expect_nrow(results$traces, n_main_iteration) + expect_ncol(results$traces, dimension + 1) + expect_nrow(results$statistics, n_main_iteration) + expect_ncol(results$statistics, 1) +}) + test_that("Constructing target distribution from log density formula works", { log_density_formula <- ~ -(x^2 + y^2 + z^2) / 2 target_distribution <- target_distribution_from_log_density_formula( @@ -117,3 +139,25 @@ test_that("Constructing target distribution from log density formula works", { ignore_attr = TRUE ) }) + +test_that("Using sample_chains with formula as target_distribution works", { + log_density_formula <- ~ -(x^2 + y^2 + z^2) / 2 + dimension <- 3 + n_warm_up_iteration <- 50 + n_main_iteration <- 50 + results <- sample_chain( + target_distribution = log_density_formula, + initial_state = rep(0, dimension), + n_warm_up_iteration = n_warm_up_iteration, + n_main_iteration = n_main_iteration + ) + expect_named( + results, + c("final_state", "traces", "statistics"), + ignore.order = TRUE, + ) + expect_nrow(results$traces, n_main_iteration) + expect_ncol(results$traces, dimension + 1) + expect_nrow(results$statistics, n_main_iteration) + expect_ncol(results$statistics, 1) +})