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, DLMtool::Cobia, reps=1000, plot=TRUE) # 1,000 log-normal samples with CV = 0.2
#'         
AvC <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected] [email protected]"
  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*.
#' 
#' @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. `[email protected]`
#' and `[email protected]_abun` contain values
#' 
#' @examples 
#' \dontrun{
#' BK(1, DLMtool::SimulatedData, reps=1000, plot=TRUE)
#' }
#' 
#' @author T. Carruthers.
#' @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.
#' @templateVar mp BK
#' @template MPtemplate
#' @template MPuses
#' @export 
#' 
#' 
BK <- function(x, Data, reps = 100, plot=FALSE) {
  # Beddington and Kirkwood life-history analysis
  dependencies = "[email protected], [email protected], [email protected]_vbLinf, [email protected], [email protected]_Abun, [email protected], [email protected]_vbK"
  Lc <- trlnorm(reps * 10, Data@LFC[x], 0.2)
  Linfc <- trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Ac <- trlnorm(reps * 10, Data@Abun[x], Data@CV_Abun[x])
  Kc <- 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 <- 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, DLMtool::SimulatedData, reps=1000, plot=TRUE)
#' }
#' 
#' @export 
BK_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  dependencies = "[email protected], [email protected], [email protected]_vbLinf, [email protected], [email protected]_vbK, [email protected], [email protected]"
  Lc <- trlnorm(reps, Data@LFC[x], 0.2)
  Linfc <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Mdb <- trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- 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)] <- 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 <- 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, DLMtool::SimulatedData, reps=100, plot=TRUE)
#' }
#' 
#' @export 
BK_ML <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected], [email protected], [email protected]_vbLinf, [email protected], [email protected]_vbK, [email protected], [email protected]"
  Lc <- trlnorm(reps * 10, Data@LFC[x], 0.2)
  Linfc <- trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- trlnorm(reps * 10, Data@vbK[x], Data@CV_vbK[x])
  Mdb <- 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 <- 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 <- 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, DLMtool::Cobia, plot=TRUE)
#' 

#### Constant Catch MPs ####
CC1 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, xx = 0) {
  dependencies = "[email protected], [email protected]_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 <- [email protected][x, (length([email protected]) - (yrsmth - 1)):length([email protected])]
  C_dat <- Data@Cat[x, yrfirst:yrlast]
  TAC <- (1 - xx) * 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 <- 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,], 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, DLMtool::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, DLMtool::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, DLMtool::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, DLMtool::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, DLMtool::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 `[email protected]_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, DLMtool::SimulatedData, plot=TRUE)
#' @family Constant Catch MPs
GB_CC <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected],[email protected]"
  Catrec <- Data@Cat[x, length(Data@Cat[x, ])]
  TAC <- 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 <- 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 <- trlnorm(1, Data@Mort[x], Data@CV_Mort)
    hc <- sample_steepness2(1, Data@steep[x], Data@CV_steep[x])
    Linfc <- trlnorm(1, Data@vbLinf[x], Data@CV_vbLinf[x])
    Kc <- trlnorm(1, Data@vbK[x], Data@CV_vbK[x])
    if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
      t0c <- -trlnorm(1, -Data@vbt0[x], Data@CV_vbt0[x])
    } else {
      t0c <- Data@vbt0[x]
    }
    t0c[!is.finite(t0c)] <- 0
    LFSc <- trlnorm(1, Data@LFS[x], Data@CV_LFS[x])
    LFCc <- trlnorm(1, Data@LFC[x], Data@CV_LFC[x])
    AMc <- trlnorm(1, iVB(Data@vbt0[x], Data@vbK[x], Data@vbLinf[x], Data@L50[x]), Data@CV_L50[x])
    ac <- trlnorm(1, Data@wla[x], Data@CV_wla[x])
    bc <- 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, DLMtool::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 = "[email protected], [email protected]_Mort, [email protected], [email protected]_vbK, [email protected], [email protected]_vbLinf, [email protected], [email protected]_vbt0, [email protected], [email protected], [email protected]_wla, [email protected], [email protected]_wlb, [email protected], [email protected]_L50, [email protected], [email protected], [email protected]_steep, [email protected], [email protected]_LFS, [email protected], [email protected]_LFC, [email protected]"
 
  runCompSRA <- CompSRA_(x, Data, reps)
  TAC <- 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, DLMtool::SimulatedData, plot=TRUE)
#' 
#' @export CompSRA4010
CompSRA4010 <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected], [email protected]_Mort, [email protected], [email protected]_vbK, [email protected], [email protected]_vbLinf, [email protected], [email protected]_vbt0, [email protected], [email protected], [email protected]_wla, [email protected], [email protected]_wlb, [email protected], [email protected]_L50, [email protected], [email protected], [email protected]_steep, [email protected], [email protected]_LFS, [email protected], [email protected]_LFC, [email protected]"
  
  runCompSRA <- CompSRA_(x, Data, reps)
  
  TAC <- 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 (NAor0(Data@BMSY_B0[x])) stop('[email protected]_B0 is NA')
  if (NAor0(Data@CV_BMSY_B0[x])) stop("[email protected]_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 <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])  # CV of 0.5 as in MacCall 2009
  FMSY_M <- 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 <- 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("[email protected]_B0 or [email protected]_BMSY_B0 do not contain values")
    return(list(dcac=rep(NA, reps), Bt_K=Bt_K))
  }
  BMSY_K <- rbeta(reps, alphaconv(Data@BMSY_B0[x], Data@BMSY_B0[x] * Data@CV_BMSY_B0[x]), 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, DLMtool::Atlantic_mackerel, plot=TRUE) 
#' 
#' @export 
DCACs <- function(x, Data, reps = 100, plot=FALSE) {
  rundcac <- DCAC_(x, Data, reps)
  TAC <- 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, DLMtool::Atlantic_mackerel, plot=TRUE) 
#' 
#' @export 
DCAC <- function(x, Data, reps = 100, plot=FALSE) {
  rundcac <- DCAC_(x, Data, reps, updateD=TRUE)
  TAC <- 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, DLMtool::Atlantic_mackerel, plot=TRUE) 
#' 
#' @export 
DCAC_40 <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected], [email protected], [email protected], [email protected]_Mort, [email protected]_M, [email protected]_FMSY_M, [email protected]_B0, [email protected]_BMSY_B0"
  
  rundcac <- DCAC_(x, Data, reps, Bt_K=0.4)
  TAC <- 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 <- DLMtool::Atlantic_mackerel
#' [email protected] <- 2005
#' DCAC4010(1, Data, plot=TRUE) 
#' 
#' @export 
DCAC4010 <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected], [email protected], [email protected], [email protected]_Mort, [email protected]_M, [email protected]_FMSY_M, [email protected], [email protected]_Dt, [email protected]_B0, [email protected]_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 <- 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, DLMtool::SimulatedData, plot=TRUE) 
#' 
#' @export 
DCAC_ML <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected], [email protected], [email protected], [email protected]_Mort, [email protected]_M, [email protected]_FMSY_M, [email protected]_B0, [email protected]_BMSY_B0, [email protected], [email protected], [email protected], [email protected]_vbLinf, [email protected], [email protected]_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 <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])  
  FMSY_M <- trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x])
  Linfc <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- 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 <- 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 <- DLMtool::Atlantic_mackerel
#' [email protected] <- 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 = "[email protected], [email protected], [email protected], [email protected]_Mort, [email protected]_M, [email protected]_FMSY_M, [email protected], [email protected]_Dt, [email protected]_B0, [email protected]_BMSY_B0"
  
  rundcac <- DCAC_(x, Data, reps, updateD = TRUE)
  TAC <- rundcac$dcac * rundcac$Bt_K/rundcac$BMSY_K
  
  TAC <- 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 <- DLMtool::Atlantic_mackerel
#' [email protected] <- 2005
#' HDAAC(1, Data, plot=TRUE) 
#' 
#' @export 
HDAAC <- function(x, Data, reps = 100, plot=FALSE) {
  dependencies = "[email protected], [email protected], [email protected], [email protected]_Mort, [email protected]_M, [email protected]_FMSY_M, [email protected], [email protected]_Dt, [email protected]_B0, [email protected]_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 <- 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 (NAor0(Data@CV_BMSY_B0[x])) stop("[email protected]_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, alphaconv(depo, min(depo * Data@CV_Dep[x], (1 - depo) * Data@CV_Dep[x])), 
                  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 <- 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 <- trlnorm(1, Data@FMSY_M[x], Data@CV_FMSY_M[x])
    BMSY_K <- rbeta(100, alphaconv(Data@BMSY_B0[x], Data@CV_BMSY_B0[x] *  Data@BMSY_B0[x]), 
                    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 `[email protected]_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, DLMtool::ourReefFish, plot=TRUE)
#' 
#' @export 
DBSRA <- function(x, Data, reps = 100, plot=FALSE) {
  runDBSRA <- DBSRA_(x,Data, reps)
  TAC <- 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, DLMtool::ourReefFish, plot=TRUE)
#' @export 
DBSRA_40 <- function(x, Data, reps = 100, plot=FALSE) {
  runDBSRA <- DBSRA_(x,Data, reps, depo=0.4)
  TAC <- 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, DLMtool::ourReefFish, plot=TRUE)
#' @export 
#' 
DBSRA4010 <- function(x, Data, reps = 100, plot=FALSE) {
  runDBSRA <- DBSRA_(x,Data, reps, hcr=c(0.4, 0.1))
  TAC <- 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 = "[email protected], [email protected], [email protected]_Mort, [email protected]_M, [email protected]_FMSY_M, [email protected]_B0, [email protected]_BMSY_B0, [email protected], [email protected], [email protected], [email protected]"
#   C_hist <- [email protected][x, ]
#   TAC <- rep(NA, reps)
#   DBSRAcount <- 1
#   maxIts <- 200
#   nIts <- 0
#   if (is.na([email protected][x]) | is.na([email protected]_Dep[x])) return(NA)
#   while (DBSRAcount < (reps + 1) & nIts < maxIts) {
#     Linfc <- trlnorm(1, [email protected][x], [email protected]_vbLinf[x])
#     Kc <- trlnorm(1, [email protected][x], [email protected]_vbK[x])
#     Mdb <- trlnorm(100, [email protected][x], [email protected]_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")
#       [email protected] <- TACfilter(rep(NA, reps))
#       return(Rec)
#     } 
#     FM <- Z - Mdb
#     FM[FM < 0] <- 0.01
#     nyears <- length([email protected])
#     Ct1 <- mean([email protected][x, 1:3])
#     Ct2 <- mean([email protected][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 <- trlnorm(1, [email protected]_M[x], [email protected]_FMSY_M[x])
#     BMSY_K <- rbeta(100, alphaconv([email protected]_B0[x], [email protected]_BMSY_B0[x] * 
#                                      [email protected]_B0[x]), betaconv([email protected]_B0[x], [email protected]_BMSY_B0[x] * 
#                                                                   [email protected]_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([email protected][x], [email protected][x], [email protected][x],  [email protected][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")
#   [email protected] <- 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 <- 1: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(alphaconv(UMSYpriorpar[1], prod(UMSYpriorpar)), 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=DLMtool::Atlantic_mackerel, plot=TRUE)
#' @family Delay-Difference MPs
#' 
#' @export 
DD <- function(x, Data, reps = 100, plot=FALSE) {
  runDD <- DD_(x, Data, reps)
  TAC <- 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=DLMtool::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 <- 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 (`[email protected]`) 
#' and an estimate of abundance (`[email protected]`) 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=DLMtool::Atlantic_mackerel, plot=TRUE)
#' @export 
DynF <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2) {
  
  dependencies = "[email protected], [email protected], [email protected], [email protected], [email protected], [email protected]_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 <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x]) * 
    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 <- 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 (`[email protected]`) 
#' and estimate of current abundance (`[email protected]`), 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=DLMtool::Atlantic_mackerel, plot=TRUE)
#' @export 
Fadapt <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 7, gg = 1) {
  
  dependencies = "[email protected], [email protected], [email protected], [email protected], [email protected], [email protected]_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 <- 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 (NAor0(Data@Mort[x])) stop("[email protected] is NA")
  if (NAor0(Data@FMSY_M[x])) stop("[email protected]_M is NA")
  if (NAor0(Data@L50[x])) stop("[email protected] is NA")
  if (NAor0(Data@vbK[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_vbK[x])) stop("[email protected]_vbK is NA")
  if (NAor0(Data@vbLinf[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_vbLinf[x])) stop("[email protected]_vbLinf is NA")
  if (NAor0(Data@vbt0[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_Mort[x])) stop("[email protected]_Mort is NA")
  if (NAor0(Data@wla[x])) stop("[email protected] is NA")
  if (NAor0(Data@wlb[x])) stop("[email protected] is NA")

  Mvec <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kc <- trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfc = 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 <- -trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
    } else {
      t0c <- trlnorm(reps, Data@vbt0[x], Data@CV_vbt0[x])
    }
  } else {
    t0c <- rep(Data@vbt0[x], reps)
  }
  t0c[!is.finite(t0c)] <- 0
  # hvec <- trlnorm(reps, [email protected][x], [email protected]_steep[x])
  hvec <- sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])
  
  if (is.null(Ac)) Ac <- 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 `[email protected]`)
#' @examples 
#' Fdem(1, DLMtool::SimulatedData, plot=TRUE)
#' @export 
Fdem <- function(x, Data, reps = 100, plot=FALSE) {
  runFdem <- Fdem_(x, Data, reps)
  TAC <- 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, DLMtool::SimulatedData, plot=TRUE)
#' @export 
Fdem_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  
  Mvec <- trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- 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 <- 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, DLMtool::SimulatedData, plot=TRUE)
#' @export 
Fdem_ML <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  Mvec <- trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])
  Kc <- trlnorm(reps * 10, Data@vbK[x], Data@CV_vbK[x])
  Linfc = trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  MuC <- Data@Cat[x, length(Data@Cat[x, ])]
  Cc <- 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 <- 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 <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x]) * trlnorm(reps, Data@FMSY_M[x], Data@CV_FMSY_M[x]) # estimate of Fmsy
  if (is.null(Abun)) Abun <- 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 `[email protected]`)
#' @family Fmsy/M methods
#' @examples 
#' Fratio(1, DLMtool::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 <- 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, DLMtool::Atlantic_mackerel, plot=TRUE)
#' @export 
Fratio4010 <- function(x, Data, reps = 100, plot=FALSE) {
  runFrat <- Fratio_(x, Data, reps)
  TAC <- TACfilter(runFrat$TAC)

  # 40-10 rule
  # Bt_K <- trlnorm(reps, [email protected][x], [email protected]_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, alphaconv(depo, min(depo * Data@CV_Dep[x], 
                                                (1 - depo) * Data@CV_Dep[x])), 
                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 <- 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, alphaconv(depo, min(depo * Data@CV_Dep[x], 
                                                (1 - depo) * Data@CV_Dep[x])), 
                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, DLMtool::SimulatedData, plot=TRUE)
#' @export
Fratio_CC <- function(x, Data, reps = 100, plot=FALSE, Fmin = 0.005) {
  # estimate abundance from average catch and F
  # MuC <- [email protected][x, length([email protected][x, ])]
  MuC <- mean(Data@Cat[x, ], na.rm=TRUE)
  Cc <- trlnorm(reps, MuC, Data@CV_Cat[x,1])
  Mdb <- 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)] <- 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 <- 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, DLMtool::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 <- trlnorm(reps * 10, MuC, Data@CV_Cat[x,1])
  Mdb <- trlnorm(reps * 10, Data@Mort[x], Data@CV_Mort[x])  # CV of 0.5 as in MacCall 2009
  Linfc <- trlnorm(reps * 10, Data@vbLinf[x], Data@CV_vbLinf[x])
  Kc <- 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 <- 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 (`[email protected]`) 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, DLMtool::SimulatedData, plot=TRUE)
#' @export 
GB_slope <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 5, lambda = 1) {
  dependencies = "[email protected], [email protected], [email protected]_Cat, [email protected]"
  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 <-  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 <- 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 (`[email protected]`),
#' *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 (`[email protected]`), 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 `[email protected]_Cref`, and \eqn{I_\textrm{target}} is assumed to be the index at MSY subject 
#' to observation error (`[email protected]_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, DLMtool::SimulatedData, plot=TRUE)
#' @family Index methods
GB_target <- function(x, Data, reps = 100, plot=FALSE, w = 0.5) {
  dependencies = "[email protected], [email protected], [email protected], [email protected]"
  Catrec <- Data@Cat[x, length(Data@Cat[x, ])]
  TACtarg <- trlnorm(reps, Data@Cref[x], Data@CV_Cref)
  Itarg <- 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 <- 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=DLMtool::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 = "[email protected], [email protected], [email protected], [email protected]"
  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 <- 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([email protected][x,],each=reps),c(reps,length([email protected][x,]))),TAC)
  # Warr<-([email protected][x]*exp([email protected][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=DLMtool::Atlantic_mackerel, plot=TRUE)
#' 
ICI <- function(x, Data, reps=100, plot=FALSE) {
  dependencies = "[email protected], [email protected]_Ind, [email protected], [email protected]_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 <- 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 <- trlnorm(reps, Cat, Data@CV_Cat[x,1])
  
  TAC <- 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=DLMtool::Atlantic_mackerel, plot=TRUE)
ICI2 <- function(x, Data, reps=100, plot=FALSE) {
  dependencies = "[email protected], [email protected]_Ind, [email protected], [email protected]_Cat"
  
  Index <- Data@Ind[x, ]
  Years <- Data@Year[!is.na(Index)]
  Index <- Index[!is.na(Index)]
  nI <- length(Index)
  Ind.samp <- 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 <- trlnorm(reps, Cat, Data@CV_Cat[x,1])
  
  TAC <- alpha * Cc
  TAC <- 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, DLMtool::Atlantic_mackerel, plot=TRUE)
#' @export
Iratio <- function(x, Data, reps=100, plot=FALSE, yrs = c(2, 5)) {
  dependencies = "[email protected], [email protected]_Ind, [email protected], [email protected]_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 <- trlnorm(reps * length(ind.num), Data@Ind[x, ind.num], Data@CV_Ind[x,1])
    I.den <- trlnorm(reps * length(ind.den), Data@Ind[x, ind.den], Data@CV_Ind[x,1])
  }
  
  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 <- trlnorm(reps, Cat, Data@CV_Cat[x,1])
  TAC <- alpha * Cc
  TAC <- 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) * 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, DLMtool::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 <- 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, DLMtool::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, DLMtool::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, DLMtool::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 = "[email protected], [email protected], [email protected]_Ind, [email protected]"
  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 * trlnorm(reps, 1, Data@CV_Ind[x,1])
  TAC <- 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 (`[email protected]`) is assumed to be a proxy for MSY. In the
#' MSE `Iref` is the index at MSY subject to observation error (`[email protected]`). 
#' 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, DLMtool::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, DLMtool::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) * 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 <- 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)], [email protected][x, (length([email protected][x,])-yrs):length([email protected][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, DLMtool::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, DLMtool::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, DLMtool::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, DLMtool::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 (`[email protected]`).
#'
#' 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=DLMtool::SimulatedData, plot=TRUE)
ITM <- function(x, Data, reps = 100, plot=FALSE) {
  
  dependencies = "[email protected], [email protected], [email protected]_Ind, [email protected], [email protected]"
  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 * trlnorm(reps, 1, Data@CV_Ind[x,1])
  TAC <- 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) * 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 <- 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=DLMtool::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=DLMtool::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=DLMtool::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=DLMtool::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=DLMtool::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=DLMtool::SimulatedData, plot=TRUE)
#' 
Lratio_BHI <- function(x, Data, reps=100, plot=FALSE, yrsmth = 3) {
  dependencies = "[email protected]_Linf, [email protected]_vbLinf, [email protected], [email protected]_Cat, [email protected], [email protected]_bins,
  [email protected], [email protected]_LFS"
  Linfc <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Lc <- 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 <- 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 <- 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=DLMtool::SimulatedData, plot=TRUE)
#' 
Lratio_BHI2 <- function(x, Data, reps=100, plot=FALSE, yrsmth = 3) {
  
  dependencies = "[email protected]_Linf, [email protected]_vbLinf, [email protected], [email protected]_Cat, [email protected], [email protected]_Mort,
  [email protected]_K, [email protected]_vbK, [email protected]_M, [email protected]_FMSY_M, [email protected], [email protected]_bins,
  [email protected], [email protected]_LFS"
  Linfc <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Lc <- trlnorm(reps, Data@LFS[x], Data@CV_LFS[x])
  
  Mvec <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  gamma <- 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 <- 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 <- 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=DLMtool::SimulatedData, plot=TRUE)
#' 
Lratio_BHI3 <- function(x, Data, reps=100, plot=FALSE, yrsmth = 3) {
  
  dependencies = "[email protected]_Linf, [email protected]_vbLinf, [email protected], [email protected]_Cat, [email protected], [email protected]_Mort,
  [email protected]_K, [email protected]_vbK, [email protected]_M, [email protected]_FMSY_M, [email protected], [email protected]_bins,
  [email protected], [email protected]_LFS"
  Linfc <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  Lc <- trlnorm(reps, Data@LFS[x], Data@CV_LFS[x])
  
  Mvec <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  gamma <- 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 <- 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 <- 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) * 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 <- 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=DLMtool::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=DLMtool::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=DLMtool::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=DLMtool::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, DLMtool::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 = "[email protected], [email protected]_Dep, [email protected]"
  if (all(is.na(Data@Cat[x, ]))) stop("No catch data")
  if (NAor0(Data@CV_Cat[x,1])) stop("[email protected]_Cat is NA")
  if (NAor0(Data@Dep[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_Dep[x])) stop("[email protected]_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, alphaconv(depo, min(depo * Data@CV_Dep[x], 
                                                (1 - depo) * Data@CV_Dep[x])), 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 <- 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, DLMtool::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 = "[email protected], [email protected]_Dep, [email protected]"
  if (all(is.na(Data@Cat[x, ]))) stop("No catch data")
  if (NAor0(Data@CV_Cat[x,1])) stop("[email protected]_Cat is NA")
  if (NAor0(Data@CV_Dep[x])) stop("[email protected]_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, alphaconv(depo, min(depo * Data@CV_Dep[x], 
                                                (1 - depo) * Data@CV_Dep[x])), 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 <- 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=DLMtool::Atlantic_mackerel, plot=TRUE)
Rcontrol <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2, glim = c(0.5, 2)) {
  dependencies = "[email protected], [email protected]_Mort, [email protected], [email protected]_vbK, [email protected], [email protected]_vbLinf, [email protected], [email protected]_vbt0 [email protected], [email protected]_steep, [email protected], [email protected], [email protected]_Dep, [email protected], [email protected]"
  if (NAor0(Data@Mort[x])) stop("[email protected] is NA")
  if (NAor0(Data@FMSY_M[x])) stop("[email protected]_M is NA")
  if (NAor0(Data@L50[x])) stop("[email protected] is NA")
  if (NAor0(Data@vbK[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_vbK[x])) stop("[email protected]_vbK is NA")
  if (NAor0(Data@vbLinf[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_vbLinf[x])) stop("[email protected]_vbLinf is NA")
  if (NAor0(Data@vbt0[x])) stop("[email protected] is NA")
  if (NAor0(Data@wla[x])) stop("[email protected] is NA")
  if (NAor0(Data@wlb[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_Mort[x])) stop("[email protected]_Mort is NA")
  
  Mvec <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfvec <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0vec <- -trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0vec <- rep(Data@vbt0[x], reps)
  }
  t0vec[!is.finite(t0vec)] <- 0
  
  # hvec <- trlnorm(reps, [email protected][x], [email protected]_steep[x])
  hvec <- 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, alphaconv(depo, min(depo * Data@CV_Dep[x], 
                                           (1 - depo) * Data@CV_Dep[x])), 
                  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 <- 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([email protected][x,],each=reps),c(reps,length([email protected][x,]))),TAC)
  # Warr<-([email protected][x]*exp([email protected][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=DLMtool::Atlantic_mackerel, plot=TRUE)
Rcontrol2 <- function(x, Data, reps = 100, plot=FALSE, yrsmth = 10, gg = 2, glim = c(0.5, 2)) {
  dependencies = "[email protected], [email protected]_Mort, [email protected], [email protected]_vbK, [email protected], [email protected]_vbLinf, [email protected], [email protected]_vbt0, [email protected], [email protected]_steep, [email protected], [email protected], [email protected]_Dep, [email protected], [email protected]"
  if (NAor0(Data@Mort[x])) stop("[email protected] is NA")
  if (NAor0(Data@FMSY_M[x])) stop("[email protected]_M is NA")
  if (NAor0(Data@L50[x])) stop("[email protected]50 is NA")
  if (NAor0(Data@vbK[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_vbK[x])) stop("[email protected]_vbK is NA")
  if (NAor0(Data@vbLinf[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_vbLinf[x])) stop("[email protected]_vbLinf is NA")
  if (NAor0(Data@vbt0[x])) stop("[email protected] is NA")
  if (NAor0(Data@CV_Mort[x])) stop("[email protected]_Mort is NA")
  if (NAor0(Data@wla[x])) stop("[email protected] is NA")
  if (NAor0(Data@wlb[x])) stop("[email protected] is NA")
  
  Mvec <- trlnorm(reps, Data@Mort[x], Data@CV_Mort[x])
  Kvec <- trlnorm(reps, Data@vbK[x], Data@CV_vbK[x])
  Linfvec <- trlnorm(reps, Data@vbLinf[x], Data@CV_vbLinf[x])
  if (Data@vbt0[x] != 0 & Data@CV_vbt0[x] != tiny) {
    t0vec <- -trlnorm(reps, -Data@vbt0[x], Data@CV_vbt0[x])
  } else {
    t0vec <- rep(Data@vbt0[x], reps)
  }
  t0vec[!is.finite(t0vec)] <- 0
  # hvec <- trlnorm(reps, [email protected][x], [email protected]_steep[x])
  hvec <- sample_steepness2(reps, Data@steep[x], Data@CV_steep[x])
  r