---
title: "Generating an Rsim ensemble with Ecosense"
author: "George A. Whitehouse"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Generating an Rsim ensemble with Ecosense}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(Rpath); library(data.table); library(viridis)
#knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
```

Rpath is an R implementation of the food web modeling program Ecopath with Ecosim
(EwE; Christensen and Pauly 1992^[Christensen V, Pauly D (1992) Ecopath II-a software 
for balancing steady-state ecosystem models and calculating network characteristics. 
Ecol Model 61:169-185. doi:10.1016/0304-3800(92)90016-8], Walters et al. 1997^[Walters 
C, Christensen V, Pauly D (1997) Structuring dynamic models of exploited ecosystems 
from trophic mass-balance assessments. Rev Fish Biol Fish 7:139-172. 
doi:10.1023/a:1018479526149]). For full documentation of the Rpath package see Lucey 
et al. (2020^[Lucey SM, Gaichas SK, Aydin KY (2020) Conducting reproducible ecosystem 
modeling using the open source mass balance model Rpath. Ecol Model 427:11. 
doi:10.1016/j.ecolmodel.2020.109057]). Ecosense is a Monte Carlo approach to generating 
an ensemble of plausible ecosystem parameter sets from a single Rpath model for 
use with `rsim.run`, the time-dynamic counterpart to Rpath. For complete documentation
of Ecosense see Whitehouse and Aydin (2020^[Whitehouse GA, Aydin KY (2020) 
Assessing the sensitivity of three Alaska marine food webs to perturbations: an 
example of Ecosim simulations using Rpath. Ecol Model 429:16. 
doi:10.1016/j.ecolmodel.2020.109074]). To use this vignette it is recommended you
are familiar with Rpath and its basic operations.

## Ecosense

This vignette includes the example Rpath models from Whitehouse and 
Aydin (2020) for the eastern Bering Sea (EBS), eastern Chukchi Sea (ECS), and Gulf 
of Alaska (GOA) marine food webs, and uses the EBS model in an example. Each of 
these models has 50–54 functional groups, including one fishery, one primary 
producer, two detrital compartments, and there are no stanzas. To use Ecosense you 
need an `rpath.params` object and the corresponding `rsim.scenario` object. First, 
load the unbalanced models, balance them, and setup the rsim scenario objects.

```{r loadbalanacemodels}
# load the unbalanced models
# load("data/Ecosense.EBS.rda")
# load("data/Ecosense.ECS.rda")
# load("data/Ecosense.GOA.rda")
# balance the models
EBS_bal <- rpath(Ecosense.EBS)
ECS_bal <- rpath(Ecosense.ECS)
GOA_bal <- rpath(Ecosense.GOA)
# create rsim scenario objects
EBS_scene <- rsim.scenario(EBS_bal, Ecosense.EBS, years=1:100)
ECS_scene <- rsim.scenario(ECS_bal, Ecosense.ECS, years=1:100)
GOA_scene <- rsim.scenario(GOA_bal, Ecosense.GOA, years=1:100)
# # source ecosense.R
# source("R/ecosense.R", local = knitr::knit_global())
# ls()
```
### Setting up Ecosense

Full sets of Ecosim parameter sets are drawn from distributions centered on the
original Rpath parameter estimates. The width of the distribution for each parameter
is defined in the data pedigree of the `rpath.params` object. The drawn parameters
include biomass, P/B, Q/B, diet composition, and natural mortality (i.e., M zero).
Additionally, vulnerability and handling time from the predator prey functional 
response can be drawn over a specified range (lower/upper bounds in log space-1) 
for each individual predator-prey interaction.

A single parameter set can be generated with a call to `rsim.sense`. The object 
returned by `rsim.sense` is equivalent to the list of parameters in an `rsim.scenario` 
object (e.g., scenario_object$params).

```{r genmod, echo=TRUE}
# One set of Ecosim parameters for the EBS model
rsim.sense(EBS_scene,Ecosense.EBS,Vvary = c(0,0), Dvary = c(0,0))
```

Each Ecosim parameter set drawn with Ecosense can be subject to an initial simulation, 
also known as a burn-in, to eliminate unstable parameter sets. This instability 
usually results from incompatible draws of parameter sets that either lead to uncontrolled
population growth or population crash. Such as predator consumption and production 
that exceeds the production of prey. The length of the burn-in period can be set
in the original `rsim.scenario` object.

```{r init, echo=TRUE}
# Setting the burn-in period in the EBS scenario object to 50 years.
EBS_scene$params$BURN_YEARS <- 50
```
During burn-in, if the biomass of a functional group exceeds 1,000x its starting
biomass or decreases to less than 1/1,000 its starting biomass, the ecosystem parameter
set is rejected and not retained for further analysis. A 50 year burn-in period 
is generally sufficient to eliminate most of the unstable configurations.

### Generating an ensemble of Ecosim parameter sets

To generate an ensemble of parameter sets you can repeat the call to `rsim.sense` 
in a loop. First, determine the number of parameter sets you wish to generate and 
create lists to store those parameter sets. Set up another vector to keep track 
of which generated parameter sets were rejected or retained.
```{r ensemble, echo=TRUE}
NUM_RUNS <- 1000 # how many ecosystem parameter sets to generate
parlist<-as.list(rep(NA,NUM_RUNS)) # create lists to store the generated parameters
kept<-rep(NA,NUM_RUNS) # object to keep track of kept systems
set.seed(666) # Optional, set seed so output can be replicated
```
In this example we will generate 1,000 parameter sets (NUM_RUNS) for the eastern 
Bering Sea food web model. parlist[[i]] is the object that will store the ith 
generated parameter set. Optionally, you can use `set.seed` to replicate these results.

The following loop generates an ensemble of Ecosim parameters for the EBS food web 
model, assigns the generated base biomass (B_BaseRef) to the starting biomass in 
the `rsim.scenario` object, sets the burn years to 50, runs a simulation with `rsim.run`, 
and rejects or retains parameter sets based on the boundaries described above.

```{r generatorloop, echo=TRUE, results='hide'}
for (i in 1:NUM_RUNS){
  EBSsense <- EBS_scene # scenario object
  # INSERT SENSE ROUTINE BELOW
  parlist[[i]]<- EBS_scene$params 		# Base ecosim params
  parlist[[i]]<- rsim.sense(EBS_scene,Ecosense.EBS,Vvary = c(-4.5,4.5), Dvary = c(0,0))	# Replace the base params with Ecosense params
  EBSsense$start_state$Biomass <- parlist[[i]]$B_BaseRef # Apply the Ecosense starting biomass
  parlist[[i]]$BURN_YEARS <- 50			# Set Burn Years to 50
  EBSsense$params <- parlist[[i]] # replace base params with the Ecosense generated params
  EBStest <- rsim.run(EBSsense, method="AB") # Run rsim with the generated system
  failList <- which(is.na(EBStest$end_state$Biomass))
  {if (length(failList)>0)
  {cat(i,": fail in year ",EBStest$crash_year,": ",EBStest$params$spname[failList],"\n"); kept[i]<-F; flush.console()}
    else 
    {cat(i,": success!\n"); kept[i]<-T;  flush.console()}} # output for the console
  parlist[[i]]$BURN_YEARS <- 1
}
```
This loop produces output in the console noting whether the generated parameter 
set was a success or failure, the year the ecosystem "crashed", and which group(s)
were responsible for the crash. Here are the results for the first two generated
parameter sets:

  1 : fail in year  1 :  Squids\
  2 : fail in year  7 :  Salmon returning\

The first two ecosystems both failed, the first in year one and the second in year
seven. Squids crashed in the first ecosystem and Salmon returning in the second. 
The 52nd generated parameter set was the first to survive the 50 year burn-in:

  52 : success!\

To determine which of the 1,000 generated systems survived burn-in, how many survived burn-in and what the rejection rate was:
```{r status, echo=TRUE}
KEPT <- which(kept==T); KEPT # the number associated with the kept system
nkept <- length(KEPT); nkept # how many were kept
1-(nkept/NUM_RUNS) # rejection rate
```
The numbers in KEPT can be used to access the retained parameter sets in parlist[[i]].
In this example, 43 parameter sets were retained and the rejection rate for generated
parameter sets was 95.7%.

All of the retained Rsim parameter sets can be run through a simulation with another
for loop. This loop subjects each of the 43 retained parameter sets in this example 
to a 100 year simulation without any additional perturbations, the first 50 years 
of which is the burn-in period.
```{r simstuff, echo=TRUE, results='hide'}
ecos <- as.list(rep(NA,length(KEPT))) # lists for simulated ecosystems
k <- 0  # counter for simulated ecosystems
for (i in KEPT) {
  EBS_scene$start_state$Biomass <- parlist[[i]]$B_BaseRef # set the starting Biomass to the generated values
  EBSsense <- EBS_scene # set up the scenario object
  EBSsense$params <- parlist[[i]] # set the params in the scenario object equal to the generated params.
  EBSsense$BURN_YEARS <- -1 # no burn-in period
  k <- k + 1 # set the number for the simulated ecosystem
  ecos[[k]] <- rsim.run(EBSsense,method='AB') # run rsim.run on the generated system
  print(c("Ecosystem no.",k,"out of",nkept)) # progress output to console
}
```

Because biomass can differ greatly between the generated parameter sets, relative
biomass can be used for plotting. This loop calculates the relative biomass of each 
group in each system, relative to their starting biomass, as drawn by Ecosense.
```{r relbio, echo=TRUE}
relB_ecos <- as.list(rep(NA,length(KEPT))) # list to output relative biomass
k <- 0
for (i in 1:nkept) {
  spname <- colnames(ecos[[i]]$out_Biomass[,2:ncol(ecos[[i]]$out_Biomass)])
  biomass <- ecos[[i]]$out_Biomass[, spname]
  n <- ncol(biomass)
  start.bio <- biomass[1, ] # the drawn starting biomass 
  start.bio[which(start.bio == 0)] <- 1
  rel.bio <- matrix(NA, dim(biomass)[1], dim(biomass)[2])
  for(isp in 1:n) rel.bio[, isp] <- biomass[, isp] / start.bio[isp] # biomass relative to biomass at t=1
  colnames(rel.bio) <- spname
  k <- k + 1
  relB_ecos[[k]] <- rel.bio
}
```
Here's an example plot of walleye pollock from all the retained Rsim parameter sets. 
First setup a matrix to store the pollock trajectories from the simulations with 
the retained parameter sets.
```{r walleye, echo=TRUE}
this_species <- "Walleye pollock"
plot_mat <- matrix(nrow=1200, ncol=nkept) # matrix of pollock trajectories from all generated systems
for (i in 1:nkept) {
  plot_mat[,i] <- relB_ecos[[i]][,this_species]
}
```
Then plot all the trajectories:
```{r plottraj, echo=TRUE}
plot_col <- viridis(nkept)
#layout(matrix(c(1,1,1,1,1,1,2,2), nrow = 1, ncol = 8, byrow = TRUE))
plot(1:1200, relB_ecos[[1]][,this_species], type='n', xlab="Months",
     ylab="Relative biomass", ylim=c(min(plot_mat),max(plot_mat)), main=this_species)
# one line for pollock in each of the generated systems
for (i in 1:nkept) {
  lines(1:1200, relB_ecos[[i]][,this_species], lwd=2, col=plot_col[i])
}
# distribution of pollock relative biomasses
boxplot(plot_mat[1200,], ylim=c(min(plot_mat),max(plot_mat)), yaxt='n')
axis(side=2, at=c(seq(0,150,50)), tick=TRUE, labels=F)
```
Each line in the left panel is the relative biomass trajectory of walleye pollock
over the 100 year simulation with the retained Rsim parameter sets.

### Run the kept systems with a perturbation

When working with an ensemble such as this, we are often interested in doing the 
same perturbation across all ensemble members to describe a range of potential outcomes.
In this example, we'll increase the fishing mortality on walleye pollock using the 
`adjust.fishing` function. The first 50 years of the simulation are the burn-in 
period. So, the perturbation is put in place from year 51 to 100, and in this example
we increase walleye pollock fishing mortality by 2X.
```{r run, echo=TRUE, results='hide'}
ecos_sp <- as.list(rep(NA,length(KEPT))) # lists for simulated ecosystems
k <- 0  # counter for simulated ecosystems
for (i in KEPT) {
  EBS_scene$start_state$Biomass <- parlist[[i]]$B_BaseRef # set the starting Biomass to the generated values
  EBSsense <- EBS_scene # set up the scenario object
  EBSsense <- adjust.fishing(EBSsense, "ForcedFRate", group=this_species, sim.year=51:100, value=2) # perturb pollock FRate
  EBSsense$params <- parlist[[i]] # set the params in the scenario object equal to the generated params.
  EBSsense$BURN_YEARS <- -1 # no burn-in period
  k <- k + 1 # set the number for the simulated ecosystem
  ecos_sp[[k]] <- rsim.run(EBSsense,method='AB') # run rsim.run on the generated system
  print(c("Ecosystem no.", k, "out of", nkept)) # progress output to console
}
```
To visualize the results, we use relative biomass again but this time relative to 
biomass at the end of the burn-in period (year 50). That is the point in time where
the perturbation began and the point from which we want to measure any response 
or change.
```{r viz, echo=TRUE}
relB_ecos_sp <- as.list(rep(NA,length(KEPT)))
k <- 0
for (i in 1:nkept) {
  spname <- colnames(ecos_sp[[i]]$out_Biomass[,2:ncol(ecos_sp[[i]]$out_Biomass)])
  biomass <- ecos_sp[[i]]$out_Biomass[, spname]
  n <- ncol(biomass)
  start.bio <- biomass[600, ] # end of burn-in
  start.bio[which(start.bio == 0)] <- 1
  rel.bio <- matrix(NA, dim(biomass)[1], dim(biomass)[2])
  for(isp in 1:n) rel.bio[, isp] <- biomass[, isp] / start.bio[isp] # biomass relative to the end of burn-in biomass
  colnames(rel.bio) <- spname
  k <- k + 1
  relB_ecos_sp[[k]] <- rel.bio
}
```
Plot the trajectories for walleye pollock from this perturbation. Note, the x-axis 
starts at the begining of the perturbation (January of year 51).
```{r walleyeplot, echo=TRUE}
plot_mat_sp <- matrix(nrow=1200, ncol=nkept)
for (i in 1:nkept) {
  plot_mat_sp[,i] <- relB_ecos_sp[[i]][,this_species]
}
#layout(matrix(c(1,1,1,1,1,1,2,2), nrow = 1, ncol = 8, byrow = TRUE))
plot(601:1200, relB_ecos_sp[[1]][601:1200,this_species], type='n', xlab="Months",
     ylab="Relative biomass", ylim=c(min(plot_mat_sp[601:1200,]),max(plot_mat_sp[601:1200,])), main=this_species)
for (i in 1:nkept) {
  lines(601:1200, relB_ecos_sp[[i]][601:1200,this_species], lwd=2, col=plot_col[i])
}
boxplot(plot_mat_sp[1200,], ylim=c(min(plot_mat_sp[601:1200,]),max(plot_mat_sp[601:1200,])))
```
Biomass for walleye pollock decreased in all 43 ecosystems in response to the perturbation.

Walleye pollock are a nodal species in the eastern Bering Sea food web and such 
a perturbation as in this example might illicit a response from many functional 
groups across the food web. Here is an example boxplot of the distribution 
of biomass outcomes for all the living functional groups relative to their biomass 
at the start of the perturbation.
```{r perturb, echo=TRUE}
sp_perturb_out <- matrix(nrow=nkept, ncol=(EBS_bal$NUM_LIVING+EBS_bal$NUM_DEAD))
for(i in 1:nkept){
  sp_perturb_out[i,] <- relB_ecos_sp[[i]][1200,]
}
colnames(sp_perturb_out) <- colnames(relB_ecos_sp[[1]])
par(mfrow=c(1,1), mar=c(9,3,0.5,0.5))
boxplot(sp_perturb_out, outline=FALSE, las=2, cex.axis=0.6)
abline(h=1, lty=2)
```