Sum the total on the axis it is plotted against

NEWS.md

@@ -159,6 +159,36 @@ 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.
+
+  `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. 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
+  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
   observed yield are then taken from those gears only, so in a model where

R/ArraySpeciesBySize-class.R

@@ -364,6 +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 <- 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)
     }
@@ -557,7 +561,13 @@ 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_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
@@ -675,7 +685,9 @@ 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_dat1 = if (total) total_contributors(x, wlim),
+                          total_dat2 = if (total) total_contributors(y, wlim))
 }
 
 check_plot2_compatible <- function(x, y, class) {
@@ -973,6 +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 <- 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)
     }
@@ -1127,11 +1143,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,23 +1176,32 @@ 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
 
     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

R/ArrayTimeBySpeciesBySize-class.R

@@ -434,12 +434,15 @@ 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
+    # 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_sums <- stats::aggregate(value ~ time + w, data = df, FUN = sum,
-                                       na.rm = TRUE)
-        total_sums$Species <- "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
@@ -452,10 +455,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 +465,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)

R/animateSpectra.R

@@ -197,17 +197,19 @@ 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()`. 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) {
-        # 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"
+        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 +222,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 +250,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) {