TrueLengthCompTest.Rmd

---
title: "Testing atlantisom: generate census length comps"
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.

Initial testing was conducted by S. Gaichas using R scripts in the R folder of this repository that are titled "PoseidonTest_[whatwastested].R". Initial tests are expanded and documented in more detail in these pages. C. Stawitz improved and streamlined the setup and intialization sections. 

## Setup

First, you will want to set up libraries and install atlantisom if you haven't already. This document assumes atlantisom is already installed. For complete setup and initialization, please see [TrueBioTest](https://sgaichas.github.io/poseidon-dev/TrueBioTest.html). 

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)
library(atlantisom)
```


## Initialize input files and directories, read in "truth"

Abbreviated here; for a full explanation please see [TrueBioTest](https://sgaichas.github.io/poseidon-dev/TrueBioTest.html). This document assumes that `atlantisom::run_truth` has already completed and stored an .RData file in the atlantis output model directory.

```{r initialize}

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

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 4: sample for length composition

This section uses the `atlantisom::create_survey()` and `atlantisom::sample_fish()` to get a biological sample dataset. From that dataset, several other functions can be run:

 *  age samples from sample_ages
 *  length samples from sample_lengths (TO BE WRITTEN)
 *  weight samples from sample_weights (TO BE WRITTEN)
 *  diet samples from sample_diet
 
Atlantis outputs numbers by cohort (stage-age) and growth informtation, but does not output size of animals directly. The function we are testing here, `atlantisom::calc_age2length` converts numbers by cohort to a length composition. The [workflow originally envisioned](https://onedrive.live.com/?authkey=%21AFQkOoKRz64TLUw&cid=59547B4CB95EF108&id=59547B4CB95EF108%216291&parId=59547B4CB95EF108%216262&o=OneUp) was to create a survey, sample fish from the survey, then apply this function. We determined that it will not work better to create a length comp for the whole population, which could then be sampled. We rethought what we needed [here](https://sgaichas.github.io/poseidon-dev/RethinkSamplingFunctions.html) and test the new functions in this document. 

### This stuff is biolerplate from other docs, repeated here for completeness: 
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 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. 

```{r stdsurveyNbased}

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

# consider saving this interim step if it takes a long time go generate


```

Next, get true annual (cohort) age comp from this census survey based on run truth. (is there a standard Atlantis output I can compare this to as we did for biomass?)

```{r truecohortagecomp}
# what is true composition? need annual by species, use code from sample_fish
# do tidyly
dat2 <- survey_testNall %>%
  group_by(species, agecl, time) %>%
  summarize(numAtAge = sum(atoutput))

totN <- dat2 %>%
  group_by(species, time) %>%
  summarize(totN = sum(numAtAge))

dat2totN <- merge(dat2, totN)

# ageclcomp <- ggplot(dat2totN, aes(x=agecl, y=numAtAge/totN, col=time)) +
#   geom_point()
# 
# ageclcomp + facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 1, scales="free")
# ageclcomp + facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 2, scales="free")
# ageclcomp + facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 3, scales="free")
# ageclcomp + facet_wrap_paginate(~species, ncol=3, nrow = 3, page = 4, scales="free")

```

Then we use the `sample_fish` function with very high effN and compare to true annual age comp calculated above as a test, which matches as we have shown [here](https://sgaichas.github.io/poseidon-dev/TrueLengthAgeCompTest.html).

However. We need both numbers and weights to get the length comp, and ultimately the weight at age for use in assessments. 

After using functions without thinking first and then [rethinking sampling flow](https://sgaichas.github.io/poseidon-dev/RethinkSamplingFunctions.html), we will now test new functons to aggregate both densities and numbers appropriately for input into the biological sampling workflow. 

What we want is to subset the nums, resn, and structn components of run_truth using the same create_survey design selecting boxes, time, and species. We need only apply survey efficiency and selectivity to nums, however. So perhaps the new aggregateDensityData function can be applied to resn and structn?

```{r aggdens}

# aggregate true resn per survey design
survey_aggresnall <- aggregateDensityData(dat = truth$resn,
                                 time = timeall,
                                 species = survspp,
                                 boxes = boxall)

# aggregate true structn per survey design
survey_aggstructnall <- aggregateDensityData(dat = truth$structn,
                                 time = timeall,
                                 species = survspp,
                                 boxes = boxall)

```

Now we should have inputs to `sample_fish` on the same scale, and they need to be aggregated across boxes into a single biological sample for the whole survey. We are not maintaining spatial structure in sampling because it isn't used in most assessments.

To do the proper aggregation and not apply the multinomial sampling to the density data, I rewrote `sample_fish` to apply the median if `sample=FALSE` in the function call. Test this:

```{r censussurvey-calcage2length, warning=FALSE, message=FALSE}

# 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)

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

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

# now cut these down to a single species for testing
# this should still represent a census but with polygon and layer aggregated

atf_numsallhigh <- numsallhigh[numsallhigh$species == "Arrowtooth_flounder",]
atf_structnall <- structnall[structnall$species == "Arrowtooth_flounder",]
atf_resnall <- resnall[resnall$species == "Arrowtooth_flounder",]

atf_length_censussurvsamp <- calc_age2length(structn = atf_structnall,
                                 resn = atf_resnall,
                                 nums = atf_numsallhigh,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)


```

Maybe this will work, we can hope. It looks like it did...

```{r atflengthsamp1-test}

lfplot <- ggplot(atf_length_censussurvsamp$natlength, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte()

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")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 4, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 5, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 6, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 7, scales="free_y")

# try dir = "v" option for vertical lf comparisons

```

## Is this working??

OK, maybe we have something here. Next I need to see how long it takes to run for all species in the survey.

```{r allspplength, eval=FALSE}

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

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

```

Answer: this ran for 5.5 hours. But it was really only 5 hours 20 minutes that was the last code block! Saving the file is highly recommended.

## References