R/MSE_functions.R

In MSEtool: Management Strategy Evaluation Toolkit

#' Check Convergence
#'
#' Have I undertaken enough simulations (nsim)? Has my MSE converged on stable
#' (reliable) performance metrics?
#'
#' Performance metrics are plotted against the number of simulations. Convergence diagnostics
#' are calculated over the last `ref.it` (default = 20) iterations. The convergence diagnostics are:
#' \enumerate{
#'   \item Is the order of the MPs stable over the last `ref.it` iterations?
#'   \item Is the average difference in performance statistic over the last `ref.it` iterations < `thresh`?
#' }
#'
#' By default three commonly used performance metrics are used:
#' \enumerate{
#'   \item Average Yield Relative to Reference Yield
#'   \item Probability Spawning Biomass is above 0.1BMSY
#'   \item Probability Average Annual Variability in Yield is < 20 per cent
#' }
#' Additional or alternative performance metrics objects can be supplied. Advanced users can develop their own performance metrics.
#'
#' @param MSEobj An MSE object of class \code{'MSE'}
#' @param PMs A character vector of names of the PM methods or a list of the PM methods
#' @param maxMP Maximum number of MPs to include in a single plot
#' @param thresh The convergence threshold. Maximum root mean square deviation over the last `ref.it` iterations
#' @param ref.it The number of iterations to calculate the convergence statistics. For example,
#' a value of 20 means convergence diagnostics are calculated over last 20 simulations
#' @param inc.leg Logical. Should the legend be displayed?
#' @param all.its Logical. Plot all iterations? Otherwise only (nsim-ref.it):nsim
#' @param nrow Numeric. Optional. Number of rows
#' @param ncol Numeric. Optional. Number of columns
#' @param silent Hide the messages printed in console?
#'
# #' @templateVar url checking-convergence
# #' @templateVar ref NULL
# #' @template userguide_link
#'
#' @return A table of convergence results for each MP
#' @examples
#' \dontrun{
#' MSE <- runMSE()
#' Converge(MSE)
#' }
#'
#' @author A. Hordyk
#' @export
#'
Converge <- function(MSEobj, PMs=c('Yield', 'P10', 'AAVY'), maxMP=15, thresh=0.5, ref.it=20,
                     inc.leg=FALSE, all.its=FALSE, nrow=NULL, ncol=NULL, silent=FALSE) {


  if(!methods::is(MSEobj, "MSE")) stop("MSEobj must be object of class 'MSE'", call.=FALSE)
  if(MSEobj@nMPs <2) stop("Converge requires more than 1 MP in the MSE object", call.=FALSE)

  if (!requireNamespace("ggrepel", quietly = TRUE)) {
    stop("Package \"ggrepel\" needed for this function to work. Please install it.",
         call. = FALSE)
  }

  nPMs <- length(PMs)

  if (mode(PMs) == 'character') {
    PMnames <- PMs
    PMs <- lapply(PMs, get)
  } else if (mode(PMs) == 'list') {
    PMnames <- 1:length(PMs)
  }

  if (is.null(ncol)) ncol <- floor(sqrt(nPMs))
  if (is.null(nrow)) nrow <- ceiling(nPMs/ncol)
  if (ncol * nrow < nPMs) stop("ncol x nrow must be > length(PMs)")

  if (MSEobj@nMPs > maxMP) {
    while (MSEobj@nMPs %% maxMP == 1) {
      maxMP <- maxMP+1
    }
    nplots <- ceiling(MSEobj@nMPs /maxMP)
  } else{
    nplots <- 1
  }

  nsim <- MSEobj@nsim
  if (nsim-ref.it < 1) {
    ref.it.new <- nsim - 1
    if (!silent) message("nsim (", nsim, ") - ref.it (", ref.it,") < 1. Setting ref.it to ", ref.it.new, "\n")
    ref.it <- ref.it.new
  }

  if (!silent) message ("Checking if order of MPs is changing in last ", ref.it, " iterations\n")
  if (!silent) message("Checking average difference in PM over last ", ref.it, " iterations is > ", thresh, "\n")

  SwitchOrd <- vector(mode = "list", length = nPMs)
  NonCon <- vector(mode = "list", length = nPMs)
  PMName <- vector("character", length=nPMs)

  getPalette <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                                   "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))
  st <- 1
  end <- min(maxMP, MSEobj@nMPs)
  it <- MP <- NULL # cran check hacks

  for (nplot in 1:nplots) {

    if (!silent) message('Plotting MPs ', st, ' - ', end)
    subMSE <- Sub(MSEobj, MPs=MSEobj@MPs[st:end])

    nMPs <-  subMSE@nMPs
    values <- rep(c("solid", "dashed", "dotted"),  nMPs)
    values <- values[1:nMPs]

    plist <- list()
    for (xx in 1:nPMs) {
      PMval <- PMs[[xx]](subMSE)
      PMName[xx] <- PMval@Name
      PMval@Prob[!is.finite(PMval@Prob)] <- 0

      cum_mean <- apply(PMval@Prob, 2, cumsum)/apply(PMval@Prob, 2, seq_along) * 100
      vals <- as.vector(cum_mean)
      mp <- rep(subMSE@MPs, each=subMSE@nsim)

      df <- data.frame(it=1:subMSE@nsim, vals, MP=mp, name=PMval@Name)
      if (!all.its) df <- subset(df, it %in% (nsim-ref.it+1):nsim)

      p <- ggplot2::ggplot(df, ggplot2::aes(x=it, y=vals, color=MP, linetype=MP)) +
        ggplot2::scale_linetype_manual(values = values) +
        ggplot2::scale_color_manual(values=getPalette(nMPs)) +
        ggplot2::geom_line() +
        ggplot2::theme_classic() +
        ggplot2::labs(x="# Simulations", y=PMval@Caption) +
        ggplot2::ggtitle(PMval@Name) +
        ggplot2::geom_vline(xintercept=MSEobj@nsim-ref.it+1, color="darkgray", linetype="dashed") +
        ggplot2::geom_vline(xintercept=MSEobj@nsim, color="darkgray", linetype="dashed") +
        ggplot2::coord_cartesian(xlim = c(min(df$it), max(df$it)))
      if (!inc.leg) p <- p + ggplot2::theme(legend.position = "none")

      # check positions & convergence
      ords <- apply(cum_mean[(MSEobj@nsim-ref.it+1):MSEobj@nsim,], 1, order, decreasing=FALSE)
      rownames(ords) <- subMSE@MPs
      mat <- matrix(subMSE@MPs[ords], nrow=subMSE@nMPs, ncol=ref.it)
      tab <- table(unlist(apply(mat, 1, unique)))
      SwitchOrd[[xx]] <- append(SwitchOrd[[xx]], rownames(tab)[which(tab > 1)])
      NonCon[[xx]] <- append(NonCon[[xx]], subMSE@MPs[apply(cum_mean, 2, Chk, MSEobj=MSEobj, thresh=thresh, ref.it)])

      noncoverg <- unique(c(SwitchOrd[[xx]], NonCon[[xx]]))
      if (length(noncoverg)>0) {
        df2 <- subset(df, MP %in% noncoverg)
        if (dim(df2)[1] > 0) {
          p <- p +
            ggrepel::geom_text_repel(
              data = subset(df2, it == max(it)),
              ggplot2::aes(label = MP),
              size = 4,
              segment.color = "grey",
              direction = "both",
              show.legend = FALSE)
        }

      }

      plist[[xx]] <- p

    }

    if (inc.leg) join_plots(plist, nrow=nrow, ncol=ncol, position="right")
    if (!inc.leg) {
      if (!requireNamespace("gridExtra", quietly = TRUE)) {
        stop("Package \"gridExtra\" needed for this function to work. Please install it.",
             call. = FALSE)
      }
      gridExtra::grid.arrange(grobs=plist, nrow=nrow, ncol=ncol)
    }
    st <- st + maxMP
    end <- min(end + maxMP, MSEobj@nMPs)

  }

  ## Report Performance
  resultsTab <- matrix(NA,nrow=2,ncol=nPMs)
  rownames(resultsTab) <- c("MPs Not Converging", "Unstable MPs")
  colnames(resultsTab) <- PMnames

  for (x in 1:nPMs) {
    SwitchOrds <- SwitchOrd[[x]]
    NonCons <- NonCon[[x]]
    # save results in table
    resultsTab[1,x] <- length(NonCons)
    resultsTab[2,x] <- length(SwitchOrds)

    if (!silent) {
      if (length(SwitchOrds)>0  | length(NonCons)>0) {
        message("\n", PMName[x], "\n")
        if (length(SwitchOrds)>0) {
          message("Order over last ", ref.it, ' iterations is not consistent for:\n ', paste(SwitchOrds, ""), ' \n')
        }
        if (length(NonCons)>0) {
          message("Mean difference over last ", ref.it, " iterations is > ", thresh, " for:\n",
                  paste(NonCons, ""), '\n')
        }
      }
    }
  }
  return(resultsTab)
}





Chk <- function(X, MSEobj, thresh, ref.it) {
  L <- length(X)
  Y <- min(MSEobj@nsim, ref.it) + 1
  x <- X[(L-Y):L]

  # root mean square deviation in last ref.it iterations is greater than thresh
  sqrt((sum((x-mean(x))^2))/(L-1)) > thresh
}


#   nm <- MSEobj@nMPs
#   nsim <- MSEobj@nsim
#   proyears <- MSEobj@proyears
#
#   Yd <- CumlYd <- array(NA, c(nm, nsim))
#   P10 <- CumlP10 <- array(NA, c(nm, nsim))
#   P50 <- CumlP50 <- array(NA, c(nm, nsim))
#   P100 <- CumlP100 <- array(NA, c(nm, nsim))
#   POF <- CumlPOF <- array(NA, c(nm, nsim))
#   yind <- max(MSEobj@proyears - 4, 1):MSEobj@proyears
#   RefYd <- MSEobj@OM$RefY
#
#   for (m in 1:nm) {
#     Yd[m, ] <- round(apply(MSEobj@C[, m, yind], 1, mean, na.rm = T)/RefYd *  100, 1)
#     POF[m, ] <- round(apply(MSEobj@F_FMSY[, m, ] >= 1, 1, sum, na.rm = T)/proyears * 100, 1)
#     P10[m, ] <- round(apply(MSEobj@B_BMSY[, m, ] <= 0.1, 1, sum, na.rm = T)/proyears * 100, 1)
#     P50[m, ] <- round(apply(MSEobj@B_BMSY[, m, ] <= 0.5, 1, sum, na.rm = T)/proyears * 100, 1)
#     P100[m, ] <- round(apply(MSEobj@B_BMSY[, m, ] <= 1, 1, sum, na.rm = T)/proyears * 100, 1)
#     CumlYd[m, ] <- cumsum(Yd[m, ])/seq_along(Yd[m, ])  #/ mean(Yd[m,], na.rm=TRUE)
#     CumlPOF[m, ] <- cumsum(POF[m, ])/seq_along(POF[m, ])  # / mean(POF[m,], na.rm=TRUE)
#     CumlP10[m, ] <- cumsum(P10[m, ])/seq_along(P10[m, ])  # / mean(P10[m,], na.rm=TRUE)
#     CumlP50[m, ] <- cumsum(P50[m, ])/seq_along(P50[m, ])  # / mean(P50[m,], na.rm=TRUE)
#     CumlP100[m, ] <- cumsum(P100[m, ])/seq_along(P100[m, ])  # / mean(P100[m,], na.rm=TRUE)
#   }
#
#
#
#   # CumlYd[is.nan(CumlYd)] <- 1 CumlPOF[is.nan(CumlPOF)] <- 1
#   # CumlP10[is.nan(CumlP10)] <- 1 CumlP50[is.nan(CumlP50)] <- 1
#   # CumlP100[is.nan(CumlP100)] <- 1
#
#   if (Plot) {
#     op <- par(mfrow = c(2, 3), cex.axis = 1.5, cex.lab = 1.7, oma = c(1, 1, 0, 0), mar = c(5, 5, 1, 1), bty = "l")
#     matplot(t(CumlYd), type = "l", xlab = "Iteration", ylab = "Rel. Yield")
#     matplot(t(CumlPOF), type = "l", xlab = "Iteration", ylab = "Prob. F/FMSY > 1")
#     matplot(t(CumlP10), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 0.1")
#     matplot(t(CumlP50), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 0.5")
#     matplot(t(CumlP100), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 1")
#
#   }
#
#   Chk <- function(X) {
#     # checks if difference in last 20 iterations is greater than thresh
#     L <- length(X)
#     Y <- 1:min(nsim, 20)
#
#     # return(mean(abs((X[L-1:10] - X[L]))/X[L], na.rm=TRUE) > thresh)
#     return(mean(abs((X[L - Y] - X[L])), na.rm = TRUE) > thresh)
#   }
#
#   NonCon <- sort(unique(c(which(apply(CumlYd, 1, Chk)),
#                           which(apply(CumlPOF, 1, Chk)),
#                           which(apply(CumlP10, 1, Chk)),
#                           which(apply(CumlP50, 1, Chk)),
#                           which(apply(CumlP100, 1, Chk)))))
#
#   if (length(NonCon) > 0) {
#     if (Plot) {
#       plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "",
#            ylab = "")
#       text(0.5, 0.5, "Some MPs have not converged", cex = 1)
#       # ev.new()
#       par(mfrow = c(2, 3), cex.axis = 1.5, cex.lab = 1.7, oma = c(1,
#                                                                   1, 0, 0), mar = c(5, 5, 1, 1), bty = "l")
#       if (length(NonCon) > 1) {
#         matplot(t(CumlYd[NonCon, ]), type = "l", xlab = "Iteration",  ylab = "Rel. Yield", lwd = 2, lty=1:length(NonCon))
#       } else matplot((CumlYd[NonCon, ]), type = "l", xlab = "Iteration",  ylab = "Rel. Yield", lwd = 2, lty=1:length(NonCon))
#       if (length(NonCon) > 1) {
#         matplot(t(CumlPOF[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. F/FMSY > 1", lwd = 2, lty=1:length(NonCon))
#       } else matplot((CumlPOF[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. F/FMSY > 1", lwd = 2, lty=1:length(NonCon))
#       if (length(NonCon) > 1) {
#         matplot(t(CumlP10[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 0.1", lwd = 2, lty=1:length(NonCon))
#       } else matplot((CumlP10[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 0.1", lwd = 2, lty=1:length(NonCon))
#       if (length(NonCon) > 1) {
#         matplot(t(CumlP50[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 0.5", lwd = 2, lty=1:length(NonCon))
#       } else matplot((CumlP50[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 0.5", lwd = 2, lty=1:length(NonCon))
#       if (length(NonCon) > 1) {
#         matplot(t(CumlP100[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 1", lwd = 2, lty=1:length(NonCon))
#       } else matplot((CumlP100[NonCon, ]), type = "l", xlab = "Iteration", ylab = "Prob. B/BMSY < 1", lwd = 2, lty=1:length(NonCon))
#
#       legend(nsim * 1.25, 50, legend = MSEobj@MPs[NonCon], col = 1:length(NonCon),
#              bty = "n", xpd = NA, lty = 1:length(NonCon), lwd = 2, cex = 1.25)
#     }
#
#     message("Some MPs may not have converged in ", nsim, " iterations (threshold = ", thresh, "%)")
#     message("MPs are: ", paste(MSEobj@MPs[NonCon], " "))
#     message("MPs #: ", paste(NonCon, " "))
#     return(data.frame(Num = NonCon, MP = MSEobj@MPs[NonCon]))
#   }
#   if (length(NonCon) == 0) {
#     if (Plot) {
#       plot(c(0, 1), c(0, 1), type = "n", axes = FALSE, xlab = "",
#            ylab = "")
#       text(0.5, 0.5, "All MPs converged", cex = 1)
#     }
#     message("All MPs appear to have converged in ", nsim, " iterations (threshold = ",
#             thresh, "%)")
#   }
#   par(op)
# }







# Kobe plot Kplot<-function(MSEobj,maxsim=60,nam=NA){
# nr<-floor((MSEobj@nMPs)^0.5) nc<-ceiling((MSEobj@nMPs)/nr)

# FMSYr<-quantile(MSEobj@F_FMSY,c(0.001,0.90),na.rm=T)
# BMSYr<-quantile(MSEobj@B_BMSY,c(0.001,0.975),na.rm=T)

# #dev.new2(width=nc*3,height=nr*3.6)
# par(mfrow=c(nr,nc),mai=c(0.45,0.45,0.45,0.01),omi=c(0.45,0.3,0.35,0.01))

# colsse<-rainbow(MSEobj@proyears,start=0.63,end=0.95)[1:MSEobj@proyears]
# colsse<-makeTransparent(colsse,95)

# for(mm in 1:MSEobj@nMPs){
# plot(c(MSEobj@B_BMSY[1,mm,1],MSEobj@B_BMSY[1,mm,2]),
# c(MSEobj@F_FMSY[1,mm,1],MSEobj@F_FMSY[1,mm,2]),xlim=BMSYr,ylim=FMSYr,
# col=colsse[1],type='l')

# OO<-round(sum(MSEobj@B_BMSY[,mm,MSEobj@proyears]<1&MSEobj@F_FMSY[,mm,MSEobj@proyears]>1,na.rm=T)/MSEobj@nsim*100,1)
# OU<-round(sum(MSEobj@B_BMSY[,mm,MSEobj@proyears]>1&MSEobj@F_FMSY[,mm,MSEobj@proyears]>1,na.rm=T)/MSEobj@nsim*100,1)
# UO<-round(sum(MSEobj@B_BMSY[,mm,MSEobj@proyears]<1&MSEobj@F_FMSY[,mm,MSEobj@proyears]<1,na.rm=T)/MSEobj@nsim*100,1)
# UU<-round(sum(MSEobj@B_BMSY[,mm,MSEobj@proyears]>1&MSEobj@F_FMSY[,mm,MSEobj@proyears]<1,na.rm=T)/MSEobj@nsim*100,1)

# #alp<-80
# #polygon(c(1,-1000,-1000,1),c(1,1,1000,1000),col=makeTransparent('orange',alp),border=makeTransparent('orange',alp))
# #polygon(c(1,1000,1000,1),c(1,1,1000,1000),col=makeTransparent('yellow',alp),border=makeTransparent('yellow',alp))
# #polygon(c(1,-1000,-1000,1),c(1,1,-1000,-1000),col=makeTransparent('yellow',alp),border=makeTransparent('yellow',alp))
# #polygon(c(1,1000,1000,1),c(1,1,-1000,-1000),col=makeTransparent('green',alp),border=makeTransparent('yellow',alp))


# abline(h=1,col='grey',lwd=3) abline(v=1,col='grey',lwd=3)
# #abline(v=c(0.1,0.5),col='grey',lwd=2) rng<-1:min(maxsim,MSEobj@nsim)
# for(i in rng){ for(y in 1:(MSEobj@proyears-1)){
# lines(c(MSEobj@B_BMSY[i,mm,y],MSEobj@B_BMSY[i,mm,y+1]),
# c(MSEobj@F_FMSY[i,mm,y],MSEobj@F_FMSY[i,mm,y+1]),
# col=colsse[y],lwd=1.6) } }

# points(MSEobj@B_BMSY[rng,mm,1],MSEobj@F_FMSY[rng,mm,1],pch=19,cex=0.8,col=colsse[1])
# points(MSEobj@B_BMSY[rng,mm,MSEobj@proyears],MSEobj@F_FMSY[rng,mm,MSEobj@proyears],pch=19,cex=0.8,col=colsse[MSEobj@proyears])

# if(mm==1)legend('right',c('Start','End'),bty='n',text.col=c(colsse[1],colsse[MSEobj@proyears]),pch=19,col=c(colsse[1],colsse[MSEobj@proyears]))
# legend('topleft',paste(OO,'%',sep=''),bty='n',text.font=2)
# legend('topright',paste(OU,'%',sep=''),bty='n',text.font=2)
# legend('bottomleft',paste(UO,'%',sep=''),bty='n',text.font=2)
# legend('bottomright',paste(UU,'%',sep=''),bty='n',text.font=2)

# mtext(MSEobj@MPs[mm],3,line=0.45) }
# mtext('B/BMSY',1,outer=T,line=1.4) mtext('F/FMSY',2,outer=T,line=0.2)
# if(is.na(nam))mtext(deparse(substitute(MSEobj)),3,outer=T,line=0.25,font=2)
# if(!is.na(nam))mtext(MSEobj@Name,3,outer=T,line=0.25,font=2) }





# Value of information analysis
# Value of information




# Manipulation of MSE Object
# --------------------------------------------------- Subset the MSE
# object by particular MPs (either MP number or name), or particular
# simulations

#' Utility functions for MSE objects
#'
#' @param MSEobj A `MSE` object
# #' @param MSEobj A MSE object. For `updateMSE`, a MSE object from a previous version of
# #' DLMtool. Also works with Stock, Fleet, Obs, Imp, and Data objects.
# #' @param MSEobjs A list of MSE objects. Must all have identical operating
# #' model and MPs. MPs which don't appear in all MSE objects will be dropped.
#' @return An object of class \code{MSE}
# #' @examples
# #' # An example of joinMSE
# #' \dontrun{
# #' OM1 <- MSEtool::testOM
# #' MSE1 <- runMSE(OM1)
# #' OM2 <- OM1
# #' OM2@seed <- OM1@seed + 1
# #' MSE2 <- runMSE(OM2)
# #' MSE <- joinMSE(list(MSE1, MSE2))
# #' MSE@nsim
# #' }
#' @author A. Hordyk
#' @describeIn checkMSE Check that an MSE object includes all slots in the latest version of DLMtool
# #' Use `updateMSE` to update the MSE object
#' @export checkMSE
checkMSE <- function(MSEobj) {
  nms <- slotNames(MSEobj)
  errs <- NULL
  for (x in seq_along(nms)) {
    chk <- try(slot(MSEobj, nms[x]), silent=TRUE)
    if ("try-error" %in% class(chk)) errs <- c(errs, x)
  }
  if (length(errs) > 0) {
    message("MSE object slots not found: ", paste(nms[errs], ""))
    stop("slot names of MSEobj don't match MSE object class. Try use `updateMSE`", call.=FALSE)
    return(FALSE)
  }
  return(TRUE)
}

#' Subset MSE object by management procedure (MP) or simulation.
#'
#' Subset the MSE object by particular MPs (either MP number or name), or
#' particular simulations, or a subset of the projection years (e.g., 1: <
#' projection years).
#'
#' @param MSEobj A MSE object.
#' @param MPs A vector MPs names or MP numbers to subset the MSE object.
#' Defaults to all MPs.
#' @param sims A vector of simulation numbers to subset the MSE object. Can
#' also be a logical vector. Defaults to all simulations.
#' @param years A numeric vector of projection years. Should start at 1 and
#' increase by one to some value equal or less than the total number of
#' projection years.
#'
# #' @templateVar url subsetting-the-mse-object
# #' @templateVar ref NULL
# #' @template userguide_link
#'
#' @author A. Hordyk
#' @examples
#' \dontrun{
#' MSE <- runMSE()
#' MSE_1 <- Sub(MSE, MPs=1:2)
#' MSE_1@MPs
#' MSE_2 <- Sub(MSE, sims=1:10)
#' MSE_2@nsim
#' }
#' @seealso \link{SubOM} for OM components and \link{SubCpars} for subsetting by simulation and projection years.
#' @export Sub
Sub <- function(MSEobj, MPs = NULL, sims = NULL, years = NULL) {

  checkMSE(MSEobj) # check that MSE object contains all slots

  Class <- class(MPs)
  if (Class == "NULL") subMPs <- MSEobj@MPs
  if (Class == "integer" | Class == "numeric") subMPs <- MSEobj@MPs[as.integer(MPs)]
  if (Class == "character") subMPs <- MPs
  if (Class == "factor") subMPs <- as.character(MPs)
  subMPs <- subMPs[!is.na(subMPs)]
  subMPs <- unique(subMPs)

  SubMPs <- match(subMPs, MSEobj@MPs)  #  which(MSEobj@MPs %in% subMPs)
  if (any(is.na(SubMPs))) {
    missing <- subMPs[is.na(SubMPs)]
    stop(paste0(missing, collapse=','), ' not found in MSE object', call.=FALSE)
  }

  not <- (subMPs %in% MSEobj@MPs)  # Check for MPs misspelled
  ind <- which(not == FALSE)
  newMPs <- MSEobj@MPs[SubMPs]
  if (length(SubMPs) < 1) stop("MPs not found")
  if (length(ind > 0)) {
    message("Warning: MPs not found - ", paste0(subMPs[ind], " "))
    message("Subsetting by MPs: ", paste0(newMPs, " "))
  }

  ClassSims <- class(sims)
  if (ClassSims == "NULL")  SubIts <- 1:MSEobj@nsim
  if (ClassSims == "integer" | ClassSims == "numeric") {
    # sims <- 1:min(MSEobj@nsim, max(sims))
    SubIts <- as.integer(sims)
  }
  if (ClassSims == "logical")  SubIts <- which(sims)
  nsim <- length(SubIts)

  ClassYrs <- class(years)
  AllNYears <- MSEobj@proyears
  if (ClassYrs == "NULL")
    Years <- 1:AllNYears
  if (ClassYrs == "integer" | ClassYrs == "numeric")
    Years <- years
  if (max(Years) > AllNYears)
    stop("years exceeds number of years in MSE")
  if (min(Years) <= 0)
    stop("years must be positive")
  if (min(Years) != 1) {
    message("Not starting from first year. Are you sure you want to do this? Probably a bad idea!")
  }
  if (!all(diff(Years) == 1))
    stop("years are not consecutive")
  if (length(Years) <= 1)
    stop("You are going to want more than 1 projection year")
  MSEobj@proyears <- max(Years)

  SubBioEco <- MSEobj@BioEco
  for (i in 1:length(SubBioEco)) {
    SubBioEco[[i]] <-  SubBioEco[[i]][SubIts, SubMPs,  Years, drop=FALSE]
  }

  SubRefPoint <- MSEobj@RefPoint
  for (i in 1:length(SubRefPoint)) {
    if ('array' %in% class(SubRefPoint[[i]])) {
      DD <- dim(SubRefPoint[[i]])
      if (length(DD) ==3) {
        SubRefPoint[[i]] <- SubRefPoint[[i]][SubIts, SubMPs, c(1:MSEobj@nyears, MSEobj@nyears+Years), drop=FALSE]
      } else if (length(DD) ==4) {
        SubRefPoint[[i]] <- SubRefPoint[[i]][SubIts, SubMPs, , c(1:MSEobj@nyears, MSEobj@nyears+Years), drop=FALSE]
      } else {
        warning('Cannot subset MSEobj@RefPoint$', names( MSEobj@RefPoint)[[i]])
        SubRefPoint[[i]] <- SubRefPoint[[i]]
      }

    } else if ('list' %in% class(SubRefPoint[[i]])) {
        for (j in names(SubRefPoint[[i]])) {
          DD <- dim(SubRefPoint[[i]][[j]])
          if (length(DD)==2) {
            SubRefPoint[[i]][[j]] <- SubRefPoint[[i]][[j]][SubIts, c(1:MSEobj@nyears, MSEobj@nyears+Years), drop=FALSE]
          } else if (length(DD)==3) {
            SubRefPoint[[i]][[j]] <- SubRefPoint[[i]][[j]][SubIts, , c(1:MSEobj@nyears, MSEobj@nyears+Years), drop=FALSE]
          } else {
            warning('Cannot subset MSEobj@RefPoint$', names( MSEobj@RefPoint)[[i]], '$', j)
            SubRefPoint[[i]][[j]] <- SubRefPoint[[i]][[j]]
          }

        }
    }

  }
  
  SubSPR <- lapply(MSEobj@SPR, function(x) x[SubIts, SubMPs, Years, drop=FALSE])
  
  subMSElist <- MSEobj@PPD[SubMPs] # doesn't subset Data by years or simulations
  
  if(length(MSEobj@Misc$extended)) {
    
    Subextended <- lapply(MSEobj@Misc$extended, function(x) {
      x[SubIts, , SubMPs, c(1:MSEobj@nyears, MSEobj@nyears+Years), , drop=FALSE]
    })

  } else {
    Subextended <- list()
  }
  
  # Numbers - dimensions changed in v3.6.2
  if (length(dim(MSEobj@N)) == 3) {
    N <- MSEobj@N[SubIts, SubMPs, Years, drop=FALSE]  
  } else {
    N <- MSEobj@N[SubIts, , SubMPs, Years, , drop=FALSE]
  }
  
  subMSE <- new("MSE",
                Name = MSEobj@Name,
                nyears=MSEobj@nyears,
                proyears=MSEobj@proyears,
                nMPs=length(newMPs),
                MPs=newMPs,
                nsim=nsim,
                OM=MSEobj@OM[SubIts,, drop=FALSE],
                Obs=MSEobj@Obs[SubIts,, drop=FALSE],
                SB_SBMSY=MSEobj@SB_SBMSY[SubIts, SubMPs,  Years, drop=FALSE],
                F_FMSY=MSEobj@F_FMSY[SubIts, SubMPs,  Years, drop=FALSE],
                N=N,
                B=MSEobj@B[SubIts, SubMPs,  Years, drop=FALSE],
                SSB=MSEobj@SSB[SubIts, SubMPs,  Years, drop=FALSE],
                VB=MSEobj@VB[SubIts, SubMPs,  Years, drop=FALSE],
                FM=MSEobj@FM[SubIts, SubMPs,  Years, drop=FALSE],
                SPR=SubSPR,
                Catch=MSEobj@Catch[SubIts, SubMPs,  Years, drop=FALSE],
                Removals=MSEobj@Removals[SubIts, SubMPs,  Years, drop=FALSE],
                Effort=MSEobj@Effort[SubIts, SubMPs,  Years, drop=FALSE],
                TAC=MSEobj@TAC[SubIts, SubMPs,  Years, drop=FALSE],
                TAE=MSEobj@TAE[SubIts, SubMPs,  Years, drop=FALSE],
                BioEco=SubBioEco,
                RefPoint=SubRefPoint,
                CB_hist=MSEobj@CB_hist[SubIts,  , drop=FALSE],
                FM_hist=MSEobj@FM_hist[SubIts,  , drop=FALSE],
                SSB_hist=MSEobj@SSB_hist[SubIts, , drop=FALSE],
                Hist=MSEobj@Hist,
                PPD=subMSElist,
                Misc=list(extended = Subextended))
  
  attr(subMSE, "version") <- packageVersion("MSEtool")
  attr(subMSE, "date") <- date()
  attr(subMSE, "R.version") <- R.version
  
  subMSE
}



# Evaluate Peformance of MPs
# --------------------------------------------------- Function examines
# how consistently an MP outperforms another.


# #' How dominant is an MP?
# #'
# #' The DOM function examines how consistently an MP outperforms another. For
# #' example DCAC might provide higher yield than AvC on average but outperforms
# #' AvC in less than half of simulations.
# #'
# #'
# #' @param MSEobj An object of class 'MSE'
# #' @param MPtg A character vector of management procedures for cross
# #' examination
# #' @return A matrix of performance comparisons length(MPtg) rows by MSE@nMPs
# #' columns
# #' @author A. Hordyk
# #' @export DOM
# DOM <- function(MSEobj, MPtg = NA) {
#   if (any(is.na(MPtg)))
#     MPtg <- MSEobj@MPs
#   proyears <- MSEobj@proyears
#   nMP <- MSEobj@nMPs
#   nsim <- MSEobj@nsim
#   ind <- which(MSEobj@MPs %in% MPtg)
#   MPr <- which(!(MSEobj@MPs %in% MPtg))
#   yind <- max(MSEobj@proyears - 4, 1):MSEobj@proyears
#   y1 <- 1:(MSEobj@proyears - 1)
#   y2 <- 2:MSEobj@proyears
#   Mat <- matrix(0, nrow = length(MPtg), ncol = nMP)
#   rownames(Mat) <- MPtg
#   colnames(Mat) <- MSEobj@MPs
#   POF <- P100 <- YieldMat <- IAVmat <- Mat
#   for (X in 1:length(MPtg)) {
#     # Overfishing (F > FMSY)
#     ind1 <- as.matrix(expand.grid(1:nsim, ind[X], 1:proyears))
#     ind2 <- as.matrix(expand.grid(1:nsim, 1:nMP, 1:proyears))
#     t1 <- apply(array(MSEobj@F_FMSY[ind1] > 1, dim = c(nsim, 1, proyears)),
#                 c(1, 2), sum, na.rm = TRUE)
#     t2 <- apply(array(MSEobj@F_FMSY[ind2] > 1, dim = c(nsim, nMP, proyears)),
#                 c(1, 2), sum, na.rm = TRUE)
#     POF[X, ] <- round(apply(matrix(rep(t1, nMP), nrow = nsim) < t2,
#                             2, sum)/nsim * 100, 0)
#     # B < BMSY
#     t1 <- apply(array(MSEobj@B_BMSY[ind1] < 1, dim = c(nsim, 1, proyears)),
#                 c(1, 2), sum, na.rm = TRUE)
#     t2 <- apply(array(MSEobj@B_BMSY[ind2] < 1, dim = c(nsim, nMP, proyears)),
#                 c(1, 2), sum, na.rm = TRUE)
#     P100[X, ] <- round(apply(matrix(rep(t1, nMP), nrow = nsim) < t2,
#                              2, sum, na.rm = TRUE)/nsim * 100, 0)
#     # Relative yield in last 5 years
#     ind1 <- as.matrix(expand.grid(1:nsim, ind[X], yind))
#     ind2 <- as.matrix(expand.grid(1:nsim, 1:nMP, yind))
#     t1 <- apply(array(MSEobj@C[ind1], dim = c(nsim, 1, length(yind))),
#                 c(1, 2), sum, na.rm = TRUE)
#     t2 <- apply(array(MSEobj@C[ind2], dim = c(nsim, nMP, length(yind))),
#                 c(1, 2), sum, na.rm = TRUE)
#     YieldMat[X, ] <- round(apply(matrix(rep(t1, nMP), nrow = nsim) >
#                                    t2, 2, sum, na.rm = TRUE)/nsim * 100, 0)
#     # interannual variation in catch
#     ind1 <- as.matrix(expand.grid(1:nsim, ind[X], y1))
#     ind2 <- as.matrix(expand.grid(1:nsim, ind[X], y2))
#     AAVY1 <- apply(array(((MSEobj@C[ind1] - MSEobj@C[ind2])^2)^0.5,
#                          dim = c(nsim, 1, length(y1))), 1, mean, na.rm = T)/apply(array(MSEobj@C[ind2],
#                                                                                         dim = c(nsim, 1, length(y1))), 1, mean, na.rm = T)
#     ind1 <- as.matrix(expand.grid(1:nsim, 1:nMP, y1))
#     ind2 <- as.matrix(expand.grid(1:nsim, 1:nMP, y2))
#     AAVY2 <- apply(array(((MSEobj@C[ind1] - MSEobj@C[ind2])^2)^0.5,
#                          dim = c(nsim, nMP, length(y1))), c(1, 2), mean, na.rm = T)/apply(array(MSEobj@C[ind2],
#                                                                                                 dim = c(nsim, nMP, length(y1))), c(1, 2), mean, na.rm = T)
#     IAVmat[X, ] <- round(apply(matrix(rep(AAVY1, nMP), nrow = nsim) <
#                                  AAVY2, 2, sum, na.rm = TRUE)/nsim * 100, 0)
#   }
#   out <- list()
#   out$POF <- POF
#   out$P100 <- P100
#   out$Yd <- YieldMat
#   out$AAVY <- IAVmat
#   return(out)
# }
#

#' Determine dominate MPs
#'
#' MPs that perform worse than comparable MPs across all performance metrics are considered 'dominated' as
#' other options are always preferable.
#'
#' The `Dom` function compares the probabilities calculated in the performance metric
#' (`PM`) functions and determines the MPs that have a lower probability across all PMs compared
#' to other MPs of the same management type (e.g., size limit, TAC, etc).
#' Consequently, it is important that all `PM` functions are constructed so that higher probabilities = better performance
#' (e.g, `PNOF` is the probability of NOT overfishing)
#'
#' @param MSEobj An object of class `MSE`
#' @param ...  Names of Performance Metrics (PMs), or other arguments to `TradePlot`.
#' First PM is recycled if number of PMs is not even
#' @param PMlist Optional list of PM names. Overrides any supplied in ... above
#' @param Refs An optional named list (matching the PM names) with numeric values to override the default `Ref` values.
#' @param Yrs An optional named list (matching the PM names) with numeric values to override the default `Yrs` values.
#'
#' @return A named list of length 2 with a character vector of non-dominated MPs in `MPs` and
#' a data.frame of dominated MPs and the names of the relevant dominated MPs in `DomMPs`
#' @export
#' @author A. Hordyk
#' @examples
#' \dontrun{
#' MSE <- runMSE(MPs=NA) # run all MPs
#' Nondom <- Dom(MSE, "P10", "LTY", "PNOF")
#' # Non-dominated MPs
#' Nondom$MPs
#'
#' # Dominated MPs
#' Nondom$DomMPs
#'
#' }
#'
Dom <- function(MSEobj, ..., PMlist=NULL, Refs=NULL, Yrs=NULL) {
  if (!methods::is(MSEobj, 'MSE')) stop("Object must be class `MSE`", call.=FALSE)
  if (is.null(PMlist)) {
    PMlist <- unlist(list(...))
  } else {
    PMlist <- unlist(PMlist)
  }

  if (!methods::is(PMlist,'character')) stop("Must provide names of PM methods")

  runPM <- vector("list", length(PMlist))
  for (X in 1:length(PMlist)) {
    ref <- Refs[[PMlist[X]]]
    yrs <- Yrs[[PMlist[X]]]
    if (is.null(ref)) {
      if (is.null(yrs)) {
        runPM[[X]] <- eval(call(PMlist[X], MSEobj))
      } else {
        runPM[[X]] <- eval(call(PMlist[X], MSEobj, Yrs=yrs))
      }
    } else {
      if (is.null(yrs)) {
        runPM[[X]] <- eval(call(PMlist[X], MSEobj, Ref=ref))
      } else {
        runPM[[X]] <- eval(call(PMlist[X], MSEobj, Ref=ref, Yrs=yrs))
      }
    }
  }

  # Create Table
  DF <- Required(MSEobj@MPs, TRUE) %>% data.frame(., stringsAsFactors = FALSE)
  DF <- dplyr::left_join(DF, MPtype(MSEobj@MPs), by="MP")

  DF3 <- lapply(runPM, slot, name='Mean') %>% do.call("cbind", .)
  colnames(DF3) <- PMlist
  DF3 <- data.frame(MP=runPM[[1]]@MPs, DF3, stringsAsFactors = FALSE)
  DF <- dplyr::left_join(DF, DF3, by="MP")

  DomList <- list()
  for (i in 1:MSEobj@nMPs) {
    # ind <- which(DF$DataClass == DF$DataClass[i] & DF$Rec == DF$Rec[i] & DF$Type ==DF$Type[i])
    ind <- which(DF$Rec == DF$Rec[i] & DF$Type ==DF$Type[i]) # group MPs by Rec Type
    ind <- ind[!ind==i]
    if (length(ind)>0) {
      df1 <- DF[i,] %>%
        dplyr::select(-MP, -Data, -DataClass, -Type, -Recs) %>% unlist()
      m1 <- matrix(df1, nrow=length(ind), ncol=length(df1), byrow=TRUE)
      df2 <- DF[ind,] %>%
        dplyr::select(-MP, -Data, -DataClass, -Type, -Recs) %>% as.matrix()
      ind2 <- which(rowSums(m1 < df2) ==0)
      if (length(ind2)>0) {
        DomList[[i]] <- data.frame(DominatedMPs=DF$MP[ind[ind2]],
                                   By=DF$MP[i], stringsAsFactors = FALSE)
      }
    }
  }
  DomDF <- do.call("rbind", DomList) %>% as.data.frame()
  DomMPs <- unique(DomDF$DominatedMPs)
  NonDom <- MSEobj@MPs[!MSEobj@MPs %in% DomMPs]
  DomList2 <- list()
  for (i in seq_along(DomMPs)) {
    By <- DomDF %>% dplyr::filter(DominatedMPs ==DomMPs[i]) %>%
      dplyr::select(By) %>% unique() %>%
      unlist(., use.names = FALSE)
    By <- By[By %in% NonDom]
    By <- By %>% paste(., collapse=", ") %>% sort()
    if (nchar(By)>0)
      DomList2[[i]] <- data.frame(MP=DomMPs[i], By=By, stringsAsFactors = FALSE)
  }

  DomDF <- do.call("rbind", DomList2) %>% as.data.frame() %>% dplyr::arrange(MP)
  NonDom <- NonDom %>% sort()

  list(MPs=NonDom, DomMPs=DomDF)
}


#' @describeIn checkMSE Adds additional MPs to an MSE object by combining
#' multiple MSE objects that have identical historical OM values but different
#' MPs.
#' @export
addMPs <- function(MSEobjs) {
  # join two or more MSE objects
  if (!inherits(MSEobjs, "list")) stop("MSEobjs must be a list")
  if (length(MSEobjs) < 2) stop("MSEobjs list doesn't contain multiple MSE objects")
  if (!all(sapply(MSEobjs, inherits, "MSE"))) stop('MSEobjs must be a list of objects of class `MSE`', call.=FALSE)

  # Check seed and hist values
  for (x in 2:length(MSEobjs)) {
    check <- all(MSEobjs[[x]]@OM == MSEobjs[[1]]@OM)
    if (!check) stop('Values for `MSEobjs[[', x, ']]@OM` are not the same as `MSEobjs[[1]]@OM`')
    
    check <- all(MSEobjs[[x]]@Obs == MSEobjs[[1]]@Obs)
    if (!check) stop('Values for `MSEobjs[[', x, ']]@Obs` are not the same as `MSEobjs[[1]]@Obs`')
  }
  
  slots_identical <- function(slotname, MSEobjs, is_logical = FALSE) {
    templist <- lapply(MSEobjs, slot, slotname)
    is_identical <- vapply(templist[-1], identical, logical(1), templist[[1]]) %>% all()
    if (is_logical) {
      return(is_identical)
    } else {
      return(templist %>% unlist() %>% unique())
    }
  }
  
  if(!slots_identical("nsim", MSEobjs, TRUE)) stop("nsim slot not identical in all MSE objects.")
  if(!slots_identical("nyears", MSEobjs, TRUE)) stop("nyears slot not identical in all MSE objects.")
  if(!slots_identical("proyears", MSEobjs, TRUE)) stop("proyears slot not identical in all MSE objects.")
  
  MSE <- new("MSE",
             Name = slots_identical("Name", MSEobjs),
             nyears = slots_identical("nyears", MSEobjs),
             proyears = slots_identical("proyears", MSEobjs),
             nMPs = sapply(MSEobjs, slot, "nMPs") %>% sum(),
             MPs = lapply(MSEobjs, slot, "MPs") %>% unlist(),
             nsim = slots_identical("nsim", MSEobjs),
             OM = MSEobjs[[1]]@OM,
             Obs = MSEobjs[[1]]@Obs,
             SB_SBMSY = join_arrays(MSEobjs, "SB_SBMSY", along = 2),
             F_FMSY = join_arrays(MSEobjs, "F_FMSY", along = 2),
             N = join_arrays(MSEobjs, "N", along = 2),
             B = join_arrays(MSEobjs, "B", along = 2),
             SSB = join_arrays(MSEobjs, "SSB", along = 2),
             VB = join_arrays(MSEobjs, "VB", along = 2),
             FM = join_arrays(MSEobjs, "FM", along = 2),
             SPR = lapply(MSEobjs, slot, "SPR") %>% join_list_of_arrays(along = 2),
             Catch = join_arrays(MSEobjs, "Catch", along = 2),
             Removals = join_arrays(MSEobjs, "Removals", along = 2),
             Effort = join_arrays(MSEobjs, "Effort", along = 2),
             TAC = join_arrays(MSEobjs, "TAC", along = 2),
             TAE = join_arrays(MSEobjs, "TAE", along = 2),
             BioEco = lapply(MSEobjs, slot, "BioEco") %>% join_list_of_arrays(along = 2),
             RefPoint = list(),
             CB_hist = MSEobjs[[1]]@CB_hist,
             FM_hist = MSEobjs[[1]]@FM_hist,
             SSB_hist = MSEobjs[[1]]@SSB_hist,
             Hist = MSEobjs[[1]]@Hist,
             PPD = lapply(MSEobjs, slot, "PPD") %>% unlist(),
             Misc = list()
  )
  
  MSE@RefPoint <- local({
    Refout <- lapply(MSEobjs, function(x) x@RefPoint[c("MSY", "FMSY", "SSBMSY", "F_SPR")]) %>% join_list_of_arrays(along = 2)
    Refout$Dynamic_Unfished <- MSEobjs[[1]]@RefPoint$Dynamic_Unfished
    Refout$ByYear <- MSEobjs[[1]]@RefPoint$ByYear
    Refout
  })
  
  if (length(MSEobjs[[1]]@Misc$extended)) {
    MSE@Misc$extended <- lapply(MSEobjs, function(x) x@Misc$extended) %>% join_list_of_arrays(along = 3)
  }
  
  attr(MSE, "version") <- packageVersion("MSEtool")
  attr(MSE, "date") <- date()
  attr(MSE, "R.version") <- R.version
  
  return(MSE)
  
}

#' @describeIn checkMSE Joins two or more MSE objects together. MSE objects must have identical
#' number of historical years, and projection years.
#' @param MSEobjs A list of MSE objects
#' @export
joinMSE <- function(MSEobjs = NULL) {
  # join two or more MSE objects
  if (!methods::is(MSEobjs, "list")) stop("MSEobjs must be a list")
  if (length(MSEobjs) < 2) stop("MSEobjs list doesn't contain multiple MSE objects")

  lapply(MSEobjs, checkMSE) # check that MSE objects contains all slots

  if (!all(unlist(lapply(MSEobjs, class))=='MSE'))
    stop('Objects in list must all be class `MSE`')

  MPNames <- lapply(MSEobjs, getElement, name = "MPs")  # MPs in each object
  allsame <- length(unique(lapply(MPNames, unique))) == 1

  if (!allsame) {
    # drop the MPs that don't appear in all MSEobjs
    mpnames <- unlist(MPNames)
    npack <- length(MSEobjs)
    tab <- table(mpnames)
    ind <- tab == npack
    commonMPs <- names(tab)[ind]
    if (length(commonMPs)<1) stop("No common MPs in MSE objects", call.=FALSE)
    MSEobjs <- lapply(MSEobjs, Sub, MPs = commonMPs)
    message("MPs not in all MSE objects:")
    message(paste(names(tab)[!ind], ""))
    message("Dropped from final MSE object.")
  }

  Nobjs <- length(MSEobjs)
  for (X in 2:Nobjs) {
    if (!all(slotNames(MSEobjs[[X]]) == slotNames(MSEobjs[[X - 1]])))
      stop("The MSE objects don't have the same slots")
    if (any(MSEobjs[[X]]@MPs != MSEobjs[[X - 1]]@MPs))
      stop("MPs must be the same for all MSE objects")
  }

  # Check that nyears and proyears are the same for all
  chkmat <- matrix(NA, nrow = Nobjs, ncol = 2)
  for (X in 1:Nobjs) {
    chkmat[X, ] <- c(MSEobjs[[X]]@nyears, MSEobjs[[X]]@proyears)
  }
  
  chk <- all(colSums(chkmat) == chkmat[1, ] * Nobjs)
  if (!chk) stop("The MSE objects have different number of nyears or proyears")
  
  # Check that MSE names are identical
  nms <- length(unique(sapply(MSEobjs, slot, "Name"))) == 1
  if (!nms) stop("MSE objects have different names")
  
  # Join them together
  MSE <- new("MSE",
             Name = MSEobjs[[1]]@Name,
             nyears = MSEobjs[[1]]@nyears,
             proyears = MSEobjs[[1]]@proyears,
             nMPs = length(MPNames[[1]]),
             MPs = MPNames[[1]],
             nsim = sapply(MSEobjs, slot, "nsim") %>% sum(),
             OM = join_rows(MSEobjs, "OM"),
             Obs = join_rows(MSEobjs, "Obs"),
             SB_SBMSY = join_arrays(MSEobjs, "SB_SBMSY"),
             F_FMSY = join_arrays(MSEobjs, "F_FMSY"),
             N = join_arrays(MSEobjs, "N"),
             B = join_arrays(MSEobjs, "B"),
             SSB = join_arrays(MSEobjs, "SSB"),
             VB = join_arrays(MSEobjs, "VB"),
             FM = join_arrays(MSEobjs, "FM"),
             SPR = lapply(MSEobjs, slot, "SPR") %>% join_list_of_arrays(),
             Catch = join_arrays(MSEobjs, "Catch"),
             Removals = join_arrays(MSEobjs, "Removals"),
             Effort = join_arrays(MSEobjs, "Effort"),
             TAC = join_arrays(MSEobjs, "TAC"),
             TAE = join_arrays(MSEobjs, "TAE"),
             BioEco = lapply(MSEobjs, slot, "BioEco") %>% join_list_of_arrays(),
             RefPoint = join_MSERefPoint(MSEobjs),
             CB_hist = join_arrays(MSEobjs, "CB_hist"),
             FM_hist = join_arrays(MSEobjs, "FM_hist"),
             SSB_hist = join_arrays(MSEobjs, "SSB_hist"),
             Hist = new("Hist"),
             PPD = join_PPD(MSEobjs),
             Misc = list()
  )
  
  MSE@Hist <- lapply(MSEobjs, slot, "Hist") %>% joinHist()
  
  if (length(MSEobjs[[1]]@Misc$extended)) {
    MSE@Misc$extended <- lapply(MSEobjs, function(x) x@Misc$extended) %>% join_list_of_arrays()
  }
  
  attr(MSE, "version") <- packageVersion("MSEtool")
  attr(MSE, "date") <- date()
  attr(MSE, "R.version") <- R.version
  
  return(MSE)
}



#' @param Hist_List A list of objects of class `Hist`
#'
#' @return A new object of class `Hist`
#' @export
#' @describeIn checkMSE Join objects of class `Hist`. Does not join slot `OM`
#' @seealso \link{joinData}
joinHist <- function(Hist_List) {
  
  OMPars <- join_rows(Hist_List, "OMPars")
  
  if (nrow(OMPars)) {
    Hist <- new("Hist",
                Data = new("Data"), 
                OMPars = OMPars,
                AtAge = lapply(Hist_List, slot, "AtAge") %>% join_list_of_arrays(),
                TSdata = join_HistTSdata(Hist_List),
                Ref = join_HistRef(Hist_List),
                SampPars = join_HistSampPars(Hist_List),
                Misc = list(),
                OM = new("OM")
    )
    
    DataList <- lapply(Hist_List, slot, "Data")
    Data <- try(joinData(DataList), silent = TRUE)
    
    if (inherits(Data, "try-error")) {
      warning("joinHist() could not join Data objects.")
    } else {
      Hist@Data <- Data
    }
    
    Hist@Misc$BioEco <- lapply(Hist_List, function(x) x@Misc$BioEco) %>% bind_rows()
    
  } else {
    Hist <- new("Hist")
  }
  
  attr(Hist, "version") <- packageVersion("MSEtool")
  attr(Hist, "date") <- date()
  attr(Hist, "R.version") <- R.version
  
  return(Hist)
}