## 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:2021			# 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 2020 
#prev.TAC <- 19184 	# 2019 TAC
#prev.TAC <- 10863 		# 2020 TAC
prev.TAC <- 10259		# 2021 TAC
TAC <- 8352

#AdvCat21 <- 5261
AdvCat22 <- 4452



#############################################################################################################################################################################
##  		**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 CM in plenary WGCSE22 5/5/22 from MixFish Advice sheet for 2022)
MixFish <-9240

# 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.571 based on scenario above and initial Fmix 0.582 to hit 9240 in 20212
Fmix<- 0.571


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








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

FC<-list()
      
      set.seed(12345)
      FC[[length(FC)+1]] <- forecast(fit, label = "MSY approach: Fmsy", fval = c(Flast,Fsq,Fmsy,Fmsy), 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,Fsq,Fpa,Fpa), 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,Fsq,Fmsyupper,Fmsyupper), 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,Fsq,Fmsylower,Fmsylower), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)
         set.seed(12345)
         FC[[length(FC)+1]] <- forecast(fit, label = "F = 0", fval = c(Flast,Fsq,0.000001, 0.000001), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)
         set.seed(12345)
         FC[[length(FC)+1]] <- forecast(fit, label = "Fpa", fval = c(Flast,Fsq,Fpa_x,Fpa_x), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)
        # 
        # #set.seed(12345)
        # #FC[[length(FC)+1]] <- forecast(fit, label = "F = Flim", fval = c(Flast,Fsq,Flim,Flim), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)
        # 
        # set.seed(12345)
        # FC[[length(FC)+1]] <- forecast(fit, label = "Fsq", fval = c(Flast,Fsq,Fsq,Fsq), 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,Fsq, NA, NA), nextssb = c(NA, NA, Blim, Blim), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)


#############################################################
#### Needed to get an F for the MixFish catch constraint......run this first!!
##############################################################
#set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label =  "Take MixFish in 2021, then Fmsy", fval = c(Flast, Fsq, Fmsy, Fmsy), ave.years=Ay, rec.years=Ry , splitLD=TRUE, nosim = nosim, #deterministic=deterministic)

#set.seed(12345)
#FC[[length(FC)+1]] <- forecast(fit, label = "New MSY approach: Fmsy", fval = c(Flast,Fsq,Fmsy,Fmsy), 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,Fsq, NA, NA), nextssb = c(NA, NA, Btrigger, Btrigger), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = #nosim,deterministic=deterministic)
      FC[[length(FC)+1]] <- forecast(fit, label =  "F = F2022", fval = c(Flast,Fsq, Fsq, Fsq), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)
      #set.seed(12345)
      #FC[[length(FC)+1]] <- forecast(fit, label =  "F=F2020", fval = c(Flast, Fsq, Fsq, Fsq), ave.years=Ay, rec.years=Ry, splitLD=TRUE, nosim = nosim,deterministic=deterministic)

	 
### Where SSB < MSY Btrigger (39301 t)
###############################################################################################
# MSY SSB in Advice Year (2023)
SSBadv <- 30343
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,Fsq,Fadj,Fmsy), 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,Fsq,Flower,Fmsy), 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,Fsq,Fupper,Fmsy), 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,Fsq, 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")