#-------------------------------------------------------------------------------
# annexIVHCR
# HCR that implements the HCR described in Annex IV of (EC COM 2010 241)
# ICES categories 6 - 9
#
# Bnow = (Iy-1 + Iy-2)/2
# Bref = (Iy-3 + Iy-4 + Iy-5)/3
#
# TAC[y+1] = gamma*TAC[y]
#
# | 1 + beta Bnow/Bref > 1 + alpha
# gamma = | f(Bnow/Bref) 1-alpha < Bnow/Bref < 1 + alpha
# | 1 - beta Bnow/Bref < 1 - alpha
#
# Base Case: alpha = 0.2 and beta = 0.15
#
# HCR 2 => f(Bnow/Bref) = 1
# HCR 4 => f(Bnow/Bref) = beta/alpha * (Bnow/Bref - 1) + 1
# Note: HCR names from De Oliveira et al. 2010. Technical Background Evaluation of Annex IV Rules.
# ICES CM 2010/ACOM:5828 pp.
#
# * A single biomass index per stock.
#
# 13/09/2011 09:47:24
# Changed: Dorleta Garcia 15/04/2013 14:23:06 generalized to work with any of
# the index contained in FLIndices objetct as it was it only
# worked with the first in the list.
#-------------------------------------------------------------------------------
#' @rdname annualTAC
#' @aliases annexIVHCR
annexIVHCR <- function(indices, advice, advice.ctrl, year, stknm,...){
Idnm <- advice.ctrl[[stknm]][['index']] # either the name or the position of the index in FLIndices object.
Id <- indices[[stknm]][[Idnm]]@index
ni <- dim(Id)[6]
# Year => Character, because the year dimension in indices does not coincide with year dimension in biol.
year.or <- year
yrnm <- dimnames(advice$TAC)[[2]][year]
year <- which(yrnm == dimnames(Id)[[2]])
Bnow <- (yearSums(Id[,(year-1):(year-2)])/2)[drop=T] # [it]
Bref <- (yearSums(Id[,(year-3):(year-5)])/3)[drop=T] # [it]
Brat <- Bnow/Bref [drop=T] # [it]
alpha <- advice.ctrl[[stknm]][['ref.pts']]['alpha',] #[it]
beta <- advice.ctrl[[stknm]][['ref.pts']]['beta',] #[it]
type <- advice.ctrl[[stknm]][['type']]
# if(Brat > rep(1+alpha,ni))
# gamma <- 1 + beta
# else{
# if(Brat < rep(1-alpha,ni))
# gamma <- 1 - beta
# else{
# if(type == 2) gamma <- 1
# else{ # type == 4
# gamma <- beta/alpha * (Brat - 1) + 1
# }
# }
# }
if(type == 2){
gamma <- ifelse(Brat > 1+alpha, 1 + beta,
ifelse(Brat < 1-alpha, 1 - beta, 1))
}
if(type == 4){
gamma <- ifelse(Brat > 1+alpha, 1 + beta,
ifelse(Brat < 1-alpha, 1 - beta, beta/alpha * (Brat - 1) + 1))
}
advice$TAC[stknm,year.or+1,] <- advice$TAC[stknm,year.or,]*gamma
return(advice)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.