From eaefc95849d97d4329e9774ec8fb0b5c09006c31 Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 07:12:49 -0800 Subject: [PATCH 1/8] [feat]: convert FIMS output to the standard format in `convert_output()` convert_output() passed FIMS output through unchanged, so `estimate` and `year` were empty and no stockplotr plot could use it. The new conout_fims() reshapes `FIMS::get_estimates()` output (or a FIMSFit object, or an .rds file with either): - Each data point (index, landings, age and length compositions) is split into "*_observed" and "*_predicted" rows, with the observation error and negative log-likelihood on the observed row, as for SS - `estimated` becomes `estimate`, and the new `start_year` argument converts FIMS year indices to years; `era` is "time", or "fore" after the last year of data. - fims_var_names.csv maps current FIMS labels (e.g., expected_recruitment, numbers_at_age, log_devs) to the standard labels. Also fixes: - convert_output() ignored a user-supplied `model` when `file` is an object. - plot_index() matched every label containing "index" (e.g., index_log_q in Rceattle output or index numbers at age in FIMS output); it now plots only index_observed and index_predicted. - plot_index() failed for a single index fleet, because the fleet facet was added back after process_data() dropped the column. --- R/conout_fims.R | 131 ++++++++++++++++++ R/convert_output.R | 37 +++-- R/plot_index.R | 10 +- inst/resources/fims_var_names.csv | 9 ++ man/conout_bam.Rd | 3 +- man/conout_rceattle.Rd | 3 +- man/conout_ss3.Rd | 3 +- man/convert_output.Rd | 25 +++- tests/testthat/fixtures/fims_output/README.md | 7 + .../fixtures/fims_output/fims_estimates.rds | Bin 0 -> 65224 bytes tests/testthat/test-convert_output_fims.R | 114 +++++++++++++++ 11 files changed, 321 insertions(+), 21 deletions(-) create mode 100644 R/conout_fims.R create mode 100644 tests/testthat/fixtures/fims_output/README.md create mode 100644 tests/testthat/fixtures/fims_output/fims_estimates.rds create mode 100644 tests/testthat/test-convert_output_fims.R diff --git a/R/conout_fims.R b/R/conout_fims.R new file mode 100644 index 00000000..33d8c3bb --- /dev/null +++ b/R/conout_fims.R @@ -0,0 +1,131 @@ +#' Convert FIMS output to the standard output format +#' +#' Reshape the estimates from a FIMS model (`FIMS::get_estimates()`) into the +#' standardized long format returned by [convert_output()]. +#' +#' @param file A tibble returned by `FIMS::get_estimates()`, a `FIMSFit` +#' object, or the path to an .rds file containing either. +#' @param out_new The empty standard output data frame from +#' [convert_output()]. +#' @param start_year Number. First year of the model. FIMS reports years as +#' indices (`year_i`, 1 = first model year), which are converted to years as +#' `year_i + start_year - 1`. NULL keeps the indices. +#' +#' @returns A data frame with the columns of `out_new`. +#' +#' @details Each data point (index, landings, and age and length compositions) +#' is one row in FIMS, with the observed value in `observed` and the expected +#' value in `estimated`. These are split into two rows, labeled +#' "*_observed" and "*_predicted" as for other models (e.g., "index_observed" +#' and "index_predicted"). Observed rows hold the observation error +#' (`exp(log_sd)`, a log-scale standard error for lognormal data) in +#' `uncertainty` and the negative log-likelihood of the observation in +#' `likelihood`. Observations coded as missing (-999) are not observations +#' and are not included. Compositions are kept as counts, as reported by FIMS. +#' +#' All other estimates keep one row each, with the value in `estimate` and +#' its standard error in `uncertainty`. Rows in years after the last year of +#' data (e.g., spawning biomass at the start of the year after the model +#' ends) have `era` "fore"; other rows with a year have `era` "time". +#' +#' @noRd +conout_fims <- function(file, out_new, start_year = NULL) { + dat <- read_fims_estimates(file) + names(dat) <- tolower(names(dat)) + if (!all(c("label", "estimated") %in% names(dat))) { + cli::cli_abort(c( + "{.arg file} is not FIMS output.", + "i" = "Supply the output of {.code FIMS::get_estimates()} or a {.cls FIMSFit} object." + )) + } + col <- function(name) { + if (name %in% names(dat)) dat[[name]] else rep(NA, nrow(dat)) + } + year_i <- suppressWarnings(as.numeric(col("year_i"))) + if (is.null(start_year)) { + cli::cli_alert_info( + "FIMS reports years as indices (1 = first model year); use {.arg start_year} to convert them to years." + ) + year <- year_i + } else { + year <- year_i + start_year - 1 + } + observed <- suppressWarnings(as.numeric(col("observed"))) + observed[observed %in% -999] <- NA + # The last year of data; later years are after the model ends + data_years <- year_i[!is.na(observed)] + last_year <- if (length(data_years) > 0) max(data_years) else max(year_i, na.rm = TRUE) + era <- dplyr::case_when( + is.na(year_i) ~ NA_character_, + year_i > last_year ~ "fore", + TRUE ~ "time" + ) + uncertainty <- suppressWarnings(as.numeric(col("uncertainty"))) + base <- data.frame( + module_name = as.character(col("module_name")), + label = as.character(dat$label), + era = era, + year = year, + initial = suppressWarnings(as.numeric(col("input"))), + estimate = suppressWarnings(as.numeric(dat$estimated)), + uncertainty = uncertainty, + uncertainty_label = ifelse(is.na(uncertainty), NA_character_, "se"), + likelihood = NA_real_, + fleet = as.character(col("fleet")), + age = suppressWarnings(as.numeric(col("age"))), + len_bins = suppressWarnings(as.numeric(col("length"))), + stringsAsFactors = FALSE + ) + + # Data points: split into observed and predicted rows + stream <- dplyr::case_when( + grepl("^index_?expected$", base$label, ignore.case = TRUE) ~ "index", + grepl("^(catch|landings)_?expected$", base$label, ignore.case = TRUE) ~ "landings", + grepl("^age_?comp(osition)?_?expected$", base$label, ignore.case = TRUE) ~ "age_composition", + grepl("^length_?comp(osition)?_?expected$", base$label, ignore.case = TRUE) ~ "length_composition", + TRUE ~ NA_character_ + ) + # Fleets with data for a stream; expected values for other fleets (e.g., + # the expected index of a fishery without an index) keep their FIMS label + has_data <- stats::ave(!is.na(observed), stream, base$fleet, FUN = any) + is_data <- !is.na(stream) & has_data %in% TRUE + obs_rows <- base[is_data & !is.na(observed), , drop = FALSE] + if (nrow(obs_rows) > 0) { + log_sd <- suppressWarnings(as.numeric(col("log_sd")))[is_data & !is.na(observed)] + obs_rows$label <- paste0(stream[is_data & !is.na(observed)], "_observed") + obs_rows$estimate <- observed[is_data & !is.na(observed)] + obs_rows$uncertainty <- exp(log_sd) + obs_rows$uncertainty_label <- ifelse(is.na(log_sd), NA_character_, "se") + obs_rows$initial <- NA_real_ + # Stock Synthesis reports the negative log-likelihood of each observation + obs_rows$likelihood <- -suppressWarnings(as.numeric(col("likelihood")))[is_data & !is.na(observed)] + } + pred_rows <- base[is_data, , drop = FALSE] + pred_rows$label <- paste0(stream[is_data], "_predicted") + out <- rbind(base[!is_data, , drop = FALSE], obs_rows, pred_rows) + + missing <- setdiff(names(out_new), names(out)) + out[missing] <- NA + out[, names(out_new)] +} + +# The estimates of a FIMS model from get_estimates() output, a FIMSFit +# object, or an .rds file with either +read_fims_estimates <- function(file) { + if (is.character(file) && length(file) == 1) { + file <- readRDS(file) + } + if (inherits(file, "FIMSFit")) { + get_estimates <- tryCatch( + getExportedValue("FIMS", "get_estimates"), + error = function(e) { + cli::cli_abort("The {.pkg FIMS} package is needed to read a {.cls FIMSFit} object.") + } + ) + file <- get_estimates(file) + } + if (!is.data.frame(file)) { + cli::cli_abort("{.arg file} must be FIMS output from {.code FIMS::get_estimates()} or a {.cls FIMSFit} object.") + } + as.data.frame(file) +} diff --git a/R/convert_output.R b/R/convert_output.R index 62bda5b8..d3b2ac99 100644 --- a/R/convert_output.R +++ b/R/convert_output.R @@ -3,7 +3,8 @@ #' Format stock assessment output file to a standardized tibble. #' #' @param file Path. Assessment model output file path. Currently accepts SS3 -#' `r4ss::ss_output()` or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object. +#' `r4ss::ss_output()` or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +#' (the tibble returned by `FIMS::get_estimates()`, a `FIMSFit` object, or an .rds file with either). #' @param model String. Assessment model used in evaluation #' #' Default: NULL @@ -20,6 +21,11 @@ #' in the file path provided. #' #' Default: NULL +#' @param start_year Number. First year of the model, for FIMS output only. +#' FIMS reports years as indices (1 = first model year), which are converted +#' to years when `start_year` is given. +#' +#' Default: NULL (years are left as indices) #' #' @author Samantha Bredeck #' @@ -30,6 +36,14 @@ #' version of a model output file. Converted data frame is always returned. #' It will also be saved if save_dir is not NULL. #' +#' For FIMS, each data point (index, landings, and age and length +#' compositions) is split into an observed row (e.g., "index_observed") and a +#' predicted row (e.g., "index_predicted"), as for the other models. Observed +#' rows hold the observation error (a log-scale standard error for lognormal +#' data) in `uncertainty` and the negative log-likelihood in `likelihood`. +#' Other estimates are relabeled with `fims_var_names.csv` (e.g., +#' "expected_recruitment" becomes "recruitment"). +#' #' #' @export #' @@ -46,7 +60,8 @@ convert_output <- function( file, model = NULL, fleet_names = NULL, - save_dir = NULL + save_dir = NULL, + start_year = NULL ) { if (is.character(file) && grepl("Report.sso", file)) { @@ -163,8 +178,12 @@ convert_output <- function( cli::cli_abort("Unknown file type. Please indicate model.") ) } - } else { + } else if (is.null(model)) { model <- switch(class(file)[1], + "FIMSFit" = { + cli::cli_alert_info("Processing FIMS output file...") + "fims" + }, "fims" = { cli::cli_alert_info("Processing FIMS output file...") "fims" @@ -210,17 +229,7 @@ convert_output <- function( cli::cli_abort("JABBA output not currently compatible.") #### FIMS #### } else if (model == "fims") { - if (is.character(file) && length(file) == 1) { - fims_output <- readRDS(file) - } else { - fims_output <- file - } - # TEMPORARY -- add in era to df - if ("era" %notin% colnames(fims_output)) { - fims_output$era <- NA - } - fims_output[setdiff(tolower(names(out_new)), tolower(names(fims_output)))] <- NA - out_new <- fims_output + out_new <- conout_fims(file = file, out_new = out_new, start_year = start_year) #### Rceattle #### } else if (model == "rceattle") { out_new <- conout_rceattle(file = file, out_new = out_new, fleet_names = fleet_names) diff --git a/R/plot_index.R b/R/plot_index.R index 4acb0a13..03552cbc 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -59,7 +59,9 @@ plot_index <- function( # Filter data prepared_data <- filter_data( dat, - label_name = "index", + # Only the observed and predicted index are plotted; other labels that + # contain "index" (e.g., catchability or index numbers at age) are not + label_name = "^index_(observed|predicted)$", era = NULL, geom = "line", # ifelse guarantees the code doesn't miss grouping when label has > 1 value @@ -90,6 +92,10 @@ plot_index <- function( if (!is.null(group)) { facet <- c(facet, group) } + # process_data() drops columns with a single value (e.g., an index from one + # fleet), so facet only by the columns that remain + facet <- intersect(facet, colnames(prepared_data)) + if (length(facet) == 0) facet <- NULL # transform prep data so group_var = label # prepared_data <- prepared_data |> # dplyr::mutate( @@ -117,7 +123,7 @@ plot_index <- function( ) ) # Overwrite facets from base plot_timeseries bc scales need to be free - if ("fleet" %in% colnames(prepared_data)) { + if (!is.null(facet) && "fleet" %in% colnames(prepared_data)) { facet <- paste("~", paste(facet, collapse = " + ")) facet_formula <- stats::reformulate(facet) plt <- plt + ggplot2::facet_wrap(facet_formula, scales = "free", labeller = function(labels) { diff --git a/inst/resources/fims_var_names.csv b/inst/resources/fims_var_names.csv index 3d295826..cc197f51 100644 --- a/inst/resources/fims_var_names.csv +++ b/inst/resources/fims_var_names.csv @@ -44,3 +44,12 @@ selectivity,slope, ,LengthCompositionProportion,length_composition_proportion ,LogRecDev,log_recruitment_deviations ,Q,catchability +Population,expected_recruitment,recruitment +Population,numbers_at_age,abundance +Population,unfished_spawning_biomass,spawning_biomass_unfished +Population,unfished_biomass,biomass_unfished +Population,mortality_M,natural_mortality +Recruitment,log_devs,log_recruitment_deviations +Recruitment,log_rzero,log_recruitment_unfished +Recruitment,log_sd,log_recruitment_sd +Recruitment,logit_steep,logit_steepness diff --git a/man/conout_bam.Rd b/man/conout_bam.Rd index 7b3ee972..8db1ca02 100644 --- a/man/conout_bam.Rd +++ b/man/conout_bam.Rd @@ -8,7 +8,8 @@ conout_bam(file, out_new, fleet_names = NULL) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{fleet_names}{Character vector. Names of fleets in the assessment model as shortened in the output file or names of fleets to replace the native fleet diff --git a/man/conout_rceattle.Rd b/man/conout_rceattle.Rd index 85612b51..e053b2a8 100644 --- a/man/conout_rceattle.Rd +++ b/man/conout_rceattle.Rd @@ -8,7 +8,8 @@ conout_rceattle(file, out_new, fleet_names = NULL) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{fleet_names}{Character vector. Names of fleets in the assessment model as shortened in the output file or names of fleets to replace the native fleet diff --git a/man/conout_ss3.Rd b/man/conout_ss3.Rd index fe5de1e8..303b7fbf 100644 --- a/man/conout_ss3.Rd +++ b/man/conout_ss3.Rd @@ -8,7 +8,8 @@ conout_ss3(file, out_new, fleet_names = NULL) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{fleet_names}{Character vector. Names of fleets in the assessment model as shortened in the output file or names of fleets to replace the native fleet diff --git a/man/convert_output.Rd b/man/convert_output.Rd index 782b50e7..b409306e 100644 --- a/man/convert_output.Rd +++ b/man/convert_output.Rd @@ -4,11 +4,18 @@ \alias{convert_output} \title{Convert Output} \usage{ -convert_output(file, model = NULL, fleet_names = NULL, save_dir = NULL) +convert_output( + file, + model = NULL, + fleet_names = NULL, + save_dir = NULL, + start_year = NULL +) } \arguments{ \item{file}{Path. Assessment model output file path. Currently accepts SS3 -\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS model object.} +\code{r4ss::ss_output()} or Report.sso file (being deprecated), BAM .rda file, Rceattle model object, and FIMS output +(the tibble returned by \code{FIMS::get_estimates()}, a \code{FIMSFit} object, or an .rds file with either).} \item{model}{String. Assessment model used in evaluation @@ -28,6 +35,12 @@ saved file name is not indicated, the file will be saved as "std_output.rda" in the file path provided. Default: NULL} + +\item{start_year}{Number. First year of the model, for FIMS output only. +FIMS reports years as indices (1 = first model year), which are converted +to years when \code{start_year} is given. + +Default: NULL (years are left as indices)} } \value{ A reformatted and standardized tibble of assessment model results @@ -41,6 +54,14 @@ Format stock assessment output file to a standardized tibble. The resulting object is a long-format, transformed and machine readable version of a model output file. Converted data frame is always returned. It will also be saved if save_dir is not NULL. + +For FIMS, each data point (index, landings, and age and length +compositions) is split into an observed row (e.g., "index_observed") and a +predicted row (e.g., "index_predicted"), as for the other models. Observed +rows hold the observation error (a log-scale standard error for lognormal +data) in \code{uncertainty} and the negative log-likelihood in \code{likelihood}. +Other estimates are relabeled with \code{fims_var_names.csv} (e.g., +"expected_recruitment" becomes "recruitment"). } \examples{ \dontrun{ diff --git a/tests/testthat/fixtures/fims_output/README.md b/tests/testthat/fixtures/fims_output/README.md new file mode 100644 index 00000000..a16f69d9 --- /dev/null +++ b/tests/testthat/fixtures/fims_output/README.md @@ -0,0 +1,7 @@ +# FIMS fixture + +`fims_estimates.rds` is a subset of `FIMS::estimates_big_optimized` (the output of +`FIMS::get_estimates()` for the `data_big` example model shipped with FIMS, +NOAA-FIMS/FIMS commit 506b290). Rows kept: every data row (non-missing +`observed` or `likelihood`), all fixed and random effects, and the spawning +biomass, recruitment, biomass, index, and catch derived quantities. diff --git a/tests/testthat/fixtures/fims_output/fims_estimates.rds b/tests/testthat/fixtures/fims_output/fims_estimates.rds new file mode 100644 index 0000000000000000000000000000000000000000..e2f03011ea8c5091265794ce3b16a13e2588b808 GIT binary patch literal 65224 zcmV(#K;*yuH+ooF0004LBHlIv03iVu0001VFXf}=0fX=VT>vZsd>Q$#^c_H`uuVy` zMMEZ)CS}@rudwH3SpP!3p@m_?scEp5S0f0yJ4)$h0 zSy8K6dNmVotlK83&-cKta-_Y4Xlpb)+%7z?;O6fn>itI6Jv&>(RVRQ~S}8Qm6_R&_ zEOaoi-ipp!70jf5?w`*8m7}YAB2}vZVIvwUnk>Pa@#Bah|tXDLFBGmoteT9=FW`=cy6?UVfDzj5fTzM?nH^U$2`6}_SD+Qy`9KYy8+tyaJyf?&95`at%)DfeM?Lm&gnFtTs_An1OVT`S z(2&o%rttP@4XKXY1s(gGxlaRA2V?oO*^x^GNG~XC4SVT)o`oP&;SCQ-04vH zIOVdU3Fj5IeU|RP^t}W>V3|C`fu!^93bNU%hF1ePoLW1iT%ianf~ylwyrUy(J$>+; zuOV1DmsI_Rbtp_PBXdg296~FiSvj!1pp{`{FiuStBcqUCcCl*X1=oiuJ&j$rKZBnQ zRC()7Sfo!3uDlg{I@BJv*a6B!rm%LS*sM@KLv3|hnP9eSY7L8oO}Em9uWM0IIn3Nr zO?aWp%H(572%#5SGR;?7G~iZ0NPcC%`&-Zc>`)n!0dcgV5@YgX9oCzm1o3iVj=K#g z&gGN`M0`s_?&6oH)^uGH|F)88`HFB1H>_kny!4oTMJ2SmcGVY8FD_phW7$Dvq@pxe z<~drJN(eC>nXQY=WZv3^#A;iK_6T7D8`Omd4;dPzRvU(DEXt_xXwP?fhU{5nzU z8t*z9t!+9D8cIDfY3Z(niDphH9%bK=deT0*x#!wrq&5`H=(zvJt#UV)98GT%`LOx& z;&waI?2vPRhGY;t)`$aaqU;$9;8mNM}`oBv1 zGy06Vfu%ptvigt}I@$}%qsYZ?pKr9Tz&e}S{sa;(5B2aOR4uuwKec00l3C~oWpmEP z5Co_I>7GXA1*wpZ)^$sohSd_axL?@2ZA4~?kl5u%3F-`5SWEQTf;@Z6~!E^)&U z3cUEO)h~}Jh|e@-5kDtj)d{Sl1oIm$q0cf?7+tGhkCvpr{*bM9M0#W_uV|%Lj>23Z z{0YM%nTFHAu)XIeTLYj~DSE_a*j}$b9?u=-Een$r##2^nWns1KK_9zHS$6s}G<~z( zqs$j)WU9uHUz*3wboMV3?R6D*rVngftQ}N9zjJ=%P;+e#m~oa%j_E zKAyKX^{e(btw2`Wn#AB6dK3j^`CUGCV^zv}Ko;a1boA)*_2G?_ISB;sZyBXKb}35Q z2qAL07qhzFXfVonHXthug&0Q)TYnS?s+%fndg<+RcEYE@b88hLzydje5;Gungex%> z4OJn=^Q)Nv$zy92Fy8||OQ~v;5fyie3qKvt)G)C;39$Lurl-iGD|WhS6KOUbx)Y#*me_vI8nsrz%;*p{Q2<{2m58wCS9se$0`A&AC5;p5v&E-N zda5E}eA@#&rMaxgLnktDu>BfbE|HoiwRmtv;b zxA11}!GCfEbMXV?QsM4D)1mRP)Tex}VM{DhAAS2e$J{C$jwe2tu{I%cH>o=(-6duv z3xa2_!CTOK>0jqHavpts8PT_TZ%^XqFr-DLPxaPJJ9uR%Q?-b>Z^I1!en0Ua`nQ?}kLOd# zA5U9#$U3(0zkXH7phy>|in+N5h`FTNT4QudR$sS{zwQ%J91mYzPg+Q%GIk$2eh&lR zUpa2RL{O*ic2=*SYm_-xSg_h->u_4f75|H$M8#78WJKmAdF z=a{Ns8hmLi47xMV*Tvi6zy$S+j#kV-@4;EnbD^?y_g~1PYJYB7PXP_?w z5K*XYb|Ht6GdWStzU)dfcWFyxPa7aQEN|RfL6wAS4djDEoPe+(QhcNwimJ}!iqx4DDAm1I zMYx)7@~4rPqwz8!6k97#?6~U{p^PWwRu0f|FLV@3vDEK9f$Bxah;IWe2H%HBMD(Zd zwSUl7soKfNop^l809!+j5zwc3@>lOu{()?+?#)cnqqCWOm5G?kb(LWJ&9N7CRG&7+ z0YE{k&e5k`$()l#K-w}hVo;xWgtq|JfjF2kVkkz#-2>DukuN-Ps7hNYpV?l?rtQ$( z&67GS)ZL;bTN_Sl%|0`r!M?wc$*-4e_e1wTgmIK^?xcK2>C*RU6Bl8we(PuvQF+gg32%%?2ASPb;QrXNS;hK^+fQV2TUX}MVyJ!s9eZmm0B5$ z*HmzVfig>2a(;>bX*mf;X4sPU9xiB*eFlq~cXgB$Pi!H*fb+_f1=KyMF<|;apY=yu z9Om*Dt(Jt22<~|bj+%jpW}wC9VkGPGhbLp1T!F>9X$p~9^N&#{EBev}03ncaR6-@V zBB8Z)umvV($rou~-l~DP%lR>WI?)v{MU4twnfC{+zJrcxZD6Z&Fjy3#Frg9R2}Y5a zNd{}mgv{rAP~rx3pIkFSiIJ3}5fuDWaJs9h3S@FU9o|a%9&m_T4a9^~X&!M=*IOCo z$GUupyW2K?9w4_DsP%})ke1{_-O2Y6yMagrSd2dK%G3Rxyrkf>OYlpTo(&AQ{O7_f zW2Fe>8pQqZhR8Gq!E8c zuiAicZ-p5^uQCACEoTVpJ69D>_egItK@+AB>l@j9y#GH_IdlK~8OS!A1p%FZ3J9Er z9YJ`%Or!gA;(RM@yM56P4T3<}zYh?p46oqwRl(M2_!U8N0OF2cJk=ahW{+MA#eB;L zAJ0$s{p{9&o0}-E!Xh!gp-0mdT_SRUq4M zboqT$#U?L*zG?9X?T}9dMu}QeAHIWDH<+FACU2+bY zPJ;JgQxQnuEhrwFQkrg-TA%8wq}iY<5^7iOt|soq$G$BSY8LQ~kcC1Ib*s6kU?8P{ znc`A@Vo?$udo1~%NoJM|JTx{H_liVnVhmP;KZlduR!I++E?k}Tx7|1?tvJeN~MXAJwXg_+(Q zVdUZdYJ4PI^msN(>oeY`Mf_4nGI>PQFFz`e3UAYEIiaNj>BtLa8{+5xL1|g?Yk7h- zm`K{d$1jsbDD@OH|E%x~cLJ=V)72bC=qfHPt*m?_yj~+ zOf|rbi$I}!5V&QuZTl=z1Z2S{*3<4ZPOREe`$%e!2WksMRm z3eX)@SR~)!^9yvWrg5V0xt@zX*R?h*?y5&76U&Dxe8( zzA}Da^wB^QwU%Dxv!}4qFly#cpj+%zk&77Uh`c7vzQwCVAYSy~JEQpAKncia z7k)w7`Yl8Ah|QHe*oC?gUGEuK6+UQbBgGV*v}3tZ@z+(VfC^`RC|IC z!?lDk zY#I90X7D|= z2f^phjRruI)_t9jRHWHs3jP2|>MgU9aWgIp!Z?%`TIq#P>}Q4dy~+%LO1}TD$n@*6 z@`5-`NlOdb3CZ@Kq}kt=Ke()D_NpW51oahK_@X2sl8jS<;g;7F%KoCccl$jAS_$el zsE?ouNNDW>?(8DdrRYI}VchMn$R8rJ)-+i6K4Qa;Spf5m(4tZ9hsrC?@fc`K7K_Wb zhtsb+gliC*H$vD4{-jdJ*S{KSCTmUB5+ju3CX>7pQ&;#f&BO)-yN!;E&N|LV%x7`v zIR@f=meTif%O1rqS-$1xbg{~4Qjuy4;?|of^ou2%1lz(Er{6zMb>0B9rk~Y6??toLndM-HxcWj&>}a{M=$tBVfAp9>C3kp&fJwHCAV&=LXB(tz zez32a=)$#$-tV8AZ5E3c6~17TVtf=zR&fJcbgPtw;Li-$gBL2tq(01mlS)$cADiFc zkWx5GGt5#Y62UJe=`NUhk3%8Lr_uUl8_)i=2D1i-$yu?f&;2qPr;DJ?SM^3)TcaJe z*cYZWjwJM0H^HgdL`>+cA_ zXiSP_Epm}r05@8ouH;a!@Dvlp@7wE_m80doV}Px-IU)v2pn%hO3JvW(jR)lL*_~M_ zZR(pLa5lBDxtds0Umqp4h7_DW9Pilf(wcdmlA;Th>O13WKBI3}Zr|}OJfzm+@=dfB zz3(#rr!6e3^htk+{c{9I_cIvA2*7Vz@EFVW4;mcPU;5tlNOer|175Kc=5x76- zA`H#b1F-v2BUuEN>n)Q(Z_915__xI^^T_|{nxjgNHh0Eoa(*vf^lpsRk@MJaz`k)eY&ZW)kC054^Z~0#+=&cmy!sK-?eyh6RYz#qBk;wM-aAWV1i_ovlIKvu|t^Lv7~zBO!5Zf8q1zQyE9Tn zo%^bdYo(ksv;9cA6@T@fp$MAEZTbWj zBbH_me|x(@z^T;R_7G5CUq4Zqv2e1j<=!dD#9|sOotqs$C+9@#pE8(d=RUOJ@v6>B zM#a+S`-2S1)jI(%E7{^+ttzTt(n(NV|hLq zUKKG%$~;hoV3Ypi8aK`1w_H`?OhB4P#0#9yqA0ARy<}|^SMZqCiUvcZG6WV@8Gj6DSyQ)Ig zihNJuDVpH846#kSZ(@EAXAD(91sLa!8Fk~0v+9$XtzhQH-OJ7ap}0y4H3(#?v^vyw zK0pgq*B^yb)x%l1I+)OG&og(WX(E4`@T+@Ob)L?LdoT5B!408(0$pH zNlM{BO0uKbcqhPZI*|z(@X(||!$E2=EoW6XN)BP+a>nvf4sX8%^elPDSMxu`y-5;i zF`hXf7NtpHgXM?duViOR#Ew8G96)P$m|FXt z0CFp-z@fGw@8{~=6CW}Ch47<$*+OpX;4gp`4=EG!gJ_ut{{#Sth*fre&H7Xv96n0Cy|J4zY>Ja@YO&bO^ZRTwg59>XzpH*k-&6ioHU0?ymY$+^Q=Tz4 zIpv4$1y4_R^o99oi-f;)saswL;{aw8)%&FzhnraSlS|rX>M(Ri9ZZH@1$8RrSnBFF zNiTYejj)9h3>A2C!W!@cg{BVYw-3sWP(=XVSvf~Sbje&Hg4@Jc-cfH4BOVo`-&V0N zX0~0EtfA&Rh`5I>yz&dfsZsZoEnLd^?WP<{>3@FlqGCP1vnm0xl@>QB_+{U679*Nl zRCh7xOBK!Z?I|XBd&)ZUr7T*r!Q{5)C!7{+srs~lGO#je*UqVX{b%zn%DIH$Ga-xR zg1uvz@swbr(F%q#u)nL3m&f@CdS1~9p@vy^p-oIU;d=0HrG*c65^vCFOs%}_O(3nm z;YM1F^+D-GHO*2asJenCyIeMm8OFewCyQq@Y2QJxaS-9~2rLslLFQ@GZMsA-2X3j* z;~n$Whz^~;ALJ_aQPaOHLl&W!1T%HbZ`uOJEVFN98_nZpYw*tfYHJ!%fMWoyXT<4b z2I-;v;VV_RVzEK3i#dyJ1^5-hHlPOG8|%r?F2&ng^&{lh&n%jQGrY;v#KF*o5yg{! zt*TrA3&Wf&mS$2{bCnokU1#3ueZwyGPBL@mh^!Zsl$CNPbb!80idJ+%&}3(H$Cnpq zjf9cgt6PU|M2X8gq`R^LD4)dUdtzvy7b>P`zz_fyPBy^jU*8&&I+eecUrQ!GsVbHxMM;Nc4zIk-~zr(P^v~&9{(mo8+$Rg_DV&^(cHv zs06{IP&6{z6Ug~mZYHCXp$M*nlh==|5;oEq9$#OKG3eEi<%gQ6_y2wwHl-}q!GjaQ6&IF2ZuwHRCIZqssH;LyPu=FjBaP?BU%B^B~>=3Vxh5b zVA$Kvm)@z`aAL!#MkbFh6F&&eE_qYR=y&Qqm`kU|2@O0JW%}RRlm9E&w!3d>9Ocv4zd?19(0p=36(&JN|I+&fhLJb)s9A4O#tP>_L7B(d&xHsMGZdGFc{ z5lvbJYw>Dn6D0VjXW6{j;pC#!_$z4rVDw(#Iou5s33p6WM4yHVXgaYU%#Ja@o5cv$ zDC32OvXT_7c~wYDOg0Z>=&q5r>qDrG<>xDw>TBrR8-q7fKH~By{!G<#0==acY0=NL z#H^KrVaoSzsmDwlQ$$q$jMQmh>xmConSy)%vrvEl8&}bwg*bdCtw%;Pl>c?c{ z@6X{(HP#IO3{RT^Rsv-r&FtHUg1RyF=%sxwb;yJ2Jh7vC)Srw-bA07`82D$&s1(oP zaT7Jt=ez8Hb0xZnDWhs=>yhUSTI)ygFJPI+Jvq{CR7+d-xUHzFfcqg ztfs@sn38Y3Slx>7Xm@@3z|6%}u1yh|s0^Q`k^EoFmlY`BzN{cq0lt84A&>V}yaJky znJ#iM-x6PB0ZPKRd`1oL?w+{k8%1oaB;YdLe92tiF8TXARcv7bqb2ewxrpnMhOyhV z^tI~87z%h|NuC{38^9LMqV3VZD!4vvlX=$OS^Z6G>pqbkZSH`_@R|9HWN=w}#;5<` zZ29Ist3h`*WIWP%1>e>}qws5~e^u^4K0Zw8GYRr`H;sG&j#f`o^XUV%T!1cvO z(g*Wg=p34C@m)F}H%WTc%d223Wb8kr7>Aa8I?+T`cg3a2#kieJAX>knPzdE)W8C^_ z-;}H_#$Ru!GXYr;ea^L!tD>f+sH0B)IupP8mEMfPP4}467zy6q_l^~K4#-qe!XdEZ zg`1+aO$zY~-LWLl0P{gc}5c!TWVKHa5Z_?m=$EDYb4#m9`<* z@E7k|EmS+Dv5(+xtN%+!uGKXobBNIjuXKCZh7s@ek6bsRx<(Ja45=mwPTS&s(rqbL zDom!-pRjWB%5jYyrk;o9OjX|f*pO$>=Qq$6ir2q6vVBb5_B7|}>%Wt{c&)79ur8M#IKGIh8a^9$byk4}`iO+nWp+&Ybt9W3lz; z6vW`;WvZct%JJ+%=mfuL@oq{BZly}G>hELxGHrT!FRUK%IZVH)(C>QutsH5H{6Foj z9N0dAmyU&7^n=LY17$eSEO1jw_Jx8A)DK~4hHdZGDnGHTDA(=cX zhvj}s^8~2pk}E8(aS@FaU3=p=9HB?+HKN32vCKk5XcWSb?i`Rd$(7nA3M)m~HMDRHTp%n4qj{TrNG z+caFY%>)!e z3`KkNVcE7~tB?_rMNJV7U6lKx40|y^16(#-9f~NE#JA~5N9olwKNm=1DXqlf&iRvu zwWkjcisn@vk}cQZm3gsKJ9DQyI8CZo!4ResmvZiENpA{L#do?-{Xv4eFnJv4!Qy5p zL;Iw$V{mg}#OUzwuH#vt3Hs%YD73rq@UOipmVb-kim6P6q9_)UzR;9sXU}|FR24M} zr&McJv|P-6ejs2-_g;chuXu>?`TQ)47z)qkd=3aKV0JEyQIt=azbDro&oUs4pS`cF za|Os|8kV7$Tr}0oLQqo8;XOkHM!mmy13)Wue3yEg&>$N8po$mA{Red%G;|8Vz>q@~ zlkM0pJ3>xjj4u{FLA{UK@2kY1)AnJgb4 zd~)Nzb|nxV>VBpB5fUCY(MFr89x=fuuEc3VEF;g>MGTtg1ws1YMn-ukAaMImr%6FF_6k6#;6vmRRIXKAlG%~L6U&u?ov;kz>Y_;uWFH_4MM z`rv0v*?jNUGB!Yh<~m?U;#%s`FONM*jE!waCO}ZP^Z4D@Kj2XTpFv=B7!5FX3+(io z&)dMv1MetSw+Y6Hc&--$LX|!li&Blnc|8#gGlSN54O_Sx7F)Ry$j0mu(9=5BUr>aL z<*<9}mRW*@c$Y5F86F7+9^GjhIP?&o-OrLdaQ@~@<-9Qmah_bdpD{qAlL z;3RHJ`ccb!IFE!Ex$MSk*5lP#av5qXdBje}+%?F_Q7LQr(aV6=yq9vQm}rL9St_sL z8-j=)WRjXS^=G?^=69H$GDhBAzZM$~nE|M!RymQmpK;zdM6EGWe0r>Fq7L>ebN5zbBD-te%DkB<+oal3CW50s}#dB zm?eWc563LWLzb|+9^_&#tK4TP4srbtKrqSB^v0noH%Cu|5GRyH=}c(Z5B&Aa z7Dxx#@9Pu2|G)zwas7(cM!CGAPJ&uXVk91jV8_$;Sr@aOV<+R{XUm=ejcA8bhic|n zs}TzeNFvf1g>|8CFRm+M1qlzNxn|he_Yk8a1HpQ*#|R#glAmA-f_vU*BdfaLTP44r z1a*Sf&n*$)FJX9k=6qmyM-kik?f=6M4kOvWy zQD$Tmlhhd&@kJ2vgL0!TLhK|e>##1}5VA?b(J&Bn!GN{!?-4 zBQCtqvOFscy`GFI&XGTZB;rd@fD_ zjOmHM2I4w*`0@(XOD_Aqk(2O=Ey8b3N306^H?G1J3KyWx4xs`_ouN?WqmrH&JAvJn zXdcj$QaG-sw|SVi#`f(PAyjIm+!@m?ILxgVCu}S6phoRLS|~K~6Bxf?g2E4*S}_%u zT~;wn>gpSv{`Or=MplAK+TQO$mjI-L*CuUAf1^sq#9XrEvDP5|EWxF=RZc&w7u9es zQ3wq~5bz&@Nx`6GF)7x*iaDz%aSU1~(#lK#)bw z?eO(M=St0{X45_a@C1`g2mOl)Lmg)mG6#-3Hlfzn$#G!~=DycNyW25Usq4-|9>?0$ z0A?n1D{*02E}3-N=f*R1Hl5hyI!2c~V?!!<_U&kY%6UFUEN`3Bz`A!JJMcMxJh+wh zcMrmJL+P3F%rs6&F236qjdd8mfn)=#TzI zWIsS6sJ?J)1c1H+Y|ilTPM$HPmD-rKQlbXgr&As|KdjWLLS?MB{|zBUq%b3v4*kg#mAJ~$;QU^^-VC8aSV2Qh<1C@s z4Ri_XpP63U%y6<;^N9@uV|d1{-Uh!ilbNO&v>h=_?&S4BYp^}N@;#0TV?UKGJ7Yzb z+Eomy5MfMQw~-?OCn)|b)Y)8H@;)}t3vIUvsyF`vx~%E)x^lN)SYx*>)-&fZwA8le}>_*E@N27OJnQ--RN8XMKmLl8U6!tu|l^MoMz7UbI}9()IYZ8 zEH>j>%B~cuAl&DpiotBjPTx!?Rap{kURJ7qbwth(oqsdk#+y*nk%{j=SU&s#5TBoP z8=-HQW@c37dtp#b;%qIKM26Wg!?zfl+a#X^`{HfW(x_xByRhpqZZJH1Dft(zn+AjA zPaW~qEcM~QdZC1ENflSYKt}SEqZpm_;ST|@jb9^&?)u#UwPfR%F4^yio3T*FwTzCB z{z@z;&7$pq3!}>GJ1t?2WdO4P2`H~F|58x*%U}~ab8=wU8o{EXb#0J=Wm&{8XsP7x zUktY&wN{$5`)>>ZgL~ma<##D`wJ4WqidWYw$p=pahH!;Dgk=%~U3>cjA{$S_7Hz5>mK28AP#ca6#0J zjmmbe``AYi6hat3-~dVLRHGbKDf$InC`0yVQw~FtPq-XOGhLzZ>9;jRzKjv}xPt(T zE;hTAmAHi?FncE;HsYh=OOa}c_58aOUj49I+>re!YJ;Cc2I#W^Nk3&CxA?mlz(RrHm{;BB@@skF00tmqc*Q83v*Mu9>ah>pWCvJ!QdUe%rII0s@Cwi)RNg zhK<7=d7{Bm%|aBdb+^tE|3T`o2i5y%zFKV3A};8e5~}SXu7Sc~|PC|-?BCJ zn+|kEXzJr^CIc;F^P(|;Rq1aN#S~Tx@fY6rS13sO)MxpuW>!h%MWpC z@p!&=%bjQqmT0X}pZC&U)buQd#A-Plk`lLqY?o@;CWiCSp@x#n3g>`Dvs(^J$rmc~ zCiK4d;Pl`L37<tju=AN(2TJRB49dDpenx!%eC?tBP$I9OwSDiwwl2R2fxCeNSZ8d~ zg-TI<7-C)@LDCNv+5zJp62QU^Kx7tZfO5o#R@b|7QuErSgH=h!!ddR9Mflzi<+GmQ znxwd8*`d_yyMOKi73~>qp2}P`z{Du#I@TT}-^)E|xOXXHxw$-La+ING3ciheDt?5u zLlD3XZFh`ewHBha~+w!VSa!;q0+Mtr%yKeocnvroM~!5 zJ8dsIm5UuDY43kgpfN4)gemY^zl97jqgijcDr8D|CCP^i)}hZ?7S@eqsMhMg7m~Wa zB`E!n8AA_Thzptmavk?T{iLF+5#amBFP6cd*ht1=(zS@E3Q=4iF*AdICAww@?w$IIsc02-~<&W3HFo#~t1L%9ROVCFUp8K41>Rj0X7^<3v#YGSg zK$+AnOatGM0l3tnXchQAknRt6$4JFDM+?w`b3Gf}MB0+gqOoT9(2bK*beqIWKW`mE z$bL3LN#4ap(ZN-twcShNqLnUSxt4K2!dv^ihO$Ctq|2s=Jeb+RUYyMXTEz5F(&{z8 zM$b7O#rznrz}J(8O407#GGvXOE)x>vbL$ua>THDoh}pn%15@npVw68nObq``IaI)S zO^H;+TaW)No{ChC#l3?VVQ%yJ$46XlQWlI93$IoKio|b=d?_%6uC7Kx)xe;yg;Vt9 z28#{Ec|KZIgi^*svux(pdB2-<(B(pxz>79Atcwo|X_W{qux0UbH0;EXynT|IJAhHwx91 zf!H*6!v41so91k8ZaGs(I*`5Z21C-yZ)-b1zI6?Iy5qpVOnr=ua=e`i7?7);!{|~u zKUvhMAX_QOZm>GPe&LG=O3cd4c=Cl4mW)>ro*8j=Rb7LEc!7EGe;|MMP`CQUq_n(* z=gGdbX%w%R$c1nrkvQ3;qf;$-#{1{TuGN)PUgs0pCauC^EAzk;$%D)KqG&1hKl__$g=@&w7}i9to}yKKQ)X4$VTF(d zN_o1HkKpmI2sh9OFD=~e3ej#2rXq4H4RQpT^hKH?7M^%N5|ZE-gIkw9x?u8(6@?{s z^At`!F5{a31ro=9TamBZ&@ z_>;1y;w!lDRWnbh1y~(hAh@M1MVBH^3~xPnIBSWmD87TzUbLFUxTj%pO2&#h{jO1X z*=c@fDQ}xRWLo|$LZ70uH`5ejF3CS8;?R^;7_-F;CCZgd?_8s=ZEeJ6L-7Sj zp`CKoB`gXL`LM@KB@R5(LuK1R2J~uBYDFXy<~5t931fwb_jNFB`Hbx%3%_TgbBBUG z=hDBr=<4vsFkj9t0yAg!&rzr6c@EV?`@EMP8W)KmKx@)l;V<-Hcg*sxBcfOL1f#^r zGnqr)7S)-F%c~QzifBFzxX;k?QLH z45|dp@Oe8<x+kO*8c7y6|Bw>8MNw)NiV2 z*9y9+w|^(gtB1FOpo9}m`7C*Cw~H`yAvi<*u~&FkX@eH00R9p40{K@Lk=NwuChMo) z4JZW5g|HunC@&s9`o2#Q?UR#;auD7c8zjY&evgL};5C;DO;>d)bAabe-j1Q=t;`YYR|Ru3agvW)}xnf#8%knhL) zZ}pX6P6Y;vFYtP=i2W;s8gBbQ0#>+nLXsvA2hjszV+Y(VaAqG5cxoT-% z>o`q?)K_EBgq^rg9mrc+^C&znHt!F_rqr-kt{FuEf<^UTm~UiNk<~fq7_&)e_;S_o zAtu+Nq2l1}pas>E0I_O{U4wr{&KXo7pefWk%x$%9`euK`{mtvGK^3TLM9Ic+DzXQSy-lT(y8X) zLG>3gy8#LkKG+sZ2rD6Iz2QfPs^K8Fe94jrwEOK&QJ6w*k=T-CjN5Dcqc!tbbhF+c zTEt1Gw*ZWsk1UAIcvN<8z-0jlgl4{W4Xb@@a{!ZmAz|ASYaoP)DBwv<27@yJp% zz8kwlwnXR7>httrMmADyI4;zcjzOQFJ79@ZM6Q63QYN259Q2pUq6!}$sX3=T@<@e zC7zAR!yS-g(VFGM0Fq0(OuCq4CWW+~#F>!+QS_WUkSTSH!L0_}3E!kKG6<;#38dPf(FIbWuf2#P+uX*?xvijqII9=jh(0Nn zB5O}}?e8fB@kvxobXjTDi_pHUV)BA(*$3yaPgQ^w_ljVFZ-1BbJ|AWqR4i4$g_rFG zB6TF~+sNOKU4bruKUxaLL(i5OVqbZSLW?Zvx-vLaK{;keI%{rSdvPjq%nEZ0L*oK6 z*rW%!pioUR$~DB#U&(8fJu66!JYi4iYLAYDq&lhLku*7|zsW1s5Y>F+eimRmY#c1Hewu8yZU5d0m}L9->P6_ zB0(`~PKb9^#V8%rl56B1=8(Id;UVjhWcFiiiRB&aD57Y{zm+lldhi67dl;@dvwbFf zc2da!=r`;>K4X5(AAxSv=MIw7TzD9 zhJBawIn*Mj{Nc~<9BfBwh}^Vjm#ajL`H>8pJLr5ieFmmi!|I|{0!7a~#*KT~u{86XW zgdyy0AOusSZGP#Z%USFD_EVCQ#c{)oPaJ1>b&}oaLM{L!p-X|lHM;Xm&56iTJ!4F2 zG=}dXPqs({t_Wj!zkRVsxkj-zzng@i6BwKGi8pdavIhCT-6t2 z6MfLoFDQh{(< zb#%fL(GK8D-*O^N0=)H?zzh~^C0F|VMmJfTaD8gxLDM&y5hC$YfgDBe6pHBd9sg-j zJaASX9OwxCOC1wp^#4k9ffS0VR^z&gH~YmO5?MrsqgAe;LTU(vX>c8dk>>AfePeLt z2XVHmWw0!qlTc6JNqzNw`iH_Mioj^N&0tGXbfHtK>mi^9|Y_}iPzDdrjR zb0l^-Xa^@x?4_;GZ_oPZOFd|}tR!1n+s5FJc*{ZDb$1t*%kHKeq}h`nd_de77`+ju z)98M-`Q*Fe%gUKaQloMsg+Zc5&~WsYepNxvXn?)?+)1FdH3(5FFK~M`rU17BHg<7&0bWS!!E3!Qqzcvew22bGl8{H&f4Hjsz8U8 zMQFu@xlJxdaa7!rV%C3bhU^HR2W-3&Jn>EK3yMdP9S(v{R9@UBR_H)G98M%yIefOM?1^Jt0a#*gs^zSX|Q{LeERBm=+n>HFes$G2C&6D&CqG zvQv#+v-MQb#eoKla^HRvC4t(Sy)27n)ZpslResdSeG!!$zPSKFK)$~vmlgq!>}`ip z@_p}TwYhCbi-~+-xz7y=4XHIC0Fs}{M=wLUxxII*Va?2dvioDM<2Jxie#Ev2Ftw*|8i;+A51K@sEV^S{uYS~McDw6 zxvLokGY&r9`}Fh=&aRqq1UmhG=yiGU`|s97s4C-nLE|Y7=Pj8%#D!C7I~Vfga-iTB zlo36lIE3vs*18Co?E&Cq5c8h`uZ*h(z$zMa`!0(larABE3(A?Yzr`zOQ&WScnztkI zb}g**WYR=8X&doW*hA`ZAIAl!6kO|tZKdl06`A zcEp#8#2hCJFroHJTrz8i6%wS&zh?aBCCZK!a9z?Y3~Y&Y2zR1j+YvePW6>-A?j%X%X4RW{I7VBJ$~FH0vONoW;k5$TL>6D)W-+hE^cp9)JMT4++Q92k0%Zkx zPRmpVmOwR=v}g9IWxvjY09XT$!iS#i)6|@{N3Unfdi{BxTT$E@wPJYdAN;Q43PSpa zu`#Lr7tHgHV>7Y{Q-P^J)Yk-FP&P_r?O`=x4}TBI3!_sWO`N`Jh@UN%{Hw@4AXqFn`Q*BW+5*AjeUVSYDfc-g` zMEqS{y;$&PfPs+h#J@q#2Sbf@!d;BB6)r&FmQNe|D-t z84(KUpXU+XG=CNVt;AyGro^c3ket@M(1^_8Y3FdV*NC-?SWKs)1BIO&NGusT{mgznz6o!0iJFIu>Ie&AY-s`bfvBp=dHudlZx#4zXWRJ#EyWI-Xl) zX6;#-FzOv;GLpUhSO~Z!*ggBQ-M=$COy^fPoBy%7W` zi+gOy819$%;+REA?DFdu=gQDGr(F`-zng2`ZLAFxOs7)!a*EuJU@^CI_t3oL7DBeteAC^{KmdiGHR%_PF z(jE2rak@o8+PL|q!nJu!*${CO2Ut%sh&wAygHf@tUed0Kqk=Tbq4);V861Tv;IlY! zL#Rhm&s#cpCfD%QLfI=TArSbxqGo~rY3Q>9=unp(!vgrGtvKGU)iKU+mATqfvNoLR zp1akyY6w8Hg&D#jb(|d$99qR1&nZ@obeffH67Zw)NQ(OvQn5)w_k?Q*dZJtjbHddU zu7__u;|fG%opjsm=EL`nMhqjc&xYIKkPYZN6+>ZA!W^~sG@jHn-EECc(#ex#4-3kU zV z601%?;VXyX`hrygZ4ZrtHzj9ba1{*W75hb`G%89QeWX#zuL5C(2g>R1Q-iHG(ztri z*^&gmdH1jOs<^@3>~fG5&Zuu>rBQ9w^bDE_pn%>_qaAsOVusE3(Eaq}1}&-@Rh zE#KWfxC=IN52u2I2UEKThiPhX?z1`0Qzg)*>BF76m!c64%NTmpg%6TX3LAZ3VMTdz{YzB!t%i#0wt^!Cbgh0`_pLA@cTApLTP(gSCUMjla6_pp|s{zydkFnxk z??YpMph3H>CN96DP=Z@4z6A^V2A<`)ID$9vlD7kaJNC?VzeDJdzB>zzG6LiY=~L+g z7z0q-9GO?!Gd*s%#hsAJ$pG~`s9-jryu9^>E2az@IbQ&qR^EiwPRLHRwaRUu&;jjT^>k*}R(sw6&HvJ?wKie{cgX@h%X6b>1RcH%7)8>x)I!+?0=@bK!~efLR9#?~~vY0I~w55I#HopDmS4e@k1?W`O~H9+K}H|6lPw)h1odj0~Mwwof6gm5)x93SoK zHQk|!(QH^OogR`Gn!b&kdfwZFu^lJbB!+q-ZKX}3!bX8a-#KHvFE|p@C9as&W|0S) z1#?$0^t772zD7Kh8uAxTBTc{#cbP`ecvb?;4AY`EH2Xr67cs%fJ<*yf?=>{g!CwzL7F zhQR@Fw_lBGO|nVw9n_A9k(R-X zogb@1&Me@YcfjdkcovWjPIcj>bhpWukDW1J`I$m|a9kzD>w^J)aHOHs!H?k|L=sW? zV+6bK`}J;))r(sg%&0`zN*uIMptlZ?tl4RM5&TVWc5Ac{L6={ciz6T4@hDxy%2su} zrcV+(*^57-Co@Hwe^y&9(SuNzB2COr%=S+3_#nfWA!lKMzItm6r6;q#$S6qK}oem~f^* zsrf<-6I-S)mkSw|Zgh3122LtgvE8i(W=+}Kli(m-a2U%pOPbgnGKxE!C@~-j?4J+p z+YSpPeK^D*ZN7PJ2E5`$i3QeJ{Zuw-bjL;@4G*WBUV!E7ixQzM6yI~>$Vud2ROqNV zCpEBV0)?O+~5I-&8jfr4)ec7 zS*I(VGDJNkeybK;Ow7`~-#l&AOr=g9Z)t|5#ig?JmbC*)({>OK@bU*)&{^CS<7)OY zoTcE}#tdM1-zY4vD@7YjyOY^$ej>vzW~>}UAENa1=jWSQ-jWh!LijBUC_gk*Iq5DE zE9RWj&ifQqHv#d=EBoUV%Rk?-aX>AwdMK*==YO_35g3@{9*L zF!%#Qde-_s1k%_jT1DO#XE0Mr6@q-;>b6z0i#9Z@HA&n}_^QOHBn>;B3P!*XjLQET zF(9(U`*a?bc_ZWG0^y98PdoM=7M9ORE(_G~z0`X#Ibn78q-i;!LOgNzWw;l|2V5Pp z1Bc8qkjR!L)yh=&|5Cf&MgXYwXs9lsN=TVZn?YP(e8Tb`NP-$$L!# zdeLpM{B&!v_VpwW>OF#1Nh#OW@^F0G39lGu!~2ql)>T79_WbdxOq!f|mtqPdrEIiK ztqf&l9(!hDQ6pWz=%vS&$UTQ?*?LQl`-AVB=?XnPYerFM>PmgWk^-(<1uYwKYTmM; zf-SW7MvgfE&<$CP(G9qCD49Y$OE(14Sj}k0UuJJ~76ZV_%4^j7f*grjFoz zE<|XFH-d2sn~J0=FgKy8wq7h$iTz3X52x;*G=QapW-S5YRsELdiu-RtiT)c8VEuw( zl)R1CnvB7gHc|VM%gO>6)QAUvdHvq7OnCc>jzz4Sp_*9pm1H78XUBMtgO_wJP=y5^ z7~7z+Ue=qpW8d^v&C`et@r#rawY56yQs2t-d+Wo+lpx6Je_acNLQujMyj24LmPtwH zU5b|lNFVr6=^I>i_pV#;h$R?S1Q|W|yP?$cB|(D{t2r9tHw#_NBf&Ma?h8K+*@37M zJBIB|fRT1`Iy^@&+n@H@4w9Dq61$iuHWsrpBawf;3$uFD4V7fyhn^U+$Q? zNELG(6)@!v3;Rx9`4-5Pbw!8ar84+45t#u}x-0vrg55oI2&2f`^q%zVyf57Kv=X*| zsG`)X@^95j-)_cs=|?sV0GI|n_MSm4`&Chx4xIr701%BROn7opPgB}+2<1OL65tVr zpFH>?j+Sqi)CZ=t>9X5y z(01WkLEa72iP1Glbv!g_`Txt@e1!0o50v8Peyk!B;LMM-3D~ZSU`_relcPMd6Q#N$ zVw!_oLafZHI%aHJ1>+B6m3q@;#77Q?%hR%V))3kcWSpYGeJt^*6|hm$#pMK~j7no{ zKQ1(^ieM{Kq27{nWcmLNCfa}D$ZNi3PTDj5wpX^+8~b_t`Bp^v7MI}b#_r_Vs^w2o z5_JfmQqL=LE)HP$wEmHA;OJFM(*45(^h4`QbzJvt6?M?Yx9$xux_Q9KsaSl+Ys7sD znGh}R4Fm+OR(xq57|v~%3hMv|cw)>mZGY3N9UbuSN4u5e?1}`WyAIK+ya;pccg&V( zUz40^-zbJl{+sp$xe1ltKI$&=0Eh-yDB$cyg+}JdWfdVg!vB_Bn%s#d_>39Z%-Tm7 z>2~$cgXMJsCn-^1ksMDV|63}2MXnRVQoW!xn}s@KGbm}}fq>Anjlh^HtcF^n?Y22+ z|7rZ7`ck|LQ|CaH$dmu+E^!Pqc_7SuJ2t~q(RhIAhkq$mNCLK!2^@pO8v~LS(Y6Ae z`9GduQ#Vlhp}nIa$}#F+PMK6LuaJ9jS&%!~<1PL)jz+{kr;|Jo%F56x<%{uQ1B0=@8epj0{%EpI@(CC4)5 zF6?nommSsp>#8LxCvav-S{CI-z~ z3(PHCD*k3tf&3u~)A?QuW!{g_Ld-M0-p9szy0yo+yIiUGEyG=MYrarjeDC%>{Bv^-XsAENOT?o6|45~ z=>S|bpWJ34YS;7@^&HS;b-`}Y)ggdz!FQwrcZg|fzG`y*4~pw+ zs6|zU&``}hdrYvxLiqn4)NsCbBP)~YC7aKQ&79VeT?FQr#z)}TsBmDp$|kO1L!t(*yznTyXLqX+%@T3?+h$311ILm{|Y=J zE2e?Stk+gdtD4j!JU&2QJw}qh{W3|^+bC%GUANWo`8+_@3}&QXtyBc0zh+CMuHcx_ zJ^Z4h&~$2mm+g+Na9%LA;o_l6uy?Qe%M&GK-nmpAZ*=cF8Anqh*29?Nd6vLio!?_r>`R{wN!e z#OUQt4#1gFaj?3bV{Y0%mUA#!F$}$5QoL3*8A$n`fH^HE08G}Lt^6_uTaz2|hmP?k zK6qc@v6@UEq59dv5dk%DkB9v#6T~E%YXW$2C$7sXkFS+}-0`K2AaEBd%6Ndce1pzd z>XTZMAsdLIb79J-*6y{2sEKSP&AkB(VgE7)#u`qzIWDjmN|@d<`3<6*luv5jy_PrN z@qhmqMv`(ButDiWjL#~+NnFkmXUkCKv6WqMlm6*=Hcp-1K{g(Ix9VcY{sSGJ(DZ&o zajQe1WE!;b-+jUkzy~YkM^s=_4GG^l=uK1zHcG!1(MgfbN<8y`AZQZIzp|W7C>=(SA zOis6|dQ4Dh@+dq?JsFj}!o}$8i|+%2XZ@%=+_!?z_11#9?VVd&}s3 z5aQgcK~!-V*+#VqyjP)=CF`wEJ&suVFQ31LHa!Nz4af3i+p8%kCSPo1jGVE?kZO!g-gAsb8*so&C#8Z#?T%Bu|Y#d4WJt z9DvMK#Iw0>DNyHaGA9CR2&v}oL-EW}%DF{Fiifj^EJR%Zaujuv$7M83{{UT^Zm(g% zo8tG`dB?}EWSam4MOxsuX{4!k6;zZp`B#ql=uTCL+s$cX4#Kp&GP^gr+ghX=DpmTn4;c6kryVqvCB<>kHlfJdkaLpWOiUK#@MLzkd8@;E z$x(+=XVS_y}+{y1{R?)=zZ@emQXP-0yQQlhO3KT zF!pk}IzGZb%_&O}+n=-4#zoPXN~VA~l&*RS{E7(S#X3P@htM=bDG2vonIcu)LvTQmPKRjFmzuMcb z#(@;0k|bP+vl{$-^yV;W1Jar`TUjd!^pBLFVripzlvY?G_J1pduh!PKs|<4$S5jpu z?G^d5JDN#HP`Gp)*6^ttnma|!MzpVpss8Cvdm9#&cT_h^yb!w0T;^u|;w}s!h8$mlpgK7$(;L7ytIN%&=>%$hr+%nBW z=WcwH*>jZ?t|YukLh=XHCXM^m2WOzu=f1&M=S#a@uAf$_6%so9F!zTZZ0ciR7m=I0 z_i(KNxU7HKbKa-y*XLhWub7ymW4kokWwn<$E1LRm5I3+y=sYG%0*^2Mk-rH-zKdWk z_6bveGLS?}IH96BT=Jl9>+NaMk10JtpmKg+y9Jk+jK8o^?_Vn9kk&E%4GH;lhVd1j zLukik-A>Am78Yf&$0Fz=vmxU4V_&8QwXmL)WDs6}f`xx|6RydXr4!wlQ*SN|i&M(N z(urlGa%2+suv6_|pAFpD#-S^ZlgJhOeQlS6P2+SUzVx&;30BtZ7<2ss+ZZ>y^_XE zcboolB2;BMAUvYk-rt%0fst~u02a^!PUbpK&_seu@EoX3{H_-(+d?m4>(1*o@#m0a zT_fbCmG@mO56=}-UYv!1rr+Z|x1GPsC$}$$Y^(TtGpvhFHk=z-6Lgc$#O~Okq!<+D zsUw)m;m(AD7r>_Frv{^r)WGM7CT7M;!S+Ap|5{jYB1D_1huQBMEa!%Zp}izehLM4} zg#*`<)en_Z?;N`o0cs~sFGnjLRX=mbV=&+~xSk5Qu|MN%*uy}?Sf>5m>|W-|5$0UFs>frTpmb9|3LV|F)1B?w z;E@@lRZG$~iH%EcB+m7r0|eG81M$+L+7fY7*vk#%A-4(1iW0<_`?@8Z?KrXaDrRj4wk#Y%X6celTm3r4 ztQ}c-560@x!iUUu9(X_*`aI0C;;i~~1i-uN(=%8H5isBT%Sl+7moZSG!b@<)&zBzC zU!70SR_c-ze+}6kX`xQ;?zmKVs+o}Uig}jkoZa(TrjSa-mKmp;1w&0kw`oOnd9*5f z<6|?t9iV%U5{h6oNk+$dL~Eq$*Tfl!yp}a^!?)ExsdQA=-$`Ze+DHSUNO+PYcFmj4 zyhiFXX8UWg%EqD(XPC}!Ji*7MM$v&!F{9?*(L#z+3<2XFi>lkMT?rE$5Uhy3LG*4! zkqTXe1&eeI!z_WetP5k5tk1WgKW-Mn0((_79top?1;HisJOKr6Bi4srNBnl~FyVg! z{Wme%c4F zu`lBF20RgijewZGE`h(wT4&&OB-Kx!)8ojyAi%}U7tJuP1)@1sOy9;{4NF!5VXVm_ zn4N4W6CV5}vyXVk?x=b5yyk-J<9!qB&m@8oGX_KGx<`e#8@2#G!u?eLqn76YmJ(k( zt(&??mHX*+9g?kkl|`nPyi_7ol`gr+w%yOTXL>{8N+`6O=^yGVbjzZ>8*9gv3^wN^ zb6cY0OS_LR@Q!W=+k5PZY*1!#0vj!iq06T~UcdsM32;ZOA#kD$g*gQSzsU7$YPTt) z4x^1=vqsEEJ*xQ4R}R~<+Cdb0%6&&^1)Nl_(q0NnmqJsF%!Q<5gN8tYR^8Cn;I=x% z;*CMbdb7=ZL0sP;<5(Cg+wVPFmvsWF(D2D^*7{Zb>fp_~z`_d}arDAY9sUL&^#6ayOm3t8GFk zbt8nqqrUVgjV)I9oA_rCfX%i(5OA#Wb1#(w-#A*+d{P3Q_qDu5z8OdzkjeZ7!y z@r*+F{NH3sa=ZruzPJ!&$v)atPOFvF?d;Zl#c7YPl@DY+lr)?w2W(t*Hy?MJom$u? z!CeUd5j3=MMEJjvAN<5v)t+p~Hww-4Xl$dj51IZ%3JTvjG`8htq20MUpHLOB8S(X^ z!|iR8yXCn}1h&~wdW(AddmOF0l*`E9xhvca#zF&pV~2JklbaT++2U4I&jzaYK357| zuWVhRM?nS7tDlA+=E2ZazAznb6STwnR!oiiZ_y%AAXbF~F!+!jv&0^${V!LJtW+h* zkg|4OGNW7(8=($n^TgNriWe33_mU8Dg6Bo1KT*Up3N;L&?;ja)Ao?`0R-Vb*#oEF< zOXV{p^0DOpp?c5!Lr~jqA)p$lQt2ECtrgAeSyx7ZDX?|_*|xNN%@WeZ>Xnxlfiu0t z{D|>o+ZZq2m!I9^UO&cKH}gKIfKu#|kzLXt@52Yj7tCM6zM=U_xrODF1HI34pEz)u zq^ipO6-gqT+E8!SUE~Q63K7A6>D*04b;*(P0A*A$HJ413+8f_I3Hng}3#5ByRiH{R z7F}>-h&YUJQEIB+@*~ADEA~Bm1@&}bx0aNwQp)sYSPz{?ptwl9CCCEqq! z#sGiv2u0pE>u>TFqZ7f_M;?hc|3QIsXH5>rYn0-Ka{qoN-wk@`0{o{n3Mw5&&zHaY z-3YMX>!H_%He+_P>qZTj)V7@%jceFMBI_YWS4T*C+%H+IEr)6$JzQh8f)4dSg^X}m z$xy(|WC4h|1RcBR^H(6oUlpTonB|O=@JV+8RUE7o49+gKjVMqrsS%L0@^mA@?X8Lt zcOn*Uthccsl5Z5o1T}SOtOl)|0;VCr)mpSQr|}s$FLfkpy*RW~$4Ig>)Y8Ce-Ct~X z4kTPf%KmFza)#XPCi5yFOJ|MHExNhE#KZs!E*3A-F!o?yv&m7Zh``n+B;URxvu-yR zZ?l_S1yRxR zbG|X+o43TchB82M*ltlr389X%?ZbOnXKbGHXEqMTaq5PH1&Bw$T(7@e`dd>n6(hw+ zoj_TpWYnrSU)m2O7loEQkIkTgT?FD(PGR=Fh_aZ)<2xIjBkybje&bwNA*pqPTzN~E z`$#0noXmsOnc5GjlmlPOyUgNTI`JJhz}K@5p)V9D;!6zj@pk%W^}-Ka-{TKiY0Ed- zE#`Mq<*I7eoCQwLmUFWK^<4B-F{GtF+XBWl_1@%fTf|#{WM?ANc&yfar2oL0A_x!q zvVBYIv;vS((cIS?pI{j;@lxhTwFRJ8dONhwo@p68Hdq9O!XB7m??dj6nuuN7X@Ikb z^ko8U91GiS9nbO{YREah+u8tF(h^oR>sZMkIyC7bD`x(eV0=YX(r*{P#9KaB+s)~9 zF8;b65ZAU^?hES{A)G83eH4sw5La=yX|naGR21IrI$z26O;x4Q)tQLEg(1?>79Sa| z?WDn|%%vOJIz0v35+0ruX7CjyZ$!qZ6NvpL5U)hj_X3mvH!3gWD7N7?loDAtw44{> z8W`rVF#ykwd9%e6LqxMfAT{yOLctvztGmSV>s}uKX9*p2%{ej`aTw>U@)kD>Ty&(! zq|uagH9AE>p*!C3Xv^3B-f4W@4GIwm_&mo_D&h-aIBG2M2UA+WG`M^ z&1RA%ugz7OV<$IEu8G zoO!xJD!*od`)Sl3*?Lo^b+=)FUMTsGupzn74qgNk3FSBmPce7jCJn51r-Il3p2%`E zTqEYwH9sF{jOh>F_N`wrp90)xmT3zP!n=}uCY)TyLUK>4~)W$y84NQ5| z1$vL}lu##jn;UP9!>>biEuL#Nn19mnbRSjgFdxF*n`*eGk!-fuQ-988>pbZu%7D92 zJ+rWsi%2O;WB6lt@j$flX8aE*B<9hepVnx@!@MDUW=Mb1WHBMDbK04{;Xkw_Ra`+p z@~2J%p6Naj*tdrWd5D`r%rppwuhIvcrCWXXbUSkls6sfx4ZinL6%X)zcQ~@wV)e)i z;(e1QNahL9{r9#l;GF`qRAM$HwCR7oN*YdWZ%a}|@VRC>9vZzP-*e>NH(`ptmQ zgb=Wj$1LVYun|_nA7A5}JnkNETE{sIqr0~*3^OKTSKeJ`k-%d&Nqut4YgE3|=!>K% ziF;k~F&Ji^i`0F4y`+YiGg=_^L;gqFU(`KsOm1E7Q{4?xX zuHF45U;rME{w_PjA_dzpR?~NFV(D@~;zyE@o@|GH$b;90(G#7#YSVO(5M#JgPpXXV zox6DDHS&yuxI(1SZEn}wd#m@^!oq$@z%RXg1o#3Tn%$m$Tn8-I%j7q#UL_dq&YOZ6 zvfKI+-N90db2FV=P_P`$WHV10qi^WIb}?k0Fj}TsHXoZbowK-Ycj{8F5Nmo~ct!Vf zIT_xJ^x>im9*L*7YLPI2J)OdFjfx8|@giv}ZJWJGl?fiucw6wT-9SR~_d5EeNTv@l zo{JqPZbfW0%NH-4WWI-DnPfB=b_jxsXW1)S(t6+T+DYn)p!$5@;7EW3oj8liZLa@_ zVMgDG?S;3lFx7Dwx>@n9`H;O-LsNqLP~Hj`2^s%gzWPm#z8o-!o|W;Q@Wlno(fPVW z7%=lJmP(_aCnZATpe!b-5?z-&0N2i+QR*mA)nKR^S+%f1=sqoHFohJ@y%{`26EF2M zFDm>J*wHVJ@Kg>J5Ri#s!4W!%7cit}&+XagA~4)4$FN`mZpsf9(L5;IHK$Xqc^Fbm zuDFXb#1{j+ieQs%yI1{jZ@V>2Jh&5c!PC*lF?4Z~G`gXO?07)N-MB;a_E#f4(oOH) zu9}6OF4eylJ$;gHmeU6dBuaA*DkZVay zVKT`yA!D^BKXBOLD@(QTl2C2~6h#BEH^2S5+CX;PJ$=pjE-*l0hrCx~mgf#~t6w%_ zTE2g%WDa50TQ7Q*Zkm6AX9Ig_T?FT`#miL|q$QCO{{RG^G#;ZIP1$90ZJ4<>8C6JL z4^UW>lyy=L=D7O{DbFU@2);BFgP-nZ4zdH*#Hl$_n)*~R_oP`3OSqlwG=L@Nr6-Yug1XqvYT6Rih;0AkxB}lY;~`i=CZt@!AL zSI}5=@($mj$FHmfQ7%Rvu`m)dN~neuXUd@WPsntoK5D5-Z|=P>FwtDRyPP$jtm&?y!_{$B{=1$1W}T~ zcN7ke?{C4q&H3>Hu=c~2UdfGYiT37}KbA1eM12z8URJ&`086O6E6QM^c5jU7n;ecCxA2cf(v&!2=*YNjSQBWO372AD$&C%1skMU7<+YJDb)cugq^Jnc?3y?BC zO}aP*{GUy(Rd5=Rql!tCK<>h1zekG{rlJsm-e5S;So zOdm-KQP+HGPZ=APM@;FES@?){`?mJ$ox&Efo0e(H^z06OpDNZl+wYnn*jKE!Sxs$_ z_MGOey74LrneZ+Rq_L8Y%(sxTLjbJeiL!dnfI9*AvM!xUA{%;8UkUwSuJ<%zFIEdq zKp~A5sHF1JMFztT|5rWTj<+pEYI9^A+BWs|ld3`g zfrfm1yq&TEl)Di8@!u-_ODuh+87NjebIo6FqN&8oEJ+} z7cKMsO!&DjvrSf(myz1gSh!IXwn*e7@;e~~6eP^=qGBz5{vbQsF@z;u=vIGzN+9=y z3qi2?tB}?6xAo{NrqFlD zdOoP|C-+cx+$eecWSq=VX+e-#9^fnMD;674m-p(j&VDGE(Vlm*?rcWuj1D{&j!u>I zvZ++pE$DOJoA~-5K`2n-V5BYgZa#`3Ro+5HlX$SX1gD6voSfZZv{(KbhE&!DRHY*# zo9-OgDlGwed+65ye)f3Poi!wO}QTt891HczvU6+j!mD~bB-zNkfC`G+C!CZU#EjaqpRf_6p zRH#{)#xfzyDkUqjaio>0eQ4$=((_IfU#D@HqffGENoJMo+2=S1{{BM{5*}~u44dKb z*yO@oS2?bGx*w-{g<%MT&MT%3%nsGUju~YbW$$mBY)~k6b#3Vm>(7Q`Goc;ILiO;R zg^QPvCu@&Nzm^l4X(fR!CJr(K!ELEU8(9Wpz0f}b6i8a6Mq4-x>{(ta)Nv#>=PP_A z8pQbQIKoT_>xM?Da-EA48GyB~5@e5Zq}Jj(6-tHnG_)Dy6Li9k`KA7-{7fF5lW8wJ zqBT_pmvAzb7-}BzvGP%sg#>S+8>wKl>Yb>hirLU%iGjcJ>(Y8LqPU+E- z5vlm-iLZzdJxi767$FhHQTU3 zNl84?%+H{5er7A0te*rVk&iT{@i@)~L51{|{DHWR00Z$v&;$6h-2UYxev6z6{<38) z_)s-5p^Om18`qpBH3=51A#r!AQlRq-(5!-MQ`lhzo(j&a6{>Z;R$+S&F2hWggKfIL zR>%+ZX8%dfp7Y>(5Idj~!w-#v!`>`}fNO4v|4pX_eR#x)?_NG~>10MKv4A$YL`Xv< z0}`CSH-E>xyK}29N}*fuNe{OYJGmR3llTDOIg5N}+CMt6H@~)r<47 z01^H`vR_Lek9dsZSy5;fM4@c3aDZ>btF;(4X<*BzNADZ%5EOh-`K>`A7w}$=bz!%R z#HZ!HxB@QJ#A$-BjKmP@=+=C?+pnje2l1kD6o>nzNnOCG0J@!v`7FbRvJI;Wg{eOyj6qEG9{u$Dm zPUt0z76m0RS@0H#UP%e^A$HA=1SlEDihgkUj@n`VbsxC>bqxOvXwVrJSh<^`q>*`+ zq<^>?*=bbodzToJd}VsYpbiIyuBI#z}bO-(Hq`mNyb&! zUT))1sJ$}pH%)E)*OVqF&H%74B>)Q!1t)&0C%=6FI@vsrslVTbvEo>}(kowb6eRNT zmG2oc<=PmnO8+sT_A;oLGq$nsYhgfEdx0-2Q`u~u4ieMyeu+JvR>G8vddi6Ed^+^w zT2SK7p8WnNjei*Yb_pYC=O?n9qJe!}fVlZ}zCPYh%v-u+I6gsCFtugTRfNZwH` zD}vJyq2E|Q8HLJqOjI0>FvnTZ1<;l;+7@NE4!TLCnP-idrZMNv;RsT+C?NgczNM%+ zF>p#6v{i%tT2*y0Bok=?OYva7{7rpeJV?P!VP5Ot?Z)w(_Q zZD#MVRy}{-+}3scTV;z%Ue?-b_n@81#4U6U#DCB5S?ibj3td_!3Hpv$NvimF5yys^q>#>;_^r*iCQUgx>T3dH3lwlqvh(1VsoxfOf|oio-SS2 zteULHOku)Hy|BUobiI^_<4v&qM_=?Sc*Zq!yOMK^a7U{9iLI#sM!L9|&p^GrD*Gs1 zX^R_syoHm7N`ku6w#%Q1>ekn}%JK8~gR$c|+Xr&%NYqi(Sh&p%z<+nTf3dOAKyj8% znS@-Zao3$6bokPa5M@?wv!=00sY|cRR3HT}3N%zI5E(u@?i&Nw=f#8Mei+G?|NIsT zGxjQ73;!zB4Wz=?IVRukyxzw-$CTSn_W~F*U{8#)U@<#AaJP`W@WMRw?H?8f>?6Z4 z*qPe&$0ISsDN_!IBvr~8iFTooUKi*=k9l=HcKA8ASFCOm@%)GJC}_I4#GtGtSiG%QGAJ1pvr^hKkdOi;bwUeoDa+Bx zAwmvdloE!%YM*TSSvn8pLFJ`hT>|@ZeZG#)9_Omd0-4tsw7^Ngbpy9sOvG?A3a<^! znBhgLWea+~ap%0aK9O^qo4MXX|T$-WA>kd|IA8k!8h0KC@VO%q*-9s0pS*~U20SJ+D z5bg<1{=Xi|SVLT2!1na(zY^Kwf0oP8_;$*`58XdT`7r@HCF$wU)XN$EcB?-ZJA;e+ zx+7!FnXtVIW92E;pU1La1iaREk*!3ub9T_@w0k{AcnOQ9#Wip#4R#A+0j6L=(mOP% zs$PG$%el=$ISaeLnHt*DIxs!!eF=e2mv1n7igw^*VnBaySD8@+(q8R2wXNv5;PpKI z*+PP~Kryn+pELkL-c8BjKe0ysM8U>5R>?I2wl44}Az1;v19z3>111d)NFy`lG>Y!< zdmZQZQZI>+|7t}osJM<-92Y@Az?`m=!q*1QavRc{1n1wrw3jps?&(~eQ_x|L7Qst< zqn6sx1>L6o-0~2TLcdO`$xn)Odrp?zbCCPHB#Ge*UElm+Jn7H}RA5%>dr*@~w+_M8 z93{}8l_guG9I-h*YlXPwdL~0Bu?VW3aKDj=RSR8>SB*DaR`UnHfp&_~A6J8R)(EaS zc&umqe!$pFTJy4`)?fZstmXeYXaxryTfZdJw+2sNog%DATj&a+4<7B?X8FP~$fg=( zEirFhRQ~+5Zlqx3=z1Uqy6&EW?r8?L=e)(m0-Z+vTkuy!!bUGu+@(6?68DC{*96PD zbdun9iF@e&YP`**sBajo8?ac+6mvDwZm*{j!i?8OREMw`bX>f)heD+hB310zkL_x5 zd}ru+k?xinh7nD3E)x*wGDo&D1&Re~M@{34oW2wH;^C0;n;K+7R^3g+w9rirk9~)Q zrsXZSW414;i3KhFXkPf4$q+&(PGSg&ma(WYqeBrg*&Lri_LUV4ZFz>Zk(_yGV%^6h zQ+z?m=8P}6l6RANFdq!8kX_3D1U}y92uONdY@+_?%44eX5~k>OSyi%M*-wt6MsBeQ z2|p09_SwQuTJjm&S zcb586m)6B*l6mv5`gGPVT^ju+oV)9cj_C$?$-mxD^B%Lo&Ij+@uk;9c?QS`QP@WvQ z>JpKg&~W^oW_eTJ2zomli@n`UT+ZLVBbSbD$Pq7Uq0|jjY$qjYP1-ffKYBA;y==eZ zDczJ^n6z^ybh&z3nsMDN54GxCYt7rBM|$DPj%$G)i%=&H3cxVH2aU5%6#_c3o%_S^ zC;aTessb*Ggs<^&Dt2OelQ16Ihc_P)1d$6N0vk?U+}Xoks_zDeNeq&u!cMhT>$EVJ z_h5GEGTp<1P?i9RfGUG_>i=XyPX~@HqtF)TVAYCY1k-~i(np32ITAw))`+r5IYSM` ztuZhi6sc?6+1o0wRC$Kh&mZp~q2>oEB1d=@1gN*$pA}CbzIruAxJ#Ku)tzVJAvHsN z$ZnC1*e^M-KQ@7l$*Vd~zugX>Bvq750IYz|0vFXn;*<;fl#e>QTu+R(`0`Io3_KJz zU(BQ}jMU0O@iWdT=nl!qEDmE4A5^x-D%5Db!8d1%l%RH35%vB<>eNCWKUHA@r4mW0 z4_A)7_BKgq7ts^AVniC8cE`AwoL2e>q1OnjCg4eCKy||eLG&+)sggJbht`&7cD9=D z+zc2ot*{yAx^L0dFDBP8zV0(ztn_=2?N7sT4Mme2I&OTfVsCb~6vT)TzTPg~OCz&D z$>AJcx&oxEaw1I`VH1Dm_`qneZ^6hjmpL>+(C9K|7I9+I)YdH<97%j-PAvB{o!Cp@ zv8Z#;=e?GQqVnUV$|P9H5s%*33?UE``p%0grPkC2WudM*!EGc6j9g)MVgv1CurRIn z{>ZP8T(26xLWig`L8t8+-j9n@^7J8o2wY}#LP}4%T}BOv5X&vHgnaj8nk<=9R9K95wO!T%12xr1 zVYITk%EN;6%XsC&XOm7f*h|DFm~MtXCxW^sL5&FBIO|37%Kzd=gPqJO3Zmx~FYaeS z-}~PMc>a0tUj0{0e%cqof5S0&->Fs~&|MOwvX?K@X@D@OGjIgGQ`k8^#zT0ysHI0O z0l7=Ubjfl|89IzFYM1x2CdxnQtIvw(>_fUhXnmG__84*TIyP*`5I~j~Scy0;6XjSs zDE1a(8$6*0WZ!i--0TN6DV*Ck%mdeISRIi4YHu!pX`wGzatSr*hCGI9GeE{!kFR%V z+b8VRX)8Tir;;|kl5fzw^u&8O(_WTFn@`mzcS577$<7Pzo`gB?+~kDpkBCr?xE}b% z)|jekI#~kUe`uPiu7xjM8-P4>o0CL{;V+j00kvQV28jvkY=(Wh+(i^*qH0dR_7WQD zONqqA*IkaP1-d0Y)=-#){x?!#MTevh686fdO4+qLR>}zPgKWmjk(bta z>^{op-EF8UNN8}b$YNY=!F4&`rCVahM&BCJlWu; zuDVRea$_{~LZrm1$_?7ORHG9#AQ3PPdrW!477+0J0y4bjqD~t2vsPwmNx@vP=fZ>UiUzV&Ug^_t(ka(v4Fb#NwUx=zT)VL3y>YN$HHnE zIfK4$5_Otf&w?9vja(LE)VQLm!9hmzS9hvIUjts!QSine6YU%O*%Ou&Qax&2Qk5FA zwZNK@D~iM>A`d&GuTls6V1UzKg9sAE`s3@eYo(4bDf#KU#)_Ih2TvmXYlvp%g+~KwIj^}}DN*sZPJM<(!_?vCWYHuIeb0h68eCjE zrj86IuLr}9jpk&&2+LkY`OYh0CEIG{tDYzKmBa{iT%0!USaL|{)Clm#$p+GuaTt!- zu!cX)Tkd?p?u-=O=&KSKW=q1C8iY~+4!VzGmi?Fesv-D&jrP?>`;+&}rpf$MzGQ%2 zZ`{(ABKTNdIJQLq0-!*1G`{_+Sh?(c0|%DWx5qESxXOR~dWII;xKnOJiWa&xWvx8m z+1OCD(GYCMbsHad{xcUfDqJz6A)Y&>RgVHXI~B*?%7*`qztX)zjf^vai@)84F$CZ|0w`o{TpSC1vb$* zSr6>+p0GXA3ekwJBdXbNmP?>YcYbBxyd0Zgy85>H3yx!T4o1T*sb8!=y_j({>EAb6$^@UcSBHGN4q#YyTr>JEFU&5#GnRuyty7E>av zyl;}rBZmTvsqk)zh?B>%VF!ke-uQ9ZLRw>D{O$;b4EMvLZ(*ktSYyJq=rR8V#xtE7 zVz6Co=}ntH)SDnUWKAS6`b!m`%BziB?6E1S>z{0Zh`K!Jc!X<8!`QZyKSiNQQm{Vm z;NYtcKKt3f`<8-;rr%f<2w3yyF|`agCx+ZqO!PrnG3C&RKVV9$iD~)~b-ie{E~PTs zarp>P-}+OD=6z`*6_%gykc+H?sM+-f^i8LO=+|)Ty6)^x)t+Rxv=S4(mL!NhP<865 zOAd$%0Q_}4zwuF7*bS>GDJgIfa<+8sm|pH*KV9Gin=Xi88!a|ZCisX&`LH1vZ+#1J zNKV4??%IO|sf3pADHW%Jnk6GL>TX@GFp8scKvDPsKjHCsx`%Usep0(d)bYzXU_Jw8 z1x{63+^hws2tI+YMY{E`J(Xp`f23YJb+{`l)O{9X#yOu59xq+zG6T;yCpI`S$i+mY zl@?u^?%rbtwSfBh=5M=j&`Fd77_luWJ$GvHK;F;rp5@7S%E}XE11=`dtSCap-0cPf zi{n1NIM7&p00|HvH3!{KTH=*jOE5{fHK6r_L4CG`pTVdfACAf$@Hq^~w@a|wKiz>` z1(x^W`%2%aFrM>)icyRNoO4r7Ur*UpE5id_WxiE=reB~6rH=p3BV6{(sw2nt)Dafj zwPy-4<{h<;W69?s&L0XOaFG4{Ohg@qOWXvEQAKpQ+(L-yWMnE~=E5D#2`g!&(w#ZXQ-V=+9vo~i zio?(2c4u1OIjMI{xv8WxKpO)9_r&mP1S%n8ni3ISB@81ohGp=#sK6jNscPIvqs>?Gs$es+mr! zDXF*GOp*WK3|VI|PosWN_dvHR$!+FrVisb`Rz|iJCZCQHR7q{J%s#^}AI2UFZ_YQ0 zeJn>5hu||*YjwOU8R)o3`p`9gsnX448cj`TM)p*2)mMdqY+wr5QO105ni-b+)zSJv zith2Mxt?}dzIxcTLj%eU)x8x)8y;3js2_4>afIlcB1)3r=EfLBD+x!fatg(&2DGN} zj+?vz>$)Far7r)Gsker`Ta^Vnja#iIqvmJK8+v4OOfTZS(8Co<58}a;4^T?^-<9un{obfpuImPOcS{gS-VHsqt0QSCq z_@n)@<$_={ZHToxUE==b++_XE2FD1r>+(&HlMr!jzMb_4MG>t`EQ1Zfw?-cVgb+?C z9kJNQB?NF%gSG{_AwP7hs~Qkw>Ye?67n2=MuoH!9rP@#7ibpA}Sd%!c#f8y2{{k(PouSgjrN0??Dl(_@l*+vmNrk9u zTh4}7HF!coqI2oM)!0k`#?g;Iz2Mu#wWbfdiyD^8=y+U(rE5i^NCStTLf`TeN=`{& ztAH%;@yXU{(F2+oSem0S&_t!&NFGSP;{>ERam^7~*xDf;a%pc@BReiD7_%{fn1u#C z8bgCmN?1!vEUPm^)5-rtP84RTem!d{8hBRtsQ)ue(9J=>?N$sESll^zYkxAg+d?sNeL_brxY&}40 zs|gsx=CAlW)cF2auTnDOYiw~XH=zd1V|l>W4OUBJW`VsS%q4P=#Sut}m_cpw9ecA_ zukt^u307hKHE(#!mGnTz()zW;<&dk@!`{d(>*lN%0V9G78)=p_p@spRz`>7Tcyh>>#$puU0maOg^_M;Tf3~-Brz&;((g;(`A*Ri&3c;eeNXA z;L{j#M-zD%ZM_k4Ih8l5wB7CvFU3kwzT3l#Av`%V0PcfyA&yAoxu@~RvFl@VRMx-#?ra3R$8OLB6xw0Af1k;%A_ztKyQq1wP?Ia+ zx!y^+3x(xiN>fmj-j=;m^Mi_Yv;T)>KjZvZA z(#!v_Cbn{G{nPx{}F>{>1^K&uXx?R6K=yH%YS`d*?msvd^)JUHv_ftA9z8 zWF+jwY;jkVM0)-yK$9towg1BP$9CgfGjxQFWi4dSyy-f~o)2c;Vt#w!`S)a?OhF&E z3NN{MI>excmtcnIj8mh+W`tj+oW2u(3N_By@y^f#oZ3%F)s>DiaHj>n<}5|JR%6Cijav+a#L_6s@Ro=F~j zVJ4Ja}f?uDQpW>rlIz5mzvkFyD-K5=>q-v*+-`weUo)Vc*mb$UnxEnj50!rv5PV`+euTLEi&`T09ho8KKlJ|)Y_R?O z_J%(l$5T+F?1R%3{*M4>HqZ3Qr^6EC&Ks;P++aBqt&YT;kEgAURn zNzW$L6Y9=X`x$QWT-2IO$k&Pu`}^7e$S(N0ByXN7cX!_2OS-S7B$6+aW(K?_LAL|S>%uW&t9@C zU#93LEbZ!?HB6-Zu!7U5oW2^G{da)*+H{0z-y-+%3gqN;uZ8q%qikNxzR^H%81`kC z%MPz+z1!B3;{9uoiK^}!e&PU&br`742m+PZD$h$`ypHlxm5vM?&q3fhqojoPO4TaX zvTZs5Nlg@oNfk&XL&^?p>d4>oOcFSeW^oU;km8E(#y=BMZ>kF(v|FVCtK!l2n4jjC zqcvPJN*SNl^d@t65#z?!v8P8C%okHbP*0+OE!K_s*;T8 zWaODpvV%~55G=}_fr`%$WVrFjPIisI=jZz;wU}2Yum1tgWJ4xen4#uGq8p-HGaiCj z)i6NUuBB{wv9F9fQ=5@U?uL?~Rpj?z8LH)E8&Sptv-mVm7_ zh#w?ry(KvOv7c`AorB;{shF#=j)w$$xNhwK(!S}RcRPZ3OHMM8X_sY+VP8Oddmns! ze6RE-o^?0#0OhUo`2L+f9M>gCnv@dJpf!y*NGL6)sX~5+%Sa{K`ReW7FFS0G-+$1S zz*$2WS5@NXQ6Qy|HbJ;aJz3n4o$0`Bn${kq%Qn@JaWbagc4BgV3WpKz;Ld=5E>Fss z<~^buR43XP1eZnO>67ZMR!eHj5O>x`mF3#*R7D~VFtLUl0q2qwOX=2#%vpSEAeOrk zZLdlKOd!f-5|o}y)WTlQA<_xNBJ=pgV|5qrO>f(ihZ>I?f0a0c8fci#8vguvJhQNj z_8kMxJj>4H?<7GEcHlvEei=knKf;_;;f>*{%uM1bbqvkU!7)D+VzbgwJhEhUBLTba zVmhWWvd9#xWQU3*86B?N$~j~%*CAp$4^Z_@C%cx_#8D(MSxBQKBUk$rtKuVh!NSa2 z>Q@FBMDs}zaEn^6Hg#!nwz$XXvMY2|?hc7MAos^NG`=g%{#2aE(s%!Pkvq(zZ49q~cNr4*m`+Ym-L%PFXVR<@7j%^;h+6mI9D ze7l}hIS#=eua|P?DIh>rP?=Dsyui0+(kkioXY}Il#xZ!jlW-q-z*iE7<&ZHPk`7t0 zlwtk{{|zxF8pnmt%sGAgG8iyF(K|fzIK2YgXtR3oZApKa$Qe0I({j^f49z=G@Imzn zw0;;0cFKvVaQpLlT|sEwdV;;Zc=f1w{_S6Xu5SF5^*|Nes!;X!obpBObGM|lo z@$@;;VGgjz=v-Z9WC}1P3-lP5xEU#{b>1v#jCSs;*8AoPP3!AL9ZuR`VU?O9LSR;4 zZprC$V2oQ}Bo*oFv0e-Y#T&RuxD1uz6>QPE0=e}o`&|j80}Z6Tfxh`S3`lX5H@U*< zy_3uBg*{gfZuu5Ldolx$;Rg=*Acb?9ePe%u$tDn(Ay1l)_o|dP)H)E5s)#|{A!_OP zU7-A0&d99sbm4~ylI;^03#7+TB%beP=Mzylat&^Z(d%w80@1T!r4=bGA2tCz*--#a zZa!r5c)tUA&50qrCGDhQE&dcqYbG!JLQlmW?)&)he#z9{@6;7E!ctsGR~_c>_69l{ z`GXIf<{{P#3$U!g*MgR%o?R^3z7wOxh8Mq-YQkq`3$AUWK04v|hyE&|?h_}`F)fgP zdFlL_4?*2BPTry+k?YI^sC151Q-gUgsE22!AtUTyTY4Nz7m3%lgYlz#j}tL;l9Ww1 zSLczFCZT5tuJiYgj_`vxO_&a>-jvdvnhX=6a5+ales|}3jEVU2;@#ug&!(gIc8ohW^&^W__Y{?AO zyVuOp(ZkHEGmGa3fDRyvf2Rgyh7XI>vQd~QF{|^(AE=SaQ(G`YNLH<0!~~G^L)wk} z(AhRx_0d{P@W%1a`~hJE3N(3m@-{}{7)DQSER@v9H}u0+)7MD`?#4ezg{mXgm8(&^ zLl{;%E2w{rFI%oZ%@rXAA?{%O3Bm{psW&gb?8QWcqRdESU=G3n$M@cc50@79+#1lW zveY2hSQ)&fZ`P5!Ab1-+Ho46qT;*wG}ZZ-Ez~Yrr#xQXMC2QVrEa6TV)Pq(-mp94_jyhNs`MV-Lj@(C^$PR4KEx+T}(RFys^9o zz6VESJ$G|Z4%@&5KYRyTV2;NocUL#Cs;j4OG8G}+L$iTMMls~qtJ_V{st(|hTngz+ z^&*!!UOi2bjsym9g||xuVOS3Ch4*jv%I?*e!$-S4vxzuChTE-;Kl2RboaZPAhaaen zC(bn{KEj~AAQS8%k3!+XvlEP$WO{3WQy+C0T{;k=i%aLLCZjJ!jC6OO!zASvF=C^) zmF^$+DK&Ot{4#%cz1K5{tX3yJI)FfGsel9Nn`PniwIb9DZ*;qrxOa&#hR{C7-kbm^ znbGhr=D0VFeb&9Fq4%B)tv+#AuQ3>x56?XTJr^v3{RB$r*6Gu_y|n(RgmCalLf%3j zH$hVzXi0bL6UfO9r6(JDMPy(8C8&-dOq&F0IKpEN&N4(%z&ub`!{YOlQl_`?KB}6V zUgPyi=CrRM1UhkBHn8z9hZLDV^1VxnL@lN8qv{v2c1D1BMIql^Cl=$&3jzGV7H(UX z=`?a*%KlyZ7~Aq^l8cf0rNz@7sjH?3V?ylh2Kf>|Xg%80%SCNs<}=&B_RsrU9hB^n zaf?nHu&BkGqBy2LnT6N)R};CdtWxgc00+?)=8O&LU-|mEJHM!tm+^8s`XEj^xp`yc z#IOZfAu>pEoU;p}p#PSh`(BHE88h#|HtxF5NCM09{CF?_(>&zOG6=Qza%s*4f6l!R z&Vi+jr2)+El^CjYQcXFcCXXTpg`7U<%pn@Q2{|-~sx1&s>?}C6i>Ao#u;X|`y-Mg? z8J1dKR9FCDJQ$D;mwqOeBMvvDr^8QjZiQ>*5+}&tP#)Fo(=(H%>+$g{3j4L}u`6E* zF(JaxVw=W7b+##xd;!;w&lbU$uIk;^BO8el>Klyy1U=Geg`r$%5SRzo-JIfjK`5U* zcM|10#o!o;r+sPr;%7m6aT2%icIjoVu%0lv=u_LdbmC8438uCMO%o5IKDMJ7T+)_t ztJ;(eR)#ioGo`x8`5Zy_11p0kz!yVI(}fg^U9;8ea`k4ZxXWR&@dv?yc9`>96MFG# zYI>W`hB|?K=!^F#s=*$wr!ADhXi0OcfooS!xN&U*8@)=36~?2FNi0~Qhbj5 zOML6km})o#=>9e3!Y*K29@Orzc{MY&WAiZgvFAeNrQPD@L{jzm{-ES0BHa629{hHQ z9*sA5h_DYwN@RCf(h7$=~(Ef6uvg7;+|cNSHA9IMch5-1l8 zbSe+(;ICS~Ii_2=T)V14NF*KPFeiL=?iV{xLW$998Vu`JG73tdxIS;4*wg zZI!#ke>yw9hGO@P44Dz0S&}k}!OpB#@fr1sU9^TWDo!YwQkNO;A|!ErJ$SAdLkLN;Yzg~=+Q4iCuL{EUOdyz4BPM8}ZY0>Tu&XZ$@~o-) z##Y*7^cG^>Bt?tmnUiB-*2eREwXzXLY;l~+D$;{du;J;9oX%}(+o)P_|JWs;_SHA5 zn>+GHM81Oz3@y4tY?-5-Xg#?NQ~zz1#Bx$4&^V!E&EeX6}cn8SJCU^5&^<$mo_U(Z8s8T8`F%%x-HQd6g;7B{>- z-etaHy@*Di+E9ERW`~Bla;B55KOg-$QT}i{d&&&8Fpo7Obn%HZ2OkR6P`b{bDK8{i zJG8BlRoxy+pF=Dbv;$h;Q)QY=Cw!^kMbQyP;VLivWEjb);rIs7*wj=!jmwyrUu0b3 zR|$e{NoOdFV)Sq-OPs@|f}8BVb+FQKgu0RHYZJc8!HMPmeKJmb z=S%ZZtBPi}!ma(a$+{_WpA1Ieyb zl;V77&Jyp=eNgE?vQBCrvclvO?A!FwplD{z1Dhd7`66$D;5l7*FGH2534r>rjd({K zH)f7Vq#{j5&__|3fZ=3=Y8nO@3JG)2^8`T9c&FeiQ(ca2XNpPdxF9_0V+;5NA`XblozSQ=r z9;IhWh_-Z)To@5-HTh9h&~P>4XWt@4$2hjRC)`jD*X-Eg&1cs|TLLI<psA?>gFVz}z(gMAfV;2lmgkIN)ieVoT1h}l zc-k(!Xtl`lQTDg46RK;L=7^20rKo3)!4P*=64wj%)FDr7u9u0_wskA%I5vZyE!Uq7L zh*N1}yanDN{jRWO{Vf35kr9P7k^woc)hw&xrjs`>EQW|pTGVKqmZ_XREm zNMl)TGn;hcTSirrj!Sk`=q$tY!0!Go2&xU)dR5P>0fj4*(NQ2u(!##67>U`{+{zm_ zc>HbQCs0m+w!@yfk7~=Tl_>+aWr{b|WS--=F}(%hg}ACN(zj-5o1U=qjujs5GLJW? z-n(_Zw>5N5!gD7}VNGZUtl$1*#wC(JoFxQR(umt4uNA@fJm!}%bMVM#Ed`P>0xI1u z`7fIGcQh4UncIvSeTDB+>V6(em@8?+E>1DK%N+mGheZ;HMN-gj z15!nkwS!OV^i>@-X!z1;5o16qY58lusN5_S$(@a&N>NUQ zn7HLvuTWLak$>@g4a1bxkf5hbhXfAVx%?~@C3vmp3ax8hwJ4b$A|JtLEUo&ab@&be zZY)G>a^k!MhJlyK{k9eJ$~XcF0NoJy5Vd6d_69*JZTOPlgJkerBVNwr7lio|PMkVs z1K!)Pw~m8xLyGLM)6AX1^k%;@n49B2H%Q)5<}^X??2(6^?<`591)19tPWZR+z^<{Q z>d?dRd%^Ie=GHJEKz(Z(5@UZb`oXj+n&}`uvT1vDmpQVXd$HRznlQ<`zz<4hYdSf) zr&MCZ)0JXdpAEfhjBrcih5Nk@z0-4=4lY`wdJ|U8vVk4XSM{@8t+`^jc8B~AT4`nV zQ#rs6+A}B|Yw>tz!zy#pL+zh$9sU<_CK-9X?fwir6?V;l*E>D$e*vjOG_WAa#Dc?! zJ|U=Gm3yb|h(wJST~AAR8;?=yZ&MKW*i}T#xw4K&-2XS4E`Gr3H#`ro*9@kov;ch{ z^1ucz{)A1oWvH_lS;gYoZe7gO`aGg@0&4Vo`9@F!Qk#>P8;>_ASfm>xRRagUw3w$A zsh6?VMk33B3I;FbdT3)6Z#hU$TJKJV)OLXz}zJq?F?&=fS}t2O3O3jK5= zSW^)<)CSv;Yr?8sE^T6)kAps-u3TjmT#s%sq$YH$4t?EDzXQx8Y0m+DXb8>H`q8%` zRjPL_u3K7O$*!2QM&Xslrrg=44EB)4T9I+JRmG2qJ`)&t7nn1Q%M0u58)%%+nYi)aU-$GVVN{r}Do8@h7N+@Dd++UXVyjDg zf&{wAyO2!vrMv-ypkGp05d)!}N8MQQb#mnxSllx&T0$Pd9x_;CNTp%-LYK4~2gXAB z@kp<30k(dBy6c7GD?3;PSqFwkB#zgHU6lZ_dt2TLz;L z+6=D^zwk#gugb(=T%!K10D{&cCtMB=e;5ZL|HdrD=kQ1)fa=nNT!Nos3z$r6qGWr5 zkjU?xM{-hu;}kKE1gP;V>ss7OmWtB|;pZ?Ktb(pllPnBtApRuPVjFC@n%)~JU1=p@WZD`aOXzFlB+1S}Dv zecM%Qp~1<%U~Y_)JqmqH>62lN9_YMTsqX>IrfCUym`RfgRKNK114syIlQ1@zw(3!! zhbvCA{fE%&jR@~+Vnji&U)m^t*m>+Y^C;}Sq)?MMoJX|ITtVH2SoB5Smt+KKso%2X zqFFl7-A|g21%ouHu*EW7H8%rIQ9B%s=|3pB;)xO2T-{SB`jjlg7XTQ^+-R$|16@WO zM){nuSD8+tB@)uD#CBf#gP)Pcf&AR1R*{?gZu~Q(V2qj$2z+s0dw7ca(J}HRREk%I z#v5-O{M2^lxQsdtlNo`?VgsNVE%;h|FW_B>o$EUr9FyNCu6>5EMdhT!$Y)I!WS_J9 z6$D6~znm|iJ63xCu=arH;=<$FH|yd&iaKrw z>0e`+HjNmkOes`MjK-?afi1syLpxK8-0={qFqXPET3a}kLe*7zRV!OO;qJ+(YFOF; zd2>g`s0ByHYf{yYEFoCdToi{TIl@p;nBl&@GlRy?RE8;}H#&gxZp%2h<-k4kiwb&C zHf5rXMj!Db$%q{kTxwVTh9{N`t4nvIix}m^C=(!C4A?B)R)v*{6q(A_*AhrrttKKb z{zf2m2UL4lUfNR1nD(&QoprfVxUzdkseU%ROz$$nv(EaXOhDHZoPqZ9tCQWkU zCmE{|y=DZQ3QsyWYtp55$bsKqjURa1sGlz{0bm-e!g}8z7KA<8bTo>O>+kYQ4 z$8W^LYOr9ZGijuzRnecnhy@wLH{GUFFNsf_6P?ilPbP?Ww&g!GtJ(Xd3iHB;1id!! zxD)yq?6)euzeLpoPs3wwF5ngj<^?ga*{8qiRTE5rxx)yJ3jtV{H{6Ykf03Bm@%or9 z6S5Hk|5${$ohzy8_=1zNh^j$tJ-|Uihe`Y)#CTJEqQrZhm5XMX!q=~?0{E(|%Motb zK(Nu;aw2;qbxPFi4Q9I;ip0b>86KViwv4@t%9)2hS1F~AO(jr0)aFI(f) zFNK^&^^uB1lmY4EApWX%*KoSa23RHK9@C}_@zcT2L~L?Sb^}T$`3Ja|cmWH_Cnk%n znrUx^x3M_LY)0 z1Kf5MmDhmz-}V5MM7;=~U!%J-txnvyXTNT9}h6~qotW!67L2T?+K zb~Acyv29jEOqks}6;q*+Um?>0rG$I4=D>jg619!vK}P7ouyygz2)x zzrmI{S{%|JK~NYI;4)C9_k5>ew`IPmt4j>5MG1Z6QTZh7h_Pcb*kwDeWuRf+q)=s);rNibaJP^+)1sr>%wnvdfX** z^Y?YeWf5LB2EAUx1jppxTTZ?0XQY;+-VPhMe+vjco0=S1Jx~)B#>Ot`VN51L3{3GA z6HZ|66A*M#qwj{=KS!RF+zjy9d>Ye|N26Jz28=z|6pl9?p1#`3I0;nG~QuA!>iZVi}a-}HApL0I*pr#)+ zh~SGQBdrw(O}$2Hc~lRy+f_F<1`MA~FqWqWOZ85cw&uC1Sv*?_Pvu4E?i3Dx7H4vx zl>-+K;)JcL3yNS?cyL5=bc8&(Kf>peOr;PL76x{ej zrI2RDYUenQ{{LpI_r9Bx&v}@*39RGt_q=q9R;&#;#q)p$zvvys*}Sb9`0x>NSh)4e z;=}cihCin@^Nn9DdiM2Rw0+!zANBmE?2Sc+pnw-hK^}q$- z=m)CC@ar>FlA1Z|WN@_-&+qi4)|~)Ap5ei8m&+l|hXAQJHTwpYEot#ZT@;_pu05SM zA&DdSa7S9y4}K~CA$H>W&7&71`Z%k!!kMk6s~mhUfD=a&MtrY$S7nRv4R+L5t&q3g zOvLMC9`e$N?PjA7nOuIlia|OGGuQ0s#Vl<9cS>Kt-}tXOVa7*>j+e}I*p(>^BDFymj)UT9YiH3=;ei}?&kQ+N0iW!g-ZKDp&57q2|aLV8S@Zv(wH z`cZm(`XFOG_+Ek*j6$CPSOA;NFbUEuFM(%g$pky>TU;0%hhGJl$bpV}*Z>Kz#-rjm zI8KrEX|^}1*cUQcgy+_U9jWhTOr+vU#EsBLyEkhN-Q0c$m;}4y3%$eD?uMCy1!V?c z&g@O1lyAl!<)Gc%2YwfHspF|4hw_E({0FgZr zb~+uDMR)UKoVK5;_BP;E%7pb}cO6Qv(C=xW1yx|;yt6b6etotiRYE&hsrY=5#SmEr zHpI})NS@CIbaOlh&Bpqm#UswN51bojyiXjwuf3;1U?MISi z8jwK1g2 zqXBMnc^c$PK@nn%1}c;Qp1iAUL33SMk&aWl@S|;M$WE0|4omI8s1FlmqHW6 zR~8C5-ff4!+0a&aHvteD48)w6cCX`FPG7!H$&-h=&RSicyvTUnA&tdtdA1s7Rkk>; zsMjB87kv1_<4YhPq+xSWIMxL1_&TgdD40S_)&R$x{OkzS7}0a%BZs1$r#wS+5dQox za|p10B0qupv2=vJJ0Xe-gA2b35CS>iB@ z9irq@xy97eP3&wZ)XH-x=?c;X!c$hi(>wazfN}pKib1DP{Ry#DGZ1nBQ{BNC%S6dy z1~Z$*9@6xa7}RpOy)APSLOi}UpS`(f3Ls328a4xZf9BWF8t&L-^ck;rhO&_1|ALWw zJrqHtP>2P_+Ky{$+c5iTQLs0ffDWjig!tGZrWg*OeW-Wz6iOd9_&24?NqWjCS%Um( zFZ>f+$E&PL;J}?f9Q2gJZJOL)gY3KkR9W)ATqVJ|2hXv2`6XQlknZR4W^_0fyoeK- z!2Q}oh@e$~RV;R3RLIKDDjgd9{iJ*6awZP4bA?hBSQYuczXhj`H~11J?d%8w>{c}2 zV%ha0YK7Fa>eg2_pM$|cr!yHPrRcmD$`5R=4)8m7!`2ieFPSJ~A%ubI!X79v zPEd?e0zd*ghdEtc+M=IXjw}p1@|-pQh=xKh7$>rw z{;yVEmEb>_PMS4oFHtvVL^Md(LK}e)i*9n~Q>2hc_q-E8CF>LDL+4NtkPj5l=)JboFJ09pteDqq_sF5sEO7_RFmaKs%YfW61j%OD6`Hyo91Hi#*fg zPO<2wf5`xqbPZyxl^L*AtQ(+pUT{Ik(@SGf3Poax?(^>j6j~eD&Vb237&iW4yCf-x zb0Cj0iUACg?!ukSM_0iR%T?8R8@14nx?i4|XJ6H-27Mff8H!6q;biam4`P`NR&X{% z?|Uo^_hJC$#nIKG0G^fSJ)@g}+?(t8vk-&q8sDTs+=srYeB1OUYnP2@OWbM@>XM1p zkQ=SOhGIA;-+)jkZuUb3hBuTRs|IkrQPN$q`($nR#ecE@fO&{{eDdD&zb4F)r&k9e zsSRmlYJJ!ql&!!H2v@reaqTAGFyLWu-abx6>Cb*lh}WpyHK*Ste0AfvK*smIir?r+ zQi3|q%7i`j)fvLP)O615#g5TTq2)RTy8Z$fU#26@uSIEtzPF-q`$$>v#1f!8YSMT+ z)oa&{SY*V(o&Aw25eV)}b6rR`Y^JKp;)A)t4!b&P^3%{phXW6P20=&xl6X8Vg=iVY zf+qaMDRK}-dALRH>Y3G5Rgw%)O6#05kMNxvXIh55a0X^GsV=XM9tiOcMvLG3F!>qa z13V&$qWC8{%_k4(@KxO(jfT2>f<}Yyn17o)+t5YBACLwhKJkB<426CDWCGr@1GE04 z2KjOb@(=I=>1s=$7hhiG=v;FT!3UKw0@`H}-tkpSk! zU+ghsO|t}Dq%>Z?Fv0=q>Wg;`gTPLNxFR)pweQ}nO?0&e!yEKYla1D)_8^T;gv;1X z#TA(@%mZs_Sb%ZQ;p5_f+5(ghH3!o_OIn3Cs`~C&rlo00JPu#2kni+G0e>!1vk@1y zGfJ#@ni}J!+&OOqTY3zNUND9b`-pXw8Lnfn?Erdw>4Uz}uUkM~)Q4fl(`l-h#8Cf6 z&{F$JwP2j4&7!$LLt?Qt%jzkFzJL=sJ+$hP9Kx7rcZ!gCpAj6PWi%>I))Fi|GX50a z8TDQPMOc#o9FcH>H;t8Kv$IYPIVEK6tp9$Xo)!am`pg%`gIp+sb-~=l8B% zk1RRTJ)(x8Tws(Tbw{+r1|zz>Zi*I&nUTNC^H#gJO8_WcT{C;O#4RF(ggu>c&QWkAR^z2JQOgu93bM;|S0Z<&#>oXZ+7eC%dS6ZevNBz5{5D*Z5YVW zA;YzFd^(wFJ-z%QT5s=@8IIBw0X06#!O-KpL-T7e7D-MR&uE%YmYAMU=S>Wk$ner& zwT^t92f6{0AXKkg81{_6{~FkbG&`5Cpf(}(DU}$|cVeUk)7q`yaIucbrVTnO1Fo2H z24#`Ga5duASYQ5#+jV)kYBn4%qQ!}RF5H_6_<)I72^7TnZa}5T5Y`7AS7*yZxXRFM z)^T}g(T;?UVW#{=feNX6_1u6vf5^fqU}$W02DT)mnTZ6(qxRq@#MeluNMeG@&SST! zeCgpjyS$YU;)cx`pE!vj;F%TtPpx>$UT0%ePvP!#QJ2f;Z~Y_ghVM zI;y2e^EuKWD~S3X#`(w3V}vVCOkrc>aC`ID`Yn2wn9>V2o#A=fWyttFn2~b@1|Y|3+c%P(a#$hllVEd0+iRLMShi4du3R&SW=(>8eyU`+?zJGzXr^sg%+{ ztg=Fr7NCO+$XP)utr8ncM6PNOv9m3j5%ATXpKcTSu%9^nvrMN#U%f*znG+?wr$%SO zZFzq90OA*w-mgRdiiQaK-{AO`rhjx>WkCYc?uQ9xkr0>g!&2P2b|4o`>zKq^X0B}G z1zUjsLp?vRyQ;C#^wGqzKw62A_ES{g5~YyRIDpm2G!_>%;}GHkYJ?u!!vq^=N2Q2~ zx}E-TtqzLMmU7D1Vz+MHWkLet@MCN8{mRAfaP2()LF+I(D3RtkZWVw3$H|R)!m$(o zujSXdlSAo?s51v#!?|dMndqP$B4gleDxtVNY;`%^Yd`}CN`|@lp7Gptq>>hB2puBo z7VlHd!v`#ayrHeSQ)5zCvvSIduC4K6NEb+J1#y8P32z@abhrao03MsZ_|Xi~#McbgG5#^qt^1Z7DAze;c*^NoykGoPJAv8A zSvKs|pXMnz$w-#eMM>e}M-H>3G1(!;>7`(hlY#>^h#jMTiAu){T)U28OpKcBWJ^U$u;^Uy?UAKw#SHi^)Wk+(}1yCLrHmvHQ9sYC3z+j2NK6-4va*0$nL zJZ|L>jCM-0bhimh_G{vqK#H~MI!|pEmyS1oJ0Yb?R+Cu3&2UUNjrX+u>)B%XEU6Nb zgiMRtS$>o{?8uEwC*k++efDL`6be8~sG2z%n)`ck$C$4lnKVLuPPyin(@|8X+RS$2 zxG97#SrQ9U%X>H9!n=sddenPxMSfA2KWhk@amZMr(J%SE>AJ2w_fY=%L$Crbvm&LV z7&J{x(D=3t)2aDSMTdRFnYZ6{9wo*-AzmNFNiW(@cnM+kg*C2MABi)i5xmg#`FZho zI*vst;MzR1D7dlvz4?6D{DmJgm~%6=Ln!wUf0^j8hO6asT);+)0K*-d5-o=ValQs;`ozZEFIBg(bDHvn>>J7b$ltYI0e zfD1H09tk)QS2-c}7CFL+aTo`v#^Rm$KfNPETI?3^^4bXWWMI(5pGmW9mX%ybT40K@zgnDv4oNuHFQQt7Lb5^L#Sg zW!?C0!PdYWw)0>o^UIaGvSDqj6xG;3^LHq|JHE z_7Dp@@yv{GD%Uh<9@TwbcS9+9Aw0h9Gw@LtS~6qcuVlJr1w0HHw7^K_tX3#x80Y7@!>oC z9F`Ix4o%4oVi>@T+=3A($x^P4i&GP1OO4|S8CZ=abMAI<8P3vtnNIy{q;%b9(9Bn1 zHjeL}O3IcMEGhGSU;yI#)0^865?)=g#BCYDugW-_#Bt#K9AH1~ zJ&UwE`=W1z5!&6?1I^j;hxTwQv0~|@$D+gt%Lm`>#$BF7btz;EBoPQ4?n+ioi(WYc zb(uViK!9~O>euzQy=-iLGAXQ>VkcsMj#l=Ptb@83JnRH5x8Pk)7nYBA{A!?T8f;harMiL zfIRcjgsHS52;2{Soc|8>Q%8r(^60Z8dZ;N+-R@nYq>d$7uVN?qE^!`eB$)O1$a%~a zbv_KT7?w6zxK$iSBh|<hIXwQ-e=YUR2dy+p-tZd*k7_AGgX4w_fVZ8ImI9SJ zrr<+7;VdC8$~^IP!Y}LvdMs!wmiR z^pTgxs08vD_leAt1O3)0A~_4m{mEeNvVwq<*iolv?O?~pAj$L=GA3Og$!(`=3*7;g zGyN9qT!il)TM}rE$oIxU%2s3P@w3MPbe=4n0meN6(~VutaerHrHFiwDW4(rwpK_Wm z^7^$3{w0H}MYRPB7Uk^|zvNeHBD+u_;}HbfI6o3`MN+v^)os7(%z=n26U0{BELKE+ zRkAj1PXC;lbRKoEk5lz;lVx!>_Ifipea}rl6>>W^oNAZW7CHAxT(1;QchjlSA>G|l z`gA=09*uH7Js)6nPL-}ZknwKk*yCib*zwke72_jbxs`;k{cK7rd|yKKo9hcAJR7&0Y|<=h(l2Kp3eajL3~?^)#hoVH)J2)G#Ei z-P&iIEF(W!X6Si@rZY-7EI^;Hi=qqCzGiPeex7S}s#(u%9Y@^hc>^9iD%jlJC0|~Y zmd0%!^jw|4c>@!A)`vLyQ~4c|XMJFrxA=jjqoU8G^1{cVta0u>&>Mk@VJZF3y>o^{ z{9vu#E!y-WT_$Dm+PI4W?{fC5%@`#g%m^Dp1< zG{c5khh>*0i*p55*&PtjUa@Jgk=?N~Kiap9+E7!W@LYOloZju@cz6df{`3*^?A5`{ z+mgcUuiIUGv#`vGcafMX+h__0%_)y^bys~5Eg1|HOhcAR9@5Vgd~)3id%cjm9UnqW z!SxXC!GfMFJ3WA=H`5Oo+^)sT3pP;7NL1(X+S^X47gtB-dB{yJtkWQi*YCB7o|r(y z23pk1*sf&+_CofJN8W<$sj6Zv$>W8#TWwTE1XIHT)3&kt(ZYQOzm#|j~-8XR!qzgVgfI&KMuWTOb!1N$nU9t zDU$T3kl~7=@EdHvL&+?^2s$5A(iaXyys0M7e@PYXTF^T4#CmN!wx_>6^)!-97NzsC z%Vln+0gV>JFvtJme5dA0lA0_{*AAZ{ugPaWuAcba`C>9iffI^=jK%pKR*a$RTRjm@ zz0SJJjt%0-2vFN>o7dsFH(0cOn^GjH&PJmbea_!%q}P{2D>c4NbJ!>8C-DhfyKe@o zd=203J5v0GM+!VZFBihicvmAfZ8w_vF6)w`57Tnt>b(dhl=GDb1EC0_z|{S*usM6c(f6U=z%lNwXt#wxK3OB<}} zLjH>ct-f$^wpcEMhq5faXvp#gGB&$w$RqxCyYe(E<Ay{C;8n^YrTY0!a|*axgu(GC_?PCNy1284`T;--&Jn>V3S zKWYVw+?Plcb5kUs-KR*my>&HxAH*AAH z+fESX+#~NJGloEJ1jyVyYsi|Z)T`0~=FJg=P1(B;qYErx-y~lhYf3_XO;jJKIw9pA zyyk#4CqsYGCu8L1+~-s(<@r@|tV->Q)3+rifDu%#q}||{%_!Nh8$HmK1n3i?!EOMf ze`0dPw1j!KciiA-EZUEf!R4R~*^Rs~^n0z~e~h>J9+LG#7sFvbRWcE+Z=CPPKqg*IdiN3y1vtqi5k9&&z((@*JctGH`fz?s&N} zUOMfo?@9!=y)Rt92QU()FywbY@3ihz>rq#z*~m|b38rS>O+j_Ag?^uJKP`j2O3GX- zTU^(=!561&JBqg`$Qqzvcf^GB_GzBRZ~KP?5a8M zpZ}@NqN1Z3Wdg)qmX$dlqer#Z-b?o@$1+c>oy==5V+E{1_bbgEw#e^s-5hfk9(wiT zk=SOSYv4r->gxRG{B}6lsGYG6wA8$7`62gua)uRINIGzeHC)2ce%AqQN14s3>D{h! zHPRBcr2i)kL)#o1)p!=l(+_Eey9kTtTeY2_OJ^o#*j?Ocr@mto)a6(#Mmb>J+);*! zlQ4ij6=|lu-R+TOK}B>K`$P9Hr$8O9d>ejn&1fAzzIU`LaWvFi>+XgmnD!pCLthoK zuEO^$?Vkx@)7fFoCIkp?vzN;7cU1bYLJMC&)as(&Q^*9nzWUV#3qfGrPvUNAz)r!j zSA&M|gbS?4Q4=8O;OtZ|_`XTqv{pa{J~uB~>kTSnq(a&OOk`C3oxR}c@{vFITUZ3~ zrxsbmMIQ8CFmA&4!}2S6gH=6G>qPVY^#^gJS7Q%n_A0jMczJOP|ysQPPpozl}xx$m?L#) zKP`)Gy<+D>z$b7z(75F+>NV2F;r&B;PB~QWQpw2q1%h994Dcz$!n7XAsp7RLQEv9w zLn2E%uHy`5ns;aRG{=Ar;%gN+&xaQ^D8PdoTzatbZJpp&1F}P;WA(V{ImFj}n+j+f zgSgN1__d#7d8Ni6H!juet6CQ?BJ+#vQw13nB8r&DEp#@psSy{)Ii@7!pkNMh5A~kAPEe6CAiea3Taa>xk-ClRuschBFQX|w@6SS$ncpX zNhei~ukDA<9>UdQOq`MJup;=I*dFt>lIK9?hl zFo=p~RyKPC7v@b}^o`9fjQcq@zJa6bb5QR!3uk%Q4TD@(<`W9uT(|Ebt|FjQeIj>0>4`*Bnaj3yW{#Pswr^?jX>tFnvEX9 z{0PZ3T>Ftp@#a29TZHM8!|3n^atLmYHw^&< zxxl=C+i7O^%=>&$Z;mKLT|r#ZbBrQiWSLHv6UA?uT*IFWfo*L*z1I9sY+xGPj_C%> zww8*@O>$|Er)_EE{mrq}FE8x`%gdL7^^w9JMmx6f3B)fiXZ5ruDs<%;@K(y7^VcSZ zLt$KLBX=2&E~S`ATo|%DOX5C2nNMf)#+9Cah!B&5n2yZ3jna2kI!;2liz@D8u?^u} zqtOutJ*9$EUY`;G!$sE*BIeN@fQ`qyLmFXoXx!hZ7yZYnJR6@Pl&p+XZ^ zb4k@X&sk9k>jhlY=p8ZMSdPoe=i~SeM7Q2ak71x;4+}(Epva3`Bu7l|t2*cw77qfx zE!LQ@doDE(o1ms$%CETiU4+B8+;C--M5`V}xo)r8A z>v7t$SXTgeVqF*{HD}T5RM}NS%r}nKTF)Vq5n>0)8)Ff<+3OL(yvpArL<*H=d>V=? zHkEF%yyGtW!I^*32gJpASHGf-_lQh(nN|al>qkHtb=lQPtC~@P!$gR5%a>-Si%XWP zIpS9aCOe^Osg!yU=}e3!27|2k(||dki3~e}IvD4HiTK?<=D;mbPZ(+o3--{kSei`w zt~ShF7xi7f{--N7OAHCaFNl3oGXmA{wd2t9x?Tpx_RQ#-;@lTqs%IF81;$3yV;*H! z+h(w*doOQCs4ur9`!&= z9KW?E?|EpJT5f#pPwT`Gjwj+=c|S&{!t%`A;#NJ*oIxP*4?&tNiqpF6(?!)Zi+*vq zxc$+rUtLvuvurRS&K<{*(BbMNCx`Wsb^LvXS5ctu(l%Tq<^0MF3%Y3!Xuwj$Xh$T3 zb3!Hio7fHqlHST&AAB-nhGA?EtaVvEv6;;F@hS~QeE4(*B$A~v>5=_<%e)LNWr%{^UKk8?ZHRN=s*-0(c zuB?!xf-H%~n^Cr#SW`a$kar{|-T`kkSLJ(i8h`4hSdRLIo8xmml76S#Sy;hRG=5+G zer7@)m>O#scpXx3r3_hG;8B^*uS*8K@c*HaCD1hI0#kPz$q_hp24zj)Kv)tK@(Q<^ z$md)oMOOx44C0_oLZ$Ar!AP+IV+S3rq(EgyLcx(`j~Yij9QfG|Sf7>~1tm#nkF z@Rk9PT*>7PHqHO>TS##tRDG-X%T29lmCj~i1{W?HbL=G~5c@R;Sgbm2DAvk|SG=D8 zNtVoqh4Jp&)aH25VIm^PUqHt>p)`wGd{}0$Av%OG*MQr^-dK<$a?!)CITkHg;+!V$ z5V_%J!={rK2js#|IanFPc55Qc0&qtA=x>g_3-C_c<}<(_pA-#AMq)cnH_>(|;5*gb zv9wHnwmv`I`Jn?35*eoeT%(@(si34S@tl_5EQ-bshm*nn)`r;DLV!6#2|f9afDM+x znh(e|U!5I_Ayd?F=lOsDqGeRu9cy2=!sFk7TbCKE7KtpSDEz1VS#&_cr|of0e5<$9+YCoPl`0&2g&nzyuyF9$%gfSPH(e)!XKv(&Xz4*jA2`G`m zP1KCPPKQ&iFg%auZ|2zX!%r}7x&?dTtBWs6&??nNxDHKT2|AY+)mab~PJyEOmmRjw zJ_o{72D=oF8sdoU$GK7n-YKd-Q-vo{)2GMni^7UUl-p7U$TEm*uWzFKd5*-&QRPUv zNA-v<=|d}PX)#Udb-!1tS$KU;pcX$aGbKI#G9vDA?4Y_D9|@}-7Z?a01y@(v0JPC! z6-_OW{z@Rc`jJjf+2V7)>;Dvw&KsL`-Wm+^!7NZ%MV@jOd122{*B?+N_^bc?R>Aww8Yy zfZm`lR#{||R(`jy_bD3t=kL!ss%v~9%(&Lt*dJhtVzpWfo6>eaPyDwHfb+dZ2sZg$ znf`R^3CQ!QKE0%K*$)w>L5$HYJK*JZ^?zP`e248qzq;$^TlGwzyA@Vkh_J(R?RYI= zXNT)WG%Px}628_W=L?-Ne?DmZI51^*nBcSBLgTSBxcOPA@1i8>q3)}^e{5KNM7K5O z7NZ-?Xf%UsGH~p62)cGtUUDeknK?w7$}D)X4c<`iVA^9E9m&QEFk7uA!|OnAXCH|u zb>_apxt1XFUYox1pGS6E5U|hlR2DaLP`pas6rIj)8<6w*J_Sj;@ zBFEcg-S0&!CDrVUd0fT%cx7zW{UJX_sN%X!@69d0FMIJ!xsk1E zE?rW6Efn3YJ%ZGcuYjc&pyKPI;gyO0LvA$~%oSo^`gs?*dXo{3{` zG^f;x8%eQZw{K$I<$y0Vx{Giqdv}N$&7Kd!77NILd?N3@Y z=f=E{8u#AdeeP$B4+Nlad1FyG;iUq+H-&txWb|2Yb7UsnqV@zb;uiA4K4%vXgZr*0 zK*Y%7zz&kGu;#731_l)d&yFTUWi-hcY8%;7dW(Sc?9oBTwcHPAcAtjU@t;3G5#toP zRGw!+qO}gW)hQykW#egrX2xCPWtlc@+*DD$t*_oR;470~aSpvM51Gv+Q*KbSc~)*h z>RirjsvcqVY?%L#>)0u}dUG#LUjBB>b>YQV_+UJ5gKJ0stL_-X?;(*g&2je{1I%op z7XCGhiW2NN!}DD#@aGOV`@_CtG*@)>)8W~ACftll zr_&l4%H)z7jeEQw6Dt0obmBfYb9i+vv<-6i=4hM4oyZK!8{P<2ilez)(s%4dllNj< z66MQ^I!LFQ2e-p{C6s}z<~1u?>0O$`vcZ_yZb%3l7aF5?N_zgZNuBytzB03g z^f060T52KtVAeHcKw2I&dSc$XqW$Qib#llZaeTQbdBSXbGsL`>AQkHt$}Y5fqt?N| zL{h7CaVtz+jTk}CFd%|gZG2FEwNBA2coc+w4uQu6Qws_G&Ym*a(=290bd9qNq7^`5 z)4OBw^a9>ZSjsj0m(b2)_1eSvO`#o*us0r}qO7hE^$s0Zzp6(^eEF_#f>|MZs=1zm zI~y9-lLoS**wk*~68^QZMQ9>=?u(Wd2cAC#tAk}(2C(<88WcEj*s^8k`6>8-=T!Bp zZ)$=B)g4;7=6bqbAl^UXIJLZ(h;-xMkR?X<^+_rdKxsh9Q_-ruGn;-~lQK*X~F_5e_>z5M#ASd{q1UTU-`o z6Xa|ad4FV7`O{-H&{b_k1)h#OZdzj?cs+Xai*c)sa)c~jqJB&I_{nWe9Z#Yfboy@k zMeiPuGI6m}nTGi}2YcIalt2p*Yb{v5BV~2x^*Seh{MNj%@k`Ea3;NY55AD=G_bG zSfZEtCz4d1^#$1E!bobgeHIBKYwZJ_JfL}uR+YHS>T`%WWQe?Ov_iPw;J`1S>ItC< zScQ$dWeElBvx62NS0m*gaKbU`BzM7*c;_xmOb^Cpi=il)JsBT1JlQ5AJ(N)$!AKr4 zBoo$WuIMKp|CH%(@vGlTWvg?{97tTMu^Ni4f3bNxljH06auQjghpq()E!=byEB?+| z-2~VbEcQh@`iSi7yx){kg*DHmo@ZSQgp3jpxBk^-wxKhRGpQBkwCWXBi$oFk!CS*_ zc;C**Q{1qgvdP!=rX6c^P+#oSY~2(V2!I6rtas!a`j~%8DAs)L0i+8m+GleN+!#A%j*;N*O&w;GEWdJ-9X62K1lpNML3Yq z5kMY1*0-~_EaEgH4$)EGs{buU8DfhGL0~^_aE!_2 zNkTyZI}cA2rN3QrtYZlt=E^}v21XC9tkJ%`2%)tL}AC}4d9M6uaQIMqMJ+p z3D%73)Q;fBm^4ROD0JM=5%ANf>B_2)0tGMtrAOyC&_sAzJPH+i4{g}ZZt>FrMaS*e zk=fHp9dXQ~W39$hJY_0$MKHxK{u}#S%eiN3xf|Qp`eNqjj_6F<{Ry&X7}koN=Ej_c z-k!WyT~P1+UJ($PHAFCv!D@D-L~Q^b87A=#z_r!kzGg8W;&_qn+UxWOdVN6Kq>$V zjbI@U1P0`Gv;GD_MOBV}-1k&)a(nCub!X0ly|!49C<9L--(=sFR&nW%Hzyzx!d zf8ElQxmYuTVn9%OeI>`Q4}Hly3)J}Yzet}RAksVlc^WW)^f1wKNFJ~6VEp6EMYg~P z5pz{USvoQSV7#)tb;DEG;E(Z{CZL_{(>UeEw@C&lG9Bh0hQn_!U*D<#)2!;e#7Dl; zbSLAY8PQ>mrd}&H*f!mVEcasS0Wa3@GYK0W7|?Eiez2H~`|r)^idbrug}+lY{;R7k z_*R_Q_Up>?&eX4Ccp;mMu|?*XuFSCE>jGjAonT~@>$ihWiQpP}7<1bhSz zg4si6S4JRL#cDh}IZswJjoPRy4XH*if5c}Zy@vaQPjQF0ObRL1kV|(*=#A4yBaP}t zMeu4h!L7{VyYbDza?j?m1NIdrP4ZKB#~JgNn1Q8S?hbTY(Sm(lJ|>~~f{ zphC)xb^Q^0>Lli;ML9W@O0NC&Fef3S4xQA)cj}yFgAMV|v8=8)X-bU`C<`(gZjx%? zobl-yK=`wJ>Bmp!=|}4O_B$cvS%;_17A*Z)m->(?1G2v%pt?>SZSwxtgig)!H1aZa zDv_|vVX*zW-LWt*e96L|3#%}y;R|nbOI2vlIlJOvx-<3I=$WT_7fly?$=By~h^Orb z-?(o^HJ#Y?CG~9H8!}_T_N5fbH2ekmSyG{l4g3*F%=wR^zm)}W)Hj0HYZLv22xDJX zco#DF)fYGbdJbG=B>Z~rO4zNZOx5TQ2z!OcYV6Dt+-iv0Hqw3DZL9#n?Pt!7D& zqW#^2cM_VuTsB>!6T+ubp62{HifjuW0=XxH=dkRb z4OrAEvX*M3Nsww{VJCZ>UnW#CvF8x`@1%u~r_-g5&`7c8+)b|;;q~P}0iTcrwO{5A zjJ=U)Yy+*AL=PF`w!HOdZ3fd*#3M`_ih`l+E8s%U9m#Ecou0AAxt4_+DnnZQEN;zS zt^3gHE}(3Iv48b`4)kHjYtg*xc+U4A^74ir1P}y6y?J2k4Y+Wtjr3g&gx}Qysx_d6 zEN`Lc$10m0Xi))emapq1xNnLtMO+uY=jqB^RhPbj&!qTK66kJwV@>{spat@{v6ufJ zX};|jc2cKJE^E;>73Es!?xAzY)PH?!r3zEg{nKilzmGQ%5oG;;=C!TDga4~rJBS&s z1G@Gg+r*mfrTa<*hYObotG}OP-2)mY) zu@2wO1fXK=UQ+QgKNuJ;N`6b=d6LXHfVnCFf{=#}g#h3fe`3z_a?Fm{2tu8?I0>Bf z#Mmun%DfzjOi}ch{fGlp)WS%a8L%mmWfjnbLbA0?Rt_lV?37^hI+>hQpP36(G;$eQ zI$Ofb=pn4i(`R8D6SBqR@i9O>7I%glO%tPB-Zje~od~n&`=nOw@<+_DY0qGs5@aW1 z?-E36huRve)zVptjlCm`f*MLfDD4-;gLN1rE}SkADptz$N+F!7Y?+qeT61#Wc|ih8 zTF!7X+1rFff4>Uj*y>_~8TdUl3I)ExR|s>Biid1C4skd!mbr$EJ)l=F0931s;lH2( zJjOXtdA@G7wQYkC{SsfszENCBPfGv-CDn;n-(ED%C4`m{hNf(#v#0z|e))?TJy;@I zLGC*qPx8sWtpMa^2&OW4YY&_%E4=K2fC-6aKiPOOHd@etx4i0vIA>;BG{qzUHy#`z zDYp9C0WM(l7TnO;o?u3m76QrV5bb+*Zgvrjr!ty!W#y@CZ&i^l@|XF7b**s z@+$o*WB4YHuI<*LV4Az~-`x&|+W2jAb3_}O1q}nWV$$As^_3s9c~KtfR$?SDbT!H$ z7OHs|^lu<^S(l9fHQOnWmOV5+K0rCG5>Qs4MMJDd&YHykK+bMb*eeYl&g$&aOT{OpgX)-ep>4XLfw z6d)0yNZp9<)xzlCbk~Q7=BZJcP=5kk;cZj1lscu)n9^EHC z>{x~@$Hws_uDrUu4}&0~H+*^negYgYAfaczX!Nno|JeN9Kr7%6_Wc8(WUkM@KWxHj z4<76kA-qjKn5UMa)~e(pjY@eJBazdU^UhuR&Y8?K6x2G{%+ZQynl-|eZOylq7@yRi zAoFZ-d4XldTh=IomNxMc#9*RWrw+oq7b!&MSCq^@N^zf0DRm6#JsIqI*ObY91)sg5 zp%&fv4?SI~QVwWBo5yJwi*fIAQDZ)RI6cdC-GD@4UX`tGJNsFMHDPeb1B#F>%ryIe zTy^JyK14%y{>ZGN*gYh=sj{96AKac@5jjD4|8>q${CO-vM|$1Ub+(#nN|P{)c3V6! zvhcfCX=*OVr>-Z5st>(M%22b`cN^mrK{eKchNB!m4*Ml`5SD6E2E8d=NUG;5HU-14 zgkqJ|(B$RjNLxLZQ8*gkt-}3OZ9uoT){nl^p2ql!o-Rzex#zt1$;EUsLr8+KAWQH8%n{G22h<9+) zL^n`zvrHzWUxBY0CWn?Wx_5iGE<9i23Ctf42@Zwh$qO`FhF=PJE0%}j{XVTmesfvz z#1LoSNBdi1e&hHB+C*vWbe!M^5dI`H)d>x>zRZ>ed=hZCkf(Ce7S3_C@OP=HadaII zdy^rd{<|0^pHZ_8>M;c*6js-gw0)G^i03$wcY0`j2bhwFk&bwvLsjURxWMnbGPSzq z@+;UEbDO`8MXz2)ky~@Kh?D;XACcfUmXgcjvM9NDkm$9 zZ){X@NseZ8v8%V8o;fIwx2`3Ial$dz8Aq}efl^)I^{o+4O(P5MxR!_v!VZSNd4%I| zdy!A53F;E8UqQdt{2a}bk5vE&wyF>%MwPR0CnwJI-b^VB?K3GRA1Nlb3|&TUMT@BN zGEItV*$3ytqsyX_*$gd5VlHy8PLG}f49I<&`l)@#=2q4=@RzZar(FKkilPvz7m%pw(P@N6mA$j- zzX^#u1kD+{9elf4Lw=0x6&}UD)=3hG>)e(~%s>i{C34UMwj-}&Rktii+9(#f_?5p$ z6jAK2uWMeIFEua>2Y@L}PPW=EAm1wQRpj};xeg|FIUx&XXDzp=j1~H44(f@C;&y!1bwZU~E{Zh?I8stpLRO!7PtV7w51J$-(Twe)XzW zT2#RU7kvh^qPWx*rd#;g{#8uc(ci(3?;y@44j3fuPTBw(`)k7&=`g3jeS!C~@&qrh z7P|^KO|IKP*V|Y2n|v&o3(u*m`gepEW)hw0633G-Okz3C1g$HTAN(!2ok0Sy;3Vq^ zK6CixZVFyw$Q2s+0$wQ|nrd3oEiIKhctjJt2*RrA9b@CpAHIK~>gn~s2!deqCsQ?F ztuymP5ML$sLr5MauB0R4lefA$-C?wpz&aLqH!EOyOTi% z1+g_$=<9)NZUWY8udQ%lDK$uC~36WwhNqsjcXwy9a_ZH9Zx7dhDjnQ+bKRz z3Pu<#tojj_2X^+ENd}f;fXB%gGTu(%pxK`Ir6n{OG=~)CJCA;QQ{zAHg;v%W#1Q6^ z9g-ZMD@EzbNky2G5kU4}N@>K7RcI`5Fx3^mW!^X`_;HPG=5o;-Io)jlM9> zfWDv1y4BgIM%Z)>gEt|n+^+FGfS=;7wecwNzBTOz!(H5-2u+H@WdYGk`qA z<~*D^PhGKR+%C0f1dK^|f~274emyVYFH4L-W50hT!ll=78jurbH=8;;+&xdpkmHsR z58MbM{E{dBW@oob@p*Y+MNt$+>7Z5xx7{zevJ+7r#=&sDs-jR*X9>D$)+?M1W747A zLqd|1@iR?7BY+YxDedYSnWB@2f~Hx{jl{+u&nd8RIB*$nr8|J7&;fPqYz{5fmv-=> zn6{0ax7Oh+{=y7pugdv(0I#_xFz$sPEWK-#B(cV2Tbp}+9Tb?!OfUYM+4h$(FW-+( zH@vK`B07AgY*}?iN2m?SS};$CE5L_=2Oz+x&bXN0LUG~~*jMybw+yQCtck422P_Q) zn8A5sKIUNu460c?l?@W|b37b9wyQ63JxkED9Hj$lO(d}rzZFlTD8L)G%)PZLzox~V zrul#eL81UP5Mh?5Ju2pl$AN*?L7M+NS(EWoGDL&vy7Z&+PILB9jm%XZbU{Ai$}H>Y zS&myen0|!9Q_sFp`(wKwXX5>9L@C#*L3CdKk26%j-A|<^0j?RgPqTj zyiceeh4nz!`Hr%lI4FT@@4xtU=(^X|@>>mhAoe%VL85hYN|xGi%Zk(o_<_HX^?^Et z3K(!+GOyzHi=56gO4<=u$w`xptntDA&d_$mQ|{MR-d#Ljg&(QA_brHX_0b*tp8Nxa zZ9RC~L9H^hT08_-n(Nd;KwefH3$j5%0ULLYdgwmmnL- Date: Fri, 6 Nov 2026 07:17:50 -0800 Subject: [PATCH 2/8] Update WORDLIST --- inst/WORDLIST | 1 + 1 file changed, 1 insertion(+) diff --git a/inst/WORDLIST b/inst/WORDLIST index c725269d..1bc1dd90 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -64,6 +64,7 @@ pkgdown preformatted rceattle rda +rds rda's rdat rowid From b6ac0de8a81cb3f14b359660a6e001c65069fa5c Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 07:43:19 -0800 Subject: [PATCH 3/8] [feat]: observation-error intervals and log scale in `plot_index()` plot_index() adds `interval`, which draws intervals around the observed index from its observation error, and `log`, which plots the index on a log scale. Intervals assume lognormal errors, as in Stock Synthesis, FIMS, BAM, and Rceattle: uncertainty labeled "se", "sd", or "log_sd" is used as the log-scale standard error, and a CV is converted with sqrt(log(1 + CV^2)). Both default to the current behavior (no intervals, linear scale). --- R/plot_index.R | 80 +++++++++++++++++++++++++++++++- man/plot_index.Rd | 30 ++++++++++++ tests/testthat/test-plot_index.R | 59 +++++++++++++++++++++++ 3 files changed, 168 insertions(+), 1 deletion(-) diff --git a/R/plot_index.R b/R/plot_index.R index 03552cbc..693d16e0 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -9,13 +9,31 @@ #' @param unit_label String. Index of abundance/CPUE units #' #' Default: "" +#' @param interval Number. Coverage of intervals around the observed index +#' based on its observation error (e.g., 0.95), drawn as vertical lines. NULL +#' draws no intervals. +#' +#' Default: NULL +#' @param log Logical. TRUE/FALSE; plot the index on a log scale. Indices +#' with lognormal errors are often easier to judge on a log scale, where the +#' intervals are symmetric. +#' +#' Default: FALSE #' @param focus Character vector. Column name or names used to filter data. #' (i.e. select names of fleets to zoom into the plot) #' #' Default: NULL #' @returns A plot showing the expected and predicted index. #' -#' @details The input is from an assessment model output file +#' @details Intervals assume lognormal observation errors, as in Stock +#' Synthesis, FIMS, BAM, and Rceattle: \eqn{o \exp(\pm z \cdot \sigma)}, +#' where \eqn{o} is the observed index and \eqn{\sigma} is its log-scale +#' standard error. Uncertainty labeled "se", "sd", or "log_sd" is used as +#' \eqn{\sigma}, and a coefficient of variation ("cv") is converted to +#' \eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an +#' uncertainty have no interval. +#' +#' The input is from an assessment model output file #' translated to a standardized output (\link[stockplotr]{convert_output}). #' There are options to return a `ggplot2` object or export an .rda object #' containing associated caption and alternative text for the figure. @@ -31,6 +49,14 @@ #' interactive = FALSE #' ) #' +#' # 95% intervals around the observed index, on a log scale +#' plot_index( +#' dat = stockplotr::example_data, +#' interval = 0.95, +#' log = TRUE, +#' interactive = FALSE +#' ) +#' plot_index <- function( dat, unit_label = "", @@ -40,6 +66,8 @@ plot_index <- function( interactive = TRUE, module = NULL, focus = NULL, + interval = NULL, + log = FALSE, make_rda = FALSE, figures_dir = getwd(), ... @@ -73,6 +101,11 @@ plot_index <- function( module = module ) + # Observation-error intervals for the observed index + if (!is.null(interval)) { + prepared_data <- index_observed_interval(prepared_data, interval) + } + # Subset data if focus if (!is.null(focus)) { prepared_data <- prepared_data |> @@ -122,6 +155,31 @@ plot_index <- function( minor.ticks = TRUE ) ) + if (!is.null(interval)) { + obs <- prepared_data[ + grepl("index_observed", prepared_data$label) & !is.na(prepared_data$estimate_lower), + ] + # Draw the intervals beneath the observed and predicted index + intervals <- ggplot2::geom_linerange( + data = obs, + ggplot2::aes( + x = .data$year, + ymin = .data$estimate_lower, + ymax = .data$estimate_upper, + color = .data$group_var + ), + alpha = 0.6, + show.legend = FALSE + ) + plt$layers <- c(list(intervals), plt$layers) + } + if (log) { + # plot_obsvpred() expands the y-axis to include zero, which has no log + plt$layers <- Filter(function(l) !inherits(l$geom, "GeomBlank"), plt$layers) + plt <- suppressMessages( + plt + ggplot2::scale_y_log10(labels = scales::label_comma()) + ) + ggplot2::labs(y = "Estimated Index (log scale)") + } # Overwrite facets from base plot_timeseries bc scales need to be free if (!is.null(facet) && "fleet" %in% colnames(prepared_data)) { facet <- paste("~", paste(facet, collapse = " + ")) @@ -178,3 +236,23 @@ plot_index <- function( # Output final plot plt } + +# Lognormal intervals around the observed index from its log-scale standard +# error (uncertainty labeled "se", "sd", or "log_sd") or its coefficient of +# variation ("cv"); other rows get no interval +index_observed_interval <- function(dat, interval) { + if (!is.numeric(interval) || length(interval) != 1 || interval <= 0 || interval >= 1) { + cli::cli_abort("{.arg interval} must be a number between 0 and 1 (e.g., 0.95).") + } + z <- stats::qnorm(1 - (1 - interval) / 2) + type <- tolower(as.character(dat$uncertainty_label)) + sigma <- dplyr::case_when( + type %in% c("se", "sd", "log_sd", "logsd", "stddev") ~ dat$uncertainty, + type == "cv" ~ sqrt(log(1 + dat$uncertainty^2)), + TRUE ~ NA_real_ + ) + observed <- grepl("index_observed", dat$label) & !is.na(sigma) + dat$estimate_lower <- ifelse(observed, dat$estimate * exp(-z * sigma), NA_real_) + dat$estimate_upper <- ifelse(observed, dat$estimate * exp(z * sigma), NA_real_) + dat +} diff --git a/man/plot_index.Rd b/man/plot_index.Rd index 19572601..bf069772 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -12,6 +12,8 @@ plot_index( interactive = TRUE, module = NULL, focus = NULL, + interval = NULL, + log = FALSE, make_rda = FALSE, figures_dir = getwd(), ... @@ -58,6 +60,18 @@ module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}} Default: NULL} +\item{interval}{Number. Coverage of intervals around the observed index +based on its observation error (e.g., 0.95), drawn as vertical lines. NULL +draws no intervals. + +Default: NULL} + +\item{log}{Logical. TRUE/FALSE; plot the index on a log scale. Indices +with lognormal errors are often easier to judge on a log scale, where the +intervals are symmetric. + +Default: FALSE} + \item{make_rda}{Logical. TRUE/FALSE; indicate whether to save the object and make an automated caption and alternative text in the form of an \code{rda} object. If TRUE, the rda will be exported to the folder indicated in the argument "figures_dir". @@ -79,6 +93,14 @@ A plot showing the expected and predicted index. Plot Index of Abundance } \details{ +Intervals assume lognormal observation errors, as in Stock +Synthesis, FIMS, BAM, and Rceattle: \eqn{o \exp(\pm z \cdot \sigma)}, +where \eqn{o} is the observed index and \eqn{\sigma} is its log-scale +standard error. Uncertainty labeled "se", "sd", or "log_sd" is used as +\eqn{\sigma}, and a coefficient of variation ("cv") is converted to +\eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an +uncertainty have no interval. + The input is from an assessment model output file translated to a standardized output (\link[stockplotr]{convert_output}). There are options to return a \code{ggplot2} object or export an .rda object @@ -91,6 +113,14 @@ plot_index( interactive = FALSE ) +# 95\% intervals around the observed index, on a log scale +plot_index( + dat = stockplotr::example_data, + interval = 0.95, + log = TRUE, + interactive = FALSE +) + } \seealso{ \code{\link[=convert_output]{convert_output()}}, \code{\link[=filter_data]{filter_data()}}, \code{\link[=process_data]{process_data()}}, \code{\link[=plot_obsvpred]{plot_obsvpred()}}, \code{\link[=export_kqs]{export_kqs()}}, \code{\link[=insert_kqs]{insert_kqs()}}, \code{\link[=create_rda]{create_rda()}} diff --git a/tests/testthat/test-plot_index.R b/tests/testthat/test-plot_index.R index 2a5239a3..eefee85a 100644 --- a/tests/testthat/test-plot_index.R +++ b/tests/testthat/test-plot_index.R @@ -30,3 +30,62 @@ test_that("plot_index generates plots without errors", { # # expect that both figures dir and the index_figure.rda file exist # check_and_remove_tmp_files("index_figure.rda") # }) + +test_that("plot_index draws lognormal intervals around the observed index", { + p <- plot_index(stockplotr::example_data, interval = 0.95, interactive = FALSE) + # The intervals are the first layer, beneath the points and lines + expect_s3_class(p$layers[[1]]$geom, "GeomLinerange") + ints <- p$layers[[1]]$data + obs <- stockplotr::example_data[ + stockplotr::example_data$label %in% "index_observed", + ] + obs <- obs[order(obs$fleet, obs$year), ] + ints <- ints[order(ints$fleet, ints$year), ] + expect_equal(nrow(ints), nrow(obs)) + expect_equal(ints$estimate_lower, obs$estimate * exp(-stats::qnorm(0.975) * obs$uncertainty)) + expect_equal(ints$estimate_upper, obs$estimate * exp(stats::qnorm(0.975) * obs$uncertainty)) + expect_no_error(ggplot2::ggplot_build(p)) + + # No intervals by default + p <- plot_index(stockplotr::example_data, interactive = FALSE) + expect_false(any(sapply(p$layers, function(l) inherits(l$geom, "GeomLinerange")))) + + expect_error(plot_index(stockplotr::example_data, interval = 95, interactive = FALSE), "between 0 and 1") +}) + +test_that("plot_index intervals use log-scale SDs and CVs", { + dat <- data.frame( + label = "index_observed", year = 1:3, fleet = "a", + estimate = c(10, 20, 30), uncertainty = c(0.2, 0.2, NA), + uncertainty_label = c("cv", "log_sd", "cv") + ) + out <- index_observed_interval(dat, interval = 0.9) + z <- stats::qnorm(0.95) + expect_equal(out$estimate_lower[1], 10 * exp(-z * sqrt(log(1 + 0.2^2)))) + expect_equal(out$estimate_upper[2], 20 * exp(z * 0.2)) + # No uncertainty, no interval + expect_true(is.na(out$estimate_lower[3])) + + # Rceattle (log_sd) and BAM (cv) output + for (dat in list(stockplotr::pollock, stockplotr::vermilion_snapper)) { + p <- suppressMessages(plot_index(dat, interval = 0.95, interactive = FALSE)) + expect_gt(nrow(p$layers[[1]]$data), 0) + expect_no_error(ggplot2::ggplot_build(p)) + } +}) + +test_that("plot_index can use a log scale", { + p <- plot_index(stockplotr::example_data, interval = 0.95, log = TRUE, interactive = FALSE) + expect_equal(p$scales$get_scales("y")$trans$name, "log-10") + expect_equal(p$labels$y, "Estimated Index (log scale)") + expect_no_error(ggplot2::ggplot_build(p)) +}) + +test_that("plot_index intervals work for FIMS output", { + est <- readRDS(test_path("fixtures", "fims_output", "fims_estimates.rds")) + fims <- suppressMessages(convert_output(est, model = "fims", start_year = 2001)) + p <- suppressMessages(plot_index(fims, interval = 0.95, interactive = FALSE)) + ints <- p$layers[[1]]$data + rows <- est[est$label == "index_expected" & !is.na(est$observed), ] + expect_equal(ints$estimate_lower, rows$observed * exp(-stats::qnorm(0.975) * exp(rows$log_sd))) +}) From 4ed72c9815d1dde2ce3ce640dd635bb4cb4c48a8 Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 6 Nov 2026 11:28:14 -0800 Subject: [PATCH 4/8] Update plot_index.R Make default interval 0.95 per @iantaylor-noaa --- R/plot_index.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/R/plot_index.R b/R/plot_index.R index 693d16e0..f1466a72 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -10,8 +10,8 @@ #' #' Default: "" #' @param interval Number. Coverage of intervals around the observed index -#' based on its observation error (e.g., 0.95), drawn as vertical lines. NULL -#' draws no intervals. +#' based on its observation error (e.g., 0.95), drawn as vertical lines. +#' Defaults to 0.95; NULL draws no intervals. #' #' Default: NULL #' @param log Logical. TRUE/FALSE; plot the index on a log scale. Indices @@ -66,7 +66,7 @@ plot_index <- function( interactive = TRUE, module = NULL, focus = NULL, - interval = NULL, + interval = 0.95, log = FALSE, make_rda = FALSE, figures_dir = getwd(), @@ -211,7 +211,7 @@ plot_index <- function( cpue.min, cpue.max, cpue.units, - dir=figures_dir + dir = figures_dir ) # Add key quantities to captions/alt text @@ -221,7 +221,7 @@ plot_index <- function( cpue.min, cpue.max, cpue.units, - dir=figures_dir + dir = figures_dir ) create_rda( From 619247bd35f685cd6a2ee47d8b9f03c82eab1500 Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 6 Nov 2026 11:47:07 -0800 Subject: [PATCH 5/8] Update plot_index.Rd redoc --- man/plot_index.Rd | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/man/plot_index.Rd b/man/plot_index.Rd index bf069772..713f4dd5 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -12,7 +12,7 @@ plot_index( interactive = TRUE, module = NULL, focus = NULL, - interval = NULL, + interval = 0.95, log = FALSE, make_rda = FALSE, figures_dir = getwd(), @@ -61,8 +61,8 @@ module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}} Default: NULL} \item{interval}{Number. Coverage of intervals around the observed index -based on its observation error (e.g., 0.95), drawn as vertical lines. NULL -draws no intervals. +based on its observation error (e.g., 0.95), drawn as vertical lines. +Defaults to 0.95; NULL draws no intervals. Default: NULL} From e0af1fa058dbd1acb3abc651350515c596f5b82e Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 17:04:47 -0800 Subject: [PATCH 6/8] Bring same files over to update per commit # 467 --- R/conout_fims.R | 38 +++++--- inst/resources/fims_var_names.csv | 114 +++++++++++----------- tests/testthat/test-convert_output_fims.R | 33 ++++++- 3 files changed, 115 insertions(+), 70 deletions(-) diff --git a/R/conout_fims.R b/R/conout_fims.R index 33d8c3bb..490000b1 100644 --- a/R/conout_fims.R +++ b/R/conout_fims.R @@ -17,14 +17,19 @@ #' is one row in FIMS, with the observed value in `observed` and the expected #' value in `estimated`. These are split into two rows, labeled #' "*_observed" and "*_predicted" as for other models (e.g., "index_observed" -#' and "index_predicted"). Observed rows hold the observation error -#' (`exp(log_sd)`, a log-scale standard error for lognormal data) in +#' and "index_predicted"). Observed rows hold the observation error in #' `uncertainty` and the negative log-likelihood of the observation in -#' `likelihood`. Observations coded as missing (-999) are not observations -#' and are not included. Compositions are kept as counts, as reported by FIMS. +#' `likelihood`. FIMS reports the observation error as `log_sd`, the log of +#' its standard deviation. For lognormal data this standard deviation +#' (\eqn{\sigma}) is on the log scale, so it is converted to the coefficient +#' of variation of the observation, \eqn{\sqrt{\exp(\sigma^2) - 1}}, and +#' labeled "cv"; for other data `exp(log_sd)` is labeled "se". Observations +#' coded as missing (-999) are not observations and are not included. +#' Compositions are kept as counts, as reported by FIMS. #' #' All other estimates keep one row each, with the value in `estimate` and -#' its standard error in `uncertainty`. Rows in years after the last year of +#' its standard error in `uncertainty`. FIMS codes estimates that could not +#' be calculated as -999; these are NA. Rows in years after the last year of #' data (e.g., spawning biomass at the start of the year after the model #' ends) have `era` "fore"; other rows with a year have `era` "time". #' @@ -60,14 +65,14 @@ conout_fims <- function(file, out_new, start_year = NULL) { year_i > last_year ~ "fore", TRUE ~ "time" ) - uncertainty <- suppressWarnings(as.numeric(col("uncertainty"))) + uncertainty <- dplyr::na_if(suppressWarnings(as.numeric(col("uncertainty"))), -999) base <- data.frame( module_name = as.character(col("module_name")), label = as.character(dat$label), era = era, year = year, initial = suppressWarnings(as.numeric(col("input"))), - estimate = suppressWarnings(as.numeric(dat$estimated)), + estimate = dplyr::na_if(suppressWarnings(as.numeric(dat$estimated)), -999), uncertainty = uncertainty, uncertainty_label = ifelse(is.na(uncertainty), NA_character_, "se"), likelihood = NA_real_, @@ -91,14 +96,21 @@ conout_fims <- function(file, out_new, start_year = NULL) { is_data <- !is.na(stream) & has_data %in% TRUE obs_rows <- base[is_data & !is.na(observed), , drop = FALSE] if (nrow(obs_rows) > 0) { - log_sd <- suppressWarnings(as.numeric(col("log_sd")))[is_data & !is.na(observed)] - obs_rows$label <- paste0(stream[is_data & !is.na(observed)], "_observed") - obs_rows$estimate <- observed[is_data & !is.na(observed)] - obs_rows$uncertainty <- exp(log_sd) - obs_rows$uncertainty_label <- ifelse(is.na(log_sd), NA_character_, "se") + is_obs <- is_data & !is.na(observed) + sigma <- exp(suppressWarnings(as.numeric(col("log_sd")))[is_obs]) + # A lognormal sigma is on the log scale; report the CV of the observation + lognormal <- grepl("^(d?lnorm|log_?normal)$", col("distribution")[is_obs], ignore.case = TRUE) + obs_rows$label <- paste0(stream[is_obs], "_observed") + obs_rows$estimate <- observed[is_obs] + obs_rows$uncertainty <- ifelse(lognormal, sqrt(exp(sigma^2) - 1), sigma) + obs_rows$uncertainty_label <- dplyr::case_when( + is.na(sigma) ~ NA_character_, + lognormal ~ "cv", + TRUE ~ "se" + ) obs_rows$initial <- NA_real_ # Stock Synthesis reports the negative log-likelihood of each observation - obs_rows$likelihood <- -suppressWarnings(as.numeric(col("likelihood")))[is_data & !is.na(observed)] + obs_rows$likelihood <- -suppressWarnings(as.numeric(col("likelihood")))[is_obs] } pred_rows <- base[is_data, , drop = FALSE] pred_rows$label <- paste0(stream[is_data], "_predicted") diff --git a/inst/resources/fims_var_names.csv b/inst/resources/fims_var_names.csv index cc197f51..7b3c2cd1 100644 --- a/inst/resources/fims_var_names.csv +++ b/inst/resources/fims_var_names.csv @@ -1,55 +1,59 @@ -module_name,label,alt_label -Fleet,age_to_length_conversion, -Fleet,agecomp_expected,age_composition_expected -Fleet,agecomp_proportion,age_composition_proportion -Fleet,index_expected, -Fleet,index_naa,index -Fleet,index_nal,index -Fleet,index_numbers, -Fleet,index_waa,index_weight -Fleet,index_weight, -Fleet,landings_expected, -Fleet,landings_naa,landings -Fleet,landings_nal,landings -Fleet,landings_numbers,landings_numbers -Fleet,landings_waa,landings_weight -Fleet,landings_weight,landings_weight -Fleet,lengthcomp_expected,length_composition_expected -Fleet,lengthcomp_proportion,length_composition_proportion -Fleet,log_Fmort,log_fishing_mortality -Fleet,log_q,log_catchability -Population,Biomass,biomass -Population,NAA,abundance -Population,Recruitment,recruitment -Population,SSB,spawning_biomass -Population,log_M,log_natural_mortality -Population,log_init_naa,log_abundance_init -maturity,inflection_point, -maturity,slope, -recruitment,log_devs,log_recruitment_deviations -recruitment,log_rzero,log_recruitment_unfished -recruitment,logit_steep, -selectivity,inflection_point, -selectivity,slope, -,AgeCompositionExpected,age_composition_expected -,AgeCompositionProportion,age_composition_proportion -,FMort,fishing_mortality -,IndexExpected,index_expected -,IndexNumberAtAge,index_numbers -,IndexNumberAtLength,index_numbers -,LandingsExpected,landings_expected -,LandingsNumberAtAge,landings_numbers -,LandingsNumberAtLength,landings_numbers -,LengthCompositionExpected,length_composition_expected -,LengthCompositionProportion,length_composition_proportion -,LogRecDev,log_recruitment_deviations -,Q,catchability -Population,expected_recruitment,recruitment -Population,numbers_at_age,abundance -Population,unfished_spawning_biomass,spawning_biomass_unfished -Population,unfished_biomass,biomass_unfished -Population,mortality_M,natural_mortality -Recruitment,log_devs,log_recruitment_deviations -Recruitment,log_rzero,log_recruitment_unfished -Recruitment,log_sd,log_recruitment_sd -Recruitment,logit_steep,logit_steepness +module_name,label,alt_label,unit +Fleet,age_to_length_conversion,, +Fleet,agecomp_expected,age_composition_expected, +Fleet,agecomp_proportion,age_composition_proportion, +Fleet,index_expected,, +Fleet,index_naa,index,numbers +Fleet,index_nal,index,numbers +Fleet,index_numbers,,numbers +Fleet,index_waa,index_weight,mt +Fleet,index_weight,,mt +Fleet,landings_expected,, +Fleet,landings_naa,landings,numbers +Fleet,landings_nal,landings,numbers +Fleet,landings_numbers,landings_numbers,numbers +Fleet,landings_waa,landings_weight,mt +Fleet,landings_weight,landings_weight,mt +Fleet,lengthcomp_expected,length_composition_expected, +Fleet,lengthcomp_proportion,length_composition_proportion, +Fleet,log_Fmort,log_fishing_mortality, +Fleet,log_q,log_catchability, +Population,Biomass,biomass,mt +Population,NAA,abundance,numbers +Population,Recruitment,recruitment,numbers +Population,SSB,spawning_biomass,mt +Population,log_M,log_natural_mortality, +Population,log_init_naa,log_abundance_init, +maturity,inflection_point,, +maturity,slope,, +recruitment,log_devs,log_recruitment_deviations, +recruitment,log_rzero,log_recruitment_unfished, +recruitment,logit_steep,, +selectivity,inflection_point,, +selectivity,slope,, +,AgeCompositionExpected,age_composition_expected, +,AgeCompositionProportion,age_composition_proportion, +,FMort,fishing_mortality, +,IndexExpected,index_expected, +,IndexNumberAtAge,index_numbers, +,IndexNumberAtLength,index_numbers, +,LandingsExpected,landings_expected, +,LandingsNumberAtAge,landings_numbers, +,LandingsNumberAtLength,landings_numbers, +,LengthCompositionExpected,length_composition_expected, +,LengthCompositionProportion,length_composition_proportion, +,LogRecDev,log_recruitment_deviations, +,Q,catchability, +Population,expected_recruitment,recruitment,numbers +Population,numbers_at_age,abundance,numbers +Population,unfished_spawning_biomass,spawning_biomass_unfished,mt +Population,unfished_biomass,biomass_unfished,mt +Population,mortality_M,natural_mortality, +Recruitment,log_devs,log_recruitment_deviations, +Recruitment,log_rzero,log_recruitment_unfished, +Recruitment,log_sd,log_recruitment_sd, +Recruitment,logit_steep,logit_steepness, +Population,biomass,,mt +Population,spawning_biomass,,mt +Fleet,landings_observed,,mt +Fleet,landings_predicted,,mt diff --git a/tests/testthat/test-convert_output_fims.R b/tests/testthat/test-convert_output_fims.R index a165b3a4..1b25ecc7 100644 --- a/tests/testthat/test-convert_output_fims.R +++ b/tests/testthat/test-convert_output_fims.R @@ -26,8 +26,9 @@ test_that("FIMS data points are split into observed and predicted rows", { expect_equal(obs$estimate, rows$observed) expect_equal(pred$estimate, rows$estimated) # Observation error and negative log-likelihood of each observation - expect_equal(obs$uncertainty, exp(rows$log_sd)) - expect_equal(unique(obs$uncertainty_label), "se") + # Lognormal observation error is reported as a CV + expect_equal(obs$uncertainty, sqrt(exp(exp(rows$log_sd)^2) - 1)) + expect_equal(unique(obs$uncertainty_label), "cv") expect_equal(obs$likelihood, -rows$likelihood) expect_equal(unique(obs$fleet), "survey1") expect_equal(obs$year, rows$year_i + 2000) @@ -55,6 +56,34 @@ test_that("FIMS observations coded as missing are not observations", { expect_equal(sum(out$label == "index_predicted"), length(idx)) }) +test_that("FIMS lognormal intervals are on the original scale", { + est <- fims_estimates() + est <- est[est$label == "landings_expected" | est$label == "catch_expected", ] + row <- which(!is.na(est$observed))[1] + est$observed[row] <- 100 + est$log_sd[row] <- log(0.2) + out <- convert_fims(est) + obs <- out[out$label == "landings_observed", ][1, ] + expect_equal(obs$uncertainty_label, "cv") + # The CV gives back the lognormal sigma, and so the lognormal interval + sigma <- sqrt(log(1 + obs$uncertainty^2)) + expect_equal(sigma, 0.2) + expect_equal(100 * exp(c(-1.96, 1.96) * sigma), c(67.57, 147.99), tolerance = 1e-4) +}) + +test_that("FIMS estimates that could not be calculated are NA", { + est <- fims_estimates() + ssb <- which(est$label == "spawning_biomass") + est$estimated[ssb[1]] <- -999 + est$uncertainty[ssb[1]] <- -999 + out <- convert_fims(est) + out_ssb <- out[out$label == "spawning_biomass", ] + # The row is kept, with a missing estimate + expect_equal(nrow(out_ssb), length(ssb)) + expect_true(is.na(out_ssb$estimate[1]) && is.na(out_ssb$uncertainty[1])) + expect_false(any(out$estimate %in% -999)) +}) + test_that("FIMS estimates are relabeled and get years and eras", { est <- fims_estimates() out <- convert_fims(est, start_year = 2001) From c8109c0fd2a1d73a63e374d6c815f4770f47ac81 Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 17:37:42 -0800 Subject: [PATCH 7/8] Forgot to change the test --- R/plot_index.R | 4 ++-- man/plot_index.Rd | 4 ++-- tests/testthat/test-plot_index.R | 6 ++++-- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/R/plot_index.R b/R/plot_index.R index b2944874..7040746c 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -11,9 +11,9 @@ #' Default: "" #' @param interval Number. Coverage of intervals around the observed index #' based on its observation error (e.g., 0.95), drawn as vertical lines. -#' Defaults to 0.95; NULL draws no intervals. +#' NULL draws no intervals. #' -#' Default: NULL +#' Default: 0.95 #' @param log Logical. TRUE/FALSE; plot the index on a log scale. Indices #' with lognormal errors are often easier to judge on a log scale, where the #' intervals are symmetric. diff --git a/man/plot_index.Rd b/man/plot_index.Rd index 4ad56611..c1902a61 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -62,9 +62,9 @@ Default: NULL} \item{interval}{Number. Coverage of intervals around the observed index based on its observation error (e.g., 0.95), drawn as vertical lines. -Defaults to 0.95; NULL draws no intervals. +NULL draws no intervals. -Default: NULL} +Default: 0.95} \item{log}{Logical. TRUE/FALSE; plot the index on a log scale. Indices with lognormal errors are often easier to judge on a log scale, where the diff --git a/tests/testthat/test-plot_index.R b/tests/testthat/test-plot_index.R index eefee85a..4d7e2299 100644 --- a/tests/testthat/test-plot_index.R +++ b/tests/testthat/test-plot_index.R @@ -46,8 +46,10 @@ test_that("plot_index draws lognormal intervals around the observed index", { expect_equal(ints$estimate_upper, obs$estimate * exp(stats::qnorm(0.975) * obs$uncertainty)) expect_no_error(ggplot2::ggplot_build(p)) - # No intervals by default - p <- plot_index(stockplotr::example_data, interactive = FALSE) + # 95% intervals by default, and none with interval = NULL + p_default <- plot_index(stockplotr::example_data, interactive = FALSE) + expect_equal(p_default$layers[[1]]$data, p$layers[[1]]$data) + p <- plot_index(stockplotr::example_data, interval = NULL, interactive = FALSE) expect_false(any(sapply(p$layers, function(l) inherits(l$geom, "GeomLinerange")))) expect_error(plot_index(stockplotr::example_data, interval = 95, interactive = FALSE), "between 0 and 1") From ac9b5139d0817d28dedda64aa35dc6edd718bfd6 Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 13 Nov 2026 10:41:25 -0800 Subject: [PATCH 8/8] Changed default from "" to NULL in documentation Fixed and better documented 3 logic pathways --- R/plot_index.R | 49 ++++++++++++++++++++++++----------------------- man/plot_index.Rd | 24 +++++++---------------- 2 files changed, 32 insertions(+), 41 deletions(-) diff --git a/R/plot_index.R b/R/plot_index.R index a7746c51..3111531a 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -1,30 +1,19 @@ #' Plot Index of Abundance #' #' @inheritParams plot_spawning_biomass -#' @param facet Character vector. Column name or names used for faceting. "fleet" is always added on to any faceting selections -#' -#' Default: NULL -#' -#' Options: May include, but is not limited to: "fleet", "year", "area", NULL -#' @param unit_label String. Index of abundance/CPUE units -#' -#' Default: "" +#' @param facet Character vector. Column name or names used for faceting. "fleet" +#' is always added on to any faceting selections. Default: NULL, but may include, +#' but is not limited to: "fleet", "year", "area", NULL +#' @param unit_label String. Index of abundance/CPUE units. Default: NULL #' @param interval Number. Coverage of intervals around the observed index #' based on its observation error (e.g., 0.95), drawn as vertical lines. -#' NULL draws no intervals. -#' -#' Default: 0.95 +#' NULL draws no intervals, default = 0.95 #' @param log Logical. TRUE/FALSE; plot the index on a log scale. Indices #' with lognormal errors are often easier to judge on a log scale, where the -#' intervals are symmetric. -#' -#' Default: FALSE +#' intervals are symmetric. Default = FALSE #' @param focus Character vector. Column name or names used to filter data. -#' (i.e. select names of fleets to zoom into the plot) -#' -#' Default: NULL +#' (i.e. select names of fleets to zoom into the plot). Default = NULL #' @returns A plot showing the expected and predicted index. -#' #' @details Intervals assume lognormal observation errors, as in Stock #' Synthesis, FIMS, BAM, and Rceattle: \eqn{o \exp(\pm z \cdot \sigma)}, #' where \eqn{o} is the observed index and \eqn{\sigma} is its log-scale @@ -93,11 +82,23 @@ plot_index <- function( interactive = interactive, module = module ) - if (!is.null(unit_label) && !is.na(unique(prepared_data$unit))) unit_label <- unique(prepared_data$unit) - if (is.na(unit_label) || is.null(unit_label)) { - cli::cli_alert_warning("No unit label found in the data. Please provide a unit label using the 'unit_label' argument.") - unit_label <- "" + + # There's 3 cases handled here + # plot_index(dat) (unit_label = NULL): Uses prepared_data$unit if present, otherwise warns / falls back to "" + # plot_index(dat, unit_label = "fish/hr"): Skips the is.null() block and keeps "fish/hr" + # plot_index(dat, unit_label = ""): Skips the is.null() block and intentionally leaves the y-axis without units ("Estimated Index") + if (is.null(unit_label)) { + data_units <- unique(stats::na.omit(prepared_data$unit)) + if (length(data_units) == 1 && nzchar(data_units)) { + unit_label <- data_units + } else { + cli::cli_alert_warning( + "No unit label found in the data. Please provide a unit label using the 'unit_label' argument." + ) + unit_label <- "" + } } + # Set cpue unit label for plot u_units <- ifelse( is.null(unit_label) | unit_label == "", @@ -237,11 +238,11 @@ plot_index <- function( dir = figures_dir ) } - + if (exists("selected_module")) { rm(selected_module, envir = .GlobalEnv) } - + plt } diff --git a/man/plot_index.Rd b/man/plot_index.Rd index f9e7786a..27152f2d 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -26,9 +26,7 @@ returned from \link[stockplotr]{convert_output}. If inputting a list of tibbles, the first tibble's reference point defined in \code{ref_line} is used to plot a reference line or calculate relative spawning biomass.} -\item{unit_label}{String. Index of abundance/CPUE units - -Default: NULL} +\item{unit_label}{String. Index of abundance/CPUE units. Default: NULL} \item{group}{String. Single column that groups the data. @@ -37,11 +35,9 @@ Set group = "none" to summarize data over all indexing values. Default: NULL Options: Including, but not limited to: "year", "area", "fleet", "sex", "none", NULL} -\item{facet}{Character vector. Column name or names used for faceting. "fleet" is always added on to any faceting selections - -Default: NULL - -Options: May include, but is not limited to: "fleet", "year", "area", NULL} +\item{facet}{Character vector. Column name or names used for faceting. "fleet" +is always added on to any faceting selections. Default: NULL, but may include, +but is not limited to: "fleet", "year", "area", NULL} \item{interactive}{Logical. TRUE/FALSE; indicate whether the environment is interactive. @@ -56,21 +52,15 @@ If the interactive and >1 module_name is found, user will select the module_name in the console. @seealso \code{\link[=filter_data]{filter_data()}}} \item{focus}{Character vector. Column name or names used to filter data. -(i.e. select names of fleets to zoom into the plot) - -Default: NULL} +(i.e. select names of fleets to zoom into the plot). Default = NULL} \item{interval}{Number. Coverage of intervals around the observed index based on its observation error (e.g., 0.95), drawn as vertical lines. -NULL draws no intervals. - -Default: 0.95} +NULL draws no intervals, default = 0.95} \item{log}{Logical. TRUE/FALSE; plot the index on a log scale. Indices with lognormal errors are often easier to judge on a log scale, where the -intervals are symmetric. - -Default: FALSE} +intervals are symmetric. Default = FALSE} \item{make_rda}{Logical. TRUE/FALSE; indicate whether to save the object and make an automated caption and alternative text in the form of an \code{rda} object. If TRUE,