R/MPs_Output.R In DLMtool: Data-Limited Methods Toolkit

## 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,

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

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 forICI 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) {
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) {
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)

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) {
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)

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