CreateStockSynthesis.Rmd

---
title: "WritingSSfiles"
author: "Christine Stawitz"
date: "May 14, 2019"
output: html_document
---

```{r setup, include=FALSE}
require(here)
require(reshape2)
require(dplyr)
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "C:/Users/chris/Documents/GitHub/atlantisom/")
```

## Introduction

This page discusses reading in an existing stock assessment file and replacing the data values with those from the `atlantisom` package. 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 packages, eval=FALSE}
devtools::install_github("r4ss/r4ss", dependencies = FALSE)
```

## Sardine example

Dummy files are in the package directory; use `r4ss` to load the dummy data file. Use the relevant species config file to load the configuration variables for each run. Note: this script will extract conditional age-at-length data, length composition data, and a CPUE in biomass for one survey during the specified time interval, and age composition data, length composition data, and catch in numbers from a single fishery. The timing and duration of the survey and fishery should be set in the config file. Unless otherwise stated, numbers going into Stock Synthesis are in thousands of fish. Weight of a single fish is in grams or kilograms, and weight of biomass is in metric tons. We assume a single sex model with a sex ratio fixed at 0.5 - different assumptions require modifying the Stock Synthesis control file. This script only runs with Stock Synthesis 3.30.   

```{r readdata, echo=FALSE}
require(r4ss)
source("./config/CCConfig.R")
source("./config/census_spec.R")
source("./config/sardine_config.R")

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

Below, we can see the structure of the catch and CPUE data
```{r}
names(stocksynthesis.data$catch) <- names(stocksynthesis.data$CPUE)

stocksynthesis.data$CPUE
```

We can next look at age composition data
```{r}
stocksynthesis.data$agecomp

stocksynthesis.data$lencomp
```

We may also need to access or modify the setup variables before the data matrices.

```{r}
#CPUE units, SD and error structure
stocksynthesis.data$CPUEinfo

#Number of length bins
stocksynthesis.data$N_lbins

#Length bin method
stocksynthesis.data$Lbin_method

#Vector of age bins
stocksynthesis.data$agebin_vector

```

## Writing index data
To edit index data, use the `atlantisom` sampling functions to extract survey biomass and numbers. These can then be written using the `SS_write_CPUE` function.

```{r}

if(!file.exists(file.path(d.name, 
                          paste0("output", scenario.name, "run_truth.RData")))){
  #Store all loaded results into an R object
  results <- run_truth(scenario = scenario.name,
                     dir = d.name,
                     file_fgs = functional.groups.file,
                     file_bgm = box.file,
                     select_groups = species,
                     file_init = initial.conditions.file,
                     file_biolprm = biol.prm.file,
                     file_runprm = run.prm.file,
                     verbose = TRUE
  )
} else {
  results <- get(load(file.path(d.name,
                              paste0("output", scenario.name, "run_truth_all.RData"))))
}
results <- get(load(atlantis_output_name))
devtools::load_all()
species=c(species)
truth <- results

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

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

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

len <- lengthwt_census_ss$natlength

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

#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_testBall, highEffN)
boxes <- unique(results$nums$polygon)

# aggregate true resn per survey design
survey_aggresnstd <- aggregateDensityData(dat = results$resn,
                                          time = unique(results$nums$time),
                                          species = species,
                                          boxes = boxes)

# aggregate true structn per survey design
survey_aggstructnstd <- aggregateDensityData(dat = results$structn,
                                             time = unique(results$nums$time),
                                             species = species,
                                             boxes = boxes)

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 = results$biolprm, fgs = results$fgs,
                                     CVlenage = CVs$lenage, remove.zeroes=TRUE)


saveRDS(ss_length_stdsurv, file.path(d.name, paste0(scenario.name, "lengthwt_census_sard.rds")))

#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_testBall,cv=cv,wtAtAge)



```


## Catch data

```{r}
fish_years <- annualmidyear[fish_years]
catch_numbers <-  create_fishery_subset(dat = results$catch,
                                         time = fish_years,
                                         species = species,
                                         boxes = boxes)


#Sample catch in numbers
catch_total <- catch_numbers %>%
  group_by(time) %>%
  summarise(catch=sum(atoutput))

#Use .txt file output biomass instead
names(results$catch_all)[2] <- "Species"
catch_total <- filter(results$catch_all, Species==species) %>%
  filter(time %in% (fish_years+8))

#Test writing CPUE in biomass
  stocksynthesis.data <- SS_write_ts(ss_data_list = stocksynthesis.data,
                ts_data = list(survObsBiom$atoutput[fish_years],
                  catch_total$catch),
                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))

```

## Writing age composition data 

```{r}
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

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

catch_numssshigh <- sample_fish(catch_numbers, effN)

#Get weights
# aggregate true resn per survey or fishery subset design
catch_aggresnss <- aggregateDensityData(dat = results$resn,
                                 time = c(0:111),
                                 species = species,
                                 boxes = boxes)

# aggregate true structn per survey or fishery subsetdesign
catch_aggstructnss <- aggregateDensityData(dat = results$structn,
                                 time = c(0:111),
                                 species = species,
                                 boxes = boxes)

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

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

str(catch_numssshigh)

catch_length <- calc_age2length(structn = catch_structnss,
                                 resn = catch_resnss,
                                 nums = catch_numssshigh,
                                 biolprm = results$biolprm, fgs = results$fgs,
                                 maxbin = maxbin,
                                 CVlenage = CVs$lenage, remove.zeroes=TRUE)

```


## Setting up EM
For the EM, we want conditional age-at-length (CAAL) and length composition data only. So we use the same function with our matrices of length and CAAL comp data to write to the data file.

```{r}

#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_length$natlength

```

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


```{r}

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

fish_len_comp_flat <- reformat_compositions(catch_length$natlength,
                                  round.places = 0,
                           comp_type="lencomp")
#remove burnin
fish_len_comp_final <- filter(fish_len_comp_flat,
                         time %in% fish_years)

fish_age_comp_flat <- reformat_compositions(catch_numssshigh,
                           comp_type="agecomp")
#remove burnin
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]])))

# 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(stocksynthesis.data, outfile = paste0(d.name, "/inst/extdata/",model_dir,
datfile_name),
                 overwrite=TRUE)
```

## Coming up with dynamic life history parameters

Below we get the weight-length relationship parameters from atlantis, calculate h and R0 from Atlantis $\alpha$ and $\beta$ parameters, use the \code{calc_Z} function to back out natural mortality (M) from Atlantis, and use the survey length-at-age from Atlantis to estimate a growth curve.
```{r}
#Get biological parameters

#Load needed inputs for biological parameters
source(here("config/CC1Config.R"))
#biological parameters
biolprm <- load_biolprm(dir = "C:/Users/chris/Documents/GitHub/atlantisom", file_biolprm = biol.prm.file)
d.name <- "C:/Users/chris/Documents/GitHub/atlantisom"
#functional groups
fgs <- results$fgs
runprm <- load_runprm(dir = d.name, file_runprm = run.prm.file)
YOY <- load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))
truenums_ss <- results$nums[results$nums$species==species,]

YOY_ss <- YOY %>%
  select(Time, "SAR.0")

fullresZ <- calc_Z(yoy = YOY_ss,
                   nums = truenums_ss,
                   fgs = fgs,
                   biolprm = biolprm,
                   toutinc = runprm$toutinc)


meanwt_spp <- ss_length_stdsurv$muweight %>% 
  filter(time>burnin) %>%
  group_by(agecl) %>%
  summarize(meanwt = mean(atoutput))
sardine.ctl <-r4ss::SS_readctl(paste0("./inst/extdata/",
model_dir,
ctlfile_name))
sardine.ctl <- SS_write_biol(sardine.ctl, biolprm, "SAR", Z = mean(fullresZ$atoutput), wtsage=meanwt_spp)

li_a_use <- biolprm$wl[match(fgs$Code[match(species,fgs$Name)],biolprm$wl[, 1]), 2]/1000
li_b_use <- biolprm$wl[match(fgs$Code[match(species,fgs$Name)],biolprm$wl[, 1]), 3]

#maturity ogive
#check what R0 is



len_nonburn <- ss_length_stdsurv$mulen %>%
  filter(time > burnin*timestep)
plot(len_nonburn$atoutput~len_nonburn$agecl, col=len_nonburn$time)

length.data <- data.frame("Year"=(len_nonburn$time-survey_sample_time)/timestep+1, length=len_nonburn$atoutput, Weight=NA, Sex="Female", age=as.integer(len_nonburn$agecl))

#require(bbmle)
#vb_estimates <-sample_fit_vbgf(length.data, 25, 45, 0.4,
#  0.1, 0.1, 20, 40, 0.05, 0.01, 0.01,
#  30, 50, 0.6, 0.3, 0.3, 0.5, 10)

```

## Running 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}
#Run stock synthesis model
run_stocksynthesis(species = species, model_dir = model_dir, show_output = TRUE)

#Add a call here to plot output 

```