R/Misc_Exported.R

In MSEtool: Management Strategy Evaluation Toolkit

get_funcs <- function(package, classy , msg) {
  pkgs <- search()
  search_package <- paste0("package:",package)
  funs <- NULL

  if (search_package %in% pkgs) {
    if (msg)
      message('Searching for objects of class ', classy, ' in package: ', package)
    funs <- ls(search_package)[vapply(ls(search_package),
                               getclass,
                               logical(1),
                               classy = classy)]
  } else {
    stop('Package ', package, ' not loaded. Use `library(', package, ')`', call. = FALSE)
  }
  funs
}

#' What objects of this class are available
#'
#' Generic class finder
#'
#' Finds objects of the specified class in the global environment or the
#' openMSE packages.
#'
#' @param classy A class of object (character string, e.g. 'Fleet')
#' @param package Optional. Names(s) of the package to search for object of class `classy`. String
#' Default is all `openMSE` packages. Always searches the global environment as well.
#' @param msg Print messages?
#' @examples
#' avail("OM", msg=FALSE)
#' @author T. Carruthers
#' @seealso \link{Can} \link{Cant} \link{avail}
#' @examples
#' Stocks <- avail("Stock")
#' Fleets <- avail("Fleet")
#' MPs <- avail("MP")
#' @export
avail <- function(classy, package=NULL, msg=TRUE) {
  temp <- try(class(classy), silent=TRUE)
  if (methods::is(temp, "try-error")) classy <- deparse(substitute(classy))
  if (temp == "function") classy <- deparse(substitute(classy))

  if (classy %in% c('Output', 'Input', "Mixed", "Reference")) {
    MPs <- avail('MP', msg=msg)
    gettype <- MPtype(MPs)
    temp <- gettype[gettype[,2] %in% classy,1]
    if (length(temp) < 1) stop("No MPs of type '", classy, "' found", call. = FALSE)
    return(temp)

  } else {

    packages <- c('MSEtool', 'SAMtool', 'DLMtool', 'MSEextra')
    if (is.null(package)) {
      package <- packages
      pkgs <- search()
      search_package <- paste0("package:",package)
      package <- package[search_package %in% pkgs]
    }

    global_funs <- ls(envir = .GlobalEnv)[vapply(ls(envir = .GlobalEnv), getclass, logical(1), classy = classy)]
    temp <- global_funs

    if ('MSEtool' %in% package) {
      MSEtool_funs <- get_funcs('MSEtool', classy, msg)
      temp <- c(temp, MSEtool_funs)
    }
    if ('SAMtool' %in% package) {
      SAMtool_funs <- get_funcs('SAMtool', classy, msg)
      temp <- c(temp, SAMtool_funs)
    }
    if ('DLMtool' %in% package) {
      DLMtool_funs <- get_funcs('DLMtool', classy, msg)
      temp <- c(temp, DLMtool_funs)
    }
    if ('MSEextra' %in% package) {
      MSEextra_funs <- get_funcs('MSEextra', classy, msg)
      temp <- c(temp, MSEextra_funs)
    }
    
    packagex <- package[!package %in% packages]
    if (length(packagex)>0) {
      other <- sapply(1:length(packagex), function(i)
        get_funcs(packagex[i], classy, msg))  
      other <- unlist(other)
      temp <- c(temp, other)
    }
    
    if (length(temp) < 1) stop("No objects of class '", classy, "' found", call. = FALSE)
    return(unique(temp))
  }
}


#' Directory of the data in the installed package on your computer
#'
#' A way of locating where the package was installed so you can find example
#' data files and code etc.
#'
#' @param stock Character string representing the name of a .csv file e.g.
#' 'Snapper', 'Rockfish'
#' @author T. Carruthers
#' @return The file path to the object
#' @examples
#' \dontrun{
#' tilefish_location <- DataDir("Gulf_blue_tilefish")
#' tilefish_Data <- new("Data", tilefish_location)
#' }
#' @export DataDir
DataDir <- function(stock = NA) {
  if (is.na(stock)) {
    file.path(system.file(package = "MSEtool"), "Data")
  } else {
    system.file(paste0('Data/', stock, ".csv"), package = "MSEtool", mustWork = TRUE)
  }
}


#' Directory of the installed package on your computer
#'
#' @param stock Character string representing the name of a .csv file e.g.
#' 'Snapper', 'Rockfish'
#'
#' @return The file path to the object
#' @export
#'
DLMDataDir <- function(stock = NA) {
  .Deprecated('DataDir')
  DataDir(stock)
}


#' Load more data from MSEextra package
#'
#' Downloads the MSEextra package from GitHub
#' @param silent Logical. Should messages to printed?
#' @param force Logical. For install from github if package is up-to-date?
#' @export
#'
MSEextra <- function(silent=FALSE, force=FALSE) {
  if (!requireNamespace("remotes", quietly = TRUE)) {
    stop("remotes is needed for this function to work. Please install it.",
         call. = FALSE)
  }

  if (!silent) message("\nDownloading 'MSEextra' from GitHub")
  remotes::install_github("blue-matter/MSEextra", quiet=FALSE, force=force)
  if (!silent) message("Use 'library(MSEextra)' to load additional data into workspace")

}



#' @rdname SubCpars
#' @export
setGeneric("SubCpars", function(x, ...) standardGeneric("SubCpars"))

#' @name SubCpars
#' @aliases SubCpars,OM-method
#' @title Subset the cpars slot in an operating model
#'
#' @description Subset the custom parameters of an operating model by simulation and projection years 
#'
#' @param x An object of class \linkS4class{OM} or \linkS4class{MOM}
#' @param sims A logical vector of length \code{x@@nsim} to either retain (TRUE) or remove (FALSE).
#' Alternatively, a numeric vector indicating which simulations (from 1 to nsim) to keep.
#' @param proyears If provided, a numeric to reduce the number of projection years (must be less than \code{x@@proyears}).
#' @param silent Logical to indicate if messages will be reported to console.
#' @param ... Arguments for method.
#' @details Useful function for running \link{multiMSE} in batches if running into memory constraints.
#' @return An object of class \linkS4class{OM} or \linkS4class{MOM} (same class as \code{x}).
#' @seealso \link{Sub} for MSE objects, \link{SubOM} for OM components.
#' @author T. Carruthers, Q. Huynh
#' @export
setMethod("SubCpars", signature(x = "OM"),
          function(x, sims = 1:x@nsim, proyears = x@proyears, silent = FALSE) {
            OM <- x
            # Reduce the number of simulations
            nsim_full <- OM@nsim
            if(is.numeric(sims)) {
              sims2 <- logical(nsim_full)
              sims2[sims] <- TRUE
            } else if(is.logical(sims) && length(sims) == nsim_full) {
              sims2 <- sims
            } else stop("Logical vector sims need to be of length ", nsim_full)
            
            if(any(!sims2) && sum(sims2) < nsim_full) {
              if (!silent) message("Removing simulations: ", paste0(which(!sims2), collapse = " "))
              OM@nsim <- sum(sims2)      
              if (!silent) message("Set OM@nsim = ", OM@nsim)
              
              if(length(OM@cpars)) {
                cpars <- OM@cpars
                OM@cpars <- lapply(names(cpars), SubCpars_sim, sims = sims2, cpars = cpars) %>% 
                  structure(names = names(cpars))
              }
            }
            
            # Reduce the number of projection years
            proyears_full <- OM@proyears
            if(proyears < proyears_full) {
              if (!silent) message("Reducing the number of projection years from ", proyears_full, " to ", proyears)
              OM@proyears <- proyears
              
              if(length(OM@cpars)) {
                cpars_p <- OM@cpars
                yr_diff <- proyears_full - proyears
                OM@cpars <- lapply(names(cpars_p), SubCpars_proyears, yr_diff = yr_diff, cpars = cpars_p) %>% 
                  structure(names = names(cpars_p))
              }
            } else if(proyears > proyears_full) {
              if (!silent) message("Number of specified projection years is greater than OM@proyears. Nothing done.")
            }
            
            return(OM)
          })

#' @name SubCpars
#' @aliases SubCpars,MOM-method
#' @export
setMethod("SubCpars", signature(x = "MOM"),
          function(x, sims = 1:x@nsim, proyears = x@proyears, silent = FALSE) {
            MOM <- x
            # Reduce the number of simulations
            nsim_full <- MOM@nsim
            if(is.numeric(sims)) {
              sims2 <- logical(nsim_full)
              sims2[sims] <- TRUE
            } else if(is.logical(sims) && length(sims) == nsim_full) {
              sims2 <- sims
            } else stop("Logical vector sims need to be of length ", nsim_full)
            
            if(any(!sims2) && sum(sims2) < nsim_full) {
              if (!silent) message("Removing simulations: ", paste0(which(!sims2), collapse = " "))
              MOM@nsim <- sum(sims2)      
              if (!silent) message("Set MOM@nsim = ", MOM@nsim)
              
              if(length(MOM@cpars)) {
                for(p in 1:length(MOM@cpars)) {
                  for(f in 1:length(MOM@cpars[[p]])) {
                    cpars <- MOM@cpars[[p]][[f]]
                    MOM@cpars[[p]][[f]] <- lapply(names(cpars), SubCpars_sim, sims = sims2, cpars = cpars) %>% 
                      structure(names = names(cpars))

                  }
                }
              }
            }
            
            # Reduce the number of projection years
            proyears_full <- MOM@proyears
            if(proyears < proyears_full) {
              if (!silent) message("Reducing the number of projection years from ", proyears_full, " to ", proyears)
              MOM@proyears <- proyears
              
              if(length(MOM@cpars)) {
                yr_diff <- proyears_full - proyears
                nms <- names(MOM@cpars)
                for(p in 1:length(MOM@cpars)) {
                  if (nms[p]=='control')
                    next()
                  for(f in 1:length(MOM@cpars[[p]])) {
                    cpars_p <- MOM@cpars[[p]][[f]]
                    MOM@cpars[[p]][[f]] <- lapply(names(cpars_p), SubCpars_proyears, yr_diff = yr_diff, cpars = cpars_p) %>% 
                      structure(names = names(cpars_p))
                  
                  }
                }
              }
            } else if(proyears > proyears_full) {
              if (!silent) message("Number of specified projection years is greater than MOM@proyears. Nothing done.")
            }
            
            return(MOM)
          })

SubCpars_sim <- function(xx, sims, cpars) {
  x <- cpars[[xx]]
  vars_nosim <- c("CAL_bins", "MPA", "plusgroup", "CAL_binsmid", "binWidth", "AddIunits", "Wa", "Wb", "Data",
                  'Real.Data.Map')
  if(any(xx == vars_nosim)) {
    return(x)
  } else if(xx == "SRR") {
    x$SRRpars <- x$SRRpars[sims, ] # data.frame
    return(x)
  } else if(is.matrix(x)) {
    return(x[sims, , drop = FALSE])
  } else if(is.array(x)) {
    if(length(dim(x)) == 3) return(x[sims, , , drop = FALSE])
    if(length(dim(x)) == 4) return(x[sims, , , , drop = FALSE])
    if(length(dim(x)) == 5) return(x[sims, , , , , drop = FALSE])
  } else if(length(x) == length(sims)) {
    return(x[sims])
  } else return(x)
}


SubCpars_proyears <- function(xx, yr_diff, cpars) {
  x <- cpars[[xx]]
  vars_noproyears <-  c("Asize", "Find", "AddIbeta", "Data", "SRR", "Real.Data.Map")
  if(xx %in% vars_noproyears) { # Matrices or arrays without projection year dimensions
    return(x)
  } else if(xx == "MPA") {
    yr_remove <- (nrow(x) - yr_diff + 1):nrow(x)
    return(x[-yr_remove, ])
  } else if(is.matrix(x)) {
    yr_remove <- (ncol(x) - yr_diff + 1):ncol(x)
    return(x[, -yr_remove])
  } else if(is.array(x)) {
    
    ldim <- length(dim(x))
    yr_remove <- (dim(x)[ldim] - yr_diff + 1):dim(x)[ldim]
    
    if(ldim == 3) return(x[, , -yr_remove, drop = FALSE])
    if(ldim == 4) return(x[, , , -yr_remove, drop = FALSE])
    if(ldim == 5) return(x[, , , , -yr_remove, drop = FALSE])
  } else {
    return(x)
  }
}
          



#' @rdname tinyErr
#' @export
setGeneric("tinyErr", function(x, ...) standardGeneric("tinyErr"))

#' @name tinyErr
#' @aliases tinyErr,OM-method
#' @title Remove observation, implementation, and process error
#'
#' @description Takes an existing OM object and converts it to one without any observation
#' error, implementation error, very little process error, and/or gradients in
#' life history parameters and catchability.
#'
#' @details Useful for debugging and testing that MPs perform as expected under perfect conditions.
#'
#' @param x An object of class `OM`
#' @param ... Arguments to generic function
#' @param obs Logical. Remove observation error? `Obs` is replaced with `Perfect_Info`
#' @param imp Logical. Remove implementation error? `Imp` is replaced with `Perfect_Imp`
#' @param proc Logical. Remove process error? All `sd` and `cv` slots in `Stock`
#' and `Fleet` object are set to 0.
#' @param grad Logical. Remove gradients? All `grad` slots in `Stock` and
#' `qinc` in `Fleet` are set to 0.
#' @param silent Logical. Display messages?
#'
#' @templateVar url modifying-the-om
#' @templateVar ref the-tinyerr-function
# #' @template userguide_link
#'
#' @return An updated object of class `OM`
#' @export
#'
#' @examples
#' OM_noErr <- tinyErr(MSEtool::testOM)
setMethod("tinyErr", signature(x = "OM"),
          function(x, obs=TRUE, imp=TRUE, proc=TRUE, grad=TRUE, silent=FALSE) {
            OM <- x
            if (length(OM@cpars)>0)
              warning("Note that this function doesn't apply to parameters in cpars.\n Must be removed manually e.g `OM@cpars$Perr_y <- NULL`")

            if (!inherits(OM, 'OM')) stop("Object must be class `OM`", call.=FALSE)
            OMperf <- new("OM", MSEtool::Albacore, MSEtool::Generic_Fleet,
                          MSEtool::Perfect_Info, MSEtool::Perfect_Imp)
            OMout <- OM

            if (obs) {
              if (!silent) message("Removing all Observation Error")
              OMout <- Replace(OMout, OMperf, "Obs", silent = TRUE)
            }
            if (imp) {
              if (!silent) message("Removing all Implementation Error")
              OMout <- Replace(OMout, OMperf, "Imp", silent = TRUE)
            }
            if (proc) {
              if (!silent) message("Removing all Process Error")
              vars <- c("cv", "sd", "Perr", "AC")
              nms <- c(slotNames('Stock'), slotNames('Fleet'))
              ind <- unique(grep(paste(vars, collapse = "|"), nms, value = FALSE))
              for (X in seq_along(ind)) {
                n <- length(slot(OMout, nms[ind[X]]))
                if (n == 0) n <- 2
                slot(OMout, nms[ind[X]]) <- rep(0, n)
              }
            }
            if (grad) {
              if (!silent)  message("Removing all Gradients")
              vars <- c("grad", "inc")
              nms <- c(slotNames('Stock'), slotNames('Fleet'))
              ind <- unique(grep(paste(vars, collapse = "|"), nms, value = FALSE))
              for (X in seq_along(ind)) {
                n <- length(slot(OMout, nms[ind[X]]))
                if (n == 0) n <- 2
                slot(OMout, nms[ind[X]]) <- rep(0, n)
              }
            }
            OMout
          })



#' Management Procedure Type
#'
#' @param MPs A vector of MP names. If none are provided function is run on all available MPs
#'
#' @return A data.frame with MP names, management type (e.g "Input", "Output") and management recommendations returned by the MP
#' (e.g, TAC (total allowable catch), TAE (total allowable effort), SL (size-selectivity), and/or or Spatial)
#' @export
#' @seealso \link{Required}
#' @examples
#' MPtype(c("AvC", "curE", "matlenlim", "MRreal", "FMSYref"))
#'
MPtype <- function(MPs=NA) {
  if(methods::is(MPs, "MP")) stop("MPs must be characters")
  availMPs <- avail("MP", msg=FALSE)
  if (any(is.na(MPs))) MPs <- availMPs
  if (!methods::is(MPs, 'character')) stop("MPs must be characters")

  existMPs <- MPs %in% availMPs

  if (any(!existMPs))
    warning(paste0('Some MPs are not found in environment: ', MPs[!existMPs], collapse=", "))

  Data <- MSEtool::SimulatedData
  runMPs <- applyMP(Data, MPs[existMPs], reps = 2, nsims=1, silent=TRUE)
  recs <- runMPs[[1]]

  type <- rep("NA", length(MPs[existMPs]))
  rec <- rep("", length(MPs[existMPs]))
  rectypes <- c("TAE", "Spatial", "Selectivity", 'Retention', "Discards")
  for (mm in seq_along(recs)) {
    Effort <- Spatial <- Selectivity <- Retention <- Discards<- FALSE
    output <- length(recs[[mm]]$TAC) > 0
    names <- names(recs[[mm]])
    names <- names[!names %in% c("TAC", "Spatial", 'type')]
    input <- sum(unlist(lapply(Map(function(x) recs[[mm]][[x]], names), length))) > 0
    if (all(!is.na(recs[[mm]]$Spatial))) input <- TRUE
    if (output) {
      type[mm] <- "Output"
      thisrec <- "TAC"
    }
    if (input) {
      # what recommendations have been made?
      if (any(is.finite(recs[[mm]]$Effort))) Effort <- TRUE
      if (any(is.finite(recs[[mm]]$Spatial))) Spatial <- TRUE
      if (any(is.finite(recs[[mm]]$LR5)) | any(is.finite(recs[[mm]]$LFR)) | any(is.finite(recs[[mm]]$HS)) |
          any(is.finite(recs[[mm]]$Rmaxlen))) Retention <- TRUE
      if (any(is.finite(recs[[mm]]$L5)) | any(is.finite(recs[[mm]]$LFS)) |
          any(is.finite(recs[[mm]]$Vmaxlen))) Selectivity <- TRUE
      if (any(is.finite(recs[[mm]]$DR)) | any(is.finite(recs[[mm]]$Fdisc))) Discards <- TRUE

      dorecs <- rectypes[c(Effort, Spatial, Selectivity, Retention, Discards)]
      thisrec <- dorecs
      type[mm] <- "Input"

    }
    if (input & output) {
      type[mm] <- "Mixed"
      thisrec <- c("TAC", thisrec)
    }
    if (length(thisrec)>1)  {
      rec[mm] <- paste(thisrec, collapse=", ")
    } else {
      rec[mm] <- thisrec
    }
  }
  type[grep("ref", MPs)] <- "Reference"

  df <- data.frame(MP=MPs[existMPs], Type=type, Recs=rec, stringsAsFactors = FALSE)
  if (sum(!existMPs)>0) {
    df_non <- data.frame(MP=MPs[!existMPs], Type='unknown', Recs='unknown', stringsAsFactors = FALSE)
    df <- rbind(df, df_non)
  }
  df <- df[order(df$Type),]
  rownames(df) <- 1:nrow(df)
  df

}

#' Is a value NA or zero.
#'
#' As title
#'
#'
#' @param x A numeric value.
#' @return TRUE or FALSE
#' @author T. Carruthers
#' @keywords internal
#' @export
NAor0 <- function(x) {
  if (length(x) == 0)
    return(TRUE)
  if (length(x) > 0)
    return(is.na(x[1]))
}


#' Print out plotting functions
#'
#' This function prints out the available plotting functions for objects of
#' class MSE or Data
#'
#' @param class Character string. Prints out the plotting functions for objects
#' of this class.
#' @param msg Logical. Should the functions be printed to screen?
#' @note Basically the function looks for any functions in the MSEtool that
#' have the word `plot` in them.  There is a chance that some plotting
#' functions are missed. Let us know if you find any and we will add them.
#' @author A. Hordyk
#' @export
plotFun <- function(class = c("MSE", "Data"), msg = TRUE) {
  class <- match.arg(class)
  tt <- lsf.str("package:MSEtool")
  p <- p2 <- rep(FALSE, length(tt))
  for (X in seq_along(tt)) {
    temp <- grep("plot", tolower(tt[[X]]))
    if (length(temp) > 0) p[X] <- TRUE
    temp2 <- grep(class, paste(format(match.fun(tt[[X]])), collapse = " "))
    if (length(temp2) > 0)  p2[X] <- TRUE
  }
  if (msg)
    message("MSEtool functions for plotting objects of class ", class,
            " are:")
  out <- sort(tt[which(p & p2)])

  out <- out[!out %in% c('plotFleet', 'plotStock', 'plotFun',
                         "COSEWIC_plot", "DFO_hist",
                         'plotOFL', "boxplot",
                         'boxplot.Data','plotDep','plotM2',
                         'plotGrowth', 'plotMat','plotRec', 'plot.OM')]

  if (class == "MSE") {
    out <- c(out, "VOI", "VOI2", "DFO_proj",
             "PWhisker")
    out <- sort(out)
  }
  if (class == "Data") {
    out <- c(out, "boxplot", "Sense")
    out <- sort(out)
  }
  if (msg) {
    if (length(out) > 5) {
      sq <- seq(from = 1, to = length(out), by = 5)
      for (x in seq_along(sq)) {
        cat(out[sq[x]:min(sq[x + 1] - 1, length(out), na.rm = TRUE)])
        cat("\n")
      }
    } else cat(out)
    cat("\n")
  }
  invisible(out)
}


#' What management procedures need what data
#'
#' A function that finds all the MPs and searches the
#' function text for slots in the Data object
#'
#' @param funcs A character vector of management procedures
#' @param noCV Logical. Should the CV slots be left out?
#' @author T. Carruthers
#' @return A matrix of MPs and their required data in terms of `slotnames('Data')`,
#' and broad Data classes for each MP
#' @seealso \link{Can} \link{Cant} \link{Needed} \link{MPtype} \linkS4class{Data}
#' @export
Required <- function(funcs = NA, noCV=FALSE) {

  if (!methods::is(funcs, 'logical') & !methods::is(funcs, "character"))
    stop("first argument must be character with MP name")

  if (all(is.na(funcs))) funcs <- avail("MP", msg=FALSE)
  for (x in 1:length(funcs)) {
    tt <- try(get(funcs[x]))
    if (!methods::is(tt,"MP")) stop(funcs[x], " is not class 'MP'")
  }

  ReqData <- MSEtool::ReqData
  builtin <- funcs[funcs %in% ReqData$MP]
  custom <- funcs[!funcs %in% ReqData$MP]

  df <- ReqData[match(builtin, ReqData$MP),]

  funcs1 <- custom

  temp <- lapply(funcs1, function(x) paste(format(match.fun(x)), collapse = " "))
  repp <- vapply(temp, match_slots, character(1))
  repp[!nzchar(repp)] <- "No data detected in this MP."

  df2 <- data.frame(MP=funcs1, Data=repp, stringsAsFactors = FALSE)

  dfout <- rbind(df, df2)

  if (noCV) {
    for (rr in 1:nrow(dfout)) {
      tt <- unlist(strsplit(dfout$Data[rr], ","))
      tt <- tt[!grepl("CV_", tt)]
      tt <- sort(trimws(sort(tt)))
      dfout$Data[rr] <- paste(tt, collapse = ", ")
    }
  }

  # Add Data Classes
  dfout$DataClass <- NULL
  for (i in 1:nrow(dfout)) {
    DataClass <- NULL

    # Catch Data
    if (grepl("LHYear", dfout$Data[i]) & grepl("Cat", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Recent Catches')
    }
    if (!grepl("LHYear", dfout$Data[i]) & grepl("Cat", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Historical Catches')
    }
    if (grepl("AvC", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Average Catch')
    }

    # Index
    if (grepl("Ind", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Index of Abundance')
    }

    # Rec
    if (grepl("Rec", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Recruitment Index')
    }

    if (grepl("Abun", dfout$Data[i]) | grepl("SpAbun", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Current Abundance')
    }

    if (grepl("Dep", dfout$Data[i]) | grepl("Dt", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Current Depletion')
    }

    if (grepl("Mort", dfout$Data[i]) | grepl("MaxAge", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Natural Mortality')
    }


    if (grepl("FMSY_M", dfout$Data[i]) | grepl("BMSY_B0", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Reference Ratios')
    }

    if (grepl("L50", dfout$Data[i]) | grepl("L95", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Maturity')
    }

    if (grepl("ML", dfout$Data[i]) | grepl("Lbar", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Mean Length')
    }

    if (grepl("Lc", dfout$Data[i]) | grepl("LFC", dfout$Data[i]) | grepl("LFS", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Selectivity')
    }

    if (grepl("CAA", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Age Composition')
    }

    if (grepl("vbK", dfout$Data[i]) | grepl("vbLinf", dfout$Data[i]) |
        grepl("vbt0", dfout$Data[i]) | grepl("LenCV", dfout$Data[i]) |
        grepl("wla", dfout$Data[i]) | grepl("wlb", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Growth Parameters')
    }

    if (grepl("steep", dfout$Data[i]) | grepl("sigmaR", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Stock-Recruitment')
    }

    if (grepl("CAL_bins", dfout$Data[i]) | grepl("CAL", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Length Composition')
    }

    if (grepl("Cref", dfout$Data[i]) | grepl("Iref", dfout$Data[i]) |
        grepl("Bref", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Reference Levels')
    }

    if (grepl("MPrec", dfout$Data[i]) | grepl("MPeff", dfout$Data[i])) {
      DataClass <- c(DataClass, 'Recent Management')
    }

    DataClass <- sort(DataClass)
    dfout$DataClass[i] <- paste0(DataClass, collapse=", ")

  }

  rownames(dfout) <- NULL

  data.frame(dfout)
}




#' Setup parallel processing
#'
#' Sets up parallel processing using the snowfall package
#'
#' @param cpus the number of CPUs to use for parallel processing. If left empty
#' all physical cores will be used, unless `logical=TRUE`, in which case both
#' physical and logical (virtual) cores will be used.
#' @param logical Use the logical cores as well? Using the virtual cores may
#' not lead to any significant decrease in run time.
#' You can test the optimal number of cores using `optCPU()`
#' @param ... other arguments passed to 'snowfall::sfInit'
#' @examples
#' \dontrun{
#' setup() # set-up the physical processors
#' setup(6) # set-up 6 processors
#' setup(logical=TRUE) # set-up physical and logical cores
#' }
#' @export
setup <- function(cpus=NULL, logical=FALSE, ...) {
  if (is.null(cpus))
    cpus <- parallel::detectCores(logical=logical)
  if(snowfall::sfIsRunning())
    snowfall::sfStop()
  snowfall::sfInit(parallel=TRUE,cpus=cpus, ...)
  sfLibrary("MSEtool", character.only = TRUE, verbose=FALSE)
  pkgs <- search()
  if ("package:DLMtool" %in% pkgs)
    sfLibrary("DLMtool", character.only = TRUE, verbose=FALSE)
  if ("package:SAMtool" %in% pkgs)
    sfLibrary("SAMtool", character.only = TRUE, verbose=FALSE)

}



# #' Open the DLMtool User Guide
# #'
# #' Opens the DLMtool User Guide website (requires internet connection)
# #'
# #' @export
# #' @examples
# #' \dontrun{
# #' userguide()
# #' }
# userguide <- function() {
#   utils::browseURL("https://dlmtool.github.io/DLMtool/userguide/introduction.html")
# }
#
# #' Opens the DLMtool Cheat-Sheets (requires internet connection)
# #'
# #' @export
# #' @examples
# #' \dontrun{
# #' cheatsheets()
# #' }
# cheatsheets <- function() {
#   # utils::browseURL("https://dlmtool.github.io/DLMtool/cheat_sheets/DLMtool_CheatSheets# .pdf")
#   utils::browseURL("https://dlmtool.github.io/DLMtool/cheat_sheets/CheatSheets.html")
# }

RepmissingVal <- function(object, name, vals=NA) {
  miss <- FALSE
  if (!name %in% slotNames(object)) return(object)
  if (!.hasSlot(object,name)) miss <- TRUE
  if (!miss) {
    if (length(slot(object, name))==0) miss <- TRUE
    if (all(is.na(slot(object, name)))) miss <- TRUE
  }
  if (miss) slot(object, name) <- vals
  return(object)
}

#' @describeIn checkMSE Updates an existing MSE object (class MSE) from a previous version of the
#' MSEtool to include slots new to the latest version. Also works with Stock,
#' Fleet, Obs, Imp, and Data objects. The new slots will be empty,
#' but avoids the 'slot doesn't exist' error that sometimes occurs.
#' Returns an object of class matching class(MSEobj)
#' @param save.name Character string. Optional file name to save the updated MSE object to disk.
#' @export
updateMSE <- function(MSEobj, save.name=NULL) {

  slots <- slotNames(MSEobj)

  if (length(slots)<1 & methods::is(MSEobj,'MSE')) {
    # incompatible version
    message('Updating MSE object from earlier version of MSEtool')
    nMSE <- new("MSE",
                Name=MSEobj@Name,
                nyears=MSEobj@nyears,
                proyears=MSEobj@proyears,
                nMPs=MSEobj@nMPs,
                MPs=MSEobj@MPs,
                nsim=MSEobj@nsim,
                OM=MSEobj@OM,
                Obs=MSEobj@Obs,
                SB_SBMSY=MSEobj@B_BMSY,
                F_FMSY=MSEobj@F_FMSY,
                N=array(),
                B=MSEobj@B,
                SSB=MSEobj@SSB,
                VB=MSEobj@VB,
                FM=MSEobj@FM,
                SPR=list(),
                Catch=MSEobj@C,
                Removals=array(),
                Effort=MSEobj@Effort,
                TAC=MSEobj@TAC,
                TAE=array(),
                BioEco=list(),
                RefPoint=list(MSEobj@Misc$MSYRefs),
                CB_hist=MSEobj@CB_hist,
                FM_hist=MSEobj@FM_hist,
                SSB_hist=MSEobj@SSB_hist,
                Hist=new('Hist'),
                PPD=MSEobj@Misc$Data,
                Misc=MSEobj@Misc
                )
    if (!is.null(save.name)) {

      message('Saving updated MSE object to: ', save.name)
      saveRDS(nMSE, save.name)
      message('Restart R session and load with `MyMSE <- readRDS(', save.name, ")`")
      return(invisible(nMSE))
    } else {
      return(nMSE)
    }

  }


  for (X in seq_along(slots)) {
    classDef <- getClassDef(class(MSEobj))
    slotTypes <- classDef@slots
    tt <- try(slot(MSEobj, slots[X]), silent = TRUE)
    if (inherits(tt, "try-error")) {
      fun <- get(as.character(slotTypes[X]))
      if(as.character(slotTypes[X]) == "vector") {
        slot(MSEobj, slots[X]) <- fun("numeric", length=0)
      } else slot(MSEobj, slots[X]) <- fun(0)
    }
  }
  MSEobj <- RepmissingVal(MSEobj, 'LenCV', c(0.08,0.15))
  MSEobj <- RepmissingVal(MSEobj, 'LR5', c(0,0))
  MSEobj <- RepmissingVal(MSEobj, 'LFR', c(0,0))
  MSEobj <- RepmissingVal(MSEobj, 'Rmaxlen', c(1,1))
  MSEobj <- RepmissingVal(MSEobj, 'DR', c(0,0))
  MSEobj <- RepmissingVal(MSEobj, 'Fdisc', c(0,0))
  MSEobj <- RepmissingVal(MSEobj, 'nareas', 2)

  MSEobj
}




#' Calculate CV from vector of values
#'
#'
#' @param x vector of numeric values
#' @author T. Carruthers
#' @return numeric
#' @keywords internal
#' @export
cv <- function(x) sd(x)/mean(x)


#' Get parameters of lognormal distribution from mean and standard deviation in normal
#' space
#'
#' @param m mean in normal space
#' @param sd standard deviation in normal space
#' @author T. Carruthers
#' @return numeric
#' @describeIn sdconv Returns sigma of lognormal distribution
#' @keywords internal
#' @export
sdconv <- function(m, sd) (log(1 + ((sd^2)/(m^2))))^0.5


#' @describeIn sdconv Returns mu of lognormal distribution
#' @export
mconv <- function(m, sd) log(m) - 0.5 * log(1 + ((sd^2)/(m^2)))

#' Calculate parameters for beta distribution from mean and standard deviation in
#' normal space
#'
#' @param m mean
#' @param sd standard deviation
#' @author T. Carruthers
#' @return numeric
#' @describeIn alphaconv Returns alpha of beta distribution
#' @keywords internal
#' @export
alphaconv <- function(m, sd) m * (((m * (1 - m))/(sd^2)) - 1)


#' @describeIn alphaconv Returns beta of beta distribution
#' @export
betaconv <- function(m, sd) (1 - m) * (((m * (1 - m))/(sd^2)) - 1)

#' Lognormal distribution 
#'
#' Variant of rlnorm which returns the mean when reps = 1.
#'
#' @param reps number of random numbers
#' @param mu mean
#' @param cv coefficient of variation
#' @param x vector
#' @author T. Carruthers
#' @return numeric
#' @describeIn trlnorm Generate log-normally distributed random numbers
#' @keywords internal
#' @export
trlnorm <- function(reps, mu, cv) {
  if (all(is.na(mu))) return(rep(NA, reps))
  if (all(is.na(cv))) return(rep(NA, reps))
  if (reps == 1)  return(mu)
  return(rlnorm(reps, mconv(mu, mu * cv), sdconv(mu, mu * cv)))
}



#' @describeIn trlnorm Calculate density of log-normally distributed random numbers
#' @export
tdlnorm <- function(x, mu, cv) dlnorm(x, mconv(mu, mu * cv), sdconv(mu, mu * cv))


#' Stock recruit parameterization
#' 
#' Convert stock recruit parameters from steepness parameterization to alpha/beta (and vice versa)
#' 
#' @param alpha Alpha parameter
#' @param beta Beta parameter
#' @param h Steepness parameter
#' @param R0 Unfished recruitment parameter
#' @param phi0 Unfished spawners per recruit
#' @param SR Stock-recruit function: (1) Beverton-Holt, or (2) Ricker
#' @param type The parameterization of the Beverton-Holt function with respect to \code{alpha} 
#' and \code{beta}. See details.
#' @details
#' The Type 1 Beverton-Holt equation is
#' \deqn{R = \alpha S/(1 + \beta S)}
#' 
#' The Type 2 Beverton-Holt equation is
#' \deqn{R = S/(\alpha + \beta S)}
#'  
#' The Ricker equation is
#' \deqn{R = \alpha S \exp(-\beta S)}
#' @return A numeric.
#' @author Q. Huynh
#' @describeIn hconv Returns steepness (h) from \code{alpha} and \code{phi0}
#' @export
hconv <- function(alpha, phi0, SR = 1, type = 1) {
  if(SR == 1) {
    type <- match.arg(as.character(type), choices = as.character(1:2))
    h <- switch(type,
                "1" = alpha*phi0/(4 + alpha*phi0),
                "2" = phi0/(4*alpha + phi0))
  } else {
    h <- 0.2 * (alpha*phi0)^0.8
  }
  return(h)
}

#' @describeIn hconv Returns unfished recruitment (R0) from \code{alpha}, \code{beta}, and \code{phi0}
#' @export
R0conv <- function(alpha, beta, phi0, SR = 1, type = 1) {
  if(SR == 1) {
    type <- match.arg(as.character(type), choices = as.character(1:2))
    R0 <- switch(type,
                "1" = (alpha*phi0 - 1)/beta/phi0,
                "2" = (phi0 - alpha)/beta/phi0)
  } else {
    R0 <- log(alpha*phi0)/beta/phi0
  }
  return(R0)
}

#' @describeIn hconv Returns \code{alpha} from \code{h} and \code{phi0}
#' @export
SRalphaconv <- function(h, phi0, SR = 1, type = 1) {
  if(SR == 1) {
    type <- match.arg(as.character(type), choices = as.character(1:2))
    alpha <- switch(type,
                    "1" = 4*h/(1-h)/phi0,
                    "2" = phi0*(1-h)/(4*h))
  } else {
    alpha <- (5*h)^1.25/phi0
  }
  return(alpha)
}

#' @describeIn hconv Returns \code{beta} from \code{h}, \code{R0}, and \code{phi0}
#' @export
SRbetaconv <- function(h, R0, phi0, SR = 1, type = 1) {
  if(SR == 1) {
    type <- match.arg(as.character(type), choices = as.character(1:2))
    beta <- switch(type,
                   "1" = (5*h-1)/(1-h)/phi0/R0,
                   "2" = (5*h-1)/(4*h)/R0)
  } else {
    beta <- log((5*h)^1.25)/phi0/R0
  }
  return(beta)
}

#' Depletion and F estimation from mean length of catches
#'
#' A highly dubious means of getting very uncertain estimates of current stock
#' biomass and (equilibrium) fishing mortality rate from growth, natural
#' mortality rate, recruitment and fishing selectivity.
#'
#' @param OM An object of class 'OM'
#' @param ML A estimate of current mean length of catches
#' @param nsim Number of simulations
#' @param ploty Produce a plot of depletion and F
#' @param Dlim Limits on the depletion that is returned as a fraction of
#' unfished biomass.
#' @return An object of class 'OM' with 'D' slot populated
#' @author T. Carruthers
#' @export
ML2D <- function(OM, ML, nsim = 100, ploty = T, Dlim = c(0.05, 0.6)) {

  nsim2<-nsim*10
  maxage <- OM@maxage
  M <- runif(nsim2, OM@M[1], OM@M[2])  # Natural mortality rate
  h <- runif(nsim2, OM@h[1], OM@h[2])  # Steepness
  Linf <- runif(nsim2, OM@Linf[1], OM@Linf[2])  # Maximum length
  K <- runif(nsim2, OM@K[1], OM@K[2])  # Maximum growth rate
  t0 <- runif(nsim2, OM@t0[1], OM@t0[2])  # Theorectical length at age zero

  if (OM@isRel == "0" | OM@isRel == "FALSE" | OM@isRel == FALSE) {
    if (max(OM@LFS) > 0) {
      LFS <- runif(nsim2*5, OM@LFS[1], OM@LFS[2])
    } else {
      LFS <- runif(nsim2*5, mean(OM@LFSLower), mean(OM@LFSUpper))
    }
  } else {
    if (max(OM@LFS) > 0) {
      LFS <- runif(nsim2*5, OM@LFS[1], OM@LFS[2]) * mean(OM@L50)
    } else {
      LFS <- runif(nsim2*5, mean(OM@LFSLower), mean(OM@LFSUpper)) *
        mean(OM@L50)
    }
  }
  LFS<-LFS[LFS<Linf][1:nsim2]
  AFS <- L2A(t0, Linf, K, LFS, maxage)
  AFS[AFS<1]<-1

  if (OM@isRel == "0" | OM@isRel == "FALSE" | OM@isRel == FALSE) {
    L5 <- runif(nsim2, OM@L5[1], OM@L5[2])
  }else{
    L5 <- runif(nsim2, OM@L5[1], OM@L5[2])* mean(OM@L50)
  }

  age05 <- L2A(t0, Linf, K, L5, maxage)
  age05[age05<0.5] <- 0.5

  Vmaxage <- runif(nsim2, OM@Vmaxlen[1], OM@Vmaxlen[2])  #runif(BT_fleet@Vmaxage[1],BT_fleet@Vmaxage[2]) # selectivity of oldest age class

  LM <- runif(nsim2*5, OM@L50[1], OM@L50[2])
  LM<-LM[LM<Linf][1:nsim2]
  AM <- L2A(t0, Linf, K, LM, maxage)
  AM[AM<1]<-1

  # age at maturity
  a <- OM@a  # length-weight parameter a
  b <- OM@b  # length-weight parameter b

  mod <- AFS  # the age at modal (or youngest max) selectivity

  # deriv <- getDNvulnS(mod, age05, Vmaxage, maxage, nsim2)  # The vulnerability schedule
  # vuln <- deriv[[1]]

  srs <- (maxage - AFS) / ((-log(Vmaxage,2))^0.5) # selectivity parameters are constant for all years
  sls <- (AFS - age05) /((-log(0.05,2))^0.5)

  vuln <- t(sapply(1:nsim2, getsel, lens=matrix(1:maxage, nrow=nsim2, ncol=maxage+1, byrow=TRUE),
                   lfs=AFS, sls=sls, srs=srs))

  Agearray <- array(rep(0:maxage, each = nsim2), c(nsim2, maxage+1))
  mat <- 1/(1 + exp((AM - (Agearray))/(AM * 0.1)))  # Maturity at age array

  nyears <- OM@nyears
  # bootfun<-function(dat,ind)mean(dat[ind]) MLo<-boot(MLt,bootfun,nsim2)
  # ML<-MLo$t
  out <- CSRA(M, h, Linf, K, t0, AM, a, b, vuln, mat, ML = rep(ML, nsim2),
              NA, NA, maxage, nyears)
  cond <- out[, 1] > Dlim[1] & out[, 1] < Dlim[2] & out[, 2] < 2.5  # Stock levels are unlikely to be above 80% unfished, F is unlikely to be above 2.5

  if (ploty & sum(cond) > 5) {
    par(mfrow = c(1, 2))
    plot(density(out[cond, 1], from = 0, adj = 0.4), main = "Depletion")
    plot(density(out[cond, 2], from = 0, adj = 0.4), main = "Fishing mortality rate")
  }
  if (sum(cond) < 5) {
    message("All estimates of Depletion outside bounds of Dlim")
    message("Operating Model object not updated")
  }

  if(sum(cond)>nsim){
    OM@cpars$D<-out[cond,1][1:nsim]
  }

  if(sum(cond) > 5) OM@D <- quantile(out[cond, 1], c(0.05, 0.95))

  OM
}

# Composition stock reduction analysis


#' Catch at size reduction analysis
#'
#' What depletion level and corresponding equilibrium F arise from data
#' regarding mean length of current catches, natural mortality rate, steepness
#' of the stock recruitment curve, maximum length, maximum growth rate, age at
#' maturity, age based vulnerability, maturity at age, maximum age and number
#' of historical years of fishing.
#'
#'
#' @param M A vector of natural mortality rate estimates
#' @param h A vector of sampled steepness (Beverton-Holt stock recruitment)
#' @param Linf A vector of maximum length (von Bertalanffy growth)
#' @param K A vector of maximum growth rate (von Bertalanffy growth)
#' @param t0 A vector of theoretical age at length zero (von Bertalanffy
#' growth)
#' @param AM A vector of age at maturity
#' @param a Length-weight conversion parameter a (W=aL^b)
#' @param b Length-weight conversion parameter b (W=aL^b)
#' @param vuln A matrix nsim x nage of the vulnerability at age (max 1) to
#' fishing.
#' @param mat A matrix nsim x nage of the maturity at age (max 1)
#' @param ML A vector of current mean length estimates
#' @param CAL A catch-at-length matrix nyears x (1 Linf unit) length bins
#' @param CAA A catch-at-age matrix nyears x maximum age
#' @param maxage Maximum age
#' @param nyears Number of historical years of fishing
#' @author T. Carruthers
#' @export CSRA
#' @keywords internal
CSRA <- function(M, h, Linf, K, t0, AM, a, b, vuln, mat, ML, CAL, CAA,
                 maxage, nyears) {
  nsim <- length(M)
  Dep <- rep(NA, nsim)
  Fm <- rep(NA, nsim)
  for (i in 1:nsim) {
    fit <- optimize(CSRAfunc, log(c(1e-04, 5)), Mc = M[i], hc = h[i],
                    maxage, nyears, Linfc = Linf[i], Kc = K[i], t0c = t0[i], AMc = AM[i],
                    ac = a, bc = b, vulnc = vuln[i, ], matc = mat[i, ], MLc = ML[i],
                    CAL = NA, CAA = NA, opt = T)


    out <- CSRAfunc(fit$minimum, Mc = M[i], hc = h[i], maxage, nyears,
                    Linfc = Linf[i], Kc = K[i], t0c = t0[i], AMc = AM[i], ac = a,
                    bc = b, vulnc = vuln[i, ], matc = mat[i, ], MLc = ML[i], CAL = NA,
                    CAA = NA, opt = 3)

    Dep[i] <- out[1]
    Fm[i] <- out[2]


  }
  cbind(Dep, Fm)
}

# The function that CSRA operates on

#' Optimization function for CSRA
#'
#' What depletion level and corresponding equilibrium F arise from data
#' regarding mean length of current catches, natural mortality rate, steepness
#' of the stock recruitment curve, maximum length, maximum growth rate, age at
#' maturity, age based vulnerability, maturity at age, maximum age and number
#' of historical years of fishing.
#'
#'
#' @param lnF A proposed value of current instantaneous fishing mortality rate
#' @param Mc Natural mortality rate estimates
#' @param hc Steepness (Beverton-Holt stock recruitment)
#' @param maxage Maximum age
#' @param nyears Number of historical years of fishing
#' @param AFSc Age at full selection
#' @param AFCc Age at first capture
#' @param Linfc Maximum length (von Bertalanffy growth)
#' @param Kc Maximum growth rate (von Bertalanffy growth)
#' @param t0c Theoretical age at length zero (von Bertalanffy growth)
#' @param AMc Age at maturity
#' @param ac Length-weight conversion parameter a (W=aL^b)
#' @param bc Length-weight conversion parameter b (W=aL^b)
#' @param vulnc A vector (nage long) of the vulnerability at age (max 1) to
#' fishing.
#' @param matc A vector (nage long) of the maturity at age (max 1)
#' @param MLc A current mean length estimates
#' @param CAL A catch-at-length matrix nyears x (1 Linf unit) length bins
#' @param CAA A catch-at-age matrix nyears x maximum age
#' @param opt Should the measure of fit be returned?
#' @param meth Are we fitting to mean length or catch composition?
#' @author T. Carruthers
#' @keywords internal
CSRAfunc <- function(lnF, Mc, hc, maxage, nyears, AFSc, AFCc, Linfc, Kc,
                     t0c, AMc, ac, bc, vulnc, matc, MLc, CAL, CAA, opt = T, meth = "ML") {

  Fm <- exp(lnF)
  Fc <- vulnc * Fm
  Lac <- Linfc * (1 - exp(-Kc * ((0:maxage) - t0c)))
  Wac <- ac * Lac^bc
  N <- exp(-Mc * ((0:maxage) - 1))
  SSN <- matc * N  # Calculate initial spawning stock numbers
  Biomass <- N * Wac
  SSB <- SSN * Wac  # Calculate spawning stock biomass

  B0 <- sum(Biomass)
  SSB0 <- sum(SSB)
  SSN0 <- SSN
  SSBpR <- sum(SSB)  # Calculate spawning stock biomass per recruit
  SSNpR <- SSN
  Zc <- Fc + Mc
  n_age <- maxage+1
  CN <- array(NA, dim = c(nyears, n_age))
  HR <- rep(0, n_age)
  pen <- 0
  for (y in 1:nyears) {
    VB <- Biomass * vulnc * exp(-Mc)
    CN[y, ] <- N * (1 - exp(-Zc)) * (Fc/Zc)
    N[2:n_age] <- N[1:(n_age - 1)] * exp(-Zc[1:(n_age - 1)])  # Total mortality
    N[1] <- (0.8 * hc * sum(SSB))/(0.2 * SSBpR * (1 - hc) + (hc - 0.2) *
                                     sum(SSB))  # Recruitment assuming regional R0 and stock wide steepness
    Biomass <- N * Wac
    SSN <- N * matc
    SSB <- SSN * Wac
  }  # end of year

  pred <- sum((CN[nyears, ] * Lac))/sum(CN[nyears, ])
  fobj <- (pred - MLc)^2  # Currently a least squares estimator. Probably not worth splitting hairs WRT likelihood functions!
  if (opt == 1) {
    return(fobj)
  } else {
    c(sum(SSB)/sum(SSB0), Fm)
  }
}


# Stochastic inverse growth curve used to back-calculate age at first
# capture from length at first capture
#' Calculate age at first capture from length at first capture and growth
#'
#' As title.
#'
#' @param t0c A vector of theoretical age at length zero (von Bertalanffy
#' growth)
#' @param Linfc A vector of maximum length (von Bertalanffy growth)
#' @param Kc A vector of maximum growth rate (von Bertalanffy growth)
#' @param LFC A vector of length at first capture
#' @param maxage Maximum age
#' @author T. Carruthers
#' @keywords internal
getAFC <- function(t0c, Linfc, Kc, LFC, maxage) {
  nsim <- length(t0c)
  agev <- c(1e-04, 0:maxage)
  agearray <- matrix(rep(agev, each = nsim), nrow = nsim)
  Larray <- Linfc * (1 - exp(-Kc * (agearray - t0c)))
  matplot(agev, t(Larray), type = "l")
  abline(h = LFC, col = "#ff000030", lwd = 2)
  AFC <- (log(1 - (LFC/Linfc))/-Kc) + t0c
  abline(v = AFC, col = "#0000ff30", lwd = 2)
  AFC
}



#' Length to age conversion
#'
#' Simple deterministic length to age conversion given inverse von Bertalanffy
#' growth.
#'
#' @param t0c Theoretical age at length zero
#' @param Linfc Maximum length
#' @param Kc Maximum growth rate
#' @param Len Length
#' @param maxage Maximum age
#' @param ploty Should a plot be included
#' @return An age (vector of ages, matrix of ages) corresponding with Len
#' @author T. Carruthers
#' @keywords internal
L2A <- function(t0c, Linfc, Kc, Len, maxage, ploty=F) {
  nsim <- length(t0c)
  agev <- c(1e-04, 0:maxage)
  agearray <- matrix(rep(agev, each = nsim), nrow = nsim)
  Larray <- Linfc * (1 - exp(-Kc * (agearray - t0c)))
  temp<-Len/Linfc
  temp[temp<0.95]<-0.95
  age <- (log(1 - (Len/Linfc))/-Kc) + t0c
  if(ploty){
    matplot(agev, t(Larray), type = "l")
    abline(h = Len, col = "#ff000030", lwd = 2)
    abline(v = age, col = "#0000ff30", lwd = 2)
  }
  age
}



#' Determine optimal number of cpus
#'
#' @param nsim Numeric. Number of simulations.
#' @param thresh Recommended n cpus is what percent of the fastest time?
#' @param plot Logical. Show the plot?
#' @param msg Logical. Should messages be printed to console?
#' @param maxn Optional. Maximum number of cpus. Used for demo purposes
#'
#' @templateVar url parallel-processing
#' @templateVar ref determining-optimal-number-of-processors
# #' @template userguide_link
#'
#' @export
#' @seealso \link{setup}
#' @examples
#' \dontrun{
#' optCPU()
#' }
#' @author A. Hordyk
optCPU <- function(nsim=96, thresh=5, plot=TRUE, msg=TRUE, maxn=NULL) {
  cpus <- 1:parallel::detectCores()
  if (!is.null(maxn)) cpus <- 1:maxn

  time <- NA
  OM <- MSEtool::testOM
  OM@nsim <- nsim
  for (n in cpus) {
    if (msg) message('Running MSE with ', nsim, ' simulations and ', n, ' of ', max(cpus), ' cpus')
    if (n == 1) {
      snowfall::sfStop()
      st <- Sys.time()
      tt <- runMSE(OM, silent = TRUE)
      time[n] <- difftime(Sys.time(), st, units='secs')
    } else{
      if (msg) {
        setup(cpus=n)
      } else {
        sink('temp')
        suppressMessages(setup(cpus=n))
        sink()
      }
      st <- Sys.time()
      tt <- runMSE(OM, silent=TRUE, parallel=TRUE)

      time[n] <- difftime(Sys.time(), st, units='secs')

    }
  }
  df <- data.frame(ncpu=cpus, time=time)
  df$time <- round(df$time,2)
  rec <- min(which(time < min(time) * (1 + thresh/100)))
  if (plot) {
    plot(df, type='b', ylab="time (seconds)", xlab= "# cpus", bty="l", lwd=2)
    points(rec, df[rec,2], cex=2, pch=16, col="blue")
  }
  return(df)
}


# #' Convert an MMSE object to an MSE object
# #'
# #' @param MMSE An object of class `MMSE`
# #' @param B Numeric. Which stock should be used to report biomass? Use NA to calculate mean for all stocks.
# #' @param C Numeric. Which fleet should be used to report catch? Use NA to calculate sum for all fleets.
# #' @param F Numeric. Which stock should be used to report fishing mortality? Use NA to calculate mean for all stocks.
# #' @param E Numeric. Which stock should be used to report effort? Use NA to calculate mean for all stocks.
# #'
# #' @return an object of class `MSE`
# #' @export
# #'
# Convert <- function(MMSE, B=1, C=NA, F=1, E=1) {
#   MPs <- MMSE@MPs
#   if (class(MPs) == 'list') {
#     # MPs by stock
#     if (class(MPs[[1]][[1]]) == 'list') {
#       # by fleet as well
#       stop("Function needs updating")
#     } else {
#       stocks <- 1:MMSE@nstocks
#       df <- data.frame(MPs)
#       colnames(df) <- stocks
#       MPnames <- rep('', nrow(df))
#       for (i in 1:nrow(df)) {
#         mps <- as.character(df[i,] )
#         MPnames[i] <- paste(unique(mps), collapse="-")
#       }
#     }
#   }
#
#   OM <- MMSE@OM[[1]][[1]]
#   ns <- MMSE@nstocks
#   nf <- MMSE@nfleets
#   nsim <- MMSE@nsim
#
#   refY <- array(NA, dim=c(nsim, ns, nf))
#   for (s in 1:ns){
#     for (f in 1:nf) {
#       refY[,s,f] <- MMSE@OM[[s]][[f]]$RefY
#     }
#   }
#
#   OM$RefY <- apply(refY, 1, sum)
#
#   if(is.na(B)) B <- 1:ns
#   if(is.na(C)) C <- 1:nf
#   if(is.na(F)) F <- 1:ns
#
#   B_BMSY <- apply(MMSE@B_BMSY[,B,,, drop=FALSE], c(1,3,4), mean)
#   F_FMSY  <- apply(MMSE@F_FMSY[,F,1:nf,, ,drop=FALSE], c(1,4,5), mean)
#   Catch <- apply(MMSE@C[,1:ns,C,, ,drop=FALSE], c(1,4,5), sum)
#   Effort <- apply(MMSE@Effort[,E,1:nf,, ,drop=FALSE], c(1,4,5), mean)
#
#   Biomass <-apply(MMSE@B[,B,,, drop=FALSE], c(1,3,4), mean)
#   SSB <- apply(MMSE@SSB[,B,,, drop=FALSE], c(1,3,4), mean)
#   VB <- apply(MMSE@VB[,B,,, drop=FALSE], c(1,3,4), mean)
#   FM <- apply(MMSE@FM[,F,1:nf,,, drop=FALSE], c(1,4,5), mean)
#   TAC <- apply(MMSE@TAC[,1:ns, C,,,drop=FALSE], c(1,4,5), sum)
#   SSB_hist <- apply(MMSE@SSB_hist[,B,,,,drop=FALSE], c(1,3,4,5), mean)
#   CB_hist <- apply(MMSE@CB_hist[,1:ns,C,,,,drop=FALSE], c(1,4,5,6), sum)
#   FM_hist <- apply(MMSE@FM_hist[,1:ns,C,,,,drop=FALSE], c(1,4,5,6), sum)
#
#   MSE <- new('MSE', MMSE@Name,
#              nyears=MMSE@nyears,
#              proyears=MMSE@proyears,
#              nMPs=MMSE@nMPs,
#              MPs=MPnames,
#              nsim=MMSE@nsim,
#              OM=OM,
#              Obs=MMSE@Obs[[1]][[1]],
#              SB_SBMSY=B_BMSY,
#              F_FMSY=F_FMSY,
#              N=array(),
#              B=Biomass,
#              SSB=SSB,
#              VB=VB,
#              FM=FM,
#              SPR=list(),
#              Catch,
#              Removals=array(),
#              TAC=TAC,
#              TAE=array(),
#              BioEco=list(),
#              RefPoint=list(),
#              Hist=new('Hist'),
#              PPD=list(),
#              Misc=list()
#   )
#
#   MSE
# }
#

#' Get help topic URL
#'
#' @param topic Name of the functions
#' @param url URL for the help documentation
#' @param nameonly Logical. Help file name only?
#'
#' @return file path to help file
#' @export
#'
#' @keywords internal
MPurl <- function(topic, url='https://dlmtool.openmse.com/reference/',
                  nameonly=FALSE) {

  paths <- file.path(.libPaths()[1], "DLMtool")

  res <- character()
  for (p in paths) {
    if (file.exists(f <- file.path(p, "help", "aliases.rds")))
      al <- readRDS(f)
    else if (file.exists(f <- file.path(p, "help", "AnIndex"))) {
      foo <- scan(f, what = list(a = "", b = ""), sep = "\t",
                  quote = "", na.strings = "", quiet = TRUE)
      al <- structure(foo$b, names = foo$a)
    }
    else next
    f <- al[topic]
    if (is.na(f))
      next
    res <- c(res, file.path(p, "help", f))

  }
  if (length(res)<1) return(NA)

  if(nameonly) {
    return(basename(res))
  } else{
    return(paste0(url, basename(res), ".html"))
  }

}


## Sampling steepness parameter (h)  ####
#
# Code from Quang Huynh that fixes the bug where h is sometimes sampled > 1 or < 0.2

#' Sample steepness given mean and cv
#'
#' @param n number of samples
#' @param mu mean h
#' @param cv cv of h
#'
#' @author Q. Huynh
#' @export
#' @keywords internal
#' @describeIn sample_steepness2 sample steepness values
sample_steepness2 <- function(n, mu, cv) {

  if(n == 1) return(mu)
  else {
    sigma <- mu * cv

    mu.beta.dist <- (mu - 0.2)/0.8
    sigma.beta.dist <- sigma/0.8

    beta.par <- derive_beta_par(mu.beta.dist, sigma.beta.dist)

    h.transformed <- rbeta(n, beta.par[1], beta.par[2])

    h <- 0.8 * h.transformed + 0.2
    h[h > 0.99] <- 0.99
    h[h < 0.2] <- 0.2

    return(h)
  }

}


#' This function reduces the CV by 5 per cent until steepness values can be sampled without error
#'
#'
#' @param mu mean h
#' @param sigma sd of h
#'
#' @describeIn sample_steepness2 derive beta parameter
#' @export
derive_beta_par <- function(mu, sigma) {

  a <- alphaconv(mu, sigma)
  b <- betaconv(mu, sigma)

  if(a <= 0 || b <= 0) {
    sigma <- 0.95 * sigma
    Recall(mu, sigma)
  }
  else return(c(a, b))

}



#' TAC Filter
#'
#' Filters vector of TAC recommendations by replacing negatives with NA and
#' and values beyond five standard deviations from the mean as NA
#'
#' @param TAC A numeric vector of TAC recommendations
#' @author T. Carruthers
#' @export
TACfilter <- function(TAC) {
  TAC[TAC < 0] <- NA  # Have to robustify due to R optmization problems.. work in progress.
  TAC[TAC > (mean(TAC, na.rm = T) + 5 * stats::sd(TAC, na.rm = T))] <- NA  # remove very large TAC samples
  return(as.numeric(TAC))
}

#' Select DataList for an MP from `MMSE@PPD`
#'
#'
#' @param PPD `PPD` slot from an MMSE object 
#' @param MP Numeric value indicating the MP to return DataList
#' @return A nested list `Data` objects (`nstock` by `nfleet`)
#' @export
select_MP <- function(PPD, MP=1) {
  DataList <- vector('list', length(PPD))
  for (s in 1:length(DataList)) {
    DataList[[s]] <- vector('list', length(PPD[[s]]))
    for (f in 1:length(DataList[[s]])) {
      DataList[[s]][[f]] <- PPD[[s]][[f]][[MP]]
    }
  }
  DataList
}