update scripts and data upload

4afc2e8f · Benoit Berges · 6047af60 · 4afc2e8f · 4afc2e8f · 4afc2e8f
Commit 4afc2e8f authored Aug 04, 2021 by Benoit Berges
--- a/02_model_SRR_segReg.R
+++ b/02_model_SRR_segReg.R
+rm(list=ls())
+
+library(FLCore)
+library(TMB)
+
+################################################
+# setup paths
+################################################
+
+data.dir        <- file.path(".","data")
+script.dir      <-  file.path(".","side_scripts")
+function.dir    <-  file.path(".","functions")
+results.dir     <-  file.path(".","results")
+figures.dir     <-  file.path(".","figures")
+
+source(file.path(function.dir,'segReg_SST.R'))
+
+runName <- '01b_SRR_fit'
+
+dir.create(file.path(results.dir,runName),showWarnings = FALSE)
+dir.create(file.path(figures.dir,runName),showWarnings = FALSE)
+
+results.dir     <-  file.path(".","results",runName)
+figures.dir     <-  file.path(".","figures",runName)
+
+################################################
+# read data
+################################################
+
+dat <- read.csv(file.path(results.dir,'..','00_NSAS_SST','SST_NSAS.csv'))
+
+SSB_pred <- seq(from=min(dat$ssb),to=max(dat$ssb),by=1000)
+SST_pred <- seq(from=min(dat$sst),to=max(dat$sst),by=0.5)#min(dat$sst)#
+
+################################################
+# segReg function
+################################################
--- a/03_setupMSE_v1.r
+++ b/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                <- 20 # number of random samples
+nits                <- 2 # number of random samples

 load(file.path(dataPath,paste0("data_sms",SMSkeyRuns,"_sf.RData")))

@@ -504,7 +504,7 @@ 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
                                                  ssb = ssb(iter(biol,its))[,ac(an(recPeriod)[1]:(max(an(recPeriod))-1))], # ssb from 1947 to 2016
                                                  model='segreg')))
-}  
+}
 for (its in itersRI){
  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
                                                 ssb = ssb(iter(biol,its))[,ac(an(recPeriod)[1]:(max(an(recPeriod))-1))], # ssb from 1947 to 2016

--- a/04_runMSE_v1.r
+++ b/04_runMSE_v1.r
-#-------------------------------------------------------------------------------
-# 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 
 #    load functions
 #-------------------------------------------------------------------------------

-#rm(list=ls())
-
-
-args=(commandArgs(TRUE))
-#args <- 'ftar=0.261_btrig=1.4e6_HCR=1_IAV=0_BB=0'
-args    <- strsplit(args,"_")
-#ftarget <- as.numeric(substr(args[[1]][1],6,9))
-#btrigger<- as.numeric(substr(args[[1]][2],7,11))
+rm(list=ls())

-ftarget   <- as.numeric(strsplit(args[[1]][1],'=')[[1]][2])
-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")
+HCR       <- 'FMSYAR'
+ftarget   <- 0.33
+btrigger  <- 1190149

 cat(ftarget,"\t",btrigger,"\t",HCR,"\n")

@@ -39,23 +16,20 @@ cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
 library(FLSAM)
 library(FLEDA)
 library(FLFleet)
+library(ggplotFL)
 library(minpack.lm)  # install.packages("minpack.lm")
 library(stats)

-# define path to directory
-path <- '/home/berge057/ICES/sanoba_nsas'
-assessment_name   <- "NSAS_sanoba"
-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,"projectNSH_FMSYAR.R"))
+source(file.path(functionPath,"MSE_assessment.r"))
+source(file.path(functionPath,"projectNSH_FMSYAR.r"))

 # ImY forecast functions
 source(file.path(functionPath,"fleet.harvest.r"))
@@ -80,19 +54,17 @@ source(file.path(functionPath,"find.FCAR.r"))
 #     - F sel: FAsel, FCsel, FBDsel
 #-------------------------------------------------------------------------------

-nits                <- 1000
+nits                <- 200
+
+load(file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))
+
 # 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 initialization for SAM
-load(file.path(outPath,paste('stkAssessment2018.init_',nits,'.RData',sep="")))
-
 # 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]),"/")
-biol@harvest.spwn[] <- 0.67
-

 strFleet    <- c('A','B','C','D')
 nFleets     <- length(strFleet)
@@ -100,6 +72,7 @@ nAges       <- dim(biol@stock.n)[1]
 surveyNames <- names(surveys)
 escapeRuns  <- numeric()

+runNameSave   <- paste0(runName,"_MSE.init_",nits)

 #-------------------------------------------------------------------------------
 # 2) ref points
@@ -107,27 +80,18 @@ escapeRuns  <- numeric()
 #
 #-------------------------------------------------------------------------------

-referencePoints <- list(Fmsy = 0.26,
+referencePoints <- list(Fmsy = 0.33,
                        Fsq  = NA,
-                        Flim = 0.34,
-                        Fpa  = 0.30,
-                        Blim = 800000,
-                        Bpa  = 900000,
-                        MSYBtrigger = 1400000,
+                        Flim = 0.42,
+                        Fpa  = 0.33,
+                        Blim = 833125,
+                        Bpa  = 1157692,
+                        MSYBtrigger = 1190149,
                        Ftarget = ftarget,
                        F01 = 0.05,
                        Btrigger = btrigger)

-managementRule  <- list(HCR = HCR,
-                        TACIAV=IAV, #"A","A+B","NULL"
-                        BB = BB)
-                        
-runName         <- paste0("NSAS_Ftar_",referencePoints$Ftarget,
-                              "_Btrig_",referencePoints$Btrigger,
-                              "_HCR_",managementRule$HCR,
-                              "_TACIAV_",paste(managementRule$TACIAV,collapse=""),
-                              "_BB_",paste(managementRule$BB,collapse=""),
-                              "_",nits,"iters")
+managementRule  <- list(HCR = HCR)

 #------------------------------------------------------------------------------#
 # 3) Housekeeping
@@ -142,10 +106,11 @@ SSBHCR                    <- FLQuant(NA, dimnames=list(age="all",year=ac(an(proj
 #------------------------------------------------------------------------------#

 start.time          <- Sys.time()
-for (iYr in an(projPeriod)){
+#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
  #----------------------------------------
@@ -153,20 +118,20 @@ for (iYr in an(projPeriod)){
  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 <- array( 0, dim=c(1,nits)) # initialize array
  
@@ -187,66 +152,68 @@ for (iYr in an(projPeriod)){
  
  recruitBio     <- recruitBio * 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(landings.n(fishery[,ac(iYr-1)]))
-  print(computeLandings(fishery[,ac(iYr-1)])/TAC[,ac(iYr-1)])
-
+  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)
-
-#  for(idxIter in 1:nits){
-#    for(idxAge in 1:nAges){
-#      stkAssessment@m[idxAge,,,,,idxIter] <- fitted(loess(data ~ year,data=data.frame(data=stkAssessment@m[idxAge,,,,,idxIter,drop=T],year=histMinYr:TaY),span=0.5))
-#    }
-#  }
-#
+  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']
+      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,
@@ -256,30 +223,37 @@ for (iYr in an(projPeriod)){
  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
  #-------------------------------------------------------------------------------
-  #(iStocks,iFishery,iYr,iTAC,iHistMaxYr,mpPoints,managementRule)
-  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
  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
@@ -288,25 +262,36 @@ for (iYr in an(projPeriod)){
  #  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); detach("package:iterators",unload=TRUE)
-
+  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)
+  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 indicates 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),"*")
+    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=paste(outPath,runName,".RData",sep=""))
+
+save.image(file=file.path(modelPath,paste0(runName,'_',nits,".RData")))


--- a/04_runMSE_v1_init.r
+++ b/04_runMSE_v1_init.r
@@ -7,8 +7,8 @@
 rm(list=ls())

 HCR       <- 'FMSYAR'
-ftarget   <- 0.35
-btrigger  <- 1.1e06
+ftarget   <- 0.33
+btrigger  <- 1190149

 cat(ftarget,"\t",btrigger,"\t",HCR,"\n")

@@ -16,6 +16,7 @@ cat(ftarget,"\t",btrigger,"\t",HCR,"\n")
 library(FLSAM)
 library(FLEDA)
 library(FLFleet)
+library(ggplotFL)
 library(minpack.lm)  # install.packages("minpack.lm")
 library(stats)

@@ -53,7 +54,7 @@ source(file.path(functionPath,"find.FCAR.r"))
 #     - F sel: FAsel, FCsel, FBDsel
 #-------------------------------------------------------------------------------

-nits                <- 20
+nits                <- 200

 # load object
 load(file.path(modelPath,paste0(runName,'_MSE_init_',ac(nits),'.RData')))
@@ -77,13 +78,13 @@ runNameSave   <- paste0(runName,"_MSE.init_",nits)
 #
 #-------------------------------------------------------------------------------

-referencePoints <- list(Fmsy = 0.35,
+referencePoints <- list(Fmsy = 0.33,
                        Fsq  = NA,
-                        Flim = 0.50,
-                        Fpa  = 0.46,
-                        Blim = 651370,
-                        Bpa  = 712681,
-                        MSYBtrigger = 1088428,
+                        Flim = 0.42,
+                        Fpa  = 0.33,
+                        Blim = 833125,
+                        Bpa  = 1157692,
+                        MSYBtrigger = 1190149,
                        Ftarget = ftarget,
                        F01 = 0.05,
                        Btrigger = btrigger)
@@ -104,7 +105,7 @@ SSBHCR                    <- FLQuant(NA, dimnames=list(age="all",year=ac(an(proj

 start.time          <- Sys.time()
 #an(projPeriod)
-for (iYr in 2033: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"))
  
@@ -165,8 +166,8 @@ for (iYr in 2033: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(landings.n(fishery[,ac(iYr-1)]))
-  print(computeLandings(fishery[,ac(iYr-1)])/TAC[,ac(iYr-1)])
+  catch.n(biol)[,ac(iYr-1)]           <- areaSums(fishery@landings.n[,ac(iYr-1)])
+  print(computeLandings(fishery)[,ac(iYr-1)]/TAC[,ac(iYr-1)])
  
  #-------------------------------------------------------------------------------
  # Assessment
@@ -178,6 +179,8 @@ for (iYr in 2033:2042){
  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)
@@ -203,12 +206,6 @@ for (iYr in 2033:2042){
  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,
@@ -227,15 +224,7 @@ for (iYr in 2033:2042){
  #-------------------------------------------------------------------------------
  # Forecast
  #-------------------------------------------------------------------------------
-  
-  iStocks <- stkAssessment
-  iFishery <- fishery
-  cYr <- iYr
-  iTAC <- TAC
-  iTAC_var <- TAC_var
-  iHistMaxYr <- histMaxYr
-  mpPoints <- referencePoints
-  
+
  projNSAS                  <- projectNSH_FMSYAR(stkAssessment,
                                                 fishery,
                                                 iYr,
@@ -243,6 +232,7 @@ for (iYr in 2033:2042){
                                                 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,19 +269,19 @@ for (iYr in 2033:2042){
  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)
+  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 indicates 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),"*")
+    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(stkAssessment.init,file=paste(outPath,runName,".RData",sep=""))
+save(stkAssessment.init,file=file.path(modelPath,paste0(runName,'_stk0_',ac(nits),'.RData')))


--- a/data/M_NSAS_smoothedSpan50_notExtrapolated_sms2010.csv
+++ b/data/M_NSAS_smoothedSpan50_notExtrapolated_sms2010.csv
+"","1974","1975","1976","1977","1978","1979","1980","1981","1982","1983","1984","1985","1986","1987","1988","1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010"
+"0",0.711909279976743,0.747328926556517,0.777391161512859,0.8014387139978,0.820693874724711,0.831917242994105,0.835844562445669,0.839510548732639,0.840188173136646,0.841517535003918,0.852189841044121,0.865641096259169,0.875573905837555,0.88941611878645,0.899808679699398,0.897286712448391,0.891004485980205,0.869649533016433,0.846416648609502,0.827275092130045,0.803056980809397,0.787268461115108,0.778731855420363,0.778134109386344,0.797300138273865,0.823124145205841,0.858105260961612,0.894146809603686,0.922662999551699,0.95204519899448,0.970772450703255,0.979333944473563,0.981569198861359,0.975915625757934,0.963107415167481,0.941047713018257,0.911031759231732
+"1",0.788498589625698,0.821369652981695,0.846036636157837,0.860595268657376,0.868431867016616,0.867567224651472,0.857891174337117,0.846692677542677,0.81744336066516,0.787606175500185,0.782841737073156,0.785405028343215,0.788463550778012,0.808672779250908,0.828348712066839,0.831710226824418,0.830327335915781,0.806623482927639,0.777913731882429,0.746765619685605,0.706400285008127,0.677820549957147,0.663080346193291,0.655998180734566,0.667003951648437,0.683514993666942,0.707638276987537,0.732264855258308,0.748882961491142,0.766162568876112,0.77018325299651,0.760907591872716,0.742914234104105,0.714271426635147,0.676228192665453,0.627138589562851,0.568003324949823
+"2",0.403498373409218,0.404304957224117,0.404692723821615,0.404799873300091,0.404386589034509,0.403482511951381,0.402096727924601,0.39989495469587,0.397726837848113,0.394669848509241,0.388671642294613,0.380842548139685,0.372934032083779,0.361348747843102,0.35024615242151,0.342939737749615,0.337747151077412,0.335783341263826,0.336103306236154,0.343426713410661,0.35531076449164,0.363814738463548,0.364256734137823,0.365474715466232,0.36886043503934,0.373927166679799,0.387391739312215,0.401121482347515,0.415458375901693,0.431562177979423,0.437906664686283,0.43500818831524,0.426326064714309,0.410047222232231,0.387174473940383,0.355549975511711,0.316401340230437
+"3",0.391417443308633,0.390227750280743,0.388806873304313,0.387369234407408,0.385596398356268,0.383536528325625,0.381201374091222,0.378083407599749,0.376063158559585,0.373062077528074,0.366026212671072,0.356959671614864,0.34700339677022,0.329725363693962,0.313088895522503,0.301654761048031,0.293128468122926,0.292541212851112,0.295188096472852,0.305852971047867,0.323133434385292,0.335597035184198,0.339600259538114,0.343133955806108,0.343460566892512,0.344281504396321,0.351400801264859,0.359163565569854,0.371581944990121,0.386613227603993,0.393664281836716,0.393123785562442,0.388212206783113,0.377310849934679,0.361341610894111,0.338355567893724,0.309415514322199
+"4",0.382507432236234,0.376592495957197,0.371195739828037,0.366545882439682,0.36233071708222,0.358568270282166,0.355275333941201,0.351975910126246,0.35099084461025,0.349728494961458,0.34618231418828,0.342149127449851,0.335489262732135,0.319453216184934,0.30321856100632,0.289236016222293,0.277628323315009,0.274122557872133,0.274010411380396,0.28366449855853,0.300296126698552,0.312769235751582,0.31901971360294,0.325133025058136,0.32739349319379,0.329947802435356,0.336931628969672,0.344258784450084,0.356167272793491,0.370325098458795,0.377426312813411,0.377876033795954,0.374419075392056,0.365573613615549,0.352185347524363,0.33247648381083,0.307414013807691
+"5",0.364598961663185,0.357016049804432,0.350062436125918,0.343888597994111,0.338306291105605,0.333324617339459,0.328939858940264,0.324845020195277,0.322678642757604,0.320547244147087,0.317526953934849,0.314699551335001,0.309468478899388,0.296304311328298,0.283016676068646,0.271355757569433,0.261735200468485,0.258075160023098,0.257384139931413,0.266157893296496,0.281175324313224,0.292900894235589,0.299761992920521,0.306798788624579,0.311010677578705,0.315592246591238,0.324054706688011,0.332657264609333,0.345071919719851,0.359419236023959,0.36617046877253,0.365678041276735,0.36079989804867,0.350053226722886,0.334297580360587,0.311786272518228,0.283466546797617
+"6",0.354383431252432,0.347074711402721,0.340418070225491,0.334582062340488,0.329286117591356,0.324732066264222,0.320889828548741,0.317063241561416,0.314437113449411,0.311676573709672,0.307476178522576,0.302896483800144,0.29701080828102,0.285172919008874,0.273576761531216,0.264345202270051,0.257312888052114,0.255865110337759,0.257192476505213,0.266213147404261,0.280549472370559,0.291866965111225,0.299053738751486,0.306406079652275,0.311294176842882,0.316301503462095,0.323830084474292,0.33128828233089,0.342593477794689,0.35563367552911,0.36147201734057,0.360468989812514,0.355195948468298,0.344181565634077,0.328315489160229,0.305969235904955,0.278019627721524
+"7",0.346146588956649,0.33913721937548,0.332700780747395,0.327021130522752,0.32172572485234,0.31711732775281,0.313167246613152,0.308959418092181,0.305331271746348,0.301410646244555,0.295142818560609,0.287647925968895,0.280207809247768,0.269095682450023,0.258845087649201,0.252540464731573,0.248777777700291,0.249754053178957,0.253321851707798,0.263223278096437,0.277812336181131,0.289429794361565,0.297363217759803,0.305261888824613,0.309936526286864,0.314574596553208,0.321621864058668,0.328590820145814,0.339861191012972,0.352935746569154,0.358547608645135,0.3570496037224,0.351125456390757,0.339265550786822,0.322385815222024,0.298803625782192,0.269414369591271
+"8",0.346146588956649,0.33913721937548,0.332700780747395,0.327021130522752,0.32172572485234,0.31711732775281,0.313167246613152,0.308959418092181,0.305331271746348,0.301410646244555,0.295142818560609,0.287647925968895,0.280207809247768,0.269095682450023,0.258845087649201,0.252540464731573,0.248777777700291,0.249754053178957,0.253321851707798,0.263223278096437,0.277812336181131,0.289429794361565,0.297363217759803,0.305261888824613,0.309936526286864,0.314574596553208,0.321621864058668,0.328590820145814,0.339861191012972,0.352935746569154,0.358547608645135,0.3570496037224,0.351125456390757,0.339265550786822,0.322385815222024,0.298803625782192,0.269414369591271
--- a/data/M_NSAS_smoothedSpan50_notExtrapolated_sms2013.csv
+++ b/data/M_NSAS_smoothedSpan50_notExtrapolated_sms2013.csv
+"","1974","1975","1976","1977","1978","1979","1980","1981","1982","1983","1984","1985","1986","1987","1988","1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013"
+"0",0.907405801131202,0.899858966040323,0.891359015574735,0.881754019480869,0.871246444973792,0.859335150192596,0.846152649833222,0.832794713983408,0.818821560416312,0.804453806659948,0.789080535747094,0.772499809979678,0.757411467981446,0.74150609840905,0.726663321737277,0.717325155843505,0.711055125889632,0.706518879153387,0.702455861389356,0.699894014370123,0.697479359104755,0.695430716523253,0.698205677811742,0.709309539905845,0.723412483878173,0.743136132380887,0.768764466061742,0.787995939305731,0.801802092550251,0.813058368532596,0.816353034265805,0.815937009834802,0.816858582682246,0.819236798961927,0.8204170055398,0.82039733896306,0.819067365017853,0.817151966073525,0.815444432074909,0.813415094003822
+"1",0.68804351171759,0.705489030452272,0.718915622125983,0.72757246373009,0.7329135018218,0.734300234800414,0.731720394979385,0.727514153746058,0.715447108223389,0.701167668553565,0.693739533014702,0.687195929480771,0.679155543261739,0.671746298108385,0.663680415406921,0.651400378780279,0.637624126393436,0.624259971822577,0.60381108471884,0.58386191251299,0.569598780101773,0.555151112920747,0.548463792023747,0.55322929577711,0.562814066435262,0.582765839549288,0.612076175699277,0.633366639314701,0.653486164325781,0.669827727075335,0.667348075534228,0.65776602338416,0.64966229723974,0.64255242170917,0.633525527216188,0.623835428471458,0.612149277072566,0.599528213687807,0.587311942061333,0.574678105817843
+"2",0.399999878350799,0.397859996818929,0.395459693786251,0.392928200683458,0.390081916399468,0.38687231700314,0.3833194263769,0.379288238361821,0.37479937259754,0.369967863567501,0.365059821040196,0.360040954930909,0.354467281566614,0.344420981727771,0.334450349983996,0.327654392392285,0.32083907460775,0.316820023038271,0.317761659470726,0.319782380492531,0.320923905144194,0.32337512948549,0.327103407654012,0.329993665222228,0.334587421008249,0.346050948546694,0.360916431021353,0.373184976912042,0.388570454506911,0.402106610655478,0.405556459595899,0.406068671830068,0.404702658075469,0.401198174841843,0.395417539393698,0.387670335493719,0.377601369215966,0.365456725877138,0.351692267287213,0.336128656595295
+"3",0.372108763284311,0.368115705230617,0.363837551278369,0.3593555906453,0.354510788945591,0.349262492995268,0.343630864874212,0.337543570432764,0.331310078712127,0.324676043820107,0.316022228227314,0.306092364259217,0.297512850148345,0.286266931327914,0.276133276334888,0.272256511559531,0.270442255023545,0.271001527062243,0.277756835167496,0.285697484315334,0.291335666537111,0.298411856985948,0.304829696906052,0.308965057987113,0.313526296090113,0.320267707479592,0.327449096061028,0.334231090796974,0.344917742316553,0.354470486006083,0.357712097075634,0.359337654769942,0.359302321239352,0.357505886870993,0.35413299660688,0.349294257907935,0.342846638331962,0.334792122006966,0.325289232429114,0.314348121590629
+"4",0.343114440418075,0.343053377312498,0.342032695960901,0.340050544991538,0.33710001794958,0.33311825167096,0.328105302890856,0.322166710261044,0.314853612660096,0.306593954550774,0.296193536048503,0.283996988246207,0.273733304638195,0.261232169497742,0.250104023581197,0.246539764978701,0.2457710998614,0.246759368569732,0.253499674965397,0.261258037616513,0.266060944416955,0.271983774564047,0.276967577413015,0.278318723977471,0.280042024628859,0.284281394076418,0.288496281643155,0.293648697641689,0.304399401413741,0.314795844796743,0.320378827982675,0.325491763281888,0.328380617928131,0.328967812001897,0.328032336631266,0.325449389619041,0.32131070091313,0.315471892536321,0.307921354714511,0.298788554813923
+"5",0.326272105287708,0.323735444114562,0.320617398570114,0.316901782870969,0.312595462027075,0.307697816795551,0.302191371425915,0.296101238056409,0.289135918710525,0.281546139197024,0.272235622031789,0.261591230224009,0.252676122329792,0.242709014650565,0.233849582887326,0.230128986570232,0.228447151244637,0.228581823588789,0.233528560280498,0.239532822286264,0.243602020809531,0.248658875037927,0.253228950059015,0.255425061328839,0.258257699688698,0.263539361753165,0.269302844811365,0.275750678040565,0.286701893889181,0.297260759921269,0.303855229252599,0.310124755630239,0.314016971851875,0.315531814213547,0.315489508965398,0.31367342408535,0.310268830532017,0.305056650505279,0.297935783029926,0.289095252664825
+"6",0.322318357600349,0.31785329691615,0.313166595154038,0.308289300249161,0.303127019912318,0.297746017561013,0.292127035574469,0.286103746096554,0.280014087133928,0.27351088178779,0.2648273158086,0.254811310572884,0.246190969305414,0.236311728588427,0.227379589585914,0.222976965836057,0.220225056380436,0.219450989571552,0.223492888896782,0.228608876626146,0.231903392991665,0.236144146315371,0.240283165258747,0.24289323227899,0.246295416411159,0.25225629809469,0.259029098149267,0.265918235673532,0.276514320710165,0.286573135855317,0.292280195851972,0.297318619106388,0.30048698904233,0.301759527426825,0.301643890894114,0.30002951336758,0.296996715454777,0.292430562906539,0.286333597564199,0.278810475870855
+"7",0.303604124401098,0.295470345888574,0.287652929959825,0.280204880277436,0.272976041863888,0.26607253988579,0.259476309732555,0.252917707046332,0.247239477102996,0.241593796409321,0.233593877488211,0.224494771372392,0.216953738315093,0.209229927632296,0.202459256275314,0.19956430405247,0.198466818241458,0.198637961752548,0.202363829202654,0.20686568106453,0.209314162532982,0.21212822214232,0.215054226595817,0.215833735062402,0.217642406733304,0.22358896147668,0.230854759186614,0.238446276638736,0.250261621168791,0.261735347867447,0.268984199452999,0.27588251910965,0.280377855944278,0.282514712335324,0.283025497827845,0.281621969703983,0.278557496309099,0.273623904094817,0.266726631317303,0.258049129156957
+"8",0.303604124401098,0.295470345888574,0.287652929959825,0.280204880277436,0.272976041863888,0.26607253988579,0.259476309732555,0.252917707046332,0.247239477102996,0.241593796409321,0.233593877488211,0.224494771372392,0.216953738315093,0.209229927632296,0.202459256275314,0.19956430405247,0.198466818241458,0.198637961752548,0.202363829202654,0.20686568106453,0.209314162532982,0.21212822214232,0.215054226595817,0.215833735062402,0.217642406733304,0.22358896147668,0.230854759186614,0.238446276638736,0.250261621168791,0.261735347867447,0.268984199452999,0.27588251910965,0.280377855944278,0.282514712335324,0.283025497827845,0.281621969703983,0.278557496309099,0.273623904094817,0.266726631317303,0.258049129156957
--- a/data/M_NSAS_smoothedSpan50_notExtrapolated_sms2016.csv
+++ b/data/M_NSAS_smoothedSpan50_notExtrapolated_sms2016.csv
+"","1974","1975","1976","1977","1978","1979","1980","1981","1982","1983","1984","1985","1986","1987","1988","1989","1990","1991","1992","1993","1994","1995","1996","1997","1998","1999","2000","2001","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016"
+"0",0.70821379746487,0.700420307927136,0.695702476022895,0.69473436677739,0.696492459647188,0.700220876839265,0.707217259803183,0.717491370939259,0.728355513880606,0.749058091796919,0.770926698351733,0.783412990390585,0.793727626549795,0.798940198889404,0.795554350289577,0.786964879180992,0.780578110949461,0.774778919418454,0.766525786608205,0.761085450199674,0.75534646001166,0.751198830654995,0.75441305732743,0.761415324438655,0.77095810067559,0.7853693292209,0.803822460413593,0.82033822367511,0.836227081824926,0.853640587643627,0.867861105494315,0.882440395298809,0.89515374395856,0.90261689066062,0.907615667174079,0.906022867204554,0.89776510080696,0.884524805844719,0.865860971132433,0.842214289148341,0.813413850806278,0.779131144689247,0.739540932550355
+"1",0.550199820855395,0.573330600342734,0.592828614845479,0.608218668036437,0.620141614965081,0.629109162368136,0.634325553076614,0.635818974946859,0.635337052253591,0.629126680764678,0.620425864900497,0.613087128610219,0.603888467931139,0.594452378248937,0.585667023624339,0.577216982448651,0.568919371603427,0.557081878497098,0.542547211653838,0.531935019628814,0.522277698970879,0.515338090790473,0.516694610358206,0.522309448213281,0.531473582656131,0.551205581321105,0.578789176111224,0.598913091729937,0.613712417437998,0.62918905928278,0.63583196714766,0.625124985760259,0.611273102517495,0.598987051532361,0.582702834973932,0.571579018813663,0.564276236969545,0.558272357747997,0.555707350838682,0.554592461905379,0.556119022463369,0.561635026008606,0.570009205963726
+"2",0.328516959703639,0.334817622797411,0.339680284629857,0.343027052803792,0.34495420344434,0.345531941554381,0.344443746425189,0.34174330365559,0.338056127709786,0.332887823186219,0.326609970800872,0.317536405272719,0.305990403732243,0.296910608437214,0.289645710621851,0.28248702628463,0.278787134829706,0.281021444330276,0.286585728592264,0.29121311561891,0.296292758630242,0.301901873663001,0.304883562286297,0.307545114429399,0.311848034860861,0.319140474490306,0.327679672354979,0.334790765414868,0.342878274840535,0.352360335127498,0.357582503223176,0.357093986034704,0.355199279356283,0.350617708989263,0.343419373004566,0.338608920084563,0.335584319175395,0.332726984624054,0.331007152405116,0.329510785627688,0.328891291777214,0.329889225859059,0.33190122788239
+"3",0.276135716202933,0.281502707091254,0.285333604582301,0.287428970600285,0.288000100884173,0.287260647852624,0.284841136867152,0.280761154887549,0.275817946729892,0.268207246986042,0.25942083277674,0.249299340513585,0.237202795345938,0.228151588989836,0.22277478871863,0.21873911102337,0.216621421446644,0.218960719193309,0.224640261800461,0.228753302309033,0.231205070539922,0.233476387469932,0.232821011514662,0.230763311033142,0.231165840715298,0.234720801948247,0.239416680112059,0.245060455886088,0.25414458145426,0.26568002094063,0.274251522934744,0.281723484275779,0.288206929687492,0.289276127868085,0.287923451552273,0.287230747211246,0.286962510961263,0.286031631376554,0.284932413613844,0.283192817825234,0.281167761661924,0.279318837649081,0.277321133840143
+"4",0.250468966257871,0.252641952620312,0.253829565321258,0.253913335936614,0.252969707402391,0.251132181243749,0.248328084598335,0.24454169242238,0.239991927968921,0.233817372360543,0.226628481801505,0.216615482044335,0.204199835416126,0.195155132654283,0.190110183550613,0.186426425241792,0.184976389482522,0.189017383757321,0.197034673510084,0.202742727434571,0.206967508928316,0.210561528505646,0.209307685067034,0.205897444872407,0.204447948329908,0.204006086292109,0.202951653484201,0.204401146123263,0.210626029747781,0.219412939696939,0.227151126888864,0.23698160917727,0.246731096854219,0.251487050548581,0.254867387224352,0.258085580959863,0.261064038445993,0.263416981770574,0.265354473805804,0.266678960747518,0.267530214754928,0.268123229401913,0.268327195001408
+"5",0.234473269801187,0.234464138753423,0.23386147145836,0.232578976038424,0.230672586324601,0.228241153123978,0.225275081115743,0.221754537506429,0.217763016894132,0.212802755318585,0.2071951729399,0.199647449282198,0.19049445725218,0.183490007702388,0.178962645295857,0.175292481544041,0.173310677572101,0.174870175290321,0.178922823809596,0.181948145056,0.185006550532776,0.187800363259124,0.187011801061815,0.184949380323408,0.184494938077612,0.184968889008533,0.1852191601773,0.187418857183152,0.193271613884862,0.20115394847414,0.208248293049634,0.217045731999811,0.225845827665294,0.230891879771927,0.234971110342937,0.23879455185128,0.242298961368831,0.245269790610803,0.247835694911712,0.249876103893813,0.251501253426681,0.25287646463801,0.253900935685525
+"6",0.226462538147919,0.224001420717782,0.221304042205047,0.21831864711188,0.215064071973297,0.21160475306262,0.208007388839063,0.204243177396947,0.200233304462258,0.19573917828349,0.190870708626046,0.184522955393262,0.17707047863342,0.171176466565288,0.167144044032085,0.163817572654042,0.161672095815447,0.161910384417184,0.163852951658839,0.165151454372747,0.16596965170766,0.166805646551369,0.165960357575744,0.164463603628726,0.164561907415313,0.166078056381,0.167968294728731,0.171179330299677,0.177158994542389,0.184987425707524,0.192059791839414,0.199895361535582,0.207695325153169,0.212850019001856,0.217261287032538,0.221614441390796,0.22583009574468,0.229732183307487,0.233475562259302,0.236915667315034,0.240143654740364,0.243292986604323,0.246277460132593
+"7",0.221422133530034,0.217948579647922,0.214359479032202,0.210619119194424,0.206735653360268,0.202754877421521,0.198736314118315,0.194656255156843,0.190435649877971,0.185868624360861,0.181063495001107,0.175046746457125,0.168154911672982,0.162745297892671,0.159091330885199,0.156145955042085,0.154285755541088,0.154660709483836,0.156657101890314,0.158005979815813,0.1587452