#' Wrapper function to run an MSE.
#'
#' Function that runs the entire MSE- generates data from OM, assesses it with EM, and forecasts. Can be run on its own, but it is intended to be run with \code{\link{run_GeMS}}
#'
#' @param out Named list output from \code{\link{ReadCTLfile}}
#' @param CTLName Name of CTL file
#' @param MSEdir Directory containing CTL files. Full absolute path preferred.
#' @param silent Logical; Showing output on console
#' @param ADoptions String of ADMB options to be passed. Always includes \code{-nohess}
#' @param ADsilent Logical; Showing ADMB output on console
#' @param echo Logical; Write to an echo file to show quantities from OM
#' @param ... WHO KNOWS
#'
#' @return GeMS Object # To be implemented
#'
#' @export
#' @examples
#' \dontrun{
#' CTLName <- "Cod_LowProd_CTL"
#' out <- ReadCTLfile(CTLName)
#' MSEdir <- "~/Examples/Cod_1_Production"
#' GeMS(out,CTLName,MSEdir,ADoptions="-gbs 1000000")
#' }
GeMS<-function(out,CTLName,MSEdir=getwd(),silent=F,ADoptions=NA,ADsilent=T,echo=F,...)
{
CurDir<-getwd()
if(!file.exists(file.path(MSEdir,paste0(CTLName,".csv")))) stop(paste0(file.path(MSEdir,paste0(CTLName,".csv")), " does not exist. Set MSEdir to directory containing CTL files."))
if(echo) echofile <- file.path(MSEdir,paste0(CTLName,"_echo.txt"))
if(echo) {if(!file.exists(echofile)) file.create(echofile)}
if(echo) cat(paste0("# Run initiated: ",Sys.time(),"\n"),file=echofile,append=F)
if(echo) cat(paste0("# Best viewed in Excel or some tab-delimiting editor.","\n"),file=echofile,append=T)
if(echo) cat("## Inputs and Derived Quantities \n",file=echofile,append=T)
if(echo) cat(CTLName,"\t # CTLName \n",file=echofile,append=T)
#===============
#==Checking OS==
#===============
os <- .Platform$OS.type
if(os != "windows") {
SimAssExec <- "SimAss"
SimAssComm <- "./SimAss -nohess"
}
if(os == "windows") {
SimAssExec <- "simass.exe"
SimAssComm <- "simass -nohess"
}
if(ADsilent) {
SimAssComm <- paste(SimAssComm, ">console.log", sep = " ") # saves console output to a file "console.log"
}
if(is.na(ADoptions[1])) {
SimAssComm <- paste(SimAssComm, ADoptions, sep = " ") # ADMB options such as "-gbs" or "-cbs" (Memory Management)
}
#=======================
#==simulation controls==
#=======================
Nsim <-out$OM$Nsim # number of simulations to do in the MSE
if(echo) cat(Nsim,"\t # Nsim \n",file=echofile,append=T)
SimYear <-out$OM$SimYear # total number of years in simulation
if(echo) cat(SimYear,"\t # SimYear \n",file=echofile,append=T)
InitYear <-out$OM$InitYear # year in which MSE starts (i.e. the number of years of data available for initial assessment)
if(echo) cat(InitYear,"\t # InitYear \n",file=echofile,append=T)
AssessmentType <-out$OM$AssessmentType # 0 = projection only, 1 = production, 2= agestructured
if(echo) cat(AssessmentType,"\t # AssessmentType \n",file=echofile,append=T)
FisheryIndepenDat <-out$OM$FisheryIndepenDat
if(echo) cat(FisheryIndepenDat,"\t # FisheryIndepenDat \n",file=echofile,append=T)
LifeHistoryPlots <-out$OM$LifeHistoryPlots
EstimationPlots <-out$OM$EstimationPlots # this plots the diagnostics for each run of the estimation model (will be a lot of plots!)
AssessmentData <-out$OM$AssessmentData # this plots the diagnostics for each run of the estimation model (will be a lot of plots!)
PlotYieldCurve <-out$OM$PlotYieldCurve
TwoPop <-out$OM$TwoPop
PlotFolder<-file.path(MSEdir, "plots")
if(!file.exists(PlotFolder)) dir.create(PlotFolder)
#==sampling uncertainty
CatchCVn <-CleanInput(out$OM$CatchCVn,SimYear)
if(echo) cat(CatchCVn,"\t # CatchCVn \n",file=echofile,append=T)
CatchCVs <-CatchCVn
if(TwoPop>0)
CatchCVs <-CleanInput(out$OM$CatchCVs,SimYear)
IndexCVn <-CleanInput(out$OM$IndexCVn,SimYear)
if(echo) cat(IndexCVn,"\t # IndexCVn \n",file=echofile,append=T)
IndexCVs <-IndexCVn
if(TwoPop>0)
IndexCVs <-CleanInput(out$OM$IndexCVs,SimYear)
LenSampleN <-CleanInput(out$OM$LenSampleN,SimYear)
if(echo) cat(LenSampleN,"\t # LenSampleN \n",file=echofile,append=T)
LenSampleS <-LenSampleN
if(TwoPop>0)
LenSampleS <-CleanInput(out$OM$LenSampleS,SimYear)
GrowthSDn <-CleanInput(out$OM$GrowthSDn,SimYear)
if(echo) cat(GrowthSDn,"\t # GrowthSDn \n",file=echofile,append=T)
GrowthSDs <-GrowthSDn
if(TwoPop>0)
GrowthSDs <-CleanInput(out$OM$GrowthSDs,SimYear)
#==========================================================
#=================population dynamics processes============
#==============================================+===========
MaxAge<-out$OM$MaxAge
if(echo) cat(MaxAge,"\t # MaxAge \n",file=echofile,append=T)
Ages<-seq(1,MaxAge)
#==natural mortality================
NatMn <-CleanInput(out$OM$NatMn,SimYear)
if(echo) cat(NatMn,"\t # NatMn \n",file=echofile,append=T)
NatMs <-NatMn
if(TwoPop>0)
NatMs <-CleanInput(out$OM$NatMn,SimYear)
#==Length at age====================
VonKn <-CleanInput(out$OM$VonKn,SimYear)
if(echo) cat(VonKn,"\t # VonKn \n",file=echofile,append=T)
LinfN <-CleanInput(out$OM$LinfN,SimYear)
if(echo) cat(LinfN,"\t # LinfN \n",file=echofile,append=T)
t0n <-CleanInput(out$OM$t0n,SimYear)
if(echo) cat(t0n,"\t # t0n \n",file=echofile,append=T)
LenAtAgeN <-matrix(nrow=SimYear,ncol=MaxAge)
LenAtAgeN[1,]<-LinfN[1]*(1-exp(-VonKn[1]*(Ages-t0n[1])))
for(i in 2:SimYear) {
LenAtAgeN[i,1]<-LinfN[i]*(1-exp(-VonKn[i]*(Ages[1]-t0n[i])))
for(j in 2:MaxAge)
LenAtAgeN[i,j]<-LenAtAgeN[i-1,j-1]+(LinfN[i]-LenAtAgeN[i-1,j-1])*(1-exp(-VonKn[i]))
}
if(echo) cat("# LenAtAgeN \n",file=echofile,append=T)
if(echo) cat("# Ages", Ages,"\n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,LenAtAgeN,file=echofile,append=T)
LenAtAgeS<-LenAtAgeN
if(TwoPop>0)
{
VonKs <-CleanInput(out$OM$VonKs,SimYear)
LinfS <-CleanInput(out$OM$LinfS,SimYear)
t0s <-CleanInput(out$OM$t0s,SimYear)
LenAtAgeS <-matrix(nrow=SimYear,ncol=MaxAge)
LenAtAgeS[1,]<-LinfS[1]*(1-exp(-VonKs[1]*(Ages-t0s[1])))
LenAtAgeS[,1]<-LenAtAgeS[1,1]
for(i in 2:SimYear) {
LenAtAgeS[i,1]<-LinfS[i]*(1-exp(-VonKs[i]*(Ages[1]-t0s[i])))
for(j in 2:MaxAge)
LenAtAgeS[i,j]<-LenAtAgeS[i-1,j-1]+(LinfS[i]-LenAtAgeS[i-1,j-1])*(1-exp(-VonKs[i]))
}
}
#==specify the number of length bins
BinWidth <-(max(LenAtAgeS,LenAtAgeN)*1.05)/(out$OM$LengthBinN)
if(echo) cat("\n",BinWidth,"\t # BinWidth \n",file=echofile,append=T)
LengthBins <-seq(0,max(LenAtAgeS,LenAtAgeN)*1.05,BinWidth)
if(echo) cat(LengthBins,"\t # LengthBins \n",file=echofile,append=T)
LengthBinsMid <-LengthBins[1:(length(LengthBins)-1)] + mean(LengthBins[1:2])
LengthBinN <-length(LengthBinsMid)
if(echo) cat(LengthBinN,"\t # LengthBinN \n",file=echofile,append=T)
#==maturity at age==========================
mat50n <-CleanInput(out$OM$mat50n,SimYear)
if(echo) cat(mat50n,"\t # mat50n \n",file=echofile,append=T)
mat95n <-CleanInput(out$OM$mat95n,SimYear)
if(echo) cat(mat95n,"\t # mat95n \n",file=echofile,append=T)
matureN <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
matureN[i,]<-1/(1+exp(-1*log(19)*(Ages-mat50n[i])/(mat95n[i]-mat50n[i])))
if(echo) cat("# matureN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,matureN,file=echofile,append=T)
matureS <-matureN
if(TwoPop>0)
{
mat50s <-CleanInput(out$OM$mat50s,SimYear)
mat95s <-CleanInput(out$OM$mat95s,SimYear)
matureS <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
matureS[i,]<-1/(1+exp(-1*log(19)*(Ages-mat50s[i])/(mat95s[i]-mat50s[i])))
}
#=====================================
#==fishery selectivitiy==============
sel50n <-CleanInput(out$OM$sel50n,SimYear)
if(echo) cat(sel50n,"\t # sel50n \n",file=echofile,append=T)
sel95n <-CleanInput(out$OM$sel95n,SimYear)
if(echo) cat(sel95n,"\t # sel95n \n",file=echofile,append=T)
vulnN <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
vulnN[i,] <-1/(1+exp(-1*log(19)*(LenAtAgeN[i,]-sel50n[i])/(sel95n[i]-sel50n[i])))
vulnN[vulnN<0.01]<-0
if(echo) cat("\n # vulnN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,vulnN,file=echofile,append=T)
vulnS <-vulnN
if(TwoPop>0)
{
sel50s <-CleanInput(out$OM$sel50s,SimYear)
sel95s <-CleanInput(out$OM$sel95s,SimYear)
vulnS <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
vulnS[i,] <-1/(1+exp(-1*log(19)*(LenAtAgeN[i,]-sel50s[i])/(sel95s[i]-sel50s[i])))
vulnS[vulnS<0.01]<-0
}
#==index selectivity=====================
surv50n <-CleanInput(out$OM$surv50n,SimYear)
if(echo) cat(surv50n,"\t # surv50n \n",file=echofile,append=T)
surv95n <-CleanInput(out$OM$surv95n,SimYear)
if(echo) cat(surv95n,"\t # surv95n \n",file=echofile,append=T)
survSelN <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
survSelN[i,] <-1/(1+exp(-1*log(19)*(LenAtAgeN[i,]-surv50n[i])/(surv95n[i]-surv50n[i])))
if(echo) cat("# survSelN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,survSelN,file=echofile,append=T)
survSelS <-survSelN
surv50s <-surv50n
surv95s <-surv95n
if(TwoPop>0)
{
surv50s <-CleanInput(out$OM$surv50s,SimYear)
surv95s <-CleanInput(out$OM$surv95s,SimYear)
survSelS <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
survSelS[i,] <-1/(1+exp(-1*log(19)*(LenAtAgeN[i,]-surv50s[i])/(surv95s[i]-surv50s[i])))
}
if(max(sel50n)>max(LinfN) | max(sel95n)>max(LinfN) | max(surv50n)>max(LinfN) | max(surv50n)>max(LinfN)) {stop("Gear is attempting to select fish larger than what's in the population (one or more of the selectivity parameters is greater than Linf)")}
if(TwoPop>0) if(max(sel50s)>max(LinfS) | max(sel95s)>max(LinfS) | max(surv50s)>max(LinfS) | max(surv50s)>max(LinfS)) {stop("Gear is attempting to select fish larger than what's in the population (one or more of the selectivity parameters is greater than Linf)")}
#==weight at age==========================
alphaN <-CleanInput(out$OM$alphaN,SimYear)
if(echo) cat(alphaN,"\t # alphaN \n",file=echofile,append=T)
betaN <-CleanInput(out$OM$betaN,SimYear)
if(echo) cat(betaN,"\t # betaN \n",file=echofile,append=T)
WeightAtAgeN <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
WeightAtAgeN[i,] <-alphaN[i]*LenAtAgeN[i,]^betaN[i]
if(echo) cat("\n # WeightAtAgeN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,WeightAtAgeN,file=echofile,append=T)
WeightAtAgeS <-WeightAtAgeN
if(TwoPop>0)
{
alphaS <-CleanInput(out$OM$alphaS,SimYear)
betaS <-CleanInput(out$OM$betaS,SimYear)
WeightAtAgeS <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
WeightAtAgeS[i,] <-alphaS[i]*LenAtAgeS[i,]^betaS[i]
}
#==movement from box to box==
MaxMovingN <-CleanInput(out$OM$MaxMovingN,SimYear)
if(echo) cat(MaxMovingN,"\t # MaxMovingN \n",file=echofile,append=T)
Move50n <-CleanInput(out$OM$Move50n,SimYear)
if(echo) cat(Move50n,"\t # Move50n \n",file=echofile,append=T)
Move95n <-CleanInput(out$OM$Move95n,SimYear)
if(echo) cat(Move95n,"\t # Move95n \n",file=echofile,append=T)
MovementN <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
MovementN[i,] <-MaxMovingN[i]/(1+exp(-1*log(19)*(Ages-Move50n[i])/(Move95n[i]-Move50n[i])))
if(echo) cat("\n # MovementN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,MovementN,file=echofile,append=T)
MovementS <-MovementN
if(TwoPop>0)
{
MaxMovingS <-CleanInput(out$OM$MaxMovingS,SimYear)
Move50s <-CleanInput(out$OM$Move50s,SimYear)
Move95s <-CleanInput(out$OM$Move95s,SimYear)
MovementS <-matrix(nrow=SimYear,ncol=MaxAge)
for(i in 1:SimYear)
MovementS[i,] <-MaxMovingS[i]/(1+exp(-1*log(19)*(Ages-Move50s[i])/(Move95s[i]-Move50s[i])))
}
if(!silent) {
if(sum(MovementS)>0 | sum(MovementN)>0)
{print("Movement is occuring between populations")}
}
#==recruitment parameters==
steepnessN <-CleanInput(out$OM$steepnessN,SimYear)
if(echo) cat(steepnessN,"\t # steepnessN \n",file=echofile,append=T)
sigmaRn <-CleanInput(out$OM$sigmaRn,SimYear)
if(echo) cat(sigmaRn,"\t # sigmaRn \n",file=echofile,append=T)
RzeroN <-CleanInput(out$OM$RzeroN,SimYear)
if(echo) cat(RzeroN,"\t # RzeroN \n",file=echofile,append=T)
steepnessS <-steepnessN
sigmaRs <-sigmaRn
RzeroS <-RzeroN
if(TwoPop>0)
{
steepnessS <-CleanInput(out$OM$steepnessS,SimYear)
sigmaRs <-CleanInput(out$OM$sigmaRs,SimYear)
RzeroS <-CleanInput(out$OM$RzeroS,SimYear)
}
RecErrN <-matrix(rnorm(SimYear*Nsim,0,sigmaRn),ncol=SimYear,byrow=T)
if(echo) cat("# RecErrN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,RecErrN,file=echofile,append=T)
RecErrS <-matrix(rnorm(SimYear*Nsim,0,sigmaRs),ncol=SimYear,byrow=T)
#==historic fishing mortality
HistoricalFn <-CleanInput(out$OM$HistoricalFn,SimYear)[1:InitYear]
if(echo) cat(HistoricalFn,"\t # HistoricalFn \n",file=echofile,append=T)
PastFsdN <-out$OM$PastFsdN
HarvestControlN <-out$OM$HarvestControlN
ConstantCatchN <-out$OM$ConstantCatchN
ConstantFn <-out$OM$ConstantFn
HCalphaN <-out$OM$HCalphaN
HCbetaN <-out$OM$HCbetaN
HistoricalFs <-HistoricalFn
PastFsdS <-PastFsdN
HarvestControlS <-HarvestControlN
ConstantCatchS <-ConstantCatchN
ConstantFs <-ConstantFn
HCalphaS <-HCalphaN
HCbetaS <-HCbetaN
if(TwoPop>0)
{
HistoricalFs <-CleanInput(out$OM$HistoricalFs,SimYear)[1:InitYear]
PastFsdS <-out$OM$PastFsdS
HarvestControlS <-out$OM$HarvestControlS
ConstantCatchS <-out$OM$ConstantCatchS
ConstantFs <-out$OM$ConstantFs
HCalphaS <-out$OM$HCalphaS
HCbetaS <-out$OM$HCbetaS
}
start_assessment<-out$OM$start_assessment
if(echo) cat(start_assessment,"\t # start_assessment \n",file=echofile,append=T)
SmallNum <-out$OM$SmalNum
InitSmooth <-out$OM$InitSmooth
FmortPen <-out$OM$FmortPen
RecruitPen <-out$OM$RecruitPen
Mpenalty <-out$OM$Mpenalty
Growthpenalty <-out$OM$Growthpenalty
SelPenalty <-out$OM$SelPenalty
EstM <-out$OM$EstM
TimeVaryM <-out$OM$TimeVaryM
EstGrowthK <-out$OM$EstGrowthK
TimeVaryGrowthK <-out$OM$TimeVaryGrowthK
EstLinf <-out$OM$EstLinf
TimeVaryLinf <-out$OM$TimeVaryLinf
TimeVarySel50 <-out$OM$TimeVarySel50
TimeVarySel95 <-out$OM$TimeVarySel95
ProjectTimeVary <-out$OM$ProjectTimeVary
InitValSig <-out$OM$InitValSig
InitBzeroMod <-out$OM$InitBzeroMod
InitGrowthRate <-out$OM$InitGrowthRate
InitShape <-out$OM$InitShape
estInit <-out$OM$estInit
InitBioProd <-out$OM$InitBioProd
#===================================================
#==Virgin numbers at age, biomass, recruitment
#===================================================
if(echo) cat("\n\n # Initialised Population \n",file=echofile,append=T)
VirInitN <-initialN(Rzero=RzeroN[1],NatM=NatMn[1],inAge=MaxAge)
if(echo) cat("# VirInitN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,VirInitN,file=echofile,append=T)
VirInitS <-initialN(Rzero=RzeroS[1],NatM=NatMs[1],inAge=MaxAge)
VirBioN <-sum(VirInitN*matureN[1,]*WeightAtAgeN[1,])
if(echo) cat("\n # VirBioN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,VirBioN,file=echofile,append=T)
VirBioS <-sum(VirInitS*matureS[1,]*WeightAtAgeS[1,])
ExploitBioN <-sum(VirInitN*vulnN[1,]*WeightAtAgeN[1,])
if(echo) cat("\n # ExploitBioN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,ExploitBioN,file=echofile,append=T)
ExploitBioS <-sum(VirInitS*vulnS[1,]*WeightAtAgeS[1,])
png(file.path(PlotFolder,paste0("LifeHistory_",CTLName,".png")),res=600,units='in',height=6,width=6)
par(mfrow=c(4,4),mar=c(3,3,0,0),oma=c(1,1,1,1))
PlotLifeHistory(LenAtAgeN,LenAtAgeS,matureN,matureS,vulnN,vulnS,survSelN,survSelS,WeightAtAgeN,
WeightAtAgeS,MovementN,MovementS,NatMs,NatMn,VirBioN,VirBioS,RzeroN,RecErrN,steepnessN,steepnessS,
RzeroS,RecErrS,sigmaRn,sigmaRs,HistoricalFn,HistoricalFs,SimYear,MaxAge)
dev.off()
if(echo) cat("\n ####### PLOTTED LIFE HISTORY ####### \n",file=echofile,append=T)
#=========================================================================
# INITALIZE THE POPULATION
#==========================================================================
tempNn <-array(dim=c(InitYear,MaxAge,Nsim))
tempNn[1,,] <-VirInitN
tempNs <-array(dim=c(InitYear,MaxAge,Nsim))
tempNs[1,,] <-VirInitS
tempCatchN <-matrix(ncol=InitYear,nrow=Nsim)
tempCatchS <-matrix(ncol=InitYear,nrow=Nsim)
tempRecN <-matrix(ncol=InitYear,nrow=Nsim)
tempRecS <-matrix(ncol=InitYear,nrow=Nsim)
tempCatchAtAgeN <-array(dim=c(InitYear,MaxAge,Nsim))
tempCatchAtAgeS <-array(dim=c(InitYear,MaxAge,Nsim))
#==Make a matrix of past Fs based on a sd and an input time series==
HistoricalFsInit <-matrix(HistoricalFs,ncol=InitYear,nrow=Nsim,byrow=T)
HistoricalFnInit <-matrix(HistoricalFn,ncol=InitYear,nrow=Nsim,byrow=T)
FerrN <-matrix(rnorm(InitYear*Nsim,1,PastFsdN),ncol=InitYear)
#cat("\n",file=echofile,append=T)
if(echo) cat("\n # FerrN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,FerrN,file=echofile,append=T)
FerrS <-matrix(rnorm(InitYear*Nsim,1,PastFsdS),ncol=InitYear)
HistoricalFsIn <-HistoricalFsInit*FerrN
#cat("\n",file=echofile,append=T)
if(echo) cat("\n # HistoricalFsIn \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,HistoricalFsIn,file=echofile,append=T)
HistoricalFnIn <-HistoricalFnInit*FerrS
for(k in 1:Nsim)
{
for (j in 2:InitYear)
{
for (i in 2:(MaxAge-1))
{
tempNn[j,i,k] <-tempNn[j-1,i-1,k]*exp(-HistoricalFnIn[k,j]*vulnN[1,i-1])*exp(-NatMn[j])
tempNs[j,i,k] <-tempNs[j-1,i-1,k]*exp(-HistoricalFsIn[k,j]*vulnS[1,i-1])*exp(-NatMs[j])
}
tempNn[j,MaxAge,k] <-(tempNn[j-1,(MaxAge-1),k])*exp(-HistoricalFnIn[k,j]*vulnN[j,MaxAge])*exp(-NatMn[j])+ tempNn[j-1,MaxAge,k]*exp(-HistoricalFnIn[k,j]*vulnN[j,MaxAge])*exp(-NatMn[j])
tempNs[j,MaxAge,k] <-(tempNs[j-1,(MaxAge-1),k])*exp(-HistoricalFsIn[k,j]*vulnS[j,MaxAge])*exp(-NatMs[j])+ tempNs[j-1,MaxAge,k]*exp(-HistoricalFsIn[k,j]*vulnS[j,MaxAge])*exp(-NatMs[j])
moveFromN <-tempNn[j,,k]*MovementN[j,]
moveFromS <-tempNs[j,,k]*MovementS[j,]
tempNn[j,,k] <-tempNn[j,,k]-moveFromN+moveFromS
tempNs[j,,k] <-tempNs[j,,k]-moveFromS+moveFromN
EggsN <-sum(tempNn[j-1,,k]*matureN[j,]*WeightAtAgeN[j,])
EggsS <-sum(tempNs[j-1,,k]*matureS[j,]*WeightAtAgeS[j,])
tempNn[j,1,k] <-Recruitment(EggsIN=EggsN,steepnessIN=steepnessN[j],RzeroIN=RzeroN[j],RecErrIN=RecErrN[k,j],recType="BH",NatMin=NatMn[j],
vulnIN=vulnN[j,],matureIN=matureN[j,],weightIN=WeightAtAgeN[j,],LenAtAgeIN=LenAtAgeN[j,],MaxAge=MaxAge,sigmaRin=sigmaRn[j])
tempNs[j,1,k] <-Recruitment(EggsIN=EggsS,steepnessIN=steepnessS[j],RzeroIN=RzeroS[j],RecErrIN=RecErrS[k,j],recType="BH",NatMin=NatMs[j],
vulnIN=vulnS[j,],matureIN=matureS[j,],weightIN=WeightAtAgeS[j,],LenAtAgeIN=LenAtAgeS[j,],MaxAge=MaxAge,sigmaRin=sigmaRs[j])
tempRecN[k,j] <-tempNn[j,1,k]
tempRecS[k,j] <-tempNs[j,1,k]
tempCatchAtAgeN[j,,k]<-((vulnN[j,]*HistoricalFnIn[k,j])/(vulnN[j,]*HistoricalFnIn[k,j]+NatMn[j])) * (1-exp(-(vulnN[j,]*HistoricalFnIn[k,j]+NatMn[j]))) * tempNn[j-1,,k]
tempCatchN[k,j] <-sum(tempCatchAtAgeN[j,,k]*WeightAtAgeN[j,])
tempCatchAtAgeS[j,,k]<-((vulnS[j,]*HistoricalFsIn[k,j])/(vulnS[j,]*HistoricalFsIn[k,j]+NatMs[j])) * (1-exp(-(vulnS[j,]*HistoricalFsIn[k,j]+NatMs[j]))) * tempNs[j-1,,k]
tempCatchS[k,j] <-sum(tempCatchAtAgeS[j,,k]*WeightAtAgeS[j,])
}
}
#===============================================================
# BEGIN SIMULATION OF ASSESSMENT AND HARVEST
#===============================================================
#==tempNn and tempNs from above are the starting points
if(echo) cat("\n\n",file=echofile,append=T)
if(echo) cat("## Full Population \n",file=echofile,append=T)
projNn <-array(dim=c(SimYear,MaxAge,Nsim))
projNs <-array(dim=c(SimYear,MaxAge,Nsim))
projCatchAtAgeN <-array(dim=c(SimYear,MaxAge,Nsim))
projCatchAtAgeS <-array(dim=c(SimYear,MaxAge,Nsim))
projCatchN <-matrix(nrow=Nsim,ncol=SimYear)
projCatchS <-matrix(nrow=Nsim,ncol=SimYear)
projSSBn <-matrix(nrow=Nsim,ncol=SimYear)
projSSBs <-matrix(nrow=Nsim,ncol=SimYear)
projExpBn <-matrix(nrow=Nsim,ncol=SimYear)
projExpBs <-matrix(nrow=Nsim,ncol=SimYear)
projSurvN <-matrix(nrow=Nsim,ncol=SimYear)
projSurvS <-matrix(nrow=Nsim,ncol=SimYear)
projRecN <-matrix(nrow=Nsim,ncol=SimYear)
projRecS <-matrix(nrow=Nsim,ncol=SimYear)
projFmortN <-matrix(nrow=Nsim,ncol=SimYear)
projFmortS <-matrix(nrow=Nsim,ncol=SimYear)
projCatLenFreqN <-array(dim=c(SimYear,LengthBinN,Nsim))
projCatLenFreqS <-array(dim=c(SimYear,LengthBinN,Nsim))
projSurvLenFreqN <-array(dim=c(SimYear,LengthBinN,Nsim))
projSurvLenFreqS <-array(dim=c(SimYear,LengthBinN,Nsim))
#=storage for the true quantities
trueRecN <-matrix(ncol=SimYear,nrow=Nsim)
trueCatchN <-matrix(ncol=SimYear,nrow=Nsim)
trueRecS <-matrix(ncol=SimYear,nrow=Nsim)
trueCatchS <-matrix(ncol=SimYear,nrow=Nsim)
trueFmortN <-matrix(ncol=SimYear,nrow=Nsim)
for(x in 1:Nsim)
trueFmortN[x,1:InitYear]<-HistoricalFn
trueFmortS <-matrix(ncol=SimYear,nrow=Nsim)
for(x in 1:Nsim)
trueFmortS[x,1:InitYear]<-HistoricalFn
trueSpbioN <-matrix(ncol=SimYear,nrow=Nsim)
trueSurvIndN <-matrix(ncol=SimYear,nrow=Nsim)
trueCPUEindN <-matrix(ncol=SimYear,nrow=Nsim)
trueSpbioS <-matrix(ncol=SimYear,nrow=Nsim)
trueSurvIndS <-matrix(ncol=SimYear,nrow=Nsim)
trueCPUEindS <-matrix(ncol=SimYear,nrow=Nsim)
trueF35 <-rep(0,SimYear)
trueSBPR35 <-rep(0,SimYear)
trueB35 <-matrix(ncol=SimYear,nrow=Nsim)
#========================================================================
# FIND REFERENCE POINTS FOR THE POPULATION
#========================================================================
if(PlotYieldCurve==1)
{
pdf(file.path(PlotFolder,paste0("Plot_init_YieldCurve_",CTLName,".pdf")))
SearchFmort <-seq(0.01,3*NatMn[1],(NatMn[1]-0.01)/100)
SearchYield <-rep(0,length(SearchFmort))
SearchBiomass <-rep(0,length(SearchFmort))
print("Plotting yield curve")
for(p in 1:length(SearchFmort))
{
tempOut<-ProjPopDym(fmortN=SearchFmort[p],MaxAge=MaxAge,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,
LenAtAgeS=LenAtAgeS,sigmaRn=sigmaRn,sigmaRs=sigmaRs,RefYear=SimYear,
MovementN=MovementN,MovementS=MovementS)
SearchYield[p]<-tempOut[[1]]
SearchBiomass[p]<-tempOut[[2]]
}
print("Yield curve plotted")
par(mfrow=c(1,2))
plot(SearchYield~SearchFmort)
plot(SearchYield~SearchBiomass)
dev.off()
}
#==============================================================
# calculate the catch proportion at (length) for assessment
#==============================================================
tempCatAtLenN<-array(dim=c(ncol(LenAtAgeN),LengthBinN,Nsim))
tempCatAtLenS<-array(dim=c(ncol(LenAtAgeN),LengthBinN,Nsim))
for(x in 1:Nsim)
for(y in 2:InitYear)
{
#==make length frequencies for catch==
for(w in 1:ncol(LenAtAgeN))
{
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeN[y,w],sd=GrowthSDn[y])
ProbN<-probtemp/sum(probtemp)
tempCatAtLenN[w,,x]<-tempCatchAtAgeN[y,w,x]*ProbN
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeS[y,w],sd=GrowthSDs[y])
ProbS<-probtemp/sum(probtemp)
tempCatAtLenS[w,,x]<-tempCatchAtAgeS[y,w,x]*ProbS
}
projCatLenFreqN[y,,x]<-apply(tempCatAtLenN[,,x],2,sum)
projCatLenFreqS[y,,x]<-apply(tempCatAtLenS[,,x],2,sum)
}
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat(paste0("# projCatLenFreqN #array! \n"),file=echofile,append=T)
if(echo) for(n in 1:Nsim) {
cat(paste0("# projCatLenFreqN sim: ",n," \n"),file=echofile,append=T)
write.table(row.names=F,col.names=F,projCatLenFreqN[,,n],file=echofile,append=T)
}
if(AssessmentData==1)
{
pdf(file.path(PlotFolder,paste0("CatchLengthProportions_",CTLName,".pdf")))
par(mfrow=c(ceiling(sqrt(InitYear)),ceiling(sqrt(InitYear))),mar=c(.1,.1,.1,.1))
for(i in 2:InitYear)
barplot(rbind(projCatLenFreqN[i,,1],projCatLenFreqS[i,,1]),beside=T,xaxt='n',yaxt='n')
dev.off()
}
#==============================================================
# calculate the survey proportion at length for assessment
#==============================================================
tempSurvAtLenN<-array(dim=c(ncol(LenAtAgeN),LengthBinN,Nsim))
tempSurvAtLenS<-array(dim=c(ncol(LenAtAgeN),LengthBinN,Nsim))
for(x in 1:Nsim)
for(y in 2:InitYear)
{
#==make length frequencies for s==
for(w in 1:ncol(LenAtAgeS))
{
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeN[y,w],sd=GrowthSDn[y])
ProbN<-probtemp/sum(probtemp)
tempSurvAtLenN[w,,x]<-tempNn[y,w,x]*survSelN[y,w]*ProbN
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeS[y,w],sd=GrowthSDs[y])
ProbS<-probtemp/sum(probtemp)
tempSurvAtLenS[w,,x]<-tempNs[y,w,x]*survSelS[y,w]*ProbS
}
projSurvLenFreqN[y,,x]<-apply(tempSurvAtLenN[,,x],2,sum)
projSurvLenFreqS[y,,x]<-apply(tempSurvAtLenS[,,x],2,sum)
}
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat(paste0("# projSurvLenFreqN #array!"," \n"),file=echofile,append=T)
if(echo) for(n in 1:Nsim) {
cat(paste0("# projSurvLenFreqN sim:",n," \n"),file=echofile,append=T)
write.table(row.names=F,col.names=F,projSurvLenFreqN[,,n],file=echofile,append=T)
}
if(AssessmentData==1)
{
pdf(file.path(PlotFolder,paste0("SurveyLengthProportions_",CTLName,".pdf")))
par(mfrow=c(ceiling(sqrt(InitYear)),ceiling(sqrt(InitYear))),mar=c(.1,.1,.1,.1))
for(i in 2:InitYear)
barplot(rbind(projSurvLenFreqN[i,,1],projSurvLenFreqS[i,,1]),beside=T,xaxt='n',yaxt='n')
dev.off()
}
#==transfer historical time series from above
for(x in 1:Nsim)
{
projNn[1:InitYear,,x] <-tempNn[,,x]
projNs[1:InitYear,,x] <-tempNs[,,x]
projCatchN[x,1:InitYear] <-tempCatchN[x,]
projCatchS[x,1:InitYear] <-tempCatchS[x,]
projRecN[x,1:InitYear] <-tempRecN[x,]
projRecS[x,1:InitYear] <-tempRecS[x,]
projFmortN[x,1:InitYear] <-HistoricalFnInit[x,]
projFmortS[x,1:InitYear] <-HistoricalFsInit[x,]
projSurvN[x,1:InitYear] <-apply(tempNn[,,x]*survSelN[1:InitYear,]*WeightAtAgeN[1:InitYear,],1,sum)
projSurvS[x,1:InitYear] <-apply(tempNs[,,x]*survSelS[1:InitYear,]*WeightAtAgeS[1:InitYear,],1,sum)
for(y in 1:InitYear)
{
projSSBn[x,y] <-sum(projNn[y,,x]*matureN[y,]*WeightAtAgeN[y,])
projSSBs[x,y] <-sum(projNs[y,,x]*matureS[y,]*WeightAtAgeS[y,])
projExpBn[x,y] <-sum(projNn[y,,x]*vulnN[y,]*WeightAtAgeN[y,])
projExpBs[x,y] <-sum(projNs[y,,x]*vulnS[y,]*WeightAtAgeS[y,])
}
}
if(echo) cat(paste0("# projNn #array!"," \n"),file=echofile,append=T)
if(echo) for(n in 1:Nsim) {
cat(paste0("# projNn sim:",n," \n"),file=echofile,append=T)
write.table(row.names=F,col.names=F,projNn[,,n],file=echofile,append=T)
}
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat("# projRecN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,projRecN,file=echofile,append=T)
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat("# projFmortN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,projFmortN,file=echofile,append=T)
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat("# projSSBn \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,projSSBn,file=echofile,append=T)
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat("# projExpBn \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,projExpBn,file=echofile,append=T)
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat("# projSurvN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,projSurvN,file=echofile,append=T)
projCatchN[,1] <-0
if(echo) cat("\n",file=echofile,append=T)
if(echo) cat("# projCatchN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,projCatchN,file=echofile,append=T)
projCatchS[,1] <-0
#==fill in true storage arrays
for(x in 1:Nsim)
{
trueRecN[x,1:InitYear] <-tempRecN[x,]
trueCatchN[x,1:InitYear] <-tempCatchN[x,]
trueSpbioN[x,1:InitYear] <-projSSBn[x,1:InitYear]
trueSurvIndN[x,1:InitYear] <-projSurvN[x,1:InitYear]
trueCPUEindN[x,1:InitYear] <-projExpBn[x,1:InitYear]
trueRecS[x,1:InitYear] <-tempRecS[x,]
trueCatchS[x,1:InitYear] <-tempCatchS[x,]
trueSpbioS[x,1:InitYear] <-projSSBs[x,1:InitYear]
trueSurvIndS[x,1:InitYear] <-projSurvS[x,1:InitYear]
trueCPUEindS[x,1:InitYear] <-projExpBs[x,1:InitYear]
}
CatchErrorN <-matrix(rnorm(Nsim*SimYear,0,CatchCVn),nrow=Nsim,ncol=SimYear)
if(echo) cat("\n # CatchErrorN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,CatchErrorN,file=echofile,append=T)
CatchErrorS <-matrix(rnorm(Nsim*SimYear,0,CatchCVs),nrow=Nsim,ncol=SimYear)
CatchAssessN <-projCatchN*exp(CatchErrorN-(CatchCVn^2/2))
if(echo) cat("\n # CatchAssessN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,CatchAssessN,file=echofile,append=T)
CatchAssessS <-projCatchS*exp(CatchErrorS-(CatchCVs^2/2))
CPUEErrorN <-matrix(rnorm(Nsim*SimYear,0,IndexCVn),nrow=Nsim,ncol=SimYear)
if(echo) cat("\n # CPUEErrorN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,CPUEErrorN,file=echofile,append=T)
CPUEErrorS <-matrix(rnorm(Nsim*SimYear,0,IndexCVs),nrow=Nsim,ncol=SimYear)
CPUEAssessN <-projExpBn*exp(CPUEErrorN-(IndexCVn^2/2))
if(echo) cat("\n # CPUEAssessN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,CPUEAssessN,file=echofile,append=T)
CPUEAssessS <-projExpBs*exp(CPUEErrorS-(IndexCVs^2/2))
SurvErrorN <-matrix(rnorm(Nsim*SimYear,0,IndexCVn),nrow=Nsim,ncol=SimYear)
if(echo) cat("\n # SurvErrorN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,SurvErrorN,file=echofile,append=T)
SurvErrorS <-matrix(rnorm(Nsim*SimYear,0,IndexCVs),nrow=Nsim,ncol=SimYear)
SurvAssessN <-projSurvN*exp(SurvErrorN-(IndexCVn^2/2))
if(echo) cat("\n # SurveAssessN \n",file=echofile,append=T)
if(echo) write.table(row.names=F,col.names=F,SurvAssessN,file=echofile,append=T)
SurvAssessS <-projSurvS*exp(SurvErrorS-(IndexCVs^2/2))
if(echo) cat(paste0("Input Ended: ",Sys.time(),"\n"),file=echofile,append=T)
#==storage for management and assessment quantities
FMSY <-matrix(nrow=Nsim,ncol=SimYear)
BMSY <-matrix(nrow=Nsim,ncol=SimYear)
TAC <-matrix(nrow=Nsim,ncol=SimYear)
trueTAC <-matrix(nrow=Nsim,ncol=SimYear)
trueOFL <-matrix(nrow=Nsim,ncol=SimYear)
CurBio <-matrix(nrow=Nsim,ncol=SimYear)
EstBio <-matrix(nrow=Nsim,ncol=SimYear)
Converge <-matrix(nrow=Nsim,ncol=SimYear)
est_init_B <-matrix(nrow=Nsim,ncol=SimYear)
#============================================
# calculate 'true' reference points
# figure out a way to determine when to do this for all sims and when for just 1
# as long as population parametesr are from the same control file, this should be the same
# for all populations regardless of fishing history or recruitment
#===================================
FindFMSYin<-1
if(FindFMSYin==1)
{
#==Find reference points for the initial year of projection
x <-0.3
FmsyOut <-nlminb(x,FindFMSY,MaxAge=MaxAge,VirInitN=VirInitN,VirInitS=VirInitS,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,LenAtAgeS=LenAtAgeS,
sigmaRn=sigmaRn,sigmaRs=sigmaRs,RefYear=InitYear,MovementN=MovementN,MovementS=MovementS)
trueFMSYbegin <-FmsyOut$par
trueUMSYbegin <- 1-(exp(-trueFMSYbegin))
trueBMSYbegin <-ProjPopDym(fmortN=trueFMSYbegin,MaxAge=MaxAge,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,LenAtAgeS=LenAtAgeS,
sigmaRn=sigmaRn,sigmaRs=sigmaRs,RefYear=InitYear,MovementN=MovementN,MovementS=MovementS)[[2]]
trueMSYbegin <- -FmsyOut$objective
#==find the reference points associated with the final year of projection
FmsyOut <-nlminb(x,FindFMSY,MaxAge=MaxAge,VirInitN=VirInitN,VirInitS=VirInitS,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,LenAtAgeS=LenAtAgeS,
sigmaRn=sigmaRn,sigmaRs=sigmaRs,RefYear=SimYear,MovementN=MovementN,MovementS=MovementS)
trueFMSYend <-FmsyOut$par
trueUMSYend <- 1-(exp(-trueFMSYend))
trueBMSYend <-ProjPopDym(fmortN=trueFMSYend,MaxAge=MaxAge,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,
LenAtAgeS=LenAtAgeS,sigmaRn=sigmaRn,sigmaRs=sigmaRs,RefYear=SimYear,
MovementN=MovementN,MovementS=MovementS)[[2]]
trueMSYend <- -FmsyOut$objective
}
#======================================
# calculate reference point proxies
# MIGHT NEED TWO OF THESE
#======================================
ConstRec <-100000
outsF35in <-FindF35(MaxAge=MaxAge,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,
LenAtAgeS=LenAtAgeS,inRec=ConstRec,RefYear=InitYear,MovementN=MovementN,MovementS=MovementS)
trueSBPR35[InitYear]<-outsF35in[[2]]
trueF35[InitYear] <-outsF35in[[1]]
if(Nsim>1) trueB35[,InitYear] <-trueSBPR35[InitYear]*apply(trueRecN[,start_assessment:SimYear],1,mean,na.rm=T)
if(Nsim==1) trueB35[,InitYear] <-trueSBPR35[InitYear]*mean(trueRecN[,start_assessment:SimYear],na.rm=T)
#=================================================================================
#==BEGIN PROJECTIONS===========================================================
#=================================================================================
for(z in 1:Nsim)
{
for(y in (InitYear+1):SimYear)
{
if(z==1)
{
# this calculates the F35 and SBPR35
# make it so this only does it if the things that change F35 and SBPR35 change
# this can be done with an if statement that calls a function that checksif the values are all the same (or columns are all the same)
outsF35in <-FindF35(MaxAge=MaxAge,vulnN=vulnN,
vulnS=vulnS,NatMn=NatMn,NatMs=NatMs,matureN=matureN,matureS=matureS,
WeightAtAgeN=WeightAtAgeN, WeightAtAgeS=WeightAtAgeS,steepnessN=steepnessN,
steepnessS=steepnessS,RzeroN=RzeroN,RzeroS=RzeroS,LenAtAgeN=LenAtAgeN,
LenAtAgeS=LenAtAgeS,inRec=ConstRec,RefYear=y,MovementN=MovementN,MovementS=MovementS)
trueSBPR35[y]<-outsF35in[[2]]
trueF35[y] <-outsF35in[[1]]
}
#==pull data for assessment
CatchDataN <- CatchAssessN[z,!is.na(CatchAssessN[z,])]
CPUEDataN <- CPUEAssessN[z,!is.na(CPUEAssessN[z,])]
SurvDataN <- SurvAssessN[z,!is.na(SurvAssessN[z,])]
CatchDataS <- CatchAssessS[z,!is.na(CatchAssessS[z,])]
CPUEDataS <- CPUEAssessS[z,!is.na(CPUEAssessS[z,])]
SurvDataS <- SurvAssessS[z,!is.na(SurvAssessS[z,])]
#=================================================================================
#===ASSESS THE STOCK, CHOOSE A METHOD IN CTL FILE=================================
#=================================================================================
CTLfolder <- file.path(MSEdir,CTLName)
if(!dir.exists(CTLfolder)) dir.create(CTLfolder)
#==Production model==
if(AssessmentType == 1)
{
if(y==(InitYear+1))
{
x <-c((InitBzeroMod*(VirBioN)),InitGrowthRate)
if(!InitShape%in%c(0,1))
x <-c(x,InitShape)
if(estInit==1)
x <-c(x,InitBioProd)
}
logx<-suppressWarnings(log(x))
if(y>(InitYear+1))
logx <-outs$par
inCatch <-CatchDataN[start_assessment:(y-1)]
inCPUE <-CPUEDataN[start_assessment:(y-1)]
outs <-suppressWarnings(nlminb(start=logx,objective=ProdMod,CatchData=inCatch,IndexData=inCPUE,estInit=estInit,InitShape=InitShape))
#outs <- optim(par=x,fn=ProdMod,CatchData=inCatch,IndexData=inCPUE,estInit=estInit)
if(sum(is.na(outs$par))>0) {stop("Production model converged on NaNs. Not really sure why. Try changing your starting values?")}
Converge[z,y]<-outs$convergence
PredBio <-ProdModPlot(exp(outs$par),inCatch,inCPUE,plots=EstimationPlots,estInit=estInit,InitShape=InitShape)
if(!InitShape%in%c(0,1)) {
FMSY[z,y] <-(exp(outs$par[2])/exp(outs$par[3]))*(1-(1/(exp(outs$par[3])+1)))
BMSY[z,y] <-exp(outs$par[1])*(exp(outs$par[3])+1)^(-1/exp(outs$par[3]))
if(estInit==1)
est_init_B[z,y]<-exp(outs$par[4])
}
if(InitShape%in%c(0,1)) {
if(InitShape==0)
InitShape<-1E-10
FMSY[z,y] <-(exp(outs$par[2])/InitShape)*(1-(1/(InitShape+1)))
BMSY[z,y] <-exp(outs$par[1])*(InitShape+1)^(-1/InitShape)
if(InitShape==1E-10)
InitShape<-0
if(estInit==1)
est_init_B[z,y]<-exp(outs$par[3])
}
MSY <-BMSY[z,y]*FMSY[z,y]
CurBio[z,y]<-PredBio[length(PredBio)-1]
inBio <-projExpBn[z,y-1]
if(CurBio[z,y]<0) CurBio[z,y]<-0
TAC[z,y] <-HarvestControlRule(FMSY=FMSY[z,y],BMSY=BMSY[z,y],ExploitBio=CurBio[z,y],SpawnBio=CurBio[z,y],
alpha=HCalphaN,beta=HCbetaN,HarvestControl=HarvestControlN,ConstantCatch=ConstantCatchN,
ConstantF=ConstantFn)
trueTAC[z,y]<-HarvestControlRule(FMSY=trueFMSYend,BMSY=trueBMSYend,ExploitBio=inBio,SpawnBio=inBio,
alpha=HCalphaN,beta=HCbetaN,HarvestControl=HarvestControlN,ConstantCatch=ConstantCatchN,
ConstantF=ConstantFn)
if(EstimationPlots==1 & AssessmentType == 1)
{
lines(BMSY[z,],col=3,lty=3)
legend("topright",bty='n',col=c(NA,1,2,3),pch=c(15,NA,NA,NA),lty=c(NA,1,2,3),
legend=c("Index obs","Index est","Catch","BMSY"))
}
if(y==SimYear)
{
EstBio[z,start_assessment:ncol(EstBio)]<-PredBio
}
if(y==SimYear & z==Nsim)
{
inFile<-file.path(CTLfolder,"ProdOutputs.csv")
file.create(inFile)
cat("#True quantities","\n",file=inFile)
cat("#","\n",file=inFile,append=TRUE)
cat("#true Catch","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(trueCatchN[m,],"\n",file=inFile,append=TRUE)
cat("#true fishing mortality","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(trueFmortN[m,],"\n",file=inFile,append=TRUE)
cat("#true spawning biomass","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(trueSpbioN[m,],"\n",file=inFile,append=TRUE)
cat("#true CPUE","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(trueCPUEindN[m,],"\n",file=inFile,append=TRUE)
cat("#true total allowable catch","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(trueTAC[m,],"\n",file=inFile,append=TRUE)
cat("#true BMSY begin","\n",file=inFile,append=TRUE)
cat(trueBMSYbegin,"\n",file=inFile,append=TRUE)
cat("#true BMSY end","\n",file=inFile,append=TRUE)
cat(trueBMSYend,"\n",file=inFile,append=TRUE)
cat("#true FMSY begin","\n",file=inFile,append=TRUE)
cat(trueFMSYbegin,"\n",file=inFile,append=TRUE)
cat("#true FMSY end","\n",file=inFile,append=TRUE)
cat(trueFMSYend,"\n",file=inFile,append=TRUE)
cat("#est cpue ind","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(EstBio[m,],"\n",file=inFile,append=TRUE)
cat("#est spawning biomass by year (terminal year)","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(CurBio[m,],"\n",file=inFile,append=TRUE)
cat("#est total allowable catch","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(TAC[m,],"\n",file=inFile,append=TRUE)
cat("#est BMSY","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(BMSY[m,],"\n",file=inFile,append=TRUE)
cat("#est FMSY","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(FMSY[m,],"\n",file=inFile,append=TRUE)
if(estInit==1)
{
cat("#est initial biomass","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(est_init_B[m,],"\n",file=inFile,append=TRUE)
}
cat("#data cpue","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(CPUEDataN,"\n",file=inFile,append=TRUE)
cat("#convergence","\n",file=inFile,append=TRUE)
for(m in 1:Nsim)
cat(Converge[m,],"\n",file=inFile,append=TRUE)
}
}
#========================================
#==Projection based model--no assessment
#=======================================
if(AssessmentType == 0)
{
#==make folder if it isn't there
if(!dir.exists(file.path(CTLfolder,z))) dir.create(file.path(CTLfolder,z))
IndSimFolder<-file.path(CTLfolder,z,y)
if(!dir.exists(IndSimFolder)) dir.create(IndSimFolder)
#==write the true values
TQfile <- file.path(IndSimFolder,"TrueQuantities.DAT")
file.create(TQfile)
#==calculate the TAC based on a SRR or a proxy, depending on data
#==allow this to be a user defined function of any of the things that can be input
TAC[z,y] <-trueSpbio[z,y-1]*(1-exp(-ConstantF))
}
#======================================================
#==Age-structured assessment model for setting the TAC
#======================================================
if(AssessmentType == 2)
{
#==make folder if it isn't there
if(!dir.exists(file.path(CTLfolder,z))) dir.create(file.path(CTLfolder,z))
IndSimFolder<-file.path(CTLfolder,z,y)
if(!dir.exists(IndSimFolder)) dir.create(IndSimFolder)
#==copy the .exe into it (where does this come from? github?)
#==Better way of doing this??
execfile <- get_exec(SimAssExec)
copied <- file.copy(from=execfile,to=IndSimFolder,overwrite=TRUE)
if(!copied) stop(paste0("Executable failed to copy from ",execfile,". Might want to check what's up."))
#==write the true values
#==Probably don't need this if it is stored in the MSE object
TQfile <- file.path(IndSimFolder,"TrueQuantities.DAT")
if(!file.exists(TQfile)) file.create(TQfile)
# TQfile now written at the end after the assessment -LQ 11/10/18
#==.CTL file !!!!!!!!!!!!!!!!!how should we deal with space?!!!!!!!!!!!!!!!!!!!!!!
CTLfile <- file.path(IndSimFolder,"SimAss.CTL")
if(!file.exists(CTLfile)) file.create(CTLfile)
cat("#Simulated assessment control file","\n",file=CTLfile)
cat("#","\n",file=CTLfile,append=TRUE)
cat("#weight parameters","\n",file=CTLfile,append=TRUE)
cat(c(alphaN[1],betaN[1]) ,"\n",file=CTLfile,append=TRUE)
cat("#how many years to average over timevarying process for projectsion","\n",file=CTLfile,append=TRUE)
cat(ProjectTimeVary,"\n",file=CTLfile,append=TRUE)
cat("#Estimate growth? if positive, indicates phase, no estimate if negative","\n",file=CTLfile,append=TRUE)
cat(EstGrowthK ,"\n",file=CTLfile,append=TRUE)
cat("#growth time-varying?","\n",file=CTLfile,append=TRUE)
cat(TimeVaryGrowthK ,"\n",file=CTLfile,append=TRUE)
cat("#Estimate growth? if positive, indicates phase, no estimate if negative","\n",file=CTLfile,append=TRUE)
cat(EstLinf ,"\n",file=CTLfile,append=TRUE)
cat("#growth time-varying?","\n",file=CTLfile,append=TRUE)
cat(TimeVaryLinf ,"\n",file=CTLfile,append=TRUE)
cat("#length parameters","\n",file=CTLfile,append=TRUE)
cat(c(VonKn[1],LinfN[1],t0n[1]) ,"\n",file=CTLfile,append=TRUE)
cat("#estimate natural mortality? if positive, number indicates phase, if negative, not estimated","\n",file=CTLfile,append=TRUE)
cat(EstM,"\n",file=CTLfile,append=TRUE)
cat("#natural mortality (used as starting value if estimated)","\n",file=CTLfile,append=TRUE)
cat(NatMn[length(NatMn)],"\n",file=CTLfile,append=TRUE)
cat("#Allow M to vary over time? 1 = yes, 0 = no","\n",file=CTLfile,append=TRUE)
cat(TimeVaryM,"\n",file=CTLfile,append=TRUE)
cat("#selectivity 50% time-varying?","\n",file=CTLfile,append=TRUE)
cat(TimeVarySel50 ,"\n",file=CTLfile,append=TRUE)
cat("#selectivity 95% time-varying?","\n",file=CTLfile,append=TRUE)
cat(TimeVarySel95 ,"\n",file=CTLfile,append=TRUE)
cat("#maturity","\n",file=CTLfile,append=TRUE)
cat(c(mat50n[length(mat50n)],mat95n[length(mat95n)]) ,"\n",file=CTLfile,append=TRUE)
cat("#steepness","\n",file=CTLfile,append=TRUE)
cat(steepnessN[length(steepnessN)],"\n",file=CTLfile,append=TRUE)
cat("#small number","\n",file=CTLfile,append=TRUE)
cat(SmallNum,"\n",file=CTLfile,append=TRUE)
cat("#initial smoothness weight","\n",file=CTLfile,append=TRUE)
cat(InitSmooth,"\n",file=CTLfile,append=TRUE)
cat("#Fishery Independent data available?","\n",file=CTLfile,append=TRUE)
cat(FisheryIndepenDat,"\n",file=CTLfile,append=TRUE)
cat("#Length at age standard deviation","\n",file=CTLfile,append=TRUE)
cat(GrowthSDn[length(steepnessN)],"\n",file=CTLfile,append=TRUE)
cat("#smoothness on fishing mortality penalty","\n",file=CTLfile,append=TRUE)
cat(FmortPen,"\n",file=CTLfile,append=TRUE)
cat("#smoothness on recruit penalty","\n",file=CTLfile,append=TRUE)
cat(RecruitPen,"\n",file=CTLfile,append=TRUE)
cat("#smoothness on time-varying natural mortality penalty","\n",file=CTLfile,append=TRUE)
cat(Mpenalty,"\n",file=CTLfile,append=TRUE)
cat("#smoothness on time-varying growth penalty","\n",file=CTLfile,append=TRUE)
cat(Growthpenalty,"\n",file=CTLfile,append=TRUE)
cat("#smoothness on time-varying selelctivity penalty","\n",file=CTLfile,append=TRUE)
cat(SelPenalty,"\n",file=CTLfile,append=TRUE)
cat("#harvest control rule selected","\n",file=CTLfile,append=TRUE)
cat(HarvestControlN,"\n",file=CTLfile,append=TRUE)
cat("#Constantcatch if harvest control ==2, set catch level","\n",file=CTLfile,append=TRUE)
cat(ConstantCatchN,"\n",file=CTLfile,append=TRUE)
cat("#Constantcatch if harvest control ==3, set fmort level","\n",file=CTLfile,append=TRUE)
cat(ConstantFn,"\n",file=CTLfile,append=TRUE)
cat("#40 10 parameters for harvest control 4","\n",file=CTLfile,append=TRUE)
cat(HCalphaN,"\n",file=CTLfile,append=TRUE)
cat("#40 10 parameters for harvest control 4","\n",file=CTLfile,append=TRUE)
cat(HCbetaN,"\n",file=CTLfile,append=TRUE)
cat("#true FMSY using final year parameters","\n",file=CTLfile,append=TRUE)
cat(trueFMSYend,"\n",file=CTLfile,append=TRUE)
#==write .PIN file to get initial values close for estimation
PINfile <- file.path(IndSimFolder,"SimAss.PIN")
if(!file.exists(PINfile)) file.create(PINfile)
cat("# Simulated assessment pin file","\n",file=PINfile)
cat("#","\n",file=PINfile,append=TRUE)
InputSurvSel50N<-selAtAgeFunc(surv50n,VonKn,LinfN,t0n)
InputSurvSel95N<-selAtAgeFunc(surv95n,VonKn,LinfN,t0n)
cat("# srv_sel50","\n",file=PINfile,append=TRUE)
cat(mean(InputSurvSel50N,na.rm=T)*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# srv_sel95","\n",file=PINfile,append=TRUE)
cat(mean(InputSurvSel95N,na.rm=T)*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# stNatLen","\n",file=PINfile,append=TRUE)
cat(log(VirInitN)*rnorm(length(log(VirInitN)),1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# log_avg_fmort_dir","\n",file=PINfile,append=TRUE)
cat(log(mean(trueFmortN,na.rm=T))*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# fmort_dir_dev","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
cat("# mean_log_rec","\n",file=PINfile,append=TRUE)
cat(log(mean(trueRecN,na.rm=T))*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# rec_dev","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
cat("# log_avg_NatM","\n",file=PINfile,append=TRUE)
if(EstM > 0) {
cat(log(mean(NatMn,na.rm=T))*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
}
if(EstM <= 0) {
cat(log(mean(NatMn,na.rm=T)),"\n",file=PINfile,append=TRUE)
}
cat("# NatM_dev","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
cat("# log_avg_GrowthK","\n",file=PINfile,append=TRUE)
if(EstGrowthK > 0) {
cat(log(mean(VonKn,na.rm=T))*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
}
if(EstGrowthK <= 0) {
cat(log(mean(VonKn,na.rm=T)),"\n",file=PINfile,append=TRUE)
}
cat("# log_avg_Linf","\n",file=PINfile,append=TRUE)
if(EstLinf > 0) {
cat(log(mean(LinfN,na.rm=T))*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
}
if(EstLinf < 0) {
cat(log(mean(LinfN,na.rm=T)),"\n",file=PINfile,append=TRUE)
}
cat("# GrowthK_dev","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
cat("# Linf_dev","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
InputSel50N<-selAtAgeFunc(sel50n,VonKn,LinfN,t0n)
InputSel95N<-selAtAgeFunc(sel95n,VonKn,LinfN,t0n)
cat("# SelPars50","\n",file=PINfile,append=TRUE)
cat(mean(InputSel50N,na.rm=T)*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# SelPars95","\n",file=PINfile,append=TRUE)
cat(mean(InputSel95N,na.rm=T)*rnorm(1,1,InitValSig),"\n",file=PINfile,append=TRUE)
cat("# SelPars_dev50","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
cat("# SelPars_dev95","\n",file=PINfile,append=TRUE)
cat(rep(0,y-start_assessment),"\n",file=PINfile,append=TRUE)
#=======================================
#==.DAT file
#=======================================
DATfile <- file.path(IndSimFolder,"SimAss.DAT")
if(!file.exists(DATfile)) file.create(DATfile)
cat("#Simulated generalized assessment","\n",file=DATfile)
cat("#","\n",file=DATfile,append=TRUE)
cat("#start/end year","\n",file=DATfile,append=TRUE)
cat(c(start_assessment,(y-1)),"\n",file=DATfile,append=TRUE)
#==aggregate the length frequencies from the two areas
cat("#number of ages","\n",file=DATfile,append=TRUE)
cat(MaxAge,"\n",file=DATfile,append=TRUE)
cat("#ages","\n",file=DATfile,append=TRUE)
cat(seq(1,MaxAge),"\n",file=DATfile,append=TRUE)
cat("#Number of length bins","\n",file=DATfile,append=TRUE)
cat(LengthBinN,"\n",file=DATfile,append=TRUE)
#==aggregate the index of abundance from two areas
cat("#CPUE years","\n",file=DATfile,append=TRUE)
cat((y-start_assessment),"\n",file=DATfile,append=TRUE)
cat("#years for CPUE","\n",file=DATfile,append=TRUE)
cat(seq(start_assessment,(y-1)),"\n",file=DATfile,append=TRUE)
cat("#survey years","\n",file=DATfile,append=TRUE)
cat((y-start_assessment),"\n",file=DATfile,append=TRUE)
cat("#years for survey","\n",file=DATfile,append=TRUE)
cat(seq(start_assessment,(y-1)),"\n",file=DATfile,append=TRUE)
cat("#length freq from catch years","\n",file=DATfile,append=TRUE)
cat((y-start_assessment),"\n",file=DATfile,append=TRUE)
cat("#years for length freq catch","\n",file=DATfile,append=TRUE)
cat(seq(start_assessment,(y-1)),"\n",file=DATfile,append=TRUE)
cat("#catch length sample sizes","\n",file=DATfile,append=TRUE)
cat(LenSampleN[seq(start_assessment,(y-1))],"\n",file=DATfile,append=TRUE)
cat("#length freq from survey years","\n",file=DATfile,append=TRUE)
cat((y-start_assessment),"\n",file=DATfile,append=TRUE)
cat("#years for length freq survey","\n",file=DATfile,append=TRUE)
cat(seq(start_assessment,(y-1)),"\n",file=DATfile,append=TRUE)
cat("#survey sample sizes","\n",file=DATfile,append=TRUE)
cat(LenSampleN[seq(start_assessment,(y-1))],"\n",file=DATfile,append=TRUE)
cat("#total survey biomass","\n",file=DATfile,append=TRUE)
cat(SurvDataN[start_assessment:(y-1)],"\n",file=DATfile,append=TRUE)
cat("#survey CV","\n",file=DATfile,append=TRUE)
cat(IndexCVn[seq(start_assessment,(y-1))],"\n",file=DATfile,append=TRUE)
cat("#total CPUE","\n",file=DATfile,append=TRUE)
cat(CPUEDataN[start_assessment:(y-1)],"\n",file=DATfile,append=TRUE)
cat("#CPUE CV","\n",file=DATfile,append=TRUE)
cat(IndexCVn[seq(start_assessment,(y-1))],"\n",file=DATfile,append=TRUE)
cat("#total catch biomass","\n",file=DATfile,append=TRUE)
cat(CatchDataN[start_assessment:(y-1)],"\n",file=DATfile,append=TRUE)
cat("#catch CV","\n",file=DATfile,append=TRUE)
cat(CatchCVn[seq(start_assessment,(y-1))],"\n",file=DATfile,append=TRUE)
cat("#Length Bins","\n",file=DATfile,append=TRUE)
cat(LengthBinsMid,"\n",file=DATfile,append=TRUE)
projCatLenFreqN[is.na(projCatLenFreqN)]<-0
projSurvLenFreqN[is.na(projSurvLenFreqN)]<-0
cat("#catch length counts","\n",file=DATfile,append=TRUE)
for(i in start_assessment:(y-1))
cat(projCatLenFreqN[i,,z],"\n",file=DATfile,append=TRUE)
cat("#survey length counts","\n",file=DATfile,append=TRUE)
for(i in start_assessment:(y-1))
cat(projSurvLenFreqN[i,,z],"\n",file=DATfile,append=TRUE)
setwd(IndSimFolder)
#==run the code
if(os == "windows") {
shell(SimAssComm)
}
if(os != "windows") {
system(SimAssComm)
}
#==delete the .exe
file.remove(file.path(IndSimFolder,SimAssExec))
setwd(CurDir)
#==================================
# PLOT THE ASSESSMENT OUTPUT (if instructed)
#==================================
REP<-readLines(file.path(IndSimFolder,"simass.rep"))
CTL<-readLines(CTLfile)
DAT<-readLines(DATfile)
data2<-scan(file.path(IndSimFolder,"simass.PAR"),what="character",quiet=T)
temp<-grep("OFL",REP)[2]
OFL<-suppressWarnings(as.numeric(unlist(strsplit(REP[temp+1],split=" "))))
#==calculate the TAC based on a SRR or a proxy, depending on data
TAC[z,y] <-OFL
}
#==================================================================================
#==find the F that removes the TAC given the fishery selectivity===================
#==this TAC is not necessarily the 'true' TAC--there is error in the biomass and FMSY
#==================================================================================
minExp<-0.00001
maxExp<-0.99999
inputFs <-0 # NEED TO BE CHANGED TO THE INPUT F FOR POP 2 OR SOME OUTPUT OF HCR
for(p in 1:30)
{
tempExp <-(maxExp+minExp)/2
tempF <- -log(1-tempExp)
tempFn <-tempF
tempFs <-inputFs
finCatLenFn <-((tempFn*vulnN[y,])/(vulnN[y,]*tempFn+NatMn[y])) * (1-exp(-((tempFn*vulnN[y,])+NatMn[y]))) * projNn[y-1,,z]
finCatFn <-sum(finCatLenFn*WeightAtAgeN[y,])
finCatLenFs <-((tempFs*vulnS[y,])/(vulnS[y,]*tempFs+NatMs[y])) * (1-exp(-((tempFs*vulnS[y,])+NatMs[y]))) * projNs[y-1,,z]
finCatFs <-sum(finCatLenFs*WeightAtAgeS[y,])
if(finCatFn<TAC[z,y])
minExp <-tempExp
if(finCatFn>TAC[z,y])
maxExp <-tempExp
initExp <-tempExp
}
ApplyFn <- -log(1-tempExp)
ApplyFs <- inputFs
trueFmortN[z,y] <-ApplyFn
#==calculate true OFL
#==find FOFL given spawning biomass
tempSpBio <-EggsN
FutMort <-trueF35[y]
if(tempSpBio<trueB35[z,y-1])
{
FutMort<-0
if(tempSpBio>HCbetaN*trueB35[z,y-1])
FutMort = trueF35[y]*(tempSpBio/trueB35[z,y-1]-HCalphaN)/(1-HCalphaN)
}
#==find the catch for the FOFL
trueCatchAtAge <-((vulnN[y,]*FutMort)/(vulnN[y,]*FutMort+NatMn[y])) * (1-exp(-(vulnN[y,]*FutMort+NatMn[y]))) * projNn[y-1,,z]
trueCatch <-sum(trueCatchAtAge*WeightAtAgeN[y,])
trueOFL[z,y] <-trueCatch
#===UPDATE DATA FOR ASSESSMENT====================================
projCatchAtAgeN[y,,z] <-((ApplyFn*vulnN[y,])/(vulnN[y,]*ApplyFn+NatMn[y])) * (1-exp(-((ApplyFn*vulnN[y,])+NatMn[y]))) * projNn[y-1,,z]
projCatchN[z,y] <-sum(projCatchAtAgeN[y,,z]*WeightAtAgeN[y,])
projCatchAtAgeS[y,,z] <-((ApplyFs*vulnS[y,])/(vulnS[y,]*ApplyFs+NatMs[y])) * (1-exp(-((ApplyFs*vulnS[y,])+NatMs[y]))) * projNs[y-1,,z]
projCatchS[z,y] <-sum(projCatchAtAgeS[y,,z]*WeightAtAgeS[y,])
#==UPDATE CATCH DATA AND INDEX (survey and CPUE) DATA
projSSBn[z,y] <-sum(projNn[y-1,,z]*matureN[y,]*WeightAtAgeN[y,])
projSSBs[z,y] <-sum(projNs[y-1,,z]*matureS[y,]*WeightAtAgeS[y,])
projExpBn[z,y] <-sum(projNn[y-1,,z]*vulnN[y,]*WeightAtAgeN[y,])
projExpBs[z,y] <-sum(projNs[y-1,,z]*vulnS[y,]*WeightAtAgeS[y,])
projSurvN[z,y] <-sum(projNn[y-1,,z]*survSelN[y,]*WeightAtAgeN[y,])
projSurvS[z,y] <-sum(projNs[y-1,,z]*survSelS[y,]*WeightAtAgeS[y,])
trueCatchN[z,y] <-projCatchN[z,y]
CatchAssessN[z,y] <-projCatchN[z,y]*exp(CatchErrorN[z,y]-(CatchCVn[y]^2/2))
CatchAssessS[z,y] <-projCatchS[z,y]*exp(CatchErrorS[z,y]-(CatchCVs[y]^2/2))
trueCPUEindN[z,y] <-projExpBn[z,y]
CPUEAssessN[z,y] <-projExpBn[z,y]*exp(CPUEErrorN[z,y]-(IndexCVn[y]^2/2))
CPUEAssessS[z,y] <-projExpBs[z,y]*exp(CPUEErrorS[z,y]-(IndexCVs[y]^2/2))
trueSurvIndN[z,y] <-projSurvN[z,y]+projSurvS[z,y]
SurvAssessN[z,y] <-projSurvN[z,y]*exp(SurvErrorN[z,y]-(IndexCVn[y]^2/2))
SurvAssessS[z,y] <-projSurvS[z,y]*exp(SurvErrorS[z,y]-(IndexCVs[y]^2/2))
trueSpbioN[z,y] <-projSSBn[z,y]
#==UPDATE CATCH LENGTH FREQS====
for(w in 1:ncol(LenAtAgeN))
{
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeN[y,w],sd=GrowthSDn[y])
ProbN<-probtemp/sum(probtemp)
tempCatAtLenN[w,,z]<-projCatchAtAgeN[y,w,z]*ProbN
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeS[y,w],sd=GrowthSDs[y])
ProbS<-probtemp/sum(probtemp)
tempCatAtLenS[w,,z]<-projCatchAtAgeS[y,w,z]*ProbS
}
projCatLenFreqN[y,,z]<-apply(tempCatAtLenN[,,z],2,sum)
projCatLenFreqS[y,,z]<-apply(tempCatAtLenS[,,z],2,sum)
#==============================================================
# calculate the survey proportion at age (length) for assessment
#==============================================================
#==make length frequencies for s==
for(w in 1:ncol(LenAtAgeS))
{
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeN[y,w],sd=GrowthSDn[y])
ProbN<-probtemp/sum(probtemp)
tempSurvAtLenN[w,,z]<-projNn[y-1,w,z]*survSelN[y,w]*ProbN
probtemp<-dnorm(LengthBinsMid,mean=LenAtAgeS[y,w],sd=GrowthSDs[y])
ProbS<-probtemp/sum(probtemp)
tempSurvAtLenS[w,,z]<-projNs[y-1,w,z]*survSelS[y,w]*ProbS
}
projSurvLenFreqN[y,,z]<-apply(tempSurvAtLenN[,,z],2,sum)
projSurvLenFreqS[y,,z]<-apply(tempSurvAtLenS[,,z],2,sum)
# sum(projSurvLenFreqN[y,,z])
# sum(projNn[y,,z])
#==apply F that would take a calculated TAC====================================================
#==update population dynamics
for (i in 2:(MaxAge-1))
{
projNn[y,i,z] <-projNn[y-1,i-1,z]*exp(-ApplyFn*vulnN[y,i-1])*exp(-NatMn[y])
projNs[y,i,z] <-projNs[y-1,i-1,z]*exp(-ApplyFs*vulnS[y,i-1])*exp(-NatMs[y])
}
projNn[y,MaxAge,z] <-(projNn[y-1,(MaxAge-1),z])*exp(-ApplyFn*vulnN[y,MaxAge])*exp(-NatMn[y])+ projNn[y-1,MaxAge,z]*exp(-ApplyFn*vulnN[y,MaxAge])*exp(-NatMn[y])
projNs[y,MaxAge,z] <-(projNs[y-1,(MaxAge-1),z])*exp(-ApplyFs*vulnS[y,MaxAge])*exp(-NatMs[y])+ projNs[y-1,MaxAge,z]*exp(-ApplyFn*vulnS[y,MaxAge])*exp(-NatMs[y])
moveFromN <-projNn[y,,z]*MovementN[y,]
moveFromS <-projNs[y,,z]*MovementS[y,]
projNn[y,,z] <-projNn[y,,z]-moveFromN+moveFromS
projNs[y,,z] <-projNs[y,,z]-moveFromS+moveFromN
EggsN <-sum(projNn[y-1,,z]*matureN[y,]*WeightAtAgeN[y,])
EggsS <-sum(projNs[y-1,,z]*matureS[y,]*WeightAtAgeS[y,])
projNn[y,1,z] <-Recruitment(EggsIN=EggsN,steepnessIN=steepnessN[y],RzeroIN=RzeroN[y],RecErrIN=RecErrN[z,y],recType="BH",NatMin=NatMn[y],
vulnIN=vulnN[y,],matureIN=matureN[y,],weightIN=WeightAtAgeN[y,],LenAtAgeIN=LenAtAgeN[y,],MaxAge=MaxAge,sigmaRin=sigmaRn[y])
projNs[y,1,z] <-Recruitment(EggsIN=EggsS,steepnessIN=steepnessS[y],RzeroIN=RzeroS[y],RecErrIN=RecErrS[z,y],recType="BH",NatMin=NatMs[y],
vulnIN=vulnS[y,],matureIN=matureS[y,],weightIN=WeightAtAgeS[y,],LenAtAgeIN=LenAtAgeS[y,],MaxAge=MaxAge,sigmaRin=sigmaRs[y])
trueRecN[z,y] <-projNn[y,1,z]
trueB35[z,y] <-trueSBPR35[y]*mean(trueRecN[z,start_assessment:SimYear],na.rm=T)
if(AssessmentType==2) {
cat("#True quantities","\n",file=TQfile)
cat("#","\n",file=TQfile,append=TRUE)
cat("#recruitment","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueRecN[q,],"\n",file=TQfile,append=TRUE)
cat("#Catch","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueCatchN[q,],"\n",file=TQfile,append=TRUE)
cat("#fishing mortality","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueFmortN[q,],"\n",file=TQfile,append=TRUE)
cat("#spawning biomass","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueSpbioN[q,],"\n",file=TQfile,append=TRUE)
cat("#survey index","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueSurvIndN[q,],"\n",file=TQfile,append=TRUE)
cat("#cpue index","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueCPUEindN[q,],"\n",file=TQfile,append=TRUE)
#cat("#total allowable catch","\n",file=TQfile,append=TRUE)
#cat(trueTAC[q,],"\n",file=TQfile,append=TRUE)
cat("#","\n",file=TQfile,append=TRUE)
cat("#BMSY begin","\n",file=TQfile,append=TRUE)
cat(trueBMSYbegin,"\n",file=TQfile,append=TRUE)
cat("#FMSY begin","\n",file=TQfile,append=TRUE)
cat(trueFMSYbegin,"\n",file=TQfile,append=TRUE)
cat("#BMSY end","\n",file=TQfile,append=TRUE)
cat(trueBMSYend,"\n",file=TQfile,append=TRUE)
cat("#FMSY end","\n",file=TQfile,append=TRUE)
cat(trueFMSYend,"\n",file=TQfile,append=TRUE)
cat("#B35in","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueB35[q,],"\n",file=TQfile,append=TRUE)
cat("#F35in","\n",file=TQfile,append=TRUE)
cat(trueF35,"\n",file=TQfile,append=TRUE)
cat("#OFL","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueOFL[q,],"\n",file=TQfile,append=TRUE)
TRU<-readLines(TQfile)
if(EstimationPlots==1 & z==Nsim & y == SimYear)
AgeAssPlot(REP,CTL,DAT,TRU,data2,MSEdir,CTLName)
}
if(AssessmentType==0) {
cat("#True quantities","\n",file=TQfile)
cat("#","\n",file=TQfile,append=TRUE)
cat("#recruitment","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueRecN[q,],"\n",file=TQfile,append=TRUE)
cat("#Catch","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueCatchN[q,],"\n",file=TQfile,append=TRUE)
cat("#fishing mortality","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueFmortN[q,],"\n",file=TQfile,append=TRUE)
cat("#spawning biomass","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueSpbioN[q,],"\n",file=TQfile,append=TRUE)
cat("#survey index","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueSurvIndN[q,],"\n",file=TQfile,append=TRUE)
cat("#cpue index","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueCPUEindN[q,],"\n",file=TQfile,append=TRUE)
#cat("#total allowable catch","\n",file=TQfile,append=TRUE)
#cat(trueTACn[q,],"\n",file=TQfile,append=TRUE)
cat("#","\n",file=TQfile,append=TRUE)
cat("#BMSY begin","\n",file=TQfile,append=TRUE)
cat(trueBMSYbegin,"\n",file=TQfile,append=TRUE)
cat("#FMSY begin","\n",file=TQfile,append=TRUE)
cat(trueFMSYbegin,"\n",file=TQfile,append=TRUE)
cat("#BMSY end","\n",file=TQfile,append=TRUE)
cat(trueBMSYend,"\n",file=TQfile,append=TRUE)
cat("#FMSY end","\n",file=TQfile,append=TRUE)
cat(trueFMSYend,"\n",file=TQfile,append=TRUE)
cat("#B35","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueB35[q,],"\n",file=TQfile,append=TRUE)
cat("#F35","\n",file=TQfile,append=TRUE)
cat(trueF35,"\n",file=TQfile,append=TRUE)
cat("#OFL","\n",file=TQfile,append=TRUE)
for(q in 1:Nsim)
cat(trueOFL[q,],"\n",file=TQfile,append=TRUE)
}
if(!silent) cat(paste("Year ",y," of ",SimYear," in simulation ",z," of ",Nsim,"\n",sep=""))
} # end y
} # end z
} # end of GeMS
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.