Some bug fixes

R/build.R

@@ -24,7 +24,7 @@ buildFLSss3 <- function(out, birthseas=out$birthseas, name=out$Control_File,
   # GET ages from catage
   ages <- getRange(out$catage)
   ages <- ac(seq(ages['min'], ages['max']))
-  dmns <- getDimnames(out, birthseas=birthseas)
+  dmns <- getDimnames(out)
   dim <- unlist(lapply(dmns, length))
 
   # EXTRACT from out
@@ -173,10 +173,6 @@ buildFLSRss3 <- function(out, ...) {
   lkhds <- out$likelihoods_used
   yrs <- recruit$Yr
 
-  # SET deviates and observed recruitment
-  recruit[, deviates := exp(dev)]
-  recruit[, obs_rec := rowSums(.SD, na.rm=TRUE), .SDcols=c("pred_recr", "deviates")]
-
   # logLik
   logLik <- lkhds[rownames(lkhds) == "Recruitment", "values"]
   class(logLik) <- "logLik"
@@ -190,6 +186,7 @@ buildFLSRss3 <- function(out, ...) {
       s=rawp[Label == "SR_BH_steep", Value],
       R0=exp(rawp[Label == "SR_LN(R0)", Value]),
       v=dquants[Label == "SSB_Virgin", Value],
+      sratio=1 / out$nsexes,
       units=c("", "1000", "t"))
     model <- "bevholtss3"
     attr(logLik, "df") <- length(rawp[!is.na(Active_Cnt), Active_Cnt])
@@ -206,19 +203,18 @@ buildFLSRss3 <- function(out, ...) {
   }
 
   # TODO SETUP for multiple recruit season/units
+
+  # LOAD FLQuants
+  dms <- list(year=yrs, season=out$spawnseas)
 
   # rec
-  rec <- FLQuant(recruit$obs_rec, dimnames=list(year=yrs, season=out$spawnseas),
-    units="1000")
+  rec <- FLQuant(recruit$pred_rec, dimnames=c(age=0, dms), units="1000")
   # ssb
-  ssb <- FLQuant(recruit$SpawnBio, dimnames=list(year=yrs, season=out$spawnseas),
-    units="t")
+  ssb <- FLQuant(recruit$SpawnBio, dimnames=c(age="all", dms), units="t")
   # fitted
-  fitted <- FLQuant(recruit$pred_recr, dimnames=list(year=yrs, season=out$spawnseas),
-    units="1000")
+  fitted <- FLQuant(recruit$bias_adjusted, dimnames=c(age=0, dms), units="1000")
   # residuals
-  residuals <- FLQuant(recruit$dev, dimnames=list(year=yrs, season=out$spawnseas),
-    units="")
+  residuals <- FLQuant(recruit$dev, dimnames=c(age=0, dms), units="")
 
   # vcov
   estpar <- parameters[grepl("SR_", Label),][!is.na(Active_Cnt),]
@@ -228,7 +224,6 @@ buildFLSRss3 <- function(out, ...) {
     vcov <- array((estpar$Parm_StDev)^2, dim=c(1,1,1),
       dimnames=list(estpar$Label, estpar$Label, iter=1))
   else if(dim(estpar)[1] > 1) {
-    # TODO
     CoVar <- data.table(out$CoVar)
     if(sum(dim(CoVar)) > 0)
       vcov <- CoVar[label.i %in% estpar$Label & label.j %in% estpar$Label,]
@@ -296,7 +291,7 @@ buildFLIBss3 <- function(out, fleets, birthseas=out$birthseas, ...) {
   index.res <- ss3index.res(cpue, fleets)
 
   # --- catch.n
-  catch <- ss3catch(catage, wtatage, dmns=getDimnames(out, birthseas=birthseas),
+  catch <- ss3catch(catage, wtatage, dmns=getDimnames(out),
     birthseas=birthseas, idx=fleets)
   catch.n <- lapply(catch, "[[", "landings.n")
 
@@ -431,7 +426,7 @@ buildFLBFss3 <- function(out, birthseas=unique(out$natage$BirthSeas)) {
   ages <- ac(seq(range['min'], range['max']))
 
   # GET dimnames
-  dmns <- getDimnames(out, birthseas=birthseas)
+  dmns <- getDimnames(out)
   dim <- unlist(lapply(dmns, length))
 
   # --- BIOL (endgrowth, natage)

R/build330.R

@@ -263,6 +263,93 @@ buildFLIBss330 <- function(out, fleets, birthseas=out$birthseas, ...) {
 
 } # }}}
 
+# buildFLIBss330 - FLIndices(FLIndexBiomass) {{{
+buildFLIss330 <- function(out, fleets, birthseas=out$birthseas, ...) {
+
+browser()
+  # SUBSET from out
+  out <- out[c("cpue", "ageselex", "endgrowth", "catage", "definitions",
+    "nsexes", "nseasons", "nareas", "birthseas")]
+
+  cpue <- data.table(out[["cpue"]])
+
+  # GET cpue fleets 
+  cpuefleets <- setNames(seq(length(unique(cpue$Fleet_name))), unique(cpue$Fleet_name))
+
+  if(missing(fleets))
+    fleets <- cpuefleets
+  else {
+    if(is.character(fleets))
+      fleets <- cpuefleets[names(cpuefleets) %in% fleets]
+    else if(is.numeric(fleets))
+      fleets <- cpuefleets[fleets]
+
+    # STOP if wrong fleets
+    if(length(fleets) == 0 | any(is.na(fleets)))
+      stop("selected fleets not found in Report.sso file")
+  }
+
+  selex <- data.table(out[["ageselex"]])
+  endgrowth <- data.table(out[["endgrowth"]])
+
+  # SET Age and unit
+  endgrowth[, Age:=int_Age]
+  endgrowth[, unit:=codeUnit(Sex, Platoon)]
+
+  wtatage <- endgrowth[,
+    c("Seas", "unit", "Age", paste0("RetWt:_", fleets)), with=FALSE]
+  catage <- data.table(out[["catage"]])
+    setkey(catage, "Area", "Fleet", "Sex", "Morph", "Yr", "Seas", "Era")
+  definitions <- data.table(out$definitions)
+
+  # --- index
+  index <- ss3index(cpue, fleets)
+
+  # --- index.q
+  index.q <- ss3index.q(cpue, fleets)
+
+  # --- sel.pattern
+  sel.pattern <- ss3sel.pattern(selex, unique(cpue$Yr), fleets,
+    morphs=unique(selex$Morph))
+
+  # --- index.var (var)
+  index.var <- ss3index.var(cpue, fleets)
+
+  # --- index.res (var)
+  index.res <- ss3index.res(cpue, fleets)
+
+  # --- catch.n
+  #catch <- ss3catch30(catage, wtatage, dmns=getDimnames(out),
+  #  birthseas=birthseas, idx=fleets)
+  #catch.n <- lapply(catch, "[[", "landings.n")
+
+  # --- FLIndices
+  cpues <- lapply(names(fleets), function(x) {
+
+    dmns <- dimnames(index[[x]])
+
+    selpat <- window(sel.pattern[[x]], start=dims(index[[x]])$minyear,
+        end=dims(index[[x]])$maxyear)
+
+    FLIndexBiomass(name=x,
+      distribution="lnorm", 
+      index=index[[x]],
+      index.q=index.q[[x]],
+      index.var=index.res[[x]],
+      # TODO How to link each cpue fleet to catch fleets for catch.n
+      # TRIM catch.n to index seasons
+   #   catch.n=unitSums(window(catch.n[[x]], start=dims(index[[x]])$minyear,
+   #     end=dims(index[[x]])$maxyear))[,,,dmns$season],
+      # NORMALIZE sel.pattern
+      sel.pattern=selpat %/% apply(selpat, 2:6, max))[,,,dmns$season]
+    })
+
+  names(cpues) <- names(fleets)
+
+  return(FLIndices(cpues))
+
+} # }}}
+
 # buildFLRPss330 - FLPar {{{
 buildFLRPss330 <- function(out, ...) {
 
@@ -315,7 +402,7 @@ buildRESss330 <- function(out, ...) {
     `Totbio_unfished`=out$derived_quants["Totbio_unfished", "Value"],
 
     # SSB_Unfished
-    `SSB_unfished`=out$derived_quants["Totbio_unfished", "Value"],
+    `SSB_unfished`=out$derived_quants["SSB_unfished", "Value"],
 
     # SSB_Virgin
     `SSB_Virgin`=out$derived_quants["SSB_Virgin", "Value"],
@@ -365,13 +452,13 @@ buildRESss330 <- function(out, ...) {
   return(res)
 } # }}}
 
-# buildFLBss330 {{{
+# buildFLBss330 - FLBiol {{{
 
-buildFLBss330 <- function(out, birthseas=out$birthseas, name=out$Control_File,
+buildFLBss330 <- function(out, morphs=out$morph_indexing$Index, name=out$Control_File,
   desc=paste(out$inputs$repfile, out$SS_versionshort, sep=" - ")) {
 
   # BUILD stk and srr
-  stk <- buildFLSss330(out, birthseas=birthseas, name=name, desc=desc)
+  stk <- buildFLSss330(out, morphs=morphs, name=name, desc=desc)
   srr <- buildFLSRss3(out)
 
   # COERCE to FLBiol
@@ -395,7 +482,7 @@ buildKobess330 <- function(out, ...) {
 
 } # }}}
 
-# buildFLBFss330 - FLStock {{{
+# buildFLBFss330 - FLBiol + FLFisheries {{{
 
 buildFLBFss330 <- function(out, morphs=out$morph_indexing$Index, name=out$Control_File,
   desc=paste(out$inputs$repfile, out$SS_versionshort, sep=" - "),
@@ -410,7 +497,8 @@ buildFLBFss330 <- function(out, morphs=out$morph_indexing$Index, name=out$Contro
     "FleetNames", "birthseas", "spawnseas", "inputs", "SS_versionshort", "parameters",
     "discard", "discard_at_age", "catch", "NatMort_option", "GrowthModel_option",
     "Maturity_option", "Fecundity_option", "Z_at_age", "M_at_age", "derived_quants",
-    "mean_body_wt", "Spawn_seas", "Spawn_timing_in_season", "morph_indexing")]
+    "mean_body_wt", "Spawn_seas", "Spawn_timing_in_season", "morph_indexing",
+    "exploitation")]
 
   # GET ages from catage
   ages <- getRange(out$catage)
@@ -491,15 +579,25 @@ buildFLBFss330 <- function(out, morphs=out$morph_indexing$Index, name=out$Contro
   selex <- lapply(lapply(fleets, function(x)
     as.FLQuant(ageselex[fleet %in% x,][, fleet:=NULL], units="NA")),
     window, start=dmns$year[1])
+
+  # effort
+  expl <- data.table(out$exploitation)
+
+  effs <- expl[Seas == 1, c("Yr", out$FleetNames[out$fleet_type == 1]), with=FALSE]
+
+  effqs <- FLQuants(lapply(colnames(fs)[-1], function(x) {
+    y <- fs[, c("Yr", x), with=FALSE]
+    setnames(y, c("year", "data"))
+    y[, effort:="all"]
+    as.FLQuant(y)
+  }))
 
-  flfs <- Map(function(c, s) {
+  flfs <- Map(function(c, s, ef) {
     ca <- FLCatch(landings.n=c$catch.n, landings.wt=c$catch.wt, catch.sel=s)
     discards.n(ca)[] <- 0
     discards.wt(ca) <- landings.wt(ca)
-    return(FLFishery(
-      effort=FLQuant(0, dimnames=list(effort='all', year=dimnames(c$catch.n)$year)),
-      A=ca))
-    }, c=catches, s=selex)
+    return(FLFishery(effort=ef, A=ca))
+    }, c=catches, s=selex, ef=effqs)
 
   # TABLE of areas and fleets
   map <- unique(catage[, .(Area, Fleet)])

R/check.R

@@ -51,6 +51,7 @@ extractSSB <- function(out) {
 }
 
 #' @rdname extractSS
+
 extractRec <- function(out) {
 
   rec <- data.table(out$derived_quants)[Label %in% paste0("Recr_",
@@ -89,6 +90,8 @@ extractFbar <- function(out) {
     year=seq(out$startyr + 1, out$endyr), unit="unique"), units="f"))
 }
 
+#' @rdname extractSS
+
 extractZatage <- function(out) {
 
   # OUTPUT years
@@ -119,5 +122,14 @@ extractZatage <- function(out) {
   return(z)
 }
 
+#' @rdname extractSS
+
+extractDevs <- function(out) {
+
+  dev <- data.table(out$recruit)[era == "Main", .(Yr, dev)]
 
+  dms <- getRange(out$catage)
+
+  return(FLQuant(dev$dev, dimnames=list(age=dms['min'], year=dev$Yr), units=""))
+}
 # }}}

R/load.R

@@ -49,7 +49,7 @@ loadOMS <- function(dir=".", subdirs=list.dirs(path=dir, recursive=FALSE),
 	# LOOP over subdirs
   out <- foreach(i=seq(length(subdirs)),
     .final = function(x) setNames(x, nm=seq(length(subdirs))),
-    .inorder=TRUE, .errorhandling="remove") %dopar% {
+    .inorder=TRUE, .errorhandling="pass") %dopar% {
 
     if(progress)
       cat("[", i, "]\n", sep="")
@@ -69,6 +69,14 @@ loadOMS <- function(dir=".", subdirs=list.dirs(path=dir, recursive=FALSE),
     run
   }
 
+  # CHECK for errors (NULL)
+  nulls <- unlist(lapply(out, function(x)
+    any(is.null(unlist(lapply(x, is.null))))))
+
+  if(any(nulls))
+    stop(paste("Some iters returned one or more NULL elements:",
+      paste(unname(which(nulls)), collapse=", ")))
+
   if(progress)
     cat("[combining now ...]\n", sep="")
 
@@ -199,7 +207,7 @@ loadRES <- function(dir=".", subdirs=list.dirs(path=dir, recursive=FALSE),
 
 	# LOOP over subdirs
   out <- foreach(i=iters, .inorder=TRUE,
-    .errorhandling="remove") %dopar% {
+    .errorhandling="pass") %dopar% {
 
     if(progress)
       cat("[", i, "]\n", sep="")
@@ -208,11 +216,17 @@ loadRES <- function(dir=".", subdirs=list.dirs(path=dir, recursive=FALSE),
     readRESss3(subdirs[i], ...)
   }
 
-  out <- rbindlist(out)
+  errs <- !unlist(lapply(out, is, "data.table"))
+
+  if(any(errs)) {
+    warning(paste("Some results could not be loaded:",
+    paste(which(errs), collapse=", ")))
+  }
+
+  out <- rbindlist(out[!errs])
 
   if(!is.null(grid)) {
-    out <- cbind(data.table(grid), out)
-    out[, iter := iters]
+    out <- cbind(data.table(grid)[!errs], out)
   }
 
   return(out)

R/read.R

@@ -73,11 +73,11 @@ readOutputss3 <- function(dir, repfile = "Report.sso",
 
   # CHECK compressed files
   idx <- file.exists(file.path(dir, paste(cfiles, compress, sep = ".")))
-  cfiles[idx] <- paste(cfiles, compress, sep = ".")
+  cfiles[idx] <- paste(cfiles[idx], compress, sep = ".")
 
   out <- SS_output(dir, verbose=FALSE, hidewarn=TRUE, warn=FALSE,
     printstats=FALSE, covarfile=cfiles["covarfile"], forecast=FALSE,
-    repfile=cfiles["repfile"], compfile=cfiles["compfile"])
+    repfile=cfiles["repfile"], compfile=cfiles["compfile"], covar=idx[3])
 
   return(out) 
 } # }}}
@@ -236,7 +236,7 @@ readFLSss3 <- function(dir, repfile="Report.sso", compfile="CompReport.sso",
 readRESss3 <- function(dir, repfile="Report.sso", compfile="CompReport.sso", ...) {
 
   # LOAD SS_output list
-  out <- readOutputss3(dir, repfile=repfile, compfile=compfile)
+  out <- readOutputss3(dir, repfile=repfile, compfile=compfile, covar=FALSE)
 
   if(out$SS_versionNumeric > 3.24)
     buildRESss330(out, ...)

R/ss3slots30.R

@@ -128,6 +128,8 @@ ss3catch30 <- function(catage, wtatage, dmns, birthseas, idx) {
   catage[, Area := if(length(unique(Area)) == 1) "unique" else Area]
 
   # MELT by Sex, BirthSeas, Yr & Seas
+  catage[, (dmns$age) := lapply(.SD, as.double), .SDcols = dmns$age]
+
 	catage <- data.table::melt(catage, id.vars=c("Area", "Fleet", "Yr", "Seas", "unit"),
     measure.vars=dmns$age, variable.name="age")
   names(catage) <- c("area", "fleet", "year", "season", "unit", "age", "data")