R/MSE_functions.R

Defines functions Dom Sub checkMSE Chk Converge

Documented in checkMSE Converge Dom Sub

#' 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(class(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
#' }
#' @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]]
          }

        }
    }

  }

  subMSElist <- MSEobj@PPD[SubMPs] # doesn't subset Data by years or simulations

  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=MSEobj@N[SubIts, SubMPs, Years, drop=FALSE],
                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=list(),
                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,  MSEobj@nyears, drop=FALSE],
                FM_hist=MSEobj@FM_hist[SubIts,  MSEobj@nyears, drop=FALSE],
                SSB_hist=MSEobj@SSB_hist[SubIts,  MSEobj@nyears, drop=FALSE],
                Hist=MSEobj@Hist,
                PPD=subMSElist,
                Misc=MSEobj@Misc)

  subMSE
}

# NOT YET WORKING FOR NEW VERSION OF MSETOOL

# #' @describeIn checkMSE Joins two or more MSE objects together. MSE objects must have identical
# #' number of historical years, and projection years. Also works for Hist objects returned
# #' by `runMSE(Hist=TRUE)`
# #' @export
# joinMSE <- function(MSEobjs = NULL) {
#   # join two or more MSE objects
#   if (class(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
#
#   ishist <- all(lapply(MSEobjs, slotNames) %>% unlist() %>% unique() %in% slotNames('Hist'))
#
#   if (ishist) {
#     out <- new("Hist")
#     sls <- slotNames('Hist')
#     nsim <- MSEobjs[[1]]@Ref$B0 %>% length()
#     for (sl in sls) {
#       obj <-lapply(MSEobjs, slot, name=sl)
#       if (sl == "Data") {
#         out@Data <- joinData(obj)
#       } else {
#         if (class(obj[[1]]) == "data.frame") {
#           slot(out, sl) <- do.call('rbind', obj)
#         }
#         if (class(obj[[1]]) == "list") {
#           out.list <- list()
#           for (nm in names(obj[[1]])) {
#             obj2 <- lapply(obj, '[[', nm)
#             ind <- which(dim(obj2[[1]]) == nsim)
#             if (length(ind)>1) ind <- ind[1]
#             if (length(ind) >0) {
#               if (class(obj2[[1]]) == "array") {
#                 tempVal <- lapply(obj2, dim)
#                 # check all dimensions the same (hack for different CAL bins)
#                 tdf <- lapply(obj2, dim) %>% unlist() %>%
#                   matrix(nrow=length(obj2), ncol=length(dim(obj2[[1]])),byrow=TRUE)
#                 nBins <- tdf[,2]
#                 Max <- max(nBins)
#                 nyrs <- max(tdf[,3])
#                 nsims <- sapply(tempVal, function(x) x[1])
#                 if (!mean(nBins) == max(nBins)) { # not all same size
#                   index <- which(nBins < Max)
#                   for (kk in index) {
#                     dif <- Max - dim(obj2[[kk]])[2]
#                     obj2[[kk]] <- abind::abind(obj2[[kk]], array(0, dim=c(nsims[kk], dif, nyrs)), along=2)
#                   }
#                 }
#               }
#               out.list[[nm]] <- abind::abind(obj2, along=ind)
#             } else {
#               out.list[[nm]] <- unlist(obj2) #  %>% unique()
#             }
#           }
#           slot(out, sl) <- out.list
#         }
#       }
#     }
#     return(out)
#   }
#
#   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 1:Nobjs) {
#     tt <- MSEobjs[[X]]
#     assign(paste0("obj", X), tt)
#     if (X > 1) {
#       tt <- MSEobjs[[X]]
#       tt2 <- MSEobjs[[X - 1]]
#       if (!all(slotNames(tt) == slotNames(tt2)))
#         stop("The MSE objects don't have the same slots")
#       if (any(tt@MPs != tt2@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)
#   nms <- NULL
#   for (X in 1:Nobjs) {
#     tt <- get(paste0("obj", X))
#     chkmat[X, ] <- c(tt@nyears, tt@proyears)
#     if (X > 1)
#       if (!any(grepl(tt@Name, nms)))
#         stop("MSE objects have different names")
#     nms <- append(nms, tt@Name)
#   }
#   chk <- all(colSums(chkmat) == chkmat[1, ] * Nobjs)
#   if (!chk) stop("The MSE objects have different number of nyears or proyears")
#
#   # Join them together
#   Allobjs <- mget(paste0("obj", 1:Nobjs))
#   sns <- slotNames(Allobjs[[1]])
#   sns<-sns[sns!="Misc"] # ignore the Misc slot
#   outlist <- vector("list", length(sns))
#   for (sn in 1:length(sns)) {
#     templs <- lapply(Allobjs, slot, name = sns[sn])
#     if (class(templs[[1]]) == "character") {
#       outlist[[sn]] <- templs[[1]]
#     }
#     if (class(templs[[1]]) == "numeric" | class(templs[[1]]) == "integer") {
#       if (sns[sn] == "CALbins") {
#         tempInd <- which.max(unlist(lapply(templs, length)))
#         CALbins <- templs[[tempInd]]
#       } else {
#         outlist[[sn]] <- do.call(c, templs)
#       }
#     }
#     if (class(templs[[1]]) == "matrix" | class(templs[[1]]) == "data.frame") {
#       outlist[[sn]] <- do.call(rbind, templs)
#     }
#     if (class(templs[[1]]) == "array") {
#       if (sns[sn] == "CAL") { # hack for different sized CAL arrays
#         tempVal <- lapply(templs, dim)
#         if (all(unlist(lapply(tempVal, length)) == 3)) {
#           nBins <- sapply(tempVal, function(x) x[3])
#           nsims <- sapply(tempVal, function(x) x[1])
#           nMPs <- sapply(tempVal, function(x) x[2])
#           if (!mean(nBins) == max(nBins)) { # not all same size
#             Max <- max(nBins)
#             index <- which(nBins < Max)
#             for (kk in index) {
#               dif <- Max - dim(templs[[kk]])[3]
#               templs[[kk]] <- abind::abind(templs[[kk]], array(0, dim=c(nsims[kk], nMPs[kk], dif)), along=3)
#             }
#           }
#           outlist[[sn]] <- abind::abind(templs, along = 1)
#         } else {
#           outlist[[sn]] <- templs[[1]]
#         }
#       } else {
#         outlist[[sn]] <- abind::abind(templs, along = 1)
#       }
#
#     }
#   }
#
#   names(outlist) <- sns
#
#   Misc<-list()
#   if (length(MSEobjs[[1]]@Misc)>0) {
#     if (!is.null(MSEobjs[[1]]@Misc$Data)) {
#       Misc$Data <- list()
#       # Posterior predicted data joining
#       for(i in 1:length(MSEobjs[[1]]@Misc$Data))
#         Misc$Data[[i]]<-joinData(lapply(MSEobjs,function(x)slot(x,"Misc")$Data[[i]]))
#     }
#
#     if (!is.null(MSEobjs[[1]]@Misc$RInd.stats)) {
#       # Error from real indices
#       nms <- unique(MSEobjs[[1]]@Misc$RInd.stats$Index) %>% as.character()
#       temp <- list()
#       for (nm in seq_along(nms)) {
#         temp1 <- list()
#         for(i in 1:length(MSEobjs)) {
#           temp1[[i]] <- MSEobjs[[i]]@Misc$RInd.stats %>% dplyr::filter(Index==nms[nm])
#         }
#         temp[[nm]] <- do.call('rbind', temp1)
#       }
#       Misc$RInd.stats <- do.call('rbind', temp)
#     }
#
#     if (!is.null(MSEobjs[[1]]@Misc$TryMP)) {
#       temp1 <- list()
#       for(i in 1:length(MSEobjs)) {
#         temp1[[i]] <- MSEobjs[[i]]@Misc$TryMP
#       }
#       Misc$TryMP <- do.call('rbind', temp1)
#     }
#
#     if (!is.null(MSEobjs[[1]]@Misc$Unfished)) {
#       Misc$Unfished <- list()
#       temp1 <- temp2 <- list()
#       for(i in 1:length(MSEobjs)) {
#         temp1[[i]] <- MSEobjs[[i]]@Misc$Unfished$Refs
#         temp2[[i]] <- MSEobjs[[i]]@Misc$Unfished$ByYear
#       }
#       Misc$Unfished$Refs <- do.call('cbind', temp1)
#       for (nm in names(temp2[[1]])) {
#         tt = lapply(temp2, "[[", nm)
#         tt <- do.call('rbind',tt)
#         Misc$Unfished$ByYear[[nm]] <- tt
#       }
#     }
#     if (!is.null(MSEobjs[[1]]@Misc$MSYRefs)) {
#       Misc$MSYRefs <- list()
#       temp1 <- temp2 <- list()
#       for(i in 1:length(MSEobjs)) {
#         temp1[[i]] <- MSEobjs[[i]]@Misc$MSYRefs$Refs
#         temp2[[i]] <- MSEobjs[[i]]@Misc$MSYRefs$ByYear
#       }
#       Misc$MSYRefs$Refs <- do.call('rbind', temp1)
#       for (nm in names(temp2[[1]])) {
#         tt <- lapply(temp2, "[[", nm)
#         if (length(dim(tt[[1]])) == 3) {
#           tt <- abind::abind(tt, along=1)
#         } else {
#           tt <- do.call('rbind',tt)
#         }
#         Misc$MSYRefs$ByYear[[nm]] <- tt
#       }
#     }
#     temp <- list()
#     nsim <- ncol(Misc$Unfished$Ref)
#     dims <- dim(MSEobjs[[1]]@Misc$LatEffort)
#     Misc$LatEffort <- array(NA, dim=c(nsim, dims[2], dims[3]))
#     Misc$Revenue <- array(NA, dim=c(nsim, dims[2], dims[3]))
#     Misc$Cost <- array(NA, dim=c(nsim, dims[2], dims[3]))
#     Misc$TAE <- array(NA, dim=c(nsim, dims[2], dims[3]))
#     st <- 1
#     for (i in 1:length(MSEobjs)) {
#       dims <- dim(MSEobjs[[i]]@Misc$LatEffort)
#       indvec <- st:(st+dims[1]-1)
#       st <- indvec[length(indvec)] + 1
#       Misc$LatEffort[indvec,,] <- MSEobjs[[i]]@Misc$LatEffort
#       Misc$Cost[indvec,,] <- MSEobjs[[i]]@Misc$Cost
#       Misc$Revenue[indvec,,] <- MSEobjs[[i]]@Misc$Revenue
#       Misc$TAE[indvec,,] <- MSEobjs[[i]]@Misc$TAE
#     }
#
#   }
#
#   # ErrList
#   nms <- names(MSEobjs[[1]]@Misc$ErrList)
#   for (nm in nms) {
#     t1 <- lapply(1:length(MSEobjs), function(i) MSEobjs[[i]]@Misc$ErrList[[nm]])
#     Misc$ErrList[[nm]] <- do.call('rbind', t1)
#   }
#
#   # Removals
#   t1 <- lapply(1:length(MSEobjs), function(i) MSEobjs[[i]]@Misc$Removals)
#   Misc$Removals <- abind::abind(t1, along=1)
#
#
#   newMSE <- new("MSE", Name = outlist$Name, nyears = unique(outlist$nyears),
#                 proyears = unique(outlist$proyears), nMP = unique(outlist$nMP),
#                 MPs = unique(outlist$MPs), nsim = sum(outlist$nsim), OM = outlist$OM,
#                 Obs = outlist$Obs, SB_SBMSY = outlist$B_BMSY, F_FMSY = outlist$F_FMSY,
#                 outlist$B, outlist$SSB, outlist$VB,
#                 outlist$FM, outlist$C, outlist$TAC, outlist$SSB_hist,
#                 outlist$CB_hist, outlist$FM_hist, outlist$Effort, outlist$PAA,
#                 outlist$CAA, outlist$CAL, CALbins, Misc=Misc)
#
#   newMSE
# }
#
# 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 (class(MSEobj) != 'MSE') stop("Object must be class `MSE`", call.=FALSE)
  if (is.null(PMlist)) {
    PMlist <- unlist(list(...))
  } else {
    PMlist <- unlist(PMlist)
  }

  if (class(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)
}

# NOT WORKING FOR NEW VERSION

# #' @describeIn checkMSE Adds additional MPs to an MSE object by combining
# #'  multiple MSE objects that have identical historical OM values but different
# #'  MPs. Note that the `Misc` slot is returned as a list of `length(MSEobjs)` with
# #'  each element containing the `Misc` object from `MSEobjs`.
# #' @export
# #'
# addMPs <- function(MSEobjs) {
#   # join two or more MSE objects
#   if (class(MSEobjs) != "list") stop("MSEobjs must be a list")
#   if (length(MSEobjs) < 2) stop("MSEobjs list doesn't contain multiple MSE objects")
#   if (!(all(lapply(MSEobjs, class) %>% unlist() =='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`')
#   }
#
#   MSEout <- MSEobjs[[1]]
#   MPs <- lapply(MSEobjs, slot, 'MPs') %>% unlist()
#   MSEout@MPs <- MPs
#   MSEout@nMPs <- length(MPs)
#
#   addmp <- function(sl='SB_SBMSY', MSEobjs, MSEout) {
#     vals <- lapply(MSEobjs, slot, sl)
#     slot(MSEout, sl) <- abind::abind(vals, along=2)
#     MSEout
#   }
#   MSEout <- addmp('SB_SBMSY', MSEobjs, MSEout)
#   MSEout <- addmp('F_FMSY', MSEobjs, MSEout)
#   MSEout <- addmp('B', MSEobjs, MSEout)
#   MSEout <- addmp('SSB', MSEobjs, MSEout)
#   MSEout <- addmp('VB', MSEobjs, MSEout)
#   MSEout <- addmp('FM', MSEobjs, MSEout)
#   MSEout <- addmp('C', MSEobjs, MSEout)
#   MSEout <- addmp('TAC', MSEobjs, MSEout)
#   MSEout <- addmp('Effort', MSEobjs, MSEout)
#   MSEout <- addmp('PAA', MSEobjs, MSEout)
#   MSEout <- addmp('CAA', MSEobjs, MSEout)
#   MSEout <- addmp('CAL', MSEobjs, MSEout)
#
#   MSEout@Misc <- lapply(MSEobjs, slot, 'Misc')
#
#   MSEout
# }
#
#

#' @describeIn checkMSE Joins two or more MSE objects together. MSE objects must have identical
#' number of historical years, and projection years. Also works for Hist objects returned
#' by `runMSE(Hist=TRUE)`
#' @param MSEobjs A list of MSE objects
#' @export
joinMSE <- function(MSEobjs = NULL) {
  # join two or more MSE objects
  if (class(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 1:Nobjs) {
    tt <- MSEobjs[[X]]
    assign(paste0("obj", X), tt)
    if (X > 1) {
      tt <- MSEobjs[[X]]
      tt2 <- MSEobjs[[X - 1]]
      if (!all(slotNames(tt) == slotNames(tt2)))
        stop("The MSE objects don't have the same slots")
      if (any(tt@MPs != tt2@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)
  nms <- NULL
  for (X in 1:Nobjs) {
    tt <- get(paste0("obj", X))
    chkmat[X, ] <- c(tt@nyears, tt@proyears)
    if (X > 1)
      if (!any(grepl(tt@Name, nms)))
        stop("MSE objects have different names")
    nms <- append(nms, tt@Name)
  }
  chk <- all(colSums(chkmat) == chkmat[1, ] * Nobjs)
  if (!chk) stop("The MSE objects have different number of nyears or proyears")

  # Join them together
  Allobjs <- mget(paste0("obj", 1:Nobjs))

  Name <- Allobjs[[1]]@Name
  nyears <- Allobjs[[1]]@nyears
  proyears <- Allobjs[[1]]@proyears
  nMPs <- length(MPNames[[1]])
  MPs <- MPNames[[1]]
  nsim <- sum(sapply(Allobjs, slot, name = 'nsim'))

  OM <- do.call('rbind', lapply(Allobjs, slot, name = 'OM'))
  Obs <- do.call('rbind', lapply(Allobjs, slot, name = 'Obs'))


  SB_SBMSY <- abind::abind(lapply(Allobjs, slot, name = 'SB_SBMSY'), along = 1)
  F_FMSY <- abind::abind(lapply(Allobjs, slot, name = 'F_FMSY'), along = 1)
  N <- abind::abind(lapply(Allobjs, slot, name = 'N'), along = 1)
  B <- abind::abind(lapply(Allobjs, slot, name = 'B'), along = 1)
  SSB <- abind::abind(lapply(Allobjs, slot, name = 'SSB'), along = 1)
  VB <- abind::abind(lapply(Allobjs, slot, name = 'VB'), along = 1)
  FM <- abind::abind(lapply(Allobjs, slot, name = 'FM'), along = 1)
  SPR <- list() # not currently returned in MSE
  Catch <- abind::abind(lapply(Allobjs, slot, name = 'Catch'), along = 1)
  Removals <- abind::abind(lapply(Allobjs, slot, name = 'Removals'), along = 1)
  Effort <- abind::abind(lapply(Allobjs, slot, name = 'Effort'), along = 1)
  TAC <- abind::abind(lapply(Allobjs, slot, name = 'TAC'), along = 1)
  TAE <- abind::abind(lapply(Allobjs, slot, name = 'TAE'), along = 1)

  BioEco <- list()
  BioEco_List <- lapply(Allobjs, slot, name = 'BioEco')
  nms <- names(Allobjs[[1]]@BioEco)
  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Allobjs)) {
      temp[[obj]] <-BioEco_List[[obj]][[nm]]
    }

    BioEco[[nm]] <- abind::abind(temp, along = 1)
  }

  RefPoint <- list()
  RefPoint_List <- lapply(Allobjs, slot, name = 'RefPoint')
  nms <- names(Allobjs[[1]]@RefPoint)
  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Allobjs)) {
      temp[[obj]] <-RefPoint_List[[obj]][[nm]]
    }
    if (class(temp[[1]]) == 'list') {
      nms2 <- names(temp[[1]])
      for (j in seq_along(nms2)) {
        temp2 <- list()
        for (k in 1:length(temp)) {
          temp2[[k]] <- temp[[k]][[j]]
        }
        nm2 <- nms2[j]
        RefPoint[[nm]][[nm2]] <- abind::abind(temp2, along = 1)
        }
    } else {
      RefPoint[[nm]] <- abind::abind(temp, along = 1)
    }
  }

  CB_hist <- abind::abind(lapply(Allobjs, slot, name = 'CB_hist'), along = 1)
  FM_hist <- abind::abind(lapply(Allobjs, slot, name = 'FM_hist'), along = 1)
  SSB_hist <- abind::abind(lapply(Allobjs, slot, name = 'SSB_hist'), along = 1)

  Hist_List <- lapply(Allobjs, slot, name = 'Hist')
  Hist <- joinHist(Hist_List)

  # PPD
  PPD <- list()
  for (m in 1:nMPs) {
    temp <- list()
    for (obj in 1:length(Allobjs)) {
      temp[[obj]] <- Allobjs[[obj]]@PPD[[m]]
    }
    if (length(temp)>0)
      PPD[[m]] <- joinData(temp)
  }


  # Misc
  # TODO
  Misc <- list()

  MSE <- new('MSE', Name, nyears, proyears, nMPs, MPs, nsim, OM, Obs,
             SB_SBMSY, F_FMSY, N, B, SSB, VB, FM, SPR, Catch, Removals, Effort,
             TAC, TAE, BioEco, RefPoint,
             CB_hist, FM_hist, SSB_hist,
             Hist, PPD, Misc)
  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`
joinHist <- function(Hist_List) {
  # TODO
  # DataList <- lapply(Hist_List, slot, name = 'Data')
  # Data <- joinData(DataList)
  Data <- new('Data')
  OMPars <- abind::abind(lapply(Hist_List, slot, name = 'OMPars'), along = 1)
  OMPars <- data.frame(OMPars)

  AtAge <- list()
  AtAge_List <- lapply(Hist_List, slot, name = 'AtAge')
  nms <- names(Hist_List[[1]]@AtAge)
  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Hist_List)) {
      temp[[obj]] <-AtAge_List[[obj]][[nm]]
    }

    AtAge[[nm]] <- abind::abind(temp, along = 1)
  }

  TSdata <- list()
  TSdata_List <- lapply(Hist_List, slot, name = 'TSdata')
  nms <- names(Hist_List[[1]]@TSdata)
  nms <- nms[nms!='Unfished_Equilibrium']

  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Hist_List)) {
      temp[[obj]] <- TSdata_List[[obj]][[nm]]
    }
    TSdata[[nm]] <- abind::abind(temp, along = 1)
  }

  TSdata$Unfished_Equilibrium <- list()
  nms <- names(Hist_List[[1]]@TSdata$Unfished_Equilibrium)
  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Hist_List)) {
      temp[[obj]] <- TSdata_List[[obj]]$Unfished_Equilibrium[[nm]]
    }
    TSdata$Unfished_Equilibrium[[nm]] <- abind::abind(temp, along = 1)
  }

  Ref <- list()
  Ref_List <- lapply(Hist_List, slot, name = 'Ref')
  Ref$ByYear <- list()
  nms <- names(Hist_List[[1]]@Ref$ByYear)[1:2]
  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Hist_List)) {
      temp[[obj]] <- Ref_List[[obj]]$ByYear[[nm]]
    }
    Ref$ByYear[[nm]] <- abind::abind(temp, along = 1)
  }

  Ref$Dynamic_Unfished <- list()
  nms <- names(Hist_List[[1]]@Ref$Dynamic_Unfished)[1:2]
  for (nm in nms) {
    temp <- list()
    for (obj in 1:length(Hist_List)) {
      temp[[obj]] <- Ref_List[[obj]]$Dynamic_Unfished[[nm]]
    }
    Ref$Dynamic_Unfished[[nm]] <- abind::abind(temp, along = 1)
  }

  temp <- list()
  for (obj in 1:length(Hist_List)) {
    temp[[obj]] <- Ref_List[[obj]]$ReferencePoints
  }
  Ref$ReferencePoints <- abind::abind(temp, along = 1)


  SampPars <- list()
  # TO DO

  Misc <- list()
  OM <- new("OM")

  sls <- slotNames("Hist")
  Hist <- new("Hist")
  for (sl in sls) {
    slot(Hist, sl) <- get(sl)
  }
  Hist
}

Try the MSEtool package in your browser

Any scripts or data that you put into this service are public.

MSEtool documentation built on July 8, 2021, 5:09 p.m.