R/MP_3b_HCRs.R
In flr/FLBEIA: Bio-Economic Impact Assessment of Management Strategies using FLR
Defines functions
annualTAC
Documented in
annualTAC
#' Harvest Control Rules (HCRs)
#'
#' There are several HCRs available in FLBEIA which are used within the main function FLBEIA to generate the management advice in each step.
#' But they can also be used independently.
#'
#' @param stocks And FLStocks object.
#' @inheritParams FLBEIA
#' @param year The position of the assessment year in the stocks and advice objects.
#' @param stknm The name of the stock for which advice is being generated.
#' @param ... Any extra arguments needed for specific HCRs.
#'
#' @return The advice input object updated with the management advice (TAC) generated by the HCR.
#'
#' @details There are two types of HCRs model-free HCRs and model-Based HCRs.
#' Model-free HCRs use abundance indices to generate the advice and hence it use FLIndices object as input data.
#' Model-based HCRs use estimates of stock abundance and stock exploitation level to generate the advice.
#'\itemize{
#' \item Model-Free HCRs: annexIVHCR, ghlHCR, little2011HCR, pidHCR, pidHCRItarg, IcesCat3HCR and IcesCat3HCR_bsafe_hrcap.
#' \item Model-Based HCRs: aneHCRE, annualTAC, CFPMSYHCR, F2CatchHCR, FroeseHCR, IcesHCR, MAPHCR, neaMAC_ltmp.
#'
#'}
#' Rules' description:
#'\itemize{
#' \item{aneHCRE:} {The HCR used in the bay of biscay anchovy long term management plan.}
#' \item{aneHCRs:} {The HCRs (escapement biomass) tested for the bay of biscay anchovy with different calendars.
#' For details see Sanchez et al. 2019. MEPS 617-618: 245-263.}
#' \item{annexIVHCR:} {The HCR used by EC and ICES to generate the TAC advice for data poor stocks.}
#' \item{annualTAC:} {A HCR that generates annual TAC advice. The HCR provides the whole flexibility of fwd.}
#' \item{CFPMSYHCR:} {HCR adapting the MAP HCR to allow flexibility in the year Fmsy is achieved.
#' The user can specify the year in which you aim to reach Fmsy, with a linear transition between
#' Fsq to Fmsy in the intervening years}
#' \item{F2CatchHCR:} {This HCR transforms the fishing mortality advice given as input data to catch advice without any other restriction.}
#' \item{FroeseHCR:} {The HCR defined in the paper by Froese, Branch et al. in Fish and Fisheries 2011.}
#' \item{ghlHCR:} {The model-free HCR used in the management of greenland-halibut}
#' \item{IcesHCR:} {The HCR used by ICES to generate TAC advice in the MSY framework.}
#' \item{IcesCat3HCR:} {HCR that implements the ICES HCR for Category 3 DLS stocks
#' (see ICES CM 2012/ACOM 68: Category 3 - Method 3.2)}
#' \item{IcesCat3HCR_bsafe_hrcap:} {Alternative approach for IcesCat3HCR with a biomass safeguard and harvest rate caps
#' (from ICES WKDLSSLS2019)}
#' \item{little2011HCR:} {The HCR defined in the paper by Little et al. in ICES Journal of Marine Science 2011.}
#' \item{MAPHRC:} {The HCR proposed by the EC in the evaluation on multi-annual management plans in 2015.}
#' \item{MultiStockHRC:} {A HCR that produces TAC advice for several stocks simultaneously. It uses a fishing mortality target and an upper bound to conciliate the TAC advices. In the case of stocks without exploitation rate estimates it uses the catch.}
#' \item{neaMAC_ltmp:} {The HCR used in the north-east atlantic mackerel long term management plan. It is a particular case of the IcesHCR.}
#' \item{pidHCR, pidHCRItarg:} {The HCRs defined in the paper by Pomaerede et al. in Aquatic Living Resources 2010.}
#' }
#' The HCRs are documented in detail in the manual of the library.
#'
#' @examples
#'\dontrun{
#' library(FLBEIA)
#' library(FLAssess) # required to use the IcesHCR. Not available for win64
#' library(FLash) # required to use the IcesHCR. Not available for win64
#' library(ggplot2)
#'
#' # Load the data to run FLBEIA in a one stock one fleet example using the HCR used by ICES
#' # in the MSY framework.
#' data(one)
#'
#' oneAdv$TAC[,ac(2009:2025)] <- NA # Put NA-s in the projection years to check how the
#' # function fills the advice object.
#'
#' res <- IcesHCR(oneSt, oneAdv, oneAdvC, 19, 'stk1')
#' # The value printed in the screen is the fishing mortality used in the advice.
#'
#' res$TAC[,'2009'] # The resulting management advice.
#' }
#-------------------------------------------------------------------------------
# HCRs
# - annualTAC.
#
# Dorleta GarcYYYa
# Created: 20/12/2010 13:26:13
# Changed: 20/12/2010 13:26:18
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# annualTAC(stocks, covars, advice, hcr.ctrl)
# year = Assessment year (POSITION).
#
# The targets and constraints will differ iteration by iteration, thus 'fwd'
# must be applied iter by iter.
#-------------------------------------------------------------------------------
# @export
annualTAC <- function(stocks, advice, advice.ctrl, year, stknm,...){
# project the stock 3 years, (current year, TAC year, TAC year + 1 for ssb or biomass constraints).
nyears <- ifelse(is.null(advice.ctrl[[stknm]][['nyears']]), 3, advice.ctrl[[stknm]][['nyears']])
wts.nyears <- ifelse(is.null(advice.ctrl[[stknm]][['wts.nyears']]), 3, advice.ctrl[[stknm]][['wts.nyears']])
fbar.nyears <- ifelse(is.null(advice.ctrl[[stknm]][['fbar.nyears']]), 3, advice.ctrl[[stknm]][['fbar.nyears']])
f.rescale <- ifelse(is.null(advice.ctrl[[stknm]][['f.rescale']]), TRUE, advice.ctrl[[stknm]][['f.rescale']])
# disc.nyears <- ifelse(is.null(advice.ctrl[[stknm]][['disc.nyears']]), wts.nyears, advice.ctrl[[stknm]][['disc.nyears']])
stk <- stocks[[stknm]]
stk@harvest[stk@harvest < 0] <- 0.00001
# if(dim(stk@m)[1] == 1){
# for(sl in c("catch.n","catch.wt","discards.n","discards.wt","landings.n",
# "landings.wt","stock.n","stock.wt","m","mat","harvest","harvest.spwn", "m.spwn"))
# {
# dimnames(slot(stk,sl))[[1]] <- 1
# }
# stk@range[6:7] <- 1
# }
#
# stk <- stf(stk, nyears = nyears, wts.nyears = wts.nyears, fbar.nyears = fbar.nyears, f.rescale = TRUE) #, disc.nyrs = disc.nyears)
#
ageStruct <- ifelse(dim(stk@m)[1] > 1, TRUE, FALSE)
if(ageStruct == TRUE)
stk <- stf(stk, nyears = nyears, wts.nyears = wts.nyears, fbar.nyears = fbar.nyears, f.rescale = f.rescale) #, disc.nyrs = disc.nyears)
else
stk <- stfBD(stk, nyears = nyears, wts.nyears = wts.nyears, fbar.nyears = fbar.nyears)
fwd.ctrl <- advice.ctrl[[stknm]]$fwd.ctrl
iter <- dim(stk@m)[6]
yrsnames <- dimnames(stk@m)[[2]]
yrsnumbs <- as.numeric(yrsnames)
assyrname <- yrsnames[year]
assyrnumb <- yrsnumbs[year]
# Refresh the years in fwd!!
fwd.ctrl@target$year <- fwd.ctrl@target$year + assyrnumb
fwd.ctrl@target$rel.year <- fwd.ctrl@target$rel.year + assyrnumb
TACvar <- FALSE
for(i in 1:iter){
stki <- iter(stk, i)
# if in <year 0> quantity = catch => set TAC in <year 0> in val
if(TACvar == TRUE | any(fwd.ctrl@target[fwd.ctrl@target$rel.year == assyrnumb,'quantity'] == 'catch')){
k <- which(fwd.ctrl@target$rel.year == assyrnumb & fwd.ctrl@target$quantity == 'catch')
fwd.ctrl@target[k,c('min','val', 'max')] <- fwd.ctrl@target[k,c('min','val', 'max')]*advice$TAC[stknm,year,,,,i]
fwd.ctrl@target[k,'rel.year'] <- NA
fwd.ctrl@trgtArray[k,c('min','val', 'max'),] <- fwd.ctrl@trgtArray[k,c('min','val', 'max'),]*c(advice$TAC[stknm,year,,,,i])
TACvar <- TRUE
}
# if(stknm == 'CMON') browser()
if(dim(stki@m)[1] > 1){
# First estimate/extract the SR model and params.
sr.pars <- advice.ctrl[[stknm]]$sr$params # sr parameters if specified.
sr.model <- advice.ctrl[[stknm]]$sr$model # sr model, mandatory.
if(is.null(sr.pars)){ # if params missing => estimate the parameters using the specified years.
if(is.null(advice.ctrl[[stknm]]$sr$years)) which(round(quantSums(stocks[[stknm]]@stock.n))!=0)[1]:(year-1)# yr0 missing => use all data years.
else{
y.rm <- as.numeric(advice.ctrl[[stknm]]$sr$years['y.rm'])
nyrs <- as.numeric(advice.ctrl[[stknm]]$sr$years['num.years'])
sr.yrs <- yrsnames[(year-y.rm-nyrs + 1):(year-y.rm)]
}
rec <- stki@stock.n[1,sr.yrs]
ssb <- ssb(stki)[,sr.yrs]
# if rec.age != 0 adjust rec and ssb.
rec.age <- as.numeric(dimnames(rec)[[1]])
if(rec.age != 0){
rec <- rec[, -(1:rec.age),]
ssb <- ssb[, 1:(dim(ssb)[2] - rec.age),]
}
if(sr.model != 'geomean') sr.pars <- try(params(fmle(FLSR(rec = rec, ssb = ssb, model = sr.model))), silent = TRUE)
if(class(sr.pars) == 'try-error' | sr.model == 'geomean'){
sr.model <- 'geomean'
sr.pars <- c(prod(c(rec))^(1/length(c(rec))))
sr.pars <- FLPar(a = ifelse(is.na(sr.pars), 0, sr.pars))
}
sr1 <- sr.pars
} else { # sr.pars not null
if(i == 1){
sr1 <- iter(sr.pars,i)
}
sr1[] <- iter(sr.pars,i)[]
}
stki <- FLash::fwd(stki, ctrl = fwd.ctrl, sr = list(model =sr.model, params = sr1))
} else {
# Extract the years to calculate the mean historical growth of the stock
if(is.null(advice.ctrl[[stknm]]$growth.years)) growth.years <- max(1,(year - 11)):(year-1)
else{
y.rm <- as.numeric(advice.ctrl[[stknm]]$growth.years['y.rm'])
nyrs <- as.numeric(advice.ctrl[[stknm]]$growth.years['num.years'])
growth.years <- yrsnames[(year-y.rm-nyrs + 1):(year-y.rm)]
}
stki <- fwdBD(stki, fwd.ctrl, growth.years)
}
Cadv <- ifelse(advice.ctrl[[stknm]][['AdvCatch']][year+1] == TRUE, 'catch', 'landings')
yy <- ifelse(slot(stki, Cadv)[,year+1] == 0, 1e-6, slot(stki, Cadv)[,year+1])
advice[['TAC']][stknm,year+1,,,,i] <- yy
# cat('---------------- HCR------------------------\n')
# cat(c(fbar(stki)[,(year-1):year]), '\n')
# cat('-------------------------------------------\n')
# save(stki, file = 'stki.RData')
}
return(advice)
}
flr/FLBEIA documentation built on July 14, 2024, 11:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions annualTAC
Documented in annualTAC
#' Harvest Control Rules (HCRs)
#'
#' There are several HCRs available in FLBEIA which are used within the main function FLBEIA to generate the management advice in each step.
#' But they can also be used independently.
#'
#' @param stocks And FLStocks object.
#' @inheritParams FLBEIA
#' @param year The position of the assessment year in the stocks and advice objects.
#' @param stknm The name of the stock for which advice is being generated.
#' @param ... Any extra arguments needed for specific HCRs.
#'
#' @return The advice input object updated with the management advice (TAC) generated by the HCR.
#'
#' @details There are two types of HCRs model-free HCRs and model-Based HCRs.
#' Model-free HCRs use abundance indices to generate the advice and hence it use FLIndices object as input data.
#' Model-based HCRs use estimates of stock abundance and stock exploitation level to generate the advice.
#'\itemize{
#' \item Model-Free HCRs: annexIVHCR, ghlHCR, little2011HCR, pidHCR, pidHCRItarg, IcesCat3HCR and IcesCat3HCR_bsafe_hrcap.
#' \item Model-Based HCRs: aneHCRE, annualTAC, CFPMSYHCR, F2CatchHCR, FroeseHCR, IcesHCR, MAPHCR, neaMAC_ltmp.
#'
#'}
#' Rules' description:
#'\itemize{
#' \item{aneHCRE:} {The HCR used in the bay of biscay anchovy long term management plan.}
#' \item{aneHCRs:} {The HCRs (escapement biomass) tested for the bay of biscay anchovy with different calendars.
#' For details see Sanchez et al. 2019. MEPS 617-618: 245-263.}
#' \item{annexIVHCR:} {The HCR used by EC and ICES to generate the TAC advice for data poor stocks.}
#' \item{annualTAC:} {A HCR that generates annual TAC advice. The HCR provides the whole flexibility of fwd.}
#' \item{CFPMSYHCR:} {HCR adapting the MAP HCR to allow flexibility in the year Fmsy is achieved.
#' The user can specify the year in which you aim to reach Fmsy, with a linear transition between
#' Fsq to Fmsy in the intervening years}
#' \item{F2CatchHCR:} {This HCR transforms the fishing mortality advice given as input data to catch advice without any other restriction.}
#' \item{FroeseHCR:} {The HCR defined in the paper by Froese, Branch et al. in Fish and Fisheries 2011.}
#' \item{ghlHCR:} {The model-free HCR used in the management of greenland-halibut}
#' \item{IcesHCR:} {The HCR used by ICES to generate TAC advice in the MSY framework.}
#' \item{IcesCat3HCR:} {HCR that implements the ICES HCR for Category 3 DLS stocks
#' (see ICES CM 2012/ACOM 68: Category 3 - Method 3.2)}
#' \item{IcesCat3HCR_bsafe_hrcap:} {Alternative approach for IcesCat3HCR with a biomass safeguard and harvest rate caps
#' (from ICES WKDLSSLS2019)}
#' \item{little2011HCR:} {The HCR defined in the paper by Little et al. in ICES Journal of Marine Science 2011.}
#' \item{MAPHRC:} {The HCR proposed by the EC in the evaluation on multi-annual management plans in 2015.}
#' \item{MultiStockHRC:} {A HCR that produces TAC advice for several stocks simultaneously. It uses a fishing mortality target and an upper bound to conciliate the TAC advices. In the case of stocks without exploitation rate estimates it uses the catch.}
#' \item{neaMAC_ltmp:} {The HCR used in the north-east atlantic mackerel long term management plan. It is a particular case of the IcesHCR.}
#' \item{pidHCR, pidHCRItarg:} {The HCRs defined in the paper by Pomaerede et al. in Aquatic Living Resources 2010.}
#' }
#' The HCRs are documented in detail in the manual of the library.
#'
#' @examples
#'\dontrun{
#' library(FLBEIA)
#' library(FLAssess) # required to use the IcesHCR. Not available for win64
#' library(FLash) # required to use the IcesHCR. Not available for win64
#' library(ggplot2)
#'
#' # Load the data to run FLBEIA in a one stock one fleet example using the HCR used by ICES
#' # in the MSY framework.
#' data(one)
#'
#' oneAdv$TAC[,ac(2009:2025)] <- NA # Put NA-s in the projection years to check how the
#' # function fills the advice object.
#'
#' res <- IcesHCR(oneSt, oneAdv, oneAdvC, 19, 'stk1')
#' # The value printed in the screen is the fishing mortality used in the advice.
#'
#' res$TAC[,'2009'] # The resulting management advice.
#' }
#-------------------------------------------------------------------------------
# HCRs
# - annualTAC.
#
# Dorleta GarcYYYa
# Created: 20/12/2010 13:26:13
# Changed: 20/12/2010 13:26:18
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# annualTAC(stocks, covars, advice, hcr.ctrl)
# year = Assessment year (POSITION).
#
# The targets and constraints will differ iteration by iteration, thus 'fwd'
# must be applied iter by iter.
#-------------------------------------------------------------------------------
# @export
annualTAC <- function(stocks, advice, advice.ctrl, year, stknm,...){
# project the stock 3 years, (current year, TAC year, TAC year + 1 for ssb or biomass constraints).
nyears <- ifelse(is.null(advice.ctrl[[stknm]][['nyears']]), 3, advice.ctrl[[stknm]][['nyears']])
wts.nyears <- ifelse(is.null(advice.ctrl[[stknm]][['wts.nyears']]), 3, advice.ctrl[[stknm]][['wts.nyears']])
fbar.nyears <- ifelse(is.null(advice.ctrl[[stknm]][['fbar.nyears']]), 3, advice.ctrl[[stknm]][['fbar.nyears']])
f.rescale <- ifelse(is.null(advice.ctrl[[stknm]][['f.rescale']]), TRUE, advice.ctrl[[stknm]][['f.rescale']])
# disc.nyears <- ifelse(is.null(advice.ctrl[[stknm]][['disc.nyears']]), wts.nyears, advice.ctrl[[stknm]][['disc.nyears']])
stk <- stocks[[stknm]]
stk@harvest[stk@harvest < 0] <- 0.00001
# if(dim(stk@m)[1] == 1){
# for(sl in c("catch.n","catch.wt","discards.n","discards.wt","landings.n",
# "landings.wt","stock.n","stock.wt","m","mat","harvest","harvest.spwn", "m.spwn"))
# {
# dimnames(slot(stk,sl))[[1]] <- 1
# }
# stk@range[6:7] <- 1
# }
#
# stk <- stf(stk, nyears = nyears, wts.nyears = wts.nyears, fbar.nyears = fbar.nyears, f.rescale = TRUE) #, disc.nyrs = disc.nyears)
#
ageStruct <- ifelse(dim(stk@m)[1] > 1, TRUE, FALSE)
if(ageStruct == TRUE)
stk <- stf(stk, nyears = nyears, wts.nyears = wts.nyears, fbar.nyears = fbar.nyears, f.rescale = f.rescale) #, disc.nyrs = disc.nyears)
else
stk <- stfBD(stk, nyears = nyears, wts.nyears = wts.nyears, fbar.nyears = fbar.nyears)
fwd.ctrl <- advice.ctrl[[stknm]]$fwd.ctrl
iter <- dim(stk@m)[6]
yrsnames <- dimnames(stk@m)[[2]]
yrsnumbs <- as.numeric(yrsnames)
assyrname <- yrsnames[year]
assyrnumb <- yrsnumbs[year]
# Refresh the years in fwd!!
fwd.ctrl@target$year <- fwd.ctrl@target$year + assyrnumb
fwd.ctrl@target$rel.year <- fwd.ctrl@target$rel.year + assyrnumb
TACvar <- FALSE
for(i in 1:iter){
stki <- iter(stk, i)
# if in <year 0> quantity = catch => set TAC in <year 0> in val
if(TACvar == TRUE | any(fwd.ctrl@target[fwd.ctrl@target$rel.year == assyrnumb,'quantity'] == 'catch')){
k <- which(fwd.ctrl@target$rel.year == assyrnumb & fwd.ctrl@target$quantity == 'catch')
fwd.ctrl@target[k,c('min','val', 'max')] <- fwd.ctrl@target[k,c('min','val', 'max')]*advice$TAC[stknm,year,,,,i]
fwd.ctrl@target[k,'rel.year'] <- NA
fwd.ctrl@trgtArray[k,c('min','val', 'max'),] <- fwd.ctrl@trgtArray[k,c('min','val', 'max'),]*c(advice$TAC[stknm,year,,,,i])
TACvar <- TRUE
}
# if(stknm == 'CMON') browser()
if(dim(stki@m)[1] > 1){
# First estimate/extract the SR model and params.
sr.pars <- advice.ctrl[[stknm]]$sr$params # sr parameters if specified.
sr.model <- advice.ctrl[[stknm]]$sr$model # sr model, mandatory.
if(is.null(sr.pars)){ # if params missing => estimate the parameters using the specified years.
if(is.null(advice.ctrl[[stknm]]$sr$years)) which(round(quantSums(stocks[[stknm]]@stock.n))!=0)[1]:(year-1)# yr0 missing => use all data years.
else{
y.rm <- as.numeric(advice.ctrl[[stknm]]$sr$years['y.rm'])
nyrs <- as.numeric(advice.ctrl[[stknm]]$sr$years['num.years'])
sr.yrs <- yrsnames[(year-y.rm-nyrs + 1):(year-y.rm)]
}
rec <- stki@stock.n[1,sr.yrs]
ssb <- ssb(stki)[,sr.yrs]
# if rec.age != 0 adjust rec and ssb.
rec.age <- as.numeric(dimnames(rec)[[1]])
if(rec.age != 0){
rec <- rec[, -(1:rec.age),]
ssb <- ssb[, 1:(dim(ssb)[2] - rec.age),]
}
if(sr.model != 'geomean') sr.pars <- try(params(fmle(FLSR(rec = rec, ssb = ssb, model = sr.model))), silent = TRUE)
if(class(sr.pars) == 'try-error' | sr.model == 'geomean'){
sr.model <- 'geomean'
sr.pars <- c(prod(c(rec))^(1/length(c(rec))))
sr.pars <- FLPar(a = ifelse(is.na(sr.pars), 0, sr.pars))
}
sr1 <- sr.pars
} else { # sr.pars not null
if(i == 1){
sr1 <- iter(sr.pars,i)
}
sr1[] <- iter(sr.pars,i)[]
}
stki <- FLash::fwd(stki, ctrl = fwd.ctrl, sr = list(model =sr.model, params = sr1))
} else {
# Extract the years to calculate the mean historical growth of the stock
if(is.null(advice.ctrl[[stknm]]$growth.years)) growth.years <- max(1,(year - 11)):(year-1)
else{
y.rm <- as.numeric(advice.ctrl[[stknm]]$growth.years['y.rm'])
nyrs <- as.numeric(advice.ctrl[[stknm]]$growth.years['num.years'])
growth.years <- yrsnames[(year-y.rm-nyrs + 1):(year-y.rm)]
}
stki <- fwdBD(stki, fwd.ctrl, growth.years)
}
Cadv <- ifelse(advice.ctrl[[stknm]][['AdvCatch']][year+1] == TRUE, 'catch', 'landings')
yy <- ifelse(slot(stki, Cadv)[,year+1] == 0, 1e-6, slot(stki, Cadv)[,year+1])
advice[['TAC']][stknm,year+1,,,,i] <- yy
# cat('---------------- HCR------------------------\n')
# cat(c(fbar(stki)[,(year-1):year]), '\n')
# cat('-------------------------------------------\n')
# save(stki, file = 'stki.RData')
}
return(advice)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.