From 7181201dffbf55091d3d4baee176745aeb1f565f Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Thu, 11 Jun 2026 18:21:45 +0100 Subject: [PATCH 1/3] Second-order plotting: draw bin averages at geometric bin centres (#382) Once a quantity is a finite-volume cell average N_j it no longer lives at the left bin edge w_j but at the geometric bin centre w*_j = sqrt(w_j w_{j+1}) = w_j sqrt(beta) (the log-midpoint, exact for the community spectrum w^-2). Plotting it at the bin boundary mis-places it by half a bin. This adds a `representation` tag ("point"/"average") to the ArraySpeciesBySize and ArrayTimeBySpeciesBySize classes, threaded through the constructors, subset/slice methods and the MizerSim rate path (sim_size_rate / get_species_size_rate_from_sim). The central size-axis helpers get_ArraySpeciesBySize_w() and the new get_ArrayTimeBySpeciesBySize_w() return the geometric bin centres (new bin_midpoints() helper) for "average" arrays, but only when the model uses second-order bin-averaging (second_order_w[["bin_average"]]) -- so default plots, as.data.frame() output and vdiffr snapshots are unchanged. The bin-averaged sinks getPredMort/getFMort/getMort/getExtMort and the reproductive investment getERepro tag themselves "average"; point-valued quantities (encounter, growth) stay "point" and on the nodes. Everything downstream (prepare_*_plot_data, as.data.frame, plot2) follows the central helpers automatically. plotSpectra/plotCDF build their own data frames and handle the w^power weight; their centre placement is the separate #383. Snapshots gained only the new benign `representation` attribute (accepted). Full suite: 2824 pass. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 10 ++++ R/ArraySpeciesBySize-class.R | 27 +++++++-- R/ArrayTimeBySpeciesBySize-class.R | 53 ++++++++++++++--- R/helpers.R | 27 +++++++++ R/plots.R | 6 +- R/rate_functions.R | 43 +++++++++----- R/setExtMort.R | 6 +- man/ArraySpeciesBySize.Rd | 14 ++++- man/ArrayTimeBySpeciesBySize.Rd | 14 ++++- man/bin_midpoints.Rd | 31 ++++++++++ man/get_ArraySpeciesBySize_w.Rd | 8 ++- man/get_ArrayTimeBySpeciesBySize_w.Rd | 23 ++++++++ man/get_species_size_rate_from_sim.Rd | 3 +- man/sim_size_rate.Rd | 1 + .../_snaps/backwards_compatibility.md | 6 ++ tests/testthat/_snaps/project_methods.md | 45 +++++++++++++++ tests/testthat/test-ArraySpeciesBySize.R | 57 +++++++++++++++++++ .../testthat/test-ArrayTimeBySpeciesBySize.R | 31 ++++++++++ tests/testthat/test-project_methods.R | 3 +- vignettes/numerical_details.qmd | 2 + 20 files changed, 375 insertions(+), 35 deletions(-) create mode 100644 man/bin_midpoints.Rd create mode 100644 man/get_ArrayTimeBySpeciesBySize_w.Rd diff --git a/NEWS.md b/NEWS.md index da2aafb09..8dabcf37f 100644 --- a/NEWS.md +++ b/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()`. diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R index c0c39976b..b28247f6b 100644 --- a/R/ArraySpeciesBySize-class.R +++ b/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,10 +39,12 @@ #' 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) } @@ -45,7 +52,8 @@ ArraySpeciesBySize <- function(x, value_name = NULL, units = NULL, 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 diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R index 8d0a8053a..7112556e9 100644 --- a/R/ArrayTimeBySpeciesBySize-class.R +++ b/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, diff --git a/R/helpers.R b/R/helpers.R index 30b657c87..88461f7b4 100644 --- a/R/helpers.R +++ b/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 diff --git a/R/plots.R b/R/plots.R index b69d81731..ca6ff584c 100644 --- a/R/plots.R +++ b/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, diff --git a/R/rate_functions.R b/R/rate_functions.R index b78098824..f3f5ca9b4 100644 --- a/R/rate_functions.R +++ b/R/rate_functions.R @@ -200,7 +200,8 @@ get_sim_rate_slice <- function(sim, time_idx) { #' @keywords internal get_species_size_rate_from_sim <- function(sim, time_range, drop, rate_fun, value_name, - units = NULL) { + units = NULL, + representation = "point") { time_elements <- get_sim_rate_time_elements(sim, time_range) # Apply the one-time rate calculation to each selected saved time. The @@ -220,11 +221,13 @@ get_species_size_rate_from_sim <- function(sim, time_range, drop, result <- ArrayTimeBySpeciesBySize(result, value_name = value_name, units = units, - params = sim@params) + params = sim@params, + representation = representation) } else if (is.matrix(result) && names(dimnames(result))[[1]] == "sp") { result <- ArraySpeciesBySize(result, value_name = value_name, - units = units, params = sim@params) + units = units, params = sim@params, + representation = representation) } result } @@ -409,7 +412,8 @@ mizer_rates_subset <- function(params, n, n_pp, n_other, t, effort, #' @keywords internal sim_size_rate <- function(sim, time_range, drop, target, slot, value_name, units = NULL, - use_sim_effort = FALSE, ...) { + use_sim_effort = FALSE, + representation = "point", ...) { params <- validParams(sim@params) rates_fns <- projectRateFunctions(params) effort <- params@initial_effort @@ -435,7 +439,8 @@ sim_size_rate <- function(sim, time_range, drop, target, slot, } m }, - value_name = value_name, units = units) + value_name = value_name, units = units, + representation = representation) } #' Build a `MizerSim` rate getter that resolves the rate functions once @@ -758,8 +763,11 @@ getPredMort.MizerParams <- function(object, n, n_pp, n_other, pred_mort <- f(params, n = n, n_pp = n_pp, n_other = n_other, t = t, pred_rate = pred_rate) } + # Predation mortality is a sink integrated against the prey density over the + # prey bin, so under second-order bin-averaging it is a bin average and is + # plotted at the geometric bin centre (see get_ArraySpeciesBySize_w()). ArraySpeciesBySize(pred_mort, value_name = "Predation mortality", - units = "1/year", params = params) + units = "1/year", params = params, representation = "average") } #' @export @@ -768,7 +776,7 @@ getPredMort.MizerSim <- function(object, n, n_pp, n_other, sim <- object sim_size_rate(sim, time_range, drop, target = "PredMort", slot = "pred_mort", value_name = "Predation mortality", - units = "1/year", ...) + units = "1/year", representation = "average", ...) } #' Alias for `getPredMort()` @@ -1119,7 +1127,8 @@ getFMort.MizerParams <- function(object, effort, time_range, drop = TRUE, n_other = n_other, time_range = t)) fmort <- do.call(f, c(args, list(...))) fmort <- ArraySpeciesBySize(fmort, value_name = "Fishing mortality", - units = "1/year", params = params) + units = "1/year", params = params, + representation = "average") return(fmort) } else if (length(effort) == no_gears) { args <- list( @@ -1131,7 +1140,8 @@ getFMort.MizerParams <- function(object, effort, time_range, drop = TRUE, n_other = n_other, time_range = t)) fmort <- do.call(f, c(args, list(...))) fmort <- ArraySpeciesBySize(fmort, value_name = "Fishing mortality", - units = "1/year", params = params) + units = "1/year", params = params, + representation = "average") return(fmort) } else { stop("Invalid effort argument") @@ -1144,7 +1154,7 @@ getFMort.MizerSim <- function(object, effort, time_range, drop = TRUE, sim <- object sim_size_rate(sim, time_range, drop, target = "FMort", slot = "f_mort", value_name = "Fishing mortality", units = "1/year", - use_sim_effort = TRUE, ...) + use_sim_effort = TRUE, representation = "average", ...) } @@ -1218,7 +1228,8 @@ getMort.MizerParams <- function(object, t = t, effort = effort, rates_fns = rates_fns, targets = "Mort", ...) return(ArraySpeciesBySize(r$mort, value_name = "Total mortality", - units = "1/year", params = params)) + units = "1/year", params = params, + representation = "average")) } #' @export @@ -1228,7 +1239,7 @@ getMort.MizerSim <- function(object, n, n_pp, n_other, effort, t, ..., if (missing(time_range) && !missing(t)) time_range <- t sim_size_rate(sim, time_range, drop, target = "Mort", slot = "mort", value_name = "Total mortality", units = "1/year", - use_sim_effort = TRUE, ...) + use_sim_effort = TRUE, representation = "average", ...) } #' Alias for `getMort()` @@ -1292,8 +1303,11 @@ getERepro.MizerParams <- function(object, n = initialN(object), erepro <- f(params, n = n, n_pp = n_pp, n_other = n_other, t = t, e = e) } + # The per-capita reproductive investment is the integrand of the + # reproduction integral, where (under second-order bin-averaging) it enters + # as a bin average, so it is plotted at the geometric bin centre. ArraySpeciesBySize(erepro, value_name = "Energy for reproduction", - units = "g/year", params = params) + units = "g/year", params = params, representation = "average") } #' @export @@ -1302,7 +1316,8 @@ getERepro.MizerSim <- function(object, n, n_pp, n_other, t, ..., sim <- object if (missing(time_range) && !missing(t)) time_range <- t sim_size_rate(sim, time_range, drop, target = "ERepro", slot = "e_repro", - value_name = "Energy for reproduction", units = "g/year", ...) + value_name = "Energy for reproduction", units = "g/year", + representation = "average", ...) } #' Alias for `getERepro()` diff --git a/R/setExtMort.R b/R/setExtMort.R index 8a4acddfe..57e6c5afe 100644 --- a/R/setExtMort.R +++ b/R/setExtMort.R @@ -178,10 +178,14 @@ getExtMort <- function(params) { } #' @export getExtMort.MizerParams <- function(params) { + # External mortality is a sink integrated against the abundance over the + # bin; under second-order bin-averaging mu_b is the exact bin average + # (see setExtMort()), so it is plotted at the geometric bin centre. ArraySpeciesBySize(params@mu_b, value_name = "External mortality", units = "1/year", - params = params) + params = params, + representation = "average") } #' @rdname setExtMort diff --git a/man/ArraySpeciesBySize.Rd b/man/ArraySpeciesBySize.Rd index 8bb7d6369..f30bdd2ac 100644 --- a/man/ArraySpeciesBySize.Rd +++ b/man/ArraySpeciesBySize.Rd @@ -4,7 +4,13 @@ \alias{ArraySpeciesBySize} \title{S3 class for species x size rate arrays} \usage{ -ArraySpeciesBySize(x, value_name = NULL, units = NULL, params = NULL) +ArraySpeciesBySize( + x, + value_name = NULL, + units = NULL, + params = NULL, + representation = c("point", "average") +) } \arguments{ \item{x}{A matrix (species x size).} @@ -15,6 +21,12 @@ ArraySpeciesBySize(x, value_name = NULL, units = NULL, params = NULL) \item{params}{A \code{MizerParams} object. Used for species colours, linetypes, and size ranges in the \code{plot()} method.} + +\item{representation}{Either \code{"point"} (the default) for a quantity sampled +at the grid nodes, or \code{"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 +(\code{second_order_w[["bin_average"]]}), so default plots are unchanged.} } \value{ An \code{ArraySpeciesBySize} object (inherits from \code{matrix} and \code{array}). diff --git a/man/ArrayTimeBySpeciesBySize.Rd b/man/ArrayTimeBySpeciesBySize.Rd index b06c11c40..ff4205f29 100644 --- a/man/ArrayTimeBySpeciesBySize.Rd +++ b/man/ArrayTimeBySpeciesBySize.Rd @@ -4,7 +4,13 @@ \alias{ArrayTimeBySpeciesBySize} \title{S3 class for time x species x size arrays} \usage{ -ArrayTimeBySpeciesBySize(x, value_name = NULL, units = NULL, params = NULL) +ArrayTimeBySpeciesBySize( + x, + value_name = NULL, + units = NULL, + params = NULL, + representation = c("point", "average") +) } \arguments{ \item{x}{A 3D array (time x species x size).} @@ -15,6 +21,12 @@ ArrayTimeBySpeciesBySize(x, value_name = NULL, units = NULL, params = NULL) \item{params}{A \code{MizerParams} object. Used for species colours, linetypes, and size ranges in the \code{plot()} and \code{animateSpectra()} methods.} + +\item{representation}{Either \code{"point"} (the default) for a quantity sampled +at the grid nodes, or \code{"average"} for a finite-volume bin average, which is +then drawn at the geometric bin centre when the model uses second-order +bin-averaging (\code{second_order_w[["bin_average"]]}). See +\code{\link[=ArraySpeciesBySize]{ArraySpeciesBySize()}}.} } \value{ An \code{ArrayTimeBySpeciesBySize} object (inherits from \code{array}). diff --git a/man/bin_midpoints.Rd b/man/bin_midpoints.Rd new file mode 100644 index 000000000..1bd5038d8 --- /dev/null +++ b/man/bin_midpoints.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{bin_midpoints} +\alias{bin_midpoints} +\title{Geometric bin centres of the size grid} +\usage{ +bin_midpoints(params, w_full = FALSE) +} +\arguments{ +\item{params}{A MizerParams object.} + +\item{w_full}{If \code{TRUE}, return the centres of the full (resource) grid +\code{params@w_full}; otherwise the consumer grid \code{params@w}.} +} +\value{ +A numeric vector of geometric bin centres, one per grid node. +} +\description{ +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 \code{w} and the full prey/resource grid +\code{w_full}. +} +\concept{helper} +\keyword{internal} diff --git a/man/get_ArraySpeciesBySize_w.Rd b/man/get_ArraySpeciesBySize_w.Rd index 8a072f195..01f8e5aa4 100644 --- a/man/get_ArraySpeciesBySize_w.Rd +++ b/man/get_ArraySpeciesBySize_w.Rd @@ -10,7 +10,13 @@ get_ArraySpeciesBySize_w(x) \item{x}{An \code{ArraySpeciesBySize} object.} } \value{ -A numeric vector giving the size represented by each column. +A numeric vector giving the size represented by each column. When the +array is tagged as a bin average (\code{representation = "average"}) \emph{and} the +model uses second-order bin-averaging (\code{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 \code{\link[=bin_midpoints]{bin_midpoints()}}). Point-valued quantities and first-order models are +unaffected, keeping default plots unchanged. } \description{ Internal helper that returns the consumer size grid \code{params@w} or the full diff --git a/man/get_ArrayTimeBySpeciesBySize_w.Rd b/man/get_ArrayTimeBySpeciesBySize_w.Rd new file mode 100644 index 000000000..185ab6d94 --- /dev/null +++ b/man/get_ArrayTimeBySpeciesBySize_w.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ArrayTimeBySpeciesBySize-class.R +\name{get_ArrayTimeBySpeciesBySize_w} +\alias{get_ArrayTimeBySpeciesBySize_w} +\title{Get the size grid for an ArrayTimeBySpeciesBySize object} +\usage{ +get_ArrayTimeBySpeciesBySize_w(x) +} +\arguments{ +\item{x}{An \code{ArrayTimeBySpeciesBySize} object.} +} +\value{ +A numeric vector giving the size represented by each size slice. +} +\description{ +Internal helper, the three-dimensional analogue of +\code{\link[=get_ArraySpeciesBySize_w]{get_ArraySpeciesBySize_w()}}. Returns the geometric bin centres (see +\code{\link[=bin_midpoints]{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 \code{params} is +attached. +} +\keyword{internal} diff --git a/man/get_species_size_rate_from_sim.Rd b/man/get_species_size_rate_from_sim.Rd index f1a492db2..78d706ad7 100644 --- a/man/get_species_size_rate_from_sim.Rd +++ b/man/get_species_size_rate_from_sim.Rd @@ -10,7 +10,8 @@ get_species_size_rate_from_sim( drop, rate_fun, value_name, - units = NULL + units = NULL, + representation = "point" ) } \arguments{ diff --git a/man/sim_size_rate.Rd b/man/sim_size_rate.Rd index 3a06319f4..1e7ee2dcc 100644 --- a/man/sim_size_rate.Rd +++ b/man/sim_size_rate.Rd @@ -13,6 +13,7 @@ sim_size_rate( value_name, units = NULL, use_sim_effort = FALSE, + representation = "point", ... ) } diff --git a/tests/testthat/_snaps/backwards_compatibility.md b/tests/testthat/_snaps/backwards_compatibility.md index c935db0a2..11aa3bc45 100644 --- a/tests/testthat/_snaps/backwards_compatibility.md +++ b/tests/testthat/_snaps/backwards_compatibility.md @@ -1991,6 +1991,8 @@ Haddock 156333250686 140753274052 126067041926 112316282311 99528337355 Cod 1859436226 1297418578 906990223 640148493 460469325 Saithe 106982064 75453886 57866920 48355351 43425648 + attr(,"representation") + [1] "point" # set_trait_model() works as in version 1 @@ -4235,6 +4237,8 @@ 9 0.01326478 0.01316974 0.01307207 0.01297159 0.01286811 0.01276144 10 0.01326478 0.01316974 0.01307207 0.01297159 0.01286811 0.01276144 11 0.01326478 0.01316974 0.01307207 0.01297159 0.01286811 0.01276144 + attr(,"representation") + [1] "point" # set_community_model() works as in version 1 @@ -4752,4 +4756,6 @@ w sp 811000 1e+06 Community 3539.996 3503.139 + attr(,"representation") + [1] "point" diff --git a/tests/testthat/_snaps/project_methods.md b/tests/testthat/_snaps/project_methods.md index b4f2d8bd8..b53f386b4 100644 --- a/tests/testthat/_snaps/project_methods.md +++ b/tests/testthat/_snaps/project_methods.md @@ -8,6 +8,11 @@ "attributes": {}, "value": [3, 100] }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["point"] + }, "dimnames": { "type": "list", "attributes": { @@ -90,6 +95,11 @@ "type": "character", "attributes": {}, "value": ["1/year"] + }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["point"] } }, "value": [1.14752828e-10, 7.69233583e-07, 5.16467603e-14, 7.46511953e-11, 1.06515634e-06, 5.53705909e-14, 1.46390086e-10, 1.47213888e-06, 2.8970455e-14, 1.37707922e-10, 2.03078552e-06, 5.39136062e-14, 1.60514742e-10, 2.79614241e-06, 2.93587784e-14, 1.53408612e-10, 3.8426749e-06, 5.70878716e-14, 1.95406033e-10, 5.27091744e-06, 4.95424293e-14, 1.56063416e-10, 7.21631356e-06, 5.76256677e-14, 1.3895232e-10, 9.86096475e-06, 2.75090857e-14, 1.31602487e-10, 0.00001345, 3.72871968e-14, 1.65180536e-10, 0.00001831, 2.55534635e-14, 1.69778775e-10, 0.00002487, 3.3636702e-14, 6.40351108e-11, 0.00003373, 1.04299851e-15, 8.17832111e-11, 0.00004565, 4.17199404e-14, 1.01841897e-10, 0.00006167, 0, 9.53648511e-11, 0.00008314, 2.39889658e-14, 6.53006502e-11, 0.00011186, 6.25799107e-15, 1.29129904e-10, 0.00015022, 1.87739732e-14, 1.04566535e-10, 0.00020132, 6.25799107e-15, 1.00751742e-10, 0.00026926, 0, 8.84694394e-11, 0.00035942, 0, 1.05453207e-10, 0.00047878, 0, 4.73014329e-11, 0.00063648, 0, 9.25684775e-11, 0.00084439, 0, 5.49533183e-11, 0.00111789, 0, 8.94385835e-11, 0.00147688, 0, 1.66739059e-11, 0.00194706, 2.50319643e-14, 3.99152586e-11, 0.00256149, 0, 0, 0.00336262, 0, 0, 0.00440481, 0, 0, 0.00575752, 0, 0, 0.00750922, 0, 0, 0.00977233, 3.33759524e-14, 0, 0.01268935, 0, 0, 0.01644035, 0, 0, 0.02125231, 0, 0, 0.02741059, 0, 0, 0.03527286, 0, 0, 0.04528627, 1.33503809e-13, 0, 0.05800822, 3.33759524e-14, 0, 0.07413156, 1.66879762e-13, 0, 0.09451507, 0, 2.1281133e-10, 0.1202201, 0, 4.58314406e-10, 0.15255454, 6.67519047e-14, 1.21104643e-09, 0.19312534, 0, 2.36645932e-09, 0.24390104, 6.67519047e-14, 4.75600646e-09, 0.30728586, 2.67007619e-13, 8.65942422e-09, 0.38620732, 0, 1.50202466e-08, 0.48421916, 2.00255714e-13, 2.50722824e-08, 0.60562208, 6.67519047e-14, 3.99336595e-08, 0.75560452, 3.33759524e-13, 6.13770413e-08, 0.9404062, 2.67007619e-13, 9.08808492e-08, 1.16750737, 0, 1.30175826e-07, 1.44584671, 0, 1.81009805e-07, 1.78607108, 2.00255714e-13, 2.45783717e-07, 2.20082051, 4.00511428e-13, 3.27090741e-07, 2.70505173, 4.67263333e-13, 4.3093748e-07, 3.31640365, 9.34526666e-13, 5.645011e-07, 4.05560814, 2.03593309e-12, 7.414353e-07, 4.94694904, 4.07186619e-12, 9.79922233e-07, 6.01877242, 8.11035642e-12, 1.30569716e-06, 7.30405027, 1.5386314e-11, 1.75204846e-06, 8.84099946, 2.91372064e-11, 2.3635429e-06, 10.6737569, 5.38020352e-11, 3.19718768e-06, 12.85311094, 9.82588037e-11, 4.33350375e-06, 15.43728769, 1.76058149e-10, 5.88787258e-06, 18.4927898, 3.11264132e-10, 8.04129544e-06, 22.09528329, 5.41341259e-10, 0.00001107, 26.33052634, 9.27183957e-10, 0.00001543, 31.29533163, 1.56369674e-09, 0.00002176, 37.09855172, 2.59906051e-09, 0.00003105, 43.86207398, 4.2583877e-09, 0.00004465, 51.7218093, 6.88116182e-09, 0.00006444, 60.82865542, 1.09713514e-08, 0.00009288, 71.34941312, 1.72695313e-08, 0.00013318, 83.46763024, 2.68508869e-08, 0.00018953, 97.38434554, 4.12632573e-08, 0.00026739, 113.31870165, 6.27138609e-08, 0.00037405, 131.50839376, 9.43279197e-08, 0.00051944, 152.2099182, 1.40499342e-07, 0.00071731, 175.69858362, 2.07369419e-07, 0.00098697, 202.2682459, 3.0347184e-07, 0.00135553, 232.23072781, 4.4060252e-07, 0.00186092, 265.91488444, 6.34979338e-07, 0.00255563, 303.66527709, 9.08780921e-07, 0.00351167, 345.84042051, 1.29217259e-06, 0.00482716, 392.81057213, 1.82595674e-06, 0.00663574, 444.95503676, 2.56502452e-06, 0.00912121, 502.65896662, 3.58284322e-06, 0.01254065, 566.30964378, 4.97730118e-06, 0.0172621, 636.29224142, 6.87837002e-06, 0.02382463, 712.98507019, 9.45826518e-06, 0.03303184, 796.75432742, 0.00001295, 0.04609141, 887.94837929, 0.00001764, 0.06481254, 986.89161935, 0.00002395, 0.09186755, 1093.87796045, 0.00003242, 0.13110992, 1209.16403126, 0.00004377, 0.18792081, 1332.96216285, 0.00005902, 0.26953291, 1465.43326394, 0.00007954, 0.38526828, 1606.67969676, 0.00010729, 0.5466472, 1756.7382768, 0.00014498, 0.76740202, 1915.57352952, 0.00019651, 1.06358703, 2083.07134522, 0.00026741, 1.45421886, 2259.03317838, 0.00036564, 1.96318886, 2443.17094032, 0.00050259, 2.6234934, 2635.10273319, 0.0006946, 3.48501062, 2834.34956926, 0.00096501, 4.62695524, 3040.33321162, 0.00134711, 6.17561973, 3252.37526058, 0.00188808, 8.32703362, 3469.69759509, 0.00265452, 11.3729761, 3691.42425933, 0.00373986, 15.72796342, 3916.58486232, 0.00527443, 21.9553048, 4144.11953295, 0.00743897, 30.79291666, 4372.8854449, 0.01048271, 43.18440517, 4601.66489546, 0.01474745, 60.32661892, 4829.17489113, 0.02069964, 83.74841602, 5054.07816002, 0.02897304, 115.43296091, 5274.99547934, 0.04042573, 157.98504791, 5490.5191744, 0.05621621, 214.82764097, 5699.22761618, 0.07790557, 290.3963241, 5899.70051719, 0.10759477, 390.30063771, 6090.53480154, 0.14810919, 521.45212851, 6270.36080551, 0.20324626, 692.22797136, 6437.85854949, 0.27810559, 912.83909206, 6591.77381288, 0.37952491, 1196.17952586, 6730.93373884, 0.51664738, 1559.51928777, 6854.26169792, 0.70164652, 2027.44995045, 6960.79114731, 0.95063099, 2636.52124534, 7049.67823629, 1.28474316, 3442.08885869, 7120.21292786, 1.73144932, 4528.1356522, 7171.82843183, 2.3259966, 6021.33198945, 7204.10877415, 3.11298132, 8111.3966555, 7216.79436134, 4.14793937, 11080.80318858, 7209.78543635, 5.49883456, 15347.75216618, 7183.14336187, 7.24729195, 21526.56272398, 7137.08970842, 9.48940885, 30508.39361717, 7072.0031666, 12.3359844, 43561.32525132, 6988.41434396, 15.9120387, 62440.96718487, 6886.99854682, 20.35555871, 89489.89205958, 6768.56668469, 25.81548637, 127686.78357711, 6634.0544688, 32.4490633, 180587.64862757, 6484.51010597, 40.41873366, 252089.36971921, 6321.08071405, 49.88885664, 345951.58779591, 6144.99770474, 61.02249581, 465048.23008866, 5957.56139347, 73.97846916, 610391.04725796, 5760.12510397, 88.90875897, 780067.29976653, 5554.07903748, 105.95611217, 968337.27059851, 5340.83417232, 125.25157978, 1165203.80726963, 5121.8064507, 146.91139245, 1356750.8934761, 4898.40149479, 171.03259004, 1526422.52167348, 4672.0000753, 197.68679332, 1657184.58052985, 4443.94453256, 226.91156566, 1734235.54691925, 4215.52632373, 258.69927939, 1747698.03680531, 3987.97484125, 292.98364531, 1694630.06174663, 3762.44761683, 329.62490861, 1579802.54040474, 3540.021994, 368.39441138, 1414987.26804181, 3321.68832026, 408.96099389, 1216898.75295298, 3108.34468091, 450.88150949, 1004298.42177957, 2900.79316354, 493.59495038, 794973.12918851, 2699.73761804, 536.42645118, 603276.6521061, 2505.78285297, 578.59845997, 438696.45434601, 2319.43517548, 619.25113233, 305575.64729172, 2141.10419099, 657.47106714, 203805.86150866, 1971.10571691, 692.33047351, 130108.52208616, 1809.66571206, 722.93067275, 79477.4937625, 1656.92505713, 748.44549138, 46440.77834305, 1512.945045, 768.17062695, 25950.62040493, 1377.71348382, 781.5619254, 13863.49191619, 1251.15115735, 788.27656983, 7078.88511422, 1133.11870644, 788.1732901, 3453.99391366, 1023.42351408, 781.34850203, 1610.07387294, 921.82673613, 768.10324261, 716.87944996, 828.05008738, 748.94914653, 304.81495148, 741.78311024, 724.56478996, 123.74702514, 662.68885363, 695.7116279, 47.95841126, 590.41092018, 663.36929537, 17.74000945, 524.57831467, 628.29563591, 0, 464.81160344, 591.44614496, 0, 410.72740327, 553.70939575, 0, 361.94252168, 515.76840044, 0, 318.07805157, 478.11871241, 0, 278.7642647, 441.27687267, 0, 243.63800016, 405.37066657, 0, 212.35173217, 371.20092364, 0, 184.57575619, 338.19830335, 0, 159.98660387, 307.05836193, 0, 138.28629587, 277.65440488, 0, 119.2111415, 249.81767463, 0, 102.48417889, 223.40944888, 0, 87.83198148, 198.65219107, 0, 0, 175.64614811, 0, 0, 154.30801511, 0, 0, 134.20480339, 0, 0, 115.90590215, 0, 0, 99.06947277, 0, 0, 83.86417027, 0, 0, 70.21045431, 0, 0, 58.47211449, 0, 0, 47.78640658, 0, 0, 38.68015927, 0, 0, 31.08724468, 0, 0, 24.57176017, 0, 0, 19.05428931, 0, 0, 14.62743547, 0, 0, 11.09310835, 0, 0, 8.12557309, 0, 0, 5.78101079, 0, 0, 4.19072833, 0, 0, 2.98005447, 0, 0, 2.05841601, 0, 0, 1.41472616, 0, 0, 0.87556236, 0, 0, 0.56659694, 0, 0, 0.27623448, 0, 0, 0.13726909, 0, 0, 0.07535991, 0, 0, 0.00128593, 0, 0, 6.14151393e-14] @@ -141,6 +151,11 @@ "type": "character", "attributes": {}, "value": ["1/year"] + }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["average"] } }, "value": [1620.54514448, 3661.04629032, 2275.82527489, 1719.21961641, 3814.15780294, 2379.51209696, 1829.29259114, 3966.98027609, 2486.1931826, 1954.47042106, 4120.12833196, 2597.19648495, 2099.41827813, 4274.59597231, 2714.30377717, 2269.88402805, 4431.81229397, 2839.81273406, 2472.97119376, 4593.75318514, 2976.66805246, 2717.76271737, 4763.18470552, 3128.75494023, 3016.55960122, 4944.13730734, 3301.47718144, 3387.03271994, 5142.72292206, 3502.75782539, 3855.60732856, 5368.41492685, 3744.61037537, 4462.45947529, 5635.93345396, 4045.45598228, 5268.68001827, 5967.94487815, 4433.44385342, 6366.52902481, 6398.92235621, 4951.20203081, 7894.2835747, 6980.73248674, 5662.71432733, 10057.89896382, 7790.78276934, 6663.3511507, 13162.34182064, 8943.80465335, 8094.37756417, 17655.62334944, 10608.41079716, 10163.34050572, 24187.65464356, 13029.22440558, 13171.31751002, 33683.21782699, 16554.31507917, 17546.69988534, 47422.61081102, 21665.51878975, 23882.5276523, 67114.14430287, 29005.69285575, 32970.05115546, 94929.9731724, 39392.18264225, 45815.31596154, 133463.22807928, 53800.69433859, 63619.3098504, 185555.60981256, 73300.45788892, 87698.13483213, 253948.77500362, 98923.13075919, 119321.59605289, 340738.59801697, 131457.57812784, 159460.52424619, 446663.19814219, 171182.14471759, 208457.13322897, 570328.37048315, 217573.38885887, 265666.39597558, 707549.55089012, 269058.63610953, 329151.37478533, 851037.95108271, 322897.94898103, 395537.89854769, 990647.41545724, 375276.94122691, 460128.84803171, 1114306.87833193, 421657.39401238, 517335.86651884, 1209596.65256261, 457369.96776551, 561409.1568167, 1265724.86703646, 478357.37940268, 587352.54203482, 1275489.94738354, 481912.33153107, 591832.06338668, 1236747.07412674, 467228.9258368, 573854.92393592, 1152975.08152589, 435615.82153971, 535031.94518009, 1032757.28889475, 390301.00812818, 479337.18655989, 888280.82881356, 335867.51171399, 412413.14163994, 733225.23567781, 277459.4506837, 340593.17049236, 580559.38963572, 219953.63438451, 269881.50816155, 440748.95181286, 167285.51903186, 205123.33150142, 320711.3575943, 122056.26643563, 149520.94427661, 223613.10023525, 85456.52682781, 104539.94006539, 149374.93780154, 57456.42112185, 70143.13363235, 95606.1241888, 37156.84814097, 45223.12927311, 58656.29975793, 23185.03508756, 28089.36457965, 34535.44992601, 14040.48049773, 16894.0749195, 19562.92181558, 8339.01387773, 9932.80762432, 10717.41072855, 4944.41392995, 5806.49006271, 5738.16004316, 3006.46599055, 3468.28624791, 3062.76702522, 1937.58941131, 2194.82622575, 1685.98499243, 1359.73855506, 1520.81397027, 1002.5541025, 1045.30888728, 1166.26522618, 670.22238186, 865.62155379, 973.06208835, 506.77497105, 752.18374209, 857.22714566, 420.78073739, 670.68927904, 776.94707081, 368.8581496, 604.83004153, 712.51119238, 327.09752245, 545.67000685, 653.34633925, 299.53010676, 494.75881605, 601.77244557, 273.34665494, 447.2566772, 552.18295719, 248.61057089, 403.11451774, 504.8142764, 225.39486001, 362.29780848, 459.9230287, 203.64309231, 324.61898718, 417.46272323, 183.53033047, 290.17322795, 377.92257889, 164.77398167, 258.51283087, 340.72132496, 147.52016629, 229.70620907, 306.27511586, 131.64441792, 203.50574034, 274.36368455, 117.0174714, 179.66936905, 244.74755871, 103.5165623, 157.9584528, 217.2055183, 91.14182445, 138.28287915, 191.80186385, 59.29457291, 71.09792507, 138.15396045, 52.09125251, 62.46069051, 121.37051478, 45.30481645, 54.32332652, 105.55839282, 39.12747897, 46.91631006, 91.16544594, 33.44384229, 40.10127195, 77.92280201, 28.31084093, 33.94648023, 65.96311611, 23.70162367, 28.41973862, 55.22382604, 19.73899851, 23.66830164, 45.99106943, 16.13172051, 19.34294826, 37.58625734, 13.05763633, 15.65692784, 30.42376535, 10.49442254, 12.58347317, 24.45158074, 8.29492728, 9.94613991, 19.32684559, 6.43234115, 7.71278191, 14.98709514, 4.93792519, 5.92088311, 11.50516628, 3.74480812, 4.49026065, 8.7252517, 2.74302848, 3.2890638, 6.39114558, 1.95155185, 2.34003351, 4.54703701, 1.41470479, 1.69632008, 3.29620503, 1.00600587, 1.20626436, 2.34395306, 0.69487945, 0.83320419, 1.61904105, 0.47758283, 0.57265187, 1.1127487, 0.29557207, 0.35440952, 0.68867099, 0.19127162, 0.22934671, 0.44565515, 0.09325115, 0.111814, 0.21727142, 0.04633926, 0.05556369, 0.1079686, 0.02543998, 0.03050413, 0.05927412, 0.0004341, 0.00052052, 0.00101145, 2.07325039e-14, 2.48595772e-14, 4.8305897e-14] @@ -250,6 +265,11 @@ "type": "character", "attributes": {}, "value": ["1/year"] + }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["average"] } }, "value": [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5] @@ -301,6 +321,11 @@ "type": "character", "attributes": {}, "value": ["1/year"] + }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["average"] } }, "value": [1620.73220406, 3661.13276768, 2275.8428408, 1719.40667599, 3814.2442803, 2379.52966286, 1829.47965071, 3967.06675345, 2486.2107485, 1954.65748064, 4120.21480932, 2597.21405085, 2099.6053377, 4274.68244967, 2714.32134307, 2270.07108763, 4431.89877133, 2839.83029996, 2473.15825333, 4593.8396625, 2976.68561836, 2717.94977694, 4763.27118288, 3128.77250613, 3016.74666079, 4944.2237847, 3301.49474734, 3387.21977951, 5142.80939942, 3502.77539129, 3855.79438813, 5368.50140421, 3744.62794127, 4462.64653486, 5636.01993132, 4045.47354818, 5268.86707784, 5968.03135551, 4433.46141933, 6366.71608439, 6399.00883357, 4951.21959671, 7894.47063427, 6980.8189641, 5662.73189323, 10058.08602339, 7790.8692467, 6663.3687166, 13162.52888021, 8943.89113071, 8094.39513007, 17655.81040901, 10608.49727453, 10163.35807162, 24187.84170314, 13029.31088294, 13171.33507592, 33683.40488657, 16554.40155653, 17546.71745124, 47422.7978706, 21665.60526711, 23882.5452182, 67114.33136244, 29005.77933311, 32970.06872136, 94930.16023198, 39392.26911961, 45815.33352744, 133463.41513885, 53800.78081595, 63619.3274163, 185555.79687214, 73300.54436628, 87698.15239804, 253948.9620632, 98923.21723655, 119321.6136188, 340738.78507655, 131457.6646052, 159460.54181209, 446663.38520176, 171182.23119495, 208457.15079487, 570328.55754272, 217573.47533623, 265666.41354148, 707549.7379497, 269058.72258689, 329151.39235123, 851038.13814228, 322898.03545839, 395537.9161136, 990647.60251681, 375277.02770427, 460128.86559761, 1114307.06539151, 421657.48048974, 517335.88408474, 1209596.83962218, 457370.05424287, 561409.17438261, 1265725.05409603, 478357.46588004, 587352.55960072, 1275490.13444311, 481912.41800843, 591832.08095258, 1236747.26118632, 467229.01231416, 573854.94150182, 1152975.26858546, 435615.90801707, 535031.96274599, 1032757.47595433, 390301.09460554, 479337.20412579, 888281.01587313, 335867.59819135, 412413.15920584, 733225.42273738, 277459.53716106, 340593.18805826, 580559.57669529, 219953.72086187, 269881.52572745, 440749.13887244, 167285.60550922, 205123.34906732, 320711.54465387, 122056.35291299, 149520.96184251, 223613.28729482, 85456.61330517, 104539.95763129, 149375.12486111, 57456.50759921, 70143.15119825, 95606.31124837, 37156.93461833, 45223.14683901, 58656.4868175, 23185.12156492, 28089.38214555, 34535.63698559, 14040.56697509, 16894.0924854, 19563.10887516, 8339.10035509, 9932.82519022, 10717.59778813, 4944.50040731, 5806.50762861, 5738.34710273, 3006.55246791, 3468.30381381, 3062.9540848, 1937.67588867, 2194.84379165, 1686.17205201, 1359.82503242, 1520.83153617, 1003.24116207, 1045.39536464, 1166.28279208, 670.90944144, 865.70803115, 973.07965425, 507.46203062, 752.27021945, 857.24471156, 421.46779696, 670.7757564, 776.96463671, 369.54520917, 604.91651889, 712.52875828, 327.78458203, 545.75648421, 653.36390516, 300.21716633, 494.84529341, 601.79001147, 274.03371452, 447.34315456, 552.20052309, 249.29763046, 403.2009951, 504.8318423, 226.08191959, 362.38428584, 459.9405946, 204.33015188, 324.70546454, 417.48028913, 184.21739004, 290.25970531, 377.94014479, 165.46104124, 259.09930823, 340.73889086, 148.20722587, 230.29268643, 306.29268176, 132.3314775, 204.0922177, 274.38125045, 117.70453097, 180.25584641, 244.76512461, 104.20362187, 158.54493016, 217.2230842, 91.82888402, 138.86935651, 191.81942975, 59.98163249, 71.68440243, 138.17152635, 52.77831208, 63.04716787, 121.38808069, 45.99187602, 54.90980388, 105.57595872, 39.81453854, 47.50278742, 91.18301184, 34.13090186, 40.68774931, 77.94036791, 28.9979005, 34.53295759, 65.98068201, 24.38868325, 29.00621598, 55.24139194, 20.42605808, 24.254779, 46.00863533, 16.81878008, 19.92942562, 37.60382324, 13.7446959, 16.2434052, 30.94133125, 11.18148211, 13.16995053, 24.96914664, 8.98198685, 10.53261727, 19.84441149, 7.11940073, 8.29925927, 15.50466104, 5.62498476, 6.50736047, 12.02273218, 4.43186769, 5.07673801, 9.2428176, 3.43008805, 3.87554116, 6.90871148, 2.63861143, 2.92651087, 5.06460291, 2.10176436, 2.28279745, 3.81377093, 1.69306545, 1.79274172, 2.86151896, 1.38193902, 1.41968156, 2.13660695, 1.1646424, 1.15912923, 1.6303146, 0.98263164, 0.94088688, 1.20623689, 0.87833119, 0.81582407, 0.96322105, 0.78031072, 0.69829136, 0.73483732, 0.73339883, 0.64204105, 0.6255345, 0.71249955, 0.61698149, 0.57684002, 0.68749368, 0.58699788, 0.51857735, 0.68705957, 0.58647736, 0.5175659] @@ -316,6 +341,11 @@ "attributes": {}, "value": [3, 100] }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["point"] + }, "dimnames": { "type": "list", "attributes": { @@ -367,6 +397,11 @@ "attributes": {}, "value": [3, 100] }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["average"] + }, "dimnames": { "type": "list", "attributes": { @@ -446,6 +481,11 @@ "attributes": {}, "value": [3, 100] }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["point"] + }, "dimnames": { "type": "list", "attributes": { @@ -497,6 +537,11 @@ "attributes": {}, "value": [3, 100] }, + "representation": { + "type": "character", + "attributes": {}, + "value": ["point"] + }, "dimnames": { "type": "list", "attributes": { diff --git a/tests/testthat/test-ArraySpeciesBySize.R b/tests/testthat/test-ArraySpeciesBySize.R index 0fbe1977f..fdb9eb526 100644 --- a/tests/testthat/test-ArraySpeciesBySize.R +++ b/tests/testthat/test-ArraySpeciesBySize.R @@ -252,3 +252,60 @@ test_that("ArraySpeciesBySize value_name attribute", { mort <- getMort(NS_params_small) expect_identical(attr(mort, "value_name"), "Total mortality") }) + +# Second-order plotting placement (#382) ----------------------------------- + +test_that("producers tag the representation honestly", { + # Bin-averaged sinks are tagged "average"; flux/point quantities "point". + expect_identical(attr(getPredMort(NS_params_small), "representation"), + "average") + expect_identical(attr(getMort(NS_params_small), "representation"), + "average") + expect_identical(attr(getFMort(NS_params_small), "representation"), + "average") + expect_identical(attr(getExtMort(NS_params_small), "representation"), + "average") + expect_identical(attr(getEncounter(NS_params_small), "representation"), + "point") + expect_identical(attr(getEGrowth(NS_params_small), "representation"), + "point") +}) + +test_that("bin_midpoints are the geometric bin centres", { + p <- NS_params_small + beta <- p@w_full[2] / p@w_full[1] + expect_equal(bin_midpoints(p), p@w * sqrt(beta)) + expect_equal(bin_midpoints(p, w_full = TRUE), p@w_full * sqrt(beta)) +}) + +test_that("the size axis honours representation only under second_order_w", { + p <- NS_params_small + # Default first-order model: averaged quantities still plot at the nodes, + # so existing plots are unchanged. + expect_equal(get_ArraySpeciesBySize_w(getPredMort(p)), p@w) + expect_equal(get_ArraySpeciesBySize_w(getEncounter(p)), p@w) + + second_order_w(p) <- c(bin_average = TRUE) + # Now averaged quantities move to the geometric bin centres ... + expect_equal(get_ArraySpeciesBySize_w(getPredMort(p)), bin_midpoints(p)) + # ... while point-valued quantities stay on the nodes. + expect_equal(get_ArraySpeciesBySize_w(getEncounter(p)), p@w) +}) + +test_that("as.data.frame x-values follow the representation tag", { + p <- NS_params_small + second_order_w(p) <- c(bin_average = TRUE) + df_avg <- as.data.frame(getPredMort(p)) + df_pt <- as.data.frame(getEncounter(p)) + expect_equal(sort(unique(df_avg$w)), sort(unname(bin_midpoints(p)))) + expect_equal(sort(unique(df_pt$w)), sort(unname(p@w))) +}) + +test_that("subsetting preserves the representation tag", { + p <- NS_params_small + second_order_w(p) <- c(bin_average = TRUE) + mort <- getPredMort(p) + sub <- mort[1:2, ] + expect_identical(attr(sub, "representation"), "average") + expect_equal(get_ArraySpeciesBySize_w(sub), bin_midpoints(p)) +}) diff --git a/tests/testthat/test-ArrayTimeBySpeciesBySize.R b/tests/testthat/test-ArrayTimeBySpeciesBySize.R index 28a6c519e..a0711a1ea 100644 --- a/tests/testthat/test-ArrayTimeBySpeciesBySize.R +++ b/tests/testthat/test-ArrayTimeBySpeciesBySize.R @@ -250,3 +250,34 @@ test_that("animateSpectra is a backward-compatible alias for animate", { expect_s3_class(animateSpectra(fmort_small), "plotly") expect_s3_class(animateSpectra(NS_sim_small), "plotly") }) + +# Second-order plotting placement (#382) ----------------------------------- + +test_that("the time-series size axis honours representation under second_order_w", { + # getFMort.MizerSim tags its result "average". + expect_identical(attr(fmort_small, "representation"), "average") + p <- NS_sim_small@params + # Default first-order: size axis stays on the nodes (df unchanged). The + # stored dimnames are rounded to 3 sig figs, so compare with tolerance. + df_off <- as.data.frame(fmort_small) + expect_equal(sort(unique(df_off$w)), sort(unname(p@w)), + tolerance = 1e-2) + + # Switch the model to second order and recompute on the sim. + sim2 <- NS_sim_small + second_order_w(sim2@params) <- c(bin_average = TRUE) + fmort2 <- getFMort(sim2) + expect_equal(get_ArrayTimeBySpeciesBySize_w(fmort2), bin_midpoints(p)) + df_on <- as.data.frame(fmort2) + expect_equal(sort(unique(df_on$w)), sort(unname(bin_midpoints(p)))) +}) + +test_that("a time slice keeps the representation tag", { + sim2 <- NS_sim_small + second_order_w(sim2@params) <- c(bin_average = TRUE) + fmort2 <- getFMort(sim2) + slice <- ArrayTimeBySpeciesBySize_slice(fmort2) + expect_s3_class(slice, "ArraySpeciesBySize") + expect_identical(attr(slice, "representation"), "average") + expect_equal(get_ArraySpeciesBySize_w(slice), bin_midpoints(sim2@params)) +}) diff --git a/tests/testthat/test-project_methods.R b/tests/testthat/test-project_methods.R index e6bc1c26c..2fcb395e8 100644 --- a/tests/testthat/test-project_methods.R +++ b/tests/testthat/test-project_methods.R @@ -213,7 +213,8 @@ test_that("species-size rate getters work for MizerSim", { test_that("getCriticalFeedingLevel matches metab over intake_max times alpha", { expected <- params@metab / params@intake_max / params@species_params$alpha expect_equal(getCriticalFeedingLevel(params), expected, - ignore_attr = c("value_name", "units", "class", "params")) + ignore_attr = c("value_name", "units", "class", "params", + "representation")) }) # getPredRate ------------------------------------------------------------- diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index 9caf3e7bc..56599edb5 100644 --- a/vignettes/numerical_details.qmd +++ b/vignettes/numerical_details.qmd @@ -106,6 +106,8 @@ In summary: The default scheme uses point values throughout, which is consistent. Replacing the bin-integrated quantities by their bin averages raises the order of those source and sink terms; the spatial order of the overall transport step is separately set by the upwind advective flux, as discussed in the Order of Accuracy section below. +**Plotting follows the same distinction.** A bin average $N_j$ does not live at the bin boundary $w_j$ but at the geometric bin centre $w^*_j=\sqrt{w_j\,w_{j+1}}=w_j\sqrt\beta$ (the log-midpoint, exact for the community spectrum $N\propto w^{-2}$). So under second-order bin-averaging mizer draws bin-averaged quantities (the abundance and the mortality/reproduction sinks) at $w^*_j$ — a uniform half-bin shift to the right on the log axis — while point-valued quantities (the encounter and growth-type rates) stay on the nodes $w_j$. The size-resolved array classes carry a `representation` tag recording which a quantity is, and the shift is applied only when `second_order_w[["bin_average"]]` is set, so default plots are unchanged. + ## Semi-Implicit Time Discretisation With the diffusion term, an explicit time discretisation would require a very small time step for stability ($\Delta t \sim \Delta w^2$). Therefore, we use a semi-implicit scheme where the densities $N$ are evaluated at time $t+1$, but the rates ($g$, $\mu$, $d$) are evaluated at time $t$. Using a fully implicit scheme would require solving a nonlinear system at each time step, which is more computationally expensive. From 330e7e9d5c44d4a367b123ce23ac6f2eb1e43511 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Thu, 11 Jun 2026 18:41:01 +0100 Subject: [PATCH 2/3] Second-order power weighting in plotSpectra/plotCDF (#383) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Refines the plotting placement of #382 for the `power` argument. For a cell average N_j the transformed density N(w) w^power is represented at the geometric bin centre w* = w sqrt(beta) by N_j (w*)^power, so using the node for both the w^power weight and the x-location is doubly off -- the location error is the half-bin shift, and the w^power error is O(power*Δx), largest for the common power = 2 Sheldon plot. plot_spectra() now evaluates both the w^power weight and the plotted location at the bin centre (species, resource, total, background) under second_order_w[["bin_average"]]; the default auto wlim follows the plotted grid so the shift does not clip the end bins. animate.MizerSim() gets the same treatment. plotlySpectra/plotlyCDF/plotSpectraRelative inherit it (the latter's w^power factor cancels in the relative ratio, so only its x-location shifts). plotCDF is the opposite case: a CDF value is cumulative up to a size, a boundary quantity, so plot_cdf() maps the density back to the bin edges (w = w*/sqrt(beta)) -- keeping the cumulative on the nodes -- while the increments stay centre-weighted, so each value*dw_k is the second-order bin integral of N w^power. getFlux() (a face/point quantity) is deliberately left untouched. Default (first-order) plots, return_data and vdiffr snapshots are unchanged. Tests cover power = 0/1/2 placement and weighting, the CDF staying on edges, and the relative-plot cancellation. Full suite: 2842 pass. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 11 ++++++ R/animateSpectra.R | 8 ++++ R/plots.R | 42 ++++++++++++++++---- tests/testthat/test-plots.R | 70 +++++++++++++++++++++++++++++++++ vignettes/numerical_details.qmd | 2 + 5 files changed, 125 insertions(+), 8 deletions(-) diff --git a/NEWS.md b/NEWS.md index 8dabcf37f..7023de650 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,16 @@ # mizer (development version) +- Under second-order bin-averaging (`second_order_w[["bin_average"]]`), the + spectrum plots `plotSpectra()`, `plotlySpectra()`, `plotSpectraRelative()` and + `animate()`/`animateSpectra()` now evaluate the `w^power` weight *and* the + marker location at the geometric bin centre `w* = w sqrt(beta)`, so each marker + is a point `(w*, N_j (w*)^power)` on the continuous `N w^power` curve rather + than being doubly mis-placed at the bin edge (the location error grows with + `power`, worst for the common `power = 2` Sheldon plot). `plotCDF()` / + `plotlyCDF()`, being cumulative integrals, keep their increments + bin-averaged but plot the cumulative on the bin **edges**, not the centres. + The default (first-order) behaviour is unchanged. (#383) + - 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 diff --git a/R/animateSpectra.R b/R/animateSpectra.R index 28b95e6e4..698bb8489 100644 --- a/R/animateSpectra.R +++ b/R/animateSpectra.R @@ -190,6 +190,14 @@ animate.MizerSim <- function(x, species = NULL, } else { y_label <- paste0("Number density * w^", power) } + # The animated spectrum is a density, so under second-order bin-averaging + # we evaluate both the w^power weight and the plotted location at the + # geometric bin centre w* = w sqrt(beta) (issue #383), matching + # plotSpectra(). Default plots are unchanged. + if (isTRUE(sim@params@second_order_w[["bin_average"]])) { + beta <- sim@params@w_full[2] / sim@params@w_full[1] + nf$w <- nf$w * sqrt(beta) + } nf <- mutate(nf, value = value * w^power) animate_plotly(nf, sim@params, log_x, log_y, y_label, wlim, llim, diff --git a/R/plots.R b/R/plots.R index ca6ff584c..4ccb75c87 100644 --- a/R/plots.R +++ b/R/plots.R @@ -1373,11 +1373,24 @@ plot_spectra <- function(params, n, n_pp, highlight, log_x, log_y, size_axis, return_data) { params <- validParams(params) size_axis <- plot_size_axis(size_axis) + # The spectrum is a density: each plotted value N_j w^power is a point on the + # continuous N(w) w^power curve, which (for the cell-average N_j) it touches + # at the geometric bin centre. Under second-order bin-averaging we therefore + # evaluate *both* the w^power weight and the plotted location at the bin + # centre w* = w sqrt(beta) (issue #383). The default keeps the grid nodes, so + # existing spectrum plots are unchanged. + second_order <- isTRUE(params@second_order_w[["bin_average"]]) + w_grid <- if (second_order) bin_midpoints(params) else params@w + w_full_grid <- if (second_order) bin_midpoints(params, w_full = TRUE) else + params@w_full + # Default size limits follow the plotted grid (the bin centres under + # second order) so the centre shift does not clip the top/bottom bins. With + # the default (first-order) nodes these reduce to the previous limits. if (is.na(wlim[1])) { - wlim[1] <- if (resource) min(params@w) / 100 else min(params@w) + wlim[1] <- if (resource) min(w_grid) / 100 else min(w_grid) } if (is.na(wlim[2])) { - wlim[2] <- max(params@w_full) + wlim[2] <- max(w_full_grid) } if (total) { @@ -1385,16 +1398,16 @@ plot_spectra <- function(params, n, n_pp, fish_idx <- (length(params@w_full) - length(params@w) + 1):length(params@w_full) total_n <- n_pp total_n[fish_idx] <- total_n[fish_idx] + colSums(n) - total_n <- total_n * params@w_full^power + total_n <- total_n * w_full_grid^power } species <- valid_species_arg(params, species) # Deal with power argument y_label <- spectra_y_label(power) - n <- sweep(n, 2, params@w^power, "*") + n <- sweep(n, 2, w_grid^power, "*") # Select only the desired species spec_n <- n[as.character(dimnames(n)[[1]]) %in% species, , drop = FALSE] # Make data.frame for plot - plot_dat <- data.frame(w = rep(params@w, + plot_dat <- data.frame(w = rep(w_grid, each = dim(spec_n)[[1]]), value = c(spec_n), Species = dimnames(spec_n)[[1]], @@ -1404,7 +1417,7 @@ plot_spectra <- function(params, n, n_pp, (params@w_full <= wlim[2]) # Do we have any resource to plot? if (sum(resource_sel) > 0) { - w_resource <- params@w_full[resource_sel] + w_resource <- w_full_grid[resource_sel] plank_n <- n_pp[resource_sel] * w_resource^power plot_dat <- rbind(plot_dat, data.frame(w = w_resource, @@ -1416,7 +1429,7 @@ plot_spectra <- function(params, n, n_pp, } if (total) { plot_dat <- rbind(plot_dat, - data.frame(w = params@w_full, + data.frame(w = w_full_grid, value = c(total_n), Species = "Total", Legend = "Total") @@ -1426,7 +1439,7 @@ plot_spectra <- function(params, n, n_pp, back_n <- n[params@species_params$is_background, , drop = FALSE] plot_dat <- rbind(plot_dat, - data.frame(w = rep(params@w, + data.frame(w = rep(w_grid, each = dim(back_n)[[1]]), value = c(back_n), Species = as.factor(dimnames(back_n)[[1]]), @@ -1660,6 +1673,19 @@ plotCDF.MizerParams <- function(object, species = NULL, plot_cdf <- function(plot_dat, params, power, normalise, log_x, log_y, wlim, llim, ylim, highlight, size_axis, return_data) { size_axis <- plot_size_axis(size_axis) + # A CDF value is cumulative *up to a size* — a boundary quantity — so it + # belongs on the bin edges, not the bin centres. Under second-order + # bin-averaging `plot_spectra()` returns its density on the geometric bin + # centres (issue #383); for the CDF we map those back to the bin edges + # (w = w* / sqrt(beta)) so the cumulative stays on the grid nodes. The + # density *values* are kept centre-weighted, so each increment + # value * dw_k = N_k (w*_k)^power dw_k is still the second-order bin integral + # of N w^power. (See get_ArraySpeciesBySize_w() and the FFT/numerical-details + # vignettes.) + if (isTRUE(params@second_order_w[["bin_average"]])) { + beta <- params@w_full[2] / params@w_full[1] + plot_dat$w <- plot_dat$w / sqrt(beta) + } if (identical(size_axis, "l")) { plot_dat_l <- convert_plot_size_axis(plot_dat, params, size_axis, drop_w = FALSE) diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index fc7caab0d..f80cfaa43 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -677,3 +677,73 @@ test_that("plotDiet works with MizerSim", { p <- plotDiet(sim, species = 2, time_range = 1:2) expect_true(is(p, "ggplot")) }) + +# Second-order power weighting in plotSpectra / plotCDF (#383) -------------- + +test_that("plotSpectra draws the spectrum at bin centres with the w^power weight there", { + p0 <- params + p1 <- params + second_order_w(p1) <- c(bin_average = TRUE) + beta <- p0@w_full[2] / p0@w_full[1] + for (pw in c(0, 1, 2)) { + d0 <- plotSpectra(p0, power = pw, return_data = TRUE) + d1 <- plotSpectra(p1, power = pw, return_data = TRUE) + d0 <- d0[order(d0$Species, d0$w), ] + d1 <- d1[order(d1$Species, d1$w), ] + expect_equal(nrow(d0), nrow(d1)) + # x moves to the geometric bin centre (a uniform sqrt(beta) shift) ... + expect_equal(unname(d1$w / d0$w), rep(sqrt(beta), nrow(d1))) + # ... and the w^power weight is evaluated there, scaling the value + # (column 2, named by the y-label) by (w*/w)^power = beta^(power/2). + expect_equal(unname(d1[[2]] / d0[[2]]), + rep(beta^(pw / 2), nrow(d1))) + } +}) + +test_that("plotSpectra default (first order) is unchanged", { + expect_identical(plotSpectra(params, power = 2, return_data = TRUE), + plotSpectra(params, power = 2, return_data = TRUE)) + # The default model never shifts: x stays on the model grid nodes. + d <- plotSpectra(params, power = 2, resource = FALSE, total = FALSE, + background = FALSE, return_data = TRUE) + expect_true(all(d$w %in% params@w)) +}) + +test_that("plotCDF keeps the cumulative on bin edges under second_order_w", { + p1 <- params + second_order_w(p1) <- c(bin_average = TRUE) + cdf <- plotCDF(p1, power = 2, return_data = TRUE) + nodes <- round(c(p1@w, p1@w_full), 6) + centres <- round(c(bin_midpoints(p1), bin_midpoints(p1, w_full = TRUE)), 6) + # The CDF x-values stay on the node (bin-edge) grid, never on the centres. + expect_true(all(round(cdf$w, 6) %in% nodes)) + expect_false(any(round(cdf$w, 6) %in% setdiff(centres, nodes))) + # Cumulative is monotonic increasing and ends at 1 (normalised). + sp1 <- cdf[cdf$Species == p1@species_params$species[1], ] + sp1 <- sp1[order(sp1$w), ] + expect_true(all(diff(sp1[[2]]) >= -1e-12)) +}) + +test_that("plotSpectraRelative shifts x to centres but cancels the power weight", { + p1a <- params + p1b <- params + p1b@initial_n <- p1b@initial_n * 1.5 + second_order_w(p1a) <- c(bin_average = TRUE) + second_order_w(p1b) <- c(bin_average = TRUE) + beta <- params@w_full[2] / params@w_full[1] + # power cancels in 2(N2-N1)/(N1+N2), so the relative value is independent of + # power; only the x-location picks up the centre shift. + p_p1 <- plotSpectraRelative(p1a, p1b, species = species, resource = FALSE, + power = 1) + p_p2 <- plotSpectraRelative(p1a, p1b, species = species, resource = FALSE, + power = 2) + d_p1 <- p_p1$data[order(p_p1$data$Species, p_p1$data$w), ] + d_p2 <- p_p2$data[order(p_p2$data$Species, p_p2$data$w), ] + expect_equal(d_p1$w, d_p2$w) + expect_equal(d_p1$rel_diff, d_p2$rel_diff) + # The x-location is the geometric bin centre, not the node. + p_node <- plotSpectraRelative(params, params, species = species, + resource = FALSE) + expect_true(all(p_node$data$w %in% params@w)) + expect_false(all(d_p1$w %in% params@w)) +}) diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index 56599edb5..15d528aa5 100644 --- a/vignettes/numerical_details.qmd +++ b/vignettes/numerical_details.qmd @@ -108,6 +108,8 @@ The default scheme uses point values throughout, which is consistent. Replacing **Plotting follows the same distinction.** A bin average $N_j$ does not live at the bin boundary $w_j$ but at the geometric bin centre $w^*_j=\sqrt{w_j\,w_{j+1}}=w_j\sqrt\beta$ (the log-midpoint, exact for the community spectrum $N\propto w^{-2}$). So under second-order bin-averaging mizer draws bin-averaged quantities (the abundance and the mortality/reproduction sinks) at $w^*_j$ — a uniform half-bin shift to the right on the log axis — while point-valued quantities (the encounter and growth-type rates) stay on the nodes $w_j$. The size-resolved array classes carry a `representation` tag recording which a quantity is, and the shift is applied only when `second_order_w[["bin_average"]]` is set, so default plots are unchanged. +For the `power`-weighted spectrum plots (`plotSpectra()` and friends) the $w^{\text{power}}$ factor must be evaluated where the density value lives, so it too is taken at the bin centre: each marker is the point $\bigl(w^*_j,\,N_j\,(w^*_j)^{\text{power}}\bigr)$ on the continuous $N(w)\,w^{\text{power}}$ curve. (Sampling the weight at the edge would mis-scale it by a factor $\beta^{\text{power}/2}$, largest for the common $\text{power}=2$ Sheldon plot.) A cumulative plot (`plotCDF()`) is the opposite case: a CDF value is cumulative *up to a size*, a boundary quantity, so its increments use the bin-averaged (centre-weighted) density but the cumulative is plotted on the bin **edges**, not the centres. + ## Semi-Implicit Time Discretisation With the diffusion term, an explicit time discretisation would require a very small time step for stability ($\Delta t \sim \Delta w^2$). Therefore, we use a semi-implicit scheme where the densities $N$ are evaluated at time $t+1$, but the rates ($g$, $\mu$, $d$) are evaluated at time $t$. Using a fully implicit scheme would require solving a nonlinear system at each time step, which is more computationally expensive. From a934b3c9816dc6af634bb039898362a26de6e91b Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Mon, 22 Jun 2026 17:33:13 +0100 Subject: [PATCH 3/3] Place plotCDF cumulative on upper bin edges (#383) The cumulative sum in plotCDF is inclusive: the sum through bin k is the integral of N w^power over all bins up to and including bin k, i.e. the CDF at that bin's upper edge w_k + dw_k, not its left edge. It was plotted at the left edge, a one-bin (O(dx)) location offset that left the CDF only first-order accurate in placement. Make the inclusive convention explicit and place each cumulative value on its bin's upper edge, in both the default and second-order schemes. Under second-order bin-averaging the CDF is then second-order accurate in its placement as well as its increments. Tests: assert the upper-edge placement in both schemes; add an analytic regression test using the community spectrum N = C w^-2 with power = 2 (integrand constant, exact CDF linear) that pins down the placement and catches the one-bin offset. Relax the length-limited CDF bound, since the top in-window point (where the normalised CDF reaches 1) now legitimately sits on the upper edge, up to one bin above llim[2]. Co-Authored-By: Claude Opus 4.8 --- NEWS.md | 15 ++++++--- R/plots.R | 29 +++++++++++++---- tests/testthat/test-plots.R | 58 +++++++++++++++++++++++++++++---- vignettes/numerical_details.qmd | 2 +- 4 files changed, 85 insertions(+), 19 deletions(-) diff --git a/NEWS.md b/NEWS.md index 7023de650..07d2517c0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -6,10 +6,17 @@ marker location at the geometric bin centre `w* = w sqrt(beta)`, so each marker is a point `(w*, N_j (w*)^power)` on the continuous `N w^power` curve rather than being doubly mis-placed at the bin edge (the location error grows with - `power`, worst for the common `power = 2` Sheldon plot). `plotCDF()` / - `plotlyCDF()`, being cumulative integrals, keep their increments - bin-averaged but plot the cumulative on the bin **edges**, not the centres. - The default (first-order) behaviour is unchanged. (#383) + `power`, worst for the common `power = 2` Sheldon plot). These + spectrum-density changes apply only under second-order bin-averaging, so + default spectrum plots are unchanged. (#383) + +- `plotCDF()` / `plotlyCDF()` now plot each cumulative value on its bin's + **upper** edge `w_k + dw_k`, following the inclusive cumulative-sum + convention (the sum through bin `k` is the integral up to that bin's upper + edge). This corrects a long-standing one-bin location offset and applies in + both the default and the second-order schemes; under second-order + bin-averaging the CDF is then second-order accurate in its placement as well + as its increments. (#383) - Size-resolved diagnostics that are finite-volume bin averages (the mortalities `getPredMort()`, `getFMort()`, `getMort()`, `getExtMort()`, and diff --git a/R/plots.R b/R/plots.R index 4ccb75c87..4de5704fb 100644 --- a/R/plots.R +++ b/R/plots.R @@ -1676,12 +1676,14 @@ plot_cdf <- function(plot_dat, params, power, normalise, log_x, log_y, wlim, lli # A CDF value is cumulative *up to a size* — a boundary quantity — so it # belongs on the bin edges, not the bin centres. Under second-order # bin-averaging `plot_spectra()` returns its density on the geometric bin - # centres (issue #383); for the CDF we map those back to the bin edges - # (w = w* / sqrt(beta)) so the cumulative stays on the grid nodes. The - # density *values* are kept centre-weighted, so each increment - # value * dw_k = N_k (w*_k)^power dw_k is still the second-order bin integral - # of N w^power. (See get_ArraySpeciesBySize_w() and the FFT/numerical-details - # vignettes.) + # centres (issue #383). Here we map those back to the grid nodes + # (w = w* / sqrt(beta)) so that `prepare_spectra_cdf_data()` can look up the + # bin widths and form the bin integrals; that helper then places the + # cumulative on the *upper* bin edges (see its comments for the inclusive + # convention). The density *values* are kept centre-weighted, so each + # increment value * dw_k = N_k (w*_k)^power dw_k is the second-order bin + # integral of N w^power. (See get_ArraySpeciesBySize_w() and the + # FFT/numerical-details vignettes.) if (isTRUE(params@second_order_w[["bin_average"]])) { beta <- params@w_full[2] / params@w_full[1] plot_dat$w <- plot_dat$w / sqrt(beta) @@ -1727,13 +1729,26 @@ prepare_spectra_cdf_data <- function(plot_dat, params, normalise = TRUE) { params <- validParams(params) y_var <- names(plot_dat)[2] plot_dat <- plot_dat[order(plot_dat$Species, plot_dat$w), ] - plot_dat[[y_var]] <- plot_dat[[y_var]] * spectra_bin_width(plot_dat$w, params) + # Each bin contributes its integral, value * dw_k, which under second-order + # bin-averaging is N_k (w*_k)^power dw_k, the second-order bin integral of + # N w^power. + widths <- spectra_bin_width(plot_dat$w, params) + plot_dat[[y_var]] <- plot_dat[[y_var]] * widths + # The cumulative sum is *inclusive*: the sum through bin k integrates + # N w^power over all bins up to and including bin k, so it equals the CDF at + # that bin's *upper* edge w_k + dw_k, not at its left edge w_k. plot_dat[[y_var]] <- stats::ave(plot_dat[[y_var]], plot_dat$Species, FUN = cumsum) if (normalise) { totals <- stats::ave(plot_dat[[y_var]], plot_dat$Species, FUN = max) plot_dat[[y_var]] <- plot_dat[[y_var]] / totals } + # Make the inclusive convention explicit by placing each cumulative value at + # its bin's upper edge w_k + dw_k, the size up to which it accumulates. This + # corrects the historical one-bin (O(dx)) location offset (the inclusive sum + # was previously plotted at the left edge) and applies in both the default + # and the second-order schemes. + plot_dat$w <- plot_dat$w + widths plot_dat } diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R index f80cfaa43..38d81c7b0 100644 --- a/tests/testthat/test-plots.R +++ b/tests/testthat/test-plots.R @@ -385,7 +385,12 @@ test_that("size-based plots support length axes", { size_axis = "l", llim = llim, return_data = TRUE) expect_true(all(cdf_l_limited$l >= llim[1])) - expect_true(all(cdf_l_limited$l <= llim[2])) + # The cumulative now sits on each bin's upper edge (#383), so the top + # in-window point — where the normalised CDF reaches 1 — may exceed llim[2] + # by up to one bin, but no more (one weight-bin ratio is a safe length bound + # since the length ratio beta^(1/b) < beta). + beta_w <- params_len@w_full[2] / params_len@w_full[1] + expect_true(all(cdf_l_limited$l <= llim[2] * beta_w)) y_var <- names(cdf_l_limited)[2] expect_equal(stats::aggregate(cdf_l_limited[[y_var]] ~ cdf_l_limited$Species, FUN = max)[[2]], @@ -709,21 +714,60 @@ test_that("plotSpectra default (first order) is unchanged", { expect_true(all(d$w %in% params@w)) }) -test_that("plotCDF keeps the cumulative on bin edges under second_order_w", { +test_that("plotCDF places the cumulative on upper bin edges (inclusive convention)", { + upper_edges <- round(params@w_full + params@dw_full, 6) + # The inclusive cumulative sum belongs on each bin's *upper* edge w_k + dw_k, + # in both the default and the second-order schemes. + cdf0 <- plotCDF(params, power = 2, return_data = TRUE) + expect_true(all(round(cdf0$w, 6) %in% upper_edges)) + # The left bin edges (nodes) are no longer used for the placement. + nodes <- round(params@w, 6) + expect_false(all(round(cdf0$w, 6) %in% nodes)) + p1 <- params second_order_w(p1) <- c(bin_average = TRUE) cdf <- plotCDF(p1, power = 2, return_data = TRUE) - nodes <- round(c(p1@w, p1@w_full), 6) centres <- round(c(bin_midpoints(p1), bin_midpoints(p1, w_full = TRUE)), 6) - # The CDF x-values stay on the node (bin-edge) grid, never on the centres. - expect_true(all(round(cdf$w, 6) %in% nodes)) - expect_false(any(round(cdf$w, 6) %in% setdiff(centres, nodes))) - # Cumulative is monotonic increasing and ends at 1 (normalised). + expect_true(all(round(cdf$w, 6) %in% upper_edges)) + # never on the geometric bin centres + expect_false(any(round(cdf$w, 6) %in% setdiff(centres, upper_edges))) + # Cumulative is monotonic increasing per species. sp1 <- cdf[cdf$Species == p1@species_params$species[1], ] sp1 <- sp1[order(sp1$w), ] expect_true(all(diff(sp1[[2]]) >= -1e-12)) }) +test_that("plotCDF cumulative is second-order and free of the one-bin offset", { + # For a community spectrum N = C w^-2 with power = 2 the integrand N w^power + # is constant, so the exact CDF is linear in w: F(w) = V (w - w_min). The + # bin-average of C w^-2 is exact at the geometric centre, so the plotted + # density value is the constant V at every bin. We can therefore check the + # cumulative analytically, which pins down the edge placement: a one-bin + # offset (plotting F at the *left* edge) would break F(w_min-bin) != 0 and + # the slope, so this is a sharp regression test for issue #383. + p <- params + second_order_w(p) <- c(bin_average = TRUE) + sp <- 1 + p@initial_n[sp, ] <- p@w^(-2) # community spectrum N proportional to w^-2 + + # The weighted density is the constant V = beta at every bin centre. + dens <- plotSpectra(p, species = sp, power = 2, resource = FALSE, + total = FALSE, background = FALSE, return_data = TRUE) + V <- mean(dens[[2]]) + expect_equal(unname(dens[[2]]), rep(V, nrow(dens))) + + cdf <- plotCDF(p, species = sp, power = 2, resource = FALSE, total = FALSE, + background = FALSE, normalise = FALSE, return_data = TRUE) + cdf <- cdf[order(cdf$w), ] + w_min <- min(p@w) + # F plotted at the upper edge equals the exact integral V (w_upper - w_min). + expect_equal(unname(cdf[[2]]), V * (cdf$w - w_min), tolerance = 1e-8) + # The smallest plotted x is the upper edge of the first bin, and its value is + # one full bin integral (not zero, and not two bins) — i.e. no offset. + expect_equal(cdf$w[1], unname(p@w[1] + p@dw_full[match(p@w[1], p@w_full)])) + expect_equal(cdf[[2]][1], V * (cdf$w[1] - w_min), tolerance = 1e-8) +}) + test_that("plotSpectraRelative shifts x to centres but cancels the power weight", { p1a <- params p1b <- params diff --git a/vignettes/numerical_details.qmd b/vignettes/numerical_details.qmd index 15d528aa5..7dd85afd7 100644 --- a/vignettes/numerical_details.qmd +++ b/vignettes/numerical_details.qmd @@ -108,7 +108,7 @@ The default scheme uses point values throughout, which is consistent. Replacing **Plotting follows the same distinction.** A bin average $N_j$ does not live at the bin boundary $w_j$ but at the geometric bin centre $w^*_j=\sqrt{w_j\,w_{j+1}}=w_j\sqrt\beta$ (the log-midpoint, exact for the community spectrum $N\propto w^{-2}$). So under second-order bin-averaging mizer draws bin-averaged quantities (the abundance and the mortality/reproduction sinks) at $w^*_j$ — a uniform half-bin shift to the right on the log axis — while point-valued quantities (the encounter and growth-type rates) stay on the nodes $w_j$. The size-resolved array classes carry a `representation` tag recording which a quantity is, and the shift is applied only when `second_order_w[["bin_average"]]` is set, so default plots are unchanged. -For the `power`-weighted spectrum plots (`plotSpectra()` and friends) the $w^{\text{power}}$ factor must be evaluated where the density value lives, so it too is taken at the bin centre: each marker is the point $\bigl(w^*_j,\,N_j\,(w^*_j)^{\text{power}}\bigr)$ on the continuous $N(w)\,w^{\text{power}}$ curve. (Sampling the weight at the edge would mis-scale it by a factor $\beta^{\text{power}/2}$, largest for the common $\text{power}=2$ Sheldon plot.) A cumulative plot (`plotCDF()`) is the opposite case: a CDF value is cumulative *up to a size*, a boundary quantity, so its increments use the bin-averaged (centre-weighted) density but the cumulative is plotted on the bin **edges**, not the centres. +For the `power`-weighted spectrum plots (`plotSpectra()` and friends) the $w^{\text{power}}$ factor must be evaluated where the density value lives, so it too is taken at the bin centre: each marker is the point $\bigl(w^*_j,\,N_j\,(w^*_j)^{\text{power}}\bigr)$ on the continuous $N(w)\,w^{\text{power}}$ curve. (Sampling the weight at the edge would mis-scale it by a factor $\beta^{\text{power}/2}$, largest for the common $\text{power}=2$ Sheldon plot.) A cumulative plot (`plotCDF()`) is the opposite case: a CDF value is cumulative *up to a size*, a boundary quantity, so its increments use the bin-averaged (centre-weighted) density but the cumulative is plotted on the bin **edges**, not the centres. Because the cumulative sum is inclusive — the sum through bin $k$ is the integral over all bins up to and including bin $k$ — each cumulative value is placed on that bin's *upper* edge $w_k+\Delta w_k$ (in both the default and second-order schemes). This makes the inclusive convention explicit and removes a one-bin offset that would otherwise leave the CDF only first-order accurate in its placement. ## Semi-Implicit Time Discretisation