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