library(stockassessment)
load("run/model.RData")

# not used
FC<-list()

#for partial mortality F_ibc
oldwd<-setwd("data")
ibf<-read.ices("whg47d_ibf.dat")

#ref points
Fmsy<-0.54  
Fpa<-0.66
Fsq<-0.085

Bpa<-210677
Blim<-150477

rec.years<-c(1983:2023)
base<-max(fit$data$years)-1

# proportion caught in area 4
propc4<-0.8023


# use recent 3 years of fishing selectivity
selYears<-max(fit$data$years) +(-3:-1)  

# to calculate partial Fibc########################
# mean recent 3 years of proportion ibc
ibc_mean<-apply(ibf[length(ibf[,1])+(-2:0),],2,mean) 

#estimated total F at age
summary_F<-summary(fit)
f_est<-summary_F[,7]

#total selectivity
fsel <- faytable(fit)/f_est

# mean recent 3 years of F, scaled t final year
catch_mean<-apply(fsel[length(fsel[,1])+(-3:-1),],2,mean)*f_est[length(f_est)-1]  
IBC_sel<-ibc_mean*catch_mean
Fibc<-mean(IBC_sel[3:6]) #age 2-5



set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1),year.base = base, rec.years=rec.years,overwriteSelYears=selYears,fval=c(Fsq,Fsq,Fsq,Fsq),splitLD =FALSE, label="FSQ, all years")


set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1),year.base = base, rec.years=rec.years,overwriteSelYears=selYears, fval=c(Fsq,Fsq,0.000001,0.000001),splitLD =FALSE, label="SQ then zero")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1),year.base = base, rec.years=rec.years,overwriteSelYears=selYears, fval=c(Fsq,Fsq,Fmsy,Fmsy),splitLD =FALSE, label="SQ then Fmsy")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), year.base = base,rec.years=rec.years,overwriteSelYears=selYears, fval=c(Fsq,Fsq,Fpa,Fpa),splitLD =FALSE, label="SQ then Fpa")


set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), year.base = base,rec.years=rec.years,overwriteSelYears=selYears, fscale=c(NA,NA,NA,NA),splitLD =FALSE, fval=c(Fsq,Fsq,0.75*Fsq,0.75*Fsq), label="SQ then 0.75*SQ")

set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1), year.base = base,rec.years=rec.years,overwriteSelYears=selYears, fscale=c(NA,NA,NA,NA), fval=c(Fsq,Fsq,1.25*Fsq,1.25*Fsq),splitLD =FALSE, label="SQ then 1.25*SQ")


set.seed(12345)
FC[[length(FC)+1]] <- forecast(fit, ave.years=max(fit$data$years) +(-3:-1),year.base = base, rec.years=rec.years,overwriteSelYears=selYears, fval=c(Fsq,Fsq,Fibc,Fibc), customSel=IBC_sel,splitLD =FALSE, label="SQ then no HC fishery") # use partial mortality for IBC of SQ


#set.seed(12345)
#FC[[length(FC)+1]]<-forecast(fit, ave.years=max(fit$data$years) +(-3:-1),rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,NA,NA,NA), #nextssb=c(NA,Blim,Blim,Blim),splitLD =FALSE, label="SQ then Blim in year 3")

#set.seed(12345)
##FC[[length(FC)+1]]<-forecast(fit,ave.years=max(fit$data$years) +(-3:-1),rec.years=2002:max(fit$data$years),overwriteSelYears=selYears, fval=c(Fsq,NA,NA,NA), nextssb=c(NA,Bpa,Bpa,Bpa),splitLD =FALSE, label="SQ then Bpa in year 3")


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