plotDiet.MizerSim ignores time-varying abundance; no time_range parameter

[gustavdelius]

Description

plotDiet.MizerSim currently discards the simulation's time-varying abundance and delegates directly to plotDiet(object@params, ...):

plotDiet.MizerSim <- function(object, species = NULL, ...) {
    plotDiet(object@params, species = species, ...)
}

This is inconsistent with every other .MizerSim plot method (plotFeedingLevel, plotPredMort, plotFMort, plotSpectra, plotCDF, …), which all accept a time_range argument and compute the plotted quantity from the simulation's actual abundance at the specified time.

The diet composition can change over the course of a simulation as species abundances shift, so ignoring the simulated n may produce misleading results when the system is not at the initial equilibrium.

Suggested fix

Add a time_range parameter to plotDiet.MizerSim. Extract the abundance at (or averaged over) the specified time step and pass it to plot_diet instead of params@initial_n. If no time is supplied, default to the final time step (matching the convention in plotFeedingLevel and plotPredMort).

[gustavdelius]

Suggested Fix for getDiet and plotDiet on MizerSim

As suggested, adding support for time-varying abundance requires introducing a getDiet.MizerSim S3 method and modifying plotDiet.MizerSim to use it.

Below are the proposed implementation details:

1. S3 Method and Generic updates for getDiet in R/summary_methods.R

First, the generic getDiet function must be extended with ... to allow methods to accept additional arguments like time_range and drop:

getDiet <- function(params,
                    n = initialN(params),
                    n_pp = initialNResource(params),
                    n_other = initialNOther(params),
                    proportion = TRUE,
                    ...) {
    UseMethod("getDiet")
}

# Update the MizerParams method to accept ...
getDiet.MizerParams <- function(params,
                    n = initialN(params),
                    n_pp = initialNResource(params),
                    n_other = initialNOther(params),
                    proportion = TRUE,
                    ...) { ... }

Next, add the S3 method for MizerSim to compute the diet over time steps/slices:

#' @rdname getDiet
#' @usage NULL
#' @export
getDiet.MizerSim <- function(params,
                             n = initialN(params),
                             n_pp = initialNResource(params),
                             n_other = initialNOther(params),
                             proportion = TRUE,
                             time_range,
                             drop = FALSE,
                             ...) {
    sim <- params
    time_elements <- get_sim_rate_time_elements(sim, time_range)
    
    diet_time <- plyr::aaply(which(time_elements), 1, function(time_idx) {
        slice <- get_sim_rate_slice(sim, time_idx)
        getDiet(sim@params, n = slice$n, n_pp = slice$n_pp,
                n_other = slice$n_other, proportion = proportion, ...)
    }, .drop = FALSE)
    names(dimnames(diet_time))[[1]] <- "time"
    
    result <- diet_time[, , , , drop = drop]
    return(result)
}

2. Update plotDiet.MizerSim in R/plots.R

Update plotDiet.MizerSim to accept time_range, average the abundance/diet over that range, and pass them to the internal plot_diet function:

#' @rdname plotDiet
#' @usage NULL
#' @export
plotDiet.MizerSim <- function(object, species = NULL,
                              wlim = c(NA, NA), llim = c(NA, NA),
                              size_axis = c("w", "l"),
                              return_data = FALSE,
                              log_x = TRUE, log_y = FALSE, log = NULL,
                              time_range, ...) {
    size_axis <- plot_size_axis(size_axis)
    log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
    
    if (missing(time_range)) {
        time_range <- max(as.numeric(dimnames(object@n)$time))
    }
    
    params <- validParams(object@params)
    diet <- getDiet(object, time_range = time_range, drop = FALSE, ...)
    if (length(dim(diet)) == 4) {
        diet <- apply(diet, c(2, 3, 4), mean)
    }
    
    time_elements <- get_time_elements(object, time_range)
    n <- apply(object@n[time_elements, , , drop = FALSE], c(2, 3), mean)
    
    plot_diet(params, n = n, diet = diet, species = species,
              log_x = log_axes$log_x, log_y = log_axes$log_y,
              wlim = wlim, llim = llim, size_axis = size_axis,
              return_data = return_data)
}

3. Namespace Registration

Register the S3 method in the NAMESPACE file:

S3method(getDiet,MizerSim)