From f53afac426e65a1388fbe65d83ef3c01d42a5f06 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Sat, 15 Aug 2026 23:50:32 +0100 Subject: [PATCH] Give the resource a weight-length relationship MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `plotSpectra(params, size_axis = "l")` used to drop the resource silently, because the length conversion needs an allometric weight-length relationship and only the species had one. `resource_params()` now carries `a` and `b`, defaulting to the equivalent spherical diameter of an organism with the density of water, `a = pi/6` and `b = 3`. That is the convention plankton ecology uses for a composite of many taxa, and on a mizer grid it puts the smallest resource sizes below a micrometre and a milligram organism at about a millimetre. The resource therefore appears on length-based spectra, and the resource array classes gain `size_axis` and `llim`, `animate()` included. The parameters feed none of the rates, so they are not arguments of `setResource()`; it fills them in as defaults so that `resource_params()` shows them, and they are also defaulted at the point of use so that models saved before they existed need no upgrade. Note that the resource then sits on the length axis at its own convention: a fish of a given weight is about 3.7 times longer than a sphere of that weight. That difference is real rather than an artefact — a 1 mg copepod really is shorter than a 1 mg fish larva — but it does mean the resource and the species are measured differently. Co-Authored-By: Claude Opus 5 --- NEWS.md | 13 +++++ R/ArrayResourceBySize-class.R | 32 ++++++------ R/ArraySpeciesBySize-class.R | 15 +++--- R/animateSpectra.R | 9 ++-- R/plots.R | 50 +++++++++++++++---- R/resource_dynamics.R | 17 +++++++ R/setResource.R | 46 +++++++++++++++++ inst/skills/analyse-and-plot/SKILL.md | 10 ++-- man/addPlot.Rd | 5 +- man/animate.Rd | 9 ++-- man/convert_plot_size_axis.Rd | 3 +- man/plot.ArrayResourceBySize.Rd | 9 ++++ man/plot2.Rd | 5 +- man/plotRelative.Rd | 5 +- man/plot_length_params.Rd | 26 ++++++++++ man/resource_length_defaults.Rd | 34 +++++++++++++ man/resource_length_params.Rd | 20 ++++++++ man/resource_params.Rd | 11 ++++ .../testthat/test-ArrayResourceBySize-class.R | 50 ++++++++++++++++++- tests/testthat/test-plots.R | 12 +++-- tests/testthat/test-setResource.R | 36 +++++++++++++ .../cheatsheet-analysis-and-plotting.Rmd | 10 ++-- 22 files changed, 362 insertions(+), 65 deletions(-) create mode 100644 man/plot_length_params.Rd create mode 100644 man/resource_length_defaults.Rd create mode 100644 man/resource_length_params.Rd diff --git a/NEWS.md b/NEWS.md index 7b2edbbdc..8561e01e2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -146,6 +146,19 @@ stability of steady states. that does not hold a density is now an error; it used to be swallowed silently by `...`. +- The resource can now be shown on length-based plots. `resource_params()` + gains the weight-length parameters `a` and `b`, defaulting to the equivalent + spherical diameter of an organism with the density of water, `a = pi/6` and + `b = 3`, which is the convention plankton ecology uses for a composite of many + taxa. `plotSpectra(params, size_axis = "l")` therefore includes the resource + spectrum, where it used to drop it silently, and the resource array plots and + `animate()` gain `size_axis` and `llim`. The parameters feed none of the + rates. Note that the resource then sits on the length axis at its own + convention: a fish of a given weight is about 3.7 times longer than a sphere + of that weight. That difference is real rather than an artefact — a 1 mg + copepod really is shorter than a 1 mg fish larva — but it does mean the + resource and the species are measured differently. + - `plotYieldObservedVsModel()` gains a `gear` argument that restricts the comparison to the catch of the selected gears. Both the model yield and the observed yield are then taken from those gears only, so in a model where diff --git a/R/ArrayResourceBySize-class.R b/R/ArrayResourceBySize-class.R index df918df59..3b548ab3b 100644 --- a/R/ArrayResourceBySize-class.R +++ b/R/ArrayResourceBySize-class.R @@ -149,9 +149,14 @@ print.summary.ArrayResourceBySize <- function(x, ...) { #' @param wlim A numeric vector of length two providing lower and upper #' limits for the weight (x) axis. Use `NA` to refer to the existing #' minimum or maximum. +#' @param llim A numeric vector of length two providing lower and upper limits +#' for the length (x) axis when `size_axis = "l"`. Use `NA` to refer to the +#' existing minimum or maximum. #' @param ylim A numeric vector of length two providing lower and upper #' limits for the value (y) axis. Use `NA` to refer to the existing #' minimum or maximum. +#' @param size_axis Whether to plot size as weight (`"w"`, default) or length +#' (`"l"`), using the weight-length relationship in [resource_params()]. #' @param per_log_size For an array that holds a density, whether to plot it #' per logarithmic size (`TRUE`) rather than per size (`FALSE`). The default, #' `NULL`, plots the density as it stands. An error for an array that does not @@ -170,9 +175,12 @@ print.summary.ArrayResourceBySize <- function(x, ...) { #' } plot.ArrayResourceBySize <- function(x, return_data = FALSE, log_x = TRUE, log_y = TRUE, log = NULL, - wlim = c(NA, NA), ylim = c(NA, NA), + wlim = c(NA, NA), llim = c(NA, NA), + ylim = c(NA, NA), + size_axis = c("w", "l"), per_log_size = NULL, y_ticks = 6, ...) { + size_axis <- plot_size_axis(size_axis) check_per_log_size(x, per_log_size) log_y <- array_log_y(x, log_y, log, !missing(log_y)) log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y) @@ -180,27 +188,29 @@ plot.ArrayResourceBySize <- function(x, return_data = FALSE, log_y <- log_axes$log_y assert_that(length(wlim) == 2, + length(llim) == 2, length(ylim) == 2) params <- attr(x, "params") plot_dat <- prepare_ArrayResourceBySize_plot_data(x, wlim = wlim) - # The resource has no weight-length relationship, so it has no length axis. - # Expressing a density per logarithmic weight needs no such relationship, - # though, so that much is available here. - plot_dat <- convert_plot_density_axis(plot_dat, params, "w", + plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis, density_wrt = array_density_wrt(x), per_log_size = per_log_size) + if (identical(size_axis, "l")) { + plot_dat <- filter_plot_length_limits(plot_dat, llim) + } if (return_data) return(plot_dat) - y_label <- array_y_label(x, default = "value", per_log_size = per_log_size) + y_label <- array_y_label(x, default = "value", size_axis = size_axis, + per_log_size = per_log_size) ylim <- array_ylim(x, ylim, log_y, plot_dat[[2]]) - plotDataFrame(plot_dat, params, xlab = "Weight (g)", + plotDataFrame(plot_dat, params, xlab = plot_size_xlab(size_axis), ylab = y_label, xtrans = if (log_x) "log10" else "identity", ytrans = if (log_y) "log10" else "identity", - xlim = wlim, ylim = ylim, + xlim = plot_size_xlim(wlim, size_axis, llim), ylim = ylim, y_ticks = y_ticks, legend_var = "Legend") } @@ -721,13 +731,7 @@ animate.ArrayTimeByResourceBySize <- function(x, species = NULL, length(wlim) == 2, length(ylim) == 2, length(tlim) == 2) warn_unused_resource_args(species, total, background) check_per_log_size(x, per_log_size) - # The length axis is derived from each species' weight-length parameters, - # and the resource is not a species. size_axis <- plot_size_axis(size_axis) - if (identical(size_axis, "l")) { - stop("A length axis is not available for resource arrays, because the ", - "weight-length relationship is a species parameter.") - } log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y) log_x <- log_axes$log_x log_y <- log_axes$log_y diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R index c265580f5..65c3f9065 100644 --- a/R/ArraySpeciesBySize-class.R +++ b/R/ArraySpeciesBySize-class.R @@ -473,9 +473,8 @@ parsePlotLog <- function(log, log_x = FALSE, log_y = FALSE) { #' limits for the length (x) axis when `size_axis = "l"`. Use `NA` to #' refer to the existing minimum or maximum.} #' \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or -#' length (`"l"`), using the allometric weight-length relationship. Not -#' available for the resource classes, because the weight-length -#' relationship is a species parameter.} +#' length (`"l"`), using the allometric weight-length relationship of +#' each species, or of the resource, see [resource_params()].} #' \item{`per_log_size`}{For an array that holds a density, whether to plot #' it per logarithmic size (`TRUE`) rather than per size (`FALSE`). The #' default, `NULL`, plots the density as it stands. Unlike `size_axis` @@ -605,9 +604,8 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second", #' limits for the length (x) axis when `size_axis = "l"`. Use `NA` to #' refer to the existing minimum or maximum.} #' \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or -#' length (`"l"`), using the allometric weight-length relationship. Not -#' available for the resource classes, because the weight-length -#' relationship is a species parameter.} +#' length (`"l"`), using the allometric weight-length relationship of +#' each species, or of the resource, see [resource_params()].} #' } #' #' **For `ArrayTimeBySpecies` methods:** @@ -895,9 +893,8 @@ check_per_log_size <- function(x, per_log_size) { #' limits for the length (x) axis when `size_axis = "l"`. Use `NA` to #' refer to the existing minimum or maximum.} #' \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or -#' length (`"l"`), using the allometric weight-length relationship. Not -#' available for the resource classes, because the weight-length -#' relationship is a species parameter.} +#' length (`"l"`), using the allometric weight-length relationship of +#' each species, or of the resource, see [resource_params()].} #' \item{`per_log_size`}{For an array that holds a density, whether to plot #' it per logarithmic size (`TRUE`) rather than per size (`FALSE`). The #' default, `NULL`, plots the density as it stands. Unlike `size_axis` diff --git a/R/animateSpectra.R b/R/animateSpectra.R index c3250fa0b..8b0517f66 100644 --- a/R/animateSpectra.R +++ b/R/animateSpectra.R @@ -24,8 +24,7 @@ #' * **`ArrayTimeByResourceBySize`** — animates the size-resolved resource #' quantity returned by [NResource()] on a `MizerSim` object. There is only a #' single resource spectrum, so `species`, `total` and `background` do nothing -#' and warn if set, and `size_axis = "l"` is an error because the -#' weight-length relationship is a species parameter. +#' and warn if set. #' #' Species linecolours and linetypes follow `params@linecolour` and #' `params@linetype`. @@ -52,9 +51,9 @@ #' it too. An error for an array that does not hold a density. The `MizerSim` #' method has its own `per_log_size`, described below. #' @param size_axis Whether to plot size as weight (`"w"`, default) or length -#' (`"l"`), using the allometric weight-length relationship. Number and -#' biomass densities are transformed to match the chosen axis. The -#' `ArrayTimeByResourceBySize` method supports only `"w"`. +#' (`"l"`), using the allometric weight-length relationship of each species, +#' or of the resource, see [resource_params()]. Number and biomass densities +#' are transformed to match the chosen axis. #' @param total A boolean value that determines whether the total over all #' selected species is plotted as an additional trace called `"Total"`. #' Default is `FALSE`. Not used by the `ArrayTimeByResourceBySize` method, diff --git a/R/plots.R b/R/plots.R index d2b80013e..43ad512e7 100644 --- a/R/plots.R +++ b/R/plots.R @@ -726,10 +726,43 @@ array_log_y <- function(x, log_y, log, given) { log_y } +#' The weight-length parameters to plot each row of plotting data with +#' +#' A length axis needs an allometric weight-length relationship for every line +#' on the plot. The species take theirs from their species parameters and the +#' resource takes its from [resource_params()], where it defaults to the +#' equivalent spherical diameter (see [resource_length_params()]). Anything else +#' — the "Total" row, for instance — has none, and is reported as `NA` so that +#' the caller can leave it out. +#' +#' @param species A vector of the species names in the plotting data. +#' @param params A MizerParams object providing the weight-length parameters. +#' @return A data frame with columns `a` and `b`, one row for each element of +#' `species`, holding `NA` where no relationship is known. +#' @keywords internal +plot_length_params <- function(species, params) { + sp <- params@species_params + if (!all(c("a", "b") %in% names(sp))) { + stop("The species parameter data frame must contain columns ", + "'a' and 'b'.") + } + species <- as.character(species) + idx <- match(species, as.character(sp$species)) + out <- data.frame(a = sp[["a"]][idx], b = sp[["b"]][idx]) + is_resource <- species == "Resource" + if (any(is_resource)) { + rp <- resource_length_params(params) + out$a[is_resource] <- rp$a + out$b[is_resource] <- rp$b + } + out +} + #' Convert plotting data from weight to length #' #' When `size_axis = "l"`, adds a length column `l` computed from the weight -#' column `w` using each species' weight-length relationship. For +#' column `w` using the weight-length relationship of each line, see +#' [plot_length_params()]. Rows with no such relationship are dropped. For #' `size_axis = "w"` the data is returned unchanged. #' #' @param plot_dat A data frame of plotting data with a `w` column and a species @@ -761,10 +794,10 @@ convert_plot_size_axis <- function(plot_dat, params, size_axis, "` column.") } - species <- as.character(plot_dat[[species_col]]) - species_idx <- match(species, as.character(params@species_params$species)) - plot_dat <- plot_dat[!is.na(species_idx), , drop = FALSE] - species_idx <- species_idx[!is.na(species_idx)] + ab <- plot_length_params(plot_dat[[species_col]], params) + known <- !is.na(ab$a) & !is.na(ab$b) + plot_dat <- plot_dat[known, , drop = FALSE] + ab <- ab[known, , drop = FALSE] if (nrow(plot_dat) == 0) { if (!drop_w) { plot_dat$l <- numeric(0) @@ -774,8 +807,7 @@ convert_plot_size_axis <- function(plot_dat, params, size_axis, } return(plot_dat[, setdiff(names(plot_dat), "w"), drop = FALSE]) } - sp <- params@species_params[species_idx, , drop = FALSE] - plot_dat$l <- w2l(plot_dat$w, sp) + plot_dat$l <- (plot_dat$w / ab$a)^(1 / ab$b) if (drop_w) { return(plot_dat[, c("l", setdiff(names(plot_dat), c("l", "w"))), drop = FALSE]) @@ -983,10 +1015,8 @@ convert_plot_density_axis <- function(plot_dat, params, size_axis, } if (needs_jacobian && nrow(plot_dat) > 0) { if (needs_length) { - species_idx <- match(as.character(plot_dat[[species_col]]), - as.character(params@species_params$species)) l <- plot_dat$l - b <- params@species_params$b[species_idx] + b <- plot_length_params(plot_dat[[species_col]], params)$b } else { # Unused by a Jacobian between two weight-based measures. Passing # NA rather than a placeholder keeps a mistake here loud. diff --git a/R/resource_dynamics.R b/R/resource_dynamics.R index 000996839..968ceeeb5 100644 --- a/R/resource_dynamics.R +++ b/R/resource_dynamics.R @@ -59,6 +59,17 @@ resource_constant <- function(params, n_pp, ...) { #' Unlike the carrying capacity, however, the initial resource abundance is #' **not** updated when you subsequently change `kappa` (or call [setResource()]). #' +#' The resource parameters `a` and `b` give the allometric weight-length +#' relationship \eqn{w = a l^b} of the resource, with \eqn{w} in grams and +#' \eqn{l} in centimetres. They feed none of the rates; they exist so that the +#' resource can be shown on the length-based plots (`size_axis = "l"`) alongside +#' the species. They default to the equivalent spherical diameter of an organism +#' with the density of water, \eqn{a = \pi/6} and \eqn{b = 3}, which is the +#' convention plankton ecology uses for a composite of many taxa. This is a +#' different convention from the one the species use, so the resource and the +#' species each sit on the length axis at their own; see +#' [resource_length_defaults]. +#' #' Assigning to `resource_params` only rebuilds the size-dependent resource rate #' and capacity arrays from these scalars (leaving any arrays you have set #' manually untouched). Changing `lambda` also recalculates any `q` and `gamma` @@ -97,6 +108,12 @@ resource_params <- function(params) { if (!is.null(value$r_pp)) { assert_that(is.number(value$r_pp), value$r_pp >= 0) } + if (!is.null(value[["a"]])) { + assert_that(is.number(value[["a"]]), value[["a"]] > 0) + } + if (!is.null(value[["b"]])) { + assert_that(is.number(value[["b"]]), value[["b"]] > 0) + } scalars <- c("kappa", "lambda", "n", "w_pp_cutoff", "r_pp") changed <- scalars[vapply(scalars, function(scalar) { diff --git a/R/setResource.R b/R/setResource.R index 01301effd..3d2e02be7 100644 --- a/R/setResource.R +++ b/R/setResource.R @@ -207,6 +207,13 @@ setResource.MizerParams <- function(params, params@resource_params[["lambda"]] <- lambda params@resource_params[["n"]] <- n params@resource_params[["w_pp_cutoff"]] <- w_pp_cutoff + # The weight-length parameters feed no rate, so they are not arguments + # here; they are filled in so that `resource_params()` shows them and the + # user can change them there. + params@resource_params[["a"]] <- + params@resource_params[["a"]] %||% resource_length_defaults$a + params@resource_params[["b"]] <- + params@resource_params[["b"]] %||% resource_length_defaults$b if (!is.null(resource_capacity) && !is.null(resource_level)) { stop("You should specify only either 'resource_level' or 'resource_capacity'.") @@ -557,6 +564,45 @@ resource_dynamics <- function(params) { setResource(params, resource_dynamics = value, balance = balance) } +#' Default weight-length parameters for the resource +#' +#' The resource is a composite of everything from bacteria to +#' macrozooplankton, so it has no taxonomic length-weight relationship. The +#' default is the geometric one that plankton ecology uses instead: the +#' **equivalent spherical diameter** of an organism with the density of water, +#' \deqn{w = \frac{\pi}{6} l^3,} +#' with \eqn{w} in grams and \eqn{l} in centimetres. On a mizer size grid this +#' puts the smallest resource sizes at a fraction of a micrometre and a +#' milligram organism at about a millimetre, which is the right order for +#' bacteria and copepods respectively. +#' +#' Note that this is a different convention from the one the species use: a fish +#' of a given weight is longer than a sphere of the same weight, by a factor +#' \eqn{(a_{fish}/a_{resource})^{-1/3}}, about 3.7 for the mizer default +#' `a = 0.01`. That difference is real rather than an artefact — a 1 mg copepod +#' really is shorter than a 1 mg fish larva — but it does mean the resource and +#' the species sit on the plot at their own conventions. +#' +#' @format A list with entries `a` and `b`. +#' @seealso [resource_params()] +#' @keywords internal +resource_length_defaults <- list(a = pi / 6, b = 3) + +#' The weight-length parameters of the resource +#' +#' Reads `a` and `b` from [resource_params()], falling back to +#' [resource_length_defaults] for a model that does not set them — which is +#' every model built before these parameters existed. +#' +#' @param params A MizerParams object. +#' @return A list with entries `a` and `b`. +#' @keywords internal +resource_length_params <- function(params) { + rp <- params@resource_params + list(a = rp[["a"]] %||% resource_length_defaults$a, + b = rp[["b"]] %||% resource_length_defaults$b) +} + # The resource capacity that the scalar resource parameters imply, or NULL if # they do not determine it. Used both to recompute a capacity that mizer # controls and to see whether a frozen one has fallen out of step with the diff --git a/inst/skills/analyse-and-plot/SKILL.md b/inst/skills/analyse-and-plot/SKILL.md index 359ddfefb..360273e4f 100644 --- a/inst/skills/analyse-and-plot/SKILL.md +++ b/inst/skills/analyse-and-plot/SKILL.md @@ -354,9 +354,13 @@ plotSpectra(sim, species = c("Cod", "Herring"), resource = FALSE) plotSpectra(sim, biomass = TRUE, size_axis = "l") # biomass density against length ``` -**The resource is dropped from a length axis.** The resource spectrum has no - length-weight relationship, so `plotSpectra(params, size_axis = "l")` shows the - species only. Use `size_axis = "w"` to see the resource. +**The resource has its own length convention.** It is a composite of many taxa, +so instead of a taxonomic weight-length relationship it uses the equivalent +spherical diameter of an organism with the density of water (`a = pi/6`, +`b = 3`, in `resource_params()`). It therefore appears on a length axis, but +measured differently from the fish: a fish of a given weight is about 3.7 times +longer than a sphere of that weight. That gap at the resource-consumer boundary +is real biology, not an artefact. ### Which density a spectrum plot shows diff --git a/man/addPlot.Rd b/man/addPlot.Rd index deb3aeec0..73eb9359a 100644 --- a/man/addPlot.Rd +++ b/man/addPlot.Rd @@ -75,9 +75,8 @@ range (\code{w_min} to \code{w_max}) are removed.} limits for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the existing minimum or maximum.} \item{\code{size_axis}}{Whether to plot size as weight (\code{"w"}, default) or -length (\code{"l"}), using the allometric weight-length relationship. Not -available for the resource classes, because the weight-length -relationship is a species parameter.} +length (\code{"l"}), using the allometric weight-length relationship of +each species, or of the resource, see \code{\link[=resource_params]{resource_params()}}.} \item{\code{per_log_size}}{For an array that holds a density, whether to plot it per logarithmic size (\code{TRUE}) rather than per size (\code{FALSE}). The default, \code{NULL}, plots the density as it stands. Unlike \code{size_axis} diff --git a/man/animate.Rd b/man/animate.Rd index 243211437..e67bd209d 100644 --- a/man/animate.Rd +++ b/man/animate.Rd @@ -65,9 +65,9 @@ for the animated time window, e.g. \code{c(1997, 2007)}. Use \code{NA} to apply limit at that end. Default is \code{c(NA, NA)}.} \item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length -(\code{"l"}), using the allometric weight-length relationship. Number and -biomass densities are transformed to match the chosen axis. The -\code{ArrayTimeByResourceBySize} method supports only \code{"w"}.} +(\code{"l"}), using the allometric weight-length relationship of each species, +or of the resource, see \code{\link[=resource_params]{resource_params()}}. Number and biomass densities +are transformed to match the chosen axis.} \item{per_log_size}{For an array that holds a density, whether to animate it per logarithmic size (\code{TRUE}) rather than per size (\code{FALSE}). The default, @@ -147,8 +147,7 @@ and \code{total} arguments. \item \strong{\code{ArrayTimeByResourceBySize}} — animates the size-resolved resource quantity returned by \code{\link[=NResource]{NResource()}} on a \code{MizerSim} object. There is only a single resource spectrum, so \code{species}, \code{total} and \code{background} do nothing -and warn if set, and \code{size_axis = "l"} is an error because the -weight-length relationship is a species parameter. +and warn if set. } Species linecolours and linetypes follow \code{params@linecolour} and diff --git a/man/convert_plot_size_axis.Rd b/man/convert_plot_size_axis.Rd index 79376f4ad..5adf0a444 100644 --- a/man/convert_plot_size_axis.Rd +++ b/man/convert_plot_size_axis.Rd @@ -32,7 +32,8 @@ dropped) when \code{size_axis = "l"}. } \description{ When \code{size_axis = "l"}, adds a length column \code{l} computed from the weight -column \code{w} using each species' weight-length relationship. For +column \code{w} using the weight-length relationship of each line, see +\code{\link[=plot_length_params]{plot_length_params()}}. Rows with no such relationship are dropped. For \code{size_axis = "w"} the data is returned unchanged. } \keyword{internal} diff --git a/man/plot.ArrayResourceBySize.Rd b/man/plot.ArrayResourceBySize.Rd index 84b7f6bfc..c36a33f65 100644 --- a/man/plot.ArrayResourceBySize.Rd +++ b/man/plot.ArrayResourceBySize.Rd @@ -11,7 +11,9 @@ log_y = TRUE, log = NULL, wlim = c(NA, NA), + llim = c(NA, NA), ylim = c(NA, NA), + size_axis = c("w", "l"), per_log_size = NULL, y_ticks = 6, ... @@ -36,10 +38,17 @@ scales, in the same form as the base \code{\link[=plot]{plot()}} argument. For e limits for the weight (x) axis. Use \code{NA} to refer to the existing minimum or maximum.} +\item{llim}{A numeric vector of length two providing lower and upper limits +for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the +existing minimum or maximum.} + \item{ylim}{A numeric vector of length two providing lower and upper limits for the value (y) axis. Use \code{NA} to refer to the existing minimum or maximum.} +\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length +(\code{"l"}), using the weight-length relationship in \code{\link[=resource_params]{resource_params()}}.} + \item{per_log_size}{For an array that holds a density, whether to plot it per logarithmic size (\code{TRUE}) rather than per size (\code{FALSE}). The default, \code{NULL}, plots the density as it stands. An error for an array that does not diff --git a/man/plot2.Rd b/man/plot2.Rd index 932ef589f..f9294f0a9 100644 --- a/man/plot2.Rd +++ b/man/plot2.Rd @@ -85,9 +85,8 @@ range (\code{w_min} to \code{w_max}) are removed.} limits for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the existing minimum or maximum.} \item{\code{size_axis}}{Whether to plot size as weight (\code{"w"}, default) or -length (\code{"l"}), using the allometric weight-length relationship. Not -available for the resource classes, because the weight-length -relationship is a species parameter.} +length (\code{"l"}), using the allometric weight-length relationship of +each species, or of the resource, see \code{\link[=resource_params]{resource_params()}}.} \item{\code{per_log_size}}{For an array that holds a density, whether to plot it per logarithmic size (\code{TRUE}) rather than per size (\code{FALSE}). The default, \code{NULL}, plots the density as it stands. Unlike \code{size_axis} diff --git a/man/plotRelative.Rd b/man/plotRelative.Rd index 578570185..02d87112c 100644 --- a/man/plotRelative.Rd +++ b/man/plotRelative.Rd @@ -67,9 +67,8 @@ range (\code{w_min} to \code{w_max}) are removed.} limits for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the existing minimum or maximum.} \item{\code{size_axis}}{Whether to plot size as weight (\code{"w"}, default) or -length (\code{"l"}), using the allometric weight-length relationship. Not -available for the resource classes, because the weight-length -relationship is a species parameter.} +length (\code{"l"}), using the allometric weight-length relationship of +each species, or of the resource, see \code{\link[=resource_params]{resource_params()}}.} } \strong{For \code{ArrayTimeBySpecies} methods:} diff --git a/man/plot_length_params.Rd b/man/plot_length_params.Rd new file mode 100644 index 000000000..861f537db --- /dev/null +++ b/man/plot_length_params.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plots.R +\name{plot_length_params} +\alias{plot_length_params} +\title{The weight-length parameters to plot each row of plotting data with} +\usage{ +plot_length_params(species, params) +} +\arguments{ +\item{species}{A vector of the species names in the plotting data.} + +\item{params}{A MizerParams object providing the weight-length parameters.} +} +\value{ +A data frame with columns \code{a} and \code{b}, one row for each element of +\code{species}, holding \code{NA} where no relationship is known. +} +\description{ +A length axis needs an allometric weight-length relationship for every line +on the plot. The species take theirs from their species parameters and the +resource takes its from \code{\link[=resource_params]{resource_params()}}, where it defaults to the +equivalent spherical diameter (see \code{\link[=resource_length_params]{resource_length_params()}}). Anything else +— the "Total" row, for instance — has none, and is reported as \code{NA} so that +the caller can leave it out. +} +\keyword{internal} diff --git a/man/resource_length_defaults.Rd b/man/resource_length_defaults.Rd new file mode 100644 index 000000000..c87561b82 --- /dev/null +++ b/man/resource_length_defaults.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/setResource.R +\name{resource_length_defaults} +\alias{resource_length_defaults} +\title{Default weight-length parameters for the resource} +\format{ +A list with entries \code{a} and \code{b}. +} +\usage{ +resource_length_defaults +} +\description{ +The resource is a composite of everything from bacteria to +macrozooplankton, so it has no taxonomic length-weight relationship. The +default is the geometric one that plankton ecology uses instead: the +\strong{equivalent spherical diameter} of an organism with the density of water, +\deqn{w = \frac{\pi}{6} l^3,} +with \eqn{w} in grams and \eqn{l} in centimetres. On a mizer size grid this +puts the smallest resource sizes at a fraction of a micrometre and a +milligram organism at about a millimetre, which is the right order for +bacteria and copepods respectively. +} +\details{ +Note that this is a different convention from the one the species use: a fish +of a given weight is longer than a sphere of the same weight, by a factor +\eqn{(a_{fish}/a_{resource})^{-1/3}}, about 3.7 for the mizer default +\code{a = 0.01}. That difference is real rather than an artefact — a 1 mg copepod +really is shorter than a 1 mg fish larva — but it does mean the resource and +the species sit on the plot at their own conventions. +} +\seealso{ +\code{\link[=resource_params]{resource_params()}} +} +\keyword{internal} diff --git a/man/resource_length_params.Rd b/man/resource_length_params.Rd new file mode 100644 index 000000000..d235f44a7 --- /dev/null +++ b/man/resource_length_params.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/setResource.R +\name{resource_length_params} +\alias{resource_length_params} +\title{The weight-length parameters of the resource} +\usage{ +resource_length_params(params) +} +\arguments{ +\item{params}{A MizerParams object.} +} +\value{ +A list with entries \code{a} and \code{b}. +} +\description{ +Reads \code{a} and \code{b} from \code{\link[=resource_params]{resource_params()}}, falling back to +\link{resource_length_defaults} for a model that does not set them — which is +every model built before these parameters existed. +} +\keyword{internal} diff --git a/man/resource_params.Rd b/man/resource_params.Rd index 3fde8a5ad..77d0fb125 100644 --- a/man/resource_params.Rd +++ b/man/resource_params.Rd @@ -48,6 +48,17 @@ created: \deqn{N_R(w) = \kappa\, w^{-\lambda}.}{N_R(w) = \kappa w^{-\lambda}.} Unlike the carrying capacity, however, the initial resource abundance is \strong{not} updated when you subsequently change \code{kappa} (or call \code{\link[=setResource]{setResource()}}). +The resource parameters \code{a} and \code{b} give the allometric weight-length +relationship \eqn{w = a l^b} of the resource, with \eqn{w} in grams and +\eqn{l} in centimetres. They feed none of the rates; they exist so that the +resource can be shown on the length-based plots (\code{size_axis = "l"}) alongside +the species. They default to the equivalent spherical diameter of an organism +with the density of water, \eqn{a = \pi/6} and \eqn{b = 3}, which is the +convention plankton ecology uses for a composite of many taxa. This is a +different convention from the one the species use, so the resource and the +species each sit on the length axis at their own; see +\link{resource_length_defaults}. + Assigning to \code{resource_params} only rebuilds the size-dependent resource rate and capacity arrays from these scalars (leaving any arrays you have set manually untouched). Changing \code{lambda} also recalculates any \code{q} and \code{gamma} diff --git a/tests/testthat/test-ArrayResourceBySize-class.R b/tests/testthat/test-ArrayResourceBySize-class.R index 033b52125..2dd215e3f 100644 --- a/tests/testthat/test-ArrayResourceBySize-class.R +++ b/tests/testthat/test-ArrayResourceBySize-class.R @@ -223,8 +223,9 @@ test_that("animate dispatches on ArrayTimeByResourceBySize", { expect_length(frames(animate(n_resource_small, tlim = c(times[2], times[3]))), 2) - expect_error(animate(n_resource_small, size_axis = "l"), - "length axis is not available") + # A length axis is now available, using the resource weight-length + # parameters rather than a species' + expect_s3_class(animate(n_resource_small, size_axis = "l"), "plotly") }) test_that("plot() can express a resource density per logarithmic size", { @@ -246,3 +247,48 @@ test_that("plot() can express a resource density per logarithmic size", { expect_error(plot(resource_level(NS_params_small), per_log_size = TRUE), "only applies to an array that holds a density") }) + +test_that("the resource can be plotted against length", { + params <- setResource(NS_params_small) + resource <- initialNResource(params) + rp <- resource_params(params) + + by_weight <- plot(resource, return_data = TRUE) + by_length <- plot(resource, size_axis = "l", return_data = TRUE) + expect_identical(names(by_length)[[1]], "l") + expect_equal(by_length$l, (by_weight$w / rp$a)^(1 / rp$b)) + + # It is a density, so the values are converted as well as the axis + expect_equal(by_length[[2]], + by_weight[[2]] * rp$b * by_weight$w / by_length$l) + p <- plot(resource, size_axis = "l") + expect_identical(p$scales$get_scales("x")$name, "Length [cm]") + expect_identical(p$scales$get_scales("y")$name, "Number density [1/cm]") + + # llim trims the length axis + expect_lt(nrow(plot(resource, size_axis = "l", llim = c(0.1, 1), + return_data = TRUE)), + nrow(by_length)) + + # A quantity that is not a density keeps its values + mort_w <- plot(getResourceMort(params), return_data = TRUE) + mort_l <- plot(getResourceMort(params), size_axis = "l", + return_data = TRUE) + expect_equal(mort_l[[2]], mort_w[[2]]) + + # and animate() no longer refuses a length axis + expect_s3_class(animate(NResource(NS_sim_small), size_axis = "l"), "plotly") +}) + +test_that("the resource joins the species on a length-axis spectrum", { + params <- setResource(NS_params_small) + on_l <- plotSpectra(params, size_axis = "l", return_data = TRUE) + expect_true("Resource" %in% on_l$Legend) + # The resource sits at its own convention, shorter than a fish of the same + # weight by the ratio of the two `a` values + on_w <- plotSpectra(params, return_data = TRUE) + res_w <- on_w[on_w$Legend == "Resource", ] + res_l <- on_l[on_l$Legend == "Resource", ] + rp <- resource_params(params) + expect_equal(res_l$l, (res_w$w / rp$a)^(1 / rp$b)) +}) diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index b09e521c5..00d8e2c82 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -508,10 +508,14 @@ test_that("size-based plots support length axes", { expect_true(all(spectra_l_limited$l >= llim[1])) expect_true(all(spectra_l_limited$l <= llim[2])) - spectra_hidden <- plotSpectra(params_len, species = species, - resource = TRUE, total = TRUE, - size_axis = "l", return_data = TRUE) - expect_false(any(spectra_hidden$Legend %in% c("Resource", "Total"))) + # The resource has its own weight-length relationship and so appears on a + # length axis; the total still does not, because it is a sum over lines + # that each have their own. + spectra_all <- plotSpectra(params_len, species = species, + resource = TRUE, total = TRUE, + size_axis = "l", return_data = TRUE) + expect_true("Resource" %in% spectra_all$Legend) + expect_false("Total" %in% spectra_all$Legend) p <- plotSpectra(params_len, species = species, resource = FALSE, size_axis = "l") diff --git a/tests/testthat/test-setResource.R b/tests/testthat/test-setResource.R index 913cae820..bddaf4dfe 100644 --- a/tests/testthat/test-setResource.R +++ b/tests/testthat/test-setResource.R @@ -418,3 +418,39 @@ test_that("explicit balancing overrides a frozen complementary array", { balance = FALSE) expect_equal(resource_rate(p3), resource_rate(p), ignore_attr = TRUE) }) + +test_that("the resource carries weight-length parameters with an ESD default", { + params <- setResource(NS_params_small) + rp <- resource_params(params) + # The default is the equivalent spherical diameter at the density of water + expect_equal(rp$a, pi / 6) + expect_equal(rp$b, 3) + expect_equal(mizer:::resource_length_params(params), list(a = pi / 6, b = 3)) + + # A model that predates the parameters falls back to the same default + old <- NS_params_small + old@resource_params$a <- NULL + old@resource_params$b <- NULL + expect_equal(mizer:::resource_length_params(old), list(a = pi / 6, b = 3)) + + # A value the user set is kept, not overwritten by a later setResource() + resource_params(params)$a <- 0.02 + resource_params(params)$b <- 2.9 + expect_equal(resource_params(setResource(params))$a, 0.02) + expect_equal(mizer:::resource_length_params(params), list(a = 0.02, b = 2.9)) + + # and is validated + expect_error({resource_params(params)$a <- -1; params}, "not greater than 0") + expect_error({resource_params(params)$b <- "3"; params}, "not a number") +}) + +test_that("the ESD default puts the resource at plausible lengths", { + # A sanity check on the convention rather than on the code: the smallest + # resource sizes should be micrometres and a milligram organism about a + # millimetre. + esd <- function(w) (w / (pi / 6))^(1 / 3) + expect_equal(esd(1e-12) * 1e4, 1.241, tolerance = 1e-3) # micrometres + expect_equal(esd(1e-3) * 10, 1.241, tolerance = 1e-3) # millimetres + # A fish of the same weight is longer, by a fixed factor + expect_equal(((1e-3 / 0.01)^(1 / 3)) / esd(1e-3), 3.741, tolerance = 1e-3) +}) diff --git a/vignettes/cheatsheet-analysis-and-plotting.Rmd b/vignettes/cheatsheet-analysis-and-plotting.Rmd index 99f8cce76..f3caf46ba 100644 --- a/vignettes/cheatsheet-analysis-and-plotting.Rmd +++ b/vignettes/cheatsheet-analysis-and-plotting.Rmd @@ -369,9 +369,13 @@ plotSpectra(sim, species = c("Cod", "Herring"), resource = FALSE) plotSpectra(sim, biomass = TRUE, size_axis = "l") # biomass density against length ``` -**The resource is dropped from a length axis.** The resource spectrum has no - length-weight relationship, so `plotSpectra(params, size_axis = "l")` shows the - species only. Use `size_axis = "w"` to see the resource. +**The resource has its own length convention.** It is a composite of many taxa, +so instead of a taxonomic weight-length relationship it uses the equivalent +spherical diameter of an organism with the density of water (`a = pi/6`, +`b = 3`, in [`resource_params()`](../reference/resource_params.html)). It therefore appears on a length axis, but +measured differently from the fish: a fish of a given weight is about 3.7 times +longer than a sphere of that weight. That gap at the resource-consumer boundary +is real biology, not an artefact. ### Which density a spectrum plot shows