New species param w_repro_max

NEWS.md

@@ -1,5 +1,7 @@
 # Development version
 
+- New species parameter `w_repro_max` giving the size at which a species 
+  invests 100% of its energy into reproduction. Set to `w_max` by default.
 - `getDiet()` now also includes the contribution of the external encounter rate
   to the diet.
 - Improved scaling of the y-axis in `plotGrowthCurves()`.

R/data.R

@@ -9,8 +9,7 @@
 #' @format A data frame with 12 rows and 7 columns. Each row is a species.
 #' \describe{
 #' \item{species}{Name of the species}
-#' \item{w_max}{The size at which the population invests 100% of its income into
-#'   reproduction so that all growth stops.}
+#' \item{w_max}{Maximum size.}
 #' \item{w_mat}{Size at maturity}
 #' \item{beta}{Size preference ratio}
 #' \item{sigma}{Width of the size-preference}
@@ -35,8 +34,7 @@
 #' @format A data frame with 12 rows and 8 columns. Each row is a species.
 #' \describe{
 #' \item{species}{Name of the species}
-#' \item{w_max}{The size at which the population invests 100% of its income into
-#'   reproduction so that all growth stops.}
+#' \item{w_max}{Maximum size.}
 #' \item{w_mat}{Size at maturity}
 #' \item{beta}{Size preference ratio}
 #' \item{sigma}{Width of the size-preference}

R/setReproduction.R

@@ -11,7 +11,11 @@
 #' that is invested into reproduction is the product of two factors: the
 #' proportion `maturity` of individuals that are mature and the proportion
 #' `repro_prop` of the energy available to a mature individual that is 
-#' invested into reproduction.
+#' invested into reproduction. There is a size `w_repro_max` at which all the
+#' energy is invested into reproduction and therefore all growth stops. There
+#' can be no fish larger than `w_repro_max`. If you have not specified the
+#' `w_repro_max` column in the species parameter data frame, then the maximum size
+#' `w_max` is used instead.
 #' 
 #' \subsection{Maturity ogive}{
 #' If the the proportion of individuals that are mature is not supplied via
@@ -29,29 +33,30 @@
 #' default it is chosen as \eqn{U = 10}, however this can be overridden by
 #' including a column \code{w_mat25} in the species parameter dataframe that
 #' specifies the weight at which 25% of individuals are mature, which sets
-#' \eqn{U = \log(3) / \log(w_{mat} / w_{25}).}{U = log(3) / log(w_mat / w_25).}
+#' \eqn{U = \log(3) / \log(w_{mat} / w_{mat25}).}{U = log(3) / log(w_mat /
+#' w_mat25).}
 #' 
-#' The sigmoidal function given above would strictly reach 1 only asymptotically.
-#' Mizer instead sets the function equal to 1 already at the species' 
-#' maximum size, taken from the compulsory `w_max` column in the
-#' species parameter data frame. Also, for computational simplicity, any 
-#' proportion smaller than `1e-8` is set to `0`.
+#' The sigmoidal function given above would strictly reach 1 only
+#' asymptotically. Mizer instead sets the function equal to 1 already at a size
+#' taken from the `w_repro_max` column in the species parameter data frame, if it
+#' exists, or otherwise from the `w_max` column. Also, for computational
+#' simplicity, any proportion smaller than `1e-8` is set to `0`.
 #' }
 #' 
 #' \subsection{Investment into reproduction}{
 #' If the the energy available to a mature individual that is 
 #' invested into reproduction is not supplied via the `repro_prop` argument,
 #' it is set to the allometric form
-#' \deqn{{\tt repro\_prop}(w) = \left(\frac{w}{w_{max}}\right)^{m-n}.}{
-#'   repro_prop(w) = (w/w_max)^(m - n).}
+#' \deqn{{\tt repro\_prop}(w) = \left(\frac{w}{w_{mat_max}}\right)^{m-n}.}{
+#'   repro_prop(w) = (w/w_repro_max)^(m - n).}
 #' Here \eqn{n} is the scaling exponent of the energy income rate. Hence
 #' the exponent \eqn{m} determines the scaling of the investment into
 #' reproduction for mature individuals. By default it is chosen to be 
 #' \eqn{m = 1} so that the rate at which energy is invested into reproduction 
 #' scales linearly with the size. This default can be overridden by including a 
-#' column `m` in the species parameter dataframe. The maximum sizes
-#' are taken from the compulsory `w_max` column in the species parameter
-#' data frame.
+#' column `m` in the species parameter dataframe. The maximum sizes are taken
+#' from the `w_repro_max` column in the species parameter data frame, if it
+#' exists, or otherwise from the `w_max` column.
 #' 
 #' The total proportion of energy invested into reproduction of an individual
 #' of size \eqn{w} is then
@@ -222,11 +227,10 @@ setReproduction <- function(params, maturity = NULL,
         maturity[] <- 
             unlist(
                 tapply(params@w, seq_along(params@w),
-                       function(wx, w_max, w_mat, w_mat25) {
+                       function(wx, w_mat, w_mat25) {
                            U <- log(3) / log(w_mat / w_mat25)
                            return((1 + (wx / w_mat)^-U)^-1)
                        },
-                       w_max = species_params$w_max,
                        w_mat = species_params$w_mat,
                        w_mat25 = species_params$w_mat25
                 )
@@ -273,27 +277,30 @@ setReproduction <- function(params, maturity = NULL,
         }
     } else {
         # Set defaults for m
-        params <- set_species_param_default(params, "m", 1)
-        if (any(params@species_params$m < params@species_params[["n"]])) {
+        species_params <- set_species_param_default(species_params, "m", 1)
+        if (any(species_params$m < species_params[["n"]])) {
             stop("The exponent `m` must not be smaller than the exponent `n`.")
         }
+        # Set defaults for w_repro_max
+        species_params <- set_species_param_default(species_params, "w_repro_max",
+                                                    params@species_params$w_max)
 
         repro_prop <- array(
             unlist(
                 tapply(params@w, seq_along(params@w),
-                       function(wx, w_max, mn) (wx / w_max)^(mn),
-                       w_max = params@species_params$w_max,
-                       mn = params@species_params[["m"]] - 
-                           params@species_params[["n"]]
+                       function(wx, w_repro_max, mn) (wx / w_repro_max)^(mn),
+                       w_repro_max = species_params$w_repro_max,
+                       mn = species_params[["m"]] - species_params[["n"]]
                 )
             ), dim = c(nrow(species_params), length(params@w)))
+        repro_prop[repro_prop > 1] <- 1
     }
 
     psi <- params@maturity * repro_prop
     # psi should never be larger than 1
-    psi[params@psi > 1] <- 1
-    # Set psi for all w > w_max to 1
-    psi[outer(species_params$w_max, params@w, "<")] <- 1
+    psi[psi > 1] <- 1
+    # Set psi for all w > w_repro_max to 1
+    psi[outer(species_params$w_repro_max, params@w, "<")] <- 1
     assert_that(all(psi >= 0 & psi <= 1))
 
     # if the slot is protected and the user did not supply a new repro_prop

R/species_params.R

@@ -33,7 +33,7 @@
 #' * `k`, `ks` and `p` are used to set activity and basic metabolic rate, 
 #'   see [setMetabolicRate()].
 #' * `z0` is used to set the external mortality rate, see [setExtMort()].
-#' * `w_mat`, `w_mat25`, `w_max` and `m` are used to set the allocation to
+#' * `w_mat`, `w_mat25`, `w_repro_max` and `m` are used to set the allocation to
 #'   reproduction, see [setReproduction()].
 #' * `pred_kernel_type` specifies the shape of the predation kernel. The default
 #'   is a "lognormal", for other options see the "Setting predation kernel"
@@ -72,8 +72,9 @@
 #'   relationship \eqn{w = a l ^ b}.
 #'   
 #' If you have supplied the `a` and `b` parameters, then you can replace weight
-#' parameters like `w_max`, `w_mat`, `w_mat25` and `w_min` by their
-#' corresponding length parameters `l_max`, `l_mat`, `l_mat25` and `l_min`.
+#' parameters like `w_max`, `w_mat`, `w_mat25`, w_repro_max and `w_min` by their
+#' corresponding length parameters `l_max`, `l_mat`, `l_mat25`, `l_mat_max` and
+#' `l_min`.
 #'   
 #' The parameters that are only used to calculate default values for other
 #' parameters are:
@@ -511,7 +512,9 @@ get_ks_default <- function(params) {
 #' inconsistent.
 #' 
 #' If a `w_inf` column is given but no `w_max` then the value from `w_inf` is
-#' used. This is for backwards compatibility.
+#' used. This is for backwards compatibility. But note that the von Bertalanffy
+#' parameter `w_inf` is not the maximum size of the largest individual, but the
+#' asymptotic size of an average individual.
 #' 
 #' Some inconsistencies in the size parameters are resolved as follows:
 #' * Any `w_mat` that is not smaller than `w_max` is set to `w_max / 4`.
@@ -543,7 +546,8 @@ validSpeciesParams <- function(species_params) {
     # Check for misspellings ----
     misspellings <- c("wmin", "wmax", "wmat", "wmat25", "w_mat_25", "Rmax",
                       "Species", "Gamma", "Beta", "Sigma", "Alpha",
-                      "W_min", "W_max", "W_mat", "e_repro", "Age_mat")
+                      "W_min", "W_max", "W_mat", "e_repro", "Age_mat",
+                      "w_max_mat")
     query <- intersect(misspellings, names(sp))
     if (length(query) > 0) {
         warning("Some column names in your species parameter data ",
@@ -575,6 +579,7 @@ validSpeciesParams <- function(species_params) {
         sp <- sp %>%
             set_species_param_from_length("w_mat", "l_mat") %>%
             set_species_param_from_length("w_mat25", "l_mat25") %>%
+            set_species_param_from_length("w_repro_max", "l_mat_max") %>%
             set_species_param_from_length("w_max", "l_max") %>%
             set_species_param_from_length("w_min", "l_min")
     }

tests/testthat/test-setReproduction.R

@@ -10,7 +10,7 @@ test_that("setReproduction works", {
                           repro_prop = p2@psi)
     comment(p3@psi) <- comment(p2@psi)
     expect_unchanged(p2, p3)
-    p3 <- setReproduction(params, repro_prop = p2@psi)
+    p3 <- setReproduction(params, repro_prop = repro_prop(params))
     comment(p3@psi) <- comment(params@psi)
     expect_unchanged(params, p3)
     expect_error(setReproduction(params, RDD = "str"),