R/deprecated.R

In mizer: Dynamic Multi-Species Size Spectrum Modelling

#' Deprecated obsolete function for setting up multispecies parameters
#'
#' @description
#' `r lifecycle::badge("deprecated")`
#'
#' This function has been deprecated in favour of the function
#' [newMultispeciesParams()] that sets better default values.
#'
#' This wrapper keeps the legacy defaults and also fills in several columns in
#' `species_params` if they are missing, using the same rules as older mizer
#' versions.
#'
#' @inheritParams newMultispeciesParams
#' @inheritParams newTraitParams
#' @param q Allometric exponent of search volume
#' @param ... Further arguments passed to [newMultispeciesParams()].
#'
#' @details
#' If `species_params` contains a `w_inf` column then it is copied to `w_max`.
#' If `max_w` is not supplied then it is set to `1.1 * max(species_params$w_max)`.
#' The supplied `min_w_pp` is shifted up by one grid step before being passed
#' to [newMultispeciesParams()] to compensate for the fact that newer mizer
#' versions extend the full size grid below `min_w_pp`.
#'
#' Missing legacy columns in `species_params` are filled as follows:
#' `gear = species`, `k = 0`, `alpha = 0.6`, `erepro = 1`,
#' `sel_func = "knife_edge"`, `knife_edge_size = w_mat` if needed,
#' `catchability = 1`, `ks = h * 0.2`, and `m = 1`.
#' If `h` is missing it is calculated from `k_vb`, `alpha`, `f0` and `w_max`.
#' If `gamma` is missing it is calculated from `f0`, `h`, `beta`, `sigma`,
#' `lambda` and `kappa`.
#' @return A MizerParams object
#' @export
#' @concept deprecated
set_multispecies_model <-
    function(
        species_params,
        interaction = matrix(1,
                             nrow = nrow(species_params),
                             ncol = nrow(species_params)),
        min_w_pp = 1e-10,
        min_w = 0.001,
        max_w = NULL,
        no_w = 100,
        n = 2 / 3,
        q = 0.8,
        f0 = 0.6,
        kappa = 1e11,
        lambda = 2 + q - n,
        r_pp = 10,
        ...) {
    lifecycle::deprecate_soft("2.0.0", "set_multispecies_model()", "newMultispeciesParams()")
    if (exists("no_w_pp")) {
        warning("New mizer code does not support the parameter no_w_pp")
    }
    # w_inf was renamed to w_max
    if ("w_inf" %in% names(species_params)) {
        species_params$w_max <- species_params$w_inf
    }
    if (is.null(max_w)) {
        max_w <- max(species_params$w_max) * 1.1
    }

    # Need to correct for the fact that new mizer extends w_full to BELOW
    # min_w_pp
        dx <- log10(max_w / min_w) / (no_w - 1)
        min_w_pp <- min_w_pp * 10 ^ dx

    species_params[["q"]] <- q
    species_params[["f0"]] <- f0
    object <- species_params
    # old code from MizerParams() in version 1.0.1

    # Set default values for column values if missing
    # If no gear_name column in object, then named after species
    if (!("gear" %in% colnames(object)))
        object$gear <- object$species
    # If no k column (activity coefficient) in object, then set to 0
    if (!("k" %in% colnames(object)))
        object$k <- 0
    # If no alpha column in object, then set to 0.6
    # Should this be a column? Or just an argument?
    if (!("alpha" %in% colnames(object)))
        object$alpha <- 0.6
    # If no erepro column in object, then set to 1
    if (!("erepro" %in% colnames(object)))
        object$erepro <- 1
    # If no sel_func column in species_params, set to 'knife_edge'
    if (!("sel_func" %in% colnames(object))) {
        message("Note: No sel_func column in species data frame. Setting selectivity to be 'knife_edge' for all species.")
        object$sel_func <- 'knife_edge'
        # Set default selectivity size
        if (!("knife_edge_size" %in% colnames(object))) {
            message("Note: No knife_edge_size column in species data frame. Setting knife edge selectivity equal to w_mat.")
            object$knife_edge_size <- object$w_mat
        }
    }
    # If no catchability column in species_params, set to 1
    if (!("catchability" %in% colnames(object)))
        object$catchability <- 1
    # Sort out h column If not passed in directly, is calculated from f0 and
    # k_vb if they are also passed in
    if (!("h" %in% colnames(object))) {
        message("Note: No h column in species data frame so using f0 and k_vb to calculate it.")
        if (!("k_vb" %in% colnames(object))) {
            stop("Except I can't because there is no k_vb column in the species data frame")
        }
        object$h <- ((3 * object$k_vb) / (object$alpha * f0)) * (object$w_max ^ (1/3))
    }
    # Sorting out gamma column
    if (!("gamma" %in% colnames(object))) {
        message("Note: No gamma column in species data frame so using f0, h, beta, sigma, lambda and kappa to calculate it.")
        ae <- sqrt(2*pi) * object$sigma * object$beta^(lambda-2) * exp((lambda-2)^2 * object$sigma^2 / 2)
        object$gamma <- (object$h / (kappa * ae)) * (f0 / (1 - f0))
    }
    # Sort out ks column
    if (!("ks" %in% colnames(object))) {
        message("Note: No ks column in species data frame so using ks = h * 0.2.")
        object$ks <- object$h * 0.2
    }

    # The m column did not exist in the old version, set it to 1
    if (!("m" %in% colnames(object))) {
        message("Note: No m column in species data frame so using m = 1.")
        object[["m"]] <- 1
    }

    return(newMultispeciesParams(object,
                                 interaction = interaction,
                                 n = n,
                                 kappa = kappa,
                                 lambda = lambda,
                                 resource_rate = r_pp,
                                 max_w = max_w,
                                 min_w = min_w,
                                 min_w_pp = min_w_pp,
                                 no_w = no_w,
                                 ...))
}

#' Alias for `set_multispecies_model()`
#'
#' `r lifecycle::badge("deprecated")`
#' An alias provided for backward compatibility with mizer version <= 1.0
#' @inherit set_multispecies_model
#' @export
#' @concept deprecated
MizerParams <- set_multispecies_model


# Copied from version 1.0.1
#' Deprecated function for setting up parameters for a trait-based model
#'
#' @description
#' `r lifecycle::badge("deprecated")`
#'
#' This function has been deprecated in favour of the function
#' [newTraitParams()] that sets better default values.
#'
#' @details
#' This functions creates a `MizerParams` object so that trait-based-type
#' models can be easily set up and run. The trait-based size spectrum model can
#' be derived as a simplification of the general size-based model used in
#' `mizer`. The species-specific parameters are the same for all species,
#' except for
#' the asymptotic size, which is considered the most important trait
#' characterizing a species. Other parameters are related to the asymptotic
#' size. For example, the size at maturity is given by \code{w_max * eta},
#' where `eta` is
#' the same for all species. For the trait-based model the number of species is
#' not important. For applications of the trait-based model see Andersen &
#' Pedersen (2010). See the `mizer` vignette for more details and examples
#' of the trait-based model.
#'
#' The function has many arguments, all of which have default values. Of
#' particular interest to the user are the number of species in the model and
#' the minimum and maximum asymptotic sizes. The asymptotic sizes of the species
#' are spread evenly on a logarithmic scale within this range.
#'
#' The stock recruitment relationship is the default Beverton-Holt style. The
#' maximum recruitment is calculated using equilibrium theory (see Andersen &
#' Pedersen, 2010) and a multiplier, `k0`. Users should adjust `k0` to
#' get the spectra they want.
#'
#' The factor for the search volume, `gamma`, is calculated using the
#' expected feeding level, `f0`.
#'
#' Fishing selectivity is modelled as a knife-edge function with one parameter,
#' `knife_edge_size`, which is the size at which species are selected. Each
#' species can either be fished by the same gear (`knife_edge_size` has a
#' length of 1) or by a different gear (the length of `knife_edge_size` has
#' the same length as the number of species and the order of selectivity size is
#' that of the asymptotic size).
#'
#' The resulting `MizerParams` object can be projected forward using
#' `project` like any other `MizerParams` object. When projecting
#' the community model it may be necessary to reduce `dt` to 0.1 to avoid
#' any instabilities with the solver. You can check this by plotting the biomass
#' or abundance through time after the projection.
#'
#' @param no_sp The number of species in the model. The default value is 10. The
#'   more species, the longer takes to run.
#' @param min_w_inf The asymptotic size of the smallest species in the
#'   community.
#' @param max_w_inf The asymptotic size of the largest species in the community.
#' @param no_w The number of size bins in the community spectrum.
#' @param min_w The smallest size of the community spectrum.
#' @param max_w Maximum size of the consumer size grid passed to
#'   [MizerParams()]. Default value is `max_w_inf * 1.1`.
#' @param min_w_pp Smallest size on the resource size grid passed to
#'   [MizerParams()]. Default value is `1e-10`.
#' @param w_pp_cutoff The cut off size of the resource spectrum. Default value
#'    is 1.
#' @param k0 Multiplier for the maximum recruitment. Default value is 50.
#' @param n Scaling of the intake. Default value is 2/3.
#' @param p Scaling of the standard metabolism. Default value is 0.75.
#' @param q Exponent of the search volume. Default value is 0.9.
#' @param eta Factor to calculate `w_mat` from asymptotic size.
#' @param r_pp Growth rate parameter for the resource spectrum. Default value is 4.
#' @param kappa Coefficient in abundance power law. Default value is
#'   0.005.
#' @param lambda Exponent of the abundance power law. Default value is (2+q-n).
#' @param alpha The assimilation efficiency of the community. The default value
#'   is 0.6
#' @param ks Standard metabolism coefficient. Default value is 4.
#' @param z0pre The coefficient of the background mortality of the community. z0
#'   = z0pre * w_inf ^ (n-1). The default value is 0.6.
#' @param h Maximum food intake rate. Default value is 30.
#' @param beta Preferred predator prey mass ratio. Default value is 100.
#' @param sigma Width of prey size preference. Default value is 1.3.
#' @param f0 Expected average feeding level. Used to set `gamma`, the
#'   factor for the search volume. The default value is 0.5.
#' @param gamma Volumetric search rate. Estimated using `h`, `f0` and
#'   `kappa` if not supplied.
#' @param knife_edge_size The minimum size at which the gear or gears select
#'   species. Must be of length 1 or no_sp.
#' @param gear_names The names of the fishing gears. A character vector, the
#'   same length as the number of species. Default is 1 - no_sp.
#' @param ... Other arguments to pass to the `MizerParams` constructor.
#'
#' @return An object of type `MizerParams`
#' @references K. H. Andersen and M. Pedersen, 2010, Damped trophic cascades
#'   driven by fishing in model marine ecosystems. Proceedings of the Royal
#'   Society V, Biological Sciences, 1682, 795-802.
#' @export
#' @concept deprecated
set_trait_model <- function(no_sp = 10,
                            min_w_inf = 10,
                            max_w_inf = 1e5,
                            no_w = 100,
                            min_w = 0.001,
                            max_w = max_w_inf * 1.1,
                            min_w_pp = 1e-10,
                            w_pp_cutoff = 1,
                            k0 = 50, # recruitment adjustment parameter
                            n = 2/3,
                            p = 0.75,
                            q = 0.9,
                            eta = 0.25,
                            r_pp = 4,
                            kappa = 0.005,
                            lambda = 2+q-n,
                            alpha = 0.6,
                            ks = 4,
                            z0pre = 0.6,
                            h = 30,
                            beta = 100,
                            sigma = 1.3,
                            f0 = 0.5,
                            gamma = NA,
                            knife_edge_size = 1000,
                            gear_names = "knife_edge_gear",
                            ...){
    lifecycle::deprecate_soft("2.0.0", "set_trait_model()", "newTraitParams()")
    if (exists("no_w_pp")) {
        warning("New mizer code does not support the parameter no_w_pp")
    }
    # If not supplied, calculate gamma using equation 2.1 in A&P 2010
    # TODO: remove this here because it is already calculated in MizerParams()
    #       Having the same code in two locations is not a good idea
    if(is.na(gamma)){
        alpha_e <- sqrt(2*pi) * sigma * beta^(lambda-2) *
            exp((lambda-2)^2 * sigma^2 / 2) # see A&P 2009
        gamma <- h * f0 / (alpha_e * kappa * (1-f0)) # see A&P 2009
    }
    w_inf <- 10^seq(from=log10(min_w_inf), to = log10(max_w_inf), length=no_sp)
    w_mat <- w_inf * eta

    # Check gears
    if (length(knife_edge_size) > no_sp){
        stop("There cannot be more gears than species in the model")
    }
    if ((length(knife_edge_size) > 1) & (length(knife_edge_size) != no_sp)){
        warning("Number of gears is less than number of species so gear information is being recycled. Is this what you want?")
    }
    if ((length(gear_names) != 1) & (length(gear_names) != no_sp)){
        stop("Length of gear_names argument must equal the number of species.")
    }

    # Make the species parameters data.frame
    trait_params_df <- data.frame(
        species = as.factor(1:no_sp),
        w_inf = w_inf,
        w_mat = w_mat,
        w_min = min_w,
        h = h, # max food intake
        gamma = gamma, # vol. search rate,
        ks = ks,# standard metabolism coefficient,
        beta = beta,
        sigma = sigma,
        z0 = z0pre * w_inf^(n-1), # background mortality
        alpha = alpha,
        #r_max = r_max,
        sel_func = "knife_edge",
        knife_edge_size = knife_edge_size,
        gear = gear_names,
        erepro = 1
    )
    # Make the MizerParams
    trait_params <-
        MizerParams(
            trait_params_df,
            min_w = min_w,
            max_w = max_w,
            no_w = no_w,
            min_w_pp = min_w_pp,
            w_pp_cutoff = w_pp_cutoff,
            n = n,
            p = p,
            q = q,
            r_pp = r_pp,
            kappa = kappa,
            lambda = lambda
        )
    # Sort out maximum recruitment - see A&P 2009 Get max flux at recruitment
    # boundary, R_max R -> | -> g0 N0 R is egg flux, in numbers per time Actual
    # flux at recruitment boundary = RDD = NDD * g0 (where g0 is growth rate) So
    # in our BH SRR we need R_max comparable to RDI (to get RDD) R_max = N0_max
    # * g0 (g0 is the average growth rate of smallest size, i.e. at f0 = 0.5) N0
    # given by Appendix A of A&P 2010 - see Ken's email 12/08/13 Taken from
    # Ken's code 12/08/13 - equation in paper is wrong!
    alpha_p <- f0 * h * beta^(2 * n - q - 1) *
        exp((2 * n * (q - 1) - q^2 + 1) * sigma^2 / 2)
    alpha_rec <- alpha_p / (alpha * h * f0 - ks)
    # Calculating dw using Ken's code - see Ken's email 12/08/13
    tmpA <- w_inf[1]
    tmpB <- (log10(w_inf[length(w_inf)]) - log10(w_inf[1])) / (no_sp - 1) # Difference between logged w_infs, fine
    dw_winf <- tmpB * tmpA *10^(tmpB*((1:no_sp)-1)) # ?
    N0_max <- k0 * w_inf^(n*2-q-3+alpha_rec) * dw_winf  # Why * dw_winf, not / ? Ken confirms * in email
    # No need to include (1 - psi) in growth equation because allocation to reproduction at this size = 0, so 1 - psi = 1
    g0 <- (alpha * f0 * h * trait_params@w[1]^n - ks * trait_params@w[1]^p)
    r_max <- N0_max * g0

    trait_params@species_params$R_max <- r_max

    return(trait_params)
}


# Copied from version 1.0.1
#' Deprecated function for setting up parameters for a community-type model
#'
#' @description
#' `r lifecycle::badge("deprecated")`
#'
#' This function has been deprecated in favour of the function
#' [newCommunityParams()] that sets better default values.
#'
#' @details
#' This functions creates a \code{\linkS4class{MizerParams}} object so that
#' community-type models can be easily set up and run. A community model has
#' several features that distinguish it from the food-web type models. Only one
#' 'species' is resolved, i.e. one 'species' is used to represent the whole
#' community. The resource spectrum only extends to the start of the community
#' spectrum. Recruitment to the smallest size in the community spectrum is
#' constant and set by the user. As recruitment is constant, the proportion of
#' energy invested in reproduction (the slot `psi` of the returned
#' `MizerParams` object) is set to 0. Standard metabolism has been turned
#' off (the parameter `ks` is set to 0). Consequently, the growth rate is
#' now determined solely by the assimilated food (see the package vignette for
#' more details).
#'
#' The function has many arguments, all of which have default values. The main
#' arguments that the users should be concerned with are `z0`,
#' `recruitment`, `alpha` and `f0` as these determine the average
#' growth rate of the community.
#'
#' Fishing selectivity is modelled as a knife-edge function with one parameter,
#' `knife_edge_size`, which determines the size at which species are
#' selected.
#'
#' The resulting `MizerParams` object can be projected forward using
#' \code{project()} like any other `MizerParams` object. When projecting
#' the community model it may be necessary to keep a small time step size
#' `dt` of around 0.1 to avoid any instabilities with the solver. You can
#' check for these numerical instabilities by plotting the biomass or abundance
#' through time after the projection.
#'
#' @param z0 The background mortality of the community. Default value is 0.1.
#' @param alpha The assimilation efficiency of the community. Default value 0.2
#' @param f0 The average feeding level of individuals who feed on a power-law
#'   spectrum. This value is used to calculate the search rate parameter
#'   `gamma` (see the package vignette). Default value is 0.7.
#' @param h The maximum food intake rate. Default value is 10.
#' @param beta The preferred predator prey mass ratio. Default value is 100.
#' @param sigma The width of the prey preference. Default value is 2.0.
#' @param q The search volume exponent. Default value is 0.8.
#' @param n The scaling of the intake. Default value is 2/3.
#' @param kappa The carrying capacity of the resource spectrum. Default value
#'   is 1000.
#' @param lambda The exponent of the resource spectrum. Default value is 2 + q - n.
#' @param r_pp Growth rate parameter for the resource spectrum. Default value is 10.
#' @param gamma Volumetric search rate. Estimated using `h`, `f0` and
#'   `kappa` if not supplied.
#' @param recruitment The constant recruitment in the smallest size class of the
#'   community spectrum. This should be set so that the community spectrum
#'   continues the resource spectrum. Default value = `kappa` *
#'   `min_w`^-`lambda`.
#' @param rec_mult Additional multiplier for the constant recruitment. Default
#'   value is 1.
#' @param knife_edge_size The size at the edge of the knife-selectivity
#'   function. Default value is 1000.
#' @param knife_is_min Is the knife-edge selectivity function selecting above
#'   (TRUE) or below (FALSE) the edge. Default is TRUE.
#' @param max_w The maximum size of the community. The `w_inf` of the
#'   species used to represent the community is set to this value. The
#'   default value is 1e6.
#' @param min_w The minimum size of the community. Default value is 1e-3.
#' @param min_w_pp The smallest size of the resource spectrum.
#' @param ... Other arguments to pass to the `MizerParams` constructor.
#' @export
#' @return An object of type \code{\linkS4class{MizerParams}}
#' @references K. H. Andersen,J. E. Beyer and P. Lundberg, 2009, Trophic and
#'   individual efficiencies of size-structured communities, Proceedings of the
#'   Royal Society, 276, 109-114
#' @concept deprecated
#' @examples
#' \donttest{
#' params <- set_community_model(f0=0.7, z0=0.2, recruitment=3e7)
#' # This is now achieved with
#' params <- newCommunityParams(f0 = 0.7, z0 = 0.2)
#' sim <- project(params, effort = 0, t_max = 100, dt=0.1)
#' plotBiomass(sim)
#' plotSpectra(sim)
#' }
set_community_model <- function(max_w = 1e6,
                                min_w = 1e-3,
                                min_w_pp = 1e-10,
                                z0 = 0.1,
                                alpha = 0.2,
                                h = 10,
                                beta = 100,
                                sigma = 2.0,
                                q = 0.8,
                                n = 2/3,
                                kappa = 1000,
                                lambda = 2+q-n,
                                f0 = 0.7,
                                r_pp = 10,
                                gamma = NA,
                                knife_edge_size = 1000,
                                knife_is_min = TRUE,
                                recruitment = kappa * min_w^-lambda,
                                rec_mult = 1,
                                ...
){
    lifecycle::deprecate_soft("2.0.0", "set_community_model()", "newCommunityParams()")
    w_inf <- max_w * 0.9
    w_pp_cutoff <- min_w
    ks <- 0 # Turn off standard metabolism
    p <- n # But not used as ks = 0
    # Estimate gamma if not supplied
    if (is.na(gamma)){
        gamma <- (f0 * h * beta^(2-lambda)) / ((1-f0)*sqrt(2*pi)*kappa*sigma)
    }
    # Make the species data.frame
    com_params_df <- data.frame(
        species = "Community",
        w_inf = w_inf,
        w_mat = w_inf * 0.25, # Has no affect as psi set to 0 but we set it to
        # something to help the constructor
        h = h, # max food intake
        gamma = gamma,# vol. search rate,
        ks = ks,# standard metabolism coefficient,
        beta = beta,
        sigma = sigma,
        z0 = z0, # background mortality
        alpha = alpha,
        erepro = 1, # not used
        sel_func = "knife_edge",
        knife_edge_size = knife_edge_size,
        knife_is_min = knife_is_min,
        constant_reproduction = recruitment * rec_mult # to be used in the SRR
    )
    com_params <- MizerParams(com_params_df, p = p, n = n, q = q, lambda = lambda,
                              kappa = kappa, min_w = min_w, max_w = max_w,
                              min_w_pp = min_w_pp,
                              w_pp_cutoff = w_pp_cutoff, r_pp = r_pp, ...)
    com_params@rates_funcs$RDD <- "constantRDD"
    com_params@psi[] <- 0 # Need to force to be 0. Can try setting w_mat but
    # due to slope still not 0
    # Set w_mat to NA for clarity - it is not actually being used
    com_params@species_params$w_mat[] <- NA
    return(com_params)
}

#### getPhiPrey ####
#' Get available energy
#'
#' @description
#' `r lifecycle::badge("deprecated")`
#'
#' This is deprecated and is no longer used by the mizer project() method.
#' Calculates the amount \eqn{E_{a,i}(w)} of food exposed to each predator as
#' a function of predator size.
#'
#' @param object An \linkS4class{MizerParams} object
#' @param n A matrix of species abundances (species x size)
#' @param n_pp A vector of the background abundance by size
#' @param ... Other arguments (currently unused)
#'
#' @return A two dimensional array (predator species x predator size)
#'   equal to `getEncounter(object, n, n_pp) / search_vol(object)`.
#' @seealso [project()]
#' @export
#' @concept deprecated
getPhiPrey <- function(object, n, n_pp, ...) {
    lifecycle::deprecate_soft("2.0.0", "getPhiPrey()", "newMultispeciesParams()")
    getEncounter(object, n, n_pp) / object@search_vol
}


#### get-prefixed rate accessors ####
#' Superseded `get`-prefixed aliases for values stored in a model
#'
#' @description
#' `r lifecycle::badge("superseded")`
#'
#' Each of these functions is an alias for a function with a shorter name that
#' returns exactly the same thing. The shorter name is the one that also has a
#' replacement function (`catchability(params) <- value` and friends), so that
#' is the name to use:
#'
#' | Superseded | Use instead |
#' |---|---|
#' | `getCatchability()` | [catchability()] |
#' | `getSelectivity()` | [selectivity()] |
#' | `getInitialEffort()` | [initial_effort()] |
#' | `getInteraction()` | [interaction_matrix()] |
#' | `getResourceDynamics()` | [resource_dynamics()] |
#' | `getResourceLevel()` | [resource_level()] |
#' | `getResourceRate()` | [resource_rate()] |
#' | `getResourceCapacity()` | [resource_capacity()] |
#' | `getPredKernel()` | [pred_kernel()] |
#' | `getSearchVolume()` | [search_vol()] |
#' | `getMaxIntakeRate()` | [intake_max()] |
#' | `getMetabolicRate()` | [metab()] |
#' | `getExtMort()` | [ext_mort()] |
#' | `getExtEncounter()` | [ext_encounter()] |
#' | `getMaturityProportion()` | [maturity()] |
#' | `getReproductionProportion()` | [repro_prop()] |
#' | `getReproductionLevel()` | [reproduction_level()] |
#'
#' The `get` prefix is reserved for the functions that *calculate* a rate from
#' the current state of a model, like [getEncounter()] or [getFMort()]. The
#' functions above only read back a value that is already stored in the
#' MizerParams object, which is what the bare names say.
#'
#' The old names are however not going away. They are plain aliases: they do
#' not warn and they will keep working, so existing code and old scripts run
#' unchanged. They are not used anywhere inside mizer and should not be used in
#' new code.
#'
#' @param params A MizerParams object
#' @return The same as the function it is an alias for.
#' @name superseded_accessors
#' @concept deprecated
NULL

#' @rdname superseded_accessors
#' @export
getCatchability <- catchability

#' @rdname superseded_accessors
#' @export
getSelectivity <- selectivity

#' @rdname superseded_accessors
#' @export
getInitialEffort <- initial_effort

#' @rdname superseded_accessors
#' @export
getInteraction <- interaction_matrix

#' @rdname superseded_accessors
#' @export
getResourceDynamics <- resource_dynamics

#' @rdname superseded_accessors
#' @export
getResourceLevel <- resource_level

#' @rdname superseded_accessors
#' @export
getResourceRate <- resource_rate

#' @rdname superseded_accessors
#' @export
getResourceCapacity <- resource_capacity

#' @rdname superseded_accessors
#' @export
getPredKernel <- pred_kernel

#' @rdname superseded_accessors
#' @export
getSearchVolume <- search_vol

#' @rdname superseded_accessors
#' @export
getMaxIntakeRate <- intake_max

#' @rdname superseded_accessors
#' @export
getMetabolicRate <- metab

#' @rdname superseded_accessors
#' @export
getExtMort <- ext_mort

#' @rdname superseded_accessors
#' @export
getExtEncounter <- ext_encounter

#' @rdname superseded_accessors
#' @export
getMaturityProportion <- maturity

#' @rdname superseded_accessors
#' @export
getReproductionProportion <- repro_prop

#' @rdname superseded_accessors
#' @export
getReproductionLevel <- reproduction_level


#### Superseded steady-state finders ####
#' Superseded names for the steady-state finders
#'
#' @description
#' `r lifecycle::badge("superseded")`
#'
#' Neither of these names says what distinguishes the two functions, and
#' `projectToSteady()` returns a different class depending on an argument. Both
#' have been replaced:
#'
#' | Superseded | Use instead |
#' |---|---|
#' | `steady()` | [tuneSteadyState()] |
#' | `projectToSteady()` | [findSteadyState()], or [projectUntilSettled()] for the trajectory |
#'
#' The distinction the new names carry is what each one holds fixed.
#' [tuneSteadyState()] holds the inputs to the fish dynamics — the reproduction
#' rate and the resource — at the values you supply while the spectra settle,
#' and then adjusts the parameters that generate them (`erepro`/`R_max` and
#' `cc_pp`) so that those held values are steady too. That is what `steady()`
#' always did. [findSteadyState()] changes no parameter and lets reproduction,
#' the resource and the spectra settle together, which is what
#' `projectToSteady()` did.
#'
#' Each of the two also gained a `solver` argument, so the same job can be done
#' either by running the dynamics (`solver = "project"`, the default and the old
#' behaviour) or by solving the steady-state equation directly with a
#' Newton-type root finder (`solver = "newton"`), which converges even at a
#' dynamically unstable steady state.
#'
#' The `return_sim` argument is gone from the new functions:
#' [tuneSteadyState()] and [findSteadyState()] always return a `MizerParams`
#' object and [projectUntilSettled()] always returns a `MizerSim`.
#'
#' The old names are however not going away. They are thin wrappers that
#' reproduce the old behaviour exactly, `return_sim` included: they do not warn
#' and they will keep working, so existing code and old scripts run unchanged.
#' They are not used anywhere inside mizer and should not be used in new code.
#'
#' @inheritParams tuneSteadyState
#' @inheritParams projectUntilSettled
#' @param tol The simulation stops when the relative change in the egg
#'   production RDI over t_per years is less than tol for every species.
#' @param t_per The interval in years at which convergence is checked, and hence
#'   also the interval at which the trajectory is saved when
#'   `return_sim = TRUE`. In [projectUntilSettled()] these two roles have been
#'   separated into `t_check` and `t_save`.
#' @param t_save Has no effect. It briefly controlled how finely the biomass
#'   series used for limit-cycle detection was sampled; that series is now
#'   sampled at every time step, which is what its default `dt` gave.
#' @param ... Further arguments will be passed on to your distance function.
#' @param return_sim If TRUE, the function returns the MizerSim object holding
#'   the result of the simulation run, saved at intervals of `t_per`. If FALSE
#'   (default) the function returns a MizerParams object with the "initial"
#'   slots set to the steady state.
#' @return A `MizerParams` object, or a `MizerSim` object if
#'   `return_sim = TRUE`, in either case carrying the `"convergence"` attribute
#'   described in [projectUntilSettled()].
#' @seealso [tuneSteadyState()], [findSteadyState()], [projectUntilSettled()]
#' @name superseded_steady
#' @concept deprecated
NULL

#' @rdname superseded_steady
#' @export
steady <- function(params, t_max = 100, t_per = 1.5, dt = 0.1, t_save = dt,
                   tol = 0.1 * dt, amplitude_tol = 0.01, amp_rel_tol = 0.01,
                   extinction_threshold = 1e-6, return_sim = FALSE,
                   preserve = c("reproduction_level", "erepro", "R_max"),
                   progress_bar = TRUE,
                   info_level = default_info_level(),
                   method = c("euler", "predictor_corrector", "tr_bdf2")) {
    UseMethod("steady")
}

#' @export
steady.MizerParams <- function(params, t_max = 100, t_per = 1.5, dt = 0.1,
                   t_save = dt,
                   tol = 0.1 * dt, amplitude_tol = 0.01, amp_rel_tol = 0.01,
                   extinction_threshold = 1e-6, return_sim = FALSE,
                   preserve = c("reproduction_level", "erepro", "R_max"),
                   progress_bar = TRUE,
                   info_level = default_info_level(),
                   method = c("euler", "predictor_corrector", "tr_bdf2")) {
    with_info_level(info_level = info_level, {
        params <- validParams(params)
        tune_steady_project(params, effort = params@initial_effort,
                            preserve = preserve,
                            t_max = t_max, t_check = t_per, dt = dt,
                            # The trajectory these wrappers have always returned
                            # is saved every `t_per`, so that is the save
                            # interval they ask for. Their own `t_save` only
                            # ever tuned the cycle detection, which no longer
                            # has a knob.
                            t_save = t_per, distance_tol = tol,
                            amplitude_tol = amplitude_tol,
                            amp_rel_tol = amp_rel_tol,
                            extinction_threshold = extinction_threshold,
                            return_sim = return_sim,
                            # Released code was written against a stopping rule
                            # that used the distance function alone, so that is
                            # what this wrapper keeps. The drift is still
                            # measured and still decides `attractor`, so the
                            # result says what it is; only the decision to keep
                            # running is left as it was.
                            require_steady = FALSE,
                            progress_bar = progress_bar,
                            info_level = info_level, method = method)
    })
}

#' @rdname superseded_steady
#' @export
projectToSteady <- function(params,
                            effort = params@initial_effort,
                            distance_func = distanceSSLogN,
                            t_per = 1.5,
                            t_max = 100,
                            dt = 0.1,
                            t_save = dt,
                            tol = 0.1 * t_per,
                            amplitude_tol = 0.01,
                            amp_rel_tol = 0.1,
                            extinction_threshold = 1e-6,
                            return_sim = FALSE,
                            progress_bar = TRUE,
                            info_level = default_info_level(),
                            method = c("euler", "predictor_corrector",
                                       "tr_bdf2"), ...) {
    UseMethod("projectToSteady")
}

#' @export
projectToSteady.MizerParams <- function(params,
                            effort = params@initial_effort,
                            distance_func = distanceSSLogN,
                            t_per = 1.5,
                            t_max = 100,
                            dt = 0.1,
                            t_save = dt,
                            tol = 0.1 * t_per,
                            amplitude_tol = 0.01,
                            amp_rel_tol = 0.1,
                            extinction_threshold = 1e-6,
                            return_sim = FALSE,
                            progress_bar = TRUE,
                            info_level = default_info_level(),
                            method = c("euler", "predictor_corrector",
                                       "tr_bdf2"), ...) {
    project_until_settled(params, effort = effort,
                          distance_func = distance_func,
                          t_check = t_per, t_max = t_max, dt = dt,
                          # As in steady() above: the trajectory keeps the
                          # `t_per` spacing it has always had.
                          t_save = t_per, distance_tol = tol,
                          amplitude_tol = amplitude_tol,
                          amp_rel_tol = amp_rel_tol,
                          extinction_threshold = extinction_threshold,
                          # As in steady() above: the old stopping rule, the
                          # new report.
                          require_steady = FALSE,
                          progress_bar = progress_bar, info_level = info_level,
                          method = method, ..., return_sim = return_sim)
}