Nothing
R/Class_definitions.R
In MSEtool: Management Strategy Evaluation Toolkit
Defines functions
show_int
show
importslot
Documented in
show
#' Label class union for performance metric objects
#'
#' @description Used internally. Nothing to see here!
#' @keywords internal
#' @export
setClassUnion(name="label.class", members=c("call", "character", "function"))
#' Prob class union for performance metric objects
#'
#' @description Used internally. Nothing to see here!
#' @keywords internal
#' @export
setClassUnion(name="prob.class", members=c("matrix", "numeric", "data.frame"))
# ---- Data Class ----
#' Class \code{'Data'}
#'
#' An object for storing fishery data for analysis
#'
#'
#' @name Data-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Data', stock)}
#' @slot Name The name of the Data object. Single value. Character string
#' @slot Common_Name Common name of the species. Character string
#' @slot Species Scientific name of the species. Genus and species name. Character string
#' @slot Region Name of the general geographic region of the fishery. Character string
#' @slot LHYear The last historical year of the simulation (before projection). Single value. Positive integer
#' @slot MPrec The previous recommendation of a management procedure. Vector of length nsim. Positive real numbers
#' @slot Units Units of the catch/absolute abundance estimates. Single value. Character string
#' @slot MPeff The current level of effort. Vector of length nsim. Positive real numbers
#' @slot nareas Number of fishing areas. Vector of length nsim. Non-negative integer
#'
#' @slot MaxAge Maximum age. Vector nsim long. Positive integer
#' @slot Mort Natural mortality rate. Vector nsim long. Positive real numbers
#' @slot CV_Mort Coefficient of variation in natural mortality rate. Vector nsim long. Positive real numbers
#' @slot vbLinf Maximum length. Vector nsim long. Positive real numbers
#' @slot CV_vbLinf Coefficient of variation in maximum length. Vector nsim long. Positive real numbers
#' @slot vbK The von Bertalanffy growth coefficient K. Vector nsim long. Positive real numbers
#' @slot CV_vbK Coefficient of variation in the von Bertalanffy K parameter. Vector nsim long. Positive real numbers
#' @slot vbt0 Theoretical age at length zero. Vector nsim long. Non-positive real numbers
#' @slot CV_vbt0 Coefficient of variation in age at length zero. Vector nsim long. Positive real numbers
#' @slot wla Weight-Length parameter alpha. Vector nsim long. Positive real numbers
#' @slot CV_wla Coefficient of variation in weight-length parameter a. Vector nsim long. Positive real numbers
#' @slot wlb Weight-Length parameter beta. Vector nsim long. Positive real numbers
#' @slot CV_wlb Coefficient of variation in weight-length parameter b. Vector nsim long. Positive real numbers
#' @slot steep Steepness of stock-recruitment relationship. Vector nsim long. Value in the range of one-fifth to 1
#' @slot CV_steep Coefficient of variation in steepness. Vector nsim long. Positive real numbers
#' @slot sigmaR Recruitment variability. Vector nsim long. Positive real numbers
#' @slot CV_sigmaR Coefficient of variation in recruitment variability. Vector nsim long. Positive real numbers
#' @slot L50 Length at 50 percent maturity. Vector nsim long. Positive real numbers
#' @slot CV_L50 Coefficient of variation in length at 50 per cent maturity. Vector nsim long. Positive real numbers
#' @slot L95 Length at 95 percent maturity. Vector nsim long. Positive real numbers
#' @slot LenCV Coefficient of variation of length-at-age (assumed constant for all age classes). Vector nsim long. Positive real numbers
#' @slot LFC Length at first capture. Vector nsim long. Positive real numbers
#' @slot CV_LFC Coefficient of variation in length at first capture. Vector nsim long. Positive real numbers
#' @slot LFS Shortest length at full selection. Vector nsim long. Positive real numbers
#' @slot CV_LFS Coefficient of variation in length at full selection. Vector nsim long. Positive real numbers
#' @slot Vmaxlen Vulnerability of individuals at asymptotic length. Vector nsim long. Real number between 0 and 1.
#'
#' @slot Year Years that corresponding to catch and relative abundance data. Vector nyears long. Positive integer
#' @slot Cat Total annual catches. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Cat Coefficient of variation in annual catches. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Cat` for all years.
#' @slot Effort Annual fishing effort. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Effort Coefficient of variation in annual effort. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Effort` for all years.
#' @slot Ind Relative total abundance index. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Ind Coefficient of variation in the relative total abundance index. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Ind` for all years
#'
#' @slot SpInd Relative spawning abundance index. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_SpInd Coefficient of variation in the relative spawning abundance index. Matrix nsim rows and either 1 or nyear columns. Positive real numbers.
#'
#' @slot VInd Relative vulnerable abundance index. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_VInd Coefficient of variation in the relative vulnerable abundance index. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers.
#'
#' @slot AddInd Optional additional indices. Array of dimensions `nsim`, n additional indices, and `nyears` (length `Year`).
#' @slot CV_AddInd Coefficient of variation for additional indices. Array of same dimensions as `AddInd`
#' @slot AddIndV Vulnerability-at-age schedules for the additional indices. Array with dimensions: `nsim`, n additional indices, `MaxAge+1`.
#' @slot AddIunits Units for the additional indices - biomass (1; default) or numbers (0). Numeric vector length n.ind.
#' @slot AddIndType Index calculated from total stock (1, default), spawning stock (2), or vulnerable stock (3). Numeric vector of length n.ind
#'
#' @slot Rec Recent recruitment strength. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Rec Log-normal CV for recent recruitment strength. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Rec` for all years.
#' @slot ML Mean length time series. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot Lc Modal length of catches. Matrix of nsim rows and nyears columns. Positive real numbers
#' @slot Lbar Mean length of catches over Lc. Matrix of nsim rows and nyears columns. Positive real numbers
#' @slot Vuln_CAA Optional vulnerability-at-age schedule for catch-at-age samples. Used to condition OM for closed-loop
#' simulation testing. Replaces the fleet selectivity schedule in the OM used to generate CAA samples. Matrix
#' with dimensions `nsim` x `MaxAge+1`.
#' @slot CAA Catch at Age data (numbers). Array of dimensions nsim x nyears x MaxAge+1. Non-negative integers
#' @slot Vuln_CAL Optional vulnerability-at-length schedule for catch-at-length samples. Used to condition OM for closed-loop
#' simulation testing. Replaces the fleet selectivity schedule in the OM used to generate CAL samples. Matrix
#' with dimensions `nsim` x `length(CAL_mids)`.
#' @slot CAL_bins The values delimiting the length bins for the catch-at-length data. Vector. Non-negative real numbers
#' @slot CAL_mids The values of the mid-points of the length bins. Optional, calculated from `CAL_bins` if not entered. Vector. Non-negative real numbers.
#' @slot CAL Catch-at-length data. An array with dimensions nsim x nyears x length(CAL_mids). Non-negative integers. By default the CAL data will be the retained lengths (i.e, not including discards). If `OM@control$CAL =="removals"` then the CAL data will include all removals (retained + discards).
#'
#' @slot Dep Stock depletion SSB(current)/SSB(unfished). Vector nsim long. Fraction.
#' @slot CV_Dep Coefficient of variation in current stock depletion. Vector nsim long. Positive real numbers
#' @slot Abun An estimate of absolute current vulnerable abundance. Vector nsim long. Positive real numbers
#' @slot CV_Abun Coefficient of variation in estimate of absolute current stock size. Vector nsim long. Positive real numbers
#' @slot SpAbun An estimate of absolute current spawning stock abundance. Vector nsim long. Positive real numbers
#' @slot CV_SpAbun Coefficient of variation in estimate of absolute spawning current stock size. Vector nsim long. Positive real numbers
#'
#' @slot FMSY_M An assumed ratio of FMSY to M. Vector nsim long. Positive real numbers
#' @slot CV_FMSY_M Coefficient of variation in the ratio in FMSY/M. Vector nsim long. Positive real numbers
#' @slot BMSY_B0 The most productive stock size relative to unfished. Vector nsim long. Fraction
#' @slot CV_BMSY_B0 Coefficient of variation in the position of the most productive stock size relative to unfished. Vector nsim long. Positive real numbers
#' @slot Cref Reference or target catch level (eg MSY). Vector of length nsim. Positive real numbers
#' @slot CV_Cref Log-normal CV for reference or target catch level. Vector of length nsim. Positive real numbers
#' @slot Bref Reference or target biomass level (eg BMSY). Vector of length nsim. Positive real numbers
#' @slot CV_Bref Log-normal CV for reference or target biomass level. Vector of length nsim. Positive real numbers
#' @slot Iref Reference or target relative abundance index level (eg BMSY / B0). Vector of length nsim. Positive real numbers
#' @slot CV_Iref Log-normalCV for reference or target relative abundance index level. Vector of length nsim. Positive real numbers
#' @slot t The number of years corresponding to AvC and Dt. Single value. Positive integer
#' @slot AvC Average catch over time t. Vector nsim long. Positive real numbers
#' @slot CV_AvC Coefficient of variation in average catches over time t. Vector nsim long. Positive real numbers
#' @slot Dt Depletion over time t SSB(now)/SSB(now-t+1). Vector nsim long. Fraction
#' @slot CV_Dt Coefficient of variation in depletion over time t. Vector nsim long. Positive real numbers
#' @slot Ref A reference management level (eg a catch limit). Single value. Positive real number
#' @slot Ref_type Type of reference management level (eg 2009 catch limit). Single value. Character string
#' @slot Log A record of events. Single value. Character string
#' @slot params A place to store estimated parameters. An object. R list
#' @slot PosMPs The methods that can be applied to these data. Vector. Character strings
#' @slot TAC The calculated catch limits (function TAC). An array with dimensions PosMPs x replicate TAC samples x nsim. Positive real numbers
#' @slot Sense The results of the sensitivity analysis (function Sense). An array with dimensions PosMPs x sensitivity increments. Positive real numbers
#' @slot MPs The methods that were applied to these data. Vector. Character strings
#' @slot OM A table of operating model conditions. R table object of nsim rows. Real numbers
#' @slot Obs A table of observation model conditions. R table object of nsim rows. Real numbers
#' @slot Misc Other information for MPs. An object. R list
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' newdata<-new('Data')
#'
setClass("Data",
representation(Name = "character", Common_Name='character',
Species='character', Region='character',
LHYear = "numeric", MPrec = "vector",
Units = "character", MPeff = "vector",
nareas = "numeric",
MaxAge = "vector", Mort = "vector", CV_Mort = "vector",
vbLinf = "vector", CV_vbLinf = "vector",
vbK = "vector", CV_vbK = "vector",
vbt0 = "vector", CV_vbt0 = "vector",
wla = "vector", CV_wla = "vector",
wlb = "vector", CV_wlb = "vector",
steep = "vector", CV_steep = "vector",
sigmaR='vector', CV_sigmaR='vector',
L50 = "vector", CV_L50 = "vector",
L95 = "vector",
LenCV="vector",
LFC = "vector", CV_LFC = "vector",
LFS = "vector", CV_LFS = "vector",
Vmaxlen = 'vector',
Year = "vector",
Cat = "matrix", CV_Cat = "matrix",
Effort = 'matrix', CV_Effort = 'matrix',
Ind = "matrix", CV_Ind = "matrix",
SpInd = "matrix", CV_SpInd = "matrix",
VInd = "matrix", CV_VInd = "matrix",
AddInd = "array", CV_AddInd = "array", AddIndV = "array",
AddIunits = 'vector', AddIndType='vector',
Rec = "matrix", CV_Rec = "matrix",
ML = "matrix", Lc = "matrix", Lbar = "matrix",
Vuln_CAA = "matrix", CAA = "array",
Vuln_CAL = 'matrix', CAL_bins = "numeric",
CAL_mids = "numeric", CAL = "array",
Dep = "vector", CV_Dep = "vector",
Abun = "vector", CV_Abun = "vector",
SpAbun= "vector", CV_SpAbun = "vector",
FMSY_M = "vector", CV_FMSY_M = "vector",
BMSY_B0 = "vector", CV_BMSY_B0 = "vector",
Cref = "vector", CV_Cref = "vector",
Bref = "vector", CV_Bref = "vector",
Iref = "vector", CV_Iref = "vector",
t = "vector", AvC = "vector", CV_AvC = "vector",
Dt = "vector", CV_Dt = "vector",
Ref = "numeric", Ref_type = "character",
Log = "list", params = "list", PosMPs = "vector",
TAC = "array", Sense = "array",
MPs = "vector", OM = "data.frame", Obs = "data.frame",
Misc = "list"))
setMethod("initialize", "Data", function(.Object, stock="nada", ...) {
if (file.exists(stock)) {
.Object <- XL2Data(stock, ...)
} else {
slots <- slotNames('Data')
for (x in seq_along(slots)) {
sl <- slots[x]
cl <- class(slot(.Object, sl))
if ("logical" %in% cl) slot(.Object, sl) <- as.numeric(NA)
if ("character" %in% cl) slot(.Object, sl) <- ''
if ("matrix" %in% cl) {
if (length(dim(slot(.Object, sl))) > 2) {
slot(.Object, sl) <- array(NA, dim=c(1,1,1))
} else {
slot(.Object, sl) <- matrix(NA)
}
}
if ("array" %in% cl) {
if (length(dim(slot(.Object, sl))) > 2) {
slot(.Object, sl) <- array(NA, dim=c(1,1,1))
} else {
slot(.Object, sl) <- matrix(NA)
}
}
if ("vector" %in% cl) slot(.Object, sl) <- NA
if ("numeric" %in% cl) slot(.Object, sl) <- as.numeric(NA)
if ("list" %in% cl) slot(.Object, sl) <- list()
}
}
# Default values
if (all(is.na(.Object@CV_Cat))) .Object@CV_Cat <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_Ind))) .Object@CV_Ind <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_SpInd))) .Object@CV_SpInd <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_VInd))) .Object@CV_VInd <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_Effort))) .Object@CV_Effort <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_Rec))) .Object@CV_Rec <- matrix(0.2, nrow=1, ncol=1)
if (NAor0(.Object@LenCV)) .Object@LenCV <- 0.1
if (NAor0(.Object@CV_Dt)) .Object@CV_Dt <- 0.25
if (NAor0(.Object@CV_AvC)) .Object@CV_AvC <- 0.2
if (NAor0(.Object@CV_Mort)) .Object@CV_Mort <- 0.2
if (NAor0(.Object@CV_FMSY_M)) .Object@CV_FMSY_M <- 0.2
if (NAor0(.Object@CV_BMSY_B0)) .Object@CV_BMSY_B0 <- 0.045
if (NAor0(.Object@CV_Cref)) .Object@CV_Cref <- 0.2
if (NAor0(.Object@CV_Bref)) .Object@CV_Bref <- 0.2
if (NAor0(.Object@CV_Iref)) .Object@CV_Iref <- 0.2
if (NAor0(.Object@CV_Dep)) .Object@CV_Dep <- 0.25
if (NAor0(.Object@CV_Abun)) .Object@CV_Abun <- 0.25
if (NAor0(.Object@CV_vbK)) .Object@CV_vbK <- 0.1
if (NAor0(.Object@CV_vbLinf)) .Object@CV_vbLinf <- 0.1
if (NAor0(.Object@CV_vbt0)) .Object@CV_vbt0 <- 0.1
if (NAor0(.Object@CV_L50)) .Object@CV_L50 <- 0.1
if (NAor0(.Object@CV_LFC)) .Object@CV_LFC <- 0.2
if (NAor0(.Object@CV_LFS)) .Object@CV_LFS <- 0.2
if (NAor0(.Object@CV_wla)) .Object@CV_wla <- 0.1
if (NAor0(.Object@CV_wlb)) .Object@CV_wlb <- 0.1
if (NAor0(.Object@CV_steep)) .Object@CV_steep <- 0.2
if (NAor0(.Object@nareas)) .Object@nareas <- 2
if (all(is.na(.Object@CAL))) .Object@CAA <- array(NA, c(1, 1, 1))
if (all(is.na(.Object@CAL))) .Object@CAL <- array(NA, c(1, 1, 1))
if (length(.Object@CAL_bins) == 0) .Object@CAL_bins <- 1
if (all(is.na(.Object@CAL))) .Object@AddInd <- array(NA, c(1, 1, 1))
if (all(is.na(.Object@CAL))) .Object@CV_AddInd <- array(NA, c(1, 1, 1))
if (all(is.na(.Object@CAL))) .Object@AddIndV <- array(NA, c(1, 1, 1))
if (length(.Object@TAC) == 0) .Object@TAC <- array(1, c(1, 1))
# if (length(.Object@TACbias) == 0) .Object@TACbias <- array(1, c(1, 1))
if (length(.Object@Sense) == 0) .Object@Sense <- array(1, c(1, 1))
if (length(.Object@ML) == 0) .Object@ML <- array(NA, c(1, 1))
if (length(.Object@Lbar) == 0) .Object@Lbar <- array(NA, c(1, 1))
if (length(.Object@Lc) == 0) .Object@Lc <- array(NA, c(1, 1))
return(.Object)
})
importslot <- function(name, length=2, Data, Names, numeric=TRUE, essential=TRUE) {
x <- Data[match(name, Names), 1:length]
if (numeric) x <- as.numeric(x)
if (essential) {
if (any(is.na(x))) warning('NAs in ', name, '. Should be length ', length)
if (length(x)<length) warning(name, ' is missing')
}
x
}
# ---- Stock Class ----
#' Class \code{'Stock'}
#'
#' An operating model component that specifies the parameters of the population
#' dynamics model
#'
#' @docType class
#'
#' @template Stock_template
#'
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Stock')}
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' showClass('Stock')
#'
setClass("Stock", representation(Name = "character",
Common_Name='character',
Species="character",
maxage = "numeric",
R0 = "numeric",
M = "numeric",
Msd = "numeric",
h = "numeric",
SRrel = "numeric",
Perr = "numeric",
AC = "numeric",
Linf = "numeric",
Linfsd = "numeric",
K = "numeric",
Ksd = "numeric",
t0 = "numeric",
LenCV="numeric",
L50 = "numeric",
L50_95 = "numeric",
D = "numeric",
a = "numeric",
b = "numeric",
Size_area_1 = "numeric",
Frac_area_1 = "numeric",
Prob_staying = "numeric",
Fdisc="numeric",
Source = "character"))
setMethod("initialize", "Stock", function(.Object, file = NA, dec=c(".", ",")) {
if (!is.na(file)) {
if (!file.exists(file)) {
stop('File not found')
} else {
# load file and populate object
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec == ".") dat <- read.csv(file, header = F,
colClasses = "character",
col.names = paste0("V", 1:Ncol))
if (dec == ",") dat <- read.csv2(file, header = F,
colClasses = "character",
col.names = paste0("V", 1:Ncol))
Names <- dat[,1]
Ncol <- ncol(dat)
Data <- dat[,2:Ncol]
.Object@Name <- importslot('Name', 1, Data, Names, FALSE)
commonname <- Data[match("Common_Name", Names), 1]
commonname1 <- Data[match("Common Name", Names), 1]
if (length(commonname)>0 && !is.na(commonname)) {
.Object@Common_Name <-commonname
} else if (length(commonname1)>0 && !is.na(commonname)) {
.Object@Common_Name <-commonname1
} else {
.Object@Common_Name <- "Not specified"
}
.Object@Species <- importslot('Species', 1, Data, Names, FALSE, FALSE)
.Object@maxage <- importslot('maxage', 1, Data, Names)
.Object@R0 <- importslot('R0', 1, Data, Names)
.Object@M <- importslot('M', 2, Data, Names)
.Object@Msd <- importslot('Msd', 2, Data, Names)
.Object@Fdisc <- importslot('Fdisc', 2, Data, Names)
.Object@h <- importslot('h', 2, Data, Names)
.Object@SRrel <- importslot('SRrel', 1, Data, Names)
.Object@Linf <- importslot('Linf', 2, Data, Names)
.Object@K <- importslot("K", 2, Data, Names)
.Object@t0 <- importslot("t0", 2, Data, Names)
.Object@LenCV <- importslot("LenCV", 2, Data, Names)
.Object@Ksd <- importslot("Ksd", 2, Data, Names)
.Object@Linfsd <- importslot("Linfsd", 2, Data, Names)
.Object@a <- importslot("a", 1, Data, Names)
.Object@b <- importslot("b", 1, Data, Names)
.Object@D <- importslot("D", 2, Data, Names)
.Object@Perr <- importslot("Perr", 2, Data, Names)
.Object@AC <- importslot("AC", 2, Data, Names)
.Object@Size_area_1 <- importslot("Size_area_1", 2, Data, Names)
.Object@Frac_area_1 <- importslot("Frac_area_1", 2, Data, Names)
.Object@Prob_staying <- importslot("Prob_staying", 2, Data, Names)
.Object@L50 <- importslot("L50", 2, Data, Names)
.Object@L50_95 <- importslot("L50_95", 2, Data, Names)
.Object@Source <- importslot("Source", 1, Data, Names, FALSE, FALSE)
}
}
.Object
})
#' Class union for isRel slot
#'
#' @description Used internally.
#' @keywords internal
#' @export
setClassUnion(name="char.log", members=c("character", "logical"))
# ---- Fleet Class -----
#' Class \code{'Fleet'}
#'
#' The component of the operating model that controls fishing dynamics
#'
#' @name Fleet-class
#' @docType class
#' @template Fleet_template
#'
#' @section Creating Object:
#' Objects can be created by calls of the form \code{new('Fleet')}
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#' showClass('Fleet')
#'
setClass("Fleet", slots = c(Name = "character",
nyears = "numeric",
CurrentYr="numeric",
EffYears = "numeric",
EffLower = "numeric",
EffUpper = "numeric",
Esd = "numeric",
qinc = "numeric",
qcv = "numeric",
L5 = "numeric",
LFS = "numeric",
Vmaxlen = "numeric",
isRel = "char.log",
LR5 = "numeric",
LFR = "numeric",
Rmaxlen = "numeric",
DR = "numeric",
Spat_targ = "numeric",
MPA='char.log',
Misc='list')
)
# initialize Fleet
setMethod("initialize", "Fleet", function(.Object, file = NA, dec=c(".", ",")) {
if (!is.na(file)) {
if (!file.exists(file)) {
stop('File not found')
} else {
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec == ".") dat <- read.csv(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol))
if (dec == ",") dat <- read.csv2(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol))
Names <- dat[, 1]
Data <- dat[, 2:ncol(dat)]
.Object@Name <- importslot('Name', 1, Data, Names, FALSE)
.Object@nyears <- importslot('nyears', 1, Data, Names)
.Object@CurrentYr <- importslot('CurrentYr', 1, Data, Names, TRUE, FALSE)
if(is.na(.Object@CurrentYr)).Object@CurrentYr<-.Object@nyears
.Object@Spat_targ <- importslot('Spat_targ', 2, Data, Names)
if (!is.na(match("Esd", Names))) {
.Object@Esd <- importslot('Esd', 2, Data, Names)
} else {
.Object@Esd <- importslot('Fsd', 2, Data, Names)
}
nEffYears <- ncol(Data[match("EffYears", Names), ])
oldw <- getOption("warn")
options(warn = -1)
chk <- as.numeric(Data[match("EffYears", Names), 1:nEffYears])
options(warn = oldw)
ind <- which(!is.na(chk))
nEffYears <- length(ind)
.Object@EffYears <- importslot('EffYears', nEffYears, Data, Names)
.Object@EffLower <- importslot('EffLower', nEffYears, Data, Names)
.Object@EffUpper <- importslot('EffUpper', nEffYears, Data, Names)
.Object@qinc <- importslot("qinc", 2, Data, Names)
.Object@qcv <- importslot("qcv", 2, Data, Names)
.Object@L5 <- importslot("L5", 2, Data, Names)
.Object@LFS <- importslot("LFS", 2, Data, Names)
.Object@Vmaxlen <- importslot("Vmaxlen", 2, Data, Names)
.Object@LR5 <- importslot("LR5", 2, Data, Names, essential = FALSE)
if (all(is.na(.Object@LR5))) .Object@LR5 <- c(0,0)
.Object@LFR <- importslot("LFR", 2, Data, Names, essential = FALSE)
if (all(is.na(.Object@LFR))) .Object@LFR <- c(0.001,0.001)
.Object@Rmaxlen <- importslot("Rmaxlen", 2, Data, Names, essential=FALSE)
if (all(is.na(.Object@Rmaxlen))) .Object@Rmaxlen <- c(1,1)
.Object@DR <- importslot("DR", 2, Data, Names)
.Object@isRel <- importslot("isRel", 1, Data, Names, FALSE, FALSE)
if (NAor0(.Object@isRel)) .Object@isRel <- "TRUE"
.Object@isRel <- as.character(.Object@isRel)
.Object@MPA <- importslot("MPA", 1, Data, Names, FALSE, FALSE)
if (NAor0(.Object@MPA)) .Object@MPA <- FALSE
.Object@MPA <- as.logical(.Object@MPA)
}
}
if (NAor0(.Object@MPA)) .Object@MPA <- FALSE
.Object
})
#' ~~ Methods for Function \code{initialize} ~~
#'
#' @name initialize-methods
#' @aliases initialize-methods initialize,Data-method
#' initialize,Fleet-method initialize,MSE-method initialize,Obs-method
#' initialize,OM-method initialize,Stock-method initialize
#' initialize,Fease-method initialize,DLM_general-method
#' @docType methods
#' @section Methods: \describe{
#'
#' \item{list('signature(.Object = \'DLM\')')}{ %% ~~describe this method
#' here~~ }
#'
#' \item{list('signature(.Object = \'Fleet\')')}{ %% ~~describe this method
#' here~~ }
#'
#' \item{list('signature(.Object = \'MSE\')')}{ %% ~~describe this method
#' here~~ }
#'
#' \item{list('signature(.Object = \'Obs\')')}{ %% ~~describe this
#' method here~~ }
#'
#' \item{list('signature(.Object = \'OM\')')}{ %% ~~describe this method here~~
#' }
#'
#' \item{list('signature(.Object = \'Stock\')')}{ %% ~~describe this method
#' here~~ }
#'
#'
#' \item{list('signature(.Object = \'Fease\')')}{ %% ~~describe this method
#' here~~ } \item{list('signature(.Object = \'DLM_general\')')}{ %% ~~describe
#' this method here~~ }
#'
#' }
#' @keywords methods ~~ other possible keyword(s) ~~
NULL
# ---- Obs Class ----
#' Class \code{'Obs'}
#'
#' An operating model component that controls the observation model
#'
#'
#' @docType class
#' @note Its questionable whether the hyperstability/hyperdepletion should be
#' categorised as an observation model characteristic as it is most often
#' driven by fleet dynamics (and therefore should be in the fleet object). Oh
#' well its here and you might want to make it hyperstable beta < 1 or
#' hyperdeplete beta > 1, only.
#'
#' @slot Name The name of the observation model object. Single value. Character string.
#' @template Obs_template
#'
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Obs')}
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' showClass('Obs')
#'
setClass("Obs", representation(Name = "character",
Cobs = "numeric",
Cbiascv = "numeric",
CAA_nsamp = "numeric",
CAA_ESS = "numeric",
CAL_nsamp = "numeric",
CAL_ESS = "numeric",
Iobs = "numeric",
Btobs = "numeric",
Btbiascv = "numeric",
beta = "numeric",
LenMbiascv = "numeric",
Mbiascv = "numeric",
Kbiascv = "numeric",
t0biascv = "numeric",
Linfbiascv = "numeric",
LFCbiascv = "numeric",
LFSbiascv = "numeric",
FMSY_Mbiascv = "numeric",
BMSY_B0biascv = "numeric",
Irefbiascv = "numeric",
Brefbiascv = "numeric",
Crefbiascv = "numeric",
Dbiascv = "numeric",
Dobs = "numeric",
hbiascv = "numeric",
Recbiascv = "numeric",
sigmaRbiascv = 'numeric',
Eobs='numeric',
Ebiascv='numeric'))
# initialize Obs
setMethod("initialize", "Obs", function(.Object, file = NA, dec=c(".", ",")) {
if (!is.na(file)) {
if (file.exists(file)) {
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec ==".") dat <- read.csv(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
if (dec ==",") dat <- read.csv2(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
dname <- dat[, 1]
dat <- dat[, 2:ncol(dat)]
.Object@Name <- dat[match("Name", dname), 1]
.Object@Cobs <- as.numeric(dat[match("Cobs", dname), 1:2])
.Object@Cbiascv <- as.numeric(dat[match("Cbiascv", dname), 1])
.Object@CAA_nsamp <- as.numeric(dat[match("CAA_nsamp", dname), 1:2])
.Object@CAA_ESS <- as.numeric(dat[match("CAA_ESS", dname), 1:2])
.Object@CAL_nsamp <- as.numeric(dat[match("CAL_nsamp", dname), 1:2])
.Object@CAL_ESS <- as.numeric(dat[match("CAL_ESS", dname), 1:2])
.Object@Iobs <- as.numeric(dat[match("Iobs", dname), 1:2])
.Object@Btobs <- as.numeric(dat[match("Btobs", dname), 1:2])
.Object@Btbiascv <- as.numeric(dat[match("Btbiascv", dname), 1])
.Object@beta <- as.numeric(dat[match("beta", dname), 1:2])
.Object@LenMbiascv <- as.numeric(dat[match("LenMbiascv", dname), 1])
.Object@Mbiascv <- as.numeric(dat[match("Mbiascv", dname), 1])
.Object@Kbiascv <- as.numeric(dat[match("Kbiascv", dname), 1])
.Object@t0biascv <- as.numeric(dat[match("t0biascv", dname), 1])
.Object@Linfbiascv <- as.numeric(dat[match("Linfbiascv", dname), 1])
.Object@LFCbiascv <- as.numeric(dat[match("LFCbiascv", dname), 1])
.Object@LFSbiascv <- as.numeric(dat[match("LFSbiascv", dname), 1])
.Object@FMSY_Mbiascv <- as.numeric(dat[match("FMSY_Mbiascv", dname), 1])
.Object@BMSY_B0biascv <- as.numeric(dat[match("BMSY_B0biascv", dname), 1])
.Object@Irefbiascv <- as.numeric(dat[match("Irefbiascv", dname), 1])
.Object@Crefbiascv <- as.numeric(dat[match("Crefbiascv", dname), 1])
.Object@Brefbiascv <- as.numeric(dat[match("Brefbiascv", dname), 1])
.Object@Dbiascv <- as.numeric(dat[match("Dbiascv", dname), 1])
.Object@Dobs <- as.numeric(dat[match("Dobs", dname), 1:2])
.Object@hbiascv <- as.numeric(dat[match("hbiascv", dname), 1])
.Object@Recbiascv <- as.numeric(dat[match("Recbiascv", dname), 1:2])
.Object@sigmaRbiascv <- as.numeric(dat[match("Recbiascv", dname), 1])
.Object@Eobs <- as.numeric(dat[match("Eobs", dname), 1:2])
.Object@Ebiascv <- as.numeric(dat[match("Ebiascv", dname), 1])
} else {
message("File doesn't exist")
}
}
.Object
})
# ---- Imp Class ----
#' Class \code{'Imp'}
#'
#' An operating model component that specifies the degree of adherence to management recommendations (Implementation error)
#'
#' @name Imp-class
#' @docType class
#' @slot Name The name of the Implementation error object. Single value. Character string.
#'
#' @template Imp_template
#'
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Imp')}#'
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' showClass('Imp')
#'
setClass("Imp", representation(Name = "character",
TACFrac = "numeric",
TACSD = "numeric",
TAEFrac = "numeric",
TAESD = "numeric",
SizeLimFrac="numeric",
SizeLimSD = "numeric"))
# initialize Imp
setMethod("initialize", "Imp", function(.Object, file = NA, dec=c(".", ",")) {
.Object@Name <- "Perfect implementation"
.Object@TACSD <- c(0,0)
.Object@TACFrac <- c(1,1)
.Object@TAESD <- c(0,0)
.Object@TAEFrac <-c(1,1)
.Object@SizeLimSD <- c(0,0)
.Object@SizeLimFrac<-c(1,1)
# .Object@Source <-"DLMtool generated"
if (!is.na(file)) {
if (file.exists(file)) {
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec ==".") dat <- read.csv(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
if (dec ==",") dat <- read.csv2(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
dname <- dat[, 1]
dat <- dat[, 2:ncol(dat)]
.Object@Name <- dat[match("Name", dname), 1]
.Object@TACSD <- as.numeric(dat[match("TACSD", dname), 1:2])
.Object@TACFrac <- as.numeric(dat[match("TACFrac", dname), 1:2])
.Object@TAESD <- as.numeric(dat[match("TAESD", dname), 1:2])
.Object@TAEFrac <- as.numeric(dat[match("TAEFrac", dname), 1:2])
.Object@SizeLimSD <- as.numeric(dat[match("SizeLimSD", dname), 1:2])
.Object@SizeLimFrac <- as.numeric(dat[match("SizeLimFrac", dname), 1:2])
# .Object@Source <- dat[match("Source", dname), 1]
} else {
message("File doesn't exist")
}
}
.Object
})
# ---- OM Class ----
#' Class \code{'OM'}
#'
#' An object containing all the parameters needed to control the MSE which can
#' be build from component Stock, Fleet, Obs, and Imp objects.
#'
#' Almost all of these inputs are a vector of length 2 which describes the upper and lower
#' bounds of a uniform distribution from which to sample the parameter.
#'
#' @name OM-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('OM', Stock, Fleet, Obs, Imp)}.
#' @slot Name Name of the operating model
#' @slot Agency Name of the agency responsible for the management of the fishery. Character string
#' @slot Region Name of the general geographic region of the fishery. Character string
#' @slot Sponsor Name of the organization who sponsored the OM. Character string
#' @slot Latitude Latitude (decimal degrees). Negative values represent the South of the Equator. Numeric. Single value
#' @slot Longitude Longitude (decimal degrees). Negative values represent the West of the Prime Meridian. Numeric. Single value
#' @slot nsim The number of simulations
#' @slot proyears The number of projected years
#' @slot interval The assessment interval - how often would you like to update the management system?
#' @slot pstar The percentile of the sample of the management recommendation for each method
#' @slot maxF Maximum instantaneous fishing mortality rate that may be simulated for any given age class
#' @slot reps Number of samples of the management recommendation for each method. Note that when this is set to 1, the mean value of
#' the data inputs is used.
#' @slot cpars A list of custom parameters. Time series are a matrix nsim rows by nyears columns. Single parameters are a vector nsim long
#' @slot seed A random seed to ensure users can reproduce results exactly
#' @slot Source A reference to a website or article from which parameters were taken to define the operating model
#'
# #' @template Stock_template
# #' @template Fleet_template
# #' @template Obs_template
# #' @template Imp_template
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#'
setClass("OM", representation(Name = "character", Agency="character",
Region="character", Sponsor="character",
Latitude="numeric", Longitude="numeric",
nsim="numeric", proyears="numeric",
interval='numeric', pstar='numeric', maxF='numeric', reps='numeric',
cpars="list",seed="numeric", Source="character"),
contains=c("Stock", "Fleet", "Obs", "Imp"))
# initialize OM
setMethod("initialize", "OM", function(.Object, Stock=NULL, Fleet=MSEtool::Generic_Fleet,
Obs=MSEtool::Generic_Obs, Imp=MSEtool::Perfect_Imp,
interval=4, pstar=0.5, maxF=0.8, reps=1, nsim=48, proyears=50, docheck=TRUE) {
if (is.null(Stock)) {
# message("No Stock object found. Returning a blank OM object")
# Check and add defaults
if (docheck)
.Object <- CheckOM(.Object, msg=FALSE, stop_if_missing=FALSE)
# Default MSE parameters
if (.hasSlot(.Object, "nsim")) .Object@nsim <- nsim
if (.hasSlot(.Object, "proyears")) .Object@proyears <- proyears
# interval, pstar, maxF, reps
if (.hasSlot(.Object, "interval")) .Object@interval <- interval
if (.hasSlot(.Object, "pstar")) .Object@pstar <- pstar
if (.hasSlot(.Object, "maxF")) .Object@maxF <- maxF
if (.hasSlot(.Object, "reps")) .Object@reps <- reps
return(.Object)
}
if (!methods::is(Stock, "Stock"))
print(paste("Could not build operating model:", deparse(substitute(Stock)), "not of class Stock"))
if (!methods::is(Fleet, "Fleet"))
print(paste("Could not build operating model:", deparse(substitute(Fleet)), "not of class Fleet"))
if (!methods::is(Obs, "Obs"))
print(paste("Could not build operating model:", deparse(substitute(Obs)), "not of class Obs"))
if (!methods::is(Imp, "Imp"))
print(paste("Could not build operating model:", deparse(substitute(Imp)), "not of class Imp"))
.Object@Name <- paste("Stock:", Stock@Name, " Fleet:", Fleet@Name, " Obs model:",
Obs@Name, " Imp model:", Imp@Name, sep = "")
# Now copy the values for stock, fleet and observation slots to same
# slots in the Sim object
Sslots <- slotNames(Stock)
for (i in 2:length(Sslots)) {
tt <- .hasSlot(Stock, Sslots[i]) # For back-compatibility
if (tt) slot(.Object, Sslots[i]) <- slot(Stock, Sslots[i])
}
Fslots <- slotNames(Fleet)
for (i in 2:length(Fslots)) {
tt <- .hasSlot(Fleet, Fslots[i])
if (tt) slot(.Object, Fslots[i]) <- slot(Fleet, Fslots[i])
}
Oslots <- slotNames(Obs)
for (i in 2:length(Oslots)) {
tt <- .hasSlot(Obs, Oslots[i])
if (tt) slot(.Object, Oslots[i]) <- slot(Obs, Oslots[i])
}
Islots <- slotNames(Imp)
for (i in 2:length(Islots)) {
tt <- .hasSlot(Imp, Islots[i])
if (tt) slot(.Object, Islots[i]) <- slot(Imp, Islots[i])
}
#
# source <- paste("Stock:", Stock@Source, "Fleet:", Fleet@Source, "Obs:", Obs@Source, "Imp:",Imp@Source)
slot(.Object, "Source") <- Stock@Source
# Default MSE parameters
if (.hasSlot(.Object, "nsim")) .Object@nsim <- nsim
if (.hasSlot(.Object, "proyears")) .Object@proyears <- proyears
# interval, pstar, maxF, reps
if (.hasSlot(.Object, "interval")) .Object@interval <- interval
if (.hasSlot(.Object, "pstar")) .Object@pstar <- pstar
if (.hasSlot(.Object, "maxF")) .Object@maxF <- maxF
if (.hasSlot(.Object, "reps")) .Object@reps <- reps
# Check and add defaults
.Object <- CheckOM(.Object, msg=FALSE, stop_if_missing=FALSE)
# if(length(.Object@LenCV) < 2) .Object@LenCV <- c(0.08,0.15)
# if(length(.Object@CurrentYr)==0).Object@CurrentYr=.Object@nyears
#
# # if(length(.Object@FecB) < 2) .Object@FecB <- c(3,3)
# # if(all(is.na(.Object@FecB))) .Object@FecB <- c(3,3)
# if(all(is.na(.Object@LenCV))) .Object@LenCV <- c(0.08,0.15)
# if(all(is.na(.Object@CurrentYr))) .Object@CurrentYr=.Object@nyears
#
# if(length(.Object@LR5) < 2) .Object@LR5 <- c(0,0)
# if(length(.Object@LFR) < 2) .Object@LFR <- c(0,0)
# if(length(.Object@Rmaxlen) < 2) .Object@Rmaxlen <- c(1,1)
# if(length(.Object@Fdisc) < 2) .Object@Fdisc <- c(0,0)
#
# if(all(is.na(.Object@LR5))) .Object@LR5 <- c(0,0)
# if(all(is.na(.Object@LFR))) .Object@LFR <- c(0,0)
# if(all(is.na(.Object@Rmaxlen))) .Object@Rmaxlen <- c(1,1)
# if(all(is.na(.Object@Fdisc))) .Object@Fdisc <- c(0,0)
#
# if (.hasSlot(.Object, "Size_area_1")) {
# if (length(.Object@Size_area_1)==0) .Object@Size_area_1 <- .Object@Frac_area_1
# if (all(is.na(.Object@Size_area_1))) .Object@Size_area_1 <- .Object@Frac_area_1
# }
#
# .Object@seed=1
.Object
})
# -- Hist Object Class ----
#' Class \code{'Hist'}
#'
#' An object for storing information generated by the end of the historical simulations
#'
#' @name Hist-class
#' @docType class
#'
#' @slot Data The Data object at the end of the historical period
#' @slot OMPars A numeric data.frame with nsim rows with sampled Stock, Fleet,
#' Obs, and Imp parameters.
#' @slot AtAge A named list with arrays of dimensions: `c(nsim, maxage+1, nyears+proyears)` or
#' `c(nsim, maxage+1, nyears, nareas)`
#' \itemize{
#' \item Length: Length-at-age for each simulation, age, and year
#' \item Weight: Weight-at-age for each simulation, age, and year
#' \item Select: Selectivity-at-age for each simulation, age, and year
#' \item Retention: Retention-at-age for each simulation, age, and year
#' \item Maturity: Maturity-at-age for each simulation, age, and year
#' \item N.Mortality: Natural mortality-at-age for each simulation, age, and year
#' \item Z.Mortality: Total mortality-at-age for each simulation, age, year and area
#' \item F.Mortality: Fishing mortality-at-age for each simulation, age, year and area
#' \item Fret.Mortality: Fishing mortality-at-age for retained fish for each
#' simulation, age, year and area
#' \item Number: Total numbers by simulation, age, year and area
#' \item Biomass: Total biomass by simulation, age, year and area
#' \item VBiomass: Vulnerable biomass by simulation, age, year and area
#' \item SBiomass: Spawning biomass by simulation, age, year and area
#' \item Removals: Removals (biomass) by simulation, age, year and area
#' \item Landings: Landings (biomass) by simulation, age, year and area
#' \item Discards: Discards (biomass) by simulation, age, year and area
#' }
#' @slot TSdata A named list with population and fleet dynamics:
#' \itemize{
#' \item Number: Total numbers; array dimensions `c(nsim, nyears, nareas)`
#' \item Biomass: Total biomass; array dimensions `c(nsim, nyears, nareas)`
#' \item VBiomass: Vulnerable biomass; array dimensions `c(nsim, nyears, nareas)`
#' \item SBiomass: Spawning Biomass; array dimensions `c(nsim, nyears, nareas)`
#' \item Removals: Removals (biomass); array dimensions `c(nsim, nyears, nareas)`
#' \item Landings: Landings (biomass); array dimensions `c(nsim, nyears, nareas)`
#' \item Discards: Discards (biomass); array dimensions `c(nsim, nyears, nareas)`
#' \item Find: Historical fishing mortality (scale-free); matrix dimensions `c(nsim, nyears)`
#' \item RecDev: Recruitment deviations (historical and projection); matrix dimensions `c(nsim, nyears+proyears+maxage)`
#' \item SPR: Named list with Equilibrium and Dynamic SPR (both matrices iwth dimensions `c(nsim, nyears)`)
#' \item Unfished_Equilibrium: A named list with unfished equilibrium numbers and biomass-at-age
#' }
#'
#' @slot Ref A named list with biological reference points:
#' \itemize{
#' \item ByYear: A named list with asymptotic reference points (i.e., calculated annually without recruitment deviations) all matrices with dimensions `nsim` by `nyears+proyears`:
#' \itemize{
#' \item N0: Asymptotic unfished total number
#' \item SN0: Asymptotic unfished spawning number
#' \item B0: Asymptotic unfished total biomass
#' \item SSB0: Asymptotic unfished spawning biomass
#' \item VB0: Asymptotic unfished vulnerable biomass
#' \item MSY: Asymptotic MSY
#' \item FMSY: Fishing mortality corresponding with asymptotic MSY
#' \item SSBMSY: Spawning stock biomass corresponding with asymptotic MSY
#' \item BMSY: total biomass corresponding with asymptotic MSY
#' \item VBMSY: Vulnerable biomass corresponding with asymptotic MSY
#' \item F01: Fishing mortality where the change in yield per recruit is 10% of that at F = 0
#' \item Fmax: Fishing mortality that maximizes yield per recruit
#' \item F_SPR: Fishing mortality corresponding to spawning potential ratio of 20 - 60% in increments of 5%; array dimensions \code{c(nsim, 9, nyears+proyears)}
#' \item Fcrash: Fishing mortality corresponding to the recruits-per-spawner at the origin of the stock-recruit relationship
#' \item Fmed: Fishing mortality corresponding to the median recruits-per-spawner in the historical period
#' \item SPRcrash: SPR corresponding to the recruits-per-spawner at the origin of the stock-recruit relationship
#' }
#' \item Dynamic_Unfished: A named list with dynamic unfished reference points for each simulation and year:
#' \itemize{
#' \item N0: Unfished total numbers
#' \item B0: Unfished total biomass
#' \item SN0: Unfished spawning numbers
#' \item SSB0: Unfished spawning biomass
#' \item VB0: Unfished vulnerable biomass
#' \item Rec: Unfished recruitment
#' }
#' \item ReferencePoints: A data.frame with `nsim` rows with with biological reference points
#' calculated as an average over age-of-maturity `ageM` years around the
#' current year (i.e. `nyears`):
#' \itemize{
#' \item N0: Average unfished numbers
#' \item B0: Average unfished biomass
#' \item SSB0: Average unfished spawning biomass (used to calculate depletion)
#' \item SSN0: Average unfished spawning numbers
#' \item VB0: Average unfished vulnerable biomass (used to calculate depletion if `cpar$control$D='VB'`)
#' \item MSY: Average maximum sustainable yield (equilibrium)
#' \item FMSY: Average fishing mortality corresponding with MSY
#' \item SSBMSY: Average spawning stock biomass corresponding with MSY
#' \item BMSY: Average total biomass corresponding with MSY
#' \item VBMSY: Average vulnerable biomass corresponding with MSY
#' \item UMSY: Average exploitation rate corresponding with MSY
#' \item FMSY_M: Average FMSY/M ratio
#' \item SSBMSY_SSB0: Average ratio of SSBMSY to SSB0
#' \item BMSY_B0: Average ratio of BMSY to B0
#' \item VBMSY_VB0: Average ratio of VBMSY to VB0
#' \item RefY: Maximum yield obtained in forward projections with a fixed F
#' }
#' }
#'
#' @slot SampPars A named list with all sampled Stock, Fleet, Obs, and Imp parameters
#' @slot OM The `OM` object (without cpars)
#' @slot Misc A list for additional information
#' @author A. Hordyk
#' @keywords classes
setClass("Hist", representation(
Data = 'Data',
OMPars = 'data.frame',
AtAge = 'list',
TSdata = 'list',
Ref = "list",
SampPars='list',
Misc = 'list',
OM = 'OM'
))
# ---- MSE Class ----
#' Class \code{'MSE'}
#'
#' A Management Strategy Evaluation object that contains information about
#' simulation conditions and performance of data-limited methods
#'
#'
#' @name MSE-class
#' @docType class
#' @template MSE_template
#'
#'
#' @author T. Carruthers and A. Hordyk
#' @keywords classes
setClass("MSE", representation(Name = "character",
nyears = "numeric",
proyears = "numeric",
nMPs = "numeric",
MPs = "character",
nsim = "numeric",
OM = "data.frame",
Obs = "data.frame",
SB_SBMSY='array',
F_FMSY='array',
N='array',
B = "array",
SSB="array",
VB="array",
FM = "array",
SPR='list',
Catch = "array",
Removals='array',
Effort='array',
TAC = "array",
TAE = 'array',
BioEco = 'list',
RefPoint = 'list',
CB_hist='array',
FM_hist='array',
SSB_hist='array',
Hist = 'Hist',
PPD = 'list',
Misc = 'list')
)
setMethod("initialize", "MSE", function(.Object, Name, nyears, proyears,
nMPs, MPs, nsim, OM, Obs,
SB_SBMSY, F_FMSY, N, B, SSB, VB, FM,
SPR, Catch, Removals, Effort, TAC, TAE,
BioEco, RefPoint,
CB_hist, FM_hist, SSB_hist,
Hist, PPD, Misc
) {
slts <- slotNames('MSE')
for (sl in slts) {
slot(.Object, sl) <- get(sl)
}
.Object
})
# ---- PMobj Class ----
#' An object for storing data for analysis using data-limited methods
#'
#' Used internally
#'
#' @name PMobj-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('PMobj')}
#' @slot Name Name of the Performance Metric. Character
#' @slot Caption A caption to be used in plots. Character, call, or function.
#' @slot Stat Statistic of interest for the PM. Dimensions: nsim, nMP, yrs. Array
#' @slot Ref Reference value to calculate probability for statistic. Numeric.
#' @slot Prob Probability (mean over years) Dimensions: nsim by MP. Matrix, numeric or data.frame
#' @slot Mean Mean probability (mean over years and simulations). Numeric. Length nMPs
#' @slot MPs Name of MPs. Single value. Character string
#' @author A. Hordyk
#' @keywords classes
setClass("PMobj", representation(Name = "character", Caption='label.class',
Stat='array', Ref='numeric', Prob='prob.class', Mean='numeric',
MPs="character"))
show <- function(object) methods::show(object)
#' Show the output of a PM
#'
#' @param object object of class MSE
#' @rdname show-MSE
#' @export
setMethod("show", signature = (object="PMobj"), function(object) {
sls <- slotNames(object)
df <- data.frame(Slot=sls, Value=NA, stringsAsFactors = FALSE)
if (any(sapply(sls, function(sl) length(slot(object, sl))) == 0)) {
cat("Incomplete PMobj\n")
for (sl in sls) {
r <- match(sl, sls)
slval <- slot(object, sl)
if ('array' %in% class(slval) & length(slval)>0 & length(dim(slval))>2) {
df[r,2] <- 'array'
} else if ('matrix' %in% class(slval) & length(slval)>0) {
df[r,2] <- 'matrix'
} else if (length(slval) > 0 & ! 'call' %in% class(slval)) {
df[r,2] <- slval
} else {
df[r, 2] <- 'not defined'
}
}
print(df)
} else {
cat(object@Name)
cat("\n", object@Caption)
cat("\n")
nMP <- length(object@MPs)
if (nMP > 1) nsim <- dim(object@Prob)[1]
if (nMP == 1) nsim <- length(object@Prob)
nprint <- min(nsim, 10)
if (nMP > 1) df <- data.frame(object@Prob[1:nprint,])
if (nMP == 1) df <- data.frame(object@Prob[1:nprint])
if (nsim > (nprint+1)) {
if (nMP > 1) lst <- object@Prob[nprint+1,]
if (nMP == 1) lst <- object@Prob[nprint+1]
} else {
if (nMP > 1) lst <- object@Prob[nprint,]
if (nMP == 1) lst <- object@Prob[nprint]
}
df <- round(df,2)
lst <- round(lst,2)
colnames(df) <- object@MPs
names(lst) <- object@MPs
if (nsim > (nprint+1)) {
df <- rbind(df,
rep(".", nMP),
rep(".", nMP),
rep(".", nMP),
lst)
rownames(df) <- c(1:(nprint+3), nsim)
}
print(df)
cat("\nMean\n")
print(round(object@Mean,2))
}
})
# ---- Summary of MSE object ----
#' Summary of MSE object
#'
#' @param object object of class MSE
#' @param ... a list of names of PM methods
#' @param silent Should summary be printed to console? Logical.
#' @param Refs An optional named list (matching the PM names) with numeric values to override the default `Ref` values. See examples.
#' @rdname summary-MSE
#' @export
setMethod('summary', signature="MSE", function(object, ..., silent=FALSE, Refs=NULL) {
PMlist <- unlist(list(...))
if(length(PMlist) == 0) PMlist <- c("PNOF", "P50", "AAVY", "LTY")
if (!methods::is(PMlist,'character')) stop("Must provide names of PM methods")
# check
for (X in seq_along(PMlist))
if (!methods::is(get(PMlist[X]), "PM")) stop(PMlist[X], " is not a valid PM method")
if (!silent) message("Calculating Performance Metrics")
storeMean <- vector('list', length(PMlist))
storeName <- vector('list', length(PMlist))
storeCap <- vector('list', length(PMlist))
storeHeading <- vector('list', length(PMlist))
storeMP <- vector('list', length(PMlist))
for (X in 1:length(PMlist)) {
ref <- Refs[[PMlist[X]]]
if (is.null(ref)) {
runPM <- eval(call(PMlist[X], object))
} else {
runPM <- eval(call(PMlist[X], object, Ref=ref))
}
storeMean[[X]] <- runPM@Mean
storeName[[X]] <- runPM@Name
storeCap[[X]] <- runPM@Caption
storeMP[[X]] <- runPM@MPs
}
df <- data.frame('MP'=storeMP[[1]], signif(do.call('cbind', storeMean),2), stringsAsFactors = FALSE)
# heading <- do.call('rbind', storeHeading)
colnames(df)[2:(length(PMlist)+1)] <- PMlist #caps # gsub(" ", "", caps)
if (!silent) {
dfprint <- data.frame('Performance Metrics' = do.call('rbind', storeName), gap="", do.call('rbind', storeCap))
names(dfprint)[2:3] <- ''
print(dfprint)
cat("\n")
cat("\nPerformance Statistics:\n")
print(df)
}
invisible(df)
})
# # ---- Plot Data Object -----
# #' Plot Data object
# #'
# #' @rdname plot-Data
# #' @param x object of class Data
# #' @param upq Upper quantile of TACs for max ylim
# #' @param lwq Lower quantile of TACs for min ylim
# #' @param outline Logical. Include outliers in plot?
# #' @param ... Optional additional arguments passed to \code{boxplot}
# #' @export
# setMethod("plot",
# signature(x = "Data"),
# function(x, upq=0.9, lwq=0.1, outline = FALSE, ...){
#
# old_par <- par(no.readonly = TRUE)
# on.exit(par(list = old_par), add = TRUE)
# boxplot.Data(x, upq, lwq, outline, ...)
# })
# # Data<-x
# if (class(Data) != "Data") stop("Must supply object of class Data")
# if (all(is.na(Data@TAC))) stop("No TAC data found")
# cols<-rep(c('black','red','green','blue','orange','brown','purple','dark grey','violet','dark red','pink','dark blue','grey'),4)
# ltys<-rep(1:4,each=13)
#
# if(is.na(funcs[1]))funcs<-Data@MPs
#
# nMPs<-length(funcs)
# nplots<-ceiling(nMPs/maxlines)
# maxl<-ceiling(nMPs/nplots)
# mbyp <- split(1:nMPs, ceiling(1:nMPs/maxl)) # assign methods to plots
#
# if(is.na(xlims[1])|length(xlims)!=2){
# xlims<-quantile(Data@TAC,c(0.005,0.95),na.rm=T)
# if(xlims[1]<0)xlims[1]<-0
# }
# if(!NAor0(Data@Ref)){
# if(xlims[1]>Data@Ref)xlims[1]<-max(0,0.98*Data@Ref)
# if(xlims[2]<Data@Ref)xlims[2]<-1.02*Data@Ref
# }
# ylims<-c(0,1)
#
# #for(m in 1:nMPs){
# # if(sum(!is.na(Data@TAC[m,,1]))>2){
# # dens<-density(Data@TAC[m,,1],na.rm=T)
# #print(quantile(dens$y,0.99,na.rm=T))
# # if(quantile(dens$y,0.9,na.rm=T)>ylims[2])ylims[2]<-quantile(dens$y,0.90,na.rm=T)
# #}
# #}
#
# #dev.new2(width=10,height=0.5+7*nplots)
# par(mfrow=c(ceiling(nplots/2),2),mai=c(0.4,0.4,0.01,0.01),omi=c(0.35,0.35,0.35,0.05))
#
# for(p in 1:nplots){
# m<-mbyp[[p]][1]
# plot(NA,NA,xlim=xlims,ylim=ylims,main="",xlab="",ylab="",col="white",lwd=3,type="l")
# abline(h=0)
# if(!NAor0(Data@Ref)){
# abline(v=Data@Ref,col="light grey",lwd=2)
# if(!NAor0(Data@Ref_type[1]))legend('right',Data@Ref_type,text.col="grey",bty='n')
# }
# #plot(density(DLM@TAC[m,,1],from=0,na.rm=T),xlim=xlims,ylim=ylims,main="",xlab="",ylab="",col=coly[m],lty=ltyy[m],type="l")
#
# if(!is.na(perc[1]))abline(v=quantile(Data@TAC[m,,1],p=perc,na.rm=T),col=cols[m],lty=ltys[m])
# #if(length(mbyp[[p]])>0){
# for(ll in 1:length(mbyp[[p]])){
# m<-mbyp[[p]][ll]
# if(sum(!is.na(Data@TAC[m,,1]))>10){ # only plot if there are sufficient non-NA TAC samples
# x<-density(Data@TAC[m,,1],from=0,na.rm=T)$x
# y<-density(Data@TAC[m,,1],from=0,na.rm=T)$y
# y<-y/max(y)
# lines(x,y,col=cols[ll])
# }else{
# print(paste("Method ",funcs[m]," produced too many NA TAC values for plotting densities",sep=""))
# }
# if(!is.na(perc[1]))abline(v=quantile(Data@TAC[m,,1],p=perc,na.rm=T),col=cols[ll],lty=2)
# }
# #}
# cind<-1:length(mbyp[[p]])
# legend('topright',funcs[mbyp[[p]]],text.col=cols[cind],col=cols[cind],lty=1,bty='n',cex=0.75)
# }
#
# mtext(paste("TAC (",Data@Units,")",sep=""),1,outer=T,line=0.5)
# mtext(paste("Standardized relative frequency",sep=""),2,outer=T,line=0.5)
# mtext(paste("TAC calculation for ",Data@Name,sep=""),3,outer=T,line=0.5)
# })
# # ---- Plot MSE object ----
# #' Plot MSE object
# #'
# #' @rdname plot-MSE
# #' @param x object of class MSE
# #' @export
# setMethod("plot",
# signature(x = "MSE"),
# function(x){
# MSEobj<-x
# Pplot(MSEobj)
# Kplot(MSEobj)
# Tplot(MSEobj)
# })
# ---- Summary of Data Object ----
#' Summary of Data object
#'
#' @rdname summary-Data
#' @param object An object of class Data
#' @param wait Logical. Wait for key press before next plot?
#' @param x iteration number for the Data object.
#' @param plots Character. What plots to show? `all`, `TS`, `CAA`, `CAL`, `PD`
#' for all plots, time-series, catch-at-age, catch-at-length, and
#' probability distributions respectively
#' @param rmd Logical. Used in a rmd file?
#' @param head Character. Heading for rmd file. Default is '##' (second level heading)
#' @param tplot Integer. Number of plots per page. Default 25
#' @export
setMethod("summary",
signature(object = "Data"),
function(object, wait=TRUE, x=1, plots='all', rmd=FALSE, head="##", tplot=25){
plots <- match.arg(plots, c('all', 'TS', 'CAA', 'CAL', 'PD'), several.ok = TRUE)
if ('all' %in% plots) plots <- c('TS', 'CAA', 'CAL', 'PD')
Freq <- n <- Var2 <- NULL # cran check
if (!methods::is(object, "Data")) stop("Object must be class `Data`", call.=FALSE)
# Time-Series
Year <- object@Year
l <- length(Year)
slts <- c("Cat", 'Ind', 'SpInd', 'VInd', 'Rec', 'ML', 'Lc')
Cat <- Ind <- SpInd <- VInd <- Rec <- ML <- Lc <- NULL # cran checks
for (sl in slts) {
tt <- slot(object, sl)[x,]
if (length(tt)!=l) tt <- c(tt, rep(NA,l-length(tt)))
assign(sl, tt)
}
Val <- c(Cat, Ind, SpInd, VInd, Rec, ML, Lc)
Var <- rep(c("Catch", "Total Index", "Spawning Index",
"Vuln. Index", "Recruitment", "Mean Length", "Mean Length above Lc"), each=length(Year))
ts.df <- data.frame(Year=object@Year, Val=Val, Var=Var, stringsAsFactors = TRUE)
# ts.df$Year <- as.factor(ts.df$Year)
ts.df$Var <- factor(ts.df$Var, levels=
c("Catch", "Total Index", "Spawning Index",
"Vuln. Index", "Recruitment", "Mean Length", "Mean Length above Lc"))
ts.df <- subset(ts.df, !is.na(Val))
if (nrow(ts.df)>0 && 'TS' %in% plots) {
P1 <- ggplot2::ggplot(ts.df, ggplot2::aes(x=Year, y=Val, group = Var)) +
ggplot2::facet_wrap(~Var, scales='free_y') + ggplot2::geom_line(linewidth=1.25) +
ggplot2::theme_classic() +
ggplot2::expand_limits(y=0) +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 1)) +
ggplot2::scale_x_continuous(breaks=pretty(rev(Year), length(Year)/5)) +
ggplot2::labs(y="")
} else {
P1 <- NULL
}
if (!is.null(P1)) {
if (rmd) {
cat('\n')
cat('\n')
cat(head, 'Time-Series')
cat('\n')
} else {
message('Plotting Time-Series')
}
print(P1)
}
if (interactive() & wait & !is.null(P1))
invisible(readline(prompt="Press [enter] to continue..."))
# Additional indices
P1a <- NULL
if(sum(dim(object@AddInd))>3) {
# additional indices exist
dd <- dim(object@AddInd[1,,])
nind <- dd[1]
nyears <- dd[2]
if (!is.null(dimnames(object@AddInd))) {
ind_names <- as.character(dimnames(object@AddInd)[[2]])
} else {
ind_names <- 1:nind
}
add.ind <- data.frame(Year=rep(Year, each=nind),
Index=ind_names,
Value=as.vector(object@AddInd[1,,]),
stringsAsFactors = TRUE)
if(nrow(add.ind)>0 && 'TS' %in% plots) {
P1a <- ggplot2::ggplot(add.ind, ggplot2::aes(x=Year, y=Value)) +
ggplot2::facet_wrap(~Index, scales='free_y') + ggplot2::geom_line(linewidth=1.25) +
ggplot2::theme_classic() +
ggplot2::expand_limits(y=0) +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 1)) +
ggplot2::scale_x_continuous(breaks=pretty(rev(Year), length(Year)/5)) +
ggplot2::labs(y="Index")
} else {
P1a <- NULL
}
if (!is.null(P1a)) {
if (rmd) {
cat('\n')
cat('\n')
cat(head, 'Additional Indices')
cat('\n')
} else {
message('Plotting Additional Indices')
}
print(P1a)
}
}
if (interactive() & wait & !is.null(P1a))
invisible(readline(prompt="Press [enter] to continue..."))
# CAA
if (all(is.na(object@CAA))) {
P2 <- FALSE
} else {
P2 <- TRUE
}
if (P2 == TRUE) {
CAA <- object@CAA[x,,]
nyrs <- nrow(CAA); maxage <- ncol(CAA)
dimnames(CAA) <- list(1:nyrs, 1:maxage)
df1 <- as.data.frame.table(CAA, stringsAsFactors = FALSE)
colnames(df1) <- c("Year", "Val", "Freq")
df1$Val <- as.numeric(df1$Val)
df1$Year <- as.numeric(df1$Year)
yr.n <- df1 %>% dplyr::group_by(Year) %>% dplyr::summarise(n=sum(Freq))
yr.ind <- yr.n %>% dplyr::filter(n>0) %>% dplyr::select(Year)
Years <- object@Year
nyears <- length(unique(df1$Year))
df1$Year_val <- (Years[(length(Years)-nyears+1):length(Years)])
df1 <- df1[!is.na(df1$Freq),] # Fliter out NA values, so we don't try to plot missing years
if (nrow(df1)>0 && 'CAA' %in% plots) {
if (rmd) {
cat('\n')
cat('\n')
cat(paste(head, 'Catch-at-Age'))
cat('\n')
} else {
message('Plotting Catch-at-Age')
}
nyears <- length(unique(df1$Year))
nbins <- length(unique(df1$Val))
if (nyears > tplot) {
npages <- ceiling(nyears/tplot)
ncol <- 5
nrow <- 5
nplot <- ncol * nrow
} else {
npages <- 1
nrow <- ceiling(sqrt(nyears))
ncol <- ceiling(nyears/nrow)
nplot <- nyears
}
pmat <- matrix(1:(nrow*ncol), nrow=nrow, ncol=ncol, byrow=TRUE)
pmat[pmat >nplot] <- NA
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
yr1 <- 1
col <- "grey"
for (pg in 1:npages) {
if(npages>1)message('Plot ', pg, ' of ', npages)
yrind <- unique(df1$Year)[yr1:(yr1+nplot-1)]
yr1 <- max(yrind) + 1
dat <- df1 %>% dplyr::filter(Year %in% yrind)
un.yrs_val <- as.numeric(unique(dat$Year_val))
un.yrs <- as.numeric(unique(dat$Year))
if (pg >1 && pg == npages) {
nplot <- nyears - (npages-1) * tplot
ncol <- ceiling(sqrt(nplot))
nrow <- ceiling(nplot/ncol)
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
pmat <- matrix(1:nplot, nrow=nrow, ncol=ncol, byrow=TRUE)
}
for (p in 1:nplot) {
pdat <- dat %>% dplyr::filter(Year==un.yrs[p])
if (nrow(pdat) > 0) {
if (p %in% pmat[nrow,]) {
barplot(pdat$Freq, names=round(pdat$Val, 2), axes=FALSE, col=col)
} else {
barplot(pdat$Freq, names=FALSE, axes=FALSE, col=col)
}
if (p %in% pmat[,1]) axis(side=2)
if (!p %in% pmat[,1]) axis(side=2, labels=TRUE)
ncount <- round(sum(pdat$Freq),0)
title(un.yrs_val[p])
text(max(pdat$Val), max(pdat$Freq), paste('n = ', ncount),
xpd=NA)
}
}
mtext(side=1, outer=TRUE, "Age", line=2, cex=1.5)
mtext(side=2, outer=TRUE, "Frequency", line=2, cex=1.5)
}
}
}
if (interactive() & wait & !is.null(P2) && P2)
invisible(readline(prompt="Press [enter] to continue..."))
# CAL
if (all(is.na(object@CAL))) {
P3 <- FALSE
} else {
P3 <- TRUE
}
if (P3 == TRUE) {
CAL <- object@CAL[x,,]
nyrs <- nrow(CAL); nbins <- length(object@CAL_bins) - 1
By <- object@CAL_bins[2] - object@CAL_bins[1]
BinsMid <- seq(object@CAL_bins[1] + 0.5*By, by=By,length.out = nbins)
dimnames(CAL) <- list(1:nyrs, BinsMid)
df1 <- as.data.frame.table(CAL, stringsAsFactors = FALSE)
colnames(df1) <- c("Year", "Val", "Freq")
df1$Val <- as.numeric(df1$Val)
df1$Year <- as.numeric(df1$Year)
yr.n <- df1 %>% dplyr::group_by(Year) %>% dplyr::summarise(n=sum(Freq, na.rm=T))
yr.ind <- yr.n %>% dplyr::filter(n>0) %>% dplyr::select(Year)
Years <- object@Year
nyears <- length(unique(df1$Year))
df1$Year_val <- (Years[(length(Years)-nyears+1):length(Years)])
if (nrow(df1) > 0 && 'CAL' %in% plots) {
if (rmd) {
cat('\n')
cat('\n')
cat(paste(head, 'Catch-at-Length'))
cat('\n')
} else {
message('Plotting Catch-at-Length')
}
nyears <- length(unique(df1$Year))
nayears <- df1 %>% group_by(Year) %>% summarize(isna=all(is.na(Freq)))
nyears <- sum(!nayears$isna)
nbins <- length(unique(df1$Val))
if (nyears > tplot) {
npages <- ceiling(nyears/tplot)
ncol <- 5
nrow <- 5
nplot <- ncol * nrow
} else {
npages <- 1
nrow <- ceiling(sqrt(nyears))
ncol <- ceiling(nyears/nrow)
nplot <- nyears
}
pmat <- matrix(1:(ncol*nrow), nrow=nrow, ncol=ncol, byrow=TRUE)
pmat[pmat>nplot] <- NA
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
yr1 <- 1
col <- "grey"
for (pg in 1:npages) {
if(npages>1)message('Plot ', pg, ' of ', npages)
if (pg ==1) {
yrind <- yr.ind$Year[1:(nplot)]
} else {
t1 <- nplot * (pg-1) +1
t2 <- t1+nplot
yrind <- yr.ind$Year[t1:t2]
}
yr1 <- max(yrind) + 1
dat <- df1 %>% dplyr::filter(Year %in% yrind)
un.yrs_val <- as.numeric(unique(dat$Year_val))
un.yrs <- as.numeric(unique(dat$Year))
if (pg >1 && pg == npages) {
nplot <- nyears - (npages-1) * tplot
ncol <- ceiling(sqrt(nplot))
nrow <- ceiling(nplot/ncol)
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
pmat <- matrix(1:(ncol*nrow), nrow=nrow, ncol=ncol, byrow=TRUE)
pmat[pmat>nplot] <- NA
}
for (p in 1:nplot) {
pdat <- dat %>% dplyr::filter(Year==un.yrs[p])
if (nrow(pdat) > 0) {
if (all(is.na(pdat$Freq))) {
} else{
if (p %in% pmat[nrow,]) {
barplot(pdat$Freq, names.arg=round(pdat$Val, 2), axes=FALSE, col=col, las=2)
} else {
barplot(pdat$Freq, names.arg=FALSE, axes=FALSE, col=col)
}
if (p %in% pmat[,1]) axis(side=2)
if (!p %in% pmat[,1]) axis(side=2, labels=TRUE)
ncount <- round(sum(pdat$Freq),0)
title(un.yrs_val[p])
text(length(unique(df1$Val)), max(pdat$Freq), paste('n = ', ncount), xpd=NA)
}
}
}
mtext(side=1, outer=TRUE, "Length", line=2, cex=1.5)
mtext(side=2, outer=TRUE, "Frequency", line=2, cex=1.5)
}
}
}
if (interactive() & wait & !is.null(P3) && P3)
invisible(readline(prompt="Press [enter] to continue..."))
# Biology & Depletion
slots<-c("Dep","Mort","FMSY_M","Dt","BMSY_B0","vbK", "vbLinf")
namey<-c("Stock depletion", "Natural Mortality rate","Ratio of FMSY to M",
"Depletion over time t","BMSY relative to unfished",
"von B. k parameter", "von B. Linf parameter")
slotsCV<-c("CV_Dep","CV_Mort","CV_FMSY_M","CV_Dt","CV_BMSY_B0","CV_vbK", "CV_vbLinf")
reps <- 5000
val <- list()
for (i in seq_along(slots)) {
mu <- slot(object, slots[i])[x]
cv <- slot(object, slotsCV[i])[x]
val[[i]] <- trlnorm(reps, mu,cv)
}
vals <- do.call("cbind", val)
colnames(vals) <- namey
df1 <- as.data.frame.table(vals, stringsAsFactors = TRUE)
df1 <- df1 %>% dplyr::filter(is.na(Freq) == FALSE)
if (nrow(df1) > 0 && 'PD' %in% plots) {
P4 <- ggplot2::ggplot(df1, ggplot2::aes(x=Freq, group=Var2)) +
ggplot2::facet_wrap(~Var2, scales="free") + ggplot2::geom_histogram(bins=30) +
ggplot2::labs(y="Frequency", x="Parameter Value") +
ggplot2::theme_classic()
} else {
P4 <- NULL
}
if (!is.null(P4)) {
if (rmd) {
cat('\n')
cat('\n')
cat(paste(head, 'Parameter Distributions'))
cat('\n')
} else {
message('Plotting Parameter Distributions')
}
print(P4)
}
})
# -- Management Recommendation Class ----
#' Class \code{'Rec'}
#'
#' An object for storing the MP recommendations
#'
#' @name Rec-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Rec')}
#' @slot TAC A numeric value with the TAC recommendation
#' @slot Effort A numeric value with the effort recommendation as a fraction of current (nyear) fishing effort
#' @slot Spatial A boolean vector of length 'nareas' specifying if area is open (1) or closed (0) to fishing
#' @slot Allocate A boolean value describing if effort should be re-allocated from close to open areas
#' @slot LR5 smallest length at 5 per cent retention - in absolute units - i.e same units as Linf and L50
#' @slot LFR smallest length at full retention - in absolute units - i.e same units as Linf and L50
#' @slot HS upper harvest slot (no retention above this) - in absolute units - i.e same units as Linf and L50
#' @slot Rmaxlen retention of the largest size class - fraction between 0 and 1
#' @slot L5 smallest length at 5 per cent selection - in absolute units - i.e same units as Linf and L50
#' @slot LFS smallest length at full selection - in absolute units - i.e same units as Linf and L50
#' @slot Vmaxlen selection of the largest size class - fraction between 0 and 1
#' @slot Fdisc fraction of discarded fish that die - fraction between 0 and 1
#' @slot DR Discard rate - the fraction of caught fish that are discarded
#' @slot Misc An empty list that can be used to store information and pass on to MPs in future
#' @author A. Hordyk
#' @keywords classes
setClass("Rec", representation(
TAC = "numeric",
Effort = "numeric",
Spatial="numeric", Allocate = "numeric",
LR5 = "numeric", LFR = "numeric", HS="numeric", Rmaxlen="numeric",
L5 = "numeric", LFS = "numeric", Vmaxlen="numeric",
Fdisc = "numeric",
DR='numeric',
Misc="list"))
#' Show the output of a single MP recommendation
#'
#' @param object object of class Rec
#' @rdname show-Rec
#' @export
setMethod("show", signature = (object="Rec"), function(object) {
Rec <- object
slots <- slotNames(Rec)
recList <- list()
perc <- 0.5
for (X in slots) { # sequence along recommendation slots
if (X == "Misc") { # convert to a list nsim by nareas
rec <- slot(Rec, X)
} else {
rec <- slot(Rec,X) # unlist(lapply(temp, slot, name=X))
}
if (X == "Spatial") { # convert to a matrix nsim by nareas
nsims <- 1
nareas <- max(2,length(rec))
rec <- matrix(rec, nareas, nsims, byrow=FALSE)
}
recList[[X]] <- rec
}
names <- c("TAC", "Effort", "LR5", "LFR", "HS", "Rmaxlen",
"L5", "LFS", 'Vmaxlen', 'Fdisc', 'DR', 'Spatial')
mat <- matrix(0, nrow=1, ncol=length(names)+nareas-1)
count <- 0
for (x in names) {
temp <- recList[[x]]
count <- count + 1
if (x!="Spatial") {
mat[,count] <- quantile(temp, probs=perc, na.rm=TRUE)
} else {
mat[,count:ncol(mat)] <- t(matrix(unlist(temp), nrow=nareas, ncol=1))
}
}
names[length(names)] <- "Area 1"
names <- c(names, paste('Area', 2:nareas))
if (perc !=0.5) names[1] <- paste0(names[1], ' (', perc, ' perc')
if (perc ==0.5) names[1] <- paste0(names[1], ' (median)')
colnames(mat) <- names
if (nrow(mat) == 1) {
mat <- as.data.frame(mat)
matout <- mat[!is.na(mat)]
names(matout) <- names[!is.na(mat)]
print(matout)
}
})
# ---- Internal show methods ----
#' @importFrom utils capture.output str
show_int <- function(object, slots_check) {
cat(paste0("S4 object of class ", dQuote(class(object)), ".\n"))
if (!missing(slots_check)) {
for(i in slots_check) {
len <- length(slot(object, i))
if (len) cat(paste0(len, " items found in slot ", dQuote(i), ".\n"))
}
}
cat("\n")
cat("Use str(), slotNames(),", dQuote("@"), "etc. to explore contents:\n\n")
txt <- capture.output(utils::str(object))
for(i in txt[1:5]) cat(i, "\n")
invisible()
}
#' @name show-MSEtool
#' @aliases show show,Data-method show,OM-method show,Hist-method show,MSE-method show,MMSE-method
#' @title Show MSEtool S4 objects
#'
#' @description Briefly prints a couple of lines from \link[utils]{str} to avoid swamping the console with
#' the contents of very large objects.
#' @param object S4 object from MSEtool
setMethod("show", "Data", function(object) show_int(object))
#' @rdname show-MSEtool
setMethod("show", "OM", function(object) show_int(object, slots_check = "cpars"))
#' @rdname show-MSEtool
setMethod("show", "Hist", function(object) show_int(object))
#' @rdname show-MSEtool
setMethod("show", "MSE", function(object) show_int(object))
Try the MSEtool package in your browser
Any scripts or data that you put into this service are public.
MSEtool documentation built on July 26, 2023, 5:21 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions show_int show importslot
Documented in show
#' Label class union for performance metric objects
#'
#' @description Used internally. Nothing to see here!
#' @keywords internal
#' @export
setClassUnion(name="label.class", members=c("call", "character", "function"))
#' Prob class union for performance metric objects
#'
#' @description Used internally. Nothing to see here!
#' @keywords internal
#' @export
setClassUnion(name="prob.class", members=c("matrix", "numeric", "data.frame"))
# ---- Data Class ----
#' Class \code{'Data'}
#'
#' An object for storing fishery data for analysis
#'
#'
#' @name Data-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Data', stock)}
#' @slot Name The name of the Data object. Single value. Character string
#' @slot Common_Name Common name of the species. Character string
#' @slot Species Scientific name of the species. Genus and species name. Character string
#' @slot Region Name of the general geographic region of the fishery. Character string
#' @slot LHYear The last historical year of the simulation (before projection). Single value. Positive integer
#' @slot MPrec The previous recommendation of a management procedure. Vector of length nsim. Positive real numbers
#' @slot Units Units of the catch/absolute abundance estimates. Single value. Character string
#' @slot MPeff The current level of effort. Vector of length nsim. Positive real numbers
#' @slot nareas Number of fishing areas. Vector of length nsim. Non-negative integer
#'
#' @slot MaxAge Maximum age. Vector nsim long. Positive integer
#' @slot Mort Natural mortality rate. Vector nsim long. Positive real numbers
#' @slot CV_Mort Coefficient of variation in natural mortality rate. Vector nsim long. Positive real numbers
#' @slot vbLinf Maximum length. Vector nsim long. Positive real numbers
#' @slot CV_vbLinf Coefficient of variation in maximum length. Vector nsim long. Positive real numbers
#' @slot vbK The von Bertalanffy growth coefficient K. Vector nsim long. Positive real numbers
#' @slot CV_vbK Coefficient of variation in the von Bertalanffy K parameter. Vector nsim long. Positive real numbers
#' @slot vbt0 Theoretical age at length zero. Vector nsim long. Non-positive real numbers
#' @slot CV_vbt0 Coefficient of variation in age at length zero. Vector nsim long. Positive real numbers
#' @slot wla Weight-Length parameter alpha. Vector nsim long. Positive real numbers
#' @slot CV_wla Coefficient of variation in weight-length parameter a. Vector nsim long. Positive real numbers
#' @slot wlb Weight-Length parameter beta. Vector nsim long. Positive real numbers
#' @slot CV_wlb Coefficient of variation in weight-length parameter b. Vector nsim long. Positive real numbers
#' @slot steep Steepness of stock-recruitment relationship. Vector nsim long. Value in the range of one-fifth to 1
#' @slot CV_steep Coefficient of variation in steepness. Vector nsim long. Positive real numbers
#' @slot sigmaR Recruitment variability. Vector nsim long. Positive real numbers
#' @slot CV_sigmaR Coefficient of variation in recruitment variability. Vector nsim long. Positive real numbers
#' @slot L50 Length at 50 percent maturity. Vector nsim long. Positive real numbers
#' @slot CV_L50 Coefficient of variation in length at 50 per cent maturity. Vector nsim long. Positive real numbers
#' @slot L95 Length at 95 percent maturity. Vector nsim long. Positive real numbers
#' @slot LenCV Coefficient of variation of length-at-age (assumed constant for all age classes). Vector nsim long. Positive real numbers
#' @slot LFC Length at first capture. Vector nsim long. Positive real numbers
#' @slot CV_LFC Coefficient of variation in length at first capture. Vector nsim long. Positive real numbers
#' @slot LFS Shortest length at full selection. Vector nsim long. Positive real numbers
#' @slot CV_LFS Coefficient of variation in length at full selection. Vector nsim long. Positive real numbers
#' @slot Vmaxlen Vulnerability of individuals at asymptotic length. Vector nsim long. Real number between 0 and 1.
#'
#' @slot Year Years that corresponding to catch and relative abundance data. Vector nyears long. Positive integer
#' @slot Cat Total annual catches. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Cat Coefficient of variation in annual catches. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Cat` for all years.
#' @slot Effort Annual fishing effort. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Effort Coefficient of variation in annual effort. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Effort` for all years.
#' @slot Ind Relative total abundance index. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Ind Coefficient of variation in the relative total abundance index. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Ind` for all years
#'
#' @slot SpInd Relative spawning abundance index. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_SpInd Coefficient of variation in the relative spawning abundance index. Matrix nsim rows and either 1 or nyear columns. Positive real numbers.
#'
#' @slot VInd Relative vulnerable abundance index. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_VInd Coefficient of variation in the relative vulnerable abundance index. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers.
#'
#' @slot AddInd Optional additional indices. Array of dimensions `nsim`, n additional indices, and `nyears` (length `Year`).
#' @slot CV_AddInd Coefficient of variation for additional indices. Array of same dimensions as `AddInd`
#' @slot AddIndV Vulnerability-at-age schedules for the additional indices. Array with dimensions: `nsim`, n additional indices, `MaxAge+1`.
#' @slot AddIunits Units for the additional indices - biomass (1; default) or numbers (0). Numeric vector length n.ind.
#' @slot AddIndType Index calculated from total stock (1, default), spawning stock (2), or vulnerable stock (3). Numeric vector of length n.ind
#'
#' @slot Rec Recent recruitment strength. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot CV_Rec Log-normal CV for recent recruitment strength. Matrix nsim rows and either 1 or nyear columns.
#' Positive real numbers. Note: built-in MPs use only the first value of `CV_Rec` for all years.
#' @slot ML Mean length time series. Matrix of nsim rows and nyears columns. Non-negative real numbers
#' @slot Lc Modal length of catches. Matrix of nsim rows and nyears columns. Positive real numbers
#' @slot Lbar Mean length of catches over Lc. Matrix of nsim rows and nyears columns. Positive real numbers
#' @slot Vuln_CAA Optional vulnerability-at-age schedule for catch-at-age samples. Used to condition OM for closed-loop
#' simulation testing. Replaces the fleet selectivity schedule in the OM used to generate CAA samples. Matrix
#' with dimensions `nsim` x `MaxAge+1`.
#' @slot CAA Catch at Age data (numbers). Array of dimensions nsim x nyears x MaxAge+1. Non-negative integers
#' @slot Vuln_CAL Optional vulnerability-at-length schedule for catch-at-length samples. Used to condition OM for closed-loop
#' simulation testing. Replaces the fleet selectivity schedule in the OM used to generate CAL samples. Matrix
#' with dimensions `nsim` x `length(CAL_mids)`.
#' @slot CAL_bins The values delimiting the length bins for the catch-at-length data. Vector. Non-negative real numbers
#' @slot CAL_mids The values of the mid-points of the length bins. Optional, calculated from `CAL_bins` if not entered. Vector. Non-negative real numbers.
#' @slot CAL Catch-at-length data. An array with dimensions nsim x nyears x length(CAL_mids). Non-negative integers. By default the CAL data will be the retained lengths (i.e, not including discards). If `OM@control$CAL =="removals"` then the CAL data will include all removals (retained + discards).
#'
#' @slot Dep Stock depletion SSB(current)/SSB(unfished). Vector nsim long. Fraction.
#' @slot CV_Dep Coefficient of variation in current stock depletion. Vector nsim long. Positive real numbers
#' @slot Abun An estimate of absolute current vulnerable abundance. Vector nsim long. Positive real numbers
#' @slot CV_Abun Coefficient of variation in estimate of absolute current stock size. Vector nsim long. Positive real numbers
#' @slot SpAbun An estimate of absolute current spawning stock abundance. Vector nsim long. Positive real numbers
#' @slot CV_SpAbun Coefficient of variation in estimate of absolute spawning current stock size. Vector nsim long. Positive real numbers
#'
#' @slot FMSY_M An assumed ratio of FMSY to M. Vector nsim long. Positive real numbers
#' @slot CV_FMSY_M Coefficient of variation in the ratio in FMSY/M. Vector nsim long. Positive real numbers
#' @slot BMSY_B0 The most productive stock size relative to unfished. Vector nsim long. Fraction
#' @slot CV_BMSY_B0 Coefficient of variation in the position of the most productive stock size relative to unfished. Vector nsim long. Positive real numbers
#' @slot Cref Reference or target catch level (eg MSY). Vector of length nsim. Positive real numbers
#' @slot CV_Cref Log-normal CV for reference or target catch level. Vector of length nsim. Positive real numbers
#' @slot Bref Reference or target biomass level (eg BMSY). Vector of length nsim. Positive real numbers
#' @slot CV_Bref Log-normal CV for reference or target biomass level. Vector of length nsim. Positive real numbers
#' @slot Iref Reference or target relative abundance index level (eg BMSY / B0). Vector of length nsim. Positive real numbers
#' @slot CV_Iref Log-normalCV for reference or target relative abundance index level. Vector of length nsim. Positive real numbers
#' @slot t The number of years corresponding to AvC and Dt. Single value. Positive integer
#' @slot AvC Average catch over time t. Vector nsim long. Positive real numbers
#' @slot CV_AvC Coefficient of variation in average catches over time t. Vector nsim long. Positive real numbers
#' @slot Dt Depletion over time t SSB(now)/SSB(now-t+1). Vector nsim long. Fraction
#' @slot CV_Dt Coefficient of variation in depletion over time t. Vector nsim long. Positive real numbers
#' @slot Ref A reference management level (eg a catch limit). Single value. Positive real number
#' @slot Ref_type Type of reference management level (eg 2009 catch limit). Single value. Character string
#' @slot Log A record of events. Single value. Character string
#' @slot params A place to store estimated parameters. An object. R list
#' @slot PosMPs The methods that can be applied to these data. Vector. Character strings
#' @slot TAC The calculated catch limits (function TAC). An array with dimensions PosMPs x replicate TAC samples x nsim. Positive real numbers
#' @slot Sense The results of the sensitivity analysis (function Sense). An array with dimensions PosMPs x sensitivity increments. Positive real numbers
#' @slot MPs The methods that were applied to these data. Vector. Character strings
#' @slot OM A table of operating model conditions. R table object of nsim rows. Real numbers
#' @slot Obs A table of observation model conditions. R table object of nsim rows. Real numbers
#' @slot Misc Other information for MPs. An object. R list
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' newdata<-new('Data')
#'
setClass("Data",
representation(Name = "character", Common_Name='character',
Species='character', Region='character',
LHYear = "numeric", MPrec = "vector",
Units = "character", MPeff = "vector",
nareas = "numeric",
MaxAge = "vector", Mort = "vector", CV_Mort = "vector",
vbLinf = "vector", CV_vbLinf = "vector",
vbK = "vector", CV_vbK = "vector",
vbt0 = "vector", CV_vbt0 = "vector",
wla = "vector", CV_wla = "vector",
wlb = "vector", CV_wlb = "vector",
steep = "vector", CV_steep = "vector",
sigmaR='vector', CV_sigmaR='vector',
L50 = "vector", CV_L50 = "vector",
L95 = "vector",
LenCV="vector",
LFC = "vector", CV_LFC = "vector",
LFS = "vector", CV_LFS = "vector",
Vmaxlen = 'vector',
Year = "vector",
Cat = "matrix", CV_Cat = "matrix",
Effort = 'matrix', CV_Effort = 'matrix',
Ind = "matrix", CV_Ind = "matrix",
SpInd = "matrix", CV_SpInd = "matrix",
VInd = "matrix", CV_VInd = "matrix",
AddInd = "array", CV_AddInd = "array", AddIndV = "array",
AddIunits = 'vector', AddIndType='vector',
Rec = "matrix", CV_Rec = "matrix",
ML = "matrix", Lc = "matrix", Lbar = "matrix",
Vuln_CAA = "matrix", CAA = "array",
Vuln_CAL = 'matrix', CAL_bins = "numeric",
CAL_mids = "numeric", CAL = "array",
Dep = "vector", CV_Dep = "vector",
Abun = "vector", CV_Abun = "vector",
SpAbun= "vector", CV_SpAbun = "vector",
FMSY_M = "vector", CV_FMSY_M = "vector",
BMSY_B0 = "vector", CV_BMSY_B0 = "vector",
Cref = "vector", CV_Cref = "vector",
Bref = "vector", CV_Bref = "vector",
Iref = "vector", CV_Iref = "vector",
t = "vector", AvC = "vector", CV_AvC = "vector",
Dt = "vector", CV_Dt = "vector",
Ref = "numeric", Ref_type = "character",
Log = "list", params = "list", PosMPs = "vector",
TAC = "array", Sense = "array",
MPs = "vector", OM = "data.frame", Obs = "data.frame",
Misc = "list"))
setMethod("initialize", "Data", function(.Object, stock="nada", ...) {
if (file.exists(stock)) {
.Object <- XL2Data(stock, ...)
} else {
slots <- slotNames('Data')
for (x in seq_along(slots)) {
sl <- slots[x]
cl <- class(slot(.Object, sl))
if ("logical" %in% cl) slot(.Object, sl) <- as.numeric(NA)
if ("character" %in% cl) slot(.Object, sl) <- ''
if ("matrix" %in% cl) {
if (length(dim(slot(.Object, sl))) > 2) {
slot(.Object, sl) <- array(NA, dim=c(1,1,1))
} else {
slot(.Object, sl) <- matrix(NA)
}
}
if ("array" %in% cl) {
if (length(dim(slot(.Object, sl))) > 2) {
slot(.Object, sl) <- array(NA, dim=c(1,1,1))
} else {
slot(.Object, sl) <- matrix(NA)
}
}
if ("vector" %in% cl) slot(.Object, sl) <- NA
if ("numeric" %in% cl) slot(.Object, sl) <- as.numeric(NA)
if ("list" %in% cl) slot(.Object, sl) <- list()
}
}
# Default values
if (all(is.na(.Object@CV_Cat))) .Object@CV_Cat <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_Ind))) .Object@CV_Ind <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_SpInd))) .Object@CV_SpInd <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_VInd))) .Object@CV_VInd <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_Effort))) .Object@CV_Effort <- matrix(0.2, nrow=1, ncol=1)
if (all(is.na(.Object@CV_Rec))) .Object@CV_Rec <- matrix(0.2, nrow=1, ncol=1)
if (NAor0(.Object@LenCV)) .Object@LenCV <- 0.1
if (NAor0(.Object@CV_Dt)) .Object@CV_Dt <- 0.25
if (NAor0(.Object@CV_AvC)) .Object@CV_AvC <- 0.2
if (NAor0(.Object@CV_Mort)) .Object@CV_Mort <- 0.2
if (NAor0(.Object@CV_FMSY_M)) .Object@CV_FMSY_M <- 0.2
if (NAor0(.Object@CV_BMSY_B0)) .Object@CV_BMSY_B0 <- 0.045
if (NAor0(.Object@CV_Cref)) .Object@CV_Cref <- 0.2
if (NAor0(.Object@CV_Bref)) .Object@CV_Bref <- 0.2
if (NAor0(.Object@CV_Iref)) .Object@CV_Iref <- 0.2
if (NAor0(.Object@CV_Dep)) .Object@CV_Dep <- 0.25
if (NAor0(.Object@CV_Abun)) .Object@CV_Abun <- 0.25
if (NAor0(.Object@CV_vbK)) .Object@CV_vbK <- 0.1
if (NAor0(.Object@CV_vbLinf)) .Object@CV_vbLinf <- 0.1
if (NAor0(.Object@CV_vbt0)) .Object@CV_vbt0 <- 0.1
if (NAor0(.Object@CV_L50)) .Object@CV_L50 <- 0.1
if (NAor0(.Object@CV_LFC)) .Object@CV_LFC <- 0.2
if (NAor0(.Object@CV_LFS)) .Object@CV_LFS <- 0.2
if (NAor0(.Object@CV_wla)) .Object@CV_wla <- 0.1
if (NAor0(.Object@CV_wlb)) .Object@CV_wlb <- 0.1
if (NAor0(.Object@CV_steep)) .Object@CV_steep <- 0.2
if (NAor0(.Object@nareas)) .Object@nareas <- 2
if (all(is.na(.Object@CAL))) .Object@CAA <- array(NA, c(1, 1, 1))
if (all(is.na(.Object@CAL))) .Object@CAL <- array(NA, c(1, 1, 1))
if (length(.Object@CAL_bins) == 0) .Object@CAL_bins <- 1
if (all(is.na(.Object@CAL))) .Object@AddInd <- array(NA, c(1, 1, 1))
if (all(is.na(.Object@CAL))) .Object@CV_AddInd <- array(NA, c(1, 1, 1))
if (all(is.na(.Object@CAL))) .Object@AddIndV <- array(NA, c(1, 1, 1))
if (length(.Object@TAC) == 0) .Object@TAC <- array(1, c(1, 1))
# if (length(.Object@TACbias) == 0) .Object@TACbias <- array(1, c(1, 1))
if (length(.Object@Sense) == 0) .Object@Sense <- array(1, c(1, 1))
if (length(.Object@ML) == 0) .Object@ML <- array(NA, c(1, 1))
if (length(.Object@Lbar) == 0) .Object@Lbar <- array(NA, c(1, 1))
if (length(.Object@Lc) == 0) .Object@Lc <- array(NA, c(1, 1))
return(.Object)
})
importslot <- function(name, length=2, Data, Names, numeric=TRUE, essential=TRUE) {
x <- Data[match(name, Names), 1:length]
if (numeric) x <- as.numeric(x)
if (essential) {
if (any(is.na(x))) warning('NAs in ', name, '. Should be length ', length)
if (length(x)<length) warning(name, ' is missing')
}
x
}
# ---- Stock Class ----
#' Class \code{'Stock'}
#'
#' An operating model component that specifies the parameters of the population
#' dynamics model
#'
#' @docType class
#'
#' @template Stock_template
#'
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Stock')}
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' showClass('Stock')
#'
setClass("Stock", representation(Name = "character",
Common_Name='character',
Species="character",
maxage = "numeric",
R0 = "numeric",
M = "numeric",
Msd = "numeric",
h = "numeric",
SRrel = "numeric",
Perr = "numeric",
AC = "numeric",
Linf = "numeric",
Linfsd = "numeric",
K = "numeric",
Ksd = "numeric",
t0 = "numeric",
LenCV="numeric",
L50 = "numeric",
L50_95 = "numeric",
D = "numeric",
a = "numeric",
b = "numeric",
Size_area_1 = "numeric",
Frac_area_1 = "numeric",
Prob_staying = "numeric",
Fdisc="numeric",
Source = "character"))
setMethod("initialize", "Stock", function(.Object, file = NA, dec=c(".", ",")) {
if (!is.na(file)) {
if (!file.exists(file)) {
stop('File not found')
} else {
# load file and populate object
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec == ".") dat <- read.csv(file, header = F,
colClasses = "character",
col.names = paste0("V", 1:Ncol))
if (dec == ",") dat <- read.csv2(file, header = F,
colClasses = "character",
col.names = paste0("V", 1:Ncol))
Names <- dat[,1]
Ncol <- ncol(dat)
Data <- dat[,2:Ncol]
.Object@Name <- importslot('Name', 1, Data, Names, FALSE)
commonname <- Data[match("Common_Name", Names), 1]
commonname1 <- Data[match("Common Name", Names), 1]
if (length(commonname)>0 && !is.na(commonname)) {
.Object@Common_Name <-commonname
} else if (length(commonname1)>0 && !is.na(commonname)) {
.Object@Common_Name <-commonname1
} else {
.Object@Common_Name <- "Not specified"
}
.Object@Species <- importslot('Species', 1, Data, Names, FALSE, FALSE)
.Object@maxage <- importslot('maxage', 1, Data, Names)
.Object@R0 <- importslot('R0', 1, Data, Names)
.Object@M <- importslot('M', 2, Data, Names)
.Object@Msd <- importslot('Msd', 2, Data, Names)
.Object@Fdisc <- importslot('Fdisc', 2, Data, Names)
.Object@h <- importslot('h', 2, Data, Names)
.Object@SRrel <- importslot('SRrel', 1, Data, Names)
.Object@Linf <- importslot('Linf', 2, Data, Names)
.Object@K <- importslot("K", 2, Data, Names)
.Object@t0 <- importslot("t0", 2, Data, Names)
.Object@LenCV <- importslot("LenCV", 2, Data, Names)
.Object@Ksd <- importslot("Ksd", 2, Data, Names)
.Object@Linfsd <- importslot("Linfsd", 2, Data, Names)
.Object@a <- importslot("a", 1, Data, Names)
.Object@b <- importslot("b", 1, Data, Names)
.Object@D <- importslot("D", 2, Data, Names)
.Object@Perr <- importslot("Perr", 2, Data, Names)
.Object@AC <- importslot("AC", 2, Data, Names)
.Object@Size_area_1 <- importslot("Size_area_1", 2, Data, Names)
.Object@Frac_area_1 <- importslot("Frac_area_1", 2, Data, Names)
.Object@Prob_staying <- importslot("Prob_staying", 2, Data, Names)
.Object@L50 <- importslot("L50", 2, Data, Names)
.Object@L50_95 <- importslot("L50_95", 2, Data, Names)
.Object@Source <- importslot("Source", 1, Data, Names, FALSE, FALSE)
}
}
.Object
})
#' Class union for isRel slot
#'
#' @description Used internally.
#' @keywords internal
#' @export
setClassUnion(name="char.log", members=c("character", "logical"))
# ---- Fleet Class -----
#' Class \code{'Fleet'}
#'
#' The component of the operating model that controls fishing dynamics
#'
#' @name Fleet-class
#' @docType class
#' @template Fleet_template
#'
#' @section Creating Object:
#' Objects can be created by calls of the form \code{new('Fleet')}
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#' showClass('Fleet')
#'
setClass("Fleet", slots = c(Name = "character",
nyears = "numeric",
CurrentYr="numeric",
EffYears = "numeric",
EffLower = "numeric",
EffUpper = "numeric",
Esd = "numeric",
qinc = "numeric",
qcv = "numeric",
L5 = "numeric",
LFS = "numeric",
Vmaxlen = "numeric",
isRel = "char.log",
LR5 = "numeric",
LFR = "numeric",
Rmaxlen = "numeric",
DR = "numeric",
Spat_targ = "numeric",
MPA='char.log',
Misc='list')
)
# initialize Fleet
setMethod("initialize", "Fleet", function(.Object, file = NA, dec=c(".", ",")) {
if (!is.na(file)) {
if (!file.exists(file)) {
stop('File not found')
} else {
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec == ".") dat <- read.csv(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol))
if (dec == ",") dat <- read.csv2(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol))
Names <- dat[, 1]
Data <- dat[, 2:ncol(dat)]
.Object@Name <- importslot('Name', 1, Data, Names, FALSE)
.Object@nyears <- importslot('nyears', 1, Data, Names)
.Object@CurrentYr <- importslot('CurrentYr', 1, Data, Names, TRUE, FALSE)
if(is.na(.Object@CurrentYr)).Object@CurrentYr<-.Object@nyears
.Object@Spat_targ <- importslot('Spat_targ', 2, Data, Names)
if (!is.na(match("Esd", Names))) {
.Object@Esd <- importslot('Esd', 2, Data, Names)
} else {
.Object@Esd <- importslot('Fsd', 2, Data, Names)
}
nEffYears <- ncol(Data[match("EffYears", Names), ])
oldw <- getOption("warn")
options(warn = -1)
chk <- as.numeric(Data[match("EffYears", Names), 1:nEffYears])
options(warn = oldw)
ind <- which(!is.na(chk))
nEffYears <- length(ind)
.Object@EffYears <- importslot('EffYears', nEffYears, Data, Names)
.Object@EffLower <- importslot('EffLower', nEffYears, Data, Names)
.Object@EffUpper <- importslot('EffUpper', nEffYears, Data, Names)
.Object@qinc <- importslot("qinc", 2, Data, Names)
.Object@qcv <- importslot("qcv", 2, Data, Names)
.Object@L5 <- importslot("L5", 2, Data, Names)
.Object@LFS <- importslot("LFS", 2, Data, Names)
.Object@Vmaxlen <- importslot("Vmaxlen", 2, Data, Names)
.Object@LR5 <- importslot("LR5", 2, Data, Names, essential = FALSE)
if (all(is.na(.Object@LR5))) .Object@LR5 <- c(0,0)
.Object@LFR <- importslot("LFR", 2, Data, Names, essential = FALSE)
if (all(is.na(.Object@LFR))) .Object@LFR <- c(0.001,0.001)
.Object@Rmaxlen <- importslot("Rmaxlen", 2, Data, Names, essential=FALSE)
if (all(is.na(.Object@Rmaxlen))) .Object@Rmaxlen <- c(1,1)
.Object@DR <- importslot("DR", 2, Data, Names)
.Object@isRel <- importslot("isRel", 1, Data, Names, FALSE, FALSE)
if (NAor0(.Object@isRel)) .Object@isRel <- "TRUE"
.Object@isRel <- as.character(.Object@isRel)
.Object@MPA <- importslot("MPA", 1, Data, Names, FALSE, FALSE)
if (NAor0(.Object@MPA)) .Object@MPA <- FALSE
.Object@MPA <- as.logical(.Object@MPA)
}
}
if (NAor0(.Object@MPA)) .Object@MPA <- FALSE
.Object
})
#' ~~ Methods for Function \code{initialize} ~~
#'
#' @name initialize-methods
#' @aliases initialize-methods initialize,Data-method
#' initialize,Fleet-method initialize,MSE-method initialize,Obs-method
#' initialize,OM-method initialize,Stock-method initialize
#' initialize,Fease-method initialize,DLM_general-method
#' @docType methods
#' @section Methods: \describe{
#'
#' \item{list('signature(.Object = \'DLM\')')}{ %% ~~describe this method
#' here~~ }
#'
#' \item{list('signature(.Object = \'Fleet\')')}{ %% ~~describe this method
#' here~~ }
#'
#' \item{list('signature(.Object = \'MSE\')')}{ %% ~~describe this method
#' here~~ }
#'
#' \item{list('signature(.Object = \'Obs\')')}{ %% ~~describe this
#' method here~~ }
#'
#' \item{list('signature(.Object = \'OM\')')}{ %% ~~describe this method here~~
#' }
#'
#' \item{list('signature(.Object = \'Stock\')')}{ %% ~~describe this method
#' here~~ }
#'
#'
#' \item{list('signature(.Object = \'Fease\')')}{ %% ~~describe this method
#' here~~ } \item{list('signature(.Object = \'DLM_general\')')}{ %% ~~describe
#' this method here~~ }
#'
#' }
#' @keywords methods ~~ other possible keyword(s) ~~
NULL
# ---- Obs Class ----
#' Class \code{'Obs'}
#'
#' An operating model component that controls the observation model
#'
#'
#' @docType class
#' @note Its questionable whether the hyperstability/hyperdepletion should be
#' categorised as an observation model characteristic as it is most often
#' driven by fleet dynamics (and therefore should be in the fleet object). Oh
#' well its here and you might want to make it hyperstable beta < 1 or
#' hyperdeplete beta > 1, only.
#'
#' @slot Name The name of the observation model object. Single value. Character string.
#' @template Obs_template
#'
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Obs')}
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' showClass('Obs')
#'
setClass("Obs", representation(Name = "character",
Cobs = "numeric",
Cbiascv = "numeric",
CAA_nsamp = "numeric",
CAA_ESS = "numeric",
CAL_nsamp = "numeric",
CAL_ESS = "numeric",
Iobs = "numeric",
Btobs = "numeric",
Btbiascv = "numeric",
beta = "numeric",
LenMbiascv = "numeric",
Mbiascv = "numeric",
Kbiascv = "numeric",
t0biascv = "numeric",
Linfbiascv = "numeric",
LFCbiascv = "numeric",
LFSbiascv = "numeric",
FMSY_Mbiascv = "numeric",
BMSY_B0biascv = "numeric",
Irefbiascv = "numeric",
Brefbiascv = "numeric",
Crefbiascv = "numeric",
Dbiascv = "numeric",
Dobs = "numeric",
hbiascv = "numeric",
Recbiascv = "numeric",
sigmaRbiascv = 'numeric',
Eobs='numeric',
Ebiascv='numeric'))
# initialize Obs
setMethod("initialize", "Obs", function(.Object, file = NA, dec=c(".", ",")) {
if (!is.na(file)) {
if (file.exists(file)) {
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec ==".") dat <- read.csv(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
if (dec ==",") dat <- read.csv2(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
dname <- dat[, 1]
dat <- dat[, 2:ncol(dat)]
.Object@Name <- dat[match("Name", dname), 1]
.Object@Cobs <- as.numeric(dat[match("Cobs", dname), 1:2])
.Object@Cbiascv <- as.numeric(dat[match("Cbiascv", dname), 1])
.Object@CAA_nsamp <- as.numeric(dat[match("CAA_nsamp", dname), 1:2])
.Object@CAA_ESS <- as.numeric(dat[match("CAA_ESS", dname), 1:2])
.Object@CAL_nsamp <- as.numeric(dat[match("CAL_nsamp", dname), 1:2])
.Object@CAL_ESS <- as.numeric(dat[match("CAL_ESS", dname), 1:2])
.Object@Iobs <- as.numeric(dat[match("Iobs", dname), 1:2])
.Object@Btobs <- as.numeric(dat[match("Btobs", dname), 1:2])
.Object@Btbiascv <- as.numeric(dat[match("Btbiascv", dname), 1])
.Object@beta <- as.numeric(dat[match("beta", dname), 1:2])
.Object@LenMbiascv <- as.numeric(dat[match("LenMbiascv", dname), 1])
.Object@Mbiascv <- as.numeric(dat[match("Mbiascv", dname), 1])
.Object@Kbiascv <- as.numeric(dat[match("Kbiascv", dname), 1])
.Object@t0biascv <- as.numeric(dat[match("t0biascv", dname), 1])
.Object@Linfbiascv <- as.numeric(dat[match("Linfbiascv", dname), 1])
.Object@LFCbiascv <- as.numeric(dat[match("LFCbiascv", dname), 1])
.Object@LFSbiascv <- as.numeric(dat[match("LFSbiascv", dname), 1])
.Object@FMSY_Mbiascv <- as.numeric(dat[match("FMSY_Mbiascv", dname), 1])
.Object@BMSY_B0biascv <- as.numeric(dat[match("BMSY_B0biascv", dname), 1])
.Object@Irefbiascv <- as.numeric(dat[match("Irefbiascv", dname), 1])
.Object@Crefbiascv <- as.numeric(dat[match("Crefbiascv", dname), 1])
.Object@Brefbiascv <- as.numeric(dat[match("Brefbiascv", dname), 1])
.Object@Dbiascv <- as.numeric(dat[match("Dbiascv", dname), 1])
.Object@Dobs <- as.numeric(dat[match("Dobs", dname), 1:2])
.Object@hbiascv <- as.numeric(dat[match("hbiascv", dname), 1])
.Object@Recbiascv <- as.numeric(dat[match("Recbiascv", dname), 1:2])
.Object@sigmaRbiascv <- as.numeric(dat[match("Recbiascv", dname), 1])
.Object@Eobs <- as.numeric(dat[match("Eobs", dname), 1:2])
.Object@Ebiascv <- as.numeric(dat[match("Ebiascv", dname), 1])
} else {
message("File doesn't exist")
}
}
.Object
})
# ---- Imp Class ----
#' Class \code{'Imp'}
#'
#' An operating model component that specifies the degree of adherence to management recommendations (Implementation error)
#'
#' @name Imp-class
#' @docType class
#' @slot Name The name of the Implementation error object. Single value. Character string.
#'
#' @template Imp_template
#'
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Imp')}#'
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#' @examples
#'
#' showClass('Imp')
#'
setClass("Imp", representation(Name = "character",
TACFrac = "numeric",
TACSD = "numeric",
TAEFrac = "numeric",
TAESD = "numeric",
SizeLimFrac="numeric",
SizeLimSD = "numeric"))
# initialize Imp
setMethod("initialize", "Imp", function(.Object, file = NA, dec=c(".", ",")) {
.Object@Name <- "Perfect implementation"
.Object@TACSD <- c(0,0)
.Object@TACFrac <- c(1,1)
.Object@TAESD <- c(0,0)
.Object@TAEFrac <-c(1,1)
.Object@SizeLimSD <- c(0,0)
.Object@SizeLimFrac<-c(1,1)
# .Object@Source <-"DLMtool generated"
if (!is.na(file)) {
if (file.exists(file)) {
dec <- match.arg(dec)
Ncol <- max(unlist(lapply(strsplit(readLines(file), ","), length)))
if (dec ==".") dat <- read.csv(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
if (dec ==",") dat <- read.csv2(file, header = F, colClasses = "character",
col.names = paste0("V", 1:Ncol)) # read 1st sheet
dname <- dat[, 1]
dat <- dat[, 2:ncol(dat)]
.Object@Name <- dat[match("Name", dname), 1]
.Object@TACSD <- as.numeric(dat[match("TACSD", dname), 1:2])
.Object@TACFrac <- as.numeric(dat[match("TACFrac", dname), 1:2])
.Object@TAESD <- as.numeric(dat[match("TAESD", dname), 1:2])
.Object@TAEFrac <- as.numeric(dat[match("TAEFrac", dname), 1:2])
.Object@SizeLimSD <- as.numeric(dat[match("SizeLimSD", dname), 1:2])
.Object@SizeLimFrac <- as.numeric(dat[match("SizeLimFrac", dname), 1:2])
# .Object@Source <- dat[match("Source", dname), 1]
} else {
message("File doesn't exist")
}
}
.Object
})
# ---- OM Class ----
#' Class \code{'OM'}
#'
#' An object containing all the parameters needed to control the MSE which can
#' be build from component Stock, Fleet, Obs, and Imp objects.
#'
#' Almost all of these inputs are a vector of length 2 which describes the upper and lower
#' bounds of a uniform distribution from which to sample the parameter.
#'
#' @name OM-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('OM', Stock, Fleet, Obs, Imp)}.
#' @slot Name Name of the operating model
#' @slot Agency Name of the agency responsible for the management of the fishery. Character string
#' @slot Region Name of the general geographic region of the fishery. Character string
#' @slot Sponsor Name of the organization who sponsored the OM. Character string
#' @slot Latitude Latitude (decimal degrees). Negative values represent the South of the Equator. Numeric. Single value
#' @slot Longitude Longitude (decimal degrees). Negative values represent the West of the Prime Meridian. Numeric. Single value
#' @slot nsim The number of simulations
#' @slot proyears The number of projected years
#' @slot interval The assessment interval - how often would you like to update the management system?
#' @slot pstar The percentile of the sample of the management recommendation for each method
#' @slot maxF Maximum instantaneous fishing mortality rate that may be simulated for any given age class
#' @slot reps Number of samples of the management recommendation for each method. Note that when this is set to 1, the mean value of
#' the data inputs is used.
#' @slot cpars A list of custom parameters. Time series are a matrix nsim rows by nyears columns. Single parameters are a vector nsim long
#' @slot seed A random seed to ensure users can reproduce results exactly
#' @slot Source A reference to a website or article from which parameters were taken to define the operating model
#'
# #' @template Stock_template
# #' @template Fleet_template
# #' @template Obs_template
# #' @template Imp_template
#'
#' @author T. Carruthers and A. Hordyk
#' @export
#' @keywords classes
#'
setClass("OM", representation(Name = "character", Agency="character",
Region="character", Sponsor="character",
Latitude="numeric", Longitude="numeric",
nsim="numeric", proyears="numeric",
interval='numeric', pstar='numeric', maxF='numeric', reps='numeric',
cpars="list",seed="numeric", Source="character"),
contains=c("Stock", "Fleet", "Obs", "Imp"))
# initialize OM
setMethod("initialize", "OM", function(.Object, Stock=NULL, Fleet=MSEtool::Generic_Fleet,
Obs=MSEtool::Generic_Obs, Imp=MSEtool::Perfect_Imp,
interval=4, pstar=0.5, maxF=0.8, reps=1, nsim=48, proyears=50, docheck=TRUE) {
if (is.null(Stock)) {
# message("No Stock object found. Returning a blank OM object")
# Check and add defaults
if (docheck)
.Object <- CheckOM(.Object, msg=FALSE, stop_if_missing=FALSE)
# Default MSE parameters
if (.hasSlot(.Object, "nsim")) .Object@nsim <- nsim
if (.hasSlot(.Object, "proyears")) .Object@proyears <- proyears
# interval, pstar, maxF, reps
if (.hasSlot(.Object, "interval")) .Object@interval <- interval
if (.hasSlot(.Object, "pstar")) .Object@pstar <- pstar
if (.hasSlot(.Object, "maxF")) .Object@maxF <- maxF
if (.hasSlot(.Object, "reps")) .Object@reps <- reps
return(.Object)
}
if (!methods::is(Stock, "Stock"))
print(paste("Could not build operating model:", deparse(substitute(Stock)), "not of class Stock"))
if (!methods::is(Fleet, "Fleet"))
print(paste("Could not build operating model:", deparse(substitute(Fleet)), "not of class Fleet"))
if (!methods::is(Obs, "Obs"))
print(paste("Could not build operating model:", deparse(substitute(Obs)), "not of class Obs"))
if (!methods::is(Imp, "Imp"))
print(paste("Could not build operating model:", deparse(substitute(Imp)), "not of class Imp"))
.Object@Name <- paste("Stock:", Stock@Name, " Fleet:", Fleet@Name, " Obs model:",
Obs@Name, " Imp model:", Imp@Name, sep = "")
# Now copy the values for stock, fleet and observation slots to same
# slots in the Sim object
Sslots <- slotNames(Stock)
for (i in 2:length(Sslots)) {
tt <- .hasSlot(Stock, Sslots[i]) # For back-compatibility
if (tt) slot(.Object, Sslots[i]) <- slot(Stock, Sslots[i])
}
Fslots <- slotNames(Fleet)
for (i in 2:length(Fslots)) {
tt <- .hasSlot(Fleet, Fslots[i])
if (tt) slot(.Object, Fslots[i]) <- slot(Fleet, Fslots[i])
}
Oslots <- slotNames(Obs)
for (i in 2:length(Oslots)) {
tt <- .hasSlot(Obs, Oslots[i])
if (tt) slot(.Object, Oslots[i]) <- slot(Obs, Oslots[i])
}
Islots <- slotNames(Imp)
for (i in 2:length(Islots)) {
tt <- .hasSlot(Imp, Islots[i])
if (tt) slot(.Object, Islots[i]) <- slot(Imp, Islots[i])
}
#
# source <- paste("Stock:", Stock@Source, "Fleet:", Fleet@Source, "Obs:", Obs@Source, "Imp:",Imp@Source)
slot(.Object, "Source") <- Stock@Source
# Default MSE parameters
if (.hasSlot(.Object, "nsim")) .Object@nsim <- nsim
if (.hasSlot(.Object, "proyears")) .Object@proyears <- proyears
# interval, pstar, maxF, reps
if (.hasSlot(.Object, "interval")) .Object@interval <- interval
if (.hasSlot(.Object, "pstar")) .Object@pstar <- pstar
if (.hasSlot(.Object, "maxF")) .Object@maxF <- maxF
if (.hasSlot(.Object, "reps")) .Object@reps <- reps
# Check and add defaults
.Object <- CheckOM(.Object, msg=FALSE, stop_if_missing=FALSE)
# if(length(.Object@LenCV) < 2) .Object@LenCV <- c(0.08,0.15)
# if(length(.Object@CurrentYr)==0).Object@CurrentYr=.Object@nyears
#
# # if(length(.Object@FecB) < 2) .Object@FecB <- c(3,3)
# # if(all(is.na(.Object@FecB))) .Object@FecB <- c(3,3)
# if(all(is.na(.Object@LenCV))) .Object@LenCV <- c(0.08,0.15)
# if(all(is.na(.Object@CurrentYr))) .Object@CurrentYr=.Object@nyears
#
# if(length(.Object@LR5) < 2) .Object@LR5 <- c(0,0)
# if(length(.Object@LFR) < 2) .Object@LFR <- c(0,0)
# if(length(.Object@Rmaxlen) < 2) .Object@Rmaxlen <- c(1,1)
# if(length(.Object@Fdisc) < 2) .Object@Fdisc <- c(0,0)
#
# if(all(is.na(.Object@LR5))) .Object@LR5 <- c(0,0)
# if(all(is.na(.Object@LFR))) .Object@LFR <- c(0,0)
# if(all(is.na(.Object@Rmaxlen))) .Object@Rmaxlen <- c(1,1)
# if(all(is.na(.Object@Fdisc))) .Object@Fdisc <- c(0,0)
#
# if (.hasSlot(.Object, "Size_area_1")) {
# if (length(.Object@Size_area_1)==0) .Object@Size_area_1 <- .Object@Frac_area_1
# if (all(is.na(.Object@Size_area_1))) .Object@Size_area_1 <- .Object@Frac_area_1
# }
#
# .Object@seed=1
.Object
})
# -- Hist Object Class ----
#' Class \code{'Hist'}
#'
#' An object for storing information generated by the end of the historical simulations
#'
#' @name Hist-class
#' @docType class
#'
#' @slot Data The Data object at the end of the historical period
#' @slot OMPars A numeric data.frame with nsim rows with sampled Stock, Fleet,
#' Obs, and Imp parameters.
#' @slot AtAge A named list with arrays of dimensions: `c(nsim, maxage+1, nyears+proyears)` or
#' `c(nsim, maxage+1, nyears, nareas)`
#' \itemize{
#' \item Length: Length-at-age for each simulation, age, and year
#' \item Weight: Weight-at-age for each simulation, age, and year
#' \item Select: Selectivity-at-age for each simulation, age, and year
#' \item Retention: Retention-at-age for each simulation, age, and year
#' \item Maturity: Maturity-at-age for each simulation, age, and year
#' \item N.Mortality: Natural mortality-at-age for each simulation, age, and year
#' \item Z.Mortality: Total mortality-at-age for each simulation, age, year and area
#' \item F.Mortality: Fishing mortality-at-age for each simulation, age, year and area
#' \item Fret.Mortality: Fishing mortality-at-age for retained fish for each
#' simulation, age, year and area
#' \item Number: Total numbers by simulation, age, year and area
#' \item Biomass: Total biomass by simulation, age, year and area
#' \item VBiomass: Vulnerable biomass by simulation, age, year and area
#' \item SBiomass: Spawning biomass by simulation, age, year and area
#' \item Removals: Removals (biomass) by simulation, age, year and area
#' \item Landings: Landings (biomass) by simulation, age, year and area
#' \item Discards: Discards (biomass) by simulation, age, year and area
#' }
#' @slot TSdata A named list with population and fleet dynamics:
#' \itemize{
#' \item Number: Total numbers; array dimensions `c(nsim, nyears, nareas)`
#' \item Biomass: Total biomass; array dimensions `c(nsim, nyears, nareas)`
#' \item VBiomass: Vulnerable biomass; array dimensions `c(nsim, nyears, nareas)`
#' \item SBiomass: Spawning Biomass; array dimensions `c(nsim, nyears, nareas)`
#' \item Removals: Removals (biomass); array dimensions `c(nsim, nyears, nareas)`
#' \item Landings: Landings (biomass); array dimensions `c(nsim, nyears, nareas)`
#' \item Discards: Discards (biomass); array dimensions `c(nsim, nyears, nareas)`
#' \item Find: Historical fishing mortality (scale-free); matrix dimensions `c(nsim, nyears)`
#' \item RecDev: Recruitment deviations (historical and projection); matrix dimensions `c(nsim, nyears+proyears+maxage)`
#' \item SPR: Named list with Equilibrium and Dynamic SPR (both matrices iwth dimensions `c(nsim, nyears)`)
#' \item Unfished_Equilibrium: A named list with unfished equilibrium numbers and biomass-at-age
#' }
#'
#' @slot Ref A named list with biological reference points:
#' \itemize{
#' \item ByYear: A named list with asymptotic reference points (i.e., calculated annually without recruitment deviations) all matrices with dimensions `nsim` by `nyears+proyears`:
#' \itemize{
#' \item N0: Asymptotic unfished total number
#' \item SN0: Asymptotic unfished spawning number
#' \item B0: Asymptotic unfished total biomass
#' \item SSB0: Asymptotic unfished spawning biomass
#' \item VB0: Asymptotic unfished vulnerable biomass
#' \item MSY: Asymptotic MSY
#' \item FMSY: Fishing mortality corresponding with asymptotic MSY
#' \item SSBMSY: Spawning stock biomass corresponding with asymptotic MSY
#' \item BMSY: total biomass corresponding with asymptotic MSY
#' \item VBMSY: Vulnerable biomass corresponding with asymptotic MSY
#' \item F01: Fishing mortality where the change in yield per recruit is 10% of that at F = 0
#' \item Fmax: Fishing mortality that maximizes yield per recruit
#' \item F_SPR: Fishing mortality corresponding to spawning potential ratio of 20 - 60% in increments of 5%; array dimensions \code{c(nsim, 9, nyears+proyears)}
#' \item Fcrash: Fishing mortality corresponding to the recruits-per-spawner at the origin of the stock-recruit relationship
#' \item Fmed: Fishing mortality corresponding to the median recruits-per-spawner in the historical period
#' \item SPRcrash: SPR corresponding to the recruits-per-spawner at the origin of the stock-recruit relationship
#' }
#' \item Dynamic_Unfished: A named list with dynamic unfished reference points for each simulation and year:
#' \itemize{
#' \item N0: Unfished total numbers
#' \item B0: Unfished total biomass
#' \item SN0: Unfished spawning numbers
#' \item SSB0: Unfished spawning biomass
#' \item VB0: Unfished vulnerable biomass
#' \item Rec: Unfished recruitment
#' }
#' \item ReferencePoints: A data.frame with `nsim` rows with with biological reference points
#' calculated as an average over age-of-maturity `ageM` years around the
#' current year (i.e. `nyears`):
#' \itemize{
#' \item N0: Average unfished numbers
#' \item B0: Average unfished biomass
#' \item SSB0: Average unfished spawning biomass (used to calculate depletion)
#' \item SSN0: Average unfished spawning numbers
#' \item VB0: Average unfished vulnerable biomass (used to calculate depletion if `cpar$control$D='VB'`)
#' \item MSY: Average maximum sustainable yield (equilibrium)
#' \item FMSY: Average fishing mortality corresponding with MSY
#' \item SSBMSY: Average spawning stock biomass corresponding with MSY
#' \item BMSY: Average total biomass corresponding with MSY
#' \item VBMSY: Average vulnerable biomass corresponding with MSY
#' \item UMSY: Average exploitation rate corresponding with MSY
#' \item FMSY_M: Average FMSY/M ratio
#' \item SSBMSY_SSB0: Average ratio of SSBMSY to SSB0
#' \item BMSY_B0: Average ratio of BMSY to B0
#' \item VBMSY_VB0: Average ratio of VBMSY to VB0
#' \item RefY: Maximum yield obtained in forward projections with a fixed F
#' }
#' }
#'
#' @slot SampPars A named list with all sampled Stock, Fleet, Obs, and Imp parameters
#' @slot OM The `OM` object (without cpars)
#' @slot Misc A list for additional information
#' @author A. Hordyk
#' @keywords classes
setClass("Hist", representation(
Data = 'Data',
OMPars = 'data.frame',
AtAge = 'list',
TSdata = 'list',
Ref = "list",
SampPars='list',
Misc = 'list',
OM = 'OM'
))
# ---- MSE Class ----
#' Class \code{'MSE'}
#'
#' A Management Strategy Evaluation object that contains information about
#' simulation conditions and performance of data-limited methods
#'
#'
#' @name MSE-class
#' @docType class
#' @template MSE_template
#'
#'
#' @author T. Carruthers and A. Hordyk
#' @keywords classes
setClass("MSE", representation(Name = "character",
nyears = "numeric",
proyears = "numeric",
nMPs = "numeric",
MPs = "character",
nsim = "numeric",
OM = "data.frame",
Obs = "data.frame",
SB_SBMSY='array',
F_FMSY='array',
N='array',
B = "array",
SSB="array",
VB="array",
FM = "array",
SPR='list',
Catch = "array",
Removals='array',
Effort='array',
TAC = "array",
TAE = 'array',
BioEco = 'list',
RefPoint = 'list',
CB_hist='array',
FM_hist='array',
SSB_hist='array',
Hist = 'Hist',
PPD = 'list',
Misc = 'list')
)
setMethod("initialize", "MSE", function(.Object, Name, nyears, proyears,
nMPs, MPs, nsim, OM, Obs,
SB_SBMSY, F_FMSY, N, B, SSB, VB, FM,
SPR, Catch, Removals, Effort, TAC, TAE,
BioEco, RefPoint,
CB_hist, FM_hist, SSB_hist,
Hist, PPD, Misc
) {
slts <- slotNames('MSE')
for (sl in slts) {
slot(.Object, sl) <- get(sl)
}
.Object
})
# ---- PMobj Class ----
#' An object for storing data for analysis using data-limited methods
#'
#' Used internally
#'
#' @name PMobj-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('PMobj')}
#' @slot Name Name of the Performance Metric. Character
#' @slot Caption A caption to be used in plots. Character, call, or function.
#' @slot Stat Statistic of interest for the PM. Dimensions: nsim, nMP, yrs. Array
#' @slot Ref Reference value to calculate probability for statistic. Numeric.
#' @slot Prob Probability (mean over years) Dimensions: nsim by MP. Matrix, numeric or data.frame
#' @slot Mean Mean probability (mean over years and simulations). Numeric. Length nMPs
#' @slot MPs Name of MPs. Single value. Character string
#' @author A. Hordyk
#' @keywords classes
setClass("PMobj", representation(Name = "character", Caption='label.class',
Stat='array', Ref='numeric', Prob='prob.class', Mean='numeric',
MPs="character"))
show <- function(object) methods::show(object)
#' Show the output of a PM
#'
#' @param object object of class MSE
#' @rdname show-MSE
#' @export
setMethod("show", signature = (object="PMobj"), function(object) {
sls <- slotNames(object)
df <- data.frame(Slot=sls, Value=NA, stringsAsFactors = FALSE)
if (any(sapply(sls, function(sl) length(slot(object, sl))) == 0)) {
cat("Incomplete PMobj\n")
for (sl in sls) {
r <- match(sl, sls)
slval <- slot(object, sl)
if ('array' %in% class(slval) & length(slval)>0 & length(dim(slval))>2) {
df[r,2] <- 'array'
} else if ('matrix' %in% class(slval) & length(slval)>0) {
df[r,2] <- 'matrix'
} else if (length(slval) > 0 & ! 'call' %in% class(slval)) {
df[r,2] <- slval
} else {
df[r, 2] <- 'not defined'
}
}
print(df)
} else {
cat(object@Name)
cat("\n", object@Caption)
cat("\n")
nMP <- length(object@MPs)
if (nMP > 1) nsim <- dim(object@Prob)[1]
if (nMP == 1) nsim <- length(object@Prob)
nprint <- min(nsim, 10)
if (nMP > 1) df <- data.frame(object@Prob[1:nprint,])
if (nMP == 1) df <- data.frame(object@Prob[1:nprint])
if (nsim > (nprint+1)) {
if (nMP > 1) lst <- object@Prob[nprint+1,]
if (nMP == 1) lst <- object@Prob[nprint+1]
} else {
if (nMP > 1) lst <- object@Prob[nprint,]
if (nMP == 1) lst <- object@Prob[nprint]
}
df <- round(df,2)
lst <- round(lst,2)
colnames(df) <- object@MPs
names(lst) <- object@MPs
if (nsim > (nprint+1)) {
df <- rbind(df,
rep(".", nMP),
rep(".", nMP),
rep(".", nMP),
lst)
rownames(df) <- c(1:(nprint+3), nsim)
}
print(df)
cat("\nMean\n")
print(round(object@Mean,2))
}
})
# ---- Summary of MSE object ----
#' Summary of MSE object
#'
#' @param object object of class MSE
#' @param ... a list of names of PM methods
#' @param silent Should summary be printed to console? Logical.
#' @param Refs An optional named list (matching the PM names) with numeric values to override the default `Ref` values. See examples.
#' @rdname summary-MSE
#' @export
setMethod('summary', signature="MSE", function(object, ..., silent=FALSE, Refs=NULL) {
PMlist <- unlist(list(...))
if(length(PMlist) == 0) PMlist <- c("PNOF", "P50", "AAVY", "LTY")
if (!methods::is(PMlist,'character')) stop("Must provide names of PM methods")
# check
for (X in seq_along(PMlist))
if (!methods::is(get(PMlist[X]), "PM")) stop(PMlist[X], " is not a valid PM method")
if (!silent) message("Calculating Performance Metrics")
storeMean <- vector('list', length(PMlist))
storeName <- vector('list', length(PMlist))
storeCap <- vector('list', length(PMlist))
storeHeading <- vector('list', length(PMlist))
storeMP <- vector('list', length(PMlist))
for (X in 1:length(PMlist)) {
ref <- Refs[[PMlist[X]]]
if (is.null(ref)) {
runPM <- eval(call(PMlist[X], object))
} else {
runPM <- eval(call(PMlist[X], object, Ref=ref))
}
storeMean[[X]] <- runPM@Mean
storeName[[X]] <- runPM@Name
storeCap[[X]] <- runPM@Caption
storeMP[[X]] <- runPM@MPs
}
df <- data.frame('MP'=storeMP[[1]], signif(do.call('cbind', storeMean),2), stringsAsFactors = FALSE)
# heading <- do.call('rbind', storeHeading)
colnames(df)[2:(length(PMlist)+1)] <- PMlist #caps # gsub(" ", "", caps)
if (!silent) {
dfprint <- data.frame('Performance Metrics' = do.call('rbind', storeName), gap="", do.call('rbind', storeCap))
names(dfprint)[2:3] <- ''
print(dfprint)
cat("\n")
cat("\nPerformance Statistics:\n")
print(df)
}
invisible(df)
})
# # ---- Plot Data Object -----
# #' Plot Data object
# #'
# #' @rdname plot-Data
# #' @param x object of class Data
# #' @param upq Upper quantile of TACs for max ylim
# #' @param lwq Lower quantile of TACs for min ylim
# #' @param outline Logical. Include outliers in plot?
# #' @param ... Optional additional arguments passed to \code{boxplot}
# #' @export
# setMethod("plot",
# signature(x = "Data"),
# function(x, upq=0.9, lwq=0.1, outline = FALSE, ...){
#
# old_par <- par(no.readonly = TRUE)
# on.exit(par(list = old_par), add = TRUE)
# boxplot.Data(x, upq, lwq, outline, ...)
# })
# # Data<-x
# if (class(Data) != "Data") stop("Must supply object of class Data")
# if (all(is.na(Data@TAC))) stop("No TAC data found")
# cols<-rep(c('black','red','green','blue','orange','brown','purple','dark grey','violet','dark red','pink','dark blue','grey'),4)
# ltys<-rep(1:4,each=13)
#
# if(is.na(funcs[1]))funcs<-Data@MPs
#
# nMPs<-length(funcs)
# nplots<-ceiling(nMPs/maxlines)
# maxl<-ceiling(nMPs/nplots)
# mbyp <- split(1:nMPs, ceiling(1:nMPs/maxl)) # assign methods to plots
#
# if(is.na(xlims[1])|length(xlims)!=2){
# xlims<-quantile(Data@TAC,c(0.005,0.95),na.rm=T)
# if(xlims[1]<0)xlims[1]<-0
# }
# if(!NAor0(Data@Ref)){
# if(xlims[1]>Data@Ref)xlims[1]<-max(0,0.98*Data@Ref)
# if(xlims[2]<Data@Ref)xlims[2]<-1.02*Data@Ref
# }
# ylims<-c(0,1)
#
# #for(m in 1:nMPs){
# # if(sum(!is.na(Data@TAC[m,,1]))>2){
# # dens<-density(Data@TAC[m,,1],na.rm=T)
# #print(quantile(dens$y,0.99,na.rm=T))
# # if(quantile(dens$y,0.9,na.rm=T)>ylims[2])ylims[2]<-quantile(dens$y,0.90,na.rm=T)
# #}
# #}
#
# #dev.new2(width=10,height=0.5+7*nplots)
# par(mfrow=c(ceiling(nplots/2),2),mai=c(0.4,0.4,0.01,0.01),omi=c(0.35,0.35,0.35,0.05))
#
# for(p in 1:nplots){
# m<-mbyp[[p]][1]
# plot(NA,NA,xlim=xlims,ylim=ylims,main="",xlab="",ylab="",col="white",lwd=3,type="l")
# abline(h=0)
# if(!NAor0(Data@Ref)){
# abline(v=Data@Ref,col="light grey",lwd=2)
# if(!NAor0(Data@Ref_type[1]))legend('right',Data@Ref_type,text.col="grey",bty='n')
# }
# #plot(density(DLM@TAC[m,,1],from=0,na.rm=T),xlim=xlims,ylim=ylims,main="",xlab="",ylab="",col=coly[m],lty=ltyy[m],type="l")
#
# if(!is.na(perc[1]))abline(v=quantile(Data@TAC[m,,1],p=perc,na.rm=T),col=cols[m],lty=ltys[m])
# #if(length(mbyp[[p]])>0){
# for(ll in 1:length(mbyp[[p]])){
# m<-mbyp[[p]][ll]
# if(sum(!is.na(Data@TAC[m,,1]))>10){ # only plot if there are sufficient non-NA TAC samples
# x<-density(Data@TAC[m,,1],from=0,na.rm=T)$x
# y<-density(Data@TAC[m,,1],from=0,na.rm=T)$y
# y<-y/max(y)
# lines(x,y,col=cols[ll])
# }else{
# print(paste("Method ",funcs[m]," produced too many NA TAC values for plotting densities",sep=""))
# }
# if(!is.na(perc[1]))abline(v=quantile(Data@TAC[m,,1],p=perc,na.rm=T),col=cols[ll],lty=2)
# }
# #}
# cind<-1:length(mbyp[[p]])
# legend('topright',funcs[mbyp[[p]]],text.col=cols[cind],col=cols[cind],lty=1,bty='n',cex=0.75)
# }
#
# mtext(paste("TAC (",Data@Units,")",sep=""),1,outer=T,line=0.5)
# mtext(paste("Standardized relative frequency",sep=""),2,outer=T,line=0.5)
# mtext(paste("TAC calculation for ",Data@Name,sep=""),3,outer=T,line=0.5)
# })
# # ---- Plot MSE object ----
# #' Plot MSE object
# #'
# #' @rdname plot-MSE
# #' @param x object of class MSE
# #' @export
# setMethod("plot",
# signature(x = "MSE"),
# function(x){
# MSEobj<-x
# Pplot(MSEobj)
# Kplot(MSEobj)
# Tplot(MSEobj)
# })
# ---- Summary of Data Object ----
#' Summary of Data object
#'
#' @rdname summary-Data
#' @param object An object of class Data
#' @param wait Logical. Wait for key press before next plot?
#' @param x iteration number for the Data object.
#' @param plots Character. What plots to show? `all`, `TS`, `CAA`, `CAL`, `PD`
#' for all plots, time-series, catch-at-age, catch-at-length, and
#' probability distributions respectively
#' @param rmd Logical. Used in a rmd file?
#' @param head Character. Heading for rmd file. Default is '##' (second level heading)
#' @param tplot Integer. Number of plots per page. Default 25
#' @export
setMethod("summary",
signature(object = "Data"),
function(object, wait=TRUE, x=1, plots='all', rmd=FALSE, head="##", tplot=25){
plots <- match.arg(plots, c('all', 'TS', 'CAA', 'CAL', 'PD'), several.ok = TRUE)
if ('all' %in% plots) plots <- c('TS', 'CAA', 'CAL', 'PD')
Freq <- n <- Var2 <- NULL # cran check
if (!methods::is(object, "Data")) stop("Object must be class `Data`", call.=FALSE)
# Time-Series
Year <- object@Year
l <- length(Year)
slts <- c("Cat", 'Ind', 'SpInd', 'VInd', 'Rec', 'ML', 'Lc')
Cat <- Ind <- SpInd <- VInd <- Rec <- ML <- Lc <- NULL # cran checks
for (sl in slts) {
tt <- slot(object, sl)[x,]
if (length(tt)!=l) tt <- c(tt, rep(NA,l-length(tt)))
assign(sl, tt)
}
Val <- c(Cat, Ind, SpInd, VInd, Rec, ML, Lc)
Var <- rep(c("Catch", "Total Index", "Spawning Index",
"Vuln. Index", "Recruitment", "Mean Length", "Mean Length above Lc"), each=length(Year))
ts.df <- data.frame(Year=object@Year, Val=Val, Var=Var, stringsAsFactors = TRUE)
# ts.df$Year <- as.factor(ts.df$Year)
ts.df$Var <- factor(ts.df$Var, levels=
c("Catch", "Total Index", "Spawning Index",
"Vuln. Index", "Recruitment", "Mean Length", "Mean Length above Lc"))
ts.df <- subset(ts.df, !is.na(Val))
if (nrow(ts.df)>0 && 'TS' %in% plots) {
P1 <- ggplot2::ggplot(ts.df, ggplot2::aes(x=Year, y=Val, group = Var)) +
ggplot2::facet_wrap(~Var, scales='free_y') + ggplot2::geom_line(linewidth=1.25) +
ggplot2::theme_classic() +
ggplot2::expand_limits(y=0) +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 1)) +
ggplot2::scale_x_continuous(breaks=pretty(rev(Year), length(Year)/5)) +
ggplot2::labs(y="")
} else {
P1 <- NULL
}
if (!is.null(P1)) {
if (rmd) {
cat('\n')
cat('\n')
cat(head, 'Time-Series')
cat('\n')
} else {
message('Plotting Time-Series')
}
print(P1)
}
if (interactive() & wait & !is.null(P1))
invisible(readline(prompt="Press [enter] to continue..."))
# Additional indices
P1a <- NULL
if(sum(dim(object@AddInd))>3) {
# additional indices exist
dd <- dim(object@AddInd[1,,])
nind <- dd[1]
nyears <- dd[2]
if (!is.null(dimnames(object@AddInd))) {
ind_names <- as.character(dimnames(object@AddInd)[[2]])
} else {
ind_names <- 1:nind
}
add.ind <- data.frame(Year=rep(Year, each=nind),
Index=ind_names,
Value=as.vector(object@AddInd[1,,]),
stringsAsFactors = TRUE)
if(nrow(add.ind)>0 && 'TS' %in% plots) {
P1a <- ggplot2::ggplot(add.ind, ggplot2::aes(x=Year, y=Value)) +
ggplot2::facet_wrap(~Index, scales='free_y') + ggplot2::geom_line(linewidth=1.25) +
ggplot2::theme_classic() +
ggplot2::expand_limits(y=0) +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, hjust = 1)) +
ggplot2::scale_x_continuous(breaks=pretty(rev(Year), length(Year)/5)) +
ggplot2::labs(y="Index")
} else {
P1a <- NULL
}
if (!is.null(P1a)) {
if (rmd) {
cat('\n')
cat('\n')
cat(head, 'Additional Indices')
cat('\n')
} else {
message('Plotting Additional Indices')
}
print(P1a)
}
}
if (interactive() & wait & !is.null(P1a))
invisible(readline(prompt="Press [enter] to continue..."))
# CAA
if (all(is.na(object@CAA))) {
P2 <- FALSE
} else {
P2 <- TRUE
}
if (P2 == TRUE) {
CAA <- object@CAA[x,,]
nyrs <- nrow(CAA); maxage <- ncol(CAA)
dimnames(CAA) <- list(1:nyrs, 1:maxage)
df1 <- as.data.frame.table(CAA, stringsAsFactors = FALSE)
colnames(df1) <- c("Year", "Val", "Freq")
df1$Val <- as.numeric(df1$Val)
df1$Year <- as.numeric(df1$Year)
yr.n <- df1 %>% dplyr::group_by(Year) %>% dplyr::summarise(n=sum(Freq))
yr.ind <- yr.n %>% dplyr::filter(n>0) %>% dplyr::select(Year)
Years <- object@Year
nyears <- length(unique(df1$Year))
df1$Year_val <- (Years[(length(Years)-nyears+1):length(Years)])
df1 <- df1[!is.na(df1$Freq),] # Fliter out NA values, so we don't try to plot missing years
if (nrow(df1)>0 && 'CAA' %in% plots) {
if (rmd) {
cat('\n')
cat('\n')
cat(paste(head, 'Catch-at-Age'))
cat('\n')
} else {
message('Plotting Catch-at-Age')
}
nyears <- length(unique(df1$Year))
nbins <- length(unique(df1$Val))
if (nyears > tplot) {
npages <- ceiling(nyears/tplot)
ncol <- 5
nrow <- 5
nplot <- ncol * nrow
} else {
npages <- 1
nrow <- ceiling(sqrt(nyears))
ncol <- ceiling(nyears/nrow)
nplot <- nyears
}
pmat <- matrix(1:(nrow*ncol), nrow=nrow, ncol=ncol, byrow=TRUE)
pmat[pmat >nplot] <- NA
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
yr1 <- 1
col <- "grey"
for (pg in 1:npages) {
if(npages>1)message('Plot ', pg, ' of ', npages)
yrind <- unique(df1$Year)[yr1:(yr1+nplot-1)]
yr1 <- max(yrind) + 1
dat <- df1 %>% dplyr::filter(Year %in% yrind)
un.yrs_val <- as.numeric(unique(dat$Year_val))
un.yrs <- as.numeric(unique(dat$Year))
if (pg >1 && pg == npages) {
nplot <- nyears - (npages-1) * tplot
ncol <- ceiling(sqrt(nplot))
nrow <- ceiling(nplot/ncol)
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
pmat <- matrix(1:nplot, nrow=nrow, ncol=ncol, byrow=TRUE)
}
for (p in 1:nplot) {
pdat <- dat %>% dplyr::filter(Year==un.yrs[p])
if (nrow(pdat) > 0) {
if (p %in% pmat[nrow,]) {
barplot(pdat$Freq, names=round(pdat$Val, 2), axes=FALSE, col=col)
} else {
barplot(pdat$Freq, names=FALSE, axes=FALSE, col=col)
}
if (p %in% pmat[,1]) axis(side=2)
if (!p %in% pmat[,1]) axis(side=2, labels=TRUE)
ncount <- round(sum(pdat$Freq),0)
title(un.yrs_val[p])
text(max(pdat$Val), max(pdat$Freq), paste('n = ', ncount),
xpd=NA)
}
}
mtext(side=1, outer=TRUE, "Age", line=2, cex=1.5)
mtext(side=2, outer=TRUE, "Frequency", line=2, cex=1.5)
}
}
}
if (interactive() & wait & !is.null(P2) && P2)
invisible(readline(prompt="Press [enter] to continue..."))
# CAL
if (all(is.na(object@CAL))) {
P3 <- FALSE
} else {
P3 <- TRUE
}
if (P3 == TRUE) {
CAL <- object@CAL[x,,]
nyrs <- nrow(CAL); nbins <- length(object@CAL_bins) - 1
By <- object@CAL_bins[2] - object@CAL_bins[1]
BinsMid <- seq(object@CAL_bins[1] + 0.5*By, by=By,length.out = nbins)
dimnames(CAL) <- list(1:nyrs, BinsMid)
df1 <- as.data.frame.table(CAL, stringsAsFactors = FALSE)
colnames(df1) <- c("Year", "Val", "Freq")
df1$Val <- as.numeric(df1$Val)
df1$Year <- as.numeric(df1$Year)
yr.n <- df1 %>% dplyr::group_by(Year) %>% dplyr::summarise(n=sum(Freq, na.rm=T))
yr.ind <- yr.n %>% dplyr::filter(n>0) %>% dplyr::select(Year)
Years <- object@Year
nyears <- length(unique(df1$Year))
df1$Year_val <- (Years[(length(Years)-nyears+1):length(Years)])
if (nrow(df1) > 0 && 'CAL' %in% plots) {
if (rmd) {
cat('\n')
cat('\n')
cat(paste(head, 'Catch-at-Length'))
cat('\n')
} else {
message('Plotting Catch-at-Length')
}
nyears <- length(unique(df1$Year))
nayears <- df1 %>% group_by(Year) %>% summarize(isna=all(is.na(Freq)))
nyears <- sum(!nayears$isna)
nbins <- length(unique(df1$Val))
if (nyears > tplot) {
npages <- ceiling(nyears/tplot)
ncol <- 5
nrow <- 5
nplot <- ncol * nrow
} else {
npages <- 1
nrow <- ceiling(sqrt(nyears))
ncol <- ceiling(nyears/nrow)
nplot <- nyears
}
pmat <- matrix(1:(ncol*nrow), nrow=nrow, ncol=ncol, byrow=TRUE)
pmat[pmat>nplot] <- NA
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
yr1 <- 1
col <- "grey"
for (pg in 1:npages) {
if(npages>1)message('Plot ', pg, ' of ', npages)
if (pg ==1) {
yrind <- yr.ind$Year[1:(nplot)]
} else {
t1 <- nplot * (pg-1) +1
t2 <- t1+nplot
yrind <- yr.ind$Year[t1:t2]
}
yr1 <- max(yrind) + 1
dat <- df1 %>% dplyr::filter(Year %in% yrind)
un.yrs_val <- as.numeric(unique(dat$Year_val))
un.yrs <- as.numeric(unique(dat$Year))
if (pg >1 && pg == npages) {
nplot <- nyears - (npages-1) * tplot
ncol <- ceiling(sqrt(nplot))
nrow <- ceiling(nplot/ncol)
op <- par(mfrow=c(nrow, ncol), no.readonly = TRUE, mar=c(2,2,2,1), oma=c(4,4,2,0))
on.exit(par(op))
pmat <- matrix(1:(ncol*nrow), nrow=nrow, ncol=ncol, byrow=TRUE)
pmat[pmat>nplot] <- NA
}
for (p in 1:nplot) {
pdat <- dat %>% dplyr::filter(Year==un.yrs[p])
if (nrow(pdat) > 0) {
if (all(is.na(pdat$Freq))) {
} else{
if (p %in% pmat[nrow,]) {
barplot(pdat$Freq, names.arg=round(pdat$Val, 2), axes=FALSE, col=col, las=2)
} else {
barplot(pdat$Freq, names.arg=FALSE, axes=FALSE, col=col)
}
if (p %in% pmat[,1]) axis(side=2)
if (!p %in% pmat[,1]) axis(side=2, labels=TRUE)
ncount <- round(sum(pdat$Freq),0)
title(un.yrs_val[p])
text(length(unique(df1$Val)), max(pdat$Freq), paste('n = ', ncount), xpd=NA)
}
}
}
mtext(side=1, outer=TRUE, "Length", line=2, cex=1.5)
mtext(side=2, outer=TRUE, "Frequency", line=2, cex=1.5)
}
}
}
if (interactive() & wait & !is.null(P3) && P3)
invisible(readline(prompt="Press [enter] to continue..."))
# Biology & Depletion
slots<-c("Dep","Mort","FMSY_M","Dt","BMSY_B0","vbK", "vbLinf")
namey<-c("Stock depletion", "Natural Mortality rate","Ratio of FMSY to M",
"Depletion over time t","BMSY relative to unfished",
"von B. k parameter", "von B. Linf parameter")
slotsCV<-c("CV_Dep","CV_Mort","CV_FMSY_M","CV_Dt","CV_BMSY_B0","CV_vbK", "CV_vbLinf")
reps <- 5000
val <- list()
for (i in seq_along(slots)) {
mu <- slot(object, slots[i])[x]
cv <- slot(object, slotsCV[i])[x]
val[[i]] <- trlnorm(reps, mu,cv)
}
vals <- do.call("cbind", val)
colnames(vals) <- namey
df1 <- as.data.frame.table(vals, stringsAsFactors = TRUE)
df1 <- df1 %>% dplyr::filter(is.na(Freq) == FALSE)
if (nrow(df1) > 0 && 'PD' %in% plots) {
P4 <- ggplot2::ggplot(df1, ggplot2::aes(x=Freq, group=Var2)) +
ggplot2::facet_wrap(~Var2, scales="free") + ggplot2::geom_histogram(bins=30) +
ggplot2::labs(y="Frequency", x="Parameter Value") +
ggplot2::theme_classic()
} else {
P4 <- NULL
}
if (!is.null(P4)) {
if (rmd) {
cat('\n')
cat('\n')
cat(paste(head, 'Parameter Distributions'))
cat('\n')
} else {
message('Plotting Parameter Distributions')
}
print(P4)
}
})
# -- Management Recommendation Class ----
#' Class \code{'Rec'}
#'
#' An object for storing the MP recommendations
#'
#' @name Rec-class
#' @docType class
#' @section Objects from the Class: Objects can be created by calls of the form
#' \code{new('Rec')}
#' @slot TAC A numeric value with the TAC recommendation
#' @slot Effort A numeric value with the effort recommendation as a fraction of current (nyear) fishing effort
#' @slot Spatial A boolean vector of length 'nareas' specifying if area is open (1) or closed (0) to fishing
#' @slot Allocate A boolean value describing if effort should be re-allocated from close to open areas
#' @slot LR5 smallest length at 5 per cent retention - in absolute units - i.e same units as Linf and L50
#' @slot LFR smallest length at full retention - in absolute units - i.e same units as Linf and L50
#' @slot HS upper harvest slot (no retention above this) - in absolute units - i.e same units as Linf and L50
#' @slot Rmaxlen retention of the largest size class - fraction between 0 and 1
#' @slot L5 smallest length at 5 per cent selection - in absolute units - i.e same units as Linf and L50
#' @slot LFS smallest length at full selection - in absolute units - i.e same units as Linf and L50
#' @slot Vmaxlen selection of the largest size class - fraction between 0 and 1
#' @slot Fdisc fraction of discarded fish that die - fraction between 0 and 1
#' @slot DR Discard rate - the fraction of caught fish that are discarded
#' @slot Misc An empty list that can be used to store information and pass on to MPs in future
#' @author A. Hordyk
#' @keywords classes
setClass("Rec", representation(
TAC = "numeric",
Effort = "numeric",
Spatial="numeric", Allocate = "numeric",
LR5 = "numeric", LFR = "numeric", HS="numeric", Rmaxlen="numeric",
L5 = "numeric", LFS = "numeric", Vmaxlen="numeric",
Fdisc = "numeric",
DR='numeric',
Misc="list"))
#' Show the output of a single MP recommendation
#'
#' @param object object of class Rec
#' @rdname show-Rec
#' @export
setMethod("show", signature = (object="Rec"), function(object) {
Rec <- object
slots <- slotNames(Rec)
recList <- list()
perc <- 0.5
for (X in slots) { # sequence along recommendation slots
if (X == "Misc") { # convert to a list nsim by nareas
rec <- slot(Rec, X)
} else {
rec <- slot(Rec,X) # unlist(lapply(temp, slot, name=X))
}
if (X == "Spatial") { # convert to a matrix nsim by nareas
nsims <- 1
nareas <- max(2,length(rec))
rec <- matrix(rec, nareas, nsims, byrow=FALSE)
}
recList[[X]] <- rec
}
names <- c("TAC", "Effort", "LR5", "LFR", "HS", "Rmaxlen",
"L5", "LFS", 'Vmaxlen', 'Fdisc', 'DR', 'Spatial')
mat <- matrix(0, nrow=1, ncol=length(names)+nareas-1)
count <- 0
for (x in names) {
temp <- recList[[x]]
count <- count + 1
if (x!="Spatial") {
mat[,count] <- quantile(temp, probs=perc, na.rm=TRUE)
} else {
mat[,count:ncol(mat)] <- t(matrix(unlist(temp), nrow=nareas, ncol=1))
}
}
names[length(names)] <- "Area 1"
names <- c(names, paste('Area', 2:nareas))
if (perc !=0.5) names[1] <- paste0(names[1], ' (', perc, ' perc')
if (perc ==0.5) names[1] <- paste0(names[1], ' (median)')
colnames(mat) <- names
if (nrow(mat) == 1) {
mat <- as.data.frame(mat)
matout <- mat[!is.na(mat)]
names(matout) <- names[!is.na(mat)]
print(matout)
}
})
# ---- Internal show methods ----
#' @importFrom utils capture.output str
show_int <- function(object, slots_check) {
cat(paste0("S4 object of class ", dQuote(class(object)), ".\n"))
if (!missing(slots_check)) {
for(i in slots_check) {
len <- length(slot(object, i))
if (len) cat(paste0(len, " items found in slot ", dQuote(i), ".\n"))
}
}
cat("\n")
cat("Use str(), slotNames(),", dQuote("@"), "etc. to explore contents:\n\n")
txt <- capture.output(utils::str(object))
for(i in txt[1:5]) cat(i, "\n")
invisible()
}
#' @name show-MSEtool
#' @aliases show show,Data-method show,OM-method show,Hist-method show,MSE-method show,MMSE-method
#' @title Show MSEtool S4 objects
#'
#' @description Briefly prints a couple of lines from \link[utils]{str} to avoid swamping the console with
#' the contents of very large objects.
#' @param object S4 object from MSEtool
setMethod("show", "Data", function(object) show_int(object))
#' @rdname show-MSEtool
setMethod("show", "OM", function(object) show_int(object, slots_check = "cpars"))
#' @rdname show-MSEtool
setMethod("show", "Hist", function(object) show_int(object))
#' @rdname show-MSEtool
setMethod("show", "MSE", function(object) show_int(object))
Try the MSEtool package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.