HW4.Rmd

---
title: "Edwards-Perez-Whelpley-HW4"
output:
  word_document: default
  html_document: default
  pdf_document: default
date: "2022/11/08"
---

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

```{r,  include=FALSE}
library(Rpath)
library(data.table)
#devtools::install_github('NOAA-EDAB/Rpath', build_vignettes = TRUE)
library(devtools)
library(ggplot2)
```

**Anchovy Bay: an imaginary ecosystem model**
\vspace{10pt}

The bay covers an area of 100 km2 and is rather isolated from other marine systems, and we can assume that the populations stay in the bay year-round. Our model will have the following 16 groups. The first 10 are living: whales, seals, cod, whiting, mackerel, anchovy, shrimp, benthos, zooplankton, and phytoplankton. There is one detrital group called detritus. In addition, there are 5 fleets: sealers, trawlers, seiners, bait boats, and shrimpers

\vspace{10pt}
**Building the parameter file for Rpath:**
\vspace{10pt}

```{r,  include=FALSE}
groups <- c('whales', 'seals', 'cod', 'whiting', 'mackerel', 'anchovy', 'shrimp',
            'benthos', 'zooplankton', 'phytoplankton', 'detritus', 'sealers', 
            'trawlers', 'seiners', 'bait boats', 'shrimpers')

types <- c(rep(0, 9), 1, 2, rep(3, 5))

AB.params <- create.rpath.params(groups, types)
```

**Fisheries Data:**
\vspace{10pt}

The sealers caught 15 seals with an average weight of 30kg. The fisheries catches were 45t of cod and 20t of whiting for the trawlers, 40t of mackerel and 120t of anchovy for the seiners, 20t of anchovy for the bait boats, and 5t of shrimp for the shrimpers. Calculate catches using the appropriate unit (t/km2/year), and enter in Rpath.

*Static parameters for this example:*
Biomass accumulation and un-assimilated production
\vspace{10pt}

```{r,  include=FALSE}
AB.params$model[Group=="seals", sealers:=0.0045] 
AB.params$model[Group=="cod", trawlers:=0.45]
AB.params$model[Group=="whiting", trawlers:=0.20]
AB.params$model[Group=="mackerel", seiners:=0.40]
AB.params$model[Group=="anchovy", seiners:=1.20]
AB.params$model[Group=="anchovy", "bait boats":=0.20]
AB.params$model[Group=="shrimp", shrimpers:=0.05]
```

**Biological Data:** 
\vspace{10pt}

```{r,  include=FALSE}
AB.params$model[Group=="whales", Biomass:=0.08] 
AB.params$model[Group=="seals", Biomass:=0.0609]
AB.params$model[Group=="cod", Biomass:=3.1]
AB.params$model[Group=="whiting", Biomass:=1.7]
AB.params$model[Group=="mackerel", Biomass:=1.22]
AB.params$model[Group=="anchovy", Biomass:=6]
AB.params$model[Group=="shrimp", Biomass:=0.8]
AB.params$model[Group=="zooplankton", Biomass:=14.8]
AB.params$model[Group=="detritus", Biomass:=10]
AB.params$model[Group=="phytoplankton", Biomass:=9]
```

P/B ratio: check!

```{r, include=FALSE}
AB.params$model[Group=="whales", PB:=0.05] 
AB.params$model[Group=="seals", PB:=0.1639]
AB.params$model[Group=="cod", PB:=0.355]
AB.params$model[Group=="whiting", PB:=0.587]
AB.params$model[Group=="mackerel", PB:=1.1478]
AB.params$model[Group=="anchovy", PB:=1.2233]
AB.params$model[Group=="shrimp", PB:=3]
AB.params$model[Group=="zooplankton", PB:=35]
AB.params$model[Group=="benthos", PB:=3]
AB.params$model[Group=="phytoplankton", PB:=240]
```

Q/B ratio: check!

```{r,  include=FALSE}
AB.params$model[Group=="whales", QB:=9]
AB.params$model[Group=="seals", QB:=15]
AB.params$model[Group=="cod", QB:=2.58]
AB.params$model[Group=="whiting", QB:=3.3]
AB.params$model[Group=="mackerel", QB:=4.4]
AB.params$model[Group=="anchovy", QB:=9.13]
```

E/E: check!

```{r,  include=FALSE}
AB.params$model[Group=="benthos", EE:=0.6]
```

```{r,  include=FALSE}
AB.params$model[Group=="shrimp", ProdCons:=0.25]
AB.params$model[Group=="benthos", ProdCons:=0.25]
AB.params$model[Group=="zooplankton", ProdCons:=0.25]
```

DIET DATA: check!
\vspace{20pt}

```{r,  include=FALSE}
AB.params$diet[Group == 'cod', whales := 0.1]
AB.params$diet[Group == 'cod', seals := 0.04]
AB.params$diet[Group == 'cod', whiting := 0.05]

AB.params$diet[Group == 'whiting', whales := 0.1]
AB.params$diet[Group == 'whiting', seals := 0.05]
AB.params$diet[Group == 'whiting', cod := 0.05]
AB.params$diet[Group == 'whiting', whiting := 0.05]

AB.params$diet[Group == 'mackerel', whales := 0.2]
AB.params$diet[Group == 'mackerel', mackerel := 0.05]

AB.params$diet[Group == 'anchovy', whales := 0.5]
AB.params$diet[Group == 'anchovy', cod := 0.10]
AB.params$diet[Group == 'anchovy', whiting := 0.45]
AB.params$diet[Group == 'anchovy', mackerel := 0.50]

AB.params$diet[Group == 'shrimp', seals := 0.01]
AB.params$diet[Group == 'shrimp', cod := 0.01]
AB.params$diet[Group == 'shrimp', whiting := 0.01]

AB.params$diet[Group == 'benthos', whales := 0.1]
AB.params$diet[Group == 'benthos', seals := 0.90]
AB.params$diet[Group == 'benthos', cod := 0.84]
AB.params$diet[Group == 'benthos', whiting := 0.44]
AB.params$diet[Group == 'benthos', shrimp := 1]
AB.params$diet[Group == 'benthos', benthos := 0.1]

AB.params$diet[Group == 'zooplankton', mackerel := .45]
AB.params$diet[Group == 'zooplankton', anchovy := 1]
AB.params$diet[Group == 'zooplankton', benthos := .1]

AB.params$diet[Group == 'phytoplankton', benthos := .1]
AB.params$diet[Group == 'phytoplankton', zooplankton := .9]

AB.params$diet[Group == 'detritus', benthos := .7]
AB.params$diet[Group == 'detritus', zooplankton := .1]
```

We now should enter the basic input parameters. Fortunately, the biologists have been busy, and we have some survey estimates of biomasses in the bay. The biomasses must be entered with the unit:t/km2. Whales: 10 individuals with an average weight of 800kg. Seals: 203 individuals with an average weight of 30kg. Cod: 310t. Whiting 170t. Mackerel: 122t. Anchovy: 600t. Shrimp: 0.8t/km2. Zooplankton: 14.8t/km2. Detritus: 10t/km2.

Next are production/biomass ratios, which with certain assumptions correspond to the total mortality, Z. The unit is year\-1, and we can often get Z from assessments. Alternatively, we have Z = F + M, so if we have the catch and the biomass, we can estimate F = C/B and add the total natural mortality to get Z. We can get estimates of M and Q/B from FishBase. Search for the species, (Gadus morhua, Merlangius merlangus, Scomber scombrus, Engraulis encrasicolus), find the life-history table, and extract the values. Estimate Z = F + M.

```{r,  include=FALSE}
AB.params$model[Type < 3, BioAcc  := 0]
AB.params$model[Type < 2, Unassim := 0.2]
AB.params$model[Type == 2, Unassim := 0]
#Detrital Fate
AB.params$model[Type < 2, detritus := 1]
AB.params$model[Type > 1, detritus := 0]
```

**Check for issues in your parameter file with**
\vspace{10pt}

```{r}
check.rpath.params(AB.params)
```

Once parameter file is built use this to run ecopath

STATUS BALANCED!

```{r, include=FALSE}
AB <- rpath(AB.params, 'Anchovy Bay')
AB
print(AB, morts = T)
webplot(AB)
```

Running Rsim: 3 step process 
Set the scene with rsim.scenario

```{r}
AB.scene <- rsim.scenario(AB, AB.params, 1:25)
```
\vspace{20pt}
**Food web plot**

```{r}
webplot(AB,
  eco.name = attr(AB, "eco.name"),  line.col = "grey",  labels = TRUE)
```

\newpage

The Production Biomass ratio (P/B) with certain assumptions correspond to the total mortality, and so the total mortality could be calculated for the remaining species to determine their PB. The total mortality (Z) was found using two methods: estimating Z using natural mortality (M) with the fishing mortality (F), and estimating Z with growth parameters using the Beverton-Holt equation.
The total mortality (Z) was calculated from the data (Z from data) using the equations: 

                                      Z=F+M and F=C/B,

where M is the natural mortality, F is the fishing mortality, C is the Catch in t/km^2, and B is the biomass in t/km^2. The natural mortality for seals was provided (M = 0.09 year\-1)). Estimates of natural mortality for cod (Gadus morhua), whiting (Merlangius merlangus), mackerel (Scomber scombrus), and anchovy (Engraulis encrasicolus) were taken from their respective life-history tables on FishBase (https://www.fishbase.se/search.php) Values of F for all species were determined using the estimates of catch and biomass of these species. Catch for each species was found by dividing the respective species fishery landings weight in tons by the area of Anchovy Bay (100 km^2). Similarly, the biomass for each species was found by dividing the respective species total survey estimate weight in tons by the area of Anchovy Bay (100 km^2). 

The total mortality (Z) was calculated from the Beverton-Holt equation (Z (B/H)) using: 

                          Z=K*((L_(inf )- L_mean ))/((L_(inf )- L')),
                          
where Lmean is the mean length of all fishes caught at sizes equal or larger than L’, which is the smallest size in the catch and here assumed to be the same as Lc, which is the mean length at entry in the fishery, assuming knife-edge selection, and thus the same as used under Yield per Recruit above.


+---------+---------+-------------+------+
| Species | Z (B&H) | Z from data | diff |
+:=======:+:=======:+:===========:+:====:+
| Cod     |  0.29   |    0.36     | -0.07|		
+---------+---------+-------------+------+		
| Whiting |  0.59   |    0.59     |  0.00|		
+---------+---------+-------------+------+		
| Mackerel|  0.63   |    1.15     | -0.51|		
+---------+---------+-------------+------+		
| Anchovy |  1.48   |    1.22     |  0.25|		
+---------+---------+-------------+------+		


The difference between the values of Z using these two different methods (diff) were examined for the fish species. Estimating Z from the data relies on estimates of catch, biomass, and natural mortality whereas estimating Z using Beverton-Holt relies on estimates of length. Estimates of Z from the data was found to be greater than those found from Beverton-Holt for cod and mackerel, while for anchovy it was lesser. For whiting, both methods resulted in the same value for Z.
\vspace{20pt}

**Exploring Ecosim with Anchovy Bay**
\vspace{10pt}
  
Scientists have noticed a gradual decline in the seal population of about 50 kg a year. The resource managers of Anchovy Bay decide that new regulations should be put in place on the seal fishery to combat this population decline. Local trawlers are concerned that an increased seal population will negatively impact their business as their chief targets are prey for seals.  

\newpage

**SCENARIOS**
\vspace{10pt}

**1) Build an rsim scenario that shows the impact of the seal decline and one that shows the impact of reducing the seal fishery in half.**
\vspace{10pt}

We examined the impact of what a reduced fishing effort on seals would have on the species targeted by the trawling fleet in Anchovy Bay. We found that a decrease in fishing effort on seals will produce higher seals biomass. We (as well as the trawling fleet in this scenario) anticipated that the increase in seal biomass would result in a decrease in the biomass of cod and whiting, their preferred prey. Instead, we found that there was an increase in cod biomass while there was a decrease in whiting biomass. Maybe due to the low diet composition value (0.04). As the trawling fleet was concerned about a decrease in the populations of their primary target species, and cod benefits from reduced fishing effort on seals, the trawling fleet will most likely not be dramatically impacted by a recovered seal population.

```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.base2 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'sealers',
                          sim.year = 1:25, value = 0.5)



AB.run2 <- rsim.run(AB.base2, years = 1:25)
#Output results
#AB.run2
#To capture results as an object use the write.Rsim function
AB.output2 <- write.Rsim(AB.run2)
```

```{r}
rsim.plot(AB.run2)
```

```{r}
seals <- extract.node(AB.run2, 'seals')
plot(seals$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Seals")

cod <- extract.node(AB.run2, 'cod')
plot(cod$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Cod")

whiting <- extract.node(AB.run2, 'whiting')
plot(whiting$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Whiting")
```

\newpage

**2) Also test the impact if seals are a strong ‘top-down’ predator. Run the simulations for 25 years.**
\vspace{10pt}

Then, seals were investigated to see what the impact would be if seals were a strong ‘top-down’ predator. To do this, we adjusted the first scenario to reflect a situation where seals had an increased effect on all of their prey.  
\vspace{10pt}

```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.base3 <- adjust.scenario(AB.base, parameter = 'VV', group = 'all', 
                           groupto = 'seals', value = 3)
```

```{r, include=FALSE}
AB.run3 <- rsim.run(AB.base3, years = 1:25)

#Output results
AB.run3
#To capture results as an object use the write.Rsim function
AB.output3 <- write.Rsim(AB.run3)
```

```{r}
#plot results
rsim.plot(AB.run3)
```

```{r}
seals <- extract.node(AB.run3, 'seals')
plot(seals$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Seals")

cod <- extract.node(AB.run3, 'cod')
plot(cod$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Cod")

whiting <- extract.node(AB.run3, 'whiting')
plot(whiting$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Whiting")
```
\vspace{10pt}

Several years later, ecotourism becomes a big part of Anchovy Bay. A portion of the trawlers shift their effort to taking customers whale watching. 

\newpage

**3) What is the impact on the groundfish species?** 
\vspace{10pt}

```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.base6 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 0.3)
```

```{r, include=FALSE}
AB.run6 <- rsim.run(AB.base6, years = 1:25)

#Output results
AB.run6
#To capture results as an object use the write.Rsim function
AB.output6 <- write.Rsim(AB.run6)
```

```{r}
#plot results
rsim.plot(AB.run6)
```

```{r}
seals <- extract.node(AB.run6, 'seals')
plot(seals$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Seals")

cod <- extract.node(AB.run6, 'cod')
plot(cod$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Cod")

whiting <- extract.node(AB.run6, 'whiting')
plot(whiting$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Whiting")
```

\newpage

**4) What if cod are a strong ‘top-down’ predator? It has been studied that whiting will increase their consumption when at lower numbers. What effect does this have on the previous scenario?**
\vspace{10pt}

```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.base5 <- adjust.scenario(AB.base, parameter = 'VV', group = 'all', 
                           groupto = 'cod', value = 3)
```

```{r, include=FALSE}
AB.run5 <- rsim.run(AB.base5, years = 1:25)

#Output results
AB.run5
#To capture results as an object use the write.Rsim function
AB.output5 <- write.Rsim(AB.run5)
```

Similar to the first scenario, cod were examined to see what the impact would be if cod were a strong ‘top-down’ predator. We adjusted the based model to reflect a situation where cod had an increased effect on all of their prey. The total biomass for cod increases, however, it increases overall to a lesser extent. For whiting, the plot of the biomass does not exhibit the initial peak as when cod is not a ‘top-down’ predator, and instead the biomass sharply declines from month 1 to about month 50, and then at a slower rate continuously declines until a very lower population equilibrium appears to be reached.

```{r}
cod <- extract.node(AB.run5, 'cod')
plot(cod$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Cod")

whiting <- extract.node(AB.run5, 'whiting')
plot(whiting$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Whiting")
```

This scenario was further built upon by including what the impact would be if the consumption by whiting increased as their numbers decreased. We found that the consumption of prey by whiting increases until a peak equilibrium appears to be reached. We also found that in response the biomass of anchovy, the main prey item for whiting, decreased sharply from month 1 to about month 50, and then a low biomass equilibrium appears to be reached. This makes sense as the biomass of whiting sharply declines during the same time period and then also appears to reach a low equilibrium.


```{r}
plot(AB.run5$annual_QB[,5], main = "Whiting consumption", xlab="years", ylab="Q/B ratio")

anchovy <- extract.node(AB.run5, 'anchovy')
plot(anchovy$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Anchovy")
```

\newpage

**5) Run scenarios varying level of trawling effort. Find the effort leading to maximum long-term yield for the trawl fishery. Plot and discuss implications for species biomass, mix of species in the trawl catch, and impact on other gears associated with fishing at maximum yield for trawl.**
\vspace{10pt}

```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.sc10 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 3.6)
AB.sc9 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 3.2)
AB.sc8 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 2.8)
AB.sc7 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 2.4)
AB.sc6 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 2.0)
AB.sc5 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 1.6)
AB.sc4 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 1.2)
AB.sc3 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 0.8)
AB.sc2 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 0.4)
AB.sc1 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'trawlers',
                          sim.year = 1:25, value = 0.1)
```

```{r, include=FALSE}
AB.run.sc10 <- rsim.run(AB.sc10, years = 1:25)
AB.run.sc9 <- rsim.run(AB.sc9, years = 1:25)
AB.run.sc8 <- rsim.run(AB.sc8, years = 1:25)
AB.run.sc7 <- rsim.run(AB.sc7, years = 1:25)
AB.run.sc6 <- rsim.run(AB.sc6, years = 1:25)
AB.run.sc5 <- rsim.run(AB.sc5, years = 1:25)
AB.run.sc4 <- rsim.run(AB.sc4, years = 1:25)
AB.run.sc3 <- rsim.run(AB.sc3, years = 1:25)
AB.run.sc2 <- rsim.run(AB.sc2, years = 1:25)
AB.run.sc1 <- rsim.run(AB.sc1, years = 1:25)
```

```{r, include=FALSE}
cod10 <- extract.node(AB.run.sc10, 'cod')
cod9 <- extract.node(AB.run.sc9, 'cod')
cod8 <- extract.node(AB.run.sc8, 'cod')
cod7 <- extract.node(AB.run.sc7, 'cod')
cod6 <- extract.node(AB.run.sc6, 'cod')
cod5 <- extract.node(AB.run.sc5, 'cod')
cod4 <- extract.node(AB.run.sc4, 'cod')
cod3 <- extract.node(AB.run.sc3, 'cod')
cod2 <- extract.node(AB.run.sc2, 'cod')
cod1 <- extract.node(AB.run.sc1, 'cod')

whit10 <- extract.node(AB.run.sc10, 'whiting')
whit9 <- extract.node(AB.run.sc9, 'whiting')
whit8 <- extract.node(AB.run.sc8, 'whiting')
whit7 <- extract.node(AB.run.sc7, 'whiting')
whit6 <- extract.node(AB.run.sc6, 'whiting')
whit5 <- extract.node(AB.run.sc5, 'whiting')
whit4 <- extract.node(AB.run.sc4, 'whiting')
whit3 <- extract.node(AB.run.sc3, 'whiting')
whit2 <- extract.node(AB.run.sc2, 'whiting')
whit1 <- extract.node(AB.run.sc1, 'whiting')
```

```{r, include=FALSE}
catch10<-last(cod10$AnnualTotalCatch)+last(whit10$AnnualTotalCatch)
catch9<-last(cod9$AnnualTotalCatch)+last(whit9$AnnualTotalCatch)
catch8<-last(cod8$AnnualTotalCatch)+last(whit8$AnnualTotalCatch)
catch7<-last(cod7$AnnualTotalCatch)+last(whit7$AnnualTotalCatch)
catch6<-last(cod6$AnnualTotalCatch)+last(whit6$AnnualTotalCatch)
catch5<-last(cod5$AnnualTotalCatch)+last(whit5$AnnualTotalCatch)
catch4<-last(cod4$AnnualTotalCatch)+last(whit4$AnnualTotalCatch)
catch3<-last(cod3$AnnualTotalCatch)+last(whit3$AnnualTotalCatch)
catch2<-last(cod2$AnnualTotalCatch)+last(whit2$AnnualTotalCatch)
catch1<-last(cod1$AnnualTotalCatch)+last(whit1$AnnualTotalCatch)
```


Scenario number 6 was the one that leads the trawl fishery to maximum long-term yield considering an effort value=2.0, obtaining a joint catch of cod and whiting of 0.7967.

```{r}
sce<-seq(1,10,1)
yields<-c(catch1,catch2,catch3,catch4,catch5,catch6,catch7, catch8,catch9,catch10)
plot(sce, yields)
```

To discuss the implications for species biomass we extract the biomass of cod and whiting and although the effort value of 2.0 maximizes the long-term yield, it produces a decrease in the biomass of the species present in the trawl fishery, therefore this strategy should not be considered.


```{r}

plot(cod6$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Cod")
plot(whit6$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Whiting")

```


**Discuss implications for species biomass, mix of species in the trawl catch, and impact on other gears associated with fishing at maximum yield for trawl.**

Regarding the impact on other gears associated with fishing at maximum yield for trawl. We use our long-term fishing effort (2.0) in the seiners fishery to see the impact on mackerel and anchovy and as a result, we obtained using the same fishing effort that leads to a long-term yield in the trawl fishery does not work for the seiner fishery, because it reduces the biomass of mackerel and anchovy to minimum levels. 


```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.base8 <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = 'seiners',
                          sim.year = 1:25, value = 2.0)
```

```{r, include=FALSE}
AB.run8 <- rsim.run(AB.base8, years = 1:25)

#Output results
AB.run8
#To capture results as an object use the write.Rsim function
AB.output8 <- write.Rsim(AB.run8)
```

```{r}
mackerel <- extract.node(AB.run8, 'mackerel')
plot(mackerel$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Mackerel")

anchovy <- extract.node(AB.run8, 'anchovy')
plot(anchovy$Biomass, xlab = 'Month', ylab = 'Biomass', main = "Anchovy")
```

\newpage

Finally, we develop three strategies to test in an MSE framework within Anchovy Bay. The objective was to maximize fishery yield but with a reference point of 0.5 BMSY as a management objective and not allowing any species to go below 0.25 BMSY as a limit reference point.


\vspace{10pt}


```{r, include=FALSE}
AB.base <- rsim.scenario(AB, AB.params, 1:25)
AB.run <- rsim.run(AB.base, years = 1:25)

#Determine reference point for cod----
#Turn off fishing to get a proxy for b0
gear <- AB.params$model[Type == 3, Group]
for(i in 1:length(gear)){
  AB.base <- adjust.fishing(AB.base, parameter = 'ForcedEffort', group = gear[i], 
                            value = 0, sim.year = 0:100)
}

AB.b0 <- rsim.run(AB.base, method = 'RK4', 1:25)
```

```{r, include=FALSE}
#Extract cod data and find b0
cod <- extract.node(AB.b0, 'cod')
cod.b0  <- max(cod$Biomass)

mackerel <- extract.node(AB.b0, 'mackerel')
mackerel.b0  <- max(mackerel$Biomass)

whiting <- extract.node(AB.b0, 'whiting')
whiting.b0  <- max(whiting$Biomass)

shrimp <- extract.node(AB.b0, 'shrimp')
shrimp.b0  <- max(shrimp$Biomass)

anchovy <- extract.node(AB.b0, 'anchovy')
anchovy.b0  <- max(anchovy$Biomass)

#Set reference point of 1/2 b0
cod.ref <- .5 * cod.b0
whiting.ref <- .5 * whiting.b0
mackerel.ref <- .5 * mackerel.b0
shrimp.ref <- .5 * shrimp.b0
anchovy.ref <- .5 * anchovy.b0

#Set limit point of 0.25 b0
cod.lim <- .25 * cod.b0
whiting.lim <- .25 * whiting.b0
mackerel.lim <- .25 * mackerel.b0
shrimp.lim <- .25 * shrimp.b0
anchovy.lim <- .25 * anchovy.b0
```

Objective values (0.5 BMSY) for:
Cod: 3.54158
Mackerel: 0.72248
Whiting: 0.978352
Shrimp: 0.45852
Anchovy: 2.9958

Limit values (0.25 BMSY) for:
Cod: 1.77
Mackerel: 0.361
Whiting: 0.489
Shrimp: 0.229
Anchovy: 1.4979

**STRATEGY 1**

First, we set our bio rule with an effort of "current.effort* 1.3" when the biomass is above the objective reference point and an effort of "current.effort * .75" when the biomass is below our limit reference point.

```{r, include=FALSE}
#control rule function----
#This modifies the effort matrix of the Rsim scenario object
bio.rule <- function(Rsim.scenario, Rsim.run, group, gear, ref.point, year){
  group.num <- which(Rsim.scenario$params$spname == group)
  gear.num  <- which(Rsim.scenario$params$spname == gear) - Rsim.scenario$params$NUM_BIO
  current.effort <- Rsim.scenario$fishing$ForcedEffort[(year - 1)*12 + 1, gear.num]
  if(Rsim.run$end_state$Biomass[group.num] > ref.point){
    Rsim.scenario <- adjust.fishing(Rsim.scenario, 'ForcedEffort', group = gear, 
                                    sim.year = year + 1, value = current.effort * 1.3)
  }
  if(Rsim.run$end_state$Biomass[group.num] < ref.point){
    Rsim.scenario <- adjust.fishing(Rsim.scenario, 'ForcedEffort', group = gear, 
                                    sim.year = year + 1, value = current.effort * .75)
  }
  return(Rsim.scenario)
}
```

```{r, include=FALSE}
#Run multistep scenario----
AB.base <- rsim.scenario(AB, AB.params, 1:100)
AB.init <- rsim.run(AB.base, method = 'AB', 1:25)
AB.full <- copy(AB.init)
for(i in 25:99){
  AB.base <- bio.rule(AB.base, AB.full, 'cod', 'trawlers', cod.ref, i)
  AB.full <- rsim.step(AB.base, AB.full, method = 'AB', i + 1)
}
```

On this scenario the cod biomass is constantly entering and leaving the target zone, but it stays far from the border zone. Mackerel remains below the target zone for much of the time series and also quite far from the limit zone while whiting, shrimp and anchovy remain above the target limit for the entire projection.

```{r}
cod <- extract.node(AB.full, 'cod')
plot(cod$Biomass, main="Cod", xlab = 'Month', ylab = 'Biomass', ylim=c(1.5,4.2))
abline(h = cod.ref, col="green")
abline(h = cod.lim, col="red")

mackerel <- extract.node(AB.full, 'mackerel')
plot(mackerel$Biomass, main="Mackerel",xlab = 'Month', ylab = 'Biomass', ylim=c(0.2,1.2))
abline(h = mackerel.ref, col="green")
abline(h = mackerel.lim, col="red")

whiting <- extract.node(AB.full, 'whiting')
plot(whiting$Biomass, main="Whiting",xlab = 'Month', ylab = 'Biomass', ylim=c(0.45,2))
abline(h = whiting.ref, col="green")
abline(h = whiting.lim, col="red")

shrimp <- extract.node(AB.full, 'shrimp')
plot(shrimp$Biomass,main="Shrimp", xlab = 'Month', ylab = 'Biomass', ylim=c(0,1.2))
abline(h = shrimp.ref, col="green")
abline(h = shrimp.lim, col="red")

anchovy <- extract.node(AB.full, 'anchovy')
plot(anchovy$Biomass, main="Anchovy",xlab = 'Month', ylab = 'Biomass', ylim=c(1,6))
abline(h = anchovy.ref, col="green")
abline(h = anchovy.lim, col="red")
```

Cod biomass last year: 3.794386
Mackerel biomass last year: 0.693222
Whiting biomass last year: 1.500186
Shrimp biomass last year: 0.897145
Anchovy biomass last year: 4.237256

Total: 11.12


**STRATEGY 2**

Then, we set our bio rule with an effort of "current.effort* 0.8" to implement a precautionary  and an effort of "current.effort * .60" when the biomass is below our limit reference point.


```{r, include=FALSE}
#control rule function----
#This modifies the effort matrix of the Rsim scenario object
bio.rule <- function(Rsim.scenario, Rsim.run, group, gear, ref.point, year){
  group.num <- which(Rsim.scenario$params$spname == group)
  gear.num  <- which(Rsim.scenario$params$spname == gear) - Rsim.scenario$params$NUM_BIO
  current.effort <- Rsim.scenario$fishing$ForcedEffort[(year - 1)*12 + 1, gear.num]
  if(Rsim.run$end_state$Biomass[group.num] > ref.point){
    Rsim.scenario <- adjust.fishing(Rsim.scenario, 'ForcedEffort', group = gear, 
                                    sim.year = year + 1, value = current.effort * 0.9)
  }
  if(Rsim.run$end_state$Biomass[group.num] < ref.point){
    Rsim.scenario <- adjust.fishing(Rsim.scenario, 'ForcedEffort', group = gear, 
                                    sim.year = year + 1, value = current.effort * .60)
  }
  return(Rsim.scenario)
}
```

```{r, include=FALSE}
#Run multistep scenario----
AB.base <- rsim.scenario(AB, AB.params, 1:100)
AB.init <- rsim.run(AB.base, method = 'AB', 1:25)
AB.full <- copy(AB.init)
for(i in 25:99){
  AB.base <- bio.rule(AB.base, AB.full, 'cod', 'trawlers', cod.ref, i)
  AB.full <- rsim.step(AB.base, AB.full, method = 'AB', i + 1)
}
```


in this scenario, cod increases its biomass and remains above the target point from month 300. Mackerel maintains its biomass below the limit point and whiting, shrimp and anchovy remain above the target point.

```{r}
cod <- extract.node(AB.full, 'cod')
plot(cod$Biomass, main="Cod", xlab = 'Month', ylab = 'Biomass', ylim=c(1.5,10.2))
abline(h = cod.ref, col="green")
abline(h = cod.lim, col="red")

mackerel <- extract.node(AB.full, 'mackerel')
plot(mackerel$Biomass, main="Mackerel",xlab = 'Month', ylab = 'Biomass', ylim=c(0.2,1.2))
abline(h = mackerel.ref, col="green")
abline(h = mackerel.lim, col="red")

whiting <- extract.node(AB.full, 'whiting')
plot(whiting$Biomass, main="Whiting",xlab = 'Month', ylab = 'Biomass', ylim=c(0.45,2))
abline(h = whiting.ref, col="green")
abline(h = whiting.lim, col="red")

shrimp <- extract.node(AB.full, 'shrimp')
plot(shrimp$Biomass,main="Shrimp", xlab = 'Month', ylab = 'Biomass', ylim=c(0,1.2))
abline(h = shrimp.ref, col="green")
abline(h = shrimp.lim, col="red")

anchovy <- extract.node(AB.full, 'anchovy')
plot(anchovy$Biomass, main="Anchovy",xlab = 'Month', ylab = 'Biomass', ylim=c(1,6))
abline(h = anchovy.ref, col="green")
abline(h = anchovy.lim, col="red")
```

Cod biomass last year: 7.531027
Mackerel biomass last year: 0.6079712
Whiting biomass last year: 1.611528
Shrimp biomass last year: 0.8271803
Anchovy biomass last year: 3.925502

Total: 14.5
 
**STRATEGY 3**

Finally, we used a value of "current.effort* 1.0" when the biomass is above the objective and  "current.effort * .80" when the biomass is below our limit reference point. Where mackerel was the only species that remained below the target point but, at the same time, well above the limit reference point.


```{r, include=FALSE}
#control rule function----
#This modifies the effort matrix of the Rsim scenario object
bio.rule <- function(Rsim.scenario, Rsim.run, group, gear, ref.point, year){
  group.num <- which(Rsim.scenario$params$spname == group)
  gear.num  <- which(Rsim.scenario$params$spname == gear) - Rsim.scenario$params$NUM_BIO
  current.effort <- Rsim.scenario$fishing$ForcedEffort[(year - 1)*12 + 1, gear.num]
  if(Rsim.run$end_state$Biomass[group.num] > ref.point){
    Rsim.scenario <- adjust.fishing(Rsim.scenario, 'ForcedEffort', group = gear, 
                                    sim.year = year + 1, value = current.effort * 1.0)
  }
  if(Rsim.run$end_state$Biomass[group.num] < ref.point){
    Rsim.scenario <- adjust.fishing(Rsim.scenario, 'ForcedEffort', group = gear, 
                                    sim.year = year + 1, value = current.effort * .80)
  }
  return(Rsim.scenario)
}
```

```{r, include=FALSE}
#Run multistep scenario----
AB.base <- rsim.scenario(AB, AB.params, 1:100)
AB.init <- rsim.run(AB.base, method = 'AB', 1:25)
AB.full <- copy(AB.init)
for(i in 25:99){
  AB.base <- bio.rule(AB.base, AB.full, 'cod', 'trawlers', cod.ref, i)
  AB.full <- rsim.step(AB.base, AB.full, method = 'AB', i + 1)
}
```

```{r}
cod <- extract.node(AB.full, 'cod')
plot(cod$Biomass, main="Cod", xlab = 'Month', ylab = 'Biomass', ylim=c(1.5,6))
abline(h = cod.ref, col="green")
abline(h = cod.lim, col="red")

mackerel <- extract.node(AB.full, 'mackerel')
plot(mackerel$Biomass, main="Mackerel",xlab = 'Month', ylab = 'Biomass', ylim=c(0.2,1.2))
abline(h = mackerel.ref, col="green")
abline(h = mackerel.lim, col="red")

whiting <- extract.node(AB.full, 'whiting')
plot(whiting$Biomass, main="Whiting",xlab = 'Month', ylab = 'Biomass', ylim=c(0.45,2))
abline(h = whiting.ref, col="green")
abline(h = whiting.lim, col="red")

shrimp <- extract.node(AB.full, 'shrimp')
plot(shrimp$Biomass,main="Shrimp", xlab = 'Month', ylab = 'Biomass', ylim=c(0,1.2))
abline(h = shrimp.ref, col="green")
abline(h = shrimp.lim, col="red")

anchovy <- extract.node(AB.full, 'anchovy')
plot(anchovy$Biomass, main="Anchovy",xlab = 'Month', ylab = 'Biomass', ylim=c(1,6))
abline(h = anchovy.ref, col="green")
abline(h = anchovy.lim, col="red")
```

Cod biomass last year: 4.50947
Mackerel biomass last year: 0.6690252
Whiting biomass last year: 1.477197
Shrimp biomass last year: 0.8782663
Anchovy biomass last year: 4.2006

Total: 11.7

**6) Show trade offs between fisheries. Assuming that cod sells for $3/pound, mackerel $2/pound, whiting $1/pound, shrimp $4/pound, and anchovies $0.50/pound is one of your strategies more economically beneficial to the community? The AB_closed_loop.R script can help you set-up the closed loop simulation.**

if we assume biomass is in pounds= 

strategy 1: 
Cod biomass last year: 3.79*3= 11.37
Mackerel biomass last year: 0.69*2= 1.38
Whiting biomass last year: 1.50*1= 1.5
Shrimp biomass last year: 0.897*4= 3.588
Anchovy biomass last year: 4.24*0.5= 2.12

Total: 19.958

strategy 2: 
Cod biomass last year: 7.531*3= 22.59
Mackerel biomass last year: 0.608*2= 1.216
Whiting biomass last year: 1.612*1= 1.612
Shrimp biomass last year: 0.827*4= 3.308
Anchovy biomass last year: 3.93*0.5= 1.965

Total: 30.691

strategy 3:
Cod biomass last year: 4.51*3= 13.53
Mackerel biomass last year: 0.669*2= 1.338
Whiting biomass last year: 1.477*1= 1.477
Shrimp biomass last year: 0.878*4= 3.512
Anchovy biomass last year: 4.2*0.5= 2.1

Total: 21.957