# Sample R codes for estimating parameters using catch and effort data, and conduct population dynamics simulations with the minimum AIC models

#--- Read functions
source("functions.R")

#--- create excutable file of ADMB
if(.Platform$OS.type=="unix"){
  system(paste("admb -r gamma_ss4"))
}
if(.Platform$OS.type=="windows"){
  shell(paste("admb -r gamma_ss4"))
}

#--- read VPA data from vpa.txt
vpa.res <- list()
vpa.res$naa <- read.table("vpa.txt",skip=1,nrows=7,header=T)
vpa.res$maa <- read.table("vpa.txt",skip=10,nrows=7,header=T)
vpa.res$M <- read.table("vpa.txt",skip=19,nrows=7,header=T)
vpa.res$saa <- read.table("vpa.txt",skip=28,nrows=7,header=T)
vpa.res$biom <- read.table("vpa.txt",skip=37,nrows=7,header=T)
vpa.res$rps <- read.table("vpa.txt",skip=46,nrows=1,header=T)
vpa.res$acatch <- read.table("vpa.txt",skip=49,nrows=1,header=T)
vpa.res$vssb <- read.table("vpa.txt",skip=52,nrows=1,header=T)
vpa.res$vbiom <- read.table("vpa.txt",skip=55,nrows=1,header=T)
vpa.res$waa <- read.table("vpa.txt",skip=58,nrows=7,header=T)
vpa.res <- lapply(vpa.res, function(x){
                   names(x) <- as.character(2003:2009); x
                 })

#--- read fishery data (dummy data)
# NOTE THAT, because this data is a dummy data created by a stochastic simulation using the estimated parameters, the estimated parameters from this data are different from results in the paper. 
ndata1 <- read.csv("ndata1.csv") #dummy data

#--- Parameter estimation (minimum AIC model)
# The explanatory variables of m1 (ship model) and m3 (catch model) were determined by manual model selection.
# The others of GLM are automatically selected by the function of 'step'.

# model structure
m1 <- expression(cbind(N1 - 1, max.N1 - N1) ~ log(biomass.c) + log(biomass.p) +  factor(fm)) # ship model (the AIC minimum model)
m2 <- expression(mcatch.wt > 0 ~ log(biomass.c)+log(biomass.p)+factor(fm)+pK+pW) # 0/1 model (model selection is done the R-function of est.coef)
m3 <- expression(mcatch.wt ~ log(N1) + log(biomass.c) + pW + factor(fm)) # catch model (the AIC minimum model)
m4 <- expression(K ~ D) # closure model A  (model selection is done the R-function of est.coef)
m5 <- expression(K ~ D+pK+X2) # closure model B  (model selection is done the R-function of est.coef)

# estimate coefficients for the above five models 
# In the main text, we used parameters from averaged models, but this R-code is for the calculation based on parameters of minimum-AIC mode for simplicity. 
ss.glm <- est.coef(ndata1,
                   models=list(formula=c(m1,m2,m3,m4,m5),
                     is.stepAIC=c(FALSE,TRUE,FALSE,TRUE,TRUE)),
                   admb.mode=list(is.admb=c(TRUE,FALSE,TRUE,FALSE,FALSE),
                     dist=rep("gamma_ss4",5),
                     phase=list(c(1,-1,1,1,1,1),NULL,c(1,1,1,1,1,1)),
                     offset=list(NULL,NULL,NULL)))

#--- simulatin starts from 2004 in fishing year
ndata2 <- subset(ndata1,fy>2003)

# the scenario with day closures
# note; The number of simulations (N) should be more than 11, because of a trivial technical reason
res1a <- make.simu2(ndata2,ss.glm,vpa.res,N=1000,float.dyn=TRUE,admb.mode=c(TRUE,FALSE,TRUE))

# the scenario without day closures
res1n <- make.simu2(ndata2,ss.glm,vpa.res,N=1000,float.dyn=TRUE,admb.mode=c(TRUE,FALSE,TRUE),type="nreg")

#--- show results
# graph
par(mfrow=c(2,2))
show.res(res1a,res1n,ndata2)
# table (corresponding to Table 3 in the main text)
make.table(res1a,res1n)