MSEtool: Management Strategy Evaluation Toolkit

Documented in make_MP

#' Make a custom management procedure (MP)
#'
#' Function operator that combines a function of class \code{Assess} and a function of
#' class \code{HCR} to create a management procedure (MP). The resulting function
#' can then be tested in closed-loop simulation via \code{\link[DLMtool]{runMSE}}.
#'
#' @param .Assess A function of class \code{Assess}.
#' @param .HCR A function of class \code{HCR}.
#' @param diagnostic A character string describing if any additional diagnostic information from the
#' assessment models will be collected during a call with \code{\link[DLMtool]{runMSE}} (\code{"none"} is the
#' default). \code{"min"} (minimal) will collect information on convergence and \code{"full"} will also collect the
#' Assessment object generated by the \code{.Assess}. This information will be written to the \code{Misc} slot in the
#' \linkS4class{MSE} object.
#' See example.
#' @param ... Additional arguments to be passed to \code{.Assess} and \code{.HCR}.
#'
#' @return A function of class \code{MP}.
#'
#' @examples
#' # A delay-difference model with a 40-10 control rule
#' DD_40_10 <- make_MP(DD_TMB, HCR40_10)
#'
#' # A delay difference model that will produce convergence diagnostics
#' DD_40_10 <- make_MP(DD_TMB, HCR40_10, diagnostic = "min")
#'
#' # MP that uses a Delay-Difference which assumes a Ricker stock-recruit function.
#' DD_Ricker <- make_MP(DD_TMB, HCR_MSY, SR = "Ricker")
#'
#' \dontrun{
#' myMSE <- DLMtool::runMSE(DLMtool::testOM, MPs = c("FMSYref", "DD_40_10"), PPD = TRUE)
#'
#' str(myMSE@@Misc)
#' diagnostic_AM(myMSE)
#' }
#'
#' @seealso \link{HCR_ramp} \link{HCR_MSY} \link{diagnostic_AM} \link{retrospective_AM}
#' @export
make_MP <- function(.Assess, .HCR, diagnostic = c("none", "min", "full"), ...) {
  diagnostic <- match.arg(diagnostic)
  if(is.character(.Assess)) {
    Assess_char <- .Assess
    .Assess <- as.symbol(.Assess)
  } else {
    .Assess <- substitute(.Assess)
    Assess_char <- as.character(.Assess)
  }
  if(is.character(.HCR)) {
    .HCR <- as.symbol(.HCR)
  } else {
    .HCR <- substitute(.HCR)
  }
  if(!inherits(eval(.Assess), "Assess")) {
    stop(paste(.Assess, "does not belong to class 'Assess'. Use: avail('Assess') to find eligible objects."))
  }
  if(!inherits(eval(.HCR), "HCR")) {
    stop(paste(.HCR, "does not belong to class 'HCR.' Use: avail('HCR') to find eligible objects."))
  }

  dots <- list(...)
  dots_in_Assess <- dots[match(names(formals(eval(.Assess))), names(dots), nomatch = 0)]
  dots_in_HCR <- dots[match(names(formals(eval(.HCR))), names(dots), nomatch = 0)]

  Assess_call <- as.call(c(.Assess, x = quote(x), Data = quote(Data), dots_in_Assess))
  HCR_call <- as.call(c(.HCR, Assessment = quote(do_Assessment), reps = quote(reps), dots_in_HCR))

  if(diagnostic == "none") {

    MP_body <- bquote({
      dependencies <- .(get_dependencies(Assess_char, dots))
      do_Assessment <- .(Assess_call)
      Rec <- .(HCR_call)
      return(Rec)
    })

  } else {

    MP_body <- bquote({
      dependencies <- .(get_dependencies(Assess_char, dots))
      do_Assessment <- .(Assess_call)
      Rec <- .(HCR_call)
      Rec@Misc <- Assess_diagnostic(x, Data, do_Assessment, include_assessment = .(diagnostic == "full"))
      return(Rec)
    })

  }

  custom_MP <- eval(call("function", as.pairlist(alist(x = , Data = , reps = 1)), MP_body))
  class(custom_MP) <- "MP"
  return(custom_MP)
}


#' Model-based management procedures
#'
#' A suite of model-based management procedures (MPs) included in the package. Additional MPs,
#' with specific model configurations (e.g., stock-recruit function or fixing certain parameters) or alternative
#' ramped harvest control rules can be created with \link{make_MP} and the available Assess and HCR objects.
#'
#' @name Model-based-MP
#' @aliases MP Data-rich-MP
#' @param x A position in the Data object.
#' @param Data An object of class Data
#' @param reps Numeric, the number of stochastic replicates for the management advice.
#'
#' @examples
#' avail("MP", all_avail = FALSE)
#'
#' \dontrun{
#' myMSE <- DLMtool::runMSE(DLMtool::testOM, MPs = c("FMSYref", "SCA_MSY", "SCA_4010"))
#' }
#' @return An object of class \linkS4class{Rec} which contains the management recommendation.
NULL




#temp <- lapply(avail("Assess"), function(x) paste(format(match.fun(x)), collapse = " "))
#Assess_dep <- lapply(temp, DLMtool:::match_slots)
#names(Assess_dep) <- avail("Assess")
Assess_dep <- list(DD_SS = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   DD_TMB = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   SCA = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@L95, Data@CAA, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   SCA2 = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@L95, Data@CAA, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   SCA_Pope = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@L95, Data@CAA, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   VPA = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@L95, Data@CAA, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   SP = "Data@Cat, Data@Ind",
                   SP_SS = "Data@Cat, Data@Ind, Data@CV_Ind",
                   cDD = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   cDD_SS = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge",
                   SSS = "Data@Cat, Data@steep, Data@Mort, Data@L50, Data@L95, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge")


get_dependencies <- function(Assess, arg = list()) {
  dep <- getElement(Assess_dep, Assess)
  formals2 <- formals(Assess)
  arg_match <- pmatch(names(arg), names(formals2))
  formals2[arg_match[!is.na(arg_match)]] <- arg[!is.na(arg_match)]

  dep_match <- pmatch(names(dep_args), names(formals2)) # From formals
  more_deps <- dep_args[do.call(c, formals2[dep_match])]

  if(length(more_deps) > 0) {
    #if(Assess == "DD_SS" && any(names(more_deps) == "fix_sigma")) more_deps$fix_sigma <- "Data@CV_Cat"
    more_deps <- paste(do.call(c, more_deps), collapse = ", ")
    dep <- paste(c(dep, more_deps), collapse = ", ")
  }
  return(dep)
}

dep_args <- list(fix_h = "Data@steep", fix_sigma = "Data@CV_Ind", fix_tau = "Data@sigmaR", fix_omega = "Data@CV_Cat",
                 use_r_prior = "Data@Mort, Data@CV_Mort, Data@steep, Data@CV_steep, Data@vbLinf, Data@vbK, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge, Data@L50, Data@L95")




#' @describeIn Model-based-MP A statistical catch-at-age model with a TAC recommendation based on fishing at FMSY,
#' and default arguments for configuring \link{SCA}.
#' @export
SCA_MSY <- make_MP(SCA, HCR_MSY, diagnostic = "min")


#' @describeIn Model-based-MP An SCA with a TAC recommendation based on fishing at 75\% of FMSY.
#' @export
SCA_75MSY <- make_MP(SCA, HCR_MSY, MSY_frac = 0.75, diagnostic = "min")


#' @describeIn Model-based-MP An SCA with a 40-10 control rule.
#' @export
SCA_4010 <- make_MP(SCA, HCR40_10, diagnostic = "min")


#' @describeIn Model-based-MP A state-space delay difference model with a TAC recommendation based on fishing at FMSY,
#' and default arguments for configuring \link{DD_SS}.
#' @export
DDSS_MSY <- make_MP(DD_SS, HCR_MSY, diagnostic = "min")


#' @describeIn Model-based-MP A state-space delay difference model with a TAC recommendation based on fishing at 75\% of FMSY.
#' @export
DDSS_75MSY <- make_MP(DD_SS, HCR_MSY, MSY_frac = 0.75, diagnostic = "min")


#' @describeIn Model-based-MP A state-space delay difference model with a 40-10 control rule.
#' @export
DDSS_4010 <- make_MP(DD_SS, HCR40_10, diagnostic = "min")


#' @describeIn Model-based-MP A surplus production model with a TAC recommendation based on fishing at FMSY,
#' and default arguments for configuring \link{SP}.
#' @export
SP_MSY <- make_MP(SP, HCR_MSY, diagnostic = "min")


#' @describeIn Model-based-MP A surplus production model with a TAC recommendation based on fishing at 75\% of FMSY.
#' @export
SP_75MSY <- make_MP(SP, HCR_MSY, MSY_frac = 0.75, diagnostic = "min")


#' @describeIn Model-based-MP A surplus production model with a 40-10 control rule.
#' @export
SP_4010 <- make_MP(SP, HCR40_10, diagnostic = "min")


#' @describeIn Model-based-MP Simple stock synthesis (terminal depletion fixed to 0.4 in \link{SSS}) with a TAC recommendation based on fishing at FMSY.
#' @export
SSS_MSY <- make_MP(SSS, HCR_MSY, diagnostic = "min")

#' @describeIn Model-based-MP Simple stock synthesis (terminal depletion fixed to 0.4) with with a TAC recommendation based on fishing at 75\% FMSY.
#' @export
SSS_75MSY <- make_MP(SSS, HCR_MSY, MSY_frac = 0.75, diagnostic = "min")

#' @describeIn Model-based-MP Simple stock synthesis (terminal depletion fixed to 0.4) with a 40-10 control rule.
#' @export
SSS_4010 <- make_MP(SSS, HCR40_10, diagnostic = "min")