R/assess_SCA.R

In SAMtool: Stock Assessment Methods Toolkit

#' Statistical catch-at-age (SCA) model
#'
#' A generic statistical catch-at-age model (single fleet, single season) that uses catch, index, and catch-at-age composition
#' data. `SCA` parameterizes R0 and steepness as leading productivity parameters in the assessment model. Recruitment is estimated
#' as deviations from the resulting stock-recruit relationship. In `SCA2`, the mean recruitment in the time series is estimated and
#' recruitment deviations around this mean are estimated as penalized parameters (`SR = "none"`, similar to Cadigan 2016). The standard deviation is set high
#' so that the recruitment is almost like free parameters. Unfished and MSY reference points are not estimated, it is recommended to use yield per recruit
#' or spawning potential ratio in harvest control rules. `SCA_Pope` is a variant of `SCA` that fixes the expected catch to the observed
#' catch, and Pope's approximation is used to calculate the annual exploitation rate (U; i.e., `catch_eq = "Pope"`).
#'
#' @aliases SCA2 SCA_Pope
#' @param x A position in the Data object (by default, equal to one for assessments).
#' @param Data An object of class Data
#' @param AddInd A vector of integers or character strings indicating the indices to be used in the model. Integers assign the index to
#' the corresponding index in Data@@AddInd, "B" (or 0) represents total biomass in Data@@Ind, "VB" represents vulnerable biomass in
#' Data@@VInd, and "SSB" represents spawning stock biomass in Data@@SpInd. Vulnerability to the survey is fixed in the model.
#' @param SR Stock-recruit function (either `"BH"` for Beverton-Holt, `"Ricker"`, or `"none"` for constant mean recruitment).
#' @param vulnerability Whether estimated vulnerability is `"logistic"` or `"dome"` (double-normal).
#' See details for parameterization.
#' @param catch_eq Whether to use the Baranov equation or Pope's approximation to calculate the predicted catch at age in the model.
#' @param CAA_dist Whether a multinomial or lognormal distribution is used for likelihood of the catch-at-age matrix. See details.
#' @param CAA_multiplier Numeric for data weighting of catch-at-age matrix if `CAA_hist = "multinomial"`. Otherwise ignored. See details.
#' @param rescale A multiplicative factor that rescales the catch in the assessment model, which
#' can improve convergence. By default, `"mean1"` scales the catch so that time series mean is 1, otherwise a numeric.
#' Output is re-converted back to original units.
#' @param max_age Integer, the maximum age (plus-group) in the model.
#' @param start Optional list of starting values. Entries can be expressions that are evaluated in the function. See details.
#' @param prior A named list for the parameters of any priors to be added to the model. See below.
#' @param fix_h Logical, whether to fix steepness to value in `Data@@steep` in the model for `SCA`. This only affects
#' calculation of reference points for `SCA2`.
#' @param fix_F_equilibrium Logical, whether the equilibrium fishing mortality prior to the first year of the model
#' is estimated. If `TRUE`, `F_equilibrium` is fixed to value provided in `start` (if provided),
#' otherwise, equal to zero (assumes unfished conditions).
#' @param fix_omega Logical, whether the standard deviation of the catch is fixed. If `TRUE`,
#' omega is fixed to value provided in `start` (if provided), otherwise, value based on `Data@@CV_Cat`.
#' @param fix_tau Logical, the standard deviation of the recruitment deviations is fixed. If `TRUE`,
#' tau is fixed to value provided in `start` (if provided), otherwise, value based on `Data@@sigmaR`.
#' @param LWT A named list (Index, CAA, Catch) of likelihood weights for the data components. For the index, a vector of length survey. For
#' CAL and Catch, a single value.
#' @param early_dev Numeric or character string describing the years for which recruitment deviations are estimated in `SCA`. By default,
#' equal to `"comp_onegen"`, where rec devs are estimated one full generation prior to the first year when catch-at-age (CAA) data are available.
#' With `"comp"`, rec devs are estimated starting in the first year with CAA. With `"all"`, rec devs start at the beginning of the model.
#' If numeric, the number of years after the first year of the model for which to start estimating rec devs. Use negative numbers for years prior to the first year.
#' @param late_dev Typically, a numeric for the number of most recent years in which recruitment deviations will
#' not be estimated in `SCA` (recruitment in these years will be based on the mean predicted by stock-recruit relationship).
#' By default, `"comp50"` uses the number of ages (smaller than the mode)
#' for which the catch-at-age matrix has less than half the abundance than that at the mode.
#' @param integrate Logical, whether the likelihood of the model integrates over the likelihood
#' of the recruitment deviations (thus, treating it as a random effects/state-space variable).
#' Otherwise, recruitment deviations are penalized parameters.
#' @param silent Logical, passed to [TMB::MakeADFun()], whether TMB
#' will print trace information during optimization. Used for diagnostics for model convergence.
#' @param opt_hess Logical, whether the hessian function will be passed to [stats::nlminb()] during optimization
#' (this generally reduces the number of iterations to convergence, but is memory and time intensive and does not guarantee an increase
#' in convergence rate). Ignored if `integrate = TRUE`.
#' @param n_restart The number of restarts (calls to [stats::nlminb()]) in the optimization procedure, so long as the model
#' hasn't converged. The optimization continues from the parameters from the previous (re)start.
#' @param control A named list of arguments for optimization to be passed to
#' [stats::nlminb()].
#' @param inner.control A named list of arguments for optimization of the random effects, which
#' is passed on to [TMB::newton()].
#' @param ... Other arguments to be passed.
#' 
#' @section Priors:
#' The following priors can be added as a named list, e.g., `prior = list(M = c(0.25, 0.15), h = c(0.7, 0.1)`. 
#' For each parameter below, provide a vector of values as described:
#' 
#' \itemize{
#' \item `R0` - A vector of length 3. The first value indicates the distribution of the prior: `1` for lognormal, `2` for uniform
#' on `log(R0)`, `3` for uniform on R0. If lognormal, the second and third values are the prior mean (in normal space) and SD (in log space).
#' Otherwise, the second and third values are the lower and upper bounds of the uniform distribution (values in normal space).
#' \item `h` - A vector of length 2 for the prior mean and SD, both in normal space. Beverton-Holt steepness uses a beta distribution, 
#' while Ricker steepness uses a normal distribution.
#' \item `M` - A vector of length 2 for the prior mean (in normal space) and SD (in log space). Lognormal prior.
#' \item `q` - A matrix for nsurvey rows and 2 columns. The first column is the prior mean (in normal space) and the second column 
#' for the SD (in log space). Use `NA` in rows corresponding to indices without priors.
#' }
#' See online documentation for more details.
#' 
#' @details
#' The basic data inputs are catch (by weight), index (by weight/biomass), and catch-at-age matrix (by numbers).
#' 
#' With `catch_eq = "Baranov"` (default in SCA and SCA2), annual F's are estimated parameters assuming continuous fishing over the year, while
#' an annual exploitation rate from pulse fishing in the middle of the year is estimated in `SCA_Pope` or `SCA(catch_eq = "Pope")`. 
#'
#' The annual sample sizes of the catch-at-age matrix is provided to the model (used in the likelihood for catch-at-age assuming
#' a multinomial distribution) and is manipulated via argument `CAA_multiplier`. This argument is
#' interpreted in two different ways depending on the value provided. If `CAA_multiplier > 1`, then this value will cap the annual sample sizes
#' to that number. If `CAA_multiplier <= 1`, then all the annual samples sizes will be re-scaled by that number, e.g. `CAA_multiplier = 0.1`
#' multiplies the sample size to 10% of the original number. By default, sample sizes are capped at 50.
#'
#' Alternatively, a lognormal distribution with inverse proportion variance can be used for the catch at age (Punt and Kennedy, 1994, as
#' cited by Maunder 2011).
#'
#' For `start` (optional), a named list of starting values of estimates can be provided for:
#' \itemize{
#' \item `R0` Unfished recruitment, except when `SR = "none"` where it is mean recruitment. 
#' By default, 150% `Data@@OM$R0[x]` is used as the start value in closed-loop simulation, and 400% of mean catch otherwise.
#' \item `h` Steepness. Otherwise, `Data@@steep[x]` is used, or 0.9 if empty.
#' \item `M` Natural mortality. Otherwise, `Data@@Mort[x]` is used.
#' \item `vul_par` Vulnerability parameters, see next paragraph.
#' \item `F` A vector of length `nyears` for year-specific fishing mortality.
#' \item `F_equilibrium` Equilibrium fishing mortality leading into first year of the model (to determine initial depletion). By default, 0.
#' \item `U_equilibrium` Same as F_equilibrium when `catch_eq = "Pope"`. By default, 0.
#' \item `omega` Lognormal SD of the catch (observation error) when `catch_eq = "Baranov"`. By default, `Data@@CV_Cat[x]`.
#' \item `tau` Lognormal SD of the recruitment deviations (process error). By default, `Data@@sigmaR[x]`.
#' }
#' 
#' Vulnerability can be specified to be either logistic or dome. If logistic, then the parameter
#' vector `vul_par` is of length 2:
#' \itemize{
#' \item `vul_par[1]` corresponds to `a_95`, the age of 95% vulnerability. `a_95` is a transformed parameter via logit transformation to constrain `a_95` to less than 75%
#' of the maximum age: `a_95 = 0.75 * max_age * plogis(x[1])`, where `x` is the estimated vector.
#' \item `vul_par[2]` corresponds to `a_50`, the age of 50% vulnerability. Estimated as an offset, i.e., `a_50 = a_95 - exp(x[2])`.
#' }
#' 
#' With dome vulnerability, a double Gaussian parameterization is used, where `vul_par`
#' is an estimated vector of length 4:
#' \itemize{
#' \item `vul_par[1]` corresponds to  `a_asc`, the first age of full vulnerability for the ascending limb. In the model, `a_asc` is estimated via logit transformation
#' to constrain `a_95` to less than 75% of the maximum age: `a_asc = 0.75 * maxage * plogis(x[1])`, where `x` is the estimated vector.
#' \item `vul_par[2]` corresponds to `a_50`, the age of 50% vulnerability for the ascending limb. Estimated as an offset, i.e.,
#' `a_50 = a_asc - exp(x[2])`.
#' \item `vul_par[3]` corresponds to `a_des`, the last age of full vulnerability (where the descending limb starts). Generated via logit transformation
#' to constrain between `a_asc` and `max_age`, i.e., `a_des = (max_age - a_asc) * plogis(x[3]) + a_asc`. By default, fixed to a small value so that the dome is effectively
#' a three-parameter function.
#' \item `vul_par[4]` corresponds to `vul_max`, the vulnerability at the maximum age. Estimated in logit space: `vul_max = plogis(x[4])`.
#' }
#' Vague priors of `vul_par[1] ~ N(0, sd = 3)`, `vul_par[2] ~ N(0, 3)`, `vul_par[3] ~ Beta(1.01, 1.01)` are used to aid convergence when parameters may not be well estimated,
#' for example, when vulnerability >> 0.5 for the youngest age class.
#' 
#' @section Online Documentation:
#' Model description and equations are available on the openMSE 
#' [website](https://openmse.com/features-assessment-models/2-sca/).
#' 
#' @references
#' Cadigan, N.G. 2016. A state-space stock assessment model for northern cod, including under-reported catches and
#' variable natural mortality rates. Canadian Journal of Fisheries and Aquatic Science 72:296-308.
#'
#' Maunder, M.N. 2011. Review and evaluation of likelihood functions for composition data in
#' stock-assessment models: Estimating the effective sample size. Fisheries Research 209:311-319.
#'
#' Punt, A.E. and Kennedy, R.B. 1997. Population modelling of Tasmanian rock lobster, Jasus edwardsii, resources. Marine and Freshwater
#' Research 48:967-980.
#'
#' @examples
#' res <- SCA(Data = MSEtool::SimulatedData)
#' res2 <- SCA2(Data = MSEtool::SimulatedData)
#' 
#' # Downweight the index
#' res3 <- SCA(Data = MSEtool::SimulatedData, LWT = list(Index = 0.1, CAA = 1))
#'
#' compare_models(res, res2)
#' @section Required Data:
#' \itemize{
#' \item `SCA`, `SCA_Pope`, and `SCA_Pope`: Cat, Ind, Mort, L50, L95, CAA, vbK, vbLinf, vbt0, wla, wlb, MaxAge
#' }
#' @section Optional Data:
#' \itemize{
#' \item `SCA`: Rec, steep, sigmaR, CV_Ind, CV_Cat
#' \item `SCA2`: Rec, steep, CV_Ind, CV_Cat
#' \item `SCA_Pope`: Rec, steep, sigmaR, CV_Ind
#' }
#' @author Q. Huynh
#' @return An object of class [Assessment-class].
#' @seealso [plot.Assessment] [summary.Assessment] [retrospective] [profile] [make_MP]
#' @export
SCA <- function(x = 1, Data, AddInd = "B", SR = c("BH", "Ricker", "none"), 
                vulnerability = c("logistic", "dome"), catch_eq = c("Baranov", "Pope"),
                CAA_dist = c("multinomial", "lognormal"),
                CAA_multiplier = 50, rescale = "mean1", max_age = Data@MaxAge,
                start = NULL, prior = list(), fix_h = TRUE, fix_F_equilibrium = TRUE, fix_omega = TRUE, fix_tau = TRUE,
                LWT = list(), early_dev = c("comp_onegen", "comp", "all"), late_dev = "comp50", integrate = FALSE,
                silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
                control = list(iter.max = 2e5, eval.max = 4e5), inner.control = list(), ...) {
  
  out <- SCA_(x = x, Data = Data, AddInd = AddInd, SR = SR, vulnerability = vulnerability, catch_eq = catch_eq, comp = "age",
              comp_dist = CAA_dist, comp_multiplier = CAA_multiplier, 
              rescale = rescale, max_age = max_age, start = start, prior = prior, fix_h = fix_h, 
              fix_F_equilibrium = fix_F_equilibrium, fix_omega = fix_omega, fix_tau = fix_tau, 
              LWT = LWT, early_dev = early_dev, late_dev = late_dev, integrate = integrate,
              silent = silent, opt_hess = opt_hess, n_restart = n_restart, control = control, 
              inner.control = inner.control, ...)
  return(out)
}
class(SCA) <- "Assess"

#' @rdname SCA
#' @param common_dev Typically, a numeric for the number of most recent years in which a common recruitment deviation will
#' be estimated (in `SCA2`, uninformative years will have a recruitment closer to the mean, which can be very misleading,
#' especially near the end of the time series). By default, `"comp50"` uses the number of ages (smaller than the mode)
#' for which the catch-at-age matrix has less than half the abundance than that at the mode.
#' @export
SCA2 <- function(x = 1, Data, AddInd = "B", vulnerability = c("logistic", "dome"), CAA_dist = c("multinomial", "lognormal"),
                 CAA_multiplier = 50, rescale = "mean1", max_age = Data@MaxAge, start = NULL, prior = list(),
                 fix_h = TRUE, fix_F_equilibrium = TRUE, fix_omega = TRUE, fix_tau = TRUE, LWT = list(),
                 common_dev = "comp50", integrate = FALSE, silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
                 control = list(iter.max = 2e5, eval.max = 4e5), inner.control = list(), ...) {
  
  out <- SCA_(x = x, Data = Data, AddInd = AddInd, SR = "none", vulnerability = vulnerability, catch_eq = "Baranov", comp = "age",
              comp_dist = CAA_dist, comp_multiplier = CAA_multiplier, rescale = rescale, max_age = max_age,
              start = start, prior = prior, fix_h = fix_h, fix_F_equilibrium = fix_F_equilibrium, fix_omega = fix_omega, fix_tau = fix_tau,
              LWT = LWT, early_dev = "all", late_dev = common_dev, integrate = integrate,
              silent = silent, opt_hess = opt_hess, n_restart = n_restart, control = control, 
              inner.control = inner.control, ...)
  return(out)
}
class(SCA2) <- "Assess"


#' @rdname SCA
#' @param fix_U_equilibrium Logical, same as `fix_F_equilibrium` for `SCA_Pope`.
#' @export
SCA_Pope <- function(x = 1, Data, AddInd = "B", SR = c("BH", "Ricker", "none"), vulnerability = c("logistic", "dome"), CAA_dist = c("multinomial", "lognormal"),
                     CAA_multiplier = 50, rescale = "mean1", max_age = Data@MaxAge, start = NULL, prior = list(),
                     fix_h = TRUE, fix_U_equilibrium = TRUE, fix_tau = TRUE, LWT = list(),
                     early_dev = c("comp_onegen", "comp", "all"), late_dev = "comp50", integrate = FALSE,
                     silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
                     control = list(iter.max = 2e5, eval.max = 4e5), inner.control = list(), ...) {
  
  out <- SCA_(x = x, Data = Data, AddInd = AddInd, SR = SR, vulnerability = vulnerability, catch_eq = "Pope", comp = "age",
              comp_dist = CAA_dist, comp_multiplier = CAA_multiplier, rescale = rescale, max_age = max_age,
              start = start, prior = prior, fix_h = fix_h, fix_F_equilibrium = fix_U_equilibrium, fix_omega = TRUE, fix_tau, 
              LWT = LWT, early_dev = early_dev, late_dev = late_dev, integrate = integrate,
              silent = silent, opt_hess = opt_hess, n_restart = n_restart, control = control, 
              inner.control = inner.control, ...)
  return(out)
}
class(SCA_Pope) <- "Assess"


#' @useDynLib SAMtool
SCA_ <- function(x = 1, Data, AddInd = "B", SR = c("BH", "Ricker", "none"), 
                 vulnerability = c("logistic", "dome"), catch_eq = c("Baranov", "Pope"), comp = c("age", "length"),
                 comp_dist = c("multinomial", "lognormal"), comp_multiplier = c(50, 50), rescale = "mean1", max_age = Data@MaxAge,
                 start = NULL, prior = list(), fix_h = TRUE, fix_F_equilibrium = TRUE, fix_omega = TRUE, fix_tau = TRUE,
                 tv_M = c("none", "walk", "DD"), M_bounds = NULL, refyear = 1,
                 LWT = list(), early_dev = c("comp_onegen", "comp", "all"), late_dev = "comp50", integrate = FALSE,
                 silent = TRUE, opt_hess = FALSE, n_restart = ifelse(opt_hess, 0, 1),
                 control = list(iter.max = 2e5, eval.max = 4e5), inner.control = list(), ...) {
  
  dependencies <- "Data@Cat, Data@Ind, Data@Mort, Data@L50, Data@L95, Data@CAA, Data@vbK, Data@vbLinf, Data@vbt0, Data@wla, Data@wlb, Data@MaxAge"

  dots <- list(...)
  start <- lapply(start, eval, envir = environment())
  
  max_age <- as.integer(min(max_age, Data@MaxAge))
  n_age <- max_age + 1
  vulnerability <- match.arg(vulnerability)
  catch_eq <- match.arg(catch_eq)
  comp <- match.arg(comp, several.ok = TRUE)
  comp_dist <- match.arg(comp_dist)
  if (length(comp_multiplier) == 1) comp_multiplier <- rep(comp_multiplier, 2)
  SR <- match.arg(SR)  
  tv_M <- match.arg(tv_M)
  
  if (is.character(early_dev)) early_dev <- match.arg(early_dev)
  if (is.numeric(early_dev)) stopifnot(early_dev < length(Data@Year))
  
  if (any(names(dots) == "yind")) {
    yind <- eval(dots$yind)
  } else {
    yind <- which(!is.na(Data@Cat[x, ]))[1]
    yind <- yind:length(Data@Cat[x, ])
  }
  Year <- Data@Year[yind]
  C_hist <- Data@Cat[x, yind]
  if (any(is.na(C_hist) | C_hist < 0)) warning("Error. Catch time series is not complete.")
  
  n_y <- length(C_hist)
  if (any(names(dots) == "M_at_age") && dots$M_at_age) {
    M <- Data@Misc$StockPars$M_ageArray[x, , n_y] * Data@Obs$Mbias[x]
    prior$M <- NULL
  } else {
    M <- rep(Data@Mort[x], n_age)
  }
  a <- Data@wla[x]
  b <- Data@wlb[x]
  Linf <- Data@vbLinf[x]
  K <- Data@vbK[x]
  t0 <- Data@vbt0[x]
  La <- Linf * (1 - exp(-K * (c(0:max_age) - t0)))
  SD_La <- La * Data@LenCV[x]
  Wa <- a * La ^ b
  A50 <- min(0.5 * max_age, iVB(t0, K, Linf, Data@L50[x]))
  A95 <- max(A50+0.5, iVB(t0, K, Linf, Data@L95[x]))
  mat_age <- c(0, 1/(1 + exp(-log(19) * (c(1:max_age) - A50)/(A95 - A50)))) # Age-0 is immature
  mat_age <- mat_age/max(mat_age)
  
  if (any(comp == "age")) {
    Data <- expand_comp_matrix(Data, "CAA") # Make sure dimensions of CAA match that in catch (nyears).
    CAA_hist <- Data@CAA[x, yind, 1:n_age]
    if (max_age < Data@MaxAge) CAA_hist[, n_age] <- rowSums(Data@CAA[x, yind, n_age:(Data@MaxAge+1)], na.rm = TRUE)
    
    if (all(is.na(CAA_hist))) warning("No age composition found in Data object", call. = FALSE)
  } else {
    CAA_hist <- matrix(0, n_y, n_age)
  }
  CAA_n_nominal <- rowSums(CAA_hist, na.rm = TRUE)
  if (comp_multiplier[1] <= 1) {
    CAA_n_rescale <- comp_multiplier[1] * CAA_n_nominal
  } else {
    CAA_n_rescale <- pmin(CAA_n_nominal, comp_multiplier[1])
  }
  
  if (any(comp == "length")) {
    CAL_hist <- Data@CAL[x, yind, ]
    CAL_sum <- colSums(CAL_hist, na.rm = TRUE) # Remove length bins with zeros for all years
    CAL_cdf <- cumsum(CAL_sum)
    CAL_ind <- which(CAL_sum > 0)[1]:which.max(CAL_cdf)[1]
    CAL_hist <- CAL_hist[, CAL_ind]
    if (all(is.na(CAL_hist))) warning("No length composition found in Data object", call. = FALSE)
    
    CAL_mids <- Data@CAL_mids[CAL_ind]
    CAL_bins <- Data@CAL_bins[CAL_ind]
    PLA <- generate_PLA(La, SD_La, CAL_bins, CAL_mids)
  } else {
    CAL_hist <- matrix(0, n_y, 1)
    PLA <- matrix(1, n_age, 1)
    CAL_bins <- CAL_mids <- 1
  }
  CAL_n_nominal <- rowSums(CAL_hist, na.rm = TRUE)
  if (comp_multiplier[2] <= 1) {
    CAL_n_rescale <- comp_multiplier[2] * CAL_n_nominal
  } else {
    CAL_n_rescale <- pmin(CAL_n_nominal, comp_multiplier[2])
  }
  
  if (early_dev == "all") {
    est_early_rec_dev <- rep(1, n_age - 1)
    est_rec_dev <- rep(1, n_y)
  } else if (early_dev == "comp") {
    est_early_rec_dev <- rep(0, n_age-1)
    if (any(comp == "age")) {
      est_rec_dev <- ifelse(1:n_y < which(CAA_n_nominal > 0)[1], 0, 1)
    } else {
      est_rec_dev <- ifelse(1:n_y < which(CAL_n_nominal > 0)[1], 0, 1)
    }
  } else if (early_dev == "comp_onegen") {
    if (any(comp == "age")) {
      istart <- which(CAA_n_nominal > 0)[1] - n_age
    } else {
      istart <- which(CAL_n_nominal > 0)[1] - n_age
    }
    if (istart < 0) {
      early_start <- n_age + istart
      est_early_rec_dev <- ifelse(2:n_age < early_start, 0, 1) %>% rev()
      est_rec_dev <- rep(1, n_y)
    } else {
      est_early_rec_dev <- rep(0, n_age - 1)
      est_rec_dev <- ifelse(1:n_y < istart, 0, 1)
    }
  } else if (is.numeric(early_dev)) {
    if (early_dev > 1) {
      est_early_rec_dev <- rep(0, n_age-1)
      est_rec_dev <- ifelse(1:n_y < early_dev, 0, 11)
    } else {
      istart <- early_dev - 1
      est_early_rec_dev <- rep(c(1, 0), c(istart, n_age - istart - 1))
      est_rec_dev <- rep(1, n_y)
    }
  }
  if (tv_M == "DD") est_early_rec_dev <- rep(0, n_age-1) # Temporary for now
  
  if (is.character(late_dev) && late_dev == "comp50") {
    if (any(comp == "age")) {
      comp_ldev <- colSums(CAA_hist, na.rm = TRUE)/max(colSums(CAA_hist, na.rm = TRUE))
    } else {
      comp_ldev <- colSums(CAL_hist, na.rm = TRUE)/max(colSums(CAL_hist, na.rm = TRUE))
    }
    comp_mode <- which.max(comp_ldev)[1]
    comp50_ind <- which(comp_ldev[1:comp_mode] <= 0.5)
    comp50_ind <- comp50_ind[length(comp50_ind)]
    if (all(comp != "age")) comp50_ind <- ceiling(LinInterp(La, 0:n_age, Data@CAL_mids[comp50_ind]))
    late_dev <- ifelse(is.na(comp50_ind), 0, comp50_ind)
  }
  if (is.numeric(late_dev) && late_dev > 0) {
    if (late_dev > length(est_rec_dev)) late_dev <- length(est_rec_dev)
    ind_late <- (length(est_rec_dev) - late_dev + 1):length(est_rec_dev)
    est_rec_dev[ind_late] <- ifelse(SR == "none", max(est_rec_dev[-ind_late]), 0)
  }
  
  if (rescale == "mean1") rescale <- 1/mean(C_hist)
  
  Ind <- lapply(AddInd, Assess_I_hist, Data = Data, x = x, yind = yind)
  I_hist <- vapply(Ind, getElement, numeric(n_y), "I_hist")
  if (is.null(I_hist) || all(is.na(I_hist))) stop("No indices found.", call. = FALSE)
  
  I_sd <- vapply(Ind, getElement, numeric(n_y), "I_sd") %>% pmax(0.05)
  I_units <- vapply(Ind, getElement, numeric(1), "I_units")
  
  I_vul <- vapply(AddInd, function(xx) {
    if (xx == "B") {
      return(rep(1, n_age))
    } else if (xx == "SSB") {
      return(mat_age)
    } else if (xx == "VB") {
      return(rep(0, n_age))
    } else {
      return(Data@AddIndV[x, suppressWarnings(as.numeric(xx)), 1:n_age])
    }
  }, numeric(n_age))
  nsurvey <- ncol(I_hist)
  
  if (!is.list(LWT)) {
    if (!is.null(LWT) && length(LWT) != nsurvey) stop("LWT needs to be a vector of length ", nsurvey)
    LWT <- list(Index = LWT)
    LWT$CAA <- LWT$CAL <- LWT$Catch <- 1
  } else {
    if (is.null(LWT$Index)) LWT$Index <- rep(1, nsurvey)
    if (is.null(LWT$CAA)) LWT$CAA <- 1
    if (is.null(LWT$CAL)) LWT$CAL <- 1
    if (is.null(LWT$Catch)) LWT$Catch <- 1 
  }
  
  # Generate priors
  prior <- make_prior(prior, nsurvey, ifelse(SR == "BH", 1, 2), msg = FALSE)
  
  # M_bounds
  if (is.null(M_bounds)) {
    if (tv_M == "none") {
      M_bounds <- c(0, 1e4)
    } else {
      M_bounds <- c(0.75, 1.25) * range(M)
    }
  }
  
  data <- list(model = "SCA", C_hist = C_hist, rescale = rescale, I_hist = I_hist,
               I_sd = I_sd, I_units = I_units, I_vul = I_vul, abs_I = rep(0, nsurvey), nsurvey = nsurvey, 
               CAA_hist = apply(CAA_hist, 1, tiny_comp) %>% t(), CAA_n = CAA_n_rescale, 
               CAL_hist = apply(CAL_hist, 1, tiny_comp) %>% t(), CAL_n = CAL_n_rescale,
               LWT = c(LWT$Index, LWT$CAA, LWT$CAL, LWT$Catch),
               n_y = n_y, n_age = n_age, n_bin = ncol(PLA), 
               M_data = 1,
               weight = Wa, PLA = PLA, mat = mat_age, vul_type = vulnerability,
               SR_type = SR, comp_dist = comp_dist, catch_eq = catch_eq,
               est_early_rec_dev = est_early_rec_dev, est_rec_dev = est_rec_dev, yindF = as.integer(0.5 * n_y),
               tv_M = tv_M, M_bounds = M_bounds, use_prior = prior$use_prior, prior_dist = prior$pr_matrix,
               sim_process_error = 0L)
  if (any(names(dots) == "M_at_age") && dots$M_at_age) data$M_data <- M
  if (data$n_bin == 1) data$CAL_hist <- t(data$CAL_hist)
  
  # Starting values
  params <- list()
  if (!is.null(start)) {
    if (!is.null(start$R0) && is.numeric(start$R0)) params$R0x <- log(start$R0[1] * rescale)
    if (!is.null(start$h) && is.numeric(start$h)) {
      if (SR == "BH") {
        h_start <- (start$h[1] - 0.2)/0.8
        params$transformed_h <- logit(h_start)
      } else if (SR == "Ricker") {
        params$transformed_h <- log(start$h[1] - 0.2)
      }
    }
    if (!is.null(start$M) && is.numeric(start$M)) params$log_M0 <- log(start$M)
    if (catch_eq == "Baranov" && !is.null(start$F_equilibrium) && is.numeric(start$F_equilibrium)) {
      params$F_equilibrium <- start$F_equilibrium
    }
    if (catch_eq == "Pope" && !is.null(start$U_equilibrium) && is.numeric(start$U_equilibrium)) {
      params$F_equilibrium <- start$U_equilibrium
    }
    if (!is.null(start$vul_par) && is.numeric(start$vul_par)) {
      if (start$vul_par[1] > 0.75 * max_age) stop("start$vul_par[1] needs to be less than 0.75 * Data@MaxAge (see help).")
      if (vulnerability == "logistic") {
        if (length(start$vul_par) < 2) stop("Two parameters needed for start$vul_par with logistic vulnerability (see help).")
        if (start$vul_par[1] <= start$vul_par[2]) stop("start$vul_par[1] needs to be greater than start$vul_par[2] (see help).")
        
        params$vul_par <- c(logit(start$vul_par[1]/max_age/0.75), log(start$vul_par[1] - start$vul_par[2]))
      }
      if (vulnerability == "dome") {
        if (length(start$vul_par) < 4) stop("Four parameters needed for start$vul_par with dome vulnerability (see help).")
        if (start$vul_par[1] <= start$vul_par[2]) stop("start$vul_par[1] needs to be greater than start$vul_par[2] (see help).")
        if (start$vul_par[3] <= start$vul_par[1] || start$vul_par[3] >= max_age) {
          stop("start$vul_par[3] needs to be between start$vul_par[1] and Data@MaxAge (see help).")
        }
        if (start$vul_par[4] <= 0 || start$vul_par[4] >= 1) stop("start$vul_par[4] needs to be between 0-1 (see help).")
        
        params$vul_par <- c(logit(start$vul_par[1]/max_age/0.75), log(start$vul_par[1] - start$vul_par[2]),
                            logit(1/(max_age - start$vul_par[1])), logit(start$vul_par[4]))
      }
    }
    if (!is.null(start$F) && is.numeric(start$F)) {
      Fstart <- numeric(n_y)
      Fstart_ind <- data$yindF + 1
      Fstart[Fstart_ind] <- log(start$F[Fstart_ind])
      Fstart[-Fstart_ind] <- log(start$F[-Fstart_ind]/Fstart[Fstart_ind])
      params$log_F_dev <- Fstart
    }
    
    if (!is.null(start$omega) && is.numeric(start$omega)) params$log_omega <- log(start$omega)
    if (!is.null(start[["tau"]]) && is.numeric(start[["tau"]])) params$log_tau <- log(start[["tau"]])
    if (!is.null(start[["tau_M"]]) && is.numeric(start[["tau_M"]])) params$log_tau_M <- log(start[["tau_M"]])
  }
  
  if (is.null(params$R0x)) {
    params$R0x <- ifelse(is.null(Data@OM$R0[x]), log(mean(data$C_hist)) + 4, log(1.5 * rescale * Data@OM$R0[x]))
  }
  if (is.null(params$transformed_h)) {
    h_start <- ifelse(!fix_h && is.na(Data@steep[x]), 0.9, Data@steep[x])
    if (SR == "BH") {
      h_start <- (h_start - 0.2)/0.8
      params$transformed_h <- logit(h_start)
    } else if (SR == "Ricker") {
      params$transformed_h <- log(h_start - 0.2)
    } else {
      params$transformed_h <- 0
    }
  }
  if (is.null(params$log_M0)) params$log_M0 <- log(M) %>% mean()
  if (is.null(params$logit_M_walk)) params$logit_M_walk <- rep(0, n_y)
  if (is.null(params$F_equilibrium)) params$F_equilibrium <- 0
  if (is.null(params$vul_par)) {
    if (any(comp == "age")) {
      comp_ldev <- colSums(CAA_hist, na.rm = TRUE)/max(colSums(CAA_hist, na.rm = TRUE))
    } else {
      comp_ldev <- colSums(CAL_hist, na.rm = TRUE)/max(colSums(CAL_hist, na.rm = TRUE))
    }
    comp_mode <- which.max(comp_ldev)[1]
    
    if (is.na(Data@LFC[x]) && is.na(Data@LFS[x]) || Data@LFC[x] > Linf || Data@LFS[x] > Linf) {
      
      if (all(comp == "length")) {
        comp_mode <- ceiling(LinInterp(La, 0:max_age, Data@CAL_mids[comp_mode]))
        if (!length(comp_mode)) comp_mode <- 1
      }
      if (vulnerability == "logistic") params$vul_par <- c(logit(comp_mode/max_age/0.75), log(1))
      if (vulnerability == "dome") {
        params$vul_par <- c(logit(comp_mode/max_age/0.75), log(1), logit(1/(max_age - comp_mode)), logit(0.5))
      }
    } else {
      A5 <- min(iVB(t0, K, Linf, Data@LFC[x]), comp_mode - 1)
      Afull <- min(iVB(t0, K, Linf, Data@LFS[x]), 0.5 * max_age)
      A5 <- min(A5, Afull - 0.5)
      A50_vul <- mean(c(A5, Afull))
      
      if (vulnerability == "logistic") params$vul_par <- c(logit(Afull/max_age/0.75), log(Afull - A50_vul))
      if (vulnerability == "dome") {
        params$vul_par <- c(logit(Afull/max_age/0.75), log(Afull - A50_vul), logit(0.1/(max_age - Afull)), logit(0.5))
      }
    }
  }
  if (is.na(params$vul_par[1])) params$vul_par[1] <- 1
  if (is.null(params$log_F_dev)) {
    Fstart <- numeric(n_y)
    Fstart[data$yindF + 1] <- log(0.75 * mean(M))
    params$log_F_dev <- Fstart
  }
  
  if (is.null(params$log_omega)) {
    sigmaC <- max(0.01, sdconv(1, Data@CV_Cat[x]), na.rm = TRUE)
    params$log_omega <- log(sigmaC)
  }
  if (is.null(params[["log_tau"]])) {
    tau_start <- ifelse(is.na(Data@sigmaR[x]), 0.6, Data@sigmaR[x])
    params$log_tau <- log(tau_start)
  }
  if (is.null(params[["log_tau_M"]])) params$log_tau_M <- log(0.05)
  
  params$log_early_rec_dev <- rep(0, n_age - 1)
  params$log_rec_dev <- rep(0, n_y)
  
  LH <- list(LAA = La, SD_LAA = SD_La, CAL_mids = CAL_mids, WAA = Wa, Linf = Linf, K = K, t0 = t0, a = a, b = b, A50 = A50, A95 = A95)
  info <- list(Year = Year, data = data, params = params, LH = LH, control = control,
               inner.control = inner.control)
  
  map <- list()
  if (catch_eq == "Baranov" && any(info$data$C_hist <= 0)) {
    ind <- info$data$C_hist <= 0
    info$params$log_F_dev[ind] <- -20
    map_logF <- length(params$log_F_dev)
    map_logF[ind] <- NA
    map_logF[!ind] <- 1:sum(!ind)
    map$log_F_dev <- factor(map_logF)
  } else if (catch_eq == "Pope") {
    map$log_F_dev <- factor(rep(NA, n_y))
  }
  if (fix_h && !prior$use_prior[2]) map$transformed_h <- factor(NA)
  if (!prior$use_prior[3]) map$log_M0 <- factor(NA)
  if (tv_M != "walk") map$logit_M_walk <- factor(rep(NA, n_y))
  if (fix_F_equilibrium) map$F_equilibrium <- factor(NA)
  if (fix_omega) map$log_omega <- factor(NA)
  if (fix_tau) map$log_tau <- factor(NA)
  map$log_tau_M <- factor(NA)
  if (any(!est_early_rec_dev)) map$log_early_rec_dev <- factor(ifelse(est_early_rec_dev, 1:sum(est_early_rec_dev), NA))
  if (any(!est_rec_dev)) map$log_rec_dev <- factor(ifelse(est_rec_dev, 1:sum(est_rec_dev), NA))
  if (vulnerability == "dome") map$vul_par <- factor(c(1, 2, NA, 3))
  
  random <- NULL
  if (integrate) random <- c("log_early_rec_dev", "log_rec_dev", "logit_M_walk")
  
  obj <- MakeADFun(data = info$data, parameters = info$params, hessian = TRUE,
                   map = map, random = random, DLL = "SAMtool", inner.control = inner.control, silent = silent)
  
  if (catch_eq == "Pope") {
    # Add starting values for rec-devs and increase R0 start value if U is too high (> 0.975)
    high_U <- try(obj$report(c(obj$par, obj$env$last.par[obj$env$random]))$penalty > 0, silent = TRUE)
    if (!is.character(high_U) && !is.na(high_U) && high_U) {
      Recruit <- try(Data@Rec[x, ], silent = TRUE)
      if (is.numeric(Recruit) && length(Recruit) == n_y && any(!is.na(Recruit))) {
        log_rec_dev <- log(Recruit/mean(Recruit, na.rm = TRUE))
        log_rec_dev[is.na(est_rec_dev) | is.na(log_rec_dev) | is.infinite(log_rec_dev)] <- 0
        info$params$log_rec_dev <- log_rec_dev
        
        obj <- MakeADFun(data = info$data, parameters = info$params, hessian = TRUE,
                         map = map, random = random, DLL = "SAMtool", inner.control = inner.control, silent = silent)
      }
      while (obj$par["R0x"] < 30 && obj$report(c(obj$par, obj$env$last.par[obj$env$random]))$penalty > 0) {
        obj$par["R0x"] <- obj$par["R0x"] + 1
      }
    }
  }
  
  mod <- optimize_TMB_model(obj, control, opt_hess, n_restart)
  opt <- mod[[1]]
  SD <- mod[[2]]
  report <- obj$report(obj$env$last.par.best)
  
  Yearplusone <- c(Year, max(Year) + 1)
  YearEarly <- (Year[1] - n_age + 1):(Year[1] - 1)
  YearDev <- c(YearEarly, Year)
  YearR <- c(YearDev, max(YearDev) + 1)
  R <- c(rev(report$R_early), report$R)
  
  Dev <- structure(c(rev(report$log_early_rec_dev), report$log_rec_dev), names = YearDev)
  report$dynamic_SSB0 <- SCA_dynamic_SSB0(obj) %>% 
    structure(names = Yearplusone)
  
  nll_report <- ifelse(is.character(opt), ifelse(integrate, NA, report$nll), opt$objective)
  Assessment <- new("Assessment", Model = "SCA", 
                    Name = Data@Name, conv = SD$pdHess,
                    B0 = report$B0, R0 = report$R0, N0 = report$N0,
                    SSB0 = report$E0, VB0 = report$VB0,
                    B = structure(report$B, names = Yearplusone),
                    B_B0 = structure(report$B/report$B0, names = Yearplusone),
                    SSB = structure(report$E, names = Yearplusone),
                    SSB_SSB0 = structure(report$E/report$E0, names = Yearplusone),
                    VB = structure(report$VB, names = Yearplusone),
                    VB_VB0 = structure(report$VB/report$VB0, names = Yearplusone),
                    R = structure(R, names = YearR),
                    N = structure(rowSums(report$N), names = Yearplusone),
                    N_at_age = report$N,
                    Selectivity = matrix(report$vul, nrow = length(Year), ncol = n_age, byrow = TRUE),
                    Obs_Catch = structure(C_hist, names = Year),
                    Obs_Index = structure(I_hist, dimnames = list(Year, paste0("Index_", 1:nsurvey))),
                    Obs_C_at_age = CAA_hist,
                    Catch = structure(report$Cpred, names = Year),
                    Index = structure(report$Ipred, dimnames = list(Year, paste0("Index_", 1:nsurvey))),
                    C_at_age = report$CAApred,
                    Dev = Dev, Dev_type = "log-Recruitment deviations",
                    NLL = structure(c(nll_report, report$nll_comp, report$prior, report$penalty),
                                    names = c("Total", paste0("Index_", 1:nsurvey), "CAA", "CAL", "Catch", "Dev", "M_dev", "Prior", "Penalty")),
                    info = info, obj = obj, opt = opt, SD = SD, TMB_report = report,
                    dependencies = dependencies)
  if (tv_M != "walk") Assessment@NLL <- Assessment@NLL[names(Assessment@NLL) != "M_dev"]
  if (all(comp != "age")) Assessment@NLL <- Assessment@NLL[names(Assessment@NLL) != "CAA"]
  if (all(comp != "length")) Assessment@NLL <- Assessment@NLL[names(Assessment@NLL) != "CAL"]
  if (catch_eq == "Pope") Assessment@NLL <- Assessment@NLL[names(Assessment@NLL) != "Catch"]
  if (SR != "none") Assessment@h <- report$h
  if (catch_eq == "Baranov") {
    Assessment@FMort <- structure(report$F, names = Year)
  } else {
    Assessment@U <- structure(report$U, names = Year)
  }
  
  if (Assessment@conv) {
    SE_Early <- as.list(SD, "Std. Error")$log_early_rec_dev %>% rev()
    SE_Main <- as.list(SD, "Std. Error")$log_rec_dev
    SE_Dev <- structure(c(SE_Early, SE_Main), names = YearDev)
    if (any(is.na(SE_Dev))) {
      Dev <- Dev[seq(which(!is.na(SE_Dev))[1], length(YearDev))]
      SE_Dev <- SE_Dev[seq(which(!is.na(SE_Dev))[1], length(YearDev))]
      SE_Dev[is.na(SE_Dev)] <- 0
    }
    
    refyear <- eval(refyear)
    
    ref_pt <- ref_pt_SCA(y = refyear, obj = obj, report = report)
    if (catch_eq == "Baranov") {
      report$FMSY <- ref_pt$FMSY
    } else {
      report$UMSY <- ref_pt$UMSY
    }
    report$MSY <- ref_pt$MSY
    report$VBMSY <- ref_pt$VBMSY
    report$RMSY <- ref_pt$RMSY
    report$BMSY <- ref_pt$BMSY
    report$EMSY <- ref_pt$EMSY
    report$refyear <- refyear
    
    if (!all(refyear == 1)) { # New reference points based on change in M
      report$new_B0 <- Assessment@B0 <- ref_pt$new_B0
      report$new_E0 <- Assessment@SSB0 <- ref_pt$new_E0
      report$new_VB0 <- Assessment@VB0 <- ref_pt$new_VB0
      report$new_R0 <- Assessment@R0 <- ref_pt$new_R0
      report$new_h <- Assessment@h <- ref_pt$new_h
      
      Assessment@B_B0 <- Assessment@B/Assessment@B0
      Assessment@SSB_SSB0 <- Assessment@SSB/Assessment@SSB0
      Assessment@VB_VB0 <- Assessment@VB/Assessment@VB0
    }
    
    if (catch_eq == "Baranov") {
      Assessment@FMSY <- report$FMSY
      Assessment@F_FMSY <- structure(report$F/Assessment@FMSY, names = Year)
    } else {
      Assessment@UMSY <- report$UMSY
      Assessment@U_UMSY <- structure(report$U/Assessment@UMSY, names = Year)
    }
    Assessment@MSY <- report$MSY
    Assessment@BMSY <- report$BMSY
    Assessment@SSBMSY <- report$EMSY
    Assessment@VBMSY <- report$VBMSY
    Assessment@B_BMSY <- structure(report$B/Assessment@BMSY, names = Yearplusone)
    Assessment@SSB_SSBMSY <- structure(report$E/Assessment@SSBMSY, names = Yearplusone)
    Assessment@VB_VBMSY <- structure(report$VB/Assessment@VBMSY, names = Yearplusone)
    Assessment@Dev <- Dev
    Assessment@SE_Dev <- SE_Dev
    Assessment@TMB_report <- report
    
    catch_eq_fn <- function(Ftarget) {
      projection_SCA(Assessment, Ftarget = Ftarget, p_years = 1, p_sim = 1, 
                     obs_error = list(array(1, c(1, 1, nsurvey)), matrix(1, 1, 1)),
                     process_error = matrix(1, 1, 1)) %>% slot("Catch") %>% as.vector()
    }
    Assessment@forecast <- list(per_recruit = ref_pt[["per_recruit"]], catch_eq = catch_eq_fn)
  }
  return(Assessment)
}

ref_pt_SCA <- function(y = 1, obj, report) {

  if (obj$env$data$SR_type == "none") {
    # Fit BH 
    R0_start <- log(mean(report$R))
    h_start <- logit((0.7 - 0.2)/0.8)
    opt <- nlminb(c(R0_start, h_start), get_SR, E = report$E, R = report$R, EPR0 = report$EPR0)
    SR_par <- get_SR(opt$par, E = report$E, R = report$R, EPR0 = report$EPR0, opt = FALSE)
    
    Arec <- SR_par$Arec
    Brec <- SR_par$Brec
    SR <- "BH"
  } else {
    SR_par <- NULL
    SR <- obj$env$data$SR_type
    Arec <- report$Arec
    Brec <- report$Brec
  }
  
  M <- apply(report$M[y, , drop = FALSE], 2, mean)
  weight <- obj$env$data$weight
  mat <- obj$env$data$mat
  vul <- report$vul
  
  catch_eq <- obj$env$data$catch_eq
  
  tv_M <- obj$env$data$tv_M
  M_bounds <- obj$env$data$M_bounds
  B0 <- report$B0
  
  max_F <- ifelse(catch_eq == "Baranov", 4, 0.99)
  opt2 <- optimize(yield_fn_SCA, interval = c(1e-4, max_F), M = M, mat = mat, weight = weight, vul = vul, 
                   SR = SR, Arec = Arec, Brec = Brec, catch_eq = catch_eq, B0 = B0, tv_M = tv_M, M_bounds = M_bounds)
  opt3 <- yield_fn_SCA(opt2$minimum, M = M, mat = mat, weight = weight, vul = vul, SR = SR, 
                       Arec = Arec, Brec = Brec, opt = FALSE, catch_eq = catch_eq, B0 = B0, tv_M = tv_M, M_bounds = M_bounds)
  
  if (catch_eq == "Baranov") {
    FMSY <- opt2$minimum
  } else {
    UMSY <- opt2$minimum
  }
  MSY <- -1 * opt2$objective
  VBMSY <- opt3["VB"]
  RMSY <- opt3["R"]
  BMSY <- opt3["B"]
  EMSY <- opt3["E"]
  
  if (catch_eq == "Baranov") {
    Fvec <- seq(0, 2.5 * FMSY, length.out = 100)
  } else {
    Fvec <- seq(0, 0.99, 0.01)
  }
  
  yield <- lapply(Fvec,
                  yield_fn_SCA, M = M, mat = mat, weight = weight, vul = vul, SR = SR, 
                  Arec = Arec, Brec = Brec, opt = FALSE, catch_eq = catch_eq, B0 = B0, tv_M = tv_M, M_bounds = M_bounds)
  EPR <- vapply(yield, getElement, numeric(1), "EPR")
  YPR <- vapply(yield, getElement, numeric(1), "YPR")
  
  new_B0 <- yield[[1]]["B"] # New due to change in M
  new_E0 <- yield[[1]]["E"]
  new_VB0 <- yield[[1]]["VB"]
  new_R0 <- yield[[1]]["R"]
  if (SR == "BH") {
    new_h <- Arec * EPR[1]/ (4 + Arec * EPR[1])
  } else {
    new_h <- 0.2 * (Arec * EPR[1])^0.8
  }
  
  if (catch_eq == "Baranov") {
    return(list(FMSY = FMSY, MSY = MSY, VBMSY = VBMSY, RMSY = RMSY, BMSY = BMSY, EMSY = EMSY,
                per_recruit = data.frame(FM = Fvec, SPR = EPR/EPR[1], YPR = YPR), SR_par = SR_par,
                new_B0 = new_B0, new_E0 = new_B0, new_VB0 = new_VB0, new_R0 = new_R0, new_h = new_h))
  } else {
    return(list(UMSY = UMSY, MSY = MSY, VBMSY = VBMSY, RMSY = RMSY, BMSY = BMSY, EMSY = EMSY,
                per_recruit = data.frame(U = Fvec, SPR = EPR/EPR[1], YPR = YPR), SR_par = SR_par,
                new_B0 = new_B0, new_E0 = new_B0, new_VB0 = new_VB0, new_R0 = new_R0, new_h = new_h))
  }
  
}

yield_fn_SCA <- function(x, M, mat, weight, fec = mat * weight, vul, SR = c("BH", "Ricker"), Arec, Brec, 
                         catch_eq = c("Baranov", "Pope"), opt = TRUE, x_transform = FALSE, B0 = 1,
                         tv_M = c("none", "walk", "DD"), M_bounds = NULL, spawn_time_frac = 0) {
  if (is.null(tv_M)) tv_M <- "none"
  tv_M <- match.arg(tv_M)
  
  if (tv_M != "DD") {
    yield_fn_SCA_int(x, M, mat, weight, fec, vul, SR, Arec, Brec, catch_eq, opt, x_transform, spawn_time_frac)
  } else {
    
    dep <- M_DD <- numeric(21)
    dep[1] <- 0.4
    for(i in 1:20) {
      M_DD[i] <- ifelse(dep[i] >= 1, M_bounds[1], 
                        ifelse(dep[i] <= 0, M_bounds[2], M_bounds[1] + (M_bounds[2] - M_bounds[1]) * (1 - dep[i])))
      out <- yield_fn_SCA_int(x, M = rep(M_DD[i], length(mat)), mat, weight, fec, vul, SR, Arec, Brec, catch_eq, 
                              opt = FALSE, x_transform = x_transform, spawn_time_frac = spawn_time_frac)
      if (abs(out["B"]/B0 - dep[i]) <= 1e-4) break
      dep[i+1] <- out["B"]/B0
    }
    
    if (opt) {
      return(-1 * out["Yield"])
    } else {
      return(out)
    }
  }
}

yield_fn_SCA_int <- function(x, M, mat, weight, fec = mat * weight, vul, SR = c("BH", "Ricker"), Arec, Brec, 
                             catch_eq = c("Baranov", "Pope"), opt = TRUE, x_transform = FALSE,
                             spawn_time_frac = 0) {
  SR <- match.arg(SR)
  catch_eq <- match.arg(catch_eq)
  if (catch_eq == "Baranov") {
    FMort <- ifelse(x_transform, exp(x), x)
    Z <- vul * FMort + M
    surv <- exp(-Z)
    spawn_surv <- exp(-spawn_time_frac * Z)
  } else {
    U <- ifelse(x_transform, ilogit(x), x)
    surv <- exp(-M) * (1 - vul * U)
    spawn_surv <- exp(-spawn_time_frac * M)
    if (spawn_time_frac > 0.5) spawn_surv <- spawn_surv * (1 - vul * U) # If tied, spawning goes first
  }
  n_age <- length(M)
  NPR <- calc_NPR(surv, n_age)
  EPR <- sum(NPR * spawn_surv * fec)
  if (SR == "BH") {
    Req <- (Arec * EPR - 1)/(Brec * EPR)
  } else if (SR == "Ricker") {
    Req <- log(Arec * EPR)/(Brec * EPR)
  }
  
  if (catch_eq == "Baranov") {
    CPR <- Baranov(vul, FMort, M, NPR)
  } else {
    CPR <- vul * U * NPR * exp(-0.5 * M)
  }
  YPR <- sum(CPR * weight)
  Yield <- YPR * Req
  if (opt) {
    return(-1 * Yield)
  } else {
    
    B <- Req * sum(NPR * weight)
    E <- Req * EPR
    VB <- Req * sum(NPR * vul * weight)
    
    return(c(EPR = EPR, Yield = Yield, YPR = YPR, B = B, E = E, VB = VB, R = Req))
  }
}

SCA_dynamic_SSB0 <- function(obj, par = obj$env$last.par.best, ...) {
  
  newdata <- obj$env$data
  if (obj$env$data$catch_eq == "Pope") {
    newdata$C_hist <- rep(1e-8, newdata$n_y)
    par[names(par) == "F_equilibrium"] <- 0
    
    obj2 <- MakeADFun(data = newdata, parameters = clean_tmb_parameters(obj), map = obj$env$map, 
                      random = obj$env$random, DLL = "SAMtool", silent = TRUE)
    out <- obj2$report(par)$E
  } else {
    par[names(par) == "log_F_dev"] <- log(1e-8)
    par[names(par) == "F_equilibrium"] <- 0
    out <- obj$report(par)$E
  }
  return(out)
}

get_SR <- function(pars, E, R, EPR0, opt = TRUE, figure = FALSE, type = c("BH", "Ricker"), fix_h = FALSE, h = NULL) {
  type <- match.arg(type)
  
  R0 <- exp(pars[1])
  E0 <- R0 * EPR0
  if (type == "BH") {
    if (!fix_h) h <- 0.2 + 0.8 * ilogit(pars[2])
    Arec <- 4*h/(1-h)/EPR0
    Brec <- (5*h-1)/(1-h)/E0
    
    Rpred <- Arec * E / (1 + Brec * E)
  } else if (type == "Ricker") {
    if (!fix_h) h <- 0.2 + exp(pars[2])
    Arec <- 1/EPR0 * (5*h)^1.25
    Brec <- 1.25 * log(5*h) / E0
    
    Rpred <- Arec * E * exp(-Brec * E)
  }
  sigmaR <- sqrt(sum((log(R/Rpred))^2)/length(R))
  
  if (opt){
    return(-sum(dnorm(log(R/Rpred), 0, sigmaR, log = TRUE)))
  } else {
    
    if (figure) {
      plot(E, R, ylim = c(0, max(R, R0)), xlim = c(0, max(E, E0)), xlab = "SSB", ylab = "Recruitment")
      
      E2 <- seq(0, E0, length.out = 500)
      if (type == "BH") Rpred2 <- Arec * E2 / (1 + Brec * E2)
      if (type == "Ricker") Rpred2 <- Arec * E2 * exp(-Brec * E2)
      
      lines(E2, Rpred2, col = "blue")
      abline(v = c(0.2 * E0, E0), h = c(h * R0, R0), lty = 2, col = "red")
      legend("topright", legend = c(paste0("h = ", round(h, 3)), paste0("log(R0) = ", round(log(R0), 3))), bty = "n")
    }
    
    return(list(Rpred = Rpred, R0 = R0, h = h, E0 = E0, sigmaR = sigmaR, Arec = Arec, Brec = Brec))
  }
}

generate_PLA <- function(LAA, SD_LAA, len_bins, len_mids) {
  n_age <- length(LAA)
  n_bin <- length(len_mids)
  
  PLA <- matrix(0, n_age, n_bin)
  PLA[, n_bin] <- 1 - pnorm(len_bins[n_bin], LAA, SD_LAA)
  PLA[, 1] <- pnorm(len_bins[2], LAA, SD_LAA)
  PLA[, 3:n_bin - 1] <- vapply(1:n_age, function(x) {
    pnorm(len_bins[3:n_bin], LAA[x], SD_LAA[x]) - pnorm(len_bins[3:n_bin - 1], LAA[x], SD_LAA[x])
  }, numeric(n_bin - 2)) %>% t()
  return(PLA)
}