# the last year of data => the first year 
library(stockassessment)
load("run/model.RData")

lastyear = max(fit$data$years)

## Recruitment years
Ry <- 2015:lastyear

## biological vars years
# use last 3 years for fishing selectivity and for average weight, maturity, M
Ay <- lastyear +(-2:0) # this is last three years

# TAC of the WG year ? TAC 2020=805 + average of the last three years of discards ... 
#TAC <- 1055

#  Catch per mixed-fisheries scenario : Haddock FMSY lower 1331
#Haddock FMSY
#MixFish <-1854


# ref pts
Fmsy <- 0.29
Flim <- 1.13
Fpa <- 0.77 # equal to F05 if equal to Flim/1.4=0.81
Bpa <- 5800
Blim <- 4200
Btrigger <- 5800
Fmsyupper <- 0.41 
Fmsylower <-0.17
F05 <- 0.77 
Fpa <- min(F05, Fpa)
#Fmsyupper <- min(Fmsyupper, F05)
#Fmsylower <- min(Fmsylower, F05)

## Get Fmix from the summary table
Fsq <- mean(tail(fbartable(fit),3)[,1])
	# F associated to the TAC constraint 
#Fsq <- 0.477

#intermidiate year being mixed fish
#Fmix <-0.735
#Fmix <-1.048

Flast <- tail(fbartable(fit),1)[1]


#Flast <- tail(fbartable(fitSAM),1)[,"Estimate"]
#Fmix <- mean(tail(fbartable(fitSAM),3)[,"Estimate"])

#or for scaled F:
#Fsq <- Flast



# main scenarios (values are for 2020:2023)

FC<-list()

set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then Fmsy", fval = c(Flast,NA, Fmsy, Fmsy), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "MSY approach: Fmsy", fval = c(Flast,Fsq,Fmsy,Fmsy), ave.years=Ay, rec.years=Ry,  splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "MSY approach: Fmsy", fval = c(Flast,Fsq,Fmsy,Fmsy), ave.years=Ay, rec.years=Ry,  splitLD=FALSE, nosim = 100001, savesim = TRUE)

set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then Fmsy upper", fval = c(Flast,NA, Fmsyupper, Fmsyupper), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "FMSY upper", fval = c(Flast,Fsq,Fmsyupper,Fmsyupper), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "FMSY upper", fval = c(Flast,Fsq,Fmsyupper,Fmsyupper), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim = 100001, savesim = TRUE)



set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then Fmsy lower", fval = c(Flast,NA, Fmsylower, Fmsylower), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "FMSY lower", fval = c(Flast,Fsq,Fmsylower,Fmsylower), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "FMSY lower", fval = c(Flast,Fsq,Fmsylower,Fmsylower), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim = 100001, savesim = TRUE)



set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then F=0", fval = c(Flast,NA, 0.000001, 0.000001), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "F = 0", fval = c(Flast,Fsq,0.000001, 0.000001), ave.years=Ay, rec.years=Ry,  splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "F = 0", fval = c(Flast,Fsq,0.000001, 0.000001), ave.years=Ay, rec.years=Ry,  splitLD=FALSE, nosim = 100001, savesim = TRUE)



set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then F=Fpa", fval = c(Flast,NA, Fpa, Fpa), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "Fpa", fval = c(Flast,Fsq,Fpa,Fpa), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "Fpa", fval = c(Flast,Fsq,Fpa,Fpa), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim = 100001, savesim = TRUE)



set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then F=Flim", fval = c(Flast,NA, Flim, Flim), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "F = Flim", fval = c(Flast,Fsq,Flim,Flim), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "F = Flim", fval = c(Flast,Fsq,Flim,Flim), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim = 100001, savesim = TRUE)

set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then F=Flim", fval = c(Flast,NA, Flim, Flim), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "F = F2021", fval = c(Flast,Fsq,Fsq,Fsq), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label = "F = F2024", fval = c(Flast,Fsq,Fsq,Fsq), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim = 100001, savesim = TRUE)



# F status quo depending on intermediate year assumption F=F2020 at wgcse 2020 
#set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then F=Fsq", fval = c(Flast,NA, Fsq, Fsq), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
##set.seed(12345)
##FC[[length(FC)+1]] <- forecast(fit, label = "Fmix", fval = c(Flast,Fmix,Fmix,Fmix), ave.years=Ay, rec.years=Ry,  splitLD=TRUE,nosim = 50001)
##set.seed(12345)
##FC[[length(FC)+1]] <- forecast(fit, label = "Fsq", fval = c(Flast,Fmix,Fsq,Fsq), ave.years=Ay, rec.years=Ry,  splitLD=TRUE,nosim = 50001)


# If archivable otherwise create an error (1/ run F=0 if ssb is higer than reference point can be run) ) 
# at WGCSE 2021 SSB 2022 for Fsq assumptio, thane F=O SSB 2022=3449	with Blim at 4200	
#set.seed(12345)
##FC[[length(FC)+1]] <- forecast(fit, label = "hit Blim", fval = c(Flast,NA, NA, NA),catchval = c(NA,TAC,NA,NA) ,nextssb = c(NA, NA, Blim, Blim), ave.years=Ay, rec.years=Ry,  splitLD=TRUE, nosim = 50001)
#FC[[length(FC)+1]] <- forecast(fit, label = "hit Blim", fval = c(Flast,Fsq, NA, NA), nextssb = c(NA, NA, Blim, Blim), ave.years=Ay, rec.years=Ry,  splitLD=TRUE, nosim = 50001)

# If archivable otherwise create an error (1/ run F=0 if ssb is higer than reference point can be run) ) 
# so not run in 2021
#set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "hit Bpa", fval = c(Flast,NA, NA, NA), catchval = c(NA,TAC,NA,NA) ,nextssb = c(NA, NA, Bpa, Bpa), ave.years=Ay, rec.years=Ry, splitLD=TRUE,  nosim = 50001)
# FC[[length(FC)+1]] <- forecast(fit, label =  "hit Bpa", fval = c(Flast,Fsq, NA, NA), nextssb = c(NA, NA, Bpa, Bpa), ave.years=Ay, rec.years=Ry, splitLD=TRUE,  nosim = 50001)





# # run in a second time 
# to modify manually : SSB in 2024 (at WGCSE 2023) -> done 
# to modify manually : SSB in 2025 (at WGCSE 2024) -> done 
# to modify manually : SSB in 2026 (at WGCSE 2025) -> done 


 SSBadv <- 371
 Fadj <- Fmsy*(SSBadv/Btrigger)
 Flower <- Fmsylower*(SSBadv/Btrigger)
 Fupper <- Fmsyupper*(SSBadv/Btrigger)
###############################################################################################
 set.seed(12345)
## FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then Fmsy HCR ", fval = c(Flast,NA, Fadj,Fmsy), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
# #FC[[length(FC)+1]] <- forecast(fit, label = " Fmsy HCR", fval = c(Flast,Fsq,Fadj,Fmsy), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = 50001)
 FC[[length(FC)+1]] <- forecast(fit, label = " Fmsy HCR", fval = c(Flast,Fsq,Fadj,Fmsy), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim = 100001, savesim = TRUE)


set.seed(12345)
# #FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then Fmsy HCR lower", fval = c(Flast,NA, Flower,Fmsy), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
# #FC[[length(FC)+1]] <- forecast(fit, label = "Fmsy HCR lower", fval = c(Flast,Fsq,Flower,Fmsy), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim =50001)
FC[[length(FC)+1]] <- forecast(fit, label = "Fmsy HCR lower", fval = c(Flast,Fsq,Flower,Fmsy), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim =100001, savesim = TRUE)



set.seed(12345)
# #FC[[length(FC)+1]] <- forecast(fit, label =  "Take TAC+Discards in 2020, then Fmsy HCR upper", fval = c(Flast,NA, Fupper,Fmsy), catchval = c(NA,TAC,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)
# #FC[[length(FC)+1]] <- forecast(fit, label = "Fmsy HCR upper", fval = c(Flast,Fsq,Fupper,Fmsy), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim =50001)
FC[[length(FC)+1]] <- forecast(fit, label = "Fmsy HCR upper", fval = c(Flast,Fsq,Fupper,Fmsy), ave.years=Ay, rec.years=Ry, splitLD=FALSE, nosim =100001, savesim = TRUE)
 
 

## new at WGCSE2021 :) 
set.seed(12345)
##FC[[length(FC)+1]] <- forecast(fit, label =  "hit Bpa", fval = c(Flast,NA, NA, NA), catchval = c(NA,TAC,NA,NA) ,nextssb = c(NA, NA, Bpa, Bpa), ave.years=Ay, rec.years=Ry, splitLD=TRUE,  nosim = 50001)
##FC[[length(FC)+1]] <- forecast(fit, label =  "Stable SSB", fval = c(Flast,Fsq, NA, NA), nextssb = c(NA, NA, SSBadv, SSBadv), ave.years=Ay, rec.years=Ry, splitLD=TRUE,  nosim = 50001)
FC[[length(FC)+1]] <- forecast(fit, label =  "Stable SSB", fval = c(Flast,Fsq, NA, NA), nextssb = c(NA, NA, SSBadv, SSBadv), ave.years=Ay, rec.years=Ry, splitLD=FALSE,  
nosim = 100001, savesim = TRUE)




# #in case one day is needed 
# # set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit, label =  "Take MixFish in 2020, then Fmsy", fval = c(Flast,NA, Fmsy, Fmsy), catchval = c(NA,MixFish,NA,NA), # # ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = 50001)



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




############################################"""
#  ICES advice rule for MSY
#MSYappr <- function(fit, Fmsy, Btrig=NULL, fscale_init=1, fval_init=NA, catchval_init=NA, label=NA){
  
  #fscale <- c(fscale_init,NA)
  #fval <- c(fval_init,Fmsy)
  #catchval <- c(catchval_init,NA)
  
  #for (i in 3:5){
    
   # 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, splitLD=TRUE, nosim=10001)
   # SSB <- attr(f1,"shorttab")[3,i]
    ##fval[i] <- min(Fmsy*SSB/Btrig, Fmsy) ## linear scaling between B=0 and B=Btrigger
    ## linear scaling between F=0 for B=Blim and F=Fmsy for B=Btrigger	
   # a <-  (Fmsy / (Btrig-Blim))
   # b <- -a*Blim
   # if(SSB < Btrig && SSB > Blim) fval[i-1] <- a*SSB + b
 
   # if(SSB < Blim) fval[i-1] <- 0.000001
  #}
  
 # set.seed(12345)
 # f_final <- forecast(fit, fscale=fscale[-5], fval=fval[-5], catchval=catchval[-5], ave.years=Ay, rec.years=Ry, label=label , splitLD=TRUE, nosim= 10001)
 # return(f_final)
#}
#FC[[length(FC)+1]] <- MSYappr(fit, Fmsy=Fmsy, Btrig=Btrigger, fscale_init = NA, fval_init=Fsq, catchval_init = NA, label="F=Fsq then Fmsy HCR")