From eaefc95849d97d4329e9774ec8fb0b5c09006c31 Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 07:12:49 -0800 Subject: [PATCH 01/10] [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 02/10] 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 03/10] [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 2f53ca08d97ebcf4640d0355a8c0b5659d281e64 Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 08:00:33 -0800 Subject: [PATCH 04/10] [feat]: plot all indices relative to their means with `plot_index(relative = TRUE)` `relative = TRUE` divides each fleet's observed and predicted index (and intervals) by the mean of its observed index and plots all fleets on one panel, colored by fleet, so that trends can be compared among indices (r4ss SSplotIndices subplot 9). Facet variables with one value for each fleet (e.g., the month of each survey) do not split the panel. `relative = TRUE` divides each fleet's observed and predicted index (and intervals) by the mean of its observed index and plots all fleets on one panel, colored by fleet, so that trends can be compared among indices (r4ss SSplotIndices subplot 9). Facet variables with one value for each fleet (e.g., the month of each survey) do not split the panel. --- R/plot_index.R | 73 ++++++++++++++++++++++++++++++-- man/plot_index.Rd | 16 +++++++ tests/testthat/test-plot_index.R | 44 +++++++++++++++++++ 3 files changed, 129 insertions(+), 4 deletions(-) diff --git a/R/plot_index.R b/R/plot_index.R index 693d16e0..ac8b7497 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -19,6 +19,14 @@ #' intervals are symmetric. #' #' Default: FALSE +#' @param relative Logical. TRUE/FALSE; divide each fleet's observed and +#' predicted index (and intervals) by the mean of its observed index, and plot +#' all fleets on one panel, colored by fleet. Indices with conflicting trends +#' (data conflicts) are a common cause of poor fits and are easiest to see +#' this way (`r4ss::SSplotIndices()` subplot 9). Any `group` other than +#' "fleet" is added to `facet`. +#' +#' 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) #' @@ -56,6 +64,13 @@ #' log = TRUE, #' interactive = FALSE #' ) +# +#' # All indices relative to their means on one panel +#' plot_index( +#' dat = stockplotr::example_data, +#' relative = TRUE, +#' interactive = FALSE +#' ) #' plot_index <- function( dat, @@ -68,6 +83,7 @@ plot_index <- function( focus = NULL, interval = NULL, log = FALSE, + relative = FALSE, make_rda = FALSE, figures_dir = getwd(), ... @@ -81,7 +97,15 @@ plot_index <- function( glue::glue("Estimated Index ({unit_label})") ) - if (is.null(facet) || "fleet" %notin% facet) { + if (relative) { + # All fleets on one panel, distinguished by color + if (!is.null(group) && group != "fleet") { + facet <- c(facet, group) + } + group <- "fleet" + facet <- setdiff(facet, "fleet") + if (length(facet) == 0) facet <- NULL + } else if (is.null(facet) || "fleet" %notin% facet) { facet <- c(facet, "fleet") } # Filter data @@ -105,6 +129,9 @@ plot_index <- function( if (!is.null(interval)) { prepared_data <- index_observed_interval(prepared_data, interval) } + if (relative) { + prepared_data <- index_relative(prepared_data) + } # Subset data if focus if (!is.null(focus)) { @@ -121,14 +148,24 @@ plot_index <- function( group <- processed_data[[2]] facet <- processed_data[[3]] - # move label to grouping and set grouping into facet - if (!is.null(group)) { + # move label to grouping and set grouping into facet (for relative + # indices, the group is shown by color on one panel) + if (!is.null(group) && !relative) { 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 + # Relative indices are on one panel; columns with one value for each fleet + # (e.g., the month of each survey) are shown by the fleet colors already + if (relative && !is.null(facet) && "fleet" %in% colnames(prepared_data)) { + nested <- vapply(facet, function(v) { + all(tapply(prepared_data[[v]], prepared_data$fleet, function(x) length(unique(x)) == 1)) + }, logical(1)) + facet <- facet[!nested] + if (length(facet) == 0) facet <- NULL + } # transform prep data so group_var = label # prepared_data <- prepared_data |> # dplyr::mutate( @@ -173,12 +210,23 @@ plot_index <- function( ) plt$layers <- c(list(intervals), plt$layers) } + if (relative) { + plt$layers <- c( + list(ggplot2::geom_hline(yintercept = 1, color = "gray70")), + plt$layers + ) + plt <- plt + ggplot2::labs(y = "Index relative to its mean") + # The model legend only has an entry for each model + if (length(unique(prepared_data$model)) == 1) { + plt <- plt + ggplot2::guides(shape = "none", linetype = "none") + } + } 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)") + ) + ggplot2::labs(y = paste(if (relative) "Index relative to its mean" else "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)) { @@ -256,3 +304,20 @@ index_observed_interval <- function(dat, interval) { dat$estimate_upper <- ifelse(observed, dat$estimate * exp(z * sigma), NA_real_) dat } + +# Divide the observed and predicted index (and intervals) of each fleet and +# model by the mean of its observed index +index_relative <- function(dat) { + keys <- intersect(c("model", "fleet"), names(dat)) + dat |> + dplyr::group_by(dplyr::across(dplyr::all_of(keys))) |> + dplyr::mutate( + index_mean = mean(.data$estimate[grepl("index_observed", .data$label)], na.rm = TRUE), + dplyr::across( + dplyr::any_of(c("estimate", "estimate_lower", "estimate_upper")), + ~ .x / .data$index_mean + ) + ) |> + dplyr::ungroup() |> + dplyr::select(-"index_mean") +} diff --git a/man/plot_index.Rd b/man/plot_index.Rd index bf069772..6ff15b28 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -14,6 +14,7 @@ plot_index( focus = NULL, interval = NULL, log = FALSE, + relative = FALSE, make_rda = FALSE, figures_dir = getwd(), ... @@ -72,6 +73,15 @@ intervals are symmetric. Default: FALSE} +\item{relative}{Logical. TRUE/FALSE; divide each fleet's observed and +predicted index (and intervals) by the mean of its observed index, and plot +all fleets on one panel, colored by fleet. Indices with conflicting trends +(data conflicts) are a common cause of poor fits and are easiest to see +this way (\code{r4ss::SSplotIndices()} subplot 9). Any \code{group} other than +"fleet" is added to \code{facet}. + +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". @@ -120,6 +130,12 @@ plot_index( log = TRUE, interactive = FALSE ) +# All indices relative to their means on one panel +plot_index( + dat = stockplotr::example_data, + relative = TRUE, + interactive = FALSE +) } \seealso{ diff --git a/tests/testthat/test-plot_index.R b/tests/testthat/test-plot_index.R index eefee85a..f0c1467f 100644 --- a/tests/testthat/test-plot_index.R +++ b/tests/testthat/test-plot_index.R @@ -89,3 +89,47 @@ test_that("plot_index intervals work for FIMS output", { 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))) }) + +test_that("plot_index can plot all indices relative to their means", { + p <- plot_index(stockplotr::example_data, relative = TRUE, interval = 0.95, interactive = FALSE) + b <- ggplot2::ggplot_build(p) + # One panel, with fleets shown by color + expect_equal(length(unique(b$layout$layout$PANEL)), 1) + expect_equal(p$labels$y, "Index relative to its mean") + expect_s3_class(p$layers[[1]]$geom, "GeomHline") + + # Observed values average one for each fleet, and the predicted index and + # intervals are scaled by the same mean + pts <- p$layers[[3]]$data + obs <- pts[pts$label == "index_observed", ] + expect_equal(as.vector(tapply(obs$estimate, obs$fleet, mean)), c(1, 1)) + raw <- stockplotr::example_data[stockplotr::example_data$label %in% c("index_observed", "index_predicted"), ] + tri <- raw[raw$fleet == "Triennial", ] + tri_mean <- mean(tri$estimate[tri$label == "index_observed"]) + pred <- p$layers[[4]]$data + expect_equal( + sort(pred$estimate[pred$fleet == "Triennial"]), + sort(tri$estimate[tri$label == "index_predicted"] / tri_mean) + ) + # Intervals keep their width relative to each observation + ints <- p$layers[[2]]$data + ints0 <- plot_index(stockplotr::example_data, interval = 0.95, interactive = FALSE)$layers[[1]]$data + ints <- ints[order(ints$fleet, ints$year), ] + ints0 <- ints0[order(ints0$fleet, ints0$year), ] + expect_equal(ints$estimate_lower / ints$estimate, ints0$estimate_lower / ints0$estimate) + expect_no_error(b) + + # Survey months that are nested within fleets do not split the panel + p <- suppressMessages(plot_index(stockplotr::pollock, relative = TRUE, interactive = FALSE)) + expect_equal(length(unique(ggplot2::ggplot_build(p)$layout$layout$PANEL)), 1) + + # Works on a log scale, with one fleet, and with FIMS output + p <- plot_index(stockplotr::example_data, relative = TRUE, log = TRUE, interactive = FALSE) + expect_equal(p$labels$y, "Index relative to its mean (log scale)") + expect_no_error(ggplot2::ggplot_build(p)) + expect_no_error(ggplot2::ggplot_build(suppressMessages(plot_index(stockplotr::ss3_simple, relative = TRUE, interactive = FALSE)))) + 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, relative = TRUE, interactive = FALSE)) + expect_equal(mean(p$layers[[2]]$data$estimate[p$layers[[2]]$data$label == "index_observed"]), 1) +}) From b663b7d4ddd0fca883f7530cea37a1269035f3a0 Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 6 Nov 2026 11:36:21 -0800 Subject: [PATCH 05/10] Update plot_index.R make sure this is consistent with other outstanding PR --- R/plot_index.R | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/R/plot_index.R b/R/plot_index.R index ac8b7497..d940a2d9 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 @@ -64,7 +64,6 @@ #' log = TRUE, #' interactive = FALSE #' ) -# #' # All indices relative to their means on one panel #' plot_index( #' dat = stockplotr::example_data, @@ -81,7 +80,7 @@ plot_index <- function( interactive = TRUE, module = NULL, focus = NULL, - interval = NULL, + interval = 0.95, log = FALSE, relative = FALSE, make_rda = FALSE, @@ -259,7 +258,7 @@ plot_index <- function( cpue.min, cpue.max, cpue.units, - dir=figures_dir + dir = figures_dir ) # Add key quantities to captions/alt text @@ -269,7 +268,7 @@ plot_index <- function( cpue.min, cpue.max, cpue.units, - dir=figures_dir + dir = figures_dir ) create_rda( From cef8f2f36576d929b9da0a2ba635a2fa396df902 Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 6 Nov 2026 11:46:12 -0800 Subject: [PATCH 06/10] 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 6ff15b28..6cf3714d 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, relative = FALSE, make_rda = FALSE, @@ -62,8 +62,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 5efd94d1f38535e18f71f4c91fdb8e7a3e648b8d Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 17:05:27 -0800 Subject: [PATCH 07/10] 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 66537498ab817ea534a8a46f6574a0e1f3492c6c Mon Sep 17 00:00:00 2001 From: ericward-noaa Date: Fri, 6 Nov 2026 18:40:25 -0800 Subject: [PATCH 08/10] update per changing defaults --- 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 64312d96..4585d736 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 ccc2ded5..ed1bc3cc 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -63,9 +63,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 f0c1467f..2ae88b18 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 b39f2e30a483ffdfbbb0423a707b8a9e1883945c Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 13 Nov 2026 13:05:24 -0800 Subject: [PATCH 09/10] Update R/plot_index.R Co-authored-by: Sam (Schiano) Bredeck <125507018+Schiano-NOAA@users.noreply.github.com> --- R/plot_index.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/plot_index.R b/R/plot_index.R index f1ddc385..e788f93b 100644 --- a/R/plot_index.R +++ b/R/plot_index.R @@ -8,7 +8,7 @@ #' Options: May include, but is not limited to: "fleet", "year", "area", NULL #' @param unit_label String. Index of abundance/CPUE units #' -#' Default: "" +#' 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. From cf5cd55c7055b6ce9b2b9924ea66a714ee4593a8 Mon Sep 17 00:00:00 2001 From: Eric Ward <5046884+ericward-noaa@users.noreply.github.com> Date: Fri, 13 Nov 2026 13:06:19 -0800 Subject: [PATCH 10/10] devtools::doc --- man/plot_index.Rd | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/man/plot_index.Rd b/man/plot_index.Rd index d2490043..20003694 100644 --- a/man/plot_index.Rd +++ b/man/plot_index.Rd @@ -27,7 +27,9 @@ 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. @@ -36,9 +38,11 @@ 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, but 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 + +Options: May include, but is not limited to: "fleet", "year", "area", NULL} \item{interactive}{Logical. TRUE/FALSE; indicate whether the environment is interactive. @@ -55,14 +59,6 @@ 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} -\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} - -\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{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. @@ -113,6 +109,14 @@ standard error. Uncertainty labeled "se", "sd", or "log_sd" is used as \eqn{\sigma = \sqrt{\log(1 + CV^2)}}. Observations without an uncertainty have no interval. +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