From d69ccd2792ea0da22f1955aee6da54058e8618ca Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Thu, 4 Apr 2024 10:24:27 -0700 Subject: [PATCH 01/19] small updates to vignettes --- R/utility.R | 5 ++++- vignettes/dsem.Rmd | 5 +++++ vignettes/simultaneous_autoregressive_process.Rmd | 2 +- 3 files changed, 10 insertions(+), 2 deletions(-) diff --git a/R/utility.R b/R/utility.R index 01a0430..c9de36b 100644 --- a/R/utility.R +++ b/R/utility.R @@ -322,7 +322,10 @@ function( x, #' #' This implementation conditions upon the maximum likelihood estimate of fixed effects #' and the empirical Bayes ("plug-in") prediction of random effects. It can -#' be described as "conditional deviance explained". +#' be described as "conditional deviance explained". A state-space model that +#' estimates measurement error variance approaching one (i.e., collapses to +#' a process-error-only model) will have a conditional deviance explained +#' that approaches 1.0 #' #' @param x output from `\code{tinyVAST()}` #' @param null_formula formula for the null model. If missing, it uses diff --git a/vignettes/dsem.Rmd b/vignettes/dsem.Rmd index b5fd1ca..d4cba7e 100644 --- a/vignettes/dsem.Rmd +++ b/vignettes/dsem.Rmd @@ -59,6 +59,11 @@ mytiny = tinyVAST( dsem = dsem, formula = logn ~ 0 + var ) mytiny +# Deviance explained relative to both intercepts +# Note that a process-error-only estimate with have devexpl -> 1 +deviance_explained( mytiny, + null_formula = logn ~ 0 + var ) + # See summary knitr::kable( summary(mytiny,"dsem"), digits=3 ) ``` diff --git a/vignettes/simultaneous_autoregressive_process.Rmd b/vignettes/simultaneous_autoregressive_process.Rmd index 1d92ae5..91c4517 100644 --- a/vignettes/simultaneous_autoregressive_process.Rmd +++ b/vignettes/simultaneous_autoregressive_process.Rmd @@ -136,7 +136,7 @@ plot( adjacency_graph, ``` We can then pass this adjacency graph to `tinyVAST` during fitting: -```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} +```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6, warning=FALSE} # Fit tinyVAST model mytiny = tinyVAST( formula = Biomass_nozeros ~ 0 + Species + Region, From 5950b9da12a40214366f1c183095931c4716c306 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Thu, 11 Apr 2024 12:33:42 -0700 Subject: [PATCH 02/19] update docs --- R/predict.R | 97 ++++++- R/utility.R | 2 +- man/deviance_explained.Rd | 5 +- man/integrate_output.Rd | 97 ++++++- vignettes/dsem.Rmd | 4 +- vignettes/multiple_data.Rmd | 2 +- vignettes/spatial.Rmd | 27 +- vignettes/web_only/VAST.Rmd | 251 ++++++++++-------- vignettes/web_only/condition.Rmd | 18 +- .../empirical_orthogonal_functions.Rmd | 1 - vignettes/web_only/figure/EOF-1.png | Bin 13110 -> 13029 bytes vignettes/web_only/figure/agecomp-1.png | Bin 23312 -> 23342 bytes .../web_only/figure/condition-maps-1.png | Bin 5131 -> 4674 bytes .../figure/condition-timeseries-1.png | Bin 5860 -> 5890 bytes 14 files changed, 345 insertions(+), 159 deletions(-) diff --git a/R/predict.R b/R/predict.R index 74c1d95..1fffdef 100644 --- a/R/predict.R +++ b/R/predict.R @@ -42,6 +42,7 @@ function( object, random = c("gamma_k","epsilon_stc","omega_sc","gamma2_k","epsilon2_stc","omega2_sc"), profile = object$internal$control$profile, DLL = "tinyVAST" ) + newobj$env$beSilent() out = newobj$report()[[what]] # Add standard errors @@ -66,16 +67,100 @@ function( object, #' #' @description Applies Monte Carlo integration to approximate area-expanded abundance #' -#' @inheritParams TMB::sdreport #' @inheritParams add_predictions #' -#' @param area value used for area-weighted expansion of estimated density surface -#' for each row of `newdata`. -#' @param V_gz Settings for expansion. -#' @param W_gz Covariates for expansion. +#' @param newdata New data-frame of independent variables used to predict the response, +#' where a total value is calculated by combining across these individual predictions. +#' If these locations are randomly drawn from a specified spatial domain, then +#' `integrate_output` applies Monte Carlo integration to approximate the +#' total over that area. If locations are drawn sysmatically from a domain, +#' then `integrate_output` is applying a midpoint approximation to the integral. +#' @param area vector of values used for area-weighted expansion of +#' estimated density surface for each row of `newdata` +#' which must have same lnegth as number of rows of `newdata`. +#' @param V_gz Integer-matrix providing settings for expansion, with two columns +#' and as many rows as `newdata`. If `V_gz` is missing, the first column is `1s`, +#' which indicates to expand each row by the first column of `W_gz`. If the +#' first column is `2`, then it weights the second +#' column of `W_gz` by the proportion at that row, e.g., to calculate +#' an abundance-weighted covariate. +#' If the first column is a `3`, then it weights that +#' row of `newdata` by a previous row of `newdata`, specified via the 2nd +#' column of `V_gz`. This 3rd option is used to weight a prediction for +#' one category based on predicted density of another category, e.g., +#' to calculate abundance-weighted condition in a bivariate model. +#' If the first column is a `0`, then that row of `newdata` is not included +#' when calculating the total across rows. Including a row of `newdata` with +#' first-column of `0` is useful, e.g., when calculating abundance at that +#' location, but where the eventual index uses abundance as weighting term +#' but without otherwise using the predicted density in calculating a total +#' value. +#' @param W_gz Matrix providing covariates for expansion, with two columns and +#' as many rows as `newdata`. If `W_gz` is missing, the first column is +#' populated from user-supplied values for `area`. +#' The second column can be used to provide a covariate +#' that are specified for use in expansion, e.g., to provide positional +#' coordinates when calculating the abundance-weighted centroid with respect +#' to that coordinate. +#' @param bias.correct logical indicating if bias correction should be applied using +#' standard methods in [TMB::sdreport()] #' @param intern Do Laplace approximation on C++ side? Passed to [TMB::MakeADFun()]. -#' @param apply.epsilon Apply epsilon bias correction? +#' @param apply.epsilon Apply epsilon bias correction using a manual calculation +#' rather than using the conventional method in [TMB::sdreport]? Using +#' `apply.epsilon` is sometimes substantially faster. +#' +#' @details +#' Analysts will often want to calculate some value by combing the predicted response at multiple +#' locations, and potentially from multiple variables in a multivariate analysis. +#' This arises in a univariate model, e.g., when calculating the integral under a predicted +#' density function, which is approximated using a midpoint or Monte Carlo approximation +#' by calculating the linear predictors at each location `newdata`, +#' applying the inverse-link-trainsformation, +#' and calling this predicted response `mu_g`. Total abundance is then be approximated +#' by multiplying `mu_g` by the area associated with each midpoint or Monte Carlo +#' approximation point (supplied by argument `area`), +#' and summing across these area-expanded values. +#' +#' In more complicated cases, an analyst can then use `V_gz` and +#' `W_gz` to calculate the weighted average +#' of a covariate (e.g., positional coordinates) for each Monte Carlo point. +#' Alternatively, an analyst fitting a multivariate model might weight one variable +#' based on another, e.g., to calculate abundance-weighted average condition, or +#' predator-expanded stomach contents. +#' +#' In practice, supplying `V_gz` and `W_gz` requires two passes through the rows of +#' `newdata` when calculating a total value. In the following, we +#' use write equations using C++ indexing conventions such that indexing starts with 0, +#' to match the way that `integrate_output` expects indices to be supplied. +#' Given `mu_g` for each row of `newdata`, +#' the first pass calculates: +#' +#' \deqn{ \nu_g = \mu_g W_{g,0} } +#' +#' where the total value from this first pass is calculated as: +#' +#' \deqn{ \nu^* = \sum_{g=0}^{G-1} \nu_g } +#' +#' The second pass then applies a further weighting, which depends upon \eqn{ V_{g,0} }, +#' and potentially upon \eqn{ V_{g,1} } and \eqn{ W_{g,1} }. +#' +#' If \eqn{V_{g,0} = 0} then \eqn{\phi_g = 0} +#' +#' If \eqn{V_{g,0} = 1} then \eqn{\phi_g = \nu_g} +#' +#' If \eqn{V_{g,0} = 2} then \eqn{\phi_g = W_{g,1} \frac{\nu_g}{\nu^*} } +#' +#' If \eqn{V_{g,0} = 3} then \eqn{\phi_g = \frac{\nu_{V_{g,1}}}{\nu^*} \mu_g } +#' +#' Finally, the total value from this second pass is calculated as: +#' +#' \deqn{ \phi^* = \sum_{g=0}^{G-1} \phi_g } #' +#' and \eqn{\phi^*} is outputted by `integrate_output`, +#' along with a standard error and potentially using +#' the epsilon bias-correction estimator to correct for skewness and retransformation +#' bias. +#' #' @export integrate_output <- function( object, diff --git a/R/utility.R b/R/utility.R index c9de36b..9d3c28a 100644 --- a/R/utility.R +++ b/R/utility.R @@ -323,7 +323,7 @@ function( x, #' This implementation conditions upon the maximum likelihood estimate of fixed effects #' and the empirical Bayes ("plug-in") prediction of random effects. It can #' be described as "conditional deviance explained". A state-space model that -#' estimates measurement error variance approaching one (i.e., collapses to +#' estimates measurement error variance approaching zero (i.e., collapses to #' a process-error-only model) will have a conditional deviance explained #' that approaches 1.0 #' diff --git a/man/deviance_explained.Rd b/man/deviance_explained.Rd index 2bdf821..3ecf76b 100644 --- a/man/deviance_explained.Rd +++ b/man/deviance_explained.Rd @@ -25,5 +25,8 @@ to calculate the proportion of deviance explained. This implementation conditions upon the maximum likelihood estimate of fixed effects and the empirical Bayes ("plug-in") prediction of random effects. It can -be described as "conditional deviance explained". +be described as "conditional deviance explained". A state-space model that +estimates measurement error variance approaching zero (i.e., collapses to +a process-error-only model) will have a conditional deviance explained +that approaches 1.0 } diff --git a/man/integrate_output.Rd b/man/integrate_output.Rd index bbcf854..5365e16 100644 --- a/man/integrate_output.Rd +++ b/man/integrate_output.Rd @@ -18,21 +18,104 @@ integrate_output( \arguments{ \item{object}{Output from \code{\link[=tinyVAST]{tinyVAST()}}.} -\item{newdata}{New data-frame of independent variables used to predict the response.} +\item{newdata}{New data-frame of independent variables used to predict the response, +where a total value is calculated by combining across these individual predictions. +If these locations are randomly drawn from a specified spatial domain, then +\code{integrate_output} applies Monte Carlo integration to approximate the +total over that area. If locations are drawn sysmatically from a domain, +then \code{integrate_output} is applying a midpoint approximation to the integral.} -\item{V_gz}{Settings for expansion.} +\item{V_gz}{Integer-matrix providing settings for expansion, with two columns +and as many rows as \code{newdata}. If \code{V_gz} is missing, the first column is \verb{1s}, +which indicates to expand each row by the first column of \code{W_gz}. If the +first column is \code{2}, then it weights the second +column of \code{W_gz} by the proportion at that row, e.g., to calculate +an abundance-weighted covariate. +If the first column is a \code{3}, then it weights that +row of \code{newdata} by a previous row of \code{newdata}, specified via the 2nd +column of \code{V_gz}. This 3rd option is used to weight a prediction for +one category based on predicted density of another category, e.g., +to calculate abundance-weighted condition in a bivariate model. +If the first column is a \code{0}, then that row of \code{newdata} is not included +when calculating the total across rows. Including a row of \code{newdata} with +first-column of \code{0} is useful, e.g., when calculating abundance at that +location, but where the eventual index uses abundance as weighting term +but without otherwise using the predicted density in calculating a total +value.} -\item{area}{value used for area-weighted expansion of estimated density surface -for each row of \code{newdata}.} +\item{area}{vector of values used for area-weighted expansion of +estimated density surface for each row of \code{newdata} +which must have same lnegth as number of rows of \code{newdata}.} -\item{W_gz}{Covariates for expansion.} +\item{W_gz}{Matrix providing covariates for expansion, with two columns and +as many rows as \code{newdata}. If \code{W_gz} is missing, the first column is +populated from user-supplied values for \code{area}. +The second column can be used to provide a covariate +that are specified for use in expansion, e.g., to provide positional +coordinates when calculating the abundance-weighted centroid with respect +to that coordinate.} -\item{bias.correct}{logical indicating if bias correction should be applied} +\item{bias.correct}{logical indicating if bias correction should be applied using +standard methods in \code{\link[TMB:sdreport]{TMB::sdreport()}}} -\item{apply.epsilon}{Apply epsilon bias correction?} +\item{apply.epsilon}{Apply epsilon bias correction using a manual calculation +rather than using the conventional method in \link[TMB:sdreport]{TMB::sdreport}? Using +\code{apply.epsilon} is sometimes substantially faster.} \item{intern}{Do Laplace approximation on C++ side? Passed to \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}.} } \description{ Applies Monte Carlo integration to approximate area-expanded abundance } +\details{ +Analysts will often want to calculate some value by combing the predicted response at multiple +locations, and potentially from multiple variables in a multivariate analysis. +This arises in a univariate model, e.g., when calculating the integral under a predicted +density function, which is approximated using a midpoint or Monte Carlo approximation +by calculating the linear predictors at each location \code{newdata}, +applying the inverse-link-trainsformation, +and calling this predicted response \code{mu_g}. Total abundance is then be approximated +by multiplying \code{mu_g} by the area associated with each midpoint or Monte Carlo +approximation point (supplied by argument \code{area}), +and summing across these area-expanded values. + +In more complicated cases, an analyst can then use \code{V_gz} and +\code{W_gz} to calculate the weighted average +of a covariate (e.g., positional coordinates) for each Monte Carlo point. +Alternatively, an analyst fitting a multivariate model might weight one variable +based on another, e.g., to calculate abundance-weighted average condition, or +predator-expanded stomach contents. + +In practice, supplying \code{V_gz} and \code{W_gz} requires two passes through the rows of +\code{newdata} when calculating a total value. In the following, we +use write equations using C++ indexing conventions such that indexing starts with 0, +to match the way that \code{integrate_output} expects indices to be supplied. +Given \code{mu_g} for each row of \code{newdata}, +the first pass calculates: + +\deqn{ \nu_g = \mu_g W_{g,0} } + +where the total value from this first pass is calculated as: + +\deqn{ \nu^* = \sum_{g=0}^{G-1} \nu_g } + +The second pass then applies a further weighting, which depends upon \eqn{ V_{g,0} }, +and potentially upon \eqn{ V_{g,1} } and \eqn{ W_{g,1} }. + +If \eqn{V_{g,0} = 0} then \eqn{\phi_g = 0} + +If \eqn{V_{g,0} = 1} then \eqn{\phi_g = \nu_g} + +If \eqn{V_{g,0} = 2} then \eqn{\phi_g = W_{g,1} \frac{\nu_g}{\nu^*} } + +If \eqn{V_{g,0} = 3} then \eqn{\phi_g = \frac{\nu_{V_{g,1}}}{\nu^*} \mu_g } + +Finally, the total value from this second pass is calculated as: + +\deqn{ \phi^* = \sum_{g=0}^{G-1} \phi_g } + +and \eqn{\phi^*} is outputted by \code{integrate_output}, +along with a standard error and potentially using +the epsilon bias-correction estimator to correct for skewness and retransformation +bias. +} diff --git a/vignettes/dsem.Rmd b/vignettes/dsem.Rmd index d4cba7e..b593496 100644 --- a/vignettes/dsem.Rmd +++ b/vignettes/dsem.Rmd @@ -33,7 +33,7 @@ options("tinyVAST.verbose" = FALSE) `tinyVAST` includes features to fit a dynamic structural equation model. We here show this using a bivariate vector autoregressive model for wolf and moose abundance on Isle Royale. -```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} +```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6, warning=FALSE} data(isle_royale, package="dsem") # Convert to long-form @@ -78,7 +78,7 @@ knitr::kable( B, digits=3) We can then compare this with package `dsem` -```{r, eval=requireNamespace("dsem"), echo=TRUE, message=FALSE, fig.width=6, fig.height=6} +```{r, eval=requireNamespace("dsem"), echo=TRUE, message=FALSE, fig.width=6, fig.height=6, warning=FALSE} library(dsem) # Keep in wide-form diff --git a/vignettes/multiple_data.Rmd b/vignettes/multiple_data.Rmd index 9af6af4..00f3284 100644 --- a/vignettes/multiple_data.Rmd +++ b/vignettes/multiple_data.Rmd @@ -30,6 +30,7 @@ knitr::opts_chunk$set( library(tinyVAST) library(sf) library(fmesher) +library(ggplot2) ``` `tinyVAST` is an R package for fitting vector autoregressive spatio-temporal (VAST) models using a minimal and user-friendly interface. @@ -144,7 +145,6 @@ plot( plot_grid, Similarly, we can plot the abundance index and standard errors ```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} -library(ggplot2) ggplot( data.frame("Year"=colnames(index),"Est"=index[1,],"SE"=index[2,]) ) + geom_point( aes(x=Year, y=Est)) + geom_errorbar( aes(x=Year, ymin=Est-SE, ymax=Est+SE) ) diff --git a/vignettes/spatial.Rmd b/vignettes/spatial.Rmd index da6c3a0..32b9e2b 100644 --- a/vignettes/spatial.Rmd +++ b/vignettes/spatial.Rmd @@ -30,6 +30,8 @@ knitr::opts_chunk$set( library(tinyVAST) library(mgcv) library(fmesher) +library(pdp) # approx = TRUE gives effects for average of other covariates +library(lattice) set.seed(101) options("tinyVAST.verbose" = FALSE) ``` @@ -67,12 +69,17 @@ plot(mesh) Finally, we can fit these data using `tinyVAST` ```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} +# Define sem, with just one variance for the single variable +sem = " + n <-> n, spatial_sd +" + # fit model out = tinyVAST( data = Data, formula = n ~ s(w), spatial_graph = mesh, control = tinyVASTcontrol(getsd=FALSE), - sem = "" ) + sem = sem) ``` We can then calculate the area-weighted total abundance: @@ -92,15 +99,8 @@ sum( Data$z ) We can compute deviance residuals and percent-deviance explained: ```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} -R1 = sum( residuals(out, type="deviance")^2 ) - -# tinyVAST null model with just a single intercept -null = tinyVAST( data = Data, - formula = n ~ 1 ) -R0 = sum( residuals(null, type="deviance")^2 ) - # Percent deviance explained -1 - R1/R0 +out$deviance_explained ``` We can then compare this with the PDE reported by `mgcv` @@ -118,8 +118,9 @@ It is then easy to confirm that mgcv and tinyVAST give (essentially) identical P ```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} out_reduced = tinyVAST( data = Data, formula = n ~ s(w) + s(x,y) ) -R1_reduced = sum( residuals(out_reduced, type="deviance")^2 ) -1 - R1_reduced/R0 + +# Extract PDE for GAM-style spatial smoother in tinyVAST +out_reduced$deviance_explained ``` # Visualize spatial response @@ -143,10 +144,6 @@ image( x=1:n_x, y=1:n_y, z=matrix(Data$z,ncol=n_y), main="True response" ) We can also compute the marginal effect of the cyclic confounder ```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} -library(pdp) # approx = TRUE gives effects for average of other covariates -library(lattice) -#library(visreg) - # compute partial dependence plot Partial = partial( object = out, pred.var = "w", diff --git a/vignettes/web_only/VAST.Rmd b/vignettes/web_only/VAST.Rmd index e4eb4e1..fc66a0e 100644 --- a/vignettes/web_only/VAST.Rmd +++ b/vignettes/web_only/VAST.Rmd @@ -62,35 +62,22 @@ mytinyVAST = tinyVAST( dsem = "logn -> logn, 1, rho", spatial_graph = mesh, family = tweedie() ) mytinyVAST -#> $call +#> Call: #> tinyVAST(formula = n ~ 0 + factor(time), data = Data, dsem = "logn -> logn, 1, rho", #> family = tweedie(), spatial_graph = mesh) #> -#> $opt -#> $opt$par -#> alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j alpha_j -#> -0.08323604 -0.13549104 -0.10579217 -0.14499114 -0.37823869 -0.21633307 -0.41489959 -0.67168423 -0.49463135 -0.13968720 0.14836187 -0.21516693 -0.20120058 0.16887045 -#> alpha_j beta_z beta_z log_sigma log_sigma log_kappa -#> 0.30040122 0.81229112 0.40988915 -0.64868475 0.04394543 0.07228542 +#> Run time: +#> Time difference of 32.35395 secs #> -#> $opt$objective -#> [1] 1717.689 +#> Family: +#> $obs #> -#> $opt$convergence -#> [1] 0 +#> Family: tweedie +#> Link function: log #> -#> $opt$iterations -#> [1] 77 #> -#> $opt$evaluations -#> function gradient -#> 107 77 -#> -#> $opt$message -#> [1] "relative convergence (4)" #> #> -#> $sdrep #> sdreport(.) result #> Estimate Std. Error #> alpha_j -0.08323604 0.15196456 @@ -115,8 +102,31 @@ mytinyVAST #> log_kappa 0.07228542 0.10755269 #> Maximum gradient component: 0.006315957 #> -#> $run_time -#> Time difference of 19.30986 secs +#> Proportion deviance explained: +#> [1] 0.6800908 +#> +#> DSEM terms: +#> heads to from parameter start lag Estimate Std_Error z_value p_value +#> 1 1 logn logn 1 1 0.8122911 0.03708631 21.90272 2.447223e-106 +#> 2 2 logn logn 2 0 0.4098891 0.03291043 12.45469 1.318626e-35 +#> +#> Fixed terms: +#> Estimate Std_Error z_value p_value +#> factor(time)1 -0.08323604 0.1519646 -0.5477333 0.583875070 +#> factor(time)2 -0.13549104 0.1867034 -0.7257021 0.468021440 +#> factor(time)3 -0.10579217 0.2052985 -0.5153090 0.606337104 +#> factor(time)4 -0.14499114 0.2178079 -0.6656835 0.505613407 +#> factor(time)5 -0.37823869 0.2269180 -1.6668518 0.095543872 +#> factor(time)6 -0.21633307 0.2302645 -0.9394981 0.347475054 +#> factor(time)7 -0.41489959 0.2345678 -1.7687834 0.076930020 +#> factor(time)8 -0.67168423 0.2383363 -2.8182199 0.004829073 +#> factor(time)9 -0.49463135 0.2386933 -2.0722464 0.038242472 +#> factor(time)10 -0.13968720 0.2373310 -0.5885756 0.556146036 +#> factor(time)11 0.14836187 0.2364020 0.6275830 0.530277181 +#> factor(time)12 -0.21516693 0.2387359 -0.9012760 0.367441605 +#> factor(time)13 -0.20120058 0.2397921 -0.8390624 0.401434274 +#> factor(time)14 0.16887045 0.2370866 0.7122734 0.476295482 +#> factor(time)15 0.30040122 0.2366004 1.2696567 0.204206943 ``` The estimated values for `beta_z` then correspond to the simulated value for `rho` and `spatial_sd`. @@ -178,21 +188,23 @@ We can then calculate the area-weighted total abundance and compare it with its # Predicted sample-weighted total (Est = sapply( seq_len(n_t), FUN=\(t) integrate_output(mytinyVAST, newdata=subset(Data,time==t)) )) -#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] -#> Estimate 97.164902 96.643634 98.362458 101.517620 84.760587 97.538110 77.520820 59.565797 75.752785 113.562128 159.66734 120.156809 137.80745 192.64174 -#> Std. Error 7.194683 7.216494 7.309313 7.572241 6.643419 7.406226 6.207919 5.090533 6.144821 8.454587 11.09945 9.107852 10.39448 13.68245 -#> Est. (bias.correct) 102.324275 102.850496 105.043004 108.373176 90.604659 104.258111 83.102278 64.065173 81.166864 120.921460 169.20673 127.519259 145.78500 203.55063 -#> Std. (bias.correct) NA NA NA NA NA NA NA NA NA NA NA NA NA NA -#> [,15] -#> Estimate 187.86973 -#> Std. Error 12.88189 -#> Est. (bias.correct) 200.54754 -#> Std. (bias.correct) NA +#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] +#> Estimate 97.164902 96.643634 98.362458 101.517620 84.760587 97.538110 77.520820 59.565797 75.752785 113.562128 159.66734 +#> Std. Error 7.194683 7.216494 7.309313 7.572241 6.643419 7.406226 6.207919 5.090533 6.144821 8.454587 11.09945 +#> Est. (bias.correct) 102.324275 102.850496 105.043004 108.373176 90.604659 104.258111 83.102278 64.065173 81.166864 120.921460 169.20673 +#> Std. (bias.correct) NA NA NA NA NA NA NA NA NA NA NA +#> [,12] [,13] [,14] [,15] +#> Estimate 120.156809 137.80745 192.64174 187.86973 +#> Std. Error 9.107852 10.39448 13.68245 12.88189 +#> Est. (bias.correct) 127.519259 145.78500 203.55063 200.54754 +#> Std. (bias.correct) NA NA NA NA # True (latent) sample-weighted total (True = tapply( Data$z, INDEX=Data$time, FUN=sum )) -#> 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 -#> 99.21643 100.10603 101.66846 109.52622 85.76973 100.97116 80.99847 68.60738 85.39974 119.62380 147.41437 122.00580 158.26179 200.56813 203.37545 +#> 1 2 3 4 5 6 7 8 9 10 11 12 13 +#> 99.21643 100.10603 101.66846 109.52622 85.76973 100.97116 80.99847 68.60738 85.39974 119.62380 147.41437 122.00580 158.26179 +#> 14 15 +#> 200.56813 203.37545 # Index = data.frame( time=seq_len(n_t), t(Est), True ) @@ -201,7 +213,6 @@ Index$high = Index[,'Est...bias.correct.'] + 1.96*Index[,'Std..Error'] # library(ggplot2) -#> Warning: package 'ggplot2' was built under R version 4.3.3 ggplot(Index, aes(time, Estimate)) + geom_ribbon(aes(ymin = low, ymax = high), # shadowing cnf intervals @@ -215,6 +226,9 @@ ggplot(Index, aes(time, Estimate)) + Next, we compare this against the current version of VAST +``` +#> Warning: package 'marginaleffects' was built under R version 4.3.3 +``` ```r settings = make_settings( purpose="index3", @@ -246,10 +260,12 @@ myVAST = fit_model( settings=settings, myVAST #> fit_model(.) result #> $par -#> beta1_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft -#> -0.58930073 0.51556564 0.46458138 0.49054919 0.46815476 0.22995172 0.39666650 0.19322490 -0.07072448 0.11680345 0.47101097 -#> beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft L_epsilon2_z logkappa2 Epsilon_rho2_f logSigmaM -#> 0.77135921 0.40831814 0.42839450 0.79998255 0.91170568 0.49266123 -4.30062279 0.85035927 0.10422305 +#> beta1_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft +#> -0.58930072 0.51556558 0.46458132 0.49054914 0.46815469 0.22995166 0.39666648 0.19322487 -0.07072451 +#> beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft beta2_ft L_epsilon2_z logkappa2 +#> 0.11680346 0.47101097 0.77135916 0.40831806 0.42839443 0.79998249 0.91170565 0.49266122 -4.30062282 +#> Epsilon_rho2_f logSigmaM +#> 0.85035928 0.10422306 #> #> $objective #> [1] 1738.337 @@ -262,10 +278,10 @@ myVAST #> 12 7 #> #> $time_for_MLE -#> Time difference of 0.9828639 secs +#> Time difference of 1.733057 secs #> #> $max_gradient -#> [1] 0.0005696325 +#> [1] 0.0005693192 #> #> $Convergence_check #> [1] "The model is likely not converged" @@ -279,57 +295,57 @@ myVAST #> #> $diagnostics #> Param starting_value Lower MLE Upper final_gradient -#> 1 beta1_ft -0.58930212 -Inf -0.58930073 Inf -2.567181e-04 -#> 2 beta2_ft 0.51556572 -Inf 0.51556564 Inf 7.261656e-06 -#> 3 beta2_ft 0.46458084 -Inf 0.46458138 Inf -1.594548e-05 -#> 4 beta2_ft 0.49055254 -Inf 0.49054919 Inf 2.243884e-04 -#> 5 beta2_ft 0.46815189 -Inf 0.46815476 Inf -2.138462e-04 -#> 6 beta2_ft 0.22994937 -Inf 0.22995172 Inf -1.954432e-04 -#> 7 beta2_ft 0.39666757 -Inf 0.39666650 Inf 7.243535e-05 -#> 8 beta2_ft 0.19322636 -Inf 0.19322490 Inf 1.265298e-04 -#> 9 beta2_ft -0.07072352 -Inf -0.07072448 Inf 8.136321e-05 -#> 10 beta2_ft 0.11680233 -Inf 0.11680345 Inf -9.143842e-05 -#> 11 beta2_ft 0.47100948 -Inf 0.47101097 Inf -1.191484e-04 -#> 12 beta2_ft 0.77135985 -Inf 0.77135921 Inf 3.399154e-05 -#> 13 beta2_ft 0.40831866 -Inf 0.40831814 Inf 3.723723e-05 -#> 14 beta2_ft 0.42839437 -Inf 0.42839450 Inf -2.700504e-05 -#> 15 beta2_ft 0.79997987 -Inf 0.79998255 Inf -1.674651e-04 -#> 16 beta2_ft 0.91170842 -Inf 0.91170568 Inf 1.970479e-04 -#> 17 L_epsilon2_z 0.49266062 -Inf 0.49266123 Inf -3.367341e-04 -#> 18 logkappa2 -4.30062172 -6.214608 -4.30062279 -3.565449 1.066371e-04 -#> 19 Epsilon_rho2_f 0.85035618 -0.990000 0.85035927 0.990000 -5.696325e-04 -#> 20 logSigmaM 0.10422369 -Inf 0.10422305 10.000000 8.907767e-05 +#> 1 beta1_ft -0.58930212 -Inf -0.58930072 Inf -2.544080e-04 +#> 2 beta2_ft 0.51556564 -Inf 0.51556558 Inf 5.379535e-06 +#> 3 beta2_ft 0.46458076 -Inf 0.46458132 Inf -1.711880e-05 +#> 4 beta2_ft 0.49055249 -Inf 0.49054914 Inf 2.249599e-04 +#> 5 beta2_ft 0.46815184 -Inf 0.46815469 Inf -2.146493e-04 +#> 6 beta2_ft 0.22994932 -Inf 0.22995166 Inf -1.963291e-04 +#> 7 beta2_ft 0.39666756 -Inf 0.39666648 Inf 7.352745e-05 +#> 8 beta2_ft 0.19322632 -Inf 0.19322487 Inf 1.263052e-04 +#> 9 beta2_ft -0.07072357 -Inf -0.07072451 Inf 8.061085e-05 +#> 10 beta2_ft 0.11680236 -Inf 0.11680346 Inf -8.984282e-05 +#> 11 beta2_ft 0.47100950 -Inf 0.47101097 Inf -1.175670e-04 +#> 12 beta2_ft 0.77135978 -Inf 0.77135916 Inf 3.317524e-05 +#> 13 beta2_ft 0.40831857 -Inf 0.40831806 Inf 3.594733e-05 +#> 14 beta2_ft 0.42839428 -Inf 0.42839443 Inf -2.828227e-05 +#> 15 beta2_ft 0.79997983 -Inf 0.79998249 Inf -1.675174e-04 +#> 16 beta2_ft 0.91170839 -Inf 0.91170565 Inf 1.978155e-04 +#> 17 L_epsilon2_z 0.49266060 -Inf 0.49266122 Inf -3.409419e-04 +#> 18 logkappa2 -4.30062178 -6.214608 -4.30062282 -3.565449 1.061017e-04 +#> 19 Epsilon_rho2_f 0.85035619 -0.990000 0.85035928 0.990000 -5.693192e-04 +#> 20 logSigmaM 0.10422369 -Inf 0.10422306 10.000000 8.657499e-05 #> #> $SD #> sdreport(.) result #> Estimate Std. Error -#> beta1_ft -0.58930073 0.05080467 -#> beta2_ft 0.51556564 0.14414517 -#> beta2_ft 0.46458138 0.17371753 -#> beta2_ft 0.49054919 0.19200615 -#> beta2_ft 0.46815476 0.20433941 -#> beta2_ft 0.22995172 0.21476798 -#> beta2_ft 0.39666650 0.21934844 -#> beta2_ft 0.19322490 0.22481460 -#> beta2_ft -0.07072448 0.22973074 -#> beta2_ft 0.11680345 0.23059098 -#> beta2_ft 0.47101097 0.22964742 -#> beta2_ft 0.77135921 0.22895295 -#> beta2_ft 0.40831814 0.23199191 -#> beta2_ft 0.42839450 0.23305156 -#> beta2_ft 0.79998255 0.23091833 -#> beta2_ft 0.91170568 0.23072494 -#> L_epsilon2_z 0.49266123 0.04727208 -#> logkappa2 -4.30062279 0.13652887 -#> Epsilon_rho2_f 0.85035927 0.03527714 -#> logSigmaM 0.10422305 0.07108160 -#> Maximum gradient component: 0.0005696325 +#> beta1_ft -0.58930072 0.05080467 +#> beta2_ft 0.51556558 0.14414517 +#> beta2_ft 0.46458132 0.17371754 +#> beta2_ft 0.49054914 0.19200615 +#> beta2_ft 0.46815469 0.20433942 +#> beta2_ft 0.22995166 0.21476799 +#> beta2_ft 0.39666648 0.21934845 +#> beta2_ft 0.19322487 0.22481460 +#> beta2_ft -0.07072451 0.22973075 +#> beta2_ft 0.11680346 0.23059099 +#> beta2_ft 0.47101097 0.22964743 +#> beta2_ft 0.77135916 0.22895295 +#> beta2_ft 0.40831806 0.23199191 +#> beta2_ft 0.42839443 0.23305157 +#> beta2_ft 0.79998249 0.23091834 +#> beta2_ft 0.91170565 0.23072495 +#> L_epsilon2_z 0.49266122 0.04727208 +#> logkappa2 -4.30062282 0.13652887 +#> Epsilon_rho2_f 0.85035928 0.03527714 +#> logSigmaM 0.10422306 0.07108160 +#> Maximum gradient component: 0.0005693192 #> #> $time_for_sdreport -#> Time difference of 4.135285 secs +#> Time difference of 6.984021 secs #> #> $time_for_run -#> Time difference of 19.02913 secs +#> Time difference of 33.26514 secs ``` Or with sdmTMB @@ -337,7 +353,6 @@ Or with sdmTMB ```r library(sdmTMB) -#> Warning: package 'sdmTMB' was built under R version 4.3.3 mesh = make_mesh(Data, c("x","y"), n_knots=n_x*n_y ) start_time = Sys.time() @@ -367,9 +382,9 @@ knitr::kable( cbind("run times (sec.)"=Times), digits=1) | | run times (sec.)| |:--------|----------------:| -|tinyVAST | 19.3| -|VAST | 22.6| -|sdmTMB | 18.1| +|tinyVAST | 32.4| +|VAST | 167.4| +|sdmTMB | 29.9| @@ -387,6 +402,7 @@ mydelta2 = tinyVAST( data = Data, delta_dsem = "logn -> logn, 1, rho"), family = delta_lognormal(type="poisson-link"), spatial_graph = mesh ) +#> Warning in deviance_explained(out): Problem detected: deviance explained should be between 0 and 1 mydelta #> Error in eval(expr, envir, enclos): object 'mydelta' not found @@ -471,49 +487,52 @@ out = tinyVAST( dsem = dsem, spatial_graph = mesh, family = tweedie() ) out -#> $call +#> Call: #> tinyVAST(formula = n ~ 0 + var, data = Data, dsem = dsem, family = tweedie(), #> spatial_graph = mesh) #> -#> $opt -#> $opt$par -#> alpha_j alpha_j beta_z beta_z beta_z beta_z beta_z log_sigma log_sigma log_kappa -#> -0.135067087 -0.050931909 0.537678546 0.471060516 -0.268664645 -0.080655519 0.355720082 -0.675048463 0.004059747 -0.566274472 -#> -#> $opt$objective -#> [1] 4300.721 +#> Run time: +#> Time difference of 3.236934 mins #> -#> $opt$convergence -#> [1] 0 +#> Family: +#> $obs #> -#> $opt$iterations -#> [1] 45 +#> Family: tweedie +#> Link function: log #> -#> $opt$evaluations -#> function gradient -#> 60 46 #> -#> $opt$message -#> [1] "relative convergence (4)" #> #> -#> $sdrep #> sdreport(.) result #> Estimate Std. Error -#> alpha_j -0.135067087 0.11243173 -#> alpha_j -0.050931909 0.09302442 -#> beta_z 0.537678546 0.06461591 +#> alpha_j -0.135067087 0.11243172 +#> alpha_j -0.050931911 0.09302442 +#> beta_z 0.537678543 0.06461591 #> beta_z 0.471060516 0.07549660 -#> beta_z -0.268664645 0.07225850 -#> beta_z -0.080655519 0.06157438 +#> beta_z -0.268664647 0.07225850 +#> beta_z -0.080655517 0.06157438 #> beta_z 0.355720082 0.01973953 -#> log_sigma -0.675048463 0.02959907 +#> log_sigma -0.675048462 0.02959907 #> log_sigma 0.004059747 0.04917005 -#> log_kappa -0.566274472 0.09138031 -#> Maximum gradient component: 0.002923767 +#> log_kappa -0.566274473 0.09138031 +#> Maximum gradient component: 0.002923743 +#> +#> Proportion deviance explained: +#> [1] 0.4630639 +#> +#> DSEM terms: +#> heads to from parameter start lag Estimate Std_Error z_value p_value +#> 1 1 d1 d1 1 1 0.53767854 0.06461591 8.321148 8.711554e-17 +#> 2 1 d2 d2 2 1 0.47106052 0.07549660 6.239493 4.389917e-10 +#> 3 1 d1 d2 3 1 -0.26866465 0.07225850 -3.718104 2.007234e-04 +#> 4 1 d2 d1 4 1 -0.08065552 0.06157438 -1.309888 1.902338e-01 +#> 5 2 d1 d1 5 0 0.35572008 0.01973953 18.020697 1.340448e-72 +#> 6 2 d2 d2 5 0 0.35572008 0.01973953 18.020697 1.340448e-72 #> -#> $run_time -#> Time difference of 2.327756 mins +#> Fixed terms: +#> Estimate Std_Error z_value p_value +#> vard1 -0.13506709 0.11243172 -1.2013254 0.2296250 +#> vard2 -0.05093191 0.09302442 -0.5475112 0.5840276 ``` The values for `beta_z` again correspond to the specified value for interaction-matrix `B` diff --git a/vignettes/web_only/condition.Rmd b/vignettes/web_only/condition.Rmd index 9677ac7..d77e99e 100644 --- a/vignettes/web_only/condition.Rmd +++ b/vignettes/web_only/condition.Rmd @@ -107,17 +107,17 @@ We can look at structural parameters using summary functions: ```r # spatial terms summary(fit, "sem") -#> heads to from parameter start Estimate Std_Error z_value p_value -#> 1 2 Biomass Biomass 1 1.4238785229 0.133047112 10.7020626 9.953750e-27 -#> 2 2 Condition Condition 2 -0.0331621242 0.004167438 -7.9574358 1.756412e-15 -#> 3 1 Condition Biomass 3 0.0000933144 0.004855041 0.0192201 9.846655e-01 +#> heads to from parameter start Estimate Std_Error z_value p_value +#> 1 2 Biomass Biomass 1 1.423877e+00 0.133046685 10.70208186 9.951682e-27 +#> 2 2 Condition Condition 2 -3.316273e-02 0.004167567 -7.95733630 1.757825e-15 +#> 3 1 Condition Biomass 3 9.292899e-05 0.004855089 0.01914053 9.847290e-01 # spatio-temporal terms summary(fit, "dsem") #> heads to from parameter start lag Estimate Std_Error z_value p_value -#> 1 2 Biomass Biomass 1 0 0.966429610 0.024274988 39.811745 0.000000e+00 -#> 2 2 Condition Condition 2 0 -0.040541568 0.002723387 -14.886454 4.036166e-50 -#> 3 1 Condition Biomass 3 0 0.008201795 0.003339302 2.456141 1.404382e-02 +#> 1 2 Biomass Biomass 1 0 0.966426764 0.024274856 39.811844 0.000000e+00 +#> 2 2 Condition Condition 2 0 -0.040541910 0.002723401 -14.886499 4.033460e-50 +#> 3 1 Condition Biomass 3 0 0.008201997 0.003339314 2.456192 1.404181e-02 ``` # Abundance-weighted expansion @@ -130,7 +130,7 @@ To explore output, we can plot output using the survey extent: region = condition_and_density$eastern_bering_sea # make extrapolation-grid -sf_grid = st_make_grid( region, cellsize=c(0.5,0.5) ) +sf_grid = st_make_grid( region, cellsize=c(0.2,0.2) ) sf_grid = st_intersection( sf_grid, region ) sf_grid = st_make_valid( sf_grid ) @@ -158,7 +158,7 @@ plot_grid = st_sf( sf_grid, "Condition.1982" = cond_1982, "Density.1982" = dens_1982 ) -plot( plot_grid ) +plot( plot_grid, border=NA ) ``` ![plot of chunk condition-maps](figure/condition-maps-1.png) diff --git a/vignettes/web_only/empirical_orthogonal_functions.Rmd b/vignettes/web_only/empirical_orthogonal_functions.Rmd index 1e119c4..f1eb01e 100644 --- a/vignettes/web_only/empirical_orthogonal_functions.Rmd +++ b/vignettes/web_only/empirical_orthogonal_functions.Rmd @@ -43,7 +43,6 @@ sf_pole = st_sfc( sf_pole, crs="+proj=longlat +datum=WGS84" ) sf_pole = st_transform( sf_pole, crs=st_crs(sf_ice) ) sf_pole = st_buffer( sf_pole, dist=3000 ) sf_ice = st_intersection( sf_ice, sf_pole ) -#> Warning: attribute variables are assumed to be spatially constant throughout all geometries Data = data.frame( st_drop_geometry(sf_ice), st_coordinates(sf_ice), diff --git a/vignettes/web_only/figure/EOF-1.png b/vignettes/web_only/figure/EOF-1.png index 32410e38fd74d840b1b45418d0b16452e98b6f70..2fd70ae43e9d6f6e7b475f218037ddd58ba192e0 100644 GIT binary patch literal 13029 zcmb7rcT`i)6K_Z$q4yFHkS<7s07?gu-b9K*=$96%(z{3^y$VPPND+`GAfXpQt~61K z0t!Jyn)ITe6s0}BzjNLn@6UJ6?%gvxJ3Hs@oo%0)n_z6DLr2X;4S_)D^mH{%ArLa+ zf0+_?QFERKA}@-+#)fybE=mYQ5dy(NAdL{nF9_rWqKJhkra`bRD&);9&@2w5t8=5Ued6TG;CU0T1#5XOxnzOHOE>r&teM}(t=fNNxMkcmJ8m@Ow0NC zn>TOzuaIWa6lc<~Gc6bCqWSr4hr;u)nDf}sbAs!6sm*!y-SZag zd8_XEu*~@<$@59E^O>~si=OtW&(AOVyO2J&F39J;cy%ftT}yumgogRQOg29I=mY}c zgXn4Ax*L+eT|iazZNQIYVYSV&dW}U%HaPac|6#O@-$S?et3n>*39;i6TH_gibr+8# z3D1AN+}NOT1KXd^6V>G8tIZef%#jm$x`qGvIIc>Et=#WaT$mc*)tDHMY}VM<#VY zWZM5$+_XV4AuEQ83eQ{)U(U(Y57Qc+ukVlRKMMQz#R;hnVEI*7pIq&Gg3yI?rEx>kA|(A2~9CA+{YiH*T7W zOe{CCvqQt1DS=Lj$j)E)NX{MkE%eqeW9mCPJB~hl0T>2)!ai=;zdm2>(j@nC2t~eDPf1hnAeSEdF7?{h1?5H{+==F->3XF2< z!4EP~8}g-wRYUnTmC{e4y>cqKsYz2VEdorl0m0HdGn);$7CG&>k(zsV-yzC|AwaKw zI?F~|h%L5(*f{39+mQ?tA$YIWIweT7vwYOfT}N#_46x~xD8;cI53QE|r;(TRAlQR# zhHPedySE}Mmz-3!|8y@cAh(lVsriv9Z-&GKbY^?>xtbTbSEE84)R&cLsJ)fRO&@1g zjk)`^jKRAm-zrNA^y6H>6>wW7NU!%k8=2;V-TA_foDcMHEDfA~FdBypsj&jym##)t zQ)2ek%s+(Adr;B3(n5kf*btE3DF{$u?m~{vWAM3@H^q!A3z!ZI#ekF&A3gCP&a3`H z0Txu#f4NH$(w{IAaasHb!}D@lKm(!^<;rB8?t6_6DZoxJ81MH;-o74UliPz2!g6n@hK$W-I&bp`P25VGFBNh)883 zIVL(*QMy_~+L(h9rR}tt;zI^%>gX)*Qxan=1>%!eU4nnEj7!~0-e@Rh$0qPemuJ90 z5x5*PlI=Dvao%gqC*)K{RYElt>#RA8v1_rj$VJlS2U1LRwsMRkJjI6=g{<$ zOG;htXe04Mv`q$GBi$*`&l)0BkjRR%(FlHe&Fhy3uab+@SzLbe(Cc0=-FSa&37`c} z)IHqdMhABNDk@okEW8EV`Z`^tI*|z&g6!BQ`AchMk!ky} zRTN$eH5@4p#8yioY6J5l(Us;NAX zP4Ec%NP{Y&A#W)jt_gL=n6paXF$_D|J|2TZm$)QM2Vl>hFZn;0Si;}`h#-tzXPh{q z4iR%!=YGP->7T*)H~SFD;~6l$#_&Tz(XlH#W9szOT0EOj`_7Bze(BpA<7GF{V%*oL zYJ%6F;E>_W9S@;Ao`KWGzHFJ&#Bi!D3*o%n`ueD?_-F82@Gz^amfE?w2y_)aYOX-ll#CQK^*v()-Pir9p!xN<=%A zinF;S=WVt+3lCybou(f`z#f<@@Sgt7$irn34sL00N(-YsB;0|kA;N|J9I~H%Dyw;=ZxRI?U;`xc;4`=fI-8L%?a{Z(3KlQQ5M^8(b*peP z%&C!P1HFJTX?0MWMTrT`z5eNSWOk2nxN5c~jV)hYlo(@VLoFI3`j~+D(yiw`Z4B6S zX6he-?3})v%WsN`89kpNF$oe)XX7S)t*iuDX~r@~T#8Hs9@NjE4^_bO5)vLs6wPbL;s6KC4E12kTAFy4Td;I;oAM*ZOCJCg20v~5p!CYd1h0j z{x{c>;=eFz(no5Nq&}aVmKVVzKQ$`!{=-nWoq=FFiWw=vX9;t^jfsx|&n6Q1Cpi+Y1#%D4$aV)v62z%2 z?I0wBAK}T9E&>4J-R7a4K7RiOnOcw^aEUAoR148Cx8@}rsy&Lhu$oYPK&vXuyOfML z7PBaCYn-!;0dB)ewZKH?X4Io6FG_T;51*CZE+!>gAiooFM3Qgf&6v&=xOrx|Fuz!Y zU>StqC2vr7cC%P_O5@+bb^td>Gfy35QaGiu79p%LWC`I8o@AZkK{BBd&6$olDc&Zz82|z?|?lHrbeDG()@Jo08jI4ft;}vWG z<3T`RUpVNwdqg<0a&*=3YX{ejSG@35{NGO_NBi`lM^S;(O^EWegq)TDq1K_Brk+5D zqjZ}@GeHqGy7~N6^Glq&ewIC|#r{z`JJu^XJ5l<*;>FA#qgUZw7-Ea#v|^$ul{Npm zFwRYoMn}_p2HXLish?Z9``4rv!mLozEy~kbYZ%1vF13fGmn|7TuZ4FtBeHPA?&^K~ z^d5RXt3^_oJT)^lqM_WMmUh-8SbX~FEipfm9*OI8NW;>r%N5 z+j#T7nw`os{ar4omlGaO>iqeXi zehWrP=xZKI=b{)^;`TPtp{7(V+cmb5w540Qaym`Mo0w4`D#L+Z@(bZe|J4Ix;?t%- zUOnJ{8VG}{h=M9VN4N|R%&w_uG%M~k|6(?DVwg}3iR zfw)h>bukCuMqGWR!pUjvL=79;b>BIF4}Imy%{kH55)Y)`_Qir~vHoa}G)~}V)L)eJ z(a=$RnzaVdvbtAhdomOn92G(a4CwAI&}Dyq+T4UL$X%$(GHqPK&6)u|nHvG1VXq|; z80~B^)UHs&2MQ2kWo+0S9|5wBT4`k|`1%PHbs*|NKSW($x_R*fsP6|2Xlwd`beNf+MP*k@T_L zmBQSikZ>nd9mwE>wp6h}P-}0T;|my(1F}ehJOb#FxQ0uRfd%ndC#RQvc+W@>I%quN zG`hnPGHXgxuYkT|Xc>$AFq9%|EU<>7hm%m-=(J$UIsR)VjsPQDPnr2Zt7qA3>++f! zI}T4aNxHfm?Aj_NV109g13DpC@!VOv^9BMg*N@LmWCIId`+7e{^_RTo0}zVC-a_bt z-XOE-94j(_l^A6N-)*Q9fhmo9E?v~Uk(ZTR&K6fA%e9&Az545KU+xC!d)8FQnw!}N7H&dTnsVYj9i zDa*L}&Y%NcwCE&+5ZXxBfUZrFs8^Tkq$S)EdCWY6CGZHrWDCMqs5z&@5IL#Jt+(!z%Rtl>ip|cL1gtp z=se@Q{Z=#+xy}2X$qAto^H5wF*v6EwLmSTb%Ai`nH?H-xZmvF$8&j>A%&r6Wp{n2G zZ`;Q?0DK6#tAGa+aQB;Ncz1_mxzBw#txFE4g45@5&yTD6{i3NbnSZflvQasd@e2bD z1IPpZ>2VU0~KD>i_<ir~<8&YPd#Hq|r5AmteYZUmYZY~oxW z*fZBL>TREl9uaznFNb*1%C^gZjzQ+t21V}bPfo6@{ZsFxz0#b}FT$X0nW4E0avU9xn6> zZ&O!(n}?n^q?#0CM)u`FUgxEsN$&j{^(G)ahC}{y@0x(uyc5k6?GKv9Po!Kb8^;L| zH{jUwwq%X&+uaP?Bzzq$x70L+BfA%lk8jsV#`8Pxa7rt(Z#d%kZ3B$P1rg>PPJ*91 zUuf!UqP=?x2X9&_0mDT4pDXLLA6iQKG#RSTeZ(0d%qR&Qt?Dh)caF9%ylN#r6%qT? zj?=nGin|_kSK^c3vHi;9B^e!gkULjh!O9jr)6+Bk6Gx({W1vWFirQYs#@i_-Kl2ft#@_MMscE3*gbh3je)bLxfn9Q>s z!5vOpS$ZVjz!cC{C?i++sDNI0kfmiV__ZUD!%vK|FMD~;A|uIU@Tg5`-0qj|~f-_({teZmNkbcmEPl`+-6EO>j7=m`a4tUL` z;LJ^AoSYo(hz5s$;_K^Iz@XIkqv8>wO_IB}nUSeKpBgwF>p*^YOKQMjLxUnzw1>`# zw+E9ht9;=i(Q9sYji*{3{PUyfox)@7_1|5!eSGZo>A1CV%2U&lX6BQ8X0=MY<}+); zlTg0?+oP+mVmxP84f@nMi<8=)+(0`!wR>I@yc?snQRiqY0$Irq@8U-j4P|p+;ugNn zJS#iS5B-cQ%OoWb=EXy-S^^%lx$H87$3K?maLsK_V9*3vucc1WR5Y2xgNg^2I`-06 zro`};HcCw_WCqxMIRP5OpToHxn(w%5wI=$RB$klp|i+p1f{LG|kz6Jh4N7Jhp^ znERgoO{Pb;xGnf(;n>#Ge@G;(uc0ir@XhaR+8x06HUsUBkHmnHn%no^X(pPqURRja zE$*y_HpAe-ZNnNe=v&HC`L3DdEp+TW^Z0v}(v--rWP|e6(uHx;FQEq`awW7RHJ+y( zM_a!|Y1wJY;iLl_Tsp1K_2?L9*^5m^N#H%Ap}L}x!qFf2L6Cfb5=t40st>ijE0mE@ zrj9eR_QV;+%!foc?>&aDHOicAEg z)bLL7b9Qrb)Q_g+a?W&*B*)I~eVe$2!^%{x4!{u!5fzbX zYY{^CP+cnAXkF~M4iiya3W{b75*^YCfh)G-a%K^)o^~p7+)v!ZR3=4=3$qtZSg$4h z%8+2PeUV6({i!THr}FFda<~(Jn=(h!DoXiTaTWqzSPI45&AB~)e`c5*)O^5g#EpTl zd>1S;5C=>zPDXe{cCWBkbuBY8^krpgSqCvI={eY7N}n9EXJVVjmau_%dl+1tJWz_L z9Ms6&sD=>Mi2gTpa|8e}H$GC4P4LyuJNg&b_s?WdFAk`|AN!r_L$uc#Vs$hPj%#{x zRAyIuk^`h1GM9a=X5@*C78%I#GB)&6iRwuD}p`oD_^ju`;E|l@_f&{OH<}@ZixCeSpgLD|F&b$;b18a)`fTaY;7q zE1!CK^8!NZTbV0?x7f>7nTRzK%IIth$3p7vRw$-YNdLjAY?B$Kf8MjM7-4t&&nJXU zg==w}kMZ}76LfWYW0>oSrognl{ol5}_}G7xvQn9|&7IO<8@UI*ERc+4+t|&7teW!6 zRCc_ttre_MMC6?2>3c?uc4Bp1KE4@(ovD3WVwbh-0L7%^Z{K(5Xw`bK<-;mvALVkA z(!*@>X_c#rhN0%AI4Vw@8uh6RKf!c$+FkNpKZ=?zd2eu!x<~fMLR*AJVAXY^+dVci z^O_-1pebRV{=@UQ);eZj`siUU-XTW5KWD&F_VT4D(bRSW@9Pm~^>5IT2@Tg7 z)ll~7)_cKRv4pmLiKB3nUHB3Z*)CBU)m^U?@E9BMO%!Q2{jYk$p(Q$1nBDxPLE4;2 zDdh^XKjdw}LTzcsU6}u^^W;>KSMhi)-Xw%b14vGTI;A5~HQ1OW} zUSu8j?c01LK4@Ir(9 z%0mKQd=ek?&6k+T$D3?2`5;>mEmR?Fo2E3$dU(B*x)18W^QJ#nC^SMALG^1G;b-d^ z(1G?%60Hn~enF6kLZM!kdV2~aKz`rIu5D?;nL+K+zwJElBG+&EEDTMy<)0xTe^A)J zU|C!28Kbg=H}x)F{r9Jr#tfvVM9RaK7G-9?vs-H#zN5B9ef*_PkCARob$&G?HqU9t`(ey*pJ)(v&#d|u1ie5 z1>6+RTP+QsT5Z?H#Y~I{5f?g!-5q``BCjP4|GpX`G4W@MH2+T2R={G(AD~I99yr-` zv{kgr(Zglg7Lc!ewSTDiaQv)vEtG5HF>U-RdlcK$@-R+os22*_V6mxzab2IIn(&SX%_0j;=SFMMymc7D3f`bdb;!&c?XH1aR<*{HPrc_>vqZA zqCUaT$}`YFF;lDvP^vjh=cVofUzd*Lk+tHbB|e}Bs`JxpKz0#vV$A<_>->L_!(>~% zmn9w@MjW0kcAqXRN}f$E)t$ZFS+}!2+4=b9toZow9fw~9Q=>nd9Tsx1zFvCSK4Sm$ zILSNmVc_oCI}ZEvUkkfWKmPkzQE^&yz87)W!!NC5{q3;({Ahi>U&HApjL4ha_v>s? z&)enWkMeI5ky(*KRUglfy(RwDoo#W{;!bzN|4#GO{!Y)u?VH2#wz{DL44oOV2d^U| z@VG7y%B|n$4>2M^(lnd-QV!2N+#K#)K0Wqufc@xoGK*6HR8)WKm|N6I6rQeEH4=jY zF;jCpV-?;1R3pDd7a}^|QV5CQ0%~^Z>??SWlPF{*i@5-IM+Q0)k$IOK{4U=5i!-cTk_l6lpJ~Cc z_50=r^kUrkhc6B;AR)-p+@mg+7gL2V$j(_a%r6X9q{h#R`p#VD^?Z7IQb7H@a^%Ws z${E+U#l_}~KR5h(qvuwdiI0Vyj3~G2w5^d88YD*HIL3UALuXX!1q$b7*L)i=1 z*%eYWk0W%-u}tN0Ov6`i#DJ@K2#viS z`w>_6Ae*x1+;8c=tw3?2IwySit2rr@8HtQ@jEpFA5t+MLrvzsY=DPB9@e*wlad6&t zAh$KLAi>t4-LFT$Gg-3km=hx2+3N#vOFjGTV78?X-Msngl_RJu8n>C)-cJiM!2RX+}D^ekqSsQCB>a9X0WF{ zuen7pJbYqCzcLi?)0lTvxU5+RQS>XOaP%M$0Tx&$v||B~6uF{muxA^4t3XQvXEDW)JSFuVH{z^XpQ}A(hld z1<~m@D&J8D$bqI&xGTc4FyiB^OC{_<;o)~3)7^oNvWGEC>8HZtXg)J&sjv+fp^B znHTq-(u#!GZe_TR7jk1{OV#aI36ktxad))re3etl2d?oeF_P#gS;DK{zt;WA!jV9+ z$w1K-z*ZhQJ3AZmN#5}quHIk}4nEMw>7gd!7hw}jvH!&a*vQ5&>!t|Hst_|CE`HGP z1cJ>Wja`yl^;_QkCdq2%&%>qAE_QOaSUrWLA_Dce%%C)=3o8`Vz?YdbrIk>!LS*_E zZZone_yb=$q-Pmd)~}p2>NmrPd-r z(yvEF-S_A*9>UzhE+8%v+0>xJmbgk|m>mb_V7#l?qgLuHn`hC+Y=k}p*#mNID6T}C zq)YYX|IASFXC#K%Rlu{Wd-d^Jt$zyN_LWnZe;n#qy>L*6gkp;?jX&2!{bHtjnLma9 z*7NHJ1&pQ#-v*@#Wpab(Q&~4d+3R}sakI}crRz6)r3s`=FtH(?!IzeN5;PcIw}WiK z6=D~XrC$f98B*3{b0EcBtbq#zL`YfnlBY5Qz%ytvAu70tgJw^vz_wP%U`vS_ za$E2(_TBw`yyJaB=zFLR2Sv5i;6;g#)!dk9NhMY<6+6naURs4VTS?sigABp+l=D$4z$LbZi80;XB> zpucw>Y5A(yrXX#79 zA9d>AzFS%>rAYVcVs2I_Z5=V(2?s&Va6n!3N;;wi0UA0-?t*@L*U5vZP{X~@+b)BG zLz|@m3f~D2>LVaqrm)7zz>1fOE{U~xm)ADWFI<^2b`Ks4?0&uR^rzZd+TIZXz|vm= zoEM8Og+?|{J?iR63Xk;#k;dR7lM4h!w)_d`qgd*G1C%QJ51bb0MNG%X9t*{)ZoXo;HkUCicGfV zfDBAy@PQ{TJ-u5D$W8YS_-KlM?~@77R_mvP9y~ix;wElFNLibcokl4!No$V#F`?tE zq>aE*kE^lCdY3DDe}MpHbq3 z6NE)1DCFEY*U_P*$l2Ed_0XJU7shie{(U*2x^HmL!G+#QGe?pKiBU5Tx~JD#N{>9) zq$37^Hw%&z3wzFZmL=Ue;fme_f2 z)G&mlgbIBNC)H;TbpF8o{h-N_P?E~w;#xp|_y9^1N5<~TKITFVpj&x}Q|nqmYf^dS z8L?o^>z39|R^s;oIP^IFs{>T_S%C`7hDtnXP5KIvwSGQrPD%_dS7;^YP?jV5@2Z!) z6rJG=$0D0ViX?NJALvv5Pyt+X+fbVBI^St!5?R5`Ci+}k1qVgQvT5$Vz~ws0$+=S3 zd^DfZkkr}Cllt5s!Yqj2ze#bclskYt>NVX@HrbZ;Aam zctHU%U`&EJJ-A7Qp1oW?@JZTZ~DeAu8~h>+q;2`5ph zxMSmfa2H&}nMVYFDoaj~9+C67{#mtW=nb=&MN%0@n_)ZWG!S0GA2bT428Th43!Xej zw2A@E*Ic~-cFf8B%h1ox3cv~Na00GfFlZ!SwB8J&j-g4raXd8?D)r<96y~kO^zDs;e>6P*VOO)jjh~)b_$8EWC#}8 z^&{R|Ybzhe)8$>cKMQItb}D6cbH40q@3rsD)uCD>E;LY2M2VoEojjqo=};$ zRk-Q&{10#L77=)WIxoZDVXi}c%EUKnon~#(gB|e$EsicF5l1D8UaV-2_{Z=S*$1pr zTAQ@E+`<)WG?S^8(vL_L$)0?Sx*;cN*DZ{seGkZw5#FzbYDZ+UqL1Oi%&%{XmoW=3%VY+r7=xL=&}GA|ALhI?ke5vzUcIMx?0VFuY;1Rwl4wtZ~;D*^YuzBSy- zo}F<2X$*&g(E~vU*1j2Yl?hljLQtBp{-S=!>g0$*qr+Z8kZ-zFstaOlAGiLb;@M%G ztiSn^9_fdti69A#cl?Hve7mPQ6R)ALj`!$Rsc`NSukjVDxWIrfJ7PH7X@`BCZBc3m zI9-}Wx#;;)!Tp@KiX;4)xB6x_J&8@9eU4vx$v+`4AJQB)uy{a;1OxN|b%yH0`)(JG zhHzN^!xR*iVOGM)uPxlY&$V@K@^@gvqpafa$+0OeG=oZzO-Tc}jM%rDZlcAFyjPQ6 zO4}*^d-Vk3nK!nq@AYZcr7SBiD>|;x*%7*2L}8{Q@Z&@VBEpA0F1*w`f-2%);7NaS z89k9+a~8d1chA1f#v9Ne^z$Ahx&fd~O2sqAB`F^t(6cr4)$DR<_!xkjx z3)2WfD)i>d-9SHC`_MIObCBw$9pvXph3*1Bzp3e>txNf(XaMZ!vXbbmJz8wFK-ffQ zrfX$*5*^w%>(OP>GcKU-N$asWptjd}@Stl?36UljPl-PMzy{|iaEcas^dw&5Hco0p zmbAh{OB-N!hsN-uHeBh|Smo^2{>~}SROqpK8%joGBwm6A>aWAe{eoE>0>6{{T$j~P zn42tHxQ^^Vo54$*EITez8aP{H#%C_%+#8`;~{cT26$)zUr!@SQ}2^&7}sP8)u z){p;kC-MQW~e_zirSRMAXw&nV(GT72^T%n8@vFp8$ zEmzMk<&l5a>y_jh?ptaJY;vZOPmV2+m*N7|y=~P2^LgEQHe|S#82Vw48sPJc0&XGH z{pw;X=lq(zauLZ{^{H=O38KJA;+xYRMoUVq3U4@IzP;kN^DsQTol2d{8 z`IWM(jZa#&oFa6{3{p>OGv$fMIxYa8B?~K8%N^|e>hEGtV|`D63uqQztn7cQP0drg zVA6*Rb(E9yhm;=M9N6Bn-ec<5cP#O0p1o-Va`UzO4*J7f$R!`;_&PSbiOlwM^K2|{3;j5X8{FLN|7f+bw70o zN7W>z?_2XuP!Aj2jU3BeL|^ALjS=47l>K2%Ijmv_deQ?LAtMy1279^Z3XV6j4&i$Jp9|TtYeYl+6nvnC#k@=Cf!1%APS%E*6Q3$K{R@}`W#)O7 znoR+F`cCa8VI^iYU&zN@j6Q&vFt+NewDcb@9UN?Cfhx!xxZ>)kh-6yCWDMvhb{kUj^i?6((WO6{ zShss$E&Q$`zfT|kF@00`VI3ojZeMl7v61=!J1O9?szvwi!pHZ_(3pSOgVQjlpmh5- zPE(Vl>O}Lm8>K#T)vi1j+TeAgBG;*rDWYvrr$yG~MM6~7FAwyEba=DXYd*K?>RXIZ zxzGa(ys3nMpW?WCC8o>;|5%kRD+7DeyJy>BEN|)05Tkc>O)mRVO8l*6+0O zo@xUA*;g(5#;`zo$8WJ8%7Skq(!Be%sp$BJ|s%fh)=cZ#=`umK1rHzDyCh} z`2C|(-BCFsb9pQdl%ZQ^VP8|#AJ|(Ix!dQ=MXN=E4+Qfp2LE5TsmGGTFd_K!xx~wFdi2EiI1}cy;omHfIdH3g4r3i|hCO+es~n zLa?huMexQAwrd=4eljZv=8+X6GCN4r5Q#&BRJ&KkBfnfZ%%|vhk|7IzV9teOylUL} z(W4qaD#{lSG6OuEFI7gE8%(nqYStFz>kP&z#Z5G4%_kC-M>9b#Arr!&5cuJh5!9eJ z@lWXSlaf_;cnvdeWFVrYR24{s0LLadfP0PM#LuRv) zdYma*(3bdXhUNgxQ#R7o`!4?v{r|EtQxO-S>Eg`g*Bl)75p%^K8)Vwe#u8>Ly3dMH zdf_6GFEiLMQ(cjtUB?>_=gxG`vhu53=|t{qeINP~qMz~ERXgK)19W0ea!(|A%kG?Ahz=|Nc6g?@rwOa!)4TaL=eFs4YQxFny+a$*I4M_Ehn%LzSGZn_ z>5k{WzvAvpQQp~*lypVdP9S`R!ijE8^S+$0-dCA+p_KnGnEd}XK1Y0YbnCH`zQgq2 P?<#s)Mw+kGfye&`hbX!T literal 13110 zcma*Oc{J2t{6GGh8OFYcktIvEY(phWRJJTb$TIeqZICT%q%u?XE&B{Y5fL$n>|XXP zC6py3St2_jiTQbd&iVa*=kxvRd(OSDbDrxx_r6}Qd(ZRneBM}Nqnk|hXXya|VA9jo zG6et_`M*w!IOR~1A=D{BZHzI~J|zJ_1psgWPy+y40I&}L{{R&nppp#W>O=vG3_wu= zD5n}KfQkx$JMAhe$tpOU3hwk-r;M4pbX0wn*Bmt-7}j7v^Vu1l_)N}i%*((jdVFW2(_sKWtuIF-8O(}t@% z6;DmoQ7Emgtpie@r;=5sl5taYr|r~&vUs)g_@^kvu$p4?lH!<0@kycFkEBF}Qj&cs zBsWT#8Kv5gQir3w)1`cpql`*Z#w00I$&^!nd-N2_sn=5~N%%gt*APDGb-7pJzww}w7KdtEhz1M!IV$owpw88!UGcVZ&)gV=3^vw~mdp#2O z1U~3zdP!NlP5pCxO+V(I`HL+NE7zTjf(O=KgiR5nc)yMNle@ zX`F2hg``WYD9(^{ni?s-#qhecQ~!;*iCZmAw)DNB{oTL<1AMHM7T) z^nbP+b+_}yMqP2>U-1cF6iE3Kej{UI^W@LoeLS1p%c}wZ3C9%I1~v(){v7WdcfQw{ z1M@*EC(kJ<3R$uP)>zuL!wdNB#(#5zR~*w;s=avYEEt@FLLcpCdA4sPM(*W(xpt2- z8x~7D(f<`)VBlb)uJd#HiDLQq@@U&cs{t+MzD0NwJ@`y^Xb3ZII$+lL%=G2_`eflP z%J=N`paGZ0L28U{3gz(!y+!pEct#ctiDF%aPbwRC4F*p<6N`6$O&C3m49wz0{k1-H z<6;tr*CosDLH+TFhwHX~+4wCh2NIyIaoLyizk;%=^tBwC9CUksb0m$rUU?;J0JG8W z{jzGx&g(Bnhp15u76EpgXQpeqxtDQ(O#A1i22p*{hEAzc+i0qBKXWj zMc=|8NyJi4y~h_c%Vhylg-S_guL!b-U74QVtDkk--A)3X`Thp2Mv&=UPU!PuJJZuR z6xhzDDv+=KkiPl7$CJO^P)SLWKXdmSUGj>dw=`L|qfF)zZ327PFKV_NMk|7Bxq@Qw zG7pu+`Ik0Q$&+S{pMm-G&qZ>^;+WVU2&+hFVef3V=6G==Z)tp@*81IRgJ)A}%28*P$3)V3mahX05ubr7N% zWWPPLQHQvwolTaFrf6kZ?T(PkC0z1Imk6v+#o5)ZS9&3t89bX;4W>wY2r!cBHIxaIErJa%9 zF`NqBC!TPzI^WVDUK#hRmsPz`K`=$`Tf>gD<8U~ywOEfXyv?L-Ry+UCl1vlu)-IIP zD~j(9iiL-z(nKL~UzbgJ4HkUwM-1EPmeZml#;kki8tdz#(H9L#ovv9ek5>M=&@jG< zm0j1o4G+490B;pSdDWSd*xvR=d}j37ENU^mhs)h;t0^le+CX({wlcti{Pl*rzpJ()B93s+yXJ5R}}OzmS6Ib)(clw?{B)`z{4Elae-RLbFA|7GpnnQue;>R)^cUC zC82(+f|2g~aW-OkvzT6er`8E?Y`+pe>%^^Vbo?d1)L3W+X`)sGO98Hfh!#d*$&iIg zW18npVEHFkfIQAj+Eu-NzD=K~0!}u)!mBY6j!oz_*FE@a6dSMG%WfT=4T%iv*awW}OHze=RL7CnVQ9 zh`8;i$+Sizo}8mnE7tisR`z|X?dJn=F49VpA9P+pNV7vEs5EwrTUQm_m{#GJ3iO~^ z1D*ef1eT&8Y1Zp>iCOOMkJ{j4^DBDCCDh>UuAY*f6>Pib&sEKNUnN^<@_3~X+2$aW z=T+g_62Qs-Q$vsm9l1K2fxXl7K6?erg2^~Bn0R#jI4i5f=?pekmIEf zQr*ESW82k1yj`xV*ArI7fi0q?kZXdj_lyAp4@?&X?6x~ z!N^mh1kRa*33hDJlu6k@g(%q?*LlhrtNugZhrck4K04saaMHAwd$JWJ64AXN+3NMd z6HTsevSv>m;nv@QTzakU?@?gdjRGU@e2eC3?VdHY(-dB<7-uC?k3E~X4!U01TiT1{ zq)iBkEf#r0ki3}!9L+^v;B4>ctSHc-BUZ1p zwL?m_;MFx6-dn$Cezt17otK%9`2t7q#7NLmcSaGWg<+fpQy2A#l6d|p#A-_40zm*9 zBYh64{Ny?=K++I9PiD2c0kkQjXV5&hZNpAq$x3n#k|)DSRsY_S^aKDuq}RrRN3kG2P~4Oa z2Sd(Tltw!Fod>?#|It9{oX|2Ws@aI-6F^9ECLI2{d!}3%RQ)MY{5aN-K!E=Q8HI4s zz0dP_TI2cz4Y^lW~bH5z)1R`tF?Am8)& z4J8JKg!0H329@@{bB{U#pHXAE3Sum$bfei|9N|+4)g34;vK z-f`q2{vE`$4Z?Td=lzgo;VpYNn!6(Ufshp_d3@+f$zRGS=}|IW{?3smb@Opsd%x60 zxClqe#V%}|5c%5QNfKC}l{d>riaY)*!^YDSy7C!vSG=ca6O=c71OG{BnQ zo!2LY9=vsOVXf}+@ooRg{|;UuCKzESHS@4&W(G;Bjk6gt2o>anm<3n4@5{leTzaOc zpzp;>Bi4QZI{!e-XWLX7hhl`Q_7JyG)EFXp#_V!|HS9S(>1+@)JSCl^`zB+Q71bju zEj0{~`u|2k`^(}Qk0hD!DfA{KqO7^G=CCHF@zyhECa*%DJ(uc|FTh^4(XiFYl#=H~FHSAic5o6Y{Azx96q z#ciFXFuz6@FjM1aRR89e?+aAeM7f7%od% zukdFG(#D#;Du)W`WpW}-aWXX&9Gyj2?^Y+Hk2k-$6pDIF(p9|DV_tpi#bJ_0rZwnkCxOJEy>$V8g2RL_}mXOt?X(usbXs~3OgjFsO^MC z&y-A@Ye+hg1gFDF^|L}jOyaVGsUmUIk>jJFoT_8HqOyn68r+F%;x;p?#e|R)2~h#K z=)Fa^*(rXbLD+na?_IRNB53rlGTe8ytB+87r_=q_7ucwX@J?GF!BAjj!rAFLRibw+ z`h>_vbltSH95k!=fPe8rK};+n5Az5r#Im{k6Ue6FrdE7Hq$RxgSlv_m&eui-Iq@I9$^`_*;ca;wCVR{U^EYfgfWCr&X5QqLEmT@C`1WGrYy|d9N0E zQ^vuHxKD~p7CJX9GRs-jSW>8Qbz)ReM}zLuojV+|(td&Anp)IGdJ&{46Q)6p(Sljh zGJn`yA!kR|*~nT2?2b@HiRiHABWDGWD}Awc{04|vD`|X7MIxV{f^k&2jL0|sF~Sa& z9bA3pOwqE*LjQ(1WeBnE3LcsuSe+S4Q*-5_oA=EHUMV`5KU;YoZx7jHlwBZ*pAP=!onkAG)Wz=yTy6@r+>wvJ3| znjSg<4f-ZGcHfUMQwydLJEbMzVdFC+ku*Iy?QRbYlv7FakE{}zMwU&`=9eLS;bZSqzzgTmskg&eBI}zF>a5(c6<67 zPK;quHAd+)#qWfXI^L4>gYYD&QK;j3>7?So;f@e`U^FOSgZK8kuQW7KSrbYO0V&rU z6Z;hH`Zoo+L!xIkVzxB8l+prinqv*A2ARy;TZtU;BlSY=YB_f=JVzg2rN2miH<3`N z4_34*Jd_g~eD$Q`^-bgkTPeq}P%(3mMgNqg^Sf`i*l6tBkGYk11?My!Xnk+st*mV$ z^+ICbKP!}PbidtRn@^35Jv4f|s%;ayy~F9uM^ck3#YlnoJ1t}H?>xTtfl6*-HOwB) zfryrnix~VS@@5_c@65+)YR{fw;MehQfp<+>zSN=G6&qAr!eX^Q$zui#GnP$wEz?QM zj(g$L^i(EJFLG|Rsuj;)x5Ermdo=pzZwaWk)tBteg-b<>WX^TxU~e-ZID~GbG@8is zGC#V53zdh5D1+`>7lEbL_N$1|q{M2pT1>$;g#6U6r`x>F@$GlSPu{gfLg`W|NTd%` z>+CHvcsy}K9`5uuFjNTIZ}mukiQG#5 zY<&tyC-vM%c8m8w7rv-4qkT^v!Y z?m@9vrKx)m@E&R;rEoE{(?6yNFFzk&F7U^F7H|3_<@L$`y%H!?{v_A zV`#t!4)Pod(=qlQ=>3kktYSiw@VlCpoUoP(ZRr+IL^jv2mF`1_^6AJfCP%foq4}Nb{Dsb?ksi!~p`H2H z4b%?5xb9!hiJUXEzSuGJiiOeFiE%6F_Y*v)MNAD)6Sm*gQ4R%5$6P4Nr)v{TDr4nE z)Px0hbvy&4qerY5_WM64(7bEH5>ue>vOk?ad6~vzwAt9m3;mT2w9WTbhCCfUN+B3G z{3D~6l%Bkx+|4s#pAfd(?R94Bx3r{=x&=!`v zbw0TzW9L`ww(rZxdN^1r=Bk4OD22Cl19he4kdlENp=C;~!w4Tq(VG_fHNP zBzt_z8*1E4^7S9?VBjV*<3f%3IhcjMH0yf?Hhr7vQ!QY>L<`cBe>a!kKD<`PFY`7?-~(aUG?{T88ARktZ1NG+*_l zB@0*yp##~FV;S)HNmRmxcEz>pcvmyJv9pOVP3{-oTjdL(o_IQv;d`UJY>NdUE`b~F-H&q|)q>cd-3 zw?0tiq3&757Bjfp=*Ft3Dr0Nn5h4SjFA3fJVdVy4;Eg$Is9Hyni!__`ONkL{90DV| zyiw}vkU+N?Msa4tM7ksZvjr1QLhSG_bX{r=XTyqF85OP&ibAVk><* zVR?j=w7kFZCY=hN=b3tcv05<6Ac8nUQ?csTK$EbFEoazmy5~ygSM+4hRo})A68|67jR{=yxfXeJN;%L2-sY8xvxd?#d93 zy3m7z2H83WSkf|J97oQ);PbC>HI266#MEDK8-axbW z)nA12_y{i%Lz$mjjk%dg67!Lz{MNeANkrnL3V8u@kc4{PD z9Hy=K2XZPC2$$jiaRQ8gSlQ|0MVTM11<8XiUQ>Hv>R?9Yw3dDj?08gjYq<`MkA?vWaxAV8OsXvs4~8md}we43PEth zQnCDHO;{bs42zU~n;mbHAX`Y#&dZeKGq`RZBUIh8`bPQ2s+528_xV4Y zi5Dh~xZMQ?pG7_~_zS!Ymf+KKppbb>2C(3wN=LT56ZrXo7FwS5oH)xun^L9o8`9JxYxZwVSS#C?1%=IvM>bOM{*>_};=drw+ZO ztVg_yUh6db#XgulXc!&Y+#l%Jg1zf+05)fS?ka#fR&cMUG9F&~)X<+9QY1Tfbv!V) zz7wJkVlyL#5>i&3-%2YCm3{96y%l~~-3glPHXIDMOMizYBtI}~rlaGJflU^GJ((&t zKEIFu{s#j!uF+B7l3#|=JiEWsIJhp_7?$qA(Q>0)jA!SPOyJ6n+uJX$bl`&!AADDv z^QX~rQ=;hWE)_O9(P*_nx@BSHYPC8zqXr*i8h4c-8ye~Uzj-1=ox&vk&HtA+`oAKV zfc=vLx3;6Yu&{mZwrfYel)qaw6pw?zp`k|yZu#Nk$5Wjb8%AR$8|dUN0!M%Lo^0la z@2^kQ$hCz(+W&rZs&a5LS9t4h{!@no%0HrSlN)o&n=?$-?hhAPKp8q32ny5~KE; z!XHK9)qH8!Mk%chVs)g*ppPW7;pMmF%5(pwaE>Ai$C{&2%AlJ1?o+RdynMBj(eaNR zUuzxhC_yahCtuzj?v9-gxF5@mJ*YciI42do<$H~?F?z)<|JW`ZiutD8>HOw)=34X~ z?X}}h;w<;tbXT7(vy1faCQvW2) z@^b+F?ym7|aAvpnyHLkDPt2nk$e1T?y77iL!8Pz{pBj~gdIKV34pXiO63!qxUF)gP z>??mV;8I*>QI-aPb5QuU@ZKk0vgr9YVW|_x2*nDWXZ#a3i(KkDw|=iuSzKV)`9*m6 zH;hN;0mY#~v78Y#>iePpljjmpheUAX#*Xe9-4r4li258IYkF7bx0=At)E<^l3fhnk ze)}4Ff6q2JcK8O^VB1sW<6%c-cNhwxz~zPc@rRd>_IEeGQHs7EacQ@~!RgDs=K^PH zdjAaQ3zTd<5H3R2NmOEZGtXxt_V-i%9sNs@{dbEJk)Q-;bSPBM=OXX!x`r;Qe4NMk z{axGL3^Zwbc>YR}UlCz@TtY}HO)(7*}Q4>oAfcATb8XHW(&uT(K; z4%*kjq*C=6oY^%HhzrM1Swu$Ng=d^rTYNeu0?5` z+!B0@d!iD+4t9W9yImE}v|ek|l+#tgXJBM4i4c`60n(Xoi@>rDpn5Rl6V&$r@+x)a+<=Hu3gO}7wf0~8 z@CwOw@arZlzUTF4Z9}{$mqiNX;vg+laaB*DJ>H3UHDHqm?;3w$L7Q+Gxi47=0Bi&Q z)dJiZ4*O&Fne(1GR4?js6W&DS=p@(KpVV@&VD8lZR}&^X6X9WBY2Tvp2~+7s53aFS z!_d!Gbh^`Db-6ylCYi}pg`TF6Y#%N&u4S-6YhNuBu@$~r;W|v@gpF<}J=6|qQm*gx z_g8AzXRTEl_A;7c7e4am3vkmj1UOi%(#%*OEY9MEmisZ5NFtNk{@w zJG31}!JJhqgp3Tc?D3o{zO?e+!=Z_-y`Yn_5yG7!3P{C(KtuWS*$?H<3 z=laNPTfEhBl-YO;wSB>jj9ryv8&qDi59sxD@-iqS2mf#+H2wkKJw;{*i4tS(JwxOK zF+~}2fyNI_J`utmZW4~QJ%EG6_wtNx>60ymE~cPsq?QODcPkCf?>Hf%*UdTca!WOeE13}MfATO3R>NWW{u!&E|ays;j7`sG!Rkc~akFY1P2aihpm>}df z{O$9blnZ_Fipw1Big=lB(|7o z0H%3#ygdF8#wAZa3#wJGO0{!@?G*!9>jOseKp0a(;o^vo@zI{F*gi>2j4)%xQMkgt zLzDdQW%)+MLa39A+I~UEa*WW|SQDaDA2|yap zG&{*pg}v3u;vfae4&@&(kA-hllh0SOy;XmJh9dmfVQ9xT?dj&~m*190>&cId?kR1Z z%d3O&BFlfW0=u(62(gy;4q(^~E<-|=GUICrHR*4YpRSrVtLqZy6M48aL)$-ARjrt* zT^v2vG$-7xJd*}hYswn;z<5*TuEa+65#$E?=L@c;zE{5vinVp_?2NvstaYet?N-xZ zLL-MHfu`Rpb7jb2p`-a5CG!z@*XsKdWfEiqsrxl*XwDyZOA8)(DMQ*R1evf;5yj0! zTE9{UdzE+@r!SB2T-yghhYA8-~6 z=MA96U68pPUxbZV$pgtDBred9bm~K&fi}3$cPuh6hpGa!}FfXJGrg=1$bGD+0-toS|txvW`MB9&|C1^Jht3rVgUj15Yo*xY9ozBX#We zp5vKBqV8xT1}$Sh8Z5*Z!lmagohOt&`E?7IMHnTc5ttY;s0)cX4auL&n{bImZ5OPS zkIsEF1n1C@)1Yu82F93Khu9ZxdC#n-<4brciC(R+v$3bg`Z|jc9$r626X1JitbTn^ zWwTqN$`t0}+ffxm;)*F#f69KdYtKd`iWd_kk2;z=ing`SSTdtz(;Ya-o~ee%`plb@ ze-#~Fd}tey5w5aw&iG|hK$^Qm7k05%R_+!T)}h=#!aj>QIA#vL9hGO#LX1()P$l|b zcM{fbI@$f?Hw15a^p1W2>37Y{KN zM$ZJW{VOqYcMr>`qf)0f`Oq&yTR5JS^%O%WWoL#Wb}V*+3!9FraoGIP+L}wGHN+TZ z2P(H|6b1cxU3DfKF4o&ItyOsGjYO4f8&pK32gN=aB-NXD88J-ozu^eds7M}`l#C*( zY`rM==R;!9l{bT+Y6jc;rSRaDxY=_`Cg7EOa8a3=2Hct(41*Cc4UfZ@&f#dJ9Yl3v zIHm8doYrDZZz#XkBE%eOe_YJGDBnF`ELByR@BOls3wuA8S>k#6PV?!mP6x6%#RF;E zn#r!k;wSH~WzX^2$^Uw?IX%lq{x&bz0Gm$^5}6|R;ZriWvyS?6*T%2#+sktm6M5Mc z&Xg1OPO&-buW$~wHj&TSvO)~q#M>O8+`o=EB!B3|c7F&X_R8MCjC#ZSaf=HMf9$o| zb(P$z@u27~T}lvJqW5-C!NWfViskeq2CW;20A6AosUV0QTF~81%)YiErt#W6`rmB% zD>tt}5d+?dnOq)6sY~P-tw0Z0!e0qVaGy#4d~h2dI*s>x=*-HvvaY4c@hfR2ht5?0 zKE5kga4R_@VhK;hbON*BSZjxoeppZI1|xb(oI%e2*|Y_)RzTKC z;I+QewE_syS3m`BIZu1Ku%KYrYXEg}Kg_~_L6ydLmIJr{G-z$}Gr-sYMXGx-dCY_& z);9iH?_Q`wf9pa9Jrf6t4L8~I@(^Rg=1s3)OPMV7qe%tEp$M)#?)69Yk1lO#@M-X9 zHE(+w#QtfD3fR9)p0}mLjK`9QKGU+pD*~hbP=u^4-Qd9s8MN+N;$5Q?)h!~8*l(#L z*`2P>Pov(EVf;Tm@OFPAVC4J<8X69spgIONm^`v7C-DIOCjW)F9(p|I*LJRz8Nacju$2KNnQrL)p*Rc^55vs{IM>^)f-w>IXzc#svp%JipF}iys_wnO zlUi<=^eUHb3fg192pXNCYjhyDKaXn-Io4ioyQ3vXm+hk5lJShk30EA^3Vyt#U#KKr znQK^QYpYs(E;r!{aU~~17#uvb0*QCqDAN-i>3XFs+&)pcHS#plB2=Us#yzX&@i*oN z;qgqTwR>wCq$_NaKR+%3Ibfu9^}&xH!cAsR!n&b7&c_Pbe$qovRsA!or?u8wK2PVod9e& z=zC267ji*slVl0oIOtfqUrLKTcr3du&yJd2q-SwUvu_|CDro0-70phiFa@01wPwV6 z9p_v{PY4(xFdP9Et`6I{x+tK@MveTqB<)(ZB3Tyuy99<#9F%5y?GP<2o4&bm(-g_M zpXV)N=fO<>JU%lOT1YB`u3| zw?X~`$-=)3u-jrfrZkZs^v*c;Dt|CYz%B?)vapd{`J6O6-L9TKQSx%@(VV!gVzLk-bsCG(FlC?%fBn51{TjZ%nu^<5c9T}&S$jaKok9GH2{kCybhp#UkEGAAgV{e9yWA8Cj36X2492{n@W zijozKZlsk8wgY>t45PQ;0*ddJd}}|$7-Vxv(tM0Ubd z9@@xAQR9*%z&t#&cyU~o^@QDj@VJ8KiG_g;n59~a6jQlGJ=p5rN*YhB=c#5#HJH=n zF{NPM{AUUDuIgy?%Ew zh6(#ytRdK@Ui?K^Vo6o1iRDKc3|W;Qt%q?t>AzRJpikO)tw1)1t!jV`sjogl#Dvhy4zPxb)UOZ5tgyMJuk1o zg7$IHEvKBObOSQevC%-H2%goc!mTM;8#0vA9L%;VOw?5EpNIfyX>K3Pw|;v)l_RJ(T1 z%iGP+6Nu@7mJf^7?!EfAa0~yOc4)2bG?Kyl?Wd_B)L0Jf1tyX&her?SA)HOcl)yno z1tcmu{#oWFxxMor1NJPI)RjQQ*gt;RFA(zl&VP`zUufQ3(c`(gd)#~V zE$Y8&O_Bi(7CSf}b1m)U9fwzu-RX?p!YqQ*N>9t93>gTKqaV!PNgcVj7NXP#E*ZSn zZ(xDqDxU7Di06{o&DX|9E{!PTFV}0&8H&Lzf*-PD&sw+)=J+}@1Ugt-z{(RWQ-T`8 zYFMnEAB~vWmdRp8qx{OMUi#nSRt`g00MfkQ%Tg~Rq^|ne%T;+&mk*X&v5}olA!5_g zL>U)ER1u6k$Vg(B4AQ&;ez7lT()*yva3#)q|3S)>7lZ`LLYE#qrXUCPgvh1VCCpR` z{rR%~ub56@k;Y)iLeiK(*B|`B)Lx+5t*&Je$vA-Ya;=8hWO+drV!Y&xKBQM#-T$Wl zx43*&{b4}x{^G$*P0dj2`ls3Iiw(DjV`qBXPAb+i_c14CqP*mau%qX4FMEC4{^k9% z^e@i4OuAYXCS^G<*-)f@(ceIX7cDb%?ppZy@YyZjzN~YUrMbx`zhm=H^o%Y@|EJ}N_RV!YaHxY&hnp~y4lCU zM1E}Zm#%e9Cc_JIX_5;c67@y}EQti`NreebuA8T~vozCwN%g@5XySEC8x4=pGlU$O z{yIxZ1)u1dgg4bVT)N%w9WS(aGrcX|%nr-0KbYk(Mf7Op^XX6~& zc{OJDRZu>OgLkS-|2QlXYH|InvOOP^Tk@S$Iic8ObsSCX0HeX$+2wXx^XoR2O|Gs@ivbbj# z1$+6XB*_rtYy#ZqY&m0z&U2k;7{^8R=q7r6`r$-z2Zf^M#U=P?FW~p?#NZ;x1Qx=>S?fYL2d<_hAWHsdw zg#{`27C74)bogmjJ~ehdK)e1z3R#knDN=CEYpg1_I8yd;W!>eFScJrkRAgeej^1jM z!N;q-cJ!ktK-Ot(>$o{Y@U2c4(~?Q?WF4LOzr!D!V)=MX&I;)l-^=)KKQV0iPXHQh!+B^+$Us>l`^!|v+)GHP}O z;o;l$I-)6e7h6(#PZIUdP>2W@+YL+aXS$3qJ+PWLjJ=KVO@xF6hsI^Rc1xpcgJq+t z_@x8v6>+_Z3SWCSmZ-MVQNgFXY?dtF(S{9IFBDTc6HdvmNj5J^)!BYgjr^20Yjs)Q z;^W6SG@?RO*~;!sj?wbZ!D-f6&h`HYh+v)J@Xl+zXmF@`CVLBM9hI|faUsG)qvI>m zQgq&D>SR8_p&@jbPSLci+IQPuzIw*IhsP)+16eT~$`|Qy+Llq|Q)=ODh>=Gq-06Zq zkMG{2DA@E4aneU0}aWAk_52Oh6XdLhq@zH6HoyC=chrSGMbX@B`CD$SKb ztDw?4mIs`Cd6IJhzt4ImR~2A^FA^2hF%~U1<1g$$Wp}$UNtm!Cgsi|H!m@KBSV2M! z4bj(m?R#&iHslQ*cL_T(a&MxlZe7wiH|-dykD3X(}dWo zbs**|b>0d(tk2uF-FkTO3)kO)d5wM^m%&TjIJMfCDRR$+q-Nd*iQeZdP>n20zEhs9 zT=Blw16@ctyfxLhQq^9O_Z@GCDKsFadc6W+9_LXxBcm}?QblIHFUcfdDCZQL>b)g{ z1uMH@8Ba``uGTxoBo9ARSF))pW2z20A+s30|4d&*um!!jok1xN5B~kQb(q)qMMxnN z=h=(=DDJX5!d2ojW`Rxsj^)~k^ILZ9T3&r99VG5jHV;i;v?(R|2$_-~t)Gu(RYAa4 z1)+v-so4>ZZ~fUS9izGUM7c;(10)RnwqmWRjr%JK8^q^PD;5Y5W*N)2dAz+U4g8gp zngZ0hMKc_M)8WI{jaTTFOF^`i1QB^yrg2WGjnun>vvpj23d20b7mcs(uVLJ0&#ypJ zeL``p?{WHI8t>_*N(-D>e(+9x?B;q855^%KLC%TpLh)2%W>+{zidjUxU_K3UEW@#4 ziN@O=L!Vq{6Wmj`a}=4%u|ztB8nOQ>Tf9aixBO2i$spT zv&cV&=;w)2dhvI&y-+Os;KH+_Ym^9i9gzm-51YycBOT8p+hs`?rHbk{hezkR2vbee z6$xKQsF90wJa-KlLG`rj&3~Z*xe>~lN?U-5B{fBnla(pCG=90zR*J{}cB4!4)4LlL zRQ{4k2)`Bd3Z`z-Hb@K-Vx8dZ%9K+=a=z8=En+~3HpAw!viCK$B~H|s!bThZw+AL7Q2Lf&yjzl zywaLxvN~^4-^_TB>x89o!?7kaM!r8&J78|*9tHjL>}X-^{A_|&$dfb&OX8>cQZXxs z-P1^*az~rSb&Sj~Ta=MW?s1u8BV&`IWGZ-X^XHC+YHaT;i3o0LS{*Jg%`mlzPA5CD z6>a1*>g4Y+!-G@yLumK%MoYDW!2HJowgGNHA$Wf5__yfW5)C~x>JbI(k++Fs|Hwu! zD44BAUeevS@VA((yJ_c0hURNT66@#Y(cjQOgN_KOxx@Zo1joTh4O1(o?~}}w#8dyKGbl$n z@&|%aknjW?KBIjCNlWsd-BD!vxdgnUc=5_)Elj15vzRHMnEv@Bd}|z zt{ZByK6M@y;0Hi=#zjqTSGpWiSs-GoOhaA$w^un^YDxLz@gwX-n%N|v6bW0{$@9h3 zD)T$}uvR>MU6mPGs!}3gBln%(Dx2+i-Q*re$94~q?N4S?idNnYTm|2@px1fin@_NL`sZODJ^SEN(dB-x} zQsn*l)rl8WGHlgz8Uy7OBHZ2>fupsjLq#n$lS^`JIh_x1sn$<<6)e6Pgk{`=sRj;Z zp1MJyujf(Qsv=sgqjOUS91ve8zbn}EEEcsEH!4?Td5b3q&67lQV1bCr1#G=sXAg+7 z`l#M<%t5z+>rHs0aUpiqpi{J=e$l`fB9UDr{GCEdT76E}4aY!3A_~5T(GsaQ*DU~w^5kGK({Pz2PcDHou z{5`bLL6Pzha`KF4*z34By{FycQ5oIon1JqQB-@LRQrBEsJqpgug>5jWCfFJ}?1CJ| zVZ8Kr;hBHsa_Ll4bk%kH^aF2UhQvllVX6?>^)Q8n>Uqi#iC)-s@2?VvkdG43tnx-9 z@Xx?n#t$IZu@TQ2n1dA{i1&Z6ZHQUTL7ZgVmfa+}HV^i^h#@MMgubZub|>qZ-5E!p zhL9xz8bR3R6O-^w^k%R;uGaDmEcVrJ0KP}zA!6&p?Uk$V{FK$PzHUgAaB;gc?P|8{ zreOyR{`%G*^R;;tGk~rz+IjXpfJ|~Y|Mic*0>b{6#;oOy(ksJZtY<*A?mo%NeAdsy zq_>$&=)4WGxf-KIYV1Rpr6N^{1$J zOs&nWH>ynM__>&!s|5q&kvsBYD$RdQMNrvO3H*j``_SzC8CdkPZ5`zhIG zo2Ty`|9G6=WI})LH26t?_i*mr{^i@LcsWhesEr@e3eRlCPp`@cN5n6~%TlOzK0iUY2 zPduB`9pTb_Z>7e@^_9_8=y^W({418}V^JchOX8{YHyfkxIE4`{A@wNaWn&H|w;OztDRD_m^hjtCW-yZlN5(BiNuA07+c2#8*&F*cKL{DG{{zRNqF@utxF)D2{Q2#*7Q+xVD&ewYGy>`;bSTL zWtUIzO*l7(bj@Yu)Q$9D)P%<2QY~^y5zFKVMGa-f_Ug@JgPLQ;yvyGVEdbGgDHqa7 z3Dl;2jhc^&BR%)SiQ&zlcGzVvm75DzZNy`(x6$rtCeZyw-xrl&DmAi31%j?2ui(|wK-%F`yY4&b&lAUC>trK%>7Uj;MlRcM#>(u;iGhfeH(Vc_i0RdO(LZ1G{>`U6 zdBgBR5Zxq`jKYcC8Jn;vE4;B8vu>==Pej}p1}nd{92cwzTpR2T-P%#M#hLZVwJ=mW zo&l>EcXHLg(5r||WivG}r!tw>RjSW0tI24^F$xZqr;<_+y>yD9y>U(zzxt$9wBeLG zUtTK{A!Nb)EV+Kbm=h{(6q=FiGTk3*+9JI&!ApX52wv?WFxV1HXL!?PRv_tZo-BTv zPwnJ6f_fvO--|m}IkF^Jg!dwnf?B-jA5bVKUt#}oX$tXtR7R7yAeK9#uT3+z5fjvt z3DOu_<}T6d`Q(xsELubI9zv;R{;w6A1H&~FJ>j%3Psi#Oo|@3uqOqbPW1z)}O6l-n zx?%md{ssk-5A@+Ce-s^9C+A8%%BRxYeK$`Zd-6D0z=uN#T%Io2RPnVCEG|cFY!UWO z6o+~o;JT?L=bn*C0DEe|^;g;IPxRF5X$RydVvtmVj>bno@=~oBNy&_aZofh9Pxw`g z{A+9TFj)iW4kVWFRNk(EFN#c7wdGQ0*+n#|6W9DrW`u&?MSwXYTJJ%wJ`03y{Zn*e`MJ0zRZmu9wxfO?arO}&0Sy3&&7DWrarst$ z3kzv*-wHv|QRcA9N}cn8^j{bW)~%WmD=?vyvxNxAQvD+1T8?|vRb$ZhSspIrnOROA zy;{^SJeBYo>b#4S($%ZQg6-t?OS8-EWai}03Ke(NWo?BmFB;vKS9IPTVr*DLuPwqB zUN%@uVMnuZ`I~dIx-53bZwzB$qmPewjR}?y5|hjJz*oroVj}J{|2P31>@ps2WfC|= zjcWNODyAXMObF&=oxFOF`Z`(fcaA$A4W4y#MmPGqhWuZ6akR$RYi3bFv9s*S1WlYd zdWs4Osdo)=#W1bRx+l$(A>^$`hD5c>m$jL>nWfg2<$q&&hDT#Wpat<1#(l%F*8D=1 zgD=oi>YRI~~MXVNP9Cn5}qR8ab8yj7(2&uk>ZwC|m`r)k3rBeeL$-7S{t5H2SC zMM7f&2n8RtqNPlsV3PT!KDKv|=e_OffF$zN8aK2j6O0&ny}UO8a-$Nsx(Ml6{oq#Z zs}`uGrpSXQw>BS$yGM2dE{cZSrw#8I3f2m%@r`)b$^R%)a+%4M6A?ynPNn|3IjbmW zoZzrcD=`m2V7izp{8beU-bJ*$Tqlm2%0sofqcW z%8W8&5S=$Lk@3maBp>Y+HKDv^eLR6fu=e&ffu#?%SnJDt{t(M1A=4;>c>-_e5MDZ1 z8#U**l=I3fm2tl3vmdn7Cgdki9BQq@1~>2kBm#VcB8G@C zJ06-k^Bod4XB0)kbI-`*QJP_qjARac?0sZ3py}3H#YK!YGlN?x+3sPEW?goQTZv3! zV8zNTYvdxq9;r?}E7!reA!Ira8ppELhBm+^f4nY_Cg|82c7Ixhb33<{f=^(BD*g|7 zCFyxnoa4O}MOwH>2d+2pVpStbz~xrl=8+;#&0cw!isrH(){j8@Cp8c=xeWW}jXU-7 zX}m?ky^-^`+SyCF%96p_is?&Tr$WTFu`a6S*FhgGX>Pjma9^~K=n>7FB5^kDm8s>1 z%X4r0C^jcYEs}R7z75gse}<@qV!k-n9~ID76yzAEjtMj$WVfiWqFewUh=Ps$w^19` zLg?_TT3>wh8?<>LL?%dr%^bKpTAg~IaZvKU@mpU>+nB0a1uvKH$>_)Rn0!uyZnvyg z`GD;4)idn>5BgW9hAKO$U(NVaelFdSl_3o!R-pA{3%>?x$#nOj>=$ZZaOur_#qs}sBqbm(f|9o>z_heQtZt3`S(2|@LJffK3gAK*}v%Vt)n1TVAHS3 zX4wPI-*odThlIEyc^eYoQ9cc7YgoM!nV~am2+#u3k5MoAXg==5<2XtlOvRJ=QlS~l zt7x?ct?yZL|GeK=E#lCx8Mnq}6sf-&;X5xO$?YzeAn?1AB zX$vlQ#9Ebbje;G?9v_cfrb?%OT2Sss=C#kyXP`T!UJ;PN@%~p$t>F5jZTh8av1W+& zNvQCqE#5njMk9)ozPYqkoB?eWv7;Q(o#jG;P|uH?V^xjdC|drpi`kcFA%y9CQW5<9 zlTh<|mTijRCdJc_;lu8~p>0@u2av6AzW&V+F6S4sK|fNJF8X%ZtVOE7MX#of=AI)|~uj zDM)Mg!`dvS>>>+5b6QiK5vjQ^!Zrf2@MX_z3i?UG1yy?tw4XMn^AK6tSE244f&~ zOl_EC((O?>_SV3ffNR~tlKvDnWHd0VBmxm0iOJTT-;j3^0aU%1TxD++W#?+aA_MbG zYiOA#Cu1|uXACN?#yQj&{ExHmFKZ#RCg&FO=qmaS821vz3kx}|+eUiDQL7HA%cny% zN|G93L)4c9Df`k=@7U=&QHmB&3Tq$@Gz-U!-91<5l_Vr2Y10kg@@<(^Bj!@~i};eh zF?0gHnmumQNUXK@?~_D`6j->6PtwOW8G^3N)EUQ{xgSzLxoPnR&n2usA?Tk7n`GDq zNmOF$%smUS5+^{brQNGuCQf?svGnjK*p{tpn+S7F5a`5Ay_={s+{F@ai$MEZo$6Ui z2|EAeUl_Zj!n!Jj>X(uzy0E(zrw1#qmi9Yb`hOige{7Q?Q9``X9%JVfEe)Yc$^g2s z9>A6xu{&8?QRjFn&sx2_G4kv#RScGHT4I*qT{ySU>DXiL+%z{gFWw8-_aIfp`G|Qd zQH?bGaTA=jzvERH`)t}v@AAYFKg_W8f*JEi6q>T;OU;h^!lLHQD`zD<^=Aj8nktEfe@xT$mW~$JCavAVjgukL}DW2b&JEWw1|z7 zu!1t4aWL;!ua{jqgX`yYZa>TnQ=kz;vCZ!td!KI>jf}2GCbIOhre{%ihO$qR!+0KB zAi`D`w$3S2;3rm1$b_AvNB> zIoM>A12PDTb&eY#h-ypcZ& zP?QIyHm=tePe;FsTrDE>2>?W89dxR1dG(tam+*e1y^!M$meR-GS6XgQqPV__KYWT2kDD6fyKxeMCzP@U9B^ve%VMx-sqCT z48$OH3{RYi!rfxdiU?;_+Q2*@5x#Gpbb$mf)=a=@!mG7?D@Z*&K7qZ?jZA;;7C$WX z0y`w1`=zst3Yv8IU~;i`8PTW2_$Xi^YYuhKQ88yT=gY;PH5MQg@Lv{yYUo~YR=|+Z z*2s(WOkTq0vVlR^=?)*7M@?)U)QmaF}bR2&Q=#+}qN3&k%fkmw9l zQy}o!;}}b@HdNVtqMn`9^OeK$AU$X>87)xUVjR5x5jwJoK3fN*AeQ{khNzwDh3J4r z62o5j!IC2!XtjL-uqDP2!DJmY=|2z$(zN$`O~nI6!VE=yF2m&P3%Uq}*PHZ~-#(-Y zNFS)Be=oNT-Kts^q)m4rI?x5nqinX@zrgN<6w3Mpw*wYfmEwZs8Up9|aUnK=TS&@I zOvCv{l%PCaI7mucr#j6U^UK_C2?`;qhP#Xa%dB!+Qg%^gLW#2$I;a?#{F}b+%HfrsperRtg#kIP-c`>J%tNFYQvO~vST(Ul~+lg3@5IoVqa@{rDp!!$%Rb_e6td5kvuRB0g?@@#s=o+61p!b{5@In zzSibG7cj9XUM{kaketSXhizWd@<(3t*<0TZ`Hy;kf0_4oH=4mwQQ8#8IqXRJee39R z^64SWWM9Ip{t(ODXA~Azgrf@*m1D@=iC+9zZOg`A&^9rW$|5W5VArro7KDL0O^|1c zKo~!tu59_1EV7TLS*s-(t4UJ+GBGWY8aU{fj**NBpYlN~`M@MP0b0!}#YS|LG;8!Z z7krfg-|xEJmGnGYZ6N0Lkn- z)D7oCeoQmy%Jw#tb9r?(TMaGlc%35q z<#KY5n>iRqLBiep)90zLA!isM(5@?=r4pw{3&Ay>nO065s&}K_h5+hKsGe@f1^!0s zXtfT8@mYxf^)eCggP{kYP)bUQQ7pscJmf7?13IY0aJtUp?-NkkYQ+;1s)e;jpuwoj zit}xhsRovM=hElQyKG^ve@SLGCH}H)Y{=D^Ir-}g_*_~)7<6Q*7ZQSlu6E}}p%IYu z@L+{(>+?ty>}nUUZ8ovser65t?^=(2a6o)U`(}Rtee)hUR%FXW>>`2zdcodNh-+L0 z(_#pNn=^E@Gp8cwYs^f>?nw?CXm#t<8XFfoo^L1M*&ZPc`D#7o_ajiB#R0G?G^y3L zrQ8(1+Ard1bT|H}qdh(2;WK)bUD-k5LQQXO}(|_;H9djRe8N9`@;ndtmXbuv-pU0%2?&;YrV zd4%cW<;_Bj<&T8^t>Gd*<5)Q|i7Cor8E>}Z--U8V*jnv)7JIvj;2k?=qyE5lAKc|+96UdLK8DkZ|tYO2=F<|s8#v5w%}K*#~*{$`cb%x zUbTH=n0xqegi9)^)s$NB#K~%V_+f$oLcXN?!LCdVd29zLXLwCh z!)Izeq-Sp^GUUoaN3^@&u?+i?tx;n)D}P63*p?r^sV20cbTTkF8)CthAAxhpdEeA8 zy5ix88Jc>&uLcT)%guOZ>@qDLa?hrUJBbQ8d-2KnkWAmC8 zTwH+YgM7*-y&2E;s_jpazLu6Zw4s~O_*hlD-L2=-lQfs%yg^d!T58G*h3ob+8*8oV zkifP_-;3QImJeU?+W$5hHx2n!&x}WLmV(X|=5X24l^zjnTFShD@<_-dkM-iZ{E_sT znQXW^VAZk=a^>saB=B1m8&m}g5cFaJH@32*M9}c1FU{si-e7JxCP?fq*8?=id3Z)D zg}15sInFgCzEKZLv+;WQDssB4>ge%hmvpHQ3O>EYrrsMJrW@R=8hxLjCW zYYOL2W)2^FWuAji5KhW|bSJjr733s&e(jvTW|`!Hqrao)-L(y zlBAj=8==_Kv&Ih4sV0v% z&kv#LL?$kB$vTmN?zxDgH;VY;o&pCxq6-@L>p*L$Yp)`WO;P?ttJ1;mn!nLo*mdg` zF0^Ws6xisEAuYYe&ViZnyUU(;F$2ukA_Z9ALFJFs~MS}X>hSGVp3Rd{r)=|>O8gB2(o-Y5W+5gS&S0+FV<0MX2S%%GO1_5DQa&eQ*o`>%Q;rR z!_DzU^t|cl^nD^#x>B{rVz2A3V$ulr`N`%OLp(NBM(++RjZMb)u%J8MoH5zg*G3pl zGYFrV^+$dhnc*LEn&?&iuF$8pt7%L*@&B)7m)B57NxI4EGQ(E6M#9{SMf@X`2SE^j%QXA*w*77H3~i%YKX~b#iC!hQt#U_dY0YH#7E*Z? z`m^!Z^x-4FLoUHH8ovtVzdPuzagQYuvDnHDCxu`326`w>LP^nqo!eL=zdB@nPjNGX| zRsbCwe?reZvAd4RR6r+8=mMeW=z6~dg!6jQk2u1Bt^kfc&`33SD$@(CdQo7D-Qs3I zrV8nt%(C$97H)9b6Zc|#6E9Y4@ptc&Qa_?!ng1Oauh*vRx7M_s#%f8B-^bl_A>e~! z*$`WpCLQm`A)7;@#lu!e_ko^)8odf{-%4JMn`bsxQw^+fak)C=Qty5DM`k<=y(D1r zMPjjku*6xJwU`b1uaDGuxZ5$^(Yy#qTh*5)7z`s4MM=wQoe({}g0!rd#CiZuT((Vt zJ^{X34TqM?_H7h1mv-bPE-#PqH%{u1zsWk}PG;E3@Hp??(%k-F)*P!6>z5=G2zNpD z^Y}m|6+McBb`ew;UN=8Jy*=KSu#E+?3~v1*VuSp`eEFGKJveG?|B23|G(kFlpuJ0! znt5|73+k@~0T4pJ(!+o1_1U6oAlpz?Yi?p7@e|c02wN)9I`vAt#n*z@Kq%Z0j=JGu z+v8XWH(iK5*ViNaDIZp4_2JBtH%e>;Zy@O|sZXd0qaS_pc$A)AO(dG@bwBjwb$b-X zyGknUGaQ79q{G}Z2v6!**W)LTNZnrPrcR`&c=y%~U#n+LK_}Ovtfj3ih+G~%RX@x6UBChvD^jy(^z_Fm zGOUw$JVsZeBzlx`SDA`wYbSOtX~}rj(LGtG+&a5DM8({vs7hJHf@21-hY^Yvv`I%! zvbTj{{2|M5&6$GWaK)Qni9G|=IJe{=^Y7K1lw1}&%Ae_qfZSv4O;Fv|8ZI|sW9l)E z%FuU`vhAosUEHZsX3C!YmR|z4A9lmXD>O^8mYw;jcd7qrFuF6=gVVVI|M1Ynf8cfv zR^6GgvFF5Qv;rTG>KexQ0s>=(oPUJV{VGf-hOju_Jo7%)1iFeuu6lp#=hb5-rrO-H zHa3*);BliUKR^f|-RJ2fUa9-DGlbdL1Wo^GU80Qi2LZ8_gc{}cT#40=l$p1To^xs1 z@dgtsar2*t$@7%G;dcvfd^HK0TxOVI>}Z&;&nNM-G50k|IsdVBnmmcyX_31D_m8iZ zr-8VF%OZ+&6TT+WxOf~duH`Jfr2NRNljWXGrGWQo@Gx~=bt&^yL1iSdKlA^@bLM5- z5T_?72Iji#6q(fYH;v+%ZMN4iayF;GZy!bTPC~5qEc$R2+g!rJ-f+0J7ITpveTpG{ zVdO;=VlsvqQ1k1bXf!&}$>T^1g?D*{-L}Xl>$po{W4XFCUw07f|B zKO?+FB{_BLyv)TbL+Lk=p;15DfQ`0(E!7sJP&X4$PlRZ-@S-8*i0>fz+`Wy z`c(B6VH_*0=J(4)wFPU=hWXWex`I)jPP$Cq1P4!vh_+x zv#J1&axJN~%qnfNN=Jxxw3Je!<{-%_;HXboO+r&e_|I{5WDj}fHMGcviM(h2ST|HvVZ&u3 zP5R_XB!uo(ht-4NLphuZMH961eI!dn!U72=p9XAeYst=070>N)*<>&tbkJa&Z%tug zC^Wg(4`BI;K}y4fHSI*50|+YU792 z6L9&@4{nEkih9khv&q)$Uiz@));MU-^}qFkU%hEfm@^*13mCxFrK$>&jFsy8aZ?@YtKqsv zy{O_~N9}u6slTfcQJQ)Y7D7FM$_aK1n^k8;B_Kf)FJemSRJR|RCI&zuZ|+F7Z{uMR zoNf=C$Ui1B2YG8ehdRGWePXWeCxLvFeQ?}U4;wv|(e;!8uE!fC6?F~H2Yvzy8v)GoRD#R=ka36Crrfz-HFoDGw#u7N1qWO4 ztz$)`Tw2zT&9RqV^k{;231COr5cBK$c z6(`Zug<;_uL=kPjn-(p6MfZ{rPc3p>~i1gsR$Y1gtOOp;f2_mzX*Z5^xobd%)+<0B^7$w9?lZT zTV`j>x4f^=cuk20!j#INoKFRnSbeOzZT6pyd8&x|{l^cEn{_vs?lw+s0dIrBClv(Y zSmUPnB4pb<7LKVKitGnNnAp`xl6teIoPJ0xK~8%9Bnh3g=jZqIfIXf&Ihe0A&EK$o&BeD77C$OW zYcsSw>IUAyPzVexEiW_A>@^%bv$slR8~4r;w<+8K&%MN1ls!}WX-8RD{G0-kuC zFqdS%R~(PgW~t{%M+vJ}KC@?a3A2u41p)no6me=dx=J`^@;t2_wt`PCR+fT~mU2}@ zjhX!_UfdnclkNWD)1+xEY{X@WA{8v1dz~xB+w~1H*1qGF-dluTUNAs9n9=`QA3px# z@(jfo+vMH@;vBJe+!K7<*?&}4V#Gfn zB^UU@=7tta<)R+3XYb<(?H044w)D17^LAc}BHSXU*CTHvL<&H~T;@tNLIFhLB&*5t zOVaZSy&5%3$T*14{qtrmb%tWPl5Zy~rXsHG`at}8-rLy2aAd^w%t*Nq^P~P~mM!D@ z_deVG9ArmId4`+@WZT=kwqI~j(wo z!9$kjhT9x=u%UAQqbQ-a!SapHLGCcC&l(MmL71^E7b#E)O>g>;e0vMV+baTg)_+>J z$@Z82;Qol5f0e8|{2D#7Nf}8^cQ0|_1ua^%uKqlUWkr3mv~JAWYxClE&YRcMAym+N zwg28i3l$W&n-#l(vap3|f|;aFjTZn+SGtXZwVK;&T5CnWZ39z(X4gM;DHGofz8lS| zNT8vrj+R^ea^+wLd-W0+O{utC*wv?=p`bUtMm>**vR93_m@S&9gj(X%#=7-#+~dX8 z@JaUa0GH%3XhDwzv~`LBr)8A}!yeLQ1Qhi)`QP~b13aJ0sV>fn`fk{sn>2lAk^g>k z{p|1Io1R5pf4{^b# zk3naOy(nLB({JNGxDgv&*C9^9d}Fx0TDaoiAao2ab6})Z*a&_G-k>(GC99tT?DpQL zPn9k}07QMTv2*4D#hb8xYptb@DE8c`jYsE%-U)|jzPCPNA6XghvQS6@I?)zNdQ20P z@nkk&y24t%kN&R1*)%_V<%2*4Ay$z|<>wr^#yat>qsiaDEp5m_zs&xkHW`cpwmxeF z(4Ue>cf;T-hEDl?5Uk%+ipMV`kt6)6u}BaNiWTTM&U}Se;4m~0r#Q0$&8WXGh_!C}Uy`bf0n4ad ziL}>mbatWZz^$nQCNKX@YHAn$50lB|OG^8aIlb)Vz%5ei!rC3Lt(V>*+*!wxOs=cf zpC%wSI#PLOhO8CJ4p834MpNl_+-E`;FI}YF8S$3yND=ve4Z0i)^ltvDpEWeVo!4!M zoCoq4WspY3j9`wLSN~)}vbH<~;AS+uEo=hX3{5Je}h(NMC@ z&1}@jOX0OpC5|r($C>3Il3}m!X;Y1i4Nb)LL%`1RUOsPzNhajrSjAauZm6mQnQkNB zQpIvBdRls4oFHw>gCPH9c8df;#v<+piI9YLuz|T&mpiZLr{{}(wxGe0*rWr_-r?1C z(127aKWlw_=`e+%b-6Vg^|w6fh)itZ^Di*bpB;^Pd23$p6M*q*i?6yy$_uJ5XXU#5 z%0Wi@E)leAlXweiRUck2r}u7GKN(WGV0dGV{cno4zzy~l9~R2@ACZ0sXb($l!wQm~ zEarVzWb0S?A2wg@>pvylp2kxYj{8np7r9|2YTe(|lqm{yT~>uFAcQ6A`NKqovVF`G zAS*i*6}Eed7r{~wMdZCzefmO;@IqGM=O^)#5{6^TJQkD0KXZwCPQknXlN(47QBtT# zBevG`;sII+kG$YY& zi#S|P7&TV_kHq4jh&k$|xJ}PZo|#H3H)5@B1J`ejolNgg z1w)uH{NvIKMMW1J9lSYjeESo%r;P1y05KML{nLY?u(dyMEu1TkDlNb5YB@t7(LZ z&g%KWYJK3qE9~u)nRV4LI$0HDM`)PD8DPMZ0^`!nj{q`>#)?>S#(z-V#Dm%&a!#zy zZU*MWWwYGGTm$c8I@;!9(hOb$uE0}U;L8XuRfB1X%bLLp!5>UM(IkaFFLdi!iXFSm z!L%nWChSEiFyXa~Zw}*n17?5PpFxb&VOw(P$&eZzo}UjN`?~9JWFjE?ks%D-Pu;e< zlQ-*58u;{_OVlD7Kjx}_NKH{y&`aSSIxGeNCn2USoEnR`%p#k5N&*2HGrv=?DVy+T zEDil=6c1Hl=j2UpioFzS*^P+>1k(OEwAzObpMx4*l@uuDC5>^1hv+B<8m_09SFWJIO_ zp8b}6dS)m+W+<)DAQ!4y7|M#S67zSPDIIAXIhA!ttw&F?6scxamWfj~0`^hU&92iX z+~Lfp`i12^^4P|oaCzUkBeylN7k2Rb^Qo$sb#?>Z42%Fs^Hshly+M^ z$M5yGW2fGiYBxQ9Dqt@?VE7=}zyH1CEKn+$zWjUPNcWjE#NW^xL7Az5wv#!lD&EbV zuP~r7xA+ow_Zd*z2b@WeJJK5YNq|BZd;V#v#AMynAePU+1i;ou-q@7{Y-zRl&nUQz zV)%J;dQW_j`!_y{kk1D%o;r<`Cm0xJ_!k_zuv&F!ZXvIQK{i1<-C;#+gC5t-WDhOJP`>oIGmX-wo#7kjS>CV=ofTZI^gNHjlk2gjSUCZ%c zC`gsH6HJHiZ`Niza-U82jgMN=m@@b2;(G*8)j;Y|ef{=h`v+s*LU!B!pLGn~V(W1? zR?iCYTzdlWXZgo2MF0Hq`iGjF)xYgRn)bObhyUdr@(7;E-2(Q|u(sEDV=*Z)>Uyk) z6~;b&I=`7b(~m9!{Q00I`3T!Ktfj_$VJ zrVBV&e|WGfWsqh0<+?~DB5$PXq*!g^yx{|IdZHz+lzd+N%A~{;u)k=o9(Tv+n%qXL zeZ~Ey>X|ORc!d@Z11ehNOh_3yUaCsAAwBNN87Ygdx-XcqsR=u79dJ$Nvm>$ zN`5+{-E?;o09>36)I1tqxFi8FcSzVud}-UYbvQX=pE*FQ7XW)&ZrFM=GPQ@1I*iH=8W=EejN_T8-^%;9=^?>I^(I72M z&|5U;j7;-yRKrAJUeRXLZ;j*fZ%PIZHW43d1<+;oRR=(k*~ERQ7=NJP|Lv4+U0~zB zX9_8*HXlc{4jt#27mXRvZ@4et?|JZCGpgaL^T9|7(m-3pN|4jI4)cqeq2aT$K-`AC zI>()I!WWa-*)~`gJg6u5&f1k=Scp6Ij+y8+P-zo+Cb}5__D7$JUpu^r59!Z9RPlA2 zRuNviF84Y@51lyBt0O~!-fui*J)yd96M(11cBb$u0$;I6(VLBuPM+tnR|RGG_XuO=S@MOBt~eJojtIHn^Sr`e6A&({`?pj{8ZZl^^@y8KS4Hrp2g)A52P zTv7FBDC3R?C(6s;2K8*6j2&+q40%dXnbH~x?K~pNZ-fmA z@opQsQXo}yV&nFCJ&mk+J};M<$I&JB-d+dti?o|Q-E~k141Ew3RHJ?jVB>_44nLo1+g?Dzn=63SrGF?eV&FG-nMSC;ZXtr_MGWL~X z8KdiLk}kuz!nrQs-gyNKYtidPngH*aRo;%Q_`%rvVDx{goZD!s5#m#gvHUf^_$f~s z2J#h;)wGx20%=&)PoZ8k{Uz|O?ApJia5;nL`|v4ZKPC5n7AV03GH2d@k}o?NNq<~- z)jXaXbp5+8$Hj-H7kPI?wZ^{}vNVboapiA$jJ#JfMyQ7m;W%3?mRIG>QBjIj;AG7A z%~V*yI9pH0*-Z505HT)#nRm>br=i2u^xA%hO&&Kq$SjVctV!jckiJ z-m#&}aOie;WvM~bnCIuXblzf(Wm`;>fr8rX=hVQIC%R&{u&*-q{+7G=RWw@meEfQ| zK@BostPaM{j{&j)sDO&uHJ6AMo=kQwhY^QM8ecDXgKh}jM&L#83=Rvko(!K3QjPSqaP2 zW1wy#Nu$62iU{5A#*Oe#zCg1IfXg|$KgTxV;GRKS4^;Tfi%X>+bWN79+r+oazIayc zSv3fALKG4)U9jneT2;a{>=js{dzZ=8?v(XQnoHYjaO70WGpN+o%(sL~%Pu9)*k2`L z2w47}v|(a)pf4L4%3;y$O6yp^pn1ba7L-7j-ZGQ9qVhC6l-WOA(c4U8N9c#i1k%8_ zEF!0j(2&l#=4Aw&y|TEjq_@fR<}oK9XoZ_dXzCGYUzk*ib6_VnkbrjS#cI##t`$&m zb71G@$Tb&q+gd%xk1+&!&FCOzYVD5g=JklJ!Aho!F!+$&SxDow3ny&0okX@%Tdto| zRG=jeV^>FK!Z0qidkmwaZcFI3a?L`pETrk~OvDOQNiSsMpTR`DqnW_>k_cDQS%B{J z4w&TP=EEcpb5u}d7-Da&ld>ILN6%hnwON)Lk0oY_a^^66xs86z4#v`)K6-JW^Pf>V zrJ`ti^HlIq(3S=K*e@OIf{~r~!y+`>f)Mn(qvk?f@w|&<83BIX|1^^h%4SS}Y|8Ks@hFt_PZP zyl+{BoJBE+}_Qd+8C?H69Lw54p=Y!Fz1e0?D5&+1}Q z#ma&B;|2=yWu;Ead{S#*uVND3##cN#tRJO@cGj1*l-`yO^#)BC-@m+iowMK=C;p+x zNn6C1F2O3^^ir}vu+7f}qC8{uk{eS&My^iTqKY|4D0^Ko?@=87v{s}{Fr1n5K`1p< zd6(j#>C}+n%&PrH$2g;=#$RmOk*A#SL^}udz*dGor^pL`SI2n-2b-+Pe39%R- zKVessHInkiK}`)*jUy0h1U(_-Q(5Y2LR$~tKWxoS1tA<8XYiGwGyE2b9_E-P_A@Hn z86{zh(!RVzT_?tYN)v$QayKzokMV7d+s26U4mUY zmbT-`GfhIO2L&unNMOpH9HDEOj{|3^bS9upJqu-3Lzoe7Nb)&n>r-WK)6u&F38gdM z4+CIo$|18IQ?n#u4r+JY84y8rk}2Q&!I#ky4s~rE+cu&ceIVDj(+bZf;N*wzDPQJK z9ZV-|BI=b~QkGYU?SI$zmhF|T({&(WYq>yw$GomjJHnXV@sVc)F#5svR6CQ@gFz!6 z^w9~l-uPC79#Okp+=t=Ep1s6S@hNEZ;2^3JP|oO^*Yg!96Wsrl0H>Ee$xZ>@UD0=O zD;i}neqOR5Vzz!KufT&BiILT{YB`{&tP%a=EwAd#+$`#36p7Ki&nEi}%@Rm~%%7Kx z`V`AX2P49P*lFnXSra~5n4OLameA$w)d4soNbzbzaGxwh8x1TL?Y{x(zBc4Z66k90 zqRzz8L+_TF0Fr*5aIO!P^HL{vH7O0{Hj+0)`R+3o*6(Gu`@wH+7w>_1=kiB=cLy;$ zGBWhCiacRl^2z2CL~kaFMh>fvq2*Vr@x9c~IvmxPB}VnhsB%^8;^P86fzFv>rn%v5 ztm1GNZ5C-#8AL2$^u~>}?j{A-yeQM}80*j#U5p)a}5ochSGBvf1-eA0@a{(Z@PWWmA-t zoG~<)y3l#SuE*2U94%T=wUIg5N##zyrWs}vx)H~X(Pdo(H>^C%X>cG(T~}B zvT2bXDj-k}KHR^zu?55YlrPsKoogDN?G zou{l1Mrqpf<(5w#x@y)dF|3~kE$W0|wL-k6)2{Bw`bRZ}1-^nSMEJ~n>w(;9qQW~j zZ2Jyyu4@IIgR1JYld3NYH78ozlyhw3IG$4&b2k!-qM#HaN30oAE#s2_WLr5>5)0vtbcdDKOLd>NLujjA+^aorhvaPCi z@e7s(1RUz}%hf3DKJF zuIxx+9RnPZ&ee}=EKh$3#18`dapk~Ibh5Td!9LA<*;?qXb3b{U)fnbSyBu{DdfB?4 ztniiN1Pq{<3rK&E$>WjKp!a)HBKm7)7wongxs?%OEka6iIz;ndsr{f73DjDkjk=mm zflrPOJq5ENe6m{tr5vJiGWT@-&M-l;7Zr~8=*i}M6x#!Txl-BZy|m=ZUtW^SgXq3! z=tbfj6+v#}pOex)>MErGQ=Y3L(V|)Mhl)OB7h$JFbUqDh(1dPXunW6+sBeeh4)YjjeSuY@juMo%x8ah5V8&OQ&kZn`* z_PiWXT9^Cs&Rxn;OwBHTfON1-UQjAaY0alnGO%L8XReh@tE>V~ZG6$4VP4%0< z&XSB@t=DFPy%J}I)qg!B;{A3~c;7N2IQuiwBbliqe!KF#l-&^!yz}!{S29)9{IIZg z=8c;{rq>*6_FcECPJ&v-iWP}KjJ!JYGtvGimi_y?o6!<>VajjnesIu>O2`jHV|b8! zq#t~W!TcL!>7=r|<5%pM@q`!H680F_*sNs_xA8ZwB%rZ=)1{)Ii=@-dK_izC^SKez zn$wo_P0v>ZtiIdoy#19_Q}XuLqw&=`n`{W{ zoS;c-EBHKWZ7CVJP7k%yA>5D=@58?NHo&KI7IUtYMX^@(pM&*R>1TmB@K8Dc8rChElKfNWC3uEhOY5&vis3SyXD%8yA^7lPe|t z!;shk?`JNu=JLpOwP~>P=_@%mlQ4MXZ)0Xzn?(q}ol!khE zu9^W3%>f!$iy#n1kbDD6C?Kx+38~joM3RL%ms55m7wsfY_L?28w+vZV|k+V{^|W+_j#B=2Z3BoqXg}Cl@H!L~qL- z$N0MJ+?sP=sd!*%pBPYjoicJ=yF2HF=S8$N)*wvUMY}d|xUym)@p@v6cZZD%-$c_C zv`9{ktjkj&Se#F{6uOY7dAu)gum<<%#KYCP@afW{qP~l$6?h%9(89&)f`l~X!+Nun^2Ru5VVHN8rY$>YneZX{Gp}X{WJ~k+{03HiOfu{P^d2$C za+9L5#on=HQepbhBy-`b_gFv$Na#U#j_gHZf4uV4-o8`Rl%V6UeYZO;Ud8`^X#bC+ z6sW}*`Pl0Q@}2RC!h30a2Svo5O4UX2N*@5dkcQrU70Gk5huZZB@7|NvB|YwWc9m%{ zD47GmPZ?UPiz^Qu1(R>{0KS?4K{bZqZhx*7kbgW6R>@=Odtineehz=X{R!VCT}jdt zp=1N0W*QgB^?O@F>m=eiOzy;&aR9F9k_VuxywS|@X}(3|Q8BRdArNqn$OKs#;KN_H zzDQzgX=8&(UkrT}h_&w(>fRU%)W`5jzw!VeN#C#HKdC8!>7qCI>5rKEpidJ|0JJ3^ zYd`zv&xsNIdTc zp8pXE+fy)PVvMpqd_jZ_c3~o4+m*}HII;BcIE~SAXt%b?3@Hvlc z0yOXIhIWft6szSIYFH69YVYt3gEGO-LqF5=`{8!$Yu+RfDN-pv4U^Y&Of*CZqA1#> z6WU=ribdK%yt#w!8<+Z#Y}f7=E!rmw&x{Nz7er8-uONNY^%}fg=o{FSYQ&1%_R7(9 z&3l4*&M%cF#Rn7l~tdp4BBVJ7k0k`1r~ogBV;*c%hF76nkPt8y=N=R zn@ebV!=aD^ledaPy4#@-HK`GS>V0$FoQ`dd8ESOq$lkEA=5p+GGjH0w^yk=P z(AQ(5U+y}VlmH$==m4GcEz2S3#T89QuUBjT9kdog;`Y-nVO5m(ZgCo~+jBBAzCptY zS&2nSP7TGhJik_5H#hDqa4s?csnRia<*%gg)R}Tln6^)bxd2jdH%sCtOd?>)(I~?rB zzgEbnR)l*(?1ud`ry+S?@!fu>s*g58v0R<5ONbv7Dj5%p}w_V_dvLn zX2qO3{yc5$IsE}R*Y2G(x(@Zc;^h!o$e=@S;AC0ZA(JQ zBD%Yi!{)+m{e0| zD&a+Yw<2gj+#iv)J#EvkpGzc2PC|rcE%Ibb-d-)Z0DvNu=$^?@yZmKsP<%)BEgTs& zs`?ikSqnM6@H2du!(i5F5pCW=7qULmJ+H|p!_RUldxn+&!I)9Ib+*P)R6xb6*;{G; z0c_6ccZJcbkgF2P;b)l_ojSBodV_@r6r5Izt(&Hsjk6(rG51tI`aCXewclrW%YO2D z z#3P}o<}-4vP5Fp0i2=CtL!&zDlS@28ZOl@n3`(uSZDA5$;}(DMS@{ZhNx2k366z>$FyISXH|3oiDKUTX1O zzQW)z8~f;!V!k3K?fsVx-x+~+Aq?c&!u1S0x$%wHFJd*KAr8;r_`*3+~V{M-_jsxGy zIBK__0vI1cZ=GO@ z$O~yD=Srxi)=Vgc^ta=O5s@hr1X%vMl4MNZDcKs(3h65e13nlxN?Wp6$Vg?}jWCvS zIz!D>HOYa?b1>~nV8@UnTM%EDBiGR@6YX9$MDgPW!-2SbnyhrpxIzRR6q?GW1@1!2 zTg4QC9kU9moy6p-^&2tO**SL>HYwlX-4AY(11sxMzbziE+SBZ7r}WOWP>vFBe!sE0 za1hx2e)xV2+j*=iq=<`=)AN3;dWmuLD~1Of;8StC@gaQ`LCNddj)>Q@0FuO>a7^== zr$6F(@MfA-v*c^LS;uMaI4pJyM?+`g*G1R4oLnOohDeO6TF`7khb<~u)?XWS^4IME zMe5j4l|0ZE4EbXa77W*#Fa9e$g^)9~mV2-`RcqOO|7)QAj z56Pj#CPt`U?r{}PE>0e^a=ocFJ_y*Vyiu>om~YTW&2$IX(j&xkpC1{aygy6QU!$i& z<{WIxvGckoA5oDvuck(x2An$xdIlNG#-&4lQl-=BF&Zwmf$($lo?SVcUsyYSWuF_cG=1N?^wai`+JJA0G)v1yK-m(_~fWFUns^SJhC9OFYET z$CPS1l&frJRYx<2MD)vz-66f!rnQi2OjaA6h%EQ)q9oTBg?#>Agw*0nz&u+b%^Gq&IuAl#3gf-O=*TyL0K}UVnN6=%5 zjT|qEvjQ-cmY$wJkP2V4`eS0F?Z3h!k9iLieKj?J+cCD6@i)!r*Kq^F(mgZRhVX*e zSxBMy)7(e*6ufMX64GLlUBASOCa&DCOk!KSp^z`zWvzBfe z2X79^r>DH`#&DR#gw6kZFk!*vOTCOYN%2*8GDSroLj2gs8ni`VbdcmG8Wk9Ahk&7e z2&lDDYezFh)hl>Ud}G6vV0`{1Y1z+h$-K5=p{~LSV zi&+5!s?*1lDs%6Ih*|b}v*<~=ivWxx9t@0$ZWVTte1-MyM0ORe_~eD_^!?TO6-emF z=7BBU(usqAS(nCSDt+Pryalz~QUZm2D=4xVJaJsJBh;FRD1kPU8`C^{&b$e=j)t};U z?r>W_kK4W?y_YE!oi}uG6_lr*`*C!Mo`+*n=J>}+M&0sj# z(+or2QaI~X%&EZ9NXxqe>lZkBD$|27f0(#XV01s9qD^A}s=7Rj_FKiYO|zU1%m&s@ zk^1emMYL{pL}2t@5pea`Myfh#&9{K{bRgo3qsCJrt@tgav*ziP4V~52>3bG{uze7^ zv~=<5@Y=A~+)Z?Sx|pR!`n?V66=-q!vZJ6USGB(q!x>ePnx5|X3Yn@ppXzuSdgYwJ zIm*8KBX_1%PcWfDQ^d=Kql#L*BIj*^maelta>%hFsItHI?X_faeP;9iM*N{EHaI<} z8%F!;=}q6;g#{@hP`bY%n1|UX!=!iM$>y#mT+nygcL7k$qsE`yNDveG+ZQK?;WA8nQokyd{YG#T$_qg1PG8f ztDlk9^%qv%e0CO1uZlqBwrJ_?`EvV! z2-BBCz#Vbg2u!|o&wbh3tVGB?_NcAZv}bM}@4&wY!eBV|Q#^<2JSf1H`&7eFc+$aQ zNw3&7T;%Ul?)%sAe0KKSg=N4)jrNeap${_h>d(e_!O9FMmGrkwwsHFwr)Yfi4WE9G zlRx`FC$y;*x`qY)Ne$&LSo&-HIIMc9yxf_N;VjVLj&0;c092JtUTujen~4&jtZ)A= zb>g)~&rn@~277Lg^#Vcs-S`lQoJl!#Qy+CAS0`t%5XJ4_Wh-rm*y=8>^3m zQ3&Fxj6!~eYElScE7sI{HSL0qda>z*%0u2r?8U63x(KMr0Yf3=xdozW((sy3{>`?= z1vvjSmSU7R$*!&og*=kRjhixwp^jRg>DqhH^-?}VIg!V|Bb{SHH#Lu$k|DlhGn@&R zPwXE!8;J8?a3RH%eM9(f4fy{snx&fMqmrh+wh}2xlXqhSng|X$c{}pRRh5T6z|j%A z!o?oV4j{B+kc0g@dT%VmaVi0FYj00LfrsjK6&qX_hB9O#&sB*J0!lrh?IF@Wx3P@e zaRgDWfB2}*M$moxGY8yW5S$nLfP;>@(fP4@9|ml5?_=~$SG@>HF0Y`Px>U0qOB8GD z_`4?9=5$^zq41*|S5Bu9;wbg>9yA7 z=lJewW3_k-oC5~(4IS`T;~8ePf_f}6m<>m%?bv22OE#oO0ImiXE>r-Ub5mpJK2ha} zo@Z^~ydgHiG)Bi2F#lXtPK~G0s)*WYPz(?nE@6|4#ga`e7J4pHPx}plj(?M1_h~_n zvRNCZT}c0RXD1O?jk{MX{lRve);tnR?+!v%@aXc{-z2DqNy>9jV(Oh}TR*6DAIngq zu`eJQs-y1v!d%pyToqvk*4s|3XfCNw{F25R{QMbN^Y^Hw>w;k_)E-^Ec>^zSr5w zw8EaV&*~pC)>hukK7caZv&53e{-=^%GaO&lMbM`n=)8g3WqIBgjz&KJ6WhSk#Wx0f z`(~i8BpTIM@{vQlGPrsVhy0Gi6r}cfd-DWppQfpsi(%^t@VW`mDr5ebN-js#WzDs{ zY#FT@g?q+--OXVxzOVI0n69y3*z5cB_qQPa%3z-K_PQnBC)d4*g@o3}YxfAbtae<^6tq=c*Q+xxH1=Y@9*@WyD1NiPT9 zdg%a4MiQLwc?#&e>x0K^?X3a-1E%Or%THPF>S_w$!9rH_^9t_P2hc{ZA<-zuuF1yt z_=d)ooh5B%&_Ych$|@RB;M0dnv}y}jnb3OZaY$;>cM29PU;39V@M`CyEl6+e z=$2QeKLwCbxnpc9X8szb8{Am6kuG{X)Pw`Du1nUKvycKOMaMEKfEg4Jl~O)6Y7a>f zrN{Ti&=v80#$`_Qb`jpfq*Vt07m19SEvX3##~M6Zx1cy5jkL3&2P95pnDvaiMPnS> zGd4-&F#_q@FR3`R@3!!3u;&W!+8PYUQ+0k$;+YAi_Tu+xRQxucrSOBnG>SF~7Nz>W zeZr~z$ko^OASJdnJu}}P4?w&hp z$W+{RxuLqu>a-$y*&P95{fvjdp2`sDj*I{i`0YnA z{9?t_{yzmO4ni6CXc;Q5DJRpfFMgGXE=C+?_+H`UeqcT{} zOTl-0p^K#a#lVBd1g&@z+#v^3lM{UEJf7QjHwl#;kr|Th{;+TKtv5j)!esk`@^FS` z2P}d!PTnv}zG2_5y~#gD2nobLe@2aqRS!{YDfPjZsE40D{4Hs_O={8efx%Llmzlqo zZTcc2x+$qx9NOfVzWV#VbhL?Gk-rEdF=uPe5|$aV=sR%|gA2i$E+Qgvfh(J=WxUoX zP6L_wrzxofs+im0-Z~jhEpcJ~pg!@YKDQ>f+{v!-ZzBK5SiG=6u?~&J)kq?zu=x0p zgq1;qUci8sV##a0tZJbNg!uz4i~-spEoJ&MN4c1cKVE!Nr15M727)s7wo468KCgsJ zX~spIR9$q~mUGdOH0YKqPBUhQqrmbtS5xsNTKq@}=IJ#54$(wxH&LgPbAX#`mQ8e) z_f(gmXy{mvt9ruE-Uj3Lp-4m2N>0tix+QfUj|LU0%$90+%@rt1Nzr>}I)7SSP@^tD zde}aaTP#>4Jw;LfDdESTGja1q@q&w>{&HtLhtjh%KgdE30(_?=C4IYL336vxs~27g zI1Qb&SPA(wtmsPxJVS{mX=UT;5u*<+<8~F_U3!DIcM#Lmg=iFF$`>AeF-HTGLyP-% zd8j7W@i>pP!uD{J?0~%JG%9>Os2`gzAjAEmpR!szC+0kH(0c&s^IqGsU)(}xx-&6J zZohxji(xeudcqertpb1AS^SJJ}sS+|46IV-otiE zug0rj<4?EJ`McS@Y+m*x7XbUwzTj8t1IL<-!Bqf5DF3MnKr;;##dh6)qH3#q?H)RQ- zGjz#^L@$A|^QFWJXatQTugz`^ih=fL{DADtMyk%d5Mq(^^E1+SVH*TO>rqSu@g1V^ zpb1A&BPuZ(B36u(vD)phs)I6{=C$XwooYHMwGiTv6V-(Mgd#OquVLQPqW@pfV;j+a zY7U_PUtM9tFH5%<4SVLMD*`iA#@`30sIvA2R4#M|QmL|vMx74j!IDBxC12gW3;5UW z`|J)Z*%WbAY+90EjeX&}SL_XP~BtWne(b{uw^+K3bf9rfBTs zB#>D7Ay&U~ixRG(%&=d%l2#}n;FKz6Zfsn=!Eh#W0k?be4_d=`%qC^I@289SwDH6+ zy3UEIwSVM!{ONtL#(U(a5>nYWxHB^$!y_f=Pt-FZ%(&r4kyc1?IyGVi%&e}R_N+fXsmy9@Im#Q$bo{VzPs zFc;=V)lozr{`^e>3^l<&f}HITu%X%;AK4VTD}g_P8$Y_k&4G=MiJDBHV{ARH9^1E% z%~gzIhoVuu?Vk1A$)sB*2kd!3P@Cl#huUe;Tb_Jp(P6D%0v{;eqce>_t6YY;UurYVbJ@zt3)D%wt;#{GwSzXdZ z+Be0TaTZZL>KaE1bKLsEgZBWF0R9ZrEU1$2Y*Gz9xT-fY^@~s{FF!z+ki3t=YE@|C zdPMq>>!;H3_H{z06H!)KyR-i7B5+t7gEF~+P^W$g#;%`(h6)VE4?GNdwRhh~d*qI^sx7P^hoQhx-?Ao~?A?1}BAos($Rle$iO8>sEqkPlEzupkMmA*>eAz$3rctFsSyYm^au1Ey|J%{v))%PV|Q)M+)MW9Ga^C&f2I{L*o!drRk z$c&{hhIq7#^^YbD53ZLD!CjW+vt8v*&A?3{v zQD6(rAAuLomghiFw^Qb@hz1aE1Rh?BWS51DNRXBjuR=R}75WjfDYnl)Xv?r~LL~ zU}@M6oP{*&98QhKnf?=nR2BvSYY^B{YRlKk(Bl2n_R_pYn^@6f#3_<>I!khA^Pfa- z(_!ED)w_lST*ROCQLjL)-?v5X&=A;H)C0vxr_RfzjaZ7R?L;*i{)td&72vo9Kkji0c7I99D4%=uYJp8!_7RR>x9+?a^cmtM% zB>`OuKJHfQJK@^kH?=Yk-&(~3_(z|A5cGfeCa&CB>iF_R$KCSdABMACH=cG=L-M)c z!U2u;PWm$j+Kw(?W@dMDX#^GK^~K?FSJQ}53CUzlU5kmV0E$AH)Cz}(9b;*qAhLV^ z2c|L0HM_Q0aXzBIv+1sryye?`6_cgez>y7++-R-Xd_6-AJ^VPcD=Udl+nNs69Dp_B zPd%U75=g0j@55>T3St?JTyIzqTGAJ%f~m8Pg9ETi)r3)mhSo@Xij;Bc=XijaND;@| z*dU)y(J+?JdTohJvb+v$sF@TjBH_(#Fq@#=B7RT|&9ETU%va~9Y-sBe$`eLY>C+6n z95QJXx6II3P}_iflAc)_n5Df``cJOj_hSexIMa9LCcv&y>f3Yd_7 zfzvOnwl;@d^Y^T329(=~H3q~9O+eDy#$sza#%tS5Q&KAPsCRHdpyKnU4{nYc(b|=a zKzX^z9>b%_L(-*X(;Wu8Cf*(_skL8#*=6-%jv^pSY11ltajVG~BTMXZ4_ZoHE6EbD zalHtMyg{HR9XmTLybBede^TDlrvBBi{G)JeJIt-D*_S)-lM^uweGRn^A{q~{Bu9do zw3Q*I9qpJBQBiBhmm{muzO=Dzs*;-NXooi!sDzZv|G-v>ZMwsSan2Pv8g%N32YeC{ z*>AW!_KtvwS`ynIJpG2*!TE@(8Qt>AqB{QfFyQT^gs z-HnVNo?hjYq4QZ1PfF*GT8nxZ)H@T7%*f`DEas|5GY?S|xh8yZ4(H7&8&$1qC|(|$ z&SS8onQtGvJ6_S@-FxpQw z0EUZ=P%G|(Nfz?Qqw#9$a#bw>k}=P!zKI3^-@5m{{P;m03Va~Z>y%EPXLEc3BSrP9 zDXg#QOapOz=N1xB^Uw4fLr!F*2wmK&6P81m2+_&FaD`{T@N*F)OAd|x@{~%{4CL18 z?R*c6OgkKt_Ru`}?|l%xWU|sz{e6#iQ(!RsC+Fbu7iSHhhQ1bo2IHR?&(Gcm&DPq# zpclt3y18C0KMw#IlCMl#F6$P%ePu2kwAhex0UhCl2db?S}ikdk(Etp>j^d5-jV zVBqHCOOI!{88WjoF!|{lgQAFpOv_niapnwn^wB8HIZBWN^`Cz`qr7#iGuy61?%y>c z*_^0bD!Ynwtbmrm#oa!Y28%U^{gHo$*!O)U?^A=1q4qNFoFsN%=Z=M&E6FF?UBZbT zXW3J~O4J$B@wEzbIK?BCUYzqs?s}b!_e0ZeQ_etjO>FfY9hC9=Y`QDAy3Wm`tn);C z_d-oyFp<{%B%Im-HyY;o^3VYaJxxablB{6%eSjW$es6^Wz{=qlaQ3W^JTv{@ zi9TYnDp+k%4-9Z>gIXFpbs``|)o9cYZQgJeqI3(Ze(BT~UHCqBdy8?m_4e$mRao`l z#1%TvMso!D^}2>J87H{Tl?4&15`QP$k{Y{ zBEFmUIFN~wb~jN>81{{GbPpMM`k6aOZKOqr!FW?(J=~Dm#4D4do-LQKOR0%weHzj0 zYxSHy9X^80`s~AsQ~0g3e(%IkqMw09@_E&$s(D$Odel1O8?ajyR~*qdG)JpFf@R7f zlBM5czcx$ZHd9Ra{D%yx+kr@n{YoLZgRKP8At7%F7KYgfgJ*`_(nHqrP*HXTJT$g* zPcl*?m9w|!@?*`J>L;Cg-~nNr%~?@AKUSCv6#&3${IuNHEK`*^E!XVO(r;Gf9@rjO z{_;gr1gzX>kr{wayUwxpz1HNM7Kbv{nDNL|$7JW=J4GgK89ICPC^bz}0T0ZYF&AZX zs$C*Hm!k}$r1fbP9uwI-_weQ1UT0^T@+x&zu?LE_D9N$MQx6dzFhFy35{H4#ziO;% zc8Dg+p|+PUS3|q$*e}`hx--&S=`AQz*?)c;QaYzHfydS#Dhw; z`X!oW|H%?9?=QP5AK|_HvAPW}I{Bi;aB#)QFjlhtVUZfS~yXTIMY`JS{9 zs>+!Xl|0krZ?m|S28CY3C>C{84qPMt1Z6(7^D++^r>r0UzLJ!dmelnc&FkqnA8goV zQ5U@M1Ni!r_DN@ujh&t237(Koj1HdjJw=xhDDEd?27NYJT3%>X?&52G z!Dl~=`~dkj{ocq!>GHkGZlOyS&{%l0&Uk0u=c~pAx`s^6FK0p^Q@*P*k%)R?h9su6 ztf?>OdGdv}ebd9jfsWr?p4E9pJ2|!n8O?JrEPY?>wDDL-%UTR7T0J7CUS2S7+>|Lz zNkT1mTHMCVmeo$c*K5o%G(S2{hqAmyg&Te>g*RvHk%3;ixoV7N!FHO*)0kP)NC`Oj z@U*WJY~;&gp>)6Eh4JHzhPGUO2`d-~mqg0v>-raHm=X;1If$)A!z|tgQVTg+7Z6W0H$A==uSx$~$-KAg3nj57sX(Kcx}zrW1s& zArm5^MxVB|;xu63@h?`u*0Pp9 z*oqlkKlc0Zo?~)WkRoAFO-cqcl*>fojDd+0mpf&1p~ie>y3v<^LaN4GXXc};7cf8K zwmRYucBzC-D;}vlY53E$#j2ss=(uNI>p!F@?-ifTWL%O^8_Afh^riP%R7r}%B8N(^$vYoOp@Y1Wz}o!nY6cQ^veNqer=@CYy7IjGbTUo~X33as&uJN;_ZY2lY1CI^f1 zj%ABqe%;Z4>}N5RsHT_HRFma7!NBf*N|}anEl0h?f6LQYKRDV@RvE`^fGJx_qK;&K zjzc3$sHx4{s0JLP+kW-an<~q23#>7y^}iBo^hlxvr{dFafIay|Z%ty6$&wCg_Y$z| z2t39OxBhhZ`;+^3wvF5ODD&HM;E~pGZC9nA{8I14-VJJok=a9(c`5*Kh-nSlSR2Ea zu>y1tEax5I-+3BWrD+OO7zB-+vspzaNLD_>vh>3j&>WR{0Is zzKSkOW1(ZhyI1g)?tooS<$gix8BODv4IYJ_rP{IJ`h1e3aC0dCCtXI?(}Xy@YFT{GKVz_sO$p z?FSHom=cfQHzs5$S~c%lU1WDXq4mv-!-IS+8LjF^eVC-YERKrS3=}zd^=wh=&pP_N zTG?GonSi6x-JrT5>{?j(R?^4%77kf;jJg%Ry-@((h7D7_>1A7nvuX@YPT{B%(b8Q0 zHcknB_@)-=1@ZP*@<1KnR7X|uN?HU;A|2-J6gW_=W8`F02&C@q?8iOWzA9+_%{ z6?GqU#XnSAd|l-NrIO%JNXwC1#^Y*SKAf|D^oxH2UOKD4`8AC?K>w|xr&s1&)H{+Y z9n_>SMUE^4zLAzv$}|3o`PnPxWe8bO7d_0T*?rrQ_CV9f&In0kUXjSp*3az<$2x_; z33MyDmEb7Y%d8TprSQQwPOWsbR9IyZRmu}BV4!iZ-u&UZ2B`Lf1G|->c^hb66}#f zyVk&{U|$CdOVlnt%U!l6|MFTB=U=zypNvg>Xah?!&R=k23uz-v%O)$MAxp5*X;kWm z@?N3`MZiT$;LF?wU|5cYzx(ZKq5OiXIex3nH+hG+#pi^40m0zA&3c08+QJ}Q=kV8A zpS|@Aq#G*4)?YQRT>Q?wZ>s!B+N2V*=f@7!Jxynd>0Ik|^)7`1w0ScOMY?2WX;Q*h zrPG~R!T&ARW>0fX;9JL*&LKwVPkpgt-q145zva(cldN>v4VbT5t7m3b@QK9evENf^A$c<* zG8z;>+aoe!tMPXPDvh7rBC$(XkynKuR-M*MT)GW?Ae|BQOe0liO-B8O$adV7rmV>V zVGBn^@u-|QEzC$^Ml8wS+5d#?)R39y=I!5q8?*vBuSmrnTy|Bz(JG$)sw)iRQOyp~ z>auNwstI8KL^E5p{3U=(8biD8rmiy+nPpn)&pZ+KufKzzS3b24TL`DgpZ%?9Lm5fHbj{9KJOJ`I zCe^0rK%$3t8ix}%9Ql;nH4~0M`UZivX5QRWH23s9SOU8_E`t^Y zIY-#WcXXHvpoZ{y9%OP_Tz-RMDmfIlGMlLgpBT!VXC9)cy2#UxlM!NX25yh1Fg?2TvKOu` zK*X@X!&JE-AWiEo!I=hbL|E&-3P|32WQJZ|-m47z{9EZn&b9mr#CEokT|`te`mwLO zbn{bMp~fQ7sFQ4Fhv~1Zx~}Lp^?3?zZSF+Uq0!CntDMX)3258dfxP2*HJ(vfT(77j z_QoFgy{fCA7WOAzXHc_y>JC=G&8$l$tuBx}A*!X<)`FuL3X%5Ax*FWpr*3fU-Q+L> z5d*H(t) zvwp1FS)!T%IcF{T234wp02*#a-mW|;GPoy9oZcwMk)Z!-%l5sH!j;O@**NK~q7k#N zt$sC>)^=W1^w}hH($XmY;i`jShCP>VE$0C*7c)xR{&Y_mKXn{XB}n!bsiOf8WQ{NK zxH9b_=i|l>3EikzBZ((1tiE5)ygtMAD;9rNSRgZVx4CpDs8zhxH4TksR-#$C8OTF8 zUd0voX2ON%di$e-j@~2{{ZGc6AZpZJ0?8W#AY#Iv1Y5OH*Q$7-$eo0n*)4a5sMl>| zmN=8QGNho&iNaEsn#$C}pPT8diGp?k5z6rjb1oAiBMD^UKS((GD*EaL6RAVGsHQh*2Sh#lF+RBbkiOX@Pi&Wi$ zFj}s2k#di}mE!Za<4fHLW24N^u5)@r7cwAVrM@vNGc#cJH6L?&pDgX9R`sIR)^&7` z9!{uL)?FC9k9@|3XCk=+;0q@1Q?G9%kMItXDfL|aY816dl6wAb=bf)VHOV`$?`NU! zbiOC`6-V&AZpJoRY3W*RKDATe8AWjF&Bmr`fFm4SYMllLFC&O$N#(Isr27&7X7yj#K#jzI1UZc<6M!Rdo`@8Deamz$ zo1V!Ky(92z{6!=auH5~9IpBA{f32LJf4Sgjcm^YGG~jc_Ae-5wNIz4KR@gVo0ctVX z`UU9#J6P@?+k5!w?T0%>A^K&6E|IQ@Y*OVxs&=%{9JSs>Ri=_X-D58KIVqpZkGl7i z*iwQu&I5?X`%>4)Gv(xue2_Dr97XqzLdpws4@YZG&6gyLH8`~pCL~#CwFyPpc!7U{ zN`KQi7w70%1WE{KVehq-O{(RceJ;Vc91wtl6`E4~i?EACj{d?U`k%wxOB!XuR_wnB z!%GN5jhVyCdPw2tqBZm`M3Hj2r1rmFJ6!V_3Hw))%cd{J_8Gh=@jtzqsX}vb@d)Wm zKgd-0u9DhI9jrK%@op^Z;%$Q#HI#MfG0zBJm{^1tomk0RIZkd1{<7Y9KshZWt^g{< z68Y_iidV*E^pQcP%2%|O!n9I)&2q{wH2PTbJ1XkEhMq?WCCD2rA1Mpx2AGnTd{=qs zSukvDVA5B+MH zIwfw4&lI2v@ydYp<@?L0B|4cxMaXA{L$rVpNyBll6W&7TsD{ouUj;7QNW?G$e>nq`zoSd>yj3FNe=W2tWn!~|5=G;xve`ci)AgD#CLLbz>XW3k%KVhPES zU3ryIB5ObKS?jw2Oj&H#5;&6(=`bX9O~HL}bHj5+RDFWkQ3<}}XJgX2a88E_cCesJ zxv;P9Je#6-3@UW3_N|{dO2EnvA0o>wNgPzp(%FszTC(hI9h9^*5sLY^k10RsP6fBg z$?U6PjhnOE96(jW{mSV$AY;IKx^;`Zw{*d*-%G8j@CM^g7H6cquD15TL$VE(hiPS4 z3QtKfAyM0WdQFdNV!+KMhG~MP_+*$6>=3n7a#Af>QtJTAYn6x5Pt|~~E{+&1K0QL| zyQlTkwy-h4DXJ%nVptQ3l1_vzGP^2cf%|EXwizY4CI2Huk~L90htDMg4LIc$>Isg@U>q`S`0R4~YJ9ofWO3)Bm6;|^F1ziyZVIUzH-cC2 zEky&Cv{=_;ReJ<%VNWYP6MKnS+2ScwRRNlN{8B!2Z+?O#5fOzm-5#m2@yF+%UMsk? zAs_RnFKNWA%X`THyyrx)dWPV9)HDs<(XQ9b{-AJ0JXLgZ72zQAACihQ< z-DkbySpV3L6vRXwD`RU*|7n5V=*snVErUc3Mp^k?_QnvZyiQnOHq{IxPjV)YYlh3B znb&O+ZE%~3x8s=QRQ17J+;60cR;%8&0V(GV9?#3b=s^VqrpuQ{mwPC+qvBKl6N(9` z#Uu6IPc0i0Qnc4sKO>zD7FrLx1i731v3oUo*`ww}=QQa4#-Zv%AcQIMDI?F&x=MP z7e65Mqr!d<^l4`4c1*V>dVfOMCuTbC$A_1`=Gke4<U9YISob|J)pa-^Ggm{^Nf01E{8LoqBmZA{9ms{(;m&VR}uS@)!rYxr+8;HqHb z9FT4wCh@ye5xR+#Xn%G6fhil5wjVmA4*GC^Rk&_B9gGyBT%%!; zu}Y^=U;gxY*+I6KV368Q8qEJVK+%ZfoY{kvX#7MC9uyg4+xA#WVf#9BuSEj`X z^e(G+c~{@rpMVPSb_ezVPoC(MsJLgjj#vW9VsS>dFnqW##q)>a(DZMRB%jU9w@6;1 zt9$QZnLc9e50$siDf?;b!!Rg$Vb`BGYveayzc2rJbz>IDEL#@+!(w4zul%uM-Q%*} zTaVZoaAvF+epa>=OaGzy`2NPj;bDS`CcsahsOFo!b8+Nc*5wH0p=x|@ z9Q(x^Vj>Ri;@bqj*VZ5M%g?NZoVRRzVV>&K^9)Q4QTmaEub454Yy?a-m*p94s;B>U zlCY%)C=u`S*+9>KX#KX+8Zm6Pd-39+K$#8e0ogyGu^I6RLtEDy|N!2I%^^g4} zCT$gzmltg@9m57Sw{n7Jw@H+xwGm7mp8(V1sesLbd<6S%{|md~C;vun#|E4Vw2);N zgb{zgtyWfwNEIb7^|)Q%A5GSkwEb!cXk7sRdd$@?DDOjUYIAIPKEnRdx6PGl-__0# zr$7@I^h1>F*IiyL3R(aO%^YnlkItA|BDIEjX=_A7zA+E@s(^#NI)nke9) z(BqOfAJs?~lc$bCZI>?%KZqjnDlzX5+c}DqX|p;fN{c-e5t*mO0EtPvCBUf~fsR~! z%U80<=UTIk%knqGrF67Q!}1?ECpXUrBd~`qC)#LPaMzvy3@h_U>Aj!bmd|Tt&4F+K zNk`)%)r|kleGbA zocB4nMIMa37(nE6ec{r^4G)QbCHmOK|z-E zix)cr>p)KPTU4lskOG)w-Q%)ArK%IV0RS$_3FL#o~f(wl5(lG!cU+Mtw9-Q0o{ z<_XP6H@O-&wDRP|DnOCovP5nys*z*JwTc3@uTpo06fxJK`~A>oqD?4D3*A7PAd+ zXaLimu9gkU2oxgt`6GDZQ9IX9j-QrAV)-6+%@1C{nY)zoPAW^C9MMOhUbmMfHe-8Vg*!Y*c_e=xC=aYPv}45p?MOmO0J?4A2e zJHF6|OKH9!%j0d8cplJLu_N#vr+3P_hx4u@d}KByh+A11+F@3V+_dEUx9F^-ND&~M z`pNocsgfwDC8+WRM(_+}HFNfd#PQY=$wWJSkJKz3us1FPV zuEF_>x%oYo23h$__so>x;%|aAVo+`*lZXzt|7x$OY?Z~43=sr$+?ICoGU@QjYK6vf zw37tz0!#ng$!L)#fk%Y(Anpa1-@(s{8!-Yi*}-d4{b61%FI3|mMF(0F0ZNYTQ);q? ze0W75OH@KInB;P{f`i8b64ebfpcX<-ol`HOJ$R3RZXOCTL#J;*+(8QF}&v`e;bVm@3jSt|8mIjLMo!tQnn2?R;o&>l0=9h<8G4S zH)FVPW~HpzjBoMN=|+Ezi%R8C3rAw@IDRj>l#MHj8>^~2TI@{gmXZaI&!GU@v>fY+ zq~LQ{pcW6`{`kdxO-f!+TZxQ7W%K1*pw*1%+F$QxLb^-@ z)=gkia$g6g8Zh}qdiQzm8gJ`?$yB7-<-_bN?&Em6d=)%*dJ;xSZDnQBDA}LRx%W=9 zE4B=V%g(l4;IE{E5FPVIBS_wm2G(=Qas&N2ks|JIw-PbymkP#S$mW;y+KZ0a7gz1& zGJGSfBF5pfX|z%aHv}mu+G!HsQUq=>S|duNWK*_a^q_{sIc=Ul@Q_2QH5>{O$f|>Ubj=&$S5d~TyXo1?)K|$1oN0w^OjZY{UlcfUJ-kdS5(&w6+eQ5mlplFeX-QGB zkP#E(!)r!oYdaeGq(tOt5t~#Aq+bh6z`Ec9OzCWM_x9rZYU(us;Z5MRaB$={V`1Ae zfkG448Axv$XTRv$|)2% zL?|<65+21J)=&5dXBcrQQ&tn>q2zd1f@r>$N~ej_K&Ws;Ar1e{IorT!`)q=V}jG zm=|2OHc{RCL9`u-t*IKeKWbh^Pkbm-t&2M>{ODZ2ZvVN-sAOyxprBnL{pmZ(uOmbG z;34<}xbMjT&WAK>OF58pUX`k2^!h;DXC1|;O4;c?Vf;R1!nNvVS*^0U!x+Kx8P6Cl zp82rVSNBsc=_l?#l`=}T*y+m)t5gl2n+6gfs8nio&7t)@8GvsZFw!XjJwi#%R0s8u zS4y)iUH~ACD&n%#8Bw_tAbS?Q8O%c1TXI`$Fj#jB*@DtUNt~Gx-0)-;_w_i`A(di z>Jd@Fib4wqVSu)CZ1`?2XLv%(xnw;cMzO@=$<7(OH>NGxBFMnhp&!Q=zF^06`|)x- zc_G$V9c4cfWY4s6&$Yf?Jy_`3%FHA(1n}`ntg9`)_UJy5cEQjjBa7rM0^fGmd2O9I z_sYD_?Q+2pnqrGw^2Mi4ibp&qgnoSb{#t3SqtmYU>N_Ly$4;AfYl;`S=i3Hs`c}h0 zx0Sv``&>)QsCsk%b~I#ckIakz4k%Xz$WemTOvKgXi<{bSp6Fh@%mX#P|D2I%`tb#o z=iJrh$g&qG+lw9A?}tnKY5B#)sU-ztx^tx#{T10S@1))%rad@^&y|SNxad~$L|c{M ztpB@G&*{?M#V6$y@QK5YDDz{o!0&*y@5J52Qa?L?)oNg7TCcmRjV7!27!%t7xW6P67CRs?>NPYmFg8TJx03Q6c!zS zsCCagK5#Zse+!YOk0~OG-;A7S>jXJO@|~VLA+jKkPXF2Ln_OI+WQ_B7WSr&-f7HDD zkUz)wO_&|va7{cV=DW-RaC{+mAXbWH_B}7HAV9O~|8c7g!xore@8UPNG~RY*CCLtY z?r9HypR8j3%leOFpA*I^vzwp)_q=G@xsbW2`4)Bc#!7@C;hMYh_58)?_!gz!zBOG9 z9Cy;-K!cr$x{mI<-@eO)f7vFN%?xXGTW1@C5ecohL)(Vt9r0SXC_w<#GsLG}`vs`g zp4T3I_jSPy{9IMorl{-qt|8aIn5ckDNO96sF^}@W1+3AlB_s?+g80^6;&ieOpg1AXeXoZo*dIKRfLF z@=!-ji?+d*4+;36%E;??Lw?0Ln`yq>2U-R23Z-c@;S>CURYBJz$#O7Xmy8{Uhk1?WeJzI)p#%83Wxq5#JH|JNrn z@z=}>p>H|GoLVwgZ!g8Pi~X)t@5nkGi}To zN`zjv#pUK-A^;gZzIj(tNmGcz<7_{DXsSk9;y@W2fi)9O zEq)Ifla;ZaE6yh~4~DdT|9r^Kh1aMmIbSvFX1ke#eWwI}}hK{@-AZu!3h#wGNla#9xkpA!@pF1s(T2b$B*k&`{(?qvEuI?^aev2e-6IiMd?wc--4rQ+2&TWx7S@iu1z&LFyAx!8#yGxgU$xs`gWrpqTAkv@fdcQ zTm=fYNom*6T=pN0(>M5*7-@}j$hkUk9`-h~)g)BbZ=SsaEp`h;PHyp%MWs=)%(2^8_yF z@wfzE7t^X{IT%JCHMoM?+-&Tz*+?QcKpeZuM9Aon&zQ-&G&1gefF7+c|LWy%d{wY( zi6>olmz2B1+FLX~IXkOyA(XxXmXq0~1&c#iqdyGA4?Saxg*q~qB58rF+^}dlrsoGb zCzQQ<>;3Ggmqz2IT|P<$gP99<{$sy`nJf}g^x8pd`pF7Hn3NtoH?W)cd>a7_Uoqe; z3<`A^6HwM1l-)*&3fW;VOH}?TS?^54TiIC5-w)~BC`t_tq}{GK{bsiv(HJ!Ng^!%M z+lmgFnn)tx@FpbuXZM-+ZW#6mShI%mJM?gEjVQ6|jcKG$^Wh!=zCY z1#Cstapraot1sKcHDRl!*Jq@&V)XM?g%@j(+^xwT0tTYh;17ZWyXf1c!kwCpPg+~~ zBnhNlLZ;W@D^i*`+j%VArU&Y`{r6DRdpF(>vfGk+oKO%9^itU^b~iCkVSZBF-tLFk z4l^0Ue(`vRP~UjkSm@F4IAR6OT3#)Nr&WBMLjiTAai|~Ri#&D2q8QwhxV!K+znP5J zOkZ5B6Sv#kazU}Y@O}@;~;-94<@f>#ZbbNe$Bj!?H4e}QhXN|I47EZO# z^k0rL8TDT57MHU$A75y@PcMnluUflI@| zQ`H-;JAtxubGx~sLXDM^&utddMjR&D*!W`zyP5L`0k9_*cV2}u)hHzyIol-gtTm;R z1OlaiiCG%19M_LAwXAHpyjdOdHkv*0a0p`#b^zL z=HNKa`bAvQMvbRv=3@Ex8arFAd5I$7GQ7m#rh-B_o}NlduR4kCkv0=pYk^b9rd&C< zKpz->t7K^H5fHFwbzbB{8FOhOW2jC*f`6WuSik<@7|5kNTNOj%k-xfM=~Hc$%a4Da zuA{SZROuX*9{e3HOrHnVa?+;iYo}_FAnyyBLg;`or<&&HDGG98T!6D{ce4Q-pQ>H1 z1A2-nfZBa}yMK(o%nja*o%D(&k8YWHlhbykO(ZEF-UJ3$=03Lf5xj>|RT&>yZ)w2y zNo0RqgIgLiS*t2B9$q;%!CMWLm5jXpofMdR))Q{|LA$l>m5g}phc@Wx6h7nqU3POB zaF#hCogJ`)qZR%1q>$oF4viOD=Q5}Y9;j7wb9)%dFa_*ijiC~pbFg8*1zFM6Gd;wd z;64w}Kph!{-Ymk!_Q>9Uc!Ct{K=&&dz+rGH$q%W@)u@;=%L1^Q|cmgbDG0iZCSs z(cxlSiyE>VY*OfPU)0dx8$0#v{w$I?$qi;rk}HhSXgfF8IerjhX7r!5uUGsA0i~sl zC5hTDsQrD`D71D-v@;O1ndtl{wL42;vMrBDSMT&PZaDr#@cG3|guj!F<~jVx>99b9 z0Aglvb9fat%)fUiWB&~?<9*kQdS)}YtA%7im1jX;Z@3M01$%mfNyvQq`&W(Jkst^A z_pWEVcu2^1)f?wQxZ|8pjZKCUxCaMA7vrVX5SWeqUr&Ji?h!NZ1kf7De-3A74~Hm- zt>8|(nCr9hg!A)=K}DaN{RPYjEjpUHoT&3CM6|#Y}>Z(9&LPw+5$6q*`>X;Wxg8wI6~-_5Yakwz1d*Gk(!g!K{V-6zak__ZcB$ zBc=<{zlcw)9_zah!i?9~m&j~xL9(`P%zqE3QXH8*+(NwP8Q`tl<}U7&P`1r}Ss zL(DK;Ln-<#WNO!}!s8v9{&{MUlM;}g@f~Q~6{<03W-lSDoeSJ;jU6Pnw=+ zyFE7n#RMc~L3X#tX&{i{WtY>e-3InFQIzG#G#Zk9IGmf%^=1c#_l?FZ5UNPEMbaUB z&5{urenlLOKv_dUk8a`smSv59|i~!7)4ZsZt4f1%@Ez;khs? z9m%gtSOv*Qw)#~LSez2mI>RdcxO9~se_fgu)5GUbrV2K2ij>kktexFfZkwr+MZA5U z-_5&-0L{uCwwI|Df}V;|s@+%Kd0eU>TL6v+L0=gibeM4RqO5V|@E;s^tA^eZCtiFK zBaneunQ9Z(YoW6Q1RA2#pr0J%&cmMT)bN8w;MteDz9;Na~3kWAI zxCuPajI%OrRydxM4l*U==_Xcoz`K&8!YEva1VXW<#Gr3#?8C7~+= z$;OrPhi2~`L|zV0N#V{MY<~4{e6dkAaJfnpbjer@>-|GEV9wrqIZ42D+P_Gk5vNuf z;;x=z)Z(?pFpv+_PZ`W6CQ`kMpl{Y@csDfGG69O+k_J4Akc=3z?_`3Z^iPA)sP*esGNFsVv7KuL7$tGEl1QS| z*^HQ;?DPfsY{s-WYuYqrP&P9Lcq&o1bXB&Mr}S0f)0b=m04aVr5@ugJCj`Mw6Hd;e z*Cx4Eoe0@*Z$su!CMPDIY`@SEnY*i~`2*osNO`xIKfD|nA#H5G3wl{WVDm$37)?FK zq24HlgPM{Lbs|>B0?w=XNA<4Y5H(VCyrJMRMi{ke;#0XQ5+nkkZsJDu4QCctKXMNe z8xhy~9ylfOb8Wu?a|my11}vBkWJREdRDDMxZ0-TcT60?qCf$O=b;z{xo|^8HW`aXi zA%sZUpvx3NP`c_5w-?ehxz1Nh<~U^i0N~ox10$Cjn#p|8$Ks#J#oE;Ox3rh0fkwso zrOy#WT)M;PiG^o1&KOUXYG;bSrRwt!#|i8M1*+Pr!!6wFRjLe_Nx$QQ6`-2JK)9x; zFi2=*{c2@%wuGz8O>RCk^wh1_oF=21lLuBqNcKUr5u@l8hv{W|Fo=ZeYm>_jaw}X5;%12MifOPbQ|F53{Sg0GF>?ju0EkeFEFOJ8QmL z?qF(_;M<`+nxTEQBkQ@#Eq1YY^n@E#CEY`sWe52*S{P0bVHScoLHb^dS)~He+Lc!N5~W)DT3rhO8o~2 z`rOBh6}vt?=+!-NwB;BOPz-Xo8DwafkIr@45*-=LeKd5;K7Ytvv4&OoX7|ypxQdu9fAT*7wNhm`3Sf(6^|0_X1m3s?fJ@FpMoisIVW zqS}vvenIX=iGfU(Q=^yQ$56KJip%sb=kD<-r*Ck}vjx=lq$N>L&L+l>H|3lFIH{kH!Waw4) Zfs!d@W;Vn6HV~Kvx~XTZTdwU8{$H`b4`=`Y diff --git a/vignettes/web_only/figure/condition-maps-1.png b/vignettes/web_only/figure/condition-maps-1.png index 9ae958edc6a505a379cdb9904501be7890923f4e..cde3f0e3acc199e82bb41b2d8405380e79c17167 100644 GIT binary patch literal 4674 zcmbuDdpK0v-@y0W5J~QrK{v;3CWa1bkWL~*j78;fRW}k_1eQA(->q1o6O)l0I)j%BW^S@lE^d~(!!uI*fa*8#%9piNC`(GLy-Xn z;zk1)G^9jK7;7eMK7$Y6+b~riJHe&_Y#R9=510+}*)aUhpo7mK^BFWg8%Yd4ALhes zJ`4|Ayjpgsf&-awa2A|+31$Yuc_7TD!E6RxXAQUPgnM;hK4YzqIY}7l4Do~E6~pdH z&e#~hd4g1{EV;^>I z*zbj9BUhinj zCx<^1M1Drgz7`M78(pU=&P7u#YHpezomVN^aQ%DoblA`8Hks>T&NX=LiW#B#XnUiE zA<1%}bJSb3c`2VDuS?I$si=?74ONfPQ;xZ}JA2)5)7LF2IBKNb%4p94%-#)W2vBH^ z(Cg&P?PLq=#D>)1>GO?^$++IwPCxbU#@Cvy3oydZ!f!`JE?yVREc<%Hc4 zH-jM<)_r1Qf0K^{WWrfe212DFR7YYd%2X>1QHa?{vrV5@h24~WqLFxl+dZBJZX}42 zXc-BJ#ZCvPijZ{Y^IXE2x%2xZMrP5}me-t5-oUt;VlV$bSQdCfN55b$&0s53-K5{; z>lP4>;)Mfvf*N#FnX_bzrp%+MgSs4bUtk=MVz#%(vm@(60PxC6qIt!l`qa7V(gQ*Ol+bjZZRUhGr1%A`gTxN4R`>i68r5i^S+P2%nFqSm8y>9>Z| zj-iD%{nzTC;w1oID8bRzMy98uz77p>o6xk{t_Ct?j_ldz0JqO5$tj;Ya!8C;=rMMA zDlO^+s1uw8P`tdsn?fu&rwVRGL-KT!Ed7T^tNUASOsYN7uILbk&$@f-sWx4`cum-o zI$MYhs!B((jP;?}@+qVBpcC>{SLR(F0-{gfO58QLOtD1mig0gfT0H6AxM_nJl#f~u zMxa7Z02r!I*#*B0){hZXvDfV<$x58>{JhaK7f>wMTcdLKyr0H5P1yeA>I3*KZRrM& zfgNwk;~HYT;wdiSZ!N>J(>wb!JUwo%KadBAkvxzkESGtv0N{UnIqFxnOiiOauYc5? z8X8EUm`4!n8t$K>l5Lk7Rl!$MtJ4k-!46H;Z97eHkbyJ&RGluJvhvQ^{mj_-mGg$) zv@+MlqTHFwNpJqnB~o%_QPeF`K?`#R^e*KehxjFJ$ck3s#9Xr?jp$tpDCWkKT!-Mf z9Wvraq@WM|!n_Q@l`vKTihsS!@+yibN4lsS?x=;zRu*S!aDve)oU0z!-6`QWt72a5 ziGzeX08t9eLWNpD!iR-}`KVqAuvu>egOjg1q8Ly`J6O#HN2+lu9B&&*ZGhs;LI$4z z#Gt2-oAntcvA}^WPeD83)Qs!sAiBDBCC#A`bn*w1i-bH&j=Igd?5oxVB|EtT-v@g` zs{$-y_(L21=I4L$M%NI58?BtSH<2C3w8PMHEvRh=h6-7({M$*ivX2zq7l;9ff z9XRjr3*H5|nKOs6%%1ZJ7d9!;Qzb{@s2k9{`ck#hYn`$wUrQGOyljd;buqEkW*cgo zYPKG|>1mRz^np6?-u>%>?Dkw%h=W|3GV`8zZE%^7+XFcx=IX3EUGavhV`YE)7hoGc z<2ARxURw-0@dagTrum8MFd0$ur2pk_Xh_pW#@0d;!1VT_cqu6AcF+A#+u+s6;NLw_ ztO6$YQA+}jsU&~!nY!b(jo_>`S1c5EZ_?~`33Y@vT5~Zx-TT)5Oi`eC5gwqpln`wV z%8)7KNp4UIZAW!W5kF-gLWhpN52FqCZ+||N4@wp3)5pdErEsO*iMI7JUA{FQ-h|4b}?>Q9yE0}hCbndhu;CT?6$45Ft3 zO40J_PAZ|V4=P1>mE1L^!aLmZJ?e_=6Xz`TulkCf7fE-|JdF8PRTlCER163yk)!eg zGR4Dhpok+tGb;_u_o&Xt0O8kmarQn5;BP>&h+@!dj7lxg{Dcg>AMdWmpm)4}HsSk) z!*Me}6dwx{F|Gk)$HkN)kCUV5;DvvOgyZ3j?@KyN=tkpbxJ=|< z@?UWq_y4Kpa(3tbR-WEac66k*>ine`-)GUa%yu~DzoUK6UV|L1{T{27-u)|_n_q)I zBs`JBJ92qZ3vqXRG!2KtjMwq>e8J!|aTM#7#XKA1%-a5;;mk>;)()AyE0Y!beAhq~ywl1g>72Dk8b!vn}S(#2@`9(6c@7J`W7zo@Ps49tB>55wb({oxH z?!SqO>H>qy&(0Rsi6U0*79ntcW3#L<3#a|W1>Czb{U&g_W4pNgX&UVN^+TqD<8LzY zwbl>rSZEr}_)8RNYb-Gy8FRi`t!Cv;ay7qJYQ>%J1Z4vE_XOwkYe%PBbk<2l^<>C2^;v4Q}{IL*YxdJ zIad7Fw;qTC8%aD>6v&!X0?GwnQ*iWMFJ8W_+c(j|GKuJw^o#iu?T_+A{`hcDT% zpsFUz)SR)c6!ujduWHLdGbKS@0m=&%%Ch570Q`1INRiW=3=kI{nG*792M@_qd|2K^ zjnr<0qg!srL0OxioB54{up(XBK!K`p1jqbU+tC$O9VdZspgNCeX&Lcqarm}1=7T>_vfY4W zABlmMQa*>wEmZ_Ucmy&0RF9gKxZnSA!g(t27DuJQ)OLQ{_UCU1?%~+AL>#nsSU~Mw zSA4}?!&7z0qcXtuy&iuZWJaLTRlcIf!=F2KY~sH<8ob< zv0Xc@vD6ZD+{kI1#|uyS7vl$;@|>&HoM8jc*Xe(VQ`7gZ=EpO+1u-4OY*$wx#PhXP zHU`qfErz?jJ=+;*{V|}t=@Gj`Je+Otb8G{>P;}K+CGpmv1hj-+4Z7&A$`a}LqlwgL zMsBSNr_Bb?E^m249QyQgk8zSW0m==H({5OaDQ{i(bM%%ChAm%mV=kyU`%^YHw=SIId61Zv>iw&s0}dHRB26;#e`x><_*^J zI?P|fnzs+09riM7`6U!qHz;za8w_sd$=9I%H(}^!_G!9Xi|R_Zc{tDpdfs;mMo*i+ z9?jG`O}M6l>`3QZV4qVB~YV$N<9la)5oi&&?qP}>LTFRXWuM%KC&OF;ErLBHUoJ3ER9*Lm>rut0S z@AhVDWZsHT@Y2MbsU>G3b($b?6$frha9%GJM{-`#{S1OyuZ(ff&UsRS-*sijok>sr z;^~*jAD-GB0nRt9GMrNpAo}w5cMSCi4}AHleXICCsIl+tB49h88QqLYIc3pHS8P?^ zcFJ3Xtqt&I zJ5iyVRpzJ=@WFp?-tR0dzace4E+tmCmo92}jDSP67Hhf zu;;0|^*xMx_;ACUF4dHkPggY#kx(QWPF#vT42xM>J^z+O6t60#u57)Puw;0|l(01l12qMx4E;)A|&+^VXBhlin&>rMWsgfstKw zqmTgT{oP|`HY6i`>bF!qEF=cK=eva71pSl15Fc-a${W-)Ww*Vj?f(i&C=A#}zcP># z=8BUCj`c`@`DWaX%iPP-6f3@Nu0Z0&_bw+@$5}C{7hg=A)EsBAQZM==obZl7icj;w zY%kHXP}b`}yl1AkJTQdZ0(8j|XTO`C1Hb}B>bI;={ucW6yVDoPlFgMl1(Ki;tvpFFj=qn!}}LO8$-}0khKXk7YLK(!DzFJHU`5A6EFlahCsuR2^cb~W#U-h(S{IM z9t?zlVYRF-g#9jL8iB@q#Jk=o@r8_mkTKf-bTG+G8kxzY5wvLp42{ex1R9M=W0Gl1 z=B!R9)3cNr^pF{YV-f?HWDJu`U^baE-=1ZTDX~2JRtAj88iVCvGMnmE!&uYe3^#X= zgh03i_AA@a)UDqTh`^+UiJ|lL0{SDKcp8~M$!@y?@y9;86d z?!Vl-%FMtEx7|Kq;UE>>VZjV&x6Cu)|E#OE2F=iLF*uS%VfQ&BWi0}i=Y|U6xq(%t zWkZbS&U(_%uS$JmR853m2Rk*E*G-lzce2XL}n8Yz35R ziSAzDQdjXF`dP^cR8@$HXeGa8)bA`>?dcRZ%gRDyUrIgXtG!&7`ME;fBQE0kbj$`_ zRATS3-tRou_$JZKuQ+m$daq}H&HR7|HP5zY^sw7jz}t_j18zBzBk_y%79YE$42ite znm_&hBMmwBW}YxFIT4!*=4SF&$6TpzEvbt`FsCnti%8FRB7$;I&u23quEdjfYFe5& z*I(|q-)PZ~jjGJ7olaUnIJB?j1iOv$CI>k@#5crv*EV~Vhxl=-ue;tVTso*qshQ}& zqvD&VVa|0C8_gaY52tcmCUOW~L@2-@{7P4;mEb?oQI8${F#61+>0yf4f&5|z;q>!H zHSf#qZev1|Z%8wD5;QUASqW*DFW_o`Jn;}xPfUp%I`SV1gEGVg75t8!)m0x`=4 zM$7t|lX;Rv1d@7`mgmYpkWLK{fCt~WI1SOzaU-KJ1QL4zQQQOvAzS?2r#o;kaVGt44Ict2vU!)d{|MagDcCl)mbZN^o>4v-E*V;Q zbsW!~e-OmJ7UPfM0<3^k=pgM8uJv{u&#wbUg(sj~w@*s3S4kluHMhARiq~U{Kd*K+ zOhpMRhJ1b6qGvb1v2g?((qsoRK9B#^a))G23fbxxbR;&2;e07Bw<@_QKaVHXKhJTw zg@3fTY=oB!Jpf8)-l8nR| zCcpZA8hN76gB&V&z1J1}o{P5%zhxn9x# zwL^M0A``kb;`|s@mKxZM+HUiM4kx%WJx^HwxsF~R7)#+(Ea~iWJ4si)8zqZT zpv$Xu?!^-Gfj_?an6GZ5_kY(c%#C-8!^y)t2;`}Scs}u3foANDuO;6y)aD~Z9*~Bo zOQwAUF6hn`Ec8hC${)kpoSuP|EQ48fPSE_nI$e)5Y-K5=7=CDeph;13$=2~kRVY`| z(heymz^)z9nf+S6S5MeJ%79wn`jp>Q56Q-vv^!>c%XP1&Qj(@2u6?zRLXQ5X7v8!( zQ#UDAAQ?R%=fl&p68@==2NXI0z%|E@VNV_erFz3kp10XYA@Qpwj#k?stafF4PT(zX zHkY95PO{X^$Yi(>M_a6oA^Hflyu9?aKn~>PU+w1^14Q2&WB#qhb4L2IB6$LQ-z7SA} zBzF>*xK@R65#50ua3XE>M0|NK6D!7s4V5{cKKXUn7QL&4waI8cZ_}p55W(6w8l-E& zzRL>QB+sM@`FfUoPpBJTT#rL3wwmn2En-y@y50nloS0|DNdsq?FWTcgYg?MghFkiq%?HetXiXXj{ z-c4LisHa-*c7GIJ>zVE(_cJpluMGAmtTDck1B{8VXxTK?XgYd_N+Dra0q2dJxs$SUI_m%0kLm&U%^PAc5Dniy3vC?|FHp>n~YqJ!LkDlTHO zL?=!#VOc6wsv0r}8_q(c&OU)cQ{<30njMmGB^gp==IwEqG9cM~wa=j$teUGO_b6&p zSL{1SY@sGHg3f!~%0?nI(>OVsRKyqaX;%CWNkH`YQ>kj7Lz_ym8ls~WcWq4qydQpVr~4D*-0#=y z9Vrj(I+`7ua?Enk3Q1mRm8wQcb17pqG_po=0!<_YQieDHPuxSH$3r`UhwsT(zx;LN zWAWKHCziOx`JS~N)eJrrai26vn{n_?zf{Kv3}HQp&p7$yT&h!ubP&&?R-mn|?Qd)) zEX})TR3T?Kkq3%FOJ2)&ai0Bw3v57{nq!QQIZ{YW=kBz0_2*;CO#yWxNL>QJQJ+t8 z{IHGB+!ac`0mAe)79DOLmOZrW14Ami5`rRRDS&4G;tm$Gab?H22gGe>#^%{H#n5>+Zy ztW&fuw=rMd(ra`5pp@lf(IZfPXSbery|aBFu*{J@^JTCpAZo+dWp7JxXdoz8QY8zY zt}@~nytFx!_@^58eBql%+|1&ip#P$Sn}96h zs()^9+&TX06i@PEeNq2iSIxS-)nvnRQsY$a0>6!14+{i?9|cLqtryu+JG<6v=N0zD3SjQLm#d1jn)&gZ{1E@ozz%J=Pv_T+7(#y)c$;-c1M*=T(*E&2Xa9=@@)m?Njw-J2B|b z9$|q*68r@YyCU^ZiaH17l6ErBQelamdBeU~nWxsegFYb-j`1U)Tnfe}-zBiPPD%cl z!vL30`rRKKEcssyb>?Sp>G)P(ksF3yc#ws8)%?Zhpq?;alc*T_iDjPD&<QbsZL`pS3+7DsLa2MM5lYlrVC!E*CVHdzmAMZ1D5Wsd8XjE z-Ash_J4XpG$_*$H9h)6bP4=rm zKxvov;6gbY?MXZ?$*;gy>kdE6bLY>dx}f9D9Kuk!$)CS@7GibOB}8#k#N|X%rTiq+ zF&cPJ<-yuC0m|ty3YlKSYz|;JjWEvm7#j=qJn)eLF#4xmrQkxm`Ka9s27_!sq(`wC zdYg1mw%~WusUK>+alU5<7IW%9ewiQdwwmV$clr#eMy!okKK~90XCQ*z4s93}zIVTC z?U8wq(na*lsP3l=U!q-1p}{qIK6kC)|AZ*q@2$K<(j%Yzyv##Z7Y=-G**l)Wb>0U{2e3DP+4du)i)Y1jVo*O51bf1eCeJqRGgcIP=_5Q|2Ci+1P;~c9TB@Lh0X0L^@>x3 zYl`JIJmWK4v3{EX6V%Bko;`KhP4d3(wV>-h$NEkxoS9@Vp{2=EucSq8l_hZ3H~ULc z3L@p#`CQO)PShC;(HdB=JAfwA_O(Qq@kDeS4>=d1{b(3(=c=I4lxDaGNwn&txG}cQIOO z*xVbYLGzFC*o?gbcah*$!1wi}w@Ay6@D+yZJ=86)OW4MYx_JLff#!a~+4KE?8di(c zcREyVcys=8G-A{4QT(cx*QA`Yiu>fof{^5|Kg---3usj8zo)Hs$#73j6z7MRv)g&z zFa|h1y;gZ@9Z{{*(^tPD)KLD-7J-bl2qUXp9Og$|ql6{BnO4X$K-cs|DLoNP30y7F z^x+0A{jRg}i9}Iq#qyw^bCI+*zr+*oznU#J9(L>S`nlwIGH^}w>YIp`O()`Up0BVE zt){$k3^}3cX)f@+{0q`0Gh<_mxMZ^tM4h4+oa}KARK zb<%_<(DyuDi0ebE#PL~%qR^aCSs8dG3~QI*b+4Iw+!ue%Z}#3k)VL`Dm`25Eh=zw$|1 z)eikS$y@i{s?Whi^%?U2K-Y3^ILll-Y;Vj$Fvq8-g!UfrhBmVSvp^OPmyA{dVQ4^+ zSR7_RByLpkS`*EbpKW!0V^G#{5Lzsm@|4msh(blA62?(Hjnkk8N%n1)6UkbXDWc^>nZY<1 z6*7}0)u7E%LJ~rj8bjIUeWr7MzxTZVz0c?KeD3GIp6j~q`&#bh`lcK_uvY>53l;{0 zDLC%4KMaGxxzbk_2~i|jd_2?$4>})lfLa)g2!l~z#4H$v0;BM-Sge_u8AKw&h(scV zNX#N;^}0x(U>BDRm*8NRV75y}E9IvQBtRj;D8wuZPXhf8cF73NXbon6%8;^VWCXJ_ zGFmgBZmT4V2+NY{qGa_-*ugIBPZAyl#-k8<5*|+?X=!O`&2VYW2yTT0S|zO#Nmo}_ zZx*pv!t0etKCy#8u`@nNBoh-8keH-3L(-cCQ6&;c6ylY#NFI{Cl1cOO;j@PSAI(+Ty@#JTg{}pwNj+(g zX(R0-ggomI{>A&~HK3Ep^X?9#j@81zJ;!asMz1mi6jk)L$O0MsCz+I=JaE4&8UQ64 zYDQt`Sh*f;peg$HtfFkswnhzMzri_p+jW1~$TwvG?w2KHTmTQ4!c_otbUWC!#~2qd zRY%+MJR2j5l4H?ccF5x^t-^uATv{m7tj0+?BxFE>F_1%w!^B*@wTioW2xVEdMg$^uN6#{vlrP56OHcSvBwgZ#Zxmw9%h`+qgUzX zwFwE))0e)#u&u8R6kbcr9GHF-?dzVKo9l6=;R)x@uO|`-%?U;d?`G@@(zhA`J%jFT zZuQoT<(Tek%WwMqtk-RA1MYNiLo>VDsvd~Vgt2@@)4XG(9gS-e0}{Jm1Sy%XNfL4j zkM}FgeP1fuQ8Ir#{Ybs%DUGL{0=}I+-0WFU7)UC=FnNf~47c(*szNmdBU&nz!Zc1H z5j>9qB!21EgJp5ug2f1AcEMH_mN}Mp2e|P-8srX`aBu_ohMz;*b@dJ(AU|ik-K!%F zhBVA$jC+d0u@mJf{~<`_tysPL`TOYi09N>$UH|y9z`Eyb?avTyWR3IVd=LiUtW>+J zs&hE#V>)q{Da`;gxK^SeNaeQ-$oD9NsXm_lGUol5$tEKxUREyLc0Z+P7bK!ZhVxt& zzm9vD-%1OtHPU5;{ymHr?SUJjWfiq9PpZoe<%Sq}x(aPFhoJc5Y0+m8EjVYJ@f%XJ zP=XxVkF93~*mEc?#Lrp5Sv!otJf<;RM7e^LPj#THZA4LY3=Tu-_+uNuCPeAqtPg!~ z%eLy*=O*|=1;^A>>H?&!)M(@6|3A5(GW-$+yWGDDvKRJOA4qirfMOeA($)Hw#ce;a=!R(7-x zobT>jc3&0`mOQNAxzq#YOz%kj{*QKzqN$=QF-r;aZHsfYB@gTBsxFB(JNo%78T@r? ze8%p$Pf9)J?k05y3f>3Sd-e!HRWIh0kCIB)9EeYTRJ=&uwiaXOB#hAQVU?{8JO@E{e_twB&0}p{N3WzJZ!kV=PwC0pU zu=8!ExDj1|z3aAe)mkh*=j|%zT|SQB9ojI)Z6k)1gO{Q8EO~%V7JRZLaKsR-00(bp z=$GgFjDaPqxXz$K?kUCwdh>ayrg9}pm>kJIAqY7Cyn1aB-QjO4mREDPPK8j^2rb}k zB31rT6Y_e}Orl2(R;qh9ct1O37xxcE=bFSAFv1r_xFyRPFj2#8G5~w&ml0RcBRT|l zDwL36vCRxk-9V+>FpPX@iW8qKLR>j2F^CRu=AlJVcbNLkD5@QI>p_hpD5@g92>Ix+ zXGNx}r9U8)-y$otY=VX#6Ng+Z#*F9=25?*JS+r=3lNWeSl&-cH6E}6*;Q$$ifZt#&+n4L7UzK)G@FSj$Q{>wz{=SLN-X-}1`#L$&_f zX;oo)8%;>34U1xcR1+2a5q83p7|_w*=R6JW!=`J_%i!wWMcPHL&2^ z)F}txR4PIpfqc+~sI0kfKzcm!cv7&o{ubu04ffp0B%b3Mj%>lw^6b3t!#ouvvgUz0 z=|*10f;;|8W};8ka7RPAT+*cN{aI7tWG;Ee<1O>Jf#*in${SxfPI+%dXSm=oSq|(> zHnoY5{raVee>Be;&b3()FVvXYtZWK*q$WLkiQoRCH}X%)jvu@pyBI=u{(o`p@bi>$sKyG&+&0Tl%%2MJ$ zy<*fRxuj{^`#7*nV^vF$C&F&7tG^lKpO*OoUCEA9Q3MVABjkByIraFrvrVbC#bOO) zmv+$765*>e#u&&aE$wol#=d*3>?{E9D&F4j-r`cxqnRBAcR8DxLRXtBJUK*6+iXH} zNNGbMF=1}N@5wJuAd-|*k9|}wS5Kart!gL}yb(OIhI~fW0jS9@?`#|@BCRKlS98A6 zZ_0D&QNUS91e>%v^}V5ZGV^5S&~+WxM_+UuaJatXab9w`i}QgO+Kg?gLmWQd&>X3Tv+Fm{Q`5pD?%a_KnXT(7JEYP{8S5!6n7zAaUXy1TCg9 z5?UeOFCmgN+)`Zc^NjBr7cm93xne%{^8G(Qq&~`yq^CC&!qqGOtns#Ysc=l}yZ{~y zI}14b(-$sA+ftrwkLEqmDPlJ_je3bI!UAJ5S>6?j`t1e>%G+hewc>kJTcIc?S=4-Ee zpvW2*Ea}Wj<@7ZxU66L$kl*Ey~15y0ufCq0M~Dsf4JX+qrZzKg-E!XKToqo~k7K3ng?aZFTN|n zenHwH(M%nmPOHw|b);px<>9Zzb899vx9qGC_~ljmvv>x?1BcW?&72D_e%CD8t_L&= zHPKyTXv;ck#PL4dL_)ymFYCB5qr(<=@bYMPwUe1E+{cTp(O_nW_pPCwoaLMS7pqg( zn*f;H3-Pt1cN^BQCQ88_J(IyTwIv#l=Zq?rZ!TwGSmH-ZTfwZv8`epeJgG4+RFy_d z=3p+w6CuSR#+!i2>GjO;$Y^iEi#6ZX^fm#f?BXponF`VSAp+cPsStu6%Y7D6G2iDm3ST`GhWpWV@ zL43t|_GwrC3>nGl!o!)Oi(o>uaDo2OFZPSgA+w0b{8brR!ypB zBp-%P^zK_Nw1_&W@!^zuY_Azk)(iMBGlCMFJ{7oKGbX4}X6c}g@X#m`43hUm%cts> zt+u-7B%#9gP_ao@$xVC;vS}n2IdLXKVkCII2E+2ZkPf!19pdLtWP4@HC1sUWJy;70 ztQ0xA)~?NZJqN8OL+Ju?=PFeI#oCKQ%kMmZfie)@NIzGM$En91y<-|Icv7Z>4moAU z%_8AQH)$KT%-A9czcz0JJg@%{!3 z(Irpa%e6SuS8rS21E;+{vPG~gW)bSGma&TT=vzTxoqf<&QlQo0+F0z! zm!TBqgVddVac)!;(#P;Sf6fP3^xBQ&g#5yq&#n8R3V0LLg#~2FV!g>nA{bwM047xS z)^;w%qJ{o{fg5HNx8_ewPhnI|co>Zx^LrfW7T@|WS%x$)ad-AyKCWvf&)LQ^2EWHC zv-AzX?3NQ9hIIIJHeWogV;4eGNReaJwt4YYm$VkZvhE%tD_&Wn%mzgXKEi(--K45D z%`R$?Pc_zHSu|J8EqXRBfle3EW`&ax0pv^aoXoS6Z?&nPS33ZHGhbH8Qm%7pB{mGY zcT#cB9x%q_Df9A&o^%Li9q0JFoJULq*===Yf(=WVOO5VO711x7@vDL|hE(32JArY#fFxa&^=L=G zF2b+o?!La6q=^a(KFU`|R1Sn_>*;`s6)@5mTvQ#c14Q=d%}H!bQXJ>E2tErkf~qwVM)cS%f_2 zJivA@($s?qamfoCCHu!^`V4G&h z`lvv$_yjZA&i!zKrGx7{?my(wZi@BLllXsss;a_mt*cP6^yY^8u@s%PYR)A5MQlSs zCLHF!5~;J=^+5Vxn^V`PmA0 z!I41;4D3c5{`+is+7`&Xoc8ZJ+qz#Cf5aAui+xzC?_n>2O6-}s+Iwv+dD8P@)|9o} z8Hqny-MaDhgj{r)6Uv&Z+;S|jwP*YhKT{nS5hZ@~D83g+-!W;u;;+4^+2y)Lbj;Cp z*=d6M4riCo`qs{Ed`{7zbN~ca9|{c)J9AgAckqV3&k^8r-|NIq1IcnELQjJF-$#IF z^Y`0LXJB^IRXXXtdJQO+6Gpz~EmUUb+G4hG$68+K^lY~2ODV-bAv%-r=hk`@q38al zZ-L-7QjN+esw;n4jmw&oC7n!opV9`ijL6quk+B}wp8FN&6{)$?-obU6>UR?nMXv;_ zg4*TOLrJ(D!P#}_r1Hc#hJo~UeQi3DoC033VDTzgkc2Oqz_a#dTt=uXK_i~@;Xgar z-+4di+wuY7>|;c2!x-m?v=8%uVhN6Ui@PmmfTm;=y z7h2X@axF+lmWd>%SQKGrBD7#4KYQ4E6JUT4}BZ_oz;5PIBD=E!+a^^F3xb z8NqW{!3n+?-lN^A@o$yIA044@8Ys}u#&A#voI{cO@A~r1g~c?@#bCI%Ss}LQKBKvC z^1!$)uZ||y!&bh$E7$fus99%5y;--?Vt8Fsy`O&-@&*KHXc6x0 z68;1?)Nyf%3_b)G$m46Ih!6w>EmaL5TJ)sYR=(iUFHP2_;QRleT#&5I*bWomb4SL$ zpAb^0RKUucTZR5A#gKk#V0|jqgb!W*2c3eNA90{b_=kTh`4?~osrJw7i8>=b5-0HRm5qoN&fhqolLvR<@_csfoNQR#Q}4vc zmawl2z0BdX$>|BD8XFT*-+2fH3x7W;j>@iB`IE5nc)a~p3c02^L&2MqBl0LXl$*|tpz%EH!1a7aq^|$I1&ztFaE@m#8IWJZ+`XZ)1 zMu|>G%U7)SrWlK#PXH7LK1&Vm=YdM|TK%(g*$uNeRM>IN)F!by-0sT8VnvVTujSqD z;=$*oc@-p!XzuvXK!Rj8bKq1}YU;JgnvuC7sQgLwA7gWpP)!(`J8I#Cv+D%9oEH+2cbE$Ck#znNOd;JjxC#wHJT>tBKDf}u?ii7V~{p*2Du@uX8bU0vN Iwi{&p4?9Q?EC2ui literal 5860 zcmai2cT`i&woeEgT4>UO(rox42nUefR1igj0Z|Epa1te0LQw>X^h6P*UQvi5pm@bk z0urRyC^=HqP^1~9$tMIM6s1bdO}O8C_q{*fdS|V(X3d`2d*J~p`j_or70z_Nubcg zZ{qX2ySsbRta|trJ$(KUBXEe3GQ{Wecszj`zbS>^lO~Yz`2tnJT;UWxpI?sD7!-Vs zWbi@H5ClScgYZSR_eHND5X#lgj(d-VKbYu~j_9dS2B6|oqXEgchue2KXJ_?nZ&9(x z?&^~KZgJ>0pWhL#y!RWz0wF1Kz$!easbW?^^rw{e_AhMb`s>Fmo%;|5(bwSwWnk+G zOcMT=nBH6C$aacti49HREq=uJ=iY@5yCOHWzU zMr)ijn~+f(+tAX1P)vUjxq59{2~HweqR)=oxPct~8_4AB9Y}z(m$)wadu=h0lL%vy z_C;?7>z17cH6YN^K@~gTzitr09;%1XPTH+d-Si?&tqPL_OAP7K&RI7T_M9Fm6)z1p z!qA3u{XrrH*qYI*!RFROovh9(F8Huy2JNI@ql7ucu0heTZXiC_h&quIFnOxWx2Ecg z&7aTBvzk3u`}z_cAyc(jv$xS#DJGmbN@&!Hm?}$#K`XDg=G25w{(V9{*s4_aWv=6+ zS7u!nFjl+&boa{Js28vKB?U)txbmT}L{I0cj}?aXVZIy7R^{Bn4n1N)Z&gikbJeKj zi2~DzF1tM2ty`Kwqin)=g`(VZA-o;*m$4H$<~Nwp@qKH({h!;Oh;PC&c83W!uPPdnNQqN^CzhYBvyKA-bkg(ZzI|iEx5hvHWl;>> z`Gi&+Lx-G^&LL1@P36(R8CfK`X-Ehe={RKO6Uhr{uCO-&*eaTz^KT2UspZPXn6&rO z*l^nWzkN3eOlDm}P!swc;LCXxEO4usloH9@ZNR`~ajGwY%9Kj4R0 z5;QI9K$8(~`!A0QBcTYkeUw8-+@_w#Ra2Ir?KH7Cf{&#MXC+b-VybEu4p4+~5$35N z?Sg#>7&DR9}2kxI4krOG#UjGF!j$1_LXX4YQe z#SU7#9{RUyD?^UT_#wf}OTwM*(eQoN?cJL}EFeHPo zvfNlT4@ZDzERGOFM*{QLb-FTkfe8K)q76-O*E8xWI3x6@jssXM+Z55deo7%xI&^c) zUqN|&6Yp(u)!8q${AJCCbm0{%@K`(bb#cxR88^A;&3Y)w73=O~LBC+K-WJ;4dw%fv z)_l{m$EYXj5TJdDIwEpL7Zw#9h$+~eppmW#ds*s$zeFuBN!9svBu_j__myWy zb-TSkw;*~N{G-&0p6~YniiYxn*G8pNm8W^N#UpPJx}RbA~-WweKG9EAiXaUmcsw%a2PK$?V)M41>a>>r;!9))W%+VwmL;e!k z4bf#eQZ$MHLkWAKZligzo~yS6NBOQSoV3nm<~zH0=nOUb?CTH4P|kspUie|)wK$|7 zuM*R8V%PPR+AT5?RT)R}>wA-xObXLRR4eZdn3NR_T)3q`C+s+>qs+OM2(u!OK(|zB zdl!$h*r@Gs%6}BggN7xK#Y=G|3fSS3%KAgNONzt1{vY?cSaq6cp0XG=dRqvnpvY_ye&FQJ1r=bs98OUoY+0Pj-%s*fIXk)AA;} zNiNK1nDH)~GGn^C{+akyo_(s%jb^Y>=`cW=-hu#ZSCyG#9q$+T5wl+h%eP^wQIkhU zp^6t)^ut)&I*1i4HonIRvM>Uq5w%wkZUw7t-|{(QkgQP@ zk|2wcyR3IP4@Y6CRr8pO8cwRd-~c07+$V27n3bmClvyzArML>O2w=Tp)C1CWVWwb? zBzim7`m{)KW5Q&4uWtu*1l9}-T8n-ug?Ur^;{`x()%gUsiE-l0G1av{j^2&FkTmcsBx_x;>my?2(u|k?$(N+8e9ZgtKM^Z zn5q1`Vf-DGeMZHhF zMP}_9Wi`AVi5>LKD2s4(DY1d&HN>pGrN5*qx0m(>H85h1r3TfzQ zT$m_33S{8yP1FU)y)LZEDkSCsfZquqB&@v+nOmtqwwRY&Sk4w5%`Dneejpb5RvhRx4LgO#P@Hru7CXAM3e&aj@h zGI*s+yk^f(O8*)&;#@A_OYA7|>nE~QyaEmG)jD{{+E9H{3+h5)3olDfHQ95fb|Srb zHT8Jt)vM?6y2&1fmwdX<4(54lJfnMWbFLw8DDXR}+_dXSl=6>Vr%-_4MDMieT9z$b zd$e|Gxgg3h!YK(y{@JTDil*_$OAmsz-tO3J0+dZ+_RVnRyIYxSH^QC|%jw)yTRZf5 z6IPUZvx^m&TjAJeRJi)cF0DLqn=U9ou)Xb`gk2(B`ShNkr8A$*s1-3=_NmU)gvWp{`UE#Z_@U;xB`3w}t8!=!bi{6n>9jx01b}D62)y8bR!IvB~ z*g1M=d-L7V9O+iPRZjREH1>t|ZXcI*x}#+!uTs$>}3b?DRTfmH{V4CCQqV$*_A4vnxfQEMf?Y zW&Sl3Egi)WEvyT;RJtFPadnGk+Jqm8IE1lH`O*| zZuq@tLrXFgP1L6Q31a{Fht@%=ks4NY2GF}3WsG_Ajn#1A267IC(z^rCO9COotE`67 zR}^UYa;`T0b+a2wG!hpERk&9~xxq?CdKeoIthFC_MboH$H&{lO#AXxI=X*Trenb1t z$mHAg#NCkvBAH^x4U?eZGr7wL2Xcm>WzAuUe4C!Rv{}6A>kuU9cWb*ymH~VTB@qjc zT&{99t-1o>B`5E|u9QY48r(_KIu<1kOl|a5R6}RViMY!~yUv+0+hs}TLgeXBF01yr z#rRAW?BBe0ZX{4YGUT zvNW!vC60zSuPxO3>#xIG2w=F&o&B9l6VN*p+1~EC1O76&`H$H5?%4Mn)`vwXT{lA# zzgbqL@BtDiyr~vHk5DEfEhCp1Zs$6`l({_F@ij5LcIiQ_F*vC;RBef^0HATcZG zEFDb&jk!SF2O`(P$^Cu?E`{=&icqAFJsV#VHSH4W#oPsZ zn0ZfgC#&EGSEgPfz|Y1XAO}e?lbaNpe%RAd!?|s!~)0N%;@6xtawa9&>vn{c$teOeOTzW>~vK z7TqD?tZzwwD3fnlCm!-V;S=ntVP+OMY(vu(ciRW$J&-;zFo^69Y=v?ti=q!7i0S+k zT;Z<^IF>QJ1uFDSY(a#SQz*=7ZC2fR=&|4%xT=6dg}-P>ae^BxYRPdkk*#A8owN`gL*bj^9L&T$Lh*k}b|&^u>XANI30GU7*S_Xn8u9~0j6qhNId?!Lh- z7U%zAh&x@~AoJ>;Qp2ehr!Tw;R;miRBH~v`b<_ZuY(BAd-oIC&1svXvwHK*2L0^<8 zZC?K`Y27pjcGTvC_1`|Ln<9&44P`&FPBX{_3e{P?0VMFrK_WVSYR{XMgHNt{@-oe} zCAuq@BgP@02?Fz5=?n18 z{g|QC`;ur=Zv5!2fS0vkOOC7OuVBGq)g0*ihT&8NbmV3SDXjVq&Sn2-4Q&~fc1=LK zxz^|T(g7pD%hFzv{KG;}CBTWszi&k9sgXa;H>m-%4{lUSwW=y7&VtgG2xt7*^B_{5 zm5{LC9i8QP9C@jcp1u3G|VuwBj;|&b|)I)Z>=a$M{!31 z>7Z7-S}6FQctvEy{5=b3`_Bd?Q)cXvnm)l!6l)n~!nwcb3hq^A^ZO4PPFFp zt!a`XA(@jzn?#SS2({-Fr3ibQT9pqt``>K;_}5p0|0w{d!p9Ck`keePVfi6l3}v6& z1g#eaP70Z}fR*}8n$r}Tgmnl%=0m_uP|_;@AJwYEt2$!9>Q^_Fbz<5nB&IkO{ldMo z!phJ(Q*EGF)poH)O5+aEt&sYYHH;MsC$1#c`Ffw|<#n8JdN4)<2&J0IKYOSY}8j*zJ3 zq~fm!jITCmc=gMi2rcD}3|a$+?_Z&Hz04db>I~j3ZBUZg{{7j2gb@FB=8(dL_ZP`SWu0Mm}=#;kn*XwU!i z(~%?RD+aQ9T4R Date: Fri, 26 Apr 2024 18:57:19 -0700 Subject: [PATCH 03/19] fix bug in integrate_output(.) --- R/fit.R | 7 +++++-- R/predict.R | 8 ++++---- R/utility.R | 3 ++- tests/testthat/test-deviance-residuals.R | 10 +++++++++- vignettes/web_only/VAST.Rmd.orig | 8 ++++++-- 5 files changed, 26 insertions(+), 10 deletions(-) diff --git a/R/fit.R b/R/fit.R index 4c1ce94..297d306 100644 --- a/R/fit.R +++ b/R/fit.R @@ -603,6 +603,9 @@ function( formula, tmb_map$log_kappa = factor(NA) } + # + tmb_random = c("gamma_k","epsilon_stc","omega_sc","gamma2_k","epsilon2_stc","omega2_sc") + ############## # Fit model ############## @@ -616,7 +619,7 @@ function( formula, obj = MakeADFun( data = tmb_data, parameters = tmb_par, map = tmb_map, - random = c("gamma_k","epsilon_stc","omega_sc","gamma2_k","epsilon2_stc","omega2_sc"), + random = tmb_random, DLL = "tinyVAST", profile = control$profile, silent = control$silent ) # @@ -689,7 +692,7 @@ function( formula, opt = opt, rep = obj$report(obj$env$last.par.best), sdrep = sdrep, - tmb_inputs = list(tmb_data=tmb_data, tmb_par=tmb_par, tmb_map=tmb_map), + tmb_inputs = list(tmb_data=tmb_data, tmb_par=tmb_par, tmb_map=tmb_map, tmb_random=tmb_random), call = match.call(), spatial_graph = spatial_graph, run_time = Sys.time() - start_time, diff --git a/R/predict.R b/R/predict.R index 1fffdef..7288a92 100644 --- a/R/predict.R +++ b/R/predict.R @@ -39,7 +39,7 @@ function( object, newobj = MakeADFun( data = tmb_data2, parameters = object$internal$parlist, map = object$tmb_inputs$tmb_map, - random = c("gamma_k","epsilon_stc","omega_sc","gamma2_k","epsilon2_stc","omega2_sc"), + random = object$tmb_inputs$tmb_random, profile = object$internal$control$profile, DLL = "tinyVAST" ) newobj$env$beSilent() @@ -177,8 +177,8 @@ function( object, if(missing(newdata)) newdata = object$data # Build new .. object$data must be same as used for fitting to get SEs / skewness of random effects tmb_data2 = add_predictions( object = object, - newdata = newdata, - remove_origdata = isFALSE(apply.epsilon) & isFALSE(bias.correct) ) + newdata = newdata ) # , + # remove_origdata = isFALSE(apply.epsilon) & isFALSE(bias.correct) ) # Area-expanded sum if(missing(W_gz)){ @@ -217,7 +217,7 @@ function( object, newobj = MakeADFun( data = tmb_data2, parameters = tmb_par2, map = object$tmb_inputs$tmb_map, - random = c("gamma_k","epsilon_stc","omega_sc"), + random = object$tmb_inputs$tmb_random, DLL = "tinyVAST", intern = intern, inner.control = inner.control, diff --git a/R/utility.R b/R/utility.R index 9d3c28a..6ecae89 100644 --- a/R/utility.R +++ b/R/utility.R @@ -190,7 +190,7 @@ function( x, # Ensure that last.par and last.par.best are right nll_new = obj$fn(x$opt$par) - if( isFALSE(nll_new==x$opt$obj) ){ + if( abs(nll_new-x$opt$obj) > 0.01 ){ stop("Model fit is not identical to recorded value: Something is not working as expected") } @@ -206,6 +206,7 @@ function( x, } } + x$obj = obj return(x) } diff --git a/tests/testthat/test-deviance-residuals.R b/tests/testthat/test-deviance-residuals.R index 1384e59..76a880f 100644 --- a/tests/testthat/test-deviance-residuals.R +++ b/tests/testthat/test-deviance-residuals.R @@ -105,7 +105,15 @@ test_that("deviance residuals for poisson works", { data = data.frame(x=x, y=y), family = poisson("log") ) resid1 = residuals( mytiny, type="deviance" ) - + # resid1^2 + # mu = mytiny$rep$mu_i; + # 2*y*log(y/mu) - (y-mu) + sign(y - mu) * pow(2*(y*log((1e-10 + y)/mu) - (y-mu)), 0.5) + # resid1 + # pow = function(a,b) a^b + # Type = function(z)z + # sign(y - mu) * pow(2*(y*log((1e-10 + y)/mu) - (y-mu)), 0.5) + # myglm = glm( y ~ 1 + x, data = data.frame(x=x, y=y), diff --git a/vignettes/web_only/VAST.Rmd.orig b/vignettes/web_only/VAST.Rmd.orig index 646a5df..87717c2 100644 --- a/vignettes/web_only/VAST.Rmd.orig +++ b/vignettes/web_only/VAST.Rmd.orig @@ -226,7 +226,7 @@ mydelta2 = tinyVAST( data = Data, family = delta_lognormal(type="poisson-link"), spatial_graph = mesh ) -mydelta +mydelta2 ``` And we can again use the `DHARMa` package to visualize simulation residuals: @@ -249,7 +249,7 @@ knitr::kable( c("Tweedie"=AIC(mytinyVAST),"delta-lognormal"=AIC(mydelta2)), digi # Bivariate spatio-temporal autoregressive model -We next highlight how to specify a bivariate spatio-temporal model with a cross-laggged (vector autoregressive) interaction. +We next highlight how to specify a bivariate spatio-temporal model with a cross-laggged (vector autoregressive) interaction. We first simulate artificial data for the sake of demonstration: ```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} # Simulate settings @@ -273,7 +273,11 @@ for( t in seq_len(n_t) ){ # Shape into longform data-frame and add error Data = data.frame( expand.grid(time=1:n_t, x=1:n_x, y=1:n_y, "var"=c("d1","d2")), z=exp(as.vector(d))) Data$n = tweedie::rtweedie( n=nrow(Data), mu=Data$z, phi=0.5, power=1.5 ) +``` + +We next set up inputs and run the model: +```{r, eval=TRUE, echo=TRUE, message=FALSE, fig.width=6, fig.height=6} # make mesh mesh = fm_mesh_2d( Data[,c('x','y')] ) From c30feee9972cafd926fc74303b53957238ca1463 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 29 Apr 2024 15:03:20 -0700 Subject: [PATCH 04/19] update `integrate_output` --- DESCRIPTION | 4 +- NAMESPACE | 2 + NEWS.md | 6 + R/fit.R | 2 +- R/predict.R | 161 +++++++++++++++-------- man/integrate_output.Rd | 81 +++++++----- tests/testthat/test-deviance-residuals.R | 2 +- 7 files changed, 160 insertions(+), 98 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index eb0b5ef..66183d6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,8 @@ Type: Package Package: tinyVAST Title: Vector Autoregressive Spatio-Temporal Model Using Minimal Feature-Set -Version: 0.5.0 -Date: 2024-04-01 +Version: 0.6.0 +Date: 2024-04-29 Authors@R: c( person(c("James", "T."), "Thorson", , "James.Thorson@noaa.gov", role = c("aut", "cre"), diff --git a/NAMESPACE b/NAMESPACE index 696f40c..06c81f2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -37,6 +37,8 @@ importFrom(TMB,sdreport) importFrom(abind,abind) importFrom(checkmate,assertClass) importFrom(checkmate,assertDataFrame) +importFrom(checkmate,checkInteger) +importFrom(checkmate,checkNumeric) importFrom(corpcor,pseudoinverse) importFrom(fmesher,fm_evaluator) importFrom(fmesher,fm_fem) diff --git a/NEWS.md b/NEWS.md index 7273433..b37d615 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# tinyVAST 0.6.0 + +* Change `integrate_output` interface by splitting `W_gz` and `V_gz` + into four vectors `area`, `type`, `covariate`, and `weighting_index` + to simplify documentations and improve naming + # tinyVAST 0.5.0 * Adding vignette showing how to fit multiple data types in an SDM diff --git a/R/fit.R b/R/fit.R index 297d306..bdf4ac9 100644 --- a/R/fit.R +++ b/R/fit.R @@ -96,7 +96,7 @@ #' model.offset model.response na.omit nlminb optimHess pnorm rnorm terms #' update.formula binomial poisson predict #' @importFrom TMB MakeADFun sdreport -#' @importFrom checkmate assertClass assertDataFrame +#' @importFrom checkmate assertClass assertDataFrame checkInteger checkNumeric #' @importFrom Matrix Cholesky solve #' @importFrom abind abind #' diff --git a/R/predict.R b/R/predict.R index 7288a92..2dae825 100644 --- a/R/predict.R +++ b/R/predict.R @@ -65,7 +65,9 @@ function( object, #' @title Monte Carlo integration for abundance #' -#' @description Applies Monte Carlo integration to approximate area-expanded abundance +#' @description +#' Calculates an estimator by summing across multiple predictions, +#' e.g., to approximate an integral when estimating area-expanded abundance. #' #' @inheritParams add_predictions #' @@ -77,31 +79,32 @@ function( object, #' then `integrate_output` is applying a midpoint approximation to the integral. #' @param area vector of values used for area-weighted expansion of #' estimated density surface for each row of `newdata` -#' which must have same lnegth as number of rows of `newdata`. -#' @param V_gz Integer-matrix providing settings for expansion, with two columns -#' and as many rows as `newdata`. If `V_gz` is missing, the first column is `1s`, -#' which indicates to expand each row by the first column of `W_gz`. If the -#' first column is `2`, then it weights the second -#' column of `W_gz` by the proportion at that row, e.g., to calculate -#' an abundance-weighted covariate. -#' If the first column is a `3`, then it weights that -#' row of `newdata` by a previous row of `newdata`, specified via the 2nd -#' column of `V_gz`. This 3rd option is used to weight a prediction for -#' one category based on predicted density of another category, e.g., -#' to calculate abundance-weighted condition in a bivariate model. -#' If the first column is a `0`, then that row of `newdata` is not included -#' when calculating the total across rows. Including a row of `newdata` with -#' first-column of `0` is useful, e.g., when calculating abundance at that -#' location, but where the eventual index uses abundance as weighting term -#' but without otherwise using the predicted density in calculating a total -#' value. -#' @param W_gz Matrix providing covariates for expansion, with two columns and -#' as many rows as `newdata`. If `W_gz` is missing, the first column is -#' populated from user-supplied values for `area`. -#' The second column can be used to provide a covariate -#' that are specified for use in expansion, e.g., to provide positional +#' with length of \code{nrow(newdata)}. +#' @param type Integer-vector indicating what type of expansion to apply to +#' each row of `newdata`, with length of \code{nrow(newdata)}. +#' \describe{ +#' \item{\code{type=1}}{Area-weighting: weight predictor by argument `area`} +#' \item{\code{type=2}}{Abundance-weighted covariate: weight `covariate` by +#' proportion of total in each row of `newdata`} +#' \item{\code{type=3}}{Abundance-weighted variable: weight predictor by +#' proportion of total in a prior row of `newdata`. +#' This option is used to weight a prediction for +#' one category based on predicted density of another category, e.g., +#' to calculate abundance-weighted condition in a bivariate model.} +#' \item{\code{type=0}}{Exclude from weighting: give weight of zero for +#' a given row of `newdata`. Including a row of `newdata` with +#' \code{type=0} is useful, e.g., when calculating abundance at that +#' location, where the eventual index uses abundance as weighting term +#' but without otherwise using the predicted density in calculating a total +#' value.} +#' } +#' @param covariate numeric-vector used to provide a covariate +#' that is used in expansion, e.g., to provide positional #' coordinates when calculating the abundance-weighted centroid with respect -#' to that coordinate. +#' to that coordinate. Only used for when \code{type=2}. +#' @param weighting_index integer-vector used to indicate a previous row +#' that is used to calculate a weighted average that is then +#' applied to the given row of `newdata`. Only used for when \code{type=3}. #' @param bias.correct logical indicating if bias correction should be applied using #' standard methods in [TMB::sdreport()] #' @param intern Do Laplace approximation on C++ side? Passed to [TMB::MakeADFun()]. @@ -110,7 +113,7 @@ function( object, #' `apply.epsilon` is sometimes substantially faster. #' #' @details -#' Analysts will often want to calculate some value by combing the predicted response at multiple +#' Analysts will often want to calculate some value by combiningg the predicted response at multiple #' locations, and potentially from multiple variables in a multivariate analysis. #' This arises in a univariate model, e.g., when calculating the integral under a predicted #' density function, which is approximated using a midpoint or Monte Carlo approximation @@ -121,36 +124,42 @@ function( object, #' approximation point (supplied by argument `area`), #' and summing across these area-expanded values. #' -#' In more complicated cases, an analyst can then use `V_gz` and -#' `W_gz` to calculate the weighted average +#' In more complicated cases, an analyst can then use `covariate` +#' to calculate the weighted average #' of a covariate (e.g., positional coordinates) for each Monte Carlo point. #' Alternatively, an analyst fitting a multivariate model might weight one variable -#' based on another, e.g., to calculate abundance-weighted average condition, or +#' based on another using `weighting_index`, e.g., +#' to calculate abundance-weighted average condition, or #' predator-expanded stomach contents. #' -#' In practice, supplying `V_gz` and `W_gz` requires two passes through the rows of +#' In practice, spatial integration in a multivariate model requires two passes through the rows of #' `newdata` when calculating a total value. In the following, we -#' use write equations using C++ indexing conventions such that indexing starts with 0, +#' write equations using C++ indexing conventions such that indexing starts with 0, #' to match the way that `integrate_output` expects indices to be supplied. -#' Given `mu_g` for each row of `newdata`, +#' Given inverse-link-transformed predictor \eqn{ \mu_g }, +#' function argument `type` as \eqn{ type_g } +#' function argument `area` as \eqn{ a_g }, +#' function argument `covariate` as \eqn{ x_g }, +#' function argument `weighting_index` as `\eqn{ h_g }` +#' function argument `weighting_index` as `\eqn{ h_g }` #' the first pass calculates: #' -#' \deqn{ \nu_g = \mu_g W_{g,0} } +#' \deqn{ \nu_g = \mu_g a_g } #' #' where the total value from this first pass is calculated as: #' #' \deqn{ \nu^* = \sum_{g=0}^{G-1} \nu_g } #' -#' The second pass then applies a further weighting, which depends upon \eqn{ V_{g,0} }, -#' and potentially upon \eqn{ V_{g,1} } and \eqn{ W_{g,1} }. +#' The second pass then applies a further weighting, which depends upon \eqn{ type_g }, +#' and potentially upon \eqn{ x_g } and \eqn{ h_g }. #' -#' If \eqn{V_{g,0} = 0} then \eqn{\phi_g = 0} +#' If \eqn{type_g = 0} then \eqn{\phi_g = 0} #' -#' If \eqn{V_{g,0} = 1} then \eqn{\phi_g = \nu_g} +#' If \eqn{type_g = 1} then \eqn{\phi_g = \nu_g} #' -#' If \eqn{V_{g,0} = 2} then \eqn{\phi_g = W_{g,1} \frac{\nu_g}{\nu^*} } +#' If \eqn{type_g = 2} then \eqn{\phi_g = x_g \frac{\nu_g}{\nu^*} } #' -#' If \eqn{V_{g,0} = 3} then \eqn{\phi_g = \frac{\nu_{V_{g,1}}}{\nu^*} \mu_g } +#' If \eqn{type_g = 3} then \eqn{\phi_g = \frac{\nu_{h_g}}{\nu^*} \mu_g } #' #' Finally, the total value from this second pass is calculated as: #' @@ -165,9 +174,10 @@ function( object, integrate_output <- function( object, newdata, - V_gz, area, - W_gz, + type = rep(1,nrow(newdata)), + weighting_index, + covariate, bias.correct = TRUE, apply.epsilon = FALSE, intern = FALSE @@ -180,25 +190,60 @@ function( object, newdata = newdata ) # , # remove_origdata = isFALSE(apply.epsilon) & isFALSE(bias.correct) ) - # Area-expanded sum - if(missing(W_gz)){ - # Default for area - if(missing(area)){ - area = rep(1, nrow(newdata)) - }else if(length(area)==1){ - area = rep(area, nrow(newdata)) - }else if( length(area)!=nrow(newdata) ){ - stop("Check length of `area`") - } - tmb_data2$W_gz = cbind(area, 0) - }else{ - tmb_data2$W_gz = W_gz + # Expansion area + if(missing(area)){ + area = rep(1, nrow(newdata)) + }else if(length(area)==1){ + area = rep(area, nrow(newdata)) + } + checkNumeric( area, lower=0, len=nrow(newdata), any.missing=FALSE ) + + # Expansion type + if(missing(type)){ + type = rep(1, nrow(newdata)) + }else if(length(type)==1){ + type = rep(type, nrow(type)) + } + checkInteger( type, lower=0, upper=3, len=nrow(newdata), any.missing=FALSE ) + + # Index for variable-weighted value + if(missing(weighting_index)){ + weighting_index = rep(0, nrow(newdata)) } - if(missing(V_gz)){ - tmb_data2$V_gz = cbind( rep(1,nrow(newdata)), 0 ) - }else{ - tmb_data2$V_gz = V_gz + if( any(weighting_index>=seq_len(nrow(newdata))) ){ + stop("Invalid `weighting_index`") + } + checkInteger( weighting_index, lower=0, len=nrow(newdata), any.missing=FALSE ) + + # + if(missing(covariate)){ + type = rep(0, nrow(newdata)) } + checkNumeric( covariate, len=nrow(newdata), any.missing=FALSE ) + + # Bundle + V_gz = cbind( type, weighting_index ) + W_gz = cbind( area, covariate ) + + # Area-expanded sum + #if(missing(W_gz)){ + # # Default for area + # if(missing(area)){ + # area = rep(1, nrow(newdata)) + # }else if(length(area)==1){ + # area = rep(area, nrow(newdata)) + # }else if( length(area)!=nrow(newdata) ){ + # stop("Check length of `area`") + # } + # tmb_data2$W_gz = cbind(area, 0) + #}else{ + # tmb_data2$W_gz = W_gz + #} + #if(missing(V_gz)){ + # tmb_data2$V_gz = cbind( rep(1,nrow(newdata)), 0 ) + #}else{ + # tmb_data2$V_gz = V_gz + #} # Abundance-weighted z #tmb_data2$W_gz = cbind( 1, newdata$x ) diff --git a/man/integrate_output.Rd b/man/integrate_output.Rd index 5365e16..c066801 100644 --- a/man/integrate_output.Rd +++ b/man/integrate_output.Rd @@ -7,9 +7,10 @@ integrate_output( object, newdata, - V_gz, area, - W_gz, + type = rep(1, nrow(newdata)), + weighting_index, + covariate, bias.correct = TRUE, apply.epsilon = FALSE, intern = FALSE @@ -25,35 +26,37 @@ If these locations are randomly drawn from a specified spatial domain, then total over that area. If locations are drawn sysmatically from a domain, then \code{integrate_output} is applying a midpoint approximation to the integral.} -\item{V_gz}{Integer-matrix providing settings for expansion, with two columns -and as many rows as \code{newdata}. If \code{V_gz} is missing, the first column is \verb{1s}, -which indicates to expand each row by the first column of \code{W_gz}. If the -first column is \code{2}, then it weights the second -column of \code{W_gz} by the proportion at that row, e.g., to calculate -an abundance-weighted covariate. -If the first column is a \code{3}, then it weights that -row of \code{newdata} by a previous row of \code{newdata}, specified via the 2nd -column of \code{V_gz}. This 3rd option is used to weight a prediction for +\item{area}{vector of values used for area-weighted expansion of +estimated density surface for each row of \code{newdata} +which must have same length as number of rows of \code{newdata}.} + +\item{type}{Integer-vector indicating what type of expansion to apply to +each row of \code{newdata}, with length of \code{nrow(newdata)}. +\describe{ +\item{\code{type=1}}{Area-weighting: weight predictor by argument \code{area}} +\item{\code{type=2}}{Abundance-weighted covariate: weight \code{covariate} by +proportion of total in each row of \code{newdata}} +\item{\code{type=3}}{Abundance-weighted variable: weight predictor by +proportion of total in a prior row of \code{newdata}. +This 3rd option is used to weight a prediction for one category based on predicted density of another category, e.g., -to calculate abundance-weighted condition in a bivariate model. -If the first column is a \code{0}, then that row of \code{newdata} is not included -when calculating the total across rows. Including a row of \code{newdata} with +to calculate abundance-weighted condition in a bivariate model.} +\item{\code{type=0}}{Exclude from weighting: give weight of zero for +a given row of \code{newdata}. Including a row of \code{newdata} with first-column of \code{0} is useful, e.g., when calculating abundance at that location, but where the eventual index uses abundance as weighting term but without otherwise using the predicted density in calculating a total value.} +}} -\item{area}{vector of values used for area-weighted expansion of -estimated density surface for each row of \code{newdata} -which must have same lnegth as number of rows of \code{newdata}.} +\item{weighting_index}{integer-vector used to indicate a previous row +that is used to calculate a weighted average that is then +applied to the given row of \code{newdata}. Only used for when \code{type=3}.} -\item{W_gz}{Matrix providing covariates for expansion, with two columns and -as many rows as \code{newdata}. If \code{W_gz} is missing, the first column is -populated from user-supplied values for \code{area}. -The second column can be used to provide a covariate +\item{covariate}{numeric-vector used to provide a covariate that are specified for use in expansion, e.g., to provide positional coordinates when calculating the abundance-weighted centroid with respect -to that coordinate.} +to that coordinate. Only used for when \code{type=2}.} \item{bias.correct}{logical indicating if bias correction should be applied using standard methods in \code{\link[TMB:sdreport]{TMB::sdreport()}}} @@ -68,7 +71,7 @@ rather than using the conventional method in \link[TMB:sdreport]{TMB::sdreport}? Applies Monte Carlo integration to approximate area-expanded abundance } \details{ -Analysts will often want to calculate some value by combing the predicted response at multiple +Analysts will often want to calculate some value by combiningg the predicted response at multiple locations, and potentially from multiple variables in a multivariate analysis. This arises in a univariate model, e.g., when calculating the integral under a predicted density function, which is approximated using a midpoint or Monte Carlo approximation @@ -79,36 +82,42 @@ by multiplying \code{mu_g} by the area associated with each midpoint or Monte Ca approximation point (supplied by argument \code{area}), and summing across these area-expanded values. -In more complicated cases, an analyst can then use \code{V_gz} and -\code{W_gz} to calculate the weighted average +In more complicated cases, an analyst can then use \code{covariate} +to calculate the weighted average of a covariate (e.g., positional coordinates) for each Monte Carlo point. Alternatively, an analyst fitting a multivariate model might weight one variable -based on another, e.g., to calculate abundance-weighted average condition, or +based on another using \code{weighting_index}, e.g., +to calculate abundance-weighted average condition, or predator-expanded stomach contents. -In practice, supplying \code{V_gz} and \code{W_gz} requires two passes through the rows of +In practice, spatial integration in a multivariate model requires two passes through the rows of \code{newdata} when calculating a total value. In the following, we -use write equations using C++ indexing conventions such that indexing starts with 0, +write equations using C++ indexing conventions such that indexing starts with 0, to match the way that \code{integrate_output} expects indices to be supplied. -Given \code{mu_g} for each row of \code{newdata}, +Given inverse-link-transformed predictor \eqn{ \mu_g }, +function argument \code{type} as \eqn{ type_g } +function argument \code{area} as \eqn{ a_g }, +function argument \code{covariate} as \eqn{ x_g }, +function argument \code{weighting_index} as \verb{\eqn{ h_g }} +function argument \code{weighting_index} as \verb{\eqn{ h_g }} the first pass calculates: -\deqn{ \nu_g = \mu_g W_{g,0} } +\deqn{ \nu_g = \mu_g a_g } where the total value from this first pass is calculated as: \deqn{ \nu^* = \sum_{g=0}^{G-1} \nu_g } -The second pass then applies a further weighting, which depends upon \eqn{ V_{g,0} }, -and potentially upon \eqn{ V_{g,1} } and \eqn{ W_{g,1} }. +The second pass then applies a further weighting, which depends upon \eqn{ type_g }, +and potentially upon \eqn{ x_g } and \eqn{ h_g }. -If \eqn{V_{g,0} = 0} then \eqn{\phi_g = 0} +If \eqn{type_g = 0} then \eqn{\phi_g = 0} -If \eqn{V_{g,0} = 1} then \eqn{\phi_g = \nu_g} +If \eqn{type_g = 1} then \eqn{\phi_g = \nu_g} -If \eqn{V_{g,0} = 2} then \eqn{\phi_g = W_{g,1} \frac{\nu_g}{\nu^*} } +If \eqn{type_g = 2} then \eqn{\phi_g = x_g \frac{\nu_g}{\nu^*} } -If \eqn{V_{g,0} = 3} then \eqn{\phi_g = \frac{\nu_{V_{g,1}}}{\nu^*} \mu_g } +If \eqn{type_g = 3} then \eqn{\phi_g = \frac{\nu_{h_g}}{\nu^*} \mu_g } Finally, the total value from this second pass is calculated as: diff --git a/tests/testthat/test-deviance-residuals.R b/tests/testthat/test-deviance-residuals.R index 76a880f..09e7de8 100644 --- a/tests/testthat/test-deviance-residuals.R +++ b/tests/testthat/test-deviance-residuals.R @@ -108,7 +108,7 @@ test_that("deviance residuals for poisson works", { # resid1^2 # mu = mytiny$rep$mu_i; # 2*y*log(y/mu) - (y-mu) - sign(y - mu) * pow(2*(y*log((1e-10 + y)/mu) - (y-mu)), 0.5) + # sign(y - mu) * pow(2*(y*log((1e-10 + y)/mu) - (y-mu)), 0.5) # resid1 # pow = function(a,b) a^b # Type = function(z)z From d26b2fd68f7604e518d69b1daef56bbcd55e7e28 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Mon, 29 Apr 2024 17:35:42 -0700 Subject: [PATCH 05/19] bug-fixes to previous commit --- R/predict.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/predict.R b/R/predict.R index 2dae825..fadf4f4 100644 --- a/R/predict.R +++ b/R/predict.R @@ -217,13 +217,13 @@ function( object, # if(missing(covariate)){ - type = rep(0, nrow(newdata)) + covariate = rep(0, nrow(newdata)) } checkNumeric( covariate, len=nrow(newdata), any.missing=FALSE ) # Bundle - V_gz = cbind( type, weighting_index ) - W_gz = cbind( area, covariate ) + tmb_data2$V_gz = cbind( type, weighting_index ) + tmb_data2$W_gz = cbind( area, covariate ) # Area-expanded sum #if(missing(W_gz)){ From 6620e52b9b40508ab340bda4916dd2e9e556dc80 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Wed, 1 May 2024 11:51:23 -0700 Subject: [PATCH 06/19] add gc(.) to integrate_output --- R/predict.R | 22 ++++++++++++++++++---- man/integrate_output.Rd | 24 ++++++++++++++++-------- 2 files changed, 34 insertions(+), 12 deletions(-) diff --git a/R/predict.R b/R/predict.R index fadf4f4..25ad320 100644 --- a/R/predict.R +++ b/R/predict.R @@ -66,7 +66,7 @@ function( object, #' @title Monte Carlo integration for abundance #' #' @description -#' Calculates an estimator by summing across multiple predictions, +#' Calculates an estimator for a derived quantity by summing across multiple predictions, #' e.g., to approximate an integral when estimating area-expanded abundance. #' #' @inheritParams add_predictions @@ -109,8 +109,8 @@ function( object, #' standard methods in [TMB::sdreport()] #' @param intern Do Laplace approximation on C++ side? Passed to [TMB::MakeADFun()]. #' @param apply.epsilon Apply epsilon bias correction using a manual calculation -#' rather than using the conventional method in [TMB::sdreport]? Using -#' `apply.epsilon` is sometimes substantially faster. +#' rather than using the conventional method in [TMB::sdreport]? +#' See details for more information. #' #' @details #' Analysts will often want to calculate some value by combiningg the predicted response at multiple @@ -168,7 +168,15 @@ function( object, #' and \eqn{\phi^*} is outputted by `integrate_output`, #' along with a standard error and potentially using #' the epsilon bias-correction estimator to correct for skewness and retransformation -#' bias. +#' bias. +#' +#' Standard bias-correction using \code{bias.correct=TRUE} can be slow, and in +#' some cases it might be faster to do \code{apply.epsilon=TRUE} and \code{intern=TRUE}. +#' However, that option is somewhat experimental, and a user might want to confirm +#' that the two options give identical results. Similarly, using \code{bias.correct=TRUE} +#' will still calculate the standard-error, whereas using +#' \code{apply.epsilon=TRUE} and \code{intern=TRUE} will not. +#' #' #' @export integrate_output <- @@ -284,6 +292,12 @@ function( object, rep = newobj$report() out = c( "Estimate"=rep$Metric, "Std. Error"=NA, "Est. (bias.correct)"=NA, "Std. (bias.correct)"=NA ) } + + # deal with memory internally + remove("newobj") + gc() + + # return stuff return(out) } diff --git a/man/integrate_output.Rd b/man/integrate_output.Rd index c066801..14b8e23 100644 --- a/man/integrate_output.Rd +++ b/man/integrate_output.Rd @@ -28,7 +28,7 @@ then \code{integrate_output} is applying a midpoint approximation to the integra \item{area}{vector of values used for area-weighted expansion of estimated density surface for each row of \code{newdata} -which must have same length as number of rows of \code{newdata}.} +with length of \code{nrow(newdata)}.} \item{type}{Integer-vector indicating what type of expansion to apply to each row of \code{newdata}, with length of \code{nrow(newdata)}. @@ -38,13 +38,13 @@ each row of \code{newdata}, with length of \code{nrow(newdata)}. proportion of total in each row of \code{newdata}} \item{\code{type=3}}{Abundance-weighted variable: weight predictor by proportion of total in a prior row of \code{newdata}. -This 3rd option is used to weight a prediction for +This option is used to weight a prediction for one category based on predicted density of another category, e.g., to calculate abundance-weighted condition in a bivariate model.} \item{\code{type=0}}{Exclude from weighting: give weight of zero for a given row of \code{newdata}. Including a row of \code{newdata} with -first-column of \code{0} is useful, e.g., when calculating abundance at that -location, but where the eventual index uses abundance as weighting term +\code{type=0} is useful, e.g., when calculating abundance at that +location, where the eventual index uses abundance as weighting term but without otherwise using the predicted density in calculating a total value.} }} @@ -54,7 +54,7 @@ that is used to calculate a weighted average that is then applied to the given row of \code{newdata}. Only used for when \code{type=3}.} \item{covariate}{numeric-vector used to provide a covariate -that are specified for use in expansion, e.g., to provide positional +that is used in expansion, e.g., to provide positional coordinates when calculating the abundance-weighted centroid with respect to that coordinate. Only used for when \code{type=2}.} @@ -62,13 +62,14 @@ to that coordinate. Only used for when \code{type=2}.} standard methods in \code{\link[TMB:sdreport]{TMB::sdreport()}}} \item{apply.epsilon}{Apply epsilon bias correction using a manual calculation -rather than using the conventional method in \link[TMB:sdreport]{TMB::sdreport}? Using -\code{apply.epsilon} is sometimes substantially faster.} +rather than using the conventional method in \link[TMB:sdreport]{TMB::sdreport}? +See details for more information.} \item{intern}{Do Laplace approximation on C++ side? Passed to \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}.} } \description{ -Applies Monte Carlo integration to approximate area-expanded abundance +Calculates an estimator for a derived quantity by summing across multiple predictions, +e.g., to approximate an integral when estimating area-expanded abundance. } \details{ Analysts will often want to calculate some value by combiningg the predicted response at multiple @@ -127,4 +128,11 @@ and \eqn{\phi^*} is outputted by \code{integrate_output}, along with a standard error and potentially using the epsilon bias-correction estimator to correct for skewness and retransformation bias. + +Standard bias-correction using \code{bias.correct=TRUE} can be slow, and in +some cases it might be faster to do \code{apply.epsilon=TRUE} and \code{intern=TRUE}. +However, that option is somewhat experimental, and a user might want to confirm +that the two options give identical results. Similarly, using \code{bias.correct=TRUE} +will still calculate the standard-error, whereas using +\code{apply.epsilon=TRUE} and \code{intern=TRUE} will not. } From ee6d2f9b316b3bdca7031271751e86c88e0b2b55 Mon Sep 17 00:00:00 2001 From: Sean Anderson Date: Thu, 2 May 2024 10:53:19 -0700 Subject: [PATCH 07/19] Fix 2 typos in docs --- R/fit.R | 4 ++-- man/tinyVAST.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/fit.R b/R/fit.R index bdf4ac9..caca3ea 100644 --- a/R/fit.R +++ b/R/fit.R @@ -57,7 +57,7 @@ #' from `data$variables`. #' @param delta_options a named list with slots for \code{delta_formula}, #' \code{delta_sem}, and \code{delta_dsem}. These follow the same format as -#' \code{family}, \code{sem}, and \code{dsem}, but specify options for the +#' \code{formula}, \code{sem}, and \code{dsem}, but specify options for the #' second linear predictor of a delta model, and are only used (or estimable) #' when a \code{\link[tinyVAST:families]{delta family}} is used for some samples. #' @param control Output from [tinyVASTcontrol()], used to define user @@ -75,7 +75,7 @@ #' #' the default `dsem=NULL` turns off all multivariate and temporal indexing, such #' that `spatial_graph` is then ignored, and the model collapses -#' to a standard model using \code{\link[mgcv]{gam}}. To specify a univeriate spatial model, +#' to a standard model using \code{\link[mgcv]{gam}}. To specify a univariate spatial model, #' the user must specify both `spatial_graph` and `dsem=""`, where the latter #' is then parsed to include a single exogenous variance for the single variable #' diff --git a/man/tinyVAST.Rd b/man/tinyVAST.Rd index fd34e7f..7f5d1f1 100644 --- a/man/tinyVAST.Rd +++ b/man/tinyVAST.Rd @@ -89,7 +89,7 @@ from \code{data$variables}.} \item{delta_options}{a named list with slots for \code{delta_formula}, \code{delta_sem}, and \code{delta_dsem}. These follow the same format as -\code{family}, \code{sem}, and \code{dsem}, but specify options for the +\code{formula}, \code{sem}, and \code{dsem}, but specify options for the second linear predictor of a delta model, and are only used (or estimable) when a \code{\link[tinyVAST:families]{delta family}} is used for some samples.} @@ -115,7 +115,7 @@ the space-variable interaction. the default \code{dsem=NULL} turns off all multivariate and temporal indexing, such that \code{spatial_graph} is then ignored, and the model collapses -to a standard model using \code{\link[mgcv]{gam}}. To specify a univeriate spatial model, +to a standard model using \code{\link[mgcv]{gam}}. To specify a univariate spatial model, the user must specify both \code{spatial_graph} and \code{dsem=""}, where the latter is then parsed to include a single exogenous variance for the single variable\tabular{ll}{ \strong{Model type} \tab \strong{How to specify} \cr From 2f73f15298aa7129432304d88349cec91559baab Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Fri, 3 May 2024 20:45:19 -0700 Subject: [PATCH 08/19] update docs --- R/predict.R | 10 +++++++--- man/integrate_output.Rd | 10 +++++++--- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/R/predict.R b/R/predict.R index 25ad320..ec5e1aa 100644 --- a/R/predict.R +++ b/R/predict.R @@ -66,8 +66,10 @@ function( object, #' @title Monte Carlo integration for abundance #' #' @description -#' Calculates an estimator for a derived quantity by summing across multiple predictions, -#' e.g., to approximate an integral when estimating area-expanded abundance. +#' Calculates an estimator for a derived quantity by summing across multiple predictions. +#' This can be used to approximate an integral when estimating area-expanded abundance, +#' abundance-weighting a covariate to calculate distribution shifts, +#' and/or weighting one model variable by another. #' #' @inheritParams add_predictions #' @@ -126,7 +128,9 @@ function( object, #' #' In more complicated cases, an analyst can then use `covariate` #' to calculate the weighted average -#' of a covariate (e.g., positional coordinates) for each Monte Carlo point. +#' of a covariate for each Monte Carlo point. For example, if the covariate is +#' positional coordinates or depth/elevation, then \code{type=2} +#' measures shifts in the average habitat utilization with respect to that covariate. #' Alternatively, an analyst fitting a multivariate model might weight one variable #' based on another using `weighting_index`, e.g., #' to calculate abundance-weighted average condition, or diff --git a/man/integrate_output.Rd b/man/integrate_output.Rd index 14b8e23..832b8f6 100644 --- a/man/integrate_output.Rd +++ b/man/integrate_output.Rd @@ -68,8 +68,10 @@ See details for more information.} \item{intern}{Do Laplace approximation on C++ side? Passed to \code{\link[TMB:MakeADFun]{TMB::MakeADFun()}}.} } \description{ -Calculates an estimator for a derived quantity by summing across multiple predictions, -e.g., to approximate an integral when estimating area-expanded abundance. +Calculates an estimator for a derived quantity by summing across multiple predictions. +This can be used to approximate an integral when estimating area-expanded abundance, +abundance-weighting a covariate to calculate distribution shifts, +and/or weighting one model variable by another. } \details{ Analysts will often want to calculate some value by combiningg the predicted response at multiple @@ -85,7 +87,9 @@ and summing across these area-expanded values. In more complicated cases, an analyst can then use \code{covariate} to calculate the weighted average -of a covariate (e.g., positional coordinates) for each Monte Carlo point. +of a covariate for each Monte Carlo point. For example, if the covariate is +positional coordinates or depth/elevation, then \code{type=2} +measures shifts in the average habitat utilization with respect to that covariate. Alternatively, an analyst fitting a multivariate model might weight one variable based on another using \code{weighting_index}, e.g., to calculate abundance-weighted average condition, or From f06c0ec13b72e3f1e65389cd01277355b13d012b Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Tue, 14 May 2024 09:25:00 -0700 Subject: [PATCH 09/19] update docs somewhat --- NEWS.md | 2 ++ R/fit.R | 7 ++++--- man/tinyVAST.Rd | 7 ++++--- src/tinyVAST.cpp | 8 ++++++++ 4 files changed, 18 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index b37d615..362a857 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,6 +3,8 @@ * Change `integrate_output` interface by splitting `W_gz` and `V_gz` into four vectors `area`, `type`, `covariate`, and `weighting_index` to simplify documentations and improve naming +* Fix bug where cloglog and logit links were not previously implemented + for use in `predict` and `integrate_output` # tinyVAST 0.5.0 diff --git a/R/fit.R b/R/fit.R index caca3ea..c52d1a2 100644 --- a/R/fit.R +++ b/R/fit.R @@ -20,11 +20,12 @@ #' `dsem=NULL` disables the space-variable interaction; see #' [make_dsem_ram()] or [make_eof_ram()]. #' @param family A function returning a class \code{family}, including [gaussian()], -#' [lognormal()], or [tweedie()]. Alternatively, can be a named list of +#' [lognormal()], [tweedie()], [binomial()], [Gamma()], or [poisson()]. +#' Alternatively, can be a named list of #' these functions, with names that match levels of #' \code{data$distribution_column} to allow different -#' families by row of data. Delta model families are possible. -#' See \code{\link[tinyVAST:families]{Families}}, +#' families by row of data. Delta model families are possible, and see +#' \code{\link[tinyVAST:families]{Families}} for delta-model options, #' @param space_columns A string or character vector that indicates #' the column(s) of `data` indicating the location of each sample. #' When `spatial_graph` is an `igraph` object, `space_columns` is a string with diff --git a/man/tinyVAST.Rd b/man/tinyVAST.Rd index 7f5d1f1..9342024 100644 --- a/man/tinyVAST.Rd +++ b/man/tinyVAST.Rd @@ -44,11 +44,12 @@ constructing a space-variable interaction. \code{\link[=make_dsem_ram]{make_dsem_ram()}} or \code{\link[=make_eof_ram]{make_eof_ram()}}.} \item{family}{A function returning a class \code{family}, including \code{\link[=gaussian]{gaussian()}}, -\code{\link[=lognormal]{lognormal()}}, or \code{\link[=tweedie]{tweedie()}}. Alternatively, can be a named list of +\code{\link[=lognormal]{lognormal()}}, \code{\link[=tweedie]{tweedie()}}, \code{\link[=binomial]{binomial()}}, \code{\link[=Gamma]{Gamma()}}, or \code{\link[=poisson]{poisson()}}. +Alternatively, can be a named list of these functions, with names that match levels of \code{data$distribution_column} to allow different -families by row of data. Delta model families are possible. -See \code{\link[tinyVAST:families]{Families}},} +families by row of data. Delta model families are possible, and see +\code{\link[tinyVAST:families]{Families}} for delta-model options,} \item{space_columns}{A string or character vector that indicates the column(s) of \code{data} indicating the location of each sample. diff --git a/src/tinyVAST.cpp b/src/tinyVAST.cpp index 9c91529..8bddfc1 100644 --- a/src/tinyVAST.cpp +++ b/src/tinyVAST.cpp @@ -759,6 +759,14 @@ Type objective_function::operator() (){ case log_link: mu_g(g) = exp(p_g(g)); break; + case logit_link: + mu_g(g) = invlogit(p_g(g)); + break; + case cloglog_link: + mu_g(g) = Type(1.0) - exp( -1*exp(p_g(g)) ); + break; + default: + error("Link not implemented."); } } if( components_e(e_g(g))==2 ){ From 920ff34149ee2b68c0b22870d20626182673ec57 Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Wed, 15 May 2024 13:29:20 -0700 Subject: [PATCH 10/19] adding informative error message ... ... when user supplies new levels for factor --- R/predict.R | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/R/predict.R b/R/predict.R index ec5e1aa..27a2323 100644 --- a/R/predict.R +++ b/R/predict.R @@ -333,9 +333,14 @@ function( object, if( !(pred %in% colnames(newdata)) ){ stop("Missing ", pred, " in newdata") } + # If predictor is a factor, make sure newdata has same levels and no new levels if( is.factor(origdata[,pred]) ){ newdata[,pred] = factor(newdata[,pred], levels=levels(origdata[,pred])) } + # Check for NAs after factor(..., levels=levels(origdata)), which converts new levels to NAs + if( any(is.na(newdata[,pred])) ){ + stop("`newdata` column ", pred, " has NAs, whether because it's a factor with new levels or otherwise.") + } } make_covariates <- @@ -473,7 +478,7 @@ function( object, stop("Check output of `add_predictions` for variables with unexpected length") } if( any( sapply(tmb_data2[c('X_gj','Z_gk','X2_gj','Z2_gk')],FUN=nrow) != nrow(newdata) ) ){ - stop("Check output of `add_predictions` for NAs") + stop("Check output of `add_predictions` for mismatch in dimensions") } return( tmb_data2 ) From e6020a79234ab24a03a23c31b8a5855d818f3b5d Mon Sep 17 00:00:00 2001 From: Jim Thorson <50178738+James-Thorson-NOAA@users.noreply.github.com> Date: Wed, 15 May 2024 14:22:38 -0700 Subject: [PATCH 11/19] further improve error messages for add_predictions(.) --- R/predict.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/predict.R b/R/predict.R index 27a2323..ee7f3fa 100644 --- a/R/predict.R +++ b/R/predict.R @@ -333,13 +333,19 @@ function( object, if( !(pred %in% colnames(newdata)) ){ stop("Missing ", pred, " in newdata") } - # If predictor is a factor, make sure newdata has same levels and no new levels + # + if( is.factor(origdata[,pred]) ){ + if( !all(unique(newdata[,pred]) %in% levels(origdata[,pred])) ){ + stop("`newdata` column ", pred, "has new levels not present in the fitted data, which is not permitted") + } + } + # If predictor is a factor, make sure newdata has same level order if( is.factor(origdata[,pred]) ){ newdata[,pred] = factor(newdata[,pred], levels=levels(origdata[,pred])) } # Check for NAs after factor(..., levels=levels(origdata)), which converts new levels to NAs if( any(is.na(newdata[,pred])) ){ - stop("`newdata` column ", pred, " has NAs, whether because it's a factor with new levels or otherwise.") + stop("`newdata` column ", pred, " has NAs, which is not permitted") } } From b33f767d05ab4f5faaaddb66f0a83b3a63b88331 Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 09:25:32 -0700 Subject: [PATCH 12/19] add warning about tibbles --- R/fit.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/R/fit.R b/R/fit.R index c52d1a2..c737da4 100644 --- a/R/fit.R +++ b/R/fit.R @@ -157,6 +157,7 @@ function( formula, #if( !is.data.frame(data) ) stop("`data` must be a data frame", call. = FALSE) assertClass(control, "tinyVASTcontrol") assertDataFrame(data) + if(inherits(data,"tbl")) stop("`data` must be a data.frame and cannot be a tibble") ############## # input telescoping @@ -206,6 +207,11 @@ function( formula, c_i = integer(0) } + # variables can't have commas, because it conflicts with how `sem` and `dsem` are parsed + if( grepl(pattern=",", x=variables, fixed=TRUE) ){ + stop("`variables` cannot include any commas") + } + ############## # DSEM RAM constructor ############## @@ -749,6 +755,7 @@ function( nlminb_loops = 1, newton_loops = 0, eval.max = 1000, iter.max = 1000, + steps = 10, getsd = TRUE, silent = getOption("tinyVAST.silent", TRUE), trace = getOption("tinyVAST.trace", 0), @@ -768,6 +775,7 @@ function( nlminb_loops = 1, newton_loops = newton_loops, eval.max = eval.max, iter.max = iter.max, + steps = steps getsd = getsd, silent = silent, trace = trace, From d7d635fd839aa8b4fbdc491a16da840de1557c4f Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 09:26:00 -0700 Subject: [PATCH 13/19] Update DESCRIPTION --- DESCRIPTION | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 66183d6..6708f58 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,8 +2,8 @@ Type: Package Package: tinyVAST Title: Vector Autoregressive Spatio-Temporal Model Using Minimal Feature-Set -Version: 0.6.0 -Date: 2024-04-29 +Version: 0.6.1 +Date: 2024-05-21 Authors@R: c( person(c("James", "T."), "Thorson", , "James.Thorson@noaa.gov", role = c("aut", "cre"), From b580166d581cfa762ed01a5f94de4c854b46a98f Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 09:27:37 -0700 Subject: [PATCH 14/19] Update DESCRIPTION --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 6708f58..8cde948 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Type: Package Package: tinyVAST Title: Vector Autoregressive Spatio-Temporal Model Using Minimal Feature-Set -Version: 0.6.1 +Version: 0.6.0 Date: 2024-05-21 Authors@R: c( person(c("James", "T."), "Thorson", , "James.Thorson@noaa.gov", From 4284cf069665b7664cafe5d33aff72d6a5463e60 Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 09:45:08 -0700 Subject: [PATCH 15/19] fix bug (typo from other project!) --- R/fit.R | 2 -- 1 file changed, 2 deletions(-) diff --git a/R/fit.R b/R/fit.R index c737da4..58fd8ad 100644 --- a/R/fit.R +++ b/R/fit.R @@ -755,7 +755,6 @@ function( nlminb_loops = 1, newton_loops = 0, eval.max = 1000, iter.max = 1000, - steps = 10, getsd = TRUE, silent = getOption("tinyVAST.silent", TRUE), trace = getOption("tinyVAST.trace", 0), @@ -775,7 +774,6 @@ function( nlminb_loops = 1, newton_loops = newton_loops, eval.max = eval.max, iter.max = iter.max, - steps = steps getsd = getsd, silent = silent, trace = trace, From da01da13f80ab14591e64a631dcd98dd07ae0eb2 Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 10:04:27 -0700 Subject: [PATCH 16/19] set calculate_deviance_explained=FALSE by default --- R/fit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/fit.R b/R/fit.R index 58fd8ad..3e66977 100644 --- a/R/fit.R +++ b/R/fit.R @@ -764,7 +764,7 @@ function( nlminb_loops = 1, gmrf_parameterization = c("separable","projection"), estimate_delta0 = FALSE, getJointPrecision = FALSE, - calculate_deviance_explained = TRUE ){ + calculate_deviance_explained = FALSE ){ gmrf_parameterization = match.arg(gmrf_parameterization) From 033a675ed599d92cd4897eca70585e014c9ede8c Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 10:10:14 -0700 Subject: [PATCH 17/19] Update fit.R --- R/fit.R | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/R/fit.R b/R/fit.R index 3e66977..c4cfd88 100644 --- a/R/fit.R +++ b/R/fit.R @@ -203,15 +203,14 @@ function( formula, } if( length(variables) > 0 ){ c_i = match( data[,variable_column], variables ) + # variables can't have commas, because it conflicts with how `sem` and `dsem` are parsed + if( grepl(pattern=",", x=variables, fixed=TRUE) ){ + stop("`variables` cannot include any commas") + } }else{ c_i = integer(0) } - # variables can't have commas, because it conflicts with how `sem` and `dsem` are parsed - if( grepl(pattern=",", x=variables, fixed=TRUE) ){ - stop("`variables` cannot include any commas") - } - ############## # DSEM RAM constructor ############## From ed9a378562973db0167a7b89461a4da41e24ff7e Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 10:31:28 -0700 Subject: [PATCH 18/19] fix bug in variable-names check --- R/fit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/fit.R b/R/fit.R index c4cfd88..759e2b5 100644 --- a/R/fit.R +++ b/R/fit.R @@ -204,7 +204,7 @@ function( formula, if( length(variables) > 0 ){ c_i = match( data[,variable_column], variables ) # variables can't have commas, because it conflicts with how `sem` and `dsem` are parsed - if( grepl(pattern=",", x=variables, fixed=TRUE) ){ + if( any(grepl(pattern=",", x=variables, fixed=TRUE)) ){ stop("`variables` cannot include any commas") } }else{ From bb64a8a32af80a1a10f0a6885d7be8fa012d06f8 Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Tue, 21 May 2024 10:47:50 -0700 Subject: [PATCH 19/19] change back to `calculate_deviance_explained = TRUE` --- R/fit.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/fit.R b/R/fit.R index 759e2b5..2a2207e 100644 --- a/R/fit.R +++ b/R/fit.R @@ -763,7 +763,7 @@ function( nlminb_loops = 1, gmrf_parameterization = c("separable","projection"), estimate_delta0 = FALSE, getJointPrecision = FALSE, - calculate_deviance_explained = FALSE ){ + calculate_deviance_explained = TRUE ){ gmrf_parameterization = match.arg(gmrf_parameterization)