# 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))
}
# }}}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.