Visualize_atlantisom.Rmd

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

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

## Relationships between `atlantsom` functions

First attempt is using the package [`DependenciesGraphs`](https://github.com/datastorm-open/DependenciesGraphs).

Other possibilities:

[`visNetwork`](https://github.com/datastorm-open/visNetwork)

[`drake`](https://docs.ropensci.org/drake/)

```{r installation, eval=FALSE}

devtools::install_github("datastorm-open/DependenciesGraphs")
devtools::install_github("datastorm-open/visNetwork") # on CRAN, this is dev version

```

## `DependenciesGraphs`

This is all the functions in `atlantisom`:

```{r depgraph, message=FALSE}

library(DependenciesGraphs)
library(atlantisom)

dep <- envirDependencies("package:atlantisom")
plot(dep)

```

This is just the `om_init` wrapper function:

```{r depgraphominit}

depominit <- funDependencies("package:atlantisom", "om_init")
plot(depominit)

```

While these map dependencies between functions, it would be nice to see what output of atlantis goes into atlantisom and what atlantisom outputs into SS3.

Try something different. Drake?

```{r atlantisomdrake}

#atlantis output files needed

#atlantisom functions (wrappers) with user-defined inputs specifying sampling frames, observation errors

#om_init oututs run_truth results and model definitions for all species

#om_species cuts full output to user selected focal species--SAVE THIS OUTPUT

#om_index takes user defined survey and fishery sampling and observation error specifications to make replicate indices, outputs both to session and saves

#atlantisom core functions called from wrappers, can be used alone

#stock assessment data rds file output

#atlantisom functions to write assessment model-specific input files

#future: wrapper functions to run and compare assessments with each other and truth

```

So I wound up doing this by hand using online flowchart software Lucidchart (see original [here](https://app.lucidchart.com/documents/edit/7af8dad3-0bc6-4000-9668-9bc1add809ba/0_0)):

Get the truth:
![](docs/images/atlantisomDataFlow_truth-link.png)

Generate the data:
![](docs/images/atlantisomDataFlow_data.png)

## Visualize results using wrappers

Meanwhile lets do final plots for the ecosystem MSE paper. Want a 3 panel plot showing atlantisom sample output (our really good survey), SS3 sample output (fit to survey), and the key feature, SS3 fit to atlantis truth:

From truth and saved wrapper output generated [here](https://sgaichas.github.io/poseidon-dev/WrapperDevelopment.html):

```{r truth-surveyind}

library(tidyverse) #fix this; build in this dependency for wrappers
library(here)
library(ggthemes)

# skip this, here for completeness, just need directory
#CC3om <- om_init(here("config/CC3config_sardLWcorr.R")) #initialize/read in truth

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

#skip this too and now just read in the saved results
#CC3om_sardine <- om_species(c("Pacific_sardine"), CC3om) #select just sardine

# TODO make option to save this output so we dont have to run om_init to get truth 
# every time and it doesnt sit around in memory. For now:

#saveRDS(CC3om_sardine, file.path(d.name, paste0(scenario.name, "omlist_ss.rds")))

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

# we have already run this so just read in the saved results, here for reference
#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_ind <- readRDS(file.path(d.name, paste0(scenario.name, "surveyB.rds")))

# TODO make plotting functions for atlantisom that do these comparisons

# needed for plotting time dimensions 
# TODO standardize time across outputs in these wrappers
omlist_ss <- CC3om_sardine
source(here("config/omdimensions.R"))

txtTrueB <- CC3om_sardine$truetotbio_ss

#WARNING! this is true bio at age so needs to be summed over ages--rename truebioage_ss
omTrueB <- CC3om_sardine$truebio_ss %>%
  group_by(species, time) %>%
  summarise(totB = sum(atoutput)) 

survObsBiomB <- CC3om_sard_ind[[1]]


plotB <-ggplot() +
  geom_point(data=txtTrueB, aes(x=time/365,y=atoutput, color="txt output true B"),
             alpha = 10/10) + 
  geom_line(data=omTrueB, aes(x=time/stepperyr,y=totB, 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")


```


From Christine's ggplot_SS.R (change to Isaac's output directory name for current best fit) we get SS3 sample output using a prettier version of r4ss plot:

```{r ggplotSS, results="hide"}

#dir <- "./inst/extdata/Sardine_SS_Files/"

dir <- "./inst/extdata/AB_Lmax25_v2_ABinControl_Mlower/"
  
#Get SS output
output_from_SS <- r4ss::SS_output(dir, covar=F, forecast=F) 

#Extract index info
index_data_fits <- output_from_SS$cpue
#Commented out - only for comparing to SS output to check
#r4ss::SS_plots(output_from_SS)

#Plot
#require(ggthemes)
#require(ggplot2)
plotB <-ggplot(index_data_fits) +
  geom_line(data=index_data_fits, aes(x=Yr,y=Exp,color="Estimated CPUE"), 
            alpha = 10/10) +
  geom_pointrange(data=index_data_fits, aes(x=Yr,y=Obs, ymin=exp(log(Obs)+qt(0.025,df=output_from_SS$survey_error[2])*SE), 
                                                                 ymax=exp(log(Obs)+qt(0.975,df=output_from_SS$survey_error[2])*SE),
                                                                          color="Input CPUE with error"),
             alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") 
plotB


```

Add the comparison of SS3 estimated total biomass with true atlantisom tot B as started [here](https://sgaichas.github.io/poseidon-dev/SkillAssessInit.html):

```{r compareB}

#point this to Isaac's output dir defined above
#replist <- r4ss::SS_output(dir = here(dir), verbose=TRUE, printstats=TRUE,covar=FALSE)

#SS_plots(replist)

#replist$timeseries
names(output_from_SS$timeseries)[19]<-"Fmort"

plotB <-ggplot() +
  geom_line(data=omTrueB, aes(x=time/stepperyr,y=totB, color="True B"), 
            alpha = 10/10) +
  geom_point(data=output_from_SS$timeseries, aes(x=Yr,y=Bio_all, color="SS3 Est B"),
             alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

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

```

Putting it together into a 3 panel figure:

```{r triplot, fig.asp=1}

require(cowplot)
require(viridis)

mycol <- viridis(4)

plotsurv <- ggplot() +
  geom_line(data=omTrueB, aes(x=time/stepperyr,y=totB, color="True total B"), 
            alpha = 10/10) +
  geom_point(data=survObsBiomB, aes(x=time/stepperyr,y=atoutput, color="Survey sample B"), 
            alpha = 10/10) +
  xlab("") + 
  ylab("") +
  scale_color_manual(values=c(mycol[2], mycol[1]),
                     guide = guide_legend(override.aes = list(
                         linetype = c("blank", "solid"),
                         shape = c(16, NA)))) + 
  xlim(20,80) +
  theme_tufte() +
  theme(legend.position = "right") +
  labs(colour="Atlantis-sardine")

plotSSfit <- ggplot() +
  geom_line(data=index_data_fits, aes(x=Yr,y=Exp,color="SS3 est. survey B"), 
            alpha = 10/10, size=2) +
  geom_pointrange(data=index_data_fits, 
                  aes(x=Yr,y=Obs,
                      ymin=exp(log(Obs)+qt(0.025,df=output_from_SS$survey_error[2])*SE),
                      ymax=exp(log(Obs)+qt(0.975,df=output_from_SS$survey_error[2])*SE),
                      fatten = .5,
                      color="Survey B + error"),
                  alpha = 10/10) + 
  scale_x_continuous(limits = c(0,120), breaks = seq(0,125,25)) +
  xlab("") + 
  ylab("Biomass, t") +
  scale_color_manual(values=c(mycol[3], mycol[2]),
                     guide = guide_legend(override.aes = list(
                         linetype = c("solid", "blank"),
                         shape = c(NA, 16)))) + 
  xlim(20,80) +
  theme_tufte() +
  theme(legend.position = "right") +
  labs(colour="SS3-survey fit")

plotSSskill <- ggplot() +
  geom_line(data=omTrueB, aes(x=time/stepperyr,y=totB, color="True total B"), 
            alpha = 10/10) +
  geom_point(data=output_from_SS$timeseries, aes(x=Yr,y=Bio_all, color="SS3 est. total B"),
             alpha = 10/10) + 
  xlab("Year") + 
  ylab("") +
  scale_color_manual(values=c(mycol[3], mycol[1]),
                     guide = guide_legend(override.aes = list(
                         linetype = c("blank", "solid"),
                         shape = c(16, NA)))) + 
  xlim(20,80) +
  theme_tufte() +
  theme(legend.position = "right") +
  labs(colour="SS3-biomass skill")

plot_grid(plotsurv + theme(legend.justification = c(0,1)), 
          plotSSfit + theme(legend.justification = c(0,1)), 
          plotSSskill + theme(legend.justification = c(0,1)), 
          ncol=1, align="v")

```