SardinesOnlyatlantisom2SStest.Rmd

---
title: "Sardines only end-to-end atlantisom test"
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)
```

```{r message=FALSE, warning=FALSE}
library(tidyr)
require(dplyr)
library(ggplot2)
library(data.table)
library(here)
library(ggforce)
library(ggthemes)
library(atlantisom)
```

## Introduction

Now we test the `atlantisom` to SS workflow with a new California Current run that includes climate forcing, recruitment variability, and the fishing scenario outlined for our project in Norway: "output_CC_2063_run11".

We will go from atlantis model output to SS input in this page for a single species, "sardine", that doesn't require any age class splitting.

As of 31 May 2019, the `run_truth` function has been modified to do the catch numbers fix for legacy models.

```{r initialize}

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/NOBAConfig.R"))

```

```{r get_names, message=FALSE, warning=FALSE}
#Load functional groups
funct.groups <- load_fgs(dir=d.name,
                         file_fgs = functional.groups.file)
#Get just the names of active functional groups
funct.group.names <- funct.groups %>% 
  filter(IsTurnedOn == 1) %>%
  select(Name) %>%
  .$Name

```

This can be run with `species_ss` in select_groups to get results only for that species. Here I'm running with all species; took about 50 minutes to return.

```{r get_truth_ss}

# 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,
                     verbose = TRUE
  )
} else {
  truth <- get(load(file.path(d.name,
                              paste0("output", scenario.name, "run_truth.RData"))))
}
```

In the big picture, we will have the true population dynamics from Atlantis which we need for reference to compare with estimated population dynamics from the stock assessment estimation models. 

Ideally we will have a function that gets the whole truth (and nothing but the truth) from Atlantis and stores it in a sensible way for later comparison and calculation of performance metrics.

Then we will have a function that implements the "survey" that produces the stock assessment inputs and sends those to the SS file writing functions.

## Truth from OM

The following settings should achieve a survey that samples all Atlantis model output timesteps, all species, and all model polygons, with perfect efficiency and full selectivity for all ages: 

```{r census-spec, message=FALSE, warning=FALSE}

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

```

Get true total biomass, total numbers. This code gets results for all species and saves the output true B and N.

```{r trueB-N, eval=FALSE}

# this uses result$biomass_ages to sample biomass directly, no need for wt@age est
survey_testBall <- create_survey(dat = truth$biomass_ages,
                                 time = timeall,
                                 species = survspp,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# make up a constant 0 cv for testing
surv_cv <- data.frame(species=survspp, cv=rep(0.0,length(survspp)))

# 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=10),
                    agecl=rep(c(1:10),length(survspp)),
                    wtAtAge=rep(1000.0,length(survspp)*10))

surveyB_frombio <- sample_survey_biomass(survey_testBall, surv_cv, wtage)

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

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

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


```

Plot true biomass and numbers for "sardines." 

```{r plot-true-B-N}

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

trueB_ss <- trueB %>%
  filter(species==species_ss)

# read Atlantis output files
atBtxt2 <- read.table(file.path(d.name, paste0("output", scenario.name, "BiomIndx.txt")), header=T)
groupslookup <- funct.groups %>%
  filter(IsTurnedOn > 0)

atBtxt2tidy <- atBtxt2 %>%
  select(Time:DIN) %>%
  #select(Time, FPL:DIN) %>%
  rename_(.dots=with(groupslookup, setNames(as.list(as.character(Code)), Name))) %>%
  gather(species, biomass, -Time) %>%
  filter(species %in% c("Pacific_sardine"))


plotB <-ggplot() +
  geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="survey census B"), 
            alpha = 10/10) +
  geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
             alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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


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

trueN_ss <- trueN %>%
  filter(species==species_ss)

plotN <-ggplot() +
  geom_line(data=trueN_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") 

```

Get true population age comp, length comp, and weight at age. Now we go to single species to speed things up.

```{r true-comps, warning=FALSE, message=FALSE, eval=FALSE}

# wrap these individual calls into a function
survey_N_ss <- create_survey(dat = truth$nums,
                                 time = timeall,
                                 species = species_ss,
                                 boxes = boxall,
                                 effic = effic1,
                                 selex = selex1)

# apply default sample fish to get numbers, effNhigh samples all, see census_spec
# NOPE! effNhigh is not big enough for census of sardines, bigger breaks function
#numssshigh <- sample_fish(survey_N_ss, effNhigh)

numsss <- aggregate(survey_N_ss$atoutput,list(survey_N_ss$species,survey_N_ss$agecl,survey_N_ss$time),sum)
names(numsss) <- c("species","agecl","time","atoutput")
numssshigh <- data.frame(species = numsss$species,
					  agecl = numsss$agecl,
					  polygon = NA,
					  layer = NA,
					  time = numsss$time,
					  atoutput = numsss$atoutput)

# save age comps
saveRDS(numssshigh, file.path(d.name, paste0(scenario.name, "natage_census_sard.rds")))

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

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

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

resnss <-  sample_fish(aggresnss, effNhigh, sample = FALSE)

# this function does the length sampling and also weight at age
length_census_ss <- calc_age2length(structn = structnss,
                                 resn = resnss,
                                 nums = numssshigh,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

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

```

Plot some of the age, length comps, and wtatage

```{r plot-true-agecomps}

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

Natageplot <- ggplot(Natage, aes(x=agecl, y=atoutput)) +
  geom_point() +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        Natage$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")

```

```{r lf-plot}

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

len <- lengthwt_census_ss$natlength

lfplot <- ggplot(len, aes(upper.bins)) +
  geom_bar(aes(weight = atoutput)) +
  theme_tufte() +
  labs(subtitle = paste(scenario.name,
                        len$species))

lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 1, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 2, scales="free_y")
lfplot + facet_wrap_paginate(~time, ncol=4, nrow = 4, page = 3, scales="free_y")


```

```{r wtage-plot}

wtage_ann <- lengthwt_census_ss$muweight %>%
  filter(time %in% annualmidyear)

# reverse to show agecl time series of wt
wageplot <- ggplot(wtage_ann, 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")

```

Get true population parameters (M in progress, S-R, growth)

```{r true-poppars}

```

Get true catch in tons (may not be different from input for our project)

```{r true-catchtot}

truecatchbio <- load_catch(d.name, file_catch = catch.file, fgs = funct.groups)

truecatchbio_ss <- truecatchbio[truecatchbio$species == species_ss,]

# note: time output of this file is in days
# what model timestep is this output matching?
# rescale catch in biomass time from days to years
truecatchbio_ss$time <- as.integer(truecatchbio_ss$time/365)

plotcatch <- ggplot() +
  geom_line(data=truecatchbio_ss, aes(x=time,y=atoutput, color="true catch bio"),
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

```

Get true catch at age and other comps (may not be different from input for our project). First see if catch output is at same timescale as population by checking toutfinc: `r fstepperyr`


```{r true-catage-lf-wtage, eval=FALSE}

# get survey nums with full (no) selectivity
catchN_ss <- create_fishery_subset(dat = truth$catch,
                                 time = timeall,
                                 species = species_ss,
                                 boxes = boxall)

# try summing this to get true catch in nums by timestep with surv_cv = 0
totcatchN <- sample_survey_numbers(catchN_ss, surv_cv)

# save true total catch in numbers by timestep, need to sum to year
saveRDS(totcatchN, file.path(d.name, paste0(scenario.name, "catchN_census_sard.rds")))

# apply default sample fish as before to get numbers
# no, dont apply sample, effN isnt high enough and can't be made high enough
#catch_numssshigh <- sample_fish(catchN_ss, effNhigh)

# just aggregte true numbers over layer and polygon
catch_numsss <- aggregate(catchN_ss$atoutput,list(catchN_ss$species,catchN_ss$agecl,catchN_ss$time),sum)
names(catch_numsss) <- c("species","agecl","time","atoutput")
catch_numssshigh <- data.frame(species = catch_numsss$species,
					  agecl = catch_numsss$agecl,
					  polygon = NA,
					  layer = NA,
					  time = catch_numsss$time,
					  atoutput = catch_numsss$atoutput)


# save true catch at age
saveRDS(catch_numssshigh, file.path(d.name, paste0(scenario.name, "catchage_census_sard.rds")))


# aggregate true resn per survey or fishery subset design
catch_aggresnss <- aggregateDensityData(dat = truth$resn,
                                 time = timeall,
                                 species = species_ss,
                                 boxes = boxall)

# aggregate true structn per survey or fishery subsetdesign
catch_aggstructnss <- aggregateDensityData(dat = truth$structn,
                                 time = timeall,
                                 species = species_ss,
                                 boxes = boxall)

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

catch_resnss <-  sample_fish(catch_aggresnss, effNhigh, sample = FALSE)

# get catch lengths and wtage at every timestep
catch_length_census_ss <- calc_age2length(structn = catch_structnss,
                                 resn = catch_resnss,
                                 nums = catch_numssshigh,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = 150,
                                 CVlenage = 0.1, remove.zeroes=TRUE)

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

```

This is true catch at age for sardine:

```{r catage1}

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

catchage <- catch_numssshigh %>%
  filter(species == "Pacific_sardine")

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

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

```

Do true catch numbers times weight at age match catch in biomass? It was about 16 times too high relative to catch.txt output, and shifted one year back (which is due to the floor function, and an easy fix). 

June 6: just corrected the catch nums correction in `run_truth` and I hope it matches now.  

```{r catchNtoB}

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

totcatchann <- totcatchN %>%
  mutate(yr = floor(time/fstepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(anncatchage = sum(atoutput))

plottotC <- ggplot(totcatchann, aes(x=yr, y=anncatchage)) +
  geom_point() +
  theme_tufte() +
  labs(subtitle = paste("Catch in N ", scenario.name,
                        totcatchann$species))

plottotC

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

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

CwtAtAge <- catch_lengthwt_census_ss$muweight %>%
  select(species, agecl, time, wtAtAge = atoutput) %>%
  mutate(wtAtAge = wtAtAge/1000)

catchbioN_ss <- merge(catch_numssshigh,CwtAtAge,by=c("species","agecl", "time"),all.x=T) %>%
  mutate(biomass = atoutput*wtAtAge/1000) %>%
  mutate(yr = floor(time/fstepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(anncatchB = sum(biomass))


plotcatch <- ggplot() +
  geom_line(data=truecatchbio_ss, aes(x=time,y=atoutput, color="true catch bio from txt"),
            alpha = 10/10) +
  geom_point(data=catchbioN_ss, aes(x=yr,y=anncatchB, color="true catch bio from N"),
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

```

Sadly, still a problem with the catch correction if it should match catch.txt when comverted to biomass. We are now off by about a factor of 2 (catch nums * wtAtAge is half true catch bio from the txt file):

```{r catchoff}

 sardcatchB <- truecatchbio_ss #%>%
# now done above
#   mutate(yr = as.integer(time/365))

sardcatchBcomp <- merge(sardcatchB, catchbioN_ss)
sardcatchBcomp <- mutate(sardcatchBcomp, mult=atoutput/anncatchB)

mean(sardcatchBcomp$mult)

summary(sardcatchBcomp$mult)

```


## “Data” from OM for EMs

Specify our survey sampling. This could come from other information such as the overlap of actual survey stations with OM polygons, experiments evaluating survey selectivity and efficiency, actual sample-based survey cv, etc. 

Here we are using Christine's specifications, which approximate the acoustic survey index used in the sardine assessment:

```{r sardine-survey-spec}
# from atlantisom/config/sardine_config.R, to be generalized here

#Set species name
species <- species_ss #"Pacific_sardine" #change to species_ss to match 

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

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

#Name of atlantis model output - this is an RData object from the runtruth() function
#atlantis_output_name <- "outputCCV3run_truth.RData"

#Efficiency
eff_case <- 0.5

#Age classes of your species
age_classes <- 1:10

#Selectivity function - a vector equivalent to length (age classes)
sel_by_age <- rep(1,length(age_classes))

#Atlantis model timestep corresponding to the true output--now from census_spec.R
timestep <- stepperyr #5

#Which atlantis timestep does the survey run in?--now from census_spec.R
# with 5 output steps per year, 0 is Jan-Feb-midMar, 1 is midMar-Apr-May, 
# 2 is June-July-midAug, 3 is midAug-Sept-Oct, 4 is Nov-Dec (ish)

survey_sample_time <- midptyr #3 this is more late summer-fall, changed to 2

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of survey times to pull
survey_sample_full <- seq(survey_sample_time,
                          total_sample, by=timestep)

#Effective sample size for composition data (annual total samples)
surveyEffN <- 1000
fisheryEffN <- 1000

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

#Number of years of data to pull
nyears <- 50

#Atlantis initialization period in years
burnin <- 30

#Maximum size of a fish in cm
maxbin <- 150

#Years to fish, assuming fishing output is annual and years index 1,2,3 etc
#fish_years <- burnin:(burnin+nyears-1)
# fishery output timestep (fstepperyr in census_spec.R) for this run is 5 per year

# here we are summing the catch in numbers for input in to SS.
# To use the total catch biomass output from catch.txt the above would work 
# because the output is annual but in units of DAYS so needs time/365

#The last timestep to sample
total_sample <- noutsteps-1 #495

#Vector of indices of catch in numbers to pull (by timestep to sum)
fish_sample_full <- c(0:total_sample)
fish_burnin <- burnin*fstepperyr+1
fish_nyears <- nyears*fstepperyr
fish_times <- fish_sample_full[fish_burnin:(fish_burnin+fish_nyears-1)]
fish_timesteps <- seq(fish_times[5], max(fish_times), by=fstepperyr) #last timestep
fish_years <- unique(floor(fish_times/fstepperyr)+1)

# potential problem: while they are from the same model year,
# fish_years dont match survey_years timesteps (survey mid year, fishery end)
# so when matching use the index, not the value itself

#Years to survey, assuming survey is once per year
survey_years <- survey_sample_full[(burnin+1):(burnin+nyears)]

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

#Efficiency of the survey
effic <- data.frame(species=species, efficiency=eff_case)

# Assume uniform selectivity, same as selex1 here
sel<- data.frame(species=rep(species,length(age_classes)),
               agecl=age_classes,
              selex=sel_by_age)

```

Get input survey biomass for SS (these have q, selectivity, cv)

```{r toSS-surveyts-comps}

#Sample survey - 3rd timestep simulates a summer survey
survey_out <- create_survey(dat=truth$nums, 
                            time=survey_sample_full, 
                            species=species, 
                            boxes=boxall, 
                            effic=effic, 
                            selex=sel)

#Try a biomass based survey for comparison
survey_outB <- create_survey(dat=truth$biomass_ages,
                            time=survey_sample_full,
                            species=species,
                            boxes=boxall,
                            effic=effic,
                            selex=sel)


#Set effective sample size for age compositions
effN <- surveyEffN
highEffN <- data.frame(species=species, effN=rep(effN, length(species)))

#Sample fish for age composition
age_comp_data <- sample_fish(survey_out, highEffN)

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

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

ss_structnstd <- sample_fish(survey_aggstructnstd,
                             effN,
                             sample=FALSE)
ss_resnstd <- sample_fish(survey_aggresnstd,
                          effN,
                          sample=FALSE)

#Extract length composition data
ss_length_stdsurv <- calc_age2length(structn = ss_structnstd,
                                     resn = ss_resnstd,
                                     nums = age_comp_data,
                                     biolprm = truth$biolprm, fgs = truth$fgs,
                                     CVlenage = CVs$lenage, remove.zeroes=TRUE)

#Need to replace with interp function
wtAtAge <- ss_length_stdsurv$muweight %>%
  select(species, agecl, time, wtAtAge = atoutput) %>%
  mutate(wtAtAge = wtAtAge/1000)

# CV function
cv <- data.frame(species="Pacific_sardine", cv=CVs$survey)

#Sample survey biomass
survObsBiom <- sample_survey_biomass(dat=survey_out,cv=cv,wtAtAge)

# check against survey with truth$biomass_ages output
wtage1 <- data.frame(species=rep(species_ss, each=10),
                    agecl=rep(c(1:10),length(species_ss)),
                    wtAtAge=rep(1000.0,length(species_ss)*10))

survObsBiomB <- sample_survey_biomass(dat=survey_outB,cv=cv,wtage1)

# survey numbers, not sure which SS needs?
survObsNum <- sample_survey_numbers(dat=survey_out,cv=cv)


```

Do some comparisons of biomass surveys with truth. 

```{r plotBseries}

plotB <-ggplot() +
  geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
             alpha = 10/10) + 
  geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="survey census B"), 
            alpha = 10/10) +
  geom_point(data=survObsBiomB, aes(x=time/stepperyr,y=atoutput, color="survey sample B from biomass_ages"), 
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

plotB <-ggplot() +
    geom_point(data=atBtxt2tidy, aes(x=Time/365,y=biomass, color="txt output true B"),
             alpha = 10/10) + 
    geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="survey census B"), 
            alpha = 10/10) +
  geom_point(data=survObsBiom, aes(x=time/stepperyr,y=atoutput, color="survey sample B from nums * wtAtAge"), 
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

Get input total catch for SS (here using numbers). SS takes a fishery cv which can be applied here to get a catch time series, or we can input true catch.

PLEASE DON'T USE THIS FOR CCA. We have determined above that catch in numbers from legacy atlantis code, even corrected, is still off by a factor of 2. We leave the info here because this code can be used for atlantis codebases after December 2015. 

```{r toSS-fisheryts}

# catch at age by area and timestep
catch_numbers <-  create_fishery_subset(dat = truth$catch,
                                         time = fish_times,
                                         species = species,
                                         boxes = boxall)

if (fstepperyr>1){
  #need to take all output time steps and sum catch over each year
  catch_numbers_area <- catch_numbers %>%
    #assign a temporary year column, aggregate nums over this year
    mutate(yr = floor(time/fstepperyr)+1) %>%
    group_by(species, agecl, polygon, yr) %>%
    summarise(anncatch_area=sum(atoutput))
  
  # catch in numbers summed over polygon 
  # and agecl after time has been aggregated into 1 yr
  catch_total <- catch_numbers_area %>%
    group_by(yr) %>%
    summarise(catch=sum(anncatch_area))
}

# catch in numbers summed over polygon and agecl if output is annual
if (fstepperyr==1){
  catch_total <- catch_numbers %>%
    group_by(time) %>%
    summarise(catch=sum(atoutput))
}

# if we want to apply the CV to this we need to write sample_total_catch

#cv <- data.frame(species="Pacific_sardine", cv=CVs$fishery)

#catchObsNum <- sample_total_catch(dat=catch_numbers, cv=cv)

```

Check this against true total catch numbers above. Should be the same (and it is, but both are incorrect in the absolute scale due to the CATCH.nc output problem in older codebases).

```{r compare-totcatch}

compareC <- ggplot() +
  geom_line(data=totcatchann, aes(x=yr, y=anncatchage, color="census catch N"), 
            alpha = 10/10) +
  geom_point(data=catch_total, aes(x=yr,y=catch, color="sample catch N"), 
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

```

Get composition inputs for SS (survey and fishery catch at age, survey and fishery lengths, survey and fishey weight at age).

Because catch composition information goes in to the assessment as a proportion, we can use fishery catch at age from this legacy codebase even with absolute catch numbers likely half what they should be overall.


```{r toSS-fisherycomps}

# Survey length comps and wtage done above to get survey ts using nums*wtage approach

#We end up using CAAL for the survey below, so let's generate fishery age comps instead
effN <- fisheryEffN/fstepperyr
effN <- data.frame(species=species, effN=effN)

#catch at age each timestep summed over polygons
catch_numsss_samp <- sample_fish(catch_numbers, effN)

#Get weights
# aggregate true resn per fishery subset design
catch_aggresnss <- aggregateDensityData(dat = truth$resn,
                                 time = fish_times,
                                 species = species,
                                 boxes = boxall)

# aggregate true structn fishery subsetdesign
catch_aggstructnss <- aggregateDensityData(dat = truth$structn,
                                 time = fish_times,
                                 species = species,
                                 boxes = boxall)

#dont sample these, just aggregate them using median
catch_structnss_samp <- sample_fish(catch_aggstructnss, effN, sample = FALSE)

catch_resnss_samp <-  sample_fish(catch_aggresnss, effN, sample = FALSE)

# these fishery lengths and weight at age are each output timestep
catch_lengthwt_samp <- calc_age2length(structn = catch_structnss_samp,
                                 resn = catch_resnss_samp,
                                 nums = catch_numsss_samp,
                                 biolprm = truth$biolprm, fgs = truth$fgs,
                                 maxbin = maxbin,
                                 CVlenage = CVs$lenage, remove.zeroes=TRUE)


```

<!--Can we match catch in weight from catch.txt? No, numbers sampled are 1e-6 compared with actual catch numbers.

```{r compareC-wt, eval=FALSE}

CswtAtAge <- catch_lengthwt_samp$muweight %>%
  select(species, agecl, time, wtAtAge = atoutput) %>%
  mutate(wtAtAge = wtAtAge/1000)

catchsampbio_ss <- merge(catch_numsss_samp,CswtAtAge,by=c("species","agecl", "time"),all.x=T) %>%
  mutate(biomass = atoutput*wtAtAge/1000) %>%
  mutate(yr = floor(time/fstepperyr)+1) %>%
  group_by(species, yr) %>%
  summarise(anncatchBsamp = sum(biomass))

plotcatch <- ggplot() +
  geom_line(data=truecatchbio_ss, aes(x=time/365,y=atoutput/1000000, color="true catch bio"),
            alpha = 10/10) +
  geom_point(data=catchsampbio_ss, aes(x=yr,y=anncatchBsamp, color="sample catch B"), 
            alpha = 10/10) +
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

```
-->

Get population parameters for SS from Atlantis

```{r toSS-poppars}

```

## Write data for Stock Synthesis

To read in the Stock Synthesis files, first ensure you have the latest version of `r4ss`. You will need to ensure the relevant Stock Synthesis files are stored in the inst/extdata/ folder, with each species example in its own folder. 

```{r readdata}
#devtools::install_github("r4ss/r4ss", dependencies = FALSE)
require(r4ss)
#source("./config/sardine_config.R") this is above for illustration

stocksynthesis.data <- r4ss::SS_readdat_3.30(paste0("./inst/extdata/",
model_dir,
datfile_name))
```

*Now changing catch units to tons!*

Write time series biomass (tons) and catch (tons) data for SS using `SS_write_ts`:

```{r SS-write-ts}


#Test writing CPUE in biomass
  stocksynthesis.data <- 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))

stocksynthesis.data$CPUE

stocksynthesis.data$catch
```

Writing age composition data 

```{r SS-write-comps}
#age_comp_data

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

age_comp_flat <- reformat_compositions(age_comp_data,
                     round.places = 4,
                     comp_type = "agecomp")

## Write age composition data for survey
stocksynthesis.data <- 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


```

Once we have the natlength matrix, we can munge the data into the proper CAAL and length bin format for SS

```{r SS-comps-formatting, eval=FALSE}

#Check length compositions match age compositions
#ss_length_stdsurv$natlength %>%
#  group_by(time,agecl) %>%
#  summarise(sum(atoutput))

len_comp_data <- ss_length_stdsurv$natlength
fish_len_comp_data <- catch_lengthwt_samp$natlength

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)) %>%
    summarise(annnatlength=sum(atoutput)) %>%
    rename(atoutput = annnatlength)
} else {
  fish_len_comp_anndata <- fish_len_comp_data
}

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

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

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

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


# below fails with two errors (only gets to first when running)

# Error in `[<-.data.frame`(`*tmp*`, indices, "Seas", value = c(1, 1, 1, : replacement has 49 rows, data has 2045

# Error in `[.data.frame`(comp_matrix[[i]], , as.character(bins[[i]])) : undefined columns selected

# both these errors come from the fishery age and length comps when I step
# through the function, survey caal and lengths read in fine

# Write CAAL and length composition data
stocksynthesis.data <- 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

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

## Run Stock Synthesis
After the data and control file modifications, ensure you copy the local SS 3.3 executable into the folder containing the model files. Then run the following chunk.

```{r runSS, eval=FALSE}
#Run stock synthesis model
run_stocksynthesis(species = species, model_dir = model_dir, show_output = TRUE)

```