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

Ry<-2020:2024     # latest 5 years as recommended from ADGdeep 2022 (wgdeep2023) 
Fsq<-0.798       # latest value from SAM 
Blim<-9340		# Blim reference point
Bpa<-11627		# Bpa, MSYBtrigger reference point
yearbase<- 2025   # assessment year-1
 # gives exactly Fsq value
set.seed(123456)
 FC[[length(FC)+1]] <- forecast(fit, fval=c(NA,NA ,0.23, 0.23, 0.23, 0.23), fscale=c(NA,NA,NA,NA,NA, NA), catchval=c(6709,5000,NA,NA,NA,NA), label="catch contraint (5000t) in 2026 therafter Fmsy", rec.years=Ry, savesim=TRUE)
set.seed(12345)
FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Fsq,Fsq,Fsq,Fsq,Fsq),rec.years=Ry, label="Rec 2020-2024, SQ all years", processNoise=FALSE,addTSB=TRUE, year.base=yearbase, savesim=TRUE) 
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Fsq,0.000001, 0.000001,0.000001, 0.000001), rec.years=Ry, label="Rec 2020-2024, SQ then zero", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Fsq,0.23, 0.23,0.23, 0.23), rec.years=Ry, label="Rec 2020-2024, SQ then Fmsy", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Fsq,0.60, 0.60,0.60, 0.60), rec.years=Ry, label="Rec 2020-2024, SQ then Fpa/Fp05", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Fsq,0.85, 0.85,0.85, 0.85), rec.years=Ry, label="Rec 2020-2024, SQ then Flim", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Fsq,0.160, 0.160,0.160, 0.160), rec.years=Ry, label="Rec 2020-2024, SQ then Fmsy*SSB2026/MSYBtrigger", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)

# IntermediateF=Fsq then SSB=Blim
set.seed(12345)
FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Fsq,NA,Fsq,Fsq,Fsq),nextssb = c(NA,NA,Blim,NA,NA,NA),label = "Rec 2020-2024, F=Fsq then SSB=Blim", rec.years=Ry, processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)

# IntermediateF=Fsq then SSB=Bpa
#set.seed(12345)
#FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Fsq,NA,Fsq,Fsq,Fsq), nextssb = c(NA,NA,Bpa,NA,NA,NA), label = "Rec 2020-2024, F=Fsq then SSB=Bpa", #rec.years=Ry, processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)

# IntermediateF=Fsq then SSB=SSB2027
ssb2027<- 6517
set.seed(12345)

### Favg NOT DEFINED
#FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Favg, NA, Favg, Favg,Favg), nextssb = c(NA,NA,ssb2027,NA,NA,NA), label = "Rec_F 2020-2024, F=Fsq then SSB=SSB2027", rec.years=Ry, processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)



Ry<-2020:2024     # latest 5 years -1 as recommended from ADGdeep 2022 (wgdeep2023) 
Fsq<-0.798       # latest value from SAM 
#Fsq<-0.607       # NOT latest value from SAM, but the F that gives catches in 2025 of 6708 tons (preliminary catch 6709) 
#Favg<-0.338        # geometric mean of years 2020-2024 latest 5 years -1 
Flow<-0.489        # Flow for Fsq2025
Blim<-9340		# Blim reference point
Bpa<-11627		# Bpa, MSYBtrigger reference point 
yearbase<- 2025   # assessment year-1

# gives exactly Fsq value
set.seed(12345)
FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Flow,Fsq,Fsq,Fsq,Fsq),rec.years=Ry, label="Rec 2020-2024, F=Flow then Fsq*", processNoise=FALSE,addTSB=TRUE, year.base=yearbase, savesim=T) 
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Flow,0.000001, 0.000001,0.000001, 0.000001), rec.years=Ry, label="Rec 2020-2024, F=Flow then zero*", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=T)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Flow,0.23, 0.23,0.23, 0.23), rec.years=Ry, label="Rec 2020-2024, F=Flow then Fmsy*", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=T)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Flow,0.60, 0.60,0.60, 0.60), rec.years=Ry, label="Rec 2020-2024, F=Flow then Fpa/Fp05*", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=T)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Flow,0.85, 0.85,0.85, 0.85), rec.years=Ry, label="Rec 2020-2024, F=Flow then Flim*", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=T)
set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq,Flow,0.160, 0.160,0.160, 0.160), rec.years=Ry, label="Rec 2020-2024, F=Flow then Fmsy*SSB2026/MSYBtrigger*", processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=T)


# IntermediateF=Flow then SSB=Blim
set.seed(12345)
FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Flow, NA, Fsq,Fsq,Fsq), nextssb = c(NA,NA,Blim,NA,NA,NA), label = "Rec 2020-2024, F=Flow then SSB=Blim*", rec.years=Ry, processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)

# IntermediateF=Flow then SSB=Bpa
set.seed(12345)
FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Flow, NA,Fsq, Fsq,Fsq), nextssb = c(NA,NA,Bpa,NA,NA,NA), label = "Rec 2020-2024, F=Flow then SSB=Bpa*", rec.years=Ry, processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)

# IntermediateF=Flow then SSB=SSB2027
ssb2027<- 8033
set.seed(12345)
FC[[length(FC)+1]]  <- forecast(fit, fscale=c(NA,NA,NA,NA,NA,NA), fval=c(Fsq, Flow, NA,Fsq, Fsq,Fsq), nextssb = c(NA,NA,ssb2027,NA,NA,NA), label = "Rec 2020-2024, F=Flow then SSB=SSB2027*", rec.years=Ry, processNoise=FALSE, addTSB=TRUE, year.base=yearbase, savesim=TRUE)




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