Round down w_repro_max to next-lower grid point for mizer > 2.5.3.9000

R/setReproduction.R

@@ -12,10 +12,11 @@
 #' 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. 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.
+#' energy is invested into reproduction and therefore all growth stops. For
+#' models created with mizer version > 2.5.3.9000, `w_repro_max` is rounded
+#' down to the nearest grid point to ensure 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
@@ -58,6 +59,12 @@
 #' from the `w_repro_max` column in the species parameter data frame, if it
 #' exists, or otherwise from the `w_max` column.
 #'
+#' For models created with mizer version > 2.5.3.9000, `w_repro_max` is 
+#' automatically rounded down to the next-lower point on the size grid. This
+#' prevents fish from growing into size classes beyond `w_repro_max` due to
+#' numerical integration behavior where a zero growth rate at the start of a
+#' size class can still lead to non-zero abundance in that class.
+#'
 #' The total proportion of energy invested into reproduction of an individual
 #' of size \eqn{w} is then
 #' \deqn{\psi(w) = {\tt maturity}(w){\tt repro\_prop}(w)}{psi(w) = maturity(w) * repro_prop(w)}
@@ -285,6 +292,22 @@ setReproduction <- function(params, maturity = NULL,
         # Set defaults for w_repro_max
         species_params <- set_species_param_default(species_params, "w_repro_max",
                                                     params@species_params$w_max)
+
+        # Round down w_repro_max to next-lower point on w grid for newer versions
+        # This prevents fish from growing into size classes beyond w_repro_max
+        # due to numerical integration behavior
+        if (params@mizer_version > "2.5.3.9000") {
+            for (i in seq_len(nrow(species_params))) {
+                # Find the largest w grid point that is <= w_repro_max
+                idx <- which(params@w <= species_params$w_repro_max[i])
+                if (length(idx) > 0) {
+                    species_params$w_repro_max[i] <- params@w[max(idx)]
+                }
+                # If w_repro_max is smaller than all grid points, leave it unchanged
+            }
+            # Save the rounded values back to params
+            params@species_params$w_repro_max <- species_params$w_repro_max
+        }
 
         repro_prop <- array(
             unlist(

tests/testthat/test-reproduction.R

@@ -55,3 +55,48 @@ test_that("noRDD returns rdi unchanged", {
     rdd <- noRDD(rdi, sp)
     expect_identical(rdd, rdi)
 })
+
+test_that("w_repro_max is rounded down to w grid for new versions", {
+    # Create a new params object with current version
+    params_new <- NS_params
+    # Set a w_repro_max that doesn't fall exactly on a grid point
+    # Use a value between two grid points
+    w_grid <- params_new@w
+    mid_idx <- length(w_grid) %/% 2
+    original_w_repro_max <- (w_grid[mid_idx] + w_grid[mid_idx + 1]) / 2
+    params_new@species_params$w_repro_max <- rep(original_w_repro_max, 
+                                                  nrow(params_new@species_params))
+
+    # Call setReproduction which should round down w_repro_max
+    params_new <- setReproduction(params_new)
+
+    # Check that w_repro_max was rounded down to a grid point
+    for (i in seq_len(nrow(params_new@species_params))) {
+        expect_true(params_new@species_params$w_repro_max[i] %in% params_new@w)
+        expect_true(params_new@species_params$w_repro_max[i] <= original_w_repro_max)
+    }
+
+    # Verify it's the largest grid point <= original value
+    expected_w_repro_max <- w_grid[mid_idx]
+    expect_equal(params_new@species_params$w_repro_max[1], expected_w_repro_max)
+})
+
+test_that("w_repro_max rounding only applies to new versions", {
+    # Create an old params object by setting mizer_version to an old value
+    params_old <- NS_params
+    params_old@mizer_version <- as.package_version("2.5.3.9000")
+
+    # Set a w_repro_max that doesn't fall exactly on a grid point
+    w_grid <- params_old@w
+    mid_idx <- length(w_grid) %/% 2
+    original_w_repro_max <- (w_grid[mid_idx] + w_grid[mid_idx + 1]) / 2
+    params_old@species_params$w_repro_max <- rep(original_w_repro_max, 
+                                                  nrow(params_old@species_params))
+
+    # Call setReproduction - should NOT round for old versions
+    params_old <- setReproduction(params_old)
+
+    # Check that w_repro_max was NOT rounded and remains as set
+    # For old versions, it should keep the off-grid value
+    expect_equal(params_old@species_params$w_repro_max[1], original_w_repro_max)
+})