flukeclimatemse: Run Climate Fluke MSE

# This script contains the operating model function: truestocksize

# To do:
  # update # Equation #Letters with equation numbers from manuscript, can search ??? for most of these equations & table references (input data)
  # make sure R distribution matches input when R turned on (may slightly underestimate recruitment relative to historic, may relate to starting catch conditions used in testing)

#' @title Generate "true" age-structured stock size
#'
#' @description This function projects the operating model recruitment and abundance forward in time and uses these values to calculate stock biomass.
#' @param Catch A vector of catches (by weight) for the current year (t) combined over states for each of 3 fleets in the order: recreational catch, recreational discards, combined commercial catch & discards, no default.
#' @param N_abund A vector of abundance at age (numbers) at the start of the year (t), no default.
#' @param meanR Mean of historical recruitments
#' @param OMR_dev Most recent year's recruitment deviation
#' @param deterministic A boolean defaulting to FALSE, if FALSE biomass and recruitment projections stochastic, if TRUE projections are deterministic & no fishing occurs
#'
#' @return A list containing:
#'  Total biomass in year t
#'  Vector of biomass-at-age in year t
#'  Projected abundance at age in year t+1
#'  Projected recruitment deviation in year t+1
#'  Projected recruitment in year t+1
#'
#' @family operating model functions
#'
#' @examples
#' # Example data: fishing fleet catch & discards, recruitment observations
#' recreational_catch <- read.csv("/Users/ahart2/Downloads/HWK4.csv", skip=1, nrow=Nyear, head=F)[,2] # Recreational catch
#' recreational_discards <- read.csv("/Users/ahart2/Downloads/HWK4.csv", skip=1, nrow=Nyear, head=F)[,2] # Recreational discards
#' commercial_catch-and-discards <- read.csv("/Users/ahart2/Downloads/HWK4.csv", skip=1, nrow=Nyear, head=F)[,2] # Combined commercial catch and discards
#' recruit_observation <- c(81955,	102427,	46954,	78263,	81397,	53988,	12474,	36963,	44019,	47704,	47264,	43928,	58403,	78348,	59520,	52374,	54518,	44100,	60551,	64979,	67860,	50131,	71270,	40634,	48153,	52646,	62460,	73747,	51331,	31296,	35187,	36719,	42271,	29833,	35853,	42415)# Recruitment timeseries (000s) from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A87, page 238, years 1982-2017.
#' # Example
#' truestocksize(Catch = c(recreational_catch[length(recreational_catch)], recreational_discards[length(recreational_discards)], commercial_catch-and-discards[length(commercial_catch-and-discards)]),
#'               N_abund = c(12000, 13000, 12900, 12800, 12500, 12000, 11000, 10000), # fish abundance for age 0-7+
#'               meanR = mean(recruit_observation),
#'               OMR_dev = 0.1) # Recruitment observation error
#'

truestocksize <- function(Catch = NULL,
                          N_abund = NULL,
                          meanR = NULL,
                          OMR_dev = NULL,
                          deterministic = FALSE){
  # Storage
  N_age_tplusone <- rep(NA, 8)
  bio_age <- rep(NA, 8)

  # Fixed/locally calculated inputs
    # Weight at age starting with age 0
    w_age <- c(0.001075, 0.022275, 0.110525, 0.282350, 0.527550, 0.824425, 1.149625, 1.483100) # in Table #H ??? See data_calculations.R for source
    # M at age
    M_age <- c(0.26, 0.26, 0.26, 0.25, 0.25, 0.25, 0.25, 0.24) # Table #H ??? See data_calculations.R for source
    # Recreational catch selectivity at age
    s_rec_cat <- c(0.02333333, 0.21666667, 0.56666667, 0.73333333, 0.87666667, 0.78333333, 0.80000000, 0.75000000) # See data_calculations.R for source
    # Recreational discards selectivity at age
    s_rec_disc <- c(0.1000000, 0.7333333, 0.7000000, 0.5600000, 0.4266667, 0.2833333, 0.2666667, 0.2666667) # See data_calculations.R for source
    # Commercial combined catch and discards selectivity at age
    s_com <- c(0.6483333, 0.5516667, 0.5850000, 0.7666667, 0.4483333, 0.4083333, 0.5550000, 0.8083333) # See data_calculations.R for source
    # Fixed input for marginal variance of recruitment deviations
    R_sigma <- 0.2 # ??? need better estimate
    # Fixed input for magnitude of recruitment autocorrelation, see data_calculations.R for details, in Table #H ???
    R_rho <- 0.414

  # Fleet exploitation rate
    # Recreational catch
    u_rec_cat <- Catch[1]/sum(w_age*s_rec_cat*N_abund*exp(-0.5*M_age)) # Equation #H ???
    # Recreational discards
    u_rec_disc <- Catch[2]/sum(w_age*s_rec_disc*N_abund*exp(-0.5*M_age)) # Equation #H ???
    # Commercial catch & discards
    u_com <- Catch[3]/sum(w_age*s_com*N_abund*exp(-0.5*M_age)) # Equation #H ???

  # Exploitation rate due to all fleets (u_exploit)
  if(deterministic == FALSE){
    u_exploit <- u_rec_cat*s_rec_cat + u_rec_disc*s_rec_disc + u_com*s_com # Equation #G ???
  } else{
    u_exploit <- rep(0, length(N_abund)) # deterministic projection with no fishing
  }

  # Recruitment
    if(deterministic == FALSE){
      delta_t <- rnorm(1, 0, (R_sigma*R_sigma)) # Equation #F ???  # 1 obs, mean=0, variance, change in recruitment
      Rdev_t <- OMR_dev # Recruitment deviation from most recent year (ε_t)
      Rdev_tplusone <- R_rho*Rdev_t + delta_t*((1-R_rho^2)^0.5) # Equation #E ???   # (ε_t+1) projected recruitment deviation
      R_tplusone <- meanR*exp(Rdev_tplusone - 0.5*(R_sigma^2)) # Equation #C ???
    } else{
      R_tplusone <- meanR
    }

  # N_age_tplusone: vector of numbers at age after natural and all fishing mortality
    N_age_tplusone[1] <- R_tplusone # Equation #B
    N_age_tplusone[2:8] <- N_abund[1:7]*exp(-1*M_age[1:7])*(1-u_exploit[1:7]) # Equation #B ???   # Fish age 0-6 in time t become 1-7 time t+1
    N_age_tplusone[8] <- N_age_tplusone[8] + N_abund[8]*exp(-1*M_age[8])*(1-u_exploit[8]) # Equation #B ???

  # Calculate total OM Biomass in year t+1
  OM_Bio <- w_age[1]*R_tplusone + sum(N_age_tplusone[2:8]*w_age[2:8]) # Equation #I ???

  # Store biomass at age
  bio_age[1] <- w_age[1]*R_tplusone # Recruitment # Breakdown of Equation #I ????
  bio_age[2:8] <- N_age_tplusone[2:8]*w_age[2:8] # Age 1-7+ # Breakdown of Equation #I ????

  return(list(OM_Bio = OM_Bio, # Total biomass in year t
              bio_at_age = bio_age, # Vector of biomass-at-age in year t
              OM_Nabundproj = N_age_tplusone, # Projected abundance at age in year t+1
              OM_Rdevproj = Rdev_tplusone, # Projected recruitment deviation in year t+1
              OM_Rproj = R_tplusone)) # Projected recruitment in year t+1
}



# # Uncomment to run examples
#
#
# ###### Real data example
# # Recreational Catch
# TotRecCatDat <- as.matrix(read.csv("/Users/ahart2/Research/Summer_Flounder_MSE/Data/RecLandings.csv"))
# # Recreational Discards
# TotRecDiscDat <- as.matrix(read.csv("/Users/ahart2/Research/Summer_Flounder_MSE/Data/RecDiscards.csv"))
# # Combined commercial catch and discards
# ComCatDiscDat <- as.matrix(read.csv("/Users/ahart2/Research/Summer_Flounder_MSE/Data/CommercialCatDisc.csv"))
#
# # Recruitment timeseries (000s) from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A87, page 238, years 1982-2017.
# RData <- c(81955,	102427,	46954,	78263,	81397,	53988,	12474,	36963,	44019,	47704,	47264,	43928,	58403,	78348,	59520,	52374,	54518,	44100,	60551,	64979,	67860,	50131,	71270,	40634,	48153,	52646,	62460,	73747,	51331,	31296,	35187,	36719,	42271,	29833,	35853,	42415)
# RMean <- mean(RData)
# Rdevs <- c(exp(4.89E-01),	exp(7.12E-01), exp(-6.81E-02),	exp(4.43E-01),	exp(4.82E-01),	exp(7.17E-02),	exp(-1.39E+00),	exp(-3.07E-01),	exp(-1.32E-01),	exp(-5.21E-02),	exp(-6.13E-02),	exp(-1.34E-01),	exp(1.51E-01),	exp(4.45E-01),	exp(1.71E-01),	exp(4.30E-02),	exp(8.36E-02),	exp(-1.28E-01),	exp(1.90E-01),	exp(2.62E-01),	exp(3.06E-01),	exp(4.73E-03),	exp(3.59E-01),	exp(-2.01E-01),	exp(-2.63E-02),	exp(5.39E-02),	exp(2.23E-01),	exp(3.88E-01),	exp(2.63E-02),	exp(-4.70E-01),	exp(-3.55E-01),	exp(-3.12E-01),	exp(-1.79E-01),	exp(-5.34E-01),	exp(-3.68E-01),	exp(-1.83E-01)) # log values of recruitment deviations from assess_cor_file.csv
# # 2018 assessment abundance estimates at age from 66th Northeast Regional Stock Assessment Workshop Assessment Report (NEFSC 2019) Table A89, pg240 2018, years 1982-2017
# AbundData <- read.csv("/Users/ahart2/Research/Summer_Flounder_MSE/Data/SourceFiles/Abund_at_age_A89_NEFSC66AssessmentReport.csv", header=TRUE)
# AbundData <- AbundData[,2:9]
# for(icol in 1:ncol(AbundData)){
#   AbundData[,icol] <- as.numeric(gsub(",","",AbundData[,icol]))
# }
# colnames(AbundData) <- c("Age0", "Age1", "Age2", "Age3", "Age4", "Age5", "Age6", "Age7plus")
#
# tempRdevs <- c(rep(0.2,length(RData))) # placeholder, assessment estimated as normally distributed and summed to 1, maybe need proxy here???
#
# # 1 year projection
# truestocksize(Catch = c(TotRecCatDat[length(TotRecCatDat)], TotRecDiscDat[length(TotRecDiscDat)], ComCatDiscDat[length(ComCatDiscDat)]), # updated
#               N_abund = unlist(AbundData[nrow(AbundData),]), # updated
#               meanR = RMean, # updated
#               OMR_dev = Rdevs[length(Rdevs)]) # log values of recruitment deviations from assess_cor_file.csv
# ### Run deterministic projection to unfished biomass
# Biomass <- NULL
# Biomass_at_age <- NULL
# RDataproj <- NULL
# tempRdevs <- NULL
# for(iyear in 1:50){
#   projResult <- truestocksize(Catch = c(0, 0, 0), # updated
#                               N_abund = unlist(AbundData[nrow(AbundData),]), # updated
#                               meanR = RMean, # updated
#                               OMR_dev = Rdevs[length(Rdevs)],# updated
#                               deterministic = TRUE)
#   AbundData <- rbind(AbundData, projResult$OM_Nabundproj)
#   RDataproj <- rbind(RDataproj, projResult$OM_Rproj)
#   tempRdevs <- rbind(tempRdevs, projResult$OM_Rdevproj) # Projection recruitment timeseries
#   Rdevs <- rbind(Rdevs, projResult$OM_Rdevproj) # Full recruitment timeseries
#   Biomass <- rbind(Biomass, projResult$OM_Bio)
#   Biomass_at_age <- rbind(Biomass_at_age, projResult$bio_at_age)
# }
# rownames(AbundData) <- c()
# Years <- seq(1982, 2017+50)
# plot(y = rowSums(AbundData)[1:47], x = Years[1:47])
# plot(y=rowSums(AbundData), x=Years)
# plot(Biomass) # plot Biomass vs. year
# plot(Biomass_at_age[,7])
# # Plot recruitment
# histogram(RDataproj)
# histogram(RData)