push changes in relation with IBP

04_runMSE_IBP_v1.r

+#-------------------------------------------------------------------------------
+# 1) load packages 
+#    setup paths 
+#    load functions
+#-------------------------------------------------------------------------------
+
+rm(list=ls())
+
+HCR       <- 'FMSYAR'
+ftarget   <- 0.31
+btrigger  <- 1232828
+
+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"
+
+# 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
+  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),,1]
+  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
+  
+  #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)
+  
+  #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(runName,'_',nits,".RData")))
+
+

04_runMSE_v1.r

@@ -169,7 +169,7 @@ for (iYr in 2021:2042){
   #- 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)])
+  #print(computeLandings(fishery)[,ac(iYr-1)]/TAC[,ac(iYr-1)])
   
   #-------------------------------------------------------------------------------
   # Assessment
@@ -218,14 +218,14 @@ for (iYr in 2021:2042){
   ret <- MSE_assessment(  stkAssessment,
                           stkAssessment.tun,
                           stkAssessment.ctrl,
-                          escapeRuns,
-                          stkAssessment.init)
+                          escapeRuns) #stkAssessment.init
+                          
   stkAssessment.init      <- ret$resInit
   escapeRuns              <- ret$escapeRuns
   stkAssessment           <- ret$stk
   
-  print(escapeRuns)
-  print(stkAssessment@stock.n[,ac(TaY)] / biol@stock.n[,ac(TaY)])
+  #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"))

03_setupMSE_v1.r → attic/03_setupMSE_v1.r

@@ -54,7 +54,7 @@ fullPeriod          <- c(histPeriod,projPeriod)
 recrPeriod          <- ac(2011:2021)
 selPeriod           <- ac(2010:2020)
 fecYears            <- ac(2010:2020)
-nits                <- 2 # number of random samples
+nits                <- 200 # number of random samples
 
 load(file.path(dataPath,paste0("data_sms",SMSkeyRuns,"_sf.RData")))
 
@@ -497,8 +497,9 @@ recPeriod <- ac(2002:2021)
 
 biol.sr <- fmle(as.FLSR(biol,model='segreg')) # just to populate the structure
 
-itersSR <- sample(1:dims(biol)$iter,round(0.15*dims(biol)$iter),replace=F)
-itersRI <- which(!(1:dims(biol)$iter) %in% itersSR)
+params(biol.sr) <- params(fmle(FLSR(rec = rec(biol)[,ac(an(recPeriod)[2]:max(an(recPeriod)))], # rec from 1948 to 2017
+                                    ssb = ssb(biol)[,ac(an(recPeriod)[1]:(max(an(recPeriod))-1))], # ssb from 1947 to 2016
+                                    model='ricker')))
 
 for (its in itersSR){
   iter(params(biol.sr),its)  <- params(fmle(FLSR( rec = rec(iter(biol,its))[,ac(an(recPeriod)[2]:max(an(recPeriod)))], # rec from 1948 to 2017