TestWeightAtStage.Rmd

---
title: "Testing atlantisom: visualize weight at (st)age output. How should we interplolate this for true ages?"
author: "Sarah Gaichas and Christine Stawitz"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
bibliography: "packages.bib"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
# automatically create a bib database for R packages
knitr::write_bib(c(
  .packages(), 'knitr', 'rmarkdown', 'tidyr', 'dplyr', 'ggplot2',
  'data.table', 'here', 'ggforce', 'ggthemes'
), 'packages.bib')
```

## Introduction

Now we delve into the calculation of true age classes and weight-at-age from standard Atlantis output age class, which is generally in 10 stages of maxage/10 true age classes for all but the newest Atlantis models. Note that modern Atlantis models can output true age classes in two output files (from the wiki):

>**ANNAGEBIO.nc**: This output provides numbers in each annual age class (so mapped from Atlantis "age class" which can contain multiple years to true annual age classes). Set flag_age_output to 1 to get this output. Tracers provided are:
>
>   + Numbers at age per species
>
>**ANNAGECATCH.nc**: This output provides numbers at annual age class (so mapped from Atlantis "age class" which can contain multiple years to true annual age classes) in the catch and discards (summed over all fleets). Set flag_age_output to 1 to get this output. Tracers provided are:
>
>   + Numbers at age per species in the catch
>   + Numbers at age per species in the discards

Future `atlantisom` users should take advantage of this direct output, and some `atlantisom` functions will have to be rewritten to use it.

However, The CCA model does not produce this output, and other legacy models may not, so we still need the `atlantisom` functions `calc_stage2age` and the associated `calc_Z`.  This function will be tested. However, it produces numbers at age but not weight. Weight-at-true-age is needed as an input to stock assessment models using empirical weight at age. See discussion [here](https://sgaichas.github.io/poseidon-dev/RethinkSamplingFunctions.html).

Therefore, our first step is understanding what the 10 class weight at age (output as `mulength` from `calc_age2length`) looks like for a variety of species.

For all setup, etc, please see previous files. Full methods are explained [here](https://sgaichas.github.io/poseidon-dev/TrueBioTest.html) and  [here](https://sgaichas.github.io/poseidon-dev/TrueLengthCompTest.html).

This page has visualizations for the NOBA model example, CERES Global Sustainability. At the end of the file we also review saved outputs from the CCA model example, Atlantis Summit Common Scenario 1.  

```{r message=FALSE, warning=FALSE}
library(tidyr)
require(dplyr)
library(ggplot2)
library(data.table)
library(here)
library(ggforce)
library(ggthemes)
library(atlantisom)
```

```{r initialize}

initCCA <- FALSE
initNEUS <- FALSE
initNOBA <- TRUE

if(initCCA) source(here("config/CCConfig.R"))

if(initNEUS) source(here("config/NEUSConfig.R"))

if(initNOBA) source(here("config/NOBAConfig.R"))

```

```{r get_names, message=FALSE, warning=FALSE}
#Load functional groups
funct.groups <- load_fgs(dir=d.name,
                         file_fgs = functional.groups.file)
#Get just the names of active functional groups
funct.group.names <- funct.groups %>% 
  filter(IsTurnedOn == 1) %>%
  select(Name) %>%
  .$Name

```

```{r load_Rdata, message=FALSE, warning=FALSE}

if(initCCA) {
  truth.file <- "outputCCV3run_truth.RData"
  load(file.path(d.name, truth.file))
  truth <- result
} 

if(initNEUS) {
  truth.file <- "outputneusDynEffort_Test1_run_truth.RData" 
  load(file.path(d.name, truth.file))
  truth <- result
}

if(initNOBA){
  truth.file <- "outputnordic_runresults_01run_truth.RData" 
  load(file.path(d.name, truth.file))
  truth <- result
}

```

## Simulate a survey part 5: weight at st(age) census

Full methods are explained [here](https://sgaichas.github.io/poseidon-dev/StdSurvLengthCompTest.html).

We will apply examples here to just a few NOBA species to examine within-year variability:
    - Cod "North_atl_cod", likely a test assessment species
    - Herring "Norwegian_ssh", likely a test assessment species
    - Greenland halibut "Green_halibut", which grows to a large size.

To create a census, the user specifies the timing of the survey, which species are captured, the spatial coverage of the survey, the species-specific survey efficiency ("q"), and the selectivity at age for each species. The following settings should achieve a survey that samples all Atlantis model output timesteps, all fish and shark species, and all model polygons, with perfect efficiency and full selectivity for all ages.


```{r census-spec, message=FALSE, warning=FALSE}

# should return a perfectly scaled survey 
effic1 <- data.frame(species=funct.group.names,
                     efficiency=rep(1.0,length(funct.group.names)))

# should return all lengths fully sampled (Atlantis output is 10 age groups per spp)
selex1 <- data.frame(species=rep(funct.group.names, each=10),
                     agecl=rep(c(1:10),length(funct.group.names)),
                     selex=rep(1.0,length(funct.group.names)*10))

# should return all model areas
boxpars <- load_box(d.name, box.file)
boxall <- c(0:(boxpars$nbox - 1))

# generalized timesteps all models
runpar <- load_runprm(d.name, run.prm.file)
noutsteps <- runpar$tstop/runpar$outputstep
stepperyr <- if(runpar$outputstepunit=="days") 365/runpar$toutinc

timeall <- c(0:noutsteps)
  
# define set of species we expect surveys to sample (e.g. fish only? vertebrates?)
# for ecosystem indicator work test all species, e.g.
survspp <- funct.group.names 

# for length and age groups lets just do fish and sharks
# NOBA model has InvertType, changed to GroupType in file, but check Atlantis default
if(initNOBA) funct.groups <- rename(funct.groups, GroupType = InvertType)

survspp <- funct.groups$Name[funct.groups$IsTurnedOn==1 &
                           funct.groups$GroupType %in% c("FISH", "SHARK")]

```

Here we use `create_survey` on the numbers output of `run_truth` to create the survey census of age composition (for just one species in this case). The `sample_fish` applies the median for aggregation and does not apply multinomial sampling if `sample=FALSE` in the function call. 

Because we don't want to wait 24 hours for this, we will look at only the first 112 time steps.

```{r stdsurveyNbased-3spp, echo=TRUE}

spp.name <- funct.group.names[funct.group.names %in% c("North_atl_cod",
                                                       "Norwegian_ssh",
                                                       "Green_halibut")]

# get survey nums with full (no) selectivity
survey_testNall <- create_survey(dat = truth$nums,
                                 time = c(0:111),
                                 species = spp.name,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# this one is high but not equal to total for numerous groups
effNhigh <- data.frame(species=survspp, effN=rep(1e+8, length(survspp)))

# apply default sample fish as before to get numbers
numsallhigh <- sample_fish(survey_testNall, effNhigh)


# aggregate true resn per survey design
aggresnall <- aggregateDensityData(dat = truth$resn,
                                 time = c(0:111),
                                 species = spp.name,
                                 boxes = boxall)

# aggregate true structn per survey design
aggstructnall <- aggregateDensityData(dat = truth$structn,
                                 time = c(0:111),
                                 species = spp.name,
                                 boxes = boxall)

#dont sample these, just aggregate them using median
structnall <- sample_fish(aggstructnall, effNhigh, sample = FALSE)

resnall <-  sample_fish(aggresnall, effNhigh, sample = FALSE)

```

Length sample with user specified max length bin (200 cm):

```{r userset-maxlen, echo=TRUE}

length_census <- calc_age2length(structn = structnall,
                                 resn = resnall,
                                 nums = numsallhigh,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 200,
                                 CVlenage = 0.1, remove.zeroes=TRUE)


```

We should get the upper end of Greenland halibut with a 200cm max length bin. That shouldnt matter for weight at age. This is all timesteps for each species by agecl.

```{r vis-wtageclass-test}

wageplot <- ggplot(length_census$muweight, aes(agecl, atoutput)) +
  geom_point(aes(colour = time)) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  scale_x_discrete(limits=c(1:10)) +
  xlab("age class") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste(scenario.name))

wageplot + facet_wrap(c("species"), scales="free_y")

```

How much does weight at age class vary over time? Here are a couple of annual cycles:

```{r vis-wtageclass-test2}

wtageyr1 <- length_census$muweight %>%
  filter(time<5)

wageplot <- ggplot(wtageyr1, aes(agecl, atoutput)) +
  geom_line(aes(colour = factor(time))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  scale_x_discrete(limits=c(1:10)) +
  xlab("age class") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " initial year"))

wageplot + facet_wrap(c("species"), scales="free_y")

```

```{r vis-wtageclass-test3}

wtageyrmid <- length_census$muweight %>%
  filter(between(time, 50,54))

wageplot <- ggplot(wtageyrmid, aes(agecl, atoutput)) +
  geom_line(aes(colour = factor(time))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  scale_x_discrete(limits=c(1:10)) +
  xlab("age class") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " year 10"))

wageplot + facet_wrap(c("species"), scales="free_y")

```

```{r vis-wtageclass-test4}

wtageyr20 <- length_census$muweight %>%
  filter(between(time, 100,104))

wageplot <- ggplot(wtageyr20, aes(agecl, atoutput)) +
  geom_line(aes(colour = factor(time))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  scale_x_discrete(limits=c(1:10)) +
  xlab("age class") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " year 20"))

wageplot + facet_wrap(c("species"), scales="free_y")

```

There are some interesting jumps between age classes. Most pronounced for herring. All of these species would be split into two true ages per age class. Interpolation could be interesting.

Change in wt at age over time for age classes using an annual mid-year snapshot (survey) (first 22+ years of NOBA model run):

```{r aggwtcomp}

# from std survey code and model step info read in above
midptyr <- round(median(seq(1,stepperyr)))
annualmidyear <- seq(midptyr, noutsteps, stepperyr)

wtage_annsurv <- length_census$muweight %>%
  filter(time %in% annualmidyear)

# reverse to show agecl time series of wt
wageplot <- ggplot(wtage_annsurv, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (5 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")

```


## Quick test with saved CCA census output

Here we check the CCA model (Atlantis Summit Common Scenario 1) saved length and weight census output for small species at Isaac's suggestion (all three of the following are true age classes, NumAgeClassSize=1 per Atlantis agecl bin:

```{r loadCCAlengthcomp, echo=TRUE}

source(here("config/CCConfig.R"))

length_censussurvsamp <- readRDS(file.path(d.name, paste0(scenario.name, "length_censussurvsamp.rds")))

```

Sardine:  
```{r plotCCAsardine}
censuswt <- length_censussurvsamp$muweight %>%
  filter(species == "Pacific_sardine")

wageplot <- ggplot(censuswt, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (1 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")
```

Anchovy:  
```{r plotCCAanchovy}
censuswt <- length_censussurvsamp$muweight %>%
  filter(species == "Anchovy")

wageplot <- ggplot(censuswt, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (1 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")
```

Herring:  
```{r plotCCAherring}
censuswt <- length_censussurvsamp$muweight %>%
  filter(species == "Herring")

wageplot <- ggplot(censuswt, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (1 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")
```

How do these three look? Herring had a rough period from model years 15-30 or so.

Now for some CCA groups with multiple true ages per stage.

Pacific hake, 2 true ages per class:  
```{r plotCCAhake}
censuswt <- length_censussurvsamp$muweight %>%
  filter(species == "Mesopel_M_Fish")

wageplot <- ggplot(censuswt, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (1 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")
```

Bocaccio rockfish, 5 true ages per age class:
```{r plotCCAbocaccio}
censuswt <- length_censussurvsamp$muweight %>%
  filter(species == "Bocaccio_rockfish")

wageplot <- ggplot(censuswt, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (1 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")
```

Yelloweye rockfish, most extreme at 12 true ages per age class:
```{r plotCCAYelloweye}
censuswt <- length_censussurvsamp$muweight %>%
  filter(species == "Yelloweye_rockfish")

wageplot <- ggplot(censuswt, aes(time, atoutput)) +
  geom_line(aes(colour = factor(agecl))) +
  theme_tufte() +
  theme(legend.position = "bottom") +
  xlab("model timestep (1 per year)") +
  ylab("average individual weight (g)") +
  labs(subtitle = paste0(scenario.name, " annual mid year survey"))

wageplot + facet_wrap(c("species"), scales="free_y")
```

I'm not sure any of this helps us decide how to interpolate, but there is clearly contrast in weight at age (class) in both models and example runs.