Merge pull request #40 from pfmc-assessments/github_actions

Add GitHub actions and vignette

.github/.gitignore

@@ -0,0 +1,2 @@
+*.html
+*.rds

.github/workflows/call-doc-and-style-r.yml

@@ -0,0 +1,12 @@
+# document and style R code using a reusable workflow
+name: call-doc-and-style-r
+# on specifies the build triggers. See more info at https://docs.github.com/en/actions/learn-github-actions/events-that-trigger-workflows
+on:
+  # workflow_dispatch requires pushing a button to run the workflow manually. uncomment the following line to add:
+  workflow_dispatch:
+  push:
+    branches: [main]
+jobs:
+  call-workflow:
+    uses: nmfs-fish-tools/ghactions4r/.github/workflows/doc-and-style-r.yml@main
+

.github/workflows/call-update-pkgdown.yml

@@ -0,0 +1,10 @@
+on:
+  workflow_dispatch:
+  push:
+    branches: [main]
+    tags: ['*']
+
+name: call-update-pkgdown
+jobs:
+  call-workflow:
+    uses: nmfs-fish-tools/ghactions4r/.github/workflows/update-pkgdown.yml@main

R/get_jitter_quants.R

@@ -20,49 +20,51 @@ get_jitter_quants <- function(mydir, model_settings, output) {
   est <- output[["est"]]
   profilesummary <- output[["profilesummary"]]
 
-  outputs <- output$profilemodels
+  outputs <- output[["profilemodels"]]
   quants <- lapply(outputs, "[[", "derived_quants")
   status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status")
-  bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")])
+  bounds <- apply(status, 2, function(x) rownames(outputs[[1]][["parameters"]])[x %in% c("LO", "HI")])
   out <- data.frame(
     "run" = gsub("replist", "", names(outputs)),
     "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_", profilesummary$endyrs[1]), "Value"),
+    "SBfinal" = sapply(quants, "[[", paste0("SSB_", profilesummary[["endyrs"]][1]), "Value"),
     "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))),
     "Lowest NLL" = ifelse(min(like) == like, "Best Fit", 0),
     stringsAsFactors = FALSE
   )
 
   # Write a md file to be included in a stock assessment document
   # Text was pirated from @chantelwetzel-noaa's 2021 dover assessment
-  file_md <- file.path(jitter_dir, "model-results-jitter.md")
-  sink(file_md)
-  on.exit(sink(), add = TRUE)
-  cat(
-    sep = "",
-    "Model convergence was in part based on starting the minimization process ",
-    "from dispersed values of the maximum likelihood estimates to determine if the ",
-    "estimation routine results in a smaller likelihood.\n",
-    "Starting parameters were jittered using the built-in functionality of ",
-    "Stock Synthesis, where you specify a jitter fraction.\n",
-    "Here we used a jitter fraction of ",
-    round(model_settings$jitter_fraction, 2), " and the jittering was repeated ",
-    xfun::numbers_to_words(model_settings$Njitter), " times.\n",
-    "A better, i.e., lower negative log-likelihood, fit was ",
-    ifelse(
-      sum(like - est < 0) == 0,
-      "not found",
-      paste0("found for ", xfun::numbers_to_words(sum(like - est < 0)), " fits")
-    ), ".\n",
-    "Several models resulted in similar log-likelihood values ",
-    "with little difference in the overall model estimates, ",
-    "indicating a relatively flat likelihood surface around the maximum likelihood estimate.\n",
-    "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.\n"
+  utils::write.csv(
+    x = data.frame(
+      caption = paste(
+        sep = "",
+        "Model convergence was in part based on starting the minimization process ",
+        "from dispersed values of the maximum likelihood estimates to determine if the ",
+        "estimation routine results in a smaller likelihood.",
+        "Starting parameters were jittered using the built-in functionality of ",
+        "Stock Synthesis, where you specify a jitter fraction.",
+        "Here we used a jitter fraction of ",
+        round(model_settings[["jitter_fraction"]], 2), " and the jittering was repeated ",
+        xfun::numbers_to_words(model_settings[["Njitter"]]), " times.",
+        "A better, i.e., lower negative log-likelihood, fit was ",
+        dplyr::if_else(
+          sum(like - est < 0) == 0,
+          true = "not found",
+          false = paste0("found for ", xfun::numbers_to_words(sum(like - est < 0)), " fits")
+        ),
+        "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."),
+      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"))
+    ),
+    file = file.path(jitter_dir, "jitterfigures4doc.csv"),
+    row.names = FALSE
   )
 
   # write tables

R/get_param_values.R

@@ -16,36 +16,36 @@
 get_param_values <- function(mydir, para = NULL, vec, summary) {
 
   x <- summary
-  n <- x$n
-  endyr <- x$endyrs[1] + 1
+  n <- x[["n"]]
+  endyr <- x[["endyrs"]][1] + 1
   out <- data.frame(
-    totlikelihood = as.numeric(x$likelihoods[x$likelihoods$Label == "TOTAL", 1:n]),
-    surveylike = as.numeric(x$likelihoods[x$likelihoods$Label == "Survey", 1:n]),
-    discardlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Discard", 1:n]),
-    lengthlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Length_comp", 1:n]),
-    agelike = as.numeric(x$likelihoods[x$likelihoods$Label == "Age_comp", 1:n]),
-    recrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Recruitment", 1:n]),
-    forerecrlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Forecast_Recruitment", 1:n]),
-    priorlike = as.numeric(x$likelihoods[x$likelihoods$Label == "Parm_priors", 1:n]),
-    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]),
-    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_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_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]),
-    cv1male = as.numeric(x$pars[grep("young_Mal_GP_1", x$pars$Label), 1:n]),
-    cv2male = as.numeric(x$pars[grep("old_Mal_GP_1", x$pars$Label), 1:n]),
+    totlikelihood = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "TOTAL", 1:n]),
+    surveylike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Survey", 1:n]),
+    discardlike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Discard", 1:n]),
+    lengthlike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Length_comp", 1:n]),
+    agelike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Age_comp", 1:n]),
+    recrlike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Recruitment", 1:n]),
+    forerecrlike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Forecast_Recruitment", 1:n]),
+    priorlike = as.numeric(x[["likelihoods"]][x[["likelihoods"]][["Label"]] == "Parm_priors", 1:n]),
+    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]),
+    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_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_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]),
+    cv1male = as.numeric(x[["pars"]][grep("young_Mal_GP_1", x[["pars"]][["Label"]]), 1:n]),
+    cv2male = as.numeric(x[["pars"]][grep("old_Mal_GP_1", x[["pars"]][["Label"]]), 1:n]),
     stringsAsFactors = FALSE
   )
 

R/get_retro_quants.R

@@ -23,7 +23,7 @@
 #'   inside of an environment with `results = "asis"`
 #'   to include a table of Mohn's rho values in a document.
 #'
-#'   `sa4ss::read_child(file.path(paste0(params$model, "_retro"), "mohnsrho.tex"))`
+#'   `sa4ss::read_child(file.path(paste0(params[["model"]], "_retro"), "mohnsrho.tex"))`
 #'
 #'
 #' @export
@@ -47,7 +47,7 @@ get_retro_quants <- function(mydir,  model_settings, output) {
       caption = paste(
         "Retrospective patterns for",
         c("spawning stock biomass (\\emph{SSB})", "fraction unfished"),
-        "when up to", xfun::numbers_to_words(max(abs(model_settings$retro_yr))),
+        "when up to", xfun::numbers_to_words(max(abs(model_settings[["retro_yrs"]]))),
         "years of data were removed from the base model.",
         "Mohn's rho (Mohn, 1999) values were",
         "recalculated for each peel given the removal of another year of data.",

R/get_settings.R

@@ -68,24 +68,24 @@ get_settings <- function(settings = NULL, verbose = FALSE) {
     subplots = c(1, 3)
   )
 
-  Settings_all$profile_details <- get_settings_profile()
+  Settings_all[["profile_details"]] <- get_settings_profile()
 
   need <- !names(Settings_all) %in% names(settings)
   Settings_all <- c(settings, Settings_all[need])
 
   # Check some items
-  if (!is.null(Settings_all$profile_details)) {
-    if (length(Settings_all$profile_details[is.na(Settings_all$profile_details)]) > 0) {
+  if (!is.null(Settings_all[["profile_details"]])) {
+    if (length(Settings_all[["profile_details"]][is.na(Settings_all[["profile_details"]])]) > 0) {
       cli::cli_abort(
         "Missing entry in the get_settings_profile data frame."
       )
     }
-    if (!is.numeric(Settings_all$profile_details$low) &
-      !is.numeric(Settings_all$profile_details$high) &
-      !is.numeric(Settings_all$profile_details$step_size)) {
+    if (!is.numeric(Settings_all[["profile_details"]][["low"]]) &
+      !is.numeric(Settings_all[["profile_details"]][["high"]]) &
+      !is.numeric(Settings_all[["profile_details"]][["step_size"]])) {
       cli::cli_abort("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) {
+    if (sum(!Settings_all[["profile_details"]][["param_space"]] %in% c("real", "relative", "multiplier")) > 0) {
       cli::cli_abort("The param_space column should be either real or relative in the get_settings_profile data frame.")
     }
   }

R/get_summary.R

@@ -22,7 +22,7 @@ get_summary <- function(mydir, para, vec, profilemodels, profilesummary) {
   outputs <- profilemodels
   quants <- lapply(outputs, "[[", "derived_quants")
   status <- sapply(sapply(outputs, "[[", "parameters", simplify = FALSE), "[[", "Status")
-  bounds <- apply(status, 2, function(x) rownames(outputs[[1]]$parameters)[x %in% c("LO", "HI")])
+  bounds <- apply(status, 2, function(x) rownames(outputs[[1]][["parameters"]])[x %in% c("LO", "HI")])
 
   out <- data.frame(
     "run" = gsub("replist", "", names(outputs)),
@@ -31,9 +31,8 @@ get_summary <- function(mydir, para, vec, 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 + 1), "Value"),
-    "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]]$endyr + 1), "Value"),
-    # "Fmsy" = sapply(quants, "[[", "annF_MSY", "Value"),
+    "SBfinal" = sapply(quants, "[[", paste0("SSB_", outputs[[1]][["endyr"]] + 1), "Value"),
+    "Deplfinal" = sapply(quants, "[[", paste0("Bratio_", outputs[[1]][["endyr"]] + 1), "Value"),
     "Nparsonbounds" = apply(status, 2, function(x) sum(x %in% c("LO", "HI"))),
     stringsAsFactors = FALSE
   )

R/plot_jitter.R

@@ -19,11 +19,11 @@ plot_jitter <- function(mydir, model_settings, output) {
   est <- output[["est"]]
   profilesummary <- output[["profilesummary"]]
 
-  ymax <- as.numeric(stats::quantile(unlist(profilesummary$likelihoods[1, keys]), 0.80))
+  ymax <- as.numeric(stats::quantile(unlist(profilesummary[["likelihoods"]][1, keys]), 0.80))
   ymin <- min(like - est) + 1
   ylab <- "Change in negative log-likelihood"
   xlab <- "Iteration"
-  pngfun(wd = jitter_dir, file = paste0("Jitter_", model_settings$jitter_fraction, ".png"), h = 12, w = 9)
+  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,
@@ -48,7 +48,7 @@ plot_jitter <- function(mydir, model_settings, output) {
   )
 
   if (ymax > 100) {
-    pngfun(wd = jitter_dir, file = paste0("Jitter_Zoomed_SubPlot_", model_settings$jitter_fraction, ".png"), h = 12, w = 9)
+    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,

R/plot_retro.R

@@ -31,18 +31,19 @@ plot_retro <- function(mydir,  model_settings, output) {
     legendlabels = c(
       "Base Model",
       sprintf("Data %.0f year%s",
-              model_settings$retro_yrs,
-              ifelse(abs(model_settings$retro_yrs) == 1, "", "s")
+              model_settings[["retro_yrs"]],
+              ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s")
       )
     ),
-    btarg = model_settings$btarg,
-    minbthresh = model_settings$minbthresh,
+    btarg = model_settings[["btarg"]],
+    minbthresh = model_settings[["minbthresh"]],
     ylimAdj = 1.2,
     plotdir = retro_dir,
     legendloc = "topright",
     print = TRUE,
     plot = FALSE,
-    pdf = FALSE
+    pdf = FALSE,
+    verbose = model_settings[["verbose"]]
   )
   savedplotinfo <- mapply(
     FUN = r4ss::SSplotComparisons,
@@ -52,9 +53,10 @@ plot_retro <- function(mydir,  model_settings, output) {
       legendloc = "topleft",
       plotdir = retro_dir,
       ylimAdj = 1.2,
-      btarg = model_settings$btarg,
-      minbthresh = model_settings$minbthresh,
-      print = TRUE, plot = FALSE, pdf = FALSE
+      btarg = model_settings[["btarg"]],
+      minbthresh = model_settings[["minbthresh"]],
+      print = TRUE, plot = FALSE, pdf = FALSE,
+      verbose = model_settings[["verbose"]]
     ),
     subplot = c(8, 10),
     legendlabels = lapply(
@@ -63,8 +65,8 @@ plot_retro <- function(mydir,  model_settings, output) {
         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"),
+                  model_settings[["retro_yrs"]],
+                  ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s"),
                   rhosall[rownames(rhosall) == x, ]
           )
         )
@@ -77,19 +79,20 @@ plot_retro <- function(mydir,  model_settings, output) {
     legendlabels = c(
       "Base Model",
       sprintf("Data %.0f year%s",
-              model_settings$retro_yrs,
-              ifelse(abs(model_settings$retro_yrs) == 1, "", "s")
+              model_settings[["retro_yrs"]],
+              ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s")
       )
     ),
-    btarg = model_settings$btarg,
-    minbthresh = model_settings$minbthresh,
+    btarg = model_settings[["btarg"]],
+    minbthresh = model_settings[["minbthresh"]],
     subplot = c(2, 4),
     ylimAdj = 1.2,
     plotdir = retro_dir,
     legendloc = "topright",
     print = TRUE,
     plot = FALSE,
-    pdf = FALSE
+    pdf = FALSE,
+    verbose = model_settings[["verbose"]]
   )
   savedplotinfo <- mapply(
     FUN = r4ss::SSplotComparisons,
@@ -99,9 +102,10 @@ plot_retro <- function(mydir,  model_settings, output) {
       legendloc = "topright",
       ylimAdj = 1.2,
       plotdir = retro_dir,
-      btarg = model_settings$btarg,
-      minbthresh = model_settings$minbthresh,
-      print = TRUE, plot = FALSE, pdf = FALSE
+      btarg = model_settings[["btarg"]],
+      minbthresh = model_settings[["minbthresh"]],
+      print = TRUE, plot = FALSE, pdf = FALSE,
+      verbose = model_settings[["verbose"]]
     ),
     subplot = c(2, 4),
     legendlabels = lapply(
@@ -110,8 +114,8 @@ plot_retro <- function(mydir,  model_settings, output) {
         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"),
+                  model_settings[["retro_yrs"]],
+                  ifelse(abs(model_settings[["retro_yrs"]]) == 1, "", "s"),
                   rhosall[rownames(rhosall) == x, ]
           )
         )
@@ -178,7 +182,7 @@ plot_retro <- function(mydir,  model_settings, output) {
   df_out <- NULL
   y <- years
   for (a in 1:n){
-    col_name <- paste0("per_diff_model", 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){
       df_out[["model"]][df_out[["model"]] == col_name] <- "Base Model"