library(stockassessment)
load("run/model.RData")
FC<-list()
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.165, .165), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="MSY approach, Fmsy lower*SSB2025/MSY Btrigger", rec.years=2021:2023, savesim=TRUE)
 
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.226, .226), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="MSY approach, Fp05 as Fmsy upper*SSB2025/MSY Btrigger", rec.years=2021:2023, , savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.226, .226), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="Fmsy *SSB2025/MSY Btrigger", rec.years=2021:2023, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.26, .26), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="Fmsy/Fpa", rec.years=2021:2023, 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(211,NA,392,NA), label="TAC2023*1.2", rec.years=2021:2023, savesim=TRUE)
 
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116 ,NA, NA), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,0,0), label="zero catch", rec.years=2021:2023, savesim=TRUE)
set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.315, .315), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="other options, Flim", rec.years=2021:2023, 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(211,NA,NA,NA),nextssb=c(NA,NA, 1850,NA),label="SSB(2026)= Blim", rec.years=2021:2023, savesim=TRUE)
set.seed(123456)
  FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,0.116, NA), fscale=c(NA,NA,NA,1), catchval=c(211,NA,NA,NA),label="other options, F=F2024", rec.years=2021:2023,       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(211,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,.116,NA,.238), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA),nextssb=c(NA,NA,2256,NA),label="SSB(2026)= SSB(2025)", rec.years=2021:2023, 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(211,NA,15,NA), label="5.0% probability of the spawning stock biomass to fall below Blim in 2026", rec.years=2021:2023, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.26, .26), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="MAP range Fupper", rec.years=2021:2023, savesim=TRUE)
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,.116,.19, .19), fscale=c(NA,NA,NA,NA), catchval=c(211,NA,NA,NA), label="MAP range Flower", rec.years=2021:2023, savesim=TRUE)
   

Blim=1850
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")