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 <chantel.wetzel@noaa.gov>
-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)
})