From a2a35b18255453631815d617e9a3b6366cc3c35b Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Sun, 16 Aug 2026 00:07:14 +0100
Subject: [PATCH 1/4] Sum the total on the axis it is plotted against
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

A total can only be formed once every line sits on the same coordinate.
On a weight axis they do: every species shares the model's weight grid.
On a length axis they do not, because each species — and, since the
previous commit, the resource — converts weight to length with its own
allometric relationship, so at a given length the lines sit at different
weights. That is why the `Total` line was dropped from length-based
plots.

It is now formed after the conversion, at equal length rather than at
equal weight, by `add_total_line()`: each series is interpolated onto the
sorted union of all the size coordinates, linearly in log size, and
contributes nothing outside its own range. Where the series already
share a grid — always on a weight axis, and on a length axis whenever the
weight-length parameters agree — the union is that grid and the
interpolation reproduces the values exactly, so nothing is approximated
in the cases where nothing needs to be. The weight-axis total is
unchanged, verified against the previous implementation for powers 0, 1
and 2.

This moves the total out of `prepare_ArraySpeciesBySize_plot_data()` and
past the conversion in every path that draws one: the array plot,
`addPlot()`, `plot2()` and `plotRelative()` via the two data-frame
helpers, `plotSpectra()`, and both `animate()` methods.

Two corrections come with it, both on the weight axis:

- `plotSpectra(total = TRUE, resource = FALSE)` now leaves the resource
  out of the total. The argument used to control only whether the
  resource was drawn, so a consumers-only total was unavailable.
- The total of a mizer array is now the sum of the lines drawn rather
  than of the whole array, so it no longer counts abundance outside a
  species' size range. Use `all.sizes = TRUE` for the old sum.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 NEWS.md                                       |  19 +++
 R/ArraySpeciesBySize-class.R                  |  31 ++--
 R/ArrayTimeBySpeciesBySize-class.R            |  18 +--
 R/animateSpectra.R                            |  49 ++++--
 R/plots.R                                     | 146 +++++++++++++++---
 inst/skills/upgrade-mizer-code/SKILL.md       |  32 ++++
 man/add_total_line.Rd                         |  48 ++++++
 man/plotComparisonDataFrame.Rd                |   6 +-
 man/plotRelativeDataFrame.Rd                  |   3 +-
 .../testthat/test-ArraySpeciesBySize-class.R  |  49 ++++++
 tests/testthat/test-plots.R                   |  96 +++++++++++-
 vignettes/upgrading.Rmd                       |  29 ++++
 12 files changed, 453 insertions(+), 73 deletions(-)
 create mode 100644 man/add_total_line.Rd

diff --git a/NEWS.md b/NEWS.md
index 8561e01e2..574a925ff 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -159,6 +159,25 @@ stability of steady states.
   copepod really is shorter than a 1 mg fish larva — but it does mean the
   resource and the species are measured differently.
 
+- The total line is now shown on length-based plots, where it used to be
+  dropped. A total can only be formed once every line sits on the same
+  coordinate, and on a length axis the lines do not: each species, and the
+  resource, converts weight to length with its own allometric relationship, so
+  at a given length they sit at different weights. The total is therefore
+  summed *after* the conversion — at equal length rather than at equal weight —
+  interpolating each series onto the union of all the size coordinates,
+  logarithmically in size, with a series contributing nothing outside its own
+  range. Where the series already share a grid, which is always the case on a
+  weight axis, the union is that grid and nothing is approximated: the
+  weight-axis total is unchanged.
+
+  Two consequences on the weight axis. `plotSpectra(total = TRUE,
+  resource = FALSE)` now leaves the resource out of the total, where the
+  `resource` argument used to control only whether the resource was drawn. And
+  the total of a mizer array is now the sum of the lines that are drawn, so it
+  no longer counts abundance outside a species' size range; pass
+  `all.sizes = TRUE` for the old sum over the whole array.
+
 - `plotYieldObservedVsModel()` gains a `gear` argument that restricts the
   comparison to the catch of the selected gears. Both the model yield and the
   observed yield are then taken from those gears only, so in a model where
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index 65c3f9065..d1769771d 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -364,6 +364,7 @@ plot.ArraySpeciesBySize <- function(x, species = NULL,
     plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
                                           density_wrt = array_density_wrt(x),
                                           per_log_size = per_log_size)
+    if (total) plot_dat <- add_total_line(plot_dat)
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -557,7 +558,8 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
                             y_ticks = y_ticks, legend_var = "Legend",
                             size_axis = size_axis,
                             density_wrt = array_density_wrt(x),
-                            per_log_size = per_log_size)
+                            per_log_size = per_log_size,
+                            total = total)
 }
 
 #' Plot relative difference between two mizer arrays
@@ -675,7 +677,8 @@ plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
                           xtrans = if (log_x) "log10" else "identity",
                           xlim = plot_size_xlim(wlim, size_axis, llim),
                           ylim = ylim,
-                          legend_var = "Legend", size_axis = size_axis)
+                          legend_var = "Legend", size_axis = size_axis,
+                          total = total)
 }
 
 check_plot2_compatible <- function(x, y, class) {
@@ -973,6 +976,7 @@ addPlot.ArraySpeciesBySize <- function(plot, x, species = NULL,
     plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
                                           density_wrt = array_density_wrt(x),
                                           per_log_size = per_log_size)
+    if (total) plot_dat <- add_total_line(plot_dat)
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -1127,11 +1131,6 @@ prepare_ArraySpeciesBySize_plot_data <- function(x, species = NULL,
     sel <- all_species %in% species
     mat <- unclass(x)[sel, , drop = FALSE]
 
-    # Compute total across all selected species before size-range trimming
-    if (total) {
-        total_row <- colSums(mat, na.rm = TRUE)
-    }
-
     plot_dat <- data.frame(
         w = rep(w, each = sum(sel)),
         value = c(mat),
@@ -1165,20 +1164,14 @@ prepare_ArraySpeciesBySize_plot_data <- function(x, species = NULL,
         }
     }
 
-    # Add total line
-    if (total) {
-        total_dat <- data.frame(
-            w = w,
-            value = total_row,
-            Species = "Total",
-            Legend = "Total"
-        )
-        total_dat <- apply_wlim(total_dat, wlim)
-        plot_dat <- rbind(plot_dat, total_dat)
-    }
-
     names(plot_dat)[2] <- value_name
 
+    # The total is deliberately not added here. It has to be formed after the
+    # size coordinate has been converted, because on a length axis the series
+    # no longer share a grid; see `add_total_line()`. The `total` argument is
+    # kept so that the callers can see what was asked for.
+    attr(plot_dat, "total") <- isTRUE(total)
+
     plot_dat
 }
 
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index ade6d64b4..4a5f1b989 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -434,13 +434,11 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                       stringsAsFactors = FALSE)
     df$value <- c(sub)
 
-    # Compute total across ALL selected species (including background) before
-    # any background filtering, matching the behaviour of plot.ArraySpeciesBySize
-    if (total) {
-        total_sums <- stats::aggregate(value ~ time + w, data = df, FUN = sum,
-                                       na.rm = TRUE)
-        total_sums$Species <- "Total"
-    }
+    # The contributors to the total are taken across ALL selected species
+    # (including background), before any background filtering, matching the
+    # behaviour of plot.ArraySpeciesBySize. They are summed only after the size
+    # axis has been converted, inside `animate_plotly()`.
+    total_dat <- if (total) df else NULL
 
     # Now handle background: exclude rows or group under "Background" legend
     df$legend_name <- df$Species
@@ -452,10 +450,7 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
         }
     }
 
-    if (total) {
-        total_sums$legend_name <- "Total"
-        df <- rbind(df, total_sums[, names(df)])
-    }
+    if (!is.null(total_dat)) total_dat$legend_name <- "Total"
 
     y_label <- array_y_label(x, default = "Value", size_axis = size_axis,
                              per_log_size = per_log_size)
@@ -465,6 +460,7 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                    size_axis = size_axis,
                    density_wrt = array_density_wrt(x),
                    per_log_size = per_log_size,
+                   total_dat = total_dat,
                    frame_duration = frame_duration,
                    transition_duration = transition_duration,
                    easing = easing)
diff --git a/R/animateSpectra.R b/R/animateSpectra.R
index 8b0517f66..9898240c1 100644
--- a/R/animateSpectra.R
+++ b/R/animateSpectra.R
@@ -197,17 +197,20 @@ animate.MizerSim <- function(x, species = NULL,
         nf <- rbind(nf, nf_back)
     }
     # Add total ----
+    # The contributors are assembled here but summed only after the size axis
+    # has been converted, inside `animate_plotly()`. Every species counts,
+    # whether or not it was selected for display, and the resource counts when
+    # it is shown, matching `plotSpectra()`.
+    total_dat <- NULL
     if (total) {
-        # Calculate total community abundance
-        fish_idx <- (length(sim@params@w_full) -
-                         length(sim@params@w) + 1):length(sim@params@w_full)
-        total_n <- sim@n_pp
-        total_n[, fish_idx] <- total_n[, fish_idx] +
-            rowSums(aperm(sim@n, c(1, 3, 2)), dims = 2)
-        nf_total <- melt(total_n[time_elements, , drop = FALSE])
-        nf_total$Species <- "Total"
-        nf_total$legend_name <- "Total"
-        nf <- rbind(nf, nf_total)
+        total_dat <- melt(sim@n[time_elements, , , drop = FALSE])
+        names(total_dat)[names(total_dat) == "sp"] <- "Species"
+        if (resource) {
+            nf_pp_total <- melt(sim@n_pp[time_elements, , drop = FALSE])
+            nf_pp_total$Species <- "Resource"
+            total_dat <- rbind(total_dat, nf_pp_total)
+        }
+        total_dat$legend_name <- "Total"
     }
 
     y_label <- spectra_y_label(power, size_axis,
@@ -220,13 +223,18 @@ animate.MizerSim <- function(x, species = NULL,
     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)
+        if (!is.null(total_dat)) total_dat$w <- total_dat$w * sqrt(beta)
     }
     nf <- mutate(nf, value = value * w^power)
+    if (!is.null(total_dat)) {
+        total_dat <- mutate(total_dat, value = value * w^power)
+    }
 
     animate_plotly(nf, sim@params, log_x, log_y, y_label, wlim, llim,
                    ylim,
                    size_axis = size_axis,
                    density_wrt = spectrum_density_wrt(spectrum$per_log_size),
+                   total_dat = total_dat,
                    frame_duration = frame_duration,
                    transition_duration = transition_duration,
                    easing = easing)
@@ -243,14 +251,27 @@ animate_plotly <- function(df, params, log_x, log_y, y_label,
                            size_axis = "w",
                            density_wrt = NA_character_,
                            per_log_size = NULL,
+                           total_dat = NULL,
                            frame_duration = 500, transition_duration = 500,
                            easing = "linear") {
     size_axis <- plot_size_axis(size_axis)
-    df <- convert_plot_density_axis(df, params, size_axis,
-                                    density_wrt = density_wrt,
-                                    per_log_size = per_log_size,
-                                    value_col = "value")
+    convert <- function(d) {
+        convert_plot_density_axis(d, params, size_axis,
+                                  density_wrt = density_wrt,
+                                  per_log_size = per_log_size,
+                                  value_col = "value")
+    }
+    df <- convert(df)
     x_var <- plot_size_x_var(size_axis)
+    # The total is summed over the converted series, so that on a length axis
+    # it is a sum at equal length rather than at equal weight.
+    if (!is.null(total_dat)) {
+        total_dat <- add_total_line(convert(total_dat), x_var, "value",
+                                    by = "time")
+        total_dat <- total_dat[total_dat$Species == "Total", ]
+        total_dat$legend_name <- "Total"
+        df <- rbind(df, total_dat[, names(df), drop = FALSE])
+    }
     legend_name_order <- intersect(names(params@linecolour),
                                    unique(df$legend_name))
     sp_order <- unlist(lapply(legend_name_order, function(ln) {
diff --git a/R/plots.R b/R/plots.R
index 43ad512e7..7c77bbd13 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -408,6 +408,8 @@ plotHover.mizer_plot <- function(x = ggplot2::last_plot(), ...,
 #'   and are left alone when the size axis changes.
 #' @param per_log_size Whether to express a density per logarithmic size.
 #'   `NULL` (the default) keeps whichever the values already are.
+#' @param total Whether to add a line for the total over the series, formed
+#'   after the size axis has been converted, see [add_total_line()].
 #' @return A `mizer_plot` (ggplot2) object.
 #' @keywords internal
 plotComparisonDataFrame <- function(frame1, frame2, params,
@@ -419,7 +421,8 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                     legend_var = "Legend",
                                     size_axis = NULL,
                                     density_wrt = NA_character_,
-                                    per_log_size = NULL) {
+                                    per_log_size = NULL,
+                                    total = FALSE) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -448,6 +451,10 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                            value_col = y_var)
         x_var <- plot_size_x_var(size_axis)
     }
+    if (total) {
+        frame <- add_total_line(frame, x_var, y_var, by = "Model")
+        frame$Model <- factor(frame$Model, levels = c(name1, name2))
+    }
 
     legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
     frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
@@ -508,7 +515,8 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
                                   ylim = c(NA, NA),
                                   highlight = NULL,
                                   legend_var = "Legend",
-                                  size_axis = NULL) {
+                                  size_axis = NULL,
+                                  total = FALSE) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -524,18 +532,28 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
              "in the data frame.")
     }
 
+    # The size axis is converted before the two frames are joined, so that a
+    # total can be formed on the axis it is plotted against. The relative
+    # difference itself is unaffected: a common Jacobian cancels.
+    if (!is.null(size_axis)) {
+        size_axis <- plot_size_axis(size_axis)
+        frame1 <- convert_plot_size_axis(frame1, params, size_axis,
+                                         species_col = group_var)
+        frame2 <- convert_plot_size_axis(frame2, params, size_axis,
+                                         species_col = group_var)
+        x_var <- plot_size_x_var(size_axis)
+    }
+    if (total) {
+        frame1 <- add_total_line(frame1, x_var, y_var)
+        frame2 <- add_total_line(frame2, x_var, y_var)
+    }
+
     by_vars <- c(x_var, group_var, legend_var)
     frame <- dplyr::inner_join(frame1, frame2, by = by_vars,
                                suffix = c(".x", ".y"))
     frame$rel_diff <- relative_difference(frame[[paste0(y_var, ".x")]],
                                           frame[[paste0(y_var, ".y")]])
     frame <- frame[is.finite(frame$rel_diff), ]
-    if (!is.null(size_axis)) {
-        size_axis <- plot_size_axis(size_axis)
-        frame <- convert_plot_size_axis(frame, params, size_axis,
-                                        species_col = group_var)
-        x_var <- plot_size_x_var(size_axis)
-    }
 
     legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
     frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
@@ -726,6 +744,72 @@ array_log_y <- function(x, log_y, log, given) {
     log_y
 }
 
+#' Add a total line to plotting data by summing over its series
+#'
+#' The total has to be formed *after* the size coordinate has been converted,
+#' not before. On a weight axis every series shares the model's weight grid, so
+#' summing at equal weight and summing at equal position are the same thing. On
+#' a length axis they are not: each species, and the resource, converts weight
+#' to length with its own allometric relationship, so at a given length the
+#' series sit at different weights and their grids no longer coincide. The sum
+#' that means something there is the sum at equal *length* — the number of
+#' organisms per unit length, whatever they are — which is what this computes.
+#'
+#' Each series is interpolated onto the sorted union of all the size
+#' coordinates, linearly in the logarithm of size, since the grid is
+#' logarithmic. A series contributes nothing outside its own range. When the
+#' series already share a grid — always on a weight axis, and on a length axis
+#' whenever the weight-length parameters agree — the union is that grid and the
+#' interpolation reproduces the values exactly, so nothing is approximated in
+#' the cases where nothing needs to be.
+#'
+#' @param plot_dat A data frame of plotting data with a size column, a value
+#'   column and a `Species` column.
+#' @param x_var Name of the size column. Defaults to the first column.
+#' @param value_col Name or index of the value column. Defaults to the second.
+#' @param by Names of further columns identifying separate plots, such as the
+#'   time of an animation frame or the model of a comparison. A total is formed
+#'   within each of their combinations.
+#' @return `plot_dat` with the total appended as a series named `"Total"`.
+#' @keywords internal
+add_total_line <- function(plot_dat, x_var = names(plot_dat)[[1]],
+                           value_col = 2, by = NULL) {
+    if (nrow(plot_dat) == 0) return(plot_dat)
+    if (is.numeric(value_col)) value_col <- names(plot_dat)[[value_col]]
+    if (length(by) > 0) {
+        groups <- split(plot_dat, plot_dat[, by, drop = FALSE], drop = TRUE)
+        out <- lapply(groups, add_total_line, x_var = x_var,
+                      value_col = value_col)
+        out <- do.call(rbind, unname(out))
+        rownames(out) <- NULL
+        return(out)
+    }
+    x_all <- sort(unique(plot_dat[[x_var]]))
+    total <- numeric(length(x_all))
+    for (series in split(plot_dat, plot_dat$Species)) {
+        series <- series[order(series[[x_var]]), , drop = FALSE]
+        x <- series[[x_var]]
+        y <- series[[value_col]]
+        if (length(x) < 2) {
+            # `approx()` needs two points; a single one contributes only where
+            # it sits.
+            hit <- match(x, x_all)
+            total[hit] <- total[hit] + y
+            next
+        }
+        contribution <- stats::approx(log(x), y, xout = log(x_all))$y
+        contribution[is.na(contribution)] <- 0
+        total <- total + contribution
+    }
+    total_dat <- plot_dat[rep(1, length(x_all)), , drop = FALSE]
+    total_dat[[x_var]] <- x_all
+    total_dat[[value_col]] <- total
+    total_dat$Species <- "Total"
+    if ("Legend" %in% names(total_dat)) total_dat$Legend <- "Total"
+    rownames(total_dat) <- NULL
+    rbind(plot_dat, total_dat)
+}
+
 #' The weight-length parameters to plot each row of plotting data with
 #'
 #' A length axis needs an allometric weight-length relationship for every line
@@ -1833,13 +1917,6 @@ plot_spectra <- function(params, n, n_pp,
         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 * w_full_grid^power
-    }
     species <- valid_species_arg(params, species)
     y_label <- spectra_y_label(power, size_axis, biomass = biomass,
                                per_log_size = per_log_size)
@@ -1867,14 +1944,6 @@ plot_spectra <- function(params, n, n_pp,
             )
         }
     }
-    if (total) {
-        plot_dat <- rbind(plot_dat,
-                          data.frame(w = w_full_grid,
-                                     value = c(total_n),
-                                     Species = "Total",
-                                     Legend = "Total")
-                          )
-    }
     if (background && any(params@species_params$is_background)) {
         back_n <- n[params@species_params$is_background, , drop = FALSE]
         plot_dat <-
@@ -1892,6 +1961,37 @@ plot_spectra <- function(params, n, n_pp,
                              (plot_dat$w <= wlim[2]), ]
     plot_dat <- convert_plot_spectrum_axis(plot_dat, params, size_axis, power,
                                            per_log_size = per_log_size)
+    if (total) {
+        # The total is summed over the series as they are plotted, which on a
+        # length axis means at equal length rather than at equal weight. It is
+        # built from every species, whether or not it was selected for display,
+        # and from the resource when the resource is shown.
+        # `n` has already been multiplied by the power of weight above.
+        total_dat <- data.frame(w = rep(w_grid, each = nrow(n)),
+                                value = c(unclass(n)),
+                                Species = rep(dimnames(n)[[1]],
+                                              times = length(w_grid)),
+                                Legend = "Total")
+        if (resource) {
+            resource_dat <- data.frame(w = w_full_grid,
+                                       value = c(n_pp * w_full_grid^power),
+                                       Species = "Resource",
+                                       Legend = "Total")
+            total_dat <- rbind(total_dat, resource_dat)
+        }
+        # Zeros are kept here, unlike in the plotted data: they are part of
+        # each series and dropping them would let the interpolation bridge a
+        # gap that is really empty.
+        total_dat <- total_dat[total_dat$w >= wlim[1] &
+                                   total_dat$w <= wlim[2], ]
+        total_dat <- convert_plot_spectrum_axis(total_dat, params, size_axis,
+                                                power,
+                                                per_log_size = per_log_size)
+        total_dat <- add_total_line(total_dat)
+        total_dat <- total_dat[total_dat$Species == "Total" &
+                                   total_dat$value > 0, ]
+        plot_dat <- rbind(plot_dat, total_dat[, names(plot_dat), drop = FALSE])
+    }
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
diff --git a/inst/skills/upgrade-mizer-code/SKILL.md b/inst/skills/upgrade-mizer-code/SKILL.md
index f45f82695..21e1dacff 100644
--- a/inst/skills/upgrade-mizer-code/SKILL.md
+++ b/inst/skills/upgrade-mizer-code/SKILL.md
@@ -52,6 +52,9 @@ plots) are in the changelog and are not repeated here.
 | A custom array plotted with `size_axis = "l"` is not transformed, or is transformed when it should not be | an array declares what it holds with `type` now | Arrays say what kind of value they hold (3.3) |
 | `array_spectrum_power()`, or `spectrum_power =` in an internal plot helper, is no longer found | replaced by the `type` metadata | Arrays say what kind of value they hold (3.3) |
 | New warning that two arrays hold "a value of type" different things in `plot2()` or `plotRelative()` | the two arrays disagree about what they hold | Arrays say what kind of value they hold (3.3) |
+| A `Total` line appears on a plot with `size_axis = "l"` where there used to be none | the total is now summed after the length conversion | The total is summed on the axis it is plotted against (3.3) |
+| `plotSpectra(total = TRUE, resource = FALSE)` gives a smaller total than before | the `resource` argument now controls the total as well as the line | The total is summed on the axis it is plotted against (3.3) |
+| `plot(<array>, total = TRUE)` gives a smaller total at some sizes | the total is now the sum of the lines drawn, not of the whole array | The total is summed on the axis it is plotted against (3.3) |
 | New warning that a change to a species or resource parameter "has not taken effect" | the rate it feeds was set by hand and is no longer calculated | A change that cannot take effect now warns (3.3) |
 | That warning appears with `given_species_params<-()` but not with `species_params<-()` | the diagnostics belong to the given species parameter setter | The two species parameter setters divide the diagnostics between them (3.3) |
 | Setting a given species parameter to `NA` now warns that the change has not taken effect | clearing a value counts as a change, and a frozen array blocks it | The two species parameter setters divide the diagnostics between them (3.3) |
@@ -239,6 +242,35 @@ single `density_wrt` argument in place of `spectrum_power` and
 `spectrum_per_log_size`, and that the internal `array_spectrum_power()` is gone.
 The `power`-based interface of `plotSpectra()` and friends is unchanged.
 
+### The total is summed on the axis it is plotted against
+
+A total can only be formed once every line sits on the same coordinate. On a
+weight axis they do: every species shares the model's weight grid. On a length
+axis they do not, because each species — and now the resource — converts weight
+to length with its own allometric relationship, so at a given length the lines
+sit at different weights. That is why the `Total` line used to be dropped from
+length-based plots.
+
+It is now summed *after* the conversion, at equal length rather than at equal
+weight, interpolating each line onto the union of all the size coordinates
+(logarithmically in size, with a line contributing nothing outside its own
+range). Where the lines already share a grid — always on a weight axis, and on
+a length axis whenever the weight-length parameters agree — the union is that
+grid and the interpolation reproduces the values exactly. **The weight-axis
+total is unchanged**, for every power.
+
+Two things do change on the weight axis, both of them corrections:
+
+- `plotSpectra(total = TRUE, resource = FALSE)` now leaves the resource out of
+  the total. The `resource` argument used to control only whether the resource
+  line was drawn; the total counted it either way, so a "consumers only" total
+  was not available.
+- `plot(<array>, total = TRUE)` now sums the lines that are drawn rather than
+  the whole array, so it no longer counts abundance outside a species' size
+  range. On the North Sea model this moves the total at about 2% of sizes, by up
+  to a factor of three where a species holds abundance beyond its `w_max`. Pass
+  `all.sizes = TRUE` to sum the whole array as before.
+
 ### Length and weight parameters follow the one you gave last
 
 A size can be given either as a weight (`w_mat`, `w_max`, …) or as the length it
diff --git a/man/add_total_line.Rd b/man/add_total_line.Rd
new file mode 100644
index 000000000..b257a4026
--- /dev/null
+++ b/man/add_total_line.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plots.R
+\name{add_total_line}
+\alias{add_total_line}
+\title{Add a total line to plotting data by summing over its series}
+\usage{
+add_total_line(
+  plot_dat,
+  x_var = names(plot_dat)[[1]],
+  value_col = 2,
+  by = NULL
+)
+}
+\arguments{
+\item{plot_dat}{A data frame of plotting data with a size column, a value
+column and a \code{Species} column.}
+
+\item{x_var}{Name of the size column. Defaults to the first column.}
+
+\item{value_col}{Name or index of the value column. Defaults to the second.}
+
+\item{by}{Names of further columns identifying separate plots, such as the
+time of an animation frame or the model of a comparison. A total is formed
+within each of their combinations.}
+}
+\value{
+\code{plot_dat} with the total appended as a series named \code{"Total"}.
+}
+\description{
+The total has to be formed \emph{after} the size coordinate has been converted,
+not before. On a weight axis every series shares the model's weight grid, so
+summing at equal weight and summing at equal position are the same thing. On
+a length axis they are not: each species, and the resource, converts weight
+to length with its own allometric relationship, so at a given length the
+series sit at different weights and their grids no longer coincide. The sum
+that means something there is the sum at equal \emph{length} — the number of
+organisms per unit length, whatever they are — which is what this computes.
+}
+\details{
+Each series is interpolated onto the sorted union of all the size
+coordinates, linearly in the logarithm of size, since the grid is
+logarithmic. A series contributes nothing outside its own range. When the
+series already share a grid — always on a weight axis, and on a length axis
+whenever the weight-length parameters agree — the union is that grid and the
+interpolation reproduces the values exactly, so nothing is approximated in
+the cases where nothing needs to be.
+}
+\keyword{internal}
diff --git a/man/plotComparisonDataFrame.Rd b/man/plotComparisonDataFrame.Rd
index aecaacaa8..eb1b3534f 100644
--- a/man/plotComparisonDataFrame.Rd
+++ b/man/plotComparisonDataFrame.Rd
@@ -21,7 +21,8 @@ plotComparisonDataFrame(
   legend_var = "Legend",
   size_axis = NULL,
   density_wrt = NA_character_,
-  per_log_size = NULL
+  per_log_size = NULL,
+  total = FALSE
 )
 }
 \arguments{
@@ -57,6 +58,9 @@ and are left alone when the size axis changes.}
 
 \item{per_log_size}{Whether to express a density per logarithmic size.
 \code{NULL} (the default) keeps whichever the values already are.}
+
+\item{total}{Whether to add a line for the total over the series, formed
+after the size axis has been converted, see \code{\link[=add_total_line]{add_total_line()}}.}
 }
 \value{
 A \code{mizer_plot} (ggplot2) object.
diff --git a/man/plotRelativeDataFrame.Rd b/man/plotRelativeDataFrame.Rd
index 1f7ff1414..5a596277d 100644
--- a/man/plotRelativeDataFrame.Rd
+++ b/man/plotRelativeDataFrame.Rd
@@ -14,7 +14,8 @@ plotRelativeDataFrame(
   ylim = c(NA, NA),
   highlight = NULL,
   legend_var = "Legend",
-  size_axis = NULL
+  size_axis = NULL,
+  total = FALSE
 )
 }
 \arguments{
diff --git a/tests/testthat/test-ArraySpeciesBySize-class.R b/tests/testthat/test-ArraySpeciesBySize-class.R
index 94e643e9d..76d9d8f2f 100644
--- a/tests/testthat/test-ArraySpeciesBySize-class.R
+++ b/tests/testthat/test-ArraySpeciesBySize-class.R
@@ -276,6 +276,55 @@ test_that("plot() can express a density per logarithmic size", {
     expect_true(any(with_total$Species == "Total"))
 })
 
+test_that("the total is shown on a length axis and matches the lines drawn", {
+    density <- initialN(NS_params_small)
+
+    on_w <- plot(density, total = TRUE, return_data = TRUE)
+    on_l <- plot(density, total = TRUE, size_axis = "l", return_data = TRUE)
+    expect_true("Total" %in% on_w$Species)
+    expect_true("Total" %in% on_l$Species)
+
+    # The species here share one weight-length relationship, so the total is
+    # the plain sum of the converted lines at each length
+    total_l <- on_l[on_l$Species == "Total", ]
+    lines_l <- on_l[on_l$Species != "Total", ]
+    summed <- vapply(total_l$l, function(len) {
+        sum(lines_l[[2]][lines_l$l == len])
+    }, numeric(1))
+    expect_equal(total_l[[2]], summed)
+
+    # The total is the sum of what is drawn, so trimming to each species' size
+    # range is respected; `all.sizes = TRUE` sums the whole array
+    all_sizes <- plot(density, total = TRUE, all.sizes = TRUE,
+                      return_data = TRUE)
+    expect_equal(all_sizes[[2]][all_sizes$Species == "Total"],
+                 unname(colSums(unclass(density))))
+})
+
+test_that("the total on a length axis interpolates onto a common grid", {
+    # Give two species their own weight-length relationships so that their
+    # length grids no longer coincide
+    params <- NS_params_small
+    species_params(params)$a <- c(0.008, 0.012, 0.01)
+    species_params(params)$b <- c(3.1, 2.9, 3)
+    density <- initialN(params)
+
+    on_w <- plot(density, total = TRUE, return_data = TRUE)
+    on_l <- plot(density, total = TRUE, size_axis = "l", return_data = TRUE)
+    # The weight axis still has one point per size, the length axis needs the
+    # union of the three grids
+    expect_identical(sum(on_w$Species == "Total"), length(params@w))
+    expect_gt(sum(on_l$Species == "Total"), length(params@w))
+
+    # At the largest length only the largest species is present, so the total
+    # is exactly that species there
+    total_l <- on_l[on_l$Species == "Total", ]
+    lines_l <- on_l[on_l$Species != "Total", ]
+    at <- max(total_l$l)
+    expect_identical(sum(lines_l$l == at), 1L)
+    expect_equal(total_l[[2]][total_l$l == at], lines_l[[2]][lines_l$l == at])
+})
+
 test_that("a proportion is plotted against the whole of [0, 1]", {
     feeding_level <- getFeedingLevel(NS_params_small)
     expect_identical(array_type(feeding_level), "proportion")
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 00d8e2c82..dc8d82053 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -508,14 +508,14 @@ test_that("size-based plots support length axes", {
     expect_true(all(spectra_l_limited$l >= llim[1]))
     expect_true(all(spectra_l_limited$l <= llim[2]))
 
-    # The resource has its own weight-length relationship and so appears on a
-    # length axis; the total still does not, because it is a sum over lines
-    # that each have their own.
+    # The resource and the total both appear on a length axis: the resource
+    # has its own weight-length relationship, and the total is summed after
+    # the conversion, at equal length.
     spectra_all <- plotSpectra(params_len, species = species,
                                resource = TRUE, total = TRUE,
                                size_axis = "l", return_data = TRUE)
     expect_true("Resource" %in% spectra_all$Legend)
-    expect_false("Total" %in% spectra_all$Legend)
+    expect_true("Total" %in% spectra_all$Legend)
 
     p <- plotSpectra(params_len, species = species, resource = FALSE,
                      size_axis = "l")
@@ -642,6 +642,94 @@ test_that("proportion_ylim widens to [0, 1] without ever hiding data", {
                      c(0, 1))
 })
 
+test_that("add_total_line sums series on a shared grid exactly", {
+    plot_dat <- data.frame(w = rep(c(1, 2, 4), 2), value = 1:6,
+                           Species = rep(c("a", "b"), each = 3))
+    out <- add_total_line(plot_dat)
+    total <- out[out$Species == "Total", ]
+    expect_identical(total$w, c(1, 2, 4))
+    expect_equal(total$value, c(1 + 4, 2 + 5, 3 + 6))
+    # The original rows are untouched
+    expect_equal(out[out$Species != "Total", names(plot_dat)], plot_dat,
+                 ignore_attr = TRUE)
+})
+
+test_that("add_total_line interpolates series that do not share a grid", {
+    # `b` lives on a shifted grid, so the total is evaluated on the union and
+    # each series is interpolated onto it, contributing nothing outside its own
+    # range.
+    plot_dat <- rbind(
+        data.frame(w = c(1, 4), value = c(1, 4), Species = "a"),
+        data.frame(w = c(2, 8), value = c(10, 40), Species = "b"))
+    total <- add_total_line(plot_dat)
+    total <- total[total$Species == "Total", ]
+    expect_identical(total$w, c(1, 2, 4, 8))
+    # At w = 1 only `a`; at w = 8 only `b`
+    expect_equal(total$value[[1]], 1)
+    expect_equal(total$value[[4]], 40)
+    # At w = 2, `a` interpolated linearly in log w between (1, 1) and (4, 4)
+    a_at_2 <- stats::approx(log(c(1, 4)), c(1, 4), xout = log(2))$y
+    expect_equal(total$value[[2]], a_at_2 + 10)
+    # At w = 4, `b` interpolated between (2, 10) and (8, 40)
+    b_at_4 <- stats::approx(log(c(2, 8)), c(10, 40), xout = log(4))$y
+    expect_equal(total$value[[3]], 4 + b_at_4)
+})
+
+test_that("add_total_line groups by the columns it is told to", {
+    plot_dat <- data.frame(w = rep(c(1, 2), 4), value = 1:8,
+                           Species = rep(c("a", "b"), each = 2),
+                           time = rep(c(1, 1, 2, 2), each = 2))
+    out <- add_total_line(plot_dat, by = "time")
+    totals <- out[out$Species == "Total", ]
+    expect_identical(nrow(totals), 4L)
+    # time 1 holds a = (1, 2) and b = (3, 4); time 2 holds a = (5, 6) and
+    # b = (7, 8)
+    expect_equal(totals$value[totals$time == 1], c(1 + 3, 2 + 4))
+    expect_equal(totals$value[totals$time == 2], c(5 + 7, 6 + 8))
+})
+
+test_that("the spectrum total is unchanged on a weight axis", {
+    # The total is now summed after the size axis is converted rather than
+    # before, which on a weight axis must make no difference at all.
+    for (power in 0:2) {
+        d <- plotSpectra(params, total = TRUE, power = power,
+                         return_data = TRUE)
+        total <- d[d$Legend == "Total", ]
+        fish_idx <- (length(params@w_full) - length(params@w) +
+                         1):length(params@w_full)
+        expected <- params@initial_n_pp
+        expected[fish_idx] <- expected[fish_idx] + colSums(params@initial_n)
+        expected <- expected * params@w_full^power
+        keep <- expected > 0 & params@w_full >= min(params@w) / 100
+        expect_equal(total[[2]], unname(expected[keep]))
+    }
+})
+
+test_that("the spectrum total appears on a length axis", {
+    on_l <- plotSpectra(params, total = TRUE, size_axis = "l",
+                        return_data = TRUE)
+    total <- on_l[on_l$Legend == "Total", ]
+    expect_gt(nrow(total), 0)
+
+    # All the species here share one weight-length relationship, so at each
+    # length the total is exactly the sum of the series drawn there
+    contributors <- on_l[on_l$Legend != "Total", ]
+    at <- total$l[[which.max(total$l)]]
+    expect_equal(total[[2]][total$l == at],
+                 sum(contributors[[2]][contributors$l == at]))
+
+    # The resource is in the total when it is shown and not when it is not,
+    # where it used to be counted either way
+    with_resource <- plotSpectra(params, total = TRUE, resource = TRUE,
+                                 return_data = TRUE)
+    without <- plotSpectra(params, total = TRUE, resource = FALSE,
+                           return_data = TRUE)
+    small <- min(params@w)
+    expect_gt(with_resource[[2]][with_resource$Legend == "Total" &
+                                     with_resource$w == small],
+              without[[2]][without$Legend == "Total" & without$w == small])
+})
+
 test_that("validate_density_wrt accepts only the known measures", {
     expect_identical(validate_density_wrt(NULL), NA_character_)
     expect_identical(validate_density_wrt(NA), NA_character_)
diff --git a/vignettes/upgrading.Rmd b/vignettes/upgrading.Rmd
index 85d24fa92..5fe105570 100644
--- a/vignettes/upgrading.Rmd
+++ b/vignettes/upgrading.Rmd
@@ -123,6 +123,35 @@ single `density_wrt` argument in place of `spectrum_power` and
 `spectrum_per_log_size`, and that the internal `array_spectrum_power()` is gone.
 The `power`-based interface of `plotSpectra()` and friends is unchanged.
 
+### The total is summed on the axis it is plotted against
+
+A total can only be formed once every line sits on the same coordinate. On a
+weight axis they do: every species shares the model's weight grid. On a length
+axis they do not, because each species — and now the resource — converts weight
+to length with its own allometric relationship, so at a given length the lines
+sit at different weights. That is why the `Total` line used to be dropped from
+length-based plots.
+
+It is now summed *after* the conversion, at equal length rather than at equal
+weight, interpolating each line onto the union of all the size coordinates
+(logarithmically in size, with a line contributing nothing outside its own
+range). Where the lines already share a grid — always on a weight axis, and on
+a length axis whenever the weight-length parameters agree — the union is that
+grid and the interpolation reproduces the values exactly. **The weight-axis
+total is unchanged**, for every power.
+
+Two things do change on the weight axis, both of them corrections:
+
+- `plotSpectra(total = TRUE, resource = FALSE)` now leaves the resource out of
+  the total. The `resource` argument used to control only whether the resource
+  line was drawn; the total counted it either way, so a "consumers only" total
+  was not available.
+- `plot(<array>, total = TRUE)` now sums the lines that are drawn rather than
+  the whole array, so it no longer counts abundance outside a species' size
+  range. On the North Sea model this moves the total at about 2% of sizes, by up
+  to a factor of three where a species holds abundance beyond its `w_max`. Pass
+  `all.sizes = TRUE` to sum the whole array as before.
+
 ### Length and weight parameters follow the one you gave last
 
 A size can be given either as a weight (`w_mat`, `w_max`, …) or as the length it

From 018dd99240e583dd0bbcacc46bc6472f0a2ce966 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Sun, 16 Aug 2026 07:49:32 +0100
Subject: [PATCH 2/4] Make `total` mean the total of everything the object
 holds

`total = TRUE` meant two different things. `plotSpectra()` summed the
resource and every species, whatever was drawn and whichever species were
selected. The array plots summed only the species selected for display,
and only the sizes inside each species' own range.

The first is the useful meaning: the total is a property of the model,
not of the plot, so that a plot of two species can be read against the
community total. Both now work that way.

- `plotSpectra()` and `animate()` on a `MizerSim` keep counting the
  resource whether or not it is drawn. The previous commit had made the
  `resource` argument exclude it from the total as well as from the plot;
  that is reverted.
- `plot()`, `addPlot()`, `plot2()`, `plotRelative()` and `animate()` on
  an array now sum the whole array. The total no longer moves when
  `species`, `all.sizes` or `background` change.

The contributors are assembled by `total_contributors()`, which is just
the usual preparation of the whole array, and summed after the size
conversion by `append_total_line()`.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 NEWS.md                                       | 12 +--
 R/ArraySpeciesBySize-class.R                  | 47 +++++++--
 R/ArrayTimeBySpeciesBySize-class.R            | 15 ++-
 R/animateSpectra.R                            | 15 ++-
 R/plots.R                                     | 95 ++++++++++++++-----
 inst/skills/analyse-and-plot/SKILL.md         |  2 +-
 inst/skills/upgrade-mizer-code/SKILL.md       | 24 +++--
 man/append_total_line.Rd                      | 38 ++++++++
 man/plotComparisonDataFrame.Rd                |  8 +-
 man/plotRelativeDataFrame.Rd                  |  3 +-
 man/total_contributors.Rd                     | 29 ++++++
 .../testthat/test-ArraySpeciesBySize-class.R  | 25 +++--
 tests/testthat/test-plots.R                   | 18 ++--
 .../cheatsheet-analysis-and-plotting.Rmd      |  2 +-
 vignettes/upgrading.Rmd                       | 21 ++--
 15 files changed, 252 insertions(+), 102 deletions(-)
 create mode 100644 man/append_total_line.Rd
 create mode 100644 man/total_contributors.Rd

diff --git a/NEWS.md b/NEWS.md
index 574a925ff..ad31d30e3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -171,12 +171,12 @@ stability of steady states.
   weight axis, the union is that grid and nothing is approximated: the
   weight-axis total is unchanged.
 
-  Two consequences on the weight axis. `plotSpectra(total = TRUE,
-  resource = FALSE)` now leaves the resource out of the total, where the
-  `resource` argument used to control only whether the resource was drawn. And
-  the total of a mizer array is now the sum of the lines that are drawn, so it
-  no longer counts abundance outside a species' size range; pass
-  `all.sizes = TRUE` for the old sum over the whole array.
+  `total = TRUE` also now means the same thing everywhere: the total of
+  everything the object holds. For `plotSpectra()` that was already so — the
+  resource and every species, whatever is drawn — and it stays so. The array
+  plots have been brought into line: `plot(<array>, total = TRUE)` used to sum
+  only the species that were selected for display, and now sums the whole
+  array, so a plot of two species can be read against the community total.
 
 - `plotYieldObservedVsModel()` gains a `gear` argument that restricts the
   comparison to the catch of the selected gears. Both the model yield and the
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index d1769771d..8d3264c0f 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -364,7 +364,10 @@ plot.ArraySpeciesBySize <- function(x, species = NULL,
     plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
                                           density_wrt = array_density_wrt(x),
                                           per_log_size = per_log_size)
-    if (total) plot_dat <- add_total_line(plot_dat)
+    if (total) {
+        plot_dat <- append_total_line(plot_dat, total_contributors(x, wlim),
+                                      params, size_axis, x, per_log_size)
+    }
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -559,7 +562,12 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
                             size_axis = size_axis,
                             density_wrt = array_density_wrt(x),
                             per_log_size = per_log_size,
-                            total = total)
+                            total_dat = if (total) {
+                                rbind(cbind(total_contributors(x, wlim),
+                                            Model = name1),
+                                      cbind(total_contributors(y, wlim),
+                                            Model = name2))
+                            })
 }
 
 #' Plot relative difference between two mizer arrays
@@ -678,7 +686,8 @@ plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
                           xlim = plot_size_xlim(wlim, size_axis, llim),
                           ylim = ylim,
                           legend_var = "Legend", size_axis = size_axis,
-                          total = total)
+                          total_dat1 = if (total) total_contributors(x, wlim),
+                          total_dat2 = if (total) total_contributors(y, wlim))
 }
 
 check_plot2_compatible <- function(x, y, class) {
@@ -976,7 +985,10 @@ addPlot.ArraySpeciesBySize <- function(plot, x, species = NULL,
     plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
                                           density_wrt = array_density_wrt(x),
                                           per_log_size = per_log_size)
-    if (total) plot_dat <- add_total_line(plot_dat)
+    if (total) {
+        plot_dat <- append_total_line(plot_dat, total_contributors(x, wlim),
+                                      params, size_axis, x, per_log_size)
+    }
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -1166,15 +1178,30 @@ prepare_ArraySpeciesBySize_plot_data <- function(x, species = NULL,
 
     names(plot_dat)[2] <- value_name
 
-    # The total is deliberately not added here. It has to be formed after the
-    # size coordinate has been converted, because on a length axis the series
-    # no longer share a grid; see `add_total_line()`. The `total` argument is
-    # kept so that the callers can see what was asked for.
-    attr(plot_dat, "total") <- isTRUE(total)
-
     plot_dat
 }
 
+#' Assemble the contributors to the total of a species-by-size array
+#'
+#' The total is the total of everything the array holds: every species, whether
+#' or not it was selected for display, and every size, whether or not it falls
+#' in a species' own size range. It is a property of the array rather than of
+#' the plot, so that a plot of two species can still be read against the
+#' community total.
+#'
+#' The rows are returned unsummed, because the sum has to be taken after the
+#' size coordinate has been converted — on a length axis the species no longer
+#' share a grid; see [add_total_line()].
+#'
+#' @param x An `ArraySpeciesBySize` object.
+#' @param wlim Numeric vector of length two giving the weight limits.
+#' @return A data frame of plotting data holding every value in the array.
+#' @keywords internal
+total_contributors <- function(x, wlim = c(NA, NA)) {
+    prepare_ArraySpeciesBySize_plot_data(x, species = NULL, all.sizes = TRUE,
+                                         wlim = wlim, background = TRUE)
+}
+
 #' @rdname plotHover
 #' @usage NULL
 #' @examples
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 4a5f1b989..50b124780 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -434,11 +434,16 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                       stringsAsFactors = FALSE)
     df$value <- c(sub)
 
-    # The contributors to the total are taken across ALL selected species
-    # (including background), before any background filtering, matching the
-    # behaviour of plot.ArraySpeciesBySize. They are summed only after the size
-    # axis has been converted, inside `animate_plotly()`.
-    total_dat <- if (total) df else NULL
+    # The total is the total of everything the array holds, so its contributors
+    # are every species rather than the selected ones, matching
+    # `plot.ArraySpeciesBySize()`. They are summed only after the size axis has
+    # been converted, inside `animate_plotly()`.
+    total_dat <- NULL
+    if (total) {
+        total_dat <- expand.grid(time = times, Species = all_species, w = w,
+                                 stringsAsFactors = FALSE)
+        total_dat$value <- c(arr)
+    }
 
     # Now handle background: exclude rows or group under "Background" legend
     df$legend_name <- df$Species
diff --git a/R/animateSpectra.R b/R/animateSpectra.R
index 9898240c1..c45e1ceaf 100644
--- a/R/animateSpectra.R
+++ b/R/animateSpectra.R
@@ -198,18 +198,17 @@ animate.MizerSim <- function(x, species = NULL,
     }
     # Add total ----
     # The contributors are assembled here but summed only after the size axis
-    # has been converted, inside `animate_plotly()`. Every species counts,
-    # whether or not it was selected for display, and the resource counts when
-    # it is shown, matching `plotSpectra()`.
+    # has been converted, inside `animate_plotly()`. The total is the total of
+    # everything the model holds: every species, whether or not it was selected
+    # for display, and the resource, whether or not it is drawn — matching
+    # `plotSpectra()`.
     total_dat <- NULL
     if (total) {
         total_dat <- melt(sim@n[time_elements, , , drop = FALSE])
         names(total_dat)[names(total_dat) == "sp"] <- "Species"
-        if (resource) {
-            nf_pp_total <- melt(sim@n_pp[time_elements, , drop = FALSE])
-            nf_pp_total$Species <- "Resource"
-            total_dat <- rbind(total_dat, nf_pp_total)
-        }
+        nf_pp_total <- melt(sim@n_pp[time_elements, , drop = FALSE])
+        nf_pp_total$Species <- "Resource"
+        total_dat <- rbind(total_dat, nf_pp_total)
         total_dat$legend_name <- "Total"
     }
 
diff --git a/R/plots.R b/R/plots.R
index 7c77bbd13..3b1e6e850 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -408,8 +408,10 @@ plotHover.mizer_plot <- function(x = ggplot2::last_plot(), ...,
 #'   and are left alone when the size axis changes.
 #' @param per_log_size Whether to express a density per logarithmic size.
 #'   `NULL` (the default) keeps whichever the values already are.
-#' @param total Whether to add a line for the total over the series, formed
-#'   after the size axis has been converted, see [add_total_line()].
+#' @param total_dat Optional data frame of the contributors to a total, with a
+#'   `Model` column identifying which of the two they belong to. The total is
+#'   summed from them after the size axis has been converted, see
+#'   [add_total_line()].
 #' @return A `mizer_plot` (ggplot2) object.
 #' @keywords internal
 plotComparisonDataFrame <- function(frame1, frame2, params,
@@ -422,7 +424,7 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                     size_axis = NULL,
                                     density_wrt = NA_character_,
                                     per_log_size = NULL,
-                                    total = FALSE) {
+                                    total_dat = NULL) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -451,8 +453,17 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                            value_col = y_var)
         x_var <- plot_size_x_var(size_axis)
     }
-    if (total) {
-        frame <- add_total_line(frame, x_var, y_var, by = "Model")
+    if (!is.null(total_dat)) {
+        # The contributors arrive unconverted and are put on the same axis as
+        # the data they join before being summed there.
+        total_dat <- convert_plot_density_axis(total_dat, params, size_axis,
+                                               density_wrt = density_wrt,
+                                               per_log_size = per_log_size,
+                                               species_col = group_var,
+                                               value_col = y_var)
+        total_dat <- add_total_line(total_dat, x_var, y_var, by = "Model")
+        total_dat <- total_dat[total_dat[[group_var]] == "Total", ]
+        frame <- rbind(frame, total_dat[, names(frame), drop = FALSE])
         frame$Model <- factor(frame$Model, levels = c(name1, name2))
     }
 
@@ -516,7 +527,7 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
                                   highlight = NULL,
                                   legend_var = "Legend",
                                   size_axis = NULL,
-                                  total = FALSE) {
+                                  total_dat1 = NULL, total_dat2 = NULL) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -543,9 +554,19 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
                                          species_col = group_var)
         x_var <- plot_size_x_var(size_axis)
     }
-    if (total) {
-        frame1 <- add_total_line(frame1, x_var, y_var)
-        frame2 <- add_total_line(frame2, x_var, y_var)
+    if (!is.null(total_dat1) && !is.null(total_dat2)) {
+        add_total <- function(frame, total_dat) {
+            if (!is.null(size_axis)) {
+                total_dat <- convert_plot_size_axis(total_dat, params,
+                                                    size_axis,
+                                                    species_col = group_var)
+            }
+            total_dat <- add_total_line(total_dat, x_var, y_var)
+            total_dat <- total_dat[total_dat[[group_var]] == "Total", ]
+            rbind(frame, total_dat[, names(frame), drop = FALSE])
+        }
+        frame1 <- add_total(frame1, total_dat1)
+        frame2 <- add_total(frame2, total_dat2)
     }
 
     by_vars <- c(x_var, group_var, legend_var)
@@ -810,6 +831,31 @@ add_total_line <- function(plot_dat, x_var = names(plot_dat)[[1]],
     rbind(plot_dat, total_dat)
 }
 
+#' Convert the contributors to a total and append the total line
+#'
+#' Wraps the two steps that every plot of an array with `total = TRUE` needs:
+#' the contributors are put on the same size axis as the data they will join,
+#' and then summed there, see [add_total_line()].
+#'
+#' @param plot_dat The converted plotting data to append the total to.
+#' @param total_dat The unconverted contributors, see [total_contributors()].
+#' @param params A MizerParams object.
+#' @param size_axis Either `"w"` (weight) or `"l"` (length).
+#' @param x The mizer array the data came from, which says whether the values
+#'   are a density.
+#' @param per_log_size Whether to express a density per logarithmic size.
+#' @return `plot_dat` with the total appended as a series named `"Total"`.
+#' @keywords internal
+append_total_line <- function(plot_dat, total_dat, params, size_axis, x,
+                              per_log_size = NULL) {
+    total_dat <- convert_plot_density_axis(total_dat, params, size_axis,
+                                           density_wrt = array_density_wrt(x),
+                                           per_log_size = per_log_size)
+    total_dat <- add_total_line(total_dat)
+    total_dat <- total_dat[total_dat$Species == "Total", , drop = FALSE]
+    rbind(plot_dat, total_dat[, names(plot_dat), drop = FALSE])
+}
+
 #' The weight-length parameters to plot each row of plotting data with
 #'
 #' A length axis needs an allometric weight-length relationship for every line
@@ -1962,23 +2008,22 @@ plot_spectra <- function(params, n, n_pp,
     plot_dat <- convert_plot_spectrum_axis(plot_dat, params, size_axis, power,
                                            per_log_size = per_log_size)
     if (total) {
-        # The total is summed over the series as they are plotted, which on a
-        # length axis means at equal length rather than at equal weight. It is
-        # built from every species, whether or not it was selected for display,
-        # and from the resource when the resource is shown.
+        # The total is the total of everything the model holds: every species,
+        # whether or not it was selected for display, and the resource, whether
+        # or not it is drawn. It is summed over the series as they are plotted,
+        # which on a length axis means at equal length rather than at equal
+        # weight.
         # `n` has already been multiplied by the power of weight above.
-        total_dat <- data.frame(w = rep(w_grid, each = nrow(n)),
-                                value = c(unclass(n)),
-                                Species = rep(dimnames(n)[[1]],
-                                              times = length(w_grid)),
-                                Legend = "Total")
-        if (resource) {
-            resource_dat <- data.frame(w = w_full_grid,
-                                       value = c(n_pp * w_full_grid^power),
-                                       Species = "Resource",
-                                       Legend = "Total")
-            total_dat <- rbind(total_dat, resource_dat)
-        }
+        resource_dat <- data.frame(w = w_full_grid,
+                                   value = c(n_pp * w_full_grid^power),
+                                   Species = "Resource",
+                                   Legend = "Total")
+        total_dat <- rbind(data.frame(w = rep(w_grid, each = nrow(n)),
+                                      value = c(unclass(n)),
+                                      Species = rep(dimnames(n)[[1]],
+                                                    times = length(w_grid)),
+                                      Legend = "Total"),
+                           resource_dat)
         # Zeros are kept here, unlike in the plotted data: they are part of
         # each series and dropping them would let the interpolation bridge a
         # gap that is really empty.
diff --git a/inst/skills/analyse-and-plot/SKILL.md b/inst/skills/analyse-and-plot/SKILL.md
index 360273e4f..73bee0d23 100644
--- a/inst/skills/analyse-and-plot/SKILL.md
+++ b/inst/skills/analyse-and-plot/SKILL.md
@@ -224,7 +224,7 @@ dedicated `plot…()` functions below — share these optional arguments:
 | `wlim`/`llim` | numeric vector `c(min, max)` — restrict the size (x) axis |
 | `ylim` | numeric vector `c(min, max)` — restrict the value (y) axis |
 | `highlight` | character vector — draw named species with thicker lines |
-| `total` | logical — add a line for the community total |
+| `total` | logical — add a line for the community total. The total of *everything the object holds*, so it does not change when you select species or hide the resource; on a length axis it is summed at equal length |
 | `log_x`, `log_y` | logical — log-scale the x or y axis |
 | `size_axis` | `"w"` (default) or `"l"` — plot against weight or against length |
 
diff --git a/inst/skills/upgrade-mizer-code/SKILL.md b/inst/skills/upgrade-mizer-code/SKILL.md
index 21e1dacff..9b1db34d9 100644
--- a/inst/skills/upgrade-mizer-code/SKILL.md
+++ b/inst/skills/upgrade-mizer-code/SKILL.md
@@ -53,8 +53,7 @@ plots) are in the changelog and are not repeated here.
 | `array_spectrum_power()`, or `spectrum_power =` in an internal plot helper, is no longer found | replaced by the `type` metadata | Arrays say what kind of value they hold (3.3) |
 | New warning that two arrays hold "a value of type" different things in `plot2()` or `plotRelative()` | the two arrays disagree about what they hold | Arrays say what kind of value they hold (3.3) |
 | A `Total` line appears on a plot with `size_axis = "l"` where there used to be none | the total is now summed after the length conversion | The total is summed on the axis it is plotted against (3.3) |
-| `plotSpectra(total = TRUE, resource = FALSE)` gives a smaller total than before | the `resource` argument now controls the total as well as the line | The total is summed on the axis it is plotted against (3.3) |
-| `plot(<array>, total = TRUE)` gives a smaller total at some sizes | the total is now the sum of the lines drawn, not of the whole array | The total is summed on the axis it is plotted against (3.3) |
+| `plot(<array>, species = ..., total = TRUE)` gives a bigger total than before | the total is the total of the whole array, not of the selected species | The total is summed on the axis it is plotted against (3.3) |
 | New warning that a change to a species or resource parameter "has not taken effect" | the rate it feeds was set by hand and is no longer calculated | A change that cannot take effect now warns (3.3) |
 | That warning appears with `given_species_params<-()` but not with `species_params<-()` | the diagnostics belong to the given species parameter setter | The two species parameter setters divide the diagnostics between them (3.3) |
 | Setting a given species parameter to `NA` now warns that the change has not taken effect | clearing a value counts as a change, and a frozen array blocks it | The two species parameter setters divide the diagnostics between them (3.3) |
@@ -259,17 +258,16 @@ a length axis whenever the weight-length parameters agree — the union is that
 grid and the interpolation reproduces the values exactly. **The weight-axis
 total is unchanged**, for every power.
 
-Two things do change on the weight axis, both of them corrections:
-
-- `plotSpectra(total = TRUE, resource = FALSE)` now leaves the resource out of
-  the total. The `resource` argument used to control only whether the resource
-  line was drawn; the total counted it either way, so a "consumers only" total
-  was not available.
-- `plot(<array>, total = TRUE)` now sums the lines that are drawn rather than
-  the whole array, so it no longer counts abundance outside a species' size
-  range. On the North Sea model this moves the total at about 2% of sizes, by up
-  to a factor of three where a species holds abundance beyond its `w_max`. Pass
-  `all.sizes = TRUE` to sum the whole array as before.
+One thing does change on the weight axis. `total = TRUE` now means the same
+thing everywhere: **the total of everything the object holds**, whatever is
+drawn. `plotSpectra()` always worked that way — the resource and every species
+count, whether or not the resource is shown and whichever species were
+selected — and it still does. The array plots did not: `plot(<array>,
+total = TRUE)` summed only the species selected for display, and only the sizes
+inside each species' own range. It now sums the whole array, so the total no
+longer moves when you change `species`, `all.sizes` or `background`, and a plot
+of two species can be read against the community total. If you were relying on
+the total of a selection, sum the selection yourself.
 
 ### Length and weight parameters follow the one you gave last
 
diff --git a/man/append_total_line.Rd b/man/append_total_line.Rd
new file mode 100644
index 000000000..425a302f2
--- /dev/null
+++ b/man/append_total_line.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plots.R
+\name{append_total_line}
+\alias{append_total_line}
+\title{Convert the contributors to a total and append the total line}
+\usage{
+append_total_line(
+  plot_dat,
+  total_dat,
+  params,
+  size_axis,
+  x,
+  per_log_size = NULL
+)
+}
+\arguments{
+\item{plot_dat}{The converted plotting data to append the total to.}
+
+\item{total_dat}{The unconverted contributors, see \code{\link[=total_contributors]{total_contributors()}}.}
+
+\item{params}{A MizerParams object.}
+
+\item{size_axis}{Either \code{"w"} (weight) or \code{"l"} (length).}
+
+\item{x}{The mizer array the data came from, which says whether the values
+are a density.}
+
+\item{per_log_size}{Whether to express a density per logarithmic size.}
+}
+\value{
+\code{plot_dat} with the total appended as a series named \code{"Total"}.
+}
+\description{
+Wraps the two steps that every plot of an array with \code{total = TRUE} needs:
+the contributors are put on the same size axis as the data they will join,
+and then summed there, see \code{\link[=add_total_line]{add_total_line()}}.
+}
+\keyword{internal}
diff --git a/man/plotComparisonDataFrame.Rd b/man/plotComparisonDataFrame.Rd
index eb1b3534f..06fd35493 100644
--- a/man/plotComparisonDataFrame.Rd
+++ b/man/plotComparisonDataFrame.Rd
@@ -22,7 +22,7 @@ plotComparisonDataFrame(
   size_axis = NULL,
   density_wrt = NA_character_,
   per_log_size = NULL,
-  total = FALSE
+  total_dat = NULL
 )
 }
 \arguments{
@@ -59,8 +59,10 @@ and are left alone when the size axis changes.}
 \item{per_log_size}{Whether to express a density per logarithmic size.
 \code{NULL} (the default) keeps whichever the values already are.}
 
-\item{total}{Whether to add a line for the total over the series, formed
-after the size axis has been converted, see \code{\link[=add_total_line]{add_total_line()}}.}
+\item{total_dat}{Optional data frame of the contributors to a total, with a
+\code{Model} column identifying which of the two they belong to. The total is
+summed from them after the size axis has been converted, see
+\code{\link[=add_total_line]{add_total_line()}}.}
 }
 \value{
 A \code{mizer_plot} (ggplot2) object.
diff --git a/man/plotRelativeDataFrame.Rd b/man/plotRelativeDataFrame.Rd
index 5a596277d..09d5f0c29 100644
--- a/man/plotRelativeDataFrame.Rd
+++ b/man/plotRelativeDataFrame.Rd
@@ -15,7 +15,8 @@ plotRelativeDataFrame(
   highlight = NULL,
   legend_var = "Legend",
   size_axis = NULL,
-  total = FALSE
+  total_dat1 = NULL,
+  total_dat2 = NULL
 )
 }
 \arguments{
diff --git a/man/total_contributors.Rd b/man/total_contributors.Rd
new file mode 100644
index 000000000..8e7d166da
--- /dev/null
+++ b/man/total_contributors.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ArraySpeciesBySize-class.R
+\name{total_contributors}
+\alias{total_contributors}
+\title{Assemble the contributors to the total of a species-by-size array}
+\usage{
+total_contributors(x, wlim = c(NA, NA))
+}
+\arguments{
+\item{x}{An \code{ArraySpeciesBySize} object.}
+
+\item{wlim}{Numeric vector of length two giving the weight limits.}
+}
+\value{
+A data frame of plotting data holding every value in the array.
+}
+\description{
+The total is the total of everything the array holds: every species, whether
+or not it was selected for display, and every size, whether or not it falls
+in a species' own size range. It is a property of the array rather than of
+the plot, so that a plot of two species can still be read against the
+community total.
+}
+\details{
+The rows are returned unsummed, because the sum has to be taken after the
+size coordinate has been converted — on a length axis the species no longer
+share a grid; see \code{\link[=add_total_line]{add_total_line()}}.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-ArraySpeciesBySize-class.R b/tests/testthat/test-ArraySpeciesBySize-class.R
index 76d9d8f2f..bfcf9650c 100644
--- a/tests/testthat/test-ArraySpeciesBySize-class.R
+++ b/tests/testthat/test-ArraySpeciesBySize-class.R
@@ -293,12 +293,16 @@ test_that("the total is shown on a length axis and matches the lines drawn", {
     }, numeric(1))
     expect_equal(total_l[[2]], summed)
 
-    # The total is the sum of what is drawn, so trimming to each species' size
-    # range is respected; `all.sizes = TRUE` sums the whole array
-    all_sizes <- plot(density, total = TRUE, all.sizes = TRUE,
-                      return_data = TRUE)
-    expect_equal(all_sizes[[2]][all_sizes$Species == "Total"],
-                 unname(colSums(unclass(density))))
+    # The total is the total of everything in the array, so it depends on
+    # neither the species selection nor the size trimming
+    full <- unname(colSums(unclass(density)))
+    expect_equal(on_w[[2]][on_w$Species == "Total"], full)
+    for (args in list(list(all.sizes = TRUE), list(species = "Cod"),
+                      list(background = FALSE))) {
+        d <- do.call(plot, c(list(density, total = TRUE, return_data = TRUE),
+                             args))
+        expect_equal(d[[2]][d$Species == "Total"], full)
+    }
 })
 
 test_that("the total on a length axis interpolates onto a common grid", {
@@ -316,10 +320,13 @@ test_that("the total on a length axis interpolates onto a common grid", {
     expect_identical(sum(on_w$Species == "Total"), length(params@w))
     expect_gt(sum(on_l$Species == "Total"), length(params@w))
 
-    # At the largest length only the largest species is present, so the total
-    # is exactly that species there
+    # At the largest length only one species reaches, so the total is exactly
+    # that species there. The total covers the whole array, so compare it
+    # against the untrimmed lines.
     total_l <- on_l[on_l$Species == "Total", ]
-    lines_l <- on_l[on_l$Species != "Total", ]
+    all_l <- plot(density, all.sizes = TRUE, size_axis = "l",
+                  return_data = TRUE)
+    lines_l <- all_l[all_l$Species != "Total", ]
     at <- max(total_l$l)
     expect_identical(sum(lines_l$l == at), 1L)
     expect_equal(total_l[[2]][total_l$l == at], lines_l[[2]][lines_l$l == at])
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index dc8d82053..14409f3e8 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -718,16 +718,16 @@ test_that("the spectrum total appears on a length axis", {
     expect_equal(total[[2]][total$l == at],
                  sum(contributors[[2]][contributors$l == at]))
 
-    # The resource is in the total when it is shown and not when it is not,
-    # where it used to be counted either way
-    with_resource <- plotSpectra(params, total = TRUE, resource = TRUE,
-                                 return_data = TRUE)
-    without <- plotSpectra(params, total = TRUE, resource = FALSE,
-                           return_data = TRUE)
+    # The total is the total of everything the model holds, so it does not
+    # depend on what is drawn: neither hiding the resource nor selecting
+    # species changes it
     small <- min(params@w)
-    expect_gt(with_resource[[2]][with_resource$Legend == "Total" &
-                                     with_resource$w == small],
-              without[[2]][without$Legend == "Total" & without$w == small])
+    total_at <- function(...) {
+        d <- plotSpectra(params, total = TRUE, return_data = TRUE, ...)
+        d[[2]][d$Legend == "Total" & d$w == small]
+    }
+    expect_equal(total_at(resource = FALSE), total_at(resource = TRUE))
+    expect_equal(total_at(species = species), total_at())
 })
 
 test_that("validate_density_wrt accepts only the known measures", {
diff --git a/vignettes/cheatsheet-analysis-and-plotting.Rmd b/vignettes/cheatsheet-analysis-and-plotting.Rmd
index f3caf46ba..2d3ac8744 100644
--- a/vignettes/cheatsheet-analysis-and-plotting.Rmd
+++ b/vignettes/cheatsheet-analysis-and-plotting.Rmd
@@ -237,7 +237,7 @@ dedicated `plot…()` functions below — share these optional arguments:
 | `wlim`/`llim` | numeric vector `c(min, max)` — restrict the size (x) axis |
 | `ylim` | numeric vector `c(min, max)` — restrict the value (y) axis |
 | `highlight` | character vector — draw named species with thicker lines |
-| `total` | logical — add a line for the community total |
+| `total` | logical — add a line for the community total. The total of *everything the object holds*, so it does not change when you select species or hide the resource; on a length axis it is summed at equal length |
 | `log_x`, `log_y` | logical — log-scale the x or y axis |
 | `size_axis` | `"w"` (default) or `"l"` — plot against weight or against length |
 
diff --git a/vignettes/upgrading.Rmd b/vignettes/upgrading.Rmd
index 5fe105570..e2783bb71 100644
--- a/vignettes/upgrading.Rmd
+++ b/vignettes/upgrading.Rmd
@@ -140,17 +140,16 @@ a length axis whenever the weight-length parameters agree — the union is that
 grid and the interpolation reproduces the values exactly. **The weight-axis
 total is unchanged**, for every power.
 
-Two things do change on the weight axis, both of them corrections:
-
-- `plotSpectra(total = TRUE, resource = FALSE)` now leaves the resource out of
-  the total. The `resource` argument used to control only whether the resource
-  line was drawn; the total counted it either way, so a "consumers only" total
-  was not available.
-- `plot(<array>, total = TRUE)` now sums the lines that are drawn rather than
-  the whole array, so it no longer counts abundance outside a species' size
-  range. On the North Sea model this moves the total at about 2% of sizes, by up
-  to a factor of three where a species holds abundance beyond its `w_max`. Pass
-  `all.sizes = TRUE` to sum the whole array as before.
+One thing does change on the weight axis. `total = TRUE` now means the same
+thing everywhere: **the total of everything the object holds**, whatever is
+drawn. `plotSpectra()` always worked that way — the resource and every species
+count, whether or not the resource is shown and whichever species were
+selected — and it still does. The array plots did not: `plot(<array>,
+total = TRUE)` summed only the species selected for display, and only the sizes
+inside each species' own range. It now sums the whole array, so the total no
+longer moves when you change `species`, `all.sizes` or `background`, and a plot
+of two species can be read against the community total. If you were relying on
+the total of a selection, sum the selection yourself.
 
 ### Length and weight parameters follow the one you gave last
 

From ba8cb0a332a1dfcbb7154f459070cfb4b3f8665b Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Sun, 16 Aug 2026 08:03:35 +0100
Subject: [PATCH 3/4] Keep the total in the comparison spectra on a length axis

`plotSpectra2()` and `plotSpectraRelative()` assembled the two spectra on
the weight axis and converted afterwards, in the data frame helper. The
total they had been handed was therefore already summed at equal weight,
and reached the conversion with no species to convert it by, so it was
silently dropped: a length-axis comparison lost its total.

They now ask `plotSpectra()` for the axis they want, so the frames they
receive are already converted and the total in them is the total on the
axis being plotted, formed the same way as in a single spectrum plot.
The data frame helpers are then called without `size_axis`, so nothing
is converted twice.

The `data_ylim` treatment is unchanged: on a length axis the inner call
still gets `c(0, NA)`, because the user's limits apply to the converted
values and are imposed at the end.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 NEWS.md                                 |  7 ++++++
 R/plots.R                               | 21 ++++++++++-------
 inst/skills/upgrade-mizer-code/SKILL.md |  8 +++++++
 tests/testthat/test-plots.R             | 31 +++++++++++++++++++++++++
 vignettes/upgrading.Rmd                 |  7 ++++++
 5 files changed, 65 insertions(+), 9 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index ad31d30e3..fbe02ce2c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -171,6 +171,13 @@ stability of steady states.
   weight axis, the union is that grid and nothing is approximated: the
   weight-axis total is unchanged.
 
+  `plotSpectra2()` and `plotSpectraRelative()` keep their total on a length
+  axis too. They used to convert the axis after assembling the two spectra, so
+  the total they had been given — summed at equal weight — arrived at the
+  conversion with no species to convert it by and was dropped. They now let
+  `plotSpectra()` convert, so the total they receive is already the total on
+  the axis being plotted.
+
   `total = TRUE` also now means the same thing everywhere: the total of
   everything the object holds. For `plotSpectra()` that was already so — the
   resource and every species, whatever is drawn — and it stays so. The array
diff --git a/R/plots.R b/R/plots.R
index 3b1e6e850..4fb19661e 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -2501,19 +2501,22 @@ plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
     log_y <- log_axes$log_y
     assert_that(length(wlim) == 2, length(llim) == 2, length(ylim) == 2)
 
+    # The size axis is converted by `plotSpectra()` rather than afterwards, so
+    # that the total it forms is the total on the axis being plotted. A total
+    # summed at equal weight cannot be converted to a length axis afterwards:
+    # it is a sum over species that each convert differently.
     data_ylim <- if (identical(size_axis, "l")) c(0, NA) else ylim
     sf1 <- plotSpectra(object1, species = species, wlim = wlim,
                        ylim = data_ylim,
                        power = power, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     sf2 <- plotSpectra(object2, species = species, wlim = wlim,
                        ylim = data_ylim,
                        power = power, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
-    density_wrt <- spectrum_density_wrt(spectrum$per_log_size)
 
     plotComparisonDataFrame(sf1, sf2, validParams(params),
                             name1 = name1, name2 = name2,
@@ -2526,9 +2529,7 @@ plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
                             ytrans = if (log_y) "log10" else "identity",
                             xlim = plot_size_xlim(wlim, size_axis, llim),
                             ylim = ylim, highlight = highlight,
-                            legend_var = "Legend",
-                            size_axis = size_axis,
-                            density_wrt = density_wrt)
+                            legend_var = "Legend")
 }
 
 #' Resolve the power of weight multiplying a spectrum
@@ -2726,13 +2727,15 @@ plotSpectraRelative <- function(object1, object2,
     size_axis <- plot_size_axis(size_axis)
     assert_that(length(wlim) == 2, length(llim) == 2, length(ylim) == 2)
 
+    # As in `plotSpectra2()`, the conversion is left to `plotSpectra()` so that
+    # the total is summed on the axis being plotted.
     sf1 <- plotSpectra(object1, species = species, wlim = wlim,
                        power = 1, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     sf2 <- plotSpectra(object2, species = species, wlim = wlim,
                        power = 1, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
@@ -2742,7 +2745,7 @@ plotSpectraRelative <- function(object1, object2,
                           xlim = plot_size_xlim(wlim, size_axis, llim),
                           ylim = ylim,
                           highlight = highlight,
-                          legend_var = "Legend", size_axis = size_axis)
+                          legend_var = "Legend")
 }
 
 #' @rdname plotSpectraRelative
diff --git a/inst/skills/upgrade-mizer-code/SKILL.md b/inst/skills/upgrade-mizer-code/SKILL.md
index 9b1db34d9..a01930fd2 100644
--- a/inst/skills/upgrade-mizer-code/SKILL.md
+++ b/inst/skills/upgrade-mizer-code/SKILL.md
@@ -54,6 +54,7 @@ plots) are in the changelog and are not repeated here.
 | New warning that two arrays hold "a value of type" different things in `plot2()` or `plotRelative()` | the two arrays disagree about what they hold | Arrays say what kind of value they hold (3.3) |
 | A `Total` line appears on a plot with `size_axis = "l"` where there used to be none | the total is now summed after the length conversion | The total is summed on the axis it is plotted against (3.3) |
 | `plot(<array>, species = ..., total = TRUE)` gives a bigger total than before | the total is the total of the whole array, not of the selected species | The total is summed on the axis it is plotted against (3.3) |
+| `plotSpectra2()` or `plotSpectraRelative()` with `size_axis = "l"` shows a `Total` line where it used to show none | the total is now formed on the axis being plotted | The total is summed on the axis it is plotted against (3.3) |
 | New warning that a change to a species or resource parameter "has not taken effect" | the rate it feeds was set by hand and is no longer calculated | A change that cannot take effect now warns (3.3) |
 | That warning appears with `given_species_params<-()` but not with `species_params<-()` | the diagnostics belong to the given species parameter setter | The two species parameter setters divide the diagnostics between them (3.3) |
 | Setting a given species parameter to `NA` now warns that the change has not taken effect | clearing a value counts as a change, and a frozen array blocks it | The two species parameter setters divide the diagnostics between them (3.3) |
@@ -258,6 +259,13 @@ a length axis whenever the weight-length parameters agree — the union is that
 grid and the interpolation reproduces the values exactly. **The weight-axis
 total is unchanged**, for every power.
 
+`plotSpectra2()` and `plotSpectraRelative()` are fixed by the same change.
+They used to convert the size axis after assembling the two spectra, so the
+total they had been handed — already summed at equal weight — reached the
+conversion with no species to convert it by and was silently dropped. They now
+let `plotSpectra()` do the conversion, so the total they receive is the total
+on the axis being plotted.
+
 One thing does change on the weight axis. `total = TRUE` now means the same
 thing everywhere: **the total of everything the object holds**, whatever is
 drawn. `plotSpectra()` always worked that way — the resource and every species
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 14409f3e8..5a80fcf4c 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -730,6 +730,37 @@ test_that("the spectrum total appears on a length axis", {
     expect_equal(total_at(species = species), total_at())
 })
 
+test_that("the comparison spectra keep their total on a length axis", {
+    params2 <- params
+    species_params(params2)$gamma <- species_params(params2)$gamma * 1.1
+
+    for (size_axis in c("w", "l")) {
+        compared <- plotSpectra2(params, params2, total = TRUE,
+                                 size_axis = size_axis)$data
+        expect_true("Total" %in% compared$Legend)
+        relative <- plotSpectraRelative(params, params2, total = TRUE,
+                                        size_axis = size_axis)$data
+        expect_true("Total" %in% relative$Legend)
+    }
+
+    # Each model in the comparison holds exactly what plotSpectra() draws for
+    # it on its own, total included, on either axis
+    for (size_axis in c("w", "l")) {
+        for (power in 0:2) {
+            compared <- plotSpectra2(params, params, total = TRUE,
+                                     power = power,
+                                     size_axis = size_axis)$data
+            single <- plotSpectra(params, total = TRUE, power = power,
+                                  size_axis = size_axis, return_data = TRUE)
+            expect_equal(sort(compared[[2]][compared$Model == "First"]),
+                         sort(single[[2]]))
+            expect_equal(sort(compared[[2]][compared$Legend == "Total" &
+                                                compared$Model == "First"]),
+                         sort(single[[2]][single$Legend == "Total"]))
+        }
+    }
+})
+
 test_that("validate_density_wrt accepts only the known measures", {
     expect_identical(validate_density_wrt(NULL), NA_character_)
     expect_identical(validate_density_wrt(NA), NA_character_)
diff --git a/vignettes/upgrading.Rmd b/vignettes/upgrading.Rmd
index e2783bb71..900fb7d99 100644
--- a/vignettes/upgrading.Rmd
+++ b/vignettes/upgrading.Rmd
@@ -140,6 +140,13 @@ a length axis whenever the weight-length parameters agree — the union is that
 grid and the interpolation reproduces the values exactly. **The weight-axis
 total is unchanged**, for every power.
 
+`plotSpectra2()` and [`plotSpectraRelative()`](../reference/plotSpectraRelative.html) are fixed by the same change.
+They used to convert the size axis after assembling the two spectra, so the
+total they had been handed — already summed at equal weight — reached the
+conversion with no species to convert it by and was silently dropped. They now
+let `plotSpectra()` do the conversion, so the total they receive is the total
+on the axis being plotted.
+
 One thing does change on the weight axis. `total = TRUE` now means the same
 thing everywhere: **the total of everything the object holds**, whatever is
 drawn. `plotSpectra()` always worked that way — the resource and every species

From 1fed9873cd2f1f10ccafadfd5b23d1f47d14e193 Mon Sep 17 00:00:00 2001
From: Gustav Delius <gustav.delius@gmail.com>
Date: Sun, 16 Aug 2026 08:14:05 +0100
Subject: [PATCH 4/4] Apply `ylim` in the comparison spectra as everywhere else
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`plotSpectra()` applies `ylim` twice over: as the axis limits and as a
filter on the data, the latter carrying a hard floor at 1e-20 that keeps
a spectrum's numerically-zero tail from stretching a log axis.

`plotSpectra2()` could not do the same on a length axis. It assembled
the two spectra on the weight axis, so the values it would have filtered
were a Jacobian away from the ones the user's limits described. It
therefore passed `c(0, NA)` to the inner calls instead, disabling the
upper cut and lowering the floor to zero, and left the limits to the
scale alone.

The previous commit moved the conversion into those inner calls, which
removes the reason: the filter now sees exactly the values the limits
describe. So `ylim` is handed straight over, and the two functions agree
on both axes — with `ylim = c(median, NA)` on a length axis they now
return the same 512 rows, where the comparison used to return 1024 of
which half were below the limit.

The plot itself is unchanged, since the axis limits hid those points
either way. What changes is `return_data = TRUE`, and a series that dips
below the lower limit mid-range now breaks into segments rather than
running off the bottom — which is already how `plotSpectra()` draws it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
---
 NEWS.md                                 |  6 +++++-
 R/plots.R                               | 11 +++++------
 inst/skills/upgrade-mizer-code/SKILL.md | 10 ++++++++++
 tests/testthat/test-plots.R             | 10 +++++++++-
 vignettes/upgrading.Rmd                 |  9 +++++++++
 5 files changed, 38 insertions(+), 8 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index fbe02ce2c..b5f3fd03c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -176,7 +176,11 @@ stability of steady states.
   the total they had been given — summed at equal weight — arrived at the
   conversion with no species to convert it by and was dropped. They now let
   `plotSpectra()` convert, so the total they receive is already the total on
-  the axis being plotted.
+  the axis being plotted. As a result `plotSpectra2()` also applies `ylim` the
+  way `plotSpectra()` does on a length axis: it could not before, because the
+  values it was filtering were a Jacobian away from the ones the limits
+  described, so with `return_data = TRUE` it returned values outside the limits
+  that a single spectrum plot would have dropped.
 
   `total = TRUE` also now means the same thing everywhere: the total of
   everything the object holds. For `plotSpectra()` that was already so — the
diff --git a/R/plots.R b/R/plots.R
index 4fb19661e..210b6d163 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -2504,15 +2504,14 @@ plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
     # The size axis is converted by `plotSpectra()` rather than afterwards, so
     # that the total it forms is the total on the axis being plotted. A total
     # summed at equal weight cannot be converted to a length axis afterwards:
-    # it is a sum over species that each convert differently.
-    data_ylim <- if (identical(size_axis, "l")) c(0, NA) else ylim
-    sf1 <- plotSpectra(object1, species = species, wlim = wlim,
-                       ylim = data_ylim,
+    # it is a sum over species that each convert differently. It also means
+    # `ylim` can be handed straight over, because the values it is applied to
+    # are then the values it describes.
+    sf1 <- plotSpectra(object1, species = species, wlim = wlim, ylim = ylim,
                        power = power, total = total, resource = resource,
                        background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
-    sf2 <- plotSpectra(object2, species = species, wlim = wlim,
-                       ylim = data_ylim,
+    sf2 <- plotSpectra(object2, species = species, wlim = wlim, ylim = ylim,
                        power = power, total = total, resource = resource,
                        background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
diff --git a/inst/skills/upgrade-mizer-code/SKILL.md b/inst/skills/upgrade-mizer-code/SKILL.md
index a01930fd2..9df6d816e 100644
--- a/inst/skills/upgrade-mizer-code/SKILL.md
+++ b/inst/skills/upgrade-mizer-code/SKILL.md
@@ -55,6 +55,7 @@ plots) are in the changelog and are not repeated here.
 | A `Total` line appears on a plot with `size_axis = "l"` where there used to be none | the total is now summed after the length conversion | The total is summed on the axis it is plotted against (3.3) |
 | `plot(<array>, species = ..., total = TRUE)` gives a bigger total than before | the total is the total of the whole array, not of the selected species | The total is summed on the axis it is plotted against (3.3) |
 | `plotSpectra2()` or `plotSpectraRelative()` with `size_axis = "l"` shows a `Total` line where it used to show none | the total is now formed on the axis being plotted | The total is summed on the axis it is plotted against (3.3) |
+| `plotSpectra2(size_axis = "l", ylim = ..., return_data = TRUE)` returns fewer rows | `ylim` now filters the data there as it does for `plotSpectra()` | The total is summed on the axis it is plotted against (3.3) |
 | New warning that a change to a species or resource parameter "has not taken effect" | the rate it feeds was set by hand and is no longer calculated | A change that cannot take effect now warns (3.3) |
 | That warning appears with `given_species_params<-()` but not with `species_params<-()` | the diagnostics belong to the given species parameter setter | The two species parameter setters divide the diagnostics between them (3.3) |
 | Setting a given species parameter to `NA` now warns that the change has not taken effect | clearing a value counts as a change, and a frozen array blocks it | The two species parameter setters divide the diagnostics between them (3.3) |
@@ -266,6 +267,15 @@ conversion with no species to convert it by and was silently dropped. They now
 let `plotSpectra()` do the conversion, so the total they receive is the total
 on the axis being plotted.
 
+That also settles what `ylim` does there. `plotSpectra()` applies it both as
+the axis limits and as a filter on the data, with a hard floor at 1e-20.
+`plotSpectra2()` could not do the same on a length axis, because the values it
+was filtering were a Jacobian away from the ones the limits described, so it
+skipped the filter. Now that it converts first, the filter applies as it does
+everywhere else. The plot is unchanged — the axis limits hid those points
+anyway — but `return_data = TRUE` no longer hands back values outside the
+limits.
+
 One thing does change on the weight axis. `total = TRUE` now means the same
 thing everywhere: **the total of everything the object holds**, whatever is
 drawn. `plotSpectra()` always worked that way — the resource and every species
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 5a80fcf4c..803c3720b 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -616,12 +616,20 @@ test_that("length-axis spectra transform their densities", {
         compared_data$l
     expect_equal(compared_data[[2]], unname(by_weight[[2]] * jacobian))
 
+    # `ylim` is applied to the converted values, not to the weight-axis ones it
+    # would be off by a Jacobian from, and it filters the data exactly as it
+    # does for a single spectrum plot.
     limit <- stats::median(compared_data[[2]])
     limited <- plotSpectra2(params, params, species = species,
                             resource = FALSE, power = 0, size_axis = "l",
                             ylim = c(limit, NA))
-    expect_true(any(limited$data[[2]] < limit))
+    expect_true(all(limited$data[[2]] >= limit))
     expect_equal(limited$scales$get_scales("y")$limits, c(log10(limit), NA))
+    single <- plotSpectra(params, species = species, resource = FALSE,
+                          power = 0, size_axis = "l", ylim = c(limit, NA),
+                          return_data = TRUE)
+    expect_equal(sort(limited$data[[2]][limited$data$Model == "First"]),
+                 sort(single[[2]]))
 })
 
 test_that("proportion_ylim widens to [0, 1] without ever hiding data", {
diff --git a/vignettes/upgrading.Rmd b/vignettes/upgrading.Rmd
index 900fb7d99..0017e5559 100644
--- a/vignettes/upgrading.Rmd
+++ b/vignettes/upgrading.Rmd
@@ -147,6 +147,15 @@ conversion with no species to convert it by and was silently dropped. They now
 let `plotSpectra()` do the conversion, so the total they receive is the total
 on the axis being plotted.
 
+That also settles what `ylim` does there. `plotSpectra()` applies it both as
+the axis limits and as a filter on the data, with a hard floor at 1e-20.
+`plotSpectra2()` could not do the same on a length axis, because the values it
+was filtering were a Jacobian away from the ones the limits described, so it
+skipped the filter. Now that it converts first, the filter applies as it does
+everywhere else. The plot is unchanged — the axis limits hid those points
+anyway — but `return_data = TRUE` no longer hands back values outside the
+limits.
+
 One thing does change on the weight axis. `total = TRUE` now means the same
 thing everywhere: **the total of everything the object holds**, whatever is
 drawn. `plotSpectra()` always worked that way — the resource and every species