library(stockassessment)
load("run/model.RData")
FC<-list()
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.074, .074), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="MSY approach, Fmsy lower*SSB2026/MSY Btrigger", rec.years=2022:2024, savesim=TRUE)
 
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.114, .114), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="MSY approach, Fmsy upper*SSB2026/MSY Btrigger", rec.years=2022:2024, , savesim=TRUE)

set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.099, .099), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="Fmsy *SSB2026/MSY Btrigger", rec.years=2022:2024, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.182, .182), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="Fmsy", rec.years=2022:2024, savesim=TRUE)
 set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.210, .210), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="Fpa", rec.years=2022:2024, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175, NA, .182), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,244,NA), label="TAC2025+20%", rec.years=2022:2024, savesim=TRUE)
 
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175 ,NA, NA), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,0,0), label="zero catch", rec.years=2022:2024, savesim=TRUE)
#set.seed(123456)
  #FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.315, .315), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="other options, Flim", rec.years=2022:2024, savesim=TRUE)
#set.seed(123456)
#  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,NA,.175), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA),nextssb=c(NA,NA, 1850,NA),label="SSB(2026)= Blim", rec.years=2022:2024, savesim=TRUE)
# this scenario is not possible - SSB(2026)=1850 t
set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,0.175, NA), fscale=c(NA,NA,NA,1), catchval=c(208,NA,NA,NA),label="other options, F=F2024", rec.years=2022:2024,       savesim=TRUE)
#set.seed(123456)
 # FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,NA,.26), fscale=c(NA,NA,NA,NA), catchval=c(2081,NA,NA,NA),nextssb=c(NA,NA,2600,NA),label="SSB(2026)= Bpa,Btrigger", rec.years=2021:2023, savesim=TRUE)
 # this scenario is not possible - SSB(2026)=2600 t
set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,NA,.175), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA),nextssb=c(NA,NA,1353,NA),label="SSB(2027)= SSB(2026)", rec.years=2022:2024, savesim=TRUE)
set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,NA,.175), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA),nextssb=c(NA,NA,1488,NA),label="SSB(2027)= SSB(2026)+10%", rec.years=2022:2024, savesim=TRUE)
set.seed(123456) 
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175 ,NA, .175), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,15,NA), label="5.0% probability of the spawning stock biomass to fall below Blim in 2026", rec.years=2022:2024, savesim=TRUE)
set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.209, .209), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="MAP range Fupper", rec.years=2022:2024, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.175,.136, .136), fscale=c(NA,NA,NA,NA), catchval=c(208,NA,NA,NA), label="MAP range Flower", rec.years=2022:2024, savesim=TRUE)
   

Blim=2094
above<-function(x){xx<-mean(x[[4]]$ssb>Blim);yr<-x[[4]]$year; round(xx*100,1)}
Pabove<-lapply(FC, above)
names(Pabove)<-lapply(FC, function(x){attr(x,"label")})
Pabove<-as.matrix(unlist(Pabove))
colnames(Pabove)<-paste0("P(SSB",FC[[1]][[4]]$year, " > ",Blim,") in %")


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