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


## Input Parameters

#  Reference points - WKNSEA 2021 Benchmark
F_MSY <- 0.21
F_MSY_lower <- 0.173
F_MSY_upper <- 0.21
Fpa <- 0.21
Flim <- 0.31
Blim <- 17286
Bpa <- 25597
MSY_Btrigger <- 25597

# TAC & advice for 2025
prev.advice <- 5116 # 3879 in 2024
prev.TAC <- 4952 # 3163 in 2024

# Current catch
#cur.cat <- 710 # 751 in 2022

# Number of simulations
numberOfSims <- 1e3


###  FORECAST ASSUMPTIONS

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

#  Start year of forecast - needs to be year before intermediate year
start.yr <-final.yr-1

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

# Use three year average for mean weights & L-D partition 
# (see above/below for age 1 & 2 fix, WGCSE 2022)
av.yrs <- (final.yr-3):(final.yr-1)

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

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

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

fit.save <- fit


## WG agreed forecast fixes
# Replace mean weight-at-age 1 for 2020 with 2019 value
fit$data$catchMeanWeight["2020", "1",1] <- fit$data$catchMeanWeight["2019", "1",1]

# 2021 catchMeanWeight & stockMeanWeight ages 0 & 1 replaced with mean of previous 
ags <- as.character(0:1)
rep.yrs <- c(2019, 2020)

# No need to adjust discMeanWeight as not used in the forecast 
# (Total disc = Catch - Land)
fit$data$catchMeanWeight["2021", ags,1] <- colMeans(fit$data$catchMeanWeight[as.character(rep.yrs),ags,1])
fit$data$stockMeanWeight["2021", ags] <- colMeans(fit$data$stockMeanWeight[as.character(rep.yrs),ags])


## 2022 catchMeanWeight & landMeanWeight replaced with 
# mean of previous 3; WGCSE 2023
ags <- as.character(0:7)
rep.yrs <- c(2019:2021)

# No need to adjust discMeanWeight as not used in the forecast 
# (Total disc = Catch - Land)
fit$data$catchMeanWeight["2022", ags,1] <- colMeans(fit$data$catchMeanWeight[as.character(rep.yrs),ags,1])
fit$data$landMeanWeight["2022", ags,1] <- colMeans(fit$data$landMeanWeight[as.character(rep.yrs),ags,1])



### FORECAST

# F scenarios
scen <- list(
  "MSY approach: F_MSY" = list(fval = c(f_sq,f_sq, F_MSY, F_MSY)),
  "Precautionary approach: Fpa" = list(fval = c(f_sq,f_sq, Fpa, Fpa)),
  "FMSY upper" = list(fval = c(f_sq,f_sq, F_MSY_upper, F_MSY_upper)),
  "FMSY lower" = list(fval = c(f_sq,f_sq, F_MSY_lower, F_MSY_lower)),
  "F = 0" = list(fval = c(f_sq,f_sq, NA, NA), fscale=c(NA,NA,0.0,0.0)),
  "F = Flim" = list(fval = c(f_sq,f_sq, Flim, Flim)),
  "Fsq"  = list(fval = c(f_sq,f_sq, f_sq, f_sq)),
  "0.05*Fsq" = list(fval = c(f_sq,f_sq, 0.05*f_sq,0.05*f_sq)),
  "0.25*Fsq" = list(fval = c(f_sq,f_sq, 0.25*f_sq,0.25*f_sq)),
  "0.5*Fsq" = list(fval = c(f_sq,f_sq, 0.5*f_sq,0.5*f_sq)),
  "0.75*Fsq" = list(fval = c(f_sq,f_sq, 0.75*f_sq,0.75*f_sq))
)

# Forecast storage list
FC <- list()

# Run F scenarios
for(i in seq(scen)){
  # RNGkind(sample.kind = "Rounding")
  set.seed(12345)
  ARGS <- scen[[i]]
  ARGS <- c(ARGS,
            list(fit = fit, 
                 ave.years = av.yrs, 
                 rec.years = R.range, 
                 year.base = start.yr,
                 label = names(scen)[i], 
                 overwriteSelYears = sel.yrs, 
                 splitLD = TRUE,
                 nosim = numberOfSims,
                 savesim = TRUE,
                 processNoiseF = FALSE))  
  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 <- rep(fscale_init,2)
  fval <- rep(fval_init,2)
  catchval <- rep(catchval_init,2)
  
  for (i in 3:4){
    
    fscale[i] <- NA
    fval[i] <- Fmsy
    catchval[i] <- NA
    
    set.seed(12345)
    f1 <- forecast(fit, year.base=start.yr,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, 
                      year.base=start.yr,
                      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"))


# SSB target
# X% increase
incr <- 10
inc <- seq(0,50,incr)/100
inc <- 0 # doesn't run for increases
SSBnext <- attributes(FC[[1]])$shorttab[3,2]

# for (i in 1:length(inc)){
#   SSBtarget <- SSBnext*(1+inc[i])
#   # RNGkind(sample.kind = "Rounding")
#   set.seed(12345)
#   FC[[length(FC)+1]] <- forecast(fit, fscale = c(NA,NA,NA), fval = c(f_sq,NA,f_sq), nextssb = c(NA,SSBtarget,NA), year.base = final.yr, label = paste0("2023F=Fsq then ",(i-1)*incr,"% SSB increase"), rec.years = R.range, ave.years = av.yrs, overwriteSelYears = sel.yrs, splitLD = TRUE, nosim = numberOfSims, savesim = TRUE)
# }


# IntermediateF=Fsq then SSB=Blim
# RNGkind(sample.kind = "Rounding")
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale = c(NA,NA,NA,NA),
                               fval = c(f_sq,f_sq,NA,f_sq),
                               nextssb = c(NA,NA,Blim,NA),
                               year.base = start.yr,
                               label = paste0(final.yr, "F=Fsq then SSB=Blim"), 
                               rec.years = R.range,
                               ave.years = av.yrs,
                               overwriteSelYears = sel.yrs,
                               nosim = numberOfSims,
                               splitLD = TRUE,
                               savesim = TRUE)


# IntermediateF=Fsq then SSB=Bpa
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale = c(NA,NA,NA,NA),
                               fval = c(f_sq,f_sq,NA,f_sq),
                               nextssb = c(NA,NA,Bpa,NA),
                               year.base = start.yr,
                               label = paste0(final.yr, "F=Fsq then SSB=Bpa"), 
                               rec.years = R.range,
                               ave.years = av.yrs,
                               overwriteSelYears = sel.yrs,
                               nosim = numberOfSims,
                               splitLD = TRUE,
                               savesim = TRUE)

# IntermediateF=Fsq then SSB2027=SSB2026
set.seed(12345)
SSB_adv.yr <-attributes(FC[[length(FC)]])$shorttab[3,3] # pulls out the SSB at the start of the advice year
FC[[length(FC)+1]] <- forecast(fit,
                               fscale = c(NA,NA,NA,NA),
                               fval = c(f_sq,f_sq,NA,f_sq),
                               nextssb = c(NA,NA,SSB_adv.yr,NA),
                               year.base = start.yr,
                               label = paste0(final.yr, "F=Fsq then SSBy+1=SSBy"), 
                               rec.years = R.range,
                               ave.years = av.yrs,
                               overwriteSelYears = sel.yrs,
                               nosim = numberOfSims,
                               splitLD = TRUE,
                               savesim = TRUE)
# RNGkind(sample.kind = "Rounding")
# set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit, fscale = c(NA,NA,NA), fval = c(f_sq,NA,f_sq), nextssb = c(NA,Bpa,NA), year.base = final.yr, label = paste0(final.yr, "F=Fsq then SSB=Bpa"), rec.years = R.range, ave.years = av.yrs, overwriteSelYears = sel.yrs, nosim = numberOfSims, splitLD = TRUE, savesim = TRUE)


# IntermediateF=Fsq then SSB=MSYBtrigger
# RNGkind(sample.kind = "Rounding")
# set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit, fscale = c(NA,NA,NA), fval = c(f_sq,NA,f_sq), nextssb = c(NA,MSY_Btrigger,NA), year.base = final.yr, label = paste0(final.yr, "F=Fsq then SSB=MSYBtrigger"), rec.years = R.range, ave.years = av.yrs, overwriteSelYears = sel.yrs, nosim = numberOfSims, splitLD = TRUE, savesim = TRUE)


# TAC reduction
# X% increase

inc <- seq(0.25,1,0.25)

for (i in 1:length(inc)){
  Catchtarget <- prev.TAC*inc[i]
  # RNGkind(sample.kind = "Rounding")
  set.seed(12345)
  FC[[length(FC)+1]] <- forecast(fit, 
                                 fval = c(f_sq,f_sq,NA,NA), 
                                 catchval = c(NA,NA,prev.TAC*inc[i],prev.TAC*inc[i]), 
                                 year.base = start.yr,
                                 label = paste0("2025F=Fsq then prev.TAC*",inc[i]), 
                                 rec.years = R.range,
                                 ave.years = av.yrs,
                                 overwriteSelYears = sel.yrs,
                                 splitLD = TRUE,
                                 nosim = numberOfSims,
                                 savesim = TRUE)
}  




# nms <- names(scen)
# LM <-length(scen)
# for (i in (LM+1):length(FC)){
#   nms <-c(nms,attr(FC[[i]],"label"))
# }
# names(FC) <-nms


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