ssmvln

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ss3om
 Title: Tools for Conditioning Fisheries Operating Models Using Stock Synthesis 3
-Version: 0.5.2.9008
+Version: 0.5.2.9010
 Authors@R: person("Iago", "Mosqueira", email = "[email protected]",
   role = c("aut", "cre"))
 Description: Tools for loading Stock Synthesis (SS3) models into FLR. Used in
@@ -18,7 +18,8 @@ Imports:
   data.table,
   foreach,
   FLFishery,
-  mse
+  mse,
+  mvtnorm
 Additional_repositories: http://flr-project.org/R
 Remotes: r4ss/r4ss
 BugReports: https://github.com/flr/ss3om/issues

NAMESPACE

@@ -6,6 +6,7 @@ import("data.table")
 import("foreach")
 import("r4ss")
 importFrom("stats", "optimise", "setNames")
+importFrom("mvtnorm", "rmvnorm")
 
 export(
   "buildFLBFss3",
@@ -42,6 +43,7 @@ export(
   "prepareRetro",
   "readCDss3",
   "readFLIBss3",
+  "readFLoemss3",
   "readFLomss3",
   "readFLRPss3",
   "readFLSss3",
@@ -54,7 +56,8 @@ export(
   "readRESIDss3",
   "readOMSss3",
   "runss3grid",
-  "runtime"
+  "runtime",
+  "ssmvln"
 )
 
 exportMethods(

R/read.R

@@ -43,23 +43,23 @@ readFLIBss3 <- function(dir, fleets, birthseas=out$birthseas,
 } # }}}
 
 # readFLomss3 {{{
-readFLomss3 <- function(dir, birthseas=out$birthseas,
-  repfile="Report.sso", compfile="CompReport.sso", ...) {
+readFLomss3 <- function(dir, repfile="Report.sso",
+  compfile="CompReport.sso", ...) {
 
   # LOAD SS_output list
   out <- readOutputss3(dir, repfile=repfile, compfile=compfile)
 
   # FLS
-  if(out$SS_versionNumeric == 3.24)
-    stk <- buildFLSss3(out, birthseas=birthseas, ...)
-  else
-    stk <- buildFLSss330(out, ...)
+  stk <- readFLSss3(dir, ...)
 
   # FLSR
   srr <- buildFLSRss3(out)
 
   # RPs
-  rps <- buildFLRPss3(out)
+  if(out$SS_versionNumeric == 3.24)
+    rps <- buildFLRPss3(out)
+  else
+    rps <- buildFLRPss330(out)
 
   return(FLom(stock=stk, sr=srr, refpts=rps))
 } # }}}

R/ssmvln.R

@@ -0,0 +1,139 @@
+# ssmvln.R - DESC
+# ss3om/R/ssmvln.R
+
+# Copyright (c) L.T. Kell, 2023.
+# Author: Laurie KELL <[email protected]>
+#
+# Distributed under the terms of the EUPL-1.2
+
+# internal {{{
+
+# sshat
+
+sshat <- function(hat) {
+
+  # SET lowercase names
+  names(hat) <- tolower(names(hat))
+
+  # EXTRACT Recr, Bratio and F
+  y <- rbind(
+    hat[grep(paste("Recr", "",  sep="_"), hat$label),],
+    hat[grep(paste("Bratio", "", sep="_"), hat$label),],
+    hat[grep(paste("F", "", sep="_"), hat$label),])
+
+  y <- y[substr(y$label, 1, 2) %in% c("F_", "Br", "Re"),]
+
+  # DROP not-annual values
+  y <- subset(y, !label %in% c("Recr_Initial","Recr_unfished",
+    "Recr_Virgin","Ret_Catch_MSY"))
+
+  # GET CV and Var
+  y <- transform(y, cv2=(y$stddev / y$value) ^ 2,
+    var=log(1 + (y$stddev / y$value) ^ 2))[, -c(4,5)]
+
+  # ASSIGN names
+  names(y)=c("label","hat","stdLog","std","cv")
+
+  return(y)
+}
+
+# sscor, returns covariances as matrix
+
+sscor <- function(covar){
+
+  # UNIFY colnames
+  setnames(covar, tolower(names(covar)))
+
+  # EXTRACT
+  flag <- unique(sort(c(
+    grep(paste("Recr",  "",sep="_"),covar$label.i),
+    grep(paste("Recr",  "",sep="_"),covar$label.j),
+    grep(paste("Bratio","",sep="_"),covar$label.i),
+    grep(paste("Bratio","",sep="_"),covar$label.j),
+    grep(paste("F",     "",sep="_"),covar$label.i),
+    grep(paste("F",     "",sep="_"),covar$label.j))))
+
+  # GET correlations for all Fs, Bratios & Recrs
+  cor <- covar[flag, c("label.i", "label.j", "corr")]
+
+  flag <- substr(cor$label.i, 1, 2) %in% c("Re","F_", "Br") & 
+    substr(cor$label.j, 1, 2) %in% c("Re","F_", "Br")
+  cor <- cor[flag,]
+
+  cor <- rbind(cor, transform(cor, label.i=label.j, label.j=label.i))
+  cor <- subset(cor, !(is.na(label.i) | is.na(label.j)))
+
+  cor <- dcast(cor, label.i ~ label.j, value.var="corr")
+  dmns <- unname(unlist(cor[,1]))
+
+  cor <- as.matrix(cor[, -1])
+  cor[is.na(cor)] <- 0
+  diag(cor) <- 1
+
+  dimnames(cor) <- list(dmns,dmns)
+
+  return(cor)
+}
+
+# }}}
+
+# ssmvln {{{
+
+#' @examples
+#' data(sma)
+#' head(ssmvln(sma$CoVar, sma$derived_quants, as.FLQuants=FALSE))
+#' ssmvln(sma)
+#' plot(ssmvln(sma))
+
+ssmvln <- function(covar, hat=NULL, mc=500, new=!FALSE, as.FLQuants=TRUE) {
+
+  # EXTRACT if covar is SS_output list
+
+  if(is.null(hat) & is.list(covar)) {
+    hat <- covar$derived_quants
+    covar <- covar$CoVar
+  }
+
+  names(hat) <- tolower(names(hat))
+  names(covar) <- tolower(names(covar))
+
+  # cov(x,y)
+  cor <- sscor(data.table(covar))
+
+  hat <- data.table(hat)
+  hat <- subset(hat,label%in%dimnames(cor)[[1]])
+  hat <- sshat(hat)
+  cor <- cor[hat$label, hat$label]
+
+  ## calculate Covariance matrix
+  if(new)
+    cvr <- cor2cov(cor, hat$cv)
+  else
+    cvr <- cor2cov(cor, hat$std^2)
+
+  cvr[!is.finite(cvr)] <- 0
+
+  mvlmu <- log(hat$hat) 
+  names(mvlmu) <- dimnames(cor)[[1]]
+
+  ## MC it all
+  rtn <- exp(rmvnorm(mc, mean = mvlmu, sigma=cvr, method = c("svd")))
+
+  ## return object
+  rtn <- data.table(rtn)
+  names(rtn) <- dimnames(cor)[[1]]
+
+  rtn[, iter:=seq(mc)]
+  rtn <- melt(rtn, id="iter")
+  rtn <- rtn[, c("variable", "year"):= tstrsplit(variable, "_")]
+
+  setnames(rtn, c("variable", "value"), c("qname", "data"))
+
+  #
+  if(as.FLQuants)
+    rtn <- as(rtn, 'FLQuants')
+
+  return(rtn)
+}
+# }}}
+

R/utilities.R

@@ -212,99 +212,3 @@ prepareRetro <- function(path, starter="starter.ss", years=5) {
 
   invisible(TRUE)
 } # }}}
-
-# sscor
-
-#' @examples
-#' data(sma)
-#' sscor(sma$CoVar)
-
-sscor<-function(covar){
-
-  # UNIFY colnames
-  setnames(covar, tolower(names(covar)))
-
-  #
-  flag <- unique(sort(c(
-    grep(paste("Bratio","",sep="_"),covar$label.i),
-    grep(paste("Bratio","",sep="_"),covar$label.j),
-    grep(paste("F",     "",sep="_"),covar$label.i),
-    grep(paste("F",     "",sep="_"),covar$label.j))))
-
-  # GET correlations for all Fs & Bratios
-  cor <- covar[flag, c("label.i", "label.j", "corr")]
-
-  flag <- substr(cor$label.i, 1, 2) %in% c("F_", "Br") & 
-    substr(cor$label.j, 1, 2) %in% c("F_", "Br")
-  cor <- cor[flag,]
-
-  cor =rbind(cor,transform(cor,label.i=label.j,label.j=label.i))
-  cor =subset(cor,!(is.na(label.i)|is.na(label.j)))
-
-  cor =dcast(cor, label.i ~ label.j, value.var="corr")
-  dmns=unname(unlist(cor[,1]))
-#  dimnames(cor)[[1]] <- dmns
-#  dmns=dimnames(cor)[[1]]
-
-  cor=as.matrix(cor[,-1])
-  cor[is.na(cor)]=0
-  diag(cor)=1
-
-  dimnames(cor)=list(dmns,dmns)
-  cor
-}
-
-#' @examples
-#' data(sma)
-#' ssmvln(sma$CoVar, sma$derived_quants)
-
-ssmvln <- function(covar,hat,mc=5000,new=!FALSE){
-
-  names(hat)  =tolower(names(hat))
-  names(covar)=tolower(names(covar))
-
-  ## Correlations on the log scale
-  ## Estimates on the untransformed scale
-
-  #cov(x,y)
-  cor <- sscor(data.table(covar))
-  hat <- data.table(hat)
-
-  if (mc<=1) return(cor)
-
-  ## Get estimates of F and Bratio
-  y=rbind(hat[grep(paste("Bratio","",sep="_"),hat$label),],
-          hat[grep(paste("F",      "",sep="_"),hat$label),])
-  y=y[substr(y$label,1,2)%in%c("F_","Br"),]
-  y=transform(y,cv2=(y$stddev/y$value)^2)
-
-  ## calculate Covariance matrix
-  if(new)
-    cvr=cor2cov(cor,y$cv2)
-  else
-    cvr=cor2cov(cor,y$var)
-
-  cvr[!is.finite(cvr)]=0
-
-  mvlmu=log(y$value) 
-  names(mvlmu)=dimnames(cor)[[1]]
-
-  ## MC it all
-  rtn=exp(mvtnorm::rmvnorm(mc, mean = mvlmu, sigma=cvr, method = c("svd")))
-  rtn=data.table(rtn)
-  names(rtn)=dimnames(cor)[[1]]
-  rtn=cbind(iter=seq(mc),rtn)
-
-  dat=melt(rtn,id="iter")
-
-  dat <- dat[, c("variable", "year"):= tstrsplit(variable, "_")]
-
-  setnames(dat, c("variable", "value"), c("qname", "data"))
-
-  res <- as(dat, 'FLQuants')
-
-  names(res)=c("stock","harvest")
-
-  return(res)
-}
-

data-raw/data.R

@@ -0,0 +1,16 @@
+# data.R - DESC
+# /home/mosqu003/FLR/pkgs/mine/ss3om/data-raw/data.R
+
+# Copyright (c) WUR, 2023.
+# Author: Iago MOSQUEIRA (WMR) <[email protected]>
+#
+# Distributed under the terms of the EUPL-1.2
+
+# XX {{{
+# }}}
+
+library(ss3om)
+
+sma <- readOutputss3('sma')
+
+save(sma, file='../data/sma.rda', compress='xz')