Put some intro-text here.
#FPRESS.run.r
#Andrew Campbell
#Marine Institite, Galway
#andrew.campbell@marine.ie
#V1 10/07/2006
#V2 28/08/2006 completed comments, some minor modifications to make
# code more readable
#V2 10/10/2006 tacopt,simtype,tacmin,tacmax,tacbias,taccv,fmin,fmax,fbarcv,
# fbias replaced with variable,resolution,varmin,varmax,varcv,varbias
#V2 16/10/2006 function changed to FPRESS.Run(), runref missing, output
# directory changed to outdata
#
# Einar Hjörleifsson
# Marines Research Institute, Reykavík
# einarhj@hafro.is
#
# FPRESS obtained from Andrew Campell on 12.8.2009
# Relative to the original run.r script, only minor changes have been made to the
# actual executable part of the program. However, major restructuring has been
# done on programming blocks and the commenting.
# Track on some critical bugs that have been discovered in the original run.r
# are kept on the DokuWiki-site fishvice.hafro.is
#
# 10 20 30 40 50 60 70 80 90
#234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234
FPRESS.Run<-function(runref="missing",batch=0,debug=0)
{
# runref - the simulation reference. Function will examine file
# options<runref>.xml
# for model parameters
# batch - batch flag. If true(1), the function will not prompt user for
# any settings and will use defaults.
# debug flag - if true(1), print debug information to the screen
if (!(is.numeric(batch)))
stop("FPRESS.Run::if batch argument is supplied it must assume a value of 0 or 1\n")
if (!(is.numeric(debug)))
stop("FPRESS.Run::if debug argument is supplied it must assume a value of 0 or 1\n")
if (!(batch==1) & !(batch==0))
stop("FPRESS.Run::if batch argument is supplied it must assume a value of 0 or 1\n")
if (!(debug==1) & !(debug==0))
stop("FPRESS.Run::if debug argument is supplied it must assume a value of 0 or 1\n")
currwd<-getwd() #save the current working directory
home <- FPRESS.Home() #get the FPRESS home directory
if (debug==1) {cat("home=",home,"\n")}
if (nchar(home)==0)
stop("Exiting - no FPRESSHome environment variable")
#prompt for the run reference (only if not in batch mode)
if (batch==0) {
if (missing(runref)) {
runref<-winDialogString("Run Reference?","")
}
} else {
if (missing(runref)) stop("No run reference")
}
if (nchar(runref)==0) stop("Exiting - no run reference provided")
##############################################################################
# READ IN THE OPTION FILE
#
# This file contains simulation parameters and controls with regards to
# the fisheries and management and user controls of output files
# The file is located in the folder options in the path given by
# FPRESS.Home. This section contains the following steps:
# 1. Check if the file exists. If not stop, otherwise read the.xml file
# 2. Assign the individual parameters in the file to a local variable name
# are assigned to a local r variable name
# 3. Generate some illegible computer ontrol with regards user requested
# output.
# 4. Some code to adapt the generic FPRESS to iCod HCR evaluation
# 5. Give some user feedback, if requested
#
# Some explanation of each variable in the control file:
# datarecyear - the final year of historical recruitment data contained in
# recfile that is to be included in the simulation (including
# the bootstrapping options). This parameter permits the use of
# only a portion of the historical recruitment data. If greater
# than the final year's data then all data in the file will
# be used
# varcv - the coefficient of variation to used when adding in noise to the
# controlling variable
# varbias - the amount of bias to apply when randomising the controlling
# variable
# varmin - the minimum (first) or only (for resolution=0) value to use
# for the controlling variable
# varmax - for simulations with resolution>0, he maximum (final) value to
# use for the controlling variable
# assessdelay - this parameter controls the timing of assessment
# events. The assessment will be based on simulated data
# corresponding to the current year less the number of years
# given in assessdelay.
# hcrchange - the proportional change to TAC of F (as a result of an HCR)
# permitted in a single year
# ssbhcr - reference SSB level used by some HCRs for comparison with
# simulated SSB levels
# fbarhcr - reference fbar value used by some HCRs for comparison with
# simulated fbar values
# minTAC - if a minimum TAC must be taken then this value will be greater
# than zero. Applies only with variable=0
# hcrperiod - (in years) between harvest control rule actions being applied
# to the assessment
# OTHER VARIABLES - see in-line comments below
# 1. Read in the whole file ##################################################
#
optionsfile = paste(FPRESS.Home(),"\\options\\","options",runref,".xml",sep="")
#check options file exists
if (!(file.exists(optionsfile))) {
stop(paste("\nOptions file",optionsfile,"doesn't exist\nUnable to continue\n",sep=" "))}
doc<-xmlTreeParse(optionsfile) # read in the file
# 2. Pick individual parameters and assign them to a local variable name #####
#
root<-xmlRoot(doc) #root xmlNode object
stockname<-nodeval(root,"stockname","") # Stock name etc.
indata.path<-paste(FPRESS.Home(),"\\indata\\",sep="") # path to input data
outdata.path<-paste(FPRESS.Home(),"\\outdata\\",sep="")# path to output
datfile<-nodeval(root,"datfile","") # age data file
recfile<-nodeval(root,"recfile","") # historical R & SBB
catchfile<-nodeval(root,"catchfile","") # historical Y & F
minage<-as.integer(nodeval(root,"minage",1)) # minimum age
maxage<-as.integer(nodeval(root,"maxage",1)) # max age, treated as a +group
fbarmin<-as.integer(nodeval(root,"fbarmin",0)) # lower age class for fbar
fbarmax<-as.integer(nodeval(root,"fbarmax",1)) # upper age class for fbar
ssbtest1<-as.integer(nodeval(root,"ssbtest1",0)) # virtual ssb test level 1 (e.g. Blim)
ssbtest2<-as.integer(nodeval(root,"ssbtest2",0)) # vitural ssb text level 2 (e.g. Bpa)
resolution<-as.integer(nodeval(root,"resolution",0)) # number of steps of F or TAC to be tested
startyear<-as.integer(nodeval(root,"startyear",1)) # simulation start year
years<-as.integer(nodeval(root,"years",1)) # number of simulation years
iterations<-as.integer(nodeval(root,"iterations",1)) # simulation number per resolution
histrec<-as.integer(nodeval(root,"histrec",0)) # flag, if 1 use historical recruits
histreccv<-as.numeric(nodeval(root,"histreccv",0)) # historical recruit cv
histcatch<-as.integer(nodeval(root,"histcatch",0)) # flag, if 1 then use historical catch/fbar
histcatchcv<-as.numeric(nodeval(root,"histcatchcv",0)) # historical catch/fbar cv
rec<-as.integer(nodeval(root,"rec",0)) # recruitment model number
recmean<-as.numeric(nodeval(root,"recmean",0)) # mean recruitment (sometimes used)
reccv<-as.numeric(nodeval(root,"reccv",0)) # recruitment cv
seggrad<-as.numeric(nodeval(root,"seggrad",1)) # slope of the segmented regression line
ssbcut<-as.numeric(nodeval(root,"ssbcut",0)) # SSB (S*) break point value in the segm. regression
datarecyear<-as.integer(nodeval(root,"datarecyear",0)) # See comments at start
variable<-as.numeric(nodeval(root,"variable",0)) # TAC(0) or F(1) based simulation
varcv<-as.numeric(nodeval(root,"varcv",0))
varbias<-as.numeric(nodeval(root,"varbias",1))
varmin<-as.numeric(nodeval(root,"varmin",1))
varmax<-as.numeric(nodeval(root,"varmax",1))
assessdelay<-as.integer(nodeval(root,"assessdelay",1))
ssbassessbias<-as.numeric(nodeval(root,"ssbassessbias",1)) # ssb assessment bias
ssbassesscv<-as.numeric(nodeval(root,"ssbassesscv",0)) # ssb assessment cv
fbarassessbias<-as.numeric(nodeval(root,"fbarassessbias",1)) # F assessment bias
fbarassesscv<-as.numeric(nodeval(root,"fbarassesscv",0)) # F assessment cv
hcrrule<-as.integer(nodeval(root,"hcrrule",0)) # hcrrule
hcrperiod<-as.integer(nodeval(root,"hcrperiod",1))
hcrchange<-as.numeric(nodeval(root,"hcrchange",1))
ssbhcr<-as.numeric(nodeval(root,"ssbhcr",0))
fbarhcr<-as.numeric(nodeval(root,"fbarhcr",0))
mintac<-as.numeric(nodeval(root,"minimumtac",0))
opftacmult<-as.numeric(nodeval(root,"opftacmult",0)) # ftac multiplier
oprecruits<-as.numeric(nodeval(root,"oprecruits",1)) # recruits
oppop<-as.numeric(nodeval(root,"oppop",0)) # population (Nay)
opftac<-as.numeric(nodeval(root,"opftac",0)) # opftac (F or TAC)
opf<-as.numeric(nodeval(root,"opf",0)) # opf
opfmult<-as.numeric(nodeval(root,"opfmult",0)) # F multiplier
opspawningmass<-as.numeric(nodeval(root,"opspawningmass",1))# SSB
opfbar<-as.numeric(nodeval(root,"opfbar",1)) # Fbar
opcatch<-as.numeric(nodeval(root,"opcatch",0)) # Catch (numbers)
opcatchweight<-as.numeric(nodeval(root,"opcatchweight",1)) # Catch (weight)
opyestock<-as.numeric(nodeval(root,"opyestock",0)) # stock in final year
opindex<-as.numeric(nodeval(root,"opindex",0)) # ?????
# 3. Code related to controling file output #################################
#
savedobjs<-c()
if (opindex>=1024){savedobjs<-c(savedobjs,"YearEndStock"); opindex<-opindex-1024}
if (opindex>=512) {savedobjs<-c(savedobjs,"CatchWeight"); opindex<-opindex-512}
if (opindex>=256) {savedobjs<-c(savedobjs,"Catch"); opindex<-opindex-256}
if (opindex>=128) {savedobjs<-c(savedobjs,"FBar"); opindex<-opindex-128}
if (opindex>=64) {savedobjs<-c(savedobjs,"NewSpawnMass"); opindex<-opindex-64}
if (opindex>=32) {savedobjs<-c(savedobjs,"FishingMortalityMult"); opindex<-opindex-32}
if (opindex>=16) {savedobjs<-c(savedobjs,"FishingMortality"); opindex<-opindex-16}
if (opindex>=8) {savedobjs<-c(savedobjs,"FTac"); opindex<-opindex-8}
if (opindex>=4) {savedobjs<-c(savedobjs,"Pop"); opindex<-opindex-4}
if (opindex>=2) {savedobjs<-c(savedobjs,"Recruits"); opindex<-opindex-2}
if (opindex>=1) {savedobjs<-c(savedobjs,"FTacMult"); opindex<-opindex-1}
if (debug==1) {cat("savedobjs=",savedobjs,"\n")}
# 4. Some shitmix related to iCod adaptation
# If variable=0 (TAC fisheries) and varmin < 1 then on assumes that
# varmin (and varmax) stands for harvest rate.
#
if(variable==0 & varmin<1) (variable<-2)
# 5. Give user some feedback, if he requested ###############################
#
if (debug==1) {cat("stockname=",stockname,"\n",
"indata.path=",indata.path,"\n",
"outdata.path=",outdata.path,"\n",
"datfile=",datfile,"\n",
"recfile=",recfile,"\n",
"catchfile=",catchfile,"\n",
"minage=",minage,"\n",
"maxage=",maxage,"\n",
"fbarmin=",fbarmin,"\n",
"fbarmax=",fbarmax,"\n",
"ssbtest1=",ssbtest1,"\n",
"ssbtest2=",ssbtest2,"\n",
"resolution=",resolution,"\n",
"startyear=",startyear,"\n",
"years=",years,"\n",
"iterations=",iterations,"\n",
"histrec=",histrec,"\n",
"histreccv=",histreccv,"\n",
"histcatch=",histcatch,"\n",
"histcatchcv=",histcatchcv,"\n",
"rec=",rec,"\n",
"recmean=",recmean,"\n",
"reccv=",reccv,"\n",
"seggrad=",seggrad,"\n",
"ssbcut=",ssbcut,"\n",
"datarecyear=",datarecyear,"\n",
"variable=",variable,"\n",
"varcv=",varcv,"\n",
"varbias=",varbias,"\n",
"varmin=",varmin,"\n",
"varmax=",varmax,"\n",
"assessdelay=",assessdelay,"\n",
"ssbassessbias=",ssbassessbias,"\n",
"ssbassesscv=",ssbassesscv,"\n",
"fbarassessbias=",fbarassessbias,"\n",
"fbarassesscv=",fbarassesscv,"\n",
"hcrrule=",hcrrule,"\n",
"hcrperiod=",hcrperiod,"\n",
"hcrchange=",hcrchange,"\n",
"ssbhcr=",ssbhcr,"\n",
"fbarhcr=",fbarhcr,"\n",
"minimumtac=",mintac,"\n")}
#
# END: READ IN THE OPTION FILE
##############################################################################
##############################################################################
# READ IN THE (AGE) DATA FILE
#
# This file contains age structured biological information about the stock
# to be used in the simulation. The code below is in following steps:
# 1. Check if the file exists, if not STOP the program
# 2. Read in the whole file. Check for consistency in the number of age
# age groups in the data file and the control file. STOP if there is
# inconsistencey
# 3. Assign individual values to local parameters
# 1. Check for existence
if (!(file.exists(datfile))) {
stop(paste("\nInitialisation file",datfile,"doesn't exist\nUnable to continue\n",sep=" "))}
# 2. Read in the whole file and check for inconsistency
indata<-matrix(scan(datfile,0,quiet=TRUE),ncol=17,byrow=TRUE) # EINAR - replaced ncol=18 with ncol=17
numage<-length(indata[,1])
if ((maxage-minage+1)!=numage) {
stop(paste("\nThe age class data supplied in the initialisation data file (",datfile,")\ndoes not match the age class range (",minage,"-",maxage,") in the options file (",optionsfile,")\n",sep=""))
}
# 3. Read in each vector file
initage<-indata[,1] #age class
initpopint<-indata[,2] #population (000s)
initpopcv<-indata[,3] #population cv
initspwnweightint<-indata[,4] #spawning weight (kg)
initspwnweightcv <-indata[,5] #spawning weight cv
initcatweightint<-indata[,6] #catch weight (kg)
initcatweightcv<-indata[,7] #catch weight cb
initmatureint<-indata[,8] #maturity
initmaturecv<-indata[,9] #maturity cv
initfint<-indata[,10] #fishing mortality
initfcv<-indata[,11] #fishing mortality cv
initfprop<-indata[,12] #proportion of fishing mort. before spawning
initdisfint<-indata[,13] #discard mortality
initdisfcv<-indata[,14] #discard mortality cv
initmint<-indata[,15] #natural mortality
initmcv<-indata[,16] #natural mortality cv
initmprop<-indata[,17] #proportio of natural mort. before spawning
#
# END: READ IN THE DATA FILE
##############################################################################
##############################################################################
# READ IN HISTORICAL RECRUITMENT, SSB, CATCH & FISHING MORTALITY
#
# 1. Read in the recruit file. The file contains year, recruitment at a
# at a specific age and spawning stock biomass. The setup of the
# external file is thus in the classical ICES summary table format
#
# 2. Read in the historical catch and and F file
# 1. Historical recruit file
input<-scan(recfile,list(0,0,0),quiet=TRUE)
recyear<-input[[1]] #recruitment data year vector
recr<-input[[2]] #recruitment data recruitment numbers vector
rssbin<-input[[3]] #recruitment data ssb vector
SS <- data.frame(year=recyear,recr=recr,ssb<-rssbin) # Einar check
# 2. Historical catch file
input<-scan(catchfile,list(0,0,0),quiet=TRUE)
catyear<-input[[1]]
ccatch<-input[[2]]
cfbar<-input[[3]]
#
# END: READ IN HISTORICAL RECRUITMENT, SSB, CATCH & FISHING MORTALITY
##############################################################################
##############################################################################
# START of INITIALISATION
#
#
# create a vector of output filenames. The number of output files
# will be resolution+1
FLRfiles<-paste("FLRObjects_",runref,"_",1:(resolution+1),".dat",sep="")
# The varmin in the options file specifies either the lowest fishing mortality
# (variable=1) or the lowest catch (variable=0) to be tested. In the running
# of the model this value gets assign to the local variable ftac.
ftac<-varmin # start with the min. value of variable (F or TAC)
ftacmult<-1 # used in controlling the fisheries via managment module
oldftacmult<-1 # initialisation
endyear<-startyear+years-1
lenrecyr<-length(recyear) # length of recruitment data set
if (histrec==1) {maxrecyr<-max(recyear)} #
if(histcatch==1){
maxcatyr<-max(catyear) #last year of catch data
lencatyr<-length(catyear)} #number of years of catch data}
else {
maxcatyr<-0} #set maximum year to zero
### Initialise arrays for storage ##########################################
# Initialise a number of arrays to STORE the various parameters of
# interest. Each array will hold all values i.e. over all age groups,
# simulations years and iterations for each resolution. Initial values are
# all zero. The arrays are designed to look like FLQuant objects (prior to
# versrion 2), hence the age, year, unit, season and area dimension names
# and some redundancy (with unit and season). The area dimension is being
# used in this context as the iteration dimension
# Note that with FLR version 2, the iterations should be moved and stored
# in the iter slot (6th dimension). This would only need minor modification
# in the run.r script. However, the FPRESS plotting routines would also
# need to be updated, hence will wait a while
# Note that newbiomassarray was put in by EINAR. The intent here is to
# store the B4+ biomass for iCod. This actually also requires a new input
# column in the datafile.
arr_dim<-c(numage,years,1,1,iterations)
arr_dimnames=list(age=minage:maxage,
year=startyear:endyear,
unit='unique',
season='all',
area=1:iterations)
catchweightarray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
catcharray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
newspawnmassarray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
newbiomassarray<-array(0,dim=arr_dim,dimnames=arr_dimnames) #EINAR - added this for iCod
yearendstockarray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
poparray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
fishingmortalityarray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
fishingmortalitymultarray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
farray<-array(0,dim=arr_dim,dimnames=arr_dimnames)
# the following arrays do not include age-disaggregated data and so
# do not expand this dimension
arr_dim_noage<-c(1,years,1,1,iterations)
arr_dimnames_noage<-list(age=0,
year=startyear:endyear,
unit='unique',
season='all',
area=1:iterations)
recruitsarray<-array(0,dim=arr_dim_noage,dimnames=arr_dimnames_noage)
fbararray<-array(0,dim=arr_dim_noage,dimnames=arr_dimnames_noage)
ftacarray<-array(0,dim=arr_dim_noage,dimnames=arr_dimnames_noage)
ftacmultarray<-array(0,dim=arr_dim_noage,dimnames=arr_dimnames_noage)
hcrsubarray<-array(0,dim=arr_dim_noage,dimnames=arr_dimnames_noage)
### Adjust spawning and recruit data ########################################
#
# Create three new vectors to hold the recruitment data offset so that
# the recruitment value corresponds to the spawning year rather than the
# year it entered the fishery. Depending on the value of minage, this is
# likely to require recruitment data for years that are after the start
# date of the projection (e.g. if the simulation starts in 2000 and
# minage = 2, then recruitment data will be required for at least 1998,
# 1999, 2000 and 2001).
#
# NOTE: The variables generated here are actually not used in the
# subsequent run. Thus this block is effectively REDUNDANT
diffxs<-max(recyear)-startyear+1 # year between start of simulation
# and final year w. availabe recruitment
spwnyr<-seq(1,1,length=lenrecyr-max(minage,diffxs))
spwnssb<-spwnyr
spwnrec<-spwnyr
for(i in 1:(lenrecyr-max(minage,diffxs))){
spwnyr[i]<-recyear[i]
spwnssb[i]<-rssbin[i]
spwnrec[i]<-recr[i+minage]
}
SSadjust1 <- data.frame(year=spwnyr,recr=spwnrec,ssb<-spwnssb) #EINAR check
### The SSB-R dataset to be used estimating SSB-R parameters ################
#
# Same procedure as above, except the vectors may be of a different
# length than those above depending on the value of datarecyear.
# This parameter gives the
# maximum year in the recfile that data is to be used from in the
# recruitment functions. This allows fictional future recruitment to be
# entered into the recruits data file, but only the real historical data to
# be used in fitting parametric functions or bootstrapping. If
# datarecyear > max(recyear) (the oldest year in the recruits data file),
# then all data in the file is used in the recruitment functions. The
# vectors spwnfnyr,spwnfnssb, spwnfnrec are used in the recruitment functions
# calculate the number of years between the last year of recruitment data
# and the final year that is to be used in the simulations
diffxt<-max(recyear)-datarecyear+1
spwnfnyr<-seq(1,1,length=lenrecyr-max(minage,diffxt))
spwnfnssb<-spwnfnyr
spwnfnrec<-spwnfnyr
for(i in 1:(lenrecyr-max(minage,diffxt))){
spwnfnyr[i]<-recyear[i]
spwnfnssb[i]<-rssbin[i]
spwnfnrec[i]<-recr[i+minage]
}
SSadjust2 <- data.frame(year=spwnfnyr,recr=spwnfnrec,ssb<-spwnfnssb) #Einar check
#
# END: INITIALISATION
##############################################################################################
################################################################################
# START OF MAIN BODY #
################################################################################
print(stockname)
print(Sys.time())
startdate<-as.POSIXlt(Sys.time(),Sys.timezone()) # Store start time
###############################################################################
# Loop over range of of F or TAC values, res is the counter variable #
###############################################################################
for(res in 0:resolution) {
if (debug==1) {cat("res=",res,"\n")}
#calculate the appropriate FTAC value for this resolution
if (resolution>0) {
ftac<-((varmax-varmin)*res/resolution)+varmin
if(ftac==0){ftac<-0.0001} #ensure we don't proceed with a 0 value
ftac2 <- ftac
# calculate the percentage of the simulation that is
# complete for reporting purpose
percent<-round((res/(resolution+1))*100,digits=2)
if (variable==1) {
cat("fval = ",ftac," ",percent,"% complete","\n",sep="")
} else {
cat("tacval = ",ftac," ",percent,"% complete","\n",sep="")
}
}
flush.console() #update the console to ouput above messages
##############################################################################
# Main iterative loop, iter is the counter variable #
##############################################################################
for(iter in 1:iterations){
#if the resolution is zero then the simulation is running for a single
#TAC or F value. In this case report progress on basis of the number of iterations
#completed. Issue report every 10%
if (resolution == 0) {
percent <- 100*iter/iterations
if ((percent%%10) == 0) {
cat(percent,"% complete","\n",sep="")
}
}
flush.console()
if (debug==1) {cat("iter=",iter,"\n")}
#############################################################################
# Year loop, year is the counter variable #
#############################################################################
for(year in 1:years){
if (debug==1) {cat("year=",year,"\n")}
simyear<-year+startyear-1 #calculate the current simulation year
if (debug==1) {cat("A&M..")}
############################################################################
# ASSESSMENT, MANAGEMENT AND F/TAC to be applied
#
# 1a. Assessment: The function (see assessment.r) calculates new
# applying user specified noise (cv) and bias to existing values
# of bio, ssb and/or fbar. This effectively emulates uncertainty
# in stock assessment.
# Output: New bio, ssb and f
# 1b. Management: The function (see managment.r) calculates the
# ftac multiplier (ftacmult) based on a user specified HCR given
# an uncertain assessment
# Output: New ftacmult
# 2 F or TAC to be applied this year (if not historical values)
# 2.a) If not historical values
# 2.b) If historical values
# Output: New Fbar or TAC (ftac)
#
# Note that the estimates are correspond to the following years
# relative to the assessment year y
# SSBy-1
# BIOy-1
# F?
# 1. Assessment model & Managment ##########################################
#
if((year-assessdelay)>0){
assres<-assessment(bio,ssb,
fbar,
ssbassessbias,
ssbassesscv,
fbarassessbias,
fbarassesscv)
bio <- assres[1]
ssb<-assres[2]
fbar<-assres[3]
man_ret<-management(bio,ssb,
fbar,
ftac,
year,
years,
assessdelay,
hcrperiod,
hcrrule,
hcrchange,
ssbhcr,
oldftacmult,
fbarhcr,
ssbtest1,
ssbtest2,
mintac)
ftacmult<-man_ret[1]
hcrsub<-man_ret[2] #...(not currently used)
} else {
ftacmult<-1 # no change is effected
hcrsub<-0
}
ftacmultarray[1,year,1,1,iter]<-ftacmult # permanent storage
hcrsubarray[1,year,1,1,iter]<-hcrsub # permanent storage
### 2a. F or TAC to be applied this year (if not historical values) ########
#
# Here the ftac parameter gets modified by the ftacmult parameter which
# came out from the management function. Note that the origftac is a
# placeholder for the ftac that is set in place at the beginning of
# each resolution-loop.
# The outputted ftac value has the following meaning depending on the
# setting of the parameter 'variable':
# variable=1: Fbar in F-based simulation
# variable=0: TAC in TAC based simulation
# variable=2: harvest rate in TAC based simulation
#
# output: ftac
#
origftac<-ftac # save the original value
# Adjust the ftac according the the management rule (ftacmult)
if ((variable==0) & (ftacmult*ftac < mintac)) ftac=mintac
else {ftac<-ftacmult*ftac}
if (debug==1) {cat("Init..")}
tmp <- ftac # EINAR: Temporary placeholder
# HARDWIRE - iCod HCR ####################################################
if(variable==2 & year>1) {
if (hcrchange == 0) {
ftac<-max(ftac*bio,mintac)
} else {
ftac<-max((ftac*bio+ftacarray[1,year-1,1,1,iter])/2,mintac)}}
if (year> 1) myText <- paste(year,"origftac=",origftac,"ftacmult=",round(ftacmult,3),"ftac1=",round(tmp,3),"ftac2=",round(ftac,0),"hr=",round(ftac/bio,2))
# 2b. F or TAC to be applied this year (if historical values) ############################
#
# Historical F og Yield values are in most cases needed , at least for the first year(s).
# One may also be intersted in running "hindcasting" from some historical starting year.
# Here the appropriate yield or fbar from the historical dataset are located and a user
# defined (histcatchcv) allocated ???? if historical values are used for catch then the
# bias and cv ssociated with F/TAC will be set to 1 and 0 respectively????
#
# Output: New Fbar or TAC (ftac)
#
newvarcv<-varcv
newvarbias<-varbias
if ((histcatch==1) & (simyear<=maxcatyr)) {
#identify the relevant year
index<-0
for (i in 1:lencatyr) {
if(catyear[i]==simyear){index<-i}}
#now perform the random draw
#F controlled fishery
if (variable==1) {
ftac<--1
#randomise
while(ftac<0){ftac<-rnorm(1,cfbar[index],cfbar[index]*histcatchcv)}
#reset cv, bias parameters
newvarcv<-0
newvarbias<-1
}
#TAC controlled fishery
if (variable %in% c(0,2)) { # EINAR - iCod adaptation
ftac<--1
#randomise
while(ftac<0){ftac<-rnorm(1,ccatch[index],ccatch[index]*histcatchcv)}
#reset cv, bias parameters
newvarcv<-0
newvarbias<-1
}
}
#update the ftac array
ftacarray[1,year,1,1,iter]<-ftac
#print(ftac)
if (debug==1) {cat("Fmult..")}
#
# END: ASSESSMENT, MANAGEMENT AND F/TAC to be applied
############################################################################
# UPDATING AGE-BASED PARAMETERS
# 1. Initialize values
# 2. Update the Nay matrix
# 2.a) Update the Nay matrix
# 2.b) Update the recruits (overrides the youngest age in step 1.a)
# 3. Update cWay, sWay, MATay, cSay, dSay and May
# 4. Update Fay
# 4.a) F-based fisheries
# 4.b) TAC-based fisheries
### 1. Initialise age vectors for start of year ############################
#
# If in the first year of the simulation, used the populations as specified
# in the input datafile. If past the first year of the simulation, use the
# population simulated at the end of the previous year.
# Thus, popcv (which conventionally would be labelled as Nay contains the
# true population numbers in the assessment year (year y). Note that this
# is different from the ssb, bio and f above, to which cv and assessment
# bias have been applied to
#
# There may be some redunance in usage here, revisit that stuff later
age<-initage #age class EINAR - looks redundant
if (year>1) {popint<-newyearplus1} else {popint<-initpopint}
popcv<-initpopcv #population cv
spwnweightint<-initspwnweightint #spawning weight
spwnweightcv <-initspwnweightcv #spawning weight cv
catweightint<-initcatweightint #catch weight
catweightcv<-initcatweightcv #catch weight cv
matureint<-initmatureint #maturity
maturecv<-initmaturecv #maturity cv
fint<-initfint #fishing mortality
fcv<-initfcv #fishing mortality cv
fprop<-initfprop #prop. of fishing mor. prior to spawning
disfint<-initdisfint #discard mortality
disfcv<-initdisfcv #discard mortality cv
mint<-initmint #natural mortality
mcv<-initmcv #natural mortality cv
mprop<-initmprop #prop. of natural mort. before spawning
### 2a. Update the Nay matrix (all except recruits) #######################
# Pop is assigned the Nay values. These numbers could either come from the
# end of the previous year (former year-loop or if this is the first year,
# from the input file.
# Adding noise / stochasticity to account for environmental/operational
# variability
# Initial population (randomised in first year of projection only).
# If year>1 the randomization has taken place at the end of the previous
# year loop.
#
pop<-popint
popintsd<-(popcv*popint)
if(year==1){
pop <- -popint
#keep drawing until a positive result is obtained
while(min(pop)<0){
for(i in 1:numage){pop[i]<-rnorm(1,popint[i],popintsd[i])}
}
}
### 2b. Updating the Nay matrix - this year recruits #######################
# Recruitments will only be supplied here if minage>0. When minage=0
# the fish enter the fishery immediately after spawning. These cases
# are dealt with in a later code (once SSBy has been calculated)
#
# Here (when minage>0) incoming recruits can come from one of two sources:
# 1. The historical recruit file
# Since these values arise from ssb values prior to the execution
# period of the program and more importantly because the estimates
# may be based on actual measurements, the recruitment values are
# stored in the historical recruitment file. An account must be taken
# of the age of recruitment in order to correctly locate the right
# recruit relative to the simulation year. Once locatated an
# historical recruitment cv value is also applied.
# As an example, if minage=3 and the model starts execution in 2000,
# the first year that the model can calculate recruitment for is 2003.
# Thus historical values are required for 2000, 2001 and 2002
# 2. Recruits generated internally from a stock-recruit function. In that
# case the values have already been generated in previous years inside
# the loop and stored in the recruitsarray
#
# Output:
#
if (debug==1) {cat("Rec after spawn..")}
if (minage>0) {
# 1. Historical recruit values, stored in the input files
if (simyear-minage <= maxrecyr) {
for (i in 1:lenrecyr) {if (recyear[i]==simyear-minage) {index<-i}}
frec<--1
while(frec<0){frec<-rnorm(1,recr[index],recr[index]*histreccv)}
recruitsarray[1,year,1,1,iter]<-frec
}
# 2. Recruit values already generated internally in prior year loops
else {frec<-as.numeric(recruitsarray[1,(year-minage),1,1,iter])}
fishrecruited<-frec
} # end if(minage>0)
if(minage>0){pop[1]<-fishrecruited}
### 3. Updating other age parameters (sWay, cWay, MATay, cSay, dSay, May #################
# The objective here is to generate numbers that are equivalent to the real numbers
# "observed" numbers (remember, we are continuously predicting in this simulation
# framework. These new numbers are then used to calculate the actual stock size, as one
# would have as measurements it in the beginning of the next year loop (end of this years
# loop). What is done here is just to supply new random numbers to the original input,
# not taking into account correlations among ages. Note, the randomization process is
# fully independent across age groups. Effectively this may mean that the sum or the
# agerage of this white noise will be zero. E.g. since SSB is the sumproduct numbers,
# weight and maturity, the white noise within each of these vectors may cancel each
# other out when it comes to calcuating SSB. Remedy: Call Höski :-)
# Additional thing - should store all these values and export
# - the stuff below calls for a seperate function
# - like myWhile(XXXint,XXXvn)
#Spawning weight
spwnweight<--spwnweightint
spwnweightsd<-spwnweightcv*spwnweightint
while(min(spwnweight)<0){
for(i in 1:numage){
spwnweight[i]<-rnorm(1,spwnweightint[i],spwnweightsd[i])}}
#Catch weight
catweight<--catweightint
catweightsd<-catweightcv*catweightint
while(min(catweight)<0){
for(i in 1:numage){
catweight[i]<-rnorm(1,catweightint[i],catweightsd[i])}}
#Maturity
mature<--matureint
maturesd<-maturecv*matureint
while(min(mature)<0){
for(i in 1:numage){
mature[i]<-min(rnorm(1,matureint[i],maturesd[i]),1)}}
#Fishing mortality - more properly referred to as selection pattern
f<--fint
fsd<-fcv*fint
while(min(f)<0){
for(i in 1:numage){
f[i]<-rnorm(1,fint[i],fsd[i])}}
#Discard mortality
disf<--disfint
disfsd<-disfcv*disfint
while(min(disf)<0){
for(i in 1:numage){
disf[i]<-rnorm(1,disfint[i],disfsd[i])}} # EINAR - BUG, removed fmult
#Natural mortality
m<--mint
msd<-mcv*mint
while(min(m)<0){
for(i in 1:numage){
m[i]<-rnorm(1,mint[i],msd[i])}}
# New bar - calculated by taking the mean of f over the age groups
# defined by the parameters fbarmin and fbarmax
newfbar<-mean(f[(fbarmin-minage+1):(fbarmax-minage+1)])
if (debug==1) {cat("ftac..")}
### 4a. Calculate Fay if F-based simulation #############################################
#
# Apply F multiplier to F vectors (only required if simulation is based on f). The
# specified level of F given by ftac is altered by possible bias (overfishing) given
# by varbias and then randomised by drawing from a Normal distribution. The randomised
# value of F (freq) is the required fishing mortality for this year. A multiplier, fmult,
# is calculated from newfbar so that when multiplied by the current F vector it will give
# the required value of freq. Finally, the F vector f is multiplied by this multiplier
# along with the discard F vector disf. In reality, discard F may not change in the same
# way as F when catches are reduced.
#
if(variable==1){
fnewval <- ftac * newvarbias # EINAR- replaced newfbias with newvarbias
nf<--1
while(nf<0){nf<-rnorm(1,fnewval,fnewval*newvarcv)} #randomise fnewval
freq<-nf #this is the required f vector
fmult<-freq/newfbar #calculate f multiplier
f<-fmult*f #apply multiplier to f
disf<-disf*fmult #apply multiplier to discard f
}
### 4a. Calculate Fay if TAC-based simulation ###########################################
#
# If using a TAC, calculate the appropriate F multiplier. The specifie level of TAC given
# by FTAC is altered by possible bias (overfishing) given by varbias and then randomised
# (this occurs each year) by drawing from a Normal distribution. The randomised value of
# TAC given by tac is the required catch for this year. Function ffactor is then used to
# calculate an F multiplier, ffac, that when applied to the current F vector f, will
# produce the required catch equal to the TAC (so that fishing mortality in the simulation
# is always actually calculated using F even when the virtual fishery uses TAC). Finally,
#the F vector f is multiplied by this multiplier along with the discard F vector disf.
if(variable %in% c(0,2)){ # EINAR - iCod adaptation
tacval<-ftac*newvarbias #apply bias
ntac<--1
while(ntac<0){ntac<-rnorm(1,tacval,tacval*newvarcv)} #randomise
tac<-ntac #the requied tac vector
ffac<-ffactor(tac,pop,f,disf,catweight,m,numage) #calculate f vector
f<-f*ffac #apply multipler to f
disf<-disf*ffac #apply multiplier to discard f
fmult<-ffac
}
#Update fishing mortality and fishing mortality multiplier arrays
fishingmortalityarray[,year,1,1,iter]<-f
fishingmortalitymultarray[,year,1,1,iter]<-fmult
#
# END: UPDATING AGE-BASED PARAMETERS
############################################################################
############################################################################
# POPULATION DYNAMICS
#
# The steps are:
# 1. Calculate population size and hence ssb at spawning time
# 2. Estimates recruitment based on a user specified recruitment model.
# 2.1 Obtain new recruits from a model (may or may not be used in the
# next step.
# 2.1 Decistion algorithm that speciefies
# a) Should recruitement be entered into the population
# immediately (only needed if fishing occurs on age 0)
# b) Should one obtain recruitment be obtained from historical
# data or should on simulate recruitment based on one
# of the user specifed recruitment above (alrady done in
# step 2.1)
# 3. Calculate the catch annual catch
# 4. Calculate the population size at the end of the simulation year
### 1. SSB - based on population numbers at spawning time ##################
if (debug==1) {cat("Dyn at spawn..")}
fbef<-(f+disf)*fprop # Fishing mortality prior to spawning
faft<-(f+disf)*(1-fprop) # Fishing mortality after spawning
mbef<-m*mprop # Natural mortality prior to spawning
maft<-m-mbef # Natural mortality after spawning
if (debug==1) {cat("fbef,mbef =",fbef,mbef,"..")}
newstock<-(pop*exp(-mbef-fbef)) # Population numbers at spawning time
newspawnpop<-(newstock*mature) # Mature population number at spawn. time
newspawnmass<-newspawnpop*spwnweight # Mature biomass at spawn. time
newspawnmassarray[,year,1,1,iter]<-newspawnmass # store biomass
ssb<-sum(newspawnmass)
### 2. RECRUITMENT #########################################################
if (debug==1) {cat("Recruit..")}
# 2.1 Recruitment model
recruitsnew <- switch(rec+1,
recruit0(recmean),
recruit1(recmean,reccv),
recruit2(spwnfnyr,spwnfnrec),
recruit3(ssb,spwnfnyr,spwnfnrec,spwnfnssb,ssbcut),
recruit4(ssb,spwnfnrec,spwnfnssb),
recruit5(ssb,spwnfnrec,spwnfnssb),
recruit6(ssb,ssbcut,recmean,seggrad),
recruit7(ssb,ssbcut,recmean,seggrad,reccv),
recruit8(ssb,reccv))
# 2.2. Entering recruitment or store for later retrieval
#
# If the minimum age is zero then the fish enter the fishery from age 0
if(minage==0){
if(simyear > maxrecyr){
# If the current simulation year is greater than the oldest year in
# the recruitment data file, then set
# the recruitment (frec) = recruitsnew
frec<-recruitsnew
} else {
#if the current simulation year is contained in the recruitment
# data file then the expected recruitment is set to the
# corresponding value given in the file
for(i in 1:lenrecyr){
if(recyear[i]==simyear){num<-i}
}
#randomise the number of recruits
frec<--1
while(frec<0){frec<-rnorm(1,recr[num],recr[num]*histreccv)}
}
#assign the new recruits to the first population vector element
newstock[1]<-frec
yearrecruits<-frec
#write data to local recruitsarray
recruitsarray[1,year,1,1,iter]<-yearrecruits
} else {
# the value recruitsnew applies to a future year write this into
# a future slot in the recruitsarray (provided there is one -
# no point in writing beyond end of simulation period)
if ((year+minage) <= years) {
recruitsarray[1,year+minage,1,1,iter]<-recruitsnew
}
}
if (debug==1) {cat("Removing Catch..")}
### 3. TAKE THE CATCH ######################################################
fbarz<-c()
for(i in (fbarmin-minage+1):(fbarmax-minage+1)){
fbarz<-c(fbarz,f[i])
}
fbar<-mean(fbarz)
#update local array
fbararray[1,year,1,1,iter]<-fbar
#calculate the catch in number
catch<-(pop*f/(f+disf+m+0.00001)*(1-exp(-f-disf-m)))
#calculate the catch weight
catchwght<-catch*catweight
#update local array
catchweightarray[,year,1,1,iter]<-catchwght
if (debug==1) {cat("Dyn After Spawn\n")}
if(year>1 & debug==2) {print(paste(myText,"fbar=",round(fbar,2)))}
### 4. POPULATION SIZE IN THE END OF THE YEAR ##############################
# calculate the year end population by applying the fishing and
# natural mortality that occurs after spawning
yearendstock<-newstock*exp(-maft-faft)
#update local arrays
yearendstockarray[,year,1,1,iter]<-yearendstock
poparray[,year,1,1,iter]<-pop
catcharray[,year,1,1,iter]<-catch
farray[,year,1,1,iter]<-f
#plus group
yearplus1<-0
for(i in 1:(numage-2)){
yearplus1<-c(yearplus1,yearendstock[i])
}
#add on the plus group to the oldest age (last spaces)
yearplus1<-c(yearplus1, yearendstock[numage-1]+yearendstock[numage])
#EINAR: Calculate the new bio
newbio <- (yearplus1*c(0,0,0,rep(1,11)))
newbiomass <- newbio*catweight
#update the local array with the new biomass for this iteration/year
newbiomassarray[,year,1,1,iter]<-newbiomass
bio <- sum(newbiomass)
#save vals in temp vars for initialisation of population in next year
newyearplus1<-yearplus1
# Note that the output has original input variances
# NOT new variances calculated from the data
ftac<-origftac
}
###########################################################
#end of year loop##########################################
###########################################################
}
##########################################################
#end of main iterative loop###############################
##########################################################
#Save data
#Create FLQuant objects
#FLQuant() is the creator function for FLQuant objects
CatchWeight<-FLQuant(catchweightarray,units='Kgs')
NewSpawnMass<-FLQuant(newspawnmassarray,units='Kgs')
FBar<-FLQuant(fbararray,units='F')
Recruits<-FLQuant(recruitsarray,units='000s')
YearEndStock<-FLQuant(yearendstockarray,units='Kgs')
Pop<-FLQuant(poparray,units='Kgs')
Catch<-FLQuant(catcharray,units='Kgs')
FishingMortality<-FLQuant(fishingmortalityarray,units='F')
FishingMortalityMult<-FLQuant(fishingmortalitymultarray,units='F')
F<-FLQuant(farray,units='F')
FTac<-FLQuant(ftacarray,units='ftac')
FTacMult<-FLQuant(ftacmultarray,units='mult')
#save specified objects
if (substr(outdata.path,nchar(outdata.path),nchar(outdata.path))=="\\") {
outfile<-paste(outdata.path,FLRfiles[res+1],sep="") } else
{
outfile<-paste(outdata.path,"\\",FLRfiles[res+1],sep="")
}
save(list=savedobjs,file=outfile)
}
############################################################
#end of resolution loop#####################################
############################################################
#some end logging
print(Sys.time())
#calculate the end time and the difference with the starttime and
#print duration of run to screen
enddate<-as.POSIXlt(Sys.time(),Sys.timezone())
diff<-difftime(enddate,startdate,tz="GMT",units="secs")
cat("Execution Time:",diff[[1]]," secs", "\n")
#restore working directory
setwd(currwd)
# EOF
}
As is always the case, the iCod dynamics and HCR has some specifics that may not be available in any off-the-shelf frameworks. So far the following have been identified, with regards to the generic FPRESS:
The key controlling variables in regards to the fisheries in FPRESS are:
if(variable==1){
#apply bias
fnewval <- ftac * newvarbias # EINAR- replaced newfbias with newvarbias
#randomise fnewval
nf<--1
while(nf<0){nf<-rnorm(1,fnewval,fnewval*newvarcv)}
#this is the required f vector
freq<-nf
#calculate f multiplier
fmult<-freq/newfbar
#apply multiplier to f
f<-fmult*f
#apply multiplier to discard f
disf<-disf*fmult
}
if(variable %in% c(0,2)){ # EINAR - iCod adaptation
#apply bias
tacval<-ftac*newvarbias #=1 no bias, <1 under, >1 over
#randomise
ntac<--1
while(ntac<0){ntac<-rnorm(1,tacval,tacval*newvarcv)}
#this is the requied tac vector
tac<-ntac
#calculate f vector
ffac<-ffactor(tac,pop,f,disf,catweight,m,numage)
#apply multipler to f
f<-f*ffac
#apply multiplier to discard f
disf<-disf*ffac
fmult<-ffac
}