diff --git a/examples/BSAI NRS/ForProjections.dat b/examples/BSAI NRS/ForProjections.dat
new file mode 100644
index 0000000..5dce63e
--- /dev/null
+++ b/examples/BSAI NRS/ForProjections.dat
@@ -0,0 +1,40 @@
+fm_projection_model_output
+0 # SSL Species???
+0 # Constant buffer of Dorn?
+1 # Number of fsheries
+2 # Number of sexes??
+0.0698939 # averagei 5yr f
+1.0 # author f
+0.4 # SPR ABC
+0.35 # SPR MSY
+2 # Spawnmo
+20 # Number of ages
+1 # Fratio
+#females first
+0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.15
+#male
+0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865 0.172865
+# Maturity Females
+0 0 0 0 0.001 0.008 0.061 0.2 0.505 0.75 0.888 0.929 0.955 0.979 0.983 0.991 0.996 0.997 0.997 0.997
+# Maturity Males same as females!!
+0 0 0 0 0.001 0.008 0.061 0.2 0.505 0.75 0.888 0.929 0.955 0.979 0.983 0.991 0.996 0.997 0.997 0.997
+# Wt spawn females
+0 0.0143 0.0426333 0.0880667 0.135333 0.249269 0.36119 0.4164 0.459727 0.488625 0.525333 0.61875 0.6755 0.7158 0.72725 0.629667 0.725308 0.7105 0.79275 0.927545
+# WtAge Females, by fishery
+0 0.01 0.070674 0.119675 0.163832 0.243534 0.309472 0.388137 0.450988 0.527343 0.546709 0.566713 0.575558 0.548343 0.568974 0.600366 0.606937 0.644775 0.604785 0.65195
+# WtAge Males, by fishery
+0 0.01 0.0706607 0.115078 0.181837 0.239772 0.280924 0.318739 0.399935 0.37515 0.421513 0.401573 0.314375 0.382681 0.415436 0.401458 0.411747 0.406325 0.449322 0.519749
+# Selectivity Females, by fishery
+0.000311118 0.000842886 0.00228149 0.00616031 0.0165244 0.0435607 0.10989 0.250739 0.475649 0.71089 0.869541 0.947554 0.97999 0.992523 0.997229 0.998976 0.999622 0.999622 0.999622 0.999622
+# Selectivity males, by fishery
+0.00020912 0.000738778 0.00260599 0.00914349 0.031494 0.10209 0.279024 0.547414 0.752056 0.840996 0.870105 0.878709 0.881174 0.881874 0.882073 0.882129 0.882144 0.882144 0.882144 0.882144
+# N at age in endyr Females, Males
+440.236 361.536 1568.42 685.803 390.204 399.757 283.75 145.383 22.4396 36.8494 23.2145 15.9373 6.9824 8.43125 9.7935 20.1012 51.2928 33.37 31.563 87.4915
+440.236 353.362 1498.29 640.245 355.799 355.337 245.094 121.789 18.1604 28.7295 17.4556 11.5895 4.92134 5.79107 6.56302 13.109 32.6172 20.979 19.8522 53.4976
+# No Recruitments
+40
+# Recruitment 1978-2004
+280.595 274.373 360.874 653.349 666.766 615.804 1010.65 874.832 835.076 1414.25 2165.87 751.508 643.514 1442.41 711.898 366.867 542.332 286.298 278.206 383.692 219.987 336.756 314.047 666.997 1037.49 1229.31 923.733 746.552 892.496 278.217 110.018 77.902 53.3355 100.487 120.679 158.777 80.5199 435.205 713.113 852.928
+# SSB
+# used only for S/R analysis
+54.1439 67.006 82.3962 94.4856 100.706 100.453 105.69 119.585 116.998 123.639 148.617 155.766 163.769 177.721 199.031 215.356 272.954 317.17 366.632 450.015 493.62 485.398 509.759 533.741 549.379 557.758 568.449 551.987 470.123 431.805 407.609 388.716 354.036 375.854 429.022 461.09 469.144 473.966 485.95 433.155
diff --git a/examples/BSAI NRS/spm.dat b/examples/BSAI NRS/spm.dat
new file mode 100644
index 0000000..7e077f8
--- /dev/null
+++ b/examples/BSAI NRS/spm.dat
@@ -0,0 +1,55 @@
+#rn
+NRS_test
+#Tier
+3
+#nalts
+7
+#alts
+1
+2
+3
+4
+5
+6
+7
+#tac_flag
+1
+#srr_type
+1
+#srr_form
+1
+#srr_conditioning
+0
+#srr_reserved
+0
+#spm_detail_flag
+1
+#nprj_yrs
+14
+#nsims
+1000
+#beg_yr
+2022
+#nyrs_fixed_catch
+2
+#nspp
+1
+#OY_min
+10
+#OY_max
+1943.248
+#data_files for each species
+../fm.prj
+#ABC_mults
+1
+#scalars
+1
+#alt4_spr
+0.6
+#nTAC_cat
+1
+#TACind
+1
+#fixed_catch
+2022 16.0143 #Beginning of October catch
+2023 40.760 #10 yr average up to 2021.
diff --git a/examples/Misc/Makefile b/examples/Misc/Makefile
new file mode 100644
index 0000000..f3d0cf2
--- /dev/null
+++ b/examples/Misc/Makefile
@@ -0,0 +1,80 @@
+run = ../src/spm -nox -nohess >/dev/null
+data = data/$(stock)_spcat.dat
+outdir = $(stock)_out
+RM = rm
+.PHONY: all run
+
+all:
+ifneq "$(wildcard $(data) )" ""
+ cp $(data) spp_catch.dat
+ @echo $(outdir)
+ $(run)
+ifneq "$(wildcard $(outdir) )" ""
+else
+ mkdir $(outdir)
+endif
+ mv *.out $(outdir)/
+ mv spm.rep $(outdir)/report.out
+ $(RM) rm eigv.rpt variance admodel.* *.r0? *.p0? fmin.log *.b0?
+else
+ # if it doesn't:
+ @echo "Oops...error, file " $(data) " appears to be missing... "
+endif
+
+
+clean:
+ $(RM) $(outdir)/*
+
+
+EXEC = spm
+DIST = ../src/
+ARGS = -nox -iprint 150
+
+ifdef ComSpec
+ RM=del /F /Q
+else
+ RM=rm -rf
+endif
+
+all: mpd $(DIST)$(EXEC).tpl
+
+$(EXEC): $(DIST)$(EXEC).tpl
+ ln -sf $(DIST)$(EXEC) $@
+ ln -sf $(DIST)$(EXEC).tpl $@.tpl
+ $(MAKE) --directory=../../src
+
+$(DIST)$(EXEC).tpl:
+ $(MAKE) --directory=../../src
+
+
+mpd: $(EXEC)
+ ./$(EXEC) $(ARGS)
+ @$(RM) $(EXEC).*[0123456789] *.rpt *.log variance gradient.* *tmp admodel.* *.eva
+
+mcmc: $(EXEC)
+ ./$(EXEC) $(ARGS) -mcmc 3000000 -mcsave 600
+ ./$(EXEC) -mceval
+
+proj: $(PROJ)
+ ./$(PROJ)
+
+debug: $(EXEC)
+ ./$(EXEC) $(ARGS)
+
+plot:
+ R CMD BATCH plot.R
+
+clean:
+ @$(RM) $(EXEC)
+ @$(RM) $(EXEC) $(EXEC).[brces]* $(EXEC).*[0123456789] *.rpt *.log variance gradient.* *tmp
+ @$(RM) admodel.*
+ @$(RM) checkfile.rep
+ @$(RM) mcout.rep
+ @$(RM) plot.Rout
+ @$(RM) Rplots.pdf
+ @$(RM) *.rep
+ @$(RM) Fprof.yld
+ @$(RM) *.prj
+ @$(RM) pm.par
+ @$(RM) SIS_out.rep
+ @$(RM) mceval.dat
diff --git a/examples/Misc/README.md b/examples/Misc/README.md
new file mode 100644
index 0000000..8c47d58
--- /dev/null
+++ b/examples/Misc/README.md
@@ -0,0 +1,17 @@
+# Example projection model run
+
+## To do list...
+
+ - [ ] get datafiles updated in examples directory
+ - [ ] make sure script runs
+ - [ ] write R script for settings and inputs and running
+ - [ ] include post-processing R script
+
+
+## Options
+
+ make mydat
+
+ run.bat mydat
+
+
diff --git a/examples/Misc/SS_example.R b/examples/Misc/SS_example.R
new file mode 100644
index 0000000..be17f39
--- /dev/null
+++ b/examples/Misc/SS_example.R
@@ -0,0 +1,222 @@
+#' # Projection model work
+#' Projection model work for demonstrating application of controls and input data
+
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+dir = "C:/WORKING_FOLDER/Model19.14.48c_T"
+Model_Name = "Model19.14.48c_T"
+
+source("../R/readData.R")
+
+#' ## Set initial "setup" parameters
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC (1 means true, otherwise false)
+ srr = 1 , #' Stock-recruitment type (1=Ricker, 2=Bholt)
+ rec_proj = 1, #' projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+ srr_cond = 0 , #' SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 1000, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+
+#' ## Set up the species specific run file
+config<-list(
+ nFixCatchYrs = 1,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote(paste0("data/",Model_Name,".dat")),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c(2019,15000.)
+)
+
+
+##
+
+write_proj<-function(data_file="Model19.14.48c_T_Proj.dat",
+ data = mod1,
+ FY=1977,
+ LY=2019,
+ RY=2,
+ fleets=3,
+ sexes=1,
+ NAGES=10,
+ SSL=0,
+ Dorn = 0,
+ AUTHOR_F = 1,
+ SPR_ABC = 0.4,
+ SPR_MSY = 0.35,
+ SPAWN_M = 1,
+ FRATIO=noquote("0.3 0.2 0.5"))
+ {
+## writing projection file
+## mean 5 year F
+Y5<-LY-5
+F_5<-mean(data$sprseries$Tot_Exploit[data$sprseries$Yr>Y5&data$sprseries$Yr<=LY])
+## population weight at age for females
+WGT<-vector("list",length=sexes)
+Nage_LY<-vector("list",length=sexes)
+M1<-vector("list",length=sexes)
+for(i in 1:sexes){
+ WGT[[i]]<-(data.table(data$endgrowth)[Sex==i]$Wt_Beg*data.table(data$endgrowth)[Sex==i]$Mat_F_Natage)[2:(NAGES+1)]
+ Nage_LY[[i]]<-subset(data$natage,data$natage[,11]=="B"&data$natage$Yr==LY&data$natage$Sex==i)
+ M1[[i]]<-as.numeric(subset(data$M_at_age,data$M_at_age$Yr==(LY-10)&data$M_at_age$Sex==i)[,4:(NAGES+3)]) ## pulling array of Ms at age by sex for LY-10
+}
+
+## selectivity at age for fishery
+sel_LY<-vector("list",length=fleets)
+wt_LY<-vector("list",length=fleets)
+for(i in 1:fleets){
+ sel_LY[[i]]<-subset(data$ageselex,data$ageselex$Fleet==i&data$ageselex$Yr==LY-RY&data$ageselex$Factor=="Asel2")
+ wt_LY[[i]]<-subset(data$ageselex,data$ageselex$Fleet==i&data$ageselex$Yr==LY-RY&data$ageselex$Factor=="bodywt")
+}
+
+## numbers at age
+Rec_1<-as.numeric(data$natage[,14][data$natage$Yr<=LY&data$natage$Yr>=FY&data$natage$Sex==1&data$natage[,11]=="B"])
+N_rec<-length(Rec_1)
+
+SSB<-as.numeric(data$sprseries$SSB[data$sprseries$Yr<=LY&data$sprseries$Yr>=FY])
+#SSB<-SSB[1:(LY-FY)]
+T1<-noquote(paste(data_file))
+write(T1,paste(data_file),ncolumns = 1 )
+T1<-noquote(paste0(SSL," # SSL Species???"))
+ write(T1,paste(data_file),ncolumns = 1,append=T)
+T1<-noquote(paste0(Dorn," # Constant Buffer Dorn?"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(fleets," # Number of fisheries"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(sexes," # Number of Sexes"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste(F_5,"# Average 5 year F"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(AUTHOR_F," # Author f"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(SPR_ABC," # SPR ABC"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(SPR_MSY," # SPR MSY"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(SPAWN_M," # Spawning month"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(NAGES," # number of ages"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote(paste0(FRATIO," # Fratio"))
+ write(T1,paste(data_file),append = T)
+T1<-noquote("# natural mortality")
+ write(T1,paste(data_file),append = T)
+ for (i in 1:sexes){
+ write(M1[[i]],paste(data_file),append = T,ncolumns = 45)
+ }
+T1<-noquote("# Maturity ")
+ write(T1,paste(data_file),append = T)
+ write(rep(1,NAGES),paste(data_file),append = T,ncolumns = 45) ## Female maturity??
+T1<-noquote("# wt spawn females")
+ write(T1,paste(data_file),append = T)
+ for(i in 1: sexes){
+ write(round(as.numeric(WGT[[i]]),4),paste(data_file),append = T,ncolumns = 45)
+ }
+T1<-noquote("# WtAge females by fishery")
+ write(T1,paste(data_file),append = T)
+ for(i in 1: fleets){
+ write(round(as.numeric(wt_LY[[i]][1,9:(NAGES+8)]),4),paste(data_file),append = T,ncolumns = 45)
+ }
+T1<-noquote("# Selectivity females by fishery")
+ write(T1,paste(data_file),append = T)
+ for(i in 1: fleets){
+ write(round(as.numeric(sel_LY[[i]][1,9:(NAGES+8)]),4),paste(data_file),append = T,ncolumns = 45)
+ }
+
+T1<-noquote(paste0("# Numbers at age females males in ",LY))
+ write(T1,paste(data_file),append = T)
+ for (i in 1:sexes){
+ write(as.numeric(Nage_LY[[i]][1,14:ncol(Nage_LY[[i]])]),paste(data_file),append = T,ncolumns = 45)
+ }
+T1<-noquote(paste0("# No Recruitments for ",FY, " to ",LY-RY))
+ write(T1,paste(data_file),append = T)
+ write((N_rec-RY),paste(data_file),append = T,ncolumns = 45)
+T1<-noquote("# Recruitment")
+ write(T1,paste(data_file),append = T)
+ write(round(Rec_1[1:(length(Rec_1)-RY)],1),paste(data_file),append = T,ncolumns = 45)
+T1<-noquote(paste("# SSB ", FY,"-",LY,sep=""))
+ write(T1,paste(data_file),append = T)
+ write(SSB,paste(data_file),append = T,ncolumns = 45)
+}
+
+
+mod1<-SS_output(dir)
+
+write_proj(data_file=paste0("data/",Model_Name,".dat"),data=mod1,
+ FY=1977,
+ LY=2020,
+ RY=2,
+ fleets=3,
+ sexes=1,
+ NAGES=10,
+ SSL=0,
+ Dorn = 0,
+ AUTHOR_F = 1,
+ SPR_ABC = 0.4,
+ SPR_MSY = 0.35,
+ SPAWN_M = 1,
+ FRATIO=noquote("0.3 0.2 0.5"))
+
+
+#' ## Save lists for running model to files expected by projection model
+library(gbm)
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,paste0("data/",Model_Name,"_spcat.dat"))
+
+file.copy(paste0("data/",Model_Name,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+
+#' ## Run projection model
+system("../src/main")
+#' ## Read in projection model mainfiles
+ .projdir= paste0(Model_Name,"/")
+ dir.create(.projdir)
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+ file.remove(list.files(getwd(), pattern="out$"))
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+#' ## Make plot of projection model simulations
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+
+
+#' ## Make tables
+ library(xtable)
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+
diff --git a/examples/Misc/admodel.cov b/examples/Misc/admodel.cov
new file mode 100644
index 0000000..410c598
Binary files /dev/null and b/examples/Misc/admodel.cov differ
diff --git a/examples/Misc/admodel.dep b/examples/Misc/admodel.dep
new file mode 100644
index 0000000..971e366
--- /dev/null
+++ b/examples/Misc/admodel.dep
@@ -0,0 +1,2 @@
+1 0
+Bzero 0
diff --git a/examples/Misc/admodel.hes b/examples/Misc/admodel.hes
new file mode 100644
index 0000000..154a7fc
Binary files /dev/null and b/examples/Misc/admodel.hes differ
diff --git a/examples/Misc/ai_re.prj b/examples/Misc/ai_re.prj
new file mode 100644
index 0000000..4894810
--- /dev/null
+++ b/examples/Misc/ai_re.prj
@@ -0,0 +1,35 @@
+Rougheye
+0 # SSL Species???
+0 # Constant buffer of Dorn
+1 # Number of fsheries
+1 # Number of sexes??
+# 0.0261738 # average 5 yr f
+# 0.0297296 # average F, 2022 and 2023
+# 0.0401244 # 2023 Fofl, start value
+# 0.028858427 # 2023 Fofl, second iteration
+ 0.031 # 2024 Fofl, start value
+# 0.0305124 # 2024 Fofl, second iteration
+1 # author f
+0.4 # ABC SPR
+0.35 # MSY SPR
+3 # Spawnmo
+43 # Number of ages
+1 # Fratio
+ # Natural mortality
+ 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307
+ # Maturity
+ 0.00340888 0.0044329 0.00576277 0.00748858 0.00972618 0.0126239 0.0163706 0.0212055 0.0274284 0.0354115 0.0456091 0.0585651 0.0749125 0.0953608 0.120663 0.151553 0.188655 0.232353 0.282646 0.339018 0.400358 0.464989 0.530817 0.59559 0.657196 0.713924 0.764628 0.808752 0.846267 0.877538 0.903175 0.923911 0.940499 0.953652 0.964008 0.972118 0.978442 0.983356 0.987164 0.99011 0.992385 0.99414 0.998196
+ # Wt Spawn
+ 55 81 112 147 188 232 280 331 384 440 497 556 615 675 735 795 854 913 970 1027 1082 1137 1189 1240 1290 1338 1384 1429 1472 1513 1553 1591 1627 1662 1695 1727 1757 1786 1814 1840 1865 1888 2006.7
+ # Wt Fish
+ 55 81 112 147 188 232 280 331 384 440 497 556 615 675 735 795 854 913 970 1027 1082 1137 1189 1240 1290 1338 1384 1429 1472 1513 1553 1591 1627 1662 1695 1727 1757 1786 1814 1840 1865 1888 2006.7
+ # selectivity
+ 0.000742953 0.00147408 0.00292258 0.0057862 0.0114235 0.0224293 0.0435709 0.0829494 0.152252 0.262859 0.414529 0.584339 0.736236 0.847145 0.916695 0.956234 0.977468 0.988524 0.994187 0.997064 0.998519 0.999254 0.999624 0.999811 0.999905 0.999952 0.999976 0.999988 0.999994 0.999997 0.999998 0.999999 1 1 1 1 1 1 1 1 1 1 1
+ # natage
+ 1.90469 1.8117 1.72321 1.48394 1.63 1.71119 1.65994 1.59983 1.32949 13.4142 1.05629 1.1426 1.22072 1.18872 0.931879 0.82203 0.933938 1.26836 0.821616 0.749373 0.705759 0.714874 0.449897 0.329574 0.251513 0.187913 0.14503 0.11613 0.0962753 0.0833661 0.0759574 0.0723329 0.0718631 0.0730712 0.0751903 0.0781308 0.0809843 0.0818001 0.0759746 0.0667667 0.0573137 0.0503982 0.792156
+ # Nrec
+35
+ # rec
+ 1.15837 1.14404 1.16416 1.21007 1.22741 1.17776 1.0395 0.895728 0.773063 0.677399 0.603982 0.554082 0.533258 0.539249 0.576319 0.645214 0.748979 0.902494 1.1235 1.36938 1.73894 2.57075 2.36166 2.3346 2.38527 3.43462 2.35989 1.93733 2.04665 2.43209 2.32951 2.04021 1.77262 21.2472 1.99402
+ # ssb
+ 4974.15 4968.42 4359.62 3715.78 3652.69 3626.45 3707.87 3827.03 3957.91 4098.91 4239.95 4368.17 4481.53 4444.2 4153.67 4180.48 3932.34 3747.46 3617.6 3551.1 3351.38 3128.31 3023.8 2959.06 2903.16 2781.58 2741.76 2727.58 2714.87 2725.5 2702.99 2691.78 2677.38 2661.55 2649.74
diff --git a/examples/Misc/ai_re_age10_m_22_1_original.dat b/examples/Misc/ai_re_age10_m_22_1_original.dat
new file mode 100644
index 0000000..4894810
--- /dev/null
+++ b/examples/Misc/ai_re_age10_m_22_1_original.dat
@@ -0,0 +1,35 @@
+Rougheye
+0 # SSL Species???
+0 # Constant buffer of Dorn
+1 # Number of fsheries
+1 # Number of sexes??
+# 0.0261738 # average 5 yr f
+# 0.0297296 # average F, 2022 and 2023
+# 0.0401244 # 2023 Fofl, start value
+# 0.028858427 # 2023 Fofl, second iteration
+ 0.031 # 2024 Fofl, start value
+# 0.0305124 # 2024 Fofl, second iteration
+1 # author f
+0.4 # ABC SPR
+0.35 # MSY SPR
+3 # Spawnmo
+43 # Number of ages
+1 # Fratio
+ # Natural mortality
+ 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307 0.0500307
+ # Maturity
+ 0.00340888 0.0044329 0.00576277 0.00748858 0.00972618 0.0126239 0.0163706 0.0212055 0.0274284 0.0354115 0.0456091 0.0585651 0.0749125 0.0953608 0.120663 0.151553 0.188655 0.232353 0.282646 0.339018 0.400358 0.464989 0.530817 0.59559 0.657196 0.713924 0.764628 0.808752 0.846267 0.877538 0.903175 0.923911 0.940499 0.953652 0.964008 0.972118 0.978442 0.983356 0.987164 0.99011 0.992385 0.99414 0.998196
+ # Wt Spawn
+ 55 81 112 147 188 232 280 331 384 440 497 556 615 675 735 795 854 913 970 1027 1082 1137 1189 1240 1290 1338 1384 1429 1472 1513 1553 1591 1627 1662 1695 1727 1757 1786 1814 1840 1865 1888 2006.7
+ # Wt Fish
+ 55 81 112 147 188 232 280 331 384 440 497 556 615 675 735 795 854 913 970 1027 1082 1137 1189 1240 1290 1338 1384 1429 1472 1513 1553 1591 1627 1662 1695 1727 1757 1786 1814 1840 1865 1888 2006.7
+ # selectivity
+ 0.000742953 0.00147408 0.00292258 0.0057862 0.0114235 0.0224293 0.0435709 0.0829494 0.152252 0.262859 0.414529 0.584339 0.736236 0.847145 0.916695 0.956234 0.977468 0.988524 0.994187 0.997064 0.998519 0.999254 0.999624 0.999811 0.999905 0.999952 0.999976 0.999988 0.999994 0.999997 0.999998 0.999999 1 1 1 1 1 1 1 1 1 1 1
+ # natage
+ 1.90469 1.8117 1.72321 1.48394 1.63 1.71119 1.65994 1.59983 1.32949 13.4142 1.05629 1.1426 1.22072 1.18872 0.931879 0.82203 0.933938 1.26836 0.821616 0.749373 0.705759 0.714874 0.449897 0.329574 0.251513 0.187913 0.14503 0.11613 0.0962753 0.0833661 0.0759574 0.0723329 0.0718631 0.0730712 0.0751903 0.0781308 0.0809843 0.0818001 0.0759746 0.0667667 0.0573137 0.0503982 0.792156
+ # Nrec
+35
+ # rec
+ 1.15837 1.14404 1.16416 1.21007 1.22741 1.17776 1.0395 0.895728 0.773063 0.677399 0.603982 0.554082 0.533258 0.539249 0.576319 0.645214 0.748979 0.902494 1.1235 1.36938 1.73894 2.57075 2.36166 2.3346 2.38527 3.43462 2.35989 1.93733 2.04665 2.43209 2.32951 2.04021 1.77262 21.2472 1.99402
+ # ssb
+ 4974.15 4968.42 4359.62 3715.78 3652.69 3626.45 3707.87 3827.03 3957.91 4098.91 4239.95 4368.17 4481.53 4444.2 4153.67 4180.48 3932.34 3747.46 3617.6 3551.1 3351.38 3128.31 3023.8 2959.06 2903.16 2781.58 2741.76 2727.58 2714.87 2725.5 2702.99 2691.78 2677.38 2661.55 2649.74
diff --git a/examples/Misc/ai_re_age10_m_22_1_original_spcat.dat b/examples/Misc/ai_re_age10_m_22_1_original_spcat.dat
new file mode 100644
index 0000000..25b398d
--- /dev/null
+++ b/examples/Misc/ai_re_age10_m_22_1_original_spcat.dat
@@ -0,0 +1,37 @@
+ #_Number_of_years with specified catch (if begin-yr = 2004, and this number is "3", then subsequent values represent catches in 2004, 05, and 06 (to evaluate alts for 2007)
+ 2
+ # Number of species
+ 1
+ # OY Minimum
+.1343248 # Note that this is for age-structured species 1330.148
+ # OY Maximum
+1943248 # Note that this is for age-structured species 1930.148
+ # data files for each species
+ # 1
+ data\ai_re_age10_m_22_1_original.dat
+ # ABC Multipliers
+ 1
+ # scalars
+ 1
+ # New Alt 4 Fabc SPRs (Rockfish = 0.75, other 0.6), Steller sea lion prey species between F40 and F60 (to meet OY Min)
+0.75
+ # Number of TAC model categories
+ 1
+ # TAC model indices (for aggregating)
+ 1
+ # Catch in each future year
+ 2022 305
+ 2023 395
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/examples/Misc/ai_re_sppcat.dat b/examples/Misc/ai_re_sppcat.dat
new file mode 100644
index 0000000..08ab02a
--- /dev/null
+++ b/examples/Misc/ai_re_sppcat.dat
@@ -0,0 +1,55 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#nyrs_fixed_catch
+2
+#nspp
+1
+#OY_min
+.1343248 # Note that this is for age-structured species 1330.148
+#OY_max
+1943248 # Note that this is for age-structured species 1930.148
+# data files for each species
+data\ai_re_age10_m_22_1_original.dat
+#ABC_mults
+1
+# scalars
+1
+# New Alt 4 Fabc SPRs (Rockfish = 0.75, other 0.6), Steller sea lion prey species between F40 and F60 (to meet OY Min)
+0.75
+# Number of TAC model categories
+1
+# TAC model indices (for aggregating)
+1
+# Catch in each future year
+ 2022 305
+ 2023 395
diff --git a/examples/Misc/ai_spm.dat b/examples/Misc/ai_spm.dat
new file mode 100644
index 0000000..9876642
--- /dev/null
+++ b/examples/Misc/ai_spm.dat
@@ -0,0 +1,55 @@
+#rn
+AI_re
+#Tier
+3
+#nalts
+7
+#alts
+1
+2
+3
+4
+5
+6
+7
+#tac_flag
+1
+#srr_type
+1
+#srr_form
+1
+#srr_conditioning
+0
+#srr_reserved
+0
+#spm_detail_flag
+1
+#nprj_yrs
+100
+#nsims
+1000
+#beg_yr
+2022
+#nyrs_fixed_catch
+2
+#nspp
+1
+#OY_min
+10
+#OY_max
+1943248 # Note that this is for age-structured species 1930.148
+#data_files
+ai_re.prj
+#ABC_mults
+1
+#scalars
+1
+#alt4_spr
+0.6
+#nTAC_cat
+1
+#TACind
+1
+#fixed_catch
+ 2022 305
+ 2023 395
diff --git a/examples/Misc/akp.R b/examples/Misc/akp.R
new file mode 100644
index 0000000..58e8ca7
--- /dev/null
+++ b/examples/Misc/akp.R
@@ -0,0 +1,127 @@
+#' # Projection model work
+#' Projection model work for demonstrating application of controls and input data
+
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+source("../R/readData.R")
+
+#' ## Set initial "setup" parameters
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC (1 means true, otherwise false)
+ srr = 1 , #' Stock-recruitment type (1=Ricker, 2=Bholt)
+ rec_proj = 1, #' projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+ srr_cond = 0 , #' SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 100, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+
+#' ## Set up the species specific run file
+config<-list(
+ nFixCatchYrs = 2,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote("data/t1.dat"),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c( 2016,55000., 2017,55000. )
+)
+
+#' ## Make list of main file w/ assessment model results
+#' E.g., "data/bsai_atka.dat"
+datfile <- list(
+ runname = noquote("M16.2"),
+ ssl_spp = 1, # SSL_spp
+ Dorn_buffer = 1, # Dorn_buffer
+ nfsh = 1, # N_fsh
+ nsex = 1, # N_sexes
+ avgF5yr = 0.0661399, # avg_5yr_F
+ F40_mult = 1, # F_40_multiplier
+ spr_abc = 0.4, # SPR_abc
+ spr_msy = 0.35, # SPR_msy
+ sp_mo = 8, # spawn_month
+ nages = 11, # N_ages
+ Frat = 1, # F_ratio
+ # M
+ M = c(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
+ # Maturity
+ pmat = c(0.005,0.037,0.224,0.688,0.944,0.992,0.999,1,1,1,1),
+ # Wt_at_age_spawners
+ wtage_sp = c(44.8,161.377,398.272,557.695,652.113,719.573,863.744,948.744,921.397,885.912,1069.87),
+ # Wt_at_age_fsh
+ wtage_fsh = c(69.3778,253.522,408.211,614.731,668.483,718.137,803.017,798.707,788.117,842.468,960.006),
+ # select
+ sel = c(0.002576427,0.040030753,0.651104228,0.768404263,0.794886081,1,0.889293108,0.604815671,0.451169778,0.403195516,0.403195516),
+ # N
+ N = c(511.179,378.528,278.443,194.385,183.423,45.5404,61.1188,19.8073,39.6285,33.6501,51.7806),
+ # Nyrs
+ nyrs = 37,
+ # recruits
+ R = c(1578.51,479.509,357.919,443.588,318.981,413.125,514.351,600.987,536.301,692.34,452.279,1618.73,702.801,372.811,597.812,1136.19,402.862,424.179,1025.36,207.05,383.695,1054.64,2224.52,1379.34,1545.81,345.556,454.884,617.194,404.876,993.497,727.754,236.795,505.948,258.653,726.627,524.198,473.54),
+ # SSB
+ SSB = c(206.391,194.569,187.097,183.296,195.289,240.774,252.143,238.163,223.199,200.776,182.378,179.671,189.193,198.242,212.123,233.484,277.891,282.004,250.709,231.816,218.403,195.275,181.628,189.976,175.912,168.624,220.206,315.124,376.621,397.171,365.476,317.159,277.777,242.719,233.415,223.636,198.117,183.537,177.91)
+ )
+
+#' ## Save lists for running model to files expected by projection model
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,"data/t1_spcat.dat")
+runfn<-"t1"
+file.copy(paste0("data/",runfn,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+
+list2dat(datfile,"data/t1.dat")
+
+#' ## Run projection model
+system("../src/main")
+#' ## Read in projection model mainfiles
+ .projdir="akp_out/"
+ dir.create(.projdir)
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+ file.remove(list.files(getwd(), pattern="out$"))
+
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+#' ## Make plot of projection model simulations
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+ library(patchwork)
+ p1/p2
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+
+
+#' ## Make tables
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+
diff --git a/examples/Misc/ebswp_spcat.dat b/examples/Misc/ebswp_spcat.dat
new file mode 100755
index 0000000..4bfc68e
--- /dev/null
+++ b/examples/Misc/ebswp_spcat.dat
@@ -0,0 +1,30 @@
+#_SETUP_FILE_FOR_THE_BERING_SEA/AI_FISHERIES
+#_Number_of_years with specified catch (if begin-yr = 2005, and this number is "3", then subsequent values represent catches in 2005, 06, and 07 (to evaluate alts for 2008)
+3
+# Number of species
+1
+# OY Minimum
+0.000000 # Note that this is for age-structured species 1330.148
+# OY Maximum
+3500.000 # Note that this is for age-structured species 1930.148
+../pm.prj
+# ABC Multipliers
+1
+# Population scalars
+1
+# New Alt 4 Fabc SPRs (Rockfish = 0.75, other 0.6), Steller sea lion prey species between F40 and F60 (to meet OY Min)
+0.6
+# Number of TAC model categories
+1
+# TAC model indices (for aggregating)
+1
+# Catch in each future year
+2019 1390
+2020 1350
+2021 1324
+2022 1341
+2023 1354
+2024 1365
+2025 1375
+2026 1377
+2017 1350
diff --git a/examples/Misc/elasticity.csv b/examples/Misc/elasticity.csv
new file mode 100644
index 0000000..e69de29
diff --git a/examples/Misc/example.R b/examples/Misc/example.R
new file mode 100644
index 0000000..b57a9c2
--- /dev/null
+++ b/examples/Misc/example.R
@@ -0,0 +1,124 @@
+#' # Projection model work
+#' Projection model work for demonstrating application of controls and input data
+
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+source("../R/readData.R")
+
+#' ## Set initial "setup" parameters
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC 1 means true, otherwise false
+ srr = 1 , #' Stock-recruitment type 1=Ricker, 2=Bholt
+ rec_proj = 1, #' projection rec form default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R
+ srr_cond = 0 , #' SR-Conditioning 0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition affects SRR fits
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file of all simulations rather than a summary, 0 means don't do it, otherwise do it Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 100, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+
+#' ## Set up the species specific run file
+config<-list(
+ nFixCatchYrs = 2,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote("data/t1.dat"),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c( 2016,55000., 2017,55000. )
+)
+
+#' ## Make list of main file w/ assessment model results
+#' E.g., "data/bsai_atka.dat"
+datfile <- list(
+ runname = noquote("M16.2"),
+ ssl_spp = 1, # SSL_spp
+ Dorn_buffer = 1, # Dorn_buffer
+ nfsh = 1, # N_fsh
+ nsex = 1, # N_sexes
+ avgF5yr = 0.0661399, # avg_5yr_F
+ F40_mult = 1, # F_40_multiplier
+ spr_abc = 0.4, # SPR_abc
+ spr_msy = 0.35, # SPR_msy
+ sp_mo = 8, # spawn_month
+ nages = 11, # N_ages
+ Frat = 1, # F_ratio
+ # M
+ M = c(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
+ # Maturity
+ pmat = c(0.005,0.037,0.224,0.688,0.944,0.992,0.999,1,1,1,1),
+ # Wt_at_age_spawners
+ wtage_sp = c(44.8,161.377,398.272,557.695,652.113,719.573,863.744,948.744,921.397,885.912,1069.87),
+ # Wt_at_age_fsh
+ wtage_fsh = c(69.3778,253.522,408.211,614.731,668.483,718.137,803.017,798.707,788.117,842.468,960.006),
+ # select
+ sel = c(0.002576427,0.040030753,0.651104228,0.768404263,0.794886081,1,0.889293108,0.604815671,0.451169778,0.403195516,0.403195516),
+ # N
+ N = c(511.179,378.528,278.443,194.385,183.423,45.5404,61.1188,19.8073,39.6285,33.6501,51.7806),
+ # Nyrs
+ nyrs = 37,
+ # recruits
+ R = c(1578.51,479.509,357.919,443.588,318.981,413.125,514.351,600.987,536.301,692.34,452.279,1618.73,702.801,372.811,597.812,1136.19,402.862,424.179,1025.36,207.05,383.695,1054.64,2224.52,1379.34,1545.81,345.556,454.884,617.194,404.876,993.497,727.754,236.795,505.948,258.653,726.627,524.198,473.54),
+ # SSB
+ SSB = c(206.391,194.569,187.097,183.296,195.289,240.774,252.143,238.163,223.199,200.776,182.378,179.671,189.193,198.242,212.123,233.484,277.891,282.004,250.709,231.816,218.403,195.275,181.628,189.976,175.912,168.624,220.206,315.124,376.621,397.171,365.476,317.159,277.777,242.719,233.415,223.636,198.117,183.537,177.91)
+ )
+
+#' ## Save lists for running model to files expected by projection model
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,"data/t1_spcat.dat")
+runfn<-"t1"
+file.copy(paste0("data/",runfn,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+
+list2dat(datfile,"data/t1.dat")
+
+#' ## Run projection model
+system("../src/main")
+#' ## Read in projection model mainfiles
+ .projdir="t1/"
+ dir.create(.projdir)
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+ file.remove(list.files(getwd(), pattern="out$"))
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+#' ## Make plot of projection model simulations
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+
+
+#' ## Make tables
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+
diff --git a/examples/Misc/goanrs_all_spcat.dat b/examples/Misc/goanrs_all_spcat.dat
new file mode 100644
index 0000000..c3c9e3e
--- /dev/null
+++ b/examples/Misc/goanrs_all_spcat.dat
@@ -0,0 +1,28 @@
+#_SETUP_FILE_FOR_THE_BERING_SEA/AI_FISHERIES
+#_Number_of_years with specified catch (if begin-yr = 2005, and this number is 3, then subsequent values represent catches in 2005, 06, and 7 (to evaluate alts for 2008
+3
+# Number of runs
+3
+# OY Minimum
+13.248 # Note that this is for age-structured species 1330.148
+# OY Maximum
+1943.248 # Note that this is for age-structured species 1930.148
+# data files for each species
+# Pollock Pacific cod sablefish Yellowfin Greenland turbot arrowtooth flounder Rock sole Flathead sole AK Plaice Pacific ocean perch Nrthrn RF Atka mackerel
+# 1 2 3 4 5 6 7 8 9 10 11 12
+# data files for each species
+data/goanrs.dat data/goanrs_western.dat data/goanrs_central.dat
+# ABC Multipliers
+1 1 1
+# Population scalars
+1 1 1
+# New Alt 4 Fabc SPRs (Rockfish = 0.75, other 0.6), Steller sea lion prey species between F40 and F60 (to meet OY Min)
+0.75 .75 .75
+# Number of TAC model categories
+1 1 1
+# TAC model indices (for aggregating)
+1 1 1
+# Catch in each future year
+2021 683.728 666 333
+2022 683.728 666 333
+2023 683.728 666 333
diff --git a/examples/Misc/goars.R b/examples/Misc/goars.R
new file mode 100644
index 0000000..8ca1232
--- /dev/null
+++ b/examples/Misc/goars.R
@@ -0,0 +1,54 @@
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+system("run.sh goanrs_central")
+system("run.sh goanrs_western")
+system("run.sh goanrs")
+
+# Or do all three above in one run
+system("run.sh goanrs_all")
+
+
+plot_ssb <- function(dirname,title=NULL){
+ df <- read_table(paste0(dirname,"_out/bigfile.out"))
+ if(is.null(title)) title=dirname
+ df %>% mutate(Alt=as.factor(Alt)) %>% group_by(Yr,Alt) %>%
+ summarise(SSB=mean(SSB),Catch=mean(Catch),ABC=mean(ABC)) %>%
+ ggplot(aes(x=Yr,y=SSB,color=Alt)) + geom_line(size=2) +
+ expand_limits(y=0) + theme_few() + ggtitle(title)
+}
+plot_ssb(dirname="goanrs_central","N rock sole, Central")
+plot_ssb(dirname="goanrs_western","N rock sole, western")
+
+plot_ssb <- function(dirname,alt=c(1,7),title=NULL){
+ df <- read_table(paste0(dirname,"_out/bigfile.out"))
+ if(is.null(title)) title=dirname
+ df %>% filter(Alt %in% alt) %>% group_by(Yr,Stock,Alt) %>%
+ summarise(ub=quantile(SSB,.9),lb=quantile(SSB,.1),
+ SSB=mean(SSB), Catch=mean(Catch),ABC=mean(ABC)) %>%
+ mutate(Stock=as.factor(Stock),Alt=as.factor(Alt)) %>%
+ ggplot(aes(x=Yr,y=SSB,ymax=ub,ymin=lb,color=Stock:Alt,fill=Stock:Alt)) + geom_line(size=2) +
+ geom_ribbon(alpha=.2) +
+ expand_limits(y=0) + theme_few() + ggtitle(title)
+}
+names(df)
+dirname <- "goanrs_all"
+dirname="goanrs_all";title="N rock sole, GOA";alt=c(1,6)
+plot_ssb(dirname="goanrs_all",title="N rock sole, GOA")
+
+plot_catch <- function(dirname,alt=c(1,7),title=NULL){
+ df <- read_table(paste0(dirname,"_out/bigfile.out"))
+ if(is.null(title)) title=dirname
+ df %>% filter(Alt %in% alt) %>% group_by(Yr,Stock,Alt) %>%
+ summarise(ub=quantile(Catch,.9),lb=quantile(Catch,.1),
+ SSB=mean(SSB), Catch=mean(Catch),ABC=mean(ABC)) %>%
+ mutate(Stock=as.factor(Stock),Alt=as.factor(Alt)) %>%
+ ggplot(aes(x=Yr,y=Catch,ymax=ub,ymin=lb,color=Stock:Alt,fill=Stock:Alt)) + geom_line(size=2) +
+ geom_ribbon(alpha=.2) +
+ expand_limits(y=0) + theme_few() + ggtitle(title)
+}
+names(df)
+plot_catch(dirname="goanrs_all",alt=1,title="N rock sole, GOA")
+plot_ssb(dirname="goanrs_all",alt=c(7,6),title="N rock sole, GOA")
+plot_catch(dirname="goanrs_all",alt=c(7,6),title="N rock sole, GOA")
+plot_ssb(dirname="goanrs_western","lN rock sole, western")
\ No newline at end of file
diff --git a/examples/Misc/overlap.csv b/examples/Misc/overlap.csv
new file mode 100644
index 0000000..cdba7ee
--- /dev/null
+++ b/examples/Misc/overlap.csv
@@ -0,0 +1,21 @@
+Years,Index,Type
+2000,-1.032574220232773,Overlap
+2001,0.95787556661264,Overlap
+2002,0.5853330121430638,Overlap
+2003,1.1939635815403216,Overlap
+2004,0.7308977125921574,Overlap
+2005,1.4097958669161312,Overlap
+2006,-0.14408474281886932,Overlap
+2007,-2.025866807486055,Overlap
+2008,-0.37906514010546777,Overlap
+2009,-1.1128692393332873,Overlap
+2010,-1.1523477685650252,Overlap
+2011,-0.24875402841591637,Overlap
+2012,-0.39238679883000743,Overlap
+2013,-1.22209318810684,Overlap
+2014,0.9169597109465762,Overlap
+2015,0.6287713981905791,Overlap
+2016,1.1721382815301298,Overlap
+2017,-0.7673162903327927,Overlap
+2018,-0.10442071014339757,Overlap
+2019,0.9860438038988355,Overlap
diff --git a/examples/Misc/projtimeline.xlsx b/examples/Misc/projtimeline.xlsx
new file mode 100644
index 0000000..3f1a8dc
Binary files /dev/null and b/examples/Misc/projtimeline.xlsx differ
diff --git a/examples/Misc/run.bat b/examples/Misc/run.bat
new file mode 100644
index 0000000..4066a28
--- /dev/null
+++ b/examples/Misc/run.bat
@@ -0,0 +1,20 @@
+:: copy data\%1_setup.dat setup.dat
+:: copy data\%1_tacpar.dat tacpar.dat
+@echo off
+
+if NOT exist data\%1_spcat.dat ( echo ---------------
+echo.
+echo Oops...Error!!!!
+echo.
+echo File data\%1_spcat.dat appears to be missing...
+echo.
+echo Exiting....
+echo ---------------
+exit /B
+)
+copy data\%1_spcat.dat spp_catch.dat
+..\src\main -nox -nohess
+if NOT exist %1_out mkdir %1_out
+copy *.out %1_out\
+copy main.rep %1_out\report.out
+copy alt2_proj.rep %1_out\alt2_proj.out
diff --git a/examples/Misc/run.sh b/examples/Misc/run.sh
new file mode 100755
index 0000000..e86afe0
--- /dev/null
+++ b/examples/Misc/run.sh
@@ -0,0 +1,10 @@
+test -f data/$1_spcat.dat || echo "Missing file..." data/$1
+test -f data/$1_spcat.dat || exit
+echo $1
+\cp data/$1_spcat.dat spp_catch.dat
+spm -nohess
+test -d $1_out || mkdir $1_out
+\mv *.out $1_out
+rm eigv.rpt variance admodel.* *.r0? *.p0? fmin.log *.b0?
+
+
diff --git a/examples/Misc/setup1.dat b/examples/Misc/setup1.dat
new file mode 100644
index 0000000..2eba790
--- /dev/null
+++ b/examples/Misc/setup1.dat
@@ -0,0 +1,19 @@
+std
+3
+7
+1
+2
+3
+4
+5
+6
+7
+1
+1
+1
+1
+0
+1
+14
+1000
+2015
diff --git a/examples/Misc/spmR_cases.xlsx b/examples/Misc/spmR_cases.xlsx
new file mode 100644
index 0000000..d257732
Binary files /dev/null and b/examples/Misc/spmR_cases.xlsx differ
diff --git a/examples/Misc/spp_catch.dat b/examples/Misc/spp_catch.dat
new file mode 100644
index 0000000..63b161f
--- /dev/null
+++ b/examples/Misc/spp_catch.dat
@@ -0,0 +1,26 @@
+# a new file
+#nFixCatchYrs
+2
+#nSpecies
+1
+#OYMin
+0.1343248
+#OYMax
+1943248
+#dataFiles
+data/t1.dat
+#ABCMult
+1
+#PoplnScalar
+1000
+#AltFabcSPR
+0.75
+#nTAC
+1
+#TACIndices
+1
+#Catch
+2016
+55000
+2017
+55000
diff --git a/examples/Misc/srecpar.dat b/examples/Misc/srecpar.dat
new file mode 100644
index 0000000..cf22bcb
--- /dev/null
+++ b/examples/Misc/srecpar.dat
@@ -0,0 +1,4 @@
+# Bzero, PhiZero, Alpha, sigmaR
+675136 0.00306701 2.34603 0.6
+# Bzero, PhiZero, Alpha, sigmaR
+675136 0.00306701 2.34603 0.6
diff --git a/examples/Misc/tacpar.dat b/examples/Misc/tacpar.dat
new file mode 100644
index 0000000..6631707
--- /dev/null
+++ b/examples/Misc/tacpar.dat
@@ -0,0 +1,10 @@
+7
+6
+ 2.60197 0.417 0.372 0.361 0.296 0.2733 0.125
+ -0.300856 -0.741664 -1.26797 -1.78614 -2.15198 -2.39595 -2.59939
+ -0.00703968 -0.956592 -1.48909 -2.09122 -2.26005 -2.31812 -2.45321
+ 0.372276 -1.35066 -1.76076 -2.22163 -2.31726 -2.28107 -2.34758
+ 0.505535 -1.19926 -1.68812 -2.39543 -2.40752 -2.5904 -2.58726
+ 0.368722 -1.59025 -1.79548 -2.67422 -2.61358 -2.41092 -2.70204
+ 0.308248 -1.78052 -2.23264 -2.85605 -3.39375 -3.05209 -2.94219
+ 0.180676 -1.7586 -1.80222 -3.09402 -2.2618 -3.43215 -3.06417
diff --git a/examples/Misc/test.R b/examples/Misc/test.R
new file mode 100644
index 0000000..aec7b47
--- /dev/null
+++ b/examples/Misc/test.R
@@ -0,0 +1,124 @@
+#' # Projection model work
+#' Projection model work for demonstrating application of controls and input data
+
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+source("../R/readData.R")
+
+#' ## Set initial "setup" parameters
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC (1 means true, otherwise false)
+ srr = 1 , #' Stock-recruitment type (1=Ricker, 2=Bholt)
+ rec_proj = 1, #' projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+ srr_cond = 0 , #' SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 100, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+
+#' ## Set up the species specific run file
+config<-list(
+ nFixCatchYrs = 2,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote("data/t1.dat"),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c( 2016,55000., 2017,55000. )
+)
+
+#' ## Make list of main file w/ assessment model results
+#' E.g., "data/bsai_atka.dat"
+datfile <- list(
+ runname = noquote("M16.2"),
+ ssl_spp = 1, # SSL_spp
+ Dorn_buffer = 1, # Dorn_buffer
+ nfsh = 1, # N_fsh
+ nsex = 1, # N_sexes
+ avgF5yr = 0.0661399, # avg_5yr_F
+ F40_mult = 1, # F_40_multiplier
+ spr_abc = 0.4, # SPR_abc
+ spr_msy = 0.35, # SPR_msy
+ sp_mo = 8, # spawn_month
+ nages = 11, # N_ages
+ Frat = 1, # F_ratio
+ # M
+ M = c(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
+ # Maturity
+ pmat = c(0.005,0.037,0.224,0.688,0.944,0.992,0.999,1,1,1,1),
+ # Wt_at_age_spawners
+ wtage_sp = c(44.8,161.377,398.272,557.695,652.113,719.573,863.744,948.744,921.397,885.912,1069.87),
+ # Wt_at_age_fsh
+ wtage_fsh = c(69.3778,253.522,408.211,614.731,668.483,718.137,803.017,798.707,788.117,842.468,960.006),
+ # select
+ sel = c(0.002576427,0.040030753,0.651104228,0.768404263,0.794886081,1,0.889293108,0.604815671,0.451169778,0.403195516,0.403195516),
+ # N
+ N = c(511.179,378.528,278.443,194.385,183.423,45.5404,61.1188,19.8073,39.6285,33.6501,51.7806),
+ # Nyrs
+ nyrs = 37,
+ # recruits
+ R = c(1578.51,479.509,357.919,443.588,318.981,413.125,514.351,600.987,536.301,692.34,452.279,1618.73,702.801,372.811,597.812,1136.19,402.862,424.179,1025.36,207.05,383.695,1054.64,2224.52,1379.34,1545.81,345.556,454.884,617.194,404.876,993.497,727.754,236.795,505.948,258.653,726.627,524.198,473.54),
+ # SSB
+ SSB = c(206.391,194.569,187.097,183.296,195.289,240.774,252.143,238.163,223.199,200.776,182.378,179.671,189.193,198.242,212.123,233.484,277.891,282.004,250.709,231.816,218.403,195.275,181.628,189.976,175.912,168.624,220.206,315.124,376.621,397.171,365.476,317.159,277.777,242.719,233.415,223.636,198.117,183.537,177.91)
+ )
+
+#' ## Save lists for running model to files expected by projection model
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,"data/t1_spcat.dat")
+runfn<-"t1"
+file.copy(paste0("data/",runfn,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+
+list2dat(datfile,"data/t1.dat")
+
+#' ## Run projection model
+system("../src/main")
+#' ## Read in projection model mainfiles
+ .projdir="t1/"
+ dir.create(.projdir)
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+ file.remove(list.files(getwd(), pattern="out$"))
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+#' ## Make plot of projection model simulations
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+
+
+#' ## Make tables
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+
diff --git a/examples/Misc/test.Rmd b/examples/Misc/test.Rmd
new file mode 100644
index 0000000..12e5e92
--- /dev/null
+++ b/examples/Misc/test.Rmd
@@ -0,0 +1,151 @@
+# Projection model work
+Projection model work for demonstrating application of controls and input data
+
+```{r }
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+source("../R/readData.R")
+```
+
+## Set initial "setup" parameters
+
+```{r }
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC (1 means true, otherwise false)
+ srr = 1 , #' Stock-recruitment type (1=Ricker, 2=Bholt)
+ rec_proj = 1, #' projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+ srr_cond = 0 , #' SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 100, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+```
+
+## Set up the species specific run file
+
+```{r }
+config<-list(
+ nFixCatchYrs = 2,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote("data/t1.dat"),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c( 2016,55000., 2017,55000. )
+)
+```
+
+## Make list of main file w/ assessment model results
+E.g., "data/bsai_atka.dat"
+
+```{r }
+datfile <- list(
+ runname = noquote("M16.2"),
+ ssl_spp = 1, # SSL_spp
+ Dorn_buffer = 1, # Dorn_buffer
+ nfsh = 1, # N_fsh
+ nsex = 1, # N_sexes
+ avgF5yr = 0.0661399, # avg_5yr_F
+ F40_mult = 1, # F_40_multiplier
+ spr_abc = 0.4, # SPR_abc
+ spr_msy = 0.35, # SPR_msy
+ sp_mo = 8, # spawn_month
+ nages = 11, # N_ages
+ Frat = 1, # F_ratio
+ # M
+ M = c(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
+ # Maturity
+ pmat = c(0.005,0.037,0.224,0.688,0.944,0.992,0.999,1,1,1,1),
+ # Wt_at_age_spawners
+ wtage_sp = c(44.8,161.377,398.272,557.695,652.113,719.573,863.744,948.744,921.397,885.912,1069.87),
+ # Wt_at_age_fsh
+ wtage_fsh = c(69.3778,253.522,408.211,614.731,668.483,718.137,803.017,798.707,788.117,842.468,960.006),
+ # select
+ sel = c(0.002576427,0.040030753,0.651104228,0.768404263,0.794886081,1,0.889293108,0.604815671,0.451169778,0.403195516,0.403195516),
+ # N
+ N = c(511.179,378.528,278.443,194.385,183.423,45.5404,61.1188,19.8073,39.6285,33.6501,51.7806),
+ # Nyrs
+ nyrs = 37,
+ # recruits
+ R = c(1578.51,479.509,357.919,443.588,318.981,413.125,514.351,600.987,536.301,692.34,452.279,1618.73,702.801,372.811,597.812,1136.19,402.862,424.179,1025.36,207.05,383.695,1054.64,2224.52,1379.34,1545.81,345.556,454.884,617.194,404.876,993.497,727.754,236.795,505.948,258.653,726.627,524.198,473.54),
+ # SSB
+ SSB = c(206.391,194.569,187.097,183.296,195.289,240.774,252.143,238.163,223.199,200.776,182.378,179.671,189.193,198.242,212.123,233.484,277.891,282.004,250.709,231.816,218.403,195.275,181.628,189.976,175.912,168.624,220.206,315.124,376.621,397.171,365.476,317.159,277.777,242.719,233.415,223.636,198.117,183.537,177.91)
+ )
+```
+
+## Save lists for running model to files expected by projection model
+
+```{r }
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,"data/t1_spcat.dat")
+runfn<-"t1"
+file.copy(paste0("data/",runfn,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+
+list2dat(datfile,"data/t1.dat")
+```
+
+## Run projection model
+
+```{r }
+system("../src/main")
+```
+
+## Read in projection model mainfiles
+
+```{r }
+ .projdir="t1/"
+ dir.create(.projdir)
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+ file.remove(list.files(getwd(), pattern="out$"))
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+```
+
+## Make plot of projection model simulations
+
+```{r }
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+```
+
+## Make tables
+
+```{r }
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+```
+
diff --git a/examples/Misc/test.html b/examples/Misc/test.html
new file mode 100644
index 0000000..48a1625
--- /dev/null
+++ b/examples/Misc/test.html
@@ -0,0 +1,455 @@
+
+
+
+
+
+Projection model work
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Projection model work
+
+Projection model work for demonstrating application of controls and input data
+
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+source("../R/readData.R")
+
+
+Set initial “setup” parameters
+
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC (1 means true, otherwise false)
+ srr = 1 , #' Stock-recruitment type (1=Ricker, 2=Bholt)
+ rec_proj = 1, #' projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+ srr_cond = 0 , #' SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 100, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+
+
+Set up the species specific run file
+
+config<-list(
+ nFixCatchYrs = 2,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote("data/t1.dat"),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c( 2016,55000., 2017,55000. )
+)
+
+
+Make list of main file w/ assessment model results
+
+E.g., “data/bsai_atka.dat”
+
+datfile <- list(
+ runname = noquote("M16.2"),
+ ssl_spp = 1, # SSL_spp
+ Dorn_buffer = 1, # Dorn_buffer
+ nfsh = 1, # N_fsh
+ nsex = 1, # N_sexes
+ avgF5yr = 0.0661399, # avg_5yr_F
+ F40_mult = 1, # F_40_multiplier
+ spr_abc = 0.4, # SPR_abc
+ spr_msy = 0.35, # SPR_msy
+ sp_mo = 8, # spawn_month
+ nages = 11, # N_ages
+ Frat = 1, # F_ratio
+ # M
+ M = c(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
+ # Maturity
+ pmat = c(0.005,0.037,0.224,0.688,0.944,0.992,0.999,1,1,1,1),
+ # Wt_at_age_spawners
+ wtage_sp = c(44.8,161.377,398.272,557.695,652.113,719.573,863.744,948.744,921.397,885.912,1069.87),
+ # Wt_at_age_fsh
+ wtage_fsh = c(69.3778,253.522,408.211,614.731,668.483,718.137,803.017,798.707,788.117,842.468,960.006),
+ # select
+ sel = c(0.002576427,0.040030753,0.651104228,0.768404263,0.794886081,1,0.889293108,0.604815671,0.451169778,0.403195516,0.403195516),
+ # N
+ N = c(511.179,378.528,278.443,194.385,183.423,45.5404,61.1188,19.8073,39.6285,33.6501,51.7806),
+ # Nyrs
+ nyrs = 37,
+ # recruits
+ R = c(1578.51,479.509,357.919,443.588,318.981,413.125,514.351,600.987,536.301,692.34,452.279,1618.73,702.801,372.811,597.812,1136.19,402.862,424.179,1025.36,207.05,383.695,1054.64,2224.52,1379.34,1545.81,345.556,454.884,617.194,404.876,993.497,727.754,236.795,505.948,258.653,726.627,524.198,473.54),
+ # SSB
+ SSB = c(206.391,194.569,187.097,183.296,195.289,240.774,252.143,238.163,223.199,200.776,182.378,179.671,189.193,198.242,212.123,233.484,277.891,282.004,250.709,231.816,218.403,195.275,181.628,189.976,175.912,168.624,220.206,315.124,376.621,397.171,365.476,317.159,277.777,242.719,233.415,223.636,198.117,183.537,177.91)
+ )
+
+
+Save lists for running model to files expected by projection model
+
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,"data/t1_spcat.dat")
+runfn<-"t1"
+file.copy(paste0("data/",runfn,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+
+
+## [1] TRUE
+
+
+list2dat(datfile,"data/t1.dat")
+
+
+Run projection model
+
+system("../src/main")
+
+
+Read in projection model mainfiles
+
+ .projdir="t1/"
+ dir.create(.projdir)
+
+
+## Warning in dir.create(.projdir): 't1' already exists
+
+
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+
+
+## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
+
+
+ file.remove(list.files(getwd(), pattern="out$"))
+
+
+## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
+
+
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+
+
+Make plot of projection model simulations
+
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+
+
+
+
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+
+
+Make tables
+
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+
+
+## <!-- html table generated in R 4.0.0 by xtable 1.8-4 package -->
+## <!-- Tue Jun 9 08:53:57 2020 -->
+## <table border=1>
+## <caption align="top"> Tier 3 projections of BSAI Atka mackerel catch for the 7 scenarios. </caption>
+## <tr> <th> Catch </th> <th> Scenario.1 </th> <th> Scenario.2 </th> <th> Scenario.3 </th> <th> Scenario.4 </th> <th> Scenario.5 </th> <th> Scenario.6 </th> <th> Scenario.7 </th> </tr>
+## <tr> <td align="right"> 2019 </td> <td align="right"> 55 </td> <td align="right"> 55 </td> <td align="right"> 55 </td> <td align="right"> 55 </td> <td align="right"> 55 </td> <td align="right"> 55 </td> <td align="right"> 55 </td> </tr>
+## <tr> <td align="right"> 2020 </td> <td align="right"> 91 </td> <td align="right"> 55 </td> <td align="right"> 21 </td> <td align="right"> 27 </td> <td align="right"> 0 </td> <td align="right"> 106 </td> <td align="right"> 91 </td> </tr>
+## <tr> <td align="right"> 2021 </td> <td align="right"> 83 </td> <td align="right"> 85 </td> <td align="right"> 22 </td> <td align="right"> 28 </td> <td align="right"> 0 </td> <td align="right"> 91 </td> <td align="right"> 83 </td> </tr>
+## <tr> <td align="right"> 2022 </td> <td align="right"> 81 </td> <td align="right"> 80 </td> <td align="right"> 25 </td> <td align="right"> 30 </td> <td align="right"> 0 </td> <td align="right"> 85 </td> <td align="right"> 93 </td> </tr>
+## <tr> <td align="right"> 2023 </td> <td align="right"> 82 </td> <td align="right"> 81 </td> <td align="right"> 27 </td> <td align="right"> 33 </td> <td align="right"> 0 </td> <td align="right"> 88 </td> <td align="right"> 91 </td> </tr>
+## <tr> <td align="right"> 2024 </td> <td align="right"> 86 </td> <td align="right"> 84 </td> <td align="right"> 29 </td> <td align="right"> 35 </td> <td align="right"> 0 </td> <td align="right"> 92 </td> <td align="right"> 93 </td> </tr>
+## <tr> <td align="right"> 2025 </td> <td align="right"> 90 </td> <td align="right"> 88 </td> <td align="right"> 30 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 96 </td> <td align="right"> 97 </td> </tr>
+## <tr> <td align="right"> 2026 </td> <td align="right"> 90 </td> <td align="right"> 88 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 96 </td> <td align="right"> 97 </td> </tr>
+## <tr> <td align="right"> 2027 </td> <td align="right"> 89 </td> <td align="right"> 86 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 94 </td> <td align="right"> 94 </td> </tr>
+## <tr> <td align="right"> 2028 </td> <td align="right"> 89 </td> <td align="right"> 86 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 95 </td> <td align="right"> 95 </td> </tr>
+## <tr> <td align="right"> 2029 </td> <td align="right"> 85 </td> <td align="right"> 82 </td> <td align="right"> 31 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 90 </td> <td align="right"> 90 </td> </tr>
+## <tr> <td align="right"> 2030 </td> <td align="right"> 84 </td> <td align="right"> 80 </td> <td align="right"> 30 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 88 </td> <td align="right"> 88 </td> </tr>
+## <tr> <td align="right"> 2031 </td> <td align="right"> 85 </td> <td align="right"> 83 </td> <td align="right"> 31 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 91 </td> <td align="right"> 91 </td> </tr>
+## <tr> <td align="right"> 2032 </td> <td align="right"> 88 </td> <td align="right"> 85 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 93 </td> <td align="right"> 93 </td> </tr>
+## </table>
+## <!-- html table generated in R 4.0.0 by xtable 1.8-4 package -->
+## <!-- Tue Jun 9 08:53:57 2020 -->
+## <table border=1>
+## <caption align="top"> Tier 3 projections of BSAI Atka mackerel ABC for the 7 scenarios. </caption>
+## <tr> <th> SSB </th> <th> Scenario.1 </th> <th> Scenario.2 </th> <th> Scenario.3 </th> <th> Scenario.4 </th> <th> Scenario.5 </th> <th> Scenario.6 </th> <th> Scenario.7 </th> </tr>
+## <tr> <td align="right"> 2019 </td> <td align="right"> 172 </td> <td align="right"> 172 </td> <td align="right"> 172 </td> <td align="right"> 172 </td> <td align="right"> 172 </td> <td align="right"> 172 </td> <td align="right"> 172 </td> </tr>
+## <tr> <td align="right"> 2020 </td> <td align="right"> 156 </td> <td align="right"> 164 </td> <td align="right"> 172 </td> <td align="right"> 171 </td> <td align="right"> 177 </td> <td align="right"> 152 </td> <td align="right"> 156 </td> </tr>
+## <tr> <td align="right"> 2021 </td> <td align="right"> 138 </td> <td align="right"> 152 </td> <td align="right"> 178 </td> <td align="right"> 175 </td> <td align="right"> 192 </td> <td align="right"> 131 </td> <td align="right"> 138 </td> </tr>
+## <tr> <td align="right"> 2022 </td> <td align="right"> 130 </td> <td align="right"> 141 </td> <td align="right"> 187 </td> <td align="right"> 182 </td> <td align="right"> 209 </td> <td align="right"> 122 </td> <td align="right"> 128 </td> </tr>
+## <tr> <td align="right"> 2023 </td> <td align="right"> 132 </td> <td align="right"> 142 </td> <td align="right"> 204 </td> <td align="right"> 197 </td> <td align="right"> 233 </td> <td align="right"> 123 </td> <td align="right"> 126 </td> </tr>
+## <tr> <td align="right"> 2024 </td> <td align="right"> 135 </td> <td align="right"> 144 </td> <td align="right"> 218 </td> <td align="right"> 210 </td> <td align="right"> 254 </td> <td align="right"> 125 </td> <td align="right"> 127 </td> </tr>
+## <tr> <td align="right"> 2025 </td> <td align="right"> 139 </td> <td align="right"> 146 </td> <td align="right"> 232 </td> <td align="right"> 223 </td> <td align="right"> 275 </td> <td align="right"> 128 </td> <td align="right"> 129 </td> </tr>
+## <tr> <td align="right"> 2026 </td> <td align="right"> 140 </td> <td align="right"> 147 </td> <td align="right"> 243 </td> <td align="right"> 232 </td> <td align="right"> 292 </td> <td align="right"> 129 </td> <td align="right"> 129 </td> </tr>
+## <tr> <td align="right"> 2027 </td> <td align="right"> 139 </td> <td align="right"> 146 </td> <td align="right"> 249 </td> <td align="right"> 237 </td> <td align="right"> 302 </td> <td align="right"> 127 </td> <td align="right"> 127 </td> </tr>
+## <tr> <td align="right"> 2028 </td> <td align="right"> 138 </td> <td align="right"> 145 </td> <td align="right"> 252 </td> <td align="right"> 240 </td> <td align="right"> 309 </td> <td align="right"> 126 </td> <td align="right"> 126 </td> </tr>
+## <tr> <td align="right"> 2029 </td> <td align="right"> 137 </td> <td align="right"> 144 </td> <td align="right"> 254 </td> <td align="right"> 241 </td> <td align="right"> 314 </td> <td align="right"> 125 </td> <td align="right"> 125 </td> </tr>
+## <tr> <td align="right"> 2030 </td> <td align="right"> 134 </td> <td align="right"> 141 </td> <td align="right"> 251 </td> <td align="right"> 238 </td> <td align="right"> 313 </td> <td align="right"> 122 </td> <td align="right"> 122 </td> </tr>
+## <tr> <td align="right"> 2031 </td> <td align="right"> 134 </td> <td align="right"> 141 </td> <td align="right"> 252 </td> <td align="right"> 238 </td> <td align="right"> 315 </td> <td align="right"> 123 </td> <td align="right"> 123 </td> </tr>
+## <tr> <td align="right"> 2032 </td> <td align="right"> 137 </td> <td align="right"> 144 </td> <td align="right"> 255 </td> <td align="right"> 242 </td> <td align="right"> 319 </td> <td align="right"> 125 </td> <td align="right"> 125 </td> </tr>
+## </table>
+## <!-- html table generated in R 4.0.0 by xtable 1.8-4 package -->
+## <!-- Tue Jun 9 08:53:57 2020 -->
+## <table border=1>
+## <caption align="top"> Tier 3 projections of BSAI Atka mackerel fishing mortality for the 7 scenarios. </caption>
+## <tr> <th> F </th> <th> Scenario.1 </th> <th> Scenario.2 </th> <th> Scenario.3 </th> <th> Scenario.4 </th> <th> Scenario.5 </th> <th> Scenario.6 </th> <th> Scenario.7 </th> </tr>
+## <tr> <td align="right"> 2019 </td> <td align="right"> 0.169 </td> <td align="right"> 0.169 </td> <td align="right"> 0.169 </td> <td align="right"> 0.169 </td> <td align="right"> 0.169 </td> <td align="right"> 0.169 </td> <td align="right"> 0.169 </td> </tr>
+## <tr> <td align="right"> 2020 </td> <td align="right"> 0.299 </td> <td align="right"> 0.172 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.353 </td> <td align="right"> 0.299 </td> </tr>
+## <tr> <td align="right"> 2021 </td> <td align="right"> 0.299 </td> <td align="right"> 0.283 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.343 </td> <td align="right"> 0.299 </td> </tr>
+## <tr> <td align="right"> 2022 </td> <td align="right"> 0.284 </td> <td align="right"> 0.263 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.315 </td> <td align="right"> 0.330 </td> </tr>
+## <tr> <td align="right"> 2023 </td> <td align="right"> 0.277 </td> <td align="right"> 0.259 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.311 </td> <td align="right"> 0.317 </td> </tr>
+## <tr> <td align="right"> 2024 </td> <td align="right"> 0.275 </td> <td align="right"> 0.258 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.310 </td> <td align="right"> 0.312 </td> </tr>
+## <tr> <td align="right"> 2025 </td> <td align="right"> 0.276 </td> <td align="right"> 0.260 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.312 </td> <td align="right"> 0.314 </td> </tr>
+## <tr> <td align="right"> 2026 </td> <td align="right"> 0.277 </td> <td align="right"> 0.261 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.313 </td> <td align="right"> 0.314 </td> </tr>
+## <tr> <td align="right"> 2027 </td> <td align="right"> 0.277 </td> <td align="right"> 0.260 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.313 </td> <td align="right"> 0.313 </td> </tr>
+## <tr> <td align="right"> 2028 </td> <td align="right"> 0.278 </td> <td align="right"> 0.259 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.312 </td> <td align="right"> 0.312 </td> </tr>
+## <tr> <td align="right"> 2029 </td> <td align="right"> 0.277 </td> <td align="right"> 0.258 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.310 </td> <td align="right"> 0.310 </td> </tr>
+## <tr> <td align="right"> 2030 </td> <td align="right"> 0.275 </td> <td align="right"> 0.254 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.306 </td> <td align="right"> 0.306 </td> </tr>
+## <tr> <td align="right"> 2031 </td> <td align="right"> 0.274 </td> <td align="right"> 0.254 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.307 </td> <td align="right"> 0.307 </td> </tr>
+## <tr> <td align="right"> 2032 </td> <td align="right"> 0.275 </td> <td align="right"> 0.256 </td> <td align="right"> 0.066 </td> <td align="right"> 0.083 </td> <td align="right"> 0.000 </td> <td align="right"> 0.309 </td> <td align="right"> 0.309 </td> </tr>
+## </table>
+## <!-- html table generated in R 4.0.0 by xtable 1.8-4 package -->
+## <!-- Tue Jun 9 08:53:57 2020 -->
+## <table border=1>
+## <caption align="top"> Tier 3 projections of BSAI Atka mackerel spawning biomass for the 7 scenarios. </caption>
+## <tr> <th> ABC </th> <th> Scenario.1 </th> <th> Scenario.2 </th> <th> Scenario.3 </th> <th> Scenario.4 </th> <th> Scenario.5 </th> <th> Scenario.6 </th> <th> Scenario.7 </th> </tr>
+## <tr> <td align="right"> 2019 </td> <td align="right"> 93 </td> <td align="right"> 93 </td> <td align="right"> 22 </td> <td align="right"> 27 </td> <td align="right"> 0 </td> <td align="right"> 107 </td> <td align="right"> 107 </td> </tr>
+## <tr> <td align="right"> 2020 </td> <td align="right"> 91 </td> <td align="right"> 89 </td> <td align="right"> 21 </td> <td align="right"> 27 </td> <td align="right"> 0 </td> <td align="right"> 106 </td> <td align="right"> 106 </td> </tr>
+## <tr> <td align="right"> 2021 </td> <td align="right"> 83 </td> <td align="right"> 85 </td> <td align="right"> 22 </td> <td align="right"> 28 </td> <td align="right"> 0 </td> <td align="right"> 91 </td> <td align="right"> 96 </td> </tr>
+## <tr> <td align="right"> 2022 </td> <td align="right"> 81 </td> <td align="right"> 80 </td> <td align="right"> 25 </td> <td align="right"> 30 </td> <td align="right"> 0 </td> <td align="right"> 85 </td> <td align="right"> 93 </td> </tr>
+## <tr> <td align="right"> 2023 </td> <td align="right"> 82 </td> <td align="right"> 81 </td> <td align="right"> 27 </td> <td align="right"> 33 </td> <td align="right"> 0 </td> <td align="right"> 88 </td> <td align="right"> 91 </td> </tr>
+## <tr> <td align="right"> 2024 </td> <td align="right"> 86 </td> <td align="right"> 84 </td> <td align="right"> 29 </td> <td align="right"> 35 </td> <td align="right"> 0 </td> <td align="right"> 92 </td> <td align="right"> 93 </td> </tr>
+## <tr> <td align="right"> 2025 </td> <td align="right"> 90 </td> <td align="right"> 88 </td> <td align="right"> 30 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 96 </td> <td align="right"> 97 </td> </tr>
+## <tr> <td align="right"> 2026 </td> <td align="right"> 90 </td> <td align="right"> 88 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 96 </td> <td align="right"> 97 </td> </tr>
+## <tr> <td align="right"> 2027 </td> <td align="right"> 89 </td> <td align="right"> 86 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 94 </td> <td align="right"> 94 </td> </tr>
+## <tr> <td align="right"> 2028 </td> <td align="right"> 89 </td> <td align="right"> 86 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 95 </td> <td align="right"> 95 </td> </tr>
+## <tr> <td align="right"> 2029 </td> <td align="right"> 85 </td> <td align="right"> 82 </td> <td align="right"> 31 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 90 </td> <td align="right"> 90 </td> </tr>
+## <tr> <td align="right"> 2030 </td> <td align="right"> 84 </td> <td align="right"> 80 </td> <td align="right"> 30 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 88 </td> <td align="right"> 88 </td> </tr>
+## <tr> <td align="right"> 2031 </td> <td align="right"> 85 </td> <td align="right"> 83 </td> <td align="right"> 31 </td> <td align="right"> 37 </td> <td align="right"> 0 </td> <td align="right"> 91 </td> <td align="right"> 91 </td> </tr>
+## <tr> <td align="right"> 2032 </td> <td align="right"> 88 </td> <td align="right"> 85 </td> <td align="right"> 31 </td> <td align="right"> 38 </td> <td align="right"> 0 </td> <td align="right"> 93 </td> <td align="right"> 93 </td> </tr>
+## </table>
+
+
+
+
+
diff --git a/examples/Misc/test.log b/examples/Misc/test.log
new file mode 100644
index 0000000..e04c112
--- /dev/null
+++ b/examples/Misc/test.log
@@ -0,0 +1,686 @@
+This is pdfTeX, Version 3.14159265-2.6-1.40.19 (TeX Live 2018/MacPorts 2018.47642_8) (preloaded format=pdflatex 2018.10.29) 7 JUN 2020 22:19
+entering extended mode
+ restricted \write18 enabled.
+ %&-line parsing enabled.
+**test.tex
+(./test.tex
+LaTeX2e <2018-04-01> patch level 2
+Babel <3.18> and hyphenation patterns for 3 language(s) loaded.
+(/opt/local/share/texmf-texlive/tex/latex/base/article.cls
+Document Class: article 2014/09/29 v1.4h Standard LaTeX document class
+(/opt/local/share/texmf-texlive/tex/latex/base/size10.clo
+File: size10.clo 2014/09/29 v1.4h Standard LaTeX file (size option)
+)
+\c@part=\count80
+\c@section=\count81
+\c@subsection=\count82
+\c@subsubsection=\count83
+\c@paragraph=\count84
+\c@subparagraph=\count85
+\c@figure=\count86
+\c@table=\count87
+\abovecaptionskip=\skip41
+\belowcaptionskip=\skip42
+\bibindent=\dimen102
+) (/opt/local/share/texmf-texlive/tex/latex/lm/lmodern.sty
+Package: lmodern 2009/10/30 v1.6 Latin Modern Fonts
+LaTeX Font Info: Overwriting symbol font `operators' in version `normal'
+(Font) OT1/cmr/m/n --> OT1/lmr/m/n on input line 22.
+LaTeX Font Info: Overwriting symbol font `letters' in version `normal'
+(Font) OML/cmm/m/it --> OML/lmm/m/it on input line 23.
+LaTeX Font Info: Overwriting symbol font `symbols' in version `normal'
+(Font) OMS/cmsy/m/n --> OMS/lmsy/m/n on input line 24.
+LaTeX Font Info: Overwriting symbol font `largesymbols' in version `normal'
+(Font) OMX/cmex/m/n --> OMX/lmex/m/n on input line 25.
+LaTeX Font Info: Overwriting symbol font `operators' in version `bold'
+(Font) OT1/cmr/bx/n --> OT1/lmr/bx/n on input line 26.
+LaTeX Font Info: Overwriting symbol font `letters' in version `bold'
+(Font) OML/cmm/b/it --> OML/lmm/b/it on input line 27.
+LaTeX Font Info: Overwriting symbol font `symbols' in version `bold'
+(Font) OMS/cmsy/b/n --> OMS/lmsy/b/n on input line 28.
+LaTeX Font Info: Overwriting symbol font `largesymbols' in version `bold'
+(Font) OMX/cmex/m/n --> OMX/lmex/m/n on input line 29.
+LaTeX Font Info: Overwriting math alphabet `\mathbf' in version `normal'
+(Font) OT1/cmr/bx/n --> OT1/lmr/bx/n on input line 31.
+LaTeX Font Info: Overwriting math alphabet `\mathsf' in version `normal'
+(Font) OT1/cmss/m/n --> OT1/lmss/m/n on input line 32.
+LaTeX Font Info: Overwriting math alphabet `\mathit' in version `normal'
+(Font) OT1/cmr/m/it --> OT1/lmr/m/it on input line 33.
+LaTeX Font Info: Overwriting math alphabet `\mathtt' in version `normal'
+(Font) OT1/cmtt/m/n --> OT1/lmtt/m/n on input line 34.
+LaTeX Font Info: Overwriting math alphabet `\mathbf' in version `bold'
+(Font) OT1/cmr/bx/n --> OT1/lmr/bx/n on input line 35.
+LaTeX Font Info: Overwriting math alphabet `\mathsf' in version `bold'
+(Font) OT1/cmss/bx/n --> OT1/lmss/bx/n on input line 36.
+LaTeX Font Info: Overwriting math alphabet `\mathit' in version `bold'
+(Font) OT1/cmr/bx/it --> OT1/lmr/bx/it on input line 37.
+LaTeX Font Info: Overwriting math alphabet `\mathtt' in version `bold'
+(Font) OT1/cmtt/m/n --> OT1/lmtt/m/n on input line 38.
+) (/opt/local/share/texmf-texlive/tex/latex/amsfonts/amssymb.sty
+Package: amssymb 2013/01/14 v3.01 AMS font symbols
+(/opt/local/share/texmf-texlive/tex/latex/amsfonts/amsfonts.sty
+Package: amsfonts 2013/01/14 v3.01 Basic AMSFonts support
+\@emptytoks=\toks14
+\symAMSa=\mathgroup4
+\symAMSb=\mathgroup5
+LaTeX Font Info: Overwriting math alphabet `\mathfrak' in version `bold'
+(Font) U/euf/m/n --> U/euf/b/n on input line 106.
+)) (/opt/local/share/texmf-texlive/tex/latex/amsmath/amsmath.sty
+Package: amsmath 2017/09/02 v2.17a AMS math features
+\@mathmargin=\skip43
+For additional information on amsmath, use the `?' option.
+(/opt/local/share/texmf-texlive/tex/latex/amsmath/amstext.sty
+Package: amstext 2000/06/29 v2.01 AMS text
+(/opt/local/share/texmf-texlive/tex/latex/amsmath/amsgen.sty
+File: amsgen.sty 1999/11/30 v2.0 generic functions
+\@emptytoks=\toks15
+\ex@=\dimen103
+)) (/opt/local/share/texmf-texlive/tex/latex/amsmath/amsbsy.sty
+Package: amsbsy 1999/11/29 v1.2d Bold Symbols
+\pmbraise@=\dimen104
+) (/opt/local/share/texmf-texlive/tex/latex/amsmath/amsopn.sty
+Package: amsopn 2016/03/08 v2.02 operator names
+)
+\inf@bad=\count88
+LaTeX Info: Redefining \frac on input line 213.
+\uproot@=\count89
+\leftroot@=\count90
+LaTeX Info: Redefining \overline on input line 375.
+\classnum@=\count91
+\DOTSCASE@=\count92
+LaTeX Info: Redefining \ldots on input line 472.
+LaTeX Info: Redefining \dots on input line 475.
+LaTeX Info: Redefining \cdots on input line 596.
+\Mathstrutbox@=\box26
+\strutbox@=\box27
+\big@size=\dimen105
+LaTeX Font Info: Redeclaring font encoding OML on input line 712.
+LaTeX Font Info: Redeclaring font encoding OMS on input line 713.
+\macc@depth=\count93
+\c@MaxMatrixCols=\count94
+\dotsspace@=\muskip10
+\c@parentequation=\count95
+\dspbrk@lvl=\count96
+\tag@help=\toks16
+\row@=\count97
+\column@=\count98
+\maxfields@=\count99
+\andhelp@=\toks17
+\eqnshift@=\dimen106
+\alignsep@=\dimen107
+\tagshift@=\dimen108
+\tagwidth@=\dimen109
+\totwidth@=\dimen110
+\lineht@=\dimen111
+\@envbody=\toks18
+\multlinegap=\skip44
+\multlinetaggap=\skip45
+\mathdisplay@stack=\toks19
+LaTeX Info: Redefining \[ on input line 2817.
+LaTeX Info: Redefining \] on input line 2818.
+) (/opt/local/share/texmf-texlive/tex/generic/ifxetex/ifxetex.sty
+Package: ifxetex 2010/09/12 v0.6 Provides ifxetex conditional
+) (/opt/local/share/texmf-texlive/tex/generic/oberdiek/ifluatex.sty
+Package: ifluatex 2016/05/16 v1.4 Provides the ifluatex switch (HO)
+Package ifluatex Info: LuaTeX not detected.
+) (/opt/local/share/texmf-texlive/tex/latex/base/fixltx2e.sty
+Package: fixltx2e 2016/12/29 v2.1a fixes to LaTeX (obsolete)
+Applying: [2015/01/01] Old fixltx2e package on input line 46.
+
+Package fixltx2e Warning: fixltx2e is not required with releases after 2015
+(fixltx2e) All fixes are now in the LaTeX kernel.
+(fixltx2e) See the latexrelease package for details.
+
+Already applied: [0000/00/00] Old fixltx2e package on input line 53.
+) (/opt/local/share/texmf-texlive/tex/latex/base/fontenc.sty
+Package: fontenc 2017/04/05 v2.0i Standard LaTeX package
+(/opt/local/share/texmf-texlive/tex/latex/base/t1enc.def
+File: t1enc.def 2017/04/05 v2.0i Standard LaTeX file
+LaTeX Font Info: Redeclaring font encoding T1 on input line 48.
+)) (/opt/local/share/texmf-texlive/tex/latex/base/inputenc.sty
+Package: inputenc 2018/04/06 v1.3b Input encoding file
+\inpenc@prehook=\toks20
+\inpenc@posthook=\toks21
+) (/opt/local/share/texmf-texlive/tex/latex/base/textcomp.sty
+Package: textcomp 2017/04/05 v2.0i Standard LaTeX package
+Package textcomp Info: Sub-encoding information:
+(textcomp) 5 = only ISO-Adobe without \textcurrency
+(textcomp) 4 = 5 + \texteuro
+(textcomp) 3 = 4 + \textohm
+(textcomp) 2 = 3 + \textestimated + \textcurrency
+(textcomp) 1 = TS1 - \textcircled - \t
+(textcomp) 0 = TS1 (full)
+(textcomp) Font families with sub-encoding setting implement
+(textcomp) only a restricted character set as indicated.
+(textcomp) Family '?' is the default used for unknown fonts.
+(textcomp) See the documentation for details.
+Package textcomp Info: Setting ? sub-encoding to TS1/1 on input line 79.
+(/opt/local/share/texmf-texlive/tex/latex/base/ts1enc.def
+File: ts1enc.def 2001/06/05 v3.0e (jk/car/fm) Standard LaTeX file
+Now handling font encoding TS1 ...
+... processing UTF-8 mapping file for font encoding TS1
+(/opt/local/share/texmf-texlive/tex/latex/base/ts1enc.dfu
+File: ts1enc.dfu 2018/04/05 v1.2c UTF-8 support for inputenc
+ defining Unicode char U+00A2 (decimal 162)
+ defining Unicode char U+00A3 (decimal 163)
+ defining Unicode char U+00A4 (decimal 164)
+ defining Unicode char U+00A5 (decimal 165)
+ defining Unicode char U+00A6 (decimal 166)
+ defining Unicode char U+00A7 (decimal 167)
+ defining Unicode char U+00A8 (decimal 168)
+ defining Unicode char U+00A9 (decimal 169)
+ defining Unicode char U+00AA (decimal 170)
+ defining Unicode char U+00AC (decimal 172)
+ defining Unicode char U+00AE (decimal 174)
+ defining Unicode char U+00AF (decimal 175)
+ defining Unicode char U+00B0 (decimal 176)
+ defining Unicode char U+00B1 (decimal 177)
+ defining Unicode char U+00B2 (decimal 178)
+ defining Unicode char U+00B3 (decimal 179)
+ defining Unicode char U+00B4 (decimal 180)
+ defining Unicode char U+00B5 (decimal 181)
+ defining Unicode char U+00B6 (decimal 182)
+ defining Unicode char U+00B7 (decimal 183)
+ defining Unicode char U+00B9 (decimal 185)
+ defining Unicode char U+00BA (decimal 186)
+ defining Unicode char U+00BC (decimal 188)
+ defining Unicode char U+00BD (decimal 189)
+ defining Unicode char U+00BE (decimal 190)
+ defining Unicode char U+00D7 (decimal 215)
+ defining Unicode char U+00F7 (decimal 247)
+ defining Unicode char U+0192 (decimal 402)
+ defining Unicode char U+02C7 (decimal 711)
+ defining Unicode char U+02D8 (decimal 728)
+ defining Unicode char U+02DD (decimal 733)
+ defining Unicode char U+0E3F (decimal 3647)
+ defining Unicode char U+2016 (decimal 8214)
+ defining Unicode char U+2020 (decimal 8224)
+ defining Unicode char U+2021 (decimal 8225)
+ defining Unicode char U+2022 (decimal 8226)
+ defining Unicode char U+2030 (decimal 8240)
+ defining Unicode char U+2031 (decimal 8241)
+ defining Unicode char U+203B (decimal 8251)
+ defining Unicode char U+203D (decimal 8253)
+ defining Unicode char U+2044 (decimal 8260)
+ defining Unicode char U+204E (decimal 8270)
+ defining Unicode char U+2052 (decimal 8274)
+ defining Unicode char U+20A1 (decimal 8353)
+ defining Unicode char U+20A4 (decimal 8356)
+ defining Unicode char U+20A6 (decimal 8358)
+ defining Unicode char U+20A9 (decimal 8361)
+ defining Unicode char U+20AB (decimal 8363)
+ defining Unicode char U+20AC (decimal 8364)
+ defining Unicode char U+20B1 (decimal 8369)
+ defining Unicode char U+2103 (decimal 8451)
+ defining Unicode char U+2116 (decimal 8470)
+ defining Unicode char U+2117 (decimal 8471)
+ defining Unicode char U+211E (decimal 8478)
+ defining Unicode char U+2120 (decimal 8480)
+ defining Unicode char U+2122 (decimal 8482)
+ defining Unicode char U+2126 (decimal 8486)
+ defining Unicode char U+2127 (decimal 8487)
+ defining Unicode char U+212E (decimal 8494)
+ defining Unicode char U+2190 (decimal 8592)
+ defining Unicode char U+2191 (decimal 8593)
+ defining Unicode char U+2192 (decimal 8594)
+ defining Unicode char U+2193 (decimal 8595)
+ defining Unicode char U+2329 (decimal 9001)
+ defining Unicode char U+232A (decimal 9002)
+ defining Unicode char U+2422 (decimal 9250)
+ defining Unicode char U+25E6 (decimal 9702)
+ defining Unicode char U+25EF (decimal 9711)
+ defining Unicode char U+266A (decimal 9834)
+ defining Unicode char U+FEFF (decimal 65279)
+))
+LaTeX Info: Redefining \oldstylenums on input line 334.
+Package textcomp Info: Setting cmr sub-encoding to TS1/0 on input line 349.
+Package textcomp Info: Setting cmss sub-encoding to TS1/0 on input line 350.
+Package textcomp Info: Setting cmtt sub-encoding to TS1/0 on input line 351.
+Package textcomp Info: Setting cmvtt sub-encoding to TS1/0 on input line 352.
+Package textcomp Info: Setting cmbr sub-encoding to TS1/0 on input line 353.
+Package textcomp Info: Setting cmtl sub-encoding to TS1/0 on input line 354.
+Package textcomp Info: Setting ccr sub-encoding to TS1/0 on input line 355.
+Package textcomp Info: Setting ptm sub-encoding to TS1/4 on input line 356.
+Package textcomp Info: Setting pcr sub-encoding to TS1/4 on input line 357.
+Package textcomp Info: Setting phv sub-encoding to TS1/4 on input line 358.
+Package textcomp Info: Setting ppl sub-encoding to TS1/3 on input line 359.
+Package textcomp Info: Setting pag sub-encoding to TS1/4 on input line 360.
+Package textcomp Info: Setting pbk sub-encoding to TS1/4 on input line 361.
+Package textcomp Info: Setting pnc sub-encoding to TS1/4 on input line 362.
+Package textcomp Info: Setting pzc sub-encoding to TS1/4 on input line 363.
+Package textcomp Info: Setting bch sub-encoding to TS1/4 on input line 364.
+Package textcomp Info: Setting put sub-encoding to TS1/5 on input line 365.
+Package textcomp Info: Setting uag sub-encoding to TS1/5 on input line 366.
+Package textcomp Info: Setting ugq sub-encoding to TS1/5 on input line 367.
+Package textcomp Info: Setting ul8 sub-encoding to TS1/4 on input line 368.
+Package textcomp Info: Setting ul9 sub-encoding to TS1/4 on input line 369.
+Package textcomp Info: Setting augie sub-encoding to TS1/5 on input line 370.
+Package textcomp Info: Setting dayrom sub-encoding to TS1/3 on input line 371.
+Package textcomp Info: Setting dayroms sub-encoding to TS1/3 on input line 372.
+
+Package textcomp Info: Setting pxr sub-encoding to TS1/0 on input line 373.
+Package textcomp Info: Setting pxss sub-encoding to TS1/0 on input line 374.
+Package textcomp Info: Setting pxtt sub-encoding to TS1/0 on input line 375.
+Package textcomp Info: Setting txr sub-encoding to TS1/0 on input line 376.
+Package textcomp Info: Setting txss sub-encoding to TS1/0 on input line 377.
+Package textcomp Info: Setting txtt sub-encoding to TS1/0 on input line 378.
+Package textcomp Info: Setting lmr sub-encoding to TS1/0 on input line 379.
+Package textcomp Info: Setting lmdh sub-encoding to TS1/0 on input line 380.
+Package textcomp Info: Setting lmss sub-encoding to TS1/0 on input line 381.
+Package textcomp Info: Setting lmssq sub-encoding to TS1/0 on input line 382.
+Package textcomp Info: Setting lmvtt sub-encoding to TS1/0 on input line 383.
+Package textcomp Info: Setting lmtt sub-encoding to TS1/0 on input line 384.
+Package textcomp Info: Setting qhv sub-encoding to TS1/0 on input line 385.
+Package textcomp Info: Setting qag sub-encoding to TS1/0 on input line 386.
+Package textcomp Info: Setting qbk sub-encoding to TS1/0 on input line 387.
+Package textcomp Info: Setting qcr sub-encoding to TS1/0 on input line 388.
+Package textcomp Info: Setting qcs sub-encoding to TS1/0 on input line 389.
+Package textcomp Info: Setting qpl sub-encoding to TS1/0 on input line 390.
+Package textcomp Info: Setting qtm sub-encoding to TS1/0 on input line 391.
+Package textcomp Info: Setting qzc sub-encoding to TS1/0 on input line 392.
+Package textcomp Info: Setting qhvc sub-encoding to TS1/0 on input line 393.
+Package textcomp Info: Setting futs sub-encoding to TS1/4 on input line 394.
+Package textcomp Info: Setting futx sub-encoding to TS1/4 on input line 395.
+Package textcomp Info: Setting futj sub-encoding to TS1/4 on input line 396.
+Package textcomp Info: Setting hlh sub-encoding to TS1/3 on input line 397.
+Package textcomp Info: Setting hls sub-encoding to TS1/3 on input line 398.
+Package textcomp Info: Setting hlst sub-encoding to TS1/3 on input line 399.
+Package textcomp Info: Setting hlct sub-encoding to TS1/5 on input line 400.
+Package textcomp Info: Setting hlx sub-encoding to TS1/5 on input line 401.
+Package textcomp Info: Setting hlce sub-encoding to TS1/5 on input line 402.
+Package textcomp Info: Setting hlcn sub-encoding to TS1/5 on input line 403.
+Package textcomp Info: Setting hlcw sub-encoding to TS1/5 on input line 404.
+Package textcomp Info: Setting hlcf sub-encoding to TS1/5 on input line 405.
+Package textcomp Info: Setting pplx sub-encoding to TS1/3 on input line 406.
+Package textcomp Info: Setting pplj sub-encoding to TS1/3 on input line 407.
+Package textcomp Info: Setting ptmx sub-encoding to TS1/4 on input line 408.
+Package textcomp Info: Setting ptmj sub-encoding to TS1/4 on input line 409.
+) (/opt/local/share/texmf-texlive/tex/latex/upquote/upquote.sty
+Package: upquote 2012/04/19 v1.3 upright-quote and grave-accent glyphs in verba
+tim
+) (/opt/local/share/texmf-texlive/tex/latex/microtype/microtype.sty
+Package: microtype 2018/01/14 v2.7a Micro-typographical refinements (RS)
+(/opt/local/share/texmf-texlive/tex/latex/graphics/keyval.sty
+Package: keyval 2014/10/28 v1.15 key=value parser (DPC)
+\KV@toks@=\toks22
+)
+\MT@toks=\toks23
+\MT@count=\count100
+LaTeX Info: Redefining \textls on input line 793.
+\MT@outer@kern=\dimen112
+LaTeX Info: Redefining \textmicrotypecontext on input line 1339.
+\MT@listname@count=\count101
+(/opt/local/share/texmf-texlive/tex/latex/microtype/microtype-pdftex.def
+File: microtype-pdftex.def 2018/01/14 v2.7a Definitions specific to pdftex (RS)
+
+LaTeX Info: Redefining \lsstyle on input line 913.
+LaTeX Info: Redefining \lslig on input line 913.
+\MT@outer@space=\skip46
+)
+Package microtype Info: Loading configuration file microtype.cfg.
+(/opt/local/share/texmf-texlive/tex/latex/microtype/microtype.cfg
+File: microtype.cfg 2018/01/14 v2.7a microtype main configuration file (RS)
+)) (/opt/local/share/texmf-texlive/tex/latex/parskip/parskip.sty
+Package: parskip 2001/04/09 non-zero parskip adjustments
+) (/opt/local/share/texmf-texlive/tex/latex/hyperref/hyperref.sty
+Package: hyperref 2018/02/06 v6.86b Hypertext links for LaTeX
+(/opt/local/share/texmf-texlive/tex/generic/oberdiek/hobsub-hyperref.sty
+Package: hobsub-hyperref 2016/05/16 v1.14 Bundle oberdiek, subset hyperref (HO)
+
+(/opt/local/share/texmf-texlive/tex/generic/oberdiek/hobsub-generic.sty
+Package: hobsub-generic 2016/05/16 v1.14 Bundle oberdiek, subset generic (HO)
+Package: hobsub 2016/05/16 v1.14 Construct package bundles (HO)
+Package: infwarerr 2016/05/16 v1.4 Providing info/warning/error messages (HO)
+Package: ltxcmds 2016/05/16 v1.23 LaTeX kernel commands for general use (HO)
+Package hobsub Info: Skipping package `ifluatex' (already loaded).
+Package: ifvtex 2016/05/16 v1.6 Detect VTeX and its facilities (HO)
+Package ifvtex Info: VTeX not detected.
+Package: intcalc 2016/05/16 v1.2 Expandable calculations with integers (HO)
+Package: ifpdf 2017/03/15 v3.2 Provides the ifpdf switch
+Package: etexcmds 2016/05/16 v1.6 Avoid name clashes with e-TeX commands (HO)
+Package etexcmds Info: Could not find \expanded.
+(etexcmds) That can mean that you are not using pdfTeX 1.50 or
+(etexcmds) that some package has redefined \expanded.
+(etexcmds) In the latter case, load this package earlier.
+Package: kvsetkeys 2016/05/16 v1.17 Key value parser (HO)
+Package: kvdefinekeys 2016/05/16 v1.4 Define keys (HO)
+Package: pdftexcmds 2018/01/30 v0.27 Utility functions of pdfTeX for LuaTeX (HO
+)
+Package pdftexcmds Info: LuaTeX not detected.
+Package pdftexcmds Info: \pdf@primitive is available.
+Package pdftexcmds Info: \pdf@ifprimitive is available.
+Package pdftexcmds Info: \pdfdraftmode found.
+Package: pdfescape 2016/05/16 v1.14 Implements pdfTeX's escape features (HO)
+Package: bigintcalc 2016/05/16 v1.4 Expandable calculations on big integers (HO
+)
+Package: bitset 2016/05/16 v1.2 Handle bit-vector datatype (HO)
+Package: uniquecounter 2016/05/16 v1.3 Provide unlimited unique counter (HO)
+)
+Package hobsub Info: Skipping package `hobsub' (already loaded).
+Package: letltxmacro 2016/05/16 v1.5 Let assignment for LaTeX macros (HO)
+Package: hopatch 2016/05/16 v1.3 Wrapper for package hooks (HO)
+Package: xcolor-patch 2016/05/16 xcolor patch
+Package: atveryend 2016/05/16 v1.9 Hooks at the very end of document (HO)
+Package atveryend Info: \enddocument detected (standard20110627).
+Package: atbegshi 2016/06/09 v1.18 At begin shipout hook (HO)
+Package: refcount 2016/05/16 v3.5 Data extraction from label references (HO)
+Package: hycolor 2016/05/16 v1.8 Color options for hyperref/bookmark (HO)
+) (/opt/local/share/texmf-texlive/tex/latex/oberdiek/auxhook.sty
+Package: auxhook 2016/05/16 v1.4 Hooks for auxiliary files (HO)
+) (/opt/local/share/texmf-texlive/tex/latex/oberdiek/kvoptions.sty
+Package: kvoptions 2016/05/16 v3.12 Key value format for package options (HO)
+)
+\@linkdim=\dimen113
+\Hy@linkcounter=\count102
+\Hy@pagecounter=\count103
+(/opt/local/share/texmf-texlive/tex/latex/hyperref/pd1enc.def
+File: pd1enc.def 2018/02/06 v6.86b Hyperref: PDFDocEncoding definition (HO)
+Now handling font encoding PD1 ...
+... no UTF-8 mapping file for font encoding PD1
+)
+\Hy@SavedSpaceFactor=\count104
+(/opt/local/share/texmf-texlive/tex/latex/latexconfig/hyperref.cfg
+File: hyperref.cfg 2002/06/06 v1.2 hyperref configuration of TeXLive
+)
+Package hyperref Info: Option `unicode' set `true' on input line 4383.
+(/opt/local/share/texmf-texlive/tex/latex/hyperref/puenc.def
+File: puenc.def 2018/02/06 v6.86b Hyperref: PDF Unicode definition (HO)
+Now handling font encoding PU ...
+... no UTF-8 mapping file for font encoding PU
+)
+Package hyperref Info: Hyper figures OFF on input line 4509.
+Package hyperref Info: Link nesting OFF on input line 4514.
+Package hyperref Info: Hyper index ON on input line 4517.
+Package hyperref Info: Plain pages OFF on input line 4524.
+Package hyperref Info: Backreferencing OFF on input line 4529.
+Package hyperref Info: Implicit mode ON; LaTeX internals redefined.
+Package hyperref Info: Bookmarks ON on input line 4762.
+\c@Hy@tempcnt=\count105
+(/opt/local/share/texmf-texlive/tex/latex/url/url.sty
+\Urlmuskip=\muskip11
+Package: url 2013/09/16 ver 3.4 Verb mode for urls, etc.
+)
+LaTeX Info: Redefining \url on input line 5115.
+\XeTeXLinkMargin=\dimen114
+\Fld@menulength=\count106
+\Field@Width=\dimen115
+\Fld@charsize=\dimen116
+Package hyperref Info: Hyper figures OFF on input line 6369.
+Package hyperref Info: Link nesting OFF on input line 6374.
+Package hyperref Info: Hyper index ON on input line 6377.
+Package hyperref Info: backreferencing OFF on input line 6384.
+Package hyperref Info: Link coloring OFF on input line 6389.
+Package hyperref Info: Link coloring with OCG OFF on input line 6394.
+Package hyperref Info: PDF/A mode OFF on input line 6399.
+LaTeX Info: Redefining \ref on input line 6439.
+LaTeX Info: Redefining \pageref on input line 6443.
+\Hy@abspage=\count107
+\c@Item=\count108
+\c@Hfootnote=\count109
+)
+Package hyperref Info: Driver (autodetected): hpdftex.
+(/opt/local/share/texmf-texlive/tex/latex/hyperref/hpdftex.def
+File: hpdftex.def 2018/02/06 v6.86b Hyperref driver for pdfTeX
+\Fld@listcount=\count110
+\c@bookmark@seq@number=\count111
+(/opt/local/share/texmf-texlive/tex/latex/oberdiek/rerunfilecheck.sty
+Package: rerunfilecheck 2016/05/16 v1.8 Rerun checks for auxiliary files (HO)
+Package uniquecounter Info: New unique counter `rerunfilecheck' on input line 2
+82.
+)
+\Hy@SectionHShift=\skip47
+)
+Package hyperref Info: Option `breaklinks' set `true' on input line 35.
+(/opt/local/share/texmf-texlive/tex/latex/geometry/geometry.sty
+Package: geometry 2018/03/24 v5.7 Page Geometry
+\Gm@cnth=\count112
+\Gm@cntv=\count113
+\c@Gm@tempcnt=\count114
+\Gm@bindingoffset=\dimen117
+\Gm@wd@mp=\dimen118
+\Gm@odd@mp=\dimen119
+\Gm@even@mp=\dimen120
+\Gm@layoutwidth=\dimen121
+\Gm@layoutheight=\dimen122
+\Gm@layouthoffset=\dimen123
+\Gm@layoutvoffset=\dimen124
+\Gm@dimlist=\toks24
+) (/opt/local/share/texmf-texlive/tex/latex/graphics/color.sty
+Package: color 2016/07/10 v1.1e Standard LaTeX Color (DPC)
+(/opt/local/share/texmf-texlive/tex/latex/graphics-cfg/color.cfg
+File: color.cfg 2016/01/02 v1.6 sample color configuration
+)
+Package color Info: Driver file: pdftex.def on input line 147.
+(/opt/local/share/texmf-texlive/tex/latex/graphics-def/pdftex.def
+File: pdftex.def 2018/01/08 v1.0l Graphics/color driver for pdftex
+)) (/opt/local/share/texmf-texlive/tex/latex/fancyvrb/fancyvrb.sty
+Package: fancyvrb 2008/02/07
+
+Style option: `fancyvrb' v2.7a, with DG/SPQR fixes, and firstline=lastline fix
+<2008/02/07> (tvz)
+\FV@CodeLineNo=\count115
+\FV@InFile=\read1
+\FV@TabBox=\box28
+\c@FancyVerbLine=\count116
+\FV@StepNumber=\count117
+\FV@OutFile=\write3
+) (/opt/local/share/texmf-texlive/tex/latex/framed/framed.sty
+Package: framed 2011/10/22 v 0.96: framed or shaded text with page breaks
+\OuterFrameSep=\skip48
+\fb@frw=\dimen125
+\fb@frh=\dimen126
+\FrameRule=\dimen127
+\FrameSep=\dimen128
+) (/opt/local/share/texmf-texlive/tex/latex/graphics/graphicx.sty
+Package: graphicx 2017/06/01 v1.1a Enhanced LaTeX Graphics (DPC,SPQR)
+(/opt/local/share/texmf-texlive/tex/latex/graphics/graphics.sty
+Package: graphics 2017/06/25 v1.2c Standard LaTeX Graphics (DPC,SPQR)
+(/opt/local/share/texmf-texlive/tex/latex/graphics/trig.sty
+Package: trig 2016/01/03 v1.10 sin cos tan (DPC)
+) (/opt/local/share/texmf-texlive/tex/latex/graphics-cfg/graphics.cfg
+File: graphics.cfg 2016/06/04 v1.11 sample graphics configuration
+)
+Package graphics Info: Driver file: pdftex.def on input line 99.
+)
+\Gin@req@height=\dimen129
+\Gin@req@width=\dimen130
+) (/opt/local/share/texmf-texlive/tex/latex/oberdiek/grffile.sty
+Package: grffile 2017/06/30 v1.18 Extended file name support for graphics (HO)
+Package grffile Info: Option `multidot' is set to `true'.
+Package grffile Info: Option `extendedchars' is set to `false'.
+Package grffile Info: Option `space' is set to `true'.
+Package grffile Info: \Gin@ii of package `graphicx' fixed on input line 494.
+)
+No file test.aux.
+\openout1 = `test.aux'.
+
+LaTeX Font Info: Checking defaults for OML/cmm/m/it on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for T1/cmr/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for OT1/cmr/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for OMS/cmsy/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for OMX/cmex/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for U/cmr/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for TS1/cmr/m/n on input line 111.
+LaTeX Font Info: Try loading font information for TS1+cmr on input line 111.
+
+(/opt/local/share/texmf-texlive/tex/latex/base/ts1cmr.fd
+File: ts1cmr.fd 2014/09/29 v2.5h Standard LaTeX font definitions
+)
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for PD1/pdf/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Checking defaults for PU/pdf/m/n on input line 111.
+LaTeX Font Info: ... okay on input line 111.
+LaTeX Font Info: Try loading font information for T1+lmr on input line 111.
+(/opt/local/share/texmf-texlive/tex/latex/lm/t1lmr.fd
+File: t1lmr.fd 2009/10/30 v1.6 Font defs for Latin Modern
+)
+LaTeX Info: Redefining \microtypecontext on input line 111.
+Package microtype Info: Generating PDF output.
+Package microtype Info: Character protrusion enabled (level 2).
+Package microtype Info: Using protrusion set `basicmath'.
+Package microtype Info: Automatic font expansion enabled (level 2),
+(microtype) stretch: 20, shrink: 20, step: 1, non-selected.
+Package microtype Info: Using default expansion set `basictext'.
+Package microtype Info: No adjustment of tracking.
+Package microtype Info: No adjustment of interword spacing.
+Package microtype Info: No adjustment of character kerning.
+(/opt/local/share/texmf-texlive/tex/latex/microtype/mt-cmr.cfg
+File: mt-cmr.cfg 2013/05/19 v2.2 microtype config. file: Computer Modern Roman
+(RS)
+)
+\AtBeginShipoutBox=\box29
+Package hyperref Info: Link coloring OFF on input line 111.
+(/opt/local/share/texmf-texlive/tex/latex/hyperref/nameref.sty
+Package: nameref 2016/05/21 v2.44 Cross-referencing by name of section
+(/opt/local/share/texmf-texlive/tex/generic/oberdiek/gettitlestring.sty
+Package: gettitlestring 2016/05/16 v1.5 Cleanup title references (HO)
+)
+\c@section@level=\count118
+)
+LaTeX Info: Redefining \ref on input line 111.
+LaTeX Info: Redefining \pageref on input line 111.
+LaTeX Info: Redefining \nameref on input line 111.
+\@outlinefile=\write4
+\openout4 = `test.out'.
+
+*geometry* driver: auto-detecting
+*geometry* detected driver: pdftex
+*geometry* verbose mode - [ preamble ] result:
+* driver: pdftex
+* paper:
+* layout:
+* layoutoffset:(h,v)=(0.0pt,0.0pt)
+* modes:
+* h-part:(L,W,R)=(72.26999pt, 469.75502pt, 72.26999pt)
+* v-part:(T,H,B)=(72.26999pt, 650.43001pt, 72.26999pt)
+* \paperwidth=614.295pt
+* \paperheight=794.96999pt
+* \textwidth=469.75502pt
+* \textheight=650.43001pt
+* \oddsidemargin=0.0pt
+* \evensidemargin=0.0pt
+* \topmargin=-37.0pt
+* \headheight=12.0pt
+* \headsep=25.0pt
+* \topskip=10.0pt
+* \footskip=30.0pt
+* \marginparwidth=65.0pt
+* \marginparsep=11.0pt
+* \columnsep=10.0pt
+* \skip\footins=9.0pt plus 4.0pt minus 2.0pt
+* \hoffset=0.0pt
+* \voffset=0.0pt
+* \mag=1000
+* \@twocolumnfalse
+* \@twosidefalse
+* \@mparswitchfalse
+* \@reversemarginfalse
+* (1in=72.27pt=25.4mm, 1cm=28.453pt)
+
+(/opt/local/share/texmf-texlive/tex/context/base/mkii/supp-pdf.mkii
+[Loading MPS to PDF converter (version 2006.09.02).]
+\scratchcounter=\count119
+\scratchdimen=\dimen131
+\scratchbox=\box30
+\nofMPsegments=\count120
+\nofMParguments=\count121
+\everyMPshowfont=\toks25
+\MPscratchCnt=\count122
+\MPscratchDim=\dimen132
+\MPnumerator=\count123
+\makeMPintoPDFobject=\count124
+\everyMPtoPDFconversion=\toks26
+) (/opt/local/share/texmf-texlive/tex/latex/oberdiek/epstopdf-base.sty
+Package: epstopdf-base 2016/05/15 v2.6 Base part for package epstopdf
+(/opt/local/share/texmf-texlive/tex/latex/oberdiek/grfext.sty
+Package: grfext 2016/05/16 v1.2 Manage graphics extensions (HO)
+)
+Package epstopdf-base Info: Redefining graphics rule for `.eps' on input line 4
+38.
+Package grfext Info: Graphics extension search list:
+(grfext) [.pdf,.png,.jpg,.mps,.jpeg,.jbig2,.jb2,.PDF,.PNG,.JPG,.JPE
+G,.JBIG2,.JB2,.eps]
+(grfext) \AppendGraphicsExtensions on input line 456.
+(/opt/local/share/texmf-texlive/tex/latex/latexconfig/epstopdf-sys.cfg
+File: epstopdf-sys.cfg 2010/07/13 v1.3 Configuration of (r)epstopdf for TeX Liv
+e
+))
+LaTeX Font Info: Try loading font information for OT1+lmr on input line 113.
+
+(/opt/local/share/texmf-texlive/tex/latex/lm/ot1lmr.fd
+File: ot1lmr.fd 2009/10/30 v1.6 Font defs for Latin Modern
+)
+LaTeX Font Info: Try loading font information for OML+lmm on input line 113.
+
+(/opt/local/share/texmf-texlive/tex/latex/lm/omllmm.fd
+File: omllmm.fd 2009/10/30 v1.6 Font defs for Latin Modern
+)
+LaTeX Font Info: Try loading font information for OMS+lmsy on input line 113
+.
+(/opt/local/share/texmf-texlive/tex/latex/lm/omslmsy.fd
+File: omslmsy.fd 2009/10/30 v1.6 Font defs for Latin Modern
+)
+LaTeX Font Info: Try loading font information for OMX+lmex on input line 113
+.
+(/opt/local/share/texmf-texlive/tex/latex/lm/omxlmex.fd
+File: omxlmex.fd 2009/10/30 v1.6 Font defs for Latin Modern
+)
+LaTeX Font Info: External font `lmex10' loaded for size
+(Font) <12> on input line 113.
+LaTeX Font Info: External font `lmex10' loaded for size
+(Font) <8> on input line 113.
+LaTeX Font Info: External font `lmex10' loaded for size
+(Font) <6> on input line 113.
+LaTeX Font Info: Try loading font information for U+msa on input line 113.
+(/opt/local/share/texmf-texlive/tex/latex/amsfonts/umsa.fd
+File: umsa.fd 2013/01/14 v3.01 AMS symbols A
+) (/opt/local/share/texmf-texlive/tex/latex/microtype/mt-msa.cfg
+File: mt-msa.cfg 2006/02/04 v1.1 microtype config. file: AMS symbols (a) (RS)
+)
+LaTeX Font Info: Try loading font information for U+msb on input line 113.
+(/opt/local/share/texmf-texlive/tex/latex/amsfonts/umsb.fd
+File: umsb.fd 2013/01/14 v3.01 AMS symbols B
+) (/opt/local/share/texmf-texlive/tex/latex/microtype/mt-msb.cfg
+File: mt-msb.cfg 2005/06/01 v1.0 microtype config. file: AMS symbols (b) (RS)
+)
+LaTeX Font Info: Try loading font information for T1+lmtt on input line 115.
+
+(/opt/local/share/texmf-texlive/tex/latex/lm/t1lmtt.fd
+File: t1lmtt.fd 2009/10/30 v1.6 Font defs for Latin Modern
+)
+LaTeX Font Info: Font shape `T1/lmtt/bx/n' in size <10> not available
+(Font) Font shape `T1/lmtt/b/n' tried instead on input line 117.
+LaTeX Font Info: Try loading font information for TS1+lmtt on input line 137
+.
+(/opt/local/share/texmf-texlive/tex/latex/lm/ts1lmtt.fd
+File: ts1lmtt.fd 2009/10/30 v1.6 Font defs for Latin Modern
+) [1
+
+{/opt/local/var/db/texmf/fonts/map/pdftex/updmap/pdftex.map}]
+
+! Package inputenc Error: Unicode character ^^[ (U+1B)
+(inputenc) not set up for use with LaTeX.
+
+See the inputenc package documentation for explanation.
+Type H for immediate help.
+ ...
+
+l.280 \end{verbatim}
+
+Here is how much of TeX's memory you used:
+ 12391 strings out of 494590
+ 181660 string characters out of 6177728
+ 292211 words of memory out of 5000000
+ 15859 multiletter control sequences out of 15000+600000
+ 27919 words of font info for 47 fonts, out of 8000000 for 9000
+ 14 hyphenation exceptions out of 8191
+ 32i,6n,38p,1296b,373s stack positions out of 5000i,500n,10000p,200000b,80000s
+
+! ==> Fatal error occurred, no output PDF file produced!
diff --git a/examples/Misc/test.md b/examples/Misc/test.md
new file mode 100644
index 0000000..08459f8
--- /dev/null
+++ b/examples/Misc/test.md
@@ -0,0 +1,275 @@
+# Projection model work
+Projection model work for demonstrating application of controls and input data
+
+
+```r
+library(tidyverse)
+library(ggplot2)
+library(ggthemes)
+library(gtable)
+source("../R/readData.R")
+```
+
+## Set initial "setup" parameters
+
+
+```r
+thisyr=2019
+setup<-list(
+ Run_name = noquote("Std"),
+ Tier = 3 ,
+ nalts = 7 ,
+ alts = c(1,2,3,4,5,6,7),
+ tac_abc = 1, #' Flag to set TAC equal to ABC (1 means true, otherwise false)
+ srr = 1 , #' Stock-recruitment type (1=Ricker, 2=Bholt)
+ rec_proj = 1, #' projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+ srr_cond = 0 , #' SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+ srr_prior = 0.0, #' Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+ write_big = 1, #' Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+ nyrs_proj = 14, #' Number of projection years
+ nsims = 100, #' Number of simulations
+ beg_yr_label = thisyr #' Begin Year
+)
+```
+
+## Set up the species specific run file
+
+
+```r
+config<-list(
+ nFixCatchYrs = 2,
+ nSpecies = 1,
+ OYMin = .1343248,
+ OYMax = 1943248,
+ dataFiles = noquote("data/t1.dat"),
+ ABCMult = 1,
+ PoplnScalar = 1000,
+ AltFabcSPR = 0.75,
+ nTAC = 1,
+ TACIndices = 1,
+ Catch = c( 2016,55000., 2017,55000. )
+)
+```
+
+## Make list of main file w/ assessment model results
+E.g., "data/bsai_atka.dat"
+
+
+```r
+datfile <- list(
+ runname = noquote("M16.2"),
+ ssl_spp = 1, # SSL_spp
+ Dorn_buffer = 1, # Dorn_buffer
+ nfsh = 1, # N_fsh
+ nsex = 1, # N_sexes
+ avgF5yr = 0.0661399, # avg_5yr_F
+ F40_mult = 1, # F_40_multiplier
+ spr_abc = 0.4, # SPR_abc
+ spr_msy = 0.35, # SPR_msy
+ sp_mo = 8, # spawn_month
+ nages = 11, # N_ages
+ Frat = 1, # F_ratio
+ # M
+ M = c(0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3,0.3),
+ # Maturity
+ pmat = c(0.005,0.037,0.224,0.688,0.944,0.992,0.999,1,1,1,1),
+ # Wt_at_age_spawners
+ wtage_sp = c(44.8,161.377,398.272,557.695,652.113,719.573,863.744,948.744,921.397,885.912,1069.87),
+ # Wt_at_age_fsh
+ wtage_fsh = c(69.3778,253.522,408.211,614.731,668.483,718.137,803.017,798.707,788.117,842.468,960.006),
+ # select
+ sel = c(0.002576427,0.040030753,0.651104228,0.768404263,0.794886081,1,0.889293108,0.604815671,0.451169778,0.403195516,0.403195516),
+ # N
+ N = c(511.179,378.528,278.443,194.385,183.423,45.5404,61.1188,19.8073,39.6285,33.6501,51.7806),
+ # Nyrs
+ nyrs = 37,
+ # recruits
+ R = c(1578.51,479.509,357.919,443.588,318.981,413.125,514.351,600.987,536.301,692.34,452.279,1618.73,702.801,372.811,597.812,1136.19,402.862,424.179,1025.36,207.05,383.695,1054.64,2224.52,1379.34,1545.81,345.556,454.884,617.194,404.876,993.497,727.754,236.795,505.948,258.653,726.627,524.198,473.54),
+ # SSB
+ SSB = c(206.391,194.569,187.097,183.296,195.289,240.774,252.143,238.163,223.199,200.776,182.378,179.671,189.193,198.242,212.123,233.484,277.891,282.004,250.709,231.816,218.403,195.275,181.628,189.976,175.912,168.624,220.206,315.124,376.621,397.171,365.476,317.159,277.777,242.719,233.415,223.636,198.117,183.537,177.91)
+ )
+```
+
+## Save lists for running model to files expected by projection model
+
+
+```r
+# Setup.dat
+list2dat(setup,"setup.dat")
+# spp_catch.dat
+list2dat(config,"data/t1_spcat.dat")
+runfn<-"t1"
+file.copy(paste0("data/",runfn,"_spcat.dat"),"spp_catch.dat",overwrite=TRUE)
+```
+
+```
+## [1] TRUE
+```
+
+```r
+list2dat(datfile,"data/t1.dat")
+```
+
+## Run projection model
+
+
+```r
+system("../src/main")
+```
+
+## Read in projection model mainfiles
+
+
+```r
+ .projdir="t1/"
+ dir.create(.projdir)
+```
+
+```
+## Warning in dir.create(.projdir): 't1' already exists
+```
+
+```r
+ file.copy(list.files(getwd(), pattern="out$"), .projdir,overwrite=TRUE)
+```
+
+```
+## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
+```
+
+```r
+ file.remove(list.files(getwd(), pattern="out$"))
+```
+
+```
+## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
+```
+
+```r
+ bf <- data.frame(read.table(paste0(.projdir,"bigfile.out"),header=TRUE,as.is=TRUE))
+ bfs <- bf %>% filter(Sim<=30)
+ #write.csv(bfs,"data/proj.csv")
+ # head(bfs)
+ bfss <- bfs %>% filter(Alt==2) %>% select(Alt,Yr,Catch,SSB,Sim)
+ pf <- data.frame(read.table(paste0(.projdir,"percentdb.out"),header=F) )
+ names(pf) <- c("stock","Alt","Yr","variable","value")
+```
+
+## Make plot of projection model simulations
+
+
+```r
+ p1 <- pf %>% filter(substr(variable,1,1)=="C",variable!="CStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=CMean),width=1.2) + geom_ribbon(aes(ymax=CUCI,ymin=CLCI),fill="goldenrod",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 ABC (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=Cabc)) + geom_line(aes(y=Cofl),linetype="dashed") + geom_line(data=bfss,aes(x=Yr,y=Catch,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ p2 <- pf %>% filter(substr(variable,1,1)=="S",variable!="SSBStdn",Alt==2) %>% select(Yr,variable,value) %>% spread(variable,value) %>%
+ ggplot(aes(x=Yr,y=SSBMean),width=1.2) + geom_ribbon(aes(ymax=SSBUCI,ymin=SSBLCI),fill="coral",alpha=.5) + theme_few() + geom_line() +
+ scale_x_continuous(breaks=seq(thisyr,thisyr+14,2)) + xlab("Year") + ylab("Tier 3 Spawning biomass (kt)") + geom_point() +
+ expand_limits(y=0) +
+ geom_line(aes(y=SSBFabc)) + geom_line(aes(y=SSBFofl),linetype="dashed")+ geom_line(data=bfss,aes(x=Yr,y=SSB,col=as.factor(Sim)))+ guides(size=FALSE,fill=FALSE,alpha=FALSE,col=FALSE)
+ t3 <- grid.arrange(p1, p2, nrow=2)
+```
+
+
+
+```r
+ ggsave(paste0(.projdir,"tier3_proj.pdf"),plot=t3,width=5.4,height=7,units="in")
+```
+
+## Make tables
+
+
+```r
+ # Stock Alt Sim Yr SSB Rec Tot_biom SPR_Implied F Ntot Catch ABC OFL AvgAge AvgAgeTot SexRatio FABC FOFL
+ bfsum <- bf %>% select(Alt,Yr,SSB,F,ABC ,Catch) %>% group_by(Alt,Yr) %>% summarise(Catch=mean(Catch),SSB=mean(SSB),F=mean(F),ABC=mean(ABC))
+ t1 <- bfsum %>% select(Alt,Yr,Catch) %>% spread(Alt,Catch)
+ names(t1) <- c("Catch","Scenario 1","Scenario 2","Scenario 3","Scenario 4","Scenario 5","Scenario 6","Scenario 7")
+
+ print_Tier3_tables(bf)
+```
+
+```
+##
+##
+##
+## Tier 3 projections of BSAI Atka mackerel catch for the 7 scenarios.
+## | Catch | Scenario.1 | Scenario.2 | Scenario.3 | Scenario.4 | Scenario.5 | Scenario.6 | Scenario.7 |
+## | 2019 | 55 | 55 | 55 | 55 | 55 | 55 | 55 |
+## | 2020 | 91 | 55 | 21 | 27 | 0 | 106 | 91 |
+## | 2021 | 83 | 85 | 22 | 28 | 0 | 91 | 83 |
+## | 2022 | 81 | 80 | 25 | 30 | 0 | 85 | 93 |
+## | 2023 | 82 | 81 | 27 | 33 | 0 | 88 | 91 |
+## | 2024 | 86 | 84 | 29 | 35 | 0 | 92 | 93 |
+## | 2025 | 90 | 88 | 30 | 37 | 0 | 96 | 97 |
+## | 2026 | 90 | 88 | 31 | 38 | 0 | 96 | 97 |
+## | 2027 | 89 | 86 | 31 | 38 | 0 | 94 | 94 |
+## | 2028 | 89 | 86 | 31 | 38 | 0 | 95 | 95 |
+## | 2029 | 85 | 82 | 31 | 37 | 0 | 90 | 90 |
+## | 2030 | 84 | 80 | 30 | 37 | 0 | 88 | 88 |
+## | 2031 | 85 | 83 | 31 | 37 | 0 | 91 | 91 |
+## | 2032 | 88 | 85 | 31 | 38 | 0 | 93 | 93 |
+##
+##
+##
+##
+## Tier 3 projections of BSAI Atka mackerel ABC for the 7 scenarios.
+## | SSB | Scenario.1 | Scenario.2 | Scenario.3 | Scenario.4 | Scenario.5 | Scenario.6 | Scenario.7 |
+## | 2019 | 172 | 172 | 172 | 172 | 172 | 172 | 172 |
+## | 2020 | 156 | 164 | 172 | 171 | 177 | 152 | 156 |
+## | 2021 | 138 | 152 | 178 | 175 | 192 | 131 | 138 |
+## | 2022 | 130 | 141 | 187 | 182 | 209 | 122 | 128 |
+## | 2023 | 132 | 142 | 204 | 197 | 233 | 123 | 126 |
+## | 2024 | 135 | 144 | 218 | 210 | 254 | 125 | 127 |
+## | 2025 | 139 | 146 | 232 | 223 | 275 | 128 | 129 |
+## | 2026 | 140 | 147 | 243 | 232 | 292 | 129 | 129 |
+## | 2027 | 139 | 146 | 249 | 237 | 302 | 127 | 127 |
+## | 2028 | 138 | 145 | 252 | 240 | 309 | 126 | 126 |
+## | 2029 | 137 | 144 | 254 | 241 | 314 | 125 | 125 |
+## | 2030 | 134 | 141 | 251 | 238 | 313 | 122 | 122 |
+## | 2031 | 134 | 141 | 252 | 238 | 315 | 123 | 123 |
+## | 2032 | 137 | 144 | 255 | 242 | 319 | 125 | 125 |
+##
+##
+##
+##
+## Tier 3 projections of BSAI Atka mackerel fishing mortality for the 7 scenarios.
+## | F | Scenario.1 | Scenario.2 | Scenario.3 | Scenario.4 | Scenario.5 | Scenario.6 | Scenario.7 |
+## | 2019 | 0.169 | 0.169 | 0.169 | 0.169 | 0.169 | 0.169 | 0.169 |
+## | 2020 | 0.299 | 0.172 | 0.066 | 0.083 | 0.000 | 0.353 | 0.299 |
+## | 2021 | 0.299 | 0.283 | 0.066 | 0.083 | 0.000 | 0.343 | 0.299 |
+## | 2022 | 0.284 | 0.263 | 0.066 | 0.083 | 0.000 | 0.315 | 0.330 |
+## | 2023 | 0.277 | 0.259 | 0.066 | 0.083 | 0.000 | 0.311 | 0.317 |
+## | 2024 | 0.275 | 0.258 | 0.066 | 0.083 | 0.000 | 0.310 | 0.312 |
+## | 2025 | 0.276 | 0.260 | 0.066 | 0.083 | 0.000 | 0.312 | 0.314 |
+## | 2026 | 0.277 | 0.261 | 0.066 | 0.083 | 0.000 | 0.313 | 0.314 |
+## | 2027 | 0.277 | 0.260 | 0.066 | 0.083 | 0.000 | 0.313 | 0.313 |
+## | 2028 | 0.278 | 0.259 | 0.066 | 0.083 | 0.000 | 0.312 | 0.312 |
+## | 2029 | 0.277 | 0.258 | 0.066 | 0.083 | 0.000 | 0.310 | 0.310 |
+## | 2030 | 0.275 | 0.254 | 0.066 | 0.083 | 0.000 | 0.306 | 0.306 |
+## | 2031 | 0.274 | 0.254 | 0.066 | 0.083 | 0.000 | 0.307 | 0.307 |
+## | 2032 | 0.275 | 0.256 | 0.066 | 0.083 | 0.000 | 0.309 | 0.309 |
+##
+##
+##
+##
+## Tier 3 projections of BSAI Atka mackerel spawning biomass for the 7 scenarios.
+## | ABC | Scenario.1 | Scenario.2 | Scenario.3 | Scenario.4 | Scenario.5 | Scenario.6 | Scenario.7 |
+## | 2019 | 93 | 93 | 22 | 27 | 0 | 107 | 107 |
+## | 2020 | 91 | 89 | 21 | 27 | 0 | 106 | 106 |
+## | 2021 | 83 | 85 | 22 | 28 | 0 | 91 | 96 |
+## | 2022 | 81 | 80 | 25 | 30 | 0 | 85 | 93 |
+## | 2023 | 82 | 81 | 27 | 33 | 0 | 88 | 91 |
+## | 2024 | 86 | 84 | 29 | 35 | 0 | 92 | 93 |
+## | 2025 | 90 | 88 | 30 | 37 | 0 | 96 | 97 |
+## | 2026 | 90 | 88 | 31 | 38 | 0 | 96 | 97 |
+## | 2027 | 89 | 86 | 31 | 38 | 0 | 94 | 94 |
+## | 2028 | 89 | 86 | 31 | 38 | 0 | 95 | 95 |
+## | 2029 | 85 | 82 | 31 | 37 | 0 | 90 | 90 |
+## | 2030 | 84 | 80 | 30 | 37 | 0 | 88 | 88 |
+## | 2031 | 85 | 83 | 31 | 37 | 0 | 91 | 91 |
+## | 2032 | 88 | 85 | 31 | 38 | 0 | 93 | 93 |
+##
+```
+
diff --git a/examples/Misc/yfs.dat b/examples/Misc/yfs.dat
new file mode 100644
index 0000000..49f8148
--- /dev/null
+++ b/examples/Misc/yfs.dat
@@ -0,0 +1,41 @@
+BSAI_YFINSS
+0 # SSL Species???
+0 # Constant buffer of Dorn?
+1 # Number of fisheries
+2 # Number of sexes??
+0.0522 # average 5yr f
+1.0 # author f
+0.4 # SPR ABC
+0.35 # SPR MSY
+5 # Spawnmo
+20 # Number of ages
+1 # Fratio
+ # 0.5 0.5 # Fratio
+#females first
+0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
+#males
+0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12 0.12
+# Maturity Females first row used in projections (average of 2012 and 1992 studies)
+0.0001 0.0003 0.0008 0.0018 0.0042 0.011 0.0382 0.097 0.2134 0.4134 0.6459 0.825 0.9195 0.9675 0.9855 0.996 0.9985 1 1 1
+# Maturity Males
+0.00026 0.0006 0.00159 0.0035 0.0084 0.012 0.0463 0.1039 0.2167 0.3967 0.6118 0.7899 0.899 0.955 0.981 0.992 0.997 1 1 1
+# Wt spawn females
+0.006 0.008 0.011 0.016 0.039 0.053 0.122 0.21 0.273 0.36 0.399 0.42 0.434 0.453 0.471 0.508 0.529 0.537 0.546 0.59
+# WtAge Females, by fishery
+0.006 0.008 0.032 0.05 0.066 0.074 0.112 0.186 0.338 0.372 0.412 0.408 0.455 0.456 0.485 0.508 0.515 0.532 0.555 0.59
+# WtAge Males, by fishery
+0.004 0.008 0.043 0.057 0.063 0.082 0.116 0.171 0.253 0.319 0.318 0.331 0.338 0.346 0.381 0.383 0.408 0.434 0.413 0.46
+# Selectivity Females, by fishery using last year in time-varying model
+0.000701093 0.00157433 0.00353139 0.00790201 0.0175864 0.0386771 0.0829259 0.168903 0.313545 0.506555 0.697632 0.838332 0.920976 0.963226 0.983297 0.992499 0.996648 0.996648 0.996648 0.996648
+# Selectivity males, by fishery using last year in time-varying model
+0.000663833 0.0014836 0.00331236 0.00737866 0.0163549 0.0358564 0.076796 0.156873 0.293871 0.482098 0.675547 0.823232 0.91241 0.958847 0.981173 0.991494 0.996179 0.996179 0.996179 0.996179
+# N at age in 2018 Females, Males
+1179.8 1033.99 934.165 1186.63 1173.63 564.941 178.599 250.413 660.076 421.127 310.649 287.582 182.132 119.927 221.512 114.413 97.0637 58.0578 75.7241 362.037
+1179.8 1033.94 934.086 1186.45 1173.45 565.078 178.989 252.484 673.115 433.036 318.252 291.14 182.187 118.758 217.846 112.506 96.4746 58.1467 75.3316 360.317
+# Number of Recruits
+33
+# Recruitment 1978-2010
+2467.15 1574.561 3038.13 2258.42 6528.73 1209.186 5401.76 1868.121 1435.065 1959.887 2685 2682.95 1339.353 1503.374 3331.54 1988.833 1679.79 1692.189 4175.32 1805.508 1499.919 1839.244 2606.42 1658.825 2284.2 2219.71 3523.92 1533.483 1844.777 2349.6 2112.732 2434.22 3272.3
+# SSB
+# used only for S/R analysis
+293.496 409.183 526.225 656.554 776.716 843.983 945.256 1024.84 1071.79 1058.17 1051.2 993.336 965.759 976.509 1057.25 1139.5 1174.45 1177.44 1176.38 1110.03 1073.52 1009.26 1000.98 987.701 983.504 981.707 990.525 1021.07 1036.84 1056.94 1061.63 1036.55 1000.41 975.522 952.413 934.539 924.131 881.156 872.721 874.66 856.504 854.096
diff --git a/examples/atka/plot.R b/examples/atka/plot.R
new file mode 100644
index 0000000..8555737
--- /dev/null
+++ b/examples/atka/plot.R
@@ -0,0 +1,39 @@
+radian
+rm(list=ls())
+ls()
+#source(paste0(here::here("R"),"/prelims.R"))
+source("../../R/prelims.R")
+#-------------------------------------------------------------------------------
+# Visual compare runs
+#-------------------------------------------------------------------------------
+library(ggridges)
+
+source("../../R/compareRuns.r")
+#--Projections---------
+pdf("proj.pdf")
+pfn <- "amak"
+pfn <- "5yr"
+pdt <- read_csv(paste0("spm_detail.csv"))
+pdt$Alternative <- as.factor(pdt$Alternative)
+names(pdt)
+
+pt <- pdt%>%filter(Yr>2021,Alternative==2) %>% group_by(Yr,Alternative) %>% summarise(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL),SSB=median(SSB) ,lb=quantile(SSB,.2) ,ub=quantile(SSB,.8) )
+pt
+ggplot(pt,aes(x=Yr,y=SSB,fill=Alternative)) + geom_line() + mytheme + ylim(c(0,340000)) + geom_ribbon(aes(ymin=lb,ymax=ub,fill=Alternative),alpha=0.25) + labs(y="Spawning biomass (kt)",x="Year") + scale_x_continuous(breaks=seq(2022,2035,2))
+c1 <- ggplot(pt,aes(x=Yr,y=Catch,color=Alternative,size=1.)) + geom_line(size=1.5) + mytheme + labs(y="Catch (kt)",x="Year") + scale_x_continuous(breaks=seq(2015,2032,2))
+c1 <- c1 + geom_line(aes(x=Yr,y=ABC),size=1)
+#c1 <- c1 + geom_line(data=pt[as.numeric(Alternative)==2,.(Yr,ABC)],aes(x=Yr,y=ABC))
+c1
+pt[as.numeric(Alternative)==2,.(Yr,ABC)]
+pt <- pdt[Yr>2018,.(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL)),.(Yr,Alternative)]
+pt
+ggplot(pt,aes(x=Yr,y=OFL,color=Alternative)) + geom_line() + mytheme
+pdt
+pdx <-rbind(pdt)
+setkey(pdx,Yr,Alternative)
+
+pt <- pdx[.(Yr>2018,(Alternative)==1),.(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL),SSB=median(SSB) ),.(Yr,config)]
+pt <- pdx[Yr>2018&Alternative==1,.(Catch=mean(Catch),ABC=mean(ABC),OFL=mean(OFL),SSB=median(SSB) ,lb=quantile(Catch,.1) ,ub=quantile(Catch,.9) ),.(Yr,config)]
+ggplot(pt,aes(x=Yr,y=Catch,color=config,shape=config)) + geom_line() + geom_point() + mytheme + geom_ribbon(aes(ymin=lb,ymax=ub,fill=config),alpha=0.25) + labs(x="Year")+ scale_x_continuous(breaks=seq(2015,2032,2)) + ylim(c(0,130))
+dev.off()
+
diff --git a/examples/atka/readme.txt b/examples/atka/readme.txt
new file mode 100644
index 0000000..0d84398
--- /dev/null
+++ b/examples/atka/readme.txt
@@ -0,0 +1 @@
+added ABC output for Tier 1 under any scenario...
diff --git a/examples/atka/setup.dat b/examples/atka/setup.dat
new file mode 100644
index 0000000..2311eae
--- /dev/null
+++ b/examples/atka/setup.dat
@@ -0,0 +1,25 @@
+#--------------------------------------------
+Std # Run name 5 scenarios
+#--------------------------------------------
+7 # Number of Alternatives
+#--------------------------------------------
+1
+2 # List of alternatives
+3
+4
+5
+6
+7
+#--------------------------------------------
+1 # Flag to set TAC equal to ABC (1 means true, otherwise false)
+#--------------------------------------------
+1 # Stock-recruitment type (1=Ricker, 2=Bholt)
+1 # projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+0 # SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+0.0 # Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+#--------------------------------------------
+1 # Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+#--------------------------------------------
+14 #_Number of projection years
+1000 #_Number of simulations
+2021 #_Begin Year
diff --git a/examples/atka/spm.dat b/examples/atka/spm.dat
new file mode 100644
index 0000000..d350203
--- /dev/null
+++ b/examples/atka/spm.dat
@@ -0,0 +1,53 @@
+#_SETUP_FILE_FOR_skates
+std # Run name 5 scenarios
+# Tier
+3
+#--------------------------------------------
+7 # Number of Alternatives
+#--------------------------------------------
+1
+2 # List of alternatives
+3
+4
+5
+6
+7
+#--------------------------------------------
+1 # Flag to set TAC equal to ABC (1 means true, otherwise false)
+#--------------------------------------------
+2 # Stock-recruitment type (1=Ricker, 2=Bholt)
+1 # projection rec form (default: 1 = use observed mean and std, option 2 = use estimated SRR and estimated sigma R)
+0 # SR-Conditioning (0 means no, 1 means use Fmsy == F35%?, 2 means Fmsy == F35% and Bmsy=B35% condition (affects SRR fits)
+0 # Condition that there is a prior that mean historical recruitment is similar to expected recruitment at half mean SSB and double mean SSB 0 means don't use, otherwise specify CV
+#--------------------------------------------
+1 # Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it) Write_Big
+#--------------------------------------------
+15 #_Number of projection years
+1000 #_Number of simulations
+2022 #_Begin Year
+#_Number_of_years with specified catch
+3
+# Number of species
+1
+# OY Min
+0
+# OY Max
+2.00E+06
+# data files for each species
+../amak.prj
+# ABC Multipliers
+1
+# scalars
+1
+# New Alt 4 Fabc SPRs (Rockfish = 0.75, other 0.6), Steller sea lion prey species between F40 and F60 (to meet OY Min)
+0.75
+# Number of TAC model categories
+1
+# TAC model indices (for aggregating)
+1
+# Catch in each future year
+# Catch in each future year
+2022 66481
+2023 83800
+2024 73495
+
diff --git a/examples/atka/spm.tpl b/examples/atka/spm.tpl
new file mode 100644
index 0000000..b791fc1
--- /dev/null
+++ b/examples/atka/spm.tpl
@@ -0,0 +1,2509 @@
+ // DONT FORGET ABOU RHO!!!
+ ////////////////////////////////////////////////////////////////////////////
+ // spm.tpl
+ // Version 0.6 Released Oct 2022
+ //
+ // NOTE: Draft subject to change
+ // Conventions:
+ // nspp, ispp = index for species/stock
+ // ngear, igear= index for gears used within age-structured species/stocks
+ //
+ // Oct 2005 add ability to compute objective function for maximizing yield while minimizing variability
+ // want to have the ability to compute the average first-diff squared and avg yield for each trajectory
+ // Alt3b_obj_fun = avg_yield - 0.5 * sqrd_1st_diff
+ //
+ ////////////////////////////////////////////////////////////////////////////
+DATA_SECTION
+ !!CLASS ofstream means_out("means.out")
+ !!CLASS ofstream alts_proj("alt_proj.out")
+ !!CLASS ofstream percent_out("percentiles.out")
+ !!CLASS ofstream percent_db("percentdb.out")
+ // !!CLASS ofstream Alt3bstuff("alt3b.out")
+ !!CLASS ofstream detail_out("spm_detail.csv")
+ !!CLASS ofstream prof_F("F_profile.out");
+ !!CLASS ofstream elasticity("elasticity.csv");
+ int condition_SR
+ int ipro
+ int isim
+ number kink_adj
+ !! condition_SR = 0;//
+ int rnseed // Random number seed
+ !! rnseed = 123;
+ ////////////////////////////////////////////////////////////////////////////
+ // !! ad_comm::change_datafile_name("setup.dat"); Jim changed to spm.dat file
+ !! *(ad_comm::global_datafile) >> run_name; // Read in the name of this run
+ init_int Tier
+ int alt ; // Alternative specfications (relic of PSEIS)
+ vector rec_vector(1,300) // storage vector for historical and future recruitment...
+ vector wtd_rec(1,25)
+ vector wtd_div(1,25) // denominator of wtd recruitment
+ init_int nalts
+ init_ivector alt_list(1,nalts)
+ init_int TAC_ABC // Flag to set TAC equal to ABC (1 means true, otherwise false)
+ init_int SrType // Type of recruitment curve (1=Ricker, 2 Bholt)
+ // Specify form of recruitment generator (1 = use observed mean and std
+ // 2 = use estimated SRR and estimated sigma R
+ // 3 = use input parameters from srecpar.dat :
+ // 4 = use input parameters from srecpar.dat :
+ init_int Rec_Gen
+ init_int Fmsy_F35 // Specify if conditioned so that Fmsy = F35 (affects SRR fitting) may need to be species specific...
+ init_number Rec_Cond // Specify prior condition that recruitment at half and double average SSB is similar to average historical Rec
+ init_int Write_Big // Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it)
+ init_int npro // Number of projection years
+ init_int nsims // Number of simulaions
+ init_int styr // First year of projection
+ !! cout<< "First year:\t"<> spp_file_name(i);
+
+ write_log(nyrs_catch_in);
+ write_log(nspp);
+ write_log(OY_min);
+ write_log(OY_max);
+ write_log(spp_file_name);
+
+ END_CALCS
+ !! cout <<"OYMax "<> bzero_in;
+ *(ad_comm::global_datafile) >> phizero_in;
+ *(ad_comm::global_datafile) >> alpha_in;
+ *(ad_comm::global_datafile) >> sigmar_in;
+ *(ad_comm::global_datafile) >> rho_in;
+ // rho_in=0.82;
+ }
+ // Open up tac-model parameters
+ ad_comm::change_datafile_name("tacpar.dat");
+ END_CALCS
+ init_int nntmp
+ init_int nnodes
+ init_vector maxabc(1,ntacspp)
+ init_matrix theta(0,nnodes,1,ntacspp)
+ !! cout<<"read tacpar"<> spname(i); // 1
+ cout<> SSL_spp(i); // 2
+ *(ad_comm::global_datafile) >> Const_Buffer(i); // 3
+ *(ad_comm::global_datafile) >> ngear(i); // 4
+ *(ad_comm::global_datafile) >> nsexes(i); // 6
+ *(ad_comm::global_datafile) >> avg_5yrF(i); // 7
+ *(ad_comm::global_datafile) >> FABC_Adj(i); // 8
+ *(ad_comm::global_datafile) >> SPR_abc(i); // 9
+ *(ad_comm::global_datafile) >> SPR_ofl(i); // 10
+ *(ad_comm::global_datafile) >> spawnmo(i); // 11
+ *(ad_comm::global_datafile) >> nages(i); // 12
+ cout<<"nages: "<> Fratiotmp(i,j); // 13
+ write_log( Fratiotmp(i));
+
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> M_Ftmp(i,k); // 14
+ if (nsexes(i)==2)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> M_Mtmp(i,k); // 15
+ write_log( M_Ftmp(i));
+
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> pmaturetmp_F(i,k); // 16
+ write_log( pmaturetmp_F(i));
+
+ if (nsexes(i)==2)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> pmaturetmp_M(i,k); // 15
+ write_log( pmaturetmp_M(i));
+ cout << "Mature: "<< pmaturetmp_F(i)(1,nages(i)) <> wt_Ftmp(i,k); // 17
+ write_log( wt_Ftmp(i));
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> wt_gear_Ftmp(i,j,k); // 18
+ write_log( wt_gear_Ftmp(i));
+
+ if (nsexes(i)==2)
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> wt_gear_Mtmp(i,j,k);// 19
+ write_log( wt_gear_Mtmp(i));
+
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> sel_Ftmp(i,j,k); // 20
+ write_log( sel_Ftmp(i));
+
+ if (nsexes(i)==2)
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> sel_Mtmp(i,j,k); // 21
+ write_log( sel_Mtmp(i));
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> n0_Ftmp(i,k); // 22
+ write_log( n0_Ftmp(i));
+
+ if (nsexes(i)==2)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> n0_Mtmp(i,k); // 23
+ write_log( n0_Mtmp(i));
+
+ cout<<"N: "<> nrec(i); // 24
+ cout<<"nrec: "<> Rtmp(i,j); // 25
+ cout<<"Rec: "<> SSBtmp(i,j); // 26
+ cout<<"SSB: "<> Rsim(ispp,i,j);
+ Rsim(ispp,i,j) *= exp(rnorms(i,j)*.375); // about 15% CV to get historical mean
+ }
+ if (nsims<=5) Rsim(ispp,i,j) = AMeanRec(ispp); // XXX constant recruitment
+ }
+ }
+ envin.close();
+ AMeanRec(ispp) *= .5; // Arithmetic mean
+ cout <<"recruits"< 0 )
+ F_begin_yr(k,ispp) = SolveF2(n0_F(ispp), n0_M(ispp), Obs_Catch(k,ispp),ispp); // F_yr_one(ispp) = SolveF2(n0_F(ispp), n0_M(ispp), yr_one_catch(ispp),ispp);
+ else
+ F_begin_yr(ispp) = 0; // F_yr_one(ispp) = 0; // cout< 1) // Debugging if LP to be done or not
+ {
+ int on=0;
+ if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nolp"))>-1)
+ dolp=0;
+ else
+ dolp=1;
+ }
+ END_CALCS
+
+
+PRELIMINARY_CALCS_SECTION
+ double tmp1;
+ write_alts_hdr();
+ get_SB100();
+ cout<<"SSB, Biomass at unfished: "<1)
+ nyrs_catch = nyrs_catch_in;
+ else
+ nyrs_catch = 1;// NOTA BUENO: this is changed under the new EIS Alternatives (May 06)
+ Do_Sims();
+ // if (alt==2) write_alts();
+ write_alts();
+ }
+
+FUNCTION compute_obj_fun
+ dvariable tmp1;
+ dvariable tmp2;
+ tmp1.initialize();
+ for (int ispp=1;ispp<=nspp;ispp++)
+ {
+ Get_Bzero(ispp);
+ // write_srec();exit(1);
+ if(Fmsy_F35>0)
+ {
+ get_msy(ispp);
+ switch (current_phase())
+ {
+ case 1 : tmp1 = 1.e1*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ case 2 : tmp1 = 1.e2*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ case 3 : tmp1 = 1.e2*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ case 4 : tmp1 = 1.e3*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ default: tmp1 = 1.e3*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ }
+ if(Fmsy_F35==2)
+ switch (current_phase())
+ {
+ case 1 : tmp1 += 1.e1*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ case 2 : tmp1 += 1.e2*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ case 3 : tmp1 += 1.e2*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ case 4 : tmp1 += 1.e3*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ default : tmp1+= 1.e3*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ }
+ obj_fun += tmp1;
+ }
+ if (Rec_Cond>0.)
+ {
+ // More conditioning here--make mean recruitment consistent with half and double
+ // recruitment at avg spawning biomass...
+ double vartmp = 2.*Rec_Cond*Rec_Cond;
+ /*
+ double ssb1 = AMeanSSB(ispp) * 0.5;
+ double ssb2 = AMeanSSB(ispp) ;
+ double ssb3 = AMeanSSB(ispp) * 2.0;
+ // double ssb4 = value(Bzero(ispp)) ;
+ tmp2.initialize();
+ tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb1, ispp)))/ vartmp;
+ tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb2, ispp)))/ vartmp;
+ tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb3, ispp)))/ vartmp;
+ // tmp2 += 24.* square(log(AMeanRec(ispp)) - log(SRecruit( ssb4, ispp)));
+ // tmp2 += 24.* square(log(AMeanSSB(ispp)) - log( ssb4 ));
+ */
+ tmp2 = square(log(Bzero(ispp)) - log(SB100(ispp)))/ vartmp;
+ obj_fun += tmp2;
+ // cout << "SSB "<OY_max)
+ Actual_Catch = TAC/sum(TAC)*OY_max;
+ else
+ Actual_Catch = TAC;
+ }
+ /* ////////////////////////////////////////////////////////////////////////////// //int use_max = 1; // Need to move this into the setup file... double max_catch=1500.;// EBS POllock special case (CHANGE THIS) //if (use_max==1) for (int ispp=1;ispp<=nspp;ispp++) Actual_Catch(ispp) = min(TAC(ispp),max_catch); // cout<<"AC: "<2)
+ write_sim("Alternative ",ispp);
+ }
+ write_by_time();
+
+
+FUNCTION Alt4_TAC
+ // Return vector of Alt4 TACs (given Alt4_Fabc, and the SSL prey condition to be at least up to SSL MaxPerm)
+ // Get all ABCs under Alt4_Fabc
+ // Sum and test if below OY_min
+ double sumabc;
+ sumabc = sum(ABC);
+ if (sumabc3)
+ {
+ Alt4_Fcalc = 1;
+ double diff;
+ double ssl_sum=0.;
+ diff = OY_min-sumabc;
+ // if so, sum up SSL ABCs
+ ssl_sum = SSL_spp*ABC;
+ for (int ispp=1;ispp<=nspp;ispp++)
+ {
+ // then apportion them to difference between current ABC and target (OY_min) difference
+ if (SSL_spp(ispp))
+ {
+ TAC(ispp) = ABC(ispp) + diff * ABC(ispp)/ssl_sum;
+ // Special for BSAI Pcod to subtract off 3%...and get the totals to match up
+ // if (ispp==2) TAC(ispp) /= 0.97; For EIS work
+ }
+ else
+ {
+ TAC(ispp) = ABC(ispp);
+ }
+ }
+ // cout << ipro<<" "<0.01) ftmp = SolveF2(N_F(ispp),N_M(ispp),TAC(ispp),ispp);
+ // cout <= 0.2 *SBzero(ispp) & SBtmp < SBKink(ispp) ) // NOTE Same as Am 56 here (until 20% bzero reached)
+ Ftmp = (Fabc(ispp)*(1/(1-alpha ))*(SBtmp/SBKink(ispp) - alpha ));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fabc(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ return(Ftmp);
+
+FUNCTION double Get_F_Am56(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ {
+ for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct
+ {
+ if (SBtmp < alpha*SBKink(ispp))
+ Ftmp = 0.;
+ if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) )
+ Ftmp = (Fabc(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fabc(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ }
+ return(Ftmp);
+
+FUNCTION double Get_Fofl_t2(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ {
+ for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct
+ {
+ if (SBtmp < alpha*SBKink(ispp))
+ Ftmp = 0.;
+ if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) )
+ Ftmp = (Fofl(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fofl(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ }
+ return(Ftmp);
+
+FUNCTION double Get_Fofl_t(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ {
+ for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct
+ {
+ if (SBtmp < alpha*SBKink(ispp))
+ Ftmp = 0.;
+ if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) )
+ Ftmp = (Fofl(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fofl(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ }
+ return(Ftmp);
+
+
+FUNCTION double Get_F_t(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ if (SSL_spp(ispp))
+ {
+ Ftmp = Get_F_SSL_prey(F_age, N_females, ispp);
+ }
+ else
+ {
+ Ftmp = Get_F_Am56(F_age, N_females, ispp);
+ }
+ return(Ftmp);
+
+FUNCTION void Project_Pops(const int& isim, const int& i)
+ double ctmp;
+ if ( i == npro && isim%int(nsims/4)==0) cout << "Year "< 0)
+ {
+ if (i <= nyrs_catch || TAC_ABC==0) // Use TAC setting algorithm for alt 2 only, for all others, set TAC==ABC
+ {
+ Ftmp = SolveF2(N_F(ispp),N_M(ispp),Actual_Catch(ispp),ispp);
+ }
+ else
+ {
+ if (alt==7 && i<=3)
+ Ftmp = Get_F(1,ispp); // Set to the F rather than solving every time...
+ else
+ if (alt==77 && i<=2)
+ Ftmp = Get_F(1,ispp); // Set to the F rather than solving every time...
+ else
+ {
+ // if (alt==2 && sum(TAC)>OY_max)
+ if ((alt==2 ||alt==98||alt==97) && sum(TAC)>OY_max)
+ {
+ Ftmp = SolveF2(N_F(ispp),N_M(ispp),OY_max*TAC(ispp)/sum(TAC),ispp);
+ // cout<<"HERE "<< Ftmp<<" "< 0.)
+ {
+ for (m=1;m<=ngear(ispp);m++)
+ {
+ Ftottmp_F = Ftmp*Frat(ispp,m)*sel_F(ispp,m);
+ Ftottmp_spr += Ftottmp_F ;
+ Ftottmp_M = Ftmp*Frat(ispp,m)*sel_M(ispp,m);
+ ctmp += (wt_gear_F(ispp,m) * elem_prod(elem_div(Ftottmp_F, Z_F(ispp)),elem_prod(1.-S_F(ispp),N_F(ispp)))); // Catch equation (vectors)
+ ctmp += (wt_gear_M(ispp,m) * elem_prod(elem_div(Ftottmp_M, Z_M(ispp)),elem_prod(1.-S_M(ispp),N_M(ispp)))); // Catch equation (vectors)
+ }
+ SPRsim(ispp,isim,i) = Implied_SPR(Ftottmp_spr,ispp);
+ // cout <NUL ");
+ ifstream tac("tac.dat");
+ tac>>TAC;
+ tac.close(); /* */
+ for (int ispp=1;ispp<=nspp;ispp++)
+ Actual_Catch(ispp) = min(TAC(ispp),ABC(ispp));
+ }
+
+
+FUNCTION Avg_Age
+ Avg_Age_End.initialize();
+ Avg_Age_Mat.initialize();
+ for (int ispp=1;ispp<=nspp;ispp++)
+ {
+ // if(!isit_const(ispp) )
+ {
+ dvector age_seq(1,nages(ispp));
+ age_seq.initialize();
+ age_seq.fill_seqadd(1,1);
+ dvector ntmp(1,nages(ispp));
+ ntmp = N_F(ispp);
+ Avg_Age_End(ispp) += age_seq * ntmp /sum(ntmp);
+ Avg_Age_sum(ispp) += Avg_Age_End(ispp) ;
+ ntmp = elem_prod(ntmp,pmature_F(ispp)) ;
+ Avg_Age_Mat(ispp) += elem_prod(ntmp,pmature_F(ispp)) * age_seq/sum(ntmp);
+ // cout< 1e-6)
+ {
+ iter++;
+ ftmp += (TACin-cc) / btmp;
+ Ftottmp_F.initialize();
+ Ftottmp_M.initialize();
+ for (m=1;m<=ngear(ispp);m++)
+ {
+ Ftottmp_F += ftmp*Fratsel_F(m);
+ Ftottmp_M += ftmp*Fratsel_M(m);
+ }
+ Z_F = Ftottmp_F + M_F(ispp);
+ Z_M = Ftottmp_M + M_M(ispp);
+ S_F = mfexp( -Z_F );
+ S_M = mfexp( -Z_M );
+ cc = 0.0;
+ for (m=1;m<=ngear(ispp);m++)
+ {
+ cc += (wt_gear_F(ispp,m) * elem_prod(elem_div(ftmp*Fratsel_F(m), Z_F),elem_prod(1.-S_F,N_F))); // Catch equation (vectors)
+ cc += (wt_gear_M(ispp,m) * elem_prod(elem_div(ftmp*Fratsel_M(m), Z_M),elem_prod(1.-S_M,N_M))); // Catch equation (vectors)
+ }
+ dd = cc / TACin - 1.;
+ //cout << ispp<<" "<< ftmp << " "<< cc << " "<100) {cerr<<"Bombed on catch solver--check scales for "<nages(ispp))
+ {
+ Tg += double(ii) * wt_mature_F(ispp,nages(ispp)) * ntmp;
+ tmp += wt_mature_F(ispp,nages(ispp)) * ntmp;
+ ntmp *= exp(-M_F(ispp,nages(ispp)));
+ }
+ else
+ {
+ Tg += double(ii) * wt_mature_F(ispp,ii) * ntmp;
+ tmp += wt_mature_F(ispp,ii) * ntmp;
+ ntmp *= exp(-M_F(ispp,ii));
+ }
+ }
+ // Tg /= value(phizero(ispp));
+ Tg /= tmp;
+ report << Tg<< " ";
+ }
+ report << endl;
+ report << "Options SR_Type Project_Recr_Assmption SR_Condition"< 1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+FUNCTION void write_TACs(const adstring& Title)
+ // This one prints out species over time (means only), but without headings
+ means_out <<"Alternative "< 1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+
+
+FUNCTION void write_sim(const adstring& Title, d3_array& Outtmp)
+ // This one prints out species over time (means only), but without headings
+ d3_array mtmp(1,nspp,1,npro,1,nsims);
+ for (int ispp=1;ispp<=nspp;ispp++) for (int i=1;i<=npro;i++) for (int k=1;k<=nsims;k++)
+ mtmp(ispp,i,k)=Outtmp(ispp,k,i);
+ means_out <<"Alternative "< 1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+FUNCTION void write_sim(const adstring& Title,const d3_array& Outtmp,const dvar_vector& bb)
+ // This one prints out species over time (means only)
+ d3_array mtmp(1,nspp,1,npro,1,nsims);
+ for (int ispp=1;ispp<=nspp;ispp++) for (int i=1;i<=npro;i++) for (int k=1;k<=nsims;k++)
+ mtmp(ispp,i,k)=Outtmp(ispp,k,i);
+ means_out <<"Alternative "<1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+FUNCTION void write_sim(const adstring& Title,const int& ispp)
+ // For each species separately
+ percent_out <<"Alternative "<
+ adstring xspname;
+ adstring_array targsppname(1,20);
+ adstring_array spp_file_name(1,20);
+ adstring_array gearname(1,8);
+ adstring_array spname(1,90);
+ adstring_array areaname(1,8);
+ adstring run_name;
+
+ ofstream write_log("Input_Log.rep");
+ #undef write_log
+ #define write_log(object) write_log << #object "\n" << object << endl;
+
+ // A routine to get transpose, sort and return a matrix ----
+ dmatrix TranSort (const dmatrix m1)
+ {
+ RETURN_ARRAYS_INCREMENT();
+ dmatrix vtmp=trans(m1);
+ int npro=m1.colmax();
+ int nsim=m1.rowmax();
+ for (int i=1;i<=npro;i++)
+ vtmp(i) = sort(vtmp(i),nsim);
+
+ RETURN_ARRAYS_DECREMENT();
+ return(vtmp);
+ }
+ // #include
+
+FUNCTION dvariable SRecruit(const dvariable& Stmp,const int& ispp)
+ RETURN_ARRAYS_INCREMENT();
+ dvariable RecTmp;
+ switch (SrType)
+ {
+ case 1:
+ RecTmp = (Stmp / phizero(ispp)) * mfexp( sr_alpha(ispp) * ( 1. - Stmp / Bzero(ispp) )) ; //Ricker form from Dorn
+ break;
+ case 2:
+ RecTmp = Stmp / ( sr_alpha(ispp) + beta(ispp) * Stmp); //Beverton-Holt form
+ break;
+ case 3:
+ // RecTmp = mfexp(log_avgrec); //Avg recruitment
+ break;
+ case 4:
+ RecTmp = Stmp * mfexp( sr_alpha(ispp) - Stmp * beta(ispp)) ; //old Ricker form
+ break;
+ }
+
+ RETURN_ARRAYS_DECREMENT();
+ return RecTmp;
+
+
+FUNCTION dvar_vector SRecruit(const dvar_vector& Stmp,const int& ispp)
+ RETURN_ARRAYS_INCREMENT();
+ dvar_vector RecTmp(Stmp.indexmin(),Stmp.indexmax());
+ switch (SrType)
+ {
+ case 1:
+ RecTmp = elem_prod((Stmp / phizero(ispp)) , mfexp( sr_alpha(ispp) * ( 1. - Stmp / Bzero(ispp) ))) ; //Ricker form from Dorn
+ break;
+ case 2:
+ RecTmp = elem_prod(Stmp , 1. / ( sr_alpha(ispp) + beta(ispp) * Stmp)); //Beverton-Holt form
+ break;
+ case 3:
+ // RecTmp = mfexp(log_avgrec); //Avg recruitment
+ break;
+ case 4:
+ RecTmp = elem_prod(Stmp ,mfexp( sr_alpha(ispp) - Stmp * beta(ispp))); //old Ricker form
+ break;
+ }
+ RETURN_ARRAYS_DECREMENT();
+ //cout <5||F1<0.01))
+ // {
+ // ii=5;
+ // F1=M(ispp,nages(ispp)); // When things bomb (F <0 or F really big then just set it to F35...)
+ // }
+ // else
+ {
+ F2 = F1 + df*.5;
+ F3 = F2 - df; //F1 = double(ii)/100;
+ yld1 = yield(F1,Stmp,Rtmp,ispp); //cout <> run_name; // Read in the name of this run
+ init_int Tier
+ int alt ; // Alternative specfications (relic of PSEIS)
+ vector rec_vector(1,300) // storage vector for historical and future recruitment...
+ vector wtd_rec(1,25)
+ vector wtd_div(1,25) // denominator of wtd recruitment
+ init_int nalts
+ init_ivector alt_list(1,nalts)
+ init_int TAC_ABC // Flag to set TAC equal to ABC (1 means true, otherwise false)
+ init_int SrType // Type of recruitment curve (1=Ricker, 2 Bholt)
+ // Specify form of recruitment generator (1 = use observed mean and std
+ // 2 = use estimated SRR and estimated sigma R
+ // 3 = use input parameters from srecpar.dat :
+ // 4 = use input parameters from srecpar.dat :
+ init_int Rec_Gen
+ init_int Fmsy_F35 // Specify if conditioned so that Fmsy = F35 (affects SRR fitting) may need to be species specific...
+ init_number Rec_Cond // Specify prior condition that recruitment at half and double average SSB is similar to average historical Rec
+ init_int Write_Big // Flag to write big file (of all simulations rather than a summary, 0 means don't do it, otherwise do it)
+ init_int npro // Number of projection years
+ init_int nsims // Number of simulaions
+ init_int styr // First year of projection
+ !! cout<< "First year:\t"<> spp_file_name(i);
+
+ write_log(nyrs_catch_in);
+ write_log(nspp);
+ write_log(OY_min);
+ write_log(OY_max);
+ write_log(spp_file_name);
+
+ END_CALCS
+ !! cout <<"OYMax "<> bzero_in;
+ *(ad_comm::global_datafile) >> phizero_in;
+ *(ad_comm::global_datafile) >> alpha_in;
+ *(ad_comm::global_datafile) >> sigmar_in;
+ *(ad_comm::global_datafile) >> rho_in;
+ // rho_in=0.82;
+ }
+ // Open up tac-model parameters
+ ad_comm::change_datafile_name("tacpar.dat");
+ END_CALCS
+ init_int nntmp
+ init_int nnodes
+ init_vector maxabc(1,ntacspp)
+ init_matrix theta(0,nnodes,1,ntacspp)
+ !! cout<<"read tacpar"<> spname(i); // 1
+ cout<> SSL_spp(i); // 2
+ *(ad_comm::global_datafile) >> Const_Buffer(i); // 3
+ *(ad_comm::global_datafile) >> ngear(i); // 4
+ *(ad_comm::global_datafile) >> nsexes(i); // 6
+ *(ad_comm::global_datafile) >> avg_5yrF(i); // 7
+ *(ad_comm::global_datafile) >> FABC_Adj(i); // 8
+ *(ad_comm::global_datafile) >> SPR_abc(i); // 9
+ *(ad_comm::global_datafile) >> SPR_ofl(i); // 10
+ *(ad_comm::global_datafile) >> spawnmo(i); // 11
+ *(ad_comm::global_datafile) >> nages(i); // 12
+ cout<<"nages: "<> Fratiotmp(i,j); // 13
+ write_log( Fratiotmp(i));
+
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> M_Ftmp(i,k); // 14
+ if (nsexes(i)==2)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> M_Mtmp(i,k); // 15
+ write_log( M_Ftmp(i));
+
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> pmaturetmp_F(i,k); // 16
+ write_log( pmaturetmp_F(i));
+
+ if (nsexes(i)==2)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> pmaturetmp_M(i,k); // 15
+ write_log( pmaturetmp_M(i));
+ cout << "Mature: "<< pmaturetmp_F(i)(1,nages(i)) <> wt_Ftmp(i,k); // 17
+ write_log( wt_Ftmp(i));
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> wt_gear_Ftmp(i,j,k); // 18
+ write_log( wt_gear_Ftmp(i));
+
+ if (nsexes(i)==2)
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> wt_gear_Mtmp(i,j,k);// 19
+ write_log( wt_gear_Mtmp(i));
+
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> sel_Ftmp(i,j,k); // 20
+ write_log( sel_Ftmp(i));
+
+ if (nsexes(i)==2)
+ for (int j=1;j<=ngear(i);j++)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> sel_Mtmp(i,j,k); // 21
+ write_log( sel_Mtmp(i));
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> n0_Ftmp(i,k); // 22
+ write_log( n0_Ftmp(i));
+
+ if (nsexes(i)==2)
+ for (int k=1;k<=nages(i);k++)
+ *(ad_comm::global_datafile) >> n0_Mtmp(i,k); // 23
+ write_log( n0_Mtmp(i));
+
+ cout<<"N: "<> nrec(i); // 24
+ cout<<"nrec: "<> Rtmp(i,j); // 25
+ cout<<"Rec: "<> SSBtmp(i,j); // 26
+ cout<<"SSB: "<> Rsim(ispp,i,j);
+ Rsim(ispp,i,j) *= exp(rnorms(i,j)*.375); // about 15% CV to get historical mean
+ }
+ if (nsims<=5) Rsim(ispp,i,j) = AMeanRec(ispp); // XXX constant recruitment
+ }
+ }
+ envin.close();
+ AMeanRec(ispp) *= .5; // Goes to females only
+ cout <<"Mean recruits "<< AMeanRec(ispp) < 0 )
+ F_begin_yr(k,ispp) = SolveF2(n0_F(ispp), n0_M(ispp), Obs_Catch(k,ispp),ispp); // F_yr_one(ispp) = SolveF2(n0_F(ispp), n0_M(ispp), yr_one_catch(ispp),ispp);
+ else
+ F_begin_yr(ispp) = 0; // F_yr_one(ispp) = 0; // cout< 1) // Debugging if LP to be done or not
+ {
+ int on=0;
+ if ( (on=option_match(ad_comm::argc,ad_comm::argv,"-nolp"))>-1)
+ dolp=0;
+ else
+ dolp=1;
+ }
+ END_CALCS
+
+
+PRELIMINARY_CALCS_SECTION
+ double tmp1;
+ write_alts_hdr();
+ get_SB100();
+ cout<<"SSB, Biomass at unfished: "<1)
+ nyrs_catch = nyrs_catch_in;
+ //else
+ // nyrs_catch = 1;// NOTA BUENO: this is changed under the new EIS Alternatives (May 06)
+ Do_Sims();
+ // if (alt==2) write_alts();
+ write_alts();
+ }
+
+FUNCTION compute_obj_fun
+ dvariable tmp1;
+ dvariable tmp2;
+ tmp1.initialize();
+ for (int ispp=1;ispp<=nspp;ispp++)
+ {
+ Get_Bzero(ispp);
+ // write_srec();exit(1);
+ if(Fmsy_F35>0)
+ {
+ get_msy(ispp);
+ switch (current_phase())
+ {
+ case 1 : tmp1 = 1.e1*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ case 2 : tmp1 = 1.e2*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ case 3 : tmp1 = 1.e2*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ case 4 : tmp1 = 1.e3*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ default: tmp1 = 1.e3*square(log(Fmsy(ispp))-log(F35(ispp))); break;
+ }
+ if(Fmsy_F35==2)
+ switch (current_phase())
+ {
+ case 1 : tmp1 += 1.e1*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ case 2 : tmp1 += 1.e2*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ case 3 : tmp1 += 1.e2*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ case 4 : tmp1 += 1.e3*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ default : tmp1+= 1.e3*square(log(Bmsy(ispp))-log(0.35 * SB100(ispp))); break;
+ }
+ obj_fun += tmp1;
+ }
+ if (Rec_Cond>0.)
+ {
+ // More conditioning here--make mean recruitment consistent with half and double
+ // recruitment at avg spawning biomass...
+ double vartmp = 2.*Rec_Cond*Rec_Cond;
+ /*
+ double ssb1 = AMeanSSB(ispp) * 0.5;
+ double ssb2 = AMeanSSB(ispp) ;
+ double ssb3 = AMeanSSB(ispp) * 2.0;
+ // double ssb4 = value(Bzero(ispp)) ;
+ tmp2.initialize();
+ tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb1, ispp)))/ vartmp;
+ tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb2, ispp)))/ vartmp;
+ tmp2 += square(log(AMeanRec(ispp)) - log(SRecruit( ssb3, ispp)))/ vartmp;
+ // tmp2 += 24.* square(log(AMeanRec(ispp)) - log(SRecruit( ssb4, ispp)));
+ // tmp2 += 24.* square(log(AMeanSSB(ispp)) - log( ssb4 ));
+ */
+ tmp2 = square(log(Bzero(ispp)) - log(SB100(ispp)))/ vartmp;
+ obj_fun += tmp2;
+ // cout << "SSB "<OY_max)
+ Actual_Catch = TAC/sum(TAC)*OY_max;
+ else
+ Actual_Catch = TAC;
+ }
+ /* ////////////////////////////////////////////////////////////////////////////// //int use_max = 1; // Need to move this into the setup file... double max_catch=1500.;// EBS POllock special case (CHANGE THIS) //if (use_max==1) for (int ispp=1;ispp<=nspp;ispp++) Actual_Catch(ispp) = min(TAC(ispp),max_catch); // cout<<"AC: "<3)
+ {
+ Alt4_Fcalc = 1;
+ double diff;
+ double ssl_sum=0.;
+ diff = OY_min-sumabc;
+ // if so, sum up SSL ABCs
+ ssl_sum = SSL_spp*ABC;
+ for (int ispp=1;ispp<=nspp;ispp++)
+ {
+ // then apportion them to difference between current ABC and target (OY_min) difference
+ if (SSL_spp(ispp))
+ {
+ TAC(ispp) = ABC(ispp) + diff * ABC(ispp)/ssl_sum;
+ // Special for BSAI Pcod to subtract off 3%...and get the totals to match up
+ // if (ispp==2) TAC(ispp) /= 0.97; For EIS work
+ }
+ else
+ {
+ TAC(ispp) = ABC(ispp);
+ }
+ }
+ // cout << ipro<<" "<0.01) ftmp = SolveF2(N_F(ispp),N_M(ispp),TAC(ispp),ispp);
+ // cout <= 0.2 *SBzero(ispp) & SBtmp < SBKink(ispp) ) // NOTE Same as Am 56 here (until 20% bzero reached)
+ Ftmp = (Fabc(ispp)*(1/(1-alpha ))*(SBtmp/SBKink(ispp) - alpha ));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fabc(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ return(Ftmp);
+
+FUNCTION double Get_F_Am56(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ {
+ for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct
+ {
+ if (SBtmp < alpha*SBKink(ispp))
+ Ftmp = 0.;
+ if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) )
+ Ftmp = (Fabc(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fabc(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ }
+ return(Ftmp);
+
+FUNCTION double Get_Fofl_t2(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ {
+ for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct
+ {
+ if (SBtmp < alpha*SBKink(ispp))
+ Ftmp = 0.;
+ if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) )
+ Ftmp = (Fofl(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fofl(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ }
+ return(Ftmp);
+
+FUNCTION double Get_Fofl_t(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ {
+ for (ii=1;ii<=3;ii++) // Iterate to get month of spawning correct
+ {
+ if (SBtmp < alpha*SBKink(ispp))
+ Ftmp = 0.;
+ if (SBtmp >= alpha*SBKink(ispp) & SBtmp < SBKink(ispp) )
+ Ftmp = (Fofl(ispp)*(1/(1-alpha))*(SBtmp/SBKink(ispp) - alpha));
+ if (SBtmp > SBKink(ispp) )
+ Ftmp = (Fofl(ispp));
+ SBtmp = N_females * elem_prod( wt_mature_F(ispp), mfexp( -yrfrac(ispp)*(M_F(ispp) + Ftmp*F_age )));
+ }
+ }
+ return(Ftmp);
+
+
+FUNCTION double Get_F_t(const dvector& F_age, const dvector& N_females, const int ispp )
+ double Ftmp; // cout<< spname(ispp)<< endl;
+ if (SSL_spp(ispp))
+ {
+ Ftmp = Get_F_SSL_prey(F_age, N_females, ispp);
+ }
+ else
+ {
+ Ftmp = Get_F_Am56(F_age, N_females, ispp);
+ }
+ return(Ftmp);
+
+FUNCTION void Project_Pops(const int& isim, const int& i)
+ double ctmp;
+ if ( i == npro && isim%int(nsims/4)==0) cout << "Year "< 0)
+ {
+ if (i <= nyrs_catch || TAC_ABC==0) // Use TAC setting algorithm for alt 2 only, for all others, set TAC==ABC
+ {
+ Ftmp = SolveF2(N_F(ispp),N_M(ispp),Actual_Catch(ispp),ispp);
+ }
+ else
+ {
+ if (alt==7 && i<=3)
+ Ftmp = Get_F(1,ispp); // Set to the F rather than solving every time...
+ else
+ if (alt==77 && i<=2)
+ Ftmp = Get_F(1,ispp); // Set to the F rather than solving every time...
+ else
+ {
+ // if (alt==2 && sum(TAC)>OY_max)
+ if ((alt==2 ||alt==98||alt==97) && sum(TAC)>OY_max)
+ {
+ Ftmp = SolveF2(N_F(ispp),N_M(ispp),OY_max*TAC(ispp)/sum(TAC),ispp);
+ // cout<<"HERE "<< Ftmp<<" "< 0.)
+ {
+ for (m=1;m<=ngear(ispp);m++)
+ {
+ Ftottmp_F = Ftmp*Frat(ispp,m)*sel_F(ispp,m);
+ Ftottmp_spr += Ftottmp_F ;
+ Ftottmp_M = Ftmp*Frat(ispp,m)*sel_M(ispp,m);
+ ctmp += (wt_gear_F(ispp,m) * elem_prod(elem_div(Ftottmp_F, Z_F(ispp)),elem_prod(1.-S_F(ispp),N_F(ispp)))); // Catch equation (vectors)
+ ctmp += (wt_gear_M(ispp,m) * elem_prod(elem_div(Ftottmp_M, Z_M(ispp)),elem_prod(1.-S_M(ispp),N_M(ispp)))); // Catch equation (vectors)
+ }
+ SPRsim(ispp,isim,i) = Implied_SPR(Ftottmp_spr,ispp);
+ // cout <NUL ");
+ ifstream tac("tac.dat");
+ tac>>TAC;
+ tac.close(); /* */
+ for (int ispp=1;ispp<=nspp;ispp++)
+ Actual_Catch(ispp) = min(TAC(ispp),ABC(ispp));
+ }
+
+
+FUNCTION Avg_Age
+ Avg_Age_End.initialize();
+ Avg_Age_Mat.initialize();
+ for (int ispp=1;ispp<=nspp;ispp++)
+ {
+ // if(!isit_const(ispp) )
+ {
+ dvector age_seq(1,nages(ispp));
+ age_seq.initialize();
+ age_seq.fill_seqadd(1,1);
+ dvector ntmp(1,nages(ispp));
+ ntmp = N_F(ispp);
+ Avg_Age_End(ispp) += age_seq * ntmp /sum(ntmp);
+ Avg_Age_sum(ispp) += Avg_Age_End(ispp) ;
+ ntmp = elem_prod(ntmp,pmature_F(ispp)) ;
+ Avg_Age_Mat(ispp) += elem_prod(ntmp,pmature_F(ispp)) * age_seq/sum(ntmp);
+ // cout< 1e-6)
+ {
+ iter++;
+ ftmp += (TACin-cc) / btmp;
+ Ftottmp_F.initialize();
+ Ftottmp_M.initialize();
+ for (m=1;m<=ngear(ispp);m++)
+ {
+ Ftottmp_F += ftmp*Fratsel_F(m);
+ Ftottmp_M += ftmp*Fratsel_M(m);
+ }
+ Z_F = Ftottmp_F + M_F(ispp);
+ Z_M = Ftottmp_M + M_M(ispp);
+ S_F = mfexp( -Z_F );
+ S_M = mfexp( -Z_M );
+ cc = 0.0;
+ for (m=1;m<=ngear(ispp);m++)
+ {
+ cc += (wt_gear_F(ispp,m) * elem_prod(elem_div(ftmp*Fratsel_F(m), Z_F),elem_prod(1.-S_F,N_F))); // Catch equation (vectors)
+ cc += (wt_gear_M(ispp,m) * elem_prod(elem_div(ftmp*Fratsel_M(m), Z_M),elem_prod(1.-S_M,N_M))); // Catch equation (vectors)
+ }
+ dd = cc / TACin - 1.;
+ //cout << ispp<<" "<< ftmp << " "<< cc << " "<100) {cerr<<"Bombed on catch solver--check scales for "<nages(ispp))
+ {
+ Tg += double(ii) * wt_mature_F(ispp,nages(ispp)) * ntmp;
+ tmp += wt_mature_F(ispp,nages(ispp)) * ntmp;
+ ntmp *= exp(-M_F(ispp,nages(ispp)));
+ }
+ else
+ {
+ Tg += double(ii) * wt_mature_F(ispp,ii) * ntmp;
+ tmp += wt_mature_F(ispp,ii) * ntmp;
+ ntmp *= exp(-M_F(ispp,ii));
+ }
+ }
+ // Tg /= value(phizero(ispp));
+ Tg /= tmp;
+ report << Tg<< " ";
+ }
+ report << endl;
+ report << "Options SR_Type Project_Recr_Assmption SR_Condition"<2)
+ write_sim("Alternative ",ispp);
+ }
+ write_by_time();
+
+ // Writing routines here ....
+FUNCTION write_by_time
+ write_sim("Catch", Csim,Cabc); // Total Catch
+ cout< 1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+FUNCTION void write_TACs(const adstring& Title)
+ // This one prints out species over time (means only), but without headings
+ means_out <<"Alternative "< 1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+
+
+FUNCTION void write_sim(const adstring& Title, d3_array& Outtmp)
+ // This one prints out species over time (means only), but without headings
+ d3_array mtmp(1,nspp,1,npro,1,nsims);
+ for (int ispp=1;ispp<=nspp;ispp++) for (int i=1;i<=npro;i++) for (int k=1;k<=nsims;k++)
+ mtmp(ispp,i,k)=Outtmp(ispp,k,i);
+ means_out <<"Alternative "< 1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+FUNCTION void write_sim(const adstring& Title,const d3_array& Outtmp,const dvar_vector& bb)
+ // This one prints out species over time (means only)
+ d3_array mtmp(1,nspp,1,npro,1,nsims);
+ for (int ispp=1;ispp<=nspp;ispp++) for (int i=1;i<=npro;i++) for (int k=1;k<=nsims;k++)
+ mtmp(ispp,i,k)=Outtmp(ispp,k,i);
+ means_out <<"Alternative "<1e-6)
+ means_out << mean_value <<" ";
+ else
+ means_out << " NA ";
+ }
+ means_out << endl;
+ }
+ means_out << endl;
+
+FUNCTION void write_sim(const adstring& Title,const int& ispp)
+ // For each species separately
+ percent_out <<"Alternative "<
+ adstring xspname;
+ adstring_array targsppname(1,20);
+ adstring_array spp_file_name(1,20);
+ adstring_array gearname(1,8);
+ adstring_array spname(1,90);
+ adstring_array areaname(1,8);
+ adstring run_name;
+
+ ofstream write_log("Input_Log.rep");
+ #undef write_log
+ #define write_log(object) write_log << #object "\n" << object << endl;
+
+ // A routine to get transpose, sort and return a matrix ----
+ dmatrix TranSort (const dmatrix m1)
+ {
+ RETURN_ARRAYS_INCREMENT();
+ dmatrix vtmp=trans(m1);
+ int npro=m1.colmax();
+ int nsim=m1.rowmax();
+ for (int i=1;i<=npro;i++)
+ vtmp(i) = sort(vtmp(i),nsim);
+
+ RETURN_ARRAYS_DECREMENT();
+ return(vtmp);
+ }
+ // #include
+
+FUNCTION dvariable SRecruit(const dvariable& Stmp,const int& ispp)
+ RETURN_ARRAYS_INCREMENT();
+ dvariable RecTmp;
+ switch (SrType)
+ {
+ case 1:
+ RecTmp = (Stmp / phizero(ispp)) * mfexp( sr_alpha(ispp) * ( 1. - Stmp / Bzero(ispp) )) ; //Ricker form from Dorn
+ break;
+ case 2:
+ RecTmp = Stmp / ( sr_alpha(ispp) + beta(ispp) * Stmp); //Beverton-Holt form
+ break;
+ case 3:
+ // RecTmp = mfexp(log_avgrec); //Avg recruitment
+ break;
+ case 4:
+ RecTmp = Stmp * mfexp( sr_alpha(ispp) - Stmp * beta(ispp)) ; //old Ricker form
+ break;
+ }
+
+ RETURN_ARRAYS_DECREMENT();
+ return RecTmp;
+
+
+FUNCTION dvar_vector SRecruit(const dvar_vector& Stmp,const int& ispp)
+ RETURN_ARRAYS_INCREMENT();
+ dvar_vector RecTmp(Stmp.indexmin(),Stmp.indexmax());
+ switch (SrType)
+ {
+ case 1:
+ RecTmp = elem_prod((Stmp / phizero(ispp)) , mfexp( sr_alpha(ispp) * ( 1. - Stmp / Bzero(ispp) ))) ; //Ricker form from Dorn
+ break;
+ case 2:
+ RecTmp = elem_prod(Stmp , 1. / ( sr_alpha(ispp) + beta(ispp) * Stmp)); //Beverton-Holt form
+ break;
+ case 3:
+ // RecTmp = mfexp(log_avgrec); //Avg recruitment
+ break;
+ case 4:
+ RecTmp = elem_prod(Stmp ,mfexp( sr_alpha(ispp) - Stmp * beta(ispp))); //old Ricker form
+ break;
+ }
+ RETURN_ARRAYS_DECREMENT();
+ //cout <5||F1<0.01))
+ // {
+ // ii=5;
+ // F1=M(ispp,nages(ispp)); // When things bomb (F <0 or F really big then just set it to F35...)
+ // }
+ // else
+ {
+ F2 = F1 + df*.5;
+ F3 = F2 - df; //F1 = double(ii)/100;
+ yld1 = yield(F1,Stmp,Rtmp,ispp); //cout <