Second-order power weighting in plotSpectra/plotCDF (#383)

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

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,

R/plots.R

@@ -1373,28 +1373,41 @@ 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) {
         # Calculate total community abundance
         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)

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))
+})

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.