Update ss3 scripts

pbs-assess · Jul 16, 2024 · 2ea50a7 · 2ea50a7
1 parent 41d273b
commit 2ea50a7
Show file tree

Hide file tree

Showing 3 changed files with 84 additions and 39 deletions.
diff --git a/ss3/02-outside-ss3-r4ss.R b/ss3/02-outside-ss3-r4ss.R
@@ -1,7 +1,7 @@
 
 model_dir <- "m37_IPHC_g4"
 ss_home <- here::here("ss3")
-ss_home <- "C:/users/qhuynh/Desktop/dogfish"
+ss_home <- "C:/users/quang/Documents/dogfish"
 
 #### Fit ss3 model version 3.30.22.1 with custom compilation that fixes lognormal prior density function
 fit_ss3 <- function(model_dir = "model1",
@@ -34,13 +34,15 @@ fit_ss3 <- function(model_dir = "model1",
 # fit_ss3(model_dir, hessian = FALSE, ss_home = ss_home, max_phase = 10L)
 
 
-mods <- c("A1", "A2_USgrowth", "A3_highmat", "A4_USgrowth_highmat", "A5_exHS", "A6_IPHC", "A7_SYN", "B1_2005")
-snowfall::sfInit(parallel = TRUE, cpus = length(mods))
+mods <- c(#"A1",
+          "A2_USgrowth", "A3_highmat", "A4_USgrowth_highmat", "A5_exHS", "A6_IPHC+CPUE", "A7_SYNonly", "A8_HBLLonly", "A9_lowM",
+          "B1_1990inc", "B2_2010step", "B3_2005step", "B4_1990inc_lowM", "B5_2010step_lowM")
+snowfall::sfInit(parallel = TRUE, cpus = 6)
 snowfall::sfLapply(mods, fit_ss3, hessian = TRUE, ss_home = ss_home)
 snowfall::sfStop()
 
 # Load r4ss list
-model_dir <- "B1_2005"
+model_dir <- "A1"
 covar <- TRUE
 replist <- r4ss::SS_output(
   file.path(ss_home, model_dir),

diff --git a/ss3/03-outside-ss3-figures.R b/ss3/03-outside-ss3-figures.R
@@ -5,56 +5,86 @@ source("ss3/ss3_functions.R")
 
 # Compare SS models
 #ss_home <- here::here("ss3")
-ss_home <- "C:/users/qhuynh/Desktop/dogfish"
+ss_home <- "C:/users/quang/Documents/dogfish"
 
-# Set A models with growth and maturity scenarios
-mods <- c("A1", "A2_USgrowth", "A3_highmat", "A4_USgrowth_highmat")
-model_name <- c("(A1) BCgrowth", "(A2) USgrowth", "(A3) BCgrowth, high mat", "(A4) USgrowth, high mat")
 
-fig_dir <- "figs/ss3/set_a_mat"
+# Specify which set of plots to generate here
+set_to_plot <- "growth"
 
-# Set A models jackknifing indices
-mods <- c("A1", "A6_IPHC+CPUE", "A7_SYNonly", "A8_HBLLonly")
-model_name <- c("(A1) BCgrowth + all indices", "(A6) IPHC + CPUE only", "(A6) SYN only", "(A7) HBLL only")
+set_to_plot <- match.arg(set_to_plot, choices = c("growth", "index", "Mchange", "lowbaseM"))
 
-fig_dir <- "figs/ss3/set_a_ind"
+if (set_to_plot == "growth") {
 
+  # Set A models with growth and maturity scenarios
+  mods <- c("A1", "A2_USgrowth", "A3_highmat", "A4_USgrowth_highmat")
+  model_name <- c("(A1) BCgrowth", "(A2) USgrowth", "(A3) BCgrowth, high mat", "(A4) USgrowth, high mat")
 
-# Set B (M change)
-mods <- c("A1", "B1_1990inc", "B2_2010step", "B3_2005step")
-model_name <- c("(A1) Constant M", "(B1) Linear increase 1990", "(B2) Step increase 2010", "(B3) Step increase 2005")
+  fig_dir <- "figs/ss3/set_a_mat"
 
-fig_dir <- "figs/ss3/set_b"
+  multi_rep <- lapply(mods, function(x) {
+    r4ss::SS_output(file.path(ss_home, x),
+                    verbose = FALSE,
+                    printstats = FALSE,
+                    hidewarn = TRUE)
+  })
+  saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_mat.rds"))
 
+} else if (set_to_plot == "index") {
 
+  # Set A models jackknifing indices
+  mods <- c("A1", "A6_IPHC+CPUE", "A7_SYNonly", "A8_HBLLonly")
+  model_name <- c("(A1) BCgrowth + all indices", "(A6) IPHC + CPUE only", "(A6) SYN only", "(A7) HBLL only")
 
+  fig_dir <- "figs/ss3/set_a_ind"
 
+  multi_rep <- lapply(mods, function(x) {
+    r4ss::SS_output(file.path(ss_home, x),
+                    verbose = FALSE,
+                    printstats = FALSE,
+                    hidewarn = TRUE)
+  })
+  saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_jackknife.rds"))
 
+} else if (set_to_plot == "Mchange") {
 
+  # Set B (M change)
+  mods <- c("A1", "B1_1990inc", "B2_2010step", "B3_2005step")
+  model_name <- c("(A1) Constant M", "(B1) Linear increase 1990", "(B2) Step increase 2010", "(B3) Step increase 2005")
 
+  fig_dir <- "figs/ss3/set_b"
 
-# Read and save r4ss report lists ----
-.ggsave <- function(filename, ...) {
-  filename <- file.path(fig_dir, filename)
-  ggsave(filename, ...)
-}
+  multi_rep <- lapply(mods, function(x) {
+    r4ss::SS_output(file.path(ss_home, x),
+                    verbose = FALSE,
+                    printstats = FALSE,
+                    hidewarn = TRUE)
+  })
+  saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_Minc.rds"))
 
-multi_rep <- lapply(mods, function(x) {
-  r4ss::SS_output(file.path(ss_home, x),
-                  verbose = FALSE,
-                  printstats = FALSE,
-                  hidewarn = TRUE)
-})
+} else if (set_to_plot == "lowbaseM") {
 
-saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_mat.rds"))
-multi_rep <- readRDS(file = file.path(ss_home, "multi_rep_mat.rds"))
+  # Compare low M (combination of A and B models)
+  mods <- c("A1", "A9_lowM", "B1_1990inc", "B4_1990inc_lowM", "B2_2010step", "B5_2010step_lowM")
+  model_name <- c("(A1) M = 0.074", "(A9) M = 0.05",
+                  "(B1) M = 0.074, inc. 1990", "(B4) M = 0.05, inc. 1990",
+                  "(B2) M = 0.074, Step increase 2010",  "(B5) M = 0.05, Step increase 2010")
 
-saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_jackknife.rds"))
-multi_rep <- readRDS(file = file.path(ss_home, "multi_rep_jackknife.rds"))
+  fig_dir <- "figs/ss3/lowM"
 
-saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_Minc.rds"))
-multi_rep <- readRDS(file = file.path(ss_home, "multi_rep_Minc.rds"))
+  multi_rep <- lapply(mods, function(x) {
+    r4ss::SS_output(file.path(ss_home, x),
+                    verbose = FALSE,
+                    printstats = FALSE,
+                    hidewarn = TRUE)
+  })
+  saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_lowM.rds"))
+}
 
+# Read and save r4ss report lists ----
+.ggsave <- function(filename, ...) {
+  filename <- file.path(fig_dir, filename)
+  ggsave(filename, ...)
+}
 
 # Check parameters and bounds
 multi_rep[[1]]$estimated_non_dev_parameters %>% View()

diff --git a/ss3/05-outside-ss3-retro.R b/ss3/05-outside-ss3-retro.R
@@ -1,11 +1,10 @@
 
 
 ss_home <- here::here("ss3")
-ss_home <- "C:/users/qhuynh/Desktop/dogfish"
-model_dir <- c("A1", "B2", "C1")
+ss_home <- "C:/users/quang/Documents/dogfish"
+model_dir <- c("A1", "A9_lowM", "B2_2010step")
 
 # Do retro
-model_dir <- "B2_2010step"
 ypeel <- seq(0, -7, -1)
 for(i in model_dir) {
 
@@ -18,13 +17,27 @@ for(i in model_dir) {
 
 }
 
+library(snowfall)
+sfInit(parallel = TRUE, cpus = length(model_dir))
+sfExport(list = c("ss_home", "ypeel"))
+sfLapply(model_dir, function(i) {
+  r4ss::retro(
+    file.path(ss_home, i),
+    years = ypeel,
+    exe = "ss",
+    extras = c("-nohess modelname ss")
+  )
+})
+sfStop()
+
+
 # Make figures
 source("ss3/ss3_functions.R")
 
 ret <- r4ss::SSgetoutput(
-  dirvec = file.path(ss_home, model_dir[1], "retrospectives", paste0("retro", ypeel))
+  dirvec = file.path(ss_home, model_dir[3], "retrospectives", paste0("retro", ypeel))
 )
 g <- SS3_retro(ret)
-#ggsave("figs/ss3/ret_A1.png", g, height = 8, width = 5)
+ggsave("figs/ss3/ret_A1.png", g, height = 8, width = 5)
 ggsave("figs/ss3/ret_B2.png", g, height = 8, width = 5)