diff --git a/NEWS.md b/NEWS.md index 79e72ed12..73b655072 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,7 @@ # development version +- Adds a warning to `project()` and `projectToSteady()` / `steady()` if the abundance of any species at its maximum size `w_max` is not negligible, indicating that the chosen `w_max` is too small and individuals are being cut off by the upper boundary condition. + - New experimental `steadyNewton()` finds a steady state by solving the steady-state equation directly with a Newton-type root finder (using the `nleqslv` package) instead of running the dynamics to convergence. Unlike diff --git a/R/helpers.R b/R/helpers.R index b0c4000ab..3fa7bca6b 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -302,3 +302,52 @@ get_steady_state_n <- function(params, g, mu, D, N0, return(zero_above_support(n, w_top)) } + +#' Warn if w_max is too small for any species +#' +#' Checks whether the total flux (advective + diffusive) of any species at its +#' maximum size `w_max` is non-negligible compared to the maximum total flux +#' across all sizes. +#' +#' It also issues a warning if species do not even reach maturity size. +#' +#' @param params A MizerParams object. +#' @param n An array (species x size) with the number density. +#' @param rates A list of rates as returned by `mizerRates()`. +#' @param threshold The threshold ratio above which a warning is issued. +#' @noRd +warn_if_w_max_violation <- function(params, n, rates, threshold = 1e-5) { + w_top <- support_top_idx(params) + no_sp <- nrow(params@species_params) + flux <- getFlux(params, n = n) + violating_species <- character() + starving_species <- character() + for (i in seq_len(no_sp)) { + w_mat_idx <- sum(params@w <= params@species_params[i, "w_mat"]) + w_max_idx <- sum(params@w <= params@species_params[i, "w_max"]) + # increase by 1 is needed in second-order scheme and it is also + # needed in first-order scheme to not warn on existing models + w_max_idx <- min(w_max_idx + 1, length(params@w)) + if (flux[i, w_mat_idx] > 0) { + ratio <- flux[i, w_top[i]] / flux[i, w_mat_idx] + if (ratio > threshold) { + violating_species <- c(violating_species, params@species_params$species[i]) + } + } else { + starving_species <- c(starving_species, params@species_params$species[i]) + } + } + if (length(violating_species) > 0) { + warning("For the following species, the abundance is not negligible at their maximum size w_max:\n", + paste(violating_species, collapse = ", "), + ".\nThis means that individuals are being cut off by the upper boundary condition. ", + "You should choose a larger w_max for these species.", + call. = FALSE) + } + if (length(starving_species) > 0) { + warning("The following species never reach their maturity size:\n", + paste(violating_species, collapse = ", "), + call. = FALSE) + } +} + diff --git a/R/project.R b/R/project.R index d4cfeb0f8..bc8d8a3b3 100644 --- a/R/project.R +++ b/R/project.R @@ -391,6 +391,8 @@ project.MizerParams <- function(object, effort, sim@n_other[i, ] <- unserialize(serialize(n_list$n_other, NULL)) } + warn_if_w_max_violation(params, n = n_list$n, rates = n_list$rates) + return(sim) } diff --git a/R/steady.R b/R/steady.R index fff46c24b..a77f41423 100644 --- a/R/steady.R +++ b/R/steady.R @@ -218,6 +218,8 @@ projectToSteady.MizerParams <- function(params, params@initial_n_pp[] <- current$n_pp params@initial_n_other[] <- current$n_other + warn_if_w_max_violation(params, n = current$n, rates = current$rates) + if (return_sim) { sim@params <- params sel <- 1:i diff --git a/tests/testthat/test-project.R b/tests/testthat/test-project.R index 4cf799b0a..ded0247f4 100644 --- a/tests/testthat/test-project.R +++ b/tests/testthat/test-project.R @@ -498,3 +498,15 @@ test_that("Predation diffusion changes simulation trajectory", { sim_base <- project(NS_params_small, t_max = 1) expect_false(identical(sim_d@n, sim_base@n)) }) + +test_that("warn_if_w_max_violation warns when w_max is too small", { + # Check that it doesn't warn on default NS_params_small + expect_warning(project(NS_params_small, t_max = 1), NA) + + # Check that it warns when w_max of Herring is set too small + params <- NS_params_small + params@species_params$w_max[2] <- 5 + params <- validParams(params) + expect_warning(project(params, t_max = 1), "abundance is not negligible at their maximum size w_max") +}) + diff --git a/tests/testthat/test-steadyNewton.R b/tests/testthat/test-steadyNewton.R index dd85d6117..395edb2ad 100644 --- a/tests/testthat/test-steadyNewton.R +++ b/tests/testthat/test-steadyNewton.R @@ -100,7 +100,8 @@ test_that("steadyNewton works with the second-order (van Leer) scheme", { sow$flux <- "van_leer" sow$bin_average <- TRUE second_order_w(p) <- sow - ps <- steady(p, t_max = 100, progress_bar = FALSE) + expect_warning(ps <- steady(p, t_max = 100, progress_bar = FALSE), + "abundance is not negligible") pn <- steadyNewton(ps) expect_s4_class(pn, "MizerParams") @@ -108,7 +109,8 @@ test_that("steadyNewton works with the second-order (van Leer) scheme", { # dynamics. The van Leer limiter is only Lipschitz, so the Newton residual # cannot be driven to machine precision, but the projected drift is the # honest test and stays small. - sim <- project(pn, t_max = 1, dt = 0.25, t_save = 1) + expect_warning(sim <- project(pn, t_max = 1, dt = 0.25, t_save = 1), + "abundance is not negligible") n0 <- pn@initial_n n1 <- finalN(sim) support <- n0 > 0