From c4f328d7dbd0d3bd3dac5400612ae56f1f53ca83 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Fri, 5 May 2023 07:53:30 -0700 Subject: [PATCH 01/28] add inputs for management target lines --- R/get_settings.R | 4 ++++ R/profile_plot.R | 22 +++++++++++++++------- R/profile_wrapper.R | 4 ++++ R/retro_wrapper.R | 4 ++++ 4 files changed, 27 insertions(+), 7 deletions(-) diff --git a/R/get_settings.R b/R/get_settings.R index 15eb553..e510666 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -42,6 +42,10 @@ get_settings <- function(settings = NULL, verbose = FALSE) { retro_yrs = -1:-5, overwrite = TRUE, show_in_console = FALSE, + + # Plot target lines + btarg = NULL, + minbthresh = NULL, # Profile Settings remove_files = TRUE, diff --git a/R/profile_plot.R b/R/profile_plot.R index f14858e..76001ed 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -18,6 +18,7 @@ #' @seealso [profile_wrapper] and [rerun_profile_vals] call `profile_plot`. profile_plot <- 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)", @@ -147,7 +148,8 @@ profile_plot <- function(mydir, rep, para, profilesummary) { # Get the relative management targets - only grab the first element since the targets should be the same btarg <- as.numeric(profilesummary$btargs[1]) - thresh <- ifelse(btarg == 0.40, 0.25, 0.125) + #thresh <- as.numeric(profilesummary$minbthreshs[1]) + thresh <- ifelse(btarg == 0.40, 0.25, ifelse(btarg == 0.25, 0.125, -1)) pngfun(wd = mydir, file = paste0("parameter_panel_", para, ".png"), h = 7, w = 7) par(mfrow = c(2, 2), mar = c(4, 4, 2, 2), oma = c(1, 1, 1, 1)) @@ -163,11 +165,14 @@ profile_plot <- function(mydir, rep, para, profilesummary) { # parameter vs. final depletion plot(x, depl, type = "l", lwd = 2, xlab = label, ylab = "Fraction of unfished", ylim = c(0, 1.2)) points(est, depl_est, pch = 21, col = "black", bg = "blue", cex = 1.5) - abline(h = c(thresh, btarg), lty = c(2, 2), col = c("red", "darkgreen")) - legend("bottomright", - legend = c("Management target", "Minimum stock size threshold"), - lty = 2, col = c("red", "darkgreen"), bty = "n" - ) + abline(h = c(btarg, thresh), lty = c(2, 2), col = c("red", "darkgreen")) + if(btarg > 0){ + legend("bottomright", + legend = c("Management target", "Minimum stock size threshold"), + lty = 2, col = c("red", "darkgreen"), bty = "n" + ) + } + # parameter vs. SB0 plot(x, sb0, type = "l", lwd = 2, xlab = label, ylab = expression(SB[0]), ylim = c(0, max(sb0))) @@ -202,7 +207,10 @@ profile_plot <- function(mydir, rep, para, profilesummary) { ifelse(est == x, " (base)", "") ), ylimAdj = 1.15, - plotdir = mydir, subplots = c(1, 3), + btarg = btarg, + minbthresh = thresh, + plotdir = mydir, + subplots = c(1, 3), pdf = FALSE, print = TRUE, plot = FALSE, filenameprefix = paste0(para, "_trajectories_") ) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 15594a8..e3bc650 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -169,6 +169,10 @@ profile_wrapper <- function(mydir, model_settings) { profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) profilesummary <- r4ss::SSsummarize(biglist = profilemodels) + if(!is.null(model_settings$btarg)){ + profilesummary$btargs <- model_settings$btarg + profilesummary$minbthreshs <- model_settings$minbthresh + } profile_output <- list() profile_output$mydir <- profile_dir diff --git a/R/retro_wrapper.R b/R/retro_wrapper.R index 9cf7268..1069804 100644 --- a/R/retro_wrapper.R +++ b/R/retro_wrapper.R @@ -145,6 +145,8 @@ retro_wrapper <- function(mydir, model_settings) { ifelse(abs(model_settings$retro_yrs) == 1, "", "s") ) ), + btarg = model_settings$btarg, + minbthresh = model_settings$minbthresh, plotdir = retro_dir, legendloc = "topright", print = TRUE, @@ -158,6 +160,8 @@ retro_wrapper <- function(mydir, model_settings) { endyrvec = endyrvec, legendloc = "topleft", plotdir = retro_dir, + btarg = model_settings$btarg, + minbthresh = model_settings$minbthresh, print = TRUE, plot = FALSE, pdf = FALSE ), subplot = c(2, 4), From 75a8d40c24ee8b4307f36458b7ad2dfe8711fa1a Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 14 Jun 2023 07:37:32 -0700 Subject: [PATCH 02/28] adjust endyr The code was defining the endyr based on the r4ss$endyr object which is equal to the final model year with data. This adjusts the endyr to be equal to the first post-data year to align with West Coast assessments --- R/get_summary.R | 8 ++++---- R/profile_plot.R | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/get_summary.R b/R/get_summary.R index 44ab6e9..9054dd9 100644 --- a/R/get_summary.R +++ b/R/get_summary.R @@ -30,8 +30,8 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) { "likelihood" = sapply(sapply(outputs, "[[", "likelihoods_used", simplify = FALSE), "[", 1, 1), "gradient" = sapply(outputs, "[[", "maximum_gradient_component"), "SB0" = sapply(quants, "[[", "SSB_Virgin", "Value"), - "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]]$endyr), "Value"), - "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr), "Value"), + "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]]$endyr + 1), "Value"), + "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr + 1), "Value"), # "Fmsy" = sapply(quants, "[[", "annF_MSY", "Value"), "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))), stringsAsFactors = FALSE @@ -57,8 +57,8 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) { parmlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_devs", 1:n]), R0 = as.numeric(x$pars[x$pars$Label == "SR_LN(R0)", 1:n]), SB0 = as.numeric(x$SpawnBio[x$SpawnBio$Label == "SSB_Virgin", 1:n]), - SBfinal = as.numeric(x$SpawnBio[x$SpawnBio$Label == paste0("SSB_", endyr), 1:n]), - deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr), 1:n]), + SBfinal = as.numeric(x$SpawnBio[x$SpawnBio$Label == paste0("SSB_", endyr + 1), 1:n]), + deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr + 1), 1:n]), yieldspr = as.numeric(x$quants[x$quants$Label == "Dead_Catch_SPR", 1:n]), steep = as.numeric(x$pars[x$pars$Label == "SR_BH_steep", 1:n]), mfem = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Fem_GP_1", 1:n]), diff --git a/R/profile_plot.R b/R/profile_plot.R index 76001ed..e9a264b 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -131,7 +131,7 @@ profile_plot <- function(mydir, rep, para, profilesummary) { dev.off() - maxyr <- min(profilesummary$endyrs) + maxyr <- min(profilesummary$endyrs + 1) minyr <- max(profilesummary$startyrs) est <- rep$parameters[rep$parameters$Label == para, "Value", 2] @@ -210,7 +210,7 @@ profile_plot <- function(mydir, rep, para, profilesummary) { btarg = btarg, minbthresh = thresh, plotdir = mydir, - subplots = c(1, 3), + subplots = c(1, 3, 9, 10, 11, 12), pdf = FALSE, print = TRUE, plot = FALSE, filenameprefix = paste0(para, "_trajectories_") ) From e42f14b80085e9c67875df9bdaa05ec97c42ecf7 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 14 Jun 2023 07:51:35 -0700 Subject: [PATCH 03/28] deprecate the prior like input Deprecate the use_prior_like function input for parameter profiles (issue #6) since r4ss is now dynamically showing this in profiles if a prior likelihood is used r4ss/r4ss#822 --- R/get_settings_profile.R | 29 ++++++++++++++++++----------- R/profile_plot.R | 1 - R/profile_wrapper.R | 5 +---- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/R/get_settings_profile.R b/R/get_settings_profile.R index 5737c34..ba131d8 100644 --- a/R/get_settings_profile.R +++ b/R/get_settings_profile.R @@ -21,7 +21,7 @@ #' @param param_space options: real, mulitplier, relative indicates how to interpret the low and high bound values. #' real indicates bounds in the parameter space, relative indicates how far to go from the base parameter, and #' multiplier indicates that low and high bounds are set at x\% above and below the base parameter. -#' @param use_prior_like options: Option to include or exclude the prior likelihoods in the likelihood +#' @param use_prior_like Deprecated Option to include or exclude the prior likelihoods in the likelihood #' profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude #' and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation #' (e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood @@ -41,11 +41,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'), -#' use_prior_like = c(1, 1, 0) +#' 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 @@ -54,8 +53,7 @@ #' low = c(0.40, 0.40), #' high = c(0.40, 0.40), #' step_size = c(0.005, 0.005), -#' param_space = c("multiplier", "multiplier"), -#' use_prior_like = c(0, 1) +#' param_space = c("multiplier", "multiplier") #' ) #' } #' @@ -64,7 +62,9 @@ get_settings_profile <- function(parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_ high = c(0.40, 1.0, 2), step_size = c(0.01, 0.05, 0.25), param_space = c("multiplier", "real", "relative"), - use_prior_like = c(0, 0, 0)) { + use_prior_like = lifecycle::deprecated()) +{ + if (length(parameters) != length(low) | length(parameters) != length(high) | length(parameters) != length(step_size) | @@ -73,13 +73,20 @@ get_settings_profile <- function(parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_ stop("Error: input vectors do match in length.") } + if (lifecycle::is_present(use_prior_like)) { + lifecycle::deprecate_warn( + #when = "1.46.0", + what = "get_settings_profile(use_prior_like)", + details = "Input 'use_prior_like' no longer needed." + ) + } + out <- data.frame( parameters = parameters, low = low, high = high, step_size = step_size, - param_space = param_space, - use_prior_like = use_prior_like + param_space = param_space ) return(out) diff --git a/R/profile_plot.R b/R/profile_plot.R index e9a264b..52613a7 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -173,7 +173,6 @@ profile_plot <- function(mydir, rep, para, profilesummary) { ) } - # parameter vs. SB0 plot(x, sb0, type = "l", lwd = 2, xlab = label, ylab = expression(SB[0]), ylim = c(0, max(sb0))) points(est, sb0_est, pch = 21, col = "black", bg = "blue", cex = 1.5) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index e3bc650..7042367 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -37,10 +37,9 @@ profile_wrapper <- function(mydir, model_settings) { for (aa in 1:N) { para <- model_settings$profile_details$parameters[aa] - prior_used <- model_settings$profile_details$use_prior_like[aa] # 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, "_prior_like_", prior_used)) + profile_dir <- file.path(mydir, paste0(model_settings$base_name, "_profile_", para)) dir.create(profile_dir, showWarnings = FALSE) # Check for existing files and delete @@ -87,8 +86,6 @@ profile_wrapper <- function(mydir, model_settings) { starter <- r4ss::SS_readstarter(file = file.path(profile_dir, "starter.ss")) starter$ctlfile <- "control_modified.ss" starter$init_values_src <- model_settings$init_values_src - # make sure the prior likelihood is calculated for non-estimated quantities - starter$prior_like <- prior_used r4ss::SS_writestarter(mylist = starter, dir = profile_dir, overwrite = TRUE) # Read in the base model From e83968e372f5de172b201665557b4dacb4f23e98 Mon Sep 17 00:00:00 2001 From: Ian Taylor Date: Fri, 16 Jun 2023 10:40:38 -0700 Subject: [PATCH 04/28] add user control of oldctlfile for profiles --- R/get_settings.R | 1 + R/profile_wrapper.R | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/R/get_settings.R b/R/get_settings.R index 15eb553..00dec79 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -45,6 +45,7 @@ get_settings <- function(settings = NULL, verbose = FALSE) { # Profile Settings remove_files = TRUE, + oldctlfile = "control.ss_new", newctlfile = "control_modified.ss", linenum = NULL, string = NULL, diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 15594a8..1dfd3bc 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -143,7 +143,7 @@ profile_wrapper <- function(mydir, model_settings) { profile <- r4ss::profile( dir = profile_dir, - oldctlfile = "control.ss_new", + oldctlfile = model_settings$oldctlfile, newctlfile = model_settings$newctlfile, linenum = model_settings$linenum, string = para, From c44d1166725a7c8880ab840b68319b37b175ad21 Mon Sep 17 00:00:00 2001 From: Ian Taylor Date: Mon, 19 Jun 2023 16:57:42 -0700 Subject: [PATCH 05/28] feature: add up-then-down approach to profile #23 - also better implement control of oldctlfile and newctlfile --- R/profile_wrapper.R | 109 ++++++++++++++++++++++++++++------------- R/rerun_profile_vals.R | 14 +++--- 2 files changed, 82 insertions(+), 41 deletions(-) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 1dfd3bc..f8dc6f4 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -56,21 +56,31 @@ profile_wrapper <- function(mydir, model_settings) { ), file = "run_diag_warning.txt") message(paste0("Running profile for ", para, ".")) - # Copy the control file to run from the copy - if (!file.exists(file.path(profile_dir, "control.ss_new"))) { - orig_dir <- getwd() - setwd(profile_dir) - r4ss::run( - dir = profile_dir, - exe = model_settings$exe, - extras = model_settings$extras - ) - setwd(orig_dir) + # check for whether oldctlfile exists + if (!file.exists(file.path(profile_dir, model_settings$oldctlfile))) { + # if the oldctlfile is control.ss_new, and doesn't exist, + # run the model to create it + if (model_settings$oldctlfile == "control.ss_new") { + if (model_settings$verbose) { + "running model to get control.ss_new file" + } + orig_dir <- getwd() + setwd(profile_dir) + r4ss::run( + dir = profile_dir, + exe = model_settings$exe, + extras = model_settings$extras + ) + setwd(orig_dir) + } else { + # stop if the file isn't control.ss_new and doesn't exist + stop("File not found ", model_settings$oldctlfile) + } } # Use the SS_parlines funtion to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( - ctlfile = "control.ss_new", + ctlfile = model_settings$oldctlfile, dir = profile_dir, verbose = FALSE, version = model_settings$version, @@ -79,13 +89,17 @@ profile_wrapper <- function(mydir, model_settings) { if (sum(check_para) == 0) { print(para) - stop(paste0("The input profile_custom does not match a parameter in the control.ss_new file.")) + stop("The input profile_custom does not match a parameter in the file ", + model_settings$oldctlfile) } - file.copy(file.path(profile_dir, "control.ss_new"), file.path(profile_dir, model_settings$newctlfile)) + # Copy oldctlfile to newctlfile before modifying it + 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")) - starter$ctlfile <- "control_modified.ss" + starter$ctlfile <- model_settings$newctlfile starter$init_values_src <- model_settings$init_values_src # make sure the prior likelihood is calculated for non-estimated quantities starter$prior_like <- prior_used @@ -137,30 +151,55 @@ profile_wrapper <- function(mydir, model_settings) { } else { high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size)) } - +browser() vec <- c(low, high) num <- sort(vec, index.return = TRUE)$ix - profile <- r4ss::profile( - dir = profile_dir, - oldctlfile = model_settings$oldctlfile, - newctlfile = model_settings$newctlfile, - linenum = model_settings$linenum, - string = para, - profilevec = vec, - usepar = model_settings$usepar, - globalpar = model_settings$globalpar, - parlinenum = model_settings$parlinenum, - parstring = model_settings$parstring, - saveoutput = model_settings$saveoutput, - overwrite = model_settings$overwrite, - whichruns = model_settings$whichruns, - prior_check = model_settings$prior_check, - read_like = model_settings$read_like, - exe = model_settings$exe, - verbose = model_settings$verbose, - extras = model_settings$extras - ) + # loop over down, then up + for (iprofile in 1:2) { + if (iprofile == 1) { + # subset of requested runs which are in the "low" vector + whichruns <- which(vec %in% low) + if (!is.null(model_settings$whichruns)) { + whichruns <- intersect(model_settings$whichruns, whichruns) + } + if (model_settings$verbose) { + message("Running profiles down from estimated value:", + paste(whichruns, sep = " ")) + } + } else { + # subset of requested runs which are in the "low" vector + whichruns <- which(vec %in% high) + if (!is.null(model_settings$whichruns)) { + whichruns <- intersect(model_settings$whichruns, whichruns) + } + if (model_settings$verbose) { + message("Running profiles down from estimated value:", + paste(whichruns, sep = " ")) + } + } + + profile <- r4ss::profile( + dir = profile_dir, + oldctlfile = model_settings$oldctlfile, + newctlfile = model_settings$newctlfile, + linenum = model_settings$linenum, + string = para, + profilevec = vec, + usepar = model_settings$usepar, + globalpar = model_settings$globalpar, + parlinenum = model_settings$parlinenum, + parstring = model_settings$parstring, + saveoutput = model_settings$saveoutput, + overwrite = model_settings$overwrite, + whichruns = whichruns, # values set above + prior_check = model_settings$prior_check, + read_like = model_settings$read_like, + exe = model_settings$exe, + verbose = model_settings$verbose, + extras = model_settings$extras + ) + } # Save the output and the summary name <- paste0("profile_", para) diff --git a/R/rerun_profile_vals.R b/R/rerun_profile_vals.R index a9efc94..e19a64d 100644 --- a/R/rerun_profile_vals.R +++ b/R/rerun_profile_vals.R @@ -15,7 +15,7 @@ #' @author Chantel Wetzel. #' @export #' -#' @example +#' @examples #' rerun_profile_vals(mydir = file.path(directory, "base_model"), #' model_settings = model_settings, #' para_name = "NatM_uniform_Fem_GP_1", @@ -46,8 +46,8 @@ rerun_profile_vals <- function(mydir, dir.create(temp_dir, showWarnings = FALSE) file.copy(file.path(profile_dir, "ss.exe"), temp_dir, overwrite = TRUE) - file.copy(file.path(profile_dir, "control_modified.ss"), temp_dir) - file.copy(file.path(profile_dir, "control.ss_new"), temp_dir) + file.copy(file.path(profile_dir, model_settings$newctlfile), temp_dir) + file.copy(file.path(profile_dir, model_settings$oldctlfile), temp_dir) file.copy(file.path(profile_dir, data_file_nm), temp_dir) file.copy(file.path(profile_dir, "starter.ss_new"), temp_dir) file.copy(file.path(profile_dir, "forecast.ss_new"), temp_dir) @@ -56,13 +56,15 @@ rerun_profile_vals <- function(mydir, # Use the SS_parlines funtion to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( + ctlfile = model_settings$oldctlfile, dir = temp_dir, verbose = FALSE, active = FALSE )$Label == para if (sum(check_para) == 0) { - stop(paste0("The input profile_custom does not match a parameter in the control.ss_new file.")) + stop("The input profile_custom does not match a parameter in the file", + model_settings$oldctlfile) } load(file.path(profile_dir, paste0(para_name, "_profile_output.Rdata"))) @@ -80,8 +82,8 @@ rerun_profile_vals <- function(mydir, for (i in run_num) { setwd(temp_dir) r4ss::SS_changepars( - ctlfile = "control_modified.ss", - newctlfile = "control_modified.ss", + ctlfile = model_settings$newctlfile, + newctlfile = model_settings$newctlfile, strings = para, newvals = vec[i], estimate = FALSE From 482a1820bffb644bed3eb0420b67525d1e16f3a7 Mon Sep 17 00:00:00 2001 From: Ian Taylor Date: Mon, 19 Jun 2023 16:57:53 -0700 Subject: [PATCH 06/28] devtools::document() --- DESCRIPTION | 2 +- man/get_settings_profile.Rd | 6 +++--- man/rerun_profile_vals.Rd | 8 ++++++++ 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 932accf..e9c8645 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,6 +24,6 @@ Suggests: devtools Remotes: github::r4ss/r4ss -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 Encoding: UTF-8 Roxygen: list(markdown = TRUE) diff --git a/man/get_settings_profile.Rd b/man/get_settings_profile.Rd index 98c17af..104036f 100644 --- a/man/get_settings_profile.Rd +++ b/man/get_settings_profile.Rd @@ -5,7 +5,7 @@ \title{Get Default Settings For Profiles} \usage{ get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), + parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), low = c(0.4, 0.25, -2), high = c(0.4, 1, 2), step_size = c(0.01, 0.05, 0.25), @@ -59,7 +59,7 @@ parameter. A user can select any of the options for specifying a parameter range # Define each parameter in real space get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), + 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), @@ -69,7 +69,7 @@ get_settings_profile( # Example 2: Run a profile for natural mortality one with the prior likelihood and one without get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "NatM_p_1_Fem_GP_1"), + parameters = c("NatM_uniform_Fem_GP_1", "NatM_uniform_Fem_GP_1"), low = c(0.40, 0.40), high = c(0.40, 0.40), step_size = c(0.005, 0.005), diff --git a/man/rerun_profile_vals.Rd b/man/rerun_profile_vals.Rd index 3d2c944..d80e4b3 100644 --- a/man/rerun_profile_vals.Rd +++ b/man/rerun_profile_vals.Rd @@ -31,6 +31,14 @@ you would like to rerun.} Generate likelihood profiles To be called from the run_diagnostics function after creating the model settings using the get_settings function. +} +\examples{ +rerun_profile_vals(mydir = file.path(directory, "base_model"), + model_settings = model_settings, + para_name = "NatM_uniform_Fem_GP_1", + run_num = c(3, 4), + data_file_nm = "data.ss") + } \author{ Chantel Wetzel. From b27cf4767189d812939171c15fe5930670fdcf36 Mon Sep 17 00:00:00 2001 From: Ian Taylor Date: Mon, 19 Jun 2023 17:16:38 -0700 Subject: [PATCH 07/28] remove browser(), oops --- R/profile_wrapper.R | 4 ++-- R/rerun_profile_vals.R | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index f8dc6f4..e05a40e 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -78,7 +78,7 @@ profile_wrapper <- function(mydir, model_settings) { } } - # Use the SS_parlines funtion to ensure that the input parameter can be found + # Use the SS_parlines function to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( ctlfile = model_settings$oldctlfile, dir = profile_dir, @@ -151,7 +151,7 @@ profile_wrapper <- function(mydir, model_settings) { } else { high <- c(seq(round_any(est, step_size, f = ceiling), range[2], step_size)) } -browser() + vec <- c(low, high) num <- sort(vec, index.return = TRUE)$ix diff --git a/R/rerun_profile_vals.R b/R/rerun_profile_vals.R index e19a64d..1ab81aa 100644 --- a/R/rerun_profile_vals.R +++ b/R/rerun_profile_vals.R @@ -54,7 +54,7 @@ rerun_profile_vals <- function(mydir, file.rename(file.path(temp_dir, "starter.ss_new"), file.path(temp_dir, "starter.ss")) file.rename(file.path(temp_dir, "forecast.ss_new"), file.path(temp_dir, "forecast.ss")) - # Use the SS_parlines funtion to ensure that the input parameter can be found + # Use the SS_parlines function to ensure that the input parameter can be found check_para <- r4ss::SS_parlines( ctlfile = model_settings$oldctlfile, dir = temp_dir, From 4f8cd88d031bb65552edf92718579092a1ae2a2b Mon Sep 17 00:00:00 2001 From: Ian Taylor Date: Mon, 19 Jun 2023 19:32:42 -0700 Subject: [PATCH 08/28] fix: backup control.ss_new for 2nd half of profile --- R/profile_wrapper.R | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index e05a40e..a6ed02b 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -155,6 +155,10 @@ profile_wrapper <- function(mydir, model_settings) { vec <- c(low, high) num <- sort(vec, index.return = TRUE)$ix + # 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")) + # loop over down, then up for (iprofile in 1:2) { if (iprofile == 1) { @@ -177,6 +181,9 @@ profile_wrapper <- function(mydir, model_settings) { message("Running profiles down from estimated value:", paste(whichruns, sep = " ")) } + # 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)) } profile <- r4ss::profile( From 2b29719775d9f0724e88e50763a857f99034f02f Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 21 Jun 2023 16:23:57 -0700 Subject: [PATCH 09/28] rm prior like check --- R/get_settings_profile.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/R/get_settings_profile.R b/R/get_settings_profile.R index ba131d8..f9b7a1a 100644 --- a/R/get_settings_profile.R +++ b/R/get_settings_profile.R @@ -68,8 +68,7 @@ get_settings_profile <- function(parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_ if (length(parameters) != length(low) | length(parameters) != length(high) | length(parameters) != length(step_size) | - length(parameters) != length(param_space) | - length(parameters) != length(use_prior_like)) { + length(parameters) != length(param_space)) { stop("Error: input vectors do match in length.") } From 9ca8d97d5ad493ff33bd69fb954938610732cd77 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Thu, 22 Jun 2023 07:59:44 -0700 Subject: [PATCH 10/28] add ability to request specific plots add adjustment depending upon prior inclusion correction to plot the base model point on the 4-panel parameter plot to be on the profile line regardless of the setting in the starter file to include the prior likelihood component. --- R/profile_plot.R | 16 ++++++++++++---- R/profile_wrapper.R | 3 ++- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/R/profile_plot.R b/R/profile_plot.R index 52613a7..9dad09b 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -133,14 +133,21 @@ profile_plot <- function(mydir, rep, para, profilesummary) { maxyr <- min(profilesummary$endyrs + 1) minyr <- max(profilesummary$startyrs) - + est <- rep$parameters[rep$parameters$Label == para, "Value", 2] sb0_est <- rep$derived_quants[rep$derived_quants$Label == "SSB_Virgin", "Value"] sbf_est <- rep$derived_quants[rep$derived_quants$Label == paste0("SSB_", maxyr), "Value"] depl_est <- rep$derived_quants[rep$derived_quants$Label == paste0("Bratio_", maxyr), "Value"] x <- as.numeric(profilesummary$pars[profilesummary$pars$Label == para, n]) - like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - rep$likelihoods_used[1, 1]) + # determine whether to include the prior likelihood component in the likelihood profile + starter <- r4ss::SS_readstarter(file = file.path(mydir, "starter.ss")) + if (starter$prior_like == 0) { + like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - + profilesummary$likelihoods[profilesummary$likelihoods$Label == "Parm_priors", n] - rep$likelihoods_used[1, 1]) + } else { + like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - rep$likelihoods_used[1, 1]) + } ylike <- c(min(like) + ifelse(min(like) != 0, -0.5, 0), max(like)) sb0 <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Label) == "SSB_Virgin", n]) sbf <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Yr) == maxyr, n]) @@ -148,7 +155,6 @@ profile_plot <- function(mydir, rep, para, profilesummary) { # Get the relative management targets - only grab the first element since the targets should be the same btarg <- as.numeric(profilesummary$btargs[1]) - #thresh <- as.numeric(profilesummary$minbthreshs[1]) thresh <- ifelse(btarg == 0.40, 0.25, ifelse(btarg == 0.25, 0.125, -1)) pngfun(wd = mydir, file = paste0("parameter_panel_", para, ".png"), h = 7, w = 7) @@ -194,6 +200,8 @@ profile_plot <- function(mydir, rep, para, profilesummary) { ) ) ) + + subplots <- profilesummary$subplots r4ss::SSplotComparisons( summaryoutput = profilesummary, @@ -209,7 +217,7 @@ profile_plot <- function(mydir, rep, para, profilesummary) { btarg = btarg, minbthresh = thresh, plotdir = mydir, - subplots = c(1, 3, 9, 10, 11, 12), + subplots = subplots, pdf = FALSE, print = TRUE, plot = FALSE, filenameprefix = paste0(para, "_trajectories_") ) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 7042367..95b6033 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -170,7 +170,8 @@ profile_wrapper <- function(mydir, model_settings) { profilesummary$btargs <- model_settings$btarg profilesummary$minbthreshs <- model_settings$minbthresh } - + profilesummary$subplots <- model_settings$subplots + profile_output <- list() profile_output$mydir <- profile_dir profile_output$para <- para From a6c66060eb5c726fc4d0118fdecadd8e5f753913 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Thu, 22 Jun 2023 08:00:27 -0700 Subject: [PATCH 11/28] typo --- R/get_summary.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/get_summary.R b/R/get_summary.R index 9054dd9..fbb6427 100644 --- a/R/get_summary.R +++ b/R/get_summary.R @@ -15,7 +15,7 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) { - # Need to identify a way to determine if a model estmates male growth parameters as offsets from females + # Need to identify a way to determine if a model estimates male growth parameters as offsets from females # get output outputs <- profilemodels From a7d15d6678c5b964b946f9352370eb6a1f0b640c Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Thu, 22 Jun 2023 08:00:42 -0700 Subject: [PATCH 12/28] add default subplot seelction --- R/get_settings.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/R/get_settings.R b/R/get_settings.R index e510666..684706e 100644 --- a/R/get_settings.R +++ b/R/get_settings.R @@ -17,6 +17,7 @@ #' settings <- list() #' get_settings <- function(settings = NULL, verbose = FALSE) { + if (is.vector(settings)) settings <- as.list(settings) Settings_all <- list( @@ -62,7 +63,8 @@ get_settings <- function(settings = NULL, verbose = FALSE) { whichruns = NULL, prior_check = FALSE, read_like = TRUE, - init_values_src = 0 + init_values_src = 0, + subplots = c(1, 3) ) Settings_all$profile_details <- get_settings_profile() @@ -84,7 +86,7 @@ get_settings <- function(settings = NULL, verbose = FALSE) { if (!is.numeric(Settings_all$profile_details$low) & !is.numeric(Settings_all$profile_details$high) & !is.numeric(Settings_all$profile_details$step_size)) { - stop("There is a non-numeric value in the low, high, or step size fiedl of the get_settings_profile data frame.") + stop("There is a non-numeric value in the low, high, or step size field of the get_settings_profile data frame.") } if (sum(!Settings_all$profile_details$param_space %in% c("real", "relative", "multiplier")) > 0) { stop("The param_space column should be either real or relative in the get_settings_profile data frame.") From 4d9a23e6f7bf37cfc3252b438a976c12ea9d213b Mon Sep 17 00:00:00 2001 From: Ian Taylor Date: Thu, 22 Jun 2023 09:49:53 -0700 Subject: [PATCH 13/28] more attempts to work on up/down profile --- R/profile_wrapper.R | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index a6ed02b..b7debd1 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -157,7 +157,13 @@ profile_wrapper <- function(mydir, model_settings) { # 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")) + 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) # loop over down, then up for (iprofile in 1:2) { @@ -169,23 +175,31 @@ profile_wrapper <- function(mydir, model_settings) { } if (model_settings$verbose) { message("Running profiles down from estimated value:", - paste(whichruns, sep = " ")) + paste(whichruns, collapse = " "), + "out of ", + length(vec) } - } else { + } + if (iprofile == 2) { # subset of requested runs which are in the "low" vector whichruns <- which(vec %in% high) if (!is.null(model_settings$whichruns)) { whichruns <- intersect(model_settings$whichruns, whichruns) } if (model_settings$verbose) { - message("Running profiles down from estimated value:", - paste(whichruns, sep = " ")) + message("Running profiles up from estimated value:", + paste(whichruns, collapse = " "), + "out of ", + length(vec)) } # 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)) + # 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) } - profile <- r4ss::profile( dir = profile_dir, oldctlfile = model_settings$oldctlfile, From 704db5b7753a5325e919bf5f257d584fbdcdcc70 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Thu, 22 Jun 2023 11:23:30 -0700 Subject: [PATCH 14/28] fix add missing ) --- R/profile_wrapper.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index b7debd1..66caf24 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -177,7 +177,7 @@ profile_wrapper <- function(mydir, model_settings) { message("Running profiles down from estimated value:", paste(whichruns, collapse = " "), "out of ", - length(vec) + length(vec)) } } if (iprofile == 2) { From e052a48c68ea633f35393458b15b7c6e5935c6e9 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Thu, 22 Jun 2023 11:26:52 -0700 Subject: [PATCH 15/28] doc --- DESCRIPTION | 2 +- man/get_settings_profile.Rd | 22 ++++++++++------------ 2 files changed, 11 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 932accf..e9c8645 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -24,6 +24,6 @@ Suggests: devtools Remotes: github::r4ss/r4ss -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 Encoding: UTF-8 Roxygen: list(markdown = TRUE) diff --git a/man/get_settings_profile.Rd b/man/get_settings_profile.Rd index 98c17af..0f95438 100644 --- a/man/get_settings_profile.Rd +++ b/man/get_settings_profile.Rd @@ -5,12 +5,12 @@ \title{Get Default Settings For Profiles} \usage{ get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), + parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_steep", "SR_LN(R0)"), low = c(0.4, 0.25, -2), high = c(0.4, 1, 2), step_size = c(0.01, 0.05, 0.25), param_space = c("multiplier", "real", "relative"), - use_prior_like = c(0, 0, 0) + use_prior_like = lifecycle::deprecated() ) } \arguments{ @@ -26,7 +26,7 @@ get_settings_profile( real indicates bounds in the parameter space, relative indicates how far to go from the base parameter, and multiplier indicates that low and high bounds are set at x\\% above and below the base parameter.} -\item{use_prior_like}{options: Option to include or exclude the prior likelihoods in the likelihood +\item{use_prior_like}{Deprecated Option to include or exclude the prior likelihoods in the likelihood profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation (e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood @@ -59,22 +59,20 @@ parameter. A user can select any of the options for specifying a parameter range # Define each parameter in real space get_settings_profile( - parameters = c("NatM_p_1_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'), - use_prior_like = c(1, 1, 0) + 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') ) # Example 2: Run a profile for natural mortality one with the prior likelihood and one without get_settings_profile( - parameters = c("NatM_p_1_Fem_GP_1", "NatM_p_1_Fem_GP_1"), + parameters = c("NatM_uniform_Fem_GP_1", "NatM_uniform_Fem_GP_1"), low = c(0.40, 0.40), high = c(0.40, 0.40), step_size = c(0.005, 0.005), - param_space = c("multiplier", "multiplier"), - use_prior_like = c(0, 1) + param_space = c("multiplier", "multiplier") ) } From 03ad041bde075feddde6234b0234b86b6d05a4f4 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Fri, 23 Jun 2023 15:45:50 -0700 Subject: [PATCH 16/28] fix M name in summary table --- R/get_summary.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/get_summary.R b/R/get_summary.R index fbb6427..510cdc6 100644 --- a/R/get_summary.R +++ b/R/get_summary.R @@ -61,13 +61,13 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) { deplfinal = as.numeric(x$Bratio[x$Bratio$Label == paste0("Bratio_", endyr + 1), 1:n]), yieldspr = as.numeric(x$quants[x$quants$Label == "Dead_Catch_SPR", 1:n]), steep = as.numeric(x$pars[x$pars$Label == "SR_BH_steep", 1:n]), - mfem = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Fem_GP_1", 1:n]), + mfem = as.numeric(x$pars[x$pars$Label == "NatM_uniform_Fem_GP_1", 1:n]), lminfem = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Fem_GP_1", 1:n]), lmaxfem = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Fem_GP_1", 1:n]), kfem = as.numeric(x$pars[x$pars$Label == "VonBert_K_Fem_GP_1", 1:n]), cv1fem = as.numeric(x$pars[grep("young_Fem_GP_1", x$pars$Label), 1:n]), cv2fem = as.numeric(x$pars[grep("old_Fem_GP_1", x$pars$Label), 1:n]), - mmale = as.numeric(x$pars[x$pars$Label == "NatM_p_1_Mal_GP_1", 1:n]), + mmale = as.numeric(x$pars[x$pars$Label == "NatM_uniform_Mal_GP_1", 1:n]), lminmale = as.numeric(x$pars[x$pars$Label == "L_at_Amin_Mal_GP_1", 1:n]), lmaxmale = as.numeric(x$pars[x$pars$Label == "L_at_Amax_Mal_GP_1", 1:n]), kmale = as.numeric(x$pars[x$pars$Label == "VonBert_K_Mal_GP_1", 1:n]), From 5c2159bb288af15985cf39dda3f7acfd8fb943a3 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Tue, 27 Jun 2023 07:38:58 -0700 Subject: [PATCH 17/28] Revise the output warning if no control file found --- R/profile_wrapper.R | 3 +-- R/rerun_profile_vals.R | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 3326057..40aabb4 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -72,8 +72,7 @@ profile_wrapper <- function(mydir, model_settings) { ) setwd(orig_dir) } else { - # stop if the file isn't control.ss_new and doesn't exist - stop("File not found ", model_settings$oldctlfile) + stop("Can not find ", model_settings$oldctlfile) } } diff --git a/R/rerun_profile_vals.R b/R/rerun_profile_vals.R index 1ab81aa..c5e9a06 100644 --- a/R/rerun_profile_vals.R +++ b/R/rerun_profile_vals.R @@ -63,7 +63,7 @@ rerun_profile_vals <- function(mydir, )$Label == para if (sum(check_para) == 0) { - stop("The input profile_custom does not match a parameter in the file", + stop("The input profile_custom does not match a parameter in the ", model_settings$oldctlfile) } From 711f4d8610349f7e286143b352bf5352ec51ba2a Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Tue, 27 Jun 2023 07:39:33 -0700 Subject: [PATCH 18/28] update the function input for use_prior_like Now describes why this function input was removed --- R/get_settings_profile.R | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/R/get_settings_profile.R b/R/get_settings_profile.R index f9b7a1a..d0a670f 100644 --- a/R/get_settings_profile.R +++ b/R/get_settings_profile.R @@ -21,11 +21,9 @@ #' @param param_space options: real, mulitplier, relative indicates how to interpret the low and high bound values. #' real indicates bounds in the parameter space, relative indicates how far to go from the base parameter, and #' multiplier indicates that low and high bounds are set at x\% above and below the base parameter. -#' @param use_prior_like Deprecated Option to include or exclude the prior likelihoods in the likelihood -#' profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude -#' and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation -#' (e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood -#' contribution may be wanted in specific instances. The default setting excludes the prior likelihood contributions. +#' @param use_prior_like Deprecated: The use_prior_like input is no longer needed since r4ss now +#' automatically plots the likelihood profile with and without any parameter prior likelihood contributions +#' regardless of the setting in the user starter file. #' #' @return A matrix of low, high, and step size values for the default parameters #' that should be profiled. The goal is to provide users with a template @@ -74,9 +72,8 @@ get_settings_profile <- function(parameters = c("NatM_uniform_Fem_GP_1", "SR_BH_ if (lifecycle::is_present(use_prior_like)) { lifecycle::deprecate_warn( - #when = "1.46.0", - what = "get_settings_profile(use_prior_like)", - details = "Input 'use_prior_like' no longer needed." + when = "1.1.2", + what = "get_settings_profile(use_prior_like)" ) } From 5ac35327b0b402864eeaf4c4bd1f4f495b6c7c10 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Tue, 27 Jun 2023 07:40:20 -0700 Subject: [PATCH 19/28] remove character return --- R/profile_plot.R | 1 - 1 file changed, 1 deletion(-) diff --git a/R/profile_plot.R b/R/profile_plot.R index 9dad09b..d38f3f3 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -133,7 +133,6 @@ profile_plot <- function(mydir, rep, para, profilesummary) { maxyr <- min(profilesummary$endyrs + 1) minyr <- max(profilesummary$startyrs) - est <- rep$parameters[rep$parameters$Label == para, "Value", 2] sb0_est <- rep$derived_quants[rep$derived_quants$Label == "SSB_Virgin", "Value"] sbf_est <- rep$derived_quants[rep$derived_quants$Label == paste0("SSB_", maxyr), "Value"] From ae6fdcdbe4ff3798d2e362a9ea0173b6453c95b2 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Tue, 27 Jun 2023 07:40:39 -0700 Subject: [PATCH 20/28] change version number --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index e9c8645..1f8bc95 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nwfscDiag Type: Package Title: Generate standard NWFSC assessment diagnostics -Version: 1.1.1 +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. From 8f107f01b24a7df7fd0afffcb968c25293f6b750 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 08:38:49 -0700 Subject: [PATCH 21/28] add space --- .github/workflows/r-cmd-check.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/r-cmd-check.yml b/.github/workflows/r-cmd-check.yml index 678ea96..944cd77 100644 --- a/.github/workflows/r-cmd-check.yml +++ b/.github/workflows/r-cmd-check.yml @@ -20,7 +20,7 @@ on: pull_request: # specifying branches means the workflow will only run when the pull request is to the merge into the branch listed, in this case, main. branches: - -main + - main # jobs specifies the jobs to run. you can have multiple jobs per github actions workflow, but in this case we have only one, which is # named run-R-CMD-check. jobs: From 6ea06cb6517f2c1c7817835bd484c9150600fe2c Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 08:39:46 -0700 Subject: [PATCH 22/28] update M parameter name --- R/get_summary.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/get_summary.R b/R/get_summary.R index 510cdc6..aa1ca7d 100644 --- a/R/get_summary.R +++ b/R/get_summary.R @@ -81,10 +81,10 @@ get_summary <- function(mydir, para, vec, name, profilemodels, profilesummary) { if (para == "SR_LN(R0)") { colnames(new_out) <- paste0("R0 ", vec) } - if (para == "NatM_p_1_Fem_GP_1") { + if (para == "NatM_uniform_Fem_GP_1") { colnames(new_out) <- paste0("M_f ", vec) } - if (para == "NatM_p_1_Mal_GP_1") { + if (para == "NatM_uniform_Mal_GP_1") { colnames(new_out) <- paste0("M_m ", vec) } if (para == "SR_BH_steep") { From acf380fa8525cc260046594bce9e4f2fcca571ff Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 08:59:05 -0700 Subject: [PATCH 23/28] corrections/revisions 1. Add text output into message() 2. Remove directory checks that were redundant with r4ss::run() 3. correct btarg/minbthresh checks 4. revise like calculation based on user settings for priors --- R/profile_plot.R | 16 ++++++++-------- R/profile_wrapper.R | 9 +++------ 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/R/profile_plot.R b/R/profile_plot.R index d38f3f3..7926072 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -141,20 +141,20 @@ profile_plot <- function(mydir, rep, para, profilesummary) { x <- as.numeric(profilesummary$pars[profilesummary$pars$Label == para, n]) # determine whether to include the prior likelihood component in the likelihood profile starter <- r4ss::SS_readstarter(file = file.path(mydir, "starter.ss")) - if (starter$prior_like == 0) { - like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - - profilesummary$likelihoods[profilesummary$likelihoods$Label == "Parm_priors", n] - rep$likelihoods_used[1, 1]) - } else { - like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - rep$likelihoods_used[1, 1]) - } + like <- as.numeric(profilesummary$likelihoods[profilesummary$likelihoods$Label == "TOTAL", n] - + ifelse(starter$prior_like == 0, + profilesummary$likelihoods[profilesummary$likelihoods$Label == "Parm_priors", n], + 0) - + rep$likelihoods_used[1, 1]) + ylike <- c(min(like) + ifelse(min(like) != 0, -0.5, 0), max(like)) sb0 <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Label) == "SSB_Virgin", n]) sbf <- as.numeric(profilesummary$SpawnBio[na.omit(profilesummary$SpawnBio$Yr) == maxyr, n]) depl <- as.numeric(profilesummary$Bratio[na.omit(profilesummary$Bratio$Yr) == maxyr, n]) # Get the relative management targets - only grab the first element since the targets should be the same - btarg <- as.numeric(profilesummary$btargs[1]) - thresh <- ifelse(btarg == 0.40, 0.25, ifelse(btarg == 0.25, 0.125, -1)) + btarg <- as.numeric(profilesummary$btarg[1]) + thresh <- as.numeric(profilesummary$minbthresh[1]) # ifelse(btarg == 0.40, 0.25, ifelse(btarg == 0.25, 0.125, -1)) pngfun(wd = mydir, file = paste0("parameter_panel_", para, ".png"), h = 7, w = 7) par(mfrow = c(2, 2), mar = c(4, 4, 2, 2), oma = c(1, 1, 1, 1)) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index 40aabb4..d263ac7 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -61,16 +61,13 @@ profile_wrapper <- function(mydir, model_settings) { # run the model to create it if (model_settings$oldctlfile == "control.ss_new") { if (model_settings$verbose) { - "running model to get control.ss_new file" + message("running model to get control.ss_new file") } - orig_dir <- getwd() - setwd(profile_dir) r4ss::run( dir = profile_dir, exe = model_settings$exe, extras = model_settings$extras ) - setwd(orig_dir) } else { stop("Can not find ", model_settings$oldctlfile) } @@ -226,8 +223,8 @@ profile_wrapper <- function(mydir, model_settings) { profilemodels <- r4ss::SSgetoutput(dirvec = profile_dir, keyvec = num) profilesummary <- r4ss::SSsummarize(biglist = profilemodels) if(!is.null(model_settings$btarg)){ - profilesummary$btargs <- model_settings$btarg - profilesummary$minbthreshs <- model_settings$minbthresh + profilesummary$btarg <- model_settings$btarg + profilesummary$minbthresh <- model_settings$minbthresh } profilesummary$subplots <- model_settings$subplots From b018ce4352a1a2fae286b841cd5998283bfeb6df Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 09:07:56 -0700 Subject: [PATCH 24/28] correct line colors --- R/profile_plot.R | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/R/profile_plot.R b/R/profile_plot.R index 7926072..53df251 100644 --- a/R/profile_plot.R +++ b/R/profile_plot.R @@ -170,11 +170,11 @@ profile_plot <- function(mydir, rep, para, profilesummary) { # parameter vs. final depletion plot(x, depl, type = "l", lwd = 2, xlab = label, ylab = "Fraction of unfished", 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("red", "darkgreen")) + abline(h = c(btarg, thresh), lty = c(2, 2), col = c("darkgreen", "red")) if(btarg > 0){ legend("bottomright", legend = c("Management target", "Minimum stock size threshold"), - lty = 2, col = c("red", "darkgreen"), bty = "n" + lty = 2, col = c("darkgreen", "red"), bty = "n" ) } @@ -199,8 +199,6 @@ profile_plot <- function(mydir, rep, para, profilesummary) { ) ) ) - - subplots <- profilesummary$subplots r4ss::SSplotComparisons( summaryoutput = profilesummary, @@ -216,7 +214,7 @@ profile_plot <- function(mydir, rep, para, profilesummary) { btarg = btarg, minbthresh = thresh, plotdir = mydir, - subplots = subplots, + subplots = profilesummary$subplots, pdf = FALSE, print = TRUE, plot = FALSE, filenameprefix = paste0(para, "_trajectories_") ) From 9e7778528bc6a231023c3a026dc4930381c9cee6 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 09:08:20 -0700 Subject: [PATCH 25/28] update version number and info on backward compatibility --- README.md | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index 9e27347..372105a 100644 --- a/README.md +++ b/README.md @@ -2,12 +2,14 @@ # nwfscDiag: Diagnostic Package for West Coast Groundfish Assessments -The package provides the functionality to conduct model diagnostics for SS models. The standard diagnostic included in this package are standard required analysis for West Coast Groundfish stock assessments. The package was designed to perform model diagnostics and create plots and tables in a standardized format. The standardized approach will facilitate the use of these outputs in the assessment template approach [sa4ss](https://github.com/pfmc-assessments/sa4ss). +The package provides the functionality to conduct model diagnostics for Stock Synthesis (SS3) models. The standard diagnostic included in this package are standard required analysis for U.S. West Coast Groundfish stock assessments managed by the Pacific Fisheries Management Council. The package was designed to perform model diagnostics and create plots and tables in a standardized format. The standardized approach will facilitate the use of these outputs in the assessment template approach [sa4ss](https://github.com/pfmc-assessments/sa4ss). The diagnostics created by the package are: -- jitter runs to ensure model convergence at the MLE -- retrospective runs to examine model sensitivity to recent data -- likelihood profiles across parameters +- jitter runs to ensure model convergence at the MLE, +- retrospective runs to examine model sensitivity to recent data, and +- likelihood profiles across parameters. + +This package does not maintain backward compatibility with previous versions of Stock Synthesis. However, if needed user can download older package versions that may work with older versions (3.30.+) of Stock Synthesis. ## Installation @@ -17,12 +19,12 @@ install.packages("remotes") remotes::install_github("pfmc-assessments/nwfscDiag") ``` ## Running the code -The package depends upon a few other packages and they should be installed upon installation of the package. The dependant packages are: +The package depends upon a few other packages and they should be installed upon installation of the package. The dependent packages are: ``` install.packages('dplyr') remotes::install_github('r4ss/r4ss') ``` -A new version of r4ss package was released on July 29, 2022. This release included some significant changes that included changes to function names and function inputs. The current version of the nwfscDiag 1.1.1 package is designed to work with the latest release of r4ss. Please see release version 1.0.1 to use earlier versions of r4ss. +A new version of r4ss package was released on July 29, 2022. This release included some significant changes that included changes to function names and function inputs. The current version of the nwfscDiag 1.1.2 package is designed to work with the latest release of r4ss. Please see release version 1.0.1 to use earlier versions of r4ss. ## Reporting problems From 459f80dc611c6c92694fd5c59bb50f990ada2c07 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 14:24:59 -0700 Subject: [PATCH 26/28] rm deprecated input --- tests/testthat/test-runs.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-runs.R b/tests/testthat/test-runs.R index 1857c2e..9ce6d53 100644 --- a/tests/testthat/test-runs.R +++ b/tests/testthat/test-runs.R @@ -52,7 +52,7 @@ test_that("Do profile using the simple model", { file.exists( file.path(path, paste0(model_settings$base_name, "_profile_", - profile_para[a], "_prior_like_", get[a, "use_prior_like"]), + profile_para[a]), paste0("profile_", profile_para[a], "_quant_table.csv"))) ) } From fe412821562a5eaa354c000dadd1a1986a7ab8be Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 14:26:06 -0700 Subject: [PATCH 27/28] update function description and rm select messages --- R/profile_wrapper.R | 72 +++++++++++++++++++++--------------------- man/profile_wrapper.Rd | 35 ++++++++++++++------ 2 files changed, 61 insertions(+), 46 deletions(-) diff --git a/R/profile_wrapper.R b/R/profile_wrapper.R index d263ac7..fe021b9 100644 --- a/R/profile_wrapper.R +++ b/R/profile_wrapper.R @@ -1,22 +1,34 @@ -#' Run [r4ss::profile] based on `model_settings` +#' Run [r4ss::profile()] based on `model_settings` #' -#' Generate likelihood profiles -#' To be called from the run_diagnostics function after creating -#' the model settings using the get_settings function. -#' -#' -#' @seealso The following functions interact with `profile_wrapper`: -#' * [run_diagnostics]: calls `profile_wrapper` -#' * [r4ss::profile]: the workhorse of `profile_wrapper` that does the parameter profiles +#' Setting up the specifications, running the profile using [r4ss::profile()], +#' and generating figures and tables can be tedious, error prone, and time +#' consuming. Thus, `profile_wrapper()` aims to further decrease the work +#' needed to generate a profile that is easily included in a assessment report +#' for the Pacific Fisheries Management Council. See the See Also section for +#' information on a workflow to use this function. #' +#' @seealso +#' The following functions interact with `profile_wrapper()`: +#' * [get_settings()]: populates a list of settings for `model_settings` +#' * [run_diagnostics()]: calls `profile_wrapper()` +#' * [r4ss::profile()]: the workhorse of `profile_wrapper()` that does the +#' parameter profiles #' #' @template mydir #' @template model_settings #' -#' @author Chantel Wetzel. +#' @author Chantel Wetzel and Ian Taylor. #' @return -#' Nothing is explicitly returned from `profile_wrapper`. -#' The following objects are saved to the disk. +#' Nothing is explicitly returned from `profile_wrapper()`. +#' The following objects are saved to the disk: +#' * `*_profile_output.Rdata` +#' * `*_trajectories_*` from [r4ss::SSplotComparisons()] +#' * `piner_panel_*.png` +#' * `parameter_panel_*.png` +#' * `run_diag_warning.txt` +#' * a copy of the control file saved to `model_settings$newctlfile` +#' * `backup_oldctlfile.ss` +#' * `backup_ss.par` #' @export profile_wrapper <- function(mydir, model_settings) { @@ -160,31 +172,19 @@ profile_wrapper <- function(mydir, model_settings) { # loop over down, then up for (iprofile in 1:2) { - if (iprofile == 1) { - # subset of requested runs which are in the "low" vector - whichruns <- which(vec %in% low) - if (!is.null(model_settings$whichruns)) { - whichruns <- intersect(model_settings$whichruns, whichruns) - } - if (model_settings$verbose) { - message("Running profiles down from estimated value:", - paste(whichruns, collapse = " "), - "out of ", - length(vec)) - } - } + whichruns <- which(vec %in% if(iprofile == 1){low} else {high}) + #if (iprofile == 1) { + # whichruns <- which(vec %in% low) + # if (!is.null(model_settings$whichruns)) { + # whichruns <- intersect(model_settings$whichruns, whichruns) + # } + #} + #if (iprofile == 2) { + # whichruns <- which(vec %in% high) + if (!is.null(model_settings$whichruns)) { + whichruns <- intersect(model_settings$whichruns, whichruns) + } if (iprofile == 2) { - # subset of requested runs which are in the "low" vector - whichruns <- which(vec %in% high) - if (!is.null(model_settings$whichruns)) { - whichruns <- intersect(model_settings$whichruns, whichruns) - } - if (model_settings$verbose) { - message("Running profiles up from estimated value:", - paste(whichruns, collapse = " "), - "out of ", - length(vec)) - } # 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)) diff --git a/man/profile_wrapper.Rd b/man/profile_wrapper.Rd index 2a208f2..37ff168 100644 --- a/man/profile_wrapper.Rd +++ b/man/profile_wrapper.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/profile_wrapper.R \name{profile_wrapper} \alias{profile_wrapper} -\title{Run \link[r4ss:profile]{r4ss::profile} based on \code{model_settings}} +\title{Run \code{\link[r4ss:profile]{r4ss::profile()}} based on \code{model_settings}} \usage{ profile_wrapper(mydir, model_settings) } @@ -17,21 +17,36 @@ Typically, created using \link{get_settings} but users can create their own list as long as it has all of the necessary components for the function at hand.} } \value{ -Nothing is explicitly returned from \code{profile_wrapper}. -The following objects are saved to the disk. +Nothing is explicitly returned from \code{profile_wrapper()}. +The following objects are saved to the disk: +\itemize{ +\item \verb{*_profile_output.Rdata} +\item \verb{*_trajectories_*} from \code{\link[r4ss:SSplotComparisons]{r4ss::SSplotComparisons()}} +\item \code{piner_panel_*.png} +\item \code{parameter_panel_*.png} +\item \code{run_diag_warning.txt} +\item a copy of the control file saved to \code{model_settings$newctlfile} +\item \code{backup_oldctlfile.ss} +\item \code{backup_ss.par} +} } \description{ -Generate likelihood profiles -To be called from the run_diagnostics function after creating -the model settings using the get_settings function. +Setting up the specifications, running the profile using \code{\link[r4ss:profile]{r4ss::profile()}}, +and generating figures and tables can be tedious, error prone, and time +consuming. Thus, \code{profile_wrapper()} aims to further decrease the work +needed to generate a profile that is easily included in a assessment report +for the Pacific Fisheries Management Council. See the See Also section for +information on a workflow to use this function. } \seealso{ -The following functions interact with \code{profile_wrapper}: +The following functions interact with \code{profile_wrapper()}: \itemize{ -\item \link{run_diagnostics}: calls \code{profile_wrapper} -\item \link[r4ss:profile]{r4ss::profile}: the workhorse of \code{profile_wrapper} that does the parameter profiles +\item \code{\link[=get_settings]{get_settings()}}: populates a list of settings for \code{model_settings} +\item \code{\link[=run_diagnostics]{run_diagnostics()}}: calls \code{profile_wrapper()} +\item \code{\link[r4ss:profile]{r4ss::profile()}}: the workhorse of \code{profile_wrapper()} that does the +parameter profiles } } \author{ -Chantel Wetzel. +Chantel Wetzel and Ian Taylor. } From b9868719da56914372ce93d2d3c1c649a40d9887 Mon Sep 17 00:00:00 2001 From: Chantel Wetzel Date: Wed, 28 Jun 2023 14:26:29 -0700 Subject: [PATCH 28/28] update deprecated input text --- man/get_settings_profile.Rd | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/man/get_settings_profile.Rd b/man/get_settings_profile.Rd index 0f95438..8808007 100644 --- a/man/get_settings_profile.Rd +++ b/man/get_settings_profile.Rd @@ -26,11 +26,9 @@ get_settings_profile( real indicates bounds in the parameter space, relative indicates how far to go from the base parameter, and multiplier indicates that low and high bounds are set at x\\% above and below the base parameter.} -\item{use_prior_like}{Deprecated Option to include or exclude the prior likelihoods in the likelihood -profiles. A value of 0 corresponds to the Stock Sythesis were this would exclude -and a value of 1 would include the prior likelihood contribution. Parameters with priors used for estimation -(e.g., natural mortality, steepness) are often profiled across and including or excluding the prior likelihood -contribution may be wanted in specific instances. The default setting excludes the prior likelihood contributions.} +\item{use_prior_like}{Deprecated: The use_prior_like input is no longer needed since r4ss now +automatically plots the likelihood profile with and without any parameter prior likelihood contributions +regardless of the setting in the user starter file.} } \value{ A matrix of low, high, and step size values for the default parameters