From 6fb2694e916a9fbaf4f884e587fcf4d442742bca Mon Sep 17 00:00:00 2001 From: "jesus.pena@bsc.es" Date: Thu, 31 Jan 2019 12:17:03 +0100 Subject: [PATCH 01/16] First versio AnoMultiOptions function --- R/AnoMultiOptions.R | 140 +++++++++++++++++++++++++++++++++++++++++ man/AnoMultiOptions.Rd | 46 ++++++++++++++ 2 files changed, 186 insertions(+) create mode 100644 R/AnoMultiOptions.R create mode 100644 man/AnoMultiOptions.Rd diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R new file mode 100644 index 00000000..054ffaf4 --- /dev/null +++ b/R/AnoMultiOptions.R @@ -0,0 +1,140 @@ +#'Anomalies relative to a climatology along sdates with or without cross-validation +#' +#'@description This function computes the anomalies relative to a climatology computed along the starting dates +#'dimension allowing the application or not of crossvalidated climatologies. The computation is carried out +#'independently for experimental and observational data products. +#' +#'Climatologies are computed using the Clim() function that allows the combination of different methods +#'(kharin and NDV flags: see documentation) +#' +#'@param var_exp a multidimensional array related to experimental simulations with the specific dimensions required by Clim(): +#' i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to +#' c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) +#' +#'@param var_obs a multidimensional array related to observations with the specific dimensions required by Clim(): +#' i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to +#' c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) +#' +#'@param cross a logical value indicating whether cross-validation should be applied or not. Default = FALSE +#' +#'@param memb a logical value indicating whether Clim() computes one climatology for each experimental data product member(TRUE) or +#'it computes one sole climatology for all members (FALSE). Default = TRUE +#' +#'@param kharin a logical value indicating whether Clim() applies Kharin method or not. Default = FALSE +#' +#'@param NDV a logical value indicating whether Clim() applies Fukhar method or not. Default = FALSE +#' +#'@return A list of length 3: +#'\itemize{ +#'\item\code{$dtr.ref}{An array with the same dimensions as the input \code{data}, but with the time dimension reduce from daily to monthly or seasonal resolution depending on the selected resolution in \code{by.season}.} +#'\item\code{$year}{A vector of the corresponding years.} +#'\item\code{$season}{A vector of the seasons or months corresponding to the resolution selected in \code{by.season}}} +#' +#' +#'@examples +#'## Example 1: +#'t <- 1 : 360 +#'obsS <- sin(t * pi / 180) +#'expS <- obsS + rnorm(360, 0., 0.1) +#'res <- AnoMultiOptions() +#'str(res) +#' +#'@export +AnoMultiOptions <- function(var_exp, var_obs, cross = FALSE, memb = TRUE, kharin = FALSE, NDV = FALSE) { + # Duplicate clim dimensions to have same dimensions as var + + + dimexp <- dim(var_exp) + dimobs <- dim(var_obs) + + #Input variables dimension checks + if (length(dimexp) < 4 | length(dimobs) < 4) { + stop("Check var_exp or var_obs dimensions, at least 4 dim needed : c(nexp/nobs, nmemb, nsdates, nltime)") + } + #Dimensions lengths for nsdates, nltime and optional nlevel,nlon and nlat should be match for exp and obs + for (jn in 3:max(length(dimexp), length(dimobs))) { + if (dimexp[jn] != dimobs[jn]) { + stop("Wrong input dimensions: Dimensions lengths for nsdates, nltime (and optionally nlevel,nlon and nlat) should match for exp and obs") + } + } + + #Flags checks + if (!is.logical(cross) | !is.logical(memb) | !is.logical(kharin) | !is.logical(NDV)) { + stop("Flags 'cross', 'memb', 'kharin' and 'NDV' should be all logical.") + } + if (length(cross) > 1) { + cross <- cross[1] + warning("Parameter ") + } + + + # + # Enlarge the number of dimensions of var_exp and var_obs to 7 if necessary + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # + + var_exp <- Enlarge(var_exp, 7) + var_obs <- Enlarge(var_obs, 7) + + enl_ano_exp <- array(dim = dim(var_exp)) + enl_ano_obs <- array(dim = dim(var_obs)) + + + # Computation of anomalies with cross-validation + if (cross) { + + for (jsdat in 1:dimexp[3]) { + tmp1 <- array(var_exp[, , jsdat, , , , ], dim = dim(var_exp)[-3]) + tmp2 <- array(var_obs[, , jsdat, , , , ], dim = dim(var_obs)[-3]) + tmp3 <- array(var_exp[, , -jsdat, , , , ], dim = c(dimexp[1:2], + dimexp[3] - 1, + dim(var_exp)[4:7])) + tmp4 <- array(var_obs[, , -jsdat, , , , ], dim = c(dimobs[1:2], + dimobs[3] - 1, + dim(var_obs)[4:7])) + tmp <- Clim(tmp3, tmp4, memb) + if (memb) { + clim_exp <- tmp$clim_exp + clim_obs <- tmp$clim_obs + } else { + clim_exp <- InsertDim(tmp$clim_exp, 2, dimexp[2]) + clim_obs <- InsertDim(tmp$clim_obs, 2, dimobs[2]) + } + enl_ano_exp[, , jsdat, , , , ] <- tmp1 - clim_exp + enl_ano_obs[, , jsdat, , , , ] <- tmp2 - clim_obs + } + + + # Computation of anomalies without cross-validation + } else { + tmp=Clim(var_exp,var_obs, memb,kharin, NDV) + + if (memb) { + clim_exp <- tmp$clim_exp + clim_obs <- tmp$clim_obs + } else { + clim_exp <- InsertDim(tmp$clim_exp, 2, dimexp[2]) + clim_obs <- InsertDim(tmp$clim_obs, 2, dimobs[2]) + } + + clim_exp <- InsertDim(clim_exp, 3, dimexp[3]) + clim_obs <- InsertDim(clim_obs, 3, dimobs[3]) + + enl_ano_exp <- var_exp - clim_exp + enl_ano_obs <- var_obs - clim_obs + + } + + + # Reduce the number of dimensions to the original one + ano_exp <- array(dim = dimexp) + ano_obs <- array(dim = dimobs) + ano_exp[] <- enl_ano_exp + ano_obs[] <- enl_ano_obs + + + + # Outputs + # ~~~~~~~~~ + invisible(list(ano_exp = ano_exp, ano_obs = ano_obs)) +} diff --git a/man/AnoMultiOptions.Rd b/man/AnoMultiOptions.Rd new file mode 100644 index 00000000..22a80c40 --- /dev/null +++ b/man/AnoMultiOptions.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/AnoMultiOptions.R +\name{AnoMultiOptions} +\alias{AnoMultiOptions} +\title{Anomalies relative to a climatology along sdates with or without cross-validation} +\usage{ +AnoMultiOptions(var_exp, var_obs, cross = FALSE, memb = TRUE, + kharin = FALSE, NDV = FALSE) +} +\arguments{ +\item{var_exp}{a multidimensional array related to experimental simulations with the specific dimensions required by Clim(): +i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to +c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} + +\item{var_obs}{a multidimensional array related to observations with the specific dimensions required by Clim(): +i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to +c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} + +\item{cross}{a logical value indicating whether cross-validation should be applied or not. Default = FALSE} + +\item{memb}{a logical value indicating whether Clim() computes one climatology for each experimental data product member(TRUE) or +it computes one sole climatology for all members (FALSE). Default = TRUE} + +\item{kharin}{a logical value indicating whether Clim() applies Kharin method or not. Default = FALSE} + +\item{NDV}{a logical value indicating whether Clim() applies Fukhar method or not. Default = FALSE} +} +\value{ +an array +} +\description{ +This function computes the anomalies relative to a climatology computed along the starting dates +dimension allowing the application or not of crossvalidated climatologies. The computation is carried out +independently for experimental and observational data products. + +Climatologies are computed using the Clim() function that allows the combination of different methods +(kharin and NDV flags: see documentation) +} +\examples{ +## Example 1: +t <- 1 : 360 +obsS <- sin(t * pi / 180) +expS <- obsS + rnorm(360, 0., 0.1) + +} + -- GitLab From b16cb2b5f725460c505127c216e3cf70f6341385 Mon Sep 17 00:00:00 2001 From: "jesus.pena@bsc.es" Date: Fri, 1 Feb 2019 18:08:58 +0100 Subject: [PATCH 02/16] Improved script formatting of AnoMultiOptions function --- R/AnoMultiOptions.R | 184 ++++++++++++++++++++--------------------- man/AnoMultiOptions.Rd | 74 +++++++++++------ 2 files changed, 141 insertions(+), 117 deletions(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 054ffaf4..3192e517 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -1,113 +1,119 @@ #'Anomalies relative to a climatology along sdates with or without cross-validation #' -#'@description This function computes the anomalies relative to a climatology computed along the starting dates -#'dimension allowing the application or not of crossvalidated climatologies. The computation is carried out -#'independently for experimental and observational data products. +#'@description This function computes the anomalies relative to a climatology computed along the +#'starting dates dimension allowing the application or not of crossvalidated climatologies. The +#'computation is carried out independently for experimental and observational data products. #' -#'Climatologies are computed using the Clim() function that allows the combination of different methods -#'(kharin and NDV flags: see documentation) +#'@param var_exp A multidimensional array related to experimental simulations with the specific +#'dimensions required by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) +#'up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) #' -#'@param var_exp a multidimensional array related to experimental simulations with the specific dimensions required by Clim(): -#' i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to +#'@param var_obs A multidimensional array related to observations with the specific dimensions required +#'by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to #' c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) #' -#'@param var_obs a multidimensional array related to observations with the specific dimensions required by Clim(): -#' i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to -#' c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) -#' -#'@param cross a logical value indicating whether cross-validation should be applied or not. Default = FALSE +#'@param cross A logical value indicating whether cross-validation should be applied or not. Default = FALSE. #' -#'@param memb a logical value indicating whether Clim() computes one climatology for each experimental data product member(TRUE) or -#'it computes one sole climatology for all members (FALSE). Default = TRUE +#'@param memb A logical value indicating whether Clim() computes one climatology for each experimental data +#'product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE. #' -#'@param kharin a logical value indicating whether Clim() applies Kharin method or not. Default = FALSE #' -#'@param NDV a logical value indicating whether Clim() applies Fukhar method or not. Default = FALSE -#' -#'@return A list of length 3: +#'@return A list of length 2: #'\itemize{ -#'\item\code{$dtr.ref}{An array with the same dimensions as the input \code{data}, but with the time dimension reduce from daily to monthly or seasonal resolution depending on the selected resolution in \code{by.season}.} -#'\item\code{$year}{A vector of the corresponding years.} -#'\item\code{$season}{A vector of the seasons or months corresponding to the resolution selected in \code{by.season}}} +#'\item\code{$var_exp}{ Array of anomalies for \code{var_exp} with same dimensions input.} +#'\item\code{$var_obs}{ Array of anomalies for \code{var_obs} with same dimensions input.}} #' #' #'@examples -#'## Example 1: -#'t <- 1 : 360 -#'obsS <- sin(t * pi / 180) -#'expS <- obsS + rnorm(360, 0., 0.1) -#'res <- AnoMultiOptions() -#'str(res) +#' +#'# Example 1: +#'Two arrays to compute anomalies: \code{var_exp} and \code{var_obs} +#' +#'dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100) +#'dim_obs <- c(mod = 1, member = 1, sdates = 2, ftime = 100) +#' +#'var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) +#'var_obs <- array(rnorm(prod(dim_obs)), dim = dim_obs) +#' +#'Computation of anomalies without cross-validation +#'anoM <- AnoMultiOptions(var_exp = var_exp, var_obs = var_obs, +#' cross = FALSE, memb = TRUE) +#'ano_exp <- anoM$var_exp +#'ano_obs <- anoM$var_obs +#' +#' +#'# Example 2: +#'One array to compute anomalie: \code{var_exp} +#' +#'dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100, nlon = 100, nlat = 50) +#' +#'var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) +#' +#'Computation of anomalies with cross-validation using kharin and NDV methods +#'anoM <- AnoMultiOptions(var_exp , +#' cross = TRUE, memb = TRUE) +#'ano_exp <- anoM$var_exp +#' +#'@seealso +#'\code{\link{Ano_CrossValid, }} +#'\code{\link{Clim.}} +#' #' #'@export -AnoMultiOptions <- function(var_exp, var_obs, cross = FALSE, memb = TRUE, kharin = FALSE, NDV = FALSE) { - # Duplicate clim dimensions to have same dimensions as var - + +AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb = TRUE) { dimexp <- dim(var_exp) dimobs <- dim(var_obs) - #Input variables dimension checks + if (is.null(var_exp) & is.null(var_obs)) { + stop("Not-NULL array var_exp or var_obs or both should be provided.") + } + if (is.null(var_exp)) { + var_exp <- var_obs + dimexp <- dim(var_exp) + warning('No experimental data (var_exp) provided.') + } + if (is.null(var_obs)) { + var_obs <- var_exp + dimobs <- dim(var_obs) + warning('No observational data provided (var_obs).') + } if (length(dimexp) < 4 | length(dimobs) < 4) { - stop("Check var_exp or var_obs dimensions, at least 4 dim needed : c(nexp/nobs, nmemb, nsdates, nltime)") + stop("Invalid input dimensions: At least 4 dimensions needed for var_exp or var_obs, + c(nexp/nobs, nmemb, nsdates, nltime).") + } + if (length(dimexp) != length(dimobs)) { + stop('Invalid input dimensions: var_exp and var_obs should have same number of dimensions.') } - #Dimensions lengths for nsdates, nltime and optional nlevel,nlon and nlat should be match for exp and obs for (jn in 3:max(length(dimexp), length(dimobs))) { - if (dimexp[jn] != dimobs[jn]) { - stop("Wrong input dimensions: Dimensions lengths for nsdates, nltime (and optionally nlevel,nlon and nlat) should match for exp and obs") + if (dimexp[jn] != dimobs[jn]) { + stop("Invalid input dimensions: var_exp and var_obs dimensions should have same length + for nsdates and nltime (and for nlevel,nlon and nlat if present).") } } - - #Flags checks - if (!is.logical(cross) | !is.logical(memb) | !is.logical(kharin) | !is.logical(NDV)) { - stop("Flags 'cross', 'memb', 'kharin' and 'NDV' should be all logical.") + if (dimexp[3] == 1 | dimobs[3] == 1) { + stop("Invalid input dimensions: sdate dimension length should be greater than 1") } - if (length(cross) > 1) { - cross <- cross[1] - warning("Parameter ") + if (!is.logical(cross) | !is.logical(memb) ) { + stop("Invalid flags used: Flags 'cross' and 'memb' should be all logical.") + } + if (length(cross) > 1 | length(memb) > 1 ) { + warning("Invalid flags used: Flags 'cross' and 'memb' should all have length 1.") } - # - # Enlarge the number of dimensions of var_exp and var_obs to 7 if necessary - # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - # - - var_exp <- Enlarge(var_exp, 7) - var_obs <- Enlarge(var_obs, 7) - - enl_ano_exp <- array(dim = dim(var_exp)) - enl_ano_obs <- array(dim = dim(var_obs)) - + # Computating anomalies + #--------------------- - # Computation of anomalies with cross-validation + # With cross-validation if (cross) { + ano <- Ano_CrossValid(var_exp, var_obs, memb = memb) + - for (jsdat in 1:dimexp[3]) { - tmp1 <- array(var_exp[, , jsdat, , , , ], dim = dim(var_exp)[-3]) - tmp2 <- array(var_obs[, , jsdat, , , , ], dim = dim(var_obs)[-3]) - tmp3 <- array(var_exp[, , -jsdat, , , , ], dim = c(dimexp[1:2], - dimexp[3] - 1, - dim(var_exp)[4:7])) - tmp4 <- array(var_obs[, , -jsdat, , , , ], dim = c(dimobs[1:2], - dimobs[3] - 1, - dim(var_obs)[4:7])) - tmp <- Clim(tmp3, tmp4, memb) - if (memb) { - clim_exp <- tmp$clim_exp - clim_obs <- tmp$clim_obs - } else { - clim_exp <- InsertDim(tmp$clim_exp, 2, dimexp[2]) - clim_obs <- InsertDim(tmp$clim_obs, 2, dimobs[2]) - } - enl_ano_exp[, , jsdat, , , , ] <- tmp1 - clim_exp - enl_ano_obs[, , jsdat, , , , ] <- tmp2 - clim_obs - } - - - # Computation of anomalies without cross-validation + # Without cross-validation } else { - tmp=Clim(var_exp,var_obs, memb,kharin, NDV) + tmp <- Clim(var_exp , var_obs , memb ) if (memb) { clim_exp <- tmp$clim_exp @@ -120,21 +126,15 @@ AnoMultiOptions <- function(var_exp, var_obs, cross = FALSE, memb = TRUE, kharin clim_exp <- InsertDim(clim_exp, 3, dimexp[3]) clim_obs <- InsertDim(clim_obs, 3, dimobs[3]) - enl_ano_exp <- var_exp - clim_exp - enl_ano_obs <- var_obs - clim_obs - + ano_exp <- var_exp - clim_exp + ano_obs <- var_obs - clim_obs + + ano <- list(ano_exp = ano_exp, ano_obs = ano_obs) } - - - # Reduce the number of dimensions to the original one - ano_exp <- array(dim = dimexp) - ano_obs <- array(dim = dimobs) - ano_exp[] <- enl_ano_exp - ano_obs[] <- enl_ano_obs - - - + # Outputs # ~~~~~~~~~ - invisible(list(ano_exp = ano_exp, ano_obs = ano_obs)) + invisible(ano) } + + \ No newline at end of file diff --git a/man/AnoMultiOptions.Rd b/man/AnoMultiOptions.Rd index 22a80c40..4cd26605 100644 --- a/man/AnoMultiOptions.Rd +++ b/man/AnoMultiOptions.Rd @@ -4,43 +4,67 @@ \alias{AnoMultiOptions} \title{Anomalies relative to a climatology along sdates with or without cross-validation} \usage{ -AnoMultiOptions(var_exp, var_obs, cross = FALSE, memb = TRUE, - kharin = FALSE, NDV = FALSE) +AnoMultiOptions(var_exp = NULL, var_obs = NULL, cross = FALSE, + memb = TRUE) } \arguments{ -\item{var_exp}{a multidimensional array related to experimental simulations with the specific dimensions required by Clim(): -i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to -c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} +\item{var_exp}{A multidimensional array related to experimental simulations with the specific +dimensions required by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) +up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} -\item{var_obs}{a multidimensional array related to observations with the specific dimensions required by Clim(): -i.e. dim(var_exp) goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to +\item{var_obs}{A multidimensional array related to observations with the specific dimensions required +by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} -\item{cross}{a logical value indicating whether cross-validation should be applied or not. Default = FALSE} - -\item{memb}{a logical value indicating whether Clim() computes one climatology for each experimental data product member(TRUE) or -it computes one sole climatology for all members (FALSE). Default = TRUE} +\item{cross}{A logical value indicating whether cross-validation should be applied or not. Default = FALSE.} -\item{kharin}{a logical value indicating whether Clim() applies Kharin method or not. Default = FALSE} - -\item{NDV}{a logical value indicating whether Clim() applies Fukhar method or not. Default = FALSE} +\item{memb}{A logical value indicating whether Clim() computes one climatology for each experimental data +product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE.} } \value{ -an array +A list of length 2: +\itemize{ +\item\code{$var_exp}{ Array of anomalies for \code{var_exp} with same dimensions input.} +\item\code{$var_obs}{ Array of anomalies for \code{var_obs} with same dimensions input.}} } \description{ -This function computes the anomalies relative to a climatology computed along the starting dates -dimension allowing the application or not of crossvalidated climatologies. The computation is carried out -independently for experimental and observational data products. - -Climatologies are computed using the Clim() function that allows the combination of different methods -(kharin and NDV flags: see documentation) +This function computes the anomalies relative to a climatology computed along the +starting dates dimension allowing the application or not of crossvalidated climatologies. The +computation is carried out independently for experimental and observational data products. } \examples{ -## Example 1: -t <- 1 : 360 -obsS <- sin(t * pi / 180) -expS <- obsS + rnorm(360, 0., 0.1) +# Example 1: +Two arrays to compute anomalies: \\code{var_exp} and \\code{var_obs} + +dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100) +dim_obs <- c(mod = 1, member = 1, sdates = 2, ftime = 100) + +var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) +var_obs <- array(rnorm(prod(dim_obs)), dim = dim_obs) + +Computation of anomalies without cross-validation +anoM <- AnoMultiOptions(var_exp = var_exp, var_obs = var_obs, + cross = FALSE, memb = TRUE) +ano_exp <- anoM$var_exp +ano_obs <- anoM$var_obs + + +# Example 2: +One array to compute anomalie: \\code{var_exp} + +dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100, nlon = 100, nlat = 50) + +var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) + +Computation of anomalies with cross-validation using kharin and NDV methods +anoM <- AnoMultiOptions(var_exp , + cross = TRUE, memb = TRUE) +ano_exp <- anoM$var_exp + +} +\seealso{ +\code{\link{Ano_CrossValid, }} +\code{\link{Clim.}} } -- GitLab From 93753324eb16530b484590742c111b779cc8cea5 Mon Sep 17 00:00:00 2001 From: "jesus.pena@bsc.es" Date: Mon, 11 Feb 2019 12:58:59 +0100 Subject: [PATCH 03/16] Possibility of 1 or 2 inputs. Plus ftime climatology option git status pwd --- R/AnoMultiOptions.R | 145 +++++++++++++++++++++++++++++++------------- 1 file changed, 104 insertions(+), 41 deletions(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 3192e517..42273436 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -1,8 +1,9 @@ -#'Anomalies relative to a climatology along sdates with or without cross-validation +#'Anomalies relative to a climatology along selected dimension with or without cross-validation #' #'@description This function computes the anomalies relative to a climatology computed along the -#'starting dates dimension allowing the application or not of crossvalidated climatologies. The -#'computation is carried out independently for experimental and observational data products. +#'selected dimension (usually starting dates or forecast time) allowing the application or not of +#'crossvalidated climatologies. The computation is carried out independently for experimental and +#'observational data products. #' #'@param var_exp A multidimensional array related to experimental simulations with the specific #'dimensions required by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) @@ -17,6 +18,9 @@ #'@param memb A logical value indicating whether Clim() computes one climatology for each experimental data #'product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE. #' +#'#'@param dim_ano An integer indicating the dimension along which the climatology will be computed. It +#'usually corresponds to 3 (sdates) or 4 (ftime). Default = 3. +#' #' #'@return A list of length 2: #'\itemize{ @@ -27,32 +31,32 @@ #'@examples #' #'# Example 1: -#'Two arrays to compute anomalies: \code{var_exp} and \code{var_obs} #' -#'dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100) -#'dim_obs <- c(mod = 1, member = 1, sdates = 2, ftime = 100) +#'dim_exp=c(mod = 1, member = 3, sdates = 2, ftime = 100) +#'dim_obs=c(mod = 1, member = 1, sdates = 2, ftime = 100) #' -#'var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) -#'var_obs <- array(rnorm(prod(dim_obs)), dim = dim_obs) +#'var_exp = array(rnorm(prod(dim_exp)), dim=dim_exp) +#'var_obs = array(rnorm(prod(dim_obs)), dim=dim_obs) #' #'Computation of anomalies without cross-validation #'anoM <- AnoMultiOptions(var_exp = var_exp, var_obs = var_obs, #' cross = FALSE, memb = TRUE) -#'ano_exp <- anoM$var_exp -#'ano_obs <- anoM$var_obs +#'str(anoM) #' #' #'# Example 2: -#'One array to compute anomalie: \code{var_exp} #' -#'dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100, nlon = 100, nlat = 50) +#'dim_exp=c(mod = 1, member = 3, sdates = 2, ftime = 100, nlon = 100, nlat = 50) +#'dim_obs=c(mod = 1, member = 1, sdates = 2, ftime = 100, nlon = 100, nlat = 50) #' -#'var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) +#'var_exp = array(rnorm(prod(dim_exp)), dim=dim_exp) +#'var_obs = array(rnorm(prod(dim_obs)), dim=dim_obs) #' #'Computation of anomalies with cross-validation using kharin and NDV methods -#'anoM <- AnoMultiOptions(var_exp , +#'anoM <- AnoMultiOptions(var_exp = var_exp, var_obs = var_obs, #' cross = TRUE, memb = TRUE) -#'ano_exp <- anoM$var_exp +#'str(anoM) +#'anoM <- AnoMultiOptions2(var_exp = var_exp) #' #'@seealso #'\code{\link{Ano_CrossValid, }} @@ -61,57 +65,87 @@ #' #'@export -AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb = TRUE) { - - dimexp <- dim(var_exp) - dimobs <- dim(var_obs) +AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb = TRUE, dim_ano = 3) { + case_exp = case_obs = 0 if (is.null(var_exp) & is.null(var_obs)) { - stop("Not-NULL array var_exp or var_obs or both should be provided.") + stop("Parameter 'var_exp' and 'var_obs' cannot be NULL.") } if (is.null(var_exp)) { var_exp <- var_obs - dimexp <- dim(var_exp) - warning('No experimental data (var_exp) provided.') + case_obs = 1 + warning("Parameter 'var_exp' is not provided and will be recycled.") } if (is.null(var_obs)) { var_obs <- var_exp - dimobs <- dim(var_obs) - warning('No observational data provided (var_obs).') + case_exp = 1 + warning("Parameter 'var_obs' is not provided and will be recycled.") } - if (length(dimexp) < 4 | length(dimobs) < 4) { + + dim_exp <- dim(var_exp) + dim_obs <- dim(var_obs) + + + if (length(dim_exp) < 4 | length(dim_obs) < 4) { stop("Invalid input dimensions: At least 4 dimensions needed for var_exp or var_obs, c(nexp/nobs, nmemb, nsdates, nltime).") } - if (length(dimexp) != length(dimobs)) { + if (length(dim_exp) != length(dim_obs)) { stop('Invalid input dimensions: var_exp and var_obs should have same number of dimensions.') } - for (jn in 3:max(length(dimexp), length(dimobs))) { - if (dimexp[jn] != dimobs[jn]) { + for (jn in 3:max(length(dim_exp), length(dim_obs))) { + if (dim_exp[jn] != dim_obs[jn]) { stop("Invalid input dimensions: var_exp and var_obs dimensions should have same length for nsdates and nltime (and for nlevel,nlon and nlat if present).") } } - if (dimexp[3] == 1 | dimobs[3] == 1) { + + if (dim_exp[3] == 1 | dim_obs[3] == 1) { stop("Invalid input dimensions: sdate dimension length should be greater than 1") } if (!is.logical(cross) | !is.logical(memb) ) { stop("Invalid flags used: Flags 'cross' and 'memb' should be all logical.") } if (length(cross) > 1 | length(memb) > 1 ) { - warning("Invalid flags used: Flags 'cross' and 'memb' should all have length 1.") + stop("Invalid flags used: Flags 'cross' and 'memb' should all have length 1.") + } + if ( !(dim_ano == 3 | dim_ano == 4)) { + warning("A not-temporal dimension selected, check 'dim_ano' parameter.") + } + if (dim_exp[dim_ano] < 2) { + stop("Invalid length for selected anomaly time dimension 'dim(var_exp[dim_ano])': length should be greater than 1.") + } + + # Selecting time dimension through dimensions permutation + if (dim_ano != 3) { + dimnames_data <- names(dim_exp) + dimperm <- 1 : length(dim_exp) + dimperm[3] <- dim_ano + dimperm[dim_ano] <- 3 + var_exp <- aperm(var_exp, perm = dimperm) + var_obs <- aperm(var_obs, perm = dimperm) + + #Updating permuted dimensions + dim_exp <- dim(var_exp) + dim_obs <- dim(var_obs) } # Computating anomalies - #--------------------- + #---------------------- # With cross-validation if (cross) { ano <- Ano_CrossValid(var_exp, var_obs, memb = memb) - - # Without cross-validation + if (case_obs == 1) { + ano <- list(ano_obs = ano$ano_obs) + } + else if (case_exp == 1) { + ano <- list(ano_exp = ano$ano_exp) + } + + # Without cross-validation } else { tmp <- Clim(var_exp , var_obs , memb ) @@ -119,22 +153,51 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb clim_exp <- tmp$clim_exp clim_obs <- tmp$clim_obs } else { - clim_exp <- InsertDim(tmp$clim_exp, 2, dimexp[2]) - clim_obs <- InsertDim(tmp$clim_obs, 2, dimobs[2]) + clim_exp <- InsertDim(tmp$clim_exp, 2, dim_exp[2]) + clim_obs <- InsertDim(tmp$clim_obs, 2, dim_obs[2]) } - clim_exp <- InsertDim(clim_exp, 3, dimexp[3]) - clim_obs <- InsertDim(clim_obs, 3, dimobs[3]) + + clim_exp <- InsertDim(clim_exp, 3, dim_exp[3]) + clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) ano_exp <- var_exp - clim_exp ano_obs <- var_obs - clim_obs - - ano <- list(ano_exp = ano_exp, ano_obs = ano_obs) + + if (case_obs == 1) { + ano <- list(ano_obs = ano_obs) + } + else if (case_exp == 1) { + ano <- list(ano_exp = ano_exp) + } + else { + ano <- list(ano_exp = ano_exp, ano_obs = ano_obs) + } } - + + # Permuting back dimensions to original order + if (dim_ano != 3) { + # var_exp <- aperm(var_exp, perm = dimperm) + # var_obs <- aperm(var_obs, perm = dimperm) + # + # attr(var_exp, 'dimensions') <- dimnames_data + # attr(var_obs, 'dimensions') <- dimnames_data + + ano$ano_exp <- aperm(ano$ano_exp, perm = dimperm) + ano$ano_obs <- aperm(ano$ano_obs, perm = dimperm) + + #Updating back permuted dimensions + dim_exp <- dim(var_exp) + dim_obs <- dim(var_obs) + } + + #Adding dimensions names + attr(ano$ano_exp, 'dimensions') <- dimnames_data + attr(ano$ano_obs, 'dimensions') <- dimnames_data + # Outputs # ~~~~~~~~~ invisible(ano) } - \ No newline at end of file + -- GitLab From 984918bb29b4af67d15e9c91e73ed58bb12852bf Mon Sep 17 00:00:00 2001 From: "jesus.pena@bsc.es" Date: Mon, 11 Feb 2019 17:29:54 +0100 Subject: [PATCH 04/16] length anomaly dimension greater than 1 --- R/AnoMultiOptions.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 42273436..219abfc3 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -7,7 +7,7 @@ #' #'@param var_exp A multidimensional array related to experimental simulations with the specific #'dimensions required by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) -#'up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon) +#'up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon). #' #'@param var_obs A multidimensional array related to observations with the specific dimensions required #'by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to @@ -100,8 +100,8 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb } } - if (dim_exp[3] == 1 | dim_obs[3] == 1) { - stop("Invalid input dimensions: sdate dimension length should be greater than 1") + if (dim_ano == dim_ano & (dim_exp[dim_ano] == 1 | dim_obs[dim_ano] == 1)) { + stop("Invalid input dimensions: selected dim_ano dimension length should be greater than 1") } if (!is.logical(cross) | !is.logical(memb) ) { stop("Invalid flags used: Flags 'cross' and 'memb' should be all logical.") -- GitLab From df37b0a4e5d107205d240ca7227dcb2709e8c5a1 Mon Sep 17 00:00:00 2001 From: jpena Date: Tue, 12 Feb 2019 16:44:51 +0100 Subject: [PATCH 05/16] New parameter for dimension selection added --- R/AnoMultiOptions.R | 21 +++++++++++++++------ man/AnoMultiOptions.Rd | 38 +++++++++++++++++++++----------------- 2 files changed, 36 insertions(+), 23 deletions(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 219abfc3..4bc70a3c 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -18,7 +18,7 @@ #'@param memb A logical value indicating whether Clim() computes one climatology for each experimental data #'product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE. #' -#'#'@param dim_ano An integer indicating the dimension along which the climatology will be computed. It +#'@param dim_ano An integer indicating the dimension along which the climatology will be computed. It #'usually corresponds to 3 (sdates) or 4 (ftime). Default = 3. #' #' @@ -122,6 +122,7 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb dimperm <- 1 : length(dim_exp) dimperm[3] <- dim_ano dimperm[dim_ano] <- 3 + var_exp <- aperm(var_exp, perm = dimperm) var_obs <- aperm(var_obs, perm = dimperm) @@ -183,8 +184,13 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb # attr(var_exp, 'dimensions') <- dimnames_data # attr(var_obs, 'dimensions') <- dimnames_data - ano$ano_exp <- aperm(ano$ano_exp, perm = dimperm) - ano$ano_obs <- aperm(ano$ano_obs, perm = dimperm) + if (case_obs == 0) { + ano$ano_exp <- aperm(ano$ano_exp, perm = dimperm) + } + if (case_exp == 0) { + ano$ano_obs <- aperm(ano$ano_obs, perm = dimperm) + } + #Updating back permuted dimensions dim_exp <- dim(var_exp) @@ -192,9 +198,12 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb } #Adding dimensions names - attr(ano$ano_exp, 'dimensions') <- dimnames_data - attr(ano$ano_obs, 'dimensions') <- dimnames_data - + if (case_obs == 0) { + attr(ano$ano_exp, 'dimensions') <- dimnames_data + } + if (case_exp == 0) { + attr(ano$ano_obs, 'dimensions') <- dimnames_data + } # Outputs # ~~~~~~~~~ invisible(ano) diff --git a/man/AnoMultiOptions.Rd b/man/AnoMultiOptions.Rd index 4cd26605..8912530c 100644 --- a/man/AnoMultiOptions.Rd +++ b/man/AnoMultiOptions.Rd @@ -2,15 +2,15 @@ % Please edit documentation in R/AnoMultiOptions.R \name{AnoMultiOptions} \alias{AnoMultiOptions} -\title{Anomalies relative to a climatology along sdates with or without cross-validation} +\title{Anomalies relative to a climatology along selected dimension with or without cross-validation} \usage{ AnoMultiOptions(var_exp = NULL, var_obs = NULL, cross = FALSE, - memb = TRUE) + memb = TRUE, dim_ano = 3) } \arguments{ \item{var_exp}{A multidimensional array related to experimental simulations with the specific dimensions required by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) -up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} +up to c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon).} \item{var_obs}{A multidimensional array related to observations with the specific dimensions required by Clim(): i.e. \code{dim(var_exp)} goes from c(nmod/nexp, nmemb/nparam, nsdates, nltime) up to @@ -20,6 +20,9 @@ c(nmod/nexp, nmemb/nparam, nsdates, nltime, nlevel, nlat, nlon)} \item{memb}{A logical value indicating whether Clim() computes one climatology for each experimental data product member(TRUE) or it computes one sole climatology for all members (FALSE). Default = TRUE.} + +\item{dim_ano}{An integer indicating the dimension along which the climatology will be computed. It +usually corresponds to 3 (sdates) or 4 (ftime). Default = 3.} } \value{ A list of length 2: @@ -29,38 +32,39 @@ A list of length 2: } \description{ This function computes the anomalies relative to a climatology computed along the -starting dates dimension allowing the application or not of crossvalidated climatologies. The -computation is carried out independently for experimental and observational data products. +selected dimension (usually starting dates or forecast time) allowing the application or not of +crossvalidated climatologies. The computation is carried out independently for experimental and +observational data products. } \examples{ # Example 1: -Two arrays to compute anomalies: \\code{var_exp} and \\code{var_obs} -dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100) -dim_obs <- c(mod = 1, member = 1, sdates = 2, ftime = 100) +dim_exp=c(mod = 1, member = 3, sdates = 2, ftime = 100) +dim_obs=c(mod = 1, member = 1, sdates = 2, ftime = 100) -var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) -var_obs <- array(rnorm(prod(dim_obs)), dim = dim_obs) +var_exp = array(rnorm(prod(dim_exp)), dim=dim_exp) +var_obs = array(rnorm(prod(dim_obs)), dim=dim_obs) Computation of anomalies without cross-validation anoM <- AnoMultiOptions(var_exp = var_exp, var_obs = var_obs, cross = FALSE, memb = TRUE) -ano_exp <- anoM$var_exp -ano_obs <- anoM$var_obs +str(anoM) # Example 2: -One array to compute anomalie: \\code{var_exp} -dim_exp <- c(mod = 1, member = 3, sdates = 2, ftime = 100, nlon = 100, nlat = 50) +dim_exp=c(mod = 1, member = 3, sdates = 2, ftime = 100, nlon = 100, nlat = 50) +dim_obs=c(mod = 1, member = 1, sdates = 2, ftime = 100, nlon = 100, nlat = 50) -var_exp <- array(rnorm(prod(dim_exp)), dim = dim_exp) +var_exp = array(rnorm(prod(dim_exp)), dim=dim_exp) +var_obs = array(rnorm(prod(dim_obs)), dim=dim_obs) Computation of anomalies with cross-validation using kharin and NDV methods -anoM <- AnoMultiOptions(var_exp , +anoM <- AnoMultiOptions(var_exp = var_exp, var_obs = var_obs, cross = TRUE, memb = TRUE) -ano_exp <- anoM$var_exp +str(anoM) +anoM <- AnoMultiOptions2(var_exp = var_exp) } \seealso{ -- GitLab From 9670dac77156e43084828a7bbd6398edd14614e8 Mon Sep 17 00:00:00 2001 From: jpena Date: Wed, 13 Feb 2019 10:30:57 +0100 Subject: [PATCH 06/16] Removing inconsistency with dimnames_data --- R/AnoMultiOptions.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 4bc70a3c..7baaf2ab 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -84,6 +84,7 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb dim_exp <- dim(var_exp) dim_obs <- dim(var_obs) + dimnames_data <- names(dim_exp) if (length(dim_exp) < 4 | length(dim_obs) < 4) { @@ -118,7 +119,6 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb # Selecting time dimension through dimensions permutation if (dim_ano != 3) { - dimnames_data <- names(dim_exp) dimperm <- 1 : length(dim_exp) dimperm[3] <- dim_ano dimperm[dim_ano] <- 3 -- GitLab From 5eb9a907d2ffa0dfab8d4a5bacc0ad5f64a0dbbd Mon Sep 17 00:00:00 2001 From: jpena Date: Wed, 13 Feb 2019 17:53:14 +0100 Subject: [PATCH 07/16] minor changes --- R/AnoMultiOptions.R | 5 ----- 1 file changed, 5 deletions(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 7baaf2ab..98a50e49 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -178,11 +178,6 @@ AnoMultiOptions <- function(var_exp = NULL, var_obs = NULL, cross = FALSE, memb # Permuting back dimensions to original order if (dim_ano != 3) { - # var_exp <- aperm(var_exp, perm = dimperm) - # var_obs <- aperm(var_obs, perm = dimperm) - # - # attr(var_exp, 'dimensions') <- dimnames_data - # attr(var_obs, 'dimensions') <- dimnames_data if (case_obs == 0) { ano$ano_exp <- aperm(ano$ano_exp, perm = dimperm) -- GitLab From bab1e4e83505a1a6222211f283596388b4b9a2b7 Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 09:49:50 +0100 Subject: [PATCH 08/16] Fixing CRAN inconsistencies --- R/AnoMultiOptions.R | 4 ++-- man/AnoMultiOptions.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/AnoMultiOptions.R b/R/AnoMultiOptions.R index 98a50e49..8ec1fb5c 100644 --- a/R/AnoMultiOptions.R +++ b/R/AnoMultiOptions.R @@ -59,8 +59,8 @@ #'anoM <- AnoMultiOptions2(var_exp = var_exp) #' #'@seealso -#'\code{\link{Ano_CrossValid, }} -#'\code{\link{Clim.}} +#'\code{\link{Ano_CrossValid}} +#'\code{\link{Clim}} #' #' #'@export diff --git a/man/AnoMultiOptions.Rd b/man/AnoMultiOptions.Rd index 8912530c..8c5fa4f5 100644 --- a/man/AnoMultiOptions.Rd +++ b/man/AnoMultiOptions.Rd @@ -68,7 +68,7 @@ anoM <- AnoMultiOptions2(var_exp = var_exp) } \seealso{ -\code{\link{Ano_CrossValid, }} -\code{\link{Clim.}} +\code{\link{Ano_CrossValid}} +\code{\link{Clim}} } -- GitLab From 69bee51e1a1f0532be3538c82884cd904a0c3288 Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 12:05:55 +0100 Subject: [PATCH 09/16] Loading vignette for AnoMultiOptions --- vignettes/Anomalies.md | 272 +++++++++++++++++++++++ vignettes/DeseasonalizedTASanomalies.png | Bin 0 -> 46373 bytes vignettes/Nino3.4TimeSeries.png | Bin 0 -> 59920 bytes 3 files changed, 272 insertions(+) create mode 100644 vignettes/Anomalies.md create mode 100644 vignettes/DeseasonalizedTASanomalies.png create mode 100644 vignettes/Nino3.4TimeSeries.png diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md new file mode 100644 index 00000000..3c9d9f84 --- /dev/null +++ b/vignettes/Anomalies.md @@ -0,0 +1,272 @@ +--- +author: "Suso" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::knitr} + %\VignetteIndexEntry{anomaly_agreement} + \usepackage[utf8]{inputenc} +--- +Anomalies +======== + +Function AnoMultiOption() allows the user different possibilities of computing anomalies for observational and modelled data. This vignette illustrates a complete example of how to load datasets using the Load() function following different cases for anomalies computations. + + +### 1- Load dependencies + +This example requires the following system libraries: + +- libssl-dev +- libnecdf-dev +- cdo + + +The **s2dverification R package** should be loaded by running the following lines in R, onces it is integrated into CRAN mirror. + +```r +library(s2dverification) +library(devtools) +``` + +You will need to install specific versions of certain R packages as follows: + +```r +source('https://earth.bsc.es/gitlab/es/s2dverification/blob/develop-anom/R/AnoMultiOptions.R') +``` + + + +### 2- Define the problem and the correspondent data and parameters + +We want to show how well seasonal forecast models perform for up to a 7-months forecast window. The idea is to compare a observational dataset represented by the ERAinterim reanalysis with the EUROSIP multi-model seasonal forecasting system. + +For more information about both datasets see the next documentation: +- [**ERAinterim**](https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim). +- [**EUROSIP system**](https://www.ecmwf.int/en/forecasts/documentation-and-support/long-range/seasonal-forecast-documentation/eurosip-user-guide/multi-model). + +Both datasets can be downloaded from the corresponding server. However in this example we will use the BSC esarchive database. Files path parameters for the Load() function are defined as follows: + +```r +#Path for experimental data files path +path_exp_EUROSIP <- list(name = 'system3_m1', + path = file.path('/esarchive/exp/EUROSIP/ecmwf', + '$EXP_NAME$/monthly_mean', + '$VAR_NAME$_f6h/$VAR_NAME$_$START_DATE$.nc')) + +#Path for observational data files +path_obs_eraI <- list(name = 'erainterim', + path = file.path('/esarchive/recon/ecmwf', + '$OBS_NAME$/monthly_mean', + '$VAR_NAME$/$VAR_NAME$_$YEAR$$MONTH$.nc')) +``` + +Our case of study will be a proxy of El Niño3.4 index defined here as the 5-month moving average temperature air surface (tas) anomalies averaged for the region comprised between longitudes [170ºW - 120ºW] and latitudes [5ºS - 5ºN]. We will focus on the 1990s and 2000s decades and especially on the strong El Niño signal during 1997/1998 and the strong La Niña during 1999/2000. Note that the accepted definition of El Niño3.4 index uses sea surface temperature (SST) anomalies instead of tas used in our proxy. Although there are some discrepancies between both, our index shows the same general pattern and it therefore suits for this practical example. + +```r +#Predefined domain for El Niño3.4 index +lonNino <- c(360 - 170, 360 - 120) +latNino <- c(-5, 5) + +#Selected domain for map visualization is defined as: +lonMap <- c(-250, -60) +latMap <- c(-30, 30) + +``` + + +### 3- Loading observational data + +As a first step we will explore the observational data for the last decades. + +In the case of using Load() for loading only observational data, all available data from the +given starting date will be loaded as a time series along the forecast time dimension. Anomalies +should therefore be computed using along this dimension using parameter `dim_ano = 4`. Parameter `'output'` +should be set to `'areave'` to compute El Niño3.4 domain averages and to `'lonlat'` to load snap-shot maps. + + +```r +#Starting dates parameter for observations +sdates_obs <- '19900101' + +#Loading data +data_Nino <- Load(var = 'tas', + obs = list(path_obs_eraI), + sdates = sdates_obs, + output = 'areave', + lonmin = lonNino[1], lonmax = lonNino[2], + latmin = latNino[1], latmax = latNino[2]) + +data_map <- Load(var = 'tas', + obs = list(path_obs_eraI), + sdates = sdates_obs, + output = 'lonlat', + lonmin = lonMap[1], lonmax = lonMap[2], + latmin = latMap[1], latmax = latMap[2]) + + +#Retrieving observational data +obs_Nino <- data_Nino$obs +obs_map <- data_map$obs + +#Selecting temporal period from 1990 to 2010 +i1990to2010 <- seq(1,12*20) +obs_Nino <- array(obs_Nino[,,,i1990to2010], dim = c(1,1,1,length(i1990to2010))) + +#Retrieving dates +dates_obs <- as.Date(data_Nino$Dates$start[i1990to2010]) + + +#Computing anomalies along forecast time dimension +ano_obs_Nino <- AnoMultiOptions(var_obs = obs_Nino, dim_ano = 4)$ano_obs +ano_obs_map <- AnoMultiOptions(var_obs = obs_map, dim_ano = 4)$ano_obs + +#Running 5-months moving window average +ano_obs_Nino_sm <- Smoothing(ano_obs_Nino, runmeanlen = 5, numdimt = 4) + +#Selecting november-1997 and november-1999 snap-shot map (index starting from january 1990) +inov97 <- 12*7 + 11 +inov99 <- 12*9 + 11 + + +#Plotting the nov1997 and nov1999 anomalies maps +taslim <- c(-2, 2) +brks <- seq(taslim[1], taslim[2], length.out = 21) +layout(matrix(c(1, 1, + 2, 3, + 4, 4), + 3, 2, byrow = TRUE), + heights = c(4, 6, 1), + widths = c(1, 1)) + +plot(x = dates_obs, y = ano_obs_Nino[,,,],type='l',col='grey', + xlab = '', ylab = '', + xlim= c(dates_obs[1], dates_obs[20*12]), ylim=taslim) +lines(x = dates_obs, y = ano_obs_Nino_sm[,,,],col='darkgreen') +lines(x = dates_obs, y = rep(0,length(dates_obs)), col = 'black') +lines( rep(dates_obs[inov97], 2), taslim, lwd = 4, col = 'darkred') +lines(rep(dates_obs[inov99],2), taslim, lwd = 4, col = 'darkblue') +grid() +title('Tas anomaly between 1990 and 2010 for Niño3.4 domain', line = 1.5) + +PlotEquiMap(ano_obs_map[1, 1, 1, inov97,,], data_map$lon, data_map$lat, + brks = brks, drawleg = F, toptitle = 'Tas anomaly (El Niño nov-97)', sizetit = 0.8) +rect(lonNino[1], latNino[1], lonNino[2], latNino[2]) + +PlotEquiMap(ano_obs_map[1, 1, 1, inov99,,], data_map$lon, data_map$lat, + brks = brks, drawleg = F, + toptitle = 'Tas anomaly (La Niña nov-99)', sizetit = 0.8) +rect(lonNino[1], latNino[1], lonNino[2], latNino[2]) + +ColorBar(brks, vert = FALSE, subsampleg = 5) + +```r + +The strong signal of El Niño / La Niña (positive / negative anomalies) for november 1997 and 1999 can be clearly seen in both time-series and snap-shot maps. However, since the anomaly is computed with the annual climatological mean, the clear seasonality presented in the time series is not removed hence obscuring the real signal of El Niño variability. We can clarify this problem computing monthly anomalies as described in the next section. + + + +### 4- Seasonal forecast versus observational data + + +In this section we will compare the tas anomalies for the Niño3.4 between the observational ERAinterim data and the EUROSIP forecast system. The latter provide a 7-month forecast window for each starting date, therefore the aim of this section will be to show how this forecast performs during the two events seen in section 3, El Niño 1997 and La Niña 1999. + +We will also compare two different methods when calculating anomalies, using or not using cross-validation during the climatology computation. + +```r +#Starting dates parameter for time series +sdates_exp_df <- seq(as.Date("1994/11/01"), as.Date("2002/11/01"), "years") +sdates_exp <- format(sdates_exp_df,'%Y%m%d') + +#Loading data +data <- Load(var = 'tas', + exp = list(path_exp_EUROSIP), + obs = list(path_obs_eraI), + sdates = sdates_exp, + output = 'areave', + lonmin = lonNino[1], lonmax = lonNino[2], + latmin = latNino[1], latmax = latNino[2]) + + +#Retrieving experimental and observational data +obs <- data$obs +exp <- data$mod + +#Retrieving dates +dates_exp <- as.Date(data$Dates$start) + +#Selecting november-1997 and november-1999 snap-shot map (index starting from january 1990) +inov97 <- 12*7 + 11 +inov99 <- 12*9 + 11 + +#Computing anomalies along the starting date dimension without cross-validation +ano_NoCV <- AnoMultiOptions(var_exp = exp, var_obs = obs, dim_ano = 3, cross = FALSE) +ano_obs_NoCV <- ano$ano_obs +ano_exp_NoCV <- ano$ano_exp + +#Computing anomalies along the starting date dimension with cross-validation +ano_CV <- AnoMultiOptions(var_exp = exp, var_obs = obs, dim_ano = 3, cross = TRUE) +ano_obs_CV <- ano$ano_obs +ano_exp_CV <- ano$ano_exp + +#Running 5-months moving window average for obs and averaging across experimental member +ano_exp_NoCV_m <- Mean1Dim(var = ano_exp_NoCV_sm, posdim = 2 ) +ano_exp_CV_m <- Mean1Dim(var = ano_exp_NoCV_sm, posdim = 2 ) + + +#Plotting +taslim <- c(-2, 2) + +layout(matrix(c(1, 1, 2, 2),4,1, byrow = TRUE)) + +plot(x = dates_obs, y = rep(0,length(dates_obs)), + type = 'l', col = 'black', + xlab = 'Forecast time', ylab = 'tas anomalies', + xaxt = 'n', + xlim= c(dates_obs[inov97-40], dates_obs[inov99+40]), ylim=taslim) +axis.Date(side = 1, at = sdates_exp_df,format = '%h%Y') + +for (isd in 1:dim(exp)[3]) { + jsd = (isd-1)*7+1 + for (im in 1:dim(exp)[2]) { + lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_NoCV[1, im, isd, ], + col = 'grey') + } + lines(x = dates_exp[jsd : (jsd+6)], y = ano_obs_NoCV[1,1,isd,], + col = 'black', lwd = 2) + lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_NoCV_m[1,isd,], + col = 'green', lwd = 2) +} +grid() +text(sdates_exp_df[5], 2, 'Without CrossValidation', font = 2, adj = 0) +title('Deseasonalized tas anomalies for observations and forecast') + + + +plot(x = dates_obs, y = rep(0,length(dates_obs)), + type = 'l', col = 'black', + xlab = 'Forecast time', ylab = 'tas anomalies', + xaxt = 'n', + xlim= c(dates_obs[inov97-40], dates_obs[inov99+40]), ylim=taslim) +axis.Date(side = 1, at = sdates_exp_df,format = '%h%Y') + +for (isd in 1:dim(exp)[3]) { + jsd = (isd-1)*7+1 + for (im in 1:dim(exp)[2]) { + lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_CV[1, im, isd, ], + col = 'grey') + } + lines(x = dates_exp[jsd : (jsd+6)], y = ano_obs_CV[1,1,isd,], + col = 'black', lwd = 2) + lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_CV_m[1,isd,], + col = 'green', lwd = 2) +} +grid() +text(sdates_exp_df[5],2,'With CrossValidation',font = 2, adj = 0) + +```r + +The figure above shows the deseasonalized tas anomalies for El Nino3.4 domain comparing observational ERAinterim data (black colour lines) with 7-months forecasts from the EUROSIP system for several starting dates (different members [grey] and ensemble mean [green]). The lower / upper panels corresponds with the application / not-application of CrossValidation in the computation of the climatology used in the anomaly calculation. It is clearly seen that there is almost no difference in these two cases due to the significantly high number of sdates (9) in the simulation. + +Overall the agreement between observations and simulations is quite good, especially because in opposition to section 3 this anomalies does not include any seasonality. Since the anomaly is computed along the starting dates dimension, i.e. from the corresponding monthly means of the two decades, it assures that the seasonality is removed. This means that changes in the deseasonalized tas anomalies represented in the figure and thus directly linked to natural variability, like El Niño / La Niña signal, are reasonably well reproduced by the EUROSIP system. + diff --git a/vignettes/DeseasonalizedTASanomalies.png b/vignettes/DeseasonalizedTASanomalies.png new file mode 100644 index 0000000000000000000000000000000000000000..b863eaacd5a3bd8113b71f4829a6b4523b61d13e GIT binary patch literal 46373 zcmcG$bySsW)IGXE6r@E;T0lfvL_iv(LsGg$q)SRtNdYMlX^;|-?vQkYfPj>AmvndE z#X0Bu-S7AB9pmmX^c*+r_kG?c)|zY1xpt7!GbtQQGE4*lfg>X=u7W_IOke%MK!ta9 zDjxhmAZQRW;-YGqijzW6%Qbd2~%)ehAa&U0y5_XEdFf=r*VvP0u z+1S`)IMRXt>lu@>%^Th}JMlI(!VuhVAXSR(v!fLPFk7uCd~36%OohiYI+TLqmVl zDy_$eNlDFi*C|s=qf4SH@J27xe-u%7it-1mWC1%xbHZAq_ z@xi6wMId+`=A}(c=JT{lmmE*dj(3enD}FT@KY#P4y4?g7QM>(hm^bRTjaV4x^5wSzKo8J9)y--2yQT<{Q@k>#iWfgW#0k1{BMcXRCsE)Ab}v>Ie9C;o8okmsS5dt^_NIr@>RxWV zMfqs;>qXAE&e+)4hRciZ=P}JHBfq_L+3{k1oe8(Mw_#o;>f9^UnCk53ehwLU>lqlN z#Kg3kmd3E^IB!ii1qTPu&dxSBH&c5ZHLX+z;!<)ZbX8YZW8ZjjxG|v~wV+#>QzlE7WsaH@Rb^M`)NkJn zL*{Z%CwAQ`I~ZNEpu5naLWNpIzrj$)wGZ#OW7rLDAP7lFkQ`Rr2{+Mx{7&Kz#J(ez z)&EpXY}l$ID{HaI7aau=6cjX8X{(Vtpx*E}p^JoF-(xY!C0otH!J*{&Y~_@f`*hAbwWjCdO#AOF90-AM5Z^zc?+BiH_~3QbK-9N2|Am27Nm z8`aSk4Uw!2L?k41(QX!UDvX3l&xTIkRVY+At*R3QW%QbrX|u-33%PC?{&(E6Q_-o&9trCAGB>a|7pmfo~QZn^{ z7n8MbA|oT~-1bosn)d=OpGXqiVvKml5h?9QgAsUoiEa`7RHsaP-RC;i4Q$dPlXe2c z+h`Qi(mOsh2oZixvr|bSX4oz%&ne|*ol8p6&LwRUE^;LL&A14H_Qp*(uywd-(|d<6 zu30BpYiep*S*_K3o;DkASJ}@c>-_gnba*xA@K8DJCSrZ1BWp)TN4WL$q}`~md&1JH zs$db$k$A#;dwUd59OmDmT-%te8~gKT4~}^K>Bje5^?W|3l|ReO^n01`&ugb~ERoYG!s4+S5B}Jp0 zpB^g0JBo_J4)blQwnNK~>*OQ$a(pUDxTXVMbYr@Es!?~a1CW7u)B>(bxdR8AQw`@Z zK{NC1QQjxZS$CfX*zqQ2smo^dJN=SQAiu1}JE^{*IZ&pUA_O;g(bu0Gvu%KQH62z< zpU)g(q%h<(X`mU=ON)3bLeu>7r+&4=^#G$F@tu|2LQc{TaF9)3y;@vY*k2h45vM&q zJ*9r)(7aNarI>0S|~=@(i!L;vaKh zS6+(bm(;l*vLOsRKC@bnwf5D1mUWu99~(bf7g)lr`1o95z+xZ5&$IMih|>N2{bSxT zEYamN%+h9V6noLk8nSY7=s0BaXKDZZwA~jS5BmlZo!Y~2Y+Y~q3km`@+CieUDQDati^oyz z;=S~XU(}cc1q8NSJv7Ms%rjITBy??9F!P0~=UzuBF@7rG&gic)_azUN43$J0a@34Q zJWE$%G~l|z7-^HYis-J^LygS$w~G5lHai|E-tP<#6kLDu&!Fv_lEUE z&43AK@A4kC_gSwYM@0L4(Cz#6kYPV_L73l%<^iA1satF+|CvD)H73O2ITul%cam%YGk4)9JA&fB)+p zhiPZ*+S*!Qp9)#%=g3H9U0r#40y;W6BFq-ZB@p|4=sP+(u&}Wy9vb_#n|v!MpyDt@ zgtWvzdHuD?@Mu0tQRrle4k6;<;bCE6vHBf>c=bE+cqWkg;$*dOb*Rw$biML$bLw)Z zgIy|=Y#>*IJ=bGns$r_opkZlg$$G4!y#5)y6cnx_68T*=%7^qJuIO)0)**lV;I^Io zKzKJ&nh;aM+}s?}a(Y@C896x$qS9ufguXP^FQbg_!M=fwWMyT2 z{`@(&?w>z@%*wSBY&U;kWvXcD=}jz{4%?;nnpxE+buYR^Z!<;1^Z(w7>P?q=_3D*D zwZlt@>mws0Lj`*7Fa}K);vn(S(pQQQg5Xyo`#(XNjYvpvhnV@;a`5Q%l!TB_00L9E z@C64cwrnOrfKlnlINbH50MOR!#owV+f_$8uef|BV_4PvB-0EyGSpGjgeIm^D2%{2| zP4E4k_(WAh1A?g#J9{BSl}|A-R|ue`r3FR-duj+`TOTW7JAY3FO_OA0Q2`9eeS8L*OT~NY=(=BV0If|A5Ay;8j(s~U?J^vb;lPvV(vdnta`opj5HuB zDvHN?w9L2#bADaI%d#5F2o|{XRkvh8v}~q&t{&vu@0;~l$fD;$uzRp@ala}tGOOpA zdEy|q06e*DjJMC&Ir2?}lJgk9xO!bfMa4%(c*9O|o^C`T4fyV~deBPkUG~iM`K??B zjY9o8yBRUkfYZI@U3i_ezEt`1;q z^Nhz@ByZ;x0+Na8Y}AyY5T~_2+L~!NKhQ}L@^m}gIES1Ct6x)|Wc?Zg2dD1!Vwf6J zS9|+S9GsUNqyqy3w2|50zI|KQqm%ql<+jh1(#`(h!S>-{gVQQ80j9Va zJUy_Usy_#q_Sv#pH#sq})(9C;{G#=vbOt5AGYjB8z_aS={(gRNCbo|~gi`PjMY4mV z^tle}=AB`0Whj<>e0;;+zaK6$zs*VX$a!rT0`O?ctJA7|`C_bkG|HyrFFppE4Qi3! zFflO`H`B)??LScRe{MH9IX>RsoJwgjgd|HV93C3FuzMWCrdRUIXyI4%Gb6&A9;6Dg{({Q2R|X6^1XNU9L|%glPyV8lYq)Z7Vm4hz3*>yQ7! zCd9pYvsm|0TqndOX2M{A^(Y7ih4{TLo+;g0*INNb<=WQm>!sa704~TMTd0|tF-A&T zC9%-cXOzpRFQiR04p!|~F8e4`;Ex{%nexXydKDHP-VF#hNtaDB25~JmYw_WJ=Xv#M z!fFCw`KvhO@*%3FvFII?VKDrWTBR1&xi1fckjKL;_9BSU^bus%?jhQS7911^?G;ud zJ_Xk6hWG5NZrPA=`@&4le#RZm4^(D5B5*`=rw@xAja?4y??y{`aKMaToBple}&vDPc#O>6gbfgxFOG0Z)^jv=? z8^44C>>tG6zJ$OB`GjlT80FHQTD60 zx1)G_P@;Gf!@KBgNCN;$@YqgL5D`6=5P}5zSM{KU!~lw5J|i7F7zt-=^E)=1<$ouN zss%+9_1*MbxG*jj)(aUeeYZkm)9u73uYVws8u{7_ON@&D>mw6xqzgXmqEbcP&6K31 z7@|sLv^DRw=m~-b^L1jPqFyijB2w#gCx-K6GJZ1FOpoi>NRPYd`|JH!M+dXyE!9P3 zjR_S&$Uu#aB2tAKKj?Aq|F;ZeFnT;Dx{Zj42;P?z*wbR#chN*Qp>zbiv9q(Y)_wfX zvo^sc?d8zd2a${<9yj?=ndofDLW$FZn~+Del9A%HT4iRpLM6w_Eo07~2xDE1+~(Zc zm;jv*CTdsIvbl?)82QIYqB9L}1#aLYh%g3T#}bT$X&>L{{I3Q0c{S?6e@3l-OQ)0N z^RGaf&svahq8`AdnAl@pcOYTEi04F(n-WI?m?;V#;Mt>F_f`*@;JfD{3pCcj_bnP8-=FA1@O?N06 zZ~qr-1xwHR_vPGY+>>f#c4A`UpFe7t{yY7O!quOi!oIyK<%3mQL^;&j7{$kxG&N&m zVq)khxc+(IXqXEVPA#^W6n%HHk53=qJorD?B53`pkPs{Xw?OYZWCImd)g;s5*RW~P zs7n-!7^-5muto7_gCL(IMMwJ^-@^FkY1Lp#j!sVc%*uMEdwNEdl*^VNEa+L^LfEW) z8H?eSm6Zixu?wILn2M>XDQ|#ZzsvGHIfAaC1O>1e6Pr6}hx}4jQhMvQhy~~9!>KF< zrvH!b?!wemBh(kD?{6?yBUldc*Nz6>(><$ zdA9*?xYbPt-4MDhKP&Q<=Dq)m{iwa~3!ZYeJsl<`5fS(R{Qpq?vsRU#dF2_R!ow9q zYVSR}^@_|2QlfT^qxq2b0xZ22lCgBeW2q;VXUgLdFYoXUQX+)-F*ocSZEEmo&8>q2D%v|@Fd=Re;%*^r0(WaD_5X(UhIlA+%U3#Ht0+QRPhQrH0a0-)*aukZ~5 z7iu6fHO}1go~P}bDbDy04%M|#ZQt-knRp}i--X|Li+zwbRe52O9vF!J_Y1!QKa!C< zfs%qEF(rka5EE+oxAf#s9O$u;E}N511#yaS$mic z;-6is2fKE9+m-L#q~1u~Th)lzSXT%tY-Rr}n=%C)PY)r3+yP<|l1HQ;pngG!z>lE{ z+!(K((|weJ%3xpi)%|eeX?$QrXU1ycsx%=}c?< zcpHTx<|2JxOq+1VjaU`r`v6%=dPi@>P8QPXeu z>$R}v@C`k>doi4x=l9qcc?<0)*_SU~q=%Fl#qd#_}*akVLG zvJ9%VsF616)79w%$_Hj-G!Y>EJDN7rP1&OEC)W^-z~ue?ope&FQUj_p z?flQu^l2(kBKpv@NsRUok&X9~AP@kd$Fa$wL|_H#(o2js@@)Tt_pN1^c|mS2P(VQK zG*~}bUB}w;HzciC^zvI-Ss9G9UVrrJHw86y3X}n(1~xiW0TL)R48Um`4``LZZ)A-OoWBqn6TgcfVnPizqbLkj^^gO z_;>xxSq-B{E7u#aJgymJs@%JGFH$;#Ehe>?5tkVQ4fQt1sm!U~oStuxot+&ARd|EL z0;cxy3uTIDf3Lrv6~Vy3kp6U!xHvfNRMKU}O{_AJ-s;uYs_JSj+DKP7H;3wRo^O)~ zL?Sl7Id{T_X1@0@ibH|4=9utdvfk8OFEaPNCR%|K7cquU?T56Ms|tHc0)!+7O;0U{6fko(8CpCvIV` z*`$b`=S|L6nb+0K%Fb?U{9CmvskguqZ=C!_Z+c*_8E{nZZnJm*Nb4o$OTrcH$ve^V zAiJcRn_L?y?K0g3z{SMO{D_ltZLjz0pzD(cT|-Q7t_l~uF%-2!2~}AA_MGtQ`n(@z z?RbS4Qf}hhR#|%Blu-j=odCnH;Q14dpt+ps?W5k=-g>vGD3tr>=jU$(CjGHUvelTd zkxE!A3-KD7pZ*_Hbk^_aS zCn(Zn;fxcgcN3Tp7am$42bS@V#XfMyf9K|&(%>^OFX4?Ayo(n z3f{%QGz#d+k z`8fFil=d6SqnbGa^k`Xs8fn_d zJkxyK=YbpHBhOVo90me*wPU|k?e&;KxgM&ljRFKyXQH|%%1-;Iy($b7uf{B2ab}4y4nWc(k zBnxGTly;?pN0>PbtBN`#ku)DLPO=pmP&!W!djX{-NQ}Qq(!Ja>84ww+>!ls*a zHgH-I5$z@m3kz;_ULG5+*%+fTQ+Eg{Ot{FbPdTp2ms(c>bOcQVAt`(KQ0>=9=J^K? z9{es_xrtr}Q5N_ptI6795JMmbI)9~1N`px@se1~%_tU4>N&}9+Z9U&sD@OYa)t_X( zI{%;F6;(;i#a`zp!k1^ac`cjQmzUd3OM!y{&49;q)N$Nd_g?yrTE_z*I@`V~I6>w8 z?+6+JqpF-|NUEoA>P}q=cya22GhdjQm|)F8cx8#{*mkXbw7x``^brC87a8l_yLW|z zk|oo-G=B5oTB(ufvi*fys+mz6g=<4hg zKkP|4c=JOP$Q`mca@`V>_Ds2$L{srs%b#EmZ0+rRYBzBeNN%0O`GWoad9&E&(-FMV zV~>k;sqn<*^a;7yD~<;n3IDwDo6Npclf_@`77@{NU}w{NiR()v4zUVESMRDOKdVc>h%oEUF zkX|h;Et9TBwZ;*QocS|2C`(TO6A&lQEypMKml+w_2pLeZaA@)G-!umcQ7M?MIanNB zHlDk8QDDVjUxkN`VPj&Bk`~`m^Yu^Q8BzEjFTenn_37Hx((yi<%`Pm& z!@>#)31MbsE#K>ocrmZ^pqLWkzesxT*FTC)*J`J3G|aVzcPu!Z%B^82@-Bf|L@n%9 zo0Ie4@#B4{T7DNd-r!@l2Ug$D4=qPcTTLxvz=EWQgab$Zx5aN@U4g}O;?o21@U=&o zy9@=KcMw6ja09PBqp(LuL&NYFg9r3&`l6ou{kqF)!eC{!)b6`hbgMR;C|(M5Y@~Gq z8M@3>d4ojyiqO!Esi>=|!Nd>~5iQKmOJ}N#{k@-OCI!dO_05}mTLb~Y0o?LtX798U z892_mMW^aLq4@a?+!9b;?d|P2I5_`}7=d^%ZY(Y??&!#cCc00wKS$Oo6CjtNwyY?8 z{^z@kZMn9drsfijrOuorft0{T4Y_jwn%G2ua>FLQ*qz0B;ODm`*N4+Yt*?xWWF@#@{Yl~r3b#-;X!0H#` zVPZlK|8hOAP_v=WG80G3h#>I1ZDD6+#p0ca;r{DLX|N{1kTHUu77l?vLl$b<4w&Ad zTGvGThY4rujzT%ReS~xdpcMN(d*xq@<5E?*Lj|xTki$Tl86L(!n3;P3c{=KkEV0F=Df@Z3{z^JNyBZ?*`HkGzo2K;DDQOiyq# z5Gm@FJ+!**DgXqg?1zHp)S*=@Ow0(?N!V{0>bc8~e35XoAWVPShYyRPPzGoYzJ!Cq z?YV8&ZiU;oZau7XJo%g5H7y859!LkiO^nRUxw*NP=H@c;@@4%26=$27ik<|E>=tXh z+|m^JMjWJ|{b^5d6$kPjR1O2|1*x`qYKl6*s9HDdOwbX3gVvrj^e5;XIyK#m%XBsj zQ9%O0GvWbvtey*{kyRUd6z`gNrV5BxCOphSDG=l$Q$obkvJSX=KJwBN$cT!Hnwm0D zh4YYwf?nNr-M8sIm`ninU}AzjLCM(I7^vIWhcDLHi#>yoKL~z)oGC{(5f>#?k%zbD zsxgrUK*15x;vjCeb4juHD#eqSfZ$_tvKJKG;7F*c;roO26pKGlS(mLcGxKGY;}QK# zM&NB~2W|(U)3s8#4+U;kui3AT4$7>C?;*`FhPyExI{apf;Z=0Bw16AiO*nA_Z)5sb zy!4TwaHhq_fYghEKtsLD$yt(}on2IRHPqO#?#guag zIIscP{xadTwcP|K%gKpbsN^nCmq3d5jE$XUrFfC?Ic5vqqb*~%*TmSu($uK_y<4(4 zReEuOybOPL4gKTQf)WxDz0{#MI~90nCE&85D=Q18pYAKUCu=r-nC@THwBYc>UtUYA z!r`0gw>_*JyLdaOO@aRAr6@3H@B(RtQG1L`BV<;IDH^h^Wcn>c&q>6ft`D88j6-h$omEaf*~rWA2yhtsOUR+2KTcg`?2yW zXITn{sGCpNHEAOoIr7wVLnYG-3k&0{y{0dx8UmRz=b=APaxq^CGJ?#XD0`b(XV z!o$0HYdx~-p0=5AaQb8WUiPu<*zU952KY!l4BiAh0a_C*_J>NxVOHxk>hM=SE zxvecbbk70V%s$Y0oz%^P#c81s0YA=+U9)kq#G)Wfq?nl4D!T$xbVbrF2QCy86o3uV zc4g19K&nmZLa|@v9v@zCl=0P3e*T<@fWY|J*3(lkymiiwmopAq;A!C_vQVZixY=3MccEkLIw}QsglyxPjcHo0|UjHEEwP|iI+@YCh&&b z9~~EW6Uh~q3Fw7XB`E7KL?c9K-iL&^uROqPt2FoZ zRcDG0rxJv1#K7K>#mybuVjR^zzqq)#zV3Ojrd6jeY}!d0Bo4>hZGYu@jZL!x4L&Hmxlm8^CylF3 zA@o@7;bYufT;ftv*AN*h=0-+$v5^qHF#UZ*XhiT$I7oq15f}Fx2!({8#zgc^La{ox9b;@*NGNyX-@YfWmic<;OR9a#`+wH)D z6WO4o^>Q$h5f{smPw>Hm<&)@DV;@(A_4xI7W8~UoY`56npemlSns>pr0F??7hqTKJ z@aKh)uuE2QlohE0R)7TpqXG^7d+{_~Hbr@PkNr6f)t5PzS$WUx)pWkFUWJ30kGx=f zctR1T%SIB?!a|G{hy;rgBr6t!dZ?wwhlXOrQqj=RXpqBf6!z}!?qF6K!)~Mx`)$lc z*755X1{$0oR1h8FLaw6>Ea}n-Cul#t40VDjvUFzi90DMo#Wo^kKolm{4Z? z$n8VzQ}IzbBD{F}7p0V{bM$0Gyheu`1np{kb#gQvih%I zzkYWnHk)?|stXL@PR>C`S5;A&RB+=_xn1p7Gf#Jvj_fh|cTVD4=O_hS@h4CzP?(yw zfCH0f0IK_SwQi6&=W4muVpKOC3)zwNe5D%YKJ$)-nWjbPsx$>bP_q?xuG+Vc9|=Gx zP-27#({2KG22}-xE}#J*=)(b*XLu@R@ZGLiAyefhnvZqn3j(Pu(AA+R2D~rL=>IGW za;fNEW>FDGv8JxFaw^XfxeqyT2^z5)0dWDD1j)RsHKhct*M9kaTKeQZdFaTi7t3e$ z;^ODcuQ%m**dLRH0({I;$pJ>8E%g&{udpEC5d_y6U(z$CXj;*e&q|*B|3La!5ctz1 zLatm7AO+&)k~xp8awhz^;izVHOs&XkR;?>!^aHJ9c%>0V$6448_2A!{VXUX;t!rjY zO$}t1+0LpNVZ3k;;g3}BogZELhDcLzU}K@7=Btz6zOAaS&lW9vj(b|X0RA@`azM@>uqfBg|2el!MMVX$waRh3 zuRz1)dV-HAqilsx5KM=31b>itje+|m-(k&Lw(P}0jf8JRoD%Ydv7q2Xasf~}Jw5&H z8tJVx00T`W=Tqz4T}qGi-uFjl{Hv2%@NgWSCwtFx)dz-#Tvi9WOgKj> zUcNaMnE7%C`zP>o9-f|&mkjsrInjqe!~n10H6LJ{{g5+47e{k9Z!>(cF|+K^-VB%b zen1$!v*Q4UzfsFV>*1m%kYVj7*fb(|BMfB$ol8E`8S~8d`C@oSYl3UyLwIO>H&GQc{v@1x!CciIQTISGH0D&ou7X?)pcO>oP?{ zDPzL*AW|A^VYq?riaV3H4kF)_?@nfW_y%YF;J|*m1!iZ^Op;Slgsp}-<|m4N;g1B( zRQt!ve|^y!W<6eI4^BOhSqL!wuUsz8LAL?aLeSrK^;m6gbMRJFRD`lgjft`sqpfB= zZ0b7V{f7_oii#V@bvHRSlcqLV zn3zV3jC@m4rhM|Qi`mWu-le9xOC6c8aNj*SF;>uC1c`v1MI7?|muUif4R}?EiJkXX zo|6=3&5``J-Gcqi#lKSQs)DdW|z!B56Wh&iSx3<}K4j?YdU zyIE;NtB`uv+EZEKZ~=Tdow|e=Ax748i)K}n;?o-`NL9N`i=rHRrUMCgwC+HGetx#w zje=19D&OALCPHkVnUyv5-FwsS_Y7XHFvSasQUd*(NITGmNaC#@Vfw>fcPuBW=y~SL z^Uo_?Cs8&Nc2h@t`vb-YN6)ub7yG)=8PO4>oy5QEolXaV%Zc1LtFRtJAi!l`4Ga?j zg*9Z!r|P7R_@Ffs1YPVGdf(Nn^S(TtYU?J#lu%bEF1dy7Kz6sU&WQr zh-ic4)3Q+iL%eU)_f&*t(Y+PJUuNj=9ukvxmufMLpXz+DLV~AP?@BqCK z*AONqCXVG#?9|llmbN=Quytgk5N@| zapbeCaI(ix@PVcQFFH;-CM%HixVY44@HM^CV3?()FFEg|bQk95OXSixihgzhYpq;* zuR!_c%Pe0nDZbtb~Py z5r{s#lIF_lwvCvih3jKBgl|M3IA*t6!no;rtNLwUWwY&N%*>|;Nb&gd2oHK$^^ilb zE^!(SCSnE_4@9a=l?Nd$D=@|}bDe}ra=l{X>nYvfL`R24knr%rygYFA!8H&56wo=^ zeOG3-;~meCiDEOT|G5DSy}f% zCgC6jUo2DfFR-|l@K=Pad5vi4XRJ%Y(TI@nG-K^DO9GjjS6qmmAtYue=EfN%4SN)l z;{!!#aOYpxu$vGh@!Jz0+aJ>Y7h>6*;L1CUE3Mi>EPOCs%8=le64^+Rx=|iU! zeOafp04h@$;|)C}t?WKxwixiLr+oZqyj${qv9ErFF9|d%DKN7$v9V2lcPU)h!*mz< zs}IgC6>aU4(Fm78Ln*bV&=%H(=-@XMm(k5aE^lJ1S)@2`p z*+`oHqd4tDvQRJ%0~G*(5~!T{ZCCKH#K>lHk*#ahN&LYJbsGIDlBAbyZ1g8v8Qz30 z4+{RlNy`1v9q|Pxz6sD2APmAegO?IJHY>$lpdNq~qY8&xIE?3Ltjke&2wj396ii6R z7pI#fEShgwW=aF0t3?7H>-=C0`$jqlMMgo;b?UeF~UN2K0;Pb-fCN(wn)So|&JM--T(8xla!R-tQ9D)Fp!P*sD-b}cB zT!S}b~CdrA%Ne{bsW5RyLfQfql(@Ls^XG5eVM-DxIUYvbt)?8Ar^ne>UMV6-q6f++|t_PF$!6 z!FoFFPz}dnJ&)93ymHyh>N8NW5bkqwbKUna{9j+oa}}EfMT_^%PR0n7N?q&+>Dh zIywNHwAnu*P<*VD)xSPke&rO`&Uc2?r&D2-1IIf@O<7slnOH+#KTJ2a%T&|{ape2+ zkbHsqm%o^itn9BV^L!=(AH#2JYYR3YIIltLg8vT)xB+Ux@^X+21!O`{Gt9r|Xp)T^ za?qrys;Q~He3=RK{?EL~c9VC>@3IN$2{9Fql9X;ePq7qcCnkh062NYd4&I4uogYpG z1_Vsjd-AihZ-%hf1E<+dIgrX`&x?})LTKZs2dZ1HXVCT&#PUos@;(KZ!~6h@stA&0 zx#dur4ngpj5|UO`(Ho~=xxdx1($s%23sp1J94xG?fH9lFGu<|~9npZ#)6QeN0&W17 zrFK~y?7P(&WE*^B2$C?GxbC9-AX`QG&DX;?4e>%W zIUgqz$~4!tiSInTysqGBb`|gfxd_DGhYuf8ws$4peqKXNelT0qxZnWjcm4fq-4xfc z@r5?#K-_eYdFI5+DhOI3WkPmtfNmNP2cZ&yh}+!Kl4+qra>^=%rIb+A^~uDvv@{DN z%cE>B9Zc(16E#k7T&!Jmq~5=O55i#^M;h$^Har7|Afy==6CtLYB`*qq7S6bDxw%2& zX#m4OJi)@kx{DFd8mHIB=pl8At9@PXg`Q zKmp(${j#kt5l+Kt^|KLL`w6Sa5wvSD{c zj?Mnlr%yDItl9H}w}90lCnM`LEp0N)&dKSuD2Gso*l)Z(>7T(YvtYh3ZAmEmFisv% zkb1|lmg6~TYq)Q=1vCm|sN|%0od9NAhoU2kJJR*ZJ>0GCX>zB$@LDl8=o#944#Qu8t0KOos%0ZPGrEeYOCv#k=8`_T)B_w@SC*{2p+@PelXm5Qp_Eo+c?% zSM@PCLraexkMDe~c=yXW9hMUKT5(#4wIjl>$8JUF>EqmZ3U;Qu7)*?e_|CTr4s@$m&u}%OtHQCt0aM48eZP$O~ zb7*g*jf5sn;b&lLhjgU(QvH-yZPt!4zYBsCfId<|H%Dn{X=i8du#PzqqM%_OQIe_!Ln_zK# zass@bDgGney(f0vbIH_eJWFGut&hTNN6T*JH{5o?a6ur>)y{p4eGv%BbbAmM0Lv~< ze=^k3*A7C83Otml@;G+jj5~u;;Qsyl(94yPYvsm~l=b`9rYye#m`S%KY->?XxUZ zM#A|Y_D^||oxsJ^WAG|=d{QO`Sn>+@L&ZV zo)+HuAV`NSnZa`QG>e6fj3dV&g}sqW_m9u!z6RW!dtFeM*D$7$3s$oKa7{yU{Qq|fEcFD ztoNMtn}E7LyB7QW1lWcRJolelS>*}dr>zG6e3c-T)!+YD@P}3d5d68@K0QM!+zZ1Y z=Z2PAz}WdPTgaJ2L3d{$suWYtxnW8aAB1Qw|l13f4&KP^;Jrs|r zI4@t|f&X(j&jHd7u%v+Qasaf-8j%Vuyk4PSrmrEU#qj^}{yhe6VER9uoRQK3twJtI z?&qu*Ikvwrn?FDk1@_8>^BBqTX3Lfse}F@D{dYoaw-C92BzCnpKBGqMqHT|C%(4=ORty=kuZ`DV~zTS-a}m! z1ig6_Oon9o2y7<6;a4SPr=vmo$p`E6tfn|Q;kJ3>PtY>^O|KkOKz@3+qGg4ulkVF>7@0Ea2 z@W5W6*Vj2F+HZr=9*ns(@B*@%mzNicY#tsB8vHe!v0UWp=V*55kIugL5|R}}D&+@U z+$?qqrL<620m-Td8a*aFKFt15$4whg{g^GJPy*N!QhUV_bezZ z1R5;xDIN(YRuDm2$|Rl-jlURlm;osVre7%gA&-tt6=eG^b)9&qUkraO20@1@+6KrI zU|OFBU>85q0c&o}s9)BWINqq9Izf#dA@p$QfzK1Tf@b4~C4TYGvUslZoG+XXtv)Wq zo5wB7@G}}V<4SuhRXAYLJ34-z89W=&v0zcoM^7(KPmpwRzSpyosxZ&cl}n}-*?Egj zu3GoEaWN41(BqRHLy0BCt?|)-@M?P0f3WmsP=%r4_<{x*njB&85X*B!C|v|a&iKI8&Jr_UcaKAnfhtf1jkBG z>QmFwH_)F2f*Wivkl@6HuJKCVm0VG}a>*L+V`yO<ajGHjqQ`UBI1ja^g`h2zzhNg80B_i=!jYHE--?|rAMZnOJb^g zVeS(4U?#{I1Qn;eyxfJ?rQ(M5I^FuCJg?X|W6)FYLhJFt!C!}J&Nvo6 zvLtKbo|c!F4rSpu%2yN%3vt82sHIz7y^GtXuJOjs&xNL2VeyaP%O8Nm1LXyBOas~l zXX+cV%ae=8pUCE;-_k4r0|nADR4ytxge}t=XEYuvUVW23b#VF_17eJnp{HvddP#Aq z1Rx4Q_u3U={Pt~wbHj^xN6U}N*m7T8fC_ptg+YGrUlylTHZn?WH*sh(l;7f9Z-L4x zDmfYIzq5f915cvefyDSH51UEBSWqVKV>C?gAXOy?lu9K=uEqe{@7^6KiA0Z?L5zZv zh$af;SE|z4!I&*{|H;eCi_7JU4eM4u!KC&N&g&(!V*j__20Skix`4re5^r^L(?DJQ zfMaDQ%PZyQgSiWx9Ru&L4L1g9w*Q)NJ^^FPNN3yQH!}2LAV!1J4gMyka&}uXuZXPc zt0I_qs30Q~^wRyM_!J=hbJnd&0?GE`&o4d^WAcB3-q$A-&!nX>;A}B5oq$@?yZ;;U z*1U)7(QVImB3RTGCyn$n-&HAlsM4Wac>@3(s>d6?SM7XYbE@>TNy$o~X1P=Q*n9QM zkti^h1#6M{t}&yb@A<|^{8&{nZ9 z1w_6L~OSV2tQh#|~c-8$Z8JNdIRx#!wXr zh0pS)0jUB5@BFsmD+@4u0CO`lGkfEnYlZiFvt1-gra$K4p=Z0k1y2J;gJ1FQ3aN-)035D4!0 zcrwxe@9>lJ$i;MaYVSK{nAkxWIfaFr1GzE2zuWBb*`E>yiqrO*amB^OWlYK3`u3jQ zlb{vdU%?=y2;I!(`dZqy*%+1ei_e@v;?ZFlILTb@a!sU44^so4+j>S$!cD+CFaXpFVy$5aU({^3LH>FS73BHwJrcbE%ecV8hCnb!{YqD*iJha>hUWigBPPSk@1B z$H$Ey<%*U~FVW0{Ksh~a2swFi7rK+7Sp6iMnKTIoU>W+4O-8S0GM^9H?t5UM6l718oH7r0pnGw(X{OD_d4kjF>Yn7NSs-qhmvLnI|rioX4QmgUtZGVMsDvOaB9G5Oh4Ms$~Q5%KwMGH}R*kZQF*I2o*w- zu|Y`YN|K=vDpSZ5B_Wxo%yWj!R6@vD6cRFJCLy6Tn5PhuDVdV#-PU#A_xC>E`##_E z7kqxd`?_5h);iaD9>=lo`?hb}wnvWG)p@;dcWj4Ufc&viXSnKMO?34|M(Z zV*Qs-8K3rE|8wuqIWhiGfhcdV#SgZ6}Y40+s2O%9gkv(+DzL3 zTSJp~Cj0MTjZe+<=Rc5o+w^|35wLeG-f^i;AkMv_`tAHbdgY>6+>4&e#<7q?%BbI| z&O1|nwD;r3-9+gITdu8S7^Laooq=v!ZtidVZ*j9%$mkWKZY9%+A?N}pZj1F79){3` z^#xzjMzuz1S~vlLR5Q3|ii^^f z&DZ*`5o#Hh83>zERSV%^BfUd^;%8dk2v+O@<6Dk_%%iS*PeZ7}?L2|z2o@#&xAZ3y z%=mR}ltp6ADwmmzZR-rwh5Ujw7$&`L^i*37t=2r+;tEzYbSVixJ+#l8W~%k{?x6Ug z4{qZwA4outGO$@gp3bnd0i490lzqePBenZv&&Zq$t6!#dd^_-U#PE3h@%R@#BnT=H zBR${=Dg0U%=CP09tIIz=Z)b3`C}tgvT%>=S0v&z|`HUDERpt%oN${#HRSj=mEqb}d)NI^iBnIm`z7pa}Ar=_otjv5HA zP?b*?e$3Aqf#e+YTv+40iR3g3Iu}NIosv0~-A5p@l%#7U|5$`2k;RJc3t)FhuvmAs zBkxk?K5_PJ&AVrmg!Y#1UAjZQos*6zI8aU2c|U9!1JxSriDZy)(#{1bz-WYNfU*qOC>5(kT{@ zZPAz+*}9*UBxP9Z3&aflU}7!qLrPmxzr)K`JhegQ9YMt`jPzy04W)#a4a1f1d%jtU zdHr=n{|kO<&~ApvosVQoqkbySm%b`lU1x|;8_Q4AqYb_>kx&Qh3da_LF7oo%^EJ5O zfD*C6`C>6E^45`jDKZsav(^z|LcNdQ=HUyB08?V?hf&8u7zds~*TIwGKG%+1Il>rr zGmoqS{kNeC7a>SQCtXNoB-l*K4h>_t*6U}V$~EfGbqaOhxKQvDF|{c zPcCw5pe_mu3aaATI!pF{K6NaZ9e>Hn!XB@X*rEg{5bp_3jAET8Nz1#?uL{|kyjHjN zRQ%l7*ePt!Gr%Bs;f_!MReQjGh11{U+gT{ZM{LB^x$dzt9SPghxC<&K^oD53 z=C96X9WqzYj6h3;HwNm3ZXN5(#9LT>3NNs`op=M;?N1ep0i5;ps%*d5?Pf7?_wL=d z_93NjY#Vru{U6Gp)c_JrwML>$J53NHMK$el>HodOXm|hxtV1&|$NnJ%0xS7%9$3|3CW~AA%=zIcHAi9dGORo0$%m_wepnRFsNQ zbgCWtK*Ci_{_hQ3Bxbc+=6t`3KqT19ED@p?j!}}l`y0H(BfX~~1TJVWdokc!EG7JZ zTpLAbI?hw2xxH!Ow4dZRVm1akW&z$1d&7HI75=-s2mW!0+-Eb9`K5V+@BGpzh^?Q@ z99u0zQrT&+yAGN7-enA$IrJmABCCOfK)HK2nB${A8wK17IJ0@d0X;qaAQ?;Pzde5l z1jaD=G_}W$s+YqrD=S2!rS-hWU**t0_s3B>%>Ul5{obvs+rzoPvEw$l>2XV=mqIGfg+BHNkO}jYM;shvX#-lu#9xq5 z<0!`Z8d3>E2FHp{DMT-i!C-|IA%EU4AsjQrJF(^zMylk7gB0}prqP~rnmGsmb;{Yf z`q`M7y>?s3jEn0+n;PCI7XG&0|2|c(mx|Gn9g1!Sw2yX-VS(;4{c_yDw}Tkj4>X=( zTqd8P0gTlOHel7&)onCj{-IB#ym`Q(49FLoyn3!SAALwGAIdx>)lwaqEKPAgsuMfm zmPlt(R1A8Y8sv-9H9s8P(Yb~_4_V-xZjhx23y7jF%;Qe4VaZozij3=EaL=chbvUc& zG9;5ggzZmGje-cY(re|8%G7Zs+&v znEo=OYv@ZDFL7)cN zL1LR|1-VC>8qed$kJ+q;i#iWtpRUIA_SLRally;KJO?|R`}aSgH=jWI@h2h=D;0UM z-bVS-+#kd%#<0M(s#R1*U7eliJoi<8nQuJ?9#wro3FE=VL?epDGCfk6=$!^1QD})t zNPxk5KJtTaKD+7Fv`8RGPyy90eemu0dQKl1I2UR)iMXT8k%;2Ljc@-q4Dt)gzb<=>i+6w z*hfW9y8I@sX9DWFYWqhBaeOo-D3VB_Z zc^=k0ym2w8l$>uv9hE=PL<|LdXNz6H*-!x4kfuRnWxr*BCNfBAxvLXT<`@x6H|Kn| zEN%A^p(TVW3w#HV0*!1_e5Z5W%V}%V%Ogh`fB{rP=J;u^kYZ!=vP!FxL{^UtW$JGi0=XHK)AnkHzIPERDYgg22haq9KrvKZH87En z!v6Vlo(iEQtXW#H{a3i|fX`jwJTqPW`;}~3aA>Fl976J_(2Iuxb;a6(eIkmUnb~|L zL;wwWq;>dOj+qyC{lFQ3m#2)B4kJaC%f`TzdWe?tt;`(~mjCWwjmUnWw1{9jwoOSa z9^F5z_+1oA8J|Jv-9~`KUexhNJM5H_B;8(_GfP12=mz$a^FydrJCUzJr35l!pSr@i zbM;`0!5Hlcx`@5%jo^yzfpCbPK3{-+nrjrLGk6MTk-iZ9{5mabuQg7^@y#CJxoZ~% z6_YSroQMQ0V22kldg~1m;1pUHmy#vRV7q|&AM7`x!y48q7!&&huhyE9{Bo^?lpl7ZPoF)z zUDn_DN-uX1XD=YwX}%}<`M*H!!zbE=k;ydAXwBA!ExJ2sXM+-Ed!CsmbQ?hw<_<8| zlz7i{)_Wc688?LZ+jvA%QvMa;AR|B_o%ib1gEmu$5zfoY->{sP^4VOExXHJqU}+*U zb4=pBnR|F!tSP6mFWhs`a)Qa&$z@}3&he}+a>6=kKYC5jxW6>0O1O8geeod_-OE77 ziB?~7a&PgbtRUz$!%;R}t@Vn<-t){@tx{N5@WumI+Dw>fnhSbo^{lT#)qh4{XJu*0 zU#m{;cMxk)nfpul*r5GGdsvF#bK9W#2Jx=*bcCstF*P&7KUO0a!%(s1zR+}CKCter z$3)F*brkdyOCE%0HtBL?snBvulM-Mus-~~MHr;cA;0yB-RaKhMbkwN$IH=RX1nJpt zL*(H~y8V8@rS!+s`^)jrbOahnB>V_`XaQB3fE&3-zxQ#oNy2=9UtD-g>1QmcaWK(T4 zks!({Dk9c^+~wfnV$1I``FT4bAD9Dz-a=%7BSGq>C-OkkEKd;fpbF5aY7!PD0lREJtZ60xvL9pp;+F_5JUUo&Whz z`k#O5mupI~>ftc_CNQsNvmB=v!S~R)yPm(Bk03qB(&T5%QE^HZv)@J##@cUmz(Pgd z5A5U@9rYbuDynu`BVic3!g|km4Ba98KLQ4hyV$rmY?0T+OkH z3}zFYt%;OtT1t>#S>a(|@%6c{uP@w2(DwmIy@}igjfh?Cx@Gu#@T<@)LaG2dRh?me zkJbGPMM#}+RyDpqqLN_JmF&A7fg2-BMZ+%sFgclqpo=Y4Vh=0ItDc_gsMyfgkizK( z3y-<#H9=8Pc=O#s`8f#UF${@`{+T!gv2mfhA~KQnnAb^toUPd@6L(HF>$Bn{{7P(a=DOv@$w6iVtI;qobp#i7il!djmfjoCY}I zKy1Vw{7zJQAiIF6EI8Rcj-{coE-%gWPKSS>O(Y=|32+q|cn#({{hA$qiqm9iWo7G| z|86{}#}6NV|M^q=_INAY&vElUz)(>0Hn>SAfLD0+g5i(y6?}vCmLC^v$G|X@l}Wa!v(OS6Y3urKs`x(snXKYDy~QV zVPQ&No|Pexz}2(aOhjBlVy?RY((o^niy%~7N)XAKe}rXh10`77PZ@Pb7K6?tm_%ZZp=XR1L1jU<>3MMml*#xSe#vzG2 zKtn^~4siyEl}Z;b@CXS}5YS=j{VlHFXeRQ0%)zqMh?-&NcA1n^F8UDSm7-F*N)ry3 zO^D@mRKsinV2TU?;2=io5H`Ki)^sDtnhX=R*H6!Qsgq--Ec z@$i`8eAeA#RG=#uhviUe*(E3`S(%-^lVAs`Ut>MXUZeA(mM*}o=c-ncOG{md8C%K? zJYm_ie4VM-gK;hFzXsnPk{&{*4et)dHE3_+tUl*p0)ecgZfg|b5UB-EDAbndwI-akKW+@bGY4E{@!jhYyc}w~A}n z4=1rtW@r(3B0a$UQ)iG@Uqq6I!~;a-x}Uf=itw?4+c*OwBbpKL_j=Yf1#gcJNfL-Y z;yOBA$Z?6#dr#p6g2Poo7nC$emI08P6zCE$!51&4(e(#E3vRVYQ_;0af%f0B^96C~ zAWcJTBGj53Hy6=^4+hG}6`Y#-`hdE@aV)av$c9)8DQK>GT3FaFfR<)%E5f9zRsLN3ED7pHAcPQ4bdmWoP*|J5REUR2(w*F4%_u)z6X`whn8M=OIKb4Z#X zj{+qbZjfU@%Jt58kGu`RV=b(_MydX;xfu=1($bx{6-lblM4wsx9*sPw>FrxjFnJK9 zu0d#wOpoaBixBpNCv+n?-3&7eAQ9U={TA@I>SHlFK0bz>F!zF^@5eK1RPvv&V{8Y? z_sV|)B@5B1f|Ki{c)&mD> zV1LR0l!w~><0^_LyKf8qFv6Hr1z&P3c+$POTK0#==0z#rXA zzh@8ZDWdP)gT}3TYjcgr&PY6!cjX)~E<{!VtYkq!a)OA&i>x=MGkX16e+YoF>JCW< z1q8gfUbEk0FOe9C$JBtmf4AK>l>HZn8Feh@UbYc7-i#PH#rK_2s7Y=(X7Y|!G;wBZ ztRq+d)V_Uj@7`UQUU>KJ9fXC*qN@s$(d65hy1w-jEHW%-Q)GmA-Nr4wA~SX5U}|b= zSa?{OFjIqU6u-0rr+soVGLcv)41?(oe|K;*GU6Z{Jap&=_AW{uId#npjgy@vj+39U zM`2>v4neSDhNu+(!Go#9&oa_yYkKa7#}|@9G1{V;xhDeOi8Pd7zq-qvaw7TyySFo> zZcFmbHJ)qky?kBEXS_8Yai8wkCB|*X{Lw9|Dkb4^sgGo1L2zfkuBd2A2|IsVtmw?GTH0H^@%zjkhGu?E8SOs&(Y0ve6K?U$P|evBC-$6ovnADI6D9wD1dc!<(*1%sIY5zAA`$%d;L46Tt5)d{ZnFo%t z?jX2M&O4wHij!AVml1gHf8K*yfTqvbn~}l~RAz{zg?V^_4-dbd?El^EcED><9mZ(b z?YT)~I|OK0N&f|uaxV(MClT`NG+Yht>+0GkbC1$ZPEexse5p2vL_dI71Mu#g#{1KA_=DnrQwmj zt$Ui%&lB`NyrlQv#jLTDHPrv*xBt@xx~kBR@}TOK6f@UE`S;L<5ga6Fc;oE5VzBUk zT?fof1c@oMBs&6Suk!k^|Cg61V$=qQZBCskz2IiDU*Zx+G=N@y z-PEorY`_Smu{_)l)^loMqcCz7GN=UY`zZD#;z;S^v%r&an?h4bSn}vjF{WxUXdNxo zeDWzED98uyd}lWnqhsFi;0>tg=2WDzQ_LPeI@-d2^eFNWj{ENhLpKBd)pelAqbn^e zy#4U_PQtJ7C1rj_QsNa5OI`>FM>Kp+03|jKHiL$xe49nx+4L>@D6c5lG|=7r>kMtx zw|1jf_V15DCgVq4XC=kja+jPy#U_@I6CWDCx{OWJe^$?@5yMYch6c85z7^id@z*(E zfre!~7y?!D8)wnPkYiK00qBWPJNi3+k~dL zI<6!CqWp(5U0-7t-9Zyyn}`OY`7e_5&d$zW1r0uIHv;t4#c`!V&b+atWy;*71n37Q z9m}C3t4uH>!{WRfj9;ql={a5QXDjFS$?u)7H6ngaMW?W2`lBHF0Pb zM?BTuCIF8*7z=7Bq){fcwKiJ)fGV0Gt)M_|pdx}FvYwL*BZbsy!=tUH-Z?P-LgudQ zUCXfxTj-Q*jS2I@Yo3cMu+DHD+BllMv}r)hQ$LxpKV+w6Xtr;zdxK*R{`ZL|w-N6q zl??ISl2k!uJEEN|ZSrpZ=jeZ4E29BEVD7N{#^Qsn5SI6eJJys^8}1Ctto+hhn-}hbj zJZ9Z1g6akbUXO-;3myn*%iRw*CoGB*6ygrU#&;!1h`4}9^w?m=2)Rvpx%;CQVOo+0 z52D$IFG=Ih{TAjueLB9*5Tp{=Fc8aYKtH@kj~1kr$NDl-+=Xw?3o&iQm@_vd&8ga< zGehDbR~^rFJ+cCrsQ}AN=|6B#y7MAVaSe|!4pU#7x)G#P#7QUsvGPe#N=6cjzvXVOx9Kgb-(ayNa~5Ix`UnezxllOC2KKOM3wtD*+n#F{^jozhM%fysHPe zwuMEumloY|=~6>`E70ZTKM~9Nv#@2SK5%9oUDr>Oj-m{Iv~gSSVC>A>N#0ramW^1| zylBM?ghwRkuPMxk_d!;l4OD{q`X@4-?J|RP-F-hMCQiYkWi6Y-#Miq|7|KNAG==vM zj{=&DIDptzn3*|FYpgV9mt!l5_i8Ad7GvrA?f7c9$Dg1+e>5^V=RMg+S_cS zu&bfmZSfD|<1uK#`n(o4i*2D8Y-Bb!)-b=IYq1`kTyhEuiFY}&X7(Z@v$KQe`-4o= z2OKLxqy&uScFcQy3nU(u<-|6Pi_2X0Y$(U`<9ax{cDkQ1k^TY9Ms8e?+j@#15Sdwog7A{p0W6Fp7xos}y_95ss*eQVR=82{)4X6i9OX9r zn}b6T>ihage_^(5i;z$vR}?Sb(ME&LQa5Ia{NB2!m;qWNPv}pgO_hm4T1kmEcXZW! zjm9N^1$M>9N_W+4Co51agBf1$=2MmRoVnsrgKH%-aie!#&zCT%?x4 zF?UX{R=`mF7_J^~^>cPM(ikR9a6xaB`P)9rZsYH#my0HE5wTQ&v>{j9I0Yu%>c7=_ zqzA}f6H%L}17}4mO!r3x{nL9gx?)F8@P$lUHCb{M6SHh3d=)dW8cxy7pnspM$xBCu z$B8m7d|M32$Bu|&JwAsxzmrVV9ojtf@=jk?+GUdBZ%r=_REZdX#PRW?kbr;+cl2oq ziNr1*!R<+^AaKoelyhy*(17QleO~;>iIqJ|ee5L_kFHMh{YFz14Y!vHPe5&bEpkJI zg^tYr-WkULno?&Hf^Yisft!!i+s#Bk(q0Qtc;YrydC^yMsJi;qtKUc;kg%P|Gfld= zqZ(Ly?BghPh@C|Wb8iwSrR?4vh2b1}gR@h_jCD|R>+|Z(W2M9#D z%3=>h4;j4vlWW0l;Pqc3Jva^6bT$hn*!cL&@=egxEQFho6QDL-b_N<lOXwWv zAT#Nq_ejc5CN^#(a*YFJ;{_5`6q%JGm=x`u zH*DbP>4|H8_Uu_+o-NpKTv4^e!*A|E!<~Z{4p<*J;=hN9)S%kae7Mlh4N<@JGu|7{ zwRhBo(=Fbbf0(3y_5h=1(5&g3??DP;3+ECz>+gPl#!gf`g<47i-Ru~}B}_MQb9IGl zBL=P&TB5>;aEs%oAxJ;JZsuaYlKoPP`>jb_+X8q%`e(x-`M;c{nR8EzpVZQSWDtXelao9xr!Tv3vFA_YS!4O@+e75UeKsqL_f{iIuB79R zW~QUNaXh{;CgOHi#mC5mnbO)hbkyMO0zSj((=h-+R!V4z&-mIaF~YdRyT9$qA|FIl zO-(b$%;l}Yd?{aIq=4ZQtsFn2&;2MCOCu{qN840W+ zpvFJx5S_xgo1c|3)>msd=S^=Q(nwTlaC1{g`sevS!j1a1vNCtM>-=PcNs=neV{aH4jex^r#hFRsjEKxxgeB-l(Bi|TR@GMS9T$jZ2^(hHrNf@c z1L%WY9n1v7Li;!ezn3?kJj=c1GCGHi2B-Z>d1GR596W|ZS(@Yi-{7@-Y7e&dZ9Y>? zdNuDm9Ybg{6oyUQC$sR75c}^GRh{}EUJJ&B9Sp`ux!bqb^p=b7X4?;qdk??auH-=4 zde@GT=x7-!b+X&=?acHv2%d@=XiXt2^MA4X-c^&F1?{swXK}%U5hO@*7iO>E1+ORA zKIB)FQV@u4!I%LXC9OQNA3(y`AZ-N129g*PCSVCT6@`dWKgnOE0b>*F*lCC5wKpsB zPvnd=J%4|{7^8SoD=F`sSeP+}$)Ug}En%^*C@+VfCC*Bu@@VJA76od4*I4wZp95?e zd?oBBF*B<4Am;6TeG+NBX?HJ=p6?1_QoWt4*xMPH{(H(bCHJ zGI^JV{qM0QV>Lo7ZyfNT8v8%1al~v|_aiYmL;Ka-%q(AF?Ar{n+(1ywo0`kKZ7gVd z#(j};NX18QJci(fd5O~l-(cg5*$yTr=&UnTh?n2d0hs~}?U0lveUj4HAL1ElBKd(w zP|*AL_YeVNkplimz_jT+;%4m!4Rh#D1=Zhg85AS&Pbe(6V~9qxn73$q{4=PE}Nrx(%tsETGV_Y{5PBpiBk$$whg}HBhkz0zI?Y9<^E!6gC*)4uA|Jb6n=iRx^CNZ9`pUlj|GqA;by3_P*8#_wHRWTz z^>;h-g6+gg&!|8_JA#w+RToPhIE)BPE&5i}SW5s5;n(>T)~xVYhaF^zx*@rBF}S-5 zj*dlg9L;TQ9cCgaVm@y=3*B1WfMz0KfjA8bZ=p&!cWi?_uxSX~VW(qQPge?-vKo)$ z?b~0bmUv@$zn(2>84z4zpISLHWT{|{9lW0>wA#B%x0xt4L;48pT0iy`0s%e{ZpDcc zmpSe~TJcV2#vLUPP^z2BI1ZIpchqOI-XjOs+K2oHwl9W{MVMnmq%* z@)uBgsX;?FwZ30ZJrB7*P|nGbkq11=BKMUy`8J5t#V~RJMK4q*nwkL>UrlM>ji3k2 zLLi(-y<>g8A4L)Bm1kejqK4ILJ3)KjOM~wv@4~ZjVx0yB7hgDjdnVu|1W(hRpw4_T ztWc=J3v?rvz^nq9swx)r`PcoIQIW!mT3p=y`9iy$yqNi!$)CkJI7UE3$?NL+3d0Sk z4|C_+7AIZ~cDK08wrhtyIf`dzQCxna;`XhJ6R(3}1#~?BfVJcZ4Gj%o74*fwFW}-=uRE@9yRFx|Q)Ywu+XY5QRn2jJK;}gx z47+v*i6_AV3awvLicw1~#Fl~rGW%LLw z<>v24uwg1#H9i$!OZXy#5wAh@sH5EAh=VWS$I8m#fh-0SAH&X!X<41kDp-1S!(Ifk zn!q;GsF!R44)pSyjBZGFu|`}yfNk?eH|NitbGUXlHl2w8!25K9aZz@q(Wim_!)8CFzF6f=EhhhcmTfm zHZlT^5Hlm=I2t7=BjA{cvlR)u{-99{iiZqQZV6Iu34q|MU zPr0|4Q&I4<(9(jj;}7QsV6f|utn(uc0>i89nT#TxK=JAAP^zP z0egNwzim^ib8-b@05-s5|3a`{qjP~;UR|9oNfj13vT5y=idasEz~n|iBoH9(MC${% zNpv*2Z}J@Xb8@Dpl6}V*gPKj=4EmshNA8chM2!J94j4I$C5kPE0|}1KfGPoaqa5CX z_7^e4@?sbh)?LIO*VPxidSicsrhJa9uN19JqJIHaCuUoreP+OI)>3)gN)oshwhPqA z;3EfbXtC`9W>4^iyi%FlxScu7K#(P4Ki3s*x=2Q#8%Kb6%5gZ&^mf$Olb$bDu`lAM z!6s5uL!C$ar($z0y{p0N0WBJDYze7D3f@n{09}kr(Wg3F|ygjef>|t zN$yXTMu#wC7F~R>-@*7^GjGprG7|xZSv1jla1t)R&`n54U{khUOg#LPv%F?>^UhI< zGi|FP2M>M&sse3>c6rX;oq;e59c}&mA}0s7o9L9lSOjLvpC)kFTVa2Fkx=2dU%9F9 zt`}oCzK@R^HGjgOlsdzrVsic$iT3FO1xi7uo)i_~Ra57i;EQud8q_}wFm(;PeoL4M z$G=-?ly*+r`}g?xRF>>GjyoKz`K{@JV@&vza3;`s+u1g=&_z?YX zN82ZMluSl-LQru#ls#@RDjpwXpd{OHU+Oz`o{fgayJXcn|#bR=Z%KKRBoLT5D^4GfEN(T zYC|auCo21Gjw&kDNvgLBbkUE)C{>S1urcAh0h)*Al`ELFNJ|D!QaotEc@9a*b>!4M zq1aw0RGZeRImpdCoLmogf8k@afPHd=0`LQq7WTXBolI;s^r<-C@;}X@^l0-i zw_)W6cb?#Rrsm~K8hLd@Fm!Y;|Gj){x6udgWWm7% z0z_~qR~h6t*4HtNCs0;~RxYm#a}ROYnjR_XAX-NwsGm>{=N#QWw9BJ%#q4759WK(8 z^yY;ESA9GK+6aUaFx z6$686!+U3ywg1BfsPluw4(^92v(bR0Hl7h1advUZ!>+nSC=%-ptuu68{S6OqjJ3zi z4)6|6or-HPO5@?B51}GIXKjtS-sp22J?on%y;m?;sqjmmi9OoKh*?}$+`1-a?Q2FG z!bfJ8298KR&Yvt);SK!xc6iSL&tEcndc8niLCySX;RJdf)KFj~Xy+^~EYyEKX1}97 zw^U3RdJ-q%RH>#YpU%-jVW)1k@yDw_);HdHx5S8y#Sl=WqpiEolQp0GHD^7D!Qx4< z27xZL!Dv3<*q#qZlIP<4ZX0(iWX7&iZ?H-}6m7ESrvc}h*eMtT^nsCOR|O~LUEybf zBo9Fy(PepM$v|i3x>x_|pPS2LuEw4_$6~5p_F33Vnn+6Sx0mdY(o*_Np5@nhMqi>;%xd;- z8SXgBXJzifq9TUAhsTA5^TE62<7*ijf)78kF#d}*$FzoI|EW2+l^+R%N(62r)qyjT z)gqg_OV7Th60~cSZ%M%ka`QORrh=tO+hZl6@b4Y!{C*UMzgE}bvck`+1x5FCPqYcg z(%Pk({EzjIz2;}fw}uZgLXv>E2hJg=diWGD{_W}_x-0?WLCy6CjRiESxfpgr<9Ojh zFkq$pDPkH~sKXY3dw~rjqh6Gr?)Py54P2g>_aG#roxdmXUb^AYWsljKANB=)SB9+| z_vf-z1(*Ejp#3ZwR!PuUD352l~*^CHNw_c5#8- z7wil^PQVMyDu|ss`Y3$}f-em2zC_RUnMI>SMsz~5gK-I-o9n0Jm1#05JWhLDc|Po9 zFnJBGlNcgvX|gaJHb$(E;PsT6Doxbs@g9$#{{HW|bvGc1~~UJN~!A(SiR++bPUi}He~8s^>$mFGxwM0}69y3<9rEx`V>{kD*U*Z>zaL~d*nozh@1m1WL_8tvt# zLoq68O>$IdEL3-a$Rpu`t95=4C;WJE&=AM2Foum`=oS-Gu!pL?z5vay8bT6=7~wce ze|~ayO2@qp^Ab~2zv0vCD&(ybGlEP8hFH@XSDt-AcG?)o)#onG zXg|$0IPDi9ah2mfGzC^7M;XJeUA>BSjN|+wX+P;@zD!l#gtwe0K-ZyL1fLzob8e7% zLzeDrEs{g4A2KY&3co!Nstqcftugbc_6}*F>@bKg_)q&W1R^B`F(7e&xXOx)lM^xu zI9jnJ`L^}hY+NDuqA(>Eyeubo7WzNev?w&rkNc%^$$U>9nHdBt=n2}3u^}^r%WKiA zon$i56@V8fc^x$UK7dT66gykw4$u6*mh_WSevzgoVaMOw(-+EOFA%bDQ4m(biF%g;e<6sa$k!}Hy9_a4};RPy_bR2faj?r1O!{ru+gYYPfN7e-f7gl9U)1}8^_TWC}zF3bysR%ZJl7(v$Ow2)y z_1RoV0kEZeM3>a+OmGaZEj4CIS$xtYlh1spdQs^TBHJM3Biji~ryWf6^jJhgrci)o ztz2HYHM^BQ9aeX@F05UiCltk1uX;%jggIqbZEW|Bx61}~{qr^2_gEvV5>sx0_AG}@ zUsv+pQ2edsBJbks7xE^g)Mf|{miw2*Y9MeV)S^}c`6iU>{Hb#}9Y!ap^swbFLiBic zbAI=<`drG^`jO=_xMf|s)WkSa>wNg_e%st2HJ-1Nj@DAd>w;@d(C(osk9{5*HjP9m z!FI@|BwJ_OB9sOoi($+Q^n;*)n>%L+7`yqFkb~ugU&E3oC*bMz=+h)5CvRd$elwVf zo<8D9T@&=Y3B3rwR%;NDAErER@s)?aW(gIAd({~}z9}>VQ8h=l@4jaWM-|+xkxk#< zXsiNem7WmBsK^?09Xz3!AR#6pu^Xfy5cA-$Yhq$3dRP$p z6u712BjKFiQ>>SSSt>D#XwZR6_sfFvdj90KJ>1G_IUTm1%ddo`YmQ(06?DK`3>1Hq z_gk>iwCZ`)@R-6O(Qw$|iERMULZb-l3SM?+H~)|hy=OHab9jP?O>{+SO!+41UY%@k zw#3E?`AAe#UyY9s*m~M8bk#L9-i$fGHUVy>Fl8nXx^Jr)dnfw8@b1wbEErQLIExWd zyJuGuG_t-*pR`))WqD35nR#~k>EC9R;kLkpK=&jJ?z*A$3|WGCaWBr~i(svaJ%)0& zYO}XKU+N z-P}(5*Fm8RW{<5z)2yz!v7(EM3xhoFxlNU|m0n^0bJ?_7M%g<1i#Aq+cMvh=yaumL z4XS(<004FJWN~GsCmP^r4^>xh;yEC+Kw1D?qpwdtOVi!eRpe#l%)y-NB{44MWuLG| ztTYZ7Tt1P|shp;^{AnfRM|xSqZLiI#!yZS^sWbGk?Da^DOGqHv`hf{rV%bI0*R?G% zxbFq8=bVbwfoSl{DtWE$XbYbDD-Ku@l4tE+zwA}H7?4_ZpD=sY-AsQ!g=7m8Oxe$R z_fMPC`CM#yRjfAytq+k19v25{=*z|dK7JC96;u|O%3AaN8P{aWvr7Vr?~U>wCnaI| z&?1zQSE!QFtiZ@+GZDDEcYh^4uef0zPKR%{#DcfYt^Rw=iHPM9!kM?v+~jwc4(cc8lT&{ZQRR zQ+9)^SK-r!x1}XJicyxf)-|X6(w{?&)YdPMc9AiMF)lXN-NOSf)6v-*#|!Bb5Hs(D7NbSdz9DqPE0fCscn_CgWp5v;vbZfcvyHSaN84%vlM0m#=y8M6>06e4qHA)M_rQncUAd)St#%^LX_gD^h{PYn+e@JwkO zs~0N*2m)@W$KiQ}9uU#tI)ugnmLHyzimEE4YiN?8alk05gPGKOfm>d5-0_S2WJ=e& z>U%pdHK{2IsCh*^DGcHK{QWVN3(6t9CB6rk1tb$yM3EQ@zH{ebTYlMm+OEnVDqa<@ zi02MIxkWUlQ2#CeJM#IS8js{pC(Vmm-WjXBd7PL_2Lk9>*I%hmpMJ-wPDMpki6j}4 zt2hCD$XYiYG=nM{Zk%W>!E?ec#n2Ztf$F)BjG{ro()Zf}B@nb?LvxI8)(O^A5nO-_p_&yA9-v+0oH2uiHWoTm8C0F7~qFJ|-~JQd%v4^K*KA zf%6UH)$rh8xXr){&)|-M)K@7nst;0jyW{SLewRODe5qPzCAp^ix>A=>WD_5)e_IKEG`0)+SDB=W%20!a9$C;xP>|ht2*~nCN#yHwpHu z!F6l)!i#>KFk zVwu-Y%DxlGdvOcRHI}Z>SLAsuXTV?v;Eky|Nva=PS}x~wG$peqvYzNtnkmk9-1{To z9MN|5yl0e}`cK!MKL*#A(MDd!fji#jF;xi6#@ut)_v!U5zP!ygfk<$84dVq~f9+%e zr5WoGY!nQo|A$MD_LtD@x$Sj<`{O(q`UMC>!|jmfdQShtr#x^YNEviIpl4Xg;C+jB0He*}Y#~;9_AQzS0zJ&@xn5o-png8>PAl2P9JZmHF z^GhW0qz`nU!~2dTLiroV_d2{3i_&wy7Hv0Df-~#9B zHu5v54N(8Qp-TL0O8a$0x6o7tjTV%TNTJ$B=dqYhI8W;yHoD4Kz+ZpbIL<=ndm-kA zh6r4+Zd_K)yDj|ZV251N%SVw_TLolboSIg!p$Uj)9{5r807>ts^8Br6nrkv^`lS_o zh+*eBSy@-fWW8L(KMY#sUzv8=T#5=iyZQX*;aoWnc)+9CnQKth+t!wQthiCuCfN2} z)^EX3!AC2of~%2oLw|kc!4#Zo>Jsr`?ozgXpq(>=GL ze@Ya!ep3>}Bx2USt?b42BzRSG{b-Vm zi-UuriVC=kU>3ssq??f}HmX2yN(iboGj~X`iuI~Y+Y(SNN3os${6PE$v<8**fd%ma zh|houirQZl@q%1^ta$zxw-NO}CQ0nuzyHIpThqI{hvfTjx-Qw-+S&>U z<<-{KqSw##(0zn%T|?Hhm*g7mQ`u3qo{^4$EYNV<8# z3Vn#gwnptGh7dV&L*4!Pl6rfpGEE|CqRwL2<;$1ezrK-_yg)W#3BXG5_TdxV%p}w# zbkmfCS_nFT{9{pAr3IaP>#?W)m_VdXa%#x&=zGD8O~ez`(-Y|k2A4wvQuyi>ay94e zGoK$6FYFv&LWG@1mi!)4d%xuu<=lJ zhg?4LH+CsEqbqwzN9lG(C({)WiQ(t0or9|Mx_FJ`w2y0%r=#q-ZD9?-F+(wp-ED!} zzYb1!#)Ul&Wo&WHFOQ}_Bv(AAb5rNmtqI_v?@eCepLzxYts}K@XGs^KTqgR+;;(Uo zpo1n7)14o791mbsz+L@la-5p^I2dY%Zp9~yPYwlnf9(2wBu2%eye&#G0~$Kl;;D&= zpsrmRHMny65DZCJbSFyeOHW&$pb`yNX1?NHNu<+a(bnpu0UL34*3RjGbdB7wo2{Uy z z9kdD{FrFXPWs3|yUL0j1 zKm3J&&DPpJ9XB3)yPkOd}fj?b@$xrZA!T) z`(X$7!xIkqA{>4cPxTt!a48mt25L>dx5>)AramQqGm+8hx{iQlj)vvT`;6I_R|h(c zoI2cpTI|!mDiqTxkgLbAo6VSPyDa_fyKLhd_e!k{z!>ZVgxOvQ)vlce6?dZ@jrw_%STP_*3w&)VUnKo^aH@$octWT$d)3qwiZ6Enl<^?>jJcTVl z^Vn~Us({TFST%41(MU|9lZ}_lR8p)O+sQ~#EWBd#(7@1pZ;*D*28`BG9I;}WZfW0_ z&!1t;|B)Qy3eMbrs^VVrbK|J}<=(!&_{aAv<{xR=zpnZQS)J7V6}cGOj;WtrPu8Cw z^`CrWQ;Oj+908A2pQ%1un4d43C||kq_SR(NVs|2(yGx;(0DgGlH`e*Z2W%KfD>dRM z9>ftXyO*Na&~1C0UD%VtLS)HcLZB$XJS0!p38A&mE4x%Zvaw#VnpJwa;mpYLva{bo z9&?_b?E5K*e^L14$ta+*`~s^|#vzSI47ZDKXQ(sybPNlQ7U<=Y$-saG2(J&@s{N6i z+YzP_7dS0!ZC#pryaH0h>`{uQ=cr2rqJ{Q2+d}nH?Y|yUy#S41Tm|%)@yf9Y2|w+` zVBX0l*cAInfZy+AiQ9{}=SSw|hG2$6ODimwpa!3C8vM)vRc+ZzCm|M zLFCjagnRGcM?TAz~jJ+v5ihWSd`BkuH2}3%XMm!TY+_E{ych+hM!@u$vM?o7&&>>!<5G zvpb#t`D$!x8g-hP9=~5I)@pj<;sca>oY#+Vu|HrgXQn)aM6|oxOf*qnUq8kDXBGw} zqvatZIB@*%@tfAK0=;*|j6T$;5P8Q~RCL?C;-%EVgS97$LmEa<%OKqhH}7(+BE|ku zyLXgP@al6sQd>K5V1?hF&1@gnhllC3N7oA1iwCybstaD_9=N(8{q+!iO9RG!GOd0{ zmi!G8?-nmGsR3S%K{X3Z_@Dj*Eym{yLx3W2H*1D21^Mmb#tR z!!vIZ-VBq}_XUgIk^XkYnmCv@ry~J=i0=m6nDfP+2l|m!0r*3pO;dfJEbPfVCLL+w z>~K?;Yp6M;Q&gbo;X?SIc!|Ffn0WYk(h&j1oI>Wh>-QWn7`Vl?RMGm}F{ES8W!`qb z_M4D96UHp%UzC&4?CR3V&;O(=8cuplVuI* z4FKdLLqae=X<^!G%w{;{(Ic@(UTz1TV^i&vG>ca*!3ky{ijFy8Qsp(!6;7ZKeZVQe z7zWZ{ZE9KBYRHjfre7^V?uLckou_hra+Y&1$EkXmtn++exh2B#PX@>sB4=|h`jE!j z8A=)&Q@`*{nr^ES?VOpGM$^JF`92MioD23^bQZgdb@{>7%grTri}1=>&yw-6QM9y~ z?$&^nG%6U(Qr~&sDA{y>m?5T50Vn(LLDhvBd~JY`Zle*W;3?t;H{)L8?z}to!*+@C zRG+@A+^MOl8L)ZgIgW{vLVN~6K@z`Sq zESTDy*rt8uO8lco?3!yS-L*)7#^prW*^`@--`_;%)EFckNU_E8arRJC zqt)7d=~sW8L7uIdncWJFQ`uGekbsX95=kF`53Y#9PwLNjqAk?yRWFCWZPR5N5CcA9 zX*oqSFj#ir%@Ny#snumMeIeog5ot=?Ynq@hot>FoQ^2JLZXNu)npWVL!%7+DyceeC z=PA8%{={Bu@jf}yk(rr^rgp`Es++R{CA?vgjH52}I|$yzUM40ntaEG@%&dB8eFkF; zpy)P?OQNHh7_T=p-AfU(!yShDoi}_*CFbJdbntJsJuFkp(&nl#g~|kfm3+{{lxqIWP6Vi~**wPyR@H zav}IZh{^saASP@y6k{79<(OfGoM--9eGaBTAu7hlpDSNEw9VaKZI}17Fv zQ$cZ$mE`{2=Xd|5dpUy7eDLD$>#;tE=x_e{UwvJ9JXQPpCi9S?YtAjs*$EjpqR6n1 zoD`yjj2mIg9Ho>Yl*&*Uq9mj^9dnT}`+`Iqev-Vnh zt@rz`_xrrh`#df&BFYV#=QZK+dSWh<@McrC=%yCI07)#lN4A3rVbB8}>p z9BSs9$4Alw+7&0_b+8Ik($eIw?xCOj zgA`M)M%~aySzPhuNUAdYQXkp(x6S*yCCuf`v1bY&QIv>8p0M+#0gl33!lUrf`J(k2 zR^l(iM+H@_#QWGf1 zhVrG1Z$$oixJ!6bkSMh-B14f}IkZ>;ckn}Hs|`-Y;OoIJEbAVn+EFXH^>}}t(Z2j; zRGKxk8O^0pm4@vDR_Uo^KU1b3Cl5@{)QiW+;3u@_hKGhK<=%rnhd}LIpJT>;!(5X? z65faBOU}Uy-MTg;*uw+1Q-T-uL78ig-!hYR`_3J}n>F(<_RYvE|M*?>_pV}G6yqa} zY9mWUI58ZCh&M_+i1?U1hs+{K0N7V9WG=90$~=>yjIAHTAQHVnSryp+J*%GFb>C&5 zo?bmhoNJIW(FQJ9q$qjv zJY&>BW3AJE7v4nw(=iUfx4ETd?yCu~Fk%f3T3ffIJ^QDJ_~h0lm{h<%K&efV2EnbT za4>w9t^=p{l&-EWv~~(~5LS8ealcjzN4>%Oe>TG1OwUp>i1WM#Yw7D)N7*k%+9B_7eoVtXw zMmN&vK~`}YYwQpb^oSbYev3Y`@`2VjFV-oR%gYOPb=pRoWNmh8G za7X=axA6%`#<8vhawonN{Z8*dM& z%+V2ISSZ>l1?9Ws;}F15SAOs|zB_@Rba?u(tpa)`#5t=5)dA!IdUVMlO2;pCcjAQy z3!EE|Cc7b^5EXtp{wMQa>b~-=q{RvKo6E#!ED)TJ1 zWzB4wXYLGh4)!KOYWu#U&RTwEWzXI=!>0L|Mx$YDBQPKUe2sFhDGL8$G@ck%Ugef~ z1~tp)_A4WdKVogh)lVW4@++uK7J_!^@xh=XwKtu{lFi5#9bql4u~X;EKm<6o=%RxE z!Ux-3@ozUDMf}g@^|s;1pd`sO-NT2(@0`aNTSR|g@3Y+fkRQ`E0}G(opsv33<2&d*3F;`^H^B8Zcjy58)Pb&n7#I50-$7WUNY(iJ~Sy) z!tj1p)-PEd+@7%EW1Txz$qJ$=2a#2&$U5KXN=N@$s*^hg{};ES0Dt{jijO@d<0C&; zd7EqJ(bAGu5vw|{?Hfl`ARNcftUNH)0!g{Xj29{-I4Hr`M{84IO^r54$Upy_Mbqx` ztm@q)E-t=-{<+;jxt*1TJ<>0H;Im_^S<4jobt25cWEcah_pRjKxp0#Qw4YIuJNg-u zTSt90dRLA-cndgFn2!%W%M?pkQD!r`;ZH>>`uqFYh}UDeI($CD8X<+S?$M_3+`t0;^&t zvSQBtCd!}1lvU`#r@HC6iezi^k7PkM(jGaX)+0nIH@8Tyx?*kh&26%A3Q$j?S3?<> zkmRh*pAX*^kEV*Ld*585&T|6~QH!G^}h^x(y}0Dx**$qlL4 z{lf(EkeQj@*AKVfdWL*R3uIjh`rL6FygbnB0&*pwDHO6izP0ku6(DvFq2JYdWuKt2 zjfi`mAC4*V+JNy5_!EG^RCS*h`Qb2png(BE%t8<2Ou)HuAmFsZf`V+D-MnOX^b#4g zcyz;A85qdYS27qG7(_T(0yhvj6Pi`Q{CPL)s4T(_0%=KE5Xg$)ToR_`&!X6ZyOsNF zwWq2W3X4CSlphhfSsT&*at}3<_wP$YIsVSa6Gx)y&(C>ozU3sJ>@0y9=aOc|#0CSy z`E?u^U<%9pzb7V8=>H1;HJna#g{W6A%T^%%RsJ7L6gNU{48B}HMFKdbOBJC!KwWV< zt~WK`w&F9$VYkhi=*%u;AJ8u_&?4qyC3nrz{uEk7rxBg1El?Tq2%>)~X;*myn!(Hc zI_(pR36?qN8ZFu}jS@F$*~%oFSz3bSvJ)>~korm=K75v{uBf6i2bW^rjh)g=r}SeC zL~lWLu(>XCa?@N$G4_`jWrvAlB%v8HxWB)eK#_$@KsOi2MEec3R_`)UK3TYg&*AjJ>#m?S+dMcR*4smI%qH*;|+TDBo~-7;#H)9ewxEYhx_2 z#~)u_))y3QAc{w;ePSx5ad3}`Nw!&2#MMmOCg`>T&^L^~31KoIZD7DTmUY7_z8~A6 zO^;PuKE{Z|61B?wcBXw7R8jg^9&W%clhgtd3+ORRmvJHo1Uoixfm2<|E-i_O6#($U z_9rDV2hpln8_>@JN+-}=pT26q$ubGel}#48$I+hz>~+`D0O+|}dNAbB7CV03+#x-T zXgylG3>f+5b}zA5;+}YYc2H9DtxGQ%h~jtzFNoTwK2e3;6F-rUNGE@a()ot@qPzF* z8OkR0p_L4%qt}_POSUU6N?I)V>Y%8yn9E7;;nLSJT?ERe%Vj|e6X=L6U4~Dcw)QQm zli}SLkO#V!l#~WB&|SNF=eEpN&ppFzKmeMfn4|&}T=40dyiH#S&5OrPBHuzLu^5ws z2(?>KNblUQb6i6_pjj!k4XhiV>Q7A&mR(Nel1{{2LvZjTT?rB0sz6Q`y2UPc&f*EZr~N8yc{eWYV=5&sS4~RB4_3612@71F8&pZs zAZFnM$qz}?Q$di=izh4u1!Gzn~e3l4o7o(%#*L$T#we+62F|LLyIrM^B!VTCoD6Ub`XGaNd%P)hz+OyFXwuPH46lVZzeqOvr|RvU z&{S2O+RTk$1A?K)4)QQOe!!7{h#;uMlan+w2m-q`wr{t%XAX%Qy1prm78r|u5kCp2 zeMd(J)8s+Di`RAlKOu`AEoR$L;O= z7?@4BIe+~;dPX}GogP(#tkgT1SN!kJPW$rI6b^Y zx(Z~Hai1}=!&dmaqORgvnhB><;2AzzG5+Ut^ItO2-~P+BcPi|00yJGeUk?2?j6BFs zYTvcggLtpj)v|)XT1ZIEXqE4q z+w&ov&V;Os3H7v$-eE9ip9=giYhXtS?(*lgj#ImN(|PWj0t02KTwYh zi?z&8!aAH`6E;p!lEg$T#FAuW+|jVzdzY?z3`!!1Lpt~Z(5WNB!?90={%o~-XExGA z?U$iwJYfv5C?wnDzozfyntO+vnhcgMqol@G2Rd>c$c!*`AMn+92r*o>9Ry092QvOg z3Cw*@a=x3Z-))kkA}PrzN%u`(dzg?cCXpJj5->RwO$CVxp^3bdhGB;aTR>JmHD9Ho z8ppraKV74oOArXml7D=dt1|X3>t=i1x|v~_3oj9!9$tmLd-ecS^W5H@rT`R58Qf^B z0&KHUN?+0en8HowlzTl7PI+4*L`e}6 z^eG5Nk)k6>+~@SzrQ_tn^h&UX#|1xgtfFJVnn&QKZQ=Cy141YTQjjPG#;2>`2DcX! z7)bE>WLCMaF_8M?gAs^egGDQG{=fW4Xd+vrrl1WosUixuvW7~;iB#|-lh2h=KWC~$ zWoqAz3G@1N?R0;C!SU+d$9~g}%Si>>*-99Qn}-L@n>Q2tmqSUko^Eb*{#$2Wb(Uim z>}j>NwMt4#(qA;Sw8FbhvIV^qWIYtiEU_T zcvU%Ycd^Sx@Y{3EpSEW5c&TYDMfdNu z^%))=9UUG%+!{)P1Yy%E!y_WnQd7?_Ed^dq$Y5e(3JVM4GaIP>k~iv`uQf9cxcdP4 z0dB+6)|LS;db{b??+Lt;l!WB^XaNqgcXX7=?@phk92Xf`R$3~nqB3a08XFt?`t@t_ zfZO*nGKdJFd~SQHT3UDrp>j$}nlY7Mw63CmqCD=Gl{@rs<8nORv`xv{QJzbl|a z4O4z*N+YiYV=C(AA1cTeDn!`Uj|!>T458yhddRS>1x)dqz4 z_)WHR)uaL*sD$iWf4W2Eb%&zAa@lIRxO~^3_y7B+d%8^hC49((mZ#@U3?YZd(R{7@ z;cSUY&dFj!MN3P|w{PR4qYf*r?v7oL+z=D+Vb0FZ;Q1FXUJ#%b=l*;A7#i9?-yV)i zNS-VY!1T^@-T5tsK`ksy_Nb$$#~<);ySKO3*Vi{atx*(bZEf8fj)sDS^eHMzUPWcD zYo+t?dP%_ZWVzW{;_N>E*VEHQ~dmZ+) zcB=`>qI@ACA-lkPW=4}H>F)OS zCkG~S*x_7_w%-O_Ig^Q-v07zG@gOoAl2tC#??LjwcC*3HkM}qD`1r0Z20T1XDJdzs zbrxvk!r&b@)AL~dhq_P-GqqzkHJQd0%A)}yNTwH8!Zo0a<;!-gbX+Bb>6WNUC zirt*9si>%Q`QA9_>t{?3c6JI|T3SMaKo~~!PIRr13HjjQAcJpv|0SVG?awzi`HD|xO7!x5E^!ldefcT0`&}B>BnkQ}>BBLYiGjlH|F;jjZ6vZ%-B} zz`??XzamG6=ydws^+)3qq7HrLsiZDUqeenj zZV|R|dPhU4a^lDE9h%m}PEu)&qp><>k+USYn>)s;cAb zP(!;Rkm;zBRAxhYMaAU~-}Ush7aRcJECT6y)ctsjDw7E%7+6K7bg^b^a)x*Yxk<>Ego9 zKMKwlIgdl8nt6lu)oUs$vxDg}uwOJ7NFZ{xp#-F)Svfg5Sy@@>>E*X8?u|>}E!EZ3 zCMPFHGx-YK82I?S*L%YLB@ZDVH#xb;VoKHewHM)n6kk+C9_wps-&xj`T7WvVT3A>( z7)L%gJsmGw2oguOu)iLCvJptBd3i81xhgrrRYfm;Po`2>zz|uX)LBjBpP!$DV2mXd zXsoKTA4+@!5}lEeQQ*T_L_TsEjuKZ+kpkF>R?fF?_kRS#tp(gyb#--RzoB1E=dw*v zq*0&bp~H#F$ysZ0-EFj62vNqUGVT&)Wpx7C46I>=Ruwob#O&rs={6v%a1z;r=Rwfc z)zyJ|S*m7U6bI$e@bbC^87<}O*S+<1knyyN;=t>OI4pjH#fCtlru_Mr?5_@H&d$z2 zHsf^a=;;oX!bO_wNC%uh*yF8MXCRd64ou4JYE3-$F)wwf2(qTHb^eO&rd@` zPoj2FObdq++0LNWOFio zP;}Uw{ZYvE`RV?2dwZLWjm_fi;!(UTb)+OLTKzbfJ|U)?09B_x);5;sApRmu1G_bAo~OUw@I$I;3}ve zE6dATHAWpDhVeAd5=JbHTiy20*S2X%D>e|h=Eq(p-=FFQH8jQ@OWnMTxqUI63YZO&emM?{(MX~;a? zg|&e2K}=gyGfU9h1yns$gw4?`0eluC$vF0@zMe1?Lj!|O=m7)va2P7%Nynqd)7d5@ zXlU@W6oJcp@>{<%shypXKv)R3JyD0}Y4g)t9M}X`5UAJJ9tB)h95K0%kB{JhkC?DR zAmA7*R@GfiO-yueo-rCm!=Pr2oSd9Ax-lWJVnr5~l?_cc5X?`{%;>fv@4$@WF!X|b zew7kI0D=kJWJ+~E1R7XEXUnkvAAim!j)eQ1oBTRyNNYpRBPe8WQ;s@O zhmo0yDS61;Di`c!q-09fyeN!3N%6(9EQu3*l_Ka=W2M!hVbjyoy!`z6Ptq}>jCjwU z(%W|2`iyuNcYg3<>@FNpE&W}1PxRPB(lgWQV6dvMum4Degb-TmBvpkQDLFYiyR*Lj zem@G;uR!P*344T>=8e9J849@79lw>}G%Ik=k*$y;#E9AQk^%oq^N8c~**+PR8~4BL}e*wLBw3UqK-n zEOp?_bWeR_l^8bF8{sTjpZ4g+s!Q4`l8bTVI?01rF6kY<8* zOw|Q?b|)vNKTf>mP4Wr~!@qtp8uh~XzcVy^WOS*ns2Cm_V_{;N5lK!>O_fVx_{1k_ zTR@9~$Eeffkfr;J3mrXTYGI*lpE+kNRhAklg!771Dv&Ho^Xtsyq}sBE;Kkl=Fg=0< z8v4t$ZULm`J$ECcY(*Llx)xdL@`{Qi`67g(?>K^FEF3JZ1M7rf?7uo$(X6Ij+ zRZ7&mm-^53;MjvfW^HW^L=?p-9(!6$zbT=HO5OO}T%L<#YOo4EM+#5~mB zaoh*0H$*W&3)9!z+uPkO*3c3kux7&b1C=v`0vy-(@83(uyzmmt2F)nAdXNMg8fMIU zKr5?3Z+Owgb*k*8?Y!gq?c&eQszt4aomJj*>P29U-_0+cS6ysw@GH~PTc@Y5AQVYM z;(4m-oEj@ysQaKU{VeBVPSu}yhq;GowRiiWaO~$YbosdlDSdDg6rm3@AC!bRlf^9T z+UshoN|e~d9UW`HaugJp$;d=_M}x{HjnOxbc7qi1B{el7A_A=BtK1T^OZv|#DR?wV z`+xqJfH@1THw?xC6n?TSWm>FA>G~{D3UD3lk{(@00q-o|si?#)S{X1meKEEftoP;H z`B0?489K1rv`akQPd+0oUnM`XmL2%|A_Sg_%?h3=S^iyD7B!Y!Qt-+Cepy+WwWFgp z>e0c0cryFrFgr@egONZvX@V%H>;?oPc`s$DVflbBJ#N69Dof+gT)z%YKRkWY-sC(v zurE`*#?Vnh5UgARa_!aib@~Lw_4PG-nxeM$Tw^2mN-ZuxT)+TSG__AgMh1WuX&ffJ z=mk3u5l2{7y^X!YLsu6U@fu9jeBR3hS?b4A8P6a>EA)s^`BnfejGr(h$lhW--TH3^ zGC{BE1afMj(=oXp*waL5Zuu^DHx)itjCPEE`kJ7Yn*aRAkTQXKj@nUiQTDa!tdgE5 z!%m&fmZ%PvH6BO-s1PBDu}i=f29+!2N|>bUE;Fp#f2 zHmjlMhz7x(?$Z6XU=(H?H*fLr@i8zk0AK=H96l0!gJHwf%GF*yb!azh+ov*8;T+?6 zEdpYUXy^tjoqdpksFZJi5XU$PI<0`fQm$d(ivLRL1#?aIH=;-I7KHEvpN9la9NW7r zdeW&_$Xja92a;#B@nZGe95*`%gE@pd&i26>f`R`1u}ootgDG{`{tu_3h@k&u zuw32ek#ftKiUCl=B}k7Nk({`vDKJ3G5`wKI=Q zf=`G9fI@Bs$S^Xp1Zhrx!$SGo`DC7X_0-e=O5wP?T%$_&#j%UsE%K8hjjztxxMRiI zN@5!~H+Q(6+O!fZXbAO7HR!9#%3dNPv)^YnI%6U7(?f4l<3zIV-X`_R(aZdN>4y4|%8C@g&!0cRp^}TadM`&y;Wu?PlRoYp?nQUMKc~O=lD|UNTuc8K|G<=C_5D}q)^i+? zXdOg6Ep(H=Y#K`Dhm9@XS)gZlvt8s14>_i;8L|2VKWrXO)R03rjZY$O{Q^8}uhZ%| z)|l?lL%Y4ay}3a3<>UnanUJl|VmU^2M!9>K+yy8!6sJTM9AYrI|-iyFYwF#PgU|mPm?b?C_9*Km|K`SI@94|fLUHy0oDFpM152zF9|7WTWhO0 zVnjgK(fyk>7Hj+l;yi;H&*{K`n+xOFhwt_WjdOz?Pv`Ne+F8a$9LX9qhOW*da~EC3 z7<7=JkMSSVtJ96LVFa{v7*n#c3~5VTN6+{~QlCifrH&ffyBC6A)ce--$| zTEM@#oUvO^P88j3Vn@dv3XE;B@)u$mwrJk?XY!qis5QKiE;YBQS7X`U!ZEzJ5P zH!rsCWX+Jb%b@B zjB}m-Sv!$tBT$r-{0Q2vc-!${Z$#zY?Ci@K>O9HKe6NPgF`+K!U(zyA4a{Md5f?- zuboH7u_z1h9yffO3ov{>Jg*%Vfo@aLB$buA01zuIJlw& zUSo2mn~)=g)%VfWF|&8`Z~2LIGQkPogLA(b6Q*jSLgFUeW-{!`(TXG%W^|3jPhW#4 z-o7hibWm<5QzG?ejSb=g7_FIG;{DeLaYZKO+`LI>X=(e-0*O^VjSorv4XdbY%S*$ z4O|7DuP3zQMsozoBYK(0+uB1{YM(QSK2=f-`$QxapgaPA}NM0X5%lPAX^Am z_sTpV{!ZXkxoph0Q$zi3_JN1dZ?gGu|=(Eogd|rI19uD25FQQ(9A!(j% zeroh$u#y}Ga%vp#+`)Fbu(W79X8CS9f9+Qw$6k?6r}Dtikld6JnP%&Kt2H%I`4ADU zGlmm7Dre5vDW1EC7t{pIQ>V0@#?xxkOC@vO7F;V_z3mbEF~$8WbP2M4Z~qq6 zSzAhBVIDZhSvt^?Wq!pLMvT5IxtOxz42k)xcWIl!>@fMoX3Wu7YoNNm==l(pf;cL+ zDb!QUKX_EUdky7oOi;Lv(&GgrNRg&sPy@>DC!z>-Qv5-1XLWi1VJ*Ttj0yUFB6HK# z0cSa0zDS2rtNJuamYV7_F}LX_z?dU1p@7c0;<~yBy_}pCp|zw$W+<93wEHU&@09_R z@A#!7-fAWux6%~4&-8mV%o~+M+f-f%#0(*BsA}%@k%tSvwHMKp^ui2Fu;{0&sE+JF zLj}D&$*)WC0&H8x3XkI>AKT#FKKTEy7Qo^J`Z#ntks}k)9C5pNci?8mYd>GV6=fcZ zS`Xldb*H8%harr#Cxs8}6&w-X4wn=oc}6`bO;&X-W|*eRy<4 zGD={yqW0)IJ62vMuPt;mWjgO1q%_S!|{08hQu2}JcWh^Xw0(e zY6k}gEg>EU>)(@guA^jzZkKA;MIY8!e0OIBoP9StTW^we=xM-YS4tu+OCp}fTEnN3 z6x;}y{7UN=DngNF3Bzu* ziucC!Z!|4_R{Og%cM1J{J=n;B4h0 zzWiZiJiN-$rm8FNVtY5J86#;TO7||=Qiqe9I+Y_sWvN>^mMaA6y{TG_eT=7SFhjc) zktPwG-GXxT15NhLV6Dw=Q3D|;A1T?jC^ZcY=jq=7a z#&T9h;At2_)ka;4`zyG30=7{@*&qQ_tU z{{H*>`#v5|2cK<3XhoDf?d|OR{QbFB{5?=n;YLhrr^zDi&CKdpP1ku|b-v3$7K>*I zd_M2*X5r`KvrZ2vQl%>mPVqF&!a{^l{C>k&Aid;~oVl1e@$?Ip6FE{c9ak77IZ^Yi zEmHyZzC)u!c@xlM0MZu;K17s4l*h`thFU%*P|?}4Hk%$Xb$FZVVX=B~qO)KVcEenB zcGPihcXwibiQo@jkP#TU;oA`AYN3wsO#fial(sg3G(~PJbrm1`*2FSHevAE~FG-8) zLe>FTxTWo|dQAMUVz1867ytHzW0%oR3#SA2d)kS-)zeDrmeH4(3mezPIYzt?!|CJe zPKjbfB@w=xt&Z3mcH!STYz|V^)@8*ak*VSBZaZg$S!?;BGu?&k$+I}7wK__HtLMoc zE;+**lM}07at6icAlNewd~^rw?#q`i zffgoH(EGefze#aS{iP-hQUAaILT~TimWQJ#yO}mFG_Eq-GWFgPA1^PfUxd;9zkqaL zmH=(yEr>B`ug_i@XR=c(!U5!Dp>v#}GL#7q$Bz|@71IuQygtu+lM@Ov^6*bjPj37x z{ryEC{{#-MQ~68dHQ^reU5sv?RSUL!cUhLd+!VlGyE_|-Fg)9;ESaf!`zz zbqjYF_{#)B|4b>58qb#G=0a9J>lhgPoAK5R13Uge+sD1HulX$!--U#R0?udx*8((8 zlnJuxr0t~b{aCf_21@2LAEfT34%Io4`PYLxdxzk(RdoEBFETNytZ0&Sjf}{MhydA9 zzU#CuH=OjYv%qce(AY`9I5;38(TKSL5c!h@1Ad{wHZ}TZ zSPZS>N;?-9mu5ip;O2-Cn=kwa_$=|j=f@jd^|yvNxVS*b1$Hw3^*?g|axbCy=U#t; zzszFA;=WhwBCqG^S{8@j=yhBp!>mkz7-a&t{bJvdnW_0#C1Reurd8k54w2^vbqb85 zXY;ORT}7HCt{KA6?shA5=(`Ha^AS(;MC>;OTO9%gm_7mGi36 zkHQAi+D20e|FrY|Boy%#08$ruM@Amp4+oZ8rUpzjup;fFngGi`?g1cpP}`vOV?dGd z^l&luvKt#cz(@5gMUt`Q>GON70yq*{^umLM>?YkL|~88Ee+^OMukW z*4Cae;@YlwtF~Y)^LGY-(!0BX*N2w0>dJWd_ys75obSMa0SB%lJgVLEtb?N~Yl7C| zd~<8xSVbH4?)*Bnt>$Pp+h%M;NG4E-{&{|GZf<6#h)8KF?PHW=AlAq_p+~xbhX)u| z+~2)RSZvc9rLF`-3Axgpsj4|1dG6q|J}orl@vLh5#tgmKn~Uo*S~Q?fmInf(D@Kh#!2fY5 zuX^KE=rpXhdb}JSp}7sYT8g_#(VdSJ*W`fC5f9U~($@wNZ$XQwJEOUg3?awF2$PWC-4IAL{3~Rx=$Q zKnx~S9Vqx=52yc!cD`|ou8X>F&(Y!$rZCkOd29;kOII6Lx2v{WHy_bj=>QWlhm?l3 z;+UwPBXIVkZ;rX^x=^o%sXs8~neez`XFc2};@ZV`rt>zMIzRj*T=i5{h@n`0BR=$m z-P3}%CeA3(WHp$st5iGR$h&>}Y-cEA#Tr0A$I|4C`7hYhjznF?16t!HX)xu(qS{I? zVcf&YFAt;9K8-jI@z9auMN8$Wp|1*g*(2lMn!mx=~L9KHMG+YYqdYWyc2JTj{o8^ zt}7u0I;PKSe4m-~rm&PKFbcVfNz?o!*~ZiM6@Or*`(~`zC(l*?DmyGkae211A&VO% z_W%2_tKfNRwxriT;*cNU>6y~UJDhBAwY;i4K;>72(bB5tOoUa*ee z^=a>}d>ycX)=raQzfo*yt2JAb4N?Sq+&>nFWtUqS_vayniHoH!{Pk+b?0$ZbSzKHM z^aeN@VPRqWx}{IKO9vaU8*21wGzb-*zcY+e$@d7XtVRK<(h&XmnzWvbgip8<2ivwNFapI$$%Y9AP{JEQylWXF+ViV#x4+h@wwI;esI4j+IxtCg!0 zIunwlm5I1qG|bHlN6Cpn7v$pP6dyQkiHk@UN50O#@icW`-c-F~7z>6_wfBJe#*02Y zIDmzP?ey8M_#MZ=(P&v{CwJ#?=4pr5R#i1Jqq$)DDYBVu#Rp!9|EAozS0<YRm5J!bR@b;UCqw4S77_m?l~H}Rf90C76a&?x*pD$ zIeSKW`kouKjyL}$>3c&xuRvK)_HJH>ym=se6wkOP@^tbeP*kNlcspPH=WMSpy2PQr}#6Q(>*7Z)RnB zCFUxXMsHF$^-xxP+4-OSP^VOTehd;7Dg9VNq>5ie0g>sMHC(GeXzRDBoi13*N|+$ z>(cDmhk-8t+M;}hm@$LnL_2+w4!zLJoya`M<5BqD&u0V@%Ic2)Ji<&N+ zCoCT!v#ppIC@PbHAPhw6N-6+-@tVDz{Q>t)&9}C|r4B)tmQig{smp&D0MMcTy|_>} zbmcScL(0-io4CEb{UXBTnb#VhDQBboe%5Uk)q~KkdLth}N~n7JI2*O7rl#hzdP)6H zcjw<5QH{t-Nz#xYEv<^Uq`2ZXysbO7@f~sUNH~aSk=mx;b3-JiSiC-N)2gzce#u{a zBBFASuSbCF{4mt`w(rHjz914!4V&X{6RtM$w>ZB-{=Qffe$etShqmWANZ?8+Ch&fO z&)2Y+C#(cA)jItgp>jWsS;OMz5ZzjxWvYYPqDA(ffpWb%s2)9kUTU+?SmC&krHZ8F zyf2Q5jaACEH!;*UR6jL_l>+~2*=x0Xwu>R?1x8AsKpE(z0q283M)qX$P!ag$IS<-3 z-+f||_|@XUmtddcd#)VID-?4DgA?zj zlXp?)`6pF6MjMuNjLRB#;hO$V5BlvdDH^(5y6>Kd)|Mw+Db|#{JWp$U|E2td4%hXp zF$@3VW))3xu-A;kU&h~hDxAP3!+l`4VHeL5k3E4ahxbO{y5>9}OdH-#sZaK5ezeUE z?KTn5lmMQP>CV^N)7z-dMCXi`K2yO!e^rb(yFHiD%U8# zazE}K9mW&v%a(bD;TTi|Dht)$kL+C7_i0>dd81I~kR!|g?kOBl zuPrS(W;EN++-X{CT8o%Js7|U;<~V;Br(C8a)3(Mr|9<#gB2#>Nwa2s;@Yub3z04Md z8ZKrYheK8G916ZvM!v-Y1>&aRhRJqfSoFe6kTOkY%`k4XYmAIdpn^$ZOakGf`1;i6 zlJ)HMs9;_a2fgw73+vq(8nKlnb=k$$Wou47P0g7Or5nH!^UsqW(p1~p+q*kYu+ipy z)TRFXv|jW#^<-vQ-Gw?{Gl}c>N|4+;e-Z4I0G?mDS#OLiM`PIC2kUqZ=q zieajErI%|3$O^9AE;9dch)vvIhnBnwF$}F$t4v(9is}pQUHAJc?DH8YBaFN~Tb7S^ zF`8Z{lO`|vk3A)&j3rU3v5V=drsE5p?vijZC&EDD7|0^ptnDaqfIa~zE3gtS9eUa# zc1S^O0?;0(ThmX`1dyFTh+nI(t1w=RFdR3FGaL7Le&;fwN*4a;8*Zs_@Kq89f)v)1 z-kH9>+N-&gsa*2vf*{2>b!}m#Y|Esvqh_UdJ?;17S9H$tU;l1YgVo(B{!*5?vi;j@ zXE-zkUuop;y5Ov{SI{U-t;D3GGEEsnAhfh4sS-v{Cu65up?gAI#|!5@ctEyu$!;6f zCq{uGNSLjnpL2k%OMzj|X-;?Wb0T1(KdT?w3;_ii#Y3H{*|)`;FDWZJBU;n(O)Aac z?SR!-%8>c5T{u8M5fc+zbEfAp2qQ-I!9a4c$}XsObH*Navo>v3XNF*sM(9KfB%x9} z(>r@H;ZK%vf64M1VFLmJ)<3VOPBAk;2W;Lq@T!tlHg@OzC-1=ZFD+$AQYI%qa&7qH zF>L(~OX^5cq*rIy8p6<(F^*{K80ZM0dma3`F4S^Nh`T9b#+C2>oMKQcyD_`mctJ8U z^wVpXk3M&+Ko_>(u)hgAI3!^d<@L6RSt_T2L_UINxZAYWHDXNYPBiv(9Ed8oNSpHT zfUF8A?skd?Tq`{g;=zsR2w7hx=QA81b@|yC|7?%n%dzr+BNe7=&+ibo26Ty^?V##xe+?E}UUp*V29vSwu{h@fwr~CHHatB1N#tyf0$@rKBogel z(~P%)lb6)vZ&i1jiYKg4WLRR(q9W&U&umX>az-8$DKJW_i`k_%MDwu9^hP<7dv|Kz ze3I-`Mu)I^{GD%foOd#8|~2Zsjn^8P(MR5T=C@%RxdPnFTFo=dS3)IJWYVA5ZZ zL=bH_upgE_+Fx|04Pv^a`(z;T!oE=b?Qm+PK@>uEyh$eRE1lUnXEn256tzp$a4}4-XH(*HL+`LJ%fPjqXZTTwjN<-93d! zwIXv+Fy`gsv_AHqvM$A5iDjb!g#|k$Ee6eu@HaF|ER4n)21a+{OamAnqMk~ za5plJ4fMgimv8N7Yw^1}>sF^$v8{wSueUK-N;r#i?`^Ahz)WDrj~>e5vw2oGAYbu8 z3^kEg^_z~Jl~rk7-33Q^V#}xVSoycwX1U}kA|xmL27wjweb|rS*7Wrk>>A78vPpps z4-C|SGS8D*D-LYudgnM-XkTjK_BNcoCKYQREN`Hv(LZi+iR>94-%Jj69Mg5go7 zE$$?z0~Jo;UxT=~i0(<0-N<$P+t#ueHSIh`y5K$8dtFT;+IOmerk`Ni(L;|l#Ax{= zz8aPf%^%1ffb&Fe$Gi!j(X)i0)zgQ|lhA=yKS$kwDLSGKQRRm^|BzJp&&B)2eco2F z2t2b|9TkPGBhz>CqhF=+Z5P@^j;R?XIybu`z`&xtMFxde;` z@KRRiZq%&U;vVB}+ixdi4?qog;)CCS9U?$QXq{vpe-l47=#4jO4mv2r*lA%E;N{w~ zU23Szbj}sX^gUjNRvhAI>v(M?Q*n}v*j%BR;CU7Rxg^jwyn9Fl^7hLrvqzl^c~tE8 z>}i$b|hy6ZA(9S1tkRxd*b z#qGDe{K+RDHZ_+GNDiu<55)ou13p#w_)Cno4jp}+OHM05Tdlfu_SH|=I6yuKFA^?HJ~y8TKVu|d z6b+~vSAc!0Uecsi=VEs(=Y|KU8-X1GDFlp*VApNia;Ns$(|&KvcuC_^bx4EcN*B&+ zmL<~~(`Q<{L!Vz#g6N)8;}-qrfT+`x<)2)ZU@-{=J_HgG)C%Qi`jE_-9CJZha;;Tp zX>PC2L}Y3k4J>P!qCfcg*0-H0GbzKp(WdU)xoPcaetmNM?BjZR)B#X1|o3T z%&CCI&o;boZfnaZbf(%VC@+_vssg_|=}Hswy}m!~ilYGfBcPo9n=irwgrJswLIsU= zcxZ2}nRWe4odIqn@+IKrxT&jfYJ%fsgA4$3481KrSi6>`@X7=TeBVtKvT?-4!BG_8 ztSYF16N~Tg<=DGqfthdezL4tX8--wT`9WHGUoQfOkJ1FYnh{B<9qm zm6a+V6pBwaxcY}| ze-j@I5TTXT)XZh)B_*i?ai}Oo?u6A?y>6iz-12U@sWP09OJQS258R7QEoW350#xaI z-==7XRUozN0ENIu-Gt|YZ#wd zY6^f+@N9&JrjJ;!5lpPUbV(yXQ87Vp3MAt)G2aG~ zs<(c)kN*-Jtl(7zUYlgmj{B39w$qiYBY}lgqnf(0tHC}^5(or1s#LzuZ4`7w%t_E% z3t}cq7XUv_gm2g`&lV6|7PWz1T*5G-T&WZ}SQH5f>03jCe=*`0IyMY#=gU_*j2`jF zxs!NY8DGJa0!sqBWhlN&j}8e64TVFG;AmUIxcFCv$C^NAB1+o~;9L-@!00gY#BJZ& z7aY~4LD`U6@Q@+CDIS2xbn3Y6IW)?fv~B z`I;L-Y2SoM%o$lI*=9ccft&FYb7@?(HX|M}F7AibIhi^*@QKp58w5J6fWgWT0xn7b zhHblyFzQH|U6n5g^mT24mph z;1Jn7B?Ty^@;dQmYS3qK`OR)u4%?{64we2QI4K(u1p>$BJzhq#Xp^Tm5|(oouiQy( zwjE_vH^xIUvh-=&)ERLZ!nEBrS~#B0QY+&6*eHP8B5l_Q;6)$@g+VM^|8=VqR*JB7 z5i9CO9cUab|EmR1uig^YV8lh~7Z(V3z@!jn?F^gBCrEsk#5Vn*GO;0~ zE(3Z#xwa6xC)%7Pj^LuF8}^~Go|r&eR(O%_9WmB_ z-!yz4ZK!leRfF$|?vZL#MER@_mCjkP>A{owvn80g(2>J{bEFrpXUL2F)fgC6zN|qD zF5{f=J5V{%qlq^~q@T`wy00)-wL122U&xc_fjO5wjqxp~TqkB(#*$JWb|gW5esDO} zpexE3Ud28}EMRnx{PYPF_=kZiV768nuAhG?4NEVhGUefyGa0q$hqo2W)TL3 zoIxZFRWwopp8(au`O@?b<^sS%9Y9)Jt5bPvp7%d~DYoc;f&WVMD}f+e0e2|vMxHR* zXA#~`n+JM$u}gza`Zy0j;YIe-eq*oX^0kx!ek%{402b)YVlpYvPg5?=%p~VC9uk#^ z+-w2&Kx2v*DS?g;_dve5KxKgf4_hG9?I*kjl^F+7 zkDfGyii~`IbtTmN*Cf#qOG+9B@<}p(S8ao}SC$sIgx|a=0&@KAU&(Vd6>Av+L7+0t zio}_ZfNCB(nlpn9G^bdR;1!}2tQ$;-hNYIj_S)9KMP&m#* zR<_~&!uzxiSQo%~uDi6TmnbXA%Kp8))QI|FQP%{#`rD2`M;?w&8b$F*lmd92>l&lM ztPZ@-|1mpqa&iK{W7%pC&~#TCM!&&eXhONsA}$+1&&p#0UM4hDLRMDR>OT*_RFstLetX~hO@}cUeANe%mOG$MWvQ_e zWa|^5sC_8Fh{Bz_!Q>$!0ls*lr@FdyB;(l64|N2*E!cZ-J)k=HdU@@oMfE0-@TrZj zGV65sNSYVEg4))4-CK2mQw9tr;(0o9{sS%U>ldWtOu@WK=Oev*jO{z7DqrvdqtjTZ=cmOL#|nr)ZQ-W>+4IOY!ly? zkhMaOHDqeK>@E`~`g>%1I2F7K?4!Mbi>pl-G)L}hcJ}8P-}#(`1R#))sru3_DJjXx z#RY6>m->+wz@QDR>$8&$!1f1x`vC7w#ImPFCXxaE6Hy8$cg?`bWRC#oeho=wVe-GU zCjFTtEj`_)3>o#4)}*qyc+X0k&9Za`N}0yN))q+hg&Y<>KCQq!Gw~wL6~2lV7=pmn z%!oXt8V4-o?1lrC?A6tkbIWo7!vv$bxw)j26mVby_ntX>q~u23Z3GzgO=}5O-dn+g z%mj=Cd#NJ<0PFJh_MVxY)hl$%Qw3k+u$W`XDva{*_z2_=*@jQRqC}Zs1t@oZGgH>D z;A;|I!0-8@7Q383@r?f;_RjLFt2TP`hwg5W?(UH8PU)0x1nKUOZjkPhZt0S4Ndf8Z zl!iIaZ`QnC@kFTj6i2&G;}@awju!qO7t;hdz%+rthxWXy)@3Jsu} zes|X#TFlhU>^~yL1zS|kfyjCy`l}woKWe4w;s*s^cfZkgDpYXsgF=SFvLUdN3EG#y z!;%Y}anL5|ZwR6~K){M1$@_+Za;*v$=*iC!2@obH-d9aw z-OH=3b08%(wW6!V)sHzuxLU)|5al^(MeEZGwQxtbxgo4xVvc|ZsPr`I^knDbdq`Fh zPW%i7%xE2|kQOgSmeS~fzGzUealYQ=^Y2tQQquIuRu6y;`3km(QyI1WLGlWG@Q6T6 zsU4E`lKA1Ayh0jNke{3bBmIyU_$b~ex8QaKXQI>V#u98x4W59-d$3uRp$OZ0W@}pr zEr9|8=kD)RiF=+nn@ zK;8$4i(oR;MXEC$j2aq}zP-Bv`90UDSOf$W|DJmOB;MLkgT2P(^!XOd6|0%granH; zoW}}_|HCnYhOE-v20<7J?|d_&5}@lp6%vIDG^p#)|NkHVFQZ|~4kZMx(-=gMm)XAj2Y)0|fPe))WnfX3%E+n{xDU=b0J!qMPV$13p-EmdQ&T3*ytl2RxiawlO?`(g*Er5Xq@tp- zR+&1u0Yh~Kwk$;=mz79i0^ ztLy1S2_p;buzv0q1kV?mWO9(*b$NL?!#NNYg?Q^_QojtkilN1L9M_QnK9BS4(a+nu zSDSrchtAsiaHZJ}6eo0bc7l3b5HzwIAS)*)2L#@41dxg)fPN~F#G=zaMsr1aF#*z( z!qfORy_;89s(H}a)ze#8pulVnq}oheaaf`&0)Fy3;#Ty70?cXv}^+JH|6O0_h{ zBprSP2SbvQ=7DmoiA-*=t`Qm$aj>^na0R*!hvG@Wq7f)>2Pq1mIuRXOPC?=3>Pol6 z-NC@%L-MFWllAQ14i5=&aj;np4jyRgTls|j>}qOC0gy2=Dv-R5G1k&8CQn*o!H22{ zKSIA+hXzTXaDsX<$pS`38unX#MpNq*ubX`cz}dwG1NuLlRP{fBi%uv2L~zNVn}cu< zAP{6uj9$o>{Qmj#XF>uxSbQNHiNvA>08lG761XR_o}r1-REnEmW?H+L5%BJ6qx8DQ zWVq6fd<(LePyvu`1R{JQYHL|xLu_nq!Rs&-05V;`(~bz|8KjfoaH?g2ZegTXE0C7L z&;Kcl?<++x8QAm%o8sW)gCHjmdjSD3={12DZtyUJ0^p(h3#Fn#1pxqnBYCww2o5O_ z2r_o09O$F|9)c7nAP_upDd4Zcx(XEF|MqA#=)J$Y3jj3+ZNdH;epzkacR#=?JSl09 z*)S^53!a|~?iLa~%@pCKUCw)mfH8RZ!3TqA6Jp>JtR^v9{VY2uL-BNh|K=i(rMe<7t3*LRJ;)yvM1p5L zh|~fDqna2GGxHf}x&{KlF$ISg2Mi8f-pycM%#PJuF=G|dZ3^05K>xLehX=G+9H|H@ zz~_BS7fc4GB^(@>0GQOx2nz+gfe3{e%-N}{*ZU190F_j@ss(%S%lzNJpd1CZI6NsuJI z&lqG>fz5s_tR2wlf*Fa;+39{p$S}V!#-TOo!E_r8I=GEyjJ2Ub|o+6zF|> zeZD&v&>9Zt{E}pXUNsjw(}MJuKa=T>2?kF(HBh z_$0DErYO?{iyZ$Hi$~X%u`~l`)l>yd_fc(zM~@`m1~;qdB%8=$DuqHiShgSt=*aAw z{Z7tCfsY5>BI~DFyg8*->{-2!*oX@^PoN~SYvcT_(-9iTP_%5;b7}hk#Kgx#^@1q~ z5Lj3qOT@b^4g3!%sabDw8CMu7_>aA?=jA4NDp*=tu7eUu84@pO(r7kzb_f6k1!b4~ zJzW%xa)0!Tmri=_svZ464qOB(@%D-OktXY@PP5)?g^J%X7#PjW%xG0|eP6D?ddMS| zieQjp0GI=Kf&abd^~~}4_3iWf(~`$uH$1PiX`y=qCfw(zrJlHHyT(7opXd)#wCIIx z`kAg#olH$WF+0Bu;@Sv*tc{Ybt-`zY=z0E|_jX$*@-!O{-b@@?|ArHmm6i(Mw@ss< zp}ijjyd8j#X9Ow{i`OuOqtKX{nE|K0F&sdY)C91$wbff*!xzooyLf-QFsok%bC%KN z)wMNb095gxwuN5pu?Dm?HeyX0P~rM6^bh<0Ys{NGxY>Bw=XGgYn^O@%0l+fvf6Od! zm_QB-=o)WzSapBAwum1NnYemDy6s0|1#bf6`@w9yyr2|%rtJ&Ld}v#too=9pHkWYZ zaTHfxwzgEU1+R@tO#r+Of+D(qOFeIQ;O2Y}f+MkZZ~)_7l~yIQ z5B3H<35ZW}^uNgpxZw@hac+Bheg+Y8ceisYM9AS`ld-c0vofTb45`z?&){!#f%;Pd z@IvDQVGL4m#?zR~snOwm?cfVXq-@=f}%Vpl74JwgdSe24p)_j%hvi^FPoZlM#OFQ^Bl&GStv zq+9cHb;aHvovIE}$z+QdF!}Vo=eC4Z)0G)1;W%oDR?RTc|MhO^<6b8NMOOp>-Yf_2v5JAPIq6Y)}{zWU0}O+ zyn^9Izs3IX?Y`$7eCF4qyf^xgehn^OY0{9s%NYr@y_?4fQvdJut}ub1&L~6C{q9(E z4KsfDJ+;?mT8P$?2Dk*VeqUBwl+3MuMpohqYZa&pg2^^;nKA7vO#|AoKqTSTMZo*B zaNb3F&&wVNSwth|4_zMt^HFHPyslrA$@2B-ZkyFMY$us!B2qFK8^pP6Mb*GCq#w%< z50}@--i&&=Nj_O(FLT!gDX#ZCPCD}43*Ss|F1t0sVuA;a0JwsjrwhX|=R+SYRdEt) zy7DU{KffMpgOP!gFm|LhJr>M=)Un-mM-r)Xml0(G0z`oMFXK_-O5Yt}7ch@fOSW~+ z9-@05qQMM>_5GUll$snw99)!%ydGKQ`MF2Zw)9_}HctPSe_yP3jSD3bMwKJnOy=uR zD9u*#NGPFj&#Mkp#}bCr6oeDn>J*fB<=d5m-)hm>vK|uJg5Tmqtsu37|5k1iR>+?i z3htTrxQ0;hw|<*}pOUZ9M^iQVfG`cCSWtU7dcd|BEKD@px8he{|M*~;Q&^$#t&XTh8ZHGOIps$QSa9;IQC0x2mcc3#l)2~bS=eIbe(5Pj)W?S{Cz-_`o145L5 zD5HYuBEdA1xBG+W`iCYl;Xxn}nuJ`!KKX|JLpd)*uc|aJr!O&JTH11ZtBSMtF-ddn z(%dG{m8v(*D!xClSOYL#qz3?@5qz5F8Louo=Y7l6Po8gOBGR!7CNr=hCZMR^q65aF zRCRd$ONamGj=>r*#h}g1g{Rhw8I;s@&*?8rf(+y#m0@O=>qD8#njgT1{yk#pc&QHw)L}OgYg9Hvk@)o?DJoKEV=XD{&g1u$kb>L;@{(d1h?44 zfZ_UOlnovDl;eWamjeTY!|C~oRH$Ajrn!pLWM*nr4k935q0+d?SdPsj6Tc(oOxz7- z6rA1>D`@hv*5ubwT5V&UAqhjzP>gRxVXK!|lL)t$;R63UUph9F;!tB5R*&ybXJLk+ zQ8E5|!aEcTE|}b9sx+&2R}X&H*a>8aUi853)=EtmjM1Acyw|j@_?bjlK7vw~v7)43$W8!7_A(ia)6{GMRvR+tTI zsMHPR#%^wdRffW-{lr!$GHE)xCkuz3PylQ8A4(0`jQHv&$e!%*f+v*4dz(6liDb6` z)g`2oF;Vqj#P3Uvo3PRH^2$Rs>YEDSiX5Rvp0E>~A%tV(59KEMDPnXKz2N8sQEl#L zR#4~+)w^bci0o2+2>{mo?oZ+FG8gAiXdgt}!cnIYR8nJ;CAR$W)bD=Daqd*~JDl3g z91`&>@d&kZ9{#Aegz`)wf?x6dfHc}5YU3&i%%T`kX@1wsD9v$sxKW&-SuxTh&Lu$cD!l}Y3q988Dd%-(F zz(&$|H0#c1TS{P_PBRr6QYBK2QVyJs0<_$gb+krL$Yk8MC2rEd_uQKekKtV202c|* zYvd1SXeAHpI;SpjQ2}-?>SYy1(WiyIX82Uv3cH!oS@D%e#UIkd@C0||ByaT<@4VIh zW_=(s;Pmy*IGWR_WWu=7Q$}dpko$WBST}TJzf(4Pv6C(ji9#pRhz%hR8bjP>w8?Jk z!f}S?WfNsOIole*h&DUKO7>r=0>BzcUe#A-Ma9TG5g>+MzZupqT$UCuJC+P{5|(}% zlRvMF%5zpxc;^(xRF{c3@#_JKZysg!datLzQr)d;6MvJy5`^eH_m`T==3E(%Q)f7Q?sm^S|z*0oA-a%arh~?5={s^-rj1UxgE{nUu6|V7gox;(YEMp(~J?&E#HTLoSoG_!_ckGlzu! zR1RVqMmQ4gFJN0oN*=Sb$?C(kOe2ezh$SE#AL^} zvSSn+uQC7_NOme-Tn$iW65AHIF zTX=xtonGgps1wg`zl|~wWLUa|HsUtX5nBQ&veJ7&l1q^rM{PQ@-bqrtMj8@LD`J%S zab6X7mJ)9kmN=6$vwFwKlBp1A8XFT7jSviD&Bz{ONa}l;5OuZGct}VM`Rd3q&}e*V zB+}UAQW1og<0mjH4MLxitr?Q6Y)XTtfO(KSH=i`CE5%_AA8W(-NOEZsiQnWm$#2g^ zDFKI-m@vNQjBY6kjmF59S)t3?OwKmpjc-pfxf~YAkFEc3#QE4aOe(5a7&NXo0$Vzb zJT3UUDpW<%88KWp^2}I_oFLpipk0&37OMjg3JOsuxg@?39?=duB^NZ5b zqOlww(QYOc){3&!`|O)Q@L{8R#uapzs(cFKX1%uw{wi^S8gJVZ#{OW?7MT?tSrQ_T z1r`|&t&l6c zoHjPr|IHH_Yy;zfZ8Vq`k@x&E7^O5c`Gn_K)ysFCBDwAnurmeqXaJ*kx+ldq=>Qco zsbqJYaq@~}6uQnJF+B$3MnYZq#3NbPTF0gpbafW`vCOyC~!&E4iR6-S%G?1sE$&n0Wm8?^! zYZ)8-FsUa`l-+KJwdk)Mqdmg*fFxB`ZPLk83suGK$uDV=O#4K@ImYlyiFFpq`=Y(2 z!aB3Q5DrMLeo!*2%KG^~f}+hz6!{n;GD z3)`vP=ee~0nKREh#HD#00o;E5U^epaej{jlEj9Wfce&{Z4{tCA5v0?e*EI5r(&_zH zOJcO$d1bpK1a4dYBq9the$CUR7elqpjvK%qui2aVcS~)B^p65th|iIv3BVX(D*&^~ z|NV{tA^2c-Wr%U2$0e#%oZXF_TC{cCiaENBdG zpGEYMw&Q>_Qdjk$w|J8vR1RDkTRb{_JzIkKbbBd6{{Yy+Co%;-W~?}F%v;PsC?q2Y z?u$`zN=qTk=_2CDd~!&YEqSyLpJI36n6JZ39EB96!Y{9^?Ga{9a}7O`FRH01$c6&w znjv0Ey8ZWlhq>UOHXwgTRR5STLhaR7tmWi|Bg&Et0p!4+Ot@|DeXQitaK$Aouf#wC z%(27P6x(N@&j`z{(RtE>#RT;AXv8BGWY!l??mYhbB-G<2DoKCiuj_7t{7B%*u9DE^ z!OkP!M`pe2MPN^8^CTxgUT1K=tQKbeiwv?&zWkdoR6KWNyjN!)R99jebX1C>$7Ar!|ML_uVWbMnEUyir4W+bJp{7F4aH= zieT&Om%SZa_};v~UFHo=P|zk)uCjjVI}!P<{OxRu7}bwMom9?>L%xWKs5a~4L94?E zDyb${jTmu)Thr-Z_!0$vCYX>#?xy`u2ukyY3aWJsWYxS7f%JG|``3tVYo~NQ;W?L@ zg*CYmzgDPLU0DGsC{H1lqCdz$#;N>gORqgQ28KDCg4i*_5F?oHNW(yH3esMKZks_+ z!y1B_$@OWn{W#0|l&QS2)od%{W(+pSwW-^V^@bp)e0%2|u_43G(2b zqb?=hiyBHrk|JH+yp~i75(DJ{AWKy|U+X$3+A$?ANJB!1T#*N>@T~J4;W%N;4ycur z5CCQcHo$ZjO)>r3z}E1E(0gF?lgT`V!{6^jJsL-^IWxtE z4Jm33wcL&x#K7R$m+BV#RUdL(m(AJRrzFVMn$Q%oWO=C{!QGq8o)pWBcIzer$Hi88 z7pvke?`{2Zvxd@t&!5jD-u1X_4|PVPNM4(~W$b2PhH=6Cp zdYd0u_xBDn#En>*=-7o{ylc@>JnBUULX&`Fae%yH!x)}HzxUym$anQ83@M#-V_(JC z=7P3U!6|3hEfED{nNu*aqHL~VaZWN%Upzep_!~CaVFPH;! zGE@Nag*qW24eE$DA&wCa{dhV0%NKcdI6y;~(P2HI?S({Gp06|H=A4;7q{i_JT%hkf z8j$iMtv0}o%*zcM2#13M@{K56c+z!oe%~rp$L=CZx!6@s3xFgY)YA)u+a3$JqQA)^ zkwco|j{FQo1uvB2-YXEIENP1C-<_ZyEQnBm>VkS3^lO(XDhgk<#0Qm8nOgKcR~KU=kgM)$m{Djt*Z@y(xI)@PYZMxayM+=7k!te)}oB!U}l% zEu$4*6%h#FwngaRu`}slvLc)XnGS7C^$#BK@*fI}(@%C|aIneZsUlefV<7imKV1?? zUHLkDcT^QHUw`hgmnAX)-c)}<&LD_w&2mP<>zZJKAAV|zT_ta{i*34@Q@-o3C0IJ) zLYhR1-!d?O{o7Cz|3a_HV)YtL$RNILF9In3>;h&At8}Qf>Hl4{b(M@F020aPu9jgd z2-H(X*H%eJzk?;Y*O2+e8`)T=&!o?PwlPqswrZI+L5WU@>&7MxRb=Tl4T8>$oi zp-^Aqh^rwi{_79S{Llo9lji_YZTtylfIm;f<8%Z`JYmo1x_9y== zCK(;4920d6hnRMPS6BF67&}3Pmsg9OxvHjD8t^ycLpS4#nBQWwGs(=&A8%?s7lzZ% zmm|@#5G|*jAtuPDPC7wOA7J|q#aebat%TfZ_VFz2?-d<{uqLO-8al}PdDi`D?H6#UEsblZVZZ! z)uTFqwTdGMAaMPTIr}LVW>;oU!GHWPBvg#g4MRVquX=9N-lQ4n(;PP8g*G7*@PP&DK{`>sK&n=bFNd@B zJRs&aPsz3lp*r>l8=No7c(jtp6RYjpWTgW-B7B^NX}5^`Exfk?6reOWEFzpKi>)r; z&Minh2d{xopUH*SXR zf4Y4X(10fV1#%F$!xi(nV+IT!Rr|(YUA`7`X${s)%i>Ig8u2$unW3jn9;C@|ZFV`| zo_WgF{Tv+&kxXDZmLf8I=_3QyI_z$&&fEU}(Rm5&e0@uN5!B!4sVS+Q?4GR2iDh+Z z6Iao@seg+@rnYKgM6t{k-^Gf{cRb%2h|tOtR0uWlgN@bM|1ly7L1w(V>99tUgez4j z!2S*E(%qr!@1L%kHV^e*7uBx03+6Mnl|^RHhQ#iSbqpeRO|9ZSc9Er;C1TL(is3$A zj7E*&pg;H8w>GR17xs-f4iD!#3ABGIPAt2*9rC*;Iz9tQ16Dov{tY4m!ctUQrb`VR z9+}IiUNK+3$)qpHJP(OGJ-`Bd0{;|Ml$@1~i8_1e15-E;zyrsyei!iU?DxJKVhRa2%Zu9Gs3Bobm%B;kEiTqsNPp&$V25-pB!>UZNYj_5n;w)b>Jew*b z+l)y`PGyL=l(Ot&2c4RK7b2p#;=dB!hfn`3eCBQn4VKOkClgKSRB&ZtVylNCL&NrV zm)nD_YlyDn{N)^|Q=QCs+wkXha0LoOGW3Pn*irW(`J?$oCNk|L@d>2v(TKo1_w0@A zQYXtdCebnA$ZQNu47KlAkk|>F*@9ifr3~rc*}k_&7RF zU(@62l0~f_Jb+MpGKnVqb~}JZx5rVL!z+YJz)@YA>`W&d2sRKNU*xa8=^4)GJMg;@ zcgdgUUJjGan!~b570}5$Zn@bXyMJ$O!CS{Jt3gJE{uEHRXp7}gUcpa=vshFKAy(K} zoC;fBr!(&y1WE4Ye#mDIxO|B!6eKr7;<5NDt>j?6ijJR=gD-qpK|ATT*LYa3mSY9p z+#0av>)XObGN3dg9Z+Pd*t|JVLfcZZ2kO=}Ptx^fV&G#FJ+Hgj;7HVSjTZTtk8p-+G`kF#&Mf4?R`-65nenfS+twQ%_q~i^ObP1R zQ(tsC_Nyur*kK1G`<|LAn{D3!EVLQTUjWJaKC zuHT8j7;spr%3s;a&Z1DG&nGkzdQCbzF|cxGJQaKfb#G8 zI!L=-QQN3I3IFv`k%3$rS5BQmE6H5BLIbjAE=mz{b;)g^xi*Kmt=$mA%Zf>mP`-;3 z7xMD9d4*)nN(PnxZcQ(mzP`PL$9B%&i-=J$(cY%Ml~BAu7T}L^TMYX^`uX5Y30C z&t^A^D5JuV%GazQcr2K-m^@L#gH5pLtJPc^nkb(r&!+zB0>n+g7mxB?8fy-WAUEp{ zwLhrrm$?w2*61JiiHrn{vY&~sd-dFki5VqDSA!b2i1k3Q#wNF{#ixksZla@nZJx}-&^a3H zsc24nr?-@Gn2-X6t(mkc(ptr-NRbu>*$UgBTBKmfA6$An?OATVdXTQ%y!a*&rrAgB+(g|Xm#&v0(B!-l)6lhh1A!)M3-vwNxVGE_NCB|NbKC6}^<$W!!>LES| z8jPCVL&YPl7@1WssV|=o&&{RUt zpg@NZV9Uz&K3ExgPMPCaXT6$V9&u90f zPTNJKZ+~`Yu?VdREsr-^hmurK0A48t_VI*95VcdNw+m&3butb z_8bT#G7ETZeDv&-&)UMqhcT>~{hz@=+#^JKZF#|N0)6!;oR@f>0|G@sx% zYtwp%@HXdP7S5h*e9mV%iFThKey8I^>G3@Aze*s&ETVPK%T>(TY^ zg3|nl+~HP4^-%4y>cO_I4w6LCdGl+T07EVQ{QPpl-0GFckWp$%15q>Jv`7pS7}4t5 zqJwghtp!*&==WS_gs7Yy8y-|8t6~F(3(;@lDGJnoZz?2Bgb73tJtf|^Z%`TN<$^T7 zB0A|3(Rj2_-p1QQ06-AYH$pQ^l&6SQ0)#s}iJqXKyr(Z6Bvz^0$dtc$v0lxwx3F#S zC8JR135a!g1R-stWdDBS+V;B6HLP~P(GRqwzy_cXOwGB)Pw)R4FBa5QMITJ*Y5ev% z@7VY=zqZm03uMQ}u?B3pJd^rXHCIL?$pWI4j8w_{zIjZEj@mqG=+-tLYYYbvBVwjR z{@8hq=X9`5VepT@nQ+HRv$S<}T-_=bi(x_qPFA-}R`UQOqcnO%$JI1l+*YuG-XZ*% z6U)fsc?pcT^w8-`yI#!_uKegp^%AwmpdFQJVnBpX^2!R~`4ueZ1MW zgP9AB7Av`{94joGAc=vB7Lg6MchlnM{e@ST8p25rS2#rVl_aX<52_2>MJTN*YWCda zw`bb2p;o<|&OcLa`Eyr#^a*ONv0Uw?8IUtL9y&P!ybMZD+d+h321Pni(9(z*(C|zA z?#zrI4KEW)jsz1vjit(n1CRwtRDATDkI6S)CXG7{6%4;-til_KQ4{YZ=X_KDvM;k( zJfub4f5%of>-ROAc$Z2T%`s##B)qb|`?uQSSbGSQGxsq;J-~^DM;ssFYs%(+>B{%3#xoNS~-eOhKcT?lPZdzfmf2)m56!IS+uEjbHw7s^9dNuRUW9O-B=tlcQ_p*)kcG(e zRoIC2hYXdT_x;L(w(D!blgvn}@)6`iqK;*sKdqT+Mw;Wv!OS41tbsZptmFYT6tY{g zQ&qRue9uXz!X&jne@W!urk|=BrX9dA?1$s>dwyaOy4=6TG~*h$9BSm3g2C9OAW#3)%fR!Eoe2IJN(?Bsv)xHJ zKvXlO-HpUIS2g2Nm*(76bc^P$%CfW? z7e}6O;F;g%{WEqu=G^NiYuLS;+!v*4@dXh(os#8FlbpuctJ7)!*@edaNI0}kkL#C; zo#7Ew7c9z)w-=cQE92smBN@r7l!oQM#U5>DGCw%dtm*dej_#By$wK5&QpJxCv&YJb zhJ^8N)f|!KJos+}vgGB_t8C_9&WfZ?q8WKp!a44~N}3K&ts2&3>`0D(D(7=5PhB@R z3Z7f2yDe|j;U|4~>ADvIfUKt$L=O^T62!1d^|l%#=42+XS@Dv>%iu+PTZGgvO2s%& z7sQ!@f($T$uXurlAxXgLpt4hBBmPnfA`V!i;af zd`-T=bdW?xc=OiG!LX=HZ;(&@MnF6gP%pie=j6H&hNM{8QsZyp}xMv*1 zY?+#tO5c>S?%NHbfgn?LB>!zq`<|&x8@j_;QPN}yR}iZYhF;T;(XhmSJd?E?IoX?Y zOgLH=y{zq|LiQGp)#W0@IoBPi4$5`@qh5=Kta)Nk)F(Ba^9z{D8Y`Ct^@CAAYhqvr za;wRHXyNY^QSMGNP%9*18{M8t2K@du%=cN61CvGK9QY|5vPbs8{acS(^zhi#eL@K8 zx0}3y&z?MfLvk{SL2=0k9I~!|<`U=ehbe>`pDkZ2ikmU6fg$C>{zQ~1b&hL%cnLS#R0^|f?2pE2#u!#*7QY7wLU=1IYHN!L* z{&w9f1G16>7<~qb+-npGV+YW72$~rBEoWpm!YpcnF$Oenfk+2aFQ;YGQ;JK4W{!Sq z{=-JyhFnTWOL^7sU)zlvP1RUML2yVgbh~)HbpQ%%er=c7&MZM<|TtFY_!WR3B zt?0B|$)ungkrDA^Vc^drh4M#pcU8jMp-y-aE2Yx7!%ZS=cq)5tmxBgqS>rmsFXBWM zf|{Nkt(ha8C{$gEET(q4#LyTVAqiBhwrxt+)2>$1=K0i+tCk+$9g zEPvs6%Pj+0-y!}g1dJKS{H|T0KAEp+_&9+6RnpvQNiU9({z`-o9|u(`N?%zjM& z{1$Xbuq-j3#O$m3rHM>TRH`0f^ds#5C2wI8DJzG#16m&~5a8=e|$yu$>^5)6>xVPAe1B z@k_t3qX~*g=YasTJFtrI5WOL1mO?QvY?_*llp_UI=A;%rJ?nA0y}gehq-OF9&Xi1MUymx5xaNf$yC%*34ELGO5C1Gg_`-J;Fw&_HW5@A@aHjgt>3V5fAHsO2vLaFtnc zi3R0jr@bWc`BEXBZa*$dEa2@X?_kKSk>>ws0qg_b9uM%`35%4W0Bou1i!<(JeauK> zwj8c2+OM@+TEU-a9PFB6r?CB`N1PZA@4aRK3_MdI|qxCaBwJ$IA{V@K) z^vZ5j3h1pVEWtqI=6&ubB0bmHaQp00u!v+MmB`)1+>4V55d=+LTqvRA^{H^=0Gd&O zftH~G9nohbNe>lP){cb0WrU?){2z8Ma>q+35c;B1J@rNI;gwU{DMY?VN{buN+RK0J z$Xxi-eSdnRve7oEOfCH3W}~$-UWf_k-L(GG#w?nPv0Td5Wa{8H4-q(JzDTg}7J4zE z&nm;?(00)a(K2rWl}{(K0VCJ_z(E&-%g!SNTSrq|rO(xkvkM!wxbr<#9!gb1t4@q$ zaAs&T<3Vgx(U$)`urbf|lT|RR%mD5!kxTsua$LzV;a9t_li059+$T$xQ@E7!F|c2# zi7Gbuc0#m2Kg_0I^u*Fn_miw;SHq!x5GA_ud=c1qoQh`cei$HqyR0b7-6kqwO5r%7 ze0E=Ek_OnU7l2KaP>9f|uM19$1?%=sJp83rh@=Q%4Zh$ztJ7(*kVFCw;epfjHKIt` zJ3`GrMSN;*URre6ts5>qiw1B2Lo6yF?!OT%Q)#jg^c(-RX-a%DB57=lHCNoIl|K(G zu6o6Waw+v(k1+<^X$Rabo5dUcbP_AWLWC2ID?jDvqTUUod()XE^o^8lLxPtB*E zcu9aspou6IOtoXd8X?Hd0N?vjfBFRp>-7vjJF{^H!t0nl0CM(V)1GRZaOjNq-n9tx zLWgRF$<|P%EZ~yj$!5{1Gr#^3LusdENyIu=JM?s=Q6=XVDy5xd;=QfH^$N00lQtwI z1Wo961})%WK-?hWl=-di-2msC6H}V3g)>g2nXS&k;fc+_sdpsWE%5yNFav`Z7{R&5N{0iCt$+Jn zkM$a^@?H{C#~hE{$T;@+4+p`~MM47{^>vPx6a?30w%0ERAqI9cE7|pbY=DGoEY!@> zY5!aQ_dO&)baeM%IvX}{YxXJ`wi(}R-(QUMW5eZ#L;L35W}A(^N#e{$rgleS5f!z0 zA);S3MT3U~4~Pv2I!NQ^Fjaf#9H8uvFNz{6f?-o3>8u( z(o4dUINw}&2xw_4L}|rG6*dH3l>JLp2N{m;W1f(HpXBs#g?s+OqIRS@$GZhzunqK{ z+<6P;8n7(Q%&30u5e+e-NiVbEPaG@Xld?zyg5H|%FgU|~Q7%?pJ_VqLA_yPg|(B=PmO* zu`U)ViP(4fMmbVW+?sydwU1pwSs&wDQZwkb$-YSz4Skp7+&hSUaB36k5xK<$mbtDg zbVvYK%ahFJllW{wS|=R%)UifD@C&b$8bd5Bz;oJlF*e%@A}6(Ssn>mDF{o^AGB z8Uaof0*R+xyWhRPfjtN^CKb)R;L8wS%q-u_BfJ8!fD=a617F!B+0=LL9oPH^TKAA0 zr~&LzS@K|c6t-{}z)zyAyS{R@#)@jEwY)f&ZUVNDZrG4tFGZhf!FuPg7kR_i;LG&I z8VUHa8JCNhVkF@B_+)jdO5}-M_bTcD=Cy>??`Y%gY@@|-BaUm3OX+jB?<3>q$FX{| zFg=%P@I3{oj8d&y`CgPBGJ<{0v0BaIu)qm_<}A`m@lw-A$BLN-!e64?dFn%)LUJx2 z&#Zp$G1Gz$EI5EmNd~Gr@d}~x54KZS7+Qr+Zjg2^?a;6uzYn@!ORHM}UXy%TD;HOK zju%QRK1WjUjr7_HWF~AHaqyJo80vLj#{~eY57yh-BJ0({1zmdb)hp)rC$s8_;|t60 zPL`61|AA&O_cOhhn`+*5XcCL7K!MAVMV2PNvrptl&05uhACvi)2QcJx3KKENW+fBg zNsl&n=huo-UB-8$sMQ!O9yzjzDU-%ar?XeA+{OK1@#`73)bQjArM#EeH7qB{D76m6 zJHUliZDV<_G6l05Kl@eAPT;&i@_s-XA-B2bWdu@w&CV9$=X<*)znR(j}|o!H*B z@8Jwl3PJHp*P%&wzL=N9bj}|0uftOO zPqpjVeny3jWs8O*0xuq?pOaN9(K1H$pE;(?R zQzzQBY`YoHFt3BG{$D}2xzGK%VWTAqfQ*VtBH$hqlkQ*xer@T@u!CV?I&8=7sef;@ z1QR%PTH$qX?ta+!@zE#2Nopqi@{uo7@1MxSGF*faY~VN^>#5@zY+#1^#oM|4+gk+5 z?9rVOZ)*TAt#xK8-SO#YGGr)TfwR}X`{*~=PO+XP>-xgHPZh{ORV;Jd6!-t4=^TSA z>$)v`VxyyuZFG{3ZFFqgwmY_M+qP|XY}>Z(dB3Xre^;G4yVhQ7t~tjT&+)c9SyPv6 zq?ijr3;c(|{cc+^*rOr<8Mmw2)dkQS{IXaU^N!QJO{Qp>GzNQ?U^mXx;&q6#{W5Ip z^S%yDv56UB0;Ud|Or`&4i1T?(@lij8(8A2IuEYNK7YmSy!Si1svc6^p=GzCO10jvj zLY--i$Aj|*(DShPogf3%*}yncXle5HIbaLc`-0pe5`-G+^y0Z?#gX%7F?K)s&sozj zK}sc;W!M7>&}hA4V$TfHyKqwbua9LTJXPwfT1*-cd((EmDn)(JVB!VSGvat_2>2IF zQ!8md&SXof4~+MyFQr9Rp^zRp-Q9n_-tu2}l6b}Q>Y93L*Tv<903`Umn9zv*(zc<4 zDmg7vf>L1ZjBU7Mgi)B4-|DK-n&p=Hgw*9lnBL(p4BCmnA(dyXgA4o&}ly)0#?og>o zM-xciFTk&!mT>-Z`mHuz&80zLpF+ItM`lD^2T&zX_<)ax$YXwvV=rb0#|fw6>XuJt zEPAs`ITSqB41ut<3-p%2;oap3(UE6tp+{9OLV)OkK*- z^ogSZ+S4Ukojj0&U`^0?2kd4EfpI3bfUmXX?%~z10D=!++kdITic%+mv_!8TpU)S| zER8}&Ujf${yt91D3JG>wIXOX*ZlbevIVlQaf{?^D{SVC%9NsJ@JD~{I|WIJug6$@byZaVP8oZ6V~Ju8U4ulJn!9QCd#g&@<0_-!Pv%X zi{tp=-6jx$`q#F{dd{Z-E%nItA&`2kVaaGVoz^JA)HBqR%k(#k`xhy(1kFs)Io;g> zM!`>jptI#c*>7)hW4}Y-KIC=X^*@spT;W`cR*xnnuk>{3 z&sa^C3(kq(chu39el^6YLKVh595Bl-MuI`DEL5)?u<51LJ`TMIgS9 zU6zx@iyD}J2b#0aTeg}v8cOv=0A9F;!Nq}vlPR*?<%)|L;LTN^S6zm{Gwuv8iro+k zVAZu>&8F*vZq3GdG2~ET zuNY_iZa=p!MlEuDGw*F>Wv;-G86XnL z_R742oq^nBS`wh;>yH-JK+v__R=iwIqlH&bYu#@B+z}Ffr}aDh|;cPQMj+URClqkkBXPg!KpEj$!X|Bm2YtwI7Ft zvLc#+b_veIXbtsdD(XR(BunSqniU{7$(Impf!18QyW}fTmoI7XZbGi0-(&d#^nD$1 z;s8buSlp}q!f2K84uv0;aiN8>1qI2|4~M#}y=I~;RYCM7F5eAM(wdzXhfh7U_TssF z;15`+lh=%1`r%fw;p>C6*6W#sVD$x7KjSR?^HO>_5gZ0qu7DFqnZ-cysbGMe6;avL z4sLV)$`J7pUj8h-R`+kzR8(#&EdFs{+xFO=V94~iQpbT|>FQk^fZd2pG7D_!Sv5+mx{xf@r-~~nCxKsd7f~o;=AfO8iA2X7Nh`Q^u2eT= zL+j~c3^O^t;!$WqH2&C^Oq)pyuJZ(;j&DFd5Mf8+#L2`Etel83M;d*vC_w|d%^M?NRr_!};$lNbJm-UY z@ZwF|@&j3Iw%Yzq4;o_a=g%CJN7qmW_O>IGB<@y_v&en(Q#`g}U8C@0f z5&y|qy|T1|;*z;Z+h3uq-<(6d*i(KR{sr`y7{rq7boIN4|IgV;!bP&oqs$Gl82X`$ z2FY_s?Oslw!!dmVjHX+%R(F$9MQFjN%!bdQ0nu^&vnJBaO_Nxzz25R`H_}4wwKSSr zRd)l0gG;dhE z=qgO~@9X{m22t_DyLGE%zU{`7AH>8&>Sfo*ZnQ4@27i>RWD~i49^NEvJ8vo%4fGpn zjvS#;i$sYMzyRXMcYxSQjHKEzcqSedZaQjx72>R+^;IWCz;AA-xi+bMN5_FIU1X*U zgk3d`J8D&wX@LcQHTWkG7J3_Z?Hp3gJ{tGS=;~m$prJbdE3d7{S^MWQTzW0c?HG!X zL#hSiFAB5iA1fm4SlAq-jqo(YpapTXwEj8Lb95R>l3%kqUYFHSd_(x*RfqiU?0%zU zgL1}p9=H)83S7b_!LgyjT0slzZ^42NeOgYxn7z6G?sB`MT@43!*G1%w_Y!Ifzp3xr zu(v7bC}gQxCMn|rZ&qDed+!(7)QZui!Pt65l_sntt zp18C+>!_pLFuw4@H@rF_`o3E51bXhTX|kt!TwPUt(b zN5K?nB0l&&XAN-pv5v1pX^YK!qzpyKi2h^*zhsi7-;Q5V)3ms!Tv9C=%k6n=fwR`J zEDpSVSh_zNH$4J&0po8zK99N&I@K4J-kHm;7u3y>1=1m(NtGM5)xPx0e&N4JQ zQZ~2TRc1*HF)5pO2@_=AF&=G@p{Mo|M#!-{D8x1SfXKBF?FU}LIew;Qrx!xs7mY1Z zCXM+eHN^e!8pJHQXai4#N1tqc;P1>dB9$j5>U5xuM27QMprn=pcxDgv=E#N3Y^Nz% z+MWIC4-=+N6x;-gwSd8p>EpYYVlwM8Rk(oQ>S>_VFH*nCT69>p5KRA!6W|-mq43R= zpkvQ@qDG2)GZ}lgabgue!9km5pag3Lfv2ez1pY#v`Nq*_ULE583AtqcvY(bnJf*H0 zl!XQZ^u#!ju>T}^#LSd{XT%O00^|L=YA%e;87Dt$G=CpNyJt2Kkbwim6US;yJi($= zRJsP|VPUXWP(bh-JD{PdeUl|rkg%a%=U?WXmO@okF!qbU2u#t87F;hkVxs(temc|_ zEJ*c0ghhb!c`{fwlr4lngWG*}a-Iy@f5-isRUVi-xQa7GJVklEK5^ zbJ$|pib2jUm;9qxJ$a^Du)Ro6N|WrY;{@KL_8YR}_T#pIL?{3dD_)F9z`@C@-q7Xg zATV!dXt&5cwmENthl`oxm(INmG{s&fzErTR!K&+H~Es8KRa@Y|GUM5 z!QcfgR^-Od6}?eRjFESQp}A0K$Ly?_#rSoQ>Llq!Kbdg`B2fP0#kOCUtxZ#AS592l zTXQ?L9UE|w>6!+S2peYW?QBV*NS+@{hDf-jEgA!Hp_V2Zbh2A-#B9dHeKbpyw^O3}Fr-pAiiwmSJ*?-*`KI9eJ`8T^>gM;&r_2}t0i4t^V_Qsr z=V)a{h->dI-nxk1Oi76s{?1vJ9CsAz!UC8aeE{_4lm)(32B<|OJ?b@VP1QirN3{RC zW=IcjcRN6Tb~=0pPh6I{GGoJw0#2^XfHxzrO#|72;H0odsV{ps7Ifm%Ete1t+s>D- z^An1AQ`3Th!g)}BSxFMWO4ZR}rUAkPR1))Y9B$d3Z^U%2>oarov>Tw;yaBuMOpckC z>3F-ri)pQ5(c)}sSckh@F!v?YB~MsD$F$5pD{w0N|IRMVgK?N)>;P``!(Sjk?77QXhI6E;cP=Q ztK4fvy!*SbvVrmcPO_|a(YCIY{|@A zABJTdQ$zn86jg-rX5h&HMvWDQIo}g}ZbH^Pxl<-eyb?H_`<`V4<$1I(VcTeB=N-S~ zHl?AAmTMhqtF(Pw8Yf#y&Wx#ksty5N(m?*j6jq$dG#ed4EBx!82zfno)dFa58eGcP zevmGNAW_YBOn-Ht*s#|)Cta95w_pL|^u!LH*)pgDVeA+y1hqcNo0Pwk@IWK>ylJ3J zRi|VH|F<1n&aE94k^CqE(x%LxAaE=8DY6x;KBtn$-)L`{Q!h}jgIfhlBzXez4tSpMtYuSTzoI8-MO^-3X_CX2Dwus zt5O^J@hZ8#tWiILRO}D0=4PxP@54&% zz*o)ffwk(yn>PJ6Q14;l1%Nc5VT5C*pA;nl3=Z_EQ8N6Llz!6vdiD7*vNC|3J9D!0P87f_GDuWukkMVihtEG4uxc;1n65|uxhmJS_y?cO*REfiEQD3mar$Z@> z!SnCVyHDPo9h2fko>(}Y;FeaVuWD=MJb~;$0BBWX&163NOo8}L7Ag~ULHP^RY%jDV z6`P0W%Vg0drz1dtX|h^lW~FX4oHvVP=_5pdOOPk2!>3`@@Dr5%NwM8ue68fiJhM{w zja|mO{p<5;`;7yT5)BA+qeSrdjDFM>&z9A}AH3}2*V{kheLi?8#%rM;=(KD#A1zE& zr*Vbgh^)XAjj9|_F;3dt>(1@qS|`HG%%(lcT3s`-v~Ap-`3#y+?Py;n`5fuBo8=+H zkJ?yjp~@Ej$9$P~&uLX*PNW5l(_PIx&hd_mfhms6Q$_cItYvYXB)tQ1wQpKzvG>xH zWcqU7y<6O3rH~-nc$RQNgiL>&ELatpdn}zZVLwn8-94-}I3hH^yV`>IQXpP@P3;<%W(bxsXwo(gyX_8& zsm^41JRp6G+cGl>D;1gqXs>XQ&&EP9zJ$9ch;(plUN8V>ufY@AAwn~0(O^P)^qLcO z5&`f*+ta|a_ArQx77$Syz8 zialShG&h4?YR|8YCLdQE%YgL08cr$zFw>NsWi1HfrfR4mFJ0O|SP!6aTBa3%-vPYd zw)Qi?mwC_JE>VeoY4Sbc$RICO*zxp0{=ret?cnCN;l0_`V6cHSer$Q$<9Qb9dVW0Z z4~IBD(Hv0>{fmgw8)F6PIt-afv0$*&maIvI@G7ggD`zkPzKogJ6Omt3Na+z`Yo`sYUtoJAq(ckOE|Ka!t8=u7)d>@8v&No83W+F#qL(89! z1&s9IWycYq+GeNG)|cs$C}Sotnmbd4(0Gu$wrhNG!Ev2}aqV{f?!M%5#d1#X+m%l* zj`8#aO}XH9yc?|1^`_Gtr(d96M9uMl?rGYuQ)?D2RYjs4x_h@&4L70Ei_SlSAon)Z`I(+fbzxhN7z_9yk?KY%j z!-~b@``>6s$*-w5A%NZ6#}E^RY+(NF?O2)5>(@6Dz-VlG>@-bv%Pqr+%NRIhsQ;d8 zE7uW5hzRxMPHSxh=JgNjf%@i_n7!>M9^)9!!@DDj5i6piTDfOurFq={vjn~4-c9H7 za*nmNMW8{rb^01NT2y!>Q4$-HK5vb)&rcxG6ewWx0BIDSa`>0-5+uVA0prMEzRzy(ou zvU4G-VLplhnaTolEctsP)MBZS5!sf+k$50Ev?|DI$EDNH2@e?3@~e5XEGcEtp_ZrE zR)NtU3JsHjg}rF2aZk(m5H{Kl?7@G!P{F8x9*q^}0>D52aB4YVd9u3Q{TCkqS#pfO zwP_lgREHe}kAUxfeWML^CG>}*@Q18`yooa?zpduHPxsST@)lx~0ueW&P@bRS;~qw5 zN+U)y=Hi7@0jAo-E(1n(F``&nWkRtBokNhQ2BNQl&EFjb*Tlzo9cRv1B3T0=$FvkR zky69}6`BL2CuauiM+*pmnG#uj(geEj)A+$gMnRq<8*H_6x0!*S0VEDQgx%rXkYUNR z!xfNnQxx)$a4yybqrp|-1d zc@G6iZqi(B%!_8F)#h|NIfk0qL5gBw{dXR?Hk3KPv6JDl<&v_Bfrf%prT%!tP1EVB zODk?qhnw;R03ce5;T(2(h^btPlZw8xSbCMg7Y1-acWLJWgoSwnbuE43AA^Tqa?Hh{ ze>HKP3snMvgK4hagcE|*tpPLOkcX?;=Jk<2J4U2HV^jLjXy_uKFJy72sT*X#%yzN1{+Y%wUTs~L zMbaHA=}t|qMQvs0%tDU)%x&Gx`$Tnkt1_gk(ZBkGY7O0I+O^B+ zAbQ@Mk?DN;wd)K7T_KwC6tZXy0&MD!meP^(1y_>4Kzn(f0@X_(eK`}RvD7bb zT*eaSQY|W|*;=@kU6@u5ducc_opoZ zQB#>G*OZTxfpyEQO3hA|kHN`X2mrhJk!ID@Nq}BewLNp=-axQweT4$B+hAs|bVN8) z^Svu1;LdvT7^S>KGLZej_oL{al8RYUlPy0{W+1sEU|p8e7vKj63SPT?TPp4#Pr4+22vJ|_ERRC9Um65gxlszGtpt7v8`i4Xf(3B zUum(eyTXTLtANOWM>~5#ft+lni3%b-@tFRio6Xdnm(`6LGprpn&yyFlh~3IU+06xV zND4QMT_fxr34amvVq?cQsma~U$zpnt6!nCa9ob@u3g^Fau|;JQ25Kc^GxLnN2K6UT5=Fk635VvXM~^ zPnbZJNbS1$OrE9RMe_HdXk@L3TXCm6(L6X*$-tK>(#GuB&JnQ#a>B92M3WNBdrWeD z^^ibwE_q8p!Ip?37l0cBhtJFsOI?Bfot$QJnE8=o=g4gSB2HDzgKyE}eR}7;2Tmq} zbI_#t^WP>cAgNEjVG%1zbSATBWn;}_;&OwXJ2I6X;$8M-vty_+elaC*7${QkKBln! zdN^bQ0X&ET^jIpe$qE!Kp#x&V!$1Io`6r{y1{HjNU5vS~qCokH(=kpVfTj1hEN#_Q zjF>`qi<{FTOv_*`I!)cfgq1_Wz?GTwX>BH&u)qzYzjPb!97PO--!TR4x(2RB)42I; zfSh)OMFDv&GNkYK*CK__=Xwi$j8@c)M{l~jwdKE4{LbfS-JZWvf-%{dch1jGyXxGm zOzuY?Lr)s(BX1f273bqfZIYJ=zWQ$!Ryx{Rlqx37?axD`=U2L1LyIZYB$u${t)^Zq-To!|`>`<@fu&m?Ax{ zbnKTS@x6y1bIeT@&tPH9esSHR_|1|hDpcV*S#_Xc;U~z;go6%_TN9s1f3^bxRIMYk@q|r7#AZ&9~jIth+oo!Rx0=4(bP1rS2hZ5H&uV4 zASjAAqrr2YipnPo04j5S4(vvx)O;K)dVRn4Pn5d0_UVE91_k}2_kNl$B3t!hpb{ZD znnw0f$u4EVF6e>;V4r0*A%;w43J@8qzdWaGiLN8CDAUQv4`eN~;Ol60GMv4QRlH|= zeP*RR;XwlaQ~BnJUu^4oaP_Nm%a6m%WZor_pfBaVCNWv#J$4LhKZSE@{tYOoYoi)2 zWyf~%#j)RAeFydZHO62OWkTFq*rEageI7aQz3o`jxz*z>XZhQRj@-y z*i8_=PCpRCKgOk3g^da_MFA^uTZ;153axyx3Z(Iswv0!OC-x$8zkhwvaM`L%11+rPc;T*gFGk9lU04m2(+PBGMIcNARG7^ zt1n6vz)BV^NO;<~xL7?{^iUb6ot<9r^I>MmKmlx0Xln8lVgF$Ll#IW7GFz}9HPnb2 z)}O)6WFSyCSzT<-bomqDtmk4g$kyVr>G4qCQiq)tOid}>ZpQEQnoq-DM{3M`;;!=7 zj@(IsEW(j9eN-@uiq%H3%dv4Q{VpR5Z$>V0j)&uQ_4^HdWYmln z<6Mo+caty+&l;GGVcUJ_(xuZ|N_;{sz`BEH$7DP<)Nw`r|4&T(M`XkZP)#2 z%yvNIMCI8OsBi2WJV5_DQKCqpLP=bsvKXESYo94kj0?|MlBbTnLFH$k!^)z9kFnE5 z@cK3JVHz&Z>E+G@3e|Ay17#U67{DV8*gP~1eKd6=apm++cNdzFLq80j`k073&Zxctj{yn-cZvgK{WjbQ6o(~FaTCZ;eJob+WeYv zM8T)?x-XJ_pXWByn~Ol*kyOR8Yww6U&e++-@h&J>I)FeS8Nl9^7Lxbi*LI z3N_@1s}|F_b|2qdK{#YJ*A_*0T1SyDXwedKD}U77%`5*7{S)ZdrFG;n1bH@Io_vui zwSR9e=7gGGMXJSB$pAbMxyY|^VMGNEG;07DC`gP_P3M%X8I`x4B0LP7oCc0SPJ^gj zpzMivk^7UO1l}QCu3*}~gnN$=WzZEHjfTk~to{&M|Hw3rKtXKiVLA5Ctblpre#OqM zYm3?+I7n5GZ6i6%*9!Z{_T|HV;d>0#4xNIp6}D6qzsUwu0vX|Y^i zK0d7k*88F4Zw&^p7J6TB3I_GoH$?aDe|cq6F@<<}SW3^ZX+s!uJlk=S#2eWoGCL&X z!g8>QnzOIlBnk1CFya$3I`NgF3dO>~kKZwnDog1E7pb458u{l(l{rL$Eejd15~5M} zA}^1Su8)~Wj2_ zS59xEsx%+QG?CUG=l1_z$Z)Ehp(Fga4*4bq7Ev7r1GtIN61VD_EyhSz4c12q6Bbh3 z-x;dFVYSAv7#+4Z3+QC=$*AeP=W!r4@aDcV9)98e3Nz^QN4b`m ztOu2m6b!@nI<&?-nKB9FNNC`G=lRt$vbZHTKi45prLu2JQPs)|0BGG#uOg(k{sooP zk5q2LO2BOd9AX`E>g)KNKbOd6ry_tOe(VMy%^-yIOO{M$ z!u=8jOx)+eeNM4-{4Gxi9Pb<&uVi#MwEi^?;P=w5V(HRxm(|=Ka4q4SY5xSN<2NTf za6i~{r=Bq$eKjWwxC~@D85$Y_5Rpoe;a^ zH3pGpTsW-7e4J0TCv(1kJ5sEbWj`!(UF3zsnzElTtpf=0QKxq6RF#AL1pq6x=J#oe z6xO0s?YuTcJ`ynU@331+BGugIkC2~qIR@M*8FXv6a3z49vm~03Em?QCXS{c4!1N-= zg@>|4_I{ZB_Xluf{j6e!Un-+Mrx1iN_ov<}1z09b$8>CDvw61=Tv}Dtl)0{e6md!=d@uv0XZuu&QZFd3tMV=6hgP#I*J2WtC(g~ zP-tBZQ7-%KK0$2&fL&?x#`j+`p!^+i5IRci(NHxPX0_lOTupjuc-8+3e~s4DPK#7l zRS2Oz7#P7+_V?i*85zBQa``i=&9_`_3`sJ1BboXit&FoLB{9SCwgn}KgGb@8YhhpS zYH}d2Ecp@A;RAP2;#xA{b|a`$?t@gWmuw+pYdknb5|W0?ti^F!n?Pw?pYLa%9w34V z$EJ1n>m9fun!CSWw(#AC$$lK_!iqK6pRh5mqku?KDjR%!1#UlV9-DZ-X%4K5x!9h5 zt*l{SD0L5e0MxqAPHV=Rqw_)HXB`UwKq(FNL)mI>9xs_^$xJrEagxR0%7ro|Z~K}} z&XcK(UlAIXw)pOeJmTx&%2Xw9x!STCgmi`mqkniUG3@%_A7y$OMBQqIk_8JOmHN_g{)_fnoC>?ocqq{f zxT!Scd|m(nC_w%UP{IyEm?rgg&-Z-~+@k4R&=(~vXkGqn_D(W)p4T_4xV1nvv?oL_Uz*01&t4=e!6aS z3}i<7&JHp^13IVbFPsQt*}IT5_$Y4WW4c;o{~G^8odZussd(g%%zq6;Ng>{i)LrXJ zF6tNL%&p}Vb@>YJ^ztA~d?L@|g*#FxzLcsVi1>>byHDB`R0fJ?*cWH_QBI$lr{Lf{ znqt)I$@eSE(}bh5FN@Z(H6HtQ5r8l|ryONG*W3O^8O%i?qSwmMt{QfI7=(ysJC_bvEC?4(%*bqEUDh zPb_2DcGuvWU+EHAv6P7wWQXuX294$r2<=f zXb~2$pHST=bd3P|u}IOZa?(UQk$n>kva=*1YEy=rCAztNaAK=bRamuH+VbF<)oKR; zfE4orHf^wf2tfhBKX`ja4y|nR?zC(|oyx}`c46Ltz%Jl9-R5-q9}E%zkAwuo*?s_# zJ5UPz|LsBKXvFV@Vc1&K#m9TuSlWhDf$jMadQ{46Fc!o_0C&_zLLP)gGyiN93v6ix zT~)5rr^-_^CZ&I0XonCDaw{+z-P{ohs+#+yj*uVS6N(KPYiIovREiD_I#?NovYI66 zxjc?ole~*f&fYNdF7r5Vm%IUfM7(Px`cV1xmnw6x`Ten6)t z-ChfJ(Ng&{QQ@O&y+CW*-C|A5%zFyz`w4}3;e8@dik$uY*Jg98!+c^ikmgV8wfrar z*4LOiQFZl+ygY>tx}3F8hh7-&PSV*=Zc)%_W5-dX#49jD3rwQvJ+^CPtK&Wv-%2fz zu%9Ju#J|@n+hIR6`KsSW>e2CgV=2H{22eNoO z6vU*c$fq`nTwSb^mK$&&Ua!-9p4@6hp-c&g${|5|xl8Gh>h9D>rPs?Q#EO{VX{q_1 z?*f=O|Hk6wUEnzfT}j^_Rm6l3=6hcd`nH~@pdn~Ahxe71p!%)6xYz)sX>q;o2E_1I zk8Vq}_6*JyWZ-XIXu3uP5$G-*pmZur#5y^_{I-o&fi;+z%m2hGpseOBcHNqB}viUH>jiQT$ZK zF(^;ofCrm!51-|66I3`*OU?Ri%W|8O;|vFpjD|Y$1DvGtr&gxA#r0K((~**qeV!*FulGhWXtpd=aELrzWk82lqIwqyE! zo@BSx#4Y~R#;>1HHJhkQkZ(3c_>rS~Xl(u#32Z6hqv%+L|XoBSSC zm5#IDFEw?m}EQ^@fVJ!`yhInOyG)v>0e=DI*J;LUB&pveId`Q6NbqasrmiD_P~VU%$xj)N=vP5l!|KIRW1tDaR) zNr{w_N#SV+-(HBh<%nMJ((N>>{_u-0F4V)wc%Y(!v z5TUCs&Q5O_h{%$5?eXdECZn?&ZQYyIKmwK*Cb;NaN%M?9>>YK%9122VVvZWcRDx&q>fL}a=S{rRum!rmLqJ(6(s z+D1|W2Q8wcmdSg25>!HB-rjwR^R0CPCdY@R7DYzH*GYf;eutbau*R z<{%dPi;cGX+YtgEpdF&wW{PzQG$~We$FV7z9J+ND@CPsU!PUN#($()diYm)?0||iI zoSw4zwHEWs5+75|(X6y>HaJlj(CV}bfjuGUo3RzHX=BQvA;V-CnY&cG_B*X6t(q2_ zbPl>PCtQ#9tE7^hjYVev9S-a(599&GO&@x~O z0s_lM@wG&30JkaW8MfGO?q=Ow%Vq8TV+yXnz>oU;AIx=N6}d@MrDX6%63eFmeHs%7 zFDkh2kXkeN8>}Bh*0~z>+lVi461#nwCE1Z*)+pHLXZ7OEoFuMAwCKCV?RJWI z0Mg586E$a;UuDbyh8D^l(n76sOw)i}pyt_oRSw1dsWj@_gSaTjZBevl4f+S5Dak22 zumHG}_LU9twV(pk6r(t|`4Q!D9{w|sT-5Z}GT<2ZrkQv z`$xNImX1$vQ0Ly(%}4YdNEwS!7t+=>!XJgVot{Q|YSzH7JdbCV>=AK*Ek&wLt*@h>=K9?)qLS2SW~G$=nSz z##!?!OoH2Y8>8H&IAlAMDr%b+AY`)u(D6G#DREZz7P!{F*+E?@6Ln_M>bJ~UKon7g zvv{k^XeCA+{7**qrd#ld?|XF9X-%vM@0IgPK!rPvJ9fdP^Ai8>X3+#&stO1DujwUr z@Mq0ARJ%CswJNPTN{Lsw6aHv}JvJCsABS4=E7Zh#_BMVm>Sdr1iWNn`-ytoxfH>;Z z>9%oE*MTh4tZecla033)I*$HKPyZe!;*T)J$3AndIW}`c{7X9;4LyyRR3{iDnPh(@ zsJEbT{{_k)3_{~WKtDaWN{LBd$WlY*IEo|5{m66;vndW-Izj24E|8FDv^RO6)v&r0 zS06r@f>PSAaM52E?q`is;pTm%m|Pgm2CIoPn16-%f|6IXavVy@u7$xh05=5yc`GFe2`}dYD9SvyS zH;N99V+UBLSX-iZ$n9~)A6osZ!Tt#Sj3|G?{CS078bn1$wS;9tge(nJZB(Lz%6`;) zwdM^Ys-y43=OKT)zb zX9~z%TNDJJH#;-;EnvDEDVN~ALe>Nt248+Uc3rOJbk-{wLJ=ACF(7-08(6?==9tRh zH*9)S`1;Z)+nQQ_hg`u#P*hdZ(4XCoX!ALW%{qtMh*AGji|li2DBqyuw<6BnyH`%pg?}IY zjaM_KKhrBFN}>DSn!EtL^gNRqCqxQ z{hFbnED#vM!hnwGp3fh#fjYPYZG_Fh2J>yVZi&R>oKPv+i-q^HjufAfGx4@3#U+3l zaRubsnKzB4Sp6wmH|;;e%VzX4$nIB;~1VDnGYJeW1Qh!L;EEb`l<* zP9L}x;wAZm6hRyc!Gk_K)2XggG$LrVi>A;zuB|#xtoSp-yoZEs zhbQZvBmptQIM+5x@A)c{oGQB4lpCLitUHd>E9Uo-h|Kzc)y_H_u$_Qubx0U@$e?wj zRmbX*7vVXIX)}4p{6nHZ;9a>wH+3ynatSuj4&a-uiW5rGWlyFRq>;77B9TR8>_^aHuF~A(^czCmIk<7Ek$zAGDyC z19)o8<(nw|m-W3fd;i8-kIvX8ZLscsK3AsC+V?R{aK(IHotb088&~{Hrp}74c4W5a zCi3#bLHx%~$aWb3UzZ!2n;Q-uB5 zt)sT5IvX!k91Hjs2VAg)X2`2b3F^*t?$c-9*hP&NG;p`1vlQR{DMUFty1D5v!w;OO zTYdjKpl~rFnGtJR|OUiagGE&b*i9aGmQ^}kkk}dLA%e1T~>7iAJW%b_)Ldah{W6ECv=ef}2KZol);g>||L^G(I9@_W8^t^SwHa}E0VY5Dk@Qt(clp3k> z>*BtYhERdps5J-$PA5EPvVq6X(|(KiOmM+82bE3m@KaPkNOafVJ>2KZiwm4NyX+*w z;Um^LsL(Y+=;RwQoUk&+2l`9v{tZ8TxL0zWv3!}B1Me?E%*#aOpph0#)OR)<<9>?- zm(NSUKhN{5F6%>-+R;x1xuJz%mxE@WhqHgLQ=nJIFu@DJUA%qR%Wz{X4ubUzd*pkx zr4iwBb%hPHBB2gwUvn1fh9K-VX!O&~C+H6mL;T+TC;!nw=!19A*UMyj+woH(!sh7- zA#grl?E8g_LrelnF#%st{Ru%}l(x@G?MCZQ8*B(haeTi`b0nH_fCmb`% zUXde34y4!0pQej3h`N2lenMIs?{aPgC{Y-R#=UzhB37r8uqP*qyCo43AUBgihF6I0 zx0u)behb<7W%l<3F-fPS=6Lu!W=a@v9vO%dwPOZ?`v^h>As%L9asCfyj{8E51iC_c zG)+a>+}gJC?-_k}Xd>E!tX1y`okA7A48?k*AFco z&Yh!^SvSq(ZCUkL6uf$8hD$czf-w_DuX-lei7bEULRVrkkI&2P&)2k#>eq}syG)3d z92{-VHsk?af>$2clj-Q<1wMe9d}bZLi%TO=6Vb_JorPvskfQAJS5xv{(dSlVI_UMB z8r<}B%!hWB7K6fH^26tau92w9Y(hDoEHwwwapBng+n+yDd*1>^%bYCz`)Q$#E;gsvlc)Waa&swwGu}AMywY}>N7l+g zy=&0H$4i6;{&RK}uIKAglOU@@E^|Ffx6SF$^Z#o5s;D@EuG@hGcY;H3w?PKC;O-7V zg9S~HAi>?;-3AE+m*5r%3=mus+&v5s9B$`Z-+j1i-IxF2pNDyvHK$K?)u~h6ReSFf zUln(ROjOulm`_L^tsPUcTp?MMdUz~4>T8hBx7IEB=j-}@DthmtDDwV5wR@gI$o`c4 z?o^Rx(I8nsklax3oy{IXFYzv`Am-8su%Eht6JjDfF7Z6SvZK1witjF|>Q#RLB4gY=*Mus37VzQ4akqU&0PLuX@K-eNm&;Vu-BZ9dN40=#O=tTHEsLw=w*s&^L61_G6Fv zs#PRa)v(`u`kFBe@7=yOfbk5Ba#vY>@`ymJUvGDjtQGy-j6IDp-XsxqTF2U$-vOK%5f^#R^ANs`QX0;zCJ2lfys=fI44Me*qlf7SaWH4xRH-WzDEp(NY2rgrV zl{yx9vFDpFTbEW?jfjH(*5K|#4)f>bSSJP}(paVja|&NHB1KV~^6ZdUNw+ahR>0g% zC{TBeV@wjwH5u)cM%JOy$|bQ;(zUt+O64z_qzsKyvIVr+zSpO0o%vDE-qM(fH_b7u ziG?42s4)`K!{vmen64^(mr$EhAx_fvtDIh9r!Ce z*=IbkciK>(qPi zk{_y_=6@YKzfq75frZAUv`;vN7@A~4<#FqcjC}7iGwWk6;x6(cmTu6d5FX4(8(iH* zzE1l~$?#K3S(Pd1s?zosbFjCJx;Zp$Q-mEU2nHO=^)Af znB*rdrq|u2fWEG&C}25@nEDRJVg(8efOJAAj;0Lzr!ZdXIyP82c*JQ)I8J*=AI>fv zNpAWxMI{Vem)KO-^n%M1ccMtfn^4|l?DQ!WX+kZfEp{}wY*vSMTOu7c!`9+ea4re= zPxTq4u4|0FH|$^)A{g`7I~k)O1rLRaP(p;>>0p1<+zio|ufs_$wJb5G|KvsZ(%{Wj z3WxK6Y9O}6$p%dAsK_HWsno(U21Xj6b zv%bFhy4kf0n^&)t%ibKV?++-Uhm(AzGtnV%=tc__@w zkMm#)E0c$1R>_Zbz37a{V%9wIcjE~z2}kI?@2e;D<_zz$cBM(*p<;Zn))lc|UA>oc zvn{CPyn1w{1Uut4W(zaZp!wb^N*K`I6Wn9o^4G`mb546RXWFxE2+91Q3wXyrHHq5C z80y_$@3?rR=O&}0=$#^61lgBtUJ2M*-W>ZG$&lCCPN8cA3f$7~LL$Z5xhuYbQy8UR zE7yFfK@OZLN5j>kcPN^Z`YxLRV|D@O>nQK8tLc0PYd?uInqI+ zpix^1F2`>#&WMj%l&xD7XF~S%N##kb1KB0i{rWGk z&b5bGlx2|M;)_yrps}rqxAsPU?^hmn8?IdM%V2I>!Jsi=y&)3_FzoF?)GMA+z8p2u zb|I`3phAlR;G&_nOyv?EKC&GV2Zd!is>?~g(^l5r+IV!S*U29x%Hp>J zz_fOX#?~4V*Yl>(8kQKts-5GGA{$`is%bpPP%% za3@il{u398X9RcS!la4fyv7}QC|Vq$ed|757&g^{aaXB9cTIvMon zU;0AB5eJb_IXJX4EUBN&o!KCknqIPrkN-LF*!hB}wTB-+JCh|yH9Q$>RD~t2&-)fZ zNQs-=Iv#UWp+^qV+9I4|wCmfTUP^$Ho}@2rK#FBlrS~Fmrr1{8M?MK-SFWz=YQKX? zV)Z>vybl71a6fVP^z0ZJG}_LO`i)QYgh!cy%yk;$=8iA&lG2pJ1_vZH_gYAn@mA7n z$+m=$0C;k{Z=+Df{UE;yyt8Q5*<0E>TK=K1b92c_W+MI_)WIZIee~^-L$b4f&oQ;_ zx=3Bws-oDbo`(k%imw{(s1$BvIqXcbo%p9^v3D|jJsC1eMg6Rhn;n^Ir<3;r6z+#H zY!Ug|vra28RmXp@={+}+AtGE3H)HQy%OGU_z4BW*YD2@)?7P%W+{80BZCsAj#Y_<@ zmZO?jg0{M(#Of`={6da!%M8`^Lahu#doBiv;wUX)G`7X$EJ56zRCz^B5I@ZsxF-4I zWcOeuw4Q3cGf zLt`+0k1EkN4iDl;u}8nsNa`7G3ebgrE*%RQ3g{v;THoaBj_h-op#D@*ae;fCj03z( zcur>Q_#bTb&H6b$aS}gXGOI)-RZ=tP9AVm?TDzU>hvtHpD?e&;ixw7WhiK zFj#pblk3W3#bO;>R{Al8DC<2EQoo3<<(h~r%dEY37OhhFk}f%4l<-oXDth1hk`bO` zSv_nSU^|0v%$Y%-bzFd%K$+fdXBg|rGVZrDzxu1H7{Oy>>eu1HY;VzKB-6Jh3^?iD z)uobn@xjzV(gK72{R9XCg4vr4s#nf68F~DAs1M==BRVKjd@eNF`Y*Y^2lEa3j+P!s zkHKhK$70+K{oou2Prq~hiF;!Lb{2gFuN07Iw6SLT&l>7ns`2p2307_ADwRlWJWt2u zVWwOHP4SlcXHWuEvu~+Y_sYO|gGydO4#yQPGYl8~#`<+_DY6*3_3Yua{uuy6_Pg?A zIA+?`IWamo8Sf#elnJq_qC_*{89fytDV4c+#-?K+bHWh*t#EVCAjjNV5xEX6)un{E zK9%(9&QB`y9AY3k&%w3ujL2h(^wMN45BsH1y_#`u9q_mlG2dBdb!cr!3P+qLy049d&p&kZ(Y7$eA6>=oWcL>|Xc{i4%N$6=DYuuGdpTz-`lUkNbiJ&<`0Oe*tX@Tz zGJuKs_MJ`!&)841LqR*4`#hvGS)kYsoE{wU1UiMa z#Obu~YzgMyDRos#LSAo|lv&XX-K4Ae+uZx29lua3ftN~GUX#>^;d=|mv!^TbGZOdF zkFbMjZ%EJ8jt}V{g)H-HKF}jGz*_nFTXiWv6K5-*O+?l-R5lh;5D8!WBsbio z<+!tcg9eT1LVePrR|zH!dqplR*n=0P)s-tAv7I=a^vCJ-H-45?I*UHs7o^vdgDKqJ zgUJN*g};oA0ScQl$uh5gaM{~esV!u3Cg0eFg>J+c_G#ou@n$pys8=boNLh!GLwUt;Qj8f@ej)L%Mlb%g?O(zJeE{xsC+lKP6PCeAqoG}E58U`cU3YE7mvu^WW>aJ#$X;m&xL-#E4 z8Ko;D-i&1I$QWh@j!Ndgn|kYLnA$*#myxHenTTeji`=J(rhQNqtO2)juD^`(Wes#q z+*?exizMGMvUo4qGN@v~L<7=U5dI+w2ZSU0LXCc5a)T-kn^Ghv-}_6c^kT&`s-20q zU4;wAg);cQCLIu(LckBl?qh!EFHT7NPEIeqh4?9nR0374Pb#VO#x1Ez;;jIzm86O> zB#}ZjM*7S~sIFhr44T7sMFQ#-;}!U7_SkGWK;Z8_?Hyu;%G)Aqb;cpHuxB77sGm}o zsE{+NP(&*>lBmc7jUi3iQ}QYj`mBr>o3}K2U>$so(3$A`-t+ zYDH`f%fo6#+~}hTaGp9)(JFHau+0cRgB^y+o~`6R^o0QD#oa$v-%qKl?(pqL z^Hd0;gK$%can)F)lfO*chGx4%sa}BgSz*Pjs0>t@rH(FclA*QDnD&H~mQA6jd>{p`q_2LcQ2!~Ic)|o>IyxFcpY8t49{=_oYDABA{wXqW zsB(-5(iuaEOEEQ>d+8lG^eO0b{8Bzd53#pFi7LKUPr|?;eJRh1#i>Z3Q7d7*2Nx+Y z%bWm2``wZ&(Z4;)0IG&CySVLFnE^Vi?k_JysK#bzHaE$GopL6rAQwP!uQv(P(%F_7@u_#YH<8O868d#zx~c8R0%)&ks{GcAtZsAfER>X{dO;T9F?snlC1GhMq~kHsl=NH z3Y!ee+=M>>4oVVk(vatn?9PSfhrz9cC^6&Y$Ui4%8cl(>^muT6dVRW}753i{4?uG0c zhs}F`*s|Cd485!kFQEKDxejwPm>+y85$rw3Otp2k{kz|0_Hc}GOA46{ML&NDzFyOD z3cEm8mw3D!1^2v#TZmqa@dEhNS|#yVs)uzH0E}We5^?@d0seNc0DQ|8jgw|;9dRzt zHetO1D$-ykK&KLx)Mh47%F^O|)nr=W@=&G+D3o0ARnAzv%@|-7Ul=bb8DKu)f)&@{ zf!M)da!l_u0VRMT+blehZyb#&4b^%Nf%{xfTM)RD3`W3sUx(AaV9*09H=%_n%1G)A zFlV4;O97rJV1wT!|9Qj)HiNEDeqGOMs~4Ski9YR@Bz*Glxcd9|bfts%j{l2-@&uJhW%VECS z`|px15{h0j@T;+|Zdm?MZ||GR@FWHLe6P)@{QUf{#;ue>uBCzl!w9SssH_vy1H=0H z7qJJU$vuz1dYtU-?+*2QqD`LwM1)iGP3zvv$N~_JT)d;a|Tj5W`bzBk_1t>0z1)tN6&pnP*6}Z=GOK| z3V_I%KDzSr^8;}S2^mr>B_gHSVy4+1SOC}!+3Z5klVw8ehy+TpXwy~KQ4Ua@?L@xW z!zS^B*fj;nt7{d|@C9j9>WGPni5E;9xVNu3Usi(N)GibV(RCOUTyo+1_XGey{QK)O z@J}xBi>!&lxDXT^HC0t$P~s)(fHqt$n7KGP`R!)eBBP*C`Tpq}%>ZXuk(ij607d~1 z=RYDMBEV+d#}~KJbmR?%@Yo_ayI8E~PMZRNoRqlf!!a^4>b@BSpZ;!3xgfgx^=oz( z=pUVZ3(?2R8MC|TjF($k=#W?$+Vi|K-@PzmQ?HA0fySs3)}LIT&Q`lVUjy3E>p>{O zA|ey{V#di9z6Z)?!d`#O)z#G<9FCTo9XN=bM^Y_e;p`Qn+TOc3bn4!mssulh=;TEy z6V%G@cLjQuCMUTdkXhhcR1_1$K8lTlBLIO|bw+A+KHVK%uXs`|bM|v>p5Fjs?PXbEgS=M%4#j_IWyNIQfk_?O>TTubXG1lmLg)9W zJ^6*uO8VM?>Rg$xn#+U= z&uoEy&S9h{f58l~Xi=wHfzjBWy9<5N#6iB@Y58aByTRw%ujTnZtOW>>n!RX7)~*#N zCqe}^F8kwqr@J;Af9(vlbSF+um@#bq%dh#D;lg4)9kum)S9?X3Jhiu7L#Y`mMwuVP z%N>HfRgMn^~Qtj5fTA3L<=FIY#0ICph1!K9emZg+6U zQ%^HBuF*!CF2?H!s{Ef-D!NwO$rA^>9~orm`b&@vcx6b#-a~gwq5fBuWz}Q$B)??8j+* z=WAf%Vz{46sZOeuUl_05&A&?!(e^bDljUY`&x32_57!gP9S#9>@zMIQjs2+- zCnqOh?Fs@kjfUcCyYK9SoZU^}3C}kPz`Xv|>H@xcQ$xzg!~`sZHnz4a&Yr*u!9~J; z65Ao2U8^EZ_q|^?v(B z5nsN1AtNWhzq`Yx5Y#xtGf~#?9*v_Na}A5${FiX1hSJ#dP5rLsfOJf>FSo!z^?dWg z!NFn8|NcY^XFI!Z>l`8=P;&D0BPeeVCBNIab<04~ry z>Cllt{Tn(QsEtdxDvZh^Bru@9P5w^mY+=A<2~|<`4}dw9rBWtO=`B+=hmDO5knNbg zGLG=!Q-)vXwQ%;9ONxp%x~^BHe`#xL3#|FNDMm=&mIPX4N=>*b*-xsAo`0Vzkp~t< z`(`6Ep93iHH1lvVF|#B@J}j?)SpPQV9~N`~M*=fqAn48ddSA$YIy91bK|(fj2i Date: Thu, 14 Feb 2019 12:10:34 +0100 Subject: [PATCH 10/16] Correction of format --- vignettes/Anomalies.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index 3c9d9f84..ce4440f2 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -160,7 +160,7 @@ rect(lonNino[1], latNino[1], lonNino[2], latNino[2]) ColorBar(brks, vert = FALSE, subsampleg = 5) -```r +``` The strong signal of El Niño / La Niña (positive / negative anomalies) for november 1997 and 1999 can be clearly seen in both time-series and snap-shot maps. However, since the anomaly is computed with the annual climatological mean, the clear seasonality presented in the time series is not removed hence obscuring the real signal of El Niño variability. We can clarify this problem computing monthly anomalies as described in the next section. @@ -264,7 +264,7 @@ for (isd in 1:dim(exp)[3]) { grid() text(sdates_exp_df[5],2,'With CrossValidation',font = 2, adj = 0) -```r +``` The figure above shows the deseasonalized tas anomalies for El Nino3.4 domain comparing observational ERAinterim data (black colour lines) with 7-months forecasts from the EUROSIP system for several starting dates (different members [grey] and ensemble mean [green]). The lower / upper panels corresponds with the application / not-application of CrossValidation in the computation of the climatology used in the anomaly calculation. It is clearly seen that there is almost no difference in these two cases due to the significantly high number of sdates (9) in the simulation. -- GitLab From d4148d28c1baccf027e5cbbbdf76afdc0b8b41ef Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 12:19:28 +0100 Subject: [PATCH 11/16] Correction of format 2 --- vignettes/Anomalies.md | 62 +++++++++++++++++++++++------------------- 1 file changed, 34 insertions(+), 28 deletions(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index ce4440f2..0a788206 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -110,8 +110,8 @@ obs_Nino <- data_Nino$obs obs_map <- data_map$obs #Selecting temporal period from 1990 to 2010 -i1990to2010 <- seq(1,12*20) -obs_Nino <- array(obs_Nino[,,,i1990to2010], dim = c(1,1,1,length(i1990to2010))) +i1990to2010 <- seq(1, 12*20) +obs_Nino <- array(obs_Nino[ , , , i1990to2010], dim = c(1, 1, 1, length(i1990to2010))) #Retrieving dates dates_obs <- as.Date(data_Nino$Dates$start[i1990to2010]) @@ -139,17 +139,17 @@ layout(matrix(c(1, 1, heights = c(4, 6, 1), widths = c(1, 1)) -plot(x = dates_obs, y = ano_obs_Nino[,,,],type='l',col='grey', - xlab = '', ylab = '', - xlim= c(dates_obs[1], dates_obs[20*12]), ylim=taslim) -lines(x = dates_obs, y = ano_obs_Nino_sm[,,,],col='darkgreen') -lines(x = dates_obs, y = rep(0,length(dates_obs)), col = 'black') -lines( rep(dates_obs[inov97], 2), taslim, lwd = 4, col = 'darkred') -lines(rep(dates_obs[inov99],2), taslim, lwd = 4, col = 'darkblue') +plot(x = dates_obs, y = ano_obs_Nino[ , , , ], type = 'l', col = 'grey', + xlab = NA, ylab = NA, + xlim= c(dates_obs[1], dates_obs[20*12]), ylim = taslim) +lines(x = dates_obs, y = ano_obs_Nino_sm[ , , , ], col = 'darkgreen') +lines(x = dates_obs, y = rep(0, length(dates_obs)), col = 'black') +lines(rep(dates_obs[inov97], 2), taslim, lwd = 4, col = 'darkred') +lines(rep(dates_obs[inov99], 2), taslim, lwd = 4, col = 'darkblue') grid() title('Tas anomaly between 1990 and 2010 for Niño3.4 domain', line = 1.5) -PlotEquiMap(ano_obs_map[1, 1, 1, inov97,,], data_map$lon, data_map$lat, +PlotEquiMap(ano_obs_map[1, 1, 1, inov97, , ], data_map$lon, data_map$lat, brks = brks, drawleg = F, toptitle = 'Tas anomaly (El Niño nov-97)', sizetit = 0.8) rect(lonNino[1], latNino[1], lonNino[2], latNino[2]) @@ -162,6 +162,9 @@ ColorBar(brks, vert = FALSE, subsampleg = 5) ``` + + + The strong signal of El Niño / La Niña (positive / negative anomalies) for november 1997 and 1999 can be clearly seen in both time-series and snap-shot maps. However, since the anomaly is computed with the annual climatological mean, the clear seasonality presented in the time series is not removed hence obscuring the real signal of El Niño variability. We can clarify this problem computing monthly anomalies as described in the next section. @@ -176,7 +179,7 @@ We will also compare two different methods when calculating anomalies, using or ```r #Starting dates parameter for time series sdates_exp_df <- seq(as.Date("1994/11/01"), as.Date("2002/11/01"), "years") -sdates_exp <- format(sdates_exp_df,'%Y%m%d') +sdates_exp <- format(sdates_exp_df, '%Y%m%d') #Loading data data <- Load(var = 'tas', @@ -217,24 +220,24 @@ ano_exp_CV_m <- Mean1Dim(var = ano_exp_NoCV_sm, posdim = 2 ) #Plotting taslim <- c(-2, 2) -layout(matrix(c(1, 1, 2, 2),4,1, byrow = TRUE)) +layout(matrix(c(1, 1, 2, 2), 4, 1, byrow = TRUE)) plot(x = dates_obs, y = rep(0,length(dates_obs)), type = 'l', col = 'black', xlab = 'Forecast time', ylab = 'tas anomalies', xaxt = 'n', - xlim= c(dates_obs[inov97-40], dates_obs[inov99+40]), ylim=taslim) -axis.Date(side = 1, at = sdates_exp_df,format = '%h%Y') + xlim= c(dates_obs[inov97-40], dates_obs[inov99+40]), ylim = taslim) +axis.Date(side = 1, at = sdates_exp_df, format = '%h%Y') -for (isd in 1:dim(exp)[3]) { - jsd = (isd-1)*7+1 - for (im in 1:dim(exp)[2]) { - lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_NoCV[1, im, isd, ], +for (isd in 1 : dim(exp)[3]) { + jsd = (isd-1)*7 + 1 + for (im in 1 : dim(exp)[2]) { + lines(x = dates_exp[jsd : (jsd + 6)], y = ano_exp_NoCV[1, im, isd, ], col = 'grey') } - lines(x = dates_exp[jsd : (jsd+6)], y = ano_obs_NoCV[1,1,isd,], + lines(x = dates_exp[jsd : (jsd + 6)], y = ano_obs_NoCV[1, 1, isd, ], col = 'black', lwd = 2) - lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_NoCV_m[1,isd,], + lines(x = dates_exp[jsd : (jsd + 6)], y = ano_exp_NoCV_m[1, isd, ], col = 'green', lwd = 2) } grid() @@ -247,24 +250,27 @@ plot(x = dates_obs, y = rep(0,length(dates_obs)), type = 'l', col = 'black', xlab = 'Forecast time', ylab = 'tas anomalies', xaxt = 'n', - xlim= c(dates_obs[inov97-40], dates_obs[inov99+40]), ylim=taslim) + xlim= c(dates_obs[inov97-40], dates_obs[inov99+40]), ylim = taslim) axis.Date(side = 1, at = sdates_exp_df,format = '%h%Y') -for (isd in 1:dim(exp)[3]) { - jsd = (isd-1)*7+1 - for (im in 1:dim(exp)[2]) { - lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_CV[1, im, isd, ], +for (isd in 1 : dim(exp)[3]) { + jsd = (isd-1)*7 + 1 + for (im in 1 : dim(exp)[2]) { + lines(x = dates_exp[jsd : (jsd + 6)], y = ano_exp_CV[1, im, isd, ], col = 'grey') } - lines(x = dates_exp[jsd : (jsd+6)], y = ano_obs_CV[1,1,isd,], + lines(x = dates_exp[jsd : (jsd + 6)], y = ano_obs_CV[1, 1, isd, ], col = 'black', lwd = 2) - lines(x = dates_exp[jsd : (jsd+6)], y = ano_exp_CV_m[1,isd,], + lines(x = dates_exp[jsd : (jsd + 6)], y = ano_exp_CV_m[1, isd, ], col = 'green', lwd = 2) } grid() -text(sdates_exp_df[5],2,'With CrossValidation',font = 2, adj = 0) +text(sdates_exp_df[5], 2, 'With CrossValidation', font = 2, adj = 0) ``` + + + The figure above shows the deseasonalized tas anomalies for El Nino3.4 domain comparing observational ERAinterim data (black colour lines) with 7-months forecasts from the EUROSIP system for several starting dates (different members [grey] and ensemble mean [green]). The lower / upper panels corresponds with the application / not-application of CrossValidation in the computation of the climatology used in the anomaly calculation. It is clearly seen that there is almost no difference in these two cases due to the significantly high number of sdates (9) in the simulation. -- GitLab From 5b5309f2b93e3ae4ab8cc5f06719e090aa86fcb0 Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 12:30:20 +0100 Subject: [PATCH 12/16] Correcting format 3 --- vignettes/Anomalies.md | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index 0a788206..b19a3c2b 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -1,11 +1,8 @@ --- -author: "Suso" -date: "`r Sys.Date()`" -output: rmarkdown::html_vignette -vignette: > - %\VignetteEngine{knitr::knitr} - %\VignetteIndexEntry{anomaly_agreement} - \usepackage[utf8]{inputenc} +author: +date: +output: +vignette: --- Anomalies ======== -- GitLab From 3b58caa82db4cd3d7a6af8ef06d8fee6be0a47aa Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 12:34:15 +0100 Subject: [PATCH 13/16] Correcting format 4 --- vignettes/Anomalies.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index b19a3c2b..3299a2f2 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -7,7 +7,7 @@ vignette: Anomalies ======== -Function AnoMultiOption() allows the user different possibilities of computing anomalies for observational and modelled data. This vignette illustrates a complete example of how to load datasets using the Load() function following different cases for anomalies computations. +Function `AnoMultiOption()` allows the user different possibilities of computing anomalies for observational and modelled data. This vignette illustrates a complete example of how to load datasets using the `Load()` function following different cases for anomalies computations. ### 1- Load dependencies @@ -42,7 +42,7 @@ For more information about both datasets see the next documentation: - [**ERAinterim**](https://www.ecmwf.int/en/forecasts/datasets/archive-datasets/reanalysis-datasets/era-interim). - [**EUROSIP system**](https://www.ecmwf.int/en/forecasts/documentation-and-support/long-range/seasonal-forecast-documentation/eurosip-user-guide/multi-model). -Both datasets can be downloaded from the corresponding server. However in this example we will use the BSC esarchive database. Files path parameters for the Load() function are defined as follows: +Both datasets can be downloaded from the corresponding server. However in this example we will use the BSC esarchive database. Files path parameters for the `Load()` function are defined as follows: ```r #Path for experimental data files path @@ -76,7 +76,7 @@ latMap <- c(-30, 30) As a first step we will explore the observational data for the last decades. -In the case of using Load() for loading only observational data, all available data from the +In the case of using `Load()` for loading only observational data, all available data from the given starting date will be loaded as a time series along the forecast time dimension. Anomalies should therefore be computed using along this dimension using parameter `dim_ano = 4`. Parameter `'output'` should be set to `'areave'` to compute El Niño3.4 domain averages and to `'lonlat'` to load snap-shot maps. -- GitLab From 444873d14acc4b137471f65962b2c2bc959d18d8 Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 12:42:54 +0100 Subject: [PATCH 14/16] Correcting format 5 --- vignettes/Anomalies.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index 3299a2f2..65a11f63 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -271,5 +271,5 @@ text(sdates_exp_df[5], 2, 'With CrossValidation', font = 2, adj = 0) The figure above shows the deseasonalized tas anomalies for El Nino3.4 domain comparing observational ERAinterim data (black colour lines) with 7-months forecasts from the EUROSIP system for several starting dates (different members [grey] and ensemble mean [green]). The lower / upper panels corresponds with the application / not-application of CrossValidation in the computation of the climatology used in the anomaly calculation. It is clearly seen that there is almost no difference in these two cases due to the significantly high number of sdates (9) in the simulation. -Overall the agreement between observations and simulations is quite good, especially because in opposition to section 3 this anomalies does not include any seasonality. Since the anomaly is computed along the starting dates dimension, i.e. from the corresponding monthly means of the two decades, it assures that the seasonality is removed. This means that changes in the deseasonalized tas anomalies represented in the figure and thus directly linked to natural variability, like El Niño / La Niña signal, are reasonably well reproduced by the EUROSIP system. +Overall the agreement between observations and simulations is quite good, especially because in opposition to section 3 this anomalies does not include any seasonality. Since the anomaly is computed along the starting dates dimension (`dim_ano = 3`), i.e. from the corresponding monthly means of the two decades, it assures that the seasonality is removed. This means that changes in the deseasonalized tas anomalies represented in the figure and thus directly linked to natural variability, like El Niño / La Niña signal, are reasonably well reproduced by the EUROSIP system. -- GitLab From 024ed7b3a92f3fba0f6d1bb97b6d167458eae065 Mon Sep 17 00:00:00 2001 From: jpena Date: Thu, 14 Feb 2019 12:47:23 +0100 Subject: [PATCH 15/16] Editing final text --- vignettes/Anomalies.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index 65a11f63..6813f732 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -271,5 +271,5 @@ text(sdates_exp_df[5], 2, 'With CrossValidation', font = 2, adj = 0) The figure above shows the deseasonalized tas anomalies for El Nino3.4 domain comparing observational ERAinterim data (black colour lines) with 7-months forecasts from the EUROSIP system for several starting dates (different members [grey] and ensemble mean [green]). The lower / upper panels corresponds with the application / not-application of CrossValidation in the computation of the climatology used in the anomaly calculation. It is clearly seen that there is almost no difference in these two cases due to the significantly high number of sdates (9) in the simulation. -Overall the agreement between observations and simulations is quite good, especially because in opposition to section 3 this anomalies does not include any seasonality. Since the anomaly is computed along the starting dates dimension (`dim_ano = 3`), i.e. from the corresponding monthly means of the two decades, it assures that the seasonality is removed. This means that changes in the deseasonalized tas anomalies represented in the figure and thus directly linked to natural variability, like El Niño / La Niña signal, are reasonably well reproduced by the EUROSIP system. +Overall the agreement between observations and simulations is quite good, especially because in opposition to section 3 this anomalies does not include any seasonality. Since the anomaly is computed along the starting dates dimension (`dim_ano = 3`), i.e. from the corresponding monthly means of the two decades, it assures that the seasonality is removed. This means that changes in the deseasonalized tas anomalies represented in the figure are linked to natural variability, like El Niño / La Niña signal. Therefore, the agreement between observations and model confirms that this natural variability is reasonably well reproduced by the EUROSIP system. -- GitLab From fdcf2c3160992376613de92be28bb42e6accf4b2 Mon Sep 17 00:00:00 2001 From: aho Date: Thu, 28 Mar 2019 12:11:14 +0100 Subject: [PATCH 16/16] Add leadtimemax, remove additional period selection. --- vignettes/Anomalies.md | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/vignettes/Anomalies.md b/vignettes/Anomalies.md index 6813f732..84b7a569 100644 --- a/vignettes/Anomalies.md +++ b/vignettes/Anomalies.md @@ -83,13 +83,15 @@ should be set to `'areave'` to compute El Niño3.4 domain averages and to `'lonl ```r -#Starting dates parameter for observations +#Starting dates and forecast period parameters for observations sdates_obs <- '19900101' +time_length <- 240 #Loading data data_Nino <- Load(var = 'tas', obs = list(path_obs_eraI), sdates = sdates_obs, + leadtimemax = time_length, output = 'areave', lonmin = lonNino[1], lonmax = lonNino[2], latmin = latNino[1], latmax = latNino[2]) @@ -97,6 +99,7 @@ data_Nino <- Load(var = 'tas', data_map <- Load(var = 'tas', obs = list(path_obs_eraI), sdates = sdates_obs, + leadtimemax = time_length, output = 'lonlat', lonmin = lonMap[1], lonmax = lonMap[2], latmin = latMap[1], latmax = latMap[2]) @@ -106,13 +109,6 @@ data_map <- Load(var = 'tas', obs_Nino <- data_Nino$obs obs_map <- data_map$obs -#Selecting temporal period from 1990 to 2010 -i1990to2010 <- seq(1, 12*20) -obs_Nino <- array(obs_Nino[ , , , i1990to2010], dim = c(1, 1, 1, length(i1990to2010))) - -#Retrieving dates -dates_obs <- as.Date(data_Nino$Dates$start[i1990to2010]) - #Computing anomalies along forecast time dimension ano_obs_Nino <- AnoMultiOptions(var_obs = obs_Nino, dim_ano = 4)$ano_obs @@ -135,6 +131,9 @@ layout(matrix(c(1, 1, 3, 2, byrow = TRUE), heights = c(4, 6, 1), widths = c(1, 1)) + +#Define the date for plotting +dates_obs <- data_Nino$Dates$start plot(x = dates_obs, y = ano_obs_Nino[ , , , ], type = 'l', col = 'grey', xlab = NA, ylab = NA, -- GitLab