added/tested PB estimation feature to match EwE

Author: kerim aydin

R/ecopath.R

@@ -18,7 +18,7 @@
 #'@export
 rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
   #Need to define variables to eliminate check() note about no visible binding
-  Type <- Group <- DetInput <- ProdCons <- PB <- QB <- noB <- noEE <- alive <- NULL
+  Type <- Group <- DetInput <- ProdCons <- PB <- QB <- noB <- noEE <- alive <- noPB <- NULL
   BEE <- Biomass <- Q <- BioAcc <- BioQB <- diag.a <- EEa <- B <- M0 <- NULL
   QBloss <- Unassim <- Ex <- NULL
 
@@ -99,11 +99,13 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
   model[, noEE  := 0]
   model[, alive := 0]
   model[, BEE   := 0]
+  model[, noPB  := 0]
   model[is.na(Biomass), noB   := 1]
   model[is.na(EE),      noEE  := 1]
   model[Type < 2,       alive := 1]
-  model[noB == 0 & noEE == 0, BEE := 1]
-  
+  model[noB == 0 & noEE == 0,  BEE := 1]
+  model[BEE == 1 & is.na(PB), noPB := 1]    
+          
   # define detritus fate matrix
   detfate <- model[, (10 + 1):(10 + ndead), with = F]
   detdetfate <- model[Type==2, (10 + 1):(10 + ndead), with = F]
@@ -120,6 +122,7 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
   #Set up A matrix
   living[noEE == 1, diag.a := Biomass * PB]
   living[noEE == 0, diag.a := PB * EE]
+  living[noPB == 1, diag.a := Biomass * EE] # this needs to be after noEE==0 case
 
   #Special case where B and EE are known then need to solve for BA
   #living[BEE == 1, b := b - (Biomass * PB * EE)]
@@ -148,6 +151,9 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
   living[, B := x * noB]
   living[is.na(Biomass), Biomass := B]
 
+  living[, PBa := x * noPB]
+  living[is.na(PB), PB := PBa]
+  
   # detritus EE calcs
   living[, M0 := PB * (1 - EE)]
   living[, QBloss := QB]
@@ -184,7 +190,7 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
   inDetB  <- model[(nliving + 1):(nliving + ndead), Biomass]
   DetPB   <- ifelse(is.na(inDetPB), Default_Detrital_PB, inDetPB)
   DetB    <- ifelse(is.na(inDetB), detinputs / DetPB, inDetB)
-  DetPB   <- detinputs / DetB
+  DetPB   <- as.numeric(detinputs) / DetB
 
   # Trophic Level calcs
   b             <- rep(1, ngroups)
@@ -219,16 +225,18 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
   #to match header file format (replacing NAs with 0.0s)
   Bplus  <- c(living[, Biomass], DetB, rep(0.0, ngear))
 
-  PBplus <- model[, PB] 
-  PBplus[(nliving + 1):(nliving + ndead)] <- DetPB
+  #PBplus <- model[, PB] 
+  #PBplus[(nliving + 1):(nliving + ndead)] <- DetPB
+  PBplus  <- c(living[,PB],DetPB , rep(0.0, ngear))
   PBplus[is.na(PBplus)] <- 0.0
 
   EEplus <- c(EE, rep(0.0, ngear))
 
   QBplus <- model[, QB]
   QBplus[is.na(QBplus)] <- 0.0
 
-  GE[is.na(GE)] <- 0.0
+  GE <- PBplus/QBplus
+  GE[is.na(GE) | is.nan(GE) | is.infinite(GE)] <- 0.0
 
   RemPlus <- model[, totcatch]
   RemPlus[is.na(RemPlus)] <- 0.0

R/param.R

@@ -216,12 +216,14 @@ check.rpath.params <- function(Rpath.params) {
     w <- w + 1
   }
   if (length(Rpath.params$model[!is.na(Biomass) &
-                                !is.na(EE) & Type < 2, Group]) > 0) {
+                                !is.na(EE) & 
+                                (!is.na(PB) | (is.na(PB) & !is.na(QB) & !is.na(ProdCons))) &
+                                               Type < 2, Group]) > 0) {
     warning(
       paste(
         Rpath.params$model[!is.na(Biomass) & !is.na(EE) & Type < 2, Group],
-        'have both Biomass and EE...Note that Rpath does not calculate BA
-               please enter a value for BA if appropriate \n',
+        'have all of Biomass, EE, and PB entered... Note that Rpath does
+        not calculate BA, please enter a value for BA if appropriate \n',
         sep = ' '
       )
     )
@@ -255,7 +257,7 @@ check.rpath.params <- function(Rpath.params) {
     warning(
       paste(
         Rpath.params$model[Type > 1 & !is.na(QB), Group],
-        'are not living and should not have a QB...set to NA \n',
+        'are not living and should not have a QB... please set to NA \n',
         sep = ' '
       )
     )
@@ -265,7 +267,7 @@ check.rpath.params <- function(Rpath.params) {
     warning(
       paste(
         Rpath.params$model[Type > 1 & !is.na(EE), Group],
-        'are not living and should not have a EE...set to NA \n',
+        'are not living and should not have a EE... please set to NA \n',
         sep = ' '
       )
     )
@@ -276,7 +278,7 @@ check.rpath.params <- function(Rpath.params) {
     warning(
       paste(
         Rpath.params$model[Type > 1 & !is.na(ProdCons), Group],
-        'are not living and should not have a ProdCons...set to NA \n',
+        'are not living and should not have a ProdCons... please set to NA \n',
         sep = ' '
       )
     )
@@ -287,12 +289,14 @@ check.rpath.params <- function(Rpath.params) {
   if (length(Rpath.params$model[Type < 2 & is.na(PB), Group]) > 0) {
     no.pb <- Rpath.params$model[Type < 2 & is.na(PB), Group]
     if (length(Rpath.params$model[Group %in% no.pb &
-                                  (is.na(QB) | is.na(ProdCons)), Group]) > 0) {
+                                  (is.na(QB) | is.na(ProdCons)) &
+                                  (is.na(Biomass) | is.na(EE)), Group]) > 0) {
       warning(
         paste(
           Rpath.params$model[Group %in% no.pb &
-                               (is.na(QB) | is.na(ProdCons)), Group],
-          'are missing a PB without a QB and PQ...set to >= 0 \n',
+                               (is.na(QB) | is.na(ProdCons)) &
+                               (is.na(Biomass) | is.na(EE)), Group],
+          'are missing a PB without either a (QB and ProdCons) or (EE and B) to estimate PB... please set to >= 0 \n',
           sep = ' '
         )
       )