From 4a7066eac9b6c88a8ac7d3f27bc65c5cc1234c19 Mon Sep 17 00:00:00 2001 From: Adam Howes Date: Wed, 18 Sep 2024 17:32:48 +0100 Subject: [PATCH] Issue #246: `dplyr` over `data.table` (#329) * Remove data.table from simulate.R * Remove data.table from postprocess.R * Remove data.table from observe.R aside from filter_obs_by_ptime * Remove data.table from roxygen in latent_individual.R * Rewrite as_latent_individual to use dplyr * Rebase * Use dplyr to subsample * Using dplyr in epidist vignette * Need to round here! * Use dplyr in epidist_diagnostics * Use dplyr in filter_obs_by_ptime * Altering creation of latent_individual to correctly work with dplyr (and uncovering bugs here that existed before) * Working through getting tests to pass * Two fixes to epidist vignette * Update FAQ vignette away from dt * Update approximate inference vignette * Remove call to arrange which creates bug (add issue for this) * Remove data.table from Imports * Rebase * Rebase * Lint * Fix to logo and improve a little * data.frame not data.table * Hexsticker fixes * Don't library data.table * Removing final uses of data.table using find in files... * Import runif * Remove another mention of data.table * Use dplyr in cmdstan check * Need to use , with data.frame * Fix to index being a factor issue * Perhaps this is needed? * Typo * Regenerate globals and namespace * Remove excess imports in line with R packages (2e) recommendations * Need higher version of R for data importing * Missed some stats:: qualifiers * Remove old test code * Somethings break when I don't import all of brms (because I'm not importing Stan). I think the saved data * Run document * Skip tests with fit on CRAN * Revert import all of brms * Resize logo --- .github/workflows/check-cmdstan.yaml | 2 +- DESCRIPTION | 4 +- NAMESPACE | 25 +----- R/defaults.R | 4 - R/diagnostics.R | 29 +++---- R/epidist-package.R | 4 +- R/globals.R | 34 ++------ R/latent_gamma.R | 16 ++-- R/latent_individual.R | 56 ++++++------- R/latent_lognormal.R | 14 ++-- R/observe.R | 77 ++++++++---------- R/postprocess.R | 23 +++--- R/prior.R | 6 +- R/simulate.R | 42 ++++------ R/utils.R | 12 +-- inst/make_hexsticker.R | 70 ++++++++-------- inst/vector-real.R | 10 --- man/as_latent_individual.Rd | 2 +- man/epidist_diagnostics.Rd | 2 +- ...pidist_family.epidist_latent_individual.Rd | 2 +- ...dist_validate.epidist_latent_individual.Rd | 2 +- man/figures/logo.png | Bin 35128 -> 22311 bytes man/is_latent_individual.Rd | 2 +- man/simulate_exponential_cases.Rd | 2 +- man/simulate_gillespie.Rd | 2 +- man/simulate_secondary.Rd | 2 +- man/simulate_uniform_cases.Rd | 2 +- tests/testthat/setup.R | 19 ++--- tests/testthat/test-latent_gamma.R | 6 ++ tests/testthat/test-latent_individual.R | 1 - tests/testthat/test-latent_lognormal.R | 6 ++ tests/testthat/test-postprocess.R | 30 +++---- vignettes/approx-inference.Rmd | 8 +- vignettes/ebola.Rmd | 13 ++- vignettes/epidist.Rmd | 29 +++---- vignettes/faq.Rmd | 8 +- 36 files changed, 240 insertions(+), 326 deletions(-) delete mode 100644 inst/vector-real.R diff --git a/.github/workflows/check-cmdstan.yaml b/.github/workflows/check-cmdstan.yaml index 57a61532c..cd6f35e43 100644 --- a/.github/workflows/check-cmdstan.yaml +++ b/.github/workflows/check-cmdstan.yaml @@ -59,7 +59,7 @@ jobs: - name: Compile model and check syntax run: | - dummy_obs <- data.table::data.table(case = 1L, ptime = 1, stime = 2, + dummy_obs <- dplyr::tibble(case = 1L, ptime = 1, stime = 2, delay_daily = 1, delay_lwr = 1, delay_upr = 2, ptime_lwr = 1, ptime_upr = 2, stime_lwr = 1, stime_upr = 2, obs_at = 100, censored = "interval", censored_obs_time = 10, ptime_daily = 1, diff --git a/DESCRIPTION b/DESCRIPTION index ddcf85fd5..5dd3aebcf 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -23,11 +23,10 @@ URL: https://epidist.epinowcast.org/, https://github.com/epinowcast/epidist/ BugReports: https://github.com/epinowcast/epidist/issues/ Depends: - R (>= 2.10) + R (>= 3.5.0) Imports: brms, cmdstanr, - data.table, ggplot2, purrr, stats, @@ -54,7 +53,6 @@ Suggests: patchwork Remotes: stan-dev/cmdstanr, - Rdatatable/data.table, paul-buerkner/brms Config/Needs/website: r-lib/pkgdown, diff --git a/NAMESPACE b/NAMESPACE index 008ba33c8..60cd41c92 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -37,29 +37,10 @@ export(simulate_exponential_cases) export(simulate_gillespie) export(simulate_secondary) export(simulate_uniform_cases) -import(brms) import(cmdstanr) -import(data.table) import(ggplot2) -importFrom(brms,brmsterms) -importFrom(checkmate,assert_data_frame) -importFrom(checkmate,assert_int) -importFrom(checkmate,assert_names) -importFrom(checkmate,assert_numeric) -importFrom(cli,cli_abort) -importFrom(dplyr,all_of) +importFrom(brms,bf) +importFrom(brms,prior) importFrom(dplyr,filter) -importFrom(dplyr,full_join) +importFrom(dplyr,mutate) importFrom(dplyr,select) -importFrom(purrr,map_vec) -importFrom(rstan,lookup) -importFrom(stats,dgamma) -importFrom(stats,dlnorm) -importFrom(stats,pgamma) -importFrom(stats,plnorm) -importFrom(stats,rexp) -importFrom(stats,rgamma) -importFrom(stats,rlnorm) -importFrom(stats,runif) -importFrom(stats,update) -importFrom(utils,capture.output) diff --git a/R/defaults.R b/R/defaults.R index 5ceca0384..5037777b5 100644 --- a/R/defaults.R +++ b/R/defaults.R @@ -3,7 +3,6 @@ #' @inheritParams epidist_validate #' @param ... Additional arguments passed to method. #' @family defaults -#' @importFrom cli cli_abort #' @export epidist_validate.default <- function(data, ...) { cli::cli_abort( @@ -17,7 +16,6 @@ epidist_validate.default <- function(data, ...) { #' @inheritParams epidist_formula #' @param ... Additional arguments passed to method. #' @family defaults -#' @importFrom cli cli_abort #' @export epidist_formula.default <- function(data, ...) { cli::cli_abort( @@ -31,7 +29,6 @@ epidist_formula.default <- function(data, ...) { #' @inheritParams epidist_family #' @param ... Additional arguments passed to method. #' @family defaults -#' @importFrom cli cli_abort #' @export epidist_family.default <- function(data, ...) { cli::cli_abort( @@ -45,7 +42,6 @@ epidist_family.default <- function(data, ...) { #' @inheritParams epidist_stancode #' @param ... Additional arguments passed to method. #' @family defaults -#' @importFrom cli cli_abort #' @export epidist_stancode.default <- function(data, ...) { cli::cli_abort( diff --git a/R/diagnostics.R b/R/diagnostics.R index 8787e80fd..987a4a976 100644 --- a/R/diagnostics.R +++ b/R/diagnostics.R @@ -2,7 +2,7 @@ #' #' This function computes diagnostics to assess the quality of a fitted model. #' When the fitting algorithm used is `"sampling"` (HMC) then the output of -#' `epidist_diagnostics` is a `data.table` containing: +#' `epidist_diagnostics` is a `data.frame` containing: #' * `time`: the total time taken to fit all chains #' * `samples`: the total number of samples across all chains #' * `max_rhat`: the highest value of the Gelman-Rubin statistic @@ -35,19 +35,20 @@ epidist_diagnostics <- function(fit) { } if (fit$algorithm == "sampling") { np <- brms::nuts_params(fit) - divergent_indices <- np$Parameter == "divergent__" - treedepth_indices <- np$Parameter == "treedepth__" - diagnostics <- data.table( - "time" = sum(rstan::get_elapsed_time(fit$fit)), - "samples" = nrow(np) / length(unique(np$Parameter)), - "max_rhat" = round(max(brms::rhat(fit), na.rm = TRUE), 3), - "divergent_transitions" = sum(np[divergent_indices, ]$Value), - "per_divergent_transitions" = mean(np[divergent_indices, ]$Value), - "max_treedepth" = max(np[treedepth_indices, ]$Value) - ) - diagnostics[, no_at_max_treedepth := - sum(np[treedepth_indices, ]$Value == max_treedepth)] - diagnostics[, per_at_max_treedepth := no_at_max_treedepth / samples] + divergent_ind <- np$Parameter == "divergent__" + treedepth_ind <- np$Parameter == "treedepth__" + diagnostics <- dplyr::tibble( + time = sum(rstan::get_elapsed_time(fit$fit)), + samples = nrow(np) / length(unique(np$Parameter)), + max_rhat = round(max(brms::rhat(fit), na.rm = TRUE), 3), + divergent_transitions = sum(np[divergent_ind, ]$Value), + per_divergent_transitions = mean(np[divergent_ind, ]$Value), + max_treedepth = max(np[treedepth_ind, ]$Value) + ) |> + mutate( + no_at_max_treedepth = sum(np[treedepth_ind, ]$Value == max_treedepth), + per_at_max_treedepth = no_at_max_treedepth / samples + ) } else { cli::cli_abort(c( "!" = paste0("Unrecognised algorithm: ", fit$algorithm) diff --git a/R/epidist-package.R b/R/epidist-package.R index 6e903f271..4b493960e 100644 --- a/R/epidist-package.R +++ b/R/epidist-package.R @@ -3,9 +3,9 @@ #' @import ggplot2 #' @import cmdstanr -#' @import brms ## usethis namespace: start -#' @import data.table +#' @importFrom dplyr filter select +#' @importFrom brms bf prior ## usethis namespace: end NULL diff --git a/R/globals.R b/R/globals.R index d7a6f93b9..6dfc6b4a0 100644 --- a/R/globals.R +++ b/R/globals.R @@ -1,57 +1,37 @@ # Generated by roxyglobals: do not edit by hand utils::globalVariables(c( - "no_at_max_treedepth", # "max_treedepth", # - "per_at_max_treedepth", # + "no_at_max_treedepth", # "samples", # - "id", # - "obs_t", # "obs_at", # "ptime_lwr", # - "pwindow", # "stime_lwr", # "ptime_upr", # "stime_upr", # - "woverlap", # - "swindow", # - "delay", # "row_id", # "woverlap", # - "row_id", # - "ptime_daily", # "ptime", # - "ptime_lwr", # - "ptime_upr", # - "stime_daily", # + "ptime_daily", # "stime", # - "stime_lwr", # - "stime_upr", # + "stime_daily", # "delay_daily", # - "delay_lwr", # - "delay_upr", # - "obs_at", # - "obs_at", # "ptime", # - "censored_obs_time", # + "obs_at", # "ptime_lwr", # - "censored", # "stime_upr", # - "censored", # "ptime_upr", # "stime_upr", # + ":=", # "ptime", # - "censored_obs_time", # "ptime_lwr", # "mu", # "sigma", # - "sd", # "mu", # - "sd", # "shape", # - "delay", # - "stime", # + "rlnorm", # "ptime", # + "delay", # "prior_old", # <.replace_prior> "prior_new", # <.replace_prior> "source_new", # <.replace_prior> diff --git a/R/latent_gamma.R b/R/latent_gamma.R index 95ddced02..ec3d2aab0 100644 --- a/R/latent_gamma.R +++ b/R/latent_gamma.R @@ -6,7 +6,6 @@ #' @param prep The result of a call to [brms::posterior_predict()] #' @param ... Additional arguments #' @autoglobal -#' @importFrom stats rgamma #' @keywords internal posterior_predict_latent_gamma <- function(i, prep, ...) { # nolint: object_length_linter mu <- brms::get_dpar(prep, "mu", i = i) @@ -20,8 +19,8 @@ posterior_predict_latent_gamma <- function(i, prep, ...) { # nolint: object_leng d_censored <- obs_t + 1 # while loop to impose the truncation while (d_censored > obs_t) { - p_latent <- runif(1, 0, 1) * pwindow - d_latent <- rgamma(1, shape = shape[s], scale = mu[s] / shape[s]) + p_latent <- stats::runif(1, 0, 1) * pwindow + d_latent <- stats::rgamma(1, shape = shape[s], scale = mu[s] / shape[s]) s_latent <- p_latent + d_latent p_censored <- .floor_mult(p_latent, pwindow) s_censored <- .floor_mult(s_latent, swindow) @@ -53,7 +52,6 @@ posterior_epred_latent_gamma <- function(prep) { # nolint: object_length_linter #' @param i The index of the observation to calculate the log likelihood of #' @param prep The result of a call to [brms::prepare_predictions()] #' @autoglobal -#' @importFrom stats dgamma pgamma #' @keywords internal log_lik_latent_gamma <- function(i, prep) { mu <- brms::get_dpar(prep, "mu", i = i) @@ -63,8 +61,8 @@ log_lik_latent_gamma <- function(i, prep) { pwindow <- prep$data$vreal2[i] swindow <- prep$data$vreal3[i] - swindow_raw <- runif(prep$ndraws) - pwindow_raw <- runif(prep$ndraws) + swindow_raw <- stats::runif(prep$ndraws) + pwindow_raw <- stats::runif(prep$ndraws) swindow <- swindow_raw * swindow @@ -77,7 +75,9 @@ log_lik_latent_gamma <- function(i, prep) { d <- y - pwindow + swindow obs_time <- obs_t - pwindow - lpdf <- dgamma(d, shape = shape, scale = mu / shape, log = TRUE) - lcdf <- pgamma(obs_time, shape = shape, scale = mu / shape, log.p = TRUE) + lpdf <- stats::dgamma(d, shape = shape, scale = mu / shape, log = TRUE) + lcdf <- stats::pgamma( + obs_time, shape = shape, scale = mu / shape, log.p = TRUE + ) return(lpdf - lcdf) } diff --git a/R/latent_individual.R b/R/latent_individual.R index 2947810dd..c69ae6609 100644 --- a/R/latent_individual.R +++ b/R/latent_individual.R @@ -35,31 +35,30 @@ assert_latent_individual_input <- function(data) { #' object, which may be passed to [epidist()] to perform inference for the #' model. #' -#' @param data A `data.frame` or `data.table` containing line list data +#' @param data A `data.frame` containing line list data #' @rdname as_latent_individual #' @method as_latent_individual data.frame #' @family latent_individual -#' @importFrom checkmate assert_data_frame assert_names assert_int -#' assert_numeric #' @autoglobal #' @export as_latent_individual.data.frame <- function(data) { assert_latent_individual_input(data) - class(data) <- c(class(data), "epidist_latent_individual") - data <- data.table::as.data.table(data) - data[, id := seq_len(.N)] - data[, obs_t := obs_at - ptime_lwr] - data[, pwindow := ifelse( - stime_lwr < ptime_upr, ## if overlap - stime_upr - ptime_lwr, - ptime_upr - ptime_lwr - )] - data[, woverlap := as.numeric(stime_lwr < ptime_upr)] - data[, swindow := stime_upr - stime_lwr] - data[, delay := stime_lwr - ptime_lwr] - data[, row_id := seq_len(.N)] + class(data) <- c("epidist_latent_individual", class(data)) + data <- data |> + mutate( + obs_t = obs_at - ptime_lwr, + pwindow = ifelse( + stime_lwr < ptime_upr, + stime_upr - ptime_lwr, + ptime_upr - ptime_lwr + ), + woverlap = as.numeric(stime_lwr < ptime_upr), + swindow = stime_upr - stime_lwr, + delay = stime_lwr - ptime_lwr, + row_id = dplyr::row_number() + ) if (nrow(data) > 1) { - data <- data[, id := as.factor(id)] + data <- mutate(data, row_id = factor(row_id)) } epidist_validate(data) return(data) @@ -72,9 +71,7 @@ as_latent_individual.data.frame <- function(data) { #' `is_latent_individual()` is true, it also checks that `data` is a #' `data.frame` with the correct columns. #' -#' @param data A `data.frame` or `data.table` containing line list data -#' @importFrom checkmate assert_data_frame assert_names assert_int -#' assert_numeric +#' @param data A `data.frame` containing line list data #' @method epidist_validate epidist_latent_individual #' @family latent_individual #' @export @@ -85,23 +82,22 @@ epidist_validate.epidist_latent_individual <- function(data) { names(data), must.include = c("case", "ptime_lwr", "ptime_upr", "stime_lwr", "stime_upr", "obs_at", - "id", "obs_t", "pwindow", "woverlap", + "obs_t", "pwindow", "woverlap", "swindow", "delay", "row_id") ) if (nrow(data) > 1) { - checkmate::assert_factor(data$id) + checkmate::assert_factor(data$row_id) } checkmate::assert_numeric(data$obs_t, lower = 0) checkmate::assert_numeric(data$pwindow, lower = 0) checkmate::assert_numeric(data$woverlap, lower = 0) checkmate::assert_numeric(data$swindow, lower = 0) checkmate::assert_numeric(data$delay, lower = 0) - checkmate::assert_integer(data$row_id, lower = 0) } #' Check if data has the `epidist_latent_individual` class #' -#' @param data A `data.frame` or `data.table` containing line list data +#' @param data A `data.frame` containing line list data #' @family latent_individual #' @export is_latent_individual <- function(data) { @@ -110,11 +106,10 @@ is_latent_individual <- function(data) { #' Check if data has the `epidist_latent_individual` class #' -#' @param data A `data.frame` or `data.table` containing line list data +#' @param data A `data.frame` containing line list data #' @param family Output of a call to `brms::brmsfamily()` #' @param ... ... #' -#' @importFrom rstan lookup #' @method epidist_family epidist_latent_individual #' @family latent_individual #' @export @@ -154,8 +149,6 @@ epidist_family.epidist_latent_individual <- function(data, #' @param ... ... #' @method epidist_formula epidist_latent_individual #' @family latent_individual -#' @importFrom brms brmsterms -#' @importFrom stats update #' @export epidist_formula.epidist_latent_individual <- function(data, family, formula, ...) { @@ -176,7 +169,6 @@ epidist_formula.epidist_latent_individual <- function(data, family, formula, #' @method epidist_stancode epidist_latent_individual #' @family latent_individual #' @autoglobal -#' @importFrom purrr map_vec #' @export epidist_stancode.epidist_latent_individual <- function(data, family = @@ -222,19 +214,19 @@ epidist_stancode.epidist_latent_individual <- function(data, stanvars_data <- brms::stanvar( block = "data", scode = "int wN;", - x = nrow(data[woverlap > 0]), + x = nrow(filter(data, woverlap > 0)), name = "wN" ) + brms::stanvar( block = "data", scode = "array[N - wN] int noverlap;", - x = data[woverlap == 0][, row_id], + x = filter(data, woverlap == 0)$row_id, name = "noverlap" ) + brms::stanvar( block = "data", scode = "array[wN] int woverlap;", - x = data[woverlap > 0][, row_id], + x = filter(data, woverlap > 0)$row_id, name = "woverlap" ) diff --git a/R/latent_lognormal.R b/R/latent_lognormal.R index a8c14a409..b098656cc 100644 --- a/R/latent_lognormal.R +++ b/R/latent_lognormal.R @@ -7,7 +7,6 @@ #' @param prep The result of a call to [brms::posterior_predict()] #' @param ... Additional arguments #' @autoglobal -#' @importFrom stats rlnorm #' @keywords internal posterior_predict_latent_lognormal <- function(i, prep, ...) { # nolint: object_length_linter mu <- brms::get_dpar(prep, "mu", i = i) @@ -21,8 +20,8 @@ posterior_predict_latent_lognormal <- function(i, prep, ...) { # nolint: object_ d_censored <- obs_t + 1 # while loop to impose the truncation while (d_censored > obs_t) { - p_latent <- runif(1, 0, 1) * pwindow - d_latent <- rlnorm(1, meanlog = mu[s], sdlog = sigma[s]) + p_latent <- stats::runif(1, 0, 1) * pwindow + d_latent <- stats::rlnorm(1, meanlog = mu[s], sdlog = sigma[s]) s_latent <- p_latent + d_latent p_censored <- .floor_mult(p_latent, pwindow) s_censored <- .floor_mult(s_latent, swindow) @@ -56,7 +55,6 @@ posterior_epred_latent_lognormal <- function(prep) { # nolint: object_length_lin #' @param i The index of the observation to calculate the log likelihood of #' @param prep The result of a call to [brms::prepare_predictions()] #' @autoglobal -#' @importFrom stats dlnorm plnorm #' @keywords internal log_lik_latent_lognormal <- function(i, prep) { mu <- brms::get_dpar(prep, "mu", i = i) @@ -69,8 +67,8 @@ log_lik_latent_lognormal <- function(i, prep) { # Generates values of the swindow_raw and pwindow_raw, but really these should # be extracted from prep or the fitted raws somehow. See: # https://github.com/epinowcast/epidist/issues/267 - swindow_raw <- runif(prep$ndraws) - pwindow_raw <- runif(prep$ndraws) + swindow_raw <- stats::runif(prep$ndraws) + pwindow_raw <- stats::runif(prep$ndraws) swindow <- swindow_raw * swindow @@ -83,7 +81,7 @@ log_lik_latent_lognormal <- function(i, prep) { d <- y - pwindow + swindow obs_time <- obs_t - pwindow - lpdf <- dlnorm(d, meanlog = mu, sdlog = sigma, log = TRUE) - lcdf <- plnorm(obs_time, meanlog = mu, sdlog = sigma, log.p = TRUE) + lpdf <- stats::dlnorm(d, meanlog = mu, sdlog = sigma, log = TRUE) + lcdf <- stats::plnorm(obs_time, meanlog = mu, sdlog = sigma, log.p = TRUE) return(lpdf - lcdf) } diff --git a/R/observe.R b/R/observe.R index 45fc049ff..7f88e7131 100644 --- a/R/observe.R +++ b/R/observe.R @@ -19,25 +19,19 @@ #' @autoglobal #' @export observe_process <- function(linelist) { - clinelist <- data.table::copy(linelist) - clinelist[, ptime_daily := floor(ptime)] - clinelist[, ptime_lwr := ptime_daily] - clinelist[, ptime_upr := ptime_daily + 1] - # How the second event would be recorded in the data - clinelist[, stime_daily := floor(stime)] - clinelist[, stime_lwr := stime_daily] - clinelist[, stime_upr := stime_daily + 1] - # How would we observe the delay distribution - # previously delay_daily would be the floor(delay) - clinelist[, delay_daily := stime_daily - ptime_daily] - clinelist[, delay_lwr := purrr::map_dbl(delay_daily, ~ max(0, . - 1))] - clinelist[, delay_upr := delay_daily + 1] - # We assume observation time is the ceiling of the maximum delay - clinelist[, obs_at := stime |> - max() |> - ceiling()] - - return(clinelist) + linelist |> + mutate( + ptime_daily = floor(ptime), + ptime_lwr = ptime_daily, + ptime_upr = ptime_daily + 1, + stime_daily = floor(stime), + stime_lwr = stime_daily, + stime_upr = stime_daily + 1, + delay_daily = stime_daily - ptime_daily, + delay_lwr = purrr::map_dbl(delay_daily, ~ max(0, . - 1)), + delay_upr = delay_daily + 1, + obs_at = ceiling(max(stime)) + ) } #' Filter observations based on a observation time of secondary events @@ -48,14 +42,14 @@ observe_process <- function(linelist) { #' @autoglobal #' @export filter_obs_by_obs_time <- function(linelist, obs_time) { - truncated_linelist <- data.table::copy(linelist) - truncated_linelist[, obs_at := obs_time] - truncated_linelist[, obs_time := obs_time - ptime] - truncated_linelist[, censored_obs_time := obs_at - ptime_lwr] - truncated_linelist[, censored := "interval"] - truncated_linelist <- truncated_linelist[stime_upr <= obs_at] - - return(truncated_linelist) + linelist |> + mutate( + obs_at = obs_time, + obs_time = obs_time - ptime, + censored_obs_time = obs_at - ptime_lwr, + censored = "interval" + ) |> + filter(stime_upr <= obs_at) } #' Filter observations based on the observation time of primary events @@ -69,27 +63,26 @@ filter_obs_by_obs_time <- function(linelist, obs_time) { filter_obs_by_ptime <- function(linelist, obs_time, obs_at = c("obs_secondary", "max_secondary")) { obs_at <- match.arg(obs_at) - pfilt_t <- obs_time - truncated_linelist <- data.table::copy(linelist) - truncated_linelist[, censored := "interval"] - truncated_linelist <- truncated_linelist[ptime_upr <= pfilt_t] - + truncated_linelist <- linelist |> + mutate(censored = "interval") |> + filter(ptime_upr <= pfilt_t) if (obs_at == "obs_secondary") { # Update observation time to be the same as the maximum secondary time - truncated_linelist[, obs_at := stime_upr] + truncated_linelist <- mutate(truncated_linelist, obs_at = stime_upr) } else if (obs_at == "max_secondary") { - truncated_linelist[, obs_at := stime_upr |> max() |> ceiling()] + truncated_linelist <- truncated_linelist |> + mutate(obs_at := stime_upr |> max() |> ceiling()) } - - # make observation time as specified - truncated_linelist[, obs_time := obs_at - ptime] - # Assuming truncation at the beginning of the censoring window - truncated_linelist[, censored_obs_time := obs_at - ptime_lwr] - - # set observation time to artifial observation time + # Make observation time as specified + truncated_linelist <- truncated_linelist |> + mutate( + obs_time = obs_at - ptime, + censored_obs_time = obs_at - ptime_lwr + ) + # Set observation time to artificial observation time if needed if (obs_at == "obs_secondary") { - truncated_linelist[, obs_at := pfilt_t] + truncated_linelist <- mutate(truncated_linelist, obs_at = pfilt_t) } return(truncated_linelist) } diff --git a/R/postprocess.R b/R/postprocess.R index 21eb2beaf..e8c506b5d 100644 --- a/R/postprocess.R +++ b/R/postprocess.R @@ -23,11 +23,10 @@ predict_delay_parameters <- function(fit, newdata = NULL, ...) { df[[dpar]] <- as.vector(lp_dpar) } class(df) <- c( - class(df), paste0(sub(".*_", "", fit$family$name), "_samples") + paste0(sub(".*_", "", fit$family$name), "_samples"), class(df) ) - dt <- as.data.table(df) - dt <- add_mean_sd(dt) - return(dt) + df <- add_mean_sd(df) + return(df) } #' @rdname predict_delay_parameters @@ -73,10 +72,10 @@ add_mean_sd.default <- function(data, ...) { #' @autoglobal #' @export add_mean_sd.lognormal_samples <- function(data, ...) { - nat_dt <- data.table::copy(data) - nat_dt <- nat_dt[, mean := exp(mu + sigma ^ 2 / 2)] - nat_dt <- nat_dt[, sd := mean * sqrt(exp(sigma ^ 2) - 1)] - return(nat_dt[]) + mutate(data, + mean = exp(mu + sigma ^ 2 / 2), + sd = mean * sqrt(exp(sigma ^ 2) - 1) + ) } #' Add natural scale mean and standard deviation parameters for a latent gamma @@ -91,8 +90,8 @@ add_mean_sd.lognormal_samples <- function(data, ...) { #' @autoglobal #' @export add_mean_sd.gamma_samples <- function(data, ...) { - nat_dt <- data.table::copy(data) - nat_dt <- nat_dt[, mean := mu] - nat_dt <- nat_dt[, sd := mu / sqrt(shape)] - return(nat_dt[]) + mutate(data, + mean = mu, + sd = mu / sqrt(shape) + ) } diff --git a/R/prior.R b/R/prior.R index 5c5d58f8c..85ee15086 100644 --- a/R/prior.R +++ b/R/prior.R @@ -89,10 +89,10 @@ epidist_family_prior.default <- function(family, formula, ...) { #' @family prior #' @export epidist_family_prior.lognormal <- function(family, formula, ...) { - prior <- brms::prior("normal(1, 1)", class = "Intercept") + prior <- prior("normal(1, 1)", class = "Intercept") if ("sigma" %in% names(formula$pforms)) { # Case with a model on sigma - sigma_prior <- brms::prior( + sigma_prior <- prior( "normal(-0.7, 0.4)", class = "Intercept", dpar = "sigma" ) } else if ("sigma" %in% names(formula$pfix)) { @@ -100,7 +100,7 @@ epidist_family_prior.lognormal <- function(family, formula, ...) { sigma_prior <- NULL } else { # Case with no model on sigma - sigma_prior <- brms::prior( + sigma_prior <- prior( "lognormal(-0.7, 0.4)", class = "sigma", lb = 0, ub = "NA" ) } diff --git a/R/simulate.R b/R/simulate.R index a9c3dc25a..e78c0be53 100644 --- a/R/simulate.R +++ b/R/simulate.R @@ -7,15 +7,14 @@ #' @param t Upper bound of the uniform distribution to generate primary event #' times. #' -#' @return A `data.table` with two columns: `case` (case number) and `ptime` +#' @return A `data.frame` with two columns: `case` (case number) and `ptime` #' (primary event time). #' #' @family simulate -#' @importFrom stats runif #' @export simulate_uniform_cases <- function(sample_size = 1000, t = 60) { - data.table::data.table( - case = 1:sample_size, ptime = runif(sample_size, 0, t) + data.frame( + case = 1:sample_size, ptime = stats::runif(sample_size, 0, t) ) } @@ -31,11 +30,10 @@ simulate_uniform_cases <- function(sample_size = 1000, t = 60) { #' @param seed The random seed to be used in the simulation process. #' @param t Upper bound of the survival time. Defaults to 30. #' -#' @return A `data.table` with two columns: `case` (case number) and `ptime` +#' @return A `data.frame` with two columns: `case` (case number) and `ptime` #' (primary event time). #' #' @family simulate -#' @importFrom stats runif #' @export simulate_exponential_cases <- function(r = 0.2, sample_size = 10000, @@ -44,7 +42,7 @@ simulate_exponential_cases <- function(r = 0.2, if (!missing(seed)) { set.seed(seed) } - quant <- runif(sample_size, 0, 1) + quant <- stats::runif(sample_size, 0, 1) if (r == 0) { ptime <- quant * t @@ -52,10 +50,7 @@ simulate_exponential_cases <- function(r = 0.2, ptime <- log(1 + quant * (exp(r * t) - 1)) / r } - cases <- data.table::data.table( - case = seq_along(ptime), - ptime = ptime - ) + cases <- data.frame(case = seq_along(ptime), ptime = ptime) return(cases) } @@ -72,11 +67,10 @@ simulate_exponential_cases <- function(r = 0.2, #' @param N The total population size. Defaults to 10000. #' @param seed The random seed to be used in the simulation process. #' -#' @return A `data.table` with two columns: `case` (case number) and `ptime` +#' @return A `data.frame` with two columns: `case` (case number) and `ptime` #' (primary event time). #' #' @family simulate -#' @importFrom stats rexp #' @export simulate_gillespie <- function(r = 0.2, gamma = 1 / 7, @@ -97,7 +91,7 @@ simulate_gillespie <- function(r = 0.2, srates <- sum(rates) if (srates > 0) { - deltat <- rexp(1, rate = srates) + deltat <- stats::rexp(1, rate = srates) t <- t + deltat wevent <- sample(seq_along(rates), size = 1, prob = rates) @@ -113,11 +107,7 @@ simulate_gillespie <- function(r = 0.2, } } - cases <- data.table::data.table( - case = seq_along(ptime), - ptime = ptime - ) - + cases <- data.frame(case = seq_along(ptime), ptime = ptime) return(cases) } @@ -131,17 +121,17 @@ simulate_gillespie <- function(r = 0.2, #' @param dist The delay distribution to be used. Defaults to [rlnorm()]. #' @param ... Arguments to be passed to the delay distribution function. #' -#' @return A `data.table` that augments `linelist` with two new columns: `delay` +#' @return A `data.frame` that augments `linelist` with two new columns: `delay` #' (secondary event latency) and `stime` (the time of the secondary event). #' #' @family simulate #' @autoglobal +#' @importFrom dplyr mutate #' @export simulate_secondary <- function(linelist, dist = rlnorm, ...) { - obs <- data.table::copy(linelist) - - obs[, delay := dist(.N, ...)] - obs[, stime := ptime + delay] - - return(obs) + linelist |> + mutate( + delay = dist(dplyr::n(), ...), + stime = ptime + delay + ) } diff --git a/R/utils.R b/R/utils.R index 2c373e3cd..15f92b6ea 100644 --- a/R/utils.R +++ b/R/utils.R @@ -36,7 +36,6 @@ #' #' @param x A number to be rounded down #' @param f A positive number specifying the multiple to be rounded down to -#' @importFrom checkmate assert_numeric #' @keywords internal .floor_mult <- function(x, f = 1) { checkmate::assert_numeric(f, lower = 0) @@ -53,9 +52,6 @@ #' #' @param old_prior One or more prior distributions in the class `brmsprior` #' @param new_prior One or more prior distributions in the class `brmsprior` -#' @importFrom cli cli_abort -#' @importFrom utils capture.output -#' @importFrom dplyr full_join filter select all_of #' @autoglobal #' @keywords internal .replace_prior <- function(old_prior, new_prior) { @@ -70,8 +66,8 @@ if (any(is.na(prior$prior_old))) { missing_prior <- utils::capture.output(print( prior |> - dplyr::filter(is.na(prior_old)) |> - dplyr::select( + filter(is.na(prior_old)) |> + select( prior = prior_new, dplyr::all_of(cols), source = source_new ) )) @@ -83,8 +79,8 @@ } prior <- prior |> - dplyr::filter(!is.na(prior_old), !is.na(prior_new)) |> - dplyr::select(prior = prior_new, dplyr::all_of(cols), source = source_new) + filter(!is.na(prior_old), !is.na(prior_new)) |> + select(prior = prior_new, dplyr::all_of(cols), source = source_new) return(prior) } diff --git a/inst/make_hexsticker.R b/inst/make_hexsticker.R index 24a542c36..745b28c43 100644 --- a/inst/make_hexsticker.R +++ b/inst/make_hexsticker.R @@ -1,6 +1,7 @@ library(hexSticker) library(sysfonts) library(ggplot2) +library(dplyr) # font setup font_add_google("Zilla Slab Highlight", "useme") @@ -8,66 +9,65 @@ font_add_google("Zilla Slab Highlight", "useme") # make standard plot outbreak <- simulate_gillespie(seed = 101) -secondary_dist <- data.table(mu = 1.8, sigma = 0.5) -class(secondary_dist) <- c(class(secondary_dist), "lognormal_samples") +secondary_dist <- data.frame(mu = 1.8, sigma = 0.5) +class(secondary_dist) <- c("lognormal_samples", class(secondary_dist)) secondary_dist <- add_mean_sd(secondary_dist) obs <- outbreak |> simulate_secondary( - meanlog = secondary_dist$meanlog[[1]], - sdlog = secondary_dist$sdlog[[1]] + meanlog = secondary_dist$mu[[1]], + sdlog = secondary_dist$sigma[[1]] ) |> observe_process() truncated_obs <- obs |> - filter_obs_by_obs_time(obs_time = 25) + filter_obs_by_obs_time(obs_time = 25) |> + slice_sample(n = 200, replace = FALSE) -truncated_obs <- truncated_obs[sample(seq_len(.N), 200, replace = FALSE)] +combined_obs <- bind_rows( + truncated_obs, + mutate(obs, obs_at = max(stime_daily)) +) |> + mutate(obs_at = factor(obs_at)) -combined_obs <- combine_obs(truncated_obs, obs) -meanlog <- secondary_dist$meanlog[[1]] -sdlog <- secondary_dist$sdlog[[1]] +meanlog <- secondary_dist$mu[[1]] +sdlog <- secondary_dist$sigma[[1]] -plot <- combined_obs |> +hex_plot <- combined_obs |> ggplot() + - aes(x = delay_daily) + + aes(x = delay_daily, fill = obs_at) + geom_histogram( - aes(y = after_stat(density), fill = obs_at), + aes(y = after_stat(density)), binwidth = 1, position = "dodge" ) + - lims(x = c(0, 18)) - -if (!missing(meanlog) && !missing(sdlog)) { - plot <- plot + - stat_function( - fun = dlnorm, args = c(meanlog, sdlog), n = 100, - col = "#696767b1" - ) -} - -# strip out most of the background -hex_plot <- plot + - scale_fill_brewer(palette = "Blues", direction = 1) + + lims(x = c(0, 18)) + + stat_function( + fun = dlnorm, args = c(meanlog, sdlog), n = 100, + col = "#696767b1" + ) + + scale_fill_brewer(palette = "Set2", direction = 1) + scale_y_continuous(breaks = NULL) + labs(x = "", y = "") + theme_void() + theme_transparent() + - theme(legend.position = "none", - panel.background = element_blank()) + theme( + legend.position = "none", + panel.background = element_blank() + ) -# make and save hexsticker +# Make and save hexsticker sticker( hex_plot, package = "epidist", p_size = 23, p_color = "#646770", - s_x = 1, - s_y = 0.85, - s_width = 1.3, - s_height = 0.75, + p_x = 1.3, + p_y = 1.15, + s_x = 0.85, + s_y = 1, + s_width = 1.2, + s_height = 1.2, h_fill = "#ffffff", h_color = "#646770", - filename = file.path("man", "figures", "logo.png"), - u_color = "#646770", - u_size = 3.5 + filename = file.path("man", "figures", "logo.png") ) diff --git a/inst/vector-real.R b/inst/vector-real.R deleted file mode 100644 index 7aa124792..000000000 --- a/inst/vector-real.R +++ /dev/null @@ -1,10 +0,0 @@ -source("~/Documents/cfa/delays/epidist/tests/testthat/setup.R", echo = TRUE) -prep_obs <- as_latent_individual(sim_obs) -set.seed(1) - -# Fails -fit <- epidist( - data = prep_obs, - formula = brms::bf(mu ~ 1), - seed = 1 -) diff --git a/man/as_latent_individual.Rd b/man/as_latent_individual.Rd index 9e7d35c9f..1c4d097d0 100644 --- a/man/as_latent_individual.Rd +++ b/man/as_latent_individual.Rd @@ -10,7 +10,7 @@ as_latent_individual(data) \method{as_latent_individual}{data.frame}(data) } \arguments{ -\item{data}{A \code{data.frame} or \code{data.table} containing line list data} +\item{data}{A \code{data.frame} containing line list data} } \description{ This function prepares data for use with the latent individual model. It does diff --git a/man/epidist_diagnostics.Rd b/man/epidist_diagnostics.Rd index 6860a82b4..7e9bf9838 100644 --- a/man/epidist_diagnostics.Rd +++ b/man/epidist_diagnostics.Rd @@ -12,7 +12,7 @@ epidist_diagnostics(fit) \description{ This function computes diagnostics to assess the quality of a fitted model. When the fitting algorithm used is \code{"sampling"} (HMC) then the output of -\code{epidist_diagnostics} is a \code{data.table} containing: +\code{epidist_diagnostics} is a \code{data.frame} containing: \itemize{ \item \code{time}: the total time taken to fit all chains \item \code{samples}: the total number of samples across all chains diff --git a/man/epidist_family.epidist_latent_individual.Rd b/man/epidist_family.epidist_latent_individual.Rd index 25e54db58..61b49a7a1 100644 --- a/man/epidist_family.epidist_latent_individual.Rd +++ b/man/epidist_family.epidist_latent_individual.Rd @@ -7,7 +7,7 @@ \method{epidist_family}{epidist_latent_individual}(data, family = "lognormal", ...) } \arguments{ -\item{data}{A \code{data.frame} or \code{data.table} containing line list data} +\item{data}{A \code{data.frame} containing line list data} \item{family}{Output of a call to \code{brms::brmsfamily()}} diff --git a/man/epidist_validate.epidist_latent_individual.Rd b/man/epidist_validate.epidist_latent_individual.Rd index b3c941093..e29b63a74 100644 --- a/man/epidist_validate.epidist_latent_individual.Rd +++ b/man/epidist_validate.epidist_latent_individual.Rd @@ -7,7 +7,7 @@ \method{epidist_validate}{epidist_latent_individual}(data) } \arguments{ -\item{data}{A \code{data.frame} or \code{data.table} containing line list data} +\item{data}{A \code{data.frame} containing line list data} } \description{ This function checks whether the provided \code{data} object is suitable for diff --git a/man/figures/logo.png b/man/figures/logo.png index b2b1e58d9f9f25cdf0a219e860c040c6a85f37cf..84a01d17b43266d97de6625b596618ad013c320a 100644 GIT binary patch literal 22311 zcmb?igIi^9w9Ylzwr$&(rkZR|oNU`}vhB&!$u=fTwr#ua{@wd8+@7cB)H$`!-utZY zUF*dc`CU;42_7FF3=9lOPF7L{c)S7w1ILDg0e;dF!EXf~NKIu`6u`i|sKLMjLcqXY zfTseEz`$Ht!N5+8z`*#^!N71FvfGpeffry*~{CO^>@=p;YVTjsm8+?M!pDgDhqGIxL`(I5M>bG zAz?3K&Iu_kRRX;oyTO^QWxj^r^$Rm8>Ym zSjIdvvYj3sHZQ%O5zUdv?-2v5VBsI8^CzxGvLJdxbPsUcT2!eY;XS7?ZSyKfo%aVn zues+7@*JQ9!L*tV%%D2p$I%R6QWYZj z-&XDe(u0J)?63BusZ$bMe3&%t8j^cSrNokoJd}Vm|L~OE_BVfcci0bZldFNs1X@9K$`RI3Ya^j7{R$1u;?h*{w z&3b-R9dW#U&O=B>!yXrak*v=0Tc${v^{;qxEkXk}uJ0kw!t&|C!9krvFCjKly1b&| z15y8edQYu>t66K1y4Exje>ZeRc)<6LcaGIo5 z3pFyfDpGfV>w z(nX!{+672Y7?rhOoJ2|6uK#N8u5%*0Z^7{fdib4>gV`l7ld^ny(I#zrb!b&>VNQ!? zGnKxn-z-HRScmyWNUT036&SqzBn2Y^_ld{zNag0zxW83DrSsC?MLGY+rX7Rd%W3;7ogH8 z<{CRMyN$8v_cr_%2>ER}p-?Qh_1N=>8W*6Iwc)|Fq41)%1NZdhb}du0hx;2_^>JYb zoR63kM@ibd-Yh2`pA#=b8jJZ&@wG?i291fSsw&6p1BuCMoRFHjI%MZWth!iB&6PxQ z^XcZ>-258O$3!Vl_Vabp0NNFX$vltYK%G)_wQ3lxAnaRK1aj9bT2WGAA4Do+?o>x> zqRr1Ao1+<9Wz^XK6(VB`$!3 zqxLQB69-8lHl18M!bf>aO0AUK`tCp~B{6X_dznM4?FaavROmbRl&@Bya&+{5+JE(< zS~i6aAL@FZFbSG99Y(t{-Ft326}vKNRFhn!^ZAC#vSwg-fmyS7 zYP=KT_*<_&M42Ld*Qs^K)(-_DkW?lVQty=2c-h z>$H9w?uqT{Z6!^aRm!XoZd2{aJ23lV*yXP?QgA_kSe?nTNvT4P?ANbfJN{l`sFxJf zN3MH6G&(N|TU%RuyF(~TCCK+kPv>E;Vqsw&7!iJ+d<>&Gh}>*^a>eDP`4iWo5kkEo z&{si8-|c^BCFBhkhTdwi`g1QYYhQDJcGcTvMMZrD*{qHh^bcFw+Mcd5sqW`Tbj97! zd@GV7SLU|!6542HHFAJjQ~GVtPE0V0IzLr;Bv0B1#f8;6jz_`7l#;1+e4*@qWy|T~LJa=*4?UGd zcVGffRqpH<&>U;5%IP}s@d=CPBOv**v_!_poq9QlK;gs=AU+dof{?O{%aGdR`UZ-15 z1Cs!+*B#|~py*}4bxL_3^YkL>2etSrwZ-LApd-{9m|(ZL5JW|kaC9B=kHGRauW0O> zSJiZGS9XR}7L%idUO54v6#KSW%M^$EuvzpvTr`(VWYEDebpk{pZ`#a8#Uwai7N})L z9-}F?UXA8Y($MBg$QuLFM1`ZL`-LiH(l?LLert7we8Cvnx%VKxse}wpP1SfH;hg}T zx3_A0BB-;*0Ej!Lur(ai70Tlk8Br0Zdmyt5_eAWoYq*GEi6uAHBoom9|Ad1@*DA2Y zd;p=XqB77xg^(filmGY_h2(%|`n!5%)z76aQxxw+8G%WSprE^^4Q;}1UUVqHitRA}tp&I7CpNof7y4NimTz##)Jzl1&`apGrXF1J6DuDV|ow6&F< zUV0CaKsq4pGi`+%jJpb1N_)O4SQH=-g5TO#Enw73cGwXizA+p*lOvupOw`tq<_!}4 zcd2;S+#<)y>UJ6%CZ`~MLSB!M)YJvx!(sSvv@V1L9lkc4cly) z!q+3#?!d6hf78<=HF0%=D=Tg7RdNILOoa|;W)xSd4h5Ey2Qyh&fjlbEkRCaXeVY6h92SGcJPL8>RjXU;%Z0*)s% zP)=dx8LTH!io4X=f`Za0`C@<{yP?msf05919NDkZF^d3K@(J*XjepCC;R2>@}9@pTy%1lpB?=5#adu*(L zb;S|LqtL{>gxlsA@Gt_Y*fUk_2`z~zQ9VQYQ}k0#3cg(i;U4amKy9L-Bs!hXnck|R z;+lJNo8r5=_0Q9#`kyS`4}|$$1@-Zi5QRN>(>wu&cv|z^cHBsHPrKk0TeTvu{`~@l zM1|7eS6m4=P(<#BP!DeDhS;}Ye4Zza|1H7#H>bRL>_BsMUsw-)#8*^ntt^uU?}!d~ ztoIx=YP9+WB6o9NfQpcImPdGf1Q3hCULE`{mr1o8%&&M%_X2(Q>s_if`y~lYNkE=q zJF%`hVxnd-*Hd`Z`H;G`u#@gPjj2T>9Q%SS!D(t0S$totbZbNN{p-J=ULx z&=cb9d*^qnm-M2_4u#efqhECkf07e7!>)py0$3P)7h zBTga2#-jQFu&^Xi$6D-1QN<9lO#sI!feIFjBZ?1YGU#%ueyo!M`y;0A$I$o4c7{=| zHTo-$VH^7GEz@5HZMG*f=E_#Ri&0$T(X=i68Hq$V8EYJy!*vYGM|e}hVBVN?W@%|4 zw-a`$-kC*xzRA6$;~doQAvL4W6MtTJ8Fq|B>C-eAf_2s0-I~b|drPM%B8+=SKsrUi zf2NpZx)a(LB#6zDJM{*PLSthw&CkeR!faQ8j~*6k z>(i5VsI`=C9dOJ*ga{ka=hoGEHwHqwo$A<6fS!?cK^Qc-32k7r^!{-`Gw~-fLmD}1 z*L)S|dvkKuYB6F)z?q=p;StC!6+k4lSv>$q(ip&%MiO`bG&MgGpF=4pchYxIk9m_; zXKlpP+Pkhnb!bu3>>Z zKKi7{69GJQk!{7=pUSU=d{y4ju)U5%czWhKRi3l$Aq3{{vcS2G9p`mB9MhpS5!(Wo1uyt}oD7&?V@57Wo7mv7l(%lUXE!u2)H0UlRy}&WCRD zQZf)0erSWVR5!EiOwO=tn`kOBz43osBa<)o59fdPf3&c2Zf>F8^pN=BOh{Iu);Kx` zVZjUk3B0eCk=KZ<^%NSWo_+=JdxvTO3h9 zqIqW)IOrSx{G872NzS0jXuLE{mR`y+ueIYB0O8gKj*6rDitFE8)OIOM0Knz8)D zU~|1!(fA0bOBrFAL0?PPrjW3u`~GL-eq`)e)?!!v2i{df5B!nP#MJ3MyVXuh4XX&S z2#rZm;b8dO&M3`v-5mB$?hOuxbF`uU5N|gJ`Ai{a(G31MZ{`-#iGQ1h_*k^Bfn@zQ zo(}b5i`E74Yh8%93h35X#u8ffK3`PJ6e$OR6@+V4ZXI7Yc#^r&WfCve`g=q=goL2Q zZUxHYXa6JdNgE2`ahLVE7etA-ndQlwSlC!}&<9UKe57o+YJABvD&iqJr-&r;kK3h2+a{;{k(FM4~Yavlf#m zBz%7!duzSBVQX|=$;%1q*fDpV;aDIVLpG|OOG--Ee;#nAvgi-wiw#Z+E^b$lqenVj zyh(AV*&a?T-@0oi9i8okxFvK4b-|p0_n|cS`xsX0082bTMjZ$26e?Hj%3WmeoK=Kp zcPQ3)|As1thy%WP_1KXp*I;F3h0kXX7Z5{sUSgw!tC;9@kZ$C+kmAj9gdxK%?eU8b^zvYQrSBh-d9&RO|meaesBAu;u z#sdUQbrB8Zk;apj4{@lzOnQZ|T09ohHxyg@;CW#H{Z?&XQ$Y1AbQlpU>PCtFx0u z?E2HPdco>YZtn-X$-{ijB!9bOi8 zMlD$h;xyoS+8Y3=yuLdt0Hu7;Bw*60zoi?_fu28(ZL{qwheW?8`xA0yJ%>`X#VzYO z#%&t%1qB9{aH~4YYYW*GP#A>|EAkR<-(@)4Fs}|hBcyp%Qx^ck=M~NLqVg#i?o9K> zS;miA`Sht;KWH00Tz*}3L-f5yTC2`|UAI=f7(Ah0|8)Z6H!Z2 za9!(1-aZxJmngEJlCJ}x9os+sszJlM<|}h@F@$_Ekg#!euCL5?8y)ja-W23S%bT1? zdnA~tfzKF5yFA_ZV17u4w_`=b@`x|J?L3vh@;^#EnWsN0ah z*I1U6n%bj|O!)vD1;Nkh5uB(y%}$+=kT4hX*vg905EK!KF~oEf_ur)?+b;KBW0YOO>&8T?P#1L2A8ke!)WS$?{ACVo#>_<< zTy!Lzt{w}V8lD@)(CYGGqXcK zt}k`}($aXC!|KOq_uh6A|$pAcC9z`v~81Idcu)4lm1VYL2%xeHVw$ zf%~3}iNm+!>L81%S2`Mj);}^jxik6xYv8KG!NH-0$M?*IAHvN+)`SpTfR~ajE`W&e z7p_bXGO}C@g%NBxA=9|IeSR`nIFQEf^g7sJhD7Qkkil z3v=9r4XiT4h7kM-)G#W{*aYgWcV#d~ot$Bw^F4jx9i{NOx1JdNIq%*Yq(u~qfp$_}jx zCkzmJiX>xo+x={a7PeDld>+=hYTyvIwjH;_A`K}a$3ZEy75Z-Xqn zUVCsy>-#-Kgy@N7rLb|5_$?_VSk zIs4J-C3ZNrpnnfw=PTy$Oa0v_@(?r-zK*mgXOP>S+o!-KUwY0+ZYlWFl)HOx`*#1> zEp`wi)*uLQk^+KMGt*TKO+3AR*4&q%TS$W*FenN=eg||-RSIVWm2`vFc$BQ!p2OJt z26!vC2&GaGstpIPrmDn|5tQkN2j@JF{b^EA+l$wUHGzqlHF(~lgq777X3b9_lW=5` z{poE}fQ+p%^Kh7dGrq{q=pK`!5=ua(-Nu0EDF@Y&x>F|q(}qyDOAPzd3FpP%vs9Zc)de=C ziEQE2Q{sx%ERcf!#<-M>FpyH;HT?VCLt5LwI6g!ZrB@v|CT+p)LlP211_&}UmYCxW zaAD$0p!Lhe#5m)|h~7@nLI1Wxn$S$S6fGvTwf41Lsh@=rf)-STm3Y?5$V1N_CxO^2 zeG~kEq#&CDUW%283}XxePofd!hWjtfa{9_B+tu=uqP)7iv*S0s2*fiompO6zlbQ0r zM=n&0L>QNy{A%jGKWXv$oo2?ti))};b&x1_d^sW&{QQ`ysH6s53il58w;u-V(O3lp zqQ!!v>VGfMk}@%2`F|b5XfEpN0#6&4l9qlk<@*f4uH%Hw2jK}ae{De0L+_Q|BKECw zdDG!|Ke}?c8&KbK;+n9De>|t;98YBigl27hg`b(!badCcVA;a9?7PFU$0wdh{iOR} z5L!w>ahsppQ@^!9BWPFsr1J{0q%r#Gp+9tGO%VhGlp2vo!wFnvZ85fQvX)et5c~5R zNpL*!Z;^fo5~JkuKD4;+xV1aBuRFH;0m$3^@t1_xzbhMHN0t)Z~0KhX>J1DqNdZ zZz`RblN0*Cg3I8|_*#c&G{6;5WtEhQid!PK`@?@F<+x65#pyOz@_AqM=~s1Bjsf!@ zuJGH(!;I}}X?=ZU6)-($YfAx2e4-^(^Qr=X`_RvxoZQ^Ldtj9X1&J6mQ%39cX2q45 zI9Tj{Z|TA!x`0iQD+J5ycf)yqj|Avyzb!E(9opBQ-syf=Sg<(_uwqZ<0qa~Yekxp? zY~y|W^=z#PM;8k%Lb?c0aIy(_?A)^A5q5{;OGrg76Y~5O7^EY(eI*eXde1mn4J z{uMHgjOSMrrczYHP%ST*uKPFIX33}3G`PMhPUd$0XJVzAxp$7J_wKWgPEMM<{22vk zYD8b(oA8vTuIBKR-!;L?DrCz#^T#zw*90(DEj53#lGtsG{T{w%~hucc-D-s~^8K*3|o{h)^el9b8l6{efqp zbrCIHY1TfKshsoh+8#bWzGY>Fgdj5~%Wh%~Re;x8G0NCL`R{$TDpn44THkN|hkz1T zv;OM&$a1M_xH&zWFe&o&+g?+XtIyT8#@Kl6Ra5N;5rqum>3w#NcO$ttC3;C|TPl;T z`%Sp}hq$*lDaeqynT3KJVEXiS8aM5m-STGpOkHn-Q;E5`Zi}temSy(rM3Y$mT>r-7 z-&E3DOc`vs`#*7{)#{yJXG|)&3Ew5k?sbJdq?)a}k)j=c9WYw!&ux2w@w*k--2S)E zfEzZmuW)eharNv&T0irh7Pqf~Nk#P?7!O>1h0$YgK45h{FYU8TB`!7$*u7RiMN|Bo zo#eiPH%kWjXgQ|lV>#E2+!eFd8hNAz{H(09va=R{6XtxMX-2j&(yAmC-q-8d(u)2b zG`(qeS2=dk3Hx+@xlzz5G(Eoy2xBTLda_+eO1xKYx~W3PlX>mbR2f02)Y#0SXe5f?l_Qn5r`ed64bu*@kn# zL=2+vXRFc%mf-jA3FGNA6!q#A6*VPI)#qJ@m1R@ehh-w4TwGjDM0sf6W{c}s^L7fw zfeQxICue7brl$M3g?aZctt<@G@`!dAp?eDV@U` zCW$(ZhJgVR6{zgFxjCo$H}){&?`mq(_s<6h*v$G}eOv#~)e3GFn8XVl77@zwNkzDW zj#2BaaQ6<(AICY8k|AJN2Yb%a!|_d~JEPltp4ow>^>+-Smf)oqfo490?U;LCY+Egh zl>t#I-e)0KuhT21`)xm%`TE7j01yQG3PYXS4+<;4>go!}#1S8^cl(G+NKk~40fKr8 z)2iD|3XAVWZvQv5g33x5bt)Q~XeHssU)mgHiK(e8Gq*Z=?S8=`UAM#*lew@aor9C> zuj%Y&iLX(pm6es-?%kglz@5|Y`s?Z0>f^JoVm_Jv z)ys=Nj7+~%U0xoc`{gf*@XcF)7_#}h{7+OgzsGQq+vdmd3(pR-Wj)637_()fMOAF( z^zpRrr)2{OwxI^g8A0+yxGp;}v}Ns#-}>)JU_rg8pVv7g4(DD{i_yD{c^^=X77CY6 z^k@Z}L4q*Z3}&pjd*)J7jEm}u_}0QN@cC|UG1m;i9p7=#6#v}Y&uXU*KIC|pS5%*# zcxnQg=b`W|A^@N*dNnmytuZk!sf-3aIKq!|9UDsYI!AZy%UlckX#+ zgExae%*~B5sl$x7u3PZ9oOli2$8cyF8ISN6xVgBTkEUnZ!~E)}1|9+~d6FV+S8Z17 zvjEU`xxWGK%G<5TyAYtDkO>GNq>Z1o7(zC$N&)UuMg|@UpX2eL;d3|BOcL98>$Cmo z<3Mt9L_EEfvL3xsxu~hdBF)PU7~*&Ww4(4F<}~2nLl;kbPkrk^FSIb(0R$(f6*ZOO zf_9uV=QB!6fe?7q-{c$FBzcwRB_&K%yK}lL+rnncrUq-~>rL8z%)AMcP@);X{Gp4% z75**%;~ptiugm9t<8h{5io^9r()@4s%HfsMbiRrCROW{1hUbMipRH~DMWTkfdjB@X zaJz`mS*!gw{1@P?hYxO5a72HZ?P50B;)|iZwhEjK>y?f_0a{HY3sEc zTqmNmt)83T?0tE>Q$t4^iYHJ24Ae!0LZI0JhJ{SmqB9P&3ILu>T55WFBQbcYuhmQ> zP=I^Q|KsX>kOS@EjIRe6EsN^w>9&JOf4L#A;Q;okmIfv}N6+aC+|(hgM4(MGJDo!7 zyrtvah;r$-9BG}EqQ^TfE?50q)_|-O^|z0B|2W&wdncbVg2R$=%{C1kfk#&nYv`rJ zdljweN02HFAOid_FP+`kh>Ji)470TH8bP{8W*pxEA}|0rf18z5$K&(C1*$lHEevIR zB14tWrrOr_zc(=d?Ct-S$r=YxS!EraAQ396ZNW>Jm;fnKLU?6$_0X@3 zxCO&Q=KlJy>B*&p));nbY7WY42a#`R>U%g#8w>-7;XEHMHP(*4_4 zGBwM1rGb~i5J+zy$#KN#%w-{>q*tq$3H5I2PzDWH>)A>6s@c&I9{!b@3&w0Qy zPOBYw2ox~n)U4ewI(^VbBo>S~R;A%EeZh2cda;0R6~SedutY0pDoM#av0fuaF;|qi zlMxZYiAuQC7%HPrbVx6fi9>}CB|gutg!1MKUmI3B?8knsufH($S8h9x-$0mk*=ASZ z9|i+e>oZiKM4M^EJ(*tMCjpL;{_AvCD>x6Zk(|#@)TQI` zZX;+G?EKYW(vS=DNiPG5F;_LDV{ljOYnnl}?!Wonc_nafM&$O?j_XA!#r2B9+LF@J z#w|o{84(Vqw&#wW35DE#F=*D#tgI~5>(nG0w)L@m+&2IA`!K}ic5&!XAl8DYDK94_ zL{3V}L&3vav#vOHm$}iZwcsdkNl>go-?{+u!Kk%CHlbrW#)(x!qBGfv3M8*V7^1l-nW0xX{8 zPH&M{WI}RGcvIZIYTfy>rPlD?5-P{3C_d7#Y4J0%@-}s7QOE0DTP-T=Pk!4M7aw7g z!^xao4fXl{RTgm_YbC?ogM-wX!E)I5V;1}qf;Pbnwcf+q#GI6S*+Wi+q&vH4zsi+a zGtR)P=Z+YBzU*4x`{>_#UB$yX{lMtT$VTR}VQN)J*@RgS;?mMz zDJk+29)7YbAQ;OZrv@7E=!_^1%l zA-c_WXmC(~Q9~{ua2Q#Q9!DgQ+E+*p^gUkhg9(W|i|=x1rqb~}Ak+SX1rYbf>!~^5 zsuSyG-An1}68QUFN&uc?JILF~sUd^U)VdpSW;l*WN>0vk(R~iEU-Le1gk4?lAGx2e z{|v~OSN~y*rXT_*%^9{ZGfVh4?iGMrzx8=>Skt2Ytfg|+CyK`URo|9r#xf*3Td%`l zDrq#O;H$X?b`rG$V86zZh&6fFH*B9#K_j0>gXL0#a6TZT!kjj`jmCzgfQnKvS z(5MA+5sN3(8IH>{nMLvcFU#T8#DpG+8YlUpKTtXL0U1uW+7LXhS@Umh?i(*U8GbA} zB225YW0L#Fo#w`;gDwhM93Dqh%w2nPw_a&95Pq_w5Z#`uXv`LLO%){-WyLK73?jf- z_dEXa%O6Vz06$e%hJzhV_M_2iM7K6p#Z1mfNLXZt_ZL^@lr#>jo4s)U%_X)XkB8-y zx9fO+GF*5va$uSmxF_JDr5!0F6^XyQ<3#Ky#e|=kUXlbRAVB&3@uQf-<_I+^8bRcS zA7Foece{ovYHFeYn_`}&L`78A)g1=72S7&G<4~~w(>ef=vD|RHW-^B(iGYWGxy@HI z?fX<-_sD$JA|3t2C44W(04Xj!;G_@467als-u)y+H3Cus+1=k2rKS0S{{Ez@;?SxW zB2EVMJMR{?l}dnHm;i~B&J~kP7=VzQUv>op>I$Mnd6oTExX25%C($V(aB8KBIV1l8 z(T2=7sPs1oX&^E%x3BFv>LS+Up3li{ppFkBr2%S?2-|(<>{)x1|W- ztzPFanzidkp5`Uve>vU!2%1FNU{%lIf{u2lD=YDV$OMGR=n=o>r+`jTl@R}D&7b1* zTA%0JWMul|+e~0mwm68XI9qqG`^h=~Z{@41^Tev6`gPgM7pS1P0Pr&qX)nmlNkw&Z z)YIv9QxT9=SCU7wml}g_ILcZK{A0GY>11c%Jix#qJ=X00iwD&eNC3_}3NX5sRa5|e z9t77#a%QX+F`CR5RvBF3D6SVc9%ol}T-7Rq-Ef54}BPwc7Cze9ZVY zlCVVYWg$$vOJPt9c>X>g&Uq#iN6Kd_jd?dW<3MS+QPQ`}!g4>T0A#xO$>+bC`Z3ZO zeivpOFyN&yz2R+-=H%?$R*LE+bE2Y7FzW(v|G_=EUF|LxG%BjcIH<>0)u*vEVHscm zghe4L0oXm@Jf<`a?4MlCT!wo(-XLg+>Oi8#{ZdX4({Ou`!Bp~eK@6!Bg7g6 z?Ip0ur-WvP3fZ-#mH>1h6h31P#Nb4VaJCCXZs^^0;Bnb2Y8CI3{Ig@@HmdeF-;R=Fwxup`#OK|2GSO zC8dZ++FEvYjcra>I+rs(5LA<7LP)b>^>8@W9|nfou{In{Rc$lIEjYZNYafW4taP8&Z=yrsWhr*)$0=kV6=>q~yi*~FmjF=)ju4~Mi;k*c!mDjhl zaRheF(aVjcNqKpT)p3bNxG0x!FBQ`;rv~y*`WaFtwoU>i=l!5at4u(nxR<|a>yq^YU*+Z;*6@ptiUO9s zB+2YAUlgcs68RLbL1Zu#?_;F3P9+rQlSS|?!}~AG*8*ld_Nu5Q*QvEh(Ww|NEAnjeW|c2!dCL+}$~46_vXo)*{6X6os>= z6n*G*pGF`tg2YS?xE$3iBFL*e`>9d6LXq9+6x%vj67Q!18!jySLYxjF4IwW^r)n=j zZ*L&%bj9lf&vh}ggKp3t8&IRdk99g5nHfHj2CT9Uhzl1^xDe%W>RD*wS+Jtr!!Oaz z7!oFU0u5_@O{a@3YJPMZO#7?>at7d#OYjQk`CkKJ1<^klzktLOW$>~p|4VWa@_pka zjBD0LhD@OhWxXL+AdcXT2Z^dAt&5peli4)EN(*8r6uAtgmjLxbA0l@BFHXn?F< zcE?^hK@?XUO^}&vF-hxK82BJ(X}%rmts_t7@=5If%}Pd&Q$jB33-U2|JE+w9^J~Cn ze9Y%hp*?b>*bTf^ge_TU(U$YB6V;@kD!3X@`IP<428U}I9E3t%$zmlXj&gY?POUGr zUIz(hl$WSFSXM{_gnTNUTkIGnAFJzKnB95Dvld%!?R>5`l=JhCn$07<+kSE;;E;do zqwb5JO&}C`lk0bffHD88zCO`#JW+JB7&2tc9#);L1Sxh#k)M4B$oDd`@>lAK#Y+qx zZdmWN>2P_K_uYy(GI9=OMF8{{se}jAbm}_^w9DbqkVp(VX6Dit4XR;xJ24g#neq6- z*m%0llzvqSCwQ%Z?@&R;bGAcXTauZd8)D|>)Uy^70kY{ux|?MeJhcV>PsloclpCdd zlzgLfMggKwebKZE&fE0^kqgt_9sRD0n6VaLo*DqOoku}})D{!12*cj7>#(M{?% z4U`_Yuql?WGi6Szb1BuctJ8QR^? z$x$61!iX;NLe?aabkKrk;BtIZR#~sx1tZ`9X#lT?l z+4+$7x>dLUbOC%Gg+*Ha4*8wfEeOLc9{H@Lt}4KK4LO>SidPH8}u4%fkXPj zg2sRu_4`{h5B1?Q-mWG-ld(mSSP%a9OWN>ZQI#((U0u0Gq+kZ_ePuJAxNoHmth3a! zYLiU}GCAqS;qjyuxd~w}UVR0ucK-0}VBSX1VI{PU9GD2ki?Pyn>H?9MgfG+;P-|;z z-{ztKP7-RBiTHlxH-&H zNfvI0WK+VU)-!<&L=;q=xeDQF5mqv3(OM6G4r`1!(aAU3T3e@J&?vH+?8%x>Wnf&c zI?0OulIVe$u}lWa>Dw%&fB*qo=8o;WP-OYwMWy_jWrZ{*V)9I{2F$2TJOL)|uBN$p zPn>AFt!}dl06LCxavk>LWWBPPWG<|z`6VULwTtrR=H`9zxzN3skgAA!(r2Ih!r=&& zL6Hp(mdyQOh=52K)$liqh(|uWT^1fZ8txtm35Z$kKE95Qj*>Z6f{Vnx#l+qPY&AhEN_aY-@uLczg3%vs3D* z_A)7_|4<0|;yZjCUZXHu}SeP_$ODQ-wH1_=hTzYes z`Gim00L|%X;Lh<0DXXZuFP9q9JK#@Xe}B@Igr%6JmDS+vQvC_R+ zT;6(U0XGOYrrR}<>T!vNIB(m;xKu$qX*5ontQo7cf?w!b_1g}DR>;thv=0+gadB}% zPLjkpm{t#Z;Au6|S#4KPMHH3zS9v1&i)%bpBe7H!dxr2##1FRv{){vh4xfPtkIguV zY}q%mKce|fDc@-&0c01hW!b>KCkPub2~_CsVMUCIKax`?8~iSuYRfT((&T#I@p{p& zw%&$ysp}Yl5y@htZfQ3Rvb}+azb<()X>!d|vB7FsA1k}D`P$)dzQVQP`=}4t7(2F^ z5`?l4;E)QovXqi;{pw4%uy>YtOG@_e8qi+_x$}rq@`cpWeJ4tAw$a@AxlfFO)|91t1{0?CKpr^ezyeEb4vq;au zXp243Ic^b&OAy#5W|IBNDhgLw*W9$%tdi+$Kd}=nle<5j^^wn(qH>oeIry+Ew3kOV zKkXB8iR?Luh%qA9Q|M%!!J2Fu*6~l(cWhN5OmJbG4<>gp7&N!ILLrf`{VN)*zUW5G znm}d_V6%#<3M_rCA?OIZWqk;>!k-}G$+oybdZQ+^X7ecfT48S$4WM{@;s91%fFzS8 zf(iVdck&)R<;xx$w7JES5B(x_Khp5emFK+6w~Q>s{eEzHlc~U9v|^@+)d2FTSt@22ZUe}TDxAqY5vA2yzrix$TRiDl{LpPiFl1VE z%^q@8#G>9M2#pt2PyiJKQh}TBZ(eydL>n}f%%2%gXnal6ODmfLYFTP|;oX+Ri*_4( zJl!8$44O4R4t{I5tWQ?niIO=Q4#NXxs}AHY6XEfCTLTwtJIzjpdeMa4_=zkN3`j{2^p zHIml#W|qpR9TZXba~Rq0e%#RMS%uGy7Y>AFuS!H5I6x#W!&poe5S5ST*Q2fAi;GDG zY)4TNd}V|Vg-(9IcGu(e3(2h7XJnMKcxcR{alM=K$U44CptM6Fj1rXxTciH*k~8^; z1~}0ZnZh}>v13L)Uz1Zyr9a$W@Vh^Z1kzXxhd|l6ym)foz@tg4#o(l;F)sst2V>$u zg*|OCKy#B!ssCGpDVx;aIuv=QBk{QczrVj%axKiRbVD+@0;#XO6T_FHyw(jqM+!%KKM$CuF|8S&k<5&sVHspaq7~e1O&!}o?}LC61jZ8F#GFt zDEQL@v1cAm&hOfu7=XjU(^>RS6M*uz2v5f^VJpddKU8ieI3_X#(tuq&&O29PHFTLGUh8A4bibz^mWTJzk$fO@;X28_`Uk(V5{RdpDrzb%C z>jzmZz|Rtos-7><9RQcG*9bN5a`ey;#29^Z<5ZXm-FjS3v%SG07fdQ?6D^O_93csj zA7w0shmExBFWbxw-oc!d=x5e4UOL#&THwc=s(5JaO4P)~MJVJGVG8^3D|3ei;L$rB zgkcJqVt&eKz=!A2X&^bEJ>w^unW<}u0(J3n(1!VMIUfiW0vdGVf7?W22|pymB>*c^ zugm3&W0s|#w9V(@aR3QCxDExGpOZ1kJuykpPQ9%GGAD-$)GOL#H4#}fUMV}<-~NFq zpD65WhpHRw{qRSp^DX4c)qXt&z!0ZO>MWZCLTXLg7H~d8iY}7_=yU99KSh#<~(wTiJD8pWZ>b&I6_-iEh(BxwxWe(?d+Q`ke*~^ zw)(yEr?Keo^wgYttlRzj_ise@`j;kuavxd$$ixICjX>7#Gn+L#9=m3OTRJQ^PW*k7 zXS4rOsB_1SA>!NZGR#r?5R(wQfDmB4kx4aB&aw*-p_PTz3hg4ar{Ay-~rQ3YD z)6+`* zDa=VBofN?!6_u*mwp;$LEo*8@IcT%CW~fo#nE6Y<%)?xZ5UAQ{aN{{D9@ z*-w_KyJ|`Mid0?pIBS_PCzVen6#Wn*zM`XRWp7V|9*IbM9VY(Y9yA3Y_W>9_fE(${ zcxQz{->v$T+oLSTxCrbzzkWNp2WLBXpHCGk^uc<LWEX@=0Fvz%;NRjLH|c+WWqlYxA8;12{tSy(zk{KR!s)Y zM$BHRhu@Fn)%h`cLy@6K|K7Q2Ja={2eTiD=yyUd(*5laCXV{{B8FkfphFI4iN2VCrfIa6y{!*c0HBdC zdalHhg#0_}F8B@pfmhIS`y)pjD3F8N%QVNGJv?Bm`^cI0tgQU$ItDP`{}-q{U91`O za?UTk#;!2}`x08f6@nGXN~|FfnztFSmw9B!MO9a8SzOlr|KXBmW-R5T>x#iZ{dmby~r{C@Fy1eD)-Q)Xu zD;fx&tr;ipcnH1q-CJpD(^~ruh5p#?@BHW_d;-LRU?T%Ye*`mih z6oWSkY^F>HLchdHQ!(K{J?6dd(Ey%%X?P%c(6JMc4w@YJva8K$O3x?Ep`kE$UNEzJ zJ_dOnFFlF?+h%ZJ++Vzs_HX_*vl0_()zMY3h5q~QJ;n+5z|@9m_L>$R3E6&+_ zCrKhC&e_Qw*=LW;LbhaQ&m;Uk*Z2EQ z1g8pxCz6JxMyT^9S9nJ@bM}?+4fD|@{Ne!(_v_`X$2H>Bx=Cl6)Vl8%?_%put4cNA zCU_z?jyDA*d}h|U1MkL%>Eo?cwu(|9VenGDs@t>?le^x7O{tmDC+^2jC&#XoC8C43 z_qj5$=wLQdMP-xmr){DU{2u&oI2uu0f)Wn&94wo;D99P?{s@-+fhk;PD0}>AarBXU zI^Fr2M!-Oj2G_2$v)j^qE0^mGDl0B4okiDkT%R#TEAm5#`*m~)bKjV*!bx)dfZt$j zJMzay>(KqzPKsk3(rgfVKME-mNHZXaoYO3TH$R*>qch64J&F{JT0a8NsV`fFdfJ_S zZH*B;^;1$)7i#CED4qiQDKHv20WyE!F(26NM>P+AVyv?OQrMZlr{WVi6rYjA%K{nfN0h>?cElCh+hoaB;S3+~7J zO^=2ts&&3K+#F0M?m6kSZXb8mv#>Dd@FBIVg9v5(HlfEk=ih&b9ZW9z_Nxqpx>KJ- zj_;WgSt-uW?o0+*uY}fjtvoiTLbBT8?|fz-3@P+hFeyrrbC&?FK|36L$8xWikw-=r z0>0ae)g8c*&GZ~HL4f;Ky+S(?BXs4btwmFH6wNSN{NkAS>zHz~hgro-BFp>ALxaV= z+F=7K0fM}~WK!Ts01QevnnOz~?tQ|wYlN>fZuNq;)lUHW-*Xzcu{BY#1x@A8-7$MT zu0I21PfcWM-q#<86_BpKa2ZQ;jBu`=`9v_M?{A|LHK#W=V&~q~WUmNkji>)+f6?&M z0qQ8ZZy${xiZnDZ=y~HNgL7$=mAM6xjw>xuV>aTN?B7jJ&;D6(QaHPmJE%ERRUE8% znK^MaPw#3I(H@j0p!Ce|s1tblT_>nb?^^bMDZa6I9kn_?#D?GpArWBI8p2Y4WN2+Y z<%B56&oB@t{p&v(wQXE#b56!2ao(kDRl|EeNhIIA4AF-3iVtj=Fpo|KdKi_0q)dc{s>uiz^Neuspp>mbH9dmdTWp?&03C`9g zkM@--bhD3lAcv8jpva82r^~l!pM&!1#jBpR-gS2jg=ooUxTmiT zPNeW zF2V9ZrY0~bFN*fHF670iE}K-low*;eNPdSbkKdCw!>=+qkle90hC4Cd_s_@O5Nz?* z$&K}vH34#Ac1ZJkw4f*md!Xhqr6)M*@Vk2fm%sr0MT>{h&e|-$9)jBXU;M<;%=9|R zKH(N!MK>i55^LN|3@Z_W)jyJeUH<5LHg$ukLcZBo?pYdV6kO@G4|_ay0N)p*ZZMSs zAsOwskPLC}9py*?X-QdrAdP871E1?$WtNsQVvp#Si*ob75d`C51*5zeOk75plA9j* zD98?zuB@AERLRECa6eus34<9$?&fQFxMzZiR{sZyECAT00>G|6nI(th4b*FC-KQ7z zXCBPfuMo60nP#*h|Jf;W5LE8z z<}%Th%sLd0RETlOK&?j}#vm-+GUNBIV~Pw612+=?@remfa6}cG!qw!d^BsC$&)V?+ zaEiVEI7P4OlKs^VKP=tqLmKhPZdLrvLAN|&ZjR1^!Qg#J``G&N%(J85x~*7*^|B+J4!E}B>sfQ?$;y1pjY|1gizADlm~H&uheF&CbTm;SZ!M_i;F zwjp9d4Jzi=wUzO-Sn1wJG-`*Ehfj&4KmY@v zzWmwj;JFCyZX4CdM(NEAPfNIY6sDR& zESW_MaJa5W&?GG|KB<32wT&tQEvxt4v|rf@0{!OWub&fg5jMsPx{7c*l8v%b*{f>|-_!97L+zijw0GXi%%WI0On{(*3|R%7WI3ws8hb?C(`Coa#>2 zFaa&IY=Y7j6I~b$v1xplj*n^d8r|3xZ#KJb+JZM&W;2($D>6zgUxA5lp$N5UqZ=I^ zeJl4u8W{8oYHNow*<|pl+T&J)URQ&?2<``7=BAImiYxPaq*{-S4f zRN6}9xS!pGS^3E1-&Nx%yebj05^9&th{`Rpbo?h~gH9LI7T~Qa-wt89N`Tbr_@;$& zpXz8;R(E!G-gHqISCXP=5X8u!igg`_et$#Xu+Iq;dZpBoUHYj08vdwS;SeJV;c$rB zba=-DYH_=cJ#L+BQJN3i3niM8x7Yx?8IH3YA=f9%`yZkMXUgd>AQ5O|Y*leN6UuI_ zP|}@Tj!%<&UM~k=;wnM;pRE?4wspBVv&&sJ6&XEgoqpi-xn~k)*$TH;-6-O_Gc~Yu zXZ*DEdF~KjuLtBj=XQXV>YqdpSu4{7|5v~Xd(w#>TZpZ#q>xdc@5_fVy)6hEcA0qM zcy)hyz@g_C5iQ+snJGO&ioVZ4#v*Zrsn?!Pg&n-At~?5Y%H3m$6pctxz?(boDcNm` zjOJYX@XG`ozGpt#_vJ*b7G8#WAcQpn91qvtanaM$fA;QCTJFzLR127Fv;NLu2NE)n zE>G;>XH>)Azr*}s?tdOIEflmI8)M7*+{V%30a)xHn|!t|XUYM=mcx(rLj*!QB{QZ8 zu_G0lFky*axv1(}dh#7X0}NDefuJqyO4`NuXyEjWirH~0X~x*p6#qcnSZ#>R7+$by z>2GazYceoB+`*xRq&CgOWh^`oZt_y?wGaPY+K!niB{ML!A1&;}&_!Rn?X%pasl-RS z{C&;*?Hq7s`I$(qTIF2-;W(d6G6z`_R`Onq`kl2h8;{Rx8>1={IUOd0BqDiC7^G+` z8)Tf;-P7aP3WA6p-A`_qW_E(V&jlM9vg~&-p-B+>prQ?GkjVGvip`CKpqfPA4Bp!U!u7>01mN4{a;X}=*G>JtOg zYMq;URs)`u(-Q``%fQQCebE;~w}|sJ1Iq-o50wqv2i1Mf>s#L9#<%}6r&Ib$?uXS9 zoE{Zk(s}H&$vlhPt<{@zMo@>QAbMq@D6I zwT(HdbI9O@0{^cd@`taBj-0IyF#*Qy)O|T5drEiJGnS7`+`LwYS^<^wc5Rd=_0_Kk zfg=f|7O*eXAN_F}FAoBA!y)8DN&ag>rqB;#Np23mH4=+I|9u7ox-OwwyX@vWUH7CVZD6#Yg zbDl=oZ8*2g*6v>PdR{Zm3BXrY_i!atdXxM3!ICPZPa029nWoWqHV_24oNV625=9? z`!6{$QD~dQV}V-~Afm<&oV=JP$X4Sa`cNTGha*%s^c@eE8ioDRUI&04V99xPuVNS;IdnNEgdUB+rGxOAG18_htl5b(mBu$J1Ki9TCso41cE zRW|3G-?X`?HR=Y~+3hsm%IDZ8QB-W}wa{5bVTNDrJ)^F(pThz5a&U8*zC=%nt9P5a zDcX)X@RiB6S2XG^Y^}(e%Ot3}aC5o%nDxbY%7!biM1H)N-!=2&&4+@B7kVYKsld*@-LY-W$|!_-I+zF90250-L!Fs3|VUs znAf7sKFE2X{}k-u-NnvI(1PtuVvL;1wsOyObK7K=ex6Wf>G%A+^IQE4+i=paaYM?# zO19x{o~#}4vqZ4rE<##C)>rjuq>lC5+pjN$kX0~YKF~TVwz*F3H2TkB zLcl*@@A%_g240`O*Fyp<-Oe$=&za;kR+J>}fg||1|9f$o#xC}kn>-#RWg-r+!98z` zi~{|nv&_{fXzL}@HZN7?aa-JyCubR`r@DumM zaLw~}k3WXRe#woGadrxYY0giz7;oNj72(RRt9Ggcmw8+qS;_6K6pz&ljukmFuH)uVeFgcMhxX_Jk^@B5Jzxe$xw$7quTIybR zzdy__9;vld#t*_Jz5&B4E?feXR9BP2U=Qd0$Pxz#Okkc@nAJ@IPv7!IwTMuN^OPv1 ze>T7#KnVaC`PaS$k$PuRu2sdwnL%e%ICLFgCWJx8Qu^f_fAvj*`HDNmlvWG3Wrme@ zF(C8LqULBWR*H}DHW{1pAxPuk|Fx$3??I!J^WW6JZmwtn(6H2dFWUy9l3pKQUmh>( zRZPQd6Jg6+ME)C@2zw|8q8%elT*mX>#4`BW+-=Q)-54fC9Na{`%?3*Jvf2VD5f)oV zC943&7r@4A0M;Iueb;6BsS@#~44OW!Mqj;I%lGl;XFKbLtbTm;7SgU}xl6Jt8Dw*2 zOEhsNu=#l!P}Qn};x-pet)QyPuNX|`n(zgZo;6E%zF>Lj^uN?eVVCj03RO?lgbg9` zpH`(q>)5xpd?wwa1x_N@bRq?^|x%uvX*-OhVmzW z=zk!pzrLT4B+w(su{5${U>w z@*Nr*D^8~QR#`*G4Q7>}&Kwjj#-vh!PvQKE)SNhyG_k?Mz9Yec1GV+W!pc}fm^Qz> zyh*4T3i(0aay9aUcDr??(0W0&CMQ@Fj={1EKfGjyF{Mx+q%4<)mPg4C_Qt|$0|Y*t zv|%A`a>cc+_<2eK_;^FpGHW(+&Jc>oLGp%zZXKG_1oZ@#a(u;$LX4m5)C6~aDtgYF zS>zQ&y+w0ivPNj+`r)F4@~1c$p?TzLWlmO?x2cpQqUz@_bvbLWUZ|wYHcyE<^PVjm?X97;`Nm m9#3(q>)+*l^wf6ulGvOftcK+GyEYJ0i>Ilor&6wD8}@%o$<*2a literal 35128 zcmce;byStn+by~Q=?3W%5Rj0P?%0Hs)FwnyNdW;V>F!hz>8?$KN_T^lh|=8+(s|bQ zcfNDZcmKQNo^c%@W9Z)dUGKY|m~%dJ5w5PLfQv2IYx(Q8`fHx>)7cv4U(^n@-4jxGg-S^h57 z_m@kt7RxJQN97=DC}aQZVc7QlhB!cK2A(n8ge93eCij0 z_eQ?-RZ_@%%16~+NSEI$YWD*tCnp22sHh2raY<2fZ0y&6e84RwIzalOzD^LHO5pxs z>LaQ#a@+u#0#psGrYI#v-1B=yw~Yr4GS^D>HO7GAB9&=pQ@%1;n)=88dQEtE_}EY79?!tGGEI34h8&-VeR`|Sh!=#Nt}~9sMpL4-Fi4-uO^?3 zbh&Sf^X8%+Qb8nOqLlTpanrj9IY+V*S8Av>Ee9I3{JBadg}9oM7o;x)CjKR}5{s3L zRKJT(xBOL{{zM>s&d=98kDPckG{)@s>B%7Ony>-_WKA}2EFN5f*8YLI8ZmnM_5$87 zVQi$}DMp4xc(6jpr|vMqDIm^lyuk#Ub2U^wF*GsnH#nCEyN3FPhcyUq54;xh9)78{ zJqlP5Rti2{+Jxid)tdQf=9JZWE@{rVQ=X?5O{adb5 zoqY7Gx46g^97auxAp`0Rnz0}Z!c{Yc-HLueIUc-(58&Y7v~?BsWGQdp5iF0VM<3$-gmT+yM`QK4xMnp+tksqOR_^$fDd+5EJJX<1ptOBOVZghDODuvZ;?2 z@_Qt5*qZC!PeS6JSb24IlVal%YsxckPl6pmdpRSc7C1dp+=6m;YNMCqE`4w|R9g*k zeUs&`@GDCHr4cNDNtB{|wpt^oP=N|a{s zr462ePi>iY>wc4o?60w4Xv7*@h){|zDyl!_qn}6R)`t}sRrF^0-*P*3nYe%OfJb=H_HZ}gDb$4aU&J!l%NkTBt(mYf9 zHu^bnPmWGm3Y=lD`eoVVH%~f&cehxQjV_u$B4cA=S9%W!r+Uq#v_EuxrTAL_ zaVL{9gk@*D*t|Z@bF`X+2gCa`2M31RuTBjr?-%J4pFOc5(90P=X^o7%ms+#NWMD&z zMABe#x0iP&qgh3UG{ek2JUvZ!HK*?I=#)M0nke}L>#>j-SlCgo zd?F#yqG--dC>{lp7Fqn+b)pDb>Kd1yJ6gouJ#&y$YfCe?>xsj~mCjK>MC_yduJ~M* zxWhGhG50J78m92SGAy}CoG~-moC#rfRA2?DU!GaUDEP$tl*hRf`(?ssNrzI+2J#ex zLPOyZrQD>Es57n8(R-q<*+qSDZhhO|odUl7A7eS$ahP)MO?Vz@DyyUtilfERL99yI z_#Qtd|HK?67kXw#D;ix@GtlSm_iaaD=`enU?RDe>0+zhq%;lAZeoYFM)%{mH#36|2 zr=5Gby@sraqnWt1WZKcUB+G4OB1s;N2CITCA1A)m`{qGoKi5*iD}WmJ?v455ven~Z zs14Zo#PMInW6~nxcaMV8P^G2}{hB>oY*Uwo_Lm zuKV-osj6H4cQ)w>6+OMZbD9)Nt7;%{228F_Hl;OL$RJi_Y|q*j=!i2YqwNhv>u0?m zm?wVx7{UKmaB;P`220DFKFD36aIKY!j;_LKf`y@Wkv3BT8ljR-GY6eqy@vF_2l&Ww z$VFD4oSdFYuTI9h=eW784btlhx*g1S!uiAN z95lAXOvE+A!f3mE!*h9<9&o7B;VRLr6`zhy`+IEc}<`dHGw<(#Elhb6^v15*OLogjSmhIzm~`C^t*c@H00DbF?32*TWNDD@~Y8C%iI?2v5&4G zhskplm8bz&X4osNNTmAwxcMyhS|i$QQrQs>i`8qyXWO&XekHp-G}Z7xA&ma2($Zq)c7W7ch>Z| zlAg3p9hS`e36n7dzavii)|8JD*8LGViPNwf&cEYut@g1)JDDTs!o@w^Z{)R^8DDVw zKZpS@)332WeaRQQfW3(N(uAe(>bCqipI7ZY`eXaafZ18G#4&}S zf+53(lUx)w+3+;XNFbkC~s#4PRW`HZ*B;?vG)$5N%JEwTzx4FZT(E@2@ z(Qw-iO#>IVwlWa$UO1B0n-X-Y*pivUd-j=NIWa%-HF4WkVlgEskVyibXST+X7B~EP z?wEen?+(;xW=uP@UaLQYP94*)jdnv*O^}|&!z)2D zE+6A9Ys{R14l^lE%jwcfanM{L4aFc&kJ=}BI>%~Xyx{p2tG~-qUjOWJ z#K11gW&8N*m4*gk%XB@NLG4~pch0etv3;n~596~fa2Z?l*WnqRdXMv8!S5bLMjTq9^5T+vAj zzzN%#7GmSQS(FkOXfPs^rCeMFXGbnzjX3CSLbpGSn zU*r?oi!zAO0~7VcnPvl)s95UgtTqTN9?JBS5BIu7JPenrmR|71*70ztvt|o99o=~# zC{qfRg5~&LFCXQ>3W^Kmoo81&jX~ocHao!_EuG+|d$C-FXEmHTo?h-rtW) zd#B!$!FtxEBU<&I-qDr6qR4fgV98NfDrAKfB|qEXuas_L&rbt*(q#1YXxkk*j~79c zm^iGPh?PFSYwOuYrHCntrl#h^8nU3XOiPoqr=sfWMJH0e6=-u!B+3$E#k>idpgHEH zfGRtGuNdwId6Fvjh~q0B^9#yQou{cK5s%>yr>2DUS{>P2X6^$A*3$WlgDO8Io8on^xg?w(r6`xwyv@FK;jHZXy^fD^m=uXn4>! zs<@*^TIR>7n1(0EOjR=^As?_veAl<$2!2fQy5S+u1h<@c%QqXxr61kipK`Ol;-0e3YqY|AWsDt4J$%==Mrtr&h zTX*PvefK9NmhPiLpL<+5$koHn?AL6s%ZM)eEP4#sbve~PybPFqAn9GKT`=rb6&-hl z>PHRstg{m(S01lDi{-k*L|uT8A%3vV=~P6gOs#9kpmrh{D-Wx0;Y|Tc+2ZlAHYKQ1 zamM2JUX^1u7rB04FcDpKdx|7}MZP@~JhaqLn#*nNoqDxC&#G(QYs+WG#l27B3ZE(D zO@MQeqOfUbYbR)69V8Pi&-YH;vmgB^?O{&b1+2TOWr zK*68p>YGpa_H5YsHULG^?$j4;(tG&pVNNN*EKE!w%u#gF<%S~163973d3go&rh8~< zLLcN!;T%B(2g!eN4!n-KuR(py*$&M7aS})I%*fUza32Md3p1yTM#G`vC6Emj!^Fo2 zC2Fzm^|IXl(a@?>siyp?A2*KoK49=`cQF{uvW-_ zqihM{V~VT4o&FwIT94!H&ea-&{3g;XPI70I@M#w}&t9$8!M9&q+d=30Q-Jvd*i}dy^D2i7S}qs?=!?TP z;c=h?y=1drbHEJ=WM!6G`pYk3JVbUcYt>1fL(St7(*wIsQymuYEX}XLdeV21fJo3*T1~BdVDq&I>jofM z-XJ^K=py|iBf`UZW}A~pnCKUf3Gx8!>%6+?UOlp!%*uc(GC0Ygq6SY^pR(`HHscGV zHms%RCzPDDy*)Z|tW5~k%s1n4kk+zn5>l6{g-B7sl9BH2+-+{2I3N5@J|_g)+V#0s z72p_YBHlXYwnz`=C#)|Iz}49ePs7=0TajKVYX3^>(cD}}Y7;`hJ|tOmd;t4UmH-B0 z`yBS_M{*5Wr3gJ?e?KeR>v!1nMurVSuKYsH)q^11kVeGM;hDp4d{8uhpqpqW&eUsp zs*p=k0?7#xuEBt#?cQ2Qf<74Y|Kw=eJbE~3dVn?T+dHjS(CQLMx z$_jvA7(5<;WZTKFxC^_YSxz42m|-N9)~)*AIrdxN>Ofa&doC{c+UNz%A)y~9oY&ji zn?Tt0amMCcI>@wmL$`fq)6*2_t2%c%b6dPTufajtYpbP>U^akw<-o;%e6@bREmLvN z?q>`i__gIsE|i5YA5D$*?wvU{8&Nz}oaS;@|1-TBSCpSWAC~7nc9g^7@xA;wIkb3- zt`-a~4pmU_TwCo=qCXV*N6)6ecCBIpdZq^G8FV!yBBJ!d?d=B$ zhK6*Z{Zj77@mduHTVU~t9u9o;Wb3bpA`M7YSmM^&TKHIQh^bU7iwU=b@EZm=`JZ>Q z9jyr1w<8EN6^yWjkkG5%-mb1EVAH8Vx}sVPd)4d2f4zM!r54z+7I!T+hO*KV<(HXR z5uwPCL-d^Ws53kF!fDRpjfHt}lmufos4}a$hHSF^JRE3=R)C>tNvW@6LwvaWDdA43 z{B}{SpBM-U3Clse_nBJMP|UdQ_QX`UcG|!`WN{Jmr|++kM4DBw;xGZ|!`2VPBn8Mz z4nGS3`2vx75Rh+sGG?5dSR+wcSs76!98M;mCsvVSTxMGwT?|H$9IO}N6ez6JK1Yp= z@F)&uF0Jw-FJK#g!)w}A=}*q<4vCV_h3hdeFmW0+1i(s&9XIU%F|_>s(=hIPuDYi8 z-DW&Rk7_K2usD8-`JgYNLeOCK;tg)?7mkdCHUT>e>oaF! zDKt5qJ{5Yu-ZU`zzvXSc{1S=!XFlUkxb3E~ar1{SG6;@Y|IxER^wZuj)4-*zeo$YF zJcPX0#dKI6FW=p2yJQmyGOe-?uB`kH7-vtSij)@qks&_w$Be%j&p!zX0ro&)M1Z}N zS9r-`aC(6C#236mRsV~fNo1Su|IEi$E5dORds2Y$ys-;#Uy0$ z+XqSN)x=}}1LI_G{1YNcXX?EH+GBdW{4>8~fEP6mcM(6C0*Q|0aaWPSX+LfyxzJ`k znGwb~3B^3xA^U$slq*U9F~-o2UWKWE2q=_OK`H!o?+XM3$+Hwwr|`xcmE$rvUeZa- zo}h^Irx_6=Ua{sj9}10o2RYP8oN0i)S4k1VGtG7Z(1)-3b?cCH1RG_nE3 zo3&E4Vr!$f*VjQ>--^H6?9kn8I_g6_(S;NQS3CGOoFV{m>vpL+UN>=Jn>@VFY7Whp zEf5qbu}TL*XE1cpIgp6W3q#Xoxv}sFvaVVGPE1Lv7={Gr6&Dxh2J8%ndK4WkusWO| zJPiJ=iv@#+`wJSg$DX6KJlm- z(wGd)LLERqF@7{d5w}I}2T*84enCDXc<_K8zXSx0=pYs6bF0tQR$)C`pkefk&s)_@ ziB<~d=};o_VzzNBKC_Eae{zvG5EbT6h+Opu!Y1Q-J+Td)Rc$Eb)MmhxY0(knJIve% zd})Rxu&C(I4Zq>iJNoWONEFg}NwUlga6SP17HZL}d{%T;*Vfnv6Wj1T z{&m&cERJ7tb*emH06mjdqb-k51VPT{l81Gyfi#0o_F)bG1(8O;mD zM)Q~NhpcW-BVr}hG%2=!%;u1*r1R8P_=Ta2O2Y~=Gr4*cVq^|7=G^I!2;)Bis2atO zx#M^?E0R-y6P(aHB_w(Mps4D*Hg zy_s+Ad6=)&OFmb#6GEZd(xp07u6t9VLk6x{m7}E`x>Eb`B0mhiKiN;HR%g+ii}3$K zNpg!#zuE}tfCYyH{{~2EWfX4y9#N#Up96{jN*wzc-G%hp@i6716_=lfGZz>5+8LjV zoF)p$#RM4}BGn$@X`TV(ySEoz^%~&7g?U7+=L#Y31p~wxiztFS2mw%6B85jiPj~)>l zm(K*rhQ=00Jq?JX+u|Fk{f|`tisIEamBaQAvKfqGV@RJr$W}2iZksL>9LmX9P^k#( zi@_Z)83X*g^P9#c8@X7slWGco-)CKEP_z-sb20IE1UJxL}MUh%I~@h42`Sv!LEjJ9^jMWqSMmV5pWo&{B4Aff8ayj~1|GGfq2|XOHU4(NGIY zM3UR*O#RaWK^>o~jGyB~<{J6kCfh$gGzUL)^xG}{H&}~*t&rtbv7grSB`l1sYc9`I zuh~L;*jQa8R-^o(`{0aKq(PlSb^|O^saPqZy{2YEVbg9}6ylSOP!mi?+gix^&7L$h z)lu^j62dtCCtbtbiBI`l&U?l7$mr$w$F7X<#k)YIr)WU>nx||s_HTO&Y0_T%{&oh7 zm?cUsefk@sw-5?<*?cqGD1TS^VdoEbDhZei?k0@pByrs-uDvG2+x`CW?3rE9>~lYN ze5oM#sECV=>2a6hSL-<$*keLMjpR4x=8TS=q`Q(1isq%bZDx_#o9CIn$Kl|ukNNmY6M>-$bsGM*bR3m~fc`@j-8U0ki*Wq00S1^D)DANHHS# z;p!GOsG==k!e+cDqY;njvQ>m=fcU<6UGW6_7C#ZBEtn7IF?WX|Tfy1*goH2<1GZi~ zJK<;^Lc;pIlPL6mLZGmRjJa)mT-eZgrin@W0hDi9d<*PZ&W_ zcO^3!p*7@@aC(eI}IJD;tdF=8j_94xX3Vk z{G_GbGmi#<;peV0|G;tfnWt#~_}jeI<+le$$R=zD7{0%XUwf)}+3?B8AoLZG`_>5Q zzI7@VNcW7Aty22s4_d~Re?-R`JN%gc9{cn^+524o1vdQx8UftK!KJekvGp(c-sK?L zY|+d9p`p$Y)tN((zHypUB1%NagQkh^I1NbM?j5PbtI1VT`CHNR&Ilhjz$v3EV)1*T zY5XN#oGb$34Ta*>y%)qR(P3F$W^(yUj`_#{Wqx27(QrWKtnVYp zVzQ-~7;Z)~m-I3&1Wd&ZYHj=hyrciBW`}5%Q2SyqWUFx7>;7$L+-d8xDb*_GDiynF zp2;bM4`BcaSj1dODzH+1?=q@v>L%eYs!|}eb3d64aU|i0y#NStCS$= ziqLi|zobGIFRvyrov47*vc)U+{TS$I=z1!2aci1bg{nDlUn?fTDXo&qb0Fl0$EJx8 z^gVv5<7_UVqoDxc!&Y`KaYKH-Kc`@1bSF!kfkAJ=N2$TZKwZ&%Y(^4YdGf>JnbDW) zhPgd3|9|>d0={xB|7HZxz@rKf)(`tfb+M!TaBRuj{j?OOAyq6>tG4 z{))uly7vp4O=BFahQDg~*06XjzpsShpE1{%r4{;p^1C;gt5-9+&f5=nWrM@!oM_V0 z{$DPBk^e{ICjd2;9v>e+ib@pu#mRn55SQT68;Zc zgH(D6l!Y(+D-^GqoKX>BF+7*2X_){g2Wkh!3%M8>6BDUoH?VQg#M}=x%ols_j5=zs zTYG^&AwD%}d(IlHFjy>~127cIi_OKI1uGEL3RV>utRLTfj2 zP@!shrhyn2mzNSw;)LmH=*Gr$d_{B`sEqIJy%u%{2M2pWDwKGTuObiJ)-ezkKRf77 zhKWB&>N+ob16Z>oU2kV=h%pR5_!VS?oPZGa8VG)w951t7UB%$P!Av|1&Z;_2(7O9R zt<|9z9xB5b)a@5mMEOy3`JHiT#}RR1iN}*wvwttlo&X4SG3)2Vd7mnnMiauGQ3=}h zmd{g6)q(dj00jmZl%Q1 zNCvlhSEbZYY1UVfs~5kq+<3&qL{hv(@#M)9zxq*RRhw=)FZ%ZOB185rbJA(HCmhh9 zbDnjXA43qANEsT%q!ANS1y%=zp%i)y{F1+Y6sNNqrr8v%ti_7YRb}n%X+gY~>agR; zGm5D&UtIr@QdfuAOu#X@UW~wL?SsK5c4f=<9`3mI2WyLShgv1>P%%6bhq@2ygFw*L^_0%4Rbe#oM;U(3?r_?JFV zlTaxr@Tolv0T+7t)hYs*A#%+l3rLQ;ua*Mj6qS_~?2W*xu%qNrn#7t3h=_uK_?xA= zLmKO4sG}nQMk3DRS`4+gzj9K=@^Jd&&=|#F0t_rHD6nvh7HOcYt-!a>_d;bAjw&X0;*ye# z_e+6l24-ehy)l^`3-|GU0DpsyD2>;8H+`CV@#e@$uS}qECPU~0+PI~>Xu?dBK5>)S z{fE$A=?U~mBymk?X?ifRLjOdi`JSKVnT%_#K1MI)!Xa?aPvWV$FG@LBK?axT>@AjQMS4sP|0CaUcRI6V}bG;xdQg|+(^6_S(} z!;#&74qTTK^Y@3utLy8qX$gS9FuC4w;ZPmF30_SdvH*;thb8i}t{&^fSj#(ngdL4s z#9gY`cxP`9w${CPlKfZi>l?Y?))#EB^=sz*E>J$H=Z(PF4|pTO!Y~vR9SFq4#ZAAx z)KFIsy9q#^Rc?tDt^E|2oOn9=R@Bwj!0h{kjohq)GTRz0HeImb)W^3#uYWHrTpt!~ zug5R}SovEIm$7QOnrggOsPgi? z#kuxdP>qU=D>!W~F5BPm=453Nuv%N!Io^_DVPVtE|7qwNvYP+OL&;a93+D+5`7Lj; z-f2vbxCVv@7r)v?r+wjomv(e)KjP+wWMpK-rJyNm5FQ^N$7f{t1F=3nRo79jhIgSB z%*Di=M+S1;fn%=gupN;&Uf#RD#TMq~E4%fPLT-1+YUvW89Wu{6PIq;d-=27g|75#R zX^)i&2xK#tT=fFOs4odbgE}`f__S}P=f11$y+kObA}+S*beX-9mgxxzQ3q$fu|WH? zJRg)-3UBpp!ScPheS3I#_~i6r?|DPgOXQcr zjt8H$?R|kPWG#ajT1cE$8L8CeOJmO_DM_~@X#3QWivEGckUJ@(0uJlxm@EiA_qBXs ztKVO%=Iq?|uL?A2qqnE?8I58x>}nt0BoH8XhTwH2Jdek-YjwFEuEwj-+?#I3*DW`y zZA>pLL|-3CjW)hBlPOk`b}>6Ro{VMsSkwR}Ia-ZpUT@E>ju6zNO1586HnD4b97SXN zK65VI55zTM9$-En%LpoG>_)~V6HZ8V`eRQgDLNqB5AVAo%y1KL1qZ&y4t&fs_ER-z zscuZ~f76k+QJ5?DmpPG}Z3CvRe%#FLK}hXhExItGlM_NCa}A0NRf38+Qn700pbHS&`J5WC$nh+0+HEG`O34c$K_pX z>z%jQBn&P`G)K^-SEE1yPyx2D#upI2azAq?Izi*M7feiZS}hB}sKGuozH^F~u3GQO z@_fq&23jzvC?DtcJX{YnF>@Vh<{E^~F1Fg~ViwCJTVtQ8d9jA~5vKLF4b3W7OD5X5 z>yh)RYBdA5i*Ylz``avPTzc;G>W?XNtu!&LmHC|-uE#YyF@#5Zz=fl>AX8hYLs!tB})k72cNytchU=!TdVGqbb2#<^%JiDtZQE&S$}Ha z+`!0ff4Ws$tHNkbi9CD<6XUQK#PCXP>`Hu+<-NuTopnkGaIG$b0(@~TiEA}ezz?c# z0N{O9e(~dA7W+(iYPhk~Uas_hzb`rI4AJO`ky+I=n@>Gm2}MDNT6at68vwk?jas3oVf72fBO3#}0K z;PcW>&dGOs;ZNvWrB7e%QE{10j{2{qaI>?sQvnNm=gzwu!baBvEuR(ST+fO55ZL&^Frw-DPA2rg z_z@hK+LUks^2O&Sh*8|4y6U?XNCb zjF`?=l)^mor?LG{biWzc7rkWNFu?>HWHa;+Ck;aFUGVjNfbJ_4wfd&FuvfzEj^WJP zjo|j`1l|@;lo(}+GaJ?Kzq4{?A@A&#oP0(n2aYBb1;r=CFeVA69a&@PK72^-N7#wJ zRQAPq<(sY%bVuhYwa?ynJJIxGT9SWv))SBEd#rqv_Fc4Lr|s9!upSQJgM0R$)`I?W z)ip3jQ<=?2J|p1Y(*?TJ)SQI)3es=}`v-C2oCwgUEU4>aGCLb{m7V{t(=^Wf!o@EG zS66eWybX@m9RJB~`gp(iZBk5_k0bMm$92-*-KWm-+VP0Aev1A%u6Zes9Yq91b=+Yo z=XNbS-4U%l&gA-jMwC~SA<=N92bbn0+RgUF4v7Go>dH{GHZO>Ltr8v8#+g^eZTL*O zZbv^!RyUHoisk;{sjEajhTe&?xAz=9-O9Rb;M!#)#Nxa5pA7wDpXzksse)Ya#g0Nf|KnKFySj) zr6R~r0Nx38@EY%_Mof)yw*BHDWC<6jE20Vyr%%t!t(T;G3Xcu=>R+!hH@=@nz?I#x zZy)8gML-KT6mh9+oEuJG5aRRZ1x8b|Do^>gk|-@Ls_3)Af+Y(2>QK-tPIsp!_Sb>{ z(*?tN7esNU5%xgtjhT@&_zyD^IapN}_nhW>TL0k#5E841seU_K zgp|T>XY1B@i;lOnW4NH*5jLy*ADgeKTjSL!<3Pf}iWCoEfbz-J{odf7i@~ISTKp|I z@&G^#-kgK3p&?+TRgRI=DTI~pGBfbRPc~aQF!BK_GCf1b-RqQ;f~Bj|R3VouYET76mlS#cVlI$4 zZ#_Ku+KLxk^F@r6h>X*)JrF(FwYC}Pj&xl<5n*8tIGcZ{fE!6p6}`eAP!=FPeQWA0 z(*$;=&U|HE$E}|*ue)=81R053S5@cA;oSaJ9Q^W|O~>Hhjo8AJA&Xg6rJnpjUB{D6 zOhD;omTsdddFoV3o*_%sSg2A@?@u0Sxu4$#lk-*5ZT~63Ha0fibVjbP|IpVz)D7sZ z!`_uDLLW0SG8)}+$JqksGG@O3Y=vD@xWskPCE~2jk>3U62}dALK<<5LMLE6-G&MCX zuI%CBzaR!K6u{1(aBx6hdZ1}){ftXJ-{hJsEwz2s;7sz=KtSv2{HplrOiK?qN6^{X z83Y_S7)$N+a^9Y@w|y;|tEn(n%zMj8gTOe(G)DG_e- z#u$Rfj4?B}j7zfOY81^5h9$>R!3or4TOFZv~(^42sS7L{ny2 z3=!msIJh1ROzgO_JbU}AQeyg>-QF+3Umra{uPFy+IiT=ug(zVc{XINf!v$>qlZuPqq zw6@fxOL`elfsc-F$Os4)jkm5{-$*l9zTu^O9^RMhUj=UOmJ@O z@^2DgZU5BC;3vOus#-IG#l=0HT{ewS6lRsx9514XV4#&gfE1At#w@g6Vu5=$EWTD% zYcVdq)}j~Vc=+zlL!!ZLM|(n6Bggua{KecmBl*bn7{g3gmc%N1MLjE;rn3QRK3nw{ z&ojUgvcENjIj(!hlJuVvlTlO4R-ZQ)T`lE5U|R?GCp_|^6R3&tm9Db+9PO#P!{Zf+ z{WCLExe9pg;2endChjK~{i8NIE;-3qIy{=rhSSve_@gX%sP-BNI7Q`=qDrqo=~!UL zM$BTk$$#B>9s^(gPQ{OIvU^1bE3@+6{y{E3i;}vn z+_1ISi!LLBc5}4lX@iT)UGDz$x3^}%<1|{LTN}DBQDF~mH3*nc7s0u!$;JanxlyxG zwje$l8d?^(c|6qy?Ww$9++LGX2q~<;|15LXF`>3r5EyO(SkO;c|Kfzk9G;%G19q_p zY+eYb)+@Hoae5#CIm{U-6oW&1oAJ&fG!4LZ05m>&(70phB!jWZDi(W&4EhEJI9)P8 zaBu&Klmh!aTBjk)_jVdXKDP>4{_Z_uVE=7LFitDN+!olI)qyn})cGhpn&~WicRkdY zLsMW3Ei>R(@aRDS+0|zI9mEznj>i`uV*V*w8ifj&=Z|lJIjto^*ap6PMA70jZiy)Q z^NxX^zf480-HH$B2Lwx|F6y=Rs>)x|qf?kKv zX}dKaG~X>8#Dn)CqudYzbl>12t05?cW19Vq-L7<%8gVpm_pVh83{7r$)=vk>@#}$Z z_$n{oauL*0M7Dz?5FB!HehJ)^s@%y@;Op?{NZP^SQ9xiI+p41H)eo!2887N{U3Yi) zN1~$W=USSG!&yKD1M`zk+wW~Ad~I*|IzX}U_eTQ7Nk-j7Dmjy}VjqMW5Q_t%qEy~k zT7u~mflkTkf^nNXN)#H33Ns^6g`1k22*}C(+93BfZtruaY!0n40 z*^OK>ffgu@|Ed20I1Yq?hZSU0QL$|bj86~2(G7jz(Umdvr&wpew$-8HE0$()9;G_r zrieZ~Iq_5L{bUA(=7SSioF#@czrmD(A;%^(RO|oE!dJx5Rv)T8HeWLOUf*RG1PS4< zFc;iW(?L>O7_4eEto1;`I1>;SlV@OX1EQ@OeN7fdzi=>a+U)Z&&`WT`2@k3yFWm{?%sESA1^Y^%SLf>~U8ikvxc`DsT}clkvy) z=A@HzoNa@8n7Q-}aQCL?H$~_FP@s7v2 z>|i-*?L#y>*xPpEXOf33^78V)ln3?0XRW??=8KTup}uqS-0osm{NQ&6aD$*zN*Rjy zh20yU?wosJne(us_Tei1%CzvC`KGWok0FB5N>pc$yZ0;%)uEXQtIk+8&m8Vg@GDn0 zK8q^XOvU-!idBEKf0So!a1s1lms^G*XB&vjX-exu%8HF%SegS0ck_dT*wK`N^3jz1 zi#t2m;8ajr8Z{#6`T0ru&-cZTU3A_feH@{d8E#H3nkx`JVUCTsc}qZ8&jD6PiEsZ3 zI=!Y9dy}p!g@L>y*nN+)hnK#8^kT}w(sE^nktN>c7^ieSbc3!>xi!`Ry;UaNRUGn5 zliwMnM<7^=%yBe(EC;Ggp{(M#>tU=FaCfiqY~mIcTdB@UIEu0tSEp1mvxHCH+iAO2UR!5_kMIl-0N1O&<_8L z*V_b-R*_~}e_Vq$t|p|PT$y)Ms!@_6-ZxL)m~v)N3srX25I~{(L75@i*?AB6Sm~}e z=I3UFntbu^u0JD>Q2H0m?Dg{|24wPn5z}aN(GB{(ht#~ZWF1-YswL`as$i#CRl-|I zFbqTm3rP;R@uI%7KbOsK>nud;5Iu!mxo5vi$AF%Y4AA}t*yw}wZ23ItQsYqjXeRY+?0{YSUdOI&mx2v{AB9#q#AJVD^~*F(oAT z`dnuGYvyf}&TY?O8HrKxYwPJ85;S}YWO^ToDJNG~0fBPk*qfUhB^rFf_-7{Ih?vXi zvi=(h#rvrduUnfe=^0KKZ$^)5ipl6=JAPtdX(@;8+Bg=ldV}L!KOgGN1+)qRzKl1$g*$?n`N)a zkdr2cybTiiM_(lx&M5pCR6ESZBkpHg-jlIf`n1W8XcG~_w7|n$t)F;DbgRYpz4>GFcRTj8k){LT4QZ8e3)&Qy zJq-HMk_AXmh(bfuP6Kgr?|Hjv$>Hdau&`~?8P0J=11AO=fzZ3Qjf9N%l1r8EN9OG} z(D;ooNjp0GR0$;dv zUFdff2`oup_GM48o!?tk)pd~i*{s+z&T?}~f1x~M0KbeQ7>Tt|l;%l`OSk!t zk%fq$&vs&bjQl|aVf!}^IsEB4SJXmjSk4auK8b_W!LH|i~h;Knny|Z2+2t{cnfw!zT3`64GbvZ4@HD+9)Z=G zI~Srb>05KYM#?6{z>GW`joSW9fUMeZ62kOo(!A;G*M~E;b}Jqz3^oaR(QPpGOB#Jm^^dW9ag;cQH66&;kgPi?p8WbXlx94YFcUg9$E1P6a<6|v$8eC zlIwmZ&R*LTwMXDpkM`D85%;w@R>l2@$u_^c=*hr)3hXF5w7@av%*;bh9?Yks^~z*P zJgiLL-_3;l)TH%+LZi{KQ8aTdGM+|22#6L2P@Wy7SzevWI)8l|SOo6LcCLw{;ck&L z0Sd)M#YXY4)!cs4-?YSQ@ju#n>#(e%=3V%qkq+qwX_OE!=tfeMMvxGsq*GcNDU}9k zBqgM~Ll6NKDQOUp4(T}a_@3W&-gCas%Re9XUTd$Zd)CZ7*YE;Z?|zZFy@war?w{2e z(mpgaWG*TqLW4y}yE9x87Di8i{5hHDJ%`ESCq1tg_3K<&4ln|wq$mYF%BM0emcmKv zUUvI%vlDB`bawdD8$57JnbI}3y-oNf*F@3sC!Uzsp$pAywKID}M1&|G3v>3QnTceK zocn!^gm+s47WgqCBOC6pFMRdGt*GMFd0_L_KVLmo)+ zw|G_0-C{z=B&5eH9{nSb?tUiR@~d%7rV=6!Dp)2@o?hwM0aP?ZMjs8^Qj_}?lgWe- zxEk0yY;EX0X>w>s{_-UsIp#qh9IbinZNr|DHx4uDg_&>cG@J}PS@@3^sMtvBzkDvMd>i-g)^9j#=n+GahjE71MDRl=AAzW3T zCqsLABR)Q!V_?7y4mu*#VrNkQaUzLusP8OqVe;_$L;z9Bx z>rIy5>$W*5q`3C-i#?+~i|9jL+3S-o+rgxvgbk7QGZnbY9~gt=Wl<65k<=v87h?QS zdWh9?PhFW~Tn;a-TocsNq z-J)lrdd{(10@4V%fIK@rUXJA`bJ9N3(zoiZwi%oMv0OVuh;6@6*+@*U@X*#dNbtq7 zuZM@7gnHEl3Fk^DyPRJdV>gDkOHxM5*Z z9>r>Qi(YQq@=$+IpF~W2!iD@H)t$1&C=B^KAFska3d6t*+YQ{MC9s{Tx+a`D`0u5l zQPwvXe~Q%-qskcX5$+7zHa%8JoGI5eyhO1-`-UhCCm*Tz<(Q%K^(bh+ zdSZMogi4>ZYy?mdzu`qb?(O@nNt*VB6|AdR2g1Rgvuc}Jdc#InXPVi1Gc2ewK&tZ| z89Ys%Ae63tBKek?>cWI+ed%AH-&9@Q^+xfIcWK9)99g)yxy6t#c0>NOwo+hPZT@w* z{qL-aV*=Ihe9y-D&@nL^A}L@$10PgHsosSBLB~Za()cOXNGzzzct4?+aSY4P&kxUU zSFW}@jEaite}40lk&{(mTkk0m{RBGOy9qF3VsA7wD0j~H?G~S_W%Nu`K415fmR%av zr5r{hqae%NOaD%gx+5&zY%^6#_wrz4P55#-U97LKZwJJu&!B=Y6RW}f4l1Hc$lT=N zSgy;}%c7(2kcO)t*mp7P{T)T4p{eVUaf3NK6lv>LXZCSRPy(>A33i8Vy#wRg{x&Jq zmeHIDSFMdQr12A6tI{j)@AlIj zGh$xjj|0N!adFWa$~GbmZq3@HgcMNH2mT_G6%Lo1eLfwdD(-&B6U`v<8+oR+;Qfx` z->vCW%64~2L7X8-)1b5Nwz(@CLvObjNedpNo14?8P!S8s^WgPG zH5cKW)Ed@u%J&K_3wL9F0nErEAN8u+)RBsple7J_+D4=MHbGcH0yy5(jW2ks4eJTD zY!#{ys}y}49)5_&z>gZ@MIQEhWyS1!Pn5%5+JDbxgJ-MufHr|xaV1P$rvcJ8T|8>h z@y$PK^75E&wzFT6Pem>vIF2N__E#e(brc&M5WRNs&yTydSC@8+I%QAV(+wMSX^BSA zz(Z3&l8@u&O-HIl;qIP%=fqp9qs56+$Mw!D4m)iqijbxg#e5_l%(-UY5cKcY}yL?Aj1jcGHua z6`i33^FIn7(cQtZrN5?H`n74{rYZ9d2Ok4tCN_U-}gMLW>$*KBoiApwPxL zs)-g??$>C{sBInoEgIU|A?Hn^iTvqqt3CjQ_+D=41_cM#i`Dq7Yd>l57(uIK&JO8X zUAy-1qR?rzkEgKYu!1I$H>bnb-KYJ1cgV=E536K2{d5juHElt@feOY`AwnFtZeb!U zQTePJsN7@(&>`#=$#ubdB0u7TZ zj}{+2nGf>a4v_D7HNzX(l{ob?r>YJcb5qR(fe>-u@Go6o9j?g1vQU;y?=iTBt{A@O zJ!Ti7QL#HLc=>n7;!>MColcY_K&o^z1uhaHpixV|JD-(I{DyruKI5U@BSe04zd27z zBeD0Z*3pPd83ibC{IOYaadtSrI+7zQp4`_9L1o;Ed1I90eAf2P{!womg@RV2!F;hy zsKe@<6bOrWP%ui?xm&`((a$zx^k2@N$xW#ltz>3BQ8G(nE3}BVU9;YDn-#UyhOK~P zVhtZNletBgP`dcw(y8m6X1~?sUvMcoh4W3fGq^Jb_upcL}!k@$FrX^2W zgKYfcvR&wfBNg4>sPt%12hR$@Yl%T6CTN{!qzkik*zw44X#DI~7aCj$1Ok$WbMC%3 zj2!UH6A~};+vlKcxi`9Fa*b&?=easMI(prv3aI}0C(qmIf2L`~y#Br-ia8VsD=>kO zi$Obm`Ek|Q8^R{=iDXghsap5>Y zv`EZ<{zKu2Y_Jd&!MCG9c5=@yf5OoRW%LFQ z4<0F-nk+P@Pnp#^J6DJbkx{L!(I-(67k0ptr@d&Tuu*L8HgbLoyyto?r<*12`Z!&nbEH4?nGg}W6k<==%`$VSN-C8L=Ovm-M9 zl3rw{209)d>)6S)UHY8GzK;(^;ca@j*?`V$UFXr=hBYyuY+S~Jj7uC^1yv3h0avg8 z!uNSfiTPN_{t`$a+xLEGE={*=222o-eq?~Vuzu|kpp6FOd8};WG5eLV3m@F^FZX?7_MR5b$=MIlYxpm52 zc1@o>d*)&6CaY=i#Cz4OfB(}ocm=1_?6QHff{uBE6F7G*55^M0)jXyk$pK-y;cC^# z5Z+CEoOqNkauQ*uQEOb!m$uPt z3bl}Fk;D#g+ni!??7uFn(8!RPw|*Zfbpwx~X!G4-02RqdGt#6P+_A5V!q^i=9SxkD zGmcE`g|F)DH}QF3pmWcB0O@6f&6xK9)rcRRR@Nw~1tw-H3mOSa0pl zbPsC1>KoB6W4^gHaUVWTeqv~88x>#qGg_>BUfVi9KK69~5LH7Aj(|nr`xD+l>OtS~ zpdfH`ve0HBZ+hMpzWTk^zWO=!_@*1}eGek#M4tA)>tmiJ0vcj-AJG|%%RVvjj#W!Y zOA6Yv#6F*(y|2`|_9JPy6t}}Pb-pjaN3nJJrRzX}6~`^Ejt}@WM4+ap;`?|*;*3Qy zV874=zJGVv#=2R`%cx%Lqdd>%?zsEq8lqI| zV;l(u?~}Ipfy5C&=?AM5q?ebM>B)R=C}88u06vzU<<>`H$aYIUc9n zO0884&GYEF-+znkh8*v|M=Jf-_hPT@rW>-(NY&eHD?E&AWM_hd0}=;R8c_`0FS^Jm z-_2!B5Q)*+NHSs)P({4@#X~#Z%5ZNqE!whn3Q7bt((b$BK)^0jyEScM5+HfdcU&(j5I|_x+qr{S!`gnL3O0V8Gx=4Xhc0B#TNh{%o2assfQ!Sq~U;TYY^7s6a52`dS>oas3 zqk>M%%sM=GK@+32K+yW+4%Gbosw=S}TH6NJke z$t}Z%Ac*~3HfQK=>Q@!iB$E<@;nt5QBkLRNcc6N0HHk{*&3}V;e+UElHb0<4V!-dYzyL!5%!a@!9C6}UfDLo1ej%cl) z)#wrt)4k?i!=O(g)|eh)eW<^7{#zvA!_6s{Qcr_lO^$Rc>KU?XhwoY1g-|RfC1%cD z4SSad5*<(W%!7i0GG%(BEyr5KMDD*;hZrGPDFWvo;vkEap)Z8bV6NJw7gpB=P6ZuVN#)~CMpq|X0!Vrko&%goX9q6lTjg(a&rF=uE~p&CrN zqwqT6=9}2qHgWtPs~}9w0H8SP-meQJ3}*|euvLP${^0=@**VWdJhDAg~ zFq6g{_*F_Q&J+a`|!Yy`Er)cN%O>t;T>#9WWThOmz`;!wU$ebHXQk2 zrlT2`$aK-8^(@y=VSW>f;u_+iw)WxXf*b;ykWc|0wa{|1@B+&k@-H6Y+$%#6Z`MUt z+fKbn$Eqr<(g$_iP?LUga&p)8&-EZ8t|Gkc%2xScyrX6UheohR9A)UA;HTX3XLC zP5eAoR-wzLI1*1N3IfjvI*WpUEM=J0H*{4~vyt2fq;T21GD_hc9n!gwwXLcWe3qL(_Qx6(5tk4*=k(wjR8BY> zJp2y%dwJ#tc&NM_)BXF@Sy@@X=pVR71WU2#=;?)%_Qo_uQVm)_d^IR@n|berG&Ud~ zCWr;68{8G_D)8^vnKQEArmNSZ;^-R){xR)-?>JfLm)D~Dm@AO)GP18QtW)7(!X-Wn zF%?yKUjYa5b@1nGq!QS%)I509^^3jH?RM^BEG<4OTYSm&1>pru!+6xp*=EFGH=?H@ z$@q0>uzhnw<=59iX(>Abm*&I#xA)(xEj@eA!hXhZ@esZ`0fD@S$K4JVW+Hm}=<;&D zppf9*Kd54YYK70>}edeO@5ld+oThVN@MHEEC_wZLuPbj(e z-@JPtRkE>xb?8_S9Y=yh!%*@E*}4|GF+u$H_I7MsLIq1p4rv(~ViM{Ic*|B;=9#I< za2TDqOWzB2>9bqv4lhtNShA_`!os4Z_QQS;2l>a5z)hEXDGu^mh|tv6*TZ2Cv(Ao=#cZ=yB5)e-@wTDoGO&D!4!Bgp*zuc5ld@F&y^b zW}z9gGtc7N7RJHdpFYk0m3Ivx)KsCPt)2bp)3@b>E<^B9Qm{YxiVq(Wp`xJ?Ux#ht(dB-H>%K)TMkxR8HOrq3?~w_Zzp8*RBpD%774~9s zW6E)QEl4gQ=Tj9n#>oltD@yr_LvH?8jjHyfjnRSSFuDrnSTUNry?b0jLdV6+P6;R_ z(nUix{(BTU3f&gdyb`TwUkGRarxrkD_Ux>4XqExO9B8v`!Woe)=d+ZBX6bT22`6Vm z8dklJ+6P;A+Q7zVhsfThBqF~a?VrSvxZ7y|+kK4oIvM$byaj!~lJY|%eY^{C8wx>A zl-j7cYO~rYjf27MBFC8&=+Ga+%T_m0BD`S!c%>Q_2ytQ&3-+`3C{@Q?HMQe% zN^-HySy?%FxtkyJ#l^Q7Y2ZxT74IG@y?E-grDM9<5WTnO6erOY9r-3lo+(SOG9^Qi zFgREW&0VMr@AaRolYW{{kzE(!ieeEMl8_bf zS1>1JTg8okk9&jAI~p;3i(in+ib~}YyZ%Yx*qujAag>;KP2~92`@v}9mb(plM3=k*c!=_Sw zKBbb^6M%B`GGty>)X|^k4U$XU>VA`E9=? zDJbaH&hq9SlvQj0&(M|MPTlpBc^OE6xrnHs4s&0yA?P zsxv;Wuq&s{-&vkMQ`yY=-}QgPa7!JJ!ag$+S$(`WU|W;@;Nf6KQZ}mUu(byYzMJ*x z^P-a>Dhr{`vJjc78xM=GEVL_$nG#a2(hTvW%I~4lAJYaBRi8G_+oN=p2whE=4z}h~ z(_d<5b|ai`rt?@=unyanOlLe^Zq4@^;;YLb@?$~0NqE_4ICbP$tQRO%>}gkhaTUM| z8Q-T{B0)f1KX&@O>&KYIXhS4l>L8LuHiW=4fL&j#Ek&z17V0{6sid)G<>azt^3|YK z`sS4vDe;I*y~i#NIj>y``(<^wS=k`P2|&`OW zUe~UoVwm3#Y$%R<)MioK(K9$z;W-ek&udZk)r<_Y?Y;8MMaE{JqnSJ-N7z{y(4G8w z2V|42r_=-5C2}SvBa1x~Q3w}AjQ60vQn;8G;ao9GFW<`_4o!S>dX!P3)}c_HSBfM` zew|p-21Sx_OV>EX`47B{Vld}Udhwm_(%yAG6pf}Ai(V~I1kS(Fa~TLr4u--0?jErY z?ACSW+%rpN*>E;I9NX`ukm#WKnu3Xw^XS5G}bnJDLP{E%q~6f%@Gucl&DD>#x)7iV$5Cu5&N^?#z$K3xI}D{ieSi1fy@{lBy|^1perS>t&Hm7xu$ae5 zI{ZoQ_CGohd}t46{b@h`>V?Q5;0dmQ_(!*d5fb|4##@7-7ZVDXVV8f3D?5WV*w8%C z{6FOoWUQTC?}?!aAI1`d+h%6#FK{7si56O6I{Gp$F0Lr@k~$@++R~_hMgiYr!{!J` zc(!aDE0XvZx|t{iRy~tc2uBe}G8L&t?-Q=e78h2#nWj7eVFjG?+`GIN zt&bl+W?26gSd@BdI=OyZK{@{zf!cq_EH5u_yq8r>@Rcb4Cx-mupSjK^_mZKC~OR%>@jZEQMw5934GIa?u<} zK2ieq0!-ycIxKf|)C^0xO7q2}Y`=C^wKMDpNaJQXy@Tw9c1 zNN#DKKyLnzro23PLT+gXP=bT%N=3}xFOvT4y({OG*FQJ8w>po2u|*-VJ9UXDyCR4Y z;N)cO`X{@b{wL zxnKP5bN@7qwp^~GLUu3XZ$Tpf^!{;>i7Sn{g4#0xWudH|-bXKz(O04zF6YMc9raH+ z@qC`)XY{hHDR!-t#` z#HXb#9d=}E-s66CI(9U09g};+c$Vqi+qWu@bn}(x#t_n4%QaThB)qCrF~0{%kaS#c z5(`VJOn?vqNl8OX;JO(bOjL;3=l**r!`V=;z!blRua28@?W1D8!?}{A3RdQ`;sBzB z)WfeiL?k68d1%DMLQXIV;O;S*cFB{_>;)3@deplphF_l4@tN0t8ovK+4mYsz?OO}? z(cm;RJXv&Cow2Z)!KIHaWHu_=GOOp$0>8#>oDlY@gP1MlgitY zlv5t~Xxb)j`xS3%+QxH1iosX05@?wvqo(=MoBXVfw18!hyVroy8KmD-v;EaX2YTG^ zgi$_I__UsE8V*hjJsm1Nu@ta?P}sMR%LZN<1O_idYCwpgsy~A;C93B`9s$8LB~090vcK$Ml9Rrq=w9vUdzm_Jl8;~Z^Q#m;EcJ)^T04Tjn-mo>E0&Pw3SMt-S5$5ExU=gl_LD&5Mk+5ir!G%k9UorBWO2gxOV9AuNg(zA8o*- zYi9L_D^1+l$lR6+-IrmBS6t=6gTjYYzlM*y8ud+0vzY{$x5zlVw-FqT6vNq`s75KG z&?Qoc5B7ALx)z`k%3 zC=*5nse=HirV6ji6*6H+dTAxk&zmTw-Ngj{hLuw^u>DSFZuB5J?PI-aMUmA~eS}|9 zE{+)d&K(Q;=X$zUX6xUA_TDd<4V4t|C;B(+rbGELTc)0L6OdptXFRW*C^z{oY&+3_ zOAsn!et=_0DsgUW6!Gm{3^hvzFW%Pa%#npGm1@YWVu01d@tKneUv;zd@i&kk&?J5! zyYvpgt;Z^MLK*08wrTuGzDwcYxovx>h*pTO1rOfVZoiGa+gOk|1qqs_`RUn-$&Wr0 z%42E2e>0k7WyR+wt7(qUKz`jWKJ<@^qU=^4kH^Y={`^^qS4G7g#04Le zs)uCMbacU-;fuq-R*%%)A}JFSVx-wX3ym&_bXA<=wHQ1s_T2x?xtY0i^LXn3QTcr6 zkhd__UV0%P5m&9h@vR<>jXhe%R)rL4y^T8U+>O_349Cl@xNU{ObfMUWpj1-$h_0md zKeX@c$Rvk>3Uu2)mbdfT+_lWDQ$i@x{m5ZxPvB@j((-`fiujT!Cw`QI0cSA5gye>W zUxV>iOpi7G9`RoxVoU#%4&HP+|KfHG3pL#vCDUL#%_au+n>mxAz_b} zdiE>`4nG^S7b3=$Qr-%SX*Nwwju)X#Iyf(W=kH61eh~PPSkC96Qr0Bi82@452%+*Y z-U^!RT>AC_$MErx*SNDWZ#$iS@dU*UFX72!_-&8F7op`&ps#Tpps%v4 zr(c3`ud1G!b4<5)WD=P8RHX8{JCe9F)VObM#YBL!_+D68T!Y6N5}!dj z)Fr-RLVB)ump|)j@%3ixIxcQt1G|xGiV!L)C#q#Jmmn~1zGIJ)ieYbXyg(rqqVzn$ z(5>;thlblyff^<0T+I}^*0a_a@kf(V;lFn2^3rdxv8Bp%L|9-m`7s3;Yto%w|3*7z-5wm46!>0+&y^N*T%HkPT$!=VrAej<~%-nWv) zKF(Tgg`)i-u_+3(-X^4<)C4y8yXT{0x;3UZad2P5Gz;Q}q-qC~e&P5mO%ONwy*~$c zys!(0FovIMeQ0fE06M8?q`YIIfxK$KRoMj_a^u`OW&b7}@uu_y`*yhXr9)LfWmnZ& zC23yVHj3p3EuTM3p3{W}!7VSr_0N62bHXTTEL2@s1tIGY>sh|J3m_J*Q6KP7%eBNi$;Zj2z zurEu&N(22ch>UGxvRU0qhy9MxqDW50|bHLi4eNasm z)3tUY@c&vQAR~hux?*|YTf`ot!?q^{R8GjKP%TZfo>LN{LiZQuJx7=Si3%UnEjCKB zAoa8hcJvSD7s*17IHVeNB_9S@zxDdkoy) zGn!ONe*+;yd$>J&PAa@%f-iM-{PTO zA)$YqR`RlnLt^puj~|+OiSz}Ra3MwiA4F&9jV(E0)K$=jeM61F$(%BT!np=bd*XaNFWbdE6BK?clB<6b_r}@d_CcK=TBkVD;8}&Fvi%L`Ot*#>owJ{MsRg!;hOa$7K%lj%9@1F9;h1M%Jr7l z`T%r4Kq<>VSN|N~Tmedh<86G%xFCD?f!9>J|EYCBQmt`2Dsty0CMI(<$dptwI*ctu zqPf>mNCPtr?4%&C2b2-7m63LFW6Do)nZKtYD9t{$=P7Kx{e@Q64+r>eOtTkTQz!|D zUP7kj7{q)&SFKm=*CbIZp6r|)h=)gHkTpgvfaYEx(VVMTIQo(rQ*|;fd2ZSxawU9v zDE>J3b2Mvxi}tHKii$?O6(1#^$1oz2mLL7n7pbl<5&xkE2JaC2u6VDrsg4bOM*psa zIhjXf#SpB&^FxbxqXhU=#+L>CcsMv+L5?iDDzh`TFq1P>K%lZ*(3R|KT63mA4c?>k zpBiKWE3bLqa>GE9MX*!7tYDj@I~I*FyS5Rs#r-dSmT(2rEk-)G{6YowHcFD_viRR;geoa~)bsH6 z16^<|>CIS>tiz0p=d=yQ^L+TvojxYu#>}-?fv4 zW+GhKe0!wxKJcg{@E8Vz$$?P=RMMRy;9C6D)zq{-)&h`4VwxYVzCR~;*Ps|6sQ`hX zNmi1eBVPK(UvyB|qF{h=iiRziz0us$*Nxx7uvfeqcZthmvv|6=KINU! zSFYH*+=wigfm3A{p7aap;Y1e^p&Yl2Gl%ZN z?>yfTH{CD&*#^pY^_7)UrnPWRxR6ij zAD{9+9?Cv9$oNDjFbWkzFLg~P4<;6KyRiIi3%LaEbCJ9w?aczZuqS?!;(5RnH=tA6 z3G7i?J@2=K~ag! zfBw-QQCy^(jVYxRp7S)L?bEVNY~s^Qfm77X=K9rlgm((?SZAR+I?luObYn*Pk2A`M zf|iA(KN-@+IA5IYd6@6vSiv+l0Ix7$BU_+G5ofn}p((aC_XAroco^{r=#Y7rmR6GY zoXl6@&oK-7g648hZ?|dxk|Rvvc@L|9sr7$QoK}rYg0yKs=hpejN`_xtF;^V)LvbGU z(Ac9-rX?_3>tjGIFX)FNXPp~rw1EW$ilFqa{qqhG1UBn7+TG^jB0Ic6GJZ1+9PE7B zB7{?0P$%tbaeAP4lNLGuIrkrpIYdMJrpH9bC&-8>P6T@BA%i08R?|VaPjV-SZB`#_ z6;3&v95Kl33ew3SQqs5sU`&O~t3W;o#SS+2uh-sr6!;eI3|I@f9U+d4uUIBJHjHgO z=3`~elQ!oYX+Dln`n#a!weq11b+7gRKsN!Dj|SB+cj-W}e_)`(?Sgz`yk6&p$g{93 zA)IrM%kwxFn9+LF-QoXd^dE3bUgpgux{crv`LnY&)6C}xfnYEv1cI8}NZi&E2qR+A zAo;2(&k*uq5YnZcfiN@Dx}(EiloRdKeoa9W{bERux`gA z#2qFFTdy`MC5qNFWe;MHiyqJ?CMDStM01?B#O@^`r2jFXQ}e3Q#{>m=XuVllv=yci zwwsI|G;oLNf`Uk+ESrVUj6vo4BA6kkQO=i%x#%^UN^CAj<;#NcHgk0wCRYy~t_SX# zDKK(yJe0d}0<8iB^x_d@+&{C*3lz-{II+Lz3;KN}fZ>5L(3K)|7c!EAGl`TT9{>9YQmC8h0&qA0|JKqk4^6D@h=>oYzxp)jp9)b=)y8wQk_n~ZxGU#o5ly&Zlv?$wi;Fwo;&v@(OQVVut>=(AesB&O6JC-K4Ptxa zgq}}i%Q;uVVBAp(c(did6hDm4Qopf7j;Jj9l`hu6;<#W51rG%gxtKlwo!=k*d*f|e zB4#+a)-NaedK?;J(4Q4G!(BFHm$$Yjh9K31k#O=r2#RV%D~?T$J0;@v4rY{^J@}*a z%zU{ZQk4`4mj46S{HgCj!M)H2wG1I55@&>qRy-9o$*XIm@>%NJMxsXbZw9`8?fzz7 z?64PyL~A5_gqV7}m#;8z`(cc>C(hHgzTz9~>@@b~*VO{s?Eh)%Cv#`-=X)MiwuKvt z7)9)4WbwMSCbw(0$8qz&7v9qCYrQAb86$CZ+$rdqC7a&c;otPFxa=zd1=~4{8e>3a zS?8?(i6sgNt=#om^9FY!F1<6$-?h^%MzcC&EKJK~V**YLFB ztsut_1`3br9LRQHe$gW&gACHhe*5+f1~;KCu&bKMT&8Zueoh6a{*)1!VtDk28RK7e zUp5H%w<}h`NR~=iGVbl+mfN0C)82de1e7&Mu~1ON_sEaIN3%pAseoyrnCc|5tLvGw zgU39f;M_`@%I7Z#1a`p#St(5tZ0-1|EFT)=Jr4+irrbA~Pk`S%2MBB@Khgca_>I|m z4*EBSlJC1(kWzNrGwpI+uW|P+^?Mpga~l=-Ik=b$k)-SAqx?Ju6_~-neus_W!|f!0rNR!mJ1HIJDk~56n9lPaOi;rwL?PQahei7~-sysIOeD zRl(`WzYazjkvN&5v1g-9J8eX?z^Dvh+egmb6f&hezqrusKE;EWFijN5ezjTd9Mu4X z^Y^eBVdu1IGm8vw0aITe3kza|meR#RO*#`%H?!`jZCm9g_# zEl|n^g%(EG9;B(UyKPUzK%OS6>4}Lt^M4~a=It(6e1egWnH$Y~+%_5|oZ2%180RB} zgmFCNl||$7Ya8{r$n^^d#QP?-V$8jV1J}bVWnnBLgSJFMt3NikL$NKLZ#-kyO3un+ z0{&yXmlC6$E73Y%KQ~7~811-_+pmi9>5=4_nR)$ujNirjnK@2icvHH4Gb%CXTV%3j zfGg=%?xBq&7vx79T$P&yJw^2eb=-2L1GO~t!trTpOxF*mt*4mKlygURa@4Q(+VR0v zK#B=v#O(NJ+fSmy^bWf+)6>w7!>CS-{mCXDy8dPIv|yHcY|%Zy>k+L`#bdhAD=_ zh=xwPT9^vC9&~F)zivkc?&XN#<=OZ5_h=*{`R^C|VufuN_FtEtJ3+bNvp8@-#}xc^ z_!aUF_g>NtBtTkUDyx5KpJZ}FUE;w+E?->O4xScZ%1`EAfuz_DCrxJ=8uLYiL;%8dj1l?i2$HMX*UWtvh1&qSpA+))Iz$fzkO z0+D&qe>g6tq^=H~TgSc)vRs<|yrl^}ve9xtHv>2$4o!lQ*BYn8fh8>k_Vfb*_fNV27Q_1!)0 zcF@>AbZ~HlDW|6!m&1SS{a3TNZUg-}u+Yl3wluHXdYQ>vkSBa7?SKnuxRQG?*P2!{ z9LL+DSgJ{-B3lj4qOX&LloY?+WOtU7+7~OFL#uu;m31zUb8+N@tFqP`4qo2g6xY-W z=dr|LS3|v9{#PTa(HY){D-zEk9?uSiClmxIDxGSCU9ZV_V9W-FcI!e@T_p{0oJm$p zTGsmg!woV2Q5)gxDiuR@=4-UCh$N0(eEUtzvt$o+bt#F2_`ZK=o=~&iLz9et6mKPv zyt7aaDtmyDG)Oo$b8;T6j1=2%w*|CRv10rDV70L+S$wC_oPT5DeMOr z`XBzRtgOi{$E^fRz8G^O+ZexcR8zH&zru9OpIVECPT3CQk=N-y>PT&`S38CC8&W5LJ7BFLZs(z9<62|G&_i2EcCrf?eXr zPg8}AP#Lf@gY(%RS99f6-cntgrT zav0wpnSbz99|{hm$4CBZEB@*@=I)>!WZe9Y<^*q{PHxC$c-OTDIx{k3z~Oj?=jJA> ziv%Ze`{%Fkj*&lL^3_VC*h~EJW2xxnd9$ph0`f<1;50V1E!tZpBixyYbV6H|P<1V@ z_hbxZT=}ZeUD7ibJqOp&G)!aPM6A2acSVV6k5s5VjVGwxP7x&`qhLm zJBlk#Jvxl~Iw-R^Ge7?b5iRi@EyY~}#GC3f&fQjr`@WwnM-*!RxH!rFFh;(KpSM+% z`^U3{)>!E;2;-8&2dGcikjcPjUjt6x(-C`8BVQIDcc6`Bcy|V+>jpw<7H=DE$-72| wZO;%kNZba9(KYv;}zl_*@G(-r9kw;;Oe9p1pIp-rzTq}ZT#wg0k@z4=Kufz diff --git a/man/is_latent_individual.Rd b/man/is_latent_individual.Rd index 0d6df1666..0240eadb5 100644 --- a/man/is_latent_individual.Rd +++ b/man/is_latent_individual.Rd @@ -7,7 +7,7 @@ is_latent_individual(data) } \arguments{ -\item{data}{A \code{data.frame} or \code{data.table} containing line list data} +\item{data}{A \code{data.frame} containing line list data} } \description{ Check if data has the \code{epidist_latent_individual} class diff --git a/man/simulate_exponential_cases.Rd b/man/simulate_exponential_cases.Rd index 3a51cd6ac..f5ad55564 100644 --- a/man/simulate_exponential_cases.Rd +++ b/man/simulate_exponential_cases.Rd @@ -16,7 +16,7 @@ simulate_exponential_cases(r = 0.2, sample_size = 10000, seed, t = 30) \item{t}{Upper bound of the survival time. Defaults to 30.} } \value{ -A \code{data.table} with two columns: \code{case} (case number) and \code{ptime} +A \code{data.frame} with two columns: \code{case} (case number) and \code{ptime} (primary event time). } \description{ diff --git a/man/simulate_gillespie.Rd b/man/simulate_gillespie.Rd index c7582c485..b53340d53 100644 --- a/man/simulate_gillespie.Rd +++ b/man/simulate_gillespie.Rd @@ -18,7 +18,7 @@ simulate_gillespie(r = 0.2, gamma = 1/7, I0 = 50, N = 10000, seed) \item{seed}{The random seed to be used in the simulation process.} } \value{ -A \code{data.table} with two columns: \code{case} (case number) and \code{ptime} +A \code{data.frame} with two columns: \code{case} (case number) and \code{ptime} (primary event time). } \description{ diff --git a/man/simulate_secondary.Rd b/man/simulate_secondary.Rd index f066194ff..e0e097972 100644 --- a/man/simulate_secondary.Rd +++ b/man/simulate_secondary.Rd @@ -14,7 +14,7 @@ simulate_secondary(linelist, dist = rlnorm, ...) \item{...}{Arguments to be passed to the delay distribution function.} } \value{ -A \code{data.table} that augments \code{linelist} with two new columns: \code{delay} +A \code{data.frame} that augments \code{linelist} with two new columns: \code{delay} (secondary event latency) and \code{stime} (the time of the secondary event). } \description{ diff --git a/man/simulate_uniform_cases.Rd b/man/simulate_uniform_cases.Rd index ddc6efbfb..8058f0c17 100644 --- a/man/simulate_uniform_cases.Rd +++ b/man/simulate_uniform_cases.Rd @@ -13,7 +13,7 @@ simulate_uniform_cases(sample_size = 1000, t = 60) times.} } \value{ -A \code{data.table} with two columns: \code{case} (case number) and \code{ptime} +A \code{data.frame} with two columns: \code{case} (case number) and \code{ptime} (primary event time). } \description{ diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index 46d569963..c3d6abed2 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -15,9 +15,8 @@ sim_obs <- simulate_gillespie() |> sdlog = sdlog ) |> observe_process() |> - filter_obs_by_obs_time(obs_time = obs_time) - -sim_obs <- sim_obs[sample(seq_len(.N), sample_size, replace = FALSE)] + filter_obs_by_obs_time(obs_time = obs_time) |> + dplyr::slice_sample(n = sample_size, replace = FALSE) set.seed(101) @@ -34,10 +33,8 @@ sim_obs_gamma <- simulate_gillespie() |> rate = rate ) |> observe_process() |> - filter_obs_by_obs_time(obs_time = obs_time) - -sim_obs_gamma <- - sim_obs_gamma[sample(seq_len(.N), sample_size, replace = FALSE)] + filter_obs_by_obs_time(obs_time = obs_time) |> + dplyr::slice_sample(n = sample_size, replace = FALSE) # Data with a sex difference @@ -65,10 +62,6 @@ sim_obs_sex_f <- dplyr::filter(sim_obs_sex, sex == 1) |> ) sim_obs_sex <- dplyr::bind_rows(sim_obs_sex_m, sim_obs_sex_f) |> - dplyr::arrange(case) - -sim_obs_sex <- sim_obs_sex |> observe_process() |> - filter_obs_by_obs_time(obs_time = obs_time) - -sim_obs_sex <- sim_obs_sex[sample(seq_len(.N), sample_size, replace = FALSE)] + filter_obs_by_obs_time(obs_time = obs_time) |> + dplyr::slice_sample(n = sample_size, replace = FALSE) diff --git a/tests/testthat/test-latent_gamma.R b/tests/testthat/test-latent_gamma.R index 775026e43..d567de880 100644 --- a/tests/testthat/test-latent_gamma.R +++ b/tests/testthat/test-latent_gamma.R @@ -1,4 +1,5 @@ test_that("posterior_predict_latent_gamma outputs positive integers with length equal to draws", { # nolint: line_length_linter. + skip_on_cran() fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) @@ -11,6 +12,7 @@ test_that("posterior_predict_latent_gamma outputs positive integers with length }) test_that("posterior_predict_latent_gamma errors for i out of bounds", { # nolint: line_length_linter. + skip_on_cran() fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) @@ -20,6 +22,7 @@ test_that("posterior_predict_latent_gamma errors for i out of bounds", { # nolin }) test_that("posterior_predict_latent_gamma can generate predictions with no censoring", { # nolint: line_length_linter. + skip_on_cran() fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) @@ -32,6 +35,7 @@ test_that("posterior_predict_latent_gamma can generate predictions with no censo }) test_that("posterior_predict_latent_gamma predicts delays for which the data is in the 95% credible interval", { # nolint: line_length_linter. + skip_on_cran() fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) @@ -52,6 +56,7 @@ test_that("posterior_predict_latent_gamma predicts delays for which the data is }) test_that("posterior_epred_latent_gamma creates a array of non-negative numbers with the correct dimensions", { # nolint: line_length_linter. + skip_on_cran() fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) @@ -64,6 +69,7 @@ test_that("posterior_epred_latent_gamma creates a array of non-negative numbers }) test_that("log_lik_latent_gamma produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter. + skip_on_cran() fit_gamma <- readRDS( system.file("extdata/fit_gamma.rds", package = "epidist") ) diff --git a/tests/testthat/test-latent_individual.R b/tests/testthat/test-latent_individual.R index 5f7a55da6..910f73c80 100644 --- a/tests/testthat/test-latent_individual.R +++ b/tests/testthat/test-latent_individual.R @@ -6,7 +6,6 @@ as_string_formula <- function(formula) { test_that("as_latent_individual.data.frame with default settings an object with the correct classes", { # nolint: line_length_linter. prep_obs <- as_latent_individual(sim_obs) - expect_s3_class(prep_obs, "data.table") expect_s3_class(prep_obs, "data.frame") expect_s3_class(prep_obs, "epidist_latent_individual") }) diff --git a/tests/testthat/test-latent_lognormal.R b/tests/testthat/test-latent_lognormal.R index 1e2efdb45..cc8e963a9 100644 --- a/tests/testthat/test-latent_lognormal.R +++ b/tests/testthat/test-latent_lognormal.R @@ -1,4 +1,5 @@ test_that("posterior_predict_latent_lognormal outputs positive integers with length equal to draws", { # nolint: line_length_linter. + skip_on_cran() fit <- readRDS( system.file("extdata/fit.rds", package = "epidist") ) @@ -11,6 +12,7 @@ test_that("posterior_predict_latent_lognormal outputs positive integers with len }) test_that("posterior_predict_latent_lognormal errors for i out of bounds", { # nolint: line_length_linter. + skip_on_cran() fit <- readRDS( system.file("extdata/fit.rds", package = "epidist") ) @@ -20,6 +22,7 @@ test_that("posterior_predict_latent_lognormal errors for i out of bounds", { # n }) test_that("posterior_predict_latent_lognormal can generate predictions with no censoring", { # nolint: line_length_linter. + skip_on_cran() fit <- readRDS( system.file("extdata/fit.rds", package = "epidist") ) @@ -32,6 +35,7 @@ test_that("posterior_predict_latent_lognormal can generate predictions with no c }) test_that("posterior_predict_latent_lognormal predicts delays for which the data is in the 95% credible interval", { # nolint: line_length_linter. + skip_on_cran() fit <- readRDS( system.file("extdata/fit.rds", package = "epidist") ) @@ -52,6 +56,7 @@ test_that("posterior_predict_latent_lognormal predicts delays for which the data }) test_that("posterior_epred_latent_lognormal creates a array of non-negative numbers with the correct dimensions", { # nolint: line_length_linter. + skip_on_cran() fit <- readRDS( system.file("extdata/fit.rds", package = "epidist") ) @@ -64,6 +69,7 @@ test_that("posterior_epred_latent_lognormal creates a array of non-negative numb }) test_that("log_lik_latent_lognormal produces a vector with length ndraws of finite non-NA numbers", { # nolint: line_length_linter. + skip_on_cran() fit <- readRDS( system.file("extdata/fit.rds", package = "epidist") ) diff --git a/tests/testthat/test-postprocess.R b/tests/testthat/test-postprocess.R index d17d05a99..5ccc13dfc 100644 --- a/tests/testthat/test-postprocess.R +++ b/tests/testthat/test-postprocess.R @@ -9,7 +9,8 @@ test_that("predict_delay_parameters works with NULL newdata and the latent logno output_dir = fs::dir_create(tempfile()) ) pred <- predict_delay_parameters(fit) - expect_s3_class(pred, "data.table") + expect_s3_class(pred, "lognormal_samples") + expect_s3_class(pred, "data.frame") expect_named(pred, c("draw", "index", "mu", "sigma", "mean", "sd")) expect_true(all(pred$mean > 0)) expect_true(all(pred$sd > 0)) @@ -25,10 +26,12 @@ test_that("predict_delay_parameters accepts newdata arguments and prediction by data = prep_obs_sex, formula = brms::bf(mu ~ 1 + sex, sigma ~ 1 + sex), seed = 1, - silent = 2 + silent = 2, + output_dir = fs::dir_create(tempfile()) ) pred_sex <- predict_delay_parameters(fit_sex, prep_obs_sex) - expect_s3_class(pred_sex, "data.table") + expect_s3_class(pred_sex, "lognormal_samples") + expect_s3_class(pred_sex, "data.frame") expect_named(pred_sex, c("draw", "index", "mu", "sigma", "mean", "sd")) expect_true(all(pred_sex$mean > 0)) expect_true(all(pred_sex$sd > 0)) @@ -36,8 +39,9 @@ test_that("predict_delay_parameters accepts newdata arguments and prediction by expect_equal(length(unique(pred_sex$draw)), summary(fit_sex)$total_ndraws) pred_sex_summary <- pred_sex |> + dplyr::mutate(index = as.factor(index)) |> dplyr::left_join( - dplyr::select(data.frame(prep_obs_sex), index = row_id, sex), + dplyr::select(prep_obs_sex, index = row_id, sex), by = "index" ) |> dplyr::group_by(sex) |> @@ -63,13 +67,12 @@ test_that("predict_delay_parameters accepts newdata arguments and prediction by test_that("add_mean_sd.lognormal_samples works with simulated lognormal distribution parameter data", { # nolint: line_length_linter. set.seed(1) - dt <- data.table( + df <- dplyr::tibble( mu = rnorm(n = 100, mean = 1.8, sd = 0.1), sigma = rnorm(n = 100, mean = 0.5, sd = 0.05) ) - class(dt) <- c(class(dt), "lognormal_samples") - x <- add_mean_sd(dt) - expect_s3_class(x, "data.table") + class(df) <- c("lognormal_samples", class(df)) + x <- add_mean_sd(df) expect_named(x, c("mu", "sigma", "mean", "sd")) expect_true(all(x$mean > 0)) expect_true(all(x$sd > 0)) @@ -77,14 +80,13 @@ test_that("add_mean_sd.lognormal_samples works with simulated lognormal distribu test_that("add_mean_sd.gamma_samples works with simulated gamma distribution parameter data", { # nolint: line_length_linter. set.seed(1) - dt <- data.table( + df <- dplyr::tibble( shape = rnorm(n = 100, mean = 2, sd = 0.1), rate = rnorm(n = 100, mean = 3, sd = 0.2) - ) - dt[, mu := shape / rate] - class(dt) <- c(class(dt), "gamma_samples") - x <- add_mean_sd(dt) - expect_s3_class(x, "data.table") + ) |> + dplyr::mutate(mu = shape / rate) + class(df) <- c("gamma_samples", class(df)) + x <- add_mean_sd(df) expect_named(x, c("shape", "rate", "mu", "mean", "sd")) expect_true(all(x$mean > 0)) expect_true(all(x$sd > 0)) diff --git a/vignettes/approx-inference.Rmd b/vignettes/approx-inference.Rmd index 9700e4b8e..d0cd6a102 100644 --- a/vignettes/approx-inference.Rmd +++ b/vignettes/approx-inference.Rmd @@ -110,16 +110,14 @@ sdlog <- 0.5 obs_time <- 25 sample_size <- 200 -obs_cens_trunc <- simulate_gillespie(seed = 101) |> +obs_cens_trunc_samp <- simulate_gillespie(seed = 101) |> simulate_secondary( meanlog = meanlog, sdlog = sdlog ) |> observe_process() |> - filter_obs_by_obs_time(obs_time = obs_time) - -obs_cens_trunc_samp <- - obs_cens_trunc[sample(seq_len(.N), sample_size, replace = FALSE)] + filter_obs_by_obs_time(obs_time = obs_time) |> + slice_sample(n = sample_size, replace = FALSE) ``` We now prepare the data for fitting with the latent individual model, and perform inference with HMC: diff --git a/vignettes/ebola.Rmd b/vignettes/ebola.Rmd index cde97afd4..917bfe6b6 100644 --- a/vignettes/ebola.Rmd +++ b/vignettes/ebola.Rmd @@ -45,7 +45,6 @@ set.seed(123) library(epidist) library(brms) -library(data.table) library(dplyr) library(purrr) library(ggplot2) @@ -151,8 +150,7 @@ That is, $\mu$ and $\sigma$ such that when $x \sim \mathcal{N}(\mu, \sigma)$ the ## Data preparation To prepare the data, we begin by transforming the date columns to `ptime` and `stime` columns for the times of the primary and secondary events respectively. -Both of these columns are relative to the first date of symptom onset in the data. -We also transform the `data.frame` to a `data.table`: +Both of these columns are relative to the first date of symptom onset in the data: ```{r} sierra_leone_ebola_data <- sierra_leone_ebola_data |> @@ -165,8 +163,7 @@ sierra_leone_ebola_data <- sierra_leone_ebola_data |> ptime = as.numeric(date_of_symptom_onset - min(date_of_symptom_onset)), stime = as.numeric(date_of_sample_tested - min(date_of_symptom_onset)) ) |> - select(case, ptime, stime, age, sex, district) |> - as.data.table() + select(case, ptime, stime, age, sex, district) head(sierra_leone_ebola_data) ``` @@ -182,7 +179,7 @@ For the time being, we filter the data to only complete cases (i.e. rows of the ```{r} n <- nrow(obs_cens) -obs_cens <- obs_cens[complete.cases(obs_cens)] +obs_cens <- obs_cens[complete.cases(obs_cens), ] n_complete <- nrow(obs_cens) ``` @@ -196,8 +193,8 @@ Additionally, to speed up computation, we take a random `r 100 * subsample`% sub (In a real analysis, we'd recommend using all the available data). ```{r} -obs_cens <- - obs_cens[sample(seq_len(.N), n_complete * subsample, replace = FALSE)] +obs_cens <- obs_cens |> + slice_sample(n = round(n_complete * subsample), replace = FALSE) ``` ## Model fitting diff --git a/vignettes/epidist.Rmd b/vignettes/epidist.Rmd index 824d1c74f..837892e08 100644 --- a/vignettes/epidist.Rmd +++ b/vignettes/epidist.Rmd @@ -53,11 +53,10 @@ Then, in Section \@ref(fit), we show how `epidist` can be used to accurately est If you would like more technical details, the `epidist` package implements models following best practices as described in @park2024estimating and @charniga2024best. -To run this vignette yourself, as well as the `epidist` package, you will need the `data.table`^[Note that to work with outputs from `epidist` you do not need to use `data.table`: any tool of your preference is suitable!], `purrr`, `ggplot2`, `gt`, and `dplyr` packages installed. +To run this vignette yourself, along with the `epidist` package, you will need the following packages: ```{r load-requirements} library(epidist) -library(data.table) library(purrr) library(ggplot2) library(gt) @@ -66,7 +65,7 @@ library(dplyr) # Example data {#data} -Data should be formatted as a [`data.table`](https://cran.r-project.org/web/packages/data.table/index.html) with the following columns for use within `epidist`: +Data should be formatted as a `data.frame` with the following columns for use within `epidist`: * `case`: The unique case ID. * `ptime`: The time of the primary event. @@ -84,20 +83,21 @@ outbreak <- simulate_gillespie(seed = 101) (ref:outbreak) Early on in the epidemic, there is a high rate of growth in new cases. As more people are infected, the rate of growth slows. (Only every 50th case is shown to avoid over-plotting.) ```{r outbreak, fig.cap="(ref:outbreak)"} -outbreak[case %% 50 == 0, ] |> +outbreak |> + filter(case %% 50 == 0) |> ggplot(aes(x = ptime, y = case)) + geom_point(col = "#56B4E9") + labs(x = "Primary event time (day)", y = "Case number") + theme_minimal() ``` -`outbreak` is a [`data.table`](https://cran.r-project.org/web/packages/data.table/index.html) with the columns `case` and `ptime`. +`outbreak` is a `data.frame` with the columns `case` and `ptime`. Now, to generate secondary events, we will use a lognormal distribution (Figure \@ref(fig:lognormal)) for the delay between primary and secondary events: ```{r} -secondary_dist <- data.table(mu = 1.8, sigma = 0.5) -class(secondary_dist) <- c(class(secondary_dist), "lognormal_samples") +secondary_dist <- data.frame(mu = 1.8, sigma = 0.5) +class(secondary_dist) <- c("lognormal_samples", class(secondary_dist)) secondary_dist <- add_mean_sd(secondary_dist) ``` @@ -130,7 +130,8 @@ obs <- outbreak |> (ref:delay) Secondary events (in green) occur with a delay drawn from the lognormal distribution (Figure \@ref(fig:lognormal)). As with Figure \@ref(fig:outbreak), to make this figure easier to read, only every 50th case is shown. ```{r delay, fig.cap="(ref:delay)"} -obs[case %% 50 == 0, ] |> +obs |> + filter(case %% 50 == 0) |> ggplot(aes(y = case)) + geom_segment( aes(x = ptime, xend = stime, y = case, yend = case), col = "grey" @@ -141,7 +142,7 @@ obs[case %% 50 == 0, ] |> theme_minimal() ``` -`obs` is now a [`data.table`](https://cran.r-project.org/web/packages/data.table/index.html) object with further columns for `delay` and `stime`. +`obs` is now a `data.frame` with further columns for `delay` and `stime`. The secondary event time is simply the primary event time plus the delay: ```{r} @@ -192,8 +193,8 @@ sample_size <- 200 This sample size corresponds to `r 100 * round(sample_size / nrow(obs_cens_trunc), 3)`% of the data. ```{r} -obs_cens_trunc_samp <- - obs_cens_trunc[sample(seq_len(.N), sample_size, replace = FALSE)] +obs_cens_trunc_samp <- obs_cens_trunc |> + slice_sample(n = sample_size, replace = FALSE) ``` Another issue, which `epidist` currently does not account for, is that sometimes only the secondary event might be observed, and not the primary event. @@ -297,9 +298,9 @@ ggplot() + aes(x = x), fun = dlnorm, args = list( - mu = secondary_dist[["mu"]], - sigma = secondary_dist[["sigma"]] - ), + meanlog = secondary_dist[["mu"]], + sdlog = secondary_dist[["sigma"]] + ) ) + geom_function( data = data.frame(x = c(0, 30)), diff --git a/vignettes/faq.Rmd b/vignettes/faq.Rmd index 3ca8185fb..097d7e78f 100644 --- a/vignettes/faq.Rmd +++ b/vignettes/faq.Rmd @@ -37,16 +37,14 @@ sdlog <- 0.5 obs_time <- 25 sample_size <- 200 -obs_cens_trunc <- simulate_gillespie(seed = 101) |> +obs_cens_trunc_samp <- simulate_gillespie(seed = 101) |> simulate_secondary( meanlog = meanlog, sdlog = sdlog ) |> observe_process() |> - filter_obs_by_obs_time(obs_time = obs_time) - -obs_cens_trunc_samp <- - obs_cens_trunc[sample(seq_len(.N), sample_size, replace = FALSE)] + filter_obs_by_obs_time(obs_time = obs_time) |> + slice_sample(n = sample_size, replace = FALSE) data <- as_latent_individual(obs_cens_trunc_samp) fit <- epidist(