## Forecast code for WGCSE2022
## DS 04/05/2022


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


nosim <-50001 #number of simulations (i think 1001 is defaut, 10001 better, 101 for quick check)
# Note - the number of iterations 10001, up from 1001. All are stable if seed is set - standard "12345" used.
# set to 101 for quick checking
# set to 1001 for routine checking
# set to 50001 preference
deterministic <- F # Leave deterministic = F to default to stochastic forecast


###############################
## Reference points 

Fmsy <- 0.375			# Median point estimate of (F05) EqSim with segreg S/R model
Fmsyupper <- 0.375 	## ICES Advice rule (capped @ Fmsy)
Fmsylower <- 0.315 	## ICES Advice rule

Flim <- 0.64 		# Can't be calculated so using Fmsy instead	## Chk
Fpa <- 0.375 		# Can't be calculated so ICES uses Fmsy	## Chk
Bpa <- 50818		# Btrigger
Blim <- 36571		# Bpa minus erro
Btrigger <- 50818	# Bpa

F05 <- 0.375			## Without Btrigger rule
Fpa <- min(F05, Fpa)
Fpa_x <- 0.375

## Get F status quo
Fstatquo <- tail(fbartable(fit),0)

## Fsq for last 3 yrs 
## FLR version: F_SQ <- mean(fbar(stock)[,as.character(stock@range['maxyear']-2:0)]) # last 3 years - normal default
Flast <- tail(fbartable(fit),1)[,"Estimate"]     # F in the last year of the assessment
Fsq <- mean(tail(fbartable(fit),3)[,"Estimate"]) # Assumed F in the intermediate year
#Fsq <- Flast # scale to the last year


###############################
## Forecast Assumptions

lastyear = max(fit$data$years)	# Final year in assessment, first year forecast (intermediate year)
## Recruitment years (vector covering full survey series) 
Ry <- 2015:2023			# Resample R from previous years (random walk), start of the time series to end of time series (median is calculated)

## biological vars years
# use last 3 years for fishing selectivity and for average weight, maturity, M
Ay <- lastyear +(-2:0) # this is last three years
###############################



###############################
## TAC & advice for 2024 
#prev.TAC <- 19184 	# 2019 TAC
#prev.TAC <- 10863 		# 2020 TAC
#prev.TAC <- 10259		# 2021 TAC
#prev.TAC <- 8352		# 2022 TAC
prev.TAC <- 3877		# 2023 TAC

TAC_24 <- 4718			# 2024 TAC
F_TAC24	<- 0.268		# This is the F needed to hit the TAC in 2024



#AdvCat21 <- 5261
#AdvCat22 <- 4452
AdvCat23 <- 0



#############################################################################################################################################################################
##  		**MixedFish pre-Scenario to get Fmix  for inetrmediate year by setting MixedFish Catch constraint for the intermediate year. Run the Forecast and then check the  F from the Forecast for the Fsq Advice year and update Fmix if required. This may not be necessary now we are using catchval.exact function in the forecast, but keep checking it anyway**
#############################################################################################################################################################################

#  Catch from mixed-fisheries scenario 
# Predicted Whiting catch under Haddock FMSY scenario (updated by from MixFish 2022 during WGCSE23 6/5/23 from MixFish Advice sheet for 2023)
# MixFish <-8859

# set.seed(12345)
# FC[[length(FC)+1]] <- forecast(fit, label =  "Take MixFish in 2021, then Fmsy", catchval.exact = c(Flast, NA, Fmsy, Fmsy), catchval = c(NA, MixFish,NA,NA), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = nosim, nosim = # nosim,deterministic=deterministic)

#Fmix<- 0.582  ## F from MixFish Advice cheet for the current working group year (i.e. intermediate year)
##  Revised to 0.656 based on MSY scenario and initial Fmix 0.57 to hit 8859 above in 2023
# Fmix<- 0.656

###############################################################################################################################################################################



## forecast names and structure part 1
## Forecasts
# main scenarios

FC<-list()
      

         set.seed(12345)
        FC[[length(FC)+1]] <- forecast(fit, label = "F = 0", fval = c(Flast,NA,0.000001, 0.000001), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)


       
#############################################################
set.seed(12345)
      FC[[length(FC)+1]] <- forecast(fit, label = "MSY approach: Fmsy", fval = c(Flast,NA,Fmsy,Fmsy), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)


set.seed(12345)
       FC[[length(FC)+1]] <- forecast(fit, label = "FMSY upper", fval = c(Flast,NA,Fmsyupper,Fmsyupper), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)

set.seed(12345)
       FC[[length(FC)+1]] <- forecast(fit, label = "FMSY lower", fval = c(Flast,NA,Fmsylower,Fmsylower), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, 		splitLD=TRUE, nosim = nosim,deterministic=deterministic)

###############################################################


###############################################################################################
# MSY SSB in Advice Year (2024)
 SSBadv <- 21210
 Fadj <- Fmsy*(SSBadv/Btrigger)
 Flower <- Fmsylower*(SSBadv/Btrigger)
 Fupper <- Fmsyupper*(SSBadv/Btrigger)
###############################################################################################

set.seed(12345)
       FC[[length(FC)+1]] <- forecast(fit, label = "ICES Advice Basis", fval = c(Flast,NA,Fadj,Fmsy), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)

set.seed(12345)
       FC[[length(FC)+1]] <- forecast(fit, label = "ICES Advice lower", fval = c(Flast,NA,Flower,Fmsy), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)

set.seed(12345)
       FC[[length(FC)+1]] <- forecast(fit, label = "ICES Advice upper", fval = c(Flast,NA,Fupper,Fmsy), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)

#################################################################################################

set.seed(12345)
         FC[[length(FC)+1]] <- forecast(fit, label = "Precautionary approach: Fpa", fval = c(Flast,NA,Fpa,Fpa), catchval.exact = c(NA, TAC_24, NA, NA), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim =  nosim,deterministic=deterministic)


#set.seed(12345)
#        FC[[length(FC)+1]] <- forecast(fit, label = "hit Blim", fval = c(Flast,NA, NA, NA), catchval.exact = c(NA, TAC_24, NA, NA), nextssb = c(NA, NA, Blim, Blim), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)


##############################################
### Using advised catch here instead of TAC 
##############################################
#      set.seed(12345)
      #FC[[length(FC)+1]] <- forecast(fit, label =  "hit MSY Btrigger", fval = c(Flast,Fmix, NA, NA), nextssb = c(NA, NA, Btrigger, Btrigger), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = #nosim,deterministic=deterministic)


set.seed(12345)
	FC[[length(FC)+1]] <- forecast(fit, label =  "TAC_2024", fval = c(Flast,NA, NA, NA), catchval.exact = c(NA,TAC_24,TAC_24,TAC_24), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)
    

set.seed(12345)
   	FC[[length(FC)+1]] <- forecast(fit, label =  "F = F2024", fval = c(Flast,F_TAC24, F_TAC24, NA), catchval.exact = c(NA,NA,NA,NA), fscale=c(NA,NA,NA,1), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)


set.seed(12345)
         FC[[length(FC)+1]] <- forecast(fit, label =  "SSB Y+2 = SSB Y+1", fval = c(Flast,NA, NA, NA), catchval.exact = c(NA,TAC_24,NA,NA), nextssb = c(NA,NA,SSBadv,SSBadv), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim, deterministic=deterministic)




#  ICES advice rule for MSY
# MSYappr <- function(fit, Fmsy=NULL, Btrig=NULL, fscale_init=1, fval_init=NA, catchval_init=NA, label=NA){
#   
#   fscale <- fscale_init
#   fval <- fval_init
#   catchval <- catchval_init
#   
#   for (i in 2:4){
#     
#     fscale[i] <- NA
#     fval[i] <- Fmsy
#     catchval[i] <- NA
#     
#     set.seed(12345)
#     f1 <- forecast(fit, fscale=fscale, fval=fval, catchval=catchval, ave.years=Ay, rec.years=Ry, nosim=1001)
#     SSB <- attr(f1,"shorttab")[3,i]
#     ##fval[i] <- min(Fmsy*SSB/Btrig, Fmsy) ## linear scaling between B=0 and B=Btrigger
#     ## linear scaling between F=0 for B=Blim and F=Fmsy for B=Btrigger	
#     a <-  (Fmsy / (Btrig-Blim))
#     b <- -a*Blim
#     if(SSB < Btrig && SSB > Blim) fval[i] <- a*SSB + b
#     
#     if(SSB < Blim) fval[i] <- 0.000001
#   }
#   
#   set.seed(12345)
#   f_final <- forecast(fit, fscale=fscale, fval=fval, catchval=catchval, ave.years=Ay, rec.years=Ry, label=label, nosim= 1001)
#   return(f_final)
# }
# FC[[length(FC)+1]] <- MSYappr(fit, Fmsy=Fmsy, Btrig=Btrigger, fscale_init = NA, fval_init=Fsq, catchval_init = NA, label="F=Fsq then Fmsy HCR")


for(i in 1:length(FC)){
    tmp = FC[[i]]
    maxy =  max( unlist(lapply(tmp,function(xx) xx$year)))
    FC[[i]] = list(  list(year=maxy) )
    attributes(FC[[i]]) <- attributes(tmp)
}
save(FC, file="run/forecast.RData")