a_seasonal_operating_model.Rmd

---
title: "Seasonal Operating Model"
output:
github_document:
  mathjax: TRUE
  html_document:
  fig_width: 6 
  fig_height: 4 
  toc: yes
bibliography: bibliography.bib
---

```{r knitr_init, echo=FALSE, results="hide"}
library(knitr)   
## Global options
opts_chunk$set(echo      =FALSE, 
               eval      =TRUE,
               prompt    =FALSE,
               comment   =NA,
               message   =FALSE,
               warning   =FALSE,
               tidy      =TRUE,
               fig.height=6,
               fig.width =8,
               cache     =TRUE,
               cache.path="../cache/class/seasonal/")

options(digits=3)

iFig=0
```

This vignette shows how to condition a seasonal Operating Model, i.e. one with less than annual time steps that are less than annual, based on life-history theory.

[Packages](#Packages)

[Operating Model Conditioning](#OM%5D)

[More Information](#More)

[References](#References)

# Packages {#Packages}

Some 'FLR' packages are required

```{r, pkgs3, echo=TRUE}
library(FLCore)
library(ggplotFL)

library(FLBRP)
library(FLasher)

library(mydas)

library(FLRef)
```

as well as `R` packages for plotting

```{r, pkgs1, echo=TRUE}
library(ggplot2)
library(GGally)
library(ggpubr)
```

and data manipulation

```{r, pkgs2, echo=TRUE}
library(reshape)
library(plyr)
library(dplyr)
library(reshape)
```

## Quick Start {#QuickStart}

The `FLife` packge is required for modelling life-histories, and 'mydas' for estimation and developing harvest control rules

The simplest way to obtain the `FLR` pacakges are to install them from the `FLR` repository

```{r install,echo=TRUE,eval=FALSE}
install.packages("FLife", repos = "http://flr-project.org/R")
```

See `help(install.packages)` for more details.

After installation, they can be loaded

```{r lib,echo=TRUE}
library(FLife)
library(mydas)
```

[Back to Top](#top)

# Operating Model Conditioning

## Annual model

Conditioned an Operating Models on life-history parameters allows more flexibility than using a data-rich stock assessment, that may never have been formally validated.

There are relationships between the life-history parameters, e.g. the von Bertalanfy growth parameters $L_{\infty}$, $k$ and $t_0$, length at which $50\%$ of individuals reach maturity ($l_{50}$), and of the length weight relationship $W=aL^b$.

There is an example data set with `teleost` life-history parameters in `FLife`

```{r}
data(teleost)

teleost
```

```{r}
ggpairs(model.frame(teleost)[,c("linf","k","l50","a","b")])+
  scale_x_log10()+
  scale_y_log10()
```

**Figure** `r iFig=iFig+1; iFig` relationships between life-history parameters

If only the maximum size of a species is known then you can derive the other parameters from the life-history relationships.

```{r, param, echo=TRUE}
par=lhPar(FLPar(c(linf= 59.1),units="NA"))

par
```

### Equilibrium dynamics

The life-history parameters can be used to create the equilibrium dynamics using the `FLBRP` class.

```{r, eq, echo=TRUE}
eq=lhEql(par)

plot(eq,refpts="msy")
```

**Figure `` `r iFig+1; iFig` ``** Equilibrium dynamics with $MSY$ reference points.

### Time series dynamics

For the Operating Model, the equilibrium `FLBRP` object needs to be converted into an `FLStock` object.

As an example, time series are simulated for fishing mortality, where an unexploited stock is fished until it becomes over exploited, after which a recovery plan is implemented.

```{r, om-3, echo=TRUE}
fbar(eq)=refpts(eq)["msy","harvest"]%*%FLQuant(c(rep(0.1,60),
                                               seq(0.1,2,length.out=40)[-40],
                                               seq(2,0.7,length.out=11),rep(0.7,20)))[,-1]
```

Create the `FLStock`

```{r, om-1, echo=TRUE}
om=as(eq,"FLStock")
```

Turn it into a stochastic object so that Monte Carlo simulations can be performed

```{r, om-2, echo=TRUE}
om=propagate(om,100)
```

Then perform a stochastic projection

```{r, om, echo=TRUE}
set.seed(234)
srDev =rlnoise(100,fbar(eq)%=%0,0.3,0)

## project for a target fishing mortality
f=fbar(eq)[,-1]

om=ffwd(om,fbar=f,sr=eq,deviances=srDev)

## get rid of burn in period
om=window(om,start=41)
```

```{r, fig-om, fig.width=10,fig.height=6, caption="Annual Projection", echo=TRUE}
mets=list(RMSY =function(x, y) rec(  x)/refpts(y)["msy","rec"],
          YMSY =function(x, y) catch(x)/refpts(y)["msy","yield"],
          SBMSY=function(x, y) ssb(  x)/ sbmsy(y),
          FMSY =function(x, y) fbar( x)/  fmsy(y))

fqs=FLQuants(lapply(mets, function(m) m(om, eq)))

plot(fqs,iter=1)+ 
  ylim(c(0, NA))+
  geom_hline(yintercept=1, linetype=2)+
  scale_x_continuous(limits=c(50,120))+
  theme_bw(16)+
  theme(legend.position="none")
```

**Figure `r iFig=iFig+1; iFig`** Operating Model with $MSY$ reference points.

## Seasonal Operating Model

```{r, seasonalise}
source("~/Desktop/flr/FLCandy/R/seasonalise.R")
seasonalise<-function(object, season=1:4){
  
  ## Stock and recruit                                   ###
  ## set expected to 1 and model variability as deviates ###
  sr=as.FLSR(object,model="geomean")
  params(sr)=FLPar(1,dimnames=list(params="a",iter=1))
  
  recs=FLCore:::expand(rec(object),season=season)
  
  ## Add seasons                                         ###
  object=FLCore:::expand(object,season=season)
  
  ## Divide up mortality by season                       ###
  m(      object)=m(object)/dim(object)[4]
  harvest(object)=harvest(object)/dim(object)[4]
  
  ## Seasonal growth                                     ###
  stock.wt(   object)[,-dim(stock.wt(   object))[2]]=wtInterp(stock.wt(   object))
  catch.wt(   object)[,-dim(catch.wt(   object))[2]]=wtInterp(catch.wt(   object))
  landings.wt(object)[,-dim(landings.wt(object))[2]]=wtInterp(landings.wt(object))
  discards.wt(object)[,-dim(discards.wt(object))[2]]=wtInterp(discards.wt(object))

  object=adjust(object)

  ## Project for historic F                              ###
  #fbar=as(FLQuants("fbar"=fbar(object)[,-1]),"fwdControl")
  #object=fwd(object,control=fbar,sr=sr,residuals=recs)
  
  object}

wtInterp<-function(wt){
  
  mlt=wt[,-dim(wt)[2]]
  mlt=FLQuant(rep(seq(0,(dim(mlt)[4]-1)/dim(mlt)[4],1/dim(mlt)[4]),each=max(cumprod(dim(mlt)[1:3]))),
              dimnames=dimnames(mlt))[-dim(mlt)[1]]
  
  incmt=wt[,,,1]
  incmt=-incmt[-dim(incmt)[1],-dim(incmt)[2]]+incmt[-1,-1]
  
  wt[dimnames(incmt)$age,dimnames(incmt)$year]=
    wt[dimnames(incmt)$age,dimnames(incmt)$year]+
    incmt%+%(mlt%*%incmt)
  
  wt[,-dim(wt)[2]]}

annualise<-function(x) {
  
  # ADD slots (catch/landings.discards, m)
  res <- qapply(x, function(s) unitSums(seasonSums(s)))
  
  # MEAN wts
  catch.wt(res)[] <- unitSums(seasonSums(catch.wt(x) * catch.n(x))) / 
    unitSums(seasonSums(catch.n(x)))
  
  landings.wt(res)[] <- unitSums(seasonSums(landings.wt(x) * landings.n(x))) / 
    unitSums(seasonSums(landings.n(x)))
  
  discards.wt(res)[] <- unitSums(seasonSums(discards.wt(x) * discards.n(x))) / 
    unitSums(seasonSums(discards.n(x)))
  discards.wt(res)[is.na(discards.wt(res))]=seasonMeans(discards.wt(x))[is.na(discards.wt(res))]

  stock.wt(res)[] <- unitSums(seasonSums(stock.wt(x) * stock.n(x))) / 
    unitSums(seasonSums(stock.n(x)))
  
  
  # RECONSTRUCT N: N0 = N1 / exp(-M - F)
  
  stkn <- unitSums(stock.n(x)[,,,4])
  m(res) <- unitMeans(seasonSums(m(x)))
  harvest(res) <- unitMeans(seasonSums(harvest(x)))
  
  stock.n(res)[] <- stkn / exp(-m(res) - harvest(res))
  
  # mat
  
  mat(res)[] <- unitSums(seasonSums(mat(x) * stock.n(x))) / 
    unitSums(seasonSums(stock.n(x)))
  mat(res) <- mat(res) %/% apply(mat(res), c(1,3:6), max)
  mat(res)[is.na(mat(res))] <- 0
  
  # spwn
  
  m.spwn(res) <- 0.5
  harvest.spwn(res) <- 0.5
  
  # totals
  
  catch(res) <- computeCatch(res)
  landings(res) <- computeLandings(res)
  discards(res) <- computeDiscards(res)
  stock(res) <- computeStock(res)
  
  return(res)}

```

Take an annual `FLStock` object and *"seasonalise"* it, i.e. interpolate growth, and split natural mortality by season.

```{r, om4, echo=TRUE}
om4=seasonalise(window(iter(om,3), end=81))
```

If spawning is not year round, as is the case in temperate species, it is necessary to assign the stock-recruit relationship parameters to the spawning season

```{r, sr, echo=TRUE}
sr.par=FLPar(NA, dimnames=list(params=dimnames(params(eq))$params, season=1:4, iter=1))
sr.par[,2]=params(eq)[,1]
sr.par[1,2]=sr.par[1,2]
sr.par[2,2]=sr.par[2,2]
sr4=predictModel(model=model(eq), params=sr.par)

sr4
```

Since the selection pattern is asymptotic, set fbar range to be oldest true age, also make years correspond to current time period.

As there is no spawning outside of $2^{nd}$ season so set these to 0

```{r, om4-tidy, echo=TRUE}
range(om4)[c("minfbar","maxfbar")]=39

dimnames(om4) <- list(year=an(dimnames(om4)$year)+1950)

mat(om4)[,,,c(1,3:4)]=0
```

Run a deterministic projection

```{r, om4-1, echo=TRUE}
control=as(FLQuants("f"=fbar(om4)[,-1]*0.5),"fwdControl")  
om4=fwd(om4,control=control,sr=sr4)

plot(window(om4,start=2010,end=2030)) 
```

**Figure `r iFig=iFig+1; iFig`** Seasonal Operating Model

Simulate a seasonal fishery

```{r, om4-2, echo=TRUE}
f=fbar(om4)[,ac(2021:2030)]*4
f[,ac(2021:2030),,1:3]=0 

control=as(FLQuants("f"=f),"fwdControl")  
om4.2=window(fwd(om4,control=control,sr=sr4),end=2030)

plot(window(FLStocks("All year round"=om4,"Seasonal\nFishery"=om4.2),start=2010,end=2030))
```

**Figure `r iFig=iFig+1; iFig`** Seasonal Operating Model and fishery

```{r, om4-3, echo=TRUE}
om=annualise(om4)
om=ffwd(om, sr=eq, fbar=fbar(om)[,-1])

plot(window(FLStocks("Annual"=om,"Seasonal"=om4),start=2010,end=2030))
```

**Figure `r iFig=iFig+1; iFig`** Convert seasonal to an annual Operating Model

[Back to Top](#top)

# More Information {#More}

-   You can submit bug reports, questions or suggestions on `FLife` at the `FLife` issue page [^1], or on the *FLR* mailing list.
-   Or send a pull request to <https://github.com/lauriekell/FLife/>
-   For more information on the FLR Project for Quantitative Fisheries Science in R, visit the FLR webpage [^2].
-   The latest version of `FLife` can always be installed using the `devtools` package, by calling

[^1]: <https://github.com/lauriekell/FLife/issues>

[^2]: <http://flr-project.org>

```{r, devtools, echo=TRUE, eval=FALSE}
	library(devtools)
	install_github("flr/FLife")
```

\`

## Software Versions

-   `r version$version.string`
-   FLCore: `r packageVersion('FLCore')`
-   FLPKG: `r # packageVersion('FLPKG')`
-   **Compiled**: `r date()`
-   **Git Hash**: `r system("git log --pretty=format:'%h' -n 1", intern=TRUE)`

## Author information

**Laurence KELL**. [laurie\@seaplusplus.co.uk](mailto:laurie@seaplusplus.co.uk){.email}

## Acknowledgements

This vignette and the methods documented in it were developed under the MyDas project funded by the Irish exchequer and EMFF 2014-2020. The overall aim of MyDas is to develop and test a range of assessment models and methods to establish Maximum Sustainable Yield (MSY) reference points (or proxy MSY reference points) across the spectrum of data-limited stocks.

# References {#References}

[Back to Top](#top)