## Forecast code for WGCSE2025
## CM 01/05/2026


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

N.SIM <-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  
#= WKCELTIC 2020
F_MSY <- 0.375        #
F_MSY_lower <- 0.315  #
F_MSY_upper <- 0.375  #
Fpa <- 0.375          #
#Flim <- 1.40          #
Blim <- 36571          #
Bpa <- 50818          # 
MSY_Btrigger <- 50818 # 

###############################
## TAC & advice for previous year

prev.advice <- 0
#prev.TAC <- 6353

## catch in last data year
cur.cat <- rev(catchtable(fit,obs.show=T))[1]


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

first.yr <-1999 # min(fit$data$years)  # First year in the assessment (intermediate year) for help estimating intermediate year Recruitment 
final.yr <-max(fit$data$years)  # Final year in the assessment (intermediate year) for help estimating intermediate year Recruitment 
av.yrs <- final.yr+(-2:0)  # final.yr   # Use last three years for mean weights & L-D partition  ##  final.yr+(-2:0) # this is last three years, so five years would be: Ay <- lastyear +(-4:0)
R.range.yrs <-seq(2022,final.yr,1)  # Resample R from previous years, during WGCSE2026 change to the point of decline in recent time series (median is calculated)

F_last <- tail(fbartable(fit),1)[,"Estimate"]  
## Fsq can be mean of last three years:
F_sq <- mean(tail(fbartable(fit),3)[,"Estimate"])
## or for scaled F (use selection pattern for last 3 years but rescale to the Fbar from last year_
#F_sq <- F_last

################################
## initial run to get SSB at the start of the advice year
## include the next line on pc to make the random number generator consistent with stockassessment.org
#RNGkind(sample.kind = "Rounding") 
set.seed(12345)
FC0 <- forecast(fit,fval=c(F_last,F_sq,F_MSY, F_MSY),ave.years=av.yrs, rec.years=R.range.yrs, splitLD=TRUE, nosim = N.SIM)
SSBNext <- attr(FC0,"shorttab")['ssb',3]

F_MSY_scaled <- F_MSY*pmin(1,SSBNext/MSY_Btrigger)
F_MSY_lower_scaled <- F_MSY_lower*pmin(1,SSBNext/MSY_Btrigger)

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

scen <- list(
  "FMSYxSSB/Btig" =    list(fval = c(F_last, F_sq, F_MSY_scaled, F_MSY) ,rec.years=R.range.yrs),
  "FMSY" =             list(fval = c(F_last, F_sq, F_MSY, F_MSY) ,rec.years=R.range.yrs),
  "FlowerxSSB/Btrig" = list(fval = c(F_last, F_sq, F_MSY_lower_scaled, F_MSY) ,rec.years=R.range.yrs),
  "F = 0" =            list(fval = c(F_last, F_sq, 0.000001 ,0.000001) ,rec.years=R.range.yrs),
  "Flower" =           list(fval = c(F_last, F_sq, F_MSY_lower, F_MSY) ,rec.years=R.range.yrs),
  "Fupper" =           list(fval = c(F_last, F_sq, F_MSY_upper, F_MSY) ,rec.years=R.range.yrs),
  "Fpa" =              list(fval = c(F_last, F_sq, Fpa, F_MSY) ,rec.years=R.range.yrs),
  "Blim" =             list(fval = c(F_last, F_sq, NA,F_MSY), nextssb = c(NA, NA, Blim, NA) ,rec.years=R.range.yrs),
  "Bpa, Btrigger" =    list(fval = c(F_last, F_sq, NA, F_MSY), nextssb = c(NA, NA, Bpa, NA) ,rec.years=R.range.yrs),
  "Stable SSB"  =      list(fval = c(F_last, F_sq, NA, F_MSY), nextssb = c(NA, NA, SSBNext, NA) ,rec.years=R.range.yrs),
  "Fsq"  =             list(fval = c(F_last, F_sq, F_sq, F_MSY) ,rec.years=R.range.yrs)
 )

if(SSBNext>MSY_Btrigger) scen[["FMSY x SSB / Blim"]] <- NULL # this scenario is obsolete when ssb>Btrigger
if(SSBNext>MSY_Btrigger) scen[["FMSY lower x SSB / Blim"]] <- NULL # this scenario is obsolete when ssb>Btrigger
if(SSBNext<MSY_Btrigger) scen[["FMSY upper"]] <- NULL # this scenario is not used when ssb<Btrigger
if(F_sq==F_last) scen[["Fsq"]] <- NULL # if they are the same then this is obsolete

FC <- ""
FC <- vector("list", length(scen))
names(FC) <- names(scen)


for(i in seq(scen)){
  set.seed(12345)
  ARGS <- scen[[i]]
  ARGS <- c(ARGS,
            list(fit = fit, ave.years=av.yrs, #rec.years=R.range.yrs, 
            label=names(scen)[i], splitLD=TRUE, nosim = N.SIM,savesim=T))  
  FC[[i]] <- try(do.call(forecast, ARGS))
  if(class(FC[[i]])=='try-error') FC[[i]] <- FC[['F = 0']]
  
  print(paste0("forecast : ", "'", names(scen)[i], "'", " is complete"))
}



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

###########
### some extra tables to get the data formatted for advice sheet

### interrim year assumptions
iy <- data.frame(Fy=F_sq
                 ,SSBy1=SSBNext
                 ,Ry=attr(FC0,"shorttab")['rec',2]
                 ,Ry1=attr(FC0,"shorttab")['rec',3]
                 ,Cy=attr(FC0,"shorttab")['catch',2]
                 ,Ly=attr(FC0,"shorttab")['Land',2]
                 ,Dy=attr(FC0,"shorttab")['Discard',2])
names(iy) <- c(paste0('F',final.yr+1),
  paste0('SSB',final.yr+2),
  paste0('Rec',final.yr+1),
  paste0('Rec',final.yr+2),
  paste0('Cat',final.yr+1),
  paste0('Lan',final.yr+1),
  paste0('Dis',final.yr+1))

### risk to blim

cs <- data.frame(catchTACyr = NA,
                 Landings = NA,
                 Discards = NA,
                 Ftotal = NA,
                 Fland = NA,
                 Fdiscard = NA,
                 ssbTACyrplus1 = NA,
                 ssbChange = NA,
                 adviceChange = NA,
                 Risk = NA)
                 
names(cs)[1:7] <- c(paste0('Cat',final.yr+2),
	paste0('Lan',final.yr+2),
	paste0('Dis',final.yr+2),
    paste0('Ftot',final.yr+2),
    paste0('Flan',final.yr+2),
    paste0('Fdis',final.yr+2),
    paste0('SSB',final.yr+3))

# the 3rd year in the forecast is the advice year because in this stock the last year with data is included in the forecast
for (i in c(1:length(FC))){
  cs[i,1] <- attr(FC[[i]],"shorttab")[4,3]
  cs[i,2] <- attr(FC[[i]],"shorttab")[7,3]
  cs[i,3] <- attr(FC[[i]],"shorttab")[8,3]
  cs[i,4] <- attr(FC[[i]],"shorttab")[1,3]
  cs[i,5] <- attr(FC[[i]],"shorttab")[5,3]
  cs[i,6] <- attr(FC[[i]],"shorttab")[6,3]
  cs[i,7] <- attr(FC[[i]],"shorttab")[3,4]
  cs[i,8] <- 100*attr(FC[[i]],"shorttab")[3,4]/SSBNext-100
  cs[i,9] <- 100*attr(FC[[i]],"shorttab")[4,3]/prev.advice-100
}
cs[,10] <- 100*unlist(lapply(FC, function(x) sum(x[[4]]$ssb < Blim) / length(x[[4]]$ssb)))
rownames(cs) <- names(scen)
save(FC,iy,cs, file="run/forecast.RData")

# reference points (to make sure you are using the correct ones)
rp <- data.frame(Btrigger=MSY_Btrigger,Fmsy=F_MSY,Blim,Bpa=MSY_Btrigger,Fpa,Btrigger=MSY_Btrigger,Bpa,Fmsy=F_MSY,Flower=F_MSY_lower,Fupper=F_MSY_upper)
save(FC,iy,cs,rp, file="run/forecast.RData")