R/ArraySpeciesBySize-class.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

# ArraySpeciesBySize S3 class for species x size arrays
#
# Copyright 2026 Gustav Delius.
# Distributed under the GPL 3 or later.

#' S3 class for species x size rate arrays
#'
#' Many functions in mizer return two-dimensional arrays (species x size)
#' holding rates like encounter rate, feeding level, growth rate, mortality etc.
#' The `ArraySpeciesBySize` class wraps these arrays to provide convenient
#' `print()`, `summary()`, `plot()`, and `as.data.frame()` methods.
#'
#' An `ArraySpeciesBySize` object behaves just like a regular matrix for
#' arithmetic operations and subsetting. It carries a few lightweight
#' attributes:
#' \itemize{
#'   \item `value_name` – a human-readable name for the value
#'       (e.g. "Encounter rate").
#'   \item `units` – the units of the rate (e.g. "g/year").
#'   \item `type` – the kind of quantity the values are.
#' }
#'
#' @param x A matrix (species x size). For `is.ArraySpeciesBySize()`, any
#'   object to test.
#' @param value_name A string giving the human-readable name for the value.
#' @param units A string giving the units (e.g. "g/year", "1/year").
#' @param type The kind of quantity the values are, see [array_types]:
#'   `"value"` (the default) for a rate or an amount, `"density"` for an amount
#'   per gram of body weight, `"proportion"` for a fraction. This is what tells
#'   `plot()` to multiply a density by the appropriate Jacobian when it is
#'   plotted against a length axis (`size_axis = "l"`), and to show a proportion
#'   against the whole of the interval from 0 to 1. The default, `NULL`, treats
#'   a `value_name` of `"Number density"` or units of `"1/g"` as a density, the
#'   way mizer recognised one before this attribute existed.
#' @param params A `MizerParams` object. Used for species colours, linetypes,
#'   and size ranges in the `plot()` method.
#' @param representation Either `"point"` (the default) for a quantity sampled
#'   at the grid nodes, or `"average"` for a finite-volume bin average. A
#'   bin-averaged quantity is drawn at the geometric bin centre rather than the
#'   left bin edge, but only when the model uses second-order bin-averaging
#'   (`second_order_w[["bin_average"]]`), so default plots are unchanged.
#'
#' @return An `ArraySpeciesBySize` object (inherits from `matrix` and `array`).
#' @seealso [print()], [summary()], [as.data.frame()], [plot()]
#' @export
#' @examples
#' \donttest{
#' enc <- getEncounter(NS_params)
#' is.ArraySpeciesBySize(enc)
#' summary(enc)
#' }
ArraySpeciesBySize <- function(x, value_name = NULL, units = NULL,
                               type = NULL,
                               params = NULL,
                               representation = c("point", "average")) {
    if (!is.matrix(x)) {
        stop("`x` must be a matrix.")
    }
    representation <- match.arg(representation)
    type <- resolve_array_type(type, value_name, units)
    if (!is.null(params) && identical(dim(x), dim(params@metab))) {
        dimnames(x) <- dimnames(params@metab)
    }
    structure(x,
        class = c("ArraySpeciesBySize", "matrix", "array"),
        value_name = value_name,
        units = units,
        type = type,
        params = params,
        representation = representation
    )
}

#' @rdname ArraySpeciesBySize
#' @return `is.ArraySpeciesBySize()` returns `TRUE` if `x` is an
#'   `ArraySpeciesBySize` object, `FALSE` otherwise.
#' @export
is.ArraySpeciesBySize <- function(x) {
    inherits(x, "ArraySpeciesBySize")
}

#' @export
print.ArraySpeciesBySize <- function(x, ...) {
    value_name <- attr(x, "value_name") %||% "ArraySpeciesBySize"
    units_str <- attr(x, "units")
    dims <- dim(x)
    header <- paste0(value_name, " (", dims[1], " species x ", dims[2],
                     " sizes)")
    if (!is.null(units_str)) {
        header <- paste0(header, " [", units_str, "]")
    }
    cat(header, "\n")
    print_ArraySpeciesBySize_body(x)
    invisible(x)
}

# Print the (possibly truncated) species x size body of an ArraySpeciesBySize
# object, without the header line. Species are truncated to a leading subset,
# sizes to an evenly log-spaced, width-fitted sample. Shared with
# print.ArrayTimeBySpeciesBySize, which prints this body for a single time
# slice.
print_ArraySpeciesBySize_body <- function(x) {
    w <- get_ArraySpeciesBySize_w(x)
    mat <- unclass_rate(x)
    n_species <- nrow(mat)
    n_sizes <- ncol(mat)

    sp_idx <- pick_head_indices(n_species, mizer_print_defaults$species_head,
                                mizer_print_defaults$species_threshold)

    size_k <- fit_log_spaced_k(
        n_sizes, mizer_print_defaults$size_max, mizer_print_defaults$size_min,
        width_fn = function(k) {
            sz_idx <- pick_log_spaced_indices(n_sizes, k)
            matrix_display_width(mat[sp_idx, sz_idx, drop = FALSE])
        })
    sz_idx <- pick_log_spaced_indices(n_sizes, size_k)

    print(mat[sp_idx, sz_idx, drop = FALSE])

    if (length(sp_idx) < n_species) {
        omitted <- setdiff(rownames(mat), rownames(mat)[sp_idx])
        detail <- paste(utils::head(omitted, 5), collapse = ", ")
        if (length(omitted) > 5) detail <- paste0(detail, ", ...")
        cat(format_truncation_note(length(sp_idx), n_species, "species", detail), "\n")
    }
    if (length(sz_idx) < n_sizes) {
        cat(format_truncation_note(length(sz_idx), n_sizes, "sizes",
                                   format_size_range_detail(w)), "\n")
    }
    invisible(NULL)
}

#' @export
summary.ArraySpeciesBySize <- function(object, ...) {
    value_name <- attr(object, "value_name") %||% "ArraySpeciesBySize"
    units_str <- attr(object, "units")
    sp_names <- rownames(object)
    mat <- unclass(object)

    df <- data.frame(
        Species = sp_names,
        Min = apply(mat, 1, min, na.rm = TRUE),
        Mean = apply(mat, 1, mean, na.rm = TRUE),
        Max = apply(mat, 1, max, na.rm = TRUE),
        row.names = NULL,
        stringsAsFactors = FALSE
    )

    result <- list(
        value_name = value_name,
        units = units_str,
        dims = dim(object),
        per_species = df
    )
    class(result) <- "summary.ArraySpeciesBySize"
    result
}

#' @export
print.summary.ArraySpeciesBySize <- function(x, ...) {
    header <- x$value_name
    if (!is.null(x$units)) {
        header <- paste0(header, " [", x$units, "]")
    }
    cat(header, "\n")
    cat(x$dims[1], "species x", x$dims[2], "sizes\n\n")
    print(x$per_species, row.names = FALSE)
    invisible(x)
}

#' Plot mizer arrays
#'
#' Many mizer functions return values that depend on species and either size or
#' time. `plot()` creates a ggplot2 figure with one line for each species
#' showing the values against size or against time (depending on the type of
#' output). [plotHover()] creates an interactive version of the same figure.
#'
#' This works because the mizer functions that give values that depend on
#' species and size return an `ArraySpeciesBySize` object and those that
#' give values that depend on species and time return an `ArrayTimeBySpecies`
#' object. These objects have attributes that store the name of the value,
#' its units, and a reference to the `MizerParams` object that the value was
#' computed from. This allows the plots to be automatically labelled and
#' coloured appropriately.
#'
#' The resource classes `ArrayResourceBySize` and `ArrayTimeByResourceBySize`
#' work the same way, except that they hold a single spectrum rather than one
#' per species.
#'
#' To compare two mizer arrays in a single plot, use [plot2()]. To show the
#' relative difference between two arrays, use [plotRelative()]. To add an array
#' to an existing plot, use [addPlot()]. All three, and [animate()], have
#' methods for every mizer array class.
#'
#' All methods return a ggplot2 object, unless `return_data = TRUE`, in
#' which case they return the underlying data frame instead. [plotHover()]
#' returns a plotly object.
#'
#' Arguments used by all methods:
#' \describe{
#'   \item{`species`}{Character vector of species to include. `NULL`
#'     (default) means all species.}
#'   \item{`highlight`}{Name or vector of names of the species to be
#'     highlighted.}
#'   \item{`total`}{A boolean value that determines whether the total is
#'     plotted as well. The total is the total of everything the array holds,
#'     every species and every size, whatever is drawn. Default is `FALSE`.}
#'   \item{`background`}{A boolean value that determines whether background
#'     species are included. Ignored if the model does not contain background
#'     species. Default is `TRUE`.}
#'   \item{`return_data`}{If `TRUE`, return the data frame instead of the
#'     plot.}
#'   \item{`log_x`}{If `TRUE`, use a log10 x-axis. The default depends on the
#'     method; see its own help page.}
#'   \item{`log_y`}{If `TRUE`, use a log10 y-axis. The default depends on the
#'     method; see its own help page.}
#'   \item{`log`}{Character string specifying which axes should use log10
#'     scales, in the same form as the base [plot()] argument. For example,
#'     `"x"`, `"y"`, `"xy"` or `""`. If supplied, this overrides `log_x` and
#'     `log_y`.}
#'   \item{`ylim`}{A numeric vector of length two providing lower and upper
#'     limits for the value (y) axis. Use `NA` to refer to the existing
#'     minimum or maximum.}
#'   \item{`y_ticks`}{The approximate number of ticks desired on the y axis.}
#' }
#'
#' Additional arguments for [plot.ArraySpeciesBySize()] and
#' [plot.ArrayTimeBySpeciesBySize()]:
#' \describe{
#'   \item{`all.sizes`}{If `FALSE` (default), values outside a species' size
#'     range (`w_min` to `w_max`) are removed.}
#'   \item{`wlim`}{A numeric vector of length two providing lower and upper
#'     limits for the weight (x) axis. Use `NA` to refer to the existing
#'     minimum or maximum.}
#'   \item{`llim`}{A numeric vector of length two providing lower and upper
#'     limits for the length (x) axis when `size_axis = "l"`. Use `NA` to
#'     refer to the existing minimum or maximum.}
#'   \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or
#'     length (`"l"`), using the allometric weight-length relationship.
#'     Densities are transformed to match the chosen axis.}
#'   \item{`per_log_size`}{For an array that holds a density, whether to plot it
#'     per logarithmic size (`TRUE`) rather than per size (`FALSE`). The
#'     default, `NULL`, plots the density as it stands. Unlike `size_axis` this
#'     needs no weight-length relationship, so it is available for the resource
#'     classes too. An error for an array that does not hold a density.}
#' }
#'
#' Additional argument for [plot.ArrayTimeBySpecies()]:
#' \describe{
#'   \item{`tlim`}{A numeric vector of length two providing lower and upper
#'     limits for the time axis, e.g. `c(1980, 2000)`. Use `NA` to apply no
#'     limit at that end. Default is `c(NA, NA)`.}
#' }
#'
#' Additional argument for [plot.ArrayTimeBySpeciesBySize()] and
#' [plot.ArrayTimeByResourceBySize()]:
#' \describe{
#'   \item{`time`}{The time to display. Default (`NULL`) is the final time
#'     step.}
#' }
#'
#' See the individual method help pages for each method's exact arguments and
#' defaults: [plot.ArraySpeciesBySize()], [plot.ArrayTimeBySpecies()],
#' [plot.ArrayTimeBySpeciesBySize()], [plot.ArrayResourceBySize()],
#' [plot.ArrayTimeByResourceBySize()].
#'
#' @name plot
#' @family plotting functions
#' @examples
#' \donttest{
#' plot(getEncounter(NS_params))
#' plot(getFeedingLevel(NS_params), species = c("Cod", "Herring"))
#' plot(getPredMort(NS_params), species = c("Cod", "Herring"),
#'      size_axis = "l")
#' plot(getBiomass(NS_sim))
#' plot(getBiomass(NS_sim), species = c("Cod", "Herring"), total = TRUE)
#' plot(getYield(NS_sim), species = c("Cod", "Herring"))
#' plot(getFMort(NS_sim), time = 2010)
#' plot(getResourceMort(NS_params))
#' plot(initialNResource(NS_params))
#' plot(NResource(NS_sim))
#' }
NULL

#' Plot method for `ArraySpeciesBySize` objects
#'
#' See [plot()] for an overview of the mizer plotting system and the
#' arguments shared by all of its methods.
#'
#' @param x An `ArraySpeciesBySize` object.
#' @param species Character vector of species to include. `NULL`
#'   (default) means all species.
#' @param all.sizes If `FALSE` (default), values outside a species' size
#'   range (`w_min` to `w_max`) are removed.
#' @param highlight Name or vector of names of the species to be
#'   highlighted.
#' @param return_data If `TRUE`, return the data frame instead of the
#'   plot.
#' @param log_x If `TRUE`, use a log10 x-axis. Default is `TRUE`.
#' @param log_y If `TRUE`, use a log10 y-axis. Default is `FALSE`.
#' @param log Character string specifying which axes should use log10
#'   scales, in the same form as the base [plot()] argument. For example,
#'   `"x"`, `"y"`, `"xy"` or `""`. If supplied, this overrides `log_x` and
#'   `log_y`.
#' @param wlim A numeric vector of length two providing lower and upper
#'   limits for the weight (x) axis. Use `NA` to refer to the existing
#'   minimum or maximum.
#' @param llim A numeric vector of length two providing lower and upper
#'   limits for the length (x) axis when `size_axis = "l"`. Use `NA` to
#'   refer to the existing minimum or maximum.
#' @param ylim A numeric vector of length two providing lower and upper
#'   limits for the value (y) axis. Use `NA` to refer to the existing
#'   minimum or maximum.
#' @param size_axis Whether to plot size as weight (`"w"`, default) or
#'   length (`"l"`), using the allometric weight-length relationship.
#' @param per_log_size For an array that holds a density, whether to plot it per
#'   logarithmic size (`TRUE`) rather than per size (`FALSE`). The default,
#'   `NULL`, plots the density as it stands. Unlike `size_axis` this needs no
#'   weight-length relationship, so it is available for the resource classes
#'   too. An error for an array that does not hold a density.
#' @param total A boolean value that determines whether the total is plotted
#'   as well. The total is the total of everything the array holds, every
#'   species and every size, whatever is drawn. Default is `FALSE`.
#' @param background A boolean value that determines whether background
#'   species are included. Ignored if the model does not contain background
#'   species. Default is `TRUE`.
#' @param y_ticks The approximate number of ticks desired on the y axis.
#' @param ... Unused.
#'
#' @return A ggplot2 object, unless `return_data = TRUE`, in which case a
#'   data frame is returned.
#' @keywords internal
#' @export
#' @examples
#' \donttest{
#' plot(getEncounter(NS_params))
#' plot(getFeedingLevel(NS_params), species = c("Cod", "Herring"))
#' plot(getPredMort(NS_params), species = c("Cod", "Herring"),
#'      size_axis = "l")
#' }
plot.ArraySpeciesBySize <- function(x, species = NULL,
                            all.sizes = FALSE, highlight = NULL,
                            return_data = FALSE, log_x = TRUE, log_y = FALSE,
                            log = NULL,
                            wlim = c(NA, NA), llim = c(NA, NA),
                            ylim = c(NA, NA),
                            size_axis = c("w", "l"),
                            per_log_size = NULL,
                            total = FALSE, background = TRUE,
                            y_ticks = 6, ...) {
    size_axis <- plot_size_axis(size_axis)
    check_per_log_size(x, per_log_size)
    log_y <- array_log_y(x, log_y, log, !missing(log_y))
    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
    log_x <- log_axes$log_x
    log_y <- log_axes$log_y

    assert_that(length(wlim) == 2,
                length(llim) == 2,
                length(ylim) == 2)
    params <- attr(x, "params")
    plot_dat <- ArraySpeciesBySize_plot_data(
        x, species = species, all.sizes = all.sizes, wlim = wlim, llim = llim,
        total = total, background = background, size_axis = size_axis,
        per_log_size = per_log_size)

    if (return_data) return(plot_dat)

    ylim <- array_ylim(x, ylim, log_y, plot_dat[[2]])
    y_label <- array_y_label(x, default = "Rate", size_axis = size_axis,
                             per_log_size = per_log_size)

    plotDataFrame(plot_dat, params, xlab = plot_size_xlab(size_axis),
                  ylab = y_label,
                  xtrans = if (log_x) "log10" else "identity",
                  ytrans = if (log_y) "log10" else "identity",
                  xlim = plot_size_xlim(wlim, size_axis, llim), ylim = ylim,
                  highlight = highlight, y_ticks = y_ticks,
                  legend_var = "Legend")
}

#' Parse the log-axis arguments of a mizer plot function
#'
#' Internal helper that resolves the various ways of specifying which axes
#' should use a logarithmic scale into a consistent pair of logical flags. It
#' is exported so that extension packages (such as mizerMR) can reuse it in
#' their own array `plot()` methods.
#'
#' @param log Either `NULL`, a single logical (legacy form, toggling only the
#'   y-axis), or a character string containing only the letters `"x"` and/or
#'   `"y"` to indicate which axes should be logarithmic.
#' @param log_x,log_y Default logical flags used when `log` is `NULL`.
#'
#' @return A list with logical components `log_x` and `log_y`.
#' @keywords internal
#' @export
parsePlotLog <- function(log, log_x = FALSE, log_y = FALSE) {
    if (is.null(log)) {
        return(list(log_x = log_x, log_y = log_y))
    }
    # Backward compatibility: legacy logical `log` toggles only the y-axis.
    if (is.logical(log)) {
        if (length(log) != 1 || is.na(log)) {
            stop("`log` must be a single logical value or a character string ",
                 "containing only \"x\" and/or \"y\".")
        }
        return(list(log_x = FALSE, log_y = isTRUE(log)))
    }
    if (!is.character(log) || length(log) != 1 || is.na(log) ||
        grepl("[^xy]", log)) {
        stop("`log` must be a single logical value or a character string ",
             "containing only \"x\" and/or \"y\".")
    }
    list(
        log_x = grepl("x", log, fixed = TRUE),
        log_y = grepl("y", log, fixed = TRUE)
    )
}

#' Compare two mizer arrays in a single plot
#'
#' `plot2()` compares two compatible mizer array objects in a single ggplot.
#' Colours identify species or groups, and linetype identifies which object
#' the values came from.
#'
#' @param x The first of two compatible mizer array objects to compare.
#'   Can be an `ArraySpeciesBySize`, `ArrayTimeBySpecies`,
#'   `ArrayTimeBySpeciesBySize`, `ArrayResourceBySize` or
#'   `ArrayTimeByResourceBySize` object.
#' @param y The second mizer array object, compatible with `x`.
#' @param name1,name2 Labels for the two objects, used in the linetype legend.
#' @param species Character vector of species to include. `NULL` (default) means
#'   all species. A resource array holds a single spectrum, so this argument is
#'   not used by the resource methods, which warn if it is set.
#' @param log_x If `TRUE`, use a log10 x-axis. Default is `TRUE` for size
#'   spectra and `FALSE` for time series.
#' @param log_y If `TRUE`, use a log10 y-axis. Default is `FALSE` for
#'   `ArraySpeciesBySize` and `TRUE` for `ArrayTimeBySpecies` and for the
#'   resource classes.
#' @param log Character string specifying which axes should use log10 scales,
#'   in the same form as the base [plot()] argument. For example, `"x"`,
#'   `"y"`, `"xy"` or `""`. If supplied, this overrides `log_x` and `log_y`.
#' @param ylim A numeric vector of length two providing lower and upper limits
#'   for the value (y) axis. Use `NA` to refer to the existing minimum or
#'   maximum.
#' @param total A boolean value that determines whether the total is plotted
#'   as well. The total is the total of everything the array holds, every
#'   species and every size, whatever is drawn. Default is `FALSE`. Not used by
#'   the resource methods, which warn if it is set.
#' @param background A boolean value that determines whether background species
#'   are included. Ignored if the model does not contain background species.
#'   Default is `TRUE`. Not used by the resource methods, which warn if it is
#'   set.
#' @param highlight Name or vector of names of the species to be highlighted
#'   with a thicker line.
#' @param y_ticks The approximate number of ticks desired on the y axis.
#' @param ... Further arguments used by only some of the methods:
#'
#'   **For the `ArraySpeciesBySize`, `ArrayTimeBySpeciesBySize`,
#'   `ArrayResourceBySize` and `ArrayTimeByResourceBySize` methods:**
#'   \describe{
#'     \item{`wlim`}{A numeric vector of length two providing lower and upper
#'       limits for the weight (x) axis. Use `NA` to refer to the existing
#'       minimum or maximum.}
#'   }
#'
#'   **For the `ArraySpeciesBySize` and `ArrayTimeBySpeciesBySize` methods:**
#'   \describe{
#'     \item{`all.sizes`}{If `FALSE` (default), values outside a species' size
#'       range (`w_min` to `w_max`) are removed.}
#'     \item{`llim`}{A numeric vector of length two providing lower and upper
#'       limits for the length (x) axis when `size_axis = "l"`. Use `NA` to
#'       refer to the existing minimum or maximum.}
#'     \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or
#'       length (`"l"`), using the allometric weight-length relationship of
#'       each species, or of the resource, see [resource_params()].}
#'     \item{`per_log_size`}{For an array that holds a density, whether to plot
#'       it per logarithmic size (`TRUE`) rather than per size (`FALSE`). The
#'       default, `NULL`, plots the density as it stands. Unlike `size_axis`
#'       this needs no weight-length relationship, so it is available for the
#'       resource classes too. An error for an array that does not hold a
#'       density.}
#'   }
#'
#'   **For `ArrayTimeBySpecies` methods:**
#'   \describe{
#'     \item{`tlim`}{A numeric vector of length two providing lower and upper
#'       limits for the time axis, e.g. `c(1980, 2000)`. Use `NA` to apply no
#'       limit at that end. Default is `c(NA, NA)`.}
#'   }
#'
#'   **For the `ArrayTimeBySpeciesBySize` and `ArrayTimeByResourceBySize`
#'   methods:**
#'   \describe{
#'     \item{`time`}{The time to display. Default (`NULL`) is the final time
#'       step.}
#'   }
#'
#' @return A ggplot2 object.
#'
#' @family plotting functions
#' @export
#' @examples
#' \donttest{
#' plot2(getEncounter(NS_params), getEncounter(NS_params))
#' plot2(getResourceMort(NS_params), getResourceMort(NS_params))
#' }
plot2 <- function(x, y, name1 = "First", name2 = "Second",
                  species = NULL, log_x, log_y, log = NULL,
                  ylim = c(NA, NA), total = FALSE, background = TRUE,
                  highlight = NULL, y_ticks = 6, ...) {
    UseMethod("plot2", x)
}

#' @rdname plot2
#' @usage NULL
#' @export
plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
                                     species = NULL,
                                     log_x = TRUE, log_y = FALSE, log = NULL,
                                     ylim = c(NA, NA),
                                     total = FALSE, background = TRUE,
                                     highlight = NULL,
                                     y_ticks = 6,
                                     all.sizes = FALSE,
                                     wlim = c(NA, NA), llim = c(NA, NA),
                                     size_axis = c("w", "l"),
                                     per_log_size = NULL, ...) {
    check_plot2_compatible(x, y, "ArraySpeciesBySize")
    compare_array_metadata(x, y)
    size_axis <- plot_size_axis(size_axis)
    check_per_log_size(x, per_log_size)
    log_y <- array_log_y(x, log_y, log, !missing(log_y))
    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
    log_x <- log_axes$log_x
    log_y <- log_axes$log_y
    assert_that(length(wlim) == 2,
                length(llim) == 2,
                length(ylim) == 2)

    params <- attr(x, "params")
    y_label <- array_y_label(x, default = "Rate", size_axis = size_axis,
                             per_log_size = per_log_size)
    # Each array is prepared with its own model, so that a length axis and a
    # density Jacobian use the weight-length relationship the values belong to.
    plot_dat1 <- ArraySpeciesBySize_plot_data(
        x, species = species, all.sizes = all.sizes, wlim = wlim, llim = llim,
        total = total, background = background, size_axis = size_axis,
        per_log_size = per_log_size)
    plot_dat2 <- ArraySpeciesBySize_plot_data(
        y, species = species, all.sizes = all.sizes, wlim = wlim, llim = llim,
        total = total, background = background, size_axis = size_axis,
        per_log_size = per_log_size)

    ylim <- array_ylim(x, ylim, log_y, c(plot_dat1[[2]], plot_dat2[[2]]))

    plotComparisonDataFrame(plot_dat1, plot_dat2, params,
                            name1 = name1, name2 = name2,
                            xlab = plot_size_xlab(size_axis), ylab = y_label,
                            xtrans = if (log_x) "log10" else "identity",
                            ytrans = if (log_y) "log10" else "identity",
                            xlim = plot_size_xlim(wlim, size_axis, llim),
                            ylim = ylim, highlight = highlight,
                            y_ticks = y_ticks, legend_var = "Legend")
}

#' Plot relative difference between two mizer arrays
#'
#' `plotRelative()` plots the difference between two compatible mizer array
#' objects relative to their average. If the values in the first object are
#' \eqn{N_1} and the values in the second are \eqn{N_2}, it plots
#' \deqn{2 (N_2 - N_1) / (N_1 + N_2).}
#'
#' @param x The first of two compatible mizer array objects to compare.
#'   Can be an `ArraySpeciesBySize`, `ArrayTimeBySpecies`,
#'   `ArrayTimeBySpeciesBySize`, `ArrayResourceBySize` or
#'   `ArrayTimeByResourceBySize` object.
#' @param y The second mizer array object, compatible with `x`.
#' @param species Character vector of species to include. `NULL` (default) means
#'   all species. A resource array holds a single spectrum, so this argument is
#'   not used by the resource methods, which warn if it is set.
#' @param log_x If `TRUE`, use a log10 x-axis. Default is `TRUE` for size
#'   spectra and `FALSE` for time series.
#' @param ylim A numeric vector of length two providing lower and upper limits
#'   for the value (y) axis.
#' @param total A boolean value that determines whether the total is plotted
#'   as well. The total is the total of everything the array holds, every
#'   species and every size, whatever is drawn. Default is `FALSE`. Not used by
#'   the resource methods, which warn if it is set.
#' @param background A boolean value that determines whether background species
#'   are included. Ignored if the model does not contain background species.
#'   Default is `TRUE`. Not used by the resource methods, which warn if it is
#'   set.
#' @param highlight Name or vector of names of the species to be highlighted
#'   with a thicker line.
#' @param ... Further arguments used by only some of the methods:
#'
#'   **For the `ArraySpeciesBySize`, `ArrayTimeBySpeciesBySize`,
#'   `ArrayResourceBySize` and `ArrayTimeByResourceBySize` methods:**
#'   \describe{
#'     \item{`wlim`}{A numeric vector of length two providing lower and upper
#'       limits for the weight (x) axis. Use `NA` to refer to the existing
#'       minimum or maximum.}
#'   }
#'
#'   **For the `ArraySpeciesBySize` and `ArrayTimeBySpeciesBySize` methods:**
#'   \describe{
#'     \item{`all.sizes`}{If `FALSE` (default), values outside a species' size
#'       range (`w_min` to `w_max`) are removed.}
#'     \item{`llim`}{A numeric vector of length two providing lower and upper
#'       limits for the length (x) axis when `size_axis = "l"`. Use `NA` to
#'       refer to the existing minimum or maximum.}
#'     \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or
#'       length (`"l"`), using the allometric weight-length relationship of
#'       each species, or of the resource, see [resource_params()].}
#'     \item{`per_log_size`}{For an array that holds a density, whether to
#'       express it per logarithmic size (`TRUE`) rather than per size
#'       (`FALSE`). The default, `NULL`, leaves the density as it stands. An
#'       error for an array that does not hold a density.}
#'   }
#'
#'   **For `ArrayTimeBySpecies` methods:**
#'   \describe{
#'     \item{`tlim`}{A numeric vector of length two providing lower and upper
#'       limits for the time axis, e.g. `c(1980, 2000)`. Use `NA` to apply no
#'       limit at that end. Default is `c(NA, NA)`.}
#'   }
#'
#'   **For the `ArrayTimeBySpeciesBySize` and `ArrayTimeByResourceBySize`
#'   methods:**
#'   \describe{
#'     \item{`time`}{The time to display. Default (`NULL`) is the final time
#'       step.}
#'   }
#'
#' @return A ggplot2 object.
#'
#' @family plotting functions
#' @export
#' @examples
#' \donttest{
#' params <- NS_params
#' given_species_params(params)["Cod", "w_mat"] <- 1200
#' plotRelative(getEGrowth(NS_params), getEGrowth(params),
#'              wlim = c(500, 2000), log_x = FALSE, species = "Cod")
#'
#' # The same works for the resource
#' params2 <- setResource(NS_params,
#'                        resource_capacity = 2 * resource_capacity(NS_params))
#' plotRelative(resource_capacity(NS_params), resource_capacity(params2))
#' }
plotRelative <- function(x, y, species = NULL, log_x,
                         ylim = c(NA, NA), total = FALSE,
                         background = TRUE, highlight = NULL, ...) {
    UseMethod("plotRelative", x)
}

#' @rdname plotRelative
#' @usage NULL
#' @export
plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
                                            log_x = TRUE,
                                            ylim = c(NA, NA),
                                            total = FALSE,
                                            background = TRUE,
                                            highlight = NULL,
                                            all.sizes = FALSE,
                                            wlim = c(NA, NA),
                                            llim = c(NA, NA),
                                            size_axis = c("w", "l"),
                                            per_log_size = NULL, ...) {
    check_plot2_compatible(x, y, "ArraySpeciesBySize")
    compare_array_metadata(x, y)
    size_axis <- plot_size_axis(size_axis)
    check_per_log_size(x, per_log_size)
    assert_that(length(wlim) == 2,
                length(llim) == 2,
                length(ylim) == 2)
    params <- attr(x, "params")
    # As in `plot2()`, each array is prepared with its own model. The Jacobian
    # of a density no longer cancels out of the ratio when the two models
    # convert weight to length differently, and the two grids no longer
    # coincide, which is why the comparison is interpolated.
    plot_dat1 <- ArraySpeciesBySize_plot_data(
        x, species = species, all.sizes = all.sizes, wlim = wlim, llim = llim,
        total = total, background = background, size_axis = size_axis,
        per_log_size = per_log_size)
    plot_dat2 <- ArraySpeciesBySize_plot_data(
        y, species = species, all.sizes = all.sizes, wlim = wlim, llim = llim,
        total = total, background = background, size_axis = size_axis,
        per_log_size = per_log_size)

    plotRelativeDataFrame(plot_dat1, plot_dat2, params,
                          xlab = plot_size_xlab(size_axis),
                          xtrans = if (log_x) "log10" else "identity",
                          xlim = plot_size_xlim(wlim, size_axis, llim),
                          ylim = ylim, highlight = highlight,
                          legend_var = "Legend", interpolate = TRUE)
}

check_plot2_compatible <- function(x, y, class) {
    if (!inherits(y, class)) {
        stop("Both objects must be of class `", class, "`.")
    }
}

compare_array_metadata <- function(x, y) {
    value_name1 <- attr(x, "value_name")
    value_name2 <- attr(y, "value_name")
    if (!is.null(value_name1) && !is.null(value_name2) &&
            !identical(value_name1, value_name2)) {
        warning("The first array has value name `", value_name1,
                "`, but the second array has value name `", value_name2, "`.")
    }
    units1 <- attr(x, "units")
    units2 <- attr(y, "units")
    if (!is.null(units1) && !is.null(units2) &&
            nzchar(units1) && nzchar(units2) &&
            !identical(units1, units2)) {
        warning("The first array has y units `", units1,
                "`, but the second array has y units `", units2, "`.")
    }
    type1 <- array_type(x)
    type2 <- array_type(y)
    if (!identical(type1, type2)) {
        stop("The first array holds a value of type `", type1,
             "`, but the second array holds a value of type `", type2,
             "`. The type decides how the values are transformed and how the ",
             "axes are scaled, so two arrays of different types cannot be ",
             "drawn on one pair of axes.")
    }
}

#' Kinds of quantity a mizer array can hold
#'
#' Mizer arrays record what kind of quantity their values are in their `type`
#' attribute, because some kinds need handling that the numbers alone do not
#' reveal:
#' \describe{
#'   \item{`"value"`}{the default: a rate, an amount, anything that needs no
#'     special handling.}
#'   \item{`"density"`}{an amount per gram of body weight, like a number
#'     density. Plotting a density against a length axis restates it per
#'     centimetre, which changes the values and not just the axis.}
#'   \item{`"proportion"`}{a fraction, like the feeding level. Plotted on a
#'     linear y axis showing the whole of the interval from 0 to 1, so that the
#'     value can be read against the scale it belongs to.}
#' }
#'
#' A `"proportion"` is not *restricted* to the interval from 0 to 1: the
#' critical feeding level and the resource level can both exceed 1, and their
#' plots show it. The type is a statement about what the number means, not a
#' bound that mizer enforces.
#'
#' @format A character vector of the three types.
#' @keywords internal
array_types <- c("value", "density", "proportion")

#' Validate the type of a mizer array
#'
#' @param type One of [array_types].
#' @return The validated type.
#' @keywords internal
validate_array_type <- function(type) {
    if (!is.character(type) || length(type) != 1 || is.na(type) ||
            !type %in% array_types) {
        stop("`type` must be one of ",
             paste0("\"", array_types, "\"", collapse = ", "), ".")
    }
    type
}

#' Resolve the type of a mizer array
#'
#' Called by the array constructors. An explicit `type` is validated and used as
#' given; `NULL` means the constructor was called without the argument, in which
#' case a density is recognised from the other metadata, the way mizer
#' recognised one before the `type` attribute existed. That keeps arrays built
#' by extension packages, and arrays saved by earlier versions, behaving as they
#' did.
#'
#' @param type The type supplied to the constructor, or `NULL`.
#' @param value_name The `value_name` of the array.
#' @param units The `units` of the array.
#' @return One of [array_types].
#' @keywords internal
resolve_array_type <- function(type, value_name = NULL, units = NULL) {
    if (!is.null(type)) {
        return(validate_array_type(type))
    }
    if (identical(value_name, "Number density") || identical(units, "1/g")) {
        return("density")
    }
    "value"
}

#' The type of a mizer array
#'
#' @param x A mizer array object.
#' @return One of [array_types].
#' @keywords internal
array_type <- function(x) {
    resolve_array_type(attr(x, "type"),
                       attr(x, "value_name"), attr(x, "units"))
}

#' The density measure of a mizer array
#'
#' The bridge from the array metadata into the density machinery of the plots.
#' Mizer arrays are indexed by the model's weight grid, so a stored density is
#' always a density with respect to weight; the other measures in
#' [density_measures] arise only for quantities that the spectrum plots compute
#' on the fly, such as a density per logarithmic weight.
#'
#' @param x A mizer array object.
#' @return `"w"` if the array holds a density, otherwise `NA_character_`.
#' @keywords internal
array_density_wrt <- function(x) {
    if (identical(array_type(x), "density")) "w" else NA_character_
}

array_units <- function(x, size_axis = "w", per_log_size = NULL) {
    density_wrt <- array_density_wrt(x)
    target <- density_target_measure(density_wrt, size_axis, per_log_size)
    convert_density_units(attr(x, "units"), density_wrt, target)
}

array_y_label <- function(x, default = "Value", size_axis = "w",
                          per_log_size = NULL) {
    value_name <- attr(x, "value_name") %||% default
    # A density per logarithmic size is a different quantity from the density
    # itself and has to say so, since its units no longer distinguish it.
    if (isTRUE(per_log_size) && !is.na(array_density_wrt(x))) {
        value_name <- paste0(value_name,
                             if (identical(plot_size_axis(size_axis), "l"))
                                 " in log length" else " in log weight")
    }
    label_with_units(value_name, array_units(x, size_axis, per_log_size))
}

#' Check that `per_log_size` applies to a mizer array
#'
#' Expressing values per logarithmic size only means anything for a density,
#' so asking for it on anything else is an argument error rather than something
#' to be quietly ignored — which is what `...` used to do with it.
#'
#' @param x A mizer array object.
#' @param per_log_size The `per_log_size` argument of the plot method.
#' @return `per_log_size`, invisibly, if it applies.
#' @keywords internal
check_per_log_size <- function(x, per_log_size) {
    if (!is.null(per_log_size)) {
        assert_that(is.flag(per_log_size), noNA(per_log_size))
        if (is.na(array_density_wrt(x))) {
            stop("`per_log_size` only applies to an array that holds a ",
                 "density, but this one holds a value of type `",
                 array_type(x), "`.")
        }
    }
    invisible(per_log_size)
}

#' Add lines to an existing plot
#'
#' `r lifecycle::badge("experimental")`
#' `addPlot()` adds another set of values to an existing ggplot, for example to
#' compare the same rate before and after a model change. There are methods for
#' all the mizer array classes. Each checks whether the existing plot uses a
#' compatible x variable, and warns if the y variable or y-axis units appear to
#' differ.
#'
#' @param plot A ggplot2 object to which the new values should be added.
#' @param x An object containing the values to add. Can be an
#'   `ArraySpeciesBySize`, `ArrayTimeBySpecies`, `ArrayTimeBySpeciesBySize`,
#'   `ArrayResourceBySize` or `ArrayTimeByResourceBySize` object.
#' @param species Character vector of species to include. `NULL` (default) means
#'   all species. A resource array holds a single spectrum, so this argument is
#'   not used by the resource methods, which warn if it is set.
#' @param total A boolean value that determines whether the total is plotted
#'   as well. The total is the total of everything the array holds, every
#'   species and every size, whatever is drawn. Default is `FALSE`. Not used by
#'   the resource methods, which warn if it is set.
#' @param background A boolean value that determines whether background species
#'   are included. Ignored if the model does not contain background species.
#'   Default is `TRUE`. Not used by the resource methods, which warn if it is
#'   set.
#' @param colour Optional fixed colour for the added lines. If `NULL`, the
#'   species colours from the existing plot are used. Because a resource array
#'   is a single line whose "Resource" level may be missing from the existing
#'   plot's colour scale, the resource methods instead default to the fixed
#'   resource colour from `getColours()`.
#' @param linetype Optional fixed line type for the added lines. If `NULL`, the
#'   species line types from the existing plot are used.
#' @param linewidth Width of the added lines.
#' @param alpha Transparency of the added lines.
#' @param ... Further arguments used by only some of the methods:
#'
#'   **For the `ArraySpeciesBySize`, `ArrayTimeBySpeciesBySize`,
#'   `ArrayResourceBySize` and `ArrayTimeByResourceBySize` methods:**
#'   \describe{
#'     \item{`wlim`}{A numeric vector of length two providing lower and upper
#'       limits for the weight (x) axis. Use `NA` to refer to the existing
#'       minimum or maximum.}
#'   }
#'
#'   **For the `ArraySpeciesBySize` and `ArrayTimeBySpeciesBySize` methods:**
#'   \describe{
#'     \item{`all.sizes`}{If `FALSE` (default), values outside a species' size
#'       range (`w_min` to `w_max`) are removed.}
#'     \item{`llim`}{A numeric vector of length two providing lower and upper
#'       limits for the length (x) axis when `size_axis = "l"`. Use `NA` to
#'       refer to the existing minimum or maximum.}
#'     \item{`size_axis`}{Whether to plot size as weight (`"w"`, default) or
#'       length (`"l"`), using the allometric weight-length relationship of
#'       each species, or of the resource, see [resource_params()].}
#'     \item{`per_log_size`}{For an array that holds a density, whether to plot
#'       it per logarithmic size (`TRUE`) rather than per size (`FALSE`). The
#'       default, `NULL`, plots the density as it stands. Unlike `size_axis`
#'       this needs no weight-length relationship, so it is available for the
#'       resource classes too. An error for an array that does not hold a
#'       density.}
#'   }
#'
#'   **For `ArrayTimeBySpecies` methods:**
#'   \describe{
#'     \item{`tlim`}{A numeric vector of length two providing lower and upper
#'       limits for the time axis, e.g. `c(1980, 2000)`. Use `NA` to apply no
#'       limit at that end. Default is `c(NA, NA)`.}
#'     \item{`ylim`}{A numeric vector of length two providing lower and upper
#'       limits for the value (y) axis.}
#'   }
#'
#'   **For the `ArrayTimeBySpeciesBySize` and `ArrayTimeByResourceBySize`
#'   methods:**
#'   \describe{
#'     \item{`time`}{The time to display. Default (`NULL`) is the final time
#'       step.}
#'   }
#'
#' @return A ggplot2 object.
#' @export
#' @family plotting functions
#'
#' @examples
#' \donttest{
#' p <- plot(getEncounter(NS_params), species = "Cod")
#' addPlot(p, getEncounter(NS_params), species = "Cod")
#'
#' pr <- plot(getResourceMort(NS_params))
#' addPlot(pr, getResourceMort(NS_params))
#' }
addPlot <- function(plot, x, species = NULL, total = FALSE,
                    background = TRUE, colour = NULL, linetype = "dashed",
                    linewidth = 0.8, alpha = 1, ...) {
    UseMethod("addPlot", x)
}

#' @rdname addPlot
#' @usage NULL
#' @export
addPlot.ArraySpeciesBySize <- function(plot, x, species = NULL,
                                       total = FALSE,
                                       background = TRUE,
                                       colour = NULL,
                                       linetype = "dashed",
                                       linewidth = 0.8,
                                       alpha = 1,
                                       all.sizes = FALSE,
                                       wlim = c(NA, NA),
                                       llim = c(NA, NA),
                                       size_axis = c("w", "l"),
                                       per_log_size = NULL,
                                       ...) {
    if (!inherits(plot, "ggplot")) {
        stop("The `plot` argument must be a ggplot object.")
    }
    assert_that(is.number(linewidth),
                is.number(alpha),
                alpha >= 0,
                alpha <= 1)
    size_axis <- plot_size_axis(size_axis)
    check_per_log_size(x, per_log_size)
    assert_that(length(wlim) == 2,
                length(llim) == 2)

    plot <- deep_copy(plot)
    plot_dat <- ArraySpeciesBySize_plot_data(
        x, species = species, all.sizes = all.sizes, wlim = wlim, llim = llim,
        total = total, background = background, size_axis = size_axis,
        per_log_size = per_log_size)
    x_var <- plot_size_x_var(size_axis)
    y_var <- names(plot_dat)[2]
    check_addPlot_compatible(plot, x_var = x_var, y_var = y_var,
                             units = array_units(x, size_axis, per_log_size))

    mapping <- aes(x = .data[[x_var]], y = .data[[y_var]],
                   group = .data[["Species"]])
    if (is.null(colour)) {
        mapping$colour <- rlang::quo(.data[["Legend"]])
    }
    if (is.null(linetype)) {
        mapping$linetype <- rlang::quo(.data[["Legend"]])
    }

    layer_args <- list(
        data = plot_dat,
        mapping = mapping,
        linewidth = linewidth,
        alpha = alpha,
        inherit.aes = FALSE
    )
    if (!is.null(colour)) {
        layer_args$colour <- colour
    }
    if (!is.null(linetype)) {
        layer_args$linetype <- linetype
    }

    plot + do.call(geom_line, layer_args)
}

deep_copy <- function(x) {
    unserialize(serialize(x, NULL))
}

check_addPlot_compatible <- function(plot, x_var, y_var, units = NULL) {
    mapping <- plot$mapping
    if (length(plot$layers) > 0) {
        layer_mapping <- plot$layers[[1]]$mapping
        if (!is.null(layer_mapping$x)) {
            mapping$x <- layer_mapping$x
        }
        if (!is.null(layer_mapping$y)) {
            mapping$y <- layer_mapping$y
        }
    }

    plot_x_var <- plot_mapping_var(mapping$x)
    if (!is.null(plot_x_var) && !identical(plot_x_var, x_var)) {
        stop("The data can only be added to a plot with x variable `", x_var,
             "`. The existing plot uses x variable `", plot_x_var, "`.")
    }

    plot_y_var <- plot_mapping_var(mapping$y)
    if (!is.null(plot_y_var) && !identical(plot_y_var, y_var)) {
        warning("The existing plot appears to use y variable `", plot_y_var,
                "`, but the added data uses `", y_var, "`.")
    }

    plot_units <- plot_y_units(plot)
    if (!is.null(plot_units) && !is.null(units) &&
            nzchar(plot_units) && nzchar(units) &&
            !identical(plot_units, units)) {
        warning("The existing plot appears to use y units `", plot_units,
                "`, but the added data uses `", units, "`.")
    }
}

plot_mapping_var <- function(mapping) {
    if (is.null(mapping)) {
        return(NULL)
    }

    label <- rlang::as_label(mapping)
    match <- regmatches(label, regexec("\\.data\\[\\[\"([^\"]+)\"\\]\\]", label))[[1]]
    if (length(match) == 2) {
        return(match[[2]])
    }
    match <- regmatches(label, regexec("\\.data\\$([^ ]+)$", label))[[1]]
    if (length(match) == 2) {
        return(match[[2]])
    }
    if (grepl("^[[:alnum:]_.]+$", label)) {
        return(label)
    }

    NULL
}

plot_y_units <- function(plot) {
    scales <- plot$scales$scales
    for (scale in scales) {
        if ("y" %in% scale$aesthetics && is.character(scale$name)) {
            return(label_units(scale$name))
        }
    }
    label_units(plot$labels$y)
}

label_units <- function(label) {
    if (is.null(label) || !is.character(label) || !nzchar(label)) {
        return(NULL)
    }
    match <- regmatches(label, regexec("\\[([^]]+)\\]\\s*$", label))[[1]]
    if (length(match) == 2) {
        return(match[[2]])
    }
    NULL
}

#' Restrict plot data to a range of weights
#'
#' Internal helper that filters a plot data frame to the weight range given by
#' `wlim`. It is exported so that extension packages (such as mizerMR) can reuse
#' it in their own array `plot()` methods.
#'
#' @param data A data frame with a numeric `w` column.
#' @param wlim A length-2 numeric vector giving the lower and upper weight
#'   limits. Either entry may be `NA` to leave that side unrestricted.
#'
#' @return The subset of `data` with `w` inside `wlim`.
#' @keywords internal
#' @export
apply_wlim <- function(data, wlim) {
    if (!is.na(wlim[1])) data <- data[data$w >= wlim[1], ]
    if (!is.na(wlim[2])) data <- data[data$w <= wlim[2], ]
    data
}

#' The complete plotting data of a species-by-size array
#'
#' Everything a plot of an `ArraySpeciesBySize` needs, prepared once: the
#' species selection, the background grouping, the masking of sizes outside a
#' species' own range, the weight limits, the conversion of the values and of
#' the size coordinate onto the requested axis, the total line, and the length
#' limits.
#'
#' All of it uses the array's *own* `params`. That matters for the comparison
#' plots, where the two operands may come from different models: a length axis
#' and a density Jacobian are both built from the weight-length relationship of
#' the model the values came from, so preparing the second array with the first
#' one's parameters would put it in the wrong place on the axis. Each operand is
#' therefore prepared here, on its own, and the comparison renderers receive
#' data that is already on the axis it will be drawn against.
#'
#' @param x An `ArraySpeciesBySize` object.
#' @param species Character vector of species to include, or `NULL` for all.
#' @param all.sizes If `FALSE`, values outside a species' size range are
#'   removed.
#' @param wlim Numeric vector of length two giving the weight limits.
#' @param llim Numeric vector of length two giving the length limits, applied
#'   only on a length axis.
#' @param total Whether to append the total line, see [total_contributors()].
#' @param background Whether background species are included.
#' @param size_axis Either `"w"` (weight) or `"l"` (length).
#' @param per_log_size Whether to express a density per logarithmic size.
#' @return A data frame with the size coordinate in its first column, the values
#'   in its second, and `Species` and `Legend` columns.
#' @keywords internal
ArraySpeciesBySize_plot_data <- function(x, species = NULL,
                                         all.sizes = FALSE,
                                         wlim = c(NA, NA), llim = c(NA, NA),
                                         total = FALSE, background = TRUE,
                                         size_axis = "w",
                                         per_log_size = NULL) {
    params <- attr(x, "params")
    size_axis <- plot_size_axis(size_axis)
    plot_dat <- prepare_ArraySpeciesBySize_plot_data(
        x, species = species, all.sizes = all.sizes, wlim = wlim,
        total = total, background = background)
    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)
    }
    plot_dat
}

prepare_ArraySpeciesBySize_plot_data <- function(x, species = NULL,
                                                 all.sizes = FALSE,
                                                 wlim = c(NA, NA),
                                                 total = FALSE,
                                                 background = TRUE) {
    params <- attr(x, "params")
    w <- get_ArraySpeciesBySize_w(x)

    all_species <- rownames(x)
    if (is.null(species)) {
        species <- all_species
    } else {
        species <- intersect(species, all_species)
        if (length(species) == 0) {
            stop("None of the selected species are in the rate array.")
        }
    }

    value_name <- attr(x, "value_name") %||% "value"
    sel <- all_species %in% species
    mat <- unclass(x)[sel, , drop = FALSE]

    plot_dat <- data.frame(
        w = rep(w, each = sum(sel)),
        value = c(mat),
        Species = rownames(mat)
    )

    if (!all.sizes && !is.null(params)) {
        sp_params <- params@species_params
        for (sp in species) {
            if (sp %in% sp_params$species) {
                sp_row <- sp_params[sp_params$species == sp, ]
                plot_dat$value[plot_dat$Species == sp &
                                   (plot_dat$w < sp_row$w_min[1] |
                                        plot_dat$w > sp_row$w_max[1])] <- NA
            }
        }
        plot_dat <- plot_dat[complete.cases(plot_dat), ]
    }

    plot_dat <- apply_wlim(plot_dat, wlim)

    plot_dat$Legend <- plot_dat$Species

    # Handle background species
    if (!is.null(params) && isTRUE(any(params@species_params$is_background))) {
        bkgrd_sp <- params@species_params$species[params@species_params$is_background]
        if (background) {
            plot_dat$Legend[plot_dat$Species %in% bkgrd_sp] <- "Background"
        } else {
            plot_dat <- plot_dat[!plot_dat$Species %in% bkgrd_sp, ]
        }
    }

    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
#' \donttest{
#' plotHover(getEncounter(NS_params))
#' }
#' @export
plotHover.ArraySpeciesBySize <- function(x, ...) {
    plotHover(plot(x, ...), ...)
}

#' @export
as.data.frame.ArraySpeciesBySize <- function(x, row.names = NULL,
                                     optional = FALSE, ...) {
    w <- get_ArraySpeciesBySize_w(x)
    sp_names <- rownames(x)
    mat <- unclass(x)
    data.frame(
        w = rep(w, each = nrow(mat)),
        value = c(mat),
        Species = sp_names,
        row.names = row.names,
        check.names = !optional,
        stringsAsFactors = FALSE
    )
}

#' Get the size grid for an ArraySpeciesBySize object
#'
#' Internal helper that returns the consumer size grid `params@w` or the full
#' prey/resource size grid `params@w_full`, depending on the number of columns
#' in the array.
#'
#' @param x An `ArraySpeciesBySize` object.
#'
#' @return A numeric vector giving the size represented by each column. When the
#'   array is tagged as a bin average (`representation = "average"`) *and* the
#'   model uses second-order bin-averaging (`second_order_w[["bin_average"]]`),
#'   the geometric bin centres are returned instead of the left bin edges, so
#'   that bin-averaged quantities are drawn at the size where they actually live
#'   (see [bin_midpoints()]). Point-valued quantities and first-order models are
#'   unaffected, keeping default plots unchanged.
#' @keywords internal
get_ArraySpeciesBySize_w <- function(x) {
    params <- attr(x, "params")
    if (is.null(params)) {
        w <- as.numeric(colnames(x))
        if (any(is.na(w))) {
            w <- seq_len(ncol(x))
        }
        return(w)
    }
    average <- identical(attr(x, "representation"), "average") &&
        isTRUE(params@second_order_w[["bin_average"]])
    if (ncol(x) == length(params@w)) {
        return(if (average) bin_midpoints(params) else params@w)
    }
    if (ncol(x) == length(params@w_full)) {
        return(if (average) bin_midpoints(params, w_full = TRUE) else params@w_full)
    }
    stop("Can not determine the size grid for this ArraySpeciesBySize object. ",
         "The number of columns is ", ncol(x), ", but the params object has ",
         length(params@w), " consumer sizes and ", length(params@w_full),
         " full-spectrum sizes.")
}

#' @export
`[.ArraySpeciesBySize` <- function(x, i, j, ..., drop = TRUE) {
    result <- NextMethod()
    # Preserve class only if result is still a 2D matrix
    if (is.matrix(result) && length(dim(result)) == 2) {
        attr(result, "value_name") <- attr(x, "value_name")
        attr(result, "units") <- attr(x, "units")
        attr(result, "type") <- attr(x, "type")
        attr(result, "params") <- attr(x, "params")
        attr(result, "representation") <- attr(x, "representation")
        class(result) <- c("ArraySpeciesBySize", "matrix", "array")
    }
    result
}

#' @export
Ops.ArraySpeciesBySize <- function(e1, e2) {
    # Strip ArraySpeciesBySize class so that arithmetic returns a plain matrix.
    # We unclass both operands and call the generic directly.
    if (is.ArraySpeciesBySize(e1)) e1 <- unclass_rate(e1)
    if (!missing(e2) && is.ArraySpeciesBySize(e2)) e2 <- unclass_rate(e2)
    op <- match.fun(.Generic)
    if (missing(e2)) op(e1) else op(e1, e2)
}

# Helper to strip all ArraySpeciesBySize attributes
unclass_rate <- function(x) {
    x <- unclass(x)
    attr(x, "value_name") <- NULL
    attr(x, "units") <- NULL
    attr(x, "type") <- NULL
    attr(x, "params") <- NULL
    attr(x, "representation") <- NULL
    x
}

# Strip the `params` back-reference (a whole MizerParams) before calling the
# default str(), otherwise it would dump the entire model. We restore a plain
# array class for a normal summary, relabel line 1 with the real class, and
# append a one-line summary of the params attribute.
#' @export
str.ArraySpeciesBySize <- function(object, ...) {
    params <- attr(object, "params")
    attr(object, "params") <- NULL
    class(object) <- c("matrix", "array")
    out <- utils::capture.output(utils::str(object, ...))
    out[1] <- paste0(" 'ArraySpeciesBySize' ", sub("^ ", "", out[1]))
    cat(paste0(out, collapse = "\n"), "\n", sep = "")
    if (!is.null(params)) {
        cat(" - attr(*, \"params\")=Formal class 'MizerParams' [package \"mizer\"] with ",
            length(slotNames(params)), " slots\n", sep = "")
    }
    invisible(NULL)
}