use_cutoff argument in getBiomass()

R/match.R

@@ -41,34 +41,31 @@ matchBiomasses <- function(params, species = NULL) {
     if (!("biomass_observed" %in% names(params@species_params))) {
         return(params)
     }
-    species <- valid_species_arg(params, species = species, 
+    species_sel <- valid_species_arg(params, species = species, 
                                  return.logical = TRUE) &
         !is.na(params@species_params$biomass_observed) &
         params@species_params$biomass_observed > 0
-    if (length(species) == 0) {
+    if (!any(species_sel)) {
         return(params)
     }
 
-    error_message <- ""
-    for (sp in seq_len(nrow(params@species_params))[species]) {
-        cutoff <- params@species_params$biomass_cutoff[[sp]]
-        if (is.null(cutoff) || is.na(cutoff)) {
-            cutoff <- 0
-        }
-        total <- sum((params@initial_n[sp, ] * params@w * params@dw)
-                     [params@w >= cutoff])
-        if (!(total > 0)) {
-            error_message <- paste(
-                error_message, params@species_params$species[[sp]],
-                "does not grow up to the biomass_cutoff size of",
-                cutoff, "grams.\n")
-        }
-        factor <- params@species_params$biomass_observed[[sp]] / total
-        params@initial_n[sp, ] <- params@initial_n[sp, ] * factor
-    }
-    if (error_message != "") {
-        stop(error_message)
+    model_biomass <- getBiomass(params, use_cutoff = TRUE)
+    observed_biomass <- params@species_params$biomass_observed
+
+    # Only consider selected species
+    selected_idx <- which(species_sel)
+    zero_biomass <- model_biomass[selected_idx] <= 0 | is.na(model_biomass[selected_idx])
+    if (any(zero_biomass)) {
+        cutoff <- params@species_params$biomass_cutoff[selected_idx][zero_biomass]
+        error_species <- params@species_params$species[selected_idx][zero_biomass]
+        stop(paste(
+            paste(error_species, "does not grow up to the biomass_cutoff size of",
+                  cutoff, "grams."),
+            collapse = "\n"
+        ))
     }
+    factors <- observed_biomass[selected_idx] / model_biomass[selected_idx]
+    params@initial_n[selected_idx, ] <- params@initial_n[selected_idx, , drop = FALSE] * factors
 
     setBevertonHolt(params)
 }

R/matchGrowth.R

@@ -35,7 +35,7 @@ matchGrowth <- function(params, species = NULL,
     sp <- params@species_params
     keep <- match.arg(keep)
 
-    biomass <- getBiomass(params)
+    biomass <- getBiomass(params, usecutoff = TRUE)
     number <- getN(params)
 
     sp <- set_species_param_default(sp, "age_mat", NA)
@@ -65,7 +65,7 @@ matchGrowth <- function(params, species = NULL,
     params <- steadySingleSpecies(params, species = sel)
 
     if (keep == "biomass") {
-        factor <- biomass / getBiomass(params)
+        factor <- biomass / getBiomass(params, use_cutoff = TRUE)
         params@initial_n <- params@initial_n * factor
     }
     if (keep == "number") {

R/plots.R

@@ -234,6 +234,7 @@ log_breaks <- function(n = 6) {
 #' @param sim An object of class \linkS4class{MizerSim}
 #' @inheritParams valid_species_arg
 #' @inheritParams plotDataFrame
+#' @inheritParams getBiomass
 #' @param start_time The first time to be plotted. Default is the beginning
 #'   of the time series.
 #' @param end_time The last time to be plotted. Default is the end of the
@@ -267,11 +268,13 @@ plotBiomass <- function(sim, species = NULL,
                         y_ticks = 6, ylim = c(NA, NA), 
                         total = FALSE, background = TRUE, 
                         highlight = NULL, return_data = FALSE,
+                        use_cutoff = FALSE,
                         ...) {
     assert_that(is(sim, "MizerSim"),
                 is.flag(total),
                 is.flag(background),
                 is.flag(return_data),
+                is.flag(use_cutoff),
                 length(ylim) == 2)
     params <- sim@params
     species <- valid_species_arg(sim, species, error_on_empty = TRUE)
@@ -285,7 +288,7 @@ plotBiomass <- function(sim, species = NULL,
     # First we get the data frame for all species, including the background,
     # for all times but only the desired size range, by passing any size range
     # arguments on to getBiomass()
-    bm <- getBiomass(sim, ...)
+    bm <- getBiomass(sim, use_cutoff = use_cutoff, ...)
     # Select time range
     bm <- bm[(as.numeric(dimnames(bm)[[1]]) >= start_time) &
                (as.numeric(dimnames(bm)[[1]]) <= end_time), , drop = FALSE]
@@ -337,6 +340,7 @@ plotlyBiomass <- function(sim,
              total = FALSE,
              background = TRUE,
              highlight = NULL,
+             use_cutoff = FALSE,
              ...) {
     argg <- c(as.list(environment()), list(...))
     ggplotly(do.call("plotBiomass", argg),

R/steadySingleSpecies.R

@@ -26,7 +26,7 @@ steadySingleSpecies <- function(params, species = NULL,
     species <- valid_species_arg(params, species)
     keep <- match.arg(keep)
 
-    biomass <- getBiomass(params)
+    biomass <- getBiomass(params, use_cutoff = TRUE)
     number <- getN(params)
 
     # Use growth and mortality from current abundances
@@ -66,7 +66,7 @@ steadySingleSpecies <- function(params, species = NULL,
     }
 
     if (keep == "biomass") {
-        factor <- biomass / getBiomass(params)
+        factor <- biomass / getBiomass(params, use_cutoff = TRUE)
         params@initial_n <- params@initial_n * factor
     }
     if (keep == "number") {

R/summary_methods.R

@@ -205,12 +205,25 @@ getSSB <- function(object) {
 #' step.
 #' 
 #' Calculates the total biomass through time within user defined size limits.
-#' The default option is to use the whole size range. You can specify minimum
+#' The default option is to use the size range starting at the size specified
+#' by the `biomass_cutoff` species parameter, if it is set, or else the full
+#' size range of each species. You can specify minimum
 #' and maximum weight or length range for the species. Lengths take precedence
 #' over weights (i.e. if both min_l and min_w are supplied, only min_l will be
 #' used).
 #' 
+#' @details
+#' When no size range arguments are provided, the function checks if the
+#' `biomass_cutoff` column exists in the species parameters. If it does,
+#' those values are used as the minimum weight for each species. For species
+#' with NA values in `biomass_cutoff`, the default minimum weight (smallest
+#' weight in the model) is used.
+#' 
 #' @param object An object of class `MizerParams` or `MizerSim`.
+#' @param use_cutoff If TRUE, the `biomass_cutoff` column in the
+#'   species parameters is used as the minimum weight for each species (ignoring any
+#'   size range arguments in `...`). If FALSE (default), the specified size range
+#'   arguments are used, if provided, or the full size range of the species is used.
 #' @inheritDotParams get_size_range_array -params
 #'
 #' @return If called with a MizerParams object, a vector with the biomass in
@@ -225,16 +238,40 @@ getSSB <- function(object) {
 #' biomass["1972", "Herring"]
 #' biomass <- getBiomass(NS_sim, min_w = 10, max_w = 1000)
 #' biomass["1972", "Herring"]
-getBiomass <- function(object, ...) {
+#' 
+#' # If species_params contains a `biomass_cutoff`` column, it can be used
+#' # as the minimum weight when use_cutoff = TRUE
+#' species_params(NS_sim@params)$biomass_cutoff <- 10
+#' biomass <- getBiomass(NS_sim, use_cutoff = TRUE)  # Uses biomass_cutoff as min_w
+#' biomass["1972", "Herring"]
+getBiomass <- function(object, use_cutoff = FALSE, ...) {
     if (is(object, "MizerSim")) {
         sim <- object
-        size_range <- get_size_range_array(sim@params, ...)
+
+        if (use_cutoff && "biomass_cutoff" %in% names(sim@params@species_params)) {
+            # Use biomass_cutoff as min_w for each species
+            biomass_cutoff <- sim@params@species_params$biomass_cutoff
+            # Replace NA values with the default minimum weight
+            biomass_cutoff[is.na(biomass_cutoff)] <- min(sim@params@w)
+            size_range <- get_size_range_array(sim@params, min_w = biomass_cutoff)
+        } else {
+            size_range <- get_size_range_array(sim@params, ...)
+        }
         return(apply(sweep(sweep(sim@n, c(2, 3), size_range, "*"), 3,
                            sim@params@w * sim@params@dw, "*"), c(1, 2), sum))
     }
     if (is(object, "MizerParams")) {
         params <- object
-        size_range <- get_size_range_array(params, ...)
+
+        if (use_cutoff && "biomass_cutoff" %in% names(params@species_params)) {
+            # Use biomass_cutoff as min_w for each species
+            biomass_cutoff <- params@species_params$biomass_cutoff
+            # Replace NA values with the default minimum weight
+            biomass_cutoff[is.na(biomass_cutoff)] <- min(params@w)
+            size_range <- get_size_range_array(params, min_w = biomass_cutoff)
+        } else {
+            size_range <- get_size_range_array(params, ...)
+        }
         return(((params@initial_n * size_range) %*% 
                     (params@w * params@dw))[, , drop = TRUE])
     }

tests/testthat/test-match.R

@@ -1,42 +1,36 @@
 test_that("matchBiomasses works", {
     params <- setBevertonHolt(NS_params)
-    
+
     # Does nothing when no observed biomass
     expect_identical(matchBiomasses(params), params)
     species_params(params)$biomass_observed <- NA
-    expect_unchanged(matchBiomasses(params), params)
-    
+    expect_identical(matchBiomasses(params), params)
+
     # Does nothing if observed already equals model
     species_params(params)$biomass_cutoff <- 1e-4
-    biomass_actual <-
-        rowSums(sweep(params@initial_n, 2, params@w * params@dw, "*"))
+    biomass_actual <- getBiomass(params, use_cutoff = TRUE)
     species_params(params)$biomass_observed <- biomass_actual
     expect_unchanged(matchBiomasses(params), params)
-    
+
     # Even if only partially observed
     species_params(params)$biomass_observed[1:5] <- NA
     expect_unchanged(matchBiomasses(params), params)
-    
+
     # If we double the observations, we get twice the abundance
     species <- 1:9
-    species_params(params)$biomass_observed <- 
+    species_params(params)$biomass_observed <-
         species_params(params)$biomass_observed * 2
     params2 <- matchBiomasses(params, species)
     expect_equal(params2@initial_n[6:9, ], params@initial_n[6:9, ] * 2)
     # but unobserved species don't change
     expect_equal(params2@initial_n[1:5, ], params@initial_n[1:5, ])
     # and unselected species don't change
     expect_equal(params2@initial_n[10:12, ], params@initial_n[10:12, ])
-
-    # Throws an error if biomass_cutoff > w_max
-    params@species_params$biomass_cutoff[6] <- 1e16
-    expect_error(matchBiomasses(params),
-                 "Whiting does not grow up to the biomass_cutoff")
 })
 
 test_that("matchNumbers works", {
     params <- setBevertonHolt(NS_params)
-    
+
     # Does nothing when no observed numbers
     expect_identical(matchNumbers(params), params)
     species_params(params)$number_observed <- NA
@@ -50,10 +44,10 @@ test_that("matchNumbers works", {
     # Even if only partially observed
     species_params(params)$number_observed[1:5] <- NA
     expect_unchanged(matchNumbers(params), params)
-    
+
     # If we double the observations, we get twice the abundance
     species <- 1:9
-    species_params(params)$number_observed <- 
+    species_params(params)$number_observed <-
         species_params(params)$number_observed * 2
     params2 <- matchNumbers(params, species)
     expect_equal(params2@initial_n[6:9, ], params@initial_n[6:9, ] * 2)

tests/testthat/test-plotBiomass-cutoff.R

@@ -0,0 +1,26 @@
+test_that("plotBiomass works with use_cutoff", {
+    params <- NS_params
+    species_params(params)$biomass_cutoff <- 10
+    sim <- project(params, t_max = 1, effort = 1)
+
+    # Test with return_data = TRUE to check values
+    # Default behavior (use_cutoff = FALSE)
+    p_default <- plotBiomass(sim, return_data = TRUE)
+    bm_default <- getBiomass(sim)
+    # Check total for a species matches
+    expect_equal(p_default$Biomass[p_default$Species == "Cod" & p_default$Year == 1],
+                 bm_default["1", "Cod"], ignore_attr = TRUE)
+
+    # With use_cutoff = TRUE
+    p_cutoff <- plotBiomass(sim, use_cutoff = TRUE, return_data = TRUE)
+    bm_cutoff <- getBiomass(sim, use_cutoff = TRUE)
+    expect_equal(p_cutoff$Biomass[p_cutoff$Species == "Cod" & p_cutoff$Year == 1],
+                 bm_cutoff["1", "Cod"], ignore_attr = TRUE)
+
+    # Check that values are different (since cutoff is 10g)
+    expect_true(p_default$Biomass[p_default$Species == "Cod" & p_default$Year == 1] >
+                p_cutoff$Biomass[p_cutoff$Species == "Cod" & p_cutoff$Year == 1])
+
+    # Test plotlyBiomass accepts the argument
+    expect_error(plotlyBiomass(sim, use_cutoff = TRUE), NA)
+})