library(stockassessment)
load("run/model.RData")
FC<-list()
Ry <- 2013:2022 # years to sample for recruitment, 10 years and exclude last 2 years
#BioY <- BioY # years for bio data e.g. M but here is as default
CIntY <- 5223 # this is the TAC for 2025 agreed - If reissued advice were followed for 2025 this would be 3194
Fmsy <- 0.25
# Flim <- NA
Fpa <- 0.25
Blim  <- 7352
Bpa  <- 11793
BTrigger  <- 11793
thisY <- max(fit$data$years)+1 # Assessment year

### Helper functions needed when optimization
## Squared difference of SSB at the end of the forecast to a target SSB
ophcr <- function(Fm) {
  args <- c(list(args = I(list(fval = c(NA,NA,Fm,Fm))), nm = paste("F =", Fm)))
  res <- do.call(run_one_scenario, args)
  ssbadvY <- attr(res, "tab")[as.character(thisY),"ssb:median"]
  target <- Fmsy*ssbadvY/BTrigger
  (attr(res, "tab")[as.character(thisY),"fbar:median"] - target)^2
}
## Squared difference of SSB at the end of the forecast to a target SSB
op <- function(Fm, target) {
  args <- c(list(args = I(list(fval = c(NA,NA,Fm,Fm))), nm = paste("F =", Fm)))
  res <- do.call(run_one_scenario, args)
  (attr(res, "tab")[as.character(thisY+1),"ssb:median"] - target)^2
}
## Run one scenario
run_one_scenario <- function(args, nm="") {
  set.seed(123456)
  ARGS <- args
  ARGS <- c(ARGS,
            list(fit = fit, fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label=nm, rec.years=Ry, savesim=TRUE, processNoiseF=FALSE))
  res <- do.call(forecast, ARGS)
  print(paste0("forecast : ", "'", nm, "'", " is complete"))
  res
}
tol <- 0.0001 # tolerance for optimization


# 1st run with F=Fmsy
set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,Fmsy,Fmsy), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="Fmsy", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE) # fish at FMSY from 2026

# ICES MSY approach
SSB_advYMSY <- attributes(FC[[1]])$tab[as.character(thisY),"ssb:median"]
if (SSB_advYMSY<BTrigger) {
#  F_hcr <- optimise(ophcr, c(0, 1), tol = tol)$minimum # estimate F in advice year
  F_hcr<-Fmsy*SSB_advYMSY/BTrigger
  set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,F_hcr,F_hcr), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="MSY approach, ICES HCR", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)
} else {
  set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,Fmsy,Fmsy), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="MSY approach, F=Fmsy", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE) # fish at FMSY from 2025
}

#set.seed(123456)
#FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,Flim,Flim), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="Flim", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

#set.seed(123456)
#FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,Fpa,Fpa), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="Fpa", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

Fsq <- attributes(FC[[1]])$tab[as.character(thisY),1] # F in intermediate year
set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,Fsq,Fsq), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="Fsq", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,NA,NA), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,0,0),  label="zero catch", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)


## Find fishing mortality that leads to Blim, Bpa, Btrigger
F_Blim <- optimise(op, c(0, 1.5), target = Blim, tol = tol)$minimum
F_Bpa <- optimise(op, c(0, 1.5), target = Bpa, tol = tol)$minimum
F_Btr <- optimise(op, c(0, 1.5), target = BTrigger, tol = tol)$minimum

set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,F_Blim,F_Blim), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="SSB_advY=Blim", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,F_Bpa,F_Bpa), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="SSB_advY=Bpa", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,F_Btr,F_Btr), fscale=c(1,NA,NA,NA), catchval=c(NA,CIntY,NA,NA), label="SSB_advY=MSYBtrigger", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

### Extra scenario for 2025 Advice sheet - headline HCR advice but assuming Int Year catch from Reissued advice

  set.seed(123456)
  
### Extra scenario for 2025 Advice sheet - headline HCR advice but assuming Int Year catch from Reissued advice 
set.seed(123456) 
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,Fmsy,Fmsy), fscale=c(1,NA,NA,NA), catchval=c(NA,3206,NA,NA), label="TMP TO GET SSB NOT TO USE", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE) 


SSB_advYMSY <- attributes(FC[[length(FC)]])$tab[as.character(thisY),"ssb:median"] 
F_hcr_int <- Fmsy*SSB_advYMSY/BTrigger 

set.seed(123456) 
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,F_hcr_int,F_hcr_int), fscale=c(1,NA,NA,NA), catchval=c(NA,3206,NA,NA), label="MSY approach, ICES HCR - WITH INT YEAR CATCH FROM REISSUE 2024", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE) 


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