From 90db368d190d05284b01431803934e25766c1bdb Mon Sep 17 00:00:00 2001 From: chantelwetzel-noaa Date: Wed, 2 Oct 2024 18:40:00 +0000 Subject: [PATCH] style and docs: run devtools::document() and styler::style_pkg() --- DESCRIPTION | 9 +-- R/check_profile_range.R | 86 ++++++++++++------------ R/get_jitter_quants.R | 7 +- R/get_param_values.R | 29 ++++----- R/get_retro_quants.R | 15 +++-- R/get_settings.R | 6 +- R/get_settings_profile.R | 9 ++- R/get_summary.R | 1 - R/jitter_wrapper.R | 3 +- R/plot_jitter.R | 16 ++--- R/plot_retro.R | 72 ++++++++++++--------- R/profile_plot.R | 26 +++++--- R/profile_wrapper.R | 1 - R/rerun_profile_vals.R | 12 ++-- R/retro_wrapper.R | 4 +- R/run_diagnostics.R | 3 +- R/run_profile.R | 38 +++++++---- R/run_retro.R | 11 ++-- man/get_settings_profile.Rd | 8 +-- man/rerun_profile_vals.Rd | 4 +- tests/testthat.R | 2 +- tests/testthat/test-runs.R | 126 ++++++++++++++++++++---------------- 22 files changed, 261 insertions(+), 227 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 016ceb1..c6602b5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,14 +1,14 @@ -Package: nwfscDiag Type: Package +Package: nwfscDiag Title: Generate Standard NWFSC Assessment Diagnostics Version: 1.1.2 Author: Chantel Wetzel Maintainer: Chantel Wetzel -Description: Package that can automates diagnositics for SS3 models by running jitters, retrospective, and profiles. +Description: Package that can automates diagnositics for SS3 models by + running jitters, retrospective, and profiles. License: GPL (>=3) URL: https://github.com/pfmc-assessments/nwfscDiag BugReports: https://github.com/pfmc-assessments/nwfscDiag/issues -LazyData: true Depends: R (>= 3.5) Imports: @@ -28,6 +28,7 @@ Suggests: testthat Remotes: github::r4ss/r4ss -RoxygenNote: 7.3.2 Encoding: UTF-8 +LazyData: true Roxygen: list(markdown = TRUE) +RoxygenNote: 7.3.2 diff --git a/R/check_profile_range.R b/R/check_profile_range.R index 1dc6b35..b55cd32 100644 --- a/R/check_profile_range.R +++ b/R/check_profile_range.R @@ -10,7 +10,7 @@ #' @export #' #' -check_profile_range <- function(mydir, model_settings){ +check_profile_range <- function(mydir, model_settings) { # Read in the base model rep <- r4ss::SS_output( file.path(mydir, model_settings$base_name), @@ -20,54 +20,54 @@ check_profile_range <- function(mydir, model_settings){ ) N <- nrow(model_settings$profile_details) - for (aa in 1:N){ + for (aa in 1:N) { profile_details <- model_settings[["profile_details"]][aa, ] para <- profile_details[, "parameters"] est <- rep$parameters[rep$parameters$Label == para, "Value"] - # Determine the parameter range - if (profile_details$param_space == "relative") { - range <- c( - est + profile_details$low, - est + profile_details$high - ) - } - if (profile_details$param_space == "multiplier") { - range <- c( - est - est * profile_details$low, - est + est * profile_details$high - ) - } - if (profile_details$param_space == "real") { - range <- c( - profile_details$low, - profile_details$high - ) - } - step_size <- profile_details$step_size + # Determine the parameter range + if (profile_details$param_space == "relative") { + range <- c( + est + profile_details$low, + est + profile_details$high + ) + } + if (profile_details$param_space == "multiplier") { + range <- c( + est - est * profile_details$low, + est + est * profile_details$high + ) + } + if (profile_details$param_space == "real") { + range <- c( + profile_details$low, + profile_details$high + ) + } + step_size <- profile_details$step_size - # Create parameter vect from base down and the base up - if (est != round_any(est, step_size, f = floor)) { - low <- rev(seq( - round_any(range[1], step_size, f = ceiling), - round_any(est, step_size, f = floor), step_size - )) - } else { - low <- rev(seq( - round_any(range[1], step_size, f = ceiling), - round_any(est, step_size, f = floor) - step_size, step_size - )) - } + # Create parameter vect from base down and the base up + if (est != round_any(est, step_size, f = floor)) { + low <- rev(seq( + round_any(range[1], step_size, f = ceiling), + round_any(est, step_size, f = floor), step_size + )) + } else { + low <- rev(seq( + round_any(range[1], step_size, f = ceiling), + round_any(est, step_size, f = floor) - step_size, step_size + )) + } - if (est != round_any(est, step_size, f = ceiling)) { - high <- c(est, seq(round_any(est, step_size, f = ceiling), range[2], step_size)) - } else { - high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size)) - } + if (est != round_any(est, step_size, f = ceiling)) { + high <- c(est, seq(round_any(est, step_size, f = ceiling), range[2], step_size)) + } else { + high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size)) + } - vec <- c(low, high) - cli::cli_inform( - "Profiling over {para} across values of {vec}." - ) + vec <- c(low, high) + cli::cli_inform( + "Profiling over {para} across values of {vec}." + ) } } diff --git a/R/get_jitter_quants.R b/R/get_jitter_quants.R index bb7d54d..ab628ff 100644 --- a/R/get_jitter_quants.R +++ b/R/get_jitter_quants.R @@ -58,10 +58,11 @@ get_jitter_quants <- function(mydir, model_settings, output) { "Through the jittering analysis performed here and ", "the estimation of likelihood profiles, ", "we are confident that the base model as presented represents the ", - "best fit to the data given the assumptions made."), + "best fit to the data given the assumptions made." + ), alt_caption = "Comparison of the negative log-likelihood across jitter runs", - label = c("jitter", "jitter-zoomed"), - filein = file.path("..", jitter_dir, c("jitter.png", "jitter_zoomed.png")) + label = c("jitter", "jitter-zoomed"), + filein = file.path("..", jitter_dir, c("jitter.png", "jitter_zoomed.png")) ), file = file.path(jitter_dir, "jitterfigures4doc.csv"), row.names = FALSE diff --git a/R/get_param_values.R b/R/get_param_values.R index d1d2e6b..6738f9e 100644 --- a/R/get_param_values.R +++ b/R/get_param_values.R @@ -14,7 +14,6 @@ #' @export get_param_values <- function(mydir, para = NULL, vec, summary) { - x <- summary n <- x[["n"]] endyr <- x[["endyrs"]][1] + 1 @@ -52,20 +51,20 @@ get_param_values <- function(mydir, para = NULL, vec, summary) { out <- t(out) colnames(out) <- vec - if(!is.null(para)) { - name <- para - if (para == "SR_LN(R0)") { - colnames(out) <- paste0("R0 ", vec) - } - if (para == "NatM_uniform_Fem_GP_1") { - colnames(out) <- paste0("M_f ", vec) - } - if (para == "NatM_uniform_Mal_GP_1") { - colnames(out) <- paste0("M_m ", vec) - } - if (para == "SR_BH_steep") { - colnames(out) <- paste0("h ", vec) - } + if (!is.null(para)) { + name <- para + if (para == "SR_LN(R0)") { + colnames(out) <- paste0("R0 ", vec) + } + if (para == "NatM_uniform_Fem_GP_1") { + colnames(out) <- paste0("M_f ", vec) + } + if (para == "NatM_uniform_Mal_GP_1") { + colnames(out) <- paste0("M_m ", vec) + } + if (para == "SR_BH_steep") { + colnames(out) <- paste0("h ", vec) + } } utils::write.csv(x = out, file = file.path(mydir, paste0(name, "_quant_table.csv")), row.names = TRUE) diff --git a/R/get_retro_quants.R b/R/get_retro_quants.R index 0378462..46eb1da 100644 --- a/R/get_retro_quants.R +++ b/R/get_retro_quants.R @@ -28,7 +28,7 @@ #' #' @export -get_retro_quants <- function(mydir, model_settings, output) { +get_retro_quants <- function(mydir, model_settings, output) { retro_dir <- output[["plotdir"]] endyrvec <- output[["endyrvec"]] retroSummary <- output[["retroSummary"]] @@ -38,7 +38,7 @@ get_retro_quants <- function(mydir, model_settings, output) { get_param_values( mydir = retro_dir, para = "retro", - vec = c("Base Model", paste0("Retro -", 1:(length(endyrvec)-1))), + vec = c("Base Model", paste0("Retro -", 1:(length(endyrvec) - 1))), summary = retroSummary ) @@ -53,11 +53,12 @@ get_retro_quants <- function(mydir, model_settings, output) { "recalculated for each peel given the removal of another year of data.", "See Table \\ref{tab:RetroMohnsrho} for other derivations of Mohn's rho." ), - alt_caption = sprintf("Each successive peel of data led to a Mohn's rho of %s for %s.", - lapply(c("SSB", "Bratio"), function(x) { - knitr::combine_words(sprintf("%.2f", (rhosall[rownames(rhosall) == x, ]))) - }), - c("SSB", "fraction unfished") + alt_caption = sprintf( + "Each successive peel of data led to a Mohn's rho of %s for %s.", + lapply(c("SSB", "Bratio"), function(x) { + knitr::combine_words(sprintf("%.2f", (rhosall[rownames(rhosall) == x, ]))) + }), + c("SSB", "fraction unfished") ), label = c("RetroSsb", "RetroFractionunfished"), filein = file.path("..", retro_dir, c("compare2_spawnbio_uncertainty.png", "compare4_Bratio_uncertainty.png")) diff --git a/R/get_settings.R b/R/get_settings.R index c4926ce..13a2b94 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -17,7 +17,6 @@ #' get_settings(list("Njitter" = 10)) #' get_settings <- function(settings = NULL, verbose = FALSE) { - if (is.vector(settings)) settings <- as.list(settings) Settings_all <- list( @@ -91,10 +90,11 @@ get_settings <- function(settings = NULL, verbose = FALSE) { } if ("profile" %in% Settings_all[["run"]]) { - if (Settings_all[["verbose"]]){ + if (Settings_all[["verbose"]]) { check_profile_range( mydir = mydir, - model_settings = Settings_all) + model_settings = Settings_all + ) } } diff --git a/R/get_settings_profile.R b/R/get_settings_profile.R index 2f8196e..6b01c59 100644 --- a/R/get_settings_profile.R +++ b/R/get_settings_profile.R @@ -61,10 +61,10 @@ #' # Define each parameter in real space #' get_settings_profile( #' parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), -#' low = c(0.02, 0.25, 8), -#' high = c(0.07, 1.0, 11), -#' step_size = c(0.005, 0.05, 0.25), -#' param_space = c('real', 'real', 'real') +#' low = c(0.02, 0.25, 8), +#' high = c(0.07, 1.0, 11), +#' step_size = c(0.005, 0.05, 0.25), +#' param_space = c("real", "real", "real") #' ) #' #' # Example 2: Run a profile for natural mortality one with the prior likelihood and one without @@ -83,7 +83,6 @@ get_settings_profile <- function(parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_ step_size = c(0.01, 0.05, 0.25), param_space = c("multiplier", "real", "relative"), use_prior_like = lifecycle::deprecated()) { - if (length(parameters) != length(low) | length(parameters) != length(high) | length(parameters) != length(step_size) | diff --git a/R/get_summary.R b/R/get_summary.R index e80847c..72e5377 100644 --- a/R/get_summary.R +++ b/R/get_summary.R @@ -15,7 +15,6 @@ #' @export get_summary <- function(mydir, para, vec, profilemodels, profilesummary) { - # Need to identify a way to determine if a model estimates male growth parameters as offsets from females # get output diff --git a/R/jitter_wrapper.R b/R/jitter_wrapper.R index f114113..562e1a2 100644 --- a/R/jitter_wrapper.R +++ b/R/jitter_wrapper.R @@ -20,7 +20,8 @@ jitter_wrapper <- function(mydir, model_settings) { output <- run_jitter( mydir = mydir, - model_settings = model_settings) + model_settings = model_settings + ) plot_jitter( mydir = mydir, model_settings = model_settings, diff --git a/R/plot_jitter.R b/R/plot_jitter.R index aef6e20..4312907 100644 --- a/R/plot_jitter.R +++ b/R/plot_jitter.R @@ -26,8 +26,8 @@ plot_jitter <- function(mydir, model_settings, output) { pngfun(wd = jitter_dir, file = "jitter.png", h = 12, w = 9) on.exit(grDevices::dev.off(), add = TRUE) plot(keys, like - est, - ylim = c(ymin, ymax), cex.axis = 1.25, cex.lab = 1.25, - ylab = ylab, xlab = xlab + ylim = c(ymin, ymax), cex.axis = 1.25, cex.lab = 1.25, + ylab = ylab, xlab = xlab ) graphics::abline(h = 0, col = "darkgrey", lwd = 2) find <- which(est == like) @@ -43,16 +43,16 @@ plot_jitter <- function(mydir, model_settings, output) { ) } graphics::legend("topleft", - legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), - bty = "n", pch = 16, col = c("green3", "blue", "red") + legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), + bty = "n", pch = 16, col = c("green3", "blue", "red") ) if (ymax > 100) { pngfun(wd = jitter_dir, file = "jitter_zoomed.png", h = 12, w = 9) on.exit(grDevices::dev.off(), add = TRUE) plot(keys, like - est, - ylim = c(ymin, 100), cex.axis = 1.25, cex.lab = 1.25, - ylab = ylab, xlab = xlab + ylim = c(ymin, 100), cex.axis = 1.25, cex.lab = 1.25, + ylab = ylab, xlab = xlab ) graphics::abline(h = 0, col = "darkgrey", lwd = 2) find <- which(est == like) @@ -68,8 +68,8 @@ plot_jitter <- function(mydir, model_settings, output) { ) } graphics::legend("topleft", - legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), - bty = "n", pch = 16, col = c("green3", "blue", "red") + legend = c("Base Model Likelihood", "Higher Likelihood", "Lower Likelihood"), + bty = "n", pch = 16, col = c("green3", "blue", "red") ) } } diff --git a/R/plot_retro.R b/R/plot_retro.R index a7960d3..b960337 100644 --- a/R/plot_retro.R +++ b/R/plot_retro.R @@ -17,7 +17,7 @@ #' #' @export -plot_retro <- function(mydir, model_settings, output) { +plot_retro <- function(mydir, model_settings, output) { retro_dir <- output[["plotdir"]] endyrvec <- output[["endyrvec"]] retroSummary <- output[["retroSummary"]] @@ -30,9 +30,10 @@ plot_retro <- function(mydir, model_settings, output) { endyrvec = endyrvec, legendlabels = c( "Base Model", - sprintf("Data %.0f year%s", - model_settings[["retro_yrs"]], - ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s") + sprintf( + "Data %.0f year%s", + model_settings[["retro_yrs"]], + ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s") ) ), btarg = model_settings[["btarg"]], @@ -64,13 +65,15 @@ plot_retro <- function(mydir, model_settings, output) { function(x) { c( "Base Model", - sprintf("Data %.0f year%s (Revised Mohn's rho %.2f)", - model_settings[["retro_yrs"]], - ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s"), - rhosall[rownames(rhosall) == x, ] + sprintf( + "Data %.0f year%s (Revised Mohn's rho %.2f)", + model_settings[["retro_yrs"]], + ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s"), + rhosall[rownames(rhosall) == x, ] ) ) - }) + } + ) ) r4ss::SSplotComparisons( @@ -78,9 +81,10 @@ plot_retro <- function(mydir, model_settings, output) { endyrvec = endyrvec, legendlabels = c( "Base Model", - sprintf("Data %.0f year%s", - model_settings[["retro_yrs"]], - ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s") + sprintf( + "Data %.0f year%s", + model_settings[["retro_yrs"]], + ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s") ) ), btarg = model_settings[["btarg"]], @@ -113,20 +117,22 @@ plot_retro <- function(mydir, model_settings, output) { function(x) { c( "Base Model", - sprintf("Data %.0f year%s (Revised Mohn's rho %.2f)", - model_settings[["retro_yrs"]], - ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s"), - rhosall[rownames(rhosall) == x, ] + sprintf( + "Data %.0f year%s (Revised Mohn's rho %.2f)", + model_settings[["retro_yrs"]], + ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s"), + rhosall[rownames(rhosall) == x, ] ) ) - }) + } + ) ) label <- ifelse( retroSummary[["SpawnOutputUnits"]][1] == "biomass", "Spawning Biomass", "Spawning Output" - ) + ) n <- ncol(retroSummary[["SpawnBio"]]) - 2 years <- retroSummary[["startyrs"]][1]:retroSummary[["endyrs"]][1] + 1 denom <- paste0("model", 1:n) @@ -169,10 +175,11 @@ plot_retro <- function(mydir, model_settings, output) { col_names <- c("Yr", "Reference_Point", paste0("per_diff_model", 1:n)) df <- rbind( - sb[ , colnames(sb) %in% col_names], + sb[, colnames(sb) %in% col_names], bratio[, colnames(bratio) %in% col_names], - f[ , colnames(f) %in% col_names], - rec[, colnames(rec) %in% col_names]) |> + f[, colnames(f) %in% col_names], + rec[, colnames(rec) %in% col_names] + ) |> tidyr::pivot_longer( cols = starts_with("per_diff_model"), names_to = "model", @@ -181,18 +188,18 @@ plot_retro <- function(mydir, model_settings, output) { df_out <- NULL y <- years - for (a in 1:n){ + for (a in 1:n) { col_name <- paste0("per_diff_model", a) df_out <- rbind(df_out, df[df[["Yr"]] %in% y & df[["model"]] %in% col_name, ]) - if (a == 1){ + if (a == 1) { df_out[["model"]][df_out[["model"]] == col_name] <- "Base Model" } else { - df_out[["model"]][df_out[["model"]] == col_name] <- paste0("Retro -", a-1) + df_out[["model"]][df_out[["model"]] == col_name] <- paste0("Retro -", a - 1) } y <- y - 1 } colnames(df_out)[colnames(df_out) == "model"] <- "Run" - leg_order <- c("Base Model", paste0("Retro -", 1:(length(endyrvec) -1))) + leg_order <- c("Base Model", paste0("Retro -", 1:(length(endyrvec) - 1))) df_out[["Run"]] <- factor(df_out[["Run"]], levels = leg_order) xrange <- c(ifelse(min(df_out[["Yr"]]) < 1980, 1980, min(df_out[["Yr"]])), max(df_out[["Yr"]])) find <- df_out |> @@ -205,27 +212,28 @@ plot_retro <- function(mydir, model_settings, output) { ggplot2::ggplot(df_out, ggplot2::aes(x = Yr, y = diff, col = Run)) + ggplot2::geom_line() + ggplot2::geom_point() + - #ylim(yrange) + - ggplot2::scale_x_continuous(limits = xrange, expand = c(0,0)) + + # ylim(yrange) + + ggplot2::scale_x_continuous(limits = xrange, expand = c(0, 0)) + ggplot2::ylab("% Differece from Base Model") + ggplot2::xlab("Year") + ggplot2::scale_colour_viridis_d() + ggplot2::theme_bw(base_size = 15) + - ggplot2::facet_wrap("Reference_Point", scales = 'free_y') + ggplot2::facet_wrap("Reference_Point", scales = "free_y") ggplot2::ggsave(filename = file.path(retro_dir, "retro_percent_difference_4_panel.png"), width = 10, height = 12) - sub_out <- df_out[df_out[["Reference_Point"]] != "Recruits",] + sub_out <- df_out[df_out[["Reference_Point"]] != "Recruits", ] find <- sub_out |> dplyr::filter(Yr >= xrange[1]) |> dplyr::summarize( - bound = abs(max(diff))) + bound = abs(max(diff)) + ) yrange <- c(-1 * find[["bound"]] - 5, find[["bound"]] + 5) ggplot2::ggplot(sub_out, ggplot2::aes(x = Yr, y = diff, col = Run)) + ggplot2::geom_line() + ggplot2::geom_point() + ggplot2::ylim(yrange) + - ggplot2::scale_x_continuous(limits = xrange, expand = c(0,0)) + + ggplot2::scale_x_continuous(limits = xrange, expand = c(0, 0)) + ggplot2::ylab("% Differece from Base Model") + ggplot2::xlab("Year") + ggplot2::scale_colour_viridis_d() + @@ -233,7 +241,7 @@ plot_retro <- function(mydir, model_settings, output) { ggplot2::facet_wrap("Reference_Point", ncol = 1, nrow = 3) ggplot2::ggsave(filename = file.path(retro_dir, "retro_percent_difference_3_panel.png"), width = 10, height = 12) - endyrvec_long <- (min(endyrvec)-5):max(endyrvec) + endyrvec_long <- (min(endyrvec) - 5):max(endyrvec) pngfun(wd = retro_dir, file = "recruitment_retrospective_squid.png", h = 7, w = 7) r4ss::SSplotRetroRecruits( retroSummary = retroSummary, diff --git a/R/profile_plot.R b/R/profile_plot.R index 22f7c13..bfa7232 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -19,7 +19,6 @@ #' @seealso [profile_wrapper] and [rerun_profile_vals] call `plot_profile`. plot_profile <- function(mydir, rep, para, profilesummary) { - label <- ifelse(para == "SR_LN(R0)", expression(log(italic(R)[0])), ifelse(para %in% c("NatM_p_1_Fem_GP_1", "NatM_uniform_Fem_GP_1"), "Natural Mortality (female)", ifelse(para %in% c("NatM_p_1_Mal_GP_1", "NatM_uniform_Mal_GP_1"), "Natural Mortality (male)", @@ -46,7 +45,7 @@ plot_profile <- function(mydir, rep, para, profilesummary) { like_comp <- unique(profilesummary[["likelihoods_by_fleet"]][["Label"]][ c( -grep("_lambda", profilesummary[["likelihoods_by_fleet"]][["Label"]]), - -grep("_N_use", profilesummary[["likelihoods_by_fleet"]][["Label"]]), + -grep("_N_use", profilesummary[["likelihoods_by_fleet"]][["Label"]]), -grep("_N_skip", profilesummary[["likelihoods_by_fleet"]][["Label"]]) ) ]) @@ -144,7 +143,8 @@ plot_profile <- function(mydir, rep, para, profilesummary) { like <- as.numeric(profilesummary[["likelihoods"]][profilesummary[["likelihoods"]][["Label"]] == "TOTAL", n] - ifelse(starter[["prior_like"]] == 0, profilesummary[["likelihoods"]][profilesummary[["likelihoods"]][["Label"]] == "Parm_priors", n], - 0) - + 0 + ) - rep[["likelihoods_used"]][1, 1]) ylike <- c(min(like) + ifelse(min(like) != 0, -0.5, 0), max(like)) @@ -172,23 +172,29 @@ plot_profile <- function(mydir, rep, para, profilesummary) { plot(x, depl, type = "l", lwd = 2, xlab = label, ylab = expression("Fraction Unfished"[final]), ylim = c(0, 1.2)) points(est, depl_est, pch = 21, col = "black", bg = "blue", cex = 1.5) abline(h = c(btarg, thresh), lty = c(2, 2), col = c("darkgreen", "red")) - if(btarg > 0){ + if (btarg > 0) { graphics::legend("bottomright", - legend = c("Management target", "Minimum stock size threshold"), - lty = 2, col = c("darkgreen", "red"), bty = "n" + legend = c("Management target", "Minimum stock size threshold"), + lty = 2, col = c("darkgreen", "red"), bty = "n" ) } # parameter vs. SB0 - plot(x, sb0, type = "l", lwd = 2, xlab = label, + plot(x, sb0, + type = "l", lwd = 2, xlab = label, ylab = ifelse(profilesummary[["SpawnOutputUnits"]][1] == "numbers", - expression(SO[0]), expression(SB[0])), ylim = c(0, max(sb0))) + expression(SO[0]), expression(SB[0]) + ), ylim = c(0, max(sb0)) + ) points(est, sb0_est, pch = 21, col = "black", bg = "blue", cex = 1.5) # parameter vs. SBfinal - plot(x, sbf, type = "l", lwd = 2, xlab = label, + plot(x, sbf, + type = "l", lwd = 2, xlab = label, ylab = ifelse(profilesummary[["SpawnOutputUnits"]][1] == "numbers", - expression(SO[final]), expression(SB[final])), ylim = c(0, max(sbf))) + expression(SO[final]), expression(SB[final]) + ), ylim = c(0, max(sbf)) + ) points(est, sbf_est, pch = 21, col = "black", bg = "blue", cex = 1.5) # Create the sb and depl trajectories plot diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 0b3ea26..0109ce3 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -32,7 +32,6 @@ #' @export profile_wrapper <- function(mydir, model_settings) { - N <- nrow(model_settings$profile_details) for (aa in 1:N) { diff --git a/R/rerun_profile_vals.R b/R/rerun_profile_vals.R index fdb2433..0718c59 100644 --- a/R/rerun_profile_vals.R +++ b/R/rerun_profile_vals.R @@ -42,7 +42,7 @@ #' run = c("profile"), #' profile_details = get_settings_profile( #' parameters = c("NatM_uniform_Fem_GP_1"), -#' low = c(0.20), +#' low = c(0.20), #' high = c(0.25), #' step_size = c(0.02), #' param_space = c("multiplier") @@ -56,11 +56,11 @@ #' rerun_profile_vals( #' mydir = temp_profile_dir, #' model_settings = model_settings, -#' para_name = "NatM_uniform_Fem_GP_1", +#' para_name = "NatM_uniform_Fem_GP_1", #' run_num = c(1), #' data_file_nm = "data.ss" #' ) -#'} +#' } rerun_profile_vals <- function(mydir, para_name, model_settings, @@ -79,7 +79,7 @@ rerun_profile_vals <- function(mydir, } para <- para_name - profile_dir <- paste(mydir, "profile", para_name, sep = "_") + profile_dir <- paste(mydir, "profile", para_name, sep = "_") temp_dir <- file.path(profile_dir, "temp") dir.create(temp_dir, showWarnings = FALSE) @@ -99,7 +99,7 @@ rerun_profile_vals <- function(mydir, # Use the SS_parlines function to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( - ctlfile = model_settings$oldctlfile, + ctlfile = model_settings$oldctlfile, dir = temp_dir, verbose = FALSE, active = FALSE @@ -155,7 +155,7 @@ rerun_profile_vals <- function(mydir, if (like >= like_check[i]) { for (ii in 1:5) { starter <- r4ss::SS_readstarter(file = file.path(temp_dir, "starter.ss")) - if (ii == 1){ + if (ii == 1) { starter$jitter_fraction <- 0.01 } else { starter$jitter_fraction <- add + starter$jitter_fraction diff --git a/R/retro_wrapper.R b/R/retro_wrapper.R index 882b434..edb489f 100644 --- a/R/retro_wrapper.R +++ b/R/retro_wrapper.R @@ -40,7 +40,7 @@ #' #' @export -retro_wrapper <- function(mydir, model_settings) { +retro_wrapper <- function(mydir, model_settings) { output <- run_retro( mydir = mydir, model_settings = model_settings @@ -55,5 +55,5 @@ retro_wrapper <- function(mydir, model_settings) { model_settings = model_settings, output = output ) - cli::cli_inform("Finished retrospective runs.") + cli::cli_inform("Finished retrospective runs.") } diff --git a/R/run_diagnostics.R b/R/run_diagnostics.R index dd2d4fe..758d051 100644 --- a/R/run_diagnostics.R +++ b/R/run_diagnostics.R @@ -12,10 +12,9 @@ #' @export run_diagnostics <- function(mydir, model_settings) { - exe <- r4ss::check_exe(exe = model_settings$exe, dir = file.path(mydir, model_settings[["base_name"]]))[["exe"]] model_settings[["exe"]] <- exe - '%>%' <- magrittr::'%>%' + "%>%" <- magrittr::"%>%" # Check for Report file model_dir <- file.path(mydir, paste0(model_settings[["base_name"]])) diff --git a/R/run_profile.R b/R/run_profile.R index d0916a8..4bc8550 100644 --- a/R/run_profile.R +++ b/R/run_profile.R @@ -35,7 +35,6 @@ #' @export run_profile <- function(mydir, model_settings, para) { - # Create a profile folder with the same naming structure as the base model # Add a label to show if prior was used or not profile_dir <- file.path(mydir, paste0(model_settings[["base_name"]], "_profile_", para)) @@ -77,7 +76,7 @@ run_profile <- function(mydir, model_settings, para) { # Use the SS_parlines function to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( - ctlfile = model_settings[["oldctlfile"]], + ctlfile = model_settings[["oldctlfile"]], dir = profile_dir, verbose = FALSE, version = model_settings[["version"]], @@ -90,8 +89,10 @@ run_profile <- function(mydir, model_settings, para) { } # Copy oldctlfile to newctlfile before modifying it - file.copy(file.path(profile_dir, model_settings[["oldctlfile"]]), - file.path(profile_dir, model_settings[["newctlfile"]])) + file.copy( + file.path(profile_dir, model_settings[["oldctlfile"]]), + file.path(profile_dir, model_settings[["newctlfile"]]) + ) # Change the control file name in the starter file starter <- r4ss::SS_readstarter(file = file.path(profile_dir, "starter.ss")) @@ -153,28 +154,37 @@ run_profile <- function(mydir, model_settings, para) { # backup original control.ss_new file for use in second half of profile file.copy(file.path(profile_dir, model_settings[["oldctlfile"]]), - file.path(profile_dir, "backup_oldctlfile.ss"), - overwrite = model_settings$overwrite) + file.path(profile_dir, "backup_oldctlfile.ss"), + overwrite = model_settings$overwrite + ) # backup original par file for use in second half of profile # if usepar = TRUE file.copy(file.path(profile_dir, "ss.par"), - file.path(profile_dir, "backup_ss.par"), - overwrite = model_settings[["overwrite"]]) + file.path(profile_dir, "backup_ss.par"), + overwrite = model_settings[["overwrite"]] + ) # loop over down, then up for (iprofile in 1:2) { - whichruns <- which(vec %in% if(iprofile == 1){low} else {high}) + whichruns <- which(vec %in% if (iprofile == 1) { + low + } else { + high + }) if (!is.null(model_settings[["whichruns"]])) { whichruns <- intersect(model_settings[["whichruns"]], whichruns) } if (iprofile == 2) { # copy backup back to use in second half of profile - file.copy(file.path(profile_dir, "backup_oldctlfile.ss"), - file.path(profile_dir, model_settings[["oldctlfile"]])) + file.copy( + file.path(profile_dir, "backup_oldctlfile.ss"), + file.path(profile_dir, model_settings[["oldctlfile"]]) + ) # copy backup back to use in second half of profile file.copy(file.path(profile_dir, "backup_ss.par"), - file.path(profile_dir, "ss.par"), - overwrite = model_settings[["overwrite"]]) + file.path(profile_dir, "ss.par"), + overwrite = model_settings[["overwrite"]] + ) } profile <- r4ss::profile( dir = profile_dir, @@ -204,7 +214,7 @@ run_profile <- function(mydir, model_settings, para) { profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) profilesummary <- r4ss::SSsummarize(biglist = profilemodels) - if(!is.null(model_settings[["btarg"]])){ + if (!is.null(model_settings[["btarg"]])) { profilesummary[["btarg"]] <- model_settings[["btarg"]] profilesummary[["minbthresh"]] <- model_settings[["minbthresh"]] } diff --git a/R/run_retro.R b/R/run_retro.R index 9f8a475..a955b2b 100644 --- a/R/run_retro.R +++ b/R/run_retro.R @@ -40,9 +40,8 @@ #' #' @export -run_retro <- function(mydir, model_settings) { - - if(!file.exists(file.path(mydir, model_settings[["base_name"]], "Report.sso"))) { +run_retro <- function(mydir, model_settings) { + if (!file.exists(file.path(mydir, model_settings[["base_name"]], "Report.sso"))) { base <- model_settings[["base_name"]] cli::cli_abort("There is no Report.sso file in the base model directory {file.path(mydir, base}") } @@ -50,7 +49,7 @@ run_retro <- function(mydir, model_settings) { # Create a jitter folder with the same naming structure as the base model retro_dir <- file.path(mydir, paste0(model_settings[["base_name"]], "_retro_", length(model_settings[["retro_yrs"]]), "_yr_peel")) dir.create(retro_dir, showWarnings = FALSE) - all_files = list.files(file.path(mydir, model_settings[["base_name"]])) + all_files <- list.files(file.path(mydir, model_settings[["base_name"]])) ignore <- file.copy( from = file.path(mydir, model_settings[["base_name"]], all_files), to = retro_dir, @@ -73,11 +72,11 @@ run_retro <- function(mydir, model_settings) { ignore <- file.remove(from = file.path(retro_dir, all_files)) runs <- list() - for(aa in 1:(length(model_settings[["retro_yrs"]]) + 1)) { + for (aa in 1:(length(model_settings[["retro_yrs"]]) + 1)) { if (aa == 1) { runs[[aa]] <- r4ss::SS_output(dir = file.path(mydir, model_settings[["base_name"]]), verbose = FALSE, printstats = FALSE) } else { - tmp = file.path(retro_dir, model_settings[["newsubdir"]], paste0("retro", model_settings[["retro_yrs"]][aa-1])) + tmp <- file.path(retro_dir, model_settings[["newsubdir"]], paste0("retro", model_settings[["retro_yrs"]][aa - 1])) runs[[aa]] <- r4ss::SS_output(dir = tmp, verbose = FALSE, printstats = FALSE) } } diff --git a/man/get_settings_profile.Rd b/man/get_settings_profile.Rd index 4b42811..12dc120 100644 --- a/man/get_settings_profile.Rd +++ b/man/get_settings_profile.Rd @@ -81,10 +81,10 @@ any of the options for specifying a parameter range for any parameter. # Define each parameter in real space get_settings_profile( parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), - low = c(0.02, 0.25, 8), - high = c(0.07, 1.0, 11), - step_size = c(0.005, 0.05, 0.25), - param_space = c('real', 'real', 'real') + low = c(0.02, 0.25, 8), + high = c(0.07, 1.0, 11), + step_size = c(0.005, 0.05, 0.25), + param_space = c("real", "real", "real") ) # Example 2: Run a profile for natural mortality one with the prior likelihood and one without diff --git a/man/rerun_profile_vals.Rd b/man/rerun_profile_vals.Rd index a35f2dc..d6552d5 100644 --- a/man/rerun_profile_vals.Rd +++ b/man/rerun_profile_vals.Rd @@ -55,7 +55,7 @@ model_settings <- get_settings( run = c("profile"), profile_details = get_settings_profile( parameters = c("NatM_uniform_Fem_GP_1"), - low = c(0.20), + low = c(0.20), high = c(0.25), step_size = c(0.02), param_space = c("multiplier") @@ -69,7 +69,7 @@ run_diagnostics( rerun_profile_vals( mydir = temp_profile_dir, model_settings = model_settings, - para_name = "NatM_uniform_Fem_GP_1", + para_name = "NatM_uniform_Fem_GP_1", run_num = c(1), data_file_nm = "data.ss" ) diff --git a/tests/testthat.R b/tests/testthat.R index 03c1bc6..b284012 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,4 +1,4 @@ library(testthat) library(nwfscDiag) -test_check("nwfscDiag") \ No newline at end of file +test_check("nwfscDiag") diff --git a/tests/testthat/test-runs.R b/tests/testthat/test-runs.R index 621e185..7ae6717 100644 --- a/tests/testthat/test-runs.R +++ b/tests/testthat/test-runs.R @@ -15,84 +15,96 @@ on.exit(unlink(tmp_path, recursive = TRUE)) test_path <- file.path(runs_path, "simple_small") skip_test <- TRUE -if(.Platform$OS.type == "windows") { - if (file.exists(file.path(test_path, "ss3.exe"))) { - skip_test <- FALSE - } +if (.Platform$OS.type == "windows") { + if (file.exists(file.path(test_path, "ss3.exe"))) { + skip_test <- FALSE + } } -if(.Platform$OS.type == "unix") { - if (file.exists(file.path(test_path, "ss3"))) { - skip_test <- FALSE - } +if (.Platform$OS.type == "unix") { + if (file.exists(file.path(test_path, "ss3"))) { + skip_test <- FALSE + } } test_that("Do profile using the simple model", { - path <- file.path(runs_path) - get <- get_settings_profile( - parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), - low = c(0.20, 0.40, -0.5), - high = c(0.20, 1.0, 0.5), - step_size = c(0.01, 0.2, 0.25), - param_space = c('multiplier', 'real', 'relative')) - - model_settings <- get_settings( - settings = list( - base_name = "simple_small", - run = c("profile"), - profile_details = get ) - ) - - run_diagnostics(mydir = path, model_settings = model_settings) - - profile_para <- c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)") - check <- 0 - for (a in 1:length(profile_para)){ - check <- check + - as.numeric( - file.exists( - file.path(path, - paste0(model_settings$base_name, "_profile_", - profile_para[a]), - paste0(profile_para[a], "_quant_table.csv"))) - ) - } - expect_true(check == length(profile_para)) + get <- get_settings_profile( + parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), + low = c(0.20, 0.40, -0.5), + high = c(0.20, 1.0, 0.5), + step_size = c(0.01, 0.2, 0.25), + param_space = c("multiplier", "real", "relative") + ) + + model_settings <- get_settings( + settings = list( + base_name = "simple_small", + run = c("profile"), + profile_details = get + ) + ) + + run_diagnostics(mydir = path, model_settings = model_settings) + + profile_para <- c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)") + check <- 0 + for (a in 1:length(profile_para)) { + check <- check + + as.numeric( + file.exists( + file.path( + path, + paste0( + model_settings$base_name, "_profile_", + profile_para[a] + ), + paste0(profile_para[a], "_quant_table.csv") + ) + ) + ) + } + expect_true(check == length(profile_para)) }) test_that("Do jitters using the simple model", { - path <- file.path(runs_path) - model_settings <- get_settings( - settings = list(base_name = "simple_small", - run = c("jitter"), - Njitter = 3)) + model_settings <- get_settings( + settings = list( + base_name = "simple_small", + run = c("jitter"), + Njitter = 3 + ) + ) - run_diagnostics(mydir = path, model_settings = model_settings) + run_diagnostics(mydir = path, model_settings = model_settings) - jitter_folder <- file.path( - path, - paste0(model_settings$base_name, "_jitter_", model_settings$jitter_fraction)) + jitter_folder <- file.path( + path, + paste0(model_settings$base_name, "_jitter_", model_settings$jitter_fraction) + ) - expect_true(file.exists(file.path(jitter_folder, "jitter_results.csv"))) + expect_true(file.exists(file.path(jitter_folder, "jitter_results.csv"))) }) test_that("Do retrospectives using the simple model", { - path <- file.path(runs_path) - model_settings <- get_settings( - settings = list(base_name = "simple_small", - run = c("retro"), - retro_yrs = -1:-2)) + model_settings <- get_settings( + settings = list( + base_name = "simple_small", + run = c("retro"), + retro_yrs = -1:-2 + ) + ) - run_diagnostics(mydir = path, model_settings = model_settings) + run_diagnostics(mydir = path, model_settings = model_settings) n <- length(model_settings$retro_yrs) - check <- file.exists( - file.path(path, paste0(model_settings$base_name, "_retro_", n, "_yr_peel"), "mohnsrho.csv")) - expect_true(check) + check <- file.exists( + file.path(path, paste0(model_settings$base_name, "_retro_", n, "_yr_peel"), "mohnsrho.csv") + ) + expect_true(check) })