library(stockassessment)
load("run/model.RData")
FC<-list()

## Establish some common settings
DataYear <- as.numeric(format(Sys.time(), "%Y"))-1
nruns <- 5000 #for exploratory runs, reduce the nbr of simulations (~5000). For final results, increase the nbr (min 50,000)

## Reference points
MSYBtrigger <- 13460
Bpa <- 13460
Blim <- 11119
SSB_5pcMSY <- 11251.37
Fpa <- 0.149
Fmsy <- 0.149
Fmsy_lower <- 0.137
Fmsy_upper <- 0.149
Fmsy_unconstrained <- 0.46
LastIcesAdvice <- 25365
DiscardSurvLastYear <- 0.22256482

## Settings
bioAveYears <- 3
bay <- bioAveYears-2
recyears <- 2002:DataYear

## Run Forecasts for various scenarios
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,1,1),
                               label="SQ all years",
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)  

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,NA,NA),
                               fval=c(NA,NA,Fmsy,Fmsy),
                               label="Fmsy",
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=50000,
                               savesim=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,NA,NA),
                               fval=c(NA,NA,0.000001,0.000001),
                               label="F=0",
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,1,1),
                               fval=c(NA,NA,NA,NA),
                               label=paste0("F=F",DataYear),
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,NA,NA),
                               fval=c(NA,NA,Fpa,Fpa),
                               label="Fpa=Fp05",
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)
set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit,
# 								fscale=c(1,1,NA,NA),
# 								fval=c(NA,NA,1.00,1.00),
# 								label="Flim",
# 								ave.years=(DataYear-bay):(DataYear+1),
# 								rec.years=2002:DataYear,
# 								year.base=DataYear,
# 								splitLD=TRUE,
#								nosim=nruns,
#                               savesim=TRUE)

set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit,
# 								fscale=c(1,1,NA,NA),
# 								nextssb=c(NA,NA,MSYBtrigger,MSYBtrigger),
# 								label=paste0("SSB(",DataYear+3,")=MSYBtrigger"),
# 								ave.years=(DataYear-bay):(DataYear+1),
# 								rec.years=2002:DataYear,
# 								year.base=DataYear,
#								nosim=nruns,
# 								splitLD=TRUE,
#                               savesim=TRUE)

set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit,
# 								fscale=c(1,1,NA,NA),
# 								nextssb=c(NA,NA,13460,13460),
# 								label=paste0("SSB(",DataYear+3,")=Bpa"),
# 								ave.years=(DataYear-bay):(DataYear+1),rec.years=2002:DataYear,
# 								year.base=DataYear,nosim=nruns,
# 								splitLD=TRUE,
#                               savesim=TRUE)

set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit,
# 								fscale=c(1,1,NA,NA),
# 								nextssb=c(NA,NA,Blim,Blim),
# 								label=paste0("SSB(",DataYear+3,")=Blim"),
# 								ave.years=(DataYear-bay):(DataYear+1),
# 								rec.years=2002:DataYear,
# 								year.base=DataYear,nosim=nruns,
# 								splitLD=TRUE,
#                               savesim=TRUE)

set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit,
#                                fscale=c(1,1,NA,NA),
#                                nextssb=c(NA,NA,
#                                          attr(FC[[2]], "tab")[rownames(attr(FC[[2]], "tab")) == (as.character(DataYear+2)), "ssb:median"],
#                                          attr(FC[[2]], "tab")[rownames(attr(FC[[2]], "tab")) == (as.character(DataYear+2)), "ssb:median"]),
#                                label=paste0("SSB(",DataYear+3,")= SSB(", DataYear+2,")"),
#                                ave.years=(DataYear-bay):(DataYear+1),
#                                rec.years=recyears,
#                                year.base=DataYear,
#                                nosim=nruns,
#                                splitLD=TRUE,
#                               savesim=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,NA,NA),
                               fval=c(NA,NA,Fmsy_lower,Fmsy_lower),
                               label="Fmsy (lower)",
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,NA,NA),
                               fval=c(NA,NA,Fmsy_upper,Fmsy_upper),
                               label="Fmsy (upper)",
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit,
                               fscale=c(1,1,NA,NA),
                               catchval=c(NA,NA,
                                          (LastIcesAdvice-(LastIcesAdvice*DiscardSurvLastYear))*1.2,
                                          (LastIcesAdvice-(LastIcesAdvice*DiscardSurvLastYear))*1.2),
                               label=paste0("Catch(", DataYear+2, ")<=20% Advice(", DataYear+1, ")"),
                               ave.years=(DataYear-bay):(DataYear+1),
                               rec.years=recyears,
                               year.base=DataYear,
                               splitLD=TRUE,
                               nosim=nruns,
                               savesim=TRUE)

## Prepare probabilities of falling below Blim
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 results
save(FC, Pabove, file="run/forecast.RData")