Commit de31e9c5 authored by Berges, Benoit's avatar Berges, Benoit
Browse files

update scripts

parent 38b5097d
......@@ -54,8 +54,9 @@ myFLSR <- fmle(myFLSR)
rickerParams <- myFLSR@params
model(myFLSR) <- rickerCa
#model(myFLSR) <- ricker
myFLSR <- fmle(myFLSR,fixed=list(a=47.6))
myFLSR <- fmle(myFLSR)#,fixed=list(a=47.6)
rickerCaParams <- myFLSR@params
......
......@@ -55,7 +55,7 @@ fullPeriod <- c(histPeriod,projPeriod)
recrPeriod <- ac(2011:2021)
selPeriod <- ac(2010:2020)
fecYears <- ac(2010:2020)
nits <- 200 # number of random samples
nits <- 10 # number of random samples
# reading the raw M and applying plus group
......
......@@ -54,7 +54,7 @@ fullPeriod <- c(histPeriod,projPeriod)
recrPeriod <- ac(2011:2021)
selPeriod <- ac(2010:2020)
fecYears <- ac(2010:2020)
nits <- 200 # number of random samples
nits <- 10 # number of random samples
load(file.path(dataPath,paste0("data_sms",SMSkeyRuns,"_sf.RData")))
......
rm(list=ls())
library(icesTAF)
library(ggplotFL)
load('model/NSAS_sanoba_200.RData')
mkdir('figures')
metricsPeriod <- projPeriod[(length(projPeriod)-10):(length(projPeriod)-1)]
compObj <- new('FLStocks')
f26 <- ac(2:6)
f01 <- ac(0:1)
load('model/NSAS_sanoba_option3_1000.RData')
biol@catch <- computeCatch(biol)
biol@stock <- computeStock(biol)
biol@landings <- computeLandings(biol)
compObj[['transfer']] <- biol
load('model/NSAS_sanoba_notrans_1000.RData')
biol@catch <- computeCatch(biol)
biol@stock <- computeStock(biol)
biol@landings <- computeLandings(biol)
compObj[['no_transfer']] <- biol
plot(window(compObj,start=2000))
#metricsPeriod <- projPeriod[(length(projPeriod)-10):(length(projPeriod)-1)]
metricsPeriod <- ac(2031:2039)
f26 <- ac(2:6)
f01 <- ac(0:1)
# risk of SSB < Blim
SSB <- ssb(biol[,metricsPeriod])
......@@ -23,10 +37,10 @@ SSB_bool <- array(FALSE,dim=c(1,nits))
for(idxIter in 1:nits){
# store value per year
SSB_riskMat[which(SSB[,idxIter] <- referencePoints$Blim),idxIter] <- TRUE
SSB_riskMat[which(SSB[,idxIter] < referencePoints$Blim),idxIter] <- TRUE
# TRUE/FALSE for each iteration
if(length(which(SSB[,idxIter] <- referencePoints$Blim)!=0))
if(length(which(SSB[,idxIter] < referencePoints$Blim)!=0))
SSB_bool[idxIter] <- TRUE
}
......@@ -73,7 +87,7 @@ perfMetrics <- cbind(LTR1,LTR3,LTY,LTF01,LTF26)
write.csv(x=perfMetrics,file = 'temp.csv')
# plotting
taf.png("NSAS_risk_ssb")
taf.png("figures/NSAS_risk_ssb")
# plot
par(mfrow=c(2,1))
plot(metricsPeriod,SSBTail[2,],type='l',
......@@ -91,12 +105,10 @@ plot(metricsPeriod,SSB_prob,type='l',
lines(c(2010,2040),c(0.05,0.05),lty=2,col='green')
dev.off()
taf.png("NSAS_risk_ssb")
taf.png("MSE_projections")
taf.png("figures/MSE_projections")
plot(biol)
dev.off()
taf.png("MSE_projections_zoomed")
taf.png("figures/MSE_projections_zoomed")
plot(window(biol,start=2000))
dev.off()
\ No newline at end of file
......@@ -54,7 +54,7 @@ source(file.path(functionPath,"find.FCAR.r"))
# - F sel: FAsel, FCsel, FBDsel
#-------------------------------------------------------------------------------
nits <- 200
nits <- 10
#load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
......@@ -107,7 +107,7 @@ SSBHCR <- FLQuant(NA, dimnames=list(age="all",year=ac(an(proj
start.time <- Sys.time()
#an(projPeriod)
for (iYr in 2032:2042){
for (iYr in 2021:2042){
cat(iYr,"\n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
......
btriggerMat <- seq(from=0.9*1e6,to=1.5*1e6,by=1e5)
ftargetMat <- seq(from=0.18,to=0.22,by=0.01)
for(btrigger in btriggerMat ){
for(ftarget in ftargetMat){
rm(list=setdiff(ls(), c("btriggerMat",'btrigger','ftargetMat','ftarget')))
#-------------------------------------------------------------------------------
# 1) load packages
# setup paths
# load functions
#-------------------------------------------------------------------------------
HCR <- 'FMSYAR'
cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
library(FLSAM)
library(FLEDA)
library(FLFleet)
library(ggplotFL)
library(minpack.lm) # install.packages("minpack.lm")
library(stats)
runName <- "NSAS_sanoba"
saveName <- paste0(runName,'_',ftarget,'_',btrigger,'_',HCR)
# paths to different subfolders
dataPath <- file.path(".","data/")
modelPath <- file.path(".","model/")
scriptPath <- file.path(".","side_scripts/")
functionPath <- file.path(".","functions/")
source(file.path(functionPath,"MSE_assessment.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
#
# Define MSE parameters,
# load objects initialized previously
# * Biology
# - biol
# - surveys
# * Fisheries
# - FCProp
# - catchD
# - F sel: FAsel, FCsel, FBDsel
#-------------------------------------------------------------------------------
nits <- 200
#load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
# load object
load(file.path(modelPath,paste0(runName,'_IBP_MSE_init_',ac(nits),'.RData')))
stkAssessment.ctrl <- NSH.ctrl
# load MSE parameters
load(file.path(modelPath,paste0(runName,'_IBP_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')
nFleets <- length(strFleet)
nAges <- dim(biol@stock.n)[1]
surveyNames <- names(surveys)
escapeRuns <- numeric()
runNameSave <- paste0(runName,"_MSE_",nits)
#-------------------------------------------------------------------------------
# 2) ref points
# Ftarget = FMSY and Btrigger = MSYBtrigger for now
#
#-------------------------------------------------------------------------------
referencePoints <- list(Fmsy = 0.31,
Fsq = NA,
Flim = 0.39,
Fpa = 0.31,
Blim = 874198,
Bpa = 956483,
MSYBtrigger = 1232828,
Ftarget = ftarget,
F01 = 0.05,
Btrigger = btrigger)
managementRule <- list(HCR = HCR)
#------------------------------------------------------------------------------#
# 3) Housekeeping
#------------------------------------------------------------------------------#
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
#------------------------------------------------------------------------------#
start.time <- Sys.time()
#an(projPeriod)
for (iYr in 2021:2042){
cat(iYr,"\n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#----------------------------------------
# define year names
#----------------------------------------
TaY <- ac(iYr-1) # terminal year in the assessment
ImY <- ac(iYr) # intermediate year in the short term forecast, ie, current year
FcY <- ac(iYr+1) # year for which the advice is given, ie, forecast two years ahead the last year in the assessment
FuY <- c(ImY,FcY) # combination of future years
#----------------------------------------
# update the biol number at age in ImY
#----------------------------------------
#- Define mortality rates for iYr-1 to calculate survivors to iYr
m <- m(biol)[,ac(iYr-1),,]
z <- areaSums(landings.sel(fishery)[,ac(iYr-1),,,,]) + m
#- Update biological model to iYr
survivors <- stock.n(biol)[,ac(iYr-1)] * exp(-z)
stock.n(biol)[ac((range(biol,"min")+1):range(biol,"max")),ac(iYr),,] <- survivors[-dim(survivors)[1],,,,,]@.Data
biol@harvest[,ac(iYr-1)] <- areaSums(landings.sel(fishery)[,ac(iYr-1),,,,])
#- Update recruitment
#biol.sr@rec <- rec(biol[,ac(iYr-1)])
biol.sr@ssb[,ac(iYr-1)] <- ssb(biol[,ac(iYr-1)])
recruitBio <- predict(biol.sr[,ac(iYr-1)]) * exp(sr.res[,ac(iYr),drop=T])
stock.n(biol)[1,ac(iYr)] <- recruitBio
#- Plusgroup
if (!is.na(range(biol,"plusgroup"))){
stock.n(biol)[ac(range(biol,"max")),ac(iYr),] <- stock.n(biol)[ac(range(biol,"max")),ac(iYr),] + survivors[ac(range(biol,"max"))]
}
#- Apply process error per cohort age 1 to 8
for(idxAge in 2:nAges){
biol@stock.n[idxAge,ac(iYr)] <- biol@stock.n[idxAge,ac(iYr)]*varProccError[idxAge-1,ac(an(ac(iYr))-(idxAge-1))]
}
cat("\n Finished biology \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#- Update fishery to year iYr-1
landings.n(fishery)[,ac(iYr-1)] <- sweep(sweep(landings.sel(fishery)[,ac(iYr-1),,,,],c(1:4,6),z,"/"),c(1:4,6),stock.n(biol)[,ac(iYr-1)]*(1-exp(-z)),"*")
catch.n(biol)[,ac(iYr-1)] <- areaSums(fishery@landings.n[,ac(iYr-1)])
#print(computeLandings(fishery)[,ac(iYr-1)]/TAC[,ac(iYr-1)])
#-------------------------------------------------------------------------------
# Assessment
#-------------------------------------------------------------------------------
# filter stock object up to intermediate year to include biological variables
stkAssessment <- window(biol,end=an(TaY))
stkAssessment@catch.n <- stkAssessment@catch.n * catchVar[,ac(histMinYr:TaY),,'residuals']
stkAssessment@landings.n <- stkAssessment@catch.n
stkAssessment@landings <- computeLandings(stkAssessment)
#stkAssessment <- window(biol,end=an(TaY))
# smooth M prior to running the assessment, median filter of order 5
require(doParallel); ncores <- detectCores()-1; ncores <- ifelse(nits<ncores,nits,ncores);cl <- makeCluster(ncores); registerDoParallel(cl)
dat <- as.data.frame(stkAssessment@m)
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)
# update surveys
for(idxSurvey in surveyNames){
agesSurvey <- an(rownames(surveys[[idxSurvey]]@index))
yearSurvey <- an(colnames(surveys[[idxSurvey]]@index))
surveyProp <- mean(c(surveys[[idxSurvey]]@range[6],surveys[[idxSurvey]]@range[7]))
surveys[[idxSurvey]]@index[,TaY] <- surveyVars[ac(agesSurvey),TaY,idxSurvey,'catchabilities']*
exp(-z[ac(agesSurvey),TaY]*surveyProp)*
biol@stock.n[ac(agesSurvey),TaY]*surveyVars[ac(agesSurvey),TaY,idxSurvey,'residuals']
}
stkAssessment.tun <- window(surveys,end=an(TaY)+1)
stkAssessment.tun[["HERAS"]] <- window(surveys[["HERAS"]],end=an(TaY))
stkAssessment.tun[["IBTS-Q3"]] <- window(surveys[["IBTS-Q3"]],end=an(TaY))
stkAssessment.ctrl@range[5] <- an(TaY)+1
stkAssessment.ctrl@residuals <- F
#stk0 <- stkAssessment
#idx0 <- stkAssessment.tun
#sam0.ctrl <- stkAssessment.ctrl
#escapeRuns <- escapeRuns
#resInit <- stkAssessment.init
# Run stock assessment
ret <- MSE_assessment( stkAssessment,
stkAssessment.tun,
stkAssessment.ctrl,
escapeRuns) #stkAssessment.init
stkAssessment.init <- ret$resInit
escapeRuns <- ret$escapeRuns
stkAssessment <- ret$stk
#compObj <- new('FLStocks')
#compObj[['assessment']] <- stkAssessment
#compObj[['biol']] <- biol
#windows()
#plot(compObj)
#print(escapeRuns)
#print(stkAssessment@stock.n[,ac(TaY)] / biol@stock.n[,ac(TaY)])
cat("\n Finished stock assessment \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#-------------------------------------------------------------------------------
# Forecast
#-------------------------------------------------------------------------------
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
SSBHCR[,FcY,,"CtY"] <- projNSAS$SSB$CtY #store HCR SSB in the continuation year
cat("\n Finished forecast \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
#-------------------------------------------------------------------------------
# TAC to real F again
#-------------------------------------------------------------------------------
#-Calculate effort accordingly (assuming constant catchability)
#- Get a good starting condition
#mults <- matrix(NA,nrow=nits,ncol=4)
#for(idxIter in 1:nits)
# mults[idxIter,] <- optim(par=runif(4),fn=TAC2sel_V2,iYr=ImY,iBiol=biol[,ImY],iFishery=fishery[,ImY],iTAC=TAC[,ImY],catchVar=catchVar,TAC_var=TAC_var,iTer=idxIter,control=list(maxit=1000),lower=rep(1e-8,4),method="L-BFGS-B")$par
CATCH[,ImY,,,"A"] <- TAC[,ImY,,,"A"] + TAC_var[ImY,,'Ctransfer'] * TAC[,ImY,,,"C"]
CATCH[,ImY,,,"B"] <- TAC[,ImY,,,"B",drop=T] * TAC_var[ImY,,'Buptake',drop=T]
require(doParallel); ncores <- detectCores()-1; ncores <- ifelse(nits<ncores,nits,ncores);cl <- makeCluster(ncores); clusterEvalQ(cl,library(FLCore)); clusterEvalQ(cl,library(minpack.lm)); registerDoParallel(cl)
multso <- do.call(rbind,foreach(idxIter = 1:nits) %dopar% nls.lm(par=runif(4),
fn=TAC2sel,
iYr=ImY,
iBiol=biol[,ImY],
iFishery=fishery[,ImY],
iTAC=CATCH[,ImY],
catchVar=catchVar,
TAC_var=TAC_var,
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)
#catchVar[1,ac(iYr),,"FCprop",,]
#Check for very high F
idx <- which(quantMeans(sweep(landings.sel(fishery)[,ac(ImY)],3:6,t(multso),"*")[ac(2:6),,,,"A"]) > 5)
if(length(idx)>0){
print(idx)
fishery@landings.sel[,ac(ImY),,,,-idx] <- sweep(landings.sel(fishery[,ac(ImY),,,,-idx]),3:6,t(multso)[,-idx],"*")
fishery@landings.sel[,ac(ImY),,,,idx] <- landings.sel(fishery[,ac(an(ImY)-1),,,,idx])
} else {
#When setting a very high F this indicate s that the optimiser doesn't work well, so replace with last year F
fishery@landings.sel[,ac(ImY)] <- sweep(landings.sel(fishery)[,ac(ImY)],3:6,t(multso),"*")
}
cat("\n Finished effort calc \n")
cat(paste("\n Time running",round(difftime(Sys.time(),start.time,unit="mins"),0),"minutes \n"))
}
save.image(file=file.path(modelPath,paste0(saveName,'_',nits,".RData")))
}
}
\ No newline at end of file
......@@ -73,7 +73,7 @@ projectNSH_FMSYAR_no_transfer <- function(iStocks,
iCATCH <- iTAC
iCATCH[,ImY,,,"A"] <- iTAC[,ImY,,,"A"]# + iTAC_var[ImY,,'Ctransfer'] * iTAC[,ImY,,,"C"]
iCATCH[,ImY,,,"B"] <- iTAC[,ImY,,,"B",drop=T]*iTAC_var[ImY,,'Buptake',drop=T]
iCATCH[,ImY,,,"C"] <- (1-iTAC_var[ImY,,'Ctransfer'])#*iTAC[,ImY,,,"C"]*iTAC_var[ImY,,'Csplit']
iCATCH[,ImY,,,"C"] <- iTAC[,ImY,,,"C"]*iTAC_var[ImY,,'Csplit']#*(1-iTAC_var[ImY,,'Ctransfer'])
iCATCH[,ImY,,,"D"] <- iTAC[,ImY,,,"D",drop=T]*iTAC_var[ImY,,'Duptake',drop=T]*iTAC_var[ImY,,'Dsplit',drop=T]
# computing F
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment