From ffd2805740be70bfce324c28a774d13900443665 Mon Sep 17 00:00:00 2001 From: Gustav Delius Date: Mon, 5 Jul 2021 16:05:00 +0100 Subject: [PATCH] Improvements to `setBevertonHolt()` to make it easier to set parameters for only some species. More robust code. Added tests and examples. --- R/setBevertonHolt.R | 116 ++++++++++++++++---------- docs/reference/setBevertonHolt.html | 59 ++++++++++++- man/setBevertonHolt.Rd | 26 +++++- tests/testthat/test-setBevertonHolt.R | 83 ++++++++++++++++-- 4 files changed, 225 insertions(+), 59 deletions(-) diff --git a/R/setBevertonHolt.R b/R/setBevertonHolt.R index 888bd245..00b932d6 100644 --- a/R/setBevertonHolt.R +++ b/R/setBevertonHolt.R @@ -89,9 +89,12 @@ #' ratio between the density-dependent reproduction rate \eqn{R_{dd}}{R_dd} and #' the maximal reproduction rate \eqn{R_{max}}{R_max}. #' -#' The parameter you provide can be either a vector, with one value for each -#' species, or a single value that is then used for all species. The -#' values for `R_max` must be larger than \eqn{R_{dd}}{R_dd} and can range +#' The parameter you provide can be either a vector with one value for each +#' species, or a named vector where the names determine which species are +#' affected, or a single unnamed value that is then used for all species. Any +#' species for which the given value is `NA` will remain unaffected. +#' +#' The values for `R_max` must be larger than \eqn{R_{dd}}{R_dd} and can range #' up to `Inf`. The values for the `reproduction_level` must be positive and #' less than 1. The values for `erepro` must be large enough to allow the #' required reproduction rate. If a smaller value is requested a warning is @@ -116,6 +119,22 @@ #' `reproduction_level = 1 / R_factor` instead. #' #' @return A MizerParams object +#' @examples +#' params <- NS_params +#' species_params(params)$erepro +#' # Attempting to set the same erepro for all species +#' params <- setBevertonHolt(params, erepro = 0.1) +#' t(species_params(params)[, c("erepro", "R_max")]) +#' # Setting erepro for some species +#' params <- setBevertonHolt(params, erepro = c("Gurnard" = 0.6, "Plaice" = 0.95)) +#' t(species_params(params)[, c("erepro", "R_max")]) +#' # Setting R_max +#' R_max <- 1e17 * species_params(params)$w_inf^-1 +#' params <- setBevertonHolt(NS_params, R_max = R_max) +#' t(species_params(params)[, c("erepro", "R_max")]) +#' # Setting reproduction_level +#' params <- setBevertonHolt(params, reproduction_level = 0.3) +#' t(species_params(params)[, c("erepro", "R_max")]) #' @export setBevertonHolt <- function(params, erepro, R_max, reproduction_level, R_factor) { @@ -126,86 +145,93 @@ setBevertonHolt <- function(params, erepro, R_max, reproduction_level, hasArg("R_factor") != 1) { stop("You should only provide `params` and one other argument.") } + # No matter which argument is given, I want to manipulate the values + if (hasArg("erepro")) values <- erepro + if (hasArg("R_max")) values <- R_max + if (hasArg("reproduction_level")) values <- reproduction_level + if (hasArg("R_factor")) values <- R_factor - rdd <- getRDD(params) - rdi <- getRDI(params) - if (any(rdi == 0)) { + if (length(values) == 1 && is.null(names(values))) { + values <- rep(values, no_sp) + } + if (length(values) != no_sp && is.null(names(values))) { + stop("You need to supply a vector of length ", no_sp, + " or a single number or a named vector.") + } + if (is.null(names(values))) { + names(values) <- params@species_params$species + } + values <- values[!is.na(values)] + if (length(values) == 0) return(params) # Nothing to do + if (!all(is.numeric(values))) { + stop("You provided invalid non-numeric values.") + } + # select the species that are affected + species <- valid_species_arg(params, names(values)) + sp_idx <- match(species, params@species_params$species) + + rdd <- getRDD(params)[species] + rdi <- getRDI(params)[species] + if (any(rdi == 0)) { # This should never happen, but did happen in the past. stop("Some species have no reproduction.") } params@rates_funcs$RDD <- "BevertonHoltRDD" - # Calculate required rdd - mumu <- getMort(params) - gg <- getEGrowth(params) - rdd_new <- getRequiredRDD(params) + rdd_new <- getRequiredRDD(params)[species] if (!missing(erepro)) { - if (length(erepro) != 1 && length(erepro) != no_sp) { - stop("`erepro` needs to be a vector of length ", no_sp, - " or a single number.") - } - rdi_new <- rdi * erepro / params@species_params$erepro + erepro_new <- values + erepro_old <- params@species_params$erepro[sp_idx] + rdi_new <- rdi * erepro_new / erepro_old wrong <- rdi_new < rdd_new if (any(wrong)) { warning("For the following species the requested `erepro` ", "was too small and has been increased to the smallest ", "possible value: ", - paste(params@species_params$species[wrong], collapse = ", ")) + paste(species[wrong], collapse = ", ")) rdi_new[wrong] <- rdd_new[wrong] - erepro[wrong] <- (params@species_params$erepro * rdi_new / rdi)[wrong] + erepro_new[wrong] <- (erepro_old * rdi_new / rdi)[wrong] } r_max_new <- rdi_new * rdd_new / (rdi_new - rdd_new) r_max_new[is.nan(r_max_new)] <- Inf - params@species_params$erepro <- erepro - params@species_params$R_max <- r_max_new + params@species_params$erepro[sp_idx] <- erepro_new + params@species_params$R_max[sp_idx] <- r_max_new return(params) } # We now know that we are setting R_max, which however can have been # specified in different ways. if (!missing(reproduction_level)) { - if (length(reproduction_level) != 1 && length(reproduction_level) != no_sp) { - stop("`reproduction_level` needs to be a vector of length ", no_sp, - " or a single number.") - } - if (!all(reproduction_level > 0 & reproduction_level < 1)) { + if (!all(values > 0 & values < 1)) { stop("The reproduction level must be a number strictly between 0 and 1.") } - r_max_new <- rdd_new / reproduction_level + r_max_new <- rdd_new / values } if (!missing(R_factor)) { - if (length(R_factor) != 1 && length(R_factor) != no_sp) { - stop("`R_factor` needs to be a vector of length ", no_sp, - " or a single number.") - } - if (!all(R_factor > 1)) { + if (!all(values > 1)) { stop("The R_factor must be greater than 1.") } - r_max_new <- rdd_new * R_factor + r_max_new <- rdd_new * values } if (!missing(R_max)) { - if (length(R_max) != 1 && length(R_max) != no_sp) { - stop("`R_max` needs to be a vector of length ", no_sp, - " or a single number.") - } - wrong <- R_max < rdd_new + wrong <- values < rdd_new if (any(wrong)) { - stop("For the following species the `R_max` you have provided is - too small: ", - paste(params@species_params$species[wrong], collapse = ", ")) + stop("For the following species the `R_max` you have provided is ", + "too small: ", + paste(species[wrong], collapse = ", ")) } - r_max_new <- R_max + r_max_new <- values } rdi_new <- rdd_new / (1 - rdd_new / r_max_new) - params@species_params$R_max <- r_max_new - params@species_params$erepro <- params@species_params$erepro * - rdi_new / rdi - wrong <- params@species_params$erepro > 1 + params@species_params$R_max[sp_idx] <- r_max_new + params@species_params$erepro[sp_idx] <- + params@species_params$erepro[sp_idx] * rdi_new / rdi + wrong <- params@species_params$erepro[sp_idx] > 1 if (any(wrong)) { warning("The following species require an unrealistic reproductive ", "efficiency greater than 1: ", - paste(params@species_params$species[wrong], collapse = ", ")) + paste(species[wrong], collapse = ", ")) } return(params) diff --git a/docs/reference/setBevertonHolt.html b/docs/reference/setBevertonHolt.html index f8fbdb27..29051915 100644 --- a/docs/reference/setBevertonHolt.html +++ b/docs/reference/setBevertonHolt.html @@ -69,6 +69,8 @@ + + @@ -233,9 +235,11 @@

Details R_max you can alternatively specify the reproduction_level which is the ratio between the density-dependent reproduction rate \(R_{dd}\) and the maximal reproduction rate \(R_{max}\).

-

The parameter you provide can be either a vector, with one value for each -species, or a single value that is then used for all species. The -values for R_max must be larger than \(R_{dd}\) and can range +

The parameter you provide can be either a vector with one value for each +species, or a named vector where the names determine which species are +affected, or a single unnamed value that is then used for all species. Any +species for which the given value is NA will remain unaffected.

+

The values for R_max must be larger than \(R_{dd}\) and can range up to Inf. The values for the reproduction_level must be positive and less than 1. The values for erepro must be large enough to allow the required reproduction rate. If a smaller value is requested a warning is @@ -249,6 +253,55 @@

Details energy income. As a result the species will also be less sensitive to fishing, leading to a higher F_MSY.

+

Examples

+ 
params <- NS_params
+species_params(params)$erepro
+
#>  [1] 1 1 1 1 1 1 1 1 1 1 1 1
# Attempting to set the same erepro for all species
+params <- setBevertonHolt(params, erepro = 0.1)
+
#> Warning: For the following species the requested `erepro` was too small and has been increased to the smallest possible value: Gurnard, Plaice
t(species_params(params)[, c("erepro", "R_max")])
+
#>               Sprat      Sandeel       N.pout      Herring          Dab
+#> erepro 1.000000e-01 1.000000e-01 1.000000e-01 1.000000e-01 1.000000e-01
+#> R_max  8.071481e+11 4.112049e+11 3.472063e+13 1.197577e+12 1.167176e+10
+#>            Whiting         Sole   Gurnard    Plaice      Haddock          Cod
+#> erepro 1.00000e-01 1.000000e-01 0.5582259 0.9212325 1.000000e-01          0.1
+#> R_max  6.22081e+11 4.007876e+10       Inf       Inf 3.929056e+12 8280106764.0
+#>              Saithe
+#> erepro 1.000000e-01
+#> R_max  1.145835e+11
# Setting erepro for some species
+params <- setBevertonHolt(params, erepro = c("Gurnard" = 0.6, "Plaice" = 0.95))
+t(species_params(params)[, c("erepro", "R_max")])
+
#>               Sprat      Sandeel       N.pout      Herring          Dab
+#> erepro 1.000000e-01 1.000000e-01 1.000000e-01 1.000000e-01 1.000000e-01
+#> R_max  8.071481e+11 4.112049e+11 3.472063e+13 1.197577e+12 1.167176e+10
+#>            Whiting         Sole      Gurnard       Plaice      Haddock
+#> erepro 1.00000e-01 1.000000e-01 6.000000e-01 9.500000e-01 1.000000e-01
+#> R_max  6.22081e+11 4.007876e+10 1.047481e+13 1.082568e+15 3.929056e+12
+#>                 Cod       Saithe
+#> erepro          0.1 1.000000e-01
+#> R_max  8280106764.0 1.145835e+11
# Setting R_max
+R_max <- 1e17 * species_params(params)$w_inf^-1
+params <- setBevertonHolt(NS_params, R_max = R_max)
+
#> Warning: The following species require an unrealistic reproductive efficiency greater than 1: Plaice
t(species_params(params)[, c("erepro", "R_max")])
+
#>               Sprat      Sandeel       N.pout      Herring          Dab
+#> erepro 9.274305e-03 1.297184e-04 7.257409e-02 8.045063e-03 4.224791e-03
+#> R_max  3.030303e+15 2.777778e+15 1.000000e+15 2.994012e+14 3.086420e+14
+#>             Whiting         Sole      Gurnard       Plaice      Haddock
+#> erepro 1.292557e-02 3.571380e-03 5.609587e-01 3.773957e+01 6.020565e-02
+#> R_max  8.389262e+13 1.154734e+14 1.497006e+14 3.360215e+13 2.316692e+13
+#>                 Cod       Saithe
+#> erepro 6.375069e-05 2.433424e-03
+#> R_max  2.509328e+12 2.521521e+12
# Setting reproduction_level
+params <- setBevertonHolt(params, reproduction_level = 0.3)
+
#> Warning: The following species require an unrealistic reproductive efficiency greater than 1: Plaice
t(species_params(params)[, c("erepro", "R_max")])
+
#>               Sprat      Sandeel       N.pout      Herring          Dab
+#> erepro 1.324581e-02 1.852847e-04 1.026645e-01 1.145066e-02 6.035197e-03
+#> R_max  2.441029e+12 1.368905e+12 3.256200e+13 3.671952e+12 3.726224e+10
+#>             Whiting         Sole      Gurnard       Plaice      Haddock
+#> erepro 1.834577e-02 5.100263e-03 7.974655e-01 1.316046e+00 7.954325e-02
+#> R_max  1.807310e+12 1.288262e+11 2.430980e+12 1.092730e+14 5.804489e+12
+#>                 Cod       Saithe
+#> erepro 9.077210e-05 3.322021e-03
+#> R_max  2.758282e+10 3.730633e+11