R/z_MPs.R

In SAMtool: Stock Assessment Methods Toolkit

#' Make a custom management procedure (MP)
#'
#' @description
#' Function operator that creates a management procedure (MP) by combining an assessment model (function of class `Assess`) with
#' a harvest control rule (function of class `HCR`). The resulting function can then be tested in closed-loop simulation via 
#' [MSEtool::runMSE()]. 
#' 
#' \itemize{
#' \item Use `make_MP` to specify constant TAC between assessments; the frequency of
#' assessments is specified in `OM@@interval`. 
#' \item Use `make_projection_MP` to set catches according to a schedule set by projections,
#' specify assessment frequency in argument `assessment_interval` and ensure that `OM@@interval <- 1`. 
#' \item Use `make_interim_MP` to use an interim procedure to adjust the TAC between 
#' assessments using an index (Huynh et al. 2020), with the frequency of assessments specified in argument `assessment_interval` when 
#' making the MP; ensure that `OM@@interval <- 1`.
#' }
#' 
#' @param .Assess Assessment model, a function of class `Assess`.
#' @param .HCR Harvest control rule, a function of class `HCR`. Currently not used in projection MPs.
#' @param diagnostic A character string describing if any additional diagnostic information from the
#' assessment models will be collected during the closed-loop simulation. `"min"` (minimal) will collect information 
#' on convergence (default) and `"full"` will also collect the model estimates of biomass and F generated by `.Assess`.
#' `"none"` skips this step. 
#' @param ... Additional arguments to be passed to `.Assess` and `.HCR`.
#'
#' @return A function of class `MP`.
#'
#' @examples
#' # A statistical catch-at-age model with a 40-10 control rule
#' SCA_40_10 <- make_MP(SCA, HCR40_10)
#'
#' # An SCA that will produce convergence diagnostics
#' SCA_40_10 <- make_MP(SCA, HCR40_10, diagnostic = "min")
#'
#' # MP with an SCA that uses a Ricker stock-recruit function.
#' SCA_Ricker <- make_MP(SCA, HCR_MSY, SR = "Ricker")
#' show(SCA_Ricker)
#' 
#' \donttest{
#' # TAC is calculated annually from triennial assessments with projections between
#' # assessments with F = 0.75 FMSY
#' # Projections by default assume no process error.
#' OM <- MSEtool::testOM
#' OM@interval <- 1
#' pMP <- make_projection_MP(SCA, assessment_interval = 3, 
#'                           Ftarget = expression(0.75 * Assessment@FMSY),
#'                           proj_args = list(process_error = 1))
#' }
#'
#' @seealso [HCR_ramp] [HCR_MSY] [diagnostic] [retrospective_AM]
#' @export
make_MP <- function(.Assess, .HCR, diagnostic = c("min", "full", "none"), ...) {
  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_char <- .HCR
    .HCR <- as.symbol(.HCR)
  } else {
    .HCR <- substitute(.HCR)
    HCR_char <- as.character(.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)]
  dots_lost <- dots[!names(dots) %in% c(names(formals(eval(.Assess))), names(formals(eval(.HCR))))]
  
  #if (length(dots_in_Assess)) message_info("Arguments passed to ", Assess_char, ": ", names(dots_in_Assess) %>% paste0(collapse = ","))
  #if (length(dots_in_HCR)) message_info("Arguments passed to ", HCR_char, ": ", names(dots_in_HCR) %>% paste0(collapse = ","))
  if (length(dots_lost)) warning("Unused arguments: ", names(dots_lost) %>% paste0(collapse = ", "))

  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))

  MP_body <- bquote({
    dependencies <- .(get_dependencies(Assess_char, dots))
    do_Assessment <- .(Assess_call)
    Rec <- .(HCR_call)
    if (diagnostic != "none") Rec@Misc <- Assess_diagnostic(x, Data, do_Assessment, include_assessment = .(diagnostic == "full"))
    if (!is.null(do_Assessment@info$Misc)) Rec@Misc <- c(Rec@Misc, do_Assessment@info$Misc)
    return(Rec)
  })

  custom_MP <- eval(call("function", as.pairlist(alist(x = , Data = , reps = 1)), MP_body))
  formals(custom_MP)$diagnostic <- diagnostic
  return(structure(custom_MP, class = "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 [make_MP] and the available Assess and HCR objects with constant
#' TAC between assessment years. 
#'
#' @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.
#' @param diagnostic Character string describing the assessment diagnostic to save, see [make_MP].
#'
#' @examples
#' MSEtool::avail("MP", package = "SAMtool")
#'
#' \donttest{
#' myMSE <- MSEtool::runMSE(MSEtool::testOM, MPs = c("FMSYref", "SCA_4010"))
#' }
#' @return An object of class [Rec-class] which contains the management recommendation.
NULL





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",
                   SCA_RWM = "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@L95, Data@CAA, 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, Data@LFC, Data@LFS",
                   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",
                   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")


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) {
    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 [SCA].
#' @export
SCA_MSY <- make_MP(SCA, HCR_MSY)


#' @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)


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


#' @describeIn Model-based-MP A state-space delay difference model with a TAC recommendation based on fishing at FMSY,
#' and default arguments for configuring [DD_SS].
#' @export
DDSS_MSY <- make_MP(DD_SS, HCR_MSY)


#' @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)


#' @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)


#' @describeIn Model-based-MP A surplus production model with a TAC recommendation based on fishing at FMSY,
#' and default arguments for configuring [SP].
#' @export
SP_MSY <- make_MP(SP, HCR_MSY)


#' @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)


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


#' @describeIn Model-based-MP Simple stock synthesis (terminal depletion fixed to 0.4 in [SSS]) with a TAC recommendation based on fishing at FMSY.
#' @export
SSS_MSY <- make_MP(SSS, HCR_MSY)

#' @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)

#' @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)