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


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

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,38181,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="MSY approach, Fmsy lower", rec.years=rec.ys, savesim=TRUE, sampleFirstR=TRUE, year.base=2024)

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,38181,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="MSY approach, Fmsy upper", rec.years=rec.ys, savesim=TRUE, sampleFirstR=TRUE, year.base=2024)

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

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

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

### '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,38181,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), rec.years=rec.ys, savesim=TRUE, sampleFirstR=TRUE, year.base=2024)
  attributes(fc)$shorttab[3,"2027"]
}
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,38181,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), label="SSB(2027)=Blim", rec.years=rec.ys, savesim=TRUE, sampleFirstR=TRUE, year.base=2024)
  attributes(fc)$shorttab[3,"2027"]-attributes(fc)$shorttab[3,"2026"]
}
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,38181,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), label="SSB(2027)=Blim", rec.years=rec.ys, savesim=TRUE, sampleFirstR=TRUE, year.base=2024)

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

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


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

set.seed(123456)
FC[[length(FC)+1]] <- adhocforecast(fit, fval=c(NA,NA,0.523,0.523), fscale=c(1,NA,NA,NA), catchval=c(NA,38181,NA,NA), ave.years=max(fit$data$years) + ((-2):(0)), 
label="p(SSB < Blim) = 5%", rec.years=rec.ys, savesim=TRUE, sampleFirstR=TRUE, year.base=2024)

Blim=52076
year=2027
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")