msy: Estimation of Equilibrium Reference Points for Fisheries

Documented in eqsim_run

#' Simulates the Equilibrium Results for a Population.
#'
#' Simulate a fish stock forward in time given biological parameters, fishery
#' parameters and advice parameters.
#'
#' @param fit A list returned from the function fitModels
#' @param bio.years The years to sample maturity, weights and M from, given as
#'                  a vector of length 2, i.e. c(2010, 2015) select from the
#'                  years 2010 to 2015 inclusive.
#' @param bio.const A flag (default FALSE), if TRUE mean of the biological values from the
#'                  years selected are used
#' @param sel.years The years to sample the selection patterns from, given as
#'                  a vector of length 2, i.e. c(2010, 2015) select from the
#'                  years 2010 to 2015 inclusive.
#' @param sel.const A flag (default FALSE), if TRUE mean of the selection patterns from the
#'                  years selected are used
#' @param Fscan F values to scan over, i.e. seq(0, 2, by = 0.05)
#' @param Fcv Assessment error in the advisory year
#' @param Fphi Autocorrelation in assessment error in the advisory year
#' @param SSBcv Spawning stock biomass error in the advisory year
#' @param rhologRec A flag for recruitment autocorrelation, default (TRUE), or a
#'                  vector of numeric values specifcying the autocorrelation
#'                  parameter for the residuals for each SR model.
#' @param Blim SSB limit reference point
#' @param Bpa SSB precuationary reference point
#' @param recruitment.trim A numeric vector with two log-value clipping the
#'        extreme recruitment values from a continuous lognormal distribution.
#'        The values must be set as c("high","low").
#' @param Btrigger If other than 0 (default) the target F applied is reduced by
#'                 SSB/Btrigger. This is the "ICES Advice Rule".
#' @param Nrun The number of years to run in total (the last 50 years from that
#'             will be retained to compute equilibrium values from)
#' @param process.error Use stochastic recruitment or mean recruitment?
#'                      TRUE (default) uses the predictive distribution of recruitment,
#'                      model estimate of recruitment + simulated observation
#'                      error.  FALSE uses model prediction of recruitment with
#'                      no observation error.
#' @param verbose Flag, if TRUE (default) indication of the progress of the
#'        simulation is provided in the console. Useful to turn to FALSE when
#'        knitting documents.
#' @param extreme.trim a pair of quantiles (low, high) which are used to trim
#'                     the equilibrium catch values, across simulations within
#'                     an F scenario, when calculating the mean catch and
#'                     landings for that F scenario.  These mean values
#'                     calculated accross simulations within an F scenario
#'                     are used to find which F scenario gave the maximum catch.
#'                     \code{extreme.trim} can therefore be used to stablise the
#'                     estimate of mean equilibrium catch and landings by F
#'                     scenario.  The default is c(0, 1) which includes all the
#'                     data and is effectively an untrimmed mean.
#' @param R.initial Initial recruitment for the simulations.  This is common
#'                  accross all simulations. Default = mean of all recruitments
#'                  in the series.
#' @param keep.sims Flag, if TRUE returns a matrix of population tragectories
#'                  for each value of F in Fscan (see examples).
#' @return
#' A list containing the results from the forward simulation and the reference
#' points calculated from it.
#'
#' @details
#' Details of the steps required to evaluate reference points are given in
#' ICES (2017).  WHile, details of the calculation of MSY ranges is given in
#' ICES (2015).
#'
#' @references
#' ICES (2015) Report of the Workshop to consider F MSY ranges for stocks in
#' ICES categories 1 and 2 in Western Waters (WKMSYREF4).
#' \href{http://ices.dk/sites/pub/Publication\%20Reports/Expert\%20Group\%20Report/acom/2015/WKMSYREF4/01\%20WKMSYREF4\%20Report.pdf}{01
#' WKMSYREF4 Report.pdf}
#'
#' ICES (2017) ICES fisheries management reference points for category 1 and 2
#' stocks.
#' DOI: \href{https://doi.org/10.17895/ices.pub.3036}{10.17895/ices.pub.3036}
#'
#' @seealso
#'
#' \code{\link{eqsr_fit}} fits multiple stock recruitment models to a data set.
#'
#' \code{\link{eqsr_plot}} plots the results from eqsr_fit.
#'
#' \code{\link{eqsim_plot}} summary plot of the forward simulation showing estimates
#'   of various reference points.
#'
#' \code{\link{eqsim_plot_range}} summary plots of the forward simulation showing
#'   the estimates of MSY ranges (ICES, 2015)
#'
#' \code{\link{msy-package}} gives an overview of the package.
#'
#' @examples
#' \dontrun{
#'
#' data(icesStocks)
#'
#' FIT <- eqsr_fit(icesStocks$saiNS,
#'   nsamp = 1000,
#'   models = c("Ricker", "Segreg")
#' )
#' SIM <-
#'   eqsim_run(
#'     FIT,
#'     bio.years = c(2004, 2013),
#'     sel.years = c(2004, 2013),
#'     Fcv = 0.24,
#'     Fphi = 0.42,
#'     Blim = 106000,
#'     Bpa = 200000,
#'     Fscan = seq(0, 1.2, len = 40)
#'   )
#'
#' # extract tragectories
#' ssbsim <- SIM$rbya$ssb
#' years <- SIM$rbya$simyears
#' models <- SIM$rbya$srmodels$model
#' Ftarget <- SIM$rbya$Ftarget
#'
#' Fval <- which(Ftarget == 0)
#' Fval <- which(Ftarget > .3)[1]
#' x <- ssbsim[Fval, , ]
#' df <- data.frame(
#'   year = 1:nrow(x),
#'   ssb = c(x),
#'   sim = rep(1:ncol(x), each = nrow(x)),
#'   model = rep(models, each = nrow(x))
#' )
#' xyplot(ssb ~ year | model, groups = sim, data = df, type = "l", col = grey(0.5, alpha = 0.5))
#'
#' fit <- density(x[x > 1e-3], from = 0)
#' plot(fit$x, fit$y * mean(x > 1e-3), col = "red", type = "l")
#' lines(x = 0, y = mean(x <= 1e-3), type = "h", lwd = 3)
#' }
#'
#' @importFrom TAF msg
#' @importFrom stats rnorm loess predict approx cor
#' @importFrom stats quantile median density complete.cases fitted dnorm
#'
#' @importFrom FLCore window catch.n catch.wt harvest
#' @importFrom FLCore harvest.spwn landings.n landings.wt m m.spwn mat stock.wt
#'
#' @export
eqsim_run <- function(fit,
                      bio.years = c(-5, -1) + dims(fit$stk)$maxyear, # years sample weights, M and mat
                      bio.const = FALSE,
                      sel.years = c(-5, -1) + dims(fit$stk)$maxyear, # years sample sel and discard proportion by number from
                      sel.const = FALSE,
                      Fscan = seq(0, 2, len = 40), # F values to scan over
                      Fcv = 0,
                      Fphi = 0,
                      SSBcv = 0,
                      rhologRec = TRUE,
                      Blim,
                      Bpa,
                      recruitment.trim = c(3, -3),
                      Btrigger = 0,
                      Nrun = 200, # number of years to run in total
                      process.error = TRUE, # use predictive recruitment or mean recruitment? (TRUE = predictive)
                      verbose = TRUE,
                      extreme.trim = c(0, 1),
                      R.initial = mean(fit$rby$rec),
                      keep.sims = FALSE) {
  if (abs(Fphi) >= 1) stop("Fphi, the autocorelation parameter for log F should be between (-1, 1)")
  if (diff(recruitment.trim) > 0) stop("recruitment truncation must be given as c(high, low)")
  # commented out as above line is a better check
  # if ((recruitment.trim[1] + recruitment.trim[2]) > 0) stop("recruitment truncation must be between a high - low range")

  if (verbose) msg("Setting up...")

  if (length(bio.years) > 2) {
    stop("bio.years must be given as a length two vector: c(first, last)")
  }
  if (length(sel.years) > 2) {
    stop("sel.years must be given as a length two vector: c(first, last)")
  }

  btyr1 <- bio.years[1]
  btyr2 <- bio.years[2]
  slyr1 <- sel.years[1]
  slyr2 <- sel.years[2]
  # Keep at most 50 simulation years (which will be the last 50 of the Nrun
  #  forward simulated years)
  keep <- min(Nrun, 50)

  SR <- fit$ sr.sto
  data <- fit$ rby[, c("rec", "ssb", "year")]
  stk <- fit$ stk

  # forecast settings (mean wt etc)
  stk.win <- window(stk, start = btyr1, end = btyr2)
  stk.winsel <- window(stk, start = slyr1, end = slyr2)

  littleHelper <- function(x, i) {
    x2 <- x
    x2[i] <- NA
    x2[] <- apply(x2, 1, mean, na.rm = TRUE)
    x[i] <- x2[i]
    return(x)
  }

  west <- matrix(stock.wt(stk.win), ncol = btyr2 - btyr1 + 1)
  i <- west == 0
  if (any(i)) west <- littleHelper(west, i)
  weca <- matrix(catch.wt(stk.win), ncol = btyr2 - btyr1 + 1)
  i <- weca == 0
  if (any(i)) weca <- littleHelper(weca, i)
  wela <- matrix(landings.wt(stk.win), ncol = btyr2 - btyr1 + 1)
  if (any(i)) wela <- littleHelper(wela, i)

  Mat <- matrix(mat(stk.win), ncol = btyr2 - btyr1 + 1)
  M <- matrix(m(stk.win), ncol = btyr2 - btyr1 + 1)
  landings <- matrix(landings.n(stk.winsel), ncol = slyr2 - slyr1 + 1)
  # if zero, use 0.10 of minimum value

  catch <- matrix(catch.n(stk.winsel), ncol = slyr2 - slyr1 + 1)
  sel <- matrix(harvest(stk.winsel), ncol = slyr2 - slyr1 + 1)
  Fbar <- matrix(fbar(stk.winsel), ncol = slyr2 - slyr1 + 1)
  sel <- sweep(sel, 2, Fbar, "/")

  if (sel.const == TRUE) { # take means of selection
    sel[] <- apply(sel, 1, mean)
    landings[] <- apply(landings, 1, mean)
    catch[] <- apply(catch, 1, mean)
  }

  # 22.2.2014 Added weight of landings per comment from Carmen
  if (bio.const == TRUE) { # take means of wts Mat and M and ratio of landings to catch
    west[] <- apply(west, 1, mean)
    weca[] <- apply(weca, 1, mean)
    wela[] <- apply(wela, 1, mean)
    Mat[] <- apply(Mat, 1, mean)
    M[] <- apply(M, 1, mean) # me
  }
  land.cat <- landings / catch # ratio of number of landings to catch

  # TODO: Check if this is sensible
  i <- is.na(land.cat)
  if (any(i)) land.cat[i] <- 1

  Fprop <- apply(harvest.spwn(stk.winsel), 1, mean)[drop = TRUE] # vmean(harvest.spwn(stk.win))
  Mprop <- apply(m.spwn(stk.win), 1, mean)[drop = TRUE] # mean(m.spwn(stk.win))

  # get ready for the simulations
  Nmod <- nrow(SR)
  NF <- length(Fscan)
  ages <- dims(stk)$age
  ssb_lag <- fit$rby$ssb_lag[1]

  ssby <- Ferr <- array(0, c(Nrun, Nmod), dimnames = list(year = 1:Nrun, iter = 1:Nmod))
  Ny <- Fy <- WSy <- WCy <- Cy <- Wy <- Wl <- Ry <-
    array(0, c(ages, Nrun, Nmod),
      dimnames = list(
        age = (range(stk)[1]:range(stk)[2]),
        year = 1:Nrun,
        iter = 1:Nmod
      )
    )
  # TODO per note from Carmen:
  #  NOTE: If we want Ferr to be a stationary AR(1) process, it would make
  #        more sense to initialise Ferr as a Normal dist with zero mean and
  #        standard deviation of AR(1) marginal distribution, i.e. standard
  #        deviation of initial Ferr = Fcv/sqrt(1- Fphi^2), instead of just
  #        initialising Ferr=0
  #  2014-03-12: Changed per note form Carmen/John
  Ferr[1, ] <-    rnorm(n = Nmod, mean = 0, sd = 1) * Fcv / sqrt(1 - Fphi^2)
  for (j in 2:Nrun) {
    Ferr[j, ] <- Fphi * Ferr[j - 1, ] + Fcv *    rnorm(n = Nmod, mean = 0, sd = 1)
  }

  # 2014-03-12: Changed per note form Carmen/John
  #  Errors in SSB: this is used when the ICES MSY HCR is applied for F
  SSBerr <- matrix(   rnorm(n = Nrun * Nmod, mean = 0, sd = 1), ncol = Nmod) * SSBcv

  rsam <- array(sample(1:ncol(weca), Nrun * Nmod, TRUE), c(Nrun, Nmod))
  rsamsel <- array(sample(1:ncol(sel), Nrun * Nmod, TRUE), c(Nrun, Nmod))
  Wy[] <- c(weca[, c(rsam)])
  Wl[] <- c(wela[, c(rsam)])
  Ry[] <- c(land.cat[, c(rsamsel)])

  # initial recruitment
  R <- R.initial

  # set up arrays to contain simulations
  ssbs <- cats <- lans <- recs <- array(0, c(7, NF))
  ferr <- ssbsa <- catsa <- lansa <- recsa <- array(0, c(NF, keep, Nmod))
  if (keep.sims) {
    # keep all ssbs
    ssbsall <- catsall <- lansall <- recsall <- array(0, c(NF, Nrun, Nmod))
  }
  begin <- Nrun - keep + 1

  # New from Simmonds' 29.1.2014
  #   Residuals of SR fits (1 value per SR fit and per simulation year
  #     but the same residual value for all Fscan values):
  resids <- array(   rnorm(Nmod * (Nrun + 1), 0, SR$cv), c(Nmod, Nrun + 1))

  # 2014-03-12: Changed per note form Carmen/John
  #  Autocorrelation in Recruitment Residuals:
  if (rhologRec == TRUE) {
    fittedlogRec <- do.call(cbind, lapply(c(1:nrow(fit$sr.sto)), function(i) {
      FUN <- match.fun(fit$sr.sto$model[i])
      FUN(fit$sr.sto[i, ], fit$rby$ssb)
    }))
    # Calculate lag 1 autocorrelation of residuals:
    rhologRec <- apply(log(fit$rby$rec) - fittedlogRec, 2, function(x) {
         cor(x[-length(x)], x[-1])
    })
  }
  if (is.numeric(rhologRec)) {
    # Draw residuals according to AR(1) process:
    for (j in 2:(Nrun + 1)) {
      resids[, j] <- rhologRec * resids[, j - 1] + resids[, j] * sqrt(1 - rhologRec^2)
    }
  }


  # Limit how extreme the Rec residuals can get:
  lims <- t(array(SR$cv, c(Nmod, 2))) * recruitment.trim
  for (k in 1:Nmod) {
    resids[k, resids[k, ] > lims[1, k]] <- lims[1, k]
  }
  for (k in 1:Nmod) {
    resids[k, resids[k, ] < lims[2, k]] <- lims[2, k]
  }
  # end New from Simmonds 29.1.2014

  if (verbose) msg("Running forward simulations.")
  if (verbose) loader(0)

  # Looping over each F value in Fscan. For each of the Nmod SR fits
  # (replicates), do a forward simulation during Nrun years
  # There are Rec residuals for each SR fit and year, which take the same
  # values for all Fscan
  for (i in 1:NF) {
    # The F value to test
    Fbar <- Fscan[i]

    ############################################################################
    # Population in simulation year 1:

    # Zpre: Z that occurs before spawning
    Zpre <- Fbar * sel[, rsamsel[1, ]] * Fprop + M[, rsam[1, ]] * Mprop

    # Zpos: Z that occurs after spawning
    # Zpos not used anywhere
    Zpos <- Fbar * (1 - Fprop) * sel[, rsamsel[1, ]] + M[, rsam[1, ]] * (1 - Mprop)

    # run Z out to age 50 for plus group...
    # TODO:
    # Comments from Carmen: Zcum is a cumulative sum:
    #  There is a matrix of F-at-age and a matrix of M-at-age (each has 49 ages, Nmod replicates)
    #  The F and M matrices are summed, giving Z-at-age (49 ages, Nmod replicates)
    Ztot <- Fbar * sel[c(1:ages, rep(ages, 49 - ages)), rsamsel[1, ]] + M[c(1:ages, rep(ages, 49 - ages)), rsam[1, ]]
    Zcum <- apply(Ztot, 2, function(x) c(0, cumsum(x)))
    # create initial population structure
    N1 <- R * exp(-unname(Zcum))

    # set up age structure out to age 50 in first years for all simulations
    Ny[, 1, ] <- rbind(N1[1:(ages - 1), ], colSums(N1[ages:50, ]))

    # calculate ssb in first year using a different stock.wt and Mat selection and M for each simulation
    ssby[1, ] <- colSums(Mat[, rsam[1, ]] * Ny[, 1, ] * west[, rsam[1, ]] / exp(Zpre))

    # if rec recruiting year class comes from previous years ssb, as in fish recruiting
    # at age 2 or winter ring herring ageing then run some more initial years
    # using the same intial population
    # - NOTE roll forward one year incase ssb_lag is 0 so that we always have a year j-1.
    for (j in 2:pmax(2, ssb_lag + 2)) {
      Ny[, j, ] <- rbind(N1[1:(ages - 1), ], colSums(N1[ages:50, ]))
      ssby[j, ] <- colSums(Mat[, rsam[j - 1, ]] * Ny[, 1, ] * west[, rsam[j - 1, ]] / exp(Zpre))
    }

    # Years (2 + ssb_lag) to Nrun:
    for (j in (2 + ssb_lag):Nrun) {
      #  year j is the projection year,
      #  year j-1 is where fishing is going to take place
      #  conceptually the assessment takes place in year j-2

      # apply HCR
      # (intended) Fbar to be applied in year j-1 (depends on SSB in year j-1):
      # 2014-03-12: Changed per note form Carmen/John
      # Fnext <- Fbar * pmin(1, SSB/Btrigger)
      Fnext <- Fbar * pmin(1, ssby[j - 1, ] * exp(SSBerr[j - 1, ]) / Btrigger)

      # apply some noise to the F
      # Notes from Carmen:
      #  Assessment and/or implementation error (modifies intended F to get
      #  realised F)
      #  Error: AR(1) process on log(F) with autocorrel = Fphi, and
      #  conditional stand deviation = Fcv
      #  Might make more sense to have the "Ferr" matrix calculated before
      #  the Fscan loop starts so that the same errors in F are applied to
      #  all Fscan values ???? (as for Rec residuals)

      # Outcommented 2014-03-12 because F-error already been drawn outside the
      #   loop, so this line here is no longer needed:
      # Ferr[j,] <- Fphi * Ferr[j-1,] + rnorm(Nmod, 0, Fcv)

      # realised Fbar in year j-1:
      Fnext <- exp(Ferr[j, ]) * Fnext

      # get a selection pattern for each simulation and apply this to get N
      Zpre <- rep(Fnext, each = length(Fprop)) * Fprop * sel[, rsamsel[j - 1, ]] + M[, rsam[j - 1, ]] * Mprop

      # get Fy
      Fy[, j - 1, ] <- rep(Fnext, each = ages) * sel[, rsamsel[j - 1, ]]

      # roll population one year forward having decided in the F value
      Ny[-1, j, ] <- Ny[1:(ages - 1), j - 1, ] * exp(-Fy[1:(ages - 1), j - 1, ] - M[1:(ages - 1), rsam[j - 1, ]])
      # calculate plus group
      Ny[ages, j, ] <- Ny[ages, j, ] + Ny[ages, j - 1, ] * exp(-Fy[ages, j - 1, ] - M[ages, rsam[j - 1, ]])

      if (ssb_lag == 0) {
        # calculate ssb ignores contribution of recruiting age 0 fish
        ssby[j, ] <- apply(array(Mat[, rsam[j, ]] * Ny[, j, ] * west[, rsam[j, ]] / exp(Zpre), c(ages, Nmod)), 2, sum)
      }

      # simulate recruitment in year j
      # get ssb from appropriate year, if ssb_lag is zero, then current year ssb is used
      SSBforRec <- ssby[j - ssb_lag, ]

      # predict recruitment using various models
      if (process.error) {
        # Changes 29.1.2014
        # new random draws each time
        # allrecs <- sapply(unique(SR $ mod), function(mod) exp(match.fun(mod) (SR, SSB) + rnorm(Nmod, 0, SR $ cv)))
        # same random draws used for each F
        ###### 2014-03-13  TMP COMMENT - ERROR OCCURS HERE
        allrecs <- sapply(unique(SR$mod), function(mod) exp(match.fun(mod)(SR, SSBforRec) + resids[, j]))
        # end Changes 29.1.2014
      } else {
        allrecs <- sapply(unique(SR$mod), function(mod) exp(match.fun(mod)(SR, SSBforRec)))
      }

      # Comment from Carmen:
      #  For each of the Nmod replicates, this selects the appropriate SR model
      #   type to use in that replicate
      #  Note that the order of SR model types that comes out in "select" is
      #   not necessarily the same order in which the SR model types were
      #   entered as inputs -- I presume the **next 2 lines** of code have
      #   been checked to avoid potential bugs due to this reordering  ????
      select <- cbind(seq(Nmod), as.numeric(factor(SR$mod, levels = unique(SR$mod))))
      Ny[1, j, ] <- allrecs[select]

      # calculate ssb now we have the recruiting age class abundance
      ssby[j, ] <- apply(array(Mat[, rsam[j, ]] * Ny[, j, ] * west[, rsam[j, ]] / exp(Zpre), c(ages, Nmod)), 2, sum)

      # calculate catch.n (should this be j-1?  does it matter?)
      Cy[, j, ] <- Ny[, j - 1, ] * Fy[, j - 1, ] / (Fy[, j - 1, ] + M[, rsam[j - 1, ]]) * (1 - exp(-Fy[, j - 1, ] - M[, rsam[j - 1, ]]))
    }

    # convert to catch weight
    Cw <- Cy * Wy # catch Numbers *catch wts
    land <- Cy * Ry * Wl # catch Numbers * Fraction (in number) landed and landed wts
    Lan <- apply(land, 2:3, sum)
    Cat <- apply(Cw, 2:3, sum)

    # summarise everything and spit out!
    ferr[i, , ] <- Ferr[begin:Nrun, ]
    ssbsa[i, , ] <- ssby[begin:Nrun, ]
    catsa[i, , ] <- Cat[begin:Nrun, ]
    lansa[i, , ] <- Lan[begin:Nrun, ]
    recsa[i, , ] <- Ny[1, begin:Nrun, ]

    # store quantiles
    quants <- c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975)
    ssbs[, i] <-    quantile(ssbsa[i, , ], quants)
    cats[, i] <-    quantile(catsa[i, , ], quants)
    lans[, i] <-    quantile(lansa[i, , ], quants)
    recs[, i] <-    quantile(recsa[i, , ], quants)

    # if user has requested full simulations
    if (keep.sims) {
      ssbsall[i, , ] <- ssby
      catsall[i, , ] <- Cat
      lansall[i, , ] <- Lan
      recsall[i, , ] <- Ny[1, , ]
    }

    if (verbose) loader(i / NF)
  }

  if (verbose) msg("Summarising simulations")

  dimnames(ssbs) <- dimnames(cats) <-
    dimnames(lans) <- dimnames(recs) <-
    list(
      quants = c("p025", "p05", "p25", "p50", "p75", "p95", "p975"),
      fmort = Fscan
    )

  rbp2dataframe <- function(x, variable) {
    x <- data.frame(t(x))
    x$variable <- variable
    x$Ftarget <- as.numeric(row.names(x))
    rownames(x) <- NULL
    return(x)
  }
  rbp <- rbind(
    rbp2dataframe(recs, "Recruitment"),
    rbp2dataframe(ssbs, "Spawning stock biomass"),
    rbp2dataframe(cats, "Catch"),
    rbp2dataframe(lans, "Landings")
  )
  rbp <- rbp[, c(9, 8, 1:7)]

  # STOCK REFERENCE POINTS

  FCrash05 <- Fscan[which.max(cats[2, ]):NF][which(cats[2, which.max(cats[2, ]):NF] < 0.05 * max(cats[2, ]))[1]]
  FCrash50 <- Fscan[which.max(cats[4, ]):NF][which(cats[4, which.max(cats[4, ]):NF] < 0.05 * max(cats[4, ]))[1]]


  # Einar amended 30.1.2014
  if (missing(extreme.trim)) {
    catm <- apply(catsa, 1, mean)
    lanm <- apply(lansa, 1, mean)
  } else {
    # 2014-03-12 Outcommented per note from Carmen/John - see below
    # x <- catsa
    # i <- x > quantile(x,extreme.trim[2]) |
    #  x < quantile(x,extreme.trim[1])
    # x[i] <- NA
    # catm <- apply(x, 1, mean, na.rm=TRUE)
    #
    # x <- lansa
    # i <- x > quantile(x,extreme.trim[2]) |
    #  x < quantile(x,extreme.trim[1])
    # x[i] <- NA
    # lanm <- apply(x, 1, mean, na.rm=TRUE)

    # 2014-03-12: Above replaced with the following per note from Carmen/John
    #  If we want to remove whole SR models, we could use the following code. But it is too extreme, it ends up getting rid of most models:
    # auxi2 <- array( apply(catsa, 1, function(x){auxi<-rep(TRUE,Nmod); auxi[x > quantile(x, extreme.trim[2]) | x < quantile(x, extreme.trim[1])] <- FALSE; x <- auxi } ), dim=c(keep,Nmod,NF))
    # auxi2 <- (1:Nmod)[apply(auxi2, 2, function(x){length(unique(as.vector(x)))})==1]
    # apply(catsa[,,auxi2],1,mean)

    # So I think the alternative is not to get rid of whole SR models, but of different SR models depending on the value of F:
    catm <- apply(catsa, 1, function(x) {
      mean(x[x <=    quantile(x, extreme.trim[2]) & x >=    quantile(x, extreme.trim[1])])
    })
    lanm <- apply(lansa, 1, function(x) {
      mean(x[x <=    quantile(x, extreme.trim[2]) & x >=    quantile(x, extreme.trim[1])])
    })
  }

  # end Einar amended 30.1.2014

  maxcatm <- which.max(catm)
  maxlanm <- which.max(lanm)

  # Einar added 29.1.2014
  rbp$Mean <- NA
  rbp$Mean[rbp$variable == "Catch"] <- catm
  rbp$Mean[rbp$variable == "Landings"] <- lanm
  # end Einar added 29.1.2014


  catsam <- apply(catsa, c(1, 3), mean)
  lansam <- apply(lansa, c(1, 3), mean)
  maxpf <- apply(catsam, 2, which.max)
  maxpfl <- apply(lansam, 2, which.max)

  FmsyLan <- Fscan[maxpfl]
  msymLan <- mean(FmsyLan)
  vcumLan <-    median(FmsyLan)
  fmsy.densLan <-    density(FmsyLan)
  vmodeLan <- fmsy.densLan$x[which.max(fmsy.densLan$y)]

  FmsyCat <- Fscan[maxpf]
  msymCat <- mean(FmsyCat)
  vcumCat <-    median(FmsyCat)
  fmsy.densCat <-    density(FmsyCat)
  vmodeCat <- fmsy.densCat$x[which.max(fmsy.densCat$y)]

  pFmsyCat <- data.frame(
    Ftarget = fmsy.densCat$x,
    value = cumsum(fmsy.densCat$y * diff(fmsy.densCat$x)[1]),
    variable = "pFmsyCatch"
  )
  pFmsyLan <- data.frame(
    Ftarget = fmsy.densLan$x,
    value = cumsum(fmsy.densLan$y * diff(fmsy.densLan$x)[1]),
    variable = "pFmsyLandings"
  )
  pProfile <- rbind(pFmsyCat, pFmsyLan)

  # PA REFERENCE POINTS
  if (!missing(Blim)) {
    pBlim <- apply(ssbsa > Blim, 1, mean)

    i <- max(which(pBlim > .95))
    grad <- diff(Fscan[i + 0:1]) / diff(pBlim[i + 0:1])
    flim <- Fscan[i] + grad * (0.95 - pBlim[i]) # linear interpolation i think..

    i <- max(which(pBlim > .90))
    grad <- diff(Fscan[i + 0:1]) / diff(pBlim[i + 0:1])
    flim10 <- Fscan[i] + grad * (0.9 - pBlim[i]) # linear interpolation i think..

    i <- max(which(pBlim > .50))
    grad <- diff(Fscan[i + 0:1]) / diff(pBlim[i + 0:1])
    flim50 <- Fscan[i] + grad * (0.5 - pBlim[i]) # linear interpolation i think..

    pBlim <- data.frame(Ftarget = Fscan, value = 1 - pBlim, variable = "Blim")
    pProfile <- rbind(pProfile, pBlim)
  } else {
    flim <- flim10 <- flim50 <- Blim <- NA
  }

  if (!missing(Bpa)) {
    pBpa <- apply(ssbsa > Bpa, 1, mean)
    pBpa <- data.frame(Ftarget = Fscan, value = 1 - pBpa, variable = "Bpa")
    pProfile <- rbind(pProfile, pBpa)
  } else {
    Bpa <- NA
  }

  # GENERATE REF-TABLE
  catF <- c(flim, flim10, flim50, vcumCat, Fscan[maxcatm], FCrash05, FCrash50)
  lanF <- c(NA, NA, NA, vcumLan, Fscan[maxlanm], NA, NA)
  catC <-    approx(Fscan, cats[4, ], xout = catF)$y
  lanC <-    approx(Fscan, lans[4, ], xout = lanF)$y
  catB <-    approx(Fscan, ssbs[4, ], xout = catF)$y
  lanB <-    approx(Fscan, ssbs[4, ], xout = lanF)$y

  Refs <- rbind(catF, lanF, catC, lanC, catB, lanB)
  rownames(Refs) <- c("catF", "lanF", "catch", "landings", "catB", "lanB")
  colnames(Refs) <- c("F05", "F10", "F50", "medianMSY", "meanMSY", "FCrash05", "FCrash50")

  # TODO: id.sim - user specified.

  # 2014-03-12 Ammendments per note from Carmen/John
  # CALCULATIONS:

  # Fmsy: value that maximises median LT catch or median LT landings
  auxi <-    approx(Fscan, cats[4, ], xout = seq(min(Fscan), max(Fscan), length = 200))
  FmsyMedianC <- auxi$x[which.max(auxi$y)]
  MSYMedianC <- max(auxi$y)
  # Value of F that corresponds to 0.95*MSY:
  FmsylowerMedianC <- auxi$x[min((1:length(auxi$y))[auxi$y / MSYMedianC >= 0.95])]
  FmsyupperMedianC <- auxi$x[max((1:length(auxi$y))[auxi$y / MSYMedianC >= 0.95])]

  auxi <-    approx(Fscan, lans[4, ], xout = seq(min(Fscan), max(Fscan), length = 200))
  FmsyMedianL <- auxi$x[which.max(auxi$y)]
  MSYMedianL <- max(auxi$y)

  # Value of F that corresponds to 0.95*MSY:
  FmsylowerMedianL <- auxi$x[min((1:length(auxi$y))[auxi$y / MSYMedianL >= 0.95])]
  FmsyupperMedianL <- auxi$x[max((1:length(auxi$y))[auxi$y / MSYMedianL >= 0.95])]

  F5percRiskBlim <- flim

  refs_interval <- data.frame(
    FmsyMedianC = FmsyMedianC,
    FmsylowerMedianC = FmsylowerMedianC,
    FmsyupperMedianC = FmsyupperMedianC,
    FmsyMedianL = FmsyMedianL,
    FmsylowerMedianL = FmsylowerMedianL,
    FmsyupperMedianL = FmsyupperMedianL,
    F5percRiskBlim = F5percRiskBlim,
    Btrigger = Btrigger
  )

  # END 2014-03-12 Ammendments per note from Carmen/John

  sim <- list(
    ibya = list(
      Mat = Mat, M = M, Fprop = Fprop, Mprop = Mprop,
      west = west, weca = weca, sel = sel
    ),
    rbya = list(
      ferr = ferr, ssb = ssbsa, catch = catsa,
      landings = lansa, rec = recsa,
      srmodels = SR, Ftarget = Fscan, simyears = begin:Nrun
    ),
    rby = fit$rby,
    rbp = rbp,
    Blim = Blim,
    Bpa = Bpa,
    Refs = Refs,
    pProfile = pProfile,
    id.sim = fit$id.sr,
    refs_interval = refs_interval,
    rhologRec = rhologRec
  )

  if (keep.sims) {
    sim$rbya_all <- list(ssb = ssbsall, catch = catsall, landings = lansall, rec = recsall)
  }

  if (verbose) msg("Calculating MSY range values")

  sim <- eqsim_range(sim)

  return(sim)
}