library(stockassessment)
load("run/model.RData")


## Input Parameters

#  Reference points WKBNSCS
F_MSY <- 0.21
F_MSY_lower <- 0.16
F_MSY_upper <- 0.314
Fpa <- 0.78
Flim <- 1.062
Blim <- 1670
Bpa <- 2322
MSY_Btrigger <- 2322

# TAC & advice for previous year
prev.advice <- 0 
prev.TAC <- 721

# Number of simulations
numberOfSims <- 1e4



##  FORECAST assumptions ####

# Final year in the assessment
final.yr <- max(fit$data$years)

# Resample R from previous 10 years Rec
R.range <-(final.yr-19):(final.yr-1)

# Use five year average for mean weights & L-D partition 
av.yrs <- (final.yr-4):(final.yr)

# Use 5 year average selection patter
sel.yrs <- (final.yr-4):(final.yr)

# F status quo is last observed value
f_sq <- fbartable(fit)[as.character(final.yr), "Estimate"]

# F mean over selected years
f_mean <- mean(fbartable(fit)[as.character(sel.yrs),"Estimate"])

# Final assessment year F - needed to run the forecast through next year to get numbers/SSB at start of intermediate year
f_final <- f_sq



## dummy forecast
# 2022F=Fsq then F_MSY - run a dummy forecast to get the SSB at start of TAC year
#RNGkind(sample.kind = "Rounding")
set.seed(1337)
dummy.forecast <- forecast(fit,
                           fscale = c(NA,NA,NA,NA),
                           fval = c(f_final,f_sq,F_MSY,F_MSY),
                           year.base = final.yr,
                           label = "",
                           rec.years = R.range,
                           ave.years = av.yrs,
                           processNoiseF = FALSE,
                           overwriteSelYears = sel.yrs,
                           splitLD = TRUE,
                           nosim = numberOfSims,
                           savesim = TRUE)

# forecast proper
SSBnext <- attributes(dummy.forecast)$shorttab[3,3]

scen <- list(
  "MSY approach: F_MSY" = list(fval = c(f_final, f_sq, F_MSY, F_MSY)),
  "Precautionary approach: Fpa" = list(fval = c(f_final, f_sq, Fpa, Fpa)),
  "FMSY upper" = list(fval = c(f_final, f_sq, F_MSY_upper, F_MSY_upper)),
  "FMSY lower" = list(fval = c(f_final, f_sq, F_MSY_lower, F_MSY_lower)),
  "F = 0" = list(fval = c(f_final, f_sq, NA, NA), fscale=c(NA,NA,0.0,0.0)),
  "F = Flim" = list(fval = c(f_final, f_sq, Flim, Flim)),
  "Fsq"  = list(fval = c(f_final, f_sq, f_sq, f_sq)),
  "0.05*Fsq" = list(fval = c(f_final, f_sq, 0.05*f_sq,0.05*f_sq)),
  "0.25*Fsq" = list(fval = c(f_final, f_sq, 0.25*f_sq,0.25*f_sq)),
  "0.5*Fsq" = list(fval = c(f_final, f_sq, 0.5*f_sq,0.5*f_sq)),
  "0.75*Fsq" = list(fval = c(f_final, f_sq, 0.75*f_sq,0.75*f_sq)),
  "hit Blim" = list(fval = c(f_final, f_sq, NA,NA), nextssb = c(NA, NA, Blim, Blim)),
  # "hit Bpa" = list(fval = c(f_final, f_sq, NA,NA), nextssb = c(NA, NA, Bpa, Bpa)),
  # "hit MSY Btrigger" = list(fval = c(f_final, f_sq, NA, NA), nextssb = c(NA, NA, MSY_Btrigger, MSY_Btrigger)),
  "SSB Y+2 = SSB Y+1"  = list(fval = c(f_final, f_sq, NA, NA), nextssb = c(NA, NA, SSBnext, SSBnext))
  
)

# Forecast storage list
FC <- list()

# RUN F scenarios ####
for(i in seq(scen)){
  
  #RNGkind(sample.kind = "Rounding")
  set.seed(1337)
  ARGS <- scen[[i]]
  ARGS <- c(ARGS,
            list(fit = fit, 
                 ave.years = av.yrs, 
                 rec.years = R.range, 
                 year.base = final.yr,
                 label = names(scen)[i], 
                 overwriteSelYears = sel.yrs, 
                 splitLD = TRUE,
                 processNoiseF = FALSE,
                 nosim = numberOfSims,
                 savesim = TRUE))  
  try(FC[[i]] <- do.call(forecast, ARGS))
  
  print(paste0("forecast : ", "'", names(scen)[i], "'", " is complete"))
}


# ICES advice rule for MSY
# 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:3){
#     
#     fscale[i] <- NA
#     fval[i] <- Fmsy
#     catchval[i] <- NA
#     
#     set.seed(12345)
#     f1 <- forecast(fit, fscale = fscale, fval = fval, catchval = catchval, ave.years = av.yrs, rec.years = R.range, overwriteSelYears = sel.yrs)
#     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 = av.yrs, 
#                       rec.years = R.range, 
#                       label = label, 
#                       overwriteSelYears = sel.yrs, 
#                       splitLD = TRUE, 
#                       nosim = numberOfSims)
#   return(f_final)
# }
# FC[[length(FC)+1]] <- MSYappr(fit, Fmsy = F_MSY, Btrig = MSY_Btrigger, fscale_init = NA, fval_init = f_sq, catchval_init = NA, label = paste0(final.yr, "F=Fsq then Fmsy HCR"))
# 
# FC[[length(FC)+1]] <- MSYappr(fit, Fmsy = F_MSY_lower, Btrig = MSY_Btrigger, fscale_init = NA, fval_init = f_sq, catchval_init = NA, label = paste0(final.yr, "F=Fsq then Fmsy HCR lower"))
# 
# FC[[length(FC)+1]] <- MSYappr(fit, Fmsy = F_MSY_upper, Btrig = MSY_Btrigger, fscale_init = NA, fval_init = f_sq, catchval_init = NA, label = paste0(final.yr, "F=Fsq then Fmsy HCR upper"))


# Save list of forecasts
save(FC, file = "run/forecast.RData")