# Spring SAM forecast (code updated April 2020)

library(stockassessment)

# HC catch TAC for area 4 in  2020, UPDATE
TAC<-17158

# proportion caught in area 4 (HC), UPDATE
propc4<-0.81


#year<-2019
#setwd(paste0("N:\\1_Tanja\\1_MARLAB\\WGNSSK\\Whiting",year,"\\AssessmentNSWhiting\\output_SAM\\NSwhiting_",year,"_update_FC","\\"))
load("run/model.RData")

base.y<-max(fit$data$years)-1  # in spring

Fsq <- mean(faytable(fit)[length(summary(fit)[,1])-1,3:7])
Nsq<-ntable(fit)
Nsq<-Nsq[rownames(Nsq)==base.y,] # N for intermediate year


#ref points
Fmsy     <-0.172  
Fmsylower<-0.158
Flim     <-0.458
Fpa      <-0.330
Fmanage  <-0.15

BMSY <-166708
Bpa  <-166708
Blim <-119970

startR<-2002  # sampling from recruitment series, start year



myforecast <- function(fit, fscale=NULL,Nbase=Nsq, catchval=NULL, catchval.exact=NULL, fval=NULL, nextssb=NULL, landval=NULL, cwF=NULL, nosim=5000, year.base=max(fit$data$years)-1, ave.years=max(fit$data$years)+(-4:0), rec.years=max(fit$data$years)+(-9:0), label=NULL, overwriteSelYears=NULL, deterministic=FALSE, customWeights=NULL, customSel=NULL, lagR=FALSE, splitLD=FALSE, addTSB=FALSE, baselineF=NULL){
  
  resample <- function(x, ...){
    if(deterministic){
      ret <- mean(x)
    }else{
      ret <- x[sample.int(length(x), ...)]
    }
    return(ret)
  }
  ns<-max(length(fscale), length(catchval), length(catchval.exact), length(fval), length(nextssb), length(cwF))    
  if(missing(fscale)&missing(fval)&missing(catchval)&missing(catchval.exact)&missing(nextssb)&missing(cwF))stop("No scenario is specified")    
  if(missing(fscale)) fscale <- rep(NA,ns)
  if(missing(fval)) fval <- rep(NA,ns)
  if(missing(catchval)) catchval <- rep(NA,ns)
  if(missing(catchval.exact)) catchval.exact <- rep(NA,ns)    
  if(missing(nextssb)) nextssb <-rep(NA,ns)
  if(missing(landval)) landval <-rep(NA,ns)  
  if(missing(cwF)) cwF <-rep(NA,ns)  
  if(missing(baselineF)) baselineF <- rep(0,ncol(faytable(fit)))
  
  if(!all(rowSums(!is.na(cbind(fscale, catchval, catchval.exact, fval, nextssb, landval, cwF)))==1)){
    stop("For each forecast year exactly one of fscale, catchval or fval must be specified (all others must be set to NA)")
  }
  
  if(!is.null(overwriteSelYears)){
    fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)  
    Ftab <- faytable(fit)
    fixedsel <- colMeans(Ftab[as.integer(rownames(Ftab))%in%overwriteSelYears,,drop=FALSE])
    fixedsel <- fixedsel/mean(fixedsel[fromto[1]:fromto[2]])
  }
  
  if(!is.null(customSel)){
    fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)  
    customSel <- customSel/mean(customSel[fromto[1]:fromto[2]])
  }
  
  
  getF <- function(x, allowSelOverwrite=FALSE){
    idx <- fit$conf$keyLogFsta[1,]+1
    nsize <- length(idx)
    ret <- exp(x[nsize+idx])
    ret[idx==0] <- 0
    if(allowSelOverwrite){
      if(!is.null(overwriteSelYears)){
        fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)    
        thisfbar<-mean(ret[fromto[1]:fromto[2]])
        ret<-fixedsel*thisfbar
      }
      if(!is.null(customSel)){
        fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)    
        thisfbar<-mean(ret[fromto[1]:fromto[2]])
        ret<-customSel*thisfbar
      }
    }
    ret
  }
  
  fbar<-function(x){
    fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)  
    mean(getF(x)[fromto[1]:fromto[2]])
  }
  
  fbarI<-function(x,baselineF){
    fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)  
    mean(baselineF[fromto[1]:fromto[2]])
  }
  
  
  fbarFrac<-function(x, lf,baselineF){
    fromto <- fit$conf$fbarRange-(fit$conf$minAge-1)  
    mean((lf*(getF(x)-baselineF))[fromto[1]:fromto[2]])
  }    
  
  getCWF<-function(x,w){
    sum(getF(x)*w)
  }    
  
  getN <- function(x){
    idx <- fit$conf$keyLogFsta[1,]+1
    nsize <- length(idx)
    ret <- exp(x[1:nsize])
    ret
  }
  
  getState <- function(N,F){
    k <- fit$conf$keyLogFsta[1,]
    F <- F[k>=0]
    k <- unique(k[k>=0]+1)
    x <- log(c(N,F[k]))
    x
  }    
  
  getProcessVar <- function(fit){
    cof <- coef(fit)
    sdN <- exp(cof[names(cof)=="logSdLogN"][fit$conf$keyVarLogN+1])
    sdN[1]<-0
    nN <- length(sdN)
    sdF <- exp(cof[names(cof)=="logSdLogFsta"][fit$conf$keyVarF+1])
    k<-fit$conf$keyLogFsta[1,]
    sdF <- sdF[k>=0]
    k <- unique(k[k >= 0] + 1)
    sdF <-sdF[k]
    nF <- length(sdF)
    if(fit$conf$corFlag==0){
      corr <- diag(nF)    
    }
    if(fit$conf$corFlag==1){
      y <- cof[names(cof)=="itrans_rho"]
      rho <- 2/(1+exp(-2*y))-1
      corr <- matrix(rho, nrow=nF, ncol=nF)
      diag(corr) <- 1
    }
    if(fit$conf$corFlag==2){
      y <- cof[names(cof)=="itrans_rho"]
      rho <- 2/(1+exp(-2*y))-1
      corr <- diag(sdF)
      corr <- rho^abs(row(corr)-col(corr))
    }
    cov <- matrix(0,nrow=nN+nF,ncol=nN+nF)
    cov[1:nN,1:nN] <- diag(sdN^2)
    cov[nN+1:nF,nN+1:nF] <- (sdF%*%t(sdF))*corr
    cov
  }
  
  step <- function(x, nm, recpool, scale, inyear=FALSE){
    F <- getF(x, allowSelOverwrite=!inyear)
    N <- getN(x)
    if(!inyear){
      Z <- F+nm
      n <- length(N)
      N <- c(resample(recpool,1),N[-n]*exp(-Z[-n])+c(rep(0,n-2),N[n]*exp(-Z[n])))
    }
    F <- F*scale
    xx <- getState(N,F)
    return(xx)
  }
  
  scaleF <- function(x, scale){
    oldF <- getF(x)
    oldFbar <- fbar(x)
    N <- getN(x)
    baseFbar <- fbar(getState(N,baselineF))
    alpha <- oldFbar-baseFbar
    deltaFbar <- (scale-1)*oldFbar
    ss <- 1+deltaFbar/alpha
    F <- baselineF + ss*(oldF-baselineF)
    if(scale==0)F<-baselineF
    F[F<1.0e-9] <- 1.0e-9  
    xx <- getState(N,F)
    return(xx)
  }
  
  sel<-function(x){
    getF(x)/fbar(x)
  }
  
  catch <- function(x, nm, cw){
    F <- getF(x)
    Z <- F+nm
    N <- getN(x)
    C <- F/Z*(1-exp(-Z))*N
    return(sum(cw*C))
  }
  
  catchIBC <- function(x,baselineF, nm, cw){
    F <- getF(x)
    Z <- F+nm
    N <- getN(x)
    C <- baselineF/Z*(1-exp(-Z))*N
    return(sum(cw*C))
  }
  
  catchHC <- function(x,baselineF, nm, cw){
    F <- getF(x)
    Z <- F+nm
    N <- getN(x)
    C <- (F-baselineF)/Z*(1-exp(-Z))*N
    return(sum(cw*C))
  }
  
  catchFrac <- function(x, nm, w, frac, baselineF){
    F <- getF(x)
    Z <- F+nm
    N <- getN(x)
    C <- (F-baselineF)/Z*(1-exp(-Z))*N
    return(sum(frac*w*C))
  }
  
  catchatage <- function(x, nm){
    F <- getF(x)
    Z <- F+nm
    N <- getN(x)
    C <- F/Z*(1-exp(-Z))*N
    return(C)
  }
  
  ssb <- function(x, nm, sw, mo, pm, pf){
    F <- getF(x)
    N <- getN(x)*exp(-pm*nm-pf*F)
    return(sum(N*mo*sw))
  }
  
  tsb <- function(x, sw){
    F <- getF(x)
    N <- getN(x)
    return(sum(N*sw))
  }
  
  rectab<-rectable(fit)
  recpool<-rectab[rownames(rectab)%in%rec.years,1]
  
  # Get final state
  if(year.base==max(fit$data$years)){
    est <- fit$sdrep$estY
    cov <- fit$sdrep$covY
  }
  if(year.base==(max(fit$data$years)-1)){  
    est <- fit$sdrep$estYm1
    cov <- fit$sdrep$covYm1
  }
  if(year.base<(max(fit$data$years)-1)){
    stop("State not saved, so cannot proceed from this year")
  }
  if(deterministic)cov<-cov*0
  sim<-rmvnorm(nosim, mu=est, Sigma=cov)
  sim[,1:9]<-l<-sapply(log(Nbase), rep, nosim)  # overwriting terminal N at age
  
  if(is.null(overwriteSelYears) & is.null(customSel)){  
    if(!all.equal(est,getState(getN(est),getF(est))))stop("Sorry something is wrong here (check code for getN, getF, and getState)")  
  }
  
  doAve  <- function(x,y)colMeans(x[rownames(x)%in%ave.years,,drop=FALSE]) 
  ave.sw <- doAve(fit$data$stockMeanWeight)
  ave.cw <- doAve(fit$data$catchMeanWeight)
  ave.mo <- doAve(fit$data$propMat)
  ave.nm <- doAve(fit$data$natMor)
  ave.lf <- doAve(fit$data$landFrac)
  ave.lw <- doAve(fit$data$landMeanWeight)
  ave.pm <- doAve(fit$data$propM)
  ave.pf <- doAve(fit$data$propF)
  getThisOrAve <- function(x,y, ave){
    if(y %in% rownames(x)){
      ret <- x[which(rownames(x)==y),]
    }else{
      ret <- ave
    }
    ret
  }
  procVar<-getProcessVar(fit)  
  simlist<-list()
  for(i in 0:(length(fscale)-1)){
    y<-year.base+i  
    
    sw<-getThisOrAve(fit$data$stockMeanWeight, y, ave.sw)
    cw<-getThisOrAve(fit$data$catchMeanWeight, y, ave.cw)
    mo<-getThisOrAve(fit$data$propMat, y, ave.mo)
    nm<-getThisOrAve(fit$data$natMor, y, ave.nm)
    lf<-getThisOrAve(fit$data$landFrac, y, ave.lf) 
    lw<-getThisOrAve(fit$data$landMeanWeight, y, ave.lw)
    pm<-getThisOrAve(fit$data$propM, y, ave.pm)
    pf<-getThisOrAve(fit$data$propF, y, ave.pf)
    
    set.seed(1999)
    sim <- t(apply(sim, 1, function(s)step(s, nm=nm, recpool=recpool, scale=1, inyear=(i==0))))
    if(i!=0){
      if(deterministic){procVar<-procVar*0}
      if(!is.null(overwriteSelYears)){nn<-length(fit$conf$keyLogFsta[1,]); procVar[-c(1:nn),-c(1:nn)] <- 0}
      if(!is.null(customSel)){nn<-length(fit$conf$keyLogFsta[1,]); procVar[-c(1:nn),-c(1:nn)] <- 0}
      sim <- sim + rmvnorm(nosim, mu=rep(0,nrow(procVar)), Sigma=procVar)
    }
    
    if(!is.na(fscale[i+1])){
      sim<-t(apply(sim, 1, scaleF, scale=fscale[i+1]))    
    }
    
    if(!is.na(fval[i+1])){
      curfbar<-median(apply(sim, 1, fbar))
      adj<-fval[i+1]/curfbar
      sim<-t(apply(sim, 1, scaleF, scale=adj))    
    }
    
    if(!is.na(cwF[i+1])){
      if(missing(customWeights))stop("customWeights must be supplied when using the cwF option")  
      curcwF<-median(apply(sim, 1, getCWF, w=customWeights))
      adj<-cwF[i+1]/curcwF
      sim<-t(apply(sim, 1, scaleF, scale=adj))    
    }
    
    if(!is.na(catchval[i+1])){
      simtmp<-NA
      fun<-function(s){
        simtmp<<-t(apply(sim, 1, scaleF, scale=s))
        simcat<-apply(simtmp, 1, catchHC, nm=nm, cw=cw, baselineF=baselineF)
        return(catchval[i+1]-median(simcat))
      }
      ff <- uniroot(fun, c(0,100))$root
      sim <- simtmp
    }
    
    if(!is.na(catchval.exact[i+1])){
      simtmp<-sim
      funk<-function(k){
        one<-function(s){
          simtmp[k,]<<-scaleF(sim[k,],s)
          simcat<-catchHC(simtmp[k,], nm=nm, cw=cw, baselineF=baselineF)
          return(catchval.exact[i+1]-simcat)
        }
        ff <- uniroot(one, c(0,100))$root
      }
      dd  <- sapply(1:nrow(sim),funk)
      sim <- simtmp
    }
    
    if(!is.na(landval[i+1])){
      simtmp<-NA
      fun<-function(s){
        simtmp<<-t(apply(sim, 1, scaleF, scale=s))
        simcat<-apply(simtmp, 1, catchFrac, nm=nm, w=lw, frac=lf)
        return(landval[i+1]-median(simcat))
      }
      ff  <- uniroot(fun, c(0,100))$root
      sim <- simtmp
    }
    
    if(!is.na(nextssb[i+1])){
      if((sum(pm)!=0) | (sum(pf)!=0))stop("nextssb option only available if SSB is calculated in the beginning of the year")  
      simtmp<-NA
      
      fun<-function(s){
        simtmp<<-t(apply(sim, 1, scaleF, scale=s))
        set.seed(1999) 
        simsim<-t(apply(simtmp, 1, function(s)step(s, nm=nm, recpool=recpool, scale=1, inyear=(i==0))))
        if(i!=0){
          if(deterministic)procVar<-procVar*0  
          simsim <- simsim + rmvnorm(nosim, mu=rep(0,nrow(procVar)), Sigma=procVar)
        }
        simnextssb<-apply(simsim, 1, ssb, nm=nm, sw=sw, mo=mo, pm=pm, pf=pf)
        return(nextssb[i+1]-median(simnextssb))
      }
      ff  <- uniroot(fun, c(0,100))$root
      sim <- simtmp
    }
    
    fbarsim  <- apply(sim, 1, fbar)
    fbarLsim <- apply(sim, 1, fbarFrac, lf=lf, baselineF=F_ibc)
    fbarIsim <- apply(sim, 1, fbarI, baselineF=F_ibc)
    catchsim <- apply(sim, 1, catch, nm=nm, cw=cw)
    IBCsim   <- apply(sim, 1, catchIBC, nm=nm, cw=cw, baselineF=F_ibc)
    landsim  <- apply(sim, 1, catchFrac, nm=nm, w=lw, frac=lf, baselineF=F_ibc)
    catchatagesim <- apply(sim, 1, catchatage, nm=nm)
    ssbsim   <- apply(sim, 1, ssb, nm=nm, sw=sw, mo=mo, pm=pm, pf=pf)
    tsbsim   <- apply(sim, 1, tsb, sw=sw)
    if(lagR){
      recsim <- exp(sim[,2])
    }else{
      recsim <- exp(sim[,1])
    }
    cwFsim <- rep(NA,nrow(sim))
    if(!missing(customWeights)){
      cwFsim <- apply(sim, 1, getCWF, w=customWeights)
    }
    simlist[[i+1]] <- list(sim=sim, fbar=fbarsim, catch=catchsim,ibc=IBCsim, ssb=ssbsim, rec=recsim,
                           cwF=cwFsim, catchatage=catchatagesim, land=landsim, fbarL=fbarLsim, fbarI=fbarIsim, tsb=tsbsim, year=y)
  }
  
  attr(simlist, "fit")<-fit
  
  collect <- function(x){
    quan <- quantile(x, c(.50,.025,.975))
    c(median=quan[1], low=quan[2], high=quan[3])
  }
  
  fbar  <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$fbar))),5)
  fbarL <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$fbarL))),5)
  fbarI <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$fbarI))),5)
  rec   <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$rec))),1)
  ssb   <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$ssb))),1)
  tsb   <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$tsb))),1)
  catch <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$catch))),1)
  ibc   <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$ibc))),1)
  land  <- round(do.call(rbind, lapply(simlist, function(xx)collect(xx$land))),1)
  
  
  caytable<-round(do.call(rbind, lapply(simlist, function(xx)apply(xx$catchatage,1,collect)))) 
  tab <- cbind(fbar,fbarI,rec,ssb,catch, ibc)
  if(splitLD){
    tab<-cbind(tab,fbarL,fbar-fbarL-fbarI,land,catch-land-ibc)
  }
  if(addTSB){
    tab<-cbind(tab,tsb)
  }
  if(!missing(customWeights)) tab <- cbind(tab,cwF=round(do.call(rbind, lapply(simlist, function(xx)collect(xx$cwF))),3))
  rownames(tab) <- unlist(lapply(simlist, function(xx)xx$year))
  nam <- c("median","low","high")
  basename<-c("fbar:","fbarI:","rec:","ssb:","catch:","ibc:")
  if(splitLD){
    basename<-c(basename,"fbarL:","fbarD:", "Land:","Discard:")    
  }
  if(addTSB){
    basename<-c(basename,"tsb:")    
  }
  if(!missing(customWeights)){
    basename<-c(basename,"cwF:")    
  }
  colnames(tab)<-paste0(rep(basename, each=length(nam)), nam)
  attr(simlist, "tab")<-tab
  shorttab<-t(tab[,grep("median",colnames(tab))])
  rownames(shorttab)<-sub(":median","",paste0(label,if(!is.null(label))":",rownames(shorttab)))
  attr(simlist, "shorttab")<-shorttab
  attr(simlist, "label")   <- label
  attr(simlist, "caytable")<-caytable
  attr(sim,"sim")<-sim
  class(simlist) <- "samforecast"
  simlist
}




#NOTICE MYFORECAST NOT FORECAST


load("run/model.RData")

FC<-list()

#for partial mortality F_ibc

ibcf<-read.ices("data/whg47d_ibf.dat")


# use recent 3 years of fishing selectivity
selYears<-max(fit$data$years) +(-3:-1)  


# to calculate partial Fibc########################
# mean recent 3 years of proportion ibc
ibc_mean<-apply(ibcf[length(ibcf[,1])+(-2:0),],2,mean) 

#estimated total F at age
summary_F<-summary(fit)
f_est    <-summary_F[,7]

#total selectivity
fsel <- faytable(fit)/f_est

# mean recent 3 years of F, scaled to final catch data year
catch_mean<-apply(fsel[length(fsel[,1])+(-3:-1),],2,mean)*f_est[length(f_est)-1] 

IBC_f <-ibc_mean*catch_mean
Fibc  <-mean(IBC_f[3:7])   #age 2-6

F_ibc <-c(IBC_f[[1]],IBC_f[[2]],IBC_f[[3]],IBC_f[[4]],IBC_f[[5]],IBC_f[[6]],IBC_f[[7]],IBC_f[[8]],IBC_f[[9]])



FC<-list()
Fo<-Fibc

#catch option table
for(i in 1:30){
  Fo<-Fo+0.01
  set.seed(12345)
  FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears, addTSB=TRUE, deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fo, Fo),baselineF=F_ibc,label=paste("option",round(Fo,3)))
  
}

save(FC, file="run/forecast_options.RData")





FC<-list()

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fsq,Fsq),baselineF=F_ibc,label="Fsq")                                   
SSBnext<-attributes(FC[[1]])$tab[,"ssb:median"][[3]] # check if above Bmsy

FC<-list()

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fmsy*SSBnext/BMSY,Fmsy*SSBnext/BMSY),baselineF=F_ibc,label="Fmsy*SSB(2021)/MSYBtrigger")                                   

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE, deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fmsy,Fmsy),baselineF=F_ibc,label="Fmsyupper=Fmsy")                                   

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE, deterministic=FALSE,splitLD =TRUE, fval=c(Fsq,Fsq,0,0),baselineF=F_ibc,label="No HC fishery")                                   

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fsq,Fsq), baselineF=F_ibc,label="FSQ")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, catchval.exact=c(NA,NA,TAC/propc4,TAC/propc4), fval=c(Fsq,Fsq,NA, NA), baselineF=F_ibc,label="Roll-over TAC")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, catchval.exact=c(NA,NA,TAC/propc4*0.85,TAC/propc4*0.85), fval=c(Fsq,Fsq,NA, NA), baselineF=F_ibc,label="15%TAC decrease (27.4)")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, catchval.exact=c(NA,NA,TAC/propc4*1.15,TAC/propc4*1.15), fval=c(Fsq,Fsq,NA, NA), baselineF=F_ibc,label="15%TAC increase (27.4)")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE, deterministic=FALSE,splitLD =TRUE, fval=c(Fsq,Fsq,0.75*Fsq, 0.75*Fsq), baselineF=F_ibc,label="0.75*FSQ")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,1.25*Fsq, 1.25*Fsq), baselineF=F_ibc,label="1.25*FSQ")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fpa, Fpa), baselineF=F_ibc,label="Fpa")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Flim, Flim), baselineF=F_ibc,label="Flim")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, nextssb=c(NA,NA,Bpa,Bpa), fval=c(Fsq,Fsq,NA, NA),baselineF=F_ibc,label="Bpa MSYBtrigger")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, nextssb=c(NA,NA,Blim,Blim), fval=c(Fsq,Fsq,NA, NA), baselineF=F_ibc,label="Blim")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fmanage, Fmanage),baselineF=F_ibc,label="Fmanagement")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,Fmsylower,Fmsylower), baselineF=F_ibc,label="Fmsylower")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE, deterministic=FALSE,splitLD =TRUE, fval=c(Fsq,Fsq,Fmsylower*SSBnext/BMSY,Fmsylower*SSBnext/BMSY), baselineF=F_ibc,label="Fmsylower*SSB(2021)/MSYBtrigger")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,0.14*SSBnext/220000, 0.14*SSBnext/220000),baselineF=F_ibc,label="A")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,0.16*SSBnext/200000, 0.16*SSBnext/200000),baselineF=F_ibc,label="B")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE, deterministic=FALSE,splitLD =TRUE, fval=c(Fsq,Fsq,0.14*SSBnext/220000, 0.14*SSBnext/220000),baselineF=F_ibc,label="C")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,0.16*SSBnext/250000, 0.16*SSBnext/250000),baselineF=F_ibc,label="A+D")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,0.16*SSBnext/210000, 0.16*SSBnext/210000),baselineF=F_ibc,label="B+E")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, fval=c(Fsq,Fsq,0.15*SSBnext/230000, 0.15*SSBnext/230000),baselineF=F_ibc,label="C+E")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, catchval.exact=c(NA,NA,TAC/propc4*0.8,TAC/propc4*0.8), fval=c(Fsq,Fsq,NA, NA), baselineF=F_ibc,label="20%TAC decrease (27.4)")

set.seed(12345)
FC[[length(FC)+1]] <- myforecast(fit,year.base=base.y,ave.years=max(fit$data$years) +(-3:-1), rec.years=startR:max(fit$data$years),overwriteSelYears=selYears,addTSB=TRUE,deterministic=FALSE, splitLD =TRUE, catchval.exact=c(NA,NA,TAC/propc4*1.25,TAC/propc4*1.25), fval=c(Fsq,Fsq,NA, NA), baselineF=F_ibc,label="25%TAC increase (27.4)")


save(FC, file="run/forecast.RData")