library(stockassessment)
load("run/model.RData")
FC<-list()


TAC <- 17679

Fmsy <- 0.31
Fmsy_lo <- 0.198
Btrig <- Bpa <- 150000
Blim <- 107000

Ay<-2017:2020
Ry<-1998:2019
Sy<-2017:2019

# Overwrite first two ave.years in landFrac so only the last value is taken forward in forecasts
fit$data$landFrac[(nrow(fit$data$landFrac)-2),] <- fit$data$landFrac[nrow(fit$data$landFrac),]
fit$data$landFrac[(nrow(fit$data$landFrac)-1),] <- fit$data$landFrac[nrow(fit$data$landFrac),]

MSYappr <- function(fit, Fmsy=NULL, Btrig=NULL, fscale_init=1, fval_init=NA, catchval_init=NA, label=NA){
  
  fscale <- fscale_init
  fval <- fval_init
  catchval <- catchval_init
  
  for (i in 2:4){
    
    fscale[i] <- NA
    fval[i] <- Fmsy
    catchval[i] <- NA
    
    set.seed(12345)
    f1 <- forecast(fit, fscale=fscale, fval=fval, catchval=catchval, ave.years=Ay, rec.years=Ry, overwriteSelYears=Sy)
    SSB <- attr(f1,"shorttab")[3,i]
    fval[i] <- min(Fmsy*SSB/Btrig, Fmsy)
    
  }
  
  set.seed(12345)
  f_final <- forecast(fit, fscale=fscale, fval=fval, catchval=catchval, ave.years=Ay, rec.years=Ry, label=label, overwriteSelYears=Sy, splitLD=TRUE)
  return(f_final)
}
FC[[length(FC)+1]] <- MSYappr(fit, Fmsy=Fmsy, Btrig=Btrig, fscale_init = NA, catchval_init = TAC, label="TAC, Fmsy HCR")

FC[[length(FC)+1]] <- MSYappr(fit, F=Fmsy_lo, Btrig=Btrig, fscale_init = NA, catchval_init = TAC, label="TAC, MAP")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), fval=c(NA,0.0000001,0.0000001,0.0000001), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, zero catch", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), fval=c(NA,0.39,0.39,0.39), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, Fpa", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), fval=c(NA,0.54,0.54,0.54), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, Flim", overwriteSelYears=Sy, splitLD=TRUE)

optim_B <- function (f2target, target, fscale_init, fval_init, catchval_init){
  set.seed(12345)
  f1 <- forecast(fit, fscale=c(fscale_init,NA,1,1), fval = c(fval_init,f2target,NA,NA), catchval = c(catchval_init,NA,NA,NA),
                 ave.years=Ay, rec.years=Ry, overwriteSelYears=Sy)
  ssb <- attr(f1,"shorttab")[3,3]
  out <- (ssb - target)^2
  return(out)
}
val <- optim(par=0.5, optim_B, target=Blim, fscale_init=NA, fval_init=NA, catchval_init=TAC,
             method= "Brent", lower=0, upper=1)[[1]]
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,1,1), fval = c(NA,val,NA,NA), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, SSB(2022)=Blim", overwriteSelYears=Sy, splitLD=TRUE)

val <- optim(par=0.2, optim_B, target=Btrig, fscale_init=NA, fval_init=NA, catchval_init=TAC, 
             method= "L-BFGS-B", lower=0, upper=1)[[1]]
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,1,1), fval = c(NA,val,NA,NA), catchval = c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, SSB(2022)=Bpa", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,1,1), fval = c(NA,val,NA,NA), catchval = c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, SSB(2022)=Btrigger", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,0.8*TAC,0.8*TAC,0.8*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)-20%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,0.85*TAC,0.85*TAC,0.85*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)-15%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,0.9*TAC,0.9*TAC,0.9*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)-10%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,0.95*TAC,0.95*TAC,0.95*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)-5%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,TAC,TAC,TAC), ave.years=Ay, rec.years=Ry, label="TAC, constant TAC", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,1.05*TAC,1.05*TAC,1.05*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)+5%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,1.1*TAC,1.1*TAC,1.1*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)+10%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,1.15*TAC,1.15*TAC,1.15*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)+15%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), catchval=c(TAC,1.2*TAC,1.2*TAC,1.2*TAC), ave.years=Ay, rec.years=Ry, label="TAC, TAC(20)+20%", overwriteSelYears=Sy, splitLD=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,1,1,1), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, Fsq", overwriteSelYears=Sy, splitLD=TRUE) 

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), fval = c(NA,0.198,0.198,0.198), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, Fmsy lower", overwriteSelYears=Sy, splitLD=TRUE) 

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), fval = c(NA,0.31,0.31,0.31), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, Fmsy", overwriteSelYears=Sy, splitLD=TRUE)

#set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA), fval = c(NA,0.46,0.46,0.46), catchval=c(TAC,NA,NA,NA), ave.years=Ay, rec.years=Ry, label="TAC, Fmsy upper", overwriteSelYears=Sy, splitLD=TRUE)



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