R/MPs_Output.R

Defines functions YPR_ML YPR_CC YPR YPRopt YPR_ SPSRA_ML SPSRA SPSRAopt SPSRA_ SPMSY SPslope SPmod SBT2 SBT1 Rcontrol2 Rcontrol MCD4010 MCD LstepCC1 LstepCC_ Lratio_BHI3 Lratio_BHI2 Lratio_BHI L95target Ltarget1 Ltarget_ ITM Itarget1 Itarget_ IT5 IT_ Islope1 Islope_ Iratio ICI2 ICI Gcontrol GB_target GB_slope Fratio_ML Fratio_CC DepF Fratio4010 Fratio Fratio_ Fdem_ML Fdem_CC Fdem Fdem_ demofn demographic2 getr Fadapt DynF DD4010 DD DD_R DD_ DBSRA4010 DBSRA_40 DBSRA DBSRA_ DBSRAopt HDAAC DAAC DCAC_ML DCAC4010 DCAC_40 DCAC DCACs DCAC_ CompSRA4010 CompSRA CompSRA_ GB_CC CC1 BK_ML BK_CC BK AvC

Documented in AvC BK BK_CC BK_ML CC1 CompSRA CompSRA_ CompSRA4010 DAAC DBSRA DBSRA_ DBSRA_40 DBSRA4010 DBSRAopt DCAC DCAC_ DCAC_40 DCAC4010 DCAC_ML DCACs DD DD_ DD4010 DD_R demofn demographic2 DepF DynF Fadapt Fdem Fdem_ Fdem_CC Fdem_ML Fratio Fratio_ Fratio4010 Fratio_CC Fratio_ML GB_CC GB_slope GB_target Gcontrol getr HDAAC ICI ICI2 Iratio Islope_ Islope1 IT_ IT5 Itarget_ Itarget1 ITM L95target Lratio_BHI Lratio_BHI2 Lratio_BHI3 LstepCC1 Ltarget1 MCD MCD4010 Rcontrol Rcontrol2 SBT1 SBT2 SPmod SPMSY SPslope SPSRA SPSRA_ SPSRA_ML SPSRAopt YPR YPR_ YPR_CC YPR_ML YPRopt

## Output Control MPs ####



# ----- Average Catch ----

#' Average Catch
#'
#' A simple average catch MP that is included to demonstrate a 'status quo' management option
#'
#' The average catch method is very simple. The mean historical catch is calculated and used to set
#' a constant catch limit (TAC). If `reps` > 1 then the `reps` samples are drawn from a log-normal
#' distribution with mean `TAC` and standard deviation (in log-space) of 0.2.
#'
#' For completeness, the TAC is calculated by:
#'
#' \deqn{\textrm{TAC} =\frac{\sum_{y=1}^{\textrm{n}}{C_y}}{\textrm{n}}}
#'
#' where \eqn{\textrm{TAC}} is the the mean catch recommendation, \eqn{n} is the number of historical years, and
#' \eqn{C_y}  is the catch in historical year \eqn{y}
#'
#' @templateVar mp AvC
#' @template MPtemplate
#' @template MPuses
#' @export
#' @author T. Carruthers
#'
#' @family Average Catch MPs
#'
#' @examples
#' Rec <- AvC(1, MSEtool::Cobia, reps=1000, plot=TRUE) # 1,000 log-normal samples with CV = 0.2
#'
AvC <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@Cat Data@LHYear"
  if (length(Data@Year)<1 | is.na(Data@LHYear[1])) {
    Rec <- new("Rec")
    Rec@TAC <- rep(as.numeric(NA), reps)
    return(Rec)
  }
  yrs <- min(Data@Year):(Data@Year[Data@Year==Data@LHYear[1]])
  yr.ind <- match(yrs, Data@Year)
  histCatch <- Data@Cat[x, yr.ind]
  meanC <- mean(histCatch, na.rm = T)
  if (reps >1) {
    TAC <- rlnorm(reps, log(meanC), 0.2)
  } else {
    TAC <- meanC
  }
  Rec <- new("Rec")
  Rec@TAC <- TAC

  if (plot) AvC_plot(x, Data, Rec, meanC, histCatch, yr.ind, lwd=3, cex.lab=1.25)
  Rec
}
class(AvC) <- "MP"


#### Beddington-Kirkwood Fmax estimation ####

#' Beddington and Kirkwood life-history MP
#'
#' Family of management procedures that sets the TAC by approximation of Fmax
#' based on the length at first capture relative to asymptotic length and the
#' von Bertalanffy growth parameter *K*.
#'
#' @templateVar mp BK
#' @template MPtemplate
#' @template MPuses
#' @export
#' @author T. Carruthers.
#'
#' @details The TAC is calculated as:
#' \deqn{\textrm{TAC} = A F_{\textrm{max}}}
#' where \eqn{A} is (vulnerable) stock abundance, and \eqn{F_{\textrm{max}}} is calculated as:
#' \deqn{F_{\textrm{max}} = \frac{0.6K}{0.67-L_c/L_\infty}}
#' where \eqn{K} is the von Bertalanffy growth coefficient, \eqn{L_c} is the
#' length at first capture, and \eqn{L_\infty} is the von Bertalanffy asymptotic length
#'
#'
#' Abundance (*A*) is either assumed known (`BK`) or estimated (`BK_CC` and `BK_ML`):
#' \deqn{A = \frac{\bar{C}}{\left(1-e^{-F}\right)}}
#' where \eqn{\bar{C}} is the mean catch, and *F* is estimated.
#' See Functions section below for the estimation of *F*.
#'
#' @note
#' Note that the Beddington-Kirkwood method is designed to estimate \eqn{F_\textrm{max}},
#' that is, the fishing mortality that produces the maximum yield *assuming constant
#' recruitment independent of spawning biomass*.
#'
#' Beddington and Kirkwood (2005) recommend estimating *F* using other methods
#' (e.g., a catch curve) and comparing the estimated *F* to the estimated
#' \eqn{F_\textrm{max}} and adjusting exploitation accordingly.
#' These MPs have not been implemented that way.
#'
#' @describeIn BK Assumes that abundance is known, i.e. `Data@Abun`
#' and `Data@CV_abun` contain values
#'

#' @references Beddington, J.R., Kirkwood, G.P., 2005. The estimation of
#' potential yield and stock status using life history parameters. Philos.
#' Trans. R. Soc. Lond. B Biol. Sci. 360, 163-170.
#'
#'
#' @examples
#' \dontrun{
#' BK(1, MSEtool::SimulatedData, reps=1000, plot=TRUE)
#' }
#'
BK <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@LFC, Data@vbLinf, Data@CV_vbLinf, Data@Abun, Data@CV_Abun, Data@vbK, Data@CV_vbK"
  Lc <- MSEtool::trlnorm(reps * 10, Data@LFC[x], 0.2)
  Linfc <- MSEtool::trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Ac <- MSEtool::trlnorm(reps * 10, Data@Abun[x], Data@CV_Abun[x])
  Kc <- MSEtool::trlnorm(reps * 10, Data@vbK[x], Data@CV_vbK[x])
  Fmax <- (0.6 * Kc)/(0.67 - (Lc/Linfc))
  TAC <- Ac * Fmax
  ind <- TAC > 0 & Fmax > 0
  TAC <- MSEtool::TACfilter(TAC[ind][1:reps])  # Interval censor only those positive catch recommendations
  Rec <- new("Rec")
  Rec@TAC <- TAC

  if (plot) {
    DF <- data.frame(vals=c(Lc[ind][1:reps]/Linfc[ind][1:reps], Kc[ind][1:reps],
                            Fmax[ind][1:reps], Ac[ind][1:reps], TAC[ind][1:reps]),
                     vars=rep(c('Lc/Linf', 'K', 'Fmax', 'Abundance', 'TAC'), each=reps),
                     stringsAsFactors = FALSE)
    DF <- DF[!is.na(DF$vals),]
    DF$vars <- factor(DF$vars, levels=c('Lc/Linf', 'K', 'Fmax', 'Abundance', 'TAC'), ordered=TRUE)
    BK_plot(DF)
  }

  Rec

}
class(BK) <- "MP"

#' @param Fmin The minimum fishing mortality rate that is derived from the
#' catch-curve (interval censor).
#' @describeIn BK Abundance is estimated using an age-based catch curve
#' to estimate *Z* and *F*, and abundance estimated from recent catches and *F*.
#' @templateVar mp BK_CC
#' @template MPuses
#' @examples
#' \dontrun{
#' BK_CC(1, MSEtool::SimulatedData, reps=1000, plot=TRUE)
#' }
#'
#' @export
BK_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  dependencies = "Data@LFC, Data@vbLinf, Data@CV_vbLinf, Data@vbK, Data@CV_vbK, Data@CAA, Data@Mort"
  Lc <- MSEtool::trlnorm(reps, Data@LFC[x], 0.2)
  Linfc <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Mdb <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, MuC, Data@CV_Cat[x,1])
  Zdb <- CC(x, Data, reps = reps * 10)
  Fdb <- Zdb - Mdb
  Ac <- Cc/(1 - exp(-Fdb))
  Fmax <- (0.6 * Kc)/(0.67 - (Lc/Linfc))

  ind <- (1:(reps * 10))[Fdb > Fmin][(1:reps)*5]
  Fdb <- Fdb[ind]
  Mdb <- Mdb[ind]
  SM <- sum(is.na(ind))
  if (SM > 0) {
    Mdb[is.na(ind)] <- MSEtool::trlnorm(SM, Data@Mort[x], Data@CV_Mort[x])
    Fdb[is.na(ind)] <- Fmin
  }
  Ac <- Cc/(1 - exp(-Fdb))

  TAC <- Ac * Fmax
  ind <- TAC > 0 & Fmax > 0
  TAC <- MSEtool::TACfilter(TAC[ind][1:reps])

  Rec <- new("Rec")
  Rec@TAC <- TAC

  if (plot) {
    DF <- data.frame(vals=c(Lc[ind][1:reps]/Linfc[ind][1:reps], Kc[ind][1:reps],
                            Fmax[ind][1:reps], Ac[ind][1:reps], TAC[ind][1:reps]),
                     vars=rep(c('Lc/Linf', 'K', 'Fmax', 'Abundance', 'TAC'), each=reps),
                     stringsAsFactors = FALSE)
    DF <- DF[!is.na(DF$vals),]

    DF$vars <- factor(DF$vars, levels=c('Lc/Linf', 'K', 'Fmax', 'Abundance', 'TAC'), ordered=TRUE)
    BK_plot(DF)

  }
  Rec
}
class(BK_CC) <- "MP"

#' @describeIn BK Abundance is estimated using mean length
#' to estimate *Z* and *F*, and abundance estimated from recent catches and *F*.
#' @templateVar mp BK_ML
#' @template MPuses
#' @examples
#' \dontrun{
#' BK_ML(1, MSEtool::SimulatedData, reps=100, plot=TRUE)
#' }
#'
#' @export
BK_ML <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@LFC, Data@vbLinf, Data@CV_vbLinf, Data@vbK, Data@CV_vbK, Data@CAL, Data@Mort"
  Lc <- MSEtool::trlnorm(reps * 10, Data@LFC[x], 0.2)
  Linfc <- MSEtool::trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- MSEtool::trlnorm(reps * 10, Data@vbK[x], Data@CV_vbK[x])
  Mdb <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  Z <- MLne(x, Data, Linfc = Linfc, Kc = Kc, ML_reps = reps * 10, MLtype = "F")
  if (all(is.na(Z))) {
    Rec <- new("Rec")
    Rec@TAC <- rep(as.numeric(NA), reps)
    return(Rec)
  }
  FM <- Z - Mdb

  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps * 10, MuC, Data@CV_Cat[x,1])
  Ac <- Cc/(1 - exp(-FM))
  ind <- Ac>0

  Fmax <- (0.6 * Kc)/(0.67 - (Lc/Linfc))  # robustifying for use in MSE
  TAC <- Ac[ind] * Fmax[ind]
  TAC <- MSEtool::TACfilter(TAC[1:reps])

  Rec <- new("Rec")
  Rec@TAC <- TAC

  if (plot) {
    DF <- data.frame(vals=c(Lc[ind][1:reps]/Linfc[ind][1:reps], Kc[ind][1:reps],
                            Fmax[ind][1:reps], Ac[ind][1:reps], TAC[ind][1:reps]),
                     vars=rep(c('Lc/Linf', 'K', 'Fmax', 'Abundance', 'TAC'), each=reps),
                     stringsAsFactors = FALSE)
    DF <- DF[!is.na(DF$vals),]
    DF$vars <- factor(DF$vars, levels=c('Lc/Linf', 'K', 'Fmax', 'Abundance', 'TAC'), ordered=TRUE)
    BK_plot(DF)
  }

  Rec
}
class(BK_ML) <- "MP"







#' Geromont and Butterworth (2015) Constant Catch
#'
#' The TAC is the average historical catch over the last `yrsmth` (default 5) years,
#' multiplied by (1-`xx`)
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = (1-x)C_{\textrm{ave}}}
#' where *x* lies between 0 and 1, and \eqn{C_{\textrm{ave}}} is average historical
#' catch over the previous `yrsmth` years.
#'
#' The TAC is constant for all future projections.
#'
#' @param yrsmth Years over which to calculate mean catches
#' @param xx Parameter controlling the TAC. Mean catches are multiplied by
#' (1-`xx`)
#'
#' @template MPtemplate
#' @templateVar mp CC1
#' @template MPuses
#' @export
#' @family Constant Catch MPs
#' @author T. Carruthers
#' @references Geromont, H. F., and D. S. Butterworth. 2015.
#' Generic Management Procedures for Data-Poor Fisheries: Forecasting with Few Data.
#'  ICES Journal of Marine Science: Journal Du Conseil 72 (1). 251-61.
#'
#' @describeIn CC1 TAC is average historical catch from recent `yrsmth` years
#' @examples
#' CC1(1, MSEtool::Cobia, plot=TRUE)
#'

#### Constant Catch MPs ####
CC1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0) {
  dependencies = "Data@Cat, Data@CV_Cat"
  if (length(Data@Year)<1) {
    Rec <- new("Rec")
    Rec@TAC <- rep(as.numeric(NA), reps)
    return(Rec)
  }

  yrlast <- match(Data@LHYear[1], Data@Year)
  yrfirst <- yrlast - yrsmth + 1
  # C_dat <- Data@Cat[x, (length(Data@Year) - (yrsmth - 1)):length(Data@Year)]
  C_dat <- Data@Cat[x, yrfirst:yrlast]
  TAC <- (1 - xx) * MSEtool::trlnorm(reps, mean(C_dat, na.rm=TRUE), Data@CV_Cat[x,1]/(yrsmth^0.5))  # mean catches over the interval
  Rec <- new("Rec")
  Rec@TAC <- MSEtool::TACfilter(TAC)
  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    par(mfrow=c(1,1))
    ylim <- c(0, max(c(Data@Cat[x,], MSEtool::TACfilter(TAC))))
    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x,],NA), type="l", lwd=2, las=1, bty="l",
         xlab="Year", ylab=paste0("Catch (", Data@Units, ")"),
         cex.lab=1.5, cex.axis=1.25, ylim=ylim)
    abline(v=Data@LHYear[1], lty=3, col="darkgray")
    lines(Data@Year[yrfirst:yrlast], rep(mean(C_dat, na.rm=TRUE), yrsmth), col="blue", lwd=3)
    boxplot(Rec@TAC, at=max(Data@Year)+1, add=TRUE, axes=FALSE)
  }

  Rec
}
class(CC1) <- "MP"

#' @describeIn CC1 TAC is average historical catch from recent `yrsmth` years reduced by 10\%.
#' @examples
#' CC2(1, MSEtool::Cobia, plot=TRUE)
#'
#' @export
CC2 <- CC1
formals(CC2)$xx <- 0.1
class(CC2) <- "MP"

#' @describeIn CC1 TAC is average historical catch from recent `yrsmth` years reduced by 20\%.
#' @examples
#' CC3(1, MSEtool::Cobia, plot=TRUE)
#'
#' @export
CC3 <- CC1
formals(CC3)$xx <- 0.2
class(CC3) <- "MP"


#' @describeIn CC1 TAC is average historical catch from recent `yrsmth` years reduced by 30\%.
#' @examples
#' CC4(1, MSEtool::Cobia, plot=TRUE)
#'
#' @export
CC4 <- CC1
formals(CC4)$xx <- 0.3
class(CC4) <- "MP"

#' @describeIn CC1 TAC is average historical catch from recent `yrsmth` years reduced by 40\%.
#' @examples
#' CC5(1, MSEtool::Cobia, plot=TRUE)
#'
#' @export
CC5 <- CC1
formals(CC5)$xx <- 0.4
class(CC5) <- "MP"

#' @describeIn CC1 TAC is fixed at last historical catch
#' @examples
#' CurC(1, MSEtool::Cobia, plot=TRUE)
#'
#' @export
CurC <- CC1
formals(CurC)$yrsmth <- 1
class(CurC) <- "MP"


#' Geromont and Butterworth Constant Catch Harvest Control Rule
#'
#' A simple MP that aims for a reference catch (as a proxy for MSY)
#' subject to imperfect information.
#'
#' Note that this is my interpretation of their MP and is now stochastic.
#' Currently it is generalized and is not 'tuned' to more detailed assessment
#' data which might explain why in some cases it leads to stock declines.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = C_\textrm{ref}}
#' where \eqn{C_\textrm{ref}} is a reference catch assumed to be a proxy for MSY.
#' In the MSE \eqn{C_\textrm{ref}} is the calculated MSY subject to observation error
#' defined in `Obs@CV_Cref`.
#'
#' The TAC is subject to the following conditions:
#' 1. if next TAC > 1.2 last catch, then TAC = 1.2 last catch
#' 2. if next TAC < 0.8 last catch, then TAC = 0.8 last catch
#'
#' @template MPtemplate
#' @templateVar mp GB_CC
#' @template MPuses
#'
#' @author T. Carruthers
#' @references Geromont, H.F. and Butterworth, D.S. 2014. Complex assessment or
#' simple management procedures for efficient fisheries management: a
#' comparative study. ICES J. Mar. Sci. doi:10.1093/icesjms/fsu017
#' @export
#' @examples
#' GB_CC(1, MSEtool::SimulatedData, plot=TRUE)
#' @family Constant Catch MPs
GB_CC <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@Cref,Data@Cat"
  Catrec <- Data@Cat[x, length(Data@Cat[x, ])]
  TAC <- MSEtool::trlnorm(reps, Data@Cref[x], Data@CV_Cref)
  TAC[TAC > (1.2 * Catrec)] <- 1.2 * Catrec
  TAC[TAC < (0.8 * Catrec)] <- 0.8 * Catrec
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) GB_CC_plot(x, Catrec, TAC, Data)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(GB_CC) <- "MP"



#### Age-Comp SRA ####
#' Internal function for CompSRA MP
#'
#' @param x Simulation number
#' @param Data An object of class Data
#' @param reps Number of reps
#'
#' @return A list
#' @keywords internal
#' @export
#'
CompSRA_ <- function(x, Data, reps=100) {
 maxage <- Data@MaxAge
 TAC <- Bt_K <- FMSY <- Ac <- rep(NA, reps)
 predout <- list()
 for (i in 1:reps) {
   Mc <- MSEtool::trlnorm(1, Data@Mort[x], Data@CV_Mort)
   hc <- MSEtool::sample_steepness2(1, Data@steep[x], Data@CV_steep[x])
   Linfc <- MSEtool::trlnorm(1, Data@vbLinf[x], Data@CV_vbLinf[x])
   Kc <- MSEtool::trlnorm(1, Data@vbK[x], Data@CV_vbK[x])
   if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
     t0c <- -MSEtool::trlnorm(1, -Data@vbt0[x], Data@CV_vbt0[x])
   } else {
     t0c <- Data@vbt0[x]
   }
   t0c[!is.finite(t0c)] <- 0
   LFSc <- MSEtool::trlnorm(1, Data@LFS[x], Data@CV_LFS[x])
   LFCc <- MSEtool::trlnorm(1, Data@LFC[x], Data@CV_LFC[x])
   AMc <- MSEtool::trlnorm(1, iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x], Data@L50[x]), Data@CV_L50[x])
   ac <- MSEtool::trlnorm(1, Data@wla[x], Data@CV_wla[x])
   bc <- MSEtool::trlnorm(1, Data@wlb[x], Data@CV_wlb[x])
   Catch <- Data@Cat[x, ]
   ny <- length(Catch)
   nyCAA <- dim(Data@CAA)[2]
   CAA <- Data@CAA[x, max(nyCAA - 2, 1):nyCAA, ]  # takes last three years as the sample (or last year if there is only one)
   Nac <- exp(-Mc * ((1:maxage) - 1))  # put a rough range on estimate of R0 assuming a mean harvest rate of 10%
   Lac <- Linfc * (1 - exp(-Kc * ((1:maxage) - t0c)))
   Wac <- ac * Lac^bc
   AFC <- log(1 - min(0.99, LFCc/Linfc))/-Kc + t0c
   AFS <- log(1 - min(0.99, LFSc/Linfc))/-Kc + t0c
   if (AFC >= 0.7 * maxage) AFC <- 0.7 * maxage
   if (AFS >= 0.9 * maxage) AFS <- 0.9 * maxage
   KES <- max(2, ceiling(mean(c(AFC, AFS))))
   pred <- Nac * Wac
   pred[1:(KES - 1)] <- 0
   pred <- pred/sum(pred)
   pred <- ((mean(Catch, na.rm=TRUE)/0.1) * pred/Wac)/exp(-(1:maxage) * Mc)
   pred <- pred[pred > 0]
   R0range <- c(mean(pred)/1000, mean(pred) * 1000)
   fit <- optimize(SRAfunc, log(R0range), Mc, hc, maxage, LFSc, LFCc, Linfc,
                   Kc, t0c, AMc, ac, bc, Catch, CAA)
   getvals <- SRAfunc(fit$minimum, Mc, hc, maxage, LFSc, LFCc, Linfc, Kc, t0c,
                      AMc, ac, bc, Catch, CAA, opt = 2)
   Ac[i] <- getvals$B
   Bt_K[i] <-  getvals$D
   predout[[i]] <- getvals$pred
   fit2 <- optimize(SRAFMSY, log(c(1e-04, 3)), Mc, hc, maxage, LFSc, LFCc, Linfc, Kc, t0c, AMc, ac, bc)
   FMSY[i] <- exp(fit2$minimum)
   if ((FMSY[i]/Mc) > 3) FMSY[i] <- 3 * Mc
   TAC[i] <- Ac[i] * FMSY[i]
 }
 return(list(TAC=TAC, Bt_K=Bt_K, FMSY=FMSY, Ac=Ac, pred=predout, CAA=CAA))
}

#' Age-Composition Stock-Reduction Analysis
#'
#' A stock reduction analysis (SRA) model is fitted to the age-composition
#' from the last 3 years (or less if fewer data are available)
#'
#' A stock reduction analysis (SRA) model is fitted to the age-composition
#' from the last 3 years (or less if fewer data are available) assuming a constant
#' total mortality rate (*Z*) and used to estimate current stock depletion (*D*),
#'  \eqn{F_\textrm{MSY}}, and stock abundance (*A*).
#'
#'
#' Abundance is estimated in the SRA. \eqn{F_{\textrm{MSY}}} is calculated assuming
#' knife-edge vulnerability at the age of full selection.
#'
#' The TAC is calculated as \eqn{F_\textrm{MSY} A}. `CompSRA4010` uses a 40-10
#' harvest control rule to reduce TAC at low biomass.
#'
#' @templateVar mp CompSRA
#' @template MPtemplate
#' @template MPuses
#'
#' @describeIn CompSRA TAC is FMSY x Abundance
#' @author T. Carruthers
#'
#' @examples
#' \dontrun{
#' CompSRA(1, MSEtool::SimulatedData, plot=TRUE)
#' }
#'
#' @export
CompSRA <- function(x, Data, reps = 100, plot=FALSE) {
 # optimize for fixed F to get you to current depletion C/Fcur =
 # abundance
 dependencies = "Data@Mort, Data@CV_Mort, Data@vbK, Data@CV_vbK, Data@vbLinf, Data@CV_vbLinf, Data@vbt0, Data@CV_vbt0, Data@MaxAge, Data@wla, Data@CV_wla, Data@wlb, Data@CV_wlb, Data@L50, Data@CV_L50, Data@CAA, Data@steep, Data@CV_steep, Data@LFS, Data@CV_LFS, Data@LFC, Data@CV_LFC, Data@Cat"
 runCompSRA <- CompSRA_(x, Data, reps)
 TAC <- MSEtool::TACfilter(runCompSRA$TAC)
 if (plot) CompSRA_plot(runCompSRA, TAC)
 Rec <- new("Rec")
 Rec@TAC <- TAC
 Rec
}
class(CompSRA) <- "MP"
#' @templateVar mp CompSRA4010
#' @template MPuses
#' @describeIn CompSRA With a 40-10 control rule based on estimated depletion
#' @examples
#' CompSRA4010(1, MSEtool::SimulatedData, plot=TRUE)
#'
#' @export CompSRA4010
CompSRA4010 <- function(x, Data, reps = 100, plot=FALSE) {
 dependencies = "Data@Mort, Data@CV_Mort, Data@vbK, Data@CV_vbK, Data@vbLinf, Data@CV_vbLinf, Data@vbt0, Data@CV_vbt0, Data@MaxAge, Data@wla, Data@CV_wla, Data@wlb, Data@CV_wlb, Data@L50, Data@CV_L50, Data@CAA, Data@steep, Data@CV_steep, Data@LFS, Data@CV_LFS, Data@LFC, Data@CV_LFC, Data@Cat"
 runCompSRA <- CompSRA_(x, Data, reps)
 TAC <- MSEtool::TACfilter(runCompSRA$TAC)
 Bt_K <- runCompSRA$Bt_K
 # 40-10 rule
 cond1 <- Bt_K < 0.4 & Bt_K > 0.1
 cond2 <- Bt_K < 0.1
 TAC[cond1] <- TAC[cond1] * (Bt_K[cond1] - 0.1)/0.3
 TAC[cond2] <- TAC[cond2] * tiny  # this has to still be stochastic albeit very small
 if (plot) CompSRA_plot(runCompSRA, TAC)
 Rec <- new("Rec")
 Rec@TAC <- TAC
 Rec
}
class(CompSRA4010) <- "MP"
# ---- Depletion Corrected Average Catch MPs -----

#' Depletion Corrected Average Catch
#'
#' Internal code to calculate dcac
#'
#' @param x Simulation number
#' @param Data An object of class `Data`
#' @param reps The number of reps
#' @param Bt_K Optional value of depletion (otherwise taken from Data object)
#' @param updateD Logical. Should depletion be updated in projection years?
#'
#' @return A list with dcac, Bt_K, and BMSY_K
#' @keywords internal
#' @export
#'
DCAC_ <- function(x, Data, reps=100, Bt_K=NULL, updateD=FALSE) {
  if (length(Data@Year)<1) return(list(dcac=NA, Bt_K=NA, BMSY_K=NA))
  if (MSEtool::NAor0(Data@BMSY_B0[x])) stop('Data@BMSY_B0 is NA')
  if (MSEtool::NAor0(Data@CV_BMSY_B0[x])) stop("Data@CV_BMSY_B0 is NA")
  yr.lst <- match(Data@LHYear[1], Data@Year)
  yrs <- 1:yr.lst
  C_tot <-  Data@AvC[x] * Data@t[x]
  Mdb <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])  # CV of 0.5 as in MacCall 2009
  FMSY_M <- MSEtool::trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x])  # standard deviation of 0.2 - referred to as 'standard error' in MacCall 2009
  if (is.null(Bt_K))  Bt_K <- MSEtool::trlnorm(reps, Data@Dt[x], Data@CV_Dt[x])
  if (!updateD) {
    if (Data@LHYear[1] != max(Data@Year)) {
      dcac <- rep(Data@MPrec[x], reps) # catch limit is static for future
      return(list(dcac=dcac, Bt_K=Bt_K))
    }
  }

  Bt_K[Bt_K>1] <-1
  Bt_K[Bt_K<0] <-0
  if (any(is.na(c(Data@BMSY_B0[x], Data@CV_BMSY_B0[x])))) {
    warning("Data@BMSY_B0 or Data@CV_BMSY_B0 do not contain values")
    return(list(dcac=rep(NA, reps), Bt_K=Bt_K))
  }
  BMSY_K <- rbeta(reps, MSEtool::alphaconv(Data@BMSY_B0[x], Data@BMSY_B0[x] * Data@CV_BMSY_B0[x]), MSEtool::betaconv(Data@BMSY_B0[x], Data@BMSY_B0[x] * Data@CV_BMSY_B0[x]))

  dcac <- C_tot/(yr.lst + ((1 - Bt_K)/(BMSY_K * FMSY_M * Mdb)))

  return(list(dcac=dcac, Bt_K=Bt_K, BMSY_K=BMSY_K))

}


#' Depletion Corrected Average Catch
#'
#' This group of MPs calculates a catch limit (*dcac*; intended as an MSY proxy) based on
#' average historical catch while accounting for the windfall catch that got the stock down to
#' its current depletion level (*D*).
#'
#' The method calculates the depletion-corrected average catch (*dcac*) as:
#'
#' \deqn{\textrm{dcac} = \frac{\sum_{y=1}^{n}{C_y}}{n+(1-D)/Y_{\textrm{pot}}}}
#'
#' where
#' \deqn{Y_{\textrm{pot}} = \frac{B_{\textrm{MSY}}}{B_0}\frac{F_{\textrm{MSY}}}{M}M}
#' and \eqn{C} is the historical catches; i.e \eqn{C} does not change in the future projections in the MSE
#'
#' The methods differ in the assumptions of current depletion (*D*). See the **Functions** section below for details.
#'
#' @templateVar mp DCAC
#' @template MPtemplate
#'
#'
#' @note It's probably worth noting that DCAC TAC recommendations do not tend
#' to zero as depletion tends to zero. It adjusts for depletion only in
#' calculating historical average catch. It follows that at stock levels much
#' below BMSY, DCAC tends to chronically overfish.
#'
#' @references
#' MacCall, A.D., 2009. Depletion-corrected average catch: a simple
#' formula for estimating sustainable yields in data-poor situations. ICES J.
#' Mar. Sci. 66, 2267-2271.
#'
#' @author T. Carruthers
#' @family Average Catch MPs
#' @describeIn DCAC Depletion is not updated in the future projections.
#' The TAC is static and not updated in the future years. This represents
#' an application of the DCAC method where a catch limit is calculated based on
#' current estimate of depletion and time-series of catch from the beginning
#' of the fishery, and the TAC is fixed at this level for all future projections.
#'
#' @templateVar mp DCACs
#' @template MPuses
#' @examples
#' DCACs(1, MSEtool::Atlantic_mackerel, plot=TRUE)
#'
#' @export
DCACs <- function(x, Data, reps = 100, plot=FALSE) {
  rundcac <- DCAC_(x, Data, reps)
  TAC <- MSEtool::TACfilter(rundcac$dcac)
  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DCACs) <- "MP"


#' @describeIn DCAC Depletion is estimated each management interval and used to
#' update the catch limit recommendation based on the historical catch (which is
#' not updated in the future projections).
#' @templateVar mp DCAC
#' @template MPuses
#' @examples
#' DCAC(1, MSEtool::Atlantic_mackerel, plot=TRUE)
#'
#' @export
DCAC <- function(x, Data, reps = 100, plot=FALSE) {
  rundcac <- DCAC_(x, Data, reps, updateD=TRUE)
  TAC <- MSEtool::TACfilter(rundcac$dcac)

  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DCAC) <- "MP"



#' @describeIn DCAC Current stock biomass is assumed to be exactly at 40 per cent
#' of unfished levels. The 40 percent depletion assumption may not
#' really affect DCAC that much as it already makes TAC recommendations that are
#' quite MSY-like.
#' @templateVar mp DCAC_40
#' @template MPuses
#' @examples
#' DCAC_40(1, MSEtool::Atlantic_mackerel, plot=TRUE)
#'
#' @export
DCAC_40 <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@AvC, Data@t, Data@Mort, Data@CV_Mort, Data@FMSY_M, Data@CV_FMSY_M, Data@BMSY_B0, Data@CV_BMSY_B0"

  rundcac <- DCAC_(x, Data, reps, Bt_K=0.4)
  TAC <- MSEtool::TACfilter(rundcac$dcac)
  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DCAC_40) <- "MP"


#' @describeIn DCAC The dynamic DCAC (depletion is updated) is paired with the 40-10 rule that throttles
#' back the OFL to zero at 10 percent of unfished stock size (the OFL is not subject to downward
#' adjustment above 40 percent unfished). DCAC can overfish below BMSY levels. The 40-10
#' harvest control rule largely resolves this problem providing an MP with surprisingly good
#' performance even at low stock levels.
#' @templateVar mp DCAC4010
#' @template MPuses
#' @examples
#' Data <- MSEtool::Atlantic_mackerel
#' Data@LHYear <- 2005
#' DCAC4010(1, Data, plot=TRUE)
#'
#' @export
DCAC4010 <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@AvC, Data@t, Data@Mort, Data@CV_Mort, Data@FMSY_M, Data@CV_FMSY_M, Data@Dt, Data@CV_Dt, Data@BMSY_B0, Data@CV_BMSY_B0"
  rundcac <- DCAC_(x, Data, reps, updateD=TRUE)
  dcac <- rundcac$dcac
  Bt_K <- rundcac$Bt_K
  TAC <- dcac

  # 40-10 rule
  cond1 <- Bt_K < 0.4 & Bt_K > 0.1
  cond2 <- Bt_K < 0.1
  if (length(cond1) > 0) TAC[cond1] <- TAC[cond1] * (Bt_K[cond1] - 0.1)/0.3
  if (length(cond2) > 0)  TAC[cond2] <- TAC[cond2] * tiny  # this has to still be stochastic albeit very small
  if (length(cond1) < 1 & length(cond2) < 1)  return(NA)
  TAC <- MSEtool::TACfilter(TAC)

  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DCAC4010) <- "MP"



#' @describeIn DCAC This variant uses the mean length estimator to calculate current stock
#' depletion. The mean length extension was programmed by Gary Nelson as part of his
#' excellent R package 'fishmethods'.
#' @templateVar mp DCAC_ML
#' @template MPuses
#' @examples
#' DCAC_ML(1, MSEtool::SimulatedData, plot=TRUE)
#'
#' @export
DCAC_ML <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@AvC, Data@t, Data@Mort, Data@CV_Mort, Data@FMSY_M, Data@CV_FMSY_M, Data@BMSY_B0, Data@CV_BMSY_B0, Data@Year, Data@CAL, Data@vbLinf, Data@CV_vbLinf, Data@vbK, Data@CV_vbK"
  if (is.na(Data@BMSY_B0[x]) | is.na(Data@CV_BMSY_B0[x])) return(NA)
  if (is.na(Data@FMSY_M[x]) | is.na(Data@CV_FMSY_M[x])) return(NA)

  Mdb <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  FMSY_M <- MSEtool::trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x])
  Linfc <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Z <- MLne(x, Data, Linfc = Linfc, Kc = Kc, ML_reps = reps, MLtype = "dep")
  if (all(is.na(Z))) {
    Rec <- new("Rec")
    Rec@TAC <- rep(as.numeric(NA), reps)
    return(Rec)
  }
  FM <- Z - Mdb
  FM[FM<0.05] <- 0.05

  nyears <- length(Data@Year)
  Ct1 <- mean(Data@Cat[x, 1:3], na.rm=TRUE)
  Ct2 <- mean(Data@Cat[x, (nyears - 2):nyears], na.rm=TRUE)
  dep <- rep(c(Ct1, Ct2), each = reps)/(1 - exp(-FM))
  if (reps == 1) {
    Bt_K <- dep[2]/dep[1]
    Bt_K <- min(Bt_K, 0.8) # make depletion < 0.8 # sometimes estimated dep >> 1
  }
  if (reps > 1) {
    Bt_K <- dep[, 2]/dep[, 1]
    for (rr in 1:reps) Bt_K[rr] <- min(Bt_K[rr], 0.8)
  }


  rundcac <- DCAC_(x, Data, reps, Bt_K=Bt_K)
  TAC <- MSEtool::TACfilter(rundcac$dcac)

  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DCAC_ML) <- "MP"


#' @describeIn DCAC Depletion Adjusted Average Catch: essentially DCAC (with updated Depletion)
#'  divided by BMSY/B0 (Bpeak) (Harford and Carruthers, 2017).
#' @references
#' Harford W. and Carruthers, T. 2017. Interim and long-term performance of
#' static and adaptive management procedures. Fish. Res. 190, 84-94.
#' @templateVar mp DAAC
#' @template MPuses
#'
#' @examples
#' Data <- MSEtool::Atlantic_mackerel
#' Data@LHYear <- 2005
#' DAAC(1, Data, plot=TRUE)
#'
#' @export
DAAC <- function(x, Data, reps = 100, plot=FALSE) {
  # extended depletion-corrected average catch (Harford and Carruthers 2017)
  dependencies = "Data@AvC, Data@t, Data@Mort, Data@CV_Mort, Data@FMSY_M, Data@CV_FMSY_M, Data@Dt, Data@CV_Dt, Data@BMSY_B0, Data@CV_BMSY_B0"

  rundcac <- DCAC_(x, Data, reps, updateD = TRUE)
  TAC <- rundcac$dcac * rundcac$Bt_K/rundcac$BMSY_K

  TAC <- MSEtool::TACfilter(TAC)
  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DAAC) <- "MP"



#' @describeIn DCAC Hybrid Depletion Adjusted Average Catch: essentially DCAC  (with updated Depletion)
#' divided by BMSY/B0 (Bpeak) when below BMSY,
#' and DCAC above BMSY (Harford and Carruthers 2017).
#' @templateVar mp HDAAC
#' @template MPuses
#' @examples
#' Data <- MSEtool::Atlantic_mackerel
#' Data@LHYear <- 2005
#' HDAAC(1, Data, plot=TRUE)
#'
#' @export
HDAAC <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@AvC, Data@t, Data@Mort, Data@CV_Mort, Data@FMSY_M, Data@CV_FMSY_M, Data@Dt, Data@CV_Dt, Data@BMSY_B0, Data@CV_BMSY_B0"

  rundcac <- DCAC_(x, Data, reps, updateD = TRUE)
  TAC <- rundcac$dcac
  ddcac <- rundcac$dcac * rundcac$Bt_K/rundcac$BMSY_K
  TAC[rundcac$Bt_K < rundcac$BMSY_K] <- ddcac[rundcac$Bt_K < rundcac$BMSY_K]

  TAC <- MSEtool::TACfilter(TAC)

  if (plot)  DCAC_plot(x, Data, dcac=rundcac$dcac, TAC, Bt_K=rundcac$Bt_K, yrs=1:length(Data@Year))

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec

}
class(HDAAC) <- "MP"


#### Depletion Based SRA ####

#' Internal optimization function
#'
#' @param lnK log unfished biomass
#' @param C_hist Historical catches
#' @param nys number of years (length of C_hist)
#' @param Mdb M
#' @param FMSY_M  FMSY/M
#' @param BMSY_K BMSY/K
#' @param Bt_K Current depletion
#' @param adelay age of maturity
#' @param opt Numeric. 1 for optimization. 2 to return Biomass trend
#'
#' @keywords internal
#' @export
#'
DBSRAopt <- function(lnK, C_hist, nys, Mdb, FMSY_M, BMSY_K, Bt_K, adelay, opt=1) {
  # the optimization for B0 given DBSRA assumptions
  Kc <- exp(lnK)
  fn <- function(n, BMSY_K) {
    # optimizer to find parameter n according to sampled BMSY/B0 (theta)
    thetapred <- n^(-1/(n - 1))
    (BMSY_K - thetapred)^2
  }
  n <- optimize(fn, c(0.01, 6), BMSY_K = BMSY_K)$minimum  #get the optimum getn(BMSY_K)
  g <-  (n^(n/(n - 1)))/(n - 1)# gety(n) get the y parameter for n
  FMSY <- FMSY_M * Mdb
  UMSY <- (FMSY/(FMSY + Mdb)) * (1 - exp(-(FMSY + Mdb)))
  MSY <- Kc * BMSY_K * UMSY
  # Bjoin rules from Dick & MacCall 2011
  Bjoin_K <- 0.5
  if (BMSY_K < 0.3) Bjoin_K <- 0.5 * BMSY_K
  if (BMSY_K > 0.3 & BMSY_K < 0.5) Bjoin_K <- 0.75 * BMSY_K - 0.075
  Bjoin <- Bjoin_K * Kc
  prodPTF <- function(depletion, n, MSY) {
    # Pella-Tomlinson production function required for DB-SRA
    y <- (n^(n/(n - 1)))/(n - 1)
    MSY * y * depletion - MSY * y * depletion^n
  }
  PBjoin <- prodPTF(Bjoin_K, n, MSY)
  cp <- (1 - n) * g * MSY * (Bjoin^(n - 2)) * Kc^-n
  Bc <- rep(NA, nys)
  Bc[1] <- Kc
  obj <- 0
  for (yr in 2:nys) {
    yref <- max(1, yr - adelay)
    if (Bc[yref] > Bjoin | BMSY_K > 0.5) {
      Bc[yr] <- Bc[yr - 1] + g * MSY * (Bc[yref]/Kc) - g * MSY *
        (Bc[yref]/Kc)^n - C_hist[yr - 1]
    } else {
      Bc[yr] <- Bc[yr - 1] + Bc[yref] * ((PBjoin/Bjoin) + cp * (Bc[yref] - Bjoin)) - C_hist[yr - 1]
    }
    if (Bc[yr] < 0)
      obj <- obj + log(-Bc[yr])
    Bc[yr] <- max(1e-06, Bc[yr])
  }
  if (opt ==1) return(obj + ((Bc[nys]/Kc) - Bt_K)^2)
  if (opt ==2) return(Bc)

}  # end of DBSRA optimization function

#' Depletion-based SRA internal function
#'
#' @param x iteration number
#' @param Data object of class Data
#' @param reps number of reps
#' @param depo Optional fixed depletion (single value)
#' @param hcr Optional harvest control rule for throttling catch as a function of B/B0.
#' Numeric vector of length 2 specifying HCR break points - e.g c(0.4, 0.1) for 40-10 HCR
#'
#' @export
#'
#' @keywords internal
#'
DBSRA_ <- function(x, Data, reps = 100, depo=NULL, hcr=NULL) {
  # returns a vector of DBSRA estimates of the TAC for a particular
  # simulation x for(x in 1:nsim){
  if (MSEtool::NAor0(Data@CV_BMSY_B0[x])) stop("Data@CV_BMSY_B0 is NA")
  C_hist <- Data@Cat[x, ]
  TAC <- rep(NA, reps)
  Btrend <- matrix(NA, nrow=reps, ncol=length(C_hist))
  Bt_Kstore <- FMSY_Mstore <- BMSY_K_Mstore <- rep(NA, reps)
  DBSRAcount <- 1
  if (is.null(depo)) {
    if (is.na(Data@Dep[x]) | is.na(Data@CV_Dep[x])) {
       out <- new("Rec")
       out@TAC <- rep(as.numeric(NA), reps)
       return(out)
     }
  } else {
    Data@CV_Dep[x] <- tiny
  }
  if (is.na(Data@BMSY_B0[x])){
    out <- new("Rec")
    out@TAC <- rep(as.numeric(NA), reps)
    return(out)
  }

  while (DBSRAcount < (reps + 1)) {
    if (is.null(depo)) depo <- max(0.01, min(0.99, Data@Dep[x]))  # known depletion is between 1% and 99% - needed to generalise the Dick and MacCall method to extreme depletion scenarios
    Bt_K <- rbeta(100, MSEtool::alphaconv(depo, min(depo * Data@CV_Dep[x], (1 - depo) * Data@CV_Dep[x])),
                  MSEtool::betaconv(depo, min(depo * Data@CV_Dep[x], (1 - depo) * Data@CV_Dep[x])))  # CV 0.25 is the default for Dick and MacCall mu=0.4, sd =0.1
    Bt_K <- Bt_K[Bt_K > 0.00999 & Bt_K < 0.99001][1]  # interval censor (0.01,0.99)  as in Dick and MacCall 2011

    Mdb <- MSEtool::trlnorm(100, Data@Mort[x], Data@CV_Mort[x])
    Mdb <- Mdb[Mdb < 0.9][1]  # !!!! maximum M is 0.9   interval censor
    if (is.na(Mdb)) Mdb <- 0.9  # !!!! maximum M is 0.9   absolute limit
    FMSY_M <- MSEtool::trlnorm(1, Data@FMSY_M[x], Data@CV_FMSY_M[x])
    BMSY_K <- rbeta(100, MSEtool::alphaconv(Data@BMSY_B0[x], Data@CV_BMSY_B0[x] *  Data@BMSY_B0[x]),
                    MSEtool::betaconv(Data@BMSY_B0[x], Data@CV_BMSY_B0[x] * Data@BMSY_B0[x]))  #0.045 corresponds with mu=0.4 and quantile(BMSY_K,c(0.025,0.975)) =c(0.31,0.49) as in Dick and MacCall 2011
    tryBMSY_K <- BMSY_K[BMSY_K > 0.05 & BMSY_K < 0.95][1]  # interval censor (0.05,0.95) as in Dick and MacCall, 2011
    if (is.na(tryBMSY_K)) {
      Min <- min(BMSY_K, na.rm = TRUE)
      Max <- max(BMSY_K, na.rm = TRUE)
      if (Max <= 0.05) BMSY_K <- 0.05
      if (Min >= 0.95) BMSY_K <- 0.95
    }
    if (!is.na(tryBMSY_K))  BMSY_K <- tryBMSY_K

    Bt_Kstore[DBSRAcount] <- Bt_K
    FMSY_Mstore[DBSRAcount] <- FMSY_M
    BMSY_K_Mstore[DBSRAcount] <- BMSY_K

    adelay <- max(floor(iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x],  Data@L50[x])), 1)

    # scale catches for optimization
    scaler <- 1000/mean(C_hist, na.rm=TRUE)
    C_hist2 <- scaler * C_hist
    opt <- optimize(DBSRAopt, log(c(0.01 * mean(C_hist2, na.rm=TRUE), 1000 * mean(C_hist2, na.rm=TRUE))), C_hist = C_hist2,
                    nys = length(C_hist2), Mdb = Mdb, FMSY_M = FMSY_M, BMSY_K = BMSY_K,
                    Bt_K = Bt_K, adelay = adelay, tol = 0.01)

    Bctemp <- DBSRAopt(opt$minimum,  C_hist = C_hist2,
                   nys = length(C_hist2), Mdb = Mdb, FMSY_M = FMSY_M, BMSY_K = BMSY_K,
                   Bt_K = Bt_K, adelay = adelay, opt=2)

    Btrend[DBSRAcount,] <- Bctemp/scaler
    Kc <- exp(opt$minimum) / scaler
    BMSYc <- Kc * BMSY_K
    FMSYc <- Mdb * FMSY_M
    UMSYc <- (FMSYc/(FMSYc + Mdb)) * (1 - exp(-(FMSYc + Mdb)))
    MSYc <- Kc * BMSY_K * UMSYc
    TAC[DBSRAcount] <- UMSYc * Kc * Bt_K

    if(!is.null(hcr)) {
      if (length(hcr)!=2) stop("hcr must be numeric vector of length 2")
      # 40-10 rule
      if (Bt_K < hcr[1] & Bt_K > hcr[2])  TAC[DBSRAcount] <- TAC[DBSRAcount] * (Bt_K - hcr[2])/(hcr[1]-hcr[2])
      if (Bt_K < hcr[2]) TAC[DBSRAcount] <- TAC[DBSRAcount] * tiny  # this has to still be a numeric value, albeit very small
    }
    DBSRAcount <- DBSRAcount + 1

  }  # end of reps
  list(TAC=TAC, Btrend=Btrend, C_hist=C_hist, Bt_Kstore=Bt_Kstore, FMSY_Mstore=FMSY_Mstore,
       BMSY_K_Mstore=BMSY_K_Mstore, hcr=hcr)

}


#' Depletion-Based Stock Reduction Analysis
#'
#' Depletion-Based Stock Reduction Analysis (DB-SRA) is a method designed for
#' determining a catch limit and management reference points for data-limited
#' fisheries where catches are known from the beginning of exploitation.
#' User prescribed BMSY/B0, M, FMSY/M are used to find B0 and therefore the a catch limit
#' by back-constructing the stock to match a user specified level of stock
#' depletion.
#'
#'
#' DB-SRA assumes that a complete time-series of catch from the beginning of
#' exploitation is available. Users prescribe estimates of current depletion \eqn{(D)},
#' biomass at MSY relative to unfished \eqn{\left(\frac{B_\textrm{MSY}}{B_0}\right)},
#' the natural mortality rate \eqn{(M)}, and the ratio fishing mortality at
#' MSY to M \eqn{\left(\frac{F_{\textrm{MSY}}}{M}\right)}.
#'
#'
#' You may have noticed that you -the user- specify three of the factors that
#' make the quota recommendation. So this can be quite a subjective method.
#' In the MSE the MSY reference points (e.g., \eqn{\left(\frac{F_\textrm{MSY}}{M}\right)})
#' are taken as the true value calculate in the MSE with added uncertainty specified
#' in the Obs object (e.g `Obs@FMSY_Mbiascv`).
#'
#' The catch limit, for the Base Version, is calculated as:
#' \deqn{\textrm{TAC} = M \frac{F_{\textrm{MSY}}}{M} D B_0}
#'
#'
#' @templateVar mp DBSRA
#' @template MPtemplate
#' @template MPuses
#'
#' @note The DB-SRA method of this package isn't exactly the same as the
#' original method of Dick and MacCall (2011) because it has to work for
#' simulated depletions above BMSY/B0 and even on occasion over B0. It also
#' doesn't have the modification for flatfish life histories that has
#' previously been applied by Dick and MacCall (2011).
#'
#' @author T. Carruthers
#' @references
#' Dick, E.J., MacCall, A.D., 2010. Estimates of sustainable yield for 50 data-poor
#' stocks in the Pacific Coast groundfish fishery management plan. Technical memorandum.
#' Southwest fisheries Science Centre, Santa Cruz, CA. National Marine Fisheries Service,
#' National Oceanic and Atmospheric Administration of the U.S. Department of Commerce.
#' NOAA-TM-NMFS-SWFSC-460.
#'
#' Dick, E.J., MacCall, A.D., 2011. Depletion-Based Stock Reduction
#' Analysis: A catch-based method for determining sustainable yields for
#' data-poor fish stocks. Fish. Res. 110, 331-341.
#'
#' @describeIn DBSRA Base Version. TAC is calculated assumed MSY harvest rate
#' multiplied by the estimated current abundance (estimated B0 x Depletion)
#' @examples
#' DBSRA(1, MSEtool::ourReefFish, plot=TRUE)
#'
#' @export
DBSRA <- function(x, Data, reps = 100, plot=FALSE) {
  runDBSRA <- DBSRA_(x,Data, reps)
  TAC <- MSEtool::TACfilter(runDBSRA[[1]])
  if (plot) DBSRA_plot(runDBSRA, Data, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DBSRA) <- "MP"


#' @templateVar mp DBSRA_40
#' @template MPuses
#' @describeIn DBSRA Same as the Base Version but assumes 40 percent current depletion (Bcurrent/B0 = 0.4), which is
#' more or less the most optimistic state for a stock (ie very close to BMSY/B0 for many stocks).
#' @examples
#' DBSRA_40(1, MSEtool::ourReefFish, plot=TRUE)
#' @export
DBSRA_40 <- function(x, Data, reps = 100, plot=FALSE) {
  runDBSRA <- DBSRA_(x,Data, reps, depo=0.4)
  TAC <- MSEtool::TACfilter(runDBSRA[[1]])
  if (plot)  DBSRA_plot(runDBSRA, Data, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}  # end of DBSRA_apply
class(DBSRA_40) <- "MP"


#' @templateVar mp DBSRA4010
#' @template MPuses
#' @describeIn DBSRA Base version paired with the 40-10 rule that linearly throttles
#' back the TAC when depletion is below 0.4 down to zero at 10 percent of unfished biomass.
#' @examples
#' DBSRA4010(1, MSEtool::ourReefFish, plot=TRUE)
#' @export
#'
DBSRA4010 <- function(x, Data, reps = 100, plot=FALSE) {
  runDBSRA <- DBSRA_(x,Data, reps, hcr=c(0.4, 0.1))
  TAC <- MSEtool::TACfilter(runDBSRA[[1]])
  if (plot) DBSRA_plot(runDBSRA, Data, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DBSRA4010) <- "MP"



# #' Depletion-Based Stock Reduction Analysis using mean length estimator of
# #' stock depletion
# #'
# #' DBSRA using the mean length estimator to calculate current stock depletion.
# #'
# #'
# #' @param x A position in a data-limited methods data object
# #' @param Data A data-limited methods data object
# #' @param reps The number of stochastic samples of the quota recommendation
# #' @note The mean length extension was programmed by Gary Nelson as part of his
# #' excellent R package 'fishmethods'
# #' @author T. Carruthers
# #' @references Dick, E.J., MacCall, A.D., 2011. Depletion-Based Stock Reduction
# #' Analysis: A catch-based method for determining sustainable yields for
# #' data-poor fish stocks. Fish. Res. 110, 331-341.
# #' @export
# DBSRA_ML <- function(x, Data, reps = 100) {
#   dependencies = "Data@Cat, Data@Mort, Data@CV_Mort, Data@FMSY_M, Data@CV_FMSY_M, Data@BMSY_B0, Data@CV_BMSY_B0, Data@L50, Data@CAL, Data@Year, Data@Cat"
#   C_hist <- Data@Cat[x, ]
#   TAC <- rep(NA, reps)
#   DBSRAcount <- 1
#   maxIts <- 200
#   nIts <- 0
#   if (is.na(Data@Dep[x]) | is.na(Data@CV_Dep[x])) return(NA)
#   while (DBSRAcount < (reps + 1) & nIts < maxIts) {
#     Linfc <- MSEtool::trlnorm(1, Data@vbLinf[x], Data@CV_vbLinf[x])
#     Kc <- MSEtool::trlnorm(1, Data@vbK[x], Data@CV_vbK[x])
#     Mdb <- MSEtool::trlnorm(100, Data@Mort[x], Data@CV_Mort[x])
#     Mdb <- Mdb[Mdb < 0.9][1]  # !!!! maximum M is 0.9   interval censor
#     if (is.na(Mdb)) Mdb <- 0.9  # !!!! maximum M is 0.9   absolute limit
#     Z <- MLne(x, Data, Linfc = Linfc, Kc = Kc, ML_reps = 1, MLtype = "dep")
#     if (all(is.na(Z))) {
#       Rec <- new("Rec")
#       Rec@TAC <- MSEtool::TACfilter(rep(NA, reps))
#       return(Rec)
#     }
#     FM <- Z - Mdb
#     FM[FM < 0] <- 0.01
#     nyears <- length(Data@Year)
#     Ct1 <- mean(Data@Cat[x, 1:3])
#     Ct2 <- mean(Data@Cat[x, (nyears - 2):nyears])
#     # dep<-c(Ct1,Ct2)/(1-exp(-FM[,c(1,2)]))
#     dep <- c(Ct1, Ct2)/(1 - exp(-FM))
#     Bt_K <- dep[2]/dep[1]
#
#     if (Bt_K < 0.01) Bt_K <- 0.01  # interval censor / temporary hack to avoid doing multiple depletion estimates that would take far too long
#     if (Bt_K > 0.99) Bt_K <- 0.99  # interval censor / temporary hack to avoid doing multiple depletion estimates that would take far too long
#
#
#     FMSY_M <- MSEtool::trlnorm(1, Data@FMSY_M[x], Data@CV_FMSY_M[x])
#     BMSY_K <- rbeta(100, MSEtool::alphaconv(Data@BMSY_B0[x], Data@CV_BMSY_B0[x] *
#                                      Data@BMSY_B0[x]), MSEtool::betaconv(Data@BMSY_B0[x], Data@CV_BMSY_B0[x] *
#                                                                   Data@BMSY_B0[x]))  #0.045 corresponds with mu=0.4 and quantile(BMSY_K,c(0.025,0.975)) =c(0.31,0.49) as in Dick and MacCall 2011
#     tryBMSY_K <- BMSY_K[BMSY_K > 0.05 & BMSY_K < 0.95][1]  # interval censor (0.05,0.95) as in Dick and MacCall, 2011
#
#     if (is.na(tryBMSY_K)) {
#       Min <- min(BMSY_K, na.rm = TRUE)
#       Max <- max(BMSY_K, na.rm = TRUE)
#       if (Max <= 0.05)
#         BMSY_K <- 0.05
#       if (Min >= 0.95)
#         BMSY_K <- 0.95
#     }
#     if (!is.na(tryBMSY_K))  BMSY_K <- tryBMSY_K
#     if (all(is.na(BMSY_K))) return(rep(NA, reps))
#     adelay <- max(floor(iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x],  Data@L50[x])), 1)
#     opt <- optimize(DBSRAopt, log(c(0.1 * mean(C_hist), 1000 * mean(C_hist))),
#                     C_hist = C_hist, nys = length(C_hist), Mdb = Mdb, FMSY_M = FMSY_M,
#                     BMSY_K = BMSY_K, Bt_K = Bt_K, adelay = adelay, tol = 0.01)
#     nIts <- nIts + 1
#     if (opt$objective < 0.1) {
#       Kc <- exp(opt$minimum)
#       BMSYc <- Kc * BMSY_K
#       FMSYc <- Mdb * FMSY_M
#       UMSYc <- (FMSYc/(FMSYc + Mdb)) * (1 - exp(-(FMSYc + Mdb)))
#       MSYc <- Kc * BMSY_K * UMSYc
#       TAC[DBSRAcount] <- UMSYc * Kc * Bt_K
#       DBSRAcount <- DBSRAcount + 1
#       nIts <- 0
#     }
#   }  # end of reps
#   Rec <- new("Rec")
#   Rec@TAC <- MSEtool::TACfilter(TAC)
#   Rec
# }
# class(DBSRA_ML) <- "MP"
#




## Delay-Difference MP ####

#' Delay-Difference Internal Function
#'
#' @param x iteration number
#' @param Data object of class Data
#' @param reps number of reps
#' @param hcr Optional harvest control rule for throttling catch as a function of B/B0.
#' Numeric vector of length 2 specifying HCR break points - e.g c(0.4, 0.1) for 40-10 HCR
#' @export
#' @keywords internal
DD_ <- function(x, Data, reps = 100, hcr=NULL) {

  Winf <- Data@wla[x] * Data@vbLinf[x]^Data@wlb[x]
  age <- 0:Data@MaxAge
  la <- Data@vbLinf[x] * (1 - exp(-Data@vbK[x] * ((age - Data@vbt0[x]))))
  wa <- Data@wla[x] * la^Data@wlb[x]
  a50V <- iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x],  Data@L50[x]) # assume vulnerability = maturity
  a50V <- max(a50V, 1)

  yind <- which(!is.na(Data@Cat[x, ] + Data@Ind[x, ]))[1] # First year with both catch and index
  yind <- yind:length(Data@Cat[x, ])
  Year <- Data@Year[yind]
  C_hist <- Data@Cat[x, yind]
  I_hist <- Data@Ind[x, yind]
  if(any(is.na(C_hist))) { # Linear interpolation of any missing catch
    C.xind <- 1:length(C_hist)
    C_hist <- approx(x = C.xind[!is.na(C_hist)], y = C_hist[!is.na(C_hist)], n = length(C.xind))$y
  }
  if(any(is.na(I_hist))) { # Linear interpolation of any missing index
    I.xind <- 1:length(I_hist)
    I_hist <- approx(x = I.xind[!is.na(I_hist)], y = I_hist[!is.na(I_hist)], n = length(I.xind))$y
  }
  E_hist <- C_hist/I_hist
  E_hist <- E_hist/mean(E_hist, na.rm=TRUE)
  ny_DD <- length(C_hist)
  k_DD <- ceiling(a50V)  # get age nearest to 50% vulnerability (ascending limb)
  k_DD[k_DD > Data@MaxAge/2] <- ceiling(Data@MaxAge/2)  # to stop stupidly high estimates of age at 50% vulnerability
  Rho_DD <- (wa[k_DD + 2] - Winf)/(wa[k_DD + 1] - Winf) # Brody growth coefficient
  Alpha_DD <- Winf * (1 - Rho_DD)
  So_DD <- exp(-Data@Mort[x])  # get So survival rate
  wa_DD <- wa[k_DD]
  UMSYpriorpar <- c(1 - exp(-Data@Mort[x] * 0.5), 0.3) # Prior for UMSY is that corresponding to F = 0.5 M with CV = 0.3
  UMSYprior <- c(MSEtool::alphaconv(UMSYpriorpar[1], prod(UMSYpriorpar)), MSEtool::betaconv(UMSYpriorpar[1], prod(UMSYpriorpar))) # Convert to beta parameters
  params <- log(c(UMSYpriorpar[1]/(1 - UMSYpriorpar[1]), 3*mean(C_hist, na.rm = T), Data@Mort[x]))
  opt <- optim(params, DD_R, opty = 1, So_DD = So_DD, Alpha_DD = Alpha_DD,
               Rho_DD = Rho_DD, ny_DD = ny_DD, k_DD = k_DD, wa_DD = wa_DD, E_hist = E_hist,
               C_hist = C_hist, UMSYprior = UMSYprior, method = "BFGS", hessian = TRUE)

  if (reps > 1) {
    if (requireNamespace("mvtnorm", quietly = TRUE)) {
      samps <- mvtnorm::rmvnorm(reps,opt$par,solve(opt$hessian)) # assuming log
    } else {
      warning('package `mvtnorm` needs to be installed for reps>1. Using MLE estimates')
      samps <- matrix(c(opt$par[1], opt$par[2], opt$par[3]), nrow = reps)
    }

  } else {
    samps <- matrix(c(opt$par[1], opt$par[2], opt$par[3]), nrow = 1)
  }

  #U_hist <- 1 - exp(-exp(opt$par[3]) * E_hist)
  UMSY <- 1/(1 + exp(-opt$par[1]))
  EMSY <- -log(1 - UMSY)/exp(opt$par[3]) # q in denominator
  eff <- EMSY/E_hist[Data@LHYear[1] - Year[1] + 1] # Effort advice is ratio of EMSY and obs. Eff in LHYear
  eff[!is.finite(eff)] <- 0.01
  eff[eff > 1e+05] <- 0.01

  # parameters are multivariate normal hessian approximation
  # samps <- cbind(rnorm(reps, opt$par[1], ((opt$par[1])^2)^0.5 * 0.1),
  #                rnorm(reps, opt$par[2], ((opt$par[2])^2)^0.5 * 0.1),
  #                rnorm(reps, opt$par[3], ((opt$par[3])^2)^0.5 * 0.1))

  getVals <- sapply(1:reps, function(i) DD_R(params=samps[i, ], opty = 2, So_DD = So_DD,
                                             Alpha_DD = Alpha_DD, Rho_DD = Rho_DD, ny_DD = ny_DD, k_DD = k_DD,
                                             wa_DD = wa_DD, E_hist = E_hist, C_hist = C_hist, UMSYprior = UMSYprior))

  TAC <- unlist(getVals[1,])
  dep <- unlist(getVals[2,])
  Cpredict <- do.call("cbind", getVals[3,])
  B_DD <- do.call("cbind", getVals[4,])

  if (!is.null(hcr)) {
    cond1 <- !is.na(dep) & dep < hcr[1] & dep > hcr[2]
    cond2 <- !is.na(dep) & dep < hcr[2]
    TAC[cond1] <- TAC[cond1] * (dep[cond1] - hcr[2])/(hcr[1]-hcr[2])
    TAC[cond2] <- TAC[cond2] * tiny  # this has to still be stochastic albeit very small
  }


  return(list(TAC=TAC, C_hist=C_hist, I_hist=I_hist, E_hist=E_hist, dep=dep,
              Cpredict=Cpredict, B_DD=B_DD, hcr=hcr, Year=Data@Year[yind], eff=eff))
}


#' Delay-Difference Internal Function
#'
#' @param params internal parameter
#' @param opty  internal parameter
#' @param So_DD  internal parameter
#' @param Alpha_DD  internal parameter
#' @param Rho_DD  internal parameter
#' @param ny_DD  internal parameter
#' @param k_DD internal parameter
#' @param wa_DD  internal parameter
#' @param E_hist  internal parameter
#' @param C_hist internal parameter
#' @param UMSYprior internal parameter
#'
#' @export
#'
#' @keywords internal
DD_R <- function(params, opty, So_DD, Alpha_DD, Rho_DD, ny_DD, k_DD, wa_DD, E_hist,
                 C_hist, UMSYprior) {
  UMSY_DD = 1/(1 + exp(-params[1])) # Logit transform to constrain u between 0-1
  MSY_DD = exp(params[2])
  q_DD = exp(params[3])
  SS_DD = So_DD * (1 - UMSY_DD)  # Initialise for UMSY, MSY and q leading.
  Spr_DD = (SS_DD * Alpha_DD/(1 - SS_DD) + wa_DD)/(1 - Rho_DD * SS_DD)
  DsprDu_DD = ((Alpha_DD + Spr_DD * (1 + Rho_DD - 2 * Rho_DD * SS_DD))/((1 - Rho_DD * SS_DD) * (1 - SS_DD)) +
                 Alpha_DD * SS_DD/((1 - Rho_DD * SS_DD) * (1 - SS_DD)^2) - Spr_DD/(1 - SS_DD)) * -So_DD
  Arec_DD = 1/(((1 - UMSY_DD)^2) * (Spr_DD + UMSY_DD * DsprDu_DD))
  Brec_DD = UMSY_DD * (Arec_DD * Spr_DD - 1/(1 - UMSY_DD))/MSY_DD
  Spr0_DD = (So_DD * Alpha_DD/(1 - So_DD) + wa_DD)/(1 - Rho_DD * So_DD)
  Ro_DD = (Arec_DD * Spr0_DD - 1)/(Brec_DD * Spr0_DD)
  Bo_DD = Ro_DD * Spr0_DD
  No_DD = Ro_DD/(1 - So_DD)

  B_DD <- rep(NA, ny_DD + 1)
  N_DD <- rep(NA, ny_DD + 1)
  R_DD <- rep(NA, ny_DD + k_DD)
  Cpred_DD <- rep(NA, ny_DD)

  B_DD[1] = Bo_DD
  N_DD[1] = No_DD
  R_DD[1:k_DD] = Ro_DD

  for (tt in 1:ny_DD) {
    Surv_DD = So_DD * exp(-q_DD * E_hist[tt])
    Cpred_DD[tt] = B_DD[tt] * (1 - exp(-q_DD * E_hist[tt]))
    Sp_DD = B_DD[tt] - Cpred_DD[tt]
    R_DD[tt + k_DD] = Arec_DD * Sp_DD/(1 + Brec_DD * Sp_DD)
    B_DD[tt + 1] = Surv_DD * (Alpha_DD * N_DD[tt] + Rho_DD * B_DD[tt]) +
      wa_DD * R_DD[tt + 1]
    N_DD[tt + 1] = Surv_DD * N_DD[tt] + R_DD[tt + 1]
  }
  tiny <- 1e-15
  Cpred_DD[Cpred_DD < tiny] <- tiny

  if (opty == 1) {
    # The following conditions must be met for positive values
    # of Arec_DD and Brec_DD, respectively:
    # umsy * DsprDu + Spr_DD > 0 and Arec_DD * Spr_DD * (1 - UMSY_DD) - 1 > 0
    # Thus, create a likelihood penalty of 100 if either condition is not met
    umsy_penalty <- ifelse(Spr_DD + UMSY_DD * DsprDu_DD > 0, 0, UMSY_DD * 100)
    alpha_penalty <- ifelse(Arec_DD * Spr_DD * (1 - UMSY_DD) - 1 > 0, 0, UMSY_DD * 100)

    sigma <- sqrt(sum((log(C_hist) - log(Cpred_DD))^2)/ny_DD) # Analytical solution

    test <- dnorm(log(C_hist), log(Cpred_DD), sigma, log = T)

    test2 <- dbeta(UMSY_DD, UMSYprior[1], UMSYprior[2], log = T)
    test[is.na(test)] <- -1000
    test[test == (-Inf)] <- -1000
    if (is.na(test2) | test2 == -Inf | test2 == Inf)  test2 <- 1000
    return(-sum(test, test2) + umsy_penalty + alpha_penalty)  # return objective function
  } else {
    return(list(TAC= UMSY_DD * B_DD[ny_DD + 1], dep=B_DD[tt + 1]/Bo_DD,
                Cpred_DD=Cpred_DD,B_DD=B_DD))
  }
  # } if (opty == 2) {
  #   # return MLE TAC estimate
  #
  # } else if (opty == 3) {
  #
  # } else {
  #   cbind(C_hist, Cpred_DD)  # return observations vs predictions
  # }
}

#' Delay - Difference Stock Assessment
#'
#' A simple delay-difference assessment with UMSY and MSY as leading parameters
#' that estimates the TAC using a time-series of catches and a relative
#'
#' This DD model is observation error only and has does not estimate
#' process error (recruitment deviations). Assumption is that knife-edge
#' selectivity occurs at the age of 50% maturity. Similar to many other assessment
#' models it depends on a whole host of dubious assumptions such as temporally
#' stationary productivity and proportionality between the abundance index and
#' real abundance. Unsurprisingly the extent to which these assumptions are
#' violated tends to be the biggest driver of performance for this method.
#'
#' The method is conditioned on effort and estimates catch. The effort is calculated
#' as the ratio of catch and index. Thus, to get a complete effort time series, a full
#' time series of catch and index is also needed. Missing values are linearly interpolated.
#'
#' A detailed description of the delay-difference model can be found in Chapter 9 of Hilborn
#' and Walters (1992).
#'
#' @templateVar mp DD
#' @template MPtemplate
#' @template MPuses
#'
#' @author T. Carruthers
#' @references
#' Carruthers, T, Walters, C.J,, and McAllister, M.K. 2012. Evaluating methods that classify
#' fisheries stock status using only fisheries catch data. Fisheries Research 119-120:66-79.
#'
#' Hilborn, R., and Walters, C. 1992. Quantitative Fisheries Stock Assessment: Choice,
#' Dynamics and Uncertainty. Chapman and Hall, New York.
#' @describeIn DD Base version where the TAC = UMSY * Current Biomass.
#'
#' @examples
#' DD(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
#' @family Delay-Difference MPs
#'
#' @export
DD <- function(x, Data, reps = 100, plot=FALSE) {
  runDD <- DD_(x, Data, reps)
  TAC <- MSEtool::TACfilter(runDD$TAC)

  if (plot) DD_plot(x, runDD, Data, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DD) <- "MP"



#' @templateVar mp DD4010
#' @template MPuses
#' @examples
#' DD4010(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
#' @describeIn DD A 40-10 rule is imposed over the TAC recommendation.
#' @export
DD4010 <- function(x, Data, reps = 100, plot=FALSE) {

  runDD <- DD_(x, Data, reps, hcr=c(.4, .1))
  TAC <- MSEtool::TACfilter(runDD$TAC)

  if (plot) DD_plot(x, runDD, Data, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DD4010) <- "MP"







## Dynamic Fratio ####






#' Dynamic Fratio MP
#'
#' The Fratio MP with a controller that changes the level of F according to the
#' estimated relationship between surplus production and biomass. Ie lower F when dSP/dB
#' is positive and higher F when dSP/dB is negative.
#'
#' The method smoothes historical catches and biomass and then infers the
#' relationship between surplus production and biomass (as suggested by Mark
#' Maunder and Carl Walters). The approach then regulates a F based policy
#' according to this gradient in which F may range between two different
#' fractions of natural mortality rate.
#'
#' The core advantage is the TAC(t) is not strongly determined by TAC(t-1) and
#' therefore errors are not as readily propagated. The result is method that
#' tends to perform alarmingly well and therefore requires debunking ASAP.
#'
#' The catch limit (TAC) is calculated as:
#' \deqn{\textrm{TAC}=F  B}
#' where \eqn{F} is fishing mortality and \eqn{B} is the estimated current biomass.
#'
#' \eqn{F} is calculated as:
#' \deqn{F = F_{\textrm{MSY}} \exp{-gG}}
#' where \eqn{F_{\textrm{MSY}}} is calculated from assumed values of \eqn{\frac{F_{\textrm{MSY}}}{M}} and
#' \eqn{M}, *g* is a gain parameter and *G* is the estimated gradient in surplus
#' production (*SP*) as a function of biomass (*B*). Surplus production for year *y* is calculated as:
#' \deqn{SP_y = B_{y+1} - B_y + C_y}
#' Trends in historical catch (*C*) and biomass (*B*) are both estimated using a loess smoother, over the last `yrsmth` years,
#' of available catch and a time-series of abundance, calculated from an index of abundance (`Data@Ind`)
#' and an estimate of abundance (`Data@Abun`) for the current year.
#'
#'
#' @templateVar mp DynF
#' @template MPtemplate
#' @template MPuses
#'
#' @param gg A gain parameter that modifies F according to the gradient in
#' surplus production with biomass
#' @param yrsmth The number of historical recent years used for smoothing catch
#' and biomass data
#'
#' @author T. Carruthers
#' @references Made-up for this package.
#' @family Fmsy/M methods
#' @examples
#' DynF(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
#' @export
DynF <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2) {

  dependencies = "Data@Year, Data@Cat, Data@Ind, Data@Abun, Data@Mort, Data@FMSY_M"
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  years <- Data@Year[ind]
  C_dat <- log(Data@Cat[x, ind])
  C_dat[C_dat == -Inf] <- 0
  B_dat <- log(Data@Ind[x, ind]/Data@Ind[x, ind[yrsmth]] * Data@Abun[x])
  B_dat[B_dat == -Inf] <- 0
  C_hist <- exp(predict(loess(C_dat ~ ind, degree = 1)))
  B_hist <- exp(predict(loess(B_dat ~ ind, degree = 1)))


  ind <- 2:yrsmth
  ind1 <- 1:(yrsmth - 1)
  SP_hist <- B_hist[ind] - B_hist[ind1] + C_hist[ind1]

  Frat <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x]) *
    MSEtool::trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x])
  Flim <- matrix(NA, nrow = 2, ncol = reps)
  Flim[1, ] <- Frat * 0.5
  Flim[2, ] <- Frat * 2

  # yind <- 1:length(SP_hist)
  # SP_mu <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1))
  # SP_se <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1), se = T)$se.fit
  # SP_new <- rnorm(reps, SP_mu, SP_se/2)
  Glm <- summary(lm(SP_hist ~ B_hist[ind1]))$coefficients[2, 1:2]  # gradient in SP as function of biomass
  G_new <- rnorm(reps, Glm[1], Glm[2]/2) # sample gradients
  # G_new[G_new>2*Frat]<-2*Frat[G_new<(2*Frat)]
  # G_new[G_new<(-2*Frat)]<--2*Frat[G_new<(-2*Frat)]
  G_new[G_new > 0] <- G_new[G_new > 0] * 3
  newF <- Frat * exp(-G_new * gg)
  newF[newF < Flim[1]] <- Flim[1]
  newF[newF > Flim[2]] <- Flim[2]

  TAC <- MSEtool::TACfilter(newF * B_hist[yrsmth])

  if (plot) DynF_plot(C_dat, C_hist, TAC, yrsmth, B_dat, B_hist, Data, SP_hist,
                      ind, ind1, G_new, Frat,newF, years)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec

}
class(DynF) <- "MP"


#' Adaptive Fratio
#'
#' An adaptive MP that uses trajectory in inferred suplus production and
#' fishing mortality rate to update a TAC
#'
#' Fishing rate is modified each year according to the gradient of surplus
#' production with biomass (aims for zero).  F is bounded by FMSY/2 and 2FMSY
#' and walks in the logit space according to dSP/dB. This is derived from the
#' theory of Maunder 2014.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC}_y= F_y B_{y-1}}
#' where \eqn{B_{y-1}} is the most recent biomass, estimated with a loess smoother
#' of the most recent `yrsmth` years from the index of abundance (`Data@Ind`)
#' and estimate of current abundance (`Data@Abun`), and
#'
#' \deqn{F_y = F_{\textrm{lim}_1} + \left(\frac{\exp^{F_{\textrm{mod}_2}}} {1 + \exp^{F_{\textrm{mod}_2}}} F_{\textrm{lim}_3} \right)   }
#'
#' where \eqn{F_{\textrm{lim}_1} = 0.5 \frac{F_\textrm{MSY}}{M}M},
#' \eqn{F_{\textrm{lim}_2} = 2 \frac{F_\textrm{MSY}}{M}M}, \eqn{F_{\textrm{lim}_3}} is \eqn{F_{\textrm{lim}_2} - F_{\textrm{lim}_1}},
#' \eqn{F_{\textrm{mod}_2}} is
#' \deqn{F_{\textrm{mod}_1} + g -G}
#' where \eqn{g} is gain parameter `gg`, `G` is the predicted surplus production given current abundance,
#' and:
#'   \deqn{F_{\textrm{mod}_1} =
#'              \left\{\begin{array}{ll}
#'              -2 & \textrm{if } F_\textrm{old} < F_{\textrm{lim}_1} \\
#'              2 & \textrm{if }  F_\textrm{old} > F_{\textrm{lim}_2} \\
#'              \log{\frac{F_\textrm{frac}}{1-F_\textrm{frac}}} & \textrm{if } F_{\textrm{lim}_1} \leq F_\textrm{old} \leq F_{\textrm{lim}_2} \\
#'              \end{array}\right.
#'            }
#' where \eqn{-F_{\textrm{frac}} = \frac{F_{\textrm{old}} - F_{\textrm{lim}_1}}{F_{\textrm{lim}_3}} },
#' \eqn{F_\textrm{old} = \sum{\frac{C_\textrm{hist}}{B_\textrm{hist}}}/n}
#' where \eqn{C_\textrm{hist}} and \eqn{B_\textrm{hist}} are smooth catch and biomass over last `yrsmth`,
#' and \eqn{n} is `yrsmth`.
#'
#'
#'
#' Tested in Carruthers et al. 2015.
#'
#' @templateVar mp Fadapt
#' @template MPtemplate
#' @template MPuses
#' @param yrsmth Years over which to smooth recent estimates of surplus
#' production
#' @param gg A gain parameter controlling the speed in update in TAC.
#' @return A numeric vector of quota recommendations
#' @author T. Carruthers
#' @references Carruthers et al. 2015. Performance evaluation of simple
#' management procedures. ICES J. Mar Sci. 73, 464-482.
#'
#' Maunder, M. 2014.
#' http://www.iattc.org/Meetings/Meetings2014/MAYSAC/PDFs/SAC-05-10b-Management-Strategy-Evaluation.pdf
#' @family Fmsy/M methods
#' @family Surplus production MPs
#' @examples
#' Fadapt(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
#' @export
Fadapt <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 7, gg = 1) {

  dependencies = "Data@Year, Data@Cat, Data@Ind, Data@Abun, Data@Mort, Data@FMSY_M"
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  years <- Data@Year[ind]
  C_dat <- log(Data@Cat[x, ind])
  C_dat[C_dat == -Inf] <- 0
  B_dat <- log(Data@Ind[x, ind]/Data@Ind[x, ind[yrsmth]] * Data@Abun[x])
  B_dat[B_dat == -Inf] <- 0
  C_hist <- exp(predict(loess(C_dat ~ ind, degree = 1)))
  B_hist <- exp(predict(loess(B_dat ~ ind, degree = 1)))

  ind <- 2:yrsmth
  ind1 <- 1:(yrsmth - 1)
  SP_hist <- B_hist[ind] - B_hist[ind1] + C_hist[ind1]

  Frat <- Data@Mort[x] * Data@FMSY_M[x]
  Flim <- Frat * c(0.5, 2)
  Flimr <- Flim[2] - Flim[1]

  yind <- 1:length(SP_hist)
  SP_mu <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1))
  SP_se <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1), se = T)$se.fit
  SP_new <- rnorm(reps, SP_mu, SP_se/2)
  Glm <- summary(lm(SP_hist ~ B_hist[ind1]))$coefficients[2, 1:2]  # plot(B_hist[ind1],SP_hist) # points(B_hist[ind1],SP_hist,col='green')
  if (reps > 1) {
    G_new <- rnorm(reps, Glm[1], Glm[2])
  } else {
    G_new <- Glm[1]
  }

  Fold <- mean(C_hist/B_hist, na.rm=TRUE)

  if (Fold < Flim[1]) Fmod1 <- (-2)
  if (Fold > Flim[2]) Fmod1 <- 2
  if (Fold > Flim[1] & Fold < Flim[2]) {
    Ffrac <- (Fold - Flim[1])/Flimr
    Fmod1 <- log(Ffrac/(1 - Ffrac))
  }
  Fmod2 <- Fmod1 + gg * -G_new
  newF <- Flim[1] + (exp(Fmod2)/(1 + exp(Fmod2))) * Flimr
  TAC <- MSEtool::TACfilter(newF * B_hist[yrsmth])

  if (plot) Fadapt_plot(C_dat, C_hist, TAC, yrsmth, B_dat, B_hist, Data, SP_hist,
                      ind, ind1, G_new, Frat,newF, years)


  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fadapt) <- "MP"



#' Internal function to estimate r
#'
#' @param x iteration number
#' @param Data Object of class `Data`
#' @param Mvec M
#' @param Kvec K
#' @param Linfvec Linf
#' @param t0vec t0
#' @param hvec steepness
#' @param maxage maximum age
#' @param r_reps number of reps
#'
#' @export
#'
#' @keywords internal
getr <- function(x, Data, Mvec, Kvec, Linfvec, t0vec, hvec, maxage,
                 r_reps = 100) {
  r <- rep(NA, r_reps)
  for (i in 1:r_reps) {
    log.r <- log(0.3)
    opt <- optimize(demofn, lower = log(1e-04), upper = log(1.4), M = Mvec[i],
                   amat = iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x],
                              Data@L50[x]), sigma = 0.2, K = Kvec[i], Linf = Linfvec[i],
                   to = t0vec[i], hR = hvec[i], maxage = maxage, a = Data@wla[x],
                   b = Data@wlb[x])
    # demographic2(opt$minimum,M[x],ageM[x],0.2,K[x],Linf,t0,steepness[x],maxage,a,b)$r
    r[i] <- exp(opt$minimum)
  }
  r
}



#' Internal demographic function
#'
#' @param log.r internal
#' @param M internal
#' @param amat internal
#' @param sigma internal
#' @param K internal
#' @param Linf internal
#' @param to internal
#' @param hR internal
#' @param maxage internal
#' @param a internal
#' @param b internal
#'
#' @export
#'
#' @keywords internal
demographic2 <- function(log.r, M, amat, sigma, K, Linf, to, hR, maxage, a, b) {
  # switch on and off to use either S or m in MC simulations
  r = exp(log.r)
  lx = exp(-M)^((1:maxage) - 1)  #survivorship
  logNormDensity = (dnorm(x = log((1:maxage)), mean = log(amat), sd = sigma))/(1:maxage)  #Maturity ogive calculation
  logNormDensity[1] = 0
  sumlogNormDen = sum(logNormDensity)
  NormalisedMaturity = logNormDensity/sumlogNormDen
  proportionMat <- vector()
  proportionMat[1] = NormalisedMaturity[1]
  for (i in 2:maxage) proportionMat[i] = proportionMat[i - 1] + NormalisedMaturity[i]
  TL = Linf * (1 - exp(-K * ((1:maxage) - to)))  #length at age
  Wa = a * TL^b  #wegith at age
  SurvWeiMat = lx * Wa * proportionMat  #survivorship X weight X maturity
  SBPR = sum(SurvWeiMat)  #Spawner biomass per recruit
  RPS = 1/(SBPR * (1 - hR)/(4 * hR))  # Beverton Holt
  # RPS=(5*hR)^(5/4)/SBPR # Ricker Recruitment per spawner biomass
  RPF = Wa * proportionMat * RPS  #Recruits per female
  Lotka = lx * RPF * exp(-(1:maxage) * r)
  sumLotka = sum(Lotka)
  epsilon = (1 - sumLotka)^2  #objective function
  return(list(epsilon = epsilon, r = r))
}

#' @describeIn demographic2 Internal function
#' @export
demofn <- function(log.r, M, amat, sigma, K, Linf, to, hR, maxage, a, b) {
  demographic2(log.r, M, amat, sigma, K, Linf, to, hR, maxage = maxage, a, b)$epsilon
}


#' Internal function for Fdem MP
#'
#' @param x internal
#' @param Data internal
#' @param reps internal
#' @param Ac Optional abundance value
#'
#' @export
#' @keywords internal
Fdem_ <- function(x, Data, reps, Ac=NULL) {
  # Demographic FMSY estimate (FMSY=r/2)
  if (MSEtool::NAor0(Data@Mort[x])) stop("Data@Mort is NA")
  if (MSEtool::NAor0(Data@FMSY_M[x])) stop("Data@FMSY_M is NA")
  if (MSEtool::NAor0(Data@L50[x])) stop("Data@L50 is NA")
  if (MSEtool::NAor0(Data@vbK[x])) stop("Data@vbK is NA")
  if (MSEtool::NAor0(Data@CV_vbK[x])) stop("Data@CV_vbK is NA")
  if (MSEtool::NAor0(Data@vbLinf[x])) stop("Data@vbLinf is NA")
  if (MSEtool::NAor0(Data@CV_vbLinf[x])) stop("Data@CV_vbLinf is NA")
  if (MSEtool::NAor0(Data@vbt0[x])) stop("Data@vbt0 is NA")
  if (MSEtool::NAor0(Data@CV_Mort[x])) stop("Data@CV_Mort is NA")
  if (MSEtool::NAor0(Data@wla[x])) stop("Data@wla is NA")
  if (MSEtool::NAor0(Data@wlb[x])) stop("Data@wlb is NA")

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kc <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfc = MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    if (Data@vbt0[x]<0) {
      t0c <- -MSEtool::trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
    } else {
      t0c <- MSEtool::trlnorm(reps, Data@vbt0[x], Data@CV_vbt0[x])
    }
  } else {
    t0c <- rep(Data@vbt0[x], reps)
  }
  t0c[!is.finite(t0c)] <- 0
  # hvec <- MSEtool::trlnorm(reps, Data@steep[x], Data@CV_steep[x])
  hvec <- MSEtool::sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])

  if (is.null(Ac)) Ac <- MSEtool::trlnorm(reps, Data@Abun[x], Data@CV_Abun[x])
  FMSY <- getr(x, Data, Mvec, Kc, Linfc, t0c, hvec, maxage = Data@MaxAge, r_reps = reps)/2
  TAC <- FMSY * Ac
  return(list(TAC=TAC, Ac=Ac, FMSY=FMSY))
}


#' Demographic FMSY method
#'
#' FMSY is calculated as r/2 where r is calculated from a demographic approach
#' (inc steepness). Coupled with an estimate of current abundance that gives
#' you the OFL.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = F_{\textrm{MSY}} A}
#' where *A* is an estimate of current abundance, and \eqn{F_{\textrm{MSY}}} is estimated
#' as \eqn{r/2}, where \eqn{r} is the intrinsic rate of population growth, estimated
#' from the life-history parameters using the methods of McAllister et al. (2001).
#'
#' @templateVar mp Fdem
#' @template MPtemplate
#' @template MPuses
#'
#' @author T. Carruthers
#' @references McAllister, M.K., Pikitch, E.K., and Babcock, E.A. 2001. Using
#' demographic methods to construct Bayesian priors for the intrinsic rate of
#' increase in the Schaefer model and implications for stock rebuilding. Can.
#' J. Fish. Aquat. Sci. 58: 1871-1890.
#' @describeIn Fdem Current abundance is assumed to be known (i.e `Data@Abun`)
#' @examples
#' Fdem(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
Fdem <- function(x, Data, reps = 100, plot=FALSE) {
  runFdem <- Fdem_(x, Data, reps)
  TAC <- MSEtool::TACfilter(runFdem$TAC)

  if (plot) Fdem_plot(runFdem, Data)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fdem) <- "MP"




#' @templateVar mp Fdem_CC
#' @template MPuses
#'
#' @param Fmin The minimum fishing mortality rate derived from the catch-curve
#' analysis
#'
#' @describeIn Fdem Current abundance is estimated from catch curve analysis
#' @examples
#' Fdem_CC(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
Fdem_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {

  Mvec <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, MuC, Data@CV_Cat[x,1])
  Zdb <- CC(x, Data, reps = reps * 10)
  Fdb <- Zdb - Mvec
  ind <- (1:(reps * 10))[Fdb > Fmin][1:reps]

  Fdb <- Fdb[ind]
  SM <- sum(is.na(ind))
  if (SM > 0) {
    Fdb[is.na(ind)] <- Fmin
  }

  Ac <- Cc/(1 - exp(-Fdb))

  runFdem <- Fdem_(x, Data, reps, Ac=Ac)
  TAC <- MSEtool::TACfilter(runFdem$TAC)

  if (plot) Fdem_plot(runFdem, Data)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fdem_CC) <- "MP"

#' @templateVar mp Fdem_ML
#' @template MPuses
#' @describeIn Fdem Current abundance is estimated from mean length
#' @examples
#' Fdem_ML(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
Fdem_ML <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  Mvec <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  Kc <- MSEtool::trlnorm(reps * 10, Data@vbK[x], Data@CV_vbK[x])
  Linfc = MSEtool::trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps * 10, MuC, Data@CV_Cat[x,1])
  Z <- MLne(x, Data, Linfc = Linfc, Kc = Kc, ML_reps = reps * 10, MLtype = "F")

  if (all(is.na(Z))) {
    Rec <- new("Rec")
    Rec@TAC <- rep(as.numeric(NA), reps)
    return(Rec)
  }

  ind <- !is.na(Z)
  FM <- Z[ind] - Mvec[ind]
  ind <- (1:(reps * 10))[FM > Fmin][1:reps]
  FM <- FM[ind]
  SM <- sum(is.na(ind))
  if (SM > 0) {
    FM[is.na(ind)] <- Fmin
  }

  Ac <- Cc[ind]/(1 - exp(-FM))

  runFdem <- Fdem_(x, Data, reps, Ac=Ac)
  TAC <- MSEtool::TACfilter(runFdem$TAC)

  if (plot) Fdem_plot(runFdem, Data)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fdem_ML) <- "MP"



##### Fratio methods ####
#' Fratio internal function
#'
#' @param x Iteration number
#' @param Data Object of class Data
#' @param reps Number of reps
#' @param Abun Optional estimate of abundance
#'
#' @export
#'
#' @keywords internal
Fratio_ <- function(x, Data, reps=100, Abun=NULL) {
  Frat <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x]) * MSEtool::trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x]) # estimate of Fmsy
  if (is.null(Abun)) Abun <- MSEtool::trlnorm(reps, Data@Abun[x], Data@CV_Abun[x])

  TAC <- Frat * Abun

  return(list(TAC=TAC, Abun=Abun, Frat=Frat))

}


#' FMSY/M ratio methods
#'
#' Calculates the OFL based on a fixed ratio of FMSY to M multiplied by a
#' current estimate of abundance.
#'
#' A simple method that tends to outperform many other approaches alarmingly
#' often even when current biomass is relatively poorly known. The low stock
#' crash potential is largely due to the quite large difference between Fmax
#' and FMSY for most stocks.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = F_{\textrm{MSY}} A}
#' where \eqn{F_{\textrm{MSY}}} is calculated as \eqn{\frac{F_\textrm{MSY}}{M} M} and
#' *A* is estimate of current abundance.
#'
#' The MP variants differ in the assumption of current abundance (see Functions section below)
#'
#' @templateVar mp Fratio
#' @template MPtemplate
#' @template MPuses
#'
#' @author T. Carruthers
#' @references Gulland, J.A., 1971. The fish resources of the ocean. Fishing
#' News Books, West Byfleet, UK.
#'
#' Martell, S., Froese, R., 2012. A simple method for estimating MSY from catch
#' and resilience. Fish Fish. doi: 10.1111/j.1467-2979.2012.00485.x.
#'
#' @describeIn Fratio Requires an estimate of current abundance (i.e `Data@Abun`)
#' @family Fmsy/M methods
#' @examples
#' Fratio(1, MSEtool::Atlantic_mackerel, plot=TRUE)
#' @export
Fratio <- function(x, Data, reps = 100, plot=FALSE) {
  # FMSY / M ratio method e.g. Gulland

  runFrat <- Fratio_(x, Data, reps)
  TAC <- MSEtool::TACfilter(runFrat$TAC)

  if (plot) Fratio_plot(x, Data, TAC, runFrat)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fratio) <- "MP"

#' @templateVar mp Fratio4010
#' @template MPuses
#' @describeIn Fratio Paired with the 40-10 rule that throttles back the OFL to zero at 10 percent of unfished biomass. Requires an estimate of current depletion.
#' @examples
#' Fratio4010(1, MSEtool::Atlantic_mackerel, plot=TRUE)
#' @export
Fratio4010 <- function(x, Data, reps = 100, plot=FALSE) {
  runFrat <- Fratio_(x, Data, reps)
  TAC <- MSEtool::TACfilter(runFrat$TAC)

  # 40-10 rule
  # Bt_K <- MSEtool::trlnorm(reps, Data@Dt[x], Data@CV_Dt[x])
  if (is.na(Data@Dep[x]) | is.na(Data@CV_Dep[x])) {
    out <- new("Rec")
    out@TAC <- rep(as.numeric(NA), reps)
    return(out)
  }
  depo <- max(0.01, min(0.99, Data@Dep[x]))
  Bt_K <- rbeta(reps * 100, MSEtool::alphaconv(depo, min(depo * Data@CV_Dep[x],
                                                (1 - depo) * Data@CV_Dep[x])),
                MSEtool::betaconv(depo, min(depo * Data@CV_Dep[x], (1 - depo) * Data@CV_Dep[x])))
  Bt_K <- Bt_K[Bt_K >= 0.01 & Bt_K <= 0.99][1:reps]  # interval censor (0.01,0.99)  as in Dick and MacCall 2011

  cond1 <- Bt_K < 0.4 & Bt_K > 0.1
  cond2 <- Bt_K < 0.1
  TAC[cond1] <- TAC[cond1] * (Bt_K[cond1] - 0.1)/0.3
  TAC[cond2] <- TAC[cond2] * tiny  # this has to still be stochastic albeit very small

  runFrat$Bt_K <- Bt_K
  if (plot) Fratio_plot(x, Data, TAC, runFrat)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fratio4010) <- "MP"

#' @templateVar mp DepF
#' @template MPuses
#'
#' @describeIn Fratio Depletion Corrected Fratio: the Fratio MP with a harvest control
#' rule that reduces F according to the production curve given an estimate of current
#' stock depletion (made-up for this package).
#' @export
DepF <- function(x, Data, reps = 100, plot=FALSE) {

  runFrat <- Fratio_(x, Data, reps)
  TAC <- MSEtool::TACfilter(runFrat$TAC)

  if (is.na(Data@Dep[x]) | is.na(Data@CV_Dep[x])) {
    out <- new("Rec")
    out@TAC <- rep(as.numeric(NA), reps)
    return(out)
  }

  depo <- max(0.01, min(0.99, Data@Dep[x]))
  Bt_K <- rbeta(reps * 100, MSEtool::alphaconv(depo, min(depo * Data@CV_Dep[x],
                                                (1 - depo) * Data@CV_Dep[x])),
                MSEtool::betaconv(depo, min(depo * Data@CV_Dep[x], (1 - depo) * Data@CV_Dep[x])))
  Bt_K <- Bt_K[Bt_K >= 0.01 & Bt_K <= 0.99][1:reps]  # interval censor (0.01,0.99)  as in Dick and MacCall 2011


  adj <- Bt_K * (1 - Bt_K) * 4
  adj[Bt_K > 0.5] <- 1

  TAC <- TAC * adj
  runFrat$Bt_K <- Bt_K

  if (plot) Fratio_plot(x, Data, TAC, runFrat)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(DepF) <- "MP"


#' @templateVar mp Fratio_CC
#' @template MPuses
#' @param Fmin Minimum current fishing mortality rate for the catch-curve analysis
#' @describeIn Fratio Current abundance is estimated using average catch and estimate of F from an age-based catch curve
#' @examples
#' Fratio_CC(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
Fratio_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  # estimate abundance from average catch and F
  # MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  MuC <- mean(Data@Cat[x, ], na.rm=TRUE)
  Cc <- MSEtool::trlnorm(reps, MuC, Data@CV_Cat[x,1])
  Mdb <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])  # CV of 0.5 as in MacCall 2009
  Zdb <- CC(x, Data, reps = reps * 10)
  Fdb <- Zdb - Mdb
  ind <- (1:(reps * 10))[Fdb > 0.005][1:reps]

  Fdb <- Fdb[ind]
  Mdb <- Mdb[ind]
  SM <- sum(is.na(ind))
  if (SM > 0) {
    Mdb[is.na(ind)] <- MSEtool::trlnorm(SM, Data@Mort[x], Data@CV_Mort[x])
    Fdb[is.na(ind)] <- Fmin
  }

  Ac <- Cc/(1 - exp(-Fdb))
  runFrat <- Fratio_(x, Data, reps, Abun=Ac)

  TAC <- MSEtool::TACfilter(runFrat$TAC)

  if (plot) Fratio_plot(x, Data, TAC, runFrat)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fratio_CC) <- "MP"

#' @templateVar mp Fratio_ML
#' @template MPuses
#' @describeIn Fratio Current abundance is estimated using average catch and estimate of F from mean lengths
#' @examples
#' Fratio_ML(1, MSEtool::SimulatedData, plot=TRUE)
#'
#' @export
Fratio_ML <- function(x, Data, reps = 100, plot=FALSE) {
  # estimate abundance from average catch and F
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps * 10, MuC, Data@CV_Cat[x,1])
  Mdb <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])  # CV of 0.5 as in MacCall 2009
  Linfc <- MSEtool::trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- MSEtool::trlnorm(reps * 10, Data@vbK[x], Data@CV_vbK[x])
  Z <- MLne(x, Data, Linfc = Linfc, Kc = Kc, ML_reps = reps * 10, MLtype = "F")
  if (all(is.na(Z))) {
    out <- new("Rec")
    out@TAC <- rep(as.numeric(NA), reps)
    return(out)
  }

  FM <- Z - Mdb
  ind <- which(FM>0)[1:reps]
  FM <- FM[ind]
  Ac <- Cc[ind]/(1 - exp(-FM))

  runFrat <- Fratio_(x, Data, reps, Abun=Ac)
  TAC <- MSEtool::TACfilter(runFrat$TAC)

  if (plot) Fratio_plot(x, Data, TAC, runFrat)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Fratio_ML) <- "MP"



### Geromont and Butterworth index slope Harvest Control Rule ####

#' Geromont and Butterworth index slope Harvest Control Rule
#'
#' An MP similar to SBT1 that modifies a time-series of catch recommendations
#' and aims for a stable catch rates.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC}_y= C_{y-1} \left(1+\lambda I\right)}
#' where \eqn{C_{y-1}} is catch from the previous year, \eqn{\lambda} is a gain parameter, and \eqn{I} is
#' the slope of the linear regression of log Index (`Data@Ind`) over the last
#' `yrsmth` years.
#'
#' The TAC is subject to the following conditions:
#' 1. if next TAC > 1.2 last catch, then TAC = 1.2 last catch
#' 2. if next TAC < 0.8 last catch, then TAC = 0.8 last catch
#'
#' Note that this is my interpretation of their approach and is now stochastic.
#' Currently it is generalized and is not 'tuned' to more detailed assessment
#' data which might explain why in some cases it leads to stock declines.
#'
#' @templateVar mp GB_slope
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth Number of years for evaluating slope in relative abundance
#' index
#' @param lambda A gain parameter
#' @author T. Carruthers
#' @references Geromont, H.F. and Butterworth, D.S. 2014. Complex assessment or
#' simple management procedures for efficient fisheries management: a
#' comparative study. ICES J. Mar. Sci. doi:10.1093/icesjms/fsu017
#' @family Index methods
#' @examples
#' GB_slope(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
GB_slope <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, lambda = 1) {
  dependencies = "Data@Year, Data@Cat, Data@CV_Cat, Data@Ind"
  Catrec <- Data@Cat[x, length(Data@Cat[x, ])]
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  I_hist <- Data@Ind[x, ind]
  yind <- 1:yrsmth
  # slppar <- summary(lm(I_hist ~ yind))$coefficients[2, 1:2]
  slppar <- summary(lm(log(I_hist) ~ yind))$coefficients[2, 1:2]

  if (reps >1) {
    Islp <- rnorm(reps, slppar[1], slppar[2])
  } else {
    Islp <- slppar[1]
  }

  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <-  MSEtool::trlnorm(reps, MuC, Data@CV_Cat[x,1])

  TAC <- Cc * (1 + lambda * Islp)
  TAC[TAC > (1.2 * Catrec)] <- 1.2 * Catrec
  TAC[TAC < (0.8 * Catrec)] <- 0.8 * Catrec

  TAC <- MSEtool::TACfilter(TAC)
  if (plot) GB_slope_plot(Data, ind, I_hist, MuC, TAC, Islp)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(GB_slope) <- "MP"


#' Geromont and Butterworth target CPUE and catch MP
#'
#' An MP similar to SBT2 that modifies a time-series of catch recommendations
#' and aims for target catch rate and catch level based on BMSY/B0 and MSY,
#' respectively.
#'
#' The TAC is calculated as:
#' If  \eqn{I_\textrm{recent} \geq I_0}:
#' \deqn{\textrm{TAC}= C_\textrm{ref}  \left(w + (1-w)\frac{I_\textrm{rec}-I_0}{I_\textrm{target}-I_0} \right) }
#'
#' else:
#' \deqn{\textrm{TAC}= wC_\textrm{ref}  \frac{I_\textrm{rec}}{I_0}^2}
#'
#' where \eqn{C_\textrm{ref}} is a reference catch assumed to be a proxy for MSY (`Data@Cref`),
#' *w* is a gain parameter,
#' \eqn{I_\textrm{rec}} is the average index over the last 4 years,
#' \eqn{I_\textrm{target}} is the target Index (`Data@Iref`), and
#' \eqn{I_0} is 0.2 x the average index over the past 5 years.
#'
#' In the MSE \eqn{C_\textrm{ref}} is the calculated MSY subject to observation error
#' defined in `Obs@CV_Cref`, and \eqn{I_\textrm{target}} is assumed to be the index at MSY subject
#' to observation error (`Obs@CV_Iref`). Consequently, the performance of this method in the MSE
#' is strongly determined by the specified uncertainty for these parameters.
#'
#' The TAC is subject to the following conditions:
#' 1. if next TAC > 1.2 last catch, then TAC = 1.2 last catch
#' 2. if next TAC < 0.8 last catch, then TAC = 0.8 last catch
#'
#'
#' @templateVar mp GB_target
#' @template MPtemplate
#' @template MPuses
#'
#' @param w A gain parameter
#' @author T. Carruthers
#' @references Geromont, H.F. and Butterworth, D.S. 2014. Complex assessment or
#' simple management procedures for efficient fisheries management: a
#' comparative study. ICES J. Mar. Sci. doi:10.1093/icesjms/fsu017
#' @export
#' @examples
#'  GB_target(1, MSEtool::SimulatedData, plot=TRUE)
#' @family Index methods
GB_target <- function(x, Data, reps = 100, plot=FALSE, w = 0.5) {
  dependencies = "Data@Cat, Data@Cref, Data@Iref, Data@Ind"
  Catrec <- Data@Cat[x, length(Data@Cat[x, ])]
  TACtarg <- MSEtool::trlnorm(reps, Data@Cref[x], Data@CV_Cref)
  Itarg <- MSEtool::trlnorm(reps, Data@Iref[x], Data@CV_Iref)
  Iav <- mean(Data@Ind[x, (length(Data@Ind[x, ]) - 4):length(Data@Ind[x, ])], na.rm = T)
  Irec <- mean(Data@Ind[x, (length(Data@Ind[x, ]) - 3):length(Data@Ind[x, ])], na.rm = T)
  I0 <- 0.2 * Iav
  TAC <- rep(NA, reps)
  if (Irec > I0) TAC <- TACtarg * (w + (1 - w) * ((Irec - I0)/(Itarg - I0)))
  if (Irec < I0) TAC <- TACtarg * w * (Irec/I0)^2
  TAC[TAC > (1.2 * Catrec)] <- 1.2 * Catrec
  TAC[TAC < (0.8 * Catrec)] <- 0.8 * Catrec
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) GB_target_plot(Itarg, Irec, I0, Data, Catrec, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(GB_target) <- "MP"




#### G Control #####

#' G-control MP
#'
#' A harvest control rule proposed by Carl Walters that uses trajectory in
#' inferred surplus production to make upward/downward adjustments to TAC
#' recommendations
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = \textrm{SP} \left(1-gG\right)}
#' where \eqn{\textrm{SP}} is the predicted surplus production for the next year,
#' *g* is a gain parameter, and *G* is the slope of surplus production
#' as a function of biomass over the last `yrsmth` years.
#'
#' The change in TAC is bounded by the `glim` argument, which by default does not allow
#' the TAC to decrease by more than half or increase more than twice the last annual catch.
#'
#' @templateVar mp Gcontrol
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth The number of years over which to smooth catch and biomass
#' data
#' @param gg A gain parameter
#' @param glim A constraint limiting the maximum level of change in quota
#' recommendations
#'
#' @author C. Walters and T. Carruthers
#' @references
#' Carruthers et al. 2015. Performance evaluation of simple
#' management procedures. ICES J. Mar Sci. 73, 464-482.
#' @export
#' @examples
#' Gcontrol(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
#' @family Index methods
Gcontrol <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2, glim = c(0.5, 2)) {
  dependencies = "Data@Year, Data@Cat, Data@Ind, Data@Abun"
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  years <- Data@Year[ind]
  C_dat <- log(Data@Cat[x, ind])
  C_dat[C_dat == -Inf] <- 0
  B_dat <- log(Data@Ind[x, ind]/Data@Ind[x, ind[yrsmth]] * Data@Abun[x])
  B_dat[B_dat == -Inf] <- 0
  C_hist <- exp(predict(loess(C_dat ~ ind, degree = 1)))
  B_hist <- exp(predict(loess(B_dat ~ ind, degree = 1)))
  ind <- 2:yrsmth
  ind1 <- 1:(yrsmth - 1)
  SP_hist <- B_hist[ind] - B_hist[ind1] + C_hist[ind1]
  yind <- 1:length(SP_hist)
  SP_mu <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1))
  SP_se <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1), se = T)$se.fit
  Glm <- summary(lm(SP_hist ~ B_hist[ind1]))$coefficients[2, 1:2]
  if (reps >1) {
    SP_new <- rnorm(reps, SP_mu, SP_se/2)
    G_new <- rnorm(reps, Glm[1], Glm[2]/2)
  } else {
    SP_new <- SP_mu
    G_new <- Glm[1]
  }

  TAC <- SP_new * (1 - gg * G_new)
  TAC[TAC < glim[1] * C_hist[yrsmth]] <- glim[1] * C_hist[yrsmth]
  TAC[TAC > glim[2] * C_hist[yrsmth]] <- glim[2] * C_hist[yrsmth]
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) Gcontrol_plot(years, ind1, yrsmth, SP_new, SP_hist, B_dat, B_hist, C_dat, C_hist, TAC, Data)
  # Carr<-cbind(array(rep(Data@Cat[x,],each=reps),c(reps,length(Data@Cat[x,]))),TAC)
  # Warr<-(Data@Mort[x]*exp(-Data@Mort[x]*(1:ncol(Carr))))[ncol(Carr):1]
  # Warr<-Warr/sum(Warr)
  # TAC<-apply(t(matrix(Warr,nrow=ncol(Carr),ncol=reps))*Carr,1,sum)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Gcontrol) <- "MP"



#### Index-Based Methods  ####

#' Index Confidence Interval
#'
#' The MP adjusts catch based on the value of the index in the current year relative to the
#' time series mean and standard error.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC}_y=C_{y-1} \alpha}
#' where \eqn{C_{y-1}} is the catch from the previous year, and \eqn{\alpha} is
#' defined as:
#'   \deqn{\alpha =
#'              \left\{\begin{array}{ll}
#'              d & \textrm{if } I < \textrm{CI}_L \\
#'              u & \textrm{if } I > \textrm{CI}_H \\
#'              1 & \textrm{if } \textrm{CI}_L \leq I \leq \textrm{CI}_H \\
#'              \end{array}\right.
#'            }
#'  where \eqn{I} is the index in the most recent year, \eqn{d} is 0.75 for `ICI` and `ICI2`,
#'  \eqn{u} is 1.05 and 1.25 for`ICI` and `ICI2` respectively, and \eqn{\textrm{CI}_L}
#'  and \eqn{\textrm{CI}_L} are the lower and upper bound of the confidence
#'  interval of mean historical index. The confidence interval is
#'  calculated using Z-scores described in the Functions section below.
#'
#'
#' @templateVar mp ICI
#' @template MPtemplate
#' @template MPuses
#'
#' @describeIn ICI The catch is reduced by 0.75 if the Z-score of the current year's index is less than -0.44.
#' The catch is increased by 1.05 if the Z-score of the current year's index is greater than 1.96. Otherwise, the catch is held constant.
#'
#' @author Coded by Q. Huynh. Developed by Jardim et al. (2015)
#' @references Ernesto Jardim, Manuela Azevedo, Nuno M. Brites, Harvest control rules for
#' data limited stocks using length-based reference points and survey biomass indices,
#' Fisheries Research, Volume 171, November 2015, Pages 12-19, ISSN 0165-7836,
#' https://doi.org/10.1016/j.fishres.2014.11.013
#' @family Index methods
#' @export
#' @examples
#' ICI(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
#'
ICI <- function(x, Data, reps=100, plot=FALSE) {
  dependencies = "Data@Ind, Data@CV_Ind, Data@Cat, Data@CV_Cat"

  Index <- Data@Ind[x, ]
  Years <- Data@Year[!is.na(Index)]
  Index <- Index[!is.na(Index)]

  nI <- length(Index)
  if (reps ==1) {
    Ind.samp <-Index
  }else {
    Ind.samp <- MSEtool::trlnorm(reps * nI, Index, Data@CV_Ind[x,1])
  }
  Ind.samp <- matrix(Ind.samp, ncol = reps)

  muI <- apply(Ind.samp, 2, mean, na.rm = TRUE)
  sigmaI <- apply(Ind.samp, 2, sd, na.rm = TRUE)

  Ind <- Ind.samp[nI, ]
  z.low <- -0.44 #qnorm(0.33)
  z.upp <- 1.96 #qnorm(0.974)

  ci.low <- muI + z.low * sigmaI / sqrt(nI)
  ci.high <- muI + z.upp * sigmaI / sqrt(nI)

  alpha <- rep(NA, reps)
  alpha[Ind < ci.low] <- 0.75
  alpha[Ind > ci.high] <- 1.05
  alpha[Ind >= ci.low & Ind <= ci.high] <- 1

  Cat <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, Cat, Data@CV_Cat[x,1])

  TAC <- MSEtool::TACfilter(alpha * Cc)
  if (plot) ICI_plot(Years, Index, ci.low, ci.high, TAC, Cat, Data)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(ICI) <- "MP"


#' @templateVar mp ICI2
#' @template MPuses
#' @describeIn ICI This method is less precautionary of the two ICI MPs by allowing for a
#' larger increase in TAC and a lower threshold of the index to decrease the TAC. The catch is reduced by 0.75
#' if the Z-score of the current year's index is less than -1.96. The catch is increased by 1.25
#' if the Z-score of the current year's index is greater than 1.96. Otherwise, the catch is held constant.
#' @export
#' @examples
#' ICI2(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
ICI2 <- function(x, Data, reps=100, plot=FALSE) {
  dependencies = "Data@Ind, Data@CV_Ind, Data@Cat, Data@CV_Cat"

  Index <- Data@Ind[x, ]
  Years <- Data@Year[!is.na(Index)]
  Index <- Index[!is.na(Index)]
  nI <- length(Index)
  Ind.samp <- MSEtool::trlnorm(reps * nI, Index, Data@CV_Ind[x,1])
  Ind.samp <- matrix(Ind.samp, ncol = reps)

  muI <- apply(Ind.samp, 2, mean, na.rm = TRUE)
  sigmaI <- apply(Ind.samp, 2, sd, na.rm = TRUE)

  Ind <- Ind.samp[nI, ]
  z.low <- -1.96 #qnorm(0.025)
  z.upp <- 1.96 #qnorm(0.975)

  ci.low <- muI + z.low * sigmaI / sqrt(nI)
  ci.high <- muI + z.upp * sigmaI / sqrt(nI)

  alpha <- rep(NA, reps)
  alpha[Ind < ci.low] <- 0.75
  alpha[Ind > ci.high] <- 1.25
  alpha[Ind >= ci.low & Ind <= ci.high] <- 1

  Cat <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, Cat, Data@CV_Cat[x,1])

  TAC <- alpha * Cc
  TAC <- MSEtool::TACfilter(alpha * Cc)
  if (plot) ICI_plot(Years, Index, ci.low, ci.high, TAC, Cat, Data)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(ICI2) <- "MP"


#' Mean Index Ratio
#'
#' The TAC is adjusted by the ratio alpha, where the numerator
#' being the mean index in the most recent two years of the time series and the denominator
#' being the mean index in the three years prior to those in the numerator. This MP is the
#' stochastic version of Method 3.2 used by ICES for Data-Limited Stocks (ICES 2012).
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC}_y = \alpha C_{y-1}}
#' where \eqn{C_{y-1}} is the catch from the previous year, and \eqn{\alpha} is
#' the ratio of the mean index in the most recent two years of the time series and the
#' mean index in 3-5 years before current time (reference years are specified as
#' `yrs` argument.
#'
#' @templateVar mp Iratio
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrs Vector of length 2 specifying the reference years
#' @author Coded by Q. Huynh. Developed by Jardim et al. (2015)
#' @references Ernesto Jardim, Manuela Azevedo, Nuno M. Brites, Harvest control rules for
#' data limited stocks using length-based reference points and survey biomass indices,
#' Fisheries Research, Volume 171, November 2015, Pages 12-19, ISSN 0165-7836,
#' https://doi.org/10.1016/j.fishres.2014.11.013
#'
#' ICES. 2012. ICES Implementation of Advice for Data-limited Stocks in 2012 in its 2012
#' Advice. ICES CM 2012/ACOM 68. 42 pp.
#' @family Index methods
#' @examples
#' Iratio(1, MSEtool::Atlantic_mackerel, plot=TRUE)
#' @export
Iratio <- function(x, Data, reps=100, plot=FALSE, yrs = c(2, 5)) {
  dependencies = "Data@Ind, Data@CV_Ind, Data@Cat, Data@CV_Cat"

  ind.num <- (length(Data@Year) - yrs[1]+1):length(Data@Year)
  ind.den <- (length(Data@Year) - yrs[2]+1):(length(Data@Year) - yrs[1])

  if (reps ==1 ) {
    I.num <- Data@Ind[x, ind.num]
    I.den <-  Data@Ind[x, ind.den]
  } else {
    I.num <- MSEtool::trlnorm(reps * length(ind.num), Data@Ind[x, ind.num], Data@CV_Ind[x,ind.num])
    I.den <- MSEtool::trlnorm(reps * length(ind.den), Data@Ind[x, ind.den], Data@CV_Ind[x,ind.den])
  }

  I.num <- matrix(I.num, ncol = reps)
  I.den <- matrix(I.den, ncol = reps)

  alpha <- apply(I.num, 2, mean, na.rm = TRUE)/apply(I.den, 2, mean, na.rm = TRUE)
  Cat <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, Cat, Data@CV_Cat[x,1])
  TAC <- alpha * Cc
  TAC <- MSEtool::TACfilter(TAC)

  if(plot) Iratio_plot(Data, I.num, ind.num, I.den, ind.den, alpha, TAC, Cat)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Iratio) <- "MP"







#' Index Slope Internal Function
#'
#' @param x Iteration number
#' @param Data An object of class `Data`
#' @param reps Number of samples of the TAC
#' @param yrsmth Years over which to calculate index
#' @param lambda  A gain parameter controlling the speed in update in TAC.
#' @param xx Parameter controlling the fraction of mean catch to start using in
#' first year
#'
#' @export
#'
#' @keywords internal
Islope_ <- function(x, Data, reps = 100, yrsmth = 5, lambda = 0.4,xx = 0.2) {

  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  Years <- Data@Year[ind]
  ylast <- (Data@LHYear[1] - Data@Year[1]) + 1  #last historical year
  C_dat <- Data@Cat[x, ind]
  if (is.na(Data@MPrec[x]) || length(Data@Year) == ylast + 1) {
    TACstar <- (1 - xx) * MSEtool::trlnorm(reps, mean(C_dat, na.rm=TRUE), Data@CV_Cat/(yrsmth^0.5))
  } else {
    TACstar <- rep(Data@MPrec[x], reps)
  }

  I_hist <- Data@Ind[x, ind]
  yind <- 1:yrsmth
  # slppar <- summary(lm(I_hist ~ yind))$coefficients[2, 1:2]
  slppar <- summary(lm(log(I_hist) ~ yind))$coefficients[2, 1:2]
  if (reps >1) {
    Islp <- rnorm(reps, slppar[1], slppar[2])
  } else {
    Islp <- slppar[1]
  }
  TAC <- TACstar * (1 + lambda * Islp)

  list(TAC=TAC, TACstar=TACstar, I_hist=I_hist, Islp=Islp, C_dat=C_dat, Data=Data, Years=Years)

}


#' Index Slope Tracking MP
#'
#'
#' A management procedure that incrementally adjusts the TAC to maintain a
#' constant CPUE or relative abundance index.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = \textrm{TAC}^* \left(1+\lambda I \right)}
#' where \eqn{\textrm{TAC}^*} is \eqn{1-xx} multiplied by the mean catch from the past `yrsmth` years for the
#' first year and catch from the previous year in projection years,
#' \eqn{\lambda} is a gain parameter, and \eqn{I} is the slope of log index over the past `yrsmth` years.
#'
#' @templateVar mp Islope1
#' @template MPtemplate
#' @template MPuses
#' @param yrsmth Years over which to calculate index
#' @param lambda A gain parameter controlling the speed in update in TAC.
#' @param xx Parameter controlling the fraction of mean catch to start using in
#' first year
#'
#' @author T. Carruthers
#' @references Carruthers et al. 2015. Performance evaluation of simple
#' management procedures. ICES J. Mar Sci. 73, 464-482.
#'
#' Geromont, H.F., Butterworth, D.S. 2014. Generic management procedures for
#' data-poor fisheries; forecasting with few data. ICES J. Mar. Sci.
#' doi:10.1093/icesjms/fst232
#'
#' @describeIn Islope1 The least biologically precautionary of the Islope methods
#' @export
#' @examples
#' Islope1(1, MSEtool::SimulatedData, plot=TRUE)
#' @family Index methods
Islope1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, lambda = 0.4,xx = 0.2) {

  runIslope <- Islope_(x, Data, reps, yrsmth, lambda, xx)
  TAC <- MSEtool::TACfilter(runIslope$TAC)
  runIslope$TAC <- TAC

  if(plot) Islope_plot(runIslope, Data)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec

}
class(Islope1) <- "MP"




#' @describeIn Islope1 More biologically precautionary. Reference TAC is 0.7 average catch
#' @export
#' @examples
#' Islope2(1, MSEtool::SimulatedData, plot=TRUE)
Islope2 <- Islope1
formals(Islope2)$lambda <- 0.4
formals(Islope2)$xx <- 0.3
class(Islope2) <- "MP"


#' @describeIn Islope1 More biologically precautionary. Reference TAC is 0.6 average catch
#' @export
#' @examples
#' Islope3(1, MSEtool::SimulatedData, plot=TRUE)
Islope3 <- Islope1
formals(Islope3)$lambda <- 0.4
formals(Islope3)$xx <- 0.4
class(Islope3) <- "MP"


#' @describeIn Islope1 The most biologically precautionary of the Islope methods.
#' Reference TAC is 0.6 average catch and gain parameter is 0.2
#' @export
#' @examples
#' Islope4(1, MSEtool::SimulatedData, plot=TRUE)
Islope4 <- Islope1
formals(Islope4)$lambda <- 0.2
formals(Islope4)$xx <- 0.4
class(Islope4) <- "MP"






#' Index Target Internal Function
#'
#' @param x Iteration number
#' @param Data An object of class `Data`
#' @param reps Number of replicates
#' @param plot Logical. Show the plot?
#' @param yrsmth The number of historical years over which to average the index
#' @param mc The maximum fractional change in the TAC among years.
#'
#' @return A named list with TAC recommendations of length `reps`
#' @export
#'
#' @keywords internal
IT_ <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, mc = 0.05) {

  dependencies = "Data@Ind, Data@MPrec, Data@CV_Ind, Data@Iref"
  ind <- max(1, (length(Data@Year) - yrsmth + 1)):length(Data@Year)
  if(is.na(Data@Iref[x])) return(list(TAC=rep(as.numeric(NA), reps)))
  deltaI <- mean(Data@Ind[x, ind], na.rm=TRUE)/Data@Iref[x]
  if (deltaI < (1 - mc)) deltaI <- 1 - mc
  if (deltaI > (1 + mc)) deltaI <- 1 + mc
  TAC <- Data@MPrec[x] * deltaI * MSEtool::trlnorm(reps, 1, Data@CV_Ind[x,1])
  TAC <- MSEtool::TACfilter(TAC)
  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(op)
    par(mfrow=c(1,2), oma=c(1,1,1,1), mar=c(5,4,1,4))
    ylim <- range(c(Data@Ind[x,ind], Data@Iref[x]))
    plot(Data@Year[ind], Data@Ind[x,ind], xlab="Year",
         ylab= paste0("Index (previous ", yrsmth, "years)"), bty="l", type="l",
         lwd=2, ylim=ylim)
    lines(Data@Year[ind], rep(mean(Data@Ind[x, ind], na.rm=TRUE), length(ind)), lty=2)
    text(quantile(Data@Year[ind],0.15), mean(Data@Ind[x, ind], na.rm=TRUE), "Mean Index", pos=3)
    lines(Data@Year[ind], rep(mean(Data@Iref[x], na.rm=TRUE), length(ind)), lty=3)
    text(quantile(Data@Year[ind],0.15), Data@Iref[x], "Reference Index", pos=3)

    boxplot(TAC, ylab=paste0("TAC (", Data@Units, ")"))
    points(1, Data@MPrec[x], cex=2, col="blue", pch=16)
    text(1, Data@MPrec[x], cex=1, col="blue", "Last Catch", pos=1)
  }
  list(TAC=TAC)

}

#' Iterative Index Target MP
#'
#' An index target MP where the TAC is modified according to current index
#' levels (mean index over last 5 years) relative to a target level.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC}_y = C_{y-1} I_\delta}
#' where \eqn{C_{y-1}} is the catch from the previous year and \eqn{I_\delta} is
#' the ratio of the mean index over the past `yrsmth` years to a reference index level.
#' The maximum allowable change in TAC is determined by `mc`: e.g `mc=0.05`
#' means that the maximum change in TAC from the previous
#' catch is 5%.
#'
#' The reference index level (`Data@Iref`) is assumed to be a proxy for MSY. In the
#' MSE `Iref` is the index at MSY subject to observation error (`Obs@Irefbiascv`).
#' Consequently the performance of these methods in MSE is strongly determined by the
#' uncertainty the in reference index.
#'
#' @templateVar mp IT5
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth The number of historical years over which to average the index
#' @param mc The maximum fractional change in the TAC among years.
#'
#' @author T. Carruthers
#' @describeIn IT5 Maximum annual changes in TAC are 5 per cent.
#' @export
#' @examples
#' IT5(1, MSEtool::SimulatedData, plot=TRUE)
IT5 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, mc = 0.05) {
  runIT <- IT_(x, Data, reps, plot, yrsmth, mc)
  Rec <- new("Rec")
  Rec@TAC <- runIT$TAC
  Rec
}
class(IT5) <- "MP"

#' @templateVar mp IT10
#' @template MPuses
#'
#' @describeIn IT5 Maximum annual changes are 10 per cent.
#' @export
#' @examples
#' IT10(1, MSEtool::SimulatedData, plot=TRUE)
IT10 <- IT5
formals(IT10)$mc <- 0.1
class(IT10) <- "MP"



#' Incremental Index Target MP Internal Function
#'
#' @param x Iteration number
#' @param Data An object of class `Data`
#' @param reps Number of replicates
#' @param plot Logical. Show the plot?
#' @param yrsmth Years over which the average index is calculated.
#' @param xx Parameter controlling the fraction of mean catch to start using in
#' first year
#' @param Imulti Parameter controlling how much larger target CPUE / index is
#' compared with recent levels.
#'
#' @return A named list with TAC of length `reps`
#' @export
#'
#' @keywords internal
Itarget_ <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, Imulti = 1.5) {
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)  # recent 5 years
  ylast <- (Data@LHYear[1] - Data@Year[1]) + 1  #last historical year
  ind2 <- ((ylast - (yrsmth - 1)):ylast)  # historical 5 pre-projection years
  ind3 <- ((ylast - (yrsmth * 2 - 1)):ylast)  # historical 10 pre-projection years
  C_dat <- Data@Cat[x, ind2]
  TACstar <- (1 - xx) * MSEtool::trlnorm(reps, mean(C_dat, na.rm=TRUE), Data@CV_Cat[x,1]/(yrsmth^0.5))
  Irecent <- mean(Data@Ind[x, ind], na.rm=TRUE)
  Iave <- mean(Data@Ind[x, ind3], na.rm=TRUE)
  Itarget <- Iave * Imulti
  I0 <- 0.8 * Iave
  if (Irecent > I0) {
    TAC <- 0.5 * TACstar * (1 + ((Irecent - I0)/(Itarget - I0)))
  } else {
    TAC <- 0.5 * TACstar * (Irecent/I0)^2
  }
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(op)
    par(mfrow=c(1,2))

    ylim <- range(c(Data@Ind[x, ], Itarget, I0))
    plot(Data@Year, Data@Ind[x, ], type="l", lwd=2, bty="l",
         xlab="Year", ylab="Index", ylim=ylim)
    points(max(Data@Year), mean(Data@Ind[x, ind], na.rm=TRUE), cex=2, pch=16,col='blue')
    text(max(Data@Year), mean(Data@Ind[x, ind], na.rm=TRUE), cex=1, 'Irecent', pos=3, col='blue', xpd=NA)

    lines(Data@Year[ind3], rep(mean(Data@Ind[x, ind3], na.rm=TRUE), length(ind3)), lty=2, col="orange")
    text(mean(Data@Year[ind3]), mean(Data@Ind[x, ind3], na.rm=TRUE), "Iave", col="orange", pos=1)

    points(max(Data@Year), Itarget, cex=2, pch=16,col='green')
    text(max(Data@Year), Itarget, cex=1, 'Itarget', pos=3, col='green', xpd=NA)

    points(max(Data@Year), I0, cex=2, pch=16,col='red')
    text(max(Data@Year), I0, cex=1, 'I0', pos=3, col='red', xpd=NA)

    ylim <- range(c(C_dat, TACstar, TAC))
    Years <- Data@Year[ind2]
    if (max(Years) != max(Data@Year)) {
      Years <- c(Years,  (max(Data@Year[ind2])+1):max(Data@Year))
    }

    yrs <- length(Years)-length(C_dat)
    Cdat <- c(C_dat, rep(NA, yrs))
    plot(c(Years, max(Years)+1), c(Cdat, NA), type="l", lwd=2, bty="l",
         xlab="Year", ylab=paste0("Catch (",  Data@Units, ")"), ylim=ylim)
    abline(v=max(Data@Year[ind2]), col="gray", lty=3)

    # lines(Years[(length(Years)-yrs):length(Years)], Data@Cat[x, (length(Data@Cat[x,])-yrs):length(Data@Cat[x,])])

    points(max(Data@Year[ind2]), mean(TACstar, na.rm=TRUE), cex=2, col="orange", pch=16)
    text(max(Data@Year[ind2]), mean(TACstar, na.rm=TRUE), "mean TAC*", pos=2, xpd=NA, col="orange")

    boxplot(TAC, at=max(Years)+1, add=TRUE, col="gray", axes=FALSE)

  }

  list(TAC=TAC )
}




#' Incremental Index Target MP
#'
#'
#' A management procedure that incrementally adjusts the TAC (starting from
#' reference level that is a fraction of mean recent catches)
#' to reach a target CPUE / relative abundance index
#'
#' Four index/CPUE target MPs proposed by Geromont and Butterworth 2014.
#' Tested by Carruthers et al. 2015.
#'
#' The TAC is calculated as:
#' If  \eqn{I_\textrm{recent} \geq I_0}:
#' \deqn{\textrm{TAC}= 0.5 \textrm{TAC}^* \left[1+\left(\frac{I_\textrm{recent} - I_0}{I_\textrm{target} - I_0}\right)\right]}
#'
#' else:
#' \deqn{\textrm{TAC}= 0.5 \textrm{TAC}^* \left[\frac{I_\textrm{recent}}{I_0}^2\right]}
#'
#' where \eqn{I_0} is \eqn{0.8 I_{\textrm{ave}}} (the average index over the 2 x `yrsmth` years prior to the projection period),
#' \eqn{I_\textrm{recent}} is the average index over the past `yrsmth` years, and
#' \eqn{I_\textrm{target}} is `Imulti` times \eqn{I_{\textrm{ave}}},
#' and \eqn{\textrm{TAC}^*} is:
#' \deqn{(1-x)C}
#' where \eqn{x} is argument `xx` and C is the average catch over the last 5 years of the historical period.
#'
#' @templateVar mp Itarget1
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth Years over which the average index is calculated.
#' @param xx Parameter controlling the fraction of mean catch to start using in
#' first year
#' @param Imulti Parameter controlling how much larger target CPUE / index is
#' compared with recent levels.

#' @author T. Carruthers
#' @references Carruthers et al. 2015. Performance evaluation of simple
#' management procedures. ICES J. Mar Sci. 73, 464-482.
#'
#' Geromont, H.F., Butterworth, D.S. 2014. Generic management procedures for
#' data-poor fisheries; forecasting with few data. ICES J. Mar. Sci. 72, 251-261.
#' doi:10.1093/icesjms/fst232
#' @describeIn Itarget1 The less precautionary TAC-based MP
#' @family Index methods
#' @export
#' @examples
#' Itarget1(1, MSEtool::Atlantic_mackerel, plot=TRUE)
Itarget1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, Imulti = 1.5) {
  runItarget <- Itarget_(x, Data, reps, plot, yrsmth, xx, Imulti)
  Rec <- new("Rec")
  Rec@TAC <- runItarget$TAC
  Rec
}
class(Itarget1) <- "MP"

#' @describeIn Itarget1 Increasing biologically precautionary TAC-based MP
#' @export
#' @examples
#' Itarget2(1, MSEtool::Atlantic_mackerel, plot=TRUE)
Itarget2 <- Itarget1
formals(Itarget2)$Imulti <- 2
class(Itarget2) <- "MP"

#' @describeIn Itarget1 Increasing biologically precautionary TAC-based MP
#' @export
#' @examples
#' Itarget3(1, MSEtool::Atlantic_mackerel, plot=TRUE)
Itarget3 <- Itarget1
formals(Itarget3)$Imulti <- 2.5
class(Itarget3) <- "MP"

#' @describeIn Itarget1 The most biologically precautionary TAC-based MP
#' @export
#' @examples
#' Itarget4(1, MSEtool::Atlantic_mackerel, plot=TRUE)
Itarget4 <- Itarget1
formals(Itarget4)$xx <- 0.3
formals(Itarget4)$Imulti <- 2.5
class(Itarget4) <- "MP"



#' Index Target based on natural mortality rate
#'
#' An index target MP where the TAC is modified according to current index
#' levels (mean index over last number of years determined by natural mortality
#' (*M*)) relative to a target level.
#'
#' The TAC is caluclated as:
#' \deqn{\textrm{TAC}_y = \textrm{TAC}_{y-1} \delta I}
#' where \eqn{\delta I} is the ratio of the  mean index over \eqn{4\frac{1}{M}^{1/4}} years
#' to the reference index (`Data@Iref`).
#'
#' The maximum fractional change in TAC is determined by \eqn{mc}, defined as
#' \eqn{mc = \textrm{max}\left(\frac{5 + 25M}{100}, 0.2\right)}
#'
#' @templateVar mp ITM
#' @template MPtemplate
#' @template MPuses
#'
#' @author T. Carruthers
#' @export
#' @examples
#' ITM(1, Data=MSEtool::SimulatedData, plot=TRUE)
ITM <- function(x, Data, reps = 100, plot=FALSE) {

  dependencies = "Data@Ind, Data@Cat, Data@CV_Ind, Data@Iref, Data@Mort"
  mc <- (5 + Data@Mort[x] * 25)/100
  if (mc > 0.2)   mc <- 0.2
  yrsmth <- floor(4 * (1/Data@Mort[x])^(1/4))
  ind <- max(1, (length(Data@Year) - yrsmth + 1)):length(Data@Year)

  deltaI <- mean(Data@Ind[x, ind], na.rm=TRUE)/Data@Iref[x]
  if (deltaI < (1 - mc)) deltaI <- 1 - mc
  if (deltaI > (1 + mc)) deltaI <- 1 + mc

  TAC <- Data@MPrec[x] * deltaI * MSEtool::trlnorm(reps, 1, Data@CV_Ind[x,1])
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(op)
    par(mfrow=c(1,2))
    ylim <- range(c(Data@Ind[x,ind], Data@Iref[x]))
    plot(Data@Year[ind], Data@Ind[x,ind], type="l", xlab="Year", ylab="Index",
         lwd=2,  bty="l", ylim=ylim)
    lines(Data@Year[ind], rep(Data@Iref[x], length(ind)), lty=2, col="darkgray")
    text(quantile(Data@Year[ind],0.15), Data@Iref[x], "Iref", pos=3)

    ylim <- range(c(TAC,Data@MPrec[x] ))
    boxplot(TAC, col="gray", ylab=paste0("TAC (", Data@Units, ")"), ylim=ylim)
    points(1, Data@MPrec[x], cex=2, pch=16, col="green")
    text(1, Data@MPrec[x], "Last TAC", col="green", pos=3)
  }

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(ITM) <- "MP"



### Length Target Methods #####
Ltarget_ <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, xL = 1.05,
                     Ltarget=NULL, L0=NULL) {
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)  # recent 5 years
  ylast <- (Data@LHYear[1] - Data@Year[1]) + 1  #last historical year
  ind2 <- ((ylast - (yrsmth - 1)):ylast)  # historical 5 pre-projection years
  ind3 <- ((ylast - (yrsmth * 2 - 1)):ylast)  # historical 10 pre-projection years
  C_dat <- Data@Cat[x, ind2]
  TACstar <- (1 - xx) * MSEtool::trlnorm(reps, mean(C_dat, na.rm=TRUE), Data@CV_Cat/(yrsmth^0.5))
  Lrecent <- mean(Data@ML[x,ind], na.rm=TRUE)
  Lave <- mean(Data@ML[x,ind3], na.rm=TRUE)
  if (is.null(L0)) L0 <- 0.9 * Lave
  if (is.null(Ltarget)) Ltarget <- xL * Lave
  if (Lrecent >= L0) {
    TAC <- 0.5 * TACstar * (1 + ((Lrecent - L0)/(Ltarget - L0)))
  } else {
    TAC <- 0.5 * TACstar * (Lrecent/L0)^2
  }

  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(op)
    par(mfrow=c(1,2))
    ylim <- range(c(Data@ML[x,], Lrecent, Lave, Ltarget, L0))
    plot(Data@Year, Data@ML[x,], type="l", xlab="Year", ylab="Mean Length",
         lwd=2,  bty="l", ylim=ylim)
    lines(Data@Year[ind], rep(Lrecent, length(ind)), lty=2, col="green")
    text(quantile(Data@Year[ind],0.15),Lrecent, "Lrecent", pos=3, col='green')

    lines(Data@Year[ind3], rep(Lave, length(ind3)), lty=2, col="orange")
    text(quantile(Data@Year[ind3],0.15),Lave, "Lave", pos=3, col='orange')

    points(max(Data@Year[ind3]), Ltarget, lty=2, col="blue", pch=16)
    text(max(Data@Year[ind3]), Ltarget, "Ltarget", pos=3, col='blue', xpd=NA)

    points(max(Data@Year[ind3]), L0, lty=2, col="red", pch=16)
    text(max(Data@Year[ind3]), L0, "L0", pos=3, col='red', xpd=NA)


    ylim <- range(c(Data@Cat[x,], TAC, TACstar ))

    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x,],NA), type="l", xlab="Year",
    ylab=paste0("Catch (", Data@Units, ")"),
         lwd=2,  bty="l", ylim=ylim)

    boxplot(TAC, col="blue", axes=FALSE, at=max(Data@Year)+1, add=TRUE)
    text(max(Data@Year)+1, quantile(TAC, 0.95, na.rm=TRUE), "TAC", col="blue")
    points(max(Data@Year[ind3]), mean(TACstar, na.rm=TRUE), cex=2, pch=16, col="green")
    text(max(Data@Year[ind3]), mean(TACstar, na.rm=TRUE), "TAC*", col="green", pos=1)
  }

  list(TAC=TAC)

}


#' Length Target TAC MP
#'
#' A management procedure that incrementally adjusts the TAC to reach
#' a target mean length in catches.
#'
#' Four target length MPs proposed by Geromont and Butterworth 2014.
#' Tested by Carruthers et al. 2015.
#'
#' The TAC is calculated as:
#'
#' If \eqn{L_\textrm{recent} \geq L_0}:
#' \deqn{\textrm{TAC} = 0.5 \textrm{TAC}^* \left[1+\left(\frac{L_\textrm{recent}-L_0}{L_\textrm{target}-L_0}\right)\right]  }
#'
#' else:
#' \deqn{\textrm{TAC} = 0.5 \textrm{TAC}^* \left[\frac{L_\textrm{recent}}{L_0}^2\right] }
#'
#' where \eqn{\textrm{TAC}^*} is (1 - `xx`) mean catches from the last `yrsmth` historical years (pre-projection),
#' \eqn{L_\textrm{recent}} is mean length in last `yrmsth` years, \eqn{L_0} is (except for `L95target`) 0.9 average catch in the last
#' 2 x `yrsmth` historical (pre-projection years) (\eqn{L_\textrm{ave}}), and \eqn{L_\textrm{target}} is
#' (except for `L95target`) `xL` \eqn{L_\textrm{ave}}.
#'
#' @templateVar mp Ltarget1
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth Years over which to calculate mean length.
#' @param xx Parameter controlling the fraction of mean catch to start using in
#' first year
#' @param xL Parameter controlling the magnitude of the target mean length of
#' catches relative to average length in catches.
#' @author T. Carruthers
#' @references Carruthers et al. 2015. Performance evaluation of simple
#' management procedures. ICES J. Mar Sci. 73, 464-482.
#'
#' Geromont, H.F., Butterworth, D.S. 2014. Generic management procedures for
#' data-poor fisheries; forecasting with few data. ICES J. Mar. Sci.
#' doi:10.1093/icesjms/fst232
#' @describeIn Ltarget1 The least biologically precautionary TAC-based MP.
#' @family Length target MPs
#' @export
#' @examples
#' Ltarget1(1, Data=MSEtool::SimulatedData, plot=TRUE)
Ltarget1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, xL = 1.05) {
  runLtarget <- Ltarget_(x, Data, reps, plot, yrsmth, xx, xL)
  Rec <- new("Rec")
  Rec@TAC <- runLtarget$TAC
  Rec
}
class(Ltarget1) <- "MP"

#' @templateVar mp Ltarget2
#' @template MPuses
#' @describeIn Ltarget1 Increasingly biologically precautionary (`xL` = 1.1).
#' @export
#' @examples
#' Ltarget2(1, Data=MSEtool::SimulatedData, plot=TRUE)
Ltarget2 <- Ltarget1
formals(Ltarget2)$xL <- 1.1
class(Ltarget2) <- "MP"

#' @templateVar mp Ltarget3
#' @template MPuses
#' @describeIn Ltarget1 Increasingly biologically precautionary (`xL` = 1.1).
#' @export
#' @examples
#' Ltarget3(1, Data=MSEtool::SimulatedData, plot=TRUE)
Ltarget3 <- Ltarget1
formals(Ltarget3)$xL <- 1.15
class(Ltarget3) <- "MP"

#' @templateVar mp Ltarget4
#' @template MPuses
#' @describeIn Ltarget1 The most biologically precautionary TAC-based MP (`xL` = 1.1, `xx`=0.2).
#' @export
#' @examples
#' Ltarget4(1, Data=MSEtool::SimulatedData, plot=TRUE)
Ltarget4 <- Ltarget1
formals(Ltarget4)$xx <- 0.2
formals(Ltarget4)$xL <- 1.15
class(Ltarget4) <- "MP"


#' @templateVar mp L95target
#' @template MPuses
#' @describeIn Ltarget1 Same as Ltarget1 but here the target and limit
#' mean lengths are based on the length at maturity distribution rather
#' than an arbitrary multiplicative of the mean length
#' @export
#' @examples
#' L95target(1, Data=MSEtool::SimulatedData, plot=TRUE)
L95target <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, xL=1.05) {
  runLtarget <- Ltarget_(x, Data, reps, plot, yrsmth, xx, xL, Ltarget = Data@L95[x],
                         L0=0.9 * Data@L50[x])
  Rec <- new("Rec")
  Rec@TAC <- runLtarget$TAC
  Rec
}
class(L95target) <- "MP"





## Mean Length-Based Indicator ####

#' Mean length-based indicator MP of Jardim et al. 2015
#'
#' The TAC is calculated as the most recent catch, modified by the ratio alpha,
#' where the numerator is the mean length of the catch (of lengths larger than Lc) and
#' the denominator is the mean length expected at MSY. Here, Lc is the length at
#' full selection (LFS).
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC}_y = C_{y-1} \frac{L}{L_\textrm{ref}}}
#' where \eqn{C_{y-1}} is the catch from the previous year, \eqn{L} is the mean
#' length of the catch over the last `yrsmth` years (of lengths larger than Lc) and \eqn{L_\textrm{ref}} is the
#' mean length expected at MSY. Here, Lc is the length at full selection (LFS).
#'
#'
#' @templateVar mp Lratio_BHI
#' @template MPtemplate
#' @template MPuses

#'
#' @param yrsmth The most recent years of data to smooth the calculation of the mean length
#' @author Coded by Q. Huynh. Developed by Jardim et al. (2015)
#' @references Ernesto Jardim, Manuela Azevedo, Nuno M. Brites, Harvest control rules for
#' data limited stocks using length-based reference points and survey biomass indices,
#' Fisheries Research, Volume 171, November 2015, Pages 12-19, ISSN 0165-7836,
#' https://doi.org/10.1016/j.fishres.2014.11.013
#'
#' @describeIn Lratio_BHI Assumes M/K = 1.5 and FMSY/M = 1. Natural mortality M and von Bertalanffy
#' K are not used in this MP (see Appendix A of Jardim et al. 2015).
#' @family Length target MPs
#' @export
#' @examples
#' Lratio_BHI(1, Data=MSEtool::SimulatedData, plot=TRUE)
#'
Lratio_BHI <- function(x, Data, reps=100, plot=FALSE, yrsmth = 3) {
  dependencies = "Data@vb_Linf, Data@CV_vbLinf, Data@Cat, Data@CV_Cat, Data@CAL, Data@CAL_bins,
  Data@LFS, Data@CV_LFS"
  Linfc <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Lc <- MSEtool::trlnorm(reps, Data@LFS[x], Data@CV_LFS[x])
  Lref <- 0.75 * Lc + 0.25 * Linfc

  Cat <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, Cat, Data@CV_Cat[x,1])

  LYear <- dim(Data@CAL)[2]
  nlbin <- ncol(Data@CAL[x,,])
  mlbin <- 0.5 * (Data@CAL_bins[1:nlbin] + Data@CAL_bins[2:(nlbin + 1)])

  ind.year <- (LYear - yrsmth + 1):LYear
  CAL <- colSums(Data@CAL[x, ind.year, ])

  LSQ <- rep(NA, reps)

  for(i in 1:reps) {
    lensamp <- sample(mlbin, 0.5*sum(CAL), replace = T, prob = CAL)
    LSQ[i] <- mean(lensamp, na.rm=TRUE)
  }

  TAC <- MSEtool::TACfilter((LSQ/Lref) * Cc)

  if (plot) Lratio_BHI_plot(mlbin, CAL, LSQ, Lref, Data, x, TAC, Cc, yrsmth)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Lratio_BHI) <- "MP"



#' @templateVar mp Lratio_BHI2
#' @template MPuses
#'
#' @describeIn Lratio_BHI More general version that calculates the reference mean length as a
#' function of M, K, and presumed FMSY/M.
#' @export
#' @examples
#' Lratio_BHI2(1, Data=MSEtool::SimulatedData, plot=TRUE)
#'
Lratio_BHI2 <- function(x, Data, reps=100, plot=FALSE, yrsmth = 3) {

  dependencies = "Data@vb_Linf, Data@CV_vbLinf, Data@Cat, Data@CV_Cat, Data@Mort, Data@CV_Mort,
  Data@vb_K, Data@CV_vbK, Data@FMSY_M, Data@CV_FMSY_M, Data@CAL, Data@CAL_bins,
  Data@LFS, Data@CV_LFS"
  Linfc <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Lc <- MSEtool::trlnorm(reps, Data@LFS[x], Data@CV_LFS[x])

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  gamma <- MSEtool::trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x])
  theta <- Kvec/Mvec

  Lref <- (theta * Linfc + Lc * (gamma + 1)) / (gamma + theta + 1)

  Cat <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, Cat, Data@CV_Cat[x,1])

  LYear <- dim(Data@CAL)[2]
  nlbin <- ncol(Data@CAL[x,,])
  mlbin <- 0.5 * (Data@CAL_bins[1:nlbin] + Data@CAL_bins[2:(nlbin + 1)])

  ind.year <- (LYear - yrsmth + 1):LYear
  CAL <- colSums(Data@CAL[x, ind.year, ])

  LSQ <- rep(NA, reps)
  for(i in 1:reps) {
    lensamp <- sample(mlbin, 0.5*sum(CAL), replace = T, prob = CAL)
    LSQ[i] <- mean(lensamp, na.rm=TRUE)
  }

  TAC <- MSEtool::TACfilter((LSQ/Lref) * Cc)

  if (plot) Lratio_BHI_plot(mlbin, CAL, LSQ, Lref, Data, x, TAC, Cc, yrsmth)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Lratio_BHI2) <- "MP"


#' @templateVar mp Lratio_BHI3
#' @template MPuses
#'
#' @describeIn Lratio_BHI A modified version of Lratio_BHI2 where mean length
#' is calculated for lengths > modal length (Lc)
#' @export
#' @examples
#' Lratio_BHI3(1, Data=MSEtool::SimulatedData, plot=TRUE)
#'
Lratio_BHI3 <- function(x, Data, reps=100, plot=FALSE, yrsmth = 3) {

  dependencies = "Data@vb_Linf, Data@CV_vbLinf, Data@Cat, Data@CV_Cat, Data@Mort, Data@CV_Mort,
  Data@vb_K, Data@CV_vbK, Data@FMSY_M, Data@CV_FMSY_M, Data@CAL, Data@CAL_bins,
  Data@LFS, Data@CV_LFS"
  Linfc <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Lc <- MSEtool::trlnorm(reps, Data@LFS[x], Data@CV_LFS[x])

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  gamma <- MSEtool::trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x])
  theta <- Kvec/Mvec

  Lref <- (theta * Linfc + Lc * (gamma + 1)) / (gamma + theta + 1)

  Cat <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps, Cat, Data@CV_Cat[x,1])

  LYear <- dim(Data@CAL)[2]
  nlbin <- ncol(Data@CAL[x,,])
  mlbin <- 0.5 * (Data@CAL_bins[1:nlbin] + Data@CAL_bins[2:(nlbin + 1)])

  ind.year <- (LYear - yrsmth + 1):LYear
  CAL <- colSums(Data@CAL[x, ind.year, ])

  CAL <- CAL[mlbin >= mean(Lc, na.rm=TRUE)]
  mlbin <- mlbin[mlbin >= mean(Lc, na.rm=TRUE)]

  LSQ <- rep(NA, reps)
  for(i in 1:reps) {
    lensamp <- sample(mlbin, 0.5*sum(CAL), replace = T, prob = CAL)
    LSQ[i] <- mean(lensamp, na.rm=TRUE)
  }

  TAC <- MSEtool::TACfilter((LSQ/Lref) * Cc)

  if (plot) Lratio_BHI_plot(mlbin, CAL, LSQ, Lref, Data, x, TAC, Cc, yrsmth)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Lratio_BHI3) <- "MP"




#### Stepwise Constant Catch ####

LstepCC_ <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, stepsz = 0.05,
                     llim = c(0.96, 0.98, 1.05)) {
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)  # recent 5 years
  ylast <- (Data@LHYear[1] - Data@Year[1]) + 1  #last historical year
  ind2 <- ((ylast - (yrsmth - 1)):ylast)  # historical 5 pre-projection years
  ind3 <- ((ylast - (yrsmth * 2 - 1)):ylast)  # historical 10 pre-projection years
  C_dat <- Data@Cat[x, ind2]
  if (is.na(Data@MPrec[x]) || length(Data@Year) == ylast) {
    TACstar <- (1 - xx) * MSEtool::trlnorm(reps, mean(C_dat), Data@CV_Cat/(yrsmth^0.5))
  } else {
    TACstar <- rep(Data@MPrec[x], reps)
  }
  step <- stepsz * TACstar
  Lrecent <- mean(Data@ML[x,ind], na.rm=TRUE)
  Lave <- mean(Data@ML[x,ind3], na.rm=TRUE)
  rat <- Lrecent/Lave
  if (rat < llim[1]) {
    TAC <- TACstar - 2 * step
  } else if (rat < llim[2]) {
    TAC <- TACstar - step
  } else if (rat > llim[3]) {
    TAC <- TACstar + step
  } else {
    TAC <- TACstar
  }

  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(op)
    par(mfrow=c(1,2))
    ylim <- range(c(Data@ML[x,]))
    plot(Data@Year, Data@ML[x,], ylab="Mean length", xlab="Year", bty="l", lwd=2,
         type="l", ylim=ylim)
    lines(Data@Year[ind], rep(Lrecent, length(ind)), lty=2)
    text(max(Data@Year[ind]), Lrecent, "Lrecent", xpd=NA, pos=3)

    lines(Data@Year[ind3], rep(Lave, length(ind3)), lty=2, col='blue')
    text(max(Data@Year[ind3]), Lave, "Lave", xpd=NA, pos=3, col='blue')


    ylim <- range(c(Data@Cat[x,], TAC, TACstar ))

    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x,],NA), type="l", xlab="Year",
         ylab=paste0("Catch (", Data@Units, ")"),
         lwd=2,  bty="l", ylim=ylim)

    boxplot(TAC, col="blue", axes=FALSE, at=max(Data@Year)+1, add=TRUE)
    text(max(Data@Year)+1, quantile(TAC, 0.95, na.rm=TRUE), "TAC", col="blue", pos=4, xpd=NA)
    points(max(Data@Year), mean(TACstar, na.rm=TRUE), cex=1, pch=16, col="green")
    text(max(Data@Year), mean(TACstar, na.rm=TRUE), "TAC*", col="green", pos=1, xpd=NA)
  }
  list(TAC=TAC)
}

#' Step-wise Constant Catch
#'
#' A management procedure that incrementally adjusts the TAC according to the
#' mean length of recent catches.
#'
#' The TAC is calculated as:
#'   \deqn{\textrm{TAC} =
#'              \left\{\begin{array}{ll}
#'              \textrm{TAC}^* - 2 S\textrm{TAC}^* & \textrm{if } r < 0.96  \\
#'              \textrm{TAC}^* - S \textrm{TAC}^* & \textrm{if } r < 0.98 \\
#'              \textrm{TAC}^* & \textrm{if } > 1.058 \\
#'              \end{array}\right.
#'            }
#' where \eqn{\textrm{TAC}^*} is (1-`xx`) times average catch in the first year,
#' and previous catch in all projection years, \eqn{S} is step-size determined by `stepsz`,
#' and \eqn{r} is the ratio of \eqn{L_\textrm{recent}} and \eqn{L_\textrm{ave}}
#' which are mean length over the most recent `yrsmth`  years and 2 x `yrsmth` historical
#' years respectively.
#'
#' The conditions are specified in the `llim` argument to the function.
#'
#' @templateVar mp LstepCC1
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth Years over which to calculate mean length.
#' @param xx Parameter controlling the fraction of mean catch to start using in
#' first year
#' @param stepsz Parameter controlling the size of update increment in TAC or effort.
#' @param llim A vector of length reference points that determine the
#' conditions for increasing, maintaining or reducing the TAC or effort.
#'
#' @author T. Carruthers
#' @references Carruthers et al. 2015. Performance evaluation of simple
#' management procedures. ICES J. Mar Sci. 73, 464-482.
#'
#' Geromont, H.F., Butterworth, D.S. 2014. Generic management procedures for
#' data-poor fisheries; forecasting with few data. ICES J. Mar. Sci.
#' doi:10.1093/icesjms/fst232
#' @describeIn LstepCC1 The least biologically precautionary TAC-based MP.
#' @export
#' @examples
#' LstepCC1(1, Data=MSEtool::SimulatedData, plot=TRUE)
#'
LstepCC1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0, stepsz = 0.05,
                     llim = c(0.96, 0.98, 1.05)) {
  runLstep <- LstepCC_(x, Data, reps, plot, yrsmth, xx, stepsz, llim)
  Rec <- new("Rec")
  Rec@TAC <- runLstep$TAC
  Rec
}
class(LstepCC1) <- "MP"


#' @templateVar mp LstepCC2
#' @template MPuses
#' @describeIn LstepCC1 More biologically precautionary than `LstepCC1` (`xx` = 0.1)
#' @export
#' @examples
#' LstepCC2(1, Data=MSEtool::SimulatedData, plot=TRUE)
LstepCC2 <- LstepCC1
formals(LstepCC2)$xx <- 0.1
class(LstepCC2) <- "MP"

#' @templateVar mp LstepCC3
#' @template MPuses
#' @describeIn LstepCC1 More biologically precautionary than `LstepCC2` (`xx` = 0.2)
#' @export
#' @examples
#' LstepCC3(1, Data=MSEtool::SimulatedData, plot=TRUE)
LstepCC3 <- LstepCC1
formals(LstepCC3)$xx <- 0.2
class(LstepCC3) <- "MP"

#' @templateVar mp LstepCC4
#' @template MPuses
#' @describeIn LstepCC1 The most precautionary TAC-based MP.
#' @export
#' @examples
#' LstepCC4(1, Data=MSEtool::SimulatedData, plot=TRUE)
LstepCC4 <- LstepCC1
formals(LstepCC4)$xx <- 0.3
class(LstepCC4) <- "MP"



## Mean Catch Depletion ####

#' Mean Catch Depletion
#'
#' A simple average catch-depletion MP that was included to demonstrate just
#' how informative an estimate of current stock depletion can be.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = 2 \bar{C} D }
#' where \eqn{\bar{C}} is mean historical catch, and \eqn{D} is estimate of current
#' depletion.
#'
#' The TAC is modified by a harvest control rule in `MCD4010`.
#'
#' @templateVar mp MCD
#' @template MPtemplate
#' @template MPuses
#'
#' @describeIn MCD The calculated TAC = 2 \* depletion \* AvC
#' @author T. Carruthers
#' @export
#' @examples
#' MCD(1, MSEtool::Atlantic_mackerel, plot=TRUE)
MCD <- function(x, Data, reps = 100, plot=FALSE) {
  # Daft method to demonstrate the relative value of information of
  # current depletion
  dependencies = "Data@Dep, Data@CV_Dep, Data@Cat"
  if (all(is.na(Data@Cat[x, ]))) stop("No catch data")
  if (MSEtool::NAor0(Data@CV_Cat[x,1])) stop("Data@CV_Cat is NA")
  if (MSEtool::NAor0(Data@Dep[x])) stop("Data@Dep is NA")
  if (MSEtool::NAor0(Data@CV_Dep[x])) stop("Data@CV_Dep is NA")
  depo <- max(0.01, min(0.99, Data@Dep[x]))  # known depletion is between 1% and 99% - needed to generalise the Dick and MacCall method to extreme depletion scenarios
  Bt_K <- rbeta(reps * 100, MSEtool::alphaconv(depo, min(depo * Data@CV_Dep[x],
                                                (1 - depo) * Data@CV_Dep[x])), MSEtool::betaconv(depo, min(depo * Data@CV_Dep[x],
                                                                                                  (1 - depo) * Data@CV_Dep[x])))  # CV 0.25 is the default for Dick and MacCall mu=0.4, sd =0.1
  Bt_K <- Bt_K[Bt_K > 0.00999 & Bt_K < 0.99001][1:reps]  # interval censor (0.01,0.99)  as in Dick and MacCall 2011

  if (reps > 1) {
    AvC <- stats::rlnorm(reps, log(mean(Data@Cat[x, ], na.rm = T)), Data@CV_Cat[x,1])
  } else {
    AvC <- mean(Data@Cat[x, ], na.rm = T)
  }
  TAC <- MSEtool::TACfilter(AvC * 2 * Bt_K)

  if (plot) MCD_plot(Data, AvC, Bt_K, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(MCD) <- "MP"


#' @templateVar mp MCD4010
#' @template MPuses
#' @describeIn MCD Linked to a 40-10 harvest control rule
#' @author T. Carruthers
#' @export
#' @examples
#' MCD4010(1, MSEtool::Atlantic_mackerel, plot=TRUE)
MCD4010 <- function(x, Data, reps = 100, plot=FALSE) {
  # Daft method to demonstrate the relative value of information of
  # current depletion
  dependencies = "Data@Dep, Data@CV_Dep, Data@Cat"
  if (all(is.na(Data@Cat[x, ]))) stop("No catch data")
  if (MSEtool::NAor0(Data@CV_Cat[x,1])) stop("Data@CV_Cat is NA")
  if (MSEtool::NAor0(Data@CV_Dep[x])) stop("Data@CV_Dep is NA")
  if (is.na(Data@Dep[x])) {
    Rec <- new("Rec")
    Rec@TAC <- rep(NA, reps)
    return(Rec)
  }
  depo <- max(0.01, min(0.99, Data@Dep[x]))  # known depletion is between 1% and 99% - needed to generalise the Dick and MacCall method to extreme depletion scenarios
  Bt_K <- rbeta(reps * 100, MSEtool::alphaconv(depo, min(depo * Data@CV_Dep[x],
                                                (1 - depo) * Data@CV_Dep[x])), MSEtool::betaconv(depo, min(depo * Data@CV_Dep[x],
                                                                                                  (1 - depo) * Data@CV_Dep[x])))  # CV 0.25 is the default for Dick and MacCall mu=0.4, sd =0.1
  Bt_K <- Bt_K[Bt_K > 0.00999 & Bt_K < 0.99001][1:reps]  # interval censor (0.01,0.99)  as in Dick and MacCall 2011
  if (reps > 1) {
    AvC <- stats::rlnorm(reps, log(mean(Data@Cat[x, ], na.rm = T)), Data@CV_Cat[x,1])
  } else {
    AvC <- mean(Data@Cat[x, ], na.rm = T)
  }
  TAC <- AvC * 2 * Bt_K

  # 40-10 HCR
  cond1 <- Bt_K < 0.4 & Bt_K > 0.1
  cond2 <- Bt_K < 0.1
  TAC[cond1] <- TAC[cond1] * (Bt_K[cond1] - 0.1)/0.3
  TAC[cond2] <- TAC[cond2] * tiny  # this has to still be stochastic albeit very small
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) MCD_plot(Data, AvC, Bt_K, TAC)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(MCD4010) <- "MP"




#### Intrinsic Rate of Increase ####

#' Intrinsic rate of Increase MP
#'
#' An MP proposed by Carl Walters that modifies the TAC according to trends in
#' apparent surplus production that includes information from a demographically
#' derived prior for intrinsic rate of increase
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = \textrm{SP} (1-gG)}
#' where \eqn{g} is a gain parameter, \eqn{\textrm{SP}} is estimated surplus production,
#' and \eqn{G} is:
#' For `Rcontrol`:  \eqn{G = r (1-2D)} where \eqn{r} is the estimated intrinsic rate
#' of increase, and \eqn{D} is assumed depletion.
#'
#' For `Rcontrol2`: \eqn{G = r - 2bB_\textrm{hist}} where \eqn{B_\textrm{hist}}
#' is the smoothed biomass overlast `yrsmth` years and:
#' \deqn{b = \sum{\frac{\textrm{SP}}{B_\textrm{hist}} - r} \frac{\sum{B_\textrm{hist}}}{\sum{B_\textrm{hist}^2}}   }.
#'
#'
#' The TAC is subject to conditions limit the maximum change from the smoothed catch
#' over the last `yrsmth` years by the `glim` argument, e.g, default values of `glim = c(0.5, 2)`
#' means that maximum decrease in TAC is 50% of average catch and maximum increase
#'  is 2 x average catch.
#'
#' @templateVar mp Rcontrol
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth The number of years for smoothing catch and biomass data
#' @param gg A gain parameters
#' @param glim Limits for the change in TAC among years
#'
#' @author C. Walters and T. Carruthers
#' @references Made-up for this package.
#' @export
#' @family Surplus production MPs
#' @describeIn Rcontrol Base version `Rcontrol`
#' @examples
#' Rcontrol(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
Rcontrol <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2, glim = c(0.5, 2)) {
  dependencies = "Data@Mort, Data@CV_Mort, Data@vbK, Data@CV_vbK, Data@vbLinf, Data@CV_vbLinf, Data@vbt0, Data@CV_vbt0 Data@steep, Data@CV_steep, Data@MaxAge, Data@Dep, Data@CV_Dep, Data@Cat, Data@Ind"
  if (MSEtool::NAor0(Data@Mort[x])) stop("Data@Mort is NA")
  if (MSEtool::NAor0(Data@FMSY_M[x])) stop("Data@FMSY_M is NA")
  if (MSEtool::NAor0(Data@L50[x])) stop("Data@L50 is NA")
  if (MSEtool::NAor0(Data@vbK[x])) stop("Data@vbK is NA")
  if (MSEtool::NAor0(Data@CV_vbK[x])) stop("Data@CV_vbK is NA")
  if (MSEtool::NAor0(Data@vbLinf[x])) stop("Data@vbLinf is NA")
  if (MSEtool::NAor0(Data@CV_vbLinf[x])) stop("Data@CV_vbLinf is NA")
  if (MSEtool::NAor0(Data@vbt0[x])) stop("Data@vbt0 is NA")
  if (MSEtool::NAor0(Data@wla[x])) stop("Data@wla is NA")
  if (MSEtool::NAor0(Data@wlb[x])) stop("Data@wlb is NA")
  if (MSEtool::NAor0(Data@CV_Mort[x])) stop("Data@CV_Mort is NA")

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfvec <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0vec <- -MSEtool::trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0vec <- rep(Data@vbt0[x], reps)
  }
  t0vec[!is.finite(t0vec)] <- 0

  # hvec <- MSEtool::trlnorm(reps, Data@steep[x], Data@CV_steep[x])
  hvec <- MSEtool::sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])
  rsamp <- getr(x, Data, Mvec, Kvec, Linfvec, t0vec, hvec, maxage = Data@MaxAge,
                r_reps = reps)

  depo <- max(0.01, min(0.99, Data@Dep[x]))  # known depletion is between 1% and 99% - needed to generalise the Dick and MacCall method to extreme depletion scenarios
  if (any(is.na(c(Data@Dep[x], Data@CV_Dep[x])))) {
    rec <- new("Rec")
    rec@TAC <- rep(as.numeric(NA), reps)
    return(rec)
  }

  if (reps > 1) {
    Bt_K <- rbeta(100, MSEtool::alphaconv(depo, min(depo * Data@CV_Dep[x],
                                           (1 - depo) * Data@CV_Dep[x])),
                  MSEtool::betaconv(depo, min(depo * Data@CV_Dep[x], (1 - depo) * Data@CV_Dep[x])))  # CV 0.25 is the default for Dick and MacCall mu=0.4, sd =0.1
  } else {
    Bt_K <- depo
  }
  Bt_K <- Bt_K[Bt_K > 0.01 & Bt_K < 0.99][1:reps]  # interval censor (0.01,0.99)  as in Dick and MacCall 2011

  G_new <- rsamp * (1 - 2 * Bt_K)  # here is a big difference from SPHCR

  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  C_dat <- log(Data@Cat[x, ind])
  C_dat[C_dat == -Inf] <- 0
  B_dat <- log(Data@Ind[x, ind]/Data@Ind[x, ind[yrsmth]] * Data@Abun[x])
  B_dat[B_dat == -Inf] <- 0
  C_hist <- exp(predict(loess(C_dat ~ ind, degree = 1)))
  B_hist <- exp(predict(loess(B_dat ~ ind, degree = 1)))
  inda <- 2:yrsmth
  ind1 <- 1:(yrsmth - 1)
  SP_hist <- B_hist[inda] - B_hist[ind1] + C_hist[ind1]
  yind <- 1:length(SP_hist)
  SP_mu <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1))
  SP_se <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) + 1), se = T)$se.fit
  if (reps >1) {
    SP_new <- rnorm(reps, SP_mu, SP_se/2)
  } else {
    SP_new <- SP_mu
  }
  TACa <- SP_new * (1 - gg * G_new)
  TAC <- TACa
  TAC[TAC < glim[1] * C_hist[yrsmth]] <- glim[1] * C_hist[yrsmth]
  TAC[TAC > glim[2] * C_hist[yrsmth]] <- glim[2] * C_hist[yrsmth]

  TAC <- MSEtool::TACfilter(TAC)
  if (plot) Rcontrol_plot(rsamp, ind, G_new, B_hist, SP_hist, SP_mu, B_dat, Data, yind, C_hist, TAC, TACa)

    # Carr<-cbind(array(rep(Data@Cat[x,],each=reps),c(reps,length(Data@Cat[x,]))),TAC)
  # Warr<-(Data@Mort[x]*exp(-Data@Mort[x]*(1:ncol(Carr))))[ncol(Carr):1]
  # Warr<-Warr/sum(Warr)
  # TAC<-apply(t(matrix(Warr,nrow=ncol(Carr),ncol=reps))*Carr,1,sum)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Rcontrol) <- "MP"


#' @describeIn Rcontrol This is different from
#' `Rcontrol` because it includes a quadratic approximation of recent trend in
#' surplus production given biomass
#' @export
#' @examples
#' Rcontrol2(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
Rcontrol2 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2, glim = c(0.5, 2)) {
  dependencies = "Data@Mort, Data@CV_Mort, Data@vbK, Data@CV_vbK, Data@vbLinf, Data@CV_vbLinf, Data@vbt0, Data@CV_vbt0, Data@steep, Data@CV_steep, Data@MaxAge, Data@Dep, Data@CV_Dep, Data@Cat, Data@Ind"
  if (MSEtool::NAor0(Data@Mort[x])) stop("Data@Mort is NA")
  if (MSEtool::NAor0(Data@FMSY_M[x])) stop("Data@FMSY_M is NA")
  if (MSEtool::NAor0(Data@L50[x])) stop("Data@L50 is NA")
  if (MSEtool::NAor0(Data@vbK[x])) stop("Data@vbK is NA")
  if (MSEtool::NAor0(Data@CV_vbK[x])) stop("Data@CV_vbK is NA")
  if (MSEtool::NAor0(Data@vbLinf[x])) stop("Data@vbLinf is NA")
  if (MSEtool::NAor0(Data@CV_vbLinf[x])) stop("Data@CV_vbLinf is NA")
  if (MSEtool::NAor0(Data@vbt0[x])) stop("Data@vbt0 is NA")
  if (MSEtool::NAor0(Data@CV_Mort[x])) stop("Data@CV_Mort is NA")
  if (MSEtool::NAor0(Data@wla[x])) stop("Data@wla is NA")
  if (MSEtool::NAor0(Data@wlb[x])) stop("Data@wlb is NA")

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfvec <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0vec <- -MSEtool::trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0vec <- rep(Data@vbt0[x], reps)
  }
  t0vec[!is.finite(t0vec)] <- 0
  # hvec <- MSEtool::trlnorm(reps, Data@steep[x], Data@CV_steep[x])
  hvec <- MSEtool::sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])
  rsamp <- getr(x, Data, Mvec, Kvec, Linfvec, t0vec, hvec, maxage = Data@MaxAge,
                r_reps = reps)

  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  C_dat <- log(Data@Cat[x, ind])
  C_dat[C_dat == -Inf] <- 0
  B_dat <- log(Data@Ind[x, ind]/Data@Ind[x, ind[yrsmth]] * Data@Abun[x])
  B_dat[B_dat == -Inf] <- 0
  C_hist <- exp(predict(loess(C_dat ~ ind, degree = 1)))
  B_hist <- exp(predict(loess(B_dat ~ ind, degree = 1)))
  inda <- 2:yrsmth
  ind1 <- 1:(yrsmth - 1)
  SP_hist <- B_hist[inda] - B_hist[ind1] + C_hist[ind1]
  yind <- 1:length(SP_hist)
  SP_mu <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) +
                                                       1))
  SP_se <- predict(lm(SP_hist ~ yind), newdat = list(yind = length(SP_hist) +
                                                       1), se = T)$se.fit
  if (reps >1) {
    SP_new <- rnorm(reps, SP_mu, SP_se/2)
  } else {
    SP_new <- SP_mu
  }
  SParr <- array(rep(SP_hist, each = reps), dim = c(reps, yrsmth - 1))
  Barr <- array(rep(B_hist[inda], each = reps), dim = c(reps, yrsmth - 1))
  rarr <- array(rep(rsamp, yrsmth - 1), dim = c(reps, yrsmth - 1))
  b2 <- apply(SParr/Barr - rarr, 1, sum) * apply(Barr, 1, sum)/apply(Barr^2, 1, sum)
  G_new <- rsamp - 2 * b2 * B_hist[yrsmth]

  TACa <- SP_new * (1 - gg * G_new)
  TAC <- TACa
  TAC[TAC < glim[1] * C_hist[yrsmth]] <- glim[1] * C_hist[yrsmth]
  TAC[TAC > glim[2] * C_hist[yrsmth]] <- glim[2] * C_hist[yrsmth]
  # Carr<-cbind(array(rep(Data@Cat[x,],each=reps),c(reps,length(Data@Cat[x,]))),TAC)
  # Warr<-(Data@Mort[x]*exp(-Data@Mort[x]*(1:ncol(Carr))))[ncol(Carr):1]
  # Warr<-Warr/sum(Warr)
  # TAC<-apply(t(matrix(Warr,nrow=ncol(Carr),ncol=reps))*Carr,1,sum)
  TAC <- MSEtool::TACfilter(TAC)
  if (plot) Rcontrol_plot(rsamp, ind, G_new, B_hist, SP_hist, SP_mu, B_dat, Data, yind, C_hist, TAC, TACa)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(Rcontrol2) <- "MP"



#### SBT ####

#' SBT simple MP
#'
#' An MP that makes incremental adjustments to TAC recommendations based on the
#' apparent trend in CPUE, a an MP that makes incremental adjustments to TAC
#' recommendations based on index levels relative to target levels (BMSY/B0)
#' and catch levels relative to target levels (MSY).
#'
#' For `SBT1` the TAC is calculated as:
#'   \deqn{\textrm{TAC}_y =
#'              \left\{\begin{array}{ll}
#'               C_{y-1} (1+K_2 \left| \lambda \right| ) & \textrm{if } \lambda \geq 0 \\
#'               C_{y-1} (1-K_1 \left| \lambda \right| ^\gamma) & \textrm{if } \lambda < 0\\
#'              \end{array}\right.
#'            }
#'  where \eqn{\lambda} is the slope of index over the last `yrmsth` years, and
#'  \eqn{K_1}, \eqn{K_2}, and \eqn{\gamma} are arguments to the MP.
#'
#' For `SBT2` the TAC is calculated as:
#' \deqn{\textrm{TAC}_y = 0.5 (C_{y-1} + C_\textrm{targ}\delta)}
#' where \eqn{C_{y-1}} is catch in the previous year, \eqn{C_{\textrm{targ}}}
#' is a target catch (`Data@Cref`), and :
#'   \deqn{\delta=
#'              \left\{\begin{array}{ll}
#'              R^{1-\textrm{epsR}} & \textrm{if } R \geq 1 \\
#'              R^{1+\textrm{epsR}} & \textrm{if } R < 1 \\
#'              \end{array}\right.
#'            }
#' where \eqn{\textrm{epsR}} is a control parameter and:
#' \eqn{R = \frac{\bar{r}}{\phi}}
#' where \eqn{\bar{r}} is mean recruitment over last `tauR` years and \eqn{\phi}
#' is mean recruitment over last 10 years.
#'
#' This isn't exactly the same as the proposed methods and is stochastic in
#' this implementation. The method doesn't tend to work too well under many
#' circumstances possibly due to the lack of 'tuning' that occurs in the real
#' SBT assessment environment. You could try asking Rich Hillary at CSIRO about
#' this approach.
#'
#' @templateVar mp SBT1
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth The number of years for evaluating trend in relative abundance
#' indices
#' @param k1 Control parameter
#' @param k2 Control parameter
#' @param gamma Control parameter
#'
#' @author T. Carruthers
#' @references http://www.ccsbt.org/site/recent_assessment.php
#' @family SBT MPs
#' @export
#' @describeIn SBT1 Simple SBT MP
#' @examples
#' SBT1(1, Data=MSEtool::SimulatedData, plot=TRUE)
SBT1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, k1 = 1.5, k2 = 3, gamma = 1) {
  dependencies = "Data@Cat, Data@Year, Data@Ind"
  Cr <- length(Data@Cat[x, ])
  cct <- MSEtool::trlnorm(reps, Data@Cat[x, Cr], Data@CV_Cat[x,1])
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  I_hist <- Data@Ind[x, ind]
  test <- summary(lm(I_hist ~ ind))$coefficients[2, 1:2]
  if (reps > 1) {
    lambda <- rnorm(reps, test[1], test[2])
  } else {
    lambda <- test[1]
  }
  # TAC <- cct * 1 + k2 * lambda
  # see https://github.com/DLMtool/DLMtool/issues/17
  TAC <- cct * (1 + k2 * lambda)
  cond <- lambda < 0
  # TAC[cond] <- cct[cond] * 1 - k1 * -lambda[cond]^gamma
  TAC[cond] <- cct[cond] * (1 - k1 * -lambda[cond]^gamma)
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    par(mfrow=c(1,2))
    plot(Data@Year[ind], I_hist, bty="l", xlab="Year", ylab="Index", type="l", lwd=2)
    lines(Data@Year[ind], predict(lm(I_hist ~ ind)), lty=2)

    ylim <- range(c(Data@Cat[x, ind], TAC))
    plot(c(Data@Year[ind], max(Data@Year[ind])+1), c(Data@Cat[x, ind], NA), bty="l", xlab="Year",
         ylab=paste0("Catch (", Data@Units, ")"), type="l", lwd=2, ylim=ylim)
    points(max(Data@Year[ind]), Data@Cat[x, max(ind)], pch=16, cex=1.5)
    boxplot(TAC, at=max(Data@Year[ind])+1, add=TRUE, axes=FALSE, col="blue")

  }


  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SBT1) <- "MP"







#' @templateVar mp SBT2
#' @template MPuses
#'
#' @param epsR Control parameter
#' @param tauR Control parameter
#'
#' @export
#' @describeIn SBT1 Complex SBT MP
#' @examples
#' SBT2(1, Data=MSEtool::SimulatedData, plot=TRUE)
SBT2 <- function(x, Data, reps = 100, plot=FALSE, epsR = 0.75, tauR = 5,
                 gamma = 1) {
  dependencies = "Data@Cref, Data@Rec, Data@Cat"

  Ctarg <- MSEtool::trlnorm(reps, Data@Cref[x], Data@CV_Cref)
  muR <- mean(Data@Rec[x, (length(Data@Rec[x, ]) - tauR + 1):length(Data@Rec[x, ])], na.rm=TRUE)
  phi <- mean(Data@Rec[x, (length(Data@Rec[x, ]) - 9):length(Data@Rec[x,])], na.rm=TRUE)
  Rrat <- muR/phi
  deltaR <- rep(NA, reps)
  deltaR[Rrat > 1] <- Rrat[Rrat > 1]^(1 - epsR)
  deltaR[Rrat < 1] <- Rrat[Rrat < 1]^(1 + epsR)
  TAC <- 0.5 * (Data@Cat[x, length(Data@Cat[x, ])] + Ctarg *  deltaR)

  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    par(mfrow=c(1,2))
    plot(Data@Year, Data@Rec[x,], bty="l", xlab="Year",
         ylab="Recruitment", type="l", lwd=2)
    nyrs <- length(Data@Year)
    nyrs2 <- length((nyrs - tauR + 1) : nyrs)
    lines(Data@Year[(nyrs - tauR + 1) : nyrs], rep(muR, nyrs2), lty=2, col="blue")
    text(mean(Data@Year[(nyrs - tauR + 1) : nyrs]), muR, col="blue", 'muR', pos=4, xpd=NA)

    nyrs2 <- length((nyrs - 9) : nyrs)
    lines(Data@Year[(nyrs - 9 ) : nyrs], rep(phi, nyrs2), lty=2, col="red")
    text(mean(Data@Year[(nyrs - 9) : nyrs]), phi, col="red", 'phi', pos=2, xpd=NA)

    ylim <- range(c(Data@Cat[x, ], TAC, Ctarg))
    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x, ], NA), bty="l", xlab="Year",
         ylab=paste0("Catch (", Data@Units, ")"), type="l", lwd=2, ylim=ylim)
    abline(h=Ctarg, lty=2, col="darkgray")
    boxplot(TAC, at=max(Data@Year)+1, add=TRUE, axes=FALSE, col="blue")
  }

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SBT2) <- "MP"




##### Surplus Production based MP ####

#' Surplus production based catch-limit modifier
#'
#' An MP that makes incremental adjustments to TAC recommendations based on the
#' apparent trend in surplus production. Based on the theory of Mark Maunder
#' (IATTC)
#'
#' Note that this isn't exactly what Mark has previously suggested and is
#' stochastic in this implementation.
#'
#' The TAC is calculated as:
#'   \deqn{\textrm{TAC}_y =
#'              \left\{\begin{array}{ll}
#'               C_{y-1} \textrm{bet}_1 & \textrm{if } r < \alpha_1 \\
#'               C_{y-1} & \textrm{if } \alpha_1 < r < \alpha_2 \\
#'               \textrm{bet}_2 (b_2 - b_1 + C_{y-2} )  & \textrm{if } r > \alpha_2 \\
#'              \end{array}\right.
#'            }
#'  where \eqn{\textrm{bet}_1} and \eqn{\textrm{bet}_2} are elements in `bet`,
#'  \eqn{r} is the ratio of the index in the most recent two years, \eqn{C_{y-1}}
#'   is catch in the previous year, \eqn{b_1} and \eqn{b_2} are ratio of index
#'   in \eqn{y-2} and \eqn{y-1} over the estimate of catchability \eqn{\left(\frac{I}{A}\right)},
#'   and \eqn{\alpha_1}, \eqn{\alpha_2}, and \eqn{\alpha_3} are specified in argument
#'   `alp`.
#'
#'
#' @templateVar mp SPmod
#' @template MPtemplate
#' @template MPuses
#'
#' @param alp Condition for modifying the TAC (bounds on change in abundance)
#' @param bet Limits for how much the TAC can change among years
#' @return A numeric vector of TAC recommendations
#' @author T. Carruthers
#' @references
#' http://www.iattc.org/Meetings/Meetings2014/MAYSAC/PDFs/SAC-05-10b-Management-Strategy-Evaluation.pdf
#' @family Surplus production MPs
#' @export
#' @examples
#' SPmod(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
SPmod <- function(x, Data, reps = 100, plot=FALSE, alp = c(0.8, 1.2), bet = c(0.8, 1.2)) {
  dependencies = "Data@Cat, Data@Ind, Data@Abun, Data@CV_Ind, Data@CV_Cat,  Data@CV_Abun"
  Ir <- length(Data@Ind[x, ])
  Cr <- length(Data@Cat[x, ])
  rat <- MSEtool::trlnorm(reps, Data@Ind[x, Ir], Data@CV_Ind[x,1])/MSEtool::trlnorm(reps, Data@Ind[x, Ir - 1], Data@CV_Ind[x,1])
  cct <- MSEtool::trlnorm(reps, Data@Cat[x, Cr], Data@CV_Cat[x,1])
  Abun <- MSEtool::trlnorm(reps, Data@Abun[x], Data@CV_Abun[x])
  TAC <- rep(NA, reps)
  TAC[rat < alp[1]] <- cct[rat < alp[1]] * bet[1]
  TAC[rat > alp[1] & rat < alp[2]] <- cct[rat > alp[1] & rat < alp[2]]

  cond <- rat > alp[2]
  reps2 <- sum(cond)
  if (reps2 > 0) {
    qq1 <- MSEtool::trlnorm(reps2, Data@Ind[x, Ir]/Abun[cond], Data@CV_Ind[x,1])
    bio1 <- Data@Ind[x, Ir - 1]/qq1
    bio2 <- Data@Ind[x, Ir]/qq1
    cct1 <- MSEtool::trlnorm(reps2, Data@Cat[x, Cr - 1], Data@CV_Cat[x,1])
    PP <- bio2 - bio1 + cct1
    TAC[cond] <- bet[2] * PP
  }
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    par(mfrow=c(1,2))
    plot(Data@Year[(Ir-1):Ir], Data@Ind[x,(Ir-1):Ir], bty="l", xlab="Year",
         ylab="Index", type="l", lwd=2)

    ylim <- range(c(Data@Cat[x, ], TAC), na.rm=TRUE)
    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x, ], NA), bty="l", xlab="Year",
         ylab=paste0("Catch (", Data@Units, ")"), type="l", lwd=2, ylim=ylim)
    boxplot(TAC, at=max(Data@Year)+1, add=TRUE, axes=FALSE, col="blue")
  }

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SPmod) <- "MP"




#' Slope in surplus production MP
#'
#' A management procedure that makes incremental adjustments to TAC
#' recommendations based on the apparent trend in recent surplus production.
#' Based on the theory of Mark Maunder (IATTC)
#'
#' Note that this isn't exactly what Mark has previously suggested and is
#' stochastic in this implementation.
#'
#' The TAC is calculated as:
#'   \deqn{\textrm{TAC}_y =
#'              \left\{\begin{array}{ll}
#'               M \bar{C}  & \textrm{if } r < \alpha_1 \\
#'               \bar{C} & \textrm{if } \alpha_1 < r < \alpha_2 \\
#'               \textrm{bet}_2 \textrm{SP}  & \textrm{if } r > \alpha_2 \\
#'              \end{array}\right.
#'            }
#'  where \eqn{r} is the ratio of predicted biomass in next year to biomass in
#'  current year \eqn{\bar{C}} is the mean catch over the last `yrmsth` years, \eqn{\alpha_1}
#'  and  \eqn{\alpha_2} are specified in `alp`, \eqn{\textrm{bet}_1} and \eqn{\textrm{bet}_2}
#'  are specified in `bet`, \eqn{\textrm{SP}} is estimated surplus production in most recent year,
#'  and:
#'  \deqn{M = 1-\textrm{bet}_1 \frac{B_y - \tilde{B}_y}{B_y}}
#'  where \eqn{B_y} is the most recent estimate of biomass and \eqn{\tilde{B}}
#'  is the predicted biomass in the next year.
#'
#' @templateVar mp SPslope
#' @template MPtemplate
#' @template MPuses
#'
#' @param yrsmth Years over which to smooth recent estimates of surplus
#' production
#' @param alp Condition for modifying the Data (bounds on change in
#' abundance)
#' @param bet Limits for how much the Data can change among years
#' @author T. Carruthers
#' @family Surplus production MPs
#' @references
#' http://www.iattc.org/Meetings/Meetings2014/MAYSAC/PDFs/SAC-05-10b-Management-Strategy-Evaluation.pdf
#' @export
#' @examples
#' SPslope(1, Data=MSEtool::Atlantic_mackerel, plot=TRUE)
SPslope <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 4, alp = c(0.9, 1.1),
                    bet = c(1.5, 0.9)) {

  dependencies = "Data@Year, Data@Cat, Data@Ind, Data@Abun"
  ind <- (length(Data@Year) - (yrsmth - 1)):length(Data@Year)
  yind <- 1:yrsmth
  C_dat <- Data@Cat[x, ind]
  B_dat <- Data@Ind[x, ind]/Data@Ind[x, ind[yrsmth]] * Data@Abun[x]
  Pt_mu <- max(B_dat[yrsmth] - B_dat[yrsmth - 1] + C_dat[yrsmth - 1], tiny)
  Pt_1 <- MSEtool::trlnorm(reps, Pt_mu, Data@CV_Cat[x,1])
  It <- exp(predict(lm(log(B_dat) ~ yind), newdat = list(yind = yrsmth + 1)))
  Ilast <- B_dat[yrsmth]
  MC <- max(mean(C_dat, na.rm=TRUE), tiny)
  Ct_1 <- MSEtool::trlnorm(reps, MC, Data@CV_Cat[x,1]/(yrsmth^0.5))  # mean catches over the interval

  rat <- It/Ilast

  mult <- max((1 - bet[1] * (Ilast - It)/Ilast), tiny)
  if (rat < alp[1]) TAC <- mult * Ct_1
  if (rat > alp[1] & rat < alp[2]) TAC <- Ct_1
  if (rat > alp[2]) TAC <- bet[2] * Pt_1
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(par(op))
    par(mfrow=c(1,2))
    ylim <- range(c(B_dat, It), na.rm=TRUE)
    plot(c(Data@Year[ind], max(Data@Year[ind])+1), c(B_dat, NA), bty="l", xlab="Year",
         ylab="Estimated Biomass", type="l", lwd=2, ylim=ylim)
    points(max(Data@Year)+1, mean(It, na.rm=TRUE), col="blue", pch=16)

    ylim <- range(c(Data@Cat[x, ], TAC), na.rm=TRUE)
    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x, ], NA), bty="l", xlab="Year",
         ylab=paste0("Catch (", Data@Units, ")"), type="l", lwd=2, ylim=ylim)
    boxplot(TAC, at=max(Data@Year)+1, add=TRUE, axes=FALSE, col="blue")
  }

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SPslope) <- "MP"



#' Catch trend Surplus Production MSY MP
#'
#' An MP that uses Martell and Froese (2012) method for estimating MSY to
#' determine the OFL. Since their approach estimates stock trajectories based
#' on catches and a rule for intrinsic rate of increase it also returns
#' depletion. Given their surplus production model predicts K, r and depletion
#' it is straight forward to calculate the OFL based on the Schaefer productivity
#' curve.
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = D K \frac{r}{2}}
#' where \eqn{D} is depletion, \eqn{K} is unfished biomass, and \eqn{r} is
#' intrinsic rate of increasase, all estimated internally by the method based
#' on trends in the catch data and life-history information.
#'
#' Requires the assumption that catch is proportional to abundance, and a catch
#' time-series from the beginning of exploitation.
#'
#' Occasionally the rule that limits r and K ranges does not allow r-K pairs to
#' be found that lead to the depletion inferred by the catch trajectories. In
#' this case this method widens the search.
#'
#' @templateVar mp SPMSY
#' @template MPtemplate
#' @template MPuses
#'
#' @author T. Carruthers
#' @references Martell, S. and Froese, R. 2012. A simple method for estimating
#' MSY from catch and resilience. Fish and Fisheries. DOI:
#' 10.1111/j.1467-2979.2012.00485.x
#' @family Surplus production MPs
#' @export
#' @examples
#' SPMSY(1, Data=MSEtool::SimulatedData, plot=TRUE)
SPMSY <- function(x, Data, reps = 100, plot=FALSE) {
  # Martell and Froese 2012 Schaefer SP estimate of MSY given priors on
  # r, k and depletion for(x in 1:100){
  dependencies = "Data@MaxAge, Data@vbK, Data@L50, Data@Cat"
  if (all(is.na(Data@Cat[x, ]))) stop("No catch data")
  maxage <- Data@MaxAge
  nsamp <- reps * 200

  # Froese 2012 http://www.fishbase.de/rfroese/Catch-MSY_SummaryFinal.doc
  rule <- rep(4, 3)

  if (Data@vbK[x] > 0.3) {
    # K rules
    rule[1] <- 1
  } else if (Data@vbK[x] < 0.3 & Data@vbK[x] > 0.16) {
    rule[1] <- 2
  } else if (Data@vbK[x] < 0.16 & Data@vbK[x] > 0.05) {
    rule[1] <- 3
  }
  AM <- iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x], Data@L50[x])
  if (AM < 1.5) {
    # Age at maturity rules
    rule[2] <- 1
  } else if (AM < 4.5 & AM > 1.5) {
    rule[2] <- 2
  } else if (AM < 10 & AM > 4.5) {
    rule[2] <- 3
  }

  if (Data@MaxAge < 4) {
    # Maximum age rules
    rule[3] <- 1
  } else if (Data@MaxAge < 11 & Data@MaxAge > 3) {
    rule[3] <- 2
  } else if (Data@MaxAge < 31 & Data@MaxAge > 10) {
    rule[3] <- 3
  }

  if (mean(rule) < 1.5)   rsamp <- runif(nsamp, 0.6, 1.5)
  if (mean(rule) > 1.5 & mean(rule) < 2.5)   rsamp <- runif(nsamp, 0.2, 1)
  if (mean(rule) > 2.5 & mean(rule) < 3.5)  rsamp <- runif(nsamp, 0.05, 0.5)
  if (mean(rule) > 3.5) rsamp <- runif(nsamp, 0.015, 0.1)

  Ksamp <- runif(nsamp, mean(Data@Cat[x, ], na.rm=TRUE)/rsamp, (10 * mean(Data@Cat[x, ], na.rm=TRUE))/rsamp)
  nyears <- length(Data@Cat[x, ])
  B <- array(NA, dim = c(nsamp, nyears))

  if (Data@Cat[x, 1] < (0.5 * max(Data@Cat[x, ]))) {
    # Martell and Froese decision rules (makes absolutely no sense to me!)
    B[, 1] <- Ksamp * runif(nsamp, 0.5, 0.9)
  } else {
    B[, 1] <- Ksamp * runif(nsamp, 0.3, 0.6)
  }

  if (Data@Cat[x, nyears] < (0.5 * max(Data@Cat[x, ]))) {
    # Martell and Froese decision rules (makes absolutely no sense to me!)
    LB <- 0.01
    UB <- 0.4
  } else {
    LB <- 0.3
    UB <- 0.7
  }

  for (i in 2:nyears) {
    B[, i] <- B[, i - 1] - Data@Cat[x, i - 1]
    B[, i] <- B[, i] + rsamp * B[, i] * (1 - B[, i]/Ksamp)
  }
  B <- B/rep(Ksamp, nyears)
  cond <- (B[, nyears] >= LB) & (B[, nyears] <= UB)
  if (sum(cond) <= 1) {
    B[B[, nyears] >= UB, nyears] <- UB
    B[B[, nyears] <= LB, nyears] <- LB
    cond <- (B[, nyears] >= LB) & (B[, nyears] <= UB)
  }
  Best <- B[cond,]
  if (nrow(Best)>reps) {
    Best <- Best[1:reps,]
  }

  dep <- B[cond, nyears][1:reps]
  MSY <- rsamp[cond][1:reps] * Ksamp[cond][1:reps]/4
  Kc <- Ksamp[cond][1:reps]
  rc <- rsamp[cond][1:reps]
  TAC <- Kc * dep * rc/2

  if (sum(!is.na(TAC)) < ceiling(reps/10)) {
    # a fudge of the original method that widens current depletion to the
    # lowest and higest bounds to get an TAC sample
    cond <- (B[, nyears] >= 0.01) & (B[, nyears] <= 0.7)
    dep <- B[cond, nyears][1:reps]
    MSY <- rsamp[cond][1:reps] * Ksamp[cond][1:reps]/4
    Kc <- Ksamp[cond][1:reps]
    rc <- rsamp[cond][1:reps]
    TAC <- Kc * dep * rc/2

  }
  # }
  TAC <- MSEtool::TACfilter(TAC)

  if (plot) {
    op <- par(no.readonly = TRUE)
    on.exit(op)
    par(mfrow=c(2,3), mar=c(4,4,1,1), oma=c(1,1,1,1))

    boxplot(Kc, ylab=paste0('Unfished Biomass (', Data@Units, ')'))
    boxplot(rc, ylab='Intrinsic rate of increase')
    boxplot(dep, ylab='Depletion')
    boxplot(MSY, ylab='MSY')

    matplot(Data@Year, t(Best), xlab="Year", ylab='Relative Biomass',
         bty="l", type="l", lwd=2)

    ylim <- range(c(Data@Cat[x,], TAC), na.rm=TRUE)
    plot(c(Data@Year, max(Data@Year)+1), c(Data@Cat[x,], NA), xlab="Year", ylab=paste0('Catch (', Data@Units, ')'),
         bty="l", type="l", lwd=2, ylim=ylim)
    boxplot(TAC, axes=FALSE, add=TRUE, at=max(Data@Year)+1, col="blue", width=2)

  }
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SPMSY) <- "MP"



#### Surplus Production Stock Reduction Analysis ####

#' Surplus Production Stock Reduction Analysis Internal Function
#'
#' @param x A position in a data-limited methods data object
#' @param Data A data-limited methods data object (class DLM)
#' @param reps The number of samples of the TAC taken for the calculation of
#' the quota
#' @param dep Optional numeric vector of depletion
#' @export
#' @keywords internal
#' @describeIn SPSRA_ internal function
SPSRA_ <- function(x, Data, reps = 100, dep=NULL) {
  if (MSEtool::NAor0(Data@Mort[x])) stop("Data@Mort is NA")
  if (MSEtool::NAor0(Data@FMSY_M[x])) stop("Data@FMSY_M is NA")
  if (MSEtool::NAor0(Data@L50[x])) stop("Data@L50 is NA")
  if (MSEtool::NAor0(Data@Dep[x])) stop("Data@Dep is NA")
  if (MSEtool::NAor0(Data@vbK[x])) stop("Data@vbK is NA")
  if (MSEtool::NAor0(Data@CV_vbK[x])) stop("Data@CV_vbK is NA")
  if (MSEtool::NAor0(Data@vbLinf[x])) stop("Data@vbLinf is NA")
  if (MSEtool::NAor0(Data@CV_vbLinf[x])) stop("Data@CV_vbLinf is NA")
  if (MSEtool::NAor0(Data@vbt0[x])) stop("Data@vbt0 is NA")
  if (MSEtool::NAor0(Data@CV_Mort[x])) stop("Data@CV_Mort is NA")
  if (MSEtool::NAor0(Data@wla[x])) stop("Data@wla is NA")
  if (MSEtool::NAor0(Data@wlb[x])) stop("Data@wlb is NA")
  if (MSEtool::NAor0(Data@CV_Dep[x])) stop("Data@CV_Dep is NA")

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfvec <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0vec <- -MSEtool::trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0vec <- rep(Data@vbt0[x], reps)
  }
  t0vec[!is.finite(t0vec)] <- 0
  if (all(is.nan(t0vec)))
    t0vec <- rep(0, reps)
  # hvec <- MSEtool::trlnorm(reps, Data@steep[x], Data@CV_steep[x])
  hvec <- MSEtool::sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])
  if (all(!is.finite(hvec)))
    return(list(TAC=rep(as.numeric(NA), reps), Ksamp=rep(as.numeric(NA), reps),
                dep=rep(as.numeric(NA), reps), rsamp=rep(as.numeric(NA), reps),
                MSY=rep(as.numeric(NA), reps)))

  rsamp <- getr(x, Data, Mvec, Kvec, Linfvec, t0vec, hvec, maxage = Data@MaxAge,
                r_reps = reps)
  if (is.null(dep)) dep <- MSEtool::trlnorm(reps, Data@Dep[x], Data@CV_Dep[x])
  Ct <- Data@Cat[x, ]
  Csamp <- array(rep(Ct, each = reps) * MSEtool::trlnorm(length(Ct) * reps, 1,
                                                Data@CV_Cat[x,1]), dim = c(reps, length(Ct)))
  Psamp <- array(MSEtool::trlnorm(length(Ct) * reps, 1, 0.1), dim = c(reps, length(Ct)))
  Ksamp <- rep(NA, reps)
  for (i in 1:reps) Ksamp[i] <- exp(optimize(SPSRAopt, log(c(mean(Csamp[i, ], na.rm=TRUE),
                                                             1000 * mean(Csamp[i, ], na.rm=TRUE))),
                                             dep = dep[i], r = rsamp[i], Ct = Csamp[i,],
                                             PE = Psamp[i, ])$minimum)
  MSY <- Ksamp * rsamp/4
  TAC <- MSEtool::TACfilter(Ksamp * dep * rsamp/2)
  return(list(TAC=TAC, Ksamp=Ksamp, dep=dep, rsamp=rsamp, MSY=MSY))
}


#' SPSRA Optimizer
#'
#' @param lnK internal parameter
#' @param dep internal parameter
#' @param r internal parameter
#' @param Ct internal parameter
#' @param PE internal parameter
#'
#' @export
#'
#' @describeIn SPSRA_ internal function
SPSRAopt <- function(lnK, dep, r, Ct, PE) {
  nyears <- length(Ct)
  B <- rep(NA, nyears)
  B[1] <- exp(lnK)
  OBJ <- 0
  for (y in 2:nyears) {
    if ((B[y - 1] - Ct[y - 1]) < 0)
      OBJ <- OBJ + (B[y - 1] - Ct[y - 1])^2
    B[y] <- max(0.01, B[y - 1] - Ct[y - 1])
    B[y] <- B[y] + r * B[y] * (1 - B[y]/B[1]) * PE[y]
  }
  return(OBJ + ((B[nyears]/B[1]) - dep)^2)
}


#' Surplus Production Stock Reduction Analysis
#'
#' A surplus production equivalent of DB-SRA that uses a demographically
#' derived prior for intrinsic rate of increase (McAllister method, below)
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = K D \frac{r}{2}}
#' where \eqn{K} is estimated unfished biomass, \eqn{D} is depletion, and \eqn{r}
#' is the estimated intrinsic rate of increase.
#'
#' Like all SRA methods, this MP requires a time-series of catch extending from the
#' beginning of exploitation.
#'
#' @templateVar mp SPSRA
#' @template MPtemplate
#' @template MPuses
#'
#' @author T. Carruthers
#' @references McAllister, M.K., Pikitch, E.K., and Babcock, E.A. 2001. Using
#' demographic methods to construct Bayesian priors for the intrinsic rate of
#' increase in the Schaefer model and implications for stock rebuilding. Can.
#' J. Fish. Aquat. Sci. 58: 1871-1890.
#' @family Surplus production MPs
#' @describeIn SPSRA Base version. Requires an estimate of current depletion
#' @export
#' @examples
#' SPSRA(1, MSEtool::SimulatedData, plot=TRUE)
SPSRA <- function(x, Data, reps = 100, plot=FALSE) {
  # Surplus productin stock reduction analysis T.Carruthers - basically
  # an SP version of DBSRA
  runSPSRA <- SPSRA_(x, Data, reps)
  TAC <- runSPSRA$TAC

  if (plot) SPSRA_plot(runSPSRA, Data, x)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SPSRA) <- "MP"

#' @templateVar mp SPSRA_ML
#' @template MPuses
#' @describeIn SPSRA Variant that uses a mean-length mortality estimator to obtain
#' a prior for current stock depletion.
#' @export
#' @examples
#' SPSRA_ML(1, MSEtool::SimulatedData, plot=TRUE)
SPSRA_ML <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "Data@Mort, Data@CV_Mort, Data@vbK, Data@CV_vbK, Data@vbLinf, Data@CV_vbLinf, Data@vbt0, Data@CV_vbt0, Data@CAL, Data@Cat, Data@steep"
  if (MSEtool::NAor0(Data@Mort[x])) stop("Data@Mort is NA")
  if (MSEtool::NAor0(Data@FMSY_M[x])) stop("Data@FMSY_M is NA")
  if (MSEtool::NAor0(Data@L50[x])) stop("Data@L50 is NA")
  if (MSEtool::NAor0(Data@vbK[x])) stop("Data@vbK is NA")
  if (MSEtool::NAor0(Data@CV_vbK[x])) stop("Data@CV_vbK is NA")
  if (MSEtool::NAor0(Data@vbLinf[x])) stop("Data@vbLinf is NA")
  if (MSEtool::NAor0(Data@CV_vbLinf[x])) stop("Data@CV_vbLinf is NA")
  if (MSEtool::NAor0(Data@vbt0[x])) stop("Data@vbt0 is NA")
  if (MSEtool::NAor0(Data@CV_Mort[x])) stop("Data@CV_Mort is NA")
  if (MSEtool::NAor0(Data@wla[x])) stop("Data@wla is NA")
  if (MSEtool::NAor0(Data@wlb[x])) stop("Data@wlb is NA")

  Mvec <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfvec = MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0vec <- -MSEtool::trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0vec <- rep(Data@vbt0[x], reps)
  }
  t0vec[!is.finite(t0vec)] <- 0
  # hvec <- MSEtool::trlnorm(reps, Data@steep[x], Data@CV_steep[x])
  hvec <- MSEtool::sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])
  rsamp <- getr(x, Data, Mvec, Kvec, Linfvec, t0vec, hvec, maxage = Data@MaxAge,
                r_reps = reps)
  Z <- MLne(x, Data, Linfc = Linfvec, Kc = Kvec, ML_reps = reps,  MLtype = "dep")
  if (all(is.na(Z))) {
    Rec <- new("Rec")
    Rec@TAC <- MSEtool::TACfilter(rep(NA, reps))
    return(Rec)
  }
  FM <- Z - Mvec
  FM[FM<0.05] <- 0.05

  nyears <- length(Data@Year)
  Ct1 <- mean(Data@Cat[x, 1:3], na.rm=TRUE)
  Ct2 <- mean(Data@Cat[x, (nyears - 2):nyears], na.rm=TRUE)
  dep <- rep(c(Ct1, Ct2), each = reps)/(1 - exp(-FM))
  if (reps == 1) {
    dep <- dep[2]/dep[1]
    dep <- min(dep, 0.8) # make depletion < 0.8 # sometimes estimated dep >> 1
  }
  if (reps > 1) {
    dep <- dep[, 2]/dep[, 1]
    for (rr in 1:reps) dep[rr] <- min(dep[rr], 0.8)
  }



  runSPSRA <- SPSRA_(x, Data, reps, dep=dep)
  TAC <- runSPSRA$TAC

  if (plot) SPSRA_plot(runSPSRA, Data, x)

  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(SPSRA_ML) <- "MP"



# A generic VPA (Walters and Licandeo UBC)
# VPA <- function(x, Data, reps = reps) {
#
#   # now do optimization for FMSY
#   dependencies = "Data@Mort, Data@CV_Mort, Data@vbK, Data@CV_vbK, Data@vbLinf, Data@CV_vbLinf, Data@vbt0, Data@CV_vbt0, Data@MaxAge, Data@wla, Data@CV_wla, Data@wlb, Data@CV_wlb, Data@L50, Data@CV_L50, Data@CAA, Data@steep, Data@CV_steep, Data@LFS, Data@CV_LFS, Data@LFC, Data@CV_LFC, Data@Cat"
#   CAAind <- (Data@CAA[x, , ] == 0) * array(rep(1:Data@MaxAge,
#                                                each = length(Data@CAA[x, , 1])), dim(Data@CAA[x, , ]))
#   maxage <- min(CAAind[CAAind != 0])
#   maxage <- which.min(abs(cumsum(apply(Data@CAA[x, , ], 2, sum))/sum(Data@CAA[x,
#                                                                               , ]) - 0.75))
#   CAAv <- Data@CAA[x, , 1:maxage]
#   CAAv[, maxage] <- CAAv[, maxage] + apply(Data@CAA[x, , (maxage +
#                                                             1):length(Data@CAA[x, 1, ])], 1, sum)
#
#   TAC <- Bt_K <- rep(NA, reps)
#
#   for (i in 1:reps) {
#
#     Mc <- MSEtool::trlnorm(1, Data@Mort[x], Data@CV_Mort[x])
#     # hc <- MSEtool::trlnorm(1, Data@steep[x], Data@CV_steep[x])
#     hc <- sample_steepness2(1, Data@steep[x], Data@CV_steep[x])
#     Linfc <- MSEtool::trlnorm(1, Data@vbLinf[x], Data@CV_vbLinf[x])
#     Kc <- MSEtool::trlnorm(1, Data@vbK[x], Data@CV_vbK[x])
#     t0c <- -MSEtool::trlnorm(1, -Data@vbt0[x], Data@CV_vbt0[x])
#     LFSc <- MSEtool::trlnorm(1, Data@LFS[x], Data@CV_LFS[x])
#     LFCc <- MSEtool::trlnorm(1, Data@LFC[x], Data@CV_LFC[x])
#     AMc <- MSEtool::trlnorm(1, iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x],
#                           Data@L50[x]), Data@CV_L50[x])
#     ac <- MSEtool::trlnorm(1, Data@wla[x], Data@CV_wla[x])
#     bc <- MSEtool::trlnorm(1, Data@wlb[x], Data@CV_wlb[x])
#
#     pmat <- rep(1, maxage)
#     pmat[1:ceiling(AMc)] <- 0
#     age <- 1:maxage
#     la <- Data@vbLinf[x] * (1 - exp(-Data@vbK[x] * ((age -
#                                                        Data@vbt0[x]))))
#     wa <- ac * la^bc
#
#     Cat <- Data@Cat[x, ]
#     Cat[1] <- Cat[2]  # temporary fix until effort simulation function gets sorted
#
#     CAAv[, maxage][CAAv[, maxage] == 0] <- 1
#
#     opt = optim(c(-3, -2), VPAopt, Cat = CAAv, yt = Data@Ind[x, ],
#                 S = exp(-Mc), maxage = maxage, wa = wa, pmat = pmat, method = "L-BFGS-B",
#                 lower = c(-5, -5), upper = c(5, 5))
#     out = VPAopt(opt$par, Cat = CAAv, yt = Data@Ind[x, ], S = exp(-Data@Mort[x]),
#                  maxage = maxage, wa = wa, pmat = pmat, opt = F)
#
#     fit2 <- optimize(VPAFMSY, log(c(1e-04, 3)), Mc = Mc, hc = hc, maxage = maxage,
#                      vul = out$va, Linfc = Linfc, Kc = Kc, t0c = t0c, AMc = AMc,
#                      ac = ac, bc = bc)
#     FMSY <- VPAFMSY(fit2$minimum, Mc = Mc, hc = hc, maxage = maxage,
#                     vul = out$va, Linfc = Linfc, Kc = Kc, t0c = t0c, AMc = AMc,
#                     ac = ac, bc = bc, opt = F)
#     if ((FMSY/Mc) > 3) FMSY <- 3 * Mc
#     TAC[i] <- out$bt[length(out$bt)] * FMSY
#   }
#
#   Rec <- new("Rec")
#   Rec@TAC <- MSEtool::TACfilter(TAC)
#   Rec
#
# }
# class(VPA) <- "MP"
#

#### Yield-per-recruit MPs ####

#' Internal function for YPR MPs
#'
#' @param x Iteration number
#' @param Data Object of class Data
#' @param reps Number of reps
#' @param Abun Optional numeric of length `reps` of abundance
#'
#' @export
#' @describeIn YPR_ Internal function for YPR MPs
#'
#' @keywords internal
YPR_ <- function(x, Data, reps = 100, Abun=NULL) {
  Linfc <- MSEtool::trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- MSEtool::trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0c <- -MSEtool::trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0c <- rep(Data@vbt0[x], reps)
  }
  t0c[!is.finite(t0c)] <- 0
  Mdb <- MSEtool::trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  LFS <- MSEtool::trlnorm(reps, Data@LFS[x], Data@CV_LFS[x])
  a <- Data@wla[x]
  b <- Data@wlb[x]
  runYPR <- YPRopt(Linfc, Kc, t0c, Mdb, a, b, LFS, maxage=Data@MaxAge, reps)
  if (is.null(Abun)) Abun <- MSEtool::trlnorm(reps, Data@Abun[x], Data@CV_Abun[x])
  TAC <- Abun * runYPR$F0.1
  runYPR$TAC <- TAC
  runYPR$Ac <- Abun
  return(runYPR)

}



#' @describeIn YPR_ YPR calculation function
#' @note Yield per recruit estimate of FMSY Meaghan Bryan 2013
#' @param Linfc von Bertalanffy  Linf
#' @param Kc von Bertalanffy K
#' @param t0c von Bertalanffy t0
#' @param a alpha parameter of L-W relationship
#' @param b beta parameter of L-W relationship
#' @param LFS Length at full selection
#' @param maxage Maximum age
#'
#' @export
YPRopt <- function(Linfc, Kc, t0c, Mdb, a, b, LFS, maxage, reps = 100) {

  nf <- 200
  frates <- seq(0, 3, length.out = nf)
  rat <- LFS/Linfc
  rat[rat > 0.8] <- 0.8  # need to robustify this for occasionally very high samples of LFS
  tc = log(1 - rat)/-Kc + t0c
  tc = round(tc, 0)
  tc[tc < 1] <- 1
  tc[tc > maxage] <- maxage

  vul <- array(0, dim = c(reps, maxage))
  mat <- array(0, dim = c(reps, maxage))
  lx <- array(NA, dim = c(reps, maxage))
  lxo <- array(NA, dim = c(reps, maxage))

  ypr <- array(NA, dim = c(reps, nf))
  sbpr <- array(NA, dim = c(reps, nf))
  sbpr.ratio <- array(NA, dim = c(reps, nf))
  sbpr.dif <- array(NA, dim = c(reps, nf))

  # average weight at age - follow von Bertalanffy growth
  age <- array(rep(1:maxage, each = reps), dim = c(reps, maxage))
  la <- Linfc * (1 - exp(-Kc * ((age - t0c))))
  wa <- a * la^b

  # vulnerability schedule - assumes knife-edge vulnerability, where all
  # individuals age tc to maxage are fully vulnerbale all individulas
  # less than age tc are not vulnerable
  for (i in 1:reps) {
    if (tc[i] > 0) vul[i, tc[i]:maxage] <- 1
    if (tc[i] > 1) mat[i, max(1, tc[i] - 1):maxage] <- 1
  }

  lx[, 1] <- 1
  lxo[, 1] <- 1
  for (k in 1:nf) {
    for (i in 2:maxage) {
      lx[, i] = lx[, i - 1] * exp(-(Mdb + vul[, i - 1] * frates[k]))
      lxo[, i] = lx[, i] * exp(-Mdb)
    }
    phi_vb = apply(lx * wa * vul, 1, sum)
    sbpro = apply(lxo * wa * mat, 1, sum)

    ypr[, k] = (1 - exp(-frates[k])) * phi_vb
    sbpr[, k] = apply(lx * wa * mat, 1, sum)
    sbpr.ratio[, k] = sbpr[, k]/sbpro
    sbpr.dif[, k] = abs(sbpr.ratio[, k] - 0.3)  #hard code comparison ratio
  }

  # frates[apply(ypr,1,which.max)] Fmaxypr

  # More code that derived F0.1 in 'per recruit analysis.R' (Meaghan
  # Bryan)
  slope.origin = (ypr[, 2] - ypr[, 1])/(frates[2] - frates[1])
  slope.10 = round(0.1 * slope.origin, 2)

  slope = array(NA, dim = dim(ypr))  #vector(length=length(ypr))
  slope[, 1] = slope.origin
  for (i in 3:ncol(ypr)) {
    slope[, i - 1] = round((ypr[, i] - ypr[, i - 1])/(frates[i] - frates[i - 1]), 2)
  }
  dif = abs(slope - slope.10)
  dif[is.na(dif)] <- 1e+11
  F0.1 <- frates[apply(dif, 1, which.min)]  #frates[which.min(dif)]

  return(list(F0.1=F0.1, frates=frates, ypr=ypr, dif=dif))
}




#' Yield Per Recruit analysis to get FMSY proxy F01
#'
#' A simple yield per recruit approximation to FMSY (F01) which is the position
#' of the ascending YPR curve for which dYPR/dF = 0.1(dYPR/d0)
#'
#' The TAC is calculated as:
#' \deqn{\textrm{TAC} = F_{0.1} A}
#' where \eqn{F_{0.1}} is the fishing mortality (*F*) where the slope of the yield-per-recruit
#' (YPR) curve is 10\% of the slope at the origin, and *A* is an estimate of current abundance.
#'
#' The YPR curve is calculated using an equilibrium age-structured model with life-history and
#' selectivity parameters sampled from the `Data` object.
#'
#' The variants of the YPR MP differ in the method to estimate current abundance (see Functions section below). #'
#' @templateVar mp YPR
#' @template MPtemplate
#' @template MPuses
#'
#' @note Based on the code of Meaghan Bryan
#' @author Meaghan Bryan and Tom Carruthers
#' @references Beverton and Holt. 1954.
#' @describeIn YPR Requires an external estimate of abundance.
#' @examples
#' YPR(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
YPR <- function(x, Data, reps = 100, plot=FALSE) {
  runYPR <- YPR_(x, Data, reps = reps, Abun=NULL)
  TAC <- MSEtool::TACfilter(runYPR$TAC)

  if (plot) YPR_plot(runYPR, Data, reps)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec

}
class(YPR) <- "MP"


#' @templateVar mp YPR_CC
#' @template MPuses
#' @param Fmin The minimum fishing mortality rate inferred from the catch-curve
#' analysis
#' @describeIn YPR A catch-curve analysis is used to determine recent Z which given M (Mort)
#' gives F and thus abundance = Ct/(1-exp(-F))
#' @examples
#' YPR_CC(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
YPR_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {

  # MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  MuC <- mean(Data@Cat[x, ], na.rm=TRUE)
  Cc <- MSEtool::trlnorm(reps, MuC, Data@CV_Cat[x,1])

  Mdb <- MSEtool::trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  Zdb <- CC(x, Data, reps = reps * 10)
  Fdb <- Zdb - Mdb
  ind <- (1:(reps * 10))[Fdb > Fmin][1:reps]

  Fdb <- Fdb[ind]
  Mdb <- Mdb[ind]
  SM <- sum(is.na(ind))
  if (SM > 0) {
    Mdb[is.na(ind)] <- MSEtool::trlnorm(SM, Data@Mort[x], Data@CV_Mort[x])
    Fdb[is.na(ind)] <- Fmin
  }

  Ac <- Cc/(1 - exp(-Fdb))

  runYPR <- YPR_(x, Data, reps = reps, Abun=Ac)
  TAC <- MSEtool::TACfilter(runYPR$TAC)

  if (plot) YPR_plot(runYPR, Data, reps)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(YPR_CC) <- "MP"


#' @templateVar mp YPR_ML
#' @template MPuses
#' @describeIn YPR A mean-length estimate of recent Z is used to infer current
#' abundance.
#' @examples
#' YPR_ML(1, MSEtool::SimulatedData, plot=TRUE)
#' @export
YPR_ML <- function(x, Data, reps = 100, plot=FALSE) {

  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- MSEtool::trlnorm(reps*10, MuC, Data@CV_Cat[x,1])
  Linfc <- MSEtool::trlnorm(reps*10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- MSEtool::trlnorm(reps*10, Data@vbK[x], Data@CV_vbK[x])
  Mdb <- MSEtool::trlnorm(reps*10, Data@Mort[x], Data@CV_Mort[x])
  Z <- MLne(x, Data, Linfc = Linfc, Kc = Kc, ML_reps = reps * 10, MLtype = "F")

  if (all(is.na(Z))) {
    out <- new("Rec")
    out@TAC <- rep(as.numeric(NA), reps)
    return(out)
  }

  FM <- Z - Mdb
  ind <- which(FM>0)[1:reps]
  Ac <- Cc[ind]/(1 - exp(-FM[ind]))

  runYPR <- YPR_(x, Data, reps = reps, Abun=Ac)
  TAC <- MSEtool::TACfilter(runYPR$TAC)

  if (plot) YPR_plot(runYPR, Data, reps)
  Rec <- new("Rec")
  Rec@TAC <- TAC
  Rec
}
class(YPR_ML) <- "MP"

Try the DLMtool package in your browser

Any scripts or data that you put into this service are public.

DLMtool documentation built on June 20, 2022, 5:14 p.m.