R/SSS_code_newSRs.r
In shcaba/SSS: Simple Stock Synthesis
Defines functions
SSS
Documented in
SSS
##' This file contails all the SSS code
##' @param filepath location where the SSS will look for model files and run the executable
##' @param file.name vector of file names for the data and control file where the expected input is c("data file", "control file")
##' @param reps number of random draws to perform
##' @param seed.in seed number to fix where random sample are drawn.
##' @param Dep.in vector defining distribution, mean, sd, and bounds for depletion prior. Expected input is c(distribution shape, mean, sd). The distribution options are 2 = 1 - beta, 4 = uniform, 10 = truncated normal.
##' @param M.in vector defining natural mortality distribuition, mean, and . Expected input is c(distrbution shape for females, mean for females, sd for females, distribution shape for males, mean for males, sd for males). The distibution options are 0 = normal, 3 = lognormal, and 4 = uniform.
##' @param SR_type The shape of the stock-recruitment curve. Options are based on SS stock-recruit options. Option 3 = Beverton-holt, 8 = Shepherd 3-parameter, 9 = Ricker 3-parameter
##' @param h.in vector defining the steepness distribution, mean, and sd. Expected input is c(distribution, mean, sd). Distribution options are 1 = truncated beta, 2 = beta, 10 = truncated normal, 30 = truncated lognormal, 4 = uniform.
##' @param FMSY_M.in vector defining the Fmsy/M ratio distribution, mean, and sd. Expected input is c(distribution, mean, sd). Distribution options are; negative value = ?, 2 = truncated beta, 10 = truncated normal, 30 = truncated lognormal, 4 = uniform.
##' @param BMSY_B0.in vector defining the Bmsy/B0 ratio distribution, mean, and sd. Expected input is c(distribution, mean, sd). Distribution options are; negative value = ?, 2 = truncated beta, 10 = truncated normal, 30 = truncated lognormal, 4 = uniform.
##' @param L1.in vector defining the minimum length for a given age. Setting this to 0 in the control file makes it the VBGF t0. A vector of zeros will not draw values for this value. Expected input values are c(female mean, female sd, male mean, male sd)
##' @param Linf.k.cor input defining the correlation between Linf and k when using the multivariate normal distribution
##' @param Linf.in vector defining the maximum length. This is an optional feature. Expected input values are c(prior type (only -1 and 1 (normal) options),female mean, female sd, prior type (only -1 and 1 (normal) options),male mean, male sd)
##' @param k.in vector defining the growth coefficient k. This is an optional feature. Expected input values are c(prior type (only -1 and 1 (normal) options), female mean, female sd, prior type(only -1 and 1 (normal) options),male mean, male sd)
##' @param t0.in vector defining the t0 parameter to calculate the L1 paremter for SS. Expected input values are c(prior type (only -1 and 1 (normal) options), female mean, female sd, prior type (only -1 and 1 (normal) options),male mean, male sd)
##' @param Zfrac.Beta.in is the Zfrac beta for stock recruit function 7 in Stock Synthesis (Survivorship function). Sometimes used with elasmobranchs. Inputs are prior type (- values skips the draws) and prior type inputs.
##' @param R_start vector allowing the user to control the starting R0 value in the control file. Expected value is c( switch option, input value) where the switch optionas are 1= draw from a random draw from a truncated lognormal distribution based on the input value and 0 = start from the input value.
##' @param doR0.loop allows for a profile over initial R0 values in order to find a converged model. It will stop the profile once a converged model is found. Inputs are feature on/off, staring profile value, ending profile, profile step. A 0 for the first value means you will only consider models that start at the given intial R0. Highly recommended to keep this as TRUE, though it can take more computational time.
##' @param sum_age summary age for total biomass calculation used by SS
##' @param ts_yrs start and end years of model year
##' @param ofl_yrs the years for which OFL values should be calculated
##' @param sexes TRUE/FALSE allows for the user to specify whether or not sexeses should have the same values from drawn parameters (e.g., natural mortality, L1, Linf)
##' @param BH_FMSY_comp
##' @param OStype distinguishes operating system being used. "Windows" or "OSX_Linux"
##' @author Jason Cope and Chantel Wetzel
##' @export
##' @seealso \code{\link{Opt_s_prof}}, \code{\link{Change_SSfiles}}, \code{\link{rbeta}}, \code{\link{rbeta_solve}}, \code{\link{rtlnorm}}, \code{\link{Run_SS}}, \code{\link{loadnnames}},
##' @import msm
##' @import EnvStats
##' @import r4ss
SSS<-function(filepath,
file.name,
reps=1000,
seed.in=19,
Dep.in=c(2,0.4,0.1),
M.in=c(3,0.1,0.4,3,0.1,0.4),
SR_type=3,
h.in=c(1,0.6,0.2),
FMSY_M.in=c(-1,0.5,0.1),
BMSY_B0.in=c(-1,0.5,0.1),
Linf.k.cor=-0.9,
Linf.in=c(-1,0,0,-1,0,0),
k.in=c(-1,0,0,-1,0,0),
t0.in=c(-1,0,0,-1,0,0),
Zfrac.Beta.in=c(-99,0.2,0.6,-99,0.5,2),
R_start=c(0,8),
doR0.loop=c(1,4.1,12.1,0.5),
sum_age=0,
ts_yrs=NA,
pop.ltbins=NA,
sexes=F,
BH_FMSY_comp=F,
OStype="Windows")
{
require(msm)
require(EnvStats)
require(r4ss)
require(tmvtnorm)
VBGF<-function(Linf, k, t0, ages){
Linf * (1 - exp(-k * (ages - t0)))
}
VBGF.age<-function(Linf,k,t0,lt){
t0 - (log(1 - (lt / Linf)) / k)
}
set.seed(seed.in)
start.time<-Sys.time()
SSS.output.list<-list()
Spp.quant.out<-list()
OFL.out<-ABC.out<-list()
Input.draws<-as.data.frame(matrix(NA,nrow=reps,ncol=10))
if(SR_type==7){Input.draws<-as.data.frame(matrix(NA,nrow=reps,ncol=11))}
if(SR_type>=8 | h.in[1]==99 | BH_FMSY_comp==T)
{
SR_expo.out<-SRMQ_par<-rep(NA,reps)
names(SR_expo.out)<-"gRicker_Beta"
Input.draws.MQs<-as.data.frame(matrix(NA,nrow=reps,ncol=3))
colnames(Input.draws.MQs)<-c("FMSY/M","BMSY/B0","FMSY")
}
sb.years<-c(ts_yrs[1]:ts_yrs[2])
Quant.out<-as.data.frame(matrix(NA,nrow=reps,ncol=29))
Quant.out.bad<-as.data.frame(matrix(NA,1,ncol=29))
#Input the file names and the summary age into the starter file
starter.new <- SS_readstarter(file.path(filepath, 'starter.ss'),verbose=FALSE)
starter.new$datfile <- file.name[1]
starter.new$ctlfile <- file.name[2]
starter.new$min_age_summary_bio <- sum_age
SS_writestarter(starter.new, filepath, overwrite=TRUE,verbose=FALSE)
#Input the starting and ending years of the model
if(all(!is.na(ts_yrs)))
{
# dat.new<-readLines(paste(filepath,"/",file.name[1],sep=""))
# Styr.line<-strsplit(dat.new[grep("styr",dat.new)], " ")[[1]]
# Styr.line[1]<-ts_yrs[1]
# Endyr.line<-strsplit(dat.new[grep("endyr",dat.new)], " ")[[1]]
# Endyr.line[1]<-ts_yrs[2]
# dat.new[grep("styr",dat.new)]<-paste(Styr.line,collapse=" ")
# dat.new[grep("endyr",dat.new)]<-paste(Endyr.line,collapse=" ")
# write(dat.new,paste(filepath,"/",file.name[1],sep=""))
dat.yrs <- SS_readdat(file.path(filepath, file.name[1]),verbose=FALSE)
dat.yrs$styr<-ts_yrs[1]
dat.yrs$endyr<-ts_yrs[2]
SS_writedat(dat.yrs, file.path(filepath, file.name[1]), overwrite=TRUE,verbose=FALSE)
}
# starter.new<-readLines(paste(filepath,"/starter.ss",sep=""))
# sum_age_line<-strsplit(starter.new[grep("summary biomass",starter.new)], " ")[[1]]
# sum_age_line[1]<-sum_age
# starter.new[grep("summary biomass",starter.new)]<-paste(sum_age_line, collapse=" ")
# write(starter.new,paste(filepath,"/starter.ss",sep=""))
i<-1
ii<-1
xxx<-1
n.bad<-0
while(i<reps+1)
{
print(paste0("Attempt ",ii))
if(ii%%10==0){print(file.name[1])}
if(file.exists(paste(filepath,"/Report.sso",sep=""))){file.remove(paste(filepath,"/Report.sso",sep=""))}
### Draw random values ###
#- values = no draws taken; one value fixed in the model
#0 = normal
#10 = truncated normal
#1 = symmetric beta
#2 = beta (rbeta)
#3 = lognormal
#30 = truncated lognormal
#4 = uniform
#99 = used only for the steepness parameter. Indicates h will come from FMSY/M prior
# dat.new<-readLines(paste(filepath,"/",file.name[1],sep=""))
# ctl.new<-readLines(paste(filepath,"/",file.name[2],sep=""))
dat.new<-SS_readdat(file.path(filepath, file.name[1]),verbose=FALSE)
ctl.new<-SS_readctl(file.path(filepath, file.name[2]),use_datlist = TRUE, datlist=dat.new,verbose=FALSE)
if(length(Dep.in)==3)
{
if(Dep.in[1]==2){Dep.draw<-round(1-rbeta.ab(1,1-Dep.in[2],Dep.in[3],0.05,0.95),2)}
if(Dep.in[1]==3){Dep.draw<-round(rlnorm(1,log(Dep.in[2]),Dep.in[3]),2)}
if(Dep.in[1]==4){Dep.draw<-round(runif(1,Dep.in[2],Dep.in[3]),2)}
if(Dep.in[1]==10){Dep.draw<-round(rtnorm(1,Dep.in[2],Dep.in[3],0.01,1),2)}
Input.draws[i,2]<-Dep.draw
}
if(length(Dep.in)>3)
{
Input.draws[i,2]<-Dep.draw<-Dep.in[i]
}
#Natural mortality
if(FMSY_M.in[1]<0 & BMSY_B0.in[1]<0)
{
if(M.in[1]>=0 & length(M.in)<=6)
{
M.draw<-0
while(M.draw<=0)
{
if(M.in[1]==0){M.draw<-round(rnorm(1,M.in[2],M.in[3]),3)}
if(M.in[1]==3){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
if(M.in[1]==4){M.draw<-round(runif(1,M.in[2],M.in[3]),3)}
Input.draws[i,3]<-M.draw
}
}
else{Input.draws[i,3]<-M.draw<-M.in[i,1]}
if(sexes==T)
{
if(length(M.in)==6)
{
M.draw.M<--1
while(M.draw.M<0)
{
if(M.in[4]==0){M.draw.M<-round(rnorm(1,M.in[5],M.in[6]),3)}
if(M.in[4]==3){M.draw.M<-round(rlnorm(1,log(M.in[5]),M.in[6]),3)}
if(M.in[4]==4){M.draw.M<-round(runif(1,M.in[5],M.in[6]),3)}
Input.draws[i,7]<-M.draw.M
}
}
if(length(M.in)>6){Input.draws[i,7]<-M.draw.M<-M.in[i,2]}
}
}
#Growth parameters
# if(Linf.in[1]>-1&k.in[1]>-1)
# {
Covar_Linf_k<-Linf.k.cor*((Linf.in[3]+0.000001)*(k.in[3]+0.000001)) #caluclate covariance matrix for multivariate sampling
Linf_k_sigma <- matrix(c((Linf.in[3]+0.000001)^2,Covar_Linf_k,Covar_Linf_k,(k.in[3]+0.000001)^2),2,2) #Set-up the variance-covariance matrix for multivariate sampling
Linf_k_samps<-rtmvnorm(1,c(Linf.in[2],k.in[2]),Linf_k_sigma,lower=rep(0,length(c(Linf.in[2],k.in[2]))))
Linf.draw<-Input.draws[i,5]<-Linf_k_samps[1]
k.draw<-Input.draws[i,6]<-Linf_k_samps[2]
#if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
t0.draw<-round(rnorm(1,t0.in[2],t0.in[3]),3)
L1.draw<-Input.draws[i,4]<-Input.draws[i,8]<-0
# }
# if(any(Linf.in[1]<1&k.in[1]<1)&!all(Linf.in[1]<1&k.in[1]<1))
# {
# # if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
# Linf.draw<-Input.draws[i,5]<-round(rnorm(1,Linf.in[2],Linf.in[3]),2)
# k.draw<-Input.draws[i,6]<-round(rnorm(1,k.in[2],k.in[3]),2)
# t0.draw<-round(rnorm(1,t0.in[2],t0.in[3]),3)
# L1.draw<-Input.draws[i,4]<-VBGF(Linf.draw,k.draw,t0.draw,0)
# }
#Male draws
if(sexes==T)
{
# if(sum(L1.in[3:4])>0) {
# L1.draw.M <- round(rnorm(1,L1.in[3], L1.in[4]), 2)
# }
# if(Linf.in[4]>-1&k.in[4]>-1)
# {
if(ctl.new$parameter_offset_approach==1)
{
Covar_Linf_k<-Linf.k.cor*((Linf.in[6]+0.000001)*(k.in[6]+0.000001)) #caluclate covariance matrix for multivariate sampling
Linf_k_sigma <- matrix(c((Linf.in[6]+0.000001)^2,Covar_Linf_k,Covar_Linf_k,(k.in[6]+0.000001)^2),2,2) #Set-up the variance-covariance matrix for multivariate sampling
Linf_k_samps<-rtmvnorm(1,c(Linf.in[5],k.in[5]),Linf_k_sigma,lower=rep(0,length(c(Linf.in[5],k.in[5]))))
Linf.draw.M<-Input.draws[i,9]<-Linf_k_samps[1]
k.draw.M<-Input.draws[i,10]<-Linf_k_samps[2]
#if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
t0.draw.M<-round(rnorm(1,t0.in[5],t0.in[6]),3)
L1.draw.M<-Input.draws[i,8]<-VBGF(Linf.draw.M,k.draw.M,t0.draw.M,t0.draw)
# }
}
if(ctl.new$parameter_offset_approach==2)
{
Linf.draw.M<-Input.draws[i,9]<-Linf.in[5]
k.draw.M<-Input.draws[i,10]<-k.in[5]
#if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
t0.draw.M<-round(rnorm(1,t0.in[5],t0.in[6]),3)
L1.draw.M<-Input.draws[i,8]<-VBGF(Linf.draw.M,k.draw.M,t0.draw.M,t0.draw)
}
#Change to offset approach
# if(any(Linf.in[4]<1&k.in[4]<1)&!all(Linf.in[4]<1&k.in[4]<1))
# {
# if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
# Linf.draw.M<-Input.draws.M[i,3]<-round(rnorm(1,Linf.in[5],Linf.in[6]),2)
# k.draw.M<-Input.draws.M[i,4]<-round(rnorm(1,k.in[5],k.in[6]),2)
# t0.draw>M<-round(rnorm(1,t0.in[5],t0.in[6]),3)
# L1.draw>M<-Input.draws.M[i,2]<-VBGF(Linf.draw.M,k.draw.M,t0.draw.M,0)
# }
}
#Steenpess
if(h.in[1]>=0 & length(h.in)==3 & FMSY_M.in[1]<0 & BMSY_B0.in[1]<0)
{
if(h.in[1]==1){h.draw<-round(rbeta.ab(1,h.in[2],h.in[3],0.21,0.99),3)}
if(h.in[1]==2){h.draw<-round(rbeta(1,h.in[2],h.in[3]),3)}
if(h.in[1]==10){h.draw<-round(rtnorm(1,h.in[2],h.in[3],0.21,0.99),3)}
if(h.in[1]==30){h.draw<-round(rlnormTrunc(1,log(h.in[2]),h.in[3],0.21,0.99),2)}
if(h.in[1]==4){h.draw<-round(runif(1,h.in[2],h.in[3]),2)}
Input.draws[i,1]<-h.draw
}
if(length(h.in)>3 & FMSY_M.in[1]<0 & BMSY_B0.in[1]<0){Input.draws[i,1]<-h.draw<-h.in[i]}
#Get steepness for BH from FMSY
if(h.in[1]==99 | BH_FMSY_comp==T & FMSY_M.in[1]>=0)
{
SRdat.new<-readLines(paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
x<-1
while(x==1)
{
#Draw Ms
if(M.in[1]>=0 & length(M.in)==6)
{
if(M.in[1]==0){M.draw<-round(rnorm(1,M.in[2],M.in[3]),3)}
if(M.in[1]==3){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
if(M.in[1]==4){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
Input.draws[i,3]<-M.draw
}
else{Input.draws[i,3]<-M.draw<-M.in[i,1]}
#Male draws
if(sexes==T)
{
if(length(M.in)==6)
{
if(M.in[4]==0){M.draw.M<-round(rnorm(1,M.in[5],M.in[6]),3)}
if(M.in[4]==3){M.draw.M<-round(rlnorm(1,log(M.in[5]),M.in[6]),3)}
if(M.in[4]==4){M.draw.M<-round(runif(1,M.in[5],M.in[6]),3)}
Input.draws[i,7]<-M.draw.M
}
if(length(M.in)>6){Input.draws[i,7]<-M.draw.M<-M.in[i,2]}
}
if(FMSY_M.in[1]==2){FMSY_M.draw<-round(rbeta.ab(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==10){FMSY_M.draw<-round(rtnorm(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==30){FMSY_M.draw<-round(rlnormTrunc(1,log(FMSY_M.in[2]),FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==4){FMSY_M.draw<-round(runif(1,FMSY_M.in[2],FMSY_M.in[3]),3)}
if(length(FMSY_M.in)>3){FMSY_M.draw<-FMSY_M.in[i]}
Input.draws.MQs[i,1]<-FMSY_M.draw
if(sexes==F)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw),sep=" ")
M_FMSY<-M.draw
}
if(sexes==T)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw.M),sep=" ")
M_FMSY<-mean(c(M.draw,M.draw.M))
}
SRdat.new[17]<-as.character(FMSY_M.draw*M_FMSY)
Input.draws.MQs[i,3]<-FMSY_M.draw*M_FMSY
SRdat.new[21]<-"-0.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
if(file.exists(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))==TRUE)
{
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
SRdat.new[21]<-"-0.1"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
SRdat.new[21]<-"-1.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
}
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
SRdat.new[21]<-"-0.51"
SRdat.new[23]<-"-0.31"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
}
}
if(steep_expo.out[1]>=0.25&steep_expo.out[1]<1&steep_expo.out[2]<=10&SRMQ_par[i]<0.1){x<-2}
}
}
Input.draws[i,1]<-h.draw<-steep_expo.out[1]
SR_expo.out[i]<-steep_expo.out[2]
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep=""))
}
if(FMSY_M.in[1]>=0 & BMSY_B0.in[1]>=0)
{
SRdat.new<-readLines(paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
#SRdat.new.mod<-strsplit(SRdat.new[grep("Amax",SRdat.new)+1][1]," ")[[1]]
#Amax.line<-strsplit(dat.new[grep("Number of ages",dat.new)], " ")[[1]]
#k.line<-strsplit(ctl.new[grep("VonBert_K_Fem_GP_1",ctl.new)]," ")[[1]]
#ltwt.exp.line<-strsplit(ctl.new[grep("Wtlen_2_Fem",ctl.new)]," ")[[1]]
x<-1
while(x==1)
{
#Draw Ms
if(M.in[1]>=0 & length(M.in)==6)
{
if(M.in[1]==0){M.draw<-round(rnorm(1,M.in[2],M.in[3]),3)}
if(M.in[1]==3){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
if(M.in[1]==4){M.draw<-round(runif(1,M.in[2],M.in[3]),3)}
Input.draws[i,3]<-M.draw
}
else{Input.draws[i,3]<-M.draw<-M.in[i,1]}
#Male draws
if(sexes==T)
{
if(length(M.in)==6)
{
if(M.in[4]==0){M.draw.M<-round(rnorm(1,M.in[5],M.in[6]),3)}
if(M.in[4]==3){M.draw.M<-round(rlnorm(1,log(M.in[5]),M.in[6]),3)}
if(M.in[4]==4){M.draw.M<-round(runif(1,M.in[5],M.in[6]),3)}
Input.draws[i,7]<-M.draw.M
}
if(length(M.in)>6){Input.draws[i,7]<-M.draw.M<-M.in[i,2]}
}
if(FMSY_M.in[1]==2){FMSY_M.draw<-round(rbeta.ab(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==10){FMSY_M.draw<-round(rtnorm(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==30){FMSY_M.draw<-round(rlnormTrunc(1,log(FMSY_M.in[2]),FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==4){FMSY_M.draw<-round(runif(1,FMSY_M.in[2],FMSY_M.in[3]),3)}
if(length(FMSY_M.in)>3){FMSY_M.draw<-FMSY_M.in[i]}
if(BMSY_B0.in[1]==2){BMSY_B0.draw<-round(rbeta.ab(1,BMSY_B0.in[2],BMSY_B0.in[3],0.01,0.99),2)*100}
if(BMSY_B0.in[1]==10){BMSY_B0.draw<-round(rtnorm(1,BMSY_B0.in[2],BMSY_B0.in[3],0.01,0.99),2)*100}
if(BMSY_B0.in[1]==30){BMSY_B0.draw<-round(rlnormTrunc(1,log(BMSY_B0.in[2]),BMSY_B0.in[3],0.01,0.99),2)*100}
if(BMSY_B0.in[1]==4){BMSY_B0.draw<-round(runif(1,BMSY_B0.in[2],BMSY_B0.in[3]),2)*100}
if(length(BMSY_B0.in)>3){BMSY_B0.draw<-BMSY_B0.in[i]}
Input.draws.MQs[i,1]<-FMSY_M.draw
Input.draws.MQs[i,2]<-BMSY_B0.draw
#SRdat.new.mod[1]<-as.numeric(Amax.line[1])
if(sexes==F)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw),sep=" ")
M_FMSY<-M.draw
}
if(sexes==T)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw.M),sep=" ")
#M_FMSY<-mean(c(M.draw,M.draw.M))
M_FMSY<-M.draw
}
#print(c(M.draw,M.draw.M,M_FMSY,FMSY_M.draw,BMSY_B0.draw))
#SRdat.new.mod[3]<-Winf
#SRdat.new.mod[4]<-as.numeric(k.line[3])
#SRdat.new.mod[5]<-as.numeric(ltwt.exp.line[3])
#SRdat.new.mod[6]<-SRdat.new.mod[7]<-Amat
#SRdat.new.mod[8]<-SR_type-7 #-7 converts from SS SR_type to MQ SR_type
SRdat.new[17]<-as.character(FMSY_M.draw*M_FMSY)
SRdat.new[19]<-as.character(BMSY_B0.draw)
Input.draws.MQs[i,3]<-FMSY_M.draw*M_FMSY
#SRdat.new[grep("Amax",SRdat.new)+1][1]<-paste(SRdat.new.mod,collapse=" ")
SRdat.new[21]<-"-0.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
if(file.exists(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))==TRUE)
{
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
if(steep_expo.out[1]<0.2|steep_expo.out[1]>=2|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
#print("retry 1")
SRdat.new[21]<-"-0.1"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
#print("retry 2")
SRdat.new[21]<-"-1.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
}
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
#print("retry 3")
SRdat.new[21]<-"-0.51"
SRdat.new[23]<-"-0.31"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
}
}
if(steep_expo.out[1]>=0.2&steep_expo.out[1]<2&steep_expo.out[2]<=10&SRMQ_par[i]<0.1){x<-2}
}
}
Input.draws[i,1]<-h.draw<-steep_expo.out[1]
SR_expo.out[i]<-steep_expo.out[2]
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep=""))
}
#Three parameter S-R model
if(Zfrac.Beta.in[1]>=0)
{
if(Zfrac.Beta.in[1]==10){Zfrac.draw<-round(rtnorm(1,Zfrac.Beta.in[2],Zfrac.Beta.in[3],0,1),2)}
if(Zfrac.Beta.in[1]==30){Zfrac.draw<-round(rlnormTrunc(1,log(Zfrac.Beta.in[2]),Zfrac.Beta.in[3],0,1),2)}
if(Zfrac.Beta.in[1]==4){Zfrac.draw<-round(runif(1,Zfrac.Beta.in[2],Zfrac.Beta.in[3]),2)}
Input.draws[i,1]<-Zfrac.draw
}
if(Zfrac.Beta.in[4]>=0)
{
if(Zfrac.Beta.in[4]==10){Beta.draw<-round(rtnorm(1,Zfrac.Beta.in[5],Zfrac.Beta.in[6],0,10),2)}
if(Zfrac.Beta.in[4]==30){Beta.draw<-round(rlnormTrunc(1,log(Zfrac.Beta.in[5]),Zfrac.Beta.in[6],0,10),2)}
if(Zfrac.Beta.in[4]==4){Beta.draw<-round(runif(1,Zfrac.Beta.in[5],Zfrac.Beta.in[6]),2)}
Input.draws[i,7]<-Beta.draw
}
### Change DAT and CTL inputs ###
#change Depletion
if(Dep.in[1]>=0)
{
dat.new$CPUE[2,4]<-Dep.draw
SS_writedat(dat.new, paste0(filepath,"/",file.name[1]), overwrite=TRUE,verbose=FALSE)
}
#change M
Sys.sleep(1)
if(M.in[1]>=0){ctl.new$MG_parms[1,3]<-M.draw}
#change growth parameters
ctl.new$Growth_Age_for_L1<-t0.draw
ctl.new$MG_parms[2,3:4]<-ctl.new$MG_parms[14,3]<-0
ctl.new$MG_parms[3,3:4]<-Linf.draw
ctl.new$MG_parms[4,3:4]<-k.draw
#change male pararmeters
if(sexes==T)
{
#change M
if(M.in[4]>=0){ctl.new$MG_parms[13,3]<-M.draw.M}
#change growth parameters
ctl.new$MG_parms[14,3]<-L1.draw.M
ctl.new$MG_parms[15,3]<-Linf.draw.M
ctl.new$MG_parms[16,3]<-k.draw.M
}
ctl.new$SR_function<-SR_type
# SRtype.line<-strsplit(ctl.new[grep("SR_function",ctl.new)], " ")[[1]]
# SRtype.line[1]<-SR_type
# ctl.new[grep("SR_function",ctl.new)]<-paste(SRtype.line,collapse=" ")
#change R0
if(R_start[1]==0)
{
R0.draw<-R_start[2]
ctl.new$SR_parms[1,3:4]<-R0.draw
# R0.line<-strsplit(ctl.new[grep("R0",ctl.new)], " ")[[1]]
# R0.line[c(3,4)]<-R0.draw
# ctl.new[grep("R0",ctl.new)]<-paste(R0.line,collapse=" ")
}
if(R_start[1]==1)
{
R0.draw<-round(rlnormTrunc(1,log(R_start[2]),0.5,3,15),2)
ctl.new$SR_parms[1,3:4]<-R0.draw
# R0.line<-strsplit(ctl.new[grep("R0",ctl.new)], " ")[[1]]
# R0.line[c(3,4)]<-R0.draw
# ctl.new[grep("R0",ctl.new)]<-paste(R0.line,collapse=" ")
}
#change h
if(SR_type==3)
{
ctl.new$SR_parms[2,3:4]<-h.draw
# h.line<-strsplit(ctl.new[grep("SR_BH_steep",ctl.new)], " ")[[1]]
# h.line[c(3,4)]<-h.draw
# ctl.new[grep("SR_BH_steep",ctl.new)]<-paste(h.line,collapse=" ")
}
if(SR_type==9)
{
ctl.new$SR_parms[2,3:4]<-h.draw
ctl.new$SR_parms[3,3:4]<-SR_expo.out[i]
h.line<-strsplit(ctl.new[grep("SR_RkrPower_steep",ctl.new)], " ")[[1]]
# h.line[c(3,4)]<-h.draw
# ctl.new[grep("SR_RkrPower_steep",ctl.new)]<-paste(h.line,collapse=" ")
# SRexpo.line<-strsplit(ctl.new[grep("SR_RkrPower_gamma",ctl.new)], " ")[[1]]
# SRexpo.line[c(3,4)]<-SR_expo.out[i]
# ctl.new[grep("SR_RkrPower_gamma",ctl.new)]<-paste(SRexpo.line,collapse=" ")
}
#change Sfrac and Beta
if(SR_type==7 & Zfrac.Beta.in[1]>=0)
{
ctl.new$SR_parms[2,3:4]<-Zfrac.draw
# Zfrac.line<-strsplit(ctl.new[grep("Zfrac",ctl.new)], " ")[[1]]
# Zfrac.line[c(3,4)]<-Zfrac.draw
# ctl.new[grep("Zfrac",ctl.new)]<-paste(Zfrac.line,collapse=" ")
}
if(SR_type==7 & Zfrac.Beta.in[4]>=0)
{
ctl.new$SR_parms[3,3:4]<-Beta.draw
# Beta.line<-strsplit(ctl.new[grep("Beta",ctl.new)], " ")[[1]]
# Beta.line[c(3,4)]<-Beta.draw
# ctl.new[grep("Beta",ctl.new)]<-paste(Beta.line,collapse=" ")
}
#write(ctl.new,paste(filepath,"/",file.name[2],sep=""))
SS_writectl(ctl.new, paste0(filepath,"/",file.name[2]), overwrite=TRUE,verbose=FALSE)
#Run model
RUN.SS(paste(filepath,"/",sep=""),ss.cmd=" -nohess -nox > out.txt 2>&1",OS.in=OStype,show_in_console = FALSE)
#if(OStype=="Windows"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
#if(OStype=="OSX_Linux"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="./ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
#Evaluate convergence and record values
rep.new<-readLines(paste(filepath,"/Report.sso",sep=""))
Sys.sleep(0.5)
#The R0 loop. Makes sure R0 is changing from input and that the depletion match is happening. Ensures good models are kept.
#Profiles over R0 when intial run fails.
if(doR0.loop[1]>0)
{
xx<-1
R0.explore<-seq(doR0.loop[2],doR0.loop[3],doR0.loop[4])
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==TRUE|as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])==0){Dep.out.testR0<-10}
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==FALSE&as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])>0){Dep.out.testR0<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
# while(abs(as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][3])-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][8]))==0)
while(abs(Dep.out.testR0-Dep.draw)>0.01)
{
print(paste0("IN THE LOOP: R0= ",R0.explore[xx]))
ctl.new$SR_parms[1,3:4]<-R0.explore[xx]
SS_writectl(ctl.new, paste0(filepath,"/",file.name[2]), overwrite=TRUE,verbose=FALSE)
# ctl.new<-readLines(paste(filepath,"/",file.name[2],sep=""))
# R0.line<-strsplit(ctl.new[grep("R0",ctl.new)], " ")[[1]]
# R0.line[c(3,4)]<-R0.explore[xx]
# ctl.new[grep("R0",ctl.new)]<-paste(R0.line,collapse=" ")
# write(ctl.new,paste(filepath,"/",file.name[2],sep=""))
RUN.SS(paste(filepath,"/",sep=""),ss.cmd=" -nohess -nox > out.txt 2>&1",OS.in=OStype)
# if(OStype=="Windows"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
# if(OStype=="OSX_Linux"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="./ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
rep.new<-readLines(paste(filepath,"/Report.sso",sep=""))
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==TRUE|as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])==0){Dep.out.testR0<-10}
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==FALSE&as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])>0){Dep.out.testR0<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
print(paste0("Is depletion difference < 0.01? ",(abs(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])-Dep.draw)<0.01)))
xx<-xx+1
if(xx==length(R0.explore)){break}
}
}
SSS.output.list[[i]]<-SS_output(paste(filepath,"/",sep=""), covar=FALSE,verbose=FALSE,printstats=FALSE)
#Begin extracting info from run
#forecast.file<-readLines(paste(filepath,"/Forecast-report.sso",sep=""))
# for(iii in 1:length(ts_yrs[1]:ts_yrs[2]))
# {
# SB.out[i,iii]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",sb.years[iii],sep=""),rep.new)], " ")[[1]][3])
# TB.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][5])
# SumAge.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][6])
# SPR.out[i,iii]<-1-as.numeric(strsplit(rep.new[grep(paste("SPRratio_",sb.years[iii]+1,sep=""),rep.new)], " ")[[1]][3])
# # B_BMSY.out[i,iii]<-as.numeric(strsplit(rep.new[grep("Yr B/Bmsy F/Fmsy",rep.new)+iii], " ")[[1]][2])
# }
if(i==1)
{
SB.out<-TB.out<-SumAge.out<-SPR.out<-as.data.frame(matrix(NA,nrow=reps,ncol=length(SSS.output.list[[i]]$timeseries$SpawnBio)))
colnames(SB.out)<-colnames(TB.out)<-colnames(SumAge.out)<-colnames(SPR.out)<-sb.years
}
SB.out[i,]<-SSS.output.list[[i]]$timeseries$SpawnBio
TB.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_all
SumAge.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_smry
#SPR.out[i,]<-1-as.numeric(strsplit(rep.new[grep(paste("SPRratio_",sb.years[iii]+1,sep=""),rep.new)], " ")[[1]][3])
nforecasts<- SSS.output.list[[i]]$nforecastyears
OFL.out[[i]]<-SSS.output.list[[i]]$derived_quants[paste0("OFLCatch_", c((ts_yrs[2]+1):(ts_yrs[2]+nforecasts))), "Value"]
ABC.out[[i]]<-SSS.output.list[[i]]$derived_quants[paste0("ForeCatch_", c((ts_yrs[2]+1):(ts_yrs[2]+nforecasts))), "Value"]
Dep.series.out<-SB.out/SB.out[,1]
#Quant.out[i,1]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,2]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,3]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,4]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Fem_GP_1",rep.new)], " ")[[1]][3])
Quant.out[i,1]<-SSS.output.list[[i]]$parameters[1,3]
Quant.out[i,2]<-SSS.output.list[[i]]$parameters[2,3]
Quant.out[i,3]<-SSS.output.list[[i]]$parameters[3,3]
Quant.out[i,4]<-SSS.output.list[[i]]$parameters[4,3]
if(sexes==T)
{
# Quant.out[i,5]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,6]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,7]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,8]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Mal_GP_1",rep.new)], " ")[[1]][3])
Quant.out[i,5]<-SSS.output.list[[i]]$parameters[13,3]
Quant.out[i,6]<-SSS.output.list[[i]]$parameters[14,3]
Quant.out[i,7]<-SSS.output.list[[i]]$parameters[15,3]
Quant.out[i,8]<-SSS.output.list[[i]]$parameters[16,3]
}
# if(SR_type==3){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_BH_steep",rep.new)], " ")[[1]][3])}
if(SR_type==3){Quant.out[i,9]<-SSS.output.list[[i]]$parameters[24,3]}
if(SR_type==7){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)-1], " ")[[1]][2])}
if(SR_type==7){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)], " ")[[1]][2])}
if(SR_type==8){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_steep",rep.new)], " ")[[1]][2])}
if(SR_type==8){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_gamma",rep.new)], " ")[[1]][2])}
#Quant.out[i,10]<-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][3])
Quant.out[i,10]<-SSS.output.list[[i]]$parameters[23,3]
#Quant.out[i,11]<-as.numeric(strsplit(rep.new[grep("SSB_Virgin",rep.new)], " ")[[1]][3])
Quant.out[i,11]<-SSS.output.list[[i]]$SBzero
Quant.out[i,12]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("SSB_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,13]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("Bratio_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,14]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("SSB_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,15]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("B_MSY/SSB_unfished",SSS.output.list[[i]]$derived_quants[,1]),][2])
if(!is.na(as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2]))){Quant.out[i,16]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])}
Quant.out[i,17]<-SSS.output.list[[i]]$likelihoods_used[1,1]
Quant.out[i,18]<-SSS.output.list[[i]]$likelihoods_used[4,1]
Quant.out[i,19]<-SSS.output.list[[i]]$maximum_gradient_component
if(Dep.in[1]>=0){Quant.out[i,20]<-as.numeric(dat.new$CPUE[2,4])}
if(Dep.in[1]>=0){Quant.out[i,21]<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
Quant.out[i,22]<-SSS.output.list[[i]]$parameters[23,8]
Quant.out[i,23]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("F_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])/as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,24]<-SSS.output.list[[i]]$likelihoods_used[11,1]
warnings.out<-readLines(paste(filepath,"/warning.sso",sep=""))
num.warning<-as.numeric(strsplit(warnings.out[grep("N warnings",warnings.out)[1]], " ")[[1]][4])
#Quant.out[i,12]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])
#Quant.out[i,13]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("SSB_Virgin",rep.new)], " ")[[1]][3])
#Quant.out[i,14]<-as.numeric(strsplit(rep.new[grep("SSB_MSY",rep.new)], " ")[[1]][2])/as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
#Quant.out[i,15]<-as.numeric(strsplit(rep.new[grep("Bmsy/Bzero",rep.new)], " ")[[1]][1])
#Quant.out[i,17]<-as.numeric(strsplit(forecast.file[grep("calculate_FMSY",forecast.file)+13],split="[[:blank:]]+")[[1]][2])/as.numeric(strsplit(forecast.file[grep("BIO_Smry_unfished",forecast.file)],split="[[:blank:]]+")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2]))){Quant.out[i,16]<-as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])}
# Quant.out[i,17]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)], " ")[[1]][2])
# Quant.out[i,18]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)+2], " ")[[1]][2])
# Quant.out[i,19]<-as.numeric(strsplit(rep.new[grep("Convergence",rep.new)], " ")[[1]][2])
# Quant.out[i,23]<-as.numeric(strsplit(rep.new[grep(paste("F_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,26]<-as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,27]<-as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
#Quant.out[i,24]<-as.numeric(strsplit(rep.new[grep("Crash_Pen",rep.new)], " ")[[1]][2])
if(Dep.in[1]>=0)
{
if(FMSY_M.in[1]<0 & BMSY_B0.in[1]<0&abs(Quant.out[i,20]-Quant.out[i,21])<0.01&Quant.out[i,24]<0.0001){i<-i+1}
if(FMSY_M.in[1]>=0 & BMSY_B0.in[1]>=0 | h.in[1]==99 | BH_FMSY_comp==T)
{
if(h.in[1]==99)
{
if(abs(Quant.out[i,20]-Quant.out[i,21])<0.01&Quant.out[i,24]<0.0001&(abs(Quant.out[i,14]-Input.draws.MQs[i,3])/Input.draws.MQs[i,3])<0.6){i<-i+1}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
if(length(M.in)>6|length(Dep.in)>3|length(h.in)>3|length(FMSY_M.in)>3|length(BMSY_B0.in)>3)
{
Quant.out[i,]<-NA
i<-i+1
}
}
}
else
{
if(abs(Quant.out[i,20]-Quant.out[i,21])<0.01&Quant.out[i,24]<0.0001&(abs(Quant.out[i,14]-(Input.draws.MQs[i,2])/100)/(Input.draws.MQs[i,2]/100))<0.5){i<-i+1}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
if(length(M.in)>6|length(Dep.in)>3|length(h.in)>3|length(FMSY_M.in)>3|length(BMSY_B0.in)>3)
{
Quant.out[i,]<-NA
i<-i+1
}
}
}
}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
if(length(M.in)>6|length(Dep.in)>3|length(h.in)>3|length(FMSY_M.in)>3|length(BMSY_B0.in)>3)
{
Quant.out[i,]<-NA
i<-i+1
}
}
# }
}
if(Dep.in[1]<0)
{
if(Quant.out[i,24]<0.0001){i<-i+1}
else
{
#Run using .par file if above conditions not met
print("*** RUNNING WITH .PAR ***")
starter.new<-readLines(paste(filepath,"/starter.ss",sep=""))
par_line<-strsplit(starter.new[grep("ss.par",starter.new)], " ")[[1]]
par_line[1]<-1
starter.new[grep("ss.par",starter.new)]<-paste(par_line, collapse=" ")
write(starter.new,paste(filepath,"/starter.ss",sep=""))
RUN.SS(paste(filepath,"/",sep=""),ss.cmd=" -nohess -nox > out.txt 2>&1",OS.in=OStype)
# if(OStype=="Windows"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
# if(OStype=="OSX_Linux"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="./ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
#RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")
rep.new<-readLines(paste(filepath,"/Report.sso",sep=""))
SSS.output.list[[i]]<-SS_output(paste(filepath,"/",sep=""), covar=FALSE,verbose=FALSE,printstats=FALSE)
#forecast.file<-readLines(paste(filepath,"/Forecast-report.sso",sep=""))
# for(iii in 1:length(ts_yrs[1]:ts_yrs[2]))
# {
# SB.out[i,iii]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",sb.years[iii],sep=""),rep.new)], " ")[[1]][3])
# TB.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][5])
# SumAge.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][6])
# SPR.out[i,iii]<-1-as.numeric(strsplit(rep.new[grep(paste("SPRratio_",sb.years[iii]+1,sep=""),rep.new)], " ")[[1]][3])
# # B_BMSY.out[i,iii]<-as.numeric(strsplit(rep.new[grep("Yr B/Bmsy F/Fmsy",rep.new)+iii], " ")[[1]][2])
# }
SB.out[i,]<-SSS.output.list[[i]]$timeseries$SpawnBio
TB.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_all
SumAge.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_smry
Dep.series.out<-SB.out/SB.out[,1]
Quant.out[i,1]<-SSS.output.list[[i]]$parameters[1,3]
Quant.out[i,2]<-SSS.output.list[[i]]$parameters[2,3]
Quant.out[i,3]<-SSS.output.list[[i]]$parameters[3,3]
Quant.out[i,4]<-SSS.output.list[[i]]$parameters[4,3]
# Quant.out[i,1]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,2]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,3]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,4]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Fem_GP_1",rep.new)], " ")[[1]][3])
if(sexes==T)
{
Quant.out[i,5]<-SSS.output.list[[i]]$parameters[13,3]
Quant.out[i,6]<-SSS.output.list[[i]]$parameters[14,3]
Quant.out[i,7]<-SSS.output.list[[i]]$parameters[15,3]
Quant.out[i,8]<-SSS.output.list[[i]]$parameters[16,3]
# Quant.out[i,5]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,6]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,7]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,8]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Mal_GP_1",rep.new)], " ")[[1]][3])
}
#if(SR_type==3){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_BH_steep",rep.new)], " ")[[1]][3])}
if(SR_type==3){Quant.out[i,9]<-SSS.output.list[[i]]$parameters[24,3]}
if(SR_type==7){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)-1], " ")[[1]][3])}
if(SR_type==7){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)], " ")[[1]][3])}
if(SR_type==8){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_steep",rep.new)], " ")[[1]][3])}
if(SR_type==8){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_gamma",rep.new)], " ")[[1]][3])}
Quant.out[i,10]<-SSS.output.list[[i]]$parameters[23,3]
Quant.out[i,11]<-SSS.output.list[[i]]$SBzero
Quant.out[i,12]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("SSB_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,13]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("Bratio_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,14]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("SSB_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,15]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("B_MSY/SSB_unfished",SSS.output.list[[i]]$derived_quants[,1]),][2])
if(!is.na(as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2]))){Quant.out[i,16]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])}
Quant.out[i,17]<-SSS.output.list[[i]]$likelihoods_used[1,1]
Quant.out[i,18]<-SSS.output.list[[i]]$likelihoods_used[4,1]
Quant.out[i,19]<-SSS.output.list[[i]]$maximum_gradient_component
if(Dep.in[1]>=0){Quant.out[i,20]<-as.numeric(dat.new$CPUE[2,4])}
if(Dep.in[1]>=0){Quant.out[i,21]<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
Quant.out[i,22]<-SSS.output.list[[i]]$parameters[23,8]
Quant.out[i,23]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("F_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])/as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,24]<-SSS.output.list[[i]]$likelihoods_used[11,1]
# Quant.out[i,10]<-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][3])
# Quant.out[i,11]<-as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
# Quant.out[i,12]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])
# Quant.out[i,13]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
# #Quant.out[i,14]<-as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+1,sep=""),rep.new)], " ")[[1]][2])
# #Quant.out[i,15]<-as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+1,sep=""),rep.new)], " ")[[1]][2])
# Quant.out[i,14]<-as.numeric(strsplit(rep.new[grep("SSB_MSY",rep.new)], " ")[[1]][2])/as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
# Quant.out[i,15]<-as.numeric(strsplit(rep.new[grep("Bmsy/Bzero",rep.new)], " ")[[1]][1])
# #Quant.out[i,17]<-as.numeric(strsplit(forecast.file[grep("Maximum_Sustainable_Yield;_where_Yield_is_Dead_Catch",forecast.file)+14],split="[[:blank:]]+")[[1]][2])/as.numeric(strsplit(forecast.file[grep("BIO_Smry_unfished",forecast.file)],split="[[:blank:]]+")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2]))){Quant.out[i,16]<-as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])}
# Quant.out[i,17]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)], " ")[[1]][2])
# Quant.out[i,18]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)+2], " ")[[1]][2])
# Quant.out[i,19]<-as.numeric(strsplit(rep.new[grep("Convergence",rep.new)], " ")[[1]][2])
# Quant.out[i,20]<-as.numeric(dat.new$CPUE[2,4])
# Quant.out[i,21]<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][3])
# Quant.out[i,22]<-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][8])
# Quant.out[i,23]<-as.numeric(strsplit(rep.new[grep(paste("F_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,26]<-as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,27]<-as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
# Quant.out[i,24]<-as.numeric(strsplit(rep.new[grep("Crash_Pen",rep.new)], " ")[[1]][2])
starter.new<-readLines(paste(filepath,"/starter.ss",sep=""))
par_line<-strsplit(starter.new[grep("ss.par",starter.new)], " ")[[1]]
par_line[1]<-0
starter.new[grep("ss.par",starter.new)]<-paste(par_line, collapse=" ")
write(starter.new,paste(filepath,"/starter.ss",sep=""))
if(Quant.out[i,24]<0.0001){i<-i+1}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
}
}
}
print(paste0("Kept ",i-1))
ii<-ii+1
if(file.exists(paste(filepath,"/Report.sso",sep=""))){file.remove(paste(filepath,"/Forecast-report.sso",sep=""))}
if(file.exists(paste(filepath,"/Report.sso",sep=""))){file.remove(paste(filepath,"/Report.sso",sep=""))}
colnames(Quant.out)<-colnames(Quant.out.bad)<-c("M_f","L1_f","Linf_f","k_f","M_m","L1_m","Linf_m","k_m","h","R0","SB0",paste("SSB_",ts_yrs[2],sep=""),"Term_Yr_Dep","SBMSY","SBMSY/SB0","FMSY","-lnL","LL_survey","Gradient","Dep.Obs","Dep.Exp","R0_init",paste("F_",ts_yrs[2],"/FMSY",sep=""),"Crash_penalty","Rick_gamma")
if(ncol(Input.draws)==10){colnames(Input.draws)<-c("h","Dep","M_f","L1_f","Linf_f","k_f","M_m","L1_m","Linf_m","k_m")}
if(ncol(Input.draws)==11){colnames(Input.draws)<-c("Sfrac","Dep","M_f","L1_f","Linf_f","k_f","Beta","M_m","L1_m","Linf_m","k_m")}
if(SR_type>=8 | h.in[1]==99 | BH_FMSY_comp==T)
{
Input.draws<-cbind(Input.draws,Input.draws.MQs,SR_expo.out,SRMQ_par)
ltcolnames<-length(colnames(Input.draws))
colnames(Input.draws)[c(ltcolnames-1,ltcolnames)]<-c("Beta","Obj_fxn")
}
end.time<-Sys.time()
OFL.out.save<-do.call("rbind",OFL.out)
ABC.out.save<-do.call("rbind",ABC.out)
colnames(OFL.out.save)<-colnames(ABC.out.save)<-c((ts_yrs[2]+1):(ts_yrs[2]+nforecasts))
Spp.quant.out<-list(Input.draws,Quant.out,SB.out,Dep.series.out,TB.out,SumAge.out,OFL.out.save,ABC.out.save,Quant.out.bad,ii-1,(as.numeric(end.time)-as.numeric(start.time))/60)
names(Spp.quant.out)<-c("Priors","Posteriors","SB_series","Rel_Stock_status_series","Total_Biomass","Summary_Biomass","OFL","ABC","Rejected_draws","Total draws","Runtime_minutes")
SSS.out<-Spp.quant.out
save(SSS.output.list,file=paste(filepath,"/SSSoutput.DMP",sep=""))
save(SSS.out,file=paste(filepath,"/SSS_out.DMP",sep=""))
}
return(Spp.quant.out)
}
shcaba/SSS documentation built on Jan. 26, 2024, 5:10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions SSS
Documented in SSS
##' This file contails all the SSS code
##' @param filepath location where the SSS will look for model files and run the executable
##' @param file.name vector of file names for the data and control file where the expected input is c("data file", "control file")
##' @param reps number of random draws to perform
##' @param seed.in seed number to fix where random sample are drawn.
##' @param Dep.in vector defining distribution, mean, sd, and bounds for depletion prior. Expected input is c(distribution shape, mean, sd). The distribution options are 2 = 1 - beta, 4 = uniform, 10 = truncated normal.
##' @param M.in vector defining natural mortality distribuition, mean, and . Expected input is c(distrbution shape for females, mean for females, sd for females, distribution shape for males, mean for males, sd for males). The distibution options are 0 = normal, 3 = lognormal, and 4 = uniform.
##' @param SR_type The shape of the stock-recruitment curve. Options are based on SS stock-recruit options. Option 3 = Beverton-holt, 8 = Shepherd 3-parameter, 9 = Ricker 3-parameter
##' @param h.in vector defining the steepness distribution, mean, and sd. Expected input is c(distribution, mean, sd). Distribution options are 1 = truncated beta, 2 = beta, 10 = truncated normal, 30 = truncated lognormal, 4 = uniform.
##' @param FMSY_M.in vector defining the Fmsy/M ratio distribution, mean, and sd. Expected input is c(distribution, mean, sd). Distribution options are; negative value = ?, 2 = truncated beta, 10 = truncated normal, 30 = truncated lognormal, 4 = uniform.
##' @param BMSY_B0.in vector defining the Bmsy/B0 ratio distribution, mean, and sd. Expected input is c(distribution, mean, sd). Distribution options are; negative value = ?, 2 = truncated beta, 10 = truncated normal, 30 = truncated lognormal, 4 = uniform.
##' @param L1.in vector defining the minimum length for a given age. Setting this to 0 in the control file makes it the VBGF t0. A vector of zeros will not draw values for this value. Expected input values are c(female mean, female sd, male mean, male sd)
##' @param Linf.k.cor input defining the correlation between Linf and k when using the multivariate normal distribution
##' @param Linf.in vector defining the maximum length. This is an optional feature. Expected input values are c(prior type (only -1 and 1 (normal) options),female mean, female sd, prior type (only -1 and 1 (normal) options),male mean, male sd)
##' @param k.in vector defining the growth coefficient k. This is an optional feature. Expected input values are c(prior type (only -1 and 1 (normal) options), female mean, female sd, prior type(only -1 and 1 (normal) options),male mean, male sd)
##' @param t0.in vector defining the t0 parameter to calculate the L1 paremter for SS. Expected input values are c(prior type (only -1 and 1 (normal) options), female mean, female sd, prior type (only -1 and 1 (normal) options),male mean, male sd)
##' @param Zfrac.Beta.in is the Zfrac beta for stock recruit function 7 in Stock Synthesis (Survivorship function). Sometimes used with elasmobranchs. Inputs are prior type (- values skips the draws) and prior type inputs.
##' @param R_start vector allowing the user to control the starting R0 value in the control file. Expected value is c( switch option, input value) where the switch optionas are 1= draw from a random draw from a truncated lognormal distribution based on the input value and 0 = start from the input value.
##' @param doR0.loop allows for a profile over initial R0 values in order to find a converged model. It will stop the profile once a converged model is found. Inputs are feature on/off, staring profile value, ending profile, profile step. A 0 for the first value means you will only consider models that start at the given intial R0. Highly recommended to keep this as TRUE, though it can take more computational time.
##' @param sum_age summary age for total biomass calculation used by SS
##' @param ts_yrs start and end years of model year
##' @param ofl_yrs the years for which OFL values should be calculated
##' @param sexes TRUE/FALSE allows for the user to specify whether or not sexeses should have the same values from drawn parameters (e.g., natural mortality, L1, Linf)
##' @param BH_FMSY_comp
##' @param OStype distinguishes operating system being used. "Windows" or "OSX_Linux"
##' @author Jason Cope and Chantel Wetzel
##' @export
##' @seealso \code{\link{Opt_s_prof}}, \code{\link{Change_SSfiles}}, \code{\link{rbeta}}, \code{\link{rbeta_solve}}, \code{\link{rtlnorm}}, \code{\link{Run_SS}}, \code{\link{loadnnames}},
##' @import msm
##' @import EnvStats
##' @import r4ss
SSS<-function(filepath,
file.name,
reps=1000,
seed.in=19,
Dep.in=c(2,0.4,0.1),
M.in=c(3,0.1,0.4,3,0.1,0.4),
SR_type=3,
h.in=c(1,0.6,0.2),
FMSY_M.in=c(-1,0.5,0.1),
BMSY_B0.in=c(-1,0.5,0.1),
Linf.k.cor=-0.9,
Linf.in=c(-1,0,0,-1,0,0),
k.in=c(-1,0,0,-1,0,0),
t0.in=c(-1,0,0,-1,0,0),
Zfrac.Beta.in=c(-99,0.2,0.6,-99,0.5,2),
R_start=c(0,8),
doR0.loop=c(1,4.1,12.1,0.5),
sum_age=0,
ts_yrs=NA,
pop.ltbins=NA,
sexes=F,
BH_FMSY_comp=F,
OStype="Windows")
{
require(msm)
require(EnvStats)
require(r4ss)
require(tmvtnorm)
VBGF<-function(Linf, k, t0, ages){
Linf * (1 - exp(-k * (ages - t0)))
}
VBGF.age<-function(Linf,k,t0,lt){
t0 - (log(1 - (lt / Linf)) / k)
}
set.seed(seed.in)
start.time<-Sys.time()
SSS.output.list<-list()
Spp.quant.out<-list()
OFL.out<-ABC.out<-list()
Input.draws<-as.data.frame(matrix(NA,nrow=reps,ncol=10))
if(SR_type==7){Input.draws<-as.data.frame(matrix(NA,nrow=reps,ncol=11))}
if(SR_type>=8 | h.in[1]==99 | BH_FMSY_comp==T)
{
SR_expo.out<-SRMQ_par<-rep(NA,reps)
names(SR_expo.out)<-"gRicker_Beta"
Input.draws.MQs<-as.data.frame(matrix(NA,nrow=reps,ncol=3))
colnames(Input.draws.MQs)<-c("FMSY/M","BMSY/B0","FMSY")
}
sb.years<-c(ts_yrs[1]:ts_yrs[2])
Quant.out<-as.data.frame(matrix(NA,nrow=reps,ncol=29))
Quant.out.bad<-as.data.frame(matrix(NA,1,ncol=29))
#Input the file names and the summary age into the starter file
starter.new <- SS_readstarter(file.path(filepath, 'starter.ss'),verbose=FALSE)
starter.new$datfile <- file.name[1]
starter.new$ctlfile <- file.name[2]
starter.new$min_age_summary_bio <- sum_age
SS_writestarter(starter.new, filepath, overwrite=TRUE,verbose=FALSE)
#Input the starting and ending years of the model
if(all(!is.na(ts_yrs)))
{
# dat.new<-readLines(paste(filepath,"/",file.name[1],sep=""))
# Styr.line<-strsplit(dat.new[grep("styr",dat.new)], " ")[[1]]
# Styr.line[1]<-ts_yrs[1]
# Endyr.line<-strsplit(dat.new[grep("endyr",dat.new)], " ")[[1]]
# Endyr.line[1]<-ts_yrs[2]
# dat.new[grep("styr",dat.new)]<-paste(Styr.line,collapse=" ")
# dat.new[grep("endyr",dat.new)]<-paste(Endyr.line,collapse=" ")
# write(dat.new,paste(filepath,"/",file.name[1],sep=""))
dat.yrs <- SS_readdat(file.path(filepath, file.name[1]),verbose=FALSE)
dat.yrs$styr<-ts_yrs[1]
dat.yrs$endyr<-ts_yrs[2]
SS_writedat(dat.yrs, file.path(filepath, file.name[1]), overwrite=TRUE,verbose=FALSE)
}
# starter.new<-readLines(paste(filepath,"/starter.ss",sep=""))
# sum_age_line<-strsplit(starter.new[grep("summary biomass",starter.new)], " ")[[1]]
# sum_age_line[1]<-sum_age
# starter.new[grep("summary biomass",starter.new)]<-paste(sum_age_line, collapse=" ")
# write(starter.new,paste(filepath,"/starter.ss",sep=""))
i<-1
ii<-1
xxx<-1
n.bad<-0
while(i<reps+1)
{
print(paste0("Attempt ",ii))
if(ii%%10==0){print(file.name[1])}
if(file.exists(paste(filepath,"/Report.sso",sep=""))){file.remove(paste(filepath,"/Report.sso",sep=""))}
### Draw random values ###
#- values = no draws taken; one value fixed in the model
#0 = normal
#10 = truncated normal
#1 = symmetric beta
#2 = beta (rbeta)
#3 = lognormal
#30 = truncated lognormal
#4 = uniform
#99 = used only for the steepness parameter. Indicates h will come from FMSY/M prior
# dat.new<-readLines(paste(filepath,"/",file.name[1],sep=""))
# ctl.new<-readLines(paste(filepath,"/",file.name[2],sep=""))
dat.new<-SS_readdat(file.path(filepath, file.name[1]),verbose=FALSE)
ctl.new<-SS_readctl(file.path(filepath, file.name[2]),use_datlist = TRUE, datlist=dat.new,verbose=FALSE)
if(length(Dep.in)==3)
{
if(Dep.in[1]==2){Dep.draw<-round(1-rbeta.ab(1,1-Dep.in[2],Dep.in[3],0.05,0.95),2)}
if(Dep.in[1]==3){Dep.draw<-round(rlnorm(1,log(Dep.in[2]),Dep.in[3]),2)}
if(Dep.in[1]==4){Dep.draw<-round(runif(1,Dep.in[2],Dep.in[3]),2)}
if(Dep.in[1]==10){Dep.draw<-round(rtnorm(1,Dep.in[2],Dep.in[3],0.01,1),2)}
Input.draws[i,2]<-Dep.draw
}
if(length(Dep.in)>3)
{
Input.draws[i,2]<-Dep.draw<-Dep.in[i]
}
#Natural mortality
if(FMSY_M.in[1]<0 & BMSY_B0.in[1]<0)
{
if(M.in[1]>=0 & length(M.in)<=6)
{
M.draw<-0
while(M.draw<=0)
{
if(M.in[1]==0){M.draw<-round(rnorm(1,M.in[2],M.in[3]),3)}
if(M.in[1]==3){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
if(M.in[1]==4){M.draw<-round(runif(1,M.in[2],M.in[3]),3)}
Input.draws[i,3]<-M.draw
}
}
else{Input.draws[i,3]<-M.draw<-M.in[i,1]}
if(sexes==T)
{
if(length(M.in)==6)
{
M.draw.M<--1
while(M.draw.M<0)
{
if(M.in[4]==0){M.draw.M<-round(rnorm(1,M.in[5],M.in[6]),3)}
if(M.in[4]==3){M.draw.M<-round(rlnorm(1,log(M.in[5]),M.in[6]),3)}
if(M.in[4]==4){M.draw.M<-round(runif(1,M.in[5],M.in[6]),3)}
Input.draws[i,7]<-M.draw.M
}
}
if(length(M.in)>6){Input.draws[i,7]<-M.draw.M<-M.in[i,2]}
}
}
#Growth parameters
# if(Linf.in[1]>-1&k.in[1]>-1)
# {
Covar_Linf_k<-Linf.k.cor*((Linf.in[3]+0.000001)*(k.in[3]+0.000001)) #caluclate covariance matrix for multivariate sampling
Linf_k_sigma <- matrix(c((Linf.in[3]+0.000001)^2,Covar_Linf_k,Covar_Linf_k,(k.in[3]+0.000001)^2),2,2) #Set-up the variance-covariance matrix for multivariate sampling
Linf_k_samps<-rtmvnorm(1,c(Linf.in[2],k.in[2]),Linf_k_sigma,lower=rep(0,length(c(Linf.in[2],k.in[2]))))
Linf.draw<-Input.draws[i,5]<-Linf_k_samps[1]
k.draw<-Input.draws[i,6]<-Linf_k_samps[2]
#if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
t0.draw<-round(rnorm(1,t0.in[2],t0.in[3]),3)
L1.draw<-Input.draws[i,4]<-Input.draws[i,8]<-0
# }
# if(any(Linf.in[1]<1&k.in[1]<1)&!all(Linf.in[1]<1&k.in[1]<1))
# {
# # if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
# Linf.draw<-Input.draws[i,5]<-round(rnorm(1,Linf.in[2],Linf.in[3]),2)
# k.draw<-Input.draws[i,6]<-round(rnorm(1,k.in[2],k.in[3]),2)
# t0.draw<-round(rnorm(1,t0.in[2],t0.in[3]),3)
# L1.draw<-Input.draws[i,4]<-VBGF(Linf.draw,k.draw,t0.draw,0)
# }
#Male draws
if(sexes==T)
{
# if(sum(L1.in[3:4])>0) {
# L1.draw.M <- round(rnorm(1,L1.in[3], L1.in[4]), 2)
# }
# if(Linf.in[4]>-1&k.in[4]>-1)
# {
if(ctl.new$parameter_offset_approach==1)
{
Covar_Linf_k<-Linf.k.cor*((Linf.in[6]+0.000001)*(k.in[6]+0.000001)) #caluclate covariance matrix for multivariate sampling
Linf_k_sigma <- matrix(c((Linf.in[6]+0.000001)^2,Covar_Linf_k,Covar_Linf_k,(k.in[6]+0.000001)^2),2,2) #Set-up the variance-covariance matrix for multivariate sampling
Linf_k_samps<-rtmvnorm(1,c(Linf.in[5],k.in[5]),Linf_k_sigma,lower=rep(0,length(c(Linf.in[5],k.in[5]))))
Linf.draw.M<-Input.draws[i,9]<-Linf_k_samps[1]
k.draw.M<-Input.draws[i,10]<-Linf_k_samps[2]
#if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
t0.draw.M<-round(rnorm(1,t0.in[5],t0.in[6]),3)
L1.draw.M<-Input.draws[i,8]<-VBGF(Linf.draw.M,k.draw.M,t0.draw.M,t0.draw)
# }
}
if(ctl.new$parameter_offset_approach==2)
{
Linf.draw.M<-Input.draws[i,9]<-Linf.in[5]
k.draw.M<-Input.draws[i,10]<-k.in[5]
#if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
t0.draw.M<-round(rnorm(1,t0.in[5],t0.in[6]),3)
L1.draw.M<-Input.draws[i,8]<-VBGF(Linf.draw.M,k.draw.M,t0.draw.M,t0.draw)
}
#Change to offset approach
# if(any(Linf.in[4]<1&k.in[4]<1)&!all(Linf.in[4]<1&k.in[4]<1))
# {
# if(sum(L1.in[1:2])>0){L1.draw<-round(rnorm(1,L1.in[2],L1.in[3]),2); Input.draws[i,4]<-L1.draw}
# Linf.draw.M<-Input.draws.M[i,3]<-round(rnorm(1,Linf.in[5],Linf.in[6]),2)
# k.draw.M<-Input.draws.M[i,4]<-round(rnorm(1,k.in[5],k.in[6]),2)
# t0.draw>M<-round(rnorm(1,t0.in[5],t0.in[6]),3)
# L1.draw>M<-Input.draws.M[i,2]<-VBGF(Linf.draw.M,k.draw.M,t0.draw.M,0)
# }
}
#Steenpess
if(h.in[1]>=0 & length(h.in)==3 & FMSY_M.in[1]<0 & BMSY_B0.in[1]<0)
{
if(h.in[1]==1){h.draw<-round(rbeta.ab(1,h.in[2],h.in[3],0.21,0.99),3)}
if(h.in[1]==2){h.draw<-round(rbeta(1,h.in[2],h.in[3]),3)}
if(h.in[1]==10){h.draw<-round(rtnorm(1,h.in[2],h.in[3],0.21,0.99),3)}
if(h.in[1]==30){h.draw<-round(rlnormTrunc(1,log(h.in[2]),h.in[3],0.21,0.99),2)}
if(h.in[1]==4){h.draw<-round(runif(1,h.in[2],h.in[3]),2)}
Input.draws[i,1]<-h.draw
}
if(length(h.in)>3 & FMSY_M.in[1]<0 & BMSY_B0.in[1]<0){Input.draws[i,1]<-h.draw<-h.in[i]}
#Get steepness for BH from FMSY
if(h.in[1]==99 | BH_FMSY_comp==T & FMSY_M.in[1]>=0)
{
SRdat.new<-readLines(paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
x<-1
while(x==1)
{
#Draw Ms
if(M.in[1]>=0 & length(M.in)==6)
{
if(M.in[1]==0){M.draw<-round(rnorm(1,M.in[2],M.in[3]),3)}
if(M.in[1]==3){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
if(M.in[1]==4){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
Input.draws[i,3]<-M.draw
}
else{Input.draws[i,3]<-M.draw<-M.in[i,1]}
#Male draws
if(sexes==T)
{
if(length(M.in)==6)
{
if(M.in[4]==0){M.draw.M<-round(rnorm(1,M.in[5],M.in[6]),3)}
if(M.in[4]==3){M.draw.M<-round(rlnorm(1,log(M.in[5]),M.in[6]),3)}
if(M.in[4]==4){M.draw.M<-round(runif(1,M.in[5],M.in[6]),3)}
Input.draws[i,7]<-M.draw.M
}
if(length(M.in)>6){Input.draws[i,7]<-M.draw.M<-M.in[i,2]}
}
if(FMSY_M.in[1]==2){FMSY_M.draw<-round(rbeta.ab(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==10){FMSY_M.draw<-round(rtnorm(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==30){FMSY_M.draw<-round(rlnormTrunc(1,log(FMSY_M.in[2]),FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==4){FMSY_M.draw<-round(runif(1,FMSY_M.in[2],FMSY_M.in[3]),3)}
if(length(FMSY_M.in)>3){FMSY_M.draw<-FMSY_M.in[i]}
Input.draws.MQs[i,1]<-FMSY_M.draw
if(sexes==F)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw),sep=" ")
M_FMSY<-M.draw
}
if(sexes==T)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw.M),sep=" ")
M_FMSY<-mean(c(M.draw,M.draw.M))
}
SRdat.new[17]<-as.character(FMSY_M.draw*M_FMSY)
Input.draws.MQs[i,3]<-FMSY_M.draw*M_FMSY
SRdat.new[21]<-"-0.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
if(file.exists(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))==TRUE)
{
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
SRdat.new[21]<-"-0.1"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
SRdat.new[21]<-"-1.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
}
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
SRdat.new[21]<-"-0.51"
SRdat.new[23]<-"-0.31"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
}
}
if(steep_expo.out[1]>=0.25&steep_expo.out[1]<1&steep_expo.out[2]<=10&SRMQ_par[i]<0.1){x<-2}
}
}
Input.draws[i,1]<-h.draw<-steep_expo.out[1]
SR_expo.out[i]<-steep_expo.out[2]
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep=""))
}
if(FMSY_M.in[1]>=0 & BMSY_B0.in[1]>=0)
{
SRdat.new<-readLines(paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
#SRdat.new.mod<-strsplit(SRdat.new[grep("Amax",SRdat.new)+1][1]," ")[[1]]
#Amax.line<-strsplit(dat.new[grep("Number of ages",dat.new)], " ")[[1]]
#k.line<-strsplit(ctl.new[grep("VonBert_K_Fem_GP_1",ctl.new)]," ")[[1]]
#ltwt.exp.line<-strsplit(ctl.new[grep("Wtlen_2_Fem",ctl.new)]," ")[[1]]
x<-1
while(x==1)
{
#Draw Ms
if(M.in[1]>=0 & length(M.in)==6)
{
if(M.in[1]==0){M.draw<-round(rnorm(1,M.in[2],M.in[3]),3)}
if(M.in[1]==3){M.draw<-round(rlnorm(1,log(M.in[2]),M.in[3]),3)}
if(M.in[1]==4){M.draw<-round(runif(1,M.in[2],M.in[3]),3)}
Input.draws[i,3]<-M.draw
}
else{Input.draws[i,3]<-M.draw<-M.in[i,1]}
#Male draws
if(sexes==T)
{
if(length(M.in)==6)
{
if(M.in[4]==0){M.draw.M<-round(rnorm(1,M.in[5],M.in[6]),3)}
if(M.in[4]==3){M.draw.M<-round(rlnorm(1,log(M.in[5]),M.in[6]),3)}
if(M.in[4]==4){M.draw.M<-round(runif(1,M.in[5],M.in[6]),3)}
Input.draws[i,7]<-M.draw.M
}
if(length(M.in)>6){Input.draws[i,7]<-M.draw.M<-M.in[i,2]}
}
if(FMSY_M.in[1]==2){FMSY_M.draw<-round(rbeta.ab(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==10){FMSY_M.draw<-round(rtnorm(1,FMSY_M.in[2],FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==30){FMSY_M.draw<-round(rlnormTrunc(1,log(FMSY_M.in[2]),FMSY_M.in[3],0.0001,10),3)}
if(FMSY_M.in[1]==4){FMSY_M.draw<-round(runif(1,FMSY_M.in[2],FMSY_M.in[3]),3)}
if(length(FMSY_M.in)>3){FMSY_M.draw<-FMSY_M.in[i]}
if(BMSY_B0.in[1]==2){BMSY_B0.draw<-round(rbeta.ab(1,BMSY_B0.in[2],BMSY_B0.in[3],0.01,0.99),2)*100}
if(BMSY_B0.in[1]==10){BMSY_B0.draw<-round(rtnorm(1,BMSY_B0.in[2],BMSY_B0.in[3],0.01,0.99),2)*100}
if(BMSY_B0.in[1]==30){BMSY_B0.draw<-round(rlnormTrunc(1,log(BMSY_B0.in[2]),BMSY_B0.in[3],0.01,0.99),2)*100}
if(BMSY_B0.in[1]==4){BMSY_B0.draw<-round(runif(1,BMSY_B0.in[2],BMSY_B0.in[3]),2)*100}
if(length(BMSY_B0.in)>3){BMSY_B0.draw<-BMSY_B0.in[i]}
Input.draws.MQs[i,1]<-FMSY_M.draw
Input.draws.MQs[i,2]<-BMSY_B0.draw
#SRdat.new.mod[1]<-as.numeric(Amax.line[1])
if(sexes==F)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw),sep=" ")
M_FMSY<-M.draw
}
if(sexes==T)
{
SRdat.new[5]<-paste(as.character(M.draw),as.character(M.draw.M),sep=" ")
#M_FMSY<-mean(c(M.draw,M.draw.M))
M_FMSY<-M.draw
}
#print(c(M.draw,M.draw.M,M_FMSY,FMSY_M.draw,BMSY_B0.draw))
#SRdat.new.mod[3]<-Winf
#SRdat.new.mod[4]<-as.numeric(k.line[3])
#SRdat.new.mod[5]<-as.numeric(ltwt.exp.line[3])
#SRdat.new.mod[6]<-SRdat.new.mod[7]<-Amat
#SRdat.new.mod[8]<-SR_type-7 #-7 converts from SS SR_type to MQ SR_type
SRdat.new[17]<-as.character(FMSY_M.draw*M_FMSY)
SRdat.new[19]<-as.character(BMSY_B0.draw)
Input.draws.MQs[i,3]<-FMSY_M.draw*M_FMSY
#SRdat.new[grep("Amax",SRdat.new)+1][1]<-paste(SRdat.new.mod,collapse=" ")
SRdat.new[21]<-"-0.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
if(file.exists(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))==TRUE)
{
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
if(steep_expo.out[1]<0.2|steep_expo.out[1]>=2|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
#print("retry 1")
SRdat.new[21]<-"-0.1"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
#print("retry 2")
SRdat.new[21]<-"-1.51"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
}
if(steep_expo.out[1]<0.25|steep_expo.out[1]>=1|steep_expo.out[2]>10|SRMQ_par[i]>0.1)
{
#print("retry 3")
SRdat.new[21]<-"-0.51"
SRdat.new[23]<-"-0.31"
write(SRdat.new,paste(filepath,"/h_BMSY_FMSY/BMSYB0.dat",sep=""))
RUN.SS(paste(filepath,"/h_BMSY_FMSY/",sep=""), ss.exe="BMSYB0",ss.cmd=" -nox > out.txt")
SRMQ_rep.new<-readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
SRMQ_par[i]<-as.numeric(strsplit(readLines(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep="")),split=" ")[[1]][11])
steep_expo.out<-as.numeric(strsplit(SRMQ_rep.new[1]," ")[[1]])
}
}
if(steep_expo.out[1]>=0.2&steep_expo.out[1]<2&steep_expo.out[2]<=10&SRMQ_par[i]<0.1){x<-2}
}
}
Input.draws[i,1]<-h.draw<-steep_expo.out[1]
SR_expo.out[i]<-steep_expo.out[2]
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.rep",sep=""))
file.remove(paste(filepath,"/h_BMSY_FMSY/bmsyb0.par",sep=""))
}
#Three parameter S-R model
if(Zfrac.Beta.in[1]>=0)
{
if(Zfrac.Beta.in[1]==10){Zfrac.draw<-round(rtnorm(1,Zfrac.Beta.in[2],Zfrac.Beta.in[3],0,1),2)}
if(Zfrac.Beta.in[1]==30){Zfrac.draw<-round(rlnormTrunc(1,log(Zfrac.Beta.in[2]),Zfrac.Beta.in[3],0,1),2)}
if(Zfrac.Beta.in[1]==4){Zfrac.draw<-round(runif(1,Zfrac.Beta.in[2],Zfrac.Beta.in[3]),2)}
Input.draws[i,1]<-Zfrac.draw
}
if(Zfrac.Beta.in[4]>=0)
{
if(Zfrac.Beta.in[4]==10){Beta.draw<-round(rtnorm(1,Zfrac.Beta.in[5],Zfrac.Beta.in[6],0,10),2)}
if(Zfrac.Beta.in[4]==30){Beta.draw<-round(rlnormTrunc(1,log(Zfrac.Beta.in[5]),Zfrac.Beta.in[6],0,10),2)}
if(Zfrac.Beta.in[4]==4){Beta.draw<-round(runif(1,Zfrac.Beta.in[5],Zfrac.Beta.in[6]),2)}
Input.draws[i,7]<-Beta.draw
}
### Change DAT and CTL inputs ###
#change Depletion
if(Dep.in[1]>=0)
{
dat.new$CPUE[2,4]<-Dep.draw
SS_writedat(dat.new, paste0(filepath,"/",file.name[1]), overwrite=TRUE,verbose=FALSE)
}
#change M
Sys.sleep(1)
if(M.in[1]>=0){ctl.new$MG_parms[1,3]<-M.draw}
#change growth parameters
ctl.new$Growth_Age_for_L1<-t0.draw
ctl.new$MG_parms[2,3:4]<-ctl.new$MG_parms[14,3]<-0
ctl.new$MG_parms[3,3:4]<-Linf.draw
ctl.new$MG_parms[4,3:4]<-k.draw
#change male pararmeters
if(sexes==T)
{
#change M
if(M.in[4]>=0){ctl.new$MG_parms[13,3]<-M.draw.M}
#change growth parameters
ctl.new$MG_parms[14,3]<-L1.draw.M
ctl.new$MG_parms[15,3]<-Linf.draw.M
ctl.new$MG_parms[16,3]<-k.draw.M
}
ctl.new$SR_function<-SR_type
# SRtype.line<-strsplit(ctl.new[grep("SR_function",ctl.new)], " ")[[1]]
# SRtype.line[1]<-SR_type
# ctl.new[grep("SR_function",ctl.new)]<-paste(SRtype.line,collapse=" ")
#change R0
if(R_start[1]==0)
{
R0.draw<-R_start[2]
ctl.new$SR_parms[1,3:4]<-R0.draw
# R0.line<-strsplit(ctl.new[grep("R0",ctl.new)], " ")[[1]]
# R0.line[c(3,4)]<-R0.draw
# ctl.new[grep("R0",ctl.new)]<-paste(R0.line,collapse=" ")
}
if(R_start[1]==1)
{
R0.draw<-round(rlnormTrunc(1,log(R_start[2]),0.5,3,15),2)
ctl.new$SR_parms[1,3:4]<-R0.draw
# R0.line<-strsplit(ctl.new[grep("R0",ctl.new)], " ")[[1]]
# R0.line[c(3,4)]<-R0.draw
# ctl.new[grep("R0",ctl.new)]<-paste(R0.line,collapse=" ")
}
#change h
if(SR_type==3)
{
ctl.new$SR_parms[2,3:4]<-h.draw
# h.line<-strsplit(ctl.new[grep("SR_BH_steep",ctl.new)], " ")[[1]]
# h.line[c(3,4)]<-h.draw
# ctl.new[grep("SR_BH_steep",ctl.new)]<-paste(h.line,collapse=" ")
}
if(SR_type==9)
{
ctl.new$SR_parms[2,3:4]<-h.draw
ctl.new$SR_parms[3,3:4]<-SR_expo.out[i]
h.line<-strsplit(ctl.new[grep("SR_RkrPower_steep",ctl.new)], " ")[[1]]
# h.line[c(3,4)]<-h.draw
# ctl.new[grep("SR_RkrPower_steep",ctl.new)]<-paste(h.line,collapse=" ")
# SRexpo.line<-strsplit(ctl.new[grep("SR_RkrPower_gamma",ctl.new)], " ")[[1]]
# SRexpo.line[c(3,4)]<-SR_expo.out[i]
# ctl.new[grep("SR_RkrPower_gamma",ctl.new)]<-paste(SRexpo.line,collapse=" ")
}
#change Sfrac and Beta
if(SR_type==7 & Zfrac.Beta.in[1]>=0)
{
ctl.new$SR_parms[2,3:4]<-Zfrac.draw
# Zfrac.line<-strsplit(ctl.new[grep("Zfrac",ctl.new)], " ")[[1]]
# Zfrac.line[c(3,4)]<-Zfrac.draw
# ctl.new[grep("Zfrac",ctl.new)]<-paste(Zfrac.line,collapse=" ")
}
if(SR_type==7 & Zfrac.Beta.in[4]>=0)
{
ctl.new$SR_parms[3,3:4]<-Beta.draw
# Beta.line<-strsplit(ctl.new[grep("Beta",ctl.new)], " ")[[1]]
# Beta.line[c(3,4)]<-Beta.draw
# ctl.new[grep("Beta",ctl.new)]<-paste(Beta.line,collapse=" ")
}
#write(ctl.new,paste(filepath,"/",file.name[2],sep=""))
SS_writectl(ctl.new, paste0(filepath,"/",file.name[2]), overwrite=TRUE,verbose=FALSE)
#Run model
RUN.SS(paste(filepath,"/",sep=""),ss.cmd=" -nohess -nox > out.txt 2>&1",OS.in=OStype,show_in_console = FALSE)
#if(OStype=="Windows"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
#if(OStype=="OSX_Linux"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="./ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
#Evaluate convergence and record values
rep.new<-readLines(paste(filepath,"/Report.sso",sep=""))
Sys.sleep(0.5)
#The R0 loop. Makes sure R0 is changing from input and that the depletion match is happening. Ensures good models are kept.
#Profiles over R0 when intial run fails.
if(doR0.loop[1]>0)
{
xx<-1
R0.explore<-seq(doR0.loop[2],doR0.loop[3],doR0.loop[4])
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==TRUE|as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])==0){Dep.out.testR0<-10}
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==FALSE&as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])>0){Dep.out.testR0<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
# while(abs(as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][3])-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][8]))==0)
while(abs(Dep.out.testR0-Dep.draw)>0.01)
{
print(paste0("IN THE LOOP: R0= ",R0.explore[xx]))
ctl.new$SR_parms[1,3:4]<-R0.explore[xx]
SS_writectl(ctl.new, paste0(filepath,"/",file.name[2]), overwrite=TRUE,verbose=FALSE)
# ctl.new<-readLines(paste(filepath,"/",file.name[2],sep=""))
# R0.line<-strsplit(ctl.new[grep("R0",ctl.new)], " ")[[1]]
# R0.line[c(3,4)]<-R0.explore[xx]
# ctl.new[grep("R0",ctl.new)]<-paste(R0.line,collapse=" ")
# write(ctl.new,paste(filepath,"/",file.name[2],sep=""))
RUN.SS(paste(filepath,"/",sep=""),ss.cmd=" -nohess -nox > out.txt 2>&1",OS.in=OStype)
# if(OStype=="Windows"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
# if(OStype=="OSX_Linux"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="./ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
rep.new<-readLines(paste(filepath,"/Report.sso",sep=""))
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==TRUE|as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])==0){Dep.out.testR0<-10}
if(is.na(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2]))==FALSE&as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])>0){Dep.out.testR0<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
print(paste0("Is depletion difference < 0.01? ",(abs(as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])-Dep.draw)<0.01)))
xx<-xx+1
if(xx==length(R0.explore)){break}
}
}
SSS.output.list[[i]]<-SS_output(paste(filepath,"/",sep=""), covar=FALSE,verbose=FALSE,printstats=FALSE)
#Begin extracting info from run
#forecast.file<-readLines(paste(filepath,"/Forecast-report.sso",sep=""))
# for(iii in 1:length(ts_yrs[1]:ts_yrs[2]))
# {
# SB.out[i,iii]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",sb.years[iii],sep=""),rep.new)], " ")[[1]][3])
# TB.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][5])
# SumAge.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][6])
# SPR.out[i,iii]<-1-as.numeric(strsplit(rep.new[grep(paste("SPRratio_",sb.years[iii]+1,sep=""),rep.new)], " ")[[1]][3])
# # B_BMSY.out[i,iii]<-as.numeric(strsplit(rep.new[grep("Yr B/Bmsy F/Fmsy",rep.new)+iii], " ")[[1]][2])
# }
if(i==1)
{
SB.out<-TB.out<-SumAge.out<-SPR.out<-as.data.frame(matrix(NA,nrow=reps,ncol=length(SSS.output.list[[i]]$timeseries$SpawnBio)))
colnames(SB.out)<-colnames(TB.out)<-colnames(SumAge.out)<-colnames(SPR.out)<-sb.years
}
SB.out[i,]<-SSS.output.list[[i]]$timeseries$SpawnBio
TB.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_all
SumAge.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_smry
#SPR.out[i,]<-1-as.numeric(strsplit(rep.new[grep(paste("SPRratio_",sb.years[iii]+1,sep=""),rep.new)], " ")[[1]][3])
nforecasts<- SSS.output.list[[i]]$nforecastyears
OFL.out[[i]]<-SSS.output.list[[i]]$derived_quants[paste0("OFLCatch_", c((ts_yrs[2]+1):(ts_yrs[2]+nforecasts))), "Value"]
ABC.out[[i]]<-SSS.output.list[[i]]$derived_quants[paste0("ForeCatch_", c((ts_yrs[2]+1):(ts_yrs[2]+nforecasts))), "Value"]
Dep.series.out<-SB.out/SB.out[,1]
#Quant.out[i,1]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,2]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,3]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,4]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Fem_GP_1",rep.new)], " ")[[1]][3])
Quant.out[i,1]<-SSS.output.list[[i]]$parameters[1,3]
Quant.out[i,2]<-SSS.output.list[[i]]$parameters[2,3]
Quant.out[i,3]<-SSS.output.list[[i]]$parameters[3,3]
Quant.out[i,4]<-SSS.output.list[[i]]$parameters[4,3]
if(sexes==T)
{
# Quant.out[i,5]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,6]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,7]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,8]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Mal_GP_1",rep.new)], " ")[[1]][3])
Quant.out[i,5]<-SSS.output.list[[i]]$parameters[13,3]
Quant.out[i,6]<-SSS.output.list[[i]]$parameters[14,3]
Quant.out[i,7]<-SSS.output.list[[i]]$parameters[15,3]
Quant.out[i,8]<-SSS.output.list[[i]]$parameters[16,3]
}
# if(SR_type==3){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_BH_steep",rep.new)], " ")[[1]][3])}
if(SR_type==3){Quant.out[i,9]<-SSS.output.list[[i]]$parameters[24,3]}
if(SR_type==7){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)-1], " ")[[1]][2])}
if(SR_type==7){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)], " ")[[1]][2])}
if(SR_type==8){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_steep",rep.new)], " ")[[1]][2])}
if(SR_type==8){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_gamma",rep.new)], " ")[[1]][2])}
#Quant.out[i,10]<-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][3])
Quant.out[i,10]<-SSS.output.list[[i]]$parameters[23,3]
#Quant.out[i,11]<-as.numeric(strsplit(rep.new[grep("SSB_Virgin",rep.new)], " ")[[1]][3])
Quant.out[i,11]<-SSS.output.list[[i]]$SBzero
Quant.out[i,12]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("SSB_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,13]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("Bratio_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,14]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("SSB_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,15]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("B_MSY/SSB_unfished",SSS.output.list[[i]]$derived_quants[,1]),][2])
if(!is.na(as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2]))){Quant.out[i,16]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])}
Quant.out[i,17]<-SSS.output.list[[i]]$likelihoods_used[1,1]
Quant.out[i,18]<-SSS.output.list[[i]]$likelihoods_used[4,1]
Quant.out[i,19]<-SSS.output.list[[i]]$maximum_gradient_component
if(Dep.in[1]>=0){Quant.out[i,20]<-as.numeric(dat.new$CPUE[2,4])}
if(Dep.in[1]>=0){Quant.out[i,21]<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
Quant.out[i,22]<-SSS.output.list[[i]]$parameters[23,8]
Quant.out[i,23]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("F_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])/as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,24]<-SSS.output.list[[i]]$likelihoods_used[11,1]
warnings.out<-readLines(paste(filepath,"/warning.sso",sep=""))
num.warning<-as.numeric(strsplit(warnings.out[grep("N warnings",warnings.out)[1]], " ")[[1]][4])
#Quant.out[i,12]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])
#Quant.out[i,13]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("SSB_Virgin",rep.new)], " ")[[1]][3])
#Quant.out[i,14]<-as.numeric(strsplit(rep.new[grep("SSB_MSY",rep.new)], " ")[[1]][2])/as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
#Quant.out[i,15]<-as.numeric(strsplit(rep.new[grep("Bmsy/Bzero",rep.new)], " ")[[1]][1])
#Quant.out[i,17]<-as.numeric(strsplit(forecast.file[grep("calculate_FMSY",forecast.file)+13],split="[[:blank:]]+")[[1]][2])/as.numeric(strsplit(forecast.file[grep("BIO_Smry_unfished",forecast.file)],split="[[:blank:]]+")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2]))){Quant.out[i,16]<-as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])}
# Quant.out[i,17]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)], " ")[[1]][2])
# Quant.out[i,18]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)+2], " ")[[1]][2])
# Quant.out[i,19]<-as.numeric(strsplit(rep.new[grep("Convergence",rep.new)], " ")[[1]][2])
# Quant.out[i,23]<-as.numeric(strsplit(rep.new[grep(paste("F_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,26]<-as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,27]<-as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
#Quant.out[i,24]<-as.numeric(strsplit(rep.new[grep("Crash_Pen",rep.new)], " ")[[1]][2])
if(Dep.in[1]>=0)
{
if(FMSY_M.in[1]<0 & BMSY_B0.in[1]<0&abs(Quant.out[i,20]-Quant.out[i,21])<0.01&Quant.out[i,24]<0.0001){i<-i+1}
if(FMSY_M.in[1]>=0 & BMSY_B0.in[1]>=0 | h.in[1]==99 | BH_FMSY_comp==T)
{
if(h.in[1]==99)
{
if(abs(Quant.out[i,20]-Quant.out[i,21])<0.01&Quant.out[i,24]<0.0001&(abs(Quant.out[i,14]-Input.draws.MQs[i,3])/Input.draws.MQs[i,3])<0.6){i<-i+1}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
if(length(M.in)>6|length(Dep.in)>3|length(h.in)>3|length(FMSY_M.in)>3|length(BMSY_B0.in)>3)
{
Quant.out[i,]<-NA
i<-i+1
}
}
}
else
{
if(abs(Quant.out[i,20]-Quant.out[i,21])<0.01&Quant.out[i,24]<0.0001&(abs(Quant.out[i,14]-(Input.draws.MQs[i,2])/100)/(Input.draws.MQs[i,2]/100))<0.5){i<-i+1}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
if(length(M.in)>6|length(Dep.in)>3|length(h.in)>3|length(FMSY_M.in)>3|length(BMSY_B0.in)>3)
{
Quant.out[i,]<-NA
i<-i+1
}
}
}
}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
if(length(M.in)>6|length(Dep.in)>3|length(h.in)>3|length(FMSY_M.in)>3|length(BMSY_B0.in)>3)
{
Quant.out[i,]<-NA
i<-i+1
}
}
# }
}
if(Dep.in[1]<0)
{
if(Quant.out[i,24]<0.0001){i<-i+1}
else
{
#Run using .par file if above conditions not met
print("*** RUNNING WITH .PAR ***")
starter.new<-readLines(paste(filepath,"/starter.ss",sep=""))
par_line<-strsplit(starter.new[grep("ss.par",starter.new)], " ")[[1]]
par_line[1]<-1
starter.new[grep("ss.par",starter.new)]<-paste(par_line, collapse=" ")
write(starter.new,paste(filepath,"/starter.ss",sep=""))
RUN.SS(paste(filepath,"/",sep=""),ss.cmd=" -nohess -nox > out.txt 2>&1",OS.in=OStype)
# if(OStype=="Windows"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
# if(OStype=="OSX_Linux"){RUN.SS(paste(filepath,"/",sep=""), ss.exe="./ss",ss.cmd=" -nohess -nox > out.txt 2>&1")}
#RUN.SS(paste(filepath,"/",sep=""), ss.exe="ss",ss.cmd=" -nohess -nox > out.txt 2>&1")
rep.new<-readLines(paste(filepath,"/Report.sso",sep=""))
SSS.output.list[[i]]<-SS_output(paste(filepath,"/",sep=""), covar=FALSE,verbose=FALSE,printstats=FALSE)
#forecast.file<-readLines(paste(filepath,"/Forecast-report.sso",sep=""))
# for(iii in 1:length(ts_yrs[1]:ts_yrs[2]))
# {
# SB.out[i,iii]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",sb.years[iii],sep=""),rep.new)], " ")[[1]][3])
# TB.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][5])
# SumAge.out[i,iii]<-as.numeric(strsplit(rep.new[grep("TIME_SERIES",rep.new)+3+iii], " ")[[2]][6])
# SPR.out[i,iii]<-1-as.numeric(strsplit(rep.new[grep(paste("SPRratio_",sb.years[iii]+1,sep=""),rep.new)], " ")[[1]][3])
# # B_BMSY.out[i,iii]<-as.numeric(strsplit(rep.new[grep("Yr B/Bmsy F/Fmsy",rep.new)+iii], " ")[[1]][2])
# }
SB.out[i,]<-SSS.output.list[[i]]$timeseries$SpawnBio
TB.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_all
SumAge.out[i,]<-SSS.output.list[[i]]$timeseries$Bio_smry
Dep.series.out<-SB.out/SB.out[,1]
Quant.out[i,1]<-SSS.output.list[[i]]$parameters[1,3]
Quant.out[i,2]<-SSS.output.list[[i]]$parameters[2,3]
Quant.out[i,3]<-SSS.output.list[[i]]$parameters[3,3]
Quant.out[i,4]<-SSS.output.list[[i]]$parameters[4,3]
# Quant.out[i,1]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,2]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,3]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Fem_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,4]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Fem_GP_1",rep.new)], " ")[[1]][3])
if(sexes==T)
{
Quant.out[i,5]<-SSS.output.list[[i]]$parameters[13,3]
Quant.out[i,6]<-SSS.output.list[[i]]$parameters[14,3]
Quant.out[i,7]<-SSS.output.list[[i]]$parameters[15,3]
Quant.out[i,8]<-SSS.output.list[[i]]$parameters[16,3]
# Quant.out[i,5]<-as.numeric(strsplit(rep.new[grep("NatM_p_1_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,6]<-as.numeric(strsplit(rep.new[grep("L_at_Amin_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,7]<-as.numeric(strsplit(rep.new[grep("L_at_Amax_Mal_GP_1",rep.new)], " ")[[1]][3])
# Quant.out[i,8]<-as.numeric(strsplit(rep.new[grep("VonBert_K_Mal_GP_1",rep.new)], " ")[[1]][3])
}
#if(SR_type==3){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_BH_steep",rep.new)], " ")[[1]][3])}
if(SR_type==3){Quant.out[i,9]<-SSS.output.list[[i]]$parameters[24,3]}
if(SR_type==7){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)-1], " ")[[1]][3])}
if(SR_type==7){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_surv_Beta",rep.new)], " ")[[1]][3])}
if(SR_type==8){Quant.out[i,9]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_steep",rep.new)], " ")[[1]][3])}
if(SR_type==8){Quant.out[i,29]<-as.numeric(strsplit(rep.new[grep("SR_RkrPower_gamma",rep.new)], " ")[[1]][3])}
Quant.out[i,10]<-SSS.output.list[[i]]$parameters[23,3]
Quant.out[i,11]<-SSS.output.list[[i]]$SBzero
Quant.out[i,12]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("SSB_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,13]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("Bratio_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,14]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("SSB_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,15]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("B_MSY/SSB_unfished",SSS.output.list[[i]]$derived_quants[,1]),][2])
if(!is.na(as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2]))){Quant.out[i,16]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])}
Quant.out[i,17]<-SSS.output.list[[i]]$likelihoods_used[1,1]
Quant.out[i,18]<-SSS.output.list[[i]]$likelihoods_used[4,1]
Quant.out[i,19]<-SSS.output.list[[i]]$maximum_gradient_component
if(Dep.in[1]>=0){Quant.out[i,20]<-as.numeric(dat.new$CPUE[2,4])}
if(Dep.in[1]>=0){Quant.out[i,21]<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][2])}
Quant.out[i,22]<-SSS.output.list[[i]]$parameters[23,8]
Quant.out[i,23]<-as.numeric(SSS.output.list[[i]]$derived_quants[grep(paste("F_",ts_yrs[2],sep=""),SSS.output.list[[i]]$derived_quants[,1]),][2])/as.numeric(SSS.output.list[[i]]$derived_quants[grep("annF_MSY",SSS.output.list[[i]]$derived_quants[,1]),][2])
Quant.out[i,24]<-SSS.output.list[[i]]$likelihoods_used[11,1]
# Quant.out[i,10]<-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][3])
# Quant.out[i,11]<-as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
# Quant.out[i,12]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])
# Quant.out[i,13]<-as.numeric(strsplit(rep.new[grep(paste("SSB_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
# #Quant.out[i,14]<-as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+1,sep=""),rep.new)], " ")[[1]][2])
# #Quant.out[i,15]<-as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+1,sep=""),rep.new)], " ")[[1]][2])
# Quant.out[i,14]<-as.numeric(strsplit(rep.new[grep("SSB_MSY",rep.new)], " ")[[1]][2])/as.numeric(strsplit(rep.new[grep("SSB_Initial",rep.new)], " ")[[1]][3])
# Quant.out[i,15]<-as.numeric(strsplit(rep.new[grep("Bmsy/Bzero",rep.new)], " ")[[1]][1])
# #Quant.out[i,17]<-as.numeric(strsplit(forecast.file[grep("Maximum_Sustainable_Yield;_where_Yield_is_Dead_Catch",forecast.file)+14],split="[[:blank:]]+")[[1]][2])/as.numeric(strsplit(forecast.file[grep("BIO_Smry_unfished",forecast.file)],split="[[:blank:]]+")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2]))){Quant.out[i,16]<-as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])}
# Quant.out[i,17]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)], " ")[[1]][2])
# Quant.out[i,18]<-as.numeric(strsplit(rep.new[grep("TOTAL",rep.new)+2], " ")[[1]][2])
# Quant.out[i,19]<-as.numeric(strsplit(rep.new[grep("Convergence",rep.new)], " ")[[1]][2])
# Quant.out[i,20]<-as.numeric(dat.new$CPUE[2,4])
# Quant.out[i,21]<-as.numeric(strsplit(rep.new[grep(paste("Bratio_",as.numeric(dat.new$CPUE[2,1]),sep=""),rep.new)], " ")[[1]][3])
# Quant.out[i,22]<-as.numeric(strsplit(rep.new[grep("R0",rep.new)], " ")[[1]][8])
# Quant.out[i,23]<-as.numeric(strsplit(rep.new[grep(paste("F_",ts_yrs[2],sep=""),rep.new)], " ")[[1]][3])/as.numeric(strsplit(rep.new[grep("annF_MSY",rep.new)], " ")[[1]][2])
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,26]<-as.numeric(strsplit(rep.new[grep(paste("OFLCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
# if(!is.na(as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2]))){Quant.out[i,27]<-as.numeric(strsplit(rep.new[grep(paste("ForeCatch_",ts_yrs[2]+2,sep=""),rep.new)], " ")[[1]][2])}
# Quant.out[i,24]<-as.numeric(strsplit(rep.new[grep("Crash_Pen",rep.new)], " ")[[1]][2])
starter.new<-readLines(paste(filepath,"/starter.ss",sep=""))
par_line<-strsplit(starter.new[grep("ss.par",starter.new)], " ")[[1]]
par_line[1]<-0
starter.new[grep("ss.par",starter.new)]<-paste(par_line, collapse=" ")
write(starter.new,paste(filepath,"/starter.ss",sep=""))
if(Quant.out[i,24]<0.0001){i<-i+1}
else
{
if(n.bad==0){Quant.out.bad<-Quant.out[i,]}
if(n.bad>0){Quant.out.bad<-rbind(Quant.out.bad,Quant.out[i,])}
n.bad<-n.bad+1
}
}
}
print(paste0("Kept ",i-1))
ii<-ii+1
if(file.exists(paste(filepath,"/Report.sso",sep=""))){file.remove(paste(filepath,"/Forecast-report.sso",sep=""))}
if(file.exists(paste(filepath,"/Report.sso",sep=""))){file.remove(paste(filepath,"/Report.sso",sep=""))}
colnames(Quant.out)<-colnames(Quant.out.bad)<-c("M_f","L1_f","Linf_f","k_f","M_m","L1_m","Linf_m","k_m","h","R0","SB0",paste("SSB_",ts_yrs[2],sep=""),"Term_Yr_Dep","SBMSY","SBMSY/SB0","FMSY","-lnL","LL_survey","Gradient","Dep.Obs","Dep.Exp","R0_init",paste("F_",ts_yrs[2],"/FMSY",sep=""),"Crash_penalty","Rick_gamma")
if(ncol(Input.draws)==10){colnames(Input.draws)<-c("h","Dep","M_f","L1_f","Linf_f","k_f","M_m","L1_m","Linf_m","k_m")}
if(ncol(Input.draws)==11){colnames(Input.draws)<-c("Sfrac","Dep","M_f","L1_f","Linf_f","k_f","Beta","M_m","L1_m","Linf_m","k_m")}
if(SR_type>=8 | h.in[1]==99 | BH_FMSY_comp==T)
{
Input.draws<-cbind(Input.draws,Input.draws.MQs,SR_expo.out,SRMQ_par)
ltcolnames<-length(colnames(Input.draws))
colnames(Input.draws)[c(ltcolnames-1,ltcolnames)]<-c("Beta","Obj_fxn")
}
end.time<-Sys.time()
OFL.out.save<-do.call("rbind",OFL.out)
ABC.out.save<-do.call("rbind",ABC.out)
colnames(OFL.out.save)<-colnames(ABC.out.save)<-c((ts_yrs[2]+1):(ts_yrs[2]+nforecasts))
Spp.quant.out<-list(Input.draws,Quant.out,SB.out,Dep.series.out,TB.out,SumAge.out,OFL.out.save,ABC.out.save,Quant.out.bad,ii-1,(as.numeric(end.time)-as.numeric(start.time))/60)
names(Spp.quant.out)<-c("Priors","Posteriors","SB_series","Rel_Stock_status_series","Total_Biomass","Summary_Biomass","OFL","ABC","Rejected_draws","Total draws","Runtime_minutes")
SSS.out<-Spp.quant.out
save(SSS.output.list,file=paste(filepath,"/SSSoutput.DMP",sep=""))
save(SSS.out,file=paste(filepath,"/SSS_out.DMP",sep=""))
}
return(Spp.quant.out)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.