Commit 500f22e2 authored by Berges, Benoit's avatar Berges, Benoit
Browse files

update MSE scripts

parent 3eb93d46
......@@ -19,7 +19,7 @@ mkdir(file.path("model","packages"))
data.source <- file.path(".", "data")
model.save <- file.path(".",'model')
runName <- '01_baseline'
runName <- 'NSAS_sanoba'
SMSkeyRuns <- 2019
source(file.path('functions','config_baseline_sf.R'))
......@@ -31,6 +31,9 @@ source(file.path('functions','config_baseline_mf.R'))
load(file.path(data.source,paste0("data_sms",SMSkeyRuns,"_sf.RData"))) # NSH, NSH.tun
addM <- 0.06
NSH@m <- NSH@m + addM
# create control
pg <- range(NSH)[2]
......@@ -43,28 +46,31 @@ NSH@harvest <- NSH.sam@harvest[,ac(range(NSH)["minyear"]:range(NSH)["maxyear"])]
save(NSH, NSH.tun, NSH.ctrl, NSH.sam,
file=file.path(model.save,
paste0("results_",runName,"_sf.RData")))
paste0(runName,"_assessment_results_sf.RData")))
### ============================================================================
### multi fleet
### ============================================================================
#load(file.path(data.source,'results_mf_2020.Rdata'))
#NSH3f.sam.stk0 <- NSH3f.sam
load(file.path(data.source,paste0("data_sms",SMSkeyRuns,"_mf.RData")))
NSH3.ctrl <- config_baseline_mf(NSHs3,NSH.tun)
addM <- 0.06
NSHs3$residual@m <- NSHs3$residual@m + addM
NSHs3$sum@m <- NSHs3$sum@m + addM
load(file.path(data.source,'results_mf_2020.Rdata'))
NSH3f.sam.stk0 <- NSH3f.sam
NSH3.ctrl <- config_baseline_mf(NSHs3,NSH.tun)
# run model
NSH3f.sam <- FLSAM(NSHs3,
NSH.tun,
NSH3.ctrl,
starting.values = NSH3f.sam.stk0)
NSH3.ctrl)#,starting.values = NSH3f.sam.stk0
save(NSHs3,
NSH.tun,
NSH3.ctrl,
NSH3f.sam,
file=file.path(model.save,
paste0("results_",runName,"_mf.RData")))
paste0(runName,"_assessment_results_mf.RData")))
#-------------------------------------------------------------------------------
# WKNSMSE
#
# Author: Benoit Berges
# WMR, The Netherland
# email: benoit.berges@wur.nl
#
# MSE of North Sea Herring
#
# Date: 2018/11/18
#
# Build for R3.5.1, 64bits
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 1) load packages
# setup paths
......@@ -24,22 +10,14 @@ library(FLSAM)
library(FLEDA)
library(FLFleet)
# define path to directory
#path <- "D:/Work/Herring MSE/NSAS/"
#path <- "D:/git/wk_WKNSMSE_her.27.3a47d/R/"
#path <- "F:/WKNSMSE/wk_WKNSMSE_her.27.3a47d/R"
#path <- 'E:/wk_WKNSMSE_her.27.3a47d/R'
path <- 'E:/git/wk_WKNSMSE_her.27.3a47d/R'
#path <- '/home/berge057/ICES/wk_WKNSMSE_her.27.3a47d/R/'
assessment_name <- "NSAS_WKNSMSE2018"
try(setwd(path),silent=TRUE)
# paths to different subfolders
dataPath <- file.path(".","data/")
outPath <- file.path(".","results/")
scriptPath <- file.path(".","side_scripts/")
modelPath <- file.path(".","model/")
functionPath <- file.path(".","functions/")
runName <- "NSAS_sanoba"
SMSkeyRuns <- 2019
# loading function
source(file.path(functionPath,"randBlocks.R"))
source(file.path(functionPath,"randNums.R"))
......@@ -59,15 +37,13 @@ source(file.path(functionPath,"make.arma.resid.lst.R"))
# complicated to implement
#-------------------------------------------------------------------------------
#- Load single fleet and multi fleet assessment objects - using assessments without the LAI
#load(file.path(outPath,paste0(assessment_name,'_mf.Rdata')))
#load(file.path(outPath,paste0(assessment_name,'_sf.Rdata')))
load(file.path(outPath,paste0(assessment_name,'_mf_noLAI.Rdata')))
load(file.path(outPath,paste0(assessment_name,'_sf_noLAI.Rdata')))
load(file.path(modelPath,
paste0(runName,"_assessment_results_sf.RData")))
load(file.path(modelPath,
paste0(runName,"_assessment_results_mf.RData")))
# parameters
n.retro.years <- 7 # Number of years for which to run the retrospective
nFutureyrs <- 40 + 3
nFutureyrs <- 20 + 3
histMinYr <- dims(NSH)$minyear
histMaxYr <- dims(NSH)$maxyear
yearCurrent <- histMinYr:histMaxYr # vector the years
......@@ -75,21 +51,25 @@ futureMaxYr <- histMaxYr + nFutureyrs
histPeriod <- ac(histMinYr:histMaxYr)
projPeriod <- ac((histMaxYr+1):futureMaxYr)
fullPeriod <- c(histPeriod,projPeriod)
recrPeriod <- ac(2007:2017)
selPeriod <- ac(2007:2017)
fecYears <- ac(2007:2017)
nits <- 2000 # number of random samples
recrPeriod <- ac(2011:2021)
selPeriod <- ac(2010:2020)
fecYears <- ac(2010:2020)
nits <- 10 # number of random samples
load(file.path(dataPath,paste0("data_sms",SMSkeyRuns,"_sf.RData")))
# reading the raw M and applying plus group
#raw_M <- read.csv(file.path(dataPath,"Smoothed_span50_M_NotExtrapolated_NSASSMS2016.csv"),header=TRUE)
raw_M <- read.csv(file.path(dataPath,"Raw_NotExtrapolated_NSAS_SMS2016.csv"),header=TRUE)
raw_M <- read.csv(file.path(dataPath,
paste0("M_NSAS_smoothedSpan50_notExtrapolated_sms",ac(SMSkeyRuns),".csv")),
header=TRUE)
colnames(raw_M) <- sub("X","",colnames(raw_M))
rownames(raw_M) <- raw_M[,1]
raw_M <- raw_M[,-1]# Trim off first column as it contains 'ages'
raw_M <- cbind(replicate(as.numeric(colnames(raw_M)[1])-histMinYr,raw_M[,1]), raw_M)
raw_M <- cbind(raw_M,raw_M[,dim(raw_M)[2]])
colnames(raw_M) <- histMinYr:histMaxYr
raw_M <- raw_M[1:9,] + 0.11 # trim age 9 and add 0.11 (assessment profiling from WKPELA)
raw_M <- raw_M+0.06
#-------------------------------------------------------------------------------
# 3) create random samples using variance/covariance matrix
......@@ -137,7 +117,7 @@ for(idxIter in 1:nits){
}
biol <- window(window(biol,end=histMaxYr+1),start=histMinYr,end=futureMaxYr) # extend the FLStock object to the full projection period
biol@m.spwn[,ac(2018:2040)] <- 0.67
biol@m.spwn[,ac((histMaxYr+1):futureMaxYr)] <- 0.67
#-------------------------------------------------------------------------------
# 4) create FLStocks object using random samples (with future years as NA)
......@@ -263,7 +243,7 @@ biol@m[,histPeriod] <- as.matrix(raw_M)
#
#-------------------------------------------------------------------------------
#load(file.path(outPath,paste0(assessment_name,'init_MSE_5.RData')))
#load(file.path(modelPath,paste0(assessment_name,'init_MSE_5.RData')))
surveys <- lapply(NSH.tun,propagate,iter=nits)
surveys <- window(window(surveys,end=histMaxYr+1),start=histMinYr,end=futureMaxYr) # extend the FLStock object to the full projection period
......@@ -513,7 +493,7 @@ for(idxFleet in 1:length(fleets)){
# 10) Future recruitment
#-------------------------------------------------------------------------------
recPeriod <- ac(2002:2016)
recPeriod <- ac(2002:2019)
biol.sr <- fmle(as.FLSR(biol,model='segreg')) # just to populate the structure
......@@ -621,50 +601,56 @@ TAC <- FLQuant(NA,dimnames=list(age='all',
area=c('A','B','C','D'),
iter=1:nits))
#- TAC for A fleet in NS
TAC_A <- read.table(file.path(dataPath,'TAC_A.csv'),sep = ",")
TAC[,ac(TAC_A[,1]),,,"A"] <- TAC_A[,2]
#- TAC for B fleet in NS
TAC_B <- read.table(file.path(dataPath,'TAC_B.csv'),sep = ",")
TAC[,ac(TAC_B[,1]),,,"B"] <- TAC_B[,2]
#- TAC for D fleet in WB
TAC_D <- read.table(file.path(dataPath,'TAC_D.csv'),sep = ",")
TAC[,ac(TAC_D[,1]),,,"D"] <- TAC_D[,2]
TAC[,ac((max(TAC_D[,1])+1):(futureMaxYr+3)),,,"D"] <- TAC[,ac(max(TAC_D[,1])),,,"D"] # TAC is fixed for the D fleet
#- TAC for C fleet in WB
TAC_C <- read.table(file.path(dataPath,'TAC_C.csv'),sep = ",")
TAC[,ac(TAC_C[,1]),,,"C"] <- TAC_C[,2]
#- Fixed TAC C in WB (used for transfer to the A fleet)
TAC[,ac((max(TAC_C[,1])+1):(futureMaxYr+3)),,,"C"] <- TAC_C[dim(TAC_C)[1],2]
TACTab <- read.table(file.path(dataPath,'TAC_var','NSAS_TAC.csv'),sep = ",",header=T)
rownames(TACTab) <- TACTab[,1]
# A fleet
TAC[,ac(TACTab$year),,,"A"] <- TACTab$A
# B fleet
TAC[,ac(TACTab$year),,,"B"] <- TACTab$B
# C fleet, fixed TAC
TAC[,ac(TACTab$year),,,"C"] <- TACTab$C
TAC[,ac((max(TACTab$year)+1):(futureMaxYr+3)),,,"C"] <- TAC[,ac((max(TACTab$year))),,,"C"]
# D fleet, fixed TAC
TAC[,ac(TACTab$year),,,"D"] <- TACTab$D
TAC[,ac((max(TACTab$year)+1):(futureMaxYr+3)),,,"D"] <- TAC[,ac((max(TACTab$year))),,,"D"]
#- setup variables for transfer, uptake and split
uptakeFleets <- read.table(file.path(dataPath,'over_underfishing2017.csv'),sep = ",")
uptakeTab <- read.table(file.path(dataPath,'TAC_var','NSAS_uptake.csv'),sep = ",",header=T)
splitTab <- read.table(file.path(dataPath,'TAC_var','NSAS_split.csv'),sep = ",",header=T)
#- Transfer from C fleet TAC to fleet A
Ctransfer <- matrix(runif((length(projPeriod)+3)*nits,min=0.4, max=0.5),nrow=nits,ncol=length(projPeriod)+3) # Transfer of TAC from IIIa to IVa for C fleet in assessment year. Set between 0.4 and 0.5
# update for the D fleets
Duptake <- matrix(1,nrow=nits,ncol=length(projPeriod)+3) # assume full uptake for the D fleet
DSplitHist <- read.table(file.path(dataPath,'D_split.csv'),sep = ",") # get mean and sd from historical data for NSAS/WBSS split for the D fleet
Dsplit <- matrix(rnorm((length(projPeriod)+3)*nits,
mean=mean(DSplitHist$V2),
sd=sd(DSplitHist$V2)/2),
# Transfer of TAC from IIIa to IVa for C fleet in assessment year. Set between 0.4 and 0.5
Ctransfer <- matrix(runif((length(projPeriod)+3)*nits,min=0.4, max=0.5),nrow=nits,ncol=length(projPeriod)+3)
Csplit <- matrix(rnorm((length(projPeriod)+3)*nits,
mean=mean(splitTab$C,na.rm=TRUE),
sd=sd(splitTab$C,na.rm=TRUE)/2),
nrow=nits,ncol=length(projPeriod)+3)
# uptake and split for the D fleets
Duptake <- matrix(1,nrow=nits,ncol=length(projPeriod)+3) # assume full uptake for the D fleet
Dsplit <- matrix(rnorm((length(projPeriod)+3)*nits,
mean=mean(splitTab$D,na.rm=TRUE),
sd=sd(splitTab$D,na.rm=TRUE)/2),
nrow=nits,ncol=length(projPeriod)+3)
#Dsplit <- matrix(rnorm((length(projPeriod)+3)*nits,mean=0.6,sd=0.1),nrow=nits,ncol=length(projPeriod)+3)
# update for the B fleet
# uptake B fleet
Buptake <- matrix(rnorm ((length(projPeriod)+3)*nits,
mean(an(as.vector(uptakeFleets[2:16,3])),na.rm=TRUE), # mean over available historical values
sd(an(as.vector(uptakeFleets[2:16,3])),na.rm=TRUE)/2), # sd over available historical values
mean(an(as.vector(uptakeTab$B)),na.rm=TRUE), # mean over available historical values
sd(an(as.vector(uptakeTab$B)),na.rm=TRUE)/2), # sd over available historical values
nrow=nits,ncol=length(projPeriod)+3)
# store TAC vars
TAC_var <- array(NA,
dim=c(length(projPeriod)+3,nits,4),
dim=c(length(projPeriod)+3,nits,5),
dimnames=list('years' = ac(an(projPeriod)[1]:(an(projPeriod)[length(projPeriod)]+3)),
'iter' = 1:nits,
'var' = c('Ctransfer','Duptake','Dsplit','Buptake')))
'var' = c('Ctransfer','Csplit','Duptake','Dsplit','Buptake')))
TAC_var[,,'Ctransfer'] <- t(Ctransfer)
TAC_var[,,'Duptake'] <- t(Duptake)
TAC_var[,,'Dsplit'] <- t(Dsplit)
TAC_var[,,'Buptake'] <- t(Buptake)
TAC_var[,,'Csplit'] <- t(Csplit)
#-------------------------------------------------------------------------------
......@@ -697,7 +683,7 @@ save( biol, # biology object
NSH.ctrl, # SAM control object
TAC,
TAC_var,
file=file.path(outPath,paste0(assessment_name,'_init_MSE_',ac(nits),'.RData')))
file=file.path(modelPath,paste0(runName,'_MSE_init_',ac(nits),'.RData')))
# resetting parameters
nFutureyrs <- nFutureyrs - 3
......@@ -706,11 +692,10 @@ futureMaxYr <- histMaxYr + nFutureyrs
projPeriod <- ac((histMaxYr+1):futureMaxYr)
fullPeriod <- c(histPeriod,projPeriod)
nFutureyrs <- 40
nFutureyrs <- 20
# save parameters
save(n.retro.years,
nFutureyrs,
save(nFutureyrs,
histMinYr,
histMaxYr,
yearCurrent,
......@@ -722,7 +707,7 @@ save(n.retro.years,
selPeriod,
fecYears,
nits,
file=file.path(outPath,paste0(assessment_name,'_parameters_MSE_',ac(nits),'.RData')))
file=file.path(modelPath,paste0(runName,'_MSE_params_',ac(nits),'.RData')))
......@@ -33,47 +33,7 @@ btrigger <- as.numeric(strsplit(args[[1]][2],'=')[[1]][2])
# HCR
HCR <- ifelse(as.numeric(substr(args[[1]][3],5,nchar(args[[1]][3])))==1,"A","B")
# IAV
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 0){
IAV <- 0
}
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 1){
IAV <- 'A'
}
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 2){
IAV <- 'B'
}
if(as.numeric(substr(args[[1]][4],5,nchar(args[[1]][4]))) == 3){
IAV <- 'E'
}
# BB
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][4]))) == 0){
BB <- 0
}
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][5]))) == 1){
BB <- 'A'
}
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][5]))) == 2){
BB <- 'B'
}
if(as.numeric(substr(args[[1]][5],4,nchar(args[[1]][5]))) == 3){
BB <- 'E'
}
if(IAV == 0){
IAV <- NULL
}else if(IAV == "B"){
IAV <- c("A","B")
}
if(BB == 0){
BB <- NULL
}else if(BB == "B"){
BB <- c("A","B")
}
cat(ftarget,"\t",btrigger,"\t",HCR,"\t",IAV,"\t",BB,"\n")
cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
library(FLSAM)
......@@ -83,15 +43,8 @@ library(minpack.lm) # install.packages("minpack.lm")
library(stats)
# define path to directory
#path <- "D:/Work/Herring MSE/NSAS/"
#path <- "D:/git/wk_WKNSMSE_her.27.3a47d/R/"
#path <- "F:/WKNSMSE/wk_WKNSMSE_her.27.3a47d/R"
#path <- 'E:/wk_WKNSMSE_her.27.3a47d/R'
#path <- 'D:/Repository/NSAS_MSE/wk_WKNSMSE_her.27.3a47d/R/'
#path <- "/home/hintz001/wk_WKNSMSE_her.27.3a47d/R"
path <- '/home/berge057/ICES/wk_WKNSMSE_her.27.3a47d/R/'
#path <- 'E:/git/wk_WKNSMSE_her.27.3a47d/R'
assessment_name <- "NSAS_WKNSMSE2018"
path <- '/home/berge057/ICES/sanoba_nsas'
assessment_name <- "NSAS_sanoba"
try(setwd(path),silent=TRUE)
# paths to different subfolders
......@@ -100,9 +53,18 @@ outPath <- file.path(".","results/")
scriptPath <- file.path(".","side_scripts/")
functionPath <- file.path(".","functions/")
source(file.path(functionPath,"MSE_assessment.R"))
source(file.path(functionPath,"forecastScenarios.r"))
source(file.path(functionPath,"forecastFunctions.r"))
source(file.path(functionPath,"projectNSH_FMSYAR.R"))
# ImY forecast functions
source(file.path(functionPath,"fleet.harvest.r"))
source(file.path(functionPath,"rescaleF.r"))
source(file.path(functionPath,"TAC2sel.r"))
# fmsyAR function
source(file.path(functionPath,"fmsyAR_fun.r"))
source(file.path(functionPath,"find.FCAR.r"))
#-------------------------------------------------------------------------------
# 2) Initialize
......
#-------------------------------------------------------------------------------
# WKNSMSE
#
# Author: Benoit Berges
# WMR, The Netherland
# email: benoit.berges@wur.nl
#
# MSE of North Sea Herring
#
# Date: 2018/11/18
#
# Build for R3.5.1, 64bits
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
# 1) load packages
# setup paths
......@@ -20,13 +6,11 @@
rm(list=ls())
IAV <- NULL
BB <- NULL
HCR <- 'A'
HCR <- 'FMSYAR'
ftarget <- 0.26
btrigger <- 1.4e06
cat(ftarget,"\t",btrigger,"\t",HCR,"\t",IAV,"\t",BB,"\n")
cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
library(FLSAM)
......@@ -35,27 +19,25 @@ library(FLFleet)
library(minpack.lm) # install.packages("minpack.lm")
library(stats)
# define path to directory
#path <- "D:/Work/Herring MSE/NSAS/"
#path <- "D:/git/wk_WKNSMSE_her.27.3a47d/R/"
#path <- "F:/WKNSMSE/wk_WKNSMSE_her.27.3a47d/R"
#path <- 'E:/git/wk_WKNSMSE_her.27.3a47d/R'
#path <- 'D:/Repository/NSAS_MSE/wk_WKNSMSE_her.27.3a47d/R/'
#path <- "/home/hintz001/wk_WKNSMSE_her.27.3a47d/R"
path <- '/home/berge057/ICES/wk_WKNSMSE_her.27.3a47d/R/'
#path <- 'E:/git/wk_WKNSMSE_her.27.3a47d/R'
assessment_name <- "NSAS_WKNSMSE2018"
try(setwd(path),silent=TRUE)
runName <- "NSAS_sanoba"
# paths to different subfolders
dataPath <- file.path(".","data/")
outPath <- file.path(".","results/")
modelPath <- file.path(".","model/")
scriptPath <- file.path(".","side_scripts/")
functionPath <- file.path(".","functions/")
source(file.path(functionPath,"MSE_assessment.R"))
source(file.path(functionPath,"forecastScenarios.r"))
source(file.path(functionPath,"forecastFunctions.r"))
source(file.path(functionPath,"projectNSH_FMSYAR.R"))
# ImY forecast functions
source(file.path(functionPath,"fleet.harvest.r"))
source(file.path(functionPath,"rescaleF.r"))
source(file.path(functionPath,"TAC2sel.r"))
# fmsyAR function
source(file.path(functionPath,"fmsyAR_fun.r"))
source(file.path(functionPath,"find.FCAR.r"))
#-------------------------------------------------------------------------------
# 2) Initialize
......@@ -71,13 +53,14 @@ source(file.path(functionPath,"forecastFunctions.r"))
# - F sel: FAsel, FCsel, FBDsel
#-------------------------------------------------------------------------------
nits <- 1000
nits <- 500
# load object
load(file.path(outPath,paste0(assessment_name,'_init_MSE_',ac(nits),'.RData')))
load(file.path(modelPath,paste0(runName,'_MSE_init_',ac(nits),'.RData')))
stkAssessment.ctrl <- NSH.ctrl
# load MSE parameters
load(file.path(outPath,paste0(assessment_name,'_parameters_MSE_',ac(nits),'.RData')))
load(file.path(modelPath,paste0(runName,'_MSE_params_',ac(nits),'.RData')))
fishery@landings.sel[,projPeriod] <- sweep(fishery@landings.sel[,projPeriod],c(2:6),quantMeans(fishery@landings.sel[,projPeriod]),"/")
strFleet <- c('A','B','C','D')
......@@ -86,6 +69,8 @@ nAges <- dim(biol@stock.n)[1]
surveyNames <- names(surveys)
escapeRuns <- numeric()
runNameSave <- paste0(runName,"_MSE.init_",nits)
#-------------------------------------------------------------------------------
# 2) ref points
# Ftarget = FMSY and Btrigger = MSYBtrigger for now
......@@ -107,8 +92,6 @@ managementRule <- list(HCR = HCR,
TACIAV=IAV, #"A","A+B","NULL"
BB = BB)
runName <- paste0('stkAssessment2018.init_',nits)
#------------------------------------------------------------------------------#
# 3) Housekeeping
#------------------------------------------------------------------------------#
......@@ -117,7 +100,6 @@ CATCH <- TAC
FHCR <- FLQuant(NA, dimnames=list(age=0:8,year=ac(an(projPeriod[1]):(an(projPeriod[length(projPeriod)])+3)),unit="unique",season="all",area="unique",iter=1:nits))
SSBHCR <- FLQuant(NA, dimnames=list(age="all",year=ac(an(projPeriod[1]):(an(projPeriod[length(projPeriod)])+3)),unit="unique",season=c("FcY","CtY"),area="unique",iter=1:nits))
#------------------------------------------------------------------------------#
# 4) Start running the MSE
#------------------------------------------------------------------------------#
......@@ -203,7 +185,7 @@ for (iYr in an(projPeriod)){
datS <- split(dat,as.factor(paste(dat$age,dat$iter)))
res <- foreach(i = 1:length(datS)) %dopar% fitted(loess(data ~ year,data=datS[[i]],span=0.5))
stkAssessment@m <- FLQuant(c(aperm(array(unlist(res),dim=c(length(1947:TaY),nits,nAges)),c(3,1,2))),dimnames=dimnames(stkAssessment@m))
stopCluster(cl); detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE); detach("package:iterators",unload=TRUE)
stopCluster(cl)#; detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE); #detach("package:iterators",unload=TRUE)
# update surveys
for(idxSurvey in surveyNames){
......@@ -222,6 +204,12 @@ for (iYr in an(projPeriod)){
stkAssessment.ctrl@range[5] <- an(TaY)+1
stkAssessment.ctrl@residuals <- F
stk0 <- stkAssessment
idx0 <- stkAssessment.tun
sam0.ctrl <- stkAssessment.ctrl
escapeRuns <- escapeRuns
resInit <- NULL
# Run stock assessment
ret <- MSE_assessment( stkAssessment,
stkAssessment.tun,
......@@ -240,12 +228,14 @@ for (iYr in an(projPeriod)){
#-------------------------------------------------------------------------------
# Forecast
#-------------------------------------------------------------------------------
#(iStocks,iFishery,iYr,iTAC,iHistMaxYr,mpPoints,managementRule)
source(file.path(functionPath,"forecastScenarios.r"))
source(file.path(functionPath,"forecastFunctions.r"))
projNSAS <- projectNSH(stkAssessment,fishery,iYr,TAC,histMaxYr,referencePoints,managementRule)
projNSAS <- projectNSH_FMSYAR(stkAssessment,
fishery,
iYr,
TAC,
TAC_var,
histMaxYr,
referencePoints)
TAC[,FcY,,,c("A","B")] <- projNSAS$TAC[,,,,c("A","B")]
FHCR[,FcY] <- projNSAS$Fbar
SSBHCR[,FcY,,"FcY"] <- projNSAS$SSB$FcY #store HCR SSB in the forecast year
......@@ -279,7 +269,7 @@ for (iYr in an(projPeriod)){
iTer=idxIter,
control=nls.lm.control(maxiter=1000),
lower=rep(1e-8,4))$par)
stopCluster(cl); detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE); detach("package:iterators",unload=TRUE)
stopCluster(cl); detach("package:doParallel",unload=TRUE); detach("package:foreach",unload=TRUE)
#Check for very high F
idx <- which(quantMeans(sweep(landings.sel(fishery[,ac(ImY)]),3:6,t(multso),"*")[ac(2:6),,,,"A"]) > 5)
......
# 04 - Short Term Forecasting with FLasher/04-tutorial.R
# Copyright Iago MOSQUEIRA (WMR), 2021
# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
#
# Distributed under the terms of the EUPL-1.2
# install.packages("FLasher", repos=structure(c(CRAN="https://cloud.r-project.org/",
# FLR="https://flr-project.org/R")))
rm(list=ls())
# Load required packages
library(FLasher)
library(ggplotFL)
# Load example data: ple4
data(ple4)
# --- PREPARE for forecasting
# CALL stf()
# - 3 years ahead.
# - wts, fbar disc ratio: means of 3 last years
fut <- stf(ple4, nyears=3, wts.nyears=3, fbar.nyears=3, disc.nyears=3)
# CHECK averages
yearMeans(stock.wt(ple4)[, ac(2015:2017)]) %/% stock.wt(fut)[, ac(2018:2020)]
# TODO: ALTER future stock.wt with a decreasing 10% trend
stock.wt(fut)[, ac(2018)] <- yearMeans(stock.wt(ple4)[, ac(2015:2017)])*1
stock.wt(fut)[, ac(2019)] <- yearMeans(stock.wt(ple4)[, ac(2015:2017)])*0.9
stock.wt(fut)[, ac(2020)] <- yearMeans(stock.wt(ple4)[, ac(2015:2017)])*0.8
yearMeans(stock.wt(ple4)[, ac(2015:2017)]) %/% stock.wt(fut)[, ac(2018:2020)]
stock.wt(fut)
# 04-tutorial.R - DESC
# 04 - Short Term Forecasting with FLasher/04-tutorial.R
# Copyright Iago MOSQUEIRA (WMR), 2021
# Author: Iago MOSQUEIRA (WMR) <iago.mosqueira@wur.nl>
#
# Distributed under the terms of the EUPL-1.2
# install.packages("FLasher", repos=structure(c(CRAN="https://cloud.r-project.org/",
# FLR="https://flr-project.org/R")))
rm(list=ls())
# Load required packages
library(FLasher)
library(ggplotFL)
# Load example data: ple4
data(ple4)
# --- PREPARE for forecasting
# CALL stf()
# - 3 years ahead.
# - wts, fbar disc ratio: means of 3 last years
fut <- stf(ple4, nyears=3, wts.nyears=3, fbar.nyears=3, disc.nyears=3)
# CHECK averages
yearMeans(stock.wt(ple4)[, ac(2015:2017)]) %/% stock.wt(fut)[, ac(2018:2020)]
# TODO: ALTER future stock.wt with a decreasing 10% trend
stock.wt(fut)[, ac(2018)] <- yearMeans(stock.wt(ple4)[, ac(2015:2017)])*1
stock.wt(fut)[, ac(2019)] <- yearMeans(stock.wt(ple4)[,