From b1161a8337b6431632e95056cd1593ca6b20b096 Mon Sep 17 00:00:00 2001
From: LaTreese Denson <70522445+latreesedenson-NOAA@users.noreply.github.com>
Date: Mon, 20 May 2024 14:17:47 -0400
Subject: [PATCH 01/10] Fixes discards for multi areas

-added which() functions to subset the discard fleets
-summed over areas to remove zeros in discard timeseries
---
 R/MS_functions.R |  7 +++++++
 R/RatioBiasEM.R  | 21 ++++++++++++++-------
 R/parse_MS.R     |  3 +--
 R/runSSMSE.R     |  3 +--
 4 files changed, 23 insertions(+), 11 deletions(-)

diff --git a/R/MS_functions.R b/R/MS_functions.R
index 1284545..1bb43f0 100644
--- a/R/MS_functions.R
+++ b/R/MS_functions.R
@@ -331,6 +331,13 @@ get_EM_catch_df <- function(EM_dir, dat) {
       }
     }
     dis_df <- do.call("rbind", dis_df_list)
+    dis_df <- dis_df %>%
+      dplyr::group_by(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]]) %>%
+      dplyr::summarise(catch = sum(.data[["Discard"]])) %>%
+      merge(se, all.x = TRUE, all.y = FALSE) %>%
+      dplyr::ungroup() %>%
+      dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
+    dis_df <- as.data.frame(dis_df)
   } else {
     dis_df <- NULL
   }
diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index ea4aee0..6021e33 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -268,6 +268,13 @@ get_RatioEM_catch_df<-function(EM_dir, dat, dat_yrs,
       }
     }
     dis_df <- do.call("rbind", dis_df_list)
+    dis_df <- dis_df %>%
+      dplyr::group_by(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]]) %>%
+      dplyr::summarise(catch = sum(.data[["Discard"]])) %>%
+      merge(se, all.x = TRUE, all.y = FALSE) %>%
+      dplyr::ungroup() %>%
+      dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
+    dis_df <- as.data.frame(dis_df)
   } else {
     dis_df <- NULL
   }
@@ -437,12 +444,12 @@ RatioBiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = F
     tmp_discards <- base::merge(new_OM_catch_list$discards, sample_struct$EM2OMdiscard_bias, all.x=TRUE)
     new_OM_catch_list$discards$Discard <- new_OM_catch_list$discards$Discard * tmp_discards$bias 
   }
-  
   if(2 %in% OM_dat$fleetinfo$type){ # if bycatch fleet -- remove from catch list and keep only bycatch fleets in catch_F list
     byc_f <- as.numeric(row.names(OM_dat$fleetinfo[which(OM_dat$fleetinfo$type==2),]))
-    new_OM_catch_list$catch<- new_OM_catch_list$catch[new_OM_catch_list$catch$fleet!=byc_f,]
-    new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[new_OM_catch_list$catch$fleet!=byc_f,]
-    new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[new_OM_catch_list$catch_F$fleet==byc_f,]
+    
+    new_OM_catch_list$catch<- new_OM_catch_list$catch[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),]
+    new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),]
+    new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[which(is.element(new_OM_catch_list$catch_F$fleet,byc_f)),]
   } else{
     new_OM_catch_list$catch_F <- NULL
   }
@@ -714,9 +721,9 @@ BiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE,
   
   if(2 %in% OM_dat$fleetinfo$type){ # if bycatch fleet -- remove from catch list and keep only bycatch fleets in catch_F list
     byc_f <- as.numeric(row.names(OM_dat$fleetinfo[which(OM_dat$fleetinfo$type==2),]))
-    new_OM_catch_list$catch<- new_OM_catch_list$catch[new_OM_catch_list$catch$fleet!=byc_f,]
-    new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[new_OM_catch_list$catch$fleet!=byc_f,]
-    new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[new_OM_catch_list$catch_F$fleet==byc_f,]
+    new_OM_catch_list$catch<- new_OM_catch_list$catch[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),]
+    new_OM_catch_list$catch_bio<- new_OM_catch_list$catch_bio[which(!is.element(new_OM_catch_list$catch$fleet,byc_f)),]
+    new_OM_catch_list$catch_F<- new_OM_catch_list$catch_F[which(is.element(new_OM_catch_list$catch_F$fleet,byc_f)),]
   } else{
     new_OM_catch_list$catch_F <- NULL
   }
diff --git a/R/parse_MS.R b/R/parse_MS.R
index adbdf08..f6bfb75 100644
--- a/R/parse_MS.R
+++ b/R/parse_MS.R
@@ -66,9 +66,8 @@ parse_MS <- function(MS, EM_out_dir = NULL, EM_init_dir = NULL,
     interim_struct = interim_struct,
     seed = seed
   )
-  
   new_catch_list <- do.call(MS, args = pars_list)
-  
+
   # to do: need better checks on function name? Maybe be more explicit on
   # which environment the function is in?
   # check output before returning
diff --git a/R/runSSMSE.R b/R/runSSMSE.R
index cddaf57..76c51e0 100644
--- a/R/runSSMSE.R
+++ b/R/runSSMSE.R
@@ -96,7 +96,7 @@ run_SSMSE <- function(scen_name_vec,
   if (!is.null(custom_MS_source)) {
     source(custom_MS_source)
   }
-  
+
   # input checks
   if (!all(MS_vec %in% c("EM", "no_catch", "Interim"))) {
     invalid_MS <- MS_vec[unlist(lapply(MS_vec, function(x) !exists(x)))]
@@ -170,7 +170,6 @@ run_SSMSE <- function(scen_name_vec,
     scen_list[[i]][["scen_seed"]] <- scen_seed
   }
   
-  
   if (run_parallel) {
     if (!is.null(n_cores)) {
       n_cores <- min(max(n_cores, 1), (parallel::detectCores() - 1))

From dd93fd7f670a39269f3e10e788543893a401d0f6 Mon Sep 17 00:00:00 2001
From: LaTreese Denson <70522445+latreesedenson-NOAA@users.noreply.github.com>
Date: Mon, 20 May 2024 14:23:33 -0400
Subject: [PATCH 02/10] fixed se name

---
 R/MS_functions.R | 2 +-
 R/RatioBiasEM.R  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/MS_functions.R b/R/MS_functions.R
index 1bb43f0..2c3c3c4 100644
--- a/R/MS_functions.R
+++ b/R/MS_functions.R
@@ -334,7 +334,7 @@ get_EM_catch_df <- function(EM_dir, dat) {
     dis_df <- dis_df %>%
       dplyr::group_by(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]]) %>%
       dplyr::summarise(catch = sum(.data[["Discard"]])) %>%
-      merge(se, all.x = TRUE, all.y = FALSE) %>%
+      merge(se_dis, all.x = TRUE, all.y = FALSE) %>%
       dplyr::ungroup() %>%
       dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
     dis_df <- as.data.frame(dis_df)
diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index 6021e33..9e17506 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -271,7 +271,7 @@ get_RatioEM_catch_df<-function(EM_dir, dat, dat_yrs,
     dis_df <- dis_df %>%
       dplyr::group_by(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]]) %>%
       dplyr::summarise(catch = sum(.data[["Discard"]])) %>%
-      merge(se, all.x = TRUE, all.y = FALSE) %>%
+      merge(se_dis, all.x = TRUE, all.y = FALSE) %>%
       dplyr::ungroup() %>%
       dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
     dis_df <- as.data.frame(dis_df)

From 6f3069d01ff8cd1417ac1ac9b105df2934378f7c Mon Sep 17 00:00:00 2001
From: LaTreese Denson <70522445+latreesedenson-NOAA@users.noreply.github.com>
Date: Tue, 21 May 2024 16:11:40 -0400
Subject: [PATCH 03/10] edit typo in summation

---
 R/MS_functions.R | 2 +-
 R/RatioBiasEM.R  | 2 +-
 R/parse_MS.R     | 1 +
 3 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/R/MS_functions.R b/R/MS_functions.R
index 2c3c3c4..8909468 100644
--- a/R/MS_functions.R
+++ b/R/MS_functions.R
@@ -333,7 +333,7 @@ get_EM_catch_df <- function(EM_dir, dat) {
     dis_df <- do.call("rbind", dis_df_list)
     dis_df <- dis_df %>%
       dplyr::group_by(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]]) %>%
-      dplyr::summarise(catch = sum(.data[["Discard"]])) %>%
+      dplyr::summarise(Discard = sum(.data[["Discard"]])) %>%
       merge(se_dis, all.x = TRUE, all.y = FALSE) %>%
       dplyr::ungroup() %>%
       dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index 9e17506..d108cf6 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -270,7 +270,7 @@ get_RatioEM_catch_df<-function(EM_dir, dat, dat_yrs,
     dis_df <- do.call("rbind", dis_df_list)
     dis_df <- dis_df %>%
       dplyr::group_by(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]]) %>%
-      dplyr::summarise(catch = sum(.data[["Discard"]])) %>%
+      dplyr::summarise(Discard = sum(.data[["Discard"]])) %>%
       merge(se_dis, all.x = TRUE, all.y = FALSE) %>%
       dplyr::ungroup() %>%
       dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
diff --git a/R/parse_MS.R b/R/parse_MS.R
index f6bfb75..aabeabd 100644
--- a/R/parse_MS.R
+++ b/R/parse_MS.R
@@ -66,6 +66,7 @@ parse_MS <- function(MS, EM_out_dir = NULL, EM_init_dir = NULL,
     interim_struct = interim_struct,
     seed = seed
   )
+  
   new_catch_list <- do.call(MS, args = pars_list)
 
   # to do: need better checks on function name? Maybe be more explicit on

From 6ee08202ef900ead3e03fe7d091d9f007bc4a2c9 Mon Sep 17 00:00:00 2001
From: LaTreese Denson <70522445+latreesedenson-NOAA@users.noreply.github.com>
Date: Wed, 22 May 2024 17:50:23 -0400
Subject: [PATCH 04/10] Input bread crumbs for Nathan to follow for fixing the
 code for RS

---
 R/RatioBiasEM.R | 2 +-
 R/parse_MS.R    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index d108cf6..f68b257 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -715,7 +715,7 @@ BiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE,
   # new_OM_catch_list$catch_bio <- NULL
   
   if(!is.null(new_OM_catch_list$discards)){
-    tmp_discards <- base::merge(new_OM_catch_list$discards, sample_struct$EM2OMdiscard_bias, all.x=TRUE)
+    tmp_discards <- base::merge(abs(new_OM_catch_list$discards), abs(sample_struct$EM2OMdiscard_bias), all.x=TRUE) #need to sort to figure this out
     new_OM_catch_list$discards$Discard <- new_OM_catch_list$discards$Discard * tmp_discards$bias 
   }
   
diff --git a/R/parse_MS.R b/R/parse_MS.R
index aabeabd..0c6a007 100644
--- a/R/parse_MS.R
+++ b/R/parse_MS.R
@@ -66,7 +66,7 @@ parse_MS <- function(MS, EM_out_dir = NULL, EM_init_dir = NULL,
     interim_struct = interim_struct,
     seed = seed
   )
-  
+  # browser()
   new_catch_list <- do.call(MS, args = pars_list)
 
   # to do: need better checks on function name? Maybe be more explicit on

From 40af79510c0ace030e33901e014fe7c0366b0c66 Mon Sep 17 00:00:00 2001
From: Nathan Vaughan <nathan.vaughan@noaa.gov>
Date: Thu, 23 May 2024 16:36:21 -0400
Subject: [PATCH 05/10] Red snapper fixes

---
 R/MS_functions.R |  8 ++++++++
 R/RatioBiasEM.R  | 13 +++++++++----
 R/parse_MS.R     |  6 ++----
 R/runSSMSE.R     |  1 +
 4 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/R/MS_functions.R b/R/MS_functions.R
index 8909468..5e569bc 100644
--- a/R/MS_functions.R
+++ b/R/MS_functions.R
@@ -267,6 +267,11 @@ get_EM_catch_df <- function(EM_dir, dat) {
   catch_df <- as.data.frame(catch_df)
   catch_bio_df <- as.data.frame(catch_bio_df)
   catch_F_df <- as.data.frame(catch_F_df)
+  
+  #sort data frames by year, season, fleet
+  catch_df <- catch_df[order(abs(catch_df[["fleet"]]),abs(catch_df[["year"]]),abs(catch_df[["seas"]])),]
+  catch_bio_df <- catch_bio_df[order(abs(catch_bio_df[["fleet"]]),abs(catch_bio_df[["year"]]),abs(catch_bio_df[["seas"]])),]
+  catch_F_df <- catch_F_df[order(abs(catch_F_df[["fleet"]]),abs(catch_F_df[["year"]]),abs(catch_F_df[["seas"]])),]
   # get discard, if necessary
   if (dat[["N_discard_fleets"]] > 0) {
     # discard units: 1, biomass/number according to set in catch
@@ -338,6 +343,9 @@ get_EM_catch_df <- function(EM_dir, dat) {
       dplyr::ungroup() %>%
       dplyr::select(.data[["Yr"]], .data[["Seas"]], .data[["Flt"]], .data[["Discard"]], .data[["Std_in"]])
     dis_df <- as.data.frame(dis_df)
+    #sort data frame by year, season, fleet
+    dis_df <- dis_df[order(abs(dis_df[["Flt"]]),abs(dis_df[["Yr"]]),abs(dis_df[["Seas"]])),]
+    
   } else {
     dis_df <- NULL
   }
diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index f68b257..dc15612 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -703,19 +703,24 @@ BiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE,
   # For the simple approach, we can just apply a series of scalars from the EM catch list to create an OM catch list
   new_OM_catch_list = new_EM_catch_list
   
-  
+  sample_struct$EM2OMcatch_bias <- sample_struct$EM2OMcatch_bias[order(abs(sample_struct$EM2OMcatch_bias$fleet),abs(sample_struct$EM2OMcatch_bias$year),abs(sample_struct$EM2OMcatch_bias$seas)),]
+  sample_struct$EM2OMcatch_bias <- sample_struct$EM2OMcatch_bias[!duplicated(sample_struct$EM2OMcatch_bias),]
   if(!is.null(new_OM_catch_list$catch)){
-    tmp_catch <- base::merge(new_OM_catch_list$catch, sample_struct$EM2OMcatch_bias, all.x=TRUE, all.y=FALSE)
+    tmp_catch <- base::merge(abs(new_OM_catch_list$catch), abs(sample_struct$EM2OMcatch_bias), all.x=TRUE, all.y=FALSE)
+    tmp_catch <- tmp_catch[order(abs(tmp_catch$fleet),abs(tmp_catch$year),abs(tmp_catch$seas)),]
     new_OM_catch_list$catch$catch <- new_OM_catch_list$catch$catch * tmp_catch$bias 
   }
   if(!is.null(new_OM_catch_list$catch_bio)){
-    tmp_catch_bio <- base::merge(new_OM_catch_list$catch_bio, sample_struct$EM2OMcatch_bias, all.x=TRUE)
+    tmp_catch_bio <- base::merge(abs(new_OM_catch_list$catch_bio), abs(sample_struct$EM2OMcatch_bias), all.x=TRUE)
+    tmp_catch_bio <- tmp_catch_bio[order(abs(tmp_catch_bio$fleet),abs(tmp_catch_bio$year),abs(tmp_catch_bio$seas)),]
     new_OM_catch_list$catch_bio$catch <- new_OM_catch_list$catch_bio$catch * tmp_catch_bio$bias
   }
   # new_OM_catch_list$catch_bio <- NULL
-  
+  sample_struct$EM2OMdiscard_bias <- sample_struct$EM2OMdiscard_bias[order(abs(sample_struct$EM2OMdiscard_bias$Flt),abs(sample_struct$EM2OMdiscard_bias$Yr),abs(sample_struct$EM2OMdiscard_bias$Seas)),]
+  sample_struct$EM2OMdiscard_bias <- sample_struct$EM2OMdiscard_bias[!duplicated(sample_struct$EM2OMdiscard_bias),]
   if(!is.null(new_OM_catch_list$discards)){
     tmp_discards <- base::merge(abs(new_OM_catch_list$discards), abs(sample_struct$EM2OMdiscard_bias), all.x=TRUE) #need to sort to figure this out
+    tmp_discards <- tmp_discards[order(abs(tmp_discards$Flt),abs(tmp_discards$Yr),abs(tmp_discards$Seas)),]
     new_OM_catch_list$discards$Discard <- new_OM_catch_list$discards$Discard * tmp_discards$bias 
   }
   
diff --git a/R/parse_MS.R b/R/parse_MS.R
index 0c6a007..30dbdf7 100644
--- a/R/parse_MS.R
+++ b/R/parse_MS.R
@@ -66,16 +66,14 @@ parse_MS <- function(MS, EM_out_dir = NULL, EM_init_dir = NULL,
     interim_struct = interim_struct,
     seed = seed
   )
-  # browser()
+  
   new_catch_list <- do.call(MS, args = pars_list)
 
   # to do: need better checks on function name? Maybe be more explicit on
   # which environment the function is in?
   # check output before returning
   check_catch_df(new_catch_list[["catch"]])
-  if (isTRUE(!is.null(new_catch_list[["discards"]]))) {
-    warning("Discards are not added into the OM for SSMSE currently.")
-  }
+  
   new_catch_list
 }
 
diff --git a/R/runSSMSE.R b/R/runSSMSE.R
index 76c51e0..9f96f77 100644
--- a/R/runSSMSE.R
+++ b/R/runSSMSE.R
@@ -731,6 +731,7 @@ run_SSMSE_iter <- function(out_dir = NULL,
       "Finished running and sampling OM through year ", max(new_catch_list[["catch"]][, "year"]),
       "."
     )
+
     if (run_EM_last_yr == FALSE && isTRUE(yr == test_run_EM_yr)) {
       skip_EM_run <- TRUE
     } else {

From bef3746861f41210c070ee1a66fe07dada420277 Mon Sep 17 00:00:00 2001
From: Nathan Vaughan <nathan.vaughan@noaa.gov>
Date: Mon, 10 Jun 2024 08:20:11 -0400
Subject: [PATCH 06/10] Intitial sampling

Fixed bugs in initial sampling to ensure that bootstrapped values are used rather than exact OM data. Also added catch and discard sampling.
---
 R/initOM.R        |  2 +-
 R/manipulate_EM.R | 24 +++++++++++++++++++++---
 R/runSSMSE.R      |  2 ++
 3 files changed, 24 insertions(+), 4 deletions(-)

diff --git a/R/initOM.R b/R/initOM.R
index 5ec7b05..b0c0f23 100644
--- a/R/initOM.R
+++ b/R/initOM.R
@@ -471,7 +471,7 @@ run_OM <- function(OM_dir,
   if (is.null(seed)) {
     seed <- stats::runif(1, 1, 9999999)
   }
-
+  
   start <- r4ss::SS_readstarter(file.path(OM_dir, "starter.ss"),
     verbose = FALSE
   )
diff --git a/R/manipulate_EM.R b/R/manipulate_EM.R
index c2c035f..8c54147 100644
--- a/R/manipulate_EM.R
+++ b/R/manipulate_EM.R
@@ -66,6 +66,23 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
     check_OM_dat(OM_dat, EM_dat)
   }
   dat <- list(OM_dat = OM_dat, EM_dat = EM_dat)
+  
+  Catches <- lapply(dat, function(x) {
+    tmp <- combine_cols(x, "catch", c("year", "seas", "fleet"))
+  })
+  # match 1 way: match each EM obs with an OM obs. extract only these OM obs.
+  matches <- which(Catches[[1]][, "combo"] %in% Catches[[2]][, "combo"])
+  # extract only the rows of interest and get rid of the "combo" column
+  new_dat[["catch"]] <- Catches[[1]][matches, -ncol(Catches[[1]])]
+  
+  Discards <- lapply(dat, function(x) {
+    tmp <- combine_cols(x, "discard_data", c("Yr", "Seas", "Flt"))
+  })
+  # match 1 way: match each EM obs with an OM obs. extract only these OM obs.
+  matches <- which(Discards[[1]][, "combo"] %in% Discards[[2]][, "combo"])
+  # extract only the rows of interest and get rid of the "combo" column
+  new_dat[["discard_data"]] <- Discards[[1]][matches, -ncol(Discards[[1]])]
+  
   CPUEs <- lapply(dat, function(x) {
     tmp <- combine_cols(x, "CPUE", c("year", "seas", "index"))
   })
@@ -73,6 +90,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   matches <- which(CPUEs[[1]][, "combo"] %in% CPUEs[[2]][, "combo"])
   # extract only the rows of interest and get rid of the "combo" column
   new_dat[["CPUE"]] <- CPUEs[[1]][matches, -ncol(CPUEs[[1]])]
+ 
   # add in lcomps
   if (OM_dat[["use_lencomp"]] == 1) {
     lcomps <- lapply(dat, function(x) {
@@ -85,7 +103,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
     new_dat[["lencomp"]] <- lcomps[[1]][matches_l, -ncol(lcomps[[1]])]
   }
   # add in age comps
-  if (!is.null(dat[["agecomp"]])) {
+  if (!is.null(OM_dat[["agecomp"]])) {
     acomps <- lapply(dat, function(x) {
       tmp <- combine_cols(
         x, "agecomp",
@@ -97,7 +115,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   }
   # TODO: check this for other types of data, esp. mean size at age, k
   # and mean size.
-  if (!is.null(dat[["meanbodywt"]])) {
+  if (!is.null(OM_dat[["meanbodywt"]])) {
     meansize <- lapply(dat, function(x) {
       tmp <- combine_cols(
         x, "meanbodywt",
@@ -107,7 +125,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
     matches_meansize <- which(meansize[[1]][, "combo"] %in% meansize[[2]][, "combo"])
     new_dat[["meanbodywt"]] <- meansize[[1]][matches_meansize, -ncol(meansize[[1]])]
   }
-  if (!is.null(dat[["MeanSize_at_Age_obs"]])) {
+  if (!is.null(OM_dat[["MeanSize_at_Age_obs"]])) {
     size_at_age <- lapply(dat, function(x) {
       tmp <- combine_cols(
         x, "MeanSize_at_Age_obs",
diff --git a/R/runSSMSE.R b/R/runSSMSE.R
index 9f96f77..2156649 100644
--- a/R/runSSMSE.R
+++ b/R/runSSMSE.R
@@ -611,6 +611,7 @@ run_SSMSE_iter <- function(out_dir = NULL,
       verbose = verbose, init_run = TRUE, seed = (iter_seed[["iter"]][1] + 12345)
     )
   }
+  
   if (use_SS_boot == FALSE) {
     stop(
       "Currently, only sampling can be done using the bootstrapping ",
@@ -648,6 +649,7 @@ run_SSMSE_iter <- function(out_dir = NULL,
     seed = (iter_seed[["iter"]][1] + 123456),
     sample_struct = sample_struct # add for bias
   )
+ 
   message(
     "Finished getting catch (years ",
     min(new_catch_list[["catch"]][, "year"]), " to ", max(new_catch_list[["catch"]][, "year"]),

From 479b27c0e201e2e69fd7e0767e064b49684cdb4e Mon Sep 17 00:00:00 2001
From: Nathan Vaughan <nathan.vaughan@noaa.gov>
Date: Mon, 10 Jun 2024 13:56:57 -0400
Subject: [PATCH 07/10] Bug: Missing discards, Mean Size

Fixed bugs in naming of mean size at age columns and a check for discard data before sampling.
---
 R/manipulate_EM.R | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/R/manipulate_EM.R b/R/manipulate_EM.R
index 8c54147..ea38693 100644
--- a/R/manipulate_EM.R
+++ b/R/manipulate_EM.R
@@ -67,6 +67,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   }
   dat <- list(OM_dat = OM_dat, EM_dat = EM_dat)
   
+  if (OM_dat[["discard_data"]] == 1) {
   Catches <- lapply(dat, function(x) {
     tmp <- combine_cols(x, "catch", c("year", "seas", "fleet"))
   })
@@ -74,15 +75,19 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   matches <- which(Catches[[1]][, "combo"] %in% Catches[[2]][, "combo"])
   # extract only the rows of interest and get rid of the "combo" column
   new_dat[["catch"]] <- Catches[[1]][matches, -ncol(Catches[[1]])]
+  }
   
-  Discards <- lapply(dat, function(x) {
-    tmp <- combine_cols(x, "discard_data", c("Yr", "Seas", "Flt"))
-  })
-  # match 1 way: match each EM obs with an OM obs. extract only these OM obs.
-  matches <- which(Discards[[1]][, "combo"] %in% Discards[[2]][, "combo"])
-  # extract only the rows of interest and get rid of the "combo" column
-  new_dat[["discard_data"]] <- Discards[[1]][matches, -ncol(Discards[[1]])]
+  if (OM_dat[["discard_data"]] == 1) {
+    Discards <- lapply(dat, function(x) {
+      tmp <- combine_cols(x, "discard_data", c("Yr", "Seas", "Flt"))
+    })
+    # match 1 way: match each EM obs with an OM obs. extract only these OM obs.
+    matches <- which(Discards[[1]][, "combo"] %in% Discards[[2]][, "combo"])
+    # extract only the rows of interest and get rid of the "combo" column
+    new_dat[["discard_data"]] <- Discards[[1]][matches, -ncol(Discards[[1]])]
+  }
   
+  if (OM_dat[["CPUE"]] == 1) {
   CPUEs <- lapply(dat, function(x) {
     tmp <- combine_cols(x, "CPUE", c("year", "seas", "index"))
   })
@@ -90,7 +95,8 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   matches <- which(CPUEs[[1]][, "combo"] %in% CPUEs[[2]][, "combo"])
   # extract only the rows of interest and get rid of the "combo" column
   new_dat[["CPUE"]] <- CPUEs[[1]][matches, -ncol(CPUEs[[1]])]
- 
+  }
+  
   # add in lcomps
   if (OM_dat[["use_lencomp"]] == 1) {
     lcomps <- lapply(dat, function(x) {
@@ -125,11 +131,12 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
     matches_meansize <- which(meansize[[1]][, "combo"] %in% meansize[[2]][, "combo"])
     new_dat[["meanbodywt"]] <- meansize[[1]][matches_meansize, -ncol(meansize[[1]])]
   }
+  
   if (!is.null(OM_dat[["MeanSize_at_Age_obs"]])) {
     size_at_age <- lapply(dat, function(x) {
       tmp <- combine_cols(
         x, "MeanSize_at_Age_obs",
-        c("Yr", "Seas", "FltSvy", "Gender", "Part", "Ageerr")
+        c("Yr", "Seas", "FltSvy", "Gender", "Part", "AgeErr")
       )
     })
     matches_size_at_age <- which(size_at_age[[1]][, "combo"] %in% size_at_age[[2]][, "combo"])

From 5c846f30d3fb4765ab4b04fe15e24e1d0f7f8335 Mon Sep 17 00:00:00 2001
From: Nathan Vaughan <nathan.vaughan@noaa.gov>
Date: Mon, 10 Jun 2024 14:03:32 -0400
Subject: [PATCH 08/10] bug: fix the last bug fix

Actually fix the input checks :)
---
 R/manipulate_EM.R | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/R/manipulate_EM.R b/R/manipulate_EM.R
index ea38693..efd64cb 100644
--- a/R/manipulate_EM.R
+++ b/R/manipulate_EM.R
@@ -67,7 +67,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   }
   dat <- list(OM_dat = OM_dat, EM_dat = EM_dat)
   
-  if (OM_dat[["discard_data"]] == 1) {
+  if (!is.null(OM_dat[["catch"]])) {
   Catches <- lapply(dat, function(x) {
     tmp <- combine_cols(x, "catch", c("year", "seas", "fleet"))
   })
@@ -77,7 +77,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
   new_dat[["catch"]] <- Catches[[1]][matches, -ncol(Catches[[1]])]
   }
   
-  if (OM_dat[["discard_data"]] == 1) {
+  if (!is.null(OM_dat[["discard_data"]])) {
     Discards <- lapply(dat, function(x) {
       tmp <- combine_cols(x, "discard_data", c("Yr", "Seas", "Flt"))
     })
@@ -87,7 +87,7 @@ get_EM_dat <- function(OM_dat, EM_dat, do_checks = TRUE) {
     new_dat[["discard_data"]] <- Discards[[1]][matches, -ncol(Discards[[1]])]
   }
   
-  if (OM_dat[["CPUE"]] == 1) {
+  if (!is.null(OM_dat[["CPUE"]])) {
   CPUEs <- lapply(dat, function(x) {
     tmp <- combine_cols(x, "CPUE", c("year", "seas", "index"))
   })

From c4e1ce9dbc26bc743470099b237f16628121bb98 Mon Sep 17 00:00:00 2001
From: Nathan Vaughan <nathan.vaughan@noaa.gov>
Date: Fri, 14 Jun 2024 10:41:04 -0400
Subject: [PATCH 09/10] Bug: Discards sort

Correction for discard values being returned in wrong order due to fleet sort missmatch
---
 R/RatioBiasEM.R | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index dc15612..56a324d 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -526,14 +526,17 @@ add_new_dat_BIAS<- function (OM_dat, EM_datfile, sample_struct, EM_dir, nyrs_ass
   SIMPLIFY = FALSE, USE.NAMES = TRUE)
   
   #extracted_dat takes observations from the OM; then below we use the EM2OM multiplier to put catch back into EM units and remove the EM2OMcatch_basis element so that it doesn't get added to the datafile. 
-  
+  extracted_dat$catch <- extracted_dat$catch[order(abs(extracted_dat$catch$fleet),abs(extracted_dat$catch$year),abs(extracted_dat$catch$seas)),]
   if(!is.null(extracted_dat$catch)){
     tmp_catch<- merge(extracted_dat$catch, sample_struct$EM2OMcatch_bias)
+    tmp_catch <- tmp_catch[order(abs(tmp_catch$fleet),abs(tmp_catch$year),abs(tmp_catch$seas)),]
     extracted_dat$catch$catch<- tmp_catch$catch / tmp_catch$bias ### EM2OM edits to get EM catch back to OM
   }
   
+  extracted_dat$discard_data <- extracted_dat$discard_data[order(abs(extracted_dat$discard_data$Flt),abs(extracted_dat$discard_data$Yr),abs(extracted_dat$discard_data$Seas)),]
   if(!is.null(extracted_dat$discard_data)){
     tmp_discard<- merge(extracted_dat$discard_data, sample_struct$EM2OMdiscard_bias)
+    tmp_discard <- tmp_discard[order(abs(tmp_discard$Flt),abs(tmp_discard$Yr),abs(tmp_discard$Seas)),]
     extracted_dat$discard_data$Discard<- tmp_discard$Discard / tmp_discard$bias
   }
   

From 5c549b472f743608a825dd2f62abd9b2d7cc89b8 Mon Sep 17 00:00:00 2001
From: Nathan Vaughan <nathan.vaughan@noaa.gov>
Date: Fri, 14 Jun 2024 11:47:10 -0400
Subject: [PATCH 10/10] Feature: Extend main rec dev phase in em

Update EM functions to increment the end year of the main rec dev phase so that
---
 R/MS_functions.R | 10 ++++++++++
 R/RatioBiasEM.R  | 11 ++++++++++-
 2 files changed, 20 insertions(+), 1 deletion(-)

diff --git a/R/MS_functions.R b/R/MS_functions.R
index 5e569bc..8d1a5a2 100644
--- a/R/MS_functions.R
+++ b/R/MS_functions.R
@@ -91,6 +91,16 @@ EM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE,
       new_datfile_name = new_datfile_name,
       verbose = verbose
     )
+    
+    #Increment main recruitment phase end year by number of assessment years
+    #So the model can continue to estimate rec devs 
+    ctl <- SS_readctl(file.path(EM_out_dir, start[["ctlfile"]]),
+                      datlist = new_EM_dat
+    )
+    ctl$MainRdevYrLast <- ctl$MainRdevYrLast + nyrs_assess
+    r4ss::SS_writectl(ctl, file.path(EM_out_dir, start[["ctlfile"]]),
+                      overwrite = TRUE
+    )
   }
   # Update SS random seed
   start <- SS_readstarter(file.path(EM_out_dir, "starter.ss"),
diff --git a/R/RatioBiasEM.R b/R/RatioBiasEM.R
index 56a324d..b459201 100644
--- a/R/RatioBiasEM.R
+++ b/R/RatioBiasEM.R
@@ -300,7 +300,7 @@ get_RatioEM_catch_df<-function(EM_dir, dat, dat_yrs,
 #'  init_loop is TRUE.
 #' @template OM_out_dir
 #' @template sample_struct
-#' @template seed
+#' @template seed 
 #' @param ... Any additional parameters
 
 #For example, if EM has a positive bias (e.g., EM catch = 1 when true OM catch = 0.5), the EM2OM multiplier should be less than 1 (EM2OM = 0.5)
@@ -662,6 +662,15 @@ BiasEM <- function(EM_out_dir = NULL, init_loop = TRUE, OM_dat, verbose = FALSE,
       verbose = verbose
     )
     
+    #Increment main recruitment phase end year by number of assessment years
+    #So the model can continue to estimate rec devs 
+    ctl <- SS_readctl(file.path(EM_out_dir, start[["ctlfile"]]),
+                      datlist = new_EM_dat
+    )
+    ctl$MainRdevYrLast <- ctl$MainRdevYrLast + nyrs_assess
+    r4ss::SS_writectl(ctl, file.path(EM_out_dir, start[["ctlfile"]]),
+                      overwrite = TRUE
+    )
   } # end else not first iteration
   
   # Update SS random seed