library(stockassessment)
load("run/model.RData")
FC<-list()
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.120, .120), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="Fmsy *SSB2027/MSY Btrigger with catch constraint", rec.years=2023:2025, savesim=TRUE)

set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.090, .090), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="MSY approach, Fmsy lower*SSB2027/MSY Btrigger", rec.years=2023:2025, savesim=TRUE)
 
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.138, .138), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="MSY approach, Fmsy upper*SSB2026/MSY Btrigger", rec.years=2023:2025, , savesim=TRUE)

set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.120, .120), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="Fmsy *SSB2026/MSY Btrigger", rec.years=2023:2025, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.182, .182), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="Fmsy", rec.years=2023:2025, savesim=TRUE)
 set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.210, .210), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="Fpa", rec.years=2023:2025, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA, NA, .182), fscale=c(NA,NA,NA,NA), catchval=c(143,115,138,NA), label="TAC2026+20%", rec.years=2023:2025, savesim=TRUE)
 
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA ,NA, NA), fscale=c(NA,NA,NA,NA), catchval=c(143,115,0,0), label="zero catch", rec.years=2023:2025, savesim=TRUE)
#set.seed(123456)
  #FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.123,.315, .315), fscale=c(NA,NA,NA,NA), catchval=c(143,NA,NA,NA), label="other options, Flim", rec.years=2023:2025, savesim=TRUE)
#set.seed(123456)
#  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.17235,NA,.175), fscale=c(NA,NA,NA,NA), catchval=c(1438,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,NA,0.123, NA), fscale=c(NA,NA,NA,1), catchval=c(143,115,NA,NA),label="other options, F2027=F2026", rec.years=2023:2025,       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,NA,NA,.123), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA),nextssb=c(NA,NA,1636,NA),label="SSB(2028)= SSB(2027)", rec.years=2023:2025, savesim=TRUE)
set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,NA,.175), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA),nextssb=c(NA,NA,1800,NA),label="SSB(2028)= SSB(2027)+10%", rec.years=2023:2025, savesim=TRUE)
set.seed(123456) 
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,NA, .123), fscale=c(NA,NA,NA,NA), catchval=c(143,115,2,NA), label="5.0% probability of the spawning stock biomass to fall below Blim in 2027", rec.years=2023:2025, savesim=TRUE)
set.seed(123456)
FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.209, .209), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="MAP range Fupper", rec.years=2023:2025, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA,.136, .136), fscale=c(NA,NA,NA,NA), catchval=c(143,115,NA,NA), label="MAP range Flower", rec.years=2023:2025, 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")