As a part of the evaluation of the iCod HCR it was decided, in addition to running an in-house framework (ADMODEL builder), to make independent evaluation using a framework that had already been used within the ICES forum to evaluate HCR for other stocks. The FPRESS is currently such a candidate. It's suitability with regards to evaluation of iCod is currently (fall 2009) under investigation.
FPRESS was conceived by Cirian Kelly and Andrew Campbell at Marine Institute, Galway. The generic source code was obtained from Andrew Campell, 12.8.2009.
It should be noted that this generic code is relatively old and has not been maintained by the original team. Or as stated by Andrew Campell in e-mail correpondence: ”This version of FPRESS is rather old at this stage. Despite my best intentions, it has been very difficult to find the means and the time to maintain a generic version of the model for distribution. Work has tended to focus on developing it for specific applications such as the NEA Mackerel LTM (as you mentioned below) and the Western Horse Mackerel Management Plan development. It has also been used for Irish Sea Cod and Sardine work.” This means that numerous “bugs” that may have been discovered in later works have not been implemented/corrected in the generic code. Given this the intent is to document separately:
The following is a list of changes made to the generic code.
The following changes are absolutely necessary (has been confirmed by Andrew), otherwise the script will halt:
The changed script (with some additional “out-commenting” of redundant code) can be found here.
Additional error codes not yet in the changed script above (will do once confirmed with author):
The following changes are needed if stochastic hockey-stick is used as a recruitment model
The changed script can be found here.
In generic FPRESS code provides ample commentaries in individual places. What is however lacking, at least to one poor sod that is trying to understand the generic codes are broader commentaries related to the block of codes at each stage. What is provided here is a major revamping of the structure of the code and commentaries to the original run.r code, with (hopefully) just a minimal changes to the actual executable part [sorry, did also some iCod specifics in the generic code - will to a retro on that soon].
#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
}