From 4379f9fe41d80282bece2d027d46103baae081fd Mon Sep 17 00:00:00 2001 From: SkylarMarvel <57526164+SkylarMarvel@users.noreply.github.com> Date: Thu, 8 Feb 2024 15:16:40 -0700 Subject: [PATCH] Function wrappers, new package data, new vignettes (#32) * Function wrappers, new package data, new vignettes * More merge changes --- .gitignore | 1 + DESCRIPTION | 7 + GeoTox.Rproj | 22 - NAMESPACE | 3 +- R/calc_concentration_response.R | 54 +- R/calc_internal_dose.R | 35 +- R/calc_invitro_concentration.R | 13 +- R/data.R | 28 +- R/extract_hill_params.R | 45 -- R/fit_hill.R | 117 +++- R/sample_Css.R | 40 ++ R/simulate_age.R | 53 +- R/simulate_exposure.R | 84 ++- R/simulate_inhalation_rate.R | 59 +- R/simulate_obesity.R | 47 ++ README.Rmd | 131 ----- README.html | 553 ------------------ README.md | 151 +---- data-raw/DATASET.R | 373 +++++++++--- data/geo_tox_data.rda | Bin 413846 -> 191211 bytes man/calc_concentration_response.Rd | 9 +- man/calc_internal_dose.Rd | 25 +- man/calc_invitro_concentration.Rd | 2 - man/extract_hill_params.Rd | 29 - man/fit_hill.Rd | 31 +- man/geo_tox_data.Rd | 28 +- man/sample_Css.Rd | 21 + man/simulate_age.Rd | 19 +- man/simulate_exposure.Rd | 26 +- man/simulate_inhalation_rate.Rd | 30 +- man/simulate_obesity.Rd | 41 ++ ...-extract_hill_params.R => test-fit_hill.R} | 5 +- vignettes/dev-conc-resp.Rmd | 84 --- vignettes/dev-sensitivity.Rmd | 120 ---- vignettes/introduction.Rmd | 231 ++++++++ vignettes/package_data.Rmd | 347 +++++++++++ 36 files changed, 1493 insertions(+), 1371 deletions(-) delete mode 100644 GeoTox.Rproj delete mode 100644 R/extract_hill_params.R create mode 100644 R/sample_Css.R create mode 100644 R/simulate_obesity.R delete mode 100644 README.html delete mode 100644 man/extract_hill_params.Rd create mode 100644 man/sample_Css.Rd create mode 100644 man/simulate_obesity.Rd rename tests/testthat/{test-extract_hill_params.R => test-fit_hill.R} (70%) delete mode 100644 vignettes/dev-conc-resp.Rmd delete mode 100644 vignettes/dev-sensitivity.Rmd create mode 100644 vignettes/introduction.Rmd create mode 100644 vignettes/package_data.Rmd diff --git a/.gitignore b/.gitignore index 074696c..66f70ef 100644 --- a/.gitignore +++ b/.gitignore @@ -22,6 +22,7 @@ /*.Rcheck/ # RStudio files +*.Rproj .Rproj.user/ # produced vignettes diff --git a/DESCRIPTION b/DESCRIPTION index aab4347..ac1f3a6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -11,9 +11,15 @@ License: MIT + file LICENSE Suggests: dplyr, ggplot2, + ggpubr, + httk, knitr, + readr, + readxl, rmarkdown, scales, + sf, + stringr, testthat (>= 3.0.0), tidyr Config/testthat/edition: 3 @@ -21,6 +27,7 @@ Encoding: UTF-8 Roxygen: list(markdown = TRUE) RoxygenNote: 7.2.3 Imports: + methods, truncnorm VignetteBuilder: knitr Depends: diff --git a/GeoTox.Rproj b/GeoTox.Rproj deleted file mode 100644 index 69fafd4..0000000 --- a/GeoTox.Rproj +++ /dev/null @@ -1,22 +0,0 @@ -Version: 1.0 - -RestoreWorkspace: No -SaveWorkspace: No -AlwaysSaveHistory: Default - -EnableCodeIndexing: Yes -UseSpacesForTab: Yes -NumSpacesForTab: 2 -Encoding: UTF-8 - -RnwWeave: Sweave -LaTeX: pdfLaTeX - -AutoAppendNewline: Yes -StripTrailingWhitespace: Yes -LineEndingConversion: Posix - -BuildType: Package -PackageUseDevtools: Yes -PackageInstallArgs: --no-multiarch --with-keep.source -PackageRoxygenize: rd,collate,namespace diff --git a/NAMESPACE b/NAMESPACE index 75dd971..d5e0b62 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,8 +3,9 @@ export(calc_concentration_response) export(calc_internal_dose) export(calc_invitro_concentration) -export(extract_hill_params) export(fit_hill) +export(sample_Css) export(simulate_age) export(simulate_exposure) export(simulate_inhalation_rate) +export(simulate_obesity) diff --git a/R/calc_concentration_response.R b/R/calc_concentration_response.R index f98eafe..3490887 100644 --- a/R/calc_concentration_response.R +++ b/R/calc_concentration_response.R @@ -1,9 +1,8 @@ #' Calculate the mixture response from one of three different approaches: #' IA, GCA ,or Hazard Quotient #' -#' @param resp data frame with columns "tp", "tp.sd", "logAC50", "logAC50.sd", -#' "resp_max", "logc_min", "logc_max". -#' @param concentration concentration +#' @param C_invitro in vitro concentrations +#' @param hill_params output from `fit_hill()` #' @param tp_b_mult upper bound multiplier for tp rtruncnorm #' @param fixed if TRUE, sd = 0 #' @@ -11,44 +10,61 @@ #' Calculate the combined response of multiple chemicals. It calculates the #' generalized concentration addition response, the independent action #' response, and a hazard quotient +#' #' @return data frame +#' #' @export calc_concentration_response <- function( - resp, concentration, tp_b_mult = 1.5, fixed = FALSE + C_invitro, hill_params, tp_b_mult = 1.5, fixed = FALSE ) { - # TODO make inputs more general - # TODO add input for linear/log concentration - # TODO maybe add option for ln(resp) or not + if (methods::is(C_invitro, "matrix")) { + .calc_concentration_response(C_invitro, hill_params, tp_b_mult, fixed) + } else { + mapply( + .calc_concentration_response, + C_invitro = C_invitro, + hill_params = list(hill_params), + tp_b_mult = tp_b_mult, + fixed = fixed, + SIMPLIFY = FALSE + ) + } + +} + +.calc_concentration_response <- function( + C_invitro, hill_params, tp_b_mult, fixed +) { interval <- c(-50,50) # TODO value of b not consistent # grep "tp.sim <-" ~/github/GeoToxMIE/*.R - tp <- t(sapply(1:nrow(concentration), function(x) { + tp <- t(sapply(1:nrow(C_invitro), function(x) { truncnorm::rtruncnorm( 1, a = 0, - b = resp$resp_max * tp_b_mult, - mean = resp$tp, - sd = if (fixed) 0 else resp$tp.sd + b = hill_params$resp_max * tp_b_mult, + mean = hill_params$tp, + sd = if (fixed) 0 else hill_params$tp.sd ) })) - logAC50 <- t(sapply(1:nrow(concentration), function(x) { + logAC50 <- t(sapply(1:nrow(C_invitro), function(x) { truncnorm::rtruncnorm( 1, - a = resp$logc_min - 2, - b = resp$logc_max + 0.5, - mean = resp$logAC50, - sd = if (fixed) 0 else resp$logAC50.sd + a = hill_params$logc_min - 2, + b = hill_params$logc_max + 0.5, + mean = hill_params$logAC50, + sd = if (fixed) 0 else hill_params$logAC50.sd ) })) - GCA.eff <- IA.eff <- GCA.HQ.10 <- IA.HQ.10 <- rep(NA, nrow(concentration)) - for (i in 1:nrow(concentration)) { + GCA.eff <- IA.eff <- GCA.HQ.10 <- IA.HQ.10 <- rep(NA, nrow(C_invitro)) + for (i in 1:nrow(C_invitro)) { - C_i <- concentration[i, ] + C_i <- C_invitro[i, ] tp_i <- tp[i, ] AC50_i <- 10^logAC50[i, ] diff --git a/R/calc_internal_dose.R b/R/calc_internal_dose.R index e2fc786..3250418 100644 --- a/R/calc_internal_dose.R +++ b/R/calc_internal_dose.R @@ -1,24 +1,53 @@ #' Calculate internal chemical dose #' #' @description -#' Estimate the internal dose from inhalation of a chemical given inhalation rate, time, and body weight +#' Estimate the internal dose from inhalation of a chemical given inhalation +#' rate, time, and body weight #' #' @param C_ext ambient chemical concentration in \eqn{\frac{mg}{m^3}} #' @param IR inhalation rate in \eqn{\frac{m^3}{day}} #' @param time total time in \eqn{days} #' @param BW body weight in \eqn{kg} +#' @param scaling scaling factor encompassing any required unit adjustments #' #' @details #' TODO Additional details... #' \deqn{D_{int} = \frac{C_{ext} \,\times\, IR \,\times\, time}{BW}} #' #' @return internal chemical dose in \eqn{\frac{mg}{kg}} +#' +#' @examples +#' n_chem <- 3 +#' n_sample <- 5 +#' +#' # Single population +#' C_ext <- matrix(runif(n_sample * n_chem), ncol = n_chem) +#' IR <- runif(n_sample) +#' calc_internal_dose(C_ext, IR) +#' +#' # Multiple populations +#' C_ext <- list( +#' "a" = matrix(runif(n_sample * n_chem), ncol = n_chem), +#' "b" = matrix(runif(n_sample * n_chem), ncol = n_chem) +#' ) +#' IR <- list(runif(n_sample), runif(n_sample)) +#' calc_internal_dose(C_ext, IR) +#' #' @export -calc_internal_dose <- function(C_ext, IR, time = 1, BW = 1) { +calc_internal_dose <- function(C_ext, IR, time = 1, BW = 1, scaling = 1) { # TODO How to handle inputs with different units? # e.g. simulated inhalation rate is in m^3/(day * kg), so BW isn't needed # TODO paper states t = 365 in section 2.3, also states that C_ss achieved # in 1 day and repeated exposure accumulates additively. Computation done # with t = 1, is that correct? - C_ext * IR * time / BW + + if (methods::is(C_ext, "matrix")) { + .calc_internal_dose(C_ext, IR, time, BW, scaling) + } else { + mapply(.calc_internal_dose, C_ext, IR, time, BW, scaling, SIMPLIFY = FALSE) + } +} + +.calc_internal_dose <- function(C_ext, IR, time, BW, scaling) { + C_ext * IR * time / BW * scaling } diff --git a/R/calc_invitro_concentration.R b/R/calc_invitro_concentration.R index 61a0d40..cf78728 100644 --- a/R/calc_invitro_concentration.R +++ b/R/calc_invitro_concentration.R @@ -10,16 +10,25 @@ #' TODO Additional details... #' \deqn{C_{plasma} = C_{ss} \,\times\, D_{int}} #' -#' TODO If C_ss is not provided, uses \emph{httk} to create... -#' #' @return \emph{in vitro} equivalent plasma concentration in \eqn{\mu M} #' @export calc_invitro_concentration <- function(D_int, C_ss = NULL) { + if (is.null(C_ss)) { # TODO add real-time computation of Css values stop("real-time computation of Css values has not been implemented") } + # TODO the current C_ss data passed into this for step 01-Sensitivity.R # doesn't match the ages that were simulated? + + if (methods::is(D_int, "matrix")) { + .calc_invitro_concentration(D_int, C_ss) + } else { + mapply(.calc_invitro_concentration, D_int, C_ss, SIMPLIFY = FALSE) + } +} + +.calc_invitro_concentration <- function(D_int, C_ss) { D_int * C_ss } diff --git a/R/data.R b/R/data.R index f5501ff..4cab00a 100644 --- a/R/data.R +++ b/R/data.R @@ -1,21 +1,19 @@ #' Geo Tox data #' -#' A subset of data from the GeoToxMIE source scripts +#' Sample data for use in vignettes and function examples. See the +#' "package_data" vignette for details on how this data was gathered. #' -#' @format A list with 4 items: +#' @format A list with items: #' \describe{ -#' \item{ice}{Description} -#' \item{MC_iter}{Description} -#' \item{cyp1a1}{Description} -#' \item{age}{Description} -#' \item{Css}{Description} -#' } -#' Css is a list with 5 items: -#' \describe{ -#' \item{Css$age}{Description} -#' \item{Css$obesity}{Description} -#' \item{Css$httk}{Description} -#' \item{Css$'dose-resp'}{Description} -#' \item{Css$'ext-conc'}{Description} +#' \item{exposure}{2019 AirToxScreen exposure concentrations for a subset of +#' chemicals in North Carolina counties} +#' \item{dose_response}{A subset of data from InvitroDB for the +#' "LTEA_HepaRG_CYP1A1_up" assay} +#' \item{age}{County population estimates for 7/1/2019 in North Carolina} +#' \item{obesity}{CDC PLACES obesity data for North Carolina counties in +#' 2020} +#' \item{simulated_css}{Simulated steady-state plasma concentrations for +#' various age groups and obesity status combinations} +#' \item{boundaries}{County and state boundaries for North Carolina in 2019} #' } "geo_tox_data" diff --git a/R/extract_hill_params.R b/R/extract_hill_params.R deleted file mode 100644 index 87d641a..0000000 --- a/R/extract_hill_params.R +++ /dev/null @@ -1,45 +0,0 @@ -#' Extract results from 2-parameter fit_hill outputs -#' -#' @param fits list of outputs from fit_hill(..., fixed_slope = TRUE) -#' -#' @return data frame -#' @export -#' -#' @examples -#' df <- data.frame( -#' logc = -3:3, -#' resp = 5 / (1 + 10^(1.2 * (0.4 - rep(-3:3, each = 3)))) + rnorm(21) -#' ) -#' fit <- fit_hill(df$logc, df$resp) -#' extract_hill_params(list(fit)) -#' -#' df_list <- list(df, df) -#' fits <- lapply(df_list, function(df) {fit_hill(df$logc, df$resp)}) -#' extract_hill_params(fits) -extract_hill_params <- function(fits) { - - # TODO add ability to extract from 3-param fits? - - # Convert list into data frame - df_params <- do.call(rbind, lapply(fits, function(fit) { - fit_params <- as.data.frame(t(unlist(fit))) - fit_params <- fit_params[c( - "par.tp", "sds.tp", "par.logAC50", "sds.logAC50", - "logc_min", "logc_max", "resp_min", "resp_max", "AIC" - )] - names(fit_params)[c(1:4)] <- c("tp", "tp.sd", "logAC50", "logAC50.sd") - fit_params - })) - - # Impute sd NAs with mean - df_params$tp.sd.imputed <- is.na(df_params$tp.sd) - idx <- df_params$tp.sd.imputed - df_params$tp.sd[idx] <- df_params$tp[idx] - - df_params$logAC50.sd.imputed <- is.na(df_params$logAC50.sd) - idx <- df_params$logAC50.sd.imputed - df_params$logAC50.sd[idx] <- df_params$logAC50[idx] - - # Return - df_params -} diff --git a/R/fit_hill.R b/R/fit_hill.R index d257c02..1620376 100644 --- a/R/fit_hill.R +++ b/R/fit_hill.R @@ -1,12 +1,67 @@ #' Fit 2- or 3-parameter Hill model #' -#' @param log10_conc base-10 log scale concentration -#' @param resp response -#' @param fixed_slope if TRUE, slope is fixed at 1 +#' @param x data frame or list of data frames. +#' @param conc column name of base-10 log scaled concentration. +#' @param resp column name of response. +#' @param fixed_slope if TRUE, slope is fixed at 1. +#' +#' @return Fit parameters and other stats. "tp" is the top asymptote and +#' "logAC50" is the 50% response concentration. If the computation of the +#' standard deviations of these two parameters fails, then the standard +#' deviation is set equal to the parameter estimate and is indicated by +#' the respective imputed flag being TRUE. +#' +#' @examples +#' # Single chemical +#' data <- geo_tox_data$dose_response[ +#' geo_tox_data$dose_response$casn == "1582-09-8", +#' ] +#' fit_hill(data) +#' +#' # Multiple chemicals +#' data <- split(geo_tox_data$dose_response, ~casn) +#' # 2-param +#' fit_hill(data) +#' # 3-param +#' fit_hill(data, fixed_slope = FALSE) #' -#' @return fit and other stats #' @export -fit_hill <- function(log10_conc, resp, fixed_slope = TRUE) { +fit_hill <- function(x, conc = "logc", resp = "resp", fixed_slope = TRUE) { + + if (!any(methods::is(x, "data.frame"), methods::is(x, "list"))) { + stop("x must be a data.frame or list") + } + + if (methods::is(x, "data.frame")) { + + if (!all(c(conc, resp) %in% names(x))) { + stop("x must contain columns named by 'conc' and 'resp' inputs") + } + + fit <- .fit_hill(x, conc, resp, fixed_slope) + .extract_hill_params(list(fit)) + + } else { + + if (!all(sapply(x, function(df) all(c(conc, resp) %in% names(df))))) { + stop("x data frames must contain columns named by 'conc' and 'resp' inputs") + } + + fit <- lapply(x, function(df) .fit_hill(df, conc, resp, fixed_slope)) + out <- .extract_hill_params(fit) + out <- cbind("name" = rownames(out), out) + rownames(out) <- NULL + # Have consistent output order + out <- out[order(out$name), , drop = FALSE] + out + + } +} + +.fit_hill <- function(x, conc, resp, fixed_slope) { + + log10_conc <- x[[conc]] + resp <- x[[resp]] # Compute initial values resp_medians <- tapply(resp, log10_conc, stats::median) @@ -57,10 +112,18 @@ fit_hill <- function(log10_conc, resp, fixed_slope = TRUE) { ) ) + sds <- suppressWarnings(sqrt(diag(solve(-fit$hessian)))) + if (fixed_slope) { + par <- c(fit$par[1:2], 1, fit$par[3]) + sds <- c(sds[1:2], 0, sds) + } else { + par <- fit$par + } + # Return results out <- list( - par = fit$par, - sds = sqrt(diag(solve(-fit$hessian))), + par = par, + sds = sds, val = fit$value, convergence = fit$convergence, AIC = 2 * length(fit$par) - 2 * fit$value, @@ -69,11 +132,43 @@ fit_hill <- function(log10_conc, resp, fixed_slope = TRUE) { resp_max = resp_max, resp_min = resp_min ) - if (fixed_slope) { - names(out$par) <- names(out$sds) <- c("tp", "logAC50", "t-error") + names(out$par) <- names(out$sds) <- c("tp", "logAC50", "slope", "t-error") + + out +} + +.extract_hill_params <- function(fits) { + + # Convert list into data frame + df_params <- do.call(rbind, lapply(fits, function(fit) { + fit_params <- as.data.frame(t(unlist(fit))) + fit_params <- fit_params[c( + "par.tp", "sds.tp", "par.logAC50", "sds.logAC50", "par.slope", + "sds.slope", "logc_min", "logc_max", "resp_min", "resp_max", "AIC" + )] + names(fit_params)[c(1:6)] <- c( + "tp", "tp.sd", "logAC50", "logAC50.sd", "slope", "slope.sd" + ) + fit_params + })) + + # Impute sd NAs with mean for tp and logAC50 + if (all(is.na(df_params$tp.sd))) { + df_params$tp.sd.imputed <- FALSE } else { - names(out$par) <- names(out$sds) <- c("tp", "logAC50", "slope", "t-error") + df_params$tp.sd.imputed <- is.na(df_params$tp.sd) + idx <- df_params$tp.sd.imputed + df_params$tp.sd[idx] <- df_params$tp[idx] } - out + if (all(is.na(df_params$logAC50.sd))) { + df_params$logAC50.sd.imputed <- FALSE + } else { + df_params$logAC50.sd.imputed <- is.na(df_params$logAC50.sd) + idx <- df_params$logAC50.sd.imputed + df_params$logAC50.sd[idx] <- df_params$logAC50[idx] + } + + # Return + df_params } diff --git a/R/sample_Css.R b/R/sample_Css.R new file mode 100644 index 0000000..c9f0db5 --- /dev/null +++ b/R/sample_Css.R @@ -0,0 +1,40 @@ +#' Title +#' +#' @param simulated_css x +#' @param age x +#' @param obesity x +#' +#' @return x +#' @export +sample_Css <- function(simulated_css, age, obesity) { + + if (is.array(age)) age <- list(age) + if (is.array(obesity)) obesity <- list(obesity) + + mapply( + function(age, obesity) { + out <- matrix(NA, nrow = length(age), ncol = length(simulated_css)) + for (i in 1:length(simulated_css)) { + df <- simulated_css[[i]] + row_idx <- mapply( + function(age, obesity) { + idx <- utils::tail(which(df$age_min <= age & df$weight == obesity), 1) + ifelse(length(idx) == 1, idx, NA) + }, + age, + obesity + ) + for (j in sort(unique(row_idx))) { + idx <- which(row_idx == j) + out[idx, i] <- sample(df$css[[j]], length(idx), replace = TRUE) + } + } + colnames(out) <- names(simulated_css) + # Have consistent output order + out[, order(colnames(out)), drop = FALSE] + }, + age, + obesity, + SIMPLIFY = FALSE + ) +} diff --git a/R/simulate_age.R b/R/simulate_age.R index 133ef83..f9a96d5 100644 --- a/R/simulate_age.R +++ b/R/simulate_age.R @@ -1,26 +1,52 @@ -#' Simulate age data +#' Simulate ages #' -#' @param x data frame containing population data for age groups. +#' @param x data frame or list of data frames containing population data for age +#' groups. Each data frame must contain columns "AGEGRP" and "TOT_POP". #' @param n simulated sample size. #' -#' @return array of simulated ages. -#' @export +#' @return Array or list of arrays containing simulated ages. #' #' @examples -#' x <- data.frame(AGEGRP = 0:18, TOT_POP = c(sum(1:18), 1:18)) -#' simulate_age(x, 10) +#' # Single data frame +#' data <- geo_tox_data$age[geo_tox_data$age$FIPS == 37001, ] +#' ages <- simulate_age(data, 10) +#' +#' # List of data frames +#' data <- split(geo_tox_data$age, ~FIPS) +#' ages <- simulate_age(data, 10) +#' +#' @export simulate_age <- function(x, n = 1e3) { - # Check columns - if (!all(c("AGEGRP", "TOT_POP") %in% names(x))) { - stop("x must contain columns 'AGEGRP' and 'TOT_POP'") + if (!any(methods::is(x, "data.frame"), methods::is(x, "list"))) { + stop("x must be a data.frame or list") } - # Enforce order and check values of age group - x <- x[order(x$AGEGRP), ] - if (!(nrow(x) == 19 && all(x$AGEGRP == 0:18))) { - stop("x$AGEGRP must contain values 0:18") + if (methods::is(x, "data.frame")) { + + if (!all(c("AGEGRP", "TOT_POP") %in% names(x))) { + stop("x must contain columns 'AGEGRP' and 'TOT_POP'") + } + + .simulate_age(x, n) + + } else { + + if ( + !all(sapply(x, function(df) all(c("AGEGRP", "TOT_POP") %in% names(df)))) + ) { + stop("x data frames must contain columns 'AGEGRP' and 'TOT_POP'") + } + + lapply(x, function(df) .simulate_age(df, n)) + } +} + +.simulate_age <- function(x, n) { + + # Ensure order for AGEGRP + x <- x[order(x$AGEGRP), ] # Probability of each age group prob <- x$TOT_POP[-1] / x$TOT_POP[1] @@ -30,5 +56,4 @@ simulate_age <- function(x, n = 1e3) { # Sample ages sample(0:89, size = n, prob = prob, replace = TRUE) - } diff --git a/R/simulate_exposure.R b/R/simulate_exposure.R index ed83913..656b468 100644 --- a/R/simulate_exposure.R +++ b/R/simulate_exposure.R @@ -1,21 +1,85 @@ #' Simulate external exposure #' -#' @param mean array of concentration mean values. -#' @param sd array of concentration standard deviations, leave as 0 for fixed -#' external concentrations. +#' @param x data frame or list of data frames. +#' @param mean column name of mean values. +#' @param sd column name of standard deviations. +#' @param label column name of labeling term, required if `x` has more than one +#' row. #' @param n simulated sample size. #' -#' @return matrix \eqn{_{n \,\times\, length(mean)}} of inhalation rates. -#' @export +#' @return A matrix or list of matrices containing inhalation rates. Matrix +#' columns are named using the values in the `label` column for more than one +#' data frame row. #' #' @examples -#' simulate_exposure(mean = 1:3, n = 10) -#' simulate_exposure(mean = c(1, 10, 100), sd = c(0, 1, 5), n = 10) -simulate_exposure <- function(mean, sd = 0, n = 1e3) { +#' data <- split(geo_tox_data$exposure, ~FIPS) +#' +#' # Single data frame +#' simulate_exposure(data[[1]], n = 5) +#' +#' # List of data frames +#' simulate_exposure(data[1:3], n = 5) +#' +#' @export +simulate_exposure <- function( + x, mean = "mean", sd = "sd", label = "casn", n = 1e3 +) { + + if (!any(methods::is(x, "data.frame"), methods::is(x, "list"))) { + stop("x must be a data.frame or list") + } + + if (methods::is(x, "data.frame")) { + + if (!all(c(mean, sd) %in% names(x))) { + stop("x must contain columns named by 'mean' and 'sd' inputs") + } + + if (nrow(x) > 1 & !(label %in% names(x))) { + stop("x must contain a column named by the 'label' input") + } + + out <- .simulate_exposure(x, mean, sd, n) + if (nrow(x) > 1) { + colnames(out) <- x[[label]] + # Have consistent output order + out <- out[order(colnames(out))] + } + out + + } else { + + if (!all(sapply(x, function(df) all(c(mean, sd) %in% names(df))))) { + stop("x data frames must contain columns named by 'mean' and 'sd' inputs") + } + + nrow <- sapply(x, nrow) + label_exists <- sapply(x, function(df) label %in% names(df)) + + if (!all(nrow > 1 & label_exists)) { + stop("x data frames must contain a column named by the 'label' input") + } + + lapply(x, function(df) { + out <- .simulate_exposure(df, mean, sd, n) + if (nrow(df) > 1) { + colnames(out) <- df[[label]] + # Have consistent output order + out <- out[, order(colnames(out)), drop = FALSE] + } + out + }) + + } +} + +.simulate_exposure <- function(x, mean, sd, n) { + + mean <- x[[mean]] + sd <- x[[sd]] + if (length(mean) == 0) { matrix(0, nrow = n, ncol = 0) - } else if (length(sd) == 1 && sd == 0) { - t(replicate(n, mean)) } else { mapply( function(mean, sd, n) { diff --git a/R/simulate_inhalation_rate.R b/R/simulate_inhalation_rate.R index 6cd10f0..c2506c4 100644 --- a/R/simulate_inhalation_rate.R +++ b/R/simulate_inhalation_rate.R @@ -1,21 +1,38 @@ #' Simulate inhalation rates #' -#' @param ages array of ages. -#' @param params data frame with columns "age", "mean" and "sd". The age column -#' should be in ascending order and represent the lower value of age groups for -#' the corresponding mean and sd values. +#' @param x array or list of arrays containing ages. +#' @param params (optional) data frame with columns "age", "mean" and "sd". See +#' details for more information. #' -#' @return array of inhalation rates. -#' @export +#' @details +#' The age column of the optional `params` data frame should be in ascending +#' order and represent the lower value of age groups for the corresponding mean +#' and sd values. When not provided, the default values will come from Table 6.7 +#' of EPA's 2011 Exposure Factors Handbook using the mean of male and female +#' values. +#' +#' @return Array or list of arrays containing inhalation rates. #' #' @examples -#' simulate_inhalation_rate(c(1, 6, 20)) -simulate_inhalation_rate <- function(ages, params = NULL) { +#' # Single array +#' ages <- sample(1:100, 6, replace = TRUE) +#' simulate_inhalation_rate(ages) +#' +#' # List of arrays +#' ages <- list( +#' sample(1:100, 5, replace = TRUE), +#' sample(1:100, 3, replace = TRUE) +#' ) +#' simulate_inhalation_rate(ages) +#' +#' @export +simulate_inhalation_rate <- function(x, params = NULL) { if (is.null(params)) { # Data comes from https://www.epa.gov/sites/default/files/2015-09/documents/efh-chapter06.pdf - # Table 6.7 Distribution percentiles of physiological daily inhalation rates per unit - # body weight (m3/kg-day) for free living normal weight males and females aged 2 months to 96 years + # Table 6.7 Distribution percentiles of physiological daily inhalation rates + # per unit body weight (m3/kg-day) for free living normal weight males and + # females aged 2 months to 96 years params <- as.data.frame(rbind( c( 0, 0.495, 0.08, 0.48, 0.075), c( 1, 0.48, 0.06, 0.45, 0.08), @@ -38,18 +55,28 @@ simulate_inhalation_rate <- function(ages, params = NULL) { } } - inhalation_rate <- rep(NA, length(ages)) - # TODO keep hard-coded max age? - age_idx <- is.numeric(ages) & (ages >= 0 & ages < 100) + if (methods::is(x, "list")) { + lapply(x, function(y) .simulate_inhalation_rate(y, params)) + } else { + .simulate_inhalation_rate(x, params) + } +} + +.simulate_inhalation_rate <- function(x, params) { + + out <- rep(NA, length(x)) + + age_idx <- is.numeric(x) & (x >= 0 & x < 100) + if (any(age_idx)) { param_idx <- sapply( - ages[age_idx], + x[age_idx], function(age) max(which(age >= params$age)) ) - inhalation_rate[age_idx] <- truncnorm::rtruncnorm( + out[age_idx] <- truncnorm::rtruncnorm( 1, 0, Inf, params$mean[param_idx], params$sd[param_idx] ) } - inhalation_rate + out } diff --git a/R/simulate_obesity.R b/R/simulate_obesity.R new file mode 100644 index 0000000..3bc81fb --- /dev/null +++ b/R/simulate_obesity.R @@ -0,0 +1,47 @@ +#' Simulate obesity +#' +#' @param x data frame containing obesity data as a percentage from 0 to 100. +#' @param prev column name of prevalence. +#' @param sd column name of standard deviation. +#' @param label column name of labeling term, required if `x` has more than one +#' row. +#' @param n simulated sample size. +#' +#' @return Array or named list of arrays containing simulated obesity status. +#' +#' @examples +#' simulate_obesity(geo_tox_data$obesity[1, ], n = 4) +#' +#' simulate_obesity(geo_tox_data$obesity[1:5, ], n = 10) +#' +#' df <- data.frame(label = letters[1:3], prev = c(20, 50, 80), sd = c(5, 5, 5)) +#' simulate_obesity(df, label = "label", prev = "prev", sd = "sd", n = 10) +#' +#' @export +simulate_obesity <- function( + x, prev = "OBESITY_CrudePrev", sd = "OBESITY_SD", label = "FIPS", n = 1e3 +) { + + if (!all(c(prev, sd) %in% names(x))) { + stop("x must contain columns named by 'prev' and 'sd' inputs") + } + + if (nrow(x) > 1 & !(label %in% names(x))) { + stop("x must contain a column named by the 'label' input") + } + + out <- mapply( + function(mean, sd) { + p <- stats::rnorm(n, mean, sd) / 100 + p[p < 0] <- 0 + p[p > 1] <- 1 + status <- stats::rbinom(n, 1, p) + ifelse(status, "Obese", "Normal") + }, + mean = x[[prev]], + sd = x[[sd]], + SIMPLIFY = FALSE + ) + + if (nrow(x) > 1) stats::setNames(out, x[[label]]) else out[[1]] +} diff --git a/README.Rmd b/README.Rmd index b940cdd..45703c9 100644 --- a/README.Rmd +++ b/README.Rmd @@ -34,134 +34,3 @@ You can install the development version of GeoTox from [GitHub](https://github.c # install.packages("devtools") devtools::install_github("Spatiotemporal-Exposures-and-Toxicology/GeoTox") ``` - -## Example - -```{r example} -library(GeoTox) -library(dplyr, warn.conflicts = FALSE) -``` - -### Estimate chemical concentration-response curves - -```{r chem} -conc <- 10^rep(-2:2, each = 3) -tp <- 100 # top asymptote -ga <- 1.6 # AC50 -gw <- 1.2 # slope -resp <- tp / (1 + (ga / conc)^gw) + rnorm(length(conc), sd = 5) - -fit_2param <- fit_hill(log10(conc), resp) # slope fixed at 1 -fit_3param <- fit_hill(log10(conc), resp, fixed_slope = FALSE) - -rbind( - "inputs" = c(tp, log10(ga), gw, NA), - "3-param" = c(fit_3param$par), - "2-param" = c(fit_2param$par[1:2], 1, fit_2param$par[3]) -) -``` - -### Estimate population dose-response - -Input data - -```{r input_data} -# Number of samples to simulate -MC_iter <- 10 - -# Number of chemicals to simulate -n_chem <- 4 - -# Create age groups and group sizes -age <- data.frame( - AGEGRP = 0:18, - TOT_POP = c(0, round(runif(18, max = 1000))) -) -age$TOT_POP[1] <- sum(age$TOT_POP[-1]) - -# Create chemical exposure mean and sd -exposure <- data.frame( - mean = (1 + runif(n_chem))*1e-6, - sd = (1 + runif(n_chem))*1e-7 -) - -# Create chemical concentration-response data -conc_resp <- lapply(1:n_chem, function(idx) { - conc <- 10^rep(-2:2, each = 3) - tp <- 100 + rnorm(1, sd = 15) - ga <- 10^(2 * runif(1) - 1) - gw <- 1 + rnorm(1)/5 - resp <- tp / (1 + (ga / conc)^gw) + rnorm(length(conc)) - resp[resp < 0] <- 0 - data.frame( - logc = log10(conc), - resp = resp - ) -}) -fits <- lapply(conc_resp, function(df) { - fit_hill(df$logc, df$resp) -}) -chem_params <- do.call( - rbind, - lapply(fits, function(fit) { - as_tibble(t(unlist(fit))) %>% - rename( - tp = par.tp, - tp.sd = sds.tp, - logAC50 = par.logAC50, - logAC50.sd = sds.logAC50 - ) %>% - select( - tp, tp.sd, logAC50, logAC50.sd, - logc_min, logc_max, resp_min, resp_max, AIC - ) %>% - mutate(across(tp:AIC, ~ as.numeric(.x))) - }) -) - -# Steady-state concentration (will be generated from httk) -C_ss <- matrix(runif(n_chem * MC_iter), nrow = MC_iter) -``` - -Simulate data - -```{r simulate} -# Simulate age based on relative population group sizes -simulated_age <- simulate_age( - age, - n = MC_iter -) - -# Simulate inhalation rate using default params -simulated_IR <- simulate_inhalation_rate( - simulated_age -) - -# Simulate external exposure -simulated_exposure <- simulate_exposure( - exposure$mean, - exposure$sd, - n = MC_iter -) -``` - -Computations - -```{r calculate} -internal_dose <- calc_internal_dose( - C_ext = simulated_exposure, - IR = simulated_IR -) - -invitro_concentration <- calc_invitro_concentration( - D_int = internal_dose, - C_ss = C_ss -) - -concentration_response <- calc_concentration_response( - resp = chem_params, - concentration = invitro_concentration -) - -concentration_response -``` diff --git a/README.html b/README.html deleted file mode 100644 index bcee9a1..0000000 --- a/README.html +++ /dev/null @@ -1,553 +0,0 @@ - - - - -
- - - - - - - - -You can install the development version of GeoTox from GitHub with:
-# install.packages("devtools")
-devtools::install_github("Spatiotemporal-Exposures-and-Toxicology/GeoTox")
-library(GeoTox)
-library(dplyr, warn.conflicts = FALSE)
-conc <- 10^rep(-2:2, each = 3)
-tp <- 100 # top asymptote
-ga <- 1.6 # AC50
-gw <- 1.2 # slope
-resp <- tp / (1 + (ga / conc)^gw) + rnorm(length(conc), sd = 5)
-
-fit_2param <- fit_hill(log10(conc), resp) # slope fixed at 1
-fit_3param <- fit_hill(log10(conc), resp, fixed_slope = FALSE)
-
-rbind(
- "inputs" = c(tp, log10(ga), gw, NA),
- "3-param" = c(fit_3param$par),
- "2-param" = c(fit_2param$par[1:2], 1, fit_2param$par[3])
-)
-#> tp logAC50 slope t-error
-#> inputs 100.00000 0.2041200 1.200000 NA
-#> 3-param 97.63578 0.1516374 1.485236 1.177891
-#> 2-param 101.53494 0.2360917 1.000000 1.417567
-Input data
-# Number of samples to simulate
-MC_iter <- 10
-
-# Number of chemicals to simulate
-n_chem <- 4
-
-# Create age groups and group sizes
-age <- data.frame(
- AGEGRP = 0:18,
- TOT_POP = c(0, round(runif(18, max = 1000)))
-)
-age$TOT_POP[1] <- sum(age$TOT_POP[-1])
-
-# Create chemical exposure mean and sd
-exposure <- data.frame(
- mean = (1 + runif(n_chem))*1e-6,
- sd = (1 + runif(n_chem))*1e-7
-)
-
-# Create chemical concentration-response data
-conc_resp <- lapply(1:n_chem, function(idx) {
- conc <- 10^rep(-2:2, each = 3)
- tp <- 100 + rnorm(1, sd = 15)
- ga <- 10^(2 * runif(1) - 1)
- gw <- 1 + rnorm(1)/5
- resp <- tp / (1 + (ga / conc)^gw) + rnorm(length(conc))
- resp[resp < 0] <- 0
- data.frame(
- logc = log10(conc),
- resp = resp
- )
-})
-fits <- lapply(conc_resp, function(df) {
- fit_hill(df$logc, df$resp)
-})
-chem_params <- do.call(
- rbind,
- lapply(fits, function(fit) {
- as_tibble(t(unlist(fit))) %>%
- rename(
- tp = par.tp,
- tp.sd = sds.tp,
- logAC50 = par.logAC50,
- logAC50.sd = sds.logAC50
- ) %>%
- select(
- tp, tp.sd, logAC50, logAC50.sd,
- logc_min, logc_max, resp_min, resp_max, AIC
- ) %>%
- mutate(across(tp:AIC, ~ as.numeric(.x)))
- })
-)
-
-# Steady-state concentration (will be generated from httk)
-C_ss <- matrix(runif(n_chem * MC_iter), nrow = MC_iter)
-Simulate data
-# Simulate age based on relative population group sizes
-simulated_age <- simulate_age(
- age,
- n = MC_iter
-)
-
-# Simulate inhalation rate using default params
-simulated_IR <- simulate_inhalation_rate(
- simulated_age
-)
-
-# Simulate external exposure
-simulated_exposure <- simulate_exposure(
- exposure$mean,
- exposure$sd,
- n = MC_iter
-)
-Computations
-internal_dose <- calc_internal_dose(
- C_ext = simulated_exposure,
- IR = simulated_IR
-)
-
-invitro_concentration <- calc_invitro_concentration(
- D_int = internal_dose,
- C_ss = C_ss
-)
-
-concentration_response <- calc_concentration_response(
- resp = chem_params,
- concentration = invitro_concentration
-)
-
-concentration_response
-#> GCA.Eff IA.eff GCA.HQ.10 IA.HQ.10
-#> 1 1.316500e-04 1.316519e-04 1.197871e-05 1.197851e-05
-#> 2 1.928066e-04 1.928077e-04 1.756148e-05 1.756120e-05
-#> 3 1.509427e-04 1.509448e-04 1.339909e-05 1.339880e-05
-#> 4 1.931396e-04 1.931407e-04 1.742350e-05 1.742319e-05
-#> 5 2.000661e-04 2.000678e-04 1.748640e-05 1.748673e-05
-#> 6 9.090331e-05 9.090415e-05 8.634445e-06 8.634337e-06
-#> 7 6.336773e-05 6.336758e-05 5.870044e-06 5.869961e-06
-#> 8 1.501012e-04 1.501024e-04 1.365014e-05 1.364991e-05
-#> 9 7.009742e-05 7.009701e-05 6.962928e-06 6.962838e-06
-#> 10 7.310190e-05 7.310279e-05 6.848001e-06 6.847910e-06
-vB9#ew-{kuC#E4_cig2&vW&Ebz<%YQP?pf;H#@VgIGcl#2!%N|LU(wQ-%9gxH4( zJ(;?o@!U9SeA2M~&4kK6?B05D!{zy5SmV**>G{P)gMoL~7IL&}XML&HG4<0nxQ{TJ zFz#22WFNOfjE)wm#L3zyowv6`Ovu`v+lQD?nS=@Tg_wZ*TiZ=p;aur;+(<13w@zij zerh3nK`H)#y|%&q)VpAwwQP%X?Wp}`<)ruZ^)t4m%E@EL)OamhwDC#+pe{Z=I9T5h zS|nS8+?gz{S-W3zu&&+rBUNdsuXEBCK%GQX{Ww5euD4T}Tqu>dako=h?0ZvtiCPGz ztbFst0G@?CLm41cdNZ88cBWC9 zxL*NUtrxT>|CKeXJV0Y^X^&P`6s}ByoC{a@B^8Iw^Z5%$z vwI+g`&tqW;K x*~fSs3Q~V~kw< z#oQ5J-O5(w3NLy1X{P-w4Nlk?KTWvj(S>OH`9f&~tz>s6Eq=Eg9_BJ51cj=*RpL&3 zg2bImIaD!bN&L#*^B)_Rw;GM944QvDr@4++u}kj7`DXOPBZK)PSMO|$i{cH1a5{R- zTG^WXlB}U9W@;&g=taRrmRp1t1=~X;sNW(O&%nqHZx3>9-yNDdyLd*!?ivm0{jn9@ z5x@AkcKhgEojxaN@XF%$+1Z;1Uc3$L@9I)XgFmIrprEzNZj*1jg&oz&>|S__H(ahR z?Hr4~Q5+q7&{Z+@qv4I+lRH^U3m5K_m#6KE+ec?!95TV5W=_L|+NtahX`Ws%P9o05 z+ZRJjLQM8AC=scORlAiRR}hsSf7IEBafg9EZFpQi)=mccZW#BW= Tj{#9~oQ=$4Yyge& zb0*w8QTsNG%{J~`Uqc2%AS!Wp*?rkep;Q%kj{CY4sN7;BfcePLN69;iW9Nq2Xy+nB z)zvLHYtNS10B|@05R^b2BSF g8?|n!7)?2nKNowGgO|7j0J~1ZrZ>`$W6*;*aoD z0gn^1wp){cc$j$jq0V{pQC+DV3koPId-rDkolRp1D+nR1`QYwdlM$5w3VY?`w7{gm zWaXF+q7uDNv{OxdaFOt=7K4^%?8FmP=k$tTVI_)U>6!V}QUIExNd9^*GHAWHCTOi@ zq ;ZRc)AWfSV;ZzRzk$-lF4 z2Y#ip&DTCySEScBK^!(An+9$eYGP>u^x4mk-?*oKG6{i!OvFm>kBNGS*aGaeUFcD1 z1fX7~de_qVi>%Wj{3 Y0rjWZ`_3QB?5%fu8hCjJu6H1c>!Pg<13^t|Z3tG|QrV z+xU&yMSk%i>4y%dRfec?^f5jtHGpW*rPFg4)Sm@-W3JwXL3=%;@{+QR%dMS?8kel? z%6Cpcc(o7vMJE^7u5d7|!SV6)gF==L?V8Y| Q#`H=c3)YYSIMYs!5IdxTTyavmTr7hu}6=6Ft zF$*KQv*>u1ol~m8Jap3MaaEHD7U|juBWijbCJ~i0pg>!P`=fdqSEb3m{>fphO|8KK zST8J+jfU8VhzEjEqzFb)B^X5qh~P1R0IvUhRCJgn8{6l~@2;17jaz&zrH}Qw+TL1~ zD`*}KeEfo7&^+JIeL(xL0BA0teZHgWHvv)kyUW$fd5-Yov*%K2Xq&y0zzw^dX8K&Q zI{Ae*5HEjiBbtCw5HjsQ{mpG&&A$r)DAQ1EZ=kRO=rMh%ynVbZW_}T9{gZ|$h#Ael zQ!;sY-$Db+-}#Am=pUl}eFtlZ1F`98ok@to@$T &*)zXA=0>es{0Lr`hw zQ9WsAJsYjfT?B~S43!sGcK8UIxlQQfHdu|-ASw}6h| dxCy1xV1z@Ez(pyT>WjsrChwb34iLaoy?ggqxq!E6GP0DcS8iE3;_KGb2q zI4kJX^i{kyzfNW@*M5&Kbq+LN8~1 CtPCjv(w4Jn@Gy?x56y0T3h!(buoAgR`??P0$_pGz~sM!L5Aa6}^4H8*g zvz}XA@cXn5aK|17b+}nsGUBJqI$P z`T1|i@~RwBiJ>7n0GEEO!(2I=GXtX2sfE}?mm$d{NIp7HLiYP 2l_CdRkQVPz&iOXd(wQ;W%9-;$-j-IKjAQL1z0`(wRrh2 z!K0CawmW|c$QtF4e_!X(7`puzW9J9(GkkwvEiX?Zd;obw<)LC_Pvu4>$t1#;w~r^L zgEy`NQMG!QQu*t0W-sTTIK#!~e<-s5RA@f=Yd`Z}g!KT&(>kNS6tVw4`uwIGfQTCY zrtxf?`L0mDeUsc Q;BlQlN43{@ zUnSV3fF&!MvraA8KnH^jbnx2ANcn7~XEKNh|CIby)o&9$Zh%J1KPeQBbP(lplLt?B zyrce3hX0HQ`L M;CVY|J=%bQuzjU3D>k-=tw2 !6oqpp1Z^7bET3OHAVep@*F67inlWcplxBI^GU)<~)AO#TWq zV43)9_$KfA&wH4?GO%_2$5O8G-KT#k6aU>auivMQN+tWxyPzwlRayd$UQD?)FM6>y z&*37$XPk||bInh$J~hv~X>88@pq@K4xgNqv^|a0&d{sEDJ6h}uyW1wb9cFq)e|GnC z^pE&??u)u4YNf~{_}62;a7LA{-PfbJZ+16-?WAeVIERT3)g@9ZG=KfPV=kH~6U}|B zK(pQsw#i;s|N8>IhjF=q4 P(11K*c2!<;0?oAIfat5NbQK58M_LYFnRH*}Qrb zc~bZ%+rHmMDGBf=OJA>1tl#UuaLfNyT>Yuo`h&3C?{oD;@sCE)Z*$+P)!!ks=x+ia zsRE*@cRPB#XkOZEpFN)y?A{ga<0V9eh!CYiXq^Y}@$o_wkVvUJbue@WZgd&S1-CbM z1J67!IK!>m!;Ra+pXgNtQ>}uDLs}UXOA266uiH7;+GFQS-Dp^Mb#>**&oA(T!xg3v z>$mJq)p{NJ0RRTbvE=7EPa9BgI^|VPLUK{i-vR{r`tHxO-~eEVlLyb7p0@M>7y$ZM zxca@#l5p31>XV2_b-+*A2Mf+A+Mt^618`j%${(=Xd@lhI_z=MI^KZLqwAa37K~#Ru z2PBy!(DUQp7%;;n05UXdxRs!*HZ6c9fyDra->NNGobHVkm^6h0v@x`(aQQyYnEW*0 zCQ1d%`{ddpJODinJt{x{mEQ37<_atsQO4qf2g5^Q$y$l;gcsSU#$++;!9U#c!KyKO z37D)I`SxGe_Z-!8X7=1KveDW3J&`)j2+~N$JahlMq*4a10TZ#>iqj)V096@P4i0}% zymoZxn59@bPAq}RfC K`rj^6N@gYzs%h89{+c=ex0CR{vGlm;*L}#GusYU+xtP1IGKxw? z>nV}Wv~{0qd$BE~qR9HHhL8{YLj}kGO=|mZso4IX0{>*m`(L9#Z;05(4X?gZw2lAS z0{*{dH~&Vd?CFDXfKqYrKhU%Mk8Sx^-GQg?LFV-d&;0*W`@TP^!>yu2@H0+QxhMz! z|L?Ibjt{p+ceYCO&o3^&nphGV)4j(kFBR|=a7#~?NXA5<+L9KfDf;cd$12sY^c8@5 zOrVeBum3!6_W8>%j4hNB$&mZtL?k$=jP~b?pJuzSL` m4@80~h3V-xm}yvJK>OCp1kl zlIm(j@RsOqA9(X>A-uYu`R~qxt6j-=MN%fZJ90$>68 J3m5# z?%Vi^$&XdZU%bg6e2HoLDfrqgI)73wP8C#6_4-c8lTg0Yyecfd zL^MT4zA<-HTLLSgYYBsLmuBbIunPzOBn7Kf?Me;*0=djfziK3M{vDH+cs07gGfZSy zZD~IDs5(*O5k0 5x*Mu9{fOY3PK2}y{?m66NPPb~jMHY?AgJ`hoGu B{MH$o0Gvbqv zV`BR9;^QIxjPWrw6~He)#6 CNH?b6Ar} Bisz`}q6^yA?KBtcePH3{-%V?{oig(WpqZOnOZ*ekJW&IeT_`H}>lta)+ zj6evawMiro3Z;|tyFN&N$2 HnsgnI^+b_V#1bSB=5nbnLY$l-Pux{4q*+bjG2UF>Shxxl8BH@ zl$sXi;g4hyDF}x-#%oF2la6#NKk3k6a2DTFguco7vPlq>*38KM$n`@#oET;%#2P^N zMH*(9JiuIE_~k5KnuIZE$VTm)F*8fhF7S+{UO}O8z&Pp)6$?xz>l_k~Ed{}g c$By5dT1BamOVUqL(`_mEWa&nlCou|Y?Q+0yGwA~ z*N%r;WmSc!7d7?;0m`k$2FuhJo+dx+JGD&WhSgIWzXfS3BaveK%5%xjgJ6#bZ7ZnB zRiJaaEODlVvq84IG;{)tt&hAeOx3PsS;tXe(0JIW$hYX<5IiOsC?2pQgjUDJ#U{lT zTjz_!X68Zb)YIwkol?K<1;yeo9ErP3>fDY?&kK7(E@Y+9-Q7BjLau1%FLe(qm@0es zixYWxM)-HPhW2-pBDb17YaWMQV67rh(P 8lhs5?pBH?2Bfs^wysMk>UbUA7Om{PglZ;V# z15(un>l{Dx5V6CVERfC`W=~-7Vsy*7Vc|JYuA{%pFP#rF<4KQnt5be^jQ8p*mnk0< zLI;J?B7L4V(@99;Nukp!x8q|$QShMP0l`$pOUDF%xn&L=9S0u^8y}k-{pPD%n96Ci z7 &OMn56I|l3_$b5Il4`20SU>_rCN 6Bv 6`t9aSDeq28-6}g<{j7Z}tE^ zZ6Xq?sw~NiDJkbK&4A9#4P?e7_(7skcrg`P37PPD{2|ep$RHT)3j*t&b8l(4PsHk- zT`cvxjN)t6Z3=UUMCfEm_&7^wz(=79I(H2#n9kjBt2LPM$G+KE!zVDE?7Isb#PB@z zyk{!rOmGVRBvS>)TAR-FTBiZ6$woN^`cI7xS)|tS`CECb i9`Y9SUmehZq z;E>NN$+OouHD}fAh6Ybv@>F_$ww>CmN9`2%UNSmmWXe5PZ|2A=Sr5`wI=2bDqo!YX ziftOxpGWj1VYzc~ h6u^P^?sDAhbu zWpZ=3sI0n!O|QsM*l`o7e5^;xH>FKg63DoDNgP!^?4ieD5zD2W9*f{=)f-8o^d<3a z+md39+($RiWXqCUhHq6ANt<<6iNZ4LvV~U;syBzfZg_=z)pH+8UX&aiA@eyYFd+Dt z`H 69jLA)uAX(&)x@bJa zw_d~v-E%5xl&yW$$1B>6te2uA9dF8+ky4Q|<9=-QV4pB>N|Qg =bwg$T|6ZRp9 z<;iv3EJE|gyZDj)>d;Kbfq=OLHU(-v1{^)EpoINb-9f{%LtjU2T3IUa&0mK(AQuI6 z%Zg6&p_BrQF>t6PA&e0UpHuc@FVxJBH1>xuTG69-!d393iB#~Q!(a}UR*$WKLYVng zYgDainG$G&AhHniJcSl~LTNZ2L3)HA1564trpXT-FAa&r14ASVU4rvVuD+M^>Dqbz z;e-XnwRI6#9D7pdlE+g8n<1@Fg?YRPIjqD4NQ@62PSz20z?c4|FEus_?zuV?OJX^T znOiPMjwk`y!ASSwwP#lQW&z1NX&QEcuZpm5C+gWwu{N*rXNe2DSUXDBVH-oL0(>o) z< Z}N G5HJ gPpKkIxuA4rD)a0wed{Baw9U0)8>ThCSR$FQkKT|lm}x$443^&(pDNB zI`4uCrD3i;=HxWyK^u2lPG>w-`XMIb< (${-sp99{1$^;6FG9tE#b^bpw`9L&9J`@(J38 -S=lmeDL z54}IwQJVW9p%%Mu8nmUcwAYWq&on5~y`ia88HWrJpIwM0@kvE*IqgS1KARU*iY?`r zSu0{zF)>L9FIml5*c?#hWRA>LxmZzSWUclHch<&KKaYV*LL;fBFHXNYoFjU2@j@fH zpfu{~5FB=X#wFIvd Bk$1ukT&Up^OXQy}v82ZBei*FilWi+xztl z&Re^_?a;3vn(gp`d!<7?LsZ32F(64{XyFZ3yN0pGhrqJel{V|lI*Zb*$@8q8CfV|? zw)qr$k 6v(~49VNX` zVLNa#J|qqA`H`r|c$P)h*0L9jrD5XVGFR|588t*(%Lz6wWJbxG2K!J{R 78ldhb@?TI0$m5?cTDXn*{zIndXLh2GV|2EKkpF1MX7c| zh1unbN)+`{xkB}p1 ?ty#hiBxPWG19e#} zB_1!1T=*w@;)hf_xH&IW+jOkmboK+nl|R5?Oe=SJ-Z 0^=n<-!~isvN4q*KZyvWW3mToTC$cb##xIS@W2e z2_*`P>V)J6vj_#sP8H?=;O`%T-!&EBB|nRN+s(&}L=tkNORHKs%uhM;u_VXdv9Q?K zJ}c0C)LpBJ{=5UJ9vmRpX8Jy#Z{bcw!mC&J*aob#L@(k?2MTqb)_C&jBMvN9-sS2$ zg{qEw_(OFcqJY1Iq7t7H9)QyWSCdVlbjTv8l7G)wZN%-QDJ$`CgPvCM$s$KAJ`fy} zPK)o*WgRO)&L}KJyObG_*>q-37ZxMgd7?wCn1dx3)7li%DS7~}sa~n;VJaO=9rn^6 zKYi%YRp+4b22K`aYaSI9?$DBI#8snGC<=FOU$s|_?)>N?n^Pb{S;ww!RY(x2KAUbb zd2iyO$_2i3;sdI^(X!mK6cd+B(N7B``xL0UEYI_#-Tk7*n4My@<%F%hLz0b$HL9@A z8ArEn5_=zHKCV|B@MzPL9-2|DQLq`7q|WS}@Zj88y<3KFO8D!Nkb4qm(ROy_!_3jwl%CEw zR#pY?1^nMPywAN{y0i8%Z{PzYmLW=|+4%M%?v49XS}Lv^=`W`dib9<&Ay#6#IyG7H zaHg;QZ2N-W#})6H=;=JGP<5$9I3yH*xHPpiAcpV>t{V5wza4&>g8G=DU7T1D`tiaO zCINo2bHi8J;#7Cvk%MZ*Tw$Pltx(~x;$f%`qux54Vlh;;X-v9nih_Si!0!1>#`r9s z4ZLr%PDXx0re?jDDzwCRg~|t4@+|cTRh9*-6CjZ8M`0tt#`7m}o)B?8HI!B{8h_uf z^0m=89)TL#fYtQPEodEOB&wycMmn;8gmR(^tKy7B&yJ!ZpjK-;7I7-l=t(H1`eb1e zoj2nM6QXx>My(^qT0YfqL4CP!u<=PpP?_Snp&y>F2lxYtf<+xG(m{+ZTEq2Xd?j+< zS4`foZre@z m zq??B>rP_8kcd@iIC~ThHLtEZNnj>c`IVHZDbjhBhr-hhI{1$S7>>wg2gEsK&>wYw8 z<}?GxTH>Ir`x+#@cS;QwyUnun9ipzFFfTOW)$+`yC1o8?`i=>LMN?^L2j49(!M~IG z_WjHTfez^%gL$@22KAH7(D25UUG ObuB!YA)tJo1Z#ap*Tx0ZO0?&mLoDm#J1Q+~#fDBm4J zs=65N7x|J?BPe8mgFtmdeCZ%JPoTcuuW-?1i*JUPzs WIs#j#%i9knd}K6{ZQPuJ}0-JN)g<{&zyt20`y+; zle?kYcC)j|_bd=$;_UmvgQ(}zrndPU+O>7Yq&6-G?;GdH%Qw|-@aBk)UOO2PSm!HC z_GN-{(vN=bpQ`X+XGmuhBQ@46_#8i%@I|{UnOS(M&GYS%YoB9$w^HKUb&bI9rQ*d8 zcf%F3+?7VLg18Fr>nqR~(d(4#sV9hxaInMo sxgtgDkFr=4l~K0ZBGB_q4KMrk17i@7LyC~avD8uurWoyd<3H8RQf)(CB5 zvNZzjtgOGYS*{OC>6Wo;W$?d}7$zaj>UZpA>=J_=`z24>F6v;zSz|tE#AHy!C#Z(_ zG027LwJ`f8tr#YJaJzX!&S4#w!zpuFrb52E9Z0TLGMy%dDM1yN-ZDERt7Z2HqOrcY z;D6ApART{unj?mo_RtW2|Ikyn8#{&6mWukL$WFoQ3YN#Wa8=Q+Uq|if?d>VKb6xh^ zzH4ND^TOtjop+)_$(tR~jOJ_9G}#u(=RST<1Mh@jLGM;+?_oz_vB@&b3HfFEc_tg} zv%JiN4IrZBoICOv)rjS5by(og_qT7{&(cj4=w@fniDG
sxq@iH%~5y?N-XOsV|AT*ez#-H^wt7lh3E#LUbg z=yuj5yEMUv91QTg^UN$aEVA@;k~d+wxqbeH6*c06;@m+QZ>3`)M-!nF`Gho-a;_{a z9LL(8-=f*V??s-wkOjFjn*|yLEh*ewi-sK(Qq8Wuw@oVUGSD-Rn<#25Tb@Y|liGbR zNYxSMt#fUUAzVKJLz{>GgLhp2NwEeF@oPRAss4$H-c;U_GFuZR9p1{3O8X+CqLNT{ zUnBAul>m8qnZg6}NaBYKLg%R=(keS5gDeFRuLTg}&J3dJmDjWfmdoiEuw{apm>n7i z*Od0?j{}2Nu8Z+H$=Ma~Vm0Q_2@I}eY3CbjSU)Erh?G8jzN#52W*h#YHdbHb94$fV z_Ut=@F*nTfiYOES1y(AS9p72r@o3P=z~+>Lv1jj=4%PRCsrv`cckRL&+Shwui_}eN zq~3>lq-}4An^&$uzYi56BnK *4`Oqwm-CUq#bbxw_0;PcYn?6kVHjnh$!y~|$py4nYA-r*Ph3Ej=CJa3{BoQC2V z2eLfjWQ?~T+_^|R>Ec+Wc)<2$WcyOf94q+o@@QqyfbOZx)3^+k)O8NUfNY#SRng;S zMu9ymr%*GtkCQI@cjonAkvoc7(!#6r19|Lk1sE+>a_hhY_?H?lCs#t}a(o*%sPiJ ?nZ7X%Cd0;49kdQZQLR4Fof%w4}tKx0dM>7e0 zZTy53M)$wZCcmag(nvI?jcPF~fW4FpeCsrSTNs{j_MwYqU9M+d-$lPq>-!}Iu2#OF zY&XXkdE(bW3i%J#S^F$QfxSvU#E0+TZ7(o6>tk$RL 6&2 zF;`=X4SBqR#?76lj>$#@8IV+bhfK#%o5}-wgSfPS_eCYq4yg?98D!@Ot(cQ>;Y4To zD!Q-5i^zu`V_DL7>fRYOqUq6-(|YtJ@M4z~;g+eXq|kG(7D4%LFRRkw_6Y87c6gEg zNthCUm{`Ds~MQeMpsaxE?WjeA~?)%R{^_dZ;E7=HKe9ke?ZW8ckOc+SV45*_LY zuRd_M)ZQtp=vTg-%Kw&!@=Pl)o5GXW&S_77UTSjJD6X+SxJ7IC1s@!{+O)_5Uw*H? zH9qj$M0t%A%i>@J@0uQ?c;2Z1zhaj7x9Bs*qbc<;>Bg2f+AT@QWS0MrqpOZ 6<*0JH#*OmnS^7#&%RC>J%e!b zqy_(lJq+c3_B2T<_75@h^M$}M$Kd()h0J@~Q^XEW=vH%3&}rE2&i8#rMezL}F d>D4|4Ojd*c`McBm z57YNz28HNfPpEf@-g cTT!sm(6qx8JkbN75I2UbFvetslT!!% zlf!tH#csj*c^V^ xvnwC~%`;`5Kd2jkHQEIya6o+|}mhN;W02C-3`LCK$ zH^K&ad@(mm7{O5QX61|WCgO@-%lHTbRQ>0hB#Ns;-cK7lL1L=-e*Kssyy<+5K&bJE z&{%LeKM0wMi}avKJTJ}{`Rt xTYSg*l`ti`ng$2aNv`f27i=`$9Cl3)t3YgP@(f;g zeN?ESv6)Z%K2Ii_`uP>OgySitRe^QRot>TIH3|_~!(wMV$hFc 3ISG|$a^h@&C^Nfhn;2uUly8f*!(eN~hZ76`KOWZ29V>lZ+nxyZx z@2ECTv=JmU&2iVSMV0;L&D9$dpRi@?ZMX248WIwg%=$@6e?AgI---rlU(JDe-~1g9 z7b$ssw=)|1tZ)h?Tf=ziB|Fh7;1T6stY}uVqBBTsmx7a%rzYryusQkKV8EAGaqo`X zCD9@}M#9|4jiEI_SXV>mmVdi54$b|7q7;uvP1RS{Hf$mB5jdua?`&N_FY>Xp+VTKA zj|800E9$IGF}6jlQtPgK>rpc(BA*sm)ud%qKxQrdSxR-5E!OFgkF#|E$6>n4 R zgHnsJqXsU2=W#;ajDIl`q~-zx)DC_HS#87EN=2~jJX1%h%$~b-0y1glpM|K~%qk9J z0i (bi8`c?ZNNbzSuMHM$${ZChxTgJl 1mjTOQbv)p zrp;ST_ToRyX0y%8V9n0i+1cM_Gt>GA=XK3g*@f*Cj{@qBPus7;O@sQdqKZ-`?c5g4 zcf!XEza}!Q!tIb0OP~F*iDAW8Ph @_3xjx2Njw6pVV*fkblX=-}d@8?b9Ilu39Hwl>zIuI(g*mkf*8ivE&Z7VW3Du!a zexchTQ@_S%j|vKt*Htu5zZ$k0ip?`Ep^QE+Tudvn7N0vjq${u|7#*)X&T3JY3<`7| z`Ze&9b-?sibi!F#iA~>Ei8xk_r=R`l6q7 z6f1pUvPe!ANyY?e$tlEG=`J> Y%Y$l! z{rqb~Ko*lP#NcNV>_*ZdAsH P69;)9w)XblImN%R_ zr6F)Vn&ChNY!q!SI^}b|_D+*hs`OOI>+0j{=HEE}yJ5 |V4MIV_q)Gi5- KHs>H2Ac9_I4ZE 2 zLPZYgTLNa47zD{nS5!Oi+g}_uC)&*s1{Uk@GFyKqnQ+6)HC)r`6SR#}%_N(CsyHyr z_cjsW6~7C z4*fKxsMDWTN%X>mmgBllT+jPFJ~tK9TE<;kr7yHPe4vSTi#w_>jgMeB=EU$nsW)Y^ z42^?|o_^+|aaxWUPUf}aSf8PcTRDAOh(B3Q%CG)ZD;|1n` O6PDvl (B4|0P GFq;bV05BrW5x8O$GOc@g5NlSfAu7&2yvvlvL}KJg9Z*dqUGdBF*ZeoS4^c0}&(B zwP1=f44}xh!n6D}Hy{RNDBF&jm0flxL9?!VyNqa%-dJ-3I=n3Xm4^*mKuHo? zU4TFgB<8Uk`uxRB(L{ CYRpwY8=l#*!e=bh%Ae#>Bh z(23*%UiAlb0(VLXgky}D7rqTe(Er$APHe^UNF*6FOq%1g-=s|FRrxwtfCCo^%Uu9S zDb`09_pW!??~gs=r;%vny9PUsYFGP?$T#)K>A@{Yz{IW+Cp`XFJTSTsGhNCQI9vU& zh0H?BTmDS7iS|FGK;g41e*`5tbiF#9EpHWm&%I}8dKcRu#Khy>QPY{KuQ~EgK;l_I z1Kl^^V;?|x{EQY8kq@bbZyXd!%8MsCJK;czxp~%DP+npuo*N-__X}lhnCydFbW*2^ zq&eZa(Mr5H!b#ZQVfl^3?ioG^=d_utD{>iTgM+L!0+Hzpu^4OA4patXU4deGcp C*m1z432AR2TK{Wu5~3@8kuXTwI~Q=4KblS_ z)kMt%?Y5Ja(Z7#dJI8(ZPPDK>gXx`7;7iDN$T!M3xu)4qOvUNel@Kk`!~^n`2e%%D z8#Th8dOnH~YxL6yUQ3vrbgm@rIC%Fv-}0Td{y8NUy%B1ns9jUvv|zG&KxkcG{5!oP z=^as<{ACR0NE<}5cS^XSe*K_>A3~=>OD6;Uayc6FQ)_%fSNA?Rq#9loJu?JsKj>8P z2rm)Z?*)t?sfx51LC*N;zF7!JgF&m0o|3(tT_DqB+W>3@{LY1Zu35K;rxc08393(2 zTn0n3b @1S( z5DzjtY-z#dvs40B8(pKmYeQOKaH?;A`a52yuc7TrMS(VhnXby3wkf}80o=^eE95-u z&3AXN`>Pd)=gs3@dk*^N7Yz;v=A*Hm(IPpNH^J+>zpc;LHbqjsz4%p>h82ScKJU(? zgd@Y(d04E2nU^CKa6amY%!DxhK{ }7`6bdo zgV;Z*xDpK647bG;37;`gvL)RnG+MT!I(XsaP7J=F-F|ZCE`DO#j_1K>KNzwsXC!ld zxIH7JbJ)5xCGH81nM49UA85YPKw!B3qlBEMB!)`*UJj_SaE#1Ue2b?)_?BV8z@7_- zpEm67sqP#Fa4eT~RtD$F%T1ElU^eBzPXP7g8Vv@8OtzdK2J2#Hkze}-;un#BKv_L8 zBl*%o1Ba(K^jnLAXJpj{*z*S{8bOrnk&X7%i`i^1*vzc$K-9TjO}+!uIAL^=v$M6e z&I!r%i{KYlGsy&1&Dq_1sWP5fZp&kO-wC$Meziw)?P)fuIiF;netN7V+X4Gl`TEQA z``?~^%3xRLb=Lgc#Rb=8k91qX3AsK3p%}iExoK;Ew!k3?hzo}M_%g|D9fsDjAS}dA z7@|bS%)Nj`$i%GkYG?G6@Y&-}IDvo&8nPs ;% z%;C01|3Fy^ `j_!}%`;`1;HW-spVozDyq4%a&HcCLDnvuUTS#HEHG%~GJ` zFUgn?Y|0#gaW &X5>VX2P*A(p#~A1M}px1+)U4XGaKxwpW$UK z78-wXdf_HZ9AdS+#TPaP0O0w0BOizgeU4M7*l#Xyj-%6Dk^!`}^|oo1e5Gkk;?aar zibQ;ZSVgQr0cbMQ5q=D5SytOOHp}_x4;FwN&Ms(zgKciIRjzY{dv8!j1UXSXd+k7E zud^5?HzDK)zIwf){N{LO_Jrt2;rv*9!NcTDlm0}|`rvzr)IY5qs9;z?3z#wBIlmDy zRV*X {5Bssz z^gr)j&A(d1WqGV8S_jC))$!ycR{h)KH@U%gvVp%^wek=1jdpf=vs;YH*KzmW)2VLJ zfrbZRRy{x?KZVZI*T@2Ygh`6BYcy-IfKDOtH%aqg04q(jaU66AS%8+F<$??$o~@#R z6gF$T#-5n3%GC}a9SW1sgjX)H99OFjpy>743Jp+ym|_rB7Y!SvsyerLurzDN1O}KC zGBK@RUOQ_ nJG)%dij$6zu9*F>&UnFXQ|DRXyq^3ORzXkbXM>^RZtC z33xi-=B1e;&t0g?1(a${pkil0ZeWzRPB2ET;a%@KlOO|$yVG1TEVVp4*g {^>;5QVIwEe;z`RcfkPHh0J7b?NaI&+@%aRcBhm)FE%HuQXI`P&c=Hu@FlElo7i$Pv(I9kOX)p<2Tf)uohSeWKDh9TB%p1H%t7{ zcv@QEk%L{33`bjHMZ^c|Y+FI`eGaxZR`wdOR2(SUDTADT2CaxjD3Z;x4K-_mGU-nV z*j$S;c3-jtLcLZr7-hW8E}MI{d@iAOt5liO3^3eJ2YNaDK9;D?7tt%@;h0fAb3Z$7 zv_4H9?X|Tc!JCY~W0D}GzJ4B3J#jp?(ROR T?Q`pD z+c#fV_mLF@*B3wfQ7%KVgN9Hk8=JM32DdgNNuaL8T112?fxZRJBbSPGH!~~!&HDM3 z3BXgK&P0;t*ro19$_>xxf{n$s(}O#0&=Z$Ye~W`FwhL#y^lse~yz};m|9&nB;HO;} zvl@CRuP_$+st0I?waI;3`=ZP|@$BuID9-b`7i73$#1smt>8$cq@y8_*2jcH6-ymuS zX(xDnUp6Q%?xk6&2h$mP-_?eSB)5Qpk`R00A6eQTH0{N-rqQL773>lRveD#($~n&2 zG^z*~4zq)HVwV(C_HNtQbWI{q6AkVH3#(b#`S4{Fge>T}kU=i|Rp+LQw)DT&f;=jj zkYjZ2Z$I_hYQL5Lo79`Kx?}6u?^Z8-`<%`CO>9Piaf<7Mm0?w> z0c4LxX-Uvj&IL-w6f}5nG&%VnREFR1lF&=rO{C2eF3l=eGN(XtX~E(YZD57=m~>Eu zyW#x_(v#`{B(gH`N5q-fQOsl*Z)M4mZ_UMA`dELR+=&ounj>&iq1w@yIdg?4_sW0Q zF_(dwQz$(%u1t~fQw3&K>CO?)E-5^;d#bnB($dl@RHA6MXxmXcE8)V`-1Yy<*O>ow z&3AD*Y0Y4KiR(tAQc$4!>HS-_H?DwgyVBgCIi!PV;^>pQC` (FewOnDWzp%V)R7g&0ZVqGItc@vsNu<0{UW;9pBJ4hPlCCzQ(R6j$2DqTE{_ zZwULEGQlIvK4 -8euO~ZSZTor^P=!3%wr9tRz_zG z{h~jsaIeHkjdDOpT=Fx;8425m43>Hxg1cXyvUT$o4e-vQWJs7Nb0V7@FO2nwj42z0 z-6vJvq-64|dvvbAI6JlP0XXfz#QSr017~TVi>`}&!ppY;vbGX;UfLm&vU!1OukDU# z8>!f)TdwgzKuPCWxftIqFMr2AChtwU=6#ONcgBU$-l~afPnu?fURTL&uciqytQ~@t zVx-xl-`(zTB4C7mU&n{_<(P>^1z&XT>Q8HfU_DuGsYSV?W(M&}7~xs7@sR0!eobeF zQmtDvbebt8Dj2a^U;rhCVBw{XGVGligA f&B~|{m=8tM(>S}ODTT&m?)`)+3`upX4 z^40}RnV>B7t@``3`327q4}xu$6|U)DAC(;v5e9&gQ{N zS(>2+ $(0T Ja|sWPKTxAm8UBzL)y@}(}s2w`q7L4sKZA61`oRu?Cbn o{7DgZ1g_Cu1hzfl(CPn@!JQYK%jd)B`$ z*iNjDooSTd;jMt!DYDkF0qz3!>v0c3>lNq+Vc8R8IN38SsnR^YCr|r5434NR6tQ2? zHo;D_(bh2c=o8-PmW$U(GGn5VO}6~igkq9ggUU5~>2{T~HGaYLK*)>RPduWKV3(wC zaSz|N*1}Gso<;jaFasc|Jvv=P{w~$J_$bN3*WZf5d4B*dVymDPf!MdB$bf`H_zuHC z?}$C4pJ1znqXA#p`k2q&VtF5QqQ=v>1)z-_rh!Lhq9JAKv93a`dnL`ciujFW6M5sN zQAvM+Zngb?;t( tM-JPlHs-dge>^31|FypzVqxO&cf%m z8(908xx-ARlc}xOdRNLJgtuxn$%?zFs3@1$3bvl2Z!*H-Q9hu9@3-qqa$nhRFHza4 z&>Qf6XT<&KM}pJsN4kx5AjM;5#@^WG5*vX UIu_qm#R(t`z2y zh0}Damz3BpGUG;@`OiarAm^zkj=xSPS~ub6?hoqR{Q7*Yd~ADdsUs40;}ttTDiuo| zH8~94BwiqkCi|+{>}; #|!DQ7P@Q{?izQl_RFb9 ziPHTJm+slm7@<~2_2ki^vyg$<1caGhNH0!`Mpt-BSd*yRj9x~B1d}!B)5j|*xgXu( zS?)zF%B|I8jpVr_V|>kdYZE$JUp$-g+Qyr>K3bE#J-id%Cdzk5=Cs;)A+6irRe;Fv zqQhss0vCe0sv!eZuf+DKmuS9@6M E?3RxE0I{S2f=d#Go>R=v1C`^9k0Bw%R`ikI?42 z2`~G~b)5YN7=Z1t-gh(#%p0-QB}CqNlWRzns6Ya8$i#+{ll=;*oyPCALKgM+wk9K2 z?$5{jmMI3Wt9rh-#6hM%OI8H@b&%q|eVR>YSDq*r^atlGtos*FJPwL~{(Og#v aFeC2vn~zw AU+StJ~f80KOpHR<#TNDWStia#-%9 z6lRAi4(fOeB)g2W>V0ac&mGuD_?~iG;2l?SB`Xa@w3z&vNFzlF{jhiMpa_$(Tg^Tz zFk4!s-SZW2v<#k%5^fo=kRO*}9G0r|6hR#S$R{Vy43&pF_F+`i#j88Ll>Av|lO)`* z87uyn7J#o7wo$o)mI~|JS9j7tM8lJnS3F0-W|a-1?Oho;x&Rw6>D=G`gCiefF_GKT zS!D+48;sNl20ENOF4nWO2IHmNB`v&5|2QS +5v3Oro&b=>${W+oBe&k>M SY$V;xjf(b~JO+Z(SInF-(dtXM}f$}7*q#cIhb;z}V*TVU~3ZJ8W4vajInPQI1 zKxAR~Sc0@dZ_XYAof!h-45N1hM)awW9P$OMZD!+IF4?&G0Xy;S)iVn*yP$pb($Ogf zBbPM?bxC`_RLeM_jugmIkn>QQbovXepDvA*L_b++sq6D;29M#T6^V>E1z17B)I_De z!&<-E?@ja(iszU3je%fP|10d9r$kqg!x|jjk^;tEUCw#&h0xmO3&H)hpv9>C`iU zT!xBwhOk;I9ptd{-N?@jA<*Xo8OURf3=@~!tktZ6aUT@&r)CgcXLq-9hwzK!ACGKs zjtVhV(<^DwA7$G}x<|jZs=1H6gTUKjX67EaqxYFX>BoS;)Ej1mNhQ-MJKgyntzXBW zurZU^L`H3Ohr69XW`_5BoIUyj=9(zqTm=)|`ATNW0GrMJ;7^gkJcnq6p61kM;$A}A zPC@0voUg6KKztIBliVa9uR}KMo=cz6xiF!+WJ+>;0I7)Q1Nm@l>?`)bSHr%Yk7eIM zK2H?5jW1}bpRd&ONDzK?s9J7e4eIkcMULIZ82Og32PrEl{b#Sn@2!cx&iRq=?mex$ zpZ1(oV9BD)R$ms5Rq&rj@RnXNVL>J0(!AEA3@uJG%_{CcO4u+2?jt>W53aDBn<@+4 zT=a`qp2wIfzM$ckmkHpx6MKMo?H6Ab2kdH~Pn|U%)w5k@LMNvOgSmFyJHd!LcF~g- zWXD*xgeBKQ6qLCC`S^Lhd5{;7u{fzu?#7)7d(5pSusFy{M*?90Wty 4%1>{j(;FBtfnn7n#tG|f5RL{!O7>*>n*52E$?Y#<;kmrl~gMuX01 zz&3To^P2fW08VVasjqxnx7EGdsP}GYK<90yk`6Nl j0sjn^%h {;HFqn8%C_H6F%VE=OaT#idz3k z!^?pE+;L?S!QjKGwc3RQw_9pg!M}yyU7`8Q--y)ZyT1lkSpLD`!&{sd3Vu!Tx4KEz zp> wCd6Ko?+k%_7_G14RCk4cks$FR2y zbchR*gXj0Je;LLM`Uf*GhYL(7dFUw+nNssIm@CyKFsfLV8v`!umdX`}wbWeijGz76 z52-AZ-feZ34P2;4VWgLuejVXax;yQgifwZatA`sk`K#5NG}rD%i#-v&DKd%VS}>@W zkPtNpl}QM{!SAji?`frZnP>bM2q=ApqPktaNY8Ki-n%A1GRU9ECWwSb`O*8$R&yKA zMS-PuFRJ GPgz17Lw-Q9}@KH9Kz{pN3Q{_H5g*@%%wSS@6v zM9Y}R#evZ;wx?jl0J5(2g~B`W37h<@{Q@}V#4-Bs2+BdKBRsAm>0X L&zBO&3n>gv+DO+7 z16Tuo!5cLMp92;ARNRK}r6xBHs&JoxnY+X`SwtxKjl0kRN17tLEYg;~c13S>i7JCs zUUE=jzW*Mraafke!%YQ}N&m%k97XiO+rBsb`HeZ=!IonZ-txl4WpyehQzh>no0zSi z%(TPHiQTE>5_g|U%~m^TmYt_Y{PN&k)=sfKu1nz=ow`Pw15h