R/oem.R

In flr/mse: Tools for Running Management Strategy Evaluations using FLR

# oem.R - DESC
# mse/R/oem.R

# Copyright European Union, 2018
# Author: Ernesto Jardim (EC JRC) <ernesto.jardim@ec.europa.eu>
#         Iago Mosqueira (EC JRC) <iago.mosqueira@ec.europa.eu>
#
# Distributed under the terms of the European Union Public Licence (EUPL) V.1.1.


# perfect.oem {{{

#' A perfect observation of catch and abundances-at-age.
#'
#' This observation error model function generates a set of perfect observations
#' on catches, biology and abundance. Direct observations are made of the stock,
#' while a single age-structured index of abundance, in numbers, is created with
#' a fixed catchability of 0.01. *deviances* on either *stk$catch.n* or
#' *idx$index*, if given, are applied.
#'
#' This *oem* function generates a full observation time series every time step,
#' and does not append them to existing objects in *observations*.
#'
#' @param om An operating model, class *FLom* or *FLombf*.
#' @param observations A list of past observations, extended to the end of *om*, class *list*.
#' @param deviances A named list of observation deviances, class *list*.
#' @param args The mp dimensions arguments, as generated by `mp`, class *list*.
#' @param tracking Object to track module decisions and outputs, class *FLQuant*.
#'
#' @return A named *list* with elements *stk* (*FLStock*), *idx* (*FLIndices*), *deviances*, *observations* and *tracking*.
#'
#' @examples
#' # On FLom
#' data(sol274)
#' obs <- perfect.oem(stock(om), deviances=NULL, observations=NULL,
#'   args=list(y0=1957, dy=2021), tracking=FLQuant())

perfect.oem <- function(stk, deviances, observations, args, tracking,
  biomass=FALSE, ...) {

  # DIMENSIONS
  y0 <- ac(args$y0)
  dy <- ac(args$dy)

  # GET perfect stock
	stk <- window(stk, start=y0, end=dy, extend=FALSE)
  
  # SIMPLIFY to match observations$stk
  dio <- dim(observations$stk)
  dis <- dim(stk)

  if(!is.null(dio)) {
    if(dio[3] > dis[3])
      stk <- nounit(stk)
    if(dio[4] > dis[4])
      stk <- noseason(stk)
  }

  # SET perfect FLIndex per stock
  if(biomass) {
    abu <- catch(stk) / fbar(stk)
    idx <- FLIndices(A=survey(stk, index.q=0.01, biomass=TRUE))
  } else {
    idx <- FLIndices(A=survey(stk, index.q=0.01))
  }

  # STORE observations
  observations$stk <- stk
  observations$idx <- idx

	list(stk=stk, idx=idx, observations=observations, tracking=tracking)

} # }}}

# shortcut.oem {{{

shortcut.oem <- function(stk, deviances, observations, args, tracking, ...) {

  # DIMENSIONS
  y0 <- ac(args$y0)
  dy <- ac(args$dy)
  dyrs <- ac(seq(args$dy - args$frq + 1, args$dy))

  # GET perfect stock
	stk <- window(stk, start=y0, end=dy, extend=FALSE)
  
  # SIMPLIFY to match observations$stk
  dio <- dim(observations$stk)
  dis <- dim(stk)

  if(!is.null(dio)) {
    if(dio[3] > dis[3])
      stk <- nounit(stk)
    if(dio[4] > dis[4])
      stk <- noseason(stk)
  }

  # ADD deviances
  if(!is.null(deviances$stk)) {

    # APPLY deviances and ASSIGN to stk slots in dyrs
    for(i in names(deviances$stk)) {
      slot(stk, i)[, dyrs] <-
      do.call(i, list(object=stk))[, dyrs] * deviances$stk[[i]][, dyrs] + 0.001
    }

    # COMPUTE aggregated slots
    landings(stk)[, dyrs] <- computeLandings(stk[, dyrs])
    discards(stk)[, dyrs] <- computeDiscards(stk[, dyrs])
    catch(stk)[, dyrs] <- computeCatch(stk[, dyrs])
  }

  # STORE observations
  observations$stk <- stk

	list(stk=stk, idx=FLIndices(), observations=observations, tracking=tracking)

} # }}}

# sampling.oem {{{

#' Samples from an operating model to obtain catch, biology and abundance data
#'
#' This observation error model (OEM) function mimics the most common data
#' collection regime, in which catch-at-age and biology is sampled from the
#' population, and one or more indices of abundance are derived from surveys
#' or CPUE data.
#'
#' The FLStock object passed to *sampling.oem* by the *mp* function is
#' simplified to match the dimensions of that present in the *observations*
#' slot.
#'
#' @param stk An FLStock object as obtained by the call to *stock(om)*.
#' @param deviances A named list of deviances, see Details.
#' @param observations A named list of observations, see Details.
#' @param args Options and arguments passed on by *mp()*.
#' @param tracking The tracking object.
#' @return A named list with elements *stk*, *idx*, *observations* and *tracking*.
#' @author Iago Mosqueira (WUR) & Ernesto Jardim (MSC).
#' @seealso \link{mp}
#' @keywords function
#' @examples
#' data(sol274)
#' # Generate samples from year 2000:2016
#' sampling.oem(stock(om), deviances=deviances(oem),
#'   observations=observations(oem),
#'   args=list(y0=2000, dy=2021, frq=1), tracking=FLQuant())

sampling.oem <- function(stk, deviances, observations, stability=1,
  args, tracking) {

  # - sampling.oem

  # DIMENSIONS
  y0 <- ac(args$y0)
  dy <- ac(args$dy)
  dyrs <- ac(seq(args$dy - args$frq + 1, args$dy))

  # CHECK inputs
  # if(!any(!c("stk", "idx") %in% names(deviances)))
  #  stop("deviances(oem) must have elements 'stk' or 'idx'.")

  if(any(!c("stk", "idx") %in% names(observations)))
    stop("observations(oem) must have elements 'stk' and 'idx'.")

  # TODO: GENERATE length samples from OM: catch and indices, BUT needs selex
  # - invALK
  # - lenSamples

  # CHECK dimensions to simplify, on catch.n for multi-fleet FLStock
  simp <- (dim(catch.n(observations$stk))[c(3,4,5)] == 1) +
    (dim(catch.n(stk))[c(3,4,5)] == 1) < 2

  # SIMPLIFY stk to match dimensions of observations$stk
  if(any(simp))
    stk <- simplify(stk, c("unit", "season", "area")[simp], harvest=FALSE)
  
  # SUBSET year range
	stk <- window(stk, start=y0, end=dy, extend=FALSE)
 
  # --- STK

  if(!is.null(deviances$stk)) {

    # APPLY deviances and ASSIGN to stk slots in dyrs
    for(i in names(deviances$stk)) {
      slot(stk, i)[, dyrs] <-
      do.call(i, list(object=stk))[, dyrs] * deviances$stk[[i]][, dyrs] + 0.001
    }

    # COMPUTE aggregated slots
    landings(stk)[, dyrs] <- computeLandings(stk[, dyrs])
    discards(stk)[, dyrs] <- computeDiscards(stk[, dyrs])
    catch(stk)[, dyrs] <- computeCatch(stk[, dyrs])

  # STORE for shortcut 
  observations$stk[, dyrs] <- stk[, dyrs]

  }

  # --- IDX
  idx <- observations$idx

  # CHOOSE indices to be updated (maxyear >= dy)
  upi <- unlist(lapply(idx, function(x) unname(dims(x)$maxyear) > args$dy))

  if(is.null(deviances$idx) | length(deviances$idx) == 0) {
    deviances$idx <- lapply(observations$idx, function(x) index.q(x) %=% 1)
  }

  # APPLY survey() with deviances$idx on top of index.q

  idx[upi] <- Map(function(x, y, z) {

    dyrs <- intersect(dyrs, dimnames(y)$year)

    # CREATE survey obs
    res <- survey(stk[, dyrs], x[, dyrs], sel=sel.pattern(x)[, dyrs],
      index.q=index.q(x)[, dyrs] * y[, dyrs], stability=z)

    # ENSURE no zeroes coming, maybe from high Fs
    if(sum(index(res)[, dyrs]) == 0)
      index(res)[, dyrs] <- sqrt(.Machine$double.eps)

    # SET 0s to min / 2
    index(res)[index(res) == 0] <- c(min(index(res)[index(res) > 0] / 2))
    
    # ASSIGN index observation
    index(x)[, dyrs] <- index(res)

    return(window(x, end=dy))

  }, x=idx[upi], y=deviances$idx[upi], z=rep(stability, length(idx))[upi])

  for(i in seq(idx[upi])) {
    yrs <- intersect(dyrs, dimnames(idx[upi][[i]])$year)
    observations$idx[upi][[i]][, dyrs]<- idx[upi][[i]][, dyrs]
  }

  list(stk=stk, idx=idx, observations=observations, tracking=tracking)

} # }}}

# default.oem {{{
default.oem <- function(om) {
 
  stk <- stock(om)
  
  # observations match OM
  obs <- list(stk=stk, idx=FLIndices(A=as(stk, 'FLIndex')))

  # deviances are NULL
  devs <- list(idx=FLQuants(A=index.q(obs$idx$A) %=% 1),
    stk=FLQuants(catch.n=catch.n(obs$stk) %=% 1))

  # method is perfect.oem

  return(FLoem(method=perfect.oem, observations=obs))
}
# }}}