diff --git a/.Rbuildignore b/.Rbuildignore index 525d30c5..a27befab 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,3 +9,7 @@ ^doc$ ^Meta$ +# files/folders +CODE_OF_CONDUCT.md +CONTRIBUTING.md +testthat \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 18cc70b0..bab99b35 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: Rpath Type: Package Title: R implementation of Ecopath with Ecosim -Version: 1.0.0 +Version: 1.1.0 Authors@R: c( person("Kerim", "Aydin", email = "kerim.aydin@noaa.gov", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-3792-9828")), @@ -36,7 +36,9 @@ Imports: stats, utils, ggrepel, - ggplot2 + ggplot2, + janitor, + xml2 LinkingTo: Rcpp Suggests: here, @@ -51,6 +53,7 @@ Suggests: usethis, dplyr, generics, + DT, testthat (>= 3.0.0) VignetteBuilder: knitr RoxygenNote: 7.3.2 diff --git a/NAMESPACE b/NAMESPACE index d3fc4a57..88c58368 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -11,6 +11,7 @@ export(adjust.fishing) export(adjust.forcing) export(adjust.scenario) export(check.rpath.params) +export(create.rpath.from.eiixml) export(create.rpath.params) export(extract.node) export(frate.table) @@ -20,6 +21,7 @@ export(get.rsim.params) export(get.rsim.stanzas) export(get.rsim.start_state) export(ggwebplot) +export(import.eiixml) export(read.rpath.params) export(rpath) export(rpath.consumers) diff --git a/NEWS.md b/NEWS.md index 6292a7fa..8b38a701 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,20 @@ +# Rpath 1.1.0 + +- Added ability to import .eiixml files via the `create.rpath.from.eiixml()` function +- Added balance functionality to estimate P/B from input Biomass and EE, also +allowing a missing Q/B to be estimated from input PC and the estimated P/B +- Corrected accounting for single-stage interdetrital flows during balance +- Behavior change: if all of P/B, Q/B, and PC are supplied as inputs for a group, recalculate PC during balance to ensure internal model consistency (instead of leaving PC unchanged and possibly inaccurate) +- Improved error messages in `rpath()` balance routine to help diagnose models that are missing parameters +- Improved warning message content in `check.rpath.params()` to aid diagnosis +- Modified `rpath.stanzas()` so it no longer produces errors when called with a model that has no multistanza groups (instead it returns the model unchanged) +- Added Western Bering Sea model .eiixml and EwE output csv files as an import example +- Added `Convert_EwE_to_Rpath` vignette to describe .eiixml import process +- Tested eiixml input routines on over 150 models available via EcoBase with most producing consistent results between Rpath and EwE; some remaining balancing differences between Rpath and EwE are noted in the `Convert_EwE_to_Rpath` vignette +- Fixed some table ordering issues with stanza inputs +- Documentation fixes + + # Rpath 1.0.0 - Fix group parameter reference in `adjust.scenario` diff --git a/R/Adjustments.R b/R/Adjustments.R index ac43826d..0b6200e2 100644 --- a/R/Adjustments.R +++ b/R/Adjustments.R @@ -1,6 +1,6 @@ #'Fishing Mortality Table #' -#'Creates a table of fishing mortalities by species group and gear for an +#'Creates a table of fishing mortalities by species group and gear for an #'\code{rsim.scenario()} object. #' #'@family Rpath functions @@ -24,27 +24,30 @@ #' #'@export - - -frate.table <- function(Rsim.scenario){ +frate.table <- function(Rsim.scenario) { #Need to define variables to eliminate check() note about no visible binding Group <- Gear <- Q <- NULL - - fish <- data.table(Group = Rsim.scenario$params$FishFrom, - Gear = Rsim.scenario$params$FishThrough, - Q = Rsim.scenario$params$FishQ) + + fish <- data.table( + Group = Rsim.scenario$params$FishFrom, + Gear = Rsim.scenario$params$FishThrough, + Q = Rsim.scenario$params$FishQ + ) group <- unique(fish[Group > 0, Group]) - gear <- unique(fish[Gear > 0, Gear]) - + gear <- unique(fish[Gear > 0, Gear]) + group.name <- Rsim.scenario$params$spname[unique(fish[Group > 0, Group]) + 1] - gear.name <- Rsim.scenario$params$spname[unique(fish[Gear > 0, Gear ]) + 1] - + gear.name <- Rsim.scenario$params$spname[unique(fish[Gear > 0, Gear]) + 1] + fish.out <- c() - for(i in 1:length(group)){ + for (i in 1:length(group)) { fish.group <- fish[Group == group[i], ] fish.all.gear <- data.table(Group = group.name[i]) - for(j in 1:length(gear)){ - f.gear <- data.table(Group = group.name[i], V1 = fish.group[Gear == gear[j], sum(Q)]) + for (j in 1:length(gear)) { + f.gear <- data.table( + Group = group.name[i], + V1 = fish.group[Gear == gear[j], sum(Q)] + ) setnames(f.gear, 'V1', gear.name[j]) fish.all.gear <- merge(fish.all.gear, f.gear, by = 'Group') } @@ -57,9 +60,9 @@ frate.table <- function(Rsim.scenario){ #'Adjust Fishing Mortality #' -#'Modifies the fishing mortality value for a species by a particular gear. +#'Modifies the fishing mortality value for a species by a particular gear. #'Parameters that can be adjusted using this function are: \emph{ForcedEffort}, -#'\emph{ForcedFRate}, or \emph{ForcedCatch}. +#'\emph{ForcedFRate}, or \emph{ForcedCatch}. #' #'@family Adjust functions #' @@ -72,68 +75,78 @@ frate.table <- function(Rsim.scenario){ #' the year are modified. #'@param value New value for the parameter. #' -#'@return Returns an \code{Rsim.scenario()} object with the new fishing parameter +#'@return Returns an \code{Rsim.scenario()} object with the new fishing parameter #' values. -#' +#' #'@examples #' # Read in Rpath parameter file and generate balanced model #' Rpath <- rpath(AB.params) #' # Create a 50 yr Rsim scenario #' Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) #' # Change value of forcedFRate for Squids in years 3 through 5 to the value of 2 (for all months) -#' Rsim.scenario.adjusted.fishing <- adjust.fishing(Rsim.scenario,parameter="ForcedFRate",group="cod",sim.year=3:5,value = 2) -#' head(Rsim.scenario.adjusted.fishing$fishing$ForcedFRate) -#' -#' +#' Rsim.scenario.adjusted.fishing <- adjust.fishing(Rsim.scenario, parameter = "ForcedFRate", group = "cod", sim.year = 3:5, value = 2) +#' head(Rsim.scenario.adjusted.fishing$fishing$ForcedFRate) +#' +#' #'@export -adjust.fishing <- function(Rsim.scenario, parameter, group = NA, sim.year = 1, - sim.month = 0, value){ +adjust.fishing <- function( + Rsim.scenario, + parameter, + group = NA, + sim.year = 1, + sim.month = 0, + value +) { #Check that parameter and group exist - if(!parameter %in% c('ForcedEffort', 'ForcedFRate', 'ForcedCatch')){stop("Fishing parameter not found")} - - if(!all(group %in% Rsim.scenario$params$spname)){ - stop("Groups not found:",group[!(group %in% Rsim.scenario$params$spname)]) + if (!parameter %in% c('ForcedEffort', 'ForcedFRate', 'ForcedCatch')) { + stop("Fishing parameter not found") + } + + if (!all(group %in% Rsim.scenario$params$spname)) { + stop("Groups not found:", group[!(group %in% Rsim.scenario$params$spname)]) } #if(!group %in% Rsim.scenario$params$spname){stop("Group not found")} - + #Create index in case the number of values is equal to length(sim.year) * length(sim.month) - ivalue <- 0 - + ivalue <- 0 + #Loop over years if more than 1 sim.year provided - for(iyear in seq_along(sim.year)){ - for(imonth in seq_along(sim.month)){ + for (iyear in seq_along(sim.year)) { + for (imonth in seq_along(sim.month)) { ivalue <- ivalue + 1 #look-up what rows correspond to the year #Regex used to account for year.month row names in Effort Matrix - year.row <- which(gsub("\\..*", "", rownames(Rsim.scenario$fishing[[parameter]])) - == sim.year[iyear]) - + year.row <- which( + gsub("\\..*", "", rownames(Rsim.scenario$fishing[[parameter]])) == + sim.year[iyear] + ) + #identify what rows correspond to the sim.months - 0 indicates the whole year - if(sim.month[1] != 0){ + if (sim.month[1] != 0) { year.row <- year.row[1:12 %in% sim.month[imonth]] - } - + } + #Apply the value to the correct row - if(length(value) == 1){ + if (length(value) == 1) { #If only 1 value is supplied for multiple years need to only point to that value Rsim.scenario$fishing[[parameter]][year.row, group] <- value - }else if(sim.month[1] == 0){ - Rsim.scenario$fishing[[parameter]][year.row, group] <- value[iyear] - }else if(length(value) == length(sim.month)){ - Rsim.scenario$fishing[[parameter]][year.row, group] <- value[imonth] - }else { - Rsim.scenario$fishing[[parameter]][year.row, group] <- value[ivalue] - } + } else if (sim.month[1] == 0) { + Rsim.scenario$fishing[[parameter]][year.row, group] <- value[iyear] + } else if (length(value) == length(sim.month)) { + Rsim.scenario$fishing[[parameter]][year.row, group] <- value[imonth] + } else { + Rsim.scenario$fishing[[parameter]][year.row, group] <- value[ivalue] } } + } return(Rsim.scenario) } - + #'Adjust Rsim.scenario parameters #' -#'Modifies the various parameters of the \code{rsim.scenario()} object. Parameters that can be adjusted using this function are: +#'Modifies the various parameters of the \code{rsim.scenario()} object. Parameters that can be adjusted using this function are: #'\emph{B_BaseRef}, \emph{MzeroMort},\emph{UnassimRespFrac}, \emph{ActiveRespFrac}, \emph{FtimeAdj}, #'\emph{FtimeQBOpt}, \emph{PBopt}, \emph{NoIntegrate},\emph{HandleSelf}, \emph{ScrambleSelf}, \emph{QQ}, #' \emph{DD}, \emph{VV}, \emph{HandleSwitch}, \emph{PredPredWeight}, \emph{PreyPreyWeight} @@ -142,10 +155,10 @@ adjust.fishing <- function(Rsim.scenario, parameter, group = NA, sim.year = 1, #' #'@inheritParams adjust.fishing #' -#'@param parameter Parameters to be modified (Choose from: \code{B_BaseRef, MzeroMort, +#'@param parameter Parameters to be modified (Choose from: \code{B_BaseRef, MzeroMort, #' UnassimRespFrac, ActiveRespFrac, FtimeAdj, FtimeQBOpt, PBopt, NoIntegrate, #' HandleSelf, ScrambleSelf, QQ, DD, VV, HandleSwitch, PredPredWeight, PreyPreyWeight}) -#'@param group The model group that the parameter change will affect. Note that +#'@param group The model group that the parameter change will affect. Note that #' a value of \emph{'all'} will affect all groups associated with the `groupto` #' variable. Valid values are found in the `Group` field of the object created #' from running \code{rpath()} @@ -160,38 +173,63 @@ adjust.fishing <- function(Rsim.scenario, parameter, group = NA, sim.year = 1, #' # Create a 50 yr Rsim scenario #' Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) #' # Adjust the PBopt parameter for cod. Set to value = 2 -#' Rsim.scenario.adjusted <- adjust.scenario(Rsim.scenario, parameter="PBopt",group = "cod", groupto = "all", value = 2) +#' Rsim.scenario.adjusted <- adjust.scenario(Rsim.scenario, parameter="PBopt",group = "cod", groupto = "all", value = 2) #' #' -#'@export +#'@export -adjust.scenario <- function(Rsim.scenario, parameter, group, groupto = NA, value){ +adjust.scenario <- function( + Rsim.scenario, + parameter, + group, + groupto = NA, + value +) { #Lookup group numbers - if(group[1] == 'all'){ + if (group[1] == 'all') { groupnum <- 0:Rsim.scenario$params$NUM_GROUPS } else { - groupnum <- Rsim.scenario$params$spnum[which(Rsim.scenario$params$spname - %in% group)] + groupnum <- Rsim.scenario$params$spnum[which( + Rsim.scenario$params$spname %in% group + )] } - if(!is.na(groupto)){ - groupnumto <- Rsim.scenario$params$spnum[which(Rsim.scenario$params$spname - %in% groupto)] + if (!is.na(groupto)) { + groupnumto <- Rsim.scenario$params$spnum[which( + Rsim.scenario$params$spname %in% groupto + )] } - + #Lookup parameter number param.num <- which(names(Rsim.scenario$params) == parameter) - + #Modify parameter - if(parameter %in% c('B_BaseRef', 'MzeroMort', 'UnassimRespFrac', 'ActiveRespFrac', - 'FtimeAdj', 'FtimeQBOpt', 'PBopt', - 'NoIntegrate', 'HandleSelf', 'ScrambleSelf')){ + if ( + parameter %in% + c( + 'B_BaseRef', + 'MzeroMort', + 'UnassimRespFrac', + 'ActiveRespFrac', + 'FtimeAdj', + 'FtimeQBOpt', + 'PBopt', + 'NoIntegrate', + 'HandleSelf', + 'ScrambleSelf' + ) + ) { Rsim.scenario$params[[param.num]][groupnum + 1] <- value } - - if(parameter %in% c('QQ', 'DD', 'VV', 'HandleSwitch', 'PredPredWeight', - 'PreyPreyWeight')){ - linknum <- which(Rsim.scenario$params$PreyFrom %in% groupnum & - Rsim.scenario$params$PreyTo == groupnumto) + + if ( + parameter %in% + c('QQ', 'DD', 'VV', 'HandleSwitch', 'PredPredWeight', 'PreyPreyWeight') + ) { + linknum <- which( + Rsim.scenario$params$PreyFrom %in% + groupnum & + Rsim.scenario$params$PreyTo == groupnumto + ) Rsim.scenario$params[[param.num]][linknum] <- value } return(Rsim.scenario) @@ -207,7 +245,7 @@ adjust.scenario <- function(Rsim.scenario, parameter, group, groupto = NA, value #' #'@param bymonth Boolean value that denotes whether to use sim.year/sim.month combo #' or just sim.month as a sequential vector starting at 1. -#' +#' #'@return Returns an Rsim.scenario object with the new parameter. #' #' @@ -217,58 +255,78 @@ adjust.scenario <- function(Rsim.scenario, parameter, group, groupto = NA, value #' # Create a 50 yr Rsim scenario #' Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) #' # Adjust the ForcedPrey parameter for cod in year 1 for all months. Change the value to 10 -#' Rsim.scenario.adjusted <- adjust.forcing(Rsim.scenario, parameter="ForcedPrey",group = "cod", sim.year = 1, sim.month=0,value=10) +#' Rsim.scenario.adjusted <- adjust.forcing(Rsim.scenario, parameter="ForcedPrey",group = "cod", sim.year = 1, sim.month=0,value=10) #' head(Rsim.scenario.adjusted$forcing$ForcedPrey) #' #' #' -#'@export -adjust.forcing <- function(Rsim.scenario, parameter, group, sim.year = 1, sim.month = 0, - bymonth = F, value){ +#'@export +adjust.forcing <- function( + Rsim.scenario, + parameter, + group, + sim.year = 1, + sim.month = 0, + bymonth = F, + value +) { #Check that parameter and group exist - if(!parameter %in% c('ForcedPrey', 'ForcedMort', 'ForcedRecs', 'ForcedSearch', 'ForcedActresp', - 'ForcedMigrate', 'ForcedBio')){stop("Forcing parameter not found")} - if(!all(group %in% Rsim.scenario$params$spname)){ - stop("Groups not found:",group[!(group %in% Rsim.scenario$params$spname)]) + if ( + !parameter %in% + c( + 'ForcedPrey', + 'ForcedMort', + 'ForcedRecs', + 'ForcedSearch', + 'ForcedActresp', + 'ForcedMigrate', + 'ForcedBio' + ) + ) { + stop("Forcing parameter not found") + } + if (!all(group %in% Rsim.scenario$params$spname)) { + stop("Groups not found:", group[!(group %in% Rsim.scenario$params$spname)]) } #if(!group %in% Rsim.scenario$params$spname){stop("Group not found")} - - if(bymonth){ + + if (bymonth) { Rsim.scenario$forcing[[parameter]][sim.month, group] <- value - }else { - - #Create index in case the number of values is equal to length(sim.year) * length(sim.month) - ivalue <- 0 - - #Loop over years if more than 1 sim.year provided - for(iyear in seq_along(sim.year)){ - for(imonth in seq_along(sim.month)){ - ivalue <- ivalue + 1 - #look-up what rows correspond to the year - #Regex used to account for year.month row names in Effort Matrix - year.row <- which(gsub("\\..*", "", rownames(Rsim.scenario$forcing[[parameter]])) - == sim.year[iyear]) - - #identify what rows correspond to the sim.months - 0 indicates the whole year - if(sim.month[1] != 0){ - year.row <- year.row[1:12 %in% sim.month[imonth]] - } - - #Apply the value to the correct row - if(length(value) == 1){ - #If only 1 value is supplied for multiple years need to only point to that value - Rsim.scenario$forcing[[parameter]][year.row, group] <- value - }else if(sim.month[1] == 0){ - Rsim.scenario$forcing[[parameter]][year.row, group] <- value[iyear] - }else if(length(value) == length(sim.month)){ - Rsim.scenario$forcing[[parameter]][year.row, group] <- value[imonth] - }else { - Rsim.scenario$forcing[[parameter]][year.row, group] <- value[ivalue] + } else { + #Create index in case the number of values is equal to length(sim.year) * length(sim.month) + ivalue <- 0 + + #Loop over years if more than 1 sim.year provided + for (iyear in seq_along(sim.year)) { + for (imonth in seq_along(sim.month)) { + ivalue <- ivalue + 1 + #look-up what rows correspond to the year + #Regex used to account for year.month row names in Effort Matrix + year.row <- which( + gsub("\\..*", "", rownames(Rsim.scenario$forcing[[parameter]])) == + sim.year[iyear] + ) + + #identify what rows correspond to the sim.months - 0 indicates the whole year + if (sim.month[1] != 0) { + year.row <- year.row[1:12 %in% sim.month[imonth]] + } + + #Apply the value to the correct row + if (length(value) == 1) { + #If only 1 value is supplied for multiple years need to only point to that value + Rsim.scenario$forcing[[parameter]][year.row, group] <- value + } else if (sim.month[1] == 0) { + Rsim.scenario$forcing[[parameter]][year.row, group] <- value[iyear] + } else if (length(value) == length(sim.month)) { + Rsim.scenario$forcing[[parameter]][year.row, group] <- value[imonth] + } else { + Rsim.scenario$forcing[[parameter]][year.row, group] <- value[ivalue] + } } } } - } - + return(Rsim.scenario) } @@ -278,7 +336,7 @@ adjust.forcing <- function(Rsim.scenario, parameter, group, sim.year = 1, sim.mo #'Set Rsim.scenario parameters #' -#'Modifies the various parameters of the \code{rsim.scenario()} object. Parameters +#'Modifies the various parameters of the \code{rsim.scenario()} object. Parameters #'that can be adjusted using this function are: \code{params},\code{start_state}, #'\code{forcing},\code{fishing},\code{stanzas} #' @@ -287,13 +345,13 @@ adjust.forcing <- function(Rsim.scenario, parameter, group, sim.year = 1, sim.mo #'@inheritParams rsim.run #'@inheritParams rsim.fishing #'@param start_state Rsim starting values object generated by \code{rsim.state()} -#'@param forcing Rsim forcing matrix object generated by \code{rsim.forcing()} +#'@param forcing Rsim forcing matrix object generated by \code{rsim.forcing()} #'@param fishing Rsim fishing matrix object generated by \code{rsim.fishing()} #'@param stanzas Rsim stanza parameters object generated by \code{rsim.stanzas()} #' #'@return Returns an \code{Rsim.scenario} object with the new parameter. #' -#'@examples +#'@examples #' # Read in Rpath parameter file and generate balanced model #' Rpath <- rpath(AB.params) #' # Create a 50 yr Rsim scenario @@ -305,24 +363,51 @@ adjust.forcing <- function(Rsim.scenario, parameter, group, sim.year = 1, sim.mo #' #' #'@export -#' -set.rsim.scene<-function(Rsim.scenario,params=NULL,start_state=NULL,forcing=NULL,fishing=NULL,stanzas=NULL){ +#' +set.rsim.scene <- function( + Rsim.scenario, + params = NULL, + start_state = NULL, + forcing = NULL, + fishing = NULL, + stanzas = NULL +) { rsim <- list() class(rsim) <- 'Rsim.scenario' attr(rsim, 'eco.name') <- attr(Rsim.scenario, 'eco.name') # can add type checks later - if (!is.null(params)) {rsim$params <- params } else {rsim$params <- Rsim.scenario$params } - if (!is.null(start_state)){rsim$start_state <- start_state } else {rsim$start_state <- Rsim.scenario$start_state } - if (!is.null(forcing)) {rsim$forcing <- forcing } else {rsim$forcing <- Rsim.scenario$forcing } - if (!is.null(fishing)) {rsim$fishing <- fishing } else {rsim$fishing <- Rsim.scenario$fishing } - if (!is.null(stanzas)) {rsim$stanzas <- stanzas } else {rsim$stanzas <- Rsim.scenario$stanzas } + if (!is.null(params)) { + rsim$params <- params + } else { + rsim$params <- Rsim.scenario$params + } + if (!is.null(start_state)) { + rsim$start_state <- start_state + } else { + rsim$start_state <- Rsim.scenario$start_state + } + if (!is.null(forcing)) { + rsim$forcing <- forcing + } else { + rsim$forcing <- Rsim.scenario$forcing + } + if (!is.null(fishing)) { + rsim$fishing <- fishing + } else { + rsim$fishing <- Rsim.scenario$fishing + } + if (!is.null(stanzas)) { + rsim$stanzas <- stanzas + } else { + rsim$stanzas <- Rsim.scenario$stanzas + } return(rsim) } #'Retrieve parameters from an Rsim scenario #' -#'Helper function that will retrieve the parameters that were used -#'in an Rsim scenario +#'Helper function that will retrieve the parameters that were used +#'in an Rsim scenario #' #'@family Get functions #' @@ -337,20 +422,20 @@ set.rsim.scene<-function(Rsim.scenario,params=NULL,start_state=NULL,forcing=NULL #' Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) #' params <- get.rsim.params(Rsim.scenario) #' names(params) -#' +#' #' #' #'@export -#' -get.rsim.params<-function(Rsim.scenario){ +#' +get.rsim.params <- function(Rsim.scenario) { return(Rsim.scenario$params) } #'Retrieve starting state values from an Rsim scenario #' -#'Helper function that will retrieve the starting state values that were used -#'in an Rsim scenario +#'Helper function that will retrieve the starting state values that were used +#'in an Rsim scenario #' #'@family Get functions #' @@ -364,20 +449,20 @@ get.rsim.params<-function(Rsim.scenario){ #' # Create a 50 yr Rsim scenario #' Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) #' params <- get.rsim.start_state(Rsim.scenario) -#' names(params) +#' names(params) #' #' #' #'@export -#' -get.rsim.start_state<-function(Rsim.scenario){ +#' +get.rsim.start_state <- function(Rsim.scenario) { return(Rsim.scenario$start_state) } #'Retrieve forcing parameters from an Rsim scenario #' -#'Helper function that will retrieve the forcing parameters that were used in an -#'Rsim scenario +#'Helper function that will retrieve the forcing parameters that were used in an +#'Rsim scenario #' #'@family Get functions #' @@ -385,7 +470,7 @@ get.rsim.start_state<-function(Rsim.scenario){ #' #'@return Returns a `forcing` object. #' -#'@examples +#'@examples #' # Read in Rpath parameter file and generate balanced model #' Rpath <- rpath(AB.params) #' # Create a 50 yr Rsim scenario @@ -395,15 +480,15 @@ get.rsim.start_state<-function(Rsim.scenario){ #' #' #'@export -#' -get.rsim.forcing<-function(Rsim.scenario){ +#' +get.rsim.forcing <- function(Rsim.scenario) { return(Rsim.scenario$forcing) } #'Retrieve fishing forcing parameters from an Rsim scenario #' -#'Helper function that will retrieve the fishing forcing parameters that were used in an -#'Rsim scenario +#'Helper function that will retrieve the fishing forcing parameters that were used in an +#'Rsim scenario #' #'@family Get functions #' @@ -411,7 +496,7 @@ get.rsim.forcing<-function(Rsim.scenario){ #' #'@return Returns a `fishing` object. #' -#'@examples +#'@examples #' # Read in Rpath parameter file and generate balanced model #' Rpath <- rpath(AB.params) #' # Create a 50 yr Rsim scenario @@ -422,15 +507,15 @@ get.rsim.forcing<-function(Rsim.scenario){ #' #' #'@export -#' -get.rsim.fishing<-function(Rsim.scenario){ +#' +get.rsim.fishing <- function(Rsim.scenario) { return(Rsim.scenario$fishing) } #'Retrieve stanza parameters from an Rsim scenario #' -#'Helper function that will retrieve the stanza parameters that were used in an -#'Rsim scenario +#'Helper function that will retrieve the stanza parameters that were used in an +#'Rsim scenario #' #'@family Get functions #' @@ -438,7 +523,7 @@ get.rsim.fishing<-function(Rsim.scenario){ #' #'@return Returns a `stanzas` object. #' -#'@examples +#'@examples #' # Read in Rpath parameter file and generate balanced model #' Rpath <- rpath(AB.params) #' # Create a 50 yr Rsim scenario @@ -449,14 +534,7 @@ get.rsim.fishing<-function(Rsim.scenario){ #' #' #'@export -#' -get.rsim.stanzas<-function(Rsim.scenario){ +#' +get.rsim.stanzas <- function(Rsim.scenario) { return(Rsim.scenario$stanzas) } - - - - - - - diff --git a/R/Rpath-package.R b/R/Rpath-package.R index 20d4d66b..775c3d7f 100644 --- a/R/Rpath-package.R +++ b/R/Rpath-package.R @@ -1,12 +1,13 @@ -#' Rpath: A package implementing mass balance algorithms designed to work with +#' Rpath: A package implementing mass balance algorithms designed to work with #' fisheries data sources. #' #' The Rpath package provides two categories of important functions: #' rpath and rsim. -#' +#' +#' @keywords internal #' @section Rpath functions: #' The Rpath functions generate the balanced snap shot of energy flow through the #' system. This is done by solving a series of linear equations for unknown biomass #' or ecotrophic efficiency. #' -"_PACKAGE" \ No newline at end of file +"_PACKAGE" diff --git a/R/Rpath_pkg.R b/R/Rpath_pkg.R deleted file mode 100644 index 5bb2f078..00000000 --- a/R/Rpath_pkg.R +++ /dev/null @@ -1,12 +0,0 @@ -#' Rpath: A package implementing mass balance algorithms designed to work with -#' fisheries data sources. -#' -#' The Rpath package provides two categories of important functions: -#' rpath and rsim. -#' -#' @section rpath functions: -#' The rpath functions generate the balanced snap shot of energy flow through the -#' system. This is done by solving a series of linear equations for unknown biomass -#' or ecotrophic efficiency. -#' -"_PACKAGE" \ No newline at end of file diff --git a/R/ecopath.R b/R/ecopath.R index 603906ea..b8cc9564 100644 --- a/R/ecopath.R +++ b/R/ecopath.R @@ -1,6 +1,6 @@ #'Ecopath module of Rpath #' -#'Performs initial mass balance using a \code{\link{Rpath.params}()} file +#'Performs initial mass balance using a \code{Rpath.params()} file #' #'@family Rpath functions #' @@ -18,166 +18,210 @@ #'@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 - + # Model Parameters - Basic parameters, detritus fate, catch, discards in that order model <- copy(Rpath.params$model) - + #Diet Parameters - diet matrix, predators as columns, prey as rows - include #producers as predators even though they do not consume any groups diet <- copy(Rpath.params$diet) - + #Check that all columns of model are numeric and not logical - if(length(which(sapply(model, class) == 'logical')) > 0){ + if (length(which(sapply(model, class) == 'logical')) > 0) { logic.col <- which(sapply(model, class) == 'logical') - for(i in 1:length(logic.col)){ + for (i in 1:length(logic.col)) { set(model, j = logic.col[i], value = as.numeric(model[[logic.col[i]]])) } } - + #Remove first column if names (factor or character) - if(sapply(diet, class)[1] == 'factor') diet[, 1 := NULL] - if(sapply(diet, class)[1] == 'character') diet[, 1 := NULL] + if (sapply(diet, class)[1] == 'factor') { + diet[, 1 := NULL] + } + if (sapply(diet, class)[1] == 'character') { + diet[, 1 := NULL] + } #Adjust diet comp of mixotrophs mixotrophs <- which(model[, Type] > 0 & model[, Type] < 1) mix.Q <- 1 - model[mixotrophs, Type] - for(i in seq_along(mixotrophs)){ + for (i in seq_along(mixotrophs)) { new.dc <- diet[, mixotrophs[i], with = F] * mix.Q[i] diet[, mixotrophs[i] := new.dc] } - + #Convert NAs to zero in diet matrix diet[is.na(diet)] <- 0 - + # Get number of groups, living, dead, and gear ngroups <- nrow(model) - nliving <- nrow(model[Type < 2, ]) - ndead <- nrow(model[Type == 2, ]) - ngear <- nrow(model[Type == 3, ]) + nliving <- nrow(model[Type < 2, ]) + ndead <- nrow(model[Type == 2, ]) + ngear <- nrow(model[Type == 3, ]) nodetrdiet <- diet[1:nliving, ] model[is.na(DetInput), DetInput := 0] # fill in GE(PQ), QB, or PB from other inputs - GE <- ifelse(is.na(model[, ProdCons]), model[, PB / QB], model[, ProdCons]) - QB.1 <- ifelse(is.na(model[, QB]), model[, PB / GE], model[, QB]) - PB.1 <- ifelse(is.na(model[, PB]), model[, ProdCons * QB], model[, PB]) + # KYA Aug 2025 - changed this logic so PC would be recalculated if PB and QB supplied + #GE <- ifelse(is.na(model[, ProdCons]), model[, PB / QB], model[, ProdCons]) + GE <- ifelse( + !is.na(model[, QB]) & !is.na(model[, PB]), + model[, PB / QB], + model[, ProdCons] + ) + QB.1 <- ifelse(is.na(model[, QB]), model[, PB / GE], model[, QB]) + PB.1 <- ifelse(is.na(model[, PB]), model[, ProdCons * QB], model[, PB]) model[, QB := QB.1] model[, PB := PB.1] - + # define landings, discards, necessary sums - landmat <- model[, (10 + ndead + 1):(10 + ndead + ngear), with = F] - discardmat <- model[, (10 + ndead + 1 + ngear):(10 + ndead + (2 * ngear)), with = F] + landmat <- model[, (10 + ndead + 1):(10 + ndead + ngear), with = F] + discardmat <- model[, + (10 + ndead + 1 + ngear):(10 + ndead + (2 * ngear)), + with = F + ] totcatchmat <- landmat + discardmat - - if (is.data.frame(totcatchmat)){ + + if (is.data.frame(totcatchmat)) { totcatch <- rowSums(totcatchmat) - landings <- rowSums(landmat) - discards <- rowSums(discardmat) - gearland <- colSums(landmat, na.rm = T) + landings <- rowSums(landmat) + discards <- rowSums(discardmat) + gearland <- colSums(landmat, na.rm = T) geardisc <- colSums(discardmat, na.rm = T) - }else{ + } else { totcatch <- totcatchmat - landings <- landmat - discards <- discardmat - gearland <- sum(landmat, na.rm = T) - geardisc <- sum(discardmat, na.rm = T) - } - + landings <- landmat + discards <- discardmat + gearland <- sum(landmat, na.rm = T) + geardisc <- sum(discardmat, na.rm = T) + } + geartot <- gearland + geardisc model[, landings := landings] model[, discards := discards] model[, totcatch := totcatch] # flag missing pars and subset for estimation - model[, noB := 0] - model[, noEE := 0] + model[, noB := 0] + model[, noEE := 0] model[, alive := 0] - model[, BEE := 0] - model[is.na(Biomass), noB := 1] - model[is.na(EE), noEE := 1] - model[Type < 2, alive := 1] + 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[BEE == 1 & is.na(PB), noPB := 1] + + if (any(model$Type == 0 & is.na(model$QB) & is.na(model$ProdCons))) { + stop( + "A consumer is missing both QB and ProdCons - balance failed. Use check.rpath.params() to diagnose." + ) + } # define detritus fate matrix detfate <- model[, (10 + 1):(10 + ndead), with = F] + detdetfate <- model[Type == 2, (10 + 1):(10 + ndead), with = F] # set up and solve the system of equations for living group B or EE - living <- model[alive == 1, ] - + living <- model[alive == 1, ] + #Set up right hand side b living[, Ex := totcatch + BioAcc] living[, BioQB := Biomass * QB] - cons <- as.matrix(nodetrdiet) * living$BioQB[col(as.matrix(nodetrdiet))] - living[, b := Ex + rowSums(cons, na.rm = T)] - + cons <- as.matrix(nodetrdiet) * living$BioQB[col(as.matrix(nodetrdiet))] + living[, b := Ex + rowSums(cons, na.rm = T)] + #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)] #living[BEE == 1, diag.a := 0] #Need to work on this solution - - A <- matrix(0, nliving, nliving) + + A <- matrix(0, nliving, nliving) diag(A) <- living[, diag.a] - QBDC <- as.matrix(nodetrdiet) * living$QB[col(as.matrix(nodetrdiet))] + QBDC <- as.matrix(nodetrdiet) * living$QB[col(as.matrix(nodetrdiet))] dimnames(QBDC) <- list(NULL, NULL) QBDC[is.na(QBDC)] <- 0 #Flip noB flag for known B and EE #living[BEE == 1, noB := 1] QBDCa <- as.matrix(QBDC) * living$noB[col(as.matrix(QBDC))] - A <- A - QBDCa + A <- A - QBDCa #Switch flag back #living[BEE == 1, noB := 0] - + + # Check for any missing info that will prevent solving + if (any(is.na(A))) { + stop( + "Model is missing parameters - can't be balanced. Use check.rpath.params() to diagnose." + ) + } + # Generalized inverse does the actual solving #Invert A and multiple by b to get x (unknowns) x <- MASS::ginv(A, tol = .Machine$double.eps) %*% living[, b] - + #Assign unknown values living[, EEa := x * noEE] living[is.na(EE), EE := EEa] - + 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] living[is.na(QBloss), QBloss := 0] - loss <- c((living[, M0] * living[, Biomass]) + - (living[, Biomass] * living[, QBloss] * living[, Unassim]), - model[Type ==2, DetInput], - geardisc) - detinputs <- colSums(loss * detfate) - detdiet <- diet[(nliving + 1):(nliving + ndead), ] - BQB <- living[, Biomass * QB] - detcons <- as.matrix(detdiet) * BQB[col(as.matrix(detdiet))] + #KYA fix Aug 2025 + #loss <- c((living[, M0] * living[, Biomass]) + + # (living[, Biomass] * living[, QBloss] * living[, Unassim]), + # model[Type ==2, DetInput], + # geardisc) + #detinputs1 <- colSums(loss * detfate) + #detinputs1 is "first pass" at det inputs, final detinputs is after initial EE + loss <- c( + (living[, M0] * living[, Biomass]) + + (living[, Biomass] * living[, QBloss] * living[, Unassim]), + rep(0, ndead), + geardisc + ) + detinputs1 <- colSums(loss * detfate + model[, DetInput]) + ## end fix + detdiet <- diet[(nliving + 1):(nliving + ndead), ] + BQB <- living[, Biomass * QB] + detcons <- as.matrix(detdiet) * BQB[col(as.matrix(detdiet))] detoutputs <- rowSums(detcons, na.rm = T) - EE <- c(living[, EE], as.vector(detoutputs / detinputs)) + det_unused <- ifelse(detinputs1 > detoutputs, detinputs1 - detoutputs, 0.0) + detinputs <- detinputs1 + colSums(det_unused * detdetfate) + EE <- c(living[, EE], as.vector(detoutputs / detinputs)) # added by kya - # if a detritus biomass is put into the spreadsheet, use that and - # calculate PB. If no biomass, but a PB, use that pb with inflow to - # calculate biomass. If neither, use default PB=0.5, Bio = inflow/PB + # if a detritus biomass is put into the spreadsheet, use that and + # calculate PB. If no biomass, but a PB, use that pb with inflow to + # calculate biomass. If neither, use default PB=0.5, Bio = inflow/PB # This is done because Ecosim requires a detrital biomass. - - Default_Detrital_PB <- 0.5 - inDetPB <- model[(nliving + 1):(nliving + ndead), PB] - 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 + + Default_Detrital_PB <- 0.5 + inDetPB <- model[(nliving + 1):(nliving + ndead), PB] + 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 <- as.numeric(detinputs) / DetB # Trophic Level calcs - b <- rep(1, ngroups) - TLcoeff <- matrix(0, ngroups, ngroups) + b <- rep(1, ngroups) + TLcoeff <- matrix(0, ngroups, ngroups) diag(TLcoeff) <- rep(1, ngroups) - gearcons <- as.matrix(totcatchmat) / geartot[col(as.matrix(totcatchmat))] + gearcons <- as.matrix(totcatchmat) / geartot[col(as.matrix(totcatchmat))] dimnames(gearcons) <- list(NULL, NULL) gearcons[is.na(gearcons)] <- 0 dietplus <- as.matrix(diet) @@ -192,7 +236,7 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) { # } #Adjust for diet import (Consumption outside model) import <- which(dietplus[nrow(diet), ] > 0) - for(i in seq_along(import)){ + for (i in seq_along(import)) { import.denom <- 1 - dietplus[nrow(diet), import[i]] dietplus[, import[i]] <- dietplus[, import[i]] / import.denom } @@ -200,108 +244,136 @@ rpath <- function(Rpath.params, eco.name = NA, eco.area = 1) { dietplus <- rbind(dietplus, matrix(0, ngear, nliving)) dietplus <- cbind(dietplus, matrix(0, ngroups, ndead), gearcons) TLcoeffA <- TLcoeff - dietplus - TL <- solve(t(TLcoeffA), b) + TL <- solve(t(TLcoeffA), b) #kya changed these following four lines for detritus, and removing NAs #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 + Bplus <- c(living[, Biomass], DetB, rep(0.0, ngear)) + + #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) & PBplus > 0.0 & !(is.na(GE) | is.nan(GE) | is.infinite(GE)) + ] <- + (PBplus / GE)[ + is.na(QBplus) & PBplus > 0.0 & !(is.na(GE) | is.nan(GE) | is.infinite(GE)) + ] 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 - - balanced <- list(Group = model[, Group], - TL = TL, - Biomass = Bplus, - PB = PBplus, - QB = QBplus, - EE = EEplus, - GE = GE, - Removals = RemPlus) - - M0plus <- c(living[, M0], as.vector(detoutputs / detinputs)) - gearF <- as.matrix(totcatchmat) / living[, Biomass][row(as.matrix(totcatchmat))] - newcons <- as.matrix(nodetrdiet) * BQB[col(as.matrix(nodetrdiet))] - predM <- as.matrix(newcons) / living[, Biomass][row(as.matrix(newcons))] - predM <- rbind(predM, detcons) - morts <- list(Group = model[Type < 3, Group], - PB = model[Type < 3, PB], - M0 = M0plus, - F = gearF[1:(nliving + ndead), ], - M2 = predM) - + + balanced <- list( + Group = model[, Group], + TL = TL, + Biomass = Bplus, + PB = PBplus, + QB = QBplus, + EE = EEplus, + GE = GE, + Removals = RemPlus + ) + + M0plus <- c(living[, M0], as.vector(detoutputs / detinputs)) + gearF <- as.matrix(totcatchmat) / + living[, Biomass][row(as.matrix(totcatchmat))] + newcons <- as.matrix(nodetrdiet) * BQB[col(as.matrix(nodetrdiet))] + predM <- as.matrix(newcons) / living[, Biomass][row(as.matrix(newcons))] + predM <- rbind(predM, detcons) + morts <- list( + Group = model[Type < 3, Group], + PB = model[Type < 3, PB], + M0 = M0plus, + F = gearF[1:(nliving + ndead), ], + M2 = predM + ) + # convert from levels to characters gnames <- as.character(balanced$Group) - + # cleanup before sending to sim -- C code wants 0 as missing value, not NA balanced$Biomass[is.na(balanced$Biomass)] <- 0 - balanced$PB[is.na(balanced$PB)] <- 0 - balanced$QB[is.na(balanced$QB)] <- 0 - balanced$EE[is.na(balanced$EE)] <- 0 - balanced$GE[is.na(balanced$GE)] <- 0 - model$BioAcc[is.na(model$BioAcc)] <- 0 + balanced$PB[is.na(balanced$PB)] <- 0 + balanced$QB[is.na(balanced$QB)] <- 0 + balanced$EE[is.na(balanced$EE)] <- 0 + balanced$GE[is.na(balanced$GE)] <- 0 + model$BioAcc[is.na(model$BioAcc)] <- 0 model$Unassim[is.na(model$Unassim)] <- 0 - dietm <- as.matrix(diet) - dimnames(dietm) <- list(c(gnames[1:(nliving+ndead)],"Import"), gnames[1:nliving]) - dietm[is.na(dietm)] <- 0 - landmatm <- as.matrix(landmat) - dimnames(landmatm) <- list(gnames, gnames[(ngroups-ngear+1):ngroups]) - landmatm[is.na(landmatm)] <- 0 - discardmatm <- as.matrix(discardmat) - dimnames(discardmatm) <- list(gnames, gnames[(ngroups-ngear+1):ngroups]) - discardmatm[is.na(discardmatm)] <- 0 - detfatem <- as.matrix(detfate) - dimnames(detfatem) <- list(gnames, gnames[(nliving+1):(nliving+ndead)]) - detfatem[is.na(detfatem)] <- 0 + dietm <- as.matrix(diet) + dimnames(dietm) <- list( + c(gnames[1:(nliving + ndead)], "Import"), + gnames[1:nliving] + ) + dietm[is.na(dietm)] <- 0 + landmatm <- as.matrix(landmat) + dimnames(landmatm) <- list(gnames, gnames[(ngroups - ngear + 1):ngroups]) + landmatm[is.na(landmatm)] <- 0 + discardmatm <- as.matrix(discardmat) + dimnames(discardmatm) <- list(gnames, gnames[(ngroups - ngear + 1):ngroups]) + discardmatm[is.na(discardmatm)] <- 0 + detfatem <- as.matrix(detfate) + dimnames(detfatem) <- list(gnames, gnames[(nliving + 1):(nliving + ndead)]) + detfatem[is.na(detfatem)] <- 0 # Add names for output list - out.Group <- gnames; names(out.Group) <- gnames - out.type <- model[, Type]; names(out.type) <- gnames - out.TL <- TL; names(out.TL) <- gnames - out.Biomass <- balanced$Biomass; names(out.Biomass) <- gnames - out.PB <- balanced$PB; names(out.PB) <- gnames - out.QB <- balanced$QB; names(out.QB) <- gnames - out.EE <- balanced$EE; names(out.EE) <- gnames - out.BA <- model[, BioAcc]; names(out.BA) <- gnames - out.Unassim <- model[, Unassim]; names(out.Unassim) <- gnames - out.GE <- balanced$GE; names(out.GE) <- gnames - + out.Group <- gnames + names(out.Group) <- gnames + out.type <- model[, Type] + names(out.type) <- gnames + out.TL <- TL + names(out.TL) <- gnames + out.Biomass <- balanced$Biomass + names(out.Biomass) <- gnames + out.PB <- balanced$PB + names(out.PB) <- gnames + out.QB <- balanced$QB + names(out.QB) <- gnames + out.EE <- balanced$EE + names(out.EE) <- gnames + out.BA <- model[, BioAcc] + names(out.BA) <- gnames + out.Unassim <- model[, Unassim] + names(out.Unassim) <- gnames + out.GE <- balanced$GE + names(out.GE) <- gnames + # list structure for sim inputs - path.model <- list(NUM_GROUPS = ngroups, - NUM_LIVING = nliving, - NUM_DEAD = ndead, - NUM_GEARS = ngear, - Group = out.Group, - type = out.type, - TL = out.TL, - Biomass = out.Biomass, - PB = out.PB, - QB = out.QB, - EE = out.EE, - BA = out.BA, - Unassim = out.Unassim, - GE = out.GE, - DC = dietm, - DetFate = detfatem, - Landings = landmatm, - Discards = discardmatm) - -#Define class of output -class(path.model) <- 'Rpath' -attr(path.model, 'eco.name') <- eco.name -attr(path.model, 'eco.area') <- eco.area - -return(path.model) + path.model <- list( + NUM_GROUPS = ngroups, + NUM_LIVING = nliving, + NUM_DEAD = ndead, + NUM_GEARS = ngear, + Group = out.Group, + type = out.type, + TL = out.TL, + Biomass = out.Biomass, + PB = out.PB, + QB = out.QB, + EE = out.EE, + BA = out.BA, + Unassim = out.Unassim, + GE = out.GE, + DC = dietm, + DetFate = detfatem, + Landings = landmatm, + Discards = discardmatm + ) + + #Define class of output + class(path.model) <- 'Rpath' + attr(path.model, 'eco.name') <- eco.name + attr(path.model, 'eco.area') <- eco.area + + return(path.model) } @@ -314,32 +386,37 @@ return(path.model) #' #'@inheritParams rpath #' -#'@return Calculates and adds biomass and consumption for trailing stanza groups. +#'@return Calculates and adds biomass and consumption for trailing stanza groups. #' Also adds weight at age and number at age for multi-staza groups. -#' +#' #'@import data.table -#'@export -rpath.stanzas <- function(Rpath.params){ +#'@export +rpath.stanzas <- function(Rpath.params) { #Need to define variables to eliminate check() note about no visible binding StGroupNum <- First <- StanzaNum <- VBGF_d <- VBGF_Ksp <- Last <- GroupNum <- NULL WageS <- age <- QageS <- Survive <- Z <- survive_L <- bs.num <- qs.num <- Leading <- Oldest <- NULL Group <- Biomass <- R <- NageS <- bs.denom <- bs <- qs.denom <- qs <- Cons <- NULL QB <- BAB <- Ex <- NULL + # Added Aug 2025 - if no stanzas, silently return original (prob no warning needed?) + if (Rpath.params$stanza$NStanzaGroups == 0) { + return(Rpath.params) + } + #Determine the total number of groups with multistanzas - Nsplit <- Rpath.params$stanza$NStanzaGroups - groupfile <- Rpath.params$stanza$stgroups + Nsplit <- Rpath.params$stanza$NStanzaGroups + groupfile <- Rpath.params$stanza$stgroups stanzafile <- Rpath.params$stanza$stindiv - + # Add Oldest column so age groups will be properly calculated stanzafile[Last == max(Last), Oldest := T, by = StGroupNum] stanzafile[is.na(Oldest), Oldest := F] - + #Need to add vector of stanza number - lastmonth <- rep(NA,Nsplit) - for(isp in 1:Nsplit){ + lastmonth <- rep(NA, Nsplit) + for (isp in 1:Nsplit) { #Put the stanzas in order for each split species - stnum <- order(stanzafile[StGroupNum == isp, First]) + stnum <- order(stanzafile[StGroupNum == isp, First]) stanzafile[StGroupNum == isp, StanzaNum := stnum] #Calculate the last month for the final ("leading") stanza @@ -348,129 +425,147 @@ rpath.stanzas <- function(Rpath.params){ # with rapid growth but low mortality (e.g. marine mammals). # So instead, calculate biomass out for a very long time and # taking 0.99999 of cumulative biomass as a cutoff. - + #this selects all of the stanza lines, then picks the last one #(maybe data table has a better way...) - + stmax <- max(stanzafile[StGroupNum == isp, StanzaNum]) - st <- stanzafile[StGroupNum == isp & StanzaNum==stmax,] + st <- stanzafile[StGroupNum == isp & StanzaNum == stmax, ] - gp <- groupfile[isp,] + gp <- groupfile[isp, ] #Max age class in months should be one less than a multiple of 12 #(trying 5999 - probably overkill but for safety) AGE <- st$First:5999 - mz <- (st$Z + gp$BAB)/12 - k <- gp$VBGF_Ksp - d <- gp$VBGF_d - NN <- shift(cumprod(rep(exp(-1*mz),length(AGE))),1,1.0) - BB <- NN * (1 - exp(-k * (1 - d) * (AGE))) ^ (1 / (1 - d)) - BBcum <- cumsum(BB)/sum(BB) + mz <- (st$Z + gp$BAB) / 12 + k <- gp$VBGF_Ksp + d <- gp$VBGF_d + NN <- shift(cumprod(rep(exp(-1 * mz), length(AGE))), 1, 1.0) + BB <- NN * (1 - exp(-k * (1 - d) * (AGE)))^(1 / (1 - d)) + BBcum <- cumsum(BB) / sum(BB) #Age at which 0.99999 of cumulative leading stanza biomass is represented, #rounded to nearest higher multiple of 12 (-1 since index starts at 0) - lastmonth[isp] <- ceiling(AGE[min(which(BBcum>0.99999))]/12) * 12 - 1 + lastmonth[isp] <- ceiling(AGE[min(which(BBcum > 0.99999))] / 12) * 12 - 1 } - + #Save the maximum month vector in the table groupfile[, last := lastmonth] - - for(isp in 1:Nsplit){ + + for (isp in 1:Nsplit) { nstanzas <- groupfile[ - StGroupNum == isp, nstanzas] + StGroupNum == isp, + nstanzas + ] stanzafile[StGroupNum == isp & Oldest, Last := lastmonth[isp]] - + #Grab ecopath group codes group.codes <- stanzafile[StGroupNum == isp, GroupNum] - + #Grab index for first and last months for stanzas - first <- stanzafile[StGroupNum == isp, First] - second <- stanzafile[StGroupNum == isp, Last] + first <- stanzafile[StGroupNum == isp, First] + second <- stanzafile[StGroupNum == isp, Last] - #Calculate weight and consumption at age - StGroup <- data.table(age = stanzafile[StGroupNum == isp & StanzaNum == 1, First]: - lastmonth[isp]) + #Calculate weight and consumption at age + StGroup <- data.table( + age = stanzafile[StGroupNum == isp & StanzaNum == 1, First]:lastmonth[isp] + ) #Calculate monthly generalized k: (Ksp * 3) / 12 k <- (groupfile[StGroupNum == isp, VBGF_Ksp] * 3) / 12 - d <- groupfile[StGroupNum == isp, VBGF_d] - StGroup[, WageS := (1 - exp(-k * (1 - d) * (age))) ^ (1 / (1 - d))] - StGroup[, QageS := WageS ^ d] - + d <- groupfile[StGroupNum == isp, VBGF_d] + StGroup[, WageS := (1 - exp(-k * (1 - d) * (age)))^(1 / (1 - d))] + StGroup[, QageS := WageS^d] + #Calculate the relative number of animals at age a #Vector of survival rates from 1 stanza to the next - + #Unwind the by-stanza mortality rates into by-month survival rates #by looping through the stanzas - survive_L <- rep(NA,length(StGroup$age)) - for(ist in 1:nstanzas){ + survive_L <- rep(NA, length(StGroup$age)) + for (ist in 1:nstanzas) { #Convert Z to a monthly Z - month.z <- (stanzafile[StGroupNum == isp & StanzaNum == ist, Z] + - groupfile[StGroupNum == isp, BAB]) / 12 - survive_L[which(StGroup$age %in% first[ist]:second[ist])] <- exp(-1*month.z) - } - + month.z <- (stanzafile[StGroupNum == isp & StanzaNum == ist, Z] + + groupfile[StGroupNum == isp, BAB]) / + 12 + survive_L[which(StGroup$age %in% first[ist]:second[ist])] <- exp( + -1 * month.z + ) + } + #Shift survival rates forward one - first survival is 1.0 - survive_L <- shift(survive_L,1,1.0) + survive_L <- shift(survive_L, 1, 1.0) #Use cumulative product to get overall survival to each month StGroup[, Survive := cumprod(survive_L)] - + StGroup[, B := Survive * WageS] StGroup[, Q := Survive * QageS] - for(ist in 1:nstanzas){ + for (ist in 1:nstanzas) { #Numerator for the relative biomass/consumption calculations b.num <- StGroup[age %in% first[ist]:second[ist], sum(B)] q.num <- StGroup[age %in% first[ist]:second[ist], sum(Q)] - + stanzafile[StGroupNum == isp & StanzaNum == ist, bs.num := b.num] stanzafile[StGroupNum == isp & StanzaNum == ist, qs.num := q.num] } - + #Scale numbers up to total recruits BaseStanza <- stanzafile[StGroupNum == isp & Leading == T, ] - BioPerEgg <- StGroup[age %in% BaseStanza[, First]:BaseStanza[, Last], sum(B)] - recruits <- Rpath.params$model[Group == BaseStanza[, Group], Biomass] / BioPerEgg + BioPerEgg <- StGroup[ + age %in% BaseStanza[, First]:BaseStanza[, Last], + sum(B) + ] + recruits <- Rpath.params$model[Group == BaseStanza[, Group], Biomass] / + BioPerEgg #Save recruits groupfile[StGroupNum == isp, R := recruits] #Numbers at age S StGroup[, NageS := Survive * recruits] - + #Calculate relative biomass stanzafile[StGroupNum == isp, bs.denom := sum(bs.num)] stanzafile[StGroupNum == isp, bs := bs.num / bs.denom] - + #Calculate relative consumption stanzafile[StGroupNum == isp, qs.denom := sum(qs.num)] stanzafile[StGroupNum == isp, qs := qs.num / qs.denom] - + #Use leading group to calculate other biomasses - stanzafile[StGroupNum == isp & Leading == T, - Biomass := Rpath.params$model[Group == BaseStanza[, Group], Biomass]] + stanzafile[ + StGroupNum == isp & Leading == T, + Biomass := Rpath.params$model[Group == BaseStanza[, Group], Biomass] + ] B <- stanzafile[StGroupNum == isp & Leading == T, Biomass / bs] stanzafile[StGroupNum == isp, Biomass := bs * B] - + #Use leading group to calculate other consumption - stanzafile[StGroupNum == isp & Leading == T, - Cons := Rpath.params$model[Group == BaseStanza[, Group], QB] * - Rpath.params$model[Group == BaseStanza[, Group], Biomass]] + stanzafile[ + StGroupNum == isp & Leading == T, + Cons := Rpath.params$model[Group == BaseStanza[, Group], QB] * + Rpath.params$model[Group == BaseStanza[, Group], Biomass] + ] Q <- stanzafile[StGroupNum == isp & Leading == T, Cons / qs] stanzafile[StGroupNum == isp, Cons := qs * Q] stanzafile[, QB := Cons / Biomass] - - Rpath.params$stanzas$StGroup[[isp]] <- StGroup - + + Rpath.params$stanzas$StGroup[[isp]] <- StGroup } - + #Drop extra columns - stanzafile[, c('bs.num', 'bs.denom', 'bs', 'qs.num', 'qs.denom', 'qs') := NULL] - + stanzafile[, + c('bs.num', 'bs.denom', 'bs', 'qs.num', 'qs.denom', 'qs') := NULL + ] + #Push biomass to modfile - for(i in 1:nrow(stanzafile)){ - Rpath.params$model[Group == stanzafile[i, Group], Biomass := stanzafile[i, Biomass]] + for (i in 1:nrow(stanzafile)) { + Rpath.params$model[ + Group == stanzafile[i, Group], + Biomass := stanzafile[i, Biomass] + ] } - + #Push consumption to modfile - for(i in 1:nrow(stanzafile)){ + for (i in 1:nrow(stanzafile)) { Rpath.params$model[Group == stanzafile[i, Group], QB := stanzafile[i, QB]] } - + return(Rpath.params) -} +} diff --git a/R/param.R b/R/param.R index 2d701b17..d722edd0 100644 --- a/R/param.R +++ b/R/param.R @@ -14,7 +14,7 @@ #' #'@return Outputs a list object of \code{Rpath.params} which are populated with values #' of NA or logical default values. Values can then be filled in using -#' R. Use \code{\link{check.rpath.params}()} to ensure parameter files are filled out +#' R. Use \code{check.rpath.params()} to ensure parameter files are filled out #' correctly (NOTE: This does not ensure data is correct just that it is #' in the right places). #'@import data.table @@ -22,46 +22,56 @@ create.rpath.params <- function(group, type, stgroup = NA) { #Need to define variables to eliminate check() note about no visible binding Group <- DetInput <- V1 <- StGroupNum <- NULL - + Rpath.params <- list() - - pred.group <- group[which(type < 2)] - prey.group <- group[which(type < 3)] - det.group <- group[which(type == 2)] + + pred.group <- group[which(type < 2)] + prey.group <- group[which(type < 3)] + det.group <- group[which(type == 2)] fleet.group <- group[which(type == 3)] - + #Model parameters model <- data.table( - Group = group, - Type = type, - Biomass = as.numeric(NA), - PB = as.numeric(NA), - QB = as.numeric(NA), - EE = as.numeric(NA), - ProdCons = as.numeric(NA), - BioAcc = as.numeric(NA), - Unassim = as.numeric(NA), - DetInput = as.numeric(NA) + Group = group, + Type = type, + Biomass = as.numeric(NA), + PB = as.numeric(NA), + QB = as.numeric(NA), + EE = as.numeric(NA), + ProdCons = as.numeric(NA), + BioAcc = as.numeric(NA), + Unassim = as.numeric(NA), + DetInput = as.numeric(NA) ) - + #Add detritial groups for (i in 1:length(det.group)) { model[Group %in% det.group, DetInput := 0] model[, V1 := as.numeric(NA)] setnames(model, "V1", det.group[i]) } - + #Add fleets twice - Landings and Discards for (i in 1:length(fleet.group)) { - model[, V1 := c(rep(0, length(group) - length(fleet.group)), rep(NA, length(fleet.group)))] + model[, + V1 := c( + rep(0, length(group) - length(fleet.group)), + rep(NA, length(fleet.group)) + ) + ] setnames(model, "V1", fleet.group[i]) } for (i in 1:length(fleet.group)) { - model[, V1 := c(rep(0, length(group) - length(fleet.group)), rep(NA, length(fleet.group)))] + model[, + V1 := c( + rep(0, length(group) - length(fleet.group)), + rep(NA, length(fleet.group)) + ) + ] setnames(model, "V1", paste(fleet.group[i], '.disc', sep = '')) } Rpath.params$model <- model - + #Diet matrix diet <- data.table(Group = c(prey.group, 'Import')) for (i in 1:length(pred.group)) { @@ -69,75 +79,81 @@ create.rpath.params <- function(group, type, stgroup = NA) { setnames(diet, "V1", pred.group[i]) } Rpath.params$diet <- diet - + #Multistanza parameters if (length(stgroup) > 1) { #Group Parameters - StanzaGroups <- unique(stgroup[!is.na(stgroup)]) - nstanzas <- as.vector(table(stgroup)[StanzaGroups]) + StanzaGroups <- unique(stgroup[!is.na(stgroup)]) + nstanzas <- as.vector(table(stgroup)[StanzaGroups]) NStanzaGroups <- length(StanzaGroups) Rpath.params$stanzas$NStanzaGroups <- NStanzaGroups - + stgroups <- data.table( - StGroupNum = 1:NStanzaGroups, + StGroupNum = 1:NStanzaGroups, StanzaGroup = StanzaGroups, - nstanzas = nstanzas, - VBGF_Ksp = NA, - VBGF_d = 0.66667, - Wmat = NA, - BAB = 0, - RecPower = 1 + nstanzas = nstanzas, + VBGF_Ksp = NA, + VBGF_d = 0.66667, + Wmat = NA, + BAB = 0, + RecPower = 1 ) - + #Individual Stanza Parameters ind.stanza.group <- model[!is.na(stgroup), Group] ieco <- which(!is.na(stgroup)) + # KYA Aug 2025 - The StGroupNum =rep() line below assumes that the + # each stanza grouping is "together" not interleaved - this is not + # guaranteed. The following lookup fixes that (data.frame used because + # data.table seems intractable to this). + stframe <- data.frame(stgroups) + row.names(stframe) <- stframe$StanzaGroup + gnum <- stframe[stgroup[!is.na(stgroup)], "StGroupNum"] stindiv <- data.table( - StGroupNum = rep(stgroups[, StGroupNum], stgroups[, nstanzas]), - StanzaNum = as.integer(0), - GroupNum = ieco, - Group = ind.stanza.group, - First = NA, - Last = NA, - Z = NA, - Leading = NA + StGroupNum = gnum, #rep(stgroups[, StGroupNum], stgroups[, nstanzas]), + StanzaNum = as.integer(0), + GroupNum = ieco, + Group = ind.stanza.group, + First = NA, + Last = NA, + Z = NA, + Leading = NA ) - } else { Rpath.params$stanzas$NStanzaGroups <- 0 - + stgroups <- data.table( - StGroupNum = NA, + StGroupNum = NA, StanzaGroup = NA, - nstanzas = NA, - VBGF_Ksp = NA, - VBGF_d = NA, - Wmat = NA, - RecPower = NA + nstanzas = NA, + VBGF_Ksp = NA, + VBGF_d = NA, + Wmat = NA, + RecPower = NA ) - + stindiv <- data.table( StGroupNum = NA, - StanzaNum = NA, - GroupNum = NA, - Group = NA, - First = NA, - Last = NA, - Z = NA, - Leading = NA + StanzaNum = NA, + GroupNum = NA, + Group = NA, + First = NA, + Last = NA, + Z = NA, + Leading = NA ) } - + Rpath.params$stanzas$stgroups <- stgroups - Rpath.params$stanzas$stindiv <- stindiv - + Rpath.params$stanzas$stindiv <- stindiv + #Pedigree pedigree <- data.table( - Group = group, + Group = group, Biomass = 1, - PB = 1, - QB = 1, - Diet = 1 + PB = 1, + QB = 1, + Diet = 1 ) #Add fleet pedigree for (i in 1:length(fleet.group)) { @@ -161,17 +177,18 @@ create.rpath.params <- function(group, type, stgroup = NA) { #' #'@return Checks Rpath parameter files for consistency. An error message will be produced if one of #' the logical checks fails. Checks include: verification that all types are represented (e.g. consumer, producer, detrital, and fleet); -#' check if input parameters were entered correctly; check if Diet columns sum to 1. +#' check if input parameters were entered correctly; check if Diet columns sum to 1. #' (NOTE: This does not ensure data is correct just that it is in the right places). -#' +#' #'@import data.table #'@export #' check.rpath.params <- function(Rpath.params) { #Need to define variables to eliminate check() note about no visible binding Type <- Group <- Biomass <- EE <- PB <- QB <- ProdCons <- BioAcc <- Unassim <- DetInput <- NULL - + w <- 0 #warning counter + c <- 0 #Balance change counter #Check to make sure all types are represented if (nrow(Rpath.params$model[Type == 0, ]) == 0) { warning('Model must contain at least 1 consumer') @@ -188,12 +205,12 @@ check.rpath.params <- function(Rpath.params) { if (nrow(Rpath.params$model[Type == 3, ]) == 0) { warning('Model must contain at least 1 fleet') } - + #Check that there is the proper number of columns n.groups <- nrow(Rpath.params$model) n.living <- length(Rpath.params$model[Type <= 1, Group]) - n.dead <- length(Rpath.params$model[Type == 2, Group]) - n.fleet <- length(Rpath.params$model[Type == 3, Group]) + n.dead <- length(Rpath.params$model[Type == 2, Group]) + n.fleet <- length(Rpath.params$model[Type == 3, Group]) if (ncol(Rpath.params$model) != 10 + n.dead + 2 * n.fleet) { warning( 'Model does not have the correct number of column. There should be 10 @@ -202,10 +219,17 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - + #Check that either biomass or EE is entered and not both - if (length(Rpath.params$model[is.na(Biomass) & - is.na(EE) & Type < 2, Group]) > 0) { + if ( + length(Rpath.params$model[ + is.na(Biomass) & + is.na(EE) & + Type < 2, + Group + ]) > + 0 + ) { warning( paste( Rpath.params$model[is.na(Biomass) & is.na(EE) & Type < 2, Group], @@ -215,22 +239,42 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - if (length(Rpath.params$model[!is.na(Biomass) & - !is.na(EE) & Type < 2, Group]) > 0) { + if ( + length(Rpath.params$model[ + !is.na(Biomass) & + !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', + Rpath.params$model[ + !is.na(Biomass) & + !is.na(EE) & + (!is.na(PB) | (is.na(PB) & !is.na(QB) & !is.na(ProdCons))) & + Type < 2, + Group + ], + 'have all of Biomass, EE, and PB(or QB and ProdCons) entered... Note that Rpath does + not calculate BA, please enter a value for BA if appropriate \n', sep = ' ' ) ) w <- w + 1 } - + #Check that Biomass / PB / QB / EE / ProdCons is not entered for types 3 - if (length(Rpath.params$model[Type == 3 & - !is.na(Biomass), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type == 3 & + !is.na(Biomass), + Group + ]) > + 0 + ) { warning( paste( Rpath.params$model[Type == 3 & !is.na(Biomass), Group], @@ -240,8 +284,14 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - if (length(Rpath.params$model[Type == 3 & - !is.na(PB), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type == 3 & + !is.na(PB), + Group + ]) > + 0 + ) { warning( paste( Rpath.params$model[Type == 3 & !is.na(PB), Group], @@ -255,7 +305,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,44 +315,70 @@ 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 = ' ' ) ) w <- w + 1 } - if (length(Rpath.params$model[Type > 1 & - !is.na(ProdCons), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type > 1 & + !is.na(ProdCons), + Group + ]) > + 0 + ) { 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 = ' ' ) ) w <- w + 1 } - + #Check that types 0 and 1 have a PB unless QB and ProdCons are entered 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) { + if ( + length(Rpath.params$model[ + Group %in% + no.pb & + (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', + Rpath.params$model[ + Group %in% + no.pb & + (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 = ' ' ) ) w <- w + 1 } } - + #Check that consumers have a QB or ProdCons but not both unless missing PB - if (length(Rpath.params$model[is.na(QB) & - is.na(ProdCons) & Type < 1, Group]) > 0) { + if ( + length(Rpath.params$model[ + is.na(QB) & + is.na(ProdCons) & + Type < 1, + Group + ]) > + 0 + ) { warning( paste( Rpath.params$model[is.na(QB) & is.na(ProdCons) & Type < 1, Group], @@ -312,26 +388,43 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - if (length(Rpath.params$model[!is.na(QB) & - !is.na(ProdCons) & Type < 1, Group]) > 0) { - both <- Rpath.params$model[!is.na(QB) & - !is.na(ProdCons) & Type < 1, Group] - if (length(Rpath.params$model[Group %in% both & - !is.na(PB), Group]) > 0) { + if ( + length(Rpath.params$model[ + !is.na(QB) & + !is.na(ProdCons) & + Type < 1, + Group + ]) > + 0 + ) { + both <- Rpath.params$model[ + !is.na(QB) & + !is.na(ProdCons) & + Type < 1, + Group + ] + if (length(Rpath.params$model[Group %in% both & !is.na(PB), Group]) > 0) { warning( paste( Rpath.params$model[Group %in% both & !is.na(PB), Group], - 'have PB, QB, and ProdCons...only two should be entered \n', + 'have PB, QB, and ProdCons... ProdCons will be recalculated during balancing \n', sep = ' ' ) ) w <- w + 1 + c <- c + 1 } } - + #Check that BioAcc / Unassim is NA for fleets and numeric for types < 3 - if (length(Rpath.params$model[Type == 3 & - !is.na(BioAcc), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type == 3 & + !is.na(BioAcc), + Group + ]) > + 0 + ) { warning( paste( Rpath.params$model[Type == 3 & !is.na(BioAcc), Group], @@ -341,8 +434,14 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - if (length(Rpath.params$model[Type == 3 & - !is.na(Unassim), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type == 3 & + !is.na(Unassim), + Group + ]) > + 0 + ) { warning( paste( Rpath.params$model[Type == 3 & !is.na(Unassim), Group], @@ -352,8 +451,14 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - if (length(Rpath.params$model[Type != 3 & - is.na(BioAcc), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type != 3 & + is.na(BioAcc), + Group + ]) > + 0 + ) { warning(paste( Rpath.params$model[Type != 3 & is.na(BioAcc), Group], 'must have a number for BioAcc...set to >= 0 \n', @@ -361,8 +466,14 @@ check.rpath.params <- function(Rpath.params) { )) w <- w + 1 } - if (length(Rpath.params$model[Type != 3 & - is.na(Unassim), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type != 3 & + is.na(Unassim), + Group + ]) > + 0 + ) { warning(paste( Rpath.params$model[Type != 3 & is.na(Unassim), Group], 'must have a number for Unassim...set to >= 0 \n', @@ -370,10 +481,16 @@ check.rpath.params <- function(Rpath.params) { )) w <- w + 1 } - + #Check that only Type 2 has DetInput set - if (length(Rpath.params$model[Type != 2 & - !is.na(DetInput), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type != 2 & + !is.na(DetInput), + Group + ]) > + 0 + ) { warning(paste( Rpath.params$model[Type != 2 & !is.na(DetInput), Group], 'are not detritus...set DetInput to NA \n', @@ -381,8 +498,14 @@ check.rpath.params <- function(Rpath.params) { )) w <- w + 1 } - if (length(Rpath.params$model[Type == 2 & - is.na(DetInput), Group]) > 0) { + if ( + length(Rpath.params$model[ + Type == 2 & + is.na(DetInput), + Group + ]) > + 0 + ) { warning(paste( Rpath.params$model[Type == 2 & is.na(DetInput), Group], 'are detritus...set DetInput to 0 \n', @@ -390,24 +513,31 @@ check.rpath.params <- function(Rpath.params) { )) w <- w + 1 } - + #Check detritus fate is numeric and sum to 1 det.matrix <- Rpath.params$model[, 11:(10 + n.dead), with = F] - test.rows <- rowSums(det.matrix) - if (length(setdiff(which(Rpath.params$model[, Type] == 2), which(test.rows != 1))) > 0) { - warning(paste( - Rpath.params$model[, Group][setdiff(which(Rpath.params$model[, Type] == 2), which(test.rows != 1))], - 'detrital fate does not sum to 1 \n', - sep = ' ' - )) - w <- w + 1 - } + + # KYA August '25 - I'm not sure why this part of this test should be done - having + # detrital groups' detrital fate not sum to 1 is appropriate if detritus + # is being exported. I think silence here is more appropriate than a warning. + #test.rows <- rowSums(det.matrix) + #if (length(setdiff(which(Rpath.params$model[, Type] == 2), which(test.rows != 1))) > 0) { + # warning(paste( + # Rpath.params$model[, Group][setdiff(which(Rpath.params$model[, Type] == 2), which(test.rows != 1))], + # 'detrital fate does not sum to 1 \n', + # sep = ' ' + # )) + # w <- w + 1 + #} + if (length(which(is.na(det.matrix))) > 0) { na.group <- which(is.na(det.matrix)) - for (i in 1:length(na.group)) - while (na.group[i] > n.groups) + for (i in 1:length(na.group)) { + while (na.group[i] > n.groups) { na.group[i] <- na.group[i] - n.groups + } + } na.group <- unique(na.group) warning( paste( @@ -418,15 +548,21 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - + #Check that landings and discards are numbers for type < 3 - fleet.matrix <- Rpath.params$model[1:(n.groups - n.fleet), (11 + n.dead):ncol(Rpath.params$model), with = F] + fleet.matrix <- Rpath.params$model[ + 1:(n.groups - n.fleet), + (11 + n.dead):ncol(Rpath.params$model), + with = F + ] if (length(which(is.na(fleet.matrix))) > 0) { na.group <- which(is.na(fleet.matrix)) - for (i in 1:length(na.group)) - while (na.group[i] > n.groups) + for (i in 1:length(na.group)) { + while (na.group[i] > n.groups) { na.group[i] <- na.group[i] - n.groups + } + } na.group <- unique(na.group) warning(paste( Rpath.params$model[na.group, Group], @@ -435,9 +571,13 @@ check.rpath.params <- function(Rpath.params) { )) w <- w + 1 } - + #Check that fleets aren't catching other fleets - fleet.matrix.2 <- Rpath.params$model[(n.groups - n.fleet + 1):n.groups, (11 + n.dead):ncol(Rpath.params$model), with = F] + fleet.matrix.2 <- Rpath.params$model[ + (n.groups - n.fleet + 1):n.groups, + (11 + n.dead):ncol(Rpath.params$model), + with = F + ] if (length(which(!is.na(fleet.matrix.2))) > 0) { not.na.group <- which(!is.na(fleet.matrix.2)) for (i in 1:length(not.na.group)) { @@ -453,12 +593,15 @@ check.rpath.params <- function(Rpath.params) { )) w <- w + 1 } - + #Diet #Check that columns sum to 1 col.names <- names(Rpath.params$diet)[2:ncol(Rpath.params$diet)] - col.sums <- Rpath.params$diet[, lapply(.SD, sum, na.rm = T), .SDcols = col.names] - + col.sums <- Rpath.params$diet[, + lapply(.SD, sum, na.rm = T), + .SDcols = col.names + ] + #Check types (>0 & <=1 are primary producers) types <- Rpath.params$model[Type < 2, Type] dctype <- round(col.sums + types, 3) @@ -475,7 +618,7 @@ check.rpath.params <- function(Rpath.params) { w <- w + 1 } } - + #Check number of columns dietcol <- ncol(Rpath.params$diet) if (dietcol != (n.living + 1)) { @@ -489,29 +632,39 @@ check.rpath.params <- function(Rpath.params) { ) w <- w + 1 } - + #Check that final row of diet is "Import" - if (!Rpath.params$diet[nrow(Rpath.params$diet), 1] == 'Import' & - !Rpath.params$diet[nrow(Rpath.params$diet), 1] == 'import') { + if ( + !Rpath.params$diet[nrow(Rpath.params$diet), 1] == 'Import' & + !Rpath.params$diet[nrow(Rpath.params$diet), 1] == 'import' + ) { warning( 'Diet matrix is missing the import row. Please add "Import" as the final row. All entries can be 0 or NA.' ) w <- w + 1 } - + # Check if any diet value is < 0 dietDF <- data.frame(Rpath.params$diet) dietDF[is.na(dietDF)] <- 0 if (any(dietDF < 0)) { - warning('Found a negative diet value. Please make sure no diet values are negative.') + warning( + 'Found a negative diet value. Please make sure no diet values are negative.' + ) w <- w + 1 } - + if (w == 0) { cat('Rpath parameter file is functional. \n') } else { - cat('Rpath parameter file needs attention! \n') + if (w == c) { + cat( + 'Rpath parameters functional, though some may be recalculated during balance. \n' + ) + } else { + cat('Rpath parameter file needs attention! \n') + } } } @@ -531,50 +684,55 @@ check.rpath.params <- function(Rpath.params) { #'@param pedfile file location of the flat file containing the \emph{pedgigree} parameters. #'@return Outputs an \code{Rpath.param} object that can be used for Rpath and subsequently #' Rsim. (NOTE: This does function does not ensure data is correct or in the -#' correct locations.Run \code{\link{check.rpath.param}()} to ensure the appropriate columns are +#' correct locations.Run \code{check.rpath.params()} to ensure the appropriate columns are #' present). #'@export -read.rpath.params <- function(modfile, - dietfile, - pedfile = NA, - stanzagroupfile = NA, - stanzafile = NA) { +read.rpath.params <- function( + modfile, + dietfile, + pedfile = NA, + stanzagroupfile = NA, + stanzafile = NA +) { #Need to define variables to eliminate check() note about no visible binding Type <- Group <- V1 <- NULL - + Rpath.params <- list() Rpath.params$model <- as.data.table(read.csv(modfile, header = T)) - Rpath.params$diet <- as.data.table(read.csv(dietfile, header = T)) - + Rpath.params$diet <- as.data.table(read.csv(dietfile, header = T)) + if (!is.na(stanzagroupfile)) { stanzagroup <- as.data.table(read.csv(stanzagroupfile, header = T)) Rpath.params$stanzas$NStanzaGroups <- nrow(stanzagroup) - Rpath.params$stanzas$stgroups <- stanzagroup - Rpath.params$stanzas$stindiv <- as.data.table(read.csv(stanzafile, header = T)) + Rpath.params$stanzas$stgroups <- stanzagroup + Rpath.params$stanzas$stindiv <- as.data.table(read.csv( + stanzafile, + header = T + )) } else { Rpath.params$stanzas$NStanzaGroups <- 0 - Rpath.params$stanzas$stgroups <- data.table( - StGroupNum = NA, + Rpath.params$stanzas$stgroups <- data.table( + StGroupNum = NA, StanzaGroup = NA, - nstanzas = NA, - VBGF_Ksp = NA, - VBGF_d = NA, - Wmat = NA, - RecPower = NA, - Wmat001 = NA, - Wmat50 = NA, - Amat001 = NA, - Amat50 = NA + nstanzas = NA, + VBGF_Ksp = NA, + VBGF_d = NA, + Wmat = NA, + RecPower = NA, + Wmat001 = NA, + Wmat50 = NA, + Amat001 = NA, + Amat50 = NA ) - Rpath.params$stanzas$stindiv <- data.table( - StGroupNum = NA, - StanzaNum = NA, - GroupNum = NA, - Group = NA, - First = NA, - Last = NA, - Z = NA, - Leading = NA + Rpath.params$stanzas$stindiv <- data.table( + StGroupNum = NA, + StanzaNum = NA, + GroupNum = NA, + Group = NA, + First = NA, + Last = NA, + Z = NA, + Leading = NA ) } if (!is.na(pedfile)) { @@ -582,10 +740,10 @@ read.rpath.params <- function(modfile, } else { Rpath.params$pedigree <- data.table( Group = Rpath.params$model$Group, - B = 1, - PB = 1, - QB = 1, - Diet = 1 + B = 1, + PB = 1, + QB = 1, + Diet = 1 ) fleets <- as.character(Rpath.params$model[Type == 3, Group]) for (i in 1:length(fleets)) { @@ -612,31 +770,37 @@ read.rpath.params <- function(modfile, #' "eco.name_model.csv". #'@export write.rpath.params <- function(Rpath.params, eco.name, path = '') { - write.csv(Rpath.params$model, - file = file.path(path, paste(eco.name, '_model.csv', sep = '')), - row.names = F) - - write.csv(Rpath.params$diet, - file = file.path(path, paste(eco.name, '_diet.csv', sep = '')), - row.names = F) - - write.csv(Rpath.params$pedigree, - file = file.path(path, paste(eco.name, '_pedigree.csv', sep = '')), - row.names = F) - + write.csv( + Rpath.params$model, + file = file.path(path, paste(eco.name, '_model.csv', sep = '')), + row.names = F + ) + + write.csv( + Rpath.params$diet, + file = file.path(path, paste(eco.name, '_diet.csv', sep = '')), + row.names = F + ) + + write.csv( + Rpath.params$pedigree, + file = file.path(path, paste(eco.name, '_pedigree.csv', sep = '')), + row.names = F + ) + #Multistanza parameters are in several different files write.csv( Rpath.params$stanzas$stgroups, file = file.path(path, paste(eco.name, '_stanza_groups.csv', sep = '')), row.names = F ) - + write.csv( Rpath.params$stanzas$stindiv, file = file.path(path, paste(eco.name, '_stanzas.csv', sep = '')), row.names = F ) - + if (Rpath.params$stanzas$NStanzaGroups > 0) { for (isp in 1:Rpath.params$stanzas$NStanzaGroups) { write.csv( diff --git a/R/xml_convert.r b/R/xml_convert.r new file mode 100644 index 00000000..66c4db26 --- /dev/null +++ b/R/xml_convert.r @@ -0,0 +1,574 @@ +#' Creates an Rpath object from an EwE exported model (EIIXML format) +#' +#' @description +#' Parses an eiixml file exported using the Ecopath with Ecosim (EwE) GUI +#' into an unbalanced Rpath model object, which can then be balanced. +#' This function was tested on files exported by EwE version 6.7. +#' +#' @param eiifile Path to exported EwE XML file +#' @param verbose Logical. Use for debugging. If TRUE, prints out a list of parsed tables and number of rows read from the XML file +#' +#' @return An Rpath object (list) with the following components: +#' \item{stanzas}{Stanza parameters} +#' \item{pedigree}{Pedigree parameters} +#' \item{diet}{Diet parameters} +#' \item{model}{Model parameters} +#' +#' @family xml +#' +#' @examples +#' # Import an eiixml file previously exported from the EwE GUI, check that the +#' # unbalanced model is functional, recalculate age-structured groups based on +#' # imported stanza parameters (if the imported model has stanzas), and balance +#' # the model. +#' eiixml <- system.file("extdata/xml/Western_Bering_Sea.eiixml", package = "Rpath") +#' rpath_object <- create.rpath.from.eiixml(eiixml) +#' check.rpath.params(rpath_object) +#' rpath_object <- rpath.stanzas(rpath_object) +#' rpath_balanced <- rpath(rpath_object) +#' +#' @export + +create.rpath.from.eiixml <- function( + eiifile, + verbose = FALSE +) { + # Import the xml file and parse it into a list of data frames + parsed_object <- import.eiixml(eiifile, verbose) + + # Extract and order group, gear and stanza names to create Rpath object--------- + + # Order of groups on spreadsheets/in Rpath is listed in the + # Sequence column. This is different than the GroupID order + # so a lookup is needed. GroupID is the key used for From and To + # flows so a lookup table is needed + seq_look <- parsed_object$ewe_EcopathGroup$Sequence + ord <- order(parsed_object$ewe_EcopathGroup$Sequence) + ford <- order(parsed_object$ewe_EcopathFleet$Sequence) + + # Sort data tables by sequence order + ordgroups <- parsed_object$ewe_EcopathGroup[ord, ] + ordfleets <- parsed_object$ewe_EcopathFleet[ford, ] + + # Get and clean names and Types for living groups and gear + bio_names <- janitor::make_clean_names(ordgroups$GroupName) + gear_names <- janitor::make_clean_names(ordfleets$FleetName) + + if (sum(gear_names %in% bio_names) > 0) { + warning( + "There are Fleets with the same name as a bio group. + These fleets will be renamed to avoid confusion. + The string `_fleet` will be appended.\n + Fleet names that match bio group names: ", + paste0(gear_names[gear_names %in% bio_names], collapse = ",") + ) + } + gear_names <- ifelse( + gear_names %in% bio_names, + paste(gear_names, "fleet", sep = '_'), + gear_names + ) + det_names <- janitor::make_clean_names(ordgroups$GroupName[ + ordgroups$Type == 2 + ]) + live_names <- janitor::make_clean_names(ordgroups$GroupName[ + ordgroups$Type != 2 + ]) + + g_names <- c(bio_names, gear_names) + g_types <- c(ordgroups$Type, rep(3, length(gear_names))) + row.names(ordgroups) <- bio_names + + # Also make a couple of named vectors associating names with groupID + # (lookups for species and gear, using ID# as a character index) + pnames <- bio_names + names(pnames) <- ordgroups$GroupID + gnames <- gear_names + names(gnames) <- ordfleets$FleetID + + # Create the stanza table only if the model has stanzas defined in the xml + + if (!(identical(parsed_object$ewe_Stanza, NA))) { + stanza_info <- make_stanza_table( + parsed_object$ewe_Stanza, + parsed_object$ewe_StanzaLifeStage, + pnames, + gnames + ) + stanza_name_only <- stanza_info$stanza_name_only + ordstanzas <- stanza_info$ordstanzas + ordstages <- stanza_info$ordstages + } else { + # create a default stanza object + stanza_name_only <- NULL + } + + #------------------------------------------------------------------------------- + # Create the Unbalanced Rpath object here + unbal <- create.rpath.params( + group = g_names, + type = g_types, + stgroup = stanza_name_only + ) + + # return(list( + # g_names = g_names, + # g_types = g_types, + # stanza_names_only = stanza_name_only + # )) + # ------------------------------------------------------------------------------ + # Populate unbal Rpath object with values + + # Fill group vectors + gear_na <- rep(NA, length(gear_names)) + unbal$model$Biomass <- as.numeric(c(vec_na(ordgroups$Biomass), gear_na)) + unbal$model$PB <- as.numeric(c(vec_na(ordgroups$ProdBiom), gear_na)) + unbal$model$QB <- as.numeric(c(vec_na(ordgroups$ConsBiom), gear_na)) + unbal$model$EE <- as.numeric(c(vec_na(ordgroups$EcoEfficiency), gear_na)) + unbal$model$ProdCons <- as.numeric(c(vec_na(ordgroups$ProdCons), gear_na)) + unbal$model$BioAcc <- as.numeric(c(vec_na(ordgroups$BiomAcc), gear_na)) + unbal$model$Unassim <- as.numeric(c(vec_na(ordgroups$Unassim), gear_na)) + unbal$model$DetInput <- as.numeric(c(vec_na(ordgroups$DtImports), gear_na)) + + #EwE output seems to have a few 0s or other numbers saved that should + #be NAs, specifically with detritus. Ensuring those don't sneak through. + unbal$model$DetInput[ordgroups$Type != 2] <- NA + unbal$model$EE[ordgroups$Type == 2] <- NA + + # DIET TABLE--------------------------------------- + # TODO: where are diet imports in EwE XML data? + diet_table <- parsed_object$ewe_EcopathDietComp + diet_table$pred_name <- pnames[as.character( + parsed_object$ewe_EcopathDietComp$PredID + )] + diet_table$prey_name <- pnames[as.character( + parsed_object$ewe_EcopathDietComp$PreyID + )] + + ppmat <- data.frame(unbal$diet[, -1]) + row.names(ppmat) <- unbal$diet$Group + for (i in 1:length(diet_table$Diet)) { + #cat(i,diet_table$prey_name[i], diet_table$pred_name[i],"\n") + if ( + diet_table$prey_name[i] %in% + rownames(ppmat) & + diet_table$pred_name[i] %in% colnames(ppmat) + ) { + ppmat[ + diet_table$prey_name[i], + diet_table$pred_name[i] + ] <- diet_table$Diet[i] + } + } + # need to convert matrices to data frames or data.table is unhappy + unbal$diet[, 2:ncol(unbal$diet)] <- ppmat + + # add diet import + impdiet <- ordgroups$ImpVar + names(impdiet) <- rownames(ordgroups) + predlist <- names(unbal$diet)[2:ncol(unbal$diet)] + for (p in predlist) { + unbal$diet[Group == "Import", p] <- impdiet[p] + } + + # DETRITUS FATE FOR NON-GEAR ------------------------------------------------- + detframe <- data.frame(unbal$model)[, det_names] + if (is.null(dim(detframe))) { + detframe <- as.data.frame(detframe) + colnames(detframe) <- det_names + } + row.names(detframe) <- g_names + for (i in 1:length(diet_table$Diet)) { + if (diet_table$prey_name[i] %in% det_names) { + detframe[ + diet_table$pred_name[i], + diet_table$prey_name[i] + ] <- diet_table$DetritusFate[i] + } + } + + # eiixmls have some detritus fates NA/blank, replace with 0s + detframe[is.na(detframe)] <- 0 + + unbal$model[, det_names] <- detframe + + # CATCH AND DISCARDS -------------------------------------------------------- + if (any(is.na(parsed_object$ewe_EcopathCatch))) { + catch_table <- data.frame( + DiscardMortality = 1, + Discards = 0.0, + FleetID = 1, + GroupID = seq(1, (length(pnames))), + Landing = 0.0, + Price = 0.0 + ) + } else { + catch_table <- parsed_object$ewe_EcopathCatch + } + catch_table$gear_name <- gnames[as.character(catch_table$FleetID)] + catch_table$group_name <- pnames[as.character(catch_table$GroupID)] + + discard_names <- paste(gear_names, "disc", sep = '.') + + fishframe <- data.frame(unbal$model)[, gear_names] + if (is.null(dim(fishframe))) { + fishframe <- as.data.frame(fishframe) + colnames(fishframe) <- gear_names + } + + row.names(fishframe) <- unbal$model$Group + discardframe <- data.frame(unbal$model)[, discard_names] + if (is.null(dim(discardframe))) { + discardframe <- as.data.frame(discardframe) + colnames(discardframe) <- discard_names + } + row.names(discardframe) <- unbal$model$Group + + for (i in 1:length(catch_table$Landing)) { + if ( + catch_table$gear_name[i] %in% + colnames(fishframe) & + catch_table$group_name[i] %in% rownames(fishframe) + ) { + fishframe[ + catch_table$group_name[i], + catch_table$gear_name[i] + ] <- catch_table$Landing[i] + discardframe[ + catch_table$group_name[i], + paste(catch_table$gear_name[i], "disc", sep = '.') + ] <- catch_table$Discards[i] + } + } + + unbal$model[, gear_names] <- fishframe + unbal$model[, discard_names] <- discardframe + + # DETRITUS FATE FOR GEAR ----------------------------------------------------- + if (any(is.na(parsed_object$ewe_EcopathDiscardFate))) { + # If no fleets/catch, route fleet1 detritus to last detritus group + last_detritus <- as.numeric(names(which( + pnames == det_names[length(det_names)] + ))) + fate_table <- data.frame( + DiscardFate = 1.0, + FleetID = 1, + GroupID = last_detritus + ) + } else { + fate_table <- parsed_object$ewe_EcopathDiscardFate + } + fate_table$gear_name <- gnames[as.character(fate_table$FleetID)] + fate_table$group_name <- pnames[as.character(fate_table$GroupID)] + + fateframe <- data.frame(unbal$model)[unbal$model$Group %in% gnames, det_names] + if (is.null(dim(fateframe))) { + fateframe <- as.data.frame(fateframe) + colnames(fateframe) <- det_names + } + + row.names(fateframe) <- gnames + for (i in 1:length(fate_table$DiscardFate)) { + if ( + fate_table$gear_name[i] %in% + rownames(fateframe) & + fate_table$group_name[i] %in% colnames(fateframe) + ) { + fateframe[ + fate_table$gear_name[i], + fate_table$group_name[i] + ] <- fate_table$DiscardFate[i] + } + } + + unbal$model[ + (length(bio_names) + 1):(length(bio_names) + length(gear_names)), + det_names + ] <- fateframe + + # Configure the stanzas object + if (!(identical(parsed_object$ewe_Stanza, NA))) { + # STANZAS ---------------------------------------------------------------------- + + unbal$stanzas$stgroups$Wmat <- ordstanzas[ + unbal$stanzas$stgroups$StanzaGroup, + "WmatWinf" + ] + unbal$stanzas$stgroups$BAB <- ordstanzas[ + unbal$stanzas$stgroups$StanzaGroup, + "BABsplit" + ] + unbal$stanzas$stgroups$RecPower <- ordstanzas[ + unbal$stanzas$stgroups$StanzaGroup, + "RecPower" + ] + + row.names(ordstages) <- ordstages$gname + unbal$stanzas$stindiv$StanzaNum <- ordstages[ + unbal$stanzas$stindiv$Group, + "Sequence" + ] + unbal$stanzas$stindiv$Leading <- ifelse( + ordstages[unbal$stanzas$stindiv$Group, "LeadingLifeStage"] == + ordstages[unbal$stanzas$stindiv$Group, "Sequence"], + TRUE, + FALSE + ) + unbal$stanzas$stindiv$First <- ordstages[ + unbal$stanzas$stindiv$Group, + "AgeStart" + ] + + # Not sure how Mort is stored outside of PB, or if it's different. + pb <- unbal$model$PB + names(pb) <- unbal$model$Group + unbal$stanzas$stindiv$Z <- pb[unbal$stanzas$stindiv$Group] + + # Looping through in this weird way doesn't change the data.tables data type + # (default logical) for the Last column. (more foolishness) + unbal$stanzas$stindiv[, Last := as.numeric(Last)] + unbal$stanzas$stgroups[, VBGF_Ksp := as.numeric(VBGF_Ksp)] + + for (gg in 1:max(unbal$stanzas$stindiv$StGroupNum)) { + # Finding vonBK is a mess involving all three tables (stored on group table) + # Look up Leading Group and Stanza Group Number to get name on main table. + unbal$stanzas$stgroups[StGroupNum == gg, "VBGF_Ksp"] <- + ordgroups[ + as.character(unbal$stanzas$stindiv[ + StGroupNum == gg & Leading == T, + "Group" + ]), + ]$vbK + + # Now loop through to add Last Month to first month + szs <- unbal$stanzas$stindiv[StGroupNum == gg, ] + for (ss in 1:(max(szs$StanzaNum) - 1)) { + unbal$stanzas$stindiv[ + StGroupNum == gg & StanzaNum == ss, + "Last" + ] <- szs[szs$StanzaNum == ss + 1, "First"] - 1 + } + ss <- ss + 1 + unbal$stanzas$stindiv[ + StGroupNum == gg & StanzaNum == ss, + "Last" + ] <- 999 + } + + # TODO RPATH - why does order matter here? + unbal$stanzas$stindiv <- unbal$stanzas$stindiv[order(StGroupNum, StanzaNum)] + #Not sure why Ksp is stored on the indiv table not the main group table + #assuming the value for the leading stanzas is correct. + } + if (0) { + cat("Created unbal (unbalanced ecopath) from ", eiifile, "\n") + } + return(unbal) +} + +############################################################################ +############################################################################ +############################################################################ +############################################################################ +#' Reads in EwE exported XML file and parses into data frames +#' +#' @description +#' Parses an eiixml file exported using the Ecopath with Ecosim (EwE) GUI into a list of data frames, +#' one frame for each table in the exported XML file. This function is usually called by the function \code{create.rpath.from.eiixml()} that +#' converts these tables into an unbalanced rpath model object. However import.xml can be used on its own to examine the full set +#' of tables exported by EwE, including tables not currently imported into Rpath objects, such as Ecosim runs or model metadata. +#' This function was tested on files exported by EwE version 6.7. +#' +#' @inheritParams create.rpath.from.eiixml +#' +#' @return A list of data frames, one data frame for each node (exported EwE table) in the XML file. Each table +#' has the naming convention ewe_[table name] where [table name] is the name of the table provided by EwE. +#' +#' @family xml +#' +#' @examples +#' # Import an eiixml file previously exported from the EwE GUI into a list of +#' # data frames containing the model data +#' eiixml <- system.file("extdata/xml","Western_Bering_Sea.eiixml", package = "Rpath") +#' xml_data <- import.eiixml(eiixml) +#' +#' +#' @export + +import.eiixml <- function(eiifile, verbose = FALSE) { + # Warn as you go, not at the end + options(warn = 1) + + # create list object used to store output + eweobject <- list() + + # Read and parse xml file + dat <- xml2::read_xml(eiifile) + + # Create a list of Tables (main eiixml output data structure is "Table") + tables <- xml2::xml_find_all(dat, ".//Table") + + # Loop through the nodes (Tables) and read into data frames + for (node in tables) { + # how to select a named node from the list of tables: + node_name <- xml2::xml_attr(node, "Name") + #cat(node_name,""); #flush.console() + # Parse the Columns attribute of each Table which is var name/var type. + # This messy command splits column names by , and by : then discarding + # the : which is variable type. I'm not really sure how the `[[` works but + # it's choosing the 1st element from each list. + cols <- unlist(lapply( + strsplit(unlist(strsplit(xml2::xml_attr(node, "Columns"), ",")), ":"), + `[[`, + 1 + )) + # Get rows of data and split into a data frame with that data + rawdat <- xml2::xml_text(xml2::xml_find_all(node, "Row")) + + # special reformatting of some nodes node to avoid comma issues + # We can't split by commas, like we can for other nodes, because + # there are some fields with commas in the text (within quotes), + # so we need to grep these lines, + # and replace the commas with whitespace for all occurrences + + if ( + node_name %in% + c( + "EcopathModel", + "Pedigree", + "EcopathGroup", + "EcopathFleet", + "Auxillary", + "UpdateLog" + ) && + (length(rawdat) > 0) + ) { + input_string <- rawdat + for (irow in 1:length(input_string)) { + row_string <- input_string[irow] + match_list <- gregexpr("\"[^\"]*\"", row_string) + quoted_content <- regmatches(row_string, match_list) + modified_quoted_content <- lapply(quoted_content, function(x) { + gsub(",", " ", x, fixed = TRUE) + }) + modified_string <- row_string + regmatches(modified_string, match_list) <- modified_quoted_content + row_string <- paste0(modified_string, " ") + rawdat[irow] <- row_string + } + } + + if (length(rawdat) > 0) { + if (verbose) { + cat(node_name, "") + } + #dtable <- data.frame(matrix(as.numeric(unlist(strsplit(rawdat,","))), ncol=length(cols), byrow=T)) + unlisted_dat <- (unlist(strsplit(rawdat, ","))) + if (length(unlisted_dat) %% length(cols) != 0) { + warning("Matrix alignment (comma issue?) in table ", node_name, ":") + } + dtable <- type.convert( + data.frame(matrix(unlisted_dat, ncol = length(cols), byrow = T)), + as.is = T + ) + if (verbose) { + cat(length(rawdat), "rows\n") + } + names(dtable) <- cols + eweobject[[paste0("ewe_", node_name)]] <- dtable + #assign(paste0("eweobject$ewe_", node_name), dtable) + } else { + dtable <- data.frame(matrix(ncol = length(cols), byrow = T))[-1, ] + # Assign default value of NA to all variables that dont exist for model + eweobject[[paste0("ewe_", node_name)]] <- NA + #assign(paste0("eweobject$ewe_", node_name), NULL) + #cat("0\n") + } + # Use these here if you want to make empty tables where there's no data + #names(dtable) <- cols + #assign(paste("ewe_",node_name,sep="") , dtable) + } # end node loop + + ## Create dummy values for missing components + if (identical(eweobject$ewe_EcopathFleet, NA)) { + warning("No Fleets are present. A dummy fleet will be created") + eweobject$ewe_EcopathFleet <- data.frame( + FixedCost = 0, + FleetID = 1, + FleetName = "Fleet1", + NominalEffort = NA, + PoolColor = 0, + SailingCost = 1, + Sequence = 1, + VariableCost = 0 + ) + } + + return(eweobject) + + ## ------ ALL RAW DATA from XML should be in R data frames at this point ------- +} + +############################################################################ +############################################################################ +############################################################################ +############################################################################ +#' Function to replace negative values with NAs +#' +#' Utility function +#' @noRd +vec_na <- function(vec) { + return(ifelse(vec < (-0.1), NA, vec)) +} +############################################################################ +############################################################################ +############################################################################ +############################################################################ +#' Create the Stanza table +#' +#' @description +#' Creates a table of stanzas and their life stages from the EwE Ecobase eiixml +#' +#' @param ewe_Stanza Stanza data from XML node +#' @param ewe_StanzaLifeStage Stanza life stage data from XML node +#' @param pnames Names of biological groups in the model +#' @param gnames Names of gear groups in the model +#' +#' @return list +#' \item{ordstanzas}{} +#' \item{stanza_name_only}{} +#' \item{ordstages}{} +#' +#' @family xml +#' +#' @noRd +make_stanza_table <- function(ewe_Stanza, ewe_StanzaLifeStage, pnames, gnames) { + # Get and clean stanza names, order stanza table + sord <- order(ewe_Stanza$StanzaID) + ordstanzas <- ewe_Stanza[sord, ] + ordstanzas$stanza_names <- janitor::make_clean_names(ordstanzas$StanzaName) + row.names(ordstanzas) <- as.character(ordstanzas$StanzaID) + + # Combine the life stages and stanza info into one table (making some duplication + # as stanza parameters will have multiple copies). + ordstages <- cbind( + ewe_StanzaLifeStage, + ordstanzas[as.character(ewe_StanzaLifeStage$StanzaID), ] + ) + + ordstages$gname <- pnames[as.character(ordstages$GroupID)] + + stanza_column <- rep(NA, length(pnames)) + names(stanza_column) <- names(pnames) + stanza_column[as.character(ordstages$GroupID)] <- ordstages$stanza_names + stanza_name_only <- c(as.character(stanza_column), rep(NA, length(gnames))) + + # Rename ordered tables to use cleaned names, not GroupIDs + row.names(ordstanzas) <- ordstanzas$stanza_names + + return(list( + ordstanzas = ordstanzas, + stanza_name_only = stanza_name_only, + ordstages = ordstages + )) +} diff --git a/_pkgdown.yml b/_pkgdown.yml index d7e7f5fe..10bced1e 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -66,11 +66,12 @@ reference: - write.Rsim - frate.table - extract.node - - Rpath-package - title: "Auxiliary Tools" desc: "Additional tools" - contents: - MTI + - import.eiixml + - create.rpath.from.eiixml navbar: @@ -86,5 +87,7 @@ navbar: href: articles/RunRsim.html - text: Create an ensemble href: articles/ecosense.html + - text: Convert EwE to Rpath + href: articles/Convert_EwE_to_Rpath.html - text: Unit testing - href: articles/UnitTests.html \ No newline at end of file + href: articles/UnitTests.html diff --git a/inst/extdata/xml/Western Bering Sea-Basic estimates.csv b/inst/extdata/xml/Western Bering Sea-Basic estimates.csv new file mode 100644 index 00000000..1092a420 --- /dev/null +++ b/inst/extdata/xml/Western Bering Sea-Basic estimates.csv @@ -0,0 +1,37 @@ +,"Group name","Trophic level","Hab area (proportion)","Biomass in habitat area (t/km^2)","Biomass (t/km^2)","Production / biomass (/year)","Consumption / biomass (/year)","Ecotrophic Efficiency","Production / consumption (/year)","Biomass accumulation (t/km^2)","BA rate (/year)" +1,"Baleen whales",3.77382302,1.00000000,0.39100000,0.39100000,0.020000000,8.31200027,0.98913050,0.00240616,, +2,"Toothed whales",4.15545177,1.00000000,0.041999999,0.041999999,0.020000000,17.5000000,0.00000000,0.00114286,, +3,"Sperm whales",4.70328903,1.00000000,0.020000000,0.020000000,0.020000000,9.80000019,0.00000000,0.00204082,, +4,"Walrus& bearded seals",3.15126300,1.00000000,0.26199999,0.26199999,0.059999999,13.8999996,0.42843521,0.00431655,, +5,Seals,3.91453004,1.00000000,0.097000003,0.097000003,0.059999999,13.3400002,0.59622002,0.00449775,, +6,"Steller sea lions",3.99049592,1.00000000,0.035000000,0.035000000,0.059999999,18.0000000,0.35000011,0.00333333,, +7,Seabirds,3.95844197,1.00000000,0.0099999998,0.0099999998,0.80000001,87.0000000,0.00000000,0.00919540,, +8,"Adult pollock",3.41674399,1.00000000,15.0000000,15.0000000,0.50000000,10.00000000,0.94974858,0.05000000,, +9,"Juvenile pollock",3.35486889,1.00000000,3.75699997,3.75699997,2.50000000,13.0000000,0.48001531,0.19230770,, +10,"Pacific cod",4.02525377,1.00000000,3.18700004,3.18700004,0.52100003,3.29999995,0.35416141,0.15787880,, +11,P.halibut,4.55041885,1.00000000,0.082999997,0.082999997,0.25000000,3.50000000,0.98781377,0.07142857,, +12,"Greenland turbot",4.45009899,1.00000000,0.057999998,0.057999998,0.20000000,3.59999990,0.94655168,0.05555556,, +13,"Arrowtooth flounder",4.27181721,1.00000000,0.052000001,0.052000001,0.31999999,4.50000000,0.88955140,0.07111111,, +14,"Small flatfish",3.20804310,1.00000000,0.99199998,0.99199998,0.28799999,6.84899998,0.73974621,0.04204993,, +15,Skates,4.34995079,1.00000000,0.27100000,0.27100000,0.40000001,4.00000000,0.37822881,0.10000000,, +16,"Sculpins& Rockfish",3.75946808,1.00000000,0.67699999,0.67699999,0.40000001,3.50000000,0.94447428,0.11428570,, +17,Macrouridae,3.90579796,1.00000000,1.15600002,1.15600002,0.30000001,3.70000005,0.67480099,0.08108108,, +18,Zoarcidae,4.08010292,1.00000000,0.89999998,0.89999998,0.30000001,2.51999998,0.98806512,0.11904760,, +19,"Tanner crab",2.99084091,1.00000000,0.082999997,0.082999997,0.80000001,5.00000000,0.39797691,0.16000000,, +20,"Snow crab",2.99084091,1.00000000,0.24900000,0.24900000,0.80000001,5.00000000,0.87123322,0.16000000,, +21,"King crab",2.99432397,1.00000000,0.11900000,0.11900000,0.60000002,2.70000005,0.82682788,0.22222219,, +22,Shrimp,2.40576601,1.00000000,2.10575294,2.10575294,2.03999996,10.1999998,0.90000021,0.20000000,, +23,Epifauna,2.15079594,1.00000000,114.961998,114.961998,1.15900004,5.08699989,0.24500690,0.22783570,, +24,Infauna,2.00000000,1.00000000,125.686996,125.686996,1.97000003,12.0000000,0.27408451,0.16416670,, +25,"Benthic amphipods",2.00502491,1.00000000,13.8120003,13.8120003,2.50000000,14.0000000,0.99449003,0.17857140,, +26,"Pacific herring",3.33645201,1.00000000,0.78700000,0.78700000,0.69999999,14.6000004,0.82282239,0.04794521,, +27,"Pacific salmon",3.67482901,1.00000000,0.039000001,0.039000001,4.00000000,16.0000000,0.39076331,0.25000000,, +28,Cephalopods,3.67573190,1.00000000,4.82999992,4.82999992,3.20000005,10.6700001,0.91090471,0.29990631,, +29,"Forage fish",3.28703189,1.00000000,19.41937065,19.41937065,0.94999999,3.50000000,0.89999992,0.27142861,, +30,Jellyfish,3.14681292,1.00000000,1.39999998,1.39999998,1.50000000,3.00000000,0.01646648,0.50000000,, +31,"Large Zooplankton",2.63895106,1.00000000,120.739998,120.739998,4.39900017,14.4569998,0.55201441,0.30428171,, +32,Copepods,2.05263209,1.00000000,122.620003,122.620003,9.50000000,26.2000008,0.97425568,0.36259541,, +33,Phytoplankton,1.00000000,1.00000000,15.0000000,15.0000000,234.000000,,0.73821777,,, +34,"Pelagic detritus",1.00000000,1.00000000,,,,,,,, +35,"Benthic detritus",1.00000000,1.00000000,,,,,,,, +36,"Detrital pool",1.00000000,1.00000000,,,,,0.00015627,,, diff --git a/inst/extdata/xml/Western_Bering_Sea.eiixml b/inst/extdata/xml/Western_Bering_Sea.eiixml new file mode 100644 index 00000000..80b3c9fc --- /dev/null +++ b/inst/extdata/xml/Western_Bering_Sea.eiixml @@ -0,0 +1,422 @@ + + + +
+ 1,0,1,1,0.007,0 + 1,0,1,2,0,0 + 1,0,1,3,0,0 + 1,0,1,4,0.006,0 + 1,0,1,5,0.002,0 + 1,0,1,6,0,0 + 1,0,1,7,0,0 + 1,0,1,8,1.051,0 + 1,0,1,9,0,0 + 1,0,1,10,0.226,0 + 1,0,1,11,0.01,0 + 1,0,1,12,0.01,0 + 1,0,1,13,0.007,0 + 1,0,1,14,0.041,0 + 1,0,1,15,0.041,0 + 1,0,1,16,0.07,0 + 1,0,1,17,0.041,0 + 1,0,1,18,0.039,0 + 1,0,1,19,0,0 + 1,0,1,20,0.006,0 + 1,0,1,21,0.008,0 + 1,0,1,22,0.002,0 + 1,0,1,23,0.004,0 + 1,0,1,24,0,0 + 1,0,1,25,0,0 + 1,0,1,26,0.055,0 + 1,0,1,27,0.012,0 + 1,0,1,28,0.02,0 + 1,0,1,29,0.001,0 + 1,0,1,30,0,0 + 1,0,1,31,0,0 + 1,0,1,32,0,0 + 1,0,1,33,0,0 + 1,0,1,34,0,0 + 1,0,1,35,0,0 + 1,0,1,36,0,0 +
+ +
+ 0,0.01,0,0,8,28 + 0,0.445,0,0,8,32 + 1,0,0,0,8,36 + 0,0.04,0,0,8,29 + 0,0.005,0,0,8,9 + 0,0.495,0,0,8,31 + 0,0.005,0,0,8,22 + 0,0.5,0,0,13,8 + 0,0.2,0,0,13,28 + 1,0,0,0,13,36 + 0,0.1,0,0,13,9 + 0,0.05,0,0,13,31 + 0,0.15,0,0,13,22 + 0,0.04,0,0,1,8 + 0,0.1,0,0,1,28 + 0,0.27,0,0,1,32 + 1,0,0,0,1,36 + 0,0.13,0,0,1,29 + 0,0.06,0,0,1,9 + 0,0.315,0,0,1,31 + 0,0.02,0,0,1,10 + 0,0.02,0,0,1,26 + 0,0.045,0,0,1,22 + 0,0.005,0,0,25,25 + 0,0.995,0,0,25,35 + 1,0,0,0,25,36 + 0,0.1114488,0,0,28,25 + 0,0.1236069,0,0,28,28 + 0,0.2066869,0,0,28,32 + 1,0,0,0,28,36 + 0,0.1550152,0,0,28,29 + 0,0.4032421,0,0,28,31 + 0,0.05,0,0,32,32 + 1,0,0,0,32,36 + 0,0.25,0,0,32,34 + 0,0.7,0,0,32,33 + 0,0.001996008,0,0,23,25 + 0,0.8562874,0,0,23,35 + 1,0,0,0,23,36 + 0,0.04690619,0,0,23,23 + 0,0.09481038,0,0,23,24 + 0,0.306245,0,0,29,25 + 0,0.01421137,0,0,29,28 + 0,0.2942354,0,0,29,32 + 1,0,0,0,29,36 + 0,0.3853082,0,0,29,31 + 0,0.45,0,0,12,8 + 0,0.3,0,0,12,28 + 1,0,0,0,12,36 + 0,0.1,0,0,12,29 + 0,0.01,0,0,12,17 + 0,0.01,0,0,12,16 + 0,0.07,0,0,12,22 + 0,0.01,0,0,12,14 + 0,0.05,0,0,12,18 + 0,1,0,0,24,35 + 1,0,0,0,24,36 + 0,0.6,0,0,30,32 + 1,0,0,0,30,36 + 0,0.01,0,0,30,9 + 0,0.3,0,0,30,31 + 0,0.05,0,0,30,34 + 0,0.04,0,0,30,33 + 0,0.02,0,0,9,28 + 0,0.45,0,0,9,32 + 1,0,0,0,9,36 + 0,0.5050001,0,0,9,31 + 0,0.025,0,0,9,34 + 0,0.1,0,0,21,25 + 0,0.1,0,0,21,35 + 1,0,0,0,21,36 + 0,0.38,0,0,21,23 + 0,0.33,0,0,21,24 + 0,0.09,0,0,21,22 + 0,0.484,0,0,31,32 + 1,0,0,0,31,36 + 0,0.079,0,0,31,31 + 0,0.241,0,0,31,34 + 0,0.196,0,0,31,33 + 0,0.05,0,0,17,25 + 0,0.3,0,0,17,28 + 1,0,0,0,17,36 + 0,0.1,0,0,17,23 + 0,0.16,0,0,17,29 + 0,0.2,0,0,17,24 + 0,0.12,0,0,17,31 + 0,0.04,0,0,17,17 + 0,0.02,0,0,17,22 + 0,0.01,0,0,17,18 + 0,0.386,0,0,10,8 + 0,0.04,0,0,10,25 + 0,0.03,0,0,10,28 + 1,0,0,0,10,36 + 0,0.126,0,0,10,23 + 0,0.021,0,0,10,29 + 0,0.07500001,0,0,10,24 + 0,0.2,0,0,10,9 + 0,0.001,0,0,10,21 + 0,0.02,0,0,10,31 + 0,0.005,0,0,10,26 + 0,0.015,0,0,10,16 + 0,0.05,0,0,10,22 + 0,0.01,0,0,10,14 + 0,0.01,0,0,10,20 + 0,0.001,0,0,10,19 + 0,0.01,0,0,10,18 + 0,0.15,0,0,26,25 + 0,0.35,0,0,26,32 + 1,0,0,0,26,36 + 0,0.03,0,0,26,23 + 0,0.03,0,0,26,29 + 0,0.41,0,0,26,31 + 0,0.03,0,0,26,22 + 0,0.095,0,0,27,25 + 0,0.095,0,0,27,28 + 0,0.1,0,0,27,32 + 1,0,0,0,27,36 + 0,0.025,0,0,27,23 + 0,0.12,0,0,27,29 + 0,0.525,0,0,27,31 + 0,0.04,0,0,27,22 + 0,0.19,0,0,11,8 + 0,0.5,0,0,11,28 + 1,0,0,0,11,36 + 0,0.025,0,0,11,29 + 0,0.17,0,0,11,9 + 0,0.06,0,0,11,10 + 0,0.005,0,0,11,16 + 0,0.015,0,0,11,22 + 0,0.02,0,0,11,14 + 0,0.005,0,0,11,20 + 0,0.005,0,0,11,19 + 0,0.005,0,0,11,18 + 1,0,0,0,33,36 + 0,0.2531468,0,0,16,8 + 0,0.05012908,0,0,16,25 + 0,0.0100058,0,0,16,28 + 1,0,0,0,16,36 + 0,0.2111225,0,0,16,23 + 0,0.005012908,0,0,16,29 + 0,0.1410818,0,0,16,24 + 0,0.131076,0,0,16,9 + 0,0.00100058,0,0,16,21 + 0,0.04512618,0,0,16,31 + 0,0.02001161,0,0,16,10 + 0,0.02001161,0,0,16,26 + 0,0.005012908,0,0,16,16 + 0,0.08124714,0,0,16,22 + 0,0.01500871,0,0,16,14 + 0,0.0100058,0,0,16,20 + 0,0.00100058,0,0,16,19 + 0,0.05,0,0,7,28 + 0,0.1,0,0,7,32 + 1,0,0,0,7,36 + 0,0.19,0,0,7,29 + 0,0.03,0,0,7,24 + 0,0.01,0,0,7,30 + 0,0.25,0,0,7,9 + 0,0.32,0,0,7,31 + 0,0.04,0,0,7,26 + 0,0.01,0,0,7,27 + 0,0.18,0,0,5,8 + 0,0.2,0,0,5,28 + 1,0.228,0,0,5,36 + 0,0.07500001,0,0,5,29 + 0,0.02,0,0,5,30 + 0,0.15,0,0,5,9 + 0,0.01,0,0,5,31 + 0,0.05,0,0,5,10 + 0,0.07000001,0,0,5,26 + 0,0.017,0,0,5,27 + 0,0.1,0,0,22,25 + 0,0.6,0,0,22,35 + 0,0.1,0,0,22,32 + 1,0,0,0,22,36 + 0,0.2,0,0,22,24 + 0,0.37,0,0,15,8 + 0,0.2,0,0,15,28 + 1,0,0,0,15,36 + 0,0.06,0,0,15,23 + 0,0.1,0,0,15,29 + 0,0.01,0,0,15,24 + 0,0.2,0,0,15,9 + 0,0.033,0,0,15,10 + 0,0.02,0,0,15,22 + 0,0.005,0,0,15,20 + 0,0.002,0,0,15,19 + 0,0.059,0,0,14,25 + 0,0.01,0,0,14,28 + 0,0.049,0,0,14,32 + 1,0,0,0,14,36 + 0,0.07000001,0,0,14,23 + 0,0.006000001,0,0,14,29 + 0,0.534,0,0,14,24 + 0,0.021,0,0,14,9 + 0,0.094,0,0,14,31 + 0,0.145,0,0,14,22 + 0,0.002,0,0,14,20 + 0,0.01,0,0,14,18 + 0,0.1,0,0,20,25 + 0,0.1,0,0,20,35 + 1,0,0,0,20,36 + 0,0.33,0,0,20,23 + 0,0.37,0,0,20,24 + 0,0.1,0,0,20,22 + 0,0.005,0,0,3,13 + 0,0.9,0,0,3,28 + 1,0,0,0,3,36 + 0,0.005,0,0,3,12 + 0,0.09,0,0,3,17 + 0,0.1698454,0,0,6,8 + 0,0.004995454,0,0,6,13 + 0,0.2497727,0,0,6,28 + 1,0.2058127,0,0,6,36 + 0,0.009960935,0,0,6,29 + 0,0.1998182,0,0,6,9 + 0,0.01998181,0,0,6,21 + 0,0.02997272,0,0,6,10 + 0,0.04995454,0,0,6,26 + 0,0.01498636,0,0,6,27 + 0,0.004995454,0,0,6,11 + 0,0.01998181,0,0,6,16 + 0,0.009960935,0,0,6,14 + 0,0.009960935,0,0,6,19 + 0,0.1,0,0,19,25 + 0,0.1,0,0,19,35 + 1,0,0,0,19,36 + 0,0.33,0,0,19,23 + 0,0.37,0,0,19,24 + 0,0.1,0,0,19,22 + 0,0.156,0,0,2,8 + 0,0.005,0,0,2,13 + 0,0.001,0,0,2,1 + 0,0.26,0,0,2,28 + 1,0.159,0,0,2,36 + 0,0.047,0,0,2,29 + 0,0.104,0,0,2,9 + 0,0.01,0,0,2,21 + 0,0.003,0,0,2,17 + 0,0.104,0,0,2,10 + 0,0.104,0,0,2,26 + 0,0.012,0,0,2,27 + 0,0.01,0,0,2,11 + 0,0.002,0,0,2,5 + 0,0.021,0,0,2,14 + 0,0.001,0,0,2,6 + 0,0.001,0,0,2,4 + 0,0.01,0,0,4,8 + 0,0.1,0,0,4,25 + 1,0,0,0,4,36 + 0,0.319,0,0,4,23 + 0,0.435,0,0,4,24 + 0,0.01,0,0,4,9 + 0,0.005,0,0,4,21 + 0,0.005,0,0,4,31 + 0,0.01,0,0,4,10 + 0,0.1,0,0,4,22 + 0,0.005,0,0,4,20 + 0,0.001,0,0,4,19 + 0,0.05487256,0,0,18,8 + 0,0.3848076,0,0,18,28 + 1,0,0,0,18,36 + 0,0.09875062,0,0,18,23 + 0,0.1649175,0,0,18,29 + 0,0.07686156,0,0,18,24 + 0,0.0109945,0,0,18,9 + 0,0.131934,0,0,18,31 + 0,0.07686156,0,0,18,22 +
+ +
+
+ 0,1,Fleet1,,0,0,1,0 +
+ + -1,1,-1,0,0,15,0,,10,0,0,-9999,0,0,1.051,8,True,True,Adult pollock,0,0,-1,0,-9999,0,0.5,-9999,0,0,,,0,0,112.5,8,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.052,0,,4.5,0,0,-9999,0,0,0.007,13,True,True,Arrowtooth flounder,0,0,-1,0,-9999,0,0.32,-9999,0,0,,,0,0,0.17056,13,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.391,0,,8.312,0,0,-9999,0,0,0.007,1,True,True,Baleen whales,0,0,-1,0,-9999,0,0.02,-9999,0,0,,,0,0,2.592174,1,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,13.812,0,,14,0,0,-9999,0,0,0,25,True,True,Benthic amphipods,0,0,-1,0,-9999,0,2.5,-9999,0,0,,,0,0,81.49081,25,-9999,-1,-1,0,0.4,0,0,-1 + -1,1,-1,0,0,-9999,0,,-9999,0,0,-9999,0,0,-2214.497,35,True,True,Benthic detritus,0,0,-1,0,-9999,0,-9999,-9999,0,0,,,0,0,0,35,-9999,-1,-1,2,0.2,0,-1,-1 + -1,1,-1,0,0,4.83,0,,10.67,0,0,-9999,0,0,0.02,28,True,True,Cephalopods,0,0,-1,0,-9999,0,3.2,-9999,0,0,,,0,0,25.77288,28,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,122.62,0,,26.2,0,0,-9999,0,0,0,32,True,True,Copepods,0,0,-1,0,-9999,0,9.5,-9999,0,0,,,0,0,1405.225,32,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,-9999,0,,-9999,0,0,0,0,0,3464.95,36,True,True,Detrital pool,0,0,-1,0,-9999,0,-9999,-9999,0,0,,,0,0,0,36,-9999,-1,-1,2,0,0,-1,-1 + -1,1,-1,0,0,114.962,0,,5.087,0,0,-9999,0,0,0.004,23,True,True,Epifauna,0,0,-1,0,-9999,0,1.159,-9999,0,0,,,0,0,217.646,23,-9999,-1,-1,0,0.4,0,0,-1 + -1,1,-1,0,0,-9999,0,,3.5,0,0,0.9,0,0,0.001,29,True,True,Forage fish,0,0,-1,0,-9999,0,0.95,-9999,0,0,,,0,0,35.92584,29,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.058,0,,3.6,0,0,-9999,0,0,0.01,12,True,True,Greenland turbot,0,0,-1,0,-9999,0,0.2,-9999,0,0,,,0,0,0.15544,12,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,125.687,0,,12,0,0,-9999,0,0,0,24,True,True,Infauna,0,0,-1,0,-9999,0,1.97,-9999,0,0,,,0,0,657.343,24,-9999,-1,-1,0,0.4,0,0,-1 + -1,1,-1,0,0,1.4,0,,3,0,0,-9999,0,0,0,30,True,True,Jellyfish,0,0,-1,0,-9999,0,1.5,-9999,0,0,,,0,0,1.26,30,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,3.757,0,,13,0,0,-9999,0,0,0,9,True,True,Juvenile pollock,0,0,-1,0,-9999,0,2.5,-9999,0,0,,,0,0,29.6803,9,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.119,0,,2.7,0,0,-9999,0,0,0.008,21,True,True,King crab,0,0,-1,0,-9999,0,0.6,-9999,0,0,,,0,0,0.18564,21,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,120.74,0,,14.457,0,0,-9999,0,0,0,31,True,True,Large Zooplankton,0,0,-1,0,-9999,0,4.399,-9999,0,0,,,0,0,865.2952,31,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,1.156,0,,3.7,0,0,-9999,0,0,0.041,17,True,True,Macrouridae,0,0,-1,0,-9999,0,0.3,-9999,0,0,,,0,0,3.07496,17,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,3.187,0,,3.3,0,0,-9999,0,0,0.226,10,True,True,Pacific cod,0,0,-1,0,-9999,0,0.521,-9999,0,0,,,0,0,6.753253,10,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.787,0,,14.6,0,0,-9999,0,0,0.055,26,True,True,Pacific herring,0,0,-1,0,-9999,0,0.7,-9999,0,0,,,0,0,8.64126,26,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.039,0,,16,0,0,-9999,0,0,0.012,27,True,True,Pacific salmon,0,0,-1,0,-9999,0,4,-9999,0,0,,,0,0,0.3432,27,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,-9999,0,,-9999,0,0,-9999,0,0,-1225.267,34,True,True,Pelagic detritus,0,0,-1,0,-9999,0,-9999,-9999,0,0,,,0,0,0,34,-9999,-1,-1,2,0.2,0,-1,-1 + -1,1,-1,0,0,0.083,0,,3.5,0,0,-9999,0,0,0.01,11,True,True,P.halibut,0,0,-1,0,-9999,0,0.25,-9999,0,0,,,0,0,0.21165,11,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,15,0,,-9999,0,0,-9999,0,0,0,33,True,True,Phytoplankton,0,0,-1,0,-9999,0,234,-9999,0,0,,,0,0,0,33,-9999,-1,-1,1,0,0,0,-1 + -1,1,-1,0,0,0.677,0,,3.5,0,0,-9999,0,0,0.07,16,True,True,Sculpins& Rockfish,0,0,-1,0,-9999,0,0.4,-9999,0,0,,,0,0,1.6248,16,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.01,0,,87,0,0,-9999,0,0,0,7,True,True,Seabirds,0,0,-1,0,-9999,0,0.8,-9999,0,0,,,0,0,0.688,7,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.097,0,,13.34,0,0,-9999,0,0,0.002,5,True,True,Seals,0,0,-1,0,-9999,0,0.06,-9999,0,0,,,0,0,1.029364,5,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,-9999,0,,10.2,0,0,0.9,0,0,0.002,22,True,True,Shrimp,0,0,-1,0,-9999,0,2.04,-9999,0,0,,,0,0,12.88721,22,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.271,0,,4,0,0,-9999,0,0,0.041,15,True,True,Skates,0,0,-1,0,-9999,0,0.4,-9999,0,0,,,0,0,0.7588,15,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.992,0,,6.849,0,0,-9999,0,0,0.041,14,True,True,Small flatfish,0,0,-1,0,-9999,0,0.288,-9999,0,0,,,0,0,5.149671,14,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.249,0,,5,0,0,-9999,0,0,0.006,20,True,True,Snow crab,0,0,-1,0,-9999,0,0.8,-9999,0,0,,,0,0,0.7968,20,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.02,0,,9.8,0,0,-9999,0,0,0,3,True,True,Sperm whales,0,0,-1,0,-9999,0,0.02,-9999,0,0,,,0,0,0.1564,3,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.035,0,,18,0,0,-9999,0,0,0,6,True,True,Steller sea lions,0,0,-1,0,-9999,0,0.06,-9999,0,0,,,0,0,0.5019,6,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.083,0,,5,0,0,-9999,0,0,0,19,True,True,Tanner crab,0,0,-1,0,-9999,0,0.8,-9999,0,0,,,0,0,0.2656,19,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.042,0,,17.5,0,0,-9999,0,0,0,2,True,True,Toothed whales,0,0,-1,0,-9999,0,0.02,-9999,0,0,,,0,0,0.58716,2,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.262,0,,13.9,0,0,-9999,0,0,0.006,4,True,True,Walrus& bearded seals,0,0,-1,0,-9999,0,0.06,-9999,0,0,,,0,0,2.89772,4,-9999,-1,-1,0,0.2,0,0,-1 + -1,1,-1,0,0,0.9,0,,2.52,0,0,-9999,0,0,0.039,18,True,True,Zoarcidae,0,0,-1,0,-9999,0,0.3,-9999,0,0,,,0,0,1.5444,18,-9999,-1,-1,0,0.2,0,0,-1 +
+ +
+
+
+
+ 254200,"Aydin, K.Y.",175,kerim.aydin@noaa.gov,Russian federation,"Created: 18/04/2006 10:01:15 a.m.; 18/04/2006 10:01:35 a.m.; 18/04/2006 10:06:22 a.m.; 18/04/2006 10:06:50 a.m.; 19/04/2006 09:17:49 p.m.; 17/07/2006 08:01:23 p.m.; 17/07/2006 08:04:07 p.m.; 17/07/2006 08:20:59 p.m.; 17/07/2006 08:27:02 p.m.",continental shelf,1981,0,45877.36,,56.07418,191.1522,66.66822,161.5438,1,,Western Bering Sea,3,9,,"Aydin K.Y.,Lapko V.V.,Radchenko V.I.,Livingston P.A.(2002). A comparison of the eastern Bering and western Bering Sea shelf and slope ecosystems through the use of mass-balance food web models NOAA Technical Memorandum",http://www.afsc.noaa.gov/Publications/AFSC-TM/NOAA-TM-AFSC-130.pdf,5,,,1, +
+ +
+
+
+
+ 1200,1 +
+ +
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ 4/5/2023,,Release,6.700014 +
+ \ No newline at end of file diff --git a/man/Rpath-package.Rd b/man/Rpath-package.Rd index 0ef0375e..a6ea23ba 100644 --- a/man/Rpath-package.Rd +++ b/man/Rpath-package.Rd @@ -1,14 +1,12 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Rpath-package.R, R/Rpath_pkg.R +% Please edit documentation in R/Rpath-package.R \docType{package} \name{Rpath-package} +\alias{Rpath} \alias{Rpath-package} -\title{Rpath: A package implementing mass balance algorithms designed to work with +\title{Rpath: A package implementing mass balance algorithms designed to work with fisheries data sources.} \description{ -The Rpath package provides two categories of important functions: -rpath and rsim. - The Rpath package provides two categories of important functions: rpath and rsim. } @@ -19,13 +17,6 @@ system. This is done by solving a series of linear equations for unknown biomas or ecotrophic efficiency. } -\section{rpath functions}{ - -The rpath functions generate the balanced snap shot of energy flow through the -system. This is done by solving a series of linear equations for unknown biomass -or ecotrophic efficiency. -} - \seealso{ Useful links: \itemize{ @@ -33,13 +24,6 @@ Useful links: \item Report bugs at \url{https://github.com/NOAA-EDAB/Rpath/issues} } - -Useful links: -\itemize{ - \item \url{https://noaa-edab.github.io/Rpath/,https://github.com/NOAA-EDAB/Rpath/} - \item Report bugs at \url{https://github.com/NOAA-EDAB/Rpath/issues} -} - } \author{ \strong{Maintainer}: Kerim Aydin \email{kerim.aydin@noaa.gov} (\href{https://orcid.org/0000-0003-3792-9828}{ORCID}) @@ -61,3 +45,4 @@ Other contributors: } } +\keyword{internal} diff --git a/man/adjust.fishing.Rd b/man/adjust.fishing.Rd index 743e557b..93de3b3f 100644 --- a/man/adjust.fishing.Rd +++ b/man/adjust.fishing.Rd @@ -30,11 +30,11 @@ the year are modified.} \item{value}{New value for the parameter.} } \value{ -Returns an \code{Rsim.scenario()} object with the new fishing parameter +Returns an \code{Rsim.scenario()} object with the new fishing parameter values. } \description{ -Modifies the fishing mortality value for a species by a particular gear. +Modifies the fishing mortality value for a species by a particular gear. Parameters that can be adjusted using this function are: \emph{ForcedEffort}, \emph{ForcedFRate}, or \emph{ForcedCatch}. } @@ -44,10 +44,10 @@ Rpath <- rpath(AB.params) # Create a 50 yr Rsim scenario Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) # Change value of forcedFRate for Squids in years 3 through 5 to the value of 2 (for all months) -Rsim.scenario.adjusted.fishing <- adjust.fishing(Rsim.scenario,parameter="ForcedFRate",group="cod",sim.year=3:5,value = 2) -head(Rsim.scenario.adjusted.fishing$fishing$ForcedFRate) - - +Rsim.scenario.adjusted.fishing <- adjust.fishing(Rsim.scenario, parameter = "ForcedFRate", group = "cod", sim.year = 3:5, value = 2) +head(Rsim.scenario.adjusted.fishing$fishing$ForcedFRate) + + } \seealso{ Other Adjust functions: diff --git a/man/adjust.forcing.Rd b/man/adjust.forcing.Rd index 8a76cdea..fc254c01 100644 --- a/man/adjust.forcing.Rd +++ b/man/adjust.forcing.Rd @@ -45,7 +45,7 @@ Rpath <- rpath(AB.params) # Create a 50 yr Rsim scenario Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) # Adjust the ForcedPrey parameter for cod in year 1 for all months. Change the value to 10 -Rsim.scenario.adjusted <- adjust.forcing(Rsim.scenario, parameter="ForcedPrey",group = "cod", sim.year = 1, sim.month=0,value=10) +Rsim.scenario.adjusted <- adjust.forcing(Rsim.scenario, parameter="ForcedPrey",group = "cod", sim.year = 1, sim.month=0,value=10) head(Rsim.scenario.adjusted$forcing$ForcedPrey) diff --git a/man/adjust.scenario.Rd b/man/adjust.scenario.Rd index 31b11a6a..2431e457 100644 --- a/man/adjust.scenario.Rd +++ b/man/adjust.scenario.Rd @@ -10,11 +10,11 @@ adjust.scenario(Rsim.scenario, parameter, group, groupto = NA, value) \item{Rsim.scenario}{Scenario object that contains all of the rsim rates and forcing functions generated by \code{\link{rsim.scenario}()}.} -\item{parameter}{Parameters to be modified (Choose from: \code{B_BaseRef, MzeroMort, +\item{parameter}{Parameters to be modified (Choose from: \code{B_BaseRef, MzeroMort, UnassimRespFrac, ActiveRespFrac, FtimeAdj, FtimeQBOpt, PBopt, NoIntegrate, HandleSelf, ScrambleSelf, QQ, DD, VV, HandleSwitch, PredPredWeight, PreyPreyWeight})} -\item{group}{The model group that the parameter change will affect. Note that +\item{group}{The model group that the parameter change will affect. Note that a value of \emph{'all'} will affect all groups associated with the `groupto` variable. Valid values are found in the `Group` field of the object created from running \code{rpath()}} @@ -29,7 +29,7 @@ variable. Required for parameters \code{QQ}, \code{DD}, \code{VV}, \code{HandleS Returns an \code{rsim.scenario()} object with the new parameter. } \description{ -Modifies the various parameters of the \code{rsim.scenario()} object. Parameters that can be adjusted using this function are: +Modifies the various parameters of the \code{rsim.scenario()} object. Parameters that can be adjusted using this function are: \emph{B_BaseRef}, \emph{MzeroMort},\emph{UnassimRespFrac}, \emph{ActiveRespFrac}, \emph{FtimeAdj}, \emph{FtimeQBOpt}, \emph{PBopt}, \emph{NoIntegrate},\emph{HandleSelf}, \emph{ScrambleSelf}, \emph{QQ}, \emph{DD}, \emph{VV}, \emph{HandleSwitch}, \emph{PredPredWeight}, \emph{PreyPreyWeight} @@ -40,7 +40,7 @@ Rpath <- rpath(AB.params) # Create a 50 yr Rsim scenario Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) # Adjust the PBopt parameter for cod. Set to value = 2 -Rsim.scenario.adjusted <- adjust.scenario(Rsim.scenario, parameter="PBopt",group = "cod", groupto = "all", value = 2) +Rsim.scenario.adjusted <- adjust.scenario(Rsim.scenario, parameter="PBopt",group = "cod", groupto = "all", value = 2) } diff --git a/man/check.rpath.params.Rd b/man/check.rpath.params.Rd index 73550e1f..86946c82 100644 --- a/man/check.rpath.params.Rd +++ b/man/check.rpath.params.Rd @@ -14,7 +14,7 @@ check.rpath.params(Rpath.params) \value{ Checks Rpath parameter files for consistency. An error message will be produced if one of the logical checks fails. Checks include: verification that all types are represented (e.g. consumer, producer, detrital, and fleet); - check if input parameters were entered correctly; check if Diet columns sum to 1. + check if input parameters were entered correctly; check if Diet columns sum to 1. (NOTE: This does not ensure data is correct just that it is in the right places). } \description{ diff --git a/man/create.rpath.from.eiixml.Rd b/man/create.rpath.from.eiixml.Rd new file mode 100644 index 00000000..31a5bb39 --- /dev/null +++ b/man/create.rpath.from.eiixml.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/xml_convert.r +\name{create.rpath.from.eiixml} +\alias{create.rpath.from.eiixml} +\title{Creates an Rpath object from an EwE exported model (EIIXML format)} +\usage{ +create.rpath.from.eiixml(eiifile, verbose = FALSE) +} +\arguments{ +\item{eiifile}{Path to exported EwE XML file} + +\item{verbose}{Logical. Use for debugging. If TRUE, prints out a list of parsed tables and number of rows read from the XML file} +} +\value{ +An Rpath object (list) with the following components: +\item{stanzas}{Stanza parameters} +\item{pedigree}{Pedigree parameters} +\item{diet}{Diet parameters} +\item{model}{Model parameters} +} +\description{ +Parses an eiixml file exported using the Ecopath with Ecosim (EwE) GUI +into an unbalanced Rpath model object, which can then be balanced. +This function was tested on files exported by EwE version 6.7. +} +\examples{ +# Import an eiixml file previously exported from the EwE GUI, check that the +# unbalanced model is functional, recalculate age-structured groups based on +# imported stanza parameters (if the imported model has stanzas), and balance +# the model. +eiixml <- system.file("extdata/xml/Western_Bering_Sea.eiixml", package = "Rpath") +rpath_object <- create.rpath.from.eiixml(eiixml) +check.rpath.params(rpath_object) +rpath_object <- rpath.stanzas(rpath_object) +rpath_balanced <- rpath(rpath_object) + +} +\seealso{ +Other xml: +\code{\link{import.eiixml}()} +} +\concept{xml} diff --git a/man/create.rpath.params.Rd b/man/create.rpath.params.Rd index 05960910..77615069 100644 --- a/man/create.rpath.params.Rd +++ b/man/create.rpath.params.Rd @@ -17,7 +17,7 @@ Fleet = 3.} \value{ Outputs a list object of \code{Rpath.params} which are populated with values of NA or logical default values. Values can then be filled in using - R. Use \code{\link{check.rpath.params}()} to ensure parameter files are filled out + R. Use \code{check.rpath.params()} to ensure parameter files are filled out correctly (NOTE: This does not ensure data is correct just that it is in the right places). } diff --git a/man/frate.table.Rd b/man/frate.table.Rd index 62b06f5f..cf54087f 100644 --- a/man/frate.table.Rd +++ b/man/frate.table.Rd @@ -17,7 +17,7 @@ Returns a data table of F values for each species/gear combination. \item{Total}{Sum over all gear types} } \description{ -Creates a table of fishing mortalities by species group and gear for an +Creates a table of fishing mortalities by species group and gear for an \code{rsim.scenario()} object. } \examples{ diff --git a/man/get.rsim.fishing.Rd b/man/get.rsim.fishing.Rd index c355dd60..1d1c2694 100644 --- a/man/get.rsim.fishing.Rd +++ b/man/get.rsim.fishing.Rd @@ -14,7 +14,7 @@ forcing functions generated by \code{\link{rsim.scenario}()}.} Returns a `fishing` object. } \description{ -Helper function that will retrieve the fishing forcing parameters that were used in an +Helper function that will retrieve the fishing forcing parameters that were used in an Rsim scenario } \examples{ diff --git a/man/get.rsim.forcing.Rd b/man/get.rsim.forcing.Rd index fc0873e1..cd4d5c64 100644 --- a/man/get.rsim.forcing.Rd +++ b/man/get.rsim.forcing.Rd @@ -14,7 +14,7 @@ forcing functions generated by \code{\link{rsim.scenario}()}.} Returns a `forcing` object. } \description{ -Helper function that will retrieve the forcing parameters that were used in an +Helper function that will retrieve the forcing parameters that were used in an Rsim scenario } \examples{ diff --git a/man/get.rsim.params.Rd b/man/get.rsim.params.Rd index 366e9a83..ce94b39e 100644 --- a/man/get.rsim.params.Rd +++ b/man/get.rsim.params.Rd @@ -14,7 +14,7 @@ forcing functions generated by \code{\link{rsim.scenario}()}.} Returns a `params` object. } \description{ -Helper function that will retrieve the parameters that were used +Helper function that will retrieve the parameters that were used in an Rsim scenario } \examples{ diff --git a/man/get.rsim.stanzas.Rd b/man/get.rsim.stanzas.Rd index 9c52b75b..1719eabf 100644 --- a/man/get.rsim.stanzas.Rd +++ b/man/get.rsim.stanzas.Rd @@ -14,7 +14,7 @@ forcing functions generated by \code{\link{rsim.scenario}()}.} Returns a `stanzas` object. } \description{ -Helper function that will retrieve the stanza parameters that were used in an +Helper function that will retrieve the stanza parameters that were used in an Rsim scenario } \examples{ diff --git a/man/get.rsim.start_state.Rd b/man/get.rsim.start_state.Rd index 99c1b093..0eaf945c 100644 --- a/man/get.rsim.start_state.Rd +++ b/man/get.rsim.start_state.Rd @@ -14,7 +14,7 @@ forcing functions generated by \code{\link{rsim.scenario}()}.} Returns an `start_state` object. } \description{ -Helper function that will retrieve the starting state values that were used +Helper function that will retrieve the starting state values that were used in an Rsim scenario } \examples{ @@ -23,7 +23,7 @@ Rpath <- rpath(AB.params) # Create a 50 yr Rsim scenario Rsim.scenario <- rsim.scenario(Rpath, AB.params, years = 1:50) params <- get.rsim.start_state(Rsim.scenario) -names(params) +names(params) diff --git a/man/import.eiixml.Rd b/man/import.eiixml.Rd new file mode 100644 index 00000000..b25f40d2 --- /dev/null +++ b/man/import.eiixml.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/xml_convert.r +\name{import.eiixml} +\alias{import.eiixml} +\title{Reads in EwE exported XML file and parses into data frames} +\usage{ +import.eiixml(eiifile, verbose = FALSE) +} +\arguments{ +\item{eiifile}{Path to exported EwE XML file} + +\item{verbose}{Logical. Use for debugging. If TRUE, prints out a list of parsed tables and number of rows read from the XML file} +} +\value{ +A list of data frames, one data frame for each node (exported EwE table) in the XML file. Each table +has the naming convention ewe_[table name] where [table name] is the name of the table provided by EwE. +} +\description{ +Parses an eiixml file exported using the Ecopath with Ecosim (EwE) GUI into a list of data frames, +one frame for each table in the exported XML file. This function is usually called by the function \code{create.rpath.from.eiixml()} that +converts these tables into an unbalanced rpath model object. However import.xml can be used on its own to examine the full set +of tables exported by EwE, including tables not currently imported into Rpath objects, such as Ecosim runs or model metadata. +This function was tested on files exported by EwE version 6.7. +} +\examples{ +# Import an eiixml file previously exported from the EwE GUI into a list of +# data frames containing the model data +eiixml <- system.file("extdata/xml","Western_Bering_Sea.eiixml", package = "Rpath") +xml_data <- import.eiixml(eiixml) + + +} +\seealso{ +Other xml: +\code{\link{create.rpath.from.eiixml}()} +} +\concept{xml} diff --git a/man/read.rpath.params.Rd b/man/read.rpath.params.Rd index 363af02f..88d6021f 100644 --- a/man/read.rpath.params.Rd +++ b/man/read.rpath.params.Rd @@ -29,7 +29,7 @@ be created.} \value{ Outputs an \code{Rpath.param} object that can be used for Rpath and subsequently Rsim. (NOTE: This does function does not ensure data is correct or in the - correct locations.Run \code{\link{check.rpath.param}()} to ensure the appropriate columns are + correct locations.Run \code{check.rpath.params()} to ensure the appropriate columns are present). } \description{ diff --git a/man/rpath.Rd b/man/rpath.Rd index c5a209f2..b9d55f82 100644 --- a/man/rpath.Rd +++ b/man/rpath.Rd @@ -22,7 +22,7 @@ Returns a static \code{Rpath} model that can be supplied to the \code{rsim.scenario} function. } \description{ -Performs initial mass balance using a \code{\link{Rpath.params}()} file +Performs initial mass balance using a \code{Rpath.params()} file } \seealso{ Other Rpath functions: diff --git a/man/rpath.stanzas.Rd b/man/rpath.stanzas.Rd index fc4fbfa2..82af7219 100644 --- a/man/rpath.stanzas.Rd +++ b/man/rpath.stanzas.Rd @@ -12,7 +12,7 @@ rpath.stanzas(Rpath.params) or \code{\link{read.rpath.params}()} functions.} } \value{ -Calculates and adds biomass and consumption for trailing stanza groups. +Calculates and adds biomass and consumption for trailing stanza groups. Also adds weight at age and number at age for multi-staza groups. } \description{ diff --git a/man/set.rsim.scene.Rd b/man/set.rsim.scene.Rd index 0d31cbf7..3af496e4 100644 --- a/man/set.rsim.scene.Rd +++ b/man/set.rsim.scene.Rd @@ -31,7 +31,7 @@ forcing functions generated by \code{\link{rsim.scenario}()}.} Returns an \code{Rsim.scenario} object with the new parameter. } \description{ -Modifies the various parameters of the \code{rsim.scenario()} object. Parameters +Modifies the various parameters of the \code{rsim.scenario()} object. Parameters that can be adjusted using this function are: \code{params},\code{start_state}, \code{forcing},\code{fishing},\code{stanzas} } diff --git a/vignettes/Convert_EwE_to_Rpath.Rmd b/vignettes/Convert_EwE_to_Rpath.Rmd new file mode 100644 index 00000000..66da9efb --- /dev/null +++ b/vignettes/Convert_EwE_to_Rpath.Rmd @@ -0,0 +1,136 @@ +--- +title: "Convert EwE to Rpath" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Convert EwE to Rpath} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(Rpath) +``` + +If you have developed your model in the Ecopath with Ecosim (EwE) software package and you now want to bring it into Rpath for analysis, there are tools available in `Rpath` to do just that. First, you need to export your model from EwE to XML format, using EwE's `File -> Export Model -> To XML` menu option to create a .eiixml file (all EwE functionality described here was tested using EwE version 6.7.0). + +For this example, the Western Bering Sea EwE model available for download on EcoBase (Aydin et al. 2002) was first imported into EwE, using EwE's Import Model menu option, then exported to the file `Western_Bering_Sea.eiixml`. We can read in this file and convert it to an `Rpath` object using the function `create.rpath.from.eiixml()`. + +```{r convert1, echo = TRUE, eval = TRUE} +# Read in the model +eiixml_file <- system.file("extdata/xml", "Western_Bering_Sea.eiixml", package = "Rpath") +model <- create.rpath.from.eiixml(eiixml_file) +``` + +Currently, the `create.rpath.from.eiixml()` will only import input parameters related to the unbalanced +Ecopath model, and will not import Ecosim information or other supplementary tables. However, the full +set of tables in an .eiixml file can be examined (or manually imported) using the `import.eiixml()` +helper function. + +The variable `model` is a list object (an unbalanced Ecopath model) containing the data objects required to balance your model and run simulations. There are 4 objects that comprise an `Rpath` object: + +* model - including basic inputs, detritus fate, other production, and fishery landings/discards. +* diet - a diet matrix, including diet import. +* stanzas - parameters for multi-stanza functional groups. +* pedigree - data pedigree (quality) for the model. + +Prepare your imported model for balancing by first calculating +multistanza parameters that are based on age structure (such as B and Q/B), using the +`rpath.stanzas()` function, then using the `check.rpath.params()` function to ensure all of the +model's parameters are entered (If your model has no multistanza groups, `rpath.stanzas()` +can be called but will not change the model). The `check.rpath.params()` function will +make sure the data was read in correctly, and produce warnings in situations that might +prevent the model from balancing; for example, if a functional group is missing too many +parameters to be balanced, or if diet compositions for a predator do not sum to 1. + +```{r check, echo = TRUE, eval = TRUE} +model <- rpath.stanzas(model) +check.rpath.params(model) +``` + +## Model parameters + +To compare Rpath balance results with EwE, first balance the model in EwE and use +EwE's "Save to a csv" icon on EwE's Outputs -> Basic estimates tab (upper right +corner). It is recommended that you first increase the precision (digits displayed) +in EwE to 8 digits. The following script, using the dplyr library, can then compare the +resulting balances. + +```{r compare, echo = TRUE, eval = TRUE} + +# dplyr is used to match the Rpath and EwE balance results + library(dplyr, warn.conflicts=F) + +# EwE output files - eiixml file and Basic Estimates csv + eiifile <- system.file("extdata/xml", "Western_Bering_Sea.eiixml", package = "Rpath") + csvfile <- system.file("extdata/xml", "Western Bering Sea-Basic estimates.csv", package = "Rpath") + +# Load and balance the model in Rpath + unbal <- create.rpath.from.eiixml(eiifile) + unbal <- rpath.stanzas(unbal) + check.rpath.params(unbal) + bal <- rpath(unbal) + +# Load the csv file and clean up the format + csv.out <- read.csv(csvfile) + # If there's no stanzas, the total mortality column will be missing from the CSV - add an extra column + if(length(csv.out)<13){csv.out <- data.frame(append(csv.out, list(tm=NA), after=6))} + # cleaned column names matching EwE output + names(csv.out)<- c("X", "Group.name", "Trophic.level", "Hab.area", "Biomass.in.habitat.area", + "Biomass", "Total.mortality", "Production.biomass", "Consumption.biomass", + "Ecotrophic.Efficiency", "Production.consumption", "Biomass.accumulation", + "BA.rate") + # Drop placeholder lines that head stanza groups + csv.out <- csv.out[!is.na(csv.out$X),] + # Clean names to rpath standard + csv.out$group <- janitor::make_clean_names(csv.out$Group.name) + +# Create a data frame of rpath (balanced model) output + rpath.dat <- data.frame(group=bal$Group,type=bal$type, tl=bal$TL, biomass=bal$Biomass, + pb=bal$PB,qb=bal$QB, EE=bal$EE, pc=bal$GE, + ba=bal$BA)[c(rpath.living(bal), rpath.detrital(bal)),] + +# Join rpath and EwE balance outputs, using dplyr + rpath_ewe_table <- rpath.dat %>% + dplyr::left_join(csv.out, by="group") %>% + dplyr::mutate( biomass_test = abs(biomass-Biomass)/Biomass, + ee_test = abs(EE-Ecotrophic.Efficiency)) %>% + dplyr::select(group, rpath.biomass=biomass, ewe.biomass=Biomass, prop.bio.diff=biomass_test, + rpath.EE=EE, ewe.EE=Ecotrophic.Efficiency, EE.diff=ee_test) + + rpath_ewe_table + +``` + +## Differences between Rpath and EwE balances + +Rpath was tested using the 150+ models available on EcoBase. In general, differences +between Rpath and EwE were within 1.0e-7 of each other for any given estimated EE, B, +or P/B, with differences due to +numerical precision of inputs. Stanza groups may show differences of up to 1.0e-3 +due to accumulated differences summing across age groups. However, there are +some specific differences in functionality in some models: + +* Q/B - Rpath does not currently estimate Q/B as part of the main balancing (i.e. to +estimate Q/B, both P/B and PC must be supplied). +* Interdetrital flows - For models with three or more stages of interdetrital +flows (that is, if detritus fate specifies thart Detritus A -> Detritus B -> Detritus C), +results may differ between Rpath and EwE. Two-stage detrital flows (Detritus A -> Detritus B) +give matching results. +* Rpath will not currently estimate BA. In cases where all of B, EE, and P/B are +supplied, Rpath will not estimate BA or make changes to the model, and may produce +inconsistent results if one of those paramters isn't set to NA. +* Rpath does not currently support BA rate, only total BA. If your model has a BA rate, +it is recommended that you calculate the total BA in EwE (as an output), and enter +that value into Rpath. + + + + + diff --git a/vignettes/ModelSetup.Rmd b/vignettes/ModelSetup.Rmd index 60b65389..bffe506c 100644 --- a/vignettes/ModelSetup.Rmd +++ b/vignettes/ModelSetup.Rmd @@ -33,11 +33,14 @@ Unlike the GUI based EwE software package, Rpath relies on a parameter input fil This file is actually a list of several different parameter files: model, diet, stanzas, and pedigree. Parameter files can be created outside of R and read in using the `read.rpath.params()` function. This function will merge several different flat -files into an R object of the list type. A preferred alternative is to generate +files into an R object of the list type. An alternative is to generate the list file and populate it completely within R. The function `create.rpath.params()` will generate an `Rpath.params` list object. This ensures that all of the correct columns are present in the parameter file. +Additionally, an existing model can be imported from EwE, by first exporting the +model from EwE in .eiixml format via the GUI interface, and using the `create.rpath.from.eiixml` function to read in the exported file (see the vignette `Convert_EwE_to_Rpath`). + The parameter file contains all of the information you would normally enter in the input data tabs in EwE. There are 2 necessary pieces of information to generate the parameter file: the group names and their corresponding type. The types are: