Fix figs

ss3/02-outside-ss3-r4ss.R

@@ -41,6 +41,12 @@ mods <- c("A1", "A0",
   "A15_100discard",
   "B1_1990inc", "B2_2010step", "B3_2005step", "B4_1990inc_lowM", "B5_2010step_lowM")
 
+# for (i in seq_along(mods)) {
+#   d <- r4ss::SS_readctl(paste0("ss3/", mods[i], "/control.ss_new"))
+#   print(mods[i])
+#   print(any(grepl("extraSD", row.names(d$Q_parms))))
+# }
+
 # # Make sure starter matches...
 tocopy <- seq_along(mods)[-2] # copying 2
 for (i in tocopy) {

ss3/03-outside-ss3-figures.R

@@ -27,8 +27,6 @@ dir.create("figs/ss3/set_a_ind", showWarnings = FALSE)
 # Specify which set of plots to generate here
 # set_to_plot <- c("growth", "index", "M")[3]
 
-base_model <- 2
-
 for (set_to_plot in c("growth", "index", "M", "zfrac")) {
 
   set_to_plot <- match.arg(set_to_plot, choices = c("growth", "index", "M", "zfrac"))
@@ -48,7 +46,7 @@ for (set_to_plot in c("growth", "index", "M", "zfrac")) {
   if (set_to_plot == "growth") {
 
     # Set A models with growth and maturity scenarios
-    mods <- c("A0", "A2_USgrowth", "A3_highmat", "A4_USgrowth_highmat", "A14_lowdiscard", "A5_highdiscard", "A15-100discard")
+    mods <- c("A0", "A2_USgrowth", "A3_highmat", "A4_USgrowth_highmat", "A14_lowdiscard", "A5_highdiscard", "A15_100discard")
 
     model_name <- c("(A0) BC growth\n(base)", "(A2) US growth", "(A3) BC growth,\nhigh maturity", "(A4) USgrowth,\nhigh maturity", "(A14) Low discard\nmortality", "(A5) High discard\nmortality", "(A15) 100% discard mortality")
 
@@ -63,7 +61,7 @@ for (set_to_plot in c("growth", "index", "M", "zfrac")) {
     saveRDS(multi_rep, file = file.path(ss_home, "multi_rep_mat.rds"))
 
     # Save the model parameters in model 1 for Res. Doc.
-    multi_rep[[base_model]]$parameters %>%
+    multi_rep[[1]]$parameters %>%
       mutate(`Estimated?` = ifelse(Phase > 0, "Yes", "Fixed")) %>%
       select(Label, `Estimated?`, Value, Parm_StDev) %>%
       readr::write_csv(file = "tables/ss3_par.csv")
@@ -316,7 +314,7 @@ for (set_to_plot in c("growth", "index", "M", "zfrac")) {
     labs(x = "Year", y = "Index", colour = "Model") +
     guides(colour = guide_legend(ncol = 2)) +
     theme(panel.spacing = unit(0, "in"), legend.position = "bottom")
-  .ggsave("index_fit.png", g, height = 4, width = 6)
+  .ggsave("index_fit.png", g, height = 5, width = 6)
 
   # Plot SR
   g <- SS3_SR(multi_rep, model_name) +
@@ -456,12 +454,12 @@ for (set_to_plot in c("growth", "index", "M", "zfrac")) {
 
     # Length comp residual - histogram
     g <- SS3_compresid(multi_rep[[1]], model_name[1], fleet = fleet_int, figure = "histogram") +
-      ggtitle(model_name[base_model]) +
+      ggtitle(model_name[1]) +
       theme(panel.spacing = unit(0, "in"))
     .ggsave("len_comp_hist_A1.png", g, height = 4, width = 6)
 
     # Numbers at age
-    g <- SS3_N(multi_rep[base_model], model_name[1], age = seq(0, 50, 10))
+    g <- SS3_N(multi_rep[1], model_name[1], age = seq(0, 50, 10))
     g2 <- g$data %>%
       mutate(variable = paste("Age", variable)) %>%
       #mutate(value = value/max(value), .by = c(variable, Sex, scen)) %>%
@@ -471,11 +469,11 @@ for (set_to_plot in c("growth", "index", "M", "zfrac")) {
       expand_limits(y = 0) +
       facet_grid(vars(variable), vars(Sex), scales = "free_y") +
       labs(x = "Year", y = "Estimated abundance") +
-      ggtitle(model_name[base_model])
+      ggtitle(model_name[1])
     .ggsave("N_age_A1.png", g2, height = 6, width = 5)
 
     # Compare with numbers at length
-    g <- SS3_N(multi_rep[base_model], model_name[base_model], type = "length", len = seq(50, 115, 15))
+    g <- SS3_N(multi_rep[1], model_name[1], type = "length", len = seq(50, 115, 15))
     g2 <- g$data %>%
       mutate(variable = paste(variable, "cm") %>% factor(levels = paste(seq(50, 115, 15), "cm"))) %>%
       ggplot(aes(Yr, value)) +
@@ -488,15 +486,15 @@ for (set_to_plot in c("growth", "index", "M", "zfrac")) {
     .ggsave("N_len_A1.png", g2, height = 6, width = 5)
 
     # Exploitation and apical F
-    g <- SS3_apicalF(multi_rep[[base_model]])
+    g <- SS3_apicalF(multi_rep[[1]])
     g$facet$params$ncol <- 3
     .ggsave("apicalF_fleet_A1.png", g, height = 6, width = 8)
 
-    g <- SS3_apicalF(multi_rep[[base_model]], FALSE)
+    g <- SS3_apicalF(multi_rep[[1]], FALSE)
     .ggsave("apicalF_sex_A1.png", g, height = 2, width = 4)
 
     # Plot annual selectivity
-    g <- SS3_selannual(multi_rep[[base_model]])
+    g <- SS3_selannual(multi_rep[[1]])
     .ggsave("sel_annual_A1.png", g, height = 4, width = 6)
 
   }

ss3/09-outside-ss3-forecast.R

@@ -136,11 +136,11 @@ mods <- mods[keep]
 model_name <- model_name[keep]
 
 length(mods)
-(tacs <- seq(0, 1500, by = 200))
+(tacs <- seq(0, 1500, by = 100))
 
 torun <- expand.grid(model = mods, catch = tacs)
 nrow(torun)
-plan(multicore, workers =  70)
+plan(multicore, workers =  60)
 
 # out3 <- filter(torun, catch %in% c(0, 200), model == "A5_highdiscard") |>
 #   purrr::pmap(run_projection, hessian = F)
@@ -215,6 +215,9 @@ temp |> filter(catch %in% tacs[seq(1, 1e2, 2)]) |>
   make_proj_by_model()
 ggsave("figs/ss3/refpts/proj-facet-model.png", width = 8.5, height = 6.5)
 
+cols <- c("grey10", RColorBrewer::brewer.pal(12L, "Paired"))
+names(cols) <- model_name[!grepl("^\\(B", model_name)]
+
 make_proj_by_catch_level <- function(dat, ylab = "S / S<sub>0</sub>") {
   dat |>
   # filter(catch %in% tacs[seq(1, 1e2, 3)]) |>
@@ -227,8 +230,8 @@ make_proj_by_catch_level <- function(dat, ylab = "S / S<sub>0</sub>") {
     geom_ribbon(fill = "grey70", alpha = 0.7, colour = NA) +
     geom_line() +
     scale_x_continuous(breaks = seq(1960, 2090, 20)) +
-    scale_colour_brewer(palette = "Paired") +
-    scale_colour_brewer(palette = "Paired") +
+    scale_colour_manual(values = cols) +
+    # scale_colour_brewer(palette = "Paired") +
     annotate(
       "rect",
       xmin = 2024, xmax = max(x$year, na.rm = TRUE),
@@ -247,7 +250,7 @@ make_proj_by_catch_level <- function(dat, ylab = "S / S<sub>0</sub>") {
     theme(axis.title = ggtext::element_markdown())
 }
 
-temp |> filter(catch %in% tacs[seq(1, 1e2, 2)]) |>
+temp |> filter(catch %in% tacs[seq(1, 12, 2)]) |>
   filter(!grepl("B", model)) |> make_proj_by_catch_level()
 ggsave("figs/ss3/refpts/proj-facet-catch.png", width = 9, height = 4.5)
 
@@ -317,9 +320,9 @@ tempF |>
   geom_ribbon(fill = "grey70", alpha = 0.7, colour = NA) +
   geom_line() +
   scale_x_continuous(breaks = seq(1960, 2090, 20)) +
-  scale_colour_brewer(palette = "Paired") +
+  scale_colour_manual(values = cols) +
+  # scale_colour_brewer(palette = "Paired") +
   # scale_colour_viridis_d(direction = -1, option = "C") +
-  scale_colour_brewer(palette = "Paired") +
   annotate(
     "rect",
     xmin = 2024, xmax = max(x$year, na.rm = TRUE),
@@ -442,4 +445,4 @@ if (FALSE) {
 
 f <- list.files("ss3", pattern = "-forecast-", full.names = T)
 f
-# x <- sapply(f, unlink, recursive = T, force = T)
+x <- sapply(f, unlink, recursive = T, force = T)

values/ref-pts.tex

@@ -11,5 +11,5 @@
 \newcommand{\A02024Yield}{193}
 \newcommand{\BaseNextYrYield}{193}
 \newcommand{\BaseNextYrYieldCI}{174--212}
-\newcommand{\EnsNextYrYield}{133}
+\newcommand{\EnsNextYrYield}{139}
 \newcommand{\EnsNextYrYieldCI}{10--307}