Give the resource a weight-length relationship

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

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,37 +175,42 @@ 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)
     log_x <- log_axes$log_x
     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

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`

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,

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.

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) {

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

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