R/setResource.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

#' Set resource dynamics
#'
#' Sets the intrinsic resource birth rate and the intrinsic resource carrying
#' capacity as well as the name of the function used to simulate the resource
#' dynamics. By default, the birth rate and the carrying capacity are changed
#' together in such a way that the resource replenishes at the same rate at
#' which it is consumed. So you should only provide either the
#' `resource_rate` or the `resource_capacity` (or `resource_level`) because
#' the other is determined by the requirement that the resource replenishes
#' at the same rate at which it is consumed.
#'
#' You would usually set the resource dynamics only after having finished the
#' calibration of the steady state. Then setting the resource dynamics with
#' this function will preserve that steady state, unless you explicitly
#' choose to set `balance = FALSE`. Your choice of the resource dynamics only
#' affects the dynamics around the steady state. The higher the resource rate
#' or the lower the resource capacity the less sensitive the model will be to
#' changes in the competition for resource.
#'
#' If you provide the `resource_level` then that sets the `resource_capacity`
#' to the current resource number density divided by the resource level. So
#' in that case you should not specify `resource_capacity` as well.
#'
#' If you provide none of the arguments `resource_level`, `resource_rate` or
#' `resource_capacity`, and you do not change any of the resource parameters,
#' then the resource rate is kept at its previous value and, when balancing, the
#' capacity is recalculated from it. If instead you change one of the resource
#' parameters (`kappa`, `lambda`, `n` or `w_pp_cutoff`) or set `reset = TRUE`,
#' the rate and capacity are recalculated from the resource parameters (and then
#' balanced, unless `balance = FALSE`).
#'
#' @section Setting resource dynamics:
#'
#' The `resource_dynamics` argument allows you to choose the resource dynamics
#' function. By default, mizer uses a semichemostat model to describe the
#' resource dynamics in each size class independently. This semichemostat
#' dynamics is implemented by the function [resource_semichemostat()]. You can
#' change that to use a logistic model implemented by [resource_logistic()] or
#' you can use [resource_constant()] which keeps the resource constant or you
#' can write your own function.
#'
#' Both the [resource_semichemostat()] and the [resource_logistic()] dynamics
#' are parametrised in terms of a size-dependent birth rate \eqn{r_R(w)} and a
#' size-dependent capacity \eqn{c_R}. The help pages of these functions give
#' the details.
#'
#' The `resource_rate` argument can be a vector (with the same length as
#' `w_full(params)`) specifying the intrinsic resource birth rate for each size
#' class. Alternatively it can be a single number that is used as the
#' coefficient in a power law: then the intrinsic birth rate \eqn{r_R(w)} at
#' size \eqn{w} is set to
#' \deqn{r_R(w) = r_R w^{n-1}.}
#' The power-law exponent \eqn{n} is taken from the `n` argument.
#'
#' The `resource_capacity` argument can be a vector specifying the intrinsic
#' resource carrying capacity for each size class. Alternatively it can be a
#' single number that is used as the coefficient in a truncated power
#' law: then the intrinsic carrying capacity \eqn{c_R(w)} at size \eqn{w}
#' is set to
#' \deqn{c_R(w) = c_R\, w^{-\lambda}}{c_R(w) = c_R w^{-\lambda}}
#' for all \eqn{w} less than `w_pp_cutoff` and zero for larger sizes.
#' The power-law exponent \eqn{\lambda} is taken from the `lambda` argument.
#'
#' The values for `lambda`, `n` and `w_pp_cutoff` are stored in a list
#' in the `resource_params` slot of the MizerParams object so that they can be
#' re-used automatically in the future. If you specify `resource_rate` or
#' `resource_capacity` as a single number, that coefficient is likewise stored,
#' as `r_pp` and `kappa` respectively. That list can be accessed with
#' [resource_params()].
#'
#' The resource power law also determines defaults for species search volume.
#' Changing `lambda` recalculates any `q` and `gamma` values that mizer
#' calculated, and changing `kappa` (by supplying a scalar
#' `resource_capacity`) recalculates any calculated `gamma`. Species-specific
#' values that you supplied explicitly remain unchanged.
#'
#' @param params A MizerParams object
#' @param resource_rate Optional. A vector of per-capita resource birth
#'   rate for each size class or a single number giving the coefficient in the
#'   power-law for this rate, see "Setting resource dynamics" below.
#'   Must be strictly positive.
#' @param resource_capacity Optional. Vector of resource intrinsic carrying
#'   capacities or coefficient in the power-law for the capacity, see
#'   "Setting resource dynamics" below.
#'   The resource capacity must not be smaller than the resource abundance.
#' @param resource_level Optional. The ratio between the current resource number
#'   density and the resource capacity. Either a number used at all sizes or a
#'   vector specifying a value for each size. Must be greater than 0 and at
#'   most 1,
#'   except at sizes where the resource is zero, where it can be `NaN`. This
#'   determines the resource capacity, so do not specify both this and
#'   `resource_capacity`.
#' @param resource_dynamics Optional. Name of the function that determines the
#'   resource dynamics by calculating the resource spectrum at the next time
#'   step from the current state.
#' @param balance By default, if possible, the resource parameters are
#'   set so that the resource replenishes at the same rate at which it is
#'   consumed. In this case you should only specify either the resource rate
#'   or the resource capacity (or resource level) because the other is then
#'   determined automatically. Set to FALSE if you do not want the balancing.
#' @param n Used to set power-law exponent for resource rate if the
#'   `resource_rate` argument is given as a single number.
#' @param lambda Used to set power-law exponent for resource capacity if the
#'   `resource_capacity` argument is given as a single number.
#' @param w_pp_cutoff The upper cut off size of the resource spectrum power law
#'   used when `resource_capacity` is given as a single number. When changing
#'   `w_pp_cutoff` without providing `resource_capacity`, the cutoff can only be
#'   decreased. In that case, both the carrying capacity and the initial
#'   resource abundance will be cut off at the new value. To increase the
#'   cutoff, you must also provide the `resource_capacity` for the extended
#'   range.
#' @param reset
#'   If set to TRUE, then the resource capacity and birth rate will be reset
#'   to the values calculated from the resource parameters, even if they were
#'   previously overwritten with custom values. If set to FALSE (default) then a
#'   recalculation from the resource parameters will take place only if no custom
#'   values have been set.
#' @param ... Unused
#'
#' @return `setResource`: A MizerParams object with updated resource parameters
#' @seealso [setParams()]
#' @export
setResource <- function(params, resource_rate = NULL, resource_capacity = NULL,
                        resource_level = NULL, resource_dynamics = NULL,
                        lambda = resource_params(params)[["lambda"]],
                        n = resource_params(params)[["n"]],
                        w_pp_cutoff = resource_params(params)[["w_pp_cutoff"]],
                        balance = NULL, reset = FALSE, ...) {
    UseMethod("setResource")
}
#' @export
setResource.MizerParams <- function(params,
                        resource_rate = NULL,
                        resource_capacity = NULL,
                        resource_level = NULL,
                        resource_dynamics = NULL,
                        lambda = resource_params(params)[["lambda"]],
                        n = resource_params(params)[["n"]],
                        w_pp_cutoff = resource_params(params)[["w_pp_cutoff"]],
                        balance = NULL,
                        reset = FALSE,
                        ...) {

    assert_that(is.flag(reset))
    if (reset) {
        if (!is.null(resource_capacity) || !is.null(resource_rate) || !is.null(resource_level)) {
            warning("Because you set `reset = TRUE`, the values you provided for `resource_capacity`, `resource_rate`, or `resource_level` will be ignored and values will be calculated from the resource parameters.")
            resource_capacity <- NULL
            resource_rate <- NULL
            resource_level <- NULL
        }
        comment(params@cc_pp) <- NULL
        comment(params@rr_pp) <- NULL
    }

    resource_rate_user <- resource_rate
    resource_capacity_user <- resource_capacity
    if (!is.null(resource_level)) {
        resource_capacity_user <- resource_level
    }

    args <- list(...)
    # `setResource()` declares its `...` as unused, so without this check any
    # misspelled argument would be silently ignored. Only the two deprecated
    # names and the internal flag used by `resource_params<-()` are expected
    # here.
    unknown <- setdiff(names(args),
                       c("r_pp", "kappa", "resource_params_changed"))
    if (length(unknown) > 0) {
        stop("`setResource()` does not have ",
             if (length(unknown) == 1) "an argument " else "arguments ",
             paste0("`", unknown, "`", collapse = ", "), ".")
    }

    resource_param_changes <- args[["resource_params_changed"]]
    changed_resource_params <- if (is.character(resource_param_changes)) {
        resource_param_changes
    } else {
        character()
    }
    resource_params_changed <- isTRUE(resource_param_changes) ||
        length(changed_resource_params) > 0 ||
        !missing(lambda) || !missing(n) || !missing(w_pp_cutoff) || reset

    old_lambda <- params@resource_params[["lambda"]]
    old_kappa <- params@resource_params[["kappa"]]

    if ("r_pp" %in% names(args)) {
        lifecycle::deprecate_warn("1.0.0", "setResource(r_pp)",
                                  "setResource(resource_rate)")
        resource_rate <- args[["r_pp"]]
        resource_rate_user <- resource_rate
    }
    if ("kappa" %in% names(args)) {
        lifecycle::deprecate_warn("1.0.0", "setResource(kappa)",
                                  "setResource(resource_capacity)")
        resource_capacity <- args[["kappa"]]
        resource_capacity_user <- resource_capacity
    }
    assert_that(is.number(lambda),
                is.number(w_pp_cutoff), w_pp_cutoff > 0,
                is.number(n))

    # Store the old w_pp_cutoff before updating
    old_w_pp_cutoff <- params@resource_params[["w_pp_cutoff"]]

    params@resource_params[["lambda"]] <- lambda
    params@resource_params[["n"]] <- n
    params@resource_params[["w_pp_cutoff"]] <- w_pp_cutoff
    # The weight-length parameters feed no rate, so they are not arguments
    # here; they are filled in so that `resource_params()` shows them and the
    # user can change them there.
    params@resource_params[["a"]] <-
        params@resource_params[["a"]] %||% resource_length_defaults$a
    params@resource_params[["b"]] <-
        params@resource_params[["b"]] %||% resource_length_defaults$b

    if (!is.null(resource_capacity) && !is.null(resource_level)) {
        stop("You should specify only either 'resource_level' or 'resource_capacity'.")
    }

    # Check and set dynamics function ----
    if (!is.null(resource_dynamics)) {
        assert_that(is.character(resource_dynamics))
        if (!is.function(get0(resource_dynamics))) {
            stop('The resource dynamics function "', resource_dynamics, '" is not defined.')
        }
        params@resource_dynamics <- resource_dynamics
    }

    w_full <- w_full(params)
    no_w_full <- length(w_full)
    mu <- getResourceMort(params)
    NR <- initialNResource(params)

    # Check resource level ----
    if (!is.null(resource_level)) {
        assert_that(is.numeric(resource_level))
        if (length(resource_level) != 1 && length(resource_level) != no_w_full) {
            stop("The 'resource_level' should have length 1 or length ",
                 no_w_full, ".")
        }
        # The resource level is allowed to be NaN only where the resource is 0
        if (any(NR > 0 & is.nan(resource_level))) {
            stop("The resource level must be defined everywhere where the current resource is non-vanishing.")
        }
        if (any(NR > 0 &
                (resource_level <= 0 | resource_level > 1))) {
            stop("The 'resource_level' must always be greater than 0 and at most 1.")
        }
        resource_capacity <- NR / resource_level
        resource_capacity[is.nan(resource_level)] <- 0
        if (is.null(comment(resource_level))) {
            if (is.null(comment(params@cc_pp))) {
                comment(resource_capacity) <- "set manually"
            } else {
                comment(resource_capacity) <- comment(params@cc_pp)
            }
        } else {
            comment(resource_capacity) <- comment(resource_level)
        }
    }

    # Check growth rate ----
    if (!is.null(resource_rate)) {
        assert_that(is.numeric(resource_rate))
        if (length(resource_rate) == 1) {
            params@resource_params[["r_pp"]] <- resource_rate
            co <- comment(resource_rate)
            if (isTRUE(params@second_order_w[["bin_average"]])) {
                # Exact bin average of the power law r_pp * w^(n-1) over each
                # bin, so the relaxation rate is consistent with the
                # finite-volume cell-average resource density. See the
                # "Point values and bin averages" section of the
                # numerical-details vignette.
                resource_rate <- resource_rate *
                    power_law_bin_average(w_full, params@dw_full, n - 1)
            } else {
                resource_rate <- resource_rate * w_full ^ (n - 1)
            }
            comment(resource_rate) <- co
        } else if (length(resource_rate) != no_w_full) {
            stop("The 'resource_rate' should have length 1 or length ",
                 no_w_full, ".")
        } else {
            if (is.null(comment(resource_rate))) {
                if (is.null(comment(params@rr_pp))) {
                    comment(resource_rate) <- "set manually"
                } else {
                    comment(resource_rate) <- comment(params@rr_pp)
                }
            }
        }
        if (any(resource_rate < 0)) {
            stop("The 'resource_rate' must always be non-negative.")
        }
    }

    # Check capacity ----
    if (!is.null(resource_capacity)) {
        assert_that(is.numeric(resource_capacity))
        if (length(resource_capacity) == 1) {
            params@resource_params[["kappa"]] <- resource_capacity
            co <- comment(resource_capacity)
            if (isTRUE(params@second_order_w[["bin_average"]])) {
                # Exact bin average of kappa * w^(-lambda), truncated at the
                # cutoff. The unpredated semichemostat equilibrium is N_R* = c_p,
                # so the stored capacity must be the cell average of the
                # background spectrum for it to match the bin-averaged resource
                # consumed by the (bin-integrated) encounter convolution. The bin
                # straddling w_pp_cutoff gets the partial average; bins above it
                # are zero.
                resource_capacity <- resource_capacity *
                    power_law_bin_average(
                        w_full, params@dw_full, -lambda,
                        w_max = params@resource_params$w_pp_cutoff)
            } else {
                resource_capacity <- resource_capacity * w_full ^ (-lambda)
                resource_capacity[w_full >= params@resource_params$w_pp_cutoff] <- 0
            }
            comment(resource_capacity) <- co
        } else if (length(resource_capacity) != no_w_full) {
            stop("The 'resource_capacity' should have length 1 or length ",
                 no_w_full, ".")
        } else {
            if (is.null(comment(resource_capacity))) {
                if (is.null(comment(params@cc_pp))) {
                    comment(resource_capacity) <- "set manually"
                } else {
                    comment(resource_capacity) <- comment(params@cc_pp)
                }
            }
        }
        if (any(resource_capacity < 0)) {
            stop("The 'resource_capacity' must never be negative.")
        }
    }

    # Handle w_pp_cutoff increase error when capacity is not explicitly provided ----
    if (is.null(resource_capacity_user) &&
        !is.null(old_w_pp_cutoff) && w_pp_cutoff > old_w_pp_cutoff) {
        stop("You cannot increase w_pp_cutoff without also providing the resource_capacity for the extended range.")
    }

    # Recompute capacity from stored scalar kappa if not provided and not commented ----
    if (is.null(resource_capacity) && resource_params_changed) {
        capacity <- capacity_from_resource_params(params, lambda)
        if (!is.null(comment(params@cc_pp))) {
            # cc_pp is commented (frozen), so the resource parameters no longer
            # determine it. Say so, but only if a change was actually requested.
            if (!is.null(capacity) && different(capacity, params@cc_pp)) {
                signal_not_recalculated(
                    "cc_pp", "resource capacity",
                    "setResource(params, reset = TRUE)",
                    derived_from = "resource parameters")
            }
        } else if (!is.null(capacity)) {
            resource_capacity <- capacity
            if (!is.null(old_w_pp_cutoff) && w_pp_cutoff < old_w_pp_cutoff) {
                params@initial_n_pp[w_full >= w_pp_cutoff] <- 0
                NR <- params@initial_n_pp
            }
        }
    }
    # A frozen capacity still has to follow a `w_pp_cutoff` that came down,
    # because the resource cannot extend beyond the grid the cutoff defines.
    if (is.null(resource_capacity) && !is.null(comment(params@cc_pp)) &&
        !is.null(old_w_pp_cutoff) && w_pp_cutoff < old_w_pp_cutoff) {
        params@cc_pp[w_full >= w_pp_cutoff] <- 0
        params@initial_n_pp[w_full >= w_pp_cutoff] <- 0
        NR <- params@initial_n_pp
    }

    # Recompute rate from stored scalar r_pp if not provided and not commented ----
    # This is deliberately independent of whether the capacity was recomputed:
    # each array is rebuilt from its own scalars so that rate-side parameters
    # (`r_pp`, `n`) take effect even when capacity-side parameters also changed.
    if (is.null(resource_rate) && resource_params_changed) {
        rate <- rate_from_resource_params(params, n)
        if (!is.null(comment(params@rr_pp))) {
            if (!is.null(rate) && different(rate, params@rr_pp)) {
                signal_not_recalculated(
                    "rr_pp", "resource rate",
                    "setResource(params, reset = TRUE)",
                    derived_from = "resource parameters")
            }
        } else {
            resource_rate <- rate
        }
    }

    # `gamma` is calculated against the idealised resource spectrum and `q`
    # against its slope. Rebuild from the given species parameters, just as the
    # species-parameter setters do, so that mizer-owned values are recalculated
    # while explicitly given values remain protected. Do this before balancing
    # because balancing uses predation mortality, which depends on search
    # volume.
    resource_default_changes <- unique(c(
        intersect(changed_resource_params, c("kappa", "lambda")),
        if (!identical(old_kappa,
                       params@resource_params[["kappa"]])) "kappa",
        if (!identical(old_lambda,
                       params@resource_params[["lambda"]])) "lambda"
    ))
    if (length(resource_default_changes) > 0) {
        params <- rebuild_from_given(params, params@species_params)
    }

    # Balance ----
    balance_fn <- get0(paste0("balance_", params@resource_dynamics))
    if (is.null(balance)) {
        balance <- is.function(balance_fn)
    }
    if (balance) {
        num_args_user <- (!is.null(resource_rate_user)) +
            (!is.null(resource_capacity_user))
        if (num_args_user > 1) {
            stop("You should only provide either the `resource_rate` or `resource_capacity` (or `resource_level`) because the other is determined by the requirement that the resource replenishes at the same rate at which it is consumed.")
        }

        if (num_args_user == 1) {
            if (!is.null(resource_capacity_user)) {
                resource_rate <- NULL
            } else {
                resource_capacity <- NULL
            }
        } else {
            # num_args_user == 0
            if (!is.null(resource_capacity)) {
                resource_rate <- NULL
            } else if (!is.null(resource_rate)) {
                resource_capacity <- NULL
            } else {
                # Neither changed from scalar; use existing rr_pp to balance capacity
                resource_rate <- params@rr_pp
                resource_capacity <- NULL
            }
        }

        # Whichever side is NULL going into the balancing function is the one
        # that gets derived from the other. A frozen (commented) array is only
        # protected from *incidental* balancing, i.e. when the user did not
        # explicitly supply the complementary rate/capacity/level
        # (`num_args_user == 0`). An explicit request to balance against a
        # supplied value overrides the freeze; pass `balance = FALSE` to keep a
        # manually set array in that case.
        derive_rate <- is.null(resource_rate)
        derive_capacity <- is.null(resource_capacity)

        balance_fn <- get0(paste0("balance_", params@resource_dynamics))
        if (!is.function(balance_fn)) {
            stop("There is no balancing function available for ",
                 params@resource_dynamics,
                 ". You should not set `balance = TRUE`.")
        }
        balance <- balance_fn(params,
                              resource_rate = resource_rate,
                              resource_capacity = resource_capacity)

        freeze_rate <- num_args_user == 0 && derive_rate &&
            !is.null(comment(params@rr_pp))
        freeze_capacity <- num_args_user == 0 && derive_capacity &&
            !is.null(comment(params@cc_pp))
        if (freeze_rate) {
            signal_info("rr_pp", paste0(
                "The resource rate has been set manually and so it was not ",
                "rebalanced. The resource may no longer replenish at the rate ",
                "at which it is consumed. Use `reset = TRUE` to recalculate ",
                "it from the resource parameters."),
                level = 1, severity = "warning", unhandled = "show")
            resource_rate <- NULL
        } else {
            resource_rate <- balance$resource_rate
        }
        if (freeze_capacity) {
            signal_info("cc_pp", paste0(
                "The resource capacity has been set manually and so it was ",
                "not rebalanced. The resource may no longer replenish at the ",
                "rate at which it is consumed. Use `reset = TRUE` to ",
                "recalculate it from the resource parameters."),
                level = 1, severity = "warning", unhandled = "show")
            resource_capacity <- NULL
        } else {
            resource_capacity <- balance$resource_capacity
        }
    }

    # Set rates
    if (!is.null(resource_rate)) {
        params@rr_pp[] <- resource_rate
        comment(params@rr_pp) <- comment(resource_rate)
    }
    if (!is.null(resource_capacity)) {
        params@cc_pp[] <- resource_capacity
        comment(params@cc_pp) <- comment(resource_capacity)
    }

    params@time_modified <- lubridate::now()
    return(params)
}

#' @rdname setResource
#' @return A vector with the intrinsic resource birth rate for each size class.
#' @export
resource_rate <- function(params) {
    ArrayResourceBySize(params@rr_pp, value_name = "Resource birth rate",
                        units = "1/year", params = params)
}

#' @rdname setResource
#' @param value The desired new value for the respective parameter.
#' @export
`resource_rate<-` <- function(params, balance = NULL, value) {
    setResource(params, resource_rate = value, balance = balance)
}

#' @rdname setResource
#' @return A vector with the intrinsic resource capacity for each size class.
#' @export
resource_capacity <- function(params) {
    ArrayResourceBySize(params@cc_pp, value_name = "Resource capacity",
                        units = "1/g", type = "density", params = params)
}

#' @rdname setResource
#' @export
`resource_capacity<-` <- function(params, balance = NULL, value) {
    setResource(params, resource_capacity = value, balance = balance)
}


#' @rdname setResource
#' @return A vector with the ratio between the current resource number density
#'   and the resource capacity for each size class.
#' @export
resource_level <- function(params) {
    ArrayResourceBySize(params@initial_n_pp / params@cc_pp,
                        value_name = "Resource level", units = "",
                        type = "proportion", params = params)
}

#' @rdname setResource
#' @export
`resource_level<-` <- function(params, balance = NULL, value) {
    setResource(params, resource_level = value, balance = balance)
}


#' @rdname setResource
#' @return The name of the function that determines the resource dynamics.
#' @export
resource_dynamics <- function(params) {
    params@resource_dynamics
}


#' @rdname setResource
#' @export
#' @examples
#' params <- NS_params
#' resource_dynamics(params)
#' resource_dynamics(params) <- "resource_constant"
`resource_dynamics<-` <- function(params, balance = NULL, value) {
    setResource(params, resource_dynamics = value, balance = balance)
}

#' Default weight-length parameters for the resource
#'
#' The resource is a composite of everything from bacteria to
#' macrozooplankton, so it has no taxonomic length-weight relationship. The
#' default is the geometric one that plankton ecology uses instead: the
#' **equivalent spherical diameter** of an organism with the density of water,
#' \deqn{w = \frac{\pi}{6} l^3,}
#' with \eqn{w} in grams and \eqn{l} in centimetres. On a mizer size grid this
#' puts the smallest resource sizes at a fraction of a micrometre and a
#' milligram organism at about a millimetre, which is the right order for
#' bacteria and copepods respectively.
#'
#' Note that this is a different convention from the one the species use: a fish
#' of a given weight is longer than a sphere of the same weight, by a factor
#' \eqn{(a_{fish}/a_{resource})^{-1/3}}, about 3.7 for the mizer default
#' `a = 0.01`. That difference is real rather than an artefact — a 1 mg copepod
#' really is shorter than a 1 mg fish larva — but it does mean the resource and
#' the species sit on the plot at their own conventions.
#'
#' @format A list with entries `a` and `b`.
#' @seealso [resource_params()]
#' @keywords internal
resource_length_defaults <- list(a = pi / 6, b = 3)

#' The weight-length parameters of the resource
#'
#' Reads `a` and `b` from [resource_params()], falling back to
#' [resource_length_defaults] for a model that does not set them — which is
#' every model built before these parameters existed.
#'
#' @param params A MizerParams object.
#' @return A list with entries `a` and `b`.
#' @keywords internal
resource_length_params <- function(params) {
    rp <- params@resource_params
    list(a = rp[["a"]] %||% resource_length_defaults$a,
         b = rp[["b"]] %||% resource_length_defaults$b)
}

# The resource capacity that the scalar resource parameters imply, or NULL if
# they do not determine it. Used both to recompute a capacity that mizer
# controls and to see whether a frozen one has fallen out of step with the
# parameters, see `setResource()`.
capacity_from_resource_params <- function(params, lambda) {
    kappa <- params@resource_params[["kappa"]]
    if (is.null(kappa) || !is.numeric(kappa) || length(kappa) != 1) {
        return(NULL)
    }
    w_full <- params@w_full
    w_pp_cutoff <- params@resource_params$w_pp_cutoff
    if (isTRUE(params@second_order_w[["bin_average"]])) {
        return(kappa * power_law_bin_average(w_full, params@dw_full, -lambda,
                                             w_max = w_pp_cutoff))
    }
    capacity <- kappa * w_full ^ (-lambda)
    capacity[w_full >= w_pp_cutoff] <- 0
    capacity
}

# The resource replenishment rate that the scalar resource parameters imply,
# or NULL if they do not determine it.
rate_from_resource_params <- function(params, n) {
    r_pp <- params@resource_params[["r_pp"]]
    if (is.null(r_pp) || !is.numeric(r_pp) || length(r_pp) != 1) {
        return(NULL)
    }
    if (isTRUE(params@second_order_w[["bin_average"]])) {
        return(r_pp * power_law_bin_average(params@w_full, params@dw_full,
                                            n - 1))
    }
    r_pp * params@w_full ^ (n - 1)
}