## ----include = FALSE, warning=FALSE, message=FALSE----------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
start_time = Sys.time()
# Install locally
#  devtools::install_local( R'(C:\Users\James.Thorson\Desktop\Git\ecostate)', force=TRUE )
#  devtools::install_github( 'James-Thorson-NOAA/ecostate', force=TRUE )
# Build
# setwd(R'(C:\Users\James.Thorson\Desktop\Git\ecostate)'); 
# devtools::build_rmd("vignettes/model_of_intermediate_complexity.Rmd"); rmarkdown::render( "vignettes/model_of_intermediate_complexity.Rmd", rmarkdown::pdf_document())

## ----setup, echo=TRUE, message=FALSE------------------------------------------
library(ecostate)

## ----echo=TRUE, message=FALSE, fig.width=7, fig.height=7----------------------
# load data
data(eastern_bering_sea)

# Reformat inputs
years = 1982:2021 # Catch only goes through 2021, and starting pre-data in 1982 doesn't play well with fit_B0
taxa = c( "Pollock", "Cod", "Arrowtooth", "Copepod", "Other_zoop", "Chloro", "NFS", "Krill", "Benthic_invert", "Benthos", "Detritus" )

# Define types
type_i = sapply( taxa, FUN=switch, "Detritus" = "detritus",
                                   "Chloro" = "auto",
                                   "hetero" )

# Starting values
U_i = EE_i = B_i = array( NA, dim=length(taxa), 
                    dimnames=list(names(eastern_bering_sea$P_over_B)))
B_i[c("Cod", "Arrowtooth", "NFS")] = c(1, 0.5, 0.02)
EE_i[] = 1
U_i[] = 0.2

# Define default vulnerability, except for primary producers
X_ij = array( 2, dim=c(length(taxa),length(taxa)) )
dimnames(X_ij) = list(names(B_i),names(B_i))
X_ij[,'Chloro'] = 91

## ----echo=TRUE, message=FALSE, fig.width=7, fig.height=7----------------------
# Define parameters to estimate
fit_Q = c("Pollock", "Copepod", "Chloro", "Other_zoop", "Krill")
fit_B0 = c("Pollock", "Cod", "Arrowtooth", "NFS")
fit_B = c("Cod", "Arrowtooth", "NFS")  

# Define process errors to estimate
# Only estimating Pollock to speed up demonstration
fit_eps = "Pollock"

# Which taxa to include
taxa_to_include = c( "NFS", "Pollock", "Copepod", "Chloro", "Krill" )
# To run full model use:
# taxa_to_include = taxa

# Define priors
log_prior = function(p){
  # Prior on process-error log-SD to stabilize model
  logp = sum(dnorm( p$logtau_i, mean=log(0.2), sd=1, log=TRUE ), na.rm=TRUE)
}

# Run model
out = ecostate( taxa = taxa_to_include,
                years = years,
                catch = eastern_bering_sea$Catch,
                biomass = eastern_bering_sea$Survey,
                PB = eastern_bering_sea$P_over_B,
                QB = eastern_bering_sea$Q_over_B,
                DC = eastern_bering_sea$Diet_proportions,
                B = B_i,
                EE = EE_i,
                U = U_i,
                type = type_i,
                X = X_ij,
                fit_B = fit_B,
                fit_Q = fit_Q,
                fit_eps = fit_eps,
                fit_B0 = fit_B0,
                log_prior = log_prior,
                control = ecostate_control( n_steps = 20, # using 15 by default
                                            start_tau = 0.01,
                                            tmbad.sparse_hessian_compress = 0 ))

# print output
out

## ----echo=TRUE, message=FALSE, fig.width=7, fig.height=7----------------------
# Plot foodweb at equilibrium
# using pelagic producers as x-axis and trophic level as y-axis
plot_foodweb( out$rep$out_initial$Qe_ij,  
              xtracer_i = ifelse(taxa_to_include=="Krill",1,0),
              B_i = out$rep$out_initial$B_i,
              type_i = type_i[taxa_to_include] )

## ----include = FALSE, warning=FALSE, message=FALSE----------------------------
run_time = Sys.time() - start_time