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

N.SIM <- 10001  # set the number of iterations the the stochastic forecasts   
# 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 10001 preference

###############################
## Reference points  
#= WKCELTIC 2020
F_MSY <- 0.353        
F_MSY_lower <- 0.221  
F_MSY_upper <- 0.521  
Fpa <- 0.708          
Flim <- 1.40          
Blim <- 9227          
Bpa <- 12822          
MSY_Btrigger <- 12822 

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

prev.advice <- 0
prev.TAC <- 1824

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


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

first.yr <- 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(first.yr,final.yr,1)  # Resample R from previous years, start of the time series to end of time series (median is calculated)
R.range.yrs <-seq(2022,final.yr,1)  # Resample R from previous years, start of the time series to end of time series (median is calculated)

## alternatively sample from the n lowest observed years
n <- 3 # three lowest values (ordinal logistic regression shows 66% likelihood that next recruitment will be in the bottom 4 rank)
#n <- 16-1 #in 2025 this is the 50th percentile, so the median of that is the 25th %ile; take one off to make it an odd number for consistency in median
#n <- 8 #in 2025 this is the 25th percentile (8/32)
#r <- rectable(fit)
#R.range.yrs <- sort(as.numeric(rownames(r[order(r[,1])[1:n],])))


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(
  "F = 0" =            list(fval = c(F_last, F_sq, 0.000001 ,0.000001) ,rec.years=R.range.yrs),
  "FMSYxSSB/Brig" =    list(fval = c(F_last, F_sq, F_MSY_scaled, 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),
  "FMSY" =             list(fval = c(F_last, F_sq, F_MSY, F_MSY) ,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),
#  "Flim" =              list(fval = c(F_last, F_sq, Flim, 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),
#  "TAC"  =             list(fval = c(F_last, F_sq, NA, F_MSY), catchval.exact = c(NA, NA, prev.TAC, NA) ,rec.years=R.range.yrs), #catchval.exact gives error when using large number of sims
  "TAC"  =             list(fval = c(F_last, F_sq, NA, F_MSY), catchval = c(NA, NA, prev.TAC, NA) ,rec.years=R.range.yrs)

# extra options to see effect of fixing recruitment (only resampling from one year)
#  "R2025=R2022 F=0" =  list(fval = c(F_last, F_sq, 0.000001 ,0.000001), rec.years=2022),
#  "R2025=R2022 FMSY" = list(fval = c(F_last, F_sq, F_MSY, F_MSY),rec.years=2022)
 )

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]
                 ,Ry=round(median(unique(FC[[1]][[2]]$rec))) # take the median of the recruitment values, not of the resampled recruitment values. this deals nicely with the "even years issue"
                 ,Ry1=round(median(unique(FC[[1]][[3]]$rec)))
                 ,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,
                 TacChange = 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[i,9] <- 100*attr(FC[[i]],"shorttab")[4,3]/prev.TAC-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,Flim,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")