R/SS_recdevs.R

In r4ss/r4ss: R Code for Stock Synthesis

#' Insert a vector of recruitment deviations into the control file.
#'
#' A function to insert a vector of recruitment deviations into the control
#' file for simulation studies. This function was written in 2010, long before
#' the functions to read and write the input files were created. An alternative
#' approach would be to read the control file, add the recdevs, and then write
#' it again, but this function still works so there's no immediate need to
#' streamline that alternative approach.
#'
#' @param fyr First year of the recdev vector.
#' @param lyr Last year of the recdev vector.
#' @param ctl Either NULL to read anew or an already read control file.
#' Default=NULL.
#' @param recdevs Either NULL to generate anew or an already generated vector
#' of recdevs. Default=NULL.
#' @param rescale Should the recdevs be rescaled to have mean = 0 and std.
#' deviation = sigmaR? Default=TRUE.
#' @param scaleyrs Vector of years over which rescaling (if chosen) should
#' occur.
#' @template dir
#' @param ctlfile Name of control file to modify.  Default="control.ss_new".
#' @param newctlfile Name of new file to output modified control file.
#' Default="control_modified.ss".
#' @template verbose
#' @param writectl Write new file? Default=TRUE.
#' @param returnctl Return contents ctl file as an object in the R workspace.
#' Default=FALSE.
#' @param newmaxbias Replace the maximum bias adjustment fraction with any
#' non-NULL value. Default=NULL.
#' @author Ian Taylor
#' @export
SS_recdevs <-
  function(
    fyr,
    lyr,
    ctl = NULL,
    recdevs = NULL,
    rescale = TRUE,
    scaleyrs = NULL,
    dir = getwd(),
    ctlfile = "control.ss_new",
    newctlfile = "control_modified.ss",
    verbose = TRUE,
    writectl = TRUE,
    returnctl = FALSE,
    newmaxbias = NULL
  ) {
    # Determine working directory on start and return upon exit
    current_wd <- getwd()
    on.exit(setwd(current_wd))
    setwd(dir)

    # define a general function for reading values from control file
    readfun <- function(string, maxlen = Inf) {
      line1 <- grep(string, ctl)
      if (length(line1) < 1) {
        stop("no line contains the phrase, '", string, "'", sep = "")
      }
      if (length(line1) > 1) {
        stop("more than one line contains the phrase, '", string, "'", sep = "")
      }

      # split parameter line at hash mark
      splitline <- strsplit(ctl[line1], "#")[[1]]
      # split along all blank spaces, including tabs
      vecstrings <- strsplit(splitline[1], "[[:blank:]]+")[[1]]
      # convert to numeric
      vec <- as.numeric(vecstrings[vecstrings != ""])
      # check for length
      if (length(vec) > maxlen) {
        stop(paste(
          "this line has more than ",
          maxlen,
          " value",
          c("s", "")[1 + (maxlen == 1)],
          ": ",
          ctl[line1],
          sep = ""
        ))
      }
      return(vec)
    } # end readfun

    # read control file if ctl is not supplied
    if (is.null(ctl)) {
      ctl <- readLines(ctlfile)
    }

    # get sigma R
    sigmaR <- readfun("SR_sigmaR")[3]

    # make sure model includes recdevs and get some information
    do_recdev <- readfun("do_recdev", maxlen = 1)
    if (do_recdev == 0) {
      stop("do_recdev should be set to 1 or 2")
    }
    yrs <- fyr:lyr
    Nrecdevs <- lyr - fyr + 1
    phase <- readfun("recdev phase", maxlen = 1)
    advanced <- readfun("advanced options", maxlen = 1)
    if (advanced != 1) {
      stop("advanced options must be turned on in control file")
    }
    if (phase > 0) {
      newphase <- -abs(phase)
      if (verbose) {
        message("Changing recdev phase to negative: ", newphase)
      }
      ctl[grep("recdev phase", ctl)] <- paste(newphase, "#_recdev phase")
    }

    # turn on read_recdevs
    key1 <- grep("read_recdevs", ctl)
    ctl[key1] <- paste(Nrecdevs, "#_read_recdevs")

    # check for keyword at start of following section
    key2 <- grep("Fishing Mortality info", ctl)
    if (length(key2) == 0) {
      warning(
        "The phrase 'Fishing Mortality info' does not occur after the\n",
        "recdev section; Format of control file may be messy."
      )
    } else {
      key2 == key2[1]
    }

    # generate new recdevs
    if (!is.null(recdevs)) {
      if (length(recdevs) != Nrecdevs) {
        stop(paste(
          "input 'recdevs' has length=",
          length(recdevs),
          " but Nrecdevs=lyr-fyr+1=",
          Nrecdevs,
          sep = ""
        ))
      } else {
        newdevs <- recdevs
      }
    } else {
      newdevs <- rnorm(n = Nrecdevs)
    }
    if (rescale) {
      if (is.null(scaleyrs)) {
        scaleyrs <- yrs %in% yrs
      } else {
        scaleyrs <- yrs %in% scaleyrs
      }
      if (verbose) {
        message(
          "Rescaling recdevs vector so yrs ",
          min(yrs[scaleyrs]),
          ":",
          max(yrs[scaleyrs]),
          " have mean 0 and std. dev. = sigmaR = ",
          sigmaR
        )
      }
      newdevs <- sigmaR *
        (newdevs - mean(newdevs[scaleyrs])) /
        sd(newdevs[scaleyrs])
    }
    # build new recdev section
    newsection <- c(
      "#_end of advanced SR options",
      "#",
      "# read specified recr devs",
      "#_Yr Input_value"
    )
    # newsection <-c(newsection, rep("",(key2-key1-1)-length(newsection))) # preserve length of file

    for (i in 1:Nrecdevs) {
      newsection[4 + i] <-
        paste(
          (fyr:lyr)[i],
          c("  ", " ")[1 + (newdevs[i] < 0)],
          newdevs[i],
          " #_stochastic_recdev_with_sigmaR=",
          sigmaR,
          sep = ""
        )
    }

    ctl <- c(ctl[1:key1], newsection, ctl[key2:length(ctl)])
    # ctl[(key1+1):(key2-1)] <- newsection

    # if maxbias is input, then replace
    if (!is.null(newmaxbias)) {
      ctl[grep("max_bias", ctl)] <- paste(newmaxbias, "#_max_bias_adj_in_MPD")
    }

    # write and/or return the modified control file
    if (writectl) {
      writeLines(ctl, newctlfile)
      if (verbose) {
        message("Wrote new file: ", newctlfile)
      }
    }
    # reset working directory
    setwd(current_wd)
    if (returnctl) {
      return(ctl)
    }
  } # end function