R/ArrayTimeBySpecies-class.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

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

#' S3 class for time x species arrays
#'
#' Some functions in mizer return two-dimensional arrays (time x species)
#' holding quantities like biomass, abundance, or yield rate through time.
#' The `ArrayTimeBySpecies` class wraps these arrays to provide convenient
#' `print()`, `summary()`, `plot()`, and `as.data.frame()` methods.
#'
#' An `ArrayTimeBySpecies` object behaves just like a regular matrix for
#' arithmetic operations and subsetting. It carries these lightweight attributes:
#' \itemize{
#'   \item `value_name` – a human-readable name for the value
#'       (e.g. "Biomass").
#'   \item `units` – the units of the value (e.g. "g", "g/year").
#'   \item `params` – the `MizerParams` object that created the values.
#' }
#'
#' @param x A matrix (time x species). For `is.ArrayTimeBySpecies()`, 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", "g/year").
#' @param type The kind of quantity the values are, see [ArraySpeciesBySize()]
#'   and [array_types].
#' @param params A `MizerParams` object holding the model that created the
#'   values.
#'
#' @return An `ArrayTimeBySpecies` object (inherits from `matrix` and `array`).
#'
#' @seealso [print()], [summary()], [as.data.frame()], [plot()]
#' @export
#' @examples
#' \donttest{
#' bio <- getBiomass(NS_sim)
#' is.ArrayTimeBySpecies(bio)
#' summary(bio)
#' }
ArrayTimeBySpecies <- function(x, value_name = NULL, units = NULL,
                               type = NULL, params = NULL) {
    if (!is.matrix(x)) {
        stop("`x` must be a matrix.")
    }
    type <- resolve_array_type(type, value_name, units)
    structure(x,
        class = c("ArrayTimeBySpecies", "matrix", "array"),
        value_name = value_name,
        units = units,
        type = type,
        params = params
    )
}

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

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

    mat <- unclass_time(x)
    n_time <- nrow(mat)
    n_species <- ncol(mat)
    times <- parse_numeric_labels(rownames(mat), n_time)

    sp_idx <- pick_head_indices(n_species, mizer_print_defaults$species_head,
                                mizer_print_defaults$species_threshold)
    time_idx <- pick_log_spaced_indices(n_time, mizer_print_defaults$time_max,
                                        mizer_print_defaults$time_threshold)

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

    if (length(time_idx) < n_time) {
        cat(format_truncation_note(length(time_idx), n_time, "times",
                                   format_time_range_detail(times)), "\n")
    }
    if (length(sp_idx) < n_species) {
        omitted <- setdiff(colnames(mat), colnames(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")
    }
    invisible(x)
}

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

    df <- data.frame(
        Species = sp_names,
        Min = apply(mat, 2, min, na.rm = TRUE),
        Mean = apply(mat, 2, mean, na.rm = TRUE),
        Max = apply(mat, 2, 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.ArrayTimeBySpecies"
    result
}

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

#' Plot method for `ArrayTimeBySpecies` objects
#'
#' See [plot()] for an overview of the mizer plotting system and the
#' arguments shared by all of its methods.
#'
#' @param x An `ArrayTimeBySpecies` object.
#' @param species Character vector of species to include. `NULL`
#'   (default) means all species.
#' @param 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)`.
#' @param y_ticks The approximate number of ticks desired on the y axis.
#' @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 over every species the array holds,
#'   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 highlight Name or vector of names of the species to be
#'   highlighted.
#' @param log_x If `TRUE`, use a log10 x-axis. Default is `FALSE`.
#' @param log_y If `TRUE`, use a log10 y-axis. Default is `TRUE`.
#' @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 return_data If `TRUE`, return the data frame instead of the
#'   plot.
#' @param ... Unused.
#'
#' @return A ggplot2 object, unless `return_data = TRUE`, in which case a
#'   data frame is returned.
#' @keywords internal
#' @export
#' @examples
#' \donttest{
#' plot(getBiomass(NS_sim))
#' plot(getBiomass(NS_sim), species = c("Cod", "Herring"), total = TRUE)
#' plot(getYield(NS_sim), species = c("Cod", "Herring"))
#' }
plot.ArrayTimeBySpecies <- function(x, species = NULL,
                                    tlim = c(NA, NA),
                                    y_ticks = 6, ylim = c(NA, NA),
                                    total = FALSE, background = TRUE,
                                    highlight = NULL, log_x = FALSE,
                                    log_y = TRUE, log = NULL,
                                    return_data = FALSE,
                                    ...) {
    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

    value_name <- attr(x, "value_name") %||% "Value"
    units_str <- attr(x, "units")
    params <- attr(x, "params")

    # Construct y_label
    y_label <- value_name
    if (!is.null(units_str) && nzchar(units_str)) {
        y_label <- paste0(value_name, " [", units_str, "]")
    }
    plot_dat <- prepare_ArrayTimeBySpecies_plot_data(
        x, species = species, tlim = tlim,
        ylim = ylim, total = total, background = background, log_y = log_y)

    if (return_data) return(plot_dat)

    ylim <- array_ylim(x, ylim, log_y, plot_dat[[2]])
    plotDataFrame(plot_dat, params, xlab = "Year", ylab = y_label,
                  xtrans = if (log_x) "log10" else "identity",
                  ytrans = if (log_y) "log10" else "identity",
                  ylim = ylim, y_ticks = y_ticks, highlight = highlight,
                  legend_var = "Legend")
}

#' @rdname addPlot
#' @usage NULL
#' @export
addPlot.ArrayTimeBySpecies <- function(plot, x, species = NULL,
                                       total = FALSE,
                                       background = TRUE,
                                       colour = NULL,
                                       linetype = "dashed",
                                       linewidth = 0.8,
                                       alpha = 1,
                                       tlim = c(NA, NA),
                                       ylim = c(NA, NA),
                                       ...) {
    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)

    plot <- deep_copy(plot)
    # Filter the added values the way the plot they are joining was filtered.
    plot_dat <- prepare_ArrayTimeBySpecies_plot_data(
        x, species = species, tlim = tlim,
        ylim = ylim, total = total, background = background,
        log_y = plot_y_is_log(plot))
    y_var <- names(plot_dat)[[2]]
    check_addPlot_compatible(plot, x_var = "Year", y_var = y_var,
                             units = attr(x, "units"))

    mapping <- aes(x = .data[["Year"]], 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)
}

#' @rdname plot2
#' @usage NULL
#' @export
plot2.ArrayTimeBySpecies <- function(x, y, name1 = "First", name2 = "Second",
                                     species = NULL,
                                     log_x = FALSE, log_y = TRUE,
                                     log = NULL, ylim = c(NA, NA),
                                     total = FALSE, background = TRUE,
                                     highlight = NULL,
                                     y_ticks = 6,
                                     tlim = c(NA, NA), ...) {
    check_plot2_compatible(x, y, "ArrayTimeBySpecies")
    compare_array_metadata(x, y)
    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

    params <- attr(x, "params")
    y_label <- array_y_label(x, default = "Value")
    plot_dat1 <- prepare_ArrayTimeBySpecies_plot_data(
        x, species = species, tlim = tlim,
        ylim = ylim, total = total, background = background, log_y = log_y)
    plot_dat2 <- prepare_ArrayTimeBySpecies_plot_data(
        y, species = species, tlim = tlim,
        ylim = ylim, total = total, background = background, log_y = log_y)

    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 = "Year", ylab = y_label,
                            xtrans = if (log_x) "log10" else "identity",
                            ytrans = if (log_y) "log10" else "identity",
                            ylim = ylim, y_ticks = y_ticks,
                            highlight = highlight, legend_var = "Legend")
}

#' @rdname plotRelative
#' @usage NULL
#' @export
plotRelative.ArrayTimeBySpecies <- function(x, y, species = NULL,
                                            log_x = FALSE,
                                            ylim = c(NA, NA),
                                            total = FALSE,
                                            background = TRUE,
                                            highlight = NULL,
                                            tlim = c(NA, NA), ...) {
    check_plot2_compatible(x, y, "ArrayTimeBySpecies")
    compare_array_metadata(x, y)
    params <- attr(x, "params")
    # The relative difference is drawn on a linear axis, so a value of zero or
    # less is kept: it is a real difference between the two models rather than
    # something a logarithm cannot show.
    plot_dat1 <- prepare_ArrayTimeBySpecies_plot_data(
        x, species = species, tlim = tlim,
        total = total, background = background, log_y = FALSE)
    plot_dat2 <- prepare_ArrayTimeBySpecies_plot_data(
        y, species = species, tlim = tlim,
        total = total, background = background, log_y = FALSE)

    plotRelativeDataFrame(plot_dat1, plot_dat2, params,
                          xlab = "Year",
                          xtrans = if (log_x) "log10" else "identity",
                          ylim = ylim, highlight = highlight,
                          legend_var = "Legend")
}

#' The complete plotting data of a time-by-species array
#'
#' The species selection and the background grouping are done in a single pass,
#' so that no species can be both selected under its own name and appended again
#' under the `"Background"` legend — which is what a separate appending step
#' used to do to every background species whenever `species` was left at its
#' default of all of them.
#'
#' @param x An `ArrayTimeBySpecies` object.
#' @param species Character vector of species to include, or `NULL` for all.
#' @param tlim Numeric vector of length two giving the time limits.
#' @param ylim Numeric vector of length two giving the value limits. Values
#'   outside them are dropped, as they cannot be seen anyway.
#' @param total Whether to append the total, which is the total over every
#'   species the array holds, whatever is drawn.
#' @param background Whether background species are included.
#' @param log_y Whether the values will be drawn on a logarithmic axis. Only
#'   then are non-positive values dropped: they have no place on a log axis, but
#'   on a linear one they are data like any other, and a quantity that can go
#'   negative — a rate of change, a difference between two models — would
#'   otherwise lose exactly the part of it that is interesting.
#' @return A data frame with `Year`, value, `Species` and `Legend` columns.
#' @keywords internal
prepare_ArrayTimeBySpecies_plot_data <- function(x, species = NULL,
                                                 tlim = c(NA, NA),
                                                 ylim = c(NA, NA),
                                                 total = FALSE,
                                                 background = TRUE,
                                                 log_y = TRUE) {
    value_name <- attr(x, "value_name") %||% "Value"
    params <- attr(x, "params")

    # Filter to time range
    t <- as.numeric(rownames(x))
    if (any(is.na(t))) t <- seq_len(nrow(x))
    if (!is.na(tlim[1]) && !is.na(tlim[2]) && tlim[1] >= tlim[2]) {
        stop("tlim[1] must be less than tlim[2]")
    }
    if (!is.na(tlim[1])) {
        x <- x[t >= tlim[1], , drop = FALSE]
        t <- t[t >= tlim[1]]
    }
    if (!is.na(tlim[2])) {
        x <- x[t <= tlim[2], , drop = FALSE]
    }

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

    bm <- unclass(x)

    # The total is the total over every species the array holds, so it is formed
    # before the selection narrows the columns.
    if (total) {
        bm <- cbind(bm, Total = rowSums(bm))
    }

    bm <- reshape2::melt(bm)

    # Implement ylim; bring columns in desired order
    keep <- (is.na(ylim[1]) | bm$value >= ylim[1]) &
        (is.na(ylim[2]) | bm$value <= ylim[2])
    if (isTRUE(log_y)) {
        # A hard floor rather than a plain `> 0`, so that a species that has
        # died out does not drag a logarithmic axis down by twenty decades.
        keep <- keep & bm$value >= 1e-20
    }
    bm <- bm[keep, c(1, 3, 2)]
    names(bm) <- c("Year", value_name, "Species")

    bkgrd_sp <- character(0)
    if (!is.null(params) && isTRUE(any(params@species_params$is_background))) {
        bkgrd_sp <- params@species_params$species[params@species_params$is_background]
    }

    # One pass over the rows decides both what is drawn and what it is called.
    # A background species is drawn when the selection asks for it, and is
    # always labelled "Background", the way the species-by-size plots label it.
    selected <- as.character(bm$Species) %in% c("Total", species)
    if (!background) {
        selected <- selected & !(as.character(bm$Species) %in% bkgrd_sp)
    }
    plot_dat <- bm[selected, , drop = FALSE]
    plot_dat$Legend <- as.character(plot_dat$Species)
    plot_dat$Legend[plot_dat$Legend %in% bkgrd_sp] <- "Background"

    plot_dat
}

#' @rdname plotHover
#' @examples
#' \donttest{
#' plotHover(getBiomass(NS_sim))
#' }
#' @export
plotHover.ArrayTimeBySpecies <- function(x, ...) {
    plotHover(plot(x, ...), ...)
}

#' @export
as.data.frame.ArrayTimeBySpecies <- function(x, row.names = NULL,
                                             optional = FALSE, ...) {
    t <- as.numeric(rownames(x))
    if (any(is.na(t))) {
        t <- seq_len(nrow(x))
    }
    sp_names <- colnames(x)
    mat <- unclass(x)
    data.frame(
        time = rep(t, times = ncol(mat)),
        value = c(mat),
        Species = rep(sp_names, each = nrow(mat)),
        row.names = row.names,
        check.names = !optional,
        stringsAsFactors = FALSE
    )
}

#' @export
`[.ArrayTimeBySpecies` <- 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")
        class(result) <- c("ArrayTimeBySpecies", "matrix", "array")
    }
    result
}

#' @export
Ops.ArrayTimeBySpecies <- function(e1, e2) {
    # Strip ArrayTimeBySpecies class so that arithmetic returns a plain matrix.
    if (is.ArrayTimeBySpecies(e1)) e1 <- unclass_time(e1)
    if (!missing(e2) && is.ArrayTimeBySpecies(e2)) e2 <- unclass_time(e2)
    op <- match.fun(.Generic)
    if (missing(e2)) op(e1) else op(e1, e2)
}

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

# Strip the `params` back-reference so the default str() doesn't dump the whole
# MizerParams; summarise it in a single line instead. See str.ArraySpeciesBySize.
#' @export
str.ArrayTimeBySpecies <- 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(" 'ArrayTimeBySpecies' ", 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)
}