update BAM2OM

Blue-Matter · Apr 18, 2024 · 91397a1 · 91397a1
1 parent 37b4857
commit 91397a1
Show file tree

Hide file tree

Showing 3 changed files with 128 additions and 50 deletions.
diff --git a/NAMESPACE b/NAMESPACE
@@ -22,6 +22,7 @@ export(ASAP2OM)
 export(Assess2MOM)
 export(Assess2OM)
 export(Awatea2OM)
+export(BAM2MOM)
 export(CALsimp)
 export(COSEWIC_Blow)
 export(COSEWIC_Dplot)

diff --git a/R/BAM2OM.R b/R/BAM2OM.R
@@ -1,4 +1,12 @@
 
+convert_units <- function(rdat) {
+  if(rdat$info$units.weight =='lb') {
+    rdat$a.series$weight <- rdat$parms$wgt.a*
+      rdat$a.series$length^rdat$parms$wgt.b
+  }
+  rdat
+}
+
 #' Import a multi-stock, multi-fleet OM from a BAM object
 #'
 #' @param rdat A list object from the `BAMextras` package. Use `bamExtras::standardize_rdat(rdat)`
@@ -14,51 +22,58 @@
 #' @return An object of class `MOM`
 #' @export
 BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
-                    stock_name, fleet_name, LowerTri=0, report = FALSE, ...) {
-
-  a.series <- rdat[["a.series"]] # age data frame
-  maxage <- max(a.series$age)
-  n_age <- maxage + 1
-  sage <- min(a.series$age) # First age class in BAM
-  BAM_age <- a.series$age
-
-  t.series <- rdat[["t.series"]] # time data frame
-  nyears <- nrow(t.series) - 1
-  CurrentYr <- t.series$year[nyears]
+                    stock_name=NULL, fleet_name=NULL, LowerTri=0, report = FALSE, ...) {
 
   parms <- rdat[["parms"]]
 
-  nyears <- length(parms$styr:parms$endyr)
-  CurrentYr <- parms$endyr
+  # Convert all weights to metric (kg)
+  rdat <- convert_units(rdat)
 
   # Grab name of stock
   np <- 1
   snameBAM <- paste(rdat$info$title, rdat$info$species)
-  if(missing(stock_name)) {
-    stock_name <- snameBAM
-  } else if(length(stock_name) != np) {
-    stop("stock_name is not a length ", np, " character vector")
+  if(!is.null(stock_name)) {
+    if(length(stock_name) != np) {
+      stop("stock_name is not a length ", np, " character vector")
+    }
+    snameBAM <- stock_name
   }
   message("BAM model: ", snameBAM)
 
-  # Grab name of fleets
-  fnameBAM <- names(parms)[grepl("F.prop", names(parms))] %>% vapply(function(x) strsplit(x, "F.prop.")[[1]][2], character(1))
+  # Fleet names
+  fnameBAM <- names(parms)[grepl("F.prop", names(parms))] |>
+    vapply(function(x) strsplit(x, "F.prop.")[[1]][2], character(1))
 
   nf <- length(fnameBAM)
-  if(missing(fleet_name)) {
-    fleet_name <- fnameBAM
-  } else if(length(fleet_name) != nf) {
-    stop("fleet_name is not a length ", nf, " character vector")
+  if(!is.null(fleet_name)) {
+    if(length(fleet_name) != nf) {
+      stop("fleet_name is not a length ", nf, " character vector")
+    }
+    fnameBAM <- fleet_name
   }
-  message(nf, "-fleet BAM model found:\n", paste0(fleet_name, collapse = ", "))
+  message(nf, "-fleet BAM model found:\n", paste0(fnameBAM, collapse = ", "))
+
 
+  # Time-series
+  a.series <- rdat[["a.series"]] # age data frame
+  maxage <- max(a.series$age)
+  n_age <- maxage + 1
+  sage <- min(a.series$age) # First age class in BAM
+  BAM_age <- a.series$age
+
+  t.series <- rdat[["t.series"]] # time data frame
+  nyears <- nrow(rdat[["F.age"]])
+  CurrentYr <- t.series$year[nyears]
+
+  ## M-at-Age
   maa <- local({
-    M <- array(MSEtool:::tiny + .Machine$double.eps, c(nsim, n_age, nyears, np))
-    M[, BAM_age + 1, , ] <- array(a.series$M, c(length(BAM_age), nsim, nyears, np)) %>%
+    M <- array(.Machine$double.eps, c(nsim, n_age, nyears, np))
+    M[, BAM_age + 1, , ] <- array(a.series$M, c(length(BAM_age), nsim, nyears, np)) |>
       aperm(c(2, 1, 3, 4))
     M
   })
 
+  ## F-at-Age
   faa_agg <- local({
     FM_agg <- array(0, c(nsim, n_age, nyears))
     FM_agg[, BAM_age + 1, ] <- replicate(nsim, rdat[["F.age"]]) %>% aperm(3:1)
@@ -89,16 +104,18 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
 
   # Check if faa matches faa_agg
   faa_check <- apply(faa[1, , , 1, ], 1:2, sum)
-  message("Difference in aggregate F.age and sum of fishery F-at-age: ", range(faa_agg[1, , ] - faa_check) %>%
-            round(2) %>% paste0(collapse = " to "))
+  message("Difference in aggregate F.age and sum of fishery F-at-age: ", range(faa_agg[1, , ] - faa_check) |>
+            round(2) |> paste0(collapse = " to "))
 
+  ## Weight-at-Age
   waa <- local({
     Wt <- array(0, c(nsim, n_age, nyears, np))
     Wt[, BAM_age + 1, , ] <- array(a.series$weight, c(length(BAM_age), nsim, nyears, np)) %>%
       aperm(c(2, 1, 3, 4))
     Wt
   })
 
+  ## Maturity-at-Age
   mataa <- local({
     if (is.null(a.series$mat.male)) {
       mat_vec <- a.series$mat.female
@@ -112,16 +129,20 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
     mat
   })
 
+  ## Length-at-Age
   laa <- local({
     len <- array(0, c(nsim, n_age, nyears, np))
     len[, BAM_age + 1, , ] <- array(a.series$length, c(length(BAM_age), nsim, nyears, np)) %>%
       aperm(c(2, 1, 3, 4))
     len
   })
 
+  ## N-at-Age
   naa <- local({
-    N <- array(NA_real_, c(nsim, n_age, nyears + 1, np))
-    N[, BAM_age + 1, , ] <- rdat[["N.age"]] %>% array(c(nyears + 1, length(BAM_age), nsim, np)) %>% aperm(c(3:1, 4))
+    add_year <- dim(rdat[["N.age"]])[1] - nyears
+    N <- array(NA_real_, c(nsim, n_age, nyears + add_year, np))
+    N[, BAM_age + 1, , ] <- rdat[["N.age"]] %>% array(c(nyears + add_year, length(BAM_age), nsim, np)) %>% aperm(c(3:1, 4))
+
     if(sage > 0) {
       aind_missing <- sage:1 # Missing cohorts to be filled in
       for(i in 1:length(aind_missing)) {
@@ -132,8 +153,11 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
     N[, , 1:nyears, , drop = FALSE]
   })
 
+
+
+
   if (is.null(a.series$prop.male)) {
-    ##### For red snapper
+    ##### For red snapper & black sea bass
     # a.series$reprod is the product of batch fecundity, n batches per year, proportion female
     fecaa <- local({
       bfec <- array(0, c(nsim, n_age, nyears, np))
@@ -155,9 +179,9 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
   }
 
   # Compare SSB
-  #OM_SSB <- apply(naa * fecaa * exp(-parms[["spawn.time"]] * (faa + maa)), c(1, 3), sum)
-  #plot(OM_SSB[1, ], typ = 'o')
-  #lines(t.series$SSB, pch = 16)
+  # OM_SSB <- apply(naa[,,,1] * fecaa[,,,1] * exp(-parms[["spawn.time"]] * (faa_agg + maa[,,,1])), c(1, 3), sum)
+  # plot(OM_SSB[1, ], typ = 'o', ylim=c(0, max(c(OM_SSB[1,], t.series$SSB), na.rm=TRUE)))
+  # lines(t.series$SSB, pch = 16, col='red')
 
   R0 <- parms$R.virgin.bc #  [["BH.R0"]]
 
@@ -168,7 +192,7 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
   if(is.null(phi0)) phi0 <-  parms[["Phi0"]]
 
   Perr <- parms[["R.sigma.logdevs"]]
-  AC <- acf(t.series$logR.dev, lag.max = 1, plot = FALSE)$acf[2]
+  AC <- acf(t.series$logR.dev, lag.max = 1, plot = FALSE, na.action =na.pass)$acf[2]
 
   message("Running Assess2MOM()...")
 
@@ -195,10 +219,11 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
   )
 
   # Add names
-  names(MOM@Stocks) <- stock_name
-  names(MOM@Fleets) <- stock_name
-  names(MOM@Fleets[[1]]) <- fleet_name
+  names(MOM@Stocks) <- snameBAM
+  names(MOM@Fleets) <- snameBAM
+  names(MOM@Fleets[[1]]) <- fnameBAM
 
+  # custom parameters 
   LatASD <- local({
     lsd <- array(0, c(nsim, n_age, nyears + proyears))
     lsd[, BAM_age + 1, ] <- array(a.series$length.sd, c(length(BAM_age), nsim, nyears + proyears)) %>%
@@ -209,17 +234,18 @@ BAM2MOM <- function(rdat, nsim = 48, proyears = 50, interval = 1,
   for(p in 1:np) {
     MOM@Stocks[[p]]@a <- parms$wgt.a
     MOM@Stocks[[p]]@b <- parms$wgt.b
-    MOM@Stocks[[p]]@Name <- stock_name[p]
+    MOM@Stocks[[p]]@Name <- snameBAM[p]
     for(f in 1:nf) {
-      MOM@Fleets[[p]][[f]]@Name <- fleet_name[f]
+      MOM@Fleets[[p]][[f]]@Name <- fnameBAM[f]
       MOM@cpars[[p]][[f]][["LatASD"]] <- LatASD
       MOM@cpars[[p]][[f]][["Linf"]] <- rep(parms$Linf, nsim)
       MOM@cpars[[p]][[f]][["K"]] <- rep(parms$K, nsim)
       MOM@cpars[[p]][[f]][["t0"]] <- rep(parms$t0, nsim)
     }
     MOM@cpars[[p]][[1]]$spawn_time_frac <- rep(parms$spawn.time, nsim)
   }
-  return(MOM)
+  MOM
+
 }
 
 
@@ -241,20 +267,20 @@ BAM2OM <- function(rdat, nsim = 48, proyears = 50, interval = 2, report = FALSE,
 
   maa <- local({
     M <- array(MSEtool:::tiny + .Machine$double.eps, c(nsim, n_age, nyears))
-    M[, BAM_age + 1, ] <- array(a.series$M, c(length(BAM_age), nsim, nyears)) %>%
+    M[, BAM_age + 1, ] <- array(a.series$M, c(length(BAM_age), nsim, nyears)) |>
       aperm(c(2, 1, 3))
     M
   })
 
   faa <- local({
     FM <- array(0, c(nsim, n_age, nyears))
-    FM[, BAM_age + 1, ] <- replicate(nsim, rdat[["F.age"]]) %>% aperm(3:1)
+    FM[, BAM_age + 1, ] <- replicate(nsim, rdat[["F.age"]]) |> aperm(3:1)
     FM
   })
 
   waa <- local({
     Wt <- array(0, c(nsim, n_age, nyears))
-    Wt[, BAM_age + 1, ] <- array(a.series$weight, c(length(BAM_age), nsim, nyears)) %>%
+    Wt[, BAM_age + 1, ] <- array(a.series$weight, c(length(BAM_age), nsim, nyears)) |>
       aperm(c(2, 1, 3))
     Wt
   })
@@ -267,21 +293,21 @@ BAM2OM <- function(rdat, nsim = 48, proyears = 50, interval = 2, report = FALSE,
       mat_vec <- a.series$mat.female * a.series$prop.female + a.series$mat.male * (1 - a.series$prop.female)
     }
     mat <- array(0, c(nsim, n_age, nyears))
-    mat[, BAM_age + 1, ] <- array(mat_vec, c(length(BAM_age), nsim, nyears)) %>%
+    mat[, BAM_age + 1, ] <- array(mat_vec, c(length(BAM_age), nsim, nyears)) |>
       aperm(c(2, 1, 3))
     mat
   })
 
   laa <- local({
     len <- array(0, c(nsim, n_age, nyears))
-    len[, BAM_age + 1, ] <- array(a.series$length, c(length(BAM_age), nsim, nyears)) %>%
+    len[, BAM_age + 1, ] <- array(a.series$length, c(length(BAM_age), nsim, nyears)) |>
       aperm(c(2, 1, 3))
     len
   })
 
   naa <- local({
     N <- array(NA_real_, c(nsim, n_age, nyears + 1))
-    N[, BAM_age + 1, ] <- replicate(nsim, rdat[["N.age"]]) %>% aperm(3:1)
+    N[, BAM_age + 1, ] <- replicate(nsim, rdat[["N.age"]]) |> aperm(3:1)
     if(sage > 0) {
       aind_missing <- sage:1 # Missing cohorts to be filled in
       for(i in 1:length(aind_missing)) {
@@ -298,7 +324,7 @@ BAM2OM <- function(rdat, nsim = 48, proyears = 50, interval = 2, report = FALSE,
     fecaa <- local({
       bfec <- array(0, c(nsim, n_age, nyears))
 
-      bfec[, BAM_age + 1, ] <- array(a.series$reprod, c(length(BAM_age), nsim, nyears)) %>%
+      bfec[, BAM_age + 1, ] <- array(a.series$reprod, c(length(BAM_age), nsim, nyears)) |>
         aperm(c(2, 1, 3))
       bfec * mataa
     })
@@ -308,7 +334,7 @@ BAM2OM <- function(rdat, nsim = 48, proyears = 50, interval = 2, report = FALSE,
     fecaa <- local({
       fec <- array(0, c(nsim, n_age, nyears))
 
-      fec[, BAM_age + 1, ] <- array(a.series$reprod, c(length(BAM_age), nsim, nyears)) %>%
+      fec[, BAM_age + 1, ] <- array(a.series$reprod, c(length(BAM_age), nsim, nyears)) |>
         aperm(c(2, 1, 3))
       fec
     })
@@ -353,7 +379,7 @@ BAM2OM <- function(rdat, nsim = 48, proyears = 50, interval = 2, report = FALSE,
 
   OM@cpars$LatASD <- local({
     lsd <- array(0, c(nsim, n_age, nyears + proyears))
-    lsd[, BAM_age + 1, ] <- array(a.series$length.sd, c(length(BAM_age), nsim, nyears + proyears)) %>%
+    lsd[, BAM_age + 1, ] <- array(a.series$length.sd, c(length(BAM_age), nsim, nyears + proyears)) |>
       aperm(c(2, 1, 3))
     lsd
   })

diff --git a/man/BAM2MOM.Rd b/man/BAM2MOM.Rd