WrapperDevelopment.Rmd

---
title: "Develop atlantisom wrappers"
author: "Sarah Gaichas and Christine Stawitz"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

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

## Too many steps! and too much memory!

We can run atlantisom with many functions to get stock assessment model inputs.

To make life easier, here are some wrapper functions to do the basic things we want.

Basic desired functionality:

1. initialize: locate model output (function to make a config file?), read in basic run parameters and species names, run_truth for all once and save output 
    + User to make up front decision--want all species or just a subset?
    + May want multiple species that need to be split in a later step for models
    + wraps functions 
        + sourcing model config file--add function to make config file?
      + specifies top directory for source files, names key files for use later
      + load_fgs
      + load all other needed files--YOY, catch
      + run_truth
1. get assessment input data: configure survey, run all survey functions, configure fishery, run all fishery functions (how to save output?)
    + Based on decision above, save interim steps for use later
    + Need a switch for single vs multispecies estimation models?
    + Separate index and comp functions in output
    + indices:
      + create survey config file(s), default is a census of all see census_spec.R
        + survey index in biomass or numbers? (uses truth$biomass_ages vs truth$nums)
        + wraps functions create_survey, sample_survey_biomass
      + create fishery config file(s)
    + comps:
      + select age only or age and length
        + age only wraps create_survey using nums with sample_fish to get age comps
        + age and length (with weight at age)
          + create_survey followed by sample_fish, 
          + aggregateDensityData for resn and structn followed by sample_fish
          + inputs all to calc_age2length
      + catch comps wrap create_fishery_subset with sample_fish to get age comps, same as above for lengths
1. write input data for specific model: wrap functions in development for SS3, add functions for other models (ASAP, multispecies models)
    + wrap r4ss functions to read and write dat and ctl files
    + possibly start with set of dummy SS files with package (take from ss3sim?)
      + read in dat file
      + indices wrap SS_write_ts
      + comps wrap reformat_compositions, SS_write_comps, some of the checks for bin matches
      + write full dat file
    + start with dummy ctl file (any from ss3sim?)
      + read in ctl file
      + wrap SS_write_bio using outputs of calc_Z, load_YOY, load_biolprm, load_runprm, etc
1. run specified model--automate for scenarios, save outputs 
1. skill assessment metrics: compare stored Atlantis truth with model output 
    
### atlantisom initialize function

Things to add: checks for output timestep?

```{r ominit, eval=FALSE}

library(tidyverse)

om_init <- function(config = configfile){
  
  # Where are the atlantis output files? Consider filling with shiny app in future
  source(config)
  # needs these files, for example config file CC3config.R is:
  
  # d.name <- here::here("atlantisoutput","CC_2063_OA_OFF_22")
  # functional.groups.file <- "CalCurrentV3Groups.csv"
  # biomass.pools.file <- "DIVCalCurrentV3_Biol.nc"
  # biol.prm.file <- "CalCurrentV3_Biol.prm"
  # box.file <- "CalCurrentV3_utm.bgm"
  # initial.conditions.file <- "DIVCalCurrentV3_Biol.nc"
  # run.prm.file <- "CalCurrentV3_run.xml"
  # scenario.name <- "CCV3"
  # bioind.file <- "outputCCV3BiomIndx.txt"
  # catch.file <- "outputCCV3Catch.txt"

  #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 true total biomass in tons
  truetotbio <- atlantisom::load_bioind(d.name, file_bioind = bioind.file, fgs = funct.groups)
  
  # load true catch in tons
  truecatchbio <- atlantisom::load_catch(d.name, file_catch = catch.file, fgs = funct.groups)

  # load YOY
  YOY <- atlantisom::load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))
  
  # load biol_prm
  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)

  # 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 <- 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,
                       verbose = TRUE
    )
  } else {
    truth <- get(load(file.path(d.name,
                                paste0("output", scenario.name, "run_truth.RData"))))
  }
  
  omlist <-list("funct.groups" = funct.groups, 
                "funct.group.names" = funct.group.names,
                "truetotbio" = truetotbio,
                "truecatchbio" = truecatchbio,
                "YOY" = YOY,
                "biol" = biol,
                "runpar" = runpar,
                "boxpars" = boxpars,
                "truth" = truth)
  
  return(omlist)
}
  

```

Usage: 

```{r exinit, eval=FALSE}

library(here)
library(tidyverse)
library(atlantisom)
CC3om <- om_init(here("config/CC3config.r"))

```

### Get assessment input data

Split to single species or subset of assessed species, get index data, get compositional data. By default, remove full om results and keep only species subsets to save memory.

#### Split to focal species

Is this where we should reconcile timesteps? depends what functions are expecting.

```{r inputs, eval=FALSE}

om_species <- function(species = spp, omlist, removefullom = TRUE){
  # spp format c("speciesname1", "speciesname2")
  if(!all(species %in% omlist$funct.group.names)) stop("species name not found") 
  species_ss <- species
  
  #subset species true bio
  truetotbio_ss <- omlist$truetotbio[omlist$truetotbio$species %in% species_ss,]

  #subset species true catch
  truecatchbio_ss <- omlist$truecatchbio[omlist$truecatchbio$species %in% species_ss,]

  #subset species YOY
  # get code matching species name to split YOY file
  code_ss <- omlist$funct.groups$Code[which(omlist$funct.groups$Name %in% species_ss)]
  # cut to a single species in YOY file
  YOY_ss <- omlist$YOY %>%
    select(Time, paste0(code_ss, ".0"))
  # reformat to be like all the other objects
  
  # numbers at agecl at full resolution (all polygons and layers)
  truenums_ss <- omlist$truth$nums[omlist$truth$nums$species %in% species_ss,]
  
  # biomass at agecl at full resolution (all polygons and layers)
  truebio_ss <- omlist$truth$biomass_ages[omlist$truth$biomass_ages$species %in% species_ss,]
  
  # reserve nitrogen at agecl at full resolution
  trueresn_ss <- omlist$truth$resn[omlist$truth$resn$species %in% species_ss,]
  
  # structural nitrogen at agecl at full resolution
  truestructn_ss <- omlist$truth$structn[omlist$truth$structn$species %in% species_ss,]
  
  # catch in numbers at agecl at full resoluation (all polygons, no layer in output)
  truecatchnum_ss <- omlist$truth$catch[omlist$truth$catch$species %in% species_ss,]
  
  #subset species biol parameters? no, may miss some globals that aren't by species
  
  #subset species functional group info
  funct.groups_ss <- omlist$funct.groups[omlist$funct.groups$Code %in% code_ss,]
  
  #keep the runpar and also need boxes for survey selection
  
  omlist_ss <- list("species_ss" = species_ss,
                    "code_ss" = code_ss,
                    "truetotbio_ss" = truetotbio_ss,
                    "truecatchbio_ss" = truecatchbio_ss,
                    "YOY_ss" = YOY_ss,
                    "truenums_ss" = truenums_ss,
                    "truebio_ss" = truebio_ss,
                    "trueresn_ss" = trueresn_ss,
                    "truestructn_ss" = truestructn_ss,
                    "truecatchnum_ss" = truecatchnum_ss,
                    "funct.group_ss" = funct.groups_ss,
                    "biol" = omlist$biol,
                    "boxpars" = omlist$boxpars,
                    "runpar" = omlist$runpar)
  
  if(removefullom) rm(omlist) #not removing passed data object
  
  return(omlist_ss)
}

```

Usage: 

```{r exsp, eval=FALSE}

CC3om_sardine <- om_species(c("Pacific_sardine"), CC3om)

```

#### Generate index data for assessment

The hardcoded script omdimensions.R needs to sit in the atlantisom config folder.

```{r omindices, eval=FALSE}

om_index <- function(usersurvey = usersurvey_file, 
                     userfishery = userfishery_file, 
                     omlist_ss, 
                     n_reps = n_reps,
                     save = TRUE){
  
  #one script for dimension parameters to be used in multiple functions
  source("config/omdimensions.R", local = TRUE)
  
  # user options for survey--default is a census with mid-year sample
  source(usersurvey, local = TRUE)


  #biomass based fishery independent survey index
  # this uses result$biomass_ages to sample biomass directly, no need for wt@age est
  survey_B <- atlantisom::create_survey(dat = omlist_ss$truebio_ss,
                                        time = survtime,
                                        species = survspp,
                                        boxes = survboxes,
                                        effic = surveffic,
                                        selex = survselex)
  
  # call sample_survey_biomass with a bunch of 1000s for weight at age
  # in the code it multiplies atoutput by wtatage/1000 so this allows us to use
  # biomass directly
  wtage <- data.frame(species=rep(survspp, each=max(age_classes)),
                      agecl=rep(c(1:max(age_classes)),length(survspp)),
                      wtAtAge=rep(1000.0,length(survspp)*max(age_classes)))
  
  # this is the step to repeat n_reps time if we want different realizations 
  # of the same survey design specified above; only observation error differs
  # using the census cv of 0 will produce identical reps!
  survObsBiomB <- list()
  for(i in 1:n_reps){
    survObsBiomB[[i]] <- atlantisom::sample_survey_biomass(survey_B, surv_cv, wtage)
  }
  
  #save survey indices, takes a long time to generate with lots of reps/species
  if(save){
    saveRDS(survObsBiomB, file.path(d.name, paste0(scenario.name, "surveyB.rds")))
  }
  
  #configure the fishery, a default is in config/fisherycensus.R
  #fishery configuration can specify only area and time of observation
  #fishery species inherited from omlist_ss
  source(userfishery, local = TRUE)
  
  #we are not currently subsetting fishery catch because we cannot correct catch.nc
  #  instead the catch in biomass from catch.txt is read in for the index
  #  we do not apply any cv to this, but we could this way (default cv=0)
  
  fishObsCatchB <- list()
  for(i in 1:n_reps){
    fishObsCatchB[[i]] <- atlantisom::sample_fishery_totcatch(omlist_ss$truecatchbio_ss, fish_cv)
  }

  if(save){
    saveRDS(fishObsCatchB, file.path(d.name, paste0(scenario.name, "fishCatch.rds")))
  }
 
  indices <- list("survObsBiomB" = survObsBiomB,
                  "fishObsCatchB" = fishObsCatchB)
  
  return(indices)
}

```

Usage: 

```{r exindex, eval=FALSE}

CC3om_sard_ind <- om_index(usersurvey = here("config/usersurvey.R"), 
                           userfishery = here("config/fisherycensus.R"),
                           omlist_ss = CC3om_sardine, 
                           n_reps = 5, 
                           save = TRUE)

```

#### Generate compositional data for assessment

Note that if a census is truly desired for age comps, don't use the sample fish function but just aggregate. I'll make a function for that too.

ToDo: write the atlantisom function sample_ages to subset the age comp here, that way there will be way more lengths than ages, similar to most real world sampling.

```{r omcomps, eval=FALSE}

om_comps <- function(usersurvey = usersurvey_file, 
                     userfishery = userfishery_file, 
                     omlist_ss, 
                     n_reps = n_reps,
                     save = TRUE){
  
  #one script for dimension parameters to be used in multiple functions
  source("config/omdimensions.R", local = TRUE)
  
  # user options for survey--default is a census with mid-year sample
  source(usersurvey, local = TRUE)
  
  #numbers based fishery independent survey for age and length comps
  # same user specifications as indices
  survey_N <- atlantisom::create_survey(dat = omlist_ss$truenums_ss,
                                        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, "survObsAgeComp.rds")))
  }
  
  #weights needed for weight at age and length comp calcs
  # aggregate true resn per survey design
  survey_aggresn <- atlantisom::aggregateDensityData(dat = omlist_ss$trueresn_ss,
                                                     time = survtime,
                                                     species = survspp,
                                                     boxes = survboxes)
  
  # aggregate true structn per survey design
  survey_aggstructn <- atlantisom::aggregateDensityData(dat = omlist_ss$truestructn_ss,
                                                        time = survtime,
                                                        species = survspp,
                                                        boxes = survboxes)
  
  #dont sample these, just aggregate them using median
  structnss <- atlantisom::sample_fish(survey_aggresn, surveffN, sample = FALSE)
  
  resnss <- atlantisom::sample_fish(survey_aggstructn, surveffN, sample = FALSE)
  
  #this is all input into the length function, replicates follow age comp reps
  #  separating the length comps from the weight at age here
  survey_lenwt <- list()
  survObsLenComp <- list()
  survObsWtAtAge <- list()
  
  for(i in 1:n_reps){
    survey_lenwt[[i]] <- atlantisom::calc_age2length(structn = structnss,
                                                     resn = resnss,
                                                     nums = age_comp_data[[i]],
                                                     biolprm = omlist_ss$biol, 
                                                     fgs = omlist_ss$funct.group_ss,
                                                     maxbin = maxbin,
                                                     CVlenage = lenage_cv, 
                                                     remove.zeroes=TRUE)
    
    survObsLenComp[[i]] <- survey_lenwt[[i]]$natlength
    survObsWtAtAge[[i]] <- survey_lenwt[[i]]$muweight
  }
  
  if(save){
    saveRDS(survObsLenComp, file.path(d.name, paste0(scenario.name, "survObsLenComp.rds")))
    saveRDS(survObsWtAtAge, file.path(d.name, paste0(scenario.name, "survObsWtAtAge.rds")))
  }
    
  #now do fishery comps
  # user options for fishery--default is a census with mid-year sample
  source(userfishery, local = TRUE)
  
  #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 = omlist_ss$truecatchnum_ss,
                                                      time = fishtime,
                                                      species = survspp,
                                                      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, "fishObsAgeComp.rds")))
  }

  #Get catch weights for length comp calc
  # aggregate true resn per fishery subset design
  catch_aggresnss <- atlantisom::aggregateDensityData(dat = omlist_ss$trueresn_ss,
                                                      time = fishtime,
                                                      species = survspp,
                                                      boxes = fishboxes)
  
  # aggregate true structn fishery subsetdesign
  catch_aggstructnss <- atlantisom::aggregateDensityData(dat = omlist_ss$truestructn_ss,
                                                         time = fishtime,
                                                         species = survspp,
                                                         boxes = fishboxes)
  
  #dont sample these, just aggregate them using median
  catch_structnss <- atlantisom::sample_fish(catch_aggstructnss, fisheffN, sample = FALSE)
  
  catch_resnss <- atlantisom::sample_fish(catch_aggresnss, fisheffN, sample = FALSE)
  
  # these fishery lengths and weight at age are each output timestep
  #same structure as above for surveys, replicates follow age comp reps
  #  separating the length comps from the weight at age here
  fishery_lenwt <- list()
  fishObsLenComp <- list()
  fishObsWtAtAge <- list()
  
  for(i in 1:n_reps){
    fishery_lenwt[[i]] <- atlantisom::calc_age2length(structn = catch_structnss,
                                                      resn = catch_resnss,
                                                      nums = catch_age_comp[[i]],
                                                      biolprm = omlist_ss$biol, 
                                                      fgs = omlist_ss$funct.group_ss,
                                                      maxbin = maxbin,
                                                      CVlenage = lenage_cv, 
                                                      remove.zeroes=TRUE)
    
    fishObsLenComp[[i]] <- fishery_lenwt[[i]]$natlength
    fishObsWtAtAge[[i]] <- fishery_lenwt[[i]]$muweight
  }
  
  if(save){
    saveRDS(fishObsLenComp, file.path(d.name, paste0(scenario.name, "fishObsLenComp.rds")))
    saveRDS(fishObsWtAtAge, file.path(d.name, paste0(scenario.name, "fishObsWtAtAge.rds")))
  }
    
  
  comps <- list("survObsAgeComp" = age_comp_data,
                "survObsLenComp" = survObsLenComp,
                "survObsWtAtAge" = survObsWtAtAge,
                "fishObsAgeComp" = catch_age_comp,
                "fishObsLenComp" = fishObsLenComp,
                "fishObsWtAtAge" = fishObsWtAtAge)
  
  return(comps)
}

```

Usage:
```{r compex, eval=FALSE}

CC3om_sard_comp <- om_comps(usersurvey = here("config/usersurvey.R"), 
                           userfishery = here("config/fisherycensus.R"),
                           omlist_ss = CC3om_sardine, 
                           n_reps = 1, 
                           save = TRUE)

```

Test these with updated length-weight parameters for CCA sardines:

I updated CalCurrentV3_Biol.prm with new sardine length-weight parameters from Isaac

>On the basis of the last sardine stock assessment that seems to have estimated growth, 
http://www.pcouncil.org/wp-content/uploads/2014_CPS_SAFE_Sardine_assessment_Appendix_C.pdf 
>
>I say length-weight A and B params,  for W = A*(L^B)
>
>should be A= 0.0075242
>and B = 3.2332
>
>where weight is in grams
>and length is in CM. 

and renamed it CalCurrentV3_Biol_sardLWcorrection.prm

Config files for this extraction should match previous sardine survey and fishery:

CC3Config_sardLWcorr.R
sardinesurvey.R
sardinefishery.R

```{r full-test}

library(here)
library(tidyverse)
library(atlantisom)

CC3om <- om_init(here("config/CC3config_sardLWcorr.R"))

CC3om_sardine <- om_species(c("Pacific_sardine"), CC3om)

CC3om_sard_ind <- om_index(usersurvey = here("config/sardinesurvey.R"), 
                           userfishery = here("config/sardinefishery.R"),
                           omlist_ss = CC3om_sardine, 
                           n_reps = 1, 
                           save = TRUE)

CC3om_sard_comp <- om_comps(usersurvey = here("config/sardinesurvey.R"), 
                           userfishery = here("config/sardinefishery.R"),
                           omlist_ss = CC3om_sardine, 
                           n_reps = 1, 
                           save = TRUE)



```

### Write input data for specific assessment models

Use the saved outputs of above to write a Stock Synthesis dat file
This is not a function yet, copy in Christine's

**See corrections below, in fishery time fixed and nsamp added to length comps**

```{r writedat}

require(r4ss)
omlist_ss<-CC3om_sardine
source(here("/config/CC3Config_sardLWcorr.R"))
source(here("/config/omdimensions.R"))
source(here("/config/sardinesurvey.R"))
source(here("/config/sardinefishery.R")) 

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

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

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

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

```