TestBiggerMaxLen.Rmd

---
title: "Testing atlantisom: allow users to specify different max length bin for length compositions"
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

Here we briefly test changes to the function `calc_age2length` that should allow the user to specify a different max length bin. 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. For full explanation of methods, see the file linked at the beginning of each section. 

```{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 4: sample for length composition (testing revised function)

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

We will apply examples here to only one species, Greenland halibut in NOBA, which grows to a large size.

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. 

```{r sppgroups, echo=TRUE}

# make defaults that return a standard survey, implement in standard_survey
# users need to map their species groups into these general ones
#   large pelagics/reef associated/burrowers/otherwise non-trawlable
#   pelagics
#   demersals
#   selected flatfish
if(initNOBA) funct.groups <- rename(funct.groups, GroupType = InvertType)

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

if(initCCA) { #Sarah's CCA Grouping
  nontrawl <- c("Shark_C","Yelloweye_rockfish","Benthopel_Fish","Pisciv_S_Fish",
                "Pisciv_T_Fish","Shark_D","Shark_P")
  pelagics <- c("Pisciv_V_Fish","Demersal_S_Fish","Pacific_Ocean_Perch","Mesopel_M_Fish",
                "Planktiv_L_Fish","Jack_mackerel","Planktiv_S_Fish","Pacific_sardine",
                "Anchovy","Herring","Pisciv_B_Fish")
  demersals <- c("Demersal_P_Fish","Planktiv_O_Fish","Demersal_D_Fish",
                 "Demersal_DC_Fish","Demersal_O_Fish","Darkblotched_rockfish",
                 "Demersal_F_Fish","Demersal_E_Fish","Bocaccio_rockfish",
                 "Demersal_B_Fish","Shark_R","Mesopel_N_Fish","Shark_B","Spiny_dogfish",
                 "SkateRay")
  selflats <- c("Pisciv_D_Fish", "Arrowtooth_flounder","Petrale_sole")
}

if(initNEUS) { # Sarah's NEUS Grouping
  nontrawl <- c("Pisciv_T_Fish", "Shark_D", "Shark_P", "Reptile", "Mesopel_M_Fish")
  pelagics <- c("Planktiv_L_Fish", "Planktiv_S_Fish", "Benthopel_Fish", "Pisciv_S_Fish")
  demersals <- c("Pisciv_D_Fish", "Demersal_D_Fish","Demersal_E_Fish", 
                 "Demersal_S_Fish","Demersal_B_Fish","Demersal_DC_Fish",
                 "Demersal_O_Fish","Demersal_F_Fish",
                 "Shark_B", "SkateRay")
  selflats <- c("Pisciv_B_Fish")
}

if(initNOBA) { # Sarah's NOBA Grouping
  nontrawl <- c("Sharks_other", "Pelagic_large","Mesop_fish")
  pelagics <- c("Pelagic_small","Redfish_other","Mackerel","Haddock",
                "Saithe","Redfish","Blue_whiting","Norwegian_ssh","Capelin")
  demersals <- c("Demersals_other","Demersal_large","Flatfish_other","Skates_rays",
                 "Green_halibut","North_atl_cod","Polar_cod","Snow_crab")
  selflats <- c("Long_rough_dab")
}

```

The following settings are for our example standard survey once per year, most areas, with mixed efficiency and selectivity: 

```{r stdbtsurvey-spec, message=FALSE, warning=FALSE, echo=TRUE}

# general specifications for bottom trawl survey, with items defined above commented out to avoid wasting time loading already loaded files:
#   once per year at mid year
# 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

midptyr <- round(median(seq(1,stepperyr)))

annualmidyear <- seq(midptyr, noutsteps, stepperyr)

#   ~75-80% of boxes (leave off deeper boxes?)
boxpars <- load_box(d.name, box.file)
boxsurv <- c(2:round(0.8*(boxpars$nbox - 1)))

#   define bottom trawl mixed efficiency
ef.nt <- 0.01 # for large pelagics, reef dwellers, others not in trawlable habitat
ef.pl <- 0.1  # for pelagics
ef.dm <- 0.7  # for demersals
ef.fl <- 1.1  # for selected flatfish

# bottom trawl survey efficiency specification by species group
effnontrawl <- data.frame(species=nontrawl, efficiency=rep(ef.nt,length(nontrawl)))
effpelagics <- data.frame(species=pelagics, efficiency=rep(ef.pl,length(pelagics)))
effdemersals <- data.frame(species=demersals, efficiency=rep(ef.dm,length(demersals)))
effselflats <- data.frame(species=selflats, efficiency=rep(ef.fl,length(selflats)))

efficmix <- bind_rows(effnontrawl, effpelagics, effdemersals, effselflats)

#   mixed selectivity (using 10 agecl for all species)
#     flat=1 for large pelagics, reef dwellers, others not in trawlable habitat
#     sigmoid 0 to 1 with 0.5 inflection at agecl 3 for pelagics, reaching 1 at agecl 5, flat top
#     sigmoid 0 to 1 with 0.5 inflection at agecl 5 for most demersals and flatfish, reaching 1 at agecl 7, flat top
#     dome shaped 0 to 1 at agecl 6&7 for selected demersals, falling off to 0.7 by agecl 10

sigmoid <- function(a,b,x) {
  1 / (1 + exp(-a-b*x))
}

# survey selectivity specification by species group
selnontrawl <- data.frame(species=rep(nontrawl, each=10),
                          agecl=rep(c(1:10),length(nontrawl)),
                          selex=rep(1.0,length(nontrawl)*10))
selpelagics <- data.frame(species=rep(pelagics, each=10),
                          agecl=rep(c(1:10),length(pelagics)),
                          selex=sigmoid(5,1,seq(-10,10,length.out=10)))
seldemersals <- data.frame(species=rep(demersals, each=10),
                          agecl=rep(c(1:10),length(demersals)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))
selselflats <- data.frame(species=rep(selflats, each=10),
                          agecl=rep(c(1:10),length(selflats)),
                          selex=sigmoid(1,1,seq(-10,10,length.out=10)))

selexmix <- bind_rows(selnontrawl, selpelagics, seldemersals, selselflats)

# use this constant 0 cv for testing
surv_cv_0 <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

#   define bottom trawl survey cv by group
cv.nt <- 1.0 # for large pelagics, reef dwellers, others not in trawlable habitat
cv.pl <- 0.5  # for pelagics
cv.dm <- 0.3  # for demersals
cv.fl <- 0.3  # for selected flatfish

# specify cv by species groups
surv_cv_nontrawl <- data.frame(species=nontrawl, cv=rep(cv.nt,length(nontrawl)))
surv_cv_pelagics <- data.frame(species=pelagics, cv=rep(cv.pl,length(pelagics)))
surv_cv_demersals <- data.frame(species=demersals, cv=rep(cv.dm,length(demersals)))
surv_cv_selflats <- data.frame(species=selflats, cv=rep(cv.fl,length(selflats)))

surv_cv_mix <- bind_rows(surv_cv_nontrawl, surv_cv_pelagics, surv_cv_demersals, surv_cv_selflats)

```

And the numbers to be sampled for lengths each year: 

```{r stdsurvey-lensamp, warning=FALSE, message=FALSE, echo=TRUE}

#   define n fish for biological sampling by group
#   this could easily be a vector or time series, constant here
ns.nt <- 25 # for large pelagics, reef dwellers, others not in trawlable habitat
ns.pl <- 1000  # for pelagics
ns.dm <- 1000  # for demersals
ns.fl <- 1000  # for selected flatfish

effNnontrawl <- data.frame(species=nontrawl, effN=rep(ns.nt,length(nontrawl)))
effNpelagics <- data.frame(species=pelagics, effN=rep(ns.pl,length(pelagics)))
effNdemersals <- data.frame(species=demersals, effN=rep(ns.dm,length(demersals)))
effNselflats <- data.frame(species=selflats, effN=rep(ns.fl,length(selflats)))

effNmix <- bind_rows(effNnontrawl, effNpelagics, effNdemersals, effNselflats)

```

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. 

```{r stdsurveyNbased-GHR, echo=TRUE}

ss.name <- funct.group.names[funct.group.names == "Green_halibut"]

#change back to flat selectivity to see full comp
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))

# get survey nums with full (no) selectivity
ss_survey_testNstd_nosel <- create_survey(dat = truth$nums,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv,
                                 effic = efficmix,
                                 selex = selex1)

# now sample fish nums from this
ss_numsstd_nosel <- sample_fish(ss_survey_testNstd_nosel, effNmix)

# aggregate true resn per survey design
survey_aggresnstd <- aggregateDensityData(dat = truth$resn,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv)

# aggregate true structn per survey design
survey_aggstructnstd <- aggregateDensityData(dat = truth$structn,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv)

#dont sample these, just aggregate them using median (effNmix does nothing)
ss_structnstd <- sample_fish(survey_aggstructnstd, effNmix, sample = FALSE)

ss_resnstd <-  sample_fish(survey_aggresnstd, effNmix, sample = FALSE)

```

Length sample with default max length bin (150 cm):

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

ss_length_stdsurv_nosel <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd_nosel,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)


```

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

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

ss_length_stdsurv_nosel_max <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd_nosel,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 250,
                                 CVlenage = 0.1, remove.zeroes=TRUE)


```

Plots show that default length bins are not adequate for something like Greenland halibut, which get larger than the default largest bin in `calc_age2length`, 150 cm. Both census and these sampled length comps are chopped off at 150 cm for this species:

```{r sslengthsamp2-testdefault}

lfplot <- ggplot(ss_length_stdsurv_nosel$natlength, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, ss_length_stdsurv_nosel$natlength$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")
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")

```

I changed the `calc_age2length` code to pass the upper bin of 150 as a default but allow the user to change it in the function call. This plot shows results from setting maxbin = 250 for Greenland Halibut in NOBA:

```{r sslengthsamp2-testmax250}

lfplot <- ggplot(ss_length_stdsurv_nosel_max$natlength, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        ss_length_stdsurv_nosel_max$natlength$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")
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")

```

Now we can go back and apply standard survey selectivity as in the [NOBA visualization](https://sgaichas.github.io/poseidon-dev/NOBAStdSurvLengthCompTest.html) and see how the length comps look:

```{r stdsurveyNbased-GHR2, echo=TRUE}

ss.name <- funct.group.names[funct.group.names == "Green_halibut"]

# get survey nums with standard selectivity
ss_survey_testNstd <- create_survey(dat = truth$nums,
                                 time = annualmidyear,
                                 species = ss.name,
                                 boxes = boxsurv,
                                 effic = efficmix,
                                 selex = selexmix)

# now sample fish nums from this
ss_numsstd <- sample_fish(ss_survey_testNstd, effNmix)

# structn and resn stuff is exactly the same because selectivity is irrelevant
```

Length sample from standard survey with user specified max length bin (250 cm):

```{r user-maxlen2, echo=TRUE}

ss_length_stdsurv_max <- calc_age2length(structn = ss_structnstd,
                                 resn = ss_resnstd,
                                 nums = ss_numsstd,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 250,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

```

```{r sslengthsamp3-testmax250}

lfplot <- ggplot(ss_length_stdsurv_max$natlength, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        ss_length_stdsurv_max$natlength$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")
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")

```

Bigger max bin is definitely more appropriate for this species.

I suppose we could make the maxbin a vector by species, but for now may be overkill.