Second-order plotting: draw bin averages at geometric bin centres (#382)

NEWS.md

@@ -1,5 +1,15 @@
 # mizer (development version)
 
+- Size-resolved diagnostics that are finite-volume bin averages (the
+  mortalities `getPredMort()`, `getFMort()`, `getMort()`, `getExtMort()`, and
+  the reproductive investment `getERepro()`) are now drawn at the geometric bin
+  centre `sqrt(w_j w_{j+1})` rather than the left bin edge, the location where a
+  bin average actually lives. The `ArraySpeciesBySize`/`ArrayTimeBySpeciesBySize`
+  classes carry a `representation` tag (`"point"`/`"average"`) recording this,
+  and the shift is applied only when the model uses second-order bin-averaging
+  (`second_order_w[["bin_average"]]`), so default plots are unchanged.
+  Point-valued quantities (encounter, growth) stay on the grid nodes. (#382)
+
 - `plotYield()` now uses `sim2 = NULL` instead of `missing(sim2)` to detect
   the optional second simulation argument. This is backward-compatible and
   makes the function work correctly with `do.call()`.

R/ArraySpeciesBySize-class.R

@@ -23,6 +23,11 @@
 #' @param units A string giving the units (e.g. "g/year", "1/year").
 #' @param params A `MizerParams` object. Used for species colours, linetypes,
 #'   and size ranges in the `plot()` method.
+#' @param representation Either `"point"` (the default) for a quantity sampled
+#'   at the grid nodes, or `"average"` for a finite-volume bin average. A
+#'   bin-averaged quantity is drawn at the geometric bin centre rather than the
+#'   left bin edge, but only when the model uses second-order bin-averaging
+#'   (`second_order_w[["bin_average"]]`), so default plots are unchanged.
 #'
 #' @return An `ArraySpeciesBySize` object (inherits from `matrix` and `array`).
 #' @seealso [print()], [summary()], [as.data.frame()], [plot()]
@@ -34,18 +39,21 @@
 #' summary(enc)
 #' }
 ArraySpeciesBySize <- function(x, value_name = NULL, units = NULL,
-                               params = NULL) {
+                               params = NULL,
+                               representation = c("point", "average")) {
     if (!is.matrix(x)) {
         stop("`x` must be a matrix.")
     }
+    representation <- match.arg(representation)
     if (!is.null(params) && identical(dim(x), dim(params@metab))) {
         dimnames(x) <- dimnames(params@metab)
     }
     structure(x,
         class = c("ArraySpeciesBySize", "matrix", "array"),
         value_name = value_name,
         units = units,
-        params = params
+        params = params,
+        representation = representation
     )
 }
 
@@ -860,7 +868,13 @@ as.data.frame.ArraySpeciesBySize <- function(x, row.names = NULL,
 #'
 #' @param x An `ArraySpeciesBySize` object.
 #'
-#' @return A numeric vector giving the size represented by each column.
+#' @return A numeric vector giving the size represented by each column. When the
+#'   array is tagged as a bin average (`representation = "average"`) *and* the
+#'   model uses second-order bin-averaging (`second_order_w[["bin_average"]]`),
+#'   the geometric bin centres are returned instead of the left bin edges, so
+#'   that bin-averaged quantities are drawn at the size where they actually live
+#'   (see [bin_midpoints()]). Point-valued quantities and first-order models are
+#'   unaffected, keeping default plots unchanged.
 #' @keywords internal
 get_ArraySpeciesBySize_w <- function(x) {
     params <- attr(x, "params")
@@ -871,11 +885,13 @@ get_ArraySpeciesBySize_w <- function(x) {
         }
         return(w)
     }
+    average <- identical(attr(x, "representation"), "average") &&
+        isTRUE(params@second_order_w[["bin_average"]])
     if (ncol(x) == length(params@w)) {
-        return(params@w)
+        return(if (average) bin_midpoints(params) else params@w)
     }
     if (ncol(x) == length(params@w_full)) {
-        return(params@w_full)
+        return(if (average) bin_midpoints(params, w_full = TRUE) else params@w_full)
     }
     stop("Can not determine the size grid for this ArraySpeciesBySize object. ",
          "The number of columns is ", ncol(x), ", but the params object has ",
@@ -891,6 +907,7 @@ get_ArraySpeciesBySize_w <- function(x) {
         attr(result, "value_name") <- attr(x, "value_name")
         attr(result, "units") <- attr(x, "units")
         attr(result, "params") <- attr(x, "params")
+        attr(result, "representation") <- attr(x, "representation")
         class(result) <- c("ArraySpeciesBySize", "matrix", "array")
     }
     result

R/ArrayTimeBySpeciesBySize-class.R

@@ -25,6 +25,11 @@
 #' @param units A string giving the units (e.g. "1/year").
 #' @param params A `MizerParams` object. Used for species colours, linetypes,
 #'   and size ranges in the `plot()` and `animateSpectra()` methods.
+#' @param representation Either `"point"` (the default) for a quantity sampled
+#'   at the grid nodes, or `"average"` for a finite-volume bin average, which is
+#'   then drawn at the geometric bin centre when the model uses second-order
+#'   bin-averaging (`second_order_w[["bin_average"]]`). See
+#'   [ArraySpeciesBySize()].
 #'
 #' @return An `ArrayTimeBySpeciesBySize` object (inherits from `array`).
 #' @seealso [print()], [summary()], [as.data.frame()], [plot()],
@@ -38,18 +43,46 @@
 #' plot(fmort, time = 2007)
 #' }
 ArrayTimeBySpeciesBySize <- function(x, value_name = NULL, units = NULL,
-                                     params = NULL) {
+                                     params = NULL,
+                                     representation = c("point", "average")) {
     if (!is.array(x) || length(dim(x)) != 3) {
         stop("`x` must be a 3D array.")
     }
+    representation <- match.arg(representation)
     structure(x,
         class = c("ArrayTimeBySpeciesBySize", "array"),
         value_name = value_name,
         units = units,
-        params = params
+        params = params,
+        representation = representation
     )
 }
 
+#' Get the size grid for an ArrayTimeBySpeciesBySize object
+#'
+#' Internal helper, the three-dimensional analogue of
+#' [get_ArraySpeciesBySize_w()]. Returns the geometric bin centres (see
+#' [bin_midpoints()]) when the array is tagged as a bin average and the model
+#' uses second-order bin-averaging, otherwise the grid nodes read from the size
+#' dimension names. Falls back to the dimension names when no `params` is
+#' attached.
+#'
+#' @param x An `ArrayTimeBySpeciesBySize` object.
+#' @return A numeric vector giving the size represented by each size slice.
+#' @keywords internal
+get_ArrayTimeBySpeciesBySize_w <- function(x) {
+    w <- as.numeric(dimnames(x)[[3]])
+    if (any(is.na(w))) w <- seq_len(dim(x)[3])
+    params <- attr(x, "params")
+    average <- !is.null(params) &&
+        identical(attr(x, "representation"), "average") &&
+        isTRUE(params@second_order_w[["bin_average"]])
+    if (!average) return(w)
+    if (dim(x)[3] == length(params@w)) return(bin_midpoints(params))
+    if (dim(x)[3] == length(params@w_full)) return(bin_midpoints(params, w_full = TRUE))
+    w
+}
+
 #' Test if an object is an ArrayTimeBySpeciesBySize
 #'
 #' @param x An object to test.
@@ -230,6 +263,7 @@ ArrayTimeBySpeciesBySize_slice <- function(x, time = NULL) {
     params <- attr(x, "params")
     value_name <- attr(x, "value_name")
     units <- attr(x, "units")
+    representation <- attr(x, "representation") %||% "point"
 
     times <- as.numeric(dimnames(x)[[1]])
     if (is.null(time)) {
@@ -243,7 +277,8 @@ ArrayTimeBySpeciesBySize_slice <- function(x, time = NULL) {
                     nrow = dim(arr)[2],
                     dimnames = dimnames(arr)[2:3])
     ArraySpeciesBySize(slice, value_name = value_name,
-                       units = units, params = params)
+                       units = units, params = params,
+                       representation = representation)
 }
 
 #' @rdname plotHover
@@ -316,7 +351,7 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
         times <- times[times <= tlim[2]]
     }
 
-    w <- as.numeric(dimnames(arr)[[3]])
+    w <- get_ArrayTimeBySpeciesBySize_w(x)
     sub <- arr[, species, , drop = FALSE]
 
     # Build long-format data frame; time varies fastest to match c(sub)
@@ -366,8 +401,7 @@ as.data.frame.ArrayTimeBySpeciesBySize <- function(x, row.names = NULL,
     times <- as.numeric(dimnames(x)[[1]])
     if (any(is.na(times))) times <- seq_len(dim(x)[1])
     sp_names <- dimnames(x)[[2]]
-    w <- as.numeric(dimnames(x)[[3]])
-    if (any(is.na(w))) w <- seq_len(dim(x)[3])
+    w <- get_ArrayTimeBySpeciesBySize_w(x)
 
     data.frame(
         expand.grid(time = times, Species = sp_names, w = w,
@@ -385,17 +419,20 @@ as.data.frame.ArrayTimeBySpeciesBySize <- function(x, row.names = NULL,
         attr(result, "value_name") <- attr(x, "value_name")
         attr(result, "units") <- attr(x, "units")
         attr(result, "params") <- attr(x, "params")
+        attr(result, "representation") <- attr(x, "representation")
         class(result) <- c("ArrayTimeBySpeciesBySize", "array")
     } else if (is.matrix(result)) {
         dim_names <- names(dimnames(result))
         attrs <- list(value_name = attr(x, "value_name"),
                       units = attr(x, "units"),
-                      params = attr(x, "params"))
+                      params = attr(x, "params"),
+                      representation = attr(x, "representation") %||% "point")
         if (identical(dim_names, c("sp", "w"))) {
             result <- ArraySpeciesBySize(result,
                                          value_name = attrs$value_name,
                                          units = attrs$units,
-                                         params = attrs$params)
+                                         params = attrs$params,
+                                         representation = attrs$representation)
         } else if (identical(dim_names, c("time", "sp"))) {
             result <- ArrayTimeBySpecies(result,
                                          value_name = attrs$value_name,

R/helpers.R

@@ -66,6 +66,33 @@ power_law_bin_average <- function(w, dw, d, w_max = Inf) {
     }
 }
 
+#' Geometric bin centres of the size grid
+#'
+#' Internal helper for the second-order plotting code. A finite-volume cell
+#' average \eqn{N_j = (1/\Delta w_j)\int_{w_j}^{w_{j+1}} N\,dw} does not live at
+#' the left bin edge \eqn{w_j} but at the geometric bin centre
+#' \deqn{w^*_j = \sqrt{w_j\,w_{j+1}} = w_j\sqrt{\beta},}
+#' where \eqn{\beta = w_{j+1}/w_j} is the (constant) bin ratio of the
+#' logarithmic grid. This is the log-symmetric, second-order-correct location at
+#' which to plot a bin-averaged quantity (it is exact for the community spectrum
+#' \eqn{N\propto w^{-2}}). It is a uniform half-bin shift to the right on the log
+#' axis, the same for the consumer grid `w` and the full prey/resource grid
+#' `w_full`.
+#'
+#' @param params A MizerParams object.
+#' @param w_full If `TRUE`, return the centres of the full (resource) grid
+#'   `params@w_full`; otherwise the consumer grid `params@w`.
+#' @return A numeric vector of geometric bin centres, one per grid node.
+#' @concept helper
+#' @keywords internal
+bin_midpoints <- function(params, w_full = FALSE) {
+    # The grid is geometric, so the bin ratio beta is constant across the whole
+    # (consumer and resource) grid.
+    beta <- params@w_full[2] / params@w_full[1]
+    w <- if (w_full) params@w_full else params@w
+    w * sqrt(beta)
+}
+
 #' Length-weight conversion
 #'
 #' For each species, convert between length and weight using the relationship

R/plots.R

@@ -2539,7 +2539,8 @@ plotPredMort.MizerSim <- function(object, species = NULL,
     pred_mort <- ArraySpeciesBySize(pred_mort,
                                     value_name = "Predation mortality",
                                     units = "1/year",
-                                    params = object@params)
+                                    params = object@params,
+                                    representation = "average")
     plot(pred_mort, species = species, all.sizes = all.sizes,
          highlight = highlight, wlim = wlim, llim = llim,
          log_x = log_x, log_y = log_y, log = log,
@@ -2674,7 +2675,8 @@ plotFMort.MizerSim <- function(object, species = NULL,
         f <- apply(f, c(2, 3), mean)
     }
     f <- ArraySpeciesBySize(f, value_name = "Fishing mortality",
-                            units = "1/year", params = object@params)
+                            units = "1/year", params = object@params,
+                            representation = "average")
     plot(f, species = species, all.sizes = all.sizes,
          highlight = highlight, wlim = wlim, llim = llim,
          log_x = log_x, log_y = log_y, log = log,