NOBAcod.Rmd

---
title: "Generate data for NOBA cod ensemble modeling"
author: "Sarah Gaichas"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    code_fold: hide
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(here)
```

## Standard atlantisom output test (age class not age)

Let's use the wrapper functions to get initial stock assessment model inputs.

Then we need to adjust the base atlantisom functions for true age classes, and develop options to call those within the wrappers.

We need config files for the NOBA model, our surveys and fisheries. Surveys and fisheries are based on the ones we used for sardine in the California Current as a test. 

**Special note:** *Atlantis output files formerly had the prefix "output" and `atlantisom` was coded to expect this. NOBA output files (and newer Atlantis models) do not universally use the "output" prefix naming convention. Annage branch `atlantisom` functions have been updated for user-specified output filenames.*

**Special note 2:** *NOBA's nordic_groups_v04.csv file has an inFishery column, rather than the `atlantisom` expected InFishery column. I've changed the column name in the input file, but need to add robustness to atlantisom code to remove case sensitivity.*

### Config files

`NOBA2config.R` looks like this (adjusted from Alfonso's original):

```{r, code = readLines("./config/NOBA2Config.R"), eval=F}
```

`omdimensions.R` standardizes timesteps, etc.:

```{r, code = readLines("./config/omdimensions.R"), eval=F}
```

`codsurvey.R` configures the fishery independent survey and looks like this:

```{r, code = readLines("./config/codsurvey.R"), eval=F}
```

`codfishery.R` configures the fishery dependent data and looks like this:

```{r, code = readLines("./config/codfishery.R"), eval=F}
```

### Using the wrappers (now updated for annual age)

All the config files go in the config folder and we'll try running this:

```{r full-test, message=FALSE, warning=FALSE}

#library(here) now above
library(tidyverse)
library(atlantisom)
library(ggthemes)

NOBAom <- om_init(here("config/NOBA2config.R"))

NOBAom_cod <- om_species(c("North_atl_cod"), NOBAom)

NOBAom_cod_ind <- om_index(usersurvey = here("config/codsurvey.R"), 
                           userfishery = here("config/codfishery.R"),
                           omlist_ss = NOBAom_cod, 
                           n_reps = 1, 
                           save = TRUE)

NOBAom_cod_comp <- om_comps(usersurvey = here("config/codsurvey.R"), 
                           userfishery = here("config/codfishery.R"),
                           omlist_ss = NOBAom_cod, 
                           n_reps = 1, 
                           save = TRUE)


```

### Wrapper test results {.tabset}

#### Biomass index
```{r bioind-1}

omlist_ss <- NOBAom_cod
source(here("config/omdimensions.R"))

#read time series data
survObsBiom <- readRDS(file.path(d.name, paste0(scenario.name, "surveyB.rds")))

survObsBiom <- survObsBiom[[1]]

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

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

```


#### Catch time series
```{r catchind-1}
#read time series data
catchbio_ss <- readRDS(file.path(d.name, paste0(scenario.name, "fishCatch.rds")))

catchbio_ss <- catchbio_ss[[1]]

plotC <-ggplot() +
  geom_line(data=catchbio_ss, aes(x=time/365,y=atoutput, color="observed Catch"), 
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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


```


#### Survey length composition
```{r lencomp-1}

#length comps
len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsLenComp.rds")))
fish_len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))

len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]

#add this to om_indices function so that this has years when read in
fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))

len <- filter(len_comp_data, time %in% c(55:175))

  lfplot <- ggplot(len, aes(upper.bins)) +
    geom_bar(aes(weight = atoutput)) +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          len$species))
  
  lfplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")


```


#### Survey age composition (age classes)
```{r agecomp-1}

#read in comp data
age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsAgeComp.rds")))
age_comp_data <- age_comp_data[[1]]

Natage <- filter(age_comp_data, time %in% c(150:270))

Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
    geom_point() +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          Natage$species))
  
  Natageplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
```


#### Survey weight at age (age classes)
```{r wtageclass-1}

wtage <- readRDS(file.path(d.name, paste0(scenario.name, "survObsWtAtAge.rds")))
wtage <- wtage[[1]]

wageplot <- ggplot(wtage, 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")
```


#### Fishery age composition (age classes)
```{r fishagecomp-1}

#read in comp data
fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]

#add this to om_indices function so that this has years when read in
fish_age_comp$time <- fish_age_comp$time/fstepperyr

Natage <- filter(fish_age_comp, time %in% c(30:53))

Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
    geom_point() +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          Natage$species))
  
  Natageplot + facet_wrap(~time, ncol=6, scales="free_y")
```

### {-}

Wrappers still work for NOBA cod using standard age classes. 

Visualize new annual age output for numbers at age, catch at age, and weight at age (interpolated): 

### Wrapper test results for annual ages {.tabset #annageoutput}

#### Survey age composition (annual ages)
```{r plot-annagecomp}
#read in comp data
annage_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
annage_comp_data <- annage_comp_data[[1]]

#fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
#fish_age_comp <- fish_age_comp[[1]]

#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- as.integer(floor(fish_age_comp$time/fstepperyr))

Natage <- filter(annage_comp_data, time %in% c(150:270))

Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
    geom_point() +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          Natage$species))
  
  Natageplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
```

#### Survey iterpolated weight at age (annual ages) 
```{r plot-annwtage}

wtage <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullWtAtAge.rds")))
wtage <- wtage[[1]] #this still has second list component, diagnostic plot
wtage <- wtage[[1]]

wageplot <- ggplot(wtage, 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")

```

#### Fishery catch at age (annual ages)
```{r plot-catchagecomp}

fish_annage_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_annage_comp <- fish_annage_comp[[1]]

#add this to om_indices function so that this has years when read in
fish_annage_comp$time <- fish_annage_comp$time/fstepperyr

Catage <- filter(fish_annage_comp, time %in% c(30:53))

Catageplot <- ggplot(Catage, aes(x=agecl, y=atoutput)) +
    geom_point() +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          Catage$species))
  
  Catageplot + facet_wrap(~time, ncol=6, scales="free_y")
```

#### Fishery iterpolated weight at age (annual ages) 
```{r plot-fishannwtage}

wtage <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullWtAtAge.rds")))
wtage <- wtage[[1]] #this still has second list component, diagnostic plot
wtage <- wtage[[1]]

wageplot <- ggplot(wtage, 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")

```

### {-}

Note: NOBA doesn't seem to have a selective fishery for cod. Lots of age 1 in the catch!

# (Annage development details below)

Now to modify the underlying code to get full age structure.

## Updating functions for full age structure

General strategy:

1. add annual age option to run_truth function call, default it to FALSE
  + if true, looks for ANNAGEBIO and ANNAGECATCH in designated directory, error out if missing
  + if true, read these files in addition to existing set rather than replace existing
    + read them in with new load_nc_annage function
    + add them to the output object in addition to existing set
1. write load_nc_annage
  + need to understand annage nc structures
    + if same as others just with full age structure, replace previous nc
    + only numbers output? 
      + how to match with resn and structn for lengths/wts?
      + reconcile with bio parameters such as maturity ogive?
    + only for some groups?

What is in the ANNAGEBIO and ANNAGECATCH .nc files?

```{r annagebio-contents}

# taken from existing load_nc

file.nc <- file.path(d.name, "outputnordic_runresults_01ANNAGEBIO.nc")

  # Load ATLANTIS output!
  at_out <- RNetCDF::open.nc(con = file.nc)

  # Get info from netcdf file! (Filestructure and all variable names)
  var_names_ncdf <- sapply(seq_len(RNetCDF::file.inq.nc(at_out)$nvars - 1),
    function(x) RNetCDF::var.inq.nc(at_out, x)$name)
  n_timesteps <- RNetCDF::dim.inq.nc(at_out, 0)$length
  n_boxes     <- RNetCDF::dim.inq.nc(at_out, 1)$length
  n_layers    <- RNetCDF::dim.inq.nc(at_out, 2)$length

RNetCDF::close.nc(at_out)


```

ANNAGEBIO.nc variable names (first hundred): `r var_names_ncdf[1:100]` 

ANNABEBIO.nc variable names (last hundred): `r var_names_ncdf[1590:1690]`

ANNAGEBIO.nc timesteps: `r n_timesteps`

ANNAGEBIO.nc boxes: `r n_boxes`

ANNAGEBIO.nc layers: `r n_layers`

Annual age output is in numbers and weight at each timestep, layer, and box, which is really good news. We should be able to get weight at true age more directly this way than the interpolation I was planning. Also, we may be able to do length estimation differently? 

To modify `load_nc_annage` we need to figure out how many true ages each species has from the biological parameter file using `load_biolprm`:

```{r biolprm}

biolprm <- load_biolprm(d.name, biol.prm.file)
fgs <- load_fgs(d.name, functional.groups.file)

#max age should be ages per cohort times the number of age classes (<10 for some NOBA spp)

names(biolprm$agespercohort) <- c("code", "ageperagecl")

maxage <- merge(biolprm$maturityogive, biolprm$agespercohort, all.x = T) %>%
    select(code, nagecl, ageperagecl) %>%
    mutate(maxage = nagecl * ageperagecl) %>%
    mutate(name = fgs$Name[match(code, fgs$Code)])

select_groups <- c('Mackerel', 'North_atl_cod', 'Capelin')

for (i in seq_along(select_groups)) {
  select_group_ages <- 1:maxage$maxage[match(select_groups[i], maxage$name)]
}

```

Try this `load_nc_annage` function:

```{r loadncannage, eval=FALSE}

load_nc_annage <- function(dir = getwd(), file_nc, bps, fgs, biolprm, select_groups,
  select_variable =
  c("Nums", "Weight"),
  check_acronyms = TRUE, bboxes = c(0), verbose = FALSE) {
  # NOTE: The extraction procedure may look a bit complex...
  # A different approach would be to
  # create a dataframe for each variable (e.g. GroupAge_Nums)
  # and combine all dataframes at the end.
  # This alternative approach requires a lot more storage
  # and the code wouldn't be vectorized.

  # Check input of the nc file
  if (tail(strsplit(file_nc, "\\.")[[1]], 1) != "nc") {
    stop("The argument for file_nc,", file_nc, "does not end in nc")
  }
  if (is.null(dir)) {
    file.nc <- file_nc
  } else {
    file.nc <- file.path(dir, file_nc)
  }

  # Check input of select_variable as only one value is allowed
  select_variable <- match.arg(select_variable, several.ok = FALSE)

  # Check input structure!
  if (check_acronyms) {
    active_groups <- as.vector(subset(fgs, IsTurnedOn == 1)$Name)
    inactive_groups <- select_groups[which(
      !is.element(select_groups, active_groups))]
    if (length(inactive_groups) >= 1) {
      select_groups <- select_groups[!is.element(select_groups, inactive_groups)]
      warning(paste(paste("Some selected groups are not active in the model run.",
        "Check 'IsTurnedOn' in fgs\n"),
        paste(inactive_groups, collapse = "\n")))
    }
    if (all(!is.element(select_groups, active_groups))) {
      stop(paste("None of the species selected are active in the model run.",
        "Check spelling and Check 'IsTurnedOn' in fgs"))
    }
  }

  # Deal with file structures

  # Load ATLANTIS output!
  at_out <- RNetCDF::open.nc(con = file.nc)
  on.exit(RNetCDF::close.nc(at_out))

  if (select_variable != "N" & all(is.element(select_groups, bps))) {
    stop("The only output for Biomasspools is N.")
  } else{
    print(paste("Read", file.nc, "successfully"))
  }

  # Get info from netcdf file! (Filestructure and all variable names)
  var_names_ncdf <- sapply(seq_len(RNetCDF::file.inq.nc(at_out)$nvars - 1),
    function(x) RNetCDF::var.inq.nc(at_out, x)$name)
  n_timesteps <- RNetCDF::dim.inq.nc(at_out, 0)$length
  n_boxes     <- RNetCDF::dim.inq.nc(at_out, 1)$length
  n_layers    <- RNetCDF::dim.inq.nc(at_out, 2)$length

  # Extract data from the ncdf file
  # Create a vector of all potential variable names
  # Only use names which are available in the ncdf-file as an
  # extraction of missing variables is not possible
  # Create vector of available species at the end using search_clean
  # This is needed to create species-names later on


  # CHANGED hardcoded cohort number to reflect true ages from


  # To make the creation of variables as robust as possible
  # we introduce different combinations of groups, variable, and cohort
  # Only combinations present in the ncdf are used later on
  # Loop over select_groups to use the same ordering

  names(biolprm$agespercohort) <- c("code", "ageperagecl")

  maxage <- merge(biolprm$maturityogive, biolprm$agespercohort, all.x = T) %>%
    select(code, nagecl, ageperagecl) %>%
    mutate(maxage = nagecl * ageperagecl) %>%
    mutate(name = fgs$Name[match(code, fgs$Code)])

  search <- list()
  for (i in seq_along(select_groups)) {
    select_group_ages <- 1:maxage$maxage[match(select_groups[i], maxage$name)]
    search[[i]] <- c(
      unlist(lapply(paste0(select_groups[i], select_group_ages),
                    paste0, select_variable)),           # GroupCohortVariable
      unlist(lapply(paste0(select_groups[i], select_variable),
                    paste0, select_group_ages)),                   # GroupVariableCohort
      unlist(lapply(paste0(select_groups[i], select_group_ages),
                    paste, select_variable, sep = "_")), # GroupCohort_Variable
      unlist(lapply(paste(select_groups[i], select_variable, sep = "_"),
                    paste0, select_group_ages)),                   # Group_VariableCohort
      unlist(lapply(paste(select_groups[i], select_group_ages, sep = "_"),
                    paste, select_variable, sep = "_")), # Group_Cohort_Variable
      unlist(lapply(paste(select_groups[i], select_variable, sep = "_"),
                    paste, select_group_ages, sep = "_")),         # Group_Variable_Cohort
      unlist(paste0(select_groups[i], select_variable)),                                                      # GroupVariable
      unlist(paste(select_groups[i], select_variable,
                   sep = "_"))                           # Group_Variable
      )
    search[[i]] <- search[[i]][is.element(search[[i]], var_names_ncdf)]
    search[[i]] <- unique(search[[i]])
  }
  search_clean <- do.call(c, search)
  # If the combination of select_groups and select_variable ends up not being found.
  if (length(search_clean) == 0) return(0)

  at_data <- lapply(search_clean, RNetCDF::var.get.nc, ncfile = at_out)

  # Get final species and number of ageclasses per species
  final_species <- select_groups[sapply(
    lapply(select_groups, grepl, x = search_clean), any)]
  final_agecl <- maxage$maxage[
    sapply(final_species, function(x) which(x == maxage$name))]

  num_layers <- RNetCDF::var.get.nc(ncfile = at_out, variable = "numlayers")[, 1]
  # add sediment layer!
  num_layers <- num_layers + ifelse(num_layers == 0, 0, 1)

  # Create an array of layerids.
  # Every entry in the array indicates if a layer is present (= 1) or not (= 0).
  # Boxes without layers (= islands) have only 0s as id,
  # used later on to remove data from non-existent layers!
  # By default output should be 0 in those layers.
  # Layers in boundary boxes are set to 0 if bboxes is anything other than NULL!
  # Applying a boolean array to an array results in a vector!
  for (i in seq_along(num_layers)) {
    if (i == 1) layerid <- array(dim = c(n_layers, n_boxes))
    if (num_layers[i] == 0) {
      layerid[, i] <- 0
    } else {
      if (!is.null(bboxes) & is.element((i - 1), bboxes)) {
        layerid[, i] <- 0
      } else {
        layerid[, i] <- c(rep(0, times = n_layers - num_layers[i]),
                          rep(1, times = num_layers[i]))
      }
    }
  }

  # Create vectors for polygons and layers
  # Each vector has the length equal to one time-step
  # All data from islands and non-existent layers is removed
  # Therefore the length of these
  # vectors is equal for each extracted variable
  boxes <- 0:(n_boxes - 1)
  # Remove islands and boundary boxes
  island_ids <- num_layers == 0
  if (!is.null(bboxes)) {
    boundary_ids <- is.element(boxes, bboxes)
    island_ids <- island_ids | boundary_ids
  }
  boxes <- boxes[!island_ids]
  num_layers <- num_layers[!island_ids]

  polygons <- rep(boxes, times = num_layers)
  layers <- sapply(num_layers[num_layers != 0] - 2,
    function(x) c(seq(x, from = 0, by = 1), n_layers - 1))
  if (any(sapply(layers, length) != num_layers[num_layers != 0])) {
    stop("Number of layers incorrect.")
  }
  layers <- do.call(c, layers)
  if (length(polygons) != length(layers)) {
    stop("Number of polygons and layers do not match.")
  }

  # In the following section the data is transformed to a long dataframe
  # I haven't found any solution to vectorize the creation of the dataframe
  # columns (species, age, polygons,...)
  # when data from 2d and 3d arrays
  # (e.g. select_variable = "N" all biomasspools are only present in the
  # sediment layer.) are read in simultaneously.
  # Therefore the current "messy" solution splits the data
  # in 2 subpopulations: 2d-data and 3d-data
  at_data3d <- at_data[which(sapply(at_data, function(x) length(dim(x))) == 3)]
  at_data2d <- at_data[which(sapply(at_data, function(x) length(dim(x))) == 2)]

  int_fs <- final_species
  int_fa <- final_agecl

  if (length(at_data3d) >= 1) {
    # Remove biomasspools if selected variable is "N"!
    if (select_variable == "N") {
      int_fs <- final_species[!is.element(final_species, bps)]
      int_fa <- final_agecl[!is.element(final_species, bps)]
      # Note this only works if age-structured vertebrates have 10 ageclasses
      int_fa[int_fa == 10] <- 1
    }
    for (i in seq_along(at_data3d)) {# for loop over all variables
      if (i == 1) result3d <- list()
      for (j in 1:n_timesteps) {# loop over timesteps
        if (j == 1) values <- array(dim = c(length(layers), n_timesteps))
        values[, j] <- at_data3d[[i]][,, j][layerid == 1]
      }
      result3d[[i]] <- as.vector(values)
    }
    result3d <- data.frame(
      species = unlist(sapply(
        X = mapply(FUN = rep, x = int_fs, each = int_fa, SIMPLIFY = FALSE,
        USE.NAMES = FALSE),
        FUN = rep, each = length(layers) * n_timesteps, simplify = FALSE)),
      agecl = unlist(sapply(
        X = sapply(X = int_fa, FUN = seq, from = 1, by = 1, simplify = FALSE,
        USE.NAMES = FALSE),
        FUN = rep, each = length(layers) * n_timesteps, simplify = FALSE)),
      polygon = unlist(sapply(
        X = n_timesteps * int_fa, FUN = rep, x = polygons, simplify = F,
        USE.NAMES = FALSE)),
      layer = unlist(sapply(
        X = n_timesteps * int_fa, FUN = rep, x = layers, simplify = FALSE,
        USE.NAMES = FALSE)),
      time = unlist(sapply(
        X = int_fa, FUN = rep, x = rep(0:(n_timesteps - 1), each = length(layers)),
        simplify = FALSE, USE.NAMES = FALSE)),
      atoutput = do.call(c, result3d),
      stringsAsFactors = FALSE)
  }

  if (length(at_data2d) >= 1) {
    # Only select biomasspools if selected variable is "N"!
    if (select_variable == "N") {
      int_fs <- final_species[is.element(final_species, bps)]
      int_fa <- final_agecl[is.element(final_species, bps)]
    }
    # age-structured invert groups are combined in ncdf file!
    if (select_variable == "Grazing") int_fa <- 1
    for (i in seq_along(at_data2d)) {# for loop over all variables!
      if (i == 1) result2d <- list()
      for (j in 1:n_timesteps) {# loop over timesteps
        if (j == 1) values <- array(dim = c(length(boxes), n_timesteps))
        values[, j] <- at_data2d[[i]][, j][boxes + 1]
      }
      result2d[[i]] <- as.vector(values)
    }

    # Order of the data in value column = "atoutput".
    # 1. species  --> rep each with the number of
    #                 ageclasses and n_timesteps * boxes
    # 2. age      --> rep each (1:10 for each species) with n_timesteps * boxes
    # 3. timestep --> rep each timestep (1:n_timesteps)
    #                 with the number of boxes and final_agecl
    #                 (num cohorts per species)
    # 4. polygon  --> rep boxes times n_timesteps * final_agecl
    #                 (num cohorts per species)
    result2d <- data.frame(species = unlist(sapply(
      X = mapply(FUN = rep, x = int_fs, each = int_fa, SIMPLIFY = FALSE,
        USE.NAMES = FALSE),
        FUN = rep, each = length(boxes) * n_timesteps, simplify = FALSE)),
      agecl = unlist(sapply(X = sapply(X = int_fa, FUN = seq, from = 1,
        by = 1, simplify = FALSE, USE.NAMES = FALSE),
        FUN = rep, each = length(boxes) * n_timesteps, simplify = FALSE)),
      polygon = unlist(sapply(X = n_timesteps * int_fa,
        FUN = rep, x = boxes, simplify = FALSE, USE.NAMES = FALSE)),
      time = unlist(sapply(X = int_fa, FUN = rep, x = rep(0:(n_timesteps - 1),
        each = length(boxes)), simplify = FALSE, USE.NAMES = FALSE)),
      atoutput = do.call(c, result2d),
      stringsAsFactors = F)
    if (select_variable == "N") result2d$layer <- n_layers - 1
  }

  # Combine dataframes if necessary!
  if (all(sapply(lapply(at_data, dim), length) == 3) & select_variable != "N") {
    result <- result3d
  }
  if (all(sapply(lapply(at_data, dim), length) == 2) & select_variable != "N") {
    result <- result2d
  }
  if (select_variable == "N") {
    if (length(at_data2d) >= 1 & length(at_data3d) == 0) result <- result2d
    if (length(at_data2d) == 0 & length(at_data3d) >= 1) result <- result3d
    if (length(at_data2d) >= 1 & length(at_data3d) >= 1) {
      result <- rbind(result2d, result3d)
    }
  }

  # Remove min_pools if existent (well, there always are min pools... ;)).
  min_pools <- is.element(result$atoutput, c(0, 1e-08, 1e-16))
  if (length(min_pools) > 0) {
    # exclude 1st timestep and sediment layer from calculation
    print_min_pools <- sum(min_pools) -
      length(result[min_pools & result$time == 1, 1]) -
      length(result[min_pools & result$time > 1 & result$layer == 7, 1])
    if (print_min_pools > 0 & verbose) {
      warning(paste0(round(print_min_pools/dim(result)[1] * 100),
        "% of ", select_variable, " are true min-pools (0, 1e-08, 1e-16)"))
    }
    result <- result[!min_pools, ]
  }

  # Remove non-existent layers.
  # WARNING: Biomass is build up (very few) in sediment layer for
  # NON sediment groups (e.g. baleen whales)
  # Therefore, I subset all data from that layer for non biomass groups and
  # groups which cannot penetrate into the sediment!
  # UPDATE: Doesn't work with layers as species are not distributed through
  # the whole water column and do not appear in
  # every polygon.

  # Sum up N for invert cohorts if invert cohorts are present!
  # NOTE: invert cohorts of size 10 are not considered!
  if (select_variable == "N" & any(final_agecl != 10 & final_agecl > 1)) {
    result <- result %>%
      dplyr::group_by(polygon, layer, species, time) %>%
      dplyr::summarise(atoutput = sum(atoutput))
  }

  return(result)
}

```

Does it work? This snippet will go into `run_truth` if it does:

```{r testloadnc}

file_nc <- ("outputnordic_runresults_01ANNAGEBIO.nc")
bps <- load_bps(dir = d.name, fgs = functional.groups.file, file_init = biomass.pools.file)
allboxes <- load_box(dir = d.name, file_bgm = box.file)
boxes <- get_boundary(allboxes)


numsage <- load_nc_annage(dir = d.name,
                  file_nc = file_nc,
                  bps = bps,
                  fgs = fgs,
                  biolprm = biolprm,
                  select_groups = select_groups,
                  select_variable = "Nums",
                  check_acronyms = TRUE,
                  bboxes = boxes, 
                  verbose = TRUE)
  #if(verbose) message("Numbers read in from ANNAGEBIO.")

weightage <- load_nc_annage(dir = d.name,
                  file_nc = file_nc,
                  bps = bps,
                  fgs = fgs,
                  biolprm = biolprm,
                  select_groups = select_groups,
                  select_variable = "Weight",
                  check_acronyms = TRUE,
                  bboxes = boxes, 
                  verbose = TRUE)
  #if(verbose) message("Weight read in from ANNAGEBIO.")

```

Need to aggreagate these to see if they look reasonable; this code snippet from `om_comps` with `numsage` output subsituted for `nums`. 

```{r fullagecomps-test}
  
omlist_ss<-NOBAom_cod
source(here("/config/NOBA2Config.R"))
source(here("/config/omdimensions.R"))
source(here("/config/codsurvey.R"))
source(here("/config/codfishery.R")) 
  
n_reps <- 1

  #numbers based fishery independent survey for age and length comps
  # same user specifications as indices
  survey_N <- atlantisom::create_survey(dat = numsage,
                                        time = survtime,
                                        species = survspp,
                                        boxes = survboxes,
                                        effic = surveffic,
                                        selex = survselex)
  #Sample fish for age composition
  # if we want replicates for obs error this sample function will generate them
  age_comp_data <- list()
  for(i in 1:n_reps){
    age_comp_data[[i]] <- atlantisom::sample_fish(survey_N, surveffN)
  }

    # save age comps
  #if(save){
    saveRDS(age_comp_data, file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
  #}


```

Plot full age comps, have 20 years for cod now, no longer age class but annual ages.

```{r plotfullagecomp}

#read in comp data
age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
age_comp_data <- age_comp_data[[1]]

#fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
#fish_age_comp <- fish_age_comp[[1]]

#add this to om_indices function so that this has years when read in
#fish_age_comp$time <- as.integer(floor(fish_age_comp$time/fstepperyr))

Natage <- filter(age_comp_data, time %in% c(55:175))

Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
    geom_point() +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          Natage$species))
  
  Natageplot + facet_wrap(~time/stepperyr, ncol=6, scales="free_y")
  
```

What is the weight output? Documentation only mentions numbers output. Looks like average weight at age for each time and polygon? If so, need to do weighted average weight using the numbers at age in each polygon? 

This output is described in Atlantis code, file atlantis/atphysics/atagetracerIO.c, function `writeBMAnnAgeBioData`:

```
if (propit == out_nums_id) {  // Numbers case
                                for (b = 0; b < bm->nbox; b++) {
                                    val[b] = 0.0;
                                    for (k = 0; k < bm->wcnz; k++) {  // Sum over all individuals in the box
                                        val[b] += bm->wctr[b][k][id] * FunctGroupArray[sp].boxPopRatio[b][k][cohort][ai];
                                    }
                                }
                            } else  { // Size data - get weighted average of values per layer if numbers > min_dens
                                for (b = 0; b < bm->nbox; b++) {
                                    val[b] = 0.0;
                                    tot = 0.0;
                                    for (k = 0; k < bm->wcnz; k++) {  // Sum over all individuals in the box
                                        if (bm->wctr[b][k][id] > bm->min_dens) {
                                            val[b] += bm->wctr[b][k][id] * (bm->wctr[b][k][sn] + bm->wctr[b][k][rn]) * 5.7 * mg_2_g * FunctGroupArray[sp].boxPopRatio[b][k][cohort][ai];
                                            tot += bm->wctr[b][k][id];
                                        }
                                    }
                                    if (tot > 0.0)
                                        val[b] /= tot;
                                }
                            }
```

So weight output appears to be average weight at age in grams?

```{r wtage}

annwtage <- weightage %>%
  group_by(species, agecl, time) %>%
  summarise(wtage = mean(atoutput)) %>%
  filter(time %in% (survtime))

boxwtage <- weightage %>%
  mutate(wtage = atoutput) %>%
  filter(species %in% c("North_atl_cod"),
         polygon %in% c(1,7,24,30,38,45),
         agecl %in% c(16:20))

wtagedat <- boxwtage

wageplot <- ggplot(wtagedat, aes(time, wtage)) +
  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, " "))

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


```

So, to modify run_truth we need to add the annage output but not replace the regular output in case users want output from non-age structured groups (inverts, etc.). Also need to see if lengths can be generaged from the annage outputs, if not need to keep the full set of outputs from the .nc file for all species to generate them.

UPDATE: there is a difference between the annage weight at age and that calculated using 'calc_age2length()'. While the results from the standard nc output look reasonable for NOBA cod, the annage weight at age is half that for the standard (and half what real cod would be). More worrisome, the oldest age classes have smaller average weight at age than some younger ones. I cannot think of why this would be. Until we are sure the annage weight output is correct, we will instead interpolate the weight at age output of 'calc_age2length" for use as assessment input. 

Looks like the annagecatch output is only numbers and not weights. We don't get weight at age from fisheries usually so this is fine. 

```{r loadannagecatch}
# taken from existing load_nc

file.nc <- file.path(d.name, "outputnordic_runresults_01ANNAGECATCH.nc")

  # Load ATLANTIS output!
  at_out <- RNetCDF::open.nc(con = file.nc)

  # Get info from netcdf file! (Filestructure and all variable names)
  var_names_ncdf <- sapply(seq_len(RNetCDF::file.inq.nc(at_out)$nvars - 1),
    function(x) RNetCDF::var.inq.nc(at_out, x)$name)
  n_timesteps <- RNetCDF::dim.inq.nc(at_out, 0)$length
  n_boxes     <- RNetCDF::dim.inq.nc(at_out, 1)$length
  n_layers    <- RNetCDF::dim.inq.nc(at_out, 2)$length

RNetCDF::close.nc(at_out)
```

ANNAGECATCH.nc variable names (first hundred): `r var_names_ncdf[1:100]` 

ANNAGECATCH.nc variable names (last hundred): `r tail(var_names_ncdf, 100)`

ANNAGECATCH.nc timesteps: `r n_timesteps`

ANNAGECATCH.nc boxes: `r n_boxes`

ANNAGECATCH.nc layers: `r n_layers`

Catch file contains "IsFished" species discard at age and catch at age by fishery in numbers. So now modifying load_nc_annage to get the fishery specific outputs, which requires additional atlantis fishery setup files and additions to the config files above.

New function `load_fisheries()` needed to load the fishery specification file (a .csv), similar to `load_fgs()`:

```{r load_fisheries, eval=FALSE}

load_fisheries <- function(dir = getwd(), file_fish){
  if (is.null(dir)) {
    file.fish <- file_fish
  } else {
    file.fish <- file.path(dir, file_fish)
  }
  result <- read.table(file = file.fish, sep = ",",
    header = TRUE, stringsAsFactors = FALSE)
  return(result)
}
```

Call `load_fisheries()` within `load_nc_annage()` to rename select_variable as used in ANNAGECATCH.nc:

```{r addsnippet, eval=FALSE}

select_variable <- c("Catch", "Discard")
dir <- here("atlantisoutput","NOBA_march_2020")
file_fish <- "NoBAFisheries.csv"

#starting at line 115 in current load_nc_annage

if(select_variable %in% c("Catch", "Discard")){
    #read in fleet names
    fleetnames <- load_fisheries(dir = dir, file_fish = file_fish)
    svfish <- c()
    for(i in select_variable){
      sv <-paste0(i, "_", fleetnames$Code)
      svfish <- c(svfish, sv)
    }
    select_variable <- svfish
  }

```

Also added if statements to add fleets for fisheries outputs, if statements to separate output for annage and annagecatch.

Test new `load_nc_annage()`:

```{r testfishery}

# source(here("R/load_nc_annage.R"))
# source(here("R/load_fisheries.R"))

source(here("config/NOBA2Config.R"))
NOBAom_cod <- readRDS(file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))

fgs <- NOBAom_cod$funct.group_ss
biolprm <- NOBAom_cod$biol
#select_groups <- c('Mackerel', 'North_atl_cod', 'Capelin')
select_groups <- 'North_atl_cod'

bps <- load_bps(dir = d.name, fgs = functional.groups.file, file_init = biomass.pools.file)
allboxes <- load_box(dir = d.name, file_bgm = box.file)
boxes <- get_boundary(allboxes)


file_nc <- "outputnordic_runresults_01ANNAGECATCH.nc"
file_fish <- "NoBAFisheries.csv"

#for testing line by line
bboxes <- boxes
dir <- d.name
select_variable <- "Catch"

fishcatchage <- load_nc_annage(dir = d.name,
                  file_nc = file_nc,
                  file_fish = file_fish,
                  bps = bps,
                  fgs = fgs,
                  biolprm = biolprm,
                  select_groups = select_groups,
                  select_variable = "Catch",
                  check_acronyms = TRUE,
                  bboxes = boxes, 
                  verbose = TRUE)


fishcatchage[1:100,]

fishdiscage <- load_nc_annage(dir = d.name,
                  file_nc = file_nc,
                  file_fish = file_fish,
                  bps = bps,
                  fgs = fgs,
                  biolprm = biolprm,
                  select_groups = select_groups,
                  select_variable = "Discard",
                  check_acronyms = TRUE,
                  bboxes = boxes, 
                  verbose = TRUE)

fishdiscage
```

Do we want catch by fishery? Yes. No fishery column here, check what it is doing.
fishdiscage is empty, make sure there is no discard output from NOBA for these species (confirmed).

Add a check for multiple fisheries per species, if only one then fleet irrelevant.

NOBA has multiple fleets for cod (but they are not in output!), so add fleet column to output.

Update: there is now a fleet column in output which should index correctly if a species is caught by multiple fleets. 

Plot catch at age; aggregate over polygons, apply sample_fish:

```{r catchannage}

omlist_ss<-NOBAom_cod
source(here("/config/NOBA2Config.R"))
source(here("/config/omdimensions.R"))
source(here("/config/codsurvey.R"))
source(here("/config/codfishery.R")) 
  
n_reps <- 1

  #numbers based fishery independent survey for age and length comps
  # same user specifications as indices
  #fishery catch at age each observed timestep summed over observed polygons
  # catch at age by area and timestep
  catch_numbers <-  atlantisom::create_fishery_subset(dat = fishcatchage,
                                                      time = fishtime,
                                                      species = select_groups,
                                                      boxes = fishboxes)
  
  # if we want replicates for obs error this sample function will generate them
  catch_age_comp <- list()
  for(i in 1:n_reps){
    catch_age_comp[[i]] <- atlantisom::sample_fish(catch_numbers, fisheffN)
  }
  
  # save fishery age comps
  #if(save){
    saveRDS(catch_age_comp, file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
  #}


```

NOTE that the fleet column is lost in the `atlantisom::sample_fish` function!

Look at only every 10 years output (this is a sample, effN = 1000 for the year, which is 200 per timestep here; see calc below confirming):

```{r plotfullcatchage}

#read in comp data
# age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsFullAgeComp.rds")))
# age_comp_data <- age_comp_data[[1]]

fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsFullAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]

#add this to om_indices function so that this has years when read in
fish_age_comp$time <- fish_age_comp$time/fstepperyr

Catage <- filter(fish_age_comp, time %in% c(30.0, 40.0, 50.0, 60.0, 70.0, 79.8))

Catageplot <- ggplot(Catage, aes(x=agecl, y=atoutput)) +
    geom_point() +
    theme_tufte() +
    labs(subtitle = paste(scenario.name,
                          Catage$species))
  
  Catageplot + facet_wrap(~time, ncol=3, scales="free_y")

```

```{r testeffn}

fish_age_comp %>%
  filter(time<31) %>%
  summarise(sum(atoutput))

```

Appears to be working for single fleet per species now.

Quick test of `run_truth()` with annage option turned on for NOBA:

```{r testruntruth, eval=FALSE}

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

  #Load functional groups
  funct.groups <- atlantisom::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

  # # load biolprm file
  #biol <- atlantisom::load_biolprm(d.name, biol.prm.file)
  # 
  # # load run_prm
  # runpar <- atlantisom::load_runprm(d.name, run.prm.file)
  # 
  # # load box info file
  # boxpars <- atlantisom::load_box(d.name, box.file)


testtruth <- atlantisom::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,
                       file_fish = fisheries.file,
                       verbose = TRUE,
                       annage = TRUE,
                       save = TRUE
    )

```

After much debugging to properly index multiple fleets per species in catch at annual age output, this produced a run_truth.Rdata file and took about 40 minutes to run.

Now all of the wrappers need to read from annage to produce composition data. `om_init()`  now takes annage flag from the config file. `om_species()` looks for annage bio and catch output in the truth file and includes it in the species file (adding NULL if it isn't there). `om_comps()` looks for annage bio and catch in the species file and adds them to the comps list if there (adding NULL if they aren't there). 

Also still need the weight at age interpolation function. Incorporate this at the `om_comps()` step. Next rmd!

Running `om_init()`, `om_species()`, `om_index()`, `om_comps()` in chunk above to make sure the original output is preserved, then plot the annage output [above](#annageoutput). 


## Write input data for specific assessment models

We will need to alter this for SAM, SS3, etc. Right now just makes SS3.

**THIS IS STILL THE CALIFORNIA CURRENT INFO, TO BE UPDATED**

```{r writedat, eval=FALSE}

require(r4ss)
#omlist_ss<- ?
source(here("config/omdimensions.R"))
source(here("/config/codsurvey.R"))
source(here("/config/codfishery.R")) 

#Directory with SS files
model_dir <- "NOBA_cod_files/"

#Name of SS data file
datfile_name <- "ss3.dat"

#CVs for length at age, catch, and survey
CVs <- list("lenage"=0.1, "fishery"=0.01, "survey"=0.1)

#Load in model object from SAM
noba_cod_out <- system.file("inst","extdata","NOBA_cod_files","NEAcod-2020","run","model.RData",
                            package = "atlantisom") 
#This is the SAM object
noba_fit <- noba_cod_out$fit

properties <- c("minAge","maxAge", "maxAgePlusGroup")

age_info <- lapply(properties, get, noba_fit$conf)
names(age_info) <- properties

#install SAM
remotes::install_github("fishfollower/SAM/stockassessment")
require(stockassessment)




#Month of survey/fishing
#survey_month <- 7
#fishing_month <- 1

#Years to survey, assuming survey is once per year
#survey_years <- survey_sample_full[(burnin+1):(burnin+nyears)] # my original
survey_years <- survey_sample_full[burnin:(burnin+nyears-1)] #from Christine's new sardine_config.R


#read time series data
survObsBiom <- readRDS(file.path(d.name, paste0(scenario.name, "surveyB.rds")))
truecatchbio_ss <- readRDS(file.path(d.name, paste0(scenario.name, "fishCatch.rds")))

survObsBiom <- survObsBiom[[1]]
truecatchbio_ss <- truecatchbio_ss[[1]]

#add this to om_indices function so that this has years when read in
truecatchbio_ss$time <- as.integer(truecatchbio_ss$time/365)


#load dummy dat file
stocksynthesis.data <- r4ss::SS_readdat_3.30(paste0("./inst/extdata/",
model_dir,
datfile_name))

#Test writing CPUE in biomass
  stocksynthesis.data <- atlantisom::SS_write_ts(ss_data_list = stocksynthesis.data,
                ts_data = list(survObsBiom$atoutput[fish_years],
                  truecatchbio_ss$atoutput[truecatchbio_ss$time %in% fish_years]),
                CVs = c(CVs$survey,
                        CVs$fishery),
                data_years = list((survObsBiom$time[fish_years]-survey_sample_time)/timestep+1,             fish_years),
            sampling_month = list(rep(survey_month,nyears),
                                rep(fishing_month,nyears)),
                units = c("biomass","biomass"),
                data_type=c("CPUE","catch"),
                fleets = c(2,1))

#test  
stocksynthesis.data$CPUE

stocksynthesis.data$catch

#read in comp data
age_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsAgeComp.rds")))
age_comp_data <- age_comp_data[[1]]

fish_age_comp <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsAgeComp.rds")))
fish_age_comp <- fish_age_comp[[1]]


#Get the age bins
age_bin_names <- names(stocksynthesis.data$agecomp)[10:length(names(stocksynthesis.data$agecomp))]
age_bins <- sub("a","",age_bin_names)

require(maditr)
# add dependency on maditr::dcast to atlantisom
age_comp_flat <- atlantisom::reformat_compositions(age_comp_data,
                     round.places = 4,
                     comp_type = "agecomp")

## Write age composition data for survey
stocksynthesis.data <- atlantisom::SS_write_comps(ss_data_list = stocksynthesis.data,
               comp_matrix = list(age_comp_flat[burnin:(burnin+nyears-1),]),
               data_rows = list(stocksynthesis.data$styr:(stocksynthesis.data$styr+nyears-1)),
               sampling_month = list(rep(survey_month,nyears)),
               data_type = c("agecomp"),
               fleet_number = c(1),
               bins = list(age_bins),
               caal_bool = c(FALSE))
stocksynthesis.data$agecomp

#length comps
len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "survObsLenComp.rds")))
fish_len_comp_data <- readRDS(file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))

len_comp_data <- len_comp_data[[1]]
fish_len_comp_data <- fish_len_comp_data[[1]]

#add this to om_indices function so that this has years when read in
fish_len_comp_data$time <- as.integer(floor(fish_len_comp_data$time/fstepperyr))
fish_age_comp$time <- as.integer(floor(fish_age_comp$time/fstepperyr))


if(fstepperyr>1){
  fish_len_comp_anndata <- fish_len_comp_data %>%
    #mutate(yr = floor(time/fstepperyr)) %>%
    #group_by(species, agecl, lower.bins, upper.bins, time=as.integer(yr)) %>%
    group_by(species, agecl, lower.bins, upper.bins, time) %>% 
    summarise(annnatlength=sum(atoutput)) %>%
    rename(atoutput = annnatlength)
} else {
  fish_len_comp_anndata <- fish_len_comp_data
}

caal_comp_flat <- atlantisom::reformat_compositions(len_comp_data,                                round.places=4,
                comp_type="caalcomp")


#remove burnin
caal_comp_final <- filter(caal_comp_flat,
                         time %in% survey_years)


#Add over age classes to get sample size
len_comp_flat <- atlantisom::reformat_compositions(len_comp_data,
                                  round.places = 0,
                           comp_type="lencomp")
#remove burnin
len_comp_final <- filter(len_comp_flat,
                         time %in% survey_years)

length_bins <- as.integer(names(len_comp_final))
length_bins <- length_bins[!is.na(length_bins)]

# fishery length comps are still 5 timesteps per year
# need to aggregate to annual (done above)
# also,  make effN annual goal/fstepsperyr (done above)
fish_len_comp_flat <- atlantisom::reformat_compositions(fish_len_comp_anndata,
                                  round.places = 0,
                           comp_type="lencomp")

#remove burnin works after adjustment above
fish_len_comp_final <- filter(fish_len_comp_flat,
                         time %in% fish_years)

notbins <- c("time", "nsamp")

# fish_length_bins <- as.integer(names(fish_len_comp_final))
# fish_length_bins <- fish_length_bins[!is.na(fish_length_bins)]

# need to fill empty length bins with 0s to have same bins as survey for SS_write_comps
missing.lengths <- setdiff(length_bins, names(fish_len_comp_final)[!names(fish_len_comp_final) %in% notbins])
fish_len_comp_final[as.character(missing.lengths)] <- 0                    # Add them, filled with '0's
fish_len_comp_final <- fish_len_comp_final[c("time", length_bins, "nsamp")] #


# fishery age comps also 5 timesteps per year
if(fstepperyr>1){
  fish_age_comp_anndata <- fish_age_comp %>%
    #mutate(yr = floor(time/fstepperyr)) %>%
    #group_by(species, agecl, time=as.integer(yr)) %>%
    group_by(species, agecl, time) %>%
    summarise(annnatage=sum(atoutput)) %>%
    rename(atoutput = annnatage)
} else {
  fish_age_comp_anndata <- fish_age_comp
}

fish_age_comp_flat <- atlantisom::reformat_compositions(fish_age_comp_anndata,
                           comp_type="agecomp")

#remove burnin (not necessary?fish comps made with fish_years only) 
fish_age_comp_final <- filter(fish_age_comp_flat,
                         time %in% fish_years)

# #SS_write_comps breaking because fishery age bins start with 2 not 1; extracting bins from fish file may help?
# fish_age_bins <- names(fish_age_comp_flat)[!names(fish_age_comp_flat) %in% notbins]

# that leaves an empty column in data file, so instead fill with 0s
missing.ages <- setdiff(age_bins, names(fish_age_comp_final)[!names(fish_age_comp_final) %in% notbins])
fish_age_comp_final[missing.ages] <- 0                    # Add them, filled with '0's
fish_age_comp_final <- fish_age_comp_final[c("time", age_bins, "nsamp")]

comp_list <- list(caal_comp_final,len_comp_final, fish_age_comp_final, fish_len_comp_final)

apply_month <- list(rep(survey_month, nrow(comp_list[[1]])), 
                    rep(survey_month, nrow(comp_list[[2]])),
                    rep(fishing_month,nrow(comp_list[[3]])),
                    rep(fishing_month,nrow(comp_list[[4]])))


# This now runs by ensuring that survey and fishery compositions have the same bins 
# (filled with 0s for missing bins in fishery relative to survey)

# Write CAAL and length composition data
stocksynthesis.data <- atlantisom::SS_write_comps(ss_data_list = stocksynthesis.data, 
                                      comp_matrix = comp_list, 
                                      data_rows = list((comp_list[[1]]$time-survey_sample_time)/timestep + 1 , (survey_years-survey_sample_time)/timestep + 1,fish_years,fish_years),
                                      sampling_month = apply_month, 
                                      data_type = rep(c("agecomp", "lencomp"),2), 
                                      fleet_number = c(2,2,1,1),  
                                      bins = list(age_bins, 
                                                  length_bins, 
                                                  age_bins, 
                                                  length_bins), 
                                      caal_bool = c(TRUE, rep(FALSE,3)))

head(stocksynthesis.data$lencomp)
head(stocksynthesis.data$agecomp)

#Change length bin structure to match atlantis data
stocksynthesis.data$lbin_vector <- length_bins

#Get correct number of length bins
stocksynthesis.data$N_lbins <- length(length_bins)

#Set lbin_method to 1 - this makes the population length bins match the data bins 
#When lbin_method==1, we just comment out the binwidth, minimum, and maximum size arguments since they aren't used
stocksynthesis.data$lbin_method <- 1
stocksynthesis.data$binwidth <- "#"
stocksynthesis.data$minimum_size <- "#"
stocksynthesis.data$maximum_size <- "#"

#Change minimum sample size to 0.001 for CAAL data (SS won't let it go lower than this)
stocksynthesis.data$age_info$minsamplesize <- rep(0.001,2)

SS_writedat_3.30(stocksynthesis.data, outfile = paste0("./inst/extdata/",model_dir,
datfile_name),
                 overwrite=TRUE)

```