Skip to content

Commit 1b2224f

Browse files
andybeetandybeet
authored
Merge pull request #162 from NOAA-EDAB/pb_fixes
Pb fixes
2 parents effaa04 + de75fe1 commit 1b2224f

2 files changed

Lines changed: 55 additions & 21 deletions

File tree

‎R/ecopath.R‎

Lines changed: 31 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
#'@export
1919
rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
2020
#Need to define variables to eliminate check() note about no visible binding
21-
Type <- Group <- DetInput <- ProdCons <- PB <- QB <- noB <- noEE <- alive <- NULL
21+
Type <- Group <- DetInput <- ProdCons <- PB <- QB <- noB <- noEE <- alive <- noPB <- NULL
2222
BEE <- Biomass <- Q <- BioAcc <- BioQB <- diag.a <- EEa <- B <- M0 <- NULL
2323
QBloss <- Unassim <- Ex <- NULL
2424

@@ -62,7 +62,9 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
6262
model[is.na(DetInput), DetInput := 0]
6363

6464
# fill in GE(PQ), QB, or PB from other inputs
65-
GE <- ifelse(is.na(model[, ProdCons]), model[, PB / QB], model[, ProdCons])
65+
# KYA Aug 2025 - changed this logic so PC would be recalculated if PB and QB supplied
66+
#GE <- ifelse(is.na(model[, ProdCons]), model[, PB / QB], model[, ProdCons])
67+
GE <- ifelse(!is.na(model[, QB]) & !is.na(model[, PB]), model[, PB / QB], model[, ProdCons])
6668
QB.1 <- ifelse(is.na(model[, QB]), model[, PB / GE], model[, QB])
6769
PB.1 <- ifelse(is.na(model[, PB]), model[, ProdCons * QB], model[, PB])
6870
model[, QB := QB.1]
@@ -97,11 +99,16 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
9799
model[, noEE := 0]
98100
model[, alive := 0]
99101
model[, BEE := 0]
102+
model[, noPB := 0]
100103
model[is.na(Biomass), noB := 1]
101104
model[is.na(EE), noEE := 1]
102105
model[Type < 2, alive := 1]
103-
model[noB == 0 & noEE == 0, BEE := 1]
104-
106+
model[noB == 0 & noEE == 0, BEE := 1]
107+
model[BEE == 1 & is.na(PB), noPB := 1]
108+
109+
if (any(model$Type==0 & is.na(model$QB) & is.na(model$ProdCons))){
110+
stop("A consumer is missing both QB and ProdCons - balance failed. Use check.rpath.params() to diagnose.")
111+
}
105112
# define detritus fate matrix
106113
detfate <- model[, (10 + 1):(10 + ndead), with = F]
107114
detdetfate <- model[Type==2, (10 + 1):(10 + ndead), with = F]
@@ -118,6 +125,7 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
118125
#Set up A matrix
119126
living[noEE == 1, diag.a := Biomass * PB]
120127
living[noEE == 0, diag.a := PB * EE]
128+
living[noPB == 1, diag.a := Biomass * EE] # this needs to be after noEE==0 case
121129

122130
#Special case where B and EE are known then need to solve for BA
123131
#living[BEE == 1, b := b - (Biomass * PB * EE)]
@@ -134,6 +142,11 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
134142
A <- A - QBDCa
135143
#Switch flag back
136144
#living[BEE == 1, noB := 0]
145+
146+
# Check for any missing info that will prevent solving
147+
if (any(is.na(A))){
148+
stop("Model is missing parameters - can't be balanced. Use check.rpath.params() to diagnose.")
149+
}
137150

138151
# Generalized inverse does the actual solving
139152
#Invert A and multiple by b to get x (unknowns)
@@ -146,6 +159,9 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
146159
living[, B := x * noB]
147160
living[is.na(Biomass), Biomass := B]
148161

162+
living[, PBa := x * noPB]
163+
living[is.na(PB), PB := PBa]
164+
149165
# detritus EE calcs
150166
living[, M0 := PB * (1 - EE)]
151167
living[, QBloss := QB]
@@ -182,7 +198,7 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
182198
inDetB <- model[(nliving + 1):(nliving + ndead), Biomass]
183199
DetPB <- ifelse(is.na(inDetPB), Default_Detrital_PB, inDetPB)
184200
DetB <- ifelse(is.na(inDetB), detinputs / DetPB, inDetB)
185-
DetPB <- detinputs / DetB
201+
DetPB <- as.numeric(detinputs) / DetB
186202

187203
# Trophic Level calcs
188204
b <- rep(1, ngroups)
@@ -217,16 +233,20 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) {
217233
#to match header file format (replacing NAs with 0.0s)
218234
Bplus <- c(living[, Biomass], DetB, rep(0.0, ngear))
219235

220-
PBplus <- model[, PB]
221-
PBplus[(nliving + 1):(nliving + ndead)] <- DetPB
236+
#PBplus <- model[, PB]
237+
#PBplus[(nliving + 1):(nliving + ndead)] <- DetPB
238+
PBplus <- c(living[,PB],DetPB , rep(0.0, ngear))
222239
PBplus[is.na(PBplus)] <- 0.0
223240

224241
EEplus <- c(EE, rep(0.0, ngear))
225242

226243
QBplus <- model[, QB]
244+
QBplus[is.na(QBplus) & PBplus>0.0 & !(is.na(GE) | is.nan(GE) | is.infinite(GE))] <-
245+
(PBplus/GE)[is.na(QBplus) & PBplus>0.0 & !(is.na(GE) | is.nan(GE) | is.infinite(GE))]
227246
QBplus[is.na(QBplus)] <- 0.0
228247

229-
GE[is.na(GE)] <- 0.0
248+
GE <- PBplus/QBplus
249+
GE[is.na(GE) | is.nan(GE) | is.infinite(GE)] <- 0.0
230250

231251
RemPlus <- model[, totcatch]
232252
RemPlus[is.na(RemPlus)] <- 0.0
@@ -337,6 +357,9 @@ rpath.stanzas <- function(Rpath.params){
337357
Group <- Biomass <- R <- NageS <- bs.denom <- bs <- qs.denom <- qs <- Cons <- NULL
338358
QB <- BAB <- Ex <- NULL
339359

360+
# Added Aug 2025 - if no stanzas, silently return original (prob no warning needed?)
361+
if(Rpath.params$stanza$NStanzaGroups==0){return(Rpath.params)}
362+
340363
#Determine the total number of groups with multistanzas
341364
Nsplit <- Rpath.params$stanza$NStanzaGroups
342365
groupfile <- Rpath.params$stanza$stgroups

‎R/param.R‎

Lines changed: 24 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,7 @@ check.rpath.params <- function(Rpath.params) {
172172
Type <- Group <- Biomass <- EE <- PB <- QB <- ProdCons <- BioAcc <- Unassim <- DetInput <- NULL
173173

174174
w <- 0 #warning counter
175+
c <- 0 #Balance change counter
175176
#Check to make sure all types are represented
176177
if (nrow(Rpath.params$model[Type == 0, ]) == 0) {
177178
warning('Model must contain at least 1 consumer')
@@ -215,13 +216,16 @@ check.rpath.params <- function(Rpath.params) {
215216
)
216217
w <- w + 1
217218
}
218-
if (length(Rpath.params$model[!is.na(Biomass) &
219-
!is.na(EE) & Type < 2, Group]) > 0) {
219+
if (length(Rpath.params$model[!is.na(Biomass) & !is.na(EE) &
220+
(!is.na(PB) | (is.na(PB) & !is.na(QB) & !is.na(ProdCons))) &
221+
Type < 2, Group]) > 0) {
220222
warning(
221223
paste(
222-
Rpath.params$model[!is.na(Biomass) & !is.na(EE) & Type < 2, Group],
223-
'have both Biomass and EE...Note that Rpath does not calculate BA
224-
please enter a value for BA if appropriate \n',
224+
Rpath.params$model[!is.na(Biomass) & !is.na(EE) &
225+
(!is.na(PB) | (is.na(PB) & !is.na(QB) & !is.na(ProdCons))) &
226+
Type < 2, Group],
227+
'have all of Biomass, EE, and PB(or QB and ProdCons) entered... Note that Rpath does
228+
not calculate BA, please enter a value for BA if appropriate \n',
225229
sep = ' '
226230
)
227231
)
@@ -255,7 +259,7 @@ check.rpath.params <- function(Rpath.params) {
255259
warning(
256260
paste(
257261
Rpath.params$model[Type > 1 & !is.na(QB), Group],
258-
'are not living and should not have a QB...set to NA \n',
262+
'are not living and should not have a QB... please set to NA \n',
259263
sep = ' '
260264
)
261265
)
@@ -265,7 +269,7 @@ check.rpath.params <- function(Rpath.params) {
265269
warning(
266270
paste(
267271
Rpath.params$model[Type > 1 & !is.na(EE), Group],
268-
'are not living and should not have a EE...set to NA \n',
272+
'are not living and should not have a EE... please set to NA \n',
269273
sep = ' '
270274
)
271275
)
@@ -276,7 +280,7 @@ check.rpath.params <- function(Rpath.params) {
276280
warning(
277281
paste(
278282
Rpath.params$model[Type > 1 & !is.na(ProdCons), Group],
279-
'are not living and should not have a ProdCons...set to NA \n',
283+
'are not living and should not have a ProdCons... please set to NA \n',
280284
sep = ' '
281285
)
282286
)
@@ -287,12 +291,14 @@ check.rpath.params <- function(Rpath.params) {
287291
if (length(Rpath.params$model[Type < 2 & is.na(PB), Group]) > 0) {
288292
no.pb <- Rpath.params$model[Type < 2 & is.na(PB), Group]
289293
if (length(Rpath.params$model[Group %in% no.pb &
290-
(is.na(QB) | is.na(ProdCons)), Group]) > 0) {
294+
(is.na(QB) | is.na(ProdCons)) &
295+
(is.na(Biomass) | is.na(EE)), Group]) > 0) {
291296
warning(
292297
paste(
293298
Rpath.params$model[Group %in% no.pb &
294-
(is.na(QB) | is.na(ProdCons)), Group],
295-
'are missing a PB without a QB and PQ...set to >= 0 \n',
299+
(is.na(QB) | is.na(ProdCons)) &
300+
(is.na(Biomass) | is.na(EE)), Group],
301+
'are missing a PB without either a (QB and ProdCons) or (EE and B) to estimate PB... please set to >= 0 \n',
296302
sep = ' '
297303
)
298304
)
@@ -321,11 +327,12 @@ check.rpath.params <- function(Rpath.params) {
321327
warning(
322328
paste(
323329
Rpath.params$model[Group %in% both & !is.na(PB), Group],
324-
'have PB, QB, and ProdCons...only two should be entered \n',
330+
'have PB, QB, and ProdCons... ProdCons will be recalculated during balancing \n',
325331
sep = ' '
326332
)
327333
)
328334
w <- w + 1
335+
c <- c + 1
329336
}
330337
}
331338

@@ -516,7 +523,11 @@ check.rpath.params <- function(Rpath.params) {
516523
if (w == 0) {
517524
cat('Rpath parameter file is functional. \n')
518525
} else {
519-
cat('Rpath parameter file needs attention! \n')
526+
if (w==c){
527+
cat('Rpath parameters functional, though some may be recalculated during balance. \n')
528+
} else {
529+
cat('Rpath parameter file needs attention! \n')
530+
}
520531
}
521532
}
522533

0 commit comments

Comments
 (0)