Introduction

This page documents initial testing of the atlantisom package in development at https://github.com/r4atlantis/atlantisom using three different Atlantis output datasets. Development of atlantisom began at the 2015 Atlantis Summit in Honolulu, Hawaii, USA. On this page we demonstrate use of atlantisom on both California Current (CCA) and Norwegian Barents Sea Atlantis (NOBA) output files to test function calc_Z calculating total mortality. This is used both to derive “natural mortality” for comparison with stock assessments and to split age classes into true ages using the function calc_stage2age.

NOBA has output files in true ages that we can compare results with. These files are in the folder NOBAwithAnnAgeOutput (not kept on github due to large filesize).

We can also compare both models with the output[scenario.name]Mort.txt file, which reports annual M and F. The user’s manual notes: >This file is currently only useful for looking at relative M vs F values for a species, as itdoes not give accurate mortalities

So we can assume that the scale of Z from this file will not be the same as our calculation, but it should be useful for splitting out F to get M, and for comparing overall patterns.

We will first use CCA here because the contrasting F scenario has been implemented. We will test with both sardines (M needed, but not age splitting) and hake (M and age splitting needed).

All model setup and configuration is described here.

library(tidyr)
require(dplyr)
library(ggplot2)
library(data.table)
library(here)
library(ggforce)
library(ggthemes)
library(atlantisom)
initCCA <- TRUE
initNEUS <- FALSE
initNOBA <- FALSE

species_ss <- c("Pacific_sardine")

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

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

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

if(initNOBA) source(here("config/NOBAaaConfig.R"))
#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
# 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"))))
}

Test calc Z function

This section tests the atlantisom::calc_Z() function step by step. Ultimately this output goes to calc_stage2age to produce annual numbers at true age, which we can compare with the ANNAGEBIO.nc output produced by NOBA. For now we will just compare to the Mort.txt output.

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

Output timestep toutinc for the population is 73,

so steps per year in run_truth output is 5

and the number of output steps in run_truth output is 500.

True numbers total (annual)

This gives us the true total numbers on an annual basis. We can compare this with age specific outputs.

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

# save this, needed as input to calc_Z?
saveRDS(survey_testNall, file.path(d.name, paste0(scenario.name, "survey_testNall.rds")))

# 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")))
surveyN <- readRDS(file.path(d.name, paste0(scenario.name, "surveyNcensus.rds")))

surveyN_ss <- surveyN[surveyN$species == species_ss,]

plotN <-ggplot() +
  geom_line(data=surveyN_ss, aes(x=time/stepperyr,y=atoutput, color="survey census N"), 
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

What level of aggregation works for calc_Z?

The atlantisom function calc_stage2age calls the calc_Z function, applying a total mortality rate to estimate the numbers at true age within an age class, but it is unclear whether to use this on the aggregated numbers or the full resolution numbers at age. Here we test outputs of Z for several levels of aggregation:

# add YOY file to the config files
YOY <- load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))

# load biolprm in some initialize file?
biol <- load_biolprm(d.name, biol.prm.file)
## Warning in `[<-.data.frame`(`*tmp*`, , 3:12, value = structure(list(`1` =
## structure(c(469L, : provided 13 variables to replace 10 variables
# get code matching species name to split YOY file
code_ss <- funct.groups$Code[which(funct.groups$Name == species_ss)]

# cut to a single species in YOY file
YOY_ss <- YOY %>%
  select(Time, paste0(code_ss, ".0"))

# numbers at agecl at full resolution (all polygons and layers)
truenums_ss <- truth$nums[truth$nums$species == species_ss,]

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

# numbers at agecl aggregated over layers, retain polygons (output of create_survey)
survey_testN_ss <- survey_testNall[survey_testNall$species == species_ss,] 

# numbers at agecl aggregated over layers and polygons
Natage_ss <- readRDS(file.path(d.name, paste0(scenario.name, "natage_census_sard.rds")))

#calc_Z <- function(yoy, nums, fgs, biolprm, toutinc)

fullresZ <- calc_Z(yoy = YOY_ss,
                   nums = truenums_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

surveyresZ <- calc_Z(yoy = YOY_ss,
                   nums = survey_testN_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

sampleresZ <- calc_Z(yoy = YOY_ss,
                   nums = Natage_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

Fixed a bug where grep was used if a # appeared. This ran for the survey resolution numbers. I then went back and ran the two other levels above after fixing calc_Z.

Both true full resolution numbers and survey resolution numbers (aggregated over polygon) returned nonzero Z values. As with previous the NOBA test, they are identical aside from eight values with no apparent pattern (rows 100, 105, 180, 200, 220, 240, 355, 400) out of 500 that differ at e-16. Same story with difference from sample resolution; this does not matter to output.

fullresZ$atoutput - surveyresZ$atoutput
##   [1]  0.000000e+00  0.000000e+00 -1.114560e-16  1.114560e-16  0.000000e+00
##   [6]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [11]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [16]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [21]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [26]  0.000000e+00  0.000000e+00  0.000000e+00  2.229120e-16  0.000000e+00
##  [31]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [36]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [41]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [46]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [51]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [56]  0.000000e+00 -1.118897e-16  1.249001e-16  0.000000e+00  0.000000e+00
##  [61]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [66]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [71]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [76]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [81]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [86]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [91]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [96]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  2.229120e-16
## [101] -1.114560e-16  0.000000e+00  0.000000e+00 -2.229120e-16  2.237793e-16
## [106]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [111]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [116]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [121]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [126]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [131]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [136]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [141]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [146]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [151]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [156]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [161]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [166]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [171]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [176]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  2.289835e-16
## [181] -1.179612e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [186]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [191]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [196]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -2.289835e-16
## [201]  2.289835e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [206]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [211]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [216]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -1.110223e-16
## [221]  2.289835e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [226]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [231]  0.000000e+00  0.000000e+00 -1.249001e-16  1.110223e-16  0.000000e+00
## [236]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  2.289835e-16
## [241] -2.289835e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [246]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [251]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [256]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -1.179612e-16
## [261]  1.144917e-16  0.000000e+00  0.000000e+00  0.000000e+00  1.110223e-16
## [266] -1.110223e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [271]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [276]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [281]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [286]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [291]  0.000000e+00  2.255141e-16 -2.359224e-16  0.000000e+00  0.000000e+00
## [296]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [301]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [306]  0.000000e+00 -2.237793e-16  1.179612e-16  0.000000e+00  0.000000e+00
## [311]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [316]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [321]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [326]  0.000000e+00  0.000000e+00  2.498002e-16 -1.144917e-16  0.000000e+00
## [331]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [336]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [341]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [346]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [351]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.127570e-16
## [356] -1.127570e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [361]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [366] -2.255141e-16  1.127570e-16  0.000000e+00  0.000000e+00  0.000000e+00
## [371]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [376]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [381]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [386]  0.000000e+00  0.000000e+00  0.000000e+00 -2.289835e-16  1.127570e-16
## [391]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [396]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.110223e-16
## [401] -2.289835e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [406]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [411]  1.110223e-16 -1.144917e-16  0.000000e+00  0.000000e+00  0.000000e+00
## [416]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [421]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [426]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [431]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [436]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [441] -2.237793e-16  4.510281e-16 -2.359224e-16  0.000000e+00  0.000000e+00
## [446]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [451]  0.000000e+00  0.000000e+00  0.000000e+00  1.110223e-16 -2.255141e-16
## [456]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [461]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [466]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [471]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [476]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [481]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [486]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [491]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [496]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [501]  0.000000e+00
fullresZ$atoutput - sampleresZ$atoutput
##   [1]  0.000000e+00  0.000000e+00 -1.114560e-16  1.114560e-16  0.000000e+00
##   [6]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [11]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [16]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00 -1.114560e-16
##  [21]  2.229120e-16  0.000000e+00  0.000000e+00  0.000000e+00  1.114560e-16
##  [26] -1.114560e-16  0.000000e+00  0.000000e+00  2.229120e-16 -1.110223e-16
##  [31]  3.348016e-16 -1.118897e-16  1.144917e-16 -1.110223e-16  0.000000e+00
##  [36]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [41]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [46]  0.000000e+00  0.000000e+00  0.000000e+00 -1.110223e-16  1.118897e-16
##  [51]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [56]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [61]  0.000000e+00  0.000000e+00  0.000000e+00  1.114560e-16 -1.114560e-16
##  [66]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.114560e-16
##  [71] -1.114560e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [76]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [81]  1.110223e-16 -1.110223e-16  0.000000e+00  0.000000e+00  0.000000e+00
##  [86]  0.000000e+00  1.114560e-16 -1.110223e-16  0.000000e+00  0.000000e+00
##  [91]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
##  [96]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  2.229120e-16
## [101] -1.114560e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [106]  0.000000e+00 -1.118897e-16  2.498002e-16  2.229120e-16 -1.118897e-16
## [111]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [116]  0.000000e+00  0.000000e+00  1.179612e-16 -1.118897e-16  0.000000e+00
## [121]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [126]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [131]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [136]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.118897e-16
## [141] -2.229120e-16  0.000000e+00 -2.498002e-16  2.229120e-16  0.000000e+00
## [146]  0.000000e+00 -2.229120e-16  1.144917e-16  0.000000e+00  0.000000e+00
## [151]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [156]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [161]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [166]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [171]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [176]  0.000000e+00  0.000000e+00  0.000000e+00 -1.249001e-16  2.289835e-16
## [181] -1.179612e-16  1.179612e-16  0.000000e+00  0.000000e+00  0.000000e+00
## [186] -2.289835e-16  2.289835e-16  0.000000e+00  1.110223e-16 -2.255141e-16
## [191]  0.000000e+00  1.110223e-16 -1.249001e-16  0.000000e+00  0.000000e+00
## [196]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [201] -1.144917e-16  2.289835e-16  0.000000e+00 -1.249001e-16  1.110223e-16
## [206]  0.000000e+00 -1.110223e-16  1.249001e-16  0.000000e+00  0.000000e+00
## [211]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [216]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [221]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [226]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [231]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [236]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  2.289835e-16
## [241] -2.289835e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [246]  0.000000e+00  0.000000e+00  0.000000e+00  2.498002e-16 -1.144917e-16
## [251]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [256]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [261]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.110223e-16
## [266] -1.110223e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [271] -2.237793e-16  2.255141e-16  0.000000e+00  0.000000e+00  0.000000e+00
## [276]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [281]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [286]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [291]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [296]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [301]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [306]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [311]  0.000000e+00  0.000000e+00  0.000000e+00  1.179612e-16 -2.255141e-16
## [316]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [321]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [326]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [331]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [336]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [341]  0.000000e+00  0.000000e+00  0.000000e+00 -1.179612e-16  0.000000e+00
## [346]  2.255141e-16  0.000000e+00  0.000000e+00  2.289835e-16 -2.255141e-16
## [351]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.127570e-16
## [356]  0.000000e+00 -1.127570e-16  0.000000e+00  0.000000e+00  0.000000e+00
## [361]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [366] -2.255141e-16  1.127570e-16  0.000000e+00  0.000000e+00  2.289835e-16
## [371] -1.127570e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [376]  0.000000e+00  0.000000e+00  0.000000e+00 -2.359224e-16  2.255141e-16
## [381]  0.000000e+00  1.127570e-16 -1.249001e-16  0.000000e+00  1.127570e-16
## [386] -1.127570e-16  0.000000e+00  0.000000e+00 -2.289835e-16  1.127570e-16
## [391]  0.000000e+00  0.000000e+00 -2.289835e-16  1.179612e-16  0.000000e+00
## [396]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  1.110223e-16
## [401] -2.289835e-16  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [406]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [411]  1.110223e-16 -1.144917e-16  2.498002e-16 -1.179612e-16  0.000000e+00
## [416]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [421]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [426]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [431]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [436]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [441]  0.000000e+00  2.255141e-16 -2.359224e-16  0.000000e+00  0.000000e+00
## [446]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [451]  0.000000e+00  0.000000e+00  0.000000e+00  1.110223e-16 -2.255141e-16
## [456]  0.000000e+00  0.000000e+00  1.110223e-16 -1.110223e-16  0.000000e+00
## [461]  0.000000e+00  0.000000e+00  0.000000e+00  2.498002e-16 -1.144917e-16
## [466]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [471]  0.000000e+00  0.000000e+00  1.249001e-16 -1.249001e-16  0.000000e+00
## [476]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [481]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [486]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
## [491] -1.144917e-16  1.179612e-16  0.000000e+00  0.000000e+00  0.000000e+00
## [496]  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  2.289835e-16
## [501]  2.289835e-16

I would be happer if the differences were all exactly 0, but I suppose we can attribute that to rounding error.

What is F at each timestep?

A relative annual F is found in outputMort.txt!

file.mort <- file.path(d.name, paste0("output", scenario.name, "Mort.txt"))

mortish <- read.table(file.mort, header = TRUE)

relF_ss <- mortish %>%
  select(Time, relF = paste0(code_ss, ".F"))

relM_ss <- mortish %>%
  select(Time, relM = paste0(code_ss, ".M"))

plotF <- ggplot(data=relF_ss, aes(x=Time/365, y=relF)) +
  geom_line() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotF

This is consistent with the scenario implemented. So far so good.

This probably depends on how the model is set up, but for models forced with F by species we ought to be able to read this out from an input file. If it is a daily rate then we mulitply by toutinc to get the rate at each of these output timestep Z values.

Quick test to see if Z is right?

Compare with outputMort.txt!

Zish <- merge(relF_ss,relM_ss) %>%
  mutate(relZ = relF + relM)

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ)) +
  geom_point(data=fullresZ, aes(x=(time*73)/365, y=atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ

After several iterations of debugging, we are still not matching the pattern, but improved in that it isn’t a constant punctuated by occasional high. Variation still doesn’t look anything like the txt file either. The 5x yr timesteps are probably messing this up?

I’ve commented out the detailed debugging block below, but check out the rmd file if you are interested in the sausage making. Changes have been implemented in calc_Z so that it now apportions annual recruitment to timesteps and recruitment is correctly aligned with truth nums output.

However, it is not clear that subtracting off “recruitment” that has not been subject to mortality is appropriate when the output numbers at each timestep have been subject to mortality? I guess that was the input in each timestep as best we know it, but this is at best an approximation.

How to turn timestep Z into annual Z for proper comparison?

Two approaches: annual nums each year, calculate survival between years (minus annual recruitment), turn into annual Z, or add timestep Zs to get annual?

First an annual estimate based on end year numbers snapshots:

totnums <- aggregate(atoutput ~ species + time, data = truenums_ss, sum) %>%
    mutate(time.days = (time+1)*runpar$toutinc) %>% #makes time 0 into days 0->73, etc
    mutate(yr = ceiling(time.days/365))  # yr 1 is 0:stepsperyr to match recruits yr1
    #mutate(yr = floor(time.days/365)) # includes 0 value, yr 1 is 

#only want numbers at end of year, snapshot
totnumsann <- totnums %>%
  filter(time.days %in% seq(365, max(time.days), by=365))

# mg carbon converted to wet weight in tonnes
k_wetdry <- biol$kgw2d / 1000000000

# WARNING only works for CCA because YOY.txt rows 2:end are already in numbers
# and we are merging out the incorrect and irrelevant YOY row 1 (Time=0)
# also hardcoded for sardine example
recnums <- YOY_ss %>%
  mutate(yr = as.integer(round(YOY_ss$Time)/365)) %>%
  mutate(recnums = SAR.0/k_wetdry/biol$redfieldcn)

annZ1 <- merge(totnumsann, recnums) %>%
  mutate(recnums = replace_na(recnums, 0)) %>%
  mutate(numslessrec = atoutput - recnums) %>%
  mutate(surv = numslessrec/lag(atoutput, default = first(numslessrec))) %>%
  mutate(Z = -1 * log(surv))

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ)) +
  geom_point(data=annZ1, aes(x=Time/365, y=Z)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(0, 2.0)

This looks better; we are starting to see a pattern.

Second summing within year Z from calc_Z:

annZ2 <- fullresZ %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(Z = sum(atoutput))

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ)) +
  geom_point(data=annZ2, aes(x=yr, y=Z)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(0, 2.0)

Also looks better. Not sure which is more correct.

Compare methods:

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_point(data=annZ1, aes(x=Time/365, y=Z, color="endyr Z")) +
  geom_point(data=annZ2, aes(x=yr, y=Z, color="sum timestep Z")) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(0, 2.0)

Both seem to be tracking each other reasonably, and also general trends from mort.txt output. Discussion with Beth suggests that mort.txt output is an initial calculation combining deaths due to predation with deaths due to fishing, but that there is rescaling of predation (and possibly fishing?) afterward this output stage such that mort.txt does not provide final mortality. If the rescaling varies by timestep, then matching general patterns rather than the full internannual variability is probably the best we can do, since we don’t know the full interannual variability.

We can probably conclude that the calc_Z function is working as correctly as possible now, given that it is an approximation no matter what, and the differences we see here are likley from temporal resolution of the Z estimate.

So what is M? This is an approximate annual M as Z-F from mort.txt (which assumes F is not rescaled!). If we need annual and age specific M then we need a new Z calc function that is cohort specific I think.

annM1 <- merge(annZ1, relF_ss, all.x = T) %>%
  mutate(M = Z-relF)

annZ2_Time <- annZ2 %>%
  mutate(Time = yr*365)

annM2 <- merge(annZ2_Time, relF_ss) %>%
  mutate(M = Z-relF)

plotM <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_line(data=Zish, aes(x=Time/365, y=relM, color="mort.txt M")) +
  geom_point(data=annM1, aes(x=Time/365, y=M, color="M = endyr Z-F")) +
  geom_point(data=annM2, aes(x=yr, y=M, color="M = sum timestep Z-F")) +
  theme_tufte() +
  theme(legend.position = "bottom", legend.box = "horizontal") +
  scale_color_discrete(NULL) +
  guides(color = guide_legend(nrow = 1)) +
  labs(subtitle = paste(scenario.name, species_ss))

plotM + ylim(0, 2.0)

So the problem here is we get negative M if we use our calculated Z from annual numbers at age or timestep numbers at age and subtract the F from mort.txt. (Negative M does not show up as points above because I restrict the y axis scale to 0-2.0.) The F from mort.txt is actually higher than our estimated Z. So Atlantis must rescale F too???

Need to investigate backing out F…

Let’s see how it performs with another model.

Test calc_Z NOBA

Might be a good time to switch back to NOBA for comparison with true annual ages.

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

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

species_ss <- "North_atl_cod"

#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

#Get true NOBAaa
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"))))
}

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

Lets try just the Z part for NOBA:

# make a function for this
# add YOY file to the config files
YOY <- load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))

# load biolprm in some initialize file?
biol <- load_biolprm(d.name, biol.prm.file)
## Warning in `[<-.data.frame`(`*tmp*`, , 3:12, value = structure(list(`1` =
## structure(c(124L, : provided 20 variables to replace 10 variables
# get code matching species name to split YOY file
code_ss <- funct.groups$Code[which(funct.groups$Name == species_ss)]

# cut to a single species in YOY file
YOY_ss <- YOY %>%
  select(Time, paste0(code_ss, ".0"))

# numbers at agecl at full resolution (all polygons and layers)
truenums_ss <- truth$nums[truth$nums$species == species_ss,]

#calc_Z <- function(yoy, nums, fgs, biolprm, toutinc)

#oops, need to generalize calc_Z for subannual timesteps in YOY!
#or input YOY only at 0, 365, etc since the numbers repeat for timesteps 
YOY_ss <- YOY_ss %>%
  filter(Time %in% seq(0, max(Time), by=365))

fullresZ <- calc_Z(yoy = YOY_ss,
                   nums = truenums_ss,
                   fgs = funct.groups,
                   biolprm = biol,
                   toutinc = runpar$toutinc)

# compare as above with mort.txt output
file.mort <- file.path(d.name, paste0("output", scenario.name, "Mort.txt"))

mortish <- read.table(file.mort, header = TRUE)

relF_ss <- mortish %>%
  select(Time, relF = paste0(code_ss, ".F"))

relM_ss <- mortish %>%
  select(Time, relM = paste0(code_ss, ".M"))

Zish <- merge(relF_ss,relM_ss) %>%
  mutate(relZ = relF + relM)

annZ2 <- fullresZ %>%
  mutate(yr = floor(time/stepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(Z = sum(atoutput))

plotZ <- ggplot() +
  geom_line(data=Zish, aes(x=Time/365, y=relZ, color="mort.txt Z")) +
  geom_point(data=fullresZ, aes(x=time/stepperyr, y=atoutput, color="each timestep Z")) +
  geom_point(data=annZ2, aes(x=yr-1, y=Z, color="sum timestep Z")) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name, species_ss))

plotZ + ylim(-1, 10.0)
## Warning: Removed 1 rows containing missing values (geom_point).

This Z is a bad match. The Z values in NOBA’s mort.txt are huge. Not sure what to make of this.

The above figure is now the best it is going to get after a deeper dive back into calc_Z for NOBA before we proceed to true ages. Found one bug hidden by sardine’s equal recruit_period and recend, another where the recruitment period is all within one output timestep that resulted in proportions over 1 (inflating the recruitment and then subtracting it off). Further, it appears that timesteps for recruitment were still not correctly aligned, so trying a different alignment that leaves out the timestep 0 from the truth. Not sure the latter is correct, or if we need to rethink the calculation of survival based on when recruits actually enter the model.

Update: the new code actually makes the CCA sardine Z estimates more in line with each other, so the updates are correct and we are keeping them in calc_Z.

Most important: Beth says that the M in mort.txt can really be as far off as we are seeing it here. The green line in the above NOBA Z plot reflects the initial estimate of what dies by predation mortality, which is a gigantic overesimate because many availability parameters are set to 1. After this output step, the predation mortality is rescaled within the model to account for the fact that there really aren’t that many to be eaten.

Our estimate of Z is from the actual numbers alive at the end of each timestep and the known input recruitment prior to any mortality, so I think it really is the best estimate we have. More debugging code commented out below but see the rmd file for the gory and redundant details.

Biological sampling: true numbers at age class

Number of true age classes is stored in the […]groups[…].csv file, which is read in using load_fgs. We read this in as funct.groups above. The column NumAgeClassSize stores the number of true ages:

Number of age classes for Cod: 2

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

Natage_ss <- Natage[Natage$species == species_ss,]

Natageplot <- ggplot(Natage_ss, aes(x=agecl, y=atoutput)) +
  geom_point() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        Natage_ss$species))

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

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

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

Natageplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 4, scales="free")

This really is still overly optimistic, but now that I think calc_Z is working, if we feed these into calc_stage2age and get a result that matches the NOBA output ANNAGEBIO.nc, then we can assume that function is correct.

These are estimated true age comps for Cod in annual ages:

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

# numbers at ageclaggregated over layers, retain polygons (output of create_survey)
survey_testN_ss <- survey_testNall[survey_testNall$species == species_ss,]


# for single species, Z was calculated above as surveyresZ
trueages_ss <- calc_stage2age(nums_data = survey_testN_ss,
                                 biolprm = biol, 
                                 yoy = YOY_ss,
                                 fgs = funct.groups, 
                                 runprm = runpar)

# run sample fish on this for completeness? yes to aggregate
trueNagesamp_ss <- sample_fish(trueages_ss, effNhigh)

# plot it

Natageplot <- ggplot(trueNagesamp_ss, aes(x=agecl, y=atoutput)) +
  geom_point() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        trueNagesamp_ss$species))

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

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

Natageplot + facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 3, scales="free")

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

Still not looking great, but at least not constant. I am now assuming we have the wrong calculations in calc_stage2age. Load true ages from ANNAGEBIO.nc:

annage_nc <- paste0("output", scenario.name, "ANNAGEBIO.nc")
bboxes <- get_boundary(boxpars)

# need to change load_nc, at line 117 hardcoded for 10 cohorts!
trueage_atout <- load_nc(dir = d.name,
                         file_nc = annage_nc,
                         bps = boxpars,
                         fgs = funct.groups,
                         select_groups = survspp,
                         select_variable = "Nums",
                         check_acronyms = TRUE,
                         bboxes = bboxes)
## [1] "Read /Users/sgaichas/Documents/0_Data/Atlantis/Poseidon/poseidon-dev/atlantisoutput/NOBAwithAnnAgeOutput/outputnordic_runresults_01ANNAGEBIO.nc successfully"
trueage_atout_ss <- trueage_atout[trueage_atout$species == species_ss,]

#survey census to agg over layers
survey_trueage_atout_ss <- create_survey(dat = trueage_atout_ss,
                                 time = timeall,
                                 species = survspp,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# sample fish to agg over polygons
sample_trueage_atout_ss <- sample_fish(survey_trueage_atout_ss, effNhigh)

Compare true age output from atlantis with our estimate (and with aggregated age comps):

compareAge <- ggplot() +
  geom_point(data=trueNagesamp_ss, aes(x=agecl,y=atoutput, color="atlantisom"), alpha = 10/10) +
  geom_point(data=sample_trueage_atout_ss, aes(x=agecl,y=atoutput, color="atout true"), alpha = 10/10) + 
  geom_point(data=Natage_ss, aes(x=agecl*2, y=atoutput, color="cohort true"), alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(subtitle = paste(scenario.name,
                        trueNagesamp_ss$species))

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

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

compareAge + 
  facet_wrap_paginate(~time, ncol=3, nrow = 3, page = 3, scales="free")

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

We only have the first 10 cohorts of true age output due to hardcoding in load_nc. I may write another load function that changes that. Either way our early ages don’t match. Need to step through calc_stage2age next!