FLBEIA_ShortLivedSmallPelagic_MSE.Rmd

---
title: "Conditioning FLBEIA based on life-history traits: data-limited short-lived small pelagics example"
# author: "Sonia Sanchez-Maroño and FLBEIA team"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
github_document:
  mathjax: TRUE
pdf_document:
  fig_width: 6 
  fig_height: 4 
  toc: yes
tags: [FLBEIA data limited short-lived small pelagics]
license: Creative Commons Attribution-ShareAlike 4.0 International Public License
bibliography: bibliography.bib
---

```{r, ini, echo=FALSE, results='hide', message=FALSE}
# This chunk set the document environment, so it is hidden
library(knitr)
source("R/ini.R")
set.seed(1423)
```

```{r echo=FALSE, out.width='20%'}
include_graphics('images/FLBEIA_logo.png')
```



$~$

# Introduction

In this document we show how to:

1. Condition **FLBEIA** using life-history traits and an assumption on initial explotation level.
2. Test the performance of a specific HCR under different management calendars using **FLBEIA** framework. 

The first part of the tutorial deals with how to condition FLBEIA based on life-history traits and 
the exploitation level at the begining of the simulation period. We will consider a historical period
of 30 years, starting with a virgin population of 100,000 tons and assuming an initial exploitation as
follows:

* first 10 years: linear increase of fishing mortality levels from 0 to $F_{MSY}$

* next 20 years: constant fishing mortality at $F_{MSY}$

Next `FLBEIA` is run using the ICES Annex IV rule for category 3 stocks (i.e. stocks for which 
survey-based assessments indicate trends). 
In the Annex IV rule, the advice is based on the most recent advice. Generally, catch or landings data are 
adjusted to change in the abundance index for the two most recent values relative to the three preceding
values (known as the 2-over-3 rule). Other reference years may be used, based on the knowledge of 
the biology of the stock (e.g. species with a relatively large longevity) or the quality of the data.

In present tutorial, we will specifically use the rule that compares the last available index 
with the mean of the two previous ones (known as the 1-over-2 rule), as it has been assessed 
to be more appropriate for short-lived stocks [@wkdlssls2019, @wkdlssls2020].

The rule is going to be applied under three alternative management calendars with the following formulations:

<!-- ```{r fig:calendars, echo=FALSE, fig.cap="Alternative management calendars: depending on the timings of the survey, stock assessment and management. Source: @sanchez2021."} -->
<!-- include_graphics('images/flbeia_slsp/calendars.jpg') -->
<!-- ``` -->

* Interim year calendar: 
This is the usual management calendar, where the TAC for year y+1 (January-December) is based on 
an abundance index in the interim-year y (in the middle of the year: $1^{st}$ July).
So that there is no indication on age 1 abundance in the TAC year, which might be the bulk of 
the population in the case of the short-lived species.

$$TAC_{(Jan-Dec)_{y+1}} = TAC_{(Jan-Dec)_{y}} \cdot \frac{I_y}{\frac{\sum_{i=y-2}^{y-1} I_i}{2}}$$

* In-year calendar:
In this case, the calendar period is changed and the TAC is set for the period July y to June y+1,
based to an abundance index up to year y. So that we have indications on age 1 abundance during the
first half of the management period (the $2^{nd}$ semester of year y), but not during the second one 
(the $1^{st}$ semester of year y+1).

$$TAC_{Jul_{y}-Jun_{y+1}} = TAC_{Jul_{y-1}-Jun_{y}} \cdot \frac{I_y}{\frac{\sum_{i=y-2}^{y-1} I_i}{2}}$$

* Full population advice: 
For this advice type, management calendar goes from January to December, and we have available 
an abundance index up to year y+1 (specifically B1+ at the beginning of the year), 
which provides information on all the exploitable age classes.

$$TAC_{(Jan-Dec)_{y+1}} = TAC_{(Jan-Dec)_{y}} \cdot \frac{I_{y+1}}{\frac{\sum_{i=y-1}^{y} I_i}{2}}$$

In order to cover this variety of calendars and given the fact that we are going to assume that 
spawning occurs in the middle of the year, we will simulate the stock dynamics in two seasons 
within a year. As a consequence, we need to do an assumption on how the TAC is shared by semester. 
In this case, we will assume that 50% of the catches will be taken in each semester. 



$~$

# Required packages to run this tutorial

To start the R session, we need to load FLBEIA package and some others. 

```{r libs}
library(dplyr)
library(FLBEIA)
library(FLBEIAshiny)
library(FLife)
library(fishmethods)
# library(gtools)      # for rdirichlet function
# library(ggplotFL)
# library(FLBRP)
library(dplyr)       # for manipulating data.frames
library(tidyr)
```



$~$

# Case study definition

As an example, we will condition the Operating Model (OM) for an anchovy-like stock, characterised by 
high natural mortality, full maturity at age 1 and large interannual fluctuations.

The OM is going to be seasonal (in particular by semester), age-structured (ages 0 to 6+) and 
with spawning occurring at the begining of the second semester, so that the recruits (age 0)
turn age 1 at the beginning of the next year.

Firstly, we need define some general parameters, such as the number of iterations, and others related to the 
stock and the spawning moment, that will be required afterwards.

```{r pars}
nit <- 2 # Number of iterations

stkn <- "ANE"               # Stock name
ages <- 0:6                 # age classes
fbar.ages <- c(min=1,max=3) # ages for Fbar calculation

spwn <- 0.5 # percentage of year before spawning (i.e. spawning time: 1st July)

# number of years in the historic period
init.nyr <- 30
```

The specific life-history parameters adopted will be detailed in the next sections.


$~$

## Life-history parameters

### Growth

The individual growth is determined by the growth in length and by the weight-at-age relationship.

Growth is modelled using the von Bertalanffy growth model, with parameters extracted from [@bellido2000]:

$$L_{t} = L_{\infty} \cdot \left( 1 - e^{-K \cdot (t-t0)} \right),$$

where $L_{t}$ is the length at age (t, years); $L_{\infty}$ is the asymptotic average length 
(i.e., the average size at the maximum age); K is the Brody growth rate coefficient (units: $yr^{−1}$); 
and $t0$ represents the age when the average size is zero.

<!-- #! +++ TO BE DECIDED if including "vonB" and "invVonB" FUNCTIONS IN FLBEIA +++ -->

```{r vonB}
# Von Bertalanffy model: linf; k; t0
# - length-at-age
vonB    <- function(age, linf, k, t0) return(linf*(1-exp(-k*(age-t0))))
# - age for specific length
invVonB <- function (L, linf, k, t0)  return(t0 - 1/k * log(1-L/linf)) 

linf <- 18.69
k    <- 0.89
t0   <- -0.02
```

Next, lengths are converted to weight-at-age using the traditional length-weight model:
$$W = a\cdot L^b,$$ 
where $W$ represents the weight, $L$ the length and $a$ and $b$ are the parameters
that need to be estimated. In our case, we took them from `teleost` object in the FLife library [@flife].

```{r lwrel}
# Length-weight (a*L^b): a, b parameters
lw_a <- 0.004799048 # or teleost["a","Engraulis encrasicolus",drop=TRUE]
lw_b <- 3.134380952 # or teleost["b","Engraulis encrasicolus",drop=TRUE]
```

Mean weights-at-age are calculated for different moments in the year:

* at the beggining of each season, for setting the weights-at-age in the stock in each semester; and

* at the middle of each season, for setting the weights-at-age in the catches in each semester.

Therefore, to calculate the mean weights-at-age, first we calculate the mean length-at-age using the 
Von Bertalanffy model and then we transform them into weights using the length-weight relationship.


```{r wage}
# Mean weight at age: from Von Bertalanffy + length-weight relationship

# - stock
#   1st semester  (weights at age at the begging of the year)
ages_stks1 <- ages
mla_stks1 <- vonB(ages_stks1, linf, k, t0)
mwa_stks1 <- lw_a*mla_stks1^lw_b
#   2nd semester (weights at age at spawning time)
ages_stks2 <- ages + spwn
mla_stks2 <- vonB(ages_stks2, linf, k, t0)
mwa_stks2 <- lw_a*mla_stks2^lw_b

# - catch
#   1st semester (weights at age in the middle of the season)
ages_cats1 <- ages + 0.25
mla_cats1 <- vonB(ages_cats1, linf, k, t0)
mwa_cats1 <- lw_a*mla_cats1^lw_b
#   2nd semester (weights at age in the middle of the season)
ages_cats2 <- ages + 0.75
mla_cats2 <- vonB(ages_cats2, linf, k, t0)
mwa_cats2 <- lw_a*mla_cats2^lw_b

# All together
dat <- rbind( data.frame( type = "stock", season = 1, age=ages, mwa=mwa_stks1), 
              data.frame( type = "stock", season = 2, age=ages, mwa=mwa_stks2),
              data.frame( type = "catch", season = 1, age=ages, mwa=mwa_cats1), 
              data.frame( type = "catch", season = 2, age=ages, mwa=mwa_cats2))

ggplot(dat, aes( age, mwa, colour=type)) +
  geom_line() + 
  facet_grid(season ~.) + 
  labs(y = "mwa (gr)")
```

In case of lacking information on $K$, it can be aproximated based on $L_{\infty}$, 
whereas $t0$ can be aproximated based on $L_{\infty}$ and $K$ values, as follows:

* @gislason2008 provides 3 alternatives: 
  + for all combined species: $K = 3.15 L_{\infty}^{-0.64}$
  + for demersal species:     $K = 3.07 L_{\infty}^{-0.64}$
  + for pelagic species:      $K = 6.48 L_{\infty}^{-0.83}$

* @flife: $$t0 = - e^{- 0.3922 - 0.2752 \log(L_{\infty}) - 1.038 \log(K)}$$


### Maturity

We consider a knife-edge maturity, with full maturation at age 1 or older and model it 
using the `knife` function from FLife package [@flife].

```{r knifemat}
# Knife-edge maturity: a1
a1_mat <- 1 # age at full maturity

mata <- knife( age=FLQuant(ages,dimnames=list(age=ages)), 
                params=FLPar(a1=a1_mat))[drop=TRUE]

mata
```

Alternative formulations could be used, such as the logistic maturity ogive, available at 
FLife package [@flife] in `logistic` function. If neither $l50$ nor $a50$ are available,
they can be calculated by an approximation based on $L_\infty$ [@flife]. For example, 

```{r logmat}
# Logistic model: a50, ato95, asym
l50_mat   <- 0.72 * linf^0.93 # L50: length at which 50% of individuals are mature
a50_mat   <- invVonB(l50_mat, linf, k, t0) # a50: age at which 50% of individuals are mature
ato95_mat <- 1                # lag between a50 and first age with full maturity
asym_mat  <- 1                # fixed maturity for older ages

matb <- logistic( age=FLQuant(ages,dimnames=list(age=ages)),
                  params=FLPar(a50=a50_mat, ato95=ato95_mat, asym=asym_mat))[drop=TRUE]
matb
```

### Natural mortality

Alternative values for natural mortality can be estimated using the function `M.empirical` 
in fishmethods package [@fishmethods] based on the life-history parameters defined above and 
other parameters, such as the sea water temperature.
For this case, we estimate variable natural mortality-at-age, based on @gislason2010.

```{r mage}
# Natural mortality
# Linf  : Length-infinity value from a von Bertalanffy growth curve (total length-cm
# Kl    : growth coefficient (per year) from a von Bertalanffy growth curve for length
# Bl    : body length in cm
# method: 9 for variable M@age based on Gislason et al. (2010)

# calculate mean weights at age in the middle of the year
# (with some correction for age 0, as only existing in the second semester)
ages_yr <- ages+c(0.75,rep(0.5,length(ages)-1)) # mid-year
mla <- vonB(ages_yr, linf, k, t0)

Mage <- mla * NA
for (i in 1:length(Mage))
  Mage[i] <- M.empirical(Linf = linf, Kl = k, Bl=mla[i], method = 9)
Mage
```

However, alternative values could be estimated using other methods, as for example global values 
for all age classes. See detailed list of methods and required input parameters by typing 
`?M.empirical`. 
Different values can be estimated with the parameters available from now, so that we can take 
for example the mean of them discarded those considered irrealistic (e.g., Roff's method estimate).

```{r mempir, eval=FALSE}
# Constant M based on tmax: Hoening (1983)
# TC  : mean water temperature (Celsius) experienced by the stock
# tmax: oldest age observed for the species
M.empirical(Linf = linf, Kl = k, 
            TC = 16,   # assumption from http://www.watertemperature.org/Bay-of-Biscay--Geo.html
            tmax = ages[length(ages)],
            tm = 0.5, method = c(1, 3, 4, 5, 10, 11))
```


### Productivity

The stock productivity is defined by the stock-recruitment relationship (SRR). But fitting a SRR is 
not viable in the cases when there is not an assessment of the stock. 
Nevertheless, we can construct an scenario on SRR using the parameterization of the SRR which uses
steepness ($s$, proportion of  recruits produced by 20% of the virgin spawning stock), 
virgin biomass ($B0$) and spawning per recruit ($spr0$) parameters.

In this case, we assumed:

* steepness of 0.75 (that is, medium productivity, see @jardim2015), as high steepness value is indicative of a resilient population (REF Subbey) as the anchovy stocks seem to be

* virgin biomass: 10,0000 tonnes, as in this case we were interested in comparing the performance of different HCRs so that the absolute levels were not important.

* spawning per recruit: was calculated by estimating R0 as a function of B0 (as R0/B is constant).

```{r srrpars}
# Beverton-Holt SR
# s    : steepness (proportion of  recruits produced by 20% of the virgin spawning stock)
#        High steepness --> resilient population.
# v    : virgin biomass (e.g 20 times maximum observed catches)
# spr0 : spawners per recruit at F=0
s  <- 0.75
B0 <- 100000
```

In order to estimate the spawning per recruit, we create a function that calculates expected numbers-at-age
under equilibrium (i.e. F=0), based on assummed recruitment in numbers-at-age.

<!-- #! +++ TO BE DECIDED if including "nage" FUNCTION IN FLBEIA +++ -->

```{r nagefun}
# Function to estimate biomass at age given:
# - ages: ages values (numeric vector)
# - Mage: maturity at age (numeric vector)
# - R0  : initial recruits in mass (numeric)
nage <- function(ages, Mage, R) {
  
  N <- numeric(length(ages))
  
  # - age 0
  N[1] <- R
  
  # - intermediate ages
  for (a in 2:length(ages))
    N[a] <- N[a-1] * exp(-Mage[a-1])
  
  # - plusgroup
  N[length(ages)] <- sum(N[length(ages)] * 
                           exp(-(0:(100-ages[length(ages)]))*Mage[length(ages)]))
  
  return(N)
  
}
```

On the one hand, we need to calculate the spawning per recruit, which is a constant factor 
corresponding to the ratio SSB0/R0.
Based on this value, then we can obtain the R0 that leads us to the selected B0 value. 
Assuming an initial R0 of 1000 tons, we get the R0 corresponding to a virgin biomass at 100,000 tonnes.

```{r spr}
mwa_spwn <- mwa_stks2/1000 # spawning at the beginnig of the 2nd semester

# Estimate R0 as a function of B0 (as R0/B is constant)
R0ini  <- 1000 # should be at spawning time (i.e. 1st July, as spwn=0.5) 
               # --> Some alternatives:
               #      - a = sr_params["a"] * exp(-Mage[1]/2) to move to spawning time; or 
               #      - Mage[1] = Mage[1]/2 for generating Nini
Mage[1] <- Mage[1]/2
Nini    <- nage(ages, Mage = Mage, R = R0ini) # expected numbers under equilibrium, given R0 (thousands)
SSBini  <- sum(Nini * exp(-spwn*Mage) * mata * mwa_spwn) # ssb in tons
spr0    <- SSBini/R0ini # ratio R0/B is constant
R0      <- B0 * 1/spr0  # as our initial SSB is B0
```

Now we transform this parameters (`s`, `v` and `spr0`) to the classical parameterization of 
Beverton & Holt SR model (with `a` and `b` parameters, $rec = a * \frac{ssb}{b + ssb}$ 
using `abPars` function in FLCore package [@flcore].
And we save this information in an `FLSR` object, that will be afterwards used to build the 
`FLSRsim` class used in `FLBEIA` to simulate the recruitment in age structured stocks.

```{r srr}
sr_model  <- "bevholt"

# Use 'abPars' function in FLCore to obtain the parameters of the 
# classical parameterization of beverton & holt SR model.
sr_params <- unlist(abPars(s = s, v = B0, spr0 = spr0, model = sr_model))

# Create FLSR object
srm <- FLSR(name = stkn, params = FLPar(unlist(sr_params)), model = sr_model)
```

To show graphically the fitted SR:

```{r srrPlot}
sr.dat <- data.frame( ssb = seq(0,1e+05,1000)) %>% 
  mutate(rec = sr_params[["a"]] * ssb / (sr_params[["b"]] + ssb))

ggplot(sr.dat, aes(ssb, rec)) + 
  geom_line() + 
  geom_hline(yintercept=R0, linetype = "dashed", color = "red") + 
  geom_text(aes(0,R0,label = "R0", vjust = 1.5), color = "red")

```


## Selectivity

As we have no information on fishing patterns, selectivity was assumed to be equal to maturity in this case.

```{r sel}
sela <- mata
```

Below, several examples on how to estimate alternative selectivity patterns given life-history parameters.

```{r selalt}
# double normal selectivity ogive (from FLife - dnormal function)
a1_sel <- 4     # age for full selectiviy
sl_sel <- 2     # age lag between a50 and age at full selectivity
sr_sel <- 50000 # shaping parameter
sela1 <- dnormal( age=FLQuant(ages,dimnames=list(age=ages)),
                 params=FLPar(a1=a1_sel, sl=sl_sel, sr=sr_sel))[drop=TRUE]
sela1

# Knife-Edge
sela2 <- knife( age=FLQuant(ages,dimnames=list(age=ages)),
                params=FLPar(a1=round(a50_mat,0)))[drop=TRUE] #! to be decided a1 value
sela2

# Logistic (same as maturity pattern)
sela3 <- logistic( age=FLQuant(ages,dimnames=list(age=ages)),
                   params=FLPar(a50=a50_mat, ato95=ato95_mat, asym=asym_mat))[drop=TRUE]
sela3
```

For checking documentation on the different functions, use the `?` command followed by the name of the
function (e.g. `?dnormalFn`)


$~$

## Initial population

In order to estimate initial numbers-at-age, we need to calculate the pristine biomass using the selected
SRR and the natural mortality.

```{r pristineB}
pristineN <- nage(ages, Mage = Mage, R = R0) # R = sr_params["a"] aproximation

# CHECK
B0                                                 # virgin biomass
sum(pristineN * exp(-spwn*Mage) * mwa_spwn * mata) # pristine SSB in tons (if numbers in thousands)
```

At that moment, we have enough information to start conditioning the population. 
So, we generate an `FLStock` object with the previous estimates, that will serve to complete the 
different FLBEIA input objects.


```{r flstock}
# with 2 seasons and 2 iterations
stk <- FLStock( name = stkn, 
                stock.wt = FLQuant(NA, dim = c(length(ages),init.nyr,1,2,1,nit), 
                                   dimnames = list(age=ages,year=1:init.nyr,iter=1:nit)))

# nit iterations
stk <- propagate( stk, nit)

units(harvest(stk))  <- "f"
units(stock.wt(stk)) <- "kg"

# Mean weight at age by semester (in kg)
stock.wt(stk)[,,,1,] <- mwa_stks1/1000
stock.wt(stk)[,,,2,] <- mwa_stks2/1000

# F range to estimate mean F
range(stk)[c("minfbar","maxfbar")] <- fbar.ages

# Prop of spawn before M and F
ssb.ss <- rec.ss <- 2
# (assumed spawning in the middle of the year)
if (spwn <= 0.5) {
  harvest.spwn(stk)[,,,1,] <- m.spwn(stk)[,,,1,] <- spwn/0.5
  harvest.spwn(stk)[,,,2,] <- m.spwn(stk)[,,,2,] <- 0
}

# Selectivity
harvest(stk) <- sela
harvest(stk)[,1,] <- 0  # no catches in pristine status

# Natural mortality and maturity
m(stk)     <- Mage/2
m(stk)[1,,,1,] <- 0    # different for age 0, as rec occurs mid-year (not 1st Jan)
m(stk)[1,,,2,] <- Mage[1]
mat(stk) <- mata

# Initial numbers: virgin biomass
stock.n(stk)[ ,1,,1,] <- pristineN
stock.n(stk)[ ,1,,2,] <- stock.n(stk)[,1,,1,] * exp(-m(stk)[,1,,1,])
stock.n(stk)[1,1,,1,] <- 0 # recruitment occurs in the 2nd semester
units(stock.n(stk)) <- "1000"
stock(stk) <- quantSums(stock.n(stk) * stock.wt(stk))

# Catches: 0 in the 1st simulation year
catch.wt(stk)[,,,1,] <- landings.wt(stk)[,,,1,] <- discards.wt(stk)[,,,1,] <- mwa_cats1/1000
catch.wt(stk)[,,,2,] <- landings.wt(stk)[,,,2,] <- discards.wt(stk)[,,,2,] <- mwa_cats2/1000
```


$~$

## Reference points

As we have adopted the assumption that catches are equally shared between semesters (i.e., 50% each), 
we need to estimate the reference points as a function of the percentage of F in each semester.
So that we can select those according to a 50% share.

```{r msyruns}
# take 1st year with information on selectivity (i.e. year=2)
stk1 <- iter(window(stk, start = 2, end = 2), 1)

# BRPS estimation
#------------------

msyruns <- brps <- data.frame()

# Loop by percentage of catches in each season

for (p in seq(0.1,0.9,0.1)) {
  
  print(paste("p =",p))
  
  # Estimate reference points
  fruns <- brpsson( stk1, B0=B0, R0=R0, rec.ss=rec.ss, ssb.ss=ssb.ss, 
                    sr_model="bevholt", sr_params=sr_params, 
                    Fprop = c(p, 1-p))
  
  msyruns <- bind_rows(msyruns, fruns$runs)     # simulations
  brps <- bind_rows(brps, fruns$refpts) # brps values
  
}
```

We get two types of outputs, the simulation values `msyruns` and the reference points' estimates
`brps`. There is the possibility to see the simulated values graphically by using `plotBRPsson` function.

```{r msyrunsfig, eval=FALSE}
plotBRPsson( msyruns, pdfnm="ane_Fbar_vs_SPR.pdf")
```

Which generates a pdf file with YPR plots by percentage of F value in each semester:

```{r echo=FALSE, fig.show = "hold", out.width = "50%", fig.align = "default"}
include_graphics('images/flbeia_slsp/ypr1.png')
include_graphics('images/flbeia_slsp/ypr2.png')
```

Regarding the reference points, we obtain a data.frame with one line per F reference point type and F proportion.
Therefore, we need to filter those cases that are more in line with the assumption that 50% of catches are obtained 
in each semester.

```{r brpprint}
brps$refpt <- factor( brps$refpt,
                      levels = c("virgin", "msy", "f0.1", "spr.20", "spr.30", "spr.35", "spr.40", "spr.50",
                                 "F.20B0", "F.30B0", "F.35B0", "F.40B0", "F.50B0", "F.50R0", "F.90msy"))

brps <- brps %>%
  group_by(refpt) %>%
  filter(abs(pY_s1 - 0.5) == min(abs(pY_s1 - 0.5), na.rm = TRUE)) %>%
  arrange(refpt) %>%
  filter(refpt != "virgin")

as.data.frame(brps)
```

Now we need to select our $F_{MSY}$ proxy, as $F_{MSY}$ is not considered value for short-lived stocks
(as usually allows very high exploitation near to $F_{20\%B0}$,  which might be used as $B_{lim}$).
A proxy for $B_{lim}$ could be the biomass that leads to 0.7R0 according to the SR model [wkmsycat342017, wklife2018], 
in our case both approaches are in agreement, as $F_{20\%B0}$ is the biomass that leads to 0.75R0.

For this case, we decided to base our $F_{MSY}$ proxy on $F_{40\%B0}$ [@punt2014] (that is, 
the fishing mortality rate associated with a biomass of 40% B0 at equilibrium), 
which corresponds to $F=1.18$.

```{r brp}
Blim <- 0.20 * B0

F40B0 <- brps %>% filter(refpt=="F.40B0")
F40B0$fbar # 1.18

ref.pts <- c( B0 = B0, R0 = R0, Blim = Blim, Fmsy = F40B0$fbar, Bmsy=F40B0$ssb, MSY=F40B0$yield)

# F proportion in each semester (for 50% catches in each semester)
fprop_s1 <- F40B0$fprop_s1
```


$~$

# FLBEIA conditioning

With the information generated up to now (mainly the `FLStock` object), we are going to construct 
the different objects that will be used as FLBEIA inputs.

Firstly, we define general parameters. These are: simulation years (with years 1-30 for the 
historical period and 30 additional years, 31-60, for the projection period), ages in the stock, 
number of seasons and iterations.

```{r dims}
# years
hist.yrs  <- dimnames(stk)$year
hist.nyr  <- length(hist.yrs)

first.yr <- as.numeric(hist.yrs[1])
proj.yr  <- as.numeric(hist.yrs[hist.nyr])+1
proj.nyr <- 30
last.yr <- proj.yr+(proj.nyr-1)

hist.yrs <- first.yr:(proj.yr-1)
proj.yrs <- proj.yr:last.yr
yrs <- first.yr:last.yr
nyr <- length(yrs)

# ages
ages <- dimnames(stk@stock.n)$age
na   <- length(ages)

# seasons
ns  <- dim(stk)[4]

# iterations
nit  <- dim(stk)[6]
```

Next, we expand the `FLStock` object to cover, not only the historical period, but also the projection
years.

```{r stkyrs}
stk <- window( stk, start=first.yr, end=last.yr)
```

And finally we will create all the arguments necessary to run FLBEIA. We can see the name of the objects 
needed to run it using the function `args`. We can obtain further information on the objects using the 
FLBEIA help page (`?FLBEIA`).

```{r args}
args(FLBEIA)
```


$~$

## Data objects

To start with, we create some empty `FLQuants` (the basic FLR data structure) with the dimensions to be used 
in the case study (CS) that will aid the conditioning process.
The first two have no age dimension, while the last ones do.

```{r flqs}
flq  <- FLQuant(1, dim = c(1,nyr,1,ns,1,nit), 
                dimnames = list(quant = "all", year = yrs, iter = 1:nit))
flq0 <- FLQuant(0, dim = c(1,nyr,1,ns,1,nit), 
                dimnames = list(quant = "all", year = yrs, iter = 1:nit))
flqa  <- FLQuant(1, dim = c(na,nyr,1,ns,1,nit), 
                 dimnames = list(age = ages, year = yrs, iter = 1:nit))
flqa0 <- FLQuant(0, dim = c(na,nyr,1,ns,1,nit), 
                 dimnames = list(age = ages, year = yrs, iter = 1:nit))
flq1s <- FLQuant(1, dim = c(1,nyr,1,1,1,nit), 
                 dimnames = list(quant = "all", year = yrs, iter = 1:nit))
```


### biols (OM: "real" stock)

The `FLBiols` object is a named list of `FLBiol` objects with the name of the stocks represented in each component.
Each FLBiol object represents a different stock simulated in the biological OM, corresponding with the "true"
population of the MSE simulations.

We will fill the object using the information from the `FLStock` previously created.


```{r biols}
# age structured stock
minage  <- as.numeric(ages[1])
maxage  <- as.numeric(ages[na])
minfbar <- range(stk)[["minfbar"]]
maxfbar <- range(stk)[["maxfbar"]]

biols <- FLBiols( stk = FLBiol(name = name(stk), 
                               desc = "",
                               range = c(min = minage, max = maxage, plusgroup = maxage, 
                                         minyear = first.yr, maxyear = last.yr, 
                                         minfbar = minfbar, maxfbar = maxfbar),
                               n    = stk@stock.n, 
                               wt   = stk@stock.wt, 
                               fec  = predictModel(fec = flqa, model = ~ fec), 
                               mat  = predictModel(mat = stk@mat, model = ~ mat), 
                               m    = stk@m,
                               spwn = stk@m.spwn
))

names(biols) <- name(stk)

# projection period
m(biols[[1]])[,ac(proj.yrs)]   <- m(biols[[1]])[,ac(proj.yr-1)]
fec(biols[[1]])[,ac(proj.yrs)] <- fec(biols[[1]])[,ac(proj.yr-1)]
mat(biols[[1]])[,ac(proj.yrs)] <- mat(biols[[1]])[,ac(proj.yr-1)]
wt(biols[[1]])[,ac(proj.yrs)]  <- wt(biols[[1]])[,ac(proj.yr-1)]

spwn(biols[[1]])[,ac(proj.yrs)]<- spwn(biols[[1]])[,ac(proj.yr-1)]

# replace 0 to a very small value to avoid NaN when estimating F at C=0
biols[[1]]@m[1,,,1,] <- 1e-08
```


### SRs (OM: stock-recruitment)

The `FLSRsim` object is used to store the information necessary to simulate recruitment for age-structured
populations. We will introduce uncertainty in the generated recruitments (both in historical and projection 
periods) with a log-normal distribution with median equal to 1 and a standard deviation of 0.75.
Autocorrelation in residuals was not included.


```{r SRs}
SRs <- list(stk = FLSRsim( name = name(stk), desc = "", 
                           ssb=flq, model = sr_model))
names(SRs) <- name(stk)

# both spawning and recruitment assumed at the beginning of the second semester
ssb.ss <- 2
rec.ss <- 2

# - data

SRs[[1]]@ssb[]            <- ssb(stk)
SRs[[1]]@ssb[,,,-ssb.ss,] <- 0

SRs[[1]]@rec[]                    <- stock.n(stk)[1,]
SRs[[1]]@rec[,,,c(1:(rec.ss-1)),] <- 0

# - rectruiment distribution

SRs[[1]]@proportion[,,,-rec.ss,] <- 0
SRs[[1]]@proportion[,,, rec.ss,] <- 1 # recruits generated only in one season

# - time lag

SRs[[1]]@timelag['year',]         <- dims(biols[[1]])$min  # =0, because age 0 recuitment
SRs[[1]]@timelag['season',rec.ss] <- ssb.ss

# - parameters

SRs[[1]]@params[] <- sr_params


# - uncertainty

sigR <- 0.75
yr <- ac(hist.yrs[-1],proj.yrs)
SRs[[1]]@uncertainty[,yr,,rec.ss,] <- rlnorm( length(yr) * nit, 0, sigR)
```


### fleets (OM: "real" fleet)

The `FLFleetsExt` object is a named list of `FLFleetsExt` objects with the name of the fleets represented 
in each component. Each FLFleetsExt object represents a different fleet simulated in the OM, corresponding 
with the "true" dynamics of the fleets in the MSE simulations.

For present study, we will consider that there is only one fleet with one métier. 
We do it in two steps, by firstly creating the `FLCatchExt` object (containing information on méiter's catches) 
and next the `FLFleetsExt` (which contains information on effort, capacity, costs...). 
However, we are not going to include economic information, as it is not under the scope of present anaysis.

```{r fleets}
# catch object

catch.obj <- FLCatchExt( name = name(stk), 
                         alpha = flqa, beta = flqa, landings.n = stk@catch.n,
                         landings = stk@catch, landings.n = stk@catch.n, 
                         landings.wt = stk@landings.wt, 
                         discards = flq0, discards.n = flqa0, 
                         discards.wt = stk@discards.wt, 
                         landings.sel = flqa, discards.sel = flqa0)

landings.wt(catch.obj)[,ac(proj.yrs)]  <- landings.wt(catch.obj)[,ac(proj.yr-1)]
discards.wt(catch.obj)[,ac(proj.yrs)]  <- discards.wt(catch.obj)[,ac(proj.yr-1)]

catches.obj <- FLCatchesExt(stk = catch.obj)
names(catches.obj) <- name(stk)


# fleets

fleets <- FLFleetsExt(fl = FLFleetExt(name = "fl", effort= flq, capacity = flq*1e12,
                                      metiers =  FLMetiersExt(mt = FLMetierExt(name = "mt", 
                                                                               effshare = flq, 
                                                                               catches = catches.obj))))

fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,ac(hist.yrs),,,,] <- stk@harvest[,ac(hist.yrs)]

fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,ac(proj.yrs)] <- 
  yearMeans(fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,ac(hist.yrs)]) 
```


### covars (OM: other variables of interest)

In this CS, no covariate will be considered.

```{r covars}
covars <- NULL
```


### indices (MP: observation model)

The `FLIndices` object is a named list of `FLIndex` with the indices that will be used to generate advice. 
Each FLFleetsExt object represents a different observed index. The model implemented in FLBEIA to simulate 
abundance indices is the classical linear model with a multiplicative error. 

In present CS, we will have a biomass index observed in the middle of the year and multiplicative error
will be generated using a log-normal distribution with median equal to one and coefficient of variation equal 
to 20\%. This index will be available only up to 10 years prior to first projection year.

```{r indices}
idxyr <- (proj.yr-10):last.yr

flqidx <- window(trim(flq, season=ssb.ss), start=idxyr[1], end=idxyr[length(idxyr)])

indices <- list( stk = FLIndices(B1plusIdx = FLIndex(name = name(stk), catch.wt = flqidx, 
                                                     effort = flqidx, index = flqidx)))

names(indices) <- name(stk)

indices[[1]]$B1plusIdx@type <- "biomass"

yr <- idxyr[idxyr %in% hist.yrs]

cvIDX <- 0.25

indices[[1]][[1]]@index.q[] <- 1 * rlnorm(dim(flqidx)[2]*nit, 0, sqrt(log(cvIDX^2+1)))

#! BORRAR???
indices[[1]][[1]]@index[,1,] <- indices[[1]][[1]]@index.q[,1,] * 
  quantSums(wt(biols[[1]][-which(dimnames(biols[[1]])$age == "0"),1,,ssb.ss,]) * 
              n(biols[[1]][-which(dimnames(biols[[1]])$age == "0"),1,,ssb.ss,]))

indices[[1]][[1]]@range[c("startf","endf")] <- 0.5 # B1+ in the middle of the year

```

Next we need to generate the index observation in the last 10 years of the historical period.
However, we fistly we need to simulate the historical population and afterwards index observation will be done 
in a different way depending on the calendar (int & iny calendars: B1+ in the middle of the interim year; 
fpa: B1+ at the beggining of the following year).


### advice (MP: advice rule)

The advice object is a list with information on TACs and fleets' quota shares. However, in this case 
the quota share will be equal to 1 as we have assumed that there is only one fleet exploiting the stock.

```{r advice}
advice <- list(TAC = flq1s,  quota.share = list(stk = flq1s))

dimnames(advice$TAC)$quant <- names(advice$quota.share) <- name(stk)
dimnames(advice$quota.share[[1]])$quant <- names(fleets)

quant(advice$TAC) <- "stock"
quant(advice$quota.share[[1]]) <- "fleet"
```


$~$

## Control objects 

The control objects are lists used to store the values that control how each part of the simulation is carried out.
There is one control object per data object.

### main.ctrl

The `main.ctrl` object declares the initial and final year of the simulations.
We want to set the first TAC in the the first proj.yr (year 31), consequently we have to start projections 
one year before. 
Because FLBEIA in the first projection year (year 30, in this example), just projecs the population forward 
taking into account the TAC set for this year and afterwards (based on stock status) estimates the TAC 
for the following year.

```{r mctrl}
main.ctrl  <- list(sim.years = c(initial = proj.yrs[1]-1, final = proj.yrs[proj.nyr]))
```

### biols.ctrl

In the `biol.ctrl` object we declare the model to be used to carry the population forward in the simulation.
In his case we used `ASPG_Baranov` which stands for Age-Structured Populaton Growth with Baranov catch 
equation.

```{r bctrl}
biols.ctrl <- create.biols.ctrl(stksnames = name(stk), growth.model = "ASPG_Baranov")
```

In case of interested in using Cobb-Douglas equation for generating catches-at-age, then we should select
`ASPG` which assumes that catches are taken in the middle of the season, instead of the asssumption on 
continuous catches along the season in the Baranov equation.


### fleets.ctrl

The `fleets.ctrl` object controls the four processes simulated in the fleet operating model, the effort allocation, the catch production, the price formation and the capital dynamics. 

Catch production is simulated using Baranov catch equation, `Baranov`, (an alternative is to use Cobb-Douglas 
function by age class, `CobbDouglasAge`). 
When simulating  #! A CAMBIAR
Effort allocation is simulated using `SMFB` model, but as there is only one stock and one fleet with a single metier 
the model calculates just the effort that produces exactly the TAC advice for the stock. 
The price and the capital dynamics are maintained fixed. 

We create the object using `create.fleets.ctrl` which by default already sets the seasonal share of the TAC to 
50%, same as we intended.


```{r flctrl}
fleets.ctrl <- create.fleets.ctrl(fls = "fl", fls.stksnames = list(fl = name(stk)), flq = flq, 
                                  effort.models = c(fl = "SMFB"), n.fls.stks = c(fl = 1), 
                                  capital.models = c(fl = "fixedCapital"), 
                                  price.models = c(fl = "fixedPrice"), 
                                  catch.models = c("Baranov"))
```

### covars.ctrl

As there are no covariates, the control parameter is not required.

```{r cvctrl}
covars.ctl <- NULL
```

### obs.ctrl

The `obs.ctrl` object comprises the neccesary information to simulate the observed data used in the management procedure 
(MP) to generate the management advice. In this object we need to define the settings for the observation of the stock 
and the index. 

In this CS, we have a data-limited stock for which we only have an abundance index available.
Consequently, we do not need to observe stock data, and we will use `NoObsStock` option.
The abundance index is generated using the `bio1plusInd` function which updates the index slot on the object with 
the most recent abundance data (in this specific case, focusing on biomass of age 1 and older).


```{r obsctrl}
obs.ctrl <- create.obs.ctrl(stksnames = name(stk), n.stks.inds = c(1), 
                            stks.indsnames = names(indices[[1]]), 
                            stkObs.models = "NoObsStock", 
                            indObs.models = "bio1plusInd")
```

### assess.ctrl

The assess.ctrl object declares the assessment model to be used for each of the stocks and, if required, 
the additional settings needed to run the models. In this case no assessment model is used.

```{r assctrl}
assess.ctrl <- create.assess.ctrl(stksnames = name(stk), assess.models = "NoAssessment")
```

### advice.ctrl

Now we create the advice control object corresponding to the HCR used by ICES for category 3 stock in the 
data limited framework (DLS). The HCR, known as n-over-m rule, is based on the tendency of an abundance index 
in the last $n+m$ years. So that the previous year TAC is corrected by the ratio between the mean of the index 
in the last $n$ years and the mean of the index in the last $m$ years. Plus a maximum change limit in the TAC 
change can be set, both for decrease (lower uncentainty cap) and increase (upper uncentainty cap). 
And the corresponding FLBEIA function is `IcesCat3HCR_bsafe_hrcap` which additionally 
includes biomass safeguards or harvest rate limits.
See details on the function in @wkdlssls2020. 

This rule requires information on: 

* How to initialize the rule: the first year when the rule is applied (`cref.year`) and 
                              the reference TAC set for this year (`cref.value`).

* The years to be used in the ratio (`refnyrs`): $n$ and $m$ values.

* Lower and upper uncentainty caps' (`unccap.low` and `unccap.upp`) as percentages.

* Precautionary buffer (`prec.buffer`): with information on the percentage reduction (`value`) and the 
  years in which it must be applied (`year`).
  
We will also include here the values for the reference points, although these are not required for this 
specific HCR.

```{r advctrl}
advice.ctrl <- list()

# HCR and input index
advice.ctrl[[stkn]] <- list(HCR.model = "IcesCat3HCR_bsafe_hrcap", index = names(indices[[1]]))

# 1-over-2 rule
advice.ctrl[[1]][['refnyrs']] <- c(now=1, past=2)

# uncertainty caps: 80% symmetric
advice.ctrl[[1]]$unccap.low <- 0.80
advice.ctrl[[1]]$unccap.upp <- 0.80

# precautionay buffer: 20% in the 1st simulation year
advice.ctrl[[1]]$prec.buffer <- list( value = 0.20, # if 0 --> no buffer
                                      year = as.character(main.ctrl$sim.years[["initial"]]))

# reference points
advice.ctrl[[1]][['ref.pts']] <- matrix( ref.pts, nrow = length(ref.pts), ncol = nit, 
                                         dimnames = list(names(ref.pts),1:nit))
```

If aiming at includding also a biomass safeguard in the rule, we should also define an extra variable 
which sets the reference level (`Itrigger`). But for setting the trigger level, we previously need to 
generate the index observations for the historical period.

# ```{r advctrlBsafe, eval=FALSE}
# advice.ctrl_bsafe <- advice.ctrl
# 
# # trigger point: e.g. minimun of observed indices
# advice.ctrl_bsafe[[1]]$Itrigger <- 
#   apply(indices[[1]]$B1plusIdx@index, c(1,3:6), min, na.rm=TRUE)[drop=TRUE]
# ```

We are pending to set how to initialize the rule, but previously we need to simulate the historical part.
As the first year we apply the rule we will take as reference TAC the mean catches over a range of the 
last simulated historical years. 



$~$

# Initial exploitation level

We simulated the initial exploitation of the stock, by projecting the population from a virgin 
levels, with increasing exploitation levels during 10 years to reach $F_{MSY}$ 
(i.e. $F_{MSY_{proxy}} = F_{40\%B0}$) and maintaining this MSY level during 20 years more. 
Instead of setting fix F values, variability was included through a lon-normal distribution with CV of 10%.

```{r histF}
# Initial exploitation at Fmsy
Ftgt <- advice.ctrl[[1]][['ref.pts']]["Fmsy",1]

# Introducing uncertainty in the Ftgt (same F at age in both seasons)
cvfh <- 0.10

fyr <- FLQuant( c(seq(0,Ftgt,length.out = 10), rep(Ftgt, length(hist.yrs)-10)), 
                dimnames = list(age="all", year=hist.yrs, iter=1:nit)) * 
  rlnorm(length(hist.yrs) * nit, 0, sqrt(log(cvfh^2+1))) 
```

Given the objective fishing mortality (`Ftgt`) by year and the selectivity pattern for the $F_{bar}$ range, 
we need to calculate the F multiplier for each age class to that leads to $F_{bar}=F_{tgt}$.

```{r fmult}
ages.fbar <- ac(biols[[1]]@range["minfbar"]:biols[[1]]@range["maxfbar"])

msel <- quantMeans(fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[ages.fbar,hist.yrs])

Fmult_s1 <- fyr/msel[,,,1,]; Fmult_s1[is.na(Fmult_s1)] <- 0
Fmult_s2 <- fyr/msel[,,,2,]; Fmult_s2[is.na(Fmult_s2)] <- 0
```

As we have two semesters in the simulations, we needed to make an assumption on the percentage of 
fishing mortality in each semester. We assumed that this percentage was kept constant at the value
leading to 50% of the catches in each semester (0.3 in this case, which has been calculated 
in conjuction with the reference points; see `fprop_s1` value).

```{r fsem}
fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,hist.yrs,,1,] <- 
  sweep(fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,hist.yrs,,1,], 2:6, Fmult_s1, "*") * fprop_s1
fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,hist.yrs,,2,] <- 
  sweep(fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,hist.yrs,,2,], 2:6, Fmult_s2, "*") * (1-fprop_s1)
```

Now we are ready to project the stock and the fleets forward in the historical period, 
based on the expected numbers at age at equilibrium and the simulated fishing mortalities 
for each year and season.

```{r histPop}
for(yr in 1:length(hist.yrs)) {
  
  for (s in 1:ns) {
    
    if (yr > 1 | s > 1) {
      
      yy <- ifelse( s==1, yr-1,  yr)
      ss <- ifelse( s==1,   ns, s-1)
      
      f0 <- fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,yy,,ss,]
      z0 <- biols[[1]]@m[,yy,,ss,] + f0
      
      # nage: adults
      
      if (s==1) { aa <- ages[-c(na-1,na)]
      } else aa <- ages[-c(1,na)]
      
      biols[[1]]@n[-c(1,na),yr,,s,] <- biols[[1]]@n[aa,yy,,ss,]*exp(-z0[aa,])
      
      if (s==1) {
        biols[[1]]@n[na,yr,,s,]     <- biols[[1]]@n[na-1,yy,,ss,]*exp(-z0[na-1,]) + 
                                        biols[[1]]@n[na,yy,,ss,]*exp(-z0[na,])
      } else 
        biols[[1]]@n[na,yr,,s,]     <- biols[[1]]@n[na,yy,,ss,]*exp(-z0[na,])
      
      # nage: recruits
      
      yr.rec <- yr-SRs[[1]]@timelag["year",s]
      ss.ssb <- SRs[[1]]@timelag["season",s]
      
      SRs[[1]]@ssb[,yr.rec,,ss.ssb,] <- 
        quantSums(n(biols[[1]]) * wt(biols[[1]]) * fec(biols[[1]])*mat(biols[[1]]) * 
                    exp(-biols[[1]]@m*spwn(biols[[1]])), na.rm=TRUE)[,yr.rec,,ss.ssb,]
      
      biols[[1]]@n[1,yr,,s] <- 
        predict( predictModel(model=SRs[[1]]@model, params=FLPar(SRs[[1]]@params[,yr,s,])), 
                 ssb=SRs[[1]]@ssb[,yr.rec,,ss.ssb,]) * 
        SRs[[1]]@proportion[,yr,,s,] * SRs[[1]]@uncertainty[,yr,,s,]
      
    }
    
    # cage
    
    f1 <- fleets[[1]]@metiers[[1]]@catches[[1]]@catch.q[,yr,,s,]
    z1 <- biols[[1]]@m[,yr,,s,] + f1
    
    fleets[[1]]@metiers[[1]]@catches[[1]]@landings.n[,yr,,s,] <-
      (f1/z1) * (1-exp(-z1)) * biols[[1]]@n[,yr,,s,]
    
  }
  
}

# landings = sum_age(landings.n * landings.wt)

fleets[[1]]@metiers[[1]]@catches[[1]]@landings[,hist.yrs,] <- 
  quantSums(fleets[[1]]@metiers[[1]]@catches[[1]]@landings.n[,hist.yrs,] * 
              fleets[[1]]@metiers[[1]]@catches[[1]]@landings.wt[,hist.yrs,])
```

We can print the total landings by semester, to check if the condition of 50% of the catches in each
semester is fullfilled.

#! ++++ COMENTAR UNA VEZ CHEQUEADO +++ , eval=FALSE

```{r catchperc}
sweep( fleets[[1]]@metiers[[1]]@catches[[1]]@landings[,hist.yrs,,,], c(1:3,5:6), 
       seasonSums(fleets[[1]]@metiers[[1]]@catches[[1]]@landings[,hist.yrs,]), "/")
```

We assume that historical TACs are equal to yearly catches.

```{r histTAC}
advice$TAC[,hist.yrs,] <- seasonSums(quantSums(catchWStock(fleets, stock = stkn)))[,hist.yrs,]
```

Now we are ready to initialize the rule, for example, we can take the mean catches of the last m 
years for the n-over-m rule (i.e. being m the number of years in the denominator of the HCR).
In our case, the mean of the last 2 years (for the 1-over-2 rule).

```{r hcrInit}
advice.ctrl[[1]]$cref.year <- main.ctrl$sim.years["initial"]
lyr <- main.ctrl$sim.years["initial"] - 1

nyr <- advice.ctrl[[1]][['refnyrs']]["past"] # number of years in rule denominator

advice.ctrl[[1]]$cref.value <- 
  yearMeans(quantSums(seasonSums(catchWStock(fleets, stock = stkn)))[,lyr - c(0:(nyr-1)),])
```


Finally, we estimate the catchability for the Cobb-Douglas equation, based on the simulated catches-at-age. 
In FLBEIA, we have `calculate.q.sel.flrObjs` function that calculates catchability values in a different way 
depending on the catch model we have decided to use for each fleet and stock 
(`fleets.ctrl[[fl]][[st]]$catch.model`):

* if `catch.model=="CobbDouglasAge"`: $$q_a = \frac {C_a} {E^{\alpha} \cdot {B_a}^{\beta}}$$

* if `catch.model=="Baranov"`: $$q_a = \frac{\frac{C_a}{C} \cdot F_a}{E^{\alpha}}$$

Additionally, this function calculates catchability for the simulation period (`sim.yrs`) 
based on the mean of a specific range of years (`mean.yrs`). 


```{r catchq}
hl5.yr <- hist.yrs[(length(hist.yrs)-5+1):length(hist.yrs)] # last 5 years in historical period

fleets <- calculate.q.sel.flrObjs( biols = biols, fleets = fleets, fleets.ctrl = fleets.ctrl, 
                                   mean.yrs = hl5.yr, sim.yrs = proj.yrs)
```



$~$

# Management under alternative calendars

Previosly defined settings are common for all calendars. 
Now we have to define other settings which are specific to the calendar, such as: 

* The TAC for the first projection year, that we will assume to be equal to the simulated catches 
in the period covered by the specific calendar (i.e. January-December or July-June).

* The observed indices in the historical part, as they will be different depending on the calendar:
  + B1+ index in the middle of the interim year (for int and iny calendars); or
  + B1+ index at the beggining of the TAC year (for fpa calendar).

* Observation settings
  + Observation model selected (`obs.ctrl[[st]]$indObs[[idx]]$indObs.model`), which value depends on the 
    selected calendar:
      + `bio1plusInd`, which corresponds to the B1+ at the beginnig of the advice season, should be used for 
       for the int calendar (where the index corresponds to B1+ at the beginning of the second semester, 
       while the advice is given at the end of the year); or
      + `bio1plusfwdInd`, which corresponds to the B1+ at the beginnig of season following to the advice one, 
       for the iny and fpa calendars.
       should be used for iny (where the index corresponds to B1+ at the beginning of the second semester, 
       while the advice is given at the end of the first semester) and for fpa (where the index corresponds to 
       B1+ at the beginning of the first semester next year, while the advice is given at the end of the year) 
       calendars.
  + If there is observation of the index in the assessment year (`obs.ctrl[[st]]$obs.curryr <- TRUE`), 
    as by default indices are observed up to the year previous to the assessment. 
    This is the case for iny and fpa calendars.

* Advice settings
  + The season when the advice is provided (`advice.ctrl[[st]]$adv.season`), that value should be defined 
    in case it is different to the end of the last season of a year.
  + In case there is observation of the index in the assessment year, we need to detail if this last index 
    is going to be used to give advice (`advice.ctrl[[st]]$inyr.idx <- TRUE`).
    This is the case for iny and fpa calendars.

Now we are going to define these specific settings for each calendar, as explained above:

* Interim year calendar

```{r intset}
# Variables to be changed for specific calendar
indices_int  <- indices
advice_int   <- advice
obs.ctrl_int <- obs.ctrl

# Observed indices (in the middle of the interim year)
indices_int[[1]][[1]]@range[c("startf","endf")] <- 0.5 # B1+ in the middle of the year

idxyr <- dimnames(indices_int[[1]][[1]])$year
yr    <- idxyr[idxyr %in% hist.yrs[-length(hist.yrs)]]


a1plus <- which(dimnames(biols[[1]]@n)$age=='1'):dim(biols[[1]]@n)[1]

idxss <- findInterval( mean(range(indices_int[[1]][[1]])[c("startf","endf")]),
                       seq(0,1,length = ns+1))

indices_int[[1]][[1]]@index[,yr,] <- indices_int[[1]][[1]]@index.q[,yr,] *
  quantSums( n(biols[[1]])[a1plus,yr,,idxss,] * 
               wt(biols[[1]])[a1plus,yr,,idxss,])  # B1+ in the middle of the interim year

indices_int[[1]][[1]]@index[,ac(c(main.ctrl$sim.years[["initial"]],proj.yrs)),] <- NA

# Required 1st TAC
advice_int$TAC[,main.ctrl$sim.years["initial"],] <- 
  seasonSums(quantSums(catchWStock(fleets, stock = stkn)))[,main.ctrl$sim.years["initial"],]

# Observation of indices up to TAC year
obs.ctrl_int[[1]]$obs.curryr <- FALSE
```


* In-year calendar

```{r inyset}
# Variables to be changed for specific calendar
indices_iny     <- indices_int # Same index as for the interim-year calendar
advice_iny      <- advice
fleets.ctrl_iny <- fleets.ctrl
obs.ctrl_iny    <- obs.ctrl
advice.ctrl_iny <- advice.ctrl

# Required 1st TAC (July-June catches)
c2s.yr1 <- quantSums(catchWStock(fleets, stock = stkn))[,main.ctrl$sim.years["initial"] -1,,2,]
C1s.yr2 <- seasonMeans(quantSums(catchWStock(fleets, stock = stkn)))[,main.ctrl$sim.years["initial"],,,]
advice_iny$TAC[,main.ctrl$sim.years["initial"] -1,] <- c2s.yr1 + C1s.yr2

# Changed calendar --> need to recalculate historical seasonal shares
fleets.ctrl_iny$seasonal.share[[1]][,main.ctrl$sim.years["initial"] -1,,2,] <- 
  c2s.yr1 / advice_iny$TAC[,main.ctrl$sim.years["initial"] -1,]
fleets.ctrl_iny$seasonal.share[[1]][,main.ctrl$sim.years["initial"],,1,]   <- 
  C1s.yr2 / advice_iny$TAC[,main.ctrl$sim.years["initial"] -1,]

# Observation of indices up to TAC year
obs.ctrl_iny[[1]]$obs.curryr  <- TRUE
obs.ctrl_iny[[1]]$indObs[[1]]$indObs.model <- "bio1plusfwdInd"

# TAC calendar (July y - June y+1)
# -> advice at the end of the 1st semester (not at the end of the year, as default)
advice.ctrl_iny[[1]]$adv.season <- 1

# Use of interim year index in the advice rule
advice.ctrl_iny[[1]]$inyr.idx <- TRUE
```

* Full-population advice calendar

```{r fpaset}
# Variables to be changed for specific calendar
indices_fpa     <- indices
advice_fpa      <- advice
obs.ctrl_fpa    <- obs.ctrl
advice.ctrl_fpa <- advice.ctrl

# Observed indices at the beggining of TAC year
# (but assigned to interim year, as it is the expected B1+ at the begining of the following year)
indices_fpa[[1]][[1]] <- # lag one year back
  window(indices[[1]][[1]], 
         start=indices[[1]][[1]]@range[["minyear"]]-1, 
         end=indices[[1]][[1]]@range[["maxyear"]]-1)

indices_fpa[[1]][[1]]@range[c("startf","endf")] <- 0 # B1+ in 1st January

idxyr  <- dimnames(indices_fpa[[1]][[1]])$year
idxyr1 <- as.character(as.numeric(idxyr)+1)
yr     <- idxyr[idxyr %in% hist.yrs[-length(hist.yrs)]]
yr1    <- as.character(as.numeric(yr)+1)
idxss  <- 1

a1plus <- which(dimnames(biols[[1]]@n)$age=='1'):dim(biols[[1]]@n)[1]

indices_fpa[[1]][[1]]@index.q[,idxyr,] <- indices[[1]][[1]]@index.q[,idxyr1,] # lag one year back

indices_fpa[[1]][[1]]@index[,yr,] <- indices_fpa[[1]][[1]]@index.q[,yr,] * 
  quantSums( n(biols[[1]])[a1plus,yr1,,idxss,] * 
               wt(biols[[1]])[a1plus,yr1,,idxss,]) # B1+ at the beginning of the following year

indices_fpa[[1]][[1]]@index[,ac(proj.yrs-1),] <- NA

# Required 1st TAC
advice_fpa$TAC[,main.ctrl$sim.years["initial"],] <- 
  seasonSums(quantSums(catchWStock(fleets, stock = stkn)))[,main.ctrl$sim.years["initial"],]

# Observation of indices up to TAC year
obs.ctrl_fpa[[1]]$obs.curryr  <- TRUE
obs.ctrl_fpa[[1]]$indObs[[1]]$indObs.model <- "bio1plusfwdInd"

# Use of interim year index in the advice rule
advice.ctrl_fpa[[1]]$inyr.idx <- TRUE
```



$~$

# Run FLBEIA

We run FLBEIA in the 3 scenarios defined, which only differ in the management calendar used.
The rest of the settings are the same.

```{r runs, results='hide'}
# Interim year calendar

proj.int <- FLBEIA(biols, SRs, BDs = NULL, fleets, covars = NULL,
                   indices = indices_int, advice = advice_int,
                   main.ctrl, biols.ctrl, fleets.ctrl,
                   covars.ctrl = NULL, obs.ctrl_int, assess.ctrl, advice.ctrl)

# In-year calendar

proj.iny <- FLBEIA(biols, SRs, BDs = NULL, fleets, covars = NULL,
                   indices = indices_iny, advice = advice_iny,
                   main.ctrl, biols.ctrl, fleets.ctrl_iny,
                   covars.ctrl = NULL, obs.ctrl_iny, assess.ctrl, advice.ctrl_iny)

# Full-population advice calendar

proj.fpa <- FLBEIA(biols, SRs, BDs = NULL, fleets, covars = NULL,
                   indices = indices_fpa, advice = advice_fpa,
                   main.ctrl, biols.ctrl, fleets.ctrl,
                   covars.ctrl = NULL, obs.ctrl_fpa, assess.ctrl, advice.ctrl_fpa)
```

And 3 additional scenarios, with the same calendars, but using the rule with a biomass safeguard, 
using for example as reference level the minimum of observed indices in the historical period.

```{r runsbsafe, results='hide'}
# Interim year calendar with biomass safeguard (Imin)

advice.ctrl_intbsafe <- advice.ctrl
advice.ctrl_intbsafe[[1]]$Itrigger <- 
  apply(indices_int[[1]]$B1plusIdx@index, c(1,3:6), min, na.rm=TRUE)[drop=TRUE]

proj.int_Imin <- FLBEIA(biols, SRs, BDs = NULL, fleets, covars = NULL,
                   indices = indices_int, advice = advice_int,
                   main.ctrl, biols.ctrl, fleets.ctrl,
                   covars.ctrl = NULL, obs.ctrl_int, assess.ctrl, advice.ctrl_intbsafe)

# In-year calendar with biomass safeguard (Imin)

advice.ctrl_inybsafe <- advice.ctrl_iny
advice.ctrl_inybsafe[[1]]$Itrigger <- 
  apply(indices_iny[[1]]$B1plusIdx@index, c(1,3:6), min, na.rm=TRUE)[drop=TRUE]

proj.iny_Imin <- FLBEIA(biols, SRs, BDs = NULL, fleets, covars = NULL,
                   indices = indices_iny, advice = advice_iny,
                   main.ctrl, biols.ctrl, fleets.ctrl_iny,
                   covars.ctrl = NULL, obs.ctrl_iny, assess.ctrl, advice.ctrl_inybsafe)

# Full-population advice calendar with biomass safeguard (Imin)

advice.ctrl_fpabsafe <- advice.ctrl_fpa
advice.ctrl_fpabsafe[[1]]$Itrigger <- 
  apply(indices_fpa[[1]]$B1plusIdx@index, c(1,3:6), min, na.rm=TRUE)[drop=TRUE]

proj.fpa_Imin <- FLBEIA(biols, SRs, BDs = NULL, fleets, covars = NULL,
                   indices = indices_fpa, advice = advice_fpa,
                   main.ctrl, biols.ctrl, fleets.ctrl,
                   covars.ctrl = NULL, obs.ctrl_fpa, assess.ctrl, advice.ctrl_fpabsafe)
```

Next, we summarise all the results.

```{r sumstats}
scenarios <- c("int", "iny", "fpa", "int_Imin", "iny_Imin", "fpa_Imin")

# In this case same reference points for all the scenarios and iterations
stk.brp <- extractBRP(advice.ctrl, stkn)

biosem <- bio <- adv <- NULL

# summary statistics by iteration

for(sc in scenarios){
  
  res_sc <- get(paste0("proj.",sc))
  biosem <- rbind( biosem, bioSum(res_sc, scenario = sc, byyear=F))   # by semester
  bio    <- rbind( bio, bioSum(res_sc, scenario = sc, ssb_season = 2, brp = stk.brp))
  adv    <- rbind( adv, advSum(res_sc, scenario = sc))
}

# quantiles

bioQ <- bioSumQ(bio)
advQ <- advSumQ(adv)
```

To end with, we can generate several plot for comparing results.

```{r plots}
# reshape to long format for ggplot
dat_bioQ <- bioQ %>%
  gather("var_q", "value", -c("scenario","stock","year")) %>%
  separate("var_q", into = c("indicator", "quantile"), sep = "_") %>%
  tidyr::spread("quantile", "value")

dat_bioQ <- dat_bioQ %>% 
  mutate(calendar = unlist(strsplit(scenario, split = "_"))[1], 
         bsafe = case_when( scenario %in% c("int", "iny", "fpa") ~ "none",
                            TRUE ~ unlist(strsplit(scenario, split = "_"))[2]), 
         bsafe = factor(bsafe, levels = c("none","Imin")))

dat_brefs <- dat_bioQ %>% 
  filter(year == 1 & indicator %in% c("ssb","catch")) %>% 
  select(scenario, stock, calendar, bsafe)

dat_brefs <- bind_rows( data.frame(dat_brefs, indicator = "ssb", 
                                   Bref = c(Blim=advice.ctrl[[1]]$ref.pts["Blim",1])), 
                        data.frame(dat_brefs, indicator = "ssb", 
                                   Bref = c(Bcollapse=0.1*advice.ctrl[[1]]$ref.pts["B0",1])), 
                        data.frame(dat_brefs, indicator = "catch", 
                                   Bref = c(MSY=advice.ctrl[[1]]$ref.pts["MSY",1])))

# trajectories: biomass and catch
ggplot(data = dat_bioQ %>% filter(indicator %in% c("ssb","catch")), 
              aes(x = year, y = q50, group = calendar, col = calendar, fill = calendar)) + 
  geom_line() + 
  geom_ribbon(aes(x = year, ymin = q05, ymax = q95), alpha = 0.35) + 
  facet_grid(indicator ~ bsafe) + 
  expand_limits(y=0) +
  geom_vline(xintercept = proj.yr - 0.5, linetype = "longdash") +
  geom_hline(lty = 2, data = dat_brefs, 
             aes(yintercept = Bref), color = rep(c("orange", "red", "green"),each=nrow(dat_brefs)/3)) +
  theme_bw() + 
  theme(text = element_text(size = 20), 
        title = element_text(size = 16, face = "bold"), 
        strip.text = element_text(size = 20)) + 
  ylab("tonnes") + 
  theme(plot.title = element_text(hjust = 0.5))
```


<!-- Plot the biomass for the three scenarios. The best results are obtained with the little HCR. In this scenario the uncertainty is quite high but none of the interations fall below Blim (horizontal black line). At the start of the simulation the biomass increase but then it decrease again. The trend in msy scenario is similar to the trend in the dls3 scenario but the biomass level is lower. In fact, in this scenario the probability of being below Blim is higher than 50\%. In the dls3 scenario the biomass creased below blim in the initial part of the simulations and remained quite constant in the whole period. -->


<!-- ```{r,  fig.height= 10/2.54, fig.width =  16/2.54, fig.cap = 'Biomass time series obtained in each of the scenarios. The shaded area correspond with the 90% confidence interval.'} -->
<!-- id <- 'biomass' -->

<!-- p <- ggplot(subset(bioQ, indicator == id), aes(x=year, y=q50, ymin=q05, ymax=q95, group = scenario)) + -->
<!--   geom_ribbon(aes(fill = scenario, alpha=0.3)) + -->
<!--   geom_line(aes(color=scenario), lwd = 1) + -->
<!--   ggtitle(id) + -->
<!--   scale_y_continuous(name="tonnes") + -->
<!--   geom_hline(yintercept= advice.ctrl.msy$mur$ref.pts['Btrigger',1]) -->
<!-- print(p) -->
<!-- ``` -->


<!-- The catch in the dls3 scenario decreases year by year in the whole simulation. In the other two sencearios afeter a period of zero catch the catch increased sharply. In the case of msy scenario decreased sharply again in the final years of the simulation. In the little scenario in remained constant. -->


<!-- ```{r,  fig.height= 10/2.54, fig.width =  16/2.54, fig.cap = 'Catch time series obtained in each of the scenarios. The shaded area correspond with the 90% confidence interval.'} -->
<!-- id <- 'catch' -->

<!-- p <- ggplot(subset(bioQ, indicator == id), aes(x=year, y=q50, ymin=q05, ymax=q95, group = scenario)) + -->
<!--   geom_ribbon(aes(fill = scenario, alpha=0.3)) + -->
<!--   geom_line(aes(color=scenario), lwd = 1) + -->
<!--   ggtitle(id) + -->
<!--   scale_y_continuous(name="tonnes") + -->
<!--   geom_hline(yintercept= advice.ctrl.msy$mur$ref.pts['Btrigger',1]) -->
<!-- print(p) -->
<!-- ``` -->

<!-- ```{r,  fig.height= 10/2.54, fig.width =  16/2.54, fig.cap = 'TAC advice time series obtained in each of the scenarios. The shaded area correspond with the 90% confidence interval.'} -->
<!-- id <- 'tac' -->

<!-- p <- ggplot(subset(advQ, indicator == id), aes(x=year, y=q50, ymin=q05, ymax=q95, group = scenario)) + -->
<!--   geom_ribbon(aes(fill = scenario, alpha=0.3)) + -->
<!--   geom_line(aes(color=scenario), lwd = 1) + -->
<!--   ggtitle(id) + -->
<!--   scale_y_continuous(name="tonnes") + -->
<!--   geom_hline(yintercept= advice.ctrl.msy$mur$ref.pts['Btrigger',1]) -->
<!-- print(p) -->
<!-- ``` -->



$~$

# Shiny App

Another FLBEIA feature can be used to generate a Shiny Web App with the simulation results for all 
the scenarios tested.


<!-- ```{r shiny} -->
<!-- # Reference points (same for all scenarios and iterations) -->
<!-- refptn <- c("Bmsy","Fmsy", "Bpa", "Blim", "Fpa", "Flim") -->
<!-- refpt0 <- dimnames(advice.ctrl[[1]]$ref.pts)[[1]] -->
<!-- RefPts <- data.frame(stock = stkn, scenario = rep(scenarios, each = length(refptn)),  -->
<!--                      refpoint = rep(refptn, length(scenarios)), -->
<!--                      value = NA) -->

<!-- for (rp in refptn[refptn %in% refpt0]) -->
<!--   RefPts[RefPts$refpoint == rp, "value"] <- advice.ctrl[[1]]$ref.pts[rp,1] -->


<!-- bio_long <- bioQ -->

<!-- # dat_bioQ <- dat_bioQ %>% -->
<!-- #   gather("var_q", "value", -c("scenario","stock","year",OM:BSAFE,"SCnam","LHnam")) %>% -->
<!-- #   separate("var_q", into = c("indicator", "quantile"), sep = "_") %>% -->
<!-- #   tidyr::spread("quantile", "value") -->


<!-- flbeiaApp(RefPts = RefPts, bio = bioQ, adv = advQ)#,  -->
<!--           # fltStk = fltStkQ, mt = mtQ, mtStk = mtStkQ, risk = risk, -->
<!--           # years = as.character(2010:2024),  -->
<!--           # calculate_npv = FALSE, npv =  NULL, npv.y0 = NULL, npv.yrs = NULL)  -->


<!-- ``` -->



$~$

# Acknowledgments

This tutorial has been built with the financial support of the Dirección de Pesca del Gobierno Vasco 
(Basque Government - Spain). 

Special thanks also to the members of the ICES Working Workshop on Data-limited Stocks of Short-lived Species
(WKDLSSLS) for their usefull discussions.



$~$

# More information

* You can submit bug reports, questions or suggestions on this tutorial at <https://github.com/flr/doc/issues>.
* Alternatively, send a pull request to <https://github.com/flr/doc/>.
* For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage: <http://flr-project.org>.
* You can submit bug reports, questions or suggestions specific to **FLBEIA** to <flbeia@azti.es>.


$~$

## Software Versions

* `r version$version.string`
* FLCore: `r packageVersion('FLCore')`
* FLBEIA: `r packageVersion('FLBEIA')`
* spict: `r packageVersion('spict')`
* fishmethods: `r packageVersion('fishmethods')`
* **Compiled**: `r date()`


$~$

## License

This document is licensed under the [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0) license.

This code forms part of the work for the publication Sánchez 




$~$

## Author information

**Sonia Sanchez-Maroño**. AZTI, Marine Research Unit. Herrera Kaia, Portualdea z/g, 20110, Pasaia, Gipuzkoa, Spain. https://www.azti.es/.

**FLBEIA team**. AZTI. Marine Reserach Unit. Txatxarramendi Ugartea z/g, 48395, Sukarrieta, Basque Country, Spain.
**Mail** flbeia@azti.es