R/project.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

# Project function for the size based modelling package mizer

# Copyright 2012 Finlay Scott and Julia Blanchard.
# Copyright 2018 Gustav Delius and Richard Southwell.
# Distributed under the GPL 3 or later
# Maintainer: Gustav Delius, University of York, <gustav.delius@york.ac.uk>

#' @useDynLib mizer
#' @importFrom Rcpp sourceCpp
NULL


#' Project size spectrum forward in time
#'
#' Runs the size spectrum model simulation.
#' The function returns an object of type
#' \linkS4class{MizerSim} that can then be explored with a range of
#' [summary_functions], [indicator_functions] and
#' [plotting_functions].
#'
#' @param object Either a \linkS4class{MizerParams} object or a
#'   \linkS4class{MizerSim} object (which contains a `MizerParams` object).
#' @param effort The effort of each fishing gear through time. See notes below.
#' @param t_max The number of years the projection runs for. The default value
#'   is 100. When an effort array is supplied, this argument can be used to
#'   extend the simulation beyond the times specified in the effort array.
#'   See notes below.
#' @param dt Time step of the solver. The default value is 0.1. When `object`
#'   is a `MizerSim`, defaults to the value used to produce that simulation.
#' @param t_save The frequency with which the output is stored. The default
#'   value is 1. See notes below.
#' @param t_start The the year of the start of the simulation. The simulation
#'   will cover the period from `t_start` to \code{t_start + t_max}.
#'   Defaults to 0. Ignored if an array is used for the `effort`
#'   argument or a `MizerSim` for the `object` argument.
#' @param initial_n `r lifecycle::badge("deprecated")` The initial abundances of
#'   species. Instead of using this argument you should set `initialN(params)`
#'   to the desired value.
#' @param initial_n_pp `r lifecycle::badge("deprecated")` The initial abundances
#'   of resource. Instead of using this argument you should set
#'   `initialNResource(params)` to the desired value.
#' @param append A boolean that determines whether the new simulation results
#'   are appended to the previous ones. Only relevant if `object` is a
#'   `MizerSim` object. Default = TRUE.
#' @param progress_bar Either a boolean value to determine whether a progress
#'   bar should be shown in the console, or a shiny Progress object to implement
#'   a progress bar in a shiny app.
#' @param callback A function to be called at each saved time step of the
#'   simulation. The callback function is called with the `MizerSim` object,
#'   the current time index, and any additional arguments passed to `project()`
#'   via `...`.
#' @param check_steady `r lifecycle::badge("experimental")`
#'   If `TRUE`, warn when the model is not at its steady state before the
#'   projection starts, which catches the common mistake of forgetting to re-run
#'   [tuneSteadyState()] after a `match…`/`calibrate…` step. Default `FALSE`,
#'   because
#'   projecting a model away from its steady state is a perfectly normal thing
#'   to do. It is meant for a `MizerParams` object; when continuing from a
#'   `MizerSim` the starting state is deliberately wherever the previous run
#'   ended, so there is nothing to check. The check is made at the effort stored in the params
#'   object rather than at the `effort` supplied here, so that running a fishing
#'   scenario at a new effort — which legitimately starts away from the steady
#'   state of that new effort — does not warn. See [getSteadyResidual()].
#' @param method The numerical method to use for the consumer density update.
#'
#'   `"euler"` uses the first-order semi-implicit Euler update.
#'
#'   `"tr_bdf2"` uses the L-stable, second-order TR-BDF2 method and is the
#'   recommended second-order method. `"second_order"` is accepted as an
#'   easier-to-remember alias for it, matching the name of the spatial
#'   second-order option [second_order_w()].
#'
#'   `"predictor_corrector"` is superseded by `"tr_bdf2"`. It is also second
#'   order and costs the same, but its Crank-Nicolson corrector is only
#'   A-stable, so it rings at large time steps where TR-BDF2 does not. It is
#'   retained for backwards compatibility and for comparison.
#'   `"predictor-corrector"` is accepted as an alias for it.
#'
#'   The two second-order methods treat the nonlinear rates identically. The
#'   Crank-Nicolson update in `"predictor_corrector"` is only the *corrector*,
#'   applied to the transport operator with those rates frozen, so the method
#'   does not inherit Crank-Nicolson's stability properties for the full
#'   nonlinear dynamics. The same holds for the L-stability of `"tr_bdf2"`; see
#'   the section "Custom rates must depend continuously on abundance" below.
#'
#'   Aliases are resolved to the canonical name, so [getSimParams()] reports
#'   `"tr_bdf2"` even when `"second_order"` was supplied.
#'
#'   When `object` is a `MizerSim`, defaults to the
#'   value used to produce that simulation. A warning is issued if
#'   `append = TRUE` and the supplied value differs from the stored one.
#' @param ... Other arguments will be passed to rate functions.
#'
#' @section Advective flux scheme:
#' The spatial discretisation of the growth (advection) term is controlled by
#' the `flux` entry of the `second_order_w` slot of the params object, not by
#' an argument to `project()`. With the default (`"upwind"`) the first-order
#' upwind flux is used. Setting it to `"van_leer"` (for example with
#' `second_order_w(params) <- TRUE`) switches on a flux-limited (van Leer, TVD)
#' deferred correction that removes the leading numerical diffusion
#' \eqn{\approx g\,w\,\log\beta} of the upwind flux while keeping the density
#' update a tridiagonal solve and preserving positivity. The correction is most
#' useful on coarse logarithmic grids and pairs naturally with the second-order
#' time methods. Because it changes the discrete steady state, the choice lives
#' in the params object alongside the steady state rather than being a per-run
#' argument. See [second_order_w()].
#'
#' @section Custom rates must depend continuously on abundance:
#' All three methods are *semi-implicit*: the densities are solved for
#' implicitly, but the rates that build the transport operator are frozen at
#' values computed from earlier states. The second-order methods gain their
#' extra order by evaluating the rates twice — at the start of the step and from
#' a provisional prediction of its end — and averaging. That average is only
#' second order if the rates vary smoothly along the trajectory.
#'
#' A custom rate function registered with [setRateFunction()] that depends
#' discontinuously on the abundances therefore defeats all three methods alike,
#' including the L-stable `"tr_bdf2"`, whose damping applies to the frozen
#' linear operator and not to the rates. The result is a trajectory that keeps
#' changing as `dt` is refined. See the [Discontinuous rate
#' functions](https://sizespectrum.org/mizer/articles/discontinuous_rates.html)
#' article for the symptoms and the remedy.
#'
#' @note The `effort` argument specifies the level of fishing effort during the
#'   simulation. If it is not supplied, the initial effort stored in the params
#'   object is used. The effort can be specified in four different ways:
#' \itemize{
#' \item A single numeric value. This specifies the effort of all fishing gears
#' which is constant through time (i.e. all the gears have the same constant
#' effort).
#' \item A named vector whose names match with existing gear names.
#'  The values in the vector specify the constant fishing effort for those
#'  fishing gears, i.e. the effort is constant through time. The
#'   effort for gears that are not included in the effort vector is set to the
#'   default effort value, which is 1 in defaults edition 2 and later and 0 in
#'   earlier defaults editions. Missing (`NA`) effort entries are replaced in
#'   the same way.
#' \item A numerical vector which has the same length as the number of fishing
#' gears. The values in the vector specify the
#' constant fishing effort of each of the fishing gears, with the ordering
#' assumed to be the same as in the MizerParams object.
#' \item A numerical array with dimensions time x gear. This specifies the
#' fishing effort of each gear at each time step.  The first dimension, time,
#' must be named numerically and increasing. The second dimension of the array
#' must be named and the names must correspond to the gear names in the
#' `MizerParams` object. The value for the effort for a particular time
#' is used during the interval from that time to the next time in the array.
#' }
#'
#' If effort is specified as an array then the smallest time in the array is
#' used as the initial time for the simulation. Otherwise the initial time is
#' set to the final time of the previous simulation if `object` is a
#' `MizerSim` object or to `t_start` otherwise.
#'
#' When an effort array is provided, the `t_max` argument can be used to
#' extend the simulation beyond the last time specified in the effort array.
#' In this case, the effort values from the last time in the array will be
#' used for the extended period. The `t_save` argument can be used to specify
#' the frequency at which simulation results are saved. If `t_save` is not
#' supplied, the results will be saved at the times specified in the effort
#' array. If both `t_max` and `t_save` are provided with an effort array,
#' effort values will be interpolated (using step function) or extrapolated
#' (using the last known value) as needed for the new time points. The
#' `t_start` argument continues to be ignored when an effort array is supplied.
#'
#' Note that if `t_max` or `t_save` are specified, the time grid for the
#' simulation is resampled based on `t_save`. This means that if the time
#' points in the `effort` array are irregular and do not align with the new
#' grid, those specific time points may be lost and the effort values at the
#' new grid points will be calculated via interpolation.
#'
#' If the `object` argument is of class `MizerSim` then the initial
#' values for the simulation are taken from the final values in the
#' `MizerSim` object and the corresponding arguments to this function will
#' be ignored.
#'
#' @return An object of class \linkS4class{MizerSim}.
#'
#' @export
#' @examples
#' \donttest{
#' params <- NS_params
#' # With constant fishing effort for all gears for 20 time steps
#' sim <- project(params, t_max = 20, effort = 0.5)
#' # With constant fishing effort which is different for each gear
#' effort <- c(Industrial = 0, Pelagic = 1, Beam = 0.5, Otter = 0.5)
#' sim <- project(params, t_max = 20, effort = effort)
#' # With fishing effort that varies through time for each gear
#' gear_names <- c("Industrial", "Pelagic", "Beam", "Otter")
#' times <- seq(from = 1, to = 10, by = 1)
#' effort_array <- array(NA,
#'     dim = c(length(times), length(gear_names)),
#'     dimnames = list(time = times, gear = gear_names)
#' )
#' effort_array[, "Industrial"] <- 0.5
#' effort_array[, "Pelagic"] <- seq(from = 1, to = 2, length = length(times))
#' effort_array[, "Beam"] <- seq(from = 1, to = 0, length = length(times))
#' effort_array[, "Otter"] <- seq(from = 1, to = 0.5, length = length(times))
#' sim <- project(params, effort = effort_array)
#' # Extend a simulation beyond the effort array times
#' # Effort values from the final time are used for the extension
#' sim <- project(params, effort = effort_array, t_max = 15)
#' # Control save times with an effort array using t_save
#' sim <- project(params, effort = effort_array, t_save = 2)
#' }
#' @rdname project
#' @export
project <- function(object, effort,
                    t_max = 100, dt = 0.1, t_save = 1, t_start = 0,
                    initial_n, initial_n_pp,
                    append = TRUE,
                    progress_bar = TRUE,
                    callback = NULL,
                    method = c("euler", "predictor_corrector", "tr_bdf2"),
                    check_steady = FALSE,
                    ...) {
    UseMethod("project")
}

# Aliases accepted by the `method` argument of project() and friends, mapped to
# the canonical name that is stored on the MizerSim and used for dispatch.
# `"second_order"` is the memorable name for the recommended second-order
# method; `"predictor-corrector"` is the hyphenated spelling.
project_method_aliases <- c("predictor-corrector" = "predictor_corrector",
                            "second_order"        = "tr_bdf2")

normalise_project_method <- function(method) {
    if (length(method) != 1) {
        method <- method[[1]]
    }
    method <- match.arg(method, c("euler", "predictor_corrector", "tr_bdf2",
                                  names(project_method_aliases)))
    if (method %in% names(project_method_aliases)) {
        method <- project_method_aliases[[method]]
    }
    method
}

#' @rdname project
#' @usage NULL
#' @export
project.MizerParams <- function(object, effort,
                                t_max = 100, dt = 0.1, t_save = 1, t_start = 0,
                                initial_n, initial_n_pp,
                                append = TRUE,
                                progress_bar = TRUE,
                                callback = NULL,
                                method = c("euler", "predictor_corrector", "tr_bdf2"),
                                check_steady = FALSE,
                                ...) {
    params <- validParams(object)
    method <- normalise_project_method(method)
    assert_that(is.flag(check_steady))
    if (check_steady) {
        # Judged at the model's own effort, not at the `effort` supplied here:
        # see the argument documentation.
        warn_if_not_steady(
            params,
            paste("Projecting from a state that is not a fixed point means",
                  "the results include that drift."))
    }
    # The advective-flux scheme is a property of the model, held in the
    # second_order_w slot, not a per-run choice.
    flux_limiter <- flux_limiter_scheme(params)
    # Set and check initial values ----
    assert_that(t_max > 0)
    if (!missing(initial_n)) params@initial_n[] <- initial_n
    if (!missing(initial_n_pp)) params@initial_n_pp[] <- initial_n_pp

    initial_n <- params@initial_n
    initial_n_pp <- params@initial_n_pp
    initial_n_other <- params@initial_n_other

    no_sp <- length(params@w_min_idx)
    assert_that(
        is.array(initial_n),
        is.numeric(initial_n),
        are_equal(dim(initial_n), c(no_sp, length(params@w)))
    )
    assert_that(
        is.numeric(initial_n_pp),
        length(initial_n_pp) == length(params@w_full)
    )

    assert_that(is.null(initial_n_other) || is.list(initial_n_other))
    other_names <- names(params@other_dynamics)
    if (length(other_names) > 0) {
        if (is.null(names(initial_n_other))) {
            stop("The initial_n_other needs to be a named list")
        }
        if (!setequal(names(initial_n_other), other_names)) {
            stop(
                "The names of the entries in initial_n_other do not match ",
                "the names of the other components of the model."
            )
        }
    }

    # Set effort array ----
    if (missing(effort)) effort <- params@initial_effort

    if (is.null(dim(effort))) { # effort is a vector or scalar
        # Set up the effort array transposed so we can use the recycling rules
        # no point running a simulation with no saved results
        if (t_max < t_save) {
            t_save <- t_max
        }
        times <- seq(t_start, t_start + t_max, by = t_save)
        t_end <- t_start + t_max
        if (!isTRUE(all.equal(times[length(times)], t_end))) {
            warning("`t_max` (", t_max, ") is not a multiple of `t_save` (",
                    t_save, "). Results will be saved at `t_max` but the ",
                    "final save interval will be shorter than `t_save`.",
                    call. = FALSE)
            times <- c(times, t_end)
        }
        effort <- validEffortVector(effort, params)
        effort <- t(array(effort,
            dim = c(length(effort), length(times)),
            dimnames = list(
                gear = names(effort),
                time = times
            )
        ))
    } else {
        effort <- validEffortArray(effort, params)

        # When an effort array is supplied, we now respect t_max and t_save
        # if they are explicitly provided by the user
        effort_times <- as.numeric(dimnames(effort)[[1]])
        t_start_effort <- effort_times[1]
        t_end_effort <- effort_times[length(effort_times)]

        # Check if t_max or t_save were explicitly provided by the user
        t_max_provided <- !missing(t_max)
        t_save_provided <- !missing(t_save)

        if (t_max_provided || t_save_provided) {
            # Determine the end time
            if (t_max_provided) {
                t_end <- t_start_effort + t_max
            } else {
                # Keep the end time from the effort array
                t_end <- t_end_effort
            }

            # Determine the save frequency
            if (t_save_provided) {
                # Use the provided t_save
                save_freq <- t_save
            } else {
                # Infer from effort array times or use default
                if (length(effort_times) > 1) {
                    # Use the first interval as the save frequency
                    save_freq <- effort_times[2] - effort_times[1]
                } else {
                    save_freq <- 1 # default
                }
            }

            # Ensure we don't use a t_save larger than the duration
            duration <- t_end - t_start_effort
            if (duration < save_freq) {
                save_freq <- duration
            }

            # Generate new time points on the requested save grid.
            times <- seq(t_start_effort, t_end, by = save_freq)

            # Interpolate effort values for the new time points
            gear_names <- dimnames(effort)[[2]]
            new_effort <- array(NA,
                dim = c(length(times), length(gear_names)),
                dimnames = list(time = times, gear = gear_names)
            )

            for (g in seq_along(gear_names)) {
                # Use approx with rule=2 for constant extrapolation beyond range
                # method="constant" means the effort value applies from that time
                # until the next time point (step function)
                new_effort[, g] <- stats::approx(
                    x = effort_times,
                    y = effort[, g],
                    xout = times,
                    method = "constant",
                    rule = 2,
                    f = 0
                )$y
            }
            effort <- new_effort
        }
    }

    times <- as.numeric(dimnames(effort)[[1]])

    # Make the MizerSim object with the right size ----
    # We only save every t_save years
    sim <- MizerSim(params, t_dimnames = times)
    sim@sim_params <- list(method = method, dt = dt,
                           flux_limiter = flux_limiter)
    # Set initial population and effort
    sim@n[1, , ] <- initial_n
    sim@n_pp[1, ] <- initial_n_pp
    sim@n_other[1, ] <- initial_n_other
    sim@effort <- effort

    ## Initialise ----
    # get functions
    resource_dynamics_fn <- get(params@resource_dynamics)
    other_dynamics_fns <- lapply(params@other_dynamics, get)
    rates_fns <- projectRateFunctions(params)

    # Set up progress bar
    if (is(progress_bar, "Progress")) {
        # We have been passed a shiny progress object
        progress_bar$set(message = "Running simulation", value = 0)
        proginc <- 1 / length(times)
    } else if (progress_bar == TRUE) {
        pb <- progress::progress_bar$new(
            format = "[:bar] :percent ETA: :eta",
            total = length(times), width = 60
        )
        pb$tick(0)
    }

    # Call callback on initial state
    if (!is.null(callback)) {
        callback(sim, t_idx = 1, ...)
    }

    n_list <- list(
        n = initial_n, n_pp = initial_n_pp,
        n_other = unserialize(serialize(initial_n_other, NULL))
    )
    t <- times[[1]]

    ## Loop over time ----
    for (i in 2:length(times)) {
        # number of time steps between saved times
        steps <- round((times[[i]] - t) / dt)
        # advance to next saved time
        n_list <- project_simple(
            params,
            n = n_list$n, n_pp = n_list$n_pp, n_other = n_list$n_other,
            t = t, dt = dt, steps = steps,
            effort = effort[i - 1, ],
            resource_dynamics_fn = resource_dynamics_fn,
            other_dynamics_fns = other_dynamics_fns,
            rates_fns = rates_fns,
            method = method, flux_limiter = flux_limiter, ...
        )
        # Calculate start time for next iteration
        # The reason we don't simply use the next entry in `times` is that
        # those entries may not be separated by exact multiples of dt.
        t <- t + steps * dt
        # Advance progress bar
        if (is(progress_bar, "Progress")) {
            progress_bar$inc(amount = proginc)
        } else if (progress_bar == TRUE) {
            pb$tick()
        }

        # Store result
        sim@n[i, , ] <- n_list$n
        sim@n_pp[i, ] <- n_list$n_pp
        sim@n_other[i, ] <- unserialize(serialize(n_list$n_other, NULL))

        if (!is.null(callback)) {
            callback(sim, t_idx = i, ...)
        }
    }

    return(sim)
}

#' @rdname project
#' @usage NULL
#' @export
project.MizerSim <- function(object, effort,
                             t_max = 100, dt = 0.1, t_save = 1, t_start = 0,
                             initial_n, initial_n_pp,
                             append = TRUE,
                             progress_bar = TRUE,
                             callback = NULL,
                             method = c("euler", "predictor_corrector", "tr_bdf2"),
                             ...) {
    validObject(object)
    stored <- object@sim_params
    dt_provided <- !missing(dt)
    method_provided <- !missing(method)
    # Default dt and method to values used in the existing simulation. The
    # advective-flux scheme is a property of the params, so it is carried over
    # automatically.
    if (!dt_provided && !is.null(stored$dt)) {
        dt <- stored$dt
    }
    if (!method_provided && !is.null(stored$method)) {
        method <- stored$method
    }
    method <- normalise_project_method(method)
    # Warn when appending with different numerical parameters
    if (append && length(stored) > 0) {
        if (!is.null(stored$dt) && !isTRUE(all.equal(dt, stored$dt))) {
            warning("Appending a simulation run with dt = ", dt,
                    " to one that used dt = ", stored$dt, ".")
        }
        if (!is.null(stored$method) &&
                method != normalise_project_method(stored$method)) {
            warning("Appending a simulation run with method = '", method,
                    "' to one that used method = '", stored$method, "'.")
        }
    }
    params <- finalParams(object)
    t_start <- getTimes(object)[idxFinalT(object)]

    sim <- project(params,
        effort = effort, t_max = t_max, dt = dt,
        t_save = t_save, t_start = t_start,
        initial_n = initial_n, initial_n_pp = initial_n_pp,
        progress_bar = progress_bar,
        callback = callback,
        method = method, ...
    )

    if (append) {
        times <- as.numeric(dimnames(sim@n)[[1]])
        no_t_old <- dim(object@n)[1]
        no_t <- length(times)
        new_t_dimnames <- c(
            as.numeric(dimnames(object@n)[[1]]),
            times[2:no_t]
        )
        new_sim <- MizerSim(params, t_dimnames = new_t_dimnames)
        old_indices <- 1:no_t_old
        new_indices <- seq(from = no_t_old + 1, length.out = no_t - 1)
        new_sim@n[old_indices, , ] <- object@n
        new_sim@n[new_indices, , ] <- sim@n[2:no_t, , ]
        new_sim@n_pp[old_indices, ] <- object@n_pp
        new_sim@n_pp[new_indices, ] <- sim@n_pp[2:no_t, ]
        new_sim@n_other[old_indices, ] <- object@n_other
        new_sim@n_other[new_indices, ] <- sim@n_other[2:no_t, ]
        new_sim@effort[old_indices, ] <- object@effort
        new_sim@effort[new_indices, ] <- sim@effort[2:no_t, ]
        new_sim@sim_params <- sim@sim_params
        return(new_sim)
    }
    return(sim)
}

#' Project abundances by a given number of time steps into the future
#'
#' This is an internal function used by the user-facing `project()` function.
#' It is of potential interest only to mizer extension authors.
#'
#' The function does not check its arguments because it is meant to be as fast
#' as possible to allow it to be used in a loop. For example, it is called in
#' `project()` once for every saved value. The function also does not save its
#' intermediate results but only returns the result at time `t + dt * steps`.
#' During this time it uses the constant fishing effort `effort`.
#'
#' The functional arguments can be calculated from slots in the `params` object
#' with
#' ```
#' resource_dynamics_fn <- get(params@resource_dynamics)
#' other_dynamics_fns <- lapply(params@other_dynamics, get)
#' rates_fns <- lapply(params@rates_funcs, get)
#' ```
#' The reason the function does not do that itself is to shave 20 microseconds
#' of its running time, which pays when the function is called hundreds of
#' times in a row.
#'
#' This function is also used by the steady-state finders. In between calls to
#' `project_simple()` the `steady()` function checks whether the values are
#' still changing significantly, so that it can stop when a steady state has
#' been approached. Mizer extension packages might have a similar need to run
#' a simulation repeatedly for short periods to run some other code in
#' between. Because this code may want to use the values of the rates from the
#' final update step, these too are included in the returned list.
#'
#' @param params A MizerParams object.
#' @param n An array (species x size) with the number density at start of
#'   simulation.
#' @param n_pp A vector (size) with the resource number density at start of
#'   simulation.
#' @param n_other A named list with the abundances of other components at start
#'   of simulation.
#' @param t Time at the start of the simulation.
#' @param dt Size of time step.
#' @param steps The number of time steps by which to project.
#' @param effort The fishing effort to be used throughout the simulation. This
#'   must be a vector or list with one named entry per fishing gear.
#' @param resource_dynamics_fn The function for the resource
#'   dynamics. See Details.
#' @param other_dynamics_fns List with the functions for the
#'   dynamics of the other components. See Details.
#' @param rates_fns List with the functions for calculating
#'   the rates. See Details.
#' @param method The numerical method to use for the consumer density update.
#'   See [project()].
#' @param ... Other arguments that are passed on to the rate functions.
#' @return List with the final values of `n`, `n_pp`, and `n_other`, together
#'   with `rates`, the rates calculated at the start of the final update step.
#'
#' @export
#' @concept helper
project_simple <- function(params, n, n_pp, n_other, effort, t, dt, steps,
                           resource_dynamics_fn, other_dynamics_fns,
                           rates_fns,
                           method = c("euler", "predictor_corrector", "tr_bdf2"),
                           ...) {
    UseMethod("project_simple")
}

#' @export
project_simple.MizerParams <-
    function(params,
             n = params@initial_n,
             n_pp = params@initial_n_pp,
             n_other = params@initial_n_other,
             effort = params@initial_effort,
             t = 0, dt = 0.1, steps,
             resource_dynamics_fn = get(params@resource_dynamics),
             other_dynamics_fns = lapply(params@other_dynamics, get),
             rates_fns = projectRateFunctions(params),
             method = c("euler", "predictor_corrector", "tr_bdf2"), ...) {
        method <- normalise_project_method(method)
        # Advective-flux scheme is read from the model's second_order_w slot.
        flux_limiter <- flux_limiter_scheme(params)
        # Handy things ----
        no_sp <- nrow(params@species_params) # number of species
        no_w <- length(params@w) # number of fish size bins
        idx <- 2:no_w
        # Hacky shortcut to access the correct element of a 2D array using 1D
        # notation
        # This references the egg size bracket for all species, so for example
        # n[w_min_idx_array_ref] = n[,w_min_idx]
        w_min_idx_array_ref <- (params@w_min_idx - 1) * no_sp + (1:no_sp)
        # Matrices for solver
        a <- matrix(0, nrow = no_sp, ncol = no_w)
        b <- matrix(0, nrow = no_sp, ncol = no_w)
        c <- matrix(0, nrow = no_sp, ncol = no_w)
        S <- matrix(0, nrow = no_sp, ncol = no_w)

        # Advance the other components by one Euler step from the start-of-step
        # state (`n`, `n_pp`, `n_other`) using the supplied rates. Used for the
        # predictor (with the start-of-step rates `r`) and, for the
        # second-order methods, the corrector (with the midpoint rates
        # `r_mid`), so that the other components are integrated to the same
        # order of accuracy as the resource and consumer spectra.
        step_other_components <- function(rates) {
            out <- list()
            for (component in names(params@other_dynamics)) {
                out[[component]] <-
                    other_dynamics_fns[[component]](
                        params,
                        n = n, n_pp = n_pp, n_other = n_other,
                        rates = rates, t = t, dt = dt,
                        component = component, ...
                    )
            }
            out
        }

        # Loop over time steps ----
        for (i_time in 1:steps) {
            r <- rates_fns$Rates(
                params,
                n = n, n_pp = n_pp, n_other = n_other,
                t = t, effort = effort, rates_fns = rates_fns, ...
            )

            # * Update other components (predictor) ----
            # First-order Euler update using the start-of-step rates `r`. For
            # the Euler method this is the final value for the step. The
            # second-order methods use `n_other_hat` in the end-of-step rates
            # (like the resource's `n_pp_hat`) and then overwrite `n_other_new`
            # with a corrector below.
            n_other_hat <- step_other_components(r)
            n_other_new <- n_other_hat

            # * Update resource and species ----
            if (method == "predictor_corrector" || method == "tr_bdf2") {
                # Second-order methods: compute midpoint rates once and use them
                # for both the resource and the consumer update.
                # Predictor: advance resource and consumer with start-of-step
                # rates r.
                n_pp_hat <- resource_dynamics_fn(params,
                    n = n, n_pp = n_pp,
                    n_other = n_other, rates = r,
                    t = t, dt = dt,
                    resource_rate = params@rr_pp,
                    resource_capacity = params@cc_pp, ...
                )
                n_hat <- project_n(params, r, n, dt, a, b, c, S, idx,
                                   w_min_idx_array_ref, no_sp, no_w,
                                   flux_limiter = flux_limiter)
                # End-of-step rates from the prediction, then midpoint rates
                r_hat <- rates_fns$Rates(
                    params,
                    n = n_hat, n_pp = n_pp_hat, n_other = n_other_hat,
                    t = t + dt, effort = effort, rates_fns = rates_fns, ...
                )
                r_mid <- average_rates(r, r_hat)

                # Corrector: other components with the midpoint rates, so they
                # reach the same second-order accuracy as the resource and
                # consumer. Computed from the start-of-step state before the
                # resource corrector below overwrites `n_pp`.
                n_other_new <- step_other_components(r_mid)

                # Corrector: resource with midpoint resource mortality
                n_pp <- resource_dynamics_fn(params,
                    n = n, n_pp = n_pp,
                    n_other = n_other, rates = r_mid,
                    t = t, dt = dt,
                    resource_rate = params@rr_pp,
                    resource_capacity = params@cc_pp, ...
                )
                # Corrector: consumer with the same midpoint rates. Supplying
                # r_mid makes the stepper skip its internal predictor.
                if (method == "predictor_corrector") {
                    n <- project_n_2(params, r, n, dt, a, b, c, S, idx,
                                     w_min_idx_array_ref, no_sp, no_w,
                                     r_mid = r_mid, n_hat = n_hat,
                                     flux_limiter = flux_limiter)
                } else {
                    n <- project_n_tr_bdf2(params, r, n, dt, a, b, c, S, idx,
                                           w_min_idx_array_ref, no_sp, no_w,
                                           r_mid = r_mid, n_hat = n_hat,
                                           flux_limiter = flux_limiter)
                }
            } else {
                # First-order Euler method.
                n_pp <- resource_dynamics_fn(params,
                    n = n, n_pp = n_pp,
                    n_other = n_other, rates = r,
                    t = t, dt = dt,
                    resource_rate = params@rr_pp,
                    resource_capacity = params@cc_pp, ...
                )
                if (any(r$diffusion > 0) || flux_limiter != "none") {
                    n <- project_n(params, r, n, dt, a, b, c, S, idx,
                                   w_min_idx_array_ref, no_sp, no_w,
                                   flux_limiter = flux_limiter)
                } else {
                    n <- project_n_no_diffusion(params, r, n, dt, a, b, S, idx,
                                                w_min_idx_array_ref, no_sp,
                                                no_w)
                }
            }

            # * Carry the updated other components into the next step ----
            # They were computed into `n_other_new` so that the resource and
            # consumer updates above still saw the start-of-step value. Now that
            # those updates are done, advance `n_other` so that the next time
            # step (and the returned value) uses the new abundances.
            n_other <- n_other_new

            # * Update time ----
            t <- t + dt
        }

        return(list(n = n, n_pp = n_pp, n_other = n_other, rates = r))
    }

validEffortArray <- function(effort, params) {
    # Check that number and names of gears in effort array is same as in
    # MizerParams object
    no_gears <- dim(params@catchability)[1]
    if (dim(effort)[2] != no_gears) {
        stop(
            "The number of gears in the effort array (length of the second dimension = ",
            dim(effort)[2],
            ") does not equal the number of gears in the MizerParams object (",
            no_gears, ")."
        )
    }
    gear_names <- dimnames(params@catchability)[[1]]
    if (!all(gear_names %in% dimnames(effort)[[2]])) {
        stop(
            "Gear names in the MizerParams object (",
            paste(gear_names, collapse = ", "),
            ") do not match those in the effort array."
        )
    }
    # Sort effort array to match order in MizerParams
    effort <- effort[, gear_names, drop = FALSE]

    if (is.null(dimnames(effort)[[1]])) {
        stop("The time dimname of the effort argument must be numeric.")
    }
    time_effort <- as.numeric(dimnames(effort)[[1]])
    if (any(is.na(time_effort))) {
        stop("The time dimname of the effort argument must be numeric.")
    }
    if (is.unsorted(time_effort)) {
        stop("The time dimname of the effort argument should be increasing.")
    }

    # Replace any NA's with default value
    effort_default <- ifelse(defaults_edition() < 2, 0, 1)
    effort[is.na(effort)] <- effort_default

    names(dimnames(effort)) <- c("time", "gear")
    effort
}