R/SS2OM.R

In Blue-Matter/MSEtool: Management Strategy Evaluation Toolkit

#' @rdname SS2MOM
#' @param import_mov Logical. Import movement matrix?
#' @param seasons_to_years Logical, when season is the time step, whether to convert OM from a seasonal model to annual model.
#' @param model_discards Logical, how to simplify a multi-fleet SS model to an OM object. If TRUE, OM will still model discards using
#' the mean retention across fleets (weighted by fleet F). Otherwise, no discards are modeled and all fishing removals are calculated in the OM
#' from the SS F-at-age matrix.
#' @param Author Who did the assessment
#' @param report Logical, if TRUE, the function will run \link{runMSE} to generate the Hist object from the operating model
#' to compare against SS output. A markdown report will be generated.
#' @param filename If \code{report = TRUE}, character string for the name of the markdown and HTML files.
#' @param dir If \code{report = TRUE}, the directory in which the markdown and HTML files will be saved.
#' @param open_file If \code{report = TRUE}, whether the HTML document is opened after it is rendered.
#' @param gender An integer that indexes the sex for importing life history parameters (1 = usually female, 2 = usually male, 1:2 = mean across both sexes).
#' Only used for \code{SS2OM} only in a 2-sex model.
#' @export
SS2OM <- function(SSdir, nsim = 48, proyears = 50, reps = 1, maxF = 3, seed = 1, interval = 1, pstar = 0.5,
                  Obs = MSEtool::Generic_Obs, Imp = MSEtool::Perfect_Imp,
                  import_mov = TRUE, gender = 1:2, seasons_to_years = TRUE,
                  model_discards = TRUE, silent = FALSE,
                  Name = "OM generated by SS2OM function",
                  Source = "No source provided", Author = "No author provided",
                  report = FALSE, filename = "SS2OM", dir = tempdir(), open_file = TRUE, ...) {

  if(is.list(SSdir)) {
    replist <- SSdir
  } else {
    replist <- SS_import(SSdir, silent, ...)
  }

  if(replist$nsexes == 1) gender <- 1

  if(!silent) message("Converting SS output to MOM...")
  MOM <- SS2MOM(replist, nsim = nsim, proyears = proyears, reps = reps, maxF = maxF, seed = seed,
                interval = interval, pstar = pstar, Obs = Obs, Imp = Imp, silent = silent,
                Name = Name, Source = Source)

  if(!silent) message("Converting MOM to OM...")
  OM <- SSMOM2OM(MOM, replist, gender, import_mov, seed, silent, model_discards)
  
  if(replist$nseasons == 1 && replist$seasduration < 1 && seasons_to_years) {
    message("Model with season as years found. Will convert to annual time step.")
    OM <- SS_seasonalyears_to_annual(OM, replist)
  }
  if(report) plot_SS2OM(OM, replist, gender, filename, dir, open_file, silent)

  return(OM)
}


mean_vector <- function(tt, cpars, gender) {
  out <- lapply(cpars[gender], function(x) x[[1]][[tt]])
  if(length(out) == 2 && is.null(out[[2]])) {
    res <- out[[1]]
  } else {
    res <- simplify2array(out) %>% apply(1, mean,na.rm=TRUE)
  }
  return(res)
}

# This function grabs array tt from the cpars of the first fleet of each stock, and averages across genders
mean_array <- function(tt, cpars, gender) {
  xx <- lapply(cpars[gender], function(x) x[[1]][[tt]])
  Reduce("+", xx) / length(xx)
}

mean_array2 <- function(list) {
  Reduce("+", list) / length(list)
}

#' @rdname SS2MOM
#' @param MOM MOM object
#' @export
SSMOM2OM <- function(MOM, SSdir, gender = 1:2, import_mov = TRUE, seed = 1, silent = FALSE, model_discards = TRUE) {
  Factor <- NULL # variable for binding check

  if(!requireNamespace("reshape2", quietly = TRUE)) {
    stop("Package `reshape2` is required for this function. Install with `install.packages('reshape2')`", call. = FALSE)
  }

  if(is.list(SSdir)) {
    replist <- SSdir
  } else {
    replist <- SS_import(SSdir, silent)
  }

  if(replist$nsexes == 1) gender <- 1

  Stocks <- MOM@Stocks
  Fleets <- MOM@Fleets
  Obs <- MOM@Obs[[1]][[1]]
  Imp <- MOM@Imps[[1]][[1]]
  cpars <- MOM@cpars
  nsim <- MOM@nsim
  nyears <- Fleets[[1]][[1]]@nyears
  proyears <- MOM@proyears
  mainyrs <- replist$startyr:replist$endyr

  Stock <- new("Stock")
  Stock@maxage <- vapply(Stocks, slot, numeric(1), "maxage") %>% unique()
  Stock@R0 <- vapply(Stocks, slot, numeric(1), "R0") %>% sum()

  Stock@SRrel <- Stocks[[1]]@SRrel
  Stock@h <- Stocks[[1]]@h
  Stock@D <- rep(replist$current_depletion, 2)
  Stock@AC <- Stocks[[1]]@AC
  Stock@Perr <- Stocks[[1]]@Perr

  Stock@Msd <- Stock@Linfsd <- Stock@Ksd <- Stock@L50 <- Stock@L50_95 <- c(0, 0)
  Stock@Size_area_1 <- Stock@Frac_area_1 <- Stock@Prob_staying <- c(0.5, 0.5)

  # cpars
  cpars_out <- list()

  cpars_out$M_ageArray <- mean_array("M_ageArray", cpars, gender)
  cpars_out$Wt_age <- mean_array("Wt_age", cpars, gender)
  cpars_out$Len_age <- mean_array("Len_age", cpars, gender)
  cpars_out$LatASD <- mean_array("LatASD", cpars, gender)
  cpars_out$Fec_age <- mean_array("Fec_age", cpars, gender=1)
  if(!length(cpars_out$Fec_age)) cpars_out$Fec_age <- NULL

  cpars_out$Linf <- mean_vector("Linf", cpars, gender)
  cpars_out$K <- mean_vector("K", cpars, gender)
  cpars_out$t0 <- mean_vector("t0", cpars, gender)

  # empirical weight-at-age for catch - weighted average over fleets 
  # and averaged over stocks
  wt_list <- list()
  for (st in 1:length(MOM@Stocks)) {
    wt_list[[st]] <- calc_weightedmean_c(MOM@cpars[[st]])
  }
  cpars_out$Wt_age_C <- mean_array2(wt_list) # empirical weight-at-age for catch
  if(length(cpars_out$Wt_age_C)==0) cpars_out$Wt_age_C <- NULL

  # Stock placeholders (overriden by cpars mean_arrays or mean_vectors above)
  Stock@M <- vapply(Stocks[gender], slot, numeric(2), "M") %>% apply(1, mean)
  Stock@LenCV <- vapply(Stocks[gender], slot, numeric(2), "LenCV") %>% apply(1, mean)

  Stock@Linf <- cpars_out$Linf %>% range()
  Stock@K <- cpars_out$K %>% range()
  Stock@t0 <- cpars_out$t0 %>% range()

  slot2 <- function(x, y) ifelse(length(slot(x, y)), slot(x, y), NA_real_)
  Stock@a <- vapply(Stocks[gender], slot2, numeric(1), "a") %>% mean(na.rm = TRUE)
  Stock@b <- vapply(Stocks[gender], slot2, numeric(1), "b") %>% mean(na.rm = TRUE)

  # cpars for the first gender, first fleet
  .cpars <- cpars[[1]][[1]]
  cpars_out$hs <- .cpars$hs
  #cpars_out$binWidth <- .cpars$binWidth
  cpars_out$CAL_bins <- .cpars$CAL_bins
  cpars_out$CAL_binsmid <- .cpars$CAL_binsmid
  cpars_out$Mat_age <- .cpars$Mat_age

  cpars_out$Perr_y <- .cpars$Perr_y

  # Do movement
  if(import_mov && !is.null(replist$movement) && nrow(replist$movement)) {
    movement <- replist$movement[replist$movement$Seas == 1 & replist$movement$Gpattern == 1, ]
    if(!nrow(movement)) movement <- replist$movement[replist$movement$Seas == 1 & replist$movement$GP == 1, ]

    nareas <- length(unique(movement$Source_area))
    if(!silent) message(nareas, " area model found. Parameterizing movement matrix.")

    full_movement <- vapply(0:Stock@maxage, function(x) movement[, paste0("age", x)], numeric(nrow(movement)))

    cpars_out$mov <- array(0, c(nsim, Stock@maxage + 1, nareas, nareas))

    for(i in 1:nrow(full_movement)) {
      from <- movement$Source_area[i]
      to <- movement$Dest_area[i]
      for(j in 1:ncol(full_movement)) cpars_out$mov[1:nsim, j, from, to] <- full_movement[i, j]
    }
    #mov[is.na(mov)] <- 0
    #cpars_out$mov <- mov
  }

  # Fleet
  Fleet <- new("Fleet")
  Fleet@Name <- names(Fleets[[1]]) %>% paste(collapse = ", ")
  Fleet@nyears <- nyears
  Fleet@CurrentYr <- max(mainyrs)

  # Placeholders
  Fleet@Spat_targ <- c(1, 1)
  Fleet@Esd <- Fleet@qinc <- Fleet@qcv <- Fleet@L5 <- Fleet@LFS <- Fleet@Vmaxlen <- c(0, 0)
  Fleet@isRel <- FALSE
  Fleet@DR <- Stock@Fdisc <- c(0, 0) # No discards
  Fleet@MPA <- FALSE

  if(model_discards) { # Grab sel and retention (averaged over all fleets) and take the mean between sexes
    sel_par <- lapply(cpars[gender], calculate_single_fleet_dynamics)
    
    cpars_out$Find <- lapply(sel_par, getElement, "Find") %>% mean_array2()
    cpars_out$Fdisc <- lapply(sel_par, getElement, "Fdisc") %>% mean_array2()
    cpars_out$Fdisc_array1 <- lapply(sel_par, getElement, "Fdisc_array1") %>% mean_array2()
    cpars_out$Fdisc_array2 <- lapply(sel_par, getElement, "Fdisc_array2") %>% mean_array2()
  
    cpars_out$retA <- lapply(sel_par, getElement, "retA_real") %>% mean_array2()
    cpars_out$V <- lapply(sel_par, getElement, "V_real") %>% mean_array2()
    cpars_out$retL <- lapply(sel_par, getElement, "retL_real") %>% mean_array2()
    cpars_out$SLarray <- lapply(sel_par, getElement, "SLarray_real") %>% mean_array2()
    
  } else { # Calc V from the F-at-age matrix (no discards modeled in the OM object)

    FF <- dplyr::filter(replist$ageselex, Factor == "F", Yr %in% mainyrs)
    if(nrow(FF)) { # Search for F-at-age in replist$ageselex otherwise go to replist$exploitation
      # Need to average over seasons
      # Selectivity = all dead catch
      # No discards are modeled
      FF <- dplyr::filter(replist$ageselex, Factor == "F", Yr %in% mainyrs)[, -c(1, 4, 6, 7)] %>%
        reshape2::melt(list("Fleet", "Yr", "Sex"), variable.name = "Age", value.name = "F") %>%
        group_by(Yr, Sex, Age) %>% summarise(F = sum(F)) %>% # Sum over fleets
        group_by(Yr, Age) %>% summarise(F = mean(F)) %>% # Mean over sexes
        group_by(Yr) %>% mutate(F_denom = pmax(F, 1e-8), V = F/max(F_denom)) # Get vulnerability, F_denom avoids divide-by-zero
      Find <- group_by(FF, Yr) %>% summarise(F = max(F)) %>% getElement("F")
      V_temp <- reshape2::acast(FF, Age ~ Yr, value.var = "V")
    } else {
      # Selectivity = All mortality
      Vfleet <- lapply(c(1:replist$nfleets)[replist$IsFishFleet], get_V_from_Asel2, i = gender,
                       replist = replist, mainyrs = mainyrs, maxage = Stock@maxage)

      Fapic <- replist$exploitation[, replist$FleetNames[replist$IsFishFleet] %in% colnames(replist$exploitation)][, -3] %>%
        dplyr::filter(Yr %in% mainyrs) %>% structure(names = c("Yr", "Seas", c(1:replist$nfleets)[replist$IsFishFleet])) %>%
        reshape2::melt(value.name = "FF", variable.name = "Fleet", id.vars = c("Yr", "Seas")) %>% group_by(Yr, Fleet) %>%
        summarise(FF = sum(FF)) %>% reshape2::acast(list("Yr", "Fleet"), value.var = "FF")

      F_at_age <- lapply(1:length(Vfleet), function(x) t(Vfleet[[x]]) * Fapic[, x]) %>%
        simplify2array() %>% apply(1:2, sum)

      Find <- apply(F_at_age, 1, max)
      F_denom <- pmax(Find, 1e-8)
      V_temp <- t(F_at_age/F_denom)
    }
    V <- vapply(1:ncol(V_temp), function(x) { # Grab selectivity from neighboring years if all zeros
      x <- min(x, ncol(V_temp))
      if(!sum(V_temp[, x], na.rm = TRUE) && x < ncol(V_temp)) {
        Recall(x+1)
      } else {
        return(V_temp[, x])
      }
    }, numeric(nrow(V_temp)))

    Vpro <- V[, ncol(V)] %>% matrix(nrow(V), proyears)
    cpars_out$V <- cbind(V, Vpro) %>% array(c(Stock@maxage + 1, nyears + proyears, nsim)) %>% aperm(c(3, 1 ,2))
    cpars_out$Find <- matrix(Find, nsim, nyears, byrow = TRUE)
  }
  
  cpars_out$qs <- rep(1, nsim)

  # Place holders
  Fleet@EffYears <- 1:nyears
  Fleet@EffLower <- Fleet@EffUpper <- apply(cpars_out$Find, 2, mean)

  OM <- suppressMessages(new("OM", Stock = Stock, Fleet = Fleet, Obs = Obs, Imp = Imp))
  OM@Name <- MOM@Name
  OM@nsim <- MOM@nsim
  OM@proyears <- MOM@proyears
  OM@reps <- MOM@reps
  OM@maxF <- MOM@maxF
  OM@seed <- MOM@seed
  OM@interval <- MOM@interval
  OM@pstar <- MOM@pstar
  OM@cpars <- cpars_out
  return(OM)
}

#' @rdname SS2MOM
#' @param x For \code{plot_SS2OM}, an object of either class \linkS4class{OM} or \linkS4class{Hist}.
#'  For \code{plot_SS2MOM}, an object of either class \linkS4class{MOM} or `multiHist`.
#' @export
plot_SS2OM <- function(x, SSdir, gender = 1:2,
                       filename = "SS2OM", dir = tempdir(), open_file = TRUE, silent = FALSE, ...) {
  if(missing(SSdir)) stop("SSdir not found.")

  if(inherits(x, "OM")) {
    if(!silent) message("Generating Hist object from OM...")
    OM <- x
    Hist <- runMSE(OM, Hist = TRUE, silent = silent)
  } else if(inherits(x, "Hist")) {
    Hist <- x
  } else {
    stop("Neither Hist nor OM object was found.", call. = FALSE)
  }

  if(is.list(SSdir)) {
    replist <- SSdir
  } else {
    replist <- SS_import(SSdir, silent, ...)
  }

  if(replist$nsexes == 1) gender <- 1

  rmd_file <- file.path(system.file(package = "MSEtool"), "Rmd", "SS", "SS2OM.Rmd")
  rmd <- readLines(rmd_file)

  write(rmd, file = file.path(dir, paste0(filename, ".rmd")))

  if(!silent) message("Rendering markdown file to HTML: ", file.path(dir, paste0(filename, ".html")))

  out <- rmarkdown::render(file.path(dir, paste0(filename, ".rmd")), "html_document", paste0(filename, ".html"), dir,
                           output_options = list(df_print = "paged"), quiet = TRUE)
  message("Rendering complete.")

  if(open_file) browseURL(out)
  return(invisible(out))
}


calculate_single_fleet_dynamics <- function(x) {
  # x is a list of cpars by fleet for a single sex
  # If nfleet > 1, selectivity and retention calculated from weighted sum of fleet F. Otherwise just grab from the first fleet
  # Use mean discard mortality

  if(length(x) > 1) {
    SLarray <- sapply(x, getElement, "SLarray", simplify = "array")
    retL <- sapply(x, getElement, "retL", simplify = "array")
    V <- sapply(x, getElement, "V", simplify = "array")
    retA <- sapply(x, getElement, "retA", simplify = "array")
    
    Find <- sapply(x, function(xx) xx$qs * xx$Find, simplify = "array")
    Fdisc1 <- sapply(x, getElement, "Fdisc_array1", simplify = "array")
    Fdisc2 <- sapply(x, getElement, "Fdisc_array2", simplify = "array")

    nsim <- dim(Find)[1]
    nyears <- dim(Find)[2]
    proyears <- dim(V)[3] - nyears

    Fdisc_avg <- sapply(x, getElement, "Fdisc", simplify = "array") %>% 
      apply(1, mean, na.rm = TRUE)
    
    # Age-based arrays
    F_at_age <- single_fleet_F_at_age(Find, V, retA, Fdisc1, nsim, nyears)
    Fdisc1_avg <- single_fleet_Fdisc(Find, V, retA, Fdisc1, nsim, nyears)
    retA_avg <- single_fleet_retA(Find, V, retA, Fdisc1_avg, F_at_age, nsim, nyears) 
    
    # Find_out <- apply(F_at_age[, -1, ]/retA_avg[, -1, ], c(1, 3), max) # Ignore age-0 to avoid likely outliers for apical Find 
    Find_out <- apply(F_at_age, c(1, 3), max) # F_at_age already includes discards
    V_avg <- single_fleet_V(F_at_age, retA = retA_avg, Find = Find_out, nsim = nsim)
    
    retA_avg <- retA_avg*V_avg
    
    # maxA <- apply(retA_avg, c(1,3), max)
    # nage <- dim(retA_avg)[2]
    # maxA <- aperm(replicate(nage, maxA), c(1,3,2))
    # retA_avg <- retA_avg/maxA 

    # Length-based arrays
    F_at_length <- single_fleet_F_at_age(Find, V = SLarray, retA = retL, 
                                         Fdisc1 = Fdisc2, nsim = nsim, nyears = nyears)
    Fdisc2_avg <- single_fleet_Fdisc(Find = Find, V = SLarray, retA = retL, 
                                     Fdisc1 = Fdisc2, nsim = nsim, nyears = nyears) 
    retL_avg <- single_fleet_retA(Find = Find, V = SLarray, retA = retL, 
                                  Fdisc1_avg = Fdisc2_avg, 
                                  F_at_age = F_at_length, nsim = nsim, 
                                  nyears = nyears)
    
    # Find_out_L <- apply(F_at_length/retL_avg, c(1, 3), max) # Should be almost identical to Find_out
    Find_out_L <- apply(F_at_length, c(1, 3), max) 
    SL_avg <- single_fleet_V(F_at_age = F_at_length, retA = retL_avg, 
                             Find = Find_out_L, nsim = nsim)
    out <- list(Find = Find_out, 
                V_real = abind::abind(V_avg, replicate(proyears, V_avg[, , nyears]), along = 3),
                retA_real = abind::abind(retA_avg, replicate(proyears, retA_avg[, , nyears]), along = 3),
                Fdisc = Fdisc_avg,
                SLarray_real = abind::abind(SL_avg, replicate(proyears, SL_avg[, , nyears]), along = 3),
                retL_real = abind::abind(retL_avg, replicate(proyears, retL_avg[, , nyears]), along = 3),
                Fdisc_array1 = abind::abind(Fdisc1_avg, replicate(proyears, Fdisc1_avg[, , nyears]), along = 3),
                Fdisc_array2 = abind::abind(Fdisc2_avg, replicate(proyears, Fdisc2_avg[, , nyears]), along = 3))

  } else {
    fvar <- c("Find", "V_real", "retA_real", "Fdisc", "SLarray_real", "retL_real", "Fdisc_array1", "Fdisc_array2")
    out <- lapply(fvar, function(xx) getElement(x[[1]], xx)) %>% structure(names = fvar)
  }
  return(out)
}

single_fleet_F_at_age <- function(Find, V, retA, Fdisc1, nsim, nyears) {
  sapply(1:nsim, function(i) { # Sum across fleets in sim i, year j
    sapply(1:nyears, function(j) { 
      Fretain <- Find[i, j, ] * t(V[i, , j, ] * retA[i, , j, ])
      Fdiscard <- Find[i, j, ] * t(V[i, , j, ] * (1 - retA[i, , j, ]) * Fdisc1[i, , j, ])
      colSums(Fretain + Fdiscard)
    })
  }, simplify = "array") %>% aperm(c(3, 1, 2))
}

single_fleet_Fret_at_age <- function(Find, V, retA, Fdisc1, nsim, nyears) {
  sapply(1:nsim, function(i) { # Sum across fleets in sim i, year j
    sapply(1:nyears, function(j) { 
      Fretain <- Find[i, j, ] * t(V[i, , j, ] * retA[i, , j, ])
      Fdiscard <- Find[i, j, ] * t(V[i, , j, ] * (1 - retA[i, , j, ]) * Fdisc1[i, , j, ])
      colSums(Fretain)
    })
  }, simplify = "array") %>% aperm(c(3, 1, 2))
}

single_fleet_Fdisc <- function(Find, V, retA, Fdisc1, nsim, nyears) {
  sapply(1:nsim, function(i) { # Weighted mean across fleets in sim i, year j
    res <- sapply(1:nyears, function(j) {
      num <- Find[i, j, ] * t(V[i, , j, ] * (1 - retA[i, , j, ]) * Fdisc1[i, , j, ])
      den <- Find[i, j, ] * t(V[i, , j, ] * (1 - retA[i, , j, ]))
      out <- colSums(num)/colSums(den)
      out[is.na(out)] <- mean(out, na.rm = TRUE) # Denominator is zero because retention = 1
      return(out)
    })
    if(any(is.na(res))) { # There are years with entirely NA vectors, return mean of time series
      m <- apply(res, 1, mean, na.rm = TRUE)
      res <- apply(res, 2, function(jj) {
        if(all(is.na(jj))) {
          m
        } else {
          jj
        }
      })
    }
    return(res)
  }, simplify = "array") %>% aperm(c(3, 1, 2))
}




single_fleet_retA <- function(Find, V, retA, Fdisc1_avg, F_at_age, nsim, nyears) {
  
  Fretain <- sapply(1:nsim, function(i) {
  sapply(1:nyears, function(j) colSums(Find[i, j, ] * t(V[i, , j, ] * retA[i, , j, ])))
  }, simplify = "array") %>% aperm(c(3, 1, 2))

  Fdiscard <- F_at_age - Fretain
  
  sapply(1:nsim, function(i) {
    res <- sapply(1:nyears, function(j) {
      out <- Fretain[i, , j] * Fdisc1_avg[i, , j]/(Fdiscard[i, , j] + Fretain[i, , j] * Fdisc1_avg[i, , j])
      out[is.na(out)] <- mean(out, na.rm = TRUE) # Denominator is zero because F = 0, should ensure that: all(out > 0)
      return(out)
    })
    if(any(is.na(res))) { # There are years with entirely NA vectors, return mean of time series
      m <- apply(res, 1, mean, na.rm = TRUE)
      res <- apply(res, 2, function(jj) {
        if(all(is.na(jj))) {
          return(m)
        } else {
          return(jj)
        }
      })
    }
    return(res)
  }, simplify = "array") %>% aperm(c(3, 1, 2))
  
}

single_fleet_V <- function(F_at_age, retA, Find, nsim) {
  sapply(1:nsim, function(i) {
    # res <- t(F_at_age[i, , ]/retA[i, , ])/Find[i, ]
    res <- t(F_at_age[i, , ])/Find[i, ]
    if(any(is.na(res))) { # There are years with entirely NA vectors (Find = 0), return mean of time series
      m <- apply(res, 2, mean, na.rm = TRUE)
      res <- apply(res, 1, function(jj) {
        if(all(is.na(jj))) {
          return(m)
        } else {
          return(jj)
        }
      }) %>% t()
    }
    # Make robust by fixing ratio of V in age-0/age-1 to that in F_at_age (see BET example)
    Fratio <- F_at_age[i, 1, ] / F_at_age[i, 2, ]
    Fratio[is.na(Fratio)] <- mean(Fratio, na.rm = TRUE)
    res[, 1] <- res[, 2] * Fratio
    return(res)
  }, simplify = "array") %>% aperm(3:1)
}

calc_weightedmean_c <- function(l) {
  W <- fm <- relF <- NULL
  Wt_age_C <- lapply(l, function(x) x$Wt_age_C)
  
  Find <- lapply(l, function(x) x$Find)
  V <- lapply(l, function(x) x$V)
  nyears <- ncol(Find[[1]])
  totyears <- dim(Wt_age_C[[1]])[3]
  if (is.null(totyears)) 
    return(NULL)
  nage <- dim(Wt_age_C[[1]])[2]
  proyears <- totyears-nyears
  nsim <- nrow(Find[[1]])
  nfleets <- length(Wt_age_C)
  
  wtlist <- list()
  for (fl in 1:nfleets) {
    
    Fdf <- data.frame(Sim=rep(1:nsim), 
                      Yr=rep(1:nyears,each=nsim),
                      fm=as.vector(Find[[fl]]))
    
    Vdf <- data.frame(Sim=rep(1:nsim), 
                      Age=rep(0:(nage-1), each=nsim),
                      Yr=rep(1:nyears,each=nage*nsim),
                      V=as.vector(V[[fl]][,,1:nyears]))
                       
    
    df <- left_join(Fdf, Vdf, by = c("Sim", "Yr"))
    df <- df %>% mutate(relF=fm*V)
    
    
  
    wtdf <- data.frame(Sim=rep(1:nsim),
                       Yr=rep(1:nyears,each=nage*nsim),
                       W=as.vector(Wt_age_C[[fl]][,,1:nyears]),
                       Fl=fl,
                       Age=rep(0:(nage-1), each=nsim))
    df <- left_join(df, wtdf, by = c("Sim", "Yr", "Age"))
    df$relF[is.na(df$relF)] <- 1
    wtlist[[fl]] <- df
  }
  df <- do.call('rbind', wtlist)
  
  wdf_hist<- df %>% group_by(Sim, Yr, Age) %>% 
    summarize(W=weighted.mean(x=W, w=relF), .groups='drop')
  
  
  # add projection years 
  lastYr <- wdf_hist %>% dplyr::ungroup() %>% dplyr::filter(Yr==max(Yr))
  proj_list <- list()
  for (y in 1:proyears) {
    temp <- lastYr
    temp$Yr <- nyears+y
    proj_list[[y]] <- temp
  }
  proj_array <- do.call('rbind', proj_list)
  wdf <- bind_rows(wdf_hist, proj_array) %>% 
    arrange(Sim, Yr, Age)
  
  wt <- array(NA, dim=c(nage, totyears, nsim))
  wt[] <- wdf$W

  out <- aperm(wt, c(3,1,2))
  out
  
}

#' @describeIn SS2MOM Aggregate all fleets in an MOM object.
#' @export
MOM_agg_fleets <- function(MOM) {
  nsim <- MOM@nsim
  MOM_aggfleet <- MOM
  n.stock <- length(MOM@Stocks)
  sel_par <- lapply(MOM@cpars, calculate_single_fleet_dynamics)
  
  fleet <- MOM@Fleets[[1]][[1]]
  MOM_aggfleet@Fleets[[1]] <- list()
  MOM_aggfleet@Fleets[[2]] <- list()
  MOM_aggfleet@Fleets[[1]][[1]] <- fleet
  MOM_aggfleet@Fleets[[1]][[1]]@Name <- 'Aggregated Female'
  MOM_aggfleet@Fleets[[2]][[1]] <- fleet
  MOM_aggfleet@Fleets[[2]][[1]]@Name <- 'Aggregated Male'
  MOM_aggfleet@cpars[[1]] <- list()
  MOM_aggfleet@cpars[[2]] <- list()
  
  MOM_aggfleet@cpars[[1]][[1]] <- MOM@cpars[[1]][[1]]
  MOM_aggfleet@cpars[[2]][[1]] <- MOM@cpars[[2]][[1]]
  
  MOM_aggfleet@CatchFrac <- vector('list', 2)
  # update cpars 
  for (st in 1:2) {
    # update aggegrated selectivity
    for (nm in names(sel_par[[st]])) {
      MOM_aggfleet@cpars[[st]][[1]][[nm]] <- sel_par[[st]][[nm]]
    }
    # weighted average empirical catch-weight over fleets
    MOM_aggfleet@cpars[[st]][[1]]$Wt_age_C  <- calc_weightedmean_c(MOM@cpars[[st]])
    MOM_aggfleet@CatchFrac[[st]] <- matrix(apply(MOM@CatchFrac[[st]], 1, sum), nrow=nsim, ncol=1)
    
  }
  
  MOM_aggfleet
}