R/newSingleSpeciesParams.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

#' Set up parameters for a single species in a power-law background
#'
#' @description
#' `r lifecycle::badge("experimental")`
#'
#' This functions creates a \code{MizerParams} object with a single
#' species. This species is embedded in a fixed power-law community spectrum
#' \deqn{N_c(w) = \kappa w^{-\lambda}}
#' This community provides the food income for the species. Cannibalism is
#' switched off. The predation mortality arises only from the predators in the
#' power-law community and it is assumed that the predators in the community
#' have the same feeding parameters as the foreground species. The function has
#' many arguments, all of which have default values.
#'
#' @details
#' In addition to setting up the parameters, this function also sets up an
#' initial condition that is close to steady state, under the assumption of
#' no fishing.
#'
#' The function rounds `no_w` to the nearest integer and increases it if
#' necessary so that there are at least 5 size bins per factor 10 in body
#' size. It requires `w_min < w_mat < w_max`, `ext_mort_prop` in `[0, 1)`,
#' positive values for `n`, `lambda`, `kappa`, `alpha`, `h`, `beta`, `sigma`
#' and `f0`, and `fc` between 0 and `f0` if `fc` is supplied. If `gamma` is
#' supplied then `f0` is ignored after its value has been validated. The
#' function stops if the resulting feeding level is not sufficient to maintain
#' the species.
#'
#' The returned model has a single foreground species with cannibalism switched
#' off and a fixed power-law background community that provides both food and
#' predation mortality. The initial species spectrum is scaled so that its
#' maximum abundance is half the background abundance at the corresponding
#' size, and `erepro` is then adjusted so the initial state satisfies the egg
#' boundary condition.
#'
#' The diffusion rate is set to `0`. Because growth is therefore deterministic,
#' no individual grows beyond `w_repro_max`, the size at which all available
#' energy is invested into reproduction. The upper boundary of the size grid is
#' therefore placed at that size, so that `w_max = w_repro_max`, instead of the
#' `1.5 * w_repro_max` headroom that [newMultispeciesParams()] leaves to
#' accommodate the stochastic growth produced by diffusion. This choice will be
#' revisited once these constructors gain a diffusion parameter, see
#' \url{https://github.com/sizespectrum/mizer/issues/339}.
#'
#' @param species_name A string with a name for the species. Will be used in
#'   plot legends.
#' @param w_max Maximum size of species
#' @param w_min Egg size of species
#' @param eta Ratio between maturity size \code{w_mat} and maximum size
#'   \code{w_max}. Default is 10^(-0.6), approximately 1/4. Ignored if
#'   \code{w_mat} is supplied explicitly.
#' @param w_mat Maturity size of species. Default value is
#'   \code{eta * w_max}.
#' @param k_vb `r lifecycle::badge("deprecated")` The von Bertalanffy growth
#'   parameter.
#' @param w_inf `r lifecycle::badge("deprecated")` The argument has been
#'   renamed to `w_max`.
#' @inheritParams newTraitParams
#' @export
#' @importFrom lifecycle deprecated
#' @return An object of type \code{MizerParams}
#' @family functions for setting up models
#' @examples
#' params <- newSingleSpeciesParams()
#' sim <- project(params, t_max = 5, effort = 0)
#' plotSpectra(sim)
newSingleSpeciesParams <-
    function(species_name = "Target species",
             w_max = 100,
             w_min = 0.001,
             eta = 10^(-0.6),
             w_mat = w_max * eta,
             no_w = log10(w_max / w_min) * 20 + 1,
             n = 3 / 4,
             p = n,
             lambda = 2.05,
             kappa = 0.005,
             alpha = 0.4,
             h = 30,
             beta = 100,
             sigma = 1.3,
             f0 = 0.6,
             fc = 0.25,
             ks = NA,
             gamma = NA,
             ext_mort_prop = 0,
             reproduction_level = 0,
             second_order_w = FALSE,
             info_level = default_info_level(),
             R_factor = deprecated(),
             w_inf = deprecated(),
             k_vb = deprecated()) {
        with_info_level(info_level = info_level, {
        if (lifecycle::is_present(w_inf)) {
            lifecycle::deprecate_warn(
                when = "2.4.0.0",
                what = "newSingleSpeciesParams(w_inf)",
                with = "newSingleSpeciesParams(w_max)"
            )
            w_max <- w_inf
            w_mat <- w_max * eta
            no_w <- log10(w_max / w_min) * 20 + 1
        }
        if (lifecycle::is_present(k_vb)) {
            lifecycle::deprecate_warn(
                when = "2.4.0.0",
                what = "newSingleSpeciesParams(k_vb)",
                details = "Specify 'h' instead."
            )
        }
    assert_that(is.string(species_name), length(species_name) == 1)

    ## Much of the following code is copied from newTraitParams

    ## Check validity of parameters ----
    if (ext_mort_prop >= 1 || ext_mort_prop < 0) {
        stop("ext_mort_prop can not take the value ", ext_mort_prop,
             " because it should be the proportion of the total mortality",
             " coming from sources other than predation.")
    }
    if (w_min <= 0) {
        stop("The egg size w_min must be greater than zero.")
    }
    no_w <- round(no_w)
    if (no_w < 1) {
        stop("The number of size bins no_w must be a positive integer")
    }
    if (no_w < log10(w_max/w_min)*5) {
        no_w <- round(log10(w_max / w_min) * 5 + 1)
        signal_info("no_w", paste("Increased no_w to", no_w,
                                  "so that there are 5 bins ",
                                  "for an interval from w and 10w."),
                    level = 1, unhandled = "show")
    }
    if (no_w > 10000) {
        signal_info("no_w", paste0("Running a simulation with ", no_w,
                                   " size bins is going to be very slow."),
                    level = 1, unhandled = "show")
    }
    if (w_min >= w_mat) {
        stop("The egg size w_min must be smaller than ",
             "the maturity size w_mat")
    }
    if (w_mat >= w_max) {
        stop("The maturity size w_mat must be ",
             "smaller the maximum size w_max")
    }
    check_f0(f0) # nolint: object_usage_linter.
    if (!all(c(n, lambda, kappa, alpha, h, beta, sigma, f0) > 0)) {
        stop("The parameters n, lambda, kappa, alpha, h, beta, sigma ",
             "and f0, if supplied, need to be positive.")
    }
    if (!is.na(fc) && (fc < 0 || fc > f0)) {
        stop("The critical feeding level must lie between 0 and f0")
    }
    if (!is.na(gamma)) {  # If gamma is supplied, f0 is ignored
        f0 <- NA
    }
        if (lifecycle::is_present(R_factor)) {
        lifecycle::deprecate_soft("2.3.0", "newSingleSpeciesParams(R_factor)",
                                  "newSingleSpeciesParams(reproduction_level)",
                                  "Set `reproduction_level = 1 / R_factor`.")
        reproduction_level <- 1 / R_factor
    }

    ## Build Params Object ----
    erepro <- 0.1  # Will be changed later to achieve coexistence
    species_params <- data.frame(
        species = species_name,
        w_min = w_min,
        # Without diffusion nothing grows beyond w_repro_max (= w_inf here), so
        # the computational grid boundary w_max is placed there too, with no
        # headroom above the maximum size. To be revisited with diffusion (#339).
        w_inf = w_max,
        w_max = w_max,
        w_mat = w_mat,
        w_min_idx = 1,
        h = h,
        gamma = gamma,
        ks = ks,
        f0 = f0,
        fc = fc,
        beta = beta,
        sigma = sigma,
        z0 = 0,
        alpha = alpha,
        erepro = erepro,
        stringsAsFactors = FALSE
    )
    # Build with the requested bin-averaging but the robust upwind flux; the
    # chosen flux scheme is for projection only and is activated at the end,
    # after the construction-time steady-state solve below.
    target_sow <- resolve_second_order_w(second_order_w)
    params <-
        newMultispeciesParams(
            species_params,
            min_w = w_min,
            no_w = no_w,
            max_w = w_max,
            w_pp_cutoff = Inf,
            lambda = lambda,
            kappa = kappa,
            n = n,
            p = p,
            resource_dynamics = "resource_constant",
            second_order_w = c(bin_average = target_sow[["bin_average"]]),
            info_level = 0
        )
    # No cannibalism
    params@interaction[] <- 0
    # No fishing
    params@initial_effort[] <- 0

    w <- params@w
    dw <- params@dw
    h <- params@species_params[["h"]]
    ks <- params@species_params$ks
    f0 <- get_f0_default(params)

    ## Construct steady state solution ----

    # Get constants for steady-state solution
    # Predation mortality rate coefficient
    mu0 <- get_power_law_mort(params)
    # Add backgound mortality rate
    mu0 <- mu0 / (1 - ext_mort_prop)
    hbar <- alpha * h * f0 - ks
    if (hbar < 0) {
        stop("The feeding level is not sufficient to maintain the fish.")
    }
    pow <- mu0 / hbar / (1 - n)
    if (pow < 1) {
        signal_info("mu0", paste0(
            "The ratio of death rate to growth rate is too small, leading to ",
            "an accumulation of fish at their largest size."),
            level = 1, unhandled = "show")
    }

    initial_n <- params@psi  # get array with correct dimensions and names
    initial_n[, ] <- 0
    params@species_params$z_ext <- mu0
    params@given_species_params$z_ext <- mu0
    params <- setExtMort(params)
    gg <- hbar * w^n * (1 - params@psi[1, ])  # Growth rate
    # Steady state solution of the upwind-difference scheme used in project
    growth_matrix <- matrix(gg, nrow = 1)
    mort_matrix <- params@mu_b
    diffusion_matrix <- growth_matrix
    diffusion_matrix[] <- 0
    n_exact <- get_steady_state_n(params, growth_matrix, mort_matrix,
                                  diffusion_matrix, c(1))
    initial_n[1, ] <- n_exact[1, ]

    # The resource was already set up by newMultispeciesParams()
    initial_n_pp <- params@initial_n_pp

    # Normalise abundance so that the maximum ratio between the species
    # abundance and community abundance is 1/2
    fish <- (length(params@w_full) - length(params@w) + 1):length(params@w_full)
    # index in fish spectrum
    imax <- which.max(initial_n[1, ] / initial_n_pp[fish])
    # corresponding resource index
    pmax <- imax + length(params@w_full) - length(params@w)
    initial_n <- initial_n / initial_n[1, imax] * initial_n_pp[pmax] / 2
    params@initial_n <- initial_n

    ## Set reproduction to meet boundary condition ----
    params@species_params$erepro <- params@species_params$erepro *
        getRequiredRDD(params) / getRDI(params)
    params@given_species_params$erepro <- params@species_params$erepro

    reproduction_level(params) <- reproduction_level

    # Activate the chosen advective-flux scheme now construction is done.
    params@second_order_w[["flux"]] <- target_sow[["flux"]]

    return(params)
    })
}