R/solve_numbers.R

In nwfsc-assess/PEPtools: Tools for Stock Assessments used by the Population Ecology Program

#' Function to solve out numbers of fish in the forecast to equal
#' match management ABC values. If a fleet in your model has removals
#' specified as numbers of fish (1,000s), SS will expect numbers of fish
#' to be used in the forecast file for that fleet is there are pre-specified
#' forecast removals.  This creates challenges when attempting to determine the
#' right number of fish that will meet the harvest specified ABC values. This
#' function does a bisection approach to iteratively determine the number of fish
#' that will match the fleet specified removals in terms of biomass.
#'
#' @param mod_dir the directory of your model - all runs will be conducted in this folder
#' make sure you are alright if the SS files are changed.
#' @param fore_yrs a vector of forecast years where removals are pre-specified
#' @param fleet_abc a vector of fleet specific abc values for the fleet that is currently
#' in terms of numbers of fish.
#' @param fleet fleet number within SS3 of the fleet that has removals in terms of numbers of
#' fish.
#' @param thresh Percent difference that controls if the model should be run again. The
#' default value of 0.01 results in the code identifying a solution as correct that is within
#' 1 percent above or below the input fleet_abc value. The model will only be rerun once. This
#' threshold may need to be increased if the catch values are very small.
#' @param exe The executable name for Stock Synthesis to be passed to the r4ss::run function.
#' @author Chantel Wetzel
#' @export
#'
#' @examples
#' \dontrun{
#' solve_numbers(
#'   mod_dir = "C:/Models/my_model",
#'   fore_yrs = 2021:2022,
#'   fleet_abc = c(5.5, 5),
#'   fleet = 4
#' )
#' }
#'
solve_numbers <- function(mod_dir, fore_yrs, fleet_abc, fleet = NULL, exe = "ss", thresh = 0.01) {
  # Set the wd to run the model
  setwd(mod_dir)

  if (!file.exists("ss.par")) {
    stop("There is no par file in mod_dir. Please run the model.")
  }

  if (length(fore_yrs) != length(fleet_abc)) {
    stop("The length of fore_yrs and fleet_abc need to be the same.")
  }

  yrs <- fore_yrs
  abc <- fleet_abc

  if (is.null(fleet)) {
    dat <- r4ss::SS_readdat("data.ss_new")
    fleet <- which(dat$fleetinfo$units == 2)

    if (length(fleet) > 1) {
      stop("There appears to be more than one fleet with catches in numbers. \n
			  Function currently only does one fleet at a time.  \n
			  Specify the fleet input in the function to specify which fleet to do.")
    }
  }

  # Turn max phase to 0 and read from par file
  starter <- r4ss::SS_readstarter(file.path(mod_dir, "starter.ss"))
  use_par <- starter$init_values_src
  starter$init_values_src <- 1
  max_phase <- starter$last_estimation_phase
  starter$last_estimation_phase <- 0
  r4ss::SS_writestarter(starter, dir = mod_dir, overwrite = TRUE, verbose = FALSE)
  r4ss::run(exe = exe, extras = "-nohess -maxfun 0", skipfinished = FALSE)

  for (i in 1:length(yrs)) {
    fore <- r4ss::SS_readforecast(file.path(mod_dir, "forecast.ss"), verbose = FALSE)
    find <- which(fore$ForeCatch$Fleet == fleet & fore$ForeCatch$Year == yrs[i])
    fore$ForeCatch[find, "Catch or F"] <- abc[i]
    r4ss::SS_writeforecast(fore, dir = mod_dir, overwrite = TRUE)
    r4ss::run(exe = exe, extras = "-nohess -maxfun 0", skipfinished = FALSE)

    rep <- r4ss::SS_output(mod_dir, printstats = FALSE, verbose = FALSE, covar = FALSE)
    bio <- rep$timeseries[rep$timeseries$Yr == yrs[i], paste0("dead(B):_", fleet)]
    num <- rep$timeseries[rep$timeseries$Yr == yrs[i], paste0("dead(N):_", fleet)]
    wt <- bio / num
    input_num <- abc[i] / wt

    fore <- r4ss::SS_readforecast(file.path(mod_dir, "forecast.ss"), verbose = FALSE)
    find <- which(fore$ForeCatch$Fleet == fleet & fore$ForeCatch$Year == yrs[i])

    # Write out the solved input numbers based on the average weight and the abc
    fore$ForeCatch[find, "Catch or F"] <- input_num
    r4ss::SS_writeforecast(fore, dir = mod_dir, overwrite = TRUE)
    # Rerun the model
    r4ss::run(exe = exe, extras = "-nohess -maxfun 0", skipfinished = FALSE)

    rep <- r4ss::SS_output(mod_dir, covar = FALSE, verbose = FALSE, printstats = FALSE)
    bio <- rep$timeseries[rep$timeseries$Yr == yrs[i], paste0("dead(B):_", fleet)]

    print(paste0("!!!!!!!!!!!!!!!!! Biomass = ", bio, " ABC = ", abc[i], "!!!!!!!!!!!!!!!!"))
    if (bio / abc[i] < (1 - thresh) | bio / abc[i] > (1 + thres)) {
      bio <- rep$timeseries[rep$timeseries$Yr == yrs[i], paste0("dead(B):_", fleet)]
      num <- rep$timeseries[rep$timeseries$Yr == yrs[i], paste0("dead(N):_", fleet)]
      wt <- bio / num
      input_num <- abc[i] / wt

      fore <- r4ss::SS_readforecast(file.path(mod_dir, "forecast.ss"), verbose = FALSE)
      find <- which(fore$ForeCatch$Fleet == fleet & fore$ForeCatch$Year == yrs[i])

      # Write out the solved input numbers based on the average weight and the abc
      fore$ForeCatch[find, "Catch or F"] <- input_num
      r4ss::SS_writeforecast(fore, dir = mod_dir, overwrite = TRUE)
      # Rerun the model
      r4ss::run(exe = exe, extras = "-nohess -maxfun 0", skipfinished = FALSE)
      rep <- r4ss::SS_output(mod_dir, covar = FALSE, verbose = FALSE, printstats = FALSE)
      bio <- rep$timeseries[rep$timeseries$Yr == yrs[i], paste0("dead(B):_", fleet)]

      print(paste0("!!!!!!!!!!!!!!!!! Biomass = ", bio, " ABC = ", abc[i], "!!!!!!!!!!!!!!!!"))
    }
  }

  # Change the max phase back to the original value
  starter <- r4ss::SS_readstarter(file.path(mod_dir, "starter.ss"))
  starter$init_values_src <- use_par
  starter$last_estimation_phase <- max_phase
  r4ss::SS_writestarter(starter, dir = mod_dir, overwrite = TRUE, verbose = FALSE)
}