library(stockassessment)
load("run/model.RData")
source("src/adhocforecast.R")
FC<-list()


set.seed(123456)
rec.ys <- 1989:2024
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,0.28,.28), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="Fmsy, base 2025", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,0.21,.21), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="MSY approach, Fmsy lower", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,0.33,.33), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="MSY approach, Fmsy upper", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,NA,NA), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,0,0), ave.years=max(fit$data$years) + ((-2):(0)), 
label="zero catch", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

#set.seed(123456)
#FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,.49,.49), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
#label="Other options, Flim", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,.35,.35), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="Other options, Fpa", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

### 'manual' calculations 
blim=52076
obj<-function(FF){
  set.seed(123456)
  fc<-adhocforecast(fit, fval=c(NA,NA,FF,FF), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)
  attributes(fc)$shorttab[3,"2028"]
}
FFblim<-uniroot(function(FF)obj(FF)-blim, interval=c(0.001, 2))$root

Bpa=72907
FFbpa<-uniroot(function(FF)obj(FF)-Bpa, interval=c(0.001, 2))$root

obj<-function(FF){
  set.seed(123456)
  fc<-adhocforecast(fit, fval=c(NA,NA,FF,FF), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), label="SSB(2028)=Blim", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)
  attributes(fc)$shorttab[3,"2028"]-attributes(fc)$shorttab[3,"2027"]
}
FF24<-uniroot(obj, interval=c(0.001, 2))$root
                                  
set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,FFblim,FFblim), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), label="SSB(2028)=Blim", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,FFbpa,FFbpa), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), label="SSB(2028)=Bpa", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,FF24,FF24), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), label="SSB(2027)=SSB(2028)", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)


set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,0.230,.28), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="Other options, F=F(2026)", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

###Have to check this one!!!!!!
set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,0.607,0.607), fscale=c(1,NA,NA,NA), catchval=c(NA,29678,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="p(SSB < Blim) = 5%", rec.years=rec.ys, savesim=TRUE, sampleFirstR=FALSE, year.base=2025)

Blim=52076
year=2028
above<-function(x){
  idx<-which(lapply(x, function(xx)xx$year)==year)
  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",year, " > ",Blim,") in %")



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