Tutorial_2.Rmd

---
title: "Rfishpop (Second tutorial)"
author: "Marta Cousido Rocha, Santiago Cerviño López, Maria Grazia Pennino"
date: "29/7/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Population dynamic tools in support of fisheries managment

After describing the operating model used to generate "true" ecosystem dynamics including the
natural variations in the system. It is important to address this tutorial which focus on  a set of methods contained in the package to estimate Maximum Sustainable Yield (MSY) reference
points. These allow to identify management targets in terms of fishing intensity, population status and yield.


To estimate these reference points several fuctions are required: biomass-per-Recruit, yield-per-Recruit,  total Yield, Biomass and Recruitment in equilibrium, and Reference Points.

To finish the tutorial a  description and an application of the Distribution.length function  which allow us to compute the stock length or catches length  of the population simulated previously using the OM are provided.


### <span style="color:blue">Biomass-per-Recruit (BPR function)</span> 


The function *BPR* returns biomass-per-Recruit for each iteration. The arguments of the function are the following.  


- *Pop.Mod*	is a list containing the components  returned by Population.Modeling function (main function).

- *f.grid*	is a sequence of *f's* which are the annual components of fishing mortality $F = f * SEL$.

- *Fish.years*	is the number of recent years to estimate the mean of SEL (selectivity, see information about such element in Sum.Pop.Mod function).

- *Bio.years*	is the number of recent years to estimate the mean of M, W and Mat (natural mortality, weight and maturity, see information about such elements in Population.Modeling function).

- *plot*	is a vector of two elements. The first one is a logical parameter. By default is equal to TRUE, which means that a biomass per recruit graph is done. The second element refers to which iteration must be plotted.

- *Method*	is the procedure to obtain the age vector of weight, natural mortality, selectivity and maturity for each iteration. By default is "mean" which means that the mean of the last "Bio.years or "Fish.years" is used. The alternative option is "own", the user can introduce these vectors in a matrix.

- *par*	If Method="own" it is a list containing the matrices whose columns report for each iteration the age vector of weight, natural mortality, selectivity and maturity. In other case is equal to NULL.

The result is an array whose third dimension corresponds to the iterations. For each iteration the array contains a matrix reporting the biomass-per-recruit for a range of overall fishing mortalities.

First of all we need to simulate our population using Population.Modeling function. 



```{r}
library(Rfishpop)
ctrPop<-list(years=seq(1980,2020,by=1),niter=2,N0=15000,ages=0:15,minFage=2,
maxFage=5,tc=0.5,seed=NULL)
number_ages<-length(ctrPop$ages);number_years<-length(ctrPop$years)

Mvec=c(1,0.6,0.5,0.4,0.35,0.35,0.3,rep(0.3,9))
M<-matrix(rep(Mvec,number_years),ncol = number_years)
colnames(M)<-ctrPop$years
rownames(M)<-ctrPop$ages

ctrBio<-list(M=M,CV_M=0.2, L_inf=20, t0=-0.25, k=0.3, CV_L=0, CV_LC=0, a=6*10^(-6), b=3,
           a50_Mat=1, ad_Mat=-0.5,CV_Mat=0)

ctrSEL<-list(type="cte", par=list(cte=0.5),CV_SEL=0)

f=matrix(rep(0.5,number_years),ncol=number_years,nrow=2,byrow=TRUE)

ctrFish<-list(f=f,ctrSEL=ctrSEL)

a_BH=15000; b_BH=50; CV_REC_BH=0

SR<-list(type="BH",par=c(a_BH,b_BH,CV_REC_BH))

Pop.Mod<-Population.Modeling(ctrPop=ctrPop,ctrBio=ctrBio,ctrFish=ctrFish,SR=SR)
```

Now, we can use the *BPR* function to obtain the  biomass-per-recruit for a range of overall fishing mortalities.


```{r}
f.grid<-seq(0.00,0.5,by=0.01)
bpr<-BPR(Pop.Mod,f.grid,Bio.years=3,Fish.years=3,plot=c(TRUE,1),Method="mean",par=NULL)
head(bpr[,,1])
```

Note that f.grid is a sequence of fishing efforts that the user must be defined. In this case we have introduced a grid from 0 to 0.5 such that two consecutive points are separated by 0.01.

If we want to use *Method="own"* we need to specify argument *par*. Assuming that W, M, Mat and SEL are equal to their values in the last year of the population. We have the following example.
 
```{r}

W=Pop.Mod$Matrices$W[,41,]

SEL=Sum.Pop.Mod(Pop.Mod,"SEL" )$SEL[,41,]

Mat=Pop.Mod$Matrices$Mat[,41,]

M=Pop.Mod$Matrices$M[,41,]

par=list(); par$W<-W; par$SEL<-SEL; par$Mat<-Mat; par$M<-M

bpr=BPR(Pop.Mod,f.grid,plot=c(TRUE,1),Method="own",par=par)
head(bpr[,,1])
head(bpr[,,2])
```

### <span style="color:blue">Yield-per-Recruit (YPR function) </span> 



The function *YPR* returns yield-per-Recruit for each iteration. Almost all of the arguments of this function have been explained above in *BPR* function. Except the following argument which must be described again due to a slightly adjustment.

- *Bio.years* is the number of recent years to estimate the mean of M and WC (natural mortality and catches weight, see information about such elements in Population.Modeling function and Sum.Pop.Mod function, respectively).

The result is an array whose third dimension corresponds to the iterations. For each iteration the array contains a matrix reporting the yield-per-recruit for a range of overall fishing mortalities.


The following lines of code allow  us to compute the yield-per-Recruit for each iteration. 


```{r}
ypr<-YPR(Pop.Mod,f.grid,3,3,plot=c(TRUE,1), Method="mean",par=NULL)
head(ypr[,,1])
```


If we want to use *Method="own"* the same example explained above is suitable.
 
```{r}
WC=Sum.Pop.Mod(Pop.Mod,"WC" )$WC[,41,]

SEL=Sum.Pop.Mod(Pop.Mod,"SEL" )$SEL[,41,]

M=Pop.Mod$Matrices$M[,41,]

par=list(); par$WC<-WC; par$SEL<-SEL; par$M<-M
ypr=YPR(Pop.Mod,f.grid,plot=c(TRUE,1),Method="own",par=par)
head(ypr[,,1])
head(ypr[,,2])
```

### <span style="color:blue">Total Yield, Biomass and Recruiment in Equilibrium (BYR.eq function)</span> 


The BYR.eq function returns total Yield, Biomass and Recruiment in equilibrium for the corresponding fishing efforts and fishing mortalities. For each case the function also returns the corresponding biomass-per-recruit and yield-per-recruit. Furthermore, the function also reports the population size in equilibrium for each age and fishing mortality.


The arguments of this function have been explained above in the YPR and BPR functions. The argument *Bio.years* is adjusted according to its use in the BYR.eq function.

- *Bio.years*	is the number of recent years to estimate the mean of M, Mat, WC, and W (natural mortality, maturity, stock weight and catches weight, see information about such elements in Population.Modeling function and Sum.Pop.Mod function).


The result is a list of two elements:

- *N* 	an array whose third dimension is the number of iterations. For each iteration it reports a matrix containing the population size in equilibrium for each age.

- *DPM*	an array whose third dimension is the number of iterations. For each iteration it reports a matrix containing the total biomass in equilibrium, total yield in equilibrium, total recruitment in equilibrium, biomass-per-recruit, and yield-per-recruit for a range of overall fishing mortalities.



The following lines of code show how to use the BYR.eq function.


```{r}
RE<-BYR.eq(Pop.Mod,f.grid,3,3,FALSE,Method="mean",par=NULL)
N_eq<-RE$N
DPM<-RE$DPM
### First iteration
head(DPM[,,1])
### Second iteration
head(DPM[,,2])

```

If we want to use *Method="own"* the same example explained above is suitable.

```{r}
par=list(); par$W<-W; par$SEL<-SEL; par$Mat<-Mat; par$M<-M; par$WC=WC
RE=BYR.eq(Pop.Mod,f.grid,plot=FALSE,Method="own",par=par)
N_eq<-RE$N
DPM<-RE$DPM
### First iteration
head(DPM[,,1])
### Second iteration
head(DPM[,,2])

```

### <span style="color:blue">Reference Points (RF function)</span> 


The function *RF* returns the reference fishery mortality which produces maximum YPR
(FM_type="F_max"), the reference fishery mortality at which the slope of the YPR curve is reduced to 0.1 of that estimated at the origin (FM_type="F_0.1"), the reference fishery mortality at which the MSY is attained (FM_type="F_msy"), the reference fishery mortality which will drive the stock to extinction (FM_type="F_Crash"), and the reference fishery mortality at which the BPR is equal to prop multiplied by BPR0 where prop is any proportion that the user can fix and BPR0 is the BPR at virgin biomass (FM_type="F_BPR").
Furthermore for each of these fishery mortalities the corresponding effort of fishing, YPR, BYR, B, Y and R in equilibrium is also returned. 

It is important to mention that if the fishing effort is equal to 10 can be that the optimize process had not found the correct value in the default sequence (except for the reference fishery mortality which will drive the stock to extinction, in this case the sequence finishes at 60).


Some arguments coincide with those of BPR, YPR and BYR.eq functions, to avoid redundancy we focus on the new arguments. 

- *FM_type*	a vector containing which of the five reference fishery mortalities must be computed. The possibilities have been described above:  FM_type="F_max", FM_type="F_0.1", FM_type="F_msy", FM_type="F_Crash" and FM_type="F_BPR".

- *iters* is a  vector containing the iterations for which the reference fishery mortalities must be computed.


- *plot* is a logical parameter equal to TRUE if the user desires to obtain several informative plots (a deep explanation above).

- *prop*	is the proportion of BPR at virgin biomass at which FM_type="F_BPR" refers. By default is 0.3. If the user desires the result for different proportions a vector containing the different values must be introduced using prop argument.


The function returns a list reporting arrays (third dimension *iters*) containing the corresponding elements of the following list depending the above selection (*FM_type* argument).


- *F_max*	is the value of F that produces maximum YPR with the corresponding effort of fishing, YPR, BYR, B, Y and R in equilibrium.

- *F_0.1* is	the value of F at which the slope of the YPR curve is reduced to 0.1 of that estimated at the origin with the corresponding effort of fishing, YPR, BYR, B, Y and R in equilibrium.

- *F_msy* is  the value of F at which the MSY is attained with the corresponding effort of fishing, YPR, BYR, B, Y and R in equilibrium.

- *F_Crash*	is the value of F which will drive the stock to extinction with the corresponding effort of fishing, YPR, BYR, B, Y and R in equilibrium.


- *F_BPR* is the value of F at which at which the BPR is equal to prop multiplied by BPR0, where prob is any proportion that the user can fix and BPR0 is the BPR at virgin biomass, with the corresponding effort of fishing, YPR, BYR, B, Y and R in equilibrium. To access to this element use F_BPRprop following by symbol of percentage (without spaces and between quotation marks) where prop is the value that the user specifies. If prop vector length is greater than 1 the function will return several slots in the list reporting the corresponding values for each value of prop. The user can access to the elements as we mentioned before replacing prop for any of the value in such vector.

It is worth to mention that if plot=TRUE, the function reports the plots: equilibrium biomass v. fishing mortality, equilibrium yield v. fishing mortality, equilibrium recruitment v. biomass, and equilibrium yield v. biomass. All the plots are done using the information of the first iteration in the vector of parameter *iters*. The corresponding reference fishery mortalities are also plotted in the previous curves. Note that only the fishery mortalities required by the argument FM_type are plotted. If prop is a vector of length greater than 1 only the values corresponding to the first element of the vector will be represented in the plot for simplicity. On the other hand, note that F_msy and F_max coincide when the recruitment relationship is constant and hence only one of both appears in the plot. The reference fishery mortalities F_max and F_Crash can be overlapped in equilibrium recruitment v. biomass and equilibrium yield v. biomass plots then only one can be shown in the plot.

Note that F_Crash does not exist when the recruitment model is constant.

Below, we provide some examples of the use of RF function.

1. We compute the fishery mortality which produces maximum YPR for the two iterations of the Pop.Mod object and we plot the first iteration in *iters* parameter.

```{r}
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type="F_max",iters=1:2,plot=TRUE)
```

2.  We compute the fishery mortality at which the slope of the YPR curve is reduced to 0.1 of that estimated at the origin for the two iterations of the Pop.Mod object and we plot the first iteration in *iters* parameter

```{r}
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type="F_0.1",iters=1:2,plot=TRUE)
```


3. We compute the fishery mortality at which the MSY is attained  for the two iterations of the Pop.Mod object and we plot the first iteration in *iters* parameter. Using the second commented line it is possible to compute such information only for the second iteration (note that in such case the plot refers to the second iteration).

```{r}
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type="F_msy",iters=1:2,plot=TRUE)
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type="F_msy",iters=2,plot=TRUE)
```

4. We compute the reference fishery mortality which will drive the stock to extinction for the two iterations of the Pop.Mod object and we plot the first iteration in *iters* parameter.

```{r}
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type="F_Crash",iters=1:2,plot=TRUE)
```

5. We compute the reference fishery mortality at which the BPR is equal to prop=0.4 multiplied by BPR0 where BPR0 is the BPR at virgin biomass for the two iterations of the Pop.Mod object and we plot the first iteration in *iters* parameter.

```{r}
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type="F_BPR",iters=1:2,plot=TRUE,prop=0.4)
```

Below, we provide some examples where several fishery mortalities are computed at the same time.

```{r}
RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type=c("F_Crash","F_msy"),iters=1:2,plot=TRUE)

RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type=c("F_Crash","F_msy","F_0.1","F_max","F_BPR"),iters=1:2,plot=TRUE)

RF(Pop.Mod, 3,3,Method="mean",par=NULL,FM_type=c("F_Crash","F_BPR"),iters=1:2,plot=TRUE,prop = c(0.2,0.4))
```

NOTE (remember) that if we want to use *Method="own"* we need to specify argument *par*, the following example is similar to the explained previously for the other functions.
 
```{r} 
par=list(); par$W<-W; par$SEL<-SEL; par$Mat<-Mat; par$M<-M; par$WC=WC
RF.Result=RF(Pop.Mod,Method="own",par=par,FM_type=c("F_Crash","F_BPR"),iters=1:2,plot=TRUE)
RF.Result
```


- <span style="color:red"> This function is time demanding due to the necessary complex optimization process. We are working in it. While it is possible to use it for relatively small populations.</span>



### <span style="color:blue">Reference Points plots (plotRF function)</span> 

All the plots provided by the previous function RF are done using the information of the first iteration in the vector of parameter *iters*. For this reason the following function is useful, since it allows to plot the information of any of the *iters*.


The arguments of the function are the following.  

- *Pop.Mod*	 is a list containing the components returned by Population.Modeling function (main function).

- *RF.result*	is a list containing the components returned by RF function (Reference Points).
- *iter* is the iteration for which the plot must be carried out.

Below, we can see an example of the use of this simple function.

```{r}
plotRF(Pop.Mod,RF.result=RF.Result, iter=2)
```


### <span style="color:blue">Distribution length (Distribution.length function)</span> 
The 
*Distribution.length* function allow us to explore our population providing stock or length catches distribution for each year and iteration. 

The arguments of this function are:

- *Pop.Mod* is the object returned by Population.Modeling function.


- *CV* 	is the coefficient of variation associated to the log-normal distribution (see explanation above)


- *Type* is an indicator of which distribution length must be computed, length stock distribution (Type="LengthS") or length catches distribution (Type="LengthC").

- *RF.value* is the number of values generated for each age (given a year and an iteration) from the log-normal distribution (see the explanation above). By default RF.value=1000.

The function returns the stochastic length distribution of the stock (Type="LengthS") or length catches distribution (Type="LengthC") for each year and iteration. In the case of the stock length distribution it is computed generating for each age, year and iteration RF.value random values from a log-normal distribution centered in the corresponding stock length and whose variability comes from the given CV. For the catches length distribution the mean of the log-normal distribution is given by the corresponding catch length. For the stock length distribution the distribution obtain for each age (given a year and an iteration) is scaled using the corresponding stock number (N matrix), whereas in the catch distribution this role is for the catch matrix C.

The result is an array whose third dimension is the number of iterations, and the second one is the different years. Hence each column contains the distribution length (stock or catches) for each year.

We obtain the stock distribution length as follows.


```{r}
L.D<-Distribution.length(Pop.Mod,CV=0.2,Type="LengthS")
```


We can check the distribution length for each iteration and year. For example, below we explore iteration 1 and year 1980.

```{r}
L.D[,,1][,1]
plot(L.D[,,1][,1], type="b", pch=19, col="red", xlab="", ylab="",main = "Distribution of stock length year 1980 iteration 1")
LS<-Sum.Pop.Mod(Pop.Mod,c("LS"))
```

Using the above lines of code only changing Type="LengthS" by Type="LengthC" we obtain the distribution of capture length.

```{r}
L.D<-Distribution.length(Pop.Mod,CV=0.2,Type="LengthC")
```
```{r}
L.D[,,1][,1]
plot(L.D[,,1][,1], type="b", pch=19, col="red", xlab="", ylab="",main = "Distribution of capture length year 1980 iteration 1")
LC<-Sum.Pop.Mod(Pop.Mod,c("LC"))
```