Atlantisom2SSworkflowtest.Rmd

---
title: "Testing atlantisom: end to end truth for California Current"
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 = TRUE)
# 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

This page documents initial testing of the atlantisom package in development at https://github.com/r4atlantis/atlantisom using three different [Atlantis](https://research.csiro.au/atlantis/) output datasets. Development of atlantisom began at the [2015 Atlantis Summit](https://research.csiro.au/atlantis/atlantis-summit/) in Honolulu, Hawaii, USA. 

The purpose of atlantisom is to use existing Atlantis model output to generate input datasets for a variety of models, so that the performance of these models can be evaluated against known (simulated) ecosystem dynamics. Atlantis models can be run using different climate forcing, fishing, and other scenarios. Users of atlantisom will be able to specify fishery independent and fishery dependent sampling in space and time, as well as species-specific catchability, selectivty, and other observation processes for any Atlantis scenario. Internally consistent multispecies and ecosystem datasets with known observation error characteristics will be the atlantisom outputs, for use in individual model performance testing, comparing performance of alternative models, and performance testing of model ensembles against "true" Atlantis outputs.

On this page we demonstrate use of atlantisom on the California Current output files to provide inputs for a stock assessment model implemented in Stock Synthesis [insert SSREF]. 

## Setup

First, you will want to set up libraries and install atlantisom if you haven't already. This document is written in in R Markdown [@R-rmarkdown], and we use several packages to produce the outputs [@R-tidyr; @R-dplyr; @R-ggplot2; @R-here; @R-ggforce; @R-ggthemes]. 

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

If you want to load the local of the atlantisom R package, use devtools. The R package here is helpful for installing from your local atlantisom directory, as it avoids hardcoding the specific location. The code below is not running from our local atlantisom directory and is not evaluated:

```{r eval=FALSE, message=FALSE, warning=FALSE}
require(devtools)
package.dir <- here()
devtools::load_all(package.dir)
```

Or you can install directly from the Github repository.

```{r eval=FALSE, message=FALSE, warning=FALSE}
devtools::install_github("r4atlantis\atlantisom")
```

Next, load the package. 

```{r message=FALSE, warning=FALSE}
library(atlantisom)
```

## Initializing input files and directories

You will first need to tell `atlantisom` to know where to look for the output and input files from your atlantis model run. Here, we will give the directory where the Atlantis inputs and outputs are stored `d.name`, the location of the functional groups file `functional.group.file` (.csv), the biomass pools file `biomass.pools.file` (.nc), the box locations file `box.file` (.bgm), an `initial.conditions.file` (.nc), the biology .prm file `biol.prm.file` (.prm), and the run .prm file `run.prm.file` (.prm). You will also need to specify a scenario name, which will be used to define the output files (i.e. output is stored in a number of netCDF files of the format: output<scenario><value>.nc). All of these files should be stored in `d.name`.

In general, Atlantis model output files that are sufficiently detailed to mimic fishery sampling in space and time are too large to store on GitHub, so are not included as examples with the atlantisom code. 

In this example, we have a local folder "atlantisoutput" and the output of interest in a subdirectory:

  * "CalCurrent2013_OA_off": California Current Atlantis model output from Isaac Kaplan's google drive \outputFolderAllOff2013Oceanography using the Hodgson et al. 2018 configuration. Corresponds to CC1Config.R file.

Our directory structure is set up to take advantage of `here()` to allow setup on a different computer.

```{r initialize}

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

#function to make a config file? need one for each atlantis run

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

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

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

```

## Getting the "true" operating model values

There are a number of functions in the package that begin with the prefix `load` that load various files. See documentation if you'd only like to load one file. The `atlantisom::run_truth()` function uses the above file definitions and calls a number of the `load` functions to read in all of the atlantis output. Note: this call reads in a number of large .nc files, so it will take a few minutes to return.

```{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 get_truth, message=FALSE, warning=FALSE}

# default run_truth setup will save the file, so check for that first

if(!file.exists(file.path(d.name, 
                          paste0("output", scenario.name, "run_truth.RData")))){
  #Store all loaded results into an R object
  truth <- run_truth(scenario = scenario.name,
                     dir = d.name,
                     file_fgs = functional.groups.file,
                     file_bgm = box.file,
                     select_groups = funct.group.names,
                     file_init = initial.conditions.file,
                     file_biolprm = biol.prm.file,
                     file_runprm = run.prm.file
  )
} else{
  truth <- get(load(file.path(d.name,
                              paste0("output", scenario.name, "run_truth.RData"))))
}
```

Now the R object `truth` with the comprehensive results from the Atlantis model has been read in. It is also saved as an .RData file titled "output[scenario.name]run_truth.RData" in the directory with the model output, so that later analyses can use `base::load()` instead of taking the time to rerun `atlantisom::run_truth()`. 

Note: on Sarah's laptop, CCA took ~2.5 hours to return from `run_truth` and NOBA took an hour and 26 minutes to return. Therefore, reading the .RData file in as below is advised after the first run to make the plots comparing systems later in this document.

## Test survey functions: census to compare with truth

This section tests the `atlantisom::create_survey()` and `atlantisom::sample_survey_biomass()` functions by comparing a census (survey sampling everything) with the results generated by `atlantisom::run_truth()` above. 

To create a survey, 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 species, and all model polygons, with perfect efficiency and full selectivity for all ages: 

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

# make a function for this
source(here("config/census_spec.R"))

```

### True biomass survey I

Because the results of `run_truth` provide both biomass at age and numbers at age, we can use `create_survey` on both with some modifications. The following uses the biomass output of `run_truth` to create the survey, so the call to `sample_survey_biomass` will require a weight at age argument that is filled with 1's because no conversion from numbers to weight is necessary. Because we are making a census for testing, the survey cv argument is set to 0. 

```{r surveyBbased, eval=F}

# this uses result$biomass_ages to sample biomass directly, no need for wt@age est

survey_testBall <- create_survey(dat = truth$biomass_ages,
                                 time = timeall,
                                 species = survspp,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# make up a constant 0 cv for testing
surv_cv <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

# call sample_survey_biomass with a bunch of 1s for weight at age
# in the code it multiplies atoutput by wtatage so this allows us to use
# biomass directly
wtage <- data.frame(species=rep(survspp, each=10),
                    agecl=rep(c(1:10),length(survspp)),
                    wtAtAge=rep(1.0,length(survspp)*10))

surveyB_frombio <- sample_survey_biomass(survey_testBall, surv_cv, wtage)

#save for later use, takes a long time to generate
saveRDS(surveyB_frombio, file.path(d.name, paste0(scenario.name, "surveyBcensus.rds")))

```

Comparing our (census) survey based on true biomass from above with the Atlantis output file "[modelscenario]BiomIndx.txt" should give us a perfect match. Note that the our (census) survey may have more sampling in time than the Atlantis output file.

```{r matchB, fig.cap="Testing whether the survey census gives the same results as the Atlantis output biomass index file; first 36 species.", message=FALSE, warning=FALSE}
# plot some comparisons with Atlantis output

# read Atlantis output files
atBtxt2 <- read.table(file.path(d.name, paste0("output", scenario.name, "BiomIndx.txt")), header=T)

surveyB_frombio <- readRDS(file.path(d.name, paste0(scenario.name, "surveyBcensus.rds")))
  
# lookup the matching names, put in time, species, biomass column format
# WARNING hardcoded for output with last species group as DIN
groupslookup <- funct.groups %>%
  filter(IsTurnedOn > 0)

atBtxt2tidy <- atBtxt2 %>%
  select(Time:DIN) %>%
  #select(Time, FPL:DIN) %>%
  rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, biomass, -Time) %>%
  filter(species %in% levels(surveyB_frombio$species))

#all species comparison, time intervals hardcoded for 5 steps per year
compareB <-ggplot() +
  geom_line(data=surveyB_frombio, aes(x=time/5,y=atoutput, color="survey census B"), 
            alpha = 10/10) +
  geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
             alpha = 1/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

compareB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 1, scales="free") 
compareB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 2, scales="free") 
compareB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 3, scales="free") 
compareB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 4, scales="free") 

```

After comparing the survey census based on the biomass output of `run_truth`, we now look at a biomass index that is estimated based on numbers in the survey census for assessments that take numbers input directly. We need to skip the average weight part of `sample_survey_biomass` but keep the rest. Hence the new function `sample_survey_numbers`.

### True numbers survey

```{r surveyNbased, eval=FALSE}

# this uses result$nums and a new function to get survey index in numbers (abundance)

survey_testNall <- create_survey(dat = truth$nums,
                                 time = timeall,
                                 species = survspp,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# as above, make up a constant 0 cv for testing
surv_cv <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

surveyN <- sample_survey_numbers(survey_testNall, surv_cv)

#save for later use, takes a long time to generate
saveRDS(surveyN, file.path(d.name, paste0(scenario.name, "surveyNcensus.rds")))

```

### Biological sampling: true numbers at age class

We get true numbers at age at each output step by appling `sample_fish` to the survey_testNall result if the model group has one true age per output age class:

```{r truenatagecl, eval=FALSE}
# We need an effective N for the sample fish function, setting equal to actual N
# this one is high but not equal to total for numerous groups
effNhigh <- data.frame(species=survspp, effN=rep(1e+8, length(survspp)))

numsallhigh <- sample_fish(survey_testNall, effNhigh)

#save for later use, takes a long time to generate
saveRDS(numsallhigh, file.path(d.name, paste0(scenario.name, "Natageclcensus.rds")))

```

### Biological sampling: setup

The next steps of sampling can get us a proper weight at age. Here again a census to generate length comps, which need an average weight at age and carry it forward. First we `sample_fish` (nums and weight components).

On a full model run this takes far too long with all species, so we cut it down to a few assessment species here:

```{r ss-biolsampling, warning=FALSE, message=FALSE}

# get only assessed species, hardcoded here for CCA sardine, hake, bocaccio
spp.name <- funct.group.names[funct.group.names %in% c("Pacific_sardine",
                                                       "Mesopel_M_Fish",
                                                       "Bocaccio_rockfish")]

survey_ssN <- create_survey(dat = truth$nums,
                                 time = timeall,
                                 species = spp.name,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# We need an effective N for the sample fish function, setting equal to actual N
# 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
numssshigh <- sample_fish(survey_ssN, effNhigh)

# aggregate true resn per survey design
aggresnss <- aggregateDensityData(dat = truth$resn,
                                 time = timeall,
                                 species = spp.name,
                                 boxes = boxall)

# aggregate true structn per survey design
aggstructnss <- aggregateDensityData(dat = truth$structn,
                                 time = timeall,
                                 species = spp.name,
                                 boxes = boxall)

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

resnss <-  sample_fish(aggresnss, effNhigh, sample = FALSE)


```

### Biological sampling: true lengths and weight-at-age

Then we calculate lengths with `calc_age2length`. This will probably not work for the full run with all species outputs 5x per year, but should be ok with 3 species. Still very long runtime! Saved output.

```{r lengthcomp-wtage, eval=FALSE}

length_census_ss <- calc_age2length(structn = structnss,
                                 resn = resnss,
                                 nums = numssshigh,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

#save for later use, takes a long time to generate
saveRDS(length_census_ss, file.path(d.name, paste0(scenario.name, "length_census_sardhakeboca.rds")))

```

### True biomass survey II

Now we can test the full `sample_survey_biomass` from nums to biomass index with these three species (required a revision of the function to allow time varying weight at age, and a clarification of units):

```{r compare-surveys}
# apply the proper weight at age (from saved census calc_age2length) to N survey

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

# weight at age output of calc_age2length is in g
# output of survey should be in t, sample_survey_biomass expects kg wt@age
# therefore need to divide by 1000 to get wt@age in kg

wtage <- length_census_ss$muweight %>%
  select(species, agecl, time, wtAtAge = atoutput) %>%
  mutate(wtAtAge = wtAtAge/1000)

surveyB_fromN <- sample_survey_biomass(survey_ssN, surv_cv, wtage)

```

Now we plot the results of the numbers-based (census) survey against the same Atlantis output file "[modelscenario]BiomIndx.txt" as above. We have a match. 

To get back true biomass we need the full detailed time varying weight at age output, which is fairly cumbersome. The fixed weight at age may be a decent approximation but does result in sigificant misses for some species [as seen here](https://sgaichas.github.io/poseidon-dev/CCAExamples.html). 

```{r matchB_Nsurveys }

atBtxt2tidy <- atBtxt2 %>%
  select(Time:DIN) %>%
  #select(Time, FPL:DIN) %>%
  rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, biomass, -Time) %>%
  filter(species %in% levels(surveyB_fromN$species))

#all species comparison, time intervals hardcoded for 5 steps per year
compareB_N <-ggplot() +
  geom_line(data=surveyB_fromN, aes(x=time/5,y=atoutput, color="survey census N->B"), 
            alpha = 10/10) +
  geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
             alpha = 1/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

compareB_N + 
  facet_wrap(~species, scales="free") 


```


## Test fishery functions: compare census with truth

Now for the fishery outputs. We will need total catch in tons for each species, and catch composition data (age and length). First we will test whether the catchbio output in `truth` matches the atlantis output file output[scenario.name]Catch.txt:

```{r testcatchbio1}

#just try the three species for now
catchbio_census_ss <- create_fishery_subset(dat = truth$catchbio,
                                         time = timeall,
                                         species = spp.name,
                                         boxes = boxall)

catchbio_census_ss_agg <- aggregate(atoutput ~ species + time,
                                    data=catchbio_census_ss, sum)

# learned the hard way this can be different from ecosystem outputs
fstepperyr <- if(runpar$outputstepunit=="days") 365/runpar$toutfinc

# what does the output look like?
catchB_outstep <- ggplot() +
  geom_line(data=catchbio_census_ss_agg, aes(x=time/fstepperyr,y=atoutput, color="catch atlantisom"),
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

catchB_outstep + 
  facet_wrap(~species, scales="free") 


```

It doesnt. A hack gets it closer, suggesting that the problem is a bug in the output for models using codebases prior to late 2015 (CCA is mid-2015). From Beth Fulton's note in the wiki:

>Note that in CATCH.NC the numbers at age had issues in old versions. For versions pre-dating December 2015 you need to divide the numbers at age by (86400 * number_water_column_layers_in_the_box). New versions post December 2015 are now in numbers without needing further adjustment.

```{r testcatchbio2}

# should it be aggregated? I just did this above
catchbio_series <- sample_survey_numbers(catchbio_census_ss, surv_cv)

# if it is the pre-2015 bug mentioned in atlantis wiki, this may get it close
catchbio_series_adjust <- catchbio_series %>%
  mutate(atoutput = atoutput/(86400*4))

#does it match catch.txt (which should be in tons)?
atCtxt <- read.table(file.path(d.name, catch.file), header=T)
atCtxt <- atCtxt[, -grep("TsAct", colnames(atCtxt))] #hardcoded for catch.txt output 

fishedlookup <- funct.groups %>%
  filter(IsFished > 0)

atCtxttidy <- atCtxt %>%
  rename_(.dots=with(fishedlookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, catchbio, -Time) %>%
  filter(species %in% levels(catchbio_series$species))

compare_catchB <-ggplot() +
  geom_line(data=catchbio_series_adjust, aes(x=time,y=atoutput, color="catch atlantisom"), 
            alpha = 10/10) +
  geom_point(data=atCtxttidy, aes(x=Time/365,y=catchbio, color="txt output catch bio"),
             alpha = 1/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

compare_catchB + 
  facet_wrap(~species, scales="free") 

```

Test this with NOBA model, which is in a newer codebase that doesn't have this problem with the CATCH.nc file. 

Update: NOBA has correct catch at age but the derived total catch biomass [still doesnt match catch.txt](https://sgaichas.github.io/poseidon-dev/TestCatchNOBA.html). I believe this is because we cannot multiply the correct weight at the layer the fishery catch is taken, which clearly does not match the median structn and resn across all model layers (as applied in `calc_biomass-age`. This may be buried deeply in the atlantis code, initialization, and prm files, but for now we will simply use the total catch weight out of atlantis catch.txt output.

## Catch output correction--apply to all species

Moving on to correcting catch at age for CCA:

```{r correctcatage}

# new atlantisom function to apply catage correction, or put inside run_truth

# first check that codebase is prior to corrected--to implement in atlantisom
#   read code version or date from log file
#   check against update date (Dec 12, 2015)
# grep("Atlantis SVN Last Change Date", readLines(file.path(d.name, "log.txt")), value = TRUE)

# alternative: less overhead if the system does it? log file is huge
# logfile <- paste0(file.path(d.name, "log.txt"))
# codedate <- system(paste0("grep 'Atlantis SVN' ", logfile))
# compare date in codedate in an if statement--fix in function

# read catch in nums output from truth$catch, or
# read in CATCH.nc (this from run_truth, modified with config files)

nc_catch <- paste0("output", scenario.name, 'CATCH.nc')
bboxes <- get_boundary(boxpars)

catch <- load_nc(dir = d.name,
                 file_nc = nc_catch,
                 bps = boxpars, #defined in census_spec.R
                 fgs = funct.groups,
                 select_groups = survspp,
                 select_variable = "Catch",
                 check_acronyms = TRUE,
                 bboxes = bboxes)

# find number of layers in each box
# The number of water column layers per box is in the initial conditions NC file:
# DIVinit_CalCurrentV3_Biol.nc : numlayers =

# read in initial conditions NC file
# this is already done for the other variables, modify from existing
# line 152 of load_nc
# num_layers <- RNetCDF::var.get.nc(ncfile = at_out, variable = "numlayers")[, 1]

# initial conditions file has different structure, so drop [,1]
at_init <- RNetCDF::open.nc(con = file.path(d.name, initial.conditions.file))

# Get info from netcdf file! (Filestructure and all variable names)
  var_names_ncdf <- sapply(seq_len(RNetCDF::file.inq.nc(at_init)$nvars - 1),
    function(x) RNetCDF::var.inq.nc(at_init, x)$name)
  numlayers <- RNetCDF::var.get.nc(ncfile = at_init, variable = "numlayers")
  
  RNetCDF::close.nc(at_init)
  
  # are these in box order??? if so make a box-numlayer lookup
  layerbox.lookup <- data.frame(polygon=boxall, numlayers)
  
  catch.tmp <- merge(catch, layerbox.lookup)

# divide the numbers at age by (86400 * number_water_column_layers_in_the_box)
# replace truth$catch atoutput with correction 
  catch <- catch.tmp %>%
    mutate(atoutput = atoutput / 86400 * numlayers) %>%
    select(species, agecl, polygon, time, atoutput)
  
  truth$catch <- catch
  
#save for later use, ar (better) resave truth file
#saveRDS(catch, file.path(d.name, paste0(scenario.name, "catchnums_corrected.rds")))
  save(truth,
      file = file.path(d.name, paste0("output", scenario.name, "run_truth.RData")))
```

Did the correction work? Not sure what to compare this to. First put it back through `calc_biomass_age` even though this won't match either?

```{r correctcatage-tobio, eval=FALSE}

biol <- load_biolprm(dir = d.name, file_biolprm = biol.prm.file)

# call within run_truth
catchbio <- calc_biomass_age(nums = catch,
    resn = truth$resn, structn = truth$structn, biolprm = biol)

#save for later use, even if wrong? takes a long time to generate
saveRDS(catchbio, file.path(d.name, paste0(scenario.name, "catchbio_ccatage.rds")))

```

It still doesn't match, similar to NOBA, but in the right ballpark now so the error is the median structn and resn called in the calc function above. This suggests that the catch at age in numbers is now correct, however.

```{r plotcorrected-bio}

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

catchbio_corrected_agg <- aggregate(atoutput ~ species + time,
                                    data=catchbio, sum)

compare_catchB <-ggplot() +
  geom_line(data=catchbio_corrected_agg, aes(x=time,y=atoutput, color="catch atlantisom"), 
            alpha = 10/10) +
  geom_point(data=atCtxttidy, aes(x=Time/365,y=catchbio, color="txt output catch bio"),
             alpha = 1/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

compare_catchB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 1, scales="free")
compare_catchB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 2, scales="free")
compare_catchB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 3, scales="free")
compare_catchB + 
  facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 4, scales="free")

```

## Fishery biological sampling

Now we can sample fishery catch for lengths, ages, weights:

Here we use `create_fishery_subset` on the numbers output of `run_truth` to create the survey census of age composition (for just our three 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 our three speces and the first 112 time steps.

```{r catchage-len-wt-3spp, echo=TRUE, warning=FALSE, message=FALSE}

# get survey nums with full (no) selectivity
catch_testNss <- create_fishery_subset(dat = truth$catch,
                                 time = c(0:111),
                                 species = spp.name,
                                 boxes = boxall)

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


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

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

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

catch_resnss <-  sample_fish(catch_aggresnss, effNhigh, sample = FALSE)

```

This is true catch at age for sardine:

```{r catage1}

catchage <- catch_numssshigh %>%
  filter(species == "Pacific_sardine")
  
catageplot <- ggplot(catchage, aes(x=agecl, y=atoutput)) +
  geom_point() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        catchage$species))

catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 1, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 2, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 3, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 4, scales="free")

```

This is true catch at age(class) for hake:

```{r catage2}

catchage <- catch_numssshigh %>%
  filter(species == "Mesopel_M_Fish")
  
catageplot <- ggplot(catchage, aes(x=agecl, y=atoutput)) +
  geom_point() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        catchage$species))

catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 1, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 2, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 3, scales="free")
catageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 4, scales="free")

```

Hake will need `calc_stage2age` applied because it has 2 true ages per age class.

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

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

catch_length_census_ss <- calc_age2length(structn = catch_structnss,
                                 resn = catch_resnss,
                                 nums = catch_numssshigh,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 200,
                                 CVlenage = 0.1, remove.zeroes=TRUE)


```

We should get the upper end of anything with a 200cm max length bin.

Sardine catch lengths:

```{r vis-fishery-lf-test-1}

catchlen <- catch_length_census_ss$natlength %>%
  filter(species == "Pacific_sardine")

lfplot <- ggplot(catchlen, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        catchlen$species))

lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 1, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 2, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 3, scales="free_y")

```

Hake catch lengths:

```{r vis-fishery-lf-test-2}

catchlen <- catch_length_census_ss$natlength %>%
  filter(species == "Mesopel_M_Fish")

lfplot <- ggplot(catchlen, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        catchlen$species))

lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 1, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 2, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 3, scales="free_y")

```

Fishery weight at (st)age:

```{r vis-fishery-wtageclass-test}

wageplot <- ggplot(catch_length_census_ss$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")

```

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

```{r aggwtcomp}

wtage_annsurv <- catch_length_census_ss$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 sample"))

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

```

## References