library(stockassessment)
load("run/model.RData")
FC<-list()
Ry <- 2011:2021 # 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 <- 7369 # catch in interim year, since we know by now catch for 2023
Fmsy <- 0.35
Flim <- 0.45
Fpa <- 0.35
Blim  <- 6080
Bpa  <- 9064
BTrigger  <- 11208
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,Fm,Fm,Fm))), nm = paste("F =", Fm)))
  res <- do.call(run_one_scenario, args)
  ssbadvY <- attr(res, "tab")[as.character(thisY-1),"ssb:median"]
  target <- Fmsy*ssbadvY/BTrigger
  (attr(res, "tab")[as.character(thisY-1),"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,Fm,Fm,Fm))), nm = paste("F =", Fm)))
  res <- do.call(run_one_scenario, args)
  (attr(res, "tab")[as.character(thisY),"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(NA,NA,NA,NA), catchval=c(CIntY,NA,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,Fmsy,Fmsy,Fmsy), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,NA,NA), label="Fmsy", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE) # fish at FMSY from 2024

# ICES MSY approach
#SSB_advYMSY <- attributes(FC[[1]])$tab[as.character(thisY-1),"ssb:median"]
SSB_advYMSY <- attributes(FC[[1]])$tab[as.character(thisY-1),"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,F_hcr,F_hcr,F_hcr), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,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,Fmsy,Fmsy,Fmsy), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,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,Flim,Flim,Flim), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,NA,NA), label="Flim", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)

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

Fsq <- attributes(FC[[1]])$tab[as.character(thisY-1),1] # F in intermediate year
set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,Fsq,Fsq,Fsq), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,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(NA,NA,NA,NA), catchval=c(CIntY,0,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,F_Blim,F_Blim,F_Blim), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,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,F_Bpa,F_Bpa,F_Bpa), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,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,F_Btr,F_Btr,F_Btr), fscale=c(NA,NA,NA,NA), catchval=c(CIntY,NA,NA,NA), label="SSB_advY=MSYBtrigger", rec.years=Ry, savesim=TRUE, processNoiseF=FALSE)



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