R/SS2MOM.R

In MSEtool: Management Strategy Evaluation Toolkit

#' Import Stock Synthesis to MOM (2-sex multi-fleet) or OM (single-sex, single-fleet)
#'
#' @description Functions that uses the file location or the r4ss output list of a fitted SS3 model including input files to populate the
#' various slots of an \linkS4class{MOM} or \linkS4class{OM} object. \code{SS2MOM} and \code{SS2OM} mainly populates the Stock and Fleet components
#' components of the operating model. \code{SS2MOM} creates a 2-sex model and multiple fleets with discarding behavior.
#' \code{SS2OM} returns a single sex (either male, female, or averaged biological parameters) and single fleet (aggregate selectivity and mortality,
#' no explicit discarding modeled). For either, the user still needs to parameterize most of the observation and
#' implementation portions. \code{SSMOM2OM} is the internal function that simplifies the MOM object to an OM object.
#' \code{plot_SS2OM} generates a markdown report to compare the OM and SS output.
#'
#' @aliases SS2OM SSMOM2OM
#' @param SSdir A folder with Stock Synthesis input and output files in it.
#' @param nsim The number of simulations to take for parameters with uncertainty (for OM@@cpars custom parameters).
#' @param proyears The number of projection years for MSE
#' @param reps The number of stochastic replicates within each simulation in the operating model.
#' @param maxF The maximum allowable F in the operating model.
#' @param seed The random seed for the operating model.
#' @param interval The interval at which management procedures will update the management advice in \link{multiMSE}, e.g., 1 = annual updates.
#' @param pstar The percentile of the sample of the management recommendation for the MP/MMP.
#' @param Obs The observation model (class Obs). These functions do not update implementation parameters.
#' @param Imp The implementation model (class Imp). These functions do not update implementation parameters.
#' @param silent Whether to silence messages to the console.
#' @param Name The name of the operating model
#' @param Source Reference to assessment documentation e.g. a url
#' @param ... Arguments to pass to \link[r4ss]{SS_output}.
#' @note Currently tested on r4ss version 1.38.1-40.0 and SS 3.30.14.
#' @return SS2MOM returns an object of class \linkS4class{MOM}. SS2OM returns an object of class \linkS4class{OM}.
#' @author Q. Huynh
#' @export
#' @seealso \link{SS2Data} \link{SS2DataMOM}
SS2MOM <- 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, silent = FALSE,
                   Name = "MOM generated by SS2MOM", Source = "No Source provided", ...) {

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

  # Create MOM object
  MOM <- suppressMessages(new("MOM"))
  MOM@Name <- Name
  MOM@Source <- Source
  MOM@reps <- reps
  MOM@seed <- seed
  MOM@maxF <- maxF
  MOM@interval <- interval

  if(!silent) message(replist$nsexes, "- sex and ", replist$nfishfleets, "- fleet model detected.")

  mainyrs <- replist$startyr:replist$endyr
  nyears <- length(mainyrs)
  MOM@proyears <- proyears
  allyears <- nyears + MOM@proyears
  MOM@nsim <- nsim
  MOM@pstar <- pstar

  output <- lapply(seq_len(replist$nsexes), SS_stock, replist = replist, mainyrs = mainyrs, nyears = nyears,
                   proyears = proyears, nsim = nsim, single_sex = length(replist$nsexes) == 1)

  MOM@Stocks <- lapply(output, getElement, "Stock") %>% structure(names = c("Female", "Male")[seq_len(replist$nsexes)])
  MOM@Fleets <- lapply(output, getElement, "Fleet") %>% structure(names = c("Female", "Male")[seq_len(replist$nsexes)])
  MOM@cpars <- lapply(output, getElement, "cpars") %>% structure(names = c("Female", "Male")[seq_len(replist$nsexes)])

  # Sample future recruitment
  Perr_proj <- exp(sample_recruitment(log(MOM@cpars[[1]][[1]]$Perr_y), proyears,
                                      replist$sigma_R_in, MOM@Stocks[[1]]@AC[1], seed))
  MOM@cpars <- lapply(MOM@cpars, function(x) lapply(x, function(xx) {xx$Perr_y <- cbind(xx$Perr_y, Perr_proj); return(xx)}))

  MOM@Obs <- lapply(seq_len(replist$nsexes), function(x) lapply(1:length(MOM@Fleets[[1]]), function(xx) return(Obs)) %>%
                      structure(names = names(MOM@Fleets[[1]]))) %>% structure(names = c("Female", "Male")[seq_len(replist$nsexes)])
  MOM@Imps <- lapply(seq_len(replist$nsexes), function(x) lapply(1:length(MOM@Fleets[[1]]), function(xx) return(Imp)) %>%
                       structure(names = names(MOM@Fleets[[1]]))) %>% structure(names = c("Female", "Male")[seq_len(replist$nsexes)])

  if(replist$nsexes == 2) {
    MOM@SexPars <- list(SSBfrom = matrix(c(1, 0), 2, 2, byrow = TRUE))
    MOM@Complexes <- list(c(1, 2))
  }

  # Catch Fracs
  CatchFrac <- lapply(MOM@cpars, function(x) vapply(x, function(xx) xx$Data@Cat[1, ncol(xx$Data@Cat)], numeric(1)))
  MOM@CatchFrac <- lapply(CatchFrac, function(x) matrix(x/sum(x), nrow = MOM@nsim, ncol = length(x), byrow = TRUE) %>%
                            structure(dimnames = list(NULL, names(MOM@Fleets[[1]])))) %>%
    structure(names = c("Female", "Male")[seq_len(replist$nsexes)])

  return(MOM)
}


#' Reads data Stock Synthesis file structure into a nested Data object analogous with multiMSE
#'
#' @description A function that uses the file location of a fitted SS3 model including input files to population
#' the various slots of an Data object.
#' @param SSdir A folder with Stock Synthesis input and output files in it. Alternatively,
#' @param age_M A vector of ages to average across to calculate a single value of natural mortality.
#' Currently, the Data object supports a single value of M for all ages. By default, \code{NULL} averages over all ages.
#' @param comp_partition Integer vector for selecting length/age observations that are retained (2), discarded (1), or both (0). By default, only retained
#' comps are used. If multiple codes are used, then comp matrix is the sum over all codes.
#' @param silent Logical. Suppress messages?
#' @param ... Arguments to pass to \link[r4ss]{SS_output}
#' @return A nested list of Data objects, with the first index by stock/sex and the second index by fleet.
#' @note Currently tested on r4ss version 1.38.1-41 and SS 3.30.14.
#'
#' Catches in \code{Data@@Cat} are the predicted sex-specific catch calculated from the SS output.
#' @author Q. Huynh
#' @export
#' @seealso \link{SS2MOM}
SS2DataMOM <- function(SSdir, age_M = NULL, comp_partition = 2, silent = FALSE, ...) {

  if(is.list(SSdir)) {
    replist <- SSdir
  } else {
    replist <- SS_import(SSdir, silent, ...)
  }
  if(replist$nseasons > 1) warning("Currently only supporting one season per year.")

  if(!silent) message(replist$nsexes, "- sex and ", replist$nfishfleets, "- fleet model detected.")

  mainyrs <- replist$startyr:replist$endyr
  nyears <- length(mainyrs)

  output <- lapply(seq_len(replist$nsexes), SS_stock, replist = replist, mainyrs = mainyrs, nyears = nyears,
                   proyears = 0, nsim = 1, single_sex = TRUE, partition = comp_partition, age_M = age_M)

  Data <- lapply(output, function(x) lapply(x$cpars, getElement, "Data")) %>%
    structure(names = c("Female", "Male")[seq_len(replist$nsexes)])

  return(Data)
}


#' @rdname SS2MOM
#' @export
plot_SS2MOM <- function(x, SSdir, gender = 1:2,
                       filename = "SS2MOM", dir = tempdir(), open_file = TRUE, silent = FALSE, ...) {
  if(missing(SSdir)) stop("SSdir not found.")
  
  if(inherits(x, "MOM")) {
    if(!silent) message("Generating multiHist object from MOM...")
    MOM <- x
    multiHist <- SimulateMOM(MOM, silent = silent)
  } else if(inherits(x, "multiHist")) {
    multiHist <- x
  } else {
    stop("Neither multiHist nor MOM 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", "SS2MOM.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))
}